Perovskite photovoltaics is a rapidly advancing area of research. The most efficient perovskite solar cells (PSCs) have already surpassed 26% power conversion efficiency [1], [2], which rivals the leading silicon-based alternatives. Moreover, tandem cells that pair perovskites with silicon provide even higher efficiencies [3]. Despite notable advancements, the commercial deployment of PSCs continues to face significant obstacles, primarily due to unresolved issues surrounding long-term operational stability[4], [5] and concerns over lead toxicity in the most efficient perovskite compositions [6], [7].

The properties of perovskite materials can be systematically tailored through compositional engineering [8]. Among these, alloys derived from the hybrid perovskite (FA = formamidinium, ) have emerged as promising candidates for enhancing both the stability and environmental compatibility of PSCs [9]. These materials typically adopt the form , where partial substitution of FA with Cs has improved structural and thermal stability [10], replacement of Pb with Sn at the \(B\)-site offers a pathway to reduced toxicity, and halide mixing at the \(X\)-site enables precise band gap tuning [11]. Despite these advances, the fundamental mechanisms governing the stability of such quaternary alloys remain incompletely understood.

The alloys have so far received the most attention [12]–[15]. Wang et al., for example, conducted the most comprehensive exploration of the two-dimensional compositional space, synthesizing and characterizing 81 distinct compositions [14]. In contrast, the lead-free counterpart remains comparatively underexplored. To date, investigations have been limited to binary systems such as , , and [11], [16]–[18], leaving much of the quaternary compositional landscape uncharted.

Computational materials science offers a powerful route to accelerate the discovery of improved perovskite materials by elucidating the fundamental mechanisms underlying their intrinsic instabilities and enabling predictions for compositions yet to be explored experimentally. However, traditional computational methods face limitations due to the inherent complexity of these multicomponent systems. While density functional theory (DFT), the cornerstone of atomistic modeling, provides accurate insight at the electronic and structural levels, its application becomes computationally prohibitive when navigating the vast configurational space of multicomponent perovskite alloys.

To address the computational cost of DFT in alloy studies, cluster expansions have been widely used to rapidly predict energies by fitting models to a limited set of DFT calculations [19]. The cluster expansion has proven effective for analyzing the stability of binary halide perovskite alloys with \(X\)-site [20], [21] or \(B\)-site [22] substitutions. However, as an on-lattice approach, it cannot capture the distinct orientations of molecular cations in hybrid perovskites. An alternative is the special quasirandom structures (SQS) method [23], which approximates alloy properties using a small number of representative configurations. SQS has been applied, for example, to compute Gibbs free energy curves for binary hybrid perovskites such as and [24]. Yet, it lacks an inherent treatment of finite-temperature effects, and entropy contributions are typically introduced via analytical approximations. Moreover, SQS’s reliance on a few configurations introduces significant biases when modeling complex alloy behavior.

Machine learning (ML) offers a powerful alternative to traditional alloy approaches by enabling accurate predictions for structurally complex materials at a fraction of DFT’s computational cost. Atomistic ML models are now widely used for both inorganic [25], [26] and hybrid [27]–[30] single-component perovskites, with recent extensions to alloy systems. In prior work, we introduced a method using ML-based structure relaxations to sample alloy convex hulls, applying it to identify stable compositions in [31]. Other studies have followed similar paths: one used graph neural networks to model with Cd, Zn, and Br substitutions [32], while others have employed neural network potentials to simulate phase segregation in [33] and [34]. However, these approaches generally neglect finite-temperature entropy effects, limiting their ability to fully capture thermodynamic stability.

In this study, we explore the thermodynamic phase stability of the hybrid perovskite alloys and by constructing their free energy landscapes using ML (workflow in Fig. 1). Building on our previous work for inorganic perovskites, we extend our ML framework to hybrid systems, where orientable organic cations introduce added complexity. We train MACE interatomic potentials [35], [36] on DFT data to enable efficient structure relaxations, supported by an optimized training workflow developed earlier [37]. Finite-temperature effects are incorporated via the Wang-Landau algorithm [38]. To accelerate the extensive sampling required, we introduce secondary ML models trained on MACE-relaxed structures to predict relaxed energies directly from unrelaxed configurations. The workflow is rigorously validated at each stage and benchmarked against available experimental data [14], [16], [17].

