Phase transitions in inorganic halide perovskites from machine learning potentials

The atomic scale dynamics of halide perovskites have a direct impact not only on their thermal stability but their optoelectronic properties. Progress in machine learned potentials has only recently enabled modeling the finite temperature behavior of these material using fully atomistic methods with near first-principles accuracy. Here, we systematically analyze the impact of heating and cooling rate, simulation size, model uncertainty, and the role of the underlying exchange-correlation functional on the phase behavior of CsPbX3 with X=Cl, Br, and I, including both the perovskite and the delta-phases. We show that rates below approximately 30 K/ns and system sizes of at least a few ten thousand atoms are indicated to achieve convergence with regard to these parameters. By controlling these factors and constructing models that are specific for different exchange-correlation functionals we then show that the semi-local functionals considered in this work (SCAN, vdW-DF-cx, PBEsol, and PBE) systematically underestimate the transition temperatures separating the perovskite phases while overestimating the lattice parameters. Among the considered functionals the vdW-DF-cx functional yields the closest agreement with experiment, followed by SCAN, PBEsol, and PBE. Our work provides guidelines for the systematic analysis of dynamics and phase transitions in inorganic halide perovskites and similar systems. It also serves as a benchmark for the further development of machine-learned potentials as well as exchange-correlation functionals.


INTRODUCTION
Halide perovskites are among the most intensively studied materials of the last decade due to their attractive properties for applications in, for example, solar energy harvesting and lighting [1][2][3][4].Similar to their oxide counterparts many of these materials exhibit several different phases that are connected through soft modes and continuous or weak first-order phase transitions [5,6].This complex dynamic behavior turns out to be intimately connected to their remarkable optoelectronic properties.
The accuracy of such simulations is, however, limited by several factors, most notably (1) sampling time, (2) system size, (3) the quality of the exchange-correlation (XC) functional, and, in the case of MLPs, (4) the model uncertainty.Sampling time and system size are in particular a problem for ab-initio MD simulations, which are typically limited to time scales of a few ten picoseconds and system sizes on the order of a 1000 atoms.Previous MLP based studies were able to extend these ranges * erhart@chalmers.se to total run times of a few nanoseconds using about a thousand atoms [19,20,26].
Here, we carry out a systematic analysis of the four factors described above.We consider CsPbX 3 with X=Cl, Br, and I and the strongly constrained and appropriately normed (SCAN) [30], van-der-Waals-density functional with consistent exchange (vdW-DF-cx) [31], PBEsol [32], and PBE [33] XC functionals.The potential energy surface (PES) is mapped using third-generation neuroevolution potentials (NEPs) models and sampled using the gpumd package.The latter provides an efficient NEP implementation that enables us to routinely sample systems comprising on the order of 60 000 atoms for 100 or more nanoseconds.
We show that well converged results can be achieved using systems containing several ten thousand atoms and heating/cooling rates on the order of 30 K/ns or lower.Using bootstrapping and ensembles of models we are able to readily generate accurate NEP models with an uncertainty that is comparable or lower than the training errors.
By controlling rate and size effects as well as model errors, we are able to isolate the impact of the underlying XC functionals and thus to assess quantitatively the quality of different XC functionals for the description of phase transitions and finite temperature properties of halide perovskites.We find the vdW-DF-cx functional to perform the best among the XC functionals considered here when comparing transition temperatures and lattice constants to experimental data.
In the following section we analyze in order rate and size effects (Sect.Rate and size effects), mode uncertainty (Sect.Model uncertainty), the impact of the XC functional (Sect.Impact of XC functional and extension to other halides), and finally the transition temperature between the δ and perovskite phases (Sect.Transition to δ-phase).We then summarize and discuss the outcome of this analysis (Sect.Discussion).

