^{1}and Bernhardt L. Trout

^{1,a)}

### Abstract

The process of water's evaporation at its liquid/air interface has proven challenging to study experimentally and, because it constitutes a rare event on molecular time scales, presents a challenge for computer simulations as well. In this work, we simulated water's evaporation using the classical extended simple point charge model water model, and identified a minimum free energy path for this process in terms of 10 descriptive order parameters. The measured free energy change was 7.4 kcal/mol at 298 K, in reasonable agreement with the experimental value of 6.3 kcal/mol, and the mean first-passage time was 1375 ns for a single molecule, corresponding to an evaporation coefficient of 0.25. In the observed minimum free energy process, the water molecule diffuses to the surface, and tends to rotate so that its dipole and one O–H bond are oriented outward as it crosses the Gibbs dividing surface. As the water molecule moves further outward through the interfacial region, its local density is higher than the time-averaged density, indicating a local solvation shell that protrudes from the interface. The water molecule loses donor and acceptor hydrogen bonds, and then, with its dipole nearly normal to the interface, stops donating its remaining hydrogen bond. At that point, when the final, accepted hydrogen bond is broken, the water molecule is free. We also analyzed which order parameters are most important in the process and in reactive trajectories, and found that the relative orientation of water molecules near the evaporating molecule, and the number of accepted hydrogen bonds, were important variables in reactive trajectories and in kinetic descriptions of the process.

The authors would like to thank Dr. Erik E. Santiso for guidance in modifying NAMD, for sharing a software utility to effect SMCV iterations, and for many helpful technical discussions. Molecular graphics were produced using VMD.^{74}

The authors would like to kindly acknowledge support from the DuPont-MIT Alliance, and from the NIH/NIGMS Biotechnology Training Program at MIT.

I. INTRODUCTION

II. ORDER PARAMETERS AND SIMULATION SETUP

III. RESULTS

A. Minimum free energy path

B. Evaporation thermodynamics and kinetics

1. Comparison to experimental results

IV. IDENTIFICATION OF IMPORTANT ORDER PARAMETERS

A. Principal component analysis

B. Directional analysis

C. Examining MFPT as a function of order parameters

V. CONCLUSIONS AND OUTLOOK

A. Implications for additive design

### Key Topics

- Evaporation
- 50.0
- Hydrogen bonding
- 40.0
- Free energy
- 29.0
- Gas liquid interfaces
- 13.0
- Liquid surfaces
- 11.0

##### B01D1/00

## Figures

Rendering of 1025 water molecules in the 31 × 31 × (4 × 31 Å) unit cell. The *z*-axis is in the horizontal direction.

Rendering of 1025 water molecules in the 31 × 31 × (4 × 31 Å) unit cell. The *z*-axis is in the horizontal direction.

Time-averaged density profiles from four 2.0-ns simulations, with frames recorded every 5 ps.

Time-averaged density profiles from four 2.0-ns simulations, with frames recorded every 5 ps.

The dipole-dipole angle η and its distribution in bulk water. (a) The angle η between dipole vectors (blue cylinders) of two water molecules; for clarity, the dipole vectors are translated and reproduced. (b) Distribution of cosine of dipole-dipole angle η for water molecules with indicated O–O separation distances in a 1.0-ns bulk SPC/E water simulation.

The dipole-dipole angle η and its distribution in bulk water. (a) The angle η between dipole vectors (blue cylinders) of two water molecules; for clarity, the dipole vectors are translated and reproduced. (b) Distribution of cosine of dipole-dipole angle η for water molecules with indicated O–O separation distances in a 1.0-ns bulk SPC/E water simulation.

The two “absolute” orientation variables θ and ω, used to define *q* _{4} = cos (θ) and *q* _{5} = cos ^{2}ω. The two angles are defined in relation to the interfacial normal vector and the evaporating molecule's dipole vector and the molecular normal, respectively. (a) Schematic illustration of the two “absolute orientation” angles. (b) Rendering of the dipole vector μ (blue cylinder), the molecular normal vector ν (yellow cylinder), and the interfacial normal (gray cylinder).

