Data-driven enhancement of cubic phase stability in mixed-cation perovskites

The A-cation of an ABX₃ perovskite, which locates within the cavity of the BX₆ network, can be chosen from the alkali metals, Li⁺, Na⁺, K⁺, Rb⁺ and Cs⁺, or among monovalent extended molecules. In this work, for instance, we have investigated some of the most popular organic cation choices, namely, hydronium (HY⁺), ammonium (AM⁺), sulfonium (SF⁺), phosphonium (PH⁺), hydroxylammonium (HA⁺), methylammonium (MA⁺), hydrazinium (HZ⁺), formamidinium (FA⁺), formamide (FO⁺), ethylammonium (EA⁺), dimethylammonium (DMA⁺), and guanidinium (GUA⁺). In general the A-cation is expected to induce a local deformation through the interaction with the ions of BX₆ octahedra. The deformations are specific to a particular cation, especially when it is non-spherical, and determine both the stability and the equilibrium structure at a given composition. As such, it is important that the dataset used to construct the ML models contains a sufficiently large number of structures presenting distortions. For this reason, we have designed a rigorous workflow for our DFT calculations (see figure 1(a)), which enables us to cover as much distortion space as possible as a function of the $A/A^{^{\prime}}$-cation content.

Zoom In Zoom Out Reset image size

First, we prepare unit cells, where the organic cation is randomly oriented. These initial coordinates are fully relaxed and the various parameters describing the structure are extracted [21, 22]. Then we construct prototypes of two-cation perovskites, namely, $A_xA^{^{\prime}}_{1-x}BX_3$ where 0 $\leqslant$  x  $\leqslant$ 1, by using 2 × 2 × 2 supercells containing 40 ionic sites (see table S1 in the supporting information (SI) for a list of such two-cation pairs) (available online at stacks.iop.org/MLST/2/025030/mmedia). These initial structures are also relaxed, and their structural information extracted, as described below.

Such a procedure provided us a total of 1864 single- and two-cation entries. However, we found that some compounds relaxed far from the perovskite crystal structure. Thus, we excluded these non-perovskite compounds from the train and test dataset. The 1780 remaining halide perovskites, which kept the 3-dimensional (3D) connectivity of BX₆ octahedra, were included in the ML train/test dataset. Among the 1780 perovskites, there were 504 mixed-cation perovskites and 1276 pure perovskite; we carried out the perovskites' optimization with six different initial orientation to mimic the cations random arrangements.

All calculations are performed at the generalized-gradient approximation level using the Perdew–Burke–Ernzerhof (PBE) [28, 29] functional, including Tkatchenko-Scheffler dispersive-interaction corrections [30]. A plane-waves energy cutoff of 520 eV is employed in the projector-augmented wave (PAW) [31, 32] method implemented in the VASP code [33–35], alongside with a 2 × 2 × 2 k-points Monkhorst–Pack mesh. All compounds are relaxed until the energies and forces are converged within 10⁻⁷ eV/atom and 0.01 eV/Å, respectively. We quantify the octahedral deformations of the BX₆ octahedra by measuring the octahedral quadratic elongation, λ, the angle variance, σ², and the tilt, θ², of the relaxed structures, as the schematic list is shown in figure 1(b). One can find the equations as well as the additional information on how these quantities are extracted in the section of structural description and figure S1, in the SI.

In addition to the structural parameters describing the local lattice deformations, we have introduced a related thermodynamical quantity, named the distortion enthalpy of the mixture, $\Delta H _\mathrm{c,mix}$. We define $\Delta H _\mathrm{c,mix}$ as

$$\begin{equation} \Delta H _\mathrm{c,mix}(A_xA^{^{\prime}}_{1-x}{BX}_3) = x \Delta H_\mathrm{c}({ABX}_3) + (1-x) \Delta H _\mathrm{c}({A}^{^{\prime}}{BX}_3)+ \Delta H_\mathrm{mix}({A}_x{A}^{^{\prime}} _{1-x}{BX}_3)\:, \end{equation} \tag{ 1 }$$

where $\Delta H_\mathrm{c}({ABX}_3)$ and $\Delta H_\mathrm{c}({A}^{^{\prime}}{BX}_3)$ are the relaxation enthalpies of the single-cation perovskites ABX₃ and ${AB^{^{\prime}}X_3}$, respectively, where $\Delta H _\mathrm{c}$ denotes the enthalpy difference between the cubic structure and the fully relaxed one [22]:

$$\begin{equation} \Delta H_\mathrm{c} = H_\mathrm{cubic}-H_\mathrm{relax}\:. \end{equation} \tag{ 2 }$$

Finally, the enthalpy of mixing of ABX₃ and ${A^{^{\prime}}BX_3}$ in their equilibrium (distorted) structure, $\Delta H _\mathrm{mix}$, takes the usual definition:

