Machine learning potentials for multicomponent systems: The Ti-Al binary system

Given a cutoff radius $r_{c}$, the short-range part of the potential energy for a structure, $E$, may be decomposed as

where $E^{(i)}$ denotes the contribution of atom $i$ within cutoff radius $r_{c}$. The atomic contribution to the potential energy can be referred to as the atomic energy. The atomic energy is then assumed to be expressed in a functional form of the neighboring atomic densities. Figure 1 shows a schematic illustration of the neighboring atomic density around atom $i$ within cutoff radius $r_{c}$ and its decomposition into the neighboring atomic densities of the elements. In a multicomponent system composed of elements $\{\rm A,\rm B,\cdots\}$, the atomic energy is written using functional $\mathcal{F}$ dependent on the element of atom $i$ as

where $\rho^{(i)}_{(s_{i},s)}$ denotes the neighboring atomic density of element $s$ ($s\in\{{\rm A},{\rm B},\cdots\}$) around atom $i$ of element $s_{i}$.

Subsequently, the neighboring atomic density of element $s$ around atom $i$ is expanded in terms of a basis set, because the expansion enables the functional form to be replaced with a function of its expansion coefficients. For a given basis set $\{b_{n}\}$, the neighboring atomic densities can be expanded as

where $a_{n,(s_{i},s)}^{(i)}$ denotes an order parameter characterizing the neighboring atomic density of element $s$ around atom $i$ of element $s_{i}$. Using the order parameters, the atomic energy may be rewritten as

Although function $F^{\prime}$ of Eqn. (4) depends on the element of atom $i$, it is convenient to introduce a unified function that is independent of the element by combining functions $F^{\prime}$ for all elements. In this study, a unified function for the atomic energy is formulated using order parameters defined for unordered pairs of elements. This means that order parameter $a_{n,(s_{i},s_{j})}^{(i)}$ of atom $i$ and its swapped order parameter $a_{n,(s_{j},s_{i})}^{(j)}$ of atom $j$ are considered as the same variable in the unified function. They are represented by $a_{n,\{s_{i},s_{j}\}}^{(i)}$ and $a_{n,\{s_{i},s_{j}\}}^{(j)}$, respectively. Defining order parameter $a_{n,\{s_{1},s_{2}\}}^{(i)}$ to be zero if $s_{i}$ is not included in $\{s_{1},s_{2}\}$ ($s_{i}\notin\{s_{1},s_{2}\}$), the atomic energy is written in an independent form of the element of atom $i$ as

where all combinations with the replacement of elements are enumerated for each $n$.

Moreover, an arbitrary rotation leaves the atomic energy invariant, although it generally changes the neighboring atomic densities and their order parameters. Therefore, the atomic energy is required to be a function of O(3) invariants $\{d_{n^{\prime}}^{(i)}\}$ as

where the invariants are derived from the order parameters $\{a_{n,\{s_{1},s_{2}\}}^{(i)}\}$. Hereafter, the invariants representing the neighboring atomic density are referred to as “structural features”. The present formulation is useful for deriving potential energy models, the accuracy and computational efficiency of which can be controlled by the selections of structural features and function $F$.

In this study, the neighboring atomic densities are expanded in terms of a basis set composed of radial functions or a basis set composed of products of radial functions and spherical harmonics Bartók et al. (2013); Seko et al. (2019), although it is also possible to use other basis sets in principle. When the neighboring atomic density is expanded in terms of radial functions $\{f_{n}\}$, the neighboring atomic density around atom $i$ is expressed as

where $r$ denotes the distance from the position of atom $i$. Since order parameter $a^{(i)}_{n,\{s_{i},s\}}$ is invariant for the O(3) group, it can be a pairwise structural feature denoted as

where $t$ identifies the unordered pair of elements, i.e., $t\in\{\{A,A\},\{A,B\},\cdots\}$.

When the neighboring atomic density is expanded in terms of products of radial functions $\{f_{n}\}$ and spherical harmonics $\{Y_{lm}\}$, the neighboring atomic density of element $s$ at a position $(r,\theta,\phi)$ in spherical coordinates centered at the position of atom $i$ is expressed as

where order parameter $a^{(i)}_{nlm,\{s_{i},s\}}$ is component $nlm$ of the neighboring atomic density. Since the order parameters are not generally invariant for the O(3) group, polynomial invariants of the order parameters are adopted as structural features. A $p$th-order polynomial invariant for a radial index $n$ and a set of pairs composed of the angular number and the element unordered pair $\{(l_{1},t_{1}),(l_{2},t_{2}),\cdots,(l_{p},t_{p})\}$ is defined as a linear combination of products of $p$ order parameters, expressed as

