††thanks: These two authors contribute equally.††thanks: These two authors
contribute equally.

# Effects of model size in density-functional-theory study of alloys: A case
study of $\text{CsPbBr}_{2}\text{Cl}$

Fang Pan  State Key Laboratory for Manufacturing Systems Engineering;
Electronic Materials Research Laboratory, Key Laboratory of the Ministry of
Education, School of Electronic Science and Engineering, Xi’an Jiaotong
University, Xi’an 710049, China  Lin Yang  State Key Laboratory for
Manufacturing Systems Engineering; Electronic Materials Research Laboratory,
Key Laboratory of the Ministry of Education, School of Electronic Science and
Engineering, Xi’an Jiaotong University, Xi’an 710049, China  Zhuangde Jiang
State Key Laboratory for Manufacturing Systems Engineering & International
Joint Laboratory for Micro/Nano Manufacturing and Measurement Technology,
Xi’an Jiaotong University, Xi’an 710049, China  Wei Ren  State Key Laboratory
for Manufacturing Systems Engineering; Electronic Materials Research
Laboratory, Key Laboratory of the Ministry of Education, School of Electronic
Science and Engineering, Xi’an Jiaotong University, Xi’an 710049, China  Zuo-
Guang Ye  Department of Chemistry and 4D LABS, Simon Fraser University,
Burnaby, British Columbia V5A 1A6, Canada  Jingrui Li<EMAIL_ADDRESS>State Key Laboratory for Manufacturing Systems Engineering; Electronic
Materials Research Laboratory, Key Laboratory of the Ministry of Education,
School of Electronic Science and Engineering, Xi’an Jiaotong University, Xi’an
710049, China

###### Abstract

The primary challenge of density-functional-theory exploration of alloy
systems concerns the size of computational model. Small alloy models can
hardly exhibit the chemical disorder properly, while large models induce
difficulty in sampling the alignments within the massive material space. We
study this problem with the $\upgamma$ phase of the mixed halide inorganic
perovskite alloy $\text{CsPbBr}_{2}\text{Cl}$. The distribution of alloy
formation energy becomes narrower when the size of the model system increases
along $\,\sqrt[]{2}\times\,\sqrt[]{2}\times 2$, $2\times 2\times 2$, and
$2\,\sqrt[]{2}\times 2\,\sqrt[]{2}\times 2$ models. This is primarily because
the distribution of Br distribution parameters, which plays a leading role in
determining the formation energy range, is more narrow for larger models. As a
result, larger entropy stability effect can be observed with larger models
especially at high temperatures, for which the approximation using mixing
entropy based on the ideal solution model becomes better.

##  I Introduction

Halide perovskites are promising candidates for next-generation
optoelectronics owing to their outstanding electronic and optical properties,
as well as the common starting materials and low synthesis cost [1]. For
example, the record power conversion efficiency of single-junction perovskite
solar cells has broken the $26\%$ mark [2, 3], catching up with the
conventional, expensive, and still market-dominating single-crystalline
silicon devices [4, 5]; the conversion efficiency against indoor white light
has reached $\sim 45\%$ [6, 7]; high-brightness, high external quantum
efficiency, and excellent monochromaticity have been achieved with perovskite
light-emitting diodes (PeLEDs) [8, 9, 10]; and perovskite-based X-ray
detectors recently exhibit high detection sensitivity, low detection limit,
and low attenuation capability [11, 12, 13]. These emergent techniques are on
their way to industrialization and commercialization, especially providing
that the major barriers being overcome or circumvented, primary the
instability of perovskite materials and devices against heat, moisture, and
oxygen.

Alloying strategy is important and common in designing halide perovskite
materials, as the many isostructural $\text{ABX}_{3}$ (A, B, and X denote the
monovalent cation, bivalent cation, and halide anion, respectively) members
offer a large flexibility in engineering the perovskite alloy composition
toward desired materials properties [14, 15, 8]. Specifically,
$\text{FA}_{1-x-y}\text{MA}_{x}\text{Cs}_{y}\text{Pb}(\text{I}_{1-z}\text{Br}_{z})_{3}$
($\text{FA}=(\text{H}_{2}\text{N})_{2}\text{CH}^{+}$, formamidinium; and
$\text{MA}=\text{CH}_{3}\text{NH}_{3}^{+}$, methylammonium) with mixing at
both A and X sites is one of today’s most preferred core material for
perovskite solar cells [16, 17]; B-site substitution (usually Sn for Pb) is a
common mean to reduce the environmental hazard of these emergent techniques
because of lead; X-site alloying is the primary measure to tune the materials
band gap, which is important for the desired emission wavelengths of PeLEDs
[18, 19], and against different artificial light sources for indoor
photovoltaics [20, 21, 22]. Instability issues especially phase segregation
are thus induced into the perovskite materials by alloying [15, 8, 23, 9].

First-principles calculations primarily using density functional theory (DFT)
play an important role in studying the properties of perovskites including
alloys. There are different approaches in previous studies for this purpose,
mostly due to different sizes of model systems. Traverse DFT studies over all
possible configurations (i.e., alignments of mixed anions) can be conducted
for small model systems, such as the $8$-perovskite-unit model (i.e.,
$(\text{ABX}_{3})_{8}$) for $\text{FA}_{1-x}\text{Cs}_{x}\text{SnI}_{3}$ ($8$
A sites) [24] and the $4$-unit models for
$\text{CsPb}(\text{I}_{1-x}\text{Br}_{x})_{3}$ and
$\text{CsPb}(\text{Br}_{1-x}\text{Cl}_{x})_{3}$ ($12$ X sites) [25]. Energy-
minimizing algorithms such as Monte Carlo simulated annealing were employed to
search for the minimal energy configurations of large model systems such as
the $12$-unit model for $\text{FA}_{1-x}\text{Cs}_{x}\text{PbI}_{3}$ ($12$ A
sites) [26] and the $8$-unit model for
$\text{CsPb}(\text{I}_{1-x}\text{Br}_{x})_{3}$ and
$\text{CsPb}(\text{Br}_{1-x}\text{Cl}_{x})_{3}$ ($24$ A sites) [27].
Alternatively, DFT calculations data of sampled configurations were used to
train cluster-expansion models ($4$-unit models for binary mixed-halide
perovskites [28, 29]) or machine-learning (ML) models ($8$-unit models for
binary mixed-halide perovskites [27]). Notably, our group studied all four
phases of both binary alloys $\text{CsPb}(\text{I}_{1-x}\text{Br}_{x})_{3}$
and $\text{CsPb}(\text{Br}_{1-x}\text{Cl}_{x})_{3}$ with traverse DFT
calculations on the $4$-unit model systems (which are the smallest for both
$\upgamma$ and $\updelta$ phases) [25]. For the first time, we have evaluated
the thermodynamic state functions including Helmholtz free energy, internal
energy, and configurational entropy for each phases and composition at
different temperatures based on DFT data, with (i) correcting the conventional
approach using on the mixing entropy based on the ideal solution model, and
(ii) construct the temperature vs. composition phase diagrams for these
alloys.

The biggest challenges for (DFT-based) computational study of alloys concerns
the size and number of the model systems. In principle, we need large enough
computational models that can properly exhibit the long-range order of alloys,
whose short-range order is missing due to the ion or atom mixing (i.e.,
chemical disorder). And we need many configurations of alloyed ions to present
the disorder character. This is, in principle, infeasible at least with pure
DFT approach because of the massive computational resource demanded (note that
the size of materials space grows with the combinatorial number related to the
size of model system). Possible solutions should include: (i) to find the
proper size of computational model system that can balance the computational
cost and the reliability of results, and (ii) to develop advanced property-
evaluating methods (such as ML based on DFT training data) whose accuracy is
at the DFT level. To this end, size effects study of computational model
systems for alloys is important. Note: The special quasirandom structures
method [30, 31] is widely used for modeling alloy materials. It uses Monte
Carlo simulated annealing to minimize the difference between the cluster
correlation functions of the constructed supercell configuration and those of
the random structure. This method pays particular attention to the random
alignment of alloyed ions or atoms. It also encounters the size-effects
problem: the figure of merit of the best matching configuration decreases with
the size of model system.

