A machine-learning framework for accelerating spin-lattice relaxation simulations

The study of spin dynamics of coordination compounds is central to the most recent advancements in different fields such as information storage, sensing, information processing, spintronics and paramagnetic resonance for biosystems^1,2,3. Understanding how spin relaxes towards equilibrium in molecules represents a long-standing puzzle whose solution would lead to a better interpretation of experiments as well as powering computational design of new magnetic materials⁴.

Ab initio spin dynamics simulations based on open quantum systems theory, density functional theory (DFT) and multiconfigurational quantum chemistry techniques are nowadays considered a gold standard tool to perform spin relaxation simulations⁴. This methodology has now been successfully applied to a large number of magnetic systems, including spin-1/2^5,6,7, single-molecule magnets^8,9,10,11, solid-state defects^12,13, and surface-adsorbed atoms¹⁴. However, these simulations require the numerical calculation of phonons for hundreds of atoms and derivatives with multi-reference ab initio methods, making predictions extremely expensive⁴. Moreover, whilst full ab initio calculation of a few systems at the time is possible through modern-day high-performance computing resources, the discovery of new materials with long relaxation times will inevitably require the characterization of many compounds^15,16, making a brute force ab initio approach unsuitable.

Machine learning (ML) has gained increasing popularity in science over the last years as a tool to overcome known computational bottlenecks of ab initio methods or to extract hidden patterns in big data. Relevant to this work is its application in the generation of force fields (FF) from a dataset of ab initio data, referred to as a training set. A MLFF is specified by a parametric function that links atomic coordinates and the total energy of a system. The MLFF’s parameters are determined by minimization of the differences between the model’s predictions and the full quantum mechanical calculations. Once the parameters of the model are determined, i.e., the model is trained, the MLFF can rapidly output the energy and atomic forces for any molecular configuration, thus avoiding the expensive step of solving the Schrödinger equation. This field has extended considerably in the last years and a plethora of different MLFFs are now available^{17,18,19,20,21,22,23,24,25,26,27,28,29,30} and capable of reaching chemical accuracy in their predictions of several material properties^31,32,33,34. Importantly, MLFFs have been used to predict phonon spectra of a large variety of compounds such as molecular crystals³⁵, organic molecules³⁶ and inorganic solid compounds³⁷, and new schemes are currently under study to achieve high transferability^38,39.

While the accuracy that the state-of-the-art ML models can achieve has been largely shown⁴⁰, their applicability in real-world scenarios has been explored to a lesser extent. Only recently, studies have been conducted to investigate the possibility of performing long-scale molecular dynamics (MD) simulation with state-of-the-art ML models and to find key metrics beyond accuracy on benchmark sets^41,42. Arguably, the full potential of MLFFs in addressing practical chemistry or physics problems has yet to be fully unleashed and their successful integration to tackle them is still far from trivial.

We aim to fill this gap in the field of spin relaxation by presenting a general ML workflow for the efficient and robust geometrical optimization and prediction of harmonic molecular vibrations for gas-phase molecular systems, a necessary first step towards a future application in condensed matter. Our work builds on previous efforts to use equivariant ML models to predict spin Hamiltonian parameters^43,44 and molecular potential energy surfaces⁴⁵. The molecular potential energy surface is approximated as linear in the bispectrum components of the atomic environments and the training is performed with an on-the-fly active learning (AL) strategy requiring minimal user intervention⁴⁶. The method is applied to three compounds of interest for the molecular magnetism community, chosen for their challenging electronic structure and diversified complexity of the potential energy surface. We show that this scheme only requires up to 20% of the total number of calculations needed by a full ab initio simulation of spin relaxation. Excellent results on predictions of spin-phonon relaxation times are achieved by using these phonon spectra and only a minimal accuracy loss is observed when employing ML for both molecular vibrations and spin-phonon coupling. Finally, we show that the same framework applies to MD simulations, thus allowing us to open a window on spin relaxation beyond the conventional equilibrium harmonic lattice.

Spin-phonon relaxation simulations

Magnetic properties of single-ion coordination compounds can be described by mapping their electronic structure onto an effective Hamiltonian. For compounds with quenched orbital angular momentum, as most transition metal complexes, this effective Hamiltonian often takes the form