Our work is guided by the hypothesis that the reduced thermodynamic stability of tin-based hybrid perovskites arises from fundamental differences in their free energy landscapes compared to their lead-based counterparts. Accordingly, our primary objective is to compare the thermodynamic behavior of and to uncover the origins of this disparity. In doing so, we also aim to address the broader question of whether entropy-driven stabilization can enhance the thermodynamic stability of hybrid perovskite alloys.

1 Results↩︎

Our solution to the long-standing sampling challenge that has prevented computational investigations of quaternary hybrid perovskite alloys in the past follows a multi-stage workflow, summarized in Fig. 1. We first accelerated perovskite structure relaxations by training MACE graph neural network interatomic potentials [35], [36] on DFT data. Although this approach speeds up computations by approximately six orders of magnitude compared to DFT, MACE-based relaxations alone are not fast enough for free energy calculations. They did, however, allow us to generate extensive relaxation data that we then used to fit a secondary set of ML models. These models employ Kernel Ridge regression (KRR) to map unrelaxed atomic structures directly to relaxed energy values, bypassing the need for explicit structure relaxation. This second step accelerated the computations by an additional three orders of magnitude, which finally allowed us to perform the energy sampling necessary for the free energy calculations with the Wang-Landau algorithm [38].

In this section, we first present results from model training, assessing the accuracy of structural relaxations with the fitted MACE potentials, and comparing the predictions of the direct relaxation models to MACE-based relaxations. Then, we present the computed free energy landscapes for various perovskite phases of both and , identifying energetically favorable regions of the composition spaces. We also show the convex hulls and curvatures of the free energy surfaces to analyze thermodynamic stability against phase separation.

1.1 Model training↩︎

The MACE potentials were trained on DFT calculations to accelerate structure relaxations, following the efficient data generation workflow we developed earlier [37] (see Methods for details). The final models were fitted using training sets comprising 2000 DFT calculations for both the Pb- and Sn-based alloys. The data included \(2\times 2\times 2\) perovskite supercell structures from four different phases: \(Pm\overline{3}m,\) \(P4/mbm,\) \(I4/mcm,\) and \(Pnma\) (illustrated in Fig. 2).

Ideally, the relaxation accuracy of these models would be evaluated by relaxing a representative set of structures using both DFT and MACE, followed by a direct comparison of the results. However, this approach was not feasible due to the substantial computational cost associated with DFT relaxations – performing them on a statistically meaningful test set would require an order of magnitude more DFT calculations than those used for training the MACE model. Instead, we evaluated the relaxation accuracy by relaxing 100 random alloy configurations with both MACE models. We then performed single-point DFT calculations on the MACE-relaxed structures and compared the DFT energies to the MACE predictions. The corresponding mean absolute errors (MAEs) are listed in the first column of Table 1 for each phase. The overall MAEs per perovskite formula unit (f.u.) are 2.81 meV/ and 3.07 meV/ for and , respectively. Visualization of the same comparison (Fig. 3) shows that the model accuracy for the relaxed structures remains consistent across the whole energy range.

Next we trained the direct relaxations models on MACE relaxation data. To assess their accuracy, an initial test was conducted following initial data generation. The initial datasets – consisting of relaxation data obtained through clustering and Monte Carlo (MC) sampling – were divided into 80% training and 20% test data. The resulting prediction MAEs on the test set, reported in the second column of Table 1, were 10.0 meV/ on average, which is considerably higher than those observed for the MACE models.

To ensure reliable free energy calculations, it was essential to improve the accuracy of the direct relaxation models, especially in the low-energy regime. To achieve this, we applied an active learning scheme that iteratively reduces prediction errors on low-energy alloy configurations through repeated energy minimization MC simulations. To asses the final models, we utilized them in performing energy minimization simulations at all alloy compositions permitted by the 2\(\times\)​2\(\times\)​2 supercell and compared their energy predictions on the obtained minimum alloy configurations to the corresponding MACE relaxation energies. To guarantee unbiased test results, we made sure that no minimum-energy configurations included in the training sets of the models were used in the evaluation. The results of the test (third column of Table 1) reveal that the MAEs of the direct relaxation models on minimum-energy configurations now average 3.2 meV/ – a considerable improvement compared to the initial errors.

