Constrained DFT-based magnetic machine-learning potentials for magnetic alloys: a case study of Fe–Al

We propose a machine-learning interatomic potential for multi-component magnetic materials. In this potential we consider magnetic moments as degrees of freedom (features) along with atomic positions, atomic types, and lattice vectors. We create a training set with constrained DFT (cDFT) that allows us to calculate energies of configurations with non-equilibrium (excited) magnetic moments and, thus, it is possible to construct the training set in a wide configuration space with great variety of non-equilibrium atomic positions, magnetic moments, and lattice vectors. Such a training set makes possible to fit reliable potentials that will allow us to predict properties of configurations in the excited states (including the ones with non-equilibrium magnetic moments). We verify the trained potentials on the system of bcc Fe–Al with different concentrations of Al and Fe and different ways Al and Fe atoms occupy the supercell sites. Here, we show that the formation energies, the equilibrium lattice parameters, and the total magnetic moments of the unit cell for different Fe–Al structures calculated with machine-learning potentials are in good correspondence with the ones obtained with DFT. We also demonstrate that the theoretical calculations conducted in this study qualitatively reproduce the experimentally-observed anomalous volume-composition dependence in the Fe–Al system.


www.nature.com/scientificreports/
To overcome this problem constraints are introduced in the Kohn-Sham equations, either as soft constraints contributing a penalty function to the total energy [16][17][18] or, as recently proposed, hard constraints solved selfconsistently using the Lagrangian multiplier method 19 .
Despite the success of different flavours of DFT calculations in investigating different materials during the last 30 years they are still computationally expensive even when used them on modern supercomputers.Moreover, in computational studies of magnetic materials, including magnetic moments as an additional degree of freedom in the calculation of DFT energy, the computational time increases drastically.To overcome the limitations of DFT in size and time, requires alternative approaches.Effective interaction model (EIM) was proposed in Refs. 20,21nd applied to the Fe-Ni system.One approach that has gained massive interest in the computational materials community over the last 15 years are machine-learning interatomic potentials (MLIPs) [22][23][24][25][26][27][28][29][30][31][32] .The main idea behind MLIPs is their ability to avoid running the full DFT simulation by interpolating between a relatively small training set of carefully selected single-point DFT calculations.A MLIP that has been trained on configurations that are representative for the entire configurational space appearing in a simulation then approximates local DFT energies and forces with, in principle, arbitrary accuracy-contrary to (semi-)empirical interatomic potentials.3][24][25][26][27][28][29][30][31][32] ) have been developed for non-or ferromagnetic materials, which do not require taking into account magnetic degrees of freedom explicitly in the functional of MLIPs.In order to approximate DFT energies that strongly depend on the magnitude and direction of the magnetic field we must enrich the functional form of MLIPs with magnetic moment features which has been proposed in a number of works [33][34][35][36][37][38] .The open problem, however, is the generation of suitable training sets.We emphasize that in addition to the configurations with non-equilibrium atomic positions and lattice parameters we also need configurations with non-equilibrium magnetic moments for the proper fitting of magnetic MLIP.We discuss the creation of such a training set in the present work.
In this paper we generalize a single-component magnetic Moment Tensor Potential 35 (one of the recently developed MLIPs with magnetic degrees of freedom, mMTP) to the case of multi-component magnetic materials.
Here, we use the recently developed cDFT approach 19 to evaluate the magnetism of Fe-Al system by first-principles calculations being such approach seen to efficiently describe magnetism in considering a potential-based formulation of the self-consistency problem 19 .These cDFT data will be subsequently used to generate a training set for fitting an mMTP, i.e. in addition to the configurations with equilibrium magnetic moments we also include the ones with non-equilibrium (or, perturbed) magnetic moments in the training set.An mMTP fitted to such a training set allows for equilibrating magnetic moments of a configuration starting from the perturbed ones.In other words, we consider magnetic moments on the same grounds as atomic positions: we generate a training set with both relaxed and perturbed atomic positions, fit a potential, and, finally, relax both (magnetic and "positional") degrees of freedom of the structures with the trained potential.
We test our multi-component mMTPs on the system of Fe-Al with different concentrations and positions of Al and Fe.We demonstrate that the formation energies, the equilibrium lattice parameters, and the total magnetic moments of the unit cell predicted with the fitted mMTPs for various compounds of the Fe-Al alloy are in a good correspondence with the DFT ones.We also demonstrate that the theoretical calculations conducted in this paper reproduce the anomalous volume-composition dependence in the Fe-Al system shown in Ref. 13 and experimentally observed 39 .

