Structure and Binding in Halide Perovskites: Analysis of Static and Dynamic Effects from Dispersion-Corrected Density Functional Theory

TL;DR

This study evaluates how different density functional theory approximations, especially dispersive corrections, influence the structural and binding property calculations of halide perovskite MAPbI3, emphasizing computational efficiency and dynamic effects.

Contribution

It demonstrates that dispersive corrections are essential for accurate structural predictions, while hybrid functionals offer limited improvements, supporting the use of simpler methods for complex calculations.

Findings

Dispersive corrections are crucial for accurate structural properties.

Hybrid density functional theory has minor impact on key observables.

Molecular dynamics reveal the role of molecular rotation in MAPbI3 structure.

Abstract

We investigate the impact of various levels of approximation in density functional theory calculations for the structural and binding properties of the prototypical halide perovskite MAPbI$_3$. Specifically, we test how the inclusion of different correction schemes for including dispersive interactions, and how in addition using hybrid density functional theory, affects the results for pertinent structural observables by means of comparison to experimental data. In particular, the impact of finite temperature on the lattice constants and bulk modulus, and the role of dispersive interactions in calculating them, is examined by using molecular dynamics based on density functional theory. Our findings confirm previous theoretical work showing that including dispersive corrections is crucial for accurate calculation of structural and binding properties of MAPbI$_3$. They furthermore highlight that using a computationally much more expensive hybrid density functional has only minor consequences for these observables. This allows for suggesting the use of semilocal density functional theory, augmented by pairwise dispersive corrections, as a reasonable choice for structurally more complicated calculations of halide perovskites. Using this method, we perform molecular dynamics calculations and discuss the dynamic effect of molecular rotation on the structure of and binding in MAPbI$_3$, which allowed for rationalizing microscopically the simultaneous occurrence of cubic octahedral symmetry and MA disorder.

Figures

Tables

Table 1. Table 1: Volume of the primitive unit cell, V 0 subscript 𝑉 0 V_{0} , and bulk modulus, B 0 subscript 𝐵 0 B_{0} , obtained by fitting the DFT-calculated data with Eq.  1 , using various methods applied to the primitive unit-cell of MAPbI 3 subscript MAPbI 3 \mathrm{MAPbI_{3}} .

	PBE	PBE+TS	PBE+MBD	PBE+HI	HSE	HSE+TS	HSE+MBD	HSE+HI	PBE+SOC	PBE+TS+SOC	Experiment
ine $V_0$ [Å3] 111Errors are between 0.1 and 0.4 Å3	272.9	256.2	257.1	262.0	266.6	252.9	252.1	258.1	274.4	256.3	247-253 Stoumpos, Malliakas, and Kanatzidis (2013); Baikie et al. (2013); Feng (2014)
$B_0$ [GPa] 222Errors are between 0.2 and 0.4 GPa	10.6	15.7	13.6	14.3	11.1	16.4	14.7	14.2	9.6	15.0	12-16 Ferreira et al. (2018); Rakita et al. (2015)

Table 2. Table 2: Cell volume, V 0 subscript 𝑉 0 V_{0} , and bulk modulus, B 0 subscript 𝐵 0 B_{0} , obtained by fitting the DFT-calculated data with Eq.  1 , using PBE and PBE+TS applied to a 2x2x2 supercell of MAPbI 3 subscript MAPbI 3 \mathrm{MAPbI_{3}} with differently orientated MA molecules. Note that in order to improve the comparison, we here report V 0 / 8 subscript 𝑉 0 8 V_{0}/8 .

	parallel MA	anti-parallel MA
ine	PBE
$V_0$ [Å3] 111Errors are between 0.1 and 0.3 Å3	272.9	267.5
$B_0$ [GPa] 222Errors are between 0.4 and 0.4 GPa	10.8	11.2
	PBE+TS
$V_0$ [Å3] 111Errors are between 0.1 and 0.3 Å3	256.3	253.8
$B_0$ [GPa] 222Errors are between 0.4 and 0.4 GPa	15.7	14.7

Table 3. Table 3: Cell volume, V 0 subscript 𝑉 0 V_{0} , and bulk modulus, B 0 subscript 𝐵 0 B_{0} , obtained by fitting the DFT-calculated data with Eq.  1 , using PBE and PBE+TS applied to a 2x2x2 supercell of MAPbI 3 subscript MAPbI 3 \mathrm{MAPbI_{3}} computed along an MD trajectory of 20 ps. Note that in order to improve the comparison, we here report V 0 / 8 subscript 𝑉 0 8 V_{0}/8 .

	PBE	PBE+TS
ine $V_0$ [Å3] 111Errors are 0.1 Å3	267.3	248.1
$B_0$ [GPa] 111Errors are 0.1 Å3	8.8	13.0

Equations

$$E(V)=E_0+\frac{B_{0}V}{B_0^{\prime }}\left(\frac{(V_0/V{)}^{B_0^{\prime }}}{B_0^{\prime }-1}+1\right)-\frac{V_0{B}_0}{B_0^{\prime }-1}.$$

Full text

††preprint: AIP/123-QED

I Introduction

In recent years, the hope of cheap, solution-based solar cells with efficiencies above 20% has fueled the scientific research in halide perovskites (HaPs). Park (2013); Snaith (2013); Kazim et al. (2014); Stranks and Snaith (2015); Zheng et al. (2015); Brenner et al. (2016); Li et al. (2017) 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_{3}$, where the $B$-site is a heavy metallic element such as lead, and halide ions ($X$) form octahedra centered around the $B$-site. Saparov and Mitzi (2016) 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 $\mathrm{MAPbI_{3}}$ 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. Poglitsch and Weber (1987); Stoumpos, Malliakas, and Kanatzidis (2013) Importantly, HaPs are compared to other well established semiconductors Egger et al. (2018) and perovskite materials very soft mechanically, Feng (2014); Sun et al. (2015); Rakita et al. (2015); Létoublon et al. (2016); Kabakova et al. (2018); Elbaz et al. (2017) which is interesting fundamentally and 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. Capitani et al. (2016); Flores-Livas et al. (2018); Zhang, Zeng, and Wang (2017) 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. Feng (2014); Wang et al. (2013); Egger and Kronik (2014); Motta et al. (2015); Li and Rinke (2016); Faghihnasiri, Izadifard, and Ghazi (2017); Kabakova et al. (2018); Motta, El-Mellouhi, and Sanvito (2016); Egger (2018) 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.Kronik and Tkatchenko (2014); Hermann, DiStasio, and Tkatchenko (2017) In particular, a recent study has concluded that correcting for dispersive interactions in calculations of HaPs does not improve their description. Bokdam et al. (2017) 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, Perdew and Levy (1983); Sham and Schlüter (1983) which can be improved by applying computationally more costly hybrid functionals in the generalized Kohn-Sham framework. Seidl et al. (1996); Zhang et al. (2011) For classical inorganic semiconductors, screened hybrid functionals were shown to improve also the description of lattice constants compared to GGA-based approaches. Heyd and Scuseria (2004); Paier et al. (2006); Marsman et al. (2008) 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.Bokdam et al. (2017) Third, it is well established that spin-orbit coupling (SOC) due to the presence of the lead atom in $\mathrm{MAPbI_{3}}$ 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. Even et al. (2013); Whalley et al. (2017) Lastly, since HaPs are mechanically soft, they exhibit large structural dynamical effects, such as massive ionic displacements and molecular rotation around room temperature, Frost and Walsh (2016); Egger, Rappe, and Kronik (2016) 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 $\mathrm{MAPbI_{3}}$. 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 $\mathrm{MAPbI_{3}}$.

II Methods and Computational Setup

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) Tkatchenko and Scheffler (2009) and with an iterative Hirshfeld partitioning (HI) Bučko et al. (2014, 2013) as well as the many-body dispersion (MBD) method. Tkatchenko et al. (2012); Ambrosetti et al. (2014); Bučko et al. (2016) 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. Bučko et al. (2014) 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. Tkatchenko et al. (2012)

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 Perdew, Burke, and Ernzerhof (1996) and HSE functionals. Heyd, Scuseria, and Ernzerhof (2003, 2006) In the HSE functional, the short-range exchange energy is calculated as a mix of PBE and exact exchange, Heyd, Scuseria, and Ernzerhof (2003, 2006) which was often found to improve the description of electronic-structure and structural properties of semiconductors. Janesko, Henderson, and Scuseria (2009); Heyd and Scuseria (2004); Paier et al. (2006); Marsman et al. (2008) 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, Hobbs, Kresse, and Hafner (2000) 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). Birch (1947); Murnaghan (1944); Fu and Ho (1983) It describes the free energy of a solid as a function of the volume of the unit cell, $E(V)$, at isothermal conditions as:

[TABLE]

$E_{0}$ is the free energy at the equilibrium volume, $V_{0}$, and $B_{0}$ and $B^{\prime}_{0}$ 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 Stoumpos, Malliakas, and Kanatzidis (2013) 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_{0}$ and $B^{\prime}_{0}$ 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 $\mathrm{MAPbI_{3}}$ as has been reported experimentally, see Fig. 1a. Stoumpos, Malliakas, and Kanatzidis (2013) 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 2x2x2 pattern (see Fig. 1a), 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 $\mathrm{MAPbI_{3}}$. 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^{\prime}$ were quite large in the case of the MD calculations, in fitting Eq. 1 we set $B^{\prime}$ 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, Kresse and Joubert (1999) as implemented in the VASP code. Kresse and Furthmüller (1996) For unit cell calculations, a 6x6x6 grid of k-points centered around the $\Gamma$-point was used, for static and dynamic calculations that considered the supercell, the grid was reduced to a sampling of 3x3x3. 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, Bučko, Hafner, and Ángyán (2005) 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 one the atoms were below 10 meV per Å. In the canonical (NVT) MD simulations, a $\mathrm{MAPbI_{3}}$ 2x2x2 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. Lahnsteiner et al. (2016)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 Lahnsteiner et al. (2016) and trajectory length. The latter was checked explicitly by adding another 10 ps to the simulation, which resulted in only insignificant changes of the calculated observables. Schematic representations of the structures were generated with VESTA. Momma and Izumi (2008)

III Results

Fig. 2 shows the energy change due to volume change, i.e., $\Delta E(V)=E(V)-E_{0}$, for the static calculations of the $\mathrm{MAPbI_{3}}$ unit cell, calculated by using the PBE and HSE functionals augmented by different dispersion corrections. The optimized volume ($V_{0}$) and bulk modulus ($B_{0}$), obtained from the fit, are listed in Table 1. 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_{0}$ of 256.2 Å3 (PBE+TS) and 257.1 Å3 (PBE+MBD) which are close to the experimental range of 247-253 Å3. In contrast, the bare PBE calculations are far off (272.9 Å3) the experimental range, in agreement with literature findings, Feng (2014); Wang et al. (2013); Egger and Kronik (2014); Motta et al. (2015); Li and Rinke (2016); Faghihnasiri, Izadifard, and Ghazi (2017) and the PBE+HI (262.0 Å3) results basically lie in between the results of bare PBE and PBE+TS/PBE+MBD. The calculated $B_{0}$ data show similar trends, i.e., the experimentally reported values (12-16 GPa) 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 a 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\approx$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.

Fig. 2 also shows the results for regular and TS-corrected calculations using the HSE functional (see Table 1 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_{0}$ is still quite off the experimental result, in stark contrast to findings reported for conventional inorganic semiconductors. Heyd and Scuseria (2004); Paier et al. (2006); Marsman et al. (2008) 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_{0}$ value using the PBE+TS, PBE+MBD, HSE+TS, and HSE+MBD are all within $<5$ Å3, which is smaller than the range of experimentally reported values. Hence, we find that the improvement due to 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 $\mathrm{MAPbI_{3}}$. From these findings, we argue that using PBE+TS provides a reasonable choice for calculating structural and binding properties of HaPs such as $\mathrm{MAPbI_{3}}$. For completeness, we note that calculations including SOC showed very little difference to those without, as can be seen in Table 1.

It is well-known that the electronic interaction between MA and the inorganic ions in $\mathrm{MAPbI_{3}}$ 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 $\mathrm{MAPbI_{3}}$. In order to test this, we varied the MA orientation in a 2x2x2 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 $\mathrm{MAPbI_{3}}$ containing one MA unit. When calculating the structural parameters of the 2x2x2 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 computational cost (see above and Table 1). 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 2 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 2) 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 octahdra retain 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_{0}/8$ changes from 272.7Å3, for the parallel orientation, to 265.8Å3, 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 $\mathrm{MAPbI_{3}}$. Since our DFT calculations show that the antiparallel orientation is preferred in  $\mathrm{MAPbI_{3}}$, 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 $\mathrm{MAPbI_{3}}$, 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, Chen et al. (2015); Li et al. (2018) 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 $\mathrm{MAPbI_{3}}$. Fig. 4 shows the entire distribution as well as the mean value of the PBE- and PBE+TS-calculated change in free energy of $\mathrm{MAPbI_{3}}$, $\Delta 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_{0}$ parameters provided in Table 3. Furthermore, a slight decrease of $V_{0}$ 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_{0}$, while in the case of PBE it is similar to the 0 K calculation of the primitive unit cell, for PBE+TS $B_{0}$ 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. Fig. 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 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.

