The crucial role of density functional nonlocality and on-axis CH3NH3 rotation induced I2 formation in hybrid organic-inorganic CH3NH3PbI3 cubic perovskite

Effects of electronic nonlocality in density functional theory study of structural and energetic properties of a pseudocubic CH3NH3PbI3 are investigated by considering coherent rotation around C–N axis of a CH3NH3 cation. A number of truly non-local and semi-local exchange correlation density functionals are examined by comparing calculated structural parameters with experimental results. The vdW-DF-cx which takes into account the non-local van der Waals correlation and consistent exchange shows the best overall performance for density functional theory study of this system. Remarkable distinctions between results from vdW-DF-cx and those from PBEsol exchange correlation functionals are observed and indicate the need of including the non-local interaction in the study of this system, especially its dynamical properties. The obtained rotational barriers are 18.56 meV/formula and 27.71 meV/formula which correspond to rotational frequencies of 3.71 THz and 2.60 THz for vdW-DF-cx and PBEsol calculations, respectively. Interestingly, the maximally localised Wannier function analysis shows the hydrogen bonding assisted covalent character of two iodide anions at a moderate rotational angle which can lead to I2 formation and then material degradation.

. The first type is the reorientation of the C-N axis with estimated 3-14 ps time scale 16,17,19,23 . The second type is the much faster methyl and/or ammonium groups rotation about C-N axis with sub-picosecond time scale 25 that has recently been estimated at ~30 times faster 26 . Such fast-rotation cannot be precisely probed by experiments 23 . Moreover, the interplay between octahedron tilting and the MA dynamics still needs a clearer understanding, as it has recently been dubbed as the chicken-and-egg paradox by Li and Rinke 27 .
Since the emergence of HOIPs as photovoltaic materials in 2009 2 , density functional theory (DFT) based calculations have been adopted to study their properties and underlying physics. Most of the studies are based on the standard semi-local generalised gradient approximation which cannot retain the non-local nature of electron correlation. However, many suggest the importance of including the non-local van der Waals (vdW) interaction in the calculations to recover such many-body effect 24,[27][28][29][30][31] . Nevertheless, there are a number of methods for including the vdW-interaction implemented in the DFT based software packages. Each of them yields different levels of accuracy, especially, for molecular, sparse and hydrogen-bond materials. The detailed discussions on the methodology and accuracy of the vdW density functionals (vdW-DFs) can be found elsewhere [32][33][34] .
In this work, we examine the application of various vdW-DFs and discuss in detail the equilibrium structural configurations corresponding to the on-axis MA rotation obtained by structural relaxations. As a result, vdW-DF-cx non-local exchange correlation (xc) functional shows the best overall performance for this system with some remarkable differences in the atomic configurations comparing with those of PBEsol semi-local xc-functional. By adding the non-local effect, lowering and broadening of the rotational energy barrier are observed. Furthermore, maximally localised Wannier function (MLWF) analysis illustrates the formation of hydrogen bond assisted I-I covalent bonding at a specific MA rotational angle which could lead to the formation of neutral I 2 defect and so the material degradation.

Results and Discussion
Structural parameters. First, we examine the xc-functionals by calculating structural parameters of the pseudocubic MAPbI 3 unit cell. Since the assigned Pm 3 m symmetry of MAPbI 3 was derived from the average structure of the material at high temperature. Therefore, the average values of such calculated parameters, i.e. average lattice constant (a), average Pb-I bond length over three inequivalent I-atoms (PbI), average C-H bond length (CH), average N-H bond length (NH), and C-N bond length (CN) obtained by performing full structural relaxation, there is no any constraint, where the C-N axis is at the lowest enthalpy arrangement 14 , are used for the comparison with Weller et. al. 's accurate powder neutron diffraction experimental results 15 .
The structural parameters are grouped into two subgroups according to magnitudes of their absolute errors, i.e. a and PbI as one subgroup and CN, CH and NH as the other. The mean absolute error (MAE) of each subgroup is given in Table 1. The parameters a and PbI are in excellent agreement with experimental results (MAE(1) < 1%) for most of xc-functionals except for vdW-DF, vdW-DF2, and PBE. Notably, the vdW-DF-cx shows the best performance in predicting the lattice constant and the Pb-I distance, see MAE(1) in Table 1. This is in good agreement with the vdW-DF testing performed by Berland and Hyldgaard. They examined the accuracy of a number of vdW-DFs in predicting the binding energy and equilibrium separation in H-bond and dispersive-force dominated systems and found that vdW-DF-cx shows the best overall performance 32 . However, MAE(2) in Table 1 shows that errors in CN, CH, and NH are in approximately the same order for all xc-functionals, with an error of ~10% for CN and CH, and ~5% for NH. The disordered C-N axis orientation and the fast on-axis rotation of H-atoms could probably be responsible for the errors in structural parameters. The divergence from the experimental data might originate from the indistinguishability of the atomic position reported by the experiment. Since the position of C, N, H C and H N atoms cannot be experimentally captured due to the MA dynamics, their positions are approximately assigned to Wyckoff sites 4MM(100) and 1 of the Pm 3 m space group 15 . From the simulation viewpoint, by using only one pseudocubic primitive unit cell, the disordered nature of the organic cations is neglected, giving rise to the divergence from averaged value of the experimental data. Moreover, the effects of temperature which drive atoms from their zero-Kelvin equilibrium configurations are also not included. Although the simulations give large MAE(2), the lattice constant and the Pb-I bond length are in good agreement with those of the experimental results, thanks to their low dynamics. The latter validates the use of computational model for investigating the local structural and energetic responses to the rotation of H atoms about the C-N axis.