Training set
The training set for magnetic Moment Tensor Potential (mMTP) fitting consists of different configurations of bcc Fe-Al in which the concentration of Al varies from 0 to 50% .We additionally include pure bcc-Al in the training set.The configurations in the training set consist of 16 atoms ( 2 × 2 × 2 conventional bcc cell).There are altogether 2012 configurations in the training set constructed as detailed below.
For constructing the training set we start from the 23 "parent" configurations taken from Ref. 13 .These configurations differ the way Al occupies the supercell sites as given in Table 1.The notations for the supercell sites (as used in the first column of Table 1 where they denote the locations of the Al atoms) are shown in Figure 1.For generating the configurations for the training set we first conduct full geometrical optimization (relaxation) of the structures (i.e., we minimize energies with respect to atomic positions, lattice vectors, and magnetic moments) from the left column in Table 1.We took the Fe magnetic moments of 3 µ B and the Al magnetic moments of 0 µ B as an initial guess.After the minimization we added the configurations from the relaxation path to the training set (see Step 1 in Table 1).After that we compress and extend lattice vectors of equilibrium configurations by 1 % .Next, we take each of these configurations, and we apply random displacements to each coordinate of each atomic position in the range from −0.1 to 0.1 Å; we do it five times and obtain five configura- tions from each equilibrium, compressed and extended configuration.For the resulting displaced configurations we optimize magnetic moments (see Step 2 in Table 1).In Step 3, as opposed to Steps 1 and 2 in which we conduct DFT calculations without any constraints, we impose hard constraints on magnetic moments to obtain the configurations with non-equilibrium magnetic moments.In Step 3.1 we take each configuration from the previous step and five times randomly perturb each equilibrium atomic magnetic moment by increasing or decreasing this value by at most 15% , and calculate these configurations with cDFT (see Step 3.1 in Table 1).As a result, the maximum deviation from the equilibrium magnetic moment for the Fe atom could reach 0.4 µ B , for the Al atom the maximum deviation could be 0.01 µ B (in the configurations with both Fe and Al atom).In the case of pure Al we randomly perturb magnetic moments by the values from the range (−0.03, 0.03) µ B .Finally, we randomly perturb equilibrium magnetic moments in the configurations obtained in Step 1 by increasing or decreasing their value by at most 50% and also conduct cDFT calculations for the configurations with non- equilibrium magnetic moments, but with the equilibrium positions and lattice vectors (see Step 3.2 in Table 1).As a result, the maximum deviation from the equilibrium magnetic moment for the Fe atom could reach 1.3 µ B , for the Al atom the maximum deviation could be 0.04 µ B (in the configurations with both Fe and Al atom).In the configurations with pure Al we randomly perturb magnetic moments by the values from the range (−0.1, 0.1) µ B .Thus, we obtain the training set including fully equilibrium configurations and the ones with perturbed atomic positions, lattice vectors, and magnetic moments.Number of configurations converged for each step and each structure is given in Table 1.The training set contains the total of 2012 configurations.We note that around 80 % of Fe-Al configurations converged during the DFT (cDFT) calculations with the ABINIT code.We also emphasize that all the cDFT calculations with perturbed magnetic moments and equilibrium atomic positions and lattice parameters (at Step 3.2) were converged.
Table 1.Construction of the training set.We first conduct full geometrical optimization (relaxation) of structures and get fully equilibrium structures and the ones from the relaxation path (Step 1).Next we displace atomic positions and lattice vectors and run DFT calculations with the optimization of magnetic moments (Step 2).Finally, we perturb magnetic moments in the configurations obtained after Steps 1 and 2 and get the configurations with non-equilibrium magnetic moments (Steps 3.1 and 3.2).The numbers of converged configurations for each step are given in the table.The total number of converged configurations is 2012.Number of configurations in the training set: 2012.

