Evaluation of the Performance of the B3LYP, PBE0 ... - ACS Publications

Jan 25, 2016 - provides the best results when used with LOC corrections. (although other functionals ..... keyword in Jaguar, based on ligand-field th...
0 downloads 0 Views 308KB Size
Download PDF
Article pubs.acs.org/JCTC

Evaluation of the Performance of the B3LYP, PBE0, and M06 DFT Functionals, and DBLOC-Corrected Versions, in the Calculation of Redox Potentials and Spin Splittings for Transition Metal Containing Systems Dilek Coskun,† Steven V. Jerome,‡ and Richard A. Friesner*,† †

Department of Chemistry, Columbia University, New York, New York 10027, United States Schrödinger, Inc., New York, New York 10036, United States

‡

S Supporting Information *

ABSTRACT: We have evaluated the performance of the M06 and PBE0 functionals in their ability to calculate spin splittings and redox potentials for octahedral complexes containing a first transition metal series atom. The mean unsigned errors (MUEs) for these two functionals are similar to those obtained for B3LYP using the same data sets. We then apply our localized orbital correction approach for transition metals, DBLOC, in an effort to improve the results obtained with both functionals. The PBE0DBLOC results are remarkably close in both MUE and parameter values to those obtained for the B3LYP-DBLOC method. The M06-DBLOC results are less accurate, but the parameter values and trends are still qualitatively very similar. These results demonstrate that DBLOC corrected methods are substantially more accurate for these systems than any of the uncorrected functionals we have tested and that the deviations between hybrid DFT methods and experiment for transition metal containing systems exhibit striking physically based regularities which are very similar for the three functionals that we have examined, despite significant differences in the details of each model.

I. INTRODUCTION The calculation of redox potentials and spin state splittings for transition metal containing species is an essential task in quantum chemical modeling of materials and metal containing biological systems. These quantities play a critical role in a wide range of transition metal chemistry and physics, including catalysis, electron transfer, and conductivity. However, it is quite challenging to compute accurate results with currently available quantum chemical approaches. High level methods such as CCSD(T) have been applied to very small transition metal containing species5,10 (although there is considerable debate about how accurate such methods are, even for these small systems) but to date have not been extensively employed in studying larger systems that are more relevant to realistic problems in biology and materials science. DFT calculations are tractable for a much larger number of atoms (300−500), but existing DFT functionals do not yet produce robust results for © XXXX American Chemical Society

transition metal containing systems. Hence, new approaches are needed for both interpretation of experimental data and design of improved molecules and materials. In a series of previous papers,11−13,15 we have introduced a new approach to improving the energetics of DFT calculations, based on the use of localized orbital corrections (LOCs) for specific chemical bonds, lone pairs, atomic hybridization states, and a number of other features of molecular structure. For organic systems, we have found that the B3LYP functional provides the best results when used with LOC corrections (although other functionals benefit considerably as well). The B3LYP-LOC method achieves a mean unsigned error (MUE) of 0.03 eV in predicting atomization energies of the G3 set of benchmark molecules; it provides substantial improvements in Received: August 14, 2015

A

DOI: 10.1021/acs.jctc.5b00782 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

either B3LYP or M06 (although they share some parametrization data in common). Neither functional has been evaluated for transition metal performance using data sets of the type we employ here; these results are of interest in and of themselves. The amenability of the two functionals to DBLOC corrections yields some surprising results, as will be shown and discussed below. The paper is organized as follows. In section II, we briefly review the DBLOC models for calculating spin state splittings and redox potentials and give descriptions of the benchmark data sets used in the previous (and current) studies. Section III outlines the computational methods employed in the calculations. Section IV presents results obtained for the B3LYP, M06, and PBE0 functionals for both data sets, as well as the DBLOC-corrected results for all three functionals. In section V we compare the results obtained for all three functionals, with and without DBLOC corrections, and discuss the implications of these findings with regard to the systematic nature of the errors in hybrid DFT as applied to transition metal containing systems. In section VI, the Conclusion, we briefly summarize our results.

other quantities, including activation barriers, ionization potentials, and electron affinities, as well.12,15 Over the past several years, we have extended the B3LYPLOC approach to transition metal containing systems; we refer to this model as the DBLOC model.1,7,20,21 As a starting point, we have assembled large data sets (consisting of simple but diverse organometallic complexes of first row transition metals) to evaluate the performance of DFT methods for the computation of both spin state splittings and redox potentials. Uncorrected B3LYP has a large MUE versus experiment for both data sets, on the order of 0.4 eV. These errors are sufficiently large so as to potentially render DFT computations of both quantities highly problematic for chemical applications (although differences in redox potentials or spin state splittings for closely related molecules may be computed with much better accuracy, due to cancellation of error). In this endeavor, we have found that descriptors based on ligand field theory, in particular the position of ligands in the spectrochemical series, provide a very useful basis for the construction of an LOC type model for transition metals. The error in the position of the eg orbital of the complex compared to the t2g orbital is proportional to the strength of the metal− ligand interaction (which is correlated to the position of the ligand in the spectrochemical series), and this in turn has a significant effect on spin state splittings and redox energies. Other systematic errors, such as the relative energies of singly and doubly occupied d orbitals, can be understood from studying DFT calculations of ionization potentials and spin state transitions for transition metal atoms, which exhibit remarkable regularities in deviation from experiment across the entire set of the first transition metal series, as is shown in detail in ref 4. These errors in coordinatively saturated complexes (we focus our attention here exclusively on octahedral complexes) are not identical to those in atoms, but the numerical values obtained from fitting are quite similar to the atomic results and provide similarly efficacious improvement in bringing the energetics of DFT into agreement with experiment. The B3LYP-DBLOC model achieves a MUE of ∼0.1 eV for both spin state splittings and redox potentials, which is arguably within range of the experimental errors in measurement for these quantities. In our development of the B3LYP-LOC methodology for organic systems, we investigated a wide range of DFT functionals, making the performance of the LOC approach clear in the context of the available options for DFT calculations and determining the efficacy of LOC corrections for each functional (even the LDA is vastly improved, with MUE for the G2 set reduced from 3.9 eV to ∼0.4 eV).11 The objective of the present paper is to carry out a similar set of comparisons for the DBLOC approach. Because the DFT calculations for our transition metal containing data sets are much more computationally intensive and difficult to converge (for example, in many cases, one has to try a large number of initial guesses to obtain the correct spin state for a given complex), we restrict our investigations here to two alternative functionals; the M06 functional of Truhlar and co-workers9 and the PBE0 functional of by Adamo and Barone.8 M06 is one of a new generation of functionals which has gained considerable popularity over the past decade; it has been specifically optimized for transition metal performance (although the training set consisted primarily of very small, coordinatively unsaturated species). The PBE0 functional was also widely used and was developed employing a very different approach than

II. DBLOC MODELS AND DATA SETS FOR SPIN SPLITTING AND REDOX POTENTIALS In this section, we briefly review the underlying rationale, and practical implementation, of the DBLOC methodology as it pertains to redox potentials and spin splittings. The theoretical basis of the general LOC approach is discussed in refs 11 and 12; the basic ideas are outlined below. DBLOC has been presented previously in considerable detail in refs 1 and 7. The localized orbital correction (LOC) model is based on the idea that the dominant nondynamical correlation energy in a molecule arises from two electrons in an orbital (either a bond or a lone pair) avoiding each other. For example, in the H2 molecule, the Hartree−Fock solution completely fails to capture left−right correlation, which plays a critical role in the substantial bond energy of H2. Furthermore, this type of correlation is handled with quantitative accuracy by the local density approximation (LDA) which is based on the uniform electron gas, a model lacking in concentrated centers of positive charge. In ref 11, we started from the hypothesis (which had been previously suggested by other workers, for example in refs 22−24) that modern hybrid DFT functionals like B3LYP approximately model the nondynamical pair correlation energy via their self-interaction “error” term. The addition of an exact exchange component (which has no self-interaction error at all) enables the magnitude of the self-interaction term arising from the Becke exchange component17 to be optimized to fit experimental atomization data; this is how the 20% exact exchange fraction in B3LYP is obtained. A detailed exposition of this argument can be found in ref 11. The optimization of the exchange fraction in B3LYP represents an attempt to model the pair nondynamical correlation energy averaged over all bonds and lone pairs. However, it is to be expected that such a representation will not provide perfect accuracy for every electron pair; some will be closer to the average than others. In this picture, pairs which deviate from the average substantially, such as those in an ionic bond (e.g., in the NaCl bond), will therefore often exhibit substantial, and systematic, deviations from experimental results. A simple empirical correction scheme, in which it assumed that the error for a given type of electron pair is B

DOI: 10.1021/acs.jctc.5b00782 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

reoptimizing the fraction of exact exchange). As we discuss in detail in ref 1, B3LYP* was calibrated primarily using spin crossover complexes, which for the great majority of cases have six ligands falling in the middle of the spectrochemical series. DBLOC predicts that the errors in this case will be correlated for all complexes of this type, thus explaining the success of B3LYP* using a single parameter adjustment. However, when complexes with a different distribution of ligands are examined, the errors in B3LYP* are quite large, in fact often larger than those of B3LYP. In contrast, DBLOC displays a uniformly good performance regardless of the distribution of ligands across the spectrochemical series. In ref 7, the DBLOC method was extended to the calculation of redox potentials. 95 small first row transition metal complexes were extracted from the literature, and the redox potential was computed using a continuum solvation model and a thermodynamic cycle. To treat redox potentials, 3 parameters describing the error in the energy required to bring an electron from infinity to reduce the complex are defined, again corresponding to the left, middle, and right of the spectrochemical series. Physically, the redox process is more complex than the transition to an excited spin state, since the total charge on the complex is altered by the redox event. Nevertheless, it is reasonable to postulate that the errors again depend upon the degree of overlap, and energy gap, between the ligand lone pair orbitals and metal d orbitals. The improvement in MUE demonstrated in ref 7, from 0.40 to 0.12 eV, bears out the utility of this hypothesis. A number of unusual ligands require specific corrections based on their differentiated electronic structure. These include radicals, metallocenes, and anionic nitrogens. Finally, a special correction is applied for Ni containing complexes based on their possessing significantly lower multireference character than the wave functions of other transition metal complexes. An extensive discussion of the rationale for each of the spin splitting/redox DBLOC parameters is found in refs 1 and 7. The ratio of training set data points to adjustable parameters in DBLOC development is substantial, which provides some protection against overfitting. However, the data is not sufficiently extensive to have a high degree of confidence in all of the parameters. Consequently, a great deal depends upon the soundness of the physical reasoning that is used to justify the parameters. An investigation of the performance of alternative functionals for these complexes can provide important insights into the validity of the model in this regard, and this is a major reason for carrying out the present study. We will return to this question after presenting the results for the M06 and PBE0 calculations. Finally, we briefly review the composition and analysis of the data sets that will be used below to benchmark performance in computing redox potentials and spin state splittings for first row transition metal complexes. The spin splitting data set was extracted from optical spectra of transition metal complexes displaying a shoulder indicative of a spin forbidden transition which gained intensity via borrowing of oscillator strength through for example a spin−orbit coupling mechanism. As discussed in detail in ref 1, the processing of the raw data to obtain estimates of the spin state splitting is nontrivial but is validated by the success of the DBLOC model in obtaining good agreement between theory and experiment with a very limited number of parameters. Analysis of the M06 and PBE0 results will provide further evidence (one way or another) as to the soundness of the estimated spin state splitting energies. The

