Intra-molecular origin of the spin-phonon coupling in slow-relaxing molecular magnets

The design of slow relaxing magnetic molecules requires the optimization of internal molecular vibrations to reduce spin-phonon coupling.

where E el is the adiabatic electronic energy and ∆X is = X is − X • is are the 3N cartesian displacements with respect to the equilibrium X • is configuration, where i spans the atom indexes and s=x,y,z. This system can still be mapped on a set of 3N decoupled 1D harmonic oscillators by introducing the normal mode of vibration. We first start defining mass-weighted cartesian coordinates u a = √ m i ∆X is , where the single index a=3(i-1)+s runs over the 3N degrees of freedom. Diagonalizing the force-constant matrix of the energy second-order derivatives H it is possible to define normal mode of vibrations q a and their frequencies ω a : where given H eigenvectors L ab and eigenvalues diag(H) aa The inverse transformation that defines the cartesian displacement associated to a unit-less normal modē q a amount of displacement is therefore Cartesian displacements so obtained can now be decomposed in three contributions: translational, rotational and internal displacements. To do that we followed the procedure outlined by Neto et al. 1 , where Eckart-Sayvets conditions are imposed on molecular atomic displacements through a self-consistent numerical procedure. The first step of this procedure consists in the definition of translational (T α ) and rotational (θ α ) displacements. This new set of coordinates however makes the metric tensor different from the identity matrix and both covariant and controvariant coordinates must be employed. The Einstein convention on repeated indexes will be used in the rest of the section.
where (I αβ ) 0 is inertia tensor calculated in the centre of mass cartesian reference system and M α , with α = 1 − 3, are infinitesimal rotation matrices around the cartesian axis α: In this framework, atomic positions with respect to the centre of mass x is = X is − X Bs can be calculated from the initial ones x • is = X • is − X • Bs by summing the pure internal contributions ∆X int is and rotating the resulting coordinates: where Λ express a general rotation of an angle ψ around a vector ξ . The amount of rotation in terms of Λ and translations can be determined by imposing the Eckart-Sayvets conditions to the internal displacements, that read This set of equations impose the independence of internal and external degrees of freedom, requiring that the projection of one set of coordinates onto the other is null. Indeed the term in parenthesis is nothing but the Jacobian matrix that transforms the complete set of 3N coordinates into the 6 external coordinates. Eqs.7 can then be used to compute the translational and rotational contribution of the full ∆X is displacements, then the internal ones can be computed by difference through Eq. 9: where Λ can be calculated employing Eq. 10 and the relations This process should be cycled using the calculated ∆X int is as new input displacements in Eq. 11 until the equation get verified. This procedure has been used to compute the amplitude of local translation, local rotation and internal vibrations of a single [(tpa Ph )Fe] −1 molecule inside its crystal, as showed in the main text. To do so, we inserted in Eq. 11 only the 3N cartesian displacements corresponding to this specific molecule coming from normal modes, which instead describe the vibrational motion of all the primitive-cell atoms. These quantities differ from the usually defined acoustic, librational and optical modes of a lattice as the latter are defined as translations, rotation and internal vibrations of all the atoms inside the unit-cell, with no particular distinction among inter/intra molecular motions. To summarize, one starts with the normal modes provided by the DFT calculation. Next, by applying the projector technique to the normal mode cartesian displacements ∆X tot is of Eq. 5, one ends up with three additional kind of cartesian displacements: ∆X int is , ∆X rot is , and ∆X trasl is . By separately plugging these cartesian displacements in Eq. 4 one can define the normal modes projection onto the translational, rotational and intra-molecular space For completeness we here also report in Fig. 1 the contributions of acoustic, librational and optic motions as function of frequency. As expected, acoustic motions are absent at any frequency as we are looking only at Γ-point vibrations.
Spin-phonon coupling coefficients are defined as the first order derivatives of the spin Hamiltonian parameters with respect to the normal mode of vibrations ∂ D i j /∂ q α 2 . The strategy we employed to compute them starts with the evaluation of the numerical D tensor derivatives with respect to the cartesian coordinates of [(tpa Ph )Fe] −1 . To do that we scanned one coordinate X is + δ X is at the time with steps δ X is equal to ±0.0025, ±0.002, ±0.0015, ±0.001 and ±0.00075. The resulting points have then been interpolated with a second order polynomial expression in order to estimate the linear terms, that correspond to the ∂ D i j /∂ X ks coefficients. Here we would like to report a few examples to show the details of the method. Fig. 4 reports the scanning of the six independent D elements along the x component of the iron atom.  It is interesting to note how, in the case of iron coordinates derivatives, the second order component of the polynomial expression is fundamental to obtain a good fitting. This is not true however for all the other elements and for atoms outside the first coordination shell a linear expression is enough to obtain a good regression. The breaking down of a first order perturbation theory for large displacements of the first coordination shell clearly comes from the fact that spin-phonon coupling interaction is magnified in proximity of the iron d-shell electrons. For the same reason, increasing the distance between the displaced atom and iron, the value of spin-phonon coupling decreases and its evaluation become more affected by numerical noise. We decided to exclude D derivatives with an error on the linear term higher than 5%. As 1-8 | 5 an extreme example of excluded contribution, Fig. 5 reports the scanning of D elements along the y coordinate of a pyrrolide carbon not directly bounded to the nitro atom. Finally, the cartesian derivatives of the anisotropy tensor D have then been projected on the basis of normal modes by means of their cartesian coordinates definition in Eq. 5.
This procedure for the computation of ∂ D i j /∂ q a has a double advantage. The first one regards the number of CASSCF calculations needed. Assuming the electronic structure of a single [(tpa Ph )Fe] − units to be independent from the other molecules in the crystal, the calculation of the D tensor derivatives can be done directly on a single molecule without periodic boundary conditions 2 . The number of independent displacements for a given molecule corresponds to its 3N degrees of freedom but the number of normal modes instead grow like the number of the primitive cell degrees of freedom, which is generally much larger than 3N. Moreover, Eq.17 can be used also to compute all the separate contributions to spin-phonon coupling coefficients by simply replacingL ba with its own projection on specific displacements, calculated as explained in the previous section. This last step can be taken as many times as one wants as there are no additional computational costs beyond the ∂ D i j /∂ X is coefficients and phonons calculation already done. Fig. 6 reports the comparison between the spin-phonon coupling coefficients computed with the just outlined procedure and those obtained by the direct differentiation with respect to the unit-cell normal modes 2 . The two calculations show only minor differences, validating the consistency of both calculations. One final remark regarding the possibility to project the spin-phonon coupling coefficients in the basis set of internal coordinates. To accomplish that the weighting coefficients of Eq. 17 must be chosen as (∂ X is /∂ q n ), where q n s represent controvariant internal coordinates such as stretchings bendings and torsions. The q n s are defined as function of cartesian coordinates and therefore the coefficients (∂ X is /∂ q n ) can be otbained by inversion of the Wilson's matrix G n,is = (∂ q n /∂ X is ). Although these coefficients are mathematically well defined this procedure hides a potential inconsistency. Indeed, due to the non-orthogonality of internal coordinates, this coefficients depends on the definition of the basis itself. For instance the computed derivative of a bending angle coordinate with respect to cartesian coordinates will 1-8 | 6 be different depending on the set of 3N-6 internal parameters used to define the intra-molecular basis set.
A possible solution to this problem regards the factorization of the metric tensor of the internal coordinates in order to define a new set of internal orthogonal parameter 3 . However, doing so the chemically intuitive picture in terms of stretchings and bending would be lost as one would be scanning the D tensor along a non-trivial combination of them. Consequently, we decided to keep working we the usual internal molecular coordinates and proceed as follows. Since we are interested in the Iron first coordination shell we selected three atoms at the time and compute the (∂ X is /∂ q n ) coefficients for the two stretchings and one bending. Doing so we obtain one ∂ D i j /∂ q n coefficients for each bending and three coefficients for each stretching, which have been avaraged among them. Using the smallest possible number of internal vectors at the time is then possible to reduce the uncerteinity on the computed coefficients due to the non-orthogonality issue. It must be stressed that the order of magnitude of the ∂ D i j /∂ q n coefficients is conserved within different basis sets. However, the correctness of this approach should be found in the coherence of the data presented in the main text.