In this work, we study the inorganic mixed-halide perovskite alloy
$\text{CsPbBr}_{2}\text{Cl}$ in its room-temperature $\upgamma$ (space group
$Pnma$) phase. This is a typical blue-light emitting material for PeLEDs [32,
33] and closely related to candidate materials for harvesting ultraviolet
light in indoor scenarios [20, 34]. Previous computational studies indicate
that it is located at the formation free energy convex hull within a large
temperature range [25], and the fully regular configuration has the lowest
alloy formation energy through the whole
$\text{CsPb}(\text{Br}_{1-x}\text{Cl}_{x})_{3}$ space [28, 29, 27, 25]. To
study the effects of model size, we collect DFT calculation data from randomly
sampled configurations of $2\times 2\times 2$ ($8$-unit) and
$2\,\sqrt[]{2}\times 2\,\sqrt[]{2}\times 2$ ($16$-unit) model systems and
compare them with the traverse DFT data of the smallest
$\,\sqrt[]{2}\times\,\sqrt[]{2}\times 2$ ($4$-unit) model. To estimate the
impact of model size, we pay particular attention to the distribution of alloy
formation energies, and eventually use this information to evaluate the
thermodynamic state functions of each model.

The remainder of this article is organized as follows. In Sec. II, we briefly
describe the investigated model systems and introduce some distribution
parameters that distinguish the configurations. Section III presents the
computational results, based on which a detailed discussion on the model size
effects are conducted in Sec. IV. Finally, Sec. V concludes with a summary.

##  II Computational models of the $\text{CsPbBr}_{2}\text{Cl}$ alloy

###  II.1 Structure models of the $\text{CsPbBr}_{2}\text{Cl}$ alloy

In this paper we focus on the alloy formation energy which is defined as

$\displaystyle E_{i}$
$\displaystyle=E_{\text{total}}(i)-\frac{2}{3}E_{\text{total}}(\text{CsPbBr}_{3})-\frac{1}{3}E_{\text{total}}(\text{CsPbCl}_{3})$
(1)

with $i$ indexing the $\text{CsPbBr}_{2}\text{Cl}$ alloy configuration and
$E_{\text{total}}$ the DFT calculated total energy (per one perovskite unit,
i.e., $5$ atoms). We considered three different models for the $\upgamma$
phase of $\text{CsPbBr}_{2}\text{Cl}$, including the minimal one
$\,\sqrt[]{2}\times\,\sqrt[]{2}\times 2$ and two larger ones, $2\times 2\times
2$ and $2\,\sqrt[]{2}\times 2\,\sqrt[]{2}\times 2$ (Fig. 1). We must notice
that in this context, the “$\upgamma$ phase” means that the configurations
were constructed according to the following protocol: (i) a regular $\upgamma$
phase $\text{CsPbX}_{3}$ structure was generated, which can be characterized
by the $a^{-}b^{-}c^{+}$ Glazer tilting of the halide octahedra
($\text{X}_{6}$); (ii) $\text{Br}^{-}$ and $\text{Cl}^{-}$ ions were randomly
occupying the X sites with the ratio $2:1$. In most cases, these structures do
not exhibit the $Pnma$ symmetry, because (a) for the equivalent sites in
$\text{CsPbX}_{3}$, their local chemical environments are different in
$\text{CsPbBr}_{2}\text{Cl}$ because of the random Br/Cl alignment, and as a
result, (b) the octahedral tilting in the DFT relaxed structure is not as
regular as in the pure $\text{CsPbBr}_{3}$ or $\text{CsPbCl}_{3}\,$.
Nevertheless, they can still be categorized in the $\upgamma$ phase as the DFT
relaxed structures exhibit $a^{-}b^{-}c^{+}$-like octahedral tilting.

FIG 1: Three structure models (from left to right:
$\,\sqrt[]{2}\times\,\sqrt[]{2}\times 2$, $2\times 2\times 2$, and
$2\,\sqrt[]{2}\times 2\,\sqrt[]{2}\times 2$) for the $\upgamma$ phase of
$\text{CsPbBr}_{2}\text{Cl}$. Cs, Pb, in-plane halides, and out-of-plane
halides are colored in green, dark gray, brown, and red, respectively.

The relaxed structure of each configuration and thus its properties are
uniquely determined by the halide alignment. The first parameter is the Br
proportion in the in-plane (IP) halides. Here, in-plane means that the halides
(brown in Fig. 1, labeled by $\bm{ab}$ hereafter) are approximately located in
the plane spanned by the lattice vectors, around which the halide octahedral
tilt out-of-phase. Out-of-plane (OP) halides (red in Fig. 1) are located
between two neighboring $\bm{ab}$ planes. The possible numbers of IP Br atoms
and the corresponding numbers of possible alloy configurations of different
models are listed in Table 1. For example, the total number of Br/Cl
alignments with $6$ IP Br in the $\,\sqrt[]{2}\times\,\sqrt[]{2}\times 2$ is
$C_{8}^{6}C_{4}^{2}=168$.

TABLE 1: IP and OP halide distributions in the three considered models. For each model, listed are possible numbers of IP $\text{Br}^{-}$ ions, the number of configurations at each IP Br number, and its proportion in the total configurations. $\,\sqrt[]{2}\times\,\sqrt[]{2}\times 2$ |  $2\times 2\times 2$ | $2\,\sqrt[]{2}\times 2\,\sqrt[]{2}\times 2$
---|---|---
Number of | Number of |  | Number of | Number of |  | Number of | Number of |
IP $\text{Br}^{-}$ ions | configurations |  Proportion | IP $\text{Br}^{-}$ ions | configurations |  Proportion | IP $\text{Br}^{-}$ ions | configurations | Proportion
$4$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}70$ |  $0.1414$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}8$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}12870$ |  $0.0175$ | $16$ | $6.01\times 10^{8{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}}$ | $0.0003$
$5$ | $224$ |  $0.4525$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}9$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}91520$ |  $0.1244$ | $17$ | $9.05\times 10^{9{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}}$ | $0.0040$
$6$ | $168$ |  $0.3394$ | $10$ | $224224$ |  $0.3049$ | $18$ | $5.66\times 10^{10}$ | $0.0251$
$7$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}32$ |  $0.0646$ | $11$ | $244608$ |  $0.3326$ | $19$ | $1.95\times 10^{11}$ | $0.0863$
$8$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{00}}1$ |  $0.0020$ | $12$ | $127400$ |  $0.1732$ | $20$ | $4.11\times 10^{11}$ | $0.1822$
|  |  | $13$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}31360$ |  $0.0426$ | $21$ | $5.64\times 10^{11}$ | $0.2499$
|  |  | $14$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{00}}3360$ |  $0.0046$ | $22$ | $5.17\times 10^{11}$ | $0.2291$
|  |  | $15$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{000}}128$ |  $0.0002$ | $23$ | $3.21\times 10^{11}$ | $0.1423$
|  |  | $16$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{00000}}1$ |  $1.36\times 10^{-6}$ | $24$ | $1.35\times 10^{11}$ | $0.0600$
|  |  |  |  |  | $25$ | $3.85\times 10^{10}$ | $0.0171$
|  |  |  |  |  | $26$ | $7.26\times 10^{9{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}}$ | $0.0032$
|  |  |  |  |  | $27$ | $8.80\times 10^{8{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}}$ | $0.0004$
|  |  |  |  |  | $28$ | $6.54\times 10^{7{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}}$ | $2.92\times 10^{-5{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}}$
|  |  |  |  |  | $29$ | $2.78\times 10^{6{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}}$ | $1.23\times 10^{-6{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}}$
|  |  |  |  |  | $30$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{00}}59520$ | $2.64\times 10^{-8{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}}$
|  |  |  |  |  | $31$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0000}}512$ | $2.27\times 10^{-10}$
|  |  |  |  |  | $32$ | ${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{000000}}1$ | $4.43\times 10^{-13}$
In total
| $495$ |  $1$ |  | $735471$ |  $1$ |  | $2.25\times 10^{12}$ | $1$