roughly the same across many different molecules in which it is embedded, can then readily be implemented. Reference 11 demonstrates that such a scheme is remarkably successful; the mean unsigned error over the G3 data set is reduced to less than 1 kcal/mol (from a B3LYP value of 4.8 kcal/mol), and large outliers are essentially eliminated. Previous papers12−14 have extended the B3LYP-LOC model to the computation of ionization potentials and electron affinities, transition states, and reaction energies, with similar MUEs across substantial data sets. Over the past several years, we have extended the ideas of the LOC approach to address first row transition metal containing systems. In initial work, systematic regularities in the B3LYP errors for transition metal atom ionization potentials and spin state changes were observed. Most prominently, a singly occupied 3d orbital was shown to be substantially overbound, by on the order of 0.40 eV. In a singly occupied orbital, there is no intraorbital nondynamical correlation energy to represent, so the self-interaction term leads to larger errors. A second regularity in the error is in the representation of interaction between d electrons with parallel unpaired spins; in the atom, this error is on the order of 0.1 eV per pair. A similar error was seen in organic atoms, as noted in ref 1. In ref 1 we constructed the initial components of the DBLOC model, studying spin splittings in 57 octahedral transition metal complexes, with the experimental results obtained from analyzing optical transitions with shoulders indicating a spin forbidden component, combined with 14 spin crossover complexes where the splitting can be measured by EPR experiments. The model contains five parameters; the first two follow the observations for atoms outlined above. Moving an electron within either the t2g or eg manifold from a singly occupied d orbital into another singly occupied d orbital, to create a doubly occupied orbital (an excited state), is misestimated by B3LYP by ∼0.40 eV, due to the overstabilization of the singly occupied states noted above for atoms. Interactions between two unpaired spins manifest a similar, although smaller, error, to that seen in atoms, −0.04 eV. The remaining three parameters address the B3LYP error in calculating the splitting between the t2g and eg orbitals. This splitting arises predominantly from the interaction of the eg orbitals with the ligand lone pair orbitals that are facing the metal center. In the simplest picture, six ligand orbitals and one transition metal orbital combine to produce seven new orbitals, six of which are lower in energy than the original ligand orbitals (thus representing bonding orbitals), and the seventh of which (the eg orbital) is higher than the original metal d orbital and thus represents an antibonding orbital. Thus, the error in the splitting between the t2g and eg orbitals will be proportional to the strength of the interaction between the ligand and metal orbitals from which the eg orbital is formed. The strength of this interaction, according to ligand field theory, depends upon the position of the ligand in the spectrochemical series. We postulate that the total error is additive across the 6 ligands in the octahedral complex and define three correction parameters, corresponding to the left, middle, and right locations in the spectrochemical series. These parameters are determined using least-squares fitting across the entire data set. Results for correcting B3LYP were presented in ref 1, demonstrating a large reduction in the MUE as compared to experiment, from 0.4 to 0.1 eV kcal/mol. We also presented in ref 1 a comparison with B3LYP*, an alternative approach to improving the calculation of spin splittings in B3LYP (by C

DOI: 10.1021/acs.jctc.5b00782 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation redox potential data set was taken from the literature, principally reports of cyclic voltammetry experiments, after careful analysis of to eliminate systems with complications (e.g., reactions coupled to the redox process). We believe that this data set provides a good set of benchmarks, although the numerical precision of the data is difficult to estimate due to the heterogeneous nature of the data set (which contains results from many different laboratories) and the lack of direct reporting of error bars in the publications.

Table 2. Complete Set of DBLOC Parameters Optimized for the B3LYP, M06, and PBE0 Functionals DBLOC parameters symbol p ss exlss exmss

III. COMPUTATIONAL METHODS Geometry optimizations and single point calculations using the three functionals under study (B3LYP, M06, and PBE0) were performed using Jaguar version 8.32,14 with the relativistic effective core potential LACV3P, which is a triple-ζ contraction of the LACVP basis set, for metal centers and 6-311G for the rest of the atoms. This basis set was chosen because it was successfully used in developing a localized orbital correction (LOC) model for transition metal atoms and has been shown in the literature to work well in many standard DFT applications to metal containing systems.1 While it cannot be definitively stated that this basis set represents the DFT basis set limit, it has considerable flexibility on both the metal and the surrounding ligands and in our view represents a reasonable compromise between quality and required computational effort. The Poisson−Boltzmann solver in Jaguar2 has been used to provide a continuum dielectric model of the solvents used in the experiments. Starting geometries were taken from the Cambridge Structural Database.3 The amount of exact exchange in the three hybrid functionals studied in this work is 20% for B3LYP, 25% for PBE0, and 27% for M06. DFT calculations on the sorts of complexes investigated in this work can frequently be highly sensitive to the choice of the initial guess wave function. To cope with this possibility, we systematically constructed initial guesses, using the isite keyword in Jaguar, based on ligand-field theory. Additionally, Mulliken population analysis was used to verify the spin states of the metal centers.

exrss relss remss rerss diff rad N−1 Ni

Table 1. Mean Unsigned Errors of Uncorrected and DBLOC Corrected Redox Potentials and Spin Transitionsa

spin forbidden spin cross over redox potentials a

M06

PBE0

B3LYPDBLOC

M06DBLOC

PBE0DBLOC

0.44

0.51

0.41

0.09

0.17

0.09

0.23

0.95

0.54

0.14

0.21

0.14

0.39

0.34

0.40

0.12

0.19

0.14

description

B3LYP

M06

PBE0

create a doubly occupied orbital create a singly occupied-singly occupied interaction t2g → eg excitation (left spectrochemical series) t2g → eg excitation (middle spectrochemical series) t2g → eg excitation (right spectrochemical series) ∞ → t2g reduction (left spectrochemical series) ∞ → t2g reduction (middle spectrochemical series) ∞ → t2g reduction right spectrochemical series) diffuse character of reduced orbital radical character of coordinating atoms anionic nitrogen character of coordinating atoms d7/d8 strong leading order determinants

−0.44 0.05

−0.48 0.02

−0.39 0.04

−0.08

−0.14

−0.11

−0.11

−0.18

−0.13

−0.23

−0.28

−0.21

−0.04

−0.02

−0.04

0.08

0.06

0.08

0.12

0.09

0.12

−0.10 −0.45

−0.07 −0.21

−0.09 −0.39

0.16

0.07

0.11

−0.37

−0.29

−0.49

In what follows, we first present in more detail the results for the spin splitting calculations, followed by a similar presentation of the results for redox potentials. B. Spin Splitting Results. In this study, spin forbidden transition energies of 55 octahedral first row transition metal complexes were computed with M06 and PBE0. The values are compared with a database of experimental spectra collected from the literature and with the results from ref 1 as calculated with B3LYP. In ref 1, a total of 57 complexes were investigated; however, for two of these complexes, we were unable to converge the relevant excited spin states with M06 and/or PBE0. As a simple expedient, we eliminated these two complexes from the data set, leaving 55 complexes in all. We note that for both of the removed complexes, the prior B3LYP results agreed very well with experiment (better than 0.1 eV in both cases), so the excited spin states located by B3LYP are very likely correct. However, as the goal of the paper is to investigate the performance of M06 and PBE0, and their DBLOC-corrected versions, for the prediction of spin splitting and redox energies, we have not further investigated the convergence failure for these two molecules. Detailed results for all 55 cases are presented in the Supporting Information (Table S1). Molecular orbital diagrams constructed in ref 1 show that an electron pairing takes places in t2g → t2g and eg → eg spin forbidden transitions from ground state to excited state. Localized orbital correction terms are applied for doubly occupied orbital and singly occupied-singly occupied interaction. The electronic structure of t2g → eg excitations and eg → t2g de-excitations indicates that a single electron is moved to or from an antibonding eg orbital as opposed to moving an electron within nonbonding t2g or antibonding eg. The energy difference between the t2g and eg sets of orbitals, Δo (ligand field splitting energy), has been found to empirically depend on the coordinating ligands according to the spectrochemical series. Therefore, besides pairing parameters, the localized orbital correction terms for Δo are applied according to the

IV. RESULTS A. Overview. A summary of the performance of the DBLOC method for redox potentials and spin splitting energies is presented in Table 1. Optimized DBLOC parameters for all

B3LYP

values (V)

All values are in eV.

three functionals investigated in this work are presented in Table 2. Optimization of parameters was performed via leastsquares fitting using the data sets described below. Results for individual calculations are presented in the Supporting Information as discussed in what follows. D

DOI: 10.1021/acs.jctc.5b00782 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

excitation and eg → t2g de-excitation databases increases from 0.08 to 0.10 and finally to 0.22 eV for B3LYP as we move to the right of the spectrochemical series corresponding to a larger Δo and thus more strongly interacting ligands. In the same way, the localized orbital correction terms for Δo according to the coordinating ligands of the spectrochemical series from left to right are 0.14, 0.18, and 0.27 eV for M06 and 0.11, 0.12, and 0.20 eV for PBE0. The signs and magnitudes of these corrections, when interpreted in a simple ligand field picture, are in remarkably good correspondence with the analogous range of correction parameters analyzed for metal atoms4 and organic molecules.5 For spin crossover (SCO) complexes, computational results are reported and compared with experiments in Table S2. While pure functionals tend to favor low-spin states, hybrid functionals like B3LYP usually give the high-spin state as the lowest-energy state if the low spin/high spin energy splitting ΔEe is not too large.6 Three hybrid functionals studied in this work, B3LYP, M06, and PBE0, underestimate the spincrossover energies, and the M06 functional is the worst performer. The computed MAEs for the explored functionals are 0.23 eV for B3LYP, 0.95 eV for M06, and 0.54 eV for PBE0. However, the application of DBLOC parameters to SCO complexes gave rise to a huge improvement in the mean unsigned errors for M06 and PBE0 functionals on the basis of the mean unsigned errors after DBLOC, B3LYP, and PBE0 reproduce the experimental values the most accurately. Ground and excited state multiplicities of spin complexes were computed by M06 and PBE0 and compared with B3LYP. The spin multiplicities calculated by the three functionals are generally in very good agreement. The only exception is the excited state multiplicity of the Ni(en)3 complex. While B3LYP and PBE0 calculate the excited state multiplicity 0.71 and 0.61, M06 computed it as 0.00. All three functionals predict the correct ordering of ground and excited state multiplicities in almost every case. C. Redox Potential Results. In ref 7, single-electron reduction half potentials of 95 octahedral first transition metal series complexes binding a diverse set of ligands were calculated at the unrestricted pseudospectral B3LYP/LACV3P level of theory in a continuum solvent. For the present paper, we investigated these same complexes using the M06 and PBE0 functionals. We had difficulty converging three of the complexes to appropriate states with all three functionals, so these three complexes were eliminated from this study, leaving a total of 92 complexes in all. A computational protocol similar to our previous work with spin complexes was used. Redox potentials were determined from adiabatic ΔSCF methods using fully unrestricted pseudospectral M06 and PBE0 free energies of the oxidized and reduced complexes. Experimental and calculated total standard reduction half potentials of all 92 complexes, presented in Table S4, are reported relative to the experimental normal hydrogen electrode (NHE) reference potential of 4.44 V. The previous five spin-splitting parameters (applied without any further parameter adjustment) and a new set of seven redox parameters (fit to the current 92 complex database) are given in Table 2. This provides a single functional form that can be applied to the calculated energies of reactants, and products can differ in the arrangement of d-orbital electrons, as in a spin-splitting reaction, or additionally can differ in the number of electrons, as in redox reaction. First, spin-splitting parameters will be applied to correct the spin-crossover components of the spin-

positions of ligands in the spectrochemical series, as discussed above. In Table 2 above, parameters p, ss, exlss, exmss, and exrss are the DBLOC model parameters required to correct the spin splitting predictions. Sign conventions for the parameters are as shown for either creating an electron pair (p), creating a spin− spin interaction (ss), or for moving an electron from t2g → eg (exlss, exmss, and exrss, depending upon the ligand atoms coordinated to the metal) as the complex transitions from the ground to the excited state. Note that the parameter p is applied per pair, the parameter ss is applied per interaction, while the parameters exlss, exmss, and exrss are applied per coordinating ligand. Oppositely signed parameters must be used if the ground to excited state transition involves the opposite of the listed physical process.1 B.1. t2g → t2g Energetics. For most of the test cases in the data set, two low-lying states are present in the excited state manifold; the high energy state is a more restrictive spin pure state, while the low energy state is spin contaminated with the high-spin ground state. The low energy state corresponds to the state accessed in the experiments, as discussed in detail in ref 1. Mean unsigned errors across the entire data set are 0.49 eV for B3LYP, 0.35 eV for M06, and 0.42 eV for PBE0, respectively. After DBLOC correction, B3LYP and PBE0 have close mean unsigned errors which are 0.10 and 0.11 eV, respectively, while M06 has a larger mean unsigned error, 0.23 eV. B. 2. eg → eg Energetics. Mean unsigned errors are 0.52 eV for B3LYP, 0.61 eV for M06, and 0.48 eV for PBE0. After DBLOC correction, B3LYP and PBE0 have close mean unsigned errors which are 0.08 and 0.08 eV, respectively, while M06 has a large mean unsigned error, 0.14 eV. The average errors for the spin-forbidden transition energies are around 0.43 eV, as seen in Table S1, a value consistent with the t2g → t2g transitions for B3LYP and PBE0. However, the average error for eg → eg transitions is about 0.30 eV more than t2g → t2g transitions for M06. B.3. t2g → eg and eg → t2g Energetics. Mean unsigned errors are 0.30 eV for B3LYP, 0.51 eV for M06, and 0.34 eV for PBE0. After DBLOC correction, B3LYP and PBE0 have very similar mean unsigned errors, which are 0.09 and 0.10 eV, respectively, while M06 has a larger mean absolute error, 0.15 eV. B.4. Summary of the Spin Splitting Results. Overall mean unsigned errors (MUE) for spin forbidden transitions are 0.44 eV for B3LYP, 0.51 eV for M06, and 0.41 eV for PBE0. When DBLOC correction is applied, mean unsigned errors decrease to 0.09 eV for B3LYP, 0.17 eV for M06, and 0.09 eV. Based on the results, mean unsigned errors calculated by B3LYP and PBE0 with respect to experimental values are close to each other in terms of both “overall spin forbidden energies” and “t2g → t2g”, “eg → eg”, “t2g → eg and eg → t2g” spin-splitting energetics. DBLOC parameters, p, ss, exlss, and exrss in Table 2, are found by linear regression. The exmss parameter is found by the application of DBLOC to spin-crossover complexes in Table 2 since there in no middle spectrochemical series ligand in the spin complexes database. The magnitude of parameter, p, which is 0.44 eV for B3LYP, 0.47 eV for M06, and 0.39 eV for PBE0, signifies the importance of errors in the description of nondynamical correlation of occupied d orbitals in metal complexes. The interaction between two singly occupied orbitals gives an estimated correction value of 0.04 eV for B3LYP and PBE0 and a corresponding value of 0.02 eV for M06. The error in t2g → eg E

