Calculation of the relative energies in a Spin Crossover material

Question

Spin Crossover (SCO) complexes are a particular type of molecular entities which are usually formed by a metal ion (general $\ce{Fe(II)}$ or $\ce{Co(III)}$) complexed by several ligands with donor N atoms. The d orbitals splitting ($\Delta$ or 10Dq) produced in certain octahedral coordination geometries allows two spin configurations which are metastable and accesible with a stimuli such as temperature, light, pressure and electric field among others, producing a well-known hysteretic behavior.

Unfortunately, the relative energies of HS and LS in a SCO complex are extremely difficult to tackle in a DFT molecular calculation. Although some elaborated approaches have been considered, such as considering the effect of the dipolar interactions or using intrincated exchange-interaction funcionals, no evident recipe has been proposed to correctly simulate such materials.

In example, in the previous work by Cirera et al., they presented a selection of SCO complexes and calculated an estimation of the transition temperature by using Gibbs free energy expression:

$$\tag{1}\Delta G = \Delta H - T\Delta S $$

Since at equilibrium, $\Delta G$ vanishes, the transition temperature ($T_{1/2}$) is determined by:

$$\tag{2}T_{1/2} = {{\Delta H} \over {\Delta S}}$$

From this equation (more details on the article), the key contribution to obtain the transition temperature is the enthalpy difference between the high and low spin state, which can be computed using DFT.

To illustrate the complexity of this procedure, I have chosen the molecule $\ce{Fe(SCN)_2(Phen)_2}$, considered the Droshophila of SCO materials, which has been modeled with TPSSh Funcional. The results obtained are:

Experimental $T_{1/2}$ = 176.5 K

Theoretical (TPSSh Functional) $T_{1/2}$ = 454 K (basis set 1) and 237 K (basis set 2)

As can be easily seen, we are still far from the accurate result.

My question is: Up to now, which is the most correct way of simulating the SCO transition? How realistic is it to reproduce the hysteresis on the magnetization properties in a temperature sweep?

Answer 1

The question highlights the difficulty in calculating an ab initio value for $T_{1/2}$ for a molecule like $\textrm{Fe(Phen)}_2\textrm{(SCN)}_2$, and specifically the fact that a highly cited 2018 paper predicts values that are between 60K and 278K higher than the experimental value of 176.5K.

The paper highlights a major danger of using DFT, and in using theory to try to predict something that is out of reach with today's computing technology and algorithms.

The paper uses many approximations to calculate $T_{1/2}$ for a molecule. They begin with the approximation that it can be calculated as follows:

\begin{align}T_{1/2} &= \frac{\Delta H}{\Delta S} \\ &\approx \frac{\Delta E_\textrm{electronic}+\Delta E_\textrm{vibrational}}{\Delta S}, \end{align}

where only $\Delta E_\textrm{electronic}$ is attempted accurately because they claim near their Eq. 3 that "the harmonic approximation" is good enough for $\Delta E_\textrm{vibrational}$ and $\Delta S$ (see their Eq. 8 for more details). In general if you want to calculate $\Delta S$ and $\Delta H$ ab initio, resources that might help you are this paper for $\Delta S$, and this step-by-step guide for $\Delta S$, but as the authors claim that the harmonic approximation is good enough for everything except $\Delta E_\textrm{electronic}$, let's zoom into how $\Delta E_\textrm{electronic}$ is calculated.

The system mentioned in the question is labelled S11 in the table below, and you can see that with 8 different functionals, values for the HS vs LS energy gap vary from -14.6 to +9.6 kcal/mol (a 25 kcal/mol range for a number estimated to be at most 15 kcal/mol in magnitude), with only 3 out of 8 functionals even giving the correct sign. For reference, the term "chemical accuracy" means an accuracy of +/- 1 kcal/mol. As someone who mainly works in the "high-precision ab initio" field, I would stop right here and start working on a different project, however the brave scientists noticed that the TPSSh functional predicted the sign correctly for all 20 systems:

I lack experience with DFT, so I personally would try to develop more experience to decide how significant it is to me that one functional is getting all the signs correctly (I think I personally would try more functionals until I find another one that gets all signs correctly, and then compare the magnitudes, since what i see in the table indicates that the magnitudes vary wildly from what people call "chemical precision"). However the authors have more experience with DFT than me, and they have decided to use the TPSSh functional for the rest of their calculations. The $\Delta E_\textrm{electronic}$ could have been 1/3 of what it was for TPSSh if they used B3LYP* instead, and could have been the opposite sign and almost 2x larger in magnitude if they used M06 (which I understand is reputed to be good for transition-metal containing systems like the system in question), so you can see that we can get a very wide range of $T_{1/2}$ values if we wish.

After deciding that the TPSSh functional was the only one that got the signs correctly for all 20 systems studied in the above table, they increased the size of the basis set from what they call BS1 to what they call BS2.

Your question notes that there is a 217K difference between the result obtained with BS2 vs BS1, a different 2 times larger than the amount needed to boil water.

Furthermore, the difference between the two basis sets is really not very big:

• BS1 = QZVP for the Fe atom and TZV for all other atoms
• BS2 = TZVP for all atoms (this may have been a typo, since it would be strange that they choose not to use QZVP for the Fe atom in the bigger basis set, but maybe that would make the overall basis set too big, but I find that hard to believe and it doesn't matter much because it's just one of the dozens of overall atoms).

This tiny difference in the basis set (keeping $\zeta$ the same but just increasing polarization functions, and perhaps also reducing $\zeta$ for the Fe atom) caused a difference of 217K in $T_{1/2}$, so again I would probably stop here and try to do a basis set convergence study (for example, looking at cc-pVDZ, cc-pVTZ, and an extrapolation, and comparing it with another 2Z/3Z extrapolation). Personally I probably would have tried to study 5 molecules more thoroughly (with a basis set convergence study and more supporting evidence for my choice of functional) rather than 20 molecules the way they did, but I also do not work in this field and cannot fully appreciate the motivations behind the study.

---

So how can the prediction for $T_{1/2}$ be calculated more accurately? I would recommend the following things:

• This study showed that a minor difference in basis set resulted in a 217K change in $T_{1/2}$, so I would do a basis set extrapolation with at least two different basis set series and see if there is some reasonable number on which the different series agree. If the basis set extrapolations give numbers that are vastly different, and I cannot afford to use a larger basis set or to incorporate explicit correlation through something like an F12 method, then there is no shame in saying that the problem has defeated me. With so many atoms, it is understood that the complete basis set is not approachable easily (if at all).

• This study shows that different functionals give vastly different results. The number of atoms is rather large, but not out of reach for CCSD and CCSD(T). Transition metal systems can be challenging, especially for multi-reference spin states, but this system does not look impossible for coupled cluster. A comparison between MP2, CCSD, and CCSD(T) may show energies that are closer together than the different functionals used in this paper, and the agreement between CCSD and CCSD(T) can be an inidication of the overall accuracy of the calculation. Another option would be CASSCF, RASSCF or GASSCF and CASPT2/RASPT2 with OpenMOLCAS, which may allow for better accuracy but will cost more human time since these are not black box methods.

• The DKH relativistic treatment in the paper is most likely okay, but for the same cost they could have done X2C which is equivalent to "infinite order DKH".