New ab initio VB interaction potential for molecular dynamics simulation of liquid water

New ab initio VB interaction potential for molecular dynamics simulation of liquid water.

M. RAIMONDI, A. FAMULARI, E. GIANINETTI, M. SIRONI, R SPECCHIO, AND I. VANDONI Dipartimento di Chiniica Fisica ed Elettrochiinica and Centro CNR CSRSRC, Universitd di Milano. via Golg 19, 20133 Milano, Italy.

Abstract

SCF-MI (Self Consistent Field for Molecular Interactions) and non orthogonal CI were used to determine a water-water interaction potential, from which BSSE is excluded in an a priori fashion. The new potential has been employed in molecular dynamics simulation of liquid water at 25°C. The simulations were performed using MOTECC suite of programs. The results were compared with experimental data for water in the liquid phase, and good accordance was found, both in radial distribution hnctions and thermodynamic properties, as well as in geometric parameters.

I. Introduction II. The water-water interaction potential 111. BSSE free ab initio intermolecular interaction potentials IV. Results and discussion V. Ab initio calculations VI.Molecular dynamics simulation VII. Simulation results VIII. IR and Raman spectra IX. Conclusions References

ADVANCES IN QUANrVM CHEMISTRY, VOLUME 32 Copynght b 1999 by Academic Press. All rights of reproduction in any f o n resewed

0065-3276/99530.W

2 3 5

9

11 12 17 20

263

264

M. Raimondi eta/.

I. Introduction

