Molecular dynamics simulation of the density and

THE JOURNAL OF CHEMICAL PHYSICS 124, 204715 共2006兲

Molecular dynamics simulation of the density and surface tension of water by particle-particle particle-mesh method Bo Shi, Shashank Sinha, and Vijay K. Dhira兲 Department of Mechanical and Aerospace Engineering, University of California, Los Angeles, California 90095

共Received 9 December 2005; accepted 30 March 2006; published online 31 May 2006兲 In this work, molecular dynamics simulation is performed to study the density and surface tension of water for a range of temperatures from 300 to 600 K. The extended simple point charge interaction potential for water is used. The particle-particle particle-mesh method, which automatically includes untruncated long-range terms, is used for the Lennard-Jones and the Coulombic terms. The results show that the long-range correction for the Lennard-Jones term is very important for the calculation of surface tension. It is found that the calculated density and surface tension of water fit well with experimental data for temperatures less than 500 K. Near the critical temperature, the simulation results are off from the experimental data. © 2006 American Institute of Physics. 关DOI: 10.1063/1.2199849兴 I. INTRODUCTION

Due to its ubiquity in our environment, water is certainly the liquid most investigated through molecular dynamics, since computer simulation was introduced in the 1960s. After the Bernal-Fowler 共BF兲 model,1 a variety of interaction potentials for the simulation of water have been proposed. By using these potentials, many studies have been conducted to predict the density, surface tension, and other physical properties of water. Guillot2 presented a review of the density calculations of water using molecular dynamics. For ambient conditions, several researchers have successfully calculated the density by using some water potentials, such as simple point charge 共SPC兲,3 TIP3P, TIP4P,4 TIP5P, extended simple point charge 共SPC/E兲,5,6 and so on. However, all these popular models failed when the system temperature was raised to near the critical temperature. Liquid density was found to be lower and vapor density to be higher than the experimental values.7 The reason for the error is that most of the interaction potentials for water are built empirically for the ambient conditions and are not correct at higher temperatures. There is a large variation in prediction for surface tension of water. In 1988, Matsumoto and Kataoka8 reported a value of 30.5 mN/ m at 300 K, which is much smaller than the experimental value of 72 mN/ m. In 1987, Wilson et al.9 also calculated the surface tension for TIP4P model of water at 325 K. They reported a surface tension value of 132 mN/ m compared to the experimental value of 68 mN/ m. Later, they claimed that their calculations had errors and reported a corrected value of 59 mN/ m. Another result reported by Lie et al.,10 also underestimated to onethird of the experimental value at 300 K. Some researchers have reported higher values. Zhu et al.11 calculated the surface tension of SPC water at 298 K and found a value of 123.6 mN/ m, about 70% higher than the experimental value a兲

Electronic mail: [email protected]

0021-9606/2006/124共20兲/204715/7/$23.00

of 72.75 mN/ m. They claimed that their calculation had a large statistical uncertainty and the error may be caused by the SPC model itself. In 1997, Zakharov et al.12 studied the surface tension for a water droplet using TIP4P potential model. They evaluated surface tension as 54 mN/ m at 300 K. This suggested that the error in surface tension did not only exist in the plane surface, but also in other geometric shapes. In 1995, while using SPC/E model to simulate water, Alejandre et al.13 found that the long-range interaction is very important for surface tension. They introduced the Ewald sum to correct the long-range interaction of the Coulombic term and found that their results were closer to experimental values. However, this method is only applied for the Coulombic term but not for the Lennard-Jones term in the SPC/E model. To compensate the discrepancies caused by the long-range interaction of the Lennard-Jones term, they introduced an empirical tail correction to their calculation. Recently, while simulating an argon system, Sinha et al.14,15 found that truncating long-range terms of the Lennard-Jones potential function at 4.5␴ would cause errors as high as 15% in surface tension, which can be eliminated by using particle-particle particle-mesh 共PPPM兲 method for long-range terms. They also suggested that the long-range interaction of the Lennard-Jones term in SPC/E model of water is not negligible. In this paper, the PPPM method is applied to minimize the error in long-range terms for both Coulombic and Lennard-Jones terms and no tail correction is needed.

II. THEORY AND APPROACH A. PPPM method

