Insights into phonons and spin–lattice relaxation in copper(II) and vanadyl(IV) porphyrin metal–organic frameworks from density functional theory

Abstract

Employing metal–organic frameworks (MOFs) has been shown to be an effective strategy for extending spin–lattice relaxation times (T₁) by modifying the vibrational properties of molecular spin qubits. Although a previous study employed terahertz spectroscopy to probe vibrational properties, the specific types of vibrational modes that affect T₁ have remained unclear. In this study, we use periodic density functional theory calculations to investigate the vibrational properties of the MOF [{M(TCPP)}Zn₂(H₂O)₂] (M = Cu, VO, and TCPP = tetrakis(4-carboxyphenyl)porphyrin) and the corresponding molecular crystals MTPP (TPP = tetraphenylporphyrin). Although the longer T₁ of the vanadyl MOF compared to the VOTPP crystal has been attributed to the absence of low-frequency modes, our combined experimental Raman spectra and phonon simulations show that the mere presence of low-frequency modes does not necessarily lead to faster relaxation times. To rationalize the differences in T₁ between the MOFs and the corresponding molecular crystals, we calculated the spin–phonon coupling (SPC) for each Γ-point phonon mode. Furthermore, we analysed correlations with the magnitude of the SPCs, the symmetry of the modes, and the in-plane and out-of-plane distortion of the porphyrin framework. Our results reveal that no single descriptor (frequency, symmetry, or distortion magnitude) can reliably predict the SPC strength, highlighting the multifactorial nature of the SPCs in these systems. This complexity underscores the importance of explicit computational treatments for identifying the key phonon modes that drive spin–lattice relaxation, while spectroscopic techniques such as low-frequency vibrational spectroscopy can provide complementary validation and qualitative insights.

---

Introduction

Electronic spins in molecular materials are gaining attention as qubits for quantum computing and quantum sensing due to their chemical tunability and the potential for spatial control through molecular design.^1–4 Despite these advantages, such systems still suffer from relatively short coherence times (T₂),⁴ and only a few have exhibited noticeable coherence at room temperature.^5–7 T₂ represents the duration over which spins can maintain their coherent quantum superposition states and determines the maximum time window for performing quantum operations.⁴ Decoherence is typically driven by interactions with nearby spins at low temperatures (below ∼15 K), while at higher temperatures, spin–lattice relaxation becomes the major limiting factor.⁸ This process, characterized by the spin–lattice relaxation time (T₁), governs how efficiently an excited spin state exchanges energy with its surrounding lattice environment. This spin–lattice relaxation time sets a fundamental upper limit on coherence, as 2T₁ ≧ T₂.⁹ Consequently, understanding and controlling T₁ relaxation pathways is essential for improving the coherence times of molecular spin qubits, particularly at higher temperatures.

Spin–lattice relaxation occurs mainly via two mechanisms: phonon-mediated processes and processes arising from spin–spin interactions (cross relaxation). In a magnetically well-diluted system and at higher temperatures,⁸ the phonon-mediated process becomes the primary relaxation mechanism. This is because the interactions between paramagnetic centres become negligible, and the cross relaxation is temperature-independent.¹⁰ In the phonon-mediated mechanism, different processes dominate in different temperature ranges. At low temperatures (below 20 K),¹¹ the one-phonon-mediated direct process is typically dominant. In this process, the energy of the spin is directly exchanged with the phonon energy, which requires a precise energy match of the spin and the phonon. In typical electron paramagnetic resonance (EPR) experiments performed at X or Q-band, the energy difference of the total spin quantum number S = 1/2 (with magnetic sublevels m_s = ±1/2) is less than a few cm⁻¹, therefore largely limiting the available phonons for energy exchange. As the temperature increases, the two-phonon Raman process becomes dominant. The Raman process does not require a direct match of the phonon energy and the spin energy; only the difference between two phonon modes and the spin energy must match. This significantly loosens the requirement for spin–lattice relaxation and allows phonons with frequencies larger than a few cm⁻¹ to participate in the process. Combined with the increased population of higher-energy phonons with increasing temperature, the two-phonon processes become the dominant relaxation mechanism typically above 20 K.¹² Another crucial relaxation mechanism is the local process. This process is also classified as a two-phonon process, similar to the Raman process,¹¹ but involves optical phonons or vibrations with localized character and is typically identified as the major relaxation mechanism at higher temperatures.^13,14

Given that these phonon-mediated processes, particularly the Raman and local mode processes, make dominant contributions to the spin–lattice relaxation, tuning the vibrational spectrum of molecular qubits can be an effective strategy for suppressing spin–lattice relaxation at elevated temperatures. Indeed, reducing low-frequency vibrations has been shown to be an effective strategy to mitigate spin–lattice relaxation. A common approach to suppress such low-frequency vibrations is to create a rigid ligand environment, thereby shifting the vibrational modes around the spin centre to higher energy.^13,15 In addition to rigidifying the ligand environment, the vibrational spectrum can also be tuned by embedding molecular units into an extended framework. This approach, which is the focus of this study, allows control over lattice vibrations through the surrounding host structure. In this context, a previous study⁷ demonstrated that incorporating molecular units into the lattice of a metal–organic framework (MOF) can prolong the spin–lattice relaxation times, particularly at elevated temperatures. Terahertz (THz) spectroscopy suggested that the MOF lattice shifts optical phonons in the corresponding molecular crystal to higher energies, thereby reducing the efficiency of the spin–lattice relaxation. While the study by Yamabayashi et al.⁷ established the beneficial effect of framework incorporation on spin–lattice relaxation, the underlying vibrational contributions remained unexplored. Computational studies^16,17 on spin–phonon coupling in vanadyl and copper molecular spin qubits in the gas phase suggest that only a few symmetry-allowed key vibrations have a significant impact on the spin–lattice relaxation time. Therefore, understanding how the MOF framework modifies these key vibrations in molecular solids is crucial for the rational design of MOFs with improved spin–lattice relaxation times.

In order to better understand the difference between the vibrational properties of qubits in MOFs and molecular crystals, we investigated prototypical MOFs with copper(II) and vanadyl(IV) porphyrins ([{Cu(TCPP)}Zn₂(H₂O)₂] (TCPP = tetrakis(4-carboxyphenyl)porphyrin) (1) and [{VO(TCPP)}Zn₂(H₂O)₂] (2)). These MOFs have previously been shown to have longer T₁ than their corresponding molecular crystals of CuTPP (TPP = tetraphenylporphyrin) (3) and VOTPP (4)^7,18–20 (Fig. S1) and serve therefore as an ideal model system for a comparative analysis (Fig. 1). The spin–lattice relaxation process is strongly influenced by both the electronic structure and the phonons. However, several studies^7,18,19 have shown that the g factors are very similar between the MOF-based porphyrins and their molecular crystal analogues (Tables S1 and S2). The g factor describes how an electron's magnetic moment couples to an external magnetic field and is closely related to the electronic structure. The deviation of the g factor from the free-electron value (g_e ∼ 2.0023) occurs due to spin–orbit coupling and the ligand environment. The magnitude of deviation of g factors from the free electron value is known to be a useful proxy for the size of the spin–orbit coupling.¹¹ Thus, similar g factors suggest that phonons are the primary source of their distinct spin–lattice relaxation behaviours rather than the difference in their electronic structures.

Fig. 1 The crystal structures of the copper(II) and vanadyl(IV) porphyrin systems studied in this work visualized using Mercury.²² View along the a axis of ([{Cu(TCPP)}Zn₂(H₂O)₂]) (P2/m) (1) and [{VO(TCPP)}Zn₂(H₂O)₂] (2). View along the c axis of CuTPP (TPP = tetraphenyl porphyrin) (I4̄2d) (3) and VOTPP (I4) (4). The unit cells are outlined with thin lines. Colour code: Cu – orange, V – silver, Zn – purple, N – light purple, O – red, C – grey, H – omitted for clarity.

