21 Dcccmber 1984
CHEMICAL PHYSICS LETTERS
Volume 112, number 5
MONTE CARLO LIQUID WATER SIMULATION WITH FOUR-BODY INTERACTXONS INCLUDED J. DETRICH, G. CORONGIU and E. CLEMENT1 IB;II Cixporatiorz, DSD. Department 48B/Mail Statiorl 428. Neighborhood Received
21 August
1984:
in timal form
13 October
Road.
Kingston,
New
Sork 12401,
US.4
1984
A four-body interaction potential for water molecules is derived. The new terms are combined with previously dcvclopcd two- and three-body POtCntht tcrrnsond upplied in n hlctropolis-Monte Carlo simulation at 298 K for nn (N, I’,?“) ensemble with 513 water molecules. Improvements, relative to the results using only two-body, or two- and three-body interxtions, are rcportcd for the correhttion functiong(O-O), for the S-ray and neutron-beam suttering intensities, and for the entllalpy.
1. Introduction For sonic time now we have been deriving interpotentials from ab initio computations and
moleculx
usiug tlicse potentials in statistical mechanical computer espcrimcnts on pure solvents or solutions [l], with the goal of deriving an ordered set of approsimations hding to an increasingly realistic description of the liquid. An hportant step along this road is the recent wurk of Clementi and Corongiu [?I. which obtains a h-cc-body interaction potential for water from ab initio calculationsand applies thisin Metropolis-Monte Carlo simulations of liquid water. The work reported here is a direct estcasion of their work to yield the ucst devclop~~~ent in this sequence, namely the derivation aud application of an interaction potential for wstcr as two-
nwlecuIes and
that
includesfour-bod_v
three-body
tcrnts
as well
tcrnrs.
2. Potential with four-body interaction temls For an N-body systtnl intermolecular interaction scrics of the forin
we can represent the total enel-gy, U. in terms of a
(1) 426
where the summations range from 1 to N, and the various Uck)arc defined so that they are necessarily zero in caseN
(3 here we use 6 = (Y,, - (Ye, where (Ye and (Y,, denote. respectively, the polarizabihty constants perpendicular and parallel to the water OH bonds, i is the index denoting the ith water nlolccule, b is the index used to distinguish the two bonds in a given molecule, ebi’is the unit vector in the direction of the bond bi, and El,i is the electric field at the midpoint of the bond bi. Note that hyperpolarizability terms [3] of quadratic aud higher order in the electric field strength are not considered in this work. The electric field Ebi originates front the assumed permanent electrostatic moments of the surrounding water n~oleculcs; it can be expressed as a sum Ebi = zEbi j+i
j .
’
where E&i is the electric field component at the bond midpoint bi originating from the permanent electrostatic monients of thejth water molecuJe_ As in earlier work [z], we account for thesi: moments in ierrns of charges situated on the 0 and H nuclei, so we have (4)
where Rbi . Ai is the distance vector pointing from the nucleus X (0 or H) in molecule j to the bond midpoint bi, and qA denotes the charge situated on nucleus h. To first-order in the polarizability, this model yields the energy given by
If we substitute eqs. (2) and (3) in this espression, we obtain a triple sum which shows the three-body character of our interaction energy_ Continuing now to second order in the polarizability, we find m inter8ction energy which is conveniently expressed in terms of dipole-dipole interactions. namely
-
21 December 1984
CHEWCAL PHYSICS LEXTERS
Volume 112, number 5
3(rbi,b5.’pbi)(‘bi_b’j’pb'j)/'&'jl
(6)
*
where ‘bi,b5’ is the distance vector pointing from bond midpoint bi to bond midpoint by_ We cm again substitute eqs. (2) and (3) here to obtain an expression with a four-fold summation which clearly shows the four-body character of the interaction energy in eq.
Analogous calculations to support our expression for the four-body interadtion would be an ovenvhelming task, since the additional water’molecule adds more degrees of freedom to be considered in the geometry scan, and each calculation would involve a basis set augmented to accommodate the additional water molecule_ On the other hand, we have seen that the fourbody term we propose is already implicit in the model adopted by Clementi and Corongiu, so their fitted parameters must also be carried over to the new inreraction term. The only question concerns their shortrange correction, since our argument is based on classical electrostatics, and short-range effects most likely are due to non-classical effects such as exchange interaction_ To settle this question, and also to verify that our model does in fact account for the four-body interaction, we have carried out several ab initio SCF calculations for-water tetramers, using a large basis set and correcting for superposition error by the counterpoise technique_ Some of the results are given in table 1. it is seen that our model (without the short-range correction) does indeed provide a good description of the four-body interaction, and also helps somewhat in improving the description of the three-body interaction. For the two-body interaction, we adopt the Matsuoka-Clementi-Yoshinine (IMCY) ab initio potential [5]. This of course facilitates comparison with previous results. In addition, the considerations that motivated its adoption by Clementi and Corongiu are still valid, in our opinion_ It should be noted that
Table 1
(6). The actual parameters Corongiu
in realizing
used by Clementi
their expression
and
for the tliree-
body interaction term arc derived from a fit to ab initio SCF calculations on the water trimer at 173 different g&ometries; to ensure accuracy, each of these calculations uses a large basis set and corrects for superposition error by means of the counterpoise technique [4] _Their fit yields the parameters ai = 3.9 1 au, 6 = 1.42 au, q. = -0.6436 au, and qH = 0.32 18 au. In order to account for short-range effects, they modify eq, (4) by replacing the denominator Rzihj by R& Ai - C,, and-use this modification in eqS: (2) and (5) as Wek they adopt-the Values CH = 0.0 and co = 3.038 A_
Comprison of tbrce- and four-body inreraction energies estracred from model interaction L’p and Up:! with results from ab initio SCF calculations for selected warcr tetramer geometries, in kcd/mol %
qN -!- $2
SCF
three-body four-body
-1.26s -
-1.665 -0.552
-1.713 -OS00
three-body four-body
-0.925 -
-1.113 -0.383
-1297 -0.361
three-body iour-body
-0.690 -
-0.892 -0.270
-0.966 -0.260
three-body four-body
-0.524 -
-0.668 -0.192
-0.762 -0.186 427
Volume 112, number 5
21 December 1984
CHEhIlCAL PHYSICS LEI-lXRS
since we do not expect two-body interaction terms occurring in UP1 and UP2 to improve the description of the two-body interaction provided by the MCY potential, these terms are omitted from both UpI and UP1 in our calculations.
L exp _______ MCY
26 -t 21
t
Ii
3. Monte Carlo simulation
The Metropolis-Monte
Carlo simulation [6] was
carried out with an (N. V.T) ensemble for T= 298 K and IV = 5 11 water ~~~oleculcsat an experimental dcnsity of 0.998 g/cm3. Cubic boundary conditions with a minimum image cutoff have been used [7]. To hasten equilibration, we started the Monte Carlo simulation from a configuration obtained from an earlier Monte Carlo simulation which was equilibrated without the new interaction terms discussed above. We
2.8 2.1 0
6
0.7
sys1c111 1IO 1. WC prcscnt our fig. I md tltblc 2. functiong(O--0j (designated hlCY
results in figs. I-3 and table 2. In WC compare tlx pair correlation ohtaincd using our interactions + CC + DC) with earlier results [2]
obrainod wirh the two-body hlCY interactions, and with the two- and threebody interactions (designated hlCY and hlCY + CC, respectively). These results are ~llw compared with those obtained by Nnrtcn [ 1 I] mudcling S-ray and neutron-beam diffraction data. 111
a; >
t
1, i
1
i r
I _.. y-~=L__
c----
-
wh
0.oL-J
retained the optimized step size of 0.30 a. and the optimized angle of 20’. In all. lo6 steps wcrc carried out: rhc first half of the run was attributed to equilihtxion, and our statistics are based on tlic last half of the run. For a detailed description of the techniques used in our Monte Carlo programs we refer the reader clscwlrerc [S] _Tile c&xAations reported hcrc were run o11m cspcrimental system consisting of a11IBM 133 i. IWO IBM 4341s and ten FI’S-164 attached processors. in a configuration so that parallel processing cm bc performed using as many as ten attached processors performing simultaneously on one run [9]_ Tlik system gives sul~crcornputer-tyI,e performance WI large runs, which is important, since four-body crtlcuiations of ~llc type described iicrc impose very Itcavy demwds on computational rcsourccs, and would bc virtually impossible without 311extremely powcful computer. WC have already reported the migration of our hlon tc Carlo codes to this parallel
1.4
_ exp ________ MCY +CC
t r
1
1
0
2.8
’
’
1
1
-
exp ___--- MCY+CC+DC
2.1
0.0
I
t .I
9
I
3.0
I
6.0
0
.
.
/
9.0
12.0
Fis. 1. Comparison of pair corrclatim functionsg(O-0) obtained by modeling X-ray diffraction data (Narten) with simulations using the hlCY potential (top), the (hlCY + CC) potential (middle), and the (h1CY + CC+ DC) potential (bottom).
figs. 7 and 3 we present X-ray and neutron-beam the corresponding defined by H(s) = [I(s) - @)I whcrc
more direct diffraction
structure
functions
comparison with data in terms of H(s).
/G))
s is given in terms
Thcsc are (7)
of the scattering
angle 29 and
the wavclcngth h by s = 4~ sin 6/X, I(s) is the coherent scattering intensity, and @> is the scattering intensity for one independent molecule averaged over all orientations_ Details of these computations can be found in previous work [ 12]_
Volume 112, number 5
21 December
CHEMICAL PHYSJCS LE’ITERS
1984
--___ ---, ___ __ __, exp _
MCY
!._!,
I
‘.i
1:
:
.p.
,,.--.
_....._____
~_$=------
;
B s r
0.6 -
2 f =, 02 s5 = -4x2 3 2
;
exp
MCY
_--
n
4
cc
B ._
-0.6 t -1.0 0.6
:_ :! :,
! -r’
S (A-‘)
t
f
I
.
I
.
I
, I
?:z-:I: exp MCY-CC
t
-0.6
_I
Z exp _____ MCW
-0.6
_
__--.---- Iexp
CC+DC
-
1
I
4.0
I
I
8.0
. S(A-‘) . . 12.0
16.0
Fis. 2. Compxison of X-=3x structure functions H(s) from cspcrimcnts (Nartcn) with simuhtions using tile MCY POtential (top), the (MCY + CC) potential (middle), and the (WX + CC + DC) potentiril (bottom).
2.0
N&en 111 j MCY 1131 @ICY f CC) 121 @KY + CC + DC) a)
I'irSl
First
maximum
r~i~irnuI1~
4.0
6.0
8.0
10.0
Fig. 3. Comparison of neutron structure functions H(s) from eqxriments (Nsrten) with simulations winy the %CY porential (top), the (&ICY + CC) potenrial (middle), nnd the (&KY + CC + DC) pOtC,,tiai (bottom)_
Table 2 Intensity I at the radial distance R (ii A) of the mlisima -ad minima of the pair Source
MCY+CC+DC
corrdation
function g(O-0)
Second msximum
Second I~~iIn~~u~l
R
I
R
I
R
I
R
I
2.83 2.83
2.31 2.46 2.66 3.50
3.45 3.53 3.37 3.34
0.85 0.94 0.75 0.58
4.53 4.25 4.50 4.32
l-12 1.08 1.17 126
5.6 5.6 5.65 5.61
0.86 0x9 0.82 0.81
2.Sl 2.72
a) This work. 429
Volume
112. number
CHEMICAL
5
PHYSICS
4. Discussion
i~ilcrvals bctwcc~~peaks in I-l(s) for X-ray scattering to Ic~~gthen. ;1 trend
that does not enbancc agreement wiib 111ecxpcrimental results for large s. For snxdls, on !iic orhcr hand, tbc increased structure leads to Inorc rcalisrn for tlic niodcl; this is Inost clearly seen by noting the initial slope and first split peak in H(s) for X-lay scattering. While hlCY + CC gives a definite improvcnlcnt over MCY, MCY + CC + DC is bcttcr
WC obtain 311avcragcintewal
Ilie amplitude energy
of -10.73
of f
I kcal/niol,
yielding 311cnlhalpy of -S.95 2 0.01 kc;~l/~~wlc. This shoulii be coniparcd with the value of -4.6 2 0.2 kcal/~:~:~l for MCY [ 121 and -7.70 i 0.05 kcal/~~l for MCY + CC [II]. TIK cspcriulct~titl 0.0
21 Dcccmbcr
1984
References
Frond fig. 1, it CJII be seen that progression from SKY. ~hrougl~ hlCY + CC, to MCY + CC + DC corresponds to a trend toward more structure ing(O-0). Tbc amplitude of the lirst peak also increases in this progression (as can bc seen from table 2 as well), and this is rcllcctcd in a corresponding tendency for the
still: ir is tlie only result niatcliiug boIlI I~alvcs01‘ this s@it peak.
LETTERS
v3luc is 4. I kcal/niol, so we can cstimatc that elTects wglcctcd iu tllis work give ;I contribution of about 0.8 kc;d/nlol. One such rffccl ori@natcs froni ~nolecular vibrations; this should also help by broadening (and lowering) tlw peaks prcscntcd in fig. 1 and table 2.
[ 1] E. Clcmenti. Computational
nspccts for lar_ee chemical systems, Lecture Notes in Chemistry. Vol. 19 (Springer, Berlin. 1980). [2] E. Clemcnti and G. Corongiu. Inrcrn. J. Quantum Cbem. 10 (1983) 31. [3] A.D. Buckingham, Advan. Chcm. Phys. 12 (1967) 107. [4] S.F. Boys nnd F. Bernardi, hlol. Phys. 19 (1970) 553. [S] 0. hlatsuoka, E. Clcmcnti and hf. Yosbiminc, J. Cbem. Phys. 64 (1976) 1351. [6] N. hletropolis, A.W. Rosenblutb. A.H. Teller and E. Teller, J. Cbem. Pbys. 21 (1953) 1078. [7 ] J.P. Vnllenu and S.G. Wbirtington, in: Statisti= mccbnnits, Part A. Equilibrium tecbniqucs, cd. B.J. Berne (Plenum Press, New York, 1977) p. 150. [8] J. Fromm, E. Clemcnti and R. Watts, IBM Research Report. San Jose, California (1973); R. Barsotti and E. Clementi, hlontcdison Technical Report. lstituto G. Donugani, Novara. Italy (1976); S. Roman0 and E. Clementi. Gazz. Cbim. Ital. 108 (1978) 319.
191 G. Corongiu and J.H. Detrich, to be published. [IO]
J. Derrich, G. Corongiu nnd E. Clementi, Intern. J. Quantum Cilcm S, to be publisbcd . [I I] A.II. Nartcn and lI.A. Levy. J. Clrem. Pbys. 55 (1971) 2263; A.H. Narten, J. Chcrn. Phys. 56 (1972) 5661; A.H. Narten, 0RNL-4578 Report, Oak Ridge Natl. Lab., Cbcmistry Div. Oak Ridge. Tennessee (1970). 1121 G.C. Lit. E. Clemcnli and hl. Yosbiminc, J. Chem. Pbys. 64 (1976) 2314.