The simulation of the second-order derivatives poses serious challenges from a computational point of view, as the number of required simulations scales quadratically with the number of atoms and the mesh used for the polynomial interpolation. In Ref. [1]}, a solution to this problem based on machine learning was brought forward. The idea is to first learn a general function \(f_{ML}(X_{is},\alpha _{i})\)  that outputs the spin Hamiltonian parameters as function of molecular coordinates, \(X_{is}\) , and then use it to predict all the data necessary to numerically estimate the second-order derivatives of Eq. REF . The learning of such a function is done by optimizing the built-in coefficients \(\alpha _{i}\)  in order to reproduce the values of spin Hamiltonian of a set of reference calculations[2]}, [3]}. It should be noted that the method provides an advantage only if the number of calculations needed to compute the training set of \(f_{ML}(X_{is},\alpha _{i})\)  is smaller than the number of calculations required for a brute-force finite difference calculation. There is in fact no conceptual difference from using a polynomial function to perform a regression of values of spin Hamiltonian as function of molecular distortions or the use of machine learning. The crucial difference lies in the way the two approaches accounts for non-linearity and their reciprocal learning rate. In the case of polynomial functions this is handled by increasing the polynomial order, which rapidly leads to an increase of coefficients that needs to be determined. Machine learning offers a few advantages: i) symmetries are automatically included in the formalisms, ii) molecular structure can be described in terms of very compact descriptors, such as bispectrum components[4]}, which reduces the number of independent variables that needs to be determined, and iii) the function is represented by a (shallow) neural network able to optimally fit any function. Two flavours of machine learning methods were recently adopted to predict the spin-phonon coupling coefficients of magnetic molecules, both in the context of solid-state perturbation theory[1]} and molecular dynamics of molecules in liquid solution[6]}. A schematic representation of the former method, based on Ridge regression and bispectrum components as molecular fingerprints is reported in Fig. REF . In a nutshell, a supercell approach to the calculation of phonons is used in combination with DFT. The optimized molecular structure is then distorted in the range of 1000 times by applying small random perturbations of the atomic positions. The spin Hamiltonian coefficients of each distorted structure are determined with electronic structure methods and form the basis for the generation of a ML model. A part of this set is usually not used for training the ML model but just for testing purposes. The ML model used in Ref. [1]} predicted the spin Hamiltonian coefficients by means of the knowledge of each atom's coordination environment, as graphically sketched in Fig. REF B. After a successful training of the ML, it is possible to very quickly scan the spin Hamiltonian coefficients along any pair of Cartesian directions[3]}. Fig REF D reports just an example out the many 2D plots needed to computed all the derivatives needed in Eq. REF .
