A Monte Carlo study of water clusters

A Monte Carlo study of water clusters

Volume 48, number 1 15 Bfay I977 CHEMICAL PHYSICS LETTERS A MONTE CARLO STUDY OF WATER CLUSTERS Michael R. MRUZIK Demment of Mmd& Science nnd Em?ne...

427KB Sizes 3 Downloads 201 Views

Volume 48, number 1

15 Bfay I977

CHEMICAL PHYSICS LETTERS

A MONTE CARLO STUDY OF WATER CLUSTERS Michael R. MRUZIK Demment of Mmd& Science nnd Em?neering. Massachusetts hstitute Cambridge. Massachusetts 02139, USA

of Technoology,

Received 24 May 1976 Revised manuscript received 13 April 1977

The free energy, internal energy and entropy for water clusters of 8 and 64 molecules have been calculated using the Monte Carlo method and both the Lie-Clementi (C-XII) and Stillinger (ST2) potential functions. Detailed structural data for a cluster of 64 water molecules was also obtained by averaging over a large number of confgurations. Average dipole orientations and radial density distributions indicated that individual water molecules were not particularly restricted to

futed orientations and that general ordering tendencies were small.

1. Introduction The absence of direct experimental information on the thermodynamics and structure of small uncharged liquid clusters containing only a few molecules is a major impediment to a full understanding of processes which occur in these systems. In addition, the standard theoretical techniques, such as liquid droplet theory [ 1,2], or the microcrystallite approach, employ somewhat arbitrary structural models which become increasingly unrealistic as cluster size decreases. Direct computer simulation using either the molecular dynamics [3] or Monte Carlo [4] methods is not subject to these restrictions and depends only on a function describing the potential energy between molecules. In this study, a Monte Carlo technique is employed to make the first computer simulation calculation of the free energy of formation of water clusters containing 8 and 64 molecules. Internal energies and entropies are also reported_ These calculations are performed using the empirically derived potential function of Stillinger [5] (ST2) as well as the LieClementi [6] (C-XII) potential function derived from ab initio Hartree-Fock studies. Structural data for a cluster of 64 water molecules is obtained from the average of a large number of configurations using the C-XII potential function.

Average dipole orientations and hydrogen and oxygen radial density distributions indicate that the individual water moIecuIes are free to assume a large number of nearly equal energy configurations and that tgeneraI ordering tendencies are slight. A discussion of the potential functions employed are given in section 2. The ComputationaI method is described in section 3 and the results are given in section 4.

2.

Potential functions

Since Monte Carlo results are a direct consequence of a given potential energy function, the nature and accuracy of these functions is a primary consideration. The Lie-Clementi [6] (Cc-Xii) potential function for the water dimer is given by +exn

=

Q2W,,