PPPM or P3M algorithms are a class of hybrid algorithms developed by Hockney and Eastwood.17 The basic 124, 204715-1

© 2006 American Institute of Physics

Downloaded 05 Jun 2006 to 128.97.3.167. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp

204715-2

J. Chem. Phys. 124, 204715 共2006兲

Shi, Sinha, and Dhir TABLE I. f p and g p in Ref. 16. p

f p共x兲

g p共x兲

1 2 3 4 5 6 7 8

exp共−x2兲 / 共冑␲x2兲共冑␲ / x兲erfc共x兲 2⌫关0 , x2兴 / 冑␲ 2关exp共−x2兲 − x冑␲ erfc共x兲兴共4 / 3冑␲兲关exp共−x2兲 − x2⌫关0 , x2兴兴 1 关共1 − 2x2兲exp共−x2兲 + 2x3冑␲ erfc共x兲兴 3 共16x4 / 15冑␲兲关共1 − x2兲exp共−x2兲 + x4⌫关0 , x2兴兴 1 2 2 4 冑 5 45 关exp共−x 兲共3 − 2x + 4x 兲 − 4 ␲x erfc共x兲兴

erfc共x兲 ⌫关0 , x2兴 erfc共x兲 + 关2x exp共−x2兲兴 / 冑␲ exp共−x2兲共1 + x2兲 2 exp共−x 兲共6x + 4x2兲 + 3冑␲ erfc共x兲 exp共−x2兲共1 + x2 + x4 / 2兲 erfc共x兲 + 关exp共−x2兲共30x + 20x3 + 8x5兲兴 / 15冑␲ exp共−x2兲共1 + x2 + x4 / 2 + x6 / 6兲

idea of the PPPM method is to split the slow decaying potential function such as the 1 / rn term into two parts by the following identity: 1 f共r兲 1 − f共␶兲 = n + . rn r rn

共1兲

There are two requirements for the splitting. 共1兲 The first part f共r兲 / rn, called the short-range term, should decay exponentially with r increasing and become negligible 共or even zero兲 beyond some small cutoff radius rc, thus the particleparticle 共PP兲 method with the neighbor-list method could then be used without introducing large cutoff errors. 共2兲 The second part 关1 − f共r兲兴 / rn, called the long-range term, should be a slowly varying function for all r, so that its Fourier transform can be represented by only a few k vectors with 兩k兩艋 kmax in Fourier space. This permits an efficient and accurate calculation by using the particle-mesh 共PM兲 method without introducing large aliasing errors. Since the field equations are linear, the sum of these two parts gives the solution for the potential function of the original problem. However, these two requirements on the function f mentioned above leave a large freedom of choice. One choice was proposed by Essmann et al.16 Using the definition of the Euler gamma function 共⌫兲 and the Fourier integral expansion of the Gaussian, the following equation is derived: ⌫共p/2兲 = rp

冕

␤2

t

p/2−1

2

exp共− r t兲dt +

冕

⬁

t

p/2−1

2

exp共− r t兲dt.

␤2

0

共2兲 Here ␤ is an arbitrary positive number and p is a positive integer. Following the derivation of Essmann et al. we get 1 = ␲3/2␤ p−3 rp

冕

R3

f p共␲兩u兩/␤兲exp共− 2␲iu · r兲d3u +

g p共␤r兲 , rp 共3兲

where u is the coordinates of points in R3 such that −1 / 2 艋 ai · u 艋 1 / 2 and ai is the lattice vectors of the unit cell. Functions f p and g p are given by f p共x兲 =

2x p−3 ⌫共p/2兲

冕

⬁

x

s2−p exp共− s2兲ds,

共4兲

g p共x兲 =

2 ⌫共p/2兲

冕

⬁

s p−1 exp共− s2兲ds.

共5兲

x

In particular, for the cases p = 1 – 8, we have tabulated f p共x兲 and g p共x兲 in Table I. After splitting, the short-range part can be calculated using particle-particle method while employing a neighbor-list method. The long-range part can be computed by an optimized particle-mesh method following these four steps: 共1兲共2兲共3兲共4兲

Assign charge to the mesh. Solve the field equation on the mesh. Calculate the mesh-defined force field. Interpolate to find forces on the particles.