$$\begin{equation} \Delta H_\mathrm{mix} \left({A}_x{A}^{^{\prime}}_{1-x}{BX}_3 \right) = \Delta H \left({A}_x{A}^{^{\prime}}_{1-x}{BX}_3 \right) - \left\lbrace x \Delta H \left({A}_x{BX}_3 \right) + (1-x) \Delta H\left({A}^{^{\prime}}_{1-x}{BX}_3 \right) \right\rbrace\:, \end{equation} \tag{ 3 }$$

where ΔH(α) is enthalpy of formation for the α compound. Thus, $\Delta H _\mathrm{c,mix}$ can be explicitly written as

$$\begin{align*} &\Delta H _\mathrm{c,mix}(A_xA^{^{\prime}}_{1-x}{BX}_3) = \\ & = H_\mathrm{relax} ({A}_x{A}^{^{\prime}}_{1-x}{BX}_3)- \{xH_\mathrm{relax}({ABX}_3)+(1-x)H_\mathrm{relax}({A^{^{\prime}}BX}_3) \} +\\ &+x\{H_\mathrm{cubic}({ABX}_3)-H_\mathrm{relax}({ABX}_3)\} +(1-x)\{H_\mathrm{cubic}({A^{^{\prime}}BX}_3)-H_\mathrm{relax}({A^{^{\prime}}BX}_3)\}\:, \end{align*}$$

namely, it is the sum of the enthalpy of mixing of the two perovskites ABX₃ and ${A^{^{\prime}}BX}_3$, and the weighted average of the relaxation enthalpies of their single cation form. Clearly if both ABX₃ and ${A^{^{\prime}}BX}_3$ have a cubic structure, then $\Delta H_\mathrm{c,mix} = \Delta H_\mathrm{mix}$. We note that $\Delta H_\mathrm{c,mix}$ is minimized (becomes more negative) when the two-cation mixture has a lower enthalpy than the components and there is no enthalpy gain in relaxing the pure structure from the cubic ones. If one assumes that there is no elastic interaction between the octahedra of the two-cation perovskite, then $\Delta H_\mathrm{c,mix}$ will describe the tendency of forming a stable cubic mixture. We then expect that a small $\Delta H _\mathrm{c,mix}$ will characterize homogeneous two-cation cubic perovskites, while large values will be indicative of possible phase segregation, as usually observed for multi-cation compounds with high Cs concentrations [16, 36]. Throughout this work, we follow the rather standard procedure in high-throughput electronic structure theory to replace the enthalpy with the total DFT energy [37].

We built the ML models, now widely used in materials science [38–40], by using the deeplearning and xgboost modules implemented in the 'h2o' library for R [41]. 'h2o.deeplearning' is based on a multi-layer feed-forward artificial neural network, and it is denoted as DL (deep learning) within this article. DL is able to learn through non-linear mapping functions that take as input a set of features. These are then combined via a fully-connected network. In contrast, XGBoost (h2o.xgboost) is based on the tree-boosting method. The learning procedure is determined by the split points and the corresponding ensemble of the classification labels, while this algorithm creates the varying decision trees. A regularized model formalization of XGBoost controls the problem of over-fitting, helping to generate a better predictive model.

The choice of selecting DL and XGBoost is motivated by our previous experience and published work on a dataset of pure (single cation) perovskites. In that work, we compared the performances of various ML algorithms, such as a generalized linear model (GLM), random forest (RF), gradient-boosting machine (GBM), XGBoost, deep learning (DL; feed-forward neural network), in predicting the structural deformation parameters and the material stability [22]. In that case the data set consisted of ABQ₃ chalcogenide (I–V–VI₃) and ABX₃ halide (I–II–VII₃) perovskites, where charge neutrality was achieved by choosing the following elements: B = V^5 +, Nb^5 +, Ta^5 +, Ge^3 +, Sn^3 + and Pb^3 +; Q = O^2 −, S^2 −, Se^2 −, Te^2 −; and X = F⁻, Cl⁻, Br⁻ and I⁻. We found that the linear regression method results in weak models compared to neuron- and tree-based methods due to the non-linear dependence of the chemical compositions' target properties.

In the present work, the model inputs consist of 24 chemical and structural features extracted from tabulated elemental properties, while the target quantities are λ, σ², θ², and $\Delta H _\mathrm{c,mix}$. Only these target variables are obtained from the DFT calculations for each composition. Some additional features describe the organic cations, namely their effective ionic radius and the number of lone pairs. Finally, in the case of two-cation perovskites, $A_xA^{^{\prime}}_{1-x}BX_3$, we also consider the relative $A/A^{^{\prime}}$ concentration, the effective-ionic-radius difference between the A and $A^{^{\prime}}$ cations, and the weight-averaged effective ionic radius, $r^*_{A,\mathrm{avg}} = xr_A+(1-x)r_{A^{^{\prime}}}$. In contrast, each elemental ion is represented by the period/group number, the ionization energy, the electron affinity of B and X, and the electronegativity difference between B and X. Also, we include as structural features the supercell size, the octahedral and the tolerance factor. As a consequence of our selection, all the input features are costless elemental features obtained from the chemical elements' tabulated properties. These properties have been assessed in terms of their correlation and importance for the machine learning training and the impact of any less important features have been verified (see figures S5 and S6 in the SI).

