Structure and binding in halide perovskites: Analysis of static and dynamic effects from dispersion-corrected density functional theory

In recent years, the hope of cheap, solution-based solar cells with efficiencies above 20% has fueled the scientific research in halide perovskites (HaPs).^1–7 HaPs exhibit many interesting optoelectronic properties that make them ideal candidates for future semiconductor based applications, yet some of their structural and dynamic properties remain of central scientific interest since they also present major stumbling blocks in the development of stable devices. Fundamentally, HaPs can crystallize in the typical perovskite structure ABX₃, where the B-site is a heavy metallic element such as lead, and halide ions (X) form octahedra centered around the B-site.⁸ In hybrid organic-inorganic HaPs, the voids between the inorganic scaffolds are occupied by the organic A-site [in this work methylammonium (MA)]. As is typical for perovskites, HaPs undergo phase transitions with changes in temperature. The prototypical variant MAPbI₃ has an orthorhombic symmetry up to 162 K, where it changes to a tetragonal structure. The second phase transition occurs at 327 K, where it switches to a cubic crystal structure.^9,10 Importantly, HaPs are, compared to other well established semiconductors¹¹ and perovskite materials,^12–17 mechanically soft, which could become problematic when they are exposed to the operating conditions of commercially used solar cells. Furthermore, experimental and theoretical studies have also shown that the structure and binding in HaPs can be modified already by applying relatively low pressure.^18–20 Therefore, it is interesting to study the interrelation of binding and structural modifications in HaPs, which can be induced by a change in temperature or applied pressure, and may have important consequences for their optoelectronic properties.

Such insight can be generated from first-principles based computations. However, from a theoretical point of view, there are a number of physical effects and properties, all of which could play a role regarding the structure of and binding in HaPs. First, hybrid HaPs consist of highly polarizable atoms such as iodine and lead and contain hydrogen atoms bound to electronegative species (e.g., nitrogen); hence, from simple chemical arguments, dispersive interactions, such as van-der-Waals (vdW) and hydrogen bonding, can be expected to play an important role in HaPs. Several computational studies have shown that dispersive interactions are crucial in static calculations of the lattice constants and mechanical properties of HaPs.^{12,16,21–27} However, as to how their treatment for HaPs should best be implemented in density functional theory (DFT), the most widely used method for HaPs, is still an open question.^28,29 In particular, a recent study has concluded that correcting for dispersive interactions in calculations of HaPs does not improve their description.³⁰ Second, Kohn-Sham DFT functionals, such as the frequently used generalized gradient approximation (GGA), cannot accurately describe the band gap of semiconductors even in principle,^31,32 which can be improved by applying computationally more costly hybrid functionals in the generalized Kohn-Sham framework.³³ For classical inorganic semiconductors, screened hybrid functionals were shown to improve also the description of lattice constants compared to GGA-based approaches.^34–36 The question that naturally arises then is whether a similar improvement is found when using a hybrid functional for calculating the structure of HaPs, as has been suggested recently.³⁰ Third, it is well established that spin-orbit coupling (SOC) due to the presence of the lead atom in MAPbI₃ leads to strong modifications in the electronic structures, i.e., it lifts the degeneracy of the conduction and valence band and lowers the band gap.^37,38 Finally, since HaPs are mechanically soft, they exhibit large structural dynamical effects, such as massive ionic displacements and molecular rotation around room temperature,^39,40 which is the relevant scenario for device applications. Therefore, it is important to investigate the consequences of these unusual structural dynamical effects on the pertinent structural properties of HaPs. To the best of our knowledge, the impact of finite temperature on the lattice constants and mechanical properties, and the role of dispersive interactions in calculating them, has not been addressed by means of DFT-based molecular dynamics (MD).

In this study, we first investigate how the choice of vdW correction and DFT functional affects the results of calculations for structure and binding in the prototype MAPbI₃. To this end, we compare data obtained in static DFT calculations for a primitive cubic unit cell to experimental ones, using various computational approaches. The most promising choice of method is then used to explore how different MA orientations impact structural properties in static and finite-temperature MD calculations of MAPbI₃.

