# Atomic Representations of Local and Global Chemistry in Complex Alloys

M. J. McCarthy<EMAIL_ADDRESS>Center for Computing Research, Sandia
National Laboratories, Albuquerque, New Mexico 87185, USA  J. Startt
<EMAIL_ADDRESS>Center for Integrated Nanotechnologies, Sandia National
Laboratories, Albuquerque, New Mexico 87185, USA  R. Dingreville
<EMAIL_ADDRESS>Center for Integrated Nanotechnologies, Sandia National
Laboratories, Albuquerque, New Mexico 87185, USA  A. P. Thompson
<EMAIL_ADDRESS>Center for Computing Research, Sandia National
Laboratories, Albuquerque, New Mexico 87185, USA  M. A. Wood
<EMAIL_ADDRESS>Center for Computing Research, Sandia National
Laboratories, Albuquerque, New Mexico 87185, USA

###### Abstract

The exceptional properties observed in complex concentrated alloys (CCAs)
arise from the interplay between crystalline order and chemical disorder at
the atomic scale, complicating a unique determination of properties. In
contrast to conventional alloys, CCA properties emerge as distributions due to
varying local chemical environments and the specific scale of measurement.
Currently there are few ways to quantitatively define, track, and compare
‘local’ alloy compositions (versus a ’global’ label, i.e. equiatomic)
contained in a CCA. Molecular dynamics is used here to build descriptive
metrics that connect a global alloy composition to the diverse local alloy
compositions that define it. A machine-learned interatomic potential for
MoNbTaTi is developed and we use these metrics to investigate how property
distributions change with excursions in global-local composition space. Short-
range order is examined through the lens of local chemistry for the equiatomic
composition, demonstrating stark changes in vacancy formation energy with
local chemistry evolution.

## INTRODUCTION

The basic metallurgical idea behind complex concentrated alloys (CCAs) is two-
fold: the maximization of the solid solution strengthening effect and phase
stabilization via configurational entropy versus enthalpy. Composed of three
or more metallic elements, CCAs represent a family of alloys with a large
composition space for property improvement to be explored. Many physical
metallurgy research groups have explored this space and ventured away from the
reference equiatomic composition to seek novel compositions that yield
improved properties ranging from high temperature strength 1, 2, 3, 4, 5, 6,
to wear resistance 7, 8, to corrosion and oxidation resistance 9, 10, 11, to
magnetic properties 12, 13, 14. For instance, Startt and coworkers 6
investigated how incremental changes in elemental compositions affected
thermomechanical properties in the MoNbTaTi refractory CCA, finding that
increases in the concentration of Mo lead to significantly improved yield
strength, while higher Ti concentrations lead to greater ductility. Similarly,
Zhao and colleagues 12 showed that increasing the amount of Co in a
CoCrCuFeMnNi alloy resulted in a higher magnetization saturation moment and an
improved soft magnetic response.

While most studies use only the global alloy composition, $c_{\rm global}$, to
characterize and describe an alloy and its properties, the reliance on a
single global descriptor can be problematic. Though they exist in ordered
crystalline states, CCAs typically possess a large degree of intrinsic
chemical disorder, leading to large fluctuations in micro- and nanoscale
properties. Those properties can in turn be highly sensitive to local
environment and short-range ordering effects. Such local/global duality adds a
new dimension and challenges to the atomic representation of the huge
compositional space in these multicomponent alloys. Atomistic modeling
techniques such as density functional theory (DFT) and molecular dynamics (MD)
offer ways to quantify the connections between local and global chemistry in
these complex alloys. Recent developments in the training and validation of
machine-learned interatomic potentials (ML-IAPs) have made CCAs more
accessible in MD, allowing for deeper investigation into local and global
property relationships. Several ML-IAPs have already been developed to study
MoNbTaW-based quaternary and quinary refractory BCC alloy systems and have
been used to study various CCA properties such as segregation and defect
formation 15, strengthening mechanisms 16, and dislocation mobility 17. In all
cases, both the ML-IAP training procedures as well as subsequent analysis were
centered around global composition descriptions, potentially obscuring
important trends that may arise from local chemical fluctuations. These
descriptive limitations arise from the lack of a unified language or system to
describe variability in chemistry.

In this work, we are focused on the intricate relationships between the local
and global chemistry in complex alloys. We are particularly interested in the
implications these relationships have in estimating alloy properties,
validating CCA interatomic potentials’ performance through composition space,
and in training improved CCA ML-IAPs in the future. To do this, we begin by
defining the composition space of a quaternary CCA (elements labeled A, B, C,
and D) in Fig. 1, where the concept of local and global duality is
illustrated. In Fig. 1a, the tetrahedron represents the entire composition
space spanned for our quaternary alloy. The inner volume of this tetrahedron
contains all possible quaternary compositions and the center of the
tetrahedron represents the equiatomic composition (black star). The vertices,
edges, and faces of the tetrahedron correspond to pure, binary, and ternary
compositions, respectively. From a macroscopic perspective a global
composition, $c_{\rm global}$, corresponds to a single point within this
tetrahedron for an infinitely large random solid solution structure, such as
that in Fig. 1b. Conversely, Fig. 1c illustrates the concept of a _local
composition_ , $c_{\rm local}$, which can be defined as the concentration of
elements found within a spherical volume element of radius $r_{\rm cut}$
centered on any one atom or lattice site, $i$. Mapping all local compositions
for a given $r_{\rm cut}$ from all atoms in a structure yields a point cloud
of $c_{\rm local}$ values centered around the parent structure’s $c_{\rm
global}$, similar to how averaging all atoms’ local composition yields $c_{\rm
global}$. Varying $r_{\rm cut}$ changes the number of atoms locally sampled
and thus the density of the point cloud, as shown in the two tetrahedra of
Fig. 1c. The variability in $c_{\rm local}$ depends on $r_{\rm cut}$ and will
eventually dwindle and converge to $c_{\rm global}$ as the cutoff increases to
infinity.

Figure 1:  Illustration of global and local composition spaces for a
quaternary alloy of generic elements A B C D. a Global composition
perspective. The location of a quaternary alloy’s equiatomic composition in
the phase tetrahedron is shown with a star. The vertices, edges, and faces of
the tetrahedron correspond to pure, binary, and ternary compositions,
respectively. Dotted black line tracks the changes in one element’s
concentration from pure A at the vertex, to a CCA enriched in A (square), to
one depleted of A (triangle) to a ternary alloy consisting of B, C, and D.
Gray dotted line starting at element D is included as a guide for the eye. b
An example equiatomic random solid solution structure with elements A, B, C,
and D. The mapping of local compositions $c_{\rm local}$ for each atom $i$ for
different values of $r_{\rm cut}$ is shown in panel c. c Local composition
perspective. An atom $i$’s local composition is calculated by taking the
concentration of all elements found within a sampling radius $r_{\rm cut}$
around said atom. Different values of $r_{\rm cut}$, here 4 Å (top) vs. 6 Å
(bottom), capture different numbers of neighboring atoms. Mapping all of a
structure’s atom’s local compositions results in point cloud sizes that depend
on $r_{\rm cut}$.

The relationship between these two perspectives can be quantified by three
characteristic metrics: a composition deviation, $\lambda$, which measures the
distance between two compositions (either global and/or local); a volume
fraction, $V_{\rm f}$, that quantifies the amount of composition space a point
cloud spans in the local composition tetrahedron; and the cutoff radius,
$r_{\rm cut}$, which defines a physical scale of interest as well as the
chemical resolution tying both perspectives together. The equations
corresponding to each metric are included in the Methods below.

In what follows, we demonstrate how $\lambda$ and $V_{\rm f}$ can be used to
draw connections between local and global compositions, and how material
properties in CCAs can be affected by deviations in composition at the local
scale. We start by building an ML-IAP for a four-component CCA using a
training database generated from ab-initio calculations. The MoNbTaTi alloy is
a BCC refractory CCA that has been shown to exhibit not just excellent
strength and elastic stiffness but is also prone to significant deviations in
strength when incremental changes are made to elemental compositions 6. Thus,
it is hypothesized that local composition fluctuations could prove to have
significant effects on material properties and behaviors, making this CCA an
ideal choice for this present study. We then use this trained ML-IAP to
illustrate how to sample a representative local chemistry and the implications
of this sampling on scale-dependent property measurement across composition
space. Finally, we examine how local-global composition sampling changes when
incorporating CCA chemical effects such as short-range order.