In the second step, fast Fourier transform 共FFT兲 is used to solve Poisson’s equation with O共N log N兲. The calculation is called optimized because an optimized Green function is applied to make sure that the total error is minimal. The details can be found in Hockney and Eastwood.17

B. Force Evaluation

Consider the SPC/E potential model,

␾共rij兲 = 4⑀LJ

冋冉冊冉冊册 ␴ rij

12

−

␴ rij

6

+

1 q iq j . 4␲⑀0⑀r rij

共6兲

Here i and j are atoms in different molecules, rij is the distance between i and j. The Lennard-Jones 共LJ兲 parameters are ⑀LJ = 1.0797⫻ 10−21 J for oxygen-oxygen, zero for both hydrogen-oxygen and hydrogen-hydrogen, ␴ = 3.166 ⫻ 10−10 m. In Eq. 共6兲, ⑀0 = 8.854187⫻ 10−12 F / m is the permittivity of vacuum, and ⑀r is the relative permittivity which is 1 for the vacuum case. For oxygen atom, the charge is qO = −0.8476e and for hydrogen, qH = 0.4238e, where e = 1.60219⫻ 10−19 C. This potential function has a Lennard-Jones term and a Coulombic term. The Lennard-Jones term is only effective for the oxygen-oxygen interaction and Coulombic term is effective for all the atoms. Notice that for the Coulombic term, p = 1. Following Table I, by dropping the constant 1 / 4␲⑀0⑀r, the short-range part of the Coulombic term of the potential on atom i is


204715-3

J. Chem. Phys. 124, 204715 共2006兲

Surface tension of water by PPPM method

␾C-SR,i = 兺 qiq j j

erfc共␤r兲 . r

= 兺 q iq j j

冉冑

exp共− ␤2r2兲 +

␲

冊

erfc共␤r兲 r . r r2

共8兲

This term decays exponentially with increasing r so that it can be calculated by PP method. Also following Table I, the long-range part of the Coulombic term of the potential can be expressed as

␾C-LR = ␲3/2␤−2

冕

R3

f 1共␲兩u兩/␤兲exp共− 2␲iu · r兲d3u.

共9兲

By introducing a reciprocal k vector from the discrete set 兵2␲n / L : n 苸 Z3其, k = 2␲u, Eq. 共9兲 can be written as

␾C-LR =

1 4␲ 2 2 兺 e−k /4␤ 兩˜␳共k兲兩2 , 2V k⫽0 k2

共10兲

where the Fourier transformed charge density ˜␳共k兲 is defined as N

˜␳共k兲 = 兺 q jeik·r j .

共11兲

j=1

Note that the split for p = 1 gives the well known Ewald summations. To calculate the long-range accurately, an optimized Green function is introduced to minimize the error, ˆ G opt共k兲 =

˜ 共k兲 · 兺 ˜ 2共k + 共2␲/h兲m兲R ˜ 共k + 共2␲/h兲m兲 D U m苸Z3 ˜ 共k兲兩2关兺 ˜ 2共k + 共2␲/h兲m兲兴2 兩D U m苸Z3

,

共12兲 ˜ 共k兲 is the Fourier transform of the employed differwhere D ˜ 共k兲 = W ˜ 共k兲 / V is the Fourier transform entiation operator, U c of the charge assignment function per the volume of one mesh cell, and

冉

˜ 共k兲 = h sin共kh/2兲 W kh/2

冊

P

共13兲

is the Pth order assignment scheme which assigns the ˜ 共k兲 is charges to the grid and h = L / N M is the mesh spacing. R the Fourier transform of the true reference force, given by ˜ 共k兲 = − ikg ˜ 共k兲˜␥共k兲. R

共14兲

Here ˜g共k兲 and ˜␥共k兲 depend on the choice of splitting function f. For p = 1 case, ˜g共k兲 = and

4␲ , k2

i is FC-LR,i = − h3

兺 r 苸M p

2␤

共16兲

Therefore, the final long-range Coulombic force on atom

Therefore, the short-range Coulombic force on atom i can be obtained as

⳵␾C-SR,i FC-SR,i = − ⳵r

˜␥共k兲 = exp共− k2/4␤2兲.

共7兲

共15兲