A satisfactory treatment of dispersion forces within conventional approximate Kohn-Sham DFT functionals requires dispersive correction schemes. In our calculations, we considered the Tkatchenko-Scheffler method with regular Hirshfeld partitioning (TS)⁴¹ and with iterative Hirshfeld partitioning (HI)^42,43 as well as the many-body dispersion (MBD) method.^44–46 The TS and HI methods include pairwise dispersive interactions between atoms in the crystal. The HI scheme expands the Hirshfeld partitioning with an iterative process that results in a more realistic allocation of charges for strongly ionic systems.⁴² To include the screening of dispersive interactions that occurs in a many-body system, we apply the MBD method based on the random phase expression of the correlation energy, as developed by Tkatchenko et al.⁴⁴ 

For the comparison of results obtained by using a GGA functional and a hybrid functional, we apply two very common variants in the context of solid-state calculations, namely, the PBE⁴⁷ and HSE functionals.^48,49 In the HSE functional, the short-range exchange energy is calculated as a mix of PBE and exact exchange,^48,49 which was often found to improve the description of electronic-structure and structural properties of semiconductors.^34–36,50 SOC describes the interaction of the spin angular momentum with the orbital momentum, and it plays a crucial role in systems with heavy elements such as lead. In our calculations, it was described fully self-consistently within the framework of non-collinear magnetism,⁵¹ in conjunction with the various DFT methods applied here.

In order to obtain structural and mechanical parameters, we use an equation of state that accurately describes the macroscopic properties of the crystal. One of the most frequently used ones is the Birch-Murnaghan equation of state (BMEOS).^52–54 It describes the free energy of a solid as a function of the volume of the unit cell, E(V), at isothermal conditions as

E₀ is the free energy at the equilibrium volume, V₀, and B₀ and $B0′$ are the bulk modulus and its derivative, respectively. In order to obtain these parameters, we first calculate the free energy of the static system at eight equidistant volumes between 90% and 111% of the experimental value¹⁰ and fit Eq. (1) using these eight values (tests with more data points did not show any significant improvement) with the method of least squares. Note that while the E₀ and $B0′$ are required to fit Eq. (1), they are not relevant for our discussion.

In our calculations, we first considered the cubic primitive unit cell of MAPbI₃ as has been reported experimentally; see Fig. 1(a).¹⁰ To be able to also investigate the effect of MA orientation, we used a supercell that consists of eight unit-cells stacked together in a 2 × 2 × 2 pattern [see Fig. 1(a)], with the MA ions orientated in parallel or anti-parallel to each other. In this way, we can test the effect of the assumption that MA molecules are oriented perfectly in parallel throughout the material on computing the structural and binding properties of MAPbI₃. Note that in these calculations, the angle between the MA molecules was fixed during the atomic relaxation. To fit Eq. (1) with the data obtained in the MD calculations, the average free energy along a trajectory of 20 ps was calculated for each volume point, which was used, together with the standard deviation as the uncertainty, in the fit of Eq. (1). We note that since the fit errors of B′ were quite large in the case of the MD calculations, in fitting Eq. (1) we set B′ to the value obtained in the static calculation of the unit cell.