Nonzero polynomial invariants are derived from the sets satisfying the condition that the intersection of $t_{p}$ is not the empty set,

Therefore, possible sets of pairs composed of the angular number and the element unordered pair $\{(l_{1},t_{1}),(l_{2},t_{2}),\cdots,(l_{p},t_{p})\}$ are enumerated for a given maximum angular number to obtain the entire set of polynomial invariants. An example where the intersection of element pairs becomes the empty set is the case that a polynomial invariant is composed of order parameters with $t_{1}={\rm\{A,A\}}$ and $t_{2}={\rm\{B,B\}}$. Such polynomial invariants are eliminated.

A coefficient set $\{c^{l_{1}l_{2}\cdots l_{p},(\sigma)}_{m_{1}m_{2}\cdots m_{p}}\}$ is independent of the radial index $n$ and the element unordered pair $t$. Therefore, linearly independent coefficient sets are obtained using the group-theoretical projector operation method for a given set $\{l_{1},l_{2},\cdots,l_{p}\}$, as proposed in Ref. Seko et al., 2019, ensuring that the linear combinations are invariant for arbitrary rotation. In terms of fourth- and higher-order polynomial invariants, multiple invariants are linearly independent for most of the set $\{l_{1},l_{2},\cdots,l_{p}\}$, which are distinguished by index $\sigma$ if necessary. Note that the second- and third-order invariants are equivalent to a multicomponent extension of the angular Fourier series and the bispectrum reported in the literature, respectively Kondor (2008); Bartók et al. (2013).

In this study, a finite set of Gaussian-type functions is adopted as radial functions in basis sets to expand the neighboring atomic density, expressed as

where $\beta_{n}$ and $r_{n}$ are parameters. The cutoff function $f_{c}$ is given by a cosine-based function as

The order parameters of atom $i$ and element pair $\{s_{i},s\}$ are approximately estimated from the neighboring atomic density of element $s$ around atom $i$ as

where $(r_{ij},\theta_{ij},\phi_{ij})$ denotes the spherical coordinates of neighboring atom $j$ centered at the position of atom $i$. Although the Gaussian-type radial functions are not orthonormal, such an approximation of the order parameters is acceptable in the present polynomial-based framework, as also discussed in Ref. Seko et al., 2019.

Polynomial models are here employed to represent the atomic energy as a function of structural features. Given a set of structural features $D=\set{d_{1},d_{2},\cdots}$, polynomial functions are written as

where $w$ denotes a regression coefficient. Although the polynomial functions are described by all combinations of structural features, only nonzero polynomial terms are retained, which is analogous to the enumeration of nonzero polynomial invariants. In a multicomponent system, a structural feature is composed of order parameters, each of which has an attribute on the element unordered pair $t$. Therefore, when the element pair of the $p^{\prime}$th order parameter in structure feature $d_{i^{\prime}}$ of a polynomial term is denoted by $t_{i^{\prime},p^{\prime}}$, a nonzero polynomial term satisfies the condition that the intersection of $\{t_{i^{\prime},p^{\prime}}\}$ is not the empty set:

For example, the intersection of element pairs becomes the empty set if a polynomial term is composed of structural features with $t_{i^{\prime},p_{i}^{\prime}}={\rm\{A,A\}}$ and $t_{j^{\prime},p_{j}^{\prime}}={\rm\{B,B\}}$. Such polynomial terms are eliminated from the polynomial functions.

The following polynomial models for the atomic energy are systematically applied to obtain Pareto optimal MLPs as described in Sec. IV.1. The first model is a polynomial of pairwise structural features. When a set of pairwise structural features is described as

the first model is expressed as

The first model includes the special case that only powers of the pairwise structural features are considered, which was introduced for elemental systems in Refs. Seko et al., 2014 and Seko et al., 2015. The first model can also be regarded as a straightforward extension of embedded atom method (EAM) potentials Takahashi et al. (2017).

The second model for the atomic energy is a linear polynomial of polynomial invariants given by Eqn. (10). A set of polynomial invariants is described as

where a set of $p$th-order polynomial invariants is denoted by

A second-order invariant is identified with a single $l$ value because second-order linear combinations are invariant only when $l_{1}=l_{2}$ El-Batanouny and Wooten (2008); Tolédano and Tolédano (1987). The second model is then written as

which was introduced in Ref. Seko et al., 2019 for elemental systems. Note that a linear polynomial model with up to third-order invariants is equivalent to a spectral neighbor analysis potential (SNAP) Thompson et al. (2015), expressed as