Clearly, the three features concerning only the two-cation perovskites are not defined for the single-cation ones. In order to maintain the same feature-vector size for both single- and two-cation structures, in the former case, we set the minority cation concentration to zero (x = 1), while assigning a randomly chosen cation to $A^{^{\prime}}$. In practice we represent the ABX₃ composition as $A_1A^{^{\prime}}_0BX_3$. Moreover, note also that the same structure can be represented in two equivalent ways. For instance, the mixture Cs_3/4MA_1/4PbI₃ is identical to the mixture MA_1/4Cs_3/4PbI₃. Thus, to reduce the bias on the specific definition for the cation occupation, we balance our data set in such a way to have an equal number of $A_xA^{^{\prime}}_{1-x}BX_3$ structures with $x\gt0.5$ and $x\lt\,0.5$.

We performed a careful hyper-parameter optimization to construct reliable models, and we systematically vary the hyper-parameters over a regular grid. A grid search on a wide range of hyper-parameters was used for DL and XGBoost training by combining several hyper-parameters. For the hyper-parameters of DL, we varied the number of hidden layers to be either 2 or 3, and the number of neurons in each hidden layer was either 128 or 64. The used activation function is the Rectifier activation function with the learning rate in the set of {0.0001, 0.001, 0.01}. Similarly, for XGBoost, we have varied the number of trees from the set {100, 350, 1000}, the maximum tree depth within {3, 5, 7}, and the learning rate in the range {0.001, 0.01, 0.1, 0.5}.

In general, 85% of the calculated single-cation and two-cation perovskite structures are part of the training dataset, while the remaining 15% constitute the test dataset. In order to select the optimized hyper-parameters, we perform 5-fold cross-validation for each hyper-parameter combination until the model is optimized to have deviance of less than 10⁻¹⁰ or there is no improvement in performance for three consecutive epochs.

Before setting the training set size to 85%, we have inspected the ML performance so to evaluate the bias-variance trade-off problem. We split the dataset into two groups by varying the train/test split ratio and optimize the hyper-parameters at every inspection. We find that for the training set larger than 80% of the data, the mean absolute error (MAE) has to be considered converged, as the ML performance convergence is presented in figure S2 in the SI. (Note that the size of the test set was kept constant throughout the test). The optimal training/test split was found to be 85/15.

As a further test, we have also built a more standard regression model using the GLM method. GLM is found to severely underperform ($R^2 \lt\,0.6$, MAE $\gt\,35$ meV ion⁻¹ for $\Delta H _\mathrm{c,mix}$) due to the non-linear relation between the input features and the target quantities. As one can see in figure S3 of the SI, DL and XGBoost always outperform GLM for all prepared training/test datasets in several training trials. Hence, as we move forward, we will only discuss the best-performing DL and XGBoost models and analyze their predictions.

Finally, in closing this section, we wish to revisit the general philosophy of this work. Our aim is to generate a single universal model capable of predicting structure and stability of perovskites made of Group 14 elements (Ge, Sn, Pb), halides (Cl, Br, I), and multiple cations. The variety of these building blocks makes it difficult to perform efficient training using a linear model. In particular, we aim at constructing a model capable of predicting unseen concentrations and compositions, which are unreachable by direct DFT calculations. These concentrations and compositions will be then validated by experiments. As a result, the strategy used here to build continuously improved ML predictors comprises three tightly interconnected phases, namely,

• Phase 1: Build a DFT database and train the ML models.
• Phase 2: Validate the ML predictions via lab-synthesis/characterization of perovskite films.
• Phase 3: Reinforce the DFT database and retrain the ML models by taking into account the feedback from Phase 2.

An early version of the ML models (Phase 1) has been already built where we used only pure perovskites' DFT results [22, 42], while in follow-up work, we have identified, via lab synthesis, where the models could be improved [24]. This current work reports the findings of multiple months of Phase 3 work, enabling us to revise and improve our earlier model after the prediction-to-lab validation of Phase 2. Therefore, we have performed the DFT calculations for perovskites, where the cations are mixed explicitly at different concentrations. This extended dataset enables us to predict the enthalpy of mixing and addresses the issues identified in Phase 2. We systematically improve the models and provide systematic lab-scale validation and testing.