Rate and size effects
The different perovskite phases are structurally closely related and connected through phase transitions with mixed continuous-first order character [34][35][36][37].For the CsPbX 3 (with X=Cl, Br, or I) materials considered in this study the perovskite lattice transforms with increasing temperature from an orthorhombic phase (P nma) via a tetragonal phase (P 4/mbm) to a cubic phase (P m 3m).Since these transitions do not involve a switch in the sign of the Glazer angles between the orthorhombic (a − a − c + ) and tetragonal phases (a 0 a 0 c + ), unlike, e.g., MAPbI 3 [19], it is possible to observe these transitions in MD simulations.Due to the remaining first-order character and the extreme heating/cooling rates that can be realized in MD simulations, one can, however, nonetheless anticipate some degree of hysteresis.
A further aggravating factor is the finite system size.For small supercells the fluctuations are naturally larger, which renders it more challenging to achieve converged results.In this section, we discriminate the effects of heating/cooling rate and system size by considering in detail MD simulations for CsPbI 3 using the full model based on the vdW-DF-cx functional (Sect.Sect.).

Rate effects
To separate rate from size effects, we first consider the former in the large size limit, using a supercell comprising 61 440 atoms, equivalent to 16x16x12 primitive orthorhombic (20-atom) unit cells.
On heating all simulations yield the correct (experimentally observed) sequence of phases irrespective of heating rate.On cooling, this sequence is reversed again regardless of rate.At the cubic-to-tetragonal transition, for a small number of simulations, one can, however, observe the simultaneous formation of multiple tetragonal domains with incompatible orientations, which can lead to the formation of domain boundaries.Since the moving of these boundaries involves a nucleation-and-growth mechanism, they remain in the simulated structure upon cooling.
The temperatures for the transitions between the perovskite phases can be readily obtained from the lattice parameters (Fig. 1a,b), revealing a strong dependence on the heating/cooling rate.For the rates below approximately 30 K/ns, the hysteresis between heating and cooling runs is 15 K are less and no longer vary systematically with the rate (Fig. 1c).By contrast, for the largest rate of 600 K/ns considered here, which is only slightly larger than values used in some earlier MLP studies [19,26], one observes a hysteresis of 100 K or more for both the 300 400 500 Temperature (K) lower and higher temperature transitions.In addition, the transition itself is smeared out in temperature, which is particular apparent for the orthorhombic-to-tetragonal transition (Fig. 1a).We observed similar trends also for the other materials and models studied in this work.We therefore conclude that rates below approximately 30 K/ns are recommended in order to achieve convergence of the transition temperatures for this class of materials.

Size effects
Next we examine the impact of supercell size on the temperature dependence of the lattice parameters and the transition temperatures.First, simulations were carried out at a heating/cooling rate of 6 K/ns using structures comprising 1280 atoms (4x4x4 primitive orthorhombic unit cells), 7680 atoms (8x8x6), 23 040 atoms (12x12x8) or 61 440 atoms (16x16x12).For the smallest supercell (7680 atoms) one notices a marked deviation from the (reference) lattice parameter parameter data from the largest supercell (61 440 atoms) around the cubic-tetragonal phase transition (Fig. 2a,b).This deviation is, however, absent for the next larger structure of 7680 atoms and at this size the transition temperatures are already converged to within 10 K of the results for the largest supercell (Fig. 2c).
A key characteristic of a phase transition that is at least partly continuous is a peak or kink in the heat capacity.The heat capacity can be readily extracted from the fluctuations of the energy in MD simulations.For the present purpose this is, however, impractical as the transition is very sharp in temperature and the temperature range of interest is wide.Here, we therefore compute the heat capacity instead by numerically differentiating the potential energy from heating/cooling runs, i.e., C p = dH/dT ≈ ∆H/∆T .This requires averaging of the data in order to obtain numerically well behaved results.To this end, we first apply a Hamming window of 0.6 K to the energy-vs-temperature data.The resulting data is numerically differentiated, after which the data is smoothened again using a Hamming window of 6 K.The Hamming window sizes are chosen to be sufficiently large to remove noise and small enough to avoid removing features (Fig. S4).
For the largest system size (61 440 atoms) the phase transitions are clearly visible as peaks in the temperature dependence of the heat capacity (Fig. 3).These features become, however, less distinct with decreasing system size as fluctuations increase.
The analysis presented in this section suggests that supercells with at least about 10 000 atoms can be expected to yield accurate lattice parameters and transition temperatures within about 10 K of the converged results.Extracting the temperature dependence of the heat capacity requires larger systems still.Even for the largest systems considered here (61 440 atoms) the noise level is rather high, but the data still allows one to accurately extract phase transition temperatures from the heat capacity data.