The DFT calculations were performed with a plane-wave basis (cutoff energy: 400 eV) and the projector-augmented wave (PAW) method,⁵⁷ as implemented in the VASP code.⁵⁸ For unit cell calculations, a 6 × 6 × 6 grid of k-points centered around the Γ-point was used, and for static and dynamic calculations that considered the supercell, the grid was reduced to a sampling of 3 × 3 × 3. For each static calculation, the ionic coordinates were relaxed, using the PBE functional and the respective dispersive correction scheme, with the Gadget tool in internal coordinates,⁵⁹ before calculating the energies to fit Eq. (1). In the calculations applying the HSE functional or SOC, the respective PBE-based geometry was used to reduce computational costs. For all static calculations, the geometry was considered relaxed when the forces acting on the atoms were below 10 meV/Å. In the canonical (NVT) MD simulations, a MAPbI₃ 2 × 2 × 2 supercell with randomly orientated MA was used as a starting point for the structural geometry, in line with the recommendations provided by Lahnsteiner et al.⁶⁰ It was then equilibrated for 5 ps and computed for 20 ps in time steps of 1 fs at a temperature of 400 K. This protocol is sufficient to compute the impact of the dynamic nuclear effects on the structure and binding despite the limited supercell size⁶⁰ and trajectory length. The latter was checked explicitly by adding another 10 ps to the simulation, which resulted in only insignificant changes in the calculated observables. Schematic representations of the structures were generated with VESTA.⁶¹ 

Figure 2 shows the energy change due to volume change, i.e., ΔE(V) = E(V) − E₀, for the static calculations of the MAPbI₃ unit cell, calculated by using the PBE and HSE functionals augmented by different dispersion corrections. The optimized volume (V₀) and bulk modulus (B₀), obtained from the fit, are listed in Table I. Considering the PBE data, it can readily be seen that using dispersive corrections has a large impact on the result. Furthermore, it was found that the results from the TS and MBD methods are very similar and provide a V₀ of 256.2 ų (PBE + TS) and 257.1 ų (PBE + MBD) which are close to the experimental range of 247-253 ų.^9,10,55 By contrast, the bare PBE calculations are far off (272.9 ų) the experimental range, in agreement with literature findings,^12,21–27 and the PBE + HI (262.0 ų) results lie in between the results of bare PBE and PBE + TS/PBE + MBD. The calculated B₀ data show similar trends, i.e., the experimentally reported values (12-16 GPa)^14,56 agree well with the PBE + TS, PBE + MBD, and PBE + HI results but deviate more from the bare PBE value. From these data, one can see that the TS dispersive correction scheme with regular Hirshfeld partitioning and the MBD approach perform best in reproducing structural and binding properties measured in experiments. It should be noted that the experimental data were recorded at elevated temperatures, at least at T ≈ 330 K, while the calculations were performed at 0 K and do not address thermal expansion, an effect we discuss below based on our results obtained from finite-temperature MD calculations.

Figure 2 also shows the results for regular and TS-corrected calculations using the HSE functional (see Table I for full dataset using the other correction schemes). While there is some visible difference between the PBE and HSE calculations without dispersive corrections, the bare HSE value for V₀ is still quite off the experimental result, in stark contrast to findings reported for conventional inorganic semiconductors.^34–36 Considering the dispersion-corrected HSE data, it is found that these compare almost equally well to experiment as their PBE counterparts. Importantly, the optimized V₀ value using the PBE + TS, PBE + MBD, HSE + TS, and HSE + MBD is all within 5 ų, which is smaller than the range of experimentally reported values. Hence, we find that the improvement because of using the HSE functional is minor compared to using dispersive corrections, which is thus found to be the essential computational ingredient for obtaining accurate structural and binding properties of MAPbI₃. From these findings, we argue that using PBE + TS provides a reasonable choice for calculating structural and binding properties of HaPs such as MAPbI₃. For completeness, we note that calculations including SOC showed very little difference to those without, as can be seen in Table I.

It is well-known that the electronic interaction between MA and the inorganic ions in MAPbI₃ is minor since the electronic states of MA are energetically not close to the band edges. Nevertheless, electrostatic and dispersive interactions between the inorganic ions and MA could still be important for the structure and binding of MAPbI₃. In order to test this, we varied the MA orientation in a 2 × 2 × 2 supercell, considering the extreme cases of either perfectly parallel or antiparallel MA molecules. It is noted that the former scenario is equivalent to considering a primitive cubic unit cell of MAPbI₃ containing one MA unit. When calculating the structural parameters of the 2 × 2 × 2 supercell with different MA orientations, we limit the calculations to the bare PBE and the PBE + TS approach since neither using HSE nor including SOC has shown a significant improvement that would justify the increase in the computational cost (see above and Table I). Furthermore, using the TS and MBD dispersive correction scheme provided essentially equal results, but the computational costs associated with the MBD method increase more rapidly when the number of atoms becomes larger compared to the TS scheme.

