Dissociation energy of the water dimer from Quantum Monte Carlo calculations

Dissociation energy of the water dimer from Quantum Monte Carlo calculations

Abstract

We report a study of the electronic dissociation energy of the water dimer using quantum Monte Carlo (QMC) techniques. We have performed variational quantum Monte Carlo (VMC) and diffusion quantum Monte Carlo (DMC) calculations of the electronic ground state of the water monomer and dimer using all-electron and pseudopotential approaches. We have used Slater-Jastrow trial wave functions with B3LYP-like single-particle orbitals, into which we have incorporated backflow correlations. When backflow correlations are introduced, the total energy of the water monomer decreases by about 4-5 mHa, yielding a DMC energy of 76.42830(5) Ha, which is only 10 mHa above the experimental value. In our pseudopotential DMC calculations, we have compared the total energies of the water monomer and dimer obtained using the locality approximation with those from the variational scheme recently proposed by Casula [Phys. Rev. B 74, 161102(R) (2006)]. The time step errors in the Casula scheme are larger and the extrapolation of the energy to zero time step always lies above the result obtained with the locality approximation. However, the errors cancel when energy differences are taken, yielding electronic dissociation energies within error bars of each other. The dissociation energies obtained in our various all-electron and pseudopotential calculations range between 5.03(7) and 5.47(9) kcal/mol and are in good agreement with experiment. Our calculations give monomer dipole moments which range between 1.897(2) and 1.909(4) Debye and dimer dipole moments which range between 2.628(6) and 2.672(5) Debye.

I Introduction

Water, as the main agent of all aqueous phenomena and an important component of the vast majority of all chemical and biological processes, has been the subject of many experimental and theoretical studies. The characteristic physical and chemical properties of water stem from its strong polar hydrogen bonds. In spite of the apparent simplicity of hydrogen bonding – a hydrogen atom bonded to an electronegative group and a lone electron pair on another system – understanding hydrogen bonding has proved very difficult.

One would hope to gain some insight by first investigating the water monomer and then progressing to the water dimer and larger water clusters, systematically including the many-body behaviour as the number of molecules increases. The study of the structure and energetics of assemblies of water molecules poses severe challenges to computational methods because they are held together by hydrogen bonds with binding energies of only a few kcal/mol.

The hydrogen bond in the water dimer has been the subject of many electronic structure studies Feller (1992); Mas and Szalewicz (1996); Famulari et al. (1998); Halkier et al. (1997); Klopper et al. (2000); Park et al. (2001); Xantheas et al. (2001); Tschumper et al. (2002) since it represents the prototype of all hydrogen-bonded systems. In recent years the water dimer has been studied using high-level quantum chemistry techniques, such as second-order Møller–Plesset perturbation theory (MP2) Xantheas et al. (2001) and coupled-cluster CCSD(T) methods Klopper et al. (2000). These methods can treat electron correlation effects quite accurately, but they are very expensive for larger water clusters, as the costs of MP2 and CCSD(T) calculations scale as and , respectively, where is the number of electrons. The accuracy of the results depends significantly on the quality of the basis set, and corrections are normally applied for the effects of basis incompleteness. Density-functional theory (DFT) shows a more favourable scaling with system size, allowing the use of larger basis sets resulting in smaller basis set errors. However, the calculated energies of the water monomer and dimer depend significantly on the exchange-correlation functional used.

Quantum Monte Carlo (QMC) methods Foulkes et al. (2001) represent an alternative and promising way of treating electron correlation. The cost of calculating energies with QMC scales roughly as , allowing accurate calculations for large systems. Moreover, QMC algorithms are intrinsically parallel and Monte Carlo codes are easily adapted to parallel computers. As the availability of powerful computers has increased, QMC has become a very attractive and effective method for probing the electronic structure of molecules and solids Drummond et al. (2005); Maezono et al. (2007).

In this paper we report variational and diffusion Monte Carlo (VMC and DMC) calculations of the water monomer and dimer, using all-electron (AE) and pseudopotential (PP) approaches. We have used Slater-Jastrow (SJ) and Slater-Jastrow-Backflow (SJB) wave functions, finding that the introduction of backflow correlations retrieves a substantial additional fraction of the correlation energy. Our calculations yield dissociation energies in very good agreement with experiment. We believe that our calculations form the most accurate QMC study of the water monomer and dimer performed to date.