Model uncertainty
Having established heating/cooling rates and system sizes that yield converged results for lattice parameters and transition temperatures, we can now address the model uncertainty.To this end, we resort to ensemble models.The latter comprise five separate models constructed using five distinct 90-10 splits of the training data (Sect.Model construction).All simulations in this section were carried out using supercells with 61 440 atoms and a heating/cooling rate of 6 K/ns.Once again we use CsPbI 3 and models trained using reference data generated by the vdW-DF-cx functional as a representative example.The uncertainty of the model predictions can be estimated by the considering the standard deviation over the model ensemble (Fig. 4a,b).At 100 K this approach yields an uncertainty of up to 0.02 Å depending on direction.This value diminishes, however, quickly with temperature to a level of 0.003 Å in the tetragonal and even less than 0.002 Å in the cubic phases (Fig. S5).
Away from the phase transitions the heat capacity curves from the different models agree well with each other (Fig. 4).The actual transition temperatures, corresponding to the position of the peaks in the heat capacity curves, obtained from the full model (and the model ensemble) are 340 K (341(9) K) and 487 K (496(7) K) on heating and 331 K (323(17) K) and 470 K (478(5) K) on cooling.The agreement between the results obtained using the full model and the model ensemble support the good convergence of the models with respect to training data.The remaining hysteresis can be attributed to the mixed first/continuous-order character, which is evident from the very small but non-zero latent heat associated with these transitions [34][35][36][37].
We can thus conservatively estimate the error in the transition temperatures due to model uncertainty to be on the order of 20 K.In combination with the model performance measures (Sect.Sect. ) and the very good agreement with the density functional theory (DFT) reference data, this provides strong evidence that the NEP models constructed are accurate representations of the DFT potential energy landscape in the regions of configuration space included here.They can thus be used to analyze the performance of different XC functionals with Transition temperatures from MD simulations in comparison with experiment.Data were obtained using heating (upward triangles) and cooling rates (downward triangles) of 6 ns for supercells comprising 61 440 atoms.Colored and gray symbols indicate results obtained using full models and model ensembles, respectively.In the latter case the uncertainty calculated as the standard deviation over the ensemble is indicated by vertical bars.Experimental transition temperatures taken from Refs. 9, 38, and 39 are shown by horizontal dashed lines.
respect to the finite temperature behavior of CsPbI 3 and the other materials considered here (Sect.Impact of XC functional and extension to other halides).

Impact of XC functional and extension to other halides
We can now apply the framework established above to predict transition temperatures for CsPbI 3 using different XC functionals for generating reference data, and compare the results to experimental data.To this end, we use the same computational settings in terms of system size and heating/cooling rate as in the analysis of the model uncertainty.
The results show that all XC functionals considered here systematically underestimate the transition temperatures for both the orthorhombic-tetragonal and tetragonal-cubic transitions by as much as 200 K (Fig. 5).The agreement improves in the sequence PBE, PBEsol, SCAN, and vdW-DF-cx, which mirrors the trends observed previously for transition metals [40].
Conversely for the lattice parameters one observes a systematic overestimation relative to experiment (Fig. 6) with a similar trend.The closest agreement is obtained using the vdW-DF-cx functional, for which the lattice parameters of the cubic phase are still overestimated by 0.02 to 0.04 Å in the high-temperature limit.
Finally, we extend our investigation to CsPbBr 3 and CsPbCl 3 including the SCAN and vdW-DF-cx functionals.As for CsPbI 3 the transition temperatures (Fig. 5b,c) and lattice parameters (Fig. 6b,c) are systematically underestimated and overestimated, respectively.The ther- mal expansion is, however, captured well by all functionals.
We conclude that while the XC functionals considered here yield the correct sequence of phases, none of them are quantitatively in satisfactory agreement with experiment, systematically underestimating the transition temperatures (Fig. 5) and overestimating the lattice parameters (Fig. 6).