where (overrightarrow{{bf{S}}}) is the spin and D is the zero-field splitting tensor. In the case of compounds with un-quenched angular momentum, as for the lanthanide compounds, the form of the effective Hamiltonian is instead

where ({hat{O}}_{m}^{l}) are tesseral functions of the total angular momentum J. From now on, both Hamiltonians are indicated with ({hat{H}}_{0}).

To describe spin-phonon relaxation, we need to take into account the Orbach and Raman mechanisms, involving respectively one and two contemporary phonon exchanges between the spin and lattice. The relaxation rate ({hat{W}}_{ba}) between the ({hat{H}}_{0}) eigenstates (leftvert arightrangle) and (leftvert brightrangle) due to the Orbach process is given by ref. ⁸

where ℏω_ba = (E_b−E_a), E_a and E_b being the eigenvalues of ({hat{H}}_{0}), Q_α and ω_α are the phonon modes and frequencies respectively. G^1−ph is the Fourier transform of the correlation function of the phonon bath

where ({overline{n}}_{alpha }) is the Bose-Einstein distribution population of the phonon with frequency ω_α and δ is the Dirac-delta function. For computational purposes, the Dirac-delta is approximated as a normalized Gaussian function where the smearing, i.e., the standard deviation, can be changed by the user. In this work, the smearing value was selected based on the convergence of the relaxation time as advised in previous works⁴⁷.

For the Raman relaxation process, the transition rate reads¹¹

where

The function G^2−ph takes into account all the possible processes involving two phonons: double absorption, double emission and contemporary absorption and emission. For example, for the contemporary absorption of phonon of frequency ω_α and emission of a phonon ω_β, G^2−ph takes the following form:

It is important to note that in deriving Eq. (5), it is assumed that the spectrum of ({widehat{H}}_{0}) is non-degenerate. To fulfill this condition, simulations of Raman relaxation are performed in the presence of a small magnetic field (0.1 T) along the easy axis. It is worth mentioning that no magnetic fields are applied when performing ab initio calculations. A magnetic field is applied only during the spin-phonon relaxation simulations, in order to lift the degeneracy of the ground state Kramers doublets. Spin-phonon couplings, i.e., (partial {hat{H}}_{{rm{0}}}/partial {Q}_{alpha }), are calculated as

where i runs over all the Cartesian degrees of freedom, m_i are the masses of atoms, L_αiare the eigenvectors of the Hessian matrix and X_i are the Cartesian degrees of freedom. The quantities (partial {hat{H}}_{{rm{0}}}/partial {X}_{i}) are calculated by symmetric derivative with a displacement of 0.01 Å. A comprehensive derivation and discussion of the equations reported above is available in literature⁴.

Machine learning force fields

MLFFs are deployed to substitute the expensive ab initio geometrical optimization and calculation of molecular vibrations. In this work, we use the linear ML interatomic potential named Spectral Neighbor Analysis Potential (SNAP)⁴⁸. According to this model, the total energy of a system is written as the sum of single atomic energy contributions, further expanded in a linear combination of atomic descriptors called bispectrum components:

where ({{mathcal{B}}}_{k}(i)) is the k-th bispectrum component of atom i, while N_k and N_i are the number of bispectrum components in the expansion and the number of atoms in the system, respectively. The terms ({{mathcal{B}}}_{k}(i)) encode geometrical information related to the surrounding of atom i within an user-chosen cutoff radius R_cut. N_k is set by the user through the value of 2j_max⁴⁹. The interested reader can find the rigorous definition and mathematical details of bispectrum components in ref. ⁴⁹. The coefficients c_k(α_i) depend on the atom species identified by the index α_i and are determined by minimising a loss function with a regularization term λ (the exact mathematical expression of the function is reported in ref. ⁴⁶). Atomic forces are rapidly evaluated by analytical derivation of the total energy. To fully avail of the efficiency of MLFFs, we implement an AL scheme for the construction of training sets for linear models⁴⁶. The ultimate goal of AL is to build a training set for a ML model according to two different criteria: one to include a structure in the training set (i) and one to declare the AL over (ii). In the context of this paper, including a structure in the training set is a synonym of triggering an ab initio calculation to add the energy/forces data to the learning dataset. For (i), we choose the comparison between the uncertainty s on the predicted quantity and a threshold defined by the user through δ