Statistical classical mechanics states that, in principle, it is possible to deduce macroscopic system properties from its time evolution in the phase space. This requires nuclear potential energy as a function of the atomic coordinates of the system to be known. To determine the shape of the potential energy, models are considered which acceptably reproduce the greatest number of experimental physical properties of the system. The parameters included in the analytical expression of the potential are chosen so as to be consistent with experimental data (semiempirical potentials) or with the energy value calculated by means of ub initio methods. Models of this kind are employed in molecular dynamics simulations, to determine static and dynamic properties of a given system; this permits to veri@ the consistency and reliability of the model, to explain the molecular origin of physical characteristics of a system or to calculate properties for which no experimental data are available. Great interest has been devoted, for example, to the study of the structure and dynamics of liquid water, ice and solutions. Liquid water simulations are used to represent the hydrogen bond network in the liquid phase and to explain some unusual and not yet understood properties of water, such as its high thermal capacity, in terms of its structure. Still more important is the elucidation of the molecular origin of dynamic phenomenons such as the solvations of ions in water and the interaction of macromolecules with the solvent. Rahman and Stillinger (1971) have applied to the study of liquid water a model potential, named ST2, which could anticipate the existence of the so called “high frequency sound” years before this collective excitation of water molecules was observed, by means of inelastic neutron scattering experiments (Texeira, Bellisent-Funel, Chen and Darner,( 1985). Lie and Clementi (1986) obtained a notable model potential for water from extensive ub initio calculations. This statistical classical mechanics approach gives rise to a series of problems. It is questionable whether a computational method which assumes the validity of classical mechanics can reproduce the physical properties of a system. It is well known, in fact, that the study of intramolecular vibrations requires a quantum mechanical approach. Two routes are available to obtain such quantum mechanical description: a simulation can be realised with the path integral molecular @numics method (PDdD), recently applied to liquid water (Billeter, King and van Gunsteren, 1994), or a classical dynamics simulation can be realised, with an a posteriori corrections of the results (Berens, Mackay, White and Wilson, 1983; Berens and Wilson, 1981). The latter choice has been made for the calculation of gaseous CO spectra in liquid argon (Berens and Wilson, 1981). Moreover, the phase space volume considered in the course of a simulation is truncated and some properties cannot be obtained with the desired accuracy. The

Ab lnitio VB Interaction Potential of Liquid Water

265

molecular dynamics simulation evaluates the properties that are a consequence of the motion of an isolated particle more accurately than those which originate from collective phenomenons. The former, in fact, can be averaged among all the system molecules as well as over different time intervals. Thermodynamical properties are thus determined less accurately than the internal molecular geometry or the molecular dipole moment. In the present work a model potential for water will be described which, once inserted in a molecular dynamics simulation, excluding the value of the internal pressure, satisfjringly reproduces the physical properties of the system at ambient temperature and density. A general ub initio VB approach, based on the SCF-MI (Self Consistent Field for Molecular Interactions) procedure (Gianinetti et al., 1996), was used to determine the water-water pair interaction potential, from which BSSE is excluded in an a priori fashion. The recent generalisation of the SCF-MI method (Gianinetti et ul., 1997) allowed to compute without BSSE also non-additive many body contributions to the total polarisation interaction energy. The simulations were performed using MOTECC suite of programs (Sciortino, Corongiu and Clementi. 1994). The results were compared with previous simulations and with experimental data for water in the liquid phase. We underline that this scheme allows a complete elimination of BSSE in an apriori fashion, a fact which seems of vital importance to generate a satisfyingly reliable a b initio potential.

II. The water-water interaction potential In this paragraph a short summary of NCC (Niesar, Corongiu, Clementi, Kneller and Bhattacharya, 1 990), a polarizable and flexible water-water interaction model first proposed by Niesar, Corongiu, Clementi, Kneller and Bhattacharya, (1 990), will be presented. According to this approach, the interaction energy among water molecules in the liquid phase can be expressed in terms of a many-body potential (Niesar, Corongiu, Clementi, Kneller and Bhattacharya, (1990):

u = v + v,,, (1).

The term V contains two contributions arising from intermolecular and intramolecular motions and is defined as:

M. Rairnondi eta/.

266

where i and j define all the molecules of the system, rap are the intermolecular distances and qia are the internal coordinates of the molecule. Coupling between intermolecular and intramolecular motions are considered negligible, so that the parameters ak can be regarded as constants. In the following discussion an analytical form for the additive portion of the potential will be considered, which derives from the point charge model originally proposed by Bernal and Fowler (1933). According to this model, positive charges q are located at the position of the two hydrogen atoms of each water molecule (1, 2, 3 and 4), and a negative charge -2q is positioned at a point P along the CX symmetry axis of the molecule (7, 8) at a fixed distance from the oxygen atom ( 5 , 6) (see Fig. 1). The analytical expression of the pairwise additive part of the potential is therefore (Niesar, Corongiu, Clementi, Kneller and Bhattacharya, (1 990):

I

Fig. 1. Water dimer atom labels.

The geometry of the rigid water molecule corresponds to that of an isolated water molecule in the gas phase (OH bond length = 0.9572 A; HOH angle = 104.52') (Benedict, Gailar and Plyler, 1956). The terms ApH and APO in Eq. (3) are needed to correct the model potential behaviour at short distances. As regards intramolecular motions, the ub initio D-MBPT(oo) potential is employed (Bartlett, Shavitt and Purvis, 1979), which is expressed, up to the quartic terms, as a function of the three internal coordinates of water, and is obtained from many body perturbation theory.

Ab lnitioVB Interaction Potential of Liquid Water

267

The term V,, in Eq. (1) is a polarisation contribution which includes the many-body non-additive effects. It has been demonstrated (Clementi and Corongiu, 1983) that the many-body corrections are necessary for an accurate quantitative prediction of the system properties. These effects are included by ascribing to each molecule a dipole moment and a bond polarizability, according to the classical model (Bottcher, 1973). Within this scheme the polarisation energy can be expressed as:

where: Pi,

= aii

eiA

and: N, is the total number of system molecules; N, is the number of induced dipole moments per molecule; pi, is the induced dipole moment on the h-th bond of the i-th molecule; tZiAis the polarizability of the h-th bond of the i-th molecule;

Ei, is the total electric field acting on the h-th bond of the i-th molecule; Aiming to obtain a non-empirical potential to be used for liquid water simulations, the parameters of the NCC model described above are determined by fitting to extensive ah initio calculations, which will be presented in a following section. III. BSSE free ab initio intermolecular interaction potential A general ah initio VB wavehnction consisting of single and double excitations out of the SCF-MI determinant, are used to determine an ab initio interaction potential apt to the dynamics simulation of liquid water. A short summary of the theory involved will be presented. For a more detailed account see (Gianinetti et al., 1996; Famulari, 1997; Specchio, 1997). The chosen scheme permits to eliminate BSSE in an a priori fashion, both at the SCF and at the correlated level. The BSSE is a well known inconvenience which strongly s e c t s the intermolecular potential investigations of weakly interacting systems, such as van der Waals or hydrogen bonded systems. In such systems, it is common that the BSSE is of the same order of magnitude of the interaction energy involved. The

M. Raimondi eta/.

268

most frequently employed a posteriori method to correct this error is the counterpoise technique (Boys and Bernardi, 1970). Nevertheless, many authors have emphasised that the method introduced by Boys and Bernardi does not realise a clear and precise determination of the BSSE magnitude (van Duijneveldt, van Duijneveldt-van de Rijdt and van Lenthe, 1994). In addition, it is worth emphasising that the SCF-MI “a priori” strategy to avoid BSSE is always correct and not only for full CI wavefunctions as in the case of the CP a posteriori procedure (van Duijneveldt, van Duijneveldt-van de Rijdt and van Lenthe, 1994). Furthermore, the addition of the partner’s functions introduces other spurious contributions - the “secondary superposition error” - owing to the improper modification of multipole moments and polarizabilities of the monomers. This becomes particularly important in the case of anisotropic potentials where these new errors in the wavefunction introduced by the CP procedure can contribute to alter the shape of the PES and the resulting physical picture. Recently, the SCF-MI algorithm (Gianinetti et al., 1996) has been implemented (Famulari et al., 1997a) into the GAh4ESS suite of programs (Schmidt et al., 1993). Within this scheme the one determinant wavefunction of the supersystem AB is written as

..5tA (2NA)@:(2NA + 1)3:(2NA + 2)... ...mi,(2NA +2N,)]

ykB= (N!)-S 1Pt(l)G?(2).

( 6) where the molecular orbitals of each fragment are free to overlap with each other. The orbitals of fragment A and the orbitals of fragment B (@:..@:,)

are expanded in two different subsets of the total basis set

x = (xr)Ll, xA =

I$),=,

M A

M.

centred on A and

xB = (x:)~’centred on B: q=1

i = l...N,

p=l

j = l...N,

269

Ab lnifioVB Interaction Potential of Liquid Water

These restrictions obviously permit the complete elimination of BSSE in an a priori fashion, by assuming and maintaining the orbital coefficient variation matrix in block diagonal form:

(8).

From these premises it follows that the general stationary condition is equivalent to the coupled secular problems: FLT, = SaT,L, TiSLT, = I,

{

FLT, = S;T,L, TlSLT, = I,

(9) where the effective Fock and overlap matrixes F' and S' possess the correct asymptotic behaviour. By means of a general variational approach, the interaction energy is expressed as:

To take the correlation contribution to water dimer interaction energy into account, a very compact multi-structure VB - non-orthogonal CI calculation is carried out. With this aim, a method to obtain virtual orbitals corresponding to the SCF-MI wavefunction which furnishes an energetic contribution of the same accuracy as very extended CI calculations, was developed. In order to avoid that BSSE switches on, both the occupied and virtual SCF-MI orbitals are constrained to the restrictions

-

M A

= Cx;T;, p=l

Ma

= Cx;TZ, q=1

M A

0;= Cx;TPt. ,

MB

= Cx:Tq;. q=l

p=1

(10).

The above constraints imply the non-orthogonality of the orbitals. The resulting wavefunction is assumed to have a general Valence Bond form: a=1 b=l

a=l

b=l

It represents the configuration interaction between the SCF-MI wavefknction

(1 1).

270

the singly excited localised configuration state functions

M. Raimondi eta/.

(121,

and the doubly excited localised configuration state functions

obtained by simultaneous single excitation localised on A and B, where the singlet spin function for the two or four electrons involved in the single or double excitation are: The virtual orbitals @ ,: and a :, are determined according to the approximation that they maximise separately the dispersion contribution of each of the two configuration wavefunctions

where “a” labels one specific occupied MO of fragment A and ‘8” one specific occupied MO of fragment B. The dimension of the virtual space so generated is given by the product of the active occupied MOs of fragments A and B, respectively. The associated virtual orbitals are determined at the variation-perturbation level of theory by minimisation of the second order expression of the energy:

by satisfylngthe conditions

Ab lnitioVB Interaction Potential of Liquid Water

271

The minimisation strategy is the specialisation to the case of a SCF-MI zero order wavefbnction of the algorithm previously described by Raimondi et al. (1 996) in the context of the Spin Coupled theory. Acceleration of the minimisation procedure now includes the approximate evaluation of the Hessian as described by Clarke et al. (1997). It should be underlined that in this scheme the a priori exclusion of BSSE is still operative, as the virtual orbitals are expanded on the two distinct subsets centred on A and B, see Eq. (10). The virtual orbitals obtained by means of the perturbation theory approximation are then employed to construct a fbll set of all the singly- and doubly-excited configurations which provide the final VB-like wavefbnction Eq. (1 1). Finally, the multi structure VB (non-orthogonal CI) problem is set up and solved variationally according to standard VB techniques, see Raimondi et al. (1977) and Cooper et al., (1987). In the present case, for each water molecule we considered an active space of four MOs with one MO (oxygen ls2 electrons) keptpozen. The resulting dimension of the virtual space of each fragment is sixteen, corresponding to four virtual orbitals for each occupied MO. This implies a set of 32 vertical singly excited configurations. By adding all the corresponding vertical doubly excited spatial configurations the final VB-like wavefbnction was generated. By taking the dimension of the spin space into account, the size of the resulting VB matrix is 545.

Iv. Results and discussion V. Ah initio calculations

The water dimer interaction energy was investigated at the Hartree-Fock and correlated level using SCF-MI and the VB-like wavehntion described above. The 7s4p2d4s2p basis set proposed by Millot and Stone (1992) was employed in both cases. This basis set reproduces well the Hartree Fock limit results obtained by Famulari et al. (1997b). Optimal geometry, calculated employing the VB-like wavehnction Eq. (1 l), is in good accordance with the experimental data (see Table I and Fig. 2). The binding energy value of -4.67 KcaVmol falls within the experimental range (see Table 11).

M. Raimondi eta/.

272

VB

R o-o(A)

Oxygen (a) 7s4p2d

Present work 3.00

Po

134.5


2.5

Hydroge n (a) 4s2p

Total 102

Exp. (a) 2.98M.03 122.0*10 0.0f10

SAPTMBPT4

mw

mE,

Present work KcaYmol -4.67

2.953

MP4 (4 2.949

124

124.8

6.8

5.35

(b)

(b) KcaVmol -5.2rL.O.7

A h (4 Kcdmol -5.44kO.7

All the parameters which appear in the NCC potential expression are determined by the present VB ub initio calculations. In this way a new non empirical potential used in the simulation of liquid water was obtained. The calculations include the determination of the interaction energy relative to 225

Ab lnitioVB interaction Potential of Liquid Water

273

configurations of the water dimer and of the non additive three-body contributions relative to 28 trimer configurations. As regard the calculations on the trimer, the quantity in question is expressed as follows: A Enm-,dditiva (1,2,3) = A E (1,2,3) - A E( 1,2) - A E (1,3) - A E(2,3) =

E(1,2,3)-E(1,2)-E(l,3)-E(2,3)+E(l)+E(2)+E(3)

(18) where E(1,2,3) is the Hartree-Fock energy of the (1,2,3) cluster. Also in this case, the interaction energy is BSSE free, as the calculations are performed by employing the n-body extension (Gianinetti, Vandoni, Famulari and Raimondi, 1997; Famulari, 1997) of the SCF-MI method. As it has been shown elsewhere (Habitz, Bagus, Siegbahn and Clementi, 1983), the main contribution to the threebody non-additivity derives from induction effects, while corrections due to electronic correlation are mostly additive. For this reason, such contributions are determined at the Hartree-Fock level and fitted to provide the parameters in Eq. (1) defining the many body non-additive polarisation contributions to the new ab inirzo VB NCC-like potential. Geometries are generated starting from the minimum energy configuration of the trimer by varying the relative 0-0distances of the three molecules, with the hydrogens moving rigidly and maintaining their relative positions with respect the plane of the oxygen’s. The calculations on the water dimer are performed at the VB level, with single and double excitations, according to the scheme already discussed. Inclusion of the intermolecular correlation is thus guaranteed. The 225 points of the PES computed for the water dmer included 13 configurations corresponding to the “closed” or “cyclic” structures. The energies corresponding to other 22 “open” configurations, near the minimum, were also calculated. It is to be noted that, although the “bifkrcated” structures are attractive, they do not correspond to a local minimum, as sometimes reported (Muggiest, Robins and Bassez-Muguet, 1991).

VI.Molecular dynamics simulation The molecular dynamics simulation is performed using the MOTECC suite of programs (Sciortino, Corongiu and Clementi1994) in the context of microcanonical statistical ensemble. The system considered is a cube with periodic boundary conditions, which contains 343 water molecules. Compatibility of this data with the water experimental density of 0.998 g/cm’ requires a cube with a side length of 21.7446 A. In accordance with the polarizable model, a spherical cutoff with radium equal to half of the simulation cube side is imposed, together with a switching knction to suppress energy drift.

M. Raimondi eta/.

274

The interaction of a single molecule with molecules whose oxygen atoms are outside the cutoff sphere are not calculated directly from Eq. (l), but by means of Ewald sum for Coulomb interactions. The equations of motion are integrated with a Gear prediction-correction algorithm at the sixth order. During the initial equilibration time, velocities are scaled every ten steps in order to klfil the relation: l N 3 < - C m , v : >=-Nk,T 2 i 2 (19) where T is the temperature imposed on the system (this relation is used for all translational, rotational and vibrational degrees of freedom, so that there will be a temperature relative to each degree of freedom). s long. In fact, the prediction-correction algorithm is Each step is 2.0 more accurate for small time steps, and moreover, molecular vibrational frequencies are above 1014s-'. The equilibration period is 25000 time steps, that is 5ps long. During the following 10 ps, corresponding to 50000 steps, atomic position and velocities and dipole moments of the molecules are evaluated every 10 steps, and stored for the subsequent analysis of the system properties.

VII. Simulationresults During the 10 ps of simulation time, temperature was equal to 299k5.95 K. The potential energy of the system of -40.56fo.22 kJ/mol is to be compared with the experimental value for the vaporisation enthalpy of liquid water at 298 K, which is -41.5 kJ/mol. The computed molecular dipole moment in the liquid is 2.48fo.01 D, to be compared with the experimental value of 2.4-2.6, see Table IV, and a value of 1.83D predicted by the model potential for the isolated molecule. The contributions to the total energy of the molecule obtained during the simulation are shown in Table 111. TABLE ILI Contributionsi to the total potential energy of the liquid water I ~ ( ~ / m o lI ) I VKlN I 11.19M.22 I VTWGBODY VPOL

I VTOT

-36.01M.24 -9.48M.20 I -29.378M.010I

Ab lnitioVB Interaction Potential of Liquid Water

275

It should be noticed that, in contrast with previous simulation results by Niesar, Corongiu, Clementi, Kneller and Bhattacharya, 1990 - NCC-flex -, the main contribution to the total potential energy comes from the two-body term, the many-body polarisation term contributing to the total potential energy only for the 23%. Some of the properties calculated during the simulation are reported in Table IV.

TABLE IV

(a) Svishchev, Kusakik, Wang and Boyd (1 996); (b) Ahlstrom, Wallqvist, Engstrom and Jonsson (1 989); (c) Berendsen, Grigera and Straatsma (1987); (d) Billeter, King and Van Gunsteren (1 994); (e) Wallqvist and Berne (1993); ( f ) Rahman and Stilliger (1971); (g) Lie and Clementi (1986); (h) Corongiu (1 992); (i) Niesar, Corongiu, Clementi, Kneller and Bhattacharya, (1 990); (i) Jorgensen, Chandrasekhar,Madura, Impey and Klein (1983); (k) Kell, G. S. (1971). As it can be seen in Table IV, all the results of the simulation computed using the present ab initio VB potential turn out encouraging. Even the computed internal pressure, although incompatible with the experimental data, is greatly improved

276

M. Raimondi eta/

with respect to the corresponding MCYL (two body only) or NCC-flex (manybody) results obtained by Lie and Clementi (1986) and by Corongiu (1992). Addition of the many body terms gives a computed internal pressure of about 700 atm to be compared with a value of about 500 atm obtained with the two body only potential. The fact that the calculated values of the internal pressure can change much after small changes in the potential surface has been the subject of a thorough discussion by Niesar, Corongiu, Huang, Dupuis and Clementi, (1989); the variation of about 200 atm obtained in the present work when including the effect of many-body polarisation interactions, can be assumed as a very encouraging result confirming the consistency of the new ah initio VB potential. In Table V the calculated geometry of the water molecule in liquid and gas phase is reported.

TABLE V Calculated geometric parameters of water molecule in gas and liquid phase compared with tl i experimental values. G a s phase Liquid phase I I HOH HOH I Present work 10.9576 1104.59 0.972 1102.75 11.52 Exp. (a,b) 10.9572 [ 104.52 0.966k0.006 I(102.8) (c) I1.510k0.01 (a) Benedict, Gailar and Plyler (1956); (b) Thiessen and Narten (1982); (c) Calculated from and . Fig. 3-6 and Tables VI-VIII show the experimental and computed site-site pair correlation hnctions m ( r ) , g&r) and m ( r ) . The relative position and intensities of peaks are reported in Tables VI-VIII.

Ab lnifioVB Interaction Potential of Liquid Water

Fig. 3. Site-site pair correlation functions g&).

Fig. 4. Site-site pair correlation functions gOH(r).

277

M. Raimondi eta/.

278

Fig. 5. Site-site pair correlation functions gIIH(r). TABLE VI Position and intensities (in parenthesis) of the peaks of g&). 1'Max I l'min I 2ndMax 1

I

Present work Exp.(neutron die.) (a) Exp.(Xray die.) (b)

(A)

2.77 (2.93) 2.89 (3.09) 2.83 (2.31)

(Q

3.68 (0.88) 3.32 (0.73) 3.45 (0.85)

(A)

4.77 (1.04) 4.53 (1.14) 4.53 (1.12)

2ndrnin

(A) 5.6 (0.88) 5.6 (0.86)

TABLE W Position and intensities (in parenthesis) of the peaks of go&). 1 1'Max I 1'min I 2"dMax I 2"dmin

Present work Exp.(neutrondie.) (a) Exp.(X ray d i e . ) (b)

(A)

1.83 (1.37) 1.85 (1.38) 1.90 (0.80)

(b) Narten and Levy (1971).

(ri>

2.40 (0.33) 2.35 (0.26) 2.45 (0.50)

(A)

3.24 (1.55) 3.35 (1.60) 3.35 (1.70)

1

(A)

4.23 (0.92) 4.25 (0.95)

-

I

279

Ab lnitioVB Interaction Potential of Liquid Water

TABLE VIIl 1' Max

lamin

Znd Max

2ndmin

(A)

(A)

(A)

(A) -

Present work 2.40 (1.40) Exp.(neutrondifi.) (a) 2.45 (1.26) 2 35 ( I 04) Exo (X rav d i e ) (b1 (a) Soper and Phillips (1986); (b) Narten and Levy ( 197 1).

3.03 (0.83) 3.05 (0.76) 3 00 (047)

3.81 (1.13) 3.85 (1.17) 4 00 (1 081

5.45 (0.96)

-

The number of the first neighbours around a central water molecule, calculated by integration of the correlation hnction gm(r) up to the position of the first minimum, is 5 . 5 . This result is consistent with a local tetrahedron with addition of contributions from some interstitial coordination. The computed W ( r ) shows a long range structure which decays faster than the experimentally based hnction. In particular, the second peak results lower. This can be ascribed to a variety of factors, including an incomplete determination of the two body potential, due to basis and configuration spaces not large enough. In addition, it is possible that the many-body contribution to the potential, approximated with a three-body term only, could be responsible of this loss in structure for the simulated liquid. The correlation functions go&) and m ( r ) are in good accordance with experimental data. By integrating go&) up to the first minimum a value of 1.92 is obtained. This means that each water molecule forms 3.84 hydrogen bonds with its first neighbours, allowing to conclude that at room temperature, the four nearest water molecules are preferentially oriented in a tetrahedron configuration. The comparison between the simulation which includes only two-body effects and the one obtained including many body terms (see Table IV) shows a perfect equivalence as regards the calculated properties of the liquid phase. The many body contribution plays a hndamental role only in the position of the second peak of the site-site pair correlation hnction gm(r), which shows a greater accordance with the experimental. VIII. IR and Raman spectra

IR and Raman spectra are calculated by means of autocorrelation hnction of the dipole moment and of the polarizability tensor of the system. The characteristics of the spectra, shown in Fig. 6 and 7, are reported in Table IX.

M. Raimondi eta/.

280

Vibrational band (cm-')

IRandRam Infia-Red

TABLE M. I band of liquid ater. Infia-Red High Frequency Raman Present work EXP. (a)

Hindered translation v T ~ Hindered translation VT

High Frequency

Raman

Present work

Weak band

v= 65 Band appears as shoulder of vL band

Band appears as shoulder of vL band

Very broad band extending from 300 cm-'to 900 cm-'with principal maximum at:

Very broad band extending from 300 cm-' to 900 cm-' with principal maximum at:

v= 193

Libration VL

1

v= 685

v= 190

vz 651

H-0-H bending vz

v= 1645

v= 1726

Association V A

vz 2125

v= 2185

0-H stretching vs Very broad band

Very broad band with two principal with two principal maxima and a maxima: shoulder: v= 3615 vz 3280 vz 3746

v= 3490 v= 3920 (a) Eisenberg and auzmann (1969:

v= 1640

vz 1726

Very broad band with principal maximum at:

Very broad bandwith principal maximum at:

v= 3439

v= 3583

281

Ab lnitioVB Interaction Potential of Liquid Water

ntensity (arb. units)

1

'oooo

0 00

1000 00

2000 00

3000 00

4000 00

v (cm-') Fig. 6. Calculated IR spectra.

[ntensity (arb. units)

1

Fig. 7. Calculated Raman spectra

The shift in the intramolecular frequencies on passing from gaseous to liquid phase are reported in Table X.

M. Raimondi eta/.

282

TABLE X

-167

-23 1 -209

(a) Eisenberg and Kauzmann (1 969).

M.Conclusions On the basis of the results of the dynamics simulation of liquid water it turns out that the VB like wavefunction based on the SCF-MI non-orthogonal occupied and virtual orbitals, describes accurately the intermolecular potential of water. Exclusion of BSSE in an apriori fashion is ensured and geometry relaxation effects are naturally taken into account. The new ab initio results were used to determine a BSSE free water-water interaction potential which was employed in molecular dynamics simulation of the liquid phase at room temperature. Extensive calculation were performed on water dimer and trimer and the parametrization of a new NCC-like (Niesar, Corongiu, Clementi, Kneller and Bhattacharya, 1990) potential was accomplished. The results showed good accordance with the experimental data relative to radial distribution functions, thermodynamic properties and geometric parameters. The computed IR and Raman frequencies shiftsalso agree with the experimental values. The prediction of the internal pressure shows a great improvement in the computed value and an increased stability. It should be noticed that although non-additive many-body contributions are likely to be necessary for an accurate reproduction of the experimental data, two body terms appear strongly dominant. The aim of the present work was to perform a first test of the reliability of the general VB approach. The next step will be to investigate the effect of the temperature and system density variations. In particular we want to set a challenge to the present procedure to correctly describe the physical behaviour of water at the critical point, where a fundamental variation of the hydrogen bond is observed (Soper, Bruni and Ricci, 1997). Preliminary results already obtained are encouraging. Future improvements, aimed to establish the actual importance of manybody contributions, will concentrate on the extension of the basis set employed.

Ab lnitioVB Interaction Potential of Liquid Water

283

More configurations in the VB step will also be added and improvements of the simulation attempted. References Ahlstrom, P., Wallqvist, A,, Engstrom, S., and Jonsson, B. (1989). Mol. Phys. 68, 563. Bartlett, R. J., Shavitt, I., Purvis 111, G. D. (1979). J. Chem. Phys. 71,28 1. Benedict, W. S., Gailar, N., and Plyler, E. K. (1956). J. Chem. Phys. 24, 1139. Berendsen, H. J. C., Grigera, J. R., and Straatsma, T. P. (1987). J. Phys. Chem., 91, 6296. Berens, P. H., and Wilson, K. R. (1981). J. Chem. Phys. 74,4872. Berens, P. H., Mackay, D. H. J., White, G. M., and Wilson, K. R. (1983). J. Chem. Phys. 79,2375. Bernal, J. D., and Fowler, F. D. (1933). J. Chem. Phys. 1, 515. Billeter, S. R., King, P. M., and van Gunsteren, W. F. (1994). J. Chem. Phys. 100, 6692. Bottcher, C. F. J. (1973). “Theory of electric polarisation”. Elsevier, Amsterdam. VOl. 1 . Boys, S. F., and Bernardi, F. (1970). Mol.Phys. 19, 553. Clarke, N., Raimondi, M., Sironi, M., Gerratt, J., and Cooper, D. L. Theor. Chim. Acta in the press Clementi, E., and Corongiu, G. (1983). Int. J. Quantum Chem Symp. 10, 3 1 . Cooper, D. L., Gerratt, J., and Raimondi, M. (1987). Adv. Chem. Phys. 69, 3 19. Corongiu, G. (1 992). Znt. J. Quantum Chem. 44, 1209. Curtiss, L. A,, Frurip, D. J., and Blander, M. (1979). J. Chem. Phys. 71,2703. Dyke, T. R., Mack,K. M., and Muenter, J. S. (1977) J. Chem. Phys. 66,498. Eisenberg, D., and Kauzmann, W. (1969). “The structure and properties of water”. Oxford University Press, London. Famulari, A. (1997). Ph.D. Thesis, Universita’ di Milano. Famulari, A,, Gianinetti, E., Raimondi, M., and Sironi, M. (1997a). Submitted. Famulari, A., Gianinetti, E., Raimondi, M., and Sironi, M. (1997b). Submitted. Gianinetti, E., Raimondi, M., and Tornaghi, E. (1996). Znt. J. Quantum Chem. 60, 157. Gianinetti, E., Vandoni, I., Famulari, A,, and Raimondi, M. Adv. Quantum Chem. Present volume. Habitz, P., Bagus, P., Siegbahn, P., and Clementi, E. (1983) Znt. J. Quantum Chem. 23, 1803. Jorgensen, W. L., Chandrasekhar, J., Madura, J., Impey, R.W., and Hein, M.L. (1983). J. Chem. Phys.79,926. Kell, G. S. (1971). ‘Water: a comprehensive treatise”. Vol. 1, 363. Franks, F., Ed. Plenum, New York.

204

M.Raimondi eta/.

Lie, G. C., and Clementi, E. (1986). Phys. Rev. A 33,2679. Mas, E. M., and Szalewicz, K. (1996). J. Chem. Phys. 104,7606. Millot, C., and Stone, A. J. (1992). Mol. Phys. 77,439. Muguet, F. F., Robinson, G. W., and Bassez-Muguet, M. P. (1991). Int. J. Quantum Chem. 39,449. Narten, A. H., and Levy, H. A. (1971). J. Chem. Phys. 55,2263. Niesar, U . , Corongiu, G., Clementi, E., Kneller, G. R., and Bhattacharya, D. K. (1990). J. Phys. Chem. 94,7949. Niesar, U., Corongk, G., Huang, M-J., Dupuis, M., and Clementi, E., (1989). Int. J. Quantum Chem. S’p. 23,42 1. Rahrnan, A,, and Stilliger, F. H. (1971). J. Chem. Phys. 55, 3336. Raimondi, M., Campion, W.,and Karplus, M. (1 977). Mol. Phys. 34, 1483. Raimondi, M., Sironi, M., Gerratt, J., and Cooper, D. L. (1 996) Int. J. Quantum Chem. 60,225. Schmidt, M. W., Baldridge, K. K., Boatz, J. A., Elbert, S. T., Gordon, M. S., Jensen, J. H., Koseki, S., Matsunaga, N., Nguyen, K. A,, Su, S. J., Windus, T. L., Dupuis, M., Montgomery, J. A. (1993). J. Comput. Chem. 14, 1347. Sciortino, F., and Corongiu, G. (1994). “Modern Techniques in Computational Chemistry: MOTECC-94.”. Clementi, E., Ed. ESCOM, Leiden, The Netherlands. Soper, A. K., and Phillips, M.G. (1986). Chem. Phys.107,47. Soper, A. K., Bruni, F., and Ricci, M. A. (1997). J. Chem. Phys. 106,247. Specchio, R., (1997). Degree Thesis, Universita’ di Milano. Svishchev, I. M., Kusakik, P. G., Wang, J., and Boyd, R. J. (1996). J. Chem. Phys., 105,4742. Texeira, J., Bellisent-Funel, M. C., Chen, S. H., and Dorner, B. (1985). Phys. Rev. Lett. 54,268 1. Thiessen, W. E., and Narten, A. H., (1 982). J. Chem. Phys., 77,2656. van Duijneveldt, F. B., van Duijneveldt-van de Rijdt, J. G. C. M., and van Lenthe, J. H. (1994). Chem. Rev. 94, 1873. van Duijneveldt-van de Rijdt, J. G. C. M., and van Duijneveldt, F. B. (1992). J. Chem. Phys. 97, 5019. Wallqvist, A,, and Berne, B. J. (1993). J. Phys. Chem. 97, 13841.