Structure
Step For the verification of the fitted mMTPs we also generate a "verification" set.We start with the configurations obtained in Step 1 of constructing the training set, and perform Steps 2 and 3.2 (i.e., we perturb atomic positions and lattice parameters and optimize magnetic moments (Step 2) and we perturb magnetic moments at equilibrium lattice parameters and atomic positions (Step 3.2)).We thus obtain additional 336 configurations for the verification of the fitted mMTPs.Unlike for the training set, we omitted Step 3.1 for the purpose of generating the verification set (and we therefore do not call it "validation set"), because we are mostly concerned about the accuracy on equilibrium configurations.

Fitting and verification of the ensemble of mMTPs
The results of fitting and verification of the ensemble including five mMTPs are given in Table 2.It could be seen that the uncertainty in energy, force, and stress errors is small compared to the magnitude of the errors, and the fitting errors and the errors on the "verification" set are reasonable and typical for the periodic crystal systems.We use the fitted ensemble of five mMTPs for further computations of the values of interest because from Table 2 we see no overfitting of the ensemble of five mMTPs.However, it should be noted that the force error on the "verification" set is smaller than the fitting force error and the stress error on the "verification" set is two times larger than the fitting stress error.Both the difference between force and stress errors are beyond the 95 % confidence interval.The reason may be that we constructed the "verification" set without Step 3.1 as opposed to the training set.In order to check it out we carry out an additional fivefold cross-validation.

Fivefold cross-validation
For the fivefold cross validation we combine the training and "verification" sets and create the total set of 2328 configurations.Next we split this set into five non-overlapping parts and carry out cross-validation, i.e. we fit mMTPs on four parts and validate it on the fifth part.The results are given in Table 3. From the table we conclude that the fitting and validation errors are close to each other and are within 99% confidence interval, i.e., the training set was constructed correctly.

Formation energy, lattice parameter, and total magnetic moment for different Fe-Al compounds
For testing the predictive power of the ensemble of fitted mMTPs we compare the formation energies, the equilibrium lattice parameters, and the total magnetic moments of the unit cell predicted with the mMTPs and DFT (implemented in ABINIT) for different compounds of Fe-Al.Uncertainty of the errors estimation in all the mentioned quantities in figures and tables are provided with the 95 % confidence interval (i.e., 2-σ interval).
Formation energy is given in Fig. 2. The ensemble of mMTPs correctly reproduces the formation energy trend calculated with DFT: it decreases when the concentration of Al increases.The maximum error of the formation energy prediction is around 20 meV/atom for the Fe 8 Al 8 −13682'4'5'7' structure.From Fig. 2 we find the configurations with minimum energy for a given concentration of Al, i.e., closest to the convex hull: pure bcc-Fe, Fe In Fig. 3 we compare the equilibrium lattice parameters calculated with the mMTPs and DFT.As for the formation energies, the mMTP-and DFT-provided equilibrium lattice parameters are close to each other, the maximum error in their prediction is about 0.014 Å(for the Fe 11 Al 5 -12457 structure).We also observe the anomalous lattice parameter/composition dependence in the Fe-Al structures.We, first, see the nearly linear increase of the lattice parameter up to the Al concentration of 18.75 %.When the Al concentration is between Table 2. Fitting root-mean-square errors (RMSEs) and RMSEs obtained on the "verification" set for the ensemble of five mMTPs and uncertainty of the errors estimation (we provide it with the 95 % confidence interval, i.e., 2-σ interval).The obtained errors are reasonable and typical for the periodic crystal systems.The ensemble of fitted mMTPs was not overfitted.In Fig. 3 we also provide the experimental lattice parameters from the paper 39 .The experimental work 39 has somewhat different lattice parameter dependencies for the differently processed samples (cast &quenched vs crushed ones), yet they both are anomalous, i.e., there is no linear dependency at more than 18.75 % Al concentration for the crushed samples (and at more than 31.25 % Al concentration for the cast &quenched ones).).The theoretical (mMTP and DFT) lattice parameters are close to each other: the maximum error in their prediction is 0.014 Å.Both the theoretical and the experimental dependencies of lattice parameters on Al concentration are anomalous: they are nonlinear at more than 18.75 % of Al concentration (at more than 31.25 % for the crushed samples).
Compositional dependencies of the total magnetic moment of the 16-atomic unit cell divided by the number of Fe atoms obtained with the ensemble of mMTPs and DFT are shown in Figure 4.The overall agreement between the mMTPs and DFT is very good: we see that the total magnetic moment of the unit cell decreases when the concentration of Al increases.The values of the total magnetic moments obtained with the mMTPs and DFT are also close to each other except for the Fe 8 Al 8 -12345678 structure: all the fitted mMTPs in the ensemble give zero magnetic moment whereas DFT gives 0.7 µ B .Except for this discrepancy for the Fe 8 Al 8 -12345678 structure, from the results of this subsection we conclude that the ensemble of mMTPs fitted to the DFT data essentially reproduces the variations in formation energies, lattice parameters, and total magnetic moments calculated with DFT.

