Diffusion Monte Carlo studies of isotope-substituted water trimers

20 December 1996

CHEMICAL PHYSICS LETTERS

ELSEVIER

Chemical Physics Letters 263 (1996) 680-686

Diffusion Monte Carlo studies of isotope-substituted water trimers Jon M. Sorenson a,1, Jonathon K. Gregory a, David C. Clary a Department of Chemistry, Unioersity of Cambridge, Cambridge CB2 1EW, UK b Department of Chemistry, Unioersity College London, London WC1H OAJ, UK Received 6 August 1996; in final form 21 October 1996

Abstract

We report the ground-state properties of several partially deuterated water trimers calculated with the rigid-body diffusion Monte Carlo method. Rotational constants are compared with recent experiments and good agreement is found. Tunneling splittings for the hydrogen flip in (HDO) 3 are predicted. New results are predicted for the experimentally unexamined mixed trimers.

I. Introduction Water clusters have been of considerable theoretical and experimental interest recently. Studies of small water clusters offer an opportunity to rigorously test theory against experiment and obtain a thorough understanding of the intermolecular interactions in water. The recent high-resolution far-infrared vibrational rotational tunneling (FIRVRT) spectroscopy investigations of Saykally and co-workers [1-5] have uncovered rich features in the far-infrared spectra of water clusters corresponding to complex tunneling splittings. The observation of tunneling splittings in water clusters has prompted several theoretical studies. Estimates for tunneling splittings in the water dimer were given by Althorpe and Clary by diago-

i Present address: Department of Chemistry, University of California, Berkeley, CA 94720-1460, USA.

nalizing the Hamiltonian for the ASP potential energy surface using a modified discrete-variable representation [6,7]. Wales estimated splittings in the trimer using an eigenvector-following approach to locate transition states and characterize reaction paths, allowing a simple approximation to provide an order of magnitude estimate [8]. Gregory and Clary have previously applied diffusion Monte Carlo to calculate splittings in the water dimer and trimer and found quite good agreement with experiment [9]. To complement the accumulating spectroscopic information on water clusters, high-resolution studies of partially deuterated water trimers have now been conducted [10]. Observation and prediction of the structure and tunneling splittings in such mixed trimer systems provides an excellent test of our understanding of the interactions present in the water trimer. The introduction of deuterium into the pure water trimer lowers the symmetry of the cluster and removes some of the degeneracies responsible for the observed splittings. In molecular symmetry terms,

0009-2614/96/$12.00 Copyright © 1996 Elsevier Science B.V. All rights reserved PII S 0 0 0 9 - 2 6 1 4 ( 9 6 ) 0 1 2 6 4 - X

J.M. Sorenson et a l . / Chemical Physics Letters 263 (1996) 680-686

the order of the corresponding molecular symmetry (MS) group [11] is reduced, The resulting effects on tunneling splittings will be determined by how much symmetry remains in the trimer after deuterium substitution. Study of these effects can only be undertaken by a theoretical method which calculates vibrational properties. Electronic structure calculations produce equilibrium geometries and are invariant to isotope substitution. As such, they can provide little insight into the dynamical properties of mixed trimers. It is possible to estimate zero-point energies from equilibrium geometries with a harmonic approximation; however, weakly-bound clusters are often bound by very anharmonic potentials, and the resulting estimates are not very reliable. Given an accurate potential energy surface, diffusion Monte Carlo (DMC) provides a way around this limitation [ 12,13]. The application of diffusion Monte Carlo techniques to vibrational calculations has recently allowed accurate calculations of the groundstate properties of clusters, and to some extent, excited-state properties [14-16]. The chief requirement for a DMC investigation is an accurate potential for a quantum simulation. With the development of intermolecular perturbation theory and improved supermolecular approaches, reasonably accurate ab initit potentials are becoming more widely available for a variety of systems. Inversion of spectroscopic information of small clusters is also providing a rich source of potential energy surfaces [ 17]. The work described here is an extension of our earlier DMC water trimer studies [9,18]. In the first section we describe the new potentials used in this study. The resulting rotational constants and zeropoint energies are presented in the next section. We also predict the tunneling splitting due to a hydrogen flip in (HDO) 3 there. The end section concludes the results of this work and further addresses the question of calculating tunneling splittings in the mixed trimers. 2. Theory

