Rotation mechanism of methylammonium molecules in organometal halide perovskite in cubic phase: An ab initio molecular dynamics study

Organometal halide perovskites such as methylammonium (CH₃NH₃ or MA) lead iodide (MAPbI₃)¹ have a wide range of applications in solar cells and sensors because of their excellent photoexcitation properties. In particular, the large diffusion lengths of photocarriers^2–4 allow the fabrication of efficient photovoltaic devices.^5,6 It is also known that photocarriers behave as free electrons and holes,^7,8 which are generated by polarized domain structures associated with MA molecules. Previous studies have suggested that dipole moments are generated by the extensive distortion of the perovskite structure caused by the rotation of MA molecules, which in turn provide polarized local domains to induce efficient charge separation,^9–15 although the molecules themselves do not control the electronic structures of the conduction- and valence-band edges.^12,16–19

While elucidating how an MA molecule rotates in the crystalline PbI₃ cage (Fig. 1(a)) is important for the fundamental understanding of the polarized domain formation, to probe it experimentally is quite difficult. X-ray diffraction (XRD) and nuclear magnetic resonance experiments were successful in the identification of the structure of PbI₃ cage and reported that a phase transition from tetragonal to cubic occurs in MAPbI₃ bulk at among 315-331 K.^1,20,21 However, those could not identify the orientation of MA molecules. On the other hand, some of neutron scattering (NS) experiments for MAPbI₃ bulk^22,23 reported potentially preferable motions and the orientation of MA molecules, but it is still unclear why such characteristic reorientations can occur.

Under such circumstances, theoretical studies based on ab initio molecular dynamics (AIMD) simulations would provide crucial information on the dynamical behavior of MA molecules.^{10,19,22,24–27} By comparing the dynamics of MA molecules in cubic-MAPbI₃ with and without the motion of the Pb–I lattice, Mosconi et al. showed that reorientation of the MA molecule requires a strong correlation between the MA molecule and the PbI₃ cage.²⁴ On the other hand, Frost et al.^10,22 applied cubic symmetry-assisted analysis (CSAA) to their AIMD results to investigate preferable crystallographic orientation of MA molecules. It is impossible to cover all of rotating sphere of C–N bond in PbI₃ cage with sufficient statistics within the feasible AIMD calculation time. Using the O_h symmetry of cubic-MAPbI₃, however, the sphere can be folded into 1/48, and accordingly the statistics is enhanced by 48 times. Frost et al. analyzed the orientation by folding the sphere into the limited domain consisting [100], [110], and [111] directions, and then concluded that a C–N bond is preferably parallel to the [100] direction (or $⟨100⟩$ directions). Such analysis using cubic symmetry has been developed and widely applied to the orientation of MA molecules by Even et al.^28,29 This anisotropy in the direction of C–N bond would be due to hydrogen bonds between H atoms in the MA molecule and I atoms in the PbI₃ cage. As suggested from previous static ab initio calculations,^17,25,30 the H–I hydrogen bonds play a vital role in the stability of the atomic structures. Carignano et al.²⁷ suggested that MA molecules rotate with the formation and cleavage of many H–I hydrogen bonds from radial distribution functions (RDFs) obtained from their AIMD simulation.

While the AIMD simulations introduced above are of great interest to understand key factors about rotational motions of MA molecules from a dynamic perspective, the analysis methods were limited to statistical ones. By tracing the reorientation of MA molecules and the formation of associated H–I hydrogen bonds, such anisotropic dynamics of MA molecules could be revealed.

In this study, we investigate the role of H–I hydrogen bonds to guide the rotational dynamics of MA molecules in cubic-MAPbI₃ at room temperature by performing AIMD simulations, in conjunction with the statistical methods used in previous AIMD studies. The reason for studying cubic-MAPbI₃ is our interest in MAPbI₃ films. Oku et al.^31,32 demonstrated the existence of a cubic phase in MAPbI₃ film at room temperature. Also, Milot et al.³³ reported that the phase transition temperature to cubic in MAPbI₃ film becomes less than that in bulk. Since MAPbI₃ film is quite important for practical applications such as solar cells, we deal with cubic-MAPbI₃ at room temperature to obtain insights into the rotational mechanism of MA molecules in film structure. However, as mentioned below, some of the structures reported as tetragonal phase are quite close to that of cubic phase. Thus, since there are much more experiments and theoretical studies dealing with MAPbI₃ bulk than those of film, we compare our AIMD result with the previous studies of bulk to investigate whether our approach from the cubic structure is reasonable. We first show the statistical results such as RDFs and the result of CSAA to verify the structure of our MAPbI₃ and to discuss preferable orientation of MA molecules, and then explain possible atomistic mechanisms underlying these results by focusing on bonding states of H–I hydrogen bond around MA molecules. While previous studies defined H–I hydrogen bonds by imposing cutoff lengths, here, we do it using charge density between the H and I atoms obtained by population analysis,^34,35 to avoid the arbitrariness in the selection of cutoff lengths.

