(451c) 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, ε_oois 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 ρ_Bat B (and let  ρ_A  > ρ_B).
Now, we divide the phase space by interfaces, λ_iwhich
are loci of phase space points with the same value of ρ, say ρ_i.
Let the interfaces be numbered such that ρ_i > ρ_i+1.
Interface λ_Acorresponds to ρ_A
and λ_Bto ρ_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.