To incorporate dimethylammonium (DMA) with methylammonium (MA) at the cation site of the perovskite structure, we prepare two different solutions inside the glove box, namely, DMAI (0.1 g, Greatcell)/PbI₂ (0.265 g, TCI)/DMSO (450 µl, Sigma Aldrich), leading to DMAPbI₃, and MAI (0.09 g, Greatcell)/PbI₂ (0.265 g, TCI)/DMSO (450 µl, Sigma Aldrich), leading to MAPbI₃. The solutions are then mixed at different concentrations so to yield DMA_x MA_1 − x PbI₃, at the following stoichiometric fractions x = 0, 0.005, 0.01, 0.02, 0.05, 0.1 and 0.15. Then the glass substrates are inserted into the glove box for spin coating. Two pipettes are kept ready with 100 µl of the perovskite solution and 25 µl of chlorobenzene, used as an anti-solvent. The spin coater is programmed for a single run at 1000 and 4000 rpm for 10 and 30 s, respectively. Chlorobenzene is then dropped while the coater is running at a time comprised between 22 and 25 s during the second step. All the samples are then annealed at 100 ^∘C on a hot plate inside the glove box for 1 h.

The octahedral distortion and tilting in cubic perovskites are, in general, determined by the relative size of the cations with respect to the available void position and by the interaction with the surrounding ions in the BX₆ octahedral network [13, 22, 43]. Thus, displacement and tilt of the cations within the voids are necessary to stabilize the perovskite structure. While an ideal cation preserves the cubic perovskite symmetry, smaller and larger than ideal cations cause the void space to shrink or stretch.

We may relate this local structural change to the extended tolerance factor, which is defined as $\tau^* = (r^*_{A,\mathrm{avg}} + r_{X})/\sqrt{2}(r_{B} + r_{X})$. This extended definition is similar to the conventional definition of tolerance factor, with the exception that, for mixed-cation compounds, the ionic radius of the A-site is replaced by the corresponding weighted-average value taken over all the cations in the mixture, namely, $r^*_{A,\mathrm{avg}} = x r_{A} + (1-x) r_{{A}^{^{\prime}}}$. Then, τ* is used as a descriptor in the conventional way so that a composition is likely to crystallize in the cubic phase for τ* = 1, in a distorted-cubic phase for $0.8 \lt\,\tau^* \lt\,1$, while orthorhombic and hexagonal systems are expected for $\tau^* \lt\,0.8$ and $\tau^* \gt\,1$, respectively. As such, the extended tolerance factor can be used as a driving map in the synthesis of mixed-cation compounds [44], where the designing rule is simply that of fine-tuning $r^*_{A,\mathrm{avg}}$ and hence τ*.

However, unfortunately, the dependence of the final compound's stability and structure over the specific cation mixture appears more complicated than what is expected by merely averaging the cation radii, particularly for organic cations due to the non-spherical and polyatomic features. In contrast to monatomic cation, the organic cations interact with the surrounding halides by the non-spherically steric effect and the hydrogen-bonds in addition to their charge balancing role. We investigate this resulting structural deformation in detail while analyzing our DFT dataset.

Figure 2(a) helps us in understanding the influence of the cation size on the octahedral distortions by showing the dependence of λ, σ² and θ² on the extended tolerance factor, τ*. In the figure, we present results for all available $A/A^{^{\prime}}$ mixtures and the BX₃ inorganic sublattices SnBr₃, SnI₃, PbBr₃ and PbI₃. It is clear that the octahedral deformation parameters, λ and σ², depend quadratically on τ*, regardless of the relative concentration of A and $A^{^{\prime}}$. However, we also note that the minimum of such parabolic curves is found at different τ*'s for different compositions; accordingly, we find τ* = 0.90 for the low-$A^{^{\prime}}$ concentration perovskites (perovskite with a significantly more abundant cation type) and τ* = 0.92 for the high-$A^{^{\prime}}$ concentration ones (perovskites with a similar cations abundance). This observation means that the minimum of the distribution moves to higher τ* as the concentration of the second cation $A^{^{\prime}}$ gets larger (more $A^{^{\prime}}$ content is mixed). At the same time, the dependence of the tilt parameter θ² is more straightforward. It increases as τ* is reduced from τ* = 1 (the value attributed to cubic perovskites), displaying an effect associated with the effective size of the A-site cation [45].

Zoom In Zoom Out Reset image size

More details can be extracted from figures 2(b) and (c), where we present both the cell volume and the octahedral tilt, θ², as a function of the $A^{^{\prime}}$ cation fraction for MA-based perovskites, where the inorganic part is either PbI₃ or PbBr₃ and for different $A^{^{\prime}}$ cations. In general, both the volume and θ² depend linearly on (1 − x), although the sign of the linear slope is opposite for these two quantities, namely for those cations where the volume increases (decreases) with (1 − x), θ² gets smaller (larger).