⳵␳M 共r p兲关␳ M 쐓Gopt兴共r p兲, ⳵ri

共17兲

where ␳ M is the mesh based charge density which is obtained from the charge assignment. The term ␳ M 쐓Gopt can be expressed as

ឈ 关FFT ជ 关␳ 兴 ⫻ FFT ជ 关G 兴兴, ␳M 쐓Gopt = FFT M opt

共18兲

ជ is the discrete fast Fourier transform where the operator FFT ឈ is the inverse discrete fast Fourier transform. As ␳ , and FFT M

Gopt are known, the long-range force can be easily calculated. Notice that, because Gopt can be computed in advance, the PPPM method is much more efficient than the direct Ewald summations. The details of the calculation of Coulombic force by using PPPM method can be found in Deserno and Holm.18 It should be mentioned that there is an extra term called the “dipole” in Ewald summations which cannot be obtained from our splitting. This term is not large but should also be included in the calculation, Fdipole,i = −

4␲qi 兺 q jr j , 共1 + 2⑀⬘兲V j

共19兲

where ⑀⬘ is the dielectric constant, and for vacuum ⑀⬘ = 1. The Lennard-Jones term has two parts: 共1兲 the repulsive term, 1 / r12, to model the electron cloud at small distances and 共2兲 the attractive dispersive term 1 / r6. The 1 / r12 term decays fast enough to be evaluated by particle-particle method with a finite cutoff radius. The 1 / r6 term, however, decays slowly and has to be evaluated by the PPPM method. If we include the repulsive term 1 / r12, the total short-range Lennard-Jones force is FLJ-SR = −

冉冉冊冉冊

␴ d 4⑀ dr r

= 4⑀

␴ r

12

12

− 4⑀␴6

g6共␤r兲 r6

冊

6g6共␤r兲 − rg6⬘共␤r兲 12 − 4⑀␴6 , r r7

共20兲

where g6 is obtained from Table I with p = 6, and g6⬘ is the derivative of g6. The calculation of the long-range part of the LJ force is similar to that of the Coulombic force. The only difference is that because p is six, f 6 should be used, and ˜g共k兲 and ˜␥共k兲 can be derived from f 6, respectively. As the short-range parts and long-range parts of the Coulombic and Lennard-Jones terms are all obtained, the total force on each atom can be easily calculated by direct summation.

C. Surface tension evaluation

The hydrodynamic definition of surface tension is the isothermal work of formation per unit area of interface. At


204715-4

J. Chem. Phys. 124, 204715 共2006兲

Shi, Sinha, and Dhir

the atomic scale, it can be expressed as the integrated imbalance of normal and tangential pressures near the interface.

␥=

冕

phase 2

关PN共z兲 − PT共z兲兴dz.

For a plane interface perpendicular to the z axis, tangential and normal pressure components are defined as PN共z兲 = Pzz共z兲,

共22兲

PT共z兲 = 共Pxx共z兲 + Pyy共z兲兲/2.

共23兲

and

For pair potentials with fij = f ijrij / rij, the normal and tangential pressures are

冋

1 V

1 Pyy = V

册

兺 m jv2xj + 兺兺 rxij f xij , j

i⬎j

冋兺 j

j

册

共25兲

册

共26兲

j

1 V

冋

兺 rzij f zij 兺j m jv2zj + 兺 i⬎j j

.

If we divide the domain into a thin bins of thickness dz parallel to the x-y plane, V = LxLydz can be obtained where Lx and Ly are the lengths of the x-y plane, respectively. Equipartition of kinetic energy gives

冓

1 ␳共z兲兺 v2xj 2 i

冔冓 =

1 ␳共z兲兺 v2yj 2 i

冔冓 =

1 ␳共z兲兺 v2zj 2 i

冔

1 = kT具␳共z兲典V. 2

共27兲

Therefore, in terms of atomic position ri for atom i and potential ␾, the expression for surface tension is

␥=

冓

冔

x2ij + y 2ij − 2z2ij ␾⬘ij . 2Arij

兺兺 i⬎j j

冉冊

1 1 ⳵2 1 x2 = − , r8 6r6 24 ⳵x2 r4

␥LJ =

冓兺兺 j

i⬎j

␥C =

冓兺兺 i⬎j

