Hydrogen Bonds in Lead Halide Perovskites: Insights from Ab Initio Molecular Dynamics

Hydrogen bonds (HBs) play an important role in the rotational dynamics of organic cations in hybrid organic/inorganic halide perovskites, thus affecting the structural and electronic properties of the perovskites. However, the properties and even the existence of HBs in these perovskites are not well established. In this study, we investigate HBs in perovskites MAPbBr3 (MA+ = CH3NH3+), FAPbI3 (FA+ = CH(NH2)2+), and their solid solution with composition (FAPbI3)7/8(MAPbBr3)1/8, using ab initio molecular dynamics and electronic structure calculations. We consider HBs donated by X-H fragments (X = N and C) of the organic cations and accepted by the halides (Y = Br and I) and characterize their properties based on pair distribution functions and on a combined distribution function of the hydrogen–acceptor distance with the donor–hydrogen–acceptor angle. By analyzing these functions, we establish geometrical criteria for HB existence based on the hydrogen–acceptor (H–Y) distance and donor–hydrogen–acceptor angle (X–H–Y). The distance condition is defined as d(H – Y) < 3 Å for N–H-donated HBs and d(H – Y) < 4 Å for C–H-donated HBs. The angular condition is 135° < (X – H – Y) < 180° for both types of HBs. A HB is considered to be formed when both angular and distance conditions are simultaneously satisfied. At the simulated temperature (350 K), the HBs dynamically break and form. We compute the time correlation functions of HB existence and HB lifetimes, which range between 0.1 and 0.3 ps at that temperature. The analysis of HB lifetimes indicates that N–H–Br bonds are relatively stronger than N–H–I bonds, while C–H–Y bonds are weaker, with a minimal influence from the halide and cation. To evaluate the impact of HBs on the vibrational spectra, we present the power spectrum in the region of N–H and C–H stretching modes, comparing them with the normal mode frequencies of isolated cations. We show that the peaks associated with N–H stretching modes in perovskites are redshifted and asymmetrically deformed, while the C–H peaks do not exhibit these effects.


INTRODUCTION
Hybrid organic−inorganic halide perovskites (HOIHPs) constitute a versatile family of materials that have attracted much research attention during the last decade. Since first proposed in 2009 by Kojima et al., 1 perovskite solar cells have emerged as a strongly promising candidate for efficient and cheap solar cells and have attained an extraordinary 26.0% record photoconversion efficiency (PCE) as single junction solar cells and 33.7% in tandem perovskite/silicon solar cells. 2 Other promising applications of HOIHPs include X-ray detectors, 3 LED devices, 4 and water splitting photocatalysts. 5 HOIHPs are flexible in terms of morphology and can be synthetized as bulk crystals, bidimensional laminar crystals, or quantum dots.
The crystal structure of 3D halide perovskites, with the general formula ABY 3 , is well illustrated by methylammonium lead iodide (CH 3 NH 3 PbI 3 ), represented in Figure 1a. The Asite of perovskite is occupied by the organic cation , the B-site is occupied by the lead cation (Pb 2+ , gray balls), and the Y-sites are occupied by the iodide anion (I − , violet balls). Methylammonium is often abbreviated as MA or MA + , leading to the short names MAPbI 3 or MAPI for the compound CH 3 NH 3 PbI 3 . Other members of the family of 3D HOIHPs can be generated, replacing the MA by another organic cations like formamidinium (FA), CH(NH 2 ) 2 ( Figure 1b shows the structure of FAPbI 3 ), or by a large inorganic cation like cesium; the iodide anion can be also replaced by bromide or chloride, and the lead cation can be replaced by tin or germanium. The full family includes solid solutions of the pure compounds, where each crystal site can be occupied randomly by one or several of the above-described components. These mixed compounds show a remarkable improvement in their stability (the HOIHP's Achilles heel) and are related with the evolution of record PCE. 6 The perovskite structure is generally stable when the ionic radii r A , r B , r Y satisfy an empirical rule 0.8 < t < 1.0, according to the Goldschmidt tolerance factor: 9 t r r r r 2 ( ) The nature of the A-site cation has a dramatic influence on the phase stability and on crystal symmetry. 10 Substitution of the A-site cation by organic cations larger than FA, with t > 1.0, limits crystal growth and allows the formation of laminar structures (2D perovskites), 11 which are periodic in two directions and are a few unit cells wide in the third direction. Another example is CsPbI 3 , which is stable in the perovskite structure at temperatures just over 600 K, but when cooled to room temperature, it makes a fast transition to a yellow nonperovskite phase. A similar behavior is exhibited by FAPbI 3 ; however, their solid solution, Cs x FA 1−x PbI 3 , is stable at room temperature. 6 Perovskite phases show a diversity of space groups and phase transitions at different temperatures, depending on the organic cation, 10,12 as well as on the B-site cation 10 and on the halide. 13 As Figure 1 illustrates, the inorganic PbY 3 backbone presents large distortions from the ideal perovskite structure. The cells shown in this figure represent formula units of the HOIHP but do not possess the crystal symmetries inferred from diffraction techniques. Time-and space-dependent realizations of these basic cells conform the average structures that make the apparent crystal periodicity in diffraction experiments. Highertemperature perovskite phases display an apparent cubic symmetry, where the crystallographic unit cell coincides with those of Figure 1, but that represents only the average positions of the atoms, with large thermal ellipsoids. Moreover, the organic cations rotate and make orientational jumps across all the possible equivalent orientations. For the lowertemperature perovskite phases, the primitive cells are supercells of those represented in Figure 1. For example, MAPbI 3 undergoes a cubic to tetragonal transition at ∼330 K, and the primitive cell has two formula units, and a tetragonal to orthorhombic transition at ∼160 K with a primitive cell of four formula units. In the latter case, the cations do not rotate, and all atoms keep at definite positions.
Hydrogen bonds (HBs) are possible in HOIHPs between the halides Y − and the protons bound to nitrogen or carbon atoms, as indicated by dashed lines in Figure 1. HBs may influence the rotational dynamics of the organic cations and therefore play a role in the stabilization of the HOIHP, the deformations of the inorganic backbone, and indirectly on the electronic structure. 9,14−23 The understanding and characterization of HBs in perovskites have therefore attracted research attention. Figure 2 illustrates the structure of FA and MA. Methylammonium is a result of the protonation of methylamine neutral molecule (H 2 N−CH 3 ), where the proton is attracted by the lone pair of the H 2 N group. The protonated ammonium group (H 3 N) has three equivalent N−H bonds; therefore, the positive charge 24 is considered a property of the full group rather than localized at any of the hydrogen atoms. Hence, the hydrogens of the ammonium extreme are selectively attracted to the halide anions in the perovskites. FA is a result of formamidine (H 2 N−CH=NH) protonation. In formamidine, the central carbon atom is involved in a single bond with a nitrogen atom and a double bond with the other nitrogen atom, leaving a lone pair at the latter. However, upon protonation at the NH group, the structure becomes resonant between the Lewis structures H 2 N−CH=NH 2 and H 2 N=CH− NH 2 . FA has a planar structure, as shown in MD simulations, which can be attributed to the double-bond character between the carbon and both nitrogen atoms. The delocalization of the double bond toward the incoming proton can be understood as a positive charge density distributed along the N−C−N backbone, as shown in Figure 2. Hence, the attraction exerted by halide anions upon each hydrogen should be smaller than in the MA case due to this distribution of the positive charge. We consider first the HBs donated by the NH 2 groups, but later we will also consider the possibility of donation by the CH group.
Svane et al. 25 devised a method to compute the HB energies in hybrid perovskites, filtering out the electrostatic interaction of the charged organic cation with the inorganic lattice. For halide perovskites, they obtained energies in the range of 0.02 to 0.27 eV per organic cation, which is one order of magnitude smaller than in formate perovskites. In the high-temperature phases of HOIHPs, the HBs are not permanent but form and break dynamically, allowing cation rotation and determining the symmetry of the crystal. The existence of HBs in leadbased HOIHPs at room temperature has been recently questioned by Ibaceta-Janã et al. 26 based on a combined Raman and density functional theory (DFT) study.
The question of the existence of HBs at given conditions requires a careful definition of these bonds. According to IUPAC recommendations, 27,28 HBs should display some charge transfer, show directional preference, and show a bond path connecting the hydrogen with the HB acceptor. Early conceptions of HBs considered only donors like F−H, O−H, and N−H, the last one present in HOIHPs. However, the list of donors has expanded to include any molecule or molecular fragment X−H, where X is any element with electronegativity larger than H, i.e., F, N, O, C, P, S, Cl, Se, Br, and I. 28 The acceptor Y can be any of these elements. Therefore, the HB donors in HOIHPs can be the N−H and C−H groups of the organic cations, while the acceptors can be the halides Cl, Br, and I. Directionality is an important criterion; the angle X−H−Y should be close to 180°. Vibrational spectroscopy provides a signature of HBs. The X−H stretching modes typically decrease its frequency (although a blue shift is also possible), accompanied by an increase of bandwidth and change of intensity in the IR spectra. The X−H peak can even disappear from the spectrum, although this is understood as a large redshift and mixing with other modes. 28 The HBs also modify the topology of electron density, obtained either from experiments or from electronic structure calculations. In particular, the electron density must present a (3, −1) critical point along the H−Y direction, meaning that the gradient of the density is maximum (minimum) along two (one) orthogonal directions. Methods based on the electron density topology are available to identify HBs using quantum chemistry calculations. 21,29,30 Finally, HBs need to be thermally stable to have practical significance, but the absence of effects does not deny its very existence as temporary states.
In this work, we present a systematic theoretical study of the HB features for three HOIHPs, MAPbBr 3 , FAPbI 3 , and the solid solution (FAPbI 3 ) 7/8 (MAPbBr 3 ) 1/8 , using ab initio molecular dynamics (AIMD). For these compounds, we will use the short names MAPBr, FAPI, and MAFA, respectively. We will present statistical functions related to the HBs, such as radial distribution functions, combined (radial and angle) distribution functions, autocorrelation functions, and lifetimes, to address some of the open questions about HBs in halide perovskites. This work builds upon the model of MAFA solid solution presented in ref 31 and the AIMD simulations developed there for MAFA, MAPBr, and FAPI. In the present study of HBs, we have imported methods and codes from the field of ionic liquids, explaining the methodology and the numerical parameters (cutoff distances and angles) that are appropriate to define the HBs in HOIHPs.
One important precedent is ref 20, where the HBs in pure MAPBr were studied using classical molecular dynamics (MD), i.e., with parameterized force fields. The latter method, although generally less accurate than AIMD, allowed the authors to include a large number of atoms in the simulation and perform an extensive statistical study, focusing on properties different from those considered in our study, such as HB energetics, spatial correlations of cation orientation, and its influence on the electronic structure. Notably, the directionality of the HB 27,28 was not considered in ref 20, which led them to overestimate the HB population and probably its importance.
In contrast, Ibaceta-Janã et al. 26 claimed that no HBs exist at room temperature in HOIHPs, or at least they cause no observable effects in the infrared and Raman spectra. To shed some light in this issue, we also present the calculations of the power spectra that show that HBs may have observable effects in the region of highest frequencies, not covered by the experiments of ref 26. We also present the calculations of the reduced density gradient that reveal HB effects on the electronic structure of certain temporary states.

METHODS
The (FAPbI 3 ) 7/8 (MAPbBr 3 ) 1/8 solid solution is represented by a special quasi-random structure, with an ion distribution that mimics the random distribution in terms of short-range pair correlation functions, as described in previous work. 31 The coordinate files of selected configurations are available in ref 31, while the full MD ensemble is deposited in a free repository. 32 Our AIMD simulations at constant volume and temperature (NVT ensemble), with T = 350 K, were carried out using the CP2K package. 33 The ionic forces were calculated using DFT in the generalized gradient approximation as implemented in the Perdew−Burke−Ernzerhof (PBE) functional, 34 with the Grimme correction scheme DFT-D3 35 to account for the dispersion interactions. The hybrid Gaussian and plane wave method 36 implemented in the QUICKSTEP module of the CP2K package was used. The Kohn−Sham orbitals of valence electrons are expanded in a Gaussian basis set (DZVP-MOLOPT for Pb, I, Br, C, N, and H). 37 Core electrons were treated using dual-space GTH pseudopotentials. 38−40 Further details of the AIMD calculations can be found in ref 31.
The TRAVIS code 41 was used to characterize the HBs via pair distribution functions (PDFs), combined distribution functions (CDFs), and time correlation functions. 41 The PDF, also called the radial distribution function, gives the probability that two atoms of different species are separated a given distance, relative to the probability in a non-interacting gas; the distance of separation is the argument of the function. The PDF reveals the short-range order of a liquid or a disordered material; bond lengths are identified from the peaks at the shortest distances, and coordination spheres are identified from the minima of the PDF. TRAVIS also allows us to obtain other distributions like angle distribution functions and dihedral distribution functions. The CDFs are combinations of the former distributions, providing the probability of finding certain combinations of geometrical parameters, such as the pair distance and the angle between three atoms. In particular, the CDF of the halide−hydrogen distance and the halide− hydrogen−nitrogen angle allows to reveal the HBs in HOIHPs. The HB dynamics can be characterized by the time correlation functions of the pair forming the HB: where β ij (t) = 1 if there is a HB between atoms i and j (one halide and one hydrogen) or zero otherwise. The condition of existence of a HB is defined in terms of a distance and an angle, as will be discussed below. The function βĩ j (t + τ) = 1 if the HB exists for all the time between t and t + τ, and is zero if the HB breaks at any instant. This is called the continuous autocorrelation function in TRAVIS. The function C C HB (τ) decays to zero, and the lifetime of the HB can be obtained as The Journal of Physical Chemistry C pubs.acs.org/JPCC Article TRAVIS also allows us to compute intermittent autocorrelation functions, where βĩ j (t + τ) is replaced by β ij (t + τ), which equals unity if the HB exists at t + τ, regardless of whether it breaks or not at intermediate times. The intermittent autocorrelation function is useful for HBs in liquids but not in HOIHPs because the organic cations are confined in a cavity and any broken HB will form again and again, not allowing the autocorrelation function to decay to zero. However, this kind of function has been used in ref 23, allowing to extract approximate lifetimes from the initial fast decay.

HB Donated by Ammonium and Amidine
Groups. Figure 3a Figure 3c,d shows that halide−nitrogen distances are smaller than halide− carbon distances. The finding that the Y−N and Y−H−N distances are shorter than the Y−C and Y−H−C distances, respectively, has been reported for MAPI elsewhere. 15,17 This behavior can be expected from the positive charge of the NH 3 group in MA. This argument does not apply to FA, as the positive charge is delocalized along the N−C−N backbone. Instead, the elongated shape of FA favors the N atoms at the edges getting closer to the halides. Other remarkable property is that the Y−C PDFs have minima with null or almost null value at ∼6 Å, coincident with the lattice parameters (6.3620, 6.3115, and 5.9328 Å for FAPI, MAFA, and MAPBr, respectively). The corresponding minima of Y−N PDFs are not that sharp but are still well defined. In contrast, the first minima Y−H−X PDFs are still higher than 0.5. This is a consequence of H atoms being distributed in large rigid cations. However, as we will see next, sharp minima are obtained when directionality criteria are incorporated.
HBs are better resolved by means of CDFs, which are functions of the Y−H distance and the Y−H−N angle. Figure  4 shows the CDF for the same combinations as in Figure 3a.
The maximal values take place in the region of the red spots, which can be approximately defined by the simultaneous conditions d < 3 Åand 135°< (I − H − N) < 180°. Even for I−H−N(FA) in FAPI and MAFA, which do not show a peak in the PDF, the CDF shows a clear maximum at distances between 2 and 3 Å, followed by minimum between 4 and 5 Å, provided that the angle is larger than 135°. Therefore, we take this couple of geometrical conditions as the definition for the existence of a HB.
The animation videos V1 to V6 in the Supporting Information show that, with these conditions, HBs form and break dynamically during the simulation time. These animations were generated with the Visual Molecular Dynamics (VMD) graphics program, 42 where HBs are defined somewhat differently. HBs in VMD are defined by the Y−N distance instead of the Y−H distance and the same condition for the angle (Y − H − N). For these videos, the Y−N distance cutoff was set as 4 Å, i.e., 1 Å larger than the Y−H distance cutoff that we have explained above. It is tempting to derive this condition from the CDF of the Y−N distance and Y−N−H angle. This CDF, shown for FAPI and MAPBr in Figure 5a,c, shows the maxima for distances between 3 and 4 Å, but it is almost independent of the angular condition. The snapshots in Figure 5b,d show the Y−N distances smaller than 4 Å (sticks in gray color). Using just the distance condition would imply HBs with up to three halide atoms simultaneously. However, HBs defined by the combination of distance The animation video V7 in the Supporting Information complements this view, indicating that the first peaks of the Y− N PDF seen in Figure 3b,c are not explained by HB formation. Proximity of nitrogen and iodine in FAPI has instead steric origin due to the cation size and the location of N atoms at opposite sites. In the case of MA, the peak of the Br−N PDF also has only a partial contribution from HBs. The smaller cation size and the presence of only one N atom suggest that this peak is caused by the electrostatic attraction between the halide and the NH 3 + group and HB. However, the region in the CDF of Figure 5a,c, corresponding to (Y − H − N) > 135°, shows a depletion for the Y−N distance between 5 and 6 Å. This is the signature of HB in this CDF, explaining why the HB can be visualized in VMD by defining the HB through the angular condition together with the Y−N distance, although VMD displays a dotted line along the H−Y path. The same behavior can be appreciated for the solid solution MAFA in Figure S1 of the Supporting Information.
After the geometrical definition of HBs has been established, we now probe their dynamics by means of autocorrelation functions, as defined in eq 1, and the lifetimes, as defined from eq 2. Based on the autocorrelation functions ( Figure 6), we present the HB lifetimes in Table 1. The lifetimes are dependent on the limits set to the distance and angle, but given a common definition, we can observe several trends.     Table S2 of the Supporting Information. The HB with iodine has shorter lifetimes than the HB with bromine. Table 1 also includes the lifetimes of HBs with hydrogens attached to carbon atoms, which will be discussed in the next section.

HB through Carbon Atoms.
We also explore the possibility of the HB donated by the C−H groups of MA and FA. Deprotonated FA can be formamidine (H 2 N−CH=NH) or diaminocarbene (H 2 N−C−NH 2 ). The latter is highly unstable, but it has been found in mass spectroscopy experiments. 43 Regarding HOIHPs, several authors 21,44 have studied the role of Y−H−C bonds, among several noncovalent interactions, and they have pointed the importance of HBs in determining the extent of the PbI 6 octahedra tilting in the low-temperature phases. Hence, we proceed to study this possible HB with the same approach used for Y−H−N bonds.
The CDF of the I−H distance with the I−H−C angle, shown for FAPI in Figure 7a Had we set the same cutoff distance for both types of HBs, the Y−H−C HBs would have shorter lifetimes. A broader question is whether these HBs exist at all. The next section presents additional arguments.

Impact of HBs on the Power Spectrum.
It is well known that HBs modify the frequencies of N−H and C−H stretching vibrations 27,28 and the shape of the power spectrum, 45 providing experimental evidence through the IR spectra 46,47 or Raman spectrum. 20 Hence, we evaluate the associated power spectrum in FAPI and MAPBr. The power spectrum of MAFA has been shown to be well approximated by a weighted average of the spectra of pure compounds. 31 Figure 8 shows the power spectrum of FAPI (a) and MAPBr (b) in the range of 3000−3700 cm −1 , compared with the vibrational density of states (VDOS) of isolated FA and MA. The frequencies of the normal modes at the minimal energy geometry were obtained by diagonalizing the Hessian of the energy as a function of the internal coordinates. The normal mode frequencies and vibration patterns are shown in Tables S3 and S4 in the Supporting Information. This method provides the power spectrum in the harmonic approximation at zero temperature. The resulting delta-like VDOS has been broadened by means of a Lorentzian function with a full width at half-maximum of 16 cm −1 for each vibrational mode. The peak at 3162 cm −1 in Figure 8a comes from the C−H  The Journal of Physical Chemistry C pubs.acs.org/JPCC Article stretching mode, and the peaks at higher frequencies come from N−H stretching modes. Similarly, for MA, the peak at 3038 cm −1 in Figure 8b comes from the symmetric stretching of C−H bonds in CH 3 , while the peak at 3152 cm −1 comes from two asymmetric stretching modes of C−H bonds. The two peaks at higher frequencies also come from N−H stretching modes. The power spectra obtained from MD simulations by means of the Fourier transform of the velocity autocorrelation function are shown in red and blue lines for both FAPI and MAPBr. The blue lines include a frequency (wavenumber) correction intended to cancel a systematic error in the oscillator frequencies introduced by the finite time step in the Verlet algorithm. 41,48 Considering the corrected curves, the vibrational frequencies at the temperature of the MD simulation (350 K) are redshifted from the zero-temperature frequencies in all cases (even in the case of the first peak in Figure 8b, there is a small redshift of 1 cm −1 ). This redshift can be understood as caused by anharmonicity and the environment. Note that peaks of the N−H stretching modes undergo larger redshifts than the C−H stretching mode peaks but also get broadened asymmetrically toward the low-frequency side. This asymmetry of the N−H stretching mode peak is a signature of the HB. 45 The peaks of C−H stretching modes do not show asymmetric broadening, as well as none of the other peaks at frequencies under 1800 cm −1 (see Figure S2 in the Supporting Information). Another interesting aspect is that the power spectrum of MAPBr shown in Figure 8b displays only two peaks and one broad shoulder, while the harmonic spectrum displays four peaks. Figure S3b shows the deconvolution of this spectrum in four components, showing an excellent fit. We conclude that the two highest frequency peaks, those corresponding to the N−H stretching modes of isolated MA, get strongly broadened and redshifted, overlapping with the peaks from C−H stretching modes. A similar behavior can be appreciated, although less dramatically, in Figure 8a for FAPI, whose power spectrum can be fitted with five components ( Figure S3a). Furthermore, the lowest frequency peaks in Figure 8, corresponding to C−H stretching modes, are well fitted with Lorentzian functions, while the fitting functions centered at higher frequencies are Gaussians. All this supports the interpretation that C−H stretching modes of FA and MA are rather well preserved in the perovskites, but the N−H stretching modes are strongly modified. We conclude that this is an observable effect of HBs donated by the N−H groups in both FA and MA. Our analysis also suggests that the C−H−Y HBs, if they exist at all, have little effect on vibrational properties.

Electronic Structure.
HBs and other non-covalent interactions (NCIs) can be proved in electronic structure calculations by means of the reduced density gradient (RDG). 29 Varadwaj et al. 21 have shown that the RDG-NCI analysis supports the existence of HBs in MAPI. Here, we show that their results also apply to FAPI and MAPBr. The RDG is defined as s = 1/(2(3π 2 ) 1/3 )| ∇ ρ|/ρ 4/3 , where ρ is the electron density. As commented above, one of the IUPAC criteria for HBs is the existence of a (3, −1) critical point along the Y−H path. In a (3, −1) point, the RDG is null, and the Hessian matrix of the density has two negative eigenvalues and one positive eigenvalue, ordered as λ 1 ≤ λ 2 ≤ λ 3 . Hence, the negative sign of λ 2 (as for λ 1 ≤ λ 2 ) where s = 0 indicates a (3, −1) critical point. Low values of the density indicate NCIs, in contrast with covalent bonds associated with larger values of the electron density. Hence, there is a certain range of values of the pair (s, sign (λ 2 )ρ) that is typical in the vicinity of the (3, −1) critical point of HBs. This range is approximately s ∈ (0.2,0.5) and sign(λ 2 )ρ ∈ ( − 0.05, − 0.01). 21,29 Figure 9 shows the RDG isosurfaces in unit cells of FAPI (a) and MAPBr (b). The isosurfaces s(x, y, z) = 0.5 are shown, using a color scale graded from the function sign(λ 2 )ρ. The calculations were made using Gaussian09 49 for non-periodic clusters FA 19 Pb 8 I 36 and MA 19 Pb 8 Br 36 , both with charge −1 to have a closed shell electronic structure. The PBE functional was used with the LANL2DZ basis set. The cluster of FAPI was cut from a configuration of the MD ensemble that has one HB minimal I−H(N) distance (2.27 Å), while the cluster of MAPBr was cut from a relaxed periodic structure with a single unit cell, having two HBs with the Br−H(N) distance equal to 2.43 Å and (Br − H − N) = 168°. The regions shown in Figure 9 belong to the central unit cell of each cluster. For both FAPI and MAPBr, we can appreciate a small green cloud (sign(λ 2 )ρ ≈ − 0.04) in the middle of each HB (green broken lines), but no similar cloud can be seen near the other hydrogens that do not fulfill the HB geometrical conditions. It is particularly interesting that MAPBr does not show any (

DISCUSSION AND CONCLUSIONS
Let us now discuss some recent publications about HBs in HOIHPs in the context of our results. Saleh et al. 20 have conducted a classical MD study with a huge number of atoms. They concluded that the HBs control the energetics of MA orientations. As they did not consider the directional criterion 27 in their geometrical definition of the HB, they obtained HB populations that are much greater than our estimations. On the opposite side, Ibaceta-Janã et al. 26 have recently disputed the presence of HBs in halide perovskites at room temperature. They proved HBs using Raman spectroscopy, comparing signatures of HBs in water, MAI, and FAI with the spectra on MAPbY 3 and FAPbY 3 . The absence of spectroscopic evidence does not prove the inexistence of HBs but rather proves its irrelevance. However, the failure to observe HB signatures in their Raman spectra can be attributed to the fact that they centered the analysis in the range of 500− 1800 cm −1 , where our simulations also do not show appreciable effects. According to our simulations (Figure 8) 50 have compiled a long list of activation barriers and residence times associated to the movement of organic cations in different HOIHPs. Residence times associated to wobble-in-a-cone motion of MA could be related with the HB lifetimes, although an exact correspondence cannot be established. These residence times at room temperature (∼0.3 ps) are of the same order of magnitude as our HB lifetimes. These characteristic times show up in rotational time correlation functions obtained from MD simulations. 51 The difference with the lifetime of the HBs can be that the latter can break and form again before a wobbling rotation takes place. To the best of our knowledge, no direct measurements of HB lifetimes have been made.
In summary, in this study, we have studied the HBs in three halide perovskites by analyzing MD simulations. These HBs have been defined from geometrical criteria as a function of the distance and angle between the atoms. The presented geometrical criteria are based on the angle−distance CDF and the IUPAC recommendations. We have characterized the statistical properties of the HBs by means of PDF, CDF, and time correlation functions. From the latter, we obtained HB lifetimes, which are smaller for C−H−Y bonds than for N− H−Y bonds (Y = Br and I). The HB lifetimes are similar in the solid solution (FAPbI 3 ) 7/8 (MAPbBr 3 ) 1/8 and in the end compounds FAPbI 3 and MAPbBr 3 .
On the other hand, analyzing the power spectrum, it has been possible to verify that the power spectra peaks related to N−H stretching modes, especially in MA, are modified as expected from the influence of HBs. However, this is not the case for the C−H stretching modes ( Figure 8). This analysis suggests that C−H−Y HBs are irrelevant for the vibrational spectra. The HBs in these compounds are weak and possibly undetectable at the high temperature that was set in the MD simulations. HBs should manifest at lower temperatures, particularly near the phase transitions. Our simulations were performed at 350 K, but this work establishes a theoretical and computational framework for the investigation of HBs that can be applied for the analysis of MD simulations at other temperatures.