1.2 Free energy landscapes↩︎

Free energies were computed at all compositions permitted by the 2\(\times\)​2\(\times\)​2 supercell we adopted using the Wang–Landau algorithm. Over 98% of the simulations converged within 3000000 iterations. Among the converged runs, 1.2% yielded clear outlier values, which were manually removed from further analysis (detailed information in Supplementary Section S2 F). The Wang–Landau simulations were repeated with a different seed, and the mean absolute difference between the free energy values obtained from the two runs was 0.37 meV/ for the raw data, which reduced to 0.28 meV/ after outlier removal.

We assessed thermodynamic stability using the computed free energy landscapes through three complementary approaches. First, the Helmholtz free energy of mixing provides a direct measure of how energetically favorable a given alloy composition is compared to the pure, unalloyed perovskites. Lower values indicate greater stability, with negative values suggesting no tendency toward phase separation. Second, we constructed convex hulls of the free energy landscapes; compositions lying on the convex hull are predicted to be stable and resistant to decomposition. Finally, we examined the curvature of the free energy surfaces. A positive curvature indicates local convexity, which corresponds to stability against small compositional fluctuations. While all compositions on the convex hull have positive curvature, such values can also appear off the hull, suggesting potential metastability in those regions. It should be noted that the analysis presented here does not address the relative preferences over the structural phases, and thus it does not consider phase separation into other structural phases but only within a given phase.

Temperature affects the free energy landscapes through entropic contributions. Figures 4 and 5 show results obtained at 300 K, which is a reasonable temperature for storing perovskite samples and perovskite solar cells. However, higher temperatures are frequently encountered, e.g., during synthesis. Thus, we analyzed the free energy and its curvatures at the elevated temperature of 150 ° C and show them, for brevity, in Supplementary Figs. S13 and S14.

The free energy landscapes, averaged over simulations using two independent random seeds, are presented in Fig. 4. Results for the \(I4/mcm\) phase are excluded due to inconsistencies observed in structure relaxations caused by its structural proximity to the \(P4/mbm\) phase, which led to unreliable free energy values. Instead, we present the \(P4/mbm\) results here as representative of the tetragonal phases, while the results for \(I4/mcm\) are available in Supplementary Section S2 G. For the three remaining phases, both and exhibit similar overall trends with some notable differences.

In the \(Pm\overline{3}m\) phase, features a low-energy region in the middle of the 2D composition space, with a well depth of approximately 20 meV/ and a minimum located around 25% Cs and 30% Br. In contrast, for the lowest energies are reached at the Cs-free edge of the composition space without the well in the middle of the composition range observed for .

For the \(P4/mbm\) phase, the key distinction between the Pb- and Sn-based systems lies in the behavior near the binary edge: in the Pb-based system, the low-energy region extends into quaternary alloys with high Cs and low Br concentrations, indicating enhanced stability. This extension, however, is absent in the Sn-based counterpart. The free energy landscapes for \(Pnma\) resemble those of \(P4/mbm,\) with the exception that neither material space exhibits local free energy minima in the central region of the 2D composition space.

We then analyzed the thermodynamic stability further by constructing the convex hulls and curvatures of the free energy surfaces. The convex hulls for the lead- and tin-based composition spaces are presented in Fig. 5 by bright red areas (for quaternary alloys) and bold black lines (for binary alloys). In the \(Pm\overline{3}m\) phase, both material spaces exhibit a stable region on the convex hull near the lower-left corner of the composition map – corresponding to low Br and low Cs concentrations. For this stable region is broader, spanning the whole range of \(A\)-site substitutions and extending up to 60% Br concentration depending on the Cs/FA ratio. Within the \(P4/mbm\) phase, only the system features quaternary alloy compositions on the convex hull, specifically near the top-left corner of the composition space, corresponding to high Cs and low Br content. For \(Pnma,\) exhibits a small stable region at high Br content at approximately 20% Cs concentration. Stable compositions are also present in the corresponding part of the landscape, but the stable region is slightly broader, extending all the way to the corner of the composition space.

