Flexible, ab initio potential, and dipole moment surfaces ... - AIP Journals

5 downloads 86 Views 429KB Size Report
3Georgia Gwinnett College, Lawrenceville, Georgia 30043, USA. 4Division of Physical and Chemical Sciences, International Atomic Energy Agency, Vienna, ...
THE JOURNAL OF CHEMICAL PHYSICS 134, 094509 (2011)

Flexible, ab initio potential, and dipole moment surfaces for water. I. Tests and applications for clusters up to the 22-mer Yimin Wang,1 Xinchuan Huang,2,a) Benjamin C. Shepler,3 Bastiaan J. Braams,4 and Joel M. Bowman1,b) 1

Cherry L. Emerson Center for Scientific Computation and Department of Chemistry, Emory University, Atlanta, Georgia 30322, USA 2 SETI Institute, 189 Bernardo Ave, Suite 100, Mountain View, California 94043,USA 3 Georgia Gwinnett College, Lawrenceville, Georgia 30043, USA 4 Division of Physical and Chemical Sciences, International Atomic Energy Agency, Vienna, Austria

(Received 29 December 2010; accepted 26 January 2011; published online 3 March 2011) We report full-dimensional, ab initio potential energy and dipole moment surfaces, denoted PES and DMS, respectively, for arbitrary numbers of water monomers. The PES is a sum of 1-, 2-, and 3-body potentials which can also be augmented by semiempirical long-range higher-body interactions. The 1-body potential is a spectroscopically accurate monomer potential, and the 2- and 3-body potentials are permutationally invariant fits to tens of thousands of CCSD(T)/aug-cc-pVTZ and MP2/aug-ccpVTZ electronic energies, respectively. The DMS is a sum of 1- and 2-body DMS, which are covariant fits to tens of thousands MP2/aug-cc-pVTZ dipole moment data. We present the details of these new 2- and 3-body potentials and then extensive applications and tests of this PES are made to the structures, classical binding energies, and harmonic frequencies of water clusters up to the 22-mer. In addition, we report the dipole moment for these clusters at various minima and compare the results against available and new ab initio calculations. © 2011 American Institute of Physics. [doi:10.1063/1.3554905] I. INTRODUCTION

The development of an accurate, flexible, ab initio-based potential for water has been a goal of our group for several years. The motivation for such a potential is self-evident, i.e., to describe the spectroscopy, dynamics, and statistical mechanics of water. However, the difficulty in achieving this goal is perhaps equally self-evident. The basic hurdle to overcome is the very high dimensionality of such a potential. Beyond, this is the need to address the small barriers that separate the many local minima of even modest sized water clusters, as well as the long-range electrostatic and electrodynamic, i.e., induced and long-range interactions. Clearly, this calls for a highly accurate and, therefore, “expensive” ab initio method, which simply adds to the difficulty in the achieving the goal. We have recently made good progress toward achieving this goal. Specifically, we reported a full-dimensional, permutationally invariant potential energy surface for the water dimer, which, by definition, contains the 1- and 2-body interactions. This potential, which we denoted HBB2, is a precise fit to roughly 30 000 CCSD(T)/aug-cc-pVTZ (aVTZ) electronic energies.1 The HBB2 PES has the essentially exact dissociation energy, De , of 4.98 kcal/mol and highly accurate barriers to isomerization of all the known low-lying stationary points of the water dimer. (It also accurately describes the high energy chemical exchange reaction of two H atoms.2 ) This potential has also be shown to give highly accurate vibration–rotation tunneling splittings and rotation cona) Present address: MS 245-6, NASA Ames Research Center, Moffett Field,

CA 94035.

b) Electronic mail: [email protected].

0021-9606/2011/134(9)/094509/12/$30.00

stants for (H2 O)2 and (D2 O)2 .3 In order to improve the longrange dipole–dipole interaction, we also considered a hybrid dimer potential that switches to the Thole-type model modication 3 flexible monomer (TTM3-F) potential of Xantheas and Fanourgakis.4 In addition, we reported a precise fit to the flexible dimer and monomer dipole moment surfaces, based on MP2/aVTZ calculations.3 More recently, we reported a full-dimensional (21 degree-of-freedom) 3-body potential, based on permutationally invariant fitting of roughly 30 000 MP2/aVTZ 3-body energies.5, 6 We denoted full-dimensional potentials built from the 1- and 2-body potentials and 1-, 2-, and 3-body potentials as PES(1,2) and PES(1,2,3), respectively. These PESs, and especially PES(1,2,3), were shown to predict the structures, energies, and harmonic normal mode frequencies of a number of water clusters up to the hexamer, with high accuracy. Here, we present new versions of these PESs. First, we represent the 2-body potential by a new permutationally invariant fit to intrinsic CCSD(T)/aVTZ 2-body energies. For the 1-body, i.e., monomer potential we use the spectroscopically accurate, potential of Partridge and Schwenke.7 The new 2-body potential is also switched in the long-range to TTM3F, as was done previously.1 For the 3-body part of PES(1,2,3) we refit the intrinsic 3-body energies in a fashion that facilitates the use of a numerically efficient cutoff function. Finally, we augment PES(1,2,3) by inclusion of long-range 4and higher-body polarization interactions, again relying on TTM3-F for these. After giving the details of these new PESs, we present applications and tests to the structures, energies, and harmonic frequencies of all water clusters including the dimer, trimer,