where s_z is the error on the training set and δ is set by the user. δ is strictly larger or equal to 1 and the smaller the more frequently it will individuate new structures. The criterion (ii) instead is dependent on the particular task at hand. For example, in the context of AL during MD, it is standard to declare AL over if an user chosen number of MD steps is performed without including new structures. For other tasks, however, we implement different criteria that are thoroughly presented in the Results section.

We provide a complete overview of a full ML approach to predict spin-phonon relaxation by predicting also D and ({B}_{m}^{l}) with an equivariant extension of SNAP⁴⁴. We briefly present here the equation underpinning this model. Any physical property can be written as a combination of spherical tensor components ({T}_{m}^{l}) of order l and 2l + 1 components. According to this principle, D and ({B}_{m}^{l}) are decomposed in a combination of spherical tensors ({T}_{m}^{l}), furtherly written in terms of single atomic contributions ({T}_{m}^{l}(a))

where ({overline{Y}}_{m}^{l}(i)) is the spherical harmonic of the i th atomic environment and is defined as ({overline{Y}}_{m}^{l}(i)=mathop{sum }nolimits_{j}^{{N}_{j}}{Y}_{m}^{l}({hat{r}}_{ij})), where j runs over N_j neighboring atoms within R_cut from atom a. ({Y}_{m}^{l}) is the standard definition of a complex spherical harmonic and ({hat{r}}_{ij}) are the coordinates of the atom j rescaled by the coordinates of the atom i. The predictions of the spin Hamiltonian tensors are then used to evaluate (partial {hat{H}}_{{rm{0}}}/partial {X}_{i}) for the calculation of (partial {hat{H}}_{{rm{0}}}/partial {Q}_{alpha }) according to Eq. (8). By doing so, we avoid performing 6N ab initio calculations for the calculation of the symmetric derivatives. In order to build the training sets for tensorial prediction, we deploy an AL scheme where uncertainty is evaluated analogously to the case of forces in ref. ⁴⁶. For a tensor of order l, we simultaneously predict 2l + 1 components and evaluate a (2l + 1) × (2l + 1) matrix, whose diagonal elements represent the variances, i.e., uncertainties, of the predictions on the single tensorial components. The value of s to use in Eq. (10) is the square root of the largest value of its diagonal elements. The interested reader can find the specific expressions to evaluate the uncertainty on the predictions in ref. ⁴⁶.

It is possible to go beyond the harmonic approximation implied by Eq. (8) by extracting the spin-phonon couplings directly from MD simulations as shown in ref. ⁵⁰. For example, we expand at the linear order the variation of the zero field splitting tensor elements D_ij with respect to their mean values, i.e., δD_ij = D_ij − 〈D_ij〉, due to normal mode vibrations

where i and j label the Cartesian axis for tensor D. It is then possible to establish a proportionality between the half Fourier transform of the correlation functions of δD_ij and the spin-phonon couplings

Analogous expressions can be written for lanthanide complexes. Achieving an accurate evaluation of the left hand side (LHS) of Eq. (13) requires long MD simulations together with the evaluation of the tensor of interest. To address this resource-intensive process, we make use again of AL-generated ML FFs for the MD simulations and the equivariant extension of SNAP augmented with AL for fast prediction of the tensors over the MD trajectories.

Selection of molecules

Here we present the results related to the implementation of ML models for (i) the prediction of molecular vibrations and spin Hamiltonian tensors, and the profiles of spin-phonon relaxation time as a function of temperature, and (ii) the calculation of correlation functions of the spin Hamiltonian tensors over MD trajectories. Simulations are presented for three molecules whose spin relaxation profiles have been fully characterized both experimentally and with full quantum mechanical simulations⁸: [Co(C₃S₅)₂] (1)⁵¹,[CoL₂]⁵², where H₂L = 1,2-bis(methane-sulfonamido) (2) and [Dy(bbpen)Cl] (3)⁵³. The two Co²⁺ complexes exhibit a ground state with spin S = 3/2 described through Eq. (1), while the Dy³⁺ complex possesses a ground state with total angular momentum J = 15/2, described by Eq. (2). The chemical structures of 1–3 are reported in Fig. 1.