j

x2ij

y 2ij

+ − 2Arij

2z2ij

␴ rij

1 q iq j 4␲⑀0⑀r rij

12

−

冔

␴ rij

6

⬘

−2

Notice that in Eq. 共30兲 i and j are only for oxygen atoms but in Eq. 共31兲, i and j are all the atoms in different molecules. By dropping the subscripts and using the nondimensional form, Eq. 共30兲 is

.

共34兲

2 exp共␤2r2兲 6 6 共␤ r + 3␤4r4 + 6␤2r2 + 6兲 Ar8 ⫻共x2 + y 2 − 2z2兲 −

24共x2 + y 2 − 2z2兲 . Ar14

共35兲

This term decays exponentially and can be calculated by the PP method with a small cutoff radius at 3.5␴. Evaluation of the long-range part of the Lennard-Jones surface tension is similar to the force evaluation except for p = 4. However, because it has the second partial derivative of 1 / r4, the three components should be multiplied by −k2x , −k2y , and −kz2, respectively, in Fourier space. For the Coulombic part of the surface tension, Ewald sums are used to calculate the pressure tensor, which is VP␣␤ =

1 兺 q jb 兺兺 qia兺 8␲⑀0 i a j⫽i b

冋冑

2

␲

⫻exp共− ␬2r2iajb兲 + erfc共␬riajb兲 +

冋冋

冉

+

␬riajb

册

共rij兲␣共riajb兲␤ r3iajb

1 2␲ 兺 Q共h兲S共h兲S共− h兲 4␲⑀0 V h⫽0

⫻ ␦␣␤ −

,

共31兲

⳵2 1 ⳵z2 r4

By using the expression given in Table I, for p = 4, the total short-range Lennard-Jones surface tension is

共30兲

.

冋冉冊冉冊冉冊册

12共x2 + y 2 − 2z2兲 1 ⳵2 1 ⳵2 1 = + Ar8 2A ⳵x2 r4 ⳵ y 2 r4

共28兲

冋冉冊冉冊册冔

共33兲

so that the first part of Eq. 共32兲 is given as

共29兲

x2ij + y 2ij − 2z2ij 4⑀LJ 2Arij

.

In Eq. 共32兲, the 1 / r14 term decays much faster than the 1 / r8 term. So the PPPM method is only applied on the 1 / r8 term. We notice that

Since the SPC/E potential can be described as the sum of two terms, the surface tension can also be expressed as

␥ = ␥LJ + ␥C ,

册冔

共32兲

␥LJ-SR =

and Pzz =

12共x2 + y 2 − 2z2兲 24共x2 + y 2 − 2z2兲 + Ar8 Ar14

共24兲

m jv2yj + 兺兺 ryij f yij , i⬎j

冓兺冋

共21兲

phase 1

Pxx =

␥LJ =

再

2h␣h␤ h␣h␤ − h2 2␬2

冊册

1 2␲ 兺兺共Pia兲␤qiah⫽0 兺 Q共h兲ih␣关S共h兲 4␲⑀0 V i a

冎

⫻exp共− ih · r兲 − S共− h兲exp共ih · r兲兴 .

共36兲

Here S共h兲 = 兺兺 qia exp共ih · ria兲,

共37兲

Q共h兲 = exp共− h2/4␬2兲/ប2 ,

共38兲

i

a


204715-5

J. Chem. Phys. 124, 204715 共2006兲

Surface tension of water by PPPM method TABLE II. Simulation parameters.

FIG. 1. Simulation domain configuration.

riajb = ria − r jb ,

共39兲

where ria is the position of atom a in the molecule i and r jb is the position of atom b in the molecule j. Reciprocal lattice vector h is defined as h = 2␲共l/Lx,m/Ly,n/Lz兲,

共40兲

where l, m, and n are the values 0, ±1 , ± 2 , . . . , ± ⬁, ␦␣␤ is the Kronecker delta, ␣ and ␤ are x, y, and z, pia = ria − ri is a vector between atom a in the molecule i and the center of mass of molecule i. V = LxLyLz is the volume and ␬ is the Ewald parameter which turns the relative weight of the real space and reciprocal space contribution. More detail can be found in Ref. 13. Once Pxx, Pyy, and Pzz have been obtained, the Coulombic surface tension can be calculated as