134, 094509-1

© 2011 American Institute of Physics

Wang et al.

J. Chem. Phys. 134, 094509 (2011)

and hexamer and then continuously from the decamer up to the 22-mer, including many isomers. These predictions are tested against previous and, in some cases, new direct ab initio calculations. It is worth noting that this strictly ab initio approach for a flexible potential stands apart from all other model potentials for water, which make extensive use of analytical expressions for the components of the potential. A review of many of these potentials can be found in the recent frontiers article by Szalewicz et al.8 and also in a recent papers by Fanourgakis and Xantheas4 and Jordan et al.9 (The so-called CC-pol potential10 also stands apart as ab initiobased; however, it is for rigid monomers. This potential is built from a high-level 2-body potential plus a 3-body potential, obtained at the Hartree–Fock level of theory.11 ) These components largely come from molecular mechanics, e.g., distributed point charges on the monomers, simple 2-body interactions and, in the more sophisticated models, polarization terms. Typically, the parameters in these analytical models are determined empirically by comparing against bulk properties using classical statistical mechanical calculations of these properties. However, the TTM3-F potential of Xantheas and co-workers, also incorporates ab initio data to determine some parameters. The potential is arguably the best flexible water potential of its type, and so we make comparisons with it here for larger clusters. Comparisons with TTM3-F for the water trimer and hexamer, with an earlier version of the type of PES we report here, were previously reported by Wang et al.5 The paper is organized as follows. Section II presents details of the new 2-body potential. Following that, two new fits to an enlarged database of 3-body MP2/aVTZ energies are described. Then we describe how 4- and higher-body longrange polarization energies are added to the 1-, 2- and 3body representation of the potential for arbitrary numbers of monomers. A number of comparisons of the 2- and 3-body fits are made with direct ab initio calculations as well as results from TTM3-F along cuts through the potential. The fit to the dimer dipole moment, obtained from MP2/aVTZ data, is briefly reviewed and the dipole moment for arbitrary monomers is described. Details of the covariant fitting of the dipole moment are given in the Appendix. Section III presents extensive applications and tests of the new potential and dipole moment surfaces. Section IV contains a summary, conclusions, and description of future work. II. POTENTIAL AND DIPOLE MOMENT SURFACES A. 2-body potential

The 2-body potential, denoted V2-body , is defined as V2-body (1, 2) = Vdimer (1, 2) − Vmonomer (1) − Vmonomer (2), (1) where 1 and 2 denote the collective Cartesian coordinates of monomer i, Vdimer is the electronic energy of two monomers at a given configuration, and Vmonomer are the two monomer electronic energies, respectively. Each of these electronic energies was obtained using the CCSD(T)/aVTZ method/basis, as implemented in MOLPRO suite of electronic structure

CCSD(T)/aVTZ 2−body fit TTM3−F

500 0 Energy (cm−1)

094509-2

−500

6

8

10 12 14 16 18

0

−1000 −50

−1500 −2000

−100 −150

2

4

6

8 10 RO−O (Å)

12

14

FIG. 1. One-dimensional cuts of indicated 2-body potentials as a function of the OO-distance with fixed monomer geometries at a fixed orientation corresponding to the dimer global minimum. The inset shows a zoom-in of the long-range region of 2-body potentials.

programs.12 The database of electronic energies consists of roughly 30 000 intrinsic 2-body electronic energies, distributed as described previously.3, 13 The electronic energies were calculated with and without the standard counterpoise correction and then the two databases were combined in a way that produced the exact dissociation energy, De , as described elsewhere.1 As in previous work,3, 13 this 12-dimensional 2-body potential is represented as an expansion of a complete basis of permutationally invariant polynomials of Morse-type functions, e−ri j /λ , where λ is fixed at 3.0 Bohrs and ri j is the internuclear distance of atom i and j, of maximum degree seven. A total number of 5 227 linear coefficients were determined in a least-squares fashion. Since we used the same permutationally invariant polynomial bases and the same set of configurations for the fitting of the intrinsic 2-body potential as for HBB2 dimer PES, the present 2-body fit inevitably inherits from HBB2 PES the limitation in describing the long-range behavior of the 2-body interaction, as discussed in detail elsewhere.1 Figure 1 shows plots of one-dimensional cuts, described in detail in the caption, of the intrinsic 2-body potential as a function of the OO-distance from present CCSD(T)/aVTZ electronic structure calculations, the new 2-body fit, and the TTM3-F potential. As seen, over the extended range from 2.5 to 14 Å all three cuts agree well with each other and smoothly approach the correct asymptotic value, i.e., zero, as expected. However, closer examination of cuts in the range of 6–20 Å does reveal that the 2-body fit starts to deviate from the CCSD(T)/aVTZ and TTM3-F potentials at around 7 Å, and approaches zero more rapidly than the other two potentials. (This behavior is primarily due to a paucity of electronic data in the long-range region and also the nature of the Morse variables used in the fit.) By contrast, the TTM3-F potential agrees better with the CCSD(T)/aVTZ energies in this long-range region of this cut. (It should be noted that in the short-range repulsive region of the potential the present 2-body fit agrees much better with the CCSD(T)/aVTZ than the TTM3-F potential; this shortcoming of the TTM3-F potential has recently been noted by others.14 ) To improve the

Water potential and dipole moment surfaces

J. Chem. Phys. 134, 094509 (2011)

V˜2-body (i, j) = (1 − s) V2-body (i, j) + s VTTM3-F (i, j),

(2)

where VTTM3-F (i, j) is the total interaction energy between monomer i and j given by the TTM3-F potential and s could be any a suitable function that varies smoothly between zero and one. As previously,1 the switching function is a simple polynomial of the OO-distance, ROO , given by     ROO − a 3 ROO − a 4 − 15 s(ROO ) = 10 b−a b−a 5  ROO − a +6 . (3) b−a As seen, the parameters a and b determine the switching region and the values chosen are 6.5 and 7.5 Å, respectively. Note among the components of the TTM3-F 2-body potential, at ROO = 5.5 Å the electrostatic component is essentially the −3 only remaining energy contribution, and this exhibits the ROO 1 dependence of the usual static dipole–dipole interaction. Obviously, there are many other analytical expressions for this interaction in the literature,4, 15–17 and they could also be reasonable alternatives for the long-range 2-body interaction of our water potential. It is also worth noting that switching to the long-range part of TTM3-F not only “corrects” the longrange behavior of the present 2-body potential, it also requires less cpu effort to evaluate, and this is an advantage for large water clusters. Next, we show a comparison of the present 2-body fit, TTM3-F and direct CCSD(T)/aVTZ energies as functions of the free and H-bonded OH-stretches of the donor H2 O in Fig. 2 and Fig. 3, respectively. (Note, at the global minimum of the dimer the H-bonded and free OH distances are 0.96 and 0.97 Å, respectively. And for the fundamental excitation of these modes, the relevant range is roughly 0.86 to 1.15 Å.) As seen, there is a much greater variation of the 2-body poCCSD(T)/aVTZ 2−body fit TTM3−F

−1500

−1

Energy (cm )

−1300

−1700

−1900 free O−H stretch −2100 0.4

0.6

0.8

1 1.2 RO−H (Å)

1.4

1.6

1.8

FIG. 2. One-dimensional cuts of indicated 2-body potentials as a function of the OH-H distance of a free monomer O–H stretch with remaining degrees of freedom fixed at a dimer geometry corresponding to the dimer global minimum.

CCSD(T)/aVTZ 2−body fit TTM3−F

−500 −2000 −1

long-range behavior of the 2-body potential, we take an entirely practical approach, that is, switching in the long-range to the TTM3-F potential via a simple function.1 Thus, the final expression for the 2-body potential, denoted V˜2-body , is given by

Energy (cm )

094509-3

−3500 −5000

H−bonded O−H stretch

−6500

0.4

0.6

0.8

1 1.2 RO−H (Å)

1.4

1.6

1.8

FIG. 3. One-dimensional cuts of indicated 2-body potentials as a function of the OH-distance of a H-bonded monomer O–H stretch with remaining degrees of freedom fixed at a dimer geometry corresponding to the dimer global minimum.

tential for the H-bonded OH-stretch than the free OH -stretch. This is expected because the free OH potential is represented better by the 1-body monomer potential than the H-bonded OH-potential. Also, the agreement of the present 2-body fit with direct CCSD(T)/aVTZ energies is quite good and much better than the results from TTM3-F, in contrast to the longrange behavior with ROO . B. 3-body potential

In previous work, we reported a 3-body fit that was the sum of two-atom, three-atom, four-atom, and one full nineatom terms.5 Each of these multi-atom terms is a complete invariant representation in terms of corresponding Morse-type variables up to certain maximum degree. In addition, since the database (at that time) did not contain an extensive number of configurations for large OO-distances, it was necessary to multiply each multi-atom term by a damping function of the internuclear distances, which enforces the physical boundary condition that V3-body vanishes as any monomer is removed from the trimer configuration. However, by introducing additional nonlinear parameters (in the damping functions) this representation of the full 3-body potential did contain small, artificial “bumps” in the near dissociation limit. Here, we present two new fits to an enlarged database of roughly 40 000 3-body MP2/aVTZ energies obtained with, MOLPRO ,12 as described previously,5 and denoted V3-body . The enlarged database includes many configurations between the complex region and the fragment region of the dimer plus monomer or three separated monomers, i.e., geometries of one or more OO-distances in the range from 6 to 13 Å. We also added more data for OOO angles in the range of 120– 180. The fits to these energies use a single 9-body expression, without damping, and (as shown below) these fits go smoothly to zero as the OO-distance gets large. The two fits correspond to total polynomial orders of 5 and 6. The lower order fit, denoted by 5th-order fit, has 1 380 unknowns in the least-squares system corresponding to invariant polynomials up to maximum degree five. The higher-order fit, denoted by 6th-order fit, is an expansion of invariant polynomials up to

094509-4

Wang et al.

-100

MP2/aVTZ 3−body fit 5 −order 3−body fit TTM3−F th 6 −order th

200 0

-200 -300

5

−200

7

9

Energy (cm-1)

−1

Energy (cm )

J. Chem. Phys. 134, 094509 (2011)

11

20 0

−400

-500 MP2/aVTZ 6th-order 3-body fit th 5 -order 3-body fit TTM3-F

-600

−20

−600

-400

−40

-700

−60

−800

2

4

6 8 RO−O (Å)

10

-800

12

FIG. 4. One-dimensional cuts of indicated 3-body potentials as a function of the OO-distance, for the removal of one monomer of the water trimer at a fixed orientation and geometry corresponding to the trimer global minimum. The inset shows a zoom-in of the asymptotic region of 3-body potential vanishing by 6 Å.

degree six and, thus, has more free coefficients, 5 849 to be exact. The 5th-order fit runs eight times faster than the 6thorder fit. The rms fitting error over the entire data set is 51 and 15 cm−1 for the 5th-order and 6th-order fit, respectively. The smallest cluster where 3-body effects are present is the trimer and we use that cluster for comparisons of the 3-body potential fits against direct MP2/aVTZ energies. First, in Fig. 4 the 3-body potential fits are shown as onedimensional functions in the OO-distance, compared against present MP2/aVTZ calculations. This cut is for the water trimer at the global minimum. As seen, both fits agree quite well with the ab initio energies over the whole range of OOdistances up to 12 Å. Also, from this figure it is clear that the 3-body potential is essentially zero by 6 Å. By contrast, as seen in Fig. 1, the 2-body potential is of the order of −10 cm−1 even at an OO-distance of 12 Å. Thus, on this basis we observe that the 3-body potential is shorter range than the 2-body one. The rapid approach to zero of the 3-body potential enables us to use a simple cutoff function for it given by, V˜3-body (i, j, k) = (1 − s) V3-body (i, j, k)

60

80

100 120 140 θO-O-O (degree)

160

180

FIG. 5. One-dimensional cuts of indicated 3-body potentials as a function of the OOO-angle, for the opening of the angle between two fixed O–O separations with fixed monomer geometries and fixed orientations corresponding to the trimer global minimum.

3-body interaction in TTM3-F (and all other such potentials) describes the long-range induced dipole–induced dipole interaction, and thus the short-range 3-body interaction, which in this plot occurs as the OOO angle decreases, is not expected to be accurately described. The comparisons of the 3-body fits to MP2/aVTZ energies as a function of the free and H-bonded OH-stretches are given in Figs. 6 and 7, respectively. As seen, the agreement is good and, as expected, the higher-order fit shows better agreement with the ab initio energies. The results from TTM3-F, also shown, underestimate the magnitude and variation of the 3-body energy. To conclude this subsection, we comment on the contribution of the 2- and 3-body energies for the trimer. For this cluster there are three 2-body contributions and, of course, just one 3-body contribution to the energy. At the global minimum, each 2-body energy is roughly −1800 cm−1 for a total 2-body energy of roughly −4800 cm−1 compared to the −300

