# The role of disorder in the synthesis of metastable ternary nitrides

Rachel Woods-Robinson1,2,3, Vladan Stevanović4,3, Stephan Lany3, Karen
Heinselman3, Kristin A. Persson5,2, Andriy Zakutayev3 1Applied Science and
Technology Graduate Group, University of California at Berkeley, Berkeley, CA,
94720 USA, 2Materials Sciences Division, Lawrence Berkeley National
Laboratory, Berkeley, CA, 94720 USA, 3Materials and Chemistry Science and
Technology Directorate, National Renewable Energy Laboratory, Golden,
Colorado, 80401 USA 4Department of Physics, Colorado School of Mines, Golden,
Colorado, 80401 USA, 5Department of Materials Science and Engineering,
University of California at Berkeley, Berkeley, CA, 94720 USA

(August 27, 2024)

###### Abstract

In materials science, it is often assumed that the most thermodynamically
stable crystal structure is the easiest polymorph to synthesize. Ternary
nitride materials, with many possible metastable polymorphs, provide a domain
to test this assumption; for example, ZnZrN2 is predicted to have an unusual
layered “wurtsalt” crystal structure and exciting optoelectronic properties,
but can it be realized experimentally? Here, we synthesize hundreds of
Zn_$x$Zr_1-$x$N_$y$ thin film samples, and find metastable rocksalt-derived or
boron-nitride-derived structures rather than the thermodynamically stable
wurtsalt structure. Computations show that this divergence is due to both
entropic and enthalpic stabilization effects. Rocksalt- and boron-nitride-
derived structures become the most stable polymorphs in the presence of
disorder because of higher tolerances to cation cross-substitution and off-
stoichiometry than the wurtsalt structure. By understanding the role of
disorder tolerance in the synthesis of metastable polymorphs, we can enable
more accurate predictions of synthesizable crystal structures and their
achievable material properties.

## Introduction

Computational materials discovery is a rapidly progressing research field,
with the potential to revolutionize how we design and scale up materials.
However, determining whether a given predicted crystalline material is
actually experimentally synthesizable remains a key challenge. One common
assumption in research guided by density functional theory (DFT) is that the
“true” DFT-predicted ground state structure, or structures with energies near
the ground state energy, are the most likely to be experimentally realized.
Conversely, another assumption is that increased energetic instability (i.e.
formation energy farther away from the ground state energy) correlates with an
increased difficulty to synthesize. However, neither of these assumptions
necessarily hold, as demonstrated by multiple experimental and computational
studies.[1] Recent work has emerged to further explore synthesizability in
metastable materials,[2, 3, 4, 5] but so far the research community still
cannot definitively answer the simple question: “can I synthesize this
material?”[6] Thus, as materials databases grow and structure predictions
yield new predicted compounds for high-throughput screenings, it is
increasingly pertinent that the materials discovery community develops
comprehensive methods for assessing synthesizability so that we can avoid
misleading false positives and negatives.

Nitride materials provide a compelling framework through which to examine
synthesizability, in because they are more likely than any other anion class
to crystallize in metastable phases.[7, 3, 4, 8] Recent computational
predictions have yielded a multitude of new ternary nitride materials to
explore,[9, 10] yet an understanding of which polymorphs are indeed
synthesizable remains elusive. The chemical and structural richness of this
emerging class of materials, including their mixed ionic-covalent nature
compared to oxides, provides new candidates for various applications such as
hydrogen storage, photovoltaic (PV) devices, and light-emitting diodes (LEDs).
One such emerging class of ternary nitrides is the II-IV-N2 family, ternary
analogs of GaN and promising candidate for PV absorbers and green LEDs. II-
IV-N2 materials are commonly studied in two prototype classes: (1) wurtzite-
derived (WZ) structures, with four-fold coordinated cations (e.g. Zn-based
ZnSnN2, ZnGeN2, ZnSnP2, ZnSiP2),[11] and (2) rocksalt-derived structures (RS),
with six-fold coordinated cations (e.g. MgTMN2).[12] Some compounds (e.g.
MgSnN2) have been shown to co-crystallize in both of these configurations at
certain growth conditions, such as at ambient temperature at Mg-rich
stoichiometries or at increased synthesis temperature on GaN substrates.[7,
13, 14] However, these two structure classes are just a small subset of
possible structure classes in the rich space of ternary nitrides; it remains
unknown which other II-IV-N2 polymorph structures and chemistries are
stabilizable.

Figure 1: (a) Representative crystal structures for each of the six ordered
structure prototype classes, and (b) energy ordering of the predicted ordered
polymorphs, calculated by SCAN and grouped by structure prototype class, with
labeled experimentally synthesized phases from this study. Horizontal solid
lines in (b) correspond to calculated $E_{\mathrm{hull}}$ values of individual
polymorphs, and shaded regions correspond to the range of $E_{\mathrm{hull}}$
for a given prototype class. Heavily distorted versions of the prototype
structures, designated with a “d-” prefix, are plotted separately on the right
of (b).

Of particular interest to this study is the experimentally empty region of
phase space in the zinc zirconium nitride (Zn-Zr-N) material system, in
particular at its II-IV-N2 composition ZnZrN2, which serves as a case study to
gain insight for ternary nitrides as a whole. In contrast to other II-IV-N2
materials, ZnZrN2 (as well as isoelectronic ZnHfN2) has a predicted ground
state $P3m1$ (156) space group structure—a layer of Zn atoms tetrahedrally
coordinated by N (wurtzite-like), a layer of Zr atoms octahedrally coordinated
by N (rocksalt-like), and alternating Zn and Zr layers (abab)—which has been
corroborated by three different computational studies using three distinct
structure prediction algorithms.[15, 16, 17] This structure is analogous to
sulfosalt ScCuS2, though we could not locate a name for this structure;[18,
19] thus, we herein refer to this structure type as “wurtsalt” (WS), an
amalgam of wurt-zite and rock-salt, and depict this structure in the top left
of Figure 1(a), alongside other polymorphs. Despite these predictions, no
semiconducting nitride materials in the Zn_$x$Zr_1-$x$N_$y$ ternary space have
ever been confirmed experimentally, and it has not yet been confirmed whether
any other polymorphs exist.

In this study, we demonstrate that certain polymorphs can be preferentially
stabilized or destabilized due to their tolerance to disorder. To do so, we
use combinatorial sputter synthesis to explore the full cation phase space of
Zn_$x$Zr_1-$x$N_$y$ alloys, focusing on the region where $y\approx$ 0 and
$x\approx$ 0, and computationally investigate a set of 28 possible ZnZrN2
polymorphs. Under our experimental conditions, the predicted WS ground state
is not synthesizable; rather, a thermodynamically metastable rocksalt (RS)
polymorph is energetically stabilized, and a metastable hexagonal boron-
nitride-derived (BN) phase is observed at higher Zn concentrations
(x$\gtrsim$0.5). We then use a series of computational methods to explain
these findings, taking into account configurational entropy and enthalpy, as
well as tolerance to off-stoichiometry. These results suggest that it is
necessary to consider the effects of disorder tolerance on energetic
stabilization in possible polymorphs when investigating new ternary nitrides
and new computationally predicted materials in general.