- 2e2(L/r18 + 02(e-bzrl + a3(e-b3’16

+ I/r,,

f I/r23 f I!Q4)

- 4Q”lr7*

+ L/r28 * L/r37 + I/‘q7) +aIe --W56 3 + e-bzrz3

+ e-b3r26

+. e-bz’w

+ e-b3rs

+ e-bz~aj

+. e-b3Qsb

(1) L7L

Volume 48, number 1

where the rii are the distances between point charges on different molecules. The model geometry and the parameters Q and all a& b’s, and c’s are given in ref.[6] . These parameters for the C-XII function were adjusted to agree with the ab initio Hartree-Fock potential energy surface. To correct for the neglect of electron correlation effects inherent in the Hartree-Fock treatment, an analytical fit was made [6] to the dispersion data of Professor W. Kofos obtained from perturbation theory and is included as the last three terms of eq. (1). Although a rigid geometry was assumed for the water molecule, this constraint has negligible consequences for the relatively weak binding of the water dimer. The only significant error is the complete neglect of multibody interactions which are attractive representing as much as 10% of the total [6]. The Stillinger [S] (SE) potential function for the water dimer has the form

where r. is the distance between oxygen atoms, and the subscripts a and b represent point charges on opposite molecules. The switching function (SW) and the constants (E, u, Qa and Qb) are exactly as given in ref. [S] The parameters [7,8]of the ST2 function have undergone extensive refinement by comparing the corresponding moiecular dynamics calculations with a wide range of experimental data for bulk water. As a result, the ST2 function for the water dimer has been “effectively” adjusted so as to include average multi-body contributions present in the bulk liquid. This imparts a tendency to overestimate (more negative) the enthalpies and free energies of smail clusters where the surface to volume ratio is large. Potential energies, therefL>re, will be as much as 10% too positive for the C-XII function and will be somewhat too negative when the ST2 function is employed *. Structurally, the ST2 function should exhibit a greater degree of tetrahedral coordination_

3. Method of computation The Monte Carlo algorithm of Metropolis [9] was employed to evaluate the canonical average potential energy (u) as defined by

172

15 May 1977

CHEMICAL PHYSICS LETTERS

(u)

EjUeduikTdsLIJk-“‘kT

d!Z ,

(3)

where d!Y2indicates that the integral is over all of configuration space. During the calculations, the molecules were restricted to remain within a spherical volume Vc abort the center of mass of the system. It has been previously shown [l&l 11 that a realistic constraining volume is about seven times that of the same number of molecules in the bulk liquid; and radial density distributions in this work indicate that, below about 350 K, the molecules naturally restrict themselves to configurations which are not hindered by the constraining boundary. The Helmholtz free energy change AF,___ associated with N water molecules coming together at constant temperature and volume to form a cluster in a center of mass coordinate system is

where Q,,_ is the canonical partition function for N water molecules and lacks three degrees of translation freedom for the cluster as a whole. To evaluate this expression, the potential energy of the cluster UN is rewritten as

u,=

cchuq, i>j

where uii is the pairwise potential energy between molecules i and i, and h is a dimensionless coupling parameter. The free energy change is then

=


dh

(6)

1

Vc’

since Qcm = Qcm.O = 1) and Q:$ I-C= = Qc.,,,.(x = 0). The rapid increase of KJNo) as X + 0 necessitated the use of the integral

aF,_

= 4

transformation

1

iu_j&t))?~~/~ d(h1/4)

(7)

vc-