2.1. Diffusion Monte Carlo The DMC method has been successfully applied to calculations of vibrationally-averaged structures of

681

various van der Waals clusters. The molecular vibrations are simulated by approximating the ground-state nuclear wavefunction on a potential energy surface describing the electronic energy. The SchrSdinger equation for an N-body nuclear time-dependent system (in atomic units) is N 1 E 2m---'-~Vi2~b+ V~b,

O~

i--=0t

j=l

(l)

where t is the time, mj is the mass of atom j, V is the electronic potential energy surface, and ~b is the wavefunction. With the substitution r = it this equation reduces to ()lgt

N

Or

=

1 2m; ~2~b- VqJ,

(2)

which is a diffusion equation with diffusion terms and a source/sink term. This can be stochastically simulated by a random walk with a branching process. Details of the DMC procedure have been described in detail elsewhere [9,19,13]. We use the rigid-body formalism introduced by Buch [20,21] to separate intermolecular and intramolecular vibrations. Rigid-body DMC (RBDMC) allows a longer time step in the simulations and has been previously shown to significantly decrease statistical error in DMC studies of clusters [9,21,22]. Properties such as bond lengths and rotational constants are calculated with the descendant weighting method described by Kalos [23].

2.2. Potential energy surfaces The present calculations employ two recently developed potentials to describe the water-water interaction. Both potentials were developed from ab initio monomer properties using intermolecular perturbation theory. Intermolecular perturbation theory (IMPT) allows division of the intermolecular interaction energy into terms corresponding to electrostatic, induction, repulsion, and dispersion energies. To reproduce the energies of clusters with more than two molecules, non-pairwise corrections can also be added into these terms. While the electrostatic energy is pairwise additive, the remaining three terms have to be modified to account for many-body forces.

682

J.M. Sorenson et al. / Chemical Physics Letters 263 (1996) 680-686

The inclusion of three-body terms was previously found to have a significant effect on the calculated structures of water trimers [18]. The ASP-3B W2 potential used here is a recent improvement of the ASP-3B potential developed previously [18]. The ASP-3B potential is a three-body modification of the original Millot and Stone ASP-P potential [24] developed to model the water dimer. In the ASP-3B potential energy surface, a three-body induction energy and a three-body Axilrod-Teller dispersion energy were added to the ASP interaction, but short-range three-body terms were neglected. The ASP-3B W2 potential includes all of the terms present in ASP-3B with a different parametrization of the dispersion energy and damping on the induction energy. A charge-transfer term has also been added, and the distributed multipoles have been revised from the original MP2 values with multi-reference CI calculations [25]. For comparison, each system studied with the ASP-3B W2 potential was also examined with the NEMO potential developed by .~,strand et al. [26]. The NEMO potential has recently been shown to reproduce both ab initio and experimental water rimer properties well [27]. The potential was originally developed for molecular dynamics studies, and consequently is much faster to compute than the ASP potentials. Many-body effects are included with the incorporation of a three-body induction energy.

3. Results and discussion

We examined the structure of the following six mixed water trimers with RBDMC simulations using the ASP-3B W2 and NEMO potentials: (HDO)(n20) 2, (D2OXH20) 2, (HDO)2(H20), (HDO) 3, (020)2(H20), and (HDOXD20) 2. The two pure water trimers, (D20) 3 and (H20)3, were also studied for comparison with these new potentials. In the case for the mixed trimers where more than one structure is possible, only the structure with the most deuteriums bound in the plane of the oxygens was considered. This choice is justified in the next section.