Table 1: Representative ordered polymorphs from each prototype class with the
lowest $E_{\mathrm{hull}}$, calculated with SCAN (see SI for full list of
polymorphs and energies)

[HTML]EFEFEF  | Prototype
---
class
| Space
---
group
| # of
---
atoms†
| $\bm{E_{\mathrm{hull}}}$
---
(eV/at)
| $\bm{E_{\mathrm{G}}}$
---
(eV)
| $\bm{E_{\mathrm{G}}^{\mathrm{d}}}$
---
(eV)
| $\bm{m^{*}_{\mathrm{e}}}$
---
| $\bm{m^{*}_{\mathrm{h}}}$
---
WS | $P3m1$ | 4 | 0.0 | 2.47 | 3.10 | 7.30 | 1.69
LC | $Pca2_{1}$ | 16 | 0.0106 | 1.63 | 1.63 | 1.33 | 1.87
d-HX | $P2_{1}/c$ | 16 | 0.0294 | 2.62 | 2.71 | 3.88 | 2.16
d-WS | $Cm$ | 32 | 0.0312 | 2.18 | 2.18 | 1.56 | 1.25
BN | $Cm$ | 16 | 0.0327 | 2.01 | 2.01 | 1.36 | 1.49
d-RS | $Pc$ | 16 | 0.0373 | 2.22 | 2.47 | 3.41 | 2.17
WZ | $Pmc2_{1}$ | 8 | 0.0496 | 2.53 | 3.23 | 0.62 | 3.62
RS | $I4_{1}/amd$ | 16 | 0.0807 | 1.15 | 1.87 | 0.83 | 1.96
ZB | $P\bar{4}m2$ | 4 | 0.0857 | 2.04 | 3.03 | 0.52 | 1.59

† Number of atoms in primitive unit cell

## Results and Discussion

### Identifying possible ZnZrN2 polymorphs

Although only the WS phase is reported in the Materials Project database
(ZrZnN2, mp-1014244),[20, 21] this unexplored Zn-Zr-N phase space could host a
variety of different structures. Many methods exist to determine possible
polymorphs and predict synthesizable compounds,[22] ranging from simple ionic
substitution, to kinetically limited minimization (KLM)[23], random structure
searches,[24] or more expensive evolutionary and genetic algorithms.[25] Since
no single method is fully representative of configurational space, we use the
combined methods of KLM,[23] the polymorph sampler method using random
structure searching,[2] and ionic substitution of prototypes from other
ternary nitrides[14] to predict 28 unique possible ordered ZnZrN2 polymorphs
and assess their relative phase stability. We classify these polymorphs into
six distinct structure prototype classes, with representative crystals for
each of these structure prototypes depicted in Figure 1(a), and adopt a naming
convention from binary analogs as follows: rocksalt-derived (“RS”) is an fcc
anion lattice with cations in Oh-coordinated voids, wurtzite-derived (“WZ”)
exhibits a structurally face-centered tetragonal anion lattice with
tetrahedral coordinated cations, zincblende-derived (“ZB”, i.e. chalcopyrite)
is an fcc anion lattice with cations in every other tetrahedral void, wurtsalt
(“WS”) presents alternating layers of octahedrally coordinated Zr and
tetrahedrally coordinated Zn (as discussed previously), and boron-nitride-
derived (“BN”) exhibits hexagonal sheets of various stackings (similar to
graphite, but a 3D structure with M-N bonds between c-axis layers). An
additional compound, with alternating 2D layers corresponding to the mineral
litharge (PbO, with a space group $P4/nmm$) and the mineral massicot (PbO,
with a space group $P2_{1}ca$), we name with the amalgam “lithicot” (“LC”). We
also separately categorize heavily distorted versions of three of these
prototypes with the prefix “d-”, using a tolerance developed from a structure
matching algorithm (see SI for classification scheme details). The structure
class “d-HX” (HX = hexagonal) represents structures that are distortions
between BN and WZ, which are commensurate to one another. We note that Zn3N2
crystallizes in an anti-bixbyite-derived phase (“BX”), which is observed
experimentally, but this structure is not included in our set of prototypes
since deriving an analogous topotactic ZnZrN2 BX structure requires removing
atoms and is not trivial.

These candidate polymorphs are relaxed using DFT, first with a GGA functional,
then with GGA with a Hubbard +U correction of 3 eV/atom for Zr,[26] and then
with the SCAN functional to verify polymorph ordering (see SI). SCAN is used
herein, unless specified otherwise, because it has been demonstrated to more
accurately predict polymorph orderings, with the trade-off of a higher
computational cost.[27, 28] Figure 1(b) plots the resulting energy ordering of
the 28 ordered structures, with each polymorph grouped and colored according
to its prototype structure, “$E_{\mathrm{hull}}$” indicating the computed
formation energy above the convex hull, and only including structures with
$E_{\mathrm{hull}}<$ 0.160 eV/atom. Individual polymorphs are designated with
horizontal lines, and the range of $E_{\mathrm{hull}}$ for a given prototype
class is shaded. With SCAN we confirm a predicted ground state WS ($P3m1$)
that lies on the convex hull, corroborating the literature.[15, 3] Other WS
polymorphs ($P6_{3}mc$, $P\bar{3}m1$) are low in energy, ranging from 0–0.025
eV/atom, and the LC structures ($Pca2_{1}$, $Iba2$) are the next-lowest in
energy. RS polymorphs are the highest in energy, ranging from 0.080–0.155
eV/atom. GGA and GGA+U yield similar energy orderings (c.f. SI), though LC is
the predicted ground state for GGA without a Hubbard U correction (see SI).
The $E_{\mathrm{hull}}$ values of the lowest energy structure in each
prototype class, as well as their calculated band gaps ($E_{\mathrm{G}}$) and
effective masses ($m^{*}_{\mathrm{e}}$ and $m^{*}_{\mathrm{h}}$) from SCAN,
are reported in Table 1, with the full list in the SI. Optoelectronic
properties vary significantly by structure. Most polymorphs have indirect gaps
except for the LC structures, most of the BN, some distorted structures, and
one RS. The WZ $Pna2_{1}$ polymorph exhibits the largest band gap
($E_{\mathrm{G}}$ $\approx$2.99 eV with SCAN), followed by d-HX, WS, ZB, d-RS
and d-WS with $E_{\mathrm{G}}>$2 eV, while RS has among the lowest band gaps
($\sim$0–1.67 eV, depending on cation ordering). Notably, the lowest-energy WS
$P3m1$ polymorph have exceptionally low hole effective mass
$m^{*}_{\mathrm{h}}$ values ($<$2) compared to $m^{*}_{\mathrm{e}}$ values
($<$7) while retaining a wide direct band gap, $E_{\mathrm{G}}^{\mathrm{d}}>3$
eV. This combination of electronic structure properties is unique among all
the considered polymorphs, and extremely rare for other chemistries outside of
the ZnZrN2 material system.

### Synthesis of metastable phases

