Prediction and understanding of barocaloric effects in orientationally disordered materials from molecular dynamics simulations

Due to its high energy efficiency and environmental friendliness, solid-state cooling based on the barocaloric (BC) effect represents a promising alternative to traditional refrigeration technologies relying on greenhouse gases. Plastic crystals displaying orientational order-disorder solid-solid phase transitions have emerged among the most gifted materials on which to realize the full potential of BC solid-state cooling. However, a comprehensive understanding of the atomistic mechanisms on which order-disorder BC effects are sustained is still missing, and rigorous and systematic methods for quantitatively evaluating and anticipating them have not been yet established. Here, we present a computational approach for the assessment and prediction of BC effects in orientationally disordered materials that relies on atomistic molecular dynamics simulations and emulates quasi-direct calorimetric BC measurements. Remarkably, the proposed computational approach allows for a precise determination of the partial contributions to the total entropy stemming from the vibrational and molecular orientational degrees of freedom. Our BC simulation method is applied on the technologically relevant material CH$_{3}$NH$_{3}$PbI$_{3}$ (MAPI), finding giant BC isothermal entropy changes ($|\Delta S_{\rm BC}| \sim 10$ J K$^{-1}$ kg$^{-1}$) under moderate pressure shifts of $\sim 0.1$ GPa. Intriguingly, our computational analysis of MAPI reveals that changes in the vibrational degrees of freedom of the molecular cations, not their reorientational motion, have a major influence on the entropy change that accompanies the order-disorder solid-solid phase transition.


I. INTRODUCTION
Solid-state cooling represents an energy efficient and ecologically friendly solution to the environmental problems posed by conventional refrigeration technologies based on compression cycles of greenhouse gases [1][2][3][4][5][6].Upon small and moderate magnetic, electric and/or mechanical field shifts, promising caloric materials experience large adiabatic temperature variations (|∆T | ∼ 1-10 K) as a result of phase transformations entailing large isothermal entropy changes (|∆S| ∼ 10-100 J K −1 kg −1 ).Solid-state cooling relies on such caloric effects to engineer multi-step refrigeration cycles.From an applied point of view, large |∆T | and |∆S| are both necessary for substantially and efficiently removing heat from a targeted system under recurrent switching on and off of a driving field.In terms of largest |∆T | and |∆S|, mechanocaloric effects induced by uniaxial stress (elastocaloric effects) and hydrostatic pressure (barocaloric -BC-effects) emerge as particularly encouraging [2][3][4][5].
Nevertheless, despite these recent experimental advancements, a detailed atomistic understanding of colossal BC effects in plastic crystals is still missing and general theoretical approaches for assessing them have not been established.Owing to the first-order character of the involved order-disorder phase transition, a few authors have employed the Clausius-Clapeyron (CC) method in conjunction with atomistic molecular dynamics simulations to estimate ∆S BC [13,18].In the present context, however, the CC method presents several critical drawbacks: (1) the quantity that is accessed is the phase transition entropy, ∆S t , which, although related, differs from the BC descriptor ∆S BC , (2) neither ∆T BC nor the temperature span over which BC effects are operative can be directly estimated, and (3) partial entropy contributions stemming from different degrees of freedom (vibrational and orientational) cannot be determined.Alterna- ≈ 160 K, MAPI undergoes a transformation from an orthorhombic to a tetragonal phase; at T T C t ≈ 330 K, MAPI transforms into a cubic phase in which the molecular MA + ions are orientationally disordered [23].A non-negligible volume increase accompanies the order-disorder phase transition at T T C t .tively, a few researchers have resorted to oversimplified ad hoc ∆S BC analytical models with little to none predictive capability [19,20].This lack of fundamental knowledge keeps hindering the rational design of disordered materials with improved BC performances, thus delaying the deployment of commercial solid-state cooling.
In this work, we advance towards the solution of the BC modelling conundrum in plastic crystals by developing and testing a facile computational approach based on molecular dynamics simulations.Our method emulates quasi-direct calorimetry measurements [7,9,10,24] and precisely provides the vibrational and molecular orientational contributions to the entropy without resorting to ad hoc analytical models.As a case study, we apply our simulation BC approach to the technologically relevant material CH 3 NH 3 PbI 3 (MAPI, Fig. 1a), a hybrid organic-inorganic perovskite that undergoes an order-disorder solid-solid phase transition under increasing temperature (Fig. 1b) [23].In particular, we determined ∆S BC and ∆T BC under broad pressure and temperature conditions finding, for instance, a giant isothermal entropy change of 31 J K −1 kg −1 and an adiabatic temperature change of 9 K for a pressure shift of 0.2 GPa at temperatures below ambient.Intriguingly, our theoretical analysis concludes that the vibrational degrees of freedom of the molecular MA cations have a predominant role in the BC performance of MAPI, instead of the typically assumed molecular reorientations.This work establishes an insightful and predictive computational method for the estimation of BC effects in orientationally disordered systems like plastic crystals, hence it may be used to guide experiments and develop original solid-state refrigeration applications.

