Ligand Additivity and Divergent Trends in Two Types of Delocalization Errors from Approximate Density Functional Theory

Despite its widespread use, the predictive accuracy of density functional theory (DFT) is hampered by delocalization errors, especially for correlated systems such as transition-metal complexes. Two complementary tuning strategies have been developed to reduce delocalization error: eliminating the global curvature with respect to charge addition or removal, and computing a linear response Hubbard U as a measure of local curvature at the metal center at fixed charge and applying it to the transition-metal complex in a DFT+U framework. We investigate the relationship between the two measures of delocalization error as we manipulate the ligand field strength by varying the number of strong-field ligands in a series of heteroleptic complexes or by geometrically constraining the metal-ligand bond length in homoleptic octahedral complexes. We show that across these sets of complexes with varying ligand fields, an inverse relationship generally exists between global and local curvatures. We find that effects of ligand substitution on both measures of delocalization are typically additive, but the two quantities seldom coincide. The observation of ligand additivity suggests opportunities for evaluating errors on homoleptic complexes to infer corrections for lower-symmetry complexes.

local curvature. This local curvature is defined as the second derivative of the energy with respect to the removal or addition of charge from a localized subshell ( ,-. ) to the rest of the system at constant charge. That is, where ,-. is the occupation of the subshell defined by the quantum numbers , of atom . The local curvature is calculated with the motivation that its functional form should eliminate the convexity of the energy of an atom as fractional charge is added to or removed from it. The elimination of convexity should subsequently recover piecewise linearity under certain conditions 29,30 .
Although elimination of local curvature and elimination of global curvature have similar aims, the two curvatures represent different properties of the system. In fact, it has been shown that the two curvature values rarely coincide in representative homoleptic transition-metal complexes 31 and instead an inverse relationship has been suggested. 31 We now investigate these two measures of DE in TMCs to determine if similar relationships exist when the ligand field is manipulated by ligand substitution in heteroleptic complexes or by distortion of the metal-ligand bond length in homoleptic complexes. We quantify the extent of ligand additivity 32 in these energetic global and local measures of DE, which has been demonstrated to enable force fields 33 to predict DFT-derived properties. We show that while local and global measures of DE are indeed divergent, concepts of ligand additivity generally hold, suggesting the possibility of quantifying both errors and computing corrections on high-symmetry or small, representative complexes.
We first assembled octahedral complexes where we smoothly vary the ligand field from weak to strong by replacing water or ammonia ligands with carbonyl ligands in heteroleptic complexes. Specifically, we started from a hexa-aqua/hexa-amine complex for our weak-field limit and replaced the ligands one by one with carbonyls until we reached a hexa-carbonyl complex representing the strong-field limit ( Figure 1). We then computed properties of heteroleptic complexes between these two limits, with structures obtained from prior work 34 . All stoichiometries were studied except for the fac and mer complexes with three of each ligand type due to their omission from the complexes studied in prior work 34 . Linear response and curvature calculations were conducted on geometry-optimized structures in both high-spin (HS) and lowspin (LS) states for Mn(II/III) and Fe(II/III). This corresponds to a HS quintet and LS singlet for  To validate our choice of heteroleptic complexes as a way to interpolate between limits of ligand field strength, we analyzed the average metal-ligand bond length as a proxy for the shift in ligand field. As more carbonyl ligands are added to an Fe(II) complex, the LS complex bond lengths decrease, whereas for the HS complex, the average bond length increases ( Figure 2).
Although the increase in HS bond lengths is surprising, this increase in the HS complex is due to more asymmetry in the bond lengths for this complex in comparison to the LS case ( Figure 2 and Supporting Information Figure S1). Importantly, the bond lengths are smoothly varying, verifying the suitability of this model for interpolation between ligand field limits. We next computed linear response U for Fe in a +2.5 oxidation state to balance between the Fe(II) and Fe(III) limits over which global curvature is evaluated. This oxidation state choice was previously shown 31 to provide a roughly average value of those (i.e., Fe(II) and Fe(III)) two limits. We observe a monotonic relationship between nominal ligand field strength and linear response U evaluated for Fe in a +2.5 charge state, with stronger ligand fields leading to higher linear response U values, especially for the LS complex ( Figure 2). Throughout the transition from weak to strong ligand field, the HS complex U values are consistently lower than those for the LS state, which is intuitive because hybridization, which is strongly linked to delocalization, is expected to be weaker in the HS state. This effect is significant, as the HS complex U only increases by about 0.3 eV, but the LS complex linear response U increases by nearly 3 eV as the ligand field is varied from weak to strong ( Figure 2 Figure S2). It has been observed that global curvature values and linear response U values are seldom identical in transition-metal complexes 31 . If the two measures were comparable, then application of a single U in DFT+U could simultaneously eliminate both local and global curvature. A secondary concern is that when aiming to compare two different spin states, either local or global measures of curvature might differ for the states being compared, which poses challenges because energies from DFT+U calculations are only comparable when the same value of U is applied to all points. We consider our data set in light of these potential challenges. The linear response U and global curvature are most comparable across both LS and HS states in the strongfield limit, meaning that a single value for the curvature correction may correct both types of delocalization error in a balanced manner for both spin states in this regime.
An inverse relationship between global curvature and the linear response U has previously been shown 31 by comparing homoleptic complexes of strong-and weak-field ligands, i.e, for hexa-carbonyl, hexa-amine and hexa-aqua complexes. As we increase the ligand field strength in this work by adjusting the number of strong-field ligands in a series of heteroleptic complexes, we observe a consistent negative correlation between the two measures of curvature, and that curvature values for intermediate ligand fields ligands reside in between the two homoleptic endpoints. That is, from strong-to weak-field complexes the global curvature decreases while the linear response U increases ( Figure 3). Because ammonia and water are comparably weak field ligands, we find that exchanging either ligand with carbonyl has a similar effect on global curvature and linear response U (i.e., local curvature). However, we observe that the relationship between global and local curvature depends strongly on the electron configuration (i.e., spin state) ( Figure 3 and Supporting Information Figure S6). Within each spin state (i.e., high-or low-spin), increased ligand field strengths lead to decreased global curvatures and higher linear response U values. The difference between the spin states is evident from the differences in the slopes of the linear fits for each spin state ( Figure 3). In the HS state, the linear response U changes half as fast as the curvature (i.e., slope of -2.07 and R 2 of 0.97) as the field strength is tuned. In the LS state, however, the reverse is true, with a change in linear response U that changes twice as much as the curvature (i.e., slope of -0.5 and R 2 of 0.77). There is also a greater variation in the data for the HS state. The two linear fits coincide if we extrapolate past the weak-field limits defined by our homoleptic complexes. In both cases, the lines approach a low linear response U (< 3 eV) and a high curvature value (> 8 eV) at the weak-field limit. As ligand field strength increases, the differences in curvature and linear response U values for the different spin states increase ( Figure   3). This divergence can be related to differences in the strength of metal-ligand bonding, observed through the bond lengths. The bond lengths for the HS and LS states are closest in value in the case of weak crystal field, and the differences increase with increasing ligand field strength ( Figure 2). Thus, the weak-field complexes have more spin-independent bonding covalency in comparison to the stronger-field complexes, as indicated by differences in bond lengths. Our observations do not depend strongly on the metal identity, as we observe similar results for Mn compounds (Supporting Information Figure S5 and Table S1).
We next isolated ligand field strength effects from dependence on the number of atoms or chemical identity of the ligand in evaluating linear response U and the global curvature. To do so, we constructed a homoleptic iron hexa-carbonyl complex, and gradually increased the Fe-C  Figure S7). Thus, our sampled distances extend beyond the longest equilibrium bond and are comparable to the shortest strongfield Fe-C bond length. As we decrease the ligand field strength by elongating the Fe-C bond, we also decrease the linear response U for both HS and LS states, consistent with our observations on the heteroleptic compounds ( Figure 4). The LS linear response U values are again consistently higher than the HS values ( Figure 4). The LS state exhibits a greater variation in Hubbard U with Fe-C distance than the HS state; we observe a 1 eV decrease in U value as the Fe-C bond length increases from 2.2 to 2.5 Å, in comparison to a 0.5 eV increase for the equivalent HS case.
Despite the qualitative consistency between trends observed for heteroleptic complexes and complexes where the Fe-ligand bond length is manipulated, the variation in U with Fe-C bond length variations is smaller than was observed upon substitution of weak-field ligands with carbonyl ligands. Specifically, replacing all water ligands with carbonyls (i.e., x = 0 to 6 in Fe(CO)x(H2O)6-x) exhibits a change in U that is at least twice as large ( Figure 2). The difference is most likely not due to ligand field strength variation. Because these complexes maintain octahedral symmetry, we can use the d orbital splitting as a direct measure of the ligand field strength instead of relying on the proxy of the bond length. We observe this splitting to decrease by ca. 2 eV with Fe-C bond stretching, which is nearly twice that observed with heteroleptic substitutions of carbonyl for weak-field ligands (Supporting Information Figure S8). That is, while changing ligand field through chemical substitution has a significant effect on linear response U, similar changes in ligand field strength by modulating metal-ligand bond has a smaller effect on the linear response U. The linear response U is therefore not only affected by the strength of the ligand field but also by the character of the ligand such as its ability to form p bonds or its polarizability.
Variations in global curvature are smaller than variations in the linear response U as Fe-C bonds are stretched ( Figure 4). The global curvature in our model system is consistently higher in the HS state than the LS state, in agreement with observations on heteroleptic complexes.
Nevertheless, the HS state has an increase in curvature with decreasing ligand field, whereas the LS state curvature is flat or decreases slightly (Figure 4). We observe less than a 0.5 eV increase in curvature for the HS complex with stretched Fe-C bonds compared to a 1.5 eV variation in the set of Fe(CO)x(H2O)6-x complexes ( Figure 2). This is most likely the result of the frontier orbital energies, which depend on the chemical composition of the entire system, changing less when we only vary bond length and do not vary ligand chemistry. In the LS case in particular, decreasing the ligand field by stretching the Fe-C bond in our homoleptic model complex does not result in the same increase in the global curvature that we observed when we decreased the ligand field by replacing strong-field ligands with weak-field ligands. These differences can be rationalized by noting that global curvature is more sensitive to the nature of the ligands in the complex than to the strength of the metal-ligand bond. To study delocalization error, we evaluated or approximated (i.e., from eigenvalues) fractional-charge properties between the +2 and +3 metal oxidation states in this work for global curvature and linear response U. We chose high-spin and low-spin conventions in accordance with prior work as follows. The Mn(II)/Mn(III) spin multiplicity decreases from sextet to quintet in the high-spin configuration and from doublet to singlet in the low-spin case, as a majority-spin electron is removed (Supporting Information Table S2). The Fe(II)/Fe(III) spin multiplicity increases from quintet to sextet and from singlet to doublet in high-and low-spin configurations as a minority-spin electron is removed from the system (Supporting Information Table S2). The HOMO and LUMO, which are used in the average curvature expression, were chosen as the levels from which an electron is removed or to which an electron is added according to these conventions, respectively.
To compute global curvature, we used the average curvature expression as in eq.    Variation in the curvature is in this complex is small compared to other systems: 0.5 eV here vs. 1 eV in the corresponding complex with Fe for example in Figure S2. The single-point calculation for the Mn +3 (NH3)6 endpoint did not converge. As a result, the global curvature could not be calculated for this point. Therefore, this complex was omitted from the comparison. In addition, the cis Mn(CO)2(NH3)4 complex is not displayed because geometry optimizations did not succeed.   Table S2). For simplicity, only the trans structure of the Fe(CO) " (H " O) & and Fe(CO) " (NH % ) & complexes are shown in the figure. Mn hexa-aqua is not displayed in this figure due to convergence issues in the geometry optimization of the complex. Mn hexa-ammine is not shown either because the single-point calculation of the +3-oxidation state did not succeed.    ). The d-splitting energy was calculated from the eigenvalues in the LS state, i.e., the energy of the LUMO minus the energy of the HOMO. The filled circles correspond to Fe in a +2 oxidation state, and the empty circles correspond to Fe in a +3 oxidation state, the two endpoints by which the average curvature is calculated. The arrows in the figures represents the direction of increasing field strength. Table S2. Formal oxidation and spin multiplicities (2S+1) considered in this work. The first number is the multiplicity of the low-spin state and the second number is the multiplicity of the high-spin state. For each metal, global curvatures were obtained using the oxidation states as the energy deviation end points, both in high-spin and low-spin states. Linear response calculations were performed on each oxidation state in both spin states.