We performed AIMD simulations based on density functional theory (DFT).^36,37 The exchange-correlation energy was included within the generalized gradient approximation (GGA)³⁸ and the interactions between ions and electrons were described by the projector-augmented-wave (PAW) method.^39,40 Projector functions were generated for the 1s state of hydrogen (H) atoms, the 2s and 2p states of carbon (C) atoms, the 2s and 2p states of nitrogen (N) atoms, the 5s, 5p, and 5d states of iodine (I) atoms, and the 6s, 6p, and 6d states of lead (Pb) atoms. Plane-wave expansion was performed at the Γ point in the Brillouin zone. The plane-wave cutoff energy of the electronic pseudo-wave-functions was 25 Ry and that of the pseudo-charge-density was 250 Ry. The DFT-D method was used for the semiempirical correction of the van der Waals interaction.⁴¹ The energy functional was minimized using a preconditioned conjugate-gradient method.^42,43 We also used the quasi-Newton method⁴⁴ to perform structural optimization. The AIMD simulation was performed at a temperature of 300 K and ambient pressure using the NPT ensemble.^45,46 We used two systems consisting of $2×2×2$ and $3×3×3$ cubic unit cells of MAPbI₃ (a total of 96 and 324 atoms, respectively) with periodic boundary conditions. In the initial configuration, we arranged the C–N bonds of all of the MA molecules parallel to the [11$1¯$] direction in cubic boxes with dimensions of $12.66Å×12.66Å×12.66Å$ and $18.99Å×18.99Å×18.99Å$⁠, for $2×2×2$ and $3×3×3$ systems, respectively. The $3×3×3$ system is shown in Fig. 1(b). The equations of motion were numerically solved by an explicit reversible integrator⁴⁷ with a time step of Δt = 0.484 fs. The total calculation times were 48.4 ps and 19.8 ps for $2×2×2$ and $3×3×3$ systems, respectively. However, the structure of $2×2×2$ system could not reproduce those of MAPbI₃ bulk and film because the calculated total RDF did not coincide with experimental RDFs as shown in Fig. S3(a) in the supplementary material. Thus, we mainly discuss the results obtained from $3×3×3$ system.

We used electronic population analysis^34,35 to clarify the changes in the bonding properties of atoms associated with hydrogen bonds. By expanding the electronic wave functions in an atomic-orbital basis set,^48,49 we obtained the gross population Z_i(t) for each atom and the bond-overlap population (BOP), O_ij(t) for each atomic pair, which are based on a formulation generalized to the PAW method.⁵⁰ 

The Mulliken charge Q_i(t) was then obtained as the difference between the number of valence electrons of an isolated neutral atom $Zi0(t)$ and the value of the gross population Z_i(t),

We estimated the charges of atoms from Q_i(t).

Meanwhile, since BOPs are often used to define H–I hydrogen bonds in this study, we give a more detailed description of this method. O_ij(t) gives a semiquantitative estimate of the strength of covalent bonds between atoms. In the PAW method, the overlap integral Oⁿ corresponding to nth eigenstate $ψn$ is defined using an atomic-orbital basis set $ϕμ$} as follows:

with

and

where $p∼μ$ and $qμν$ represent the projector functions and the augmentation charges, respectively, in the PAW method.⁴⁰ The BOP O_ij for a pair of ith and jth atoms is obtained by summing $Oμνn$ as follows:

where f_n is the electronic occupation number of the nth eigenstate. Note that a bonding state exists between ith and jth atoms if O_ij shows a positive value. In contrast, when O_ij shows a negative value, an antibonding state exists.

To estimate the error in the projection of the electronic wave functions onto atomic-orbital basis functions, we calculated the charge spillage S_Q,⁴⁸ 