## RESULTS

### Multi-compositional refractory CCA ML-IAP for MoNbTaTi

In this section, we describe the fitting process and the performance of the
MoNbTaTi CCA ML-IAP used in the following sections for analysis. The primary
target in fitting this potential was to ensure accuracy across a wide range of
global composition space, with a particular focus on replicating DFT measured
elastic properties. Because elasticity is a property measured over an entire
simulation cell, and not locally, we did not apply the local composition
analysis techniques outlined above to optimize the fitting process. Instead,
we followed the global composition sampling scheme detailed in the Methods.
For creating fitting data, the global composition was systematically varied
via a pseudo-equiatomic composition sampling scheme by adjusting the
concentration of one element, labeled element A, against the other three
elements, labeled elements B, C, and D.

This sampling scheme follows that of the black dotted line in Fig. 1a,
illustrating that all global composition training points fall along single
lines connecting each vertex (corresponding to a pure element) to its opposing
tetrahedral face (corresponding to an equiatomic ternary alloy completely
absent of the vertex element). Each DFT training structure was averaged over
four special quasirandom structures (SQS) to approximate bulk mixing
character. Starting from the equiatomic composition in the center, each
component element A was depleted to 12.5% at minimum and enriched to 50.0% at
maximum.

Composition | Property |  DFT 6 | SNAP | Error (%)
---|---|---|---|---
Equi. |  $\mathbb{C}_{11}$ [GPa] | 239.6 | 255.3 | 6.6%
$\mathbb{C}_{12}$ [GPa] | 129.7 | 139.7 | 7.7%
$\mathbb{C}_{44}$ [GPa] | 37.8 | 39.1 | 3.4%
$B$ [GPa] | 165.7 | 178.4 | 7.7%
$G$ [GPa] | 43.6 | 42.7 | 1.0%
$E$ [GPa] | 120.1 | 126.4 | 5.2%
| $\nu$ | 0.379 | 0.382 | 0.8%
Mo 12.5 at.-% |  $\mathbb{C}_{11}$ [GPa] | 234.3 | 233.8 | 0.2%
$\mathbb{C}_{12}$ [GPa] | 134.6 | 134.3 | 0.2%
$\mathbb{C}_{44}$ [GPa] | 35.6 | 36.0 | 1.0%
$B$ [GPa] | 167.8 | 167.5 | 0.2%
$G$ [GPa] | 40.8 | 40.1 | 0.5%
$E$ [GPa] | 113.2 | 113.6 | 0.4%
| $\nu$ | 0.388 | 0.387 | 0.2%
Mo 50 at.-% |  $\mathbb{C}_{11}$ [GPa] | 334.4 | 312.5 | 6.5%
$\mathbb{C}_{12}$ [GPa] | 147.3 | 149.3 | 1.3%
$\mathbb{C}_{44}$ [GPa] | 50.331 | 50.0 | 0.6%
$B$ [GPa] | 209.7 | 203.7 | 2.9%
$G$ [GPa] | 64.7 | 60.9 | 5.8%
$E$ [GPa] | 175.9 | 166.2 | 5.5%
| $\nu$ | 0.360 | 0.364 | 1.1%
Table 1:  Elasticity value comparisons from the MoNbTaTi dataset and the SNAP
model for the equiatomic random solid solution and an example at Mo enriched
to 50 at.%. Each value was averaged over four DFT special quasirandom
structure calculations and five SNAP calculations in MD.

For the purposes of this work, we selected one well-performing ML-IAP
generated using the above parameters, though many thousand ML-IAPs were
generated and tested. A schematic of the fitting process is shown in Fig. 2a,
consisting of two major parts: a linear regression which connects the DFT data
to LAMMPS bispectrum calculations, run in the software FitSNAP 18 and a
genetic algorithm to optimize SNAP hyperparameters, run in the DAKOTA software
19. More details on each part, and the specifics of fitting to multiple
compositions, are included in the Methods. Errors on the energies and forces
for this ML-IAP are 12 meV and 176 meV/Å, respectively. Representative model
errors for elastic constants and moduli for the equiatomic and two Mo
compositions (Mo at 12.5 at.% and 50 at.%) are shown in Table. 1. All errors
for the equiatomic and for a great number of other fitted compositions remain
in general below 10% of their respective DFT values and frequently near or
below 5%. Fig. 2b plots the bulk modulus B and shear modulus G for a selection
of trained compositions for both the trained SNAP ML-IAP (curves) and DFT
values (symbols). Values from the trained model’s bulk and shear moduli ($B$
and $G$, respectively) are shown in Fig. 2b for 29 selected compositions. Data
points from the DFT training (symbols) are averaged over four SQS elasticity
calculations (taken from Startt et al. 6), and trained SNAP ML-IAP (curves)
are averaged over five MD calculations. The error bars of the SNAP curves show
the standard deviation over five modulus calculations. Note that, because the
DFT and SNAP values for the equiatomic calculations had $<1$% error, the
symbols have been left off for clarity. The horizontal axis of these plots
follows the concentration of element A, whose species assignment is indicated
by colors and symbols. In general, the ML-IAP moduli match closely to the
expected DFT values for the majority of the trained global alloy compositions.
Depletion or enrichment of Mo, Nb, and Ti from the equiatomic composition are
especially well-fit. The largest errors are found in alloys with Ta at the
highest (approx. 20% error at 50 at-% Ta) and lowest concentrations (approx.
8% error at 12.5 at-% Ta).

Figure 2:  a Schematic of the ML-IAP fitting process. The FitSNAP software 18
performs regression while LAMMPS is used to calculate descriptors. The
resulting fit is run through a series of tests, whose errors are used to
inform a single-objective genetic algorithm provided by DAKOTA. By iterating
thousands of times, the genetic algorithm is able to isolate ML-IAPs with
optimal material properties and excellent dynamic stability. b Calculations of
the bulk modulus B (top) and shear modulus G (bottom) for selected
compositions using the SNAP MoNbTaTi model. Solid lines follow the SNAP-
calculated values for enrichment or depletion of each element in the
refractory CCA where corresponding shading indicates the standard deviation.
The symbols show averaged data from the DFT training, taken from Startt et al.
6.

### Representative Sampling of Chemistry

Given the availability of an IAP, one of the key advantages of MD simulations
is their ability to scale. By dint of its larger possible size, a single MD
simulation cell ($N_{\rm MD}\leq 10^{10},\leq 1~{}{\rm\mu m}^{3}$) can sample
a vastly wider range of local compositions than even a large ensemble of DFT-
sized ones ($N_{\rm DFT}\leq 10^{4},\leq 10~{}{\rm\AA}^{3}$). This range means
that, for a given global composition $c_{\rm global}$ of a random solid
solution (i.e. MoxNbyTazTi(1-x-y-z)), a large fraction of $N_{\rm MD}$ will
contain many high probability local compositions ($c_{\rm local}$) near
$c_{\rm global}$, as well as extreme values of $c_{\rm local}$ that deviate
from $c_{\rm global}$. The first effort herein is focused on quantitatively
defining these length scales where $N_{\rm MD}$ can be considered as a
representative volume of some average chemistry, $c_{\rm global}$. Caveats to
consider are then: would it suffice to examine only the most common
environments? What represents an uncommon environment? And how different is it
from the most common ones? In this section, we illustrate how metrics of
distance and volume within composition space, represented by $\lambda$ and
$V_{\rm f}$, can be used to answer these and related questions.

The first task is to understand the distribution of $c_{\rm local}$ around a
single $c_{\rm global}$. We define the scalar distance metric $\lambda=|c_{\rm
global}$ \- $c_{\rm local}|$, taken element wise, that maps any arbitrary
local composition’s deviation from the global one. Its maximum value is
$\sqrt{2}$ (the distance between any two pure elements using the Euclidian
norm), and its minimum value is determined by $r_{\rm cut}$. The cumulative
distribution function (CDF) of $\lambda$ for three different $r_{\rm cut}$
values is shown in Fig. 3a. Reflective of the qualitative local perspective
data displayed in Fig. 1c, it is clear in Fig. 3a that the median $\lambda$ is
inversely proportional to $r_{\rm cut}$. This makes intuitive sense due to the
fact that when fewer atoms are collected in calculating $c_{\rm local}$, the
concentration will change more drastically with a single atom swap.
Additionally, these single atom composition changes also explain the more
jagged nature of the CDF for smaller $r_{\rm cut}$.