The only exception is for MA_x Cs_1 − x PbBr₃, for which both the volume and θ² decreases with (1 − x), although the variations are actually small. Such anti-correlation between the behavior of the volume and that of θ² indicates that the BX₆ framework can substitute MA with a second cation by either expanding the volume and reducing the octahedral tilt or by contracting the volume and becoming more distorted. The fact reveals that the behavior of molecular cations is not solely determined by their effective ionic radius. For instance, one can note that in MA$_xA^{^{\prime}}_{1-x}$PbI₃ the two cations HZ and Cs contract the cell volume and rotate the PbI₆ octahedra by approximately the same amount at any concentration, despite having significantly different ionic radii ($r_\mathrm{HZ} = 217$ pm and $r_\mathrm{Cs} = 188$ pm). Incidentally the ionic radius of HZ is essentially the same as that of MA ($r_\mathrm{MA} = 217$ pm), meaning that both the volume change and the distortions induced by the substitution of MA with HZ take place regardless of the fact that τ* remains constant. This behavior is attributed to the more complex non-covalent bonding and repulsion exerted by a given organic cation on the inorganic octahedral network. This fact is associated with the number of lone pairs as earlier presented regarding the organic cation's input features [22].

We now turn our attention to the analysis of the phase stability of the mixed-cation perovskites. Figure 3(a) reports the DFT-calculated $\Delta H _\mathrm{c,mix}$ (top panel) as a function of τ* for a broad range of cation compositions. The figure includes all the cation mixtures considered and the data are presented with the convention that $x\gt\,(1-x)$, thus that, for example, the perovskite MA_0.25Cs_0.75PbI₃ is represented as Cs_0.75MA_0.25PbI₃ in the plot. The lower panel of figure 3(a) displays the best-fit to the DFT data. We note that all the $\Delta H _\mathrm{c,mix}(\tau^*)$ curves have a parabolic shape with the minima located in the interval $0.9\lt\,\tau^*\lt\,1$, regardless of the relative cation concentration. Such minimum moves to higher τ* as the relative cation concentration approaches 50:50. Interestingly, $\Delta H _\mathrm{c,mix}$ appears clear that structures with an equal cation mixture tend to be less favorable, nevertheless τ* indicates a cubic structure for a given perovskite yet.

Zoom In Zoom Out Reset image size

We also observe that the minimum value of $\Delta H _\mathrm{c,mix}$ is always below the 34 meV ion⁻¹ line (dashed line in the figure) within the range of $0.85 \lt\,\tau ^* \lt\,1.00$. Furthermore, it gets lower as the relative concentration of the cations moves away from 50:50 (the minima are distributed over a range of about 10 meV ion⁻¹). Recalling the definition of $\Delta H _\mathrm{c,mix}$, the behavior observed here indicates that structures in which one cation type occupies the majority of the available sites are more enthalpically favorable but likely to be more distorted than the equal-proportion mixtures. These, in contrast, seem to prefer a cubic structure. At this point, however, it is worth noting that the spread of the DFT values is pretty large (top panel), much larger than the difference between the best-fit curves. It is then clear that many other factors, in addition to τ*, determine the stability of a given composition.

It is important to note that our dataset may include some local minima structures because we have not explored all possible molecular orientations and orders in our DFT calculations. Similar to Frost et al who performed large-scale calculations employing a 25 × 25 × 1 supercell [46], here we also find distortions of the edge-sharing BX₆ octahedra due to the electric dipole between the neighboring A-sites. However, given the relatively limited size of our supercells, our calculations cannot explore long-range order and may include local minima instead of the global ones. Note that, in general, we observe relatively small energy differences between the total energy of structures relaxed from different initial conditions, indicating that the potential energy surface as a function of the local molecule orientation is relatively flat. This fact gives us confidence that, even if local minima are considered, our data are still suitable for constructing the ML model.

As MA is the most notable organic cation in solution-processed perovskite synthesis, we now focus our analysis on MA$_xA^{^{\prime}}_{1-x}$PbI₃ (where 0 $\leqslant$  x  $\leqslant$ 1) compounds as a function of the type and concentration of $A^{^{\prime}}$. Figure 3(b) presents $\Delta H _\mathrm{c,mix}$ as a function of the concentration of the $A^{^{\prime}}$ cation for all the cations investigated. Note that $\Delta H _\mathrm{c,mix}$ for (1 − x) = 0 (beginning of the plot) corresponds to the deviation of MAPbI₃ from its cubic structure, while for (1 − x) = 1 (end of the plot) one finds the same quantity for $A^{^{\prime}}$PbI₃. In all cases $\Delta H _\mathrm{c,mix}$ at the (1 − x) = 1 endpoint is higher than at (1 − x) = 0, meaning that the deviation energy from the cubic structure of MA is less pronounced than that of all the other cations investigated. Because of this quantity definition, $\Delta H _\mathrm{c,mix}$ generally increases with (1 − x), although the rate of increase is significantly different for different $A^{^{\prime}}$.

For FA, DMA, and GUA, $\Delta H _\mathrm{c,mix}$ fluctuates only moderately as their concentration is ramped up, suggesting that such cations can be easily mixed with MA with arbitrary relative abundance. As $\Delta H _\mathrm{c,mix}$ remains close to its endpoint across the entire composition, one can conclude that the possible mixed perovskites may have a structure close to cubic. Interestingly, the ionic radii of these three cations are all significantly larger than that of MA.