To probe this, the Γ-point phonons of systems (1)–(4) were calculated and compared using periodic density functional theory (DFT), as simplified single molecule or cluster models fail to reproduce low-frequency (terahertz) vibrations because they are highly sensitive to intermolecular and lattice interactions.²¹ Since the vibrational modes that dominate spin–lattice relaxation at relevant temperatures (room temperature or below) are low-frequency modes below 400 cm⁻¹,¹⁶ where intermolecular interactions play a significant role, a precise description of these phonon modes is crucial. In addition to phonon calculations, spin–phonon coupling analysis was conducted in order to evaluate how these vibrational modes interact with the spin degrees of freedom and to understand the origin of their different spin-relaxation behaviours.

Methods

Raman spectroscopy

Unpolarized Raman spectra of the polycrystalline samples were measured using a home-built single-monochromator micro-Raman spectrometer with a 488 nm solid laser, as described elsewhere.²³ The materials were prepared using established methods.^7,18,19,24 The details are provided in the SI. The spectrometer is equipped with two volume holographic notch filters (SureBlock) from Ondax Inc. (currently Coherent Corp.)²⁵ that allow measurements down to ∼10 cm⁻¹.

Raman spectra were collected at wavelengths of 488 nm and 532 nm with an NRS-5100 spectrometer to confirm that wavelength-dependent resonance effects did not influence the relative band intensities, enabling a direct comparison with simulated non-resonant spectra. At 532 nm, the measurable range was limited to >100 cm⁻¹ due to the cutoff of the notch filter required to suppress Rayleigh scattering.

Computational methods

Periodic DFT calculation (cell optimization and Γ-point phonon simulations). Full geometry optimizations including the lattice constants as well as atomic positions, and subsequent phonon calculations at the Γ-point were performed using either CRYSTAL23²⁶ or VASP (version: 6.4.1)^27–31 in conjunction with Phonopy (version: 2.38.1).^32,33 The key parameters (basis sets, plane-wave energy cutoffs, k-point grid) that determine the quality of simulations were chosen based on convergence tests. For the periodic hybrid functional calculations, the divergence of the exact exchange term at the Γ point was treated using the standard finite-size correction schemes implemented in the respective codes. In VASP, the Coulomb singularity of the exchange interaction at the Γ point is handled by the Gygi–Baldereschi correction. In CRYSTAL23, the exchange divergence is treated analytically within the exact-exchange formalism using the auxiliary-function approach implemented for periodic boundary conditions. Details of these convergence tests and additional computational settings are described in the SI (Tables S3, S4 and Fig. S6, S7). Initial structures were taken from experimental single-crystal X-ray data reported in the Cambridge Crystallographic Data Centre (CCDC), when available. Hydrogen atoms were added where missing, while maintaining the reported space group. For system (1), hydrogen atoms were added to the oxygen atoms of the H₂O molecules coordinating to the Zn nodes, initially preserving the original space group of P2/m from the reported single-crystal X-ray structure (CCDC: 1555581).¹⁸ This resulted in an imaginary frequency mode localized on the H₂O molecule. To eliminate this mode, the symmetry was reduced to P2 by removing the mirror plane around the H₂O molecule.

For (2), no single-crystal X-ray structure has been reported, as the structure determination was hampered by their twinning or polycrystalline nature.¹⁹ However, the previous work established that this system is assumed to be isostructural with (1).¹⁹ Therefore, its structural model was constructed by substituting Cu with an oxovanadium (VO) unit in (1). Because the VO unit breaks the twofold axis present in (1), the symmetry is reduced to P1, effectively lowering the lattice system from monoclinic to triclinic. To avoid artificial distortions of the lattice during structural optimization of (2), the lattice constants were kept fixed in this case, and only the atomic positions were relaxed.

For (4), the orientation of the VO group is disordered in the experimental single-crystal structure with the I4/m space group (CCDC: 1125976³⁴ or 2290193⁷). Explicitly accounting for this disorder would require large supercells and is therefore not practical. To approximate the two limiting configurations, the unit cell was optimized in the P4/n (mirror plane lost, inversion preserved via a glide plane) and I4 (both inversion and mirror symmetry lost) space groups, as shown in as shown in Fig. S8. Geometry optimization and phonon calculation using the P4/n cell resulted in two imaginary frequency modes (∼−25 cm⁻¹). Therefore, the structure with the I4 unit cell was employed for further discussion, and the symmetry analysis was only conducted based on this space group.

For the simulations with CRYSTAL23, the frequency calculations (phonon calculations at the Γ-point) were carried out based on the fully optimized structures (including lattice relaxations). Then, the computed frequencies were checked to ensure there are no imaginary modes to confirm that the optimized structure is stable. Subsequently, the intensities of the infrared (IR) and Raman modes were calculated in order to compare the spectra with our experimental measurements. Raman intensities were computed at the temperature of measurement (298 K) and the wavelength of the incident laser (488 nm and 532 nm). For simulations with VASP, phonon calculations were carried out using the finite displacement method implemented in Phonopy, employing the primitive unit cells.^32,33 The amplitude of the atomic displacement was set to 0.01 Å. The convergence of the computed phonon modes was confirmed by the agreement of the eigenvectors (vide infra).

Spin–phonon coupling analysis. The spin–phonon coupling (SPC) coefficients were numerically calculated with ORCA (version: 6.0.0),^35,36 using the simulated phonon modes from the periodic DFT calculations. The eigenvectors obtained from the periodic DFT calculations were projected onto the optimized structures to generate distorted geometries along individual phonon modes. For modes simulated with CRYSTAL23, the displacements along the normal coordinates were performed using the implemented SCANMODE module, with step sizes of (Q = ±0.005, ±0.01, ±0.015). In the convention used in CRYSTAL23, Q is a dimensionless parameter, and a displacement of Q = 1 corresponds to the amplitude that yields a potential energy of 1/2ħω in the harmonic potential.³⁷ For phonon modes obtained from VASP, the distorted geometries were generated using the Phonopy eigenvectors, with the same step sizes and Q convention as applied in CRYSTAL23. Then, the SPC coefficients calculations were performed following the methods reported in ref. 13, 16, 17 and 38–41. For these calculations, the molecular porphyrin units (CuTPP and VOTPP) were cut out from the optimized and distorted crystal structure in (3) and (4). For (1) and (2), the MTCPP⁴⁻ (M = Cu/VO) fragment was cut out from the optimized and distorted periodic structure, similar to a previous study.³⁹ For the g factor calculations, the hybrid functional PBE0⁴² was used as this type of functionals are commonly reported to yield reliable g factors.^43,44 The relativistic corrections were included through the ZORA method.⁴⁵ The basis set ZORA-def2-VTZP⁴⁶ was used for the metal centre (Cu and V) as well as the coordinating N atoms, while the smaller basis set of ZORA-def2-SVP⁴⁶ was employed for the remaining atoms. Moreover, the SARC/J auxiliary basis set⁴⁷ were utilized, similar to previous studies.^13,16 Full computational details on the remaining settings are provided in the SI.

Normal-coordinate structural decomposition analysis. For the evaluation of out-of-plane displacements of the porphyrin ligands, normal-coordinate structural decomposition analysis (NSD)^48–50 was carried out using the PorphyStruct program.⁵¹ To analyse the out-of-plane displacements of the porphyrin ligand for each vibrational mode, displaced geometries were generated at Q = 1, which corresponds to the maximum displacement of the ground-state vibration in the harmonic potential, based on the modes computed in the phonon calculations.

Results and discussion

Raman spectra comparison

Periodic DFT calculations were carried out to simulate the Raman spectra. The 488 nm experimental measurements required two notch filters to suppress the Rayleigh scattering, which reduced the overall signal intensity and lowered the signal-to-noise ratio compared with some of the 532 nm spectra (Fig. 2 and Fig. S10–S12). Although some variations exist in the intensity depending on the wavelength, the overall spectral features are consistent between the two measurements. The calculated spectra show qualitatively good agreement with the experimentally obtained Raman spectra (Fig. 2 and Fig. S10–S12).