Figure 5 also displays the curvature landscapes of the free energy surfaces on and outside the convex hull. For the \(Pm\overline{3}m\) phase of , the highest curvature values – indicating the highest stability – are located near the and corners of the composition space, as well as along the Cs-free edge of the composition space up to 50% Br concentration. Positive curvature values are also observed at compositions outside the convex hull, most notably near the center and upper regions of the composition space (around 50% Br content), indicating possible metastability. For , the highest curvature is found near the corner – also stable according to the convex hull analysis.

In the \(P4/mbm\) phase, exhibits very high curvature values in the proximity of the corner. Moreover, both and display regions of positive curvature in the upper (Cs-rich) halves of the landscapes. In the \(Pnma\) phase, this behavior is mirrored, with metastable compositions in the Pb-based alloy space forming a band at lower Cs content of approximately 25%. For , this band is broader and covers most low-Cs compositions.

2 Discussion↩︎

To evaluate the effectiveness of our ML approach, we first examined the performance of the MACE potentials. The models demonstrated rapid convergence, requiring remarkably few data points – final training sets contained only 2000 atomic structures. The relaxation energy errors for the final MACE potentials were 2.8 meV/ and 3.1 meV/ for and , respectively (Fig. 3). The errors translate to 0.35 meV/atom on average, which compares favorably to previous ML potentials for hybrid perovskites [29], [30] as well as our earlier work on inorganic systems [31], [37]. Other metrics, such as single-point energy and force predictions (Supplementary Sections S2 A and S2 B), were equally good, further confirming the strong suitability of MACE in modeling complex perovskite systems.

By contrast, the direct relaxation models trained on the initial datasets (obtained via clustering and MC sampling) yielded higher prediction errors: 10.2 meV/ for and 9.8 meV/ for (Table 1). Increasing the training set size did not lead to substantial improvements in accuracy (see Supplementary Section S2 C). In our view, tailoring of the direct relaxation model architecture remains a key challenge for future methodological improvement.

Further training of the direct relaxation models with active learning improved predictions greatly in the low-energy regime (Table 1) – the domain most critical for free energy calculations. Our analysis (Supplementary Section S2 E) identifies two mechanisms behind the improvement: 1) enhanced accuracy on actual low-energy configurations, and 2) suppression of unrealistically low energy predictions for configurations that, in reality, lie much higher on the energy landscape. Without active learning, the high errors in the low-energy regime would have significantly degraded the quality of the free energy results, highlighting its essential role in the workflow.

Together, the two ML stages accelerated energy predictions by approximately nine orders of magnitude, enabling the evaluation of over two billion relaxation energies necessary to construct the free energy landscapes for both and . The findings of this study are intended to guide future experimental efforts in perovskite screening and compositional design. To support this goal, we evaluated how well our results align with existing experimental observations. We compared our findings for to the experimental phase mapping and photostability study by Wang et al. [14], which succeeded in sampling this composition space in full with a high-throughput setup. Our predictions for the cubic phase correspond well to the experiments, with the majority of the stable compounds concentrated in the low-Br, low-Cs region of the composition space. The match for the orthorhombic phase is inconclusive – partially because phase separation as such was not the main objective of the work done by Wang et al. Our detailed analysis (in Supplementary Section S2 I) suggests that future experimental phase mapping and phase separation analyses would be beneficial especially in high-Cs and high-Br regions to confirm if orthorhombic phase is observed. Our comparison also motivates future experimental investigations into the dependence of the material’s properties on the synthesis procedure. Understanding these effects would benefit understanding perovskite stability in both short and long-term.

Although the Sn-based system holds promise for reducing the toxicity of perovskite solar cells, experimental stability data remain limited and are currently available only for the binary alloy . Gao et al. reported successful fabrication of alloy compositions with up to 40% of Cs [16], whereas Pansa-Ngat et al. did not observe phase separation at Cs concentrations below 10% [17]. Our results are in qualitative agreement with both studies since we predict only low concentrations of Cs to be stable (up to 30% of Cs in cubic phase at 300 K, and a stability region moderately expanded from this for 150 ° C).