with

where N represents total number of electrons and $χn$ is the projected function corresponding to $ψn$⁠. In other words, S_Q represents the fraction of total electronic charge spilled when projecting the calculated eigenstates onto atomic-orbital basis functions. The resulting S_Q was only 0.1% in our simulation, indicating that our basis set was high quality.

We used our own AIMD code, which has been implemented on parallel computers using the message passing interface library for interprocessor communications.

Here, we investigate the atomic structure by radial distribution functions (RDFs) $gα−β$(r) to verify that the structure is comparable to those from the experiments and to clarify the distance between H and I atoms. The solid, dashed, and dotted curves in Fig. 2(a) are the total RDFs obtained from the simulation with X-ray form factors and those for MAPbI₃ bulk and film from an XRD experiment by Choi et al.,⁵¹ respectively. The four peak positions and their widths are in reasonably good agreement with those of both bulk and film structures, which indicates that the lattice structure is correctly described in the simulation. Since we use cubic-MAPbI₃, the corresponding Glazer notation⁵² is a⁰a⁰a⁰, i.e., perovskite’s tilting motion cannot be reproduced in our simulation. Also, the calculated distribution of I–Pb–Pb–I dihedral angle has a peak around $0°$ as shown in Fig. S1 in the supplementary material. This means that an ideal cubic phase is formed. The reason for the total RDF matches that of bulk structure as well as that of film one is that the structures used in the experiment were quite close to that of cubic. In fact, Choi et al.⁵¹ showed the refined tetragonal lattice parameters (a and c), where the reduced unit cell parameter ratios⁵³ $c/2a$ were only 1.006 (for bulk) and 1.002 (for film). For other experiments, while some of them showed the ratios more than 1% (e.g., 1.011¹ and 1.010⁵⁴), there were also works that reported the ratios less than 1%, e.g., 1.009⁵³ and 1.005.²³ Thus, although it would not be applied to all the reported MAPbI₃ bulk, it is possible that the cubic-MAPbI₃ is in good approximation to represent their structure.

The total RDFs from the XRD experiment provide only the structure consisting of Pb and I atoms. Both the peak position and width can be explained by g_Pb–I(r), g_Pb–Pb(r), and g_I–I(r) obtained from the simulation (Fig. 2(b)), indicating that the experimental RDF is determined only by the atomic structure of the inorganic Pb–I lattice. The broad peaks indicate large fluctuation in the atomic distances (Pb–I, Pb–Pb, and I–I) and angles (Pb–I–Pb and I–Pb–I) in the PbI₃ cage, which are induced by the interaction between the cage and MA molecule, as suggested in the previous studies.^{6,21,30,55,56} To clarify the atomic structure of MA molecules, the calculated partial RDFs between the PbI₃ cage and MA molecule are shown in Figs. 2(c)–2(f). Hereafter, H atoms bonded to C and N atoms are referred to as H_C and H_N atoms, respectively. From Figs. 2(c)–2(f), the shortest interatomic distance between the MA molecule and the PbI₃ cage is between H_N and I atoms. The first peak of g_HN–I(r) is located at 2.9 Å. The H_C–I interatomic distances are much longer, and the first peak of g_HC–I(r) is at 3.5 Å. Similar to the relationship between the H_N and H_C atoms, the first peak of g_N–I(r) (3.7 Å) is less than that of g_C–I(r) (4.2 Å). On the other hand, all of the interatomic distances between Pb atoms and atoms in the MA molecule are longer than those between I atoms and atoms in the MA molecule, i.e., the first peaks of g_HN–Pb(r), g_HC–Pb(r), g_N–Pb(r), and g_C–Pb(r) are located at 4.1, 4.5, 4.9, and 5.3 Å, respectively. From the positional relationship mentioned above, the most frequent interactions would occur between H_N and I atoms, which coincides with the results of previous ab initio studies.^17,27,30 As suggested in these literatures, hydrogen bonds are formed between H_N and I atoms and are involved in the stable configuration of MA molecules in PbI₃ cage. Thus, by tracing the formation process of H_N–I hydrogen bonds during our simulation, we attempt to understand the rotational motion of MA molecules. However, since the previous studies defined H_N–I hydrogen bonds by imposing cutoff lengths, there was arbitrariness in the selection of the length. Instead, we define them from the charge density existing between H_N and I atoms obtained by population analysis described in Sec. IV. Also, we define hydrogen bonds between H_C and I atoms to compare with those between H_N and I atoms there.

