Journal of Molecular Structure (Theo&m), 260 (1992) 31%322 Elsevier Science Publishers B.V., Amsterdam
“Slide-rule” basis sets: application helium-like species* Graham Doggett”
313
to the 3S state of
and Graham D. Fletcherb
“Department of Chemistry, University of York, York YOl 5DD (UK) bDepartment of Physics, University of York, York YOI 500 (UK) (Received 4 November 1991)
Abstract The energy of a single determinant wavefunction with maximum total spin is usually determined by imposing the computationally convenient requirement of orthogonality on the atomic orbital functions, which are themselves expanded over a suitable basis set in the self-consistent-field (SCF) model. An alternative to orthogonalisation, yielding the same total energies, is presented here in which, for helium-like species, the chosen basis (a scaled even-tempered basis of N 1s Slater-type orbitals) is partitioned into three disjoint subsets P, W and 9. For subsets dl and 9 containing one element, it is demonstrated that the total energy is invariant to the choice of partitioning, and that the two (non-orthogonal) atomic orbitals may be expanded over N - 1 basis functions contained in D u W and V u 9, respectively. The total energies for the % state, as a function of overall scale factor, are compared with those obtained from full configuration interaction calculations, and the radial distribution functions for the 2iS and ‘S states are presented for helium.
INTRODUCTION
Coulson [l] has explored some of the underlying principles and assumptions of the HartreeFock model as applied to the helium atom, in his paper on the saddle point character of Hartree-Fock wavefunctions. In particular, he demonstrated that the energy of the simple Hylleraas split-shell 53 wavefunction with spatial part $ll = {MW2(2)
+ 4,(W,(2))
where 41, qS2are single Slater-type orbital (STO) 1s functions, must be lower than the corresponding closed-shell wavefunction. However, in later work on the excited 3S and 2% states of helium, Sharma and Coulson [2] noted Correspondence to: G. Doggett, Department of Chemistry, University of York, York YO15DD, UK. *Dedicated to the memory of Professor Charles A. Coulson.
0166-1280/92/$05.00 0 1992 Elsevier Science Publishers
B.V. All rights reserved.
314
that the constraint of orthogonalising 41 and +2 in the 2% state leads to a loss of variational flexibility; but no such loss of quality occurs for the triplet state, as the single determinant spin eigenfunction is invariant to linear transformations on the space spanned by $1 and &. The work described in this paper is concerned with the application of the simplest example of the spin coupled model to determine the optimum orbitals for the lowest excited singlet and triplet S states. The nonorthogonality of C$J~ and & presents no problem for the singlet state but, as shown in the next section, it is necessary to use a special partition of the basis set in order to avoid the numerical instability which arises from the use of non-orthogonal orbitals r& and & for the triplet state. THEORY
The spatial part of the 93 ground or excited states may be used to form the usual spin eigenfunction: yll = d4,(W,(2)
* 3w9 2)
where S,,(l, 2) = 1/,/2[41)/?(2) - 42)/?(l)]. In the method used for determining the optimum form of the orbitals +1 and &, an expansion is made over an even-tempered basis of 1s STOs: +j
=
it,
cij5i
where the exponents for the (& > are given by hi = e/fiN-i ;
i = N, (Iv - l), . . . ) 1
/? is taken as 1.5, and a is a scale factor. Basis sets of this kind are referred to as UN, where N is the number of primitive functions, ii. In determining $1, the coefficients (cia: i = 1, . . . , N) are fixed (to define the potential experienced by the electron occupying &). ‘I’ is then expanded in the form
i=l
and the cil are determined by solving the variational problem Hc, = SC~E, where
S,, = j~L(Wz
. S&W,
- S&,%
The new c, defines the potential for the occupation of &, and an analogous
315
equation for determining c2 is obtained. The whole process is cycled to self-consistency. The scheme described here has some similarity with the procedure described by Davidson [3]; but it should be noted that we work in the mixed c, 4 basis-the computational advantages of which become apparent when dealing with three or more electrons, as only those integrals involving the reformed orbital need be updated in each self-consistent-field (SCF) cycle. In addition, unlike conventional SCF procedures based on the use of a linear combination of atomic orbitals (LCAO) expansion, the Neigenvalues obtained at each cycle correspond to total energies; furthermore, two stacks of such energy states are produced: one associated with the optimisation of & and the other with & (stacks 1 and 2, respectively). For the 2% state, for example, the lowest and second lowest roots derived from stacks 1 and 2, respectively, are selected at each cycle of the iterative procedure. Further details of the stacks, and the relation of the method to the more conventional SCF model, based on the use of effective one-electron hamiltonians, will be presented elsewhere [4]. The eigenvalue obtained for the 29 state is an upper bound to the exact energy since, although the wavefunction itself is not orthogonal to the optimum split-shell ground state wavefunction derived using the U6 basis, the methodology [4] ensures that it is orthogonal to an underlying non-optimum split-shell wavefunction. The situation for the 3S state, obtained using a U6 basis, is more interesting: the procedure outlined above for determining the optimum orbitals 41, & when using the spatial and spin functions
and $092)
= 4042)
respectively, leads to numerical instability as a result of the coalescence of & and &. However, as noted above, Y’, = d$,(l)c#3(2)S,(l, 2) is invariant to linear transformations in the space spanned by 41, &,; hence, following Sharma and Coulson [2], $1 may be orthogonalised to & without loss of generality-as is the custom when forming the unrestricted Hartree-Fock (UHF) wavefunction. Clearly, orthogonalisation is computationally desirable in some respects but, within the conventional approach, a Lagrangian multiplier has to be introduced to ensure that orthogonality is maintained during the iteration procedure-thereby avoiding the possible coalescence of 41 and &. An alternative way of avoiding possible numerical instabilities is to select basis functions for 41 and & in such a way that the respective orbital subspaces remain distinct. Davidson [5] achieved this requirement (for both the 2% and 3S states) by expanding $1 and & over completely disjoint basis
316
sets containing four and ten basis functions, respectively. Clearly, he is performing a 14 basis function calculation with reduced quality, as the maximum variational flexibility cannot be achieved with this particular partition-simply because the best wavefunction for this kind of basis can be achieved only at the single configuration level by expanding 41 and & each over the 14 basis functions, together with the loss of one degree of freedom from each orbital as a consequence of orthogonalisation. We have developed another approach to this problem of avoiding numerical instability when seeking self-consistent solutions for states with total spin S = n/2, where n is the number of electrons. The rationale for this approach arises from our work on spin coupled wavefunctions in which the orbitals +1 and & are allowed to be overlapping in the triplet state (for full variational flexibility this must also be allowed for in singlet states as discussed above). A degree of overlap can be permitted by selecting basis sets for $1 and & that are minimally disjoint in the case of two electrons: in this case, three basis orbital subsets B, Q?,and 9 can be constructed in which 99 and 9 have just one element (the choice of B and 9 is not unique). 41 and & are then expanded over the basis sets defined by $3 u %?and B u 9, respectively. All partitions of this form, involving the assignment of a single ci to both 93 and 9, will suffice as they all give rise to a Y’, with the same energy; however, different LCAO coefficients are obtained for each &. In schematic form we have [.
[-;--;-..,I
I L______--__--___,
I
which is why we refer to such basis sets as “slide-rule” basis sets. Increasing the number of basis functions included in W and/or 9, with comcomitant reduction in the size of %‘,yields a wavefunction with a variationally poorer energy. In the limiting case, as is seen in the work of Davidson [5], % is the null set and 93 and $3 contain four and ten elements, respectively-a rather drastic construction. The number of cycles necessary to achieve convergence to ten decimal places for the triplet state energy is usually between four and six; however, for particular choices of both scale factor and the subsets 93 and 9, convergence may become progressively more difficult. For example, in the case of Li’ , with 93 = {C,} and 9 = {&}, nine cycles are required for a = 5.70; but, with a incrementing in steps of 0.05, and using the self-consistent set of LCAO coefficients from the previous calculation, the number of cycles rises steadily until at a = 6.65,140 cycles of the iteration procedure fails to yield
317
convergence. At, or before, this juncture it is just necessary to redefine g = (c3), and convergence is achieved in four cycles. This problem arises periodically, and is associated with the nature of the even-tempered basis: examination of the LCAO coefficients for both orbitals shows that the same basis function has become dominant in the LCAO expansion of both atomic orbitals-thus indicating that the spatial functions are trying to coalesce, i.e. the overlap integral is too large. The remedy is simply to remove this dominant basis function from one of the orbitals by redefining 99 as described above. At the end of the calculation, the overlap integral between & and & can be obtained, but its value will show slight variations according to the partition of the basis. The results obtained with the slide-rule basis are identical with those obtained by expanding & in the subspace which forms the orthogonal complement of 41 (and vice versa), thus confirming that the constraint of orthogonalisation for the triplet state can be simulated by the use of a slide-rule basis, thereby enabling the spin coupled method to be extended to high spin states in a simple way. In contradistinction to the calculation of the triplet energies, as a function of the scale factor, the corresponding calculations of the 2’S state energies are achieved in four cycles for the range of values of a which includes the energy minima. RESULTS
The total energy (Esx) for the lowest triplet state as a function of the scale factor a is given in Figs. l-3 for the helium-like species He, Li’ , and Be2+, when using a slide-rule basis of six even-tempered basis functions (dotted lines). For comparison, the energy obtained from the full configuration interaction (CI) calculations using this basis set yields the full line in each figure. It is clear that there are three optimum U6 basis sets; furthermore, the smallest optimum value of a is very close to the value of the nuclear charge but, for the three species studied here, this value is not associated with the global energy minimum. The other interesting feature of these plots is that, for all three species, the energy difference EsR - E,, is very nearly constant within the range of scale factors that contains the energy minima. It is also useful to examine the variation in the minimum energy of the 2% state as a function of scale factor, to investigate whether the optimum basis set for the triplet state provides an adequate description for the excited singlet state. In the case of helium, the global minimum occurs for the 2% state at a = 3.56--although local minima are also seen in the vicinity of a = 2.0 and a = 2.7. Similar plots for Li’ and Be2+ yield global minima for a = 4.40 (with two other minima in the vicinity of a = 3.2 and a = 6.8) and a = 6.44 (with a local minimum at a = 10.07), respectively; but,
318
-2.17416
-2.17420
-2.17424
-2~7426:,,,,,,,,,,,,,,,,,,,,,,,,,,, 2.00
2.66
3.06
3.66
4.00
4.60
Scale Factor, a Fig. 1. Plots of the total energy for the 152s ‘5 state of helium. as a function of the basis set scale factor (a): (. . . .) ESR; (-) EC,.
TABLE 1 Summary of the overall energy results for helium State
N
Scale
E (a.u.)
C&11&>
B
D
Ref.
3
14 6 6 6
4.06 4.06 4.06
-
2.17425 2.17425 2.17425 2.17425
0.1332 0.1047 0.0650 0.1459
5, . . .5, 52
1;6. . .L
5
:;
::
14 6
3.56
- 2.14344 - 2.14346
0.0821 0.0812
6, . . -5, null
null
2lS
I,
56.. .t-14
This work This work This work
5
This work
319
-5.10914
-5.10920
i
d
B
15 G-5.10926
-5.10632
-5.10638
I,,
I
3.00
I,,‘,
,,,,,,,,,,,,,,,,,,
4.00
5.00
6.06
7.00
6.00
ScaleFactor,a Fig. 2. Plots of the total energy for the 1~2s %3state of Li’ as a function of the basis set scale factor (a): (. ’ . .) ESR;(-) EC,.
in the vicinity of the global minimum for the triplet states, the appropriate values of a are all reduced in value--corresponding to a more dilated basis set in the singlet states of the three isoelectronic species studied here. These general features of the optimum basis sets are apparent in Fig. 4, which shows plots of the radial distribution function for the ZIS and 3S states for values of a at the respective global minimum (the plots for the values of a corresponding to other local energy minima are imperceptibly different from those given in Fig. 4): furthermore, it is encouraging that the results obtained by Winkler and Porter [6], using correlated wavefunctions of Hylleraas quality, yield radial distribution functions with forms identical in appearance to those calculated here.
320
-9.29905-
-9.29510-
-9.29515-
-9.2%20-
-9.29525-
. -92952e i 6 -9.2959s $ 15
-9.29540-
-9aw-
-929590-
-9.29555-
-9.29590-
.9.23565 ;,,, 3.00
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,’,,,,,,,,
4.00
5.00
6.00
7.00
9.W
9.w
10.00
ll.w
12.00
Scale Factor, a Fig. 3. Plots of the total energy for the 152s ‘S state of Be’+ as a function of the basis set scale factor (a): (. ...)ESR. (-)E,,.
The overall energy results for helium are summarised
in Table
I.
CONCLUSION
In this work, the differences in the electronic structure between the 2% and 3S excited states for helium-like species has been explored using a method in which the two orbitals are not constrained to be orthogonal. In the case of the singlet, the relaxation of the orthogonality constraint is a necessary requirement as indicated in Sharma and Coulson’s [2] analysis. In calculating the energy of the triplet state, it is shown how the constraint of orbital orthogonality in the UHF method can be circumvented
321 1.20 -I
l.OU-
om-
I
.*
2 +i om3 u
!z
0.W
1.00
2.W
3.00
4.00
5.00
6.00 I
7.00
6.W
9.w
10.00
11.w
12.00
1a.u.
Fig. 4. Plots of the radial distribution function for the 2% (-) obtained using a U6 basis with optimum scale factor (a).
and as (. . .) states of helium
through the use of slide-rule basis sets in which, by a suitable partition of the basis functions, it is possible to simulate the consequences of orthogonality without requiring the atomic orbitals to be orthogonal. Furthermore, as further calculations show [7], this basis set partitioning procedure has useful applications to both closed and open shell systems, in that it enables programs developed for the calculation of spin coupled wavefunctions to be used, with only minor modification, for calculating SCF and UHF wavefunctions without requiring the constraint of orbital orthogonality, thus providing a unification of the SCF molecular orbital and spin coupled approaches to the determination of electronic structures. More importantly, the method permits the calculation of excited state energies and wavefunctions without any “sagging” [2], with the consequent loss of the variational upper bound. Clearly, basis set partitions of the kind described here could be contrived by carrying out linear transformations on the orbitals 4; associated with a conventional SCF or UHF wavefunction; however, no useful insight is gained by this procedure and neither can excited state energies be calculated using a single configuration. It is useful to compare the calculated triplet energies, over the range of scale factors considered here, with the results of full CI calculations using the same U6 basis. The difference in the total energies for the slide-rule and
322
full CI calculations (Figs. l-3) is found to be very nearly independent of the
choice of scale factor in the region where the total energy exhibits a local minima: a result which, if maintained for other electronic states, would suggest that the calculation of excitation energies may be more tractable than the total energies themselves. Finally, the plots of the spinless radial distribution functions for the lowest excited singlet and triplet states (Fig. 4) show an overall contraction in the valence shell for the triplet state relative to the singlet state. Further examination of the effects of spin coupling on the underlying electronic properties for these two-electron species, together with threeand four-electron species, will be reported elsewhere. ACKNOWLEDGEMENT
G.D.F. acknowledges the award of a S.E.R.C. Postgraduate Studentship. REFERENCES 1 CA. Coulson, in P.-O. Liiwdin (Ed.), Q uantum Theory of Atoms and Molecules and the Solid State, Academic Press, New York, 1966, p. 97. 2 C.S. Sharma and C.A. Coulson, Proc. Phys. Sot., 66 (1962) 81. 3 E.R. Davidson, J. Chem. Phys., 41 (1964) 656. 4 G. Doggett and G.D. Fletcher, Phys. Rev. A, submitted for publication. 5 E.R. Davidson, J. Chem. Phys., 42 (1965) 4199. 6 P. Winkler and R.N. Porter, J. Chem. Phys., 61 (1974) 2038. 7 G. Doggett and G.D. Fletcher, unpublished work.