J. Phys. Chem. 1982, 86, 1059-1064
remark that a useful extension of the present method involves the calculation of intermolecular potentials as a function of the internal molecular coordinates.N This should provide us with a reliable basis for the calculation (34) W. Rodwell, R.0. Watts, and S. F. Lam,Mol. Phys., submitted for publication.
1059
of vibrationally inelastic scattering for a larger variety of systems than available at present.
Acknowledgment. We are most grateful to Professor U. Buck and K. T. Tang for supplying us with details of their Work28B*prior to publication and to Dr*R*0.Watts and Professor D. P. Craig for helpful discussions.
Molecular Electronic Structure by Combination of Fragments Bernard Klrtman Department of Chemistry, University of Callfornla, Santa Barbara, Californla 93106 (Received: April 27, 1981; I n Final Form: August 5, 1981)
A method of combining fragments to obtain the electronic structure of a large molecule is presented. Our technique is based on an SCF treatment of the first-order density matrix, R. It is applicable within the Hartree-Fock, strongly orthogonal valence-bond, and density functional models. Primary corrections to R, which occur in the vicinity of the bond(s) between the fragments, are calculated exactly; secondary delocalization effects are determined approximately. The approximations result in a substantial saving of computer time and storage requirements. Initial tests with INDO and the unrestricted Hartree-Fock model show small errors for conformationalenergy differencesand the total density matrix. The computation time increases quadratically with the size of the molecule rather than as the third power which is the case for a conventional INDO calculation and for the related PCILO method. Unlike PCILO, however, our procedure is not inextricably bound to the zero differential overlap approximation. It can be used to substantially improve the efficiency of existing SCF treatments of large molecules and to increase their applicability as well.
I. Introduction There is much evidence to suggest that ordinary covalent molecules characteristically contain localized chemical bonds. Thus, it is natural to assume that the primary changes in electronic structure accompanying the formation of such bonds will occur within a limited region of space. Yet it is not obvious how to take advantage of this fact to simplify structure calculations. We have recently' developed a method which does so-at least for localized interactions in extended systems like those encountered in chemisorption. In the present article our approach is extended to the problem of determining molecular electronic structure by combining two or more fragments. If successful, such a treatment would have numerous applications including the study of various biological systems, interactions in aqueous solution, etc. The method presented here, in section 11, is based on the single-particle density matrix (hereafter referred to as the density matrix). It can be used in conjunction with the HartreeFock or unrestricted Hartree-Fock model, the strongly orthogonal valence-bond model, or any of the density functional approximations to the Hamiltonian. We start with the density matrices of the separate fragments that we wish to combine to form a molecule. Of course, the density matrix of the molecule is not simply the sum of the fragment density matrices (even after nonorthogonality has been taken into account). It may be assumed, however, that the primary corrections are well described by a relatively small orbital basis (= the local space) centered in or near the bonding region. But such a limited basis set will not allow for the proper transfer of charge into or out of the bonding region nor will it maintain the
molecular orbitals at unit occupancy. In other words, the density matrix will not be idempotent. A secondary correction must be included to account for this deficiency. It turns out that the idempotency requirement can be incorporated exactly without enlarging the local space by utilizing an effective Hamiltonian. Thus, our method effectively includes some variation over the entire space. Indeed, in section IV we present an example where the calculations yield exact results for the complete molecular density matrix. The major reason for introducing a density matrix formalism is the ease with which these secondary idempotency effects are taken into account. As one would expect, the local density matrix method may be simplified by introducing various common approximations such as pseudopotentials, partial or complete neglect of differential overlap, etc. In this paper we employ Pople's INDO approximation2 in conjunction with the unrestricted Hartree-Fock model to carry out some illustrative calculations on methane and ethane. In each case the molecule was divided into two fragments. If the local space is restricted to orbitals on the two atoms that connect the fragments, then the total charge densities and conformational energy differences are very well reproduced in comparison with a complete treatment. In addition, the energy of the connecting C-C or C-H bond is duplicated to within 6 kcal/mol in the worst case. In section I11 we analyze the savings in computer time that result from limiting the orbital basis to a local space. In a conventional Hartree-Fock calculation the process of diagonalizing the Hamiltonian matrix to fiid an improved set of molecular orbitals is third order in the dimension of the basis. For the INDO approximation this is the
(1) B. Kirtman and C. deMelo, J . Chem. Phys., 75, 4592 (1981).
(2) See J. A. Pople and D. L. Beveridge, 'Approximate Molecular Orbital Theory", McGraw-Hill, New York, 1971, p 81.
0022-3654/82/2086-1059$01.25/0
0 1982 American Chemical Society
Kirtman
The Journal of Physical Chemistty, Vol. 86, No. 7, 1982
1060
rate-limiting step. In the local density matrix method the corresponding step is more rapid by a factor of (N/NQ)~ where N is the average number of orbitals per fragment in the complete basis set and NQis the number in the local space. This step will remain the bottleneck if we account for an increase in the size of the molecule by increasing the number of fragments, f , while holding N fixed. An alternative strategy is to increase f and N proportionately. In that event the bottleneck lies elsewhere and the computation time depends quadratically on the size of the molecule as opposed to the third power for a conventional INDO calculation and the related PCILO m e t h ~ d .Be~ sides the saving of time there is also a substantial reduction in storage requirements. A variety of procedures have been developed for dealing with the electronic structure of biomolecules and other large systems. Within the zero differential overlap (ZDO) approximation our method is an alternative to PCILO. However, we are not bound to ZDO, or to just s and p orbitals, or to a particular zeroth-order bonding description, or to a single structure, et^.^ as PCILO is. The local space approximation is simply a technique for reducing the computations involved in setting up and solving the SCF equations regardless of where the molecular integrals come from. Thus, it can be used in conjunction with any of the existing SCF procedures. As discussed in section IV the latter are, currently (or will shortly be), limited by the SCF rather than the integration step. So there is the potential for significant impact with regard to efficiency. In addition, we show how the range of applicability can be extended by generalizing the definition of the local space. 11. Theoretical Basis It is well-known that in the Hartree-Fock and density functional models the energy is completely determined by the (first-order) density matrix, R. A similar statement4 may be made for the strongly orthogonal valence bond model as well. Within each of these models, then, the usual variation condition on the energy leads to an equation for R which, in turn, determines all the physical properties of the system. In order to be correct, however, this equation must provide for the fact that only certain variations of the density matrix are allowed. For R must be idempotent5 to ensure that it is derivable from an orbital model wave function. McWeenyG has written the idempotency condition that must be satisfied in the general form
R = R,
+ AR = (R, + v ) ( l + vtv)-'(Ro + vt)
(1)
where v = UOXR,
Uo = 1 - Ro (2) and X is completely arbitrary. Here Ro is the original (idempotent) density matrix and AR is an allowed correction to it. The latter can also be expressed as a power series in X by expanding (1 + vtv)-' to get' AR = v + + V V -~ V ~ + V ... (3) ~f
(3) (a) S. Diner, J. P. Malrieu, and P. Claverie, Theor. Chim. Acta, 13,
1 (1969); J. P. Malrieu, P. Claverie, and S. Diner, ibid., 13, 18 (1969); S.
Diner, J. P. Malrieu, F. Jordan, and M. Gilbert, ibid., 15,100 (1969). (b) For a recent review we J. P. Malrieu in 'Modern Theoretical Chemistry", Vol. 7, G. A. Segal, Ed., Plenum Press, New York, 1977. (4) B. Kirtman, work in progress. Strictly speaking we should say, in this case, that both the first-order density matrix and the energy are derivable from elementary single-particle matrices. (5) For the valence bond model it is not the denaity matrix itaelf but some closely related matrices that must be idempotent. (6) R. McWeeny, Reu. Mod.Phys., 32, 335 (1960).
This transformation shifts attention from AR to the matrix X which is freely variable. If Ro is some initial guess for the density matrix, then an approximate correction can be determined from the linear variation condition which leads to the relation &hUo + Uoh& = Roh&XUo + UoXRohRo UohUoXRo - RoXUohUo + O(X2) (4) as we have shown.'V8 In the above h = h(R) is the model Hamiltonian matrix. It may be spin dependent in which case &,Uo,and X are spin dependent as well. In addition, h is a function of R and, hence, of X. For this reason eq 4 has to be solved self-consistently starting with h = h(&), even when terms of O(X2) are ignored (as they are in our procedure). Each self-consistent iteration yields a new X or v and, from eq 1, an improved R which, then, is used to update h. The matrix v is most easily determined through the intermediary Y = RJUo + UJR,,. In terms of the latter, eq 4 becomes &hUo + Uoh& = (%h& - UohUo)Y+ Y(&h& - UohUo)+ O(Y2) (5a) and, of course v = UoYRo
(5b) By transforming eq 5a to the basis which diagonalizes &hRo - UohUowe immediately find Yij = (R,hUo + U&&)ij/(Xi + Xj) (6) where X i , Xj are eigenvalues of RohRo- UohUo. For our treatment we wish to limit the variation to some local space. This space is defined in terms of a limited set of atomic orbitals, &, centered in the region where the two fragments are connected. Since the orbitals have tails that extend to infinity they are not completely localized but they are localized in exactly the same sense as_ conventional localized molecular orbitals. The projector Q for the local space is defined as
Q
= EIdi)(S-1)ij(4jI ij
(7)
in which S-' is the inverse overlap matrix, and the projection onto this space will be indicated by a subscript &. Thus, for example XQ = QXQ (8) where Q is the matrix representation of Q. Although AR in eq 1 is idempotent its projection onto the local space, in general, will not be. In order to limit the variation, therefore, we deal directly with X. The replacement of X by XQ constitutes the local space approximation. Although XQ is local it also gives rise to secondary, or delocalized, corrections to the density matrix. The latter occur as a result of the Uo and Ro projections in eq 2. Among other things the secondary effect allows for the self-consistent flow of electronic charge into and out of the local space. Despite the fact that AR is nonlocal the linear variation condition yields a completely local analogue of eq 4, i.e. (hRu+ hUR)Q = hQRxQ(U&+ (Uo)QXQhQR hQuxQ(&)Q - (R&XQhQU + ~ ( X Q ' )(9) (7) The first-order terms in this expression are unique but the second-order terms are not (see ref 1). Naturally, this means that eq 1 is not unique either. (8) Here we have set Xt= X which is readily proved to be so. The same relation has also been derived by McWeeny in another form. See R. McWeeny, Phys. Reo., 126, 1028 (1962). eq 4.3.
The Journal of Physical Chemistry, Vol. 86, No. 7, 1982
Electronic Structure from Fragment Combination
in which
bU= (UohUo)Q (hRU+ hUR)Q (RohUo + Uoh&)Q
hQR= (RohRo)Q
(loa) (lob)
A detailed derivation of this result has been given previously.' The linear relation for XQwhich is obtained by ignoring terms of O(XQ2)in eq 9 cannot be solved directly. However, we have found that a reasonable first approximation is determined by the local analogue of eq 5a, i.e.
(hRU+ hUR)Q = (hQR - hQU)YQ + YQ(hQR - heu) (1h) with
YQ = (UO)QXQ(RJQ + (&)QXQ(UO)Q (1lb) Eguation l l a leads to an approximate YQand eq I l b yields the corresponding XQ after transforming to the basis in which (R& and (Uo)Q are simultaneously diagonal. This initial approximation to XQcan be improved by iteration but, in practice, that has not been found necessary (see section IV). It is possible to convert the analogue of eq 1into a very convenient form with the aid of the local matrices uQR,uQ", etc. defined by
AR = R0uQRRo+ ROUQ~~UO + UOUQ~~RO + UOUQ~UO (12)
Using the power series expansion of (1 + vtv)-' as an intermediate one can readily verify that = - u Q ~ ~ ( U ~ ) Q XUQQ ' = XQ(RO)QUQ~" (13a)
U Q ~
with OBRU
(13b) (uO)QxQ(&)QxQ]-'
= xQ[1Q +
and, of course, .BUR = (ugRU)'. In carrying out the self-consistent interations we have found that it is easiest to accumulate the total change in the density matrix with respect to the very first approximation, Roo, appropriate to the noninteracting fragments. For the INDO scheme Roo is simply the sum of the density matrices of the isolated fragments. In the general case the nonvanishing overlaps between orbitals of different fragments must be accounted for by symmetric orthogonalization. The cumulative change in the density matrix may be expressed as
Then, on the first iteration, 7QR= uQR,7QRU= uQRU,etc. and the new approximation to the density matrix is Rol = R,,O + AR. Given AR the corresponding total change in the Hamiltonian AhT and, particularly, the local projections [&OA~TR~]Q, [&OAhTU,O]Q,etc. are readily determined as discussed in the next section. For the second iteration (see eq 9, loa, and lob) we do not require the latter projections but rather the current ones where Roo and Uooare replaced by R,,' and Uo'. When the four 76 matrices or, better yet, the combinations KQ
+
= (Roo)Q7QR(Uo0)Q7gUR
(154
X, = 7Qu(Uoo)Q + 7QUR(h0) (15b) are used, the required transformation can be done by manipulations entirely on the local space. Thus, for example
1061
[%'AhT&'lQ = (IQ + K Q ) ~ ( [ % ~ A ~ T &+~ KQ) ] Q (+~ [boAh~Roo]~X,,) Q + ~t([u~oAh~ + u[ U~Oo~]A~~ ~T R O+~ KIQQ) ) ( (16) ~Q At this point the SCF equation (eq 9) is solved once again and a new set of UQ matrices (eq 12) is obtained. From these the new set of TQ'S is generated by transforming back from the current R,, (in this case, R,,') to R,,O utilizing the previous TQ'S. The next step is to update AhT and, at the Same time, to COnStrUCt [Ro0Ah&O]~,[R$A~TUO~]Q, etC. Finally, the iterative looping is continued in the same vein until the solution has converged. From the converged solution the Hartree-Fock energy is readily evaluated. An appropriate formula for the restricted case has been presented in a previous paper.' For the unrestricted model the energy expression is E(HF) = Eo(HF) + y2 Tr [(h, + h)ART].+ 7 2 Tr [(ho+ h)aRrlo (17) where the subscripts CY and p refer to electron spin and Eo(HF) is the Hartree-Fock energy calculated with the initial approximation to the density matrix. Likewise ho is the Hamiltonian for R = &O. The last two terms on the right-hand side of the above expression can be reduced to sums over the local space by taking advantage of eq 14 which yields Tr [(ho + h)A%1, = Tr([&(ho h)Rolg7gR+ [Uo(ho + h)ROlQ7QRU + [Uo(ho + h)uOlQ7Qu,)a (18) [Who h)UOlQ7QUR and an exactly analogous formula for spin 0. Within the INDO approximation Eo(HF)is the sum of the HartreeFock energies of the separate fragments plus an interaction term. The latter includes the interaction between atomic cores in different fragments and the effective electronnuclear attraction (or repulsion) between fragments due to the net charges on the atoms. Exactly the same terms are contained in ho as well. Although the INDO approximation is by no means necessary it does provide considerable simplification and is useful for some purposes. We employ it here for illustration and preliminary testing as well. The relevant Hamiltonian matrix element expressions are well-known2 and will not be given here. Instead, we will briefly discuss how the Hamiltonian is updated from one interation to the next. It is important to realize that one does not need to store the full Hamiltonian matrix. All that is required for haU,and (hRu+ huR)s. eq 9 and 17 are the projections bR, Although these projections involve the current density matrix they can also be calculated from projections in terms of the initial density matrix as illustrated by eq 16. The latter projections are evaluated directly in the following manner. From eq 14 the set of corrections to the density matrix (A&)rs is determined for a particular 4r and the set of orbitals 4s.Here and 4, are members of the complete orbital space of the entire molecule. This set of corrections to the density matrix yields a corresponding set of corrections to the Hamiltonian matrix, (Ah&,, as determined by the INDO approximation. Then the sum ~,(AhT)rs(&o),b is formed for each '#Ib in the local space. Because of zero differential overlap the index s includes only those orbitals associated with the same molecular fragment as I&. Finally, the product of yields the contribution to and ~s(AhT)rs(Roo)8b [RooAhTRoo],b due to the orbital &. This contribution is stored and accumulated for all r in a matrix which has the dimensions of the local space. An exactly analogous procedure is followed for the three remaining projections C#J~
1082 The Journal of Physical Chemisfry, Vol. 86, No. 7, 7982
Kirtman
obtained by replacing one or both of the Roo'swith Uoo.
about the nature of the perturbation series. Although PCILO is more rapid than a conventional CNDO calcu111. Computational Analysis lation, in either case the computation time increases with By limiting the primary variation of the density matrix the number of orbitals as NT3. to a local space we achieve a considerable reduction in In addition to the fact that PCILO is inextricably tied computer time, not to mention storage requirements. In to the zero differential overlap approximation there are order to carry out an analysis let us suppose, for sake of a number of other limitations which have been reviewed argument, that there are f fragments and that, on the by M a l r i e ~ These . ~ ~ include the convergence of the peraverage, the complete basis for each fragment contains N turbation series problem mentioned above, the arbitrariorbitals whereas the local space consists of NQ orbitals. In ness of the zeroth-order description, the inability to deal order to initiate our procedure it is necessary to have the with the dissociation of covalent bonds or with the tranresults of a conventional calculation on each of the fragsition from one bonding structure to another, and the ments. For a HartreeFock calculation the bottleneck, at restriction to s and p orbitals. In principle, at least, none least within the INDO approximation, lies in the diagoof these limitations are inherent to the local density matrix nalization of the Hamiltonian matrix. Both this step and method. However, the latter does utilize a local space the construction of the updated density matrix are N3 approximation and it does not include as much of the processes whereas the number of nonzero (INDO) integrals correlation energy'O as PCILO. Of course, the extra coris of order N2. Thus, for f fragments, the preliminary relation energy does not guarantee1' better conformational computation time will be proportional to fiV. On the other energy differences or potential energy features in general. hand, a direct calculation on the entire molecule would be The detailed analysis given above pertains specifically an PIP process. The savings factor is f which indicates to an INDO calculation. But the major time and storage that it is desirable to use as many fragments as po~sible.~ saving features will persist for any SCF treatment. For However, this must be balanced against the time spent on they derive from the fact that the Fock matrix is subthe subsequent density matrix treatment (see below) which stantially reduced in size due to the local space approxiincreases with the number of fragments for a molecule of mation. In particular, the time required for diagonalization fixed size. . significance of is decreased by a factor of ( N Q / N ) 3 The The rate-determining steps in the local density matrix this improvement is discussed in section V. calculation are updating the Hamiltonian matrix and IV. Illustrative Calculations solving the HartreeFock equation for h. In the previous As an initial test of our computer programs and of the section we described the sequence of manipulations accuracy of the local space approximation we have carried whereby the updating is accomplished. The most timeout a series of unrestricted INDO calculations on methane consuming of these is the computation of ( A b ) and and ethane. The molecule is divided into two fragments cs(AhT)rn(&o)& (or zn(AhT)rs(u$)sb)which is an P"NQ in each case although the definition of the fragments is process within the INDO approximation. In order to solve varied as are the orbitals contained in the local space. In the Hartree-Fock equation for X Q two matrix diagonalithe local density matrix calculations an iterative procedure zations and several matrix multiplications on the local occurs in two places. There are, first of all, the overall space are required. The number of individual arithmetic Hartree-Fock SCF iterations. Then, within each SCF operations is, therefore, proportional to PNQ3and does not iteration, the current Hartree-Fock equation is solved by depend at all on the average size of the fragments. The a successive approximation method. In order to increase relative timing of the two rate-determining steps will deefficiency the convergence criterion for the inner iterative pend upon the strategy that is adopted. If N is held fixed loop is adjusted from one SCF cycle to the next. As the as the molecule increases in size, then the f N Q 3process change in the density matrix gets smaller for successive will eventually become dominant. In fact, it will dominate SCF iterations this criterion is tightened proportionately. the entire calculation (including the preliminary treatment In practice we have found that the initial approximation of the separate fragments), assuming that NQ is constant for the solution of the current Hartree-Fock equation has as well. But, if the number of fragments and NQ are fixed, always been sufficient. then the preliminary f M step will eventually take over with A. Methane. In the case of methane the molecule was the updating, which is OcfWNQ)being second most imdivided for all calculations into a methyl radical and a portant. As a compromise, we may consider the case where hydrogen atom. An isolated hydrogen atom contains one both f and N are proportionately increased along with the electron of spin a,say, and zero electrons of spin /3. On size of the entire molecule. In that event the overall timing the other hand, the self-consistent electronic charge at each formula will be hydrogen atom in methane turns out to be 0.505 for both t = aNT2 i- bNT3I2 (19) spins. So the initial approximation is fairly poor in this regard. If Q is taken as the projector for the full space, where NT is the total number of orbitals and the second then our results are identical with those obtained in a term arises from the ~ N step. Q ~ Thus, the local space conventional INDO treatment just as they should be. approximation reduces the dependence upon NT from When the orbitals on the hydrogen atoms of the methyl cubic to quadratic. group are removed from the local space the energy and the It is of interest to compare our technique with the density matrix-or, at least, the local projection of the PCILO m e t h ~ d . ~As mentioned in the Introduction density matrix-are unchanged. So in this particular inPCILO was especially designed for use in conjunction with the zero differential overlap approximation. Within this (10) Some correlation energy is included in the local density matrix context it is an N T procedure ~ through second order and method when the strongly orthogonal valence bond model is employed. NT3 through third order. Normally third-order terms are Furthermore, the unrestricted Hartree-Fock treatment described below included since they may be quite substantial. Indeed, the is correlated with respect to the restricted treatment used (for zeroth order) in the PCILO method. fourth-order terms can be larger yet which raises questions (9) Of course, the fragments must remain large enough so that the different local spaces do not have any orbitals in common.
(11) Particularly since the correlation is calculated in the CNDO approximation. A case in point is the barrier to internal rotation in ethane. The PCILO value turns out to be further from experiment than the corresponding SCF value.
Electronic Structure from Fragment Combination
The Journal of Physical Chemistry, Vol. 86, No. 7, 1982
1063
TABLE I: Projected Density Matrix for Ethane Formed by Combination of Two CH, Fragmentsa 2s 2s 2PX 2PY 2PZ 2s' 2Px; 2PY, 2PZ 2s 2PX 2PY 2PZ 2s' 2Px; 2PY, 2PZ 2s 2PX 2PY 2PZ 2s' 2Px; 2PY, 2PZ
2PX
0.520 0.000 0.000 0.012 0.127
0.000
0.000 0.000
0.088 0.000
-0.232
0.000
0.521
0.000
0.000
0.475
0.000 -0.013 0.127
0.000 0.000 0.000
0.000
0.088 0.000
-
0.000 -0.231
0.472 0.000 0.000
0.000
2 PZ
2s'
A. Full Space, a or p 0.000 -0.012 0.000 0.000 0.472 0.000 0.000 0.507 0.000 0.232 0.000 0.000 0.088
0.000
0.000 -0.356
Spin 0.127 0.000 0.000 0.232 0.520 0.000 0.000 0.012
B. Local Space, a Spin 0.000 - 0.013 0.127 0.000 0.000 0.000 0.47 5 0.000 0.000 0.000 0.508 0.232 0.000 0.232 0.520 0.000 0.000 0.000
-0.232 0.000 0.000 -0.356 0.012 0.000 0.000 0.507
0.000 0.000 0.473 0.000
0.000
0.000 0.088 0.000
0.000 0.000
0.000 0.473 0.000
0.000 0.000 0.088
0.000 0.000 0.000
0.000
0.000
0.088
0.000
0.000 0.088
0.506 0.231 0.000 0.000 -0.356
0.231 0.521 0.000 0.000 0.013
0.470
0.000 0.000 0.000
0.000
0.000 0.088 0.000 0.000
0.000 0.000
0.000
0.470
0.000
0.088 0.000
0.127
0.000
0.000
0.000
0.000
0.000 0.000
0.520
2PZ'
0.470
p Spin - 0.01 1 0.000
0.000
2Py'
0.000
0.000
0.000 0.011
0.088
2PxI
0.470 0.000
0.000 -0.356
0.000 0.000 -0.011 0.127
2PY
0.000 0.000 0.000
0.000 0.000 0.000
-0.231 0.000 0.000 -0.356 0.011 0.000 0.000 0.506 -0.231 0.000 0.000 -0.356 0.013
0.475 0.000 0.000 0.088 0.000 0.475 0.000 0.232 0.000 0.000 0.000 0.000 0.508 a The local space consists of the valence shell orbitals o n both carbons. In case A the calculation is done on the full space and projected afterwards; in case B it is done on the local space in the first place.
0.000 0.000
0.088
stance the local space approximation is exact.12 A further removal of one or more of the carbon orbitals causes the agreement to vanish. For example, deletion of a 2p orbital (the 2p orbitals were taken to be perpendicular to the faces of a cube whereas the hydrogen atoms lie at the corners) causes the energy to increase by 61 kcal/mol. The major effect on the density matrix is to introduce an artificial spin splitting. For one spin the electronic charge at the hydrogen atom of the original hydrogen atom fragment is 0.582; for the other spin it is 0.422. The average of these values (0.502) is very close to the full space result quoted earlier (0.505) where the two spins are equivalent. A similar situation exists for the electronic charge in the 2s and 2p orbitals of carbon although the splitting is smaller. In other words, the total density matrix remains quite accurate even though an artificial spin splitting does appear. Of course, we could have taken advantage of the symmetry of methane to eliminate the two 2p, orbitals on carbon without affecting the density matrix at all. But the above calculation is more representative of the general situation that will be encountered. The error in our treatment is determined by how well the exact (Le., full space) orbitals can be localized to the particular subspace of interest. A quantitative measure of the degree of localization is given by the eigenvalues of the projected density matrix q.These eigenvalues represent maximum and minimum occupancies for orbitals confined to the local subspace. For the full space there are four orbitals of unit occupancy and four orbitals of zero occupancy for each spin. When one or more of the hydrogen atom orbitals of the methyl fragment are removed most of the eigenvalues deviate from zero and unity. In fact, for both spins just one totally occupied and one completely unoccupied orbital remains. But that is sufficient to make the local space approximation exact as we have seen. Any further reduction in size causes the unit (12)To four significant figures.
occupancy to disappear in one spin space and the zero occupancy to disappear in the other. At the same time the local space approximation ceases to be exact. B. Ethane. For ethane we examined the effect of the local space approximation on the barrier to internal rotation as well as the bond energy. The barrier determination represents a typical conformational problem. Calculations were done at the equilibrium geometry of ethane with two different fragment combinations CH,. CH3. and C2H,- + He. In both instances we confmed that, for the full space, our results agree with the conventional treatment. For the combination of two methyl radicals the smallest reasonable subspace (ignoring symmetry) consists of the 2s and 2p orbitals on both carbons. The density matrices computed for a and 0spin are given in Table I along with the corresponding projections of the full space calculation. A small error does arise here due to neglecting the hydrogen atom orbitals. But the individual density matrix elements differ from the full space values by no more than 0.003 and the total energy is raised by only 0.1 kcal/mol. Just about the same deviations are found in the eclipsed conformation. Thus, the discrepancy in the barrier height is even less than13 0.1 kcal/mol. For the combination of an ethyl radical with a hydrogen atom the smallest reasonable subspace consists of the valence shell orbitals of the carbon and hydrogen involved in the bond connecting the fragments. This time the size of the local space is only 5/14 that of the full space. (If we take symmetry into account the corresponding ratios are 3/4 for methane and 2/3 for CH3. + CH,..) As far as the density matrix is concerned the major effect of the approximation is, again, to introduce an artificial splitting between the CY and 0 spin. In fact, it can be seen from Table I1 that the error in the spin-averaged density matrix is less than 0.004 for each element. The energy is raised
+
(13) It is, actually, on the order of 0.01 kcal/mol.
1084
The Journal of Physical Chemistry, Vol. 86, No. 7, 1982
TABLE 11: Projected Density Matrix f o r Ethane Formed by Combination of an Ethyl Radical and a Hydrogen Atoma Is
2s
2P,
2PY
2 Pz 0.152 0.012 0.000
A. Full Space, 01 or p Is 2s 2p, 2p, 2p,
0.509 0.245 0.201 0.348 0.152
Is
0.525 0.243 0.202 0.350 0.154
0.245 0.520
0.201
0.348
0.000
0.000 0.000
0.472
0.000 0.000
0.000 0.000
0.472 0.002
0.012
0.000 0.507
B. Local Space, cy Spin 2s 2p, 2p, 2P,
0.243 0.518 -0.002 -0.003 0.011
0.202 -0.002 0.469 -0.003
0.350 -0.003 -0.003 0.466
-0.001
-0.001
0.154 0.011 -0.001
-0.001 0.505
p Spin 0.251 0.200 0.347 0.152 0.487 0.004 0.013 2s 0.251 0.521 0.002 0.005 0.002 0.002 0.477 2p, 0.200 0.483 0.003 0.004 0.005 2p, 0.341 0.003 0.510 0.013 0.002 2P, 0.152 a The local space consists of t h e valence shell orbitals o f the carbon and hydrogen involved in the bond connecting the fragments. In case A the calculation is done o n the full space and projected afterwards; in case B it is done o n the local space in the first place.
Is
by 6.0 kcal/mol with respect to the full space but, again, the eclipsed conformation suffers a similar fate and the resulting barrier height is correct to within 0.33 kcai/mol. This is a relatively small error especially considering that the choice of fragments is a poor one for the barrier property.
V. Discussion It is important to realize that our method is not bound
to the ZDO approximation nor to any particular technique for obtaining molecular integrals. What we have developed is an approximation procedure which reduces the computational effort involved in setting up and solving the SCF equations. This procedure can be employed in conjunction with any of the SCF treatments that have heretofore been proposed for large molecules. In order for the local space approximation to have a large impact the time required for integral evaluation must be a minor fraction of the total computation time. This appears to be the situation that currently exists or, at the very least, is rapidly approaching.I4 It is certainly true of the FSGO-fragment treatment devised by Christofferson et al.15 and of the Halgren and Lipscomb16PRDDO approximation scheme as well. It is (14)J. J. Kaufman, H. E. Popkie, and P. C. Hariharan, ACS Symp. Ser., No. 112,415 (1979). (15) For relative timing of the two parta of a typical FSGO-fragment calculation see R. P. Angeli, S. D. Hornung, and R. E. Christofferson in ‘Proceedings of the Third International Conference on Computers in Chemical Research, Education, and Technology”, Caracas, Venezuela, 1976,E. V. Ludena, N. H. Sabelli, and A. C. Wahl, Ed., Plenum Press, New York, 1977,p 357. In the FSGO-fragment method the fragments are ueed to determine a transferable FSGO basis which is a very different approach from ours. References for the detailed procedure are given in the above article.
Kirtman
also true whenever one is dealing with chemical or geometrical changes within a physically small region of a spatially extended molecule. For in that event most of the integrals need be computed only once. (In this connection a library of integrals could be built up for many of the common, or not so common, fragments that occur in large molecules.) From the foregoing discussion it is apparent that the method presented here will be especially advantageous for the study of conformational energy differences and of geometry variations in general. In this respect it is complementary to the non-self-consistent SAMOl’ technique which is highly inaccurate with regard to the former (i.e., differences between conformers) and relatively expensive for the latter. Of course, the accuracy of the various SCF treatments depends critically upon the basis set. The smaller the basis the more careful one must be that it is well balanced. In the FSGO-fragment approach this poses a severe difficulty in connection with conformational and other geometry problems. (A suggestion for improvement will be made below.) Note that the generalization of our method to include generalized valence bond wave functions permits the correct treatment of bond dissociation as well as those configurations where the mixing of ionic and covalent structures is important. In this paper we have focussed on the use of the local space approximation to join molecular fragments. Thus, for our examples, the local space was defined by the set of atomic orbitals on the atoms forming the bond. For comparable accuracy with larger fragments it will likely be necessary to include orbitals on neighboring atoms as well and, of course, on atoms that are sterically interacting. Indeed, the choice of the local space is completely flexible. It could be defined so as to selectively study what additions to the FSGO-fragment basis are needed to achieve a proper balance for the computation of a conformational energy difference. Alternatively, Q might span the entire basis of a relatively small fragment that is attached to a large molecular skeleton. This would be particularly appropriate for investigating a series of changes where the skeleton is invariant. And so forth. Although the applications mentioned here appear promising further calculations will be necessary to establish the accuracy and efficiency of the local space approximation.
Acknowledgment. I acknowledge the very helpful critique of this paper provided by William Palke. (16)T.A. Halgren and W. N. Lipscomb, J. Chem. Phys., 58, 1569 (1973);T.A. Halgren, D. A. Kleier, J. H. Hall, Jr., L. D. Brown, and W. N. Lipscomb, J. Am. Chem. Soc., 100,6595(1978). Actually, the statement made here is valid only up to about 95 electrons because the PRDDO method involves an NT‘ integral transformation from AO’sto orthogonalized AO’s. As indicated by Halgren and Lipscomb it is likely that the NT‘dependence can be eliminated without significant loss of accuracy. Our estimate for the break even point in the number of electrona was obtained by using (i) the ratio 5/3for the SCF w. integral times (not including the NT4term) as given for formic acid; (ii) the latest timing ~ and (iii) a ratio of 6/5 for the number formula for the N T contribution; of electrons as compared to the number of orbitals (this is the correct value for carbon). (17) For an excellent review of SAM0 which critically evaluates this approach and compares it to SCF methods see B. OLeary, B. J. Duke, and J. E. Eilers, Adu. Quantum Chem., 9,1 (1975).