The first relevant observation in the data shown in Table II is that the results of the supercell calculations considering the parallel MA orientation are, within the fitting error, identical to the results obtained with the calculations for the primitive unit cell, as expected. However, a noticeable difference occurs in both the structure (see Fig. 3) and the optimized parameters (see Table II) when the MA orientation is changed to the other extreme of perfectly antiparallel MA molecules. In the case of the parallel MA orientation, at all volumes, the octahedra retain the cubic symmetry, whereas in the antiparallel orientation, they tilt strongly (see Fig. 3). Indeed, this effect is important since for the fit of the supercell data calculated with bare PBE, we find that V₀/8 changes from 272.7 ų, for the parallel orientation, to 265.8 ų, for the antiparallel one. The same trend is confirmed in the PBE + TS calculations although the differences are smaller. Hence, the interaction between the inorganic cage and the MA molecules is important for the structure and binding in MAPbI₃. Since our DFT calculations show that the antiparallel orientation is preferred in MAPbI₃, i.e., the free energy is consistently lower for the antiparallel MA orientation at all eight volumes, this requires further investigations.

Hence, motivated by the finding that MA orientation impacts structure and binding in MAPbI₃, we performed fully unconstrained NVT MD calculations at 400 K. This is relevant, since the MA unit was found to be only weakly bound in the cubic phase,^62,63 and hence undergoes rotational and translational motion at elevated temperatures, which could impact the energetics in the material dynamically. Furthermore, the MD calculations allow for testing whether the effect of dispersive corrections seen for the static structures is still visible at elevated temperatures, a comparison that has not been attempted previously but is important for investigating dynamical effects in the structural interaction of MAPbI₃. Figure 4 shows the entire distribution as well as the mean value of the PBE- and PBE + TS-calculated change in free energy of MAPbI₃, ΔE, determined for a 20 ps MD simulation, again fitted by Eq. (1). The most apparent finding from Fig. 4 is yet again the sizable difference between the PBE and PBE + TS curves, as also quantified by the V₀ parameters provided in Table III. Furthermore, a slight decrease in V₀ is found in the MD calculations at 400 K compared to the static 0 K calculation, contrary to the expected thermal expansion. In regard to B₀, while in the case of PBE, it is similar to the 0 K calculation of the primitive unit cell, for PBE + TS, it is slightly lower than what was obtained in the static calculations.

In order to better understand the origin of these findings, the average atomic positions along the PBE + TS MD trajectories were calculated. Figure 5 shows that depending on the volume of the supercell, two different types of structures emerged. For the three smallest volumes considered in the MD, the octahedra were found to be tilted and the C-N bond of MA is still clearly visible in the average structure: successive MA molecules are oriented in parallel or antiparallel to each other along one direction and orthogonal to each other in the other two directions, aligning with the long axis of the rhombus created by the tilted octahedra. For the five larger considered volumes, on the other hand, the octahedra form an on average near perfect cubic symmetry, and the average carbon and nitrogen atoms are almost conjoined. The latter could either mean that there is absolutely no preferred direction of the MA molecules, i.e., it rotates such that it is entirely disordered over the course of the 20 ps trajectory, or that there are preferred directions which are equally occupied thermally.