In this section, we define H_N–I and H_C–I hydrogen bonds based on population analysis. Specifically, the BOP, which is described in Section II, is used. Although the BOP semiquantitatively represents covalent bond strength originally, it is valid to analyze chemical reactions involved with hydrogen bonds as demonstrated by our previous studies.^57,58

Hereafter, Subsection IV A describes the time-averaged distribution of BOPs between H_N–I and H_C–I to compare the strengths of hydrogen bonds. In addition, the time-averaged distributions of the Mulliken charges of H_N and H_C atoms are also shown to confirm that the results of the BOPs are reflected in the atomic charges. Next, the corresponding average numbers of hydrogen bonds and average bonding distances at any given time during our simulation are given in Subsection IV B.

Figure 3(a) shows the time-averaged distributions of the BOPs of H_N–I and H_C–I. Because the H_N–I interatomic distance is shorter, the BOP of H_N–I is greater than that of H_C–I. Namely, H_N atoms form more bonds with I atoms. The time-averaged distributions of the Mulliken charges of H_N and H_C atoms are shown in Fig. 3(b). The distribution of H_N atoms is slightly spread in the negative direction. H_N atoms have more positive charges than H_C atoms because of the negative charge transfer from I to H_N atoms with the formation of hydrogen bonds. The distribution of H_C atoms has a sharper peak, indicating the formation of few hydrogen bonds.

Based on the BOP, the distributions of the interatomic distances of the H_N–I and H_C–I hydrogen bonds were calculated, where two atoms are considered to be bonded when the corresponding BOP shows a positive finite value. The time-averaged distributions of the interatomic distances of H_N–I and H_C–I are shown in Fig. 4. The peaks are located at 2.9 Å, which corresponds to the most frequent bonding distance, and shorter bonds can be considered to be strong H–I hydrogen bonds. The interatomic distances of H_N–I and H_C–I are in the ranges of 2.2–3.7 and 2.5–3.2 Å, respectively. H_N atoms can form hydrogen bonds with a wide range of interatomic distances. In addition, the integrated values correspond to the time-averaged numbers of I atoms forming hydrogen bonds with H_N and H_C, which are approximately 1.3 and 0.1, respectively. Namely, all of the H_N atom bonds form at least one hydrogen bond with I atoms, while H_C atoms hardly form any hydrogen bonds. The reason why H_N atoms bond to more than one I atom is that the H_N atom instantaneously forms two hydrogen bonds when molecular rotation occurs. We therefore conclude that the rotation of MA molecules is considered to involve with only H_N–I hydrogen bonds. The information about the most frequent and average bonding distances of H_N–I hydrogen bonds is used to explain atomistic mechanisms of molecular rotation in Sec. V B.

In this section, we focus on rotational motions of MA molecules. Subsection V A first describes statistical results obtained from our AIMD simulation. After the definition of rotational angles for the reorientation of a whole molecule and intramolecular rotation around the C–N bond, we apply the CSAA to the angles and then discuss preferable orientations of a C–N bond. Also, time constants regarding the two types of rotation are calculated to compare with the previous studies. In Subsection V B, we explain the possible atomistic mechanisms to describe the results of the CSAA using the information of H_N–I hydrogen bonds obtained in Section IV.

Here, we consider two types of MA molecular rotations, i.e., the reorientation of a whole molecule and intramolecular rotation around the C–N bond, and investigate their respective rotational behavior and preferred orientations using the CSAA.

As shown in Fig. 5(a), the spherical coordinates θ and ϕ represent the direction from C to N atoms. They are used to show the time series of whole molecule reorientation, where they are defined within the ranges of $0°≤θ≤$ 180° and $0°≤ϕ≤$ 360°. We use the angle θ to represent rotation around the C–N bond, where $0°≤ψ≤360°$ and the counterclockwise direction is positive. In addition, $ψN$ and $ψC$ represent the rotation angles in NH₃ and CH₃ groups, respectively. The detailed definition of ψ is described in Section 3 in the supplementary material.

