Article pubs.acs.org/JCTC
Comparative Study of Nonhybrid Density Functional Approximations for the Prediction of 3d Transition Metal Thermochemistry John J. Determan,†,∥ Katelyn Poole,†,‡ Giovanni Scalmani,§ Michael J. Frisch,§ Benjamin G. Janesko,*,† and Angela K. Wilson⊥ †
Department of Chemistry, Texas Christian University, Fort Worth, Texas 76129, United States Department of Chemistry and Center for Advanced Scientific Computing and Modeling (CASCaM), University of North Texas, Denton, Texas 76203-5017, United States § Gaussian, Inc., 340 Quinnipiac Street, Building 40, Wallingford, Connecticut 06492, United States ⊥ Department of Chemistry, Michigan State University, East Lansing, Michigan 48824 United States ‡
S Supporting Information *
ABSTRACT: The utility of several nonhybrid density functional approximations (DFAs) is considered for the prediction of gas phase enthalpies of formation for a large set of 3d transition metal-containing molecules. Nonhybrid DFAs can model thermochemical values for 3d transition metal-containing molecules with accuracy comparable to that of hybrid functionals. The GAM-generalized gradient approximation (GGA); the TPSS, M06-L, and MN15-L meta-GGAs; and the Rung 3.5 PBE+ΠLDA(s) DFAs all give root-mean-square deviations below that of the widely used B3LYP hybrid. Modern nonhybrid DFAs continue to show utility for transition metal thermochemistry.
■
dependence on the electron density gradients ∇ρσ(r). Third rung meta-GGAs incorporate the kinetic energy density of the noninteracting Kohn−Sham reference system τ σ(r) and/or the density Lapacian. The LSDA, GGAs, and meta-GGAs are conventionally referred to as “semilocal” approximations, see for example ref 7. We will use this convention here. The exchange components of semilocal DFAs are implicitly or explicitly8 constructed from exchange hole models in which the σ-spin density matrix and exchange hole about reference point r are modeled using the density or Kohn−Sham orbitals in an infinitesimal volume element about r
INTRODUCTION Kohn−Sham density functional theory (DFT) is widely used in computational chemistry. Relatively simple density functional approximations (DFAs) for the exact exchange-correlation (XC) functional can often model a broad range of properties of atoms, molecules, solids, and surfaces at a mean-field cost.1 Density Functional Approximations. The cost and accuracy of DFT calculations largely depend on the choice of DFA, which provides many-body corrections to the underlying mean-field calculation. Most modern calculations are based on a “Jacob’s Ladder” of DFAs.2 Higher-rung DFAs incorporate additional ingredients and can satisfy more exact constraints.3 The first rung is the local spin-density approximation (LSDA) in which the XC energy density at point r is constructed from the σ = ↑,↓-spin densities ρσ(r). The LSDA is widely used for simulating properties of bulk metals, such as lattice parameters and phonon frequencies,4 but is generally too overbinding to predict realistic chemistry. As an illustration, LSDA calculations on the atomization energies of the 148 molecules of the G2 test set give a mean absolute error >80 kcal/mol.5,6 Second rung generalized gradient approximations (GGAs) incorporate a © XXXX American Chemical Society
Exσ SL = (1/2)
∫ d3rρσ (r) ∫ d3r′hSL Xσ(r, r′)/ r − r′
hXσ SL(r, r′) = −ρσ −1(r) γ SL(ρσ (r), ∇ρσ (r), τσ(r), r − r′)
(1) 2
(2) Received: July 27, 2017 Published: September 6, 2017 A
DOI: 10.1021/acs.jctc.7b00809 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation Localization of hXσSL(r,r′) about the reference point r simulates left−right correlation in covalent bonds.9 However, these DFAs tend to overdelocalize electrons giving inaccurate estimations of many properties.10 Semilocal DFAs tend to overstabilize molecules relative to separated atoms,11,12 overestimate molecular surface adsorption energies, and underestimate band gaps13 and many reaction barriers.14 Hybrid DFAs combine a fraction of eq 1 with a fraction of Hartree−Fock-like (HF) exchange obtained by substituting into eq 1 hXσ(r,r′) = −ρσ−1(r)|γσ(r,r′)|2, where γσ(r,r′) is the one-particle density matrix of the noninteracting Kohn−Sham reference system. This “exact” exchange functional provides the exact high-density limit of the exchange-correlation functional for nondegenerate ground states. “Rung 3.5” DFAs15 instead use an exchange energy density obtained by combining γσ(r,r′) with an LSDA or GGA model density matrix. For example, the PBE model density matix γPBE(ρσ(r),∇ρσ(r),r − r′) introduced in ref 16 yields the ΠPBE exchange energy Exσ ΠPBE = −1/2
formation and limited multireference character. The heats of formation of these 193 molecules form what is known as the ccCA-TM/11 set.41 A follow-up study considered the performance of fourth- and fifth-rung DFAs for this data set.23 Recently, several studies have shown the data from the ccCA-TM/11 set (or subsets that have extensively been based on the ccCA-TM/ 11 set) to be useful for assessment of various density functional methods.42−45 The present work complements these studies by considering DFAs that do not include exact exchange, including modern meta-GGAs and Rung 3.5 DFAs. This work compares the predictions of several modern GGA, meta-GGA, and Rung 3.5 DFAs for the three-dimensional (3d) transition metal heats of formation collected in ref 41. Updating ccCA-TM/11. The ccCA-TM/11 set provides a diverse assessment of 3d transition metal thermochemistry. However, several thermochemical values have been reassessed and/or reevaluated in recent experiments. Both the studies of Stanton and co-workers as well as a study by Dixon and coworkers have suggested that the values of several metal hydrides should be reassessed.44,45 Additionally, a new study by Morse and co-workers has found a higher value for the enthalpy of formation of VC than the value that was used in the original ccCA-TM/11 set (197.7 ± 2 vs 181.6 ± 15 kcal/mol).46 These newer values have been incorporated into the ccCA-TM/11 set. The updated set will be referred to as the ccCA-TM/11:17 set.
∫ d3rρσ (r) ∫ d3r′γσ(r, r′)
γ PBE(ρσ (r), ∇ρσ (r), r − r′)/|r − r′|
(3)
where −(1/2)|γPBE(ρσ(r),∇ρσ(r),r − r′)|2/ρσ(r) provides a model for the exchange hole corresponding to the Perdew− Burke−Ernzerhof (PBE) GGA for exchange. The form of γPBE is presented in ref 17. Hybrid DFAs’ HF exchange admixture tunes the balance between overdelocalization and simulation of nondynamical correlation.18 The optimum fraction of HF exchange tends to differ for different properties. This has led to a broad range of global hybrids with different constant fractions of HF exchange,19 range-separated hybrids including different fractions of HF exchange at different interelectronic separations,20 and “hyper-GGAs” incorporating position-dependent HF exchange admixture.21,22 Rung 3.5 DFAs provide an alternative admixture of nonlocal information. DFT for Transition Metal Complexes. Calculations utilizing standard DFAs are widely used in modeling transition metal complexes.23−33 New DFAs are regularly assessed for transition metal thermochemistry.34−37 Such assessments are complicated by many factors. These include limited gas-phase thermochemical data available from experiments, which can make the development of a reliable molecular properties set challenging. The diversity of transition metal chemistry means that accurate assessments require testing a variety of transition metals as well as a variety of bonding environments. Low-lying excited states, spin−orbit coupling effects, relativistic effects, enhanced electron correlation effects, metal−ligand pi backdonation, and core−core and core−valence interactions may also play significant roles in transition metal chemistry. In contrast to main-group thermochemistry, where hybrid DFAs typically outperform GGAs and meta-GGAs, previous benchmark studies have suggested that GGAs and meta-GGAs can give performance comparable to that of hybrids for transition metal complexes.24,28−30,34,38 Many of these previous studies suffered from sample sets that were relatively small and/or lacking diversity. To assess the accuracy of the correlation consistent composite approach (ccCA) developed in the Wilson research group,39,40 a large, diverse data set composed of thermochemical values for 225 transition metal species was introduced. Of these, 193 molecules were selected based on experimental heats of
■
METHODOLOGY Several GGA, meta-GGA, and Rung 3.5 DFAs are considered. The LSDA is evaluated using the Vosko−Wilk−Nusair correlation functional V.47 Tested GGAs include the nonempirical Perdew−Burke−Ernzerhof (PBE) GGA,48 the modestly empirical BPW91 GGA, 49,50 and the newer SOGGA11 and GAM.34,51 Tested meta-GGAs include the nonempirical Tao−Perdew−Staroverov−Scuseria (TPSS)5 and the empirical M06-L,30 M11-L,52 and MN15-L.38 Three different Rung 3.5 DFAs are considered. Π1PBE combines the PBE GGA for correlation, 50% of the PBE GGA for exchange, and 50% of eq 3. PBE+ΠLDA(s) combines the PBE GGA with admixture of Rung 3.5 exchange in low-gradient regions E XC = E XC(PBE) + Σσ ( −eXσ PBE(r) +
∫ d3r
exp[−ωs2 σ (r)]
∫ d3r′γσ(r, r′)γLDA(ρσ (r), r − r′)/|r − r′| (4)
γLDA(ρσ (r), r − r′) = ρσ (r)exp[−bρσ (r)|r − r′|2 ]
(5)
Here, sσ(r)=|∇ρσ(r)|/(2(6π2)1/3ρ4/3σ(r)) is the reduced density gradient. Unitless empirical parameter ω = 3. PBE+ΠPBE(s) replaces γLDA with γPBE in eq 4. The r − r′ dependence of γLDA and γPBE in eqs 3 and 4 is expanded in an auxiliary basis set of Gaussian-type functions centered at each grid point r. As these are grid point-centered rather than atom-centered auxiliary basis functions, there is no need to use different auxiliary basis sets for different atoms. Because γLDA is spherically symmetric about r − r′, the auxiliary basis set used to expand its r − r′ dependence requires only spherically symmetric (“s-type”) functions. Full details of the tested functionals and the auxiliary basis set implementation are presented in ref 53. The PBE0,54 B3LYP,55 B97-1,56 and MN1538 hybrid DFAs are included for comparison. B
DOI: 10.1021/acs.jctc.7b00809 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation Table 1. Error Statistics (kcal/mol) for the 193 ΔH0f(298 K) of the ccCA-TM/11:17 Data Set type LSDA GGA
meta-GGA
Rung 3.5
Hybrid
DFA
ME
MAE
RMSD
SVWN5 PBE SOGGA11 BPW91 GAM M11-L TPSS TPSSa M06-L MN15-L Π1PBE PBE+ΠPBE(s) PBE+ΠLDA(s) PBE+ΠLDA(s)a PBE0 B3LYPa B3LYPb B97-1 MN15
88.5 25.0 19.4 13.4 −1.9 −2.1 10.9 10.7 6.9 1.3 14.8 −13.1 −10.5 −9.2 −8.9 −15.4 −13.9 −7.2 −3.5
88.5 26.1 22.2 17.4 9.9 14.1 13.0 12.8 11.3 8.9 20.0 15.4 14.3 13.5 13.9 17.1 15.4 9.0 10.2
133.4 39.1 30.5 24.6 13.9 21.0 17.2 17.3 15.8 12.8 27.4 25.8 24.2 22.8 23.0 25.0 22.3 13.9 15.8
MAX 546.9 178.9 120.3 123.0 32.4 47.1 85.8 87.6 52.1 49.1 91.2 14.5 22.2 24.1 33.3 15.5 16.0 11.0 25.3
Sc(C5H5)3 V4O10 V4O10 V4O10 VF5 CrBr4 V4O10 V4O10 Sc2Cl6 Fe2Cl6 V4O10 ScC4 ScC4 TiC4 Sc(C5H5)3 VH VH CoH FeF2
MIN 0.8 −23.3 −53.1 −41.9 −61.7 −117.7 −33.4 −31.0 −72.3 −62.3 −63.2 −184.8 −190.8 −181.9 −158.9 −98.6 −99.4 −60.9 −120.8
Zn2 Cr2Cl4 Ti2 Cr2Cl4 Ni(PF3)4 Ni(PF3)4 Ni(PF3)4 Ni(PF3)4 Ni(PF3)4 Cr3O9 Cr2 Ni(PF3)4 Ni(PF3)4 Ni(PF3)4 Cr3O9 Cr3O9 Cr3O9 Ni(PF3)4 Cr3O9
a
Evaluated with cc-pVTZ optimized geometries, zero-point corrections, and thermal corrections. bEvaluated with cc-pVQZ basis, optimized geometries, zero-point corrections, and thermal corrections.
13.9 kcal/mol, respectively). The good performance of B97-1 is consistent with previous work.23 The results for the LSDA and GGAs are consistent with previous work. The LSDA severely overbinds, producing an RMSD value of 133.4 kcal/mol. Density gradient corrections improve the results. All tested GGAs have MAEs and RMSDs below those of LSDA. The extent of improvement depends on the functional form; for example, the BPW91 GGA gives an RMSD 15 kcal/mol smaller than that of PBE, and the recent GAM GGA gives the best performance of tested GGAs. The results for meta-GGAs are also consistent with previous work. The MN15-L meta-GGA gives the lowest overall RMSD (12.8 kcal/mol) compared with the available benchmarks. The TPSS, M06-L, and MN15-L meta-GGAs all result in RMSD values (17.2, 15.8, and 12.8 kcal/mol, respectively) lower than the widely used B3LYP global hybrid (RMSD of 25.0 kcal/ mol). Geometry optimization with TPSS gives a small change in the error statistics. These meta-GGAs also give relatively small maximum and minimum errors, indicating improved reliability. The Rung 3.5 DFAs give creditable performance for these benchmarks, though their performance is not on par with the best meta-GGAs. The PBE+ΠLDA(s) Rung 3.5 DFA gives RMSD that is comparable to that of the B3LYP global hybrid. Interestingly, full geometry optimization of the complex geometries further reduces the RMSD. Other Rung 3.5 DFAs are more problematic, with the simple ΠLDA overbinding worse than the LSDA itself. This is consistent with its performance for main-group thermochemistry.15 Importantly, global hybrids and Rung 3.5 DFAs both give mean deviations relative to experiment that are negative (underbinding), whereas the tested GGAs and meta-GGAs tend to overbind. This observation that Rung 3.5 DFAs are more like hybrid DFAs and less like semilocal DFAs is consistent with the fact that eq 3 incorporates nonlocal information similar to that of hybrid DFAs. The five species Ni(PF3)4, Cr3O9,V4O10, Fe2Cl6, and Sc2Cl6 are outliers for many of the tested DFAs, accounting for most
Calculations use the development version of the Gaussian suite of programs.57 All calculations are performed selfconsistently. Calculations expand the Kohn−Sham orbitals in the atom-centered cc-pVTZ basis set. Some test calculations use the cc-pVQZ basis set. Calculations use B3LYP/cc-pVQZ geometries, zero-point, and thermal corrections.58−61 Some test calculations obtain geometries, zero-point energies, and thermal corrections with the tested functional. Geometries are optimized with spatial symmetry constraints. All computed enthalpies of formation include spin−orbit corrections for selected diatomic molecules (selected diatomics can be found in the Supporting Information) and include spin−orbit corrections for isolated atoms. (These results differ slightly from those reported in ref 23. One may reproduce the results of ref 23 using our computed values by including spin−orbit corrections for molecules and not including spin−orbit corrections for isolated atoms.) Stability analysis is performed on all atomic wave functions and on the LSDA, BPW91, M06L, TPSS, PBE+ΠLDA(s), and PBE0 molecular calculations.62 Error statistics include mean signed errors ME (evaluated as experiment − theory), mean absolute errors (MAEs), rootmean-square deviations (RMSDs), and most positive (MAX) and most negative (MIN) deviations from experiment. Experimental enthalpies of formation, computed geometries, computed total energies, computed zero-point and thermal enthalpy corrections, and spin−orbit corrections of all atoms and all molecules evaluated at all levels of theory are included as Supporting Information.
■
RESULTS
Overall Performance. Table 1 summarizes the statistical errors for all 193 heats of formation in the ccCA-TM/11:17 set. As in previous studies of transition metal chemistry, the bestperforming DFAs do not necessarily involve HF exchange. The GAM GGA, MN15-L meta-GGA, and B97-1 global hybrid give the best performance among the tested DFAs (producing RMSD values relative to experimental values of 13.9, 12.8, and C
DOI: 10.1021/acs.jctc.7b00809 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
methods can result in large errors and unpredictable behavior when computing energies and spectroscopic properties for these molecules. Jiang and co-workers thus identified 23 molecules from the ccCA-TM/11 molecule set as having possible significant multireference character (as defined as having T1 and D1 values larger than the criteria listed above). The remaining 170 molecules can be referred to as the subset ccCA-TM/11:SR-170. Zhang and co-workers selected 70 molecules identified previously by authors of the ccCA-TM/ 11 set that had low experimental uncertainty (