The rest of the paper is organized as follows. In Section II we give short reviews of the VMC and DMC methods, the SJ and SJB wave functions, and the methods used for dealing with non-local PPs within DMC. In Section III we discuss the single-particle orbitals used and some details of the calculations. In Section IV we present and discuss our results. The first part of Section IV describes results for the water monomer, the second part gives results for the water dimer and its dissociation energy into two water molecules, and we compare results from the two different schemes used to calculate the non-local PP energy. In the third part, we report our results for the dipole moment of the monomer and dimer. Our conclusions are summarized in Section V.

Ii VMC and DMC Methods

In VMC the ground-state energy is estimated as the expectation value of the Hamiltonian with an approximate trial wave function, the integrals being evaluated by importance-sampled Monte Carlo integration. The trial wave functions contain variable parameters, whose values are obtained from an optimisation procedure formulated within VMC. There are no restrictions on the form of trial wave functions that can be used, and VMC does not suffer from a fermion sign problem. However, the choice of the approximate trial wave function is very important as it directly determines the accuracy and statistical efficiency of the calculation. Due to the difficulty of preparing trial wave functions of equivalent accuracy for different systems, the bias in the energy difference between systems is normally significant. We have used VMC methods mainly to optimize parameters in the trial wave functions, and our most accurate calculations have been performed with the DMC method.

In DMC the ground-state component of a trial wave function is projected out by evolving an ensemble of electronic configurations using the imaginary-time Schrödinger equation. The fermionic symmetry is maintained by the fixed-node (FN) approximation Anderson (1976), in which the nodal surface of the DMC wave function is constrained to equal that of the trial wave function. The FN DMC energy is higher than the exact ground-state energy, and becomes equal to it when the fixed nodal surface is exact. The dependence of the DMC energy on the trial wave function is smaller than in VMC, but in practice it is often significant. It is therefore desirable to be able to improve the nodal surfaces in order to reduce the impact of the FN approximation.

The standard SJ wave function can be written as

(1)

where R is a 3-dimensional vector denoting the position of each electron. The nodes of are defined by the Slater part of the wave function, , which takes the form , where is a Slater determinant of single particle orbitals of spin .

The Jastrow correlation factor, , contains electron-electron, electron-nucleus and electron-electron-nucleus terms, as described in Ref. Drummond et al., 2004. The Jastrow factor keeps electrons away from each other and greatly improves wave functions in general, but it does not modify the nodal surface of the wave function.

One way of reducing the FN error is to alter the nodes of the wave function by introducing backflow correlations López-Ríos et al. (2006), replacing the coordinates R in the Slater part of the wave function by the collective coordinates X, so that the SJB wave function reads

(2)

The new coordinates for each electron are given by

(3)

being the backflow displacement of particle , which depends on the position of every electron in the system. Details of the specific form of the backflow function used can be found in Ref. López-Ríos et al., 2006.

ii.1 Pseudopotential DMC calculations

We use non-local (angular-momentum dependent) PPs. Unfortunately non-local potentials cannot be used directly in DMC calculations because of the local character of the algorithm. Indeed, as we have already mentioned, the starting point of the DMC algorithm is the imaginary time Schrödinger equation. If the FN Hamiltonian, , contains a PP with a non-local part, , and is a trial energy, then the propagator for imaginary time is . This propagator contains terms of the form , which are not guaranteed to be positive for arbitrary and , and therefore cannot be interpreted as a probability density.

One way of avoiding the possible change of sign is to use the PP locality approximation (PLA) Mitas̆ et al. (1991). In this approximation the FN Hamiltonian is replaced by an effective Hamiltonian in which the non-local part operates on the trial wave function . This technique is analogous to the FN approximation. In this way the non-local part of the Hamiltonian is replaced by a local term and the PLA effective Hamiltonian reads:

(4)

where is the kinetic energy operator and is the local part of the PP. In this approximation, all of the matrix elements of the localized potential enter the branching factor of the DMC algorithm. A disadvantage of the PLA is that the energy may be higher or lower than the exact value.

Casula has recently introduced a different scheme Casula (2006) for treating non-local PPs within DMC. In the Casula scheme (CS), the FN Hamiltonian is substituted by an effective Hamiltonian

(5)

where

(6)

and

(7)

From Eqs. (5) and (6) we see that in the CS only the positive matrix elements of the non-local potential are localized, and they are absorbed into the branching factor. On the other hand, unlike in the PLA, the negative matrix elements of contribute to the drift-diffusion of the walkers. In practice, the standard DMC algorithm needs few changes: once a drift-diffusion move is proposed and accepted or rejected, and after the walker has been weighted, an extra displacement of the walker is performed according to a transition probability that depends only on the negative matrix elements of the non-local part of the PP. This method has two advantages over the PLA. First, the ground state energy of the CS Hamiltonian is an upper bound ten Haaf et al. (1995) on the ground state energy of the true FN Hamiltonian. Second, whenever the negative elements are large, the new displacement pushes the walkers away from the attractive regions of the localized potential, which avoids “population explosions” in which walkers are trapped in the attractive regions with large branching factors. We found occasional population explosions in our PLA calculations, which we dealt with by restarting the calculations at an earlier move with a new random seed. We found no such instabilities when using the CS.