Then, the CSAA is applied to our AIMD result in the same way of Frost et al.¹⁰ We first separate the surface area of the rotating sphere by whole molecule reorientation (whose origin is the position of the C atom) into equivalent 48 spherical triangles and folded these into one spherical triangle containing [100], [110], and [111] directions (referred to as 100-110-111 triangle). For example, the six spherical triangles in the first octant are shown in the left of Fig. 5(b) by solid lines on the sphere, and the 100-110-111 triangle is shown in the right. Subsequently, the 100-110-111 triangle is furthermore separated into 3 domains having similar surface areas each other using the center of gravity of the spherical triangle. We refer them to as {100}, {110}, and {111} domains,⁵⁹ where the area of the 110 domains is about 2.6% smaller than those of the {100} and {111} domains. The solid, dotted, and striped regions in the 100-110-111 triangle represent the {100}, {110}, and {111} domains, respectively. Also, it can be considered that these 3 domains consist of summations of different 6, 12, and 8 domains on the sphere. For the discussions below, we individually refer the domains to as (hkl) domain, where h, k, l ∈ {−1,0,1}, and also draw them by the solid, dotted, and striped regions. For example, the whole of striped region shown in the left of Fig. 5(b) represents (111) domain, which is one of the 8 domains contained in 111 domains.

Figure 5(c) shows the percentage of time that the C–N bond is parallel to each of the {100}, {110}, and {111} domains in the 19.8 ps simulation over all of the MA molecules to be 34.5%, 41.3%, and 24.2%, respectively, as shown by the solid bars. Namely, the times when the C–N bond is parallel to the $⟨110⟩$⁠, $⟨100⟩$⁠, and $⟨111⟩$ directions are the first, second, and third longest. The order coincides with the result of a recent NS experiment²³ of bulk structure in cubic phase performed at 350 K. However, the result is somewhat different from that by Frost et al.¹⁰ who reported that the percentages of time were 35%, 42%, and 23% for {100}, {110}, and {111} domains, respectively. Since their result was obtained using a pseudo-cubic structure, the origin of this discrepancy may be due to not considering the tilting motion in our simulation. If yes, our results could be applied only to film structures that are considered to be in cubic phase at room temperature. However, we should point out that the structure formed by I atoms in the system with a small number of atoms such as our $2×2×2$ cubic unit cell system (a total of 96 atoms) shows unusually large distortion, where the orientation of the C–N bond is fixed to $⟨100⟩$ directions due to the fact that H_N atoms form hydrogen bonds only with particular I atoms. Since 80 atoms were used in the study of Frost et al., the surrounding I atoms are likely to make MA molecules be parallel to $⟨100⟩$ directions. The detail is described in Section 2 in the supplementary material. Also, it is worth mentioning that a recent NS experiment by Ren et al.²³ showed the possible orientations of a C–N bond at room temperature estimated from the nuclear density, which were apparently inclined from the $⟨100⟩$ direction. Anyway, we will discuss the reason for our CSAA results from an atomistic point of view in Sec. V B.

To validate the reliability of our simulation, we compare the time constants involved with rotational motions experimentally obtained for MAPbI₃ bulk and film. First, we describe swinging motions of C–N bonds. The time series of the rotation angles θ, ϕ, and ψ of a selected MA molecule in 1 ps of the simulation are shown in Figs. 5(d) and 5(e), which were extracted from the original data shown in Fig. S5(a) of the supplementary material. For whole-molecule reorientation described by the time series of θ and ϕ in Fig. 5(c), ϕ vibrates around $80°$ until 450 fs while θ is nearly constant at about $90°$⁠. After 450 fs, ϕ becomes almost constant at about $180°$ and θ starts to vibrate around $80°$⁠. Using the separated domains defined above, swinging between the (011) and (111) domains is observed until 450 fs, and then the motion changes between the (⁠$1¯$00) and (⁠$1¯$01) domains. Such swinging events of the molecular C–N bond direction between neighboring domains in this way are a common feature of all MA molecules (see Fig. S5 of the supplementary material). When we evaluated the frequency of events on the 100-110-111 triangle, swinging is most likely to occur between the {100} and {110} domains. Hereafter, this event is referred to as {100} ⇌ {110}. The average number of the events per MA molecule in the 19.8 ps simulation is 225, i.e., a swinging event occurs every ∼0.1 ps, and the percentages of [100} ⇌ {110}, {110} ⇌ {111}, and {100} ⇌ {111} are 43.7%, 32.6%, and 23.7%, respectively. Thus, swinging events occur more frequently between the {100} and {110} domains, which reflects the first and second preferable orientations of $⟨110⟩$ and $⟨100⟩$ directions in our CSAA results. Bakulin et al.⁶⁰ demonstrated that two anisotropy decays of the excited-state population occurred in time scales of ∼0.3 ps and ∼3 ps in MAPbI₃ film using ultrafast 2D-IR vibrational anisotropy spectroscopy. By performing AIMD simulation, they also suggested that the faster and slower decays involve swinging motions of C–N bonds and large reorientations of MA molecules, respectively. Therefore, the frequency of swinging events estimated to be ∼0.1 ps is comparable to the faster anisotropy decay time scale. On the other hand, to compare with the time scale of slower anisotropy decay, we calculate time dependent self-correlation function C(t) of the C–N bonds,²⁷ 