Figure 3:  a CDF of $\lambda$ for three values of $r_{\rm cut}$. b Duplicate
local compositions $\delta$ versus the deviation from the global composition
$\lambda$ for each $r_{\rm cut}$. Inset is of s = 40 and $r_{\rm cut}$ = 6 Å ,
with the color representing the number of duplicates $\delta$. c Volume
fraction $V_{\rm f}$ as a function of simulation cell size, s. Dashed lines
show $V_{\rm f}$ calculated with all values of $\lambda$ (as in the top inset
of panel a), solid lines where 0.5% of the highest values of $\lambda$
excluded. Excluding 0.5% of $\lambda$’s largest values not only significantly
lowers $V_{\rm f}$, but also stabilizes it to a constant value for all $r_{\rm
cut}$ shown. d Volume fraction $V_{\rm f}=V_{\rm f}^{99.5\%}$ for a range of
compositions included in the training set (see Fig. 1a for reference), using
$s=40a_{0}$ and $r_{\rm cut}$ = 6 Å averaged over 5 runs for one example
element. Insets show the convex hull geometry of $r_{\rm cut}$ = 6 Å with all
$\lambda$ included (top), and with the largest 0.5% $\lambda$ excluded
(bottom).

From the discreteness in single atom changes to $c_{\rm local}$ and on-lattice
measurement of $\lambda$ comes another metric, the duplicate local
compositions, $\delta$, which is a count of the number of atomic environments
that have identical local composition values. Fig. 3b shows $\delta$ as a
function of $\lambda$ for different values of $r_{\rm cut}$. Projected into
the composition tetrahedron, the inset shows local composition data for
$r_{\rm cut}$ = 6 Å with a log-scaled color bar. For $r_{\rm cut}$ = 4 Å, the
smallest possible compositional change is $\sim$5.77 at-%, where this
definition of length scale corresponds to the two nearest-neighbor shells in
BCC. The smaller cutoff radius results in higher degeneracy of observed
$c_{\rm local}$, and coarser $\lambda$ spacing. Though duplicate compositions
are chemically degenerate, they are likely to be structurally unique, just as
is true for global compositions as a whole. The differences between each
duplicate environments’ atomic arrangements gives rise to a distribution of
atomic energies (as predicted by the SNAP ML-IAP developed herein) per $c_{\rm
local}$.

Due to the coarse chemical representation of short cutoff distances, the
volume contained by a bounding surface of $c_{\rm local}$, $V_{\rm f}$, shown
in Fig. 3c, also increases. As the occurrence of extrema local chemistries is
more probable at short cutoff distances, these outlier points are filtered out
by a modest $0.5\%$ rejection of the highest $\lambda$ values (Fig. 3c, solid
lines). Compositional volume fraction $V_{\rm f}$ data are plotted against the
size of a cubic simulation cell, where the number of unit cell replicas, $s$,
spans atom counts capable of being represented in DFT, 54 ($s=3$), to those
only possible in MD, $\sim 1\cdot 10^{6}$ ($s\geq 70$). Sampling of $\lambda$
is taken from 5 unique random solid solution cubes per data point with unit
side length of $s$. Shaded bands around each line and point encompass the
standard deviation for each $r_{\rm cut}$.

At small cell sizes, before $V_{\rm f}$ has been saturated, not enough unique
$c_{\rm local}$ configurations will be present for the simulation cell to be
chemically representative of $c_{\rm global}$. The apparent lack of saturation
in the data for $r_{\rm cut}$ = 4 Å is mathematically explained by the
probability for $c_{\rm local}$ to be a single element species. Element labels
drawn at random for the 15 nearest atoms equal to the central atom type is
$p=(\frac{1}{4})^{15}$, which multiplied by the $10^{7}$ atom environments
sampled for $s=70$ is not a vanishingly small probability($\approx 0.93\%$).
The same random draw where $r_{\rm cut}$ = 6 Å can be safely assumed as a
complete chemically representative volume element. Filtering out the largest
0.5% of $\lambda$ values ($V_{\rm f}^{99.5\%}$) leads to a significant drop in
volume fraction, as shown in solid lines in Fig. 3c. Fig. 3d captures
saturated $V_{\rm f}^{99.5\%}$ values while moving along a tie line of one
element A (e.g., where A = Mo, then ${\rm Mo}_{\rm x}({\rm NbTaTi})_{\rm
1-x}$). There it can be seen that the chemically representative volume element
is more compact for compositions near a dilute single component or ternary
composition of the CCA. For example, the chemically representative volume
element for random solid solution equiatomic cubes sampled with $r_{\rm cut}$
= 6 Å is $V_{\rm f}^{99.5\%}=0.15$ located directly in the center of the phase
tetrahedron. In contrast, chemically representative volume elements for near-
ternary alloys (where the concentration of element A approaches $<$ 10 at.-%
of element A) or element A-enriched alloys (increasing at-% of element A) are
comparatively smaller than the equiatomic. The same measure taken for a
$c_{\rm global}$ near a single element vertex would result in a comparatively
small $V_{\rm f}^{99.5\%}$, given the fewer number of possible local chemical
arrangements. It is also important to note that, though we have chosen to
filter the largest 0.5% values of $\lambda$, one could also choose heavier
filtering, which would decrease $V_{\rm f}$ proportionally. As a visual guide,
the two inset images of Fig. 3d show two cases of the bounding surface formed
by each of the $c_{\rm local}$ values, either where the full CDF (top inset)
is used, resulting in a rough surface, or when the lowest 99.5% values of
$\lambda$ (bottom inset) is used, yielding a smooth surface. The exact choice
of a saturation filter is less important than having the ability to make
informed decisions about the degree of chemical accuracy desired for CCA
modeling and simulation.

Taking this quantitative measure one step further, the complexity of an alloy
can be defined by its $V_{\rm f}^{99.5\%}$, which has merit for alloys of both
fewer and more elements included in them. While the visualization of the
volume contained by the collection of $c_{\rm local}$ in five or more element
composition space is challenging, using these relatively simple measures of
complexity make comparisons possible. As will be explored in an upcoming
section, complexity metrics also enable an assessment of the degree of
chemical change within a single $c_{\rm global}$. For example, where short-
range order is present in a given CCA, $\lambda$, $\delta$ and $V_{\rm f}$
will reflect the unique evolution with respect to an equivalent random solid
solution at the same $c_{\rm global}$. Furthermore, when selecting a
computational method for predicting properties of CCAs, these chemically
representative volume element measures confirm that DFT-sized simulation
cells, which have maximal atom counts on the order of $10^{2}$ ($s\leq 5$) for
a single cell, sample far too small a fraction of the chemically
representative volume element (which includes all duplicates per local
composition $\delta$) to capture global CCA properties accurately.

Having determined a reasonable simulation cell size from which to properly
sample CCA properties (i.e., $s\simeq 40$), we now turn predictions of
material properties that can be reasonably assumed to be of the same length
scale as these $r_{\rm cut}$ definitions.

### Scale-Dependent Property Measurement

Material properties in CCAs are generally expressed as distributions rather
than as single, averaged values typical of chemically simpler systems. The
question then becomes whether this property distribution is explained by the
$c_{\rm local}$ distribution defined herein. This assertion is critically
dependent on the characteristic length scale of the property of interest, and
whether that length scale demonstrates the chemical variability discussed
previously. Example properties and their associated length scales could
include: elastic moduli that encompass the entire chemically representative
volume element, dislocations that sample as many $c_{\rm local}$ as atoms
along the line, or point defects that may (depending on the strain field
generated) only reflect a single $c_{\rm local}$. In order to enable high-
throughput alloy design in CCAs, we need an ability to characterize which
properties are sensitive to the distribution of $c_{\rm local}$, or can be
assigned to $c_{\rm global}$. A material defect whose properties are highly
spatially localized is the vacancy formation energy. It has been shown that
vacancy formation and migration energies in CCAs have varied distributions 20
which can for example affect the evolution of radiation damage cascades 21,
22.

Utilizing the SNAP ML-IAP described above and in the Methods, we sample the
vacancy formation energy for each element species, measuring $\lambda$ in
order to address how properties are distributed around a given $c_{\rm
global}$. The local chemical definition coincides with the interaction range
of the ML-IAP where $r_{\rm cut}$ = 6 Å. This distance sufficiently resolves
the strain field at a vacancy, where the outermost (see Fig. 1b) neighbor
shell shows local strain values less than $6\cdot 10^{-4}$ as measured by
tools made available in OVITO 23.