Fig. 2 Experimental Raman spectra (light blue) (488 nm) and blue (532 nm) of Cu MOF (1) and simulated spectra at 488 nm using the following DFT methods: PBEsol-D3 with plane wave basis sets (energy cutoff = 850 eV) (violet), PBEsol03c/sol-def2-mSVP (magenta), and PBE-D3/VTZP (pink). The linewidths of the computed vibrational frequencies were modelled using Lorentzian broadening using a full width half maximum value of 5 cm⁻¹.

MOF are well known for their ability to adsorb solvent or gas molecules within their pores. In passing we note that in the simulations presented in this study, we used MOF models without any gas or solvent molecules in the pores. We present several justifications for this choice. First, the MOFs synthesized were solvent exchanged in acetone and then dried under vacuum. For EPR experiments to measure the relaxation time, these solid samples were sealed under high vacuum in a glass tube. This should effectively remove weakly bound molecules from the pores. Raman spectra were measured under ambient condition. In general, the phonons or vibrational modes in the THz region is very sensitive to any structural changes. Some MOFs have been reported to show noticeable difference in their vibrational spectra depending on the existence of absorbent molecules in the pore.^52,53 The agreement between the experimental spectra and the simulated spectra in this study without considering any absorbents suggest that the model in the simulation is sufficiently accurate.

Impact of the basis set and density functional on phonons

Simulations of Raman spectra were initially carried out using CRYSTAL23,²⁶ which employs atom-centred Gaussian basis sets. The extensive use of symmetry operations and the efficient handling of atom-centred Gaussian basis sets in this program make its calculations efficient.⁵⁴ However, the large variety of available basis sets makes it rather difficult to systematically improve the quality of the calculations and achieve convergence. For this reason, we also tested VASP, which uses plane-wave basis sets, as the quality of the basis set can be controlled by a single parameter – the kinetic energy cut-off (ENCUT).

For the simulations using CRYSTAL23, various basis sets were tested ranging from double-ζ-quality basis sets (sol-def2-mSVP,⁵⁵ 6-31G**^56–59) to triple-ζ basis sets (POB-TZVP,⁶⁰ VTZP⁶¹). Among these, the composite method PBEsol0-3c/sol-def2-mSVP⁵⁵ reproduced the spectral feature and lattice constants very well (Fig. 2, Fig. S10–S12 and Tables S5–S7). In passing, we note that PBEsol0-3c⁵⁵ is a solid-state composite method based on the PBEsol0 hybrid functional,⁶² combined with a compact sol-def2-mSVP⁵⁵ basis set and corrections for dispersion interactions on the D3 level⁶³ using the Becke–Johnson damping function⁶⁴ and basis-set superposition error. The method is known for its efficiency while maintaining good accuracy.⁵⁵

However, subsequent eigenvector analysis (Fig. S13, SI) revealed that the basis sets of double-ζ-quality are not sufficient in reaching convergence, especially at low-frequency modes. In this analysis, the similarity of the eigenvectors (displacement vector of each atom) of the phonon modes computed with different functionals or basis sets were compared. The similarity was evaluated by taking the dot product of normalized displacement vectors for each atom.⁶⁵ When two phonon modes from different calculations are equivalent, their normalized dot product approaches 1; conversely, values close to 0 indicate little or no similarity. Fig. S13b and d shows that the dot products with a well converged VASP calculation (PBE-D3) of the eigenvectors are less than 0.5 for various modes at low frequency, indicating that the phonon modes computed with double-ζ-quality basis sets (sol-def2-mSVP and 6-31G**) are not consistent. The convergence of VASP calculations were assumed from an excellent agreement of calculations with two different functionals (PBE-D3 and PBEsol-D3), as shown in Fig. S13a. Simulations with double-ζ-quality basis sets resulted in discrepancies in vibrational modes that show significant spin–phonon coupling coefficients in the subsequent spin–phonon coupling analysis (vide infra, Fig. S14 and S15). In order to gain consistent results across different methods, the VTZP basis set, which has previously been used for a Cu-based MOF and showed reasonable convergence,⁶⁶ was required.

Upon increasing the basis set to this level, the agreement of the eigenvector overlaps across different calculations significantly improved. As shown in Fig. S13e and f, the PBE-D3/VTZP level calculations performed using CRYSTAL23 and VASP (PBE⁶⁷ or PBEsol⁶² with ENCUT = 850 eV) gave nearly the same results for the phonon calculations and subsequent spin–phonon coupling calculations (vide infra, Fig. S17 and S18), validating the convergence of our phonon calculations. POB-TZVP⁶⁰ is also a triple-ζ-quality basis set, it is known to give inconsistent results and tends to overestimate intermolecular forces, leading to contracted lattice parameters and blue-shifted vibrational modes.⁶⁸ In line with this observation, the POB-TZVP basis set performed poorly here as well, producing the largest deviation in the lattice constants (Tables S6 and S7) and phonon eigenvector overlaps with those obtained from other methods (Fig. S13).

Spin phonon coupling calculations

Phonon calculations provide insight into the distribution of phonon modes in MOFs and molecular crystals. However, the phonon modes and their frequencies alone do not directly reveal which vibrations dominate their spin–lattice relaxation. To identify the key contributors to spin-lattice relaxation, we quantified the coupling of each phonon mode to the spin using Γ-point phonons obtained from periodic DFT calculations, following established protocols.^{13,16,17,38–41} Although all phonon modes across the Brillouin zone contribute to the spin–lattice relaxation, considering only Γ-point phonon modes has been shown to provide sufficient accuracy without significantly affecting simulation results.¹¹

Within the spin-Hamiltonian framework, several approaches exist for computing SPC coefficients, depending on whether the g-tensor (g; Zeeman interaction) or the hyperfine tensor (A; hyperfine interaction) is considered to be the leading spin–lattice relaxation mechanism, and whether first or second derivatives with respect to normal coordinates (Q) are taken. Previously reported methods include the first derivative (ref. 13, 16, 20, 38–41 and 69–72), (ref. 41 and 73), the second derivative (ref. 17), second mixed derivative (ref. 11 and 73–75) and (ref. 73 and 74). The choice of g or A to be differentiated depends on the assumed dominant mechanism of spin–lattice relaxation, which remains under debate.^11,16 Despite this ongoing debate, a previous study has shown that g and A derivatives yield consistent trends in spin–phonon coupling behaviour, giving essentially the same vibrational modes as the dominant contributors.⁷³

In this study, we employed the more widely adopted g-tensor approach to gain insights into the modes that have a substantial contribution to spin–lattice relaxation. Our primary objective is to rationalize the different relaxation behaviours of molecular crystals and MOFs bearing the same molecular units, rather than reproducing absolute T₁ values. The g-tensor framework provides a particularly transparent connection between vibrations and magnetic relaxation, as g-tensor anisotropy arises directly related to spin–orbit coupling.¹⁶ We note that a complete and quantitative description of spin–lattice relaxation generally requires consideration of both g- and A-tensor modulation. However, because the g or A-tensor approach is reported to give the same trend,⁷³ the choice between the g or A-tensor approach is unlikely to impact our qualitative final conclusions.

Results of SPC analysis

For the calculation of the g-factors, two hybrid functionals (PBE0⁴² and B3LYP^76,77), two basis sets (def2-TZVP and def2-SVP⁴⁶) as well as two relativistic correction methods (ZORA⁴⁵ and Douglas–Kroll–Hess method (DKH)^78,79) were tested to confirm that they give consistent results. The trend of g-factor changes as a function of the normal coordinate Q and was shown to be largely independent of the choice of functional, basis set, and the relativistic correction method (Fig. S16–S18). These calculations confirmed that these computational settings have only a minor influence on the results. For the preceding calculations, the functional PBE0 and the basis set def2-TZVP for the metal and coordinating nitrogen, and def2-SVP for the remaining atoms were used with the ZORA method, following a previous study.¹⁶