where $b→$ is a unit vector pointing from the C atom toward the N atom of the same MA molecule and $⟨⟩$ denotes time average. The result is shown in Fig. 5(f), and the corresponding correlation time, estimated from the exponential decay, is ∼3.6 ps. Although the rotation of MA molecules in our system might have been affected by finite size effects as suggested by Carignano et al.,²⁷ it is in reasonable agreement with the time scale of the slower decay. In addition, our result of time dependent self-correlation function is comparable to the time constants obtained in bulk structure, and the dipole rotational relaxation time of 5.37 ps in MAPbI₃ bulk from millimetre-wave spectroscopy by Poglitsch and Weber⁶¹ is also the same order of magnitude. The orientation residence time of 13.8 ps from quasi-NS experiment by Leguy et al.²² is one-order of magnitude longer than that of our result. However, it is worth noting that it coincides with the time scale that the self-correlation function reaches zero (∼12 ps).

Second, we discuss rotation around the C–N bond. We found that the average numbers of rotation for NH₃ and CH₃ groups are 4.77 and 5.33 within 19.8 ps, i.e., the times for one rotation are 4.15 and 3.71 ps, respectively. These results agree quite well with the corresponding residence time of 3.1 ps for both groups measured in bulk structure by Leguy et al.²² Interestingly, however, we observed much faster rotation than these time scales. As shown in Fig. 5(d), until 450 fs, $ψN$ shows rapid rotation, as shown by the solid part of the curve. We extracted the periods (the solid parts of the curves in Fig. 5(d)) by imposing a lower limit on the rotational speed $ψ˙$⁠. Details of the extraction are described in Section 5 in the supplementary material. The time required for one rotation is approximately 280 fs, which is much shorter than the experimental correlation time for rotation around the C–N bond of 460 fs²¹ and the average rotation times mentioned above. Such rapid rotation occurs for both of the NH₃ and CH₃ groups on average for 1.2 ps within 19.8 ps simulation and may influence the surrounding polarized structure because of a lot of H–I hydrogen bond switching events accompanied by the rotation. The time series of ψ for selected MA molecules with and almost without continuous rotation are shown in Fig. S5 of the supplementary material. On the other hand, subsequent rotation becomes slower (the dashed part of the curve). The change in rotational behavior is attributed to variation in the C–N bond after about 450 fs, as shown in Fig. 5(d), that is, rapid rotation probably depends on the orientation of the C–N bond. It is worth noting that the percentages of time of rapid rotation for $ψN$ in the {100}, {110}, and {111} domains on the 100-110-111 triangle are calculated using the CSAA to be 31.3%, 46.6%, and 22.1%, respectively, as shown by the striped bars in Fig. 5(c). Meanwhile, the percentages for $ψC$ are 28.9%, 46.7%, and 24.4%, respectively, which are almost the same as those of $ψN$⁠. As a result, we conclude that rapid rotation around the $⟨110⟩$ directions easily occurs on both sides of the MA molecule.

Although it seems that the time series of $ψC$ is different from that of $ψN$ (Fig. 5(d)), they are actually very similar (see Fig. S5 of the supplementary material) because the activation energies for independent rotation for both groups have been estimated to both be about 2 kcal/mol⁶² and it is difficult to overcome this energy barrier at 300 K. Since the H_C atom almost does not form hydrogen bonds with I atoms, it is considered that $ψC$ is followed by $ψN$⁠.