3.1. Rotational constants For direct comparison with experiments, expectation values of rotational constants were calculated. For rotational constants A and B we report here the average value, (A + B)/2. For most configurations visited in the RBDMC random walk, A and B are distinct. A calculation of the expectation values of A and B using instantaneous inertial axes will result in distinct values for A and B. It is experimentally observed [2] however, that vibrational averaging results in a symmetric top structure (A = B) for the pure water trimers. We have found that (A + B ) / 2 provides a better comparison to the experimentallydetermined rotational constants. The question of how best to calculate rotational constants in non-rigid complexes has been addressed recently by Emesti and Hutson [28]. They conclude that calculations of expectation values of rotational constants should use axes which satisfy the Eckart conditions to best separate vibrational and rotational motion. However, in the case of very large amplitude motions, as we have here in the case of a hydrogen flip, there is no definitive procedure for choosing an axis frame [29]. In the absence of a clear direction to proceed, comparison of (A + B ) / 2 with experiment is most appropriate. This choice has also been used, without discussion, in many previous DMC studies of vdW complexes [15,21,30]. Results for rotational constants calculated at the minimum energy geometry are shown in Table 1, while vibrationally-averaged rotational constants, calculated using the RBDMC method, are presented in Table 2. It is readily apparent that the ASP-3B W2 potential produces excellent agreement with the experimental rotational constants. For (D20) 3, the rotational constants are predicted within 0.2%. We would expect that such good agreement is fortuitous, as we are working within the confines of an approximate rigid-body potential. The agreement for (H20) 3 is not quite so good ( = 2%), but still better than the expected statistical error in diffusion Monte Carlo results. The ASP-3B W2 expectation value for the O - O separation was found to be 2.84 .~ in (020) 3 and 2.85 ,~ in (H20)3, also in excellent accord with the values of 2.837 and 2.840 .~, respectively, found in Liu et al.'s fit to experimental data [10]. The NEMO potential shows less good agreement with

683

J.M. Sorenson et al. / Chemical Physics Letters 263 (1996) 680-686

Table 1 Equilibrium rotational constants (MHz) for the pure and mixed water trimers (A + B ) / 2 a

(H 20)3 (HDOXH20) 2 (D2OXH 20)2 (HDO)2(H 2O) (D20)2(H 20) (HDO) 3 (HDOXD 20) 2 (D20) 3

C

ASP-3B W2 b

NEMO c

ASP-3B W2

NEMO

6807 6729 6495 6660 6212 6584 6132 5915

6718 6654 6421 6586 6134 6513 6072 5855

3477 3437 3324 3400 3196 3362 3163 3063

3433 3398 3286 3363 3160 3327 3133 3033

a See explanation in text. b Ref. [25]. c Refers to the NEMO 5-site water-water potential [26].

e x p e r i m e n t ; w e c a n trace this to the l a r g e r p r e d i c t e d O - O s e p a r a t i o n o f 2 . 9 4 .~ in ( H 2 0 ) 3 a n d 2.93 ~, in ( D 2 0 ) 3 in t h a t potential. I n p r e l i m i n a r y s t u d i e s w i t h the m i x e d t r i m e r s , it w a s f o u n d t h a t d e u t e r i u m prefers to r e m a i n b o u n d in the h y d r o g e n b o n d , and, for the a b o v e results, the g r o u n d state o f e a c h s y s t e m w a s s i m u l a t e d w i t h as m a n y d e u t e r i u m s i n v o l v e d in h y d r o g e n b o n d s as p o s s i b l e . T h i s r e s u l t is e x p e c t e d f r o m a c o n s i d e r a t i o n of zero-point energies. Deuterium bonds have lower zero-point energies than hydrogen bonds and deuterium-bonded trimers should be correspondingly m o r e stable. W e h a v e c o n d u c t e d R B D M C s i m u l a t i o n s o f ( H D O ) 3 w i t h a d e u t e r i u m in or o u t o f the

p l a n e o f the o x y g e n s a n d f o u n d that the d e u t e r i u m b o n d e n h a n c e s the stability o f the t r i m e r b y a b o u t 5 0 c m - l in the A S P - 3 B W 2 potential. T h i s is n o t a l a r g e f a c t o r b u t it b e c o m e s m o r e s i g n i f i c a n t at the l o w t e m p e r a t u r e s o f e x p e r i m e n t s , w h e r e k T -- 4 c m - l [31 ]. S p e c t r o s c o p i c studies o f i s o t o p e - s u b s t i t u t e d w a ter d i m e r s h a v e also f o u n d t h a t d e u t e r i u m - b o n d e d s t r u c t u r e s are m o r e stable t h a n t h e i r h y d r o g e n - b o n d e d c o u n t e r p a r t s [32]. 3.2. Z e r o - p o i n t energies