Surprisingly, despite the existence of at least 19 predicted polymorphs with
lower formation energy (at $T$ $=0$), a RS phase with a high
$E_{\mathrm{hull}}$ is experimentally synthesized at low growth temperatures
and ZnZrN2 stoichiometry. Using combinatorial sputter synthesis, we grow a set
of approximately 400 samples in the Zn_$x$Zr_1-$x$N_$y$ ternary alloy system,
with cation concentration ranging across the full Zn–Zr composition space and
growth temperature ranging from ambient to 500°C. Films are grown on amorphous
fused silica substrates, with Zn and Zr metallic targets and N plasma, and the
combinatorial method is used to induce composition and temperature gradients
(see SI for details).[29] Rutherford backscattering spectrometry (RBS)
corroborates the cation concentration measured by X-ray fluorescence
spectroscopy for a set of representative samples grown at ambient temperature.
As shown in Figure 2(a), depicting anion composition, $y$ $=$ (N+O)/(Zn+Zr),
as a function of cation composition, $x$ $=$ Zn/(Zn+Zr), RBS indicates N-rich
compositions at Zn-poor samples, N-poor compositions in Zn-rich samples, and
approximately stoichiometric N at the ZnZrN2 composition of interest.
Additionally, a small but nonzero presence of oxygen likely substitutes for N
and resides on the film surface ($\sim$0.3 at. % in Zn-poor samples, $\sim$5
at. % in Zn-rich samples due to reaction of zinc nitride with ambient
atmosphere; see SI). Figure 2(a) suggests our samples have crystallized along
the Zn3N2–Zr3N4 tieline, as indicated by the exponential fit to RBS referenced
to the crossed markers, with an approximate stoichiometry of
Zn_$x$Zr_1-$x$N_$y$ where $y\approx(4-2x)/3$. We note this system could
alternately be expressed as “Zn_1+$x$Zr_1-$x$N_2+$y$” to emphasize off-
stoichiometry from ZnZrN2 (see SI discussion). For simplicity and generality
we will refer to experimental alloys as “Zn_$x$Zr_1-$x$N_$y$” herein since
multiple experimental phases are observed, and focus on varying $x$ since we
are not directly varying anion concentration.

Figure 2: (a) RBS measurements of anion composition as a function of cation
composition, and extrapolated with the black line. (b) Map of over 400 samples
in Zn-Zr-N experimental phase space, with cation composition $x=$ Zn/(Zn+Zr)
on the x-axis and calibrated deposition temperature on the y-axis. (c)
Representative diffraction patterns for 3 samples marked with rectangles in
(b), with “modeled” ZnZrN2 RS and BN from ensemble averages of XRD patterns
from the polymorph sampler method for ZnZrN2 compositions, and “standard”
Zn3N2 anti-bixbyite (BX) from the Materials Project database with smearing
bandwidth of 0.05.

Using high-throughput synchrotron X-ray diffraction (XRD) analysis, and
referencing the XRD patterns from the set of predicted polymorphs, we map the
crystalline phase space by linking the structure of each of the $\sim$400
samples to its corresponding ternary composition and calibrated deposition
temperature, $T_{\mathrm{dep}}$. The resulting map is plotted in Figure 2(b)
(see SI for details about patterns and $T_{\mathrm{dep}}$ calibration). Rather
than crystallizing in its predicted ground state WS structure at and around
the ZnZrN2 ($x$ = 0.5) composition (dotted vertical line), synchrotron
diffraction demonstrates a range of alloys predominantly in an RS crystal
structure, which is energetically unstable according to Figure 1. We also
observe the transformation to a second phase at higher $x$ compositions, a
hexagonal structure likely BN or d-HX, and a third anti-bixbyite-derived phase
(“BX”) at the highest $x$ concentrations near Zn3N2 (again, simulating BX
ZnZrN2 is beyond our scope). At higher temperatures, there are regions of
mixed phases of these three polymorphs, perhaps due to miscibility or Zn
volatility. We note that we attempted to synthesize high-temperature Zn-rich
samples, but no such samples were realized due to the high vapor pressure of
Zn under our growth conditions (see SI).

Figures 2(b) and (c) represent the key structural observations in
Zn_$x$Zr_1-$x$N_$y$, where $x$ is the cation ratio Zn/(Zn+Zr). First, a RS-
derived phase dominates from $x$ = 0 (Zr3N4) up to a threshold $x$ value,
which is approximately $x\approx 0.66$ at ambient temperature growth
conditions ($T_{\mathrm{dep}}\approx$ 65°C, bottom of figure) and which drops
as $T_{\mathrm{dep}}$ increases. A diffraction pattern for a representative RS
sample of $x\approx$ 0.45 is depicted in (c) in dark teal, compared to a
modeled RS diffraction pattern in light teal, simulated as ensemble-averages
from the polymorph sampler (see next section and SI for more details). As Zn
increases, the XRD peaks around $Q$ = 2.3 and 2.65 Å , corresponding to RS
(111) and (200), respectively (indices from the ZrN RS structure; see SI),
shift to higher $Q$ values, with the former strengthening and the latter
weakening (see SI). This trade-off is likely due to shifts in texturing, as
also commonly observed in other ternary nitrides.[30]

At the threshold $x$ value ($x$ $\approx$ 0.66 at ambient temperature), there
is a phase-change to a hexagonal BN-derived structure. Figure 2(c) depicts a
representative BN diffraction pattern for a sample with $x$ $\approx$ 0.56 in
dark red, with diffraction peaks at $Q$ values of $\sim$2.10, 2.45, and 2.55 Å
corresponding to BN (001), (100), and (101) reflections, respectively (with
the modeled BN pattern also from the polymorph sampler, and indices from the
BN structure; see SI). This transformation occurs at lower $x$ values for
samples grown in the approximately range 100°C $\lesssim
T_{\mathrm{dep}}\lesssim$ 225°C, with a large region of mixed phase RS and BN
(“RS + BN”). At a second threshold composition ($x$ $\gtrsim$ 0.8 at ambient
temperature, and at lower values for high temperatures) and up to 100% Zn at a
stoichiometry of approximately Zn3N2, a second phase transition occurs to the
BX phase. The presence of a BX phase in Zn3N2 corroborates literature
reports,[31] and may be enabled by ZrZn stabilization across phase space as
indicated by a shifts to higher $Q$ with increasing $x$. There are several
regions of phase-segregated BX phases as well, in particular at high
temperatures ($T_{\mathrm{dep}}>$ 250 °C) where films are likely completely
phase-separating into binaries of RS Zr3N4 and BX Zn3N2.

Overall, the measured and simulated XRD patterns correspond very well, except
for the offsets in $Q$ that are a consequence of errors in DFT lattice
constants and/or experimental artifacts (e.g. off-stoichiometry, possible
residual strain in the films, sample misalignment in XRD etc.). This work
warrants follow-up in-depth study of this abundant alloy space and associated
materials properties, but for this analysis we focus primarily on the near-
stoichiometric regime around ZnZrN2, i.e. $x$ = 0.5 (dotted vertical line in
Figure 2b), which can be most easily probed computationally due to its simple
stoichiometry and II-IV-V_2 chemistry.

### Statistical and thermodynamic influences on polymorph stabilization