* This agrees with computer claculations of internal energy for bulk water at 298 K which give; -8.5 kcd/moIe for the ST2 function [S J, -7.2 kcal/moIe for the C-XII function [6] (without long range corrections), and -8.12 kcal/mole experimentally.

CHEMICAL

Volume 48, number 1

PHYSICS

Typically, seven to nine discrete vaIues of W,(x)> approximately equally spaced for h1j4 on [O,l] were required to reduce the integral error of eq. (7) to beIow about 2% using a closed Newton-Cotes formula. The total free energy of a cluster Fc__ is then given by F cm. = @.c.ln. + F:$$

=‘V,

v,> ,

(8)

-where Fr$ @S is the free energy of an ideal gas with respect to a’center of mass constraining vohrme and has been previously derived [9 1. The average internal energy (EcJn.) and entropy GT,,) are obtained from Wc,?

= 3kT(N -

l/2) +

(9)

-

(10)

and (Sc,?

= ((&r,.)

- Fc,_)I~

4. Results and discussion Monte Carlo values of Helmholtz free energy, internal energy, and entropy are given in table 1. As expected from the discussion of potential functions, the free energy F,,_ obtained from the ST2 potential function is uniformly more negative (about 7%) than for the C-XII potential function_ Internal energy (EC_? exhibited similar behavior and, because of surface effects was consistently less negative than the bulk experimental vahre of -8.12 kcal/mole. Entropies for the ST2 function were slightly less positive (less

15 May 1977

LETTERS

freedom for the water molecules) than for the C-XII function. Entropies and internal energies became _ more positive at higher temperatures due to a greater tendency for disorder and less average cohesion_ For a ciuster of 64 water molecules, 300 000 to SO0 000 configurations were required to obtain average potential energies with standard deviations of around 2%. The same accuracy was obtained after 50 000 to 100 000 configurations for 8 water molecules. Ah free energies had standard deviations less than about 5%. Since obtaining accurate statistics on structural features often requires very long periods of observation [12,13], 3.0 million configurations of a 64 water molecuIes cluster at 263.2 K were generated using the C-XII potentiaI function_ After equilibrium was attained around 500 000 configurations, cahzulation of the average radial density distributions was begun and the resuhs are ihustrated in fig. 1. Radial density distributions near the center of mass (R 5 3.5 A) indicate that the oxygen atoms tend to accumulate at a fried distance with the hydrogen atoms pointing both inwards and outwards_ For larger radii less than about 7.5 A, there is a slight tendency for radial dipole ordering as evidenced by a shift in the average hydrogen density towards the cluster center with respect to the oxygen density. For large radii the hydrogen density is slightiy greater. To better elucidate the microscopic structure, statistics on the cosine of the angle 0 formed by the electric dipole of a water molecule and the radial direction pointing outwards were compiled for I.0 million con-

Table 1 Free energy, internal energy, and entropy with respect to a center of mass coordinate function, temperature, and chster size as determined from hfonte Carlo calculations Potential function

Temperature

ST2 ST2 C-XII

298 320 298

ST2 ST2 ST2 ST2 C-XII C-XII C-XII

255 263.2 273 298 263.2 273 298

C-XII

350

CK)

Number of molecides 8 8 8 64 64 64 64 64 64 64 64

Mc.m.JN (kcal/mole) -4.2 -

-4.26 -4.11 -3.94 -3.56 -3.26 -3.13 -2.81 -2.25

system are given as a function of potential

i/N ~ZZ&ole)

&_m,lN (en)

-9.93 -

4.25 -3.81 -3.85

i9.L -

-11.27 -11.40 -11.57 -12.04 -10.55 -10.75 -11.28 -1253

-7.13 -6.95 -6.73 -6.21 -5.35 -5.l? 4.54 -3.86

16.2 16.9 17-7 19.6 19.8 20.6 22.6 24.8

Fc.m.IN (kcal/mole)

fj

,

i, 20

R*.

I

15 May 1977

CHEMICAL PHYSICS LETTERS

Volume 48, number 1

-

O*yqen

---

Hydrogen

L,

40

60

RADIAL

80

100

DISTANCE

,, 120

140

SQUARED

160

for radii less than about 10 A, the average value oscillates about cos 8 = 0. For Iarger radii both quantities are positive. Although these quantities indicate only a rather sligkt preference for inward dipole ordering (radius less than 10 A) with possibly skewed distributions of vaIues for cos 0, both quantities are insensitive to small amounts of preferential ordering and would miss such features as bimodal distributions about cos 0 = 0. As a result the total distribution of values for cos t? is given for seven different radii in figs. 3 and 4. Overall, the values of cos 0 are fairly evenly distri-

160

ti2,

Fig. 1. Average densities of oxygen and hydrogen atoms are given for a cluster of64 water molecules at 263.2 K using the Lie-Ctementi (C-XIlJ potential function [6]. LO-

0806-

04-

02-

8 B

AVERAGE

VALUE

0 bl443b

-02

to 626

i

‘ogt_l_lliLI -10

MOST PROBABLE 0

10 20 3040

-08

-06

-04

-02

0

02

Q4

06

06

IO

VALUE

50 6070 80 90IDDII0J201300140150160 R’.RADIlJS SOLJARED Iii21

Fig. 2. Average and most probable vahes for cos .9 are given as

a function of radius squared. cl6

26dlo?6?%

figurations. Both the average value and the most probable value? of the cosine are given as a function of radius in fig. 2. While the most probable value is negative ?’ For each confiiration the peak (most probable) value of cos e was determined. The most frequently cited value from all confiirations is referred to here simply as the “most probable” value.

174

01 -10

I -08

I -06

I -04

I -02

I

1

0 cos

02

I

I

04

06

I 08

, IO

(81

Fig. 3. Relative densities of values of cos 0 are given for different distances from the cIuster center as indicated.

Volume 48, number 1

CHEMICAL PHYSICS LETTERS

15 May 1977

viously reported [ 121, convergence of structural features (espec~aliy for small radii) is very slow. Q&ratelike structures were not observed in any detailed arrangement of water molecules studied.

0’ -10 M

-08

1 -06

I -a4

I -02

I 0 cos(6~

t

z g

AcknowIedgement I

I 02

I 04

f 06

I 00

J 10

20

!A w ‘0 2_

ii0

I

I

d

-10

-08

-06

I -04

I -02

t 02

I

I 04

I 06

I CX3

J IO

04

06

00

I IO

08

to

CihH

=20-

cl9908

to 1084

B

The author is thankful to IBM Research Laboratory, San Jose, for providing the necessary facilities to include time on an IBM 360/195 computer and to the Women7s Auxiliary of the A.I.M.E. for financial support. The advice and assistance of Dr. Farid F. Abraham of IBM San Jose and Professor G.M. Pound of Stanford University is gratefully acknowledged.

IO0 -IO

-08

-06

-04

5 ?I 0

-02

0 caste)

dl1084& -08

-Q6

-04 t

-02

-

02

References

to II 71 li 0

a2 $7

CM 06;-,

Fig. 4. Relative densities of values of cos 0 are given for different distances from the cluster center as indicated.

buted, being onIy siightiy skewed, and peak vafues do not correspond to perfect radial dipole ordering. In summary, the results are: (1) CaIcuIations which give the upper and lower limits of free energies of formation of water molecule clusters have been accomplished. Differences are primarily caused by the neglect (C-XII function) or overestimation (ST2 fonction~ of three body interactions. (2) As demonstrated by the fairly flat distribution of cos 8 values and average values of cos 0 near zero, general ordering tendencies are slight and water molecufes are able to assume a very large number of different but nearly equal energy configurations. As pre-

fl] R. Becker and W. Doring, Ann. Phys. 24 (1935) 719. [2] J. Lothe and GM. Pound, J. Chem. Phys. 36 (1962) 2080. 131 A. Rahman, Phys. Rev. 136 (1964) A40-5; L. Verlet, Phys. Rev. 159 (1967) 98. [4] W.W. Wood, in: Physics of simple liquids, eds. M-N-V. Temperley, J-S. Rowhnson and G.S. Rushbrooke (NorthHolland, Amsterdam, 1968). [ 5 ] FX. StiBiier and A. Rahman, J. Chem. Phys. 60 (19741 1545. 161 G.C. Lie and E. Clementi, 3. Chem. Phys 62 (1975) 2f95. [7] J.D. Bernal and RX. Fowler, J. Chem. Phys. l(L933) 515. IS] N. Bjerrum, Kgl. Danskk Videnskab. Selskab Math. Fys. Ilfedd. 27 (19.51) 1. [9] N. Metropolis, A.W. Rosenbtuth, M_N_Rosenblnth, A-H. Teller and E. Teller, J. Chem. Phys- 21 (1953) LO87[lo] J.K. Lee, J.A. Barker and F-F. Abraham, J. Chem. Phys 48 (1973) 3166. [ 111 M.R. Mruzik, F.F. Abraham, D-E. Schreiber and G.M. Pound, J. Chem. Phys. 64 (1976) 481. [ 121 F.F. Abraham, D.E. Schreiber and J-A. Barker, J. Chem. Phys. 62 (1975) 1958. [13] F-F. Abraham, J. Chem. Phys. 61(1974) 1221.