II. COMPUTATIONAL BC METHOD
In the present theoretical BC approach, the entropy of the low-T ordered and high-T orientationally disordered phases are determined as a function of pressure and temperature, S(p, T ), by assuming the relation: where S vib and S ori are respectively the partial entropy contributions resulting from the vibrational and molecular orientational degrees of freedom (Fig. 2), and possible vibrational-orientational molecular couplings have been disregarded.In the low-T ordered phase, S ori is null while in the high-T disordered phase is finite and positive.S vib , on the other hand, is finite and positive under any physically realizable thermodynamic conditions.Once S(p, T ) is determined for the low-T and high-T phases, the BC descriptors ∆S BC and ∆T BC are estimated just like it is done in quasi-direct calorimetric measurements [7,9,10,24] (Sec.III C).Next, we explain in detail how to calculate each entropy term appearing in Eq. (1).

A. Vibrational entropy, S vib
The vibrational density of states, ρ(ω), besides providing information on the phonon spectrum of a crystal as a function of the vibrational frequency ω, it also allows to estimate key thermodynamic quantities like the lattice free energy, F vib , and entropy, S vib .A possible way of calculating ρ(ω) from the outputs of N pT molecular dynamics simulations (N pT -MD, Methods) consists in estimating the Fourier transform of the velocity autocorrelation function (VACF) [25].
Let v i (t) be the velocity of the i-th particle expressed as a function of simulation time t.The VACF is given then by the expression: where • • • denotes statistical average in the N pT ensemble and N the number of particles in the system.Subsequently, the vibrational density of states can be estimated like [25]: In the present case, we assume the vibrational density of states to fulfill the condition: where 3N corresponds to the number of lattice vibrations with real and positively defined phonon frequencies in the system.(In the case of assuming the normalization condition ∞ 0 ρ(ω) dω = 1, an additional multiplicative factor 3N should be considered in the formulas appearing below in this subsection.) Upon determination of ρ, the vibrational free energy of the system can be straightforwardly estimated with the formula [26]: where k B is the Boltzmann's constant.Consistently, the vibrational entropy, S vib = − ∂F vib ∂T , adopts the expres-sion: It is worth noting that the dependence on pressure and temperature of the thermodynamic quantities in Eqs. ( 5)-( 6) are implicitly contained in ρ(ω), since this function is directly obtained from the outputs of atomistic N pT -MD simulations.

B. Molecular orientational entropy, Sori
For a continuous random variable x with probability density f (x), its entropy is defined as [27]: where the integral runs over all possible values of x.If x represents a microstate characterizing a particular thermodynamic macrostate, then the following Gibbs entropy can be defined for the system of interest [28]: In a orientationally disordered crystal, molecules reorient in a quasi-random fashion as shown by the time evolution of their polar (θ) and azimuthal (ϕ) angles as referred to a fixed arbitrary molecular axis (e.g., the C-N bond in the methylammonium molecule at a given time).By assuming the molecules in the crystal to be independent one from the other, one may then estimate a probability density function for their orientation, f (θ, ϕ), by considering, for instance, the atomistic trajectories generated during long N pT -MD simulations.In such a case, one can define the following three-dimensional molecular angular entropy (in which the molecules length are considered fixed): In practice, the calculation of S ang entails the construction of histograms for which the continuous polar variables are discretised, (θ, ϕ) → (θ i , ϕ j ), and bin probabilities are defined like: where ∆cos θ∆ϕ corresponds to the area of an histogram bin (assumed to be constant here).Consistently, one can then rewrite S ang in the discretised form [27]: The angular entropy defined in the equation above, however, is not necessarily equal to the entropy associated with molecular orientational disorder alone since other molecular degrees of freedom (e.g., molecular fluctuations and librations, which are of vibrational nature) are also inevitably captured by the angular probability density function f (θ, ϕ) deduced from N pT -MD simulations (Sec.III C).A practical way of getting rid of most spurious contributions to the molecular orientational entropy in the high-T disordered phase consists in offsetting S ang by the maximum value attained in the low-T ordered phase, in which molecular reorientations do not actually occur at an appreciable rate.Due to the fact that the largest molecular oscillations in the low-T ordered phase occur close to the order-disorder phase transition temperature, T t , an appropriate molecular orientational entropy then can be defined like: which is the expression that we have used throughout this work.

