Single Molecule Magnet Properties of Transition Metal Ions encapsulated in Lacunary Polyoxometalates: a Theoretical Study

Single molecule magnet (SMM) properties of transition metal complexes coordinated to lacunary polyoxometalates (POM) are studied by means of state of the art ab initio methodology. Three [M( -SiW10O36)2] (M = MnIII, FeIII, CoII) complexes synthesized by Sato et al. (Chem. Commun. 2015, 51, 4081–4084) are analyzed in detail. SMM properties for the CoII and MnIII systems can be rationalized due to the presence of low energy excitations in the case of CoII, which are much higher in energy in the case of MnIII. The magnetic behavior of both cases is consistent with simple d-orbital splitting considerations. The case of the FeIII complex is special, as it presents a sizable demagnetization barrier for a high spin d5 configuration, which should be magnetically isotropic. We conclude that a plausible explanation for this behavior is related with the presence of low lying quartet and doublet states from the iron(III) center. This scenario is supported by ab initio Ligand Field analysis based on CASSCF results, that pictures a d-orbital splitting that resembles more a square-planar geometry than an octahedron, stabilizing lower multiplicity states. This coordination environment is sustained by the rigidity of the POM ligand, that imposes a longer axial bond distance to the inner oxygen atom in comparison to the more external, equatorial donor atoms.


Introduction
Since the discovery of slow relaxation of the magnetization properties in mononuclear Tb(III)phthalocyanine compounds by Ishikawa et al., 1 the quest of new Single Molecule Magnets (SMMs) with larger demagnetization barriers shifted from the synthesis of large polymetallic complexes with the highest possible multiplicity values toward lower nuclearity compounds. The rationale for the design of large, polymetallic compounds to achieve better SMMs was the maximization of the total spin of the system to achieve higher values for the demagnetization barrier. The relaxation barrier (U) is directly proportional to the squared spin value for the ground state (U = DS 2 ), where S is the total spin of the ground state and D is the zero field splitting parameter, which quantifies the energy splitting of the lowest energy multiplet, associated with magnetic anisotropy. In this way, larger spins would be more likely to yield higher values of U.
Since 1995, when an Fe19 complex with S=33/2 was presented, 2 5 The current record value was presented in 2015 and corresponds to a Fe42 compound with S=90/2, in which 18 Fe III centers are ferromagnetically coupled and the remaining 24 iron cations are low spin iron(II). 6 Despite the success in the synthesis of molecules of increasingly higher spin, these systems did not present larger demagnetization barriers than existing, lower nuclearity compounds. [7][8][9] Furthermore, larger spin systems did not present any SMM properties in most cases. This apparent departure from the U = DS 2 relation can be understood noticing that the value of the zero field splitting parameter (D) is itself proportional to 1/S 2 , cancelling the explicit S 2 dependence of U. [10][11][12] Thus, chemical tuning of the coordination environment for the maximization of D is unavoidable to design improved SMMs, independently of their number of magnetic centers and total spin. Under this perspective, it is certainly simpler to conceive a tailored coordination environment for a single metal center in a small coordination compound than in a polynuclear compound, where the magnetic anisotropy of several centers has to be aligned and maximized simultaneously.
In this way, great attention is devoted to design new SMMs candidates using monometallic coordination compounds, where the environment for the metal center is carefully chosen to maximize magnetic anisotropy. [13][14][15][16] The employment of rigid or bulky ligands to force a particular coordination environment by geometrical constraints or steric hindrance has proven to be a successful strategy to design new SMMs with improved properties. Apart for new lanthanide phthalocyanine complexes synthesized after [TbPc2] -, 17,18 we can mention several examples of remarkable SMMs based on polydentate, rigid or bulky ligands, as Dy III sandwich complexes, 19 dicoordinate Fe I and Ni I compounds stabilized by bulky ligands 20,21 or endohedral lanthanide systems. 22 In this sense, lacunary polyoxometalates appear as promising ligands for SMM systems, given they large volume and extreme rigidity, that permits to design highly structured coordination environments for the promotion of magnetic anisotropy. Several examples of SMMs based on lanthanide magnetic centers coordinated to diamagnetic polyoxometalates exist in the literature, exhibiting large demagnetization barriers. 23,24 The influence of the polyoxometalate ligands in the energy splitting and magnetic anisotropy of the lanthanide ion has been studied in detail and has been rationalized and linked to Crystal Field Theory considerations. 25,26 In contrast, examples of polyxometalate-coordinated transition metal systems are comparatively scarce, and the relation between the observed SMMs properties and their electronic structure has not been yet described in detail.
The goal of this article is to analyze the SMMs properties of transition metal based complexes with polyoxometalate ligands by means of state of the art ab initio calculations. Sato et al. 27 reported magnetic studies about three transition metal complexes coordinated to two lacunary [ -SiW10O36] 8silicontungstates, acting as tridentate ligands. Two of the studied compounds are fieldinduced SMMs, with demagnetization barriers of 6.3 cm -1 (Fe III ) and 13.4 cm -1 (Co II ) under a static field of 0.1 T. The third compound is based on Mn III and shows temperature and frequency dependence for the ac magnetic susceptibility at very low temperatures, not allowing for the estimation of a demagnetization barrier. The coordination geometry around the transition metal centers is similar for the three complexes and can be described as distorted octahedral.
Interestingly, the Fe based complex is assigned to be high spin Fe III , with a formally isotropic d 5 configuration, being the only reported example of a high-spin d 5 system presenting SMM behavior.
The article is divided in two parts: (i) the three studied systems are analyzed in terms of ab initio electronic structure calculations to describe the origin of their observed magnetic properties in terms of excitation energies, zero-field splitting, g-factors and ab initio Ligand Field analysis. (ii) The case of high spin Fe III presenting slow relaxation of the magnetization is studied in detail by ab initio methodologies. Computational results are analyzed to describe plausible and unlikely mechanisms for the unusual SMM properties reported for this compound.

Computational Details
All electronic structure calculations were performed using the ORCA 3.0.3 package. 28,29 Multireference ab initio calculations based on the Complete Active Space Self-Consisted field (CASSCF) 30,31 method considered the five orbitals of the 3d shell and dynamic correlation was included through the N-electron Valence State Perturbation Theory (NEVPT2). 32,33 Scalar relativistic effects were accounted by the Douglas-Kroll-Hess Hamiltonian at 2 nd order (DKH2), [34][35][36][37] while state mixing due to spin orbit coupling was represented by a quasi-degenerate perturbation theory approach (QDPT). The effect of dynamic correlation was included in the QDPT step as a diagonal correction on the non-relativistic state energies. Molecular geometries for all complexes were directly obtained from X-Ray crystallographic structures. 27 As ab initio calculations are computationally demanding, truncated models of the three compounds were constructed keeping the immediate coordination environment of the transition metal ion and its second and third neighbors. Fourth neighbors were replaced with hydrogen atoms that were optimized to adjust their bond length, freezing the corresponding bond angles. In this way, the immediate coordination environment of the transition metal and its second neighbor W and Si atoms is preserved without any changes (See Figure 2). Cartesian coordinates for the molecular models are presented as Supporting Information. All calculations were performed using the DKH optimized SARC basis set for W, 38 and relativistically recontracted versions of the corresponding Ahlrichs basis set for the remaining atoms. 39,40 The COSMO solvation model with water as solvent was included in all calculations. 41 A comparison to gas-phase calculations for Fe III -POM is presented in Table S2. dorbital energies, spin orbit coupling and ligand field Hamiltonian parameters were directly fitted from CASSCF matrix elements by means of the ab initio Ligand Field (AILF) methodology. [42][43][44] In a nutshell, the AILF approach consists in the direct mapping of each matrix element of the configuration interaction matrix from a CASSCF calculation to a model Hamiltonian matrix.
Dynamical correlation can be included by replacing the CASSCF excitation energies with corrected values (i.e. energies from a NEVPT2 run). In this case, the model Hamiltonian matrix was constructed in terms of usual Ligand Field Hamiltonian parameters (i.e. Racah parameters and the general form of the ligand field matrix for d orbitals considering 15 parameters, which can be further condensed in terms of common ligand field approaches, such as Angular Overlap Model parameters). 45,46 In this way, the AILF approach treats high-and low-symmetry systems in an equal footing and can fit a larger set of parameters than a direct fit to state energies, as in considers all matrix elements of the configuration interaction matrix and not only its eigenvalues. NEVPT2 energies can eventually replace CASSCF values in the reconstruction of the ab inito CI matrix, while wavefunctions still correspond to a CASSCF calculation. Continuous shape measurements (CShM) were performed using the SHAPE 2.0 code. 47,48 CShM allows for a quantitative evaluation of the similarity between an arbitrary shape and reference shapes, discarding any effect associated with rotations and scaling. In practice, the degree of matching between target and reference shapes is expressed in terms of a number, where 0 indicates a perfect match that deteriorates with increasing values. As the studied compounds were six-coordinated, we compared the directly bonded MO6 (M = Mn, Fe, Co) fragments with the octahedron, pentagonal pyramid and trigonal prism reference shapes.

Results and Discussion
According to a continuous shape measurement (CShM) analysis, the coordination environment of the three complexes can be described as close to octahedral, with S values of 1.88, 1. The orbital splitting of all complexes is expected to be close to the classical picture of octahedral geometry, with a sizable distortion due to the longer axial distances and angular deviation.
Following the rules proposed by Gómez-Coca et al. 49 for the prediction of SMM behavior in transition metal complexes (see Figure 1), we can expect that the Mn III -POM complex would hardly present a sizable relaxation barrier, given the large Jahn-Teller like distortion of their eg orbitals and the absence of components of the angular momentum operator connecting the dx2-y2 and dz2 orbitals, in line with its experimental behavior. The case of Fe III -POM is interesting, as their high spin d 5 configuration would not be compatible with SMM behavior, conflicting with experimental results, this example will be analyzed in detail later. Finally, Co II -POM is likely to present a relaxation barrier due to the excitations associated to the doubly and singly occupied t2g orbitals. These orbitals are less affected by their coordination environment than eg orbitals and are connected by matrix elements of the angular momentum operator. This orbital set is then suitable to develop a sizable demagnetization barrier, in agreement with their ac-susceptibility measurements.

ab initio calculations:
In order to describe the magnetic anisotropy of the ground state of the three studied complexes, we performed ab initio multiconfigurational calculations (CASSCF) for the model complexes, followed by a quasi-degenerate perturbation theory (QDPT) step to account for state interaction due to Spin-Orbit coupling (SOC). The CASSCF+QDPT method is widely employed in the description of SMM complexes 15,43,[50][51][52] because it allows for a multiconfigurational description of the ground and low energy excited states, which is a key condition to describe the state mixing due to SOC properly. In the QDPT step, nonrelativistic states from CASSCF are split in their corresponding Ms components that serve as the basis of the state interaction matrix associated with the SOC operator. In this way, wavefunctions after QDPT reflect the magnetic anisotropy of the system in their structure and allow for the evaluation of Spin-Hamiltonian parameters directly from their interaction with a magnetic field, without requiring any fit to a macroscopic magnetic response.
Considering the five 3d orbitals as the active space, a CASSCF(4,5) calculation was performed for Mn III -POM considering 5 quintets, 45 triplets and 49 singlets. The lower two d orbitals are dominated by dxz and dyz contributions, with a small energy separation of 289 cm -1 , the dxy orbital has an energy of 4288 cm -1 to complete the t2g set for an ideal octahedral geometry. The dz2 orbital is a bit higher in energy (6948 cm -1 ), while the dx2-y2 is markedly higher (19903 cm -1 ) (See Figure   2). Under this scenario, the high-spin d 4 configuration for the ground state will be (dyz) 1 (dxz) 1 (dxy) 1 (dz2) 1 and the first excitation energy is as large as 13449 cm -1 , due to the large destabilization of the dx2-y2 orbital that effectively quenches state mixing associated with SOC. The inclusion of NEVPT2 does not change this picture, yielding a slightly higher excitation energy (14237 cm -1 ). The calculated splitting of the ground multiplet is 13.1 cm -1 and 13.5 cm -1 for CASSCF and NEVPT2, respectively, with a slight separation of the ML = ±2 states due to some rhombic contribution to D (E/D = 0.03 for both methods). The combination of a relatively low barrier, together with the existence of some rhombic component in the ZFS and the non-Kramers' nature of the Mn III ion are in line with the observed lack of measurable relaxation barrier for Mn III -POM. Spin Hamiltonian parameters for all calculated models are presented in Table 1.
The calculated Zero Field Splitting (ZFS) parameter D is -3.27 cm -1 and -3.37 cm -1 for CASSCF and NEVPT2, respectively. These values are lower than the experimentally fitted D value (-5.28 cm -1 ), as reflected by the comparison of the experimental and calculated susceptibility curves (see Figure 3). In the low temperature range, the decrease of the T product is more pronounced in the experimental curve, in line with the higher experimental value of D in comparison to calculations.
The origin of this discrepancy relies in the importance of quintet-triplet excitations for a quantitative recovery of magnetic anisotropy. An analysis of the contributions of each excited state to the value of D reveals that quintet-quintet excitations contribute only ~ -1.2 cm -1 to D, which is -3.19 cm -1 in the second order perturbation theory approach based on a CASSCF(4,5) calculation.
It is remarkable that several quintet-triplet excitations contribute importantly to the zero field splitting value (see Table S1 for a complete list of contributions for all quintet, triplet and singlet states). It is important to stress that zero field splitting values originating from lower multiplicity states are expected to be underestimated in CASSCF calculations, are higher multiplicity states are over stabilized under this methodology, yielding too large energy differences with respect to lower multiplicity states. This issue will manifest more clearly in the case of Fe III -POM, where magnetic anisotropy is essentially related to sextet-quartet excitations. The case of Co II -POM is different, as it presents relatively low-lying first and second excited states at 266.7 cm -1 and 726.1 cm -1 . These states cannot be described by a single d orbital filling, as they correspond to mixtures of various configurations, in which the t2g 5 eg 2 type is predominant (assuming the labelling for octahedral symmetry). In terms of d-orbital energies, we observe that  We now turn our attention to the case of Fe III -POM, that appears to conflict with common conceptions in the field of transition-metal based SMMs. It presents a high-spin d 5 configuration, which is related with a ground sextet that is likely to mix weakly with excited states of lower multiplicities. Orbital energies corresponding to a NEVPT2 corrected, CASSCF(5,5) calculation including 1 sextet, 24 doublets and 75 singlets follow the same pattern of Mn III -POM, as the first two orbitals are a mixture of the dxz and dyz orbitals, with an energy separation of 269 cm -1 . The dxy, dz2 and dx2-y2 orbitals follow at 2594 cm -1 , 9038 cm -1 and 14825 cm -1 , respectively. Again, orbital energies resemble the ordering of a square-planar complex due to the longer Fe-O distance induced by the rigidity of the POM ligand. As expected, the ground state corresponds to a sextet state, separated by 23566 cm -1 from the lowest energy quartet. The inclusion of dynamical correlation by NEVPT2 diminishes the sextet-quartet gap to 17097 cm -1 . In both cases, the energy separation is too high to allow for an efficient state mixing by Spin-Orbit coupling, resulting in too low D values (-0.32 cm -1 and -0.37 cm -1 for CASSCF and NEVPT2, respectively). The energy splitting for the ground sextet is around 2 cm -1 for both methods, too low to explain the SMM properties of the Fe III -POM complex. Comparing the experimental and simulated magnetic susceptibility curves (see Figure 3), it is clear that the calculation is not completely capturing the magnetic anisotropy of the Fe III -POM system, as the calculated decay of T at low temperatures associated with magnetic anisotropy is too steep. This difference is reflected in the considerably higher value of the experimentally fitted zero field splitting parameter (D = -1.20 cm -1 ). 27 As mentioned earlier, the origin of this discrepancy is the predominance of sextet-quartet excitations in magnetic anisotropy, associated with an underestimation of D.