Figure 3: (a) Thermodynamic density of states (TDOS) distribution of
polymorphs from the polymorph sampler using random structure searching,
demonstrating a dominant RS peak, with the crossed-out phases depicting the
absence of several predicted polymorph phases from (a). (b) Free energy as a
function of temperature, calculated from the polymorph sampler in (a).

Why are the RS-derived and BN-derived ZnZrN2 structures synthesized rather
than wurtsalt (WS), which is the predicted ground state at 0 K? To examine
this and contextualize our experimental findings, we apply a series of high-
to medium-throughput computational methods to investigate the role of cation
disorder and free energy on polymorph ordering. Figure 1(a) indicates that the
structures with lowest formation energies exhibit unique cation coordination
environments for Zr and Zn (WS, LC), while higher energy structures have
similar coordination environments for all cation sites (RS, BN, WZ, ZB). Our
hypothesis is that cation disordering during synthesis—enabled by rapid
condensation from the vapor state to the solid state in physical vapor
deposition (PVD)—favors structures with similar cation coordination
environments, thus lowering the probability of the formation of WS. To assess
this hypothesis, we (1) apply a polymorph sampler method to treat
configurational degrees of freedom, (2) develop a random disorder descriptor
to explain and understand the results of the polymorph sampler, and (3)
estimate formation energies of alloys with varied cation ratios to assess the
effects of off-stoichiometry.

#### Statistical treatment of configurational degrees of freedom

To account for the configurational degrees of freedom and the associated
entropy contributions to the free energy in various ZnZrN2 polymorph
structures, we apply the “polymorph sampler” statistical approach developed by
Stevanović and coworkers.[2, 32, 33] It has previously been shown that
treating the spectrum of structures generated by ab-initio random structure
sampling[24] as a proper statistical ensemble makes it possible to understand
the experimental realizability of metastable crystalline polymorphs phases[2,
32] as well as the structural features of glasses.[33] Namely, crystalline
structures with high ensemble probabilities — that is, high enough occurrence
in random sampling and low enough energies — were shown to correlate with
synthetically realized polymorphs in the chemical systems of MgO, ZnO, SnO2,
Si, and others. Experimental enthalpies, powder X-ray diffraction patterns,
and radial distribution functions of amorphous Si were also all reproduced
with high accuracy when compared to the corresponding computed ensemble-
averaged quantities.

Herein, the original approach is modified to include cation (lattice) disorder
in the ZnZrN2 system, the suspected phenomenon behind experimentally observed
structures. The modification pertains mainly to structure classification and
the statistical treatment (see SI for details). Specifically, we generate a
set of 5,000 random superlattice (RSL) ZnZrN2 structures with 24-atom cells
and relax them to the closest local minimum using DFT with the GGA+U
functional (note that SCAN is too computationally expensive for this many
structures). The resulting spectrum of polymorph structures is split into
prototype classes with the same underlying space group (see SI). Next, the
ensemble probability $P_{k}$ of every individual class is evaluated as:

$P_{k}=\frac{\sum_{n=m}^{m+n_{k}}\omega_{n}\mathord{\mathrm{e}}^{-\frac{\epsilon_{n}}{k_{\mathrm{B}}T}}}{\Xi}=\frac{\Xi_{k}}{\Xi}$
(1)

where $k$ represents different prototype classes, $n$ counts different
polymorph structures within a given prototype class, $\omega_{n}$ is the
frequency of occurrence of a structure $n$ belonging to the class $k$, and
$\epsilon_{n}$ is the total energy normalized per atom evaluated relative to
the ground state structure. Lastly, $\Xi_{k}$ and $\Xi$ are the partial and
the full partition functions respectively; the former pertaining only to the
$k$-class of structures with the same underlying space group, while the latter
is evaluated for the entire set of structures obtained from the random
sampling.

The thermodynamics density of states (TDOS, i.e. number of structures per
energy unit) resulting from the random structure sampling, normalized and
resolved by the underlying structure class, is shown in Figure 3(a). Two
features become immediately apparent. First, consistent with the KLM algorithm
and literature reports, the ground state WS structure is correctly identified,
but the corresponding frequency of occurrence of this structure and the
associated TDOS are very small. Second, random sampling indicates that the
class of RS structures has the largest TDOS peak (blue), indicating the
flexibility of this structure class to accommodate cation disorder in a
relatively narrow energy interval. While their energies span $\sim$0.04 –
$\sim$0.15 eV/atom range, the largest fraction of structures is concentrated
in the narrow window of 0.08–0.12 eV/atom, giving rise to the large peak shown
in Figure 3(a). Additional classes of disordered structures with more
significant occurrence in the random polymorph sampling are the BN (red) and
the two highly distorted d-RS and d-WS classes (striped patterns), with the
latter also including structures with more than one tetrahedral layer
sandwiched between the octahedrally-coordinated layers of the WS structures.
It is also important to note that none of the other previously discussed
structure types appear in the random sampling, including the LC and WZ
structure classes with T = 0 enthalpies lower than that of RS. This indicates
that all these possible structures exhibit very “narrow” local minima in
configurational energy space, leading to the very low probability of
occurrence. The same is true for a relatively large number of very low
symmetry structures (space groups $P1$ and $P\bar{1}$) typically obtained in
random structure sampling, but since none of these structures occur in large
numbers they become statistically insignificant compared to those depicted in
Figure 3.

The TDOS from Figure 3(a) allows evaluating ensemble probabilities per
Equation 1 and associated partial partition functions $\Xi_{k}$. These can
then be used to evaluate the “configurational” free energies $F_{k}$ of the
corresponding structure types using the standard statistical mechanics
equation:

$F_{k}=-k_{\mathrm{B}}T\ln\Xi_{k}$ (2)

where $k_{\mathrm{B}}$ is the Boltzmann constant and $T$ is the temperature.
The temperature dependence of the free energies displayed in Figure 3(b)
clearly shows that at low $T$ the lowest free energy structure is the ground
state WS structure, consistent with Figure 1; at $T\approx$ 1100 K the
disordered RS becomes the most favorable structure due to the large gain in
configurational entropy at higher temperatures.[34] In the temperature range
1300–1600 K, the WS structure gives way to the disordered BN as the second
most favorable structure, the order that perseveres, while at higher
temperatures the d-RS takes over. This structure, if mixed with RS, would be
virtually experimentally undetectable from RS because of their very similar
XRD patterns. In summary, at higher temperatures the polymorph sampler
procedure and the ensemble treatment suggest the following ordering of
structures according to their free energies: (1) disordered RS, (2) distorted
RS (d-RS), (3) BN, and (4) distorted WS (d-WS). This is consistent with
experimentally observed XRD patterns that are compared with the ensemble-
averaged patterns in Figure 2(c). Since previous studies have suggested that
sputter deposition occurs at effective temperatures greater than 1150 K in
ternary nitrides (note that this is not equivalent to synthesis temperature
$T_{\mathrm{dep}}$, which scales inversely with effective temperature for
$T_{\mathrm{dep}}\lesssim$ 600°C; see SI),[35, 36] these ensemble
probabilities and free energies support the observed behavior, i.e. the
stabilization of the RS phase and the destabilization of the WS phase.

Figure 4: (a) Schematic of the WS and RS derived structures with cation site
disorder. WS is an example of a structure with in-equivalent cation sites,
while RS is an example of a structure with equivalent cation sites. (b)
Comparison of formation energies of the lowest ordered polymorph structure and
the random disordered structure, derived from special quasirandom structure
(SQS) calculations, for each prototype class that emerged from the polymorph
sampler (see Figure 3). The $E_{\mathrm{rdp}}$ descriptor is denoted by the
magnitude of the dashed arrow, depicting the lower disorder tolerance for the
BN and RS structures that are observed experimentally (see SI for all
structure energies and details).

#### Configurational enthalpy of random disorder

What is the physical rationale behind the higher cation disorder tolerance in
RS ZnZrN2, and why do we also observe BN experimentally in the Zn-Zr-N phase
space? In practice, cation disorder is ubiquitous in ternary nitrides,[37, 36,
38] especially in materials synthesized at high effective temperatures (as
present in sputtering). Thus, modeling small ordered unit cells may not
adequately capture energetic information in these systems. Here, we explicitly
consider the energetic effects of random disorder in ZnZrN2 polymorph
structures. One common method to simulate disorder is the special quasirandom
structures (SQS) method, which models random atomic decorations on a lattice
in relatively small unit cells.[39] This is achieved by searching for unit
cells that reproduce or approximate pair (or higher order) correlation
functions up to a cutoff distance by minimizing an objective function (a
metric to assess the “goodness” of the disordered structure with respect to
the correlation functions). Here, we calculate a set of SQS structures with 64
atoms for each structure class using the ATAT package,[40, 41] selecting only
structures with the lowest objective functions (see SI and the python module
pymatgen.command_line.mcsqs_caller for equations and details). Each SQS
structure is relaxed, its formation enthalpy is computed with GGA and then
SCAN, and it is assigned to its closest structure prototype class via a
structure-matching algorithm (see SI) to account for any SQS structures that
relaxed to a different geometry. For each structure class $k$, the resulting
SQS formation energies all are within $\sim$10 eV of one another; these
energies are then ensemble-averaged to best represent the random disorder
formation energy, $\Delta H_{k}^{\mathrm{disordered}}$.

Formation energies of random disordered structures are presented in Figure
4(b) with thick bar markers, in comparison to those of the lowest energy
ordered structures of each prototype (thin bar markers), plotting only the
five structure classes that emerge from the polymorph sampler (see SI for
other structure classes). We find that the random disordered WS and d-WS
structures have significantly higher formation energies than their ordered
counterparts, with the disordered WS energy increasing drastically to over
0.170 eV/atom. In contrast, the BN random disordered structure has the lowest
formation energy of all disordered structures. Although RS does not have the
lowest random energy disordered structure, it is significant is that the
formation energy of the RS disordered structure is very low (within 25
meV/atom) compared to that of the RS ordered structure. To assess this
disorder tolerance for each class $k$, we introduce a new descriptor, the
“random disordered polymorph energy” $E_{{\mathrm{rdp},k}}$:

$E_{{\mathrm{rdp},k}}=\Delta H_{k}^{\mathrm{disordered}}-\Delta
H_{k}^{\mathrm{ordered}}$ (3)

which is the difference between the random disordered SQS structure’s energy
and the energy of the lowest energy ordered structure $\Delta
H_{k}^{\mathrm{ordered}}$, for a given structure class. RS and BN, which are
the structures that have been experimentally synthesized, have the lowest
$E_{\mathrm{rdp}}$ values. Physically, since the $T=0$ DFT formation energy is
an approximation of formation enthalpy, the $E_{\mathrm{rdp}}$ corresponds to
the additional configurational enthalpy that is introduced for each structure
as a result of cation disorder. Thus, we have shown that in the ZnZrN2
polymorph structures with in-equivalent cation sites (WS, LC), disordering
significantly increases enthalpy, whereas in the ZnZrN2 structures with
equivalent cation sites (RS, BN) disordering only negligibly increases
enthalpy.

Furthermore, since random disordered SQS structures effectively represent high
temperature structures, by qualitatively comparing the formation energies
$\Delta H_{k}^{\mathrm{disordered}}$ of the random disordered structures in
Figure 4(b) to the statistical ensemble free energies $F_{k}$ from Figure
3(b), we observe that the relative polymorph ordering changes: the RS phase is
far lower in free energy $F_{k}$ at high temperatures (T $\gtrsim$ 1000 K),
while the disordered BN phase is lowest in $\Delta
H_{k}^{\mathrm{disordered}}$ in Figure 4(b). Since $F_{k}$ inherently includes
configurational enthalpy and configurational entropy, and $\Delta
H_{k}^{\mathrm{disordered}}$ includes only configurational enthalpy, this
suggests that the WS phase is enthalpically destabilized, while the RS phase
is entropically stabilized compared to the BN and other phases. We acknowledge
that these SQS disordered structures are estimates (a cluster expansion should
be fit to rigorously account for local effects of cation ordering), but
ultimately these estimates support the hypothesis that disorder tolerance
stabilizes the RS and BN phases over other considered polymorphs, due to both
configurational enthalpy and configurational entropy.

Figure 5: (a) Ternary phase diagram of Zn-Zr-N and a close-up of the
Zn_$x$Zr_1-$x$N region, with computed stable and unstable compounds from the
Materials Project designated. The colored circle depict the approximate
composition and associated phases of experimental data from this work. (b)
Computed phase space across the alloy tie-line ZrN–Zn3ZrN4, with lowest
energies plotted for each computed structure. Similar structure prototypes are
connected with dotted lines for the lowest energy structures. The thick line
on the bottom indicates the approximate hull.

#### Off-stoichiometry tolerance

