J. Phys. Chem. 1996, 100, 10131-10134
10131
Combining Multiconfigurational Wave Functions with Density Functional Estimates of Dynamic Electron Correlation Nathaniel O. J. Malcolm and Joseph J. W. McDouall* Department of Chemistry, UniVersity of Manchester, Manchester M13 9PL, U.K. ReceiVed: February 19, 1996; In Final Form: April 16, 1996
X
A reliable technique for combining multiconfigurational wave functions with density functionals would provide an economic route to the correlation problem for systems that are not well described by a single configuration. However, such procedures are prone to overestimate (overcount) the correlation energy since there is no simple way to separate the dynamic and nondynamic components of the correlation. We propose a method, based on partitioning the electrons of a system into two sets and applying a wave function description to one set while obtaining the correlation of the second set and the interaction between sets from a density functional. We test our method on the electronically challenging case of trifluoride anion and a series of diatomic molecules.
Introduction Chemical systems whose electronic states are poorly described by a single orbital configuration are usually treated by multiconfigurational techniques, for example, multiconfigurational self-consistent field theory (MCSCF) or the complete active space (CAS) variant of it (CASSCF).1 In these models, the one-electron molecular orbitals and the mixing coefficients between orbital configurations are simultaneously optimized. With modern computer programs, the size of the basis set employed in the molecular orbital expansion and the number of orbital configurations that may be treated are sufficient that quite large chemical problems may be tackled efficiently.2,3 While such techniques are known to treat the nondynamic electron correlation (proper dissociation behavior, near-degeneracy, orbital relaxation) effectively, they are inappropriate for the dynamic correlation problem.4 The latter, which is essential for the accurate prediction of bond energies, activation barriers, and molecular properties, is better obtained from multireference configuration interaction (MRCI),5 perturbation theory (MRPT),6-8 or coupled cluster (MRCC)9,10 calculations. The problem arises that the realm of applicability of MRCI, MRPT, and MRCC techniques is limited to that of very small molecules. Hence, while MCSCF calculations may be carried out on reasonably large molecules, an efficient treatment of the dynamic correlation is only possible for the smallest systems. To remedy this situation a computationally inexpensive approach to the dynamic correlation problem is required. Given the current interest in density functional methods, it is natural to look for ways of combining multiconfigurational wave function based techniques with a density functional estimate of the dynamic correlation energy. In density functional theory the correlation energy, EC, is given by
EC(F) ) ∫FC(F) dV
(1)
where F is the one-electron density and C(F) is the correlation energy per particle. C(F) is usually a complicated function of F and possibly ∇F.11 For Hartree-Fock wave functions we can correct the energy, ESCF, according to12
ESCF+Corr ) ESCF + Ec(FSCF)
(2)
where FSCF is the Hartree-Fock density. We could apply the X
Abstract published in AdVance ACS Abstracts, June 1, 1996.
S0022-3654(96)00489-3 CCC: $12.00
same approach for a multiconfigurational wave function based on the multiconfigurational one-electron density, FMCSCF, i.e.
EMCSCF+Corr ) EMCSCF + EC(FMCSCF)
(3)
Since EC can be evaluated at less cost than EMCSCF, this would ostensibly provide an efficient solution. However, direct use of eq 3 can lead to an overestimation of the correlation energy which can often be numerically large. A multiconfigurational wave function already contains a component of the correlation energy. Hence, in evaluating EC(FMCSCF) we can in principle double-count the correlation contribution. Lie and Clementi13,14 proposed a way around this problem by scaling down FMCSCF. Since C(F) depends on F, it is not sensitive to the two-particle nature of the wave function and hence will overcount the correlation energy. It is possible to use an C that also depends on the two-electron density. This is the idea behind the ColleSalvetti functional,15 which yields very good correlation energies but is computationally expensive relative to the more widely used correlation functionals which depend only on F and ∇F. In CASSCF methods, given a good choice of the active space, the active space electrons are fully correlated among themselves (within their subset of orbitals). The inactive electrons are not correlated, nor is the interaction between the two-electron subspaces accounted for significantly. The main bulk of the dynamic correlation that is missing must come from these latter omissions. In the next section we propose a very simple scheme in which the active electrons are described by a wave function of the CASSCF type, while the remaining interactions are dealt with using a density correlation functional. Theory For the development in this section we shall partition the electrons of a system into two sets, which we shall term the core and valence sets. We use these terms in a rather general sense. Hence, the valence set may correspond to those electrons which are involved in bond rearrangement during a reaction (corresponding to the active/inactive partition of the CASSCF method) and thus may include only a subset of the constituent atomic valence shells. The partitioning is quite general. Given the partitioned electron space, consider a calculation of the correlation energy using, for example, a full CI-type wave function (Ψ). We can separate the total correlation energy into components due to correlating the core electrons alone, the © 1996 American Chemical Society
10132 J. Phys. Chem., Vol. 100, No. 24, 1996
Malcolm and McDouall
valence electrons alone, and a term corresponding to the correlation of the core electrons with the valence electrons
EC(Ψ) ) ECcore(Ψ) + ECvalence(Ψ) + ECcore/valence(Ψ)
method
(4)
could be obtained by exciting the core electrons into all Ecore C could be obavailable valence and virtual orbitals. Evalence C could be obtained by coupled tained similarly. Ecore/valence C excitations from the core and valence sets. In the density functional approach, the total correlation, EC(F), is obtained from a correlation functional. Again, EC(F) can be separated (at least formally) as
EC(F) ) ECcore(F) + ECvalence(F) + ECcore/valence(F)
(5)
When we combine a density functional with the energy obtained from a multiconfigurational wave function, we do not want to (F) in include the correlation due to the valence space Evalence C the functional. Writing the total density of our partitioned sets as
F ) Fcore + Fvalence
(6)
ECvalence(F) ≈ EC(Fvalence)
(7)
we can obtain
Hence, from eq 5 and eq 7
EC(F) - EC(Fvalence) ) ECcore(F) + ECcore/valence(F)
(8)
(Ψ) as estimated by any technique To this we may add Evalence C that correlates the valence electrons alone. Usually we shall (Ψ), but if a sufficiently flexible not have an exact Evalence C active space is used then a good answer should be obtained. Hence the total energy may be written
Etotal ) EMCSCF(Ψ) + EC(F) - EC(Fvalence)
(9)
In fact MCSCF/CASSCF procedures introduce part of via orbital rotations between the core and valence Ecore/valence C subspaces, which strictly should be eliminated. However, this effect is numerically small and may typically be neglected. A Difficult Test Case: F3The trifluoride anion, F3-, is an electronically challenging molecule for ab initio methods.16 The experimental evidence for F3- is incomplete, it has not been isolated so far, and its existence is inferred from vibrational spectroscopy and matrix isolation studies. When a mixture of CsF, RbF, or KF vapor and F2 is deposited with a large excess of inert gas, two bands are seen in the vibrational spectra which exhibit the infrared/ Raman mutual exclusion characteristic of a molecule with a center of symmetry.17,18 For F3- this implies a linear molecule with D∞h symmetry. Many ab initio studies have been done16,19 and also some Kohn-Sham density functional calculations.20 In this study we have investigated the geometry and stability of F3-, using a combination of a CASSCF wave function and a density functional as described above (eq 9). For the density functional part of the calculation we have used the WilsonLevy correlation functional,21
(aF + b|∇F|F-1/3) dV EC(F) ) ∫ (c + d|∇F|(F/2)-4/3 + rs)
TABLE 1: Total Energies, Bond Lengths, and Binding Energies of F2 and F3-
CASSCF/6-31+G* eq 9/6-31+G* B-LYP/6-31+G* CCSD/DZP+a CCSD(T)/DZP+a experiment CASSCF/6-31+G* eq 9/6-31+G* B-LYP/6-31+G* CCSD/DZP+a CCSD(T)/DZP+a a
energy (au) F2 -198.758 790 -199.404 282 -199.518 240 -199.109 58 -199.119 42 F3-298.215 648 -299.195 242 -299.448 658 -298.777 47 -298.808 70
Re (Å)
De (kJ mol-1)
1.503 1.421 1.445 1.428 1.447 1.412
58 106 200 104 121 160
1.812 1.713 1.781 1.723 1.773
61b 73b 201b 68b 115b
Reference 16. b Energy for the reaction F3- f F2 + F-.
rs ) (3/4πF)1/3 a ) -0.748 60 b ) 0.060 01 c ) 3.600 73 d ) 0.900 00 (10) Fuentealba and Savin22 have recently shown that the WilsonLevy functional is an excellent functional for use with HartreeFock densities. All calculations were done using the CASSCF option in our valence bond codes,23,24 which we have interfaced with the GAUSSIAN 92 suite of programs.25 The numerical integration of eq 10 was carried out using the partitioning of the integrand described by Becke26 and the schemes suggested by Handy et al.27 Starting from Hartree-Fock orbitals, the active space for the CASSCF calculations was initially chosen to be the three σ orbitals lying along the bond axis28,29 and the four electrons contained within them. These orbitals are of bonding (σu), nonbonding (σg), and antibonding (σu) character. It was considered essential to obtain at least a qualitatively correct structure at the CASSCF level (i.e., a linear geometry) before attempting to correct for the dynamic correlation. Using this active space and the 6-31G* basis set, a linear geometry was obtained. However, when the 6-31+G* basis was employed, which contains necessary diffuse functions for describing the anion, a bent geometry was obtained. We extended the number of orbitals in the active space by including higher lying σ orbitals from the virtual manifold of Hartree-Fock orbitals. Inclusion of a further three orbitals gave a linear geometry. We checked the “convergence” of this process by adding further orbitals. While the energy decreased negligibly, the geometry (obtained from analytic gradients) did not change when any further orbitals were added. Hence we adopted this four-electron/six-orbital active space in all subsequent calculations. The total density was partitioned according to eq 6 and used to calculate EC(F) from eq 10. Table 1 shows our results obtained at the CASSCF level and from eq 9. We have also performed Kohn-Sham calculations using the B-LYP30,31 functionals for comparison. We also compare our results with the high-level ab initio calculations of ref 16. The equilibrium bond length of F2 predicted at the CASSCF level is, as anticipated, too long. The application of eq 9 improves this quite significantly, the Re and De predicted are very close to CCSD results obtained in a similar quality basis (DZP+).16 The CCSD(T) results, however, are superior with respect to the estimate of De, but poorer in predicting the experimental bond length. Since we have no corresponding experimental numbers for F3-, we can only compare with the best available ab initio results (i.e., CCSD(T)). Again there is closest correspondence between the results of eq 9 and CCSD. It is interesting to note that Kohn-Sham density functional theory using the B-LYP functionals gives
Combined Wave Functions and Density Functionals
J. Phys. Chem., Vol. 100, No. 24, 1996 10133
TABLE 2: Total Energies, Bond Lengths, and Binding Energies of F2 and F3method
energy (au)
CASSCF/6-311+G(2df) eq 9/6-311+G(2df) B-LYP/6-311+G(2df) CCSD/ TZ2Pf+a CCSD(T)/ TZ2Pf+a Experiment CASSCF/6-311+G(2df) eq 9/6-311+G(2df) B-LYP/6-311+G(2df) CCSD/ TZ2Pf+a CCSD(T)/TZ2Pf+a
F2 -198.819 382 -199.465 802 -199.583 412 -199.251 82 -199.28375c F3-298.300 469 -299.278 204 -299.542 444 -298.989 66 -299.032 09
Re (Å)
De (kJ mol-1)
molecule
1.469 1.397 1.432 1.405 1.405d 1.412
66 119 203 119 146 160
1.805 1.705 1.779 1.701 1.701d
53b 59b 197b 57b 103b
Reference 16. Energy for the reaction F3 f F2 + F . The CCSD(T) energy quoted in ref 16 is incorrect. d Energy calculated at the optimized CCSD geometry. a
b
TABLE 3. Total Energies, Bond Lengths, and Binding Energies of Diatomic Molecules Obtained at the CASSCF Level and with Equation 9a
-
- c
bond lengths that are very close to the CCSD(T) values but gives binding energies that are 65-75% too large. Table 2 shows the analogous data obtained with the 6-311+G(2df) basis. The trends are very similar. The larger basis results in a shortening of the bond lengths in F2 and F3- at all levels of theory. While the binding energy of F2 is predicted to increase with the larger basis, the stability of F3- is predicted to decrease. Comparing with the CCSD and CCSD(T) calculations of ref 16, we can see that our results mimic the CCSD results quite closely but are not able to achieve the CCSD(T) results. While this is encouraging, it is useful to ask why the more accurate CCSD(T) results are not reproduced? The main point that presents itself is the inadequacy of the correlation functional of eq 10. We should stress that this functional performs relatiVely well and certainly no worse than more widely used density functionals.32 Nevertheless, the correlation energy of the F atom is significantly underestimated. Equation 10 gives a correlation energy of -0.313 au compared with the exact value of -0.322 au (∆ ≈ 24 kJ mol-1). Note that our number for the F atom is not spherically averaged (unlike the data in ref 32). The situation is considerably worse for F-: eq 10 gives -0.365 au compared with the exact value of -0.397 au (∆ ≈ 84 kJ mol-1). Hence, we may suppose that the correlation in F3- is severely underestimated and the most obvious way to improve our data would be to find a better correlation functional. We are currently investigating some recent correlation functionals suggested in the literature. Diatomic Bond Lengths and Dissociation Energies The main purpose of this paper is to suggest an economical procedure for improving results obtained at the CASSCF level of theory. Hence, we have applied eq 9 to the calculation of the bond lengths and dissociation energies of a series of diatomic molecules using large basis sets.33,34 The choice of active space in these molecules is straightforward and consists of the bonding and antibonding orbital for each two-electron bond in the molecule. The results are presented in Table 3. The mean absolute error in bond length is reduced from 0.023 Å at the CASSCF/CC-pVTZ level to 0.019 Å at the eq 9/CC-pVTZ level. CASSCF bond lengths tend to be overlong, and the inclusion of the dynamic correlation shortens them (in line with expectation). Dissociation energies (De) are much more significantly improved through the use of eq 9. The mean absolute error in De is reduced from 24.9 kcal mol-1 at the CASSCF/CC-pVTZ level to 7.3 kcal mol-1 at the eq 9/CC-pVTZ level.
N2 P2 PN Cl2 HCl HF FCl LiH
energy
Re
Re expt |∆Re|
-109.119 435 -681.571 218 -395.301 788 -919.020 112 -460.122 904 -100.081 455 -558.936 284 -7.986 162
CASSCF 1.104 1.098 1.929 1.893 1.509 1.491 2.047 1.988 1.289 1.275 0.916 0.917 1.664 1.628 1.608 1.596
-109.404 571 -682.612 765 -395.969 748 -920.471 969 -460.851 119 -100.413 142 -559.984 834 -8.050 463
Equation 9 1.097 1.098 1.852 1.893 1.493 1.491 1.938 1.988 1.269 1.275 0.912 0.917 1.603 1.628 1.616 1.596
mean N2 P2 PN Cl2 HCl HF FCl LiH
De
0.006 203.6 0.036 87.4 0.018 118.3 0.059 38.0 0.014 90.0 0.001 113.5 0.036 34.8 0.011 34.4
De expt |∆De| 228.4 117.2 148.6 58.0 106.5 141.2 61.4 58.0
0.023
mean a
0.001 0.041 0.002 0.050 0.006 0.005 0.025 0.020 0.019
24.8 29.8 30.3 20 16.5 27.7 26.6 23.6 24.9
221.6 112.9 142.9 63.6 104.2 124.9 53.9 47.9
228.4 117.2 148.6 58.0 106.5 141.2 61.4 58.0
6.8 4.3 5.7 5.6 2.3 16.3 7.5 10.1 7.3
basis33
The CC-pVTZ was used for all calculations except LiH for which the 6-311G(2df,2pd) basis34 was used. Total energies in au, Re in Å, and De in kcal mol-1.
Further reduction of the errors could be achieved by expanding the active space used in the CASSCF calculations. Clearly, however good the treatment of the Ecore and Ecore/valence is via C C (Ψ) is eqs 9 and 10, unless a good enough estimate of Evalence C obtained, the results will be poor. There is obviously a great (Ψ) deal of flexibility in how one chooses to improve Evalence C and further investigation of this is essential. Conclusion The computational scheme we have proposed (eq 9) is very simple to apply in the context of the CASSCF method. The results are quite encouraging; they mimic results obtained with the CCSD method quite closely but are not able to reproduce results obtained with the CCSD(T) method. We are investigating the use of other density functionals in an attempt to improve our approach. In the Introduction we mentioned that MCSCF/CASSCF techniques can be applied to quite large chemical problems. However, these methods become intractable as the number of electrons in the active space grows past ≈12-14. For example, a 20-electron/20-orbital CASSCF would contain >5.9 × 109 configurations and is currently beyond the ability of any CASSCF code. However, techniques such as CCSD/CCSD(T) are designed to approximate the results of full CI calculations, and they are quite successful at this. Given that a 20electron/20-orbital calculation is a trivial matter for the CCSD method, it would be possible to combine eq 9 with a CCSD reference function. This could easily be achieved by starting with a Hartree-Fock calculation, applying a window to the orbital space in which the CCSD calculation is performed (i.e., freezing low-lying occupied orbitals and high-lying virtual orbitals), and then treating the rest of the correlation problem with a density functional. This is one approach whereby a large chemical system could be treated. The orbital window would be chosen to include those orbitals which undergo rearrangement or show significant depopulation. Hence, the method described here can be extended to quite general reference wave functions, yielding a useful hybrid wave function/density functional approach to the correlation problem.
10134 J. Phys. Chem., Vol. 100, No. 24, 1996 Acknowledgment. We thank the Engineering and Physical Sciences Research Council for the award of a postdoctoral research associateship (Grant GR/K01537) during which this work was carried out. References and Notes (1) Roos, B. O. AdV. Chem. Phys. 1987, 69, 399. (2) Werner, H.-J. AdV. Chem. Phys. 1987, 69, 1. (3) Frisch, M.; Ragazos, I. N.; Robb, M. A.; Schlegel, H. B. Chem. Phys. Lett. 1992, 189, 524. (4) Knowles, P. J.; Werner, H.-J. Chem. Phys. Lett. 1985, 115, 524. (5) Knowles, P. J.; Werner, H.-J. J. Chem. Phys. 1988, 89, 5803. (6) Wolinski, K.; Sellers, H. L.; Pulay, P. Chem. Phys. Lett. 1987, 140, 225. (7) McDouall, J. J. W.; Peasley, K.; Robb, M. A. Chem. Phys. Lett. 1988, 148, 183. (8) Andersson, K.; Malmquist, P.-A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. J. Phys. Chem. 1990, 94, 5483. (9) Laidig, W. D.; Saxe, P.; Bartlett, R. J. J. Chem. Phys. 1987, 86, 887. (10) Hoffmann, M. R.; Simons, J. J. Chem. Phys. 1988, 88, 993. (11) Becke, A. D. J. Chem. Phys. 1992, 96, 2155. (12) Clementi, E.; Chakravorty, S. J. J. Chem. Phys. 1990, 93, 2591. (13) Lie, G. C.; Clementi, E. J. Chem. Phys. 1974, 60, 1275. (14) Lie, G. C.; Clementi, E. J. Chem. Phys. 1974, 60, 1288. (15) Colle, R.; Salvetti, O. J. Chem. Phys. 1983, 79, 1404. (16) Heard, G. L.; Marsden, C. J.; Scuseria, G. E. J. Phys. Chem. 1992, 96, 4359. (17) Ault, B. S.; Andrews, L. J. Am. Chem. Soc. 1976, 98, 1591.
Malcolm and McDouall (18) Ault, B. S.; Andrews, L. Inorg. Chem. 1977, 16, 2024. (19) Wright, T. G.; Lee, E. P. F. Mol. Phys. 1993, 79, 995. (20) Sosa, C.; Lee, C.; Fitzgerald, G.; Eades, R. A. Chem. Phys. Lett. 1993, 211, 265. (21) Wilson, L. C.; Levy, M. Phys. ReV. B 1990, 41, 12930. (22) Fuentealba, P.; Savin, A. Chem. Phys. Lett. 1994, 217, 566. (23) McDouall, J. J. W. Theor. Chim. Acta 1992, 83, 339. (24) Malcolm, N. O. J.; McDouall, J. J. W. J. Comput. Chem. 1994, 15, 1357. (25) Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. W.; Wong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Robb, M. A.; Replogle, E.S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzalez, C.; Martin, R. L.; Fox, D. J.; DeFrees, D. J.; Baker, J.; Stewart, J. J. P.; Pople, J.A. GAUSSIAN 91; Gaussian Inc.: Pittsburgh, PA, 1992. (26) Becke, A. D. J. Chem. Phys. 1988, 88, 2547. (27) Murray, C. W.; Handy, N. C.; Laming, G. J. Mol. Phys. 1993, 78, 997. (28) Pimentel, G. C. J. Chem. Phys. 1951, 19, 446. (29) Rundle, R. E. SurV. Prog. Chem. 1963, 1, 81. (30) Becke, A. D. Phys. ReV. B 1988, 38, 3098. (31) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (32) See, for example, the compilation in: Wilson, L. C. Chem. Phys. 1994, 181, 337. (33) Dunning, T. H., Jr. J. Chem.Phys. 1989, 90, 1007. Woon, W. E.; Dunning, T. H., Jr. J. Chem. Phys. 1993, 96, 6796. (34) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650.
JP960489B
