Multi-Configuration SCF Calculations BY R. MCWEENY Dept. of Chemistry The University Sheffield S3 7HF Received 13th September 1968 In multi-configuration self-consistent field (MC SCF) calculations a many-configuration approxi- mation to the wave function is optimized by variation of the orbitals. The equations satisfied by the optimum orbitals are formulated along with a practicable method of numerical solution. Incidental problems arising in the computer implementation of MC SCF calculations are also discussed. 1. INTRODUCTION In the configuration interaction (CI) approach to the construction of electronic wave functions the wave function is expanded in the form Y = CKCKQK (1-1) where the QK are antisymmetric many-electron functions each referring to a con- figuration of occupied orbitals (i.e to an orbital configuration).The orbitals used are drawn from say the first M members of any convenient complete set and the expansion (1. l) in principle infinite is therefore truncated accordingly ; nevertheless even with only a small number of electrons the expansion may be too long to be manageable. Thus with five electrons and ten orbitals we have twenty spin-orbitals and can set up 15,504 spin-orbital configurations if the OK were taken to be single determinants this number would be the dimension of the secular problem to determine the coefficients. It is therefore preferable to take as the mK suitable linear combinations of determinants chosen to conform to symmetry requirements in order to reduce the dimensions of the secular problem. The only type of symmetry that is always available is the spin symmetry if the QK are chosen to be spin eigenfunctions (quantum numbers S M ) we need consider only different orbital configurations (a relatively small number) the different deter- minants within each orbital conguration being combined to yield a small number of spin eigenfunctions of each (S M).When the Hamiltonian is spin free (the usual first approximation) only QK of the same (S M ) occur in the expansion which is greatly reduced in length. For five electrons and ten orbitals we can find 3,300 independent spin eigenfunctions with S = A4 = 3. In a complete CI calculation with the given one-electron basis all these functions would be admitted. No further reduction is possible unless spatial symmetry is present and although the construction of spin eigenfunctions is a necessary step in the practical application of CI methods it is not sufficient to reduce the expansion to a manageable length unless the orbital basis is quite small.In order to make progress possible it is usually necessary to truncate further the expansion (l.l) e.g. by admitting only certain types of configuration ; one might for instance take a one-determinant SCF function as the leading term and then admit all configurations formed by promoting one or two electrons into the higher-energy " virtual " orbitals obtained from the SCF calculation. Unfortunately however 7 8 MULTI-CONFIGURATION SFC CALCULATIONS the orbitals which optimise a one-determinant wave function may not be well-suited to its further improvement by CI ; they do not give the best function of given form (unless complete CI is used) and a truncated expansion can be improved greatly in the usual variational sense by variation of the orbitals.In using a basis of m orbitals (with m large) it is therefore natural to turn attention from the extremely large but linear equation system for determining the CI expansion coefficients to the m-dimen- sional but non-linear problem of determining optimum orbitals in a given truncated expansion. This latter problem is in some ways analogous to that encountered in SCF theory and several types of multi-configuration (MC) SCF theory have now been successfully used. These go back originally to Hartree et al. and later to Yutsis et aZ.,2 for atoms but calculations on molecules are relatively r e ~ e n t . ~ - ~ Usually restrictions are imposed on the types of configuration admitted but general equations for the optimum orbitals can be given and the purpose of this note is to outline a practicable method of solution and to consider some of the incidental difficulties common to all MC calculations.2.-CHOI CE OF LINEARLY INDEPENDENT S P I N EI GENFUN CTIONS For a given orbital configuration there are many methods of setting up a complete linearly independent set of spin eigengunctions of given (S M ) . The functions we shall use are of the form a = M,4Q,@,l (2.1) where R is an orbital product 0 a spin function A the antisymmetrizer and M a normalizing factor. We choose 0 to be a spin eigenfunction with M = S since this is the simplest case to consider and matrix elements between functions with M # S can always be derived easily from this case (e.g.by the Wigner-Eckart theorem) and take @,(S1 s2 . . . s,) = q S i sj)e(sk s,) . . . a(s,)a(s,) . . . (2.2) where there are g paired spins each with a factor of the form and N-2g a factors. The resultant spin function has S = M = 3(N-2g) and the prescription used is essentially the one employed in valence bond (VB) theory. FIG. 1 . Branching diagram showing pos- sible resultant spin (S) for N electrons. 0 1 2 3 4 5 6 N If i2 contains a factor u(m)u(n) (doubly occupied orbital) then 0 must contain the factor Q(s s,J and spin factors associated with doubly-occupied orbitals need not be considered further. Different spin functions correspond to different pairing schemes but the number of possibilities is much larger than the number of independent functions.One method of eliminating linear dependencies is to set down the numbers 1,2 . . . N in a ring representing each 8(si s j ) factor by an arrow i-+j and each a(sJ R . MCWEENY 9 factor by a point p . Any functions with crossed arrows are then expressible in terms of those with no crossings and may therefore be discarded. When S # 0 however more limited dependencies still remain and it seems preferable to use a set which is in one-to-one correspondence with the standard set used in the branching-diagram and group-theoretical constr~ctions.~ In the branching diagram approach the spins are coupled one at a time,? the rrth electron spin being coupled to the spin eigenfunctions for 12- 1 electrons in order to obtain those for n electrons. The resultant functions are conveniently classified using the branching diagram (fig.l) in which the routes leading from the origin to any intersection represent the coupling schemes leading to linearly independent eigenfunctions of the corresponding S value. Thus for five electrons the coupling schemes leading to S = 3 are The functions so constructed are identical with those that carry the standard irreducible representations of the spin permutation group,8 associated with standard tableaux (respectively) (the numbers in the top row correspond to up-spin couplings (/) in the branching diagram thosein the bottom row to down-spin (\)) and the five functions are in one-to-one correspondence with (though not identical with) the VB-type functions indicated by In other words the VB function associated with any given branching-diagram function is obtained by coupling the first up-spin (/) with the first down-spin (\) the second up-spin with the second down-spin and so on any remaining spins being given a fix t or s .The final result is that the branching diagram functions (which are orthogonal not only independent) can be obtained from the VB functions by Schmidt ortho- gonalization starting at the left when the functions are listed (as above) in “last letter sequence ’,.lo The whole construction is now so straightforward and easily adapted for fully automatic computation that attention may be confined to functions of VB type which turn out to be extremely easy to handle. t Spins associated with doubly occupied orbitals need not be considered. 10 MULTI-CONFIGURATION SCF CALCULATIONS The evaluation of' matrix elements betwecn spin-paired functions for one- and two-electron operators of all kinds has been discussed by Cooper and McWeeny who give general rules based upon superposition patterns.Here we summarize the rules in their simplest form for spinless operators. Let us take two spin eigenfunctions 0 and OK. draw the associated diagrams showing the paired spins and superimpose these diagrams to form a superposition pattern. For example two 6-electron triplet functions ( S = 1) and their super- position pattern would be * 6 1 6 Such patterns contain islands each formed by a closed sequence of arrows andchains each formed by an open sequence the latter are of two types 0-chains containing an odd number of centres E-chains containing an even number. The pattern above contains one island and two 0-chains.The spin-factor dependence of any matrix element is characterized by such a pattern and in particular by the quantities. g = number of paired spins in each function n K K ' = number of islands in the superposition pattern v,,) = number of arrow reversals (if any) needed to achieve head-to-head tail-to-tail matching dE = 1 no E-chains; = 0 otherwise. The dependence on orbital factors which we denote by (where 4, is simply the orbital in the nth position whatever it might be) is characterized by the following quantities CT,,. = parity of permutation needed to achieve maximum matching of the orbital products inK,mKt = numbers of doubly occupied orbitals in R, QKt mKKt = number of doubly occupied orbitals which after matching S,,. = overlap of orbital products llK,QK) (= 1 products identical; coincide,? = 0 otherwise) SKK#(i)= overlap with orbitals in the ith position removed SKK#(i,j) = overlap with orbitals in the ith andjth positions removed qij(qij) = 1 if 4i t$j (& f#$) are different orbitals = 0 if 4i + j (& 4;) are identical.t For example AABBCDEEF and AABBCDCEF each contain three doubly occupied orbitals but can give only two coincidences ( I T z K K ~ = 2) ; but AABBCCDEF and AABBKKDEF with the same numbers of doubly occupied orbitals give three. R . MCWEENY 11 In t e r m of these quantities the non-orthogonality integral and the only non-zero niatrix elements of the usual spinless Hamiltonian take the forms (3.2) and (<Dtir I <Dk) = hE~,,tA,,tS,,~ <@,PI H I @ti) = &utiK*AKKr C sKtif(N4i I 11 I 4i) + 3SC'SKKf(i,j)(+:4i I I 4i4j) N 1 i = 1 i j = 1 The matrix element vanishes identically if there are more than two E-chains.The factors CO,,. and AKKI are defined by (3.4) COKK = (- lyKK'2(mK +inK'- 2 m K K ' ) / 2 The coefficients x i j depend only on the positions of i a n d j in the superposition pattern and are listed in table 1. TABLE 1.-COEFFICIENT Xc FOR POSITIONS i j I i xij(= X j i ) no E-chains island I island I %3Plj+ 1) 0-chain J 0-chain J W P i l + 1) island I island I' (#I) - i 1 0-chain J 0-chain J' (# J) + ( p ~ + 1) or 0-chain J two E-chains E-chain K E-chain K' (#K) p~ pq = pipj where the " parity factor " pi = -+1 is assigned by giving + 1 to an arbitrary position in an island or to an end point in a chain and then proceeding along the sequence giving rtl to alternate positions.For E-chains the end point chosen is the one where the arrow points into the chain. Normally the expression (3.3) contains only a few terms vanishing altogether when QK and R,. differ in more than two orbitals; thus with one orbital difference in the pth position SKKI = 0 SKKp (i) = 0 except for i = p and SKK#(i,j) = 0 unless i = p o r i = p. Examples of the application of these results are available.12 They are readily programmed for fully automatic computation the only input data required being a specification of orbitals and coupling scheme for each <DK. 4.-THE VARIATION PROBLEM When the orbitals are selected from a complete set A 23 C . . . each matrix element takes the form Optimum values of the expansion coefficients cK for a given choice of orbitals are determined by solving the secular equation HC = EMc (4.2) 12 MULTI-CONFIGURATION S C F CALCULATIONS with matrix elements The energy may then be written The coefficients of the integrals are elements of the one- and two-electron density matrices and depend on the expansion coefficients cK thus and these coefficients determine the orbital forms of the electron density and the pair function.The basic problem of MC SCF theory is to optimize the forms of the orbitals A B C . . . which appear in (4.4) so as to obtain an improved energy. For this purpose the coefficients of the one- and two-electron integrals are regarded as numerical constants and we consider an independent variation A+A+6A B+B+6B . . . of the orbitals. A formal solution of this problem in the form of equations defining the optimum orbitals was given many years ago but no practical application of these equations appears to have been made.To obtain a practicable method we introduce an LCAO-type approximation assuming that the n orbitals referred to in the CI expansion are written as linear combinations of nz basis orbitals $1 $ 2 . . . & Thus writing A = c $ r K A (4.7) r where the coefficients TrA ( A denoting a typical orbital) may be collected into an nz x JZ matrix T the energy expression (4.4) takes the form in which the one- and two-electron integrals are now defined over basis functions and the variational problem now involves the elements of the unknown matrix T. This matrix must be chosen so as to minimize E subject to an orthonormality con- straint which (assuming an orthonormal basis) may be written in the form TfT = 1 (4.9) (the n x y1 unit matrix).The equations that determine the optimum orbitals follow on making a first-order variation T+T+GT and requiring that E be stationary the constraints being in- corporated by the method of Lagrangian multipliers. The result is hTPl +Z = TE (4.10) where E is a Hermitian matrix of mutlipliers which may be determined from the condition (4.9) on multiplying (4.10) from the left by Tt. Here h is the usual matrix of one-electron integrals with elements hrs = (r I h I s) and Z is anelectron interaction matrix whose elements are zrA = C TJs<rs I 9 I tLL)TCT~DP2CD,AB* (4.11) B,C,D s.t,u R . MCWEENY 13 Equation (4.10) reduces to the Roothaan equation l4 hFT = (h+ G)T = TL (4.12) in the one-determinant approximation in which P1 is twice the unit matrix and Z may be factorized in the form Z = 2GT; but in general Z does not factorize and the equation does not reduce to eigenvalue form.5.-METHOD OF SOLU'I'ION Instead of trying to solve (4.10) directly it appears to be simpler to start froni the energy formula (4.8) and to minimize E by a descent procedure analogous to that used in some forniulations 1 3 9 of the one-determinant SCF method. To incorporate the orthonormality constraint (4.9) we suppose T-+T+6T = UT = (l+V)T (5- 1) where U is an HZ x IIZ unitary matrix and V must therefore be to first order anti- Hermitian Vt = -v (5.2) The corresponding first-order change in energy is where 6E(') = trVt(X-Xt) X = (hTP1 + Z)Tt (5.3) (5.4) The reqiiircinent 6E(' = 0 (all V) yields incidentally an alternative stationary condition namely that the matrix X be Hermitian.Away from the stationary point however X is non-Hermitian and we may regard rs as the Hermitian scalar product of two vectors in a space of dimension 172. For a V of given magnitude this quantity takes its greatest negative value corresponding to steepest descent of the energy surface when the two vectors are anti-parallel v = -A(X-Xf) (5.5) where A is a positive number measuring the magnitude of the change. The method of reaching an energy minimum is now clear. We start from a T matrix corresponding to the orbitals of a first approximation calculate the correction matrix V with some suitable estimate * of ;1 and revise T accordingly. The process may be repeated recalculation of Z X and V leading to a new descent direction and iteration may be continued until (X-Xt) vanishes to any desired accuracy.In order to preserve orthonormality accurately (not just to first order) T may be corrected between iterations by the process T,+T,+ = +T,(31 -T,?T,) which is rapidly convergent when T (used as first approximation) has columns which are almost orthonormal. * e.g. A may be chosen so that aE/aA = 0 corresponding to the achievement of the lowest point accessible in the given descent direction.l39 l 5 14 MULTI-CONFIGURAT ION SCF CALCULATIONS 6.-CONCLU S I ON A programme has been outlined for the fully automatic calculation of mutli- configuration molecular wave functions with optimization of the orbitals and with no restrictions on the types or configuraton admitted.The main steps are (i) construc- tion of a linearly independent set of spin eigenfunctions (ii) evaluation of matrix elements and solution of the secular equations and (iii) iterative adjustment of orbitals to obtain the best CI expansion of given form. Algorithms have been devised for the computer implementation of each step and pilot calculations on small molecules are now in progress. Convergence has been found satisfactory and detailed results will be published elsewhere. D. R. Hartree W. Hartree and B. Swirles Phil. Trans. A 1939 238 229. Proc. Roy. SOC. A 1933 139 311. A. P. Yutsis Zhur. Eksp. Teor. Fiz. (Soviet Phys. JETP) 1952 23 129. G. Das and A. C. Wahl J. Chem. Physics 1966,44,87. E. Clementi and A. Veillard J. Chern. Physics 1966 44 3050. J. Hinze and C. C. J. Roothaan Suppl. Prog. Theor. Physics 1967 40 37. ti R. McWeeny Proc. Roy. SOC. A 1955,232 114. ’ M. Kotani A. Amemiya E. Ishiguro and T. Kimura Tables of Molecular Zntegrals (Maruzen. Tokyo 1966). T. Yamanouchi Proc. Phys. Math. SOC. Japan 1937 19,436. R. McWeeny 1968 unpublished. D. E. Rutherford Substitutional Analysis (Edinburgh University Press 1948). I1 I. L. Cooper and R. McWeeny J. Chem. Physics 1966,45,226. I. L. Cooper and R. McWeeny J. Chem. Physics 1968,49,3223. I3 R. McWeeny Rev. Mod. Physics 1960 32 335. I4 C. C . J. Roothaan Rev. Mod. Physics 1951 23 69. R. McWeeny Proc. Roy. Soc. A 1966 235 496.