Lastly, we have thus far suggested why RS and BN are stabilized, but why is RS
generally synthesized at Zr-rich compositions and BN at Zn-rich compositions?
Another plausible explanation for the absence of WS ZnZrN2 is that it is a
line compound, i.e. stable only in a very narrow region of ternary phase space
(see Figure 5(a)) and configurational energy space that we may have missed
with combinatorial growth. To address this, we use a simplified alloy model.
We perform cation substitution in each of the ordered ZnZrN2 polymorph
structures from Figure 1 to create a set of ordered prototypes across the
ZrN–ZnN tie-line (i.e. Zn_$x$Zr_1-$x$N_$y$ where $y$ = 1 and $x$ = 0, 0.25,
0.50, and 0.75), and then relax the structures and calculate their formation
energies (with SCAN). We note that this is a very rough approximation, since
films are N-rich for $x$ < 0.5 and N-poor for $x$ > 0.5; the experimental
alloy is closer to the Zr3N4–Zn3N2 tie-line, but this heterovalent alloy is
far trickier to model due to defect compensation and is beyond the scope here.
Figure 5(b) plots the relative phase stability for four values of $x$ across
this tie-line, normalized to the $E_{\mathrm{hull}}$ for a given composition,
with horizontal bar markers denoting the lowest energy polymorph of a given
prototype and dotted lines connecting such values (note some absent phases,
since d-RS and d-HX both relax into RS for $x$ = 0, 0.25, and WZ relaxes into
BN for $x\neq$ 0.5). This indeed suggests a narrow window of $x$ in which WS
has the lowest $E_{\mathrm{hull}}$, with RS the lowest energy polymorph for a
wide window of approximately 0 < $x$ $\lesssim$ 0.45, and BN (or d-RS) the
lowest energy polymorph at approximately $x$ $\gtrsim$ 0.55. Removing WS and
LC, following the polymorph sampler results, would suggests a phase change
from RS to BN somewhere within approximately 0.45 < $x$ < 0.55. We do not plot
$x$ > 0.75 (labeled “non-physical regime”), as ZnN is a non-physical unstable
compound since Zn is not stable in a +3 oxidation state, yet another phase
change at high Zn content is plausible since we observe Zn3N2 experimentally.
We note that in our alloy system Zn_$x$Zr_1-$x$N_$y$, the presence of a lower
density hexagonal phase (BN, here) located between two higher density cubic
phases (RS and BX, here) is indicative of a phenomenon in heterovalent
heterostructural alloys called “negative pressure” alloys.[42, 43]
Additionally, these calculations use the nominal valence of the cations,
namely, Zn2+ and Zr4+; we do not perform defect calculations nor vary cation
oxidation states. Rigorous examination of alloy phase space would require an
in-depth calculation of a temperature phase diagram, which is beyond our
scope, but this approximation supports our experimental observation of a phase
change from RS to BN as $x$ increases in Zn_$x$Zr_1-$x$N_$y$.

#### Other effects and implications

In principle, it is also necessary to acknowledge the role of dynamic and
kinetic effects in this ternary phase space. Dynamic instabilities could be
assessed by phonon calculations. As has been demonstrated in other II-IV-N2
systems, spurious oxygen incorporation from the growth chamber can influence
phase stability and result in impurity phases phases,[14] though our values
below 1% from Figure 2(a) suggest that a phase-segregated oxide is not
observed here. Additionally, it is important to note that the elemental Zr
used in synthesis contains $\sim$10 at. % Hf, which could influence relative
phase stability. Finally, surface morphology, templating during growth, and
kinetic effects could be assisting in restricting the formation of WS, and in
enabling the preferential formation of RS. Even though these films are grown
on amorphous fused silica, we also acknowledge the possibility of preferential
nucleation.

As a thought experiment to assess the role of disorder tolerance in other new
ternary pnictide systems, we can consider the set of II-IV-V_2 pnictide
semiconductors where II $=$ Ba, Mg, Ca, Be, Sr, Zn, Cd, IV $=$ Si, Ge, Sn, Pb,
Ti, Zr, Hf, and V $=$ N, P, As. Theoretically, this set includes 147 unique
compositions, though only 31 of these compositions have been confirmed
experimentally and only 43 other predicted compositions are in the Materials
Project (MP) database, leaving 73 II-IV-V_2 compositions not on MP. Within the
set of predicted but not-yet-synthesized compounds, six other systems other
than ZnZrN2 have computed ground states where cations occupy symmetrically in-
equivalent lattice sites: SrGeN2, CaTiN2, ZnHfN2, SrZrN2, SrHfN2, and CdGeN2.
Thus, it is feasible that higher-enthalpy polymorphs of these compounds, as
well as the 73 empty compositions, could be synthesized, leading to the
experimental realizability of new metastable compounds with promising
properties to be explored.

## Conclusion

In summary, we have grown the first Zn_$x$Zr_1-$x$N_$y$ alloys (a set of
approximately 400 samples), demonstrating the crystallization of metastable
rocksalt-derived (RS) and boron-nitride-derived (BN) phases rather than the
predicted “wurtsalt” (WS) ground state phase. We explain these findings using
a series of first-principles computational methods. First, through a
statistically rigorous polymorph sampler approach, we demonstrate energetic
destabilization of the predicted ground state at high temperatures, and
stabilization of RS and BN phases that support our experimental results.
Second, we propose that this stabilization is attributed to the increased
disorder tolerance in the RS and BN phases due to only minor gains in
configurational enthalpy, and suggest that the RS phase is more entropically
stabilized than BN. Third, ordered alloy calculations suggest that RS and BN
have a higher tolerance to cation off-stoichiometries compared to WS, and we
predict a phase transformation from RS to BN as $x$ increases that
corroborates our experimental findings. We acknowledge that growth conditions
matter significantly in phase stabilization, and have demonstrated the
advantages of combinatorial sputtering in accessing high effective
temperatures and quenching metastable states. Although WS is energetically
destabilized here, its realizability is not ruled out, and this case study
could provide a framework to target synthesis for WS phases (e.g. low
effective temperature, epitaxial, on-stoichiometry synthesis). The
optoelectronic properties of the synthesizable polymorphs in this system
remain to be investigated, and Table 1 indicates promising properties for
device applications such as contact materials, solar cell absorbers,
photocatalysts, piezoelectric and ferroelectric materials.[16, 44] In
particular, the synthesized BN-derived polymorph has $>$2 eV band gap and low
(<1.5) well-matched electron and hole effective masses, making it interesting
for electronic devices that can operate at elevated temperatures.
Additionally, this non-polar BN polymorph is the transition state between two
variants of the polar WZ structure, suggesting a pathway to tuning its
predicted and measured ferroelectric response.[45, 46]

A fundamental goal of materials discovery is to bridge the gap between
computationally predicted materials and materials that can actually be grown
in the laboratory with desired properties. Thus, these results are interesting
beyond this specific material system; this study points to the fact that there
may be many accessible energetic states that neither we nor nature have
realized yet, suggests that some of these accessible states may be
thermodynamically “metastable” (at $T=0$), and offers possible routes to
synthesis. Many high-throughput materials discovery studies discard materials
that have an $E_{\mathrm{hull}}$ above a certain cutoff. However, since we
have demonstrated synthesis of a polymorph phase with $E_{\mathrm{hull}}>$
0.080 eV/atom using a common physical vapor deposition technique, we must
question what other metastable phases we have been discarding. In the Zn-Zr-N
system, we have suggested tolerance to disorder and off-stoichiometry are
responsible for the realization of metastable phases, and have introduced a
new descriptor to assess disorder tolerance, but in other systems there may be
different physical mechanisms enabling synthesis of high-energy polymorphs. In
general, the materials discovery community needs to continue to redefine the
metrics by which we computationally assess phase stability in order to yield
experimentally realizable predictions that enable new functional materials.

## Acknowledgments

This work was authored in part at the National Renewable Energy Laboratory,
operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of
Energy (DOE) under Contract No. DE-AC36-08GO28308. The funding was provided by
the Office of Science (SC), Office of Basic Energy Sciences (BES), Materials
Chemistry program, as a part of the Early Career Award “Kinetic Synthesis of
Metastable Nitrides”. Rachel Woods-Robinson acknowledges financial support
from the U.C. Berkeley Chancellor’s Fellowship and the National Science
Foundation (NSF) Graduate Research Fellowship under Grant No. DGE1106400 and
DGE175814. Dr. Vladan Stevanović acknowledges financial support from NSF
Career Award No. DMR-1945010. Use of the Stanford Synchrotron Radiation
Lightsource, SLAC National Accelerator Laboratory, is supported by DOE’s
Office of Science (SC), Basic Energy Sciences (BES) under Contract No. DE-
AC02-76SF00515. The authors thank Dr. Sage Bauers, Dr. Kevin Talley, Valerie
Jacobson, and Rachel Sherbondy for experimental assistance, Dr. Matthew K.
Horton, Dr. Shyam Dwaraknath, and Matthew McDermott for assistance with
calculations, Dr. John Perkins and Dr. Apurva Mehta with characterization
assistance, and Dr. Adele Tamboli and Dr. Wenhao Sun for insightful
discussions. This work used high-performance computing resources located at
NREL and sponsored by the Office of Energy Efficiency and Renewable Energy.
The views expressed in the article do not necessarily represent the views of
the DOE or the U.S. Government.

## Supplementary materials

The following supporting information is included:

  * •

S1 Synthesis details

  * •

S2 Composition and structural characterization

  * •

S3 Polymorph and DFT ordering

  * •

S4 Computational methods

## Author contributions

Conceptualization, R.W.R., A.Z., K.A.P.; Methodology, R.W.R., V.S., S.L.,
K.H.; Computational investigation, R.W.R., V.S., S.L.; Experimental
investigation, R.W.R., K.H., A.Z.; Writing, R.W.R., V.S.; Funding Acquisition,
R.W.R., A.Z., V.S., K.A.P.; Supervision, A.Z., V.S., K.A.P.

## References

  * [1] W. Sun, A. Holder, B. Orvañanos, E. Arca, A. Zakutayev, S. Lany, and G. Ceder, “Thermodynamic routes to novel metastable nitrogen-rich nitrides,” Chemistry of Materials, vol. 29, no. 16, pp. 6936–6946, 2017\.
  * [2] V. Stevanović, “Sampling polymorphs of ionic solids using random superlattices,” Physical Review Letters, vol. 116, no. 7, p. 075503, 2016\.
  * [3] W. Sun, S. T. Dacek, S. P. Ong, G. Hautier, A. Jain, W. D. Richards, A. C. Gamst, K. A. Persson, and G. Ceder, “The thermodynamic scale of inorganic crystalline metastability,” Science Advances, vol. 2, no. 11, p. e1600225, 2016.
  * [4] M. Aykol, S. S. Dwaraknath, W. Sun, and K. A. Persson, “Thermodynamic limit for synthesis of metastable inorganic materials,” Science Advances, vol. 4, no. 4, p. eaaq0148, 2018.
  * [5] M. Aykol, V. I. Hegde, L. Hung, S. Suram, P. Herring, C. Wolverton, and J. S. Hummelshøj, “Network analysis of synthesizable materials discovery,” Nature Communications, vol. 10, no. 1, pp. 1–7, 2019.
  * [6] M. Horton, S. Dwaraknath, and K. A. Persson, “Perils and promise of materials databases,” Nature Computational Science, vol. 1, no. 1, 2021.
  * [7] P. Kroll, “Pathways to metastable nitride structures,” Journal of Solid State Chemistry, vol. 176, no. 2, pp. 530–537, 2003.
  * [8] A. L. Greenaway, C. L. Melamed, M. B. Tellekamp, R. Woods-Robinson, E. S. Toberer, J. R. Neilson, and A. C. Tamboli, “Ternary nitride materials: Fundamentals and emerging device applications,” 2020.
  * [9] Y. Hinuma, T. Hatakeyama, Y. Kumagai, L. A. Burton, H. Sato, Y. Muraba, S. Iimura, H. Hiramatsu, I. Tanaka, H. Hosono, and F. Oba, “Discovery of earth-abundant nitride semiconductors by computational screening and high-pressure synthesis,” Nature Communications, vol. 7, no. 1, p. 11962, 2016.
  * [10] W. Sun, C. J. Bartel, E. Arca, S. R. Bauers, B. Matthews, B. Orvañanos, B.-R. Chen, M. F. Toney, L. T. Schelhas, W. Tumas, J. Tate, A. Zakutayev, S. Lany, A. M. Holder, and G. Ceder, “A map of the inorganic ternary metal nitrides,” Nature Materials, vol. 18, pp. 732–739, July 2019.
  * [11] A. D. Martinez, A. N. Fioretti, E. S. Toberer, and A. C. Tamboli, “Synthesis, structure, and optoelectronic properties of II-IV-V_2 materials,” Journal of Materials Chemistry A, vol. 5, no. 23, pp. 11418–11435, 2017.
  * [12] S. R. Bauers, A. Holder, W. Sun, C. L. Melamed, R. Woods-Robinson, J. Mangum, J. Perkins, W. Tumas, B. Gorman, A. Tamboli, et al., “Ternary nitride semiconductors in the rocksalt crystal structure,” Proceedings of the National Academy of Sciences, vol. 116, no. 30, pp. 14829–14834, 2019.
  * [13] F. Kawamura, M. Imura, H. Murata, N. Yamada, and T. Taniguchi, “Synthesis of a Novel Rocksalt-Type Ternary Nitride Semiconductor MgSnN2 Using the Metathesis Reaction Under High Pressure,” European Journal of Inorganic Chemistry, vol. 2020, no. 5, pp. 446–451, 2020.
  * [14] A. L. Greenaway, A. L. Loutris, K. N. Heinselman, C. L. Melamed, R. R. Schnepf, M. B. Tellekamp, R. Woods-Robinson, R. Sherbondy, D. Bardgett, S. Bauers, et al., “Combinatorial synthesis of magnesium tin nitride semiconductors,” Journal of the American Chemical Society, vol. 142, no. 18, pp. 8421–8430, 2020.
  * [15] Y. Hinuma, T. Hatakeyama, Y. Kumagai, L. A. Burton, H. Sato, Y. Muraba, S. Iimura, H. Hiramatsu, I. Tanaka, H. Hosono, et al., “Discovery of earth-abundant nitride semiconductors by computational screening and high-pressure synthesis,” Nature Communications, vol. 7, no. 1, pp. 1–10, 2016.
  * [16] C. Tholander, C. Andersson, R. Armiento, F. Tasnadi, and B. Alling, “Strong piezoelectric response in stable TiZnN2, ZrZnN2, and HfZnN2 found by ab initio high-throughput approach,” Journal of Applied Physics, vol. 120, no. 22, p. 225102, 2016.
  * [17] W. Sun, C. J. Bartel, E. Arca, S. R. Bauers, B. Matthews, B. Orvañanos, B.-R. Chen, M. F. Toney, L. T. Schelhas, W. Tumas, et al., “A map of the inorganic ternary metal nitrides,” Nature Materials, vol. 18, no. 7, p. 732, 2019.
  * [18] J. Dismukes, R. Smith, and J. White, “Physical properties and crystal structure of a new semiconducting I-III-VI_2 compound, CuScS2,” Journal of Physics and Chemistry of Solids, vol. 32, no. 4, pp. 913–922, 1971.
  * [19] D. O. Scanlon and G. W. Watson, “Stability, geometry, and electronic structure of an alternative I-III-VI_2 material, CuScS2: A hybrid density functional theory analysis,” Applied Physics Letters, vol. 97, no. 13, p. 131904, 2010.
  * [20] Materials Project, “Materials data on zrznn2 by the Materials Project,” 5 2017.
  * [21] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, et al., “Commentary: The Materials Project: A materials genome approach to accelerating materials innovation,” Apl Materials, vol. 1, no. 1, p. 011002, 2013.
  * [22] S. M. Woodley and R. Catlow, “Crystal structure prediction from first principles,” Nature Materials, vol. 7, no. 12, pp. 937–946, 2008.
  * [23] P. P. Zawadzki, J. Perkins, and S. Lany, “Modeling amorphous thin films: Kinetically limited minimization,” Physical Review B, vol. 90, no. 9, p. 094203, 2014.
  * [24] C. J. Pickard and R. Needs, “Ab initio random structure searching,” Journal of Physics: Condensed Matter, vol. 23, no. 5, p. 053201, 2011.
  * [25] A. R. Oganov, A. O. Lyakhov, and M. Valle, “How evolutionary crystal structure prediction works and why,” Accounts of chemical research, vol. 44, no. 3, pp. 227–237, 2011.
  * [26] V. Stevanović, S. Lany, X. Zhang, and A. Zunger, “Correcting density functional theory for accurate predictions of compound enthalpies of formation: Fitted elemental-phase reference energies,” Physical Review B, vol. 85, no. 11, p. 115104, 2012.
  * [27] A. P. Bartók and J. R. Yates, “Regularized SCAN functional,” The Journal of chemical physics, vol. 150, no. 16, p. 161101, 2019.
  * [28] J. H. Yang, D. A. Kitchaev, and G. Ceder, “Rationalizing accurate structure prediction in the meta-GGA SCAN functional,” Physical Review B, vol. 100, no. 3, p. 035132, 2019.
  * [29] M. L. Green, I. Takeuchi, and J. R. Hattrick-Simpers, “Applications of high throughput (combinatorial) methodologies to electronic, magnetic, optical, and energy-related materials,” Journal of Applied Physics, vol. 113, no. 23, p. 231101, 2013.
  * [30] S. R. Bauers, D. M. Hamann, A. Patterson, J. D. Perkins, K. R. Talley, and A. Zakutayev, “Composition, structure, and semiconducting properties of Mg_xZr_2-xN_2 thin films,” Japanese Journal of Applied Physics, vol. 58, no. SC, p. SC1015, 2019.
  * [31] D. Partin, D. Williams, and M. O’Keeffe, “The crystal structures of Mg3N2 and Zn3N2,” Journal of Solid State Chemistry, vol. 132, no. 1, pp. 56–59, 1997.
  * [32] E. B. Jones and V. Stevanović, “Polymorphism in elemental silicon: Probabilistic interpretation of the realizability of metastable structures,” Physical Review B, vol. 96, no. 18, p. 184101, 2017.
  * [33] E. B. Jones and V. Stevanović, “The glassy solid as a statistical ensemble of crystalline microstates,” npj Computational Materials, vol. 6, no. 1, pp. 1–6, 2020.
  * [34] C. M. Rost, E. Sachet, T. Borman, A. Moballegh, E. C. Dickey, D. Hou, J. L. Jones, S. Curtarolo, and J.-P. Maria, “Entropy-stabilized oxides,” Nature communications, vol. 6, no. 1, pp. 1–8, 2015.
  * [35] A. N. Fioretti, J. Pan, B. R. Ortiz, C. L. Melamed, P. C. Dippo, L. T. Schelhas, J. D. Perkins, D. Kuciauskas, S. Lany, A. Zakutayev, et al., “Exciton photoluminescence and benign defect complex formation in zinc tin nitride,” Materials Horizons, vol. 5, no. 5, pp. 823–830, 2018.
  * [36] S. Lany, A. N. Fioretti, P. P. Zawadzki, L. T. Schelhas, E. S. Toberer, A. Zakutayev, and A. C. Tamboli, “Monte Carlo simulations of disorder in ZnSnN2 and the effects on the electronic structure,” Physical Review Materials, vol. 1, no. 3, p. 035401, 2017.
  * [37] P. C. Quayle, E. W. Blanton, A. Punya, G. T. Junno, K. He, L. Han, H. Zhao, J. Shan, W. R. Lambrecht, and K. Kash, “Charge-neutral disorder and polytypes in heterovalent wurtzite-based ternary semiconductors: The importance of the octet rule,” Physical Review B, vol. 91, no. 20, p. 205207, 2015.
  * [38] R. R. Schnepf, J. J. Cordell, M. B. Tellekamp, C. L. Melamed, A. L. Greenaway, A. Mis, G. L. Brennecka, S. Christensen, G. J. Tucker, E. S. Toberer, et al., “Utilizing site disorder in the development of new energy-relevant semiconductors,” ACS Energy Letters, vol. 5, no. 6, pp. 2027–2041, 2020\.
  * [39] A. Zunger, S.-H. Wei, L. Ferreira, and J. E. Bernard, “Special quasirandom structures,” Physical Review Letters, vol. 65, no. 3, p. 353, 1990.
  * [40] A. Van De Walle, “Multicomponent multisublattice alloys, nonconfigurational entropy and other additions to the Alloy Theoretic Automated Toolkit,” Calphad, vol. 33, no. 2, pp. 266–278, 2009.
  * [41] A. Van de Walle, P. Tiwary, M. De Jong, D. Olmsted, M. Asta, A. Dick, D. Shin, Y. Wang, L.-Q. Chen, and Z.-K. Liu, “Efficient stochastic generation of special quasirandom structures,” Calphad, vol. 42, pp. 13–18, 2013.
  * [42] S. Siol, A. Holder, J. Steffes, L. T. Schelhas, K. H. Stone, L. Garten, J. D. Perkins, P. A. Parilla, M. F. Toney, B. D. Huey, et al., “Negative-pressure polymorphs made by heterostructural alloying,” Science Advances, vol. 4, no. 4, p. eaaq1442, 2018.
  * [43] R. Woods-Robinson, Y. Han, J. S. Mangum, C. L. Melamed, B. P. Gorman, A. Mehta, K. A. Persson, and A. Zakutayev, “Combinatorial tuning of structural and optoelectronic properties in Cu_xZn_1-xS,” Matter, vol. 1, no. 4, pp. 862–880, 2019.
  * [44] H. Ling, S. S. Dwaraknath, and K. A. Persson, “Origin of disorder tolerance in piezoelectric materials and design of polar systems,” Chemistry of Materials, vol. 32, no. 7, pp. 2836–2842, 2020.
  * [45] C. E. Dreyer, A. Janotti, C. G. Van de Walle, and D. Vanderbilt, “Correct implementation of polarization constants in wurtzite materials and impact on iii-nitrides,” Physical Review X, vol. 6, no. 2, p. 021038, 2016.
  * [46] S. Fichtner, N. Wolff, F. Lofink, L. Kienle, and B. Wagner, “Alscn: A iii-v semiconductor based ferroelectric,” Journal of Applied Physics, vol. 125, no. 11, p. 114103, 2019.