C a l c u l a t i o n o f the e x p e c t a t i o n v a l u e o f the e n e r g y in D M C s i m u l a t i o n s a l l o w s u s to c a l c u l a t e the zero-

Table 2 Vibrationally-averaged rotational constants (MHz) for the pure and mixed water trimers. Statistical uncertainties are 0.1% for the RBDMC results (A + B ) / 2 a

(H 20)3 (HDO)(H 20)2 (D2OXH20) 2 (HDO)2(H 20) (D20)2(H zO) (HDO) 3 (HDOXD~ 0) 2 (D20) 3

C

ASP-3B W2 b

NEMO c

experiment d

ASP-3B W2

NEMO

experiment

6557 6501 6284 6470 6071 6402 6025 5805

6228 6197 5985 6152 5728 6074 5673 5494

6647 -

3349 3304 3193 3289 3118 3256 3109 2998

3127 3115 3013 3091 2893 3055 2872 2791

3088 e

6015 5796

a See explanation in text. b Ref. [25]. c Refers to the NEMO 5-site water-water potential [26]. a Ref. [10]. Liu et al. [10] comment that this number is high due to nonstructural Coriolis contamination. With a Monte Carlo fit, they find a value of C = 2999 MHz, in excellent agreement with the ASP-3B W2 result.

J.M. Sorenson et al. / Chemical Physics Letters 263 (1996) 680-686

684 2100 2000

Q,

1900

',9

'- ;

Fig. 2. The h y d r o g e n flip r e a r r a n g e m e n t in water ~ m e r s . The middle picture is the transition state.

180~ 1700

than (D20)2(H20), it is slightly less stable because it has one less deuterium.

1600 1500

(HDO)(H~O)2 (HDO)~(H~O) (HDO)3 (D~O)~ (H20)$ (D~O)(H~O)~ (D~O)~(H~O) (HDO)(D20)2 Fig. I. ASP-3B W2 zero-point energies for the pure and mixed trimers ordered b y the n u m b e r o f deuterium bonds: n o b o n d s ( O ) , one b o n d ( × ) , two b o n d s ( + ) , three b o n d s ( * ) . Energies are in e r a - ].

point energies stabilizing the trimers. The resulting bond dissociation energies and zero-point energies are tabulated in Table 3. Because we are employing the rigid-body approximation, the resulting zero-point energies do not take account of the zero-point energies present in the O - H bond. However, we can compare the relative stability of the trimers against each other. The stabilization in the mixed trimers is most dependent on the number of deuterium bonds in the complex; more deuterium bonds leads to a more stable complex. The number of unbound deuterium atoms is a secondary factor, with more deuteriums stabilizing the complex further. These effects can best be seen in Fig. 1. It is interesting to note that while (HDO) 3 possesses more deuterium bonds

3.3. Tunneling splittings In general, tunneling splittings are much harder to calculate in the mixed trimer systems. Where, before, fixed nodes could be placed at transition states by employing symmetry considerations, most of these mixed trimers lack the symmetry necessary for these assumptions. As an example, consider the flip of a hydrogen in (HDO)(H20) 2 versus the flip of a hydrogen in (H20) 3 (see Fig. 2). In the case of the pure water trimer, the flip interchanges isoenergetic minima on the potential energy surface because of the indistinguishability of the three water molecules. In contrast, the introduction of the deuterium in the case of the mixed trimer provides each water with a distinct environment and allows us to distinguish between the three. The flip no longer interchanges isoenergetic minima, and, consequently, calculation of the tunneling splitting needs to proceed differently from the standard case of degenerate tunneling. Without a symmetric transition state between the

Table 3 V i b r a t i o n a l l y - a v e r a g e d ( D o) a n d r i g i d - b o d y zero-point energies (cm - j ) for the pure a n d m i x e d w a t e r trimers. The dissociation energy c o r r e s p o n d s to the dissociation o f the trimer into its three m o n o m e r s . Statistical uncertainties are 0 . 2 % Vibrationally a v e r a g e d