First-order SPC

First-order spin–phonon coupling (SPC) governs one-phonon relaxation mechanisms, namely the Orbach and direct processes.^69,80 The Orbach process is a one-phonon relaxation process similar to the direct process, but through excited electronic states. For S = 1/2 copper(II) and vanadyl(IV) qubits, the first excited electronic state typically lies far above the thermally accessible range. Therefore, this process is considered negligible. The direct process is a one-phonon process, in which the energy of the spin is directly exchanged with the phonon. Since a precise energy match of the spin energy and the phonon energy is required, acoustic phonons are the main contributors to this process. The direct process is dominant only at extremely low temperatures (below 20 K).¹¹ Therefore, this process is outside the scope of this study, which focuses on optical phonons and spin–lattice relaxation at higher temperatures. However, first-order spin–phonon coupling coefficients are often employed to analyse multi-phonon relaxation (Raman and local relaxation processes),^13,16,20,69 under the assumption that first-order coupling strengths serve as a reliable proxy for the second-order spin–phonon coupling.^16,69

Fig. 3 presents the first-order SPC coefficients of the systems (1) and (3) determined by fitting g-values as a function of atomic displacement for each phonon mode. First-order spin–phonon coupling coefficients for the systems (2) and (4) are calculated similarly and shown in Fig. S27. Consistent with previous findings,^16,40,70 only a small fraction of vibrational modes exhibits significant spin–phonon coupling, while the majority contribute negligibly. For (1), vibrational modes around 200 cm⁻¹ and 380 cm⁻¹ show the largest first-order SPC coefficients (Fig. 3) due to both correct symmetry and large displacement of the ligand (Table S13, see Discussion). In the case of (3), the dominant contributions arise from modes around 400 cm⁻¹ similar to (1), with additional smaller contributions from around 120 cm⁻¹, 200 cm⁻¹, and 330 cm⁻¹ modes (Fig. 3). For the vanadyl systems, both (2) and (4) have major first-order spin phonon coupling coefficients around 20–30 cm⁻¹, 100 cm⁻¹, 200 cm⁻¹, 300 cm⁻¹ and 400 cm⁻¹ (Fig. S27). System (2) also exhibits additional contributions near 250 cm⁻¹ and 350 cm⁻¹. While small variations exist in the precise frequencies and number of contributing modes, the overall trend looks similar, with major contributions clustered around comparable frequency ranges (highlighted in green in Fig. 3).

Fig. 3 First-order SPC coefficients for (a) Cu MOF (1) and (b) CuTPP (3). Vibrational modes were computed using VASP with the PBE-D3 functional. Phonons that show sizable SPC coefficients at similar frequencies in (1) and (3) are marked in green, and others are marked in blue.

The magnitudes of the SPC coefficients are also comparable between the MOF and molecular crystal counterparts. For (1), the largest coefficients are approximately 4 × 10⁻⁶ (g_x and g_y) and 5 × 10⁻⁵ (g_z) (Fig. 3). For (3), the maximum values are about 1.5 × 10⁻⁶ (g_x and g_y) and 1.5 × 10⁻⁵ (g_z). A similar trend is observed for the vanadyl systems of (2) and (4): the dominant SPC coefficients are on the order of 10⁻⁸ for g_x and g_y, and 10⁻⁶ for g_z (Fig. S27). These similarities stand in contrast to the experimentally observed differences in spin–lattice relaxation between the MOF and molecular crystal systems. Typically, experimental T₁ values associated with parallel components (g_z) are larger, indicating a weaker spin–phonon coupling compared to that measured in the perpendicular (g_x and g_y) direction.⁸¹ This behaviour, however, contradicts the calculated trend, where g_z (∼10⁻⁶) show larger SPC coefficients than g_x and g_y (∼10⁻⁸), for example in (2) and (4).

SPC coefficients set the intrinsic coupling strength, but phonon populations also control each mode's contribution to the spin–lattice relaxation at finite temperatures. To incorporate phonon population effects, SPC coefficients were scaled using the Bose–Einstein distribution (Fig. S29–S32). As shown in Fig. 4a, the theoretical temperature dependence of a value proportional to T₁ derived from the summed contributions of all vibrational modes is nearly identical for (1) and (3), in apparent contrast to the experimental results (Fig. S1). For a quantitative comparison with the experiment, an empirical arbitrary factor is used to scale to the experimental value in the previous study.¹⁶ Since this method is a simplified method¹⁶ from the more complicated but rigorous method,⁸⁰ we avoid a quantitative argument of T₁ here but restrict our argument for only the qualitative trend. Nevertheless, the discrepancy between the computational and experimental results is clearly seen even at this qualitative level.

Fig. 4 Temperature dependence of T₁ in Cu MOF (1) and CuTPP (3) calculated by scaling (a) the first-order SPC coefficients with the Bose–Einstein temperature factor and (b) the second-order SPC coefficients with the Bose–Einstein temperature factor.

Second-order SPC

This discrepancy between the theoretical result based on the first-order spin–phonon coupling and experimental observation may reflect limitations of the approach. As noted, the first-order coupling formally accounts only for the first-order processes such as the direct process and Orbach process, despite the common use of this model in previous studies due to its computational simplicity.^{13,16,20,38–41,69–72} To consider a more complete picture, we also considered second derivatives of the g-tensor in addition to the first-order coupling coefficients. However, evaluating cross-coupling terms (mixed partial derivatives) presents significant computational challenges, requiring an enormous number of calculations (e.g., at least 10⁵ for VOTPP) as the number of calculations scales quadratically with the number of vibrational degrees of freedom.⁷³ Accordingly, only the non-cross terms were evaluated in this study. The non-cross terms describe a relaxation process, in which the absorbed and emitted phonons belong to the same vibrational mode or phonon branch. Although such processes would be forbidden by strict energy conservation for single-frequency optical phonons, in real crystals, optical phonons have linewidths due to finite dispersion across the Brillouin zone, and finite lifetimes. Both effects enable this relaxation channel mediated by optical phonons.^14,82

Fig. 5 presents the computed second-order (non-cross terms only) spin–phonon coupling coefficients. In contrast to the previous assumption that the first-order SPC shows a similar trend with the second-order SPC,^16,69 our results reveal markedly different distributions of coupling strengths between the first- and second-order coefficients across all systems investigated (Fig. 3, 5 and Fig. S27, S28). Moreover, contrary to the first-order SPC, where the differences between the MOF and molecular crystal systems are relatively minor, the second-order SPC reveals clear distinctions. For the Cu(II) species, (1) exhibits markedly smaller coefficients, on the order of 10⁻⁹, compared to (3), which shows coefficients on the order of 10⁻⁵ for the g_x and g_y directions (Fig. 5). Similarly, for the vanadyl case (Fig. S27 and S28), the MOF (2) displays SPC coefficients on the order of 10⁻⁸, whereas (4) reaches values around 10⁻⁶ for g_x and g_y. Reflecting these differences in SPC magnitudes, the temperature-scaled SPC coefficients (Fig. S29–S32) also differs significantly between the systems, showing better agreement with experiments: (3) and (4) exhibit faster spin–lattice relaxation than (1) and (2). Fig. 4b shows the sum of the contributions to T₁ from all vibrational modes scaled by the Bose–Einstein phonon occupation distribution. In this case, a better agreement was obtained with the experimental data that shows a less steep decrease in T₁ for (1) compared to (3) than in the case of the first-order calculation. The improved agreement obtained with the second-order approach suggests that identifying the vibrational modes relevant to T₁ at high temperatures – beyond those involved in the direct process – requires second-order spin–phonon coupling coefficients rather than only first-order terms.