Fig. 4 shows the vacancy formation energy, $E_{\rm vac}$, of removed Nb atoms
versus $\lambda$, with colors indicating CCA compositions that are Nb-depleted
and Nb-enriched in panels a and b, respectively. Details on vacancy formation
energy calculations can be found in Methods. The contour plots indicate
regions of similar data point density, calculated using a kernel density
function. The outermost contour contains 90% of data points, with each line
thereafter delineating a 20% change in point density up to the innermost line,
which encompasses the final 1% of data. This smallest contour indicates the
most probable values from the $\lambda$ and $E_{\rm vac}$ distributions. The
univariate distribution plots provide another means of interpreting point
densities with respect to its parallel axis. We keep the equiatomic
composition included in both plots as the aforementioned baseline CCA
distribution example.

Figure 4:  a-b Vacancy energy trends vs. $\lambda$, the deviation from a
chosen global composition. Each plot features distributions for Nb vacancies.
Panel a illustrates trends while depleting Nb and panel b shows those for
enriching Nb in the CCA, with each global composition indicated by different
colors. Contour lines, calculated using a kernel density estimate, enclose
(outer to inner) 90%, 70%, 50%, 30%, and 10% of data points, and the innermost
circle the most probable 1%. The dotted red lines show the average vacancy
energy in a pure Nb sample, which corresponds to $\lambda=0$. The X- and
Y-axes show the univariate distribution for the vacancy energies and
deviations from global, respectively.

In Fig. 4a, the Nb depleted compositions, which moves toward the MoTiTa
ternary face in the composition tetrahedron, show very little change in this
localized property with respect to the composition deviation and maintain a
$E_{\rm vac}$ as seen in the average equiatomic CCA (= 3.40 eV) than pure Nb
(= 2.76 eV). The converse is true for the enriched side of the tie line. The
composition deviations gradually decrease in tandem with the property
measurement, converging on the pure element’s value.

Given the locality (as determined by the rapidly diminishing strain field) of
a vacancy, it is reasonable to attribute the formation energy to the structure
change within the range of the $c_{\rm local}$ definition. Therefore, it is
incorrect to assign $E_{\rm vac}$ to $c_{\rm global}$. Rather, each $c_{\rm
local}$ carries its’ own property value and the global composition only
determines the range of compositions where that property is averaged based on
the chemically representative volume element. If a larger $r_{\rm cut}$ were
used to define the property measurement (and thus the local composition) it
would have the effect in a random solid solution to concentrate $c_{\rm
local}$ to smaller values of $\lambda$. But this effect is subverted by
allowing for short-range order to develop. In a chemically-ordered CCA
structure, more exotic local environments will be expressed in comparison to a
random solid solution. As will be shown in the following section, these
environments give even more merit to assigning properties to local chemistries
than global averaged values.

### Local composition and short-range order

As useful as a _truly random_ solid solution is for general modeling purposes,
it is unlikely to accurately reflect the equilibrium state of real CCAs. Far
more probable is the development of at least some degree of short-range order
24, 25, 26. Though small in CCAs by design, finite differences in chemical
potential between species imply that some reordering is highly probable,
especially given factors such as intrinsic defects and thermo-mechanical
processing steps that can activate such reordering. In terms of the metrics
defined above, even mild degrees of short-range order evolution will
significantly change the geometry of the local composition point cloud. The
values of $\lambda$ that are chemically representative, $V_{f}$, $\delta$, and
what constitutes an outlier will shift in turn. Critically, these all will
change even for a single global composition label. This means that, for CCAs,
global composition not only obscures important scale-dependent property
measurements for the random solid solution, but also lack descriptive power in
capturing shifts in properties due to chemical interactions. In this final
section, we will use the above-defined metrics to explore how short-range
order alters both local composition space and CCA properties. This portion of
the study is critically dependent on the availability of an accurate energy
model, even for extrapolative local chemistries.

For the MoNbTaTi system specifically, hybrid Monte Carlo MD (MC/MD) and DFT
studies of CCAs containing Mo, Ta, and Nb show that Mo and Ta strongly favor
each other as neighbors and can form B2 phases 27, 16, 15, 28. Hybrid MC/MD
simulations of the SNAP MoNbTaTi featured in this work show exactly the same
trend, which is clearly reflected in changes in the local composition point
cloud as compared to an initially randomly-ordered sample. Details of the
setup of the hybrid MC/MD simulations are included in the Methods. Fig. 5
demonstrates how chemical reordering in an equiatomic sample is reflected in
local composition space through time, using $r_{\rm cut}$ = 6 Å and a
simulation cell with 128,000 atoms ($s=40$). Initially, the structure begins
as a random solid solution (Fig. 5a), reflected in the spherical shape of the
local composition point cloud. Note that the colors of the points reflect the
species of the central atom being sampled. Fig. 5b tracks the evolution of
both $V_{f}^{99.5\%}$ and the cell’s Warren-Cowley parameters 29 for the first
nearest-neighbor shell $\alpha_{1}$. Near-zero values of alpha indicate random
ordering, negative values indicate species pair attraction, and positive
values indicate repulsion. Further details of the calculation of $\alpha_{1}$
are described in the Methods. By the final simulation step at $t=150$ ps,
shown in Fig. 5c, the structure has become clearly ordered. This change in
local chemistry is reflected in the point cloud geometry, which has narrowed
and spread into new local composition regions. The volume fraction has
significantly increased, from an initial $V_{f}^{99.5\%}=0.15$ to approx. 0.41
in the chemically-ordered cell. Especially notable is a shift towards the
binary Mo-Ta phase tie line, demonstrating these elements’ capacity to form B2
phases in CCAs.

The favored swapping towards Mo-Ta pairings deplete equiatomic local
compositions and shift those points into Ti- and Nb-richer regions of the
tetrahedron. Both of these trends correspond to the pairs with the largest
negative values of $\alpha_{1}$ at $t=150$ ps, namely Mo-Ta, Ti-Ti, and Nb-Nb.
Supplemental Fig. 1 shows the evolution of average potential energy through
time, as well as a visualization of B2 ordering.

The structural changes arising from short-range ordering will leave their mark
on the distribution of local chemical properties relative to the random solid
solution. Fig. 6 illustrates how evolved short-range order alters the vacancy
formation energy distributions of the initially-random equiatomic cell. In
both panels, dashed lines indicate the initial states of vacancy energies, and
solid lines show the result after reordering. All elements except Ti exhibit
significant shifts in vacancy formation energies, as shown in Fig. 6a. Mo
atoms undergo clear changes in the shape of the property distribution, with a
total cumulative change of -0.28 eV especially represented in the lowest 50%
of the curve. Ta undergoes an overall increase of 0.61 eV in formation
energies. All Nb vacancy formation energies decrease in the ordered state by
-0.50 eV. The highest vacancy energies of Ti increase by 0.08 eV, the smallest
change of all elements. Negative (favorable) vacancy formation energies when
removing Ti are reflective of the metastability of Ti in a BCC phase. The
unusual shift in the energy distributions of Mo vacancies is explored in more
detail in Fig. 6b, with respect to the composition deviation $\lambda$. While
there is some overlap in the contours of the random and ordered structures,
the energy values have split into two distinct regions at larger values of
$\lambda$, one centered near $E_{vac}\approx 3.25eV$ and the other much lower
at $E_{vac}\approx 1.25eV$. An analysis of the characteristic local
compositions at $\lambda>0.3$ indicates that Mo vacancies in the lower energy
region have local compositions especially rich in Ti-Nb. The higher-energy
region, which generally also overlaps the random structure’s energies, is
characterized by Mo-Ta enriched local environments.

Though further detailed analysis of trends in $E_{vac}$ is beyond the scope of
the current work, the additional information given by $\lambda$ enables
straightforward segmentation of trends in atomic-scale properties. These
results illustrate the potential power of local composition analysis,
highlighting other major benefits of MD studies and ML-IAPs for CCAs. For
example, the local composition analysis techniques developed here can be used
to uncover the degree of short-range order change and shift in chemical space
from the initial $c_{\rm global}$ not only from two single points in time, but
also throughout an entire MC/MD trajectory. Such information could be used to
correlate degrees of short-range ordering with e.g. information gathered from
advanced microscopy data 30, creating a potential bridge between the
timescales of experimental and computational ordering measures. The above
examples are viewed as avenues for future development of these alloy chemical
complexity metrics.