Transition to δ-phase
It is experimentally well established that the perovskite phases of CsPbI 3 are only metastable at low temperatures as the actual ground state structure of the material is the so-called δ-phase [41].We therefore also computed the transition temperature between the cubic perovskite and the δ-phase by free energy integration (see Fig. S6 for the free energy curves).
The transition temperature obtained by the vdW-DFcx model of 635 K with an estimated uncertainty of 20 to 40 K is in rather good agreement with experimental values of around 600 K (Table I).The SCAN and PBEsol models yield lower values that are considerably smaller than the experimental data.According to DFT (and NEP) calculations the δ-phase is the most stable structure for CsPbBr 3 .The energy difference (Table S3) is, however, much smaller than for CsPbI 3 (Table S4), leading to lower transition temperatures that are predicted to be 310 K and 151 K according to the vdW-DF-cx and SCAN models, respectively.We are unaware of experimental measurements of the transition temperature, which given the present predictions might actually be close to or below room temperature and thus difficult to observe.
For CsPbCl 3 the DFT calculations yield energy differences between perovskite and δ-phases close to zero (Table S2), suggesting that the δ-phase is actually not the most stable phase under any conditions or only at extremely low temperatures.

DISCUSSION
We have systematically analyzed four key sources of error in atomic scale simulations of phase transitions of inorganic halide perovskites, related respectively to (1) sampling time, (2) system size, (3) model uncertainty, and (4) the underlying XC functional.Based on these results, it is recommended to use heating/cooling rates of at most approximately 30 K/ns but preferably even lower and system sizes comprising at least a couple ten thousand atoms, corresponding to a few thousand primitive unit cells.We expect that these guidelines are not limited to inorganic halide perovskites but should also be heeded when modeling the dynamics of other perovskites and related systems.
The model uncertainty was assessed using ensembles of models, from which uncertainty estimates for, e.g., transition temperatures and lattice parameters can be estimated.The results show that by careful model construction the model uncertainty can be reduced to a level that allows one to quantitatively discriminate the performance of different XC functionals.
Using this approach, we found that the semi-local XC functionals considered here systematically underestimate the transition temperatures and overestimate the lattice parameters at finite temperature compared to experiment.The best overall agreement is obtained for the vdW-DF-cx functional, which also outperforms the SCAN functional.In pioneering work on the use of MLPs for probing the dynamics of halide perovskite, the latter functional had been suggested as achieving a good match with experiment [19,26,45,46].Our analysis suggests that the good agreement was likely fortuitous and likely the result of using high rates and small system sizes (see Fig. S6 for a more detailed comparison).
This raises the question how, for example, hybrid functionals or the random phase approximation would perform with regard to the finite temperature properties and dynamics of these materials [13,27,47].In the present work we included a relatively large set of structures and supercells that would pose a considerable computational challenge for either one of these methods.As noted above (Sect.Model construction), one can, however, expect that the number of training structures can be considerably reduced without a notable loss in model accuracy by using active learning.This strategy could be combined with principal component analysis to identify regions with very dense sampling [48] or entropy maximization [49] to reduce the training set size even further, eventually allowing one to build NEP or other MLP models that can represent such more accurate electronic structure methods.
Finally, we note that the accuracy of the models presented here in combination with the very high computational efficiency provide by the implementation on graphical processing units (GPUs), now enables one to sample the dynamics of these and related materials with unprecedented time resolution and accuracy.

Neuroevolution potentials
Here, we use the third generation of the NEP scheme (NEP3) [48] to build MLPs for CsPbX 3 with X=Cl, Br, and I.The NEP format employs a simple multi-layer perceptron neural network architecture with a single hidden layer [50].In NEP3 the radial part of the atomic environment descriptor is constructed from linear combinations of Chebyshev basis functions while the three-body angular part is similarly build from Legendre polynomials.Four and five-body terms of the atomic cluster expansion form [51] can be included as well but here we limit ourselves to two and three-body terms.
For the present purpose it is crucial that the NEP scheme is not only accurate but has been implemented on GPUs in the gpumd package [48].For the models described in the following this allows us to achieve a speed of 2 × 10 7 atom step/s on an NVidia A100 card, i.e., we can simulate a system of 60 000 atoms for about 150 ns per day using a time step of 5 fs.

