- Email: [email protected]
A new configuration selection method for configuration interaction calculations
9 September 1994
ELSEVIER
CHEMICAL PHYSICS LETTERS
Chemical Physics Letters 227 (1994) 327-336
A new configuration selection method for configuration interaction calculations L. Visscher a,’, H. DeRaedt b, W.C. Nieuwpoort a aDepartment of Chemical Physics and Material Science Centre, University of Groningen. Nijenborgh 4, 9747 AG Groningen, The Netherlands b Institutefor Theoretical Physics and Material Science Centre, University of Groningen, Nijenborgh 4, 9747AG Groningen, The Netherlands
Received 11 March 1994; in final form 21 June 1994
Abstract The recently proposed stochastic diagonalization method is applied to the ab initio quantum chemistry configuration interaction problem. In this context it can be viewed as a multi-reference CI method with dynamic selection of important configurations. The method is compared with other methods and tested by calculations on a number of small molecular systems for which accurate results are available. A calculation on the Cr, dimer is presented to show the capability of the algorithm to find short expansions of molecular wavefunctions.
1. Introduction
The optimal choice of (reference) configurations in truncated configuration interaction methods has been a subject of discussion over the years [ 11. The well-known configuration interaction singles doubles (CISD) method uses an extension of a given reference space that is systematic and unambiguous. The resulting algorithm can be well vectorized and is used to treat matrices of large dimensions [2]. A disadvantage is the lack of size-extensivity that can only partly be corrected by applying corrections [ 3 ] afterwards. Another disadvantage is the scaling of the method with the number of basis functions involved. Full-C1 (FCI) calculations indicate [ 2 ] that most of the configurational space does not contribute significantly to the ground state wavefunction, yet many ’ Present address: NASA Ames Research Center, MS 230/3, Moffet Field, CA 94035-l 000, USA.
configurations do arise in the standard excitation schemes. The sparsity of the CI matrix is used by methods like CIPSI [ 4 1, MRDCI [ 5 ] and the CI + PT method of Harrison [ 61 that select only the important configurations. Experience has already shown that such methods that combine perturbation theory and variational methods can give reliable results while using only a fraction of the full CI space. The aim of this Letter is therefore not to present a new efficient implementation of such a selected CI method but rather to test a method that originates from a different application field. The novel algorithm that we tested is based on the stochastic diagonalization method [ 71 that was recently developed for general quantum many-body problems. The method was used in applications on the two-dimensional Hubbard model where the size of the Hilbert space may become as large as 1035.The success of the algorithm to sample important states
0009-2614/94/$07.00 0 1994 Elsevier Science B.V. All rights reserved SSDIOOO9-2614(94)00824-S
L. Visscher et al. / Chemical Physics Letters 227 (I 994) 327-336
328
from this huge space raises the question whether the method can also be useful for large scale CI calculations. The method is based on extension of the CI space with one state at a time, judging its contribution to the ground state energy against a numerical criterion that is decreased during the process. The addition of new states is combined with a Jacobi-like diagonalization scheme making the actual diagonalization time only a small fraction of the total time involved. The method will be outlined in Section 2.
H~~m+l)=(Hji~)-H~~))Cm+,Sm+,
+Hgq&+,
-&+,)=O.
(5)
This gives the following expression for cm+1and s,+ 1:
(6)
C
with t,+ i, t
2H,Sm) m+’ = Hhm’ -Hii”’ +J($“‘) _H;i”‘)2+4H$“)2
’
(7) A plane rotation decreases the value of the diagonal element Hii by an amount A, while increasing the value of Hii by the same amount,
2. The importancesampling algorithm We will start by introducing some notations. The full CI matrix (dimension N) in a given basis is denoted H, with matrix elements HP The eigenvalues Ei of H are assumed to be ordered such that E, G Ez Q ... Q EN. These eigenvalues are obtained by transforming the matrix H by the unitary matrix U. Obviously the diagonalization problem is to find this matrix. In the classical Jacobi scheme [8] U is obtained as the product of m plane rotation matrices u(k)
U = lim U(m) ,
(1)
m-co
Uw),U(lqJ(2)***p)
(2)
, ik
uck’=uck)(ik
9
ck
,.. Sk 1 i
-Sk
... ck
jk ck Sk) = >
jk
3
1 (3) Each rotation annihilates two off-diagonal matrix elements. To find the values of cm+, and sm+, for rotation m+ 1 we will look at the intermediary matrix H(m)7 H(m)=U(m,=HU(m).
(4)
The off-diagonal element H&m’ is annihilated by solving the 2 x 2 equation
Htnt+l) ,Hlim’ -A,+, l, Him+l) cH$~’
,
+&+l,
(8)
A m+l =SZ,+~(H~im’-H~m’)+2~~+~s,+~H~mm) =t,+,H,S”‘).
(9)
The method is monotonically convergent because with each rotation the sum of the squares of diagonal elements increases and, since C&i Hhrnj2 is constant, the sum of the squares of the off-diagonal elements converges to zero. The speed of convergence is strongly dependent on the order of the rotations because with each rotation older off-diagonal zeroes may be destroyed. For good convergence it is thus expedient to annihilate the off-diagonal elements (roughly) in order of magnitude. The classical Jacobi method aims as the complete diagonalization of a matrix. In most CI calculations one is not interested in finding all eigenvalues of the matrix, but rather in finding the lowest eigenvalue. In this case we need only the first row of the transformation matrix U. This implies that we can keep i, = 1 for all m, reducing the search for large off-diagonal elements to elements in the first column of H only. Another important point in these types of calculations is the size of the off-diagonal elements. Usually these are small compared to the diagonal elements making the number of rotations necessary of order N rather than N2. These two considerations lead to the development of the stochastic diagonalization method, in this con-
L. V&her et al. 1 Chemical PhysicsLetters227 (1994) 327-336
text called the StochCI method. We combine the calculation of the first row of U with a state selection procedure. We will work in a sub block of the complete matrix that includes only those states that give rise to large off-diagonal elements with the (transformed) Hi, matrix element. When these large offdiagonal elements are annihilated the threshold for inclusion of new states is relaxed and the extended block is transformed until no more large off-diagonal elements remain. This scheme keeps the working dimension (n) of the matrices small, allowing (in core) storage of the matrix elements. The final set of states usually gives a fairly good sampling of the important states and may be used as the starting point of larger MRCI calculations or can be used directly if a compact representation of the ground state wavefunction is desired. The process is started by choosing one or more (nmf) reference determinants. We select the determinant with the lowest diagonal energy and assign it position 1 in the matrix. We then perform plane rotations with the constraint that & be smaller than nrrf and that Ak be larger than a given threshold TA. If no more rotations can be found that satisfy these conditions we generate ntrialnew determinants and check whether these give rise to rotations that do give a decrease that is larger than T,. These determinants are added one by one to the CI space, checking after each addition whether an additional rotation is necessary to keep the condition Ap T, for all states j< n. The process is repeated until no more determinants can be found that give a decrease Ak> TA. We then decrease T, by a factor TD and repeat the steps given above. The process is terminated if the number of states n reaches its maximum value N, or if TA drops below a specified value. An important point in the procedure is the generation of new trial states, since the Hamiltonian contains only one- and two-body operators we can restrict the possible trial states to single and double excitations from the determinants that span the current subspace. In general this gives such a large set of possible states that it is not feasible to calculate all Ak. We will therefore further restrict the full set of singly and doubly excited states to a smaller set. We have tested two schemes to do this. In the ‘random’ method new trial states were generated by selecting in a random manner one of the
329
determinants in the present space and exciting from one or two randomly determined occupied orbital( s) to one or two randomly chosen unoccupied orbital(s). The frequency ratio between single and double excitations was taken to be 3 : 7. After generating 100 new determinants the determinant which gave the largest value of A was chosen and added to the CI space. This process is similar to the one used in the original formulation of the stochastic diagonalization method for the two-dimensional Hubbard model [7]. In the second ‘approximation’ method we generate new states by taking all single and double excitations from statesj that satisfied the criterion TMj> TA, with TMj=Ujf, max(A,J . Thus, the maximum value of A arising from previous off-diagonal elements is used as a first estimate for the importance of states arising from reference state j. One by one the generated determinants are added to the CI space if they satisfy the criterion A> T,,
.(,. J-)-l.
(11)
To calculate A(,+, ) we need the diagonal elements HI?) and Hh’;‘,,n+, and the off-diagonal matrix element Hi’;‘,,, . Of course Hi?) is known since this is the element that we are minimizing. The new diagonal element HA’;‘,,,,, is the expectation value of a single determinant since the previous rotations did not affect this matrix element. The element HA’${,L, however, depends on the previous rotations and is given by
k=l
(12)
The calculation of this matrix element presents the most time-consuming step in the process of selecting new states. Because the ground state will be at any stage a linear combination of all determinants included in the set, an exact calculation implies the calculation of matrix elements of the new state with all previous states. To reduce this time we use a second intermediary screening before calculating the com-
L. Visscher et al. /Chemical Physics Letters 227 (1994) 327-336
330
plete matrix element, by restricting the sum over k to states with a large coefficient U,$I. A$?‘, I=
?
k:U&’
>C,
H n+dW
.
(13)
We use the approximation Ah’;‘,,, instead of exact element Him+\,, and calculate an approximate value A Only states that fulfil the conof &z+l):~(m+l). straint A$,,+,) > T, are considered for the extension of the CI space that involves the calculation of the full matrix element. For each accepted state Hi’;‘,,, is tinally calculated without approximations. Both the ‘random’ and the ‘approximation’ variants involve the storage of non-zero matrix elements Hu with i and j< n, the storage of the transformed matrix elements HI'J) , and one column of the transformation matrix U,@). For the efficient calculation of the CI matrix elements we also kept a bit representation of the selected determinants and on the MO integrals in core memory. Using 120 Mb of core memory it was possible to run calculations with up to 60000 determinants on an IBM RS/6000-320H workstation.
in the StochCI method the Jacobi procedure is stopped if no rotations that decrease the HI 1 by more than the final threshold are possible. This choice can make the diagonalization time much smaller at the expense of not giving the exact solution in the given expansion space. This is more consistent with the truncation criterion, however. Because we add the states one by one we can in principle account for the influence of the m previous rotations in the calculation of A(,+,). This matrix element is calculated exactly while this is done perturbatively in other selected MR-CI schemes. The use of the approximation Atm+l, makes this difference less prominent, however. It is also instructive to look at the correspondence of the Jacobi procedure with perturbation theory. If we look at the energy of the ground state we find that each rotation Utk) gives a contribution Ak that can be written as Ak=2HI:-1)2(Hlk-1)-HHI:-‘))-1 II
X(1+
JZ)-‘.
(14)
3. Comparison with other methods The method can obviously be classified as a selected multi-reference CI method. Of the different methods presented in the literature it bears most close resemblance to the CI + PT method of Harrison [ 6 1. In this method one also extends the initial reference space with states that lower the energy by more than a numerical criterion. This criterion is then decreased until the energy converges. Apart from the fact that the StochCI method does not employ perturbation theory to incorporate the effect of the remaining states, the differences are connected to the selection of new determinants. The random version of the algorithm is obviously different. The approximation version resembles the algorithm described by Harrison more closely. Contrary to that method, it still scans only part of the interacting space due to the restriction in Eq. ( 10) but it does this in a deterministic way. A major difference concerns the diagonalization of the final expansion space. In the MRD-CI, CIPSI and CI + PT methods the matrix is made diagonal, while
The total correlation energy is IF=, Ak which is equivalent with the second-order perturbation energy expression if the following approximations are made. First: each state is involved in only one plane rotation (m= n) and all single and double excitations from the reference are included (no screening). Second: we take Hck)nl 1 = H, 1 for all k. Third: we approximate the square root of formula ( 14) by one. We then obtain formula ( 15 ) which is equivalent to the perturbation energy expression
(15) The connection with perturbation theory may be further elaborated by using expression ( 14) to include contributions from states that are not selected. We have not tested this method but examples of similar approaches based on the above mentioned and other MRCI selection schemes can be found in the literature [ 4-691.
L. Visscheret al. /Chemical PhysicsLetters2.27 (1994) 327-336
4. Results
We have tested the method for a number of small systems for which (near) full CI results are available or easily obtainable. 4. I. Magnesium atom The first test problem was the magnesium atom. We used the same basis as was given in Ref. [ 11, with the same frozen 1s function. The initial value of T, was taken to be 0.0 1 hartree, and decreased with a TD factor of 5 during the process. The process was started with the closed shell reference determinant and two different state generation schemes (‘random’ or ‘approximation’) of extending the reference space were employed. The run was ended when the specified number of determinants (2000 or 40000) was reached (Table 1). The results for the magnesium atom show that both the random version and the approximation version give lower energies with fewer determinants than the conventional singles-doubles-triples-etc. excitation method. The one by one addition scheme used in the StochCI method leads to symmetry breaking of the wavefunction. This occurs because not all determinants that are related by a symmetry operation of the full rotational group are added simultaneously. In the calculation of (x2), (y’) and ( z2) this gives differences of 0.02% for the different directions. The problem may of course be easily solved by working
Table 1 Correlation
energies of the Mg atom calculated Method
SCF RASSCF CISD CISDT CISDTQ STOCHCI STOCHCI STOCHCI STOCHCI full CI
[2 ]
(random) (random) ( appr. ) (appr. )
with different
331
with fully symmetry adapted functions or by adding all related determinants simultaneously with the same coefficient. 4.2. Be and Be2 As a second test system we chose the Be2 problem. It is well known that this weakly bound system is difficult to describe with a single reference method that includes only single and double type of excitations. The nearly degenerate ls22s2 and ls22p2 conligurations on the Be atoms make the inclusion of triples or even quadruples necessary. Because our method should sample these important higher order corrections it is therefore a good test case. We used a rather small [ 4s2p 1d] AN0 basis [ lo] that does allow for comparison with full CI results. Be2 is not bound in this small basis but the potential curve shows a local minimum at approximately the right distance. The 1s orbitals were kept frozen. The initial value of T, was 0.01 hartree, the value of T, was 2 or 5. The StochCI process was stopped when the number of 2000 (in the Tn=2 case) or 5000 (in the T,=5 case) determinants was reached. At that stage the values of TA were 0.08 uhartree and 0.005 uhartree respectively. The potential curves calculated with different methods are shown in Fig. 1. We see that most methods do not reproduce the shallow minimum. The single-reference with doubles methods CLSD and CC-SD give a purely repulsive curve. The CI-SDT method gives a minimum but much too deep compared to the full CI result. The
methods
Energy (hartree)
Corr. energy
- 199.585214 -199.615701 - 199.721386 - 199.722039 - 199.726256 - 199.722921 - 199.726113 - 199.724237 - 199.726164
0.0000 -0.0305 -0.1362 -0.1368 -0.1410 -0.1377 -0.1390 -0.1409 -0.1410
- 199.7263
-0.1411
(hartree )
No. of determinants (in Dzh symmetry) 1 4 2960 102928 1964232 2000 40000 2000 40000 2538603250
L. Visscheret al. / ChemicalPhysicstitters 227 (1994) 327-336
332
-
CI-SD (575)
-
CI-SDT
(5159)
culation. The standard deviation of the differences of the StochCI (2000) energies with the full CI energies is 0.24 mhartree, while this is 0.10 mhartree in the MRCI calculation. In the 5000 determinant calculation the differences with the full CI fall below 0.13 mhartree with a standard deviation of only 0.04 mhartree. In the analysis of the wavefunction single, triple and quadruple excitations were found to be equally contributing, confirming the importance of the higher excitations.
A-
4.3. Hydrogen jluoride
-29.236
Hydrogen fluoride represents one of the simplest molecules that possesses a significant permanent dipole. The calculation of the dipole moment is a probe for the ability of the method to sample the important determinants for calculation of properties other than the energy. We used a reasonably large [ 6s5p3d2f], [ 6s4p3d] basis [ IO]. In such a basis full CI results are hard to obtain. We therefore compare our results with the large MRCI calculations of Bauschlicher et al. [ 13 ] and with experiment [ 14 1. The starting value of TA was taken as 0.01 hartree with T,=5 for the random algorithm and T,, = 2 for the approximation scheme. The results are given in Table 2. The random sampling algorithm (selecting the most important of 200 randomly generated determinants) gives an irregular curve in the vicinity of the minimum. The exact coincidence of the lowest minimum with the experimental value of r. is therefore not significant. The approximation version gives a lower total energy and a regular curve. In the analysis
-29.238 4.0
5.0
6.0
7.0
8.0
9.0
R (Bohr)
Fig. 1. Potential energy curves of the ‘Z: ground state of Bea, calculated with different methods. When appropriate the numbers in parentheses give the number of determinants in the wavefunction. The 1s electrons are left uncorrelated.
non-iterative inclusion [ 111 of triples in the CCSD(T) method makes the curve better but still slightly repulsive. This result is consistent with the observation of Sosa et al. [ 12 ] that the full inclusion of triples is necessary for this case. The StochCI calculation can be compared with a MRCI calculation where all single and double excitations from a CAS reference space (four electrons in the eight Be 2s and 2p orbitals) were taken. StochCI gives faster energy convergence, but the shape of the potential curve is better represented in the MRCI calculation than in the 2000 determinant StochCI calTable 2 Spectroscopic constants for the ‘Z ground state of hydrogen fluoride Method
No. of determinants
Eo
0, (eV)
r0 (A)
P (D)
SCF CISD STOCHCI (random) STOCHCI (appr. ) MRCI [ 131
1 34927 5000 5000
-
4.35 7.24 5.94 6.08 6.07
0.900 0.910 0.917” 0.922 0.919
1.922 1.814 1.787 1.802
0.917
1.826
experiment [ 141
100.07077 100.35174 100.30480 100.31062
(4.35) (5.80) (5.14) (5.25) (6.07)
6.123
The figures in parentheses give the dissociation energy corrected for the basis set superposition error using the counterpoise method t151. ’ No single minimum found due to irregularities in the curve. Value given corresponds to the lowest minimum.
L. Visscheret al. / ChemicalPhysicsLetters227 (1994) 327-336
of the latter wavefunction the 5000 sampled determinants were all found to be double excitations from the reference determinant, indicating that a single reference approach would be appropriate in this case. We estimated the basis set superposition error (BSSE) by calculating the energy of the monomers in the dimer basis at the experimental equilibrium distance. This error was found to be of the order of 1 electron volt, slightly smaller than the value for the CISD calculation ( 1.44 eV), but it clearly shows that the good agreement of the uncorrected De value with the experimental result is fortuitous. The large value of the correction arises in the last stage of the calculation of the fluor monomer energy. The calculations with and without the hydrogen ghost basis were both stopped after 5000 determinants were selected. In the calculation without the ghost basis 3438 singly and doubly excited determinants are selected, after which the weakly contributing quadruple corrections come in. In the calculation with the ghost basis only singly and doubly excited states come in. This indicates that the basis set is not fully saturated at this level of theory. The dipole moment is calculated using the calculated densities at the experimental geometry. The deviation from experiment is about twice as large as Table 3 Energies of the Hz0 molecule Method
SCF FCI [ 161 CI ]61
StochCI
Dimension
(hartree)
Threshold
with different
in the CI-SD calculation. From these results it is clear that in this case the simple CLSD scheme works better than the StochCI selection schemes. 4.4. H,O The water molecule is one of the most studied examples in chemistry. We compare our results with the benchmark full CI result of Bauschlicher [ 16 ] and with the selected CI scheme of Harrison [ 61. The basis set employed is of DZP quality and is given in Ref. [ 161. The 1s orbital is frozen after the SCF step, the other orbitals are obtained from a CASSCF calculation that has an active space of eight electrons in the 3a,4a, lb12b, lb22b2 orbitals. To compare with the results of Harrison we started with a T, value of 0.0128 hartreeandusedaTpvalueof2 (Table3). Although the thresholds used are the same, the results of both methods are still difficult to compare since different expansions are used, namely CSFs in Harrison’s method and determinants in StochCI. At first sight the results indicate that the method of Harrison converges far more rapidly. This is, however largely due to the fact that the contribution of a single determinant to the correlation energy will be smaller
methods
dim.
energy
28233466 1029 1512 2185 3210 5314 9219 15023 24337 764 1001 1269 1662 2050 2830 4106 5991
- 76.040542 - 76.256624 - 76.244338 - 76.247649 -76.249863 -76.251418 - 76.252759 -76.253829 - 76.254593 - 76.255215 - 76.234613 - 76.239674 - 76.243301 -76.246518 -76.248108 - 76.249624 -76.250907 -76.251886
given in determinants
and thresholds ISR,
R.
2.5 x lo-’ 1.25x 10-r 6.25x IO-’ 3.125~10-~ 1.5625x lO-6 7.8125x lo-’ 3.90625 x IO-’ 1.953125x lo-’ 2.5x lO-5 1.25x 10-s 6.25x lO-6 3.125x lO-6 1.5625x lO-6 7.8125~ lo-’ 3.90625x lO-7 1.953125x lo-’ of the expansion
calculated
333
2R.
dim.
energy
1187 1720 2587 3901 5825 8901 13280 20242 630 948 1479 2212 3126 4278 5852 7893
- 75.800494 -76.071405 - 76.058030 -76.062216 - 76.065023 - 76.066929 -76.068219 -76.069191 -76.069835 - 76.070352 - 76.044469 -76.051102 - 76.056805 - 76.060861 - 76.063487 - 76.065279 - 76.066542 - 76.067428
for the StochCI and the full CI calculations
dim.
energy
1090 1612 2297 3416 5116 7475 10711 16014 725 1145 1600 2290 3124 4115 5386 7166
- 75.582286 - 75.952269 -75.941241 -75.945007 -75.947190 - 75.948760 -75.949953 - 75.950692 -75.951198 - 75.95 1577 - 75.926657 - 75.934895 - 75.939507 -75.943138 - 75.945442 -75.946830 -75.947819 -75.948755
and in CSFs for the selected CI treatment.
334
L. Visscher et al. /Chemical Physics Letters 227 (1994) 327-336
than the contribution of a CSF which means that a smaller space will be selected if one uses the same criterion. At the threshold of 1.953125~ lo-’ the number of determinants selected by the StochCI scheme is 5991, 7893 and 7166 for R,= 1, 1.5 and 2 respectively, while the number of CSFs selected is 24337, 20242 and 16014. Assuming that the ratio between the number of determinants and CSFs in this expansion is the same (4.25) as in the FCI expansion [ 161 we may compare the StochCI energies with the energies found with a CSF threshold of 1.25 x 1OF5 or 6.25~ 10e6. In that case the StochCI energies are lower than the energies found with the other method. This implies that the StochCI calculation converges faster in terms of expansion space but needs to be converged to a higher value of T, to obtain the same precision. 4.5. Cr, The bonding in the chromium dimer is still one of the unsolved problems in theoretical chemistry. The delicate balance between the magnetic coupling of the two high-spin singlet atoms and the possible five-fold bond is hard to describe with a small number of configurations. We have tried a 5000 determinant description using the natural orbitals of a much larger CASSCF ( 1072 16 determinants) calculation. The purpose of the calculation was to see which determinants are important for further use in, for instance, MRCI calculations. We used the segmented [8s6p4dlf] basis set that was used by Walch et al. [ 17 ] in their CASSCF calculations on Cr2. T, (initial) was 0.01 hat-tree, with a T,, value of 2. The process was either stopped at T,= 10 phartree or when 5000 determinants were selected. The results are given in Fig. 2. Perhaps what is surprising is that the energies of the 5000 determinant calculation are lower than the reference calculation. This can be explained by the importance of excitations to the low-lying 4p and 3d’ orbitals, that were not included in the CASSCF space consisting of the 3d and 4s orbitals only. The StochCI selection scheme picks up these contributions resulting in a better description at shorter distances. The lowering of the energy gives rise to a minimum in the potential curve at about 3.1 bohr. This is shorter than the experimental bond distance of 3.17 bohr [ 181 which is probably
-0SMJ
-0.620
I I
I I
2.8
3.0
I
I I
3.2
3.4
/ 2.6
3.6
R (Bohr)
Fig. 2. Potential energy curves of the ‘Z: ground state of Cr,, calculated with different methods. When appropriate the numbers in parentheses give the number of determinants in the wavefunction. The 3d and 4s electrons were correlated.
due to the error made in the dissociation. Due to the small number of determinants dissociation cannot be described correctly. To describe the separated highspin ‘S atoms with delocalised functions one needs in principle 2x2i2=8192 determinants, which is more than the chosen expansion space. The truncation error is calculated by performing a calculation at a bond distance of 10 bohr using the same selection criteria as at the shorter distances. This calculation gives an energy that is 39 mhartree above the CASSCF energy at that distance.
5. Discussion The results show that the dynamic selection of states indeed give a (much) shorter expansion than can be obtained using unselected methods. This feature is shared with other selected CI schemes as can be seen from the calculation on the water molecule. A difficult point is the determination of a truncation criterion. In potential energy surface calculations the convergence of energy may differ much at different points, causing errors in the calculated surface. In the systems treated here calculations with a stop value of T, = 1OF5or a fixed number of determinants of 5000 give essentially the same results. This indicates that
L. Visscheret al. /Chemical PhysicsLetters227 (1994) 327-336
10
100
1000
335
10000
Fig. 3. Convergence of A as a function of the number of determinants.
at some points contributions of a large number of weakly contributing determinants are important for a correct description. The value of A as a function of the size of the CI space is plotted for the Mg calculations. Fig. 3 shows that the most important determinants are sampled in the beginning of the iterative process. A block structure is visible in the figure, this reflects the selection process: the bottom line of a block is the threshold T, at that stage of process, while the top line is the previous threshold TD States that fall either below or above these thresholds are not yet selected or were already selected, respectively. After about 1000 states the block structure breaks down. This may be caused by important states that were not found previously because the reference state from which they are formed did not meet criterion ( 10) in the previous step. The alteration of the matrix element in formula ( 12 ) with respect to the previous test may be another source of this effect. The effect of the truncation in the sum of ( 13 ) was found to be small by comparing the results of a run on the Mg atom with the default value of Cn=O.OOl with a calculation where C, was set to zero. After sampling 2000 determinants the energy difference was found to be 23 uhartree, while
with 8500 determinants the error was 12 uhartree. Because the object was to study the energy convergence of different variants of the StochCI method we have written a modular and simple program. A disadvantage of the generality of the modules is that the computational efficiency of the present implementation is not too good. The calculation of the CI matrix elements especially could be much improved upon. This limits the significance of mentioning CPU times for the different runs. To get an impression one may use the figures for the run on the Be2 molecule. The 2000 determinant calculation took about 0.7 h on a IBM RS/6000-320H workstation, while the 5000 determinant calculation took 2.7 h. Both CPU times are dominated by evaluating formulae 12 and 13, accounting for 70%-80% of the time.
6. Conclusions Several implementations of the stochastic diagonalization method have been tested on the CI problem. Random sampling of the determinantal space was found to be less efficient than selection methods that are based on the weight of the reference contig-
336
L. V&her et al. / ChemicalPhysicsLetters227 (I 994) 327-336
uration. Alteration of the scheme along these lines makes the StochCI method more resemblant to other selected CI schemes like CIPSI, MRD-CI and CI + PT, although the distinct features of adding one state at a time using an exact criterion and the incomplete diagonalisation of the final space remain. The efficiency of the method can be increased significantly by employing approximate selection criteria. The maximum number of states that can be used in the present formulation and implementation is about 50000 which is enough to obtain a good sampling of important states but in most cases too small to approach the full CI limit to the desired (chemical) accuracy. The implementation is, however, suitable for further optimisations, in particular the calculation of the test criterion A$,,+, ) can be well parallelized. The automatic selection of important states presented here can be used in cases where a good a priori choice of reference states is difficult and may be extended with a subsequent perturbative or MRCI like calculation. Calculations on Mg, Be*, HF and HZ0 show that the method is especially suited for cases where the zeroth-order ground state cannot be described with one determinant. The results for the Cr2 dimer indicate that short expansions of the wavefunction may be able to describe the bonding in this molecule provided that the selection of determinants is done carefully.
Acknowledgement
LV was supported by Cray Research Grant 93.11. A grant of computing time from the Netherlands National Computing Facilities Foundation (NCF) is gratefully acknowledged. One of the referees is thanked for a number of constructive remarks and suggesting the calculation on the Hz0 molecule reported in Section 4.
References [I ] C.W. Bauschlicher Jr., S.R. Langhoff and P.R. Taylor, Advan. Chem. Phys. 77 (1990) 103; P.E.M. Siegbahn, Lecture Notes Chem. 58 (1992) 255. [2] J. Olsen, P. Jorgensen and J. Simons, Chem. Phys. Letters 169 (1990) 463. [3] E.R. Davidson and D.W. Silver, Chem. Phys. Letters 52 (1977) 403. [4] B. Huron, J.P. Malrieu and P. Rancurel, J. Chem. Phys. 58 (1973) 5745. [5] R.J. Buenker and SD. Peyerimhoff, Theoret. Chim. Acta 39 (1975) 217; R.J. Buenker, S.D. Peyerimhoff and W. Butscher, Mol. Phys. 35 (1978) 771. [6] R.J. Harrison, J. Chem. Phys. 94 (1991) 5021. [ 71 H. DeRaedt and W. von der Linden, Intern. J. Mod. Phys. C 3 (1992) 97; Phys. Rev. B 45 (1992) 8787; H. DeRaedt and M. Frick, Phys. Rept. 23 1 ( 1993) 107. [ 81 J.H. Wilkinson, The algebraic eigenvalue problem (Clarendon Press, Oxford, 1965 ). [9] P.J.A. Ruttink and M.M.M. van Schaik, J. Comput. Chem. 6 (1985) 88; G.A. Segal, R.W. Weunore and K. Wolf, Chem. Phys. 30 (1978) 269. [lo] P.-O. Widmark, P.-A. Malmqvist and B.O. Roos, Theoret. Chim. Acta 77 ( 1990) 29 1. [ 111 J.A. Pople, M. Head-Gordon and K Raghavachari, J. Chem. Phys. 87 (1987) 5968. [ 121 C. Sosa, J. Noga and R.J. Bartlett, J. Chem. Phys. 88 ( 1988) 5974. [ 131 C.W. Bauschlicher Jr., S.P. Walch, S.R. Langhoff, P.R. Taylor and R.L. Jaffe, J. Chem. Phys. 88 ( 1988) 1743. [ 141 K.P. Huber and G. Herzberg, Constants of diatomic molecules (Van Nostrand Reinhold, New York, 1979); R.C. Weast, ed., Handbook of chemistry and physics, 72nd Ed. (CRC Press, Boca Raton, 1992). [ 151 S.F. Boys and F. Bemardi, Mol. Phys. 19 ( 1970) 553. [ 161 C.W. Bauschlicher Jr. and P.R. Taylor, J. Chem. Phys. 85 (1986) 2779. [ 171 S.P. Walch, C.W. Bauschlicher Jr., B.O. Roos and C.J. Nelin, Chem. Phys. Letters 103 (1983) 175. [ 181 M.D. Morse, Chem. Rev. 86 (1986) 1049.
ELSEVIER
CHEMICAL PHYSICS LETTERS
Chemical Physics Letters 227 (1994) 327-336
A new configuration selection method for configuration interaction calculations L. Visscher a,’, H. DeRaedt b, W.C. Nieuwpoort a aDepartment of Chemical Physics and Material Science Centre, University of Groningen. Nijenborgh 4, 9747 AG Groningen, The Netherlands b Institutefor Theoretical Physics and Material Science Centre, University of Groningen, Nijenborgh 4, 9747AG Groningen, The Netherlands
Received 11 March 1994; in final form 21 June 1994
Abstract The recently proposed stochastic diagonalization method is applied to the ab initio quantum chemistry configuration interaction problem. In this context it can be viewed as a multi-reference CI method with dynamic selection of important configurations. The method is compared with other methods and tested by calculations on a number of small molecular systems for which accurate results are available. A calculation on the Cr, dimer is presented to show the capability of the algorithm to find short expansions of molecular wavefunctions.
1. Introduction
The optimal choice of (reference) configurations in truncated configuration interaction methods has been a subject of discussion over the years [ 11. The well-known configuration interaction singles doubles (CISD) method uses an extension of a given reference space that is systematic and unambiguous. The resulting algorithm can be well vectorized and is used to treat matrices of large dimensions [2]. A disadvantage is the lack of size-extensivity that can only partly be corrected by applying corrections [ 3 ] afterwards. Another disadvantage is the scaling of the method with the number of basis functions involved. Full-C1 (FCI) calculations indicate [ 2 ] that most of the configurational space does not contribute significantly to the ground state wavefunction, yet many ’ Present address: NASA Ames Research Center, MS 230/3, Moffet Field, CA 94035-l 000, USA.
configurations do arise in the standard excitation schemes. The sparsity of the CI matrix is used by methods like CIPSI [ 4 1, MRDCI [ 5 ] and the CI + PT method of Harrison [ 61 that select only the important configurations. Experience has already shown that such methods that combine perturbation theory and variational methods can give reliable results while using only a fraction of the full CI space. The aim of this Letter is therefore not to present a new efficient implementation of such a selected CI method but rather to test a method that originates from a different application field. The novel algorithm that we tested is based on the stochastic diagonalization method [ 71 that was recently developed for general quantum many-body problems. The method was used in applications on the two-dimensional Hubbard model where the size of the Hilbert space may become as large as 1035.The success of the algorithm to sample important states
0009-2614/94/$07.00 0 1994 Elsevier Science B.V. All rights reserved SSDIOOO9-2614(94)00824-S
L. Visscher et al. / Chemical Physics Letters 227 (I 994) 327-336
328
from this huge space raises the question whether the method can also be useful for large scale CI calculations. The method is based on extension of the CI space with one state at a time, judging its contribution to the ground state energy against a numerical criterion that is decreased during the process. The addition of new states is combined with a Jacobi-like diagonalization scheme making the actual diagonalization time only a small fraction of the total time involved. The method will be outlined in Section 2.
H~~m+l)=(Hji~)-H~~))Cm+,Sm+,
+Hgq&+,
-&+,)=O.
(5)
This gives the following expression for cm+1and s,+ 1:
(6)
C
with t,+ i, t
2H,Sm) m+’ = Hhm’ -Hii”’ +J($“‘) _H;i”‘)2+4H$“)2
’
(7) A plane rotation decreases the value of the diagonal element Hii by an amount A, while increasing the value of Hii by the same amount,
2. The importancesampling algorithm We will start by introducing some notations. The full CI matrix (dimension N) in a given basis is denoted H, with matrix elements HP The eigenvalues Ei of H are assumed to be ordered such that E, G Ez Q ... Q EN. These eigenvalues are obtained by transforming the matrix H by the unitary matrix U. Obviously the diagonalization problem is to find this matrix. In the classical Jacobi scheme [8] U is obtained as the product of m plane rotation matrices u(k)
U = lim U(m) ,
(1)
m-co
Uw),U(lqJ(2)***p)
(2)
, ik
uck’=uck)(ik
9
ck
,.. Sk 1 i
-Sk
... ck
jk ck Sk) = >
jk
3
1 (3) Each rotation annihilates two off-diagonal matrix elements. To find the values of cm+, and sm+, for rotation m+ 1 we will look at the intermediary matrix H(m)7 H(m)=U(m,=HU(m).
(4)
The off-diagonal element H&m’ is annihilated by solving the 2 x 2 equation
Htnt+l) ,Hlim’ -A,+, l, Him+l) cH$~’
,
+&+l,
(8)
A m+l =SZ,+~(H~im’-H~m’)+2~~+~s,+~H~mm) =t,+,H,S”‘).
(9)
The method is monotonically convergent because with each rotation the sum of the squares of diagonal elements increases and, since C&i Hhrnj2 is constant, the sum of the squares of the off-diagonal elements converges to zero. The speed of convergence is strongly dependent on the order of the rotations because with each rotation older off-diagonal zeroes may be destroyed. For good convergence it is thus expedient to annihilate the off-diagonal elements (roughly) in order of magnitude. The classical Jacobi method aims as the complete diagonalization of a matrix. In most CI calculations one is not interested in finding all eigenvalues of the matrix, but rather in finding the lowest eigenvalue. In this case we need only the first row of the transformation matrix U. This implies that we can keep i, = 1 for all m, reducing the search for large off-diagonal elements to elements in the first column of H only. Another important point in these types of calculations is the size of the off-diagonal elements. Usually these are small compared to the diagonal elements making the number of rotations necessary of order N rather than N2. These two considerations lead to the development of the stochastic diagonalization method, in this con-
L. V&her et al. 1 Chemical PhysicsLetters227 (1994) 327-336
text called the StochCI method. We combine the calculation of the first row of U with a state selection procedure. We will work in a sub block of the complete matrix that includes only those states that give rise to large off-diagonal elements with the (transformed) Hi, matrix element. When these large offdiagonal elements are annihilated the threshold for inclusion of new states is relaxed and the extended block is transformed until no more large off-diagonal elements remain. This scheme keeps the working dimension (n) of the matrices small, allowing (in core) storage of the matrix elements. The final set of states usually gives a fairly good sampling of the important states and may be used as the starting point of larger MRCI calculations or can be used directly if a compact representation of the ground state wavefunction is desired. The process is started by choosing one or more (nmf) reference determinants. We select the determinant with the lowest diagonal energy and assign it position 1 in the matrix. We then perform plane rotations with the constraint that & be smaller than nrrf and that Ak be larger than a given threshold TA. If no more rotations can be found that satisfy these conditions we generate ntrialnew determinants and check whether these give rise to rotations that do give a decrease that is larger than T,. These determinants are added one by one to the CI space, checking after each addition whether an additional rotation is necessary to keep the condition Ap T, for all states j< n. The process is repeated until no more determinants can be found that give a decrease Ak> TA. We then decrease T, by a factor TD and repeat the steps given above. The process is terminated if the number of states n reaches its maximum value N, or if TA drops below a specified value. An important point in the procedure is the generation of new trial states, since the Hamiltonian contains only one- and two-body operators we can restrict the possible trial states to single and double excitations from the determinants that span the current subspace. In general this gives such a large set of possible states that it is not feasible to calculate all Ak. We will therefore further restrict the full set of singly and doubly excited states to a smaller set. We have tested two schemes to do this. In the ‘random’ method new trial states were generated by selecting in a random manner one of the
329
determinants in the present space and exciting from one or two randomly determined occupied orbital( s) to one or two randomly chosen unoccupied orbital(s). The frequency ratio between single and double excitations was taken to be 3 : 7. After generating 100 new determinants the determinant which gave the largest value of A was chosen and added to the CI space. This process is similar to the one used in the original formulation of the stochastic diagonalization method for the two-dimensional Hubbard model [7]. In the second ‘approximation’ method we generate new states by taking all single and double excitations from statesj that satisfied the criterion TMj> TA, with TMj=Ujf, max(A,J . Thus, the maximum value of A arising from previous off-diagonal elements is used as a first estimate for the importance of states arising from reference state j. One by one the generated determinants are added to the CI space if they satisfy the criterion A> T,,
.(,. J-)-l.
(11)
To calculate A(,+, ) we need the diagonal elements HI?) and Hh’;‘,,n+, and the off-diagonal matrix element Hi’;‘,,, . Of course Hi?) is known since this is the element that we are minimizing. The new diagonal element HA’;‘,,,,, is the expectation value of a single determinant since the previous rotations did not affect this matrix element. The element HA’${,L, however, depends on the previous rotations and is given by
k=l
(12)
The calculation of this matrix element presents the most time-consuming step in the process of selecting new states. Because the ground state will be at any stage a linear combination of all determinants included in the set, an exact calculation implies the calculation of matrix elements of the new state with all previous states. To reduce this time we use a second intermediary screening before calculating the com-
L. Visscher et al. /Chemical Physics Letters 227 (1994) 327-336
330
plete matrix element, by restricting the sum over k to states with a large coefficient U,$I. A$?‘, I=
?
k:U&’
>C,
H n+dW
.
(13)
We use the approximation Ah’;‘,,, instead of exact element Him+\,, and calculate an approximate value A Only states that fulfil the conof &z+l):~(m+l). straint A$,,+,) > T, are considered for the extension of the CI space that involves the calculation of the full matrix element. For each accepted state Hi’;‘,,, is tinally calculated without approximations. Both the ‘random’ and the ‘approximation’ variants involve the storage of non-zero matrix elements Hu with i and j< n, the storage of the transformed matrix elements HI'J) , and one column of the transformation matrix U,@). For the efficient calculation of the CI matrix elements we also kept a bit representation of the selected determinants and on the MO integrals in core memory. Using 120 Mb of core memory it was possible to run calculations with up to 60000 determinants on an IBM RS/6000-320H workstation.
in the StochCI method the Jacobi procedure is stopped if no rotations that decrease the HI 1 by more than the final threshold are possible. This choice can make the diagonalization time much smaller at the expense of not giving the exact solution in the given expansion space. This is more consistent with the truncation criterion, however. Because we add the states one by one we can in principle account for the influence of the m previous rotations in the calculation of A(,+,). This matrix element is calculated exactly while this is done perturbatively in other selected MR-CI schemes. The use of the approximation Atm+l, makes this difference less prominent, however. It is also instructive to look at the correspondence of the Jacobi procedure with perturbation theory. If we look at the energy of the ground state we find that each rotation Utk) gives a contribution Ak that can be written as Ak=2HI:-1)2(Hlk-1)-HHI:-‘))-1 II
X(1+
JZ)-‘.
(14)
3. Comparison with other methods The method can obviously be classified as a selected multi-reference CI method. Of the different methods presented in the literature it bears most close resemblance to the CI + PT method of Harrison [ 6 1. In this method one also extends the initial reference space with states that lower the energy by more than a numerical criterion. This criterion is then decreased until the energy converges. Apart from the fact that the StochCI method does not employ perturbation theory to incorporate the effect of the remaining states, the differences are connected to the selection of new determinants. The random version of the algorithm is obviously different. The approximation version resembles the algorithm described by Harrison more closely. Contrary to that method, it still scans only part of the interacting space due to the restriction in Eq. ( 10) but it does this in a deterministic way. A major difference concerns the diagonalization of the final expansion space. In the MRD-CI, CIPSI and CI + PT methods the matrix is made diagonal, while
The total correlation energy is IF=, Ak which is equivalent with the second-order perturbation energy expression if the following approximations are made. First: each state is involved in only one plane rotation (m= n) and all single and double excitations from the reference are included (no screening). Second: we take Hck)nl 1 = H, 1 for all k. Third: we approximate the square root of formula ( 14) by one. We then obtain formula ( 15 ) which is equivalent to the perturbation energy expression
(15) The connection with perturbation theory may be further elaborated by using expression ( 14) to include contributions from states that are not selected. We have not tested this method but examples of similar approaches based on the above mentioned and other MRCI selection schemes can be found in the literature [ 4-691.
L. Visscheret al. /Chemical PhysicsLetters2.27 (1994) 327-336
4. Results
We have tested the method for a number of small systems for which (near) full CI results are available or easily obtainable. 4. I. Magnesium atom The first test problem was the magnesium atom. We used the same basis as was given in Ref. [ 11, with the same frozen 1s function. The initial value of T, was taken to be 0.0 1 hartree, and decreased with a TD factor of 5 during the process. The process was started with the closed shell reference determinant and two different state generation schemes (‘random’ or ‘approximation’) of extending the reference space were employed. The run was ended when the specified number of determinants (2000 or 40000) was reached (Table 1). The results for the magnesium atom show that both the random version and the approximation version give lower energies with fewer determinants than the conventional singles-doubles-triples-etc. excitation method. The one by one addition scheme used in the StochCI method leads to symmetry breaking of the wavefunction. This occurs because not all determinants that are related by a symmetry operation of the full rotational group are added simultaneously. In the calculation of (x2), (y’) and ( z2) this gives differences of 0.02% for the different directions. The problem may of course be easily solved by working
Table 1 Correlation
energies of the Mg atom calculated Method
SCF RASSCF CISD CISDT CISDTQ STOCHCI STOCHCI STOCHCI STOCHCI full CI
[2 ]
(random) (random) ( appr. ) (appr. )
with different
331
with fully symmetry adapted functions or by adding all related determinants simultaneously with the same coefficient. 4.2. Be and Be2 As a second test system we chose the Be2 problem. It is well known that this weakly bound system is difficult to describe with a single reference method that includes only single and double type of excitations. The nearly degenerate ls22s2 and ls22p2 conligurations on the Be atoms make the inclusion of triples or even quadruples necessary. Because our method should sample these important higher order corrections it is therefore a good test case. We used a rather small [ 4s2p 1d] AN0 basis [ lo] that does allow for comparison with full CI results. Be2 is not bound in this small basis but the potential curve shows a local minimum at approximately the right distance. The 1s orbitals were kept frozen. The initial value of T, was 0.01 hartree, the value of T, was 2 or 5. The StochCI process was stopped when the number of 2000 (in the Tn=2 case) or 5000 (in the T,=5 case) determinants was reached. At that stage the values of TA were 0.08 uhartree and 0.005 uhartree respectively. The potential curves calculated with different methods are shown in Fig. 1. We see that most methods do not reproduce the shallow minimum. The single-reference with doubles methods CLSD and CC-SD give a purely repulsive curve. The CI-SDT method gives a minimum but much too deep compared to the full CI result. The
methods
Energy (hartree)
Corr. energy
- 199.585214 -199.615701 - 199.721386 - 199.722039 - 199.726256 - 199.722921 - 199.726113 - 199.724237 - 199.726164
0.0000 -0.0305 -0.1362 -0.1368 -0.1410 -0.1377 -0.1390 -0.1409 -0.1410
- 199.7263
-0.1411
(hartree )
No. of determinants (in Dzh symmetry) 1 4 2960 102928 1964232 2000 40000 2000 40000 2538603250
L. Visscheret al. / ChemicalPhysicstitters 227 (1994) 327-336
332
-
CI-SD (575)
-
CI-SDT
(5159)
culation. The standard deviation of the differences of the StochCI (2000) energies with the full CI energies is 0.24 mhartree, while this is 0.10 mhartree in the MRCI calculation. In the 5000 determinant calculation the differences with the full CI fall below 0.13 mhartree with a standard deviation of only 0.04 mhartree. In the analysis of the wavefunction single, triple and quadruple excitations were found to be equally contributing, confirming the importance of the higher excitations.
A-
4.3. Hydrogen jluoride
-29.236
Hydrogen fluoride represents one of the simplest molecules that possesses a significant permanent dipole. The calculation of the dipole moment is a probe for the ability of the method to sample the important determinants for calculation of properties other than the energy. We used a reasonably large [ 6s5p3d2f], [ 6s4p3d] basis [ IO]. In such a basis full CI results are hard to obtain. We therefore compare our results with the large MRCI calculations of Bauschlicher et al. [ 13 ] and with experiment [ 14 1. The starting value of TA was taken as 0.01 hartree with T,=5 for the random algorithm and T,, = 2 for the approximation scheme. The results are given in Table 2. The random sampling algorithm (selecting the most important of 200 randomly generated determinants) gives an irregular curve in the vicinity of the minimum. The exact coincidence of the lowest minimum with the experimental value of r. is therefore not significant. The approximation version gives a lower total energy and a regular curve. In the analysis
-29.238 4.0
5.0
6.0
7.0
8.0
9.0
R (Bohr)
Fig. 1. Potential energy curves of the ‘Z: ground state of Bea, calculated with different methods. When appropriate the numbers in parentheses give the number of determinants in the wavefunction. The 1s electrons are left uncorrelated.
non-iterative inclusion [ 111 of triples in the CCSD(T) method makes the curve better but still slightly repulsive. This result is consistent with the observation of Sosa et al. [ 12 ] that the full inclusion of triples is necessary for this case. The StochCI calculation can be compared with a MRCI calculation where all single and double excitations from a CAS reference space (four electrons in the eight Be 2s and 2p orbitals) were taken. StochCI gives faster energy convergence, but the shape of the potential curve is better represented in the MRCI calculation than in the 2000 determinant StochCI calTable 2 Spectroscopic constants for the ‘Z ground state of hydrogen fluoride Method
No. of determinants
Eo
0, (eV)
r0 (A)
P (D)
SCF CISD STOCHCI (random) STOCHCI (appr. ) MRCI [ 131
1 34927 5000 5000
-
4.35 7.24 5.94 6.08 6.07
0.900 0.910 0.917” 0.922 0.919
1.922 1.814 1.787 1.802
0.917
1.826
experiment [ 141
100.07077 100.35174 100.30480 100.31062
(4.35) (5.80) (5.14) (5.25) (6.07)
6.123
The figures in parentheses give the dissociation energy corrected for the basis set superposition error using the counterpoise method t151. ’ No single minimum found due to irregularities in the curve. Value given corresponds to the lowest minimum.
L. Visscheret al. / ChemicalPhysicsLetters227 (1994) 327-336
of the latter wavefunction the 5000 sampled determinants were all found to be double excitations from the reference determinant, indicating that a single reference approach would be appropriate in this case. We estimated the basis set superposition error (BSSE) by calculating the energy of the monomers in the dimer basis at the experimental equilibrium distance. This error was found to be of the order of 1 electron volt, slightly smaller than the value for the CISD calculation ( 1.44 eV), but it clearly shows that the good agreement of the uncorrected De value with the experimental result is fortuitous. The large value of the correction arises in the last stage of the calculation of the fluor monomer energy. The calculations with and without the hydrogen ghost basis were both stopped after 5000 determinants were selected. In the calculation without the ghost basis 3438 singly and doubly excited determinants are selected, after which the weakly contributing quadruple corrections come in. In the calculation with the ghost basis only singly and doubly excited states come in. This indicates that the basis set is not fully saturated at this level of theory. The dipole moment is calculated using the calculated densities at the experimental geometry. The deviation from experiment is about twice as large as Table 3 Energies of the Hz0 molecule Method
SCF FCI [ 161 CI ]61
StochCI
Dimension
(hartree)
Threshold
with different
in the CI-SD calculation. From these results it is clear that in this case the simple CLSD scheme works better than the StochCI selection schemes. 4.4. H,O The water molecule is one of the most studied examples in chemistry. We compare our results with the benchmark full CI result of Bauschlicher [ 16 ] and with the selected CI scheme of Harrison [ 61. The basis set employed is of DZP quality and is given in Ref. [ 161. The 1s orbital is frozen after the SCF step, the other orbitals are obtained from a CASSCF calculation that has an active space of eight electrons in the 3a,4a, lb12b, lb22b2 orbitals. To compare with the results of Harrison we started with a T, value of 0.0128 hartreeandusedaTpvalueof2 (Table3). Although the thresholds used are the same, the results of both methods are still difficult to compare since different expansions are used, namely CSFs in Harrison’s method and determinants in StochCI. At first sight the results indicate that the method of Harrison converges far more rapidly. This is, however largely due to the fact that the contribution of a single determinant to the correlation energy will be smaller
methods
dim.
energy
28233466 1029 1512 2185 3210 5314 9219 15023 24337 764 1001 1269 1662 2050 2830 4106 5991
- 76.040542 - 76.256624 - 76.244338 - 76.247649 -76.249863 -76.251418 - 76.252759 -76.253829 - 76.254593 - 76.255215 - 76.234613 - 76.239674 - 76.243301 -76.246518 -76.248108 - 76.249624 -76.250907 -76.251886
given in determinants
and thresholds ISR,
R.
2.5 x lo-’ 1.25x 10-r 6.25x IO-’ 3.125~10-~ 1.5625x lO-6 7.8125x lo-’ 3.90625 x IO-’ 1.953125x lo-’ 2.5x lO-5 1.25x 10-s 6.25x lO-6 3.125x lO-6 1.5625x lO-6 7.8125~ lo-’ 3.90625x lO-7 1.953125x lo-’ of the expansion
calculated
333
2R.
dim.
energy
1187 1720 2587 3901 5825 8901 13280 20242 630 948 1479 2212 3126 4278 5852 7893
- 75.800494 -76.071405 - 76.058030 -76.062216 - 76.065023 - 76.066929 -76.068219 -76.069191 -76.069835 - 76.070352 - 76.044469 -76.051102 - 76.056805 - 76.060861 - 76.063487 - 76.065279 - 76.066542 - 76.067428
for the StochCI and the full CI calculations
dim.
energy
1090 1612 2297 3416 5116 7475 10711 16014 725 1145 1600 2290 3124 4115 5386 7166
- 75.582286 - 75.952269 -75.941241 -75.945007 -75.947190 - 75.948760 -75.949953 - 75.950692 -75.951198 - 75.95 1577 - 75.926657 - 75.934895 - 75.939507 -75.943138 - 75.945442 -75.946830 -75.947819 -75.948755
and in CSFs for the selected CI treatment.
334
L. Visscher et al. /Chemical Physics Letters 227 (1994) 327-336
than the contribution of a CSF which means that a smaller space will be selected if one uses the same criterion. At the threshold of 1.953125~ lo-’ the number of determinants selected by the StochCI scheme is 5991, 7893 and 7166 for R,= 1, 1.5 and 2 respectively, while the number of CSFs selected is 24337, 20242 and 16014. Assuming that the ratio between the number of determinants and CSFs in this expansion is the same (4.25) as in the FCI expansion [ 161 we may compare the StochCI energies with the energies found with a CSF threshold of 1.25 x 1OF5 or 6.25~ 10e6. In that case the StochCI energies are lower than the energies found with the other method. This implies that the StochCI calculation converges faster in terms of expansion space but needs to be converged to a higher value of T, to obtain the same precision. 4.5. Cr, The bonding in the chromium dimer is still one of the unsolved problems in theoretical chemistry. The delicate balance between the magnetic coupling of the two high-spin singlet atoms and the possible five-fold bond is hard to describe with a small number of configurations. We have tried a 5000 determinant description using the natural orbitals of a much larger CASSCF ( 1072 16 determinants) calculation. The purpose of the calculation was to see which determinants are important for further use in, for instance, MRCI calculations. We used the segmented [8s6p4dlf] basis set that was used by Walch et al. [ 17 ] in their CASSCF calculations on Cr2. T, (initial) was 0.01 hat-tree, with a T,, value of 2. The process was either stopped at T,= 10 phartree or when 5000 determinants were selected. The results are given in Fig. 2. Perhaps what is surprising is that the energies of the 5000 determinant calculation are lower than the reference calculation. This can be explained by the importance of excitations to the low-lying 4p and 3d’ orbitals, that were not included in the CASSCF space consisting of the 3d and 4s orbitals only. The StochCI selection scheme picks up these contributions resulting in a better description at shorter distances. The lowering of the energy gives rise to a minimum in the potential curve at about 3.1 bohr. This is shorter than the experimental bond distance of 3.17 bohr [ 181 which is probably
-0SMJ
-0.620
I I
I I
2.8
3.0
I
I I
3.2
3.4
/ 2.6
3.6
R (Bohr)
Fig. 2. Potential energy curves of the ‘Z: ground state of Cr,, calculated with different methods. When appropriate the numbers in parentheses give the number of determinants in the wavefunction. The 3d and 4s electrons were correlated.
due to the error made in the dissociation. Due to the small number of determinants dissociation cannot be described correctly. To describe the separated highspin ‘S atoms with delocalised functions one needs in principle 2x2i2=8192 determinants, which is more than the chosen expansion space. The truncation error is calculated by performing a calculation at a bond distance of 10 bohr using the same selection criteria as at the shorter distances. This calculation gives an energy that is 39 mhartree above the CASSCF energy at that distance.
5. Discussion The results show that the dynamic selection of states indeed give a (much) shorter expansion than can be obtained using unselected methods. This feature is shared with other selected CI schemes as can be seen from the calculation on the water molecule. A difficult point is the determination of a truncation criterion. In potential energy surface calculations the convergence of energy may differ much at different points, causing errors in the calculated surface. In the systems treated here calculations with a stop value of T, = 1OF5or a fixed number of determinants of 5000 give essentially the same results. This indicates that
L. Visscheret al. /Chemical PhysicsLetters227 (1994) 327-336
10
100
1000
335
10000
Fig. 3. Convergence of A as a function of the number of determinants.
at some points contributions of a large number of weakly contributing determinants are important for a correct description. The value of A as a function of the size of the CI space is plotted for the Mg calculations. Fig. 3 shows that the most important determinants are sampled in the beginning of the iterative process. A block structure is visible in the figure, this reflects the selection process: the bottom line of a block is the threshold T, at that stage of process, while the top line is the previous threshold TD States that fall either below or above these thresholds are not yet selected or were already selected, respectively. After about 1000 states the block structure breaks down. This may be caused by important states that were not found previously because the reference state from which they are formed did not meet criterion ( 10) in the previous step. The alteration of the matrix element in formula ( 12 ) with respect to the previous test may be another source of this effect. The effect of the truncation in the sum of ( 13 ) was found to be small by comparing the results of a run on the Mg atom with the default value of Cn=O.OOl with a calculation where C, was set to zero. After sampling 2000 determinants the energy difference was found to be 23 uhartree, while
with 8500 determinants the error was 12 uhartree. Because the object was to study the energy convergence of different variants of the StochCI method we have written a modular and simple program. A disadvantage of the generality of the modules is that the computational efficiency of the present implementation is not too good. The calculation of the CI matrix elements especially could be much improved upon. This limits the significance of mentioning CPU times for the different runs. To get an impression one may use the figures for the run on the Be2 molecule. The 2000 determinant calculation took about 0.7 h on a IBM RS/6000-320H workstation, while the 5000 determinant calculation took 2.7 h. Both CPU times are dominated by evaluating formulae 12 and 13, accounting for 70%-80% of the time.
6. Conclusions Several implementations of the stochastic diagonalization method have been tested on the CI problem. Random sampling of the determinantal space was found to be less efficient than selection methods that are based on the weight of the reference contig-
336
L. V&her et al. / ChemicalPhysicsLetters227 (I 994) 327-336
uration. Alteration of the scheme along these lines makes the StochCI method more resemblant to other selected CI schemes like CIPSI, MRD-CI and CI + PT, although the distinct features of adding one state at a time using an exact criterion and the incomplete diagonalisation of the final space remain. The efficiency of the method can be increased significantly by employing approximate selection criteria. The maximum number of states that can be used in the present formulation and implementation is about 50000 which is enough to obtain a good sampling of important states but in most cases too small to approach the full CI limit to the desired (chemical) accuracy. The implementation is, however, suitable for further optimisations, in particular the calculation of the test criterion A$,,+, ) can be well parallelized. The automatic selection of important states presented here can be used in cases where a good a priori choice of reference states is difficult and may be extended with a subsequent perturbative or MRCI like calculation. Calculations on Mg, Be*, HF and HZ0 show that the method is especially suited for cases where the zeroth-order ground state cannot be described with one determinant. The results for the Cr2 dimer indicate that short expansions of the wavefunction may be able to describe the bonding in this molecule provided that the selection of determinants is done carefully.
Acknowledgement
LV was supported by Cray Research Grant 93.11. A grant of computing time from the Netherlands National Computing Facilities Foundation (NCF) is gratefully acknowledged. One of the referees is thanked for a number of constructive remarks and suggesting the calculation on the Hz0 molecule reported in Section 4.
References [I ] C.W. Bauschlicher Jr., S.R. Langhoff and P.R. Taylor, Advan. Chem. Phys. 77 (1990) 103; P.E.M. Siegbahn, Lecture Notes Chem. 58 (1992) 255. [2] J. Olsen, P. Jorgensen and J. Simons, Chem. Phys. Letters 169 (1990) 463. [3] E.R. Davidson and D.W. Silver, Chem. Phys. Letters 52 (1977) 403. [4] B. Huron, J.P. Malrieu and P. Rancurel, J. Chem. Phys. 58 (1973) 5745. [5] R.J. Buenker and SD. Peyerimhoff, Theoret. Chim. Acta 39 (1975) 217; R.J. Buenker, S.D. Peyerimhoff and W. Butscher, Mol. Phys. 35 (1978) 771. [6] R.J. Harrison, J. Chem. Phys. 94 (1991) 5021. [ 71 H. DeRaedt and W. von der Linden, Intern. J. Mod. Phys. C 3 (1992) 97; Phys. Rev. B 45 (1992) 8787; H. DeRaedt and M. Frick, Phys. Rept. 23 1 ( 1993) 107. [ 81 J.H. Wilkinson, The algebraic eigenvalue problem (Clarendon Press, Oxford, 1965 ). [9] P.J.A. Ruttink and M.M.M. van Schaik, J. Comput. Chem. 6 (1985) 88; G.A. Segal, R.W. Weunore and K. Wolf, Chem. Phys. 30 (1978) 269. [lo] P.-O. Widmark, P.-A. Malmqvist and B.O. Roos, Theoret. Chim. Acta 77 ( 1990) 29 1. [ 111 J.A. Pople, M. Head-Gordon and K Raghavachari, J. Chem. Phys. 87 (1987) 5968. [ 121 C. Sosa, J. Noga and R.J. Bartlett, J. Chem. Phys. 88 ( 1988) 5974. [ 131 C.W. Bauschlicher Jr., S.P. Walch, S.R. Langhoff, P.R. Taylor and R.L. Jaffe, J. Chem. Phys. 88 ( 1988) 1743. [ 141 K.P. Huber and G. Herzberg, Constants of diatomic molecules (Van Nostrand Reinhold, New York, 1979); R.C. Weast, ed., Handbook of chemistry and physics, 72nd Ed. (CRC Press, Boca Raton, 1992). [ 151 S.F. Boys and F. Bemardi, Mol. Phys. 19 ( 1970) 553. [ 161 C.W. Bauschlicher Jr. and P.R. Taylor, J. Chem. Phys. 85 (1986) 2779. [ 171 S.P. Walch, C.W. Bauschlicher Jr., B.O. Roos and C.J. Nelin, Chem. Phys. Letters 103 (1983) 175. [ 181 M.D. Morse, Chem. Rev. 86 (1986) 1049.