In (a) Co²⁺ compound 1, (b) Co²⁺ compound 2 and (c) the Dy³⁺ compound 3. Color code: cobalt in pink, dysprosium in cyan, oxygen in red, carbon in dark gray, sulphur in yellow, nitrogen in blue, hydrogen in white and chlorine in green.

Full size image

Machine learning of molecular vibrations

The AL strategy is here applied to the optimization of molecular structures and calculation of molecular vibrations. The goal of the protocol is to train a MLFF able to accurately predict the forces for the configurations used to calculate the molecular vibrations and frequencies. For this purpose, we designed the following workflow:

• geometrical optimization: the initial training set is the sole molecular crystallographic structure and its energy and atomic forces are computed with electronic structure methods. A geometrical optimization process is started with on-the-fly AL;

• MD simulation at 50 K: once the optimization is complete, the training set is reset. The only structure in the new training set is the one obtained from the final geometry of the ML-aided optimization. Also here an on-the-fly AL scheme is again employed during MD;

• molecular vibration prediction: using the ML force field trained through the above steps, vibrational spectra are predicted.

In particular, we note the importance of the second step, where an MD at 50 K is performed, to achieve good accuracy in the prediction of frequencies. This step in fact guarantees a sufficient sampling of the harmonic portion of the potential, which is not automatically performed during optimization.