Here, we try to explain the atomistic mechanisms for the results of CSAA, which the $⟨110⟩$ directions are preferred for both of the two rotations, from the viewpoint of the relation between H_N–I hydrogen bonds and the orientation of an MA molecule. Figure 6 shows snapshots of the atomic configurations of the MA molecule whose time series of θ, ϕ, and ψ were discussed in Section V A. The configurations at 23, 716, and 217 fs with the surrounding PbI₃ cage are shown, which are when the C–N axes are almost parallel to the [011], [$1¯$00], and [111] directions, respectively (see Fig. 5(c)). Actual spherical coordinates (θ, ϕ) are (44°, 98°), (90°, 176°), and (38°, 52°). These are used as selected examples of MA molecules parallel to the $⟨100⟩$⁠, $⟨110⟩$⁠, and $⟨111⟩$ directions. Analysis based on such limited atomic configurations may not lead to comprehensive conclusion regarding rotational dynamics, but it certainly provides valuable insights into its mechanisms. For the discussions, we also use the corresponding optimized structures of the MA molecule in the fixed PbI₃ cages shown in Fig. S10 in the supplementary material. Note that we describe all of the switching events of hydrogen bonds within 1 ps of the simulation with additional snapshots and population analysis results in Section 6 in the supplementary material, providing the rotation process in great detail. A supplementary movie is also available. The labels in Fig. 6 are based on the order in the supplementary material.

First, the relation of the $⟨110⟩$ directions with H_N–I hydrogen bonds is discussed. When the C–N bond is parallel to the [011] direction, the second-neighbor I atoms I4 and I10 are close to the H_N atoms, as well as the first-neighbor I atoms I1, I2, I3, and I5. Here, the concepts of first- or second-neighbor I atoms are assumed for a perfect Pb–I lattice. The interatomic distances shown in Fig. 6(a) are within or close to the range of hydrogen bond distance of H_N–I of 2.2-3.7 Å. In addition, there exist short, middle, and long hydrogen bonds, among which only the short bond is a strong hydrogen bond because its distance is shorter than the most frequent H_N–I hydrogen bonding distance of 2.9 Å (see Fig. 4). This means that the H_N atoms easily form hydrogen bonds with the six I atoms arranged in a hexagon on the plain through the three H_N atoms, without being strongly bound to multiple I atoms. It is therefore considered that the geometry of the surrounding I atoms and their suitable interatomic distances to form hydrogen bonds help to maintain the direction of the C–N bond and induce rotation around the C–N bond. Note that spherical coordinates in the optimized MA molecule shown in Fig. S10(a) of the supplementary material are almost unchanged (θ = 42°, ϕ = 93°), which is maintained by three H_N–I hydrogen bonds with three of the six I atoms. Similar to the configuration before optimization, these have different bond distances and only one bond is shorter than 2.9 Å. The vibration of ϕ until 450 fs is induced due to that the H_N atoms switch the binding partners from among the six I atoms. The vibration process accompanied by rotation around the C–N bond is described in Section 6.1 in the supplementary material. On the other hand, the large reorientation after 450 fs can be explained from the formation of hydrogen bonds between the I atom labeled I8 and the H_N atoms because I8 is also close to the H_N atoms (interatomic distances are 3.6–3.7 Å). This process is described in Sections 6.1 and 6.2 in the supplementary material.

Second, the relation of the $⟨100⟩$ directions with H_N–I hydrogen bonds is discussed. When the C–N bond is parallel to the [$1¯$00] direction, as shown in Fig. 6(b), the interatomic distances between the H_N atoms and the first-neighbor I atoms I1, I2, I6, and I7 are relatively short (2.2-3.7 Å), whereas the interatomic distances between H_N atoms and the second-neighbor I atoms I4, I8, I10, and I12 are relatively long, and some are greater than 4 Å. This means that H_N atoms mainly form hydrogen bonds with the first-neighbor I atoms arranged in a square on the plain through the three H_N atoms. In this situation, two of the H_N atoms form short hydrogen bonds and the remaining one forms a long hydrogen bond. Since it is difficult for such an asymmetric bonding state to maintain the direction of C–N bond, it would be immediately inclined from the [$1¯$00] direction. In fact, the C–N bond tilts from the [$1¯$00] to the [$1¯$01] directions after 716 fs in our simulation. The spherical coordinates of the optimized structure also show a large change in θ (θ = 69°, ϕ = 179°). Meanwhile, the reason why the $⟨100⟩$ directions are the second preferable is that the interatomic distances between H_N atoms and first-neighbor I atoms would be too close. There are at least two shorter bonds than the most frequent bonding distance of 2.9 Å in both of the configuration shown in Fig. 6(b) and its optimized one. Hence, the molecules hardly go outside the influence range of the four I atoms. In addition, the situation would induce the swinging motion of C–N bonds as seen, e.g., for θ at t > 500 fs in Fig. 5(c), the process of which is described in Section 6.3 in the supplementary material.