Conclusion
In this paper we proposed the machine-learning interatomic potential with magnetic degrees of freedom (magnetic Moment Tensor Potential, mMTP) for prediction the properties of magnetic alloys.This potential was trained on data obtained with the recently developed method of cDFT calculations 19 that allows us to compute energies of configurations with non-equilibrium (excited) magnetic moments and, thus, to consider magnetic moments as degrees of freedom along with atomic positions, atomic types, and lattice vectors.We verify the developed magnetic multi-component machine-learning potentials on the Fe-Al system.We, first, created a training set including fully equilibrium atomic positions, lattice vectors, magnetic moments and the perturbed ones for different concentrations of Fe and Al in the Fe-Al system.Next, we fitted the ensemble of five mMTPs to the cDFT data and compared the dependencies of formation energies, equilibrium lattice parameters, and total magnetic moments of unit cell on the concentration of Al atoms predicted with the ensemble of mMTPs and DFT.We concluded that the mentioned mMTP and DFT differences are minor.Both mMTPs and DFT reproduced anomalous volume-composition dependence in the Fe-Al system obtained theoretically in the previous studies and has been experimentally observed.The main difference between the mMTP and DFT results was found for the Fe 8 Al 8 -12345678 structure: the ensemble of mMTPs gave the local minimum with zero magnetic moments for the Fe atoms whereas DFT predicted the minimum with magnetic moments of 0.7 µ B for the Fe atoms.Nev- ertheless, the rest of the results obtained with DFT and mMTP are in good correspondence.
In future, we are planning to develop an active learning algorithm for mMTP.Our confidence of developing an efficient active learning algorithm stems from the fact that with cDFT we are able to treat magnetic moments on the same footing as atomic positions; and for learning on atomic positions/geometries various successful active learning algorithms already exist.With active learning, we see in principle no obstacle in applying our methodology to predict material defect properties of other multi-component systems, e.g., of high-entropy alloys.In Ref. 40 we have developed an MTP-based algorithm for computing stacking fault energies, surface energies, and elastic constants, for (non-magnetic) bcc random alloys and used it to predict ductility of Mo-Nb-Ta over the entire composition space.Extending it to the case of magnetic alloys should be straightforward once we have developed an active learning algorithm, and this will allow us to screen for new materials with exotic mechanical properties over a much larger space of (magnetic and non-magnetic) metallic alloys than the space that can currently be approached with todays's state-of-the-art methods.Finally, active learning will allow us to train mMTP during molecular dynamics simulations and apply it to investigating the processes and predicting the properties of magnetic materials at finite temperature.The concept of magnetic multi-component Moment Tensor Potential (mMTP) presented in the current research is based on the previously developed non-magnetic MTP for multi-component systems 41,42 and magnetic MTP for single-component systems 35 .
The mMTP potential is local, i.e., the energy of the atomistic system is a sum of energies of individual atoms: where i stands for the individual atoms in an N a -atom system.We note that any configuration includes lattice vectors L = {l 1 , l 2 , l 3 } , atomic positions R = {r 1 , . . ., r N a } , types Z = {z 1 , . . ., z N a } (we also denote N types by the total number of atomic types in the system), and magnetic moments M = {m 1 , . . ., m N a } .The energy of the atom E i , in turn, has the form: where ξ = {ξ α } are the "linear" parameters to be optimized and B α are the so-called basis functions, which are contractions of the descriptors 25 of atomistic environment n i , yielding a scalar.The α max parameter can be changed to provide potentials with different amount of parameters 35 .
The descriptors are composed of the radial part, i.e., the scalar function depending on the interatomic distances and atomic magnetic moments, and the angular part, which is a tensor of rank ν: where n i stands for the atomic environment, including all the atoms within the R cut distance (or less) from the central atom i, µ is the number of the radial function, ν is the rank of the angular part tensor, |r ij | is the distance between the atoms i and j, z i and z j are the atomic types, m i and m j are the magnetic moments of the atoms.
The radial functions are expanded in a basis of Chebyshev polynomials: Here c = {c ζ ,β,γ µ,z i ,z j } are the "radial" parameters to be optimized, each of the functions φ ζ (|r ij |) , ψ β (m i ) , ψ γ (m i ) is a Chebyshev polynomial of order ζ , β and γ correspondingly, taking values from −1 to 1.The function φ ζ (|r ij |) yields the dependency on the distance between the atoms i and j, while the functions ψ β (m i ) and ψ γ (m j ) yield the dependency on the magnetic moments of the atoms i and j, correspondingly.The arguments of the functions φ ζ (|r ij |) are on the interval (R min , R cut ) , where R min and R cut are the minimum and maximum distance, corre- spondingly, between the interacting atoms.The functions ψ β (m i ) and ψ γ (m j ) are of the same structure, which we explain for the case of the former one.The argument of the function ψ β (m i ) is the magnetic moment of the atom i, taking the values on the (−M z i max , M z i max ) interval.The value M z i max itself depends on the type of atom z i , and is determined as the maximal absolute value of the magnetic moment for atom type z i in the training set.Similar to the conventional MTP, the term (R cut − |r ij |) 2 provides smooth fading to 0 when approaching the R cut distance, in accordance with the locality principle (1).
We note that magnetic degrees of freedom m i from (4) are collinear, i.e., they can take negative or positive values as projection onto the Z axis (though the choice of the axis is arbitrary).This way, in comparison to nonmagnetic atomistic systems with N atoms, in which the amount of degrees of freedom equals 4N (namely 3N for coordinates and N for types), for the description of magnetic systems additional N degrees of freedom are introduced, standing for the magnetic moment m i of each atom.The amount of parameters entering the radial functions (Eq.4) also increases in mMTP compared to the conventional MTP 41,42 .Namely, in MTP this number equals Thus, if we take N ψ = 2 (which is used in the current research), the amount of the parameters entering the radial functions would be four times more in mMTP then in MTP.
We denote all the mMTP parameters by θ = {ξ , c} and the total energy (1) of the atomic system by