Computational parameters
In this study we used the same hyperparameters for all models, which were chosen based on experience and pre-trials [48].The cutoffs for two and three-body interactions are 8 Å and 4 Å, respectively.There are 8 radial and 6 angular descriptor components, 8 basis functions for building both the radial and angular descriptor functions, and the angular components are expanded up to fourth order.The hidden layer contains 50 neurons.
The weights for energies, forces, and virials in the loss function were set to 1, 5, and 0.2 in gpumd units, respectively, while the weights for the 1 and 2 regularization terms were set to 0.1.The neuroevolution strategy [52] used for optimizing the parameters used a population size of 50 and was run for 200 000 generations.

Model construction
To construct NEP models we employed a bootstrapping strategy.First we identified potentially relevant phases.This included the cubic perovskite structure (P m 3m, Glazer notation a 0 a 0 a 0 ), two tetragonal structures (I4/mcm → a 0 a 0 c − , P 4/mbm → a 0 a 0 c + ), representing out-of-phase and in-phase tilts relative to the caxis, respectively, one orthorhombic structure (P nma → a − a − c + ) as well as the so-called delta-phase (P nma), which is experimentally known to be the most stable structure at least for CsPbI 3 and CsPbBr 3 .
We then calculated energy-volume curves for these five prototype structures using DFT calculations (Sect.Reference calculations) allowing both the ionic coordinates and the cell shape to relax under the constraint of constant volume until the maximum force on any atom fell below 30 meV/ Å.
Subsequently we generated supercells for each prototype with random atomic displacements using the Monte Carlo rattle procedure from the hiphive package [53] with a standard deviation of 0.04 Å.The supercell size was chosen to be between 120 and 160 atoms and the volume was varied between 85% and 110% of the respective equilibrium volume with five structures per volume and prototype.
Using these data we generated a first iteration of NEP models using the gpumd package for the optimization [48] and the calorine package for data preparation and analysis [54].One model was generated using the full data set ("full model") and five additional models ("model ensemble") were generated by using five different 90-10 splits of the available data.Using the full model we generated new structures for each prototype by running short MD simulations at pressures between −1 and 10 GPa using a temperature ramp from 20 to 620 K over 3 ns.From each trajectory we selected 12 configurations.For each of these configurations we then computed the standard deviation of the energy and forces using the model ensemble.The standard deviation over the ensemble predictions provided a measure for the uncertainty of the current model generation for the respective conditions (temperature, pressure, structure).We then computed energy and forces for the new structures using DFT calculations, added these to the training set and repeated the procedure.Typically after four generations we found that the uncertainty in the energy and forces was comparable are smaller than the respective training error indicating convergence of the model construction.We note that in principle one could have adapted an active learning strategy based on the model ensemble and only included configurations with high uncertainty as additional reference structures.Here, we decided to include rather more data in the training set but we expect that the number of structures can be reduced considerably without a notable decrease in model performance.
The final models yield RMSE scores for training and validation sets of about 2 meV/atom, 50 meV/ Å, and 15 meV/atom or better for energies, forces, and virials, respectively (Table II and Table S1).Importantly the models closely reproduce the energy differences and energy-volume curves of all the structures of interest in the present study (Fig. S1; Fig. S2; Table S2 to Table S5).The final models were subsequently used in large scale MD simulations to predict, for example, transition temperatures or lattice parameters (Sect.MD simulations).

MD simulations
All MD simulations were carried out using the gpumd code.Temperature and pressure were controlled using stochastic velocity [55] and cell rescaling [56] and the time step was 5 ps, where all simulations were run at zero pressure.
For studying the convergence with size (Sect.Size effects), we considered system sizes between 1280 and 61 440 atoms, equivalent to 4x4x4 to 16x16x12 primitive orthorhombic perovskite (20-atom) unit cells.To analyze the impact of heating and cooling rates (Sect.Rate effects) the temperature was linearly varied between 20 K and 620 K over 1 ns to 100 ns.
The production runs used to quantify model uncertainty (Sect.Model uncertainty) and the impact of the XC functional (Sect.Impact of XC functional and extension to other halides) were carried out using supercells comprising 16x16x12 primitive orthorhombic unit cells (61 440 atoms).The total simulation time was set to 100 ns and the temperature was varied over a range of 400 to 600 K corresponding to a heating/cooling rate of 4 to 6 K/ns.

Free energy calculations
For CsPbI 3 and possibly CsPbBr 3 the perovskite phases are only metastable at lower temperatures.Provided sufficient kinetic activation, below a certain temperature the perovskite structure transforms into the socalled δ-phase via a first order transition.To determine the transition temperature from the NEP models we calculated the free energies of the δ and cubic perovskite phases through thermodynamic integration using the classical method by Frenkel and Ladd [57][58][59], as implemented in ase [60].In these calculations, we used an Einstein solid as reference system, for which the free energy can be computed analytically, and used supercells containing about 1500 atoms for each phase.For each temperature the integration was carried out over 50 ps and the results were averaged over ten independent runs.

Reference calculations
DFT calculations were performed using the projector augmented-wave method [61] as implemented in the Vienna ab-initio simulation package [62,63].The exchangecorrelation contribution was represented using the vdW-DF-cx method [31], the SCAN density functional [30], the PBEsol functional [32], and the PBE functional [33].The Brillouin zone was sampled with a Γ-centered grid with a k-point density of 0.18/ Å and Gaussian smearing with a width of 0.1 eV.For the calculation of the forces a finer support grid was employed to improve their numerical accuracy.

FIG. 2 .
FIG. 2. Convergence with supercell size.Lattice parameters for CsPbI3 as a function of temperature from MD simulations during (a) heating and (b) cooling at a rate of R = 6 K/ns for different supercells.The transition temperatures extracted from these data are indicated by diamonds.A Hamming window of 0.6 K was applied.(c) Transition temperatures as a function of supercell size.All results were obtained using the full model (Sect.Model construction) based on the vdW-DF-cx exchange-correlation functional.

FIG. 3 .
FIG. 3.Convergence of heat capacity with supercell size.Isobaric heat capacity of CsPbI3 as a function of temperature for different supercell sizes obtained through numerical differentiation as detailed in the text.The phase transitions are visible as peaks in the heat capacity.All results were obtained using the full model (Sect.Model construction) based on the vdW-DF-cx exchange-correlation functional.

FIG. 4 .
FIG. 4.Assessing model uncertainty through ensemble of models.(a,b) Lattice parameters and (c,d) heat capacity as a function of temperature during (a,c) heating and (b,d) cooling from full (blue lines) and ensemble models (orange lines).All results are for CsPbI3 using models based on the vdW-DF-cx exchange-correlation functional.
FIG. 5.Transition temperatures from MD simulations in comparison with experiment.Data were obtained using heating (upward triangles) and cooling rates (downward triangles) of 6 ns for supercells comprising 61 440 atoms.Colored and gray symbols indicate results obtained using full models and model ensembles, respectively.In the latter case the uncertainty calculated as the standard deviation over the ensemble is indicated by vertical bars.Experimental transition temperatures taken from Refs. 9, 38, and 39 are shown by horizontal dashed lines.

3 FIG. 6 .
FIG.6.Lattice parameters as a function of temperature from simulation in comparison with experiment.Data were obtained using "full" NEP models trained using different XC functionals in comparison.Experimental data from Refs. 9, 38, and 39.

TABLE I .
Temperatures for the transition between the δphase and the cubic perovskite phase.Experimental data from Refs.42-44.

TABLE II .
RMSE scores for the final NEP models obtained by training against the full data set.Additional performance measures including RMSE and Pearson correlation coefficients for model ensembles can be found in TableS1.