Iii Details of the calculations

All the calculations have been performed using the casino code Needs et al. (2006). The single-particle orbitals forming the Slater determinants have been obtained from DFT calculations using the crystal98 code Saunders et al. (1998) with the Roos augmented double zeta ANO (Roos aug-DZ-ANO) Gaussian basis set Widmark et al. (1990) used in its decontracted form. Extensive tests of the basis sets listed at http://www.emsl.pnl.gov/forms/basisform.html showed that the Roos basis set gives very good DMC energies. The dependence of the DMC energies on the orbitals is generally quite weak. However, previous calculations for neon, neon radicals Gurtubay et al. (2006) and water Benedek et al. (2006) have shown that B3LYP orbitals give slightly lower DMC energies than Hartree-Fock (HF) orbitals, suggesting that the former can improve the nodal surface of the wave function. Based on this, we have investigated whether one can obtain further improvements in the DMC energies by optimising the form of the exchange-correlation (XC) functional used to generate the single-particle orbitals. Becke’s three parameter B3LYP XC functional Becke (1993) may be written as Saunders et al. (1998):

(8)

where determines the amount of Fock exchange, and and determine the amounts of non-local exchange and correlation, respectively. In Eq. (8), the Vosko-Wilk-Nusair correlation potential (VWN) is used to extract the local part of the LYP correlation potential. The usual B3LYP functional corresponds to , and . In principle we could try and minimize the DMC energy with respect to , and , but this would be very expensive. Setting in Eq. (8), we have obtained a DMC energy for the water molecule of 76.4218(1) Ha 1, which is very close to the value of 76.42205(8) Ha obtained with HF orbitals. This suggested that the contribution from the correlation part (terms involving ) does not significantly affect the DMC energy. We therefore decided to search for the amount of Fock exchange () which minimizes the DMC energy, while keeping and equal to their standard B3LYP values. The results are shown in Fig. 1. Note that the DMC energy (squares) is much less sensitive to the changes in the XC functional than the DFT energy (circles). The dashed line represents the best quadratic fit to the DMC energies for 0.6. The lowest DMC energy is obtained when 25% Fock exchange is used in the above functional. We therefore used to generate the single-particle orbitals for the Slater parts of all of the trial wave functions used in this work. In what follows, we refer to these as B3LYP-like orbitals.

Figure 1: Gurtubay et al., Journal of Chemical Physics

In our AE calculations, the single-particle molecular orbitals close to the nuclei have been corrected using the scheme described in Ref. Ma et al., 2005 so that each orbital obeys the appropriate cusp condition Kato (1957) at each nucleus, ensuring that the local energy is finite at each nucleus. In our PP calculations, we have used the two schemes described in Section II.1 with PPs generated within Hartree-Fock theory, including scalar relativistic effects Trail and Needs (2005a, b). These PPs have been found to work very well in conjuction with QMC methods Gurtubay et al. (2006). One of the advantages of using PPs is that they avoid the short-range variations in the wave function near the nuclei, and hence we can use larger time steps than is possible in AE calculations. We have performed calculations for different time steps and then extrapolated to the limit of zero time step by fitting to quadratic functions. Population-control biases have been checked by running with different population sizes, and no such bias was detected.

The parameters in the Jastrow and backflow functions were obtained by minimizing the variance of the local energy Kent et al. (1999); Drummond and Needs (2005), except for the AE calculations with SJB wave functions. In this case, the Jastrow and backflow parameters were obtained by minimizing the mean absolute deviation of the set of the local energies from the median.

Iv Results and discussion

The calculations for the water monomer were carried out in the experimental equilibrium geometry with r = r = 0.9572  and . This geometry was used in previously reported QMC calculations Lüchow et al. (1997); Benedek et al. (2006).

For the water dimer we have chosen the geometry derived by Klopper et al. Klopper et al. (2000) from CCSD(T) calculations, extrapolating to the basis set limit using MP2-R12 results. The optimized water dimer geometry (see Fig. 2) yields a final equilibrium oxygen-oxygen distance of 2.912 . When the dimer is formed, some structural deformation occurs in each of the monomers: a slight change in the angle between the inter-oxygen vector and the plane of the proton-accepting monomer occurs and the bond angle for the proton acceptor widens by about 0.5. We have neglected these small changes, assuming that their effect on the intermolecular binding energy is small. The electronic dissociation energy in the equilibrium geometry, , with respect to dissociation into two isolated monomers, is defined as

