### Formalism

Two parts of a crystal can be shifted relative to each other by a vector (overrightarrow{r}) that is parallel to a glide plane as shown in Fig. 1a or the crystal can be cleaved by a distance *δ* perpendicular to the glide plane as in Fig. 1b. The gliding of the two parts of the crystal relative to each other results in a planar fault. The energy per unit area as the two halves of a crystal are shifted relative to each other by (overrightarrow{r}) is conventionally referred to as the GSF energy. It can be defined to be the glide energy either at fixed *δ* or at a value of *δ* corresponding to zero tractions perpendicular to the glide plane. Throughout we will use periodic boundary conditions.

The GSF energy of a single component crystal is a periodic function over the space of two-dimensional glide vectors (overrightarrow{r}). Translating the two halves of a crystal by a full lattice vector recovers the unfaulted bulk crystal. Each point in the two-dimensional glide space (overrightarrow{r}) that coincides with a translation vector will therefore have the same GSF energy. Figure 2a shows a (001) glide plane of a body-centered cubic (bcc) crystal. The corresponding two-dimensional glide vector space for the (001) glide plane of bcc is shown in Fig. 2b with glide vectors that recover the bulk bcc crystal shown as squares. Glide vectors that differ from a lattice translation produce a bicrystal with a planar fault. The introduction of a planar fault usually increases the energy of the crystal.

The GSF energies of multicomponent alloys differ from those of pure elements. The periodicity of the GSF energy that characterize a single component crystal is generally broken in a multicomponent alloy. This is illustrated for a particular ordering of red and green atoms on bcc in Fig. 3a. A glide of the upper half of the crystal relative to the lower half by a [100] translation vector of the underlying bcc crystal structure results in a different ordering of red and green atoms. Though the bcc crystal structure is recovered, a planar defect referred to as an anti-phase boundary has been created. The energy of the crystal before and after the glide by a bcc translation vector is therefore no longer the same since the arrangement of red and green atoms has changed. Two lattice translations along the (001) plane are required for the example of Fig. 3a to recover the original ordering. Figure 3b shows the symmetry in the two-dimensional space of glide vectors (overrightarrow{r}), with blue squares corresponding to the original ordering and yellow squares corresponding to orderings with an anti-phase boundary. For the particular ordering in Fig. 3a, glide vectors along the [010] direction do not change the ordering and the energy remains unchanged.

The example of Fig. 3a shows that it is necessary to track the degree of order in a multicomponent alloy since glides by a translation vector of the underlying parent crystal structure can change the arrangement of chemical species in the crystal. The state of ordering in a multicomponent solid can be described mathematically by assigning occupation variables to each site of its crystal. Consider a large crystal with *N* sites (assuming periodic boundary conditions), where each site can be occupied by one of two chemical components A or B. Any ordering of A and B atoms on this crystal can be represented with an occupation vector (overrightarrow{sigma }={{sigma }_{1},{sigma }_{2},cdots ,{sigma }_{{mathrm{N}}}}), where *σ*_{i} is an occupation variable that takes the value of +1 if site *i* is occupied by A and a value of −1 otherwise. The labels *i* refer to sites in a reference crystal that is not deformed. For the purpose of tracking the configuration of a deformed crystal, we map each site of the deformed crystal onto the nearest site of the reference crystal. A particular configuration of A and B atoms, ({overrightarrow{sigma }}_{1}), for example, may then be converted into a new configuration ({overrightarrow{sigma }}_{2}) upon application of a glide that coincides with an elementary translation vector of the underlying parent crystal structure. A second glide by another translation vector may convert ({overrightarrow{sigma }}_{2}) into ({overrightarrow{sigma }}_{3}). The changes in configuration upon application of glides coinciding with parent crystal translations can be represented in the two-dimensional glide space of (overrightarrow{r}) as schematically illustrated in Fig. 4a. Since the energy of the crystal depends on how the A and B atoms are arranged, it will also vary upon the application of a glide that is equal to a parent crystal translation. This is schematically illustrated in Fig. 4b.

