J. Phys. Chem. 1995,99, 17085-17087
17085
Effective Core Potentials for DFT Calculations Thomas V. RUSSO,*Richard L. Martin, and P. Jeffrey Hay Theoretical Division, MS B268, Los Alamos National Laboratory, Los Alamos, New Mexico 87545 Received: August 9, 1995; In Final Form: October 12, 1995@
Effective core potentials have been generated for Ti and Ni using atomic density functional theory (DFT) wave functions within the local density approximation. W e find that these effective potentials give good agreement between all-electron and valence-electron calculations in TiF4 and Ni(C0)4 for both the localdensity (S-VWN) and the gradient-corrected (B-LYP) approximation. Furthermore, we demonstrate for a small but representative set of transition metal complexes that the effective core potentials generated previously from Hartree-Fock atomic calculations may be used in DFT-based methods as well.
Introduction The dynamic electron correlation effects which make accurate calculations based on the Hartree-Fock approximation so difficult for transition metal complexes appear to be treated reasonably accurately with the new generation of “gradientcorrected” density functional theories (DFT). Our recent study’ of the excitation and ionization energies of the atoms of the fisst transition series showed that the Becke exchange functional2 coupled with the Lee-Yang-Pan correlation fun~tional,~ the (B-LYP) functional, gave exceptional results for the ionization energies (mean errors of -0.1 eV) and acceptable results for excitation energies (mean errors of -0.5 eV). This can be compared with mean errors in the Hartree-Fock (HF) approximation of -1.7 and -0.8 eV, respectively. This level of agreement with experiment is consistent with the experience of other^^-^ using different gradient-corrected exchange and correlation functionals. The number of applications of these methods to inorganic complexes is increasing rapidly, and the thermochemistries reported for both “gradient-corrected” DFT methods and hybrid approaches, in which a component of the exact HF exchange interaction is included in the Hamilt~nian,~are very enco~raging.~-I~ As was observed in calculations on main group molecule^,'^ the gradient corrections serve to moderate the tendency of the local density approximation (LDA) to overbind. For example, in a study of ScF3, TiF4, VF5, and crF6, we foundI3 that the B-LYP method overestimated the average bond energy in the series by -8 kcdmol. This can be compared with an error of -20-30 kcal/mol for the LDA and errors of -30-70 kcal/mol for Hartree-Fock. Remarkably, the hybrid B3-LYP a p p r ~ a c h , ~which , ’ ~ uses a component of the exact exchange, was in error by only -1-2 kcal/mol throughout this series. Ricca and BauschlicherI2 obtained similarly accurate results with the B3-LYP approach in a study of the binding energies of a series of metal carbonyls. Although there is certainly room for improvement in the functionals, the level of accuracy achieved, when considered with the computational scaling advantages of these methods, makes them quite attractive. In order for these methods to reach their full potential in inorganic chemistry, however, the effects of relativity must be addressed. Approaches to relativistic DFT which treat all the electrons in the molecule appear quite promising.16-20 Alternatively, the use of effective core potentials (ECPs), which replace the chemically inert core electrons with an effective @
Abstract published in Advance ACS Abstracts, November 15, 1995.
0022-365419512099-17085$09.00/0
potential, can be used to incorporate the dominant effects of relativity and have the advantage that the size of the basis set needed to treat the molecule can be significantly reduced. ECPs generated from atomic LDA cores have been used quite effectively in LDA calculations in both the solid-state2’-28and quantum chemistry c o m u n i t i e ~ .Relativistic ~ ~ ~ ~ ~ ECPs generated from the HF core are also used extensively in quantum c h e m i ~ t r y ~ for ’ - ~HF and post-HF calculations. In this letter we address some issues regarding the use of ECPs in DFT calculations. In particular, one might expect that an ECP appropriate for use in a specific DFT calculation should be based on the analogous atomic DFT calculation. Given the large number of combinations of exchange and correlation functionals available in the literature, this could prove to be a formidable task. On the other hand, in our study of transitionmetal atoms, we found that a basis set in which the core orbitals were contracted with the atomic HF coefficients gave results for DFT (B-LYP) calculations in good agreement with the fully uncontracted basis. This suggests that a frozen-core approximation, with the core determined via the HF method, is an acceptable approximation in DFT work. This is implicitly assumed in studies of main group chemistry as well, where basis sets which heavily contract the core, such as the 6-31G basis, are used with good success. One might be tempted, therefore, to conclude that ECPs based on atomic HF calculations might also be utilized in DFT, and in fact, we have noticed a number of recent publications in which this approach has been used without justification. These studies are suspect, however, for a number of reasons. For example, because the exchange and correlation functionals are nonlinear in the density, the core and valence densities are not strictly separable. A frozen-core approximation is therefore somewhat different conceptually from an ECP approach which replaces the sum of the core and valence density with a pseudo “valence only” density. In addition, the relativistic influence of the core on the valence electrons can be captured in a relativistic ECP, whereas it is less clear how to accomplish this with a frozen-core approach. For these and other reasons, then, the use of ECPs generated from HF calculations must be carefully tested for their applicability in DFT work. In this work we investigate two questions. The first is the extent to which ECPs derived from LDA atomic cores can be utilized in gradient-corrected D R calculations. In order to study this question, new ECPs have been generated for Ti and Ni using the Slater exchange and Vosko-Wilks-Nusair correlation functionals (S-VWNI4) and the approaches described the sole difference being that previously for the HF 0 1995 American Chemical Society
17086 J. Phys. Chem., Vol. 99, No. 47, 1995 TABLE 1: Equilibrium Bond Lengths
Letters
(A) Calculated Using Hartree-Fock, HF
molecule TiF, Ni(C0)J
r w rco
AE 1.755 1.833 1.134
HFECP 1.756 1.838 1.133
AE 1.744
1.743 1.170
S-VWN, and B-LYP S-VWN S-VWNECP AE 1.747 1.716 1.753 1.801 1.170 1.180
B-LYP S-VWNECP 1.777 1.815 1.180
All calculations used the (6s5p2d) contraction of the all-electron basis set of Wachters.’ The column labeled AE is an all-electron calculation; the columns marked HFECP and S-VWNECP used effective core potentials generated from either atomic HF’3 or S-VWN calculations described herein.
TABLE 2: Binding Energies (kcaVmo1) Calculated Using Hartree-Fock, S-VWN, and B-LYP“ HF S-VWN molecule AE HFECP AE S-VWNECP AE TiF4 316.8 323.3 666.9 664.1 571.5 -81.9 285.5 279.0 186.9 Ni(C0)J -85.3
B-LYP S-VWNECP 564.1 173.7
‘All calculations used the (6s5p2d) contraction of the all-electron basis set of Wachters.’ The column labeled AE is an all-electron calculation; the columns marked HFECP and S-VWNECP used effective core potentials generated from either atomic HF33or the S-VWN calculations described herein. For Ni(C0)4, the binding energy is relative to Ni 4CO.
+
numerical S-VWN atomic wave functions were used as the starting point. The second question we address is the appropriateness of using HF-based E C P S in ~ ~LDA and gradientcorrected DFT calculations. We conclude that the differences between all-electron DFT calculations and those using ECPs derived from either LDA or HF atomic calculations are comparable and of the same order of magnitude as those found previously in HF and post-HF calculations.
Computational Details The DFT calculations were carried out using the Truchas suite of program^.^' The Kohn-Sham procedure used is strictly analogous to traditional Hartree-Fock calculations, with the replacement of the analytical exchange matrix with an exchangecorrelation matrix calculated by numerical integration. Analytic gradients of the energy were computed including the derivatives of the quadrature weight^;'^ the ECP calculations utilized the gradient capabilities described recently.42 On the metals we use a simple grid of 100 concentric spheres, with radii chosen according to an Euler-MacLaurin scheme,43 on which are positioned quadrature points chosen according to the scheme of L e b e d e ~ ,using ~ ~ . ~the~ 194-point (23rd-order) formulas. This 100-point radial grid for the transition metal has been shown to provide atomic energies accurate to low4au in a previous study.’ It is much larger than the SG1 grids and could be made significantly more economical with little loss in accuracy by applying a “grid pruning” technique.46-48 The “Standard Grid-1’’ of Gill et al.46is used for main group atoms; this grid prunes the number of grid points by varying the order of the Lebedev grids used based on the radius of the spherical shell. The atomic calculations were performed using a spinrestricted, open-shell1 procedure as described in a previous work.’ Our aim in this work is not to compare the results of ECP calculations with experiment, but rather to compare the predictions of ECP calculations with analogous all-electron (AE) calculations. In order not to cloud the issue, a common basis set was used for both the AE and ECP calculations. For the metal atoms, the (14s9p5d) Cartesian-Gaussian basis set of Wachters” was contracted’ to [6s5p2d] using Wachters’ Hartree-Fock coefficients and the general contraction concepts of RaffenettLso The more diffuse primitives of each space are left free. This contracted set gave B-LYP total energies for ScCr within 10-3 au of the completely uncontracted results and
excitation and ionization energies to within 0.01 eV.’ For the main group atoms, the 6-31G basis set was used.
Results and Discussion Use of Local-Density-Derived ECPs for Gradient-Corrected Calculations. The first issue we wish to address is whether an ECP derived from an atomic S-VWN calculation may be used with a gradient-corrected approximation, such as B-LYP. We shall denote such an ECP as S-VWNECP, in order to distinguish it from a potential derived from an atomic HF calculation, which we denote as HFECP. We have generated new S-VWNECPs which replace the [Ne] core in Ti and Ni and used them in calculations on TiF4 and Ni(C0)d. These two molecules were chosen as prototypical of ionic and covalent bonding, respectively. The first major column of Table 1 compares all-electron HF bond distances with the results using a HFECP, the second compares all-electron S-VWN calculations with S-VWNECP results, and the third reports all-electron B-LYP results versus B-LYP calculations using the S-VWN/ ECP. Table 2 compares the binding energies, the total molecular energy at the equilibrium geometry minus the sum of the atomic energies; for Ni(C0)I the binding energy is calculated relative to atomic nickel and optimized free CO molecules. The first set of entries, HF calculations using the previously derived HFECPs, provide a base line for what has constituted an acceptable error in the ast. Note that the distances are reproduced to within -0.01 and the binding energies to within -2 kcaVmol per bond. The differences observed between allelectron and S-VWNECP calculations are similar in magnitude. The largest difference in bond distances, 0.01 A, occurs for the Ni-C bond in Ni(C0)d. The differences in binding energies are also less than 2 kcaVmol per bond. As has been noted many times before, the S-VWN approximation binds Ni(C0)4 relative to Ni 4C0, whereas the HF approximation predicts that this molecule is unbound. Note also in Tables 1 and 2 that the results of all-electron B-LYP calculations and B-LYP calculations using the S-VWNECP are also in good agreement with one another. The errors observed when using the S-VWNECP for the gradient-corrected calculation are somewhat larger than those observed in the S-VWN calculations, but still acceptable and of the order of -0.01 A and -2 kcal/mol per bond. On the basis of this very limited study, it appears that one can expect accuracy in S-VWN-based ECPs which is comparable to that observed in HF calculations with HFECPs and that the S-VWNECPs can be used in gradient-corrected calculations as well.
1
+
Letters
J. Phys. Chem., Vol. 99, No. 47, 1995 17087
TABLE 3: Equilibrium Bond Lengths (A) Calculated Using Hartree-Fock, S-VWN,and B-LYP HF
s-VWN
B-LYP
molecule
AE
HFECP
AE
HFECP
ZE
HFECP
SCF~ TiF4 VF5
1.857 1.755 1.737 1.701 1.580 1.833 1.134 1.823
1.859 1.756 1.738 1.701 1.579 1.838 1.133 1.825 0.001
1.819 1.744 1.747 1.713 1.641 1.743 1.170 1.726
1.824 1.746 1.750 1.715 1.643 1.750 1.170 1.727 0.001
1.850 1.776 1.783 1.749 1.677 1.801 1.180 1.765
1.854 1.776 1.785 1.750 1.676 1.811 1.180 1.767 0.001
Ti0 Ni(C0)4
rax rq ~N,C
rco CuF
A
a All calculations used the (6s5p2d) contraction of the all-electron basis set of Wachters.] Columns labeled AE are all-electron calculations; columns marked H F E C P used effective core potentials generated previously from atomic HF calculation^.^^ The mean absolute deviation between AE and H F E C P calculations is given in the last row.
TABLE 4: Binding Energies (kcaymol) Calculated Using Hartree-Fock, S-VWN,and B-LYP HF
for all the individual components of the hybrid functionals suggests that they may be used in that case as well.
S-VWN
B-LYP
molecule
AE
HFECP
AE
HFECP
AE
HFECP
SCF~ TiF4 VF5 Ti0 Ni(C0)d CuF A
300.11 316.81 190.44 -21.69 -85.29 35.81
301.79 323.34 194.43 -19.51 -81.86 36.33 0.5
505.81 666.85 728.09 126.37 285.49 108.71
507.21 669.56 727.03 125.40 278.67 108.25 0.4
439.54 571.45 605.59 89.19 186.87 89.94
440.17 570.31 599.00 87.94 173.28 89.45 0.6
All calculations used the (6s5p2d) contraction of the all-electron basis set of Wachters.' Columns labeled AE are all-electron calculations; columns marked HFECP used effective core potentials generated previously from atomic HF calculation^.^^ For Ni(C0)4, the binding energy is relative to Ni 4CO. The mean absolute deviation between all-electron and ECP calculations is given in the last row (per bond).
+
Use of Hartree-Fock-Derived ECPs in Density Functional Calculations. The second issue we address is the appropriateness of HFECPs for DFT calculations. Tables 3 and 4 report the results of such an investigation for a number of inorganic complexes. Once again the first major column provides a base line in which all-electron HF calculations are compared with HFECP calculations. Examination of these tables shows that the use of HFECPs in the D I T calculations leads to errors which are of the same order of magnitude as the error in the HF calculations. The average absolute error in the bond distances for the molecules in Table 3 is 0.001 %, for HF and 0.001 %, for S - V W N and B-LYP as well. For the binding energies, the mean absolute error per bond is 0.5, 0.4, and 0.6 kcal/mol per bond for HF, S-VWN, and B-LYP, respectively. Conclusions We have shown that, for a small but representative set of test cases, the ECPs derived previously from atomic HF calculations may be used with little loss of accuracy in DFT calculations as well. It should be stressed that the comparisons here were performed using a common basis set appropriate for the all-electron case. In particular, we have not shown that the basis sets derived for use with HFECPs may also be used in DFT studies. In fact, we have observed highly erratic results if the I-IFECP basis set is employed. Preliminary results suggest that appropriate DFT basis sets may be generated from the primitive ECP sets if they are contracted according to the coefficients of an atomic S - V W N calculation in the primitive basis. Finally, although we have not examined the hybrid functionals, Le., those which combine a component of the HF exchange with DFT functionals, the fact that the HFECPs perform well
Acknowledgment. This work was sponsored by the US. Department of Energy through the Laboratory Directed Research and Development program at Los Alamos National Laboratory. References and Notes (1) Russo, T. V.; Martin, R. L.; Hay, P. J. J . Chem. Phys. 1994, 101, 7729. (2) Becke, A. D. Phys. Rev. A 1988, 38, 3098. (3) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. E 1988,37, 785. (4) Kutzler, F.; Painter, G. Phys. Rev. E 1991, 43, 6865. (5) Perdew, J. P.; Chevary, J.; Vosko, S.; Jackson, K.; Pederson, M.; Singhand, D.; Fiolhas, C. Phys. Rev. E 1992, 46, 6671. (6) Ziegler, T.; Li, J. Can. J. Chem. 1994, 72, 783. (7) Becke, A. D. J . Chem. Phys. 1993, 98, 5648. (8) Ziegler, T. Chem. Rev. 1991, 91, 651. (9) Delley, B.; Wrinn, M.; Luethi, H. P. J . Chem. Phys. 1994, 100, 5785. (10) Deeth, R. J. J . Phys. Chem. 1993, 97, 11625. (11) Li, J.; Schreckenbach, G.;Ziegler, T. J. Phys. Chem. 1994, 98, 4838. (12) Ricca, A.; Bauschlicher, C. W. J . Phys. Chem. 1994, 98, 12899. (13) Russo, T. V.; Martin, R. L.; Hay, P. J. J . Chem. Phys. 1995, 102, 8023. (14) Johnson, B. G.; Gill, P. M. W.; Pople, J. A. J . Chem. Phys. 1993, 98, 5612. (15) Stevens, P.; Devlin, F.; Chabarowski, C.; Frisch, M. J. J . Phys. Chem. 1994, 98, 11623. (16) Bastung, T.; Heinemann, D.; Sepp, W.; Kolb, D.; Fricke, B. Chem. Phys. Lett. 1993, 211, 119. (17) Ziegler, T.; Tschinke, V.; Baerends, E.; Snijders, J.; Ravenek, V. J . Chem. Phys. 1989, 93, 3050. (18) Haberlen, 0.;Roesch, N. Chem. Phys. Lett. 1992, 199, 491. (19) Haberlen, 0.; Chung, S.-C.; Roesch, N. Znt. J . Quantum Chem. Symp. 1994, 28, 595. (20) Chung, S.-C.; Krueger, S.; Pacchioni, G.; Roesch, N. J . Chem. Phys. 1995, 102, 3695. (21) Hamman, D.; Schlueter, M.; Chiang, C. Phys. Rev. Lett. 1979,43, 1494. (22) Bachelet, G.; Hamman, D.; Schlueter, M. Phys. Rev. E 1982, 26, 4199. (23) Kleinman, L.; Bylander, D. Phys. Rev. Lett. 1982, 48, 1425. (24) Hamman, D. Phys. Rev. E 1989, 40, 2980. (25) Troullier, N.; Martins, J. Phys. Rev. E 1991, 43, 1993. (26) Pickett, W. Comput. Phys. Rep. 1989, 9, 115. (27) Lee, S.; Bylander, D.; Kleinman, L. Phys. Rev. E 1988,37, 10035. (28) Martins, J. L.; Buttet, J.; Car, R. Phys. Rev. E 1985, 31, 1804. (29) Andzelm, J.; Radzio, E.; Salahub, D. J . Chem. Phys. 1985,83,4573. (30) Chen, H.; Krasowski, M.; Fitzergerald, G. Unpublished. (31) Hay, P. J.; Wadt, W. R. J . Chem. Phys. 1985, 82, 270. (32) Wadt, W. R.; Hay, P. J. J . Chem. Phys. 1985, 82, 284. (33) Hay, P. J.; Wadt, W. R. J . Chem. Phys. 1985, 82, 299. (34) Pacios, L.; Christiansen, P. J. Chem. Phys. 1985, 82, 2664. (35) Hurley, M.; Pacios, L.; Christiansen, P.; Ross, R.; Ermler, W. J . Chem. Phys. 1986, 84, 6840. (36) LaJohn, L.; Christiansen, P. A.; Ross, R.; Asharoo, T.; Ermler, W. J . Chem. Phys. 1986, 87, 2812. (37) Ross, R.; Powers, J.; Asharoo, T. J . Chem. Phys. 1990, 93, 6654. (38) Krauss, M.; Stevens, W.; Basch, H.; Jasien, P. Can. J . Chem. 1992, 70, 612. (39) Stevens, W.; Basch, H.; Krauss, M. J . Chem. Phys. 1984,81,6026. (40) Cundari, T.; Stevens, W. J . Chem. Phys. 1993, 98, 5555. (41) Martin, R. L.; Russo, T. V.; Lengsfield, B. H.; Saxe, P. W.; Tawa, G. J.; Page, M.; Hay, P. J.; Rappt, A. K. Truchas, 1994. (42) Russo, T. V.; Martin, R. L.; Hay, P. J.; RappB, A. K. J . Chem. Phys. 1995, 102, 9315. (43) Murray, C . W.; Handy, N. C.; Laming, G. J. Mol. Phys. 1993, 78, 997. (44) Lebedev, V. I . Zh. Vyschisl. Mat. Mar. Fiz. 1975, 15, 48. (45) Lebedev, V. I. Zh. Vyschisl. Mat. Mat. Fiz. 1976, 16, 293. (46) Gill, P. M. W.; Johnson, B. G.; Pople, J. A. Chem. Phys. Lett. 1993, 209, 506. (47) Perez-Jorda, J. M.; Becke, A. D.; San-Fabian, E. J . Chem. Phys. 1994, 101, 6520. (48) Baker, 3.; Andzelm, J.; Scheiner, A,; Delley, B. J . Chem. Phys. 1994, 101, 8894. (49) Wachters, A. J. H. J. Chem. Phys. 1970, 52, 1003. (50) Raffenetti, R. C. J . Chem. Phys. 1973, 58, 4452. JP952296E