DOI: 10.1021/acs.jctc.5b00782 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

optimization involving a diverse set of experimental data, including explicitly transition metal containing systems. The optimization strategy is at the other end of the spectrum as compared to PBE0, which does not fit parameters to experimental data. It is therefore not surprising that on average the M06 values for spin splittings and redox potentials differ very substantially from both B3LYP and PBE0. The present results allow us to tentatively compare, in a preliminary fashion, the effectiveness of the various strategies in obtaining transferable results for the properties considered here, spin splittings and redox potentials for coordinatively saturated organometallic species. They also enable an assessment as to the regularities in the errors of all three methods. Systematic errors are more easily correctible than those which do not manifest a discernible pattern. We first discuss the performance of the three uncorrected functionals. Overall, there is not much to choose from as M06 does better for redox potentials but correspondingly worse for spin splittings, particularly for the spin crossover systems. The performances of PBE0 and B3LYP are very comparable. One might have expected M06 to outperform B3LYP and PBE0 on these data sets, as it was trained on some transition metal cases, but the paucity of data of this type in the training set, as well as the significant differences in chemical structures between the training set and the test cases, led to only minimal improvement (perhaps manifested in the improved redox potential results). The DBLOC corrections yield very large improvements for all three functionals. B3LYP-DBLOC and PBE0-DBLOC yield an accuracy that is likely comparable to that of the experimental data, while M06 is 1.5−2× worse than this. The results for redox potentials are particularly interesting; M06 starts out 50% better but ends up 2× worse. We must conclude that the B3LYP and PBE0 functionals have errors concentrated in components that are additive with the chemical nature of the ligand, whereas M06 has a significant component of the error which is more delocalized. Given that the M06 overall functional form diverges substantially from the common framework shared by B3LYP and PBE0, such an effect is certainly plausible (although not predictable a priori, one way or the other). The most interesting results are the close correspondence of the DBLOC correction parameters for all three functionals, including M06 despite the significantly larger residual errors, major fluctuations for selected individual complexes, and very different functional form and parametrization approach. This correspondence clearly establishes that DBLOC is capturing underlying physical regularities in the errors in the hybrid DFT treatment of transition metals, which appears to be consistent across the three different functionals examined here. This confirmation is particularly important for a number of the parameters which were based on a small number of training set examples, such as the corrections for diffuse orbitals, radicals, anionic nitrogens, and nickel containing complexes. For all of these parameters, the PBE0 and B3LYP correction parameters are in good quantitative agreement (within 0.1 eV in all cases), while for M06 the agreement is more qualitative but still reasonable. For the remaining parameters, agreement is semiquantitative in all cases, with a maximum deviation of 0.07 eV and more typically agreement to within 0.02−0.03 eV. Agreement for these parameters is even closer if only PBE0 and B3LYP are considered.

crossover coupled redox reactions. Then, the new set of seven redox parameters are introduced. The relss, remss, and rerss parameters correlate with the spectrochemical series similar to the exlss, exmss, and exrss parameters and describe the functionals’ errors in HOMO/LUMO energetics. The diffuse parameter, diff, is needed when an electron is added to a highly diffuse orbital, and the radical parameter, rad, is used to describe ligands that retain significant doublet character, for example, NO• in the [Mn(II) (CN)5(NO)]−2 and [Fe(II) (CN)5(NO)]−2 complexes. The anionic nitrogen parameter, N−1, is used in complexes which coordinate a nitrogen with significant negative charge, as in the [Cr(III)(en)2(ncs)2]+1 complex. The Ni parameter is used in nickel containing redox couples with fully filled t2g manifolds without spin-crossovers. All Ni complexes in Table S4 need the Ni parameter correction term, as noted above and discussed in detail in ref 7. As in the case of the spin splitting calculations, the spin states of the initial and final states of the various redox processes enumerated in Table S4 display good agreement across all three functionals. This agreement establishes that, at least for the systems investigated here, the underlying description of the electronic structure by all three functionals is qualitatively very similar. The final DBLOC model brings the mean unsigned errors (MUEs) from 0.39 to 0.12 V for B3LYP, 0.34 to 0.19 V for M06, and 0.40 to 0.14 V for PBE0. In Table 2, all parameters are reported in the electrochemical framework and can be used directly to compute standard reduction potentials versus NHE. The values of final DBLOC parameters for B3LYP and PBE0 are very close the each other, and the M06 parameters have similar trends although the numerical values have larger differences. The implications of these results will be considered further below in the Discussion section.

