This section describes how we compiled the dataset using ab initio calculations. Then, we propose a method to determine whether the dopant atom is more likely to substitute Fe or Rh from the ab initio calculations. Lastly, we show that the favorability of substitution is independent of the dopant concentration or its position in the lattice and propose a model that relies on the dopant’s atomic radius, electronegativity, and valency.Supercells and dopantsIn this part, we describe the part of the configuration space that will be studied. We start with a standard Fe–Rh cubic lattice structure. The primitive cell of this alloy is defined with orthogonal lattice vectors equal to \(\left( l_0 \approx 2.99 {\mathring{\text{A}}} \right) \) in which one Fe atom and one Rh atom are placed with coordinates (0, 0, 0) and (0.5, 0.5, 0.5), respectively, in the lattice basis. We study doped FeRh alloys where doping occurs alongside a vacancy defect. To achieve that, we remove one atom of Fe and one atom of Rh from the lattice structure and then substitute one of the vacancies with a dopant atom. This procedure ensures that the solution’s atomic composition does not depend on the substitution type (whether the atom of Rh or the atom of Fe is substituted), which allows for the straightforward comparison of the per-atom energy.To calculate the concentration of the dopant, we use \(c = \frac{1}{N}\), where N is the total amount of atoms in a lattice structure excluding the vacancies. This work considers dopant concentration levels in the alloy of 1% to 5%. To obtain the desired concentration percentages, we expanded the original \((l_0, l_0, l_0)\) configuration by periodically repeating our initial primitive cell in space. We assume that there is only a single dopant atom and a single vacancy in each resulting configuration. For example, to obtain a dopant concentration of approximately 5% percent, we repeat the original lattice three times along one dimension and two times along two other dimensions. We call the resulting lattice \((3l_0, 2l_0, 2l_0)\) Small Supercell (SS). To get lower percentages of dopant concentration, we perform an analogous procedure and get Medium Supercell (MS) and Large Supercell (LS). The details on these configurations can be found in Table 1. Our work investigates the following \(M=16\) dopants: Ir, Co, Cr, Ru, Pd, Pt, Mn, Ni, Te, Ti, Ta, Re, Nb, Os, Sn, Sb.Table 1 Atomic composition of considered supercells.Note that the formation of interstitial solid solutions is unlikely for these transition metals since the atomic radii of these elements differ by not more than 5% to the atomic radius of Fe and Rh. Thus, according to the assumptions of Hume-Rothery rules43, all of them form substitutional solid solutions and occupy sites in the nodes of the crystalline lattice of FeRh. As indicated by previously published experimental studies, the first eight chosen components are elements that maximally preserve the significant magneto-thermal properties of binary alloys22,27,31.Fig. 13D representations of all possible non-equivalent double substitutions. Red atoms correspond to Fe; blue atoms correspond to Rh. The dotted bounding box represents the volume spanned by lattice vectors. Yellow atoms indicate pairs of atoms that are first removed and then successively replaced by the dopant atom.Non-equivalent configurationsGenerally, there are three different scenarios in double substitutions: we can either substitute two Fe atoms (Fe–Fe type), two Rh atoms (Rh–Rh type), or both Fe and Rh atoms (Fe–Rh type). For the reasons described in Section “Determining the replacement type”, we only consider Fe–Rh type substitutions in this work. Our goal is to study the thermodynamic properties of each symmetrically inequivalent realization of the structure, which can correspond to fixed atomic occupancies and, consequently, hinder applications of the virtual crystal44 and coherent potential approximations45.Additionally, we obtain different resulting conformations by varying the locations from which Fe and Rh atoms are removed. Note that many double substitutions are equivalent to each other when translations and reflections are considered. To find all possible unique double substitutions, we employ Supercell software46 to obtain all non-equivalent pairs of atoms. We have two options for placing the dopant for any such pair, so the total number of double substitutions \(K = 2 * [\text {total number of non-equivalent pairs}]\). We illustrate non-equivalent double substitutions in Fig. 1. There are \(K_{SS} = 4\) unique double substitutions for small supercell, \(K_{MS} = 6\) unique double substitutions for medium supercell, and \(K_{LS} = 8\) unique double substitutions for large supercell. The total number of double substitutions is thus \((K_{SS} + K_{MS} + K_{LS}) * M = 288\).Overall databaseTo complete the database, we extend it with additional structures. First, we include the baseline supercells without any substitutions, thus adding three additional conformations to the database. Second, we include the so-called “single substitutions,” in which a single Fe (Rh) atom is replaced with a dopant. The total number of single substitutions for all considered supercells is \(2 * 3 (\text {number of supercells}) * M = 96\). Additionally, it is interesting to consider the case when some of the Fe and Rh atoms are interchanged with each other. Since this procedure could be done for each combinatorially non-equivalent case, thus we increase our database by \(\frac{1}{2}(K_{SS} + K_{MS} + K_{LS}) = 9\) entities. Our dataset consists of \(288 + 3 + 96 + 9 = 396\) initial configurations. To get the final dataset, we obtain geometry optimization trajectories for all initial configurations (see Section “Determining the replacement type”).Fig. 2Absolute energy differences \(|\Delta E^{D, SC}|\) per atom in relaxed states for different dopants and supercell types. All the values are in milielectronvolts (meV). The red and blue colors indicate double substitutions where \(E^{D, SC} < 0\) and \(E^{D, SC} > 0\) respectively.Determining the replacement typeIn this section, we propose a method to determine whether the dopant atom is more likely to substitute Fe or Rh in case of the double vacancy defect through ab initio calculations (see Section “Additional details” for the details on calculations). For a given dopant type D and supercell type SC, we divide the set of all non-equivalent double substitutions into two subsets \(SC_{Fe}(D)\) and \(SC_{Rh}(D)\), where an atom of Iron and atom of Rh is substituted by a dopant atom, correspondingly. Next, for each configuration c, we calculate the per-atom energy \({\text {E}}(c)\) in the relaxed state. To get the relaxed structure and the potential energy for a substitution, we perform geometry relaxation using Vienna Ab initio Simulation Package(VASP)47,48,49 (see Section “Additional details”). We simultaneously perform volume relaxation to account for different sizes of dopant atoms. Finally, we compare the minimum energies of two types: \(\min \limits _{c \in SC_{Fe}(D)} {\text {E}}(c)\) and \(\min \limits _{c \in SC_{Rh}(D)} {\text {E}}(c)\). The configuration with lower energy indicates a more desirable substitution type for the dopant. Note that the direct comparison of energies for two configurations is only possible for structures with the same atomic composition. This is the main reason we only consider double substitutions with one atom of each Fe and Rh removed. Such alloys can be denoted as \(\text {Fe}_{x-1}\text {Rh}_{x-1}\text {D}_1\), where x is the total number of Fe(Rh) atoms in a supercell, and D denotes a dopant.To illustrate the preferred replacement types for various dopants, we calculate the difference in minimal total energies after the relaxation between the cases where the dopant atom replaces either Fe or Rh: \(\Delta E^{D, SC} = \min \limits _{c \in SC_{Rh}(D)} {\text {E}}(c) – \min \limits _{c \in SC_{Fe}(D)} {\text {E}}(c)\), where D stands for a dopant type, and SC stands for the supercell type. The results are summarized in Fig. 2. We color entries in Fig. 2 depending on the \(\Delta E^{D, SC}\) sign. If \(\Delta E^{D, SC} < 0\) (replacing the atom of Fe is more energetically preferable for the dopant atom), the color is set to purple; if \(\Delta E^{D, SC} > 0\) (replacing the atom of Rh is more energetically preferable for the dopant atom), the color is set to blue. Qualitatively, we can observe that Cr, Mn, Te, Ti, Ta, Re, Nb, Sn, and Sb tend to replace Fe, while Ir, Co, Ru, Pd, Pt, Ni, and Os are more likely to substitute Rh. Moreover, for each dopant, we observe that the concentration determined by the supercell type does not influence (with precision up to 1.25 meV) the preferability of the substitution.Rules for predictionFig. 3Barplots of atomic radii, electronegativities, and maximum valencies of dopants. Dopants replacing Fe, according to the proposed computational method, are colored red, while elements replacing Rh are colored blue. Dashed lines indicate the atomic radii, electronegativities, and maximum valencies of Fe and Rh.In this work, we are concerned with substitutional solid solutions of three metals. There is a basic set of rules50 that indicates when an element could dissolve in a metal, forming a substitutional solid solution. This set of rules operates with three numerical features of atoms: empirically measured covalent radius (also known as Slater radius51), electronegativity by Pauling scale52, and maximum valency. However, this set of rules was initially derived for two-element alloys, whereas our work studies three-element alloys. Moreover, these rules do not provide an intuition on this manuscript’s central question: whether the dopant atom is more likely to substitute Fe or Rh in case of the double vacancy defect in the alloy.We hypothesize that the favorability of the substitution does not depend on the concentration of the dopant or its position in the lattice. We propose using the aforementioned numerical features to infer rules to determine the favorability directly from ab initio data using classification decision trees53. To train a classification decision tree, we must first map dopant atoms to vectors in the feature space defined by atomic radius, electronegativity, and valency of the atom. This procedure results in a feature matrix \(\textbf{X} \in \textbf{R}^{M \times 3}\), which is visualized in Fig. 3.To define the targets for classification, we use ab initio calculations described in Section “Determining the replacement type” and presented in Fig. 2. We say that an element belongs to class “0” (replaces Fe) if it replaces Fe in the majority of concentrations (see Fig. 2). Analogously, we say that an element belongs to class “1” (replaces Rh) if it replaces Rh in the majority of concentrations. This results in a target vector \(\textbf{y} \in \{0, 1\}^M\).To infer the classification rules, we use the DecisionTreeClassifier from the sklean package54 and fit it using \(\textbf{X}\) as the feature matrix and \(\textbf{y}\) as the target vector. We specifically limit the maximum depth of the tree to 3 to avoid complex rules. The resulting tree is visualized in Fig. 4. Every leaf node is a terminal node with a certain class (“replaces Fe” or “replaces Rh”) assigned to it. The class of the leaf node is determined during the tree fitting as the class that the majority of the elements in the leaf have. The leaf nodes also contain the Gini impurity53 value that indicates how homogeneous the objects in the leaf are in terms of their target value. The Gini impurity of 0 means that all objects in the leaf belong to the same class. Note that the Gini impurity equals 0 in all leaf nodes, so the decision tree perfectly separates all dopants.Every non-leaf node contains a rule consisting of a feature and a threshold. For example, the root node contains the following rule: If the dopant’s electronegativity is larger than 2.15, the dopant replaces Rh (the right child of the root node is the leaf node with class “1” assigned to it); otherwise, continue with the rule in the left child. Note that all non-leaf nodes contain rules with unique features, which indicates that all three features are essential to classify all the dopants perfectly.Fig. 4A visualization of the classification decision tree. A leaf node contains the Gini impurity value, the total amount of samples from data \(\textbf{X}\) that ended up in the leaf, class values for these samples (i.e., \(\text {values} = [0, 5]\) stands for 0 samples from class “0” and 5 samples from class “1”), and the class name of the majority of samples in the leaf. A non-leaf node additionally contains a rule by which to partition the samples.Additional detailsSpin-polarized total energy DFT calculations and total energy relaxations were executed using the Vienna Ab initio Simulation Package(VASP)47,48,49. We used a fully automatic generation scheme of the Bloch vectors with \(R_k = 20 {\mathring{\text{A}}} \) where \(R_k\) defines subdivisions \(N_i\) along the reciprocal lattice vectors \(\mathbf {b_i}\) in the following way: \(N_i = int(max(1, R_k|\mathbf {b_i}| + 0.5))\). The energy cutoff for the plane-wave basis set was set at 600 eV. For self-consistent field conditions, the criteria for breaking the cycle were that the total energy change and the band structure energy change between two steps are smaller than \(10^{-6}\) eV. Gaussian smearing was used for electronic occupancies with a smearing of 0.04 eV. All calculations were considered spin-polarized ones with initial magnetic moments equal to \(2\mu _b\) for each atom. Ionic relaxation was executed via the conjugate gradient algorithm until the norms of all the forces acting on the atoms were smaller than \(10^{-2}\) eV/Å. We used the default blocked-Davidson-iteration scheme. Finally, the Projector Augmented-Wave (PAW)55 pseudopotentials were used to reproduce the atomic core effects in the electronic density of the valence electrons. In the standard mode, VASP performs a fully relativistic calculation for the core electrons and treats valence electrons in a scalar relativistic approximation.Table 2 presents the averages and standard deviations of energies and maximum forces for each system, along with the number of structures analyzed. These properties characterize the systems at the final stage of ionic optimization. Notably, the average maximum force meets the criterion for halting the ionic relaxation, as the forces are smaller than \(10^{-2}\) eV/Å.Table 2 Averages and standard deviations of energies per atom (eV) and maximum forces (eV/Å) for each system, alongside the number of analyzed structures at the final stage of ionic optimization.