(9)

where is the experimental dissociation energy and ZPE denotes the zero point energy of the nuclei.

Figure 2: Gurtubay et al, Journal of Chemical Physics

iv.1 Water monomer

Our AE SJ and SJB QMC results with B3LYP-like orbitals for the water monomer are given in Table 1, together with energies calculated using other ab initio methods. The expansion order of the polynomials and the spin dependencies used in the Jastrow factor as described in Ref. Drummond et al., 2004 were the following: N=N=6, N=3, S=1 and S=S=0. This led to 128 parameters in the Jastrow factor. When backflow correlations were incorporated the total number of parameters to be optimized increased to 338, corresponding to the following expansion order of the polynomials in the backflow term (Ref. López-Ríos et al., 2006): N=N=6, N=3, S=S=1 and S=0. We have mentioned in Section II that the VMC energy depends strongly on the trial wave function. This is illustrated by Table 1: when backflow correlations are included, the VMC energy for the water molecule [76.4034(2) Ha] decreases by nearly 26 mHa. Our DMC energy extrapolated to zero time step with the SJ wave function (DMC-SJ), 76.42376(5) Ha, is nearly 1 mHa lower than that recently reported by Benedek et al. Benedek et al. (2006) using true B3LYP orbitals and lies about 14 mHa above the experimental result of 76.438 Ha. When backflow correlations are introduced, the DMC energy decreases by a further 4.5 mHa, which is 25% of the correlation energy missing with the SJ wave function, yielding a DMC energy of 76.42830(5) Ha. Note that, as mentioned in Section II, the difference between the energies obtained with the SJ and SJB trial wave functions is smaller at the DMC level than at the VMC level.

Lüchow and Fink Lüchow and Fink (2000) performed DMC calculations using a pair natural orbital configuration-interaction (PNO-CI) trial wave function consisting of 300 determinants and a Jastrow factor, obtaining an energy of 76.429(1) Ha, which is within error bars of our single-determinant SJB DMC result. For completeness, Table 1 also gives the DMC energy obtained by Casula et al. Casula et al. (2004) using a trial wave function consisting of a product of an antisymmetrized geminal power (AGP) wave function and a Jastrow factor, which is nearly 7 mHa above our best DMC result. Table 1 also gives the DMC energies calculated with SJ wave functions and HF nodes from Refs. Gurtubay et al., 2006 and Benedek et al., 2006. As was pointed out in Section III, both of these energies are slightly higher than the DMC results obtained using B3LYP and B3LYP-like orbitals, which suggests that the hybrid DFT functional yields better nodes than the HF orbitals. The small difference between the two DMC-HF values given in Table 1 arises from the use of different basis sets. For comparison, we have also included energies from a CCSD(T)-R12 calculation Müller et al. (1997) and a recently reported Correlation Energy Extrapolation by Intrinsic Scaling (CEEIS) calculation Bytautas and Ruedenberg (2006) which, although not variational, give results very close to the “exact” one.

Method Total energy (Ha) Basis set Ref.
DFT Roos Aug-DZ-ANO This work
VMC-SJ Roos Aug-DZ-ANO This work
VMC-SJB Roos Aug-DZ-ANO This work
DMC-SJ Roos Aug-DZ-ANO This work
DMC-SJB Roos Aug-DZ-ANO This work
DMC-HF 6-311++G(2d,2p) Gurtubay et al.,2006
DMC-HF Roos Aug-DZ-ANO Benedek et al.,2006
DMC-B3LYP Roos Aug-DZ-ANO Benedek et al.,2006
DMC-PNO-CI STO Lüchow and Fink,2000
DMC-AGP STO-DZ Casula et al.,2004
CCSD(T)-R12 Dunning cc-pV6Z Müller et al.,1997
CEEIS CBS limit Bytautas and Ruedenberg,2006
Exact Feller et al.,1987
Table 1: Total energy of the water monomer calculated with different ab initio methods, and the exact non-relativistic static point nuclei value. CBS limit stands for extrapolation to the complete basis set limit.

iv.2 Water dimer and electronic dissociation energy

