# Atomistic simulation of barocaloric effects

Claudio Cazorla<EMAIL_ADDRESS>Departament de Física, Universitat
Politècnica de Catalunya, Campus Nord B4-B5, Barcelona 08034, Spain

###### Abstract

Due to critical environmental issues there is a pressing need to switch from
current refrigeration methods based on compression of greenhouse gases to
novel solid-state cooling technologies. Solid-state cooling capitalizes on the
thermal response of materials to external fields named “caloric effect”. The
barocaloric (BC) effect driven by hydrostatic pressure is particularly
promising from a technological point of view since typically presents larger
cooling potential than other caloric variants (e.g., magnetocaloric and
electrocaloric effects driven by magnetic and electric fields, respectively).
Atomistic simulation of BC effects represents an efficient and physically
insightful strategy for advancing solid-state cooling by complementing, and in
some cases even guiding, experiments. Atomistic simulation of BC effects
involves approaches ranging from computationally inexpensive force fields to
computationally very demanding, but quantitatively accurate, first-principles
methods. Here, we survey several methods and strategies involved in atomistic
simulation of BC effects like the quasi-harmonic approximation and
direct/quasi-direct estimation approaches. The Review finalizes with a
collection of case studies in which some of these methods were employed to
simulate and predict original BC effects.

##  I Introduction

When the isothermal entropy changes associated with fluctuations in the
electric, magnetic and structural properties of materials are large, $|\Delta
S_{T}|\sim 10$–$100$ kJK-1kg-1, these may render sizable caloric effects
(e.g., adiabatic temperature changes of $|\Delta T_{S}|\sim 1$–$10$ K). Solid-
state cooling is an environmentally friendly, highly energy-efficient, and
highly scalable technology that exploits caloric effects in materials for
refrigeration purposes [cazorla19b, ; moya14, ; cazorla17b, ; cazorla16, ;
scott11, ; pecharsky97, ; bonnot08, ; shirsath20, ]. By applying external
fields on caloric materials it is possible to achieve reversible temperature
shifts that can be integrated in cooling cycles and do not involve the use of
greenhouse gases. Materials undergoing abrupt structural phase transitions
under small hydrostatic pressure are specially well suited for solid-state
cooling applications based on the barocaloric (BC) effect [cazorla19b, ;
moya14, ; cazorla17b, ].

The magnitude of BC effects can be estimated accurately with theoretical and
atomistic simulation methods. Computer simulations can be used to rationalize
the origins of experimentally observed BC phenomena since the materials and
pressure-induced phase transitions of interest can be accessed at the atomic
scale under controlled conditions. Moreover, from a resources point of view
computer simulations are inexpensive. For instance, by using open-source
software and modest computational resources it is possible, in some cases, to
simulate complex BC effects and assess the magnitude of the accompanying
isothermal entropy and adiabatic temperature shifts. Thus, modeling of BC
effects can be done systematically in order to complement, and in some cases
also guide, the experiments.

The reliability of computer simulations, however, depends strongly on the
simplifications made on the adopted structural and interatomic interaction
models. Typically, increasing the reliability of the structural models comes
at the expense of reducing the accuracy in the description of the interatomic
forces (due to practical computational limitations). For instance, if the
simulation approach to be employed is accurate first-principles methods most
calculations are likely to be performed at zero temperature by considering
perfectly ordered atomic structures. Such simulation conditions obviously
differ from the actual experimental conditions. On the other hand, to simulate
BC effects directly at finite temperatures for realistic systems containing
crystalline defects and/or other inhomogeneities one should use
computationally inexpensive interatomic modeling techniques, which may suffer
from transferability issues and in general have modest predictive power.
Fortunately, there are well-established simulation strategies that allow to
achieve a suitable balance between computational accuracy and model
reliability and can be used to obtain physically meaningful results.

In this Review, we survey general computational techniques based on atomistic
simulation methods that can be used to theoretically estimate and predict BC
effects. We start by briefly describing genuine first-principles methods, like
density functional theory, and other practical approaches that allow to
describe the interactions between atoms in materials, namely, classical
interatomic and machine learning potentials. Next, several computational
techniques that allow to describe phase transitions and estimate thermodynamic
properties of materials (e.g., free energies and entropies) are summarized
(e.g., the quasi-harmonic approximation and thermodynamic integration
methods). Subsequently, different simulation strategies that can be employed
for indirect, quasi-direct and direct evaluation of BC effects are succinctly
explained. We finalize the Review by providing some representative examples in
which atomistic simulation methods have been used to predict original BC
effects in families of materials that are technologically relevant, namely,
fast-ion conductors [cazorla16, ; sagotra17, ; sagotra18, ; min20, ], complex
hydrides [sau21, ] and multiferroic perovskite oxides [menendez20, ;
machado21, ]. This Review is addressed to both experienced computational
materials scientists with an interest on caloric effects and solid-state
cooling as well as to novice researchers wanting to learn modeling techniques
of phase transitions and free-energy methods.

##  II Modeling of atomic interactions

The fundamentals of density functional theory (DFT), a first-principles
approach widely used in condensed-matter and materials science, and classical
interatomic potentials are briefly reviewed next. The latter interatomic
interactions modeling approach is approximate and typically rely on the
results of accurate first-principles methods (also called _ab initio_) like
DFT. On the other hand, the computational cost of DFT methods is several
orders of magnitude higher than that of classical interatomic potentials,
hence in many occasions first-principles methods may not be directly employed
in the study of BC effects. Machine learning potentials are also succinctly
described since they constitute a modern atomistic modeling technique in the
middle way between DFT and classical potential approaches both in terms of
accuracy and computational expense.

###  II.1 $Ab$ $initio$ methods

In solids, the dynamics of the electrons and nuclei can be decoupled to a good
approximation because their respective masses differ by several orders of
magnitude. The wave function of the corresponding many-electron system,
$\Psi({\bf r}_{1},{\bf r}_{2},...,{\bf r}_{N})$, therefore can be determined
by solving the Schrödinger equation corresponding to the non-relativistic
Born-Oppenheimer Hamiltonian:

$\displaystyle
H=-\frac{1}{2}\sum_{i}{\bf\nabla}^{2}_{i}-\sum_{I}\sum_{i}\frac{Z_{I}}{|{\bf
R}_{I}-{\bf r}_{i}|}$ $\displaystyle+\frac{1}{2}\sum_{i}\sum_{j\neq
i}\frac{1}{|{\bf r}_{i}-{\bf r}_{j}|}~{},$ (1)

where $Z_{I}$ are the nuclear charges, ${\bf r}_{i}$ the positions of the
electrons, and ${\bf R}_{I}$ the positions of the nuclei (considered to be
fixed). In real materials, $\Psi$ is a complex mathematical function that in
most cases is unknown. At the heart of any first-principles method is to find
a good approximation for $\Psi$, or an equivalent quantity (e.g., the
electronic density), that is manageable enough to perform calculations.
Examples of _ab initio_ methods include density functional theory (DFT),
Møller-Plesset perturbation theory (MP2), the coupled-cluster method with
single, double and perturbative triple excitations [CCSD(T)], and quantum
Monte Carlo (QMC), to cite just a few. Among these techniques, DFT methods are
frequently applied to the study of materials and pressure-induced phase
transitions hence for this reason we summarise their foundations in what
follows.

In 1965, Kohn and Sham developed a pioneering theory to effectively calculate
the energy and properties of many-electron systems without the need of
explicitly knowing $\Psi$ [kohn65, ; sham66, ]. The main idea underlying this
theory, called density functional theory (DFT), is that the exact ground-state
energy, $E$, and electron density, $n({\bf r})$, can be determined by solving
an effective one-electron Schrödinger equation of the form:

$H_{\rm KS}\Psi_{i\sigma}=\epsilon_{i\sigma}\Psi_{i\sigma}~{},$ (2)

where $H_{\rm KS}$ is the Kohn-Sham Hamiltonian, index $i$ labels different
one-electron orbitals and $\sigma$ different spin states. The KS Hamiltonian
can be expressed as:

$H_{\rm KS}=-\frac{1}{2}\nabla^{2}+V_{ext}({\bf r})+\int\frac{n({\bf
r^{\prime}})}{|{\bf r}-{\bf r^{\prime}}|}d{\bf r^{\prime}}+V_{xc}({\bf
r})~{},$ (3)

where

$n({\bf r})=\sum_{i\sigma}|\Psi_{i\sigma}({\bf r})|^{2}~{},$ (4)

$V_{ext}$ represents an external field and $V_{xc}({\bf r})=\delta
E_{xc}/\delta n({\bf r})$ is the exchange-correlation potential.

The exchange-correlation energy has a purely quantum mechanical origin and can
be defined as the interaction energy difference between a quantum many-
electron system and its classical counterpart. Despite $E_{xc}$ represents a
relatively small fraction of the total energy, this contribution is extremely
crucial for all materials and molecules because it acts directly on the
bonding between atoms. In general, $E_{xc}[n]$ is unknown and needs to be
approximated. This is the only source of fundamental and uncontrollable errors
in DFT methods. In standard DFT approaches, $E_{xc}[n]$ typically is
approximated with the expression:

$E_{xc}^{\rm approx}[n]=\int\epsilon_{xc}^{\rm approx}({\bf r})n({\bf r})d{\bf
r}~{},$ (5)

where $\epsilon_{xc}^{\rm approx}$ is made to depend on $n({\bf r})$, $\nabla
n({\bf r})$, and/or the electronic kinetic energy $\tau({\bf
r})=\frac{1}{2}\sum_{i\sigma}|\nabla\Psi_{i\sigma}({\bf r})|^{2}$ . Next, we
summarise the basic aspects of the most popular $E_{xc}[n]$ functionals
employed in modeling of materials and phase transitions.

In local approaches (e.g., local density approximation –LDA–), $E_{xc}^{\rm
approx}$ in Eq.(5) is calculated by considering the exchange-correlation
energy of an uniform electron gas with density $n({\bf r})$,
$\epsilon_{xc}^{unif}$, which is known exactly from quantum Monte Carlo
calculations [ceperley80, ; perdew81, ]. In order to deal with the non-
uniformity of real electronic systems, the space is partitioned into
infinitesimal volume elements that are considered to be locally uniform. In
semi-local approaches (e.g., generalized gradient approximation –GGA–),
$E_{xc}$ is approximated similarly than in local approaches but
$\epsilon_{xc}^{\rm approx}$ is made to depend also on the gradient of $n({\bf
r})$ [perdew92, ; perdew96, ]. Both local and semi-local approximations
satisfy some exact $E_{xc}$ constraints and can work notably well for systems
in which the electronic density varies slowly over the space (e.g., crystals).
An extension of the GGA approach is provided by meta-GGA functionals, in which
the non-interacting kinetic energy density is considered also as an energy
functional input. An example of this latter type of functionals is the
recently proposed meta-GGA SCAN [scan, ; zhang18, ].