␥C =

1 2

冕

⬁

Parameters

Water

Temperature Density Time step

T = 300– 600 K ␳ = 1000 kg/ m3 ⌬t = 1.00⫻ 10−15 s

tances are constrained at 1.0 Å. The HOH bond angle is ␪0 = 109.47°. The intermolecular interactions consist of three point charges located at the oxygen and hydrogen sites, and the LJ potential used for O–O interactions. The equations of motion are solved using the Verlet algorithm with a time step of 1 ⫻ 10−15 s. The constraints inside the water molecules are handled using the SHAKE method.19 The simulations are carried out in an NVT ensemble. The system was equilibrated at a desired temperature for 100 000 time steps by velocity rescaling and 100 000 time steps of nonthermostat equilibration. Statistics were sampled 共density and surface tension兲 for an additional 100 000 time steps. The density and surface tension are calculated every ten steps. Simulation parameters are given in Table II. In using of PPPM method, the arbitrary positive number ␤ is chosen as 0.9 as suggested by Deserno and Holm.18,20 IV. RESULTS AND DISCUSSION A. Density

关PN共z兲 − PT共z兲兴dz,

共41兲

−⬁

where PN共z兲 = Pzz ,

共42兲

1 PT共z兲 = 关Pxx + Pyy兴. 2

共43兲

and

III. MODELLING AND SIMULATION

In order to simulate the liquid/vapor interface of water, a rectangular box is used for the simulation domain. Its dimensions are 32, 32, and 64 Å, respectively. The x, y, and z boundary conditions are set periodically. The initial condition is a liquid slab in the middle of the z direction in the simulation domain. 800 water molecules are placed in the liquid region using a cubic crystal structure as shown in Fig. 1. The initial velocities of molecules are set randomly corresponding to the system temperature. No external force field is applied. The cutoff radius for the short-range part is chosen as 0.98 nm. For the evaluation of the long-range part, 32⫻ 32 ⫻ 64 mesh points are used for both force and surface tension evaluation. In this simulation, the SPC/E interaction potential for water is used. Water 共H2O兲 molecule has a tetrahedral shape and the two hydrogen atoms are connected to the oxygen atom through covalent bonds. In this model, OH bond dis-

Once the system has reached equilibrium, the properties are calculated. The statistics start at 200 001 time step. Figure 2 is a snapshot of the equilibrium state of the simulation system at a temperature of 572 K. To obtain the density profile, the system was divided along the z axis into 64 cells, and the number of water molecules falling within each cell is counted. The density profile at different temperatures are shown in Fig. 3. It is clear from the profiles that a well defined interface is established and pure liquid and vapor regions exist on different sides. The calculated density profile is fitted to a hyperbolic tangent function of the form,

冉冊

1 1 z − z0 ␳共z兲 = 共␳l + ␳v兲 + 共␳l − ␳v兲tanh , 2 2 d

共44兲

where the fitting parameter ␳l is the bulk liquid density, ␳v the bulk vapor density, and z0 the middle of the interface between liquid and vapor. The average thickness d is found to be 1.68 Å at a temperature of 302 K. The distribution is

FIG. 2. Snapshot of the simulation domain at a temperature of 572 K.


204715-6

J. Chem. Phys. 124, 204715 共2006兲

Shi, Sinha, and Dhir

FIG. 3. Density profiles at various temperatures.

shown in Fig. 4. For the hyperbolic tangent function the “10–90 thickness,” t, is related to d by t = 2.1972d. The thickness can be obtained from above to be 3.91 Å. This is smaller than the values of 4.1 and 8.3 Å Matsumoto8 obtained from ellipsometric and x-ray reflectivity experiments. From Fig. 3, we can see that the thickness of the liquid region is more than 20 Å, which is wide enough to make an accurate estimate of the liquid and vapor density. In Fig. 5, the calculated results are compared with the data from experiments.21 Here we can see that at low temperatures, density values agree well with experimental results. In the region near the critical temperature, the density of the liquid from simulation is lower than that of the experimental data while vapor density is higher. This shows that the SPC/E model itself is not very accurate for the high temperature simulations. Currently, no potential model accurately simulates that region.