Magnetic symmetrization of mMTP
The tensor (Eq.( 4)) includes collinear magnetic moments in its functional form.However, it is not invariant with respect to the inversion of magnetic moments, i.e., E(θ; M) � = E(θ; −M) , while both original and spin-inverted configurations should yield the same energy due to the arbitrary orientation of the projection axis, which we further call the magnetic symmetry.
We use data augmentation followed by explicit symmetrization with respect to magnetic moments to train a symmetric mMTP as we discuss below.Assume we have K configurations in the training set with DFT energies E DFT k , forces f DFT i,k , and stresses σ DFT ab,k ( a, b = 1, 2, 3 ) calculated.We find the optimal parameters θ (fit mMTP) by minimizing the objective function: where w e , w f , and w s are non-negative weights.By minimizing ( 5) we find such optimal parameters θ that yield E k ( θ ; M) ≈ E k ( θ ; −M) , k = 1, . . ., K (the same fact takes place for the mMTP forces and stresses), i.e., we symmetrize the training set to make mMTP learn the required symmetry from the data itself-this is called data augmentation.
Next, we modify mMTP to make the energy used for the simulations (e.g., relaxation of configurations) to satisfy the exact symmetry: That is, we substitute the mMTP energy (1) into ( 6) and get a functional form which satisfies the exact identity E symm ( θ; M) = E symm ( θ; −M) for any configuration.We also note that E( θ ) ≈ E symm ( θ ).

