Assessing the Calculation of Exchange Coupling Constants and Spin Crossover Gaps Using the Approximate Projection Model to Improve Density Function Calculations

Assessing the Calculation of Exchange Coupling Constants and Gaps the Approximate Improve This work evaluates the quality of exchange coupling constant and spin crossover gap calculations using density functional theory corrected by the Approximate Projection model. Results show that improvements using the Approximate Projection model range from modest to significant. This study demonstrates that, at least for the class of systems examined here, spin-projection generally improves the quality of density functional theory calculations of J-coupling constants and spin crossover gaps. Furthermore, it is shown that spin-projection can be important for both geometry optimization and energy evaluations. The Approximate Project model provides an affordable and practical approach for effectively correcting spin-contamination errors in molecular exchange coupling constant and spin crossover gap calculations. Abstract This work evaluates the quality of exchange coupling constant and spin crossover gap calculations using density functional theory corrected by the Approximate Projection model. Results show that improvements using the Approximate Projection model range from modest to signiﬁcant. This study demonstrates that, at least for the class of systems examined here, spin-projection generally improves the quality of density functional theory calculations of J -coupling constants and spin crossover gaps. Furthermore, it is shown that spin-projection can be important for both geometry optimization and energy evaluations. The Approximate Project model provides an aﬀordable and practical approach for eﬀectively correcting spin-contamination errors in molecular exchange coupling constant and spin crossover gap calculations.


Introduction
Transition metal complexes play an important role in a wide range of chemistries including catalysis, 1 magnetic materials, 2 and biological systems. 3 Over the past several decades, computational models and methods for transition metals have become increasingly efficient and accurate; meanwhile computers have also become orders of magnitude faster. All these factors have made computational study of transition metal chemistry more wide spread, reliable, and predictive. While computational study has become an important driving force in transition metal chemistry, [3][4][5] it can be quite difficult to achieve chemical accuracy (error ∼ 1 kcal/mol) using available models due to the complex nature of transition metals. As a result, routine modeling transition metal systems often employs Density Functional Theory (DFT). 6 As a single-reference method, DFT treatment of transition metals lacks a description of static correlation energy, often yielding unsatisfying results and limiting the scope of potential applications. 7 An alternative and computationally feasible approach that has seen widespread application is broken-symmetry DFT. This scheme employs an unrestricted Kohn-Sham determinant that breaks spin symmetry in the density. Analysis of the broken-symmetry determinant shows that it is a mixture of different spin-pure configurations. Errors associated with this spin contamination can significantly change the energy and distort the potential energy surface (PES). These errors are particularly prevalent in geometrically strained and spin-polarized systems, such as transition structures and transition metal clusters. [8][9][10] Nevertheless, in certain cases symmetry breaking yields an improved treatment of static correlation at single-reference cost, though static and dynamic correlation are can become poorly balanced in such cases. 11 To correct spin contamination error, spin projection methods can be used to remove the contaminating spin states from a broken symmetry determinant. A (approximate or partially) projected determinant will be free of contamination and still retain a fair degree of static correlation contributions. This scheme is often referred to as projection-after-variation (PAV), and it can recover static correlation without explicitly resorting to multi-reference methods while maintaining mean-field cost. In contrast, variation-after-projection (VAP) methods variationally optimize a projected determinant, which will obtain not only static correlation, but some dynamic correlation as well. Examples of VAP schemes include the extended Hartree-Fock method of Löwdin and, more recently, the projected Hatree-Fock (PHF) method developed recently by  et. al.. PAV and VAP results are not directly comparable and the formal weaknesses of PAV are well documented. 15 However, PAV allows a clean separation of correlation types, which may present a relatively straightforward pathway for ongoing theory and methods developments.
Approximate Projection (AP), a projection method proposed by Yamaguchi and coworkers, 16 has shown its usefulness for energy calculations of transition metal complexes and diradicals. [17][18][19] However, to the best of our knowledge, there has not been a systematic investigation on the effectiveness of AP in the study of magnetic properties in transition metal systems. To address this need, we applied AP to two kinds of problems: 1) calculation of exchange coupling constants (J-couplings) for a set of binuclear transition metal complexes and 2) spin crossover gaps for mononuclear transition metal complexes.
This work also investigates geometry optimization effects, i.e., how geometry optimization affects J-coupling prediction. It has been found that full relaxation of the crystal structure in gas phase will slightly worsen J-coupling results relative to experiment. 20 It is also expected that projection will change geometry minima since it modifies energy as well as the energy landscape. Therefore, it is beneficial to know if those changes are meaningful for prediction of physical properties. In recent years, our group has extended the use of AP by developing and implementing efficient first and second derivatives, 21,22 making it efficient to obtain optimized structures on AP-corrected potential energy surfaces.
In this report, we use AP to test how projection methods affect the prediction of magnetic properties with special focus on the importance of geometry optimization in such calculations. Specifically, we compare the AP predicted J-coupling and spin crossover gap values with those obtained by conventional broken-symmetry DFT and experiment. We also compare results obtained from different geometries and using different approximate density functionals. We show that spin projection can alter PESs and that addressing this consideration in geometry optimizations is beneficial for computational studies of systems affected by spin contamination.

Computational Details
All electronic structure calculations were carried out using a local development version of Gaussian. 23 Calculations using the AP model employed our lab's local implementations of analytic first-and second-derivative programs. 11,21,22 The examination of J-coupling calculations has been carried out using five different functionals -B3LYP, 24,25 LC-ωPBE, 26,27 CAM-B3LYP, 28 B3PW91 29 and ωB97XD. 30 Previous studies have shown that B3LYP and LC-ωPBE yield the best J-couplings among several density functionals. 31,32 All J-coupling calculations were performed with Ahlrich's all electron triple-ζ valence plus polarization for the metal centers 33 and Ahlrich's double-ζ valence plus polarization basis sets for non-metal atoms. 34 A pruned integration grid of 75 radial shells and 302 angular points per shell was used in DFT calculations (corresponding to the so-called "Ultrafine" grid option in Gaussian 16). Broken symmetry solutions were generated by using a fragment guess approach to ensure spin-localized initial Kohn-Sham (KS) determinant guesses. All converged KS solutions were verified by stability analysis. 35,36 To explore the impact of spin-projection on the types of calculations being explored here, we have chosen to use the AP method as a proxy for other spin-projection models. Indeed, recent work from our lab has demonstrated the conditions under which such approaches are rigorously equivalent and those cases for which such models are numerically similar. 11 AP calculations require both a broken-symmetry state and of the high-spin state. Those results are used to assemble an AP energy according to Subscriptions low spin (LS) and high spin (HS) refer to the corresponding ideal spin-pure states. AP geometries were obtained by standard geometry optimization methods. Similar to UDFT geometry, the phrase "AP geometry" will be used to refer to the optimized ground state structure of the complex. J was computed for each combination of compound, functional, and geometry according to Eqs. (4), (5), (6).
For each approximate density functional, three different geometries were used for calculating J-couplings. The first approach, and perhaps most widely used methodology in similar benchmark studies, employs the known crystal structures unaltered from their original literature report. The second source of geometries was the minimum energy structures on the HS or broken-symmetry PESs. The third set of geometries used minimum energy structures located on the AP PESs. In the discussion below, the phrase "UDFT geometry" refers to the lower-energy structure of the HS and broken-symmetry results for each complex, i.e., this refers to a broken-symmetry results when the ground state species is anti-ferromagnetic and a HS results when the ground state is ferromagnetic.
Our tests evaluating spin crossover gaps were carried out using a test set previously used by Hughes

Results and Discussion
As discussed above, we have explored the effect of spin projection on the quality of calculated J-coupling values and spin crossover gaps for two commonly used test sets. A particular interest in this work is the dependence of such calculations on the geometries one employs.
Specifically, we have used our in-house analytic AP first-and second-derivative codes to examine the extent to which spin contamination degrades the quality of the potential energy surface, predicted minimum energy structures, and ultimately the calculated (vertical) gaps between spin-states. This section first reports the results from our tests on a well characterized and widely used J-coupling parameters corresponding to a set of nine transition metal complexes. Next, we examine this effect on spin crossover gaps for a set of 65 complexes recently compiled and reported by Hughes and Friesner. 37 In the final part of this section we discuss general implications of these results.

J-Coupling Tests
Our examination of the role of geometry optimization and spin projection on predicting Jcouplings used the test set constructed by Peralta and Rudra. 31,40 This test set consists of  Table S1).
where J ij is the coupling constant. 42 S i and S j represent theŜ z spin operator for each spin site. Positive J corresponds to FM spin coupling and a HS ground state, while negative J indicates a preferred AF spin coupling and a LS ground state. The magnitude of J indicates the energy gap between two spin states. Experimentally, J is obtained using techniques such as electron paramagnetic resonance (EPR). 43 Using computation, J can be estimated by relating it to the energy difference between HS and LS states described by the Heisenberg Hamiltonian, according to where E B is the energy of the broken-symmetry state. This equation is commonly referred to as the non-projected (NP) approach as it does not involve spin projection. 32 The deficiency of the NP method is that the LS energy typically involves a spin contaminated KS determinant. Addressing this concern, Dai and Whangbo 44 proposed a spinprojected (SP) approach by mapping the broken-symmetry and HS states onto the Heisenberg Hamiltonian, and derived the following expression for J, where E B is the energy from a broken-symmetry calculation. Eq. (4) corrects for the spin contamination of the broken-symmetry result using an approximation for spin projection and will yield a J whose magnitude is necessarily greater than that computed using Eq. (4).
The validity of SP and NP approaches depends on the magnitude of spin coupling. In the weak coupling regime, where the overlap between two interacting magnetic orbitals a and b is close to zero, i.e. a|b ≈ 0, SP should be chosen. NP tends to perform better in the strong coupling regime, where a|b ≈ 1. However, in practice, DFT tends to overestimate the magnitude of J-couplings. 45 NP tends to underestimate the magnitude of J, and when used with DFT, gives better results even in the weak coupling regime. 46,47 Furthermore, as we will show, NP necessarily has smaller errors than SP (vide infra).
The AP model smoothly unifies the two limiting cases treated by the SP and NP approaches. The AP model employs the dominant HS contaminant, finds its weight in the broken-symmetry state, removes it from the contaminated energy, and renormalizes the en-ergy to yield a corrected low-spin energy. The exchange coupling parameter in this model is given by where E AP is the AP-corrected LS energy and Ŝ 2 LS is the ideal spin-squared expectation value of the pure LS state. In the weak coupling regime, when the broken symmetry solution exhibits a large degree of spin polarization, AP reduces to SP (Eq.5). In the strong coupling regime, S 2 B approaches S 2 LS and AP reduces to NP. Our study began by carrying geometry optimizations of the nine complexes comprising the J-coupling benchmark set used by Peralta and others. 20,31,45,48 Initial geometries were take from the crystal structures used in prior studies. Minimum energy structures were located on the HS, broken-symmetry LS, and AP corrected low-spin PESs.
To quantify the extent of geometric relaxation resulting from these optimization calculations, the root-mean-square-difference-at-maximum-overlap values for the crystal structure, HS, and broken-symmetry LS minimum energy geometries were computed relative to corresponding optimized AP geometries. These RMSD values are shown in Table 1. As one might expect, meaningful geometric relaxation is observed when optimizing structures in the gas phase from initial crystal structures. The average RMSD of the crystal structures relative to the APgeometries was 0.48Å. The largest internal coordinate displacements of bond lengths, angles, and dihedrals involving a metal center averaged over the set of nine compounds were 0.09Å, 12.6 • , and 34.2 • , measured between crystal and AP geometries. On the other hand, HS and broken-symmetry LS geometries differ with AP geometry only slightly (RMSD ∼ 0.01 − 0.02Å), which suggests that AP energy correction only slightly changes the potential energy surface, in agreement with prior findings. 49 This observation indicates that -at least for the test set considered here -the geometries are not very sensitive to the presence of spin contamination.
For example, Fig. 2 shows the geometry change of compound 3 from crystal geometry to the AP optimized geometry. Overall, the three-dimensional structure remains qualitatively Table 1: RMSD (Å) of crystal, optimized HS and optimized broken-symmetry LS structures relative to corresponding optimized AP geometries. unchanged, but evaluation of specific internal coordinates shows a few key changes. The most notable change is in the copper-nitrogen bond distance labeled in Fig. 2, which decreases during geometry optimization on the AP PES by more than 0.2 Å. The complete set of geometric parameter changes due to geometry optimization are given in Table S2). We next investigated the quality of calculated J-coupling parameters by employing optimized AP structures. As described earlier, we considered three different approaches for calculating J-couplings -NP (Eq.(4)), SP (Eq.(5)), and AP (Eq.(6)). Table 2 gives Jcouplings calculated by Eqs. (4)-(6) on AP geometries using the B3LYP and LC-ωPBE approximate density functionals. The table also includes experimentally determined J values and calculated expectation values of the S 2 operator for each compound. Results using other approximate functionals and geometries (crystal structures, and HS and broken-symmetry In all cases, calculated J-coupling values had the correct sign, indicating that hybrid functionals (with or without range separation included) predict the correct ground state. In addition, NP J-couplings are always smaller in magnitudes than those obtained by AP. This result is in line with the suggestion that NP counteracts the overestimation of J-couplings by DFT. [50][51][52] Comparing NP and AP results with experimental data, it is clear that the non-projection approach yields better agreement with experiment. Impotantly, it should be kept in mind that the magnitudes of the J-couplings in these systems are much smaller than the typical error limits associated with DFT calculations. Assuming a normal distribution for this error in DFT energy calculations, the calculated HS and LS energies can be related to the true energies according to, and E calc HS = E true HS + X HS (8) where X LS and X HS are two scalar values taken from a normal distribution with deviation σ. Invoking Eq.(4) we obtain and ∆X = X HS − X LS . ∆X is the difference between the two normally distributed variates and thus follows a normal distribution of zero mean and 2σ variance. Thus, the NP error, NP , is given by Similarly, the AP J-coupling error, AP , can determined by using Eq.(6) and is given by J true AP − J exp is the difference between the J value calculated using the correct (or true) LS and HS energies with the AP scheme and the experimentally measured J. This value is of similar magnitude to the experimental J value, if not much smaller. The same observation holds for J true NP − J exp . If one assumes these two parameters are given by the same small value ∆J, e.g. ∼ 100 cm -1 (0.286 kcal/mol), and it is further assumed that σ is 1 kcal/mol, and S 1 and S 2 are assumed to be 0.5 (as is the case for most singlet-triplet splittings) we have and Within the range ∆X < 0 ∪ ∆X > 4 3 ∆J, NP < AP is always true. In a normal distribution, the probability of ∆X being in that range is calculated to be 0.9244. For other fair assumptions of ∆X, ∆J, S 1 and S 2 values, the probability of NP < AP remains close to 1 as long as ∆J is significantly smaller than ∆X. In other words, if a calculated value is normally distributed around the true value with a deviation well beyond its magnitude, a method that consistently results in values of smaller magnitude will yield statistically better results.
The geometry dependence of J-couplings was further examined by calculating NP and AP coupling values using crystal geometries and those found from geometry optimizations on the UDFT and AP PESs. As a representative case, Table 3 reports the J values calculated by LC-ωPBE. Overall, J-coupling dependence on geometry is very small. Relative to the crystal structure results, UDFT geometries resulted in slightly improved J-couplings with the NP method but slightly degraded results using the AP method. Using AP geometries, on the other hand, yielded somewhat improved results relative to UDFT geometry by both AP and NP methods, especially when compound 6 is excluded. Figure 3 shows MAE's of calculated J-couplings for four different approximate density functionals using three difference types of geometries and evaluating J using both AP and NP approaches. In agreement with previous work by Phillips et al., we found that range separated approximate functionals are found to be superior to global hybrid functionals in predicting J-couplings. 48 We note, that that work suggests that the improvement of such results using range separated functionals (relative to global hybrids) is not an effect of range separation. Instead, it appears that the improved agreement with experiment is die to their inclusion of more HF exchange. Our results also show that using geometries optimized on AP PESs, J-coupling values calculated with range separated functionals reduced the MAE by 50-60% (excluding complex 6) with both AP and NP methods, relative to B3LYP results. The  three range-corrected functionals behave similarly. Geometry dependence, although small, can be seen clearly where AP and UDFT geometries give better agreement with experiment than crystal structure with the range separated functionals considered.

Spin Crossover Gaps
As a second test of the impact of spin-contamination on geometry and spin-state energy differences, we explored the role of spin-projection in calculating spin crossover gaps in a broad set of transition metal compounds. The set of spin crossover species was used by Hughes and Friesner in the development of their DBLOC model. 37 This test set includes 65 octahedral first-row mononuclear transition metal complexes, including V, Cr, Mn, Ni, Fe, Co with oxidation states ranging from +2 to +4, total charge ranging from -4 to +4, multiplicity ranging from 1 to 6, and total number of valence d-shell electrons ranging from 2 to 8. The set is divided into several subgroups (transition types) based on the number of valence electrons and transition type (see Table 4). For example, for each member of the t 2g -t 2g d 3 transition type group a t 2g electron is excited to another t 2g orbital.
To determine the dependence on geometry optimization in spin crossover calculations, spin crossover energies were calculated for each compound in the test set using the same three types of geometries as in the exchange coupling calculations described above: (1) crystal structures; (2) * These two groups are experimentally observed to have very small spin crossover gaps. Despite the overall improvement from UDFT to UDFT//AP, no significant difference between them is observed in the t 2g -e g (d 6 ), *t 2g -e g (d 6 ) and **t 2g -e g (d 7 ) transition types.
This can be understood by considering the AP α parameter defined by Eq. (2), which represents the degree of contamination in a broken symmetry determinant. The closer α is to one, the less contaminated a solution. Therefore, AP will have little (or no) effect on spin crossover gaps when α is close to 1.0.
On the other hand, AP correction dramatically reduced the errors of all α > 1 cases, , t 2g -t 2g (d 2 ) and e g -e g (d 8 ), from 11.24 to 3.03 kcal/mol (see last row of Table 5). Using AP these transition types all yielded relatively accurate AP behaved particularly well for group t 2g -t 2g (d 3 ), which is shown in detail in Table 6.
Crystal and UDFT values were only slightly different, meaning that geometry changes impose similar effects on HS and LS energies and thus the change in spin crossover gap is negligible.
The MAE obtained on Crystal structures is 12.83 kcal/mol, which was significantly improved to 1.21 kcal/mol by using AP. In this case, only LS states were affected by AP because HS states were not spin contaminated. Therefore AP had a substantial impact on the energy gaps. AP//AP values were calculated on AP geometry instead of UDFT geometry, but no significant improvement in spin crossover gap was seen. Both UDFT//AP and AP//AP Calculations without solvent effect were also performed and yielded very similar results.
The MAE for UDFT//AP and AP//AP model chemistries were the same within two decimal places, and the MAE of α > 1 cases of these two methods were 3.04 and 2.98 kcal/mol (See Table S9). From our results, solvent effect for this test set does not seem to be significant.
Overall, AP provides reasonably good results for spin crossover gaps predicted by UDFT.

Conclusions
This work explored the effect of spin projection on the quality of calculated J-coupling values and spin crossover gaps in transition metal complexes. Furthermore, the dependence of such calculations on geometry has been studied. A widely used test set of J-couplings for nine transition metal complexes was studied with and without applying spin-projection corrections. Our results show that spin-projection does not often result in large geometric changes in the minimum energy structure relative to optimized spin-contaminated unrestricted DFT geometries. However, the use of geometries optimized with the APmethod in J-coupling computations using either the spin-projected or non-projected approaches generally results in better agreement with experiment than using unoptimized crystal structure geometries.
A second set of benchmark tests involved a set of 65 transition metal complex spin crossover gaps. Projection methods were found to be effective for calculating spin crossover energy gaps. Focusing on a subset of 43 members of this set where the AP model is expected to improve the quality of mean-field calculations showed that AP geometry optimization and energy evaluations leads to very good agreement with experiment. Importantly, these results are similar in quality to the DBLOC model, which employs an empirical correction model involving a fitted parameter set.
Overall, the results of this work demonstrate that correcting for spin-projection can improve the performance of DFTcalculations of transition metal complex J-couplings and spin crossover gaps. However, for some systems this correction is modest. Also, as expected, the quality of calculations on systems with multiple contaminating spin-states are degraded by the use of the AP model. In such circumstances, more complicated projection models would likely improve calculated results. Nevertheless, this investigation supports the use of the AP model in calculations on molecular systems of the sort considered in this work.

Acknowledgements
The  S-1    1 Spin quantum number of each spin site 2 The first letter denotes geometry, C stands for crystal, L stands for low spin, H stands for high spin, A stands for AP; the last two letters denote method, SP stands for spin projected, NP stands for non projected, AP stands for approximate projection 3 Spin square expectation value at AP geometry for both HS and LS state S-2  1 Spin quantum number of each spin site 2 The first letter denotes geometry, C stands for crystal, L stands for low spin, H stands for high spin, A stands for AP; the last two letters denote method, SP stands for spin projected, NP stands for non projected, AP stands for approximate projection 3 Spin square expectation value at AP geometry for both HS and LS state S-3