IV Discussion

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 $\mathrm{MAPbI_{3}}$.Feng (2014); Wang et al. (2013); Egger and Kronik (2014); Motta et al. (2015); Li and Rinke (2016); Faghihnasiri, Izadifard, and Ghazi (2017) 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,Bučko et al. (2014, 2013) $\mathrm{MAPbI_{3}}$ 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,Egger and Kronik (2014) are screened by the surrounding dielectric environment. We considered a comparison of the $C_{6}$ parameters calculated from PBE+TS and PBE+MBD, which determine the dispersive correction in either case,Kronik and Tkatchenko (2014); Hermann, DiStasio, and Tkatchenko (2017) 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 $\mathrm{MAPbI_{3}}$ are such that higher-order interactions are indeed negligible compared to the pairwise contribution, and also that the polarizability is largely isotropic (see ref. Hermann, DiStasio, and Tkatchenko (2017) 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 experimentally reported range for these two important structural parameters. Since our data, presented in Table 1, 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, Egger and Kronik (2014); Li and Rinke (2016); Kabakova et al. (2018); Motta et al. (2015); Motta, El-Mellouhi, and Sanvito (2016); Egger (2018) 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. 2017 suggested the choice of DFT functional to have a bigger impact on the structural properties of $\mathrm{MAPbI_{3}}$ 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 $\mathrm{MAPbI_{3}}$. We further note that the different conclusions of our work and Bokdam et al. 2017 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. 2017 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. Toulouse et al. (2009); Hermann, DiStasio, and Tkatchenko (2017)

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 $\mathrm{MAPbI_{3}}$, 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 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 iodines, 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 cancelled in the case of parallel MA orientation, and our geometry relaxations did not automatically adapt to the more favorable antiparallel szenario. This finding is also quite interesting in view of the fact that the experimentally-determined lattice symmetry of  $\mathrm{MAPbI_{3}}$ above 327 K is almost perfectly cubic,Poglitsch and Weber (1987); Stoumpos, Malliakas, and Kanatzidis (2013) 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 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.Poglitsch and Weber (1987) Hence, the inorganic atoms essentialy 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 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 $\mathrm{MAPbI_{3}}$ 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,Li and Rinke (2016) and partially hinders free MA motion. Such tilting of the octahedra into an orthorhombic-like structure was also observed experimentally in pressure experimentsCapitani et al. (2016) for $\mathrm{MAPbI_{3}}$, and the variation in the preferred direction for MA at different volumes was discussed also in a recent study reporting MD simulations.Lahnsteiner et al. (2018) Therefore, in agreement with previous findings our study shows that even relatively mild changes in external pressure can result in large changes of the structure and binding in $\mathrm{MAPbI_{3}}$.

V Conclusion

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 $\mathrm{MAPbI_{3}}$. 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 $\mathrm{MAPbI_{3}}$, 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 $\mathrm{MAPbI_{3}}$, we analyzed the dynamic effect of molecular motion and its interplay with the structure of and binding in $\mathrm{MAPbI_{3}}$. From this analysis, we could rationalize microscopically the simultaneous occurrence of a preferred cubic octahedral symmetry and MA disorder.

VI Acknowledgements

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).

VII References