As illustrated in Fig. 2, the two hydrogen atoms of the acceptor monomer (H) are not equivalent to the donor hydrogen (H) and free hydrogen (H) atoms in the donor subunit, and the two oxygen atoms (O and O) are inequivalent. However, we have obtained slightly lower DMC energies by using Jastrow and backflow functions which treat all four hydrogen atoms and both oxygen atoms as equivalent. Treating atoms as equivalent in this way reduces the number of variable parameters, and presumably this results in the better performance of the stochastic optimization procedure. When the H and O atoms are treated as equivalent, the SJ trial wave functions for the dimer contains 128 parameters (for both the AE and PP calculations) while the SJB calculations have been performed with 390 parameters for the AE calculations and 286 parameters for the PP calculations. These number of parameters correspond to N=N=6, N=3, S=S=1 and S=0 in both the AE and PP Jastrow factors. In the backflow term, N=N=6, S=S=1, S=0 in both AE and PP calculations, and N=3 and N=2 in the AE and PP cases, respectively.

Figure 3: Gurtubay et al, Journal of Chemical Physics

Figure 3 shows the AE DMC energies of two water monomers and the water dimer as a function of time step. The zero-time-step extrapolations of the fitted curves are shown by filled triangles of the same color as the original symbols. The DMC-SJ and DMC-SJB values shown in Table 1 correspond to the black and blue triangles in this figure, i.e., the extrapolated values from the curves labelled 2HO(SJ) and 2HO(SJB), respectively. As we mentioned previously for the water monomer, the DMC energy decreases by about 4.5 mHa when backflow correlations are introduced. Figure 3 shows that the reductions in energy from including backflow correlations are essentially independent of the time step for both the monomer and dimer.

Figure 4 shows PP SJ DMC energies for two water monomers and the dimer as a function of the time step. The CS and PLA curves cross at a time step of about au. The extrapolated SJ DMC energy of two water monomers obtained within the CS scheme is less than 4 mHa above the PLA value. It is not surprising that the CS energies are above the PLA ones, because the CS scheme gives an upper bound on the true energy, while the PLA scheme does not. The energy of the water monomer as a function of the time step in both the PLA and CS schemes can be found in Table 2.

Figure 4: Gurtubay et al, Journal of Chemical Physics
Energy (Ha)
Time step (au) PLA CS
Table 2: Comparison of PP SJ DMC energies obtained within the PLA and CS schemes.
Figure 5: Gurtubay et al, Journal of Chemical Physics

In Fig. 5 we have summarized the total energies of the PP water monomer and dimer energies with SJ and SJB wave functions within the PLA. The CS yields a similar figure (not shown) with curves parallel to the red squares of Fig. 4. As in Fig. 3, including backflow correlation in the monomer and dimer results in an essentially rigid shift of the corresponding SJ curve to lower energies. Furthermore, because the black and red curve corresponding to the SJ calculations, and the blue and green curves, corresponding to the SJB calculations, have nearly the same shape, the time step errors largely cancel when computing the of the water dimer. This cancellation is illustrated in Fig. 6, where we have plotted as a function of the time step for the four combinations of SJ/SJB wave functions within the AE/PP approaches from Figs. 3 and 5.

Figure 6: Gurtubay et al, Journal of Chemical Physics

The top panel of Fig. 6 shows the AE values of . The magenta diamonds (orange triangles) show the results with SJ (SJB) wave functions as a function of time step. We fitted the energies to quadratic functions and extrapolated to zero time step, obtaining the values of given in Table 3.

Method (kcal/mol) Ref.
AE-SJ This work
AE-SJB This work
PP-SJ (PLA) This work
PP-SJ (CS) This work
PP-SJB (PLA) This work
PP-SJB (CS) This work
AE-HF-SJ Benedek et al.,2006
AE-B3LYP-SJ Benedek et al.,2006
CCSD(T) Klopper et al.,2000
Experiment + Theory Curtiss et al.,1979
Experiment + Theory Mas et al.,2000
Experiment + Theory Leforestier et al.,2002
Table 3: Top: DMC dissociation energies (kcal/mol) obtained in this work with SJ and SJB wave functions using B3LYP-like orbitals within the AE and PP (PLA and CS) approaches. Bottom: Results from other ab initio calculations and from “Experiment + Theory”. No error bar was reported in Ref. Leforestier et al., 2002.

The high degree of cancellation of the time step errors shown in Fig. 6 is remarkable. Our AE-SJ data show smaller time step errors than those reported in Ref. Benedek et al., 2006 with HF single-particle orbitals. When backflow correlations are introduced, the cancellation of time step errors at larger time steps is not as good. This can also be observed in Fig. 3 where the green line bends more than the blue one at larger time steps. Nevertheless, the extrapolated AE values of for the SJ wave function [5.16(10) kcal/mol] and the SJB wave function [5.35(10) kcal/mol] are within error bars of each other. For the PP calculations we have seen that the PLA and CS schemes give significantly different time step errors. However, the time step errors are very similar for the monomer and dimer (see Fig. 4). Therefore the values of computed within the PLA and CS schemes are very similar at each time step. We find 5.03(7) kcal/mol [5.07(7) kcal/mol] for SJ wave functions within the PLA [CS], and 5.47(9) kcal/mol [5.38(8) kcal/mol] when backflow correlations are introduced. We have also found that, in both our AE and PP calculations, backflow correlations tend to increase the value of slightly.

Our DMC results range between 5.03(7) and 5.47(9) kcal/mol and are similar to the DMC values of reported in Ref. Benedek et al., 2006 (with HF and true B3LYP orbitals) and the CCSD(T) result of Klopper et al. Klopper et al. (2000) of 5.02(5) kcal/mol.

For many years the experimental value of has been taken to be 5.440.7 kcal/mol. This value was derived from measurements of the enthalpy of dimerization of water corrected with a theoretical estimate of the ZPE Curtiss et al. (1979). Mas et al. Mas et al. (2000) have recently calculated a more accurate value of the ZPE and have estimated to be 5.000.7 kcal/mol. Leforestier et al. Leforestier et al. (2002) have reported a value of of 5.14 kcal/mol, correcting for the ZPE by determining potential energy surfaces for the monomer and dimer via direct inversion of spectroscopic data.

iv.3 QMC dipole moment of the water monomer and dimer

We have calculated the mean dipole moment, , of the water monomer and dimer. The dipole moment operator does not commute with the Hamiltonian and therefore the error in the DMC mixed estimate () is linear in the error in the wave function. This error can be reduced to a quadratic error by using the so called extrapolated estimator, , where is the VMC estimate of the dipole moment. Both VMC and DMC estimates must, of course, be calculated using the same trial wave function.

Figure 7: Gurtubay et al, Journal of Chemical Physics

Figure 7 shows our extrapolated estimator results for the dipole moment of a water molecule as a function of time step. Red squares (black triangles) show AE SJ (SJB) results while green circles (blue diamonds) represent PP SJ (SJB) results calculated with the CS. Dashed lines of the same color as the original symbols are linear fits used to find the zero-time-step dipole moment. Our results are summarized in Table 4 together with the results of other ab initio calculations. Our results for the monomer are slightly above (0.05 Debye) the experimental value, as are the other theoretical values. For comparison, the inset of Fig. 7 shows the DMC mixed estimates of the dipole moment for each of the four cases investigated. As we have mentioned above, the DMC estimate has a larger dependency on the trial wave function and this can be observed in the larger differences between the data sets displayed in the inset. However, when the extrapolated estimator is considered the error is reduced and the results come into closer agreement.

Method (Debye) Ref.
AE-SJ This work
AE-SJB This work
PP-SJ (CS) This work
PP-SJB (CS) This work
aug-cc-VDZ MP2 Gregory et al.,1997
aug-cc-pVQZ MP2 and MP4 Xantheas and Dunning,1993
aug-cc-pVDZ MP2-R12 Klopper et al.,1995
Experiment Gregory et al.,1997
Table 4: Top: Dipole moment of the water monomer, evaluated using the extrapolated estimator, and also extrapolated to zero time step. Bottom: Other ab initio results and the experimental value.

We have arrived at similar conclusions for the dipole moment of the water dimer (figure not shown). The extrapolated estimator values of the dipole moment of the water dimer, extrapolated to zero time step, are summarized in Table 5. Our AE results are lower than the experimental value Dyke et al. (1977), but within 0.02 Debye of it, and our PP results are higher than the experimental value, but within 0.03 Debye of it.

Method (Debye) Ref.
AE-SJ This work
AE-SJB This work
PP-SJ (CS) This work
PP-SJB (CS) This work
aug-cc-VDZ MP2 Gregory et al.,1997
aug-cc-pVDZ B3LYP Xantheas,1995
Experiment Dyke et al.,1977
Table 5: As in Table 4, but for the water dimer.

V Conclusions

We have performed all-electron and pseudopotential QMC calculations of the total energy of the water monomer and dimer using SJ and SJB wave functions with B3LYP-like single-particle orbitals. The introduction of backflow correlations recovers a significant amount of correlation energy and we have obtained a final DMC energy (extrapolated to zero time step) for the water molecule of 76.42830(5) Ha. We have performed a careful study of the time step errors and have shown that they largely cancel when energy differences, such as the electronic dissociation energy of the water dimer, are considered. We have compared two different schemes for treating non-local PPs within DMC: the pseudopotential localization approximation and a scheme recently proposed by Casula. Although the total energies from the two schemes differ significantly, these differences are essentially the same in the monomer and dimer, and therefore they cancel when evaluating . Our extrapolated AE and PP results for the show good agreement with each other when both SJ and SJB wave functions are used, and they are all within the error bars of the values determined from experiment. The total spread of our AE and PP values of is only 0.44 kcal/mol. We have also reported QMC estimates of the dipole moments of the water monomer and dimer. Our AE and PP extrapolated DMC results yield a dipole moment of the monomer within 0.05 Debye of the experimental value, while for the dimer the AE results are within 0.02 Debye of experiment and the PP results are within 0.03 Debye of experiment.

Acknowledgements.
I.G.G. gratefully acknowledges financial support from the Basque Government through grant BFI04/183, the EC 6 framework Network of Excellence NANOQUANTA (NMP4-CT-2004-500198), and the Ministerio de Educación y Ciencia (FIS2006-01343 and CSD2006-53). R.J.N. gratefully acknowledges support from the Engineering and Physical Sciences Research Council (EPSRC) of the United Kingdom. Computational facilities were provided by the Cambridge High Performance Computing Service and the SGI/IZO-SGIker UPV/EHU (supported by the National Program for the Promotion of Human Resources within the National Plan of Scientific Research, Development and Innovation-Fondo Social Europeo, MCyT and the Basque Government).

Figure captions

Fig. 1: DMC energy (squares) as a function of the percentage of Fock exchange in the exchange-correlation B3LYP-like functional [ in Eq. (8)] used to obtain the single-particle orbitals for the SJ wave functions. The same Jastrow factor of 128 parameters was used in all the calculations. The dashed line shows a quadratic fit to the points around the minimum. The dotted line represents the DMC energy with HF orbitals. Inset: DFT energy as a function of the percentage of Fock exchange.

Fig. 2: (Color online) The equilibrium structure of the water dimer from Klopper et al. Klopper et al. (2000) The “a” stands for acceptor, “d” for donor and “f” for free.

Fig. 3: (Color online) All-electron total energies of two water monomers and the water dimer as a function of time step, calculated with SJ (circles) and SJB (squares) wave functions. Dashed and dotted lines show quadratic fits to the monomer and dimer energies, respectively. The filled triangles represent the extrapolation to zero time step. The statistical error bars are of order  Ha.

Fig. 4: (Color online) SJ PP energies for two water monomers and for the water dimer, calculated with the PLA (black circles) and the CS (red squares). Dashed (dotted) lines are fits for the monomer (dimer). Triangles denote the energies extrapolated to zero time step.

Fig. 5: (Color online) As in Fig. 3, but using PPs within the PLA.

Fig. 6: (Color online) Dissociation energy of the water dimer as a function of time step. The curves represent time step extrapolations and filled symbols give the zero time step . Top panel: AE calculations. Bottom panel: PP calculations. For comparison, the AE values extrapolated to zero time step are also shown (with the same colors and symbols as in the top panel).

Fig. 7: (Color online) Extrapolated estimates of the dipole moment of the water molecule as a function of the time step. Filled symbols correspond to zero time step values obtained from linear fits (dashed lines). Inset: DMC mixed estimator.

Footnotes

  1. The value in brackets denotes the standard error in the mean in the last figure. In this case, for example, 76.4218(1) means 76.42180.0001.