The optimization routine is declared over when the following conditions are contemporarily fulfilled: (i) the RMS of the gradient vector is smaller then 0.01 kcal/mol/Å, (ii) the maximum value of the elements in the gradient is smaller than 0.1 kcal/mol/Å and (iii) the variation in energy between two consecutive minimization steps is smaller than 0.0001 kcal/mol. These conditions are chosen in order to closely resemble the “TightOpt” option as implemented in ORCA, therefore making it possible to measure the performance of the ML workflow compared to the full ab initio approach. Geometry optimization is performed by testing different values of δ to prove the robustness of the protocol. In the preliminary stage of this work, different values of δ have been tested in order to find a range that would balance accuracy with small training set sizes. Finally, we set δ equal to 5, 10, 20, and 30. With the ML models trained at the end of optimization, we have calculated the molecules’ vibrational modes and frequencies. To check whether the results for the root mean square error (RMSE) values on frequencies and hessian elements could be improved, we augment the training set by sampling more structures with AL during MD simulation at 50 K, thus guaranteeing a more effective sampling of structures around the energy minimum. Since all of the ML models trained after optimization lead to similar results, (see Tab. S1), only the optimized structure obtained for δ = 30 is considered for further exploration during MD at 50 K, because of the smallest final training set size. Parity plots for the comparison of the predictions between ML and the ab initio approach are shown in Figs. S1-S12. Different values of δ are also tested during MD at 50 K. In this case a range of plausible values of δ is suggested by a previous study⁴⁶, although it has been slightly modified for this work. The values of δ tested are 2.5, 2.75, 3, and 3.25. The AL-augmented MD stops when 100 ps are simulated without finding new configurations to include in the training set. As an example, we present here only the results for δ = 2.5 during MD in Fig. 2 and Tab. 1. The RMSE on vibrational frequencies and Hessian elements is always lower after MD at 50 K and equivalent results are obtained for all the values of δ (see Tab. S2 and Figs. S13-S24). Beyond accuracy, a key metric to judge the ML performance is the efficiency factor (EF), defined as [1 − (N_ML/N_DFT)] ⋅ 100, where N_ML and N_DFT are the number of ab initio calculations performed with the aid of ML and with a traditional full quantum mechanical method, respectively. N_ML is given by the sum of the final training set sizes for optimization and MD. For the three molecules, the efficiency factor is always around 80 % (see Tab. 1) and higher efficiency factors are obtained for other values of δ (see Tab. S2, clearly pointing out that ML, while preserving near-to-full ab initio accuracy, also leads to a consistent computational saving.

Comparison between molecular vibration frequencies calculated with ML and full ab initio approach. Results for (a) 1, (b) 2, and (c) 3. For clarity purposes, we report the bisector line in the graphs.

Full size image

Given a robust and accurate protocol for the prediction of vibrational spectra with ML, we proceed to test its integration within the framework of spin-phonon relaxation simulations to speed them up. Spin-phonon relaxation profiles show two different curves related to the Orbach and Raman mechanisms. In the Orbach mechanism, an initially fully polarized spin state M_S is excited to an intermediate spin state via a series of resonant phonon absorptions. Subsequently, a series of resonant phonon emissions leads the system towards the spin state—M_S. This mechanism leads to a relaxation time that exponentially increases in temperature as the phonon population grows according to the Bose-Einstein distribution. For the Raman mechanism, let us consider the case of contemporary absorption and emission of two phonons. In Raman relaxation, a fully polarized spin state M_S relaxes to the spin state—M_S mediated by virtual excited states. Low-energy phonons are usually found to drive this mechanism⁸. In general, therefore, the spin relaxation profiles show a regime at high temperatures where the Orbach process leads to relaxation and a Raman-dominated regime at low temperatures. Combining the goal of achieving a consistent acceleration of these simulations with the need for near-to-chemical accuracy for the quantities in Eqs. (3) and (5) represents a serious test-bed for any ML model. Results for δ = 2.5 during MD are reported in Fig. 3, where the discrepancy between the full quantum mechanical and ML approach is always within a factor of 10 for all the compounds, i.e., within the upper bound of the typical error that affects these simulations⁸.

Comparison between spin-phonon relaxation profiles as a function of temperature calculated with phonons from ML and with full ab initio approach. Results for (a) 1, (b) 2, and (c) 3. Solid and dotted lines refer to Orbach and Raman relaxation, respectively.

Full size image

Machine learning of spin-phonon coupling

The optimal results shown on spin relaxation predictions have motivated us to further predict the spin Hamiltonian tensors in order to realize a complete ML scheme. For this purpose, we use the training sets built during MD with δ = 2.5 as learning datasets to predict the spin Hamiltonian tensors with equivariant SNAP and evaluate the spin-phonon couplings. In Fig. 4 we report the training and test sets predictions for all compounds (see RMSE values in Tab. S3), where the test set is made of the 6N configurations used to evaluate the derivatives (partial {hat{H}}_{{rm{0}}}/partial {X}_{i}) appearing in Eq. (8). The plots relative to the fit of higher-order tensors of 3 are reported in Fig. S25.

Results for (a) 1, (b) 2 and (c) 3 (l = 2). The plots for l = 4, 6 are reported in SI. For clarity purposes, we report the bisector line in the graphs.

Full size image

The spin-phonon couplings appearing in Eq. (8) are evaluated using the predictions of equivariant SNAP and subsequently used to evaluate the Orbach and Raman relaxation times as a function of temperature. A brief comparison of the spin relaxation profiles in Fig. 5 obtained within the full ML approach and Fig. 3 shows that every compound is differently affected by using ML to obtain the spin-phonon couplings.

Comparison between spin-phonon relaxation profiles as a function of temperature calculated with a full ML and with full ab initio approach. Results for (a) 1, (b) 2, and (c) 3. Solid and dotted lines refer to Orbach and Raman relaxation respectively.

Full size image

While accuracy is maintained for 3, the discrepancy between ML and ab initio approach gets approximately 10 times worse for 1 and 2. The larger discrepancy is due to the limited accuracy of the ML model in predicting (partial {hat{H}}_{{rm{0}}}/partial {X}_{i}) with R² = 0.69 for 1, R² = 0.52 for 2 and R² = 0.32, 0.03, −7.5 × 10⁻³ for 3, order l = 2, 4, 6 respectively (see Fig. S26) and the spin-phonon density, defined as ({sum }_{alpha }{left(partial {D}_{ij}/partial {Q}_{alpha }right)}^{2}delta (omega -{omega }_{alpha })) averaged over the Cartesian degrees of freedom, i and j, for the Co²⁺ complexes and ({sum }_{alpha }{left(partial {B}_{m}^{l}/partial {Q}_{alpha }right)}^{2}delta (omega -{omega }_{alpha })) averaged over the allowed values of m for l = 2, 4, 6 for 3 (see Fig. S27). Considering the Co²⁺ complexes for example, error analysis shows that the relative error on the terms (partial {hat{H}}_{{rm{0}}}/partial {X}_{i}) is proportional to σ_D/ΔD, where σ_D is the error in the prediction of D and ΔD the variation of the tensor component at the displaced geometries when evaluating the symmetric derivatives. Going to larger values of the derivative step in order to increase ΔD may improve the fit, but it could also invalidate the linear approximation underlying the derivative calculation. Manual tuning of the model’s hyperparameters has not led to significant improvements.

Machine learning the full potential energy surface and its spin Hamiltonian

Spin-phonon couplings can be extracted through analysis of MD simulations according to Eq. (13). However, evaluation of the LHS of Eq. (13) requires long MD simulations in order to accurately calculate the averages, e.g., 〈δD_ij(0)δD_ij(t)〉. To tackle this challenge, we use the following workflow, which combines ML models for force fields and spin Hamiltonian tensors together with AL:

• to speed up the construction of the learning set, we decide to generate the ML model with the set built for molecular vibrations with δ = 2.5. This is not a strict requirement and starting from scratch is possible;

• new configurations are added to the training set during an AL-augmented MD performed at 300 K (for Co²⁺ complexes) and 100 K (for Dy³⁺ complex) with δ set to 1.75;

• using the ML force field trained through the above steps, MD trajectories approximately 40–60 ns long at 25, 50, 75 and 100 K are generated. Correlation functions for the spin Hamiltonian tensors are evaluated on these MD simulations.

Analysis of the correlation functions is done only up to 100 K because experimental spin relaxation data is commonly available only at low temperatures. AL-augmented MD is performed at a lower temperature for the Dy³⁺ compound as the FF becomes substantially less accurate in its predictions if trained at 300 K, possibly due to the limited stability of this molecule at high-temperature and gas phase. The AL is over when no new structures are included after 0.5 ns of MD. To assess the accuracy of the ML models, we benchmark them on test sets drawn from MD simulations at all 4 different temperatures. The test set is made of 50 structures sampled evenly during 2 ns of dynamics. The generated FFs achieve chemical accuracy on test sets for 1 and 3 with an upper bound on the RMSE for energies (forces) of 0.30 (0.99) and 0.57 (2.24) respectively (energies in kcal/mol and forces in kcal/mol/Å). For 2, the RMSE upper bounds on energies (forces) of 2.38 (1.68) kcal/mol (kcal/mol/Å) are instead slightly above chemical accuracy. For results at all temperatures for all compounds and on the training set, see Tabs. S4 and S5 and Figs. S28-S39. To evaluate the correlation functions of the spin Hamiltonian tensors over the MD trajectories, we use the following protocol:

• with no a priori training dataset we perform an AL (δ = 1.5) targeting D and ({B}_{m}^{l}) over 2 ns of trajectory;

• with the ML models trained through the above step, we predict the spin Hamiltonian tensors over approximately 10M structures for every compound at each temperature, hence their correlation functions along the MD trajectories.

Also in this case, we validate the accuracy of the ML model by calculating the RMSE of its predictions for the same configurations in the test set used for MD. RMSE values on test sets for spin Hamiltonian tensors for 1 range between 1.49 and 2.90 cm⁻¹, and 0.97 and 1.90 cm⁻¹ for 2. For 3, RMSE test set values depend on the order l: it ranges between 0.31 and 0.74 cm⁻¹ (l = 2), 2.7 × 10⁻³ and 5.9 × 10⁻³ cm⁻¹ (l = 4) and 5 × 10⁻⁵ and 1.2 × 10⁻⁴ cm⁻¹ (l = 6). For more details on RMSE values on the training sets and for specific RMSE values on the test set, see Tabs. S4 and S5, respectively. Parity plots related to the results on tensorial predictions are shown in Figs. S28–S39.

Figure 6, reports the plot for element D₁₁ of compound 1 at four different temperatures. As expected, the heights of the peaks directly correlate with the temperature since the proportionality factor omitted in Eq. (13) is the thermal population of phonons. As can be seen from these results, the simulation of spin-phonon coupling from MD automatically accounts for the presence and the thermal evolution of the phonons’ frequencies and linewidths. Figures 6 and S40–S43 show that such features cannot be reproduced with a simple equilibrium harmonic approximation of the lattice dynamics.

a Average of the Fourier transforms of the correlation functions of δD_ij at four different temperatures. b Average of the Fourier transforms of the correlation functions of δD_ij calculated as in LHS of Eq. (13) at T = 25 K.

Full size image

The correlation functions could in principle be used to estimate the relaxation time. However, a proper molecular dynamics spin relaxation theory is not available yet, and its derivation is beyond the scope of this work. Indeed, whilst MD is able to include non-equilibrium effects and phonons’ dissipation, it is based on a classical notion of nuclei, which is expected to be of limited accuracy at such low temperatures. Moreover, MD has only been adapted to second-order perturbation theory so far, while Raman relaxation will require the development of a fourth-order theory.

In this study, we have presented a ML workflow based on AL that leads to computational savings of at least 80% in the prediction of molecular vibrations and spin-phonon couplings keeping near-to-full ab initio accuracy. The precision achieved on vibrational spectra has been proven to be sufficient to reach excellent results also in the prediction of spin-phonon relaxation. Due to the non-trivial dependence of the spin-phonon relaxation rate on molecular vibrations, this has to be considered a remarkable result.

We have also tested an equivariant extension of SNAP to predict the Spin Hamiltonian tensors proving its successful integration within the spin relaxation framework. We notice a worse agreement between ML and ab initio calculation for the Co-based complexes compared to the Dy-based one, which appears as non-trivial and warrants further studies. This aspect of the work has room for refinement and we believe that different and complete chemical descriptors, improvable ad libitum^27,54, represent a viable solution to this problem. However, the higher number of parameters associated with these descriptors is usually linked to larger training sets, sometimes by orders of magnitude, therefore possibly reducing the computational savings offered by the presented workflow. A detailed benchmark of different descriptors and their trade-off between efficiency and accuracy is a necessary next step.

Futhermore, in this work, we have not investigated the strategies to tune the dataset in order to better describe the vibrations more relevant for the relaxation mechanism, e.g., by better fitting forces on specific atoms. We believe that this should be object of further studies in order to achieve higher accuracies.

We note that vibrational spectra could have also been obtained through analytical differentiation. However, in general, the memory requirements of the analytical approach scale unfavorably with the system size and are not necessarily more efficient than the numerical differentiation⁵⁵. In view of the applications of this method to large molecules and crystals, the favorable scaling of energy and gradients is key. Moreover, analytical frequencies are not necessarily always available for any level of theory in any software. The present workflow on the other hand can be adapted to any electronic structure method used to produce the ab initio reference data. Similarly, the analytical calculation of spin-phonon coupling coefficients has been recently proposed⁵⁶, but it lacks the inclusion of spin-orbit coupling derivatives, and a numerical differentiation remains advisable for the time being.

Aiming to go beyond the equilibrium harmonic approximation for molecular vibrations, we have shown how SNAP and its equivariant version can be deployed to study the correlation functions of the spin Hamiltonian tensors, a task demanding tens of millions of tensorial predictions and generation of MD simulations on the scale of 50 ns. While this is far beyond the reach of any ab initio approach, regardless of the availability of analytical Hessians and gradients, it has required the execution of only hundreds of first-principles calculations. In⁵⁰, only the use of ML techniques has made possible the aforementioned calculations for an organic radical in solution. Here we have made a step forward with respect to⁵⁰ by challenging the capabilities of the ML model on the dynamics of more complex systems and by deploying an AL scheme for the construction of the training dataset to generate MD simulations. These kinds of simulations represent a challenging testbed for the AL framework since even a single inaccurate prediction during MD can lead to catastrophic results. Furthermore, equivariant SNAP is combined with the AL scheme in ref. ⁴⁶ to achieve a rational construction of the training set also for the prediction of D and ({B}_{m}^{l}).

This is a proof of the possibility to successfully deploy ML models to study magneto-structural correlations, which are of paramount importance to design molecular systems with tailored properties, e.g., design of potential molecular qubits^57,58.

In conclusion, we have presented here a ML method for the prediction of vibrational spectra and spin-phonon couplings for molecules in the gas phase with computational savings of at least 80% compared to the full ab initio approach. We anticipate that this work, along with its future extensions to periodic systems, will significantly expand the scope of numerical simulations for the study of spin relaxation in condensed matter and non-equilibrium settings.

Electronic structure calculations

To calculate the magnetic properties of the compounds we use the software ORCA 4.2.1⁵⁹. Firstly, a multireference wavefunction is obtained at the level of Complete Active Space Self Consistent Field (CASSCF) theory. Spin-orbit coupling is included at the level of quasi-degenerate perturbation theory to calculate the zero-field splitting tensor D and the coefficients ({B}_{m}^{l}). For Dy³⁺ ions, we employ an active space of 7 4f orbitals with 9 electrons using all the solutions with multiplicity 6. For Co²⁺ we consider 7 electrons in 5 3d orbitals (7,5) using 10 solutions with multiplicity 4 and 40 solutions with multiplicity 2. Relativistic effects of the CASSCF calculations are treated by the Douglas-Kroll-Hess method⁶⁰.The RIJCOSX approximation for coulomb and exchange integrals is used for complexes containing both ions. DKH-def2-TZVPP basis set is used for all elements with an “AutoAux” generated auxiliary set, except for Dy for which SARC-DKH-TZVPP is used as a basis set.

To evaluate the expressions in Eqs. (3)–(8), it is necessary to calculate the phonon spectrum. For this purpose, firstly we use ORCA to optimize the geometrical configuration at hand at the level of DFT. Tight optimization settings, exchange-correlation functional BP86 with D3BJ⁶¹ corrections along with def2-TZVPP and def2/J as basis and auxilary set are used. Large integration grid settings are deployed (Grid 7, GridPruning 0), with increased precision for cobalt (SpecialGridIntAcc 10) and very tight SCF convergence criteria. For the dysprosium compounds, DFT calculations are carried out replacing Dy with Y. The diamagnetic substitution Dy-Y is performed to allow for faster convergence of DFT calculations without leading to major structural differences with respect to the original magnetic molecule. It should be noted that this approximation should be carefully checked when comparing simulations and experiments, but bears no prejudice to the quality of the present investigation as the comparison is strictly among simulated relaxation times, both obtained under the same approximation. Once the optimized structure is obtained, 6N DFT single point calculations with the same settings used for optimization are performed for configurations where atoms are individually displaced along the Cartesian axes by 0.001 Å. Atomic forces are extracted from these calculations and provided as input to MolForge to determine the molecular vibrational spectra.

Machine learning force fields

All the training and AL for the scalar version of SNAP are performed targeting both energies and forces. With reference to the notation in ref. ⁴⁶, all the hyperparameters are kept fixed across the different compounds (2j_max = 12, λ = 0.1, R_cut = 4.0 Å). Only the relative weight for the training of energy and forces depends on the specific compound and is set to 3N, where N is the number of atoms in the system in order to equally weight energies and forces during the training. Optimizations performed with MLFFs employ the Adam algorithm⁶². Using the notation of⁶², the learning rate α = 10⁻⁴, the exponential decay rates for the first and second moment estimates β₁ = 0.9 and β₂ = 0.999 and ϵ = 10⁻⁷, a “small” number that avoids division by 0 in the execution of the algorithm.

For equivariant SNAP, the hyperparameters are set to 2j_max = 8, λ = 0.1, and R_cut = 4.0 Å for all compounds.

MD trajectories for the study of the correlation functions are produced with Large-scale Atomic/Molecular Massively Parallel Simulator(LAMMPS)⁶³, setting the timestep to 1 fs and geometrical configurations are dumped every 5 fs.