Figure 5:  Evolution of the local composition tetrahedra ($r_{\rm cut}$ = 6 Å
) and physical samples of a 128,000-atom equiatomic simulation of hybrid Monte
Carlo/molecular dynamics at $T_{MD}$ = $T_{MC}$ = 300 K for $t=150$ ps
(768,000 MC atom swap attempts). The colors of the points in the tetrahedra
correspond to the species of the sampled central atom. a The state of the
composition tetrahedron at $t=0$ ps. The local composition point cloud is
spherical, reflecting the initially random atomic ordering. b Change of the
volume fraction $V_{f}^{99.5\%}$ (top panel) and the Warren-Cowley parameter
for the first nearest-neighbor shell (bottom panel) through simulation time.
Negative numbers indicate pair attraction and positive indicate pair
repulsion. The most attractive pair by this criteria is Mo-Ta, followed by Ti-
Ti and Nb-Nb. c The initially-spherical composition point cloud spreads out
into broader regions of the local composition tetrahedron, saturating parts of
the binary Mo-Ta phase tie line, centered around 50 at.% and expanding into
CCA compositions roughly spanning the entire Ti-Nb tie line.  Figure 6:  a
Cumulative distribution functions of the vacancy formation energies $E_{\rm
vac}$ for each element in the random solid solution (RSS, dashed curves) and
short-range ordered (SRO, solid curves) cells from Fig. 5, panels a and c
respectively. The shaded regions indicate the shift in energies due to
chemical ordering. b $E_{\rm vac}$ and $\lambda$ for Mo vacancies from panel
a. The contour outlines share the same parameters as those from Fig. 4. Two
regions of different average energies can be observed at values of
$\lambda\gtrsim 0.3$, which can be further analyzed to isolate local
environmental trends. The region with distinctly lower energies ($\approx
1.25eV$) is correlated with local compositions rich in Ti-Nb.

## DISCUSSION

While the details of the model development were not the focus of this work,
generating an interatomic potential for a CCA has brought about unique
perspective on measuring material properties in CCAs. What makes classical
molecular dynamics a powerful tool to study materials is the computational
efficiency derived from the key approximation that interactions between atoms
are localized within some distance. In turn, we have addressed how this
locality has derived alloy complexity measures such as the chemically
representative volume element and volume occupied in composition space for a
given global composition. Additional scrutiny is now required on our
computational methods when predicting properties of CCAs due to these measures
of chemical complexity and locality of a property, which we believe will cause
a dramatic shift in this field.

Looking forward, data of $c_{\rm local}$ and chemically representative volume
elements of CCAs should be universally used in training ML-IAPs themselves to
ensure proper chemical fidelity. Clear metrics of chemical complexity allow
confirmation that targeted regions of composition space have been adequately
represented in training, and enables the automated generation of new training
structures where alloy chemistry data is lacking. Looking beyond model
validation, these techniques could also open new possibilities in scale
bridging. In principle, one could compute material properties in a pre-defined
grid in composition space. When querying a new alloy of a desired (or
measured) $c_{\rm global}$ a property can then be inferred based on averaging
nearby $c_{\rm local}$ and observing whether those points are contained within
the surface of a chemically representative volume element. This type of
prediction would allow MD results to inform property variability at much
larger (micro-, meso-) scales as well because chemical inhomogeneities can be
defined at numerous length scales. For example, if experiments or CALPHAD
indicate that certain phases or regions of global composition should be
present in CCAs, the chemically representative volume elements for those
phases can be added to the training process to ensure models are sampling
chemical environments (and also material properties) adequately. This kind of
diverse-by-design chemical training approach could be expanded to allow CCA
models to capture important phenomena not only through composition space, but
also temperature and pressure. Such capabilities would greatly expand MD
simulation utility in becoming truly integrated into workflows of alloy design
and discovery.

All levels of property analysis involve choices about length scale - indeed,
navigating those choices is perhaps the most fundamental challenge of
materials science. Thus, what is needed to bring modeling and simulation
efforts in closer concert with experimental capability is a review of the
definition of critical length scales. When practitioners of ab initio methods
identify a material with its global composition, it should be cast as one of
many local composition measurements in MD. The same can be said for global
composition in MD - it is a very local measurement when compared to even the
most resolved experimental measurements in, for example, atom probe
tomography. If we are to achieve high-throughput, or even a baseline
capability to screen for optimal CCA compositions, these ideas of chemical
representation of length scales need to be addressed more completely. The
results herein bring quantitative measures to this challenge, and address how
scale-bridging efforts between DFT and MD efforts can resolve the atomic
length scale of chemical representation and its influence on simple material
properties.

## METHODS

### Definition of local and global characteristic metrics

For a simulation cell of arbitrary size, a single point in global composition
space $c_{\rm global}$ is equivalent to the average chemistry of all its local
atomic environments, $\bar{c}_{\rm local}$:

$c_{\rm global}=\frac{1}{N}\sum^{N}_{i=1}c^{{\rm atom},i}_{\rm
local}=\bar{c}_{\rm local}(r_{\rm cut})$ (1)

where $N$ is the total number of atoms (i.e., local chemical environments),
$c^{{\rm atom},i}_{\rm local}(r_{\rm cut})$ is the local composition of atom
$i$ given a sampling radius of $r_{\rm cut}$. This metric considers the joint
probability of chemical occupations of all neighbor shells within $r_{\rm
cut}$, though for smaller $r_{\rm cut}$ can be approximated by other well-
known chemical ordering approaches (e.g., an average of Warren-Cowley
parameters 29).

For the analysis of local composition, we define two metrics in composition
space, the composition deviation $\lambda$, and the convex hull volume
fraction $V_{\rm f}$. The former value is the magnitude of a vector mapping of
any two points in composition space:

$\lambda^{2}=\sum_{i=\\\ {\rm Mo,Nb,Ta,Ti}}(c_{\rm global}^{i}-c_{\rm
local}^{i})^{2}$ (2)

The latter value requires an extra step. We can create a surface using the
outermost points in the point clouds shown in 1b, or in other words, a convex
hull. The volume enclosed by the surface of that convex hull $V_{hull}$
encompasses a certain space in a unit tetrahedron, whose side length in 4
dimensions is $\sqrt{2}$. By dividing the hull volume by the volume of a unit
tetrahedron, we gain a fraction representing the span of a given sample in
composition space:

$V_{\rm f}=V_{\rm hull}\cdot\frac{1}{3}$ (3)

### Density functional theory calculations for training set

DFT training data for the MoNbTaTi ML-IAP developed in this work is based on
the work of Startt et al. 6. All data for that study was generated using the
Vienna Ab-Initio Simulation Package (VASP) 31, 32, 33 using plane wave basis
sets for the orbital wavefunctions and the the projector augmented wave (PAW)
method 34, 35 to describe interactions between the core and valence-state
electrons. The CCA training data in this work was constructed around two
primary composition classes: compositionally varied quaternary compounds and
pure elemental systems. Within each class, DFT calculations were performed to
extract forces from several material state conditions to be used in the
training of the ML-IAP. These states included: structures with displaced
atoms, isotropically and uniaxially strained structures, surfaces, and high-
temperature NVT ensemble states.

Regarding the simulations for the quaternary compounds, we calculated the
aforementioned properties for specific compositional ranges following the
composition nomenclature: Ax(BCD)1-x, where all four metals (Mo, Nb, Ta, and
Ti) were taken in turn as the ’A’ component. The value of x was sampled at x =
0.125, 0.1667, 0.2083, 0.25, 0.2917, 0.3333, 0.375, and 0.5 atomic %. For
quaternary compounds, we performed structural relaxations of four unique
special quasirandom structure (SQS) 36, 72-atom supercells (a $4\times 3\times
3$ multiplication of the 2-atom conventional BCC unit cell) at each
composition point. These structures were modelled using a $3\times 5\times 5$
Monkhorst-Pack k-point grid. For all four structures belonging to each
composition, an atomic displacement calculation (i.e. an elastic constant
calculation) was performed, where each atom was displaced a total of six times
(two times $\pm$ along each Cartesian coordinate). For each displacement all
atoms are fixed in place and the electronic wave-function and atomic forces
are minimized. For other material conditions, a subset of calculations were
performed over a reduced range of composition for x = 0.1667, 0.25, 0.3333
atomic %. Uniaxially and isotropically strained supercells were modelled for
this reduced composition range and over compression to elongation ranges of
0.86 to 1.06 and 0.94 to 1.03 times the equilibrium lattice spacing,
respectfully. Surface supercells were modeled for the BCC $100$ (72 atoms),
$110$ (72), $111$ (162); and $112$ (128) surfaces. Lastly, thermodynamic ab-
initio MD (AIMD) NVT simulations were performed at temperatures of 300K,
1200K, 2400K, and 3200K, and at volumes set to match the expected thermal
expansion as predicted by 6 at each temperature.

Regarding simulations for the pure elements, we simulated each elemental
component of the MoNbTaTi quaternary individually in their pure ground state
structures (i.e., Mo, Nb, Ta - BCC, Ti - HCP (hexagonally close-packed)).
These systems were modelled over a similar set of deformation, strain, and
thermodynamic conditions as the quaternary compounds, with some notable
differences. Given the simplicity of the atomic structures, structural
relaxations were performed on smaller simulation cells, utilizing far denser
k-point grids. Additionally, surface slabs were constructed and minimized for
the 8-9 most stable surfaces known for each element as listed on the materials
project website 37. AIMD NVT simulations at were carried out at 300K, 1200K
and one point at least 100K above the known melting temperature of each metal.
Volumes at each temperature were again set to match the expected lattice
constant according to known thermal expansions.

For both classes of simulations, certain simulation parameters were kept
constant. Notably, plane-wave energy cutoffs were kept at 400 eV for all
material systems. Wavefunction energy and ionic force minimization criteria
were set to 10-6 eV and 0.02 eV/Å2, respectively. Partial orbital occupancies
were set according to a Gaussian smearing scheme using a smearing width of
0.02 eV. Exchange and correlation effects were handled according the commonly
employed Perdew, Burke, and Ernzerhof formalism of the generalized gradient
approximation 38. Lastly, the effects of spin-polarization were found to be
negligible for all systems after rigorous testing and screening and so were
excluded from calculations in the dataset. For more detailed information and
an in-depth analysis of some features of the DFT training set, see Ref. 6.

### Interatomic Potential Construction

The energy model is constructed using the standard form of SNAP that is
described in more detail in earlier publications.39, 40 The total potential
energy of a configuration of atoms is written as the sum of atomic energies
combined with an additional reference potential,

$E({\bf r}^{N})=E_{ref}({\bf r}^{N})+\sum_{i=1}^{N}E_{i}({\bf r}^{N}),$ (4)

where $E$ is the total potential, ${\bf r}^{N}$ are the positions of the $N$
atoms in the configuration, $E_{ref}$ is the reference energy, and $E_{i}$ is
the atomic energy of atom $i$. The atomic energy of atom $i$ is expressed as a
sum of the bispectrum components ${\bf B}_{i}$ for that atom weighted by
regression coefficients

$\displaystyle E_{i}({\bf r}^{N})$ $\displaystyle=$
$\displaystyle{{\bm{\beta}}}_{\nu_{i}}\cdot({{\bf B}}_{i}-{{\bf
B}}_{0\nu_{i}})\,,$ (5)

where the elements of the vector ${{\bm{\beta}}}_{\nu}$ are constant linear
coefficients for atoms of element ${\nu}$ whose values are determined in
training. The vector ${{\bf B}}_{i}$ is a flattened list of bispectrum
components for atom $i$, while ${{\bf B}}_{0\nu}$ is the list of bispectrum
components for an isolated atom of type $\nu$. By construction, the energy of
an isolated atom is zero. The bispectrum components are real, rotationally
invariant triple-products of four-dimensional hyperspherical harmonics ${\bf
U}_{j}$ 41

$\displaystyle B_{j_{1}j_{2}j}$ $\displaystyle=$ $\displaystyle{\bf
U}_{j_{1}}\otimes_{j_{1}j_{2}}^{j}{\bf U}_{j_{2}}\colon{\bf U}_{j}^{*}\,,$ (6)

where symbol $\otimes_{j_{1}j_{2}}^{j}$ indicates a Clebsch-Gordan product of
two matrices of arbitrary rank, while $:$ corresponds to an element-wise
scalar product of two matrices of equal rank. The total hyperspherical
harmonics for a central atom $i$ are written as sums over neighbor
contributions,

$\displaystyle{\bf U}_{j}$ $\displaystyle=$ $\displaystyle{\bf u}_{j}({\bf
0})+\sum_{k\in\mathcal{N}(i)}~{}f_{c}(r_{ik})w_{\nu_{k}}{\bf u}_{j}({\bf
r}_{ik})\,,$ (7)

where the summation is over all neighbor atoms $k$ within a cutoff distance
$R_{\nu_{i}\nu_{k}}$ of atom $i$. Atoms of different chemical elements are
distinguished by the element weights $w_{\nu}$. The radial cutoff function
$f_{c}(r)$ ensures that atomic contributions go smoothly to zero as $r$
approaches $R$ from below. We used the 55 lowest order bispectrum components
$B_{j_{1}j_{2}j}$ with half-integral indices restricted to the range $0\leq
j_{2}\leq j_{1}\leq j\leq 4$. The SNAP element weights and cutoff distances
were optimized for each element using the genetic algorithm search described
below. To ensure strong repulsion at short separations between all pairs of
atoms, a short-ranged ZBL42 reference potential was added ($Z$ = 44.5,
$R_{cut}$ = 5.0 Å).

### Multi-compositional refractory CCA ML-IAP for MoNbTaTi

In this section, we will briefly describe the fitting process and the
performance of the MoNbTaTi refractory CCA ML-IAP. The primary target in
fitting this potential was to ensure accuracy across a wide range of global
composition space, with a particular focus on replicating DFT predicted
elastic properties. Because elasticity is a property measured over an entire
simulation cell, and not locally, we did not apply the local composition
analysis techniques outlined above to optimize the fitting process. Instead,
we followed the global composition sampling scheme detailed in the ”Density
functional theory calculations for training set” described above.

For the present the case of the MoNbTaTi CCA ML-IAP, two methodologies were
combined: linear regression to generate individual IAPs and a genetic
algorithm to explore parameter space. The first method is the core model
connecting the DFT training set and a resulting IAP. We use a linear
regression scheme whose general form is:

$\hat{\bm{\sigma}}=\underset{\bm{\beta}}{\operatorname{argmin}}(\|\bm{\epsilon}\circ(D\bm{\beta}-T)\|^{2}-\gamma_{n}~{}\|\bm{\beta}\|^{n})$
(8)

$\hat{\bm{\sigma}}$ represents regression error between the bispectrum
descriptor predictions, $D$ and the training data reference energies, $T$.
$\gamma_{n}~{}\|\bm{\beta}\|^{n}$ is a regularization penalty of order $n$ and
weight $\gamma_{n}$.

To implement this scheme, descriptors of the training structures from the DFT
must be reproduced in MD simulations. For this purpose, we used the open-
source software FitSNAP 18 (GitHub link: https://github.com/FitSNAP/FitSNAP),
which parses configurations from DFT into LAMMPS. Once each structures’
corresponding bispectrum descriptors, reference energies, and reference forces
have been calculated, FitSNAP retrieves that information from LAMMPS to form
the D matrix from Eq. 8. The generation process is completed by solving for
$\hat{\bm{\beta}}$ using singular value decomposition and outputting an IAP.
Though the first process as described above will generate an IAP for use in
LAMMPS, it does not guarantee that it will be optimal for applications of
interest. To address the issue of tuning variables, we wrap a single-objective
genetic algorithm (GA) around FitSNAP, as implemented by the DAKOTA
optimization software 19. The purpose of the GA is to take a user-generated
selection of target IAP properties and discover the optimal set of variables
to use within the core fitting calcuations in FitSNAP. To accomplish this, the
user creates a series of short simulations in LAMMPS that are run and are used
to evaluate a generated IAP’s overall quality.

Including multi-compositional data in the training set is not enough to
guarantee that an IAP will perform well across those compositions. We found
the key technique to be the setup of the simulations that feed into the GA’s
objective functions. For the purposes of this work, we aimed to optimize ML-
IAPs on the first-principles elasticity data (see “Density functional theory
calculations for training set” subsection above), which indicate that MoNbTaTi
alloys undergo linear changes in elastic moduli with the enrichment and
depletion of single elements. Thus, our post-fit testing and matching GA
objective functions were designed to make the GA especially sensitive to the
rate of change of the bulk ($B$) and shear ($G$) moduli as compositions are
varied. Fitting the moduli to slopes instead of single-composition values
forces the GA to favor not only low regression errors on single-composition
fits, but also couple those errors’ tendencies to favor accurate calculations
of $B$ and $G$. To fit these slopes, it is also necessary to create objective
functions that minimize errors on the components of the elasticity tensor,
$\mathbb{C}_{xx}$. As the stable MoNbTaTi single-phase refractory CCAs take on
BCC structure, only three of the tensor components ($\mathbb{C}_{11}$,
$\mathbb{C}_{12}$, and $\mathbb{C}_{44}$), which we will also refer to as the
elastic constants, need to be calculated per IAP and composition tested.
Having obtained those, the calculation of the bulk and shear moduli is as
follows:

$B=\frac{\mathbb{C}_{11}+2\mathbb{C}_{12}}{3}$ (9)
$G=\frac{1}{2}\left[\frac{\mathbb{C}_{11}-\mathbb{C}_{12}+3\mathbb{C}_{44}}{5}+\frac{5\mathbb{C}_{44}\left(\mathbb{C}_{11}-\mathbb{C}_{12}\right)}{4\mathbb{C}_{44}+3\left(\mathbb{C}_{11}-\mathbb{C}_{12}\right)}\right]$
(10)

Once calculated for three separate compositions, a slope for $B$ and $G$ for
the ML-IAP can be fit and tested against slopes calculated from the training
data.

### Vacancy energy simulations

In alloys, the vacancy formation energy of a given atomic species $\mu$,
$E_{\rm vac}^{\mu}$, is found by first calculating the average cohesive energy
of all atoms in a cell $E_{\rm cohesive}^{\rm all}$, removing an atom of one
species $\mu$, and then recalculating the cell’s new energy $E_{\rm
removed}^{\mu}$. The difference of the new energy and the total potential
energy corresponding to a cell one atom smaller gives the penalty for vacancy
formation for that species $\mu$:

$E_{\rm vac,f}^{\mu}=E_{\rm removed}^{\mu}-E_{\rm total}^{\rm
all}\cdot\frac{N_{\rm cell}-1}{N_{\rm cell}}$ (11)

### Hybrid Monte Carlo/molecular dynamics (MC/MD) simulations

To induce chemical reordering in the MoNbTaTi SNAP potential, we used a hybrid
Monte Carlo/molecular dynamics procedure on one equiatomic random solid
solution cube with side length of $s\cdot a_{0}=40\cdot 3.2546$ Å , totalling
128,000 atoms. The random solid solution was initially relaxed at T = 300K for
10 ps in the NPT ensemble (Nose-Hoover barostat), using a 1 fs timestep. For
the hybrid algorthim, the same settings were used for all further MD steps.
The MC element type swaps were conducted at intervals of 50 MD steps (0.050
ps) using a series of _fix atom/swap_ in LAMMPS at a MC temperature of 300K.
In total, 150 ps of hybrid MC/MD (75,000 MC steps with 768,000 atom swap
attempts, 138,748 accepted) were conducted to achieve the results found in
Fig. 5a. Warren-Cowley parameters 29 were used to calculate trends in pair
ordering for the 1st nearest-neighbor shell of central atom $i$ relative to
elements A and B in Fig. 5b:

$\alpha_{i}^{AB}=1-\frac{P_{i}^{AB}}{c_{j}}$ (12)

Details of the convergence of potential energy through MC steps can be found
in the Supplemental Information.

## DATA AVAILABILITY

The full training data sets as well as all validation and test cases are
available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

This work was supported by the U.S. Department of Energy, Office of Fusion
Energy Sciences (OFES) under Field Work Proposal Number 20-023149, and the the
Center for Integrated Nanotechnologies, an Office of Science user facility
operated for the U.S. Department of Energy. This article has been authored by
an employee of National Technology & Engineering Solutions of Sandia, LLC
under Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE). The
employee owns all right, title and interest in and to the article and is
solely responsible for its contents. The United States Government retains and
the publisher, by accepting the article for publication, acknowledges that the
United States Government retains a non-exclusive, paid-up, irrevocable, world-
wide license to publish or reproduce the published form of this article or
allow others to do so, for United States Government purposes. The DOE will
provide public access to these results of federally sponsored research in
accordance with the DOE Public Access Plan
https://www.energy.gov/downloads/doe-public-access-plan.

## AUTHOR CONTRIBUTIONS

MJM : Development of interatomic potential, performed MD simulations,
implementation of local composition analysis code, chemical complexity
analysis. JS : Performed DFT simulations, chemical complexity analysis. APT :
Implemented interatomic potential in LAMMPS, chemical complexity analysis. RD
: Validation of DFT, MD and interatomic potential. MAW : Development of
interatomic potential, chemical complexity analysis. All authors participated
in conceiving the research and writing the manuscript.

## COMPETING INTERESTS

The authors declare no competing interests.

## ADDITIONAL INFORMATION

Supplementary Information is available for this paper.

## REFERENCES

  * [1] Senkov, O., Miller, J., Miracle, D. & Woodward, C.  Accelerated exploration of multi-principal element alloys for structural applications.  _Calphad_ 50, 32–48 (2015).  URL https://doi.org/10.1016/j.calphad.2015.04.009.
  * [2] Yao, H. _et al._ MoNbTaV medium-entropy alloy.  _Entropy_ 18, 189 (2016).  URL https://doi.org/10.3390/e18050189.
  * [3] Senkov, O., Miracle, D., Chaput, K. & Couzinié, J.-P.  Development and exploration of refractory high entropy alloys – A review.  _J. Mater. Res._ 33, 3092–3128 (2018).  URL https://doi.org/10.1557/jmr.2018.153.
  * [4] Coury, F., Kaufman, M. & Clarke, A.  Solid-solution strengthening in refractory high entropy alloys.  _Acta Mater._ 175, 66–81 (2019).  URL https://doi.org/10.1016/j.actamat.2019.06.006.
  * [5] Senkov, O., Gorsse, S. & Miracle, D.  High temperature strength of refractory complex concentrated alloys.  _Acta Mater._ 175, 394–405 (2019).  URL https://doi.org/10.1016/j.actamat.2019.06.032.
  * [6] Startt, J., Kustas, A., Pegues, J., Yang, P. & Dingreville, R.  Compositional effects on the mechanical and thermal properties of MoNbTaTi refractory complex concentrated alloys.  _Mater. Des._ 213, 110311 (2022).  URL https://doi.org/10.1016/j.matdes.2021.110311.
  * [7] Cui, Y., Shen, J., Manladan, S. M., Geng, K. & Hu, S.  Wear resistance of FeCoCrNiMnAl${}_{\textrm{x}}$ high-entropy alloy coatings at high temperature.  _Appl. Surf. Sci._ 512, 145736 (2020).  URL https://doi.org/10.1016/j.apsusc.2020.145736.
  * [8] Du, Y. _et al._ Mechanical and tribological performance of CoCrNiHf${}_{\textrm{x}}$ eutectic medium-entropy alloys.  _J. Mater. Sci. Technol._ 90, 194–204 (2021).  URL https://doi.org/10.1016/j.jmst.2021.03.023.
  * [9] Liu, Y. _et al._ Effect of Al content on high temperature oxidation resistance of Al${}_{\textrm{x}}$CoCrCuFeNi high entropy alloys (x= 0, 0.5, 1, 1.5, 2).  _Vacuum_ 169, 108837 (2019).  URL https://doi.org/10.1016/j.vacuum.2019.108837.
  * [10] Esmaily, M. _et al._ High-temperature oxidation behaviour of Al${}_{\textrm{x}}$FeCrCoNi and AlTiVCr compositionally complex alloys.  _npj Mater. Degrad._ 4, 1–10 (2020).  URL https://doi.org/10.1038/s41529-020-00129-2.
  * [11] Scully, J. R. _et al._ Controlling the corrosion resistance of multi-principal element alloys.  _Scripta Mater._ 188, 96–101 (2020).  URL https://doi.org/10.1016/j.scriptamat.2020.06.065.
  * [12] Zhao, R.-F., Ren, B., Zhang, G.-P., Liu, Z.-X. & Zhang, J.-J.  Effect of Co content on the phase transition and magnetic properties of Co${}_{\textrm{x}}$CrCuFeMnNi high-entropy alloy powders.  _J. Magn. Magn. Mater._ 468, 14–24 (2018).  URL https://doi.org/10.1016/j.jmmm.2018.07.072.
  * [13] Chaudhary, V., Chaudhary, R., Banerjee, R. & Ramanujan, R.  Accelerated and conventional development of magnetic high entropy alloys.  _Mater. Today_ 49, 231–252 (2021).  URL https://doi.org/10.1016/j.mattod.2021.03.018.
  * [14] Mishra, R. K., Kumari, P., Gupta, A. K. & Shahi, R. R.  Design and development of Co35Cr5Fe${}_{20-\textrm{x}}$Ni${}_{20+\textrm{x}}$Ti${}_{2}0$ high entropy alloy with excellent magnetic softness.  _J. Alloys Compd_ 889, 161773 (2021).  URL https://doi.org/10.1016/j.jallcom.2021.161773.
  * [15] Byggmästar, J., Nordlund, K. & Djurabekova, F.  Modeling refractory high-entropy alloys with efficient machine-learned interatomic potentials: Defects and segregation.  _Phys. Rev. B_ 104, 104101 (2021).  URL https://doi.org/10.1103/PhysRevB.104.104101.
  * [16] Li, X.-G., Chen, C., Zheng, H., Zuo, Y. & Ong, S. P.  Complex strengthening mechanisms in the NbMoTaW multi-principal element alloy.  _npj Comput. Mater._ 6, 1–10 (2020).  URL https://doi.org/10.1038/s41524-020-0339-0.
  * [17] Yin, S. _et al._ Atomistic simulations of dislocation mobility in refractory high-entropy alloys and the effect of chemical short-range order.  _Nat. Commun._ 12, 1–14 (2021).  URL https://doi.org/10.1038/s41467-021-25134-0.
  * [18] Rohskopf, A. _et al._ Fitsnap: Atomistic machine learning with lammps.  _Journal of Open Source Software_ 8, 5118 (2023).  URL https://doi.org/10.21105/joss.05118.
  * [19] Adams, B. _et al._ Dakota: a multilevel parallel object-oriented framework for design optimization, parameter estimation, uncertainty quantification, and sensitivity analysis.  Tech. Rep. SAND2010-2183, Sandia National Laboratories (2009).  URL https://doi.org/10.2172/1630693.
  * [20] Zhou, X., He, S. & Marian, J.  Vacancy energetics and diffusivities in the equiatomic multielement Nb-Mo-Ta-W alloy.  _Materials_ 15, 5468 (2022).  URL https://doi.org/10.3390/ma15155468.
  * [21] Cusentino, M., Wood, M. & Dingreville, R.  Compositional and structural origins of radiation damage mitigation in high-entropy alloys.  _J. Appl. Phys._ 128, 125904 (2020).  URL https://doi.org/10.1063/5.0024014.
  * [22] Li, C. _et al._ First principle study of magnetism and vacancy energetics in a near equimolar NiFeMnCr high entropy alloy.  _J. Appl. Phys._ 125, 155103 (2019).  URL https://doi.org/10.1063/1.5086172.
  * [23] Stukowski, A. & Arsenlis, A.  On the elastic–plastic decomposition of crystal deformation at the atomic scale.  _Model. Simul. Mat. Sci. Eng._ 20, 035012 (2012).  URL https://doi.org/10.1088/0965-0393/20/3/035012.
  * [24] Ferrari, A. _et al._ Frontiers in atomistic simulations of high entropy alloys.  _J. Appl. Phys._ 128, 150901 (2020).  URL https://doi.org/10.1063/5.0025310.
  * [25] Antillon, E., Woodward, C., Rao, S. I. & Akdim, B.  Chemical short range order strengthening in BCC complex concentrated alloys.  _Acta Mater._ 215, 117012 (2021).  URL https://doi.org/10.1016/j.actamat.2021.117012.
  * [26] Miracle, D. B. & Senkov, O. N.  A critical review of high entropy alloys and related concepts.  _Acta Mater._ 122, 448–511 (2017).  URL https://doi.org/10.1016/j.actamat.2016.08.081.
  * [27] Widom, M., Huhn, W. P., Maiti, S. & Steurer, W.  Hybrid Monte Carlo/molecular dynamics simulation of a refractory metal high entropy alloy.  _Metall. Mater. Trans. A_ 45, 196–200 (2014).  URL https://doi.org/10.1038/s41524-020-0339-0.
  * [28] Körmann, F., Ruban, A. V. & Sluiter, M. H.  Long-ranged interactions in bcc NbMoTaW high-entropy alloys.  _Mater. Res. Lett._ 5, 35–40 (2017).  URL https://doi.org/10.1080/21663831.2016.1198837.
  * [29] Cowley, J. M.  An approximate theory of order in alloys.  _Phys. Rev._ 77, 669 (1950).  URL https://doi.org/10.1103/PhysRev.77.669.
  * [30] Zhang, R. _et al._ Short-range order and its impact on the CrCoNi medium-entropy alloy.  _Nature_ 581, 283–287 (2020).  URL https://doi.org/10.1038/s41586-020-2275-z.
  * [31] Kresse, G. & Hafner, J.  Ab initio molecular dynamics for liquid metals.  _Phys. Rev. B_ 47, 558–561 (1993).  URL https://doi.org/10.1103/PhysRevB.47.558.
  * [32] Kresse, G. & Hafner, J.  Ab initio molecular-dynamics simulation of the liquid-metalamorphous- semiconductor transition in germanium.  _Phys. Rev. B_ 49, 14251–14269 (1994).  URL https://doi.org/10.1103/PhysRevB.49.14251.
  * [33] Kresse, G.  Ab initio molecular dynamics for liquid metals.  _J. Non-Cryst. Solids_ 192-193, 222–229 (1995).  URL https://doi.org/10.1016/0022-3093(95)00355-X.
  * [34] Blochl, P. E.  Projector augmented-wave method.  _Phys. Rev. B_ 50, 17979 (1994).  URL https://doi.org/10.1103/PhysRevB.50.17953.
  * [35] Kresse, G. & Joubert, D.  From ultrasoft pseudopotentials to the projector augmented-wave method.  _Phys. Rev. B_ 59, 1758–1775 (1999).  URL https://doi.org/10.1103/PhysRevB.59.1758.
  * [36] Van de Walle, A. _et al._ Efficient stochastic generation of special quasirandom structures.  _Calphad_ 42, 13–18 (2013).  URL https://doi.org/10.1016/j.calphad.2013.06.006.
  * [37] Jain, A. _et al._ Commentary: The materials project: A materials genome approach to accelerating materials innovation.  _APL Mater._ 1, 011002 (2013).  URL https://doi.org/10.1063/1.4812323.
  * [38] Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized gradient approximation made simple.  _Phys. Rev. Lett._ 77, 3865–3868 (1996).  URL https://doi.org/10.1103/PhysRevLett.77.3865.
  * [39] Thompson, A. P., Swiler, L. P., Trott, C. R., Foiles, S. M. & Tucker, G. J.  Spectral neighbor analysis method for automated generation of quantum-accurate interatomic potentials.  _J. Comput. Phys._ 285, 316–330 (2015).  URL http://dx.doi.org/10.1016/j.jcp.2014.12.018.
  * [40] Wood, M. A. & Thompson, A. P.  Extending the accuracy of the SNAP interatomic potential form.  _J. Chem. Phys._ 148, 1–10 (2018).  URL http://dx.doi.org/10.1063/1.5017641.  eprint 1711.11131.
  * [41] Bartók, A. P., Payne, M. C., Kondor, R. & Csányi, G.  Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons.  _Phys. Rev. Lett._ 104, 136403 (2010).  URL https://doi.org/10.1103/PhysRevLett.104.136403.
  * [42] Ziegler, J., Biersack, J. & Littmark, U.  _The Stopping and Range of Ions in Solids_ (Pergamon, 1985).
  * [43] Stukowski, A.  Visualization and analysis of atomistic simulation data with OVITO-the Open Visualization Tool.  _Model. Simul. Mat. Sci. Eng._ 18 (2010).  URL https://doi.org/10.1088/0965-0393/18/1/015012.

SUPPLEMENTAL FIG. 1:  a Average potential energy per atom vs. Monte Carlo (MC)
steps, showing the reduction in overall system energy due to successful atom
swaps at $T=300K$. The simulations were stopped at 150 ps (150,000 MD steps,
75,000 MC steps, 768,000 total swap attempts with 138,748 successful swaps).
The inset shows 2,235 Mo (green) and Ta (yellow) atoms identified as
B2-ordered by the Ovito software’s ’Polyhedral Template Matching’ chemical
ordering modifier 43.