(451c) Liquid to Vapor TRANSITION of WATER UNDER Hydrophobic CONFINEMENT
- Conference: AIChE Annual Meeting
- Year: 2011
- Proceeding: 2011 Annual Meeting
- Group: Computational Molecular Science and Engineering Forum
- Time: Wednesday, October 19, 2011 - 9:12am-9:30am
Liquid to Vapor transition of water UNDER
Sumit Sharma, Pablo G. Debenedetti
of Chemical and Biological Engineering, Princeton University, Princeton, NJ
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 103 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.
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
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,
= (εcc . εoo)0.5 (1)
= (σcc + σoo)/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.
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.
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,
= ϕ(λ1|λA) ∏iP(λi+1|λi)
ϕ(λ1|λA) is the flux of
trajectories which leave λA and reach λ1.
P(λi+1|λi) is the
probability that a trajectory starting from λi would
reach λi+1 before reaching λA.
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 ϕ(λ1|λA) and collect N
configurations at λ1, which correspond to the first
crossing of λ1 by a trajectory coming from λA.
From each of the N configurations, C trajectories are shot (with
slightly perturbing velocities) to determine P(λ2|
λ1). The configurations at the interface λ2 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
Results and Discussion
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 109 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 ~1010 in evaporation rates. This
implies that the kinetics of evaporation of water in confined geometries is
strongly dependent on d.
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 kBT. 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 d2scaling 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.
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.
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.