The results from the static unit cell calculations confirmed previous findings that showed the importance of including dispersive interaction in DFT calculations for the structure and binding in MAPbI₃.^12,21–27 Here, we went beyond in testing a range of dispersive-correction schemes: PBE calculations corrected by the TS scheme with the regular Hirshfeld partitioning and the MBD scheme showed good agreement with experimental data. The TS method with iterative Hirshfeld partitioning performed slightly worse than the regular TS and MBD schemes. While the iterative Hirshfeld partitioning has indeed been shown to improve the description of ionic materials,^42,43 MAPbI₃ is a more complex case since it is a hybrid organic-inorganic system that contains covalent bonds (in the MA molecule) as well as partially covalent-ionic bonds (in the inorganic framework). One perhaps surprising result is that the MBD method did not improve the results significantly compared to the TS method even though one would expect that the iodine atoms, which account for most of the vdW-interactions,²² are screened by the surrounding dielectric environment. We considered a comparison of the C₆ parameters calculated from PBE + TS and PBE + MBD, which determine the dispersive correction in either case,^28,29 to further understand this result. We find that the changes are minor, on the order of 1% or smaller, and, furthermore, do not depend on the volume of the unit cell. This implies that the screening is barely affected by the changes in distance in the calculations at different volumes and thereby only leads to a constant energy shift. Furthermore, we considered the higher-order contributions to the dispersive energy calculated with PBE + MBD, to find that indeed the second-order term is dominating. This could imply that the dispersive interactions in MAPbI₃ are such that higher-order interactions are indeed negligible compared to the pairwise contribution and also that the polarizability is largely isotropic (see Ref. 29 for further discussion).

Furthermore, we found that using the hybrid functional HSE, which improves the description of the electronic structure, did not result in improvements for the optimized unit-cell volume and bulk modulus, once it was combined with dispersive corrections. Indeed, the results from the PBE + TS, PBE + MBD, HSE + TS, and HSE + MBD calculations are all within the experimentally reported range for these two important structural parameters. Since our data, presented in Table I, also confirm the effect of SOC for the structure and binding to be minor, and since the PBE + TS approach was shown to be very accurate for structural and mechanical properties of multiple HaP crystals,^{16,22–24,26,27} we conclude that using PBE + TS for investigating the structure and binding in static and dynamical calculations of HaPs is a reasonable choice.

In this context, it is worth noting that Bokdam et al.³⁰ suggested the choice of DFT functional to have a bigger impact on the structural properties of MAPbI₃ than the addition of dispersive corrections. While this is certainly true for electronic properties, and while our calculations showed some minor differences between the PBE and HSE results, we have shown that this clearly cannot diminish the important role of dispersive interactions in MAPbI₃. We further note that the different conclusions of our work and Bokdam et al.³⁰ could be related to the different points of reference used in either case: while we have chosen to consider experimental data on the lattice constants and bulk modulus, Bokdam et al.³⁰ considered energies calculated in the random-phase approximation (RPA) as a reference. Indeed, RPA energies can be very accurate for solids and are broadly applicable, but it is worth noting that “standard” RPA suffers from known deficiencies, such as potentially incorrect descriptions of short-range correlation.^29,64

Using the PBE + TS approach allowed for studying more complicated static and dynamic structural phenomena: First, we investigated the impact of MA orientation on the structural parameters of MAPbI₃, studying the extreme cases of either perfectly parallel or antiparallel MA molecules contained in a supercell. We found that only for the parallel orientation of MA do the octahedra retain the cubic symmetry, while for antiparallel MA molecules, they strongly tilt. Note that the lattice vectors were still constrained to the primitive unit cell, i.e., the volume change corresponds to a hydrostatic pressure in the system. Since the relative energy gain/cost due to these structural distortions depends on the unit-cell volume, this effect modified the obtained volume and bulk modulus. Hence, the interaction between MA and the inorganic atoms is important, especially because the calculations showed that the antiparallel orientations are actually energetically preferred.

