# Accelerated Discovery of a Large Family of Quaternary Chalcogenides with
Very Low Lattice Thermal Conductivity

Koushik Pal<EMAIL_ADDRESS>Department of Materials Science and
Engineering, Northwestern University, Evanston, IL 60208, USA  Yi Xia
Department of Materials Science and Engineering, Northwestern University,
Evanston, IL 60208, USA  Jiahong Shen  Department of Materials Science and
Engineering, Northwestern University, Evanston, IL 60208, USA  Jiangang He
Department of Materials Science and Engineering, Northwestern University,
Evanston, IL 60208, USA  Yubo Luo  Department of Chemistry, Northwestern
University, Evanston, IL 60208, USA  Mercouri G. Kanatzidis  Department of
Chemistry, Northwestern University, Evanston, IL 60208, USA  Chris Wolverton
<EMAIL_ADDRESS>Department of Materials Science and Engineering,
Northwestern University, Evanston, IL 60208, USA

###### Abstract

The development of efficient thermal energy management devices such as
thermoelectrics, barrier coatings, and thermal data-storage disks often relies
on compounds that possess very low lattice thermal conductivity
($\kappa_{l}$). Here, we present the computational prediction of a large
family of 628 thermodynamically stable quaternary chalcogenides, AMM’Q3 (A =
alkali/alkaline earth/post-transition metals; M/M’ = transition metals,
lanthanides; Q = chalcogens) using high-throughput density functional theory
(DFT) calculations. We validate the presence of low-$\kappa_{l}$ in this
family of materials by calculating $\kappa_{l}$ of several predicted stable
compounds using the Peierls-Boltzmann transport equation within a first-
principles DFT framework. Our analysis reveals that the low-$\kappa_{l}$ in
the AMM’Q3 compounds originates from the presence of either a strong lattice
anharmonicity that enhances the phonon scatterings or rattlers cations that
lead to multiple scattering channels in their crystal structures. Our
predictions suggest new experimental research opportunities in the synthesis
and characterization of these stable, low-$\kappa_{l}$ compounds.

##  I Introduction

An important focus in materials science research has been to discover novel
materials with properties that might hold the keys to solving the most
pressing problems in renewable energy, energy harvesting, or semiconductor
power electronics. The augmentation of new materials discovery and the
prediction of their properties have been accelerated by the advent of advanced
computer algorithms coupled with high-throughput (HT) screening
methodsCurtarolo  _et al._ (2013); Mounet  _et al._ (2018); Saal  _et al._
(2013); Kirklin  _et al._ (2015); Jain  _et al._ (2013, 2016); Meng and
Arroyo-de Dompablo (2009); Greeley  _et al._ (2006); Gautier  _et al._ (2015);
Zakutayev  _et al._ (2013) using accurate quantum mechanical calculations
based on density functional theory (DFT). In the recent past, several
computational predictions have led to the successful synthesis of new solid-
state compounds in a variety of chemistry and structure types in the family of
half-HeuslersGautier  _et al._ (2015); Zakutayev  _et al._ (2013), double
half-HeuslersAnand  _et al._ (2019), electridesWang  _et al._ (2017), AB2X4
based chalcogendesXi  _et al._ (2018), and rocksalt based compoundsHao  _et
al._ (2014).

Crystalline solids with extreme thermal transport properties are
technologically important for the efficient management of thermal
energySamanta  _et al._ (2020). While materials with high lattice thermal
conductivity ($\kappa_{l}$) are used in microelectronic devices for heat
dissipation, materials with low $\kappa_{l}$ are used in thermal barrier
coatingsDarolia (2013), thermal data-storage devicesMatsunaga  _et al._
(2011), and high-performance thermoelectrics (TEs)Zhao  _et al._ (2016);
Biswas  _et al._ (2012) which can convert heat into electrical energy. The
conversion efficiency of the TEs is defined by the figure-of-merit,
ZT=$S^{2}\sigma/(\kappa_{l}+\kappa_{e})$, where S, $\sigma$, $\kappa_{e}$ and
$\kappa_{l}$ are the Seebeck coefficient, electrical conductivity, electronic
thermal conductivity, and lattice thermal conductivity, respectively.
Engineering the electronic band structures of crystalline compounds that
already possess low $\kappa_{l}$ has emerged to be a very popular strategy to
increase the ZT. Therefore, crystalline semiconductors with intrinsically
low-$\kappa_{l}$ are highly sought after in thermoelectrics and other thermal
energy management devices.

In pursuit of finding new low-$\kappa_{l}$ materials, different classes of
crystalline compounds like the perovskitesvan Roekeghem  _et al._ (2016),
half-HeuslersCarrete  _et al._ (2014), full-HeuslersHe  _et al._ (2016), and
double half-HeuslersAnand  _et al._ (2019), have been explored using HT
computational methods, and subsequently, some of the predicted compounds were
experimentally synthesized successfullyGautier  _et al._ (2015); Zakutayev
_et al._ (2013); Anand  _et al._ (2019). Despite the existence of a large
number of crystalline compounds in various materials databases e.g., the
Inorganic Crystal Structure Database (ICSD)Belsky  _et al._ (2002), Open
Quantum Materials Database (OQMD)Saal  _et al._ (2013); Kirklin  _et al._
(2015), Materials ProjectJain  _et al._ (2013), and AflowlibCurtarolo  _et
al._ (2012), it is quite important to look for novel stable and metastable
compounds which might exhibit exciting physical and chemical properties. In
this work, we present the computational discovery of a large number of stable
(628) and low-energy metastable (852) quaternary chalcogenides AMM’Q3 (A =
alkali, alkaline earth, post-transition metals; M/M’ = transition metals,
lanthanides; Q = chalcogens) that span a huge chemical space across the
periodic table. Our results are based on reliable, accurate and robust HT-DFT
calculations where (a) we generated initial AMM’Q3 compositions following the
experimentally known AMM’Q3 compounds formation criteria, (b) calculated their
energetics in all known seven crystallographic prototypes that are found in
this materials family, and (c) performed thermodynamic phase stability
analysis of these compositions against all possible competing phases that are
present in the OQMD.

Figure 1: Crystallographic prototypes of the AMM’Q3 compounds. (a-g) Seven
structural prototypes in the family of experimentally known AMM’Q3 compounds.
Most of them have layered crystal structures where the layers of Am+ cations
are sandwiched between the covalently bonded layers of [MM’Q3]m- anion
sublattice, and interact with the later through electrostatic interactions.
The black solid lines in each figure indicate the conventional unit cell of
each compound. (h) Distributions of the experimentally known AMM’Q3 compounds
in the seven structure types.

About 192 quaternary chalcogenides (see Supplementary Information (SI) for a
complete list) with the generic formula AMM’Q3 have been synthesized
experimentallyKoscielski and Ibers (2012); Prakash  _et al._ (2015); Ruseikina
_et al._ (2017, 2019a, 2019b); Sikerina and Andreev (2007); Maier  _et al._
(2016); Azarapin  _et al._ (2020); Strobel and Schleid (2006) which reveal
that these compounds possess rich chemistries and structure types like the
perovskites and Heusler compounds. Koscielski et al.Koscielski and Ibers
(2012) noted that these known AMM’Q3 compounds contain no Q-Q bonds and the
elements (A, M, M’, Q) balance the charge in their crystal structures with
their expected formal oxidation states, making them charge balanced. The
AMM’Q3 compounds are further classified into three categories depending on the
nominal oxidation states of the three cations A, M, and M’, namely:

  * •

Type-I (A1+M1+M’4+Q3)

  * •

Type-II (A2+M1+M’3+Q3)

  * •

Type-III (A1+M2+M’3+Q3),

where the oxidation states of the cations are indicated with the superscripts.
In all cases, we assume that the chalcogen atoms Q (S, Se, and Te) have a
nominal 2- charge. Examining the experimentally known AMM’Q3 compounds (see SI
Figure S1), we observe the following: (a) The A site in these compounds is
always occupied by alkali, alkaline earth, or post-transition metals with the
only exception of Eu (in a 2+ charge state) which occupies the A-site of some
of the Type-II compounds. (b) Whereas only transition metals occupy the M
site, the M’ site can be filled either by the transition metals, lanthanides,
or actinides. (c) No observed AMM’Q3 compound contains more than one alkali,
alkaline earth, or post-transition metals. As detailed later, we will use
these criteria in designing our HT workflow for generating the initial crystal
structures through prototype decoration. Although the crystal chemistries of
these compounds have been characterized somewhat in detailKoscielski and Ibers
(2012); Prakash  _et al._ (2015); Ruseikina  _et al._ (2017, 2019a, 2019b);
Sikerina and Andreev (2007); Maier  _et al._ (2016); Azarapin  _et al._
(2020); Strobel and Schleid (2006), their properties have remained largely
unexplored. Recently, it was shown experimentallyHao  _et al._ (2019) and
theoreticallyHao  _et al._ (2019); Pal  _et al._ (2019a, b, c) that many known
semiconducting compounds in this crystal family exhibit ultralow $\kappa_{l}$.
In addition, some compounds are shown to possess electronic bands favorable to
support high thermoelectric performanceHao  _et al._ (2019); Pal  _et al._
(2019b, c).