(H20) 3 ( H D O ) ( H 2O)2 (D20)(H 20) 2 ( H D O ) 2 ( H 2O) (D/O)2(H 20) (HDO) 3 (HDOXD20) 2 (D20) 3

Zero-point

ASP-3B W2 a

NEMO b

ASP-3B W2 ¢

NEMO d

-

-

2068 1979 1930 1868 1741 1762 1614 1566

1789 1702 1666 1611 1538 1527 1451 1413

3820 3910 3958 ~,020 4147 4126 4275 4322

2908 2996 3032 3087 3160 3170 3247 3284

Ref. [25]. b Ref. [26]. ¢ D e is - 5 8 8 8 c m - ' for the A S P - 3 B W 2 equilibrium g e o m e t r y . d De is -- 4698 c m - ~ for the N E M O equilibrium geometry.

J.M. Sorenson et al. / Chemical Physics Letters 263 (1996) 680-686

Table 4 Predicted tunneling splittings due to hydrogen flips (in cm-t). Statistical uncertainties are in parentheses

(H 20)3 (HDO) 3 (D20) 3

NEMO

ASP-3B W2

49(2) 45(3) 26(2)

16(4) 15(5) 4(3)

minima, it is not possible to a priori apply a fixed node. As a first estimate, we would generally expect the tunneling splitting to be decreased in the case of non-degenerate tunneling due to the smaller overlap between wavefunctions localized at local minima in the trimer potential energy surface. There are several sources of tunneling splittings in the water trimers; here we concentrate on the simplest splitting caused by the hydrogen flip. For the special case of (HDO) 3, with all three deuteriums bound in the OOO plane, flips do interchange isoenergetic minima, and we can calculate the expected tunneling splitting following the method detailed in Ref. [9]. The excited tunneling state is simulated by fixing a node in the wavefunction in the middle of the flipping coordinate of a hydrogen (or deuterium, as the case may be) not bound in the plane of the oxygens. Table 4 summarizes the predicted tunneling splittings. Because the flip in (HDO) 3 involves motion of a hydrogen, and not the heavier deuterium, the resulting tunneling is very similar to the case of (H20) 3, and we find that the magnitude of the splitting is statistically indistinguishable from that in (H20) 3. In the splitting pattern for (HDO)3, we would expect to observe the quartet resulting from the hydrogen flip [8,9] with no additional splittings due to the symmetry-disallowed bifurcation tunneling. Recent high-level ab initio calculations of the transition state for hydrogen flips have shown that the ASP-3B W2 barrier height of 403 cm-1 is too high [33], and we might expect that the predicted tunneling splitting is too small compared to experiment. In this case we might put more faith in the larger tunneling splittings predicted by the NEMO potential.

685

tunneling splittings in water trimers. The mixed trimers only exhibit a subset of the predicted tunneling splittings, and thus we can confidently fix the origins of each splitting. The case of (HDO) 3 is an excellent example: the only experimental splittings which should be observed will come about from the hydrogen flip. To isolate the splitting due only to bifurcation tunneling, the trimers (D20)(H20) 2 and (D20)2(H20) should be examined - where tunneling due to flips are now symmetry-disallowed. DMC calculations of this splitting are possible [9], but this splitting is small enough to require too many DMC simulations for such calculations to be feasible. Besides offering definitive answers for the origin of tunneling splittings in water trimers, studies of the mixed trimers also provide a good series of structures for testing potential energy surfaces. Here the ability of DMC to predict vibrationally-averaged structures is a crucial ingredient in the validation process. By predicting dynamic properties, DMC allows us to verify new potential energy surfaces with the complete set of possible water trimer isotopomers, not only the one pure water trimer normally examined for validation. Through this process, we have found here that the ASP-3B W2 surface provides rotational constants that compare well with experiment for several water trimers, and we have every reason to believe that future spectroscopic studies on the mixed trimers will further confirm the results reported here.