III. RESULTS
We used the method introduced in Sec.II to estimate BC effects in bulk CH 3 NH 3 PbI 3 (MAPI, Fig. 1), a highly promising photovoltaic material with exceptional and highly tunable optoelectronic properties [32][33][34].The reasons for this material selection are several.First, at a temperature of ≈ 330 K, MAPI un-dergoes a structural transformation from a tetragonal to a cubic phase in which the methylammonium molecular cations (MA + ) are orientationally disordered [23] (Fig. 1b).MAPI, therefore, conforms to the broad definition of orientationally disordered crystal despite being generally, and more conveniently, classified as a hybrid organic-inorganic perovskite.Second, previous experimental and theoretical works have already addressed the existence of possible caloric effects in MAPI [35,36], thus our results may be compared with those of other studies.And third, a reliable force field has been already reported and tested for MAPI in the literature [29][30][31] that allows to perform thousands-of-atoms N pT -MD simulations in a very efficient manner (Methods).The organization of this section is as follows.First, we characterize the molecular order-disorder (OD) phase transition in MAPI under broad temperature and pressure conditions by means of atomistic N pT -MD simulations performed with the force field developed by Mattoni and co-workers [29][30][31].Next, we analyse the entropy change associated with the phase transition of interest by using the well-known Clausius-Clapeyron equation [13].Finally, we explain in detail the estimation of the partial entropy contributions S vib (p, T ) and S ori (p, T ) and of the barocaloric descriptors ∆S BC and ∆T BC .
A. Order-disorder phase transition in MAPI Figure 3 shows the theoretical characterization of the OD phase transition in MAPI at pressures of 0 ≤ p ≤ 0.5 GPa as obtained from N pT -MD simulations (Meth- ods).In the high-T phase the molecular MA cations are orientationally disordered while in the low-T phase they arrange orderly in an antiferroelectric pattern [37].At zero pressure, this phase transition experimentally occurs near room temperature and is accompanied by a relative volume expansion of 0.2-0.3%[38,39].The MAPI force field employed in this work [29][30][31] qualitatively mimics such a OD phase transition, however, it does not quantitatively reproduce the experimental T t and transition volume expansion value.
Figure 3a shows the evolution of the volume per formula unit calculated across the OD MAPI phase transformation at different pressure conditions.In the N pT -MD simulations, the transition temperature is identified by a sudden increase in volume that coincides with the stabilization of molecular orientational disorder (confirmed by rapidly decaying angular autocorrelation functions estimated under the same conditions [13], Supplementary Fig. 1).In our zero-pressure simulations, we calculated a relative volume expansion of 0.8% at a transition temperature of ≈ 240 K, values that are appreciably larger and smaller, respectively, than the corresponding experimental figures.Nevertheless, since the main objective of the present study is to develop and test a solid BC computational approach, such a lack of quantitative force field accuracy should not be considered as a limiting factor (Sec. III C).
Figure 3b shows the dependence of T t estimated for MAPI under pressure.When considering the numerical uncertainties in our N pT -MD simulations, a linear p-dependence of the OD transition temperature clearly emerges.In particular, a first-order polynomial fit of the form T t (p) = 239.2+63.7p, in which the temperature and pressure are respectively expressed in units of K and GPa, is found to best reproduce the computed transition temperatures.The barocaloric coefficient dT t /dp deduced from the simulations roughly amounts to 60 K GPa −1 , which in this case is in fairly good agreement with the experimental values of 76-35 K GPa −1 reported at zeropressure conditions [35].Table I summarizes the key features of the pressure-induced OD MAPI phase transition as obtained from our N pT -MD simulations, including the phase transition entropy change, ∆S t .

