- Email: [email protected]
A smooth truncation for polarizable water
Chemical Physics ELSEVIER
Chemical Physics 184 (1994) 163-169
A smooth truncation for polarizable water Bruce J. Palmer Analytic Sciences Department, Pacific Northwest Laboratory, Richland, WA99352, USA Received 16 March 1993; in final form 17 February 1994
Abstract
A systematic approach to smoothly truncatingthe Coulombinteractionin systems containing polarizable centers is presented and applied to the revised POLl model of water. The truncation leads to energy conserving trajectories without severly distorting the properties of the original model. Equilibrium properties and the diffusion constant are reported for the truncated POLl model.
1. Introduction
There has recently been a great deal of interest in developing polarizable water models for molecular dynamics simulations [l-5]. However, because the resulting induced dipoles represent a many-body interaction, the addition of polarizability greatly complicates the force calculation. The induced dipole moments associated with each of the polarizable centers are calculated via an iterative procedure and each iteration is equivalent to a complete force calculation for a nonpolarizable system. The result is that simulations of polarizable systems typically need much more CPU time than comparable simulations on unpolarizable models. To keep CPU requirements at reasonable levels, the use of a cutoff is highly desirable. Furthermore, the cutoff should lead to energy conserving trajectories so that both equilibrium and nonequilibrium properties can be obtained from the simulation. However, the use of cutoffs in polarizable systems has not been extensively discussed in the literature. Several authors have used molecule-centered cutoffs in simulations [ 1,5,6]. Dang has compared the use of a sharp molecule-centered cutoff in polarizable systems to an Ewald sum and found that the truncation gives results that are very 0301-0104/94/$07.00 0 1994 Elsevier Science B.V. All rights reserved .SSD10301-0104(94)00088-R
close to the Ewald sum. Van Belle et al. [7] investigated the use of a sharp cutoff for the induced dipole interactions, but used the full Coulomb interaction for the static charges. They found that the energy due to the induced dipoles was relatively insensitive to the cutoff distance over a wide range of values. The molecule-centered truncations appear to work well for water, but molecule-centered truncation schemes are not easily generalized to larger molecules. This paper will discuss a site-site truncation of both the static charge and induced dipole interactions in polarizable models. This truncation works fairly well and can be applied to molecules of arbitrary size. The truncation is applied to the POLl water model recently proposed by Dang [sl. Many of the properties of pure water can be reproduced using models composed of short-ranged interactions and fixed partial charges. The solvent induced changes in the dipole are incorporated in the model, in a mean field sense, by using parameters that have been optimized to reproduce the properties of the liquid. These models usually have permanent dipole moments that are substantially higher than the gas phase value. For homogeneous systems, where the average electric field vanishes, this is expected to be a good approximation. For dissolved ions, however, the assumption
164
B.J. Palmer/Chemical Physics 184 (1994) 163-169
that all water molecules experience the same environment breaks down. It is clear that water molecules close to the ion will be subjected to stronger electric fields than molecules farther away. The polarizability of water implies that the average dipole moment of the water molecules changes as a function of the ion-water distance. The accurate reproduction of this behavior requires polarizable models. Several polarizable models of water have appeared in the last few years. Besides the polarizable water model of Dang [ 51 discussed in this paper, polarizable versions of SPC [ 1 ] and TIP4P [ 21 water have been developed. Like SPC water, the polarizable versions of SPC [ 1 ] and TIP4P [ 21 water have been developed. Like SPC water, the polarizable SPC model (PSPC) gives relatively poor pair distributions. However the overall energetics and average dipole moment are fairly good [ 11. The polarizable TIP4P model of Sprik and Klein [ 21 gives pair distribution functions that arecomparable in quality to the standard TIP4P model and a value for the energy that is within a few percent of experimental results. The average dipole moment for polarizable TIWP water is in good agreement with experimental results although it is about 10% too large. The diffusion constant for polarizable TIP4P is ( 1.5 *OS) X lo-’ cm’/s [ 81, and is substantially lower than the experimental value of 2.4X 10m5 cm*/s [ 91. In contrast, the calculations presented here indicate that the truncated version of Dang’s polarizable water model is capable of accurately reproducing almost all the properties of water under standard conditions.
to produce continuous forces due to the induced dipole moments, both the first and second derivatives of the truncated Coulomb potential must be continuous. Therefore, more parameters are required in the truncation. For these reasons, the second approach of applying the truncations term by term is adopted in this paper. The formula for the electrostatic interactions in a system containing fixed charges and polarizable centers is [ 1,101
(2.1) where qi is the charge at site i, pi is the induced dipole at site i, and Ei is the static electric field at site i due to all the other fixed charges in the system. The electric fields Ei and the induced dipole moments pi are given by the equations Ei=
Ciru, j#i
(2.2)
rij
CT.*pj),
P.=ai(Er+
(2.3)
r#j
where (Yiis the polarizability associated with site i and T, is the dipole tensor (2.4) The truncation scheme used in this paper is to separately truncate the Coulomb interaction and the factor of rw3 appearing in both the electric field and dipole tensor expressions. The Coulomb interaction is smoothly truncated by replacing ri; ’ in
2. Truncating the potential
(2.5) There are two approaches to truncating the potential for a polarizable system. The first is to apply a truncation to the original Coulomb interaction and then derive expressions for the electric field, dipole tensor, and induced dipole moments based on the modified Coulomb interaction. The second approach is to start with the equations for the full Coulomb potential, electric field, and induced dipole moments and apply truncations to these terms separately. The first approach has the virtue of consistency, but errors introduced in the Coulomb interaction are magnified in the electric field and induced dipole calculations. Furthermore, in order
with the shifted force expression qb(,)=
_.L -t
+ -$(rti--rc),
rij
r,j
=
0,
rti> r,,
(2.6)
where r, is the cutoff distance. The shifted force truncation has been shown to give very good results for standard water models [ 111. The electric field interactions and the dipole tensor are truncated by replacing r V-3 with the function
B.J. Palmer/Chemical
165
Physics 184 (1994) 16.3-169
From Eq. (2.12) the gradient vector does not act on the implicit r dependence of the pi. The force on site i due to this term is then
= 0,
rii > r, .
(2.7)
(M*(VnE))i=
C [GJl’(ri/)rc(ri,*pi)
+qj$(cj)Pi
j#i
These two replacements are enough to guarantee that all interactions vanish for r > r, and that both the potential energy and the forces are continuous functions at r,. To calculate the forces and the contribution to the pressure due to the polarization energy
where
Upi= - i Cpi.Ei
Similarly, the force on site i from the dipole tensor term in Eq. (2.13) is
(2.8)
I
it is convenient to write the coordinates, induced dipoles, and the local electric fields as large single vectors, R , M , E , respectively, of dimension 3N. Similarly, the dipole tensor and polarizabilities can be thought of as 3N X 3N tensors, T and A. Using this notation, the equation for the induced dipoles can be written as M=A.[E+T*M]
,
(2.9)
-4iG’(rij)rij(rijaPj)
q(r)
=
=
+‘( r&
(2.10)
The polarization energy can then be written as C&i=--$E.M=-tE*[I-A-T]-‘.A.E.
5
[ II
--Jl(rii) : rl
(2.14)
@i l rij) @j “ij 1 -Pi ‘Pj
1
[$rij@i*rij)(p,r,) u
-Pi(pi’rij)-@i’rij)Pj
which can be formally solved for M to get M=[I-A-T]-‘*A-E.
9
d@(r)l --=(-;+;)s. dr r
xi
j#i
-4i9Nrij)Pjl
II
.
(2.15)
It remains to calculate the contribution to the pressure from the polarization energy. The pressure is given by the expression [ 12,131
(2.11)
Taking the gradient of U,, with respect to R yields -V,U,,=E.A[I-A*T]-‘.(VnE) +$E*A[I-T.A]-‘.(V,T)[I-A-T]-‘*A-E =M*(VuE)
(2.12)
+&M.(V,T)*M.
This equation follows immediately from the relations A[I-T-A]-‘=[I-A.T]-‘.A and VRB-‘=
-B-‘.(V,B)nB-’
,
where B is an arbitrary, spatially-dependent matrix. Returning to the single particle notation, the M-E term can be written as CPi*Ei= i
C C ~@(r~)~ij*Pi. i j#i
(2.13)
where V is the volume and NT is the number of independently translating centers. For systems without any rigid internal constraints, NT is equal to the total number of atoms in the system. For rigid molecule systems, NT is equal to the total number of molecules. The factor of - 1 is for periodic systems where the center-of-mass translational motion is conserved. For the rigid water model considered in this paper, care must be taken in evaluating the derivative with respect to volume [ 13,141. Ordinarily, the volume derivative is done by writing the particle coordinates as r.-Vi’3s. I11 where Si is an unscaled coordinate. If V”3~i is substituted for ri in the potential energy, then the expression for NJ/W follows immediately. For rigid molecules, this formulation is not valid because it is no longer
166
B.J. Palmer/Chemical Physics 184 (I 994) 163-l 69
meaningful to contract all distances uniformly. While it is still possible to translate molecules as a whole towards each other, the internal distances must be held fixed [ 13,151. To calculate tire volume derivatives in this case, the location of each of the interaction sites is written as
where ri, is the location of a reference site and di’i is a vector specifying the displacement of the site at ri relative to the reference site rij. Each rigid unit has one reference site associated with it so several sites ri will have the same reference site riS. Using unscaled coordinates, the ri can be written as ri=v “3Si, +d. 1’1 .*
(2.16)
For rigid systems, the displacements diri remain fixed as the reference sites are translated relative to one another. Using Eq. (2.16)) the derivative of a scalar function #I of rii with respect to the volume is
From this it is clear that tire pressure calculation is almost identical to the force calculation and gives
using an 8.5 A cutoff in a 18.64 A cubic cell. The water model used in the revised POLl model recently proposed by Dang [ 51. The portion of the potential that is independent of the induced dipoles, Upair,has the form
The Lennard-Jones terms
are located only on the oxygen atoms while the partial charges are located on all atoms. Each atom also has a finite polarizability associated with it. A complete list of the parameters for this model is given in Table 1. The Lennard-Jones potentials are truncated by replacing the standard Lennard-Jones potential with
+ai,+b&r, -ru>,
ril
=o,
rii > r, ,
and fixing aij and b, by the requirement that both the potential and its first derivative with respect to rv vanish at the cutoff. The expressions for au and b, are .“=-4qFr b,j=-
)(pj.ri7,)
)I
, (2.17)
3. Simulation details Using the truncation procedure described above, simulations of 216 water molecules were performed
_ (:T],
$[+$6(~~].
The POLl model was originally designed to reproduce the 1 atmosphere pressure of water at 0.97 g/cm3 and 298 K. THe pressure is extremely sensitive to details of the potential, and the density decreased about 4% when the original POLI potential was truncated. The model was reparameterized for the truncated potential to give approximately a 1 atmosphere presTable 1 Potential pammeters for the truncated FOLl model. The molecular geometry is ROM= I .OA and LHOH = 109.5 Atom type
u(A)
t (kcal/mol)
q
ff (A?
oxygen
3.196 0.000
0.147 0.000
- 0.730 + 0.365
0.528 0.170
hydrogen
B.J. Palmer/Chemical
sure under standard conditions. This was accomplished by adjusting the well depth of the Lennard-Jones term on the oxygen atom, which required only a modest change in the well depth from 0.160 to 0.147 kcal/mol. The pressure calculation was verified by running simulations using the Andersen-Nod algorithm for constant temperature-constant pressure dynamics and checking that the Hamiltonian associated with this algorithm was conserved [ 15-171. The pressure fluctuation are typically very large, for the simulations reported here the root mean square pressure fluctuations were about 800 atmospheres. This made it difficult to determine the pressure for the model exactly, but for these parameters the model pressure is probably within 50 atmospheres of 1 atmosphere at the experimental density. A series of six simulations were run and the averages collected. Each simulation was 25.0 ps long and consisted of 10000 steps using a 2.5 fs timestep. The temperature was maintained at 298 K using a Nose thermostat with the temperature mass set at 5 X lo3 amu A’. The time constant for oscillations in the temperature scaling variable was about 1 ps for this value of the temperature mass. The equations of motion were integrated using the velocity Verlet algorithm recast as a predictor-corrector [ 181 and the rigid water geometry was maintained using a variant of the SHAKE algorithm [ 19,203. For these simulation conditions, the ratio of the root mean square fluctuations in the Nose Hamiltonian to the average kinetic energy was 1 X 10W3,indicating excellent Hamiltonian conservation.
Table 2 Equilibrium properties of the truncated POLl model at 298 K and p=O.997 g/cm’
prvm
Simulation
(II) (kcai/mol) (U,,& (kcal/mol) P (bar)
- 9.38 f 0.03 -9.90*0.09 -2.84f0.01 - 2.96 -30*55 1.0 2.592 f 0.002 2.61 f 0.01 22.1 f 5.9
M (debye) C, (kcal/mol K)
Experiment
Dang ’
- 9.92 b
1 2.6’ 17.8 b
’ Ref. [S]. Dang’s results are for 303 K. b Ref. [21]. ’ Ref. [22].
The value obtained from the simulations is relatively high compared to experimental results [ 211. However, Cv is sensitive to the way the intramolecular vibrations are handled and also has substantial quantum corrections [ 231. It is therefore problematic whether exact agreement can be obtained using a rigid model and classical dynamics. The pair distribution functions for the polarizable water model are shown in Figs. l-3, along with the experimental curves of Soper and Phillips [ 241. The oxygen-xygen pair distribution function is almost identical to the curve reported by Dang, and also agrees quite well with experimental results. The main peak in the oxygen+xygen pair distribution function for the polarizable model is shifted about 0.15 A from the experimental curve and is slightly lower than experiment. The secondary peaks and minima are also more c
4. Results The equilibrium properties calculated from the simulations are listed in Table 2. The average potential energy is in reasonable agreement with experiment [ 211 although the agreement is not as close as for the original POLl model. Some additional reparameterization of the truncated potential will be necessary to get complete agreement with experiment. Like the original POLl model, the average dipole moment of the truncated model matches the experimental value almost exactly [ 221. The constant volume heat capacity Cv can be calculated from fluctuations in the energy using the relation [ 181
167
Physics 184 (1994) 163-169
-I
:
0.0 0
2
I
I
.
I
4 8 6 Distance (A) Fig. 1. Oxygen-oxygen pair distribution function. ( -_) work. (---) Experimental curve of Soper and Phillips.
a
10 This
B.J. Palmer/Chemical
168
Physics 184 (1994)
3.5
grated out to the first minimum at 2.4 A, they both give 1.83 hydrogen atoms in the first peak, so the differences between the simulation and experiment apparently cancel out. The diffusion constant D was calculated from the displacement formula
3.0 2.5
s
2.0
-
9
1.5
-
1.0
-
0.5
-
D= ;;
1
II
0.0 0
163-169
2
I
I
I
4 6 Distance (A)
I
I 8
Fig. 2. Oxygen-hydrogen pair distribution function. (work. (- - -) Experimental curve of Soper and Phillips.
I 10 )ThiS
3.5 3.0
([r(r)
-r(O)]‘).
The mean square displacement of the oxygen atoms was calculated as a function of time using multiple time origins. The diffusion constant was calculated from the slope of a linear least squares fit to ([r(t) - r( 0) ] ‘) on the interval between 2.0 and 4.0 ps. The calculated diffusion constant is (2.50f0.16) X 10e5 cm’/s, which reproduces the experimental value of 2.4 X 10m5 cm2/s [9] almost exactly.
2.5
5. Conclusions
0.0
3 0
2
8
10
Ltallce($ Fig. 3. Hydrogen-hydrogen pair distribution function. (work. (-- -) Experimental curve of Soper and Phillips.
) This
pronounced than experiment. The secondary peaks and minima are also more pronounced in the experimental distribution function. A very small feature in the theoretical curve can be seen at the cutoff distance of 8.5 8, and is the result of the truncation. However, the small size of this feature, plus the fact that it is nearly unobservable in the oxygenAydrogen and hydrogenhydrogen distribution functions, indicates that the cutoff is not seriously disrupting the local structure. The remaining pair distribution functions also agree closely with experimental results. The single hugest difference is the hydrogen bonding peak in the oxygenhydrogen pair distribution function. Part of the discrepancy can be accounted for by the choice of parameters for the rigid water geometry, which is slightly different than the experimental geometry of liquid water. However, if the experimental and simulation curves are inte-
A scheme for systematically truncating polarizable potentials has been presented and applied to the revised POLl model of water. The truncated potential produces trajectories that conserve energy to the same degree as nonpolarizable models. The truncation does not produce any large artifacts in the pair distribution functions at the cutoff. For POLl water, the truncation does not significantly distort the properties of the original model. Preliminary simulations on 256 molecules indicate that most properties of the model are insensitive to system size. The only exception is the pressure, which may increase to lOG200 atmospheres. The availability of smooth truncation schemes should serve to increase both the range of polarizable models that can be conveniently simulated and the number of properties, particularly dynamical quantities, that can be calculated from simulations.
Acknowledgement
The author is indebted to Liem Dang and Bruce Garrett for many useful discussions. This work was performed under the auspices of the Division of Chemical Sciences, Office of Basic Energy Sciences, U.S. Department of Energy. Pacific Northwest Laboratory
B.J. Palmer/Chemical
is operated for the U.S. Department of Energy by Battelle Memorial Institute under Contract DE-ACO676RL0 1830.
References [ 11 P. Ahlstr6m. A. Wallqvist, S. Engstriim and B. Jiinsson. Mol. Phys. 68 (1989) 563. [ 21 M. Sprik and M.L. Klein, J. Chem. Phys. 89 ( 1988) 7556. [3] A. Wallqvist, Chem. Phys. 148 (1990) 439. [4] L.X. Dang, J. Chem. Phys. 96 ( 1992) 6970. [5] L.X. Dang, J. Chem. Phys. 97 ( 1992) 2659. [6] P. Barnes, J.L. Finney, J.D. Nicholas and J.E. Quinn, Nature 282 (1979) 459. [7] D. Van Belle, 1. Couplet, M. Prevost and S.J. Wodak. J. Mol. Biol. 198 (1987) 721. [8] M. Sprik. M.L. Klein and K. Watanabe, J. Phys. Chem. 94 (1990) 6483. [9] K. Krynicki, C.D. Gteen and D.W. Sawyer, Faraday Discussions Chem. Sot. 66 ( 1978) 199.
Physics 184 (1994) 163469
169
[lo] F.J. Vesely, J. Comput. Phys. 24 (1977) 361. [ 111 C.L. Brooks, B.M. Pettitt and M. Karplus, J. Chem. Phys. 83 (1985) 5897. [ 121 D.A. McQuarrie, Statistical mechanics (Harper and Row, New York, 1976). [ 131 W. Smith, CCP5 Info. Quart. 26 (1987) 43. [ 141 M. Ferrario and J.P. Ryckaert, Mol. Phys. 54 ( 1985) 587. [ 151 G. Ciccotti, M. Ferrario and J.P. Ryckaert. Mol. Phys. 47 (1982) 1253. [ 161 H.C. Andersen, J. Chem. Phys. 72 (1980) 2384. [17] S.Nod.J.Chem.Phys.81 (1984) 511. [ 181 M.P. Allen and D.J. Tildesley, Computer simulation of liquids (Clarendon Press, Oxford, 1987). [ 191 B.J. Palmer, J. Comput. Phys. 104 (1993) 470. [20] B.J. Palmer and B.C. Garrett, J. Chem. Phys. 98 (1993) 4047. [ 2 11 D. Eisenberg and W. Kauzmann, The structure and properties of water (Oxford Univ. Press, Oxford, 1969). [22] C.A. Coulson and D. Eisenberg, Proc. Roy. Sot. A 291 ( 1966) 445. [23] P.H. Betens, D.H.J. Mackay, G.M. White and K.R. Wilson, J. Chem. Phys. 79 (1983) 2375. [24] A.K. Soper and M.G. Phillips, Chem. Phys. 107 ( 1986) 47.
Chemical Physics 184 (1994) 163-169
A smooth truncation for polarizable water Bruce J. Palmer Analytic Sciences Department, Pacific Northwest Laboratory, Richland, WA99352, USA Received 16 March 1993; in final form 17 February 1994
Abstract
A systematic approach to smoothly truncatingthe Coulombinteractionin systems containing polarizable centers is presented and applied to the revised POLl model of water. The truncation leads to energy conserving trajectories without severly distorting the properties of the original model. Equilibrium properties and the diffusion constant are reported for the truncated POLl model.
1. Introduction
There has recently been a great deal of interest in developing polarizable water models for molecular dynamics simulations [l-5]. However, because the resulting induced dipoles represent a many-body interaction, the addition of polarizability greatly complicates the force calculation. The induced dipole moments associated with each of the polarizable centers are calculated via an iterative procedure and each iteration is equivalent to a complete force calculation for a nonpolarizable system. The result is that simulations of polarizable systems typically need much more CPU time than comparable simulations on unpolarizable models. To keep CPU requirements at reasonable levels, the use of a cutoff is highly desirable. Furthermore, the cutoff should lead to energy conserving trajectories so that both equilibrium and nonequilibrium properties can be obtained from the simulation. However, the use of cutoffs in polarizable systems has not been extensively discussed in the literature. Several authors have used molecule-centered cutoffs in simulations [ 1,5,6]. Dang has compared the use of a sharp molecule-centered cutoff in polarizable systems to an Ewald sum and found that the truncation gives results that are very 0301-0104/94/$07.00 0 1994 Elsevier Science B.V. All rights reserved .SSD10301-0104(94)00088-R
close to the Ewald sum. Van Belle et al. [7] investigated the use of a sharp cutoff for the induced dipole interactions, but used the full Coulomb interaction for the static charges. They found that the energy due to the induced dipoles was relatively insensitive to the cutoff distance over a wide range of values. The molecule-centered truncations appear to work well for water, but molecule-centered truncation schemes are not easily generalized to larger molecules. This paper will discuss a site-site truncation of both the static charge and induced dipole interactions in polarizable models. This truncation works fairly well and can be applied to molecules of arbitrary size. The truncation is applied to the POLl water model recently proposed by Dang [sl. Many of the properties of pure water can be reproduced using models composed of short-ranged interactions and fixed partial charges. The solvent induced changes in the dipole are incorporated in the model, in a mean field sense, by using parameters that have been optimized to reproduce the properties of the liquid. These models usually have permanent dipole moments that are substantially higher than the gas phase value. For homogeneous systems, where the average electric field vanishes, this is expected to be a good approximation. For dissolved ions, however, the assumption
164
B.J. Palmer/Chemical Physics 184 (1994) 163-169
that all water molecules experience the same environment breaks down. It is clear that water molecules close to the ion will be subjected to stronger electric fields than molecules farther away. The polarizability of water implies that the average dipole moment of the water molecules changes as a function of the ion-water distance. The accurate reproduction of this behavior requires polarizable models. Several polarizable models of water have appeared in the last few years. Besides the polarizable water model of Dang [ 51 discussed in this paper, polarizable versions of SPC [ 1 ] and TIP4P [ 21 water have been developed. Like SPC water, the polarizable versions of SPC [ 1 ] and TIP4P [ 21 water have been developed. Like SPC water, the polarizable SPC model (PSPC) gives relatively poor pair distributions. However the overall energetics and average dipole moment are fairly good [ 11. The polarizable TIP4P model of Sprik and Klein [ 21 gives pair distribution functions that arecomparable in quality to the standard TIP4P model and a value for the energy that is within a few percent of experimental results. The average dipole moment for polarizable TIWP water is in good agreement with experimental results although it is about 10% too large. The diffusion constant for polarizable TIP4P is ( 1.5 *OS) X lo-’ cm’/s [ 81, and is substantially lower than the experimental value of 2.4X 10m5 cm*/s [ 91. In contrast, the calculations presented here indicate that the truncated version of Dang’s polarizable water model is capable of accurately reproducing almost all the properties of water under standard conditions.
to produce continuous forces due to the induced dipole moments, both the first and second derivatives of the truncated Coulomb potential must be continuous. Therefore, more parameters are required in the truncation. For these reasons, the second approach of applying the truncations term by term is adopted in this paper. The formula for the electrostatic interactions in a system containing fixed charges and polarizable centers is [ 1,101
(2.1) where qi is the charge at site i, pi is the induced dipole at site i, and Ei is the static electric field at site i due to all the other fixed charges in the system. The electric fields Ei and the induced dipole moments pi are given by the equations Ei=
Ciru, j#i
(2.2)
rij
CT.*pj),
P.=ai(Er+
(2.3)
r#j
where (Yiis the polarizability associated with site i and T, is the dipole tensor (2.4) The truncation scheme used in this paper is to separately truncate the Coulomb interaction and the factor of rw3 appearing in both the electric field and dipole tensor expressions. The Coulomb interaction is smoothly truncated by replacing ri; ’ in
2. Truncating the potential
(2.5) There are two approaches to truncating the potential for a polarizable system. The first is to apply a truncation to the original Coulomb interaction and then derive expressions for the electric field, dipole tensor, and induced dipole moments based on the modified Coulomb interaction. The second approach is to start with the equations for the full Coulomb potential, electric field, and induced dipole moments and apply truncations to these terms separately. The first approach has the virtue of consistency, but errors introduced in the Coulomb interaction are magnified in the electric field and induced dipole calculations. Furthermore, in order
with the shifted force expression qb(,)=
_.L -t
+ -$(rti--rc),
rij
r,j
=
0,
rti> r,,
(2.6)
where r, is the cutoff distance. The shifted force truncation has been shown to give very good results for standard water models [ 111. The electric field interactions and the dipole tensor are truncated by replacing r V-3 with the function
B.J. Palmer/Chemical
165
Physics 184 (1994) 16.3-169
From Eq. (2.12) the gradient vector does not act on the implicit r dependence of the pi. The force on site i due to this term is then
= 0,
rii > r, .
(2.7)
(M*(VnE))i=
C [GJl’(ri/)rc(ri,*pi)
+qj$(cj)Pi
j#i
These two replacements are enough to guarantee that all interactions vanish for r > r, and that both the potential energy and the forces are continuous functions at r,. To calculate the forces and the contribution to the pressure due to the polarization energy
where
Upi= - i Cpi.Ei
Similarly, the force on site i from the dipole tensor term in Eq. (2.13) is
(2.8)
I
it is convenient to write the coordinates, induced dipoles, and the local electric fields as large single vectors, R , M , E , respectively, of dimension 3N. Similarly, the dipole tensor and polarizabilities can be thought of as 3N X 3N tensors, T and A. Using this notation, the equation for the induced dipoles can be written as M=A.[E+T*M]
,
(2.9)
-4iG’(rij)rij(rijaPj)
q(r)
=
=
+‘( r&
(2.10)
The polarization energy can then be written as C&i=--$E.M=-tE*[I-A-T]-‘.A.E.
5
[ II
--Jl(rii) : rl
(2.14)
@i l rij) @j “ij 1 -Pi ‘Pj
1
[$rij@i*rij)(p,r,) u
-Pi(pi’rij)-@i’rij)Pj
which can be formally solved for M to get M=[I-A-T]-‘*A-E.
9
d@(r)l --=(-;+;)s. dr r
xi
j#i
-4i9Nrij)Pjl
II
.
(2.15)
It remains to calculate the contribution to the pressure from the polarization energy. The pressure is given by the expression [ 12,131
(2.11)
Taking the gradient of U,, with respect to R yields -V,U,,=E.A[I-A*T]-‘.(VnE) +$E*A[I-T.A]-‘.(V,T)[I-A-T]-‘*A-E =M*(VuE)
(2.12)
+&M.(V,T)*M.
This equation follows immediately from the relations A[I-T-A]-‘=[I-A.T]-‘.A and VRB-‘=
-B-‘.(V,B)nB-’
,
where B is an arbitrary, spatially-dependent matrix. Returning to the single particle notation, the M-E term can be written as CPi*Ei= i
C C ~@(r~)~ij*Pi. i j#i
(2.13)
where V is the volume and NT is the number of independently translating centers. For systems without any rigid internal constraints, NT is equal to the total number of atoms in the system. For rigid molecule systems, NT is equal to the total number of molecules. The factor of - 1 is for periodic systems where the center-of-mass translational motion is conserved. For the rigid water model considered in this paper, care must be taken in evaluating the derivative with respect to volume [ 13,141. Ordinarily, the volume derivative is done by writing the particle coordinates as r.-Vi’3s. I11 where Si is an unscaled coordinate. If V”3~i is substituted for ri in the potential energy, then the expression for NJ/W follows immediately. For rigid molecules, this formulation is not valid because it is no longer
166
B.J. Palmer/Chemical Physics 184 (I 994) 163-l 69
meaningful to contract all distances uniformly. While it is still possible to translate molecules as a whole towards each other, the internal distances must be held fixed [ 13,151. To calculate tire volume derivatives in this case, the location of each of the interaction sites is written as
where ri, is the location of a reference site and di’i is a vector specifying the displacement of the site at ri relative to the reference site rij. Each rigid unit has one reference site associated with it so several sites ri will have the same reference site riS. Using unscaled coordinates, the ri can be written as ri=v “3Si, +d. 1’1 .*
(2.16)
For rigid systems, the displacements diri remain fixed as the reference sites are translated relative to one another. Using Eq. (2.16)) the derivative of a scalar function #I of rii with respect to the volume is
From this it is clear that tire pressure calculation is almost identical to the force calculation and gives
using an 8.5 A cutoff in a 18.64 A cubic cell. The water model used in the revised POLl model recently proposed by Dang [ 51. The portion of the potential that is independent of the induced dipoles, Upair,has the form
The Lennard-Jones terms
are located only on the oxygen atoms while the partial charges are located on all atoms. Each atom also has a finite polarizability associated with it. A complete list of the parameters for this model is given in Table 1. The Lennard-Jones potentials are truncated by replacing the standard Lennard-Jones potential with
+ai,+b&r, -ru>,
ril
=o,
rii > r, ,
and fixing aij and b, by the requirement that both the potential and its first derivative with respect to rv vanish at the cutoff. The expressions for au and b, are .“=-4qFr b,j=-
)(pj.ri7,)
)I
, (2.17)
3. Simulation details Using the truncation procedure described above, simulations of 216 water molecules were performed
_ (:T],
$[+$6(~~].
The POLl model was originally designed to reproduce the 1 atmosphere pressure of water at 0.97 g/cm3 and 298 K. THe pressure is extremely sensitive to details of the potential, and the density decreased about 4% when the original POLI potential was truncated. The model was reparameterized for the truncated potential to give approximately a 1 atmosphere presTable 1 Potential pammeters for the truncated FOLl model. The molecular geometry is ROM= I .OA and LHOH = 109.5 Atom type
u(A)
t (kcal/mol)
q
ff (A?
oxygen
3.196 0.000
0.147 0.000
- 0.730 + 0.365
0.528 0.170
hydrogen
B.J. Palmer/Chemical
sure under standard conditions. This was accomplished by adjusting the well depth of the Lennard-Jones term on the oxygen atom, which required only a modest change in the well depth from 0.160 to 0.147 kcal/mol. The pressure calculation was verified by running simulations using the Andersen-Nod algorithm for constant temperature-constant pressure dynamics and checking that the Hamiltonian associated with this algorithm was conserved [ 15-171. The pressure fluctuation are typically very large, for the simulations reported here the root mean square pressure fluctuations were about 800 atmospheres. This made it difficult to determine the pressure for the model exactly, but for these parameters the model pressure is probably within 50 atmospheres of 1 atmosphere at the experimental density. A series of six simulations were run and the averages collected. Each simulation was 25.0 ps long and consisted of 10000 steps using a 2.5 fs timestep. The temperature was maintained at 298 K using a Nose thermostat with the temperature mass set at 5 X lo3 amu A’. The time constant for oscillations in the temperature scaling variable was about 1 ps for this value of the temperature mass. The equations of motion were integrated using the velocity Verlet algorithm recast as a predictor-corrector [ 181 and the rigid water geometry was maintained using a variant of the SHAKE algorithm [ 19,203. For these simulation conditions, the ratio of the root mean square fluctuations in the Nose Hamiltonian to the average kinetic energy was 1 X 10W3,indicating excellent Hamiltonian conservation.
Table 2 Equilibrium properties of the truncated POLl model at 298 K and p=O.997 g/cm’
prvm
Simulation
(II) (kcai/mol) (U,,& (kcal/mol) P (bar)
- 9.38 f 0.03 -9.90*0.09 -2.84f0.01 - 2.96 -30*55 1.0 2.592 f 0.002 2.61 f 0.01 22.1 f 5.9
M (debye) C, (kcal/mol K)
Experiment
Dang ’
- 9.92 b
1 2.6’ 17.8 b
’ Ref. [S]. Dang’s results are for 303 K. b Ref. [21]. ’ Ref. [22].
The value obtained from the simulations is relatively high compared to experimental results [ 211. However, Cv is sensitive to the way the intramolecular vibrations are handled and also has substantial quantum corrections [ 231. It is therefore problematic whether exact agreement can be obtained using a rigid model and classical dynamics. The pair distribution functions for the polarizable water model are shown in Figs. l-3, along with the experimental curves of Soper and Phillips [ 241. The oxygen-xygen pair distribution function is almost identical to the curve reported by Dang, and also agrees quite well with experimental results. The main peak in the oxygen+xygen pair distribution function for the polarizable model is shifted about 0.15 A from the experimental curve and is slightly lower than experiment. The secondary peaks and minima are also more c
4. Results The equilibrium properties calculated from the simulations are listed in Table 2. The average potential energy is in reasonable agreement with experiment [ 211 although the agreement is not as close as for the original POLl model. Some additional reparameterization of the truncated potential will be necessary to get complete agreement with experiment. Like the original POLl model, the average dipole moment of the truncated model matches the experimental value almost exactly [ 221. The constant volume heat capacity Cv can be calculated from fluctuations in the energy using the relation [ 181
167
Physics 184 (1994) 163-169
-I
:
0.0 0
2
I
I
.
I
4 8 6 Distance (A) Fig. 1. Oxygen-oxygen pair distribution function. ( -_) work. (---) Experimental curve of Soper and Phillips.
a
10 This
B.J. Palmer/Chemical
168
Physics 184 (1994)
3.5
grated out to the first minimum at 2.4 A, they both give 1.83 hydrogen atoms in the first peak, so the differences between the simulation and experiment apparently cancel out. The diffusion constant D was calculated from the displacement formula
3.0 2.5
s
2.0
-
9
1.5
-
1.0
-
0.5
-
D= ;;
1
II
0.0 0
163-169
2
I
I
I
4 6 Distance (A)
I
I 8
Fig. 2. Oxygen-hydrogen pair distribution function. (work. (- - -) Experimental curve of Soper and Phillips.
I 10 )ThiS
3.5 3.0
([r(r)
-r(O)]‘).
The mean square displacement of the oxygen atoms was calculated as a function of time using multiple time origins. The diffusion constant was calculated from the slope of a linear least squares fit to ([r(t) - r( 0) ] ‘) on the interval between 2.0 and 4.0 ps. The calculated diffusion constant is (2.50f0.16) X 10e5 cm’/s, which reproduces the experimental value of 2.4 X 10m5 cm2/s [9] almost exactly.
2.5
5. Conclusions
0.0
3 0
2
8
10
Ltallce($ Fig. 3. Hydrogen-hydrogen pair distribution function. (work. (-- -) Experimental curve of Soper and Phillips.
) This
pronounced than experiment. The secondary peaks and minima are also more pronounced in the experimental distribution function. A very small feature in the theoretical curve can be seen at the cutoff distance of 8.5 8, and is the result of the truncation. However, the small size of this feature, plus the fact that it is nearly unobservable in the oxygenAydrogen and hydrogenhydrogen distribution functions, indicates that the cutoff is not seriously disrupting the local structure. The remaining pair distribution functions also agree closely with experimental results. The single hugest difference is the hydrogen bonding peak in the oxygenhydrogen pair distribution function. Part of the discrepancy can be accounted for by the choice of parameters for the rigid water geometry, which is slightly different than the experimental geometry of liquid water. However, if the experimental and simulation curves are inte-
A scheme for systematically truncating polarizable potentials has been presented and applied to the revised POLl model of water. The truncated potential produces trajectories that conserve energy to the same degree as nonpolarizable models. The truncation does not produce any large artifacts in the pair distribution functions at the cutoff. For POLl water, the truncation does not significantly distort the properties of the original model. Preliminary simulations on 256 molecules indicate that most properties of the model are insensitive to system size. The only exception is the pressure, which may increase to lOG200 atmospheres. The availability of smooth truncation schemes should serve to increase both the range of polarizable models that can be conveniently simulated and the number of properties, particularly dynamical quantities, that can be calculated from simulations.
Acknowledgement
The author is indebted to Liem Dang and Bruce Garrett for many useful discussions. This work was performed under the auspices of the Division of Chemical Sciences, Office of Basic Energy Sciences, U.S. Department of Energy. Pacific Northwest Laboratory
B.J. Palmer/Chemical
is operated for the U.S. Department of Energy by Battelle Memorial Institute under Contract DE-ACO676RL0 1830.
References [ 11 P. Ahlstr6m. A. Wallqvist, S. Engstriim and B. Jiinsson. Mol. Phys. 68 (1989) 563. [ 21 M. Sprik and M.L. Klein, J. Chem. Phys. 89 ( 1988) 7556. [3] A. Wallqvist, Chem. Phys. 148 (1990) 439. [4] L.X. Dang, J. Chem. Phys. 96 ( 1992) 6970. [5] L.X. Dang, J. Chem. Phys. 97 ( 1992) 2659. [6] P. Barnes, J.L. Finney, J.D. Nicholas and J.E. Quinn, Nature 282 (1979) 459. [7] D. Van Belle, 1. Couplet, M. Prevost and S.J. Wodak. J. Mol. Biol. 198 (1987) 721. [8] M. Sprik. M.L. Klein and K. Watanabe, J. Phys. Chem. 94 (1990) 6483. [9] K. Krynicki, C.D. Gteen and D.W. Sawyer, Faraday Discussions Chem. Sot. 66 ( 1978) 199.
Physics 184 (1994) 163469
169
[lo] F.J. Vesely, J. Comput. Phys. 24 (1977) 361. [ 111 C.L. Brooks, B.M. Pettitt and M. Karplus, J. Chem. Phys. 83 (1985) 5897. [ 121 D.A. McQuarrie, Statistical mechanics (Harper and Row, New York, 1976). [ 131 W. Smith, CCP5 Info. Quart. 26 (1987) 43. [ 141 M. Ferrario and J.P. Ryckaert, Mol. Phys. 54 ( 1985) 587. [ 151 G. Ciccotti, M. Ferrario and J.P. Ryckaert. Mol. Phys. 47 (1982) 1253. [ 161 H.C. Andersen, J. Chem. Phys. 72 (1980) 2384. [17] S.Nod.J.Chem.Phys.81 (1984) 511. [ 181 M.P. Allen and D.J. Tildesley, Computer simulation of liquids (Clarendon Press, Oxford, 1987). [ 191 B.J. Palmer, J. Comput. Phys. 104 (1993) 470. [20] B.J. Palmer and B.C. Garrett, J. Chem. Phys. 98 (1993) 4047. [ 2 11 D. Eisenberg and W. Kauzmann, The structure and properties of water (Oxford Univ. Press, Oxford, 1969). [22] C.A. Coulson and D. Eisenberg, Proc. Roy. Sot. A 291 ( 1966) 445. [23] P.H. Betens, D.H.J. Mackay, G.M. White and K.R. Wilson, J. Chem. Phys. 79 (1983) 2375. [24] A.K. Soper and M.G. Phillips, Chem. Phys. 107 ( 1986) 47.