Our results provide further insight into the quaternary alloys. While its free energy landscapes closely resemble those of in both shape and magnitude, the absence of central low-energy regions indicates lower stability of Sn-based quaternary alloys. Consistently, the convex hull results (Fig. 5) exhibit a broader stability range for , suggesting fewer options for compositional engineering in the Sn-based alloy. The only exception occurs in the low-Cs, high-Br region of the composition space in the \(Pnma\) phase, where more Sn-based compounds are predicted to be stable. Experimental evidence, however, suggests that \(Pnma\) is not the dominant phase for these compositions [14]. Aside from this exception, shows no stability in regions of the composition space where is unstable, which indicates that future experimental efforts to optimize could focus on compositions that have already been identified as favorable for the Pb-based counterpart.

Based on the curvature analysis (Figure 5), at room temperature (300 K), the highest stability regions are located near the corners of the composition space, and no additional stability is gained in the center. This limits the extent to which the quaternary perovskite alloys could be stabilized through mixing entropy, which is at its highest in the middle of the composition space. Regarding the temperature dependence of our results, we observe the same overall trend at a typical perovskite synthesis temperature of 150 ° C, as shown in Supplementary Fig. S14.

In this work, the stability of the perovskite alloys was investigated in a specified scope that has implications on both comparing the experimental results to ours and on guiding future experiments. Our analysis focused on the intrinsic stability of perovskites against phase separation, without directly addressing decomposition mechanisms driven by external stresses such as moisture or air exposure. Furthermore, we computed the free energies of mixing within each phase separately, which restricts the analysis to cases of phase separation where the decomposition products retain the same lattice type as the original compound.

Direct comparison of structural phases in perovskites requires accounting for vibrational entropy contributions, as neglecting them leads to a consistent preference for the lowest-symmetry \(Pnma\) phase due to its lower internal energy. Including vibrational entropy would also enable meaningful comparisons with non-perovskite phases, such as the optoelectronically inactive \(\delta\)-phase, offering deeper insight into decomposition pathways from active to inactive materials. While ML offers a promising route for efficiently estimating vibrational entropies through accelerated phonon calculations, the interatomic potentials developed in this work were not sufficiently accurate for that purpose. We are currently working to enhance our workflow to support efficient training with higher-tier DFT calculations, thereby enabling reliable phase comparisons.

3 Methods↩︎

3.1 DFT calculations↩︎

The DFT calculations for the MACE model data generation were carried out using the all-electron code FHI-aims [39]–[42], with the Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional [43]. To account for van der Waals interactions, we employed a non-local many-body dispersion correction [44]. Brillouin-zone integrations were performed using 4\(\times\)​4\(\times\)​4 and 8\(\times\)​8\(\times\)​8 \(\Gamma\)-centered \(k\)-point meshes for the 2\(\times\)​2\(\times\)​2 and 1\(\times\)​1\(\times\)​1 perovskite supercells, respectively. All calculations used the standard tier-2 basis sets and "tight" integration grids provided by FHI-aims. Scalar relativistic effects were treated using the zeroth-order regular approximation [45].

3.2 MACE interatomic potentials↩︎

We employ MACE, an equivariant graph neural network model [35], [36], to fit ML interatomic potentials for structural relaxation of and hybrid perovskite structures. The DFT training data for the models were generated through a data generation workflow that we developed in our previous work [37]. This workflow utilizes clustering and active learning for diverse and efficient sampling of the structural space. For this study, only minor modifications were made to the workflow: treatment of orientable molecules was implemented to account for the FA cations, and an uncertainty estimation approach utilizing an ensemble of ML models was introduced into the active learning process to enable the transition from our previous regression model to MACE. See Supplementary Section S1 A for a detailed description of the entire training process.

To preserve the intended lattice type during structural relaxations we imposed phase-specific constraints on the lattice parameters and the halide positions. For the cubic \(Pm\overline{3}m\) phase, we forced the simulation cell to maintain a cubic shape by allowing only isotropic scaling, and prohibited any tilting of the coordination octahedra. The tetragonal \(P4/mbm\) and \(I4/mcm\) phases were treated with a different set of constraints: octahedral tilting was permitted only around the \(c\)-axis, and the cell was allowed to deform only by adjusting the height-to-width ratio. No constraints were applied to the orthorhombic \(Pnma\) phase. Structure relaxations with both DFT and MACE were performed using the ASE package [46], employing the Broyden–Fletcher–Goldfarb–Shanno (BFGS) minimizer [47]. The relaxations had a convergence criterion of 5 meV/Å on the maximum atomic force.

3.3 Direct relaxation models↩︎