V. DISCUSSION The raw PBE0 values for both spin splittings and redox potentials are in many cases very close (within ∼0.1 eV) to the B3LYP values. However, there are individual cases where there are quite large (absolute) differences; 0.41 eV for spin splitting (excluding spin crossover complexes), 0.34 eV for spin splitting in spin crossover complexes, and 0.59 eV for redox potentials. For M06, the raw values on average have significantly greater differences, with maximum absolute differences of 0.56, 1.10, and 0.86 eV, respectively, for the three quantities enumerated above. The closeness of the B3LYP and PBE0 results is rather surprising in view of the significant differences in the functional form and method of parameter optimization used for B3LYP and PBE0. PBE0 uses the Perdew−Wang (PW) correlation functional,18 while B3LYP employs the Lee−Yang−Parr (LYP) correlation functional;16 among other differences, the former has no adjustable parameters while the latter is fit to the correlation energy of the He atom. Similarly, for PBE0 the PBE exchange functional,18 again containing no parameters fit to experimental data, is used, whereas Becke exchange17 (parametrized to the energy of rare gas atoms) is used for B3LYP. On the other hand, the overall functional form of the two functionals has important similarities (both combine the LDA, a GGA functional, and a fraction of exact exchange, in the range of 20−25%). Overall, PBE0 has no parametrization to experimental data, whereas B3LYP is fit to atomization energies of the G2 data set of Pople and co-workers.19 M06 differs substantially from both of the above functionals in its use of meta-GGA exchange and correlation components and in its use of a large number of fitting parameters and F

DOI: 10.1021/acs.jctc.5b00782 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

■