Figure 2: Elemental distributions of the experimentally known and
theoretically decorated AMM’Q3 compounds. The AMM’Q3 compounds are categorized
into three types depending on the nominal oxidation states of the cations:
Type-I (A1+M1+M’4+Q3), Type-II (A2+M1+M’3+Q3) and Type-III (A1+M2+M’3+Q3). The
color-coded periodic table indicates the group of elements (whose oxidation
states are written in blue color) that are used for prototype decorations in
our HT-DFT calculations. Elements that are denoted with the check marks
constitute the experimentally known AMM’Q3 compounds. Salmon, cyan and yellow
colors represent elements that occupy A, M/M’ and Q sites, respectively. The
grey color represents the elements that were not used for prototype
decoration. Eu is the only element that occupies both A and M’ site. All
experimental and decorated compounds are charge-balanced.

From the distribution of the elements forming the known AMM’Q3 compounds, we
see that all the cations (A, M, M’) are in their common oxidation states of
1+, 2+, 3+ or 4+ coming from the chemical groups (i.e., alkali, alkaline
earth, transition, post-transition metals, lanthanides, actinides) that span a
large part of the periodic table. Yet, although a large number (192) of these
compounds have been reported experimentally, this number is small compared to
the vast number of possible compounds that can be obtained based on charge-
balanced combinatorial substitutions of the elements in a prototype crystal
structure of AMM’Q3. Performing this combinatorial exercise experimentally
would require a massive amount of resources and time to discover new AMM’Q3
compounds. However, computational screening can be very helpful in narrowing
down the search space of the target compounds that would have higher chances
of synthesizability in the laboratoryCurtarolo  _et al._ (2013); Greeley  _et
al._ (2006); Gautier  _et al._ (2015); Zakutayev  _et al._ (2013); Anand  _et
al._ (2019); Hautier  _et al._ (2012); Sun  _et al._ (2016); Aykol  _et al._
(2018); Jain  _et al._ (2011); Ward and Wolverton (2017); Haastrup  _et al._
(2018). Here, we have performed HT-DFT calculations followed by accurate
ground-state phase stability analysis, and suggest (T = 0 K) thermodynamically
stable and metastable AMM’Q3 compounds for experimental synthesis and
exploration of their properties. Our calculations of the thermal transport
properties of some of the predicted stable compounds using the Peierls-
Boltzmann transport equation (PBTE) show that these compounds exhibit innate
low-$\kappa_{l}$ due to the presence of strong lattice anharmonicity or
rattler cations.

##  II Results

###  II.1 Structural prototypes

The experimentally known AMM’Q3 compounds crystallize in seven structure
typesKoscielski and Ibers (2012); Prakash  _et al._ (2015); Ruseikina  _et
al._ (2017, 2019a, 2019b); Sikerina and Andreev (2007); Maier  _et al._
(2016); Azarapin  _et al._ (2020); Strobel and Schleid (2006): KCuZrSe3 (space
group (SG): Cmcm, #63), Eu2CuS3 (SG: Pnma, #62), BaCuLaS3 (SG: Pnma, #62),
Ba2MnS3 (SG: Pnma, #62), NaCuTiS3 (SG: Pnma, #62), BaAgErS3 (SG: C2/m, #12)
and TlCuTiTe3 (SG: P21/m, #11). All these structure types are visualized in
their extended unit cells in Figure 1, where the conventional unit cell is
outlined with black lines. Among these, five structure types (KCuZrSe3,
Eu2CuS3, Ba2MnS3, NaCuTiS3, TlCuTiTe3) are layered where the rows of Am+
cations stack alternatively with the [MM’Q3]m- layers and interact through
electrostatic interactionsPal  _et al._ (2019c, a). The strength of the
interactions vary with the charges on the cations (i.e., m+) and induce
modifications in the structures as well as in their propertiesPal  _et al._
(2019b). When the interactions between the layers increase significantly, the
atoms from the neighboring [MM’Q3]m- layers interact, giving rise to the
three-dimensional channel structures (BaCuLaS3 and BaAgErS3). In BaCuLaS3
(Figure 1c) and BaAgErS3 (Figure 1f), two rows of the A-site cations occupy
the empty spaces inside the channels formed by the M, M’, and Q atoms. Figure
1h shows that 71% of the known AMM’Q3 compounds crystallize in the KCuZrSe3
structure followed by 16% of the compounds crystallizing in the Eu2CuS3
structure type. The rest of the known AMM’Q3 compounds (13 %) crystallize in
the other five structure types. It is worth noting that in Eu2CuS3, the Eu
atoms have mixed oxidation states (Eu2+CuEu3+S3) and sit in two different
sites in the crystal structure. Similarly, the Ba atoms in Ba2MnS3 also occupy
two different sites. We have used all these structure types in our HT-DFT
design and discovery of new AMM’Q3 compounds. We note that the KCuZrSe3 and
TlCuTiTe3 structure types have 12 atoms in their primitive unit cells and the
rest of the structure types have 24 atoms.

Figure 3: A brief workflow and summary of high-throughout DFT calculations.
(a) A schematic flow-chart of the HT-DFT calculations for the discovery of new
AMM’Q3 compounds. (b) The distribution of 628 predicted stable compounds in
the 7 structure types. (c) 570 out of 628 compounds possess finite band gaps.
Histogram plot for the DFT-calculated band gaps of these 570 semiconductors.

###  II.2 Materials design strategy

The discovery of new compounds through HT-DFT method often starts with the
decoration of prototype crystal structures with chemically similar elements
from the periodic table to generate their initial crystal structures. DFT
calculations are then performed on these newly decorated compounds followed by
rigorous thermodynamic phase stability analysis to screen for stable and
metastable compounds. Rather than generating the input crystal structures in a
brute-force manner by substituting every element of the periodic table at all
atomic sites in the prototype structures, in this work we restricted ourselves
by the following screening criteria which are derived by examining the
experimentally known AMM’Q3 compounds: (a) We substitute alkali, alkaline
earth, or post-transition metal elements at the A site. The M and M’ sites are
populated with the transition metals and lanthanides. Three chalcogens i.e.,
S, Se, and Te are substituted at the Q site. (b) We choose A-site cations with
nominal oxidation states of 1+ or 2+, and M/M’-site cations with nominal
oxidation states of 1+, 2+, 3+, or 4+. The elements that are chosen for
substitutions at the A, M, M’ and Q sites along with their oxidation states
are shown in Figure 2. (c) We only consider compound compositions that are
charge-balanced. (d) We exclude any radioactive elements during the prototype
decorations although some of the known AMM’Q3 compounds contain them. Adhering
to these pre-conditions helps us narrow down our search space of compound
exploration, reduces the computational cost, and most importantly it increases
the success rates of stable compound prediction through HT-DFT calculations
which will be evident later.

First, we generate the crystal structures of 4659 AMM’Q3 compositions (see
Table I) using the KCuZrSe3 structure as (i) it is the most prevalent
structure type in this family and (ii) all experimentally known AMM’Q3
compounds have low energies (within 50 meV/atom above the convex hull) in this
structure type. After performing DFT calculations for all these compositions
followed by T=0 K phase stability analysis, we kept only those compounds (
$\sim$ 1700) that have an energy within 50 meV/atom of the ground state convex
hull and discarded the rest from our search space. The DFT relaxed structures
of all these 1700 compounds retain the KCuZrSe3 structure type. In the next
step, we take these 1700 compositions and regenerate their crystal structures
in each of six additional structure types to perform DFT calculations of 6
$\times$ 1,700 = 10,200 AMM’Q3 compounds. Next, we perform T = 0K
thermodynamic phase stability analysis of those 1700 compositions considering
all seven structure types and their competing phases that are available in the
OQMD. In the final step, we obtain 628 thermodynamically stable and 852
metastable hitherto unknown AMM’Q3 compounds after performing a total number
of 4,659 + 10,200 = 14,859 DFT calculations. The stable 628 compounds include
69 Type-I, 231 Type-II, and 328 Type-III compounds, and among them, a total
number of 570 compounds possess finite band gaps. A schematic of the HT-DFT
flowchart is shown in Figure 3a, and a summary is given in Table I.

Figure 4: Elemental distributions of the predicted stable AMM’Q3 compounds.
Distributions of the cations (A, M, M’) forming the 628 predicted stable
AMM’Q3 compounds in the three categories: Type-I, Type-II, and Type-III. The
bar corresponding to an element represents the number of stable compounds
containing it. Total number of sulfides, selenides and tellurides are
mentioned in each panel. These compounds exclude the experimentally known 192
AMM’Q3 compounds.

###  II.3 Phase stability analysis

We now present a detailed analysis of the T= 0 K ground state phase stability
of all (known and predicted) AMM’Q3 compounds. We begin our assessment with
the phase stability analysis of the experimentally known AMM’Q3. Out of 192
known compounds, we find only 119 compounds (Type-I: 39, Type-II: 30, and
Type-III: 50) in OQMD before we performed any new calculations from this work.
We designate these 119 compounds as Set-I. As the initial crystal structures
of the experimentally known compounds in OQMD mostly come from the ICSD, the
DFT calculations of Set-I compounds in OQMD were performed based on their
experimental crystal structures taken from the ICSD. We designate the rest of
the known 192-119 = 73 compounds (Type-I: 1, Type-II: 60, and Type-III: 12) as
Set-II, which were not present in OQMD due to (a) the absence of their
structures in the ICSD and (b) no previous HT-DFT calculations were performed
based on the prototype decorations in this AMM’Q3 family. However, our HT-DFT
calculations of all the decorated AMM’Q3 compositions include the
experimentally known AMM’Q3 compounds in Set-II. Hence, we will first analyze
the phase stability of the Set-I compounds to see if our DFT calculations are
able to correctly capture the energetics of the known compounds in Set-I and
then we will utilize the phase stability data of Set-II compounds to validate
the reliability of our approach for the discovery of new stable AMM’Q3
compounds based on prototype decoration through HT-DFT calculations.

####  II.3.1 Set-I compounds

As detailed in the Methods section and as well as in other referencesCurtarolo
_et al._ (2013); Saal  _et al._ (2013); Kirklin  _et al._ (2015); Jain  _et
al._ (2013); Curtarolo  _et al._ (2012); Jain  _et al._ (2016, 2011); Emery
_et al._ (2016) the hull distance (hd) is a metric of the thermodynamic
stability of a compound. If the formation energy of the compound breaks the
convex hull, then it is considered to be thermodynamically stable with hd = 0,
indicating the likelihood of its synthesizability. On the other hand,
compounds with a small positive hull distance (typically within a few tens of
meV/atom) are called metastable and may also be in some cases experimentally
synthesizedSun  _et al._ (2016); Aykol  _et al._ (2018). According to this
criterion, all experimentally known AMM’Q3 compounds should possess zero or
small positive hull distances. Our analysis reveals that in Set-I, all (39)
Type-I compounds and all but one (29) Type-II compounds have hd = 0, which is
in line with our expectations. The one Type-II compound that has a small
positive hull distance is Eu2+CuEu3+S3 (hd = 37 meV/atom). Also, all but three
(47) Type-III compounds of Set-I have hd = 0. These three Type-III compounds
are CsCoYbS3 (hd =192 meV/atom), CsCoYbSe3 (hd =151 meV/atom), CsZnYbSe3
(hd=76 meV/atom). So, 115 of 119 experimentally synthesized compounds in Set-I
are thermodynamically stable in the OQMD. Hence, stability is an excellent
metric for the synthesizability of the predicted compounds. From this
analysis, we also find that one of the Type-II compounds in Set-I, BaAgErS3,
possesses a small positive hd of 12 meV/atom when its calculation is performed
on the experimental crystal structure having C2/m SG (#12), indicating that it
is metastable at 0 K in this structure. Interestingly, our HT-DFT calculations
reveal that BaAgErS3 is stable T = 0 K in the KCuZrSe3 structure type (see SI
Figure S2).

####  II.3.2 Set-II compounds

Since the 73 compounds in Set-II did not exist in OQMD before our HT-DFT
calculations, we generate their crystal structures through prototype
decorations as mentioned before to perform DFT calculations and T = 0 K phase
stability analysis. As these compounds have already been synthesized
experimentally, our DFT calculations and phase stability analysis provide a
key test of our methodology. It also gives us an opportunity to examine how
reliably our calculations can predict hitherto unknown stable AMM’Q3
compounds. After performing the thermodynamic stability analysis, we found
that among the 73 compounds in Set-II (type-I: 1, type-II: 60, type-III: 12),
64 compounds have hd = 0, and only 4 compounds (SrCuYbS3: hd = 84 meV/atom,
BaCuYbTe3: hd=62 meV/atom, EuCuYbS3: hd=72 meV/atom, PbCuYbS3: hd =104
meV/atom) of Type-II and 5 compounds (CsMnYbSe3: hd=8 meV/atom, CsZnYbS3 :
hd=163 meV/atom, CsZnYbTe3: hd=71 meV/atom, RbZnYbSe3: hd=102 meV/atom,
RbZnYbTe3: hd=82 meV/atom) of Type-III have positive hull distances.

|  Type-I |  Type-II |  Type-III |  Total
---|---|---|---|---
A |  Li, Na, K, Rb, Cs, In, Tl (7) |  Mg, Ca, Sr, Ba, Ge, Sn, Pb, Eu (8) |  Li, Na, K, Rb, Cs, In, Tl (7) |  -
M |  Cu, Ag, Au (3) |  Cu, Ag, Au (3) |  Mn, Fe, Co, Ni, Zn, Cd, Hg, Eu (8) |  -
M’ |  Ti, Mn, Zr, Mo, Ru, Pd, Hf, W, Re, Os, Ir, Pt, Ce (13) |  Sc, Y, La, Ce, Pr, Nd, Sm, Eu Gd, Tb, Dy, Ho, Er, Tm, Yb, Lu (16) |  Sc, Y, La, Ce, Pr, Nd, Sm, Eu Gd, Tb, Dy, Ho, Er, Tm, Yb, Lu (16) |  -
Q |  S, Se, Te (3) |  S, Se, Te (3) |  S, Se, Te (3) |  -
Total number of initial compositions |  7 $\times$ 3 $\times$ 13 $\times$ 3 = 819 |  8 $\times$ 3 $\times$ 16 $\times$ 3 = 1152 |  7 $\times$ 8 $\times$ 16 $\times$ 3 = 2688 |  4659
Thermodynamically stable compounds |  69 |  231 |  328 |  628
Thermodynamically metastable compounds |  59 |  282 |  511 |  852
Table 1: A table summarizing the design and discovery of new materials through
prototype decoration. Elements chosen for decorating the prototype AMM’Q3
structures and results of HT-DFT calculations are shown here. The values
within the parentheses indicate the total number of elements used in each of
three categories of compounds for substitutions at the A, M, M’ and Q sites.
All generated initial compositions satisfy the charge-neutrality criteria
based on their nominal oxidation states. The numbers of predicted stable (hd =
0) and low-energy metastable (0 $<$ hd $\leq$ 50 meV/atom) compounds in the
table exclude the experimentally known 192 AMM’Q3 compounds.

Thus, our HT-DFT calculations based on the decorated structures successfully
capture the stability of all but 13 (out of 192) of the experimentally known
AMM’Q3 compounds. The 13 compounds which are experimentally observed but with
hd $>$ 0 all contain Eu and Yb, and the calculated energetics of these
compounds may originate from the choice of incorrect pseudopotentials as
discussed later. These results and analysis give us strong confidence in
designing and discovering new AMM’Q3 compounds using the HT-DFT calculations
and thermodynamic T=0 K stability as a metric for synthesizable compounds. We
have provided the phase stability data of all 192 experimentally known AMM’Q3
compounds in the SI.

###  II.4 Novel AMM’Q3 compounds

After performing the T= 0 K stability analysis of all newly decorated AMM’Q3
compounds for which HT-DFT calculations are performed in all 7 structure
types, we discover a large number of 628 (type-I: 69, type-II: 231, type-
III:328) thermodynamically stable compounds that exclude the experimentally
known 192 compounds. To put this number into perspective, the OQMD (containing
more than 900,000 entries as of September 2020) has stable (hd = 0) 1161 full-
Heuslers (SG: Fm$\bar{3}$m, #225), 618 half-Heuslers (SG: F$\bar{4}$3m, #216),
353 cubic (SG: Pm$\bar{3}$m, #221) perovskites, 242 orthorhombic (SG: Pnma,
#62) perovskites. Similar to the experimentally known AMM’Q3 compounds, the
KCuZrSe3 and Eu2CuS3 structure types are the most common, which constitute 41%
and 25% of the predicted stable compounds, respectively (Figure 3b). Among the
other structure types, NaCuTiSe3 and TlCuTiTe3 are quite common, which
constitute 19% and 9% of the predicted stable compounds. The rest 6% compounds
have the BaCuLaS3 and BaAgErS3 structure types. However, we found no new
stable compound in the Ba2MnS3 structure type. Our analysis shows that 570 out
of 628 compounds possess finite bands that range from 0.2 eV to 5.34 eV, among
which most of the compounds have band gaps within 1.5 eV (Figure 3c). This is
not surprising since all the decorated compositions are charge-balanced. In
addition, we found 852 potentially synthesizable metastable compounds (Type-I:
59, Type-II: 282, Type-III: 511) with small positive hull distances (i.e., 0
$<$ hd $\leq$ 50 meV/atom). A summary of the HT-DFT calculations is shown in
Table-I, and the lists of all predicted stable and metastable compounds are
given in the SI.

We further examine the predicted stable compounds within each type in terms of
their chemistries (sulfides, selenides, and tellurides) and the cations (A, M,
M’). The results are displayed as bar charts in Figures 4, where the bar
corresponding to an element represents the number of stable compounds that
contain it. We see that: (a) the elements that form the stable compounds
constitute almost the entire periodic table. (b) There are more Cu-containing
compounds in Type-I and Type-II categories compared to Ag or Au. In Type-III,
there are more Mn compounds than other M elements with a 2+ oxidation state.
(c) The number of stable compounds increases from Type-I (69) to Type-II (231)
to Type-III (328). This is not surprising given that the number of elements
that can occupy the M and M’ sites (satisfying the charge-neutrality criteria)
increase from Type-I to Type-III compounds (Table-I). A general trend that is
noticeable across three types is that as we go from sulfides to selenides to
tellurides, the number of compounds with smaller A cations decrease whereas
the number increases with larger A cations. For example, there are two, two,
and no Li-containing compounds in sulfides, selenides, and tellurides of
Type-I compounds, respectively. Similarly, the number of compounds that have
Sr increase from sulfides (7) to selenides (20) to tellurides (34) in Type-II.

###  II.5 Lattice thermal transport properties

We now focus on exploring the lattice thermal transport properties of the
predicted stable AMM’Q3 compounds. An accurate estimation of $\kappa_{l}$ of a
compound within a first-principles DFT framework is computationally very
expensiveLi  _et al._ (2014). Hence, the calculations of $\kappa_{l}$ for all
the predicted stable compounds would require a massive amount of computational
resources. However, to demonstrate the thermal transport properties of our
newly predicted stable compounds, we randomly select a handful of compounds
with some criteria. We first screen for non-magnetic and semiconducting
compounds, where the lattice contribution dominates the total thermal
conductivity. Next, we search for those compounds which have the KCuZrSe3
structure type as it possesses the highest crystal symmetry (SG: Cmcm, #63)
and the smallest unit cell (12 atoms). Finally, we randomly select 10
compounds for $\kappa_{l}$ calculations. The selected compounds, which include
sulfides, selenides and tellurides, are: CsCuZrS3, BaCuScSe3, BaCuScTe3,
BaCuTbSe3, BaAgGdSe3, CsZnYS3, CsZnGdS3, CsZnScSe3, CsZnScTe3, CsCdYTe3. We
note that this list also includes Type-I (the first in the list), Type-II (the
next 4) as well as Type-III (the last 5) compounds. The electronic structures
and phonon dispersions of these compounds are given in the SI.

Figure 5: Lattice thermal transport properties of the AMM’Q3 compounds. The
calculated temperature-dependent lattice thermal conductivity ($\kappa_{l}$)
of 10 predicted stable AMM’Q3 compound which are non-magnetic and
semiconducting. $\kappa_{l}^{\bot}$ (blue disk) and $\kappa_{l}^{\|}$ (orange
pentagon) are the two components of $\kappa_{l}$ that are perpendicular and
parallel to the stacking directions in the crystal structure, respectively.

We calculate the $\kappa_{l}$ of these 10 compounds using the PBTE (see
Methods section) and present the results in Figure 5. We see that all these
compounds exhibit very low $\kappa_{l}$ where the in-plane
$(\kappa_{l}^{\bot})$ and the cross-plane $(\kappa_{l}^{\|})$ components are
lower than 3 Wm-1K-1 and 1.2 Wm-1K-1, respectively, for T $\geq$ 300 K. Here,
$\kappa_{l}^{\bot}$ is perpendicular to the stacking direction of the layers
in the crystal structure of the AMM’Q3 compounds and $\kappa_{l}^{\|}$ is
parallel to it. As a reference, we compare our calculated $\kappa_{l}$ with
that of a prototypical thermoelectric material SnSe, which was experimentally
shown to possesses low-$\kappa_{l}$ that leads to a very high thermoelectric
figure-of-meritZhao  _et al._ (2014). The measuredWu  _et al._ (2017)
$\kappa_{l}^{\bot}$ and $\kappa_{l}^{\|}$ for single-crystalline
stoichiometric samples of SnSe are 1.9 Wm-1K-1 and 0.9 Wm-1K-1, respectively,
at 300 K, which become much lower in the off-stoichiometric polycrystalline
samplesZhao  _et al._ (2014). Further examination of the results in Figure 5
reveals that in terms of the anisotropy of the $\kappa_{l}^{\bot}$ and
$\kappa_{l}^{\|}$ components, Type-I and Type-III compounds are quite similar,
but Type-II compounds are different from the rest i.e.,
$\kappa_{l}^{\bot}/\kappa_{l}^{\|}$ (Type-I/III) $>$
$\kappa_{l}^{\bot}/\kappa_{l}^{\|}$ (Type-II). The difference in the
anisotropy of the properties arises from the fact that in Type-II compound the
electrostatic attractions between A2+ and [MM’Q3]2- layers are stronger than
that of between A1+ and [MM’Q3]1- layers in Type-I/III compounds. The stronger
interlayer interactions give rise to a shorter interlayer distance in Type-II,
which make $\kappa_{l}^{\bot}$ and $\kappa_{l}^{\|}$ less anisotropic.

To gain a deeper understanding, we examine the lattice dynamics and thermal
transport properties of BaCuScTe3 (Type-II) and CsCdYTe3 (Type-III) in detail.
Our analysis reveals that the underlying physical principles governing the
low-$\kappa_{l}$ in Type-II compounds are different from the Type-I/III
compounds. We start our analysis with the phonon dispersion of BaCuScTe3
(Figure 6a) which shows the presence of very low-frequency acoustic ($<$ 45
cm-1), and optical ($\sim$ 15 cm-1 along X-S-R directions) phonon modes, which
give rise to low phonon group velocities. The phonon dispersion and the phonon
density of states (Figure 6b) of BaScCuTe3 show that a strong hybridization
exists between the phonon branches up to 100 cm-1, where Ba and Te atoms have
large contributions. Soft acoustic phonon branches give rise to very low sound
velocities and a strong hybridization between the phonons at low energies
enhances the phonon scattering phase space. Both these factors help suitably
to give rise to a very low-$\kappa_{l}$ ($\kappa_{l}^{\bot}$ = 1.7 Wm-1K-1 and
$\kappa_{l}^{\|}$ = 0.76 Wm-1K-1 at 300 K). On the other hand, the phonon
dispersion of CsCdYTe3 (Figure 6d) features nearly dispersion-less optical
phonon branches along the X-S-R-Z directions in the Brillouin zone at low
energies, which are the characteristics of rattler vibrations in the crystal
structure. In addition, it also has soft acoustic phonon branches ($<$ 35
cm-1). The calculated $\kappa_{l}$ of CsCdYTe3 becomes ultralow with
$\kappa_{l}^{\bot}$ and $\kappa_{l}^{\|}$ being 0.82 Wm-1K-1 0.25 Wm-1K-1,
respectively at 300 K.

Rattler phonon modes are highly localized, which strongly inhibit the
transport of phonons, giving rise to ultralow-$\kappa_{l}$ in many crystalline
solids such as TlInTe2Jana  _et al._ (2017), CsAg5Te3Lin  _et al._ (2016),
etc. It was shown that the filler atoms in clathratesTadano  _et al._ (2015)
and skutteruditesLi and Mingo (2015) act as ideal rattlers, which give rise to
dispersion-less phonon branches where the phonon frequencies remain highly
localized having very small participation ratio (PR) values $\sim$ 0.2 (see
Methods section). The PRs of the phonon modes of BaCuScTe3 and CsCdYTe3 are
color-coded in Figures 6a and 6d, respectively. We see that while most of the
low-energy phonon modes ($<$ 100 cm-1) of BaCuScTe3 have PR values close to 1,
signifying the absence of phonon localization, the low-energy dispersion-less
phonon branches of CsCdYTe3 have small PRs ( $<$ 0.2), indicating their highly
localized nature of the phonon modes. The phonon density of states (Figure 6e)
also reveals that these localized phonons primarily arise from the Cs atoms
(confined in 25-55 cm-1) that act as rattlers. Thus, our analysis shows that
the ultralow-$\kappa_{l}$ in CsCdYTe3 is primarily caused by the localized
vibrations of the rattling phonons.

Figure 6: Harmonic and anharmonic lattice dynamical properties of the AMM’Q3
compounds. Harmonic phonon dispersions, atom-resolved phonon density of
states, anharmaonic frequency-cumulative $\kappa_{l}$, and their derivatives
with respect to the frequency for BaCuScTe3 (a-c) and CsCdYTe3 (d-f). The
phonon dispersions are color coded with the participation ratios of the phonon
modes (indicated in the color bar). $\kappa_{l}^{\bot}$ (solid blue line) and
$\kappa_{l}^{\|}$ (solid red line) are perpendicular and parallel to the
stacking directions in the crystal structure of the AMM’Q3 compound,
respectively. ${\kappa^{\prime}_{l}}^{\bot}$ (dotted blue line) and
${\kappa^{\prime}_{l}}^{\|}$ (dotted red line) are the first-order derivatives
of $\kappa_{l}^{\bot}$ and $\kappa_{l}^{\|}$, respectively, with respect to
the frequency, which are in arbitrary units in the figures.

We further investigate the origin of such poor thermal transport properties in
BaCuScTe3 in terms of a more fundamental material quantity, the lattice
anharmonicity. Strong lattice anharmonicity is one of the important factors
that induces very low-$\kappa_{l}$ in compounds like SnSeLi  _et al._ (2015),
NaSbTe2Nielsen  _et al._ (2013), etc. To estimate the strength of the
intrinsic anharmonicity, we calculate the macroscopic average Gruneisen
parameter, $\gamma=\sum_{i}{\gamma_{i}}C_{v,i}/\sum_{i}C_{v,i}$, where
$\gamma_{i}$ and $C_{v,i}$ are the Gruneisen parameter and specific heat
capacity at constant volume for the $i$-th phonon mode. The calculated
$\gamma$ of BaCuScTe3 (1.5) is larger than that of CsCdYTe3 (1.2), signifying
the presence of stronger anharmonicity in the former. These $\gamma$ values
are comparable to that of NaSbSe2 (1.7), NaSbTe2 (1.6), NaBiTe2 (1.5) etc.
compounds which are experimentallyNielsen  _et al._ (2013) shown to possess
ultralow-$\kappa_{l}$. Thus, we see that while the low-$\kappa_{l}$ in Type-II
compounds is caused by the low-sound velocities and a stronger lattice
anharmonicity, the presence of rattler cations in Type-III as well Type-I
compounds (see SI Figure S4) are primarily responsible for inducing
low-$\kappa_{l}$ in them.

To examine which phonons are primarily responsible for the conduction of heat
in BaCuScTe3 and CsCdYTe3, we plot the cumulative-$\kappa_{l}$ and their
first-order derivatives ($\kappa^{\prime}_{l}$) with respect to the frequency
at T = 300 K. We see from Figure 6c that while the acoustic and low-energy
optical phonons up to 90 cm-1 mainly contribute to $\kappa_{l}^{\bot}$,
$\kappa_{l}^{\|}$ is primarily contributed by the phonons up to 45 cm-1 in
BaCuScTe3. On the other hand, Figure 6f shows that while the acoustic as well
as the optical phonons up to 100 cm-1 have large contributions towards
$\kappa_{l}^{\bot}$, only the acoustic phonons (up to 35 cm-1) primarily carry
heat for $\kappa_{l}^{\|}$ in CsCdYTe3. We also notice that the anisotropy
$(\kappa_{l}^{\bot}/\kappa_{l}^{\|}$ = 3.3 at T = 300 K) in CsCdYTe3 is much
larger that of BaCuScTe3 ($\kappa_{l}^{\bot}/\kappa_{l}^{\|}$ = 2.2 at T = 300
K). The origin of this anisotropy can be attributed to the contrasting
interlayer and intralayer interactions in BaCuScTe3 and CsCdYTe3. For example,
the analysis of the interatomic force-constants (IFCs) reveals that interlayer
interactions in CsCdYTe3 are much weaker (IFC(Cs-Te) = -0.333 eV/$\AA^{2}$)
compared to BaCuScTe3 (IFC(Ba-Te) =-1.204 eV/$\AA^{2}$), which makes the
transport of optical phonons (above 25 cm-1) along the stacking direction
(i.e., $\kappa_{l}^{\|}$) of CsCdYTe3 less effective. On the other hand, the
intralayer interactions in CsCdYTe3 are much stronger (IFC(Cd-Te)= -3.923
eV/$\AA^{2}$, IFC(Y-Te): -2.216 eV/$\AA^{2}$) than those of BaCuScTe3 (IFC(Cu-
Te)= -2.396 eV/$\AA^{2}$, IFC(Sc-Te)= -1.890 eV/$\AA^{2}$). As a result, while
phonons up to 90 cm-1 mainly carry the heat for $\kappa_{l}^{\bot}$ in
BaCuScTe3, in CsCdYTe3 they are carried by phonons with frequencies up to 100
cm-1 very effectively.

##  III Discussion

We notice that in Set-I and Set-II, all the known compounds that have positive
hull distances have either Yb3+ or Eu3+ cations in them. Given the fact OQMD
used Yb_2 and Eu_2 PPs which are intended for compounds having Yb2+ and Eu2+
cations, those energetic results are somewhat suspect. However, we note that
our HT-DFT calculations predict the stability of other rare-earth elements
containing known compounds in this family correctly. Hence, the newly
predicted compounds containing Eu3+ and Yb3+ should be taken with caution. We
note that none of our predicted stable compounds contain these cations.
However, 11 and 7 of the predicted metastable compounds contain Eu3+ and Yb3+
cations, respectively.

Concerning the experimental validation of our prediction, we note that the
experimental synthesis of a large number of compounds is a daunting task.
While automated synthesis of a large number of compounds is now possible
through high-throughput experimental facilitiesGreen  _et al._ (2017), here,
we suggest only four compounds (Table II) that can be immediately picked up
for experimental verification of our prediction. These compounds contain no
toxic element and their calculated lattice thermal conductivity is very low.
In Table II, we provide their DFT calculated band gaps, $\gamma$, and average
$\kappa_{l}$ calculated at 300 K. These calculated quantities along with
$\kappa_{l}$ can be compared with the experimentally measured values of these
materials. Finally, we note that $\kappa_{l}$ for each compound has been
calculated using only three-phonon scattering processes. The inclusion of
additional four-phonon scattering ratesXia  _et al._ (2020a, b) and grain-
boundaryPal  _et al._ (2020) limited phonon scattering mechanisms could
further lower the calculated $\kappa_{l}$ in this family of compounds. Also,
as the electronic structure of the compounds features nearly flat bands and
multiple peaks near the valence/conduction band extrema (Figure S3, SI), some
of these compounds are expected to exhibit good potential for thermoelectric
applications as well.

Compound |  Band gap (eV) |  Average $\gamma$ at 300 K |  Average $\kappa_{l}$ (Wm-1K-1) at 300 K
---|---|---|---
BaCuScTe3 |  0.39 |  1.51 |  1.23
BaAgGdSe3 |  1.45 |  1.60 |  1.24
CsZnScTe3 |  1.20 |  1.30 |  0.79
CsCdYTe3 |  1.74 |  1.30 |  0.54
Table 2: Suggested compounds for experimental synthesis and characterization.
A list of 4 predicted stable compounds are suggested for the experimental
synthesis and measurement of their properties. The DFT-calculated band gap,
macroscopic average Gruneisen parameter ($\gamma$), the average $\kappa_{l}$
are provided to compare against experiments.

In summary, we use HT-DFT calculations to discover a large number of 628
thermodynamically stable quaternary chalcogenides (AMM’Q3). As all
compositions in this family are charge-balanced, our analysis shows that 570
of 628 compounds possess finite band gaps which vary between 0.2 and 5.34 eV.
Our calculations of the thermal transport properties show that AMM’Q3
compounds exhibit intrinsically very low $\kappa_{l}$, and the anisotropy in
$\kappa_{l}$ is much smaller in Type-II compounds compared to Type-I/III
compounds. Our analysis further reveals that low-$\kappa_{l}$ in this family
originates either due to the presence of rattling cations (in Type-I/III
compounds) or stronger lattice anharmonicity (Type-II compound). While the
rattler cations give rise to localized phonon modes that inhibit the
propagation of phonons, a stronger lattice anharmonicity enhances the phonon
scattering phase space, leading to a low-$\kappa_{l}$. In addition, there
exists a strong coupling between the acoustic and low-energy optical phonon
modes in Type-II compounds which increases the phonon scattering rates of the
heat-carrying phonons. Our work is thus interesting not just from the
perspective of new materials discovery but also for finding the presence of
low $\kappa_{l}$ in them, which hold promises for further research and
possible applications in energy materials and related devices.

##  IV Methods

###  IV.1 DFT calculations

We performed DFT calculations using the Vienna Ab-initio Simulation Package
(VASP)Kresse and Furthmüller (1996) using the projector-augmented wave
(PAW)Kresse and Joubert (1999) potentials with the Perdew-Burke-Ernzerhof
(PBE)Perdew  _et al._ (1996) generalized gradient approximation (GGA) to the
exchange-correlation functional. The atomic positions and other cell degrees
of freedom of the compounds were fully relaxed and spin-polarized calculations
were performed for compounds that contain partially filled d or f-shells
elements with a ferromagnetic arrangement of spins in accordance with the
high-throughput framework as laid out in the qmpy suite of toolsSaal  _et al._
(2013); Kirklin  _et al._ (2015). For more details on the calculation
parameters, we refer to Ref.Saal  _et al._ (2013); Kirklin  _et al._ (2015).
T=0 K phase stability analysis often serves as an excellent indicator for the
possibility of synthesizability of a predicted compound in the
laboratoryZakutayev  _et al._ (2013); Gautier  _et al._ (2015); Anand  _et
al._ (2019). To assess the T = 0 K thermodynamic stability of the compounds,
we calculate their formation energies ($\Delta$Hf) utilizing the DFT total
energy (ground state) of each compound using the formula, $\Delta$Hf($\sigma$)
= E($\sigma$) - $\sum_{i}$ n${}_{i}\mu_{i}$, where E is the DFT total energy
(at 0 K) of an AMM’Q3 compound in a crystal structure denoted by $\sigma$,
$\mu_{i}$ is the chemical potential of element $i$ with its fraction ni in
that compound. For each composition, we used a number of prototype crystal
structures, $\sigma$, based on structural prototypes of known AMM’Q3
compounds. To determine the thermodynamic stability of a compound, we need to
compare its formation energy against all its competing phases, not only
against other compounds at the same composition. To this end, we generated the
quaternary phase diagram (i.e., the T = 0 K convex hull of the A-M-M’-Q phase
space) for every AMM’Q3 compound considering all elemental, binary, ternary,
and quaternary phases present in the OQMD, which (as of September 2020)
corresponded to nearly 900,000 entries of DFT-calculated energies. The
calculated convex hull distance (hd, explained in the next section) then
serves as a metric to determine whether a compound is stable (i.e., hd = 0),
metastable (small positive hd), or unstable (large positive hd).

Figure 7: Visualization of quaternary and ternary phase diagrams. (a) The four
dimensional phase diagram (T= 0 K) of the Li-Cu-Zr-S quaternary system, which
is the isometric shot of a Gibbs’ tetrahedron. (b) One of the faces of the
tetrahedron represents the three-dimensional phase diagram of the Cu-Zr-S
ternary system which is presented as the Gibbs’ triangle. Each cyan node in
(a) and (b) represents a stable compound. For clarity, we do not show any
metastable/unstable compounds and mark only few stable compounds in this
figure. LiCuZrS3 is one of the predicted stable compounds in this work, which
and its competing phases are denoted with dashed arrows in (a).

###  IV.2 Convex hull construction

To construct the convex hull of an AMM’Q3 compound, it is necessary to
identify the set of phases (elemental, binary, ternary as well as quaternary)
in the four-dimensional composition space of A-M-M’-Q that have the lowest
formation energies at their compositions. In OQMD, we present the convex hull
of a quaternary compound through a four-dimensional phase diagram which is
represented as an isometric shot of the Gibbs’ tetrahedron (Figure 7a). Each
face of the tetrahedron represents a three-dimensional phase diagram of a
ternary composition that is represented as the Gibbs’ triangle (Figure 7b).
The vertices in Figures 7a & 7b represent the elements constituting the
quaternary and ternary phase space, respectively, and the edges connecting any
two vertices represent the binary composition axis between those two elements.
Any node within Figures 7a & 7b represents a stable compound, which is denoted
by a cyan disk. The metastable and unstable compounds are not shown in these
figures for clarity which fall off the nodes. The stability of an AMM’Q3
compound is given by the difference (which is defined as the hull distance, hd
= $\Delta$Hf \- $\Delta$He) between the calculated formation energy
($\Delta$Hf) of an AMM’Q3 compound under consideration and its hull energy
($\Delta$He) at that composition. The hull energy is defined as the energy at
the convex hull at that AMM’Q3 composition. By definition, the hd of a stable
compound is zero, whereas for metastable and unstable compounds they are real
positive numbers. In keeping with the heuristic conventions used in
literatureZakutayev  _et al._ (2013); Cerqueira  _et al._ (2015); Wu  _et al._
(2013) we term those AMM’Q3 compounds to be metastable whose hull distances
lie within 50 meV/atom above the convex hull (i.e., 0 $<$ hd $\leq$ 50
meV/atom). The metastable compounds are also potentially synthesizable in the
laboratory Sun  _et al._ (2016); Aykol  _et al._ (2018).

###  IV.3 Phonon participation ratio

Phonon dispersions have been calculated using 2 $\times$ 2 $\times$ 1
supercell of the primitive unit cell using PhonopyTogo and Tanaka (2015). The
high symmetry paths in the Brillouin zones were adopted following the
conventions used by Setyawan et al.Setyawan and Curtarolo (2010). To examine
the extent of localization of the phonon modes, we calculate their phonon
participation ratio using the formulaTadano  _et al._ (2015); Pailhès  _et
al._ (2014):
$P(\omega_{q})=(\sum_{i}^{N}\frac{{{|e}_{i}(\omega_{q})|}^{2}}{M_{i}})^{2}/\sum_{i}^{N}\frac{{{|e}_{i}(\omega_{q})|}^{4}}{M_{i}^{2}}$,
where $e_{i}(\omega_{q})$ is the eigenvector of the phonon mode at q with
frequency $\omega$, $M_{i}$ is the mass of the $i$-th atom in the unit cell
containing a total number of $N$ atoms. The value of $P(\omega_{q})$ ranges
between 0 and 1. In an ordered crystal, when $P(\omega_{q})$ becomes close to
1, it indicates that the phonon mode is propagative where all the atoms in the
unit cell participate. On the other hand, very low values of PR ($\sim$
0.2)Pailhès  _et al._ (2014); Beltukov  _et al._ (2013) indicate the strong
localization of the phonon modes (e.g., rattling phonons) where only a few
atoms in the unit cell participate in the vibrations. Examples of rattler
atoms containing compounds include filled clathrates (Ba8Si46 and
Ba8Ga16Ge30)Tadano  _et al._ (2015); Pailhès  _et al._ (2014), where the
filler atoms act as ideal rattlers that induce ultralow-$\kappa_{l}$ in them.

###  IV.4 Thermal conductivity calculations

We calculate the $\kappa_{l}$ utilizing the phonon lifetimes obtained from the
third-order interatomic force constants (IFCs)Chaput  _et al._ (2011); Togo
_et al._ (2015); Zhou  _et al._ (2014), which was shown to reproduce
$\kappa_{l}$ within 5 % of the experimentally measured $\kappa_{l}$ in this
AMM’Q3 family of compoundsHao  _et al._ (2019); Pal  _et al._ (2019a). We
constructed the third-order IFCs of each compound based on DFT calculations of
displaced supercell configurations by limiting the cut-off distance (rc) up to
the third nearest neighbor. We used 2 $\times$ 2 $\times$ 1 supercell
(containing 48 atoms) of the primitive unit cell (with 12 atoms) using
thirdorder.pyLi  _et al._ (2014) utility for the calculation of IFCs. Using
the second and third-order IFCs in the ShengBTE codeLi  _et al._ (2014), we
calculate the temperature-dependent phonon scattering rates and $\kappa_{l}$
utilizing a full iterative solution to the PBTE for phonons using a 12
$\times$ 12 $\times$ 12 q-point mesh. The calculated $\kappa_{l}$ generally
depends on the rc which accounts for the maximum range of interaction in the
third-order IFCsLi  _et al._ (2015). It was shown that good convergence of
$\kappa_{l}$ was obtained by limiting rc to the third nearest neighbor within
the crystal structure in this family of compoundsPal  _et al._ (2019a).

##  V Data availability

The data that supports the findings of the work are in the manuscript and
Supplementary Information. The structures and energetics of the predicted
compounds would be made available through the Open Quantum Materials database
(OQMD) in a future release. Additional data will be available upon reasonable
request.

##  VI Code availability

Open-source codes are used throughout this work.

##  VII Materials & Correspondence

Request for materials and correspondence should be addressed to K.P. or C.W.

##  VIII Acknowledgements

###### Acknowledgements.

K.P. and C.W. acknowledge support from the U.S. Department of Energy under
Contract No. DE-SC0014520 (thermal conductivity calculations) and the Center
for Hierarchical Materials Design (CHiMaD) and from the U.S. Department of
Commerce, National Institute of Standards and Technology under Award No.
70NANB14H012 (HT-DFT calculations). J.S. and J.H. acknowledges support from
the National Science Foundation through the MRSEC program (NSF-DMR 1720139) at
the Materials Research Center (phase stability). Y.X. acknowledges support
from Toyota Research Institute (TRI) through the Accelerated Materials Design
and Discovery program (lattice dynamics). Y.L. and M.G.K. were supported in
part by the National Science Foundation Grant DMR-2003476. K.P. sincerely
thanks Sean Griesemer for useful discussion on the abundance of various
crystallographic prototypes in the OQMD. We acknowledge the computing
resources provided by (a) the National Energy Research Scientific Computing
Center (NERSC), a U.S. Department of Energy Office of Science User Facility
operated under Contract No. DE-AC02-05CH11231, (b) Quest high-performance
computing facility at Northwestern University which is jointly supported by
the Office of the Provost, the Office for Research, and Northwestern
University Information Technology, and (c) the Extreme Science and Engineering
Discovery Environment (National Science Foundation Contract ACI-1548562).

##  IX Author contributions

K.P. conceived and designed the project. K.P. performed calculations and
analysis with help and suggestions from Y.X., J.S., J.H., Y.L, M.G.K, and C.W.
C.W. supervised the whole project. All authors discussed the results, provided
comments, and contributed to writing the manuscript.

##  X Competing interests

The authors declare no competing financial or non-financial interests.

## References

  * Curtarolo  _et al._ (2013) S. Curtarolo, G. L. Hart, M. B. Nardelli, N. Mingo, S. Sanvito, and O. Levy, The high-throughput highway to computational materials design, Nature materials 12, 191 (2013).
  * Mounet  _et al._ (2018) N. Mounet, M. Gibertini, P. Schwaller, D. Campi, A. Merkys, A. Marrazzo, T. Sohier, I. E. Castelli, A. Cepellotti, G. Pizzi, _et al._ , Two-dimensional materials from high-throughput computational exfoliation of experimentally known compounds, Nature nanotechnology 13, 246 (2018).
  * Saal  _et al._ (2013) J. E. Saal, S. Kirklin, M. Aykol, B. Meredig, and C. Wolverton, Materials design and discovery with high-throughput density functional theory: the open quantum materials database (oqmd), Jom 65, 1501 (2013).
  * Kirklin  _et al._ (2015) S. Kirklin, J. E. Saal, B. Meredig, A. Thompson, J. W. Doak, M. Aykol, S. Rühl, and C. Wolverton, The open quantum materials database (oqmd): assessing the accuracy of dft formation energies, npj Computational Materials 1, 1 (2015).
  * Jain  _et al._ (2013) 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 1, 011002 (2013).
  * Jain  _et al._ (2016) A. Jain, Y. Shin, and K. A. Persson, Computational predictions of energy materials using density functional theory, Nature Reviews Materials 1, 1 (2016).
  * Meng and Arroyo-de Dompablo (2009) Y. S. Meng and M. E. Arroyo-de Dompablo, First principles computational materials design for energy storage materials in lithium ion batteries, Energy & Environmental Science 2, 589 (2009).
  * Greeley  _et al._ (2006) J. Greeley, T. F. Jaramillo, J. Bonde, I. Chorkendorff, and J. K. Nørskov, Computational high-throughput screening of electrocatalytic materials for hydrogen evolution, Nature materials 5, 909 (2006).
  * Gautier  _et al._ (2015) R. Gautier, X. Zhang, L. Hu, L. Yu, Y. Lin, T. O. Sunde, D. Chon, K. R. Poeppelmeier, and A. Zunger, Prediction and accelerated laboratory discovery of previously unknown 18-electron abx compounds, Nature chemistry 7, 308 (2015).
  * Zakutayev  _et al._ (2013) A. Zakutayev, X. Zhang, A. Nagaraja, L. Yu, S. Lany, T. O. Mason, D. S. Ginley, and A. Zunger, Theoretical prediction and experimental realization of new stable inorganic materials using the inverse design approach, Journal of the American Chemical Society 135, 10048 (2013).
  * Anand  _et al._ (2019) S. Anand, M. Wood, Y. Xia, C. Wolverton, and G. J. Snyder, Double half-heuslers, Joule 3, 1226 (2019).
  * Wang  _et al._ (2017) J. Wang, K. Hanzawa, H. Hiramatsu, J. Kim, N. Umezawa, K. Iwanaka, T. Tada, and H. Hosono, Exploration of stable strontium phosphide-based electrides: theoretical structure prediction and experimental validation, Journal of the American Chemical Society 139, 15668 (2017).
  * Xi  _et al._ (2018) L. Xi, S. Pan, X. Li, Y. Xu, J. Ni, X. Sun, J. Yang, J. Luo, J. Xi, W. Zhu, _et al._ , Discovery of high-performance thermoelectric chalcogenides through reliable high-throughput material screening, Journal of the American Chemical Society 140, 10785 (2018).
  * Hao  _et al._ (2014) S. Hao, L.-D. Zhao, C.-Q. Chen, V. P. Dravid, M. G. Kanatzidis, and C. M. Wolverton, Theoretical prediction and experimental confirmation of unusual ternary ordered semiconductor compounds in sr–pb–s system, Journal of the American Chemical Society 136, 1628 (2014).
  * Samanta  _et al._ (2020) M. Samanta, K. Pal, U. V. Waghmare, and K. Biswas, Intrinsically low thermal conductivity and high carrier mobility in dual topological quantum material, n-type bite, Angew. Chem. 132, 4852 (2020).
  * Darolia (2013) R. Darolia, Thermal barrier coatings technology: critical review, progress update, remaining challenges and prospects, Int. Mater. Rev. 58, 315 (2013).
  * Matsunaga  _et al._ (2011) T. Matsunaga, N. Yamada, R. Kojima, S. Shamoto, M. Sato, H. Tanida, T. Uruga, S. Kohara, M. Takata, P. Zalden, G. Bruns, I. S. H. C. Wille, R. P. Hermann, and M. Wuttig, Phase-change materials: vibrational softening upon crystallization and its impact on thermal properties, Adv. Funct. Mater. 21, 2232 (2011).
  * Zhao  _et al._ (2016) L.-D. Zhao, G. Tan, S. Hao, J. He, Y. Pei, H. Chi, H. Wang, S. Gong, H. Xu, V. P. Dravid, _et al._ , Ultrahigh power factor and thermoelectric performance in hole-doped single-crystal snse, Science 351, 141 (2016).
  * Biswas  _et al._ (2012) K. Biswas, J. He, I. D. Blum, C.-I. Wu, T. P. Hogan, D. N. Seidman, V. P. Dravid, and M. G. Kanatzidis, High-performance bulk thermoelectrics with all-scale hierarchical architectures, Nature 489, 414 (2012).
  * van Roekeghem  _et al._ (2016) A. van Roekeghem, J. Carrete, C. Oses, S. Curtarolo, and N. Mingo, High-throughput computation of thermal conductivity of high-temperature solid phases: the case of oxide and fluoride perovskites, Phys. Rev. X 6, 041061 (2016).
  * Carrete  _et al._ (2014) J. Carrete, W. Li, N. Mingo, S. Wang, and S. Curtarolo, Finding unprecedentedly low-thermal-conductivity half-heusler semiconductors via high-throughput materials modeling, Physical Review X 4, 011019 (2014).
  * He  _et al._ (2016) J. He, M. Amsler, Y. Xia, S. S. Naghavi, V. I. Hegde, S. Hao, S. Goedecker, V. Ozoliņš, and C. Wolverton, Ultralow thermal conductivity in full heusler semiconductors, Phys. Rev. Lett. 117, 046602 (2016).
  * Belsky  _et al._ (2002) A. Belsky, M. Hellenbrandt, V. L. Karen, and P. Luksch, New developments in the inorganic crystal structure database (icsd): accessibility in support of materials research and design, Acta Crystallographica Section B: Structural Science 58, 364 (2002).
  * Curtarolo  _et al._ (2012) S. Curtarolo, W. Setyawan, S. Wang, J. Xue, K. Yang, R. H. Taylor, L. J. Nelson, G. L. Hart, S. Sanvito, M. Buongiorno-Nardelli, _et al._ , Aflowlib. org: A distributed materials properties repository from high-throughput ab initio calculations, Computational Materials Science 58, 227 (2012).
  * Koscielski and Ibers (2012) L. A. Koscielski and J. A. Ibers, The structural chemistry of quaternary chalcogenides of the type amm’q3, Zeitschrift für anorganische und allgemeine Chemie 638, 2585 (2012).
  * Prakash  _et al._ (2015) J. Prakash, A. Mesbah, J. C. Beard, and J. A. Ibers, Syntheses and crystal structures of baagtbs3, bacugdte3, bacutbte3, baagtbte3, and csagute3, Zeitschrift für anorganische und allgemeine Chemie 641, 1253 (2015).
  * Ruseikina  _et al._ (2017) A. Ruseikina, S. Kol’tsov, O. Andreev, and L. Pimneva, Crystal structure of eulnags 3 (ln= gd and ho) compounds, Russian Journal of Inorganic Chemistry 62, 1632 (2017).
  * Ruseikina  _et al._ (2019a) A. V. Ruseikina, L. A. Solovyov, V. . Chernyshev, A. S. Aleksandrovsky, O. V. Andreev, S. N. Krylova, A. S. Krylov, D. . Velikanov, M. S. Molokeev, N. G. Maximov, _et al._ , Synthesis, structure, and properties of euercus3, Journal of Alloys and Compounds 805, 779 (2019a).
  * Ruseikina  _et al._ (2019b) A. V. Ruseikina, L. A. Solovyov, M. V. Grigoriev, and O. V. Andreev, Crystal structure variations in the series srlncus3 (ln= la, pr, sm, gd, er and lu), Acta Crystallographica Section C: Structural Chemistry 75, 584 (2019b).
  * Sikerina and Andreev (2007) N. Sikerina and O. Andreev, Crystal structures of srlncus 3 (ln= gd, lu), Russian Journal of Inorganic Chemistry 52, 581 (2007).
  * Maier  _et al._ (2016) S. Maier, J. Prakash, D. Berthebaud, O. Perez, S. Bobev, and F. Gascoin, Crystal structures of the four new quaternary copper (i)-selenides a0. 5cuzrse3 and acuyse3 (a= sr, ba), Journal of Solid State Chemistry 242, 14 (2016).
  * Azarapin  _et al._ (2020) N. O. Azarapin, A. S. Aleksandrovsky, V. V. Atuchin, T. A. Gavrilova, A. S. Krylov, M. S. Molokeev, S. Mukherjee, A. S. Oreshonkov, and O. V. Andreev, Synthesis, structural and spectroscopic properties of orthorhombic compounds balncus3 (ln= pr, sm), Journal of Alloys and Compounds 832, 153134 (2020).
  * Strobel and Schleid (2006) S. Strobel and T. Schleid, Three structure types for strontium copper (i) lanthanide (iii) selenides srcumse3 (m= la, gd, lu), Journal of alloys and compounds 418, 80 (2006).
  * Hao  _et al._ (2019) S. Hao, L. Ward, Z. Luo, V. Ozolins, V. P. Dravid, M. G. Kanatzidis, and C. Wolverton, Design strategy for high-performance thermoelectric materials: The prediction of electron-doped kzrcuse3, Chemistry of Materials 31, 3018 (2019).
  * Pal  _et al._ (2019a) K. Pal, Y. Xia, J. He, and C. Wolverton, Intrinsically low lattice thermal conductivity derived from rattler cations in an amm′ q3 family of chalcogenides, Chem. Mater. 31, 8734 (2019a).
  * Pal  _et al._ (2019b) K. Pal, X. Hua, Y. Xia, and C. Wolverton, Unraveling the structure-valence-property relationships in amm′ q3 chalcogenides with promising thermoelectric performance, ACS Applied Energy Materials 3, 2110 (2019b).
  * Pal  _et al._ (2019c) K. Pal, Y. Xia, J. He, and C. Wolverton, High thermoelectric performance in baagyte 3 via low lattice thermal conductivity induced by bonding heterogeneity, Physical Review Materials 3, 085402 (2019c).
  * Hautier  _et al._ (2012) G. Hautier, A. Jain, and S. P. Ong, From the computer to the laboratory: materials discovery and design using first-principles calculations, Journal of Materials Science 47, 7317 (2012).
  * Sun  _et al._ (2016) 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 2, e1600225 (2016).
  * Aykol  _et al._ (2018) M. Aykol, S. S. Dwaraknath, W. Sun, and K. A. Persson, Thermodynamic limit for synthesis of metastable inorganic materials, Science advances 4, eaaq0148 (2018).
  * Jain  _et al._ (2011) A. Jain, G. Hautier, C. J. Moore, S. P. Ong, C. C. Fischer, T. Mueller, K. A. Persson, and G. Ceder, A high-throughput infrastructure for density functional theory calculations, Computational Materials Science 50, 2295 (2011).
  * Ward and Wolverton (2017) L. Ward and C. Wolverton, Atomistic calculations and materials informatics: A review, Current Opinion in Solid State and Materials Science 21, 167 (2017).
  * Haastrup  _et al._ (2018) S. Haastrup, M. Strange, M. Pandey, T. Deilmann, P. S. Schmidt, N. F. Hinsche, M. N. Gjerding, D. Torelli, P. M. Larsen, A. C. Riis-Jensen, _et al._ , The computational 2d materials database: high-throughput modeling and discovery of atomically thin crystals, 2D Materials 5, 042002 (2018).
  * Emery  _et al._ (2016) A. A. Emery, J. E. Saal, S. Kirklin, V. I. Hegde, and C. Wolverton, High-throughput computational screening of perovskites for thermochemical water splitting applications, Chemistry of Materials 28, 5621 (2016).
  * Li  _et al._ (2014) W. Li, J. Carrete, N. A. Katcho, and N. Mingo, Shengbte: A solver of the boltzmann transport equation for phonons, Comp. Phys. Comm. 185, 1747 (2014).
  * Zhao  _et al._ (2014) L.-D. Zhao, S.-H. Lo, Y. Zhang, H. Sun, G. Tan, C. Uher, C. Wolverton, V. P. Dravid, and M. G. Kanatzidis, Ultralow thermal conductivity and high thermoelectric figure of merit in snse crystals, Nature 508, 373 (2014).
  * Wu  _et al._ (2017) D. Wu, L. Wu, D. He, L.-D. Zhao, W. Li, M. Wu, M. Jin, J. Xu, J. Jiang, L. Huang, M. G. Zhu, Yimei Kanatzidis, and J. He, Direct observation of vast off-stoichiometric defects in single crystalline snse, Nano Energy 35, 321 (2017).
  * Jana  _et al._ (2017) M. K. Jana, K. Pal, A. Warankar, P. Mandal, U. V. Waghmare, and K. Biswas, Intrinsic rattler-induced low thermal conductivity in zintl type tlinte2, J. Am. Chem. Soc. 139, 4350 (2017).
  * Lin  _et al._ (2016) H. Lin, G. Tan, J.-N. Shen, S. Hao, L.-M. Wu, N. Calta, C. Malliakas, S. Wang, C. Uher, C. Wolverton, _et al._ , Concerted rattling in csag5te3 leading to ultralow thermal conductivity and high thermoelectric performance, Angewandte Chemie 128, 11603 (2016).
  * Tadano  _et al._ (2015) T. Tadano, Y. Gohda, and S. Tsuneyuki, Impact of rattlers on thermal conductivity of a thermoelectric clathrate: A first-principles study, Phys. Rev. Lett. 114, 095501 (2015).
  * Li and Mingo (2015) W. Li and N. Mingo, Ultralow lattice thermal conductivity of the fully filled skutterudite ybfe4sb12 due to the flat avoided-crossing filler modes, Phys. Rev. B 91, 144304 (2015).
  * Li  _et al._ (2015) C. W. Li, J. Hong, A. F. May, D. Bansal, S. Chi, T. Hong, G. Ehlers, and O. Delaire, Orbitally driven giant phonon anharmonicity in snse, Nat. Phys. 11, 1063 (2015).
  * Nielsen  _et al._ (2013) M. D. Nielsen, V. Ozolins, and J. P. Heremans, Lone pair electrons minimize lattice thermal conductivity, Energy Environ. Sci. 6, 570 (2013).
  * Green  _et al._ (2017) M. L. Green, C. Choi, J. Hattrick-Simpers, A. Joshi, I. Takeuchi, S. Barron, E. Campo, T. Chiang, S. Empedocles, J. Gregoire, _et al._ , Fulfilling the promise of the materials genome initiative with high-throughput experimental methodologies, Applied Physics Reviews 4, 011105 (2017).
  * Xia  _et al._ (2020a) Y. Xia, K. Pal, J. He, V. Ozoliņš, and C. Wolverton, Particlelike phonon propagation dominates ultralow lattice thermal conductivity in crystalline tl3vse4, Phys. Rev. Lett. 124, 065901 (2020a).
  * Xia  _et al._ (2020b) Y. Xia, V. I. Hegde, K. Pal, X. Hua, D. Gaines, S. Patel, J. He, M. Aykol, and C. Wolverton, High-throughput study of lattice thermal conductivity in binary rocksalt and zinc blende compounds including higher-order anharmonicity, Physical Review X 10, 041029 (2020b).
  * Pal  _et al._ (2020) K. Pal, Y. Xia, and C. Wolverton, Microscopic mechanism of unusual lattice thermal transport in tlinte $\\_2$, arXiv preprint arXiv:2012.04072 (2020).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comp. Mater. Sci. 6, 15 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Cerqueira  _et al._ (2015) T. F. Cerqueira, S. Lin, M. Amsler, S. Goedecker, S. Botti, and M. A. Marques, Identification of novel cu, ag, and au ternary oxides from global structural prediction, Chemistry of Materials 27, 4562 (2015).
  * Wu  _et al._ (2013) Y. Wu, P. Lazic, G. Hautier, K. Persson, and G. Ceder, First principles high throughput screening of oxynitrides for water-splitting photocatalysts, Energy & environmental science 6, 157 (2013).
  * Togo and Tanaka (2015) A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater. 108, 1 (2015).
  * Setyawan and Curtarolo (2010) W. Setyawan and S. Curtarolo, High-throughput electronic band structure calculations: Challenges and tools, Comp. Mater. Sci. 49, 299 (2010).
  * Pailhès  _et al._ (2014) S. Pailhès, H. Euchner, V. M. Giordano, R. Debord, A. Assy, S. Gomès, A. Bosak, D. Machon, S. Paschen, and M. De Boissieu, Localization of propagative phonons in a perfectly crystalline solid, Phys. Rev. Lett. 113, 025506 (2014).
  * Beltukov  _et al._ (2013) Y. Beltukov, V. Kozub, and D. Parshin, Ioffe-regel criterion and diffusion of vibrations in random lattices, Phys. Rev. B 87, 134203 (2013).
  * Chaput  _et al._ (2011) L. Chaput, A. Togo, I. Tanaka, and G. Hug, Phonon-phonon interactions in transition metals, Phys. Rev. B 84, 094302 (2011).
  * Togo  _et al._ (2015) A. Togo, L. Chaput, and I. Tanaka, Distributions of phonon lifetimes in brillouin zones, Phys. Rev. B 91, 094306 (2015).
  * Zhou  _et al._ (2014) F. Zhou, W. Nielson, Y. Xia, and V. Ozoliņš, Lattice anharmonicity and thermal conductivity from compressive sensing of first-principles calculations, Phys. Rev. Lett. 113, 185501 (2014).