Hybrid functionals comprise a combination of non-local exact Hartree-Fock and
local exchange energies, together with semi-local correlation energies. The
proportion in which both non-local and local exchange densities are mixed
generally relies on empirical rules. The popular B3LYP approximation [becke93,
], for instance, takes a $20$% of the exact HF exchange energy and the rest
from the GGA and LDA functionals. Other well-known hybrid functionals are PBE0
[adamo99, ] and the range-separated HSE06 proposed by Scuseria and
collaborators [hse06, ]. In contrast to local and semi-local functionals,
hybrids describe to some extent the delocalisation of the exchange-correlation
hole around an electron hence they partially correct for electronic self-
interaction errors (which are ubiquitous in standard DFT) [franchini14, ].
This technical feature is specially useful when dealing with strongly
correlated systems that contain $d$ and $f$ electronic orbitals (e.g.,
transition-metal oxide perovskites) [bilc08, ; cazorla16b, ; cazorla17c, ].

###  II.2 Classical interatomic potentials

Using first-principles methods to describe the interactions between electrons
and ions in crystals requires dedicated computational resources. In some
cases, the interatomic interactions can be approximated satisfactorily by
analytical functions known as classical interatomic potentials or force fields
and consequently the simulations can be accelerated dramatically with respect
to _ab initio_ calculations. Classical interaction models contain a number of
parameters that are adjusted to reproduce experimental or _ab initio_ data,
and their analytical expressions are constructed based on physical knowledge
and intuition. The force matching method proposed by Ercolessi and Adams
[ercolessi94, ] is an example of a classical potential fitting technique that
is widely employed in the fields of condensed matter physics and materials
science [cazorla15c, ; mattoni15, ; hata17, ; cazorla18b, ]. Nonetheless, the
ways in which classical interatomic potentials are constructed are neither
straightforward nor uniquely defined and the thermodynamic intervals over
which they remain reliable are limited.

A typical pairwise interaction model that has been widely employed to simulate
materials at finite temperatures is the Coulomb–Buckingham (CB) potential,
which adopts the simple form [tinte99, ; jackson05, ; sepliarsky01, ;
hashimoto14, ]:

$V_{\alpha\beta}(r_{ij})=A_{\alpha\beta}e^{-\frac{r_{ij}}{\rho_{\alpha\beta}}}-\frac{C_{\alpha\beta}}{r_{ij}^{6}}+\frac{Z_{\alpha}Z_{\beta}}{r_{ij}}~{},$
(6)

where subscripts $\alpha$ and $\beta$ represent atomic species in the system,
$r_{ij}$ the radial distance between a pair of $\alpha$ and $\beta$ atoms
labelled $i$ and $j$ respectively, $Z$ ionic charges, and $A$, $\rho$ and $C$
are potential parameters. The CB potential is composed of three different
energy contributions. The exponential term accounts for short-range repulsive
forces resulting from the interactions between nearby electrons; the second
term represents long-range attractive interactions arising from dispersive van
der Waals forces; the third term is the usual Coulomb interaction between
point charges. In order to describe atomic polarizability effects, “core-
shell” modeling can be performed on top of the CB potential. In standard core-
shell approaches, each atom is decomposed into a charged core, which interact
with others through $V_{ij}$’s analogous to the expression shown in Eq.(6),
and a charged shell that is harmonically bound to the core [jackson05, ;
sepliarsky01, ; hashimoto14, ; gambuzzi14, ].

Some technologically relevant materials (e.g., oxide and hybrid organic-
inorganic perovskites) are characterized by a delicate balance between short-
range and long-range forces, which may be respectively originated by
electronic orbital hybridizations and the Coulomb interactions between
permanent electric dipoles and higher-order moments [cohen92, ]. The
simplicity of the pairwise interaction model in Eq.(6) may not be adequate to
fully grasp the complexity of such variety of interatomic interactions. Bond-
valence (BV) potentials, for instance, represent an improvement with respect
to the CB model because they can mimic chemical bonding in complex materials
more precisely [haomin17, ]. (Many other interatomic interaction models going
beyond CB exist, like the embedded atom potential [taioli07, ; cazorla07, ]
and multibody force fields [cazorla15c, ], but for the sake of focus those are
not covered here.)

A general BV potential is [grinberg02, ; grinberg04, ; shin05, ; liu13, ]:

$V_{\rm BV}(r,\theta)=V_{\rm bind}(r)+V_{\rm charge}(r)+V_{\rm rep}(r)+V_{\rm
nl}(\theta)~{},$ (7)

where the first term in the right-hand side represents the bond-valence
potential energy, the second the Coulomb potential energy, the third the
repulsive potential energy, and the fourth an angle potential energy (to
prevent unphysically large distortions of covalently connected structural
units; for instance, the oxygen octahedra in oxide perovskites and molecular
cations in hybrid organic-inorganic compounds).

The bond-valence energy term generally is expressed as:

$V_{\rm bind}(r)=\sum_{\alpha=1}^{N_{s}}S_{\alpha}\sum_{i=1}^{N_{\alpha}}\mid
V_{i\alpha}(r_{i})-V_{\alpha}\mid^{\gamma_{\alpha}},$ (8)

with

$V_{i\alpha}(r_{i})=\sum_{\beta=1}^{N_{s}}\sum_{j}^{NN}\left(\frac{r_{0}^{\alpha\beta}}{r_{ij}^{\alpha\beta}}\right)^{C_{\alpha\beta}},$
(9)

where $N_{s}$ represents the number of atomic species in the system,
$S_{\alpha}$ are fitting parameters, $N_{\alpha}$ the number of $\alpha$
atoms, $V_{\alpha}$ the desired atomic valence for ion $\alpha$,
$\gamma_{\alpha}$ fitting parameters typically set to $1$, $j$ an atomic index
that runs over nearest-neighbour (NN) ions, $r_{0}^{\alpha\beta}$ and
$C_{\alpha\beta}$ parameters determined by empirical rules, and
$r_{ij}^{\alpha\beta}$ the radial distance between ions $i$ and $j$.

For the repulsive energy term, $V_{\rm rep}$, the following expression
normally is employed:

$V_{\rm
rep}(r)=\epsilon\sum_{\alpha=1}^{N_{s}}\sum_{i=1}^{N_{\alpha}}\sum_{\beta=1}^{N_{s}}\sum_{j=1}^{N_{\beta}}\left(\frac{B_{\alpha\beta}}{r_{ij}^{\alpha\beta}}\right)^{12},$
(10)

where $\epsilon$ and $B_{\alpha\beta}$ are fitting parameters. Meanwhile, an
harmonic function is used for the angle potential energy that reads:

$V_{\rm nl}(\theta)=k\sum_{i=1}^{N_{\rm
bun}}\left(\theta_{i,x}^{2}+\theta_{i,y}^{2}+\theta_{i,z}^{2}\right),$ (11)

where $k$ is a fitting parameter, $N_{\rm bun}$ the number of covalently
bonded units in the simulated system (e.g., oxygen octahedral in oxide
perovskites and molecular cations in hybrid organic-inorganic compounds), and
$\\{\theta_{i,\gamma}\\}$ the angular degrees of freedom of such entities.

Reliable BV, or formally analogous, interatomic potentials have been developed
for archetypal ferroelectric materials (e.g., BaTiO3, PbTiO3, and
PbZr0.2Ti0.8O3) [grinberg02, ; grinberg04, ; shin05, ; liu13, ; xu15, ],
hybrid organic-inorganic perovskites (e.g., MAPbI3) [mattoni15, ; hata17, ]
and complex hydrides (e.g., Li2B12H12 and LiCB11H12) [sau21b, ], among others.
Recently, the outcomes of molecular dynamics simulations performed with some
of those classical force fields have predicted the existence of novel and
giant caloric effects [sau21, ; qi18, ; liu16, ].

###  II.3 Machine learning potentials

Situations in which the use of first-principles methods is prohibitive and the
available classical potentials are not versatile enough to reproduce the phase
transition phenomena of interest, machine learning techniques may be very
useful. Machine learning (ML) is a subfield of artificial intelligence that
exploits the systematic identification of correlation in data sets to make
predictions and analysis [behler10, ; rupp15, ].

When calculations are performed in a series of similar systems or a number of
configurations involving a same system, the results contain redundant
information. An example is to run a molecular dynamics simulation in which the
total internal energy and atomic forces are calculated at each time step;
after a sufficiently long time, points which are close in configurational
space and have similar energies are visited during the sampling of the
potential-energy surface. Such a redundancy can be exploited to perform
computationally intensive calculations (i.e., of first-principles type) only
in few selected configurations and then use machine learning (ML) to
interpolate between those, thus obtaining approximate solutions for most
configurations. The success of this approach depends on a balance between
incurred errors due to interpolation and invested computational effort.

ML modeling tools can provide both the energy and atomic forces directly from
the atomic positions, hence they can be regarded as a particular class of
atomistic potential. ML potentials, however, rely on very flexible analytic
functions rather than on physically motivated functionals hence they may be
highly versatile without requiring prior detailed knowledge of the simulated
system. Promising analytic approaches that have been recently proposed to
construct ML potentials include permutation invariant polynomials, the
modified Shepard method using Taylor expansions, Gaussian processes and
artificial neural networks. Artificial neural networks, in particular, have
been demonstrated to be a class of “universal approximators” [behler15, ]
since allow to approximate unknown multidimensional functions to within
arbitrary accuracy for a given set of known function values.

ML potentials are being employed with increasing frequency to solve problems
in materials science via atomistic simulations due to their ability to reach
high accuracy levels at moderate computational expense (see work [deringer19,
] for a general review on this topic). In the context of caloric effects, only
few works relying on ML techniques can be found in the literature (see, for
instance, works [zhang20, ; gong22, ]); however, based on the recent
progression realized in other similar research fields, we foresee a surge in
the application of ML methods to the analysis and prediction of caloric
effects in the near future.

##  III Free-energy methods

Barocaloric (BC) effects associated with first-order phase transitions that
entail abrupt volume and structural changes typically are large. The two main
physical descriptors of BC effects are the isothermal entropy change, $\Delta
S_{T}$, and adiabatic temperature change, $\Delta T_{S}$, associated to the
phase transition. These two quantities are thermodynamically related, as we
explain next. Consider that the entropy of a condensed matter system depends
on temperature and pressure, $S(T,P)$; an infinitesimal entropy variation then
can be expressed as:

$dS=\left(\frac{\partial S}{\partial T}\right)_{P}dT+\left(\frac{\partial
S}{\partial P}\right)_{T}dP~{}.$ (12)

In an adiabatic process the entropy is conserved, $dS=0$, thus from the
equation above it follows:

$\frac{C_{P}}{T}~{}dT=-\left(\frac{\partial S}{\partial P}\right)_{T}dP~{},$
(13)

in which the constant-pressure heat capacity $C_{P}\equiv T\left(\partial
S/\partial T\right)_{P}$ has been introduced.

By integrating both sides of the latter equality, one obtains:

$\Delta T_{S}=-\int\frac{T}{C_{P}}~{}dS_{T}~{}.$ (14)