B. The Clausius-Clapeyron method
The well-known Clausius-Clapeyron equation, allows for the estimation of the entropy change accompanying a phase transition, ∆S t , based on the knowledge of the slope of its p-T phase boundary and the concomitant change in volume, ∆V t .By applying this relation to the results of our N pT -MD simulations, we obtained the phase transition entropy changes summarized in Table I.
The predicted entropy change approximately amounts to 28 J K −1 kg −1 at zero pressure and steadily decreases under increasing compression (e.g., ∆S t = 19.2J K −1 kg −1 at 0.2 GPa).Such a p-induced phase transition entropy decrease is due exclusively to the ∆V t reduction originated by pressure (Table I) since, to a good approximation, the slope of the p-T phase boundary in our N pT -MD simulations is constant (Fig. 3b).
In the absence of external applied pressure, the experimental ∆S t value reported for MAPI amounts to 15.65 J K −1 kg −1 [35], which is roughly two times smaller than the one estimated here with N pT -MD simulations (Table I).This noticeable difference follows from the dp/dT t and ∆V t discrepancies found between the experiments and our calculations, which have been explained in the previous paragraphs.
As it can be appreciated in Fig. 3b, the assessed T t 's come with some numerical uncertainties (i.e., error bars in the figure) that result from the fluctuations introduced by the thermostat and barostat employed in the N pT -MD simulations along with the discreteness of the simulated p-T conditions (Methods).The presence of pretransitional effects also manifests at the highest simulated pressures, as shown by the T -dependence of the volume system obtained across the phase transition at p = 0.4 and 0.5 GPa (i.e., the volume variations become less abrupt making more difficult the identification of ∆V t , Fig. 3a).These errors, which to a certain extent are unavoidable in N pT -MD simulations, lead to inaccuracies in the estimation of dp/dT t and phase transition volume changes, and inevitably propagate to ∆S t .In the following section, we present a more refined method for the calculation of phase transition entropy changes, and entropy curves in general, based on the computational strategy introduced in Sec.II, thus overcoming the usual numerical limitations of the Clausius-Clapeyron approach [13,18].

Vibrational entropy
Figure 4 shows the vibrational density of states (VDOS) and partial organic and inorganic contributions estimated for MAPI at zero pressure and temperatures slightly above and below the simulated OD transition point (i.e., 245 and 240 K).The VDOS contribution assigned to the inorganic part, namely, the PbI 3 octahedra, is clearly dominant in the low-frequency range 0 ≤ ν ≤ 5 THz (Figs. 4a-c).This result follows from the fact that the lighter atoms, which typically vibrate at higher frequencies, entirely reside in the organic molecules.Consistently, the range of moderate and high frequencies, ν ≥ 5 THz, is mostly governed by the MA cation vibrations.
Albeit the VDOS enclosed in Figs.4a,b may seem quite similar at first glance, there are significant differences among them.In particular, the high-T disordered phase accumulates more phonon modes in the low-frequency range 0 ≤ ν ≤ 2 THz than the low-T ordered phase (Fig. 4c and Supplementary Fig. 2).This effect has a strong impact on the vibrational entropy of the system, S vib , which is significantly larger in the high-T disordered phase.
Figure 4d shows the S vib estimated at zero pressure as a function of temperature.A clear surge in the vibrational entropy is observed at the OD transition point, which amounts to ∆S vib = 21.66J K −1 kg −1 .The primary contribution to such a vibrational entropy increase comes from the molecular cations, which is equal to 13.93 J K −1 kg −1 and almost twice as large as that calculated for the inorganic part (7.73 J K −1 kg −1 ).Thus, although the low-frequency range in the VDOS is dominated by the inorganic anions, the organic MA cations have a larger influence on the vibrational entropy change associated with the order-disorder phase transition, ∆S vib .
Figures 5a,b show the angular probability density associated with the orientation of a molecular MA cation in MAPI at temperatures below and above the simulated OD transition point (i.e., averaged over all the molecules in the simulation cell).The coordinates (ϕ, θ) represent the azimuthal and polar angles for a MA molecule considering its C-N bond axis and the crystallographic axes as the rest reference system.In the low-T ordered phase, the molecules do not reorient but tightly fluctuate around their preferential orientations (represented by sharp red spots in Fig. 5a), which alternate from cation to cation forming a global antiferroelectric pattern [37].In the high-T disordered phase, on the contrary, elongated regions of roughly equal probability appear that connect the preferential molecular orientations (i.e., the blurry greenish areas in Fig. 5b).These equiprobable density regions represent the actual paths through which the MA cations reorient (a three-dimensional sketch of those rotational tracks is provided in the Supplementary Fig. 3).
It is worth noting that in the high-T disordered phase the fluctuations of the MA + cations around their equilibrium orientations are very broad, as it is deduced from the enlarged areas of maximum probability in Fig. 5b (in comparison to those shown in Fig. 5a).These fluctuations conform to molecular librations with very low vibrational frequencies due to the large effective mass associated with them.Such a collective librational dynamics appears to be the main responsible for the vibrational entropy gain identified for the high-T disordered phase in the previous Sec.III C 1.
Figure 5c shows the orientational entropy of the MA + cations, S ori , estimated as a function of temperature and pressure by employing the method introduced in Sec.II B. At temperatures below the theoretical tran-sition point, S ori is null since molecular reorientations do not occur at an appreciable rate.At temperatures above T t , on the other hand, S ori is sizeable and practically independent of temperature.The orientational entropy plateau that is rapidly attained at T > T t can be understood in terms of two different facts: (1) the paths through which the molecules reorient in the disordered phase are not significantly affected by increasing thermal excitations, and (2) although molecular reorientations occur at a higher rate at higher temperatures, the probability of finding a molecule around an equilibrium orientation or in a transient orientational state is pretty independent of temperature.It is worth noting that the behaviour described here for S ori is plainly different from that of a rigid molecular free rotor [19,20], for which the entropy monotonically increases under increasing temperature (Discussion).
The pressure dependence of the molecular orientational entropy also can be appreciated in Fig. 5c.In the high-T ordered phase, compression only moderately reduces S ori .For instance, at zero pressure it approximately amounts to 11 J K −1 kg −1 while at p = 0.5 GPa is equal to ≈ 8.5 J K −1 kg −1 .In the compression interval 0.1 ≤ p ≤ 0.4 GPa, however, S ori behaves quite irregularly and adopts values of 9-10 J K −1 kg −1 .Such fluctuations are due to the numerical uncertainties in our calculations and the mild influence of pressure on S ori .We performed a conclusive test to assess the accuracy of our S ori calculation method.At zero pressure, the entropy change associated with the order-disorder phase transition can be directly obtained from the N pT -MD simulations since the Gibbs free energy equality between the two involved phases, G O (0, T t ) = G D (0, T t ), leads to the relation ∆S t = ∆U t /T t , where ∆U t ≡ U D − U O is the internal energy difference between the disordered and ordered phases at the specified p-T conditions.Then, by considering that ∆S t = ∆S vib + ∆S ori and assuming that ∆S vib is perfectly evaluated with the method described in Sec.II A, it is possible to directly compute S ori without the need of performing integrations over polar probability maps.By proceeding so, we obtained a minute discrepancy of only 0.5 J K −1 kg −1 with respect to the S ori value obtained with the method intro-duced in Sec.II B and shown in Fig. 5c.Similar excellent agreement between the two independent S ori evaluation approaches was also found for other transition points obtained under p = 0 conditions (in this latter case, ∆S t = 1/T t [∆U t + p∆V t ]).Therefore, based on these findings we may conclude that (1) for any arbitrary p-T conditions, not necessarily ascribed to a phase transition point, our S ori estimation method based on N pT -MD simulations is fully reliable, and (2) the neglect of possible vibrational-orientational molecular couplings in Eq. ( 1) appears to be reasonable (at least, for MAPI).

Total entropy curves and estimation of BC effects
We calculated the total entropy of MAPI as a function of temperature and pressure by adding up the vibrational and orientational contributions reported in the previous two sections, namely, S(p, T ) = S vib (p, T ) + S ori (p, T ).The resulting S(p, T ) curves are very well behaved, as it is explicitly shown in Fig. 6b.In analogy to quasi-direct calorimetry experiments [7,9,10,24], the estimation of the barocaloric isothermal entropy change, ∆S BC , and adiabatic temperature change, ∆T BC , now turns out to be quite straightforward, as it is schematically shown in Fig. 6a.In particular, we have that: where T f stands for the temperature fulfilling the condition S(T f , p) = S(T, 0) (Fig. 6a).
Figure 7 shows our barocaloric ∆S BC and ∆T BC results obtained in MAPI for pressure shifts of 0 ≤ ∆p ≤ 0.5 GPa.In Fig. 7a, it is clearly appreciated that under increasing compression the maximum ∆S BC value also increases.For instance, the maximum isothermal entropy change obtained for ∆p = 0.1 GPa amounts to −28.19 J K −1 kg −1 whereas for ∆p = 0.5 GPa is equal to −38.99 J K −1 kg −1 .It is worth noting that for ∆p ≥ 0.4 GPa the growth rate of ∆S BC is found to diminish drastically.The temperature span over which BC effects remain sizeable also increases under increasing pressure shifts, changing for instance from ≈ 10 K for ∆p = 0.2 GPa to ≈ 20 K for 0.5 GPa.These isothermal entropy shift values, although probably are somewhat overestimated due to the limitations of the employed force field evidenced in Secs.III A-III B, can be regarded as giant since they are of the order of 10 J K −1 kg −1 [5].
The evolution of ∆T BC under increasing ∆p is reported in Fig. 7b.As already expected from the previous isothermal entropy results, the maximum adiabatic temperature change is also significantly enhanced under increasing compression.In particular, the ∆T BC peak changes from 4.41 K for a pressure shift of ∆p = 0.1 GPa to 22.03 K for 0.5 GPa.In this case, the size of ∆T BC is not found to saturate at the highest considered pressure.It is worth noting here that at conditions other than the phase transition points the estimated BC effects are not null (Fig. 7).These results follow from the thermal expansion properties of MAPI in the ordered and disordered phases [24], not from the OD phase transition itself.
In a previous work, Liu and Cohen investigated the possible existence of elastocaloric effects in MAPI (i.e., induced by uniaxial stress, σ) based on similar molecular dynamics simulations than reported here [36] (e.g., employed also the force field developed by Mattoni and co-workers [29][30][31]).Those authors used a direct elastocaloric estimation approach that requires of very long simulation times to attain adiabatic conditions and which cannot straightforwardly provide elastocaloric ∆S results.Interestingly, Liu and Cohen reported an elastocaloric ∆T of 10.2 K for a tensile ∆σ of 0.55 GPa [36].Although these figures are not directly comparable to our barocaloric ∆T BC results (i.e., hydrostatic pressure involves isotropic compression whereas tensile uniaxial strain conveys unidirectional stretching), they are of similar size and confirm the promising mechanocaloric potential of MAPI.
Finally, we note that the barocaloric ∆S BC results shown in Fig. 7a differ greatly from the phase entropy changes, ∆S t , estimated with the Clausius-Clapeyron method (Sec.III B).Firstly, for a pressure shift of 0.1 GPa, for example, the estimated barocaloric isothermal entropy change is larger than the corresponding phase transition entropy change by about ≈ 10%.And secondly, under increasing pressure ∆S t gets progressively reduced (Table I) whereas ∆S BC steadily increases (Fig. 7a), thus rendering opposite p-induced behaviours.Therefore, the Clausius-Clapeyron method, although it may provide an approximate order of magnitude for ∆S BC at very low pressures [13,18], is not adequate for estimating BC effects in plastic crystals.(The fact that ∆S t = ∆S BC is well established in experimental works [7][8][9][10], however, it probably needs to be more emphasized in the context of theoretical calculations.)

IV. DISCUSSION
An interesting question to answer for MAPI, or any other orientationally disordered crystal, is the following: which are the partial contributions to the order-disorder phase transition entropy change stemming from the vibrational and molecular orientational degrees of freedom?The computational approach introduced in this work can provide a quantitative answer to this question, which is shown in Fig. 8.At zero pressure, it is found that ∆S vib is practically twofold larger than ∆S ori , namely, 21.66 and 11.02 J K −1 kg −1 , respectively.(It is noted that the ∆S t = ∆S vib + ∆S ori value obtained in this case is slightly larger than the corresponding entropy change estimated with the Clausius-Clapeyron method -Table I-due to the inherent numerical inaccuracies of the latter method explained in Sec.III B.) Under increasing pressure, the vibrational entropy change gets effectively reduced whereas ∆S ori remains less affected.In spite of this behaviour, at the highest pressure considered in this study ∆S vib continues being noticeably larger than the molecular orientational entropy change, namely, 12.96 versus 8.13 J K −1 kg −1 (Fig. 8).This result appears to be consistent with recent experimental findings reported for the archetypal plastic crystal adamantane [40] and the orientationally disordered ferroelectric ammonium sulfate [41], for which it has been shown that the vibrational contributions to the order-disorder phase entropy change surpass those resulting from molecular reorientational motion.
The results presented in this article were obtained from N pT -MD simulations performed with the MAPI force field developed by Mattoni and co-workers [29][30][31].To assess the accuracy of this classical interatomic potential in providing correct molecular orientational entropies, we carried out complementary ab initio molecular dynamics (AIMD) simulations based on density functional theory [26] (DFT, Methods).Supplementary Fig. 4 shows the angular probability densities obtained for the molecular MA cations at temperatures of 250 and 400 K and zero pressure from AIMD simulations.Under the same thermodynamic conditions, the S ori values obtained from N pT -MD and AIMD simulations (without offsetting) agree remarkably well, namely, within 1%.This finding, on one hand, corroborates the reliability of the molecular orientational entropy results presented above and, on the other hand, suggests that intrinsically accurate firstprinciples methods, although being computationally very expensive, may be feasibly employed for the analysis of order-disorder phase transitions in orientationally disordered crystals.
Finally, we comment on the reasonableness of employing analytical free rotor models to describe S ori in orientationally disordered molecular crystals [19,20].From textbooks, it is well known that the entropy of a nonsymmetrical molecular free rotor can be analytically expressed as [28]: where {I i } are the three principal moments of inertia of the molecule.Since the atomic structure of the methylammonium cation is well established, we can provide here a quantitative and meaningful comparison between the two entropies, S ori and S rot , estimated for MAPI.Supplementary Fig. 5 shows the size and temperature dependence of S rot obtained for a molecular MA cation, which are significantly different from the S ori results enclosed in Fig. 5c.Firstly, S rot monotonically grows under increasing temperature (Eq.16) whereas S ori remains practically constant in the disordered phase.Secondly, at room temperature, for instance, S rot is two orders of magnitude larger than S ori , namely, ∼ 10 3 and ∼ 10 1 J K −1 kg −1 , respectively.And thirdly, Eq. ( 16) cannot reproduce any entropy dependence on pressure, even if that is mild, in contrast to S ori .The huge quantitative and qualitative differences evidenced here for S rot and S ori can be understood in terms of the facts that in real materials (i) molecules reorient only along very specific and well defined paths (Fig. 5b and Supplementary Fig. 3), not along all possible directions like a molecular free rotor does, and (ii) molecules are not in a perpetual state of orientational motion but instead integrate vibrations and fluctuations around equilibrium states with the actual reorientations.Therefore, as it has been demonstrated in this section, S rot does not seem to conform to a physically well motivated model for understanding and quantifying the orientational entropy and associated changes in plastic crystals and/or hybrid organicinorganic perovskites [19,20].

V. CONCLUSIONS
We have introduced a computational approach based on molecular dynamics simulations that allows to investigate barocaloric effects in crystals resulting from molecular order-disorder phase transitions.Both the barocaloric isothermal entropy and adiabatic temperature changes can be quantified from entropy curves expressed as a function of pressure and temperature, just like it is done in quasi-direct barocaloric measurements.The presented simulation method automatically provides the partial contributions to the crystal entropy stemming from the vibrational and molecular orientational degrees of freedom.As a case study, we applied our barocaloric computational approach to MAPI, a technologically relevant perovskite compound that experimentally undergoes an order-disorder phase transition near room temperature at normal pressure.We found giant barocaloric isothermal entropy and adiabatic temperature changes of the order of ∼ 10 J K −1 kg −1 and ∼ 10 K, respectively, for moderate pressure shifts of ∼ 0.1 GPa (although these figures probably are somewhat overestimated due to the evidenced limitations of the employed force field).Interestingly, the vibrational degrees of freedom, and in particular those ascribed to the molecular cations, were found to contribute most significantly to the phase transition entropy change at all considered pressures.Furthermore, we demonstrated that the well known analytical model for a molecular free rotor may not be adequate to physically understand and quantify the entropy changes occurring in orientationally disordered materials.We expect that the simulation approach introduced in this study will be broadly adopted by researchers working in the fields of caloric effects and/or disordered materials, thus making a potential great impact in the disciplines of energy materials and condensed matter physics.

A. Classical molecular dynamics simulations
We used the LAMMPS simulation code [42] to perform systematic classical molecular dynamics simulations in the N pT ensemble for bulk MAPI by using the force field developed by Mattoni and co-workers [29][30][31].In these simulations, the temperature was steadily increased up to a targeted value during 1 ns and subsequently the system was equilibrated during an additional 1 ns under a targeted fixed pressure.The production runs then lasted for about 100 ps with ∆t = 0.5 fs, from which the velocities of the atoms and other key quantities (e.g., the potential energy and volume of the system) were gathered.The average value of the temperature and pressure were set by using Nose-Hoover thermostats and barostats with mean fluctuations of 5 K and 0.05 GPa, respectively.The simulation box contained 3072 atoms (equivalent to 256 MAPI unit cells) and periodic boundary conditions were applied along the three Cartesian directions.The long-range electrostatic interactions were calculated by using a particle-particle particle-mesh solver to compute Ewald sums up to an accuracy of 10 −4 kcal mol −1 Å−1 in the forces.The force field uses an hybrid formulation of the Lennard-Jones and Buckingham pairwise interaction models together with long-range Coulomb and harmonic forces (the latter for the molecules).The cutoff distance for the evaluation of the potential energy was set to 10 Å.To determine the p-T phase diagram of MAPI and its barocaloric performance, we carried out comprehensive N pT -MD simulations in the temperature range 180 ≤ T ≤ 340 K, taken at intervals of 10 K, and in the pressure range 0 ≤ p ≤ 0.5 GPa, taken at intervals of 0.1 GPa.

B. Ab initio molecular dynamics simulations
First-principles calculations based on density functional theory (DFT) [26] were performed to analyze the orientational properties of MAPI at zero pressure.The DFT calculations were carried out with the VASP code [43] by following the generalized gradient approximation to the exchange-correlation energy due to Perdew et al. (PBE) [44].The projector augmented-wave (PAW) method was used to represent the ionic cores [45], and the following electronic states were considered as valence: 2s 2 2p 2 for C, 2s 2 2p 3 for N, 2s 2 3p 5 for I, 1s 1 for H and 5d 10 6s 2 6p 2 for Pb.We performed ab initio MD (AIMD) simulations for large supercells containing 384 atoms and on which periodic boundary conditions were applied along the three Cartesian directions.The plane-wave basis set was truncated at 500 eV and for integrals within the first Brillouin zone we employed Γ-point sampling.The AIMD simulations typically lasted for about 200 ps and the employed time step was equal to 1.5 fs.A Nose-Hoover thermostat was used to constrain the temperature in our AIMD calculations in which we mimicked the tetragonal (I4cm) and cubic (P m 3m) phases of MAPI.The dimensions of the system were constrained to the equilibrium volume determined at zero temperature, thus thermal expansion effects were disregarded in our AIMD simulations.

FIG. 1 .
FIG. 1. Characterization of CH3NH3PbI3 (MAPI) and experimental sequence of T -induced phase transitions.(a) Representation of the cubic unit cell of MAPI.(b) At T OT t

FIG. 2 .
FIG. 2. Sketch of T -induced molecular vibrational and orientational excitations in MAPI.Vibrations entail periodic motion of the atoms around their equilibrium positions, corresponding to fluctuations around a particular local minimum in the potential energy surface.Orientational disorder, on the other hand, appears when thermal excitations are high enough for molecules to overcome the energy barriers between different local minima in the potential energy surface.

FIG. 3 .
FIG. 3. Order-disorder phase transition in MAPI characterized by N pT -MD simulations.(a) Volume per formula unit expressed as a function of temperature at pressures 0.0 ≤ p ≤ 0.5 GPa.The phase transition is characterized by a non-negligible change in volume.The lines are guides to the eye.(b) p-T order-disorder phase boundary estimated for MAPI (red line, linear fit).The error bars result from the statistical fluctuations of the barostat and thermostat employed in the N pT -MD simulations and the discreteness of the simulated p-T conditions (Methods).

FIG. 4 .
FIG. 4. Vibrational density of states and vibrational entropy, S vib , estimated for MAPI.Total VDOS and partial contributions stemming from the PbI3 and MA motifs at zero pressure and T = 240 K (a) and 245 K (b).(c) Comparison of the VDOS obtained for low vibrational frequencies at temperatures right below and above the MAPI order-disorder phase transition.(d) Vibrational entropy calculated for MAPI as a function of temperature at zero pressure; partial contributions from the PbI3 and the MA groups are also shown.

FIG. 5 .
FIG. 5. Characterization of the angular distribution of molecular cations in MAPI and resulting orientational entropy.(a) Angular probability density for the molecular MA cations in the ordered phase at T = 200 K. High-probability regions are represented with red color whereas low-probability regions with dark blue color.(b) Angular probability density for the molecular MA cations in the disordered phase at T = 300 K. (c) Molecular orientational entropy, Sori, calculated at different temperatures and pressures.The lines are guides to the eye.

FIG. 6 .FIG. 7 .
FIG. 6. Quasi-direct estimation of BC descriptors and total entropy curves calculated for MAPI.(a) The barocaloric descriptors ∆SBC and ∆TBC can be directly assessed from the S-T curves obtained at different pressures, as it is schematically shown in the figure.(b) Entropy curves calculated for MAPI from N pT -MD simulations as a function of temperature and considering different fixed pressures.

FIG. 8 .
FIG. 8. Contributions to the total entropy change associated with the order-disorder phase transition in MAPI.∆S vib and ∆Sori represent the total vibrational and molecular orientational entropies, respectively, and it follows that ∆St = ∆S vib + ∆Sori.

TABLE I .
[29][30][31]er solid-solid phase transition in MAPI characterized by N pT -MD simulations.Results were obtained with the MAPI force field developed by Mattoni and co-workers[29][30][31].The phase-transition entropy changes, ∆St, were estimated with the Clausius-Clapeyron relation.The numerical uncertainties of the reported pressures and transition temperatures amount to 0.05 GPa and 5 K, respectively.