An extension of the second model is a polynomial of polynomial invariants described as

Note that a quadratic polynomial model of polynomial invariants up to the third order is equivalent to a quadratic SNAP Wood and Thompson (2018b). Other extended models are also introduced, which are given by

They are decomposed into a linear polynomial of structural features and a polynomial of a subset of the structural features.

The atomic energy in all the models is measured from the sum of the energies of isolated atoms. Moreover, the forces acting on atoms and the stress tensors in a multicomponent system can be derived in a similar manner to those in the elemental system derived in Ref. Seko et al., 2019. Note that the above polynomial function forms were also applied to develop MLPs for elemental systems included in the Machine Learning Potential Repository Seko (2020); Mac .

The polynomial models given in Sec. II.3 are systematically applied to develop MLPs in the Ti-Al binary system. Because the accuracy and computational efficiency of the MLPs strongly depend on several input parameters, a systematic grid search is performed to find their optimal values. The parameters in the grid search are the cutoff radius, the type of structural features, the type of potential energy model, the number of radial functions, the polynomial order in the potential energy model, and the truncation of the polynomial invariants, i.e., the maximum angular numbers of spherical harmonics, $\{l_{\rm max}^{(2)},l_{\rm max}^{(3)},\cdots,l_{\rm max}^{(p_{\rm max})}\}$, and the polynomial order of the invariants, $p_{\rm max}$.

Figure 2 (a) shows the distribution of MLPs obtained from the grid search. The root mean squared (RMS) error for the test dataset is used as an estimator of the accuracy of MLPs. The computational time indicates the elapsed time for a single point calculation normalized by the number of atoms. Figure 2 (a) also shows the Pareto optimal MLPs when optimizing both the accuracy and computational efficiency simultaneously. As can be seen in Fig. 2 (a), the accuracy and computational efficiency of MLPs are conflicting properties; hence, the Pareto optimal MLPs can be optimal ones with different trade-offs.

In performing an atomistic simulation, an appropriate MLP must be chosen from the Pareto optimal ones according to its target system and purpose. Therefore, a convenient score that can estimate the computational cost performance is required to find an MLP with high computational cost performance in a simplified manner. In this study, functions $t+\Delta E$ and $10t+\Delta E$ that should be minimized are introduced, where $t$ and $\Delta E$ denote the computational time with the unit of ms/atom/step and the RMS error with the unit of meV/atom, respectively.

Figure 2 (a) shows the MLP with the lowest RMS error and two Pareto optimal MLPs showing high computational cost performance. The MLP with the lowest RMS error is denoted by “MLP3”, whereas the two MLPs showing high computational cost performance are denoted by “MLP1” and “MLP2”. MLP1 and MLP2 are obtained by minimizing the scores $10t+\Delta E$ and $t+\Delta E$, respectively. As can be seen in Fig. 2 (a), they exhibit high computational efficiency without significantly increasing the RMS error. The model parameters and RMS errors of MLP1, MLP2, and MLP3 are listed in Table 2. The RMS errors of MLP2 and MLP3 are close, 2.09 and 1.67 meV/atom, respectively. The RMS error of MLP1 is 4.51 meV/atom, which is greater than those of MLP2 and MLP3, whereas MLP1 is ten times more computationally efficient than MLP2. As described in Table 2, all the MLPs are derived from polynomial invariants. However, MLP1 is developed only using the $l=0$ component of spherical harmonics. Therefore, MLP1 can be regarded as a polynomial model of pairwise structural features.

Figure 2 (b) shows the computational time required for a single point calculation using the EAM potential Zope and Mishin (2003), MLP1, MLP2, and MLP3. The computational time is evaluated for structures with up to 32,000 atoms constructed by the expansions of the four-atom unit cell of a AuCu-type ($L1_{0}$) structure with a lattice constant of 4 Å. The computational time is measured by implementing the present MLPs Lam in the lammps code lammps code ; Plimpton (1995). As can be seen in Fig. 2 (b), the EAM potential and the MLPs show linear scaling of the computational time with respect to the number of atoms. Therefore, the computational time normalized by the number of atoms can be an estimator of the computational efficiency of MLPs, and the computational time required for a simulation of $n_{\rm step}$ steps using a structure with $n_{\rm atom}$ atoms can be estimated as $t\times n_{\rm atom}\times n_{\rm step}$.