These observations suggest two future research directions. First, they provide a solid foundation for continuing to refine and expand the DBLOC methodology. Additional benchmark data (from either theory or experiment) will be needed, but once it is available, improved accuracy and coverage of properties for transition metals should be possible. We believe that the DBLOC method, while semiempirical in character due to the required parametrization, can perform a useful function in analyzing and predicting experimental results for transition metal containing systems. At present, achieving accurate results with alternative approaches for the properties discussed here is a challenging task. Second, other approaches to the development of functionals may be able to exploit the systematic nature of the errors we have uncovered here and address them via modifications of functional form. The present results clearly identify where the weaknesses of currently available functionals lie. It is noteworthy that despite the substantial increase in complexity of the M06 functional (which has 22 adjustable parameters), the errors in predicting spin splittings and redox potentials for octahedral transition metal complexes are quantitatively and qualitatively very similar to those of the older B3LYP and PBE0 functionals, which have many fewer parameters and terms. New theoretical approaches are needed to address these deficiencies. The present work provides insight into the nature of these deficiencies across a wide range of example complexes. The data sets assembled herein also constitute test cases for new methods which propose to obtain better results. Finally, we should note that the data sets discussed here do not cover the full range of geometries of first row transition metal containing species. Tetrahedral complexes may require different parameters but would be expected to be amenable to a DBLOC correction scheme. Coordinatively unsaturated complexes, and transition states, may require a more sophisticated approach. Based on our work for organic systems, in which a modified B3LYP-LOC model displays a good ability to improve the agreement of calculated barrier heights with experiment, we are optimistic that such problems can be addressed; however, evaluation using a benchmark data set will be required to establish the validity of this supposition, one way or another. For transition metals, readily interpretable experimental data of this type is difficult to find in the literature; high level computation may solve this problem going forward.

Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.5b00782. Additional data describing experimental and calculated spin-forbidden transitions, spin-crossover transitions, redox potentials for each of the compounds studied computed at the B3LYP-DBLOC, M06-DBLOC, PBE0DBLOC levels of theory, along with ligand field diagrams used to determine the DBLOC corrections (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare the following competing financial interest(s): R. Friesner has a significant financial stake in Schrodinger, Inc., is a consultant to Schrodinger, Inc., and is on the Scientific Advisory Board of Schrodinger, Inc.

■

ACKNOWLEDGMENTS This work was supported in part by the grant from the Department of Energy program through solar photochemistry (DE-FG02-90ER-14162) to R.A.F.

■

REFERENCES

(1) Hughes, T. F.; Friesner, R. A. Correcting Systematic Errors in DFT Spin Splitting Energetics for Transition Metal Complexes. J. Chem. Theory Comput. 2011, 7, 19−32. (2) Jaguar, version 8.3; Schrödinger, LLC: New York, NY, 2014. (3) Allen, F. H. The Cambridge Structural Database: A Quarter of A Million Crystal Structures and Rising. Acta Crystallogr., Sect. B: Struct. Sci. 2002, 58 (3), 380−388. (4) Rinaldo, D.; Tian, L.; Harvey, J. N.; Friesner, R. A. Density Functional Localized Orbital Corrections For Transition Metals. J. Chem. Phys. 2008, 129 (16), 164108. (5) Xu, X. F.; Zhang, W. J.; Tang, M. S.; Truhlar, D. G. Do Practical Standard Coupled Cluster Calculations Agree Better than Kohn-Sham Calculations with Currently Available Functionals When Compared to the Best Available Experimental Data for Dissociation Energies of Bonds to 3d Transition Metals? J. Chem. Theory Comput. 2015, 11 (5), 2036−2052. (6) Reiher, M. Theoretical study of the Fe(phen)2(NCS)2 spincrossover complex with reparametrized density functionals. Inorg. Chem. 2002, 41 (25), 6928−6935. (7) Hughes, T. F.; Friesner, R. A. Development of Accurate DFT Methods for Computing Redox Potentials of Transition Metal Complexes: Results for Model Complexes and Application to Cytochrome P450. J. Chem. Theory Comput. 2012, 8, 442−459. (8) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0Model. J. Chem. Phys. 1999, 110 (13), 6158−6170. (9) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, And Transition Elements: Two New Functionals And Systematic Testing Of Four M06-Class Functionals And 12 Other Functionals. Theor. Chem. Acc. 2008, 120 (1−3), 215−241. (10) Jiang, W.; Deyonker, N. J.; Determan, J. J.; Wilson, A. K. Toward Accurate Theoretical Thermochemistry Of First Row Transition Metal Complexes. J. Phys. Chem. A 2012, 116 (2), 870− 885. (11) Friesner, R. A.; Knoll, E. H.; Yixiang, C. A Localized Orbital Analysis of the Thermochemical Errors in Hybrid Density Functional

VI. CONCLUSION We have evaluated the performance of two widely used functionals, M06 and PBE0, for predicting redox potentials and spin state splittings for large data sets of octahedral complexes of first transition metal series atoms. We find that both methods yield errors that are very similar to those seen in the case of B3LYP, on the order of 0.4 eV MUE. The DBLOC approach successfully corrects the errors in both methods. For PBE0DBLOC, the resulting MUE and parameter values are remarkably similar to those seen for B3LYP, with the MUE close to experimental error. For M06, the MUE is 50−100% larger, but the parameter values still track the B3LYP values closely. The congruence of the parameter sets demonstrates that the physical model for the errors that underlies DBLOC provides a good description of such errors for diverse hybrid functionals. G

DOI: 10.1021/acs.jctc.5b00782 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation Theory: Achieving Chemical Accuracy via a Simple Empirical Correction Scheme. J. Chem. Phys. 2006, 125 (12), 124107. (12) Knoll, E. H.; Friesner, R. A. Localized Orbital Corrections for the Calculation of Ionization Potentials and Electron Affinities in Density Functional Theory. J. Phys. Chem. B 2006, 110 (38), 18787− 18802. (13) Goldfeld, D. A.; Bochevarov, A. D.; Friesner, R. A. Localized orbital corrections applied to thermochemical errors in density functional theory: The role of basis set and application to molecular reactions. J. Chem. Phys. 2008, 129 (21), 214105. (14) Bochevarov, A. D.; Harder, E.; Hughes, T. F.; Greenwood, J. R.; Braden, D. A.; Philipp, D. M.; Rinaldo, D.; Halls, M. D.; Zhang, J.; Friesner, R. A. Jaguar: A High-Performance Quantum Chemistry Software Program With Strengths In Life And Materials Sciences. Int. J. Quantum Chem. 2013, 113 (18), 2110−2142. (15) Hall, M. L.; Goldfeld, D. A.; Bochevarov, A. D.; Friesner, R. A. Localized Orbital Corrections for the Barrier Heights in Density Functional Theory. J. Chem. Theory Comput. 2009, 5 (11), 2996− 3009. (16) Lee, C. T.; Yang, W. T.; Parr, R. G. Development of the ColleSalvetti Correlation-Energy Formula into a Functional of the ElectronDensity. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37 (2), 785− 789. (17) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic-Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38 (6), 3098−3100. (18) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77 (18), 3865. (19) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Pople, J. A. Assessment of Gaussian-2 and density functional theories for the computation of enthalpies of formation. J. Chem. Phys. 1997, 106 (3), 1063−1079. (20) Jerome, S. V.; Hughes, T. F.; Friesner, R. A. Accurate pK(a) Prediction In First-Row Hexaaqua Transition Metal Complexes Using The B3LYP-DBLOC Method. J. Phys. Chem. B 2014, 118 (28), 8008− 8016. (21) Hughes, T. F.; Harvey, J. N.; Friesner, R. A. A B3LYP-DBLOC Empirical Correction Scheme For Ligand Removal Enthalpies Of Transition Metal Complexes: Parameterization Against Experimental And CCSD(T)-F12 Heats Of Formation. Phys. Chem. Chem. Phys. 2012, 14 (21), 7724−7738. (22) Polo, V.; Kraka, E.; Cremer, D. Electron correlation and the selfinteraction error of density functional theory. Mol. Phys. 2002, 100 (11), 1771−1790. (23) Gräfenstein, J.; Kraka, E.; Cremer, D. Effect of the selfinteraction error for three-electron bonds: On the development of new exchange-correlation functionals. Phys. Chem. Chem. Phys. 2004, 6 (6), 1096−1112. (24) Gräfenstein, J.; Kraka, E.; Cremer, D. The impact of the selfinteraction error on the density functional theory description of dissociating radical cations: Ionic and covalent dissociation limits. J. Chem. Phys. 2004, 120 (2), 524−539.

H

DOI: 10.1021/acs.jctc.5b00782 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX