Exploring Librational Pathways with on-the-Fly Machine-Learning Force Fields: Methylammonium Molecules in MAPbX3 (X = I, Br, Cl) Perovskites

Two seemingly similar crystal structures of the low-temperature (∼100 K) MAPbX3 (X = I, Br, Cl) perovskites, but with different relative methylammonium (MA) ordering, have appeared as representatives of this orthorhombic phase. Distinguishing them by X-ray diffraction experiments is difficult, and conventional first-principles-based molecular dynamics approaches are often too computationally intensive to be feasible. Therefore, to determine the thermodynamically stable structure, we use a recently introduced on-the-fly machine-learning force field method, which reduces the computation time from years to days. The molecules exhibit a large degree of anharmonic motion depending on temperature: that is, rattling, twisting, and tumbling. We observe the crystal’s “librational pathways” while slowly heating it in isothermal–isobaric simulations. Marked differences in the thermal evolution of structural parameters allow us to determine the real structure of the system via a comparison with experimentally determined crystal structures.


■ INTRODUCTION
The crystal structure of hybrid halide perovskites is a topic of study that has surfaced several times in the past four decades. X-ray powder diffraction experiments of Weber et al. on methylammonium (MA)-PbX 3 with halogens X = {I, Br, Cl} have established a high-temperature cubic phase for all X. 1 A perovskite structure is formed by PbX 6 corner-sharing octahedra enclosing the MA molecules. In later years, the low-temperature phases and librational modes of the MA molecule at various temperatures were studied. 2−6 The development of high-efficiency solar cells based on these perovskites sparked a rival of interest in its structure characterization. 7−10 It was shown by high-quality powder neutron-diffraction experiments 9 that the low-temperature orthorhombic phase of MAPbI 3 actually belongs to the Pnma space group. This is the space group that was also determined for MAPbBr 3 11 and MAPbCl 3 . 12 The potential for optoelectronic applications raised new questions that are all (in)directly related to the atomic structure: the effect of MA rotation on charge dynamics, 13,14 dynamic or permanent deformations of the PbX 6 octahedra, 15−18 the extent of electron−phonon coupling, 19−21 and the Rashba effect, 22−27 to name a few. Considerable progress has been made, but consensus has not always been achieved. This is in part the result of differences in interpretation of the local microscopic structure. The disorder, be it static or dynamic, of the molecular C−N axes is a wellknown problem for diffraction techniques that makes it difficult to determine their precise orientation. 18,28 Firstprinciples (FP) methods such as density functional theory (DFT) have shown to be very useful for the determination of the crystal structure by augmenting the experimentally resolved inorganic framework with the ordering of the molecules. 29−37 However, even though commonly used density functional approximations (DFAs) have the required chemical accuracy, 38 their computational complexity prohibits the large length-andtimescale molecular dynamics (MD) calculations necessary to resolve the free energy landscape and thereby the finitetemperature crystal structure. 31 We will use the on-the-fly machine-learning force field (MLFF) method, 39,40 which makes it possible to explore the full diversity of atomic structures while going through the entropy-driven phase transformations in hybrid perovskites. This method substantially reduces the computational cost while retaining near-FP accuracy. Recently, it has been shown to be capable to resolve the orthorhombic−tetragonal (Ort−Tet) and tetragonal− cubic (Tet−Cub) phase transitions in MAPbI 3 and the inorganic halide perovskites CsPbX 3 in good agreement with experiment. 39 Furthermore, it can be systematically extended to describe mixed MA x FA 1−x PbI 3 perovskites under isothermal−isobaric conditions. 41 The starting points of our search for the low-temperature (∼100 K) orthorhombic structure of MAPbX 3 are two seemingly similar but distinctly different structures: sA and sB. They have the same lattice vectors and inorganic coordinates but a different molecular ordering pattern, as sketched in Figure 1. We have labeled the lattice vectors such that the molecules lie in the ab-plane. sB is created out of sA by an in-plane rotation of half of the molecules by 180°, as indicated by the curved arrows. Note that in both arrangements, the neighboring molecules in the c-direction (not shown in the figure) are antiparallel. These structures have been prepared in a 2 × 2 × 2 supercell such that it accommodates the ( a a a 2 , 2 , 2 p p p ) basis of X = I, Br as well as the (2a p , 2a p , 2a p ) basis of X = Cl, where a p is the pseudocubic lattice constant of the parent Pm3̅ m cell. 12 The molecules in the often-referenced experimental (Pnma) structures for X = I, Br, and Cl of refs 8 11, and 12 are arranged as in the sB, sA, and sA configuration, respectively. Other experimental works did not distinguish between these two arrangements 7,42 since refinement of the model structure by permuting N and C with respect to the measured diffraction spectra does not lead to significant improvements of the fit. 18 To date, even for the extensively studied MAPbI 3 perovskite, different studies report opposite arrangements: sA 32,43,44 and sB. 29,34,45,46 This is unexpected because if we focus only on the dipole moment of the MA molecule and compute the total electrostatic energy in the point-dipole approximation, then the sB pattern is clearly favored. This pattern shows a closer resemblance to the "head−tails" ground state of a point-dipole model. 47 In this work, we will analyze the "librational pathways" of the MA molecules and PbX 6 octahedra in MAPbX 3 and use them to identify the most representative low-temperature orthorhombic structure. We sample structures of the crystal by slowly heating up two plausible low-temperature structures (sA and sB) in isothermal−isobaric MD simulations. Structures on the explored pathways through the structural phase space are thermodynamically linked to the starting configuration and result in marked differences in lattice and order parameters that are compared to temperature-dependent diffraction studies.