Figure 3 (a) shows the distribution of the prediction errors for structures in the training and test datasets. The degree of scattering in the distribution decreases in the order of MLP1, MLP2, and MLP3, which coincides with the decreasing order of the RMS errors shown in Table 2. Figure 3 (b) shows the distribution of the prediction errors for the equilibrium structure generators listed in Table 1. The distribution of the prediction errors for the structure generators indicates the predictive power for the cohesive energy. The prediction errors for some of the structure generators are significant in MLP1. On the other hand, the prediction errors are trivial for all the structure generators in MLP2 and MLP3, which indicates that MLP2 and MLP3 should have high predictive power for the cohesive energy in a wide variety of binary ordered structures.

The formation energy of a given ordered structure is more challenging to predict accurately than its cohesive energy because high predictive power is required not only for the ordered structure but also for the reference structures of hexagonal-close-packed (hcp) Ti and face-centered-cubic (fcc) Al. Furthermore, the local geometry relaxation for the ordered structure is crucial for evaluating the formation energy. Therefore, MLPs are needed to derive an accurate potential energy surface around the initial and equilibrium structures.

Tables 3 and 4 show the formation energies of selected ordered structures predicted using the EAM potential, MLP1, MLP2, and MLP3, compared with those predicted by DFT calculation. Note that the ordered structures cover a wide range of compositions. The RMS errors of the EAM potential, MLP1, MLP2, and MLP3 for the formation energy are 68.0, 13.8, 2.0, and 1.5 meV/atom, respectively. This reveals that MLP2 and MLP3 have high predictive power for the formation energy in a wide range of structures. The RMS error of MLP1 for the formation energy is significant as a consequence of the systematic deviation of the formation energy for the overall ordered structures. As can be seen in Tables 3 and 4, the formation energies of most of the ordered structures predicted using MLP1 are approximately 10 meV/atom lower than those predicted by DFT calculation, which originates from the fact that the prediction error of MLP1 for hcp-Ti is significant ($+$27.4 meV/atom). Moreover, the prediction error of the EAM potential is much more significant, and the EAM potential fails to reconstruct the hierarchy of the formation energies predicted by DFT calculation.

Tables 5, 6, and 7 show the lattice constants and the elastic constants of nine structures, i.e., TiAl ($L1_{0}$), TiAl ($B2$), TiAl ($B8_{1}$), Ti₃Al ($D0_{19}$), Ti₃Al ($D0_{22}$), Ti₃Al ($L1_{2}$), TiAl₃ ($D0_{22}$), TiAl₃ ($D0_{19}$), and TiAl₃ ($L1_{2}$), which are predicted using the EAM potential, MLP1, MLP2, and MLP3, along with the lattice constants and the elastic constants obtained by DFT calculation. The lattice constants predicted using the EAM potential deviate from those obtained by DFT calculation in most of the structures. This deviation may arise from the fact that the EAM potential was developed by fitting to experimental lattice constants Zope and Mishin (2003), excluding the descriptive power of the EAM potential model. In contrast, the MLPs can compute the lattice constants accurately in all of the structures. Regarding the elastic constants, the elastic constants predicted using the EAM potential and MLP1 are close to those obtained by DFT calculation in many structures; however, the EAM potential and MLP1 fail to predict the elastic constants accurately in a few structures. Also, the elastic constants predicted using MLP2 and MLP3 are almost the same as those obtained by DFT calculation in all structures.

The predictive power of the MLPs for phonon properties and thermal expansion is examined. The phonon properties and thermal expansion are calculated using a finite displacement method implemented in the phonopy code Togo and Tanaka (2015). Figure 4 (a) shows the phonon density of states for 13 structures predicted using the EAM potential, MLP1, MLP2, and MLP3. The supercells required to compute the phonon properties are constructed by the expansions of the conventional unit cells of the structures. The number of atoms included in the supercells ranges from 54 to 162. The EAM potential predicts the phonon density of states well in the low-frequency region for many structures, while the deviation from the DFT phonon density of states is large in the high-frequency region. Conversely, the phonon density of states predicted using the MLPs and those predicted by DFT calculation overlap for all the structures, particularly those predicted using MLP2 and MLP3.

Figure 4 (b) shows the phonon dispersion curves for TiAl (AuCu-type, $L1_{0}$) predicted using the EAM potential and the MLPs, compared with those predicted by DFT calculation. The deviation of the EAM phonon dispersions from the DFT phonon dispersions is significant in the high-frequency region. On the other hand, the phonon dispersions of MLP2 and MLP3 are consistent with those of DFT calculation.