The origin of these distortions is related to the interactions between the MA molecules and iodine atoms since the partially positive ammonium-part of the MA interacts stronger with the partially negatively iodide ion than the methyl-side. In Fig. 6, we illustrate that the nitrogen atoms of MA interact mainly with five of the neighboring iodine atoms, which results in a distortion of the octahedra such that the distances between the nitrogen and iodine atoms are maintained between 3.64 and 3.81 Å. These interactions seem to be the main driving force behind the MA-induced distortions of the octahedra in the antiparallel case, which are energetically favorable for the system. Due to symmetry, these interactions are canceled in the case of parallel MA orientation, and our geometry relaxations did not automatically adapt to the more favorable antiparallel scenario. This finding is also quite interesting in view of the fact that the experimentally determined lattice symmetry of MAPbI₃ above 327 K is almost perfectly cubic,^9,10 despite the fact that different MA orientations can induce lower energy structures.

The implications of these findings can be fully understood by means of fully unconstrained MD calculations at 400 K since in these, the MA molecules, as well as all other atoms, are allowed to move freely. The first important finding from these data was that inclusion of dispersive corrections is equally important to obtain reasonable structural parameters in MD calculations. Second, the results obtained from the MD simulations are quite surprising at the first sight since the unit-cell volume was found to be smaller than the one obtained from the 0 K static calculation. To understand this finding, consider that the time scales associated with MA motion are much faster than the ones corresponding to the octahedral distortions of the heavy Pb-I cage, which is included in the MD but absent in the static calculations. Our analysis further showed that the average nuclear positions of MA are such that the carbon and nitrogen atoms are essentially conjoined at this temperature, meaning that MA is disordered, which is well known.⁹ Hence, the inorganic atoms essentially respond to a time-averaged volume corresponding to the moving MA molecule. Due to the MA disorder, this volume is effectively smaller at 400 K than for a fixed pattern of MA orientations, allowing for shorter I—Pb—I bonds and hence smaller crystal volumes. Last but not least, this rationale also explains the finding from the MD data showing that the time-averaged octahedral symmetry is almost perfectly cubic at 400 K in agreement with the experiment: MA disorder implies that all possible octahedral distortions induced by MA-iodine interactions are statistically equally likely. Therefore, the cubic perovskite structure is maintained as the most energetically favorable average crystal structure, exhibiting all the electronic and optical properties that render MAPbI₃ so favorable for device applications.

Finally, the data contained in Fig. 5 showed that at smaller volumes, the MD-averaged nuclear positions exhibit octahedral distortions together with a more preferred orientational order of MA. This makes sense, considering that smaller unit-cell volumes correspond to smaller voids between the octahedra, which increases the organic-inorganic interactions,²⁴ and partially hinders free MA motion. Such tilting of the octahedra into an orthorhombic-like structure was also observed experimentally in pressure experiments¹⁸ for MAPbI₃, and the variation in the preferred direction for MA at different volumes was discussed also in a recent study reporting MD simulations.⁶⁵ Therefore, in agreement with previous findings, our study shows that even relatively mild changes in external pressure can result in large changes in the structure and binding in MAPbI₃.

In summary, we investigated the impact of various levels of DFT-related approximations for calculations of the structural and binding properties of the prototypical HaP material MAPbI₃. Our tests considered the effects of including different dispersive correction schemes, applying a hybrid functional, including SOC, and also addressed the role of dynamic effects in MD calculations. The data confirmed previous theoretical work showing that dispersive corrections are important for accurate calculations of MAPbI₃ and also highlight that applying a computationally much more expensive hybrid functional improves the description of structural and mechanical properties by only a small amount. From this, we conclude that the use of a semilocal functional, augmented by pairwise dispersive interactions, is a suitable choice when computing more complicated static as well as structural dynamical phenomena in HaPs. Applying this methodology to DFT-based MD calculations of MAPbI₃, we analyzed the dynamic effect of molecular motion and its interplay with the structure of and binding in MAPbI₃. From this analysis, we could rationalize microscopically the simultaneous occurrence of a preferred cubic octahedral symmetry and MA disorder.

Funding provided by the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the German Federal Ministry of Education and Research is acknowledged. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC).