Benchmark thermodynamic analysis of methylammonium lead iodide decomposition from first principles

Hybrid organic–inorganic perovskites (HOIPs) such as methylammonium lead iodide (MAPbI3) are promising candidates for use in photovoltaic cells and other semiconductor applications, but their limited chemical stability poses obstacles to their widespread use. Ab initio modeling of finite-temperature and pressure thermodynamic equilibria of HOIPs with their decomposition products can reveal stability limits and help develop mitigation strategies. We here use a previously published experimental temperature-pressure equilibrium to benchmark and demonstrate the applicability of the harmonic and quasiharmonic approximations, combined with a simple entropy correction for the configurational freedom of methylammonium cations in solid MAPbI3 and for several density functional approximations, to the thermodynamics of MAPbI3 decomposition. We find that these approximations, together with the dispersion-corrected hybrid density functional HSE06, yield remarkably good agreement with the experimentally assessed equilibrium between T = 326 K and T = 407 K, providing a solid foundation for future broad thermodynamic assessments of HOIP stability.


Introduction
Methylammonium (MA) lead iodide (CH 3 NH 3 PbI 3 or MAPbI 3 ) is a hybrid organic-inorganic perovskite (HOIP) semiconductor that has garnered much interest as a light-harvesting material for photovoltaic cells due to its facile synthesis and high power-conversion efficiency .As a result of this development, certified power conversion efficiency of devices based on the HOIP paradigm has grown to more than 25.5% in recent years for single-junction cells [31], reaching 26.1% in July 2023 [32].Tandem cells coupling HOIP and Si devices have rapidly increased in efficiency, reaching 33.7% certified efficiency in May 2023 [32].While exact composition data are not always immediately available for the latest devices, an earlier (2022) record device reached 31.3% efficiency in a device without MAPbI 3 and 30.9% in a device containing a small amount of MAPbI 3 , respectively [33].However, stability remains an ongoing research priority for hybrid metal halide based devices [1-3, 5-7, 11-18, 21, 22, 24, 26, 27, 29, 34, 35].The prototypical compound MAPbI 3 tends to decompose when exposed to realistic operating conditions, most notably high temperatures and humidity [1, 2, 7, 10-18, 22, 24, 27, 29, 34, 35].Several approaches have been investigated to improve the lifespans of perovskite photovoltaic cells, such as encapsulating the cells with impermeable materials [1,3,5,6,24] and altering the chemical composition of the perovskites.For instance, methylammonium ions can be replaced with formamidinium (CH 5 N + 2 ) [36][37][38] or cesium [11,[38][39][40].Indeed, demonstrated device stability in methylammonium-free perovskite cells can reach 2000 hours in the current literature [6,7,16,21,24,26], but perovskite solar cell stability nevertheless remains an ongoing issue for technological viability over commercially relevant device lifetimes.Long-term stability remains a general research topic for the much broader class of metal halide HOIPs.
One step towards developing mitigation strategies against decomposition is to assess thermodynamic stability boundaries of the broader class of HOIP materials in a given environment.Specifically, for assessing the effectiveness of encapsulation-based [1,3,5,6,24], it would be useful to be able to predict temperature-pressure equilibria of HOIPs with their decomposition products from first principles.However, MAPbI 3 and other hybrid perovskites present some difficulties with respect to computational modeling due to the strongly anharmonic nature of the organic-inorganic lattice.Specifically, the polarity and rotational asymmetry of the methylammonium ions introduce complex internal degrees of freedom to the material, causing the lattice around them to distort.[4,9,13,19,20,28,30,[41][42][43][44] Depending on the temperature, these ions may be temporarily fixed by hydrogen bonds (at cryogenic temperatures [44]), they may switch between several minimum-energy orientations or they may even dynamically tumble [41,44].Additionally, the orientational preferences of these ions causes the surrounding lattice to distort into non-cubic structures (specifically, orthorhombic and tetragonal phases) at low temperatures.One option is to apply ab-initio molecular dynamics and thermodynamic integration of anharmonic terms to the problem of MAPbI 3 thermodynamics and decomposition [28,41,45], but this method is computationally demanding, with either significant computational cost (if carried out by a direct first-principles approach) or a nontrivial workflow by machine-learning an interpolated potential energy surface with lower computational cost and sufficient precision [4,27].For simpler solids, a widely employed method to capture a first approximation of temperature/entropy components is the quasiharmonic method (QHM) [15,28,45,46], which uses the locally stable equilibrium geometry and takes lattice thermal expansion into account, but otherwise remains within the bounds of the harmonic approximation.While the quasiharmonic approximation does not account for all consequences of anharmonicity [47], some of its limitations can potentially be mitigated.As a workflow, the QHM is computationally tractable and, additionally, naturally includes approximate (i.e.harmonic) nuclear quantum effects in the free energy.The QHM has previously been successfully applied to MAPbI 3 to capture properties of the material itself (e.g.lattice parameters, [20,48] elastic moduli and phase transitions [20,28]).
We here assess the quality of several density functionals in combination with the QHM and a simple entropy correction to account for MA + cation orientational disorder (similar to a correction applied and validated previously by Saidi and Choi [20]) for a somewhat more challenging purpose, i.e. a chemical reaction equilibrium.Specifically, for MAPbI 3 , we are in the fortunate position that a well established temperature-pressure equilibrium with its decomposition products is published [14,34].Ciccioli et al [14,34] determined the product gas pressures needed to bring the system into equilibrium as a function of temperature, using Knudsen effusion mass loss (KEML) and Knudsen effusion mass spectrometry (KEMS).Three candidate reactions were considered [34]: Ciccioli et al found that the pressure calculated from pathway (2) closely matched measurements of the MAPbI 3 system.In principle, the products of pathway (1) are more energetically favorable than those of pathway (2) [12,16,34,35].However, pathway (1) is kinetically hindered at moderate temperatures (326-407 K for the KEML and KEMS experiments [14]), allowing pathway (2) to dictate the equilibrium pressures in this temperature range.Therefore, pathway (2) will be considered the dominant mode of decomposition in the present work.We stress that our study thus focuses on the experimentally established thermodynamic equilibrium (2), i.e. the starting and end points of a decomposition reaction, but we do not make any claims regarding the numerous possible, detailed kinetic processes involved in the decomposition or formation of an experimental sample.The decomposition mechanism (2) is in line with the broader literature, [6,12,27,34] and others (though without the detailed analysis of Ciccioli et al) also report typical MAPbI 3 decomposition temperatures in the range of 100 • C (373 K) [8,49].It is worth noting that the equilibrium pressures associated with the gas-phase decomposition products for the equilibrium temperature range investigated by Ciccioli et al are extremely low (approximately 10 −2 -10 −3 Pa) compared to typical laboratory gas pressures, indicating that the stability of MAPbI 3 based devices against decomposition might increase significantly even in a trace gas atmosphere of its prospective decomposition products.Our analysis indicates that the QHM, coupled with the Heyd-Scuseria-Ernzerhof (HSE06) hybrid density functional, [50,51] the Tkatchenko-Scheffler (TS) dispersion correction [52] and the aforementioned entropy correction [20] reproduces the published experimental equilibrium remarkably well.Interestingly, already the harmonic approximation for fixed lattice parameters (i.e.without the QHM correction) captures the reaction equilibrium, with only small energy corrections from volume expansion in the QHM.This observation lays the foundation to use the same, tractable approach for other hybrid metal halide perovskite materials in the future.

Methodological choices
The way to perform a computational evaluation of pathway (2) for comparison with experiment is, in principle, straightforward.One must simply calculate the Gibbs free energies g i of the reactants and products i and then, for any given temperature, find the pressure at which the energy difference between the reactants and products, ∆g, is zero.Even if the product gases are initially trapped in a solid-solid interface, they may thermodynamically be considered to be in equilibrium with an external gas phase (especially for the anticipated long lifetimes of an encapsulated photovoltaic device), and therefore it is valid to treat the material as if it is in contact with a gas.
Several key decisions must be made for a practical first-principles computation of the relevant Gibbs free energy differences and the resulting temperature-pressure equilibria: (i) Selection of the underlying method to compute total energies (the potential energy surface), i.e. the density functional approximation and relativistic formalism methods to be applied.In general, the validity and limits of different methods are not a priori clear.Providing this validation for several functionals is one of the objectives of this paper.Specifically, we apply the Perdew-Burke-Ernzerhof (PBE) functional [53], the PBE-based hybrid functional known as PBE0 [54], and the HSE06 hybrid functional [50,51].In each case, we account for the energy associated with van der Waals type dispersion interactions using the pairwise TS van der Waals correction [52].Our choice of the TS approach is motivated by a significant body of experience with the TS method applied to HOIP materials in our group, showing that the TS method generally produces equilibrium geometries in good agreement with published experimental structures [4,13,20,[55][56][57][58][59][60].The supporting information of [55] also shows that predicted structure differences for MAPbI 3 between the TS method and a more sophisticated many-body dispersion approach [61] are insignificant.Finally, a simple, theoretically motivated reparameterization [62] of the atomic polarizability of Cs renders the TS method applicable to Cs-containing HOIP nanostructures as well [63].(ii) Selection of the method used to calculate the free energy components associated with temperature and entropy of the solid and gaseous species.As mentioned above, the QHM is here applied to MAPbI 3 and PbI 2 .The free energies of the gases are calculated by the harmonic oscillator / rigid rotor model.Anharmonic corrections to the free energies of the molecular species could be incorporated by more complex workflows of their own and have been demonstrated, for example, by Piccini et al for small alkanes in zeolites [64,65].In contrast, our analysis relies initially on the QHM only.The relevant equations to assess the g i in solids and gases using the harmonic approximation or variants thereof are well known and included in the appendix to keep this paper self-contained.We then incorporate a simple configurational entropy estimate to account for the most likely orientations of the MA + cations, aiming to capture the zero-order entropy difference due to MA + orientational changes between important phases of MAPbI 3 (specifically, the low-temperature orthorhombic phase and the high-temperature cubic phase).
The free energy calculations for the cubic phase of MAPbI 3 are based on a fully-relaxed (including unit cell) orthorhombic geometry.The geometry of the MA + cation dictates that a minimal unit of MAPbI 3 cannot be locally of cubic symmetry, unlike the spatial or thermal average of such a unit in a diffraction experiment.Rather, even the 'cubic' phase of MAPbI 3 in a diffraction sense will still exhibit significant local lattice distortion and tilting of octahedra [19,28,30,42].Even at temperatures in which MAPbI 3 exists in the cubic phase, the orthorhombic unit cell may therefore provide a better approximation of the instantaneously preferred configurations of the crystal on the nanoscale.
In our calculations, we neglect the pV term of the Gibbs free energy of either MAPbI 3 or PbI 2 since the partial pressures of the involved gases remain relatively low (<50 mPa [34]) We note that, in high-pressure experiments, this term would be important.The high-pressure physics of MAPbI 3 has also been studied in other works [15,66] but does not pertain to the specific decomposition equilibrium conditions considered here.

Correction for configurational entropy of MAPbI 3
The QHM is suitable for a first-order approximation of the thermodynamics of MAPbI 3 , but it does not take into account the intrinsic disorder of the MAPbI 3 lattice.At moderate to high temperatures, the methylammonium ions within the lattice are capable of realigning or even tumbling, causing distortion of the lattice around them.This overall anharmonic motion is also intricately linked to the different phases exhibited by methylammonium lead halides at different temperatures, as, e.g.some of us have shown in other work, [44] comparing temperature-dependent high-resolution x-ray diffraction and spectroscopies to ab initio molecular dynamics derived structural expectation values in MAPbBr 3 .In general, modeling this behavior adds considerable complexity to any numerical approach; molecular dynamics simulations have been used to characterize it [20], but these are computationally intensive, and may become prohibitive if one wishes to model the behavior of MAPbI 3 over a range of temperatures.In other solid systems with internal rotational degrees of freedom (such as hydrocarbons in zeolite lattices) [64], accurate calculations of the free energies associated with these degrees of freedom can be made by mapping the rotational energy well of the free species and thus calculating its rotational partition function.For example, Piccini et al quantify the net impact of the anharmonic contribution to the adsorption Gibbs free energy ∆G of an ethane molecule in a zeolite cavity to ≈0.12 eV at T = 313 K [64].However, this method is expected to be most accurate in systems where each rotating species is functionally isolated from the others, and may therefore be treated separately from them.Due to the proximity of the methylammonium ions in the MAPbI 3 lattice and the flexible nature of the surrounding structure, such a method may not be applicable.
We here adopt a simple entropy correction strategy that accounts for expected preferential orientations of MAPbI 3 as validated earlier through explicit molecular dynamic simulations by Saidi and Choi [20].If each methylammonium ion is considered to be in one of several low-energy orientations (each with probability p i ), then one may formulate a straightforward configurational entropy term from the random alignment of ions: Molecular dynamics simulations have shown that, in the cubic phase, there are six most probable orientations of the methylammonium ion in MAPbI 3 [30,41,43,67].For six equally-probable configurations per methylammonium ion (and, thus, per formula unit), the formula equation ( 4) simplifies to the following expression: It is important to note that, by assigning a probability of 1/6 to each possible orientation of each methylammonium ion, we have made the assumption that the individual ions' orientations are uncorrelated with one another.Molecular dynamics simulations [68,69] suggest that the static correlations between the orientations of neighboring MA + ions are highly significant at low temperatures, but decrease greatly as temperature increases to 300 K and above.As such, we expect this method to systematically and significantly overestimate the configurational entropy at low temperatures, but it may be approximately accurate in the temperature range (300-400 K) considered here.At T = 313 K, equation ( 5) yields Ts config = 0.05 eV.As expected, this value is somewhat lower than the ethane/zeolite estimate of 0.12 eV by Piccini et al [64].Due to the charged and polar nature of the MA + ion, it is expected that the forces between the ion and the host lattice would be significantly stronger than those between the ethane molecule and the zeolite lattice.These forces restrict the motion of the MA + ion to the extent that a reduction of configurational entropy by a factor of ∼2 compared to the ethane/zeolite system is understandable.

Computational details
The structures of MAPbI 3 , PbI 2 , methylamine (MAe, to distinguish the neutral species CH 3 NH 2 from the MA cation CH 3 NH + 3 ) and HI were obtained by relaxing them with the FHI-aims code [70,71] using the PBE [53] generalized gradient approximation functional with the TS van der Waals correction.FHI-aims is an all-electron code based on numeric atom-centered basis functions, an efficient, high-precision recipe for electronic structure calculations verified by numerous past benchmarks [72,73] and also validated in numerous studies of perovskites [55][56][57][58][59][60].As an all-electron code, FHI-aims represents the potential near the nucleus without any shape approximations, i.e. pseudopotentials or projectors are not used.Specific numerical settings used, including k-space grids for solids, are summarized in an appendix 'Computational Settings' .For static total-energy calculations, we also obtained fully relaxed structures using the HSE06+TS and PBE0+TS hybrid functionals.In addition to the atomic coordinates, the lattice vectors of MAPbI 3 and PbI 2 were relaxed as well, except for the QHM related geometries, where fixed lattice parameters are necessary to simulate thermal expansion.Note that FHI-aims can perform both periodic and non-periodic calculations on equal footing and isolated molecules were thus treated non-periodically.Also, the algorithm used to relax structures [74] is robust even for complex, subtle HOIP crystal structures and has been used in our group to validate detailed single-crystal x-ray refined structures [59,60].The starting structure for MAPbI 3 was the computationally relaxed (PBE+TS) low-temperature orthorhombic unit cell of [55].PbI 2 is the DFT-relaxed version of the P3m1 literature structure (the 2 H polytype) [75].For DFT-PBE+TS calculations, k-space sampling was carried out using Γ-centered integration grids of 6 × 6 × 6 points for MAPbI 3 and 12 × 12 × 9 points for PbI 2 .For hybrid DFT calculations (HSE06+TS and PBE0+TS), we used k-space grids of 4 × 4 × 4 points for MAPbI 3 and 10 × 10 × 8 points for PbI 2 .The parameters for the HSE06 functional were set as suggested by Krukau et al [76] (exchange mixing 0.25, screening parameter 0.11 inverse Bohr radii).Scalar-relativistic total energy calculations were calculated using the atomic zero-order regular approximation as defined in equations ( 55) and ( 56) of [71].Experimental and computationally relaxed lattice parameters are compared in table 1 and are in close agreement with one another.The harmonic vibrational free energies of the isolated molecular species MAe and HI were calculated using finite displacements of the nuclei of 0.0025 Å, using an existing script [78] implemented in the FHI-aims package, which also provides the rotational free energy contribution.The vibrational free energies of MAPbI 3 and PbI 2 were obtained using finite-difference phonon calculations as implemented in the PHONOPY package [79], with an atomic displacement parameter of 0.01 Å. Due to the large number of individual first-principles calculations necessary to determine the free energies (especially for MAPbI 3 ) in the finite-difference based QHM approach, vibrational and phonon calculations were carried out using DFT-PBE+TS, which is computationally affordable for this task.For the HSE06+TS and PBE0+TS hybrid functionals, PBE+TS vibrational energy and entropy terms were used as fixed corrections.
In order to account for the thermal expansion of MAPbI 3 and PbI 2 , the QHM was used.As mentioned above, for MAPbI 3 , we use the molecular arrangement of the low-temperature orthorhombic phase as the most plausible short-range order pattern to apply the QHM and later apply an entropy correction term to account for configurational disorder.For the QHM, we scale the computationally relaxed orthorhombic lattice (table 1) by uniform factors, i.e. we model the volume expansion while keeping the cell fixed at the local cell shape known from the orthorhombic phase).We note that the computational orthorhombic cell's lattice parameters deviate from exactly equal lattice parameters by 2.8%.However, it is important to note that local fluctuations will always occur in the soft lattice of a halide perovskite and our simple scaling of the computationally determined unit cell is therefore a justifiable approximation of experimental reality within the confines of the QHM.At each volume probed for phonon calculations of MAPbI 3 , all cell-internal atomic positions were fully relaxed.Volumetric scaling factors of 0.99 to 1.06 were used.The case of PbI 2 is somewhat more complicated, in that it has a hexagonal crystal structure with two non-equivalent axes.Therefore, we sampled free energy landscape over a two-dimensional grid of scaling factors of the two independent lattice parameters a and c.Specifically, a was scaled by length scaling factors (LFs) between 0.98 and 1.05.The overall unit cell volume V was scaled by a separate volume scaling factor (VF) between 0.99 and 1.1, defining the respective values for the c lattice parameters associated with probed each pair of (a, V).In figure 1, the resulting grids of free-energy values are shown for three different temperatures considered for the QHM.As expected, the predicted volume and lattice spacing a increase with increasing temperature.The actual minimum free-energy values of a and V as a function of temperature were determined from a smooth bivariate spline interpolation of free-energy landscapes on (a, V) grids analogous to those displayed in figure 1.

Internal energies
Using the total energies associated with the DFT-relaxed structures (no temperature or entropy corrections), decomposition energies of MAPbI 3 were calculated straightforwardly as ∆e = e (PbI 2 ) + e (HI) + e (CH 3 NH 2 ) − e (MAPbI 3 ) . ( The resulting ∆e values for the three different density functionals are listed in table 2. Clearly, MAPbI 3 is stable against decomposition absent vibrational energy and entropy corrections, consistent with its experimental stability at low temperature.The variation in predicted reaction energies between three different DFT methods is less than 0.1 eV per formula unit.The hybrid functionals (PBE0+TS and HSE06+TS) predict a slightly higher (less favorable) internal energy change for the decomposition reaction than the PBE+TS functional does.

Gibbs free energies of solid and gaseous species
We next consider the impact of finite temperature and (for gaseous species) finite partial pressure contributions to the Gibbs free energy difference associated with reaction (2).Specifically, we consider the impact of phonons at varying unit cell volumes in the solid materials and vibrational, rotational and ideal-gas translational degrees of freedom for the molecular species in equilibrium with their gas phase, as summarized in the appendix (equations (A.1)-(A.7)).The Gibbs free energy change of reaction varies with respect to three independent variables: temperature, total gas pressure p tot , and the composition of the gas (fixed by the mole fraction of hydrogen iodide, x).Thus, in order to graph the equilibrium total pressure as a function of temperature, we must choose a value of x.While it would be sensible enough to choose x = 1 2 as a default choice, this value is not optimal for the purposes of comparing our calculated pressures to experimental KEMS and KEML results, as the gas mixture inside a Knudsen cell is not equimolar, but is governed by the Knudsen effusion condition [80,81], which is given in equation ( 7): where p i , M i , and dm i dt are the partial pressure, molecular weight, and mass effusion rate of species i, A is the area of the effusion orifice, and f is a geometric correction factor.We can reformulate equation ( 7) in terms of the particle effusion rate rather than the mass effusion rate ( dn i dt = 1 M i dm i dt ) and thus obtain Because each formula unit of MAPbI 3 decomposes into one molecule each of HI and MAe (in addition to the solid PbI 2 ), we know that the particle effusion rates of these species are equal ( dnHI dt = dnMA dt ) when the Knudsen cell has reached a steady-state condition (i.e.neither gas is accumulating in the cell).As such, we may calculate the quotient of the partial pressures according to which dictates the value of x as The inclusion of the vibrational free energy contributions introduces a zero-point energy term.The zero-point energy for each reactant and product is reported in table 3. The lattice parameters of MAPbI 3 and PbI 2 are kept fixed for this step of the analysis, i.e. quasiharmonic contributions related to thermal expansion are not yet included.The zero-point energy contributions disfavor the stability of MAPbI 3 by about 0.28 eV.
In addition to zero-point energy corrections, finite-temperature free energy terms result from the vibrational/rotational and phonon calculations.These terms cause the decomposition products to be increasingly favored over MAPbI 3 as temperature increases.Figure 2 shows the reactant and product free energies vs. temperature, at four different, exemplary product gas pressures, this time also including the quasiharmonic terms for MAPbI 3 and PbI 2 .The point at which the reactant and product free energies cross, for a given gas pressure, represents equilibrium at that pressure.The equilibrium temperature depends on the partial pressures of the product gases due to the pressure dependence of the translational free energy of the gas.Above the equilibrium temperature, MAPbI 3 is unstable and is predicted to decompose.As seen in figure 2, low product gas pressures correspond to low equilibrium temperatures, while high pressures result in high equilibrium temperatures.Figure 2 indicates that the equivalent product gas pressures associated with MAPbI 3 decomposition are indeed small (less than 0.1 Pa) around 400 K (the upper limit of the range investigated by Ciccioli et al).Even modest counterpressures of the product gases MAe and HI would increase the decomposition temperature of MAPbI 3 by a substantial amount, far outside the anticipated operating temperatures of typical MAPbI 3 based devices.
Figure 2 also compares the change in the Gibbs free energy that results from including or excluding the entropy correction equation ( 5) that accounts approximately for MA + cationic configurational disorder.
Although the difference appears small on the y-axis scale (several eV) plotted in figure 2, it does cause a shift of the prospective equilibrium temperature to measurably higher values (around ∼10 K), as discussed in more detail further below.

Impact of the density functional approximation used
In the left panel of figure 3, the equilibrium pressure of MAPbI 3 and its degradation products is plotted as a function of temperature, showing different computational approaches alongside the experimental data  5), respectively) and decomposition products vs temperature, at four selected pressures (green, red, purple, and brown curves, for p = 10 −5 Pa, 0.001 Pa, 0.1 Pa, and 10 Pa, respectively).The partial pressures of HI and MAe are chosen according to equation (10), and the reported pressure p is their sum.Free energies are defined in reference to the calculated DFT energy at the Born-Oppenheimer surface (E DFT ).All values are calculated using the DFT-PBE+TS level of theory.
reported by Ciccioli et al [14,34].This plot is the central result of the present work, since it allows us to directly compare the impact of different computational approximations on the accuracy of predicted temperature-pressure equilibria.As noted above, the DFT-PBE+TS approach including the entropy correction (blue dashed curve) yields higher equilibrium temperatures than without the entropy correction (blue solid curve), shifted by about 10 K.Although this shift is small, it is significant on the scale of the data of Ciccioli et al.The dashed blue curve (DFT-PBE+TS with entropy correction) is significantly closer to the experimental data than the solid curve without the entropy correction, but it still corresponds to equilibrium temperatures that are visibly underestimated (or, correspondingly, equilibrium partial pressures of the gaseous decomposition products HI and MAe that are visibly overestimated).
As described in the methods section, we considered the impact of two other density functionals (HSE06+TS and PBE0+TS) on the predicted temperature/pressure equilibria.Specifically, we employed the total energies of the reactants and products for structures that were fully relaxed at the target density functional, but kept the vibrational and rotational free energy corrections determined at the DFT-PBE+TS level of theory.Remarkably, either the HSE06+TS hybrid density functional or the PBE0+TS hybrid density functional produce predicted temperature/pressure equilibria that are almost precisely consistent with the experimental data, as long as the entropy correction is included.The original derivation of the PBE0 hybrid density functional does suggest a generally better performance of this functional than that of the simpler (and computationally cheaper) semilocal PBE density functional [54,82].Similarly, the HSE06 density functional and its predecessors were originally derived with the aim of achieving similar accuracy to the PBE0 functional, but at anticipated lower cost of a short-range screened hybrid density functional [50,51].It is therefore not surprising that the hybrid functionals yield somewhat improved thermodynamic predictions for the stability of MAPbI 3 , compared to the simpler PBE+TS functional.In summary, our results suggest that using total energies and static structures derived from the hybrid PBE0+TS or HSE06+TS functionals, combined with quasiharmonic vibrational and rotational entropy corrections from the more affordable PBE+TS density functional, can yield accurate stability predictions for MAPbI 3 .Given the general nature of the underlying approaches, we expect that this trend will also hold for HOIPs beyond MAPbI 3 , the case for which detailed experimental reference data are available.Equilibrium total pressure vs. temperature, for various DFT methods, including harmonic contributions with and without considering the volume expansion of the solid phases (labeled 'with QHM' and 'without QHM' , respectively).The graphs look nearly identical, illustrating the relatively minor role of the volume-dependent terms for the temperature-pressure equilibrium with the gas phase species.For the calculated and KEML results, the partial pressure ratio was determined by equation (10).For the KEMS results, the partial pressures are measured independently, so no such assumption is necessary.KEML and KEMS results were taken from [14].

Impact of the quasi-harmonic correction
A final point concerns the impact of the QHM.As mentioned above, this approximation allows one to predict temperature-dependent unit cell volumes and, thus, relaxes the restriction of free-energy calculations to a fixed unit cell.The underlying free energy is, however, still computed in the harmonic approximation, since anharmonic motion of the atoms themselves is not accounted for.For the purposes of this paper, the most important role of the QHM is the overall correction that it imparts to the free energies of the solids in equation ( 2) and, thus, to the calculated temperature-pressure equilibrium curve of MAPbI 3 decomposition.As a second result, however, one can also derive the QHM-corrected unit cell volumes themselves, albeit leaving out true anharmonic effects, such as molecular rotation.For example, Saidi and Choi [20] used the QHM for this purpose for MAPbI 3 .
Figure 4 compares different temperature-dependent experimental and computational results for the unit cell volumes of MAPbI 3 (figure 4(a)) and PbI 2 (figure 4(b)).Commenting on PbI 2 first, there is a practically constant offset between the computational and experimental data, but the trend with temperature is otherwise identical within the temperature range covered by the experimental literature data shown in figure 4. The offset of a few percent between the DFT-PBE+TS data and the experiments is within the range that is typical of experiment-computation discrepancies for unit cell volumes and the temperature trend is essentially the same.The situation for MAPbI 3 is more complex and likely indicative of the limits of the QHM for anharmonic structural degrees of freedom.In the low-temperature (orthorhombic) regime below ∼160 K, the slope of the computational and experimental volumes as a function of temperature agrees.The experimental data are only separated from the computational data by different constant offsets, as in the case of PbI 2 .The offset between our computational results and those of Saidi and Choi [20] is also constant and attributed to different density functionals used (DFT-PBE+TS here vs. DFT-PBEsol by Saidi and Choi).The trend of the computational data is otherwise nearly exactly the same over the entire temperature range plotted.In contrast, the experimental unit cell volume experiences a discontinuous change at the orthorhombic-tetragonal phase transition at ∼160 K.This change cannot be captured by the QHM and is likely associated with a freeing-up of the molecular conformation to a larger ensemble in which the molecules can rotate [44].This is the very rotation that, for the Gibbs free energy, we (as well as by Saidi and Choi) capture via the qualitative entropy term of equation ( 5), but which does not enter the QHM unit cell volume determination.As expected, volume changes above the phase transition temperature of ∼160 K can therefore not be captured in their entirety by the QHM alone.
Interestingly, and as shown in table 4, the actual QHM energy contributions (i.e.those contributions that are not already captured as harmonic terms and which are genuinely attributable to the volume expansion) only make a minute difference to the actual free energy differences of relevance in the thermodynamic   temperature-pressure equilibrium between MAPbI 3 and its decomposition products in equation (2).The overall energy scale of Gibbs free energies in figure 2 is of the order of several eV, whereas the energy changes due to the QHM in table 4 are roughly two orders of magnitude smaller.In fact, the entropy correction in equation (5) itself is larger than the overall QHM correction, making this term more important than the volume-change related energy differences.Thus, even the simple harmonic approximation at static volumes together with the configurational entropy correction would, according to our results, have captured the majority of the contribution necessary to predict the temperature-pressure stability equilibrium of figure 3.
The right panel of figure 3, which omits the effect of volume expansion, shows that the resulting changes are very small.

Conclusion
In this work, we demonstrate that applying a set of computationally affordable models to the thermodynamic stability limits of the paradigmatic hybrid organic-inorganic semiconductor MAPbI 3 can model temperature-pressure equilibrium data in close agreement with experimental reference data.Specifically, remarkably good accuracy is achieved when using (i) total energies for static structures predicted at the level of the dispersion corrected hybrid density functionals PBE0+TS or HSE06+TS, (ii) harmonic free-energy corrections for gaseous and solid species at the more affordable, dispersion corrected semilocal DFT-PBE+TS level of theory, (iii) the low-temperature local order of molecular species in the orthorhombic phase of MAPbI 3 to approximate residual short-range order, combined with a simple entropy correction to account for molecular orientational degrees of freedom and (iv) quasiharmonic corrections to account for thermal expansion effects.In contrast, just using the computationally cheaper PBE+TS approximation for (i), i.e. the internal energies, would not produce the same quality of computational results.Although still more expensive than semilocal DFT calculations such as PBE, hybrid DFT calculations can now be carried out with affordable efficiency and with high precision up to very large system sizes, e.g. up to 3383 atoms in a perovskite structure in [87].Interestingly, the quasiharmonic corrections (while important to predict qualitatively correct volume expansion trends for the low-temperature phase) do not appear to impact the overall predicted temperature-pressure equilibria of MAPbI 3 and its decomposition products strongly.We have thus validated a computational approach that is expected to be helpful to predict thermodynamic stability criteria for the much larger space of other hybrid organic-inorganic semiconductors, although we note that any necessary qualitative entropy corrections for molecular rotations will be system-dependent (e.g.depend on the number and freedom of molecular rotors present).Importantly, the underlying reaction equilibria also offer an approach to actively stabilize these materials against decomposition by offering reactants such as methylamine, HI and I 2 at sufficiently high partial pressures to establish or re-establish appropriate thermodynamic equilibrium conditions that stabilize the hybrid perovskite in operando.

Figure 1 .
Figure 1.Gibbs free energy of PbI2 at three representative temperatures at fixed volume and lattice parameter a in the harmonic approximation: (a) 100 K, (b) 300 K , and (c) 500 K, plotted as a function of the in-plane length scaling factor (LF) and the volume scaling factor (VF). Energies are aligned in the plots by subtracting the minimum energy at each temperature.

Figure 2 .
Figure 2. Gibbs free energy of MAPbI3 (orange and blue curves, without and with entropy correction equation (5), respectively) and decomposition products vs temperature, at four selected pressures (green, red, purple, and brown curves, for p = 10 −5 Pa, 0.001 Pa, 0.1 Pa, and 10 Pa, respectively).The partial pressures of HI and MAe are chosen according to equation(10), and the reported pressure p is their sum.Free energies are defined in reference to the calculated DFT energy at the Born-Oppenheimer surface (E DFT ).All values are calculated using the DFT-PBE+TS level of theory.

Figure 3 .
Figure3.Equilibrium total pressure vs. temperature, for various DFT methods, including harmonic contributions with and without considering the volume expansion of the solid phases (labeled 'with QHM' and 'without QHM' , respectively).The graphs look nearly identical, illustrating the relatively minor role of the volume-dependent terms for the temperature-pressure equilibrium with the gas phase species.For the calculated and KEML results, the partial pressure ratio was determined by equation(10).For the KEMS results, the partial pressures are measured independently, so no such assumption is necessary.KEML and KEMS results were taken from[14].

Table 4 .
QHM free energy corrections per formula unit for MAPbI3 and PbI2 at various temperatures, along with the resulting reaction energy correction.

Table 2 .
Total energy change ∆e associated with the decomposition reaction equation (2) as calculated using each of three DFT methods, per MAPbI3 formula unit (12 atoms).Note that the positive sign indicates that MAPbI3 is stable against decomposition at the DFT total energy level.

Table 3 .
Zero-point vibrational energy correction calculated for each species at the DFT-PBE+TS level of theory and using fixed lattice parameters (fully relaxed at the Born-Oppenheimer surface, i.e. not yet including any corrections due to finite temperature and/or the QHA), as well as the associated zero-point energy correction to the total reaction energy.