FIG. 5. Comparison with data from literature of the predicted liquid density and vapor density as a function of temperature.

Simulation results are compared with the thermodynamic correlations based on the equation of the corresponding states for the surface tension of water in Fig. 6. Filled symbols are the simulation results and the line is emprical fit to data from experiments.22 From the above figure it is clear that the simulation closely matches the thermodynamic correlation. This agreement shows that the surface tension should be calculated without truncating the long-rang part in the LJ potential interaction. By using the PPPM method, this error can by minimized and no artificial tail correction is needed at low temperatures. At the high temperature region, because of the low density of liquid water and high density of vapor, calculated surface tension differs from the experimental data. Another reason for the difference is, because of the long-range interactions, the two interfaces cannot be independent. To obtain more accurate results, a thicker and more stable film needs to be modeled.

B. Surface Tension V. CONCLUSION

Surface tension of the liquid-vapor interface is calculated using Eq. 共28兲. Since there are two liquid-vapor interfaces in this geometry, surface tension of the interface is one-half of the total summation value.

SPC/E water is simulated by using the PPPM method and the Ewald sum for the LJ part and the Coulomb part, respectively. The results show that in the low temperature

FIG. 4. Density profile with the tangent hyperbolic correlating function at 302 K.

FIG. 6. Comparison with data from literature of the surface tension as a function of temperature.


204715-7

region, density and surface tension agree with experimental values. In the high temperature region, the model for potential is not very accurate, as a result the density and surface tension differ from experimental data. The proposed method does not require any tail corrections. ACKNOWLEDGMENT

The authors gratefully acknowledge the support this research received from NASA under the Microgravity Fluid Physics Program. J. D. Bernal and R. H. Fowler, J. Chem. Phys. 1, 515 共1933兲. B. Guillot, J. Mol. Liq. 101, 219 共2002兲. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, and H. J. Hermans, in Intermolecular Forces, edited by Pullman 共Reidel, Dordrecht, 1981兲. 4 W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 共1983兲. 5 H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 共1987兲. 6 Y. Guissani and B. Guillot, J. Chem. Phys. 98, 8221 共1993兲. 7 J. J. de Pablo, J. M. Prausnitz, H. J. Strauch, and P. T. Cummings, J. Chem. Phys. 93, 7355 共1990兲. 8 M. Matsumoto and Y. Kataoka, J. Chem. Phys. 88, 3233 共1988兲. 1 2 3

J. Chem. Phys. 124, 204715 共2006兲

Surface tension of water by PPPM method 9

M. A. Wilson, A. Pohorille, and L. R. Pratt, J. Phys. Chem. 91, 4873 共1987兲. 10 G. C. Lie, S. Grigoras, L. X. Dang, D. Y. Yang, and D. A. McLean, J. Chem. Phys. 99, 3933 共1993兲. 11 S. B. Zhu, T. G. Fillingim, and G. W. Robinson, J. Phys. Chem. 95, 1002 共1991兲. 12 V. V. Zakharov, E. N. Brodskaya, and A. Laaksonen, J. Chem. Phys. 107, 10675 共1997兲. 13 J. Alejandre, D. J. Tildesley, and G. A. Chapela, J. Chem. Phys. 102, 4574 共1995兲. 14 S. Sinha, B. Shi, V. K. Dhir, J. Freund, and E. Darve, Proceedings of 2003 ASME Summer Heat Transfer Conferences, Las Vegas, 2003. 15 S. Sinha, V. K. Dhir, J. B. Freund, and E. Darve, J. Comput. Phys. 共submitted兲. 16 U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 共1995兲. 17 R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles 共IOP, Bristol, 1988兲. 18 M. Deserno and C. Holm, J. Chem. Phys. 109, 7678 共1998兲. 19 V. Krautler, W. Van Gunsteren, and P. H. Hunenberger, J. Comput. Chem. 22, 501 共2001兲. 20 M. Deserno and C. Holm, J. Chem. Phys. 109, 7694 共1998兲. 21 W. C. Reynolds, Thermodynamic Properties in SI 共Stanford University, Stanford, 1979兲. 22 A. F. Mills, Mass Transfer 共Prentice-Hall, Englewood Cliffs, NJ, 2001兲.