Due to the high complexity of the configuration space, the MACE-based structure relaxations are not fast enough for free energy computations. To overcome this, we accelerated the energy predictions further by using the MACE potentials to generate relaxation data for training a secondary set of ML models that predict relaxed energies directly from the unrelaxed structures – bypassing the need for explicit relaxations during the free energy sampling. We carried out this structure-to-energy mapping employing an ML model architecture that first constructs a global structural representation of an atomic geometry based on the neural network node features extracted from the the fitted MACE models. The structural representation is then mapped to an energy value using Kernel Ridge regression (KRR).

The data generation strategy for these models was based on the same clustering and active learning workflow used for the MACE potentials, with two key modifications. First, during initial data generation, MC sampling was introduced alongside clustering to explore low-energy structures more efficiently. Secondly, the acquisition strategy of the active learning stage was modified to prioritize low-energy alloy configurations over structures with high prediction uncertainty. A detailed description of the direct relaxation model architecture and model training is available in Supplementary Section S1 B.

3.4 Stability analysis via free energy calculations↩︎

Using the direct relaxation models, we were able to analyze the thermodynamic stability across the composition spaces via computing the free energies and their curvatures. To access the free energies of the perovskite alloys, we employed the Wang-Landau algorithm (see Supplementary Section S1 C for details on our implementation), which is an MC method for sampling the temperature-independent density of states (DOS) of a system [38].

The algorithm estimates the DOS with a discrete grid function, \(\rho(E),\) which allows approximating the partition function by summing over the energy grid: \[\begin{align} Z \approx \sum_{i}\rho(E_i)e^{E_i/k_BT}, \end{align}\] where \(T\) is the simulated temperature. The Helmholtz free energy can then be computed as: \[\begin{align} F = -k_BT\ln Z. \end{align}\] To facilitate interpretation, we calculated the Helmholtz free energy of mixing \[\begin{align} \Delta F_{\text{mix}} =& F_{\ce{(Cs/FA)B(Br/I)3}} - xyF_{\ce{CsBBr3}}\\ &- x(1-y)F_{\ce{CsBI3}} - (1-x)yF_{\ce{FABBr3}}\\ &- (1-x)(1-y)F_{\ce{FABI3}}, \addtocounter{equation}{1}\end{align}\] where \(x\) and \(y\) denote the elemental concentrations in . We performed the Wang-Landau simulation at every allowed alloy composition within the 2\(\times\)​2\(\times\)​2 supercell to construct the free energy landscapes.

The local stability of an alloy composition against phase separation is determined by the curvature of the free energy surface at that point. We estimated the curvatures from the simulated free energy values by first fitting a cubic spline \(f(x,y)\) on the data, and then solving the eigenvalues of the Hessian matrix \[\begin{align} H_f(x, y) &= \begin{bmatrix} \frac{\partial^2 f}{\partial x^2} & \frac{\partial^2 f}{\partial x \partial y} \\ \frac{\partial^2 f}{\partial y \partial x} & \frac{\partial^2 f}{\partial y^2} \end{bmatrix}. \end{align}\] The smaller of the two eigen values of \(H_f\) corresponds to the minimum normal curvature of the free energy surface at concentration (\(x\), \(y\)), reflecting the energetic tendency of phase separation along the most favorable direction. A positive curvature indicates local convexity of the free energy surface, and thus local stability against infinitesimal perturbations.

In addition to local stability, we assessed global thermodynamic stability by determining the convex hulls of the free energy surfaces using the spline fits. A composition lying on the convex hull is thermodynamically stable against phase separation into any other configuration of the same phase. The convex hulls for the binary alloys were determined separately by only considering decomposition along a single alloy dimension. Thus, they can suggest broader stability regions than those indicated by the full two-dimensional convex hull analysis.

4 Data availability↩︎

We made the DFT data generated for MACE model fitting and testing available on NOMAD (https://doi.org/10.17172/NOMAD/2026.03.17-1). The same data in processed form as well as final fitted MACE models for both alloy systems are available through Zenodo (https://doi.org/10.5281/zenodo.18957344).

5 Code availability↩︎

The codes used for all computational steps have been uploaded to GitLab
(https://gitlab.com/cest-group/learnsolar-hybrids).

References↩︎