The two “absolute” orientation variables θ and ω, used to define *q* _{4} = cos (θ) and *q* _{5} = cos ^{2}ω. The two angles are defined in relation to the interfacial normal vector and the evaporating molecule's dipole vector and the molecular normal, respectively. (a) Schematic illustration of the two “absolute orientation” angles. (b) Rendering of the dipole vector μ (blue cylinder), the molecular normal vector ν (yellow cylinder), and the interfacial normal (gray cylinder).

Examples of restraint force and order parameter convergence from image 9 of string 4. (a) Restraint force for order parameter . (b) Value of order parameter

Examples of restraint force and order parameter convergence from image 9 of string 4. (a) Restraint force for order parameter . (b) Value of order parameter

Values of tetrahedrality order parameters in each Voronoi dynamics image. Bars indicate one semi-standard deviation over simulation trajectory.

Values of tetrahedrality order parameters in each Voronoi dynamics image. Bars indicate one semi-standard deviation over simulation trajectory.

Frechét distance from initial string, and Frechét distance from each string's predecessor. See text for notes about methodological adjustments at string 18 and string 29.

Frechét distance from initial string, and Frechét distance from each string's predecessor. See text for notes about methodological adjustments at string 18 and string 29.

Minimum free energy path for evaporation, along with Voronoi cell boundaries between images, projected onto two order parameter dimensions at a time. The point labeled “GDS” is image 9, which contains the plane *q* _{0} = *z* _{ GDS }, the Gibbs dividing surface. The free energy changed most dramatically over the images (numbers 9–14), highlighted with white centers. Note that Voronoi cell boundaries do not necessarily appear normal to the string because they respect the scaling of order parameters (see text), and because of the plots' axis scaling.

Minimum free energy path for evaporation, along with Voronoi cell boundaries between images, projected onto two order parameter dimensions at a time. The point labeled “GDS” is image 9, which contains the plane *q* _{0} = *z* _{ GDS }, the Gibbs dividing surface. The free energy changed most dramatically over the images (numbers 9–14), highlighted with white centers. Note that Voronoi cell boundaries do not necessarily appear normal to the string because they respect the scaling of order parameters (see text), and because of the plots' axis scaling.

Snapshots from the frames in images 8–13 in which the system was closest to its OP target values, as measured by minimal restraint energy. In these images, *q* _{0} = *z* varied from 13 to 21 Å. (a) Water molecule's orientation during evaporation; the molecule's position and orientation (subject to translation) from images 9–13 is shown in a single figure. The other molecules' configuration is from image 9. The blue and yellow vectors are the dipole and molecular normal vectors, respectively. The transparent surface is the water surface. (b) Hydrogen bonds that an evaporating water molecule donates (yellow) and accepts (purple) in images 8–13. The blue arrow is the dipole vector of the evaporating molecule.

Snapshots from the frames in images 8–13 in which the system was closest to its OP target values, as measured by minimal restraint energy. In these images, *q* _{0} = *z* varied from 13 to 21 Å. (a) Water molecule's orientation during evaporation; the molecule's position and orientation (subject to translation) from images 9–13 is shown in a single figure. The other molecules' configuration is from image 9. The blue and yellow vectors are the dipole and molecular normal vectors, respectively. The transparent surface is the water surface. (b) Hydrogen bonds that an evaporating water molecule donates (yellow) and accepts (purple) in images 8–13. The blue arrow is the dipole vector of the evaporating molecule.

Transition frequencies from home cell to other cells for four selected images during Voronoi dynamics simulations. Simulations for other images exhibited transitions to sequential cells, as in panels (a) and (b). (a) Image 2. (b) Image 10. (c) Image 14. (d) Image 17.

Transition frequencies from home cell to other cells for four selected images during Voronoi dynamics simulations. Simulations for other images exhibited transitions to sequential cells, as in panels (a) and (b). (a) Image 2. (b) Image 10. (c) Image 14. (d) Image 17.

Free energy measured through Voronoi milestoning, as a function of order parameter *q* _{0} = relative *z*-position (left) and order parameter *q* _{1} = local density (right). The local density achieves a minimum value in the vapor phase when only the evaporating molecule itself is contributing to the local density.

Free energy measured through Voronoi milestoning, as a function of order parameter *q* _{0} = relative *z*-position (left) and order parameter *q* _{1} = local density (right). The local density achieves a minimum value in the vapor phase when only the evaporating molecule itself is contributing to the local density.

Free energy as a function of order parameters 6 and 7.

Free energy as a function of order parameters 6 and 7.

Free energy profile, along with average system energy values. Error bars are 1.5 standard errors.

Free energy profile, along with average system energy values. Error bars are 1.5 standard errors.

Mean first passage time to the final milestone as a function of order parameter *q* _{0} = relative *z*-position (top) and as a function of *q* _{6} and *q* _{7}, the number of hydrogen bonds accepted and donated (bottom).

Mean first passage time to the final milestone as a function of order parameter *q* _{0} = relative *z*-position (top) and as a function of *q* _{6} and *q* _{7}, the number of hydrogen bonds accepted and donated (bottom).

Schematic showing contributing trajectories in both reactant-to-product and product-to-reactant directions (solid curves) and non-contributing trajectories (dashed curves).

Schematic showing contributing trajectories in both reactant-to-product and product-to-reactant directions (solid curves) and non-contributing trajectories (dashed curves).

Projection of contributing trajectory segments in the forward (evaporating) direction onto principle components. Image centers (Voronoi support points) are the black points, while the trajectories from images 10–13 are shown in alternating shades of gray. The rightmost point represents the final, vapor-phase image.

Projection of contributing trajectory segments in the forward (evaporating) direction onto principle components. Image centers (Voronoi support points) are the black points, while the trajectories from images 10–13 are shown in alternating shades of gray. The rightmost point represents the final, vapor-phase image.

Summary of PCA results. (a) Amount of variance explained by principal components. (b) Contributions of each order parameter to principal components 1 and 2.

Summary of PCA results. (a) Amount of variance explained by principal components. (b) Contributions of each order parameter to principal components 1 and 2.

Coefficients of order parameters for the five models *A*–*E* with best BIC values, and the linear model containing all order parameters. The coefficients are for the normalized order parameters, and the error bars are 95% confidence intervals.

Coefficients of order parameters for the five models *A*–*E* with best BIC values, and the linear model containing all order parameters. The coefficients are for the normalized order parameters, and the error bars are 95% confidence intervals.

Observed and fitted values of the MFPT values at milestones, plotted against (a) local density and (b) number hydrogen bonds accepted. The fitted data were from model *B* of Table V and Figure 18 . (a) MFPT and *z*-position. (b) MFPT and number hydrogen bonds accepted.

Smooth weighting functions used for calculating local density and local relative orientational order (top), and number of hydrogen bonds donated and accepted (bottom).

Smooth weighting functions used for calculating local density and local relative orientational order (top), and number of hydrogen bonds donated and accepted (bottom).

## Tables

Description of order parameters used to describe state of water molecule near interface. Order parameters 8 and 9 have definitions that do not permit the imposition of forces, but listed force constants were used to calculate distances in order parameter-space.

Description of order parameters used to describe state of water molecule near interface. Order parameters 8 and 9 have definitions that do not permit the imposition of forces, but listed force constants were used to calculate distances in order parameter-space.

Comparison of simulation measurements to experimental values for the evaporation or “desolvation” process at 298 K. All values are in kcal/mol.

Comparison of simulation measurements to experimental values for the evaporation or “desolvation” process at 298 K. All values are in kcal/mol.

Order parameter components of first and second principle components in analysis of contributing trajectories in Images 10–13. The three largest components in PC1 and the two largest in PC2 are highlighted with boldface type.

Order parameter components of first and second principle components in analysis of contributing trajectories in Images 10–13. The three largest components in PC1 and the two largest in PC2 are highlighted with boldface type.

Results of local direction analysis for forward-directed transitions in Images 8–14. The two largest components of the mean vector in each image are underlined.

Results of local direction analysis for forward-directed transitions in Images 8–14. The two largest components of the mean vector in each image are underlined.

Best models of τ^{′} = (1 − τ/τ_{ bulk }) with different numbers of order parameters used. The combinations are sorted by BIC value.

Best models of τ^{′} = (1 − τ/τ_{ bulk }) with different numbers of order parameters used. The combinations are sorted by BIC value.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content