Finally, the relation of the $⟨111⟩$ directions with H_N–I hydrogen bonds is discussed. Although the C–N bond in Fig. 6(c) is nearly parallel to the [111] direction, the interatomic distances are much longer than those for the [$1¯$00] and [011] directions. The minimum distance is 2.94 Å. Since there is no bond shorter than the most frequent bonding distance of 2.9 Å, it would be immediately inclined from the [111] direction. In fact, the spherical coordinates of the optimized structure show a large change in ϕ (θ = 30°, ϕ = 71°). Lack of driving force for maintaining the direction of the C–N bond would be reflected in the results of our CSAA.

In summary, we have performed AIMD simulations using a cubic-MAPbI₃ model at room temperature to investigate how the MA molecule rotates in the inorganic PbI₃ cage through interaction between H atoms of the molecules and surrounding I atoms from the atomistic viewpoint as well as statistical analyses. Our motivation was also to make results of this study to be helpful in the understanding of the molecular rotational motions not only in MAPbI₃ film, which is considered to be in cubic phase at around room temperature, but also in bulk structure. While we found that our total RDFs agree well with both of those measured in bulk and film structures, our result of the CSAA to clarify the preferable orientation of C–N bonds was slightly different from that of the previous AIMD simulation by Frost et al. The $⟨110⟩$ directions in the PbI₃ cage would be preferred in our study, but the $⟨100⟩$ directions would be done in the previous study. Since the discrepancy may be attributed to the neglect of tilting motions in our cubic system, our CSAA result could be applied only to film structures. However, the time constants regarding the rotational motions statistically obtained from our simulation were in reasonably accordance with the corresponding constants measured in bulk and film structures. Thus, we consider that our simulation using the cubic structure could reproduce most of the rotational behaviors in the bulk structure. Meanwhile, we found much faster rotations around C–N bonds than the average rotation time, and revealed using the CSAA that they preferably occur when the C–N bonds are parallel to the $⟨110⟩$ directions as well. Furthermore, we have tried to explain the atomistic mechanisms leading to our CSAA results by focusing on bonding states of H_N–I hydrogen bonds around MA molecules. Here, we defined the hydrogen bonds using charge density between the H_N and I atoms obtained by population analysis and found that H_N atoms always form a hydrogen bond with an I atom with a wide range of interatomic distances (2.2–3.7 Å). In addition, the most frequent bonding distance was estimated to be 2.9 Å. From the results, we gave a possible explanation of the reason why the $⟨110⟩$ directions are preferred, that is, the atomic arrangement and distances of the surrounding I atoms are more suitable to form H_N–I hydrogen bonds than those of other directions. We believe that such atomistic-level findings and the statistical results of the CSAA and population analysis provide a deeper understanding of the rotational behavior of the MA molecules and the formation of polarized local domains in MAPbI₃ bulk and film. We also hope them to be helpful in the fabrication of efficient photovoltaic devices such as perovskite solar cells.

See supplementary material for a detailed description and a movie of rotational dynamics of the selected MA molecule, including the definition of ψ and its time derivative $ψ˙$⁠. The unusual distortion of I atoms in $2×2×2$ cubic unit-cell system is also described. In addition, the optimized atomic configurations of those in Fig. 6 are shown.

This work was partly supported by a Grant-in-Aid for JSPS research fellows (No. 15J01350). The work at USC was supported by the US Department of Energy, Office of Science, Basic Energy Sciences, Materials Science and Engineering Division, Grant No. DE-FG02-04ER-46130. Computational support for this work was provided by the Research Institute for Information Technology, Kyushu University. The authors also kindly acknowledge the Supercomputer Center, Institute for Solid State Physics, the University of Tokyo, for the use of its facilities.