■ COMPUTATIONAL DETAILS
The DFT calculations are performed with the projectoraugmented wave method, 48 as implemented in the VASP code 49,50 using the metagradient-corrected SCAN 51 DFA, which has shown good performance when compared to highquality many-body perturbation theory reference calculations. 38 A plane-wave basis with a cutoff of 350 eV, Gaussian smearing with a width of 10 meV, and 4 (I, Br) or 8 (Cl) kpoints of the Γ-centered 2 × 2 × 2 Monkhorst−Pack grid are set, which suffice to obtain the required accuracy of the calculations. 36 The computed lattice parameters as a function of temperature should (qualitatively) agree with experiment over the whole temperature range. Therefore, by not limiting the study to a 0 K DFT-based relaxation of the internal energy, biases related to the chosen DFA can be detected. 38 Before starting the MD simulations, the starting structures of Figure 1 were shortly relaxed by a conjugate gradient algorithm.
MLFFs are trained during MD simulations with the VASP based on calculated total energies, forces, and stress tensors for automatically (on-the-fly) selected structures in the isothermal−isobaric ensemble. This approach is described in detail in refs 39 and 40. In short, a Bayesian error estimation of the predicted forces is used to select either DFT or MLFF forces to propagate the structure in time (t n → t n+1 ). Whenever the predicted errors exceed the threshold, a new reference structure is picked up, a DFT calculation is performed, and the coefficients of the MLFF are reoptimized. In Figure 3c,f, a "density-of-states"-like function of the temperature (note, equivalent to simulation time) shows when in the training MD most DFT calculations were performed. It is calculated where δ(T) is a Lorentzian function. This function is normalized to the total number of DFT reference structures picked up in training, N ref = ∫ ρ FP (T) dT. The automatically picked-up reference structures form a minimal training database (containing total energies, forces, stress tensors, and atomic coordinates) that is well spread over the available structural phase space. We have shared this database via the 4TU.DataBase repository 52 to encourage development of ML potentials based on minimalistic data sets.
A variant of the GAP-SOAP 53,54 method is used as a descriptor of the local atomic configuration around each atom. Within a cutoff of 7 Å, a two-body radial probability distribution ρ i (2) (r) is built, as well as three-body angular distribution ρ i (3) (r,s,θ) within a cutoff of 4 Å. The atomic coordinates are smeared in the distributions by placing Gaussians with a width of 0.5 Å. The obtained distributions are projected on a finite basis set of spherical Bessel functions multiplied with spherical harmonics. The Bessel functions are of the order 6 and 7 for the radial and angular part, respectively. Only the angular part has a maximal angular momentum of l max = 6. The coefficients of the projections are gathered in the descriptor vector X i . A kernel-based regression method 55 is applied to map the descriptor to a local atomic energy. The similarity between two local configurations is c a l c u l a t e d b y a p o l y n o m i a l k e r n e l f u n c t i o n : = · + · K X X X X X X ( , ) 1/2( ) 1/2( ) . For each of the six on-the-fly heating trajectories, a Supplementary Movie was generated. A running average (window size of 25 K) over each atomic coordinate is computed to smoothen out high frequency movements. At

■ RESULTS AND DISCUSSION
To introduce the librational pathways of MAPbX 3 , we will illustrate them by weighted sums of pair distribution functions (PDFs) in Figure 2. The PDF for the atom types α and β is defined as where r i and r j are the coordinates of the N α and N β atoms, r ij = r i − r j , V is the volume of the simulation box, and ⟨ ⟩ denotes the ensemble average. In g inorganic (r), the pairs of framework components Pb−X, X−X, Pb−Pb are included, and only the C−N pairs of the MA molecules are included in g C,N (r). For all halides, X, we see that g inorganic (r) retains most of its structure throughout the whole temperature range and that g C,N (r) shows a transition whereby part of the order is lost. The intramolecular part (r ≈ 1.5 Å) remains intact, but the intermolecular pairs show dual peaks merging into a single broad peak centered around the nearest-neighbor distances of the consecutive cubic unit cells, that is, a a a 1 , 2 , 3 p p p . This is the result of the unfreezing of the molecules, whereby they reorient (C 4 ) rapidly. 4 Around room temperature, neighboring MA molecules are still dynamically correlated to their neighbors. 34 The differences in the DSF between the halides X are, among other things, related to different Pb−X bond lengths and to the relative orientation of the molecules in the low-temperature phase. As we will show hereafter, the thermodynamically stable molecular configurations at low temperature (sA or sB) depend on the halide type X.
The starting structures ( Figure 1) are slowly heated at a rate of 2 3 K ps using DFT-based MD with a Langevin thermostat and a time step of 2 fs. The PDFs at different temperatures have been obtained by partitioning the resulting trajectory in parts of equal length. In the NPT ensemble, all lattice degrees of freedom are allowed to change, as shown for MAPbBr 3 in Figure 3. The two plots correspond to the heating trajectories starting from the sA and sB structures. Thermal fluctuations in structural parameters are smoothed by applying running averages. To accelerate the MD, a MLFF is trained on-thefly, as described in refs 39 and 40. The algorithm switches between MLFF and DFT forces based on the predicted error of the MLFF. Structural reference configurations to train the MLFF are automatically picked up and by construction lie outside the already "learned" part of the phase space. This can be seen by the sharp increase in the density of FP calculations (ρ FP (T)) shown in Figure 3c This transition temperature is expected to be retarded because the system is out of thermal equilibrium as a result of the still considerable heating rate. Even though, the agreement with the experimental lattice parameters shown by the symbols in Figure 3 is remarkable. We show that the structural transformation and the related librational pathways of the molecules and octahedra (see Supplementary Movies) are accurately described in the on-thefly heating MD. This opens up the possibility of exploring many different perovskites because only ∼1.000 out of the total of 150.000 MD steps per heating run were DFT calculations. This reduces the computing time from years to days.
Librational Pathway Analysis. Figure 3 indicates that apart from a small step in the c lattice constant in the sB case, both trajectories qualitatively and quantitatively agree with the experimental data. This means that the lattice parameters alone provide insufficient information to select either the sA or the sB as the thermodynamically stable low-temperature structure. Therefore, we analyze the structural motif presented by the atomic coordinates in the six heating trajectories. First, we solely extract the orientation of the molecular C−N axes in time. The total electrostatic energy (H lr ) corresponding to the dipole moments of the molecules is calculated. 47 In short, this point-dipole model assumes a fixed dipole moment on all molecules, no screening, and includes all dipole−dipole interactions up to the third nearest neighbor. In Figure 4a, H lr is plotted as a function of temperature for the three halides starting from the sA (black lines) and sB (red lines) configurations. The sB configuration is clearly lower in energy. With increasing temperature, the molecules flip/reorder and the stable arrangement is broken down, eventually leading to a disordered state with H lr = 0. 56 For X = I, the sA pattern flips to the stable sB pattern before reaching the Ort−Tet phase transition, which shows that the initial structure was out of equilibrium. However, for X = Br and Cl, the pattern remains largely frozen in until the phase transition. The sA configuration can only be stable at low-T if either a potential energy contribution arising from the inorganic framework or the entropy compensates this internal energy difference. The DFT/MLFF-calculated internal energy (E) is shown in Figure 4b as the energy difference (ΔE = E sA − E sB ) between the heating trajectory starting with the sA and the sB structure. At low temperature and ambient pressure, the volume and entropy contributions to the Gibbs free energy are small, and a sizable positive/negative ΔE would indicate that the sB/sA configuration is favored, respectively. For MAPbCl 3 , it is clear that the sA structure is favored even though its electrostatic energy in the dipole model was unfavored. Above ∼175 K, the difference between the two initial configurations has been lifted by thermally induced structural rearrangements. For MAPbI 3 and MAPbBr 3 , the situation is less clear; at low temperature, ΔE is positive; however, it is of the size of the fluctuations. Even for the fully DFT-relaxed (0 K) structures, ΔE values are small: 33, 3, and −72 meV/formula unit for I, Br, and Cl, respectively. Increasing the precision of the calculation by doubling the k-point grid density and applying the tetrahedron method 57 results in 33, −6, and −75 meV/f.u., respectively. This indicates that, especially for MAPbBr 3 , we cannot distinguish sA from sB based on the internal energy alone.
The changes of the structural motif as a function of the temperature are compared in Figure 5  The breakdown of the initial molecular order occurs at a lower temperature in the Br case and coincides with the Ort−Tet phase transformation. This shows that the sA and sB molecular orders are more energetically competitive in orthorhombic X = Br than in I and that an additional measure is required to determine the stable arrangement. The O values for the experimental lowtemperature structures indicated by the circles in Figure 5a provide this measure. Based on their comparison to the simulation, the Br-sA is the stable low-temperature orthorhombic structure.
To automatically classify the instantaneous crystallographic phase of all structures within the MD trajectory, we applied a new approach based on the O order parameter. This allows, for example, to still assign the orthorhombic phase to a structure that is strained in a box with a ≈ b ≈ c. Whenever the variance of the components of O is below a threshold, it is classified as Cub, and otherwise, it is Ort or Tet. We can then differentiate between the last two by counting the number (1 → Ort, 2 → Tet) of components larger than their mean value. The red line in Figure 5a shows the result of this phase classification. Note that no Tet phase in MAPbCl 3 is recognized. This could be caused by a very small temperature window in which the Tet phase is stable or that the c/a ratio is too small to be noticed in the supercell.
Lattice parameters as a function of temperature are shown in Figure 5b. The parameters are refined in the unit cell corresponding to the classified phase and converted to pseudocubic lattice parameters (a p ). Experimentally obtained parameters are shown by the circles. Surprisingly, the plots for MAPbBr 3 and MAPbCl 3 show good to very good agreement with the experiment. Especially for MAPbI 3 , we notice the   Combining all the above presented analysis leads to the favored initial configurations, I: sB, Br: sA, and Cl: sA. These configurations are highlighted by the colored rectangles in Figure 5.
Training System Size Dependence. The PDFs in Figure  2 are plotted beyond half of the simulation box width since the 2 × 2 × 2 supercell has an average width of 2a p . They are computed in 4 × 4 × 4 supercells, which are created by replicating the original 2 × 2 × 2 cell, enabling us to sample the PDF on the [0, 2a p ] domain. These results are compared to a training run performed with a 4 × 4 × 4 supercell. Specifically, we have made a test for MAPbBr 3 starting in the sA configuration and applied the same on-the-fly training directly on the 4 × 4 × 4 supercell. In this supercell, the k-points of the k-grid (applied with the 2 × 2 × 2 supercell) all fold down on the gamma point. For computational tractability, we slightly lowered the plane-wave cutoff from 350 to 300 eV and retrained the 2 × 2 × 2 cell with the same cutoff for a fair comparison. The smaller plane-wave basis results in higher Pulay stress, which slightly affects the volume but does not qualitatively change the crystal structure. 36 Figure 6 shows g inorganic (r) and g C,N (r) for the standard supercell and the 8 times larger cell. The g inorganic (r) of the two systems are almost identical. The main deviations in g C,N (r) are the result of a different temperature at which the system switches from the Ort to Tet phase. This is to be expected for simulations with finite system size. The MD trajectory is chaotic, and the transition does not occur at the same temperature even when the initial conditions are the same. This very good agreement indicates that the applied 2 × 2 × 2 supercell is large enough to capture the crystal symmetry. This is in agreement with a fully ab initio MD study of the system size dependence of MAPbI 3 going up to the 6 × 6 × 6 supercell. 34 The accuracy of the MLFF model over the entire collected data set of structures, which includes three different perovskite phases is very high. The DFT reference energy (U DFT ) and predicted MLFF energy (U MLFF ) for the test systems are plotted in Figure 7 and show a clear linear relation over a large energy range. Note that both point clouds nicely overlap, whereby the larger variations are, as would be expected, observed for the smaller supercell. The overall root-meansquare (rms) error on the energy is only 1.7 and 0.88 meV/ atom for the MLFF trained on the 2 × 2 × 2 and 4 × 4 × 4 supercells, respectively. This is small and of the same order of magnitude of state-of-the-art ML potentials (kernel-regression, 58    The Journal of Physical Chemistry C pubs.acs.org/JPCC Article octahedra are distorted. 12 This polar distortion is highlighted in Figure 1 by the green and orange lines indicating the difference between two Pb−Cl bond lengths (3.02 and 2.73 Å) in the crystallographic a-direction and is in good agreement with our simulations. However, we find no noticeable distortion of the PbI 6 and PbBr 6 octahedra above 80 K, as shown in Figure 8, again, in agreement with experiments of refs 10 11, however, opposite to the findings of refs 16 and 17 where a permanent octahedra distortion is found in MAPbBr 3 around room temperature. In Figure 8a, the distribution of the Pb−X bond lengths as a function of temperature is shown. The heating trajectories were cut in parts of equal length, and all bond lengths in the a-direction were added to the distribution. The low-temperature distribution for Pb−Cl has two peaks. This distortion is neither observed when starting from the Cl-sB structure, nor when starting from the Br-sB or I-sA structures. The distribution is well described by a combination of two Gaussian distribution functions, μ σ μ σ + ( , ) . For Cl-sA, the mean values (μ) and the standard deviations (σ) as a function of temperature are shown in Figure 8b. These values have been obtained from a separate heating trajectory with the finished MLFF on a 4 × 4 × 4 MAPbCl 3 supercell. The so-obtained distributions agree with Figure 8a and improve statistical accuracy. Experimentally determined bond lengths from refs 12 and 18 have been added to Figure 8b and agree within the standard deviation. In the simulations, a single Gaussian suffices above ∼175 K, that is, |μ 1 − μ 2 | < (σ 1 + σ 2 ). At this temperature, the octahedral polar distortion is no longer observed on time average and the crystal is in the cubic phase. As for I and Br, instantaneous distortions of the octahedra do occur at these elevated temperatures.
The following scenario now becomes plausible, the sA ordered molecules are stabilized in the MAPbCl 3 orthorhombic phase by an antiferroelectric striped pattern of dipolar octahedra in the a-direction. As argued in ref 12, the volume of the perovskite has to be sufficiently small to induce these distortions, whereby the "hard" MA deforms the "soft" octahedra. However, the striped pattern cannot be the only stabilization mechanism because no distortion is observed in Br-sA.
We would like to note that our "ensemble average" view on the structural model results in PDFs for MAPbBr 3 that qualitatively agree with those obtained from X-ray diffraction experiments of ref 17. In Figure 9a, g inorganic (r) has been plotted on the same length scale as Figure     The Journal of Physical Chemistry C pubs.acs.org/JPCC Article change in the 150−280 K temperature range apart from thermal broadening of the peaks. However, our approach classifies tetragonal and cubic structures within this range and does not indicate that an orthorhombic structure would be a better fit throughout this temperature range. Starting from sB results in structural changes (indicated by *), in disagreement with the PDF of ref 17. This is another indication that the sB pattern is not the stable low-temperature structure. The different structural interpretation of the crystal at elevated temperatures becomes apparent in the H−X bonding, as shown by g H,X (r) in Figure 9b. We find a halogen-dependent first peak position. This finding is different from the 2.5 Å peak observed in the neutron PDFs of both MAPbBr 3 and MAPbCl 3 . 18 MA C 3 Dynamics. We would like to note that training of a very accurate MLFF for the hybrid perovskites is not fully completed by the here-performed single heating run. Precise values for phase-transition temperatures, c/a ratios, and so forth were not the aim of this work but can be obtained with an accurate MLFF which enables long MD trajectories on large supercells. 39 Limits to the accuracy can be seen in the C 3 torsion/rotation degree of freedom of the molecule, for example. Figure 10 shows the NH 3 versus the CH 3 group dihedral angle (ϕ MA ) of all MAs in MAPbBr 3 in NPT ensembles with the finished MLFF. Torsion (ϕ MA < 60°) unfreezes at ∼25 K, 60 and also, rotations (ϕ MA > 60°) occur around our starting temperature (80 K) of the on-the-fly training. Figure 10b shows that ϕ MA,i (t) of a single molecule shows step-like behavior, occasionally jumping between planes separated by 120°and is superimposed by a fast oscillation. For each of the eight molecules in the supercell, a probability distribution of ϕ MA as in Figure 10c was made. We then calculate the C 3 rotational energy barrier in the two MLFFs of MAPbBr 3 (sA and sB) by a Boltzmann inversion of the distribution. A barrier was only assigned at a temperature T when the number of 120°rotations in the ∼100 ps long MD trajectories exceeded the number of molecules in the supercell. The error bars in Figure 10d correspond to ±σ, the standard deviation of the eight obtained barriers. The barrier is, within our statistical accuracy, temperature-independent and, surprisingly, different between sA and sB. It is tempting to conclude that the sA structure for Br affords a more facile C 3 rotational degree of freedom compared to the sB structure within the orthorhombic phase. However, the difference should not persist in the high-temperature cubic phase, in which the MAs are orientationally disordered. This should be a warning that training is not yet complete and the MLFF still shows a bias depending on the initial conditions.
We are able to measure the barriers of a single molecule in vacuum in the same manner. The absence of a surrounding Pb−Br framework does not destabilize the molecule. Barriers obtained in this way are slightly lower than the DFT value of 105 meV for MA in vacuum. This value was calculated as the difference in internal energies of the optimized ϕ MA = 0 and 60°molecule, for which all internal degrees of freedom were relaxed under the constraint of ϕ MA . Using these two structures, the barriers for the Br-sA and Br-sB MLFFs are 78 and 89 meV, respectively.
Since both Br-sA and Br-sB show C 3 rotations at temperatures below the Ort−Tet phase-transition temperature, we can conclude that the unfreezing of this motion does not drive it. This is in agreement with the findings of the computational study of Kieslich et al. 61 This phase transition is driven by an increase in C 4 configurational entropy of the molecules. 4 This, however, does not preclude the possibility that entropy related to C 3 dynamics is involved in this phase transition. Based on the observed initial condition dependence (sA or sB) of the barrier, we speculate that Br-sA is entropically stabilized by the C 3 torsion/rotation degree of freedom of the molecule and thereby determines the orthorhombic structure below the phase-transition temperature.

■ CONCLUSIONS
In conclusion, we have shown that on-the-fly MLFFs are a very powerful tool in determining the atomic structure in dynamic, entropically stabilized solids. Already during the training-byheating MD, important structural characteristics are qualitatively correct and even quantitatively useful. As a prime example, the low-temperature ordering pattern of MA molecules in MAPbX 3 perovskites, which is not uniquely resolved by diffraction experiments, was studied. We determined the most likely structure by slow-heating DFTbased MD and analyzing the librational pathways. By comparing this analysis with reported temperature-dependent lattice parameters and refined structures, we show that the ordering of the molecules (sA or sB) in orthorhombic phases of MAPbBr 3 and MAPbCl 3 is similar (sA), while in MAPbI 3 , they are differently ordered (sB). This is unexpected since the sA pattern is energetically unfavorable when considering solely the intrinsic dipole moment of the MA molecules. The sA order induces a permanent structural distortion of the PbCl 6 octahedra at low temperature, resulting in an antiferroelectric striped pattern in the crystallographic a-direction. In the higher-temperature cubic phase, this distortion is no longer observed in the ensemble average; instead, instantaneous dynamic distortions appear. No permanent distortion is The Journal of Physical Chemistry C pubs.acs.org/JPCC Article observed in the PbBr, nor PbI, octahedra, even down to the lowest simulated temperature of 80 K. We have presented indications that the sA order in low-temperature, orthorhombic MAPbBr 3 is stabilized by an entropic contribution to the free energy related to the C 3 dynamics of the MA molecules. We hope that this paper will stimulate combined experimental and MLFF studies of the structure of many other complex dynamic solids.
■ ASSOCIATED CONTENT * sı Supporting Information