which typically is approximated like [cazorla19b, ; moya14, ]:

$\Delta T_{S}\approx-\frac{T_{t}}{C_{P}}\cdot\Delta S_{T}~{},$ (15)

where $T_{t}$ corresponds to the phase transition temperature and $C_{P}$ to
the averaged heat capacity of the involved phases close to the transition
point.

Atomistic simulations generally are performed in the $(N,V,T)$ and $(N,P,T)$
ensembles hence the entropy and other related thermodynamic quantities can be
theoretically determined as a function of the state variables $V$, $P$, $T$
and $\rho$ ($\equiv N/V$). Concerning BC effects, once the thermodynamic
function $S(P,T)$ of the involved phases are known both the $\Delta S_{T}$ and
$\Delta T_{S}$ descriptors can be straightforwardly calculated [see
Eqs.(12)-(15) above]. Meanwhile, the entropy of a given system is related to
its Gibbs free energy, $G(P,T)$, through the expression:

$S=-\left(\frac{\partial G}{\partial T}\right)_{P}~{},$ (16)

thus it follows, without too much surprise, that the fundamental quantity
allowing for the theoretical determination of BC effects is $G(P,T)$. For this
reason, in this section we explain several computational approaches that allow
for direct estimation of free energies in crystals by using atomic interaction
modeling techniques like those introduced in the previous section (i.e., ab
initio, classical potential and machine learning methods).

###  III.1 The quasi-harmonic approximation

In the quasi-harmonic approach (QHA) [harmonic1, ; harmonic2, ; harmonic3, ]
one assumes that the potential energy of a crystal can be approximated with a
quadratic expansion around the equilibrium atomic configuration of the form:

$E_{\rm qh}=E_{\rm
eq}+\frac{1}{2}\sum_{l\kappa\alpha,l^{\prime}\kappa^{\prime}\alpha^{\prime}}\Phi_{l\kappa\alpha,l^{\prime}\kappa^{\prime}\alpha^{\prime}}u_{l\kappa\alpha}u_{l^{\prime}\kappa^{\prime}\alpha^{\prime}}~{},$
(17)

where $E_{\rm eq}$ is the total energy of the equilibrium lattice, $\bm{\Phi}$
the corresponding force-constant matrix, and $u_{l\kappa\alpha}$ the
displacement along Cartesian direction $\alpha$ of atom $\kappa$ at lattice
site $l$. The atomic displacements normally are expressed as:

$u_{l\kappa\alpha}(t)=\sum_{q}u_{q\kappa\alpha}\exp{\left[i\left(\omega
t-\bm{q}\cdot(\bm{l}+\bm{\tau}_{\kappa}\right)\right]}~{},$ (18)

where $\bm{q}$ is a wave vector in the first Brillouin zone (BZ) that is
defined by the equilibrium unit cell; $\bm{l}+\bm{\tau}_{\kappa}$ is the
vector that locates atom $\kappa$ at cell $l$ in the equilibrium structure.
The normal modes of the crystal then are found by diagonalizing the dynamical
matrix:

$\begin{split}&D_{\bm{q};\kappa\alpha,\kappa^{\prime}\alpha^{\prime}}=\\\
&\frac{1}{\sqrt{m_{\kappa}m_{\kappa^{\prime}}}}\sum_{l^{\prime}}\Phi_{0\kappa\alpha,l^{\prime}\kappa^{\prime}\alpha^{\prime}}\exp{\left[i\bm{q}\cdot(\bm{\tau}_{\kappa}-\bm{l^{\prime}}-\bm{\tau}_{\kappa^{\prime}})\right]}~{},\end{split}$
(19)

and thus the solid can be treated as a collection of non-interacting harmonic
oscillators with frequencies $\omega_{\bm{q}s}$ (positively defined and non-
zero) and energy levels:

$E_{\bm{q}s}=\hbar\omega_{\bm{q}s}\left(\frac{1}{2}+n\right)~{},$ (20)

where $n$ is the Bose-Einstein occupation number. Accordingly, the potential
energy of a crystal can be written within the QHA like:

$E_{\rm
qh}=\sum_{\bm{q}s}\hbar\omega_{\bm{q}s}\left[\frac{1}{2}+\frac{1}{\exp{\left(\hbar\omega_{\bm{q}s}/k_{B}T\right)}-1}\right]$
(21)

and the constant-volume heat capacity like:

$\displaystyle C_{V}(T)$ $\displaystyle=\left(\frac{\partial E_{\rm
qh}}{\partial T}\right)_{V}$ (22)
$\displaystyle=\sum_{\bm{q}s}\frac{\left(\hbar\omega_{\bm{q}s}\right)^{2}}{k_{B}T^{2}}\frac{\exp{\left(\hbar\omega_{\bm{q}s}/k_{B}T\right)}}{\left[\exp{\left(\hbar\omega_{\bm{q}s}/k_{B}T\right)}-1\right]^{2}}~{}.$

Likewise, the Helmholtz free energy, $F(V,T)$, is analytically expressed as:

$F(V,T)=E_{0}(V)+k_{B}T\sum_{\bm{q}s}\ln\left[2\sinh\left(\frac{\hbar\omega_{\bm{q}s}(V)}{2k_{\rm
B}T}\right)\right]~{},$ (23)

where $E_{0}$ corresponds to the energy of the static lattice and the
dependence of the vibrational frequencies on volume has been explicitly noted.
The Helmholtz and Gibbs free energies are related through the Legendre
transformation $G=F+PV$, in which the pressure is equal to $P=-\left(\partial
F/\partial V\right)_{T}$. Thus, by knowing $F(V,T)$, $C_{V}(T)$ and the
equation of state $P(V,T)$, one can deduce the value of $G$ and $C_{P}$ as a
function of $P$ and $T$.

Crystals in which anharmonic effects are important (e.g., materials in which
imaginary phonon frequencies appear in the vibrational phonon spectra
calculated at zero temperature), however, are not adequately described by the
QHA (e.g., the Helmholtz free energy expression in Eq.(23) is ill-defined and
cannot be evaluated). In such a case, one should resort instead to genuine
anharmonic free-energy approaches for the calculation of $F(V,T)$ like
molecular dynamics, thermodynamic integration from reference models [taioli07,
; cazorla07, ; cazorla12, ] and self-consistent phonon approaches [scp1, ;
scp2, ; scp3, ; scp4, ], to cite some examples. The computational load
associated with these anharmonic free-energy methods, however, typically is
orders of magnitude higher than that of quasi-harmonic approaches. Next, we
briefly describe some of them.

###  III.2 Molecular dynamics

Molecular dynamics (MD) is a computer simulation method in which the
trajectories of atoms in a many-body interacting system are determined
numerically by solving the Newton’s equations of motion [gulp, ; lammps, ].
For this purpose, one needs to know the forces acting between the particles
and their potential energies, which can be calculated by any of the methods
described in Sec.II. Thermostat and barostat techniques are employed in MD
simulations to render constant temperature and pressure conditions [frenkel-
book, ]. Temperature and anharmonic effects in materials are naturally
accounted for in MD simulations (quantum nuclear effects [harmonic1, ], on the
other hand, are systematically disregarded in this approach). Estimation of
quantities like $T$–renormalized phonon frequencies, $\omega(V,T)$, and
vibrational density of states, $g(\omega)$, can be obtained from MD
simulations by using phonon-mode decomposition [dynaphopy, ], force-constant
optimization [hiphive, ] and effective Hamiltonian fitting [hellman13, ]
techniques, to cite some examples.

When the interatomic interactions in MD simulations are described with first-
principles methods, typically DFT potentials, they are referred to as _ab
initio_ molecular dynamics (AIMD). AIMD simulations certainly are several
orders of magnitude computationally more intensive than MD simulations
performed with classical potentials. Nevertheless, thanks to the current
steady growth in computational power and improved algorithms design (e.g.,
linear scaling DFT methods [conquest, ]), reliable AIMD simulation of complex
materials currently is within reach (see, for instance, works [sagotra19, ;
miyazaki21, ; gebbia21, ; polek20, ; luo20, ]). Machine learning based
molecular dynamics, on the other hand, represents an excellent alternative to
AIMD simulations (provided that a reliable interaction model can be generated)
due to its great balance between numerical accuracy and computational expense
[jinnouchi19, ].

A standard approach to estimate $T$–renormalized $g(\omega)$’s from MD
simulations consists in calculating the Fourier transform of the velocity-
velocity autocorrelation function as follows [sagotra19, ; miyazaki21, ;
gebbia21, ]:

$g(\omega)=\frac{1}{N_{ion}}\sum_{i}^{N_{ion}}\int_{0}^{\infty}\langle{\bf
v}_{i}(t)\cdot{\bf v}_{i}(0)\rangle e^{i\omega t}dt~{},$ (24)

where ${\bf v}_{i}(t)$ represents the velocity of the atom labelled as $i$ at
time $t$, and $\langle\cdots\rangle$ denotes statistical average in the
$(N,V,T)$ ensemble. By knowing g($\omega$), one can easily estimate $F(V,T)$
from the integral expression that is equivalent to Eq.(23), namely:

$F(V,T)=E_{0}(V)+\int_{0}^{\infty}k_{B}T\ln\left[2\sinh\left(\frac{\hbar\omega}{2k_{\rm
B}T}\right)\right]g(\omega)d\omega~{}.$ (25)

Likewise, the constant-volume heat capacity can be computed by replacing the
summation in Eq.(22) with the corresponding integral expression, and the
entropy with the formula [phonopy, ]:

$\displaystyle S(V,T)$ $\displaystyle=-\left(\frac{\partial F}{\partial
T}\right)_{V}=\int_{0}^{\infty}\frac{\hbar\omega}{2T}\coth\left(\frac{\hbar\omega}{2k_{\rm
B}T}\right)-$ (26) $\displaystyle
k_{B}\ln\left[2\sinh\left(\frac{\hbar\omega}{2k_{\rm
B}T}\right)\right]g(\omega)d\omega~{}.$

Thus, just like in the QHA method, one can deduce the value of the
thermodynamic functions $G(P,T)$ and $C_{P}(T)$, and in turn of the
barocaloric descriptors $\Delta S_{T}$ and $\Delta T_{S}$, from the knowledge
of $F(V,T)$, $C_{V}(T)$ and $P(V,T)$.

###  III.3 Thermodynamic integration

More often than not, accurate estimation of the Helmholtz free energy within
an ample thermodynamic interval based on Eqs.(24)–(25) may turn out to be
computationally too intensive. In such a situation, thermodynamic integration
techniques, in their different variants, may result very useful since they
rely on the fact that thermodynamic quantities other than $F(V,T)$ (e.g., the
internal energy and pressure) can be numerically converged more efficiently in
practice. Two popular modalities of thermodynamic integration (TI) are (1)
standard TI and (2) TI from reference models [taioli07, ; cazorla07, ;
cazorla12, ; cazorla09, ].

Standard TI exploits two well-known thermodynamic relations. First, by
considering the hydrostatic pressure expression:

$\left(\frac{\partial F}{\partial V}\right)_{V}=-P~{},$ (27)

it is possible to estimate Helmholtz free energy shifts as a function of
volume at a fixed temperature like:

$F(V_{1},T)-F(V_{0},T)=\int_{V_{1}}^{V_{0}}P~{}dV~{}.$ (28)

Likewise, based on the following expression deduced from statistical
mechanics:

$\left(\frac{\partial\left[\beta F\right]}{\partial\beta}\right)_{V}=E~{},$
(29)

where $\beta\equiv 1/k_{B}T$, it is possible to estimate Helmholtz free energy
shifts as a function of temperature at a fixed volume like:

$\beta_{1}F(V,T_{1})-\beta_{0}F(V,T_{0})=\int_{\beta_{0}}^{\beta_{1}}E~{}d\beta~{}.$
(30)

In view of Eqs.(28)–(30), it is possible then to extend the value of the
Helmholtz free energy by performing MD simulations over a dense grid of
$(V,T)$ points, from which $P$ and $E$ can be determined, and subsequently
compute the value of the involved integrals numerically.

In the standard TI approach, however, only Helmholtz free energy differences
can be estimated hence for the computation of total free energies of different
phases, and therefore of barocaloric effects, the value of $F(V,T)$ should be
known at least at one thermodynamic state. TI from a reference model can be
very useful for this end. The general principle in this approach is that the
change of Helmholtz free energy is computed as the total energy function
$E_{\lambda}({\bf r}_{1},\ldots{\bf r}_{N})$ changes adiabatically from
$E_{0}$ to $E_{1}$, the free energies associated with these energy functions
being $F_{0}$ and $F_{1}$. Specifically, the formula for TI from a reference
model is [cazorla12, ; cazorla09, ]:

$F_{1}-F_{0}=\int_{0}^{1}\langle E_{1}-E_{0}\rangle_{\lambda}~{},$ (31)

where $\langle\cdots\rangle_{\lambda}$ represents the thermal average
evaluated for the system governed by the energy function
$E_{\lambda}=(1-\lambda)E_{0}+\lambda E_{1}$. In practice, one takes $E_{1}$
to be the total energy function whose free energy wish to calculate (e.g., a
specific DFT functional), and $E_{0}$ the total energy function of a
“reference” system whose free energy can be evaluated exactly. For solids, one
typically chooses the reference system to be a perfectly harmonic system, the
corresponding Helmholtz free energy expression being that in Eqs.(23)–(25).
Meanwhile, the integral in Eq. (31) should be computed numerically and for
this end one can perform a series of MD simulations in the $(N,V,T)$ ensemble
governed by the energy function $E_{\lambda}$ at different $\lambda$ values.

##  IV Estimation of BC effects

In analogy to the experiments, BC effects can be estimated with atomistic
simulations in an indirect, quasi-direct and direct fashion. In the indirect
approach, one exploits the Maxwell relation:

$-\left(\frac{\partial S}{\partial P}\right)_{T}=\frac{\partial^{2}G}{\partial
T\partial P}=\left(\frac{\partial V}{\partial T}\right)_{P}~{},$ (32)

to estimate the adiabatic temperature and isothermal entropy shifts like [see
Eq.(13)]:

$\displaystyle\Delta T_{S}=\int\frac{T}{C_{P}}~{}\left(\frac{\partial
V}{\partial T}\right)_{P}~{}dP$ $\displaystyle\Delta
S_{T}=-\int\left(\frac{\partial V}{\partial T}\right)_{P}~{}dP~{}.$ (33)

The temperature derivatives in Eq.(33), however, are ill-defined for the case
of first-order phase transitions since the volume of the system changes
discontinously with temperature during this type of transformation.
Accordingly, indirect approaches based on the Maxwell relation above are only
valid for estimating BC effects associated with second-order like phase
transitions in which the volume of the system changes continuously rather than
abruptly [cazorla16, ; sagotra17, ; sagotra18, ; min20, ].

As it was mentioned previously, large BC effects typically are associated with
first-order phase transitions. In such a case, one can approximately estimate
the size of $\Delta S_{T}$ and $\Delta T_{S}$ with the Clausius-Clapeyron
method or exactly with quasi-direct and direct techniques. We briefly describe
each of these approaches in the following sections.

###  IV.1 The Clausius-Clapeyron method

To a first approximation, the barocaloric isothermal entropy change associated
with a first-order $P$-induced phase transition can be estimated with the
Clausius-Clapeyron relation like [moya14, ]:

$\Delta S_{T}\approx\Delta S=\Delta V\cdot\frac{dP_{t}}{dT_{t}}~{},$ (34)

where $\Delta V$ is the volume change that the system experiences during the
phase transformation occurring at conditions $P_{t}(T_{t})$ and the derivative
is estimated on the corresponding phase boundary line. Meanwhile, the
corresponding adiabatic temperature change can be roughly determined with
Eq.(15).

The Clausius-Clapeyron (CC) method is not exact and entails few technical
problems in practice. First, in experiments the entropy change of the phase
transformation, $\Delta S$, is not necessarily equal to the BC descriptor
$\Delta S_{T}$ due to the presence of pervasive hysteresis and phase
coexistence effects [cazorla19b, ], which in general are not reproducible with
atomistic simulation techniques (see next Sec.IV.2). And second, the value of
the derivative and volume change entering Eq.(34) typically are accompanied by
large numerical uncertainties that are propagated into the evaluation of
$\Delta S$. Having said this, the CC method represents a fair approach to
estimate BC effects and it has been already employed in numerous experimental
and theoretical works [sau21, ; li19, ; bermudez17, ]. In order to estimate BC
effects with the highest accuracy, however, it is recommended to employ quasi-
direct and direct methods.

Figure 1: Sketch of the entropy curves $S(P,T)$ that can be theoretically
estimated with computational free-energy methods for a system undergoing a
pressure-induced first-order phase transition, $I\to II$. Hysteresis effects
have been disregarded. The BC descriptors $\Delta S_{T}$ and $\Delta T_{S}$
associated with the phase transition can be straightforwardly deduced from
these curves. The entropy changes associated with the phase transition,
$\Delta S$’s, are not necessarily equal to the BC isothermal entropy change
$\Delta S_{T}$.

###  IV.2 Quasi-direct estimation of BC effects

The theoretical quasi-direct method emulates the analogous experimental quasi-
direct approach (see, for instance, works [cazorla17b, ; lloveras19, ;
lloveras20, ]) in which essentially the pressure and temperature dependences
of the system entropy, $S(P,T)$, are determined (Fig.1). Upon derivation of
the full entropy curves (i.e., exhibiting a discontinuity, $\Delta S$, at the
first-order phase transition points) obtained at constant pressures $P_{0}$
and $P_{1}$ and expressed as a function of temperature, it is straightforward
to graphically deduce the size of the BC descriptors $\Delta S_{T}$ and
$\Delta T_{S}$ on the corresponding $S$–$T$ diagram, as it is illustrated in
Fig.1. From a practical point of view, the $S(P,T)$ curves can be calculated
numerically with any of the free-energy methods described in Sec.III, and the
phase transition points be determined under the condition
$G_{I}(P,T_{t})=G_{II}(P,T_{t})$ ($I$ and $II$ being the two phases involved
in the first-order transformation).

It is worth mentioning, however, that there are important differences between
the theoretical and experimental quasi-direct BC approaches. In experiments,
the full entropy curves, typically determined with calorimetric techniques,
inevitably include hysteresis and phase coexistence effects. Hysteresis, for
instance, originates from the fact that in practice transition paths occur
under nonequilibrium conditions; consequently, irreversible processes take
place within the material. Irreversibility may pose severe limitations to the
overall cooling performance attained during successive field-induced cycles,
hence it is an important aspect to take into consideration when envisaging
possible solid-state cooling applications [manosa17, ]. Unfortunately,
hysteresis, and in general any nonequilibrium process, cannot be reproduced
with the simulation techniques reviewed in Sec.III since perfect equilibrium
and reversibility conditions are always assumed in them. Thus, in most cases
the value of the theoretically estimated BC descriptors should be regarded as
upper bounds of the real $\Delta S_{T}$ and $\Delta T_{S}$ shifts (i.e., those
measured in the experiments).

###  IV.3 Direct estimation of BC effects

In this simulation approach, the system is first thermalized at the desired
initial temperature and pressure in the $(N,P,T)$, ensemble. After
thermalization, the simulation is switched to the isobaric-isoenthalpic
ensemble, $(N,P,H)$ [where $H$ represents the enthalpy of the system]. At this
stage the hydrostatic pressure is ramped up to the desired final value slowly
enough to guarantee adiabaticity, and the accompanying temperature change,
$\Delta T_{\rm on}$, is monitored. The system is simulated under these
conditions for some time. Subsequently, $P$ is ramped down to its initial
value slowly enough again to guarantee adiabaticity and the corresponding
temperature change, $\Delta T_{\rm off}$, and final temperature, $T_{f}$, are
monitored. Under the condition that the system remains in thermal equilibrium
during the entire described simulation cycle, it will follow that
$T_{i}=T_{f}$ and $\Delta T_{\rm on}=\Delta T_{\rm off}$ within the
corresponding statistical uncertainties. In such a case, either $\Delta T_{\rm
on}$ or $\Delta T_{\rm off}$ can be identified with the adiabatic temperature
change associated with the simulated BC process.

The main advantages of using the direct BC simulation method are the little
supervision required and possibility to combine pressure with additional
external fields (e.g., electric bias) in order to explore original
multicaloric effects [cazorla19b, ]. The disadvantages, on the other hand, are
large computational expense, no direct access to the isothermal entropy
change, and the fact that implementation of the $(N,P,H)$ ensemble is rare in
most molecular dynamics software packages. To the best of our knowledge,
direct estimation of BC effects has been only performed for the hybrid
organic-inorganic perovskite CH3NH3PbI3 [liu16, ], also known as MAPI.

##  V Representative examples

In this section, we analyze several illustrative cases in which original BC
effects have been simulated and quantified by using atomistic simulation
methods. The systems for which those BC effects have been predicted are
technologically relevant, namely, fast-ion conductors (LiN3), orientationally
disordered complex hydrides (Li2B12H12) and multiferroic perovskite oxides
(BiFe1-xCoxO3). In these cases, the involved materials modeling techniques
comprised molecular dynamics relying on classical force fields (LiN3 and
Li2B12H12) and the quasi-harmonic approximation performed with DFT methods
(BiFe1-xCoxO3). Analogous BC simulation success can be achieved also for other
families of pressure responsive materials.

Figure 2: Predicted barocaloric effects in bulk Li3N. (a) Isothermal entropy
and (b) adiabatic temperature shifts estimated for $\beta$-Li3N. (c)
Isothermal entropy and (d) adiabatic temperature shifts estimated for
$\alpha$-Li3N. Ball-stick representation of the $\beta$ and $\alpha$
polymorphs of Li3N; Li and N ions are represented with green and blue spheres,
respectively. Adapted from work [sagotra18, ].

###  V.1 Fast-ion conductor LiN3

Fast-ion conductors (FIC) are materials that exhibit high ionic conductivity
in the solid phase [hull04, ]. Examples of archetypal FIC, also known as
_superionic_ conductors, include alkali-earth metal fluorides (CaF2), oxides
(doped CeO2), and lithium-rich compounds (Li10GeP2S12). Inherent to
superionicity is a significant increase in the concentration of point defects
(e.g., Frenkel pair defects) associated to a particular sublattice of atoms in
the crystal. Due to their unique ion-transport properties, FIC are promising
materials for the realization of all-solid-state electrochemical batteries via
replacement of customary liquid electrolytes, which involves substantial
improvements with respect to standard batteries in terms of safety,
cyclability, and electrochemical performance [wang15, ]. Owing to the large
entropy change (typically of the order of $100$ J K-1 [hull04, ]) and external
tunability associated with the superionic phase transition [cazorla18b, ;
cazorla14, ; cazorla17d, ], FIC also are being established as auspicious
caloric materials for exploitation in solid-state cooling applications
[cazorla17b, ; cazorla16, ; sagotra18, ; min20, ; cazorla15, ].

Figure 2 shows the two common polymorphs of bulk Li3N, a structurally simple
and well-known lithium-based FIC. The $\alpha$ phase (hexagonal, space group
$P6/mmm$) has a layered structure composed of alternating planes of hexagonal
Li2N and pure Li+ ions. The $\beta$ phase (hexagonal, space group
$P6_{3}/mmmc$) presents an additional layer of lithium ions intercalated
between the Li2N planes that is accompanied by a doubling of the unit cell.
Exceptionally high ionic conductivities of the order of $10^{-4}$–$10^{-3}$ S
cm-1 have been experimentally observed in Li3N at room temperature li10 ;
alpen77 ; nazri94 . In molecular dynamics simulation studies based on
classical interaction potentials and first-principles methods [sagotra18, ],
it has been found that in order to reproduce room-temperature superionicity in
$\beta$–Li3N it is necessary to introduce a small concentration of extrinsic
Li+ vacancies in the system ($\sim 1$%); by contrast, stoichiometric
$\alpha$–Li3N already displays fast-ionic conductivity at $T=300$ K.

Figure 2 also encloses the theoretical estimation of BC effects in Li3N up to
$1$ GPa at room temperature and $T=400$ K [sagotra18, ]. Those calculations
relied on molecular dynamics simulations performed with classical interatomic
potentials (Sec.II.2) and the set of equations Eq.(33). In the $\beta$ phase,
a maximum room-temperature entropy shift of $-24$ JK-1kg-1 was obtained, which
slightly increases in absolute value at higher temperatures (Fig.2a). The
magnitude of such $\Delta S$ can be regarded as “giant” (i.e., $|\Delta S|>10$
JK-1kg-1). The accompanying adiabatic temperature shifts, however, are large
but not giant (that is, $|\Delta T|<10$ K, Fig.2b) owing to the huge heat
capacity of Li3N (e.g., $\sim 4\cdot 10^{3}$ JK-1kg-1 at $T=300$ K). Analogous
BC results were found in $\alpha$–Li3N (Figs.2c–d) although in this latter
case the estimated isothermal entropy and adabatic temperature changes were
slightly larger (for instance, $\Delta S=-32$ JK-1kg-1 and $\Delta T=+2.8$ K
at $T=300$ K and $\sigma=1$ GPa).

The large BC effects found in superionic Li3N are not driven by any structural
phase transformation. This assertion is confirmed, for instance, by the radial
distribution functions calculated for all pairs of ionic species in Li3N at
$T=300$ K since they present almost identical structural traits independently
of the applied pressure [sagotra18, ]. Such an absence of phase transformation
is in stark contrast to what is observed in other BC materials in which most
of the caloric response is concentrated near phase transition points [moya14,
]. Then, which is the principal mechanism behind the disclosed giant $\Delta
S$ in Li3N? It is well-known that hydrostatic pressure depletes significantly
ionic diffusivity, and therefore the entropy, in most fast-ion conductors
[cazorla18b, ; hull04, ; cazorla14, ; cazorla17d, ]. Essentially, the
available volume to interstitial ions is effectively reduced under compression
and as a consequence the kinetic barriers and formation energy of defects
governing ion migration drastically increase. Thus, increasing compression
steadily depletes ionic diffusion and in turn the entropy of the FIC crystal.

Therefore, in work [sagotra18, ] and original and rational strategy for
achieving large (and possibly also reversible) BC effects at ambient
temperature based on fast-ion conductors was introduced and illustrated for
bulk Li3N. Instead of focusing on the triggering of a structural phase
transition, which only serendipitously will occur at room temperature, the
starting point in this new scenario is a material that is already superionic
at ambient conditions. Large BC effects then can be obtained trough the
application of hydrostatic pressure since this affects significantly and
sustainedly the ionic conductivity, and in turn the entropy and volume, of
fast-ion conductors.

Figure 3: (a) Low-$T$ (ordered) and high-$T$ (disordered) phases of bulk
Li2B12H12. Li, B, and H ions are represented with red, blue, and yellow
colours, respectively. (b) Estimated Li2B12H12 phase diagram expressed as a
function of pressure and temperature. (c) Relative volume change associated
with the $T$-induced $\alpha\to\beta$ phase transition expressed as a function
of pressure (red); volume of the $\alpha$ phase per formula unit (f.u.),
$V_{\alpha}$, at the phase-transition conditions (black). (d) Isothermal
entropy (red) and adiabatic temperature (black) changes associated with the
barocaloric response of bulk Li2B12H12 expressed as a function of pressure.
Adapted from work [sau21, ].

###  V.2 Complex hydride Li2B12H12

Recently, colossal barocaloric effects (defined as $|\Delta S|\sim 100$
JK-1kg-1) have been measured in two different families of materials that
display intriguing order-disorder phase transitions [cazorla17b, ; li19, ;
lloveras19, ]. First, giant barocaloric effects have been theoretically
predicted [sagotra17, ] and experimentally observed in the archetypal
superionic compound AgI [cazorla17b, ]. AgI exhibits a first-order normal
(low-entropy) to superionic (high-entropy) phase transition that responds to
both temperature and pressure and which involves the presence of highly mobile
silver ions in the high–$T$ superionic state [hull04, ]. Likewise, the entropy
changes estimated for other similar superionic phase transitions in general
are also large [cazorla16, ; sagotra18, ; min20, ]. And second, colossal
barocaloric effects have been reported for the molecular solid neopentylglycol
[li19, ; lloveras19, ], (CH3)2C(CH2OH)2, and other plastic crystals
[lloveras20, ]. In these solids, molecules reorient almost freely around their
centers of mass, which remain localized at well-defined lattice positions.
Molecular rotations lead to orientational disorder, which renders high
entropy. By using hydrostatic pressure, it is possible to block such molecular
reorientational motion and thus induce a fully ordered state characterized by
low entropy [cazorla19a, ]. The barocaloric effects resulting from this class
of first-order order-disorder phase transition are huge and comparable in
magnitude to those achieved in conventional refrigerators with environmentally
harmful gases [li19, ; lloveras19, ; lloveras20, ].

In work [sau21, ], it has been predicted the occurrence of colossal
barocaloric effects ($|\Delta S|\sim 100$ JK-1kg-1) in the energy material
Li2B12H12 (LBH), a complex hydride that is already known from the fields of
hydrogen storage [her08, ; lai19, ; shevlin12, ] and solid-state batteries
[paskevicius13, ; mohtadi16, ]. By using molecular dynamics simulations
[sau19, ], a pressure-induced isothermal entropy change of $|\Delta S|=367$
JK-1kg-1 and an adiabatic temperature change of $|\Delta T|=43$ K at $T=480$ K
have been identified. These colossal entropy and temperature changes were
driven by small hydrostatic pressure shifts of $\sim 0.1$ GPa (Fig.3), thus
yielding gigantic barocaloric strengths of $|\Delta S|/P\sim 10^{3}$
JK-1kg-1GPa-1 and $|\Delta T|/P\sim 10^{2}$ K GPa-1. The colossal barocaloric
effects disclosed in bulk LBH are originated by simultaneous frustration of
Li+ diffusion and (BH)${}_{12}^{-2}$ icosahedra reorientational motion driven
by hydrostatic pressure.

Specifically, at ambient conditions lithium dodecahydrododecaborate
(Li2B12H12), LBH, presents an ordered cubic $Pa\overline{3}$ phase ($Z=4$),
referred to as $\alpha$, which is characterized by Li+ cations residing on
near-trigonal-planar sites surrounded by three (BH)${}_{12}^{-2}$ icosahedron
anions. In turn, each (BH)${}_{12}^{-2}$ anion resides in an octahedral cage
surrounded by six Li+ cations (Fig.3a). A symmetry preserving order-disorder
phase transition occurs at high temperatures ($\sim 600$ K) that stabilises a
disordered state, referred to as $\beta$, in which the Li+ cations are mobile
and the (BH)${}_{12}^{-2}$ anions present reorientational motion (Fig.3a). The
relative volume expansion that has been experimentally measured for such a
first-order order-disorder phase transition is $\Delta V^{\rm
expt}/V_{\alpha}^{\rm expt}\approx 8$% [paskevicius13, ] and the corresponding
phase transition enthalpy $\Delta H^{\rm expt}\approx 130$ kJ kg-1 [verdal14,
]. The huge volume variation and enthalpy associated with the
$\alpha\leftrightarrow\beta$ transformation could be propitious for
barocaloric purposes if the involved phase transition was responsive to small
hydrostatic pressure shifts of $\sim 0.1$ GPa.

Figure 3b shows the $P$–$T$ phase diagram of bulk LBH theoretically estimated
with MD techniques [sau19, ]. The coexistence line of the $\alpha$ and $\beta$
phases was determined by conducting numerous MD simulations at small $P$–$T$
shifts of $0.05$ GPa and $12.5$ K, and by monitoring the structural, Li+
diffusion, and (BH)${}_{12}^{-2}$ reorientational properties of the system. In
consistency with the experiments, a point in the $\alpha$–$\beta$ coexistence
line of Fig.3b was identified with sharp and simultaneous changes in the
volume (Fig.3c), Li+ diffusion coefficient ($D_{\rm Li}$), and
(BH)${}_{12}^{2-}$ reorientational frequency ($\lambda_{\rm B_{12}H_{12}}$) of
bulk LBH. It was found that the critical temperature of the
$\alpha\leftrightarrow\beta$ transformation can be certainly displaced by
hydrostatic pressure (i.e., roughly by $13$ K per $0.1$ GPa).

Consequently, the BC isothermal entropy and adiabatic temperature shifts
induced by pressures of $0\leq P\leq 0.4$ GPa in bulk LBH were theoretically
estimated in work [sau21, ]. For this end, and in view of the first-order
nature of the $\alpha\leftrightarrow\beta$ phase transformation paskevicius13
, the Clausius-Clapeyron method was employed (Sec.IV.1). The resulting $\Delta
S$ and $\Delta T$ values enclosed in Fig. 3d in fact render colossal
barocaloric effects. For example, at $T=480$ K and $P=0.1$ GPa an isothermal
entropy change of $-367$ JK-1kg-1 and an adiabatic temperature change of $+43$
K were estimated. The size of $\Delta S$ and $\Delta T$ were found to
gradually decrease under pressure (e.g., at $T=515$ K and $P=0.4$ GPa we
estimated $-229$ JK-1kg-1 and $+28$ K, respectively). Meanwhile, the predicted
LBH barocaloric effects were always direct, that is, $\Delta T>0$, since the
low-entropy ordered phase ($\alpha$) was stabilized over the high-entropy
disordered phase ($\beta$) under pressure ($\Delta S<0$).

The phase transition underlying the colossal barocaloric effects enclosed in
Fig.3 is quite remarkable since it combines key ingredients of fast-ion
conductors (i.e., ionic diffusion) and molecular crystals (i.e.,
reorientational motion), materials that individually have been proven to be
excellent BC materials. Thus, alkali-metal complex borohydrides
($A_{2}$B12H12, $A=$ Li, Na, K, Cs [udovic14, ; udovic20, ]) emerge as a
promising new family of caloric materials in which the salient phase-
transition features of fast-ion conductors and plastic crystals coexist and
cooperate to render colossal barocaloric effects, broadening so the range of
caloric materials with commercial potential for solid-state cooling
applications.

Figure 4: (a) Description of the two energetically competing structures in
bulk BFCO under pressure. (b) Gibbs free energy difference between the ${\cal
T}$ and ${\cal R}$ phases of bulk BFCO expressed as a function of pressure and
temperature. (c) Phase boundary delimiting the stability region of the ${\cal
T}$ and ${\cal R}$ phases of bulk BFCO. (d) Relative volume change occurring
during the pressure-induced ${\cal T}\to{\cal R}$ phase transition and
expressed as a function of temperature. (e) Isothermal entropy and (f)
adiabatic temperature changes estimated for the pressure-induced ${\cal
T}\to{\cal R}$ phase transition expressed as a function of temperature.

###  V.3 Supertetragonal oxide perovskite BiFe1-xCoxO3

Super-tetragonal (${\cal T}$) oxide perovskites comprise a family of materials
that are fundamentally intriguing and have great potential for ferroelectric,
piezoelectric, sensor, and energy conversion applications [zhang18b, ;
yamada13, ; infante11, ]. Super-tetragonal phases exhibit giant electric
polarizations of the order of $100$ $\mu$C/cm2 and may be accompanied by
magnetism. The coexistence of ferroelectricity and magnetism in crystals,
known as multiferroics, offers the possibility of controlling the
magnetization with electric fields via their cross-order coupling.
Magnetoelectric couplings can be used, for example, to design ultra efficient
logic and memory devices and realize large piezomagnetic coefficients for the
miniaturization of antennas and sensors. Furthermore, phase transitions
involving ${\cal T}$ phases typically exhibit colossal volume changes of $\sim
10$% (e.g., PbVO3 and related solid solutions), which can be exploited in
mechanical degradation [yamamoto19, ; pan19, ] applications. Examples of
${\cal T}$ multiferroic materials are bulk BiCoO3 (BCO) and BiFeO3 (BFO) thin
films [menendez20, ; belik06, ; wang03, ].

Recently, it has been shown by means of first-principles calculations based on
DFT that BiFe1-xCoxO3 solid solutions (BFCO) with $0.25\leq x\leq 0.50$
represent ideal bulk systems in which to realize the full potential of
multiferroic ${\cal T}$ materials [menendez20b, ]. Specifically, it has been
found that under moderate hydrostatic pressures of $0.1\lesssim P\lesssim 1$
GPa (depending on the relative Fe/Co content) it is possible to trigger a
phase transition from a low-$T$ rhombohedral (${\cal R}$) phase to a high-$T$
${\cal T}$ phase (Figs.4a–c) at room temperature. The disclosed pressure-
induced ${\cal R}\to{\cal T}$ phase transformation involves (i) a colossal
increase in the electric polarization of $\sim 150$%, (ii) the existence of a
robust net magnetization of $\approx 0.13$ $\mu_{B}$ per formula unit, and
(iii) a giant volume increase of $\Delta V\sim 10$% (Fig.4d). Examples of
technologies in which these multifunctional phenomena could have an impact
include pyroelectric energy harvesting [bowen14, ; hoffmann15, ] and solid-
state cooling [gottschall18, ].

Figures 4d and e show respectively the isothermal entropy and adiabatic
temperature changes induced by hydrostatic pressure in BiFe0.5Co0.5O3 solid
solutions. Those results have been obtained with the quasi-harmonic DFT method
(Sec.III.1) and the Clausius-Clapeyron indirect estimation approach (Sec.IV.1)
[menendez22, ]. First, it is noted that the ${\cal R}\to{\cal T}$ transition
temperature can be enormously shifted by small pressure changes (Fig.4c). For
instance, in BiFe0.5Co0.5O3 a compression decrease of only $0.06$ GPa moves
the transition temperature from room temperature down to about $200$ K.
Second, an isothermal entropy change of approximately $-1.5$ J K-1mol-1
(equivalent to $\approx-4.8$ J K-1kg-1) is estimated at room temperature along
with a $\Delta T$ of $\approx 5$ K. These figures, although cannot be regarded
as “giant”, are reasonably large and promising. And third, the size of these
two barocaloric descriptors significantly increase at lower temperatures. At
$T=200$ K, for instance, $\Delta S$ increases up to $-2.8$ J K-1mol-1 and
$\Delta T$ to $6.6$ K.

The driving pressures reported in Fig.4d are very large, of the order of $1$
GPa; however, in practice those values can be drastically reduced (i.e., by
one order of magnitude, $\sim 0.1$ GPa) by increasing the relative content of
Fe ions. For instance, by using analogous first-principles computational
methods a ${\cal R}\to{\cal T}$ phase transition pressure of $0.2$ GPa is
estimated at zero temperature for bulk BiFe0.75Co0.25O3 [menendez22, ]. This
theoretical outcome is in good agreement with the available experimental data
on the compositional phase diagram of BFCO solid solutions [hojo18, ]. Another
interesting aspect of BFCO solid solutions is that, besides hydrostatic
pressure, they can react to external magnetic and electric fields since these
are simultaneously magnetic and ferroelectric [menendez20b, ]. Such a
multifunctional quality converts BFCO solid solutions into a potentially ideal
playground in which to explore and take advantage of novel multicaloric
effects and cycles [cazorla19b, ; gottschall18, ]. Therefore, BFCO solid
solutions in particular and multiferroic materials exhibiting supertetragonal
phases in general, appear to be very promising materials for the development
of advanced solid-state cooling applications.

##  VI Conclusions

Atomistic computational techniques for the simulation of barocaloric effects
are already well established. The simulation methods surveyed in this Review
(e.g., quasi-harmonic and thermodynamic integration techniques) were
originally developed for the study of temperature- and pressure-induced phase
transformations in materials, and have already achieved great success in many
research disciplines (e.g., high-pressure physics). Nonetheless, the number of
atomistic simulation studies found in the literature on barocaloric effects
currently is quite reduced and mostly limited to indirect estimation
approaches. One of the main objectives of this Review is to improve this
situation by surveying key theoretical and computational concepts in the
estimation of barocaloric effects and to promote the use of alternative and
highly accurate direct and quasi-direct simulation methods.

As it has been illustrated here, atomistic simulation approaches can be used
to reproduce with reliability the barocaloric performance of complex materials
like multiferroics, in which the structural and electronic degrees of freedom
are strongly coupled, lithium-based hydrides, in which molecular ions can be
orientationally disordered, and fast-ion conductors, in which ionic diffusion
determines most of their physical properties. Analogous barocaloric simulation
success can be achieved for other families of functional materials exhibiting
complex and unconventional, as well as standard, responses to hydrostatic
pressure.

Atomistic simulation of barocaloric effects, therefore, can indeed help
enormously in developing new materials and strategies for boosting solid-state
cooling based on the application of hydrostatic pressure. Thus, the current
pressing challenge of finding new refrigeration technologies that are
environmentally friendly and also sustainable may greatly benefit from the
reliable and physically insightful computational methods surveyed in this
Review.

## References

  * (1) Cazorla, C., Novel mechanocaloric materials for solid-state cooling applications, Appl. Phys. Rev. 6, 041316 (2019).
  * (2) Moya, X., Kar-Narayan, S. $\&$ Mathur, N. D., Caloric materials near ferroic phase transitions, Nat. Mater. 13, 439 (2014).
  * (3) Aznar, A., Lloveras, P., Romanini, M., Barrio, M., Tamarit, J. Ll., Cazorla, C., Errandonea, D., Mathur, N. D., Planes, A., Moya, X. $\&$ Ma${\rm\tilde{n}}$osa, Ll., Giant barocaloric effects over a wide temperature range in superionic conductor AgI, Nat. Commun. 8, 1851 (2017).
  * (4) Cazorla, C $\&$ Errandonea, D., Giant mechanocaloric effects in fluorite-Structured superionic materials, Nano Lett. 16, 3124 (2016).
  * (5) Scott, J. F., Electrocaloric materials, Annu. Rev. Mater. Res. 41, 229 (2011).
  * (6) Pecharsky, V. K. $\&$ Gschneidner, Jr. K. A., Giant magnetocaloric effect in Gd5(Si2Ge2), Phys. Rev. Lett. 78, 4494 (1997).
  * (7) Bonnot, E., Romero, R., Ma${\rm\tilde{n}}$osa, Ll., Vives, E. $\&$ Planes, A., Elastocaloric effect associated with the martensitic transition in shape-memory alloys, Phys. Rev. Lett. 100, 125901 (2008).
  * (8) Shirsath, S. E. _et al._ , Interface-charge induced giant electrocaloric effect in lead free rerroelectric thin-film bilayers, Nano Lett. 20, 1262 (2020).
  * (9) Sagotra, A. K., Errandonea, D. $\&$ Cazorla, C., Mechanocaloric effects in superionic thin films from atomistic simulations, Nat. Commun. 8, 963 (2017).
  * (10) Sagotra, A. K., Chu, D. $\&$ Cazorla, C., Room-temperature mechanocaloric effects in lithium-based superionic materials, Nat. Commun. 9, 3337 (2018).
  * (11) Min, J., Sagotra, A. K. $\&$ Cazorla, C., Large barocaloric effects in thermoelectric superionic materials, Phys. Rev. Mater. 4, 015403 (2020).
  * (12) Sau, K., Ikeshoji, T., Takagi, S., Orimo, S.-I., Errandonea, D., Chu, D. $\&$ Cazorla, C., Colossal barocaloric effects in the complex hydride Li2B12H12, Sci. Rep. 11, 11915 (2021).
  * (13) Menéndez, C., Chu, D. $\&$ Cazorla, C., Oxygen-vacancy induced magnetic phase transitions in multiferroic thin films, npj Comput. Mater. 6, 76 (2020)
  * (14) Machado, P. _et al._ , Enhancement of phase stability and optoelectronic performance of BiFeO3 thin films via cation co-substitution, J. Mater. Chem. C 9, 330 (2021).
  * (15) Kohn, W. $\&$ Sham, L. J., Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * (16) Sham, L. J. $\&$ Kohn, W., One-particle properties of an inhomogeneous interacting electron gas, Phys. Rev. 145, 561 (1966).
  * (17) Ceperley, D. M. $\&$ Alder, B. J., Ground state of the electron gas by a stochastic method, Phys. Rev. Lett. 45, 566 (1980).
  * (18) Perdew, J. P. $\&$ Zunger, A., Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23, 5048 (1981).
  * (19) Perdew, J. P., J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh $\&$ Fiolhais, C., Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation, Phys. Rev. B 46, 6671 (1992).
  * (20) Perdew, J. P., K. Burke $\&$ Ernzerhof, M., Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * (21) Sun, J., Ruzsinszky, A. $\&$ Perdew, J. P., Strongly constrained and appropriately normed semilocal density functional, Phys. Rev. Lett. 115, 036402 (2015).
  * (22) Zhang, Y. _et al._ , Efficient first-principles prediction of solid stability: Towards chemical accuracy, npj Comp. Mater. 4, 9 (2018).
  * (23) Becke, A. D., Density‐functional thermochemistry. III. The role of exact exchange, J. Chem. Phys. 98, 5648 (1993).
  * (24) Adamo, C. $\&$ Barone, V., Toward reliable density functional methods without adjustable parameters: The PBE0 model, J. Chem. Phys. 110, 6158 (1999).
  * (25) Krukau, A. V., Vydrov, O. A., Izmaylov, A. F. $\&$ Scuseria, G. E., Influence of the exchange screening parameter on the performance of screened hybrid functionals, J. Chem. Phys. 125, 224106 (2006).
  * (26) Franchini, C., Hybrid functionals applied to perovskites, J. Phys. Condens. Matt. 26, 253202 (2014).
  * (27) Bilc, D. I., Orlando, R., Shaltaf, R., Rignanese, G.-M., ${\rm\acute{I}}$${\rm\tilde{n}}$iguez, J. $\&$ Ghosez, Ph., Hybrid exchange-correlation functional for accurate prediction of the electronic and structural properties of ferroelectric oxides, Phys. Rev. B 77, 165107 (2008).
  * (28) Rivero, P. $\&$ Cazorla, C., Revisiting the zero-temperature phase diagram of stoichiometric SrCoO3 with first-principles methods, Phys. Chem. Chem. Phys. 18, 30686 (2016).
  * (29) Cazorla, C., Lattice effects on the formation of oxygen vacancies in perovskite thin films, Phys. Rev. Appl. 7, 044025 (2017).
  * (30) Ercolessi, F. $\&$ Adams, J. B., Interatomic potentials from first-principles calculations: The force-matching method, Europhys. Lett. 26, 583 (1994).
  * (31) Cazorla, C. $\&$ Boronat, J., First-principles modeling of three-body interactions in highly compressed solid helium, Phys. Rev. B 92, 224113 (2015).
  * (32) Mattoni, A., Filipetti, A., Saba, M. I. $\&$ Delugas, P., Methylammonium rotational dynamics in lead halide perovskite by classical molecular dynamics: The role of temperature, J. Phys. Chem. C 119, 17421 (2015).
  * (33) Hata, T., Giorgi, G., Yamashita, K., Caddeo, C. $\&$ Mattoni, A., Development of a Classical Interatomic Potential for MAPbBr3, J. Phys. Chem. C 121, 3724 (2017).
  * (34) Cazorla, C., Sagotra, A. K., King, M. $\&$ Errandonea, D., High-pressure phase diagram and superionicity of alkaline earth metal difluorides, J. Phys. Chem. C 122, 1267 (2018).
  * (35) Tinte, S., Stachiotti, M. G., Sepliarsky, M., Migoni, R. L. $\&$ Rodriguez, C. O., Atomistic modelling of BaTiO3 based on first-principles calculations, J. Phys.:Condens. Matt. 11, 9679 (1999).
  * (36) Jackson, R. A. $\&$ Valerio, M. E. G., Structural phase transition and dynamical properties of PbTiO3 simulated by molecular dynamics, J. Phys.:Condens. Matt. 17, 837 (2005).
  * (37) Sepliarsky, M., Phillpot, S. R., Wolf, D., Stachiotti, M. G. $\&$ Migoni, R. L., Ferroelectric properties of KNbO3/KTaO3 superlattices by atomic-level simulation, J. Appl. Phys. 90, 4509 (2001).
  * (38) Hashimoto, T. $\&$ Moriwake, H., Dielectric properties of BaTiO3 by molecular dynamics simulations using a shell model, Mol. Simul. 41, 1074 (2015).
  * (39) Gambuzzi, E. $\&$ Pedone, A., On the structure of Ce-containing silicophosphate glasses: a core-shell molecular dynamics investigation, Phys. Chem. Chem. Phys. 16, 21645 (2014).
  * (40) Cohen, R. E., Origin of ferroelectricity in perovskite oxides, Nature 358, 136 (1992).
  * (41) Haomin, C. $\&$ Adams, S., Bond softness sensitive bond-valence parameters for crystal structure plausibility tests, IUCrJ 4, 614 (2017).
  * (42) Taioli, S., Cazorla, C., Gillan, M. J. $\&$ Alfè, D., Melting curve of tantalum from first principles, Phys. Rev. B 75, 214103 (2007).
  * (43) Cazorla, C., Gillan, M. J., Taioli, S. $\&$ Alfè, D., Ab initio melting curve of molybdenum by the phase coexistence method, J. Chem. Phys. 126, 194502 (2007).
  * (44) Grinberg, I., Cooper, V. R. $\&$ Rappe, A. M., Relationship between local structure and phase transitions of a disordered solid solution, Nature 419, 909 (2002).
  * (45) Grinberg, I., Cooper, V. R. $\&$ Rappe, A. M., Oxide chemistry and local structure of PbZrxTi1-xO3 studied by density-functional theory supercell calculations, Phys. Rev. B 69, 144118 (2004).
  * (46) Shin, Y.-H., Cooper, V. R., Grinberg, I. $\&$ Rappe, A. M., Development of a bond-valence molecular-dynamics model for complex oxides, Phys. Rev. B 71, 054104 (2005).
  * (47) Liu, S., Grinberg, I., Takenaka, H. $\&$ Rappe, A. M., Reinterpretation of the bond-valence model with bond-order formalism: An improved bond-valence-based interatomic potential for PbTiO3, Phys. Rev. B 88, 104102 (2013).
  * (48) Xu, R., Liu, S., Grinberg, I., Karthik, J., Damodaran, A. R., Rappe, A. M. $\&$ Martin, L. W., Ferroelectric polarization reversal via successive ferroelastic transitions, Nat. Mater. 14, 79 (2015).
  * (49) Sau, K., Ikeshoji, T., Kim, S., Takagi, S. $\&$ Orimo, S.-I. Comparative molecular dynamics study of the roles of anion–cation and cation–cation correlation in cation diffusion in Li2B12H12 and LiCB11H12, Chem. Mater. 33, 2357 (2021).
  * (50) Qi, Y., Liu, S., Lindenberg, A. M. $\&$ Rappe, A. M., Ultrafast electric field pulse control of giant temperature change in ferroelectrics, Phys. Rev. Lett. 120, 055901 (2015).
  * (51) Liu, S. $\&$ Cohen, R. E., Response of methylammonium lead iodide to external stimuli and caloric effects from molecular dynamics simulations, J. Phys. Chem. C 120, 17274 (2016).
  * (52) Behler, J., Neural network potential-energy surfaces for atomistic simulations, Chem. Modell. 7, 1 (2010).
  * (53) Rupp, M., Machine learning for quantum mechanics in a nutshell, Int. J. Quantum Chem. 115, 1058 (2015).
  * (54) Behler, J., Constructing high-dimensional neural network potentials: A tutorial review, Int. J. Quantum Chem. 115, 1032 (2015).
  * (55) Deringer, V. L., Caro, M. A. $\&$ Csanyi, G. Machine learning interatomic potentials as emerging tools for materials science, Adv. Mater. 31, 1902765 (2019).
  * (56) Zhang, Y. $\&$ Xu, X. Machine learning the magnetocaloric effect in manganites from compositions and structural parameters, AIP Adv. 10, 035220 (2020).
  * (57) Gong, J., Chu, S., Mehta, R. K. $\&$ McGaughey, A. J. H. XGBoost model for electrocaloric temperature change prediction in ceramics, npj Comput. Mater. 8, 140 (2022).
  * (58) Cazorla, C. $\&$ Boronat, J., Simulation and understanding of atomic and molecular quantum crystals, Rev. Mod. Phys. 89, 035003 (2017).
  * (59) Baroni, S., Giannozzi, P. $\&$ Isaev, E., Density-functional perturbation theory for quasi-harmonic calculations, Rev. Mineral. Geochem. 71, 39 (2010).
  * (60) Cazorla, C., Errandonea, D. $\&$ Sola, E., High-pressure phases, vibrational properties, and electronic structure of Ne(He)2 and Ar(He)2: A first-principles study, Phys. Rev. B 80, 064105 (2009).
  * (61) Cazorla, C., $\&$ Alfè, D. $\&$ Gillan, M. J., Constraints on the phase diagram of molybdenum from first-principles free-energy calculations, Phys. Rev. B 85, 064113 (2012).
  * (62) Tadano, T. $\&$ Tsuneyuki, S., Self-consistent phonon calculations of lattice dynamical properties in cubic SrTiO3 with first-principles anharmonic force constants, Phys. Rev. B 92, 054301 (2015).
  * (63) Souvatzis, P., Eriksson, O., Katsnelson, M. I. $\&$ Rudin, S. P., Entropy driven stabilization of energetically unstable crystal structures explained from first principles theory, Phys. Rev. Lett. 100, 095901 (2008).
  * (64) Souvatzis, P., Phonon lifetimes from first-principles self-consistent lattice dynamics, J.Phys.:Condens. Matter 23, 445401 (2011).
  * (65) Monacelli, L., Bianco, R., Cherubini, M., Calandra, M., Errea, I. $\&$ Mauri, F. The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects, J.Phys.:Condens. Matter 33, 363001 (2021).
  * (66) Gale, J. D., GULP: A computer program for the symmetry-adapted simulation of solids, J. Chem. Soc. Faraday Trans. 93, 629 (1997).
  * (67) Plimpton, S., Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comp. Phys. 117, 1 (1995).
  * (68) Frenkel, D. Understanding Molecular Simulation : from Algorithms to Applications, San Diego:Academic Press (1996).
  * (69) Carrera, A., Togo, A. $\&$ Tanaka, I. DynaPhoPy: A code for extracting phonon quasiparticles from molecular dynamics simulations, Comput. Phys. Commun. 221, 221 (2017).
  * (70) Eriksson, F., Fransson, E. $\&$ Erhart, P. The hiphive package for the extraction of high-order force constants by machine learning, Adv. Theory Simul. 2, 1800184 (2019)
  * (71) Hellman, O., Steneteg, P., Abrikosov, I. A. $\&$ Simak, S. I. Temperature dependent effective potential method for accurate free energy calculations of solids, Phys. Rev. B 87, 104111 (2013).
  * (72) Nakata, A. et al. Large scale and linear scaling DFT with the CONQUEST code, J. Chem. Phys. 152, 164112 (2020).
  * (73) Sagotra, A. K., Chu, D. $\&$ Cazorla, C. Influence of lattice dynamics on lithium-ion conductivity: A first-principles study, Phys. Rev. Mater. 3, 035405 (2019).
  * (74) Miyazaki, Y., Nakano, M., Krivchikov, A. I., Koroyuk, O. A., Gebbia, J. F., Cazorla, C. $\&$ Tamarit, J. Ll. Low-temperature heat capacity anomalies in ordered and disordered phases of normal and deuterated thiophene, J. Phys. Chem. Lett. 12, 2112 (2021).
  * (75) Szewczyk, D., Gebbia, J. F., Jezowski, A., Krivchikov, A. I., Guidi, T., Cazorla, C. $\&$ Tamarit, J. Ll. Heat capacity anomalies of the molecular crystal 1-fluoro-adamantane at low temperatures, Sci. Rep. 11, 18640 (2021).
  * (76) Polek, A., Cazorla, C. $\&$ Kundu, D. Nature of alkali ion conduction and reversible Na-ion storage in hybrid formate framework materials, J. Phys. Chem. C 124, 26714 (2020).
  * (77) Luo, X., Rawal, A., Cazorla, C. $\&$ Aguey-Zinsou, K.-F. Facile self-forming superionic conductors based on complex borohydrides surface oxidation, Adv. Sust. Sys. 4, 1900113 (2020).
  * (78) Jinnouchi, R., Lahnsteiner, J., Karsai, F., Kresse, G. $\&$ Bokdam, M. Phase transitions of hybrid perovskites simulated by machine-learning force fields trained on the fly with bayesian inference, Phys. Rev. Lett. 122, 225701 (2019).
  * (79) Togo, A. $\&$ Tanaka, I., First principles phonon calculations in materials science, Scr. Mater. 108, 1 (2015).
  * (80) Cazorla, C., Alfé, D. $\&$ Gillan, M. J. Melting properties of a simple tight-binding model of transition metals. I. The region of half-filled $d$-band, J. Chem. Phys. 130, 174707 (2009).
  * (81) Li, B. et al. Colossal barocaloric effects in plastic crystals, Nature 567, 506 (2019).
  * (82) Bermúdez-García, J. M., Sánchez-Andújar, M. $\&$ Se${\rm\tilde{n}}$arís-Rodríguez, M. A. A new playground for organic–inorganic hybrids: Barocaloric materials for pressure-induced solid-state cooling, J. Phys. Chem. Lett. 8, 4419 (2017).
  * (83) Lloveras, P. _et al._ , Colossal barocaloric effects near room temperature in plastic crystals of neopentylglycol, Nat. Commun. 10, 1803 (2019).
  * (84) Aznar, A. _et al._ , Reversible and irreversible colossal barocaloric effects in plastic crystals, J. Mater. Chem. A 8, 639 (2020).
  * (85) Ma${\rm\tilde{n}}$osa, Ll. $\&$ Planes, A. Materials with giant mechanocaloric effects: Cooling by strength, Adv. Mater. 29, 1603607 (2017).
  * (86) Hull, S. Superionics: crystal structures and conduction processes, Rep. Prog. Phys. 67, 1233 (2004).
  * (87) Wang, Y., Richards, W. D., Ong, S. P., Miara, L. J., Kim, J. C., Mo, Y. $\&$ Ceder, G. Design principles for solid-state lithium superionic conductors, Nat. Mater. 14, 1026 (2015).
  * (88) Cazorla, C. $\&$ Errandonea, D. Superionicity and polymorphism in calcium fluoride at high pressure, Phys. Rev. Lett. 113, 235902 (2014).
  * (89) Sagotra, A. K. $\&$ Cazorla, C. Stress-mediated enhancement of ionic conductivity in fast-ion conductors, ACS Appl. Mater. Interfaces 9, 38773 (2017).
  * (90) Cazorla, C. In the search of new electrocaloric materials: Fast ion conductors, Res. Phys. 5, 262 (2015).
  * (91) Li, W., Wu, G., Araújo, M., Scheicher, R. H., Blomqvist, A., Ahuja, R., Xiong, Z., Feng, Y. $\&$ Chen, P. Li+ ion conductivity and diffusion mechanism in $\alpha$-Li3N and $\beta$-Li3N, Energy Environ. Sci. 3, 1524 (2010).
  * (92) Alpen, U. V., Rabenau, A. $\&$ Talat, G. H. Ionic conductivity in Li3N single crystals, Appl. Phys. Lett. 30, 621 (1977).
  * (93) Nazri, G. A., Julien, C. $\&$ Mavi, H. S. Structure of Li3X (X = N, P, As) superionic conductors: X-ray diffraction and FTIR studies, Sol. Stat. Ionics 70, 137 (1994).
  * (94) Cazorla, C. Refrigeration based on plastic crystals, Nature 567, 470 (2019).
  * (95) Her, J.-H., Yousufuddin, M., Zhou, W., Jalisatgi, S. S., Kulleck, J. G., Zan, J. A., Hwang, S.-J., Bowman, R. C. $\&$ Udovic, T. J. Crystal structure of Li2B12H12: A possible intermediate species in the decomposition of LiBH4, Inorg. Chem. 47, 9757 (2008).
  * (96) Lai, Q., Sun, Y., Modi, P., Cazorla, C., Demirci, U. B., Ares, J. R., Leardini, F. $\&$ Aguey-Zinsou, K. F. How to design hydrogen storage materials? Fundamentals, synthesis, and storage tanks, Adv. Sustainable Syst. 3, 1900043 (2019).
  * (97) Shevlin, S. A., Cazorla, C. $\&$ Guo, Z. X. Structure and defect chemistry of low and high temperature phases of LiBH4, J. Phys. Chem. C 116, 13488 (2012).
  * (98) Paskevicius, M., Pitt, M. P., Brown, D. H., Sheppard, D. A., Chumphongphan, S. $\&$ Buckley, C. E. First-order phase transition in the Li2B12H12 system, Phys. Chem. Chem. Phys. 15, 15825 (2013).
  * (99) Mohtadi, R. $\&$ Orimo, S. I. The renaissance of hydrides as energy materials, Nat. Rev. Mater. 2, 16091 (2016).
  * (100) Sau, K., Ikeshoji, T., Kim, S., Takagi, S., Akgi, K. $\&$ Orimo, S.-I. Reorientational motion and Li+ ion transport in Li2B12H12 system: Molecular dynamics study, Phys. Rev. Mater. 3, 075402 (2019).
  * (101) Verdal, N., Her, J.-H., Stavila, V., Soloninin, A. V., Babanova, O. A., Skripov, A. V., Udovic, T. J. $\&$ Rush, J. J. Complex high-temperature phase transitions in Li2B12H12 and Na2B12H12, J. Sol. Stat. Chem. 212, 81 (2014).
  * (102) Udovic, T. J., Matsuo, M., Unemoto, A., Verdal, N., Stavila, V., Skripov, A. V., Rush, J. J., Takamura, H. $\&$ Orimo, S.-I. Sodium superionic conduction in Na2B12H12, Chem. Commun. 50, 3750 (2014).
  * (103) Jorgensen, M., Shea, P. T., Tomich, A. W., Varley, J. B., Bercx, M., Lovera, S., Cerny, R., Zhou, W., Udovic, T. J., Lavallo, V., Jensen, T. R., Wood, B. C. $\&$ Stavila, V. Understanding superionic conductivity in lithium and sodium salts of weakly coordinating closo-hexahalocarbaborate anions, Chem. Mater. 32, 1475 (2020).
  * (104) Zhang, L. et al., Giant polarization in super-tetragonal thin films through interphase strain, Science 361, 494 (2018).
  * (105) Yamada, H. et al., Giant electroresistance of super-tetragonal BiFeO3-based ferroelectric tunnel junctions, ACS Nano 7, 5385 (2013).
  * (106) Infante, I. C. et al., Multiferroic phase transition near room temperature in BiFeO3 films, Phys. Rev. Lett. 107, 237601 (2011).
  * (107) Yamamoto, H., Ogata, T., Saki, Y. $\&$ Azuma, M. Stability of polar structure in filling-controlled giant tetragonal perovskite oxide PbVO3, Inorg. Chem. 58, 2755 (2019).
  * (108) Pan, Z. et al., Large negative thermal expansion induced by synergistic effects of ferroelectrostriction and spin crossover in PbTiO3-based perovskites, Chem. Mater. 31, 1296 (2019).
  * (109) Belik, A. A. et al., Neutron powder diffraction study on the crystal and magnetic structures of BiCoO3, Chem. Mater. 18, 798 (2006).
  * (110) Wang, J. et al., Epitaxial BiFeO3 multiferroic thin film heterostructures, Science 299, 1719 (2003).
  * (111) Menéndez, C. $\&$ Cazorla, C., Giant thermal enhancement of the electric polarization in ferrimagnetic BiFe1-xCoxO3 solid solutions near room temperature, Phys. Rev. Lett. 125, 117601 (2020).
  * (112) Bowen, C. R., Taylor, J., LeBoulbar, E., Zabek, D., Chauhan, A. $\&$ Vaish, R. Pyroelectric materials and devices for energy harvesting applications, Energy Environ. Sci. 7, 3836 (2014).
  * (113) Hoffmann, M., Schroeder, U., Künneth, C., Kersch, A., Starschich, S., Böttger, U. $\&$ Mikolajick, T. Ferroelectric phase transitions in nanoscale HfO2 films enable giant pyroelectric energy conversion and highly efficient supercapacitors, Nano Energy 18, 154 (2015).
  * (114) Gottschall, T. et al., A multicaloric cooling cycle that exploits thermal hysteresis, Nat. Mater. 17, 929 (2018).
  * (115) Menéndez, C. $\&$ Cazorla, C., Colossal room-temperature electrocaloric strength aided by hydrostatic pressure in multiferroic BiFexCo1-xO3 solid solutions, To be published (2022).
  * (116) Hojo, H., Kengo, O., Keisuke, S., Yamamoto, H., Kawabe, R. $\&$ Azuma, M. Development of bismuth ferrite as a piezoelectric and multiferroic material by cobalt substitution, Adv. Mater. 30, 1705665 (2018).