The total number of configurations are $C_{12}^{8}=495$, $C_{24}^{16}=735471$,
and $C_{48}^{32}=2.25\times 10^{12}$ for $\,\sqrt[]{2}\times\,\sqrt[]{2}\times
2$, $2\times 2\times 2$, and $2\,\sqrt[]{2}\times 2\,\sqrt[]{2}\times 2$,
respectively. For $\,\sqrt[]{2}\times\,\sqrt[]{2}\times 2$, the number is
reduced to $44$ because of symmetry thus allowing for a traverse DFT
exploration [25]. While for the other two models, it is infeasible to perform
DFT calculation on each configuration even symmetry can reduce the total
number of configurations by some orders of magnitude. We sampled
configurations of these two models according to the IP-OP distribution of
halides as listed in Table 1.

###  II.2 Halide distribution parameters and machine learning

The large number of mixed halide configurations in the large $2\times 2\times
2$ and $2\,\sqrt[]{2}\times 2\,\sqrt[]{2}\times 2$ structure models motivate
to call of ML techniques. The aim is to rapidly predict the alloy formation
energies of any given mixed halide alignments, to describe which we need a
proper coding.

We can use a set of parameters that describe the halide alignment. The IP/OP
Br proportion, denoted by $p_{\text{Br}}^{\text{IP}}$ hereafter, already
introduced in Fig. 1 and Table 1 provides the information of individual
$\text{Br}^{-}$ ions’ locations. In this regard it can be considered a one-
body parameter. Two-body parameters give the information of the probability of
two specific halides being found with a certain distance. These include

  * •

$p_{\text{SS}}$: IP nearest-neighbor (NN) Br-Br proportion with both
$\text{Br}^{-}$ ions on the same side (SS, i.e., above or below the plane in
terms of $\bm{c}$) of the $\bm{ab}$ plane (Fig. 2 red). The total number of
SS-NN halide pairs of an $n$-unit model is $2n$.

  * •

$p_{\text{DS}}$: IP NN Br-Br proportion with two $\text{Br}^{-}$ ions on the
different sides (DS) of the $\bm{ab}$ plane (Fig. 2 blue). The total number of
DS-NN halide pairs of an $n$-unit model is $2n$.

  * •

$p_{\text{OP}}$: OP NN Br-Br proportion with one IP and one OP $\text{Br}^{-}$
ions (Fig. 2 green). The total number of OP-NN halide pairs of an $n$-unit
model is $8n$.

  * •

$p_{\text{CP}}$: The portion of Br-Br pairs with both IP $\text{Br}^{-}$ ions
at the neighboring $\bm{ab}$ planes (cross-plane or CP), with similar
coordinates along both $\bm{a}$ and $\bm{b}$ (Fig. 2 orange). The total number
of OP-NN halide pairs of an $n$-unit model is $8n$.

  * •

$p_{\text{SC}}$: The portion of Br-Br pairs with both OP $\text{Br}^{-}$ ions
at the same column (SC, i.e., with similar coordinates along both $\bm{a}$ and
$\bm{b}$), connected with one $\text{Pb}^{2+}$ ion (Fig. 2 yellow). The total
number of SC halide pairs of an $n$-unit model is $n$.

  * •

$p_{\text{CC}}$: The portion of Br-Br pairs with both OP $\text{Br}^{-}$ ions
at the nearest two columns (cross-column or CC), with similar coordinates
along $\bm{c}$ (Fig. 2 purple). The total number of SC halide pairs of an
$n$-unit model is $2n$.

The whole set of parameters,
$\\{p_{\text{Br}}^{\text{IP}},p_{\text{SS}},p_{\text{DS}},p_{\text{OP}},p_{\text{CP}},p_{\text{SC}},p_{\text{CC}}\\}$,
are used as the first feature for our ML study.

FIG 2: Br-Br distribution parameters: SS, DS, OP, CP, SC, and CC (colored by
red, blue, green, orange, yellow, and purple, respectively). Also shown is the
coordinates $(\bm{a},\bm{b},\bm{c})$ system.

Another feature for describing the atomic structure of
$\text{CsPbBr}_{2}\text{Cl}$ alloy is the many-body tensor representation
(MBTR) [35]. MBTR is a structural descriptor that considers structural motifs
such as elemental contents ($k=1$), interatomic distances ($k=2$), bond angles
($k=3$), etc. to form a vector representation of an atomic geometry [35, 36].
A previous study shows that the $k=1$ term does not improve the model accuracy
when higher-order terms were included [37], and it is the same for all
configurations considered in this paper. It is generally concluded that the
inclusion of the $k=3$ term produces a minimal improvement in accuracy with
significantly increasing the computational cost [37, 27]. Therefore we only
consider the $k=2$ term in this paper, similar as in previous studies [37,
27]. The second-order MBTR encodes the distances of all pairs of atoms of the
investigated structure by

$\displaystyle\bm{M}_{\text{E}_{1}\text{E}_{2}}(x)$
$\displaystyle=\sum_{l_{1},l_{2}}w_{l_{1}l_{2}}\frac{1}{\varsigma\,\sqrt[]{2\uppi}}\exp\left(-\frac{\Big{(}x-R_{l_{1}l_{2}}^{-1}\Big{)}^{2}}{2\varsigma^{2}}\right)\,,$
(2)

where $R_{l_{1}l_{2}}$ is the interatomic (cartesian) distance between the
$l_{1}$th atom of element $\text{E}_{1}$ and the $l_{2}$ atom of
$\text{E}_{2}\,$, $\varsigma$ is the Gaussian width parameter. The sum runs
over all atoms of element $\text{E}_{1}$ (indexed by $l_{1}$) and
$\text{E}_{2}$ (indexed by $l_{2}$) within a cutoff distance, except that
$l_{1}\neq l_{2}$ when $\text{E}_{1}=\text{E}_{2}\,$. The weighting function
is defined by

$\displaystyle w_{l_{1}l_{2}}$
$\displaystyle=\left\\{\begin{array}[]{ll}\text{e}^{-sR_{l_{1}l_{2}}}\,,&\quad\text{when}\quad
R_{l_{1}l_{2}}\leqslant r_{\text{cutoff}}\,,\cr 0\,,&\quad\text{when}\quad
R_{l_{1}l_{2}}>r_{\text{cutoff}}\,,\end{array}\right.$ (5)

with $s$ controling its decay rate with $R_{l_{1}l_{2}}\,$. In practice, one
can limit all atoms of element $\text{E}_{1}$ within the unit cell without
loss of generality, and search for all atoms of $\text{E}_{2}$ within the
sphere of radius $r_{\text{cutoff}}$ centered at this $\text{E}_{1}$ atom,
thus avoiding double counting.

The thus generated second-order MBTRs are continuous functions of $x$ whose
dimension is the inverse of interatomic distance. They are then discretized at
a set of equidistance grid points over a certain range of $x$ values. In this
way, the MBTRs of all or particular atom pairs construct a vector presentation
$\bm{M}$ of a given atomic structure.

MBTR has many advantages as a descriptor that maps the atomic structures to
Hilbert space elements [35]. A particular merit is that MBTR, as a continuous
function of the atomic coordinates (note:
$R_{l_{1}l_{2}}=\big{|}\left|\bm{R}_{l_{2}}-\bm{R}_{l_{1}}\right|$ with
$\bm{R}_{l}$ denoting the cartesian coordinate of the $l$th atom), allows a
direct evaluation of its gradient. Therefore one can compute the atomic forces
and stress tensor components directly from the ML model [27]. Nevertheless, we
mainly focus on the relationship of how the formation energy of an alloy
structure depends on the alignment of the alloyed element atoms in this paper.
Thus we evaluate the MBTR of each configuration at its starting structure,
which is initiate by linear interpolation between the relaxed
$\text{CsPbBr}_{3}$ and $\text{CsPbCl}_{3}$ structures according to Vegard’s
law (i.e., the starting structures of different configurations only differ in
the order of halide ions).

We correlate the ML features with the energy labels $E_{\text{label}}$ using
kernel-ridge regression (KRR) with the Gaussian kernel $k$ based on previous
studies [37, 27] and tests:

$\displaystyle E_{\text{label}}$
$\displaystyle=\sum_{i}^{N}\beta_{i}k(\bm{s},\bm{s}_{i})\,,$ (6)
$\displaystyle k(\bm{s},\bm{s}_{i})$
$\displaystyle=\exp(-\gamma||\bm{s}-\bm{s}^{\prime}||_{2}^{2})\,,$ (7)

with $\bm{s}$ denoting the feature of a particular structure, $\beta_{i}$
fitting coefficients, and $\gamma$ the width parameter controlling the kernel.
The determination of $\\{\beta_{i}\\}$ is done according to the following
matrix form

$\displaystyle\bm{\beta}$
$\displaystyle=(\bm{K}+\alpha\bm{I})^{-1}\bm{E}_{\text{label}}^{\text{ref}}$
(8)

with $\bm{K}$ the kernel matrix $K_{ij}=k(\bm{s}_{i},\bm{s}_{j})$,
$\bm{E}_{\text{label}}^{\text{ref}}$ the energy labels of the reference data
set, and $\alpha$ the regularization parameter. The hyperparameters, including
$\alpha$ and $\gamma$ and for the MBTR-ML case also $\sigma$, are optimized
during the training of ML models using the Bayesian optimization structure
search (BOSS) code [38].

##  III Results

###  III.1 Benchmarking the DS-PAW calculations

Figure 3 compares the alloy formation energies of all $44$ small-size models
of $\text{CsPbBr}_{2}\text{Cl}$, obtained from structure optimization
calculations using DS-PAW vs. FHI-aims (data from our previous study [25]).
The results by these two softwares are generally similar, with the DS-PAW
results slightly higher than the FHI-aims results. The root-mean-square error
is $2.6\leavevmode\nobreak\ \text{meV}$ (per formula unit, similarly
hereinafter), and the largest deviation is $6.9\leavevmode\nobreak\
\text{meV}$ Specifically, the portions of configurations with which the
deviation is within $[0,1]\leavevmode\nobreak\ \text{meV}$,
$(1,2]\leavevmode\nobreak\ \text{meV}$, and $(2,3]\leavevmode\nobreak\
\text{meV}$ are $\frac{17}{44}=0.39$, $\frac{10}{44}=0.23$, and
$\frac{8}{44}=0.18$, meaning that the result difference between these two
softwares are within $3\leavevmode\nobreak\ \text{meV}$ for $82\%$ of samples.
Overall, linear refinement indicates that $\text{DS-PAW
result}=0.94\leavevmode\nobreak\ \text{FHI-aims
result}+3.0\leavevmode\nobreak\ \text{meV}$, with the coefficient of
determination $R^{2}=0.95$.

FIG 3: Comparison between the alloy formation energies of all $44$ small-size
models calculated by DS-PAW and FHI-aims.

The average expended resources for the DS-PAW and FHI-aims calculations are
$56.6$ and $133.1\leavevmode\nobreak\ \text{CPU\leavevmode\nobreak\ h}$,
respectively. Based on the comparison of computational accuracy and cost
between these two softwares, we can conclude that DS-PAW performs well as it
reaches the accuracy of FHI-aims which is widely used for halide perovskites
[39, 40, 41, 27, 25, 42] but requiring about only $40\%$ CPU hours.

###  III.2 Energy distribution diagrams of model systems of different sizes

Different to that we sampled all possible configurations of the smallest model
system, the configurations of the medium and the large model systems were
randomly sampled. Nevertheless, we expect that these random samples to contain
most of the characters of the energy level distribution. To this end, we use
density of states (DOS) to describe the distribution of alloy formation
energies of the whole ensemble of $\text{CsPbBr}_{2}\text{Cl}$ configurations
for each model system size. It is defined by

$\displaystyle\rho(E)=A\sum_{i}g_{i}\text{e}^{-(E-E_{i})^{2}/\sigma^{2}},$
$\displaystyle\quad\int\rho(E)\operatorname{d}E=N$ (9)

where $A$ is the normalization factor, $g_{i}$ the degree of degeneracy and
$E_{i}$ the alloy formation energy of the $i$th configuration, $\sigma$ the
width parameter (we chose $\sigma=0.5\leavevmode\nobreak\ \text{meV}$), $N$ is
the total number of possible configurations. The smallest system is explored
by a traverse DFT study, with $g_{i}$ obtained from symmetry analysis (see
Ref. [25]). While for either the medium or the large system, the DOS is
constructed based on $100$ randomly sampled configurations. As the sampling
was performed according to a certain distribution of alloy configurations, we
have all $g_{i}=1$.

The DOS scheme (Fig. 4) converts the discrete representation of (alloy
formation) energy levels to a (quasi-)continuum of possibility distribution
over a certain range of energies. It is obvious that as the model system size
increases, the DOS distribution becomes narrower. Specific analysis is as
follows. Firstly, the signals in the range $E<6\leavevmode\nobreak\
\text{meV}$ for the smallest model system vanish in the DOS for both larger
model systems. The minimal-energy peak of the smallest model is at
$-5.2\leavevmode\nobreak\ \text{meV}$, contributed from the fully ordered
structure ($p_{\text{Br}}^{\text{IP}}=1$). The peak centered at
$5.1\leavevmode\nobreak\ \text{meV}$ for the smallest model system is
noticeable because of the high degree of degeneracy ($32$) of the
corresponding configuration where $7$ $\text{Br}^{-}$ anions occupy the IP
sites (i.e., $p_{\text{Br}}^{\text{IP}}=0.875$). However, configurations with
both parameters were not sampled for both larger model systems due to the low
possibilities (e.g., for $p_{\text{Br}}^{\text{IP}}=0.875$, the proportions
are $4.6\times 10^{-3}$ and $2.9\times 10^{-5}$ for the medium and large
models, respectively). Secondly, for the medium model system, configurations
with alloy formation energy ranging between $28$ and $33\leavevmode\nobreak\
\text{meV}$ noticeably contribute to the DOS, while these configurations are
not allowed by the smallest model system meaning that they are relatively
largely disordered. Finally, both low- and high-energy “tails” of the medium
model system DOS vanish for the large system, indicating that the
corresponding configurations were not sampled for the latter. As a result, the
arithmetic mean energy increases from the smallest ($18.0\leavevmode\nobreak\
\text{meV}$) to the medium ($20.1\leavevmode\nobreak\ \text{meV}$) model
system, then slightly decreases when the model system size further increases
($19.2\leavevmode\nobreak\ \text{meV}$).

FIG 4: Density of states (DOS) diagrams calculated with the smallest
($\,\sqrt[]{2}\times\,\sqrt[]{2}\times 2$, in black), medium ($2\times 2\times
2$, in red), and large ($2\,\sqrt[]{2}\times 2\,\sqrt[]{2}\times 2$, in blue)
models.

###  III.3 Thermodynamics properties

From the DOS schemes we can calculate the partition function $Z$ and therewith
the thermodynamic state functions, namely the Helmholtz free energy $F$, the
internal energy $U$, and the configurational entropy $S$, for each model by

$\displaystyle Z$
$\displaystyle=\int\rho(E)\text{e}^{-nE/k_{\text{B}}T}\operatorname{d}E\,,$
(10) $\displaystyle F$ $\displaystyle=-\frac{1}{n}k_{\text{B}}T\ln(Z)\,,$ (11)
$\displaystyle U$ $\displaystyle=\frac{1}{nZ}\int
E\rho(E)\text{e}^{-nE/k_{\text{B}}T}\operatorname{d}E\,,$ (12) $\displaystyle
S$ $\displaystyle=\frac{U-F}{T}$ (13)

with $n$ denoting the number of perovskite units in the model system. The
results are presented in Fig. 5.

FIG 5: Thermodynamic state functions of the smallest (black), medium (red),
and large (blue) model systems calculated from their DOSs. a Helmholtz free
energy. b configurational entropy in terms of its proportion to the value
based on the ideal solution model. c internal energy. Also shown are: in a,
Helmholtz free energies calculated by adding the mixing entropy term based on
the ideal solution model to either the arithmetic mean energy (dashed lines)
or the minimal formation energy (dotted lines); and in c, the arithmetic mean
energy (dashed lines).

The calculated Helmholtz free energy ($F$) lines (Fig. 5a) show a cross at
$T=200\leavevmode\nobreak\ \text{K}$:
$F(\text{small})<F(\text{medium})<F(\text{large})$ when
$T<200\leavevmode\nobreak\ \text{K}$, while the opposite trend is observed
when $T>200\leavevmode\nobreak\ \text{K}$. This can be rationalized by the
alloy formation energy ranges of the sampled configurations: the contribution
of the low-energy configurations of the smallest model system dominates the
DOS at relatively low temperatures, while for the larger model systems, they
were not sampled due to the low probabilities. As temperature increases, high-
energy configurations play increasingly important role. The large
($2\,\sqrt[]{2}\times 2\,\sqrt[]{2}\times 2$) model exhibits the narrowest
$\rho(E)$ distribution, thus showing the largest entropy stabilization effect
as reflected by the smallest $F$ values.

When neglecting the volume change (which is adequate for solid states), we
have $\operatorname{d}F=-S\operatorname{d}T$. At room temperature or higher
($T\geqslant 300\leavevmode\nobreak\ \text{K}$), all three models exhibit an
approximate linear dependence of $F$ on $T$ (Fig. 5a), indicating that the
configurational entropy $S$ is nearly constant within this temperature range.
Figure 5b clearly shows this character. Define the mixing entropy based on the
ideal solution model by [25]

$\displaystyle S^{\prime}$
$\displaystyle=3k_{\text{B}}\bigg{[}\frac{2}{3}\ln\bigg{(}\frac{2}{3}\bigg{)}+\frac{1}{3}\ln\bigg{(}\frac{1}{3}\bigg{)}\bigg{]}\,,$
(14)

Fig. 5b shows the evolution of $S/S^{\prime}$ with temperature. From $300$ to
$700\leavevmode\nobreak\ \text{K}$, the proportions of configurational entropy
in the ideal-solution-model-based mixing entropy of the smallest, medium, and
large models increase from $74.1\%$ to $80.0\%$, from $81.3\%$ to $86.9\%$,
and from $88.8\%$ to $92.2\%$, respectively. The configurational entropy trend
stays at $S(\text{large})>S(\text{medium})>S(\text{small})$ within the whole
temperature range, which is in accordance with the slope trend of $F$ of these
three models shown in Fig. 5. As a final note, the configurational entropy of
the smallest model is significantly smaller than the other two models at
$100\leavevmode\nobreak\ \text{K}$, indicating that only a few configurations
contribute to the thermodynamic properties of the whole ensemble of the
smallest model at low temperatures.

Figure 5a provide an approximate method to evaluate the Helmholtz free energy
of an alloy ensemble: subtracting the entropy term $TS^{\prime}$ obtained from
the ideal solution model (Eq. 14) from the arithmetic mean energy
$E_{\text{avg}}$ (dashed lines). From the definition of internal energy $U$
(Eq. 12), we know that $U$ it is always below $E_{\text{avg}}$ unless all
configurations have exactly the same energy value. $U$ approaches
$E_{\text{avg}}$ as temperature increases, as shown in Fig. 5c. As
$E_{\text{avg}}>U$ and $S^{\prime}>S$, $(E_{\text{avg}}-TS^{\prime})$ might be
a good approximation to $F=U-TS$ as both errors cancel each other to a certain
extend.

Finally, we shortly discuss another way to evaluate the entropy stabilization
effect which was widely used in computational studies [26, 24, 43]:
subtracting the entropy term $TS^{\prime}$ obtained from the ideal solution
model (Eq. 14) from the (quasi)minimal alloy formation energy $E_{\min}$,
which can be obtained from DFT calculations on a series of test configurations
and usually in combination with some minimization algorithm such as Monte
Carlo simulated annealing [28, 29, 44, 27]. According to Fig. 5a, this is a
worse approach to $F$ compared to $(E_{\text{avg}}-TS^{\prime})$. This is
because $E_{\min}<U$ and $S{\prime}>S$ therefore the errors of both terms lead
to a smaller $F$ but can never cancel each other.

##  IV Discussion

The results presented in Sec. III show noticeable size effects of the model
systems of $\text{CsPbBr}_{2}\text{Cl}$ alloy. They can be summarized as: the
larger the model system, the narrower the distribution of alloy formation
energy, the larger the configurational entropy (i.e., the closer to the mixing
entropy based on the ideal solution model), and the smaller the Helmholtz free
energy at room temperature or higher (i.e., the larger the entropy
stabilization effect). In addition, large models are inappropriate to describe
the low-temperature behavior of this alloy because the low-energy
configurations were missed during the sampling.

###  IV.1 Energy distribution

The proportion data listed in Table 1 can be converted into a quasi-continous
distribution function $\varrho(x)$, with
$x=\frac{m}{2n}=p_{\text{Br}}^{\text{IP}}$ denoting the Br portion in all IP
halides. Accordingly we have $\int\varrho(x)\operatorname{d}x=1$ with the
integration running through the whole domain
$x\in\Big{[}\frac{1}{2},1\Big{]}$. Figure 6 shows the evolution of
$\varrho(x)$ function with the model system size $z$ increases. From
$\lim_{n\to+\infty}n!=\,\sqrt[]{2\uppi n}(n/\text{e})^{n}$, we have for the
infinite model limit

$\displaystyle\lim_{n\to+\infty}\varrho(x)$
$\displaystyle=\lim_{n\to+\infty}2n\frac{C_{2n}^{m}C_{n}^{2n-m}}{C_{3n}^{2n}}$
$\displaystyle\sim\exp\Big{(}-2n\Big{[}2(1-x)\ln(1-x)+x\ln(x)+\big{(}x-\textstyle\frac{1}{2}\big{)}\ln\Big{(}x-\textstyle\frac{1}{2}\Big{)}\Big{]}\Big{)}$
(15)

(more precisely, this is available within the major part of the domain
$x\in\Big{[}\frac{1}{2},1\Big{]}$ except when $x$ is very close to the edge
$\frac{1}{2}$ or $1$). Because of the coefficient $2n$ in the exponent, the
exponential part dominates the character of $\varrho(x)$ with the following
character: (i) $\varrho(x)$ exhibits a single, Gaussian-like peak at
$x=\frac{2}{3}$, and (ii) the width of this peak decreases as $n$ increases.
As a result, in the infinite model limit, almost all configurations will be
characterized by $x=\frac{2}{3}$ which is natural because of the
$\text{Br}_{2}\text{Cl}$ stoichiometry.

FIG 6: Configuration distribution with respective to the IP Br proportion
parameter $p_{\text{Br}}^{\text{IP}}$ (labeled by $x$). Data of model systems
of different sizes $z=4,8,16,32,64$ (only the first three considered in this
paper) are colored in black, red, green, blue, and yellow, respectively.

Figure 7 shows the alloy formation energies of sampled configurations at each
IP Br proportions. A general trend that the energy decreases with
$p_{\text{Br}}^{\text{IP}}$ increases can be observed. The limit case of this
trend was reported in previous studies: the minimal energy configuration
exhibits a fully ordered $\text{Br}^{-}$/$\text{Cl}^{-}$ alignment so that all
$\text{Br}^{-}$ anions occupy the IP halide sites [29, 27, 25]; while in the
maximal energy configuration, all OP halide sites are occupied by the larger
halide $\text{Br}^{-}$, and accordingly the $\text{Br}^{-}$ proportion in IP
halide sites is $0.5$ [25].

Such a trend eases our analysis of the model-size dependence of the energy
distribution by combining Figs. 6 and 7. The probability that the
configurations with $p_{\text{Br}}^{\text{IP}}$ parameters close to
$\frac{2}{3}$ are sampled increases with the size of model system. Fig. 7
shows that the energies of these configurations are approximately distributed
within the range $[10,30]\leavevmode\nobreak\ \text{meV}$, in good accordance
with the arithmetic mean energies of all three models as alluded to. In
addition, the energy distribution is narrower for larger model system, partly
because larger model system allows configurations with
$p_{\text{Br}}^{\text{IP}}$ closer to $\frac{2}{3}$ (the
$p_{\text{Br}}^{\text{IP}}$ parameters closest to $\frac{2}{3}$ for the
smallest, medium, and large model systems are $\frac{5}{8}$, $\frac{11}{16}$,
and $\frac{21}{32}$, $0.042$, $0.021$, and $0.010$ away from $\frac{2}{3}$,
respectively).

FIG 7: Alloy formation energy distribution at each IP Br proportion. a
Smallest ($\,\sqrt[]{2}\times\,\sqrt[]{2}\times 2$) model. b Medium ($2\times
2\times 2$) model. c Large ($2\,\sqrt[]{2}\times 2\,\sqrt[]{2}\times 2$)
model. The vertical dashed lines mark
$p_{\text{Br}}^{\text{IP}}=\frac{2}{3}\,$.

In principle, additional DFT calculation data are required to improve and
eventually the DOS. However, this correspond to tremendous demand of
computational resource and is infeasible for large model systems (such as the
“medium” and “large” in this paper). Recently, the MBTR-KRR ML model has been
successfully applied to inorganic mixed halide perovskites [27]. It can
predict not only the formation energy but also the atomic forces of a $2\times
2\times 2$ supercell at the DFT accuracy ($\sim 10^{6}$ and $\sim 10^{4}$
faster, respectively), therefore it can perform structure optimization whose
results are very similar to DFT. This would therefore be a proper choice if we
want to construct the DOS of alloy formation energies in a (quasi)traverse
manner. Nevertheless, it is natural to regard that the total energy of the
relaxed structure of a particular $\text{CsPbBr}_{2}\text{Cl}$ configuration
uniquely depends on the “topological” Br/Cl alignment, i.e., not including
detailed atomic coordinate variation. We consider either the Br distribution
parameters defined in Sec. II.2 or only-halide (i.e., Br–Br, Cl–Cl, and Br–Cl)
second-order MBTR satisfy this requirement.

Based on the DOS character of all three models, we only perform ML study on
the configurations whose $p_{\text{Br}}^{\text{IP}}$ parameters are close to
$\frac{2}{3}\,$. These include: $p_{\text{Br}}^{\text{IP}}=\frac{5}{8}$ and
$\frac{6}{8}$ for the smallest ($\,\sqrt[]{2}\times\,\sqrt[]{2}\times 2$)
model, $p_{\text{Br}}^{\text{IP}}=\frac{10}{16}$ and $\frac{11}{16}$ for the
medium ($2\times 2\times 2$) model, and
$p_{\text{Br}}^{\text{IP}}=\frac{21}{32}$ and $\frac{22}{32}$ for the largest
($2\,\sqrt[]{2}\times 2\,\sqrt[]{2}\times 2$) model. The total number of data
entries, $128$, is rather small because of the limited number of overall
samples. We divided them as $100$ for training and $28$ for testing. Figure 8
shows the ML results using the Br-distribution parameters (Fig. 8a) and
halide-only MBTR (Fig. 8b). Both ML models show good accuracy as reflected by
the relatively small mean absolute error (MAE) values. Br-distribution-
parameters-based ML model has a much lower $R^{2}$ coefficient than the MBTR-
based model, in accordance with its larger MAE. Overall, the MBTR-KRR ML model
performs better, and can fulfill the aim of approximately predicting the
formation energy of the optimized structure of an alloy configuration.

FIG 8: KRR ML model of alloy formation energies of all configurations with
$p_{\text{Br}}^{\text{IP}}\approx\frac{2}{3}\,$. a Using all Br distribution
parameters defined in Sec. II.2 as feature. b Using only-halide MBTR as
feature.

###  IV.2 Thermodynamic properties

From Eqs. 10-13, we can easily understand the observed trend of size effects
on the thermodynamic properties at high temperatures. This is directly based
on the interand when calculating $Z$ using Eq. 10,
$\rho(E)\text{e}^{-nE/k_{\text{B}}T}\,$. For larger $T$ and narrower $\rho(E)$
peak, the difference of $\text{e}^{-nE/k_{\text{B}}T}$ values for energies
within the peak is smaller. In the limit of $\frac{\langle\Delta
E\rangle}{T}\to 0$ ($\langle\Delta E\rangle$ denotes the energy deviation of
the ensemble), we have

$\displaystyle Z$ $\displaystyle=\text{e}^{-n\langle
E\rangle/k_{\text{B}}T}\int\rho(E)\operatorname{d}E=\text{e}^{-n\langle
E\rangle/k_{\text{B}}T}\varOmega\,,$ (16) $\displaystyle F$
$\displaystyle=\langle E\rangle-\frac{1}{n}k_{\text{B}}T\ln(\varOmega)\,,$
(17) $\displaystyle U$ $\displaystyle=\langle E\rangle,$ (18) $\displaystyle
S$ $\displaystyle=\frac{1}{n}\ln(\varOmega)=S^{\prime}$ (19)

with $\varOmega$ denoting the number of configurations. The underlying physics
is straightforward: since most of the configurations have nearly the same
energy, their contribution to the thermodynamic properties of the whole
ensemble are approximately equal. In this case, the ideal solution model is a
good approximation.

The timeperature-evolution of the thermodynamic properties of the smallest
model system can also be analyzed based on
$\rho(E)\text{e}^{-nE/k_{\text{B}}T}$ (Fig. 9). As the structural parameter
distribution (Fig. 6) of this model system is rather flat, the $\rho(E)$ (DOS,
Fig. 4) is wide and structured. As a result, the
$\rho(E)\text{e}^{-nE/k_{\text{B}}T}$ function which exhibit the major
contribution to the partition function $Z$ shows a transition from low- to
high-energy regions as $T$ increases. When the temperature is as low as
$100\leavevmode\nobreak\ \text{K}$, most of the contribution to the ensemble
is provided by the minimal-energy configuration ($E=-5.2\leavevmode\nobreak\
\text{meV}$ or $4E=-20.8\leavevmode\nobreak\ \text{meV}$) despite the low
degree of degeneracy ($1$). Accordingly, the configurational entropy is very
small. At $200\leavevmode\nobreak\ \text{K}$, the contribution from this
minimal-energy configuration rapidly decays, while the major contribution is
from the configuration ($32$-fold degenerate, $E=5.1\leavevmode\nobreak\
\text{eV}$). As temperature further increases, configurations with higher
energies becomes increasingly important, and accordingly the configurational
entropy increases.

In sum, the high-temperature thermodynamic behavior of the inorganic halide
perovskite alloy $\text{CsPbBr}_{2}\text{Cl}$ can be properly evaluated,
because the energy distribution of most of the configurations is relatively
narrow and can already be adequately reflected by our large model that
contains $16$ perovskite units. In contrast, the low-temperature behavior is
more complicated. Minimal-energy configurations are naturally important at low
temperatures. If their proportion is low (especially in large model systems),
the evaluation of low-temperature thermodynamic properties of alloys remains
open and will be subject to future work.

FIG 9: Interand of partition function calculation (Eq. 10) of the smallest
model system at different temperatures.

##  V Conclusion

In summary, we explore the effects of model size on the computational study of
thermodynamic properties of representative all-inorganic mixed-halide alloy
$\text{CsPbBr}_{2}\text{Cl}$. While the distribution of alloy formation
energies for the smallest model is wide and structured, it becomes narrow when
the size of model system increases, with configurations within both low- and
high-energy regions being less sampled. As a result, the Helmholtz free energy
of the alloy ensemble is more negative for larger model at room temperature or
higher, indicating a stronger entropy stability effect. Overall, we observe
obvious size effects of computational model systems for the
$\text{CsPbBr}_{2}\text{Cl}$, which rapidly decay when the model system
becomes larger.

##  VI Method

For DFT structure optimization calculations, we chose the Perdew-Burke-
Ernzerhof exchange-correlation functional for solids (PBEsol) [45] implemented
in the DFT software DS-PAW, which is based on plane wave basis functions and
projected augmented wave pseudo-potentials [46] and distributed over the
Device Studio platform [47]. PBEsol has already proved its suitability in
modeling halide perovskites [48, 49, 50, 27, 42] and was used for the
calculations of $\text{CsPbBr}_{2}\text{Cl}$ alloys [25], whose data were used
as reference in this paper. The plane-wave cutoff was set to
$600\leavevmode\nobreak\ \text{eV}$ based on test calculations. Structural
relaxation was performed using the conjugate gradient method with the
convergence criteria for total energy and atomic forces set to
$10^{-4}\leavevmode\nobreak\ \text{eV}$ and $0.05\leavevmode\nobreak\
\text{eV}\cdot\text{\AA}^{-1}$, respectively. $\Gamma$-centered $8\times
8\times 6$, $6\times 6\times 6$, and $4\times 4\times 6$ $k$-point meshes were
used for the Brillouin zone integration for the
$\,\sqrt[]{2}\times\,\sqrt[]{2}\times 2$, $2\times 2\times 2$, and
$2\,\sqrt[]{2}\times 2\,\sqrt[]{2}\times 2$ model systems, respectively.

##  VII Acknowledgement

The authors thank Hua Dong, Jie Xu, Fang Yuan, and Jarno Laakso for fruitful
discussions. We gratefully acknowledge HZWTECH for providing computational
facilities. We also thank the computational resources provided by Xi’an
Jiaotong University’s HPC platform and the Computing Center in Xi’an. F.P.
thanks Shenbo Yang and Duo Zhang from HZETECH for technical assistance. This
work was supported by the Natural Science Foundation of Shaanxi Province of
China (grant No. 2023-YBGY-447) and the National Natural Science Foundation of
China (grant No. 62281330043).

## References

  * Snaith [2013] Snaith, H. J. Perovskites: The emergence of a new era for low-cost, high-efficiency solar cells. J. Phys. Chem. Lett. 4, 3623–3630 (2013).
  * NRE [2023] _https://www.nrel.gov/pv/cell-efficiency.html https://www.nrel.gov/pv/cell-efficiency.html_ (National Renewable Energy Laboratory: Best research-cell efficiencies, 2023).
  * Zhou et al. [2024] Zhou, J. et al. Highly efficient and stable perovskite solar cells via a multifunctional hole transporting material. Joule 8, 1691–1706 (2024).
  * Ma et al. [2023a] Ma, F., Zhao, Y., Qu, Z. & You, J. Developments of highly efficient perovskite solar cells. Acc. Mater. Res. 4, 716–725 (2023a).
  * Szabó et al. [2023] Szabó, G., Park, N.-G., De Angelis, F. & Kamat, P. V. Are perovskite solar cells reaching the efficiency and voltage limits? ACS Energy Lett. 8, 3829–3831 (2023).
  * Ma et al. [2023b] Ma, Q. et al. Wide-band-gap perovskite solar minimodules exceeding $43\%$ efficiency under indoor light illumination. Device 1, 100174 (2023b).
  * Ma et al. [2024] Ma, Q. et al. One-step dual-additive passivated wide-bandgap perovskites to realize $44.72\%$-efficient indoor photovoltaics. Energy Environ. Sci. 17, 1637–1644 (2024).
  * Lu et al. [2019] Lu, M. et al. Metal halide perovskite light-emitting devices: promising technology for next-generation displays Adv. Funct. Mater. 29, 1902008 (2019).
  * Liu et al. [2021a] Liu, X.-K. et al. Metal halide perovskites for light-emitting diodes. Nat. Mater. 20, 10–21 (2021a).
  * [10] Fakharuddin, A. et al. Perovskite light-emitting diodes. Nat. Electron. 5, 203–216 (2022).
  * Yakunin et al. [2016] Yakunin, S. et al. Detection of gamma photons using solution-grown single crystals of hybrid lead halide perovskites. Nature Photon. 10, 585–589 (2016).
  * Tsai et al. [2020] Tsai, H. et al. A sensitive and robust thin-film x-ray detector using 2D layered perovskite diodes. Sci. Adv. 6, eaay0815 (2020).
  * Zhang et al. [2024] Zhang, Z. et al. Dion-Jacobson 2D perovskite $2\text{AMPPbBr}_{4}$ with type II band alignment for highly sensitive X-ray detection. Adv. Mater. Technol. 9, 2400102 (2024).
  * Correa-Baena et al. [2017a] Correa-Baena, J.-P. et al. Promises and challenges of perovskite solar cells. Science 358, 739–744 (2017a).
  * Correa-Baena et al. [2019] Correa-Baena, J.-P. et al. Homogenized halides and alkali cation segregation in alloyed organic-inorganic perovskites. Science 363, 627–631 (2019).
  * Correa-Baena et al. [2017b] Correa-Baena, J.-P. et al. The rapid evolution of highly efficient perovskite solar cells. Energy Environ. Sci. 10, 710–727 (2017b).
  * Sun et al. [2021] Sun, S. et al. A data fusion approach to optimize compositional stability of halide perovskites. Matter 4, 1305–1322 (2021).
  * Chiba et al. [2018] Chiba, T. et al. Anion-exchange red perovskite quantum dots with ammonium iodine salts for highly efficient light-emitting devices. Nature Photon. 12, 681–687 (2018).
  * Karlsson et al. [2021] Karlsson, M. et al. Mixed halide perovskites for spectrally stable and high-efficiency blue light-emitting diodes. Nature Commun. 12, 361 (2021).
  * Zuo et al. [2019] Zuo, L., Shi, X., Fu, W. & Jen, A. K.-Y. Highly efficient semitransparent solar cells with selective absorption and tandem architecture. Adv. Mater. 31, 1901683 (2019).
  * Feng et al. [2021] Feng, M., Zuo, C., Xue, D.-J., Liu, X. & Ding, L. Wide-bandgap perovskites for indoor photovoltaics. Sci. Bull. 66, 2047–2049 (2021).
  * Tong et al. [2022] Tong, Y. et al. Wide-bandgap organic-inorganic lead halide perovskite solar cells. Adv. Sci. 9, 2105085 (2022).
  * Liu et al. [2021b] Liu, L. et al. A-site phase segregation in mixed cation perovskite. Mater. Rep.: Energy 1, 100064 (2021b).
  * Gao et al. [2018] Gao, W. et al. Robust stability of efficient lead-free formamidinium tin iodide perovskite solar cells realized by structural regulation. J. Phys. Chem. Lett. 9, 6999–7006 (2018).
  * Pan et al. [2024] Pan, F. et al. Mixed-halide perovskite alloys $\text{CsPb}(\text{I}_{1-x}\text{Br}_{x})_{3}$ and $\text{CsPb}(\text{Br}_{1-x}\text{Cl}_{x})_{3}\,$: New insight of configurational entropy effect from first-principles and phase diagrams. Chem. Mater. 36, 3957–3969 (2024).
  * Yi et al. [2016] Yi, C. et al. Entropic stabilization of mixed A-cation $\text{ABX}_{3}$ metal halide perovskites for high performance perovskite solar cells. Energy Environ. Sci. 9, 656–662 (2016).
  * Laakso et al. [2022] Laakso, J., Todorović, M., Li, J., Zhang, G.-X. & Rinke, P. Compositional engineering of perovskites with machine learning. Phys. Rev. Materials 6, 113801 (2022).
  * Yin et al. [2014] Yin, W.-J., Yan, Y. & Wei, S.-H. Anomalous alloy properties in mixed halide perovskites. J. Phys. Chem. Lett. 5, 3625–3631 (2014).
  * Bechtel and van der Ven [2018] Bechtel, J. S. & van der Ven, A. First-principles thermodynamics study of phase stability in inorganic halide perovskite solid solutions. Phys. Rev. Materials 2, 045401 (2018).
  * Zunger et al. [1990] Zunger, A., Wei, S.-H., Ferreira, L. G. & Bernard, J. E. Special quasirandom structures. Phys. Rev. Lett. 65, 353–356 (1990).
  * Wei et al. [1990] Wei, S.-H., Ferreira, L. G., Bernard, J. E. & Zunger, A. Electronic properties of random alloys: Special quasirandom structures. Phys. Rev. B 42, 9622–9649 (1990).
  * Wang et al. [2020] Wang, C. et al. Dimension control of in situ fabricated $\text{CsPbClBr}_{2}$ nanocrystal films toward efficient blue light-emitting diodes. Nature Commun. 11, 6428 (2020).
  * Kim et al. [2021] Kim, Y. C. et al. High-performance perovskite-based blue light-emitting diodes with operational stability by using organic ammonium cations as passivating agents. Adv. Funct. Mater. 31, 2005553 (2021).
  * Chauhan et al. [2022] Chauhan, A., Shrivastav, A. K. & Oudhia, A. Optical, morphological and electrical studies of fully doctor bladed $\text{CsPbBr}_{2}\text{Cl}$-based perovskite thin films. Microelectron. Eng. 258, 111757 (2022).
  * Huo and Rupp [2022] Huo, H. & Rupp, M. Unified representation of molecules and crystals for machine learning. Mach. Learn.: Sci. Technol. 3, 045017 (2022).
  * Himanen et al. [2020] Himanen, L. et al. DScribe: Library of descriptors for machine learning in materials science. Comput. Phys. Commun. 247, 106949 (2020).
  * Stuke et al. [2019] Stuke, A. et al. Chemical diversity in molecular orbital energy predictions with kernel ridge regression. J. Chem. Phys. 150, 204121 (2019).
  * Todorović et al. [2019] Todorović, M., Gutmann, M. U., Corander, J. & Rinke, P. Bayesian inference of atomistic structure in functional materials. npj Comput. Mater. 5, 35 (2019).
  * Li and Rinke [2016] Li, J. & Rinke, P. Atomic structure of metal-halide perovskites from first principles: The chicken-and-egg paradox of the organic-inorganic interaction. Phys. Rev. B 94, 045201 (2016).
  * Li et al. [2018a] Li, J., Järvi, J. & Rinke, P. Multiscale model for disordered hybrid perovskites: The concept of organic cation pair modes. Phys. Rev. B 98, 045201 (2018a).
  * Li et al. [2018b] Li, J. et al. Activation energy of organic cation rotation in $\text{CH}_{3}\text{NH}_{3}\text{PbI}_{3}$ and $\text{CD}_{3}\text{NH}_{3}\text{PbI}_{3}$: Quasi-elastic neutron scattering measurements and first-principles analysis including nuclear quantum effects. J. Phys. Chem. Lett. 9, 3969–3977 (2018b).
  * Li et al. [2024] Li, J. et al. Structural disorder by octahedral tilting in inorganic halide perovskites: New insight with Bayesian optimization. Small Struct. accepted (2024).
  * Wang et al. [2022] Wang, X. et al. Entropy-driven stabilization of multielement halide double-perovskite alloys. J. Phys. Chem. Lett. 13, 5017–5024 (2022).
  * Dalpian et al. [2019] Dalpian, G. M., Zhao, X.-G., Kazmerski, L., & Zunger, A. Formation and composition-dependent properties of alloys of cubic halide perovskites. Chem. Mater. 31, 2497–2506 (2019).
  * Perdew et al. [2008] Perdew, J. P. et al. Restoring the density-gradient expansion for exchange in solids and surfaces. Phys. Rev. Lett. 100, 136406 (2008).
  * Blöchl [1994] Blöchl P. E. Projector augmented-wave method. Phys. Rev. B 50, 17953–17979 (1994).
  * [47] H. Technology, _Device studio_ , URL https://cloud.hzwtech.com/.
  * Yang et al. [2017] Yang, R. X., Skelton, J. M., da Silva, E. L., Frost, J. M. & Walsh, A. Spontaneous octahedral tilting in the cubic inorganic cesium halide perovskites $\text{CsSnX}_{3}$ and $\text{CsPbX}_{3}$ ($\text{X}=\text{F},\text{Cl},\text{Br},\text{I}$). J. Phys. Chem. Lett. 8, 4720–4726 (2017).
  * Bokdam et al. [2017] Bokdam, M., Lahnsteiner, J., Ramberger, B., Schäfer, T. & Kresse, G. Assessing density functionals using many body theory for hybrid perovskites. Phys. Rev. Lett. 119, 145501 (2017).
  * Seidu et al. [2021] Seidu, A., Dvorak, M., Rinke, P. & Li, J. Atomic and electronic structure of cesium lead triiodide surfaces J. Chem. Phys. 154, 074712 (2021).