While the gliding of two halves of a crystal by a translation vector of the parent crystal changes the configuration of the alloy and therefore its energy, most of the change in ordering is restricted to the vicinity of the glide plane. The local arrangements of A and B atoms far away from the glide plane are unaffected by the glide since those regions have simply been translated rigidly. Chemical interactions in an alloy typically decay over a distance of several nanometers when maintaining the solid in a constant state of strain. The contribution to the energy of the crystal due to a particular arrangement of A and B atoms far away from the glide plane will be identical in two configurations, ({overrightarrow{sigma }}_{1}) and ({overrightarrow{sigma }}_{2}), related by a glide since the local degree of ordering at those large distances from the glide plane are identical in both configurations. It is only in regions within the chemical interaction range of the glide plane where the local degree of order is different that an energy difference arises. This motivates the separation of the GSF energy into an average configurational energy and a “configurationally-resolved planar fault energy” (CRPF) that is a local excess energy. The energy, *E*, of a bi-crystal (shown schematically in Fig. 4b with an initial ordering (overrightarrow{sigma }) that is translated by (overrightarrow{r}) and separated by a distance *δ* can then be written as

$$E(overrightarrow{r},delta ,overrightarrow{sigma })={E}^{{rm{CRPF}}}(overrightarrow{r}-{overrightarrow{r}}_{1},delta ,{overrightarrow{sigma }}_{1})+{E}^{{rm{avg}}}(overrightarrow{r},{overrightarrow{sigma }}_{1},{overrightarrow{sigma }}_{2},{overrightarrow{sigma }}_{3}).$$

(1)

In this expression, ({E}^{{rm{CRPF}}}(overrightarrow{r}-{overrightarrow{r}}_{1},delta ,{overrightarrow{sigma }}_{1})) is the CRPF energy, with *σ*_{1} being one of the three nearest orderings on the perfect crystal for the glide vector (overrightarrow{r}) as schematically illustrated in Fig. 4a. The glide vector ({overrightarrow{r}}_{1}) converts (overrightarrow{sigma }) to ({overrightarrow{sigma }}_{1}). The average configurational energy, *E*^{avg}, is related to the energy of the three nearest configurations as schematically shown in Fig. 4a, and is defined as

$${E}^{{rm{avg}}}(overrightarrow{r},{overrightarrow{sigma }}_{1},{overrightarrow{sigma }}_{2},{overrightarrow{sigma }}_{3}),=,(1,-,{w}_{2},-,{w}_{3})E(overrightarrow{{sigma }_{1}}),+,{w}_{2}E(overrightarrow{{sigma }_{2}}),+,{w}_{3}E(overrightarrow{{sigma }_{3}}),$$

(2)

where (E({overrightarrow{sigma }}_{1})), (E({overrightarrow{sigma }}_{2})) and (E({overrightarrow{sigma }}_{3})) are the energies of the ({overrightarrow{sigma }}_{1}), ({overrightarrow{sigma }}_{2}), and ({overrightarrow{sigma }}_{3}) orderings in the perfect crystal. The weights *w*_{2} and *w*_{3} are related to the glide vectors by:

$$left[{w}_{2},,{w}_{3}right],=,{left[begin{array}{c}{overrightarrow{r}}_{12}\ {overrightarrow{r}}_{13}end{array}right]}^{-1}(overrightarrow{r},-,{overrightarrow{r}}_{1}),$$

(3)

where ({overrightarrow{r}}_{12}) is the glide vector relating configurations ({overrightarrow{sigma }}_{1}) and ({overrightarrow{sigma }}_{2}) and ({overrightarrow{r}}_{13}) connects ({overrightarrow{sigma }}_{1}) to ({overrightarrow{sigma }}_{3}).

Rigorous statistical mechanics calculations of the temperature and composition dependence of the GSF energy require the evaluation of (E(overrightarrow{r},delta ,overrightarrow{sigma })) across all possible decorations of the bi-crystal. This can be computationally intractable when using quantum mechanical techniques. Surrogate models informed from a small set of quantum mechanical calculations that accurately reproduce the bulk and CRPF energies for arbitrary configurations are thus needed to bridge the gap. In the rest of this section, we review the cluster expansion formalism to describe the configurational energy of crystalline solids and subsequently extend it to describe the CRPF energy as a function of configurational ordering.

As shown by Sanchez et al.^{24}, the configurational energy (E(overrightarrow{sigma })) in a multicomponent solid with a particular crystal structure can be expanded in terms of cluster basis functions according to

$$E(overrightarrow{sigma }),=,{V}_{0},+,mathop{sum }limits_{alpha }{V}_{alpha }{phi }_{alpha }(overrightarrow{sigma }),$$

(4)

where (E(overrightarrow{sigma })) is the energy of (overrightarrow{sigma }), *V*_{α} are expansion coefficients, referred to as effective cluster interactions, and the ({phi }_{alpha }(overrightarrow{sigma })) are cluster basis functions. For a binary alloy, the cluster basis functions are defined as

$${phi }_{alpha }(overrightarrow{sigma }),=,mathop{prod }limits_{{mathrm{j}}in alpha }{sigma }_{{mathrm{j}}},$$

(5)

where *α* refers to a cluster of sites in the crystal, such as pair clusters, triplet clusters etc. The *V*_{α} are determined by the chemistry of the alloy. In most alloys, chemical interactions decay beyond a maximum length and cluster size and the cluster expansion of Eq. (4) can be truncated. Clusters related to each other through a symmetry operation in the undecorated crystal have the same expansion coefficient. Strategies that rely on genetic algorithms^{25}, cross-validation^{26}, bayesian regression^{27}, neural networks^{28}, and quadratic programming^{29} have been succesfully applied to generate high-fidelity cluster expansion models trained to first-principles calculations. The resulting lattice models are typically used in conjunction with statistical mechanics tools such as Monte-Carlo simulations to calculate temperature and composition dependent thermodynamic properties of multi-component solids.

The cluster expansion of Eq. (4) is only valid for a fixed parent crystal structure. In the context of GSF energy surfaces, it can only be used to describe the energy of the bicrystal for glide vectors (overrightarrow{r}) that recover the underlying parent crystal structure. This includes the energies of (E({overrightarrow{sigma }}_{1})), (E({overrightarrow{sigma }}_{2})), and (E({overrightarrow{sigma }}_{3})) appearing in the expression of average configurational energy *E*^{avg} as defined by Eqs. (2), (3), and appearing in Eq. (1). We next extend the cluster expansion approach to describe the CRPF energy of Eq. (1).

We first simplify the problem by exploiting well-established analytical expressions of energy-versus-separation curves to describe the dependence of the CRPF on *δ*. For most metals, the energy versus separation curve can be accurately represented with the universal binding energy relation (UBER) of Rose et al.^{30,31} according to

$${E}^{{rm{CRPF}}},=,{E}_{0}^{{rm{CRPF}}}-2kappa left[-1,+,left(1,+,frac{delta,-,{delta }_{0}}{lambda }right)exp left(-frac{delta,-,{delta }_{0}}{lambda }right)right],$$

(6)

where ({E}_{0}^{{rm{CRPF}}}) is the CRPF energy at the equilibrium separation *δ*_{0}, 2*κ* is the surface energy at infinite separation, and *λ* is related to the curvature of the energy around the equilibrium separation. The parameters ({E}_{0}^{{rm{CRPF}}},,{delta }_{0},,kappa,) and *λ* are all functions of the configuration (overrightarrow{sigma }) and glide vector (overrightarrow{r}). While Eq. (6) is that for the UBER curve, alternate functional forms such as xUBER^{32,33} may also be employed.

The dependence of ({E}_{0}^{{rm{CRPF}}}), *κ*, *δ*_{0}, and *λ* on configuration (overrightarrow{sigma }) can be expressed as a cluster expansion. For example, ({E}_{0}^{{rm{CRPF}}}) can be written as

$${E}_{0}^{{rm{CRPF}}}(overrightarrow{r},overrightarrow{sigma }),=,{Gamma }_{0}(overrightarrow{r}),+,mathop{sum }_{alpha }{Gamma }_{alpha }(overrightarrow{r}){phi }_{alpha }(overrightarrow{sigma }),$$

(7)

where the cluster basis functions, *ϕ*_{α} are the same as those defined in Eq. (5), and where Γ_{α} are expansion coefficients that are functions of the glide vector, (overrightarrow{r}). Similar to the cluster expansion of Eq. (4), the expansion coefficients in Eq. (7) obey certain symmetry properties dictated by the space group of the undecorated bicrystal having undergone a glide (overrightarrow{r}). Since a glide of a bicrystal by (overrightarrow{r}) in general breaks symmetry, far fewer expansion coefficients will be equivalent by symmetry than for the cluster expansion of the undeformed parent crystal. For example, translation symmetry in directions perpendicular to the glide plane are lost upon application of a glide (overrightarrow{r}). Point clusters that are otherwise equivalent by symmetry in the perfect crystal, are no longer equivalent if they are at different distances from the glide plane. The same holds true for multi-body clusters.

The cluster expansions of the parameters ({E}_{0}^{{rm{CRPF}}}), *κ*, *δ*_{0}, and *λ* appearing in Eq. (6) extend over all clusters of the bicrystal. However, these cluster expansions should converge rapidly and only clusters within the chemical interaction range from the glide plane are likely necessary in a truncated cluster expansion. This becomes evident when rearranging Eq. (1) to isolate the CRPF energy according to

$${E}^{{rm{CRPF}}}(overrightarrow{r},-,{overrightarrow{r}}_{1},,delta ,,{overrightarrow{sigma }}_{1}),=,E(overrightarrow{r},,delta ,,overrightarrow{sigma }),-,{E}^{{rm{avg}}}(overrightarrow{r},,{overrightarrow{sigma }}_{1},,{overrightarrow{sigma }}_{2},,{overrightarrow{sigma }}_{3}).$$

(8)

The above equation shows that the contribution to the energy of the bi-crystal from regions far away from the glide plane is removed when subtracting off the weighted average energy ({E}^{{rm{avg}}}(overrightarrow{r},,{overrightarrow{sigma }}_{1},,{overrightarrow{sigma }}_{2},,{overrightarrow{sigma }}_{3})), since the configurations *σ*_{1}, *σ*_{2}, and *σ*_{3} have chemical orderings that are identical (up to a translation vector) to that of the bicrystal beyond the chemical interaction range of the glide plane.

In summary, the parameterization of surrogate models that accurately describe the GSF energies in multicomponent alloys requires two separate cluster expansions. The first is a cluster expansion of the formation energies of orderings over the parent crystal structure. This cluster expansion is required to calculate the average configurational energy, *E*^{avg}, in Eq. (1). Methods to parameterize these models are well-established. A second cluster expansion is necessary to describe the short-range CRPF energy. A training dataset can be generated by first calculating the bi-crystal energies, (E(overrightarrow{r},,delta ,,overrightarrow{sigma })), for several symmetrically distinct orderings, *σ*, glide vectors, (overrightarrow{r}), and separation distances, *δ*. The CRPF energies for each of these configurations can then be calculated with Eq. (8). The resulting CRPF energies for a fixed chemical ordering *σ* then serve to train the adjustable parameters of Eq. (6), which can then be cluster expanded according to Eq. (7) to describe their dependence on the degree of chemical order. In most Peirls-Nabarro-type models, the bi-crystal is assumed to be under zero stress in the direction perpendicular to the fault. As a result the GSF energy must be minimized relative to the slab separation distance. Thus, the energy given by ({E}_{0}^{{rm{CRPF}}}) in Eqs. (6, 7) is the desired quantity when modeling dislocation properties with PN models that assume zero tractions perpendicular to the glide plane.

### GSF energies in Mo–Nb

In this section, we illustrate the above cluster expansion formalism by investigating the composition and configuration dependence of USF energies in the binary Mo–Nb alloy. The Mo–Nb alloy adopts the bcc crystal structure at all compositions and forms a disordered solid solution between room temperature and the melt^{34,35,36}. The Mo–Nb binary is of current interest since Mo and Nb are both components of important bcc based high entropy alloys^{37,38}. Studies of the mechanical properties of Mo–Nb alloys indicate that their strength is primarily controlled by the formation and motion of screw dislocations and to some extent edge dislocations^{39}. The screw dislocations are formed along the 〈111〉 direction, and are known to spread on the {110} planes^{23}. Much of the physics of dislocations in Mo–Nb alloys is therefore directly related to the USF energy in the {110} plane for a relative displacement along the 〈111〉 direction. The USF energy serves as a key input to model the motion and evolution of screw dislocations in bcc alloys with phase-field dislocation dynamics^{40} and classical Peirls–Nabarro type models^{14,23,41}.

The first step in describing the dependence of the GSF energy on ordering is to construct a cluster expansion for the formation energy of the binary bcc Mo–Nb alloy. Figure 5a shows the formation energies of 847 symmetrically-distinct orderings on the bcc crystal structure in the binary Mo–Nb alloy as calculated with density functional theory (DFT). The formation energies are referenced to bcc Mo and Nb at 0 K. More details about the DFT calculations and the cluster expansion that was subsequently trained to these energies can be found in “Methods” section. The convex hull is outlined in black and shows that several ordered phases are predicted to be stable at 0 K. The energies of the 847 orderings as predicted with a cluster expansion are also shown in Fig. 5a as circles. The exceptionally low training error of 0.0008 eV/atom and the excellent qualitative agreement between the ground states as predicted with the cluster expansion and those found with DFT suggest that the configurational energy of the Mo–Nb binary alloy is well-described with a truncated cluster expansion model. Figure 5b also shows the relaxed volume of all orderings relative to that of bcc molybdenum. We find them to vary almost linearly as a function of niobium composition, albeit with a slight negative deviation in close agreement with Vegards law^{42} (shown schematically by the dashed line in the figure).

The GSF energy surface of a Mo bcc bicrystal for the {110} glide plane along the 〈111〉 direction is shown in Fig. 6a. The introduction of a planar fault due to a glide results in an energy penalty. The fault energy increases until it reaches a maximum at a glide of 1/4〈111〉. This energy corresponds to the USF energy. As the glide vector approaches a full lattice translation in the {110} glide plane (corresponding to 1/2〈111〉) the energy decreases until long-range bcc order is restored, where the energy becomes equal to that of bcc Mo. In calculating the GSF as a function of the glide vector (overrightarrow{r}) of Fig. 6a, we first calculated the energy of the bicrystal as a function of *δ* along the [110] direction to generate decohesion curves for each value of (overrightarrow{r}). A particular example of such a curve is shown in Fig. 6b. The DFT energies as a function of *δ* (for fixed (overrightarrow{r})) were then fit to the UBER^{30}, and the minimum of each curve was used to construct the GSF of Fig. 6a.

A section of the GSF energy for the B2 ordering with a composition of *x*_{Nb} = 0.5 is shown in Fig. 6c. Shifting the two halves of a B2 bicrystal through a full lattice translation results in an anti-phase boundary, which is accompanied by an increase of the energy. Similar to pure molybdenum, an USF is found to exist for a glide corresponding to half a translation vector. Figure 6d collects the calculated USF energies for pure Mo, L2_{1} (Mo_{3}Nb, MoNb_{3}), B2 (MoNb), and pure Nb. We find a strong composition dependence of the USF energies, with the values changing by almost a factor of two with increasing niobium composition. Furthermore, the USF energies vary non-linearly with composition, suggesting that short and long-range order also plays a role in addition to the average concentration. Figure 6d shows the energy at equilibrium separation for the sheared bicrystal.

The CRPF energies of 514 symmetrically distinct USF as calculated with DFT is shown in Fig. 7. The CRPF values vary from ≈0.8 J/m^{2} for pure niobium to about 2.0 J/m^{2} in the binary alloy. The spectrum of values at a particular composition is also found to span a large range of values, suggesting that the state of order among Mo and Nb plays a significant role in determining the CRPF energies.

A cluster expansion was parameterized to describe the dependence of the CRPF energies on the degree of Mo–Nb ordering. A comparison of the DFT and cluster expanded CRPF energies is shown in Fig. 8a. Details about the fitting procedure and cluster expansion model are provided in the “Methods” section. The USF energies are reproduced well by the cluster expansion model with a fitting error of 0.016 eV per two-dimensional unit cell of the (110) glide plane. The CRPF energies of configurations that have compositions close to pure molybdenum or niobium have a slightly higher error than configurations with compositions closer to *x* = 1/2. We validated the model by comparing cluster expansion predictions to DFT values of CRPF energies for 38 stochastically enumerated orderings in a 16 atom supercell. Figure 8b shows a good agreement between the model predictions and the DFT calculations with a validation error of 0.013 eV per unit cell. The interactions within the CRPF cluster expansion are relatively short-range as indicated by the sharply decaying magnitude of the point correlations shown in Fig. 9.

Having fit a cluster expansion that accurately describes the USF energy in the binary Mo–Nb system, we next investigated the composition and temperature dependence of this energy. Grand-canonical Monte-Carlo simulations at temperatures above 600 K are found to be completely disordered at all compositions, in agreement with experiment^{34}. Snapshots of disordered configurations were collected from grand-canonical Monte-Carlo simulations at 600 and 1000 K. Chemical potentials were chosen such that the average composition of niobium was 0.25, 0.5, or 0.75. For each Monte-Carlo snapshot, an USF was introduced in the cell and the USF energy was evaluated with Eq. (1) using the cluster expansions for the CRPF and the formation energy of the bcc Mo–Nb alloy. An USF energy was calculated by introducing a fault in every (110) layer within the simulation cell and for every (langle 1overline{1}1rangle) direction within the plane.

Figure 10a shows a histogram of USF energies at three different niobium compositions at a temperature of 1000 K. The USF energies decrease with increasing niobium composition. This is in agreement with the general trend of CRPF energies across compositions in Fig. 7. Our results predict that the stacking fault energies vary strongly with the average composition of the alloy. The spectrum of USF energies at elevated temperatures are normally distributed. The distribution is very sharply peaked at a niobium composition of 0.25 with a slightly more broadened distribution with increasing niobium composition.

We investigate the temperature dependence of the USF energies in Fig. 10b. Grand-canonical Monte-Carlo calculations were performed at temperatures of 600 K, and 1000K at chemical potentials that corresponded to average alloy compositions of 0.5. Figure 10b shows that the distribution of USF energies continues to be normally distributed across a wide range of temperatures. The magnitude of the fault energy decreases with increasing temperature. The distribution of energies is also broadened with increasing temperature. These results suggest that the average USF energies vary strongly with niobium composition, while the distribution of energy values is sensitive to the temperature.