(4)

−400

where s, as given in Eq. (3), is a function of the maximum OO-distance in a trimer configuration. The cutoff region defined by parameters a and b, is chosen to be between 7 and 8 Å, hence the 3-body energy of any trimer configuration with maximum OO-distance greater than 8 Å will not be evaluated. By doing so, we expect that the number of non-negligible 3body interactions will be greatly reduced for any large water cluster. Next, we show a comparison of the 3-body fits to MP2/aVTZ energies as a function of the OOO angle in Fig. 5, where again very good agreement is seen. The result from TTM3-F is in good agreement with the ab initio energies for angles greater than 120 but for angles less than 120 TTM3-F underestimates the magnitude of the 3-body energy by several hundred wavenumbers. This is not too surprising because the

−500

−1

Energy (cm )

free O−H stretch

−600 −700

MP2/aVTZ 6th−order 3−body fit 5th−order 3−body fit TTM3−F

−800 −900

0.7

0.9

1.1 RO−H (Å)

1.3

1.5

FIG. 6. One-dimensional cuts of indicated 3-body potentials as a function of the OH-distance of a H-bonded monomer O–H stretch with remaining degrees of freedom fixed at a trimer geometry corresponding to the trimer global minimum.

094509-5

Water potential and dipole moment surfaces 500

J. Chem. Phys. 134, 094509 (2011)

moment vector of the water dimer is represented a vector sum of effective (not physically meaningful) charges on six atoms,

MP2/aVTZ 3−body fit 5 −order 3−body fit TTM3−F th 6 −order th

Energy (cm−1)

−500

μ  dimer (X) =

qi (X) ri ,

(6)

i=1

−1500

−2500

−3500

6 

H−bonded O−H stretch

0.7

0.9

1.1 1.3 RO−H (Å)

1.5

1.7

FIG. 7. One-dimensional cuts of indicated 3-body potentials as a function of the OH-distance of a H-bonded monomer O–H stretch with remaining degrees of freedom fixed at a trimer geometry corresponding to the trimer global minimum.

single 3-body energy of roughly −700 cm−1 , about a factor of 7 smaller than the total 2-body contribution. As discussed in Sec. III, the contribution of the 3-body potential grows for larger clusters, as expected, because the number of the 3-body interactions grows faster than the number the 2-body interactions, i.e., roughly as n 3 vs n 2 , respectively, for n monomers. C. 4- and higher-body interactions

Finally, to complete the water potential for an arbitrary number of monomers, we add 4- and higher-body long-range polarization energies in addition to 1-, 2-, 3-body interactions described above. Again, we borrow from TTM3-F for these higher-body interactions. This is easily implemented by using the following expression:   Vn -body = Uind − Uind (i, j, k) n>4

+ (M − 3)



Uind (i, j),

(5)

where M is the number of water monomers and Uind represents the polarization induced interaction and is also the only nonpairwise term in the TTM3-F potential. Basically, this expression simply excludes from the full TTM3-F polarization energy all 2- and 3-body contributions, and this is done by subtracting polarization energies for all trimers from the polarization energy for the full water cluster and adding back discounted polarization energies for all dimers by a factor of M − 3. D. Dipole moment

The dipole moment surface (DMS) for water is constructed as a sum of 1- and 2-body dipole moment surfaces, both obtained from the full-dimensional DMS for the water dimer that has been reported previously.13 Briefly, the dimer DMS is a least-squares fit to MP2/aVTZ dipole moments of configurations in the same data base used for the intrinsic 2body fit; details of that fit are given in Ref. 13. The dipole

where X = ( r1 , r2 , r3 , r4 , r5 , r6 ) denotes the Cartesian coordinates of the dimer and ri is the xyz-coordinate for each atom i. Each of these effective charges qi (X) is an invariant polynomial representation of all internuclear distances. The fit expression of the dipole moment vector must be covariant with respect to all permutations of like atoms; that is it must transform such as a vector under permtutations of like atoms. We refer readers to the Appendix for more details on the permutationally covariant dipole fitting. Since the dimer dipole moment surface dissociates correctly to two monomers dipoles, the DMS of the water monomer (μ  1-body ) is naturally included in the dimer DMS in the dissociation limit. The expression for the dipole moment of a water cluster is denoted by DMS(1,2), since it is based on 1- and 2-body dipole moments and it is given by   DMS(1, 2) = μ  1-body (k) + μ  2-body (k, l). (7) k

k