Extended active space calculations
To obtain a more accurate description of the electronic structure of the three studied complexes, we performed CASSCF(n+4,7) calculations for all models. We are especially interested in Fe III -POM and Mn III -POM, where a clear underestimation of the zero field splitting was observed in CASSCF(n,5) calculations. The extended active space now includes the bonding counterparts of the antibonding dz2 and dx2-y2 orbitals. As observed in Table 2

The case of Fe III -POM
If we consider the d-excitation energy manifold for a square-planar complex, a sizable magnetic anisotropy could be achieved if lower multiplicity states stabilize with respect to the ground sextet, yielding low energy excitations that can couple through the SOC operator. This is related with the strong destabilization of the dx2-y2 orbital, which is directly pointing towards the four equatorial ligands. In this way, S=3/2 electronic configuration will be not so disfavored by the enhanced repulsion due to electron pairing, as it will be partially compensated by the depopulation of the significantly antibonding dx2-y2 orbital. In the case of Fe III -POM, it was already commented that CASSCF results predict the lowest quartet to be 23539 cm -1 higher in energy than the ground sextet. The inclusion of NEVPT2 results in a smaller energy gap of 17063 cm -1 , which is still too high to explain the observed relaxation barrier of Fe III -POM. The calculated zero field splitting parameter modestly increased from CASSCF (-0.32 cm -1 ) to NEVPT2 (-0.37 cm -1 ).
A similar underestimation of D has been observed by Stavretis et al. for high spin Fe III porphyrin complexes coordinated to halogens. 55 The authors relate this trend to the too ionic picture of CASSCF due to its lack of dynamical correlation effects, which is only partially corrected by for Fe III -POM. 27 In our case, D values for CASSCF and NEVPT2 are similar, in contrast with results from Stavretis et al., 55 where a three-fold increase in D is observed for NEVPT2 results.
This difference is likely to be related with the larger covalency for the Fe III -porphyrin complexes in comparison with Fe III -POM. Opposed to the case of Co II -POM, we expect a smaller decrease of the relaxation barrier associated to other relaxation mechanism in Fe III -POM due to the lack of nuclear spin of the most stable isotope of iron, the negative sign of the zero field splitting parameter and its lower rhombic component. Coupling of electronic and nuclear spin has been pointed as the main factor governing magnetic relaxation in mononuclear Co compounds (I=7/2). 56 To visualize the relation between interelectronic repulsion, the sextet-quartet gap and its relation with D, a correlation diagram for Fe III -POM was constructed from ab initio ligand field parameters (see Figure 4, left). First, a CASSCF(5,5) calculation including 1 sextet, 24 quartets and 75 doublets was performed and NEVPT2 corrected energies were obtained from the former calculation. Then, Racah B and C ligand field parameters were fitted from the CI matrix  Following the evolution of the energy gap between the ground sextet and its proximal excited states ( Figure 4, left), it is clear that larger " o/B" ratio will be favorable for a lowering of this excitation energy. Furthermore, this induces a marked increase of the zero field splitting parameter when low energy excitation become proximal (Figure 4, right). In order to achieve a D value close to experiment (-1.2 cm -1 ), the ratio between ligand field and interelectronic repulsion must be larger than the obtained for NEVPT2. To put this trend in context, we can estimate the position for CASSCF calculation in the correlation diagram, being lower than the NEVPT2 value (0.82, as marked with a vertical line in Figure 4, left). Equation 1 is an empirical way to estimate the decrease in B parameter (nephelauxetic reduction) due to metal-ligand orbital mixing. B0 corresponds to the free ion value of B and h and k are empirical parameters for ligand and metal, respectively. h and k are fitted from spectroscopic data to account for the observed reduction of the Racah parameters for several complexes with different combinations of transition metals and ligands and are designed to be transferable.
Using equation 1, we estimate a B value around 700 cm -1 for Fe III -POM, considering typical values for the nephelauxetic reduction parameters (0.24 for kFeIII and 1.5 for hL, if we assume ox 2as representative of our ligand), 59 and a free ion B parameter for Fe III of 1139 cm -1 , calculated from a fit of all sextet-quartet excitations. 60 This estimated value is significantly lower than the NEVPT2 result for B (987 cm -1 ), indicating that a larger correction is required to achieve a lower Racah parameter and a higher D. In this way, it appears that NEVPT2 is correcting in the right direction, although insufficiently, from the crude overestimation of Racah parameters inherent of the CASSCF method. As discussed previously, this picture is consistent with results from Stavretis et al. for iron(III) porphyrins. 55 The behavior of CASSCF is expected due to the neglect of dynamical correlation of this method. In this way, the repulsion between electrons is overestimated and excitation energies tend to be higher than experiment. Furthermore, covalency is likely to be underestimated for the same reasons, resulting in a too low " o/B" ratio.

Conclusions
The combination of ab initio calculations and ab initio Ligand Field analysis allowed for the identification of the key features responsible for the magnetization relaxation properties of the three studied complexes. In the case of Co II -POM, SMM properties are related with the presence of a strong unquenched angular momentum in the ground state, associated to low lying d excitations. The calculated zero field splitting parameter is larger than the experimental value, in agreement with the contribution of additional relaxation mechanisms, such as the interaction with the nuclear spin (I=7/2) of the cobalt center. In the case of Mn III -POM, its lack of a measurable relaxation barrier is consistent with the absence of low lying excited states for this complex, with a relatively low zero field splitting parameter that is reasonably reproduced by ab initio methods.
Fe III -POM is the most challenging case, as it exhibits a measurable, although small, relaxation barrier for a formally isotropic d 5 configuration. ab initio results point to the axially elongated coordination environment forced by the rigidity of the lacunary POM ligand as the main factor for the unusual SMM properties of Fe III -POM. In this way, the orbital splitting pattern of the studied compounds is best described as close to square-planar instead of near octahedral. The strong destabilization of the dx2-y2 orbital results in a lowering of the sextet to quartet gap, allowing for a moderate spin mixing and the development of some magnetic anisotropy. The presented ab initio calculations underestimate the zero field splitting due to the too repulsive picture of CASSCF associated with missing dynamical correlation. NEVPT2 only partially corrects this deficiency, reflected in an overestimation of the interelectronic repulsion parameters, thus giving a lower nephelauxetic reduction. A reasonable estimation of the nephelauxetic reduction yields an energy splitting that is compatible with the presence of a small demagnetization barrier, in line with the experimentally measured value. In this way, lacunary polyoxometalates are excellent candidates for the design of highly structured and rigid coordination environments for transition metal based SMMs, favoring the development of sizable magnetic anisotropy. Concretely, POM ligands can force a specific coordination environment, overcoming coordination preferences of each metal ion.
Furthermore, its rigidity should help to restrict vibrational motion of ligands. Lacunary ligands of W POMs assure isolation of the metal ion from centers presenting nuclear spin I ≠ 0. It is likely that all this favorable features will motivate the development of new transition metal-POM SMMs. For