The Role of Anisotropic Exchange in Single Molecule Magnets: A CASSCF/NEVPT2 Study of the Fe₄ SMM Building Block [Fe₂(OCH₃)₂(dbm)₄] Dimer

Abstract

The rationalisation of single molecule magnets’ (SMMs) magnetic properties by quantum mechanical approaches represents a major task in the field of the Molecular Magnetism. The fundamental interpretative key of molecular magnetism is the phenomenological Spin Hamiltonian and the understanding of the role of its different terms by electronic structure calculations is expected to steer the rational design of new and more performing SMMs. This paper deals with the ab initio calculation of isotropic and anisotropic exchange contributions in the Fe(III) dimer [Fe${}_2$(OCH${}_3$)${}_2$(dbm)${}_4$]. This system represents the building block of one of the most studied Single Molecule Magnets ([Fe${}_4$RC(CH${}_2$O)${}_3$)${}_2$(dpm)${}_6$] where R can be an aliphatic chain or a phenyl group just to name the most common functionalization groups) and its relatively reduced size allows the use of a high computational level of theory. Calculations were performed using CASSCF and NEVPT2 approaches on the X-ray geometry as assessment of the computational protocol, which has then be used to evinced the importance of the outer coordination shell nature through organic ligand modelization. Magneto-structural correlations as function of internal degrees of freedom for isotropic and anisotropic exchange contributions are also presented, outlining, for the first time, the extremely rapidly changing nature of the anisotropic exchange coupling.

1. Introduction

The Spin Hamiltonian (SH) provided the key to access to the microscopic explanation of the magnetic properties of materials starting from the Heisenberg, Dirac, and Van Vleck Hamiltonian proposed by Heisenberg in the 1928 [1,2,3]. Although it is widespread in both chemists and physicists communities, the direct mapping of the different SH terms with quantum mechanical approaches is still an relevant research field and its integration into a solid theoretical framework has been achieved only recently [4,5,6,7,8,9,10,11]. The combination of the Spin Hamiltonian theory together with computational methods allowed to access to microscopic insights from electronic structure theories that can be used as guide both for the interpretation of phenomena and the design of new materials. The possibility to calculate Spin Hamiltonian parameters from first principles is also of particular importance for all those systems where the origin of a complex magnetic structure cannot be easily disentangled. In fact, the inclusion of several Spin Hamiltonian terms can rapidly lead to an over-parametrization problem in the fitting of the experimental results. A remarkable example of a class of systems, showing very rich flavours of complex magnetic properties, is represented by single molecule magnets (SMMs). Many hopes are set on SMMs as possible active elements in several technological areas as spintronic and quantum computing [12,13,14,15], just to name a few of them. Although many efforts have been devoted to the study of various Spin Hamiltonian terms that contribute to SMMs’ magnetism, their usual multi paramagnetic ion-nature made these attempts not fully successful. Among the Spin Hamiltonian parameters, the most elusive one is by far the anisotropic exchange interaction. Although this interaction is generally neglected in the interpretation of experiments, this might not always be a correct assumption and especially for those systems where single-ion anisotropies do not necessarily dominate the low-lying part of the spin spectrum, this interaction might play a fundamental role, as also recently reported [16]. As showed by the authors in a previous paper [17], a possible route to the computation of the anisotropic exchange has been presented but unfortunately, a conclusive claim about the goodness of its applicability has not been possible. The reason mainly lies in the lack of benchmark references, both of experimental and computational origin. In order to get some more insights on the importance of this interaction in SMMs the most valuable strategy would concern the application of a higher level of theory with respect to the usually employed DFT, but unfortunately this is not possible for systems as large as those of practical interest and the study must be restrained to its most characteristic building block fragments.

In this respect, the high spin Fe(III)-oxo dimer Fe${}_2$(OCH${}_3$)${}_2$(dbm)${}_4$ with Hdbm = dibenzoylmethane [18] (Fe${}_2$ from now on) represents a good benchmark system mimicking the principal exchange pathways occurring inside propeller shaped SMMs. The structure of Fe${}_2$ is reported in Figure 1. This class of SMMs represents one of the most investigated systems and represents a valuable ensemble of good candidates for future device fabrication. An experimental EPR study of Fe${}_2$ shows that this system behaves as two $S=5/2$ iron ions whose spectrum is well described by the Spin Hamiltonian [19].

$$H_s(s_1,s_2)=J_{12}({\overrightarrow{\mathbf{s}}}_1·{\overrightarrow{\mathbf{s}}}_2)+{\overrightarrow{\mathbf{s}}}_1·{\mathbf{D}}_{12}·{\overrightarrow{\mathbf{s}}}_2+{\overrightarrow{\mathbf{s}}}_1·{\mathbf{D}}_1·{\overrightarrow{\mathbf{s}}}_1+{\overrightarrow{\mathbf{s}}}_2·{\mathbf{D}}_2·{\overrightarrow{\mathbf{s}}}_2$$

(1)

The leading interaction is the isotropic interaction $J_{12}$ whose value is 15.4 cm${}^{-1}$. For symmetry reasons single ion anisotropic terms are equal and of the easy plane kind. The anisotropic exchange has been estimated to be of the same order of magnitude of the single ion term suggesting its importance in the final contribution to the magnetic behaviour of this paramagnetic complex. We here present the calculation of the anisotropic exchange interaction for Fe${}_2$ at post Hartree-Fock (post-HF) level of theory. An initial assessment of the theoretical scheme will be followed by an analysis of the effect of the organic ligand nature and finally the investigation of the magneto-structural correlation pathways inside the dimer will be shown. This computational study points the isotropic exchange as a robust and smoothly varying property with respect to both the second coordination shell nature and the internal geometrical parameters. In contrast, the anisotropic exchange has been found to be rapidly modulated by these structural perturbations and the possible consequences of this effect on the magnetisation relaxation properties have been discussed.

2. Results and Discussion

2.1. Method Assessment

2.1.1. X-ray Structure

The calculation of isotropic $J_{12}$ and anisotropic exchange $D_{12}$ coupling parameters of a dimer of Fe(III) ions by ab initio methods can still be considered a challenging task. Indeed, the energetics of these parameters goes from a few to a fraction of wavenumbers, demanding a very high computational accuracy. In this regard, to verify the reliability of the proposed computational protocol, we have decided not to include further sources of possible errors given by the optimisation procedure or some modelization of the geometry. For such a reason the Fe${}_2$ X-ray structure (1-Ph${}_{\mathbf{xray}}$ from now on) has been used for the first benchmark calculation at the post-HF level. Both CASSCF and NEVPT2 approaches (see Theory and Computational Methods and Computational Details sections) have been applied and results are listed in

Table 1. Calculated single ion and anisotropic exchange anisotropy tensors.

	J${}_{12}$	D${}_{\mathit{SI}}$	E${}_{\mathit{SI}}$/D${}_{\mathit{SI}}$	D${}_{\mathit{EX}}$	E${}_{\mathit{EX}}$/D${}_{\mathit{EX}}$
1-Ph${}_{\mathbf{xray}}$	5.32	0.329	0.113	−0.108	0.321
1-Me${}_{\mathbf{xray}}$	5.52	0.090	0.124	0.411	0.135
1-Me${}_{\mathbf{opt}}$	2.92	0.197	0.105	0.116	0.131
Exp. ${}^a$	15.4	0.749	0.097	−0.159	0.176

All reported values are in cm${}^{-1}$. ${}^a$ Experimental values are taken from Ref. [19].

. As expected [6,10,23,25], the inclusion of correlation only in the five 3d orbitals of the two Fe(III) leads to a very poor agreement with the experimental value. Indeed, a ferromagnetic $J_{12}$ value of −9.12 cm${}^{-1}$ was found (see Equation 5) versus an experimental antiferromagnetic value of 15.4 cm${}^{-1}$. The introduction of dynamical correlation via NEVPT2 recovered some of the important contributions given by the bridging ligand orbitals to the super-exchange mechanism and both the experimental sign and an acceptable magnitude of 5.32 cm${}^{-1}$ for $J_{12}$ have been then reached. Possible strategies to further improve this result would require the introduction of a dedicated CI procedure [33] or equivalent [34] on top of the CASSCF instead of NEVPT2. Another possible source of the discrepancy can reside in the limited flexibility of the wavefunction resulting by a CAS(5,5), i.e., the so called 3d double-shell effect [35] or the missing of the explicit inclusion of bridging ligands orbitals in the CAS [35]. Even in these cases, the dimension of the CAS should be at least doubled. However, the computational demand for such an approach would require a severe reduction in the size of the system and will not be addressed here. Assessed that the inclusion of the dynamical correlation via perturbative approach is mandatory for a better description of the electronic states energy ladder we used the same protocol for the calculation of the anisotropic Spin Hamiltonian terms. In our treatment, only the SOC contribution is included while it is worth to remember that experimental single ion values contain both spin-spin and SOC contributions. Nevertheless, on the basis of the results reported in ref. [17,36], the lacking of the spin-spin contribution could not account for the big difference between the computed and the experimental single ion anisotropy values and the simulation protocol employed leads to a slight underestimation of the single ions contribution. However, the computed anisotropic exchange parameters result in fairly good agreement with the experimental ones both in sign and magnitude [19]. In Table 1 we reported the computed values for anisotropic exchange anisotropy tensors. In conclusion, the computational set up can be considered to be reliable enough to be used for the magneto-structural analysis.

2.1.2. Simplified X-ray Structure

A magneto-structural correlation study requires a single point calculation for several different geometries corresponding to a specific path on the molecule potential energy surface expected to strongly influence the magnetic properties of the system under study. In our case, the single point calculations are made at the CASSCF/NEVPT2 level and, therefore, very demanding both from hardware and wall-time point of view. To reduce both of them, we have modelled the X-ray structure by replacing the phenyl groups of the dbm ligands by methyls groups (1-Me${}_{\mathbf{xray}}$ model). Such a choice was supposed to slightly affect the ligand field felt by the two iron ions. Indeed, the modelization consists only in the substitution of an sp${}^2$ carbon (phenyl) with an sp${}^3$ one (methyl) as functional group in alpha position to the CO groups of the dbm ligands, which is the actual moiety that binds the metals. In fact, the computed $J_{12}$ was practically unaffected by this substitution ($J_{12}$ = 5.52 cm${}^{-1}$). On the other hand, $D_{12}$ was found, unexpectedly, to change both in sign and magnitude. The change of the outer coordination shell showed to be enough to turn the easy-axis kind of the exchange anisotropy to an easy-plane one with $D_{12}$ = 0.411 cm${}^{-1}$ and E/D = 0.135 (see Table 1). In order to confirm the reliability of this result, we checked its convergence with respect to the number of excited spin multiplets introduced in the SOC diagonalization adding up to 99 quintuplets states. Although the $D_{12}$ value gets reduced to 0.344 cm${}^{-1}$, its easy-plane anisotropy was maintained. Therefore, the variation from easy-axis to easy-plane anisotropy can be ascribed to the variation of the π contribution induced by the substitution of the CO groups. Such a result is very interesting since it shed some light on how subtle the effects on the anisotropic exchange terms can be and how important the chemical tailoring in designing magnetic materials should be.

2.2. Magneto-Structural Correlations

As reported several times in the literature [17,18,37], the most interesting collective variable connected to the isotropic exchange coupling is the FeOFe angle and in this work we further extend the study of this correlation including the assessment of the anisotropic part of the exchange. In the attempt to address magneto-structural correlations by mean of solely first principle calculations, we did not use any experimental reference structure, differently from what previously done [17]. To avoid the introduction of any bias in the choice of the collective variable along which we extract the magneto-structural correlations, we optimised the 1-Me${}_{\mathbf{xray}}$ geometry in vacuum (1-Me${}_{\mathbf{opt}}$) followed by the calculation of its hessian. Both calculations have been done with DFT (see Computational Details). The root mean square displacement (RMSD) between 1-Me${}_{\mathbf{opt}}$ and 1-Me${}_{\mathbf{xray}}$ is only 0.143 Å. The optimised main geometrical parameters are reported in

Table 2. X-ray and optimized structural parameters.

	Fe–Fe (Å)	Fe–O (Å)	FeOFe (Deg.)
1-Ph${}_{\mathbf{xray}}$	3.15	1.99	100.69
1-Me${}_{\mathbf{opt}}$	3.12	2.03	102.69

with their correspondent X-ray geometrical values.

At this point, from the diagonalization of the hessian, it has been possible to obtain the normal modes of the molecule, which were used as a guide in the choice of the most natural collective variable containing the FeOFe internal displacement.

In Figure 2 the weights of the two possible FeOFe angles inside all the normal modes is reported. Due to the presence of the point inversion symmetry the two angles are linearly combined, as the absolute value of the FeOFe amplitudes in Figure 2 shows. The normal mode at frequency ∼250 cm${}^{-1}$ (q from now on) is the one showing the larger symmetric FeOFe displacement among the whole vibrational spectra and it was therefore chosen as the collective variable to be used for the magneto-structural correlation study. Figure 3 shows a schematic representation of the cartesian displacements corresponding to the normal mode of vibration q.

At this stage, it is fundamental to stress that, due to the critical dependence of magnetic interactions on nuclear displacements, we need to consider the next results as standing alone with respect to the experimentally measured interactions. Indeed, the most natural reference state will become the 1-Me${}_{\mathbf{opt}}$ model in place of the 1-Ph${}_{\mathbf{xray}}$ model and the magneto-structural correlations will be referred to this geometry. For the study of the magneto-structural correlations, we produced six distorted structures corresponding at ±1q, ±2q and ±3q displacements, expressed in terms of normal mode unit, with respect to the equilibrium geometry. The ensemble of angles spanned by these structures is rather large (94.9${}^{\circ }$÷ 106.4${}^{\circ }$) and fully includes all the experimental scenarios encountered so far in both propeller-shaped SMMs and this oxo-bridged Fe2 dimer family [18,37].

As expected, the $J_{12}$ value nicely correlates with the FeOFe angle and the reduced value calculated for the optimised geometry is fully accountable to the reduction of the FeOFe during the optimisation. Indeed, according to the plot reported in Figure 4, for the experimental FeOFe value of 102.7${}^{\circ }$ the expected $J_{12}$ would be ∼5.5 cm${}^{-1}$, in good agreement with the calculation done on the 1-Ph${}_{\mathbf{xray}}$ model.

In contrast, $D_{12}$ values do not change smoothly with the angle FeOFe and span a relatively large amplitude of values for the selected configurations (Figure 5). Interestingly, the trend presents a sharp minimum correspondent to a negative $D_{12}$ value of −0.265 cm${}^{-1}$ for a FeOFe angle of 99${}^{\circ }$. Such a result confirms that a tailored chemical design is mandatory to obtain systems with ad hoc magnetic properties and, more interestingly, that a transition from an easy-plane to easy-axis or vice versa could be possible and instantaneously achieved by perturbing selectively the system with an energy quantum correspondent to the normal mode at issue. This last point might become important when thermal effects are considered. A long lasting magnetisation vector at high temperature is the main goal to be achieved in the SMMs research field and in order to do that the control of the spin-phonon coupling interaction is necessary. Indeed, at finite temperature, the modulation of the Spin Hamiltonian parameters is the trigger for the magnetisation vector relaxation and the larger this modulation is, the faster will be the relaxation. In this respect, the presence of degrees of freedom able to vary significantly a magnetic interaction inside the metal core of the molecule is rather unwelcome. The results presented in the work highlight a possible relevant role of the anisotropic interaction in the SMMs relaxation processes. Indeed, although its relatively small intensity, this property happens to be rapidly modulated by the FeOFe degree of freedom, which has been demonstrated to vary of several degrees due to thermal fluctuations [38,39], and thus opening a possible efficient relaxation pathway.

3. Theory and Computational Methods

The calculation of the anisotropic exchange coupling tensor has already been described by the authors [17] and here will be only briefly reviewed. This procedure is based on the twofold possibility to write the Spin Hamiltonian of two coupled ions: in a multi spin fashion, as in Equation (1), and in a giant spin fashion

$${\mathbf{H}}_{\mathbf{S}}=\overrightarrow{\mathbf{S}}·{\mathbf{D}}_{\mathbf{S}}·\overrightarrow{\mathbf{S}}$$

(2)

where $H_S$, function of the total spin $\overrightarrow{\mathbf{S}}={\overrightarrow{\mathbf{s}}}_1+{\overrightarrow{\mathbf{s}}}_2$, acts as a molecular Hamiltonian resuming all the single ions information. The parameters of the two Hamiltonians (Equations (1) and (2)) are related together by some projection coefficients (see Equation (3)), which could be numerically evaluated from Clebsh-Gordan coefficients, or by direct diagonalization of the Heisenberg Hamiltonian once Js values are known [20].

$${\mathbf{D}}_{\mathbf{S}}=d_1^S{\mathbf{D}}_{\mathbf{1}}+d_2^S{\mathbf{D}}_{\mathbf{2}}+d_{12}^S{\mathbf{D}}_{\mathbf{12}}$$

(3)

Equation (3) represents the basis for the determination of $D_{12}$. Indeed, the single ion tensors of a couple can easily be obtained with the diamagnetic substitution method (i.e., substituting a Fe(III) ion in the dimer one at a time with a diamagnetic structurally equivalent Ga(III) ion) and $D_S$ can be assessed as the anisotropy $D_{HS}$ of the high spin state, i.e., $S=s_1+s_2=10$. In this framework

$${\mathbf{D}}_{\mathbf{12}}=\frac{{\mathbf{D}}_{\mathbf{HS}}-d_1^{HS}{\mathbf{D}}_{\mathbf{1}}-d_2^{HS}{\mathbf{D}}_{\mathbf{2}}}{d_{12}^{HS}}$$

(4)

where $d_1^{HS}$, $d_2^{HS}$, $d_{12}^{HS}$ are the projection coefficients that relate the various tensors to the global high spin state. This procedure can be applied with both DFT or multi-reference electronic structure methods, as long as the spin contamination of the high spin state is negligible. The used post-HF approach consisted in performing a complete active space self-consistent field (CASSCF) followed by N-electron valence perturbation theory (NEVPT2) level of theory [21]. For the CASSCF calculation of the electronic structure of the Ga(III) doped system, a (5,5) active space, i.e., including all d-like orbitals and the five unpaired electrons, has been used. The sextuplet, all triplets and all doublets CI solutions have been calculated in a state-average fashion and used for the subsequent NEVPT2 step. While this strategy for extracting transition metal single ion anisotropy has already been widely validated [22], calculations on multi-ion transition metals are limited to a few papers [23,24,25,26,27,28]. So far, only active spaces extended over the d-like MOs, i.e., (10,10) in this case, have been considered in this kind of calculations and the same approach will be followed here in order to limit the already large computational demand. Nevertheless, as already pointed out by other papers concerning the isotropic part of the exchange, the inclusion of the MOs participating to the super-exchange interaction might be fundamental [6,10]. Such extension of the active space would be unaffordable and the effect of the ligand MOs correlation is included only at the perturbative level through NEVPT2 correction. For the originally proposed DFT recipe the choice to project the $\mathbf{D}$ tensor on the high spin multiplet solution was forced by the impossibility of DFT to reproduce intermediate spin states. On the other hand, such a limitation is not present in post-HF methods for which the access to each spin multiplicity is granted. Although the choice of the spin multiplet where to project the $\mathbf{D}$ tensor is arbitrary, the high spin solution remains the most practical pathway as it is the one with the lower number of excitations. Moreover, this solution is also in principle less affected by S-mixing problems and therefore much more prone to be included in a Spin Hamiltonian mapping. In light of this discussion, the CASSCF/NEVPT2 calculations on the non-doped dimer have been restricted to the state-average optimisation of only all the possible two higher spin multiplicity excitations i.e., one endeca-plet and ninety-nine nonu-plets. Finally, the CI solutions obtained from these calculations have been used as basis to diagonalize the spin-orbit coupling operator and, in the context of quasi-degenerate perturbation theory (QDPT), the $2S+1$ solutions of interest (the $2S+1=6$ and $2S+1=11$ solutions for the doped and non-doped systems, respectively) have been used as model spaces to fit the Spin Hamiltonians [9]. The CASSCF/NEVPT2 calculation of the electronic structure makes also possible to extract the isotropic exchange coupling constant from the energy splitting between the $S=10$ and $S=9$ multiplet at the NEVPT2 level according to the formula

$$J_{12}=(E\left(S\right)-E(S-1))/S$$

(5)

Among the many $S=9$ solutions the one that should be considered for this purpose is the one corresponding to two misaligned local $s=5/2$ spins, i.e., the solution with the largest weight on the determinant with all the active orbitals singly filled.

4. Computational Details

All the calculations have been done with ORCA [29] employing a def2-TZVP [30] basis set for magnetic elements and their first neighbours (oxygen atoms), while def2-SVP basis set has been used for all the other atoms. The RI-J approximation along with the def2-TZVP/J auxiliary basis set [31] for all the elements has been used. Speaking in the ORCA notation, grids were set to 5 and VeryTightSCF convergence criteria have been used. This set up has been tested in the context of CASSCF method with respect to calculations done with the def2-TZVP basis set on every atomic kind and no significant differences have been noted on the energy ladder of the CASSCF excited states. The computation of the magnetic properties included only SOC as relativistic interactions, neglecting for instance spin-spin coupling. Optimisation and normal modes calculations have been done at the DFT level of theory employing the PBE functional [32].

5. Conclusions

We have presented the ab initio calculations of the isotropic and anisotropic exchange contributions in the dimer of formula [Fe${}_2$(OCH${}_3$)${}_2$(dbm)${}_4$] which represents the basic building block for the Fe${}_4$ SMM family. The study has been performed on both the X-ray and a modelled structure, where in the latter the phenyl groups were substituted by methyl groups. As already showed in literature, it turned out that the addition of dynamical correlation in the protocol is mandatorily needed for a nice agreement between the experimental and computed parameters. From this study turned out that even small changes in the outer coordination sphere can strongly influence the anisotropic exchange values, while leaving the isotropic coupling almost unaffected. A magneto-structural correlation study has also been presented for both $J_{12}$ and $D_{12}$ as a function of the bending of the FeOFe angle on the basis of the hessian analysis. A smooth behaviour, in line with the already reported trends in literature, was found for $J_{12}$, while a more complex behaviour was observed for $D_{12}$, where positive and negative values were spanned. Such a result, in addition to the phenyl-methyl substitution effect, is important as it opens scenarios where the easy axis or easy plane behaviours can be a priori chemically designed. Possible implications on the spin-phonon relaxation mechanisms have also been addressed.