J . Phys. Chem. 1988, 92, 3029-3033
3029
Gaussian Basis Sets for High-Quality ab Initio Calculations Jan Almlof,* Trygve Helgaker, Supercomputer Institute, University of Minnesota, Minneapolis, Minnesota 55455
and Peter R. Taylor ELORET Institute, Sunnyvale, California 94087 (Received: August 7 , 1987)
Different types of Gaussian basis sets for accurate LCAO calculations are discussed. For calculations designed to recover a substantial portion of the correlation energy, we suggest the use of basis sets comprising natural orbitals from correlated calculations on the atoms. These basis sets have proven to be very efficient in accounting for large fractions of the molecular correlation energy. For cases in which an SCF or MCSCF treatment is adequate the use of a floating basis provides a rapid convergence to the large-basis limit.
Introduction In ab initio calculations, deviations from the exact nonrelativistic limit are caused by truncations of the one- and N-particle basis sets. It is desirable to improve these approximations in a systematic way toward a well-defined limit. Due to recent full CI benchmark calculations' the truncation effects on the N-particle basis (the configuration expansion) are now rather well understood. In contrast, the construction of basis sets for LCAO calculations is an area where development has been slow. In Schaefer's words,2 the selection of an appropriate basis set is still "more art than science". In correlated calculations, the size of the final contracted basis set often constitutes a bottleneck, since this determines the amount of work required to optimize the wave function. In contrast, the integral time is of minor concern. Clearly there is a need for compact basis sets even if it means a more expensive evaluation of integrals. A configuration basis can be gradually extended, at least in principle, from one configuration to the full C I limit. Extensions can be made by adding successively less important configurations (or classes of configurations) to the basis, and truncation effects can be explored in a systematic way in a C I calculation. Similarly, a perturbation treatment of correlation energy provides a welldefined path toward the exact solution. In contrast, most studies of one-particle basis set effects have employed rather arbitrarily chosen basis sets of different size. Variations in atomic and molecular properties are erratic and difficult to interpret, since a smaller basis set does not normally span a subspace of the larger one. The basis sets used for calculations on polyatomic molecules are almost invariably Gaussians centered on the positions of the atomic nuclei: +(l,m,n,a) = (x
- x,)'a(y
- y,)ma(z - za)nne-aa(r-ra)2 (1)
culation with an uncontracted basis of fixed size, optimizing the orbital exponents either with a numerical grid-based procedure or by using analytic derivative techniques.) The coefficients for the expansion of the self-consistent AOs are then used as contraction coefficients. Krishnan et aL6 optimized the exponents and contraction coefficients of their 6-31 1G** sets in MP2 calculations. These basis sets were designed specifically for correlated calculations, which is an obvious advantage, even though a variational treatment of MP2 energies raises some formal objections. Ideally, procedures such as those outlined above would lead to basis sets made up by the optimum atomic orbitals. An S C F calculation on the atom with this contracted basis would thus reproduce the results for the (usually much larger) primitive basis. However, this is seldom the case in practice. For reasons of convenience and tradition, segmented and nodeless basis sets are often used. In a segmented basis, a primitive Gaussian can contribute to only one contracted function. Clearly, for all atoms beyond He the nodal structure of the atomic orbitals is compromised in such a basis. (The popular small STO-nG basis sets are also segmented, but they do incorporate an approximate nodal structure in basis functions with n > 1 + 1. As a result, these basis sets sometimes produce molecular results which are competitive with those from energy-optimized, nodeless basis sets.) In general, the errors introduced by compromising the nodal structure are several orders of magnitude larger than those due to inadequacies in the primitive basis. As a remedy for the node problem, Raffenetti' suggested the use of a general contraction scheme, Le., one where each primitive function can contribute to all the contracted ones. This, of course, would make it possible to use the "true" atomic S C F orbitals as a basis, which would in turn automatically provide the correct nodal structure. There are some problems with this scheme, however. One is that only a minimal basis set is defined. have been made to circumvent this problem by augmenting the basis set with virtual orbitals, using either an N o r an N - 1 electron potential, but the results were not always encouraging. Another objection is, of course, that the entire procedure is biased toward SCF calculations on atoms.
Variational optimization of Gaussian basis sets is normally restricted to orbital exponents. Recent studies3 have shown that although many basis sets in present use are incompletely optimized, those deficiencies have a minimal effect on computed atomic and molecular properties. Attempts to improve basis sets further must therefore focus on other parameters of the basis, Le., contraction schemes and/or basis set centers. Basis sets are usually optimized for the free atom (see, however, the recent work by Schlege14on exponent optimization in molecules). The conventional approachS is to carry out an SCF cal-
Polarization Functions A basis set that provides a good description of the atoms is a necessary, though by no means sufficient, condition for balanced calculations on molecules. In addition to describing the atoms well, the basis must have the necessary flexibility to account for the deformation that the atoms undergo when a molecule is
( 1 ) Bauschlicher, C. W., Taylor, P. R. J . Chem. Phys. 1986,85,2779, and references cited therein. (2) Rothenberg, S.;Schaefer, H. F., 111 J . Chem. Phys. 1971, 54, 2764. (3) Almlof, J.; Faegri, K., Jr. J . Comput. Chem. 1986, 7, 396. (4) Tonachini, G.; Schlegel, H. B. J . Chem. Phys. 1987, 87, 514.
(5) Huzinaga, S. J . Chem. Phys. 1965, 42, 1293. (6) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J . Chem. Phys. 1980, 72, 650. (7) Raffenetti, R. C. J . Chem. Phys. 1973, 58, 4452. (8) Feller, D. F.; Ruedenberg, K. Theor. Chim. Acta 1979, 52, 231.
0022-3654/88/2092-3029$01.50/0
0 1988 American Chemical Society
3030 The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 formed. For the calculation of molecular properties, the basis must also be able to account for the change of the wave function caused by an external perturbation. The required flexibility is achieved by augmenting the basis set with diffuse functions and polarization functions. Polarization functions are basis functions with angular momenta higher than those occupied in an atomic SCF calculation. Their purpose is to change the shape of the atomic orbitals from that adopted in the atom and to account for a shift of the center of gravity of the charge density from the atomic nucleus. Polarization functions have their radial maximum in the same region as the valence orbitals. The center of the basis function (ra in eq 1 ) is usually chosen to coincide with the atomic nucleus. For calculations on atoms this is a natural choice, but for molecules where the local environment of the atom has a lower symmetry, it may be expected that this position is not necessarily the optimal one. It is easily shown that the derivative of a Gaussian with respect to its center is a combination of Gaussians with the same exponent but with higher and lower I-values: d$(l,m,n,a)/dx, = l$(l-l,m,n,a) - 2cu$(l+l,m,n,cu) ( 2 ) Therefore, if the polarization functions have exponents similar to those of valence orbitals, they will describe the effect of displacing the basis functions from the nuclei-at least to first order in a perturbation or Taylor expansion. An alternative approach-one which will be discussed later-would of course be to allow the basis function to actually move off the nucleus.
Correlation Functions In order to describe angular correlation effects, basis functions with high angular momenta are necessary. Like polarization functions they need to have their radial maxima in the region of the valence electrons, although the argument for identical orbital exponents no longer applies. Often these correlation effects demand a larger and more flexible basis set than does polarization. Polarization effects are therefore usually well accounted for in beyond-Hartree-Fock calculations with a basis set adequate to describe correlation effects. Basis Set Superposition Errors For most chemical applications relative energies are more important than absolute ones, and for all types of bonding it is important that any deficiencies which may exist in the description of the monomer fragments are not accidentally removed when the bond is formed. Unfortunately, this is not always the case. Deficiencies in the monomer description are often reduced by basis functions on the other fragment. This effect, which may occur on both the one- and N-electron level, is referred to as basis set superposition error. This undesired "improvement" is in general minimized if the description of the monomer fragments is as good as possible, while keeping the basis sets small.
A N 0 Basis Sets As mentioned above, several objections can be raised against conventional basis sets: they are segmented, which compromises the nodal structure; the functions are centered on the atomic nuclei, which increases the need for polarization functions; and the contraction scheme is usually based on SCF calculations with no consideration of electron correlation (which is what ultimately places the most stringent requirement on basis sets). For highly correlated calculations, we suggest the use of a different type of basis set: the atomic natural orbital (ANO) basis,9 a basis made up by the natural orbitals with the highest occupation numbers from an atomic CI calculation in the uncontracted basis. Although it is not required for this scheme, we suggest that the uncontracted basis be chosen large enough to saturate each I-value space for practical purposes. With modern codes the time for integral evaluation is a minor part of the total time required for a high-level correlation calculation, and with a near complete set (9) Almlof, J.; Taylor, P. R. J . Chem. Phys. 1987, 86. 4070.
Almlof et al. of primitives the ANOs can adopt a more optimal form and fewer of them are needed. Such a basis has numerous advantages: (i) It accounts for electron correlation in a way that is very nearly optimal for the atom with a given number of basis functions. (ii) It has the correct nodal structure. Since the atomic fragments are so well described and the basis sets are very compact, basis set superposition errors are minimized. (iii) It can be extended in a smooth, well-defined, and physically reasonable way toward the uncontracted basis limit by lowering the occupation number cutoff. With each such extension the larger basis set continues to contain the smaller one as a subset, and the change in computed properties can be ascribed entirely to the new functions added. This is a great advantage and is in strong contrast to most basis set investigations performed with conventional basis sets, in which smaller sets seldom survive as a subset of the larger sets used. The development of A N 0 basis sets requires methods and codes that treat atoms on a level similar to the intended molecular calculation (SDCI, MRCI, etc.) while retaining full atomic symmetry in the final natural orbitals. It has recently been shown how this can be achieved with little extra computational effortlo even for very complicated atomic states and sophisticated treatment of electron correlation. Finally, of course, a general contraction scheme is needed, Le., one that allows for a primitive basis function to contribute to several contracted ones. This should not be a problem with modern codes for integral evaluation. For efficiency reasons these codes treat all integrals involving primitive functions of four different shells in parallel, in order to obtain sufficient vector lengths for modern supercomputers. This naturally produces the integrals needed for a general contraction scheme in the appropriate grouping.
More Complicated Situations It is evidence from the above that an optimal description of the atom is crucial for the A N 0 strategy. The choice of atomic state is one question which naturally arises. In many cases, the atomic ground state is representative for the situation in the molecule, and in other cases the difference seems to be insignificant. In the methylene molecule, for example, the singlet-triplet splitting is unaffected to within 0.1 kcal/mol whether an s2p2 or an sp3 optimized basis for carbon is used. There are cases, however, where such a simple approach is likely to fail. In transition-metal chemistry, for instance, it is often necessary to have a balanced description of s2dn,sld"+', and sodn+*states in the same calculation. The radial d-electron distributions for these atomic configurations are known to be quite different, and it is not likely that an A N 0 basis optimized for one particular state would have the flexibility to describe the others on an equal footing. Similar situations arise in calculations of electron affinities, where neutral and ionic states require a balanced description. For these problems, the A N 0 approach offers an excellent solution. A composite density matrix containing all the relevant states can easily be computed, either from a variational calculation on an average of atomic states or by averaging density matrices for different pure states. The eigenvectors of such a composite density matrix would constitute the appropriate A N 0 basis molecular calculations. It is observed that when such an approach is used, the A N 0 spectrum is more dense than usual and more basis functions must be included before the occupation numbers reach the desired threshold. This, of course, reflects the appropriate response to the increased demands placed on the basis set in such a complex situation. Results with A N 0 Basis Sets A N 0 basis sets have been applied and tested for a large number of different situations. The general finding is that they perform very well for energies and closely related properties, in fact, better than STO sets of similar size for the cases where comparisons have been made. For situations where a flexible description of the diffuse regions is more important than an optimal energy;it is necessary to augment the set. This can be done by allowing the ( I 0 ) Bauschlicher, C . W., Jr.; Taylor, P. R. Theor. Chim.Acta, in press.
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3031
Gaussian Basis Sets for ab Initio Calculations TABLE I: Contraction Error in phartrees on the Total Energy of HBr for Different Contraction Schemes“ AESCF
a b C
d
633 543 229 696938
AEcor
1214 3446 496 4061
AESCF
e f g
1361 631 536
AEcor
4963 1629 3596
OThe uncontracted basis set for Br is Dunning’s” 14sllp6d augmented with f and diffuse p and d to 14s13p7d4f. For H an uncontracted 8s6p4d basis was used. The uncontracted energies are EsCF= -2572.981 707 1 hartree and E,, = -0.190724hartree. (a) Contracted A N 0 [7 6 4 2/4 3 21. (b) Contracted A N 0 [6+1,5+1, 3+1, 1 + 1 / 3+1, 2+1, 1+1]. (c) Contracted A N 0 [7+1,6+1, 4+1, 2+1/4+1, 3+1, 2+1]. (d) Contracted segmented [7 6 4 2/4 3 21 (7+2+1+1+1+1+ 1, 7+2+1+1+ 1+ 1, 4+1+1+1, 3+1/5+1+1+1, 4+1+1, 3+1; contraction coefficients for f orbitals are based on natural orbitals). (e) Contracted Raffenetti [7 6 4 2/4 3 21. (f) Contracted A N 0 Br/Br- [7 6 4 2/4 3 21. (9) Contracted A N 0 Br/Br- [6+1, 5+1, 3+1, 1+1/3+1, 2+1, 1+1]. TABLE I 1 Contraction Errors (au X 10”) in the Dipole Moment of HBP
a b C
d
SCF
cor
78 151 110 60
128 -198 -23 -172
e f g
SCF
cor
63 185 145
-174 62 -176
“The uncontracted values are wsCF = 0.373687 au and ApLcor= -0.040 390 au. The notation for different contraction schemes is as in Table I. most diffuse Gaussians to enter the basis separately and/or adding more diffuse functions. Some results from applying these basis sets to the HBr molecule e; presented in Tables I and 11. The calculations were carried out at the experimental equilibrium geometry (1.4144 A). To get reference results for a comparison, SCF and CI calculations were carried out with the uncontracted basis set. The C I was performed at the SDCI level, using the SCF orbitals and correlating the eight valence electrons. The basis was then contracted, testing a number of different contraction schemes. Table I shows that all the ANO-based contractions perform well, with a total loss of less than 1 mhartree due to contraction. As expected, the segmented basis performs very poorly in this respect. The A 0 coefficients do not readily suggest a unique segmented contraction scheme, and indeed Dunning” questioned whether a reasonable segmented contraction scheme could be found for this set (and others of the same size). The one used here has an energy loss of nearly 0.7 hartree. This is due to the incorrect nodal structure of the AOs enforced by the segmented contraction. Interestingly enough, the “Raffenetti”’ type contraction based on atomic SCF orbitals is not competitive with ANOs for HBr even at the S C F level. The reason is that this scheme defines only the strongly occupied space, i.e., a minimal basis set. Additional functions have to be chosen on rather arbitrary criteria-in this case the outer primitive shells were left uncontracted. At the S C F level, the main advantage of using ANOs is the correct nodal structure of the AOs and the unique definition of the functions outside the minimal basis. For the correlation energy another interesting observation can be made. First, it should be noted that the error in correlation energy is larger than the S C F error, so that despite the fact that the contractions are optimized for correlation the higher angular momentum sets are perfectly able to describe molecular polarization effects. The contraction loss in the correlation energy for the A N 0 sets obtained by averaging Br and Br- ANOs is considerably larger than that for the Br A N 0 sets, although the S C F contraction errors are actually slightly smaller. In Table I11 the deviations from the uncontracted results for the HBr dipole moment are listed. As the total dipole (11) Dunning, T.H., Jr. J . Chem. Phys. 1977, 66, 1382.
TABLE 111: Energy Recovered (Percent of Difference between DZ and the Hartree-Fock Limit”)
fixed floating
DZ
DZ+
DZP
DZP+
0.0 27.3
8.4 36.8
51.8 56.5
59.5 64.6
“Calculated using 6-31l++G(3df,3pd) basis moment is the sum of two terms which differ in sign, where the S C F and correlation contraction losses are opposite in sign the total dipole will be in better agreement with the uncontracted result than when the signs are the same. Consequently, the total dipole moment (at the CI level) is given best by the A N 0 contractions in which the outermost primitives have been freed. However, there is little in the way of a systematic trend in the individual results: the S C F contraction error is smallest for the SCF-derived contractions (segmented and Raffenetti), but these have the largest contraction errors at the correlated level of all except the simple A N 0 set. It should be noted, however, that all errors are small, none being larger than 0.002 au. Hence if it is desired to obtain good correlation energies and good dipole moments, any of the augmented A N 0 sets would be a suitable choice.
Floating Basis Sets As noted above, atom-centered basis sets require augmentation with polarization functions to describe the deformation in atomic densities attendant on formation of molecules, and this is particularly important for properties. For example, when an external electric field is applied to a molecule, the electrons and nuclei will be displaced in opposite directions, and in an atom-centered basis this effect can only be accounted for if polarization functions are included. An alternative to polarization functions is to allow the orbitals to move off the nuclei and occupy their optimal positions in the molecule. This adds considerable flexibility to the basis set without increasing its size. Due to the high power dependence of electronic structure calculations on the basis set size, such “floating orbital” basis sets have considerable advantages. On the other hand, the optimization of the orbital positions is clearly a nontrivial task, being at least as expensive as a conventional geometry optimization. A final judgment on the merits of floating orbitals is therefore difficult, especially as new technical developments change computational efficiencies. Interest in floating functions is, of course, not new. Hurley12 investigated floating functions in the 1950’s, Frost introduced floating spherical Gaussian orbitals 20 years ago,13 and more recently HuberI4 has carried out extensive studies using his floating orbital geometry optimization (FOGO) model. However, with recent developments in second-order techniques for optimization of geometries and calculation of molecular properties as energy derivatives, the use of floating orbitals has become more attractive. There are two main steps in floating orbital calculations: optimization of the wave function and the calculation of properties. The optimization may be carried out as a conventional geometry optimization by treating the orbitals as charge-free atoms. With a second-order convergence scheme (the trust-region method) the orbital positions usually converge (to within lo-’ in the gradient) in less than five iterations, even though the initial Hessian often has several negative eigenvalues. The analytical calculation of properties can also be carried out in a straightforward manner by treating the floating orbitals as charge-free atoms. However, when second-order properties are calculated, an additional (inexpensive) computational step is required at the end of the calculation to account for the relaxation of the orbital positions after the perturbation. For example, if we calculate the perpendicular component of the polarizability of H 2 0 using a fixed double-{ basis set, we obtain 1.29 au, considerably less than the Hartree-Fock value, which is close to 8. (12) Hurley, A. C . Proc. R . SOC.London, A 1954, 226, 170. (13) Frost, A. A. J . Chem. Phys. 1967, 47, 3707. (14) Huber, H . Theor. Chim. Acta 1980, 55, 117.
3032 The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 TABLE I V Errors in Calculated Properties (%) Relative to the Hartree-Fock LimiP ~
DZ
DZ+
DZP
DZP+
1.5 0.8
1.6 0.9
0.3 0.2
0.3 0.3
4.4
5.0 0.6
0.6 0.6
0.2 0.2
4.4 I .7
4.5
I .6
0.9 0.6
0.7 0.4
21.3 3.5
28.5 3.1
7.3 4.0
9.1 I .0
35.8 7.1
21.2 I .4
28.4 7.1
15.1 1.2
69.7 22.5
73.3 9.4
12.4 15.6
19.0 2.9
3.4 1 .o
3.7 I .o
0.6 0.5
0.4 0.3
42.2 11.0
41.0 4.6
16.0
14.4
8.9
1.7
22.8 6.0
22.3 2.8
8.3 4.7
7.4
Almlof et al. TABLE V. Properties Calculated for H 2 0 DZ DZ+ DZP ‘OH,
bond distances
polarizabilities, au in-plane i in-plane
11
polarizabilities intensities
out-of-plane frequencies, cm-‘ a, str
geometrical properties
a, def b, str
electrical properties intensities, km/mol a , str
all properties a , def
I .o
b, str
”Calculated using 6-31 I++G(3df,3pd) basis. The results for floating orbitals are given in italics.
If we now optimize the orbital positions and recalculate the polarizability using the same procedure as for the fixed orbitals, we obtain 1.13-almost the same as before. The reason is that we have only taken into account the relaxation of the orbital centers in the absence of a field. An additional relaxation term is required to describe the response of the centers to a field. This “dynamic contribution” is easily calculated from the second derivatives with respect to displacements of the orbital coordinates and the dipole derivatives with respect to the orbital coordinates. We obtain 5.03 au for this term, which gives a total perpendicular polarizability of 6.16, considerably closer to the Hartree-Fock limit. The results of our calculations are summarized in Tables I11 and IV and are based on Hartree-Fock calculations on HF, H,O, NH3, CHI, CO, and H2C0. For each molecule we have calculated geometries, dipole moments, static polarizabilities, harmonic frequencies, and double-harmonic infrared intensities using four different floating basis sets. All calculations were carried out with the SIRIUS1’ and ABACUS'^ programs. The calculated numbers have been compared with those obtained with the corresponding fmed orbitals and also with results from high-quality Hartree-Fock calculations at the 6-3 11++G(3df,3pd) level. The complete results will be published elsewhere,” and we concentrate here on some general features. All properties have been calculated at the theoretical minima rather than the experimental geometries. In this way we compare results that have been obtained without any a priori knowledge about the system. Although calculations at the experimental geometry in many cases give values closer to observed properties, our purpose is to examine the convergence of the one-electron basis, and experimental values are therefore not of immediate interest. The double-( basis (DZ) is van Duijneveldt’s (8s4p/4s) set1* contracted to (4s2p/2s). DZP is the same basic set augmented with one set of polarization functions on all atoms, using the exponents of Krishnan et a].’ Finally the DZ+ and DZP+ basis (15) Jensen, H. J. Aa.;Agren, H. Chem. Phys. Lett. 1984,IIO.1 4 0 Chem. Phys. 1986,104, 229. (16) Helgaker, T. U.; Almlof, J.; Jensen. H. J. Aa.; Jsrgensen, P. J . Chem. Phys. 1986,84, 6266. Helgaker, T. U.; Jensen, H. J. Aa.; Jsrgensen, P. J . Chem. Phys. 1986,84, 6280. (17) Helgaker, T.; Almlof, J., to be published. (18) van Duijneveldt, F. B. IBM Res. Rep. 1971,RJ 945
0.949 0.947
0.950 0.948
0.941 0.942
0.941 0.942
0.941
111.9 106.1
112.9 106.6
105.6 105.5
106.3 106.2
106.2
0.977 0.735
1.016 0.785
0.835 0.731
0.874 0.777
0.774
6.56 8.27 3.70 7.09 1.29 6.16
7.34 8.57 4.87 7.71 3.86 7.06
6.74 8.12 4.90 7.08 2.58 6.10
7.26 8.41 5.78 7.78 5.22 7.26
8.61
3986 4074 1728 1799 4150 4178
3983 4055 1695 1777 4145 4159
4149 4134 1772 1779 4259 4218
4155 4124
4128
1733 1749 4268 4206
1747
2.6 6.2 124.5 96.3 50.6 58.0
3.7 12.1 148.0 96.2
18.6 9.8 87.0 97.6 52.7 60.9
28.7 15.3 90.0 98.9 91.1 89.3
dipole moment, au
frequencies dipole moments
b
LHOH, deg
bond angles
0.4
DZP+
A
81.1 84.5
7.81 7.39
4228 14.5 92.4
88.8
” T h e results for floating orbitals are given in italics. b6-311++G(3df,3pd).
sets were obtained by adding one set of diffuse orbitals to each atom (s and p on heavy atoms and s on hydrogen). The four floating basis sets are identical with the fixed basis sets except that the orbital centers are variationally optimized. For reasons explained elsewhere” we did not float the orbitals independently: the two innermost s functions on the heavy atoms were fixed at the nuclei, while the two outermost s functions were constrained to have the same (floating) coordinates. The diffuse functions in DZ+ and DZP+ are attached to these same floating centers. As a benchmark for the smaller floating and fixed calculations, we have carried out calculations using the 6-3 1 l++G(3df,3pd) basis set of Frisch et al.I9 We do not claim that this basis gives the Hartree-Fock limit for all properties considered, but it should be large enough to represent a limit on what can be obtained in routine calculations. Also, although more accurate calculations have been carried out for selected properties using larger basis sets, most calculations of near-Hartree-Fock quality have been performed at the experimental geometry. Since many properties are strongly geometry dependent, any comparison with such results would be misleading. For example, using the floating DZP+ and the large 6-31 l++G(3df,3pd) basis, we obtained 0.060 and 0.058 au for the dipole moment of CO (with carbon positively charged). From the calculated dipole derivatives we may further estimate that at the experimental geometries we would obtain 0.104 and 0.107 au, respectively. This in good agreement with numerical Hartree-Fock, which gives 0.104 at the experimental geometry.*O One way to measure the quality of a basis set is to compare its energy with that obtained at the Hartree-Fock limit. In Table I11 we have listed how much of the difference between the Hartree-Fock energy and the DZ energy (the “residual energy”) is recovered by each basis set. It is clear from Table I11 that floating functions cannot compete with the addition of polarization functions if the only criterion is energy lowering. The errors in the calculated properties are listed in Table IV. The properties fall into two classes: geometrical (geometries and ~
~
~
( I 9) Frisch, M J , Pople, J A , Binkley, J S J Chem Phys 1984,80,
3265 ( 2 0 ) Sundholm, D , Pyykko, P : Laaksonen, L Mol Phys 1985,56,1411
J. Phys. Chem. 1988, 92, 3033-3036 frequencies) and electrical (dipole moments, polarizabilities, and infrared intensities). The geometrical properties are well described a t all levels except fixed D Z and DZ+. The improvements obtained by floating are small, particularly for bond distances. Floating DZ gives results slightly inferior to fixed DZP but is still an improvement on the fixed D Z results. The situation is different for electrical properties, which are poorly described at the fixed-orbital level: the best agreement is obtained for the DZP basis dipole moment values, which differ from the Hartree-Fock limit values by some 7%. Floating produces drastic improvements in all cases. Indeed, for all the properties considered even the smallest floating basis set gives properties which are similar to or better than the best fixed basis set. Dipole moments at the fixed D Z level are in error by about 20%, compared to 3-4% when orbital floating is allowed. With addition of diffuse functions the computed values are somewhat larger than the Hartree-Fock limit, while the inclusion of both diffuse orbitals and polarization functions (floating DZP+) gives agreement with the Hartree-Fock limit to within 1%. Since experimental dipole moments are usually a few percent smaller than the Hartree-Fock limit, a floating D Z basis often gives dipole moments which are in fortuitously good agreement with experiment, as noted by Huber.21 The requirements for a fixed-orbital basis for calculating polarizabilities are quite demanding, and it is not surprising that even the largest fixed basis set (DZP+) produces polarizabilities which are 15% too low. When floating is allowed, the errors are 7% or smaller in all cases. In particular, when diffuse functions are added, the agreement with the Hartree-Fock limit is remarkably good. (21) Huber, H. J . Mol. Struct. 1981, 76, 277.
3033
Perhaps the most difficult to calculate of these properties are infrared intensities, and fixed D Z and DZ+ calculations give intensities which are completely unreliable. Again, floating results in significant improvements, the best basis set giving errors about 3%. It is clear that floating orbitals have a distinct advantage over fixed for the calculation of electrical properties. Even a modest basis set such as floating DZ+ gives errors less than 5% for electrical properties, while the floating DZP+ results are only 2% in error. The optimization of floating orbitals is little more time-consuming than the optimization of geometries. With second-order techniques the number of iterations does not increase significantly, normally by only one or two iterations. The time spent in each iteration is also not significantly increased. The time-consuming step is the calculation of second-derivative integrals, and this takes the same time for floating and fixed orbitals. However, if the solution of the response equations is carried out in an iterative fashion, the number of iterations will necessarily increase with the increased number of perturbations. Nonetheless, with an appropriately adapted code the total time for calculation of wave functions and properties using floating orbitals should not exceed twice the time needed for the corresponding fixed basis (provided a geometry optimization is carried out). Given the dramatic improvements in calculated properties, floating orbitals appear more attractive than previously thought and considerably more attractive than extending fixed basis sets with the polarization functions that would otherwise be required. Acknowledgment. This work was supported by the National Science Foundation, Grant CHE-8610809 (J.A.), the Institute of Mathematics and its Applications (T.H.), and by NASA Grant N C C 2-371 (P.R.T.).
An Open-Shell Spin-Restricted Coupled Cluster Method: Application to Ionization Potentials in N, Magnus Rittby and Rodney J. Bartlett*' Quantum Theory Project, University of Florida, Gainesville, Florida 3261 I (Received: July 15, 1987;
In Final Form: October I , 1987)
In order to circumvent the problem of spin contamination in unrestricted Hartree-Fock based coupled cluster (CC) calculations, we present a new method of calculation for certain classes of open-shell systems. The approach ensures that the proper spin component of the resulting correlated wave function is projected out in the energy evaluation by the use of a reference function constructed from suitably chosen restricted open-shell Hartree-Fock or other orbitals. This single-reference open-shell spin-restricted CC method is applied to the calculation of ionization potentials in the N2molecule, and it is shown that highly accurate results can be obtained in a 5s4pld basis. The mean error for all the principal ionization potentials of N2 compared to experiment is 0.45%.
Introduction In single-determinant reference coupled cluster (CC) 1-5 and many-body perturbation theory (MBPT) method^,^,^,' the use of an unrestricted Hartree-Fock (UHF) wave function for open-shell systems is frequently e m p l ~ y e d . ~A- ~U H F reference generally leads to a lowering of the total energy at the Hartree-Fock (HF) level since the a-spin electrons are chosen to occupy different spatial orbitals than the @-spinelectrons. This procedure introduces important spin-polarization effects on the spin density, for example, while the corresponding open-shell R H F results would erroneously give zero.8 Also, U H F functions facilitate the correct separation of a single-determinant wave function in bond breaking.2.3 However, there are two prices paid for the additional
'Guggenheim Fellow. 0022-3654/88/2092-3033$01 SO10
flexibility introduced by alU-HF function. First, the wave function is not an eigenfunction of SBS;~ and second, having to treat different (1) (a) Cizek, J. Ado. Chem. Phys. 1969, 14, 35. (b) Paldus, J.; Cizek, J. In Energy Structure and Reactivity; Smith, D. W., McRae, W., Eds.;
Wiley: New York, 1973. (2) (a) Bartlett, R. J.; Purvis, G . D. In?.J . Quantum Chem. 1978, 14, 561. (b) Bartlett, R. J.; Purvis, G. D. Phys. Scr. 1980, 21, 255. (3) Pople, J. A,; Krishnan, R.; Schlegel, H. B.; Binkley, J. S. Int. J . Quantum Chem. 1978, 14, 545. (4) Bartlett, R. J. Annu. Rev. Phys. Chem. 1981, 32, 359 and references therein. ( 5 ) Bartlett, R. J.; Dykstra, C. E.; Paldus, J. In Advanced Theories and Computational Approaches to the Electronic Structure of Molecules; Dykstra, C. E., Ed.; Reidel: Dordrecht, 1984; p 127. (6) Kelly, H. P. Ado. Chem. Phys. 1969, 14, 129. (7) (a) Bartlett, R. J.; Silver, D. M. Int. J . Quantum Chem.,Symp. 1974, 8, 271. (b) Bartlett, R. J.; Shavitt, I. Chem. Phys. Lett. 1977, 50, 190.
0 1988 American Chemical Society