Fig. 5 Second-order SPC coefficients for (a) Cu MOF (1) and (b) CuTPP (3). Vibrational modes were computed with VASP using the PBE-D3 functional.

Although this second-order treatment improves the agreement with experiment, some discrepancies remain. These differences likely arise from methodological approximations, such as neglecting phonon dispersion over the entire Brillouin zone and omitting cross terms in the second-order coupling coefficients. While more rigorous treatments including these effects have been reported,^11,73–75 they require computationally intensive methods that are beyond the scope of the present study. Nevertheless, the current approach captures the dominant experimental trends and provides valuable qualitative insight into the differences between the MOF and molecular crystal systems. Having established the overall trends in spin–phonon coupling and relaxation behaviour, we will now discuss the origin of the longer spin–lattice relaxation times observed in the MOF systems.

Discussion – origin of longer spin–lattice relaxation times in MOFs based on phonon mode analysis

Having established a DFT framework for obtaining reliable phonon modes and SPC coefficients within the constraints of our approach (i.e., considering only the non-cross terms in the second-order SPC), we now turn to a comparative analysis of the MOF and their related molecular crystal systems. In the following discussion, we examine how variations in phonon modes and the magnitudes of their associated SPC coefficients contribute to the distinct spin–lattice relaxation behaviours observed experimentally.

Dominant phonon modes and their role in spin–lattice relaxation

The experimental Raman spectrum of (1) (Fig. S33 and Fig. 2) shows the first Raman-active peak (above the measurable cutoff of ∼10 cm⁻¹) at a lower frequency than that of the molecular crystal (3). Consistently, DFT calculations predict the lowest-frequency optical phonon mode at ∼18 cm⁻¹ for (1), compared to ∼30–40 cm⁻¹ for (3), with the exact values depending on the density functional and basis set employed (Fig. S34). Despite the longer spin–lattice relaxation times observed for the MOF (1), its vibrational spectrum does not shift to higher energies relative to (3). This indicates that the presence of lower-frequency phonon modes does not necessarily translate into faster spin–lattice relaxation. This finding contrasts with previous THz spectroscopy results on another vanadyl porphyrin MOF, where the absence of vibrational bands in the THz region was attributed to longer T₁ values at high temperatures.⁷ This discrepancy suggests that the correlation between low-frequency vibrational modes and spin–lattice relaxation is not universal but depends strongly on the coupling strength and character of the individual phonon modes, rather than their frequencies alone.

Animated visualization of the phonon modes in (1) shows that the lowest-frequency optical phonon mode at 18 cm⁻¹ corresponds to a delocalized vibration involving collective motion of the entire MOF framework (Fig. 6). In this mode, the porphyrin units undergo rigid displacements while preserving their local ligand environments, resulting in only small spin–phonon coupling (Table S13; SPC coefficients: first-order ∼10⁻¹⁴, second-order ∼10⁻⁹). Although this phonon mode has a relatively small SPC coefficient (Fig. 5), its low energy leads to significant thermal population even at low temperatures, making it the dominant contributor to the spin–lattice relaxation of (1) via second-order processes (Fig. S29). However, its contribution to the T₁ process remains much smaller than that of the most dominant phonon modes in (3) occurring at 208 cm⁻¹ and 209 cm⁻¹ due to the magnitude of their SPC couplings (on the order of 10⁻⁹ for (1) and 10⁻⁵ for (3) in g_x and g_y) (Table S15 and Fig. S29, S31). This difference in the size of spin–phonon coupling in (1) and (3) explains why the overall T₁ process in (1) is slower than in (3), even though the lowest-frequency mode dominates in (1). This demonstrates that vibrational frequency shifts alone do not dictate relaxation behaviour, but the coupling strength of specific phonon modes is equally critical.

Fig. 6 The lowest-frequency optical phonon mode at the Γ point in Cu MOF (1), occurring at 18 cm⁻¹, as obtained from periodic DFT calculations. The light blue arrows indicate the atomic displacement vectors corresponding to the phonon eigenvector. Visualized using V_Sim.

In the case of vanadyl species, the lowest experimental Raman bands in (2) and (4) are observed at similar frequencies (Fig. S33). For (2), the computed lowest-frequency optical phonon mode occurs at ∼20 cm⁻¹, comparable to that of analogous (1) (Fig. S35). While the frequency of the lowest-optical phonon mode in (4) shows dependence on the DFT method, appearing at 19 cm⁻¹ with PBE-D3/VTZP and 4 cm⁻¹ with PBE-D3 (VASP), this lowest frequency mode consistently shows the largest contribution to the second-order process in (4), regardless of the computational method used. By contrast, although (2) also exhibits major contributions from very low-frequency modes in a similar range (Fig. S30; modes in the range between 16 cm⁻¹ to 37 cm⁻¹), the magnitude of the spin–phonon coupling remains smaller than (4), resulting in weaker spin–lattice relaxation contributions. Therefore, as in the Cu(II) case, the key factor distinguishing the MOFs and molecular crystal systems is not the phonon frequencies themselves, but rather the difference in the magnitude of the spin–phonon coupling.

Correlation between Γ-point phonon-mode symmetry and spin–phonon coupling strength

To gain deeper insight into the differences in spin–phonon coupling between the MOF and molecular crystal systems, we investigated the correlation between the symmetry (irreducible representation) of the phonon modes and the magnitude of the spin–phonon coupling. The symmetry of a phonon mode determines how atomic displacements transform under the point group operations of the crystal and, consequently, which terms of the spin-Hamiltonian they can modulate. Previous studies have shown, based on group-theory analysis, that the totally symmetric vibrational modes are the major contributors to the first-order SPC.^16,69,70 The correlation of the symmetry of the phonon modes and the size of first-order SPC coefficients is shown in Fig. 7 for (1) and Fig. S36–S38 for others. The size of the SPC coefficients is plotted against the irreducible representation of the underlying phonon mode. For (2), the space group has been reduced to P1, in which all phonon modes belong to the same irreducible representation (A). Therefore, a symmetry-based comparison of coupling strength is not applicable for this system. Consistent with the previous symmetry-based analyses,^16,69,70 both (1) and (4) show the dominant contributions from totally symmetric A modes. For (3), totally symmetric A₁ modes tend to show major contributions along with A₂ modes. This highlights that symmetry considerations are crucial for suppressing first-order spin–lattice relaxation processes, as the previous studies show.^16,20,69,70

Fig. 7 Distributions of first-order and second-order spin–phonon coupling coefficients (g_x) (a) for Cu MOF (1) and (b) for CuTPP (3), categorized by the irreducible representation of the vibrational modes. Vibrational modes were computed using VASP with the PBE-D3 functional.

By contrast, for the second-order spin–phonon coupling coefficients, the correlations between specific distortion types and coupling magnitudes is less clear than as the first-order SPC (Fig. 7 and Fig. S36–S38). Nevertheless, phonon modes that do not satisfy the symmetry requirements for first-order coupling (non-totally symmetric modes) tend to exhibit larger second-order SPC values. This trend is reasonable, as the first-order coupling requires the g-factor to vary linearly as a function of Q while the second-order coupling requires it to vary quadratically.^69,83 While this may seem obvious, it highlights the need for caution for the common use of the first-order SPC coefficients to identify the vibrational modes that promote the spin–lattice relaxation based on the assumption that modes with large first-order SPC coefficients also have large second-order SPC coefficients. While there is a tendency for modes with small first-order SPC coefficients to exhibit larger second-order coefficients, this trend is less clear-cut than the prediction for first-order SPC from symmetry. For example, both A and B modes show similar sizes of second-order SPC in (1), as shown in Fig. 7a. As a result, predicting which vibrational modes play a significant role cannot be done as intuitively as for the first-order spin–phonon coupling, highlighting the importance of computational approaches in thoroughly understanding spin–lattice relaxation in addition to vibrational spectroscopy or group theory.

Structural distortions and their influence on spin–phonon coupling

