Liquid to Vapor TRANSITION of WATER UNDER Hydrophobic CONFINEMENT

                Liquid to Vapor transition of water UNDER HYDROPHOBIC confinement

Sumit Sharma, Pablo G. Debenedetti

Department of Chemical and Biological Engineering, Princeton University, Princeton, NJ 08544

Introduction

                Molecular origins of phase separation of oil in water have been a topic of great interest as it is considered the basis of biological self assembly (Nicholls, Sharp et al. 1991; Sharp, Nicholls et al. 1991). Theoretical and simulation studies predict that the behavior of water molecules near hydrophobic solutes is solute-size-dependent (Lum, Chandler et al. 1999; Rajamani, Truskett et al. 2005). At ambient conditions, water molecules are able to retain their hydrogen (H) bonds around hydrophobic solutes ≤ 1 nm (Pangali, Rao et al. 1979). The solubility of small hydrophobic solutes is shown to be proportional to the free energy of creating a cavity of the size of the solute in water (Hummer, Garde et al. 1998). For larger solutes, water molecules are unable to retain their H-bonds close to the solute and hence form a soft, vapor-like interface (Wallqvist and Berne 1995). Simple thermodynamical arguments show that water confined between two large repulsive surfaces separated by distance d should vaporize even when d is as large as 10³ nm (Lum, Chandler et al. 1999). The fact that evaporation transitions at such macroscopic length scales do not occur in finite times (Christenson and Claesson 2001) indicates that kinetics of evaporation in hydrophobic confinement are of fundamental importance and merit computational scrutiny.

                In this study, we estimate rates and activation energy for liquid to vapor transition of water confined between hydrophobic surfaces using Forward Flux Sampling (FFS) (Allen, Frenkel et al. 2006). FFS allows a good sampling of the liquid to vapor phase transition, which, being a rare event in the simulation time scales, is hard to achieve in equilibrium simulations. We found that the activation energy, ΔG scales as d for the case when the size of the hydrophobic solutes is ~ 1 nm. The transition state ensemble comprises of configurations in which almost the entire confined region is devoid of water molecules.

Simulation system and methods

We use the SPC/E model of water (Berendsen, Grigera et al. 1987). The hydrophobic surfaces were constructed with Lennard Jones (LJ) atoms arranged in a hexagonal lattice. This arrangement of atoms mimics that of carbon atoms in graphene (Choudhury and Pettitt 2005). The LJ potential parameters of carbon atoms for the Amber force field (Cornell, Cieplak et al. 1995) are: ε_cc = 0.086 kcal mol^-1 and σ_cc = 3.4 Å. The carbon atoms of the surface interact with the oxygen atoms of water molecules with LJ potential. As per the Lorentz – Berthelot rules,

ε_oc = (ε_cc . ε_oo)^0.5                                                                                                                                                                                                                                           (1)

σ_oc = (σ_cc + σ_oo)/2                                                                                                                                                                        (2)

where, ε_oo is the SPC/E water oxygen – oxygen LJ ε parameter and σ_oo is the SPC/E water LJ σ parameter. For SPC/E water, equations (1) and (2) give ε_oc = 0.1156 kcal/mol and σ_oc = 3.283 Å.  This carbon – oxygen potential leads to a very strong binding of water molecules to the surfaces (Hummer, Rasaiah et al. 2001; Choudhury and Pettitt 2005).  Hence, we simulated a weaker (more hydrophobic) wall – oxygen potential with ε_ow = ε_oc / 4 and σ_ow = σ_oc.

                Molecular dynamics (MD) simulations were conducted in the isothermal-isobaric (NPT) ensemble. 2329 water molecules were taken in a cubic simulation box, periodic in the x, y and z directions. Each wall / hydrophobic surface, made up of 60 atoms, had the lateral dimensions of 9.05 x 0.05 Å. The walls were kept fixed and were were placed parallel at a distance d from each other.

Forward Flux Sampling

Suppose a system has two locally stable states designated A and B separated by a free energy barrier. Let us identify a system property, say ρ, that changes in value from ρ_A at A to ρ_B at B (and let  ρ_A  > ρ_B). Now, we divide the phase space by interfaces, λ_i which are loci of phase space points with the same value of ρ, say ρ_i. Let the interfaces be numbered such that ρ_i > ρ_i+1. Interface λ_A corresponds to ρ_A and λ_B to ρ_B. The rate of the transition from A to B will be,

Rate = ϕ(λ₁|λ_A) ∏_iP(λ_i+1|λ_i)                                                                                                                               (3)

Where ϕ(λ₁|λ_A) is the flux of trajectories which leave λ_A and reach λ₁. P(λ_i+1|λ_i) is the probability that a trajectory starting from λ_i would reach λ_i+1 before reaching λ_A.  

In order to calculate the rate of transition from liquid to vapor state for water confined between two hydrophobic surfaces, ρ is the density of the water molecules in the confined region. From an equilibrium simulation (40 ns), we determine ϕ(λ₁|λ_A) and collect N configurations at λ₁, which correspond to the first crossing of λ₁ by a trajectory coming from λ_A. From each of the N configurations, C trajectories are shot (with slightly perturbing velocities) to determine P(λ₂| λ₁). The configurations at the interface λ₂ are further propagated to the next interface. The process of propagating the trajectories from one interface to the next is continued till the vapor phase is reached.

Results and Discussion

            Evaporation rates were calculated for d =9.0 Å to 14.0 Å. For d=9.0 Å, the confined region sampled both liquid and vapor phase-space in the equilibrium simulation. For d ≥ 9.8 Å, the confined region only sampled the liquid phase in the equilibrium simulation. The evaporation rate for d=9.0 Å was found to be 1.67 x 10⁹ nm^-2s^-1, whereas for d=14.0 Å, it was found to be 0.058 nm^-2s^-1. Therefore, an increase in d less than two times the diameter of a water molecule leads to a decrease of the order ~10¹⁰ in evaporation rates. This implies that the kinetics of evaporation of water in confined geometries is strongly dependent on d.

The activation energy to the evaporation event for d=9.8 Å was determined by doing evaporation rate calculations at T=298, 348 and 398 K. The activation energy was found to be ~ 45 k_BT. Assuming a constant pre-factor, activation energies for other values of d were predicted from the calculated evaporation rates. The predicted values of the activation energies were validated by calculating the activation energy at d=11.0 Å. The activation energy to evaporation was found to scale linearly with d. The transition state ensemble comprised of configurations in which the vapor bubble is almost of the same volume as the confined region. This observation implies that the free energy barrier to evaporation should be proportional to the free energy required to create a vapor bubble of approximately the same volume as the confined region, which explains the linear dependence of activation energy on d. A lattice gas model calculation to determine free energy barriers to evaporation between infinite plates showed a d² scaling of the activation energy (Luzar 2004).  Therefore, it appears that there should be transition in the scaling factor from linear to quadratic as the size of the hydrophobic plates is increased.  

References

Allen, R. J., D. Frenkel, et al. (2006). "Simulating rare events in equilibrium or nonequilibrium stochastic systems." J. Chem. Phys. 124(2): 024102.

Berendsen, H. J. C., J. R. Grigera, et al. (1987). "The Missing Term in Effective Pair Potentials." J. Phys. Chem. 91(24): 6269-6271.

Choudhury, N. and B. M. Pettitt (2005). "On the mechanism of hydrophobic association of nanoscopic solutes." J. Am. Chem. Soc. 127(10): 3556-3567.

Christenson, H. K. and P. M. Claesson (2001). "Direct measurements of the force between hydrophobic surfaces in water." Adv. Colloid Interface Sci. 91(3): 391-436.

Cornell, W. D., P. Cieplak, et al. (1995). "A 2nd Generation Force-Field for the Simulation of Proteins, Nucleic-Acids, and Organic-Molecules." J. Am. Chem. Soc. 117(19): 5179-5197.

Hummer, G., S. Garde, et al. (1998). "Hydrophobic effects on a molecular scale." J. Phys. Chem. B 102(51): 10469-10482.

Hummer, G., J. C. Rasaiah, et al. (2001). "Water conduction through the hydrophobic channel of a carbon nanotube." Nature 414(6860): 188-190.

Lum, K., D. Chandler, et al. (1999). "Hydrophobicity at small and large length scales." J. Phys. Chem. B 103(22): 4570-4577.

Luzar, A. (2004). "Activation barrier scaling for the spontaneous evaporation of confined water." J. Phys. Chem. B 108(51): 19859-19866.

Nicholls, A., K. A. Sharp, et al. (1991). "Protein Folding and Association - Insights from the Interfacial and Thermodynamic Properties of Hydrocarbons." Proteins: Struct., Funct., Genet. 11(4): 281-296.

Pangali, C., M. Rao, et al. (1979). "Hydrophobic Hydration around a Pair of Apolar Species in Water." J. Chem. Phys. 71(7): 2982-2990.

Rajamani, S., T. M. Truskett, et al. (2005). "Hydrophobic hydration from small to large lengthscales: Understanding and manipulating the crossover." Proc. Natl. Acad. Sci. U. S. A. 102(27): 9475-9480.

Sharp, K. A., A. Nicholls, et al. (1991). "Extracting Hydrophobic Free-Energies from Experimental-Data - Relationship to Protein Folding and Theoretical-Models." Biochemistry 30(40): 9686-9697.

Wallqvist, A. and B. J. Berne (1995). "Computer simulation of hydrophobic hydration forces on stacked plates at short range." J. Phys. Chem. 99(9): 2893-2899.