Colossal Reversible Barocaloric Effects in a Plastic Crystal Mediated by Lattice Vibrations and Ion Diffusion

Abstract Solid‐state methods for cooling and heating promise a sustainable alternative to current compression cycles of greenhouse gases and inefficient fuel‐burning heaters. Barocaloric effects (BCE) driven by hydrostatic pressure (p) are especially encouraging in terms of large adiabatic temperature changes (|ΔT| ≈ 10 K) and isothermal entropy changes (|ΔS| ≈ 100 J K−1 kg−1). However, BCE typically require large pressure shifts due to irreversibility issues, and sizeable |ΔT| and |ΔS| seldom are realized in a same material. Here, the existence of colossal and reversible BCE in LiCB11H12 is demonstrated near its order‐disorder phase transition at ≈380 K. Specifically, for Δp ≈ 0.23 (0.10) GPa, |ΔS rev| = 280 (200) J K−1 kg−1 and |ΔT rev| = 32 (10) K are measured, which individually rival with state‐of‐the‐art BCE figures. Furthermore, pressure shifts of the order of 0.1 GPa yield huge reversible barocaloric strengths of ≈2 J K−1 kg−1 MPa−1. Molecular dynamics simulations are performed to quantify the role of lattice vibrations, molecular reorientations, and ion diffusion on the disclosed BCE. Interestingly, lattice vibrations are found to contribute the most to |ΔS| while the diffusion of lithium ions, despite adding up only slightly to the entropy change, is crucial in enabling the molecular order–disorder phase transition.


I. INTRODUCTION
Solid-state methods for cooling and heating are energy efficient and ecologically friendly techniques with potential for solving the environmental problems posed by conventional refrigeration and heat pump technologies relying on compression cycles of greenhouse gases and inefficient traditional fuel-burning heaters [1].Under moderate magnetic, electric or mechanical field variations, auspicious 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 ) [2,3].Solidstate cooling and heat pumping capitalize on such caloric effects for engineering refrigeration and heating cycles.From a practical point of view, large and reversible |∆T | and |∆S| are both necessary for achieving rapid and efficient devices under recursive application and removal of the driving fields.In terms of largest |∆T | and |∆S|, mechanocaloric effects induced by uniaxial stress (elastocaloric effects) and hydrostatic pressure (barocaloric effects -BCE-) are among the most promising [4][5][6].
In spite of these recent developments, finding barocaloric materials with well-balanced and suitable features for developing thermal applications, e.g., |∆T rev | ≥ FIG. 1. Sketch of the order-disorder phase transition occurring in LCBH upon increasing temperature.(a) Ballstick representation of the low-T ordered (O) and high-T disordered (D) phases.Lithium, carbon, boron and hydrogen atoms are represented with red, brown, green and blue spheres, respectively.In the high-T phase, the Li + cations diffuse throughout the crystalline matrix while the [CB11H12] − anions reorient disorderly [20]; the volume increases significantly during the Tinduced phase transition.(b) Outline of the order-disorder phase transition in terms of Gibbs free energies.The red dotted lines represent internal energies and the blue solid lines Gibbs free energies; Tt denotes the phase transition temperature.20 K and |∆S rev | ≥ 100 J K −1 kg −1 driven by ∆p 0.1 GPa, is proving extremely difficult.From the hundred of barocaloric materials known to date [6], to the best of our knowledge only four fulfill the conditions specified above, namely, the spin-crossover complex Fe 3 (bntrz) 6 (tcnset) 6 (|∆T rev | = 35 K and |∆S rev | = 120 J K −1 kg −1 for ∆p = 0.26 GPa) [21], the layered hybrid perovskite [C 10 H 21 NH 3 ] 2 MnCl 4 (|∆T rev | = 27 K and |∆S rev | = 250 J K −1 kg −1 for ∆p = 0.19 GPa) [13,14], the plastic crystal 1-Br-adamantane (|∆T rev | = 20 K and |∆S rev | = 120 J K −1 kg −1 for ∆p = 0.10 GPa) [10], and the elastomer acetoxy silicone (|∆T rev | = 22 K and |∆S rev | = 182 J K −1 kg −1 for ∆p = 0.17 GPa) [12].Moreover, studies addressing a fundamental and quantitative understanding of the atomistic mechanisms that bring on such colossal BCE are very scarce [22][23][24][25], thus hindering the rational design of disordered materials with enhanced barocaloric performances.
In this work, we experimentally and theoretically demonstrate the existence of colossal and reversible BCE in the monocarba-closo-dodecaborate LiCB 11 H 12 (LCBH) near its order-disorder phase transition occurring at T t ≈ 380 K [20].LCBH is a well-known solid electrolyte in which at temperatures above T t the lithium cations are highly mobile and the molecular anions [CB 11 H 12 ] − reorient disorderly [26,27] (Fig. 1); thus, LCBH combines phase-transition features of both plastic crystals and superionic compounds, two families of materials for which colossal and giant BCE, respectively, have been previously reported [7][8][9]16].In par-  Very interestingly, the contribution of the lattice vibrations to ∆S was found to be the dominant at all pressures, instead of the typically assumed one resulting from molecular reorientational motion [22][23][24].Our results provide new valuable insights into the physical behavior and functionality of plastic crystals and suggest that colossal BCE similar to those reported here for LCBH could also exist in other akin closo-borate materials like NaCB 11 B 12 [20,26], KCB 11 B 12 [28], and LiCB 9 H 10 [29,30].

A. LiCB11H12 general properties
In a recent X-ray powder diffraction study [20], it has been shown that at room temperature LiCB 11 H 12 (LCBH) presents an ordered orthorhombic structure (space group P ca2 1 ) in which the Li + cations reside near trigonal-planar sites surrounded by molecular [CB 11 H 12 ] − anions arranged in a cubic sublattice.An order-disorder phase transition occurs at T t ≈ 380 K that stabilizes a disordered phase in which the Li + cations are highly mobile and the molecular anions present fast reorientational motion (Fig. 1a).At normal pressure, the lithium ion conductivity measured just above T t exceeds values of 0.1 S cm −1 [20] and the reorientational motion of the molecular anions can reach frequencies of 10 11 s −1 [20,31].Meanwhile, the T -induced order-disorder phase transition is accompanied by a huge volume increase of the order of ≈ 10% [31] that, based on the Clausius-Clapeyron (CC) equation ∆S t = ∆V t dp dT , suggests great barocaloric potential.
The described order-disorder phase transition can be qualitatively understood in terms of the Gibbs free energy difference between the high-T disordered (D) and ∆T ( K ) low-T ordered (O) phases, ∆G ≡ G D − G O (Fig. 1b).This free energy difference consists of an internal energy (∆E), entropy (−T ∆S), and volume (p∆V ) terms.The internal energy remains more or less constant during the phase transition while the volume term disfavors the stabilization of the disordered phase since ∆V > 0. Thus, the LCBH order-disorder phase transition appears to be governed by the change in entropy, ∆S, which in view of the ion conductivity and molecular reorientational frequency measured above T t should be fairly large.

B. Experimental barocaloric results
Conventional X-ray powder diffraction experiments performed at normal pressure and under varying temperature confirmed the expected structures of the low-T and high-T phases (orthorhombic and cubic symmetry, respectively).Pattern matching analysis of the obtained data yielded the temperature-dependent volume of LCBH (see Fig. 2a), which shows a huge ≈ 13% relative volume increase at the endothermic transition corresponding to ∆V ≈ 12 • 10 −5 m 3 kg −1 .
High-pressure differential thermal analysis (HP-DTA) was carried out in the pressure interval 0 ≤ p ≤ 0.23 GPa (Fig. 2b).At pressures below ≈ 0.13 GPa, a single peak in the heat flow was measured corresponding to the aforementioned orthorhombic (ordered phase, II) ↔ cubic (disordered phase, I) first-order phase transition.At pressure above ≈ 0.13 GPa, the HP-DTA signals exhibit two peaks thus indicating the appearance of a new phase that we label here as III (high-pressure enantiotropy).To the best of our knowledge, phase III has not been previously reported in the literature and its specific crystalline structure remains unknown since we did not resolve it.Interestingly, a broad peak was previously detected in differential scanning calorimetry experiments [20] that hints at the stabilization of phase III.
Transition temperatures were determined from the maximum of the HP-DTA peaks (Fig. 2c) which allowed to estimate an upper threshold for the triple point at ≈ (425 K, 0.13 GPa) given the width of the peaks obtained under the chosen scanning rate.Considering only the data measured near atmospheric pressure, the pressure dependence of the II→I transition was determined to be dT dp ≈ 420 K GPa −1 , which slightly decreases under in- Compendium of experimentally measured reversible BCE.The size of the symbols represents the reversible barocaloric strength defined as the ratio of |∆Srev| by the corresponding pressure change ∆p.Material names are indicated near each symbol or in the right side of the panel.NPG: neopentylglycol; PG: pentaglycerine; NPA: Neopentyl alcohol; o-carb: orthocarborane; m-carb: metacarborane; p-carb: paracarborane; 1-Br-ada: 1-Bromoadamantane; 1-Cl-ada: 1-Chloroadamantane; 1ada-ol: 1-adamantanol; 2ada-ol: 2-adamantanol; 2m2ada-ol: 2-methyl-2-adamantanol; ASR: Acetoxy Silicone Rubber.Numerical details and references can be found in the Supplementary Table S1.
creasing pressure due to the small convexity of the coexistence line.For the II→III and III→I transitions, linear fits to the obtained coexistence lines yielded dT dp ≈ 135 K GPa −1 and dT dp ≈ 310 K GPa −1 , respectively.Phase transition entropy changes were calculated via integration of the 1 T dQ dT function after baseline subtraction.As it was already expected, the ∆S II→I values associated to the LCBH order-disorder phase transition are noticeably large, namely, ≈ 208 J K −1 kg −1 (Fig. 2d).By plugging the measured dT dp and ∆S II→I values at atmospheric pressure in the CC equation we obtain ∆V CC ≈ 9 • 10 −5 m 3 kg −1 , which is in reasonable agreement with the ∆V determined directly in the experiments.
Above p ≈ 0.13 GPa, due to the overlapping between the II↔III and III↔I peaks, the contribution associated to each phase transition was decided at the inflection point of the cumulative entropy change function T T1 1 T dQ dT dT .∆S t remains practically constant from atmospheric pressure all the way up to the triple point.At p 0.13 GPa, we obtained ∆S II→I ≈ ∆S II→III +∆S III→I , as it is required by the condition of thermodynamic equilibrium, while above the triple point ∆S II→III ≈ ∆S III→I .Splitting of the II→I phase transition into II→III and III→I might be associated to the decoupling of the dif-fusive and orientational degrees of freedom right at the stabilization of the high-T phase, although further investigations are necessary for a more conclusive assessment of phase III.
HP-DTA measurements along with experimental differential scanning calorimetry (Supplementary Fig. S1), heat capacity (Supplementary Fig. S2) and theoretical equations of state V (T, p) (i.e., obtained from molecular dynamics simulations, Sec.II C) were used to determine the isobaric entropy curves S(T, p) (Supplementary Fig. S3), from which the BC effects can be directly calculated (Methods).Figures 3a,b  Operation of solid-state cooling and heating devices requires cyclic application and removal of the driving fields, for which reversible caloric effects, |∆S rev | and |∆T rev |, must be considered.By reversible caloric effects we mean acquitted of phase transition hysteresis effects [9].The obtained results are shown in Figs.3c,d Figure 4 compares most of the experimental |∆S rev | and |∆T rev | reported thus far in the literature for barocaloric materials.Additionally, the size of the symbols therein account for the materials BC strength, which is defined as the ratio of |∆S rev | by the corresponding pressure shift ∆p.The best performing barocaloric materials, therefore, should appear in the top right side of the panel and with the largest possible symbol area.Each material has been represented with one or two points that best illustrate their overall barocaloric performance, while for LCBH we have selected a set of barocaloric measurements.
Although   ,b).At zero pressure, we estimated a huge volume increase of about 11% at the theoretical transition temperature T t ≈ 400 K, along with the order parameter changes ∆D Li = 1.13 • 10 −6 cm 2 s −1 and ∆λ CBH = 0.33 • 10 11 s −1 .It was found that the pressure dependence of the transition temperature could be precisely reproduced by the second-order polynomial curve T t (p) = 412 + 438p − 610p 2 (red line in the inset of Fig. 5b), in which the temperature and pressure are expressed in units of K and GPa, respectively.
Partial contributions to the entropy change accompanying the order-disorder phase transition in LCBH expressed as a function of pressure.Entropy changes stem from the vibrational, ∆S vib , molecular orientational, ∆Sori, and cation diffusive, ∆S diff , degrees of freedom.
Results were obtained from comprehensive molecular dynamics simulations and Gibbs free energy calculations (Methods).
The slight dT dp decrease under increasing compression is consistent with the p-induced reduction of the transition volume change since ∆S t is roughly independent of pressure, in agreement with our experiments.It is worth noting that phase-transition hysteresis effects cannot be reproduced by the equilibrium MD approach employed in this study [25].
The LCBH p-T phase diagram obtained from MD simulations (Fig. 5b) is in quantitative good agreement with the experiments performed below the triple point found at ≈ 0.13 GPa (Fig. 2b), although the transition temperatures are slightly overestimated by theory.For example, at zero pressure and p = 0.10 GPa the MD simulations yielded T t = 410 ± 15 and 440 ± 15 K (Fig. 5), respectively, to be compared with the corresponding experimental values 390 ± 10 and 410 ± 10 K (Fig. 2b).The agreement between the predicted and measured volumes for the ordered and disordered phases at zero pressure is also notable, finding only small relative discrepancies of ∼ 1% for the low-T phase (Figs.2a and 5a).Meanwhile, the triple point observed in the experiments was not reproduced by the MD simulations.It is worth noting, however, that under p = 0 conditions and close to T t we observed pre-transitional effects in our simulations consisting of few slowly diffusing Li ions in the ordered phase (Supplementary Fig. S4).
Figures 5c,d show the theoretical barocaloric |∆S| and |∆T | deduced from the entropy curves S(p, T ) enclosed in Fig. 5b, which were obtained from data generated in the MD simulations.The agreement between these theoretical results and the corresponding experimental values is remarkably good for pressures below the experi-mental triple point.For example, for a pressure shift of 0.10 GPa we estimated an isothermal entropy change of 227 J K −1 kg −1 and an adiabatic temperature change of 32 K from the MD simulations, to be compared with the corresponding experimental values 250 J K −1 kg −1 and 24 K (Fig. 3a,b).In view of such a notable agreement, we characterized with MD simulations the contributions to the phase transition entropy change stemming from the vibrational, molecular orientational and cation diffusive degrees of freedom, a highly valuable atomistic insight that in principle cannot be obtained from the experiments.
Figures 6a,b reveal synchronized surges in D Li and λ CBH at the order-disorder phase transition points.Thus, both ion diffusion and molecular anion orientational disorder (Figs.6e,f) contribute to the transition entropy change and barocaloric effects disclosed in LCBH.Nevertheless, there is a third possible source of entropy in the crystal which is related to the lattice vibrations, S vib (Supplementary Figs.S5-S6).Figures 6c,d show examples of the cumulative S vib function expressed as a function of the vibrational phonon energy, calculated for LCBH in the ordered and disordered phases at zero pressure and evaluated for each atomic species.Therein, it is appreciated that the largest contribution to the S vib difference between the order and disordered phases comes from the B atoms (followed by hydrogen).This outcome can be rationalized in terms of the relative great abundance of this species in LCBH (≈ 45%) and its larger mass as compared to that of H atoms (10 times heavier): B ions have a predominant weight on the low-frequency vibrational modes (Fig. 6c-d) that most significantly contribute to S vib near ambient temperature.
Figure 7 shows the relative contributions of the vibrational, molecular orientational and ion diffusion degrees of freedom to the phase transition entropy change estimated at different pressures with MD simulations.Interestingly, in all the analyzed cases the largest contribution stems from changes in the lattice vibrations, ∆S vib , followed by the molecular reorientations, ∆S ori , and finally ion diffusion, ∆S diff .For example, at zero pressure the vibrational, molecular orientational and ion diffusive degrees of freedom respectively contribute in ≈ 48, 32 and 20% to ∆S t .The entropy preeminence of the lattice vibrations can be rationalized in terms of (1) the huge volume expansion accompanying the order-disorder phase transition (∼ 10%, Fig. 5a), which further curtails the frequency of the low-energy phonon bands in the disordered phase (Supplementary Fig. S5), and (2) the intensification and amplitude broadening of the molecular libration modes in the disordered phase (inferred from the angular probability density variations around the equilibrium positions in Figs.6e-f).These outcomes are highly valuable and insightful since thus far molecular reorientations were thought to be the primary source of entropy variation in plastic crystals undergoing orderdisorder phase transitions [22][23][24].
The vibrational and orientational entropy changes re-main more or less constant for pressures ≤ 0.1 GPa, whereas ∆S diff significantly decreases under compression.For instance, at 0.1 GPa the diffusive degrees of freedom contribute to ∆S t in less than 4%.These outcomes can be understood in terms of the small fraction of diffusive ions in LCBH (i.e., one Li atom per formula unit) and the marked decline in D Li induced by pressure (Fig. 6a).The appearance of pre-transitional effects in our MD simulations, specially under p = 0 conditions (Supplementary Fig. S4), also contributes to the noticeable ∆S diff drop caused by compression.Nonetheless, it is worth noting that despite the relative minuteness of ∆S diff , cation disorder was found to play a critical role on triggering molecular orientational disorder, which by contrast contributes very significantly to ∆S t .In particular, we conducted constrained MD runs in which we fixed the positions of the lithium ions so that they could not diffuse.It was found then that molecular orientational disorder only emerged at temperatures well above 550 K (Supplementary Fig. S7).Therefore, it can be concluded that cation disorder crucially assists on the realization of colossal BCE through the order-disorder phase transition, a characteristic trait that differentiates LCBH from other molecular plastic crystals bearing also great barocaloric promise.

III. CONCLUSIONS
Colossal barocaloric effects (BCE) driven by pressure shifts of the order of 0.10 GPa were experimentally and theoretically disclosed in bulk LiCB 11 H 12 (LCBH), a compound that at high temperatures presents disorder features characteristic of both plastic crystals and superionic materials, namely, molecular reorientational motion and ion diffusion.Reversible peaks of |∆S rev | = 280 J K −1 kg −1 and |∆T rev | = 32 K were experimentally measured around 400 K for a pressure shift of 0.23 GPa, yielding huge and reversible barocaloric strengths of ≈ 2 J K −1 kg −1 MPa −1 over tens of degrees intervals.Likewise, for a smaller pressure shift of 0.10 GPa we obtained very promising values of |∆S rev | = 200 J K −1 kg −1 and |∆T rev | = 10 K.These results place LCBH among the best-known barocaloric materials in terms of huge and reversible isothermal entropy and adiabatic temperature changes, two quantities that rarely are found simultaneously in a same material.
Atomistic molecular dynamics simulations yielded theoretical |∆S| and |∆T | in very good agreement with the experimental values, and allowed to quantify the importance of vibrational, molecular orientational, and ion diffusive degrees of freedom on the disclosed colossal BCE.It was found that the contribution to the phase transition entropy change stemming from the lattice vibrations was the largest, followed by that of molecular reorientations and both being much superior than the entropy associated to lithium diffusion alone.Nevertheless, cationic disorder was found to have a critical influence on the stabilization of orientational disorder thus, in spite of its small contribution to ∆S t , lithium diffusion appears to be essential for the emergence of colossal BCE in bulk LCBH.These results are of high significance since reveal the preeminence of the vibrational degrees of freedom in the phase transition entropy change of a plastic crystal, and demonstrate atomistic BCE mechanisms other than molecular reorientational disorder (i.e., lattice vibrations and ion diffusion).
LCBH belongs to the family of closo-borate materials, a promising class of solid electrolytes for all-solid-state batteries.Examples of akin compounds that have been already synthesized in the laboratory and tested for electrochemical energy storage applications are NaCB 11 H 12 [20,26], KCB 11 H 12 [28], and LiCB 9 H 10 [29,30].Colossal BCE could also exist in these materials and in other similar compounds harboring both ion diffusion and molecular orientational disorder at or near room temperature.Thus, the present combined experimentaltheoretical study opens new horizons in solid-state cooling and heating and advances knowledge in the realization of colossal BCE in plastic crystals.

Experimental techniques
Materials synthesis.LiCB 11 H 12 was obtained by drying the hydrated compound LiCB 11 H 12 •xH 2 O (Katchem, Ltd.) under vacuum (< 5 × 10 −4 Pa) at 160 • C for 12 h.X-ray powder diffraction.High-resolution X-ray powder diffraction measurements were performed using the Debye-Scherrer geometry and transmission mode with a horizontally mounted cylindrical position-sensitive INEL detector (CPS-120).Monochromatic Cu-Kα 1 radiation was selected by means of a curved germanium monochromator.Temperature-dependent measurements were performed using a liquid nitrogen 700 series Oxford Cryostream Cooler.Powder samples were introduced into 0.5 mm diameter Lindemann capillaries.Volume was obtained by pattern matching procedure.Quasi-direct barocaloric measurements.A Q100 thermal analyzer (TA Instruments) was used to perform differential scanning calorimetry experiments at atmospheric pressure with ∼ 10 mg of sample hermetically encapsulated in Aluminum pans (Supplementary Fig. S1).The standard mode (at 3, 5 and 10 K min −1 ) was used to determine the transition properties whereas the modulated mode (isothermal conditions, modulation amplitude 1 • C, modulation period 120 s) was used to measure the heat capacity in each phase (Supplementary Fig. S2).
Pressure-dependent calorimetry was performed with a custom-built high-pressure differential thermal analyzer (from Irimo, Bellota Herramientas S.A.) that uses Bridgman thermocouples as thermal sensors.The nominal operational pressure range is from atmospheric to 0.3 GPa and the temperature range is from room temperature up to 473 K. Heating ramps were performed at 3 K min −1 using a resistive heater whereas cooling were carried out at ∼ −2 K min −1 by an air stream.A few hundreds of mg of LiCB 11 H 12 were mixed with an inert perfluorinated fluid (Galden Bioblock Scientist) to remove air and sealed within tin capsules.The pressure-transmitting fluid was Therm240 (Lauda).
Isobaric entropy functions S(T, p) were determined with respect to a reference temperature T 0 below the transition using the method explained in Ref. [32] (Supplementary Fig. S3).The procedure is based on the following thermodynamic equation: where dQ dT is the heat flow in temperature due to the firstorder phase transition measured by pressure-dependent calorimetry.
In each phase, C p is the corresponding heat capacity and was considered independent of pressure as indicated by the approximately linear behavior of volume with temperature obtained in the two phases from MD simulations (Fig. 5) along with the thermodynamic equation: In the transition region C p was calculated as an average weighted according to the fraction of each phase.
To take into account the dependence of the transition region with pressure, the overall C p function at atmospheric pressure obtained in each phase and across the transition was extrapolated to higher temperatures according to the experimental value of dT dp ∆p, where ∆p is the pressure change applied in each particular case.Experimental measurement of C p at atmospheric pressure and the calculated curves at different pressures are shown in Supplementary Fig. S2.
The pressure dependence of S(T, p) was evaluated using the thermodynamic equation: where p 0 was selected equal to p atm = 1 bar.Here, we make use of the approximation ∂V ∂T T,p ∂V ∂T T,p0 , which is reasonable based on the ∂V ∂T T,p data obtained from the MD simulations (Fig. 5).
Once the entropy function S(T, p) was determined for both heating and cooling runs independently (Supplementary Fig. S3), BC effects obtained upon first application or removal of the field were calculated as: ∆S(T, p 0 → p 1 ) = S(T, p 1 ) − S(T, p 0 ) and ( 4) ∆T (T s , p 0 → p 1 ) = T (S, p 1 ) − T s (T, p 0 ) , where T s is the starting temperature of the heating/cooling process.
Here, it must be considered that for materials with dT dp > 0 BC effects on compression (p 0 = p atm , p 1 > p atm ) and decompression (p 0 > p atm , p 1 = p atm ) are calculated from S(T, p) functions obtained on cooling and heating, respectively [9].In turn, BC effects obtained reversibly on cyclic compression-decompression processes were calculated from the S(T, p) curves obtained on heating at atmospheric pressure and cooling at high pressure.

Simulation techniques
Molecular dynamics simulations.Force-field based molecular dynamics (MD) simulations were performed using a previously reported interatomic potential for LCBH [31].This force field is a combination of Coulomb-Buckingham (CB), harmonic bond, and angle-type potentials, namely: where q i denotes the charge of the ion labeled i, 0 the vacuum permittivity, A ij and ρ the short-range repulsive energy and length scales for the pairs of atoms ij, and C ij the corresponding dispersion interaction coefficient.r 0 and θ 0 are an equilibrium bond distance and angle, respectively, and k r and k θ the spring constants of the harmonic bond and angle potentials.The numerical value of these potential parameters can be found in the Supplementary Table S2.
We performed N pT -MD simulations in the temperature range 325 ≤ T ≤ 525 K at intervals of 12.5 K, and pressure range 0 ≤ p ≤ 0.15 GPa at intervals of 0.025 GPa.The temperature and pressure in the system were controlled with thermostating and barostating techniques, in which some dynamic variables are coupled with the particle velocities and simulation box dimensions.The simulation supercell comprised a total of 6400 atoms.A time step of 0.5 fs was employed for integration of the atomic forces along with the velocity Verlet algorithm.A typical N pT -MD run lasted for about 2 ns and the atomic trajectories were stored at intervals of 500 fs.Detailed analyses and statistical time averages were performed over the last 1 ns of such simulations.To guarantee proper convergence of the estimated thermodynamic properties, in few instances longer simulation times of 10 ns were carried out.Periodic boundary conditions were applied along the three Cartesian directions and the Ewald summation technique was used for evaluation of the long-range Coulomb interactions with a short-range cut-off distance of 13 Å.All the N pT -MD simulations were carried out with the LAMMPS software package [33].
Density functional theory and ab initio molecular dynamics simulations.First-principles calculations based on density functional theory (DFT) were performed to analyze the energy, structural and vibrational properties of bulk LCBH.The DFT calculations were carried out with the VASP code [34] by following the generalized gradient approximation to the exchange-correlation energy due to Perdew et al. (PBE) [35].The projector augmented-wave method was used to represent the ionic cores [36], and the electronic states 1s-2s Li, 2s-2p C, 2s-2p B and 1s H were considered as valence.Wave functions were represented in a plane-wave basis truncated at 650 eV.By using these parameters and dense k-point grids for Brillouin zone integration, the resulting energies were converged to within 1 meV per formula unit.In the geometry relaxations, a tolerance of 0.005 eV Å−1 was imposed in the atomic forces.
Ab initio molecular dynamics (AIMD) simulations based on DFT were carried out to assess the reliability of the interatomic potential model employed in the MD simulations on the description of the vibrational degrees of freedom of bulk LCBH (Supplementary Fig. S6).The AIMD simulations were performed in the canonical ensemble (N, V, T ) considering constant number of particles, volume and temperature.The constrained volumes were equal to the equilibrium volumes determined at zero temperature, an approximation that has been shown to be reasonable at moderate temperatures [37].The temperature in the AIMD simulations was kept fluctuating around a set-point value by using Nose-Hoover thermostats.A large simulation box containing 800 atoms was employed in all the simulations, and periodic boundary conditions were applied along the three Cartesian directions.Newton's equations of motion were integrated by using the customary Verlet's algorithm and a time-step length of δt = 10 −3 ps.Γ-point sampling for integration within the first Brillouin zone was employed in all the AIMD simulations.The AIMD simulations comprised long simulation times of ≈ 200 ps and temperatures in the range 200 ≤ T ≤ 500 K.
Estimation of key quantities with MD simulations.The mean square displacement of the lithium ions was estimated with the formula [38]: where r i (t j ) is the position of the migrating ion i at time t j (= j • δt), τ represents a lag time, n τ = τ /δt, N ion is the total number of mobile ions, and N step the total number of time steps.The maximum n τ was chosen equal to N step /2, hence we could accumulate enough statistics to reduce significantly the fluctuations in MSD Li (τ ) at large τ 's.The diffusion coefficient of lithium ions was calculated with the Einstein's relation: by performing linear fits to the averaged MSD Li values calculated at long τ .The angular autocorrelation function of the molecular [CB 11 H 12 ] − anions was estimated using the expression [25]: where r is a unitary vector connecting the center of mass of each closoborane unit with one of its edges and • • • denotes statistical average in the (N, p, T ) ensemble considering all the molecular anions.This autocorrelation function typically decays as ∝ exp [−λ CBH • τ ], where the parameter λ CBH represents a characteristic reorientational frequency.For significant anion reorientational motion, that is, large λ CBH , the φ CBH function decreases rapidly to zero with time.
The temperature dependence of the lithium diffusion coefficient was assumed to follow an Arrhenius law at any pressure of the form: where D 0 and E a are parameters that depend on p and k B represents the Boltzmann constant.The reorientational frequency of closoborane units, λ CBH , was assumed to follow a similar dependence on temperature.The entropy of each phase was calculated as a function of temperature and pressure, S(p, T ), by fully considering the vibrational, molecular orientational and ion diffusive degrees of freedom: S(p, T ) = S vib (p, T ) + S ori (p, T ) + S diff (p, T ) .(11) In the low-T phase, S ori and S diff are null while in the high-T phase are finite and positive.
The vibrational density of states (VDOS), g(ω), was calculated via the Fourier transform of the velocityvelocity autocorrelation function obtained directly from the N pT -MD simulations, namely: where v i (t) represents the velocity of the atom labeled i at time t, and • • • denotes statistical average in the (N, p, T ) ensemble.The vibrational entropy was subsequently estimated with the formula [39]: where ĝ(ω) is the normalized vibrational density of states ( ∞ 0 ĝ(ω)dω = 3N ion ) and the dependence on pressure (and also temperature) is implicitly contained in ĝ(ω).
The ion diffusive entropy difference was estimated at the phase transition points via equalization of the Gibbs free energies of the low-T (O) and high-T (D) phases, namely, G D (p, T t ) = G O (p, T t ), thus leading to the expression: where ∆X ≡ X D − X O and E represents the internal energy of the system.For any pressure, ∆S diff was assumed to be constant at temperatures T t ≤ T .

FIG. 2 .
FIG. 2. Experimental phase diagram of bulk LCBH and corresponding phase transition entropy changes.(a) Volume per formula unit measured as a function of temperature at normal pressure.(b) Isobaric heat flow data expressed as a function of applied pressure and temperature; data collected during heating (cooling) are represented in the positive (negative) y-axis.(c) Pressure and temperature phase diagram; transition temperatures are determined from the peaks in panel (b).(d) Phase transition entropy changes as a function of pressure and transition path.∆St remains practically constant from atmospheric pressure all the way up to the triple point.At p 0.13 GPa, ∆SII→I ≈ ∆SII→III + ∆SIII→I, while above the triple point ∆SII→III ≈ ∆SIII→I.Straight lines at pressures above the triple point are linear fits to ∆SII→III + ∆SIII→I.
show representative isothermal entropy changes, |∆S|, and adiabatic temperature changes, |∆T |, obtained upon the first application and removal of the driving pressure shift.It is worth noticing that a small ∆p ≈ 0.03 GPa already produced colossal values of |∆S| = 100 J K −1 kg −1 and |∆T | = 8 K, and similarly ∆p ≈ 0.08 GPa yielded |∆S| = 250 J K −1 kg −1 and |∆T | = 16 K.For the largest pressure shift considered in this study, namely, ∆p ≈ 0.23 GPa, the resulting |∆S| and |∆T | amount to 300 J K −1 kg −1 and 40 K, respectively.

p
FIG. 5. Colossal BCE estimated for bulk LCBH with MD simulations.(a) Volume change per formula unit across the phase transition expressed as a function of temperature and pressure.(b) Total entropy curves expressed as a function of pressure and temperature.Inset: theoretically calculated p-T phase diagram.(c) Isothermal entropy and (d) adiabatic temperature changes expressed as a function of temperature and pressure.Results were obtained from N pT -MD simulations.
. Colossal |∆S rev | were already obtained for a minimum pressure shift of ≈ 0.08 GPa.For instance, under a moderate pressure change of ≈ 0.10 GPa LCBH renders |∆S rev | = 200 J K −1 kg −1 and |∆T rev | = 10 K.Meanwhile, for the largest pressure shift considered in this study we measured outstanding values of |∆S rev | = 280 J K −1 kg −1 and |∆T rev | = 32 K.

FIG. 6 .
FIG. 6. Atomistic insights into the order-disorder phase transition in LCBH from MD simulations.(a) Lithium ion diffusion coefficient, DLi.(b) Anionic reorientational frequency, λCBH.Solid lines correspond to Arrhenius law fits.(c)-(d) Cumulative function of the vibrational entropy as a function of the phonon energy and atomic species, calculated for the ordered (T = 400 K) and disordered (T = 412 K) phases at zero pressure.Dashed lines indicate analogous asymptotic values reached in the ordered phase.(e)-(f) Angular probability density function estimated for the molecular (CB11H12) − anions calculated in the ordered (T = 350 K) and disordered (T = 550 K) phases at zero pressure, expressed as a function of the polar (θ) and azimuthal (φ) angles.Dark and bright areas represent low and high probability regions, respectively.
t j + τ ) − r i (t j )| 2 , [21] is not the best performing material in terms of a single quality, it displays an unprecedentedly well-balanced and accomplished barocaloric portfolio consisting of colossal |∆S rev |, large |∆T rev | and large BC strength obtained under moderate pressure shifts of the order of 0.10 GPa.For instance, in terms of largest |∆S rev | the plastic crystal neopentylglycol (NPG) emerges as the clear winner since it holds a gigantic value of ≈ 400 J K −1 kg −1[9]; however, as regards |∆T rev | the same material becomes a poor contestant in the presence of LCBH (that is, ≈ 8 K versus 32 K).Likewise, the |∆T rev | record holder, namely, the spincrossover complex Fe 3 (bntrz) 6 (tcnset) 6[21], presents |∆S rev | and BC strength values that roughly are halves of the LCBH maxima (for instance, ≈ 120 J K −1 kg −1 versus 280 J K −1 kg −1 ).Therefore, LCBH can be deemed as one of the most thorough and promising barocaloric materials reported to date owing to its unique parity between sizable |∆S rev | and |∆T rev | obtained under moderate pressure shifts.