Faraday Symp. Chem. SOC.,1984 19 17-37 The Lennard-Jones Lecture The Value of Very-high-accuracy Calculations in Quantum Chemistry By N. C. HANDY University Chemical Laboratory Lensfield Road Cambridge CB2 1EW Receiued 14th January 1985 Five recent topics under research at Cambridge in ab initio quantum chemistry are reported (a) the latest developments in methods for full CI calculations (b) the convergence of the Marller-Plesset perturbation-theory series at the RHF and UHF level using H,O and NH as examples with energies up to 24th order (c) a new way to calculate Marller-Plesset energy gradients and second derivatives the latter needing only the solution of the first-order CPHF equations (d) some calculations using the quantum Monte Carlo method compared with standard CI calculations and (e) high-accuracy calculations on Be, Be and Be to demonstrate the value of such research.I. INTRODUCTION This paper reviews some of the progress made in the development of methods for the accurate ab initio computation of the properties of molecules by members of the Cambridge University Theoretical Chemistry Department. Some particular applica- tions of these methods are also reported. The members of the department who have been actively involved in this research are Drs R. D. Amos S. M. Colwell R. J. Harrison P. J. Knowles and H-J. Werner and Miss J. E. Rice and Miss K. Somasundram. Much stimulation for this research came from interaction with Prof. H. F. Schaefer’s group at Berkeley and as far as these particular topics are concerned helpful conversations were held with Prof.J. A. Pople R. Ahlrichs P. Siegbahn and B. Alder. In our research we must take care to keep up with and take advantage of the latest hardware advances in modern computer technology. It will be seen in the forthcoming sections that considerable emphasis is placed on the importance of developing methods which take advantage of the ‘vector’ capability of the CRAY-IS computer to which we have had access over the past three years. (In my experience there is always at least a year’s development time before researchers think in the appropriate way to take such advantages.) We also imagine that in the next few years there will be considerable emphasis on the development of programs for the smallest computers (e.g.IBM PC) which for a few thousand dollars can give the C.P.U. power of a VAXl 1/750. Since I have spent much of my research on the development of methods for highly accurate wavefunctions this paper will be largely devoted to such topics. These are now summarised. (a) As a test of the accuracy of many methods it is desirable to have an efficient program for full configuration-interaction (CI) calculations. Much progress has been made on this problem in recent years and in section I1 a new method for these calculations is introduced which takes advantage of vectorisation procedures. Besides being useful in its own right this program can be the central algorithm of any complete active-space self-consistent-field (CASSCF) program.17 THE LENNARD-JONES LECTURE (b) The wide usage of the GAUSSIAN-80 and GAUSSIAN-82 packages and their reliance on Marller-Plesset perturbation theory for the calculation of correlation energies forces one to enquire of the convergence of this perturbation series. Using the full CI program it is possible to examine the convergence for some selected cases and these new results are discussed in section IV. In this context the accuracy of single-plus-double excitation CI(SDCI) many-body perturbation theory (MBPT) and multi-reference SDCI (MRSDCI) calculations in general are also discussed in this section. (c) We have recently observed that the way in which energy gradients and higher derivatives are evaluated for energies which include electron correlation can be improved.In section V it is shown that second-order Marller-Plesset second derivatives can be evaluated without solving the 3N2 second-order CPHF equations (N is the number of nuclei). (d) The recent innovations in the area of quantum Monte Carlo calculations for the electronic energies of small systems has led to some claims that this approach may be the answer to many of our problems. In section VI this method is compared with more standard ab initio methods and future areas for development of the method are discussed. (e) Finally as a demonstration of the very-high-accuracy type of calculations that we are now able to perform results are reported for the binding energies of Be, Be and Be,. Theoreticians predicted the existence of a bound 'Cg state of Be with several vibrational levels and this has now been confirmed experimentally.11. THE FULL-CI PROCEDURE This is a very simple problem to state given a set of orthonormal orbitals 41&...4m what are the energies obtained for the lowest states when the wavefunction is expressed as a linear combination of all possible determinants which can be formed from these orbitals? The length of the expansion is reduced if configuration state functions (CSF) are used instead of determinants but the results will be the same. Full CI calculations involve differently structured programs to the usual variational programs because all the two-electron integrals may be held in core. However the length of the CI vector rapidly becomes very long.If A4 is the number of orbitals and N and Npare the number of aand /3 electrons the number of determinants is When the Sz eigenvalue is zero it is possible to show that the combinations &(@% @4J & &(@BA @%) (2) span singlet and triplet states respectively. In this notation OAand denote a string of orbitals. The following table shows approximate relationships between the number of CSF and the number of determinants which arise in typical cases szeigenvalue state expansion function no. of determinants/ no. of CSF 4 1 2 1 1.2 1 2.6 1 1.8 1 1.4 1 1.2 1 N. C. HANDY In the full-CI programs we have developed we have used determinants as the expansion functions. If the combinations (2) are used this means that the length of the determinant vector will be no more than twice as long as the CSF vector except for the doublet case.This probably means that as far as full-CI is concerned the use of determinants is not a major disadvantage. The first full-CI program which was designed for the largest cases used the combinations (2)1-3and used the Cooper-Nesbet algorithm4 to determine the lowest eigenvector. This relies upon Acz = a1/(Hzz-E) (3) as the update procedure performed sequentially. It needs only the current eigenvector in the memory. The program was very short easy to debug and has now been used by various groups. The program is not totally vectorisable because it needs random access to the integral list and is fairly slow. It relies upon a lexical ordering of the determinants; this is possible because the orbital occupancy can give the /Iand a string labels which then give the determinant label.The largest case we have considered is a calculation on HF involving 1880346 combinations,2 or 944348 CSF in C, symmetry. This is probably the best program available for large full-CI calculations. It is possible to extend the program when it is not possible to hold even the one vector in the memory by the standard ‘paging’ techniques but this makes the program slower. In principal therefore there is no limit to the size of calculation this program can perform. In a recent publication5 Siegbahn showed how it is possible to vectorise completely the full-CI algorithm. The method relies on the Davidson algorithm6 and therefore needs two vectors in the memory.For maximum efficiency paging should not be considered. Essentially his idea was to write the evaluation of the cvector as a matrix product or where the coupling coefficients r& can be written in terms of the unitary group one-particle generators :7 which may be written as r;;Ll= k x(11 Eij I K)(Kl Ekl I J>-$(I1Ei I J> (7) K on the introduction of the resolution of the identity. The difficult term is the first and its contribution to oImay be written as This formulation leads to the scheme 0= Tr(y I-D) THE LENNARD-JONES LECTURE kl ij which can be efficiently vectorised with the sum over K as the outermost loop provided the one-electron coupling coefficients yGK (given K all J) are available from a tape or can be rapidly evaluated.Siegbahn based his program on the graphical unitary group approach (GUGA) CI algorithm,8-10 which meant that the ygK had to be predetermined and held on a tape. This meant that the method rapidly became input- output bound. The largest case Siegbahn reported was one involving 30000 CSF. If instead of CSF single determinants are used the evaluation of the ygK is trivial; they each have value 0 or 1. We have recently reported the development of this determinant CI algorithm.ll The determinants are addressed cia a lexical ordering scheme which means that the whole procedure is vectorisable on the CRAY and the algorithm is one of the most ‘vector’4ntensive in existence. Full details are given in ref.(1 1). Because of the structure of the Davidson iterative algorithm if one starts with a guessed vector which has the correct space-spin symmetry (ie.a linear combination of a small number of determinants) then that symmetry will be exactly maintained throughout the procedure provided that the diagonal elements HI are replaced by the average diagonal energy l/MZyl HI,of all those A4 determinants Orwhich have the same orbitals singly occupied. This aspect is now fully vectorised as well. On the CRAY-IS with 106 words of memory it is possible to consider calculations with 300000 determinants d(@i @&) and for such a case an iteration takes 18 C.P.U.s with no input-output requirements. This program is now available from the authors.An important by-product of this algorithm is that it forms the central part of a new CASSCF12 program developed by Werner and Knowles.13 This latest multiconfigur- ation SCF (MCSCF) development has examined carefully the coupling between orbital and CT coefficients and the evidence to date is that the new program has excellent convergence features and can deal with systems involving > 100000 CSF which means that molecules involving transition-metal atoms can be reasonably treated calculations on Fe0(5A) are being reported. The conclusion of this section is that the last five years have seen the development of highly efficient vectorised codes for full-CI calculations which besides being useful on their own for theoretical investigations on convergence form the central part of good CASSCF programs.It should also be added at this stage that the GUGA CI program is suitable for full-CI calculations when the number of electrons is very small and the number of basis functions is large. Such a calculation is not really appropriate for GUGA but we have reported calculations on Be and (H,).l* 111. THE GENERATION OF TERMS IN THE M0LLER-PLESSET PERTURBATION SERIES We may use the full-CI program previously described to examine the convergence of the Msller-Plesset perturbation series.I5 In the framework when the zeroth-order hamiltonian is the sum of one-electron Fock operators and we use standard Rayleigh-Schrodinger perturbation theory then many-body perturbation theory (MBPT)l67l7and Mdler-Plesset (MP) theory involve the same series the names only distinguishing the way in which the early terms in the series are evaluated.For N. C. HANDY higher-order terms there seems only one sensible way to proceed and that is outlined by Bartlett and Brandasla and reproduced here. Ho is defined as Ho = xF(i). (14) The eigenfunctions of Hoare all the determinants OIformed from the set of all orbitals involved and the eigenvalues are the sum of the orbital energies. Intermediate normalisation is used (yoI Vk) = 0 and the perturbation equations are Ho Yo = Eo Yo (15) k (HO-EO) lyk+(H1-El) wk-1-c Er Vk-r = (16) where H = Ho+HI. The full-CI program evaluates vectors o = Hc and thus in this case it can evaluate Ok = Hvk-l (17) where the ly are expanded in terms of all the determinants V/k = xdrk(Dr* (18) Successive perturbation energies can therefore be evaluated through Ek = (YO I Ok-1) (19) and i.e.to evaluate the kth perturbation energy needs a time equivalent to k full-CI iterations. Thus this development is a trivial change to the full-CI program for closed-shell restricted-Hartree-Fock (RHF) calculations. Of course MP theory may also be used in the unrestricted-Hartree-Fock (UHF) frame~0rk.l~ The full-CI program is amended to have available the two-electron integrals (aaI aa) (aaPP) (JP I aa) and (apI PP) in the memory. There are no other changes necessary to enable the program to generate nth order MP UHF energies. Laidig et uf.,O have also reported MBPT(n) series calculations based on their GUGA programs.IV. RESULTS AND CONCLUSIONS FROM THE FULL-CI AND PERTURBATION SERIES CALCULATIONS Several full-CI calculations have now been performed on systems with four or more electrons BH (2 -[+polarisation basis) H,O (2 -< basis) NH (2 -[basis frozen core and virtual orbital) HF (2 -[+polarisation basis frozen core and virtual orbital),,T and some ionised states (,B, 2Al and ,B,) and excited states (,B1 ,B2 ,A1 and ,A,) of H,0.21 Whilst it is obviously important to have many such calculations available in particular with as many basis functions and electrons as possible so that there are many types of higher excitations this author believes that probably many of the conclusions which can be deduced from a careful examination of the results for H,O (lA,) at three C, geometries (6 = 8, r = re l.5re and 2.0rJ can be extrapolated to many other systems.The significance of these calculations is that at the full-CI level the coefficient THE LENNARD-JONES LECTURE Table 1. Calculations on H,O using a 2-[ basis set. All orbitals active. Three bond lengths. C, symmetry ~ ~~ re 1.5re 2.0re SCF" -76.009 838 -75.803 529 -75.595 180 SDCI" -76.150 015 -75.992 140 -75.844 817 (0.9 78 7) (0.9462) (0.8714) SDTQCI" -76.157 603 -76.013 418 -75.900 896 (0.9 754) (0.9247) (0.7726) full CI" -76.157 866 -76.014 521 -75.905 247 (0.9 7 5 3) (0.923 3) (0.7636) Davidson extrapolationa -76.155 912 -76.01 1 867 -75.904 874 MBPT(2)b -76.149 315 -75.994 577 -75.852 460 MBPT( 3)b -76.150 707 -75.989 39 1 -75.834 803 MBPT(4)b -76.156 876 -76.008 395 -75.888 867 M BPT( 5)" -76.157 056 -76.009 771 -75.889 199 [2 11" -76.157 358 -76.012 555 -75.905 780 CCSDb -76.156 076 -76.008 931 -75.895 913 CCSD +T(4)b -76.157 438 -76.012 858 -75.907 996 CASSCF-6C -76.026 113 -75.842 122 -75.677 112 CASSCF-7 -76.065 100 -75.927 258 -75.828 416 CASSCF-8 -76.099 399 -75.955 930 -75.845 671 CASSCF-9 -76.132 645 -75.982 052 -75.866 155 MRCISD-6" -76.152 546 -76.004 199 -75.870 295 MRCISD-7 -76.155 929 -76.012 464 -75.903 349 MRCISD-8 -76.157 488 -76.014 032 -75.904 708 MRCISD-9 -76.157 783 -76.014 349 -75.905 105 " From Handy and coworkers ref.(1)-(3). From Bartlett and coworkers ref. (17) and (20). From Shavitt and coworkers ref.(22). of the dominant determinant is 0.975 0.923 and 0.763 respectively and thus such calculations will show clearly the effects of configuration mixing which are met in so many calculations especially away from equilibrium geometry. The size of the calculations was 502355 determinant combinations (2) or 256473 CSF. In table 1 many results are collated by various methods using the same basis set and geometries. We now discuss the results under the various headings in that table. Acknowledgement is made here to Prof. Bartlett and Shavitt whose results are q~0ted.l~~ 20* 22 SCF (self-consistent field). Note that at re l.5re and 2.0r the correlation energies are 0.148028 0.210992 and 0.310067 hartree respectively. SDCI. The singles and doubles CI calculation picks up 95 89 and 80% of the correlation energy.This is clearly an inferior method for potential-energy surfaces. SDTQCI. With triples and quadruples added although all results are better than 98.5% of the correlation energy the error at 2re is 0.00435 hartree. 17678 CSF are involved. Davidson correction This extrapolation gives the best result for 2re but must be regarded as fortuitous! MBPT(4). This is the limit of general practical calculations of MP or MBPT. The errors are 0.000990 0.006 126 and 0.01638 hartree. These results demonstrate the weakness of this method for sections of surfaces involving the breaking of bonds. N. C. HANDY 0.012 0.013 0.014 0.015 0.016 0.017 0.0 18 0.019 11111111111 3 4 5 6 7 8 9 1011 121314 Fig.1. Convergence of MP-n series for H,O at Re:SCF =-75.8884 2nd order = -76.0093 exact = -76.01851 5. CCSD. The coupled-cluster SD calculations have errors of 0.00179 0.00559 and 0.009334 hartree; the errors are comparable to MBPT(4). CASSCF. Even with 9 active orbitals only 83% of the correlation energy at Y is obtained. MRCISD. These SDCI calculations using the CASSCF functions as the reference set are by far the most uniformly accurate. Even SD-CI out of CASSCF with 7 active orbitals shows a constant error of 0.001 hartree at the three points (7096 CSF in these calculations). Some of these results are given in fig. 1 which shows very clearly that the best results are achieved by MRCISD calculations; and the conclusion of this discussion must be that the very best programs are going to be those which involve SD excitations from a carefully selected CASSCF or MCSCF reference space.The most hopeful attempt in this direction to date is due to Werner and Me~er~~ in their contracted CI scheme where SD replacements of the full reference function are included in a CI scheme. This is where the future may lie. The results of the perturbation-theory calculations which are given in table 2 and THE LENNARD-JONES LECTURE Table 2. Marller-Plesset perturbation energies En using 6-21G basis set with a frozen core orbital for re lSre and 2.0re(8 = 107.6') for H20 En re 1.5re 2.0re -0.120 865 -0.166 896 -0.241 643 -0.003 303 -0.002 015 0.006 123 -0.004 849 -0.016 925 -0.046 465 -0.000 488 -0.001 352 -0.004 219 -0.000 435 -0.003 417 -0.012 369 -0.000 076 -0.000 338 -0.001 032 -0.000 048 -0.000 827 -0.004 065 -0.000 012 -0.000 032 0.002 319 -0.000 006 -0.000 189 -0.001 417 -0.000 002 0.000 010 0.002 191 -0.000 00 1 -0.000 032 -0.000 277 10-7 0.000 007 0.000 959 10-7 -0.000 001 0.000 023 10-7 0.000 003 0.000 146 10-8 0.000 002 0.000 025 10-8 0.000 001 -0.000 140 10-9 0.000 001 0.000 043 10-9 10-7 -0.000 159 10-9 1o-6 0.000 087 E21 10-10 10-7 -0.000 095 10-10 10-7 0.000 095 E22 E23 1O-IO 10-8 -0.000 037 10-l1 lo-* 0.000 064 E24 fig.2-4 show that in most cases the series will converge to the exact eigenvalue. The calculations we have performed on H20used a 6-21G basis with a frozen core orbital the latter being to make the size of calculation tractable (61441 determinants and 245025 internal IK) determinants).Indeed the usual applications of MP theory using GAUSSIAN~~~ 25 programs use frozen core orbitals. That the effect of the core orbital does not affect a general discussion of the convergence of the series is indicated from a comparison of the MP-n (n = 2-5) perturbation energies En for the 24 (no frozen core) and these 6-21G calculations the former being due to Laidig et al. for H,O at 2r 2-r 6-21G E2 -0.257 28 -0.241 64 E3 0.017 66 0.006 12 E4 -0.054 06 -0.046 47 -0.000 33 -0.004 22 E5 -0.012 37 E6 The results in table 2 clearly show that the early rapid convergence of the MP-n series is not maintained.In the most favourable case re the energies have converged to within 0.001 hartree with E4,whilst for l.5re this happens at E, and for 2r, at 0.89 0.894 \ 0.895 0.896 ! I I 0.89 7 X 0.898 \ X \ 0.899 0.900 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 Fig. 2. Convergence of MP-n series for H,O at 1.5Re SCF = -75.7072 2nd order = -75.8741 3rd order = -75.8761 exact = -75.899202. 0.789 0.790 0 -791 0.792 0.793 0.794 v, 0.795 -X I I ,I , I ,, I Fig. 3. Convergence of MP-n series for H,O at 2.0Re SCF = -75.4914 2nd order = -75.7730 3rd order = -75.269 t3 4th order =-75.7733 5th order = -75.7776 exact = cn -75.791 269. THE LENNARD-JONES LECTURE 0.01 I 0.009 -0.008 -0.007 -0.006 -0.005 -0.OOL -0.003 -0.002 -0.001 0.001 0.002 0.003r 0.OOL 0.005 Fig.4. Energies of H20(2-c basis) using various methods (a) SDCI (b) MBPT(4) (c) MBPT(5) (d)SD T QCI (e)MRCISD-7 and (f)MBPT(8). Table 3. SCF eigenvalues of the 6-2 1G H,O calculations En re 1Sr 2.0re E2 -1.329 -1.213 -1.199 E3 -0.684 -0.506 -0.504 E4 -0.529 -0.484 -0.415 E5 -0.475 -0.465 -0.409 ‘6 0.261 0.099 -0.012 E7 0.361 0.170 0.033 ‘8 1.207 1.081 1.121 E13.The fact that the convergence is so rapid at re underlines the value of the many MP-4 calculations which are carried out at equilibrium geometries. However away from equlibrium the picture changes. Convergence is reasonable at 1.5re but note that at E5the error is 0.005 hartree which means that it is impractical to achieve 0.001 hartree accuracy.At 2re at eighth order the energy is 0.0038 below the exact value whilst at 22nd order the energy is oscillating by 0.0001. These results at 2re show the dangers of using any Pade approximants.18 These approximants are designed to approximate the energy that will be obtained with a few more terms in the series. For example an approximant based on E3 E4and E5 yields an energy -75.792 704 which is 0.0013 below the exact. The [3,3] approximant yields -75.7993 which is 0.008 below the exact. In this author’s view any approximant which is claimed to be successful at these distorted geometries can only be so by chance. A reason for the slower convergence at 1.5re and 2re can be obtained by observing the SCF eigenvalues in table 3.The two lowest virtual eigenvalues E and E, become N. C. HANDY Table 4. UHF MP-n calculations on NH (,B,) using 6-31G basis (frozen core) re 13 2re -55.530 176 -55.367 729 -55.181 -55.532 247 -55.405 143 -55.381 93 1 -0.085 5 12 -0.062 116 -0.031 541 -0.009 815 -0.01 1 393 -0.006 208 -0.003 612 -0.008 695 -0.002 508 -0.001 192 -0.005 775 -0.001 235 -0.000 463 -0.004 888 -0.000 845 -0.000 200 -0.004 066 -0.000 663 E* -0.000 100 -0.003 505 -0.000 568 CE -55.633 142 -55.505 582 -55.425 497 full CI -55.633 27 -55.526 -55.44 correlation energy RHF 0.103 0.158 0.26 UHF 0.101 0.121 0.06 (9) 0.752 1.661 2.51 0.627 0.628 0.629 0.630 0.631 I 0.632 0.633 I \ 2345678 Fig.5. Convergence of MP-n series for NH at Re:SCF = -55.53225,2nd order =-55.61 776 exact = -55.63327. THE LENNARD-JONES LECTURE 0.487 0.489 0.491 0.493 0.495 0.497 0.499 0.501 0.503 0.505 0.507 0.509 0.51 1 0.512 I I I I 1 I I 3 4 5 6 7 8 9 1 0 Fig. 6. Convergence of MP-n series for NH at 1.5Re SCF = -55.40514 2nd order = -55.46726 3rd order = -55.47865 exact = -55.526. nearly degenerate and indeed E < 0 at 2re. The MP series is so dependent on the size of the denominators that this has to be the principal reason for the convergence difficulties. We deduce that MP-n [or MBPT(n)] calculations ought to be viewed with serious suspicion if any of the occupied or virtual orbitals have the wrong sign for their SCF energies.We also have some preliminary results for UHF MP-n theory. The calculations were performed for NH (2B,),for fixed bond angle 8 = 103.2" and bond lengths re 1.5re and 2re where re = 1.013 A (C2vsymmetry) with a 63 1G basis set and a frozen core orbital. The results are given in table 4 and fig. 5-7. The results are rather different to those for the closed-shell RHF case. In all three cases the convergence appears smooth but it is very slowly convergent at l.5re and 2re. It appears to be sufficiently smooth that some sort of extrapolation may be useful this is being investigated. The results may be considered curious because although UHF picks up 2 23 and 77% of the correlation energy at re 1.5re and 2re the series needs many more terms to find the remaining correlation energy at 2re than at re.The reason of course N. C. HANDY 0.420 \ 0.422 0.423 0.424 0.425 0.426 3 4 5 6 7 8 9 10 Fig. 7. Convergence of MP-n series for NH at 2.0Re SCF = -55.38193 2nd order = -55.41347 exact = -55.439. lies in the very poor spin state found by the UHF calculation as is seen by (S2) = 2.5 at 2r,. In summary preliminary conclusions on these investigations into the convergence of MP-n series are that (a) the RHF results show slow and irregular convergence features away from equilibrium and (b) the UHF results show a much smoother convergence but also a very slow convergence away from equilibrium.Further investigations will be carried out and a complete set of results and discussion given elsewhere. V. AN IMPROVED THEORY FOR THE CALCULATION OF ENERGY GRADIENTS AND HIGHER DERIVATIVES Because the calculation of energy gradients and second derivatives is becoming an important every day tool for the ab initio chemist it is useful at this stage to underline a recent development in the theory for these calculations which has particular relevance for the calculation of CI gradients and MP-2 first and second derivatives. The possibility of calculating MP-2 second derivatives analytically is important 30 THE LENNARD-JONES LECTURE .~~ because as Pople et ~1say 'MP-2 accounts for 70% of the difference between HF theory and experiment yielding frequencies which differ from experiment by only 10-90 cm-l '.Under a movement of a nucleus the orbitals and the basis-function integrals 28 4i 4 4i+a CUai4a (21) + (aBI Y@ CaB I 74+4og I 74' (22) and for SCF or MCSCF wavefunctions the energy gradient can be straightforwardly evaluated29 with the usual notations for the integrals and the basis-function density matrices y r. In the derivation of eqn (23) the uUican be eliminated because the energy E has been optimised with respect to all parameters. For the MP-2 energy26 occ virt EM, =-a z E (iJ' II +ESCF ij ab (EafEb-ei-Ej) the energy is not optimised with respect to all parameters but we do know how the orbitals change because the uai are solutions of the SCF coupled-Hartree-Fock (CPHF) equations.28 The non-redundant uai (a = virtual i = occupied) are the explicit solution of these equations A ux = b" which in detail are ((er-ei) uFi+zz [4(irI sk)-(isI rk)-(ikI rs)]u&} sk =-&fr+S& ei++I z S&[4(ir(sk)-(isIrk)-(ikIrs)] (26) sk for closed-shell systems.The remaining ux may be expressed26 in terms of the solution of the CPHF equations but it is desirable to use formulae which do not explicitly evaluate uij (i and j occupied) and Uab (a and b virtual) because these depend on (ei-ej)-l and (e -eb)-l respectively. .~~ Pople et ~1have given a formula for the first derivative of EMp2,which can be rewritten in the form where Tx depends upon derivative integrals and Waiis given by Standard GAUSSIAN programs solve the CPHF equations for the ux (3N of them N = number of nuclei) and then calculate dE/dX according to eqn (27).N. C. HANDY However eqn (27) can be rewritten as dE -= U/TuX+P dX = U/rA-lb"+ T" (31) where ATZ= w. (33) If the gradient is evaluated through eqn (32) then only one set of simultaneous eqn (33) have to be solved for 2,instead of 3N; more importantly no derivative integral has to be stored because whether it contributes to T" or hx,it can immediately be multiplied by its appropriate factor and entered directly into the gradient vector. For second derivatives the advantages of this development are even more significant. The formula for the second derivative of the MP-2 energy has to take the form To evaluate this the first term which involves the solutions of uxY of the second order CPHF equations [(3w2of them] can be reduced by a similar analysis to eqn (30)-(33) to u/TuXY = FbXY (35) where bxY is the right-hand side of the second-order CPHF equations.Therefore to evaluate a2E/aX(3Y only the first-order CPHF equations have to be solved. The evaluation of MP-2 second derivatives should therefore be a relatively straightforward matter and we are currently implementing this procedure. These arguments also apply to the calculation of CI gradients and the details may be found in ref. (30). It also applies to the calculation of molecular properties. For example for correlated wavefunction calculations where the Hellmann-Feynman theorem is not obeyed it has been suggested by Diercksen and Sadlej31 that more accurate first-order properties may be obtained by calculating the energy derivative.This may be written for the dipole moment and for MP-2 calculations as =-NIP2 --XZai(al~li)-i ~@jb~~~(kI~li)+~x C a$bag(cIxlb) (36) a1 ai ijk ab ij abc which follows from eqn (32). These are to be compared with the formula derived from expectation values by for one-electron Marller-Plesset properties which are therefore different. VI. A VIEW OF QUANTUM MONTE CARL0 AND ITS COMPARISON WITH CONFIGURATION INTERACTION Recently a novel method for the calculation of electronic energies of small molecules was introduced by the authors of ref. (33)-(35) called the quantum Monte Carlo method. A few groups are now actively pursuing research in this area and it is of value to compare the current status of these calculations with the more standard ab initiocalculations.We start by briefly referring to our recent CI calculation on LiH (1Z+),36which was an attempt to answer the following question 'Using gaussian basis sets and CI how accurate an energy can be obtained for LiH?' THE LENNARD-JONES LECTURE Table 5. Energy and dipole moment of LiH at r = 3.015 bohr method E/ hartree PP ref. SCF (STO) -7.987 3 13 -6.001 5 41 SCF (numerical) -7.987 34 -6.002 5 42 SCF (numerical) -7.987 354 -6.002 5 31 NO-CI (ETO) -8.060 62 -5.965 43 PNO-CI (GTO) -8.064 705 -5.837 40 CI (ETO) -8.065 53 44 CI (GTO) -8.069 04 -5.86 36 -exact -8.070 49 -5.83 At the SCF level it is now possible to achieve an accuracy to 5 or 6 decimal places using the numerical procedure of Laaksonen et aL3' for the solution of Poisson-like equations.(It is interesting to that there are 'spurious' nodes in these numerical Hartree-Fock orbitals at long range just as in the Hartree-Fock orbitals for atoms. These are probably due to the long-range behaviour of these orbitals which is dominated by exchange effects.39) Up to 1983 the best CI calculation at Y had an error of 0.005 hartree (Meyer and Ro~mus,~~ PNO-CI calculation using a 9s5p4d/7s2pld basis set). The energies and dipole moments of LiH are given in table 5 for various calculations. We chose the basis set by optimising it for Li+ and H- and found that with 1 Is 1 lp 8d 3f and lg for Li and lOs 5p 3dand 2f for H the energy errors were as follows Li+ (0.00055) Li (0.00058) H-(0.00028) H (0.0000035).These calculations demonstrate the well known fact that more angular functions are required to represent the cusp in Li+ than in H-as far as the energy is concerned. The CI calculation involved all single and double replacements (1 3201 5 CSF) from a CASSCF with the lo 20 30 40 and In orbitals active. The resulting CI energy of -8.06904 hartree was in error by 145 x which can be attributed as 90 x basis-set error and 55 x quadruples error. The time for this calculations was 9 h on the IBM 3081D of which the 4-index transformation took 5 h. The conclusion to be drawn from this rather tedious calculation is that the basis-set problem remains with us at the CI level in a very much worse way than at the SCF level.To give a fair comparison with the quantum Monte Carlo method one of the authors of the LiH calculation (R. J. Harrison) undertook these calculations. He used the short time GFMC algorithm of Reynolds et al. with the fixed-node approximation. The method is well described in their paper. The equation which is simulated is where f = qh,vT and where at convergence q4 = dexact.The success of the method depends upon the ability to choose a good trial function yT because (a) the nodes off are the same as those of yT in this procedure and (b) the variance is crucially dependent on ry, as is evident from the term HyT/t//T in the equation and the fact that the energy E is calculated through E= fedV 1sfdV.s VT N. C. HANDY 33 Table 6. GFMC results for He using various forms for tptrial time-method trial wavefunctiona energy points steps blocks variational 1 term Jastrowb -2.883 6 --GFMC 1 term Jastrow -2.902 81 +O.OOO 83 100 50 900 variational 6 term Jastrowb -2.901 82 GFMC 6 term Jastrow -2.904 17+0.000 21 100 50 700 variational 6 term Hylleraas" -2.903 24 GFMC 6 term Hylleraas -2.903 9 & 0.000 4 200 200 12 GFMC 6 term Hylleraas -2.903 74f0.000 05 200 200 1156 a The orbital part = exp [ -1.6875(r +r,)] for the Jastrow function. From N. C. Handy (Ph.D. Thesis). Ref. (46). At the end it is hoped that the distribution of points represents the density f=fexact vT,but one has to be very careful about how long it takes to reach this ideal distribution.Years ago one of the justifications for the trans-correlated method45 was that if CO was an exact wavefunction then HCO = ECO (39) implied that (C-lHC) O = EO (40) and thus that the use of the non-hermitian operator CfHCwas not too important. However it is probably true to say that sufficiently accurate functions CO were never obtained to realise fully the potential of the method. In the GFMC method vTis represented as CO,where O is a linear combination of determinants and C is a Jastrow-like factor. For atoms Jastrow factors are or with similar expressions for molecules. C has no nodes. In table 6 some results are given for He. The results for the Hylleraas trial function show that with a large amount of time high-accuracy results can be obtained.However the results obtained with increasingly accurate Jastrow functions only marginally reduce the variance and question whether the Jastrow form is most appropriate. The results for LiH and Be are summarised in table 7. All of these used the simple Jastrow factor.41 For LiH the variational CO (SCFO) energy gives 47% of the correlation energy and the corresponding result with a 4-term MCSCF O yields 67% of the correlation energy. The best GFMC result is twice as accurate as the CI result for approximately the same cost in computing time. The disappointing aspect is that the variance is not reduced by using an MCSCF O,and is only reduced by extending the run time i.e.the number of blocks. For Be it is not surprising that the MCSCF results are far more accurate than those achieved with an SCF @ (98% instead of 88 % of the correlation energy for the best GFMC calculations). This is due to the considerable configuration mixing and thus 2 bAR THE LENNARD-JONES LECTURE Table 7. GFMC results for LiH (re = 3.015) and Be time-method energy points steps blocks LiH (i) SCF 0 -7.986 50 C Q (variational) -8.026 4 & 0.000 4 1000 20 821 GFMC (Z = 0.01) -8.069 6 f0.000 7 1000 50 99 (ii) MCSCF Q -8.013 64 --C @ (variational) -8.042 8 f0.000 5 1000 50 657 GFMC (Z = 0.01) -8.069 73 & 0.000 26 1000 50 887 (iii) CI -8.069 0 --(iv) exact -8.070 49 Be (i) SCF -14.572 4 --C 0(variational) -14.608 9&0.000 8 800 100 135 GFMC (Z = 0.01) -14.656 4+0.000 5 800 50 347 (ii) MCSCF 0 -14.628 0 -C Q (variational) -14.647 8 +O.OOO 7 800 100 124 GFMC (Z = 0.01) -14.665 1 & 0.000 4 800 100 270 (iii) exact -14.667 3 the change in the nodal surfaces.However again we notice that the variance is not reduced. The probable conclusions of these GFMC calculations and similar calculations by other workers are as follows. (a)The dominant configurations should be included in because they can affect the nature of the nodal surfaces. Their inclusion does not increase the time cost very much. (b) The Jastrow factor is obviously very important. However it is not clear how to determine the non-linear parameters. Does the shape near the electron4ectron cusp depend upon the local electron density? Is the form given by eqn (41) appropriate? (c) The principal effort must be to reduce the variance.(d) Therefore until either the Jastrow factor can be substantially improved or the variance reduced in some other way these calculations are likely to remain both enormously expensive and restricted to a very small class of molecules. VI. CALCULATIONS OF HIGH ACCURACY ON SMALL BERYLLIUM CLUSTERS There have been many calculations in recent years on small Be clusters [see e.g.ref. (47) and (48)] and the purpose here will be to demonstrate that high-accuracy results can be achieved on Be, Be and Be which are at variance with calculations of less accuracy. Be2 There is now complete agreement between theoreticians and experimentalists on the potential-energy curve for the ground state (Xxi).Indeed it is one of the successes of high-accuracy ab initio calculations that they predicted a 2.3 kcal minimum at 4.75 bohr before the experimentalists confirmed that Be existed with several vibrational levels.The definitive calculation on Be, which used a large basis set and sufficient CSF is due to Lengsfield et aZ.49Earlier calculations all suffer from a lack of one or both N. C. HANDY of these and now become largely historical. A complete list of these earlier works may be found in ref. (14). The results of these many ab initio calculations on Be suggest the following conclusions. First the SCF curve is purely repulsive and secondly SDCI calculations predict a shallow (< 1 kcal mol-l) minimum around 8.5 bohr.These statements are independent of basis set. For a 7s+ 3p+ Id basis set a shallow minimum is also predicted by SCEP CEPA-3 CEPA-1 and also CCSD and CCSDT-1. Lee and Bartlett50 have presented a clear explanation of these results showing in particular that neglect of triple and quadruple excitations leads to too deep a well at short distances (eg. CEPA-0) whilst inclusion of only doubles and quadruples produces too repulsive a curve at short distances. Thus MBPT(4) and full CI (SDTQCI) should produce satisfactory results and indeed they both predict a well of ca.0.8 kcal mol-1 around 5.0 bohr. A 7s +4p +2d full-CI calculation yields an interaction energy of 1.06 kcal/mol-l at 4.75 bohr. Our best full-CI result (8s+5p+2d+ If basis set) gave 1.86 kcal mol-1 at 4.75 bohr and is in complete agreement with the more restricted MRSDCI results of Lengsfield et al.49We have argued that further basis-set extension must lead to a final result in the region of 2.2-2.3 kcal mol-1 for the interaction energy with g functions estimated to contribute 0.15 kcal mol-l.These agree with the experimental results of Bondybey and English,51 who predict re =4.65 bohr and D in the range 2.14-2.29 kcal mol-l. Be has therefore turned out to be an excellent molecule to demonstrate the need for high-accuracy calculations which need both a large number of CSF and a large basis set. Be3 Increasing sp hybridisation means that there is an increase in the stability for the sequence Be, Be and Be,.However because of the above arguments for Be we would expect that a large calculation involving many excitations will be necessary to achieve the correct dissociation energy for Be,. At the SCF level Be (D3h)seems not to be bound although Whiteside et al.48found D,to be 1.5 kcal mol-l in a 6-31G* basis set. The SCF curve approaching from dissociation passes through a maximum before approaching the minimum. Whiteside et al. found in a MP(4) (SDQ) calculation that the binding energy was 6.0 kcal mol-1 for r =2.23 A i.e. that electron-correlation effects are significant. The best calculation we performed on Be involved a 7s+ 3p+ 2d basis set with an initial CASSCF calculation for those orbitals correlating with 2s 2px and 2py orbitals on the atoms (z out of the plane).A reference set was then chosen as all SD CSF with orbitals correlating with 2s 2px 2p and 2p,. Finally an SDCI out of this reference set generated 784498 CSF in C, symmetry. A binding energy of 19 kcal mol-l was calculated and unlike SCF there is no maximum in the potential curve. A reasonable basis-set extension would predict a binding energy of > 25 kcal mol-l. Be* The situation for Be is rather different because SCF calculations predict Be to be bound by 36 kcal mol-l. 52 Calculations for the correlation energy which include polarisation functions in the basis set show that it has a considerable effect on D,. Whiteside et al.48give 56.0 kcal mo1-I (2.126 A G).The best result of Bauschlicher et al.52is for an SDCI calculation and gives a 64 kcal mo1-1 D,at 3.92 bohr with the Davidson correction contributing 19 kcal mol-l to this result.The best calculations performed by us5 used a TZ+DP basis set for SDCI+ Davidson correction gave a binding enery of 66.1 kcal mol-l. As a confirmation of this 2-2 THE LENNARD-JONES LECTURE result CEPA-l gave 64.8 kcal mol-1 at the same geometry G,re = 3.92 bohr. A more complete basis-set calculation could easily increase D,to 75 kcal mol-l. SUMMARY OF Be, Be AND Be It is interesting to place the results for this sequence in a table cluster DJkcal mol r,/bohr symmetry 2.3 4.75 25 4.22 75 3.92 This shows the shortening of the bond as the sp polarisation effects increase. The earlier results in the literature which we have summarised above show the great care with which these calculations have to be performed before reliable results are obtained.N. C. Handy Chem. Phys. Lett. 1980,74,280. R. J. Harrison and N. C. Handy Chem. Phys. Lett. 1983,95 386. P. Saxe H. F. Schaefer and N. C. Handy Chem. Phys. Lett. 1981,79 202. R. K. Nesbet J. Chem. Phys. 1965,43 317. P. E. M. Siegbahn Chem. Phys. Lett. 1984 109 417. E. R. Davidson J. Comput. Phys. 1975 17 84. 'I J. Paldus J. Chem. Phys. 1974 61 5321. I. Shavitt Znt. J. Quantum. Chem. Symp. 1977 11 131; 1978 12 5. O B. R. Brooks and H. F. Schaefer J. Chem. Phys. 1979,70 5092. lo P. Saxe D. J. Fox H. F. Schaefer and N. C. Handy J. Chem. Phys. 1982 77 5584. l1 P. J. Knowles and N. C. Handy Chem.Phys. Lett. 1984 111 315. l2 B. 0. Roos P. E. M. Siegbahn and P. R. Taylor Chem. Phys. 1980,48 157. l3 P. J. Knowles and H-J. Werner Chem. Phys. Lett. 1985 115 259. l4 R.J. Harrison and N. C. Handy Chem. Phys. Lett. 1983,98 97. l5 C. Mnrller and M. S. Plesset Phys. Rev. 1934 46 618. S. Wilson Theoretical Chemistry (Spec. Period. Rep. The Chemical Society London 1981) vol. 4 p. 1. R. J. Bartlett H. Sekino and G. D. Purvis Chem. Phys. Lett. 1983 98 66. R. J. Bartlett and E. Brandas J. Chem. Phys. 1972 56 5467. lo J. S. Binkley and J. A. Pople Int. J. Quantum Chem. 1975 9 229. 2o W. D. Laidig G. Fitzgerald and R. J. Bartlett Chem. Phys. Lett. 1985 113 151. 21 K. Hirao and Y. Hatano Chem. Phys. Lett. 1984 111 533. 22 F. B.Brown I. Shavitt and R. Shepard Chem. Phys. Lett. 1984 105 363. 23 H-J. Werner and E. A. Reinsch J. Chem. Phys. 1982 76 3144; see also W. Meyer in Modern Theoretical Chemistry ed. H. F. Schaefer (Plenum New York 1977). 24 J. S. Binkley R. A. Whiteside R. Krishnan R. Seeger D. J. Defrees H. B. Schlegel S. Topiol L. R. Kahn and J. A. Pople Quantum Chemistry Program Exchange 1981 13,406. 25 J. S. Binkley M. J. Frisch D. J. Defrees R. Krishnan R. A. Whiteside R. Seeger H. B. Schlegel and J. A. Pople Gaussian-82 (Carnegie-Mellon University Pittsburgh). 26 J. A. Pople R. Krishnan J. B. Schlegel and J. S. Binkley Int J. Quantum Chem. 1979 S13 225. 27 P. Pulay Mol. Phys. 1969 17 197. 28 J. Gerratt and I. M. Mills J. Chem. Phys. 1968 49 1719. 20 J. D.Goddard N. C. Handy and H. F. Schaefer J. Chem. Phys. 1979 71 1525. 30 N. C. Handy and H. F. Schaefer J. Chem. Phys. 1984,81 5031. 31 G. H. F. Diercksen and A. J. Sadlej J. Chem. Phys. 1981 75 1253. 32 R. D. Amos Chem. Phys. Lett. 1980 73 602. 33 J. B. Anderson J. Chem. Phys. 1975,63 1499. 34 J. W. Moskowitz K. E. Schmidt M. A. Lee and M. H. Kalos J. Chem. Phys. 1982 77 349. 35 P. J. Reynolds D. M. Ceperley B. J. Alder and W. A. Lester J. Chem. Phys. 1982 77 5593. 36 N. C. Handy R. J. Harrison P. J. Knowles and H. F. Schaefer J. Phys. Chem. 1984 88 4852. 37 L. Laaksonen P. Pyykko and D. Sundholm Chem. Phys. Lett. 1983 % 1. N. C. HANDY 37 38 P. Pyykko personal communication. 39 N. C. Handy M. T. Marron and H. J. Silverstone Phys. Rev. 1969 180 45.40 W. Meyer and P. Rosmus J. Chem. Phys. 1975 63 2356. L'' P. E. Cade and W. M. Huo J. Chem. Phys. 1967,47 614. 42 E. A. McCullough Chem. Phys. Lett. 1974 24 55. 43 G. F. Bender and E. R. Davidson J. Phys. Chem. 1966 70 2675. 44 D. W. Bishop and L. M. Cheung J. Chem. Phys. 1983,78 1396. 45 N. C. Handy and S. F. Boys Proc. R. London Soc. A 1969 310 43. 46 E. A. Hylleraas Z. Phys 1929 54 347; 1930 60,624; and 1930 65 209. 47 P. S. Bagus H. F. Schaefer and C. W. Bauschlicher J. Chem. Phys. 1983 78 1390. 4* R. A. Whiteside R. Krishnan J. A. Pople M. Krogh-Jesperson P. von R. Schleyer and G. Wenke J. Comput. Chem. 1980 1 307. 49 B. H. Lengsfield A. D. McClean M. Yoshimine and B. Liu J. Chem. Phys. 1983 79 1891. 50 Y.S. Lee and R. J. Bartlett J. Chem. Phys. 1984 80,4371. 5* \1. E. Bondybey and J. H. English J. Chem. Phys. 1984 80 568. 52 C. W. Bauschlicher P. S. Bagus and B. N. Cox J. Chem. Phys. 1982 77 4032. 53 R. J. Harrison and N. C. Handy to be published.