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

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 \(\left\vert a\right\rangle\) and \(\left\vert b\right\rangle\) 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.

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.

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⁸.

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.

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.

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.

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.