Energy barriers, Equilibrium atomic configurations, and H-bonds. The rotational energy barriers are
obtained by performing structural relaxation without any constraint after coherent H-atoms rotation (solid lines in Fig. 1(c)). The rotation of MA is discretised to 12 rotational steps with 10° step size owing to its C 3v symmetry, see Fig. 1 for atomistic structural model. Hereafter, only the results from vdW-DF-cx and PBEsol xc-functionals are shown as they are the representatives of non-local and semi-local xc-functionals respectively, results from all xc-functionals are provided in supplementary information. According to Fig. 1(c), the rotational potential energy, ≅ . meV/formula) and the vdW-DF-cx rotational potential energy curve is broader than that of PBEsol.
Lowering of the vdW-DF-cx energy barrier attributes to the nature of electronic interaction incorporated in the xc-functionals and the distinction of optimised structures. To understand these, we calculate θ ΔE( ) of the vdW-DF-cx-structures for θ = 0° to 60°, due to a symmetry of the energy curves, with applications of the following xc-functionals: vdW-DF-cx, PBEsol, vdW-DF-cx-ex (non-local correlation is turned off), and PBEsol-ex (PBEsol density gradient correlation is turned off) (see Fig. 2). The electronic interaction contribution can be observed by comparing , and the semi-local part of the vdW-DF-cx is chosen such that its PBEsol correlation is positive at θ = 10°, 20° while vdW-DF-cx correlation is negative at all rotational angles. Moreover, at these beginning rotational angles, the exchange parts of both semi-local and non-local xc-functionals are almost the same. The non-local correlation becomes crucially important at the moderate rotational angles (θ = 40°, 50°, 60°). The structural contribution can be perceived by comparing ΔE PBEsol with E PBEsol(vdW Str) Δ − . This contribution attributes to: i. unit cell deformation ( Fig. 3(a,b)), ii. I-displacement (octahedron distortion) (Fig. 3(c-e)), iii. MA arrangement and its internal distortion, and iv. hydrogen bonding (Fig. 4(b)). For instance, at 60 θ =°, (Fig. 2) is a consequence of: i. short PBE-optimised lattice vectors comparing to those of vdW-DF-cx-optimised structure while a difference in unit cell shearing is negligible, ii. notably large PBE-optimised I 2 atomic displacement in opposite of [100] direction comparing to that of vdW-DF-cx-optimised structures. Moreover, at this angle, only Hn 1 -I 1 (bc) , Hn 2 -I 1 (c) , and Hn 3 -I 2 (a) hydrogen bonds play the role (Fig. 4(b)), and those of PBEsol-optimised and vdW-DF-cx-optimised structures are almost the same. It is important to note here that the internal distortion of MA is negligible. Thus, Δ >Δ −

E E
PBEsol PBEsol(vdW Str) mainly attributes to unit cell and octahedron deformations. The lowering of rotational potential barrier is predominated by the non-local interaction keeping H-atoms and I-atoms further apart from each other.
The broadening of vdW-DF-cx rotational potential energy curve can be related to an abrupt change of H-atom positions, illustrated in Fig. 4(a), accompanying with I 1 and I 3 displacements in [001] and [100] directions, respectively. These differences in optimised structures reflect the noticeable distinction in force field calculated by semi-local and non-local xc-functionals. This is emphasised by the totally difference in changes of C-N axis altitudinal angle when the MA undergo on-axis rotation shown in Fig. 3(f). This indicates the need for theoretical revision of dynamical properties of this material in which the truly non-local interaction is included as its importance has recently been realised in spin-polarised system 35 and nonmagnetic transition metals 36 .
Furthermore, we examine the rotational potential energy when there are constraints imposed in structural relaxation process. The results are shown in Fig. 1(c). The vdW-DF-cx and PBEsol rotational energy barriers are 37.7 meV and 50.7 meV, respectively when the cell parameters are kept at the lowest energy configuration (B), and 50.6 meV and 47.0 meV, respectively, when both cell parameters and I-atomic positions are kept at the lowest energy configuration (C). We also calculate the free on-axis rotational frequencies corresponding to such energy barriers. By obeying the transition-state theory 30 , the frequency is given by: when ΔE is the rotational energy barrier. The obtained frequencies, at T = 330 K, for the relaxation model A, B, and C are 3.71 (2.60) THz, 1.83 (1.16) THz, and 1.16 (1.32) THz, respectively, for vdW-DF-cx (PBEsol) calculations. These rotational frequencies may help one experimentally justifies which MA rotational mechanism it is, if the local rotational frequency could be experimentally probed.
As we have seen in Fig. 4(b), even in this high symmetry phase the complexity of hydrogen bonding takes place in this material. Especially when MA undergoes rotation, the hydrogen bonding associated with a specific H-atom can be turned on/off or switches to bond with the other I-atom depending on the rotational angle. For example, Hn 1 strongly bonds with I 3 (ab) at θ ∈ [0°, 40°] then switches to strongly bond with I 1 (bc) , (a criterion used for strong H-bond is H ··· I ≤ 3 Å). To illustrate that, we perform the maximally localised Wannier function   37 . The complete MLWFs calculated at θ = 0°, 40°, 60° are given in supplementary information. Here, it is worth highlighting a MLWF that localised in vicinity of the I 1 (bc) obtained at θ = 60°. This MLWF mainly has I-p lone pair atomic orbital character and, also, shows covalent bonding character with the I 3 (ab) assisted by Hn 1 , see Fig. 5. Note that bond character identification of the hydrogen bond supported I − − I − interaction is based on the basic idea of the quantum theory of atoms in molecules formulated by RFW. Bader 38 that the I-p lone pair of the I 1 (bc) becomes delocalized into vicinity of the I 3 (ab) instead of as localized as those of the other I-atoms. This result is the first direct theoretical evidence supporting the recently experimental observation of interstitial iodide migration that leads to I 2 formation in tetragonal MAPbI 3 via e 2I I 2 2 → + − − and then material degradation with I 2 as a product 39,40 . Moreover, recent experiment has revealed that I 2 can further accelerate the degradation of the iodide perovskites 41 . Notably, our results show that I -I covalent bonding is supported by hydrogen bond formed when   the structure is at the top of rotational potential hill which mean that I 2 formation is induced by MA rotation activated by thermal energy.

Conclusions
We have performed density functional theory (DFT) based calculations to investigate the equilibrium atomic configurations of the pseudocubic MAPbI 3 hybrid perovskite with respect to on-axis coherent methylammonium rotation. First, we assess the accuracy of a number of semi-local and non-local van der Waals exchange correlation density functionals (xc-functional) by calculating the structural parameters and compare them to the experimental values. According to the results in Table 1, vdW-DF-cx xc-functional shows the best overall performance for structural property DFT calculation of this material. Second, we report details of structural and energetic responses to the MA rotation obtained through the full structural relaxations. Comparison of results calculated with application of vdW-DF-cx with those of PBEsol xc-functionals shows noticeable lowering and broadening of vdW-DF-cx rotational potential energy curve. The non-local vdW correlation becomes crucial at the rotational angles (at θ = 40°, 50°, 60°) where the potential energy is much higher than the lowest energy (at θ = 0°). That manifest an inevitability of including truly non-local xc-functional in DFT study of this system. Also, the rotational potential barriers and corresponding rotational frequencies are calculated for the cases that the unit cell is fixed, and both unit cell and I-atom positions are fixed in the relaxation processes. Finally, a maximally localised Wannier function analysis reflects that the I -I covalent bonding can be formed at θ = 60° with the assistance of hydrogen bonding. This can lead to I 2 formation and material degradation.

Methods
Static DFT calculation is performed with the pseudopotential and plane wave method implemented in Quantum-Espresso package 42 . The ion cores are represented by original scalar relativistic ultrasoft pseudopotentials taken from PSlibrary project 43 . Spin-orbit coupling (SOC) effects are not included in this work, as previous study suggested that SOC does not significantly affect ground state properties of this material 44 . The kinetic energy cut-off is set to 70 Ry. The unshifted k-point mesh is gridded into 6 × 6 × 6 by Monkhorst-Packs scheme 45 . This configuration yields convergence of total energy within 2 meV. The nonlocal vdW exchange correlation functionals considered in this work are as follows: rvv10 46 , vdW-DF 47,48 , vdW-DF2 49 , vdW-DF-ob86 50 and vdW-DF-cx 51 . The Perdew-Burke-Ernzerhof (PBE) 52 and PBEsol 53 functionals are also applied for comparison. The relaxation algorithm of atoms and lattice parameters used is Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm with a force/atom tolerance equal to 0.001 Hartree/Bohr. The initial structure is the lowest enthalpy configuration reported by Brivio et al. study 14 (available online 54 ) re-optimised with our settings. The calculations of exchange energies are achieved by setting of the Quantum-Espresso input variable, i.e. input_dft, as following: input_dft = 'sla + pw + cx13' , and input_dft = 'sla + pw + psx + nogc' for vdW-DF-cx and PBEsol exchange (with local correlation) density functionals respectively. For maximally localised Wannier function (MLWF) calculation, the k-space is meshed by 8 × 8 × 8 grid. Twenty-five Bloch wave functions are projected onto twelve atom-centred MLWFs with convergent criteria of the spread function equal to 10 −10 Å 2 .