The nature of hydrogen-bonding interaction in the prototypic hybrid halide perovskite, tetragonal CH3NH3PbI3

In spite of the key role of hydrogen bonding in the structural stabilization of the prototypic hybrid halide perovskite, CH3NH3PbI3 (MAPbI3), little progress has been made in our in-depth understanding of the hydrogen-bonding interaction between the MA+-ion and the iodide ions in the PbI6-octahedron network. Herein, we show that there exist two distinct types of the hydrogen-bonding interaction, naming α- and β-modes, in the tetragonal MAPbI3 on the basis of symmetry argument and density-functional theory calculations. The computed Kohn-Sham (K-S) energy difference between these two interaction modes is 45.14 meV per MA-site with the α-interaction mode being responsible for the stable hydrogen-bonding network. The computed bandgap (Eg) is also affected by the hydrogen-bonding mode, with Eg of the α-interaction mode (1.73 eV) being significantly narrower than that of the β-interaction mode (2.03 eV). We have further estimated the individual bonding strength for the ten relevant hydrogen bonds having a bond critical point.

Scientific RepoRts | 6:21687 | DOI: 10.1038/srep21687 controlling the degree of the PbI 6 octahedral tilting through the steric size of the molecular cation. According to these two studies 19,20 , the optical bandgap can be reduced by decreasing the degree of the octahedral tilting, which, in turn, can be achieved by adjusting the degree of the hydrogen-bonding interaction between the halides and H atoms bonded to the MA group. Several other studies [21][22][23][24][25] also indicate the important role of the MA + -ion orientation and, thus, the hydrogen-bonding interaction in controlling the core properties of the MAPbX 3 -based perovskite solar cells, which includes the enhanced carrier diffusion length 21 , the ferroelectric photovoltaic effect 22 , and the interplay of the MA-dipole orientation with the stability of perovskite structure 24,25 .
In spite of the key role of the MA-dipole orientation and consequent hydrogen-bonding interaction, little progress has been made in our systematic understanding of (i) the stable configuration of MA + -ions in the perovskite unit cell and (ii) the nature and strength of the hydrogen bonding between the MA + -ion and the halide (X) ions in the PbX 6 -octahedron network. Herein, we show that there exist two distinct types of the hydrogen-bonding interaction in the tetragonal phase which is relevant to room-temperature performance of the prototypic MAPbI 3 -based solar cells. On the basis of symmetry consideration of the PbI 6 -octahedron network, we will predict the possibility of existence of two distinct chemical environments for the MA + -ion orientation in the tetragonal phase and computationally show that one of these two is responsible for the stable hydrogen-bonding interaction between the MA + -ion and the surrounding PbI 6 -octahedron cages.