For two other cations, AM and HZ, the situation is similar, namely $\Delta H _\mathrm{c,mix}$ increases linearly with (1 − x), although the relaxation from the cubic phase has much higher energy gain in AMPbI₃ and HZPbI₃ than in MAPbI₃, so that the curves are much steeper. This outcome means that as one moves away from MAPbI₃ such two-cation compositions are likely to form heavily distorted structures. Finally, for Cs and Rb $\Delta H _\mathrm{c,mix}$ is not monotonic with (1 − x), presenting a maximum at around 1 − x = 0.8. Cs is particularly interesting because the endpoints have similar $\Delta H _\mathrm{c,mix}$. This is significantly lower than the values taken for mixed compositions, meaning that Cs-rich MA_x Cs_1 − x PbI₃ are unlikely to be formed.

However, it should be noted that this figure displays only the part of the MA$_xA^{^{\prime}}_{1-x}$PbI₃ dataset of Ge-, Sn-, Pb-halide perovskites. In fact, when one looks at all perovskites made of Group 14 heavy elements and halides, the trend no longer follows a monotonous dependence. It is, in fact, significantly more complex, showing clear non-linearities caused by relationships between the various A, B, and X elements. For this reason, a single universal predictor, capable of learning and describing several $AA^{\prime} BX_3$ element combinations, cannot be formulated by using linear regression models. Note that a non-linear dependence in the formation energy of cubic-perovskite alloys has been reported previously in literature [47], indicating that the observed deviation of $\Delta H_\mathrm{c,mix}$ from what expected from a composition-weighted average should not be surprising.

Let us now turn our attention to the enthalpy of mixing, which is presented for MA$_xA^{^{\prime}}_{1-x}$PbI₃ as a function of (1 − x) in figures 4(a) and (b), for $A^{^{\prime}}$ cations having an effective ionic radius, respectively, smaller and larger than that of MA. In the lower panels, we show the corresponding free energy of mixing also, $\Delta G_\mathrm{mix}$, at 300 K. This is defined as $\Delta G_\mathrm{mix} = \Delta H _\mathrm{mix} - T \Delta S$, where the entropy of mixing is approximated with that of a solid state solution of A and $A^{^{\prime}}$ cations, namely $\Delta S = k_{\mathrm B}[x \ln(x) + (1-x) \ln (1-x)]$, with $k_\mathrm{B}$ being the Boltzmann constant.

Zoom In Zoom Out Reset image size

For cations smaller than MA (figure 4(a)) the enthalpy of mixing is always positive, meaning that there is never an enthalpy gain in forming two-cation structures. Since the cations interact with the surrounding anions mainly by steric forces, the positive enthalpy of mixing is attributed to the local octahedral deformation introduced by the mixtures. When also entropy is considered (see $\Delta G_\mathrm{mix}$), among all the investigated $A^{^{\prime}}$ cations, only AM presents negative free energy of mixing across the entire composition range, a fact that suggests the possibility of forming various MA_x AM_1 − x PbI₃ perovskites, as demonstrated experimentally [48, 49]. There are two other cations displaying a negative $\Delta G_\mathrm{mix}$ at 300 K, namely Cs⁺ and Rb⁺, although for only a restricted number of concentrations around MA-rich compositions.

In contrast, when MA is mixed together with another ion, presenting an ionic radius larger than that of MA, $\Delta H _\mathrm{mix}$ appears in general small, it never exceeds 10 meV ion⁻¹. And, for a few cations at some concentrations, it turns negative (see figure 4(b)). As a result, the corresponding values for the free energy of mixing at 300 K are always negative, indicating that such cations mixtures are likely to be accessible by the synthesis. Notably, MA-based perovskites incorporating a second cation with an ionic radius larger than MA often crystallize into a hexagonal phase [7, 44, 50, 51].

Note that the cubic phase stability criterion, as well as $\Delta H _\mathrm{mix}$, are expressed as total energy differences and, as such, they are sensitive to the possible uncertainty over the total energy. Assessing the error is not simple because enthalpy data are available only for a handful of compounds. At the same time, a benchmark of our calculations obtained by re-calculating the entire dataset with a different functional is simply not practical, as the numerical cost is high. Here we note that typically PBE total energies are accurate for this class of compounds. The error on total-energy differences is small since the error on total energies is not particularly bias against certain compounds. Thus, in general, we expect that our DFT results provide a quantitative description of the materials trends investigated here.

We now discuss the performance of our ML models in predicting the structural properties of the two-cation-mixed perovskites, namely, in predicting λ, θ², σ², and $\Delta H _\mathrm{c,mix}$. Figures 5(a) and (b) compare the values of the target variables predicted by both DL and XGBoost against the DFT-calculated results, while we apply the best model for the test dataset.

Zoom In Zoom Out Reset image size

