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

Carlos Escorihuela-Sayalero,^1,2 Luis Carlos Pardo,^1,2 Michela Romanini,^1,2 Nicolas Obrecht,³ Sophie Loehlé,³ Pol Lloveras,^1,2 Josep-Lluís Tamarit,^1,2 Claudio Cazorla^1,2 ¹Grup de Caracterització de Materials, Departament de Física, Universitat Politècnica de Catalunya, Campus Diagonal-Besòs, Av. Eduard Maristany 10-14, 08019 Barcelona, Spain
²Research Center in Multiscale Science and Engineering, Universitat Politècnica de Catalunya, Campus Diagonal-Besòs, Av. Eduard Maristany 10-14, 08019 Barcelona, Spain
³TotalEnergies OneTech, 69360 Salize, France

Abstract

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₃NH₃PbI₃ (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 ( $|\Delta T|\sim 1$ – $10$ K) as a result of phase transformations entailing large isothermal entropy changes ( $|\Delta S|\sim 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 $|\Delta T|$ and $|\Delta 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 $|\Delta T|$ and $|\Delta S|$ , mechanocaloric effects induced by uniaxial stress (elastocaloric effects) and hydrostatic pressure (barocaloric –BC– effects) emerge as particularly encouraging [2; 3; 4; 5].

Colossal BC effects, $|\Delta S_{\rm BC}|\geq 100$ J K^-1 kg^-1, and $|\Delta T_{\rm BC}|\geq 10$ K driven by pressure shifts of the order of $0.1$ GPa have been recently measured in several families of materials displaying orientational order-disorder solid-solid phase transitions [7; 8; 9; 10; 11; 12; 13; 14; 15; 16], thus achieving a major breakthrough in the field of solid-state cooling. Plastic crystals like neopentane derivatives [7; 8; 9], adamantane derivatives [10; 11], carboranes [12] and closo-borates [13], to cite some examples, conform the most representative family of disordered materials on which such a BC revolution has been realized. Colossal BC effects in plastic crystals have been intuitively rationalized in terms of large entropy changes predominantly originated by molecular orientational disorder stabilized under increasing temperature [17; 18; 19; 20] (as it is commonly assumed in molecular liquids [21; 22]).

Refer to caption — Figure 1: Characterization of CH₃NH₃PbI₃ (MAPI) and experimental sequence of $T$ -induced phase transitions. (a) Representation of the cubic unit cell of MAPI. (b) At $T_{t}^{OT}\approx 160$ K, MAPI undergoes a transformation from an orthorhombic to a tetragonal phase; at $T_{t}^{TC}\approx 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}^{TC}$ .

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 $\Delta S_{\rm 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, $\Delta S_{t}$ , which, although related, differs from the BC descriptor $\Delta S_{\rm BC}$ , (2) neither $\Delta T_{\rm 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. Alternatively, a few researchers have resorted to oversimplified ad hoc $\Delta S_{\rm 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₃NH₃PbI₃ (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 $\Delta S_{\rm BC}$ and $\Delta T_{\rm 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:

S(p,T)=S_{\rm vib}(p,T)+S_{\rm ori}(p,T)~{},

(1)

where $S_{\rm vib}$ and $S_{\rm 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_{\rm ori}$ is null while in the high- $T$ disordered phase is finite and positive. $S_{\rm 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 $\Delta S_{\rm BC}$ and $\Delta T_{\rm BC}$ are estimated just like it is done in quasi-direct calorimetric measurements [7; 9; 10; 24] (Sec. III.3). Next, we explain in detail how to calculate each entropy term appearing in Eq. (1).

II.1 Vibrational entropy, $S_{\rm vib}$

The vibrational density of states, $\rho(\omega)$ , besides providing information on the phonon spectrum of a crystal as a function of the vibrational frequency $\omega$ , it also allows to estimate key thermodynamic quantities like the lattice free energy, $F_{\rm vib}$ , and entropy, $S_{\rm vib}$ . A possible way of calculating $\rho(\omega)$ from the outputs of $NpT$ molecular dynamics simulations ( $NpT$ –MD, Methods) consists in estimating the Fourier transform of the velocity autocorrelation function (VACF) [25].

Let $\mathbf{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:

\text{VACF}(t)=\left\langle\mathbf{v}(0)\cdot\mathbf{v}(t)\right\rangle=\frac{1}{N}\sum_{i}^{N}\mathbf{v_{i}}(0)\cdot\mathbf{v_{i}}(t)~{},

(2)

where $\langle\cdots\rangle$ denotes statistical average in the $NpT$ ensemble and $N$ the number of particles in the system. Subsequently, the vibrational density of states can be estimated like [25]:

\rho(\omega)=\int_{0}^{\infty}\text{VACF}(t)~{}e^{i\omega t}dt~{}.

(3)

In the present case, we assume the vibrational density of states to fulfill the condition:

\int_{0}^{\infty}\rho(\omega)~{}d\omega=3N~{},

(4)

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 $\int_{0}^{\infty}\rho(\omega)~{}d\omega=1$ , an additional multiplicative factor $3N$ should be considered in the formulas appearing below in this subsection.)

Upon determination of $\rho$ , the vibrational free energy of the system can be straightforwardly estimated with the formula [26]:

F_{\rm vib}(p,T)=k_{B}T\int_{0}^{\infty}\ln{\left[2\sinh{\left(\frac{\hbar\omega}{2k_{B}T}\right)}\right]}\rho(\omega)d\omega~{},

(5)

where $k_{B}$ is the Boltzmann’s constant. Consistently, the vibrational entropy, $S_{\rm vib}=-\frac{\partial F_{\rm vib}}{\partial T}$ , adopts the expression:

	$\displaystyle S_{\rm vib}(p,T)$	$\displaystyle=$	$\displaystyle-\int_{0}^{\infty}k_{B}\ln{\left[2\sinh{\left(\frac{\hbar\omega}{2k_{B}T}\right)}\right]}\rho(\omega)d\omega+$		(6)
			$\displaystyle\int_{0}^{\infty}\frac{\hbar\omega}{2T}\tanh^{-1}{\left(\frac{\hbar\omega}{2k_{B}T}\right)}\rho(\omega)d\omega~{}.$		(6)

It is worth noting that the dependence on pressure and temperature of the thermodynamic quantities in Eqs. (5)–(6) are implicitly contained in $\rho(\omega)$ , since this function is directly obtained from the outputs of atomistic $NpT$ –MD simulations.

II.2 Molecular orientational entropy, $S_{\rm ori}$

For a continuous random variable $x$ with probability density $f(x)$ , its entropy is defined as [27]:

S=-\int_{X}f(x)\log{f(x)}~{}dx~{},

(7)

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]:

S_{G}=-k_{B}\int_{X}f(x)\log{f(x)}~{}dx~{}.

(8)

In a orientationally disordered crystal, molecules reorient in a quasi-random fashion as shown by the time evolution of their polar ( $\theta$ ) and azimuthal ( $\varphi$ ) 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(\theta,\varphi)$ , by considering, for instance, the atomistic trajectories generated during long $NpT$ –MD simulations. In such a case, one can define the following three-dimensional molecular angular entropy (in which the molecules length are considered fixed):

S_{\rm ang}=-k_{B}\int_{\Theta,\Phi}f(\theta,\varphi)\log{f(\theta,\varphi)}~{}d\cos\theta d\varphi~{}.

(9)

In practice, the calculation of $S_{\rm ang}$ entails the construction of histograms for which the continuous polar variables are discretised, $(\theta,\varphi)\to(\theta_{i},\varphi_{j})$ , and bin probabilities are defined like:

p_{i}(\theta_{i},\varphi_{j})\approx f(\theta_{i},\varphi_{j})\Delta{\cos\theta}\Delta\varphi~{},

(10)

where $\Delta{\cos\theta}\Delta\varphi$ corresponds to the area of an histogram bin (assumed to be constant here). Consistently, one can then rewrite $S_{\rm ang}$ in the discretised form [27]:

S_{\rm ang}=-k_{B}\sum_{i}p_{i}\log{p_{i}}+k_{B}\log{\left(\Delta{\cos\theta}\Delta\varphi\right)}~{}.

(11)

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(\theta,\varphi)$ deduced from $NpT$ –MD simulations (Sec. III.3). 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_{\rm 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:

S_{\rm ori}(p,T)=\left\{\begin{array}[]{lr}S_{\rm ang}(p,T)-S_{\rm ang}(p,T_{t})&T\geq T_{t}\\ 0&T<T_{t}\end{array}\right.,

(12)

which is the expression that we have used throughout this work.

$p$	$T_{t}$	$dT_{t}/dp$	$\Delta V_{t}/Z$	$\Delta V_{t}/V$	$\Delta S_{t}$
(GPa)	(K)	(K GPa^-1)	(Å³)	(%)	(J K^-1 kg^-1)
0.0	240	63.7	1.86	0.8	28.4
0.1	245	63.7	1.66	0.7	25.4
0.2	250	63.7	1.37	0.6	19.2
0.3	260	63.7	1.04	0.3	15.9
0.4	265	63.7	0.90	0.4	13.8
0.5	270	63.7	0.94	0.4	14.3

Table 1: Order-disorder solid-solid phase transition in MAPI characterized by

NpT

–MD simulations. Results were obtained with the MAPI force field developed by Mattoni and co-workers [29; 30; 31]. The phase-transition entropy changes,

\Delta S_{t}

, 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.

III Results

We used the method introduced in Sec. II to estimate BC effects in bulk CH₃NH₃PbI₃ (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 $\approx 330$ K, MAPI undergoes 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 $NpT$ -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 $NpT$ –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_{\rm vib}(p,T)$ and $S_{\rm ori}(p,T)$ and of the barocaloric descriptors $\Delta S_{\rm BC}$ and $\Delta T_{\rm BC}$ .

III.1 Order-disorder phase transition in MAPI

Figure 3 shows the theoretical characterization of the OD phase transition in MAPI at pressures of $0\leq p\leq 0.5$ GPa as obtained from $NpT$ -MD simulations (Methods). 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 $NpT$ -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 $\approx 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.3).

Figure 3b shows the dependence of $T_{t}$ estimated for MAPI under pressure. When considering the numerical uncertainties in our $NpT$ -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 zero-pressure conditions [35]. Table 1 summarizes the key features of the pressure-induced OD MAPI phase transition as obtained from our $NpT$ -MD simulations, including the phase transition entropy change, $\Delta S_{t}$ .

III.2 The Clausius-Clapeyron method

The well-known Clausius-Clapeyron equation,

\Delta S_{t}=\left(\frac{dT_{t}}{dp}\right)^{-1}\Delta V_{t}~{},

(13)

allows for the estimation of the entropy change accompanying a phase transition, $\Delta S_{t}$ , based on the knowledge of the slope of its $p$ – $T$ phase boundary and the concomitant change in volume, $\Delta V_{t}$ . By applying this relation to the results of our $NpT$ -MD simulations, we obtained the phase transition entropy changes summarized in Table 1. The predicted entropy change approximately amounts to $28$ J K^-1 kg^-1 at zero pressure and steadily decreases under increasing compression (e.g., $\Delta S_{t}=19.2$ J K^-1 kg^-1 at $0.2$ GPa). Such a $p$ -induced phase transition entropy decrease is due exclusively to the $\Delta V_{t}$ reduction originated by pressure (Table 1) since, to a good approximation, the slope of the $p$ – $T$ phase boundary in our $NpT$ -MD simulations is constant (Fig. 3b).

In the absence of external applied pressure, the experimental $\Delta 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 $NpT$ -MD simulations (Table 1). This noticeable difference follows from the $dp/dT_{t}$ and $\Delta 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 $NpT$ -MD simulations along with the discreteness of the simulated $p$ – $T$ conditions (Methods). The presence of pre-transitional 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 $\Delta V_{t}$ , Fig. 3a). These errors, which to a certain extent are unavoidable in $NpT$ -MD simulations, lead to inaccuracies in the estimation of $dp/dT_{t}$ and phase transition volume changes, and inevitably propagate to $\Delta 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].

III.3 Estimation of BC effects in MAPI

III.3.1 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₃ octahedra, is clearly dominant in the low-frequency range $0\leq\nu\leq 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, $\nu\geq 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\leq\nu\leq 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_{\rm vib}$ , which is significantly larger in the high- $T$ disordered phase.

Figure 4d shows the $S_{\rm 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 $\Delta S_{\rm vib}=21.66$ J 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, $\Delta S_{\rm vib}$ .

III.3.2 Molecular orientational entropy

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 $(\varphi,\theta)$ 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.3.1.

Figure 5c shows the orientational entropy of the MA⁺ cations, $S_{\rm ori}$ , estimated as a function of temperature and pressure by employing the method introduced in Sec. II.2. At temperatures below the theoretical transition point, $S_{\rm ori}$ is null since molecular reorientations do not occur at an appreciable rate. At temperatures above $T_{t}$ , on the other hand, $S_{\rm 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_{\rm 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_{\rm 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 $\approx 8.5$ J K^-1 kg^-1. In the compression interval $0.1\leq p\leq 0.4$ GPa, however, $S_{\rm 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_{\rm ori}$ .

We performed a conclusive test to assess the accuracy of our $S_{\rm ori}$ calculation method. At zero pressure, the entropy change associated with the order-disorder phase transition can be directly obtained from the $NpT$ -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 $\Delta S_{t}=\Delta U_{t}/T_{t}$ , where $\Delta U_{t}\equiv 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 $\Delta S_{t}=\Delta S_{\rm vib}+\Delta S_{\rm ori}$ and assuming that $\Delta S_{\rm vib}$ is perfectly evaluated with the method described in Sec. II.1, it is possible to directly compute $S_{\rm 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_{\rm ori}$ value obtained with the method introduced in Sec. II.2 and shown in Fig. 5c. Similar excellent agreement between the two independent $S_{\rm ori}$ evaluation approaches was also found for other transition points obtained under $p\neq 0$ conditions (in this latter case, $\Delta S_{t}=1/T_{t}\left[\Delta U_{t}+p\Delta V_{t}\right]$ ). 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_{\rm ori}$ estimation method based on $NpT$ -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).

III.3.3 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_{\rm vib}(p,T)+S_{\rm 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, $\Delta S_{\rm BC}$ , and adiabatic temperature change, $\Delta T_{\rm BC}$ , now turns out to be quite straightforward, as it is schematically shown in Fig. 6a. In particular, we have that:

	$\displaystyle\Delta S_{\rm BC}(T,0\rightarrow p)$	$\displaystyle=$	$\displaystyle S(T,p)-S(T,0)$		(14)
	$\displaystyle\Delta T_{\rm BC}(T,0\rightarrow p)$	$\displaystyle=$	$\displaystyle T_{\rm f}(S,p)-T(S,0)~{},$		(15)

where $T_{\rm f}$ stands for the temperature fulfilling the condition $S(T_{\rm f},p)=S(T,0)$ (Fig. 6a).

Figure 7 shows our barocaloric $\Delta S_{\rm BC}$ and $\Delta T_{\rm BC}$ results obtained in MAPI for pressure shifts of $0\leq\Delta p\leq 0.5$ GPa. In Fig. 7a, it is clearly appreciated that under increasing compression the maximum $\Delta S_{\rm BC}$ value also increases. For instance, the maximum isothermal entropy change obtained for $\Delta p=0.1$ GPa amounts to $-28.19$ J K^-1 kg^-1 whereas for $\Delta p=0.5$ GPa is equal to $-38.99$ J K^-1 kg^-1. It is worth noting that for $\Delta p\geq 0.4$ GPa the growth rate of $\Delta S_{\rm 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 $\approx 10$ K for $\Delta p=0.2$ GPa to $\approx 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.1–III.2, can be regarded as giant since they are of the order of $10$ J K^-1 kg^-1 [5].

The evolution of $\Delta T_{\rm BC}$ under increasing $\Delta 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 $\Delta T_{\rm BC}$ peak changes from $4.41$ K for a pressure shift of $\Delta p=0.1$ GPa to $22.03$ K for $0.5$ GPa. In this case, the size of $\Delta T_{\rm 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, $\sigma$ ) 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 $\Delta S$ results. Interestingly, Liu and Cohen reported an elastocaloric $\Delta T$ of $10.2$ K for a tensile $\Delta\sigma$ of $0.55$ GPa [36]. Although these figures are not directly comparable to our barocaloric $\Delta T_{\rm 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 $\Delta S_{\rm BC}$ results shown in Fig. 7a differ greatly from the phase entropy changes, $\Delta S_{t}$ , estimated with the Clausius-Clapeyron method (Sec. III.2). 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 $\approx 10$ %. And secondly, under increasing pressure $\Delta S_{t}$ gets progressively reduced (Table 1) whereas $\Delta S_{\rm 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 $\Delta S_{\rm BC}$ at very low pressures [13; 18], is not adequate for estimating BC effects in plastic crystals. (The fact that $\Delta S_{t}\neq\Delta S_{\rm 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 $\Delta S_{\rm vib}$ is practically twofold larger than $\Delta S_{\rm ori}$ , namely, $21.66$ and $11.02$ J K^-1 kg^-1, respectively. (It is noted that the $\Delta S_{t}=\Delta S_{\rm vib}+\Delta S_{\rm ori}$ value obtained in this case is slightly larger than the corresponding entropy change estimated with the Clausius-Clapeyron method –Table 1– due to the inherent numerical inaccuracies of the latter method explained in Sec. III.2.) Under increasing pressure, the vibrational entropy change gets effectively reduced whereas $\Delta S_{\rm ori}$ remains less affected. In spite of this behaviour, at the highest pressure considered in this study $\Delta S_{\rm 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 $NpT$ –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_{\rm ori}$ values obtained from $NpT$ –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 first-principles 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_{\rm ori}$ in orientationally disordered molecular crystals [19; 20]. From textbooks, it is well known that the entropy of a non-symmetrical molecular free rotor can be analytically expressed as [28]:

\frac{S_{\rm rot}}{k_{B}}=\frac{3}{2}+\ln{\left[\frac{8\pi^{2}}{h^{3}}\left(2\pi k_{B}T\right)^{3/2}\left(I_{1}I_{2}I_{3}\right)^{1/2}\right]}~{},

(16)

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_{\rm ori}$ and $S_{\rm rot}$ , estimated for MAPI. Supplementary Fig. 5 shows the size and temperature dependence of $S_{\rm rot}$ obtained for a molecular MA cation, which are significantly different from the $S_{\rm ori}$ results enclosed in Fig. 5c. Firstly, $S_{\rm rot}$ monotonically grows under increasing temperature (Eq. 16) whereas $S_{\rm ori}$ remains practically constant in the disordered phase. Secondly, at room temperature, for instance, $S_{\rm rot}$ is two orders of magnitude larger than $S_{\rm ori}$ , namely, $\sim 10^{3}$ and $\sim 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_{\rm ori}$ . The huge quantitative and qualitative differences evidenced here for $S_{\rm rot}$ and $S_{\rm 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_{\rm 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 organic-inorganic 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 $\sim 10$ J K^-1 kg^-1 and $\sim 10$ K, respectively, for moderate pressure shifts of $\sim 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.

Methods

V.1 Classical molecular dynamics simulations

We used the LAMMPS simulation code [42] to perform systematic classical molecular dynamics simulations in the $NpT$ 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 $\Delta 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 $NpT$ -MD simulations in the temperature range $180\leq T\leq 340$ K, taken at intervals of $10$ K, and in the pressure range $0\leq p\leq 0.5$ GPa, taken at intervals of $0.1$ GPa.

V.2 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 $\Gamma$ -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 ( $Pm\bar{3}m$ ) 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.

References

Mañosa et al. [2013] L. Mañosa, A. Planes, and M. Acet, Journal of Materials Chemistry A 1, 4925 (2013).
Mañosa et al. [2010] L. Mañosa, D. González-Alonso, A. Planes, E. Bonnot, M. Barrio, J.-L. Tamarit, S. Aksoy, and M. Acet, Nature Materials 9, 478 (2010).
Mañosa and Planes [2017] L. Mañosa and A. Planes, Advanced Materials 29, 1603607 (2017).
Cazorla [2019a] C. Cazorla, Applied Physics Reviews 6, 041316 (2019a).
Lloveras and Tamarit [2021] P. Lloveras and J.-L. Tamarit, MRS Energy & Sustainability 8, 3 (2021).
Hou et al. [2022] H. Hou, S. Qian, and I. Takeuchi, Nature Reviews Materials 7, 633 (2022).
Lloveras et al. [2019] P. Lloveras, A. Aznar, M. Barrio, P. Negrier, C. Popescu, A. Planes, L. Mañosa, E. Stern-Taulats, A. Avramenko, N. D. Mathur, X. Moya, and J.-L. Tamarit, Nature Communications 10, 1803 (2019).
Li et al. [2019] B. Li, Y. Kawakita, S. Ohira-Kawamura, T. Sugahara, H. Wang, J. Wang, Y. Chen, S. I. Kawaguchi, S. Kawaguchi, K. Ohara, K. Li, D. Yu, R. Mole, T. Hattori, T. Kikuchi, S.-i. Yano, Z. Zhang, Z. Zhang, W. Ren, S. Lin, O. Sakata, K. Nakajima, and Z. Zhang, Nature 567, 506 (2019).
Aznar et al. [2020] A. Aznar, P. Lloveras, M. Barrio, P. Negrier, A. Planes, L. Mañosa, N. D. Mathur, X. Moya, and J.-L. Tamarit, Journal of Materials Chemistry A 8, 639 (2020).
Aznar et al. [2021] A. Aznar, P. Negrier, A. Planes, L. Mañosa, E. Stern-Taulats, X. Moya, M. Barrio, J.-L. Tamarit, and P. Lloveras, Applied Materials Today 23, 101023 (2021).
Salvatori et al. [2022] A. Salvatori, P. Negrier, A. Aznar, M. Barrio, J. L. Tamarit, and P. Lloveras, APL Materials 10, 111117 (2022).
Zhang et al. [2022] K. Zhang, R. Song, J. Qi, Z. Zhang, Z. Zhang, C. Yu, K. Li, Z. Zhang, and B. Li, Advanced Functional Materials 32, 2112622 (2022).
Sau et al. [2021] K. Sau, T. Ikeshoji, S. Takagi, S.-i. Orimo, D. Errandonea, D. Chu, and C. Cazorla, Scientific Reports 11, 11915 (2021).
Imamura et al. [2020] W. Imamura, É. O. Usuda, L. S. Paixão, N. M. Bom, A. M. Gomes, and A. M. G. Carvalho, Chinese Journal of Polymer Science 38, 999 (2020).
Li et al. [2021] J. Li, M. Barrio, D. J. Dunstan, R. Dixey, X. Lou, J.-L. Tamarit, A. E. Phillips, and P. Lloveras, Advanced Functional Materials 31, 2105154 (2021).
Seo et al. [2022] J. Seo, R. D. McGillicuddy, A. H. Slavney, S. Zhang, R. Ukani, A. A. Yakovenko, S.-L. Zheng, and J. A. Mason, Nature Communications 13, 2536 (2022).
Cazorla [2019b] C. Cazorla, Nature 567, 470 (2019b).
Li et al. [2022] F. Li, M. Li, C. Niu, and H. Wang, Applied Physics Letters 120, 073902 (2022).
Li et al. [2020] F. B. Li, M. Li, X. Xu, Z. C. Yang, H. Xu, C. K. Jia, K. Li, J. He, B. Li, and H. Wang, Nature Communications 11, 4190 (2020).
de Oliveira [2023] N. de Oliveira, Acta Materialia 246, 118657 (2023).
Henao et al. [2016] A. Henao, S. Busch, E. Guàrdia, J. L. Tamarit, and L. C. Pardo, Phys. Chem. Chem. Phys. 18, 19420 (2016).
Pardo et al. [2015] L. C. Pardo, A. Henao, and A. Vispa, J. Non-Cryst Solids 407, 220 (2015).
Shahrokhi et al. [2022] S. Shahrokhi, M. Dubajic, Z.-Z. Dai, S. Bhattacharyya, R. A. Mole, K. C. Rule, M. Bhadbhade, R. Tian, N. Mussakhanuly, X. Guan, Y. Yin, M. P. Nielsen, L. Hu, C.-H. Lin, S. L. Y. Chang, D. Wang, I. V. Kabakova, G. Conibeer, S. Bremner, X.-G. Li, C. Cazorla, and T. Wu, Small 18, 2200847 (2022).
Aznar et al. [2017] A. Aznar, P. Lloveras, M. Romanini, M. Barrio, J.-L. Tamarit, C. Cazorla, D. Errandonea, N. D. Mathur, A. Planes, X. Moya, and L. Mañosa, Nature Communications 8, 1851 (2017).
Sagotra et al. [2019] A. K. Sagotra, D. Chu, and C. Cazorla, Phys. Rev. Mater. 3, 035405 (2019).
Cazorla and Boronat [2017] C. Cazorla and J. Boronat, Rev. Mod. Phys. 89, 035003 (2017).
Cover and Thomas [2006] T. M. Cover and J. A. Thomas, Elements of Information Theory (Wiley Series in Telecommunications and Signal Processing) (Wiley-Interscience, USA, 2006).
Pathria [1972] R. Pathria, Statistical Mechanics, International series of monographs in natural philosophy (Elsevier Science & Technology Books, 1972).
Mattoni et al. [2015] A. Mattoni, A. Filippetti, M. I. Saba, and P. Delugas, The Journal of Physical Chemistry C 119, 17421 (2015).
Mattoni et al. [2016] A. Mattoni, A. Filippetti, and C. Caddeo, Journal of Physics: Condensed Matter 29, 043001 (2016).
Caddeo et al. [2020] C. Caddeo, A. Filippetti, and A. Mattoni, Nano Energy 67, 104162 (2020).
Hu et al. [2022] L. Hu, L. Duan, Y. Yao, W. Chen, Z. Zhou, C. Cazorla, C.-H. Lin, X. Guan, X. Geng, F. Wang, T. Wan, S. Wu, S. Cheong, R. D. Tilley, S. Liu, J. Yuan, D. Chu, T. Wu, and S. Huang, Advanced Science 9, 2102258 (2022).
Ren et al. [2020] L. Ren, Y. Wang, M. Wang, S. Wang, Y. Zhao, C. Cazorla, C. Chen, T. Wu, and K. Jin, The Journal of Physical Chemistry Letters 11, 2577 (2020).
Hellmann et al. [2020] T. Hellmann, C. Das, T. Abzieher, J. A. Schwenzer, M. Wussler, R. Dachauer, U. W. Paetzold, W. Jaegermann, and T. Mayer, Advanced Energy Materials 10, 2002129 (2020).
Bermúdez-García et al. [2017] J. M. Bermúdez-García, M. Sánchez-Andújar, and M. A. Señarís-Rodríguez, The Journal of Physical Chemistry Letters 8, 4419 (2017).
Liu and Cohen [2016] S. Liu and R. E. Cohen, The Journal of Physical Chemistry C 120, 17274 (2016).
Frost and Walsh [2016] J. M. Frost and A. Walsh, Accounts of Chemical Research 49, 528 (2016).
Whitfield et al. [2016] P. S. Whitfield, N. Herron, W. E. Guise, K. Page, Y. Q. Cheng, I. Milas, and M. K. Crawford, Scientific Reports 6, 35685 (2016).
Patru et al. [2021] R. E. Patru, H. Khassaf, I. Pasuk, M. Botea, L. Trupina, C.-P. Ganea, L. Pintilie, and I. Pintilie, Materials 14, 4215 (2021).
Meijer et al. [2023] B. E. Meijer, R. J. C. Dixey, F. Demmel, R. Perry, H. C. Walker, and A. E. Phillips, Phys. Chem. Chem. Phys. 25, 9282 (2023).
Yuan et al. [2022] S. Yuan, B. E. Meijer, G. Cai, R. J. C. Dixey, F. Demmel, M. T. Dove, J. Liu, H. Y. Playford, H. C. Walker, and A. E. Phillips, Advanced Functional Materials 32, 2207717 (2022).
Plimpton [1995] S. Plimpton, Journal of Computational Physics 117, 1 (1995).
Kresse and Hafner [1993] G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
Perdew et al. [1996] J. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
Blöchl [1994] P. Blöchl, Physical Review B 50, 17953 (1994).

Acknowledgements

The authors acknowledge TotalEnergies RD High Performance Computing Center in Houston and Production High Performance Computing Center in Pau for the use of our HPC. This work was also supported by the MINECO Projects PID2020-112975GB-I00 and TED2021-130265B-C22 (Spain) and by the DGU Project 2021SGR-00343 (Catalonia). C.C. acknowledges financial support from the Spanish Ministry of Science, Innovation and Universities under the “Ramón y Cajal” fellowship RYC2018-024947-I. C.E.S. acknowledges financial support from the Spanish Ministry of Science, Innovation and Universities under the ”Becas Margarita Salas para la formación de doctores jóvenes” fellowship 2021UPC-MS-67395. Additional computational support was provided by the Red Española de Supercomputación (RES) under the grants FI-2022-1-0006, FI-2022-2-0003 and FI-2022-3-0014.

Conflict of Interest

The authors declare no conflict of interest.

Data Availability Statement

The data that support the findings of this study are available from the corresponding authors upon reasonable request.