Correlation between C-H•••Br and N-H•••Br hydrogen bond formation in perovskite CH3NH3PbBr3: a study based on statistical analysis

Academic Editor: Ilya G. Shenderovich Received: 29 April 2025 Revised: 30 May 2025 Accepted: 2 June 2025 Published: 4 June 2025 Citation: Garrote-Márquez, A.; Cruz Hernández, N.; Menéndez-Proupin, E. Correlation Between C–H···Br and N–H···Br Hydrogen Bond Formation in Perovskite CH3NH3PbBr3: A Study Based on Statistical Analysis. Solids 2025,6, 29. https://doi.org/10.3390/ solids6020029 Copyright: © 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/ licenses/by/4.0/). Article Correlation Between C–H···Br and N–H···Br Hydrogen Bond Formation in Perovskite CH3NH3PbBr3: A Study Based on Statistical Analysis Alejandro Garrote-Márquez , Norge Cruz Hernández and Eduardo Menéndez-Proupin * Departamento de Física Aplicada I, Escuela Politécnica Superior, Universidad de Sevilla, E-41011 Seville, Spain; [email protected] (A.G.-M.); [email protected] (N.C.H.) *Correspondence: [email protected] Abstract: This study investigates the potential correlation between C—H ··· Br and N— H ··· Br hydrogen bonds in CH 3 NH 3 PbBr 3 over a broad temperature range (50–350 K), using a statistical analysis of molecular dynamics simulations. The analysis focused on quantifying the relationship between both hydrogen bond types via Pearson and Spearman correlation coefficients, derived from extensive datasets obtained from simulation trajectories. The results revealed a notable discrepancy between the two coefficients at low temperatures (T ≤ 125 K): While Spearman’s values suggested a strong monotonic correlation, Pearson’s values indicated a lack of linear association. Further analysis through data segmentation and block averaging demonstrated that the high Spearman coefficients at low temperatures were not statistically robust. At higher temperatures (T > 125 K), both correlation coefficients consistently exhibited low values, confirming the absence of meaningful correlation. These findings suggest that the formation of C–H ··· Br and N–H ··· Br hydrogen bonds occurs independently, with no evidence of cooperative behavior. Keywords: hydrogen bonding; halide perovskites; statistical analysis; hydrogen bond correlation coefficients 1. Introduction The study of tin-lead halide perovskites has garnered increasing attention in the scientific community, particularly in the fields of optoelectronic devices, solar cells, and LEDs. These materials have emerged as a highly promising alternative for next-generation photovoltaic applications due to their exceptional optoelectronic properties, ease of fabrication, and remarkable power conversion efficiency [ 1 – 6 ]. Since the first reports on perovskite solar cells (PSCs) [ 7 ], their efficiency has increased dramatically [ 8 – 12 ], reaching over 26%, making them comparable to conventional silicon solar cells [ 9 , 13 , 14 ]. However, large-scale commercial implementation is not yet possible due to the lack of long-term stability of perovskite materials. One of the primary challenges associated with perovskite solar cells is their sensitivity to environmental factors such as humidity, oxygen exposure, temperature fluctuations, and prolonged light exposure [ 15 , 16 ]. These external conditions can accelerate the degradation of perovskite films, significantly reducing their operational lifespan, unlike silicon-based solar cells, which exhibit durability exceeding 25 years. To mitigate these issues, extensive research has focused on various strategies, including compositional engineering, improvements in encapsulation techniques, and structural modifications aimed at enhancing the robustness of perovskite materials under adverse conditions [17–27]. Solids 2025,6, 29 https://doi.org/10.3390/solids6020029 Solids 2025,6, 29 2 of 15 Figure 1shows a unit cell of the metal–organic halide perovskite CH 3 NH 3 PbBr 3 . It can be described as a methylammonium cation, MA=CH 3 NH 3 , inside a cubo-octahedral cavity composed of 8 Pb cations and 12 Br anions. This unit cell is repeated quasi-periodically as illustrated in Figure 2. The deviations from strict periodicity, apparent in Figure 2, are due to thermal motion, as these figures are snapshots from molecular dynamics (MD) simulations. Some effects of thermal motion are particularly strong in halide perovskites, such as disorder in cation orientation, large displacements of the halide anions, and a random network of hydrogen bonds (HBs). The organic cation MA is agglutinated by covalent C–N, N–H and C–H bonds, which exhibit lengths of approximately 1.5 Å, 1.0 Å, and 1.0 Å, respectively. These bond lengths fluctuate slightly due to thermal motion, but they remain confined within a narrow distribution. A covalent bond involves the pairing of two electrons with opposite spin projections, resulting in increased electron density between two atomic nuclei. This also allows us to understand, on qualitative grounds, the positive charge of MA. The molecular energy is minimized by filling the valence shells of C and N atoms with eight electrons and the valence shell of H with two electrons. This count includes the electrons shared with the covalently bonded atom. The 5-electron valence shell of N is completed with the electrons shared by the C atom and two H atoms, leaving one lone pair. The neutral molecule is CH 3 NH 2 , but in the ionic PbBr 3 framework, charge neutrality and energy minimization are achieved by capturing a proton attracted by the lone pair of electrons on the N atom. This qualitative picture is strongly backed by numerical quantum mechanics calculations. Figure 1. A unit cell of CH 3 NH 3 PbBr 3 , showing two X–H ··· Y hydrogen bonds (N–H ··· Br and C– H ··· Br) as dotted lines, as well as the X–Y distances (dashed lines and values) and the X–H–Y angles that define them. A different matter is the interaction of the organic cation with the inorganic lead-halide framework. To a great extent, this interaction is ionic, given that the inorganic sublattice is negatively charged. However, the cation can rotate inside the cubo-octahedral cavity, except at low temperatures. In this context, non-covalent and non-ionic interactions modify the cation movement, and they also contribute to the cohesive energy of the perovskite crystal. One of the most important factors of these interactions is hydrogen bonding [ 28 , 29 ]. These bonds play a fundamental role in maintaining the structural integrity of hybrid perovskites, which consist of both organic and inorganic components. The strength and configuration of HBs directly influence the overall stability and resilience of perovskite materials against degradation factors. Hydrogen bonding is a weak interaction between a highly electronegative atom Y, called the acceptor, and a hydrogen atom that is covalently bonded to another atom X (this pair is called the donor). The HB is denoted X–H ··· Y, where the three-dot line Solids 2025,6, 29 3 of 15 indicates a weak non-covalent interaction. According to the definition of the International Union of Pure and Applied Chemistry (IUPAC) [ 30 , 31 ], both X and Y atoms must be more electronegative than H, causing some charge transfer in the heteropolar covalent bond X–H, such that H becomes slightly positively charged and undergoes electrostatic attraction with the donor Y. In a metal–organic halide perovskite like CH3NH3PbBr3, the acceptor Y = Br, while the donor can be either N–H or C–H in the organic cation. The IUPAC definition also establishes that the HB is directional, such that the angle X–H–Y must approach 180 ◦ , and the electron density must display a (3, − 1) critical point in the H ··· Y path. All these criteria are fulfilled in metal–organic halide perovskites [ 12 , 29 , 32 , 33 ]. The HBs are represented in Figures 1and 2by means of thick dotted lines. The N–H ··· Br HBs are believed to be significantly stronger than the C–H ··· Br ones as well as those involving other halides. First, the electronegativity of carbon is relatively low compared to that of nitrogen, so the C–H group has not traditionally been considered a valid donor, although this view has changed with the current IUPAC definition of HB [ 30 , 31 ]. Second, the H–Br distance is about 1 Å shorter in the case of the N–H donor than for C–H. Third, the power spectrum of the N–H vibrational stretching modes is dramatically modified by the insertion of the organic cation into the perovskite lattice, with a much smaller effect for the C–H stretching modes [29,33]. The well-known rotational motion of the organic cation in metal–organic halide perovskites implies that HBs form and break dynamically. This process can be studied by means of MD simulations and it can be quantified by means of a dimer autocorrelation function, from which a lifetime (LT) can be derived [ 34 ]. It is natural to expect that the LTs of N–H ··· Br bonds are longer than those of C–H ··· Br bonds. However, MD simulations of CH 3 NH 3 PbBr 3 [ 33 ] have shown rather similar LTs across a wide range of temperatures for both kinds of HB. Notably, the LTs are finite even for the low-temperature orthorhombic phase, where the cations do not rotate (at least the C–N axis). Moreover, the LTs’ temperature dependence is consistent with the Arrhenius law, and the activation energies are also similar for both kinds of HBs. This similarity of LTs and activation energies seems difficult to reconcile with the idea of one kind of HB being stronger than the other. The existence of C–H ··· Br bonds is supported by the reduced electron-density gradient calculations, but they do not manifest in the vibration properties or any other observable factors. The CH 3 NH 3 cation is confined in a cubo-octahedral cavity limited by Pb and Br ions. When the ammonium group NH 3 binds to one, two, or three Br ions, the opposite methyl group, CH 3 , has available Br ions to bind at the opposite side of the cavity. It might be possible that the driving force for the formation of C–H ··· Br bonds were the N–H ··· Br bonds at the opposite side of the cation. Henceforth, we present a study of the correlation between the two kinds of HBs of a single cation along the MD trajectory. Figure 2. Representation of the MAPbBr 3 perovskite structure. Hydrogen bonds are indicated by thin dashed lines in red and bluish colour for C–H ··· Br and N–H ··· Br bonds, respectively. Image created with VMD [35]. Solids 2025,6, 29 4 of 15 A deeper understanding of this possible correlation is essential for improving perovskite stability. If a feedback mechanism exists, it could be leveraged to design more robust materials by enhancing cooperative hydrogen-bonding networks. Correlation analysis is a fundamental tool in MD studies, particularly when aiming to understand the relationships between different dynamic variables in complex systems. Previous studies have demonstrated the importance of correlation analysis in identifying hidden patterns and dependencies among variables in molecular systems [ 36 – 39 ]. In this context, Pearson and Spearman correlation coefficients are commonly used to assess the strength and direction of relationships between variables [40,41]. Pearson’s correlation coefficient is a measure of the linear relationship between two continuous variables. It evaluates how changes in one variable correspond proportionally to changes in another. This coefficient assumes that the relationship between variables is linear and that both datasets follow a normal distribution. The main advantage of Pearson’s correlation is its ability to quantify the precise degree of linear dependence between variables. However, its primary limitation is its sensitivity to outliers and its inability to accurately capture non-linear relationships. On the other hand, Spearman’s correlation coefficient is a rank-based measure that assesses the strength and direction of a monotonic relationship between two variables. Unlike Pearson’s correlation, Spearman’s coefficient does not require the data to be normally distributed, making it particularly useful for datasets with skewed distributions or nonlinear relationships. It is also less sensitive to outliers, providing a more robust measure of association in cases where extreme values may distort linear analyses. 2. Materials and Methods 2.1. Molecular Dynamics Simulations and Hydrogen Bond Characterization The source of the data is the same as in Ref. [ 33 ], where the MD simulations were performed in the canonical ensemble with the Nosé thermostat [ 42 ] using the VASP code [ 43 ]. The simulation cell includes 2592 atoms, which corresponds to a 6 × 6 × 6 supercell of the primitive cell of the cubic phase. The dimensions of the simulation cell for each temperature were determined as the average values from a previous simulation in the isothermal-isobaric ensemble. The equations of motion were integrated with a time step of 0.5 fs, while the coordinates were saved every two timesteps. The interatomic forces were computed using machine-learning force fields trained as described by Liang X et al. [44]. In previous studies [ 29 , 33 ], the identification of HBs was based on two geometric criteria, namely the distance between the hydrogen atom and the bromide ion, and the angle formed by the donor-hydrogen-acceptor triplet. Specifically, N–H ··· Br interactions are defined by an H ··· Br distance shorter than 3 Å and a N–H–Br angle greater than 135 ◦ , while for C–H ··· Br interactions, the distance threshold was extended to 4 Å, maintaining the same angular condition. These cutoffs reflect the observed maxima in the combined distribution function (CDF) analysis [ 34 ], and they are consistent with the characteristic ranges found in previous computational studies of similar systems [29,45]. In this study, the Visual Molecular Dynamics (VMD) software [ 35 ] (version 1.9.4a55 OpenGL) has been used for both visualizing the HBs and for quantifying their numbers. In VMD, a X–H ··· Y HB is defined by the X–Y distance instead of the H–Y distance and the above explained condition for the angle (X–H–Y). Because of this, the N–Br distance’s upper limit was set to 4 Å, that is, 1 Å larger than the H–Br distance limit. The same procedure was followed for the C–Br distance, setting its upper limit at 5 Å. The equivalence of the HB definition in terms of X–Y or H–Y distances was shown in Ref. [29]. Solids 2025,6, 29 5 of 15 2.2. Correlation Analysis The dataset used in this study was built by tracking the trajectory of an MA molecule, along with the atoms forming its cubo-octahedral cavity. For each timestep, the number of HBs was quantified, distinguishing those formed with the nitrogen atom from those formed with the carbon atom. A dataset consists of a sequence of pairs of integer numbers—the numbers of C–H···Br and N–H···Br HBs at each MD frame—for each temperature. To explore whether a meaningful relationship exists between the behavior of C–H ··· Br and N–H ··· Br HBs, we applied Pearson’s and Spearman’s correlation coefficients. These statistical tools enabled us to quantify both linear and monotonic associations in the time series extracted from MD simulations. Pearson’s correlation coefficient (r) measures the linear correlation between two continuous variables. It quantifies the degree to which the variables move together linearly. Pearson’s coefficient ranges between − 1 and +1, where +1 indicates a perfect positive linear relationship, − 1 a perfect negative linear relationship, and 0 indicates no linear relationship. The mathematical expression for Pearson’s correlation coefficient is given by r=∑(xi−x)(yi−y)/rh∑(xi−x)2∑(yi−y)2i, (1) where xi and yi represent individual observations from each dataset and x and y are the mean values of each dataset. Spearman’s correlation coefficient ( ρ ), on the other hand, assesses the strength and direction of a monotonic relationship between two variables. It is especially useful for non-normally distributed data or data that do not necessarily exhibit linear relationships. Spearman’s coefficient also ranges between − 1 and +1, where values close to ± 1 indicate a strong monotonic relationship, and values near 0 indicate no monotonic relationship. The formulation for Spearman’s rank correlation coefficient is as follows: ρ=1−h6∑di2/nn2−1i, (2) where di is the difference between the ranks of each paired observation and nis the number of paired observations. The rank, given by the EXCEL function RANK.AVG, refers to the relative position of a data point within an ordered dataset. Instead of using the raw numerical values of the variables being compared, Spearman’s method evaluates the strength of a monotonic relationship by assigning a rank to each value based on its position in the sorted list from the smallest to the largest. By converting the data into ranks, Spearman’s correlation becomes less sensitive to the actual magnitude of the values and instead focuses on the overall order or trend between the two variables. This characteristic makes it particularly effective for detecting consistent, non-linear relationships, even in the presence of outliers or skewed distributions. To perform these correlation analyses, both datasets were initially organized and carefully prepared to ensure consistency. A descriptive statistical analysis was initially performed to understand the basic statistical characteristics of the datasets, including the calculation of mean, median, standard deviation, variance, and range, to ensure data suitability for subsequent correlation analysis. The strength of the obtained correlations was interpreted according to experience. Typically, absolute correlation coefficient values between − 0.7 and 0.7 are regarded as weak correlations, and values above ±0.7 indicate strong correlations. To ensure that the analyzed data corresponded to thermal equilibrium conditions across the temperature range (from 50 K to 350 K), a convergence analysis of the datasets was performed. Specifically, for T≤125 K, the last 50,000, 40,000, 30,000, 20,000, and Solids 2025,6, 29 6 of 15 10,000 values from each dataset as well as the full 60,000 step dataset were analyzed. For T> 125 K, the MDs contain only 20,000 frames. Hence, 20,000 and the last 15,000, 10,000, and 5000 values from each dataset were examined. Table 1shows the results of the Pearson and Spearman coefficients for each temperature and subset indicated above. Table 1. Pearson and Spearman coefficients with data analyzed in thermal equilibrium. PEARSON Temperature/Frames 10,000 20,000 30,000 40,000 50,000 60,000 50 −0.24 0.13 0.10 0.08 0.12 0.08 70 −0.05 −0.02 −0.05 −0.04 −0.05 −0.02 90 −0.02 −0.02 0.01 −0.03 −0.02 −0.02 100 0.10 0.06 0.05 0.02 0.02 0.01 110 0.08 0.04 −0.03 0.01 0.00 −0.01 125 0.12 0.14 0.16 0.15 0.14 0.15 Temperature/Frames 5000 10,000 15,000 20,000 150 0.07 0.10 0.22 0.20 175 0.19 0.20 0.13 0.15 200 0.15 0.11 0.10 0.11 215 0.01 0.12 0.12 0.13 235 0.19 0.15 0.14 0.15 250 0.12 0.20 0.18 0.15 275 0.10 0.12 0.09 0.10 300 0.08 0.09 0.10 0.11 325 0.10 0.14 0.11 0.12 350 0.13 0.16 0.13 0.14 SPEARMAN Temperature/Frames 10,000 20,000 30,000 40,000 50,000 60,000 50 0.88 0.86 0.89 0.89 0.89 0.90 70 0.84 0.83 0.80 0.80 0.82 0.84 90 0.79 0.75 0.73 0.73 0.73 0.73 100 0.57 0.61 0.64 0.62 0.64 0.62 110 0.77 0.75 0.68 0.62 0.61 0.59 125 0.46 0.48 0.52 0.53 0.48 0.46 Temperature/Frames 5000 10,000 15,000 20,000 150 0.37 0.38 0.36 0.35 175 0.40 0.30 0.31 0.32 200 0.30 0.26 0.23 0.25 215 0.31 0.32 0.34 0.31 235 0.23 0.26 0.20 0.27 250 0.29 0.30 0.31 0.30 275 0.33 0.28 0.30 0.27 300 0.34 0.30 0.28 0.27 325 0.27 0.22 0.24 0.27 350 0.24 0.21 0.27 0.27 2.3. Data Averaging and Correlation Stability Evaluation For each dataset of T ≤ 125 K, a block-averaging procedure was implemented to smooth the data without losing statistically significant information. The block-averaging procedure was performed by grouping consecutive values within the original dataset and computing the mean for each group, generating a new set composed of those mean values. The block sizes of 30, 150, 1000, 1500, 3000, and 10,000 were used in such a Solids 2025,6, 29 7 of 15 way that the new sets contained 2000, 400, 60, 40, 20, and 6 elements, respectively. This approach facilitated the creation of more manageable data subsets, preserving overall trends while minimizing random fluctuations that could introduce bias. Pearson and Spearman correlation coefficients were then computed for these reduced datasets and compared to those obtained using the complete dataset. This procedure also allowed us to solve an apparent contradiction between the Pearson and Spearman correlation coefficients for T ≤ 125 K derived from the integer nature and the prevalence of value 3 (the number of HBs of NH3and CH3groups) in the original dataset. To assess the stability of correlation as a function of the number of averaged data points, correlation coefficients were visually analyzed using scatter plots. Different block-averaged sets were represented along with their respective coefficients, allowing for an evaluation of whether the correlation remained stable as data points were reduced or if there was a significant loss of information due to averaging. Additionally, linear regressions were performed on the subsets to analyze whether the relationship between variables remained consistent across different averaging levels. This approach enabled the determination of the optimal number of data points required to obtain stable and representative correlation coefficients, minimizing statistical noise and ensuring the validity of the conclusions drawn from the analysis. Table 2shows the results of the Pearson and Spearman coefficients for each temperature as a function of the number of elements of the block-averaged set, and it will be discussed in Section 3. Table 2. Pearson and Spearman coefficients with block-averaged data. Pearson 60,000 2000 400 200 60 40 20 6 50 K 0.08 0.15 0.36 0.46 0.54 0.65 0.68 0.57 70 K −0.02 −0.04 −0.04 −0.13 −0.20 −0.38 −0.51 −0.56 90 K −0.02 −0.03 −0.01 −0.10 −0.14 −0.39 −0.56 −0.61 100 K 0.01 0.02 0.05 0.08 0.28 0.29 0.18 0.24 110 K −0.01 −0.01 0.08 0.05 0.31 0.22 0.22 −0.11 125 K 0.15 0.21 0.40 0.44 0.46 0.53 0.48 0.65 Spearman 60,000 2000 400 200 60 40 20 6 50 K 0.90 0.76 0.49 0.33 0.29 0.16 0.13 0.54 70 K 0.84 0.57 0.10 −0.10 −0.21 −0.33 −0.46 −0.43 90 K 0.73 0.38 −0.06 −0.19 −0.20 −0.44 −0.42 −0.60 100 K 0.62 0.26 −0.06 −0.15 −0.09 0.04 0.02 0.03 110 K 0.59 0.24 −0.03 −0.12 0.19 0.06 0.11 −0.03 125 K 0.46 0.27 0.28 0.36 0.41 0.54 0.47 0.77 3. Results 3.1. Correlation Analysis with the Complete Dataset Table 3shows, for each temperature, the total number of MD frames, the HB lifetimes, and the Pearson and Spearman coefficients, considering the initial (before block averaging) dataset. As can be seen in Figure 3, the trends in the Pearson and Spearman coefficients were distinct across the entire temperature range examined (50–350 K). At lower temperatures ( ≤ 125 K), the Spearman coefficients presented high values, consistently above 0.5, indicating a possible relationship between the variables analyzed. However, the Pearson coefficient remained low in this range, fluctuating between −0.02 and 0.14. Solids 2025,6, 29 8 of 15 Table 3. Pearson and Spearman coefficients, number of frames of the MD simulations, and lifetimes. The τvalues have been obtained from Table 1 of Ref. [33]. T(K) Frames Pearson Spearman τ(ps) HBs_C τ(ps) HBs_N 50 60,000 0.08 0.90 - - 70 60,000 −0.02 0.84 7.6408 6.6992 90 60,000 −0.02 0.73 3.4374 2.7417 100 60,000 0.01 0.62 2.5533 2.0081 110 60,000 −0.01 0.59 1.5975 1.4053 125 60,000 0.15 0.46 0.7176 0.9616 150 20,000 0.20 0.35 0.5085 0.6800 175 20,000 0.15 0.32 0.3878 0.5150 200 20,000 0.11 0.25 0.3409 0.4166 215 20,000 0.13 0.31 0.3044 0.3830 235 20,000 0.15 0.27 0.2464 0.3383 250 20,000 0.15 0.30 0.2349 0.3188 275 20,000 0.10 0.27 0.2082 0.2787 300 20,000 0.11 0.27 0.1933 0.2532 325 20,000 0.12 0.27 0.1774 0.2267 350 20,000 0.14 0.27 0.1635 0.2035 Figure 3. Values of the Pearson and Spearman coefficients as a function of temperature. For T > 125 K, it was observed that Spearman coefficients decreased significantly, with values fluctuating around 0.3. This trend suggests a reduction in the strength of the monotonic correlation between variables at higher temperatures. The Pearson coefficient also showed an upward trend at intermediate temperatures but remained relatively low, indicating that the relationship between the variables remained weakly linear or was influenced by additional factors. However, it can be observed in Figure 1that the Pearson and Spearman coefficients at temperatures above 125 K are proportional and both indicate low correlation. In contrast, at lower temperatures, the Pearson and Spearman coefficients are qualitatively different, indicating low and high correlations, respectively. This could indicate two possible situations: either that, at low temperatures, the system is not thermalized, or that the behavior of the data ceases to be linear at low temperatures. This latter behavior might be interpreted as non-linear correlation between both types of HBs. At sufficiently low temperatures, the structure of the system becomes more ordered with less thermal contribution, which facilitates the emergence of HB networks in which the formation or breaking of one type of bond determines the presence of the other. Solids 2025,6, 29 9 of 15 In this context, a correlation between both types of HBs—as suggested by a high Spearman coefficient—appears plausible. If the formation of an N–H ··· Br bond forces the cation to adopt a specific orientation, it is reasonable to assume that such orientation also influences the likelihood of forming a C–H ··· Br bond in the opposite direction and vice versa. This possible cooperative relationship, although not necessarily linear, could arise from the spatial constraints imposed by the crystal structure and the need to optimize interactions between the organic cation and the inorganic framework. Within the cubooctahedral cavity of the perovskite lattice, the MA cation is confined by its surrounding Pb 8 Br 12 environment, which limits its freedom of motion and imposes specific steric and electrostatic conditions. When an N–H ··· Br bond forms, the resulting orientation of the NH 3+ group inevitably repositions the CH 3 end of the cation, potentially placing one or more hydrogen atoms of the methyl group in closer proximity to bromide ions. This geometrical rearrangement may facilitate the formation of a C–H ··· Br interaction simply due to spatial alignment without requiring a direct energetic coupling between the two types of HBs. Moreover, from an energetic standpoint, the system may favor configurations that maximize the number of stabilizing interactions within the cavity, even if the individual contributions of C–H ··· Br bonds are weaker. The presence of multiple simultaneous weak interactions can contribute collectively to the stabilization of the local structure, which may manifest statistically as a monotonic correlation between bond types. Importantly, this scenario does not imply that the formation of one bond directly triggers the other but rather that both emerge as compatible outcomes within a constrained conformational landscape. This subtle form of cooperativity—driven more by spatial accommodation than by direct coupling—could explain the emergence of monotonic patterns in the HB dynamics at low temperatures. However, due to the marked difference between the Pearson and Spearman coefficients—and, more importantly, the lack of consistency in the high Spearman values (only exceeding 0.8 at 70 K and 50 K)—it becomes necessary to carry out more detailed analyses, as developed in the following sections. Sections 3.2 and 3.3 are devoted to analyzing the difference in both types of correlation coefficients. 3.2. Correlation Analysis with Fragmented Dataset To verify that our correlation results accurately reflect equilibrium behavior and are not biased by initial conditions, we recalculated the correlation coefficients using only (later) fragments of each trajectory. Table 1shows that, at all temperatures, the correlation coefficients exhibit only fluctuations when the dataset is fragmented. Focusing on the low-temperature examples, at 50 K, the Pearson correlation using the full 60,000-frame dataset is r≈ 0.08, whereas the Spearman correlation is ρ≈ 0.90. If, instead, we use only the last 10,000 frames of that trajectory, we obtain r≈ − 0.24 and ρ≈ 0.88. Using the last 20,000 frames, we get r ≈ 0.13 and ρ≈ 0.86, and using the last 50,000 frames, we get r≈0.12 and ρ≈ 0.89. These values are very close to each other, given the statistical uncertainty (on the order of a few hundredths). The small variations in Pearson’s rat 50 K (some slightly negative, some slightly positive) likely reflect random fluctuations due to the finite sample size rather than a systematic trend. Critically, the Spearman coefficient remains consistently high (0.86–0.90) regardless of segment length and the Pearson coefficient stays close to zero. Similarly, at 70 K, whether we use 10,000 frames or 60,000 frames, Spearman’s ρ remains around 0.8 (0.80–0.84) and Pearson’s r stays between –0.05 and –0.02. At 125 K, using the last 10,000 frames gives ρ≈ 0.46 and r ≈ 0.12, while using the