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

Introduction
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 biosystems1,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 materials4.
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 simulations4. This methodology has now been successfully applied to a large number of magnetic systems, including spin-1/25,6,7, single-molecule magnets8,9,10,11, solid-state defects12,13, and surface-adsorbed atoms14. However, these simulations require the numerical calculation of phonons for hundreds of atoms and derivatives with multi-reference ab initio methods, making predictions extremely expensive4. 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 compounds15,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 available17,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 properties31,32,33,34. Importantly, MLFFs have been used to predict phonon spectra of a large variety of compounds such as molecular crystals35, organic molecules36 and inorganic solid compounds37, and new schemes are currently under study to achieve high transferability38,39.
While the accuracy that the state-of-the-art ML models can achieve has been largely shown40, 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 sets41,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 parameters43,44 and molecular potential energy surfaces45. 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 intervention46. 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.
Results
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. 8
where ℏωba = (Eb−Ea), Ea and Eb being the eigenvalues of ({hat{H}}_{0}), Qα and ωα are the phonon modes and frequencies respectively. G1−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 works47.
For the Raman relaxation process, the transition rate reads11
where
The function G2−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 ωβ, G2−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, mi are the masses of atoms, Lαiare the eigenvectors of the Hessian matrix and Xi 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 literature4.
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)48. 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 Nk and Ni 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 Rcut. Nk is set by the user through the value of 2jmax49. The interested reader can find the rigorous definition and mathematical details of bispectrum components in ref. 49. The coefficients ck(α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. 46). 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 models46. 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 sz 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 SNAP44. 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 Nj neighboring atoms within Rcut 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. 46. 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. 46.
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. 50. For example, we expand at the linear order the variation of the zero field splitting tensor elements Dij with respect to their mean values, i.e., δDij = Dij − 〈Dij〉, 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 δDij 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 simulations8: [Co(C3S5)2] (1)51,[CoL2]52, where H2L = 1,2-bis(methane-sulfonamido) (2) and [Dy(bbpen)Cl] (3)53. The two Co2+ complexes exhibit a ground state with spin S = 3/2 described through Eq. (1), while the Dy3+ 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) Co2+ compound 1, (b) Co2+ compound 2 and (c) the Dy3+ 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.
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 study46, 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 − (NML/NDFT)] ⋅ 100, where NML and NDFT are the number of ab initio calculations performed with the aid of ML and with a traditional full quantum mechanical method, respectively. NML 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.
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 MS 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—MS. 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 MS relaxes to the spin state—MS mediated by virtual excited states. Low-energy phonons are usually found to drive this mechanism8. 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 simulations8.

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.
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.
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.
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 R2 = 0.69 for 1, R2 = 0.52 for 2 and R2 = 0.32, 0.03, −7.5 × 10−3 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 Co2+ 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 Co2+ 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., 〈δDij(0)δDij(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 Co2+ complexes) and 100 K (for Dy3+ 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 Dy3+ 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−1, and 0.97 and 1.90 cm−1 for 2. For 3, RMSE test set values depend on the order l: it ranges between 0.31 and 0.74 cm−1 (l = 2), 2.7 × 10−3 and 5.9 × 10−3 cm−1 (l = 4) and 5 × 10−5 and 1.2 × 10−4 cm−1 (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 D11 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 δDij at four different temperatures. b Average of the Fourier transforms of the correlation functions of δDij calculated as in LHS of Eq. (13) at T = 25 K.
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.
Discussion
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 libitum27,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 differentiation55. 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 proposed56, 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. In50, 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 to50 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. 46 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 qubits57,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.
Methods
Electronic structure calculations
To calculate the magnetic properties of the compounds we use the software ORCA 4.2.159. 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 Dy3+ ions, we employ an active space of 7 4f orbitals with 9 electrons using all the solutions with multiplicity 6. For Co2+ 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 method60.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 D3BJ61 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. 46, all the hyperparameters are kept fixed across the different compounds (2jmax = 12, λ = 0.1, Rcut = 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 algorithm62. Using the notation of62, the learning rate α = 10−4, the exponential decay rates for the first and second moment estimates β1 = 0.9 and β2 = 0.999 and ϵ = 10−7, a “small” number that avoids division by 0 in the execution of the algorithm.
For equivariant SNAP, the hyperparameters are set to 2jmax = 8, λ = 0.1, and Rcut = 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)63, setting the timestep to 1 fs and geometrical configurations are dumped every 5 fs.
Responses