A visual inspection of figures 5(a) and (b) already suggests that $\Delta H _\mathrm{c,mix}$ is the target variable better predicted by the two ML models. In fact, we find an MAE of 13.0 and 7.5 meV ion⁻¹ for DL and XGBoost, respectively, against a quantity that can take values as large as 400 meV ion⁻¹. Such attribution is confirmed by the R² coefficient, which is found to be 0.967 for XGBoost and slightly lower for DL. Interestingly, the second-best predicted variable is θ², for which the XGBoost R² is also pretty high, 0.960 (the MAE is of the order of 0.4 for a quantity that takes values as large as 12).

In order to examine the contribution of the input features to both $\Delta H _\mathrm{c,mix}$ and θ², we have estimated the feature importance of the models (see figure S5 in the SI). Although there are differences in the ranks depending on the target variable and algorithms, the feature importance estimation shows that the weight-averaged, effective ionic, and tolerance factors are commonly high in the ranks. In contrast, the difference in effective ionic radii between A and $A^{^{\prime}}$ is low. This analysis supports the empirical approaches by which the approximate properties of mixed-ion perovskite are assessed by the weight-average effective ionic radius of the mixed ions instead of the individual radius.

The performance of the two ML models in predicting λ, and σ² is less satisfactory, and we were not able to train models with R² higher than 0.8 (0.7 for σ²). In the meantime, we are probably overlooking additional features of relevance for λ, and σ². These are likely related to the non-spherical symmetry of the organic cations, a fact which is not taken into full account in the chemical input features of the models.

In general, the two models show a similar prediction capability, as confirmed by the similar R² coefficients over the entire target-variable range (see figure 6(a)), with XGBoost slightly outperforming DL. Although the difference between the two algorithms is small, we attribute the slightly better results of the XGBoost model to its bagging property because it preserves the raw dataset's configuration rather than interpolating its trend so that its predictions are discrete values (see also figure S7 in the SI). Thus, we speculate that this XGBoost's accuracy can be attributed to the presentation of the dataset values. In contrast, the DL model would be advantageous in accessing the intermediate and dilute concentrations that are inaccessible to DFT calculation due to the large supercells needed.

Zoom In Zoom Out Reset image size

A better understanding of the possible correlation between the various target variables can be obtained by looking at the Pearson linear correlation coefficients presented in figure 6(b). It is notable that the most correlated variables are the octahedral quadratic elongation, λ, and the angle variance, σ², indicating that these types of distortions are always simultaneously present in two-cation perovskites. Besides, we find that all the structural features, λ, σ² and θ², have a positive linear correlation to the distortion enthalpy of the mixture, $\Delta H_\mathrm{c,mix}$. As the radius-related features are the main contributor to the target variable rather than each element's chemical properties, the ionic size relative to the surrounding ions plays a crucial role in the structural deformation. As expected, this relationship implies that the most distorted perovskites display the largest enthalpy gain, originating from relaxing the cubic structure into a lower symmetry one.

In our previous work, we have used the relaxation enthalpy, defined in equation (2) [22], as a descriptor for the phase stability of mixed-cation perovskites and their tendency to form cubic structures. This descriptor was constructed by interpolating the DFT data obtained for only single-cation perovskites. Taking advantage of our much larger DFT dataset, including several two-cation structures, we now evaluate whether $\Delta H _\mathrm{c,mix}$ is a better quantity to describe the stability of mixed-ion perovskites. Most importantly, $\Delta H _\mathrm{c,mix}$ can now be correlated to the various structural target properties defining the compound. Our results, as predicted by DL, are presented in figure 7(a), where we show a map of $\Delta H _\mathrm{c,mix}$ as a function of the cation radius and of the $A^{^{\prime}}$ concentration for MA$_xA^{^{\prime}}_{1-x}$PbI₃.

Zoom In Zoom Out Reset image size

Such $\Delta H _\mathrm{c,mix}$ map presents a very distinctive feature. It is a relatively smooth function for $r_{A^{^{\prime}}}\gt\,2\$Å, while it drastically increases as a function of the $A^{^{\prime}}$ cation concentration for $r_{A^{^{\prime}}}\lt\,2\$ Å. The boundary is, thus, set close to the ionic radius of MA. This essentially suggests that the incorporation of a second cation together with MA is thermodynamically allowed only at low concentrations and mostly for cations with an ionic radius larger than MA. It is then expected that relatively small cations such as Cs, may incorporate in the perovskite structure, but without forming a solid-state cation solution, namely, they will segregate.

The structural characteristics of mixed MA-based perovskites can then be extracted from the remaining panels, as shown in figures 7(b)–(d). In the small $r_{A^{^{\prime}}}$ region, one finds that all structures are predicted having a relatively large octahedral quadratic elongation, λ, regardless of the secondary cation concentration. One then expects that, if the small-radius cations can be incorporated in a diluted mixture, the resulting structure will display some tetrahedral distortion. In the same region of parameters, the angle variance, σ², remains moderate, while the angle tilt parameter, θ², increases as the concentration increases, following a pattern similar to that of $\Delta H _\mathrm{c,mix}$. In contrast, for $r_{A^{^{\prime}}}\gt\,2\$ Å, one encounters a region extending to concentrations up to 10%, where the octahedral quadratic elongation remains rather close to unity. In the same region, θ² is small, while σ² is in general large. This representation means that for cations with an ionic radius larger than MA, MA$_xA^{^{\prime}}_{1-x}$PbI₃ perovskites can be stabilized in an almost cubic structure, where the strain associated with the large cations can be accommodated into local distortions of the PbI₃ octahedra.