The above calculations qualitatively reproduce the smaller T₁ values observed in (1) and (2). However, neither frequency shifts nor symmetry considerations alone fully explain why such MOFs exhibit smaller second-order SPC coefficients and consequently slower spin–lattice relaxation. Since spin–phonon coupling arises from phonon-induced perturbations to the ligand environment and the electronic structure, as reflected in changes to the g-tensor, the difference may stem from the magnitude of ligand displacements around the metal centre. In these type of MOFs, the porphyrin ligands are embedded in an extended framework through coordination to metal nodes. This may lead to the delocalization of phonon modes across the entire framework. Then, such delocalization likely reduces local distortions at the metal centre, thereby weakening the spin–phonon interaction.

To examine this hypothesis, we quantified the in-plane and out-of-plane distortions in (1)–(4). For in-plane distortions, we evaluated the average M–N (M = Cu, V) bond-length change (defined as the average distance variation between the metal ion and the four coordinating nitrogen atoms) for each phonon mode. For out-of-plane distortions, we performed normal-coordinate structural decomposition (NSD) analysis of the porphyrin units. This analysis decomposes structural distortions into five fundamental modes: saddling (A_2u), ruffling (B_2u), doming (B_1u), waving (E_g), and propelling (A_1u) in an ideal D_4h symmetric porphyrins (Fig. S39).^50,51 The total out-of-plane distortion, D_oop, can also be computed from this analysis. It should be noted that while these symmetry labels originate from the point group D_4h, the actual phonon modes of the MOFs and molecular crystals are symmetry-lowered combinations of these distortions. Thus, the symmetry label in the NSD analysis do not directly correspond to the irreducible representations of the phonon modes. Each distorted structure at Q = 1 was analysed with respect to each phonon mode to identify correlations between specific distortion types, magnitude of the out-of-distortions and spin–phonon coupling strengths.

Fig. 8 shows the average distance between the Cu and the nitrogen and the absolute change in out-of-plane deformation D_oop with respect to the equilibrium geometry (|ΔD_oop|) obtained from the NSD analysis for (1) and (3). The results for (2) and (4) are shown in Fig. S40. When comparing the MOFs and molecular crystals, specifically (1) vs. (3) and (2) vs. (4), no drastic difference is observed in the distribution of M–N distances (Table S12). For D_oop, (3) shows a smaller value than (1), while (4) shows a larger value than (2), showing no clear trend between the MOFs and molecular crystals (Table S12).

Fig. 8 The distribution of average Cu–N distance changes of phonon modes in (a) Cu MOF (1) (b) CuTPP (3) and out-of-plane distortion (|ΔD_oop|) of phonon modes in (c) (1) (d) (3).

Additionally, the correlation plots between the magnitudes of in-plane and out-of-plane distortions and their corresponding spin–phonon coupling coefficients (Fig. S41 and S42) are scattered, showing only a rather poor correlation for both the first-order and second-order SPC. This is in contrast with the generally accepted concept that, in order to achieve a longer T₁, designing a rigid ligand that suppresses the ligand vibrations around the metal centre is the key. This absence of predictable patterns suggests that the spin–phonon coupling arises from more complex mechanisms that cannot be predicted by only the magnitude of distortions of the ligand around the metal centre. As noted in the previous section, symmetry is also likely to play a role in determining which mode can have large first-order and second-order SPC.

Consistent with this hypothesis, for the first-order spin–phonon coupling coefficients, a slightly better correlation was observed between specific distortion types and coupling strength (Fig. S43–S46). The distortion mode that shows a positive correlation is system-dependent. For (3), modes with larger saddling and ruffling distortions tend to show larger SPC coefficients (Fig. S45). By contrast, in (1), both ruffling and waving modes appear relevant (Fig. S43). In (4), phonon modes involving pronounced doming distortions tend to exhibit larger coupling coefficients (Fig. S46), while no clear correlation between specific distortion types and spin–phonon coupling was observed in (2) (Fig. S44). Although the symmetry of a given distortion cannot be directly assigned to a single irreducible representation of a phonon mode since phonon modes are mixtures of multiple distortion components, this system-dependent behaviour is in line with the idea that the symmetry imposes constraints on the spin–phonon coupling, as discussed in the previous section. For the second-order SPC (Fig. S47–S50), the correlation between the specific distortion mode and the size of the SPC coefficient is barely seen. This is similar to the symmetry analysis in the previous section, where the correlation between the symmetry and the size of SPC was not as obvious as the first-order SPC. This further suggests the more complicated nature of the second-order SPC and the difficulty in intuitively predicting the types of vibrational modes that have larger second-order SPC.

The distributions of the average distance change of the metal and the nitrogen (M–N; M = Cu or V) and ΔD_oop as a function of phonon frequency are also shown in Fig. 9 for (1) and (3) and Fig. S51 for (2) and (4), and the overall patterns are quite similar across systems. In particular, the frequency-dependent profiles of M–N distance fluctuations are nearly identical between MOFs and molecular crystals. For instance, in the vanadyl systems (2) and (4), the phonon modes that perturb the M–N distance appear at nearly the same frequencies around 200, 280, 330, and 400 cm⁻¹. In the Cu(II) systems, although (1) and (3) show slightly less similarity than the vanadyl analogues, the main Cu–N perturbing modes are still clustered around the same frequencies: approximately at 200, 280, 330, and 400 cm⁻¹.

Fig. 9 The average Cu–N distance changes of phonon modes in (a) Cu MOF (1) (b) CuTPP (3) and out-of-plane distortion (ΔD_oop) of phonon modes in (c) (1) (d) (3) as a function of the phonon mode frequency.

Despite the similarity of vibrational modes, the spin–phonon coupling (SPC) coefficients can differ significantly. For example, in compound (3), the dominant contributing modes are found at 208 and 209 cm⁻¹ (PBE-D3, VASP; Fig. S31 and Table S15, B₁, B₂ mode respectively). A similar vibrational mode is present in (1) at a wavenumber of 197 cm⁻¹ based on the NSD analysis. While this mode in (1) makes a significant contribution to the first-order SPC (Fig. S29, first-order SPC coefficients ∼10⁻⁶), its second-order contribution is much smaller (second-order SPC coefficients: ∼10⁻¹⁰) (Table S13). This difference can be partially attributed to the distinct symmetry or space group of the two systems. In the MOF (1), consistent with the group-theory based analysis, the 197 cm⁻¹ mode is atotally symmetric A mode and induces a largely linear change in the g-factor as a function of Q, which effectively reduces the second-order coupling. Thus, the difference in spin–lattice relaxation behaviour between (1) and (3) can be partially explained by symmetry effects.

In the vanadyl systems, the second-order spin–phonon relaxation process in (4) is dominated by the lowest-frequency mode, located around 4 cm⁻¹ (VASP, PBE-D3; Fig. S32). This mode exhibits largely phenyl rotations but is also accompanied by large out-of-plane distortions, as shown in Fig. S52. By contrast, while the lowest-frequency mode in (2) at 16 cm⁻¹ also contributes significantly to the second-order process, primarily due to the temperature-dependent pre-factor, its actual SPC coefficient is much smaller (∼10⁻¹⁰) (Fig. S28 and Table S14). This difference is likely due to the significantly reduced ligand distortions around the vanadyl centre in (2), as illustrated in Fig. S51 and Table S14, S16. Therefore, in the case of vanadyl, the slower spin–lattice relaxation observed in the MOF (2) can be attributed to the suppression of low-frequency modes that strongly perturb the ligand environment. In (4), the molecular crystal supports a very low-frequency mode that is mainly characterized as phenyl rotation but also perturbs the ligand environment, enhancing relaxation. The MOF lattice modifies this mode and effectively suppresses the major relaxation pathway in (4). Although new low-frequency modes emerge in the MOF (2) due to the new framework structure, it perturbs the ligand environment much less, resulting in weaker spin–phonon coupling and longer T₁.

Overall, while vibrational modes with large out-of-plane or in-plane distortions can exhibit strong first- or second-order spin–phonon coupling, as exemplified by the first optical phonon mode in (4), not all such modes behave this way. Instead, symmetry appears to play a critical role, and the magnitude of SPC coefficients is likely governed by a complex interplay of symmetry, distortion amplitude, and other factors. This indicates the difficulty in intuitively identifying the vibrational mode from commonly employed simple vibrational spectroscopy to suppress in order to achieve longer T₁ and highlights the importance of computational approaches for accurately identifying the vibrational modes that strongly contribute to spin–lattice relaxation and for guiding the rational design of systems with longer relaxation times.

Conclusion

In this study, we employed density functional theory (DFT) and spin–phonon coupling (SPC) calculations to investigate the origin of slower spin–lattice relaxation times (T₁) in MOF-based copper(II) and vanadyl(IV) porphyrin systems compared to their molecular crystal analogues. Using Γ-point phonon modes from well-converged periodic calculations, the first-order and second-order (restricted to only non-cross terms) SPC coefficients were computed for each phonon mode. Based on these results, we analysed the differences between the MOF and molecular crystal systems in terms of vibrational frequency shifts, symmetry considerations, and the magnitude of ligand displacements around the metal centre.

Our analysis revealed several important insights. First of all, first-order SPC calculations, which have been used to analyse the local mode impact on T₁,^{13,16,20,39,41,69–71} did not reproduce the experimentally observed differences in T₁, emphasizing the need to consider second-order SPC coefficients. In particular, modes that are symmetry-forbidden for first-order coupling often exhibit significant second-order contributions. Secondly, our results show that the frequency alone does not explain T₁ differences. While symmetry seems to play a somewhat important role in determining the size of SPC both for the first-order and second-order, the correlation with the magnitude of the SPC coefficients and symmetry is not straightforward, especially for the second-order SPC. Predicting which vibrational modes dominate spin–lattice relaxation based solely on symmetry remains challenging. We also analysed the correlation between porphyrin ligand distortions and SPC coefficients. While modes involving large out-of-plane or in-plane distortions can yield strong SPC, not all do. The presence of such distortions alone does not guarantee strong coupling, and the contribution of each mode appears to depend on a complex interplay between symmetry, distortion amplitude, and additional factors that were not included in the present analysis.

Overall, our results demonstrate that no single descriptor, neither frequency, symmetry, nor distortion magnitude, can reliably predict the spin–phonon coupling strength. Particularly for second-order processes, these contributions seem to arise from complex mechanisms that cannot be captured through intuitive analysis alone. Therefore, while experimental characterization such as low-frequency vibrational spectroscopy is useful, computational approaches are also essential for identifying key vibrational modes responsible for spin–lattice relaxation and for guiding the rational design of ligands and frameworks that support longer T₁ times in molecular spin qubits.

Author contributions

Y. S. – conceptualization, data curation, formal analysis, funding acquisition, investigation, visualization, writing – original draft preparation, N. S. – software, formal analysis, writing – review & editing, R. S. – writing – review & editing, M. Y. – Supervision, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The datasets generated in this work and animated relevent vibrational modes are available at the TU Graz Repository: https://repository.tugraz.at/records/t33h0-dy041 (DOI: https://doi.org/10.3217/t33h0-dy041).

Supplementary information (SI) containing additional experimental details, figures, and analysis is available. See DOI: https://doi.org/10.1039/d5cp04453g.

Acknowledgements

This work was supported by the JSPS KAKENHI (no. 23KJ0169 to Y. S.) and JST-FOREST (JPMJFR203F to R. S.). Moreover, this research was funded in part by the Graz University of Technology through the Lead Project Porous Material@Work for Sustainability (LP-03). Y. S. is very grateful to Prof. Egbert Zojer from the Institute of Solid State Physics at Graz University of Technology, who significantly supported this work. Y. S. performed all calculations in his group during her research visit, which was financially supported by the COLABS program at Tohoku University and the JSPS KAKENHI. All calculations were performed using the Austrian Scientific Computing (ASC) infrastructure. Y. S. also thanks Prof. Masami Kanzaki (Institute for Planetary Materials, Okayama University) for granting us access to his low-frequency Raman spectroscopy setup and for his assistance during the measurements.

References

1. M. Atzori and R. Sessoli, J. Am. Chem. Soc., 2019, 141, 11339–11352.
2. S. von Kugelgen and D. E. Freedman, Science, 2019, 366, 1070–1071.
3. C.-J. Yu, S. von Kugelgen, D. W. Laorenza and D. E. Freedman, ACS Cent. Sci., 2021, 7, 712–723.
4. M. J. Graham, J. M. Zadrozny, M. S. Fataftah and D. E. Freedman, Chem. Mater., 2017, 29, 1885–1897.
5. M. Atzori, L. Tesi, E. Morra, M. Chiesa, L. Sorace and R. Sessoli, J. Am. Chem. Soc., 2016, 138, 2154–2157.
6. K. Bader, D. Dengler, S. Lenz, B. Endeward, S.-D. Jiang, P. Neugebauer and J. van Slageren, Nat. Commun., 2014, 5, 5304.
7. T. Yamabayashi, M. Atzori, L. Tesi, G. Cosquer, F. Santanni, M.-E. Boulon, E. Morra, S. Benci, R. Torre, M. Chiesa, L. Sorace, R. Sessoli and M. Yamashita, J. Am. Chem. Soc., 2018, 140, 12090–12101.
8. C.-J. Yu, S. von Kugelgen, M. D. Krzyaniak, W. Ji, W. R. Dichtel, M. R. Wasielewski and D. E. Freedman, Chem. Mater., 2020, 32, 10200–10206.
9. D. D. Traficante, Concepts Magn. Reson., 1991, 3, 171–177.
10. K. J. Standley and R. A. Vaughan, Electron Spin Relaxation Phenomena in Solids, 1969.
11. L. A. Mariano, V. H. A. Nguyen, J. B. Petersen, M. Björnsson, J. Bendix, G. R. Eaton, S. S. Eaton and A. Lunghi, Sci. Adv., 2025, 11, eadr0168.
12. A. J. Fielding, S. Fox, G. L. Millhauser, M. Chattopadhyay, P. M. H. Kroneck, G. Fritz, G. R. Eaton and S. S. Eaton, J. Magn. Reson., 2006, 179, 92–104.
13. M. J. Amdur, K. R. Mullin, M. J. Waters, D. Puggioni, M. K. Wojnar, M. Gu, L. Sun, P. H. Oyala, J. M. Rondinelli and D. E. Freedman, Chem. Sci., 2022, 13, 7034–7045.
14. P. G. Klemens, Phys. Rev., 1962, 125, 1795–1798.
15. M. R. Espinosa, F. Guerrero, N. P. Kazmierczak, P. H. Oyala, A. Hong, R. G. Hadt and T. Agapie, J. Am. Chem. Soc., 2025, 147, 14036–14042.
16. N. P. Kazmierczak, R. Mirzoyan and R. G. Hadt, J. Am. Chem. Soc., 2021, 143, 17305–17315.
17. L. Escalera-Moreno, N. Suaud, A. Gaita-Ariño and E. Coronado, J. Phys. Chem. Lett., 2017, 8, 1695–1700.
18. A. Urtizberea, E. Natividad, P. J. Alonso, M. A. Andrés, I. Gascón, M. Goldmann and O. Roubeau, Adv. Funct. Mater., 2018, 28, 1801695.
19. A. Urtizberea, E. Natividad, P. J. Alonso, L. Pérez-Martínez, M. A. Andrés, I. Gascón, I. Gimeno, F. Luis and O. Roubeau, Mater. Horiz., 2019, 7, 885–897.
20. N. P. Kazmierczak, N. E. Lopez, K. M. Luedecke and R. G. Hadt, Chem. Sci., 2024, 15, 2380–2390.
21. P. A. Banks, L. Burgess and M. T. Ruggiero, Phys. Chem. Chem. Phys., 2021, 23, 20038–20051.
22. C. F. Macrae, I. Sovago, S. J. Cottrell, P. T. A. Galek, P. McCabe, E. Pidcock, M. Platings, G. P. Shields, J. S. Stevens, M. Towler and P. A. Wood, J. Appl. Crystallogr., 2020, 53, 226–235.
23. M. Kanzaki, J. Miner. Pet. Sci., 2018, 113, 171219.
24. Y. Horii, M. Makino, T. Yamamoto, S. Tatsumi, H. Suzuki, M. Noguchi, T. Yoshida, T. Kajiwara, Z.-Y. Li and M. Yamashita, Inorg. Chem. Front., 2022, 9, 6271–6278.
25. S. Lebedkin, C. Blum, N. Stürzl, F. Hennrich and M. M. Kappes, Rev. Sci. Instrum., 2011, 82, 013705.
26. A. Erba, J. K. Desmarais, S. Casassa, B. Civalleri, L. Donà, I. J. Bush, B. Searle, L. Maschio, L. Edith-Daga, A. Cossard, C. Ribaldone, E. Ascrizzi, N. L. Marana, J.-P. Flament and B. Kirtman, J. Chem. Theory Comput., 2023, 19, 6891–6932.
27. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
28. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
29. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
30. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
31. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 59, 1758–1775.
32. A. Togo, L. Chaput, T. Tadano and I. Tanaka, J. Phys.: Condens. Matter, 2023, 35, 353001.
33. A. Togo, J. Phys. Soc. Jpn., 2023, 92, 012001.
34. M. G. B. Drew, P. C. H. Mitchell and C. E. Scott, Inorg. Chim. Acta, 1984, 82, 63–68.
35. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.
36. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2025, 15(2), e70019.
37. R. Dovesi, V. R. Saunders, C. Roetti, R. Orlando, C. M. Zicovich-Wilson, F. Pascale, B. Civalleri, K. Doll, N. M. Harrison, I. J. Bush, P. D’Arco, M. Llunell, M. Causà, Y. Noël, L. Maschio, A. Erba, M. Rerat, S. Casassa, B. G. Searle and J. K. Desmarais, CRYSTAL23 User's Manual, University of Torino.
38. A. Albino, S. Benci, L. Tesi, M. Atzori, R. Torre, S. Sanvito, R. Sessoli and A. Lunghi, Inorg. Chem., 2019, 58, 10260–10268.
39. P. C. Bruzzese, Y.-K. Liao, L. Donà, B. Civalleri, E. Salvadori and M. Chiesa, J. Phys. Chem. Lett., 2024, 15, 7161–7167.
40. F. Santanni, A. Albino, M. Atzori, D. Ranieri, E. Salvadori, M. Chiesa, A. Lunghi, A. Bencini, L. Sorace, F. Totti and R. Sessoli, Inorg. Chem., 2021, 60, 140–151.
41. A. Lunghi and S. Sanvito, Sci. Adv., 2019, 5, eaax7163.
42. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
43. W. M. Ames and S. C. Larsen, Phys. Chem. Chem. Phys., 2009, 11, 8266–8274.
44. S. K. Singh, M. Atanasov and F. Neese, J. Chem. Theory Comput., 2018, 14, 4662–4677.
45. E. van Lenthe, J. G. Snijders and E. J. Baerends, J. Chem. Phys., 1996, 105, 6505–6516.
46. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
47. F. Weigend, Phys. Chem. Chem. Phys., 2006, 8, 1057–1065.
48. W. Jentzen, J.-G. Ma and J. A. Shelnutt, Biophys. J., 1998, 74, 753–763.
49. W. Jentzen, X.-Z. Song and J. A. Shelnutt, J. Phys. Chem. B, 1997, 101, 1684–1699.
50. C. J. Kingsbury and M. O. Senge, Coord. Chem. Rev., 2021, 431, 213760.
51. J. Krumsieck and M. Bröring, Chem. – Eur. J., 2021, 27, 11580–11588.
52. S. A. Ajibade, L. Catalano, J. Kölbel, D. M. Mittleman and M. T. Ruggiero, J. Phys. Chem. Lett., 2024, 15, 5549–5555.
53. K. Hamilton and J. Neu, APL Mater., 2024, 12, 010903.
54. A. Erba, J. Baima, I. Bush, R. Orlando and R. Dovesi, J. Chem. Theory Comput., 2017, 13, 5019–5027.
55. L. Doná, J. G. Brandenburg and B. Civalleri, J. Chem. Phys., 2019, 151, 121101.
56. V. A. Rassolov, J. A. Pople, M. A. Ratner and T. L. Windus, J. Chem. Phys., 1998, 109, 1223–1229.
57. K. L. Schuchardt, B. T. Didier, T. Elsethagen, L. Sun, V. Gurumoorthi, J. Chase, J. Li and T. L. Windus, J. Chem. Inf. Model., 2007, 47, 1045–1052.
58. D. Feller, J. Comput. Chem., 1996, 17, 1571–1586.
59. B. P. Pritchard, D. Altarawy, B. Didier, T. D. Gibson and T. L. Windus, J. Chem. Inf. Model., 2019, 59, 4814–4820.
60. J. Laun, D. V. Oliveira and T. Bredow, J. Comput. Chem., 2018, 39, 1285–1290.
61. A. Schäfer, H. Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571–2577.
62. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
63. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
64. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
65. L. Legenstein, L. Reicht, T. Kamencek and E. Zojer, ACS Mater. Au, 2023, 3, 371–385.
66. M. Fratschko, N. Strasser, N. Taghizade, M. Linares-Moreau, J. C. Fischer, T. Zhao, I. A. Howard, P. Falcaro, E. Zojer and R. Resel, Cryst. Growth Des., 2025, 25, 3665–3679.
67. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
68. P. A. Banks, Z. Song and M. T. Ruggiero, J. Infrared, Millimeter, Terahertz Waves, 2020, 41, 1411–1429.
69. R. Mirzoyan, N. P. Kazmierczak and R. G. Hadt, Chem. – Eur. J., 2021, 27, 9482–9494.
70. R. Mirzoyan and R. G. Hadt, Phys. Chem. Chem. Phys., 2020, 22, 11249–11265.
71. N. P. Kazmierczak and R. G. Hadt, J. Am. Chem. Soc., 2022, 144, 20804–20814.
72. E. Garlatti, L. Tesi, A. Lunghi, M. Atzori, D. J. Voneshen, P. Santini, S. Sanvito, T. Guidi, R. Sessoli and S. Carretta, Nat. Commun., 2020, 11, 1751.
73. E. Garlatti, A. Albino, S. Chicco, V. H. A. Nguyen, F. Santanni, L. Paolasini, C. Mazzoli, R. Caciuffo, F. Totti, P. Santini, R. Sessoli, A. Lunghi and S. Carretta, Nat. Commun., 2023, 14, 1653.
74. A. Lunghi and S. Sanvito, J. Phys. Chem. Lett., 2020, 11, 6273–6278.
75. A. Lunghi, Sci. Adv., 2022, 8, eabn7880.
76. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
77. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
78. M. Douglas and N. M. Kroll, Ann. Phys., 1974, 82, 89–155.
79. B. A. Hess, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 32, 756–763.
80. A. Lunghi, Challenges Adv. Comput. Chem. Phys., 2023, 219–289.
81. S. S. Eaton, T. Yamabayashi, Y. Horii, M. Yamashita and G. R. Eaton, J. Am. Chem. Soc., 2025, 147, 13815–13823.
82. C.-Y. Huang, Phys. Rev., 1967, 154, 215–219.
83. N. P. Kazmierczak, R. Mirzoyan and R. G. Hadt, J. Am. Chem. Soc., 2021, 143, 17305–17315.

---