Acknowledgement We are grateful to Richard Saykally and Kun Liu for communicating their unpublished mixed trimer results and to Gunner Karlstrbm for making available his NEMO code. JMS acknowledges a scholarship from the Winston Churchill Foundation, and JKG acknowledges a grant from the Engineering and Physical Sciences Research Council (EPSRC).

References 4. Conclusions The decreased symmetry in the mixed trimers makes them good targets for the precise study of

[1] N. Pugliano and R.J. Saykally, Science 257 (1992) 1937. [2] K. Liu, J.G. Loeser, M.J. Elrod, B.C. Host, J.A. Rzepiela, N. Pugliano and R.J. Saykally, J. Am. Chem. Soc. 116 (1994) 3507.

686

J.M. Sorenson et a L / Chemical Physics Letters 263 (1996) 680-686

[3] N. Pugliano, J.D. Cruzan, J.G. Loeser and R.J. Saykally, J. Chem. Phys. 98 (1993) 6600. [4] R.J. Saykally, K. Liu and J.D. Cruzan, Science 271 (1996) 929. [5] K. Liu, M.G. Brown, C. Carter, R.J. Saykally, J.K. Gregory and D.C. Clary, Nature 381 (1996) 501. [6] S.C. Althorpe and D.C. Clary, J. Chem. Phys. 101 (1994) 3603. [7] S.C. Althorpe and D.C. Clary, J. Chem. Phys. 102 (1995) 4390. [8] D.J. Wales, J. Am. Chem. Soc. 115 (1993) 11180. [9] J.K. Gregory and D.C. Clary, J. Chem. Phys. 102 (1995) 7817. [10] K. Liu, M.G. Brown, M.R. Viant, J.D. Cruzan and R.J. Saykally, Mol. Phys., to be published. [11] PoR. Bunker, Molecular symmetry and spectroscopy (Academic Press, New York, 1979). [12] J.B. Anderson, J. Chem. Phys. 63 (1975) 1499. [13] M.A. Suhm and R.O. Watts, Phys. Rep. 204 (1991) 293. [14] R.E. Miller, D.F. Coker and R.O. Watts, J. Chem. Phys. 82 (1985) 3554. [15] H. Sun and R.O. Watts, J. Chem. Phys. 92 (1990) 603. [16] R.N. Barnett and K.B. Whaley, Phys. Rev. A 47 (1993) 4082. [17] A. van der Avoird, P.E.S. Wormer and R. Moszynski, Chem. Rev. 94 (1994) 1931.

[18] J.K. Gregory and D.C. Clary, J. Chem. Phys. 103 (1995) 8924. [19] J.B. Anderson, Int. Rev. Phys. Chem. 14 (1995) 85. [20] V. Buch, J. Chem. Phys. 97 (1992) 726. [21] P. Sandler, J.O. Jung, M.M. Szcz~niak and V. Buch, J. Chem. Phys. 101 (1994) 1378. [22] J.K. Gregory and D.C. Clary, Chem. Phys. Lett. 228 (1994) 547. [23] M.H. Kalos, Phys. Rev. A 2 (1970) 250. [24] C. Millot and A.J. Stone, Mol. Phys. 77 (1992) 439. [25] C. Millot and A.J. Stone, to be published. [26] P.-O. ,~strand, A. Wallqvist and G. Karlstr~m, J. Chem. Phys. 100 (1994) 1262. [27] P.-O. ,~strand, P. Linse and G. Karlstr~m, Chem. Phys. 191 (1995) 195. [28] A. Emesti and J.M. Hutson, Chem. Phys. Lett. 222 (1994) 257. [29] J.M. Hutson, Mol. Phys. 84 (1995) 185. [30] K.A. Franken and C.E. Dykstra, J. Chem. Phys. 100 (1994) 2865. [31] R.J. Saykally, Acc. Chem. Res. 22 (1989) 295. [32] E.N. Karyakin, G.T. Fraser, F.J. Lovas, R.D. Suenram and M. Fujitake, J. Chem. Phys. 102 (1995) 1114. [33] T.R. Walsh and D.J. Wales, J. Chem. Soc. Faraday Trans., to be published.