To validate our analysis, we have decided to synthesize a family mixed-cation perovskites, where MA was combined with a cation presenting a significantly larger ionic radius. In particular, we have found DMA to be particularly attractive ($r_{A^{^{\prime}}}$ for DMA is about 270 pm). In fact, for small DMA concentrations, MA_x DMA_1 − x PbI₃ appears to lie in a region of high stability and small distortions.

We have then explored our predictions on the structural deformation and stability of mixed-cation perovskites by fabricating MA_x DMA_1 − x PbI₃ thin films with (1 − x) = 0.0 to 0.15 (0–15 mol% DMA cation). Our thin-film fabrication results are presented in figure 8(a), where we present the XRD analysis across the composition range. According to the XRD spectra, the integrity of the 3D perovskite structure is maintained only up to MA concentrations of 2.0 mol%. Beyond such value, the broadening of the main XRD peaks indicates the presence of significant structural distortion induced by the large DMA cation (note the appearance of a robust peak at 20^∘, which is absent for both pure MAPbI₃ and DMAPbI₃). Since MA_0.85DMA_0.15PbI₃ has been successfully synthesized in powder form before [50], we speculate that the lack of crystallization of thin films with DMA content above 2 mol% is attributed to the use of a glass substrate.

Zoom In Zoom Out Reset image size

Despite the restricted range of DMA content for thin-film growth on glass, high-resolution XRD spectra enable us to trace possible structural phase transitions as a function of the DMA concentration (see figure 8(b)). In particular, we analyze in detail the fine structure of the 2θ peaks at 14.2^∘ and 28.6^∘. Their peak decomposition in the case of MAPbI₃ is illustrated in the top panel of figure 8(b), and it is characteristic of the tetragonal phase. As DMA is incorporated into the structure (lower panel of figure 8(b)), the two main peaks forming the signal at 14.2^∘ and 28.6^∘ first merge in a single peak for 0.5 and 1 mol% DMA and then separate again for 2 mol% DMA content. Moreover, as only the minor diffraction at 2θ = 23.5^∘ can be used to distinguish between those phases [52, 53], the evident disappearance of the pattern at 2 mol% DMA incorporation, as shown in figure S11 in the SI. Thus, we observe a tetragonal to cubic transition at room temperature for small DMA concentrations, followed by a new distorted phase for $x\gt\,0.02$.

This behavior is in agreement with the predictions from our ML models (DL), presented in figure 8(c). The model forecasts that the octahedral quadratic elongation, λ, first decreases at moderate DMA contents, and then increases again, thus confirming our experimental findings. Simultaneously, the angle variance, σ², and tilt, θ², slightly decreases and increases, respectively. More DMA incorporation results in the octahedral tilt reduction as the concentration (1 − x) gets larger. The behavior of $\Delta H _\mathrm{c,mix}$ with (1 − x) seems to confirm the trend described by the structural parameters. Namely, it increases by about 5 meV ion⁻¹ across an extensive composition range, with such increase being more pronounced at small x. Clearly, we do not expect the models to be quantitative capable of describing our experimental finding because they do not include any effects arising from the substrate. Note that if the extended tolerance factor is used to analyze the structure evolution as a function of x, one will predict a monotonic increase from the cubic structure as the DMA concentration is increased, in contrast with experiments. This fact highlights the need to use a range of descriptors, instead of a single one, to predict the structural stability of hybrid compounds incorporating molecular cations.

Taking the example of the MA_x DMA_{(1 − x)}PbI₃, the lowest concentration we could reasonably achieve with DFT correspond to (1 − x) = 12.5% or 1/8 of DMA. On the experimental side, the concentration of interest could be as small as (1 − x) = 2% or 1/50 requiring a simulation cell at least 6.25 times larger. Carrying the DFT calculations on such large systems is highly impractical and poses obstacles in the computing resources since DFT calculation that often scales as $\mathcal{O}(N^3)$, (computational complexity and cost for DFT calculations as a function of the number of electrons N). Hence, $\mathcal{O}(N^3)$ would mean 6.25³ (∼244) folds increase in required computing resources to obtain the enthalpies for low concentration such as 2 mol% DMA incorporation.

After a careful analysis of the correlation in our dataset, we found that the A-site cation description could be effective if one considers the ionic radius and the number of lone electron pairs as features. Certainly, this choice is not exhaustive. For instance, additional features describing the cation alignment and its dipole moment may be relevant in driving phase transitions [54] and possibly should be included. A detailed assessment of the importance of such additional descriptors will be the subject of future work.