Results and Discussion
Two Distinct Environments for the Organic-cation Orientation. Quarti et al. 24 computationally showed that a set of polar (ferroelectric-like) structures formed by a preferred MA + -ion orientation is more stable, in general, than a set of apolar (antiferroelectric-like) structures formed by an isotropic distribution of the MA dipoles, which indicates an important role of the MA + -ion orientation in the stability of the perovskite lattice. Molecular dynamics computations 26 and first-principles study 27 further showed that for both cubic and tetragonal phases, the MA + cations are oriented parallel to the facial direction of the inorganic cage. On the basis of these theoretical studies, one can describe the orientation of the organic cations (MA + ) in the tetragonal phase by 8 different initial orientations 24 . These eight orientations are described by two characteristic angles, θ and φ, and are graphically illustrated in Fig. 1a, where a, b, and c denote the three crystallographic axes of the tetragonal perovskite lattice. In the figure, θ defines the angle between the a-axis (i.e., x-direction) and the projection of the MA + -ion orientation within the ab-plane. Thus  Herein, A, B, C, and D represent the projection of the MA + -ion orientations on the a-b plane, as measured by the azimuthal angle θ. The orientation of the MA + -ion with respect to the a-b in-plane is represented by the tilting angle φ. According to Quarti et al. 24 the two optimum φ angles are ± 30°. However, we have found that the optimum φ angle depends sensitively on the hydrogen-bonding interaction mode (See the text for details). (b) Crystal structure of the high-temperature cubic Pm m 3 phase viewed from the c-axis (lefthand side). The corresponding Pb-I inorganic cage characterized by the C 4 -rotation axis and the mirror plane σ h perpendicular to the C 4 axis (right-hand side). (c) Crystal structure of the tetragonal I4/mcm phase viewed from the c-axis (left-hand side). The corresponding Pb-I inorganic cage characterized by the improper S 4rotation axis (right-hand side).
Scientific RepoRts | 6:21687 | DOI: 10.1038/srep21687 the θ angles of 45° (for A), 135° (for B), 225° (for C), and 315° (for D). With respect to the ab-plane, the MA cations have two symmetric preferred orientations, φ = ± 30° 24 , where φ is the tilting angle of the C-N bond axis with respect to the ab-in-plane (Fig. 1a). Figure 1b shows the crystal structure of the high-temperature cubic phase composed of the central PbI 6 -octahedron cage and the surrounding MA + ions. In the cubic phase which is represented by the Pm m 3 space-group symmetry, the corner-shared PbI 6 octahedral frame does not show any tendency of the octahedral tilting along all three directions, a, b, and c. Thus, the cubic phase is represented by a a a 0 0 0 in the Glazer's notation. Figure 1c shows the tilted three-dimensional structure of the room-temperature-stable tetragonal phase which belongs to the I4/mcm space group. In this tetragonal structure, the PbI 6 octahedra do not show any alternative tilting along the a-and b-axes but exhibit out-of-phase tilting along the c-axis, which is in accordance with the − a a c 0 0 tilt pattern in the Glazer's notation. Let us now consider the difference in the point-group symmetry of the PbI 6 -octahedron cage between these two relevant phases. In the cubic phase, the PbI 6 octahedral network belongs to O h point group which is characterized by the principal 4-fold rotation axis along the c-axis (C 4 ) and the mirror plane perpendicular to this C 4 axis σ ( h ; Fig. 1b). Owing to the C 4 symmetry, a set of the following four distinct orientations of the C-N bond axis is under the same chemical environment: {+ A, + B, + C, + D}. Similarly, a set of the orientations, {− A, − B, − C, − D}, at a given MA-site is chemically equivalent in the cubic phase. Owing to the σ h symmetry, however, the two orientations having the same θ angle but with two opposite φ values (e.g., + A and -A orientations; Fig. 1a) are under the same chemical environment. Thus, in the cubic phase, the PbI 6 -octahedron cage provides all eight possible orientations of MA, {± A, ± B, ± C, ± D}, with the same chemical environment. This symmetry argument is graphically illustrated in Figure S1 of the Supplementary Information.
In the room-temperature-stable tetragonal phase, on the contrary, the PbI 6 octahedral network belongs to D 2d point group owing to the  The unit-cell structure of MAPbI 3 with the marked four distinct MA-sites is depicted in Fig. 3. As displayed in Fig. 3a, the four MA dipoles (1, 2, 1′ , and 2′ ) are located at the same a-b plane. When the cell is viewed from the a-axis (Fig. 3b), the 1 st and 2 nd MA-sites are on the same a-b plane but the 3 rd and 4 th sites are located at a different a-b plane which is (c/2) away from the former a-b plane along the c-axis of the tetragonal I4/mcm cell. Thus, the distance between the 1 st and 2 nd sites (or equivalently, between the 3 rd and 4 th sites) is given by , where a is the a-axis lattice parameter.
Two Distinct Modes of Hydrogen-bonding Interaction. We have examined the above made proposition on the existence of two non-equivalent chemical environments by investigating the MA-ion orientation in the tetragonal phase on the basis of ab initio density-functional theory (DFT) calculations. We used the experimental lattice parameters ( = = . , = . ) a b c 8 85Å 12 64Å 16 as the input parameters of our DFT calculations and subsequently obtained the optimized local structures of MAPbI 3 by applying the structure relaxation (i.e., relaxation of the internal positions at a fixed unit-cell volume). However, the volume relaxation method also gives essentially the same DFT optimized results that include the Kohn-Sham (K-S) energy and the equilibrium tilting angle, φ.
We have chosen two orientations, + A and -A, to examine the existence of two non-degenerate chemical environments at a particularly chosen MA site. However, our discussion is also valid for other pairs of the MA orientations (e.g., + C & -C). The DFT optimized value of θ is 45° for both + A and -A orientations [ Fig. 1a] 24 . However, the optimum tilting angle (φ) which corresponds to the minimum in the K-S energy depends sensitively on the MA + -ion orientation: ~+ 22° for + A orientation and ~+ 5° for -A initial orientation. It is interesting to notice that the optimum relaxed tilting angle (φ) for the -A initial orientation is ~+ 5°, instead of yielding a negative value. This is quite surprising since the input φ value for the -A orientation (usually It can be shown that for the -A initial orientation, the inconsistency between the symmetry prediction and the DFT optimized result stems mainly from the hydrogen-bonding interaction between the MA + -ion and I − ions in the PbI 6 -octahedron network. More specifically, the DFT optimized result fully reflects the site-specific hydrogen-bonding effect. On the contrary, the symmetry prediction is purely based on the D 2d structural symmetry of the PbI 6 -octahedron network without considering this site-specific hydrogen-bonding interaction between the MA + -ion and I − ions. Because of this simplification, the symmetry prediction can only be used as an initial guideline. In actual DFT calculations, we have adopted the structure relaxation at a fixed unit-cell volume, instead of using the volume relaxation, by considering computational efficiency and cost. As mentioned previously, however, the structure relaxation gives essentially the same DFT optimized results as the volume relaxation method. The above described extraordinary result indicates that the particular MA-site chosen in the present DFT calculations strongly prefers the + A orientation to the -A orientation. Let us call this particular site as the 1 st MA-site, as shown in Fig. 3. Indeed, the calculated K-S energy difference between the two orientations is as large as 45.14 meV per MA-site. Thus, a set of the orientations, {− A, + B, − C, + D}, does not practically exist at the 1 st MA-site though the symmetry argument predicts its existence. Consequently, we end up with a positive equilibrium φ even if we use a negative input φ value for the -A initial orientation. In our calculations of the K-S energy for the + A orientation at the 1 st MA-site, we have chosen the site-dependent dipole configuration of [+ A,-A,+ A,-A] which denotes the MA + -ion orientations of + A, -A, + A, and -A at 1 st , 2 nd , 3 rd , and 4 th sites, respectively. It can be shown that this particular MA + -ion configuration corresponds to the symmetry-allowed lowest energy configuration (See Subsection "Remarkably Simplified Dipole Configurations by Considering Structural Symmetry. "). On the contrary, the [-A,+ A,-A,+ A] initial configuration was used to evaluate the K-S energy for the − A orientation at the same 1 st MA-site. Thus, the K-S energy difference between these two distinct dipole configurations, [+ A,-A,+ A,-A] and [-A,+ A,-A,+ A], is as high as 180.56 meV (= 45.14 × 4).
Let us define the hydrogen-bonding interaction mode that corresponds to the tilting angle (φ) of + 22° as the α-interaction mode. Similarly, let us denote the hydrogen-bonding interaction mode corresponding to the tilting angle (φ) of + 5° as the β-interaction mode. Recall that the input φ value for the β-interaction mode is negative although the relaxed value is positive, ~+ 5°. As mentioned previously, the K-S energy difference between these two tilting-angle states is 45.14 meV per MA-site (i.e., per formula unit). The α-interaction mode with φ value of ~+ 22° is structurally illustrated in Fig. 4a by showing the 1 st MA-site (at center) and the surrounding PbI 6 -octahedron cages. On the other hand, the β-interaction mode with φ value of ~+ 5° is structurally depicted in Fig. 4b. Herein, the apical (axial) iodine atoms in the PbI 6 23 . However, as indicated in Fig. 5a, the bandgap (E g ) at the zone-center Γ -point is significantly affected by the hydrogen-bonding mode. We have further examined the partial density-of-states (PDOS) to resolve the atomic-scale origin of this bonding-mode-dependent bandgap. As indicated in Fig. 5b, the conduction-band minimum (CBM) is characterized by the Pb 6p orbitals, which is irrespective of the hydrogen-bonding interaction mode. On other hand, the valence-band maximum (VBM) is featured by the Pb 6s and I 5p orbitals. A detailed analysis of the wavefunction-character indicates that the Pb 6p-I 5p* anti-bonding orbital corresponds to the overlapping at the CBM while the Pb 6s-I 5p* anti-bonding orbital represents the VBM. It is interesting to notice that in the case of the α-interaction mode, the PDOS of the Pb 6p z orbital at the CBM further penetrates into a lower energy region (down to 1.73 eV above the VBM; Fig. 5b). This lowers the CBM value with respect to the VBM, leading to the bandgap reduction in the case of the α-interaction mode.
In addition, the Pb 6p z orbital is expected to show a certain degree of the orbital overlapping with the apical I 5p* orbital under the α-interaction mode. A careful examination of the PDOS indeed shows that the PDOS for the apical I 5p* (ap) orbital is slightly higher than that for the equatorial I 5p* (eq) orbital near the CBM under the α-interaction mode (Fig. 5b). Owing to the slightly enhanced Pb 6p z -I 5p* orbital overlapping at the CBM, it is predicted that the angle between Pb-(ap)I-Pb under the α-interaction mode is closer to 180° than the corresponding angle under the β-interaction mode. Indeed, our ab initio DFT calculations showed that the Pb-(ap) I-Pb angle under the α-interaction mode (ω = 168.6°) is substantially closer to 180° than the Pb-(ap)I-Pb angle under the β-interaction mode (ω = 160.3°).   Characteristic-angle-dependent Kohn-Sham Energy. We have then examined the orientation-dependent K-S energy to assess whether the DFT optimized α-interaction mode corresponds to the most stable configuration (in a single unit cell) or not. The orientation of the C-N bond axis is determined by the following three characteristic angles: θ (azimuthal angle), φ (tilting angle), and χ (torsion angle), where χ defines the rotation angle of the C-N bond axis 17 . For a fixed MA + -ion orientation, both αand β-interaction modes have a common azimuthal θ-angle. For instance, θ = 45° for ± A-orientations (Fig. 1a). Thus, we have examined φor χ-dependent K-S energy. In Fig. 6a, the computed K-S energy is plotted as a function of the tilting angle, φ, which indicates the equilibrium φ-angle for the 1 st (or 3 rd ) site is + 22° and + 5° for αand β-interaction modes, respectively. In case of the torsion angle, the K-S energy for the β-interaction mode shows its maximum when χ is at 0° or 120° while the K-S energy shows its minimum when χ is at 60° (Fig. 6b). Contrary to this, the K-S energy shows a reverse trend for the α-interaction mode. In this case, a pronounced increase in the K-S energy occurs upon increase in the torsion angle (χ) from 0° to 60° or upon decrease in χ from 120° to 60° (Fig. 6b). This increase in the K-S energy can be understood in terms of the rupture of the relevant hydrogen bonds upon the torsion of the C-N bond axis from its equilibrium χ values, 0°, 120°, etc (Fig. 4a). For each interaction mode, the orientation-dependent energy is described by three characteristic variables, θ, φ, and χ. Under the thermodynamic equilibrium, the K-S energy should be its true minimum, simul taneously satisfying the two criteria: in the three-dimensional θ φ χ ( , , )-space. Thus, for the α-interaction mode, the equilibrium φ and χ angles deduced from Fig. 6 correspond to the most stable state in a single unit cell.
According to the computed K-S energy shown in Fig. 6b, the activation free-energy for the C-N bond rotation is 49.4 meV for the α-interaction mode while it is 16.9 meV for the β-interaction mode. This suggests that the net hydrogen-bonding strength in the α-interaction mode is much stronger than that in the β-interaction mode. We will quantitatively examine this important point in Subsection "Evaluation of Individual Hydrogen-bonding Strength". Since the room-temperature thermal energy is 25.7 meV, an effectively free torsional motion of the C-N bond axis is expected in the β-interaction mode but not in the α-interaction mode. According to the transition-state theor y 28 , the frequency of the free torsional rotation is estimated to be: for the β-interaction mode at 300 K. We are in a position to summarize the main difference in the organic MA + -ion orientation between the αand β-interaction modes: (i) The equilibrium φ-angle is + 22° in the α-interaction mode while it is + 5° in the β-interaction mode. (ii) The orientation relationship of the − NH 3 group in the β-interaction mode with χ = 60° (Fig. 6b) can be reproduced by the 180° rotation of the N-H N (3) bond axis of the − NH 3 group in the α-interaction mode (χ = 0° or 120°) along the c-axis. This can be identified by examining the two left-hand side illustrations of Fig. 4. On the other hand, both αand β-interaction modes have a common θ-angle, as mentioned previously. Let us now extend the above argument to the remaining three MA-sites in the tetragonal unit cell (Fig. 3). On the basis of the translational symmetry of the tetragonal MAPbI 3 cell, the above symmetry rule can be directly applied to the 3 rd site. In other words, the + A (+ C) orientation is much more stable than the -A (-C) orientation at the 3 rd site, regardless of the hydrogen-bonding interaction mode. Thus, the calculated φ-dependent K-S energy for the 1 st MA-site (Fig. 6a) can be extended to the 3 rd MA-site, as shown in Fig. 7a. On the contrary, the reverse is true for the 2 nd and 4 th sites: the -A (or -C) orientation with a negative tilting angle φ is much more stable than the + A (or + C) orientation at the 2 nd or 4 th site, which is regardless of the interaction mode. The computed φ-dependent K-S energy (Fig. 7b) clearly supports this conclusion. Let us extend this argument of the site-dependent MA orientation to the ± B and ± D cases. Considering the structural symmetry rule of {+ A, − B, + C, − D} for the 1 st MA-site, one can readily obtain the reverse conclusion for the ± B and ± D orientations. More specifically, the − B (or − D) orientation is much more stable than the + B (or + D) orientation at the 1 st and 3 rd sites. On the contrary, + B (or + D) orientation is much more stable than the − B (or − D) orientation at the 2 nd and 4 th sites of the tetragonal MAPbI 3 .

Remarkably Simplified Dipole Configurations by Considering Structural
If the MA dipoles are randomly oriented as in the case of the cubic phase, the number of maximum conceivable orientations of the four MA dipoles in the tetragonal unit cell is given by (2*4) 4 = 4096 with each orientation represented by characteristic θ and φ angles. Herein, '2' takes into account ± orientations for a fixed θ, '4' represents the four possible values of θ for a fixed φ, and the power-exponent, 4, takes into account the four distinct MA-sites. Owing to the above symmetry rule of dipole orientations, however, the number of possible orientations of the four MA dipoles in the tetragonal cell can be greatly simplified. To deduce this, suppose that the 1  Table 2. Similarly, we have 64 distinct dipole configurations for each occupancy of -B or + C or -D dipole at the 1 st MA-site. Thus, we have a total of 256 conceivable dipole configurations in the tetragonal unit cell under the α -hydrogen-bonding interaction mode. Exactly the same number of the dipole at all the relevant bond critical points (BCPs) by exploiting the so-called 'quantum theory of atoms in molecules (QTAIM)' 29 . In the QTAIM, the local electronic kinetic-energy density of a given quantum system can be expressed in terms of the first-order density matrix ′ ρ r r ( , ): Mata et al. 31 further correlated the hydrogen-bonding energy (E HB ) with ( ) r G at the BCP using the following relation:

HB BCP
The calculated bonding energy and length, together with the associated topological properties ρ ( ) , are listed in Table 1 for the 10 relevant H N I bonds that are directly involved in the hydrogen-bonding interaction 17 . In addition to this, all ten BCPs (five BCPs for each interaction mode) are marked with small circles in Figure S2 Table 1 practically satisfy the criteria of the hydrogen bonding. Among 10 different hydrogen bonds, three are prominent in their E HB values and bonding lengths (< 2.70Å). Previously, they are named 'the dominant hydrogen bonds' (Subsection "Two Distinct Modes of Hydrogen-bonding Interaction") and are: H N (1)I A (1) and H N (2)I A (2) in the α-interaction mode and H N (3)I E (1) in the β-interaction mode.
According to the computed results shown in Table 1, the net difference in the hydrogen-bonding energy ∆E [ ] HB between the two interaction modes is 43.87 meV (= 381.06-337. 19) per formula cell. This clearly supports the previously made conclusion that the tetragonal MAPbI 3 perovskite cell under the α-interaction mode is much more stable than the same perovskite cell under the β-interaction mode (Subsection "Two Distinct Modes of Hydrogen-bonding Interaction"). Moreover, the estimated bonding-energy difference by the QTAIM (43.87 meV) nearly coincides with the previously calculated K-S energy difference between the two interaction modes (45.14 meV). As indicated in Eqs. (2) and (3), the computed E HB value by the QTAIM depends on ρ ( ) r at the BCP. In the DFT, ρ ( ) r uniquely determines the external potential 34 , thus, the ground-state K-S energy that comprises all the interaction terms including the Hartree energy, the external interaction energy between the nucleus and electrons, and the exchange-correlation energy. Thus, the computed value of ∆E HB (43.87 meV) by applying the QTAIM can be viewed as the overall K-S energy difference between the two interaction modes (45.14 meV), rather than being interpreted as the difference in the pure hydrogen-bonding interaction energy between the two interaction modes.

Conclusion
On the basis of symmetry argument and DFT calculations, we have made the following conclusions on the tetragonal MAPbI 3 perovskite cell: (i) There exist two distinct types of the hydrogen-bonding interaction between the MA + -ion and the iodide ions in the PbI 6 -octahedron network. We named them αand β-interaction modes. (ii) The computed K-S energy difference between these two interaction modes is 45.14 meV per MA-site with the α-interaction mode being responsible for the stable hydrogen-bonding network. (iii) Based on the individual bonding-energy calculations by exploiting the QTAIM, we have shown that five distinct hydrogen bonds are effective in the tetragonal MAPbI 3 under the stable α-interaction mode. The net difference in the total hydrogen-bonding energy between these two interaction modes is 43.87 meV per MA-site, which nearly coincides with the K-S energy difference of 45.14 meV. (iv) We have further made a remarkable simplification in the   1st  2nd  3rd  4th  1st  2nd  3rd  4th  1st  2nd  3rd  4th  1st  2nd  3rd 4th MA-dipole configurations by imposing the structural symmetry rule and the tilting-angle-dependent K-S energy to the tetragonal MAPbI 3 cell.

Methods
We have performed ab initio density functional theory (DFT) calculations on the basis of the Perdew-Burke-Enzerhof generalized gradient approximation (PBE-GGA) 35,36 implemented with projector augmented-wave (PAW) pseudopotential 37 using the Vienna ab initio Simulation Package (VASP) 38,39 . To assess the effect of the van der Waals (vdW) interaction on the structure relaxation, we have performed all ab initio calculations using the internal parameters of the Grimme's DFT-D2 vdW as implemented in VASP 40 . Most of the DFT calculations were performed by adopting (i) a 4x4× 3 Monkhorst-Pack (M-P) k-point mesh 41 centered at the Γ -point and (ii) a 500-eV plane-wave cutoff energy. In the band-structure calculations, however, we have initially adopted a 3x3× 2 M-P k-point mesh to obtain a relaxed structure and subsequently used a 9x9× 6 k-point mesh to accurately assess the k-point-dependent Kohn-Sham energy. All the structural relaxations were performed with a Gaussian broadening of 0.05 eV. The ions were relaxed until the Hellmann-Feynmann forces on them were less than 0.01 eV•Å −1 . The topological analysis of electronic density contours was performed by suitably exploiting the AIM-UC program 42 .