Quantum Monte Carlo calculations on Be and LiH

Volume 113, ntimhzr 3

CHEMICAL PHYSICS Ll3TER.S

QUANTUM MONTE CARLO_CALcuLATION5

18 January 1985

ON Be AND LiH

RJ. HARRISON and NC. HANDY Received 25 September 1984;in final form 23 October 1984

Quantum Monte Carlo (QMC) methods, as recently developed for molecular systems, are applied to Be and LiH. The importance of the trial wavefunction, written as a product of a correlation factor and an orbital p-x& is emphansed. It is shown Sat dgnifieant ~provemen~ in the accuracy of the approach are aclueved ~rnulti~on~~ra~ion wavefunctions are used in preference to self-consistent tield wavefnnctdons. Vanous forms of the correlation factor are investigated.

1. Introduction It is generally recognised that the principal weakness

of ab tiitio calculations is the basrs set deficiency_ Fur~e~ore whilst it is possible to use extended basis sets such that the Hartree-Fock limit is reached in selfconsistent field (SCF) calculations, it is not possible to reach the’full. correlation energy limit because of the slowly converg%t nature of the wavefunction expansion. This is because higb+:der polar&&ion basis functions d, f, g... are required to obtdrn the true shape of the wavefunction in the region of the electroneIectron

cusp. ‘Ibis is one of the uuun reasons

there has been much interest

recently

why

in the Green’s

Monte Carlo (GFMC) methods which solve the Schrodinger equation in amore direct fashion. It is not necessary to introduce the theory of thrs approach which is thoroughly covered m papers by Reynolds, Ceperley, Alder and Lester El ] , Ceperley [2,3]. Ceperky and Kales [4] and Moskowitz. Schmidt, LRe and Kales [s]. We use here the &tort time algorithm [l J] (also tenned diffuson Monte Carlo by Ceperley [Z]) with the fmed node approximation [1,4,5], although we have aLso ~plemented the exact WMC algorithm of Ceperley 1233, again within the fixed node approximation. We focus our function

(1) Within the fixed node approximation the ultimate energy obtained corresponds to the best variational energy possible within the nodal surfaces of the trial wavefunction. (2) The variance associated with any energy estimate is directly related to the accuracy of WT: the better the t&d wavefunction the smaller the variance for equal length simulations. (3) If 3’~ is exact then the exact energy is obtained

zero variance whatever the timestep (7) in the short time (ST) approximation.

with

Thus the more accurate that %I$ is, the larger may be the time step, ieading to more efficient sampling. Calculations with small time steps are very expensive due to the increased correlation of successive energy estimates, longer runs being needed to produce reliable error estimates. This paper is concerned with the fom of the trial wave~nction and takes note of the advances made by ab initionists in recent years, especiaUy towards the goal of obWg reasonably compact wavefi.rnctions.

2. The form of the trial wavefunction

The trial wavefunction is written as a product of

attention upon the use of importance sampling [ 1,4]

correlation factor C and an orbital factor a, i.e.

through the ~trodu~tion of a trial wave~nction !I&_ This trial wave~nct~on has a huge effect on the calcu-

*~=CG.

Iations

The factor is chosen such that the correiation cusp

for several reasons_

0 009-2614/85/S ~o~~-Ho~and

03.30 0 &evier Physics oblong

!%zienceF’ubliahera B.V. Division)

a

0)

257

Volume 113,number 3

CHEMICALPHYSICSLEITERS

condition is obeyed exactly for rdl electron pairs. A form used by many workers is the Jastrow factor [6] J = ,p,

exp ( ufj)

*

(3

where in its simplest form Uij = Qv'fj/(l+ b'ij), 1 aij=u

(3)

for ij of like spin,

=:

for i,j of unhke spin _

The simplest form which may be used for the orbital part @ is a single Slater determinan t, wluch for closedshell s-stems is det( @f +$..-$z),

where III tlu

(4)

notation $, are spatial orbit&

expanded

in some basis set {x)

(5) Since we are not concerned with the spin variables the form of the orbrtal wavefunction that would actually be used is

det(D&)

det(D$)

,

(6)

where

The orbital part 9 is crucially important because it defimes the nodal surfaces. The dimension of the basis set (m) need not be too large to approach the SCF limit, double zeta perhaps being sufficient. In recent years much progress has been made beyond the SCF limit and in particular rt is possible to determine multi-configuration SCF (MC SCF) wavefunctrons. The dominant aspects of non-dynamical electron correlation can be included with a small number of configurations, and presumably, it is these which wdI determine. except in fine detail, the location of the nodal surfaces. &Iere therefore we shrill constder the effects of representing 0 as a short MC SCF expansion

18January1985

Although the calculation of aT and its denvatives accounts for most of the computer time in a QMC calculation, the expense of the calculation need not be greatly increased by using such a trial wavefunction. There will be many common orbit& and cofactors between the various dete rminants, the most time consuming part (for the four-electron systems studied here) being the evaluation of the orbitals themselves and their derivatives. In the calculation reported below the MC SCF calculations on Be took only twice as long as those using an SCF orbital function for the same number of time steps, even though four new basis functions were introduced. The QMC calculation on Be with an MC SCF orbital wavefunction using a time step of 0.0 1 h-l took approximately 6 h of IBM3081 D time. Once the form of + has been fmed, the effect of the correlation factor is only to improve the efficiency of the simulation. A better C will reduce the vanance; it cannot alter the nodal pattern ifits form is given by (2). We cor&d+r ‘he introduction of higher terms in ’ ‘Lherepresentation of Uij, similar to those included by Boys and Handy in their transcorrelated method [7] We &all also consider a different form for Cintroduced by CoIIe and SaIvetti [8,9] and shown to give good results in some non-variational calculations [8-lo], though a recent variational calculation [ 1 l] has found it to be very poor_

#(q,

exp(-_P2&) JI x [l - rr1’2B(1 + aij~~j)/(l + ~"2p)] , rj) =

P=4CP(ri)P(Q)11'6

-

Here aij as in (3); p(r) is the electron density from the orbital function 9. This is not exactly the form of Colle and Salvetti in the definition of 0, but it is simpler computationally. Only 4 has to be determined. It may also be considered necessary for e-r to obey the nuclear cusp condition, i.e. a term such as

(9) each determinant, $I~,being interpreted according to (6). 258

ought to be included, where riA =-I ri -A I and A is a nuclear coordinate. If a Slater-type basis is used in (5) then the nuclear cusp is quite well obeyed and in

Volume 113, number 3

CHEMICAL

some t&A ~vest~gations on Be no appreciable difference in the v~tio~ energy or its variance was found between calculations in which (9) was neglected, or if it was included with Z, treated as a parameter. The calcrrlations reported here therefore only employ an electron-electron c&relation factor

3. CaIcuiations We report calculations on He, Be2+, Be and L.iH. The short time algorithm of Reynolds et al. [l] is used, the only difference being that al,! of the electrons are moved at once reducing the number of times q= and its derivatives need be computed. The effect of this is

PHYSICS LETTERS

18 January 1985

slightly increase the number Of rejeCtionS wken detailed balance is imposed from approximately 1% to 5% for the time steps used and it does not appreciably affect the variances computed for a given number of time steps. To investigate the error arising from the use of a ftite time step in the short tune approximation we report c~~u~tions on He, Be3 and Be. Using a sixtwn Hylleraas 1121 trial wavefunction for He and a tnne step (T) of 0.01 h-* it is seen (table 1) that the mor in the best ST result is a factor of ten less than the error of the variational energy of eT and any systernatic error due to the use of finite ?-is certainly less than the statktrcal uncertainty in the enerw (50 & or 11 cm-l). The results in table 1 obtained with a only to

Table 1 Helium QMC results Method a)

Trial wavefunctron h)

Enerby e)

variational ST(O.L!l) varictional ST(O.01) ST(O.01)

two-term Hylleraas e) two-term Hylleraas e) s&term HylJeraas e) six-tcrrn Hylleraas e) six-term HyUeraase)

variational

SCF X J (O-5)

variational

SCF x C(1.4)

-2.8766 (35.5%) -2.904 f 0.004 (100 t 10%) -2.90324 (98.86%) -2.9039 r 0.0004 (100 _+ 1%) -2.90374 2 0.00005 (100.05 i 0.1%) -2.8836 & O-00065 (52i 15%) -2.8956 ‘I 0.0006 (80.7 * 1.5%) -2.90372

exact f)

Detisd)

(100,100,25) (200,200,12) (200,200.1156) (200~0,475f (200,50,474)

a) Variational, energy of trjal wavefunction:ST~~),ST algorithmwith *tiestep r b) l(b), cJectrorncJastrow factor wnh b as in (3); C(q). CoBe-Sdvetta factor wrtb optimal q as in (8). c) Energy with computed error and fraction of correlation energy recovered. d) flV, m, n)N = number of indnnduul configuratmns, m = number of steps per b1ock.n = number of bIocks over which overages were computed. e, Ref. 1121. f, Ref. [13].

Table 2 Be2+SCF X correlation factor QMC results

a-d.9

Method a)

Corretation factor b)

Energy c)

variational variational ST(O.01) variational exact f)

J(1.8) J(l.8) C(1.60)

-13.6113 -13.6390 -13.6549 -13.6425 -1365557

Details d, (0%) f ~.0007 (63 * 1.5%) f 0.00051(98~ * 1.2%) f 0.0006 (705 5 1.4%)

(800,100,180) (800,~00,148~ (800,10fl.191)

Footnotes as for table 1,

259

CHEMICALPHYSICSLETTERS

Volume113,number3 TabIc 3 Beatom SCFX

correIationfactorQMCresults

Method

a)

variational variational ST(O_OS) ST(OO.01) ST(O_OO5? variational exacte)

a-d) Footnotesas~~~table1.

Correlationfactorb)

Energy=)

fietdsd)

J(1.8) J(l.8) J(1.8) J(l 8) C(2.0)

-145724(0%) -14.6089* 0.0008 (38~ 0.8%) -_14.6580~ 0.0005 (PO* 0.5%) -14.6564 * 0.0006 (88 = 0.6%) -14.6567 XIZ o_ooo6 (88*0.6%) --14.6091* 0.0011(38* l-t%> -14.66733

(800,100,135~ (800,50,327)* (800~0,347~ (8GO,lOO,354) (80050.2703

e) R~F. [ls]_

less accurate two-term Hylleraas expansion [12] further emphasise the gains in efficiency that may be obtained by using accurate trial wavefunctions. Calculations on Be*+ (table 2) using !Z?= of the form (l), where ip is an SCF fun&Ion, shows that the ST(7 = 0.01) result lies 0.8 mh above the exact result. which is only just greater than the standard deviation u. It is of course possible to perform calculations for several values of 7 and extrapolate to 7 = 0, as I& been done by other workers [ 1.141_ Tables 3 and 4 contain ST results for Be obtained with two different XV= for several vales of 7. It does not seem possible to use this small amount of data to extrapolate reliably to the 7 = 0 result, but It is again clear that the systematic ST error is of the same magnitude as the statistical errors quoted when 7 = 0.01 h-l. The only calculation to be discussed below in wluch the ST error may be significant 1s the most accurate calculation in table 5 on LiH, which differs from the exact eigen-

Table4 Bcatom MCSCFf)

value by O-76 mh or three times the given error. Since at least part of tills error is attributable to the fmed node approximation, which if implemented exactly would be variational, the ST error must be less than this unless some cancellation of errors is occurrmg. The details of the Slater type bavs sets used in these calculations are given in table 6_ Comparing the results for Be in tables 3 and 4, it appears that the ftved node results are far superior for the MC SCF @ rather than the SFC ip. The MC SCF wavefunctlon was carefully .bosen to look after near sp degeneracy and to introduce radial correlation effects for both the 1s and 2s orbit&. The situation IS different for LiH where the results show that the nodal surfaces for SCF + must be very accurate, though of course this will not be true at displaced geometries. However in these calculations it is found that improving @ does not apprec&bly affect the variance obtained from a given length simulation. This is rather dlsap-

x correlation factor QMC results

Method")

CorreIatlon factorb)

variational VZ3lXitiOIld

ST(O.05) ST(O.01) ST(O.005) sr(0.001) variational exacte) a-e) Footnotesasfortable3. 260

18Januaq1985

J(l.8) J(1.8)

JU 8)

J(1.8) J(l_8) C(2.35)

Energyc)

Detailsd)

-14.6280 (5896) -14.6478 * 0.0007 (79* 0.8%) -14.6637iO.0007 (96 f 0.8%) -I4 6651= 0.0004 (97.52 0.4%) -14.66585 0.0005 (98-O* 0.4%) -14.6657* 0.0007 (98* 08%) -14.64132 0.0009(72* 1%) -14.66733

(8OO,lOO,l24> (800,50.176) (800.100.270~ (800,100.271) (800,200,278) (80050,333)

0 Q,=O.9519226llsz2sz1-O.3O5458~ls22~2l-O.O232124~2s'z~~..

_

CHEMICALPH~SICS

Volume 113. number 3 Table 5 LiH QMC results atR

18 January 1985

LEmRs

= 3.015 00

Me@od a)

Trial wavefunction b)

Energy C)

VtitiOlKd

SCF

-7.98650(0%)

variational

SCF x J(1.35) SCF X /(l-35) MC SCF(2) f) MC SCF(2) x J(1.35) MC SCF(2) x J (1.351 MC SCF(4) MC SCF(4) x I(1.35) MC SCF(4) X I(1.35)

-8.0264 -8.0696 -8.0027 -8.0375 -8 0701 -8.01364 -8 0428

ST(O.01) variational variational ST(O.01) variational variational ST(O.01) exam c)

Dettisd) + 0.0004 i 0.0007 (1%) f 0.0004 -c 0.0004 (32%) + 0.0005

(47 + 0.5%) (98.9 f 0.8%)

(1000.20,821) (1000,50,99)

(61 f 0.5%) (99.5 * 0.5%)

(1000,20,789) (1000.20.401)

(67 f 0.6%)

(100050.657)

-8.06973 f 0.00026 (99.15 O-354) (1000,50.887)

-8.OTO49

a-d) Footnotes

8s for table 1_ e) Ref. [ 16). f)MCSCF(2)=O.YY2[lc~~2cr~~-O.l259~lo~3~~,MCSCF(4)

pointing, and not fully understood, because (as observed earlier) it is not much more time consuming to evaluate an MC SCF @ than an SCF +. It will be the correlation factor which is most important in determining the efficiency of the simulation through its critical effect on the LGsmoothness” of the local energy hypersurface which controls the variance. As found previously [ 1,5,11,14] the simple Jastrow factor (2) recovers 20--50% of the missing correlation energy for these small systems in a variational calculaTable 6 Slater basis sets used for QMC orbital wavefunctions :ie 1s 2.91093,1.45363 Be*+ 1s 3.43071,5-6315.7.35143 Be SCF 1s 5.59108.3.35538 2s 1.01122.0-61 Be MC SCF 1s 559108.3.35538 2s 4.109,1.01122,0.61 2p 0.98 LIH Li 1s 4.699,2.52117 2s l-2.0.79722 2~ 2.75, 1.2,0.73691 H 1s 1.56567.0-88775 2s 2.2 2p 1.37646

=0.987711u22~21-0.127711cr2~21-0.0Y0111~lnZI.

tion, but we find that the energy and its variance (~2) are not strondy dependent upon the single vanational parameter, b in eq. (3). We attempted to unprove J by the inclusion for Be of a higher order term in U

(10) optimising y variationally, but this gave no reduction in u and only recovered a further l-3% of the correlation energy. In the tables results are also presented using the Colle-Salvetti correlation factor (8). The energy and variance were found to be strongly dependent upon the single parameter q which was easily optimised in variational calculations which formed correlated averages for several values of q simultaneously. For the twoelectron systems the factor is seen to recover a larger proportion of the correlation energy than the Jastrow factor without significantly increasing the efficiency of the calculation whereas for Be u is substantraliy increased with no improvement in the energy. The seemingly excelient results previously obtained with a function of tti type [S-lo] must be due to the non-variational nature of these calculations, whilst the poor variational result of Moskowitz [ 1 l] is probably attnbutable to using a non-optimum value of q-

4. The dipole moment of LiH The simulation, with importance sampling, yields 261

Volume 113,nurri*r 3

CHEMICAL PHYSICS LETTERS

a set oipoints which are asymptotically distributed proportional to ‘k,ipo, wlrere +o is the exact ground state wavefunction (here with the fared node constraint)_ Any expectation value of a property $ wiIl provide a mixed estimate WM = s*+,,

dR/J

*+Q,

dR.

(11)

Ceperley and Kalos [4] showed how a perturbation type extrapolated estimate of the desired expectation value could be obtamed. By writing Q,,(R

1 = *#I

+ ET;I(R1,

then the expectation

value correct

(4)x

J

= 2WM

- (4)”

w9 to order c2 is (13)

where

As an example the drpole moment of LiH at R = 3 -015 ao was calculated using the MC SCF(4) wavefunction in table 5 The electromc expectation values (m atomrc units), wrth the ongm at the bond rmdpoint and themolecule along the Z-axis, are (z>, = -0.7054 f 0.0022, (zlM = -0.6951 f 0.0027, LZ), =-O-685 f 0.005. On addition of the nuclear drpole a predicted value of 2.33 i 0.005 is obtained, which is to be compared wrth the estimated experimental value of 2.294 [ 161. The quoted error is very uncertain and may be at least 50% in error despite the great length of the cakulation. The large population of 1000 configurations was used specrfically to reduce the error in the dipole, the total time for the calculation being approximately 18 IBM3081 hours [ 171. The introduction of control or antithetic variates [17], if possible for the ST simulation as it is for the variational calculation, could be very advantageous in the calculation of such propertres.

5. Conclusions The most important conclusions appear from the comparison of tables 3 and 4. In any molecule where there is substantial configuration mixmg due to a near degeneracy of another electronic configuration it will be necessary to use MC SCF orbital wavefunctions to 262

18 Jalluary1985

obtain nodal surfaces of reasonable accuracy. Tlus wiIl certainly apply to molecules away from equilibrium geometry. Our investigations on the optimisation of the correlatmg factor were d&appointing, In principle rt is possible to write the exact fared node wavefunction as a product of an orbital function and a correlation factor but so far little progress has been made beyond using just one term in the correlating factor. It may well be that the correlating factor has to be very different in nature depending upon the position of fi and 3 in the molecule. Such an effect may be achieved by using terms similar to Boys and Handy in theu calculatrons on Ne and EH [7]. Our attempt to use the CoIIeSaIvettr factor was not very successful and it appears to be more complex and less satisfactory than the Jastrow factor. Clearly this is an area where much work remains to be done; the gains to be made are very substantial. We may also comment on the current status of ab initio and GFMC calculatrons of very high accuracy by comparing the results in table 5 and our best variational calculation on J_iH [ 161 using 168 Gaussian basis functions with an energy of -8.0690 hartree. The MC SCF(2) ST result is below this energy by approximately 3 times the quoted error in the energy and required approximately the same time as the CI calculation (roughly 9 h IBM3081 time). For further comments on this topic see ref [ 16]_ Apart from further theoretical developments in the actual GFMC algorithms, such as the recent nodal relaxation method of Ceperley [3,18], the aim must be to develop accurate compact trial wavefunctions whose variational energy achieves at least 80% of the correlation energy for these simple systems at equilibrium geometry_ If this can be achieved then the evidence of these and other calculations is that GFMC can recover the remaining correlation energy quite efficiently. So GFMC could be very useful for obtaining some benchmark results, such as the absolute energes of some transition states. It remains that the calculations are

very expensive especially when there are many core electrons. However the algorithms are extremely well suited to vector processors such’& the CFUY, indeed if optimally coded the largest calculation on LiH should require less than one hour of CPU on a CRAY-1 S [13,15].

Volume 113. number 3

Acknowledgement The authors acknowledge preprints and commurdcattons from Dr. D.M. Ceperley, and also stimulating discussrons with Professor

18 January 1985

CHEMXAL PHYSICS LElTERS

B J. Alder and Dr. B-T.

Sutcliffe. RJH acknowledges SERC, UK support. References Reynolds, D.M. Ceperky, BJ. Alder and WA. Lester, J. Chem. Phys. 77 (1982) 5593. 121D.M. Ceperley, J. Comp. Phys. 51 (19B3) 404. [31 D-M Cepcrley and B-J. Alder, J. Chem Phys., submitted for publication. 143 D M. Ceperfey and M.H. Kales, in: Monte Carlo methods ~1 statistical physics, ed. h Bmder (Springer, Berlin,

El1 P-J.

1979). ISI J.W. Moskowitz, hE. Schmidt, M.A. Lee and MR. Kales, f. Chem. Phys 77 (1982) 349. [61 R. Jastrow, Phys. Rev. 98 (1955) 1479.

[7] NC. Handy and S.F. BOYS,Proc. Roy Sot. A310 (1969) 63;X311(1969) 309. [8] R Cdle and 0. Salvetti, Theoret. Chim. Acta 37 (1975) [9]

329_ R C&e

and 0. Salvetti, Theoret.

Cbim

Acta 53 (1979)

Amoral and R. McWeeny, Theoret. Chim. Acta 64 (1983) 171. [ 111 J.W. Moskowita. KE. Schmidt, MA. Lee and M.H. [IO]

&.V.

Kales, J. Chcm. Phys. 76 (1982) 1064. 1121 E-A. Hylleraas. 2. Physik 54 (1929) 347;60 (1930) 624; 65 (1930) 209. [ 131 CL. Pekexis, Phys- Rev. 112 (1958) 1649. [ 141 F. Mentch and J.B. Anderson, J_ Chem Phys. 80 (1984)

2675. [ 151 CF. Bunge,Phys. Rev. Al4 (1976)

1965. [ 161 NC. Handy, R J. Harrison, PJ. Knowles and H.F. Schaefer III. J. Phys. Chem., to be pub~hed_ [ 171 J-M. HammersIey and D.C. Handscomb, Monte C&lo methods (Me&men, London, 1964). [ 181 D-M. Ceperley, in: Recent progress ii. many-body thcories. eds. J. Zabolitzky et al. (Springer, Berlin, 1981).

263