Constrained density functional theory calculations
We use the cDFT approach with hard constraints (i.e., Lagrange multiplier) as proposed by Gonze et al. in Ref. 19 .One way to formulate it is to first note that in a single-point DFT calculation we minimize the Kohn-Sham total energy functional E[ρ; R] with respect to the electronic density ρ = ρ(r) (here ρ combines the spin-up and spin-down electron densities), keeping the nuclei position R fixed.In other words, we solve the following minimization problem: and from the optimal ρ * = arg minE[ρ; R] we can, e.g., find magnetization m(r) = ρ * + − ρ * − , where the subscripts denote the spin-up ( + ) and spin-down (−) densities.The magnetic moment of the ith atom can be found by integrating m(r) over some (depending on the partitioning scheme) region around the atom: Since the minimizer ρ * depends on R , m i are also the functions of R.
According to the cDFT approach 19 , we now formulate the problem of minimizing E[ρ; R] in which not only R, but also ρ is allowed to change only subject to constraints (7): The algorithmic details of how this minimization problem is solved, and how the energy derivatives (forces, stresses, torques) are computed, are described in detail in Ref. 19 .

Computational details
We used the ABINIT code 43,44 for DFT (and cDFT recently developed and described in Ref. 19 ) calculations with 6 × 6 × 6 k-point mesh and cutoff energy of 25 Hartree (about 680 eV).We utilized the PAW PBE method with the generalized gradient approximation.We applied constraints on magnetic moments of all atoms during cDFT calculations.
We fitted an ensemble of five mMTPs with 415 parameters in order to quantify the uncertainty of mMTPs predictions.For each mMTP we took R min = 2.1 Å, R cut = 4.5 Å, M Al max = 0.1 µ B , and M Fe max = 3.0 µ B .The weights in the objective function (5) were w e = 1 , w f = 0.01 Å 2 , and w s = 0.001.

Figure 2 .
Figure 2. Formation energies for different compounds of Fe-Al predicted with the ensemble of mMTPs and DFT.The trend of formation energy was correctly reproduced with the ensemble of mMTPs, maximum error in formation energy prediction is 20 meV/atom for the Fe 8 Al 8 −13682'4'5'7' structure.

Figure 3 .
Figure 3. Lattice parameter dependencies on concentration of Al in the Fe-Al compounds computed with mMTP and DFT, and obtained experimentally (the experimental data are taken from Ref.39 ).The theoretical (mMTP and DFT) lattice parameters are close to each other: the maximum error in their prediction is 0.014 Å.Both the theoretical and the experimental dependencies of lattice parameters on Al concentration are anomalous: they are nonlinear at more than 18.75 % of Al concentration (at more than 31.25 % for the crushed samples).

Figure 4 .
Figure 4. Total magnetic moment of unit cell divided by the number of Fe atoms for different compounds of Fe-Al computed with the ensemble of mMTPs and DFT.Both the mMTPs and DFT predict the decrease in total magnetic moment when the concentration of Al increases.