References

  1. D. Feller, J. Chem. Phys. 96, 6104 (1992).
  2. E. M. Mas and K. Szalewicz, J. Chem. Phys. 104, 7606 (1996).
  3. A. Famulari, M. Raimondi, M. Sironi, and E. Gianinetti, Chem. Phys. 232, 275 (1998).
  4. A. Halkier, H. Koch, P. Jørgensen, O. Christiansen, I. M. B. Nielsen, and T. Helgaker, Theor. Chem. Acc. 97, 150 (1997).
  5. W. Klopper, J. G. C. M. van Duijneveldt-van de Rijdt, and F. B. van Duijneveldt, Phys. Chem. Chem. Phys. 2, 2227 (2000).
  6. C.-Y. Park, Y. Kim, and Y. Kim, J. Chem. Phys. 115, 2926 (2001).
  7. S. S. Xantheas, C. J. Burnham, and R. J. Harrison, J. Chem. Phys. 116, 1493 (2001).
  8. G. S. Tschumper, M. L. Leininger, B. C. Hoffman, E. F. Valeev, H. F. Schaefer III, and M. Quack, J. Chem. Phys. 116, 690 (2002).
  9. W. M. C. Foulkes, L. Mitas̆, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  10. N. D. Drummond, A. J. Williamson, R. J. Needs, and G. Galli, Phys. Rev. Lett. 95, 096801 (2005).
  11. R. Maezono, A. Ma, M. D. Towler, and R. J. Needs, Phys. Rev. Lett. 98, 025701 (2007).
  12. J. B. Anderson, J. Chem. Phys. 65, 4121 (1976).
  13. N. D. Drummond, M. D. Towler, and R. J. Needs, Phys. Rev. B 70, 235119 (2004).
  14. P. López-Ríos, A. Ma, N. D. Drummond, M. D. Towler, and R. J. Needs, Phys. Rev. E 74, 066701 (2006).
  15. L. Mitas̆, E. L. Shirley, and D. M. Ceperley, J. Chem. Phys. 95, 3467 (1991).
  16. M. Casula, Phys. Rev. B 74, 161102 (2006).
  17. D. F. B. ten Haaf, H. J. M. van Bemmel, J. M. J. van Leeuwen, W. van Saarloos, and D. M. Ceperley, Phys. Rev. B 51, 13039 (1995).
  18. R. J. Needs, M. D. Towler, N. D. Drummond, and P. López-Ríos, casino version 2.0.0 User Manual, University of Cambridge, Cambridge (2006).
  19. V. R. Saunders, R. Dovesi, C. Roetti, M. Causà, N. M. Harrison, R. Orlando, and C. M. Zicovich-Wilson, crystal98 Users’s Manual, Università di Torino, Torino (1998).
  20. P. O. Widmark, P. A. Malmqvist, and B. Roos, Theor. Chim. Acta 77, 291 (1990).
  21. I. G. Gurtubay, N. D. Drummond, M. D. Towler, and R. J. Needs, J. Chem. Phys. 124, 024318 (2006).
  22. N. A. Benedek, I. K. Snook, M. D. Towler, and R. J. Needs, J. Chem. Phys. 125, 104302 (2006).
  23. A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
  24. A. Ma, M. D. Towler, N. D. Drummond, and R. J. Needs, J. Chem. Phys. 122, 224322 (2005).
  25. T. Kato, Commun. Pure Appl. Math. 10, 151 (1957).
  26. J. R. Trail and R. J. Needs, J. Chem. Phys. 122, 014112 (2005a).
  27. J. R. Trail and R. J. Needs, J. Chem. Phys. 122, 174109 (2005b).
  28. P. R. C. Kent, R. J. Needs, and G. Rajagopal, Phys. Rev. B 59, 12344 (1999).
  29. N. D. Drummond and R. J. Needs, Phys. Rev. B 72, 085124 (2005).
  30. A. Lüchow, J. B. Anderson, and D. Feller, J. Chem. Phys. 106, 7706 (1997).
  31. A. Lüchow and R. F. Fink, J. Chem. Phys. 113, 8457 (2000).
  32. M. Casula, C. Attaccalite, and S. Sorella, J. Chem. Phys. 121, 7110 (2004).
  33. H. Müller, W. Kutzelnigg, and J. Noga, Molecular Physics 92, 535 (1997).
  34. L. Bytautas and K. Ruedenberg, J. Chem. Phys. 124, 174304 (2006).
  35. D. Feller, C. M. Boyle, and E. R. Davidson, J. Chem. Phys. 86, 3424 (1987).
  36. L. A. Curtiss, D. J. Frurip, and M. Blander, J. Chem. Phys. 71, 2703 (1979).
  37. E. M. Mas, R. Bukowski, K. Szalewicz, G. C. Groenenboom, P. E. S. Wormer, and A. van der Avoird, J. Chem. Phys. 113, 6687 (2000).
  38. C. Leforestier, F. Gatti, R. S. Fellers, and R. J. Saykally, J. Chem. Phys. 117, 8710 (2002).
  39. T. R. Dyke, K. M. Mack, and J. S. Muenter, J. Chem. Phys. 66, 498 (1977).
  40. J. K. Gregory, D. C. Clary, K. Liu, M. G. Brown, and R. J. Saykally, Science 275, 814 (1997).
  41. S. S. Xantheas and T. H. Dunning, Jr., J. Chem. Phys. 99, 8774 (1993).
  42. W. Klopper, M. Schütz, H. P. Lüthi, and S. Leutwyler, J. Chem. Phys. 103, 1085 (1995).
  43. S. S. Xantheas, J. Chem. Phys. 102, 4505 (1995).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minumum 40 characters
Add comment
Cancel
Loading ...
108212
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description