Figure 4 (c) shows the temperature dependence of the thermal expansion, calculated using a quasi-harmonic approximation, in Ti₃Al ($D0_{19}$), TiAl ($L1_{0}$), and TiAl₃ ($D0_{22}$), which are experimentally observed in the Ti-Al binary system. As can be seen in Fig. 4 (c), the thermal expansion of the EAM potential differs from that of the DFT calculation in all the structures. On the other hand, MLP2 and MLP3 derive the temperature dependence of the thermal expansion accurately in all the structures, even though the thermal expansion is more challenging to predict accurately than the phonon density of states and the phonon dispersion curves. The accurate prediction of the thermal expansion indicates that MLP2 and MLP3 can accurately evaluate the volume dependence of the whole range of phonon frequencies.

The energy profile along the Bain path between the $\gamma$-TiAl ($L1_{0}$) structure and the $B2$ structure is calculated using the EAM potential and the MLPs, and compared with that obtained by DFT calculation. The structures required to compute the energy profile are obtained by transforming the $c/a$ ratio while keeping their volumes fixed to that of the equilibrium $L1_{0}$ structure. The energy profile is then evaluated by single-point calculations without geometry optimization for the structures. Figure 5 shows the energy profiles along the Bain path predicted using the EAM potential, the MLPs, and DFT calculation. The energy profiles along the Bain path of the MLPs, particularly MLP2 and MLP3, are consistent with that obtained by DFT calculation. The EAM profile is also close to the DFT profile for a $c/a$ ratio of 1.0–1.5, although the $c/a$ ratio of the equilibrium $L1_{0}$ structure of the EAM potential is slightly different from that obtained by DFT calculation. The behavior of the energy profile along the Bain path predicted by the EAM and the deviation of the $c/a$ ratio are the same as those discussed by Zope and Mishin in Ref. Zope and Mishin, 2003, who developed the EAM potential.

Computational models for the superlattice intrinsic stacking fault (SISF), the antiphase boundary (APB), and the complex stacking fault (CSF) are constructed using the procedure shown in Ref. Dumitraschkewitz et al., 2017. First, the supercell is constructed by the expansion along the $\langle 111\rangle$ direction of $\gamma$-TiAl with ideal cubic lattice parameters, and the resultant supercell is composed of 24 atoms. The equilibrium structure of the supercell is then obtained by local geometry optimization. A computational model with a stacking fault is generated by introducing a tilt into the equilibrium structure. A displacement vector defining a tilt is given by a linear combination of two vectors corresponding to the $\langle\bar{1}10\rangle/2$ and $\langle 11\bar{2}\rangle/2$ directions in ideal cubic $\gamma$-TiAl. The displacement vectors for the SISF, the APB, and the CSF are expressed as

respectively. Finally, single-point calculations are performed for the tilted structures.

Table 8 summarizes the stacking fault energies computed using the EAM potential, the MLPs, and DFT calculation. The cohesive energy of the equilibrium structure of $\gamma$-TiAl is used as a reference to compute the stacking fault energy. As can be seen in Table 8, the MLPs predict the stacking fault energies of the SISF, the APB, and the CSF accurately. On the other hand, the EAM potential lacks the predictive power for the stacking fault energy of the APB.

The stacking fault energy can be defined not only for these special stacking faults but also for other general stacking faults defined by displacement vectors. A collection of the excessive energies of the general stacking faults comprises a generalized stacking fault energy (GSFE) surface. The displacement vector identifying a tilt of the supercell for a general stacking fault is given by

where $u$ and $v$ denote the fractional coordinates defined by the vectors $\langle\bar{1}10\rangle/2$ and $\langle 11\bar{2}\rangle/2$, respectively.

Figure 6 (a) shows the profiles of the stacking fault energy along the path $\gamma$-TiAl $\rightarrow$ SISF $\rightarrow$ APB $\rightarrow$ CSF $\rightarrow$ $\gamma$-TiAl $\rightarrow$ APB, predicted using the EAM potential, the MLPs, and DFT calculation. The displacement vector changes continuously along the path. The stacking fault energy profiles of MLP2 and MLP3 are close to that obtained by DFT calculation. The MLP1 profile agrees with the DFT profile along a major part of the path, while it deviates from the DFT profile along the path between the APB and the CSF. The EAM profile also deviates from the DFT profile along the path between the SISF and the CSF. Figure 6 (b) shows the GSFE surface in $\gamma$-TiAl predicted using the EAM potential, the MLPs, and DFT calculation. The GSFE surfaces predicted using the MLPs and DFT calculation are similar, while that predicted using the EAM potential is different from that predicted by the DFT calculation. Thus, the present MLPs should have high predictive power for the stacking faults and related properties.