Monte Carlo liquid water simulation with four-body interactions included

Monte Carlo liquid water simulation with four-body interactions included

21 Dcccmber 1984 CHEMICAL PHYSICS LETTERS Volume 112, number 5 MONTE CARLO LIQUID WATER SIMULATION WITH FOUR-BODY INTERACTXONS INCLUDED J. DETRICH,...

388KB Sizes 12 Downloads 80 Views

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.