# Beyond Gaussian fluctuations of quantum anharmonic nuclei.
The case of rotational degrees of freedom

Antonio Siciliano<EMAIL_ADDRESS>Dipartimento di Fisica,
Università di Roma La Sapienza, Piazzale Aldo Moro 5, 00185 Roma, Italy
Lorenzo Monacelli  Dipartimento di Fisica, Università di Roma La Sapienza,
Piazzale Aldo Moro 5, 00185 Roma, Italy  Theory and Simulation of Materials
(THEOS), and National Centre for Computational Design and Discovery of Novel
Materials (MARVEL), École Polytechnique Fédérale de Lausanne, 1015 Lausanne,
Switzerland  Francesco Mauri  Dipartimento di Fisica, Università di Roma La
Sapienza, Piazzale Aldo Moro 5, 00185 Roma, Italy

###### Abstract

The atomic motion in molecular crystals, such as high-pressure hydrogen or
hybrid organic-inorganic perovskites, is very complex due to quantum
anharmonic effects. In addition, these materials accommodate rotational
degrees of freedom. All the approximate methods that describe the nuclei
thermodynamic using Cartesian coordinates lead to an unphysical hybridization
of roto-librations with other high-energy modes. Hence, they do not accurately
account for the free energy contributions of these degrees of freedom. So, a
reliable description of a molecular crystal’s phase diagram is only possible
with Path Integral Molecular Dynamics (PIMD) at a high computational cost.
This work shows how to include roto-librational modes in the Self-Consistent
Harmonic Approximation (SCHA) framework. SCHA approximates the nuclei
Cartesian fluctuations to be Gaussian, thus neglecting curvilinear motion.
Keeping its low computational cost, we employ the generalization of SCHA,
called nonlinear SCHA (NLSCHA). Our method relies on a Gaussian ansatz for the
nuclei density matrix on a curved manifold, allowing us to map roto-librations
into harmonic modes defined on a surface. By optimizing the surface’s
curvature variationally, we minimize the free energy, allowing the spontaneous
activation of these degrees of freedom without external parameters. Notably,
in the limit of vanishing curvature, we recover the standard SCHA.

††preprint: APS/123-QED

##  I Introduction

Thanks to the recent methodological advantages [1, 2, 3, 4, 5] in the field of
computational condensed matter, the pivotal role of quantum fluctuations,
anharmonic effects, and finite temperature excitations on the equilibrium
ionic properties has been unveiled. We emphasize that a reliable free-energy
calculation should encompass all degrees of freedom. Indeed, the crystal
configuration can accommodate new types of atomic motion as it changes with
temperature and pressure.

In simplest structures, the only degrees of freedom are those we refer to as
’linear vibrations’, i.e. Gaussian fluctuations in a Cartesian space. In this
case, the lattice excitations are defined in a flat space, as breathing modes
or molecular stretching. Certain materials exhibit rotations and librations,
i.e. partial rotations, showcasing unique characteristics as atom movement is
confined to a curved surface. Indeed, the free rotation of a diatomic molecule
is the correlated motion of two atoms on a sphere, where the diameter
corresponds to the average bond length. The situation is more complex in the
case of a molecular crystal, where a group of atoms forms a rigid structure
strongly bonded. The latter displays low-energy modes as it can rotate freely
or partially without distortions. This type of motion has a low impact on the
internal energy but, on the contrary, makes a crucial contribution to the
total entropy, increasing the phase space available for the system. If we miss
these degrees of freedom, we may not detect phase transitions driven by
internal energy and entropy competition.

A prototypical example is found in the Rigid Unit Modes (RUMs) of framework
materials [6], formed by stiff connected polyhedrons of atoms, e.g. SiO4 and
AlO4 tetrahedra. RUMs have been identified as the soft modes responsible for
the displacive transitions, where the rigid units rotate and translate from
one phase to another [7]. Similar behavior is shown by the methyl group CH3 as
it can behave as a spinning top [8, 9] depending on the environment
surrounding it. This molecule is found in many pharmaceutical products, e.g.
paracetamol, and biological compounds, such as proteins and deoxyribonucleic
acid (DNA). Another interesting case is the high-pressure hydrogen phase
diagram. In phase I, the H2 centers form an h.c.p. lattice, and the quantum
distribution of protons is almost spherical, suggesting that the molecules
behave as free rotators [10]. The increasing pressure leads to a larger
intermolecular interaction, and phase II is stabilized at around $110$ GPa
[11]. Smaller molecular librations replace free rotations, hence the
orientational disorder is reduced. Rotations are also thermally activated as
in the III-IV phase transition where the molecules behave as free rotators
above $300$ K and $220$ GPa [12].

If used for a systematic and unbiased investigation of molecular crystal phase
diagrams, a reliable free energy method should accurately describe linear
vibrations and roto-librational modes. Only Molecular Dynamics (MD)
simulations can achieve this at the classical level. Indeed, MD offers the
advantage of including non-perturbative anharmonic effects while accurately
representing all the degrees of freedom, including rotations. Nevertheless, it
does not account for quantum effects, which Path-Integral Molecular Dynamics
(PIMD), the non-perturbative, reference-exact method, includes at finite
temperatures. However, PIMD carries a high computational cost, requiring the
simultaneous evolution of a number of the system’s replicas inversely
proportional to temperature. For this reason, the development of approximate
ionic free energy methods became of paramount importance.

The most commonly used and straightforward is the harmonic approximation (HA).
The HA fails when phonon interactions dominate due to substantial deviations
of atoms from equilibrium positions, and extensive regions of the BO energy
surface (BOES) are explored. Such a condition is relevant for light nuclei,
leading to sizeable zero-point motion (ZPM), near the melting point, or during
a second-order phase transition. In addition, the HA is not suited for
molecular crystals as it completely misses rotational modes. Indeed, in this
case, the small oscillations assumption of the HA breaks down, meaning that it
completely disregards the rotations of molecules or groups of atoms within a
molecular crystal.

If sufficiently small, anharmonic corrections can be incorporated via
perturbation theory starting from the HA. However, this approach becomes
cumbersome and is ineffective for hydride superconductors, ferroelectric, and
charge-density wave compounds. The growing interest in these materials
prompted the community to develop other methods, Self-Consistent Phonons (SCP)
theories [1, 2, 4], Vibrational Self-Consistent Field (VSCF) [5] and the
Temperature-Dependent Effective Potential (TDEP) method [3]. All of them, at
different levels of approximation, describe anharmonicity for linear
vibrations. However, none of these methods can effectively describe roto-
librations, as they all depend on Cartesian coordinates, which are
insufficient for capturing the atomic motion on a surface. So, the call for
approximate methods to handle rotations is more urgent than ever.

The self-consistent description of anharmonic phonons [1, 13, 14, 15, 16, 17,
18, 19] is a family of variational methods that constrain the nuclear
probability distribution to be Gaussian in Cartesian coordinates. This
approximation allows the enforcement of symmetries, which is an excellent help
for identifying phase transition, and interpolation techniques are available
to describe the thermodynamic limit, avoiding the simulation of big
supercells. Remarkably, none of these advantages are present in MD/PIMD.
Nevertheless, when atoms rotate freely or partially, the probability
distribution deviates significantly from a Gaussian shape [8, 12]. Therefore,
self-consistent phonon (SCP) theories [1, 13, 14, 15, 16] are not reliable in
such cases.

Here, we show how the Nonlinear Self-Consistent Harmonic Approximation
(NLSCHA) [20] can overcome this problem. We employ an ad-hoc change of
variables to deform the Gaussian ansatz and allow normal modes to occur in a
curved manifold, capturing the rotations of molecules and rigid body clusters.
The key aspect is that we variationally optimize the curvature, allowing the
system to spontaneously activate these modes only if the quantum free energy
is minimized. Remarkably, if the curvature is zero, the surface on which we
constrain the atomic motion becomes flat, and in this limit, we describe only
linear vibrations with the same accuracy as standard SCP theories.

In section II, we discuss the failure of SCP theories [1, 2, 4] in accounting
for roto-librational degrees of freedom with a 2D model for H2, as already
noted in Refs [8, 12]. To address such limitation, section III shows how to
incorporate these modes in the NLSCHA framework with negligible computational
cost. In the end, in sections IV-VI, we benchmark our method at zero and
finite temperature, and in section VII we show how to generalize our method to
the three-dimensional case of molecules and crystals.

##  II Failure on rotational modes

Here we show the failure of the SCP methods with molecular rotations. In
particular, among these methods, we consider the Self-Consistent Harmonic
Approximation (SCHA) [21, 1]. To this purpose, we solve the H2 molecule
rotating in two dimensions. In the center of mass reference frame, the only
degree of freedom is the relative coordinate,
$\bm{R}=\bm{R}_{1}-\bm{R}_{2}=(x,y)$, with an effective mass $m=m_{\ch{H}}/2$.
The BOES is given by a Morse potential fitted ab initio with DFT-BLYP on
$\ch{H2}$ (see appendix A) plus an empirical crystal field $E$ along the $x$
direction to control the rotational disorder

$V^{\text{(BO)}}(\bm{R})=V_{0}+d\left\\{1-\exp{-a(|\bm{R}|-R_{\text{eq}})}\right\\}+Ex$
(1)

where $|\bm{R}|=\sqrt{x^{2}+y^{2}}$

$\displaystyle V_{0}=$ $\displaystyle-1.172\text{ Ha}$ $\displaystyle d=$
$\displaystyle 0.137\text{ Ha}$ (2) $\displaystyle a=$ $\displaystyle
1.217\text{ Bohr}^{-1}$ $\displaystyle R_{\text{eq}}=$ $\displaystyle
1.393\text{ Bohr}$

Figure 1: Fig. a: exact, harmonic, and SCHA zero-point energy (ZPE), i.e. the
difference between the ground state energy and the potential minimum. Figs
b-g: exact and SCHA probability distributions (in Bohr-2) for three values of
the crystal field marked by horizontal lines in the upper panel
($E=0,0.005,0.6$ Ha/Bohr). SCHA can not describe roto-librations as its
Gaussian trial wavefunction unphysically hybridizes low-energy rotations with
high-energy vibrations.

In Fig. 1 a, we compare the exact, harmonic, and SCHA zero-point energy (ZPE)
as a function of the rotational freedom, $E$. In the harmonic approximation,
we expand $V^{\text{(BO)}}(\bm{R})$ around its minimum, and from the second-
order terms we extract the normal modes. The details of the exact
diagonalization and SCHA simulations are in appendix A and appendix B. The ZPE
is the difference between the ground state energy and the minimum of the
potential, representing the energy excess due to quantum uncertainty. In the
harmonic approximation, the ZPE is

$\text{ZPE}_{\text{harm}}=\frac{\hbar\left(\omega_{\text{harm,vib}}+\omega_{\text{harm,rot}}\right)}{2}$
(3)

where $\omega_{\text{harm,vib}},\omega_{\text{harm,rot}}$ are the frequencies
of the vibrational and rotational modes (respectively polarized along x and
y). In Fig. 1 b, d, f, we represent the probability distribution of the
relative coordinate $\bm{R}$ increasing the value of the crystal fields. In
the absence of $E$, the H2 molecule behaves like a free rotator (Fig. 1 b),
then rotations are progressively suppressed as $E$ increases (Fig. 1 d) until
the molecule is locked with a fixed orientation (Fig. 1 f).

In the harmonic approximation, $\omega_{\text{harm,rot}}$ is zero at $E=0$
Ha/Bohr, indicating the presence of a free rotator mode. So, there is a
direction along which the propagation costs zero energy, thereby undermining
the assumption of small oscillations on which the HA relies. Consequently, in
the case of full rotational invariance, the HA is not justified at finite
temperatures. In Fig. 2, we compare the exact, SCHA and harmonic free energies
at $1000$ K for low values of $E$. In the limit of vanishing crystal field,
the harmonic entropy diverges logarithmically as the rotational frequency
tends to zero. Consequently, the HA is unreliable at finite temperatures for
low values of $E$.

Figure 2: Exact, SCHA and harmonic free energies at $1000$ K for low values of
the crystal field $E$. As we reduce $E$, the rotational frequency is vanishing
$\omega_{\text{harm,rot}}\rightarrow 0$, hence the harmonic free energy
diverges as $f_{0}\log(E/E_{0})$ (solid dotted line).

In contrast, the SCHA method effectively handles large nuclear fluctuations by
adjusting the width and center of the Gaussian, minimizing the total energy.
Thus, we always get well-defined solutions. However, the SCHA method samples
the BOES using a Gaussian distribution in Cartesian coordinates, Fig. 1 c, e,
g. The approximation employed by SCHA results in an overestimation of the
rotational barrier, as it hybridizes vibrational and rotational modes, causing
an unphysical stiffening of the latter. This is evident when comparing SCHA
and exact probability distribution, Fig. 1 b-e. So a Gaussian trial density
matrix lacks the flexibility to describe free or partially rotating molecules,
thus sampling high-energy regions of the BOES. As a result, the SCHA ZPEs are
unreliable at low crystal fields. They can be worse than the predictions from
the HA, that, however provides a non-variational energy. Therefore, the SCHA
method fails to meet the gold-standard free energy method requirements when
rotations are present [8, 12]. Nevertheless, as $E$ is increased, rotations
are suppressed, and the SCHA method performs better, adeptly capturing the
anharmonic effects of the Morse potential and reducing its error in comparison
to the exact solution, as seen in Fig. 1 f-g.

Here, we demonstrated again that the SCHA is a valuable tool when linear
vibrations dominate as already noted in Refs [8, 12]. However, as soon as
rotational modes become active, the method falters due to the trial density
matrix’s lack of flexibility. All the other methods belonging to the SCP
family have the same problem. So, SCP methods can not detect the presence of a
rotational mode, leading to an overestimation of the total energy. To address
these issues, we employ the NLSCHA theory [20], which adapts a new trial
density matrix by automatically activating a roto-librational mode if it
lowers the total energy.

##  III The solution of nonlinear Self-Consistent Harmonic Approximation

The SCP approach falls short in considering roto-librations, as these modes
generate non-Gaussian fluctuations (refer to Fig. 1 b-d), surpassing the
capabilities of these methods. We employ the NLSCHA theory to modify the
Gaussian ansatz with an invertible nonlinear transformation so that it can
accommodate rotations.

###  III.1 The nonlinear change of variables

We aim to disentangle rotations from stretching modes as the SCHA leads to
unphysical hybridization. The solution of NLSCHA is to variationally optimize
a harmonic density matrix in a suitable auxiliary space that separates these
modes.

Figure 3: The stereographic transformation represents the Earth’s surface
(Fig. a) on a 2D plane (Fig. b). Each point on the Earth, identified by
$(\theta,\phi)$, is mapped to a unique point $(x_{\text{s}},y_{\text{s}})$ on
a plane tangent to the North Pole (N) (Fig. c). To build this transformation,
We connect S to a point $(\theta,\phi)$ then the stereographic projection
$(x_{\text{s}},y_{\text{s}})$ is given by the intersection of this line with
the tangent plane as Fig. c shows. This change of variables applied to the
atomic coordinates disentangles rotational modes (e.g. the motion on the
stereographic plane) from the vibrational ones (e.g. the motion perpendicular
to the surface). The plot was made with the cartopy package [22].

The most straightforward choice to describe rotations is spherical coordinates
$(r,\theta,\phi)$ [12]

$\displaystyle r$ $\displaystyle=\sqrt{x^{2}+y^{2}+z^{2}}$ (4a)
$\displaystyle\phi$ $\displaystyle=\atan\left(\frac{y}{x}\right)$ (4b)
$\displaystyle\theta$
$\displaystyle=\arccos\left(\frac{z}{\sqrt{x^{2}+y^{2}+z^{2}}}\right)$ (4c)

Note that it is possible to define small oscillations (i.e. the harmonic
approximation) only in 2D $(r,\phi)$ but not in 3D $(r,\phi,\theta)$ as small
variations of the polar angle $\theta$ around $\theta\simeq 0$ lead to large
changes of the azimuthal one $\phi$. In addition, Eqs (4) have a different
topology from the Cartesian space, indeed $(r,\phi,\theta)$ are not defined in
$\mathbb{R}^{3}$. To use the NLSCHA framework, we must work with an auxiliary
space with the same topology as the Cartesian one so that small oscillations
are well-defined and the entropy is an analytical quantity [20]. It is not a
problem working with $r$ as it represents the atomic distance and exploring
$r\rightarrow 0$ has a low probability. Indeed, $r\rightarrow 0$ has a
substantial energetic cost as it would correspond to nuclear fusion. So, we
have to solve the problem of the angular variables’ bounded range. If we adopt
the stereographic projection, the topology of $(\phi,\theta)$ becomes the same
as the Cartesian space. This technique represents the Earth on a 2D plane with
range $\mathbb{R}^{2}$. All the Earth’s surface points are projected on a
tangent plane at the North Pole (N) using the South Pole (S) as the projecting
point. Note that we can not represent S in this way and that no invertible
transformation can achieve this.

The NLSCHA auxiliary manifold we define here derives from a parameterization
of the stereographic transformation as it allows to explore more phase space,
i.e. it represents both North and South Hemispheres in cartography. For our
purposes, this property implies that we describe large angular fluctuations.
To better illustrate our approach, we consider the H2 molecule in 2D. We
define the auxiliary variables $\bm{u}=(u_{1},u_{2})$ from

$\bm{R}=\bm{\xi}(\bm{u})$ (5)

where $\bm{R}=(x,y)$ is the H2 relative coordinate in the center of mass
reference frame and $\bm{\xi}(\bm{u})$ is the stereographic projection

$\displaystyle x(\bm{u})$
$\displaystyle=x_{\text{C}}+(u_{1}+r_{0})\cos(\phi(u_{2}))$ (6a)
$\displaystyle y(\bm{u})$
$\displaystyle=y_{\text{C}}+(u_{1}+r_{0})\sin(\phi(u_{2}))$ (6b)
$\displaystyle\phi(u_{2})$
$\displaystyle=\phi_{0}+2\atan\left(\frac{u_{2}}{2r_{0}}\right)$ (6c)

Note that $\bm{\xi}(\bm{u})$ is parametrized by $x_{\text{C}}$,
$y_{\text{C}}$, $r_{0}$, and $\phi_{0}$. They represent the center of
curvature, $\bm{\mathcal{R}}_{\text{C}}$, and the curvature vector,
$\bm{\mathcal{R}}_{\text{T}}$,

$\displaystyle\bm{\mathcal{R}}_{\text{C}}$
$\displaystyle=\begin{bmatrix}x_{\text{C}}&y_{\text{C}}\end{bmatrix}$ (7a)
$\displaystyle\bm{\mathcal{R}}_{\text{T}}$
$\displaystyle=\begin{bmatrix}r_{0}\cos(\phi_{0})&r_{0}\sin(\phi_{0})\end{bmatrix}$
(7b)

where $\bm{\mathcal{R}}_{\text{T}}$ contains the information on the curvature
$\kappa$ of the nonlinear transformation

$\kappa=|\bm{\mathcal{R}}_{\text{T}}|^{-1}=r_{0}^{-1}$ (8)

Figure 4: A geometrical representation of the nonlinear transformation, Eq.
(6), with curvature $\kappa>0$. $\bm{R}=\bm{\xi}(\bm{u})$ corresponds to a
stereographic projection in 2D. The center and radius of the circle are
defined respectively by $\bm{\mathcal{R}}_{\text{C}}$ (red arrow) and
$|\bm{\mathcal{R}}_{\text{T}}|$ (green arrow).
$\bm{\mathcal{R}}_{\text{C}}+\bm{\mathcal{R}}_{\text{T}}$ indicates the
position of the tangent line, i.e. the North Pole. The stereographic
projection maps the points $u_{2}$ on angles $\phi$, with the South Pole
$\bm{\mathcal{R}}_{\text{C}}-\bm{\mathcal{R}}_{\text{T}}$ being the projecting
point.

We examine Fig. 4 to illustrate the geometrical meaning of Eq. (6) with
$\kappa>0$. The vector $\bm{\mathcal{R}}_{\text{C}}$ identifies the circle’s
center (center of the Earth) with radius $|\bm{\mathcal{R}}_{\text{T}}|$.
$u_{1}=|\bm{R}-\bm{\mathcal{R}}_{\text{C}}|-r_{0}$, the radial coordinate, and
$u_{2}$, the stereographic projection of $\phi-\phi_{0}$, identify the
position of a point $\bm{R}=(x,y)$ in the new space. We perform the
stereographic projection on the tangent plane at
$\bm{\mathcal{R}}_{\text{C}}+\bm{\mathcal{R}}_{\text{T}}$ (the North Pole)
using $\bm{\mathcal{R}}_{\text{C}}-\bm{\mathcal{R}}_{\text{T}}$ (the South
Pole) as the projecting point. Note that the North and South Pole are free
parameters of the transformation. As
$|\bm{\mathcal{R}}_{\text{T}}|\rightarrow+\infty$, the curvature
$\kappa\rightarrow 0$ so we recover a linear transformation (see section
III.2)

$\bm{R}\simeq\bm{\mathcal{R}}_{\text{C}}-\bm{\mathcal{R}}_{\text{T}}+\bm{u}=\bm{\mathcal{R}}+\bm{u}$
(9)

which is the one employed by standard SCHA [1]
($\bm{\mathcal{R}}_{\text{C}}-\bm{\mathcal{R}}_{\text{T}}=\bm{\mathcal{R}}$ is
the average atomic position in SCHA).

###  III.2 The trial density matrix

Here, we show how to incorporate the ad-hoc change of variables, Eq. (6), with
NLSCHA [20]. Within our framework, we variationally minimize the quantum free
energy with a trial density matrix corresponding to a Gaussian probability
distribution in $\bm{u}$, not in $\bm{R}$ as in the SCHA. Indeed, depending on
$\kappa$, we have a probability distribution that supports roto-librations.
The NLSCHA density operator is defined by the matrix elements (see appendix C)

$\bra{\bm{R}}\hat{\rho}_{\text{nl}}\ket{\bm{R}^{\prime}}=\frac{\overline{\rho}_{\text{nl}}(\bm{u},\bm{u}^{\prime})}{\sqrt{\mathcal{J}(\bm{u})\mathcal{J}(\bm{u}^{\prime})}}$
(10)

where $\mathcal{J}$ is the Jacobian’s determinant

$\mathcal{J}(\bm{u})=\det\left(\partialderivative{\bm{\xi}(\bm{u)}}{\bm{u}}\right)>0$
(11)

and $\overline{\rho}_{\text{nl}}(\bm{u},\bm{u}^{\prime})$ satisfies
($\beta^{-1}=k_{\text{B}}T$)

$\displaystyle\overline{\rho}_{\text{nl}}(\bm{u},\bm{u}^{\prime})$
$\displaystyle=\frac{\bra{\bm{u}}\exp{-\beta\hat{\overline{\mathcal{H}}}_{\text{nl}}}\ket{\bm{u}^{\prime}}}{\overline{\mathcal{Z}}_{\text{nl}}}$
(12a) $\displaystyle\overline{\mathcal{Z}}_{\text{nl}}$
$\displaystyle=\int_{-\infty}^{+\infty}d\bm{u}\bra{\bm{u}}\exp{-\beta\hat{\overline{\mathcal{H}}}_{\text{nl}}}\ket{\bm{u}}$
(12b)

In NLSCHA the auxiliary harmonic Hamiltonian
$\hat{\overline{\mathcal{H}}}_{\text{nl}}$ is defined in $\bm{u}$-space

$\bra{\bm{u}}\hat{\overline{\mathcal{H}}}_{\text{nl}}\ket{\bm{u}^{\prime}}=\delta(\bm{u}-\bm{u}^{\prime})\left(-\frac{\hbar^{2}}{2}\partialderivative{\bm{u}}\cdot\overset{-1}{\bm{\mathcal{M}}}\cdot\partialderivative{\bm{u}}+\frac{1}{2}\bm{u}\cdot\bm{\Phi}_{\text{nl}}\cdot\bm{u}\right)$
(13)

The variational parameters of $\hat{\rho}_{\text{nl}}$ (Eq. (10)) are the
force constant $\bm{\Phi}_{\text{nl}}$, the mass tensor $\bm{\mathcal{M}}$ and
the free parameters of the nonlinear transformation $\bm{\xi}$

$\bm{\Gamma}_{\text{nl}}=(\bm{\mathcal{R}}_{\text{C}},\bm{\mathcal{R}}_{\text{T}})$
(14)

Note that both $\bm{\Phi}_{\text{nl}}$ and $\bm{\mathcal{M}}$ are symmetric
and positive definite.

Figure 5:  Figs a-b show respectively $\overline{\rho}_{\text{nl}}(\bm{u})$
and the corresponding physical probability distribution
$\bra{\bm{R}}\hat{\rho}_{\text{nl}}\ket{\bm{R}}$ with finite curvature
$\kappa=|\bm{\mathcal{R}}_{\text{T}}|^{-1}>0$. Figs c-d report the same
quantities but with zero curvature $\kappa\rightarrow 0$. When $\kappa$ is
finite, $\bra{\bm{R}}\hat{\rho}_{\text{nl}}\ket{\bm{R}}$ describes roto-
vibrational modes, see Fig. b. Otherwise, if $\kappa\rightarrow 0$, the metric
of the nonlinear transformation is trivial so
$\bra{\bm{R}}\hat{\rho}_{\text{nl}}\ket{\bm{R}}$ becomes a SCHA-like
distribution, see Fig. d.

We prove that $\bra{\bm{R}}\hat{\rho}_{\text{nl}}\ket{\bm{R}}$ describes roto-
librational modes. In Fig. 5 we compare $\overline{\rho}_{\text{nl}}(\bm{u})$,
Eqs (12), with the corresponding physical probability distribution
$\bra{\bm{R}}\hat{\rho}_{\text{nl}}\ket{\bm{R}}$, Eq. (10), showing the effect
of the nonlinear transformation when $\kappa>0$ (Fig. 5 a, b) or
$\kappa\rightarrow 0$ (Fig. 5 c, d). Note that in Fig. 5, we keep
$\bm{\Phi}_{\text{nl}}$ and $\bm{\mathcal{M}}$ fixed. In the auxiliary
$\bm{u}$-space $\overline{\rho}_{\text{nl}}(\bm{u})$ is a Gaussian, see Fig. 5
a-c. When the nonlinear change of variables is applied, the shape of
$\bra{\bm{R}}\hat{\rho}_{\text{nl}}\ket{\bm{R}}$ changes depending on the
curvature $\kappa$, see Fig. 5 b, d. As expected from section III.1, a finite
$\kappa$ bends the probability distribution thanks to the stereographic
projection. So $\bra{\bm{R}}\hat{\rho}_{\text{nl}}\ket{\bm{R}}$ describes a
molecule that simultaneously rotates and vibrates, see Fig. 5 b. On the
contrary, when $\kappa\rightarrow 0$, we obtain a Gaussian probability
distribution, recovering the standard SCHA, see Fig. 5 d.

From Fig. 5 we deduce the newly introduced variational manifold is a superset
of the SCHA, ensuring that NLSCHA systematically outperforms the SCHA for
roto-librations. We variationally estimate the exact BO free energy (see
appendix D)

$\displaystyle F^{\text{(BO)}}$ $\displaystyle\leq F_{\text{nl}}$ (15a)
$\displaystyle F_{\text{nl}}$
$\displaystyle=\Tr\left[\hat{\rho}_{\text{nl}}\hat{H}^{\text{(BO)}}\right]+k_{\text{B}}T\Tr\left[\hat{\rho}_{\text{nl}}\log(\hat{\rho}_{\text{nl}})\right]$
(15b)

where the entropic term
$-k_{\text{B}}\Tr\left[\hat{\rho}_{\text{nl}}\log(\hat{\rho}_{\text{nl}})\right]$
has a harmonic form and coincides with the temperature derivative of
$F_{\text{nl}}$ only if we optimize all the free parameters [20]. We minimize
$F_{\text{nl}}$ with respect to the free parameters

$\displaystyle\partialderivative{F_{\text{nl}}}{\bm{\Phi}_{\text{nl}}}$
$\displaystyle=\bm{0}$ (16a)
$\displaystyle\partialderivative{F_{\text{nl}}}{\bm{\mathcal{M}}}$
$\displaystyle=\bm{0}$ (16b)
$\displaystyle\partialderivative{F_{\text{nl}}}{\bm{\Gamma}_{\text{nl}}}$
$\displaystyle=\left(\partialderivative{F_{\text{nl}}}{\bm{\mathcal{R}}_{\text{T}}},\partialderivative{F_{\text{nl}}}{\bm{\mathcal{R}}_{\text{C}}}\right)=\bm{0}$
(16c)

where the gradient in Eqs (16) depends solely on the BO energies and forces as
in the SCHA, see appendix E and Ref. [20] for details. Once the equilibrium
conditions, Eqs (16), are reached the system’s real interacting normal modes
are described, in a self-consistent framework, by the NLSCHA auxiliary phonons

$\bm{D}_{\text{nl}}=\overset{-T}{\sqrt{\bm{\mathcal{M}}}}\cdot\bm{\Phi}_{\text{nl}}\cdot\overset{-1}{\sqrt{\bm{\mathcal{M}}}}=\sum_{\mu=1}^{2}\omega_{\text{nl,}\mu}^{2}\bm{e}_{\text{nl},\mu}\bm{e}_{\text{nl},\mu}$
(17)

where $-1$ and $-T$ indicate the inverse and its transpose. The variational
optimization of its curvature $\kappa$ allows the description of both linear
vibrations (see Fig. 5 b for $\kappa\rightarrow 0$) and roto-librations (see
Fig. 5 a for $\kappa>0$) so the crystal/molecule activates the minimum-free
energy degrees of freedom without any external constraint.

##  IV Results at zero temperature

The NLSCHA results obtained solving Eqs (16) are reported in Fig. 6. In Fig. 6
a, we report the exact, SCHA and NLSCHA zero-point energy (ZPE) changing the
crystal field $E$. Fig. 6 b-i compare the exact probability distributions with
the SCHA and NLSCHA ones for three values of the crystal field $E$.

For low values of $E$, the NLSCHA outperforms the SCHA thanks to the finite
curvature $\kappa$ allowing the bending of the vibration in the angular
variable. The description of rotations is excellent, in particular, at
$E=0.005$ Ha/Bohr (see lower panels of Fig. 6) the SCHA error is $21.6$ meV
while the NLSCHA error is one order of magnitude smaller, being just $3.0$
meV. As already discussed in section II, as $E$ is increased, only linear
vibrations survive. Note that, reducing the angular fluctuation, NLSCHA
perfectly reproduces the SCHA results, Fig. 6, see appendix F.

Figure 6: Fig. a: exact, SCHA and NLSCHA zero-point energy (ZPE), i.e. the
difference between the ground state energy and the minimum of the potential.
Fig. b-e: exact, SCHA, and NLSCHA probability distributions (in Bohr-2) for
three values of the crystal field marked by horizontal lines $E=0,0.005,0.06$
Ha/Bohr. At low values of $E$, NLSCHA describes roto-librations since the
stereographic projection disentangles rotations from vibrations (Figs d-g). At
high $E$, NLSCHA reproduces the SCHA results (Figs i-l).

Figure 7: Fig. a: exact, SCHA and NLSCHA ZPE for low values of the crystal
field (from $E=0$ Ha/Bohr up to $E=5\cdot 10^{-3}$ Ha/Bohr) which do not
suppress roto-librations. Figs b-l: exact, SCHA, and NLSCHA probability
distributions (in Bohr-2) for three values of the crystal field marked by
horizontal lines in the upper panel $E=5\cdot 10^{-3},2.5\cdot 10^{-3},5\cdot
10^{-3}$ Ha/Bohr. Contrary to SHCA, NLSCHA captures the smallest changes in
the rotational degree of freedom thanks to the curvature optimization.

In Fig. 7, we compare the exact solution with SCHA and NLSCHA for extremely
low crystal field values, where the molecule rotates almost freely. Our method
detects even the subtlest alterations in the rotational degree of freedom, see
Fig. 7 b-e. Conversely, SCHA is inaccurate as it yields an almost crystal
field-independent ZPE. While NLSCHA completely captures semi-free rotations
(up to $\pi/2$) even for small values of the crystal field, it fails when
$E=0$ Ha/Bohr (see Fig. 6). Indeed, large angular fluctuations, i.e.
$\omega_{\text{nl,}\text{rot}}\rightarrow 0$, introduce spurious distortions
into the probability distribution preventing a comprehensive characterization
of free rotations, see Fig. 8.

Note that Fig. 8 has a connection with cartography. The limit
$\omega_{\text{nl,}\text{rot}}\rightarrow 0$ implies diverging angular
fluctuations, that push the probability weight at the boundaries of the
stereographic plane where the Jacobian (Eq. (11)) deforms the distribution
(Fig. 8). Similarly, the stereographic projection deforms the geographical
areas closer to the South Pole (Fig. 3). Both in cartography and NLSCHA, this
effect is due to the Jacobian (Eq. (11)) which is dominant for continents
below the Equator (Fig. 3), and large angular fluctuations (Fig. 8).
Importantly, we emphasize that a multi-peak and sharp distribution as in Fig.
8 has a very high kinetic energy, thus making it inaccessible during the
minimization of $F_{\text{nl}}$ (Eq. (15b)).

Figure 8: The NLSCHA distribution
$\bra{\bm{R}}\hat{\rho}_{\text{nl}}\ket{\bm{R}}$ (Bohr-2) for
$\omega_{\text{nl,}\text{rot}}\rightarrow 0$. The limit of vanishing angular
frequency, $\omega_{\text{nl,}\text{rot}}\rightarrow 0$, pushes the spectral
weight at the boundaries of the stereographic plane so that the NLSCHA
Jacobian (Eq. (11)) deforms the probability distribution. Hence, NLSCHA can
not describe free rotations.

##  V Results at finite temperature

Figure 9: Fig. a compares the exact free energy with the SCHA and NLSCHA
results ($E=0.01$ Ha/Bohr). Fig. b-g: exact, SCHA and NLSCHA probability
distribution (in Bohr-2) for various temperatures ($0-500-1000$ K). Finite
temperature increases the rotational degree of freedom. SCHA completely misses
the activation of this degree of freedom whereas in NLSCHA the angular
fluctuation is optimized to variationally minimize the free energy.

Temperature plays a major role in the thermodynamics of molecules, as it can
activate their rotations. Here, we investigate the thermal effect on the H2
model comparing the SCHA and NLSCHA with exact results. Fig. 9 a presents the
exact, SCHA, and NLSCHA free energies in the temperature range from $0$ K up
to $1000$ K. Furthermore, Figs 9 b-l display the exact, SCHA, and NLSCHA
probability distributions for $0$, $500$, and $1000$ K. Figs 9 b, c, d show
how the amplitude of the rotation increases upon heating.

Notably, as the temperature rises, the SCHA error also increases as the trial
probability distribution does not bend as the NLSCHA one. The NLSCHA approach
effectively captures the temperature-induced activation of rotations by
appropriately determining the optimal $\kappa$ for each temperature $T$.
Notably, at $500$ K, the NLSCHA error is $5.3$ meV while the SCHA completely
fails with an error of $23.8$ meV. The NLSCHA error grows with temperature as
the entropy is generated by a harmonic Hamiltonian and of course, there is
space for future works to relax this hypothesis. The free parameters of NLSCHA
and SCHA are reported respectively in appendix F and appendix B.

##  VI Phase transition

Figure 10: The model for III-IV phase transition of pure hydrogen. Exact,
SCHA, and NLSCHA free energy difference $\Delta
F=F_{\text{rot}}-F_{\text{vib}}+\delta$ as a function of temperature $T$.
$F_{\text{rot}}$ is the free energy for $E=0.01$ Ha/Bohr and $F_{\text{vib}}$
for $E=0.06$ Ha/Bohr. The critical temperature $T_{\text{c}}$ (defined by
$\Delta F=0$ meV) depends on the value of the shift $\delta$, which differs
for all the methods used. NLSCHA finds the correct critical temperature as it
accurately reproduces the temperature-induced activation of rotations (see
Fig. 9).

Now, we discuss a qualitative model for the III-IV phase transition of high-
pressure hydrogen. In phase III, the orientation of the H2 molecules is fixed
[23], by increasing temperature, free rotations are activated, and phase IV is
stabilized at around $300$ K [24, 25, 12].

We consider the exact, SCHA and NLSCHA free energy difference for two values
of the crystal field $E=0.01$, $0.06$ Ha/Bohr. $E=0.01$ Ha/Bohr models phase
IV as rotations are thermally activated (see Fig. 9), while $E=0.06$ Ha/Bohr
represents phase III as these modes are always locked (see appendix G). Fig.
10 shows $\Delta F=F_{\text{rot}}-F_{\text{vib}}+\delta$ where $\delta$ is a
different shift for all the methods. We define the critical temperature
$T_{\text{c}}$ for each method as the temperature where $\Delta F=0$ meV.
Remarkably, the NLSCHA $T_{\text{c}}$ is very close to the exact one with an
error of $14\%$; on the contrary, in the SCHA, there is an overestimation of
$75\%$ K due to the unphysical hybridization of the rotations with the linear
vibrations. Hence, we expect SCHA to be completely unreliable for
investigating high-pressure hydrogen at finite temperatures and should be
replaced by the more powerful NLSCHA. Remarkably, as in the SCHA, the entropy
of NLSCHA does not need any additional complex calculation, contrary to
MD/PIMD. Indeed, it is analytical and depends solely on the NLSCHA phonons
(Eq. (17)) [20].

##  VII 3D case

In this section, we extend the nonlinear change of variables discussed in
section III.1 to the case of $N$ atoms in three dimensions. Here, we employ
the stereographic projection on the sphere, so we transform the Cartesian
coordinates $\bm{R}_{i}$ of each atom $i$ in the following way

$\displaystyle\bm{R}_{i}=\bm{\mathcal{R}}_{\text{C},i}+(u_{i,1}+r_{0,i})\begin{bmatrix}\cos(\phi(\bm{u}_{i}))\sin(\theta(\bm{u}_{i}))\\\
\sin(\phi(\bm{u}_{i}))\sin(\theta(\bm{u}_{i}))\\\
\cos(\theta(\bm{u}_{i}))\end{bmatrix}$ (18a)
$\displaystyle\phi(\bm{u}_{i})=\phi_{0,i}+\atan\left(\frac{u_{i,2}}{u_{i,3}}\right)$
(18b)
$\displaystyle\theta(\bm{u}_{i})=\theta_{0,i}+2\atan\left(\frac{\sqrt{u_{i,2}^{2}+u_{i,3}^{2}}}{2r_{0,i}}\right)$
(18c)

As in section III.1, the free parameters are the center of the curvature,
$\bm{\mathcal{R}}_{\text{C},i}$, and the curvature vector,
$\bm{\mathcal{R}}_{\text{T},i}$, defined as

$\bm{\mathcal{R}}_{\text{T},i}=r_{0,i}\begin{bmatrix}\cos(\phi_{0,i})\sin(\theta_{0,i})&\sin(\phi_{0,i})\sin(\theta_{0,i})&\cos(\theta_{0,i})\end{bmatrix}$
(19)

Note that the nonlinear transformation does not mix the coordinates of
different atoms; consequently, the number of extra free parameters scales
linearly with $N$. $\bm{\mathcal{R}}_{\text{C},i}$ is a linear shift of the
positions. $\bm{\mathcal{R}}_{\text{T},i}$ defines the position of the
stereographic plane and its length gives the inverse curvature on each atom
$r_{0,i}=|\bm{\mathcal{R}}_{\text{T},i}|=\kappa_{i}^{-1}$. So $\kappa_{i}>0$
means that atom $i$ is part of a group that rotates (e.g. an organic molecule
inside the cage of a molecular perovskite) otherwise, if
$\kappa_{i}\rightarrow 0$, the atom only vibrates (e.g. the atoms of a cage in
a molecular perovskite). Note that we can not use $(\phi,\theta)$ in NLSCHA as
small oscillations are not well-defined. Indeed, $\theta\simeq 0$ imply large
$\phi$ fluctuations. On the contrary, we always define the harmonic
approximation in the stereographic coordinates (see Fig. 3). In Fig. 11, we
present the NLSCHA probability distribution for a 3D rotating diatomic
molecule in the center-of-mass reference frame $\bm{R}=\bm{R}_{1}-\bm{R}_{2}$.
Fig. 11 highlights the applicability of our approach to the study of extended
crystals and molecules.

Figure 11: 3D isosurfaces of the NLSCHA probability distribution for a
diatomic molecule’s center of mass $\bm{R}=\bm{R}_{1}-\bm{R}_{2}$. The plot
was made with mayavi package [26].

##  VIII Conclusions

We have demonstrated that NLSCHA is a very promising method for a systematic
and unbiased investigation of molecular crystals’ phase diagrams as it takes
into account linear vibrations with the same accuracy as SCHA. Still, also it
perfectly describes roto-librations, as no other currently available
approximation. The improved flexibility is due to the curvature $\kappa$,
which minimizes the free energy, allowing the system to spontaneously activate
its degrees of freedom.

It is important to emphasize that the additional computational overhead of
NLSCHA is negligible. The SCHA and NLSCHA minimize the free energy with
respect to an auxiliary force constant matrix, but in NLSCHA we also have the
mass tensor $\bm{\mathcal{M}}$. The dimension of these matrices increases
quadratically with the number of atoms, thereby being the most computationally
intensive step in the minimization. However, while the SCHA minimizes the
centroid position $\bm{\mathcal{R}}$, representing the average atomic
position, the NLSCHA optimizes two centroids $\bm{\mathcal{R}}_{\text{T}}$ and
$\bm{\mathcal{R}}_{\text{C}}$ for each atom. The dual centroid structure
arises due to the nonlinear transformation and encapsulates information about
the particle’s average position and the curvature $\kappa_{i}$. Of particular
significance is that the nonlinear transformation acts solely on single-
particle coordinates; as such, the coordinates of different atoms remain
separate (section VII). Hence, the dimensions of $\bm{\mathcal{R}}_{\text{T}}$
and $\bm{\mathcal{R}}_{\text{C}}$ scale linearly with the number of atoms.
Consequently, the NLSCHA computational cost in the thermodynamic limit scales
quadratically with the number of atoms as in the SCHA. This is primarily due
to the dominating influence of the optimization process concerning the
auxiliary force constant matrix and the mass tensor.

In conclusion, the NLSCHA equations can be solved stochastically as discussed
in Refs [20, 1]. This makes our method the most promising competitor of PIMD
as it incorporates both vibrations and roto-librational degrees of freedom in
a quantum framework, finally solving SCHA weakness pointed out in Refs [8,
12].

## Acknowledgements

A.S. and F.M. acknowledge support from European Union under project ERC-SYN
MORE-TEM (grant agreement No 951215). L.M. acknowledges funding from the
European Research Council, Marie Curie, project THERMOH.

##  Appendix A Exact diagonalization

The Morse potential for H2 is fitted using Quantum Espresso [27, 28] combined
with the BLYP exchange functional [29]. The cutoffs on plane waves and charge
density are $80$ Ry and $320$ Ry, the k-point grid is $(10,10,10)$. The
simulation box has size $20$ Bohr.

Figure 12: Fit of DFT-BLYP energy profile for H2 obtained with Quantum
Espresso.

The exact diagonalization of the Schrodinger equation is performed on a
uniform square grid in $\bm{R}$ space of size $N=1200$ between $\pm 3$ Bohr.
To get eigenvectors and eigenfunctions of $\hat{H}^{\text{(BO)}}$, we use the
Implicitly Restarted Lanczos Method [30] as implemented in the scipy function
scipy.sparse.linalg.eigsh [31]. The Lanczos algorithm applies many times the
target Hamiltonian $\hat{H}^{\text{(BO)}}$ to a starting normalized
wavefunction $\ket{\psi}_{\text{init}}$ to get the ground state and the first
excited states. To avoid the storage of the Hamiltonian as a $N^{2}\times
N^{2}$ sparse matrix we used scipy.sparse.linalg.LinearOperator which allows
to compute on-the-fly $\hat{H}^{\text{(BO)}}\ket{\psi}$ where
$\bra{\bm{R}}\ket{\psi}$ is defined on the 2D $N\times N$ grid.

##  Appendix B SCHA simulations

The SCHA simulations were performed on a uniform square grid in $\bm{R}$ space
of size $1200$ between $\pm 3$ Bohr. The conjugate gradient (CG) minimization
of the SCHA free energy was performed with the scipy [31] function
scipy.optimize.minimize setting gtol to $10^{-9}$ and maxiter $400$. We report
the SCHA free parameters’ values both at zero (appendix B.1) and finite
temperature (appendix B.2).

###  B.1 Zero temperature

In Fig. 13 we report the SCHA free parameters at equilibrium from $E=0$
Ha/Bohr to $E=0.06$ Ha/Bohr. $\mathcal{R}_{x}$ represents the $x$ component of
the SCHA centroid, $\bm{\mathcal{R}}$, while $\omega_{\text{rot}}$ and
$\omega_{\text{vib}}$ are the rotational and vibrational SCHA frequencies. We
remark that $\mathcal{R}_{y}$ is zero by symmetry.

Figure 13: The free parameters of SCHA at zero temperature from $E=0$ Ha/Bohr
to $E=0.06$ Ha/Bohr.

In Fig. 14 we report the SCHA free parameters at equilibrium at low values of
the crystal field from $E=0$ Ha/Bohr to $E=0.005$ Ha/Bohr.

Figure 14: The free parameters of SCHA at zero temperature from $E=0$ Ha/Bohr
to $E=0.005$ Ha/Bohr.

###  B.2 Finite temperature

The free parameters $\mathcal{R}_{x}$, $\omega_{\text{rot}}$ and
$\omega_{\text{vib}}$ at finite temperature for $E=0.01-0.06$ Ha/Bohr are
reported in Figs 15 16.

Figure 15: The free parameters of SCHA at finite temperature for $E=0.01$
Ha/Bohr.

Figure 16: The free parameters of SCHA at finite temperature for $E=0.06$
Ha/Bohr.

##  Appendix C The trial density matrix

According to Ref. [20], the nonlinear SCHA trial density matrix is given by
Eq. (10) where $\overline{\rho}_{\text{nl}}(\bm{u},\bm{u}^{\prime})$ is
generated by an harmonic Hamiltonian Eq. (13)

$\displaystyle\overline{\rho}_{\text{nl}}(\bm{u},\bm{u}^{\prime})=\sqrt{\det\left(\frac{\bm{\Upsilon}_{\text{nl}}}{2\pi}\right)}\exp\left\\{-\frac{1}{4}\sum_{ab=1}^{2}u_{a}\Theta_{\text{nl}}{}_{ab}u_{b}\right.$
(20)
$\displaystyle\left.-\frac{1}{4}\sum_{ab=1}^{2}u^{\prime}_{a}\Theta_{\text{nl}}{}_{ab}u^{\prime}_{b}+\sum_{ab=1}^{2}u_{a}A_{\text{nl}}{}_{ab}u^{\prime}_{b}\right\\}$

The NLSCHA tensors are related by
$\bm{\Upsilon}_{\text{nl}}=\bm{\Theta}_{\text{nl}}-2\bm{A}_{\text{nl}}$ and
are defined by

$\displaystyle\Upsilon_{\text{nl},ab}$
$\displaystyle=\sum_{ij=1}^{2}\sqrt{\mathcal{M}}_{ai}\overline{\Upsilon}_{\text{nl,ij}}\sqrt{\mathcal{M}}^{T}_{jb}$
(21a) $\displaystyle\overline{\Upsilon}_{\text{nl,ij}}$
$\displaystyle=\sum_{\mu=1}^{2}\frac{2\omega_{\text{nl,}}{}_{\mu}}{\hbar(1+2n_{\text{nl,}\mu})}e_{\text{nl},\mu}^{i}e_{\text{nl},\mu}^{j}$
(21b)

and

$\displaystyle A_{\text{nl}}{}_{,ab}$
$\displaystyle=\sum_{ij=1}^{2}\sqrt{\mathcal{M}}_{ai}\overline{A}_{\text{nl,ij}}\sqrt{\mathcal{M}}^{T}_{jb}$
(22a) $\displaystyle\overline{A}_{\text{nl,ij}}$
$\displaystyle=\sum_{\mu=1}^{2}\frac{2\omega_{\text{nl,}}{}_{\mu}n_{\text{nl,}\mu}(1+n_{\text{nl,}\mu})}{\hbar(1+2n_{\text{nl,}\mu})}e_{\text{nl},\mu}^{i}e_{\text{nl},\mu}^{j}$
(22b)

where $n_{\text{nl,}\mu}$ is the Bose-Einstein occupation number

$n_{\text{nl,}\mu}=\frac{1}{e^{\beta\hbar\omega_{\text{nl,}\mu}}-1}$ (23)

The square root of the mass tensor $\bm{\mathcal{M}}$ satisfies

$\mathcal{M}_{ab}=\sum_{i=1}^{2}\sqrt{\mathcal{M}}_{ai}\sqrt{\mathcal{M}}^{T}_{ib}$
(24)

It seems that there is a contradiction within nonlinear SCHA because
$\overline{\rho}_{\text{nl}}(\bm{u},\bm{u}^{\prime})$ is normalized as a
Gaussian [20] but the radial variable $u_{1}$ is defined between $[0,+\infty)$
(see Eq. (6)). So we approximate the normalization condition extending the
range of $u_{1}$ in $(-\infty,+\infty)$

$\displaystyle\Tr\left[\hat{\rho}_{\text{nl}}\right]$
$\displaystyle=\int_{-\infty}^{+\infty}d\bm{R}\bra{\bm{R}}\hat{\rho}_{\text{nl}}\ket{\bm{R}}$
(25)
$\displaystyle\simeq\int_{-\infty}^{+\infty}d\bm{u}\overline{\rho}_{\text{nl}}(\bm{u},\bm{u})=1$

The above assumption is justified by the following argument. In a diatomic
molecule, the linear vibration along the radial coordinate $u_{1}$ is a high
energy mode so that the bond length is always well defined hence the
corresponding component of the wave function is very localized, i.e. it decays
rapidly to zero, otherwise the atoms of the molecule will collapse one into
the other. No approximations are considered for the normalization along
$u_{2}$ as we used the stereographic projection that maps the angles
$[0,2\pi]$ in $(-\infty,+\infty)$ (see Eq. (6c)).

##  Appendix D The nonlinear SCHA free energy

In appendix D.1 we compute all the necessary quantities to get the NLSCHA
kinetic energy according to Ref. [20]. In appendix D.2 we present the full
NLSCHA free energy following Ref. [20].

###  D.1 Preliminary definitions

We define the Jacobian of $\bm{\xi}$ (Eq. (6))

$J^{i}_{j}=\partialderivative{R_{i}}{u_{j}}=\begin{bmatrix}\partialderivative{x}{u_{1}}&\partialderivative{y}{u_{1}}\\\
\partialderivative{x}{u_{2}}&\partialderivative{y}{u_{2}}\end{bmatrix}=\begin{bmatrix}\frac{x-x_{\text{C}}}{u_{1}+r_{0}}&\frac{y-y_{\text{C}}}{u_{1}+r_{0}}\\\
-\frac{y-y_{\text{C}}}{r_{0}f^{2}(u_{2})}&-\frac{x-x_{\text{C}}}{r_{0}f^{2}(u_{2})}\end{bmatrix}$
(26)

where $f(u_{2})$ is

$f(u_{2})=\sqrt{1+\left(\frac{u_{2}}{2r_{0}}\right)^{2}}$ (27)

The determinat of Eq. (26)

$\mathcal{J}=\left|\det\left(\bm{J}\right)\right|=\frac{u_{1}+r_{0}}{r_{0}f^{2}(u_{2})}$
(28)

The inverse metric tensor $\bm{g}$ is

$g{}^{ab}=\sum_{i=1}^{2}\partialderivative{u_{a}}{R_{i}}\partialderivative{u_{b}}{R_{j}}=\begin{bmatrix}1&0\\\
0&(r_{0}f(u_{2}))^{2}h(\bm{u})\end{bmatrix}$ (29)

where $h(\bm{u})$ is

$h(\bm{u})=\left(\frac{f(u_{2})}{u_{1}+r_{0}}\right)^{2}$ (30)

In addition, we define the vector $\bm{d}$ as

$\bm{d}=\frac{1}{2}\partialderivative{\log(\mathcal{J})}{\bm{u}}=\begin{bmatrix}\frac{1}{2(u_{1}+r_{0})}\\\
-\frac{u_{2}}{(2r_{0}f(u_{2}))^{2}}\end{bmatrix}$ (31)

###  D.2 Free energy calculation

According to Ref. [20] the nonlinear SCHA kinetic energy is

$\Tr\left[\frac{\hat{\bm{P}}^{2}}{2m}\hat{\rho}_{\text{nl}}\right]=\left\langle\mathcal{K}\right\rangle_{\text{nl}}=\int_{-\infty}^{+\infty}d\bm{u}\mathcal{K}\overline{\rho}_{\text{nl}}(\bm{u})$
(32)

$\mathcal{K}$ is the kinetic kernel [20]

$\mathcal{K}=\sum_{a=1}^{2}\left(\mathcal{K}^{(2)}_{a}\mathcal{L}^{(2)}_{a}+\mathcal{K}^{(1)}_{a}\mathcal{L}^{(1)}_{a}\right)+\mathcal{K}^{(0)}$
(33)

where $\bm{\mathcal{L}}^{(1)}$, $\bm{\mathcal{L}}^{(2)}$ are defined by

$\displaystyle\mathcal{L}^{(1)}_{a}$
$\displaystyle=-\partialderivative{\log\left(\overline{\rho}_{\text{nl}}(\bm{u})\right)}{\widetilde{u}_{a}}=\frac{1}{\sqrt{m}}\sum_{i=1}^{3N}\Upsilon_{\text{nl},ai}u_{i}$
(34a) $\displaystyle\mathcal{L}^{(2)}_{a}$
$\displaystyle=-\partialderivative{\log\left(\overline{\rho}_{\text{nl}}(\bm{u})\right)}{\widetilde{u}_{a}}{\widetilde{u}_{a}}$
(34b)
$\displaystyle=\frac{1}{m}\left(\Upsilon_{\text{nl},aa}-\sum_{ij=1}^{3N}\Upsilon_{\text{nl},ai}\Upsilon_{\text{nl},aj}u_{i}u_{j}\right)$

and $\mathcal{K}^{(0)}$, $\bm{\mathcal{K}}^{(1)}$, and
$\bm{\mathcal{K}}^{(2)}$ by

$\displaystyle\mathcal{K}^{(0)}$
$\displaystyle=\frac{\hbar^{2}}{2m}\left\\{\Tr\left[\bm{g}\cdot\left(\frac{\bm{\Upsilon}_{\text{nl}}}{4}+\bm{A}_{\text{nl}}\right)\right]+\bm{d}\cdot\bm{g}\cdot\bm{d}\right\\}$
$\displaystyle=\frac{\hbar^{2}}{2m}\Tr\left[\bm{g}\cdot\left(\frac{\bm{\Upsilon}_{\text{nl}}}{4}+\bm{A}_{\text{nl}}\right)\right]+\frac{\hbar^{2}h(\bm{u})}{8m}$
(35a) $\displaystyle\bm{\mathcal{K}}^{(1)}$
$\displaystyle=\frac{\hbar^{2}}{4}\bm{g}\cdot\bm{d}=\frac{\hbar^{2}}{4\sqrt{m}}\begin{bmatrix}\frac{1}{2(u_{1}+r_{0})}\\\
-\frac{1}{4}u_{2}h(\bm{u})\end{bmatrix}$ (35b)
$\displaystyle\bm{\mathcal{K}}^{(2)}$
$\displaystyle=-\frac{\hbar^{2}}{8}\text{diag}(\bm{g})=-\frac{\hbar^{2}}{8}\begin{bmatrix}1\\\
(r_{0}f(u_{2}))^{2}h(\bm{u})\end{bmatrix}$ (35c)

where $\bm{g}$ and $\bm{d}$ are defined in Eqs (29) (31). The potential energy
is

$\displaystyle\Tr\left[\hat{V}^{\text{(BO)}}\hat{\rho}_{\text{nl}}\right]$
$\displaystyle=\int_{-\infty}^{+\infty}d\bm{u}V^{\text{(BO)}}(\bm{\xi}(\bm{u}))\overline{\rho}_{\text{nl}}(\bm{u})$
(36) $\displaystyle=\left\langle V^{\text{(BO)}}\right\rangle_{\text{nl}}$

So the NLSCHA free energy is

$F_{\text{nl}}=\left\langle\mathcal{K}\right\rangle_{\text{nl}}+\left\langle
V^{\text{(BO)}}\right\rangle_{\text{nl}}-TS_{\text{nl}}$ (37)

where the entropy is harmonic as discussed in Ref. [20]

$\displaystyle
S_{\text{nl}}=\Tr\left[\hat{\rho}_{\text{nl}}\log(\hat{\rho}_{\text{nl}})\right]$
(38)
$\displaystyle=k_{\text{B}}\sum_{\mu=1}^{2}\left[(1+n_{\text{nl,}\mu})\log(1+n_{\text{nl,}\mu})-n_{\text{nl,}\mu}\log(n_{\text{nl,}\mu})\right]$

##  Appendix E Gradient of the nonlinear SCHA free energy

We start with the gradient with respect to the center of the curvature
$\bm{\mathcal{R}}_{\text{C}}$ (Eq. (7a))

$\partialderivative{F_{\text{nl}}}{\bm{\mathcal{R}}_{\text{C}}}=\left\langle\partialderivative{V^{\text{(BO)}}}{\bm{R}}\right\rangle_{\text{nl}}$
(39)

The gradient with respect to curvature vector $\bm{\mathcal{R}}_{\text{T}}$
(Eq. (7b)) is the following

$\displaystyle\partialderivative{F_{\text{nl}}}{\bm{\mathcal{R}}_{\text{T}}}=\left\langle\partialderivative{\mathcal{K}}{\bm{\mathcal{R}}_{\text{T}}}\right\rangle_{\text{nl}}+\left\langle\partialderivative{V^{\text{(BO)}}}{\bm{\mathcal{R}}_{\text{T}}}\right\rangle_{\text{nl}}$
(40)
$\displaystyle=\frac{\bm{\mathcal{R}}_{\text{T}}}{|\bm{\mathcal{R}}_{\text{T}}|}\left\langle\partialderivative{\mathcal{K}}{r_{0}}\right\rangle_{\text{nl}}+\sum_{i=1}^{2}\partialderivative{R_{i}}{\bm{\mathcal{R}}_{\text{T}}}\left\langle\partialderivative{V^{\text{(BO)}}}{{R_{i}}}\right\rangle_{\text{nl}}$

noting that $\mathcal{K}(\bm{u})$ (Eq. (33)) depends only on the curvature
$|\bm{\mathcal{R}}_{\text{T}}|=r_{0}=\kappa^{-1}$. The derivative of the
kinetic energy kernel $\mathcal{K}$ (Eq. (33)) is

$\displaystyle\partialderivative{\mathcal{K}}{r_{0}}=\sum_{a=1}^{2}\left(\partialderivative{\mathcal{K}^{(2)}_{a}}{r_{0}}\mathcal{L}^{(2)}_{a}+\partialderivative{\mathcal{K}^{(1)}_{a}}{r_{0}}\mathcal{L}^{(1)}_{a}\right)+\partialderivative{\mathcal{K}^{(0)}(\bm{u})}{r_{0}}(\bm{u})$
(41)

where $r_{0}$ the derivative of the coefficients (Eqs (35)) are

$\displaystyle\partialderivative{\mathcal{K}^{(0)}}{r_{0}}=$
$\displaystyle\frac{\hbar^{2}}{2m}\Tr\left[\partialderivative{\bm{g}}{r_{0}}\cdot\left(\frac{\bm{\Upsilon}_{\text{nl}}}{4}+\bm{A}_{\text{nl}}\right)\right]$
$\displaystyle+\frac{\hbar^{2}}{8m}\partialderivative{h(\bm{u})}{r_{0}}$ (42a)
$\displaystyle\partialderivative{\bm{\mathcal{K}}^{(1)}}{r_{0}}=$
$\displaystyle\frac{\hbar^{2}}{4\sqrt{m}}\begin{bmatrix}\partialderivative{r_{0}}\frac{1}{2(u_{1}+r_{0})}\\\
-\frac{1}{4}u_{2}\partialderivative{r_{0}}h(\bm{u})\end{bmatrix}$ (42b)
$\displaystyle\partialderivative{\bm{\mathcal{K}}^{(2)}}{r_{0}}=$
$\displaystyle-\frac{\hbar^{2}}{8}\begin{bmatrix}0\\\
\partialderivative{[(r_{0}f(u_{2}))^{2}h(\bm{u})]}{r_{0}}\end{bmatrix}$ (42c)
$\displaystyle\partialderivative{\bm{g}}{r_{0}}=$
$\displaystyle\begin{bmatrix}0&0\\\
0&\partialderivative{[(r_{0}f(u_{2}))^{2}h(\bm{u})]}{r_{0}}\end{bmatrix}$
(42d) $\displaystyle\partialderivative{f^{2}(u_{2})}{r_{0}}=$
$\displaystyle-\frac{u_{2}{}^{2}}{2r_{0}^{3}}$ (42e)
$\displaystyle\partialderivative{h(\bm{u})}{r_{0}}=$
$\displaystyle\frac{1}{(u_{1}+r_{0})^{2}}\partialderivative{f^{2}(u_{2})}{r_{0}}-\frac{2h(\bm{u})}{r_{0}+u_{1}}$
(42f)

The gradient of the Cartesian coordinates $\bm{R}$ with respect to the
curvature vector $\bm{\mathcal{R}}_{\text{T}}$ is

$\partialderivative{R_{i}}{\bm{\mathcal{R}}_{\text{T}}}=\frac{\bm{\mathcal{R}}_{\text{T}}}{|\bm{\mathcal{R}}_{\text{T}}|}\partialderivative{R_{i}}{r_{0}}+\partialderivative{\phi_{0}}{\bm{\mathcal{R}}_{\text{T}}}\partialderivative{R_{i}}{\phi_{0}}$
(43)

where the derivatives of the Cartesian coordinates $\bm{R}=(x,y)$ are

$\begin{bmatrix}\partialderivative{x}{r_{0}}&\partialderivative{y}{r_{0}}\\\
\partialderivative{x}{\phi_{0}}&\partialderivative{y}{\phi_{0}}\end{bmatrix}=\begin{bmatrix}\frac{x-x_{\text{C}}}{u_{1}+r_{0}}+\frac{u_{2}(y-y_{\text{C}})}{r_{0}^{2}f^{2}(u_{2})}&\frac{y-y_{\text{C}}}{u_{1}+r_{0}}-\frac{u_{2}(x-x_{\text{C}})}{r_{0}^{2}f^{2}(u_{2})}\\\
-(y-y_{\text{C}})&+(x-x_{\text{C}})\end{bmatrix}$ (44)

and the derivative of $\phi_{0}$ is computed from the definition of
$\bm{\mathcal{R}}_{\text{T}}$ (Eq. (7b))

$\partialderivative{\phi_{0}}{\bm{\mathcal{R}}_{\text{T}}}=-\frac{1}{|\bm{\mathcal{R}}_{\text{T}}|^{2}}\begin{bmatrix}\mathcal{R}_{\text{T},2}&-\mathcal{R}_{\text{T},1}\end{bmatrix}$
(45)

To compute the gradient with respect to $\bm{\Phi}_{\text{nl}}$ we use the
following formula introduced by Ref. [32]

$\displaystyle\partialderivative{\left\langle
O\right\rangle_{\text{nl}}}{\bm{\Phi}_{\text{nl}}}=\left\langle\partialderivative{O}{\bm{\Phi}_{\text{nl}}}\right\rangle_{\text{nl}}$
(46)
$\displaystyle+\frac{1}{2}\sum_{ijk=1}^{2}\partialderivative{\overset{-1}{\Upsilon}_{\text{nl},ij}}{\bm{\Phi}_{\text{nl}}}\Upsilon_{\text{nl},ik}\left\langle
u_{k}\partialderivative{O}{u_{j}}\right\rangle_{\text{nl}}$

The first term of Eq. (46) takes into account the explicit dependence of $O$
on $\bm{\Phi}_{\text{nl}}$ while the second considers the change in the
probability distribution $\overline{\rho}_{\text{nl}}(\bm{u})$. With Eq. (46),
the derivative of $F_{\text{nl}}$ with respect to the auxiliary force constant
matrix $\bm{\Phi}_{\text{nl}}$ is

$\displaystyle\partialderivative{F_{\text{nl}}}{\bm{\Phi}_{\text{nl}}}=\left\langle\partialderivative{\mathcal{K}(\bm{u})}{\bm{\Phi}_{\text{nl}}}\right\rangle_{\text{nl}}$
(47)
$\displaystyle+\frac{1}{2}\sum_{ijk=1}^{2}\partialderivative{\overset{-1}{\Upsilon}_{\text{nl},ij}}{\bm{\Phi}_{\text{nl}}}\Upsilon_{\text{nl},ik}\left\langle
u_{k}\left(\partialderivative{\mathcal{K}}{u_{j}}+\partialderivative{V^{\text{(BO)}}}{u_{j}}\right)\right\rangle_{\text{nl}}$
$\displaystyle-T\partialderivative{S_{\text{nl}}}{\bm{\Phi}_{\text{nl}}}$

where derivative of the entropy $S_{\text{nl}}$ with respect to
$\bm{\Phi}_{\text{nl}}$ is computed in Ref. [20].

The first term is the derivative of $\mathcal{K}(\bm{u})$ (Eq. (33)) with
respect to $\bm{\Phi}_{\text{nl}}$

$\partialderivative{\mathcal{K}}{\bm{\Phi}_{\text{nl}}}=\sum_{a=1}^{2}\left(\mathcal{K}^{(2)}_{a}\partialderivative{\mathcal{L}^{(2)}_{a}}{\bm{\Phi}_{\text{nl}}}+\mathcal{K}^{(1)}_{a}\partialderivative{\mathcal{L}^{(1)}_{a}}{\bm{\Phi}_{\text{nl}}}\right)+\partialderivative{\mathcal{K}^{(0)}}{\bm{\Phi}_{\text{nl}}}$
(48)

The derivatives of $\bm{\mathcal{L}}^{(1)}$, $\bm{\mathcal{L}}^{(2)}$ (Eqs
(34)) are

$\displaystyle\partialderivative{\mathcal{L}^{(1)}_{a}}{\bm{\Phi}_{\text{nl}}}$
$\displaystyle=\frac{1}{\sqrt{m}}\left(\sum_{i=1}^{2}\partialderivative{\Upsilon_{\text{nl},ai}}{\bm{\Phi}_{\text{nl}}}u{}_{i}\right)$
(49a)
$\displaystyle\partialderivative{\mathcal{L}^{(2)}_{a}}{\bm{\Phi}_{\text{nl}}}$
$\displaystyle=\frac{1}{m}\left(\partialderivative{\Upsilon_{\text{nl},aa}}{\bm{\Phi}_{\text{nl}}}-2\sum_{ij=1}^{2}\partialderivative{\Upsilon_{\text{nl},ai}}{\bm{\Phi}_{\text{nl}}}\Upsilon_{\text{nl},aj}u{}_{i}u{}_{j}\right)$
(49b)

The derivative of $\mathcal{K}^{(0)}$ (Eq. (35)) is

$\partialderivative{\mathcal{K}^{(0)}}{\bm{\Phi}_{\text{nl}}}=\frac{\hbar^{2}}{2m}\Tr\left[\bm{g}\cdot\left(\frac{1}{4}\partialderivative{\bm{\Upsilon}_{\text{nl}}}{\bm{\Phi}_{\text{nl}}}+\partialderivative{\bm{A}_{\text{nl}}}{\bm{\Phi}_{\text{nl}}}\right)\right]$
(50)

The derivatives of $\bm{\Upsilon}_{\text{nl}}$/$\bm{A}_{\text{nl}}$ can be
computed with the expressions of Ref. [20].

The derivative of the kinetic energy kernel $\mathcal{K}(\bm{u})$ with respect
to the auxiliary displacements $\bm{u}$ is

$\displaystyle\partialderivative{\mathcal{K}}{\bm{u}}=$
$\displaystyle\sum_{a=1}^{2}\left(\partialderivative{\mathcal{K}^{(2)}_{a}}{\bm{u}}\mathcal{L}^{(2)}_{a}+\mathcal{K}^{(2)}_{a}\partialderivative{\mathcal{L}^{(2)}_{a}}{\bm{u}}\right)$
(51)
$\displaystyle+\sum_{a=1}^{2}\left(\partialderivative{\mathcal{K}^{(1)}_{a}}{\bm{u}}\mathcal{L}^{(1)}_{a}+\mathcal{K}^{(1)}_{a}\partialderivative{\mathcal{L}^{(1)}_{a}}{\bm{u}}\right)+\partialderivative{\mathcal{K}^{(0)}}{\bm{u}}$

The derivatives with respect to $\bm{u}$ of $\mathcal{K}^{(0)}$,
$\bm{\mathcal{K}}^{(1)}$, and $\bm{\mathcal{K}}^{(2)}$ (Eqs (35)) are

$\displaystyle\partialderivative{\mathcal{K}^{(0)}}{\bm{u}}=$
$\displaystyle\frac{\hbar^{2}}{2m}\left\\{\Tr\left[\partialderivative{\bm{g}}{\bm{u}}\cdot\left(\frac{\bm{\Upsilon}_{\text{nl}}}{4}+\bm{A}_{\text{nl}}\right)\right]\right.$
$\displaystyle\left.+\frac{1}{4}\partialderivative{h(\bm{u})}{\bm{u}}\right\\}$
(52a) $\displaystyle\partialderivative{\bm{\mathcal{K}}^{(1)}}{\bm{u}}=$
$\displaystyle\frac{\hbar^{2}}{4\sqrt{m}}\begin{bmatrix}\partialderivative{\bm{u}}\frac{1}{2(u_{1}+r_{0})}\\\
-\frac{1}{4}u_{2}\partialderivative{\bm{u}}h(\bm{u})\end{bmatrix}$ (52b)
$\displaystyle\partialderivative{\bm{\mathcal{K}}^{(2)}}{\bm{u}}=$
$\displaystyle-\frac{\hbar^{2}}{8}\begin{bmatrix}0\\\
r_{0}^{2}\partialderivative{(f^{2}(u_{2})h(\bm{u}))}{\bm{u}}\end{bmatrix}$
(52c) $\displaystyle\partialderivative{\bm{g}}{\bm{u}}=$
$\displaystyle\begin{bmatrix}0&0\\\
0&r_{0}^{2}\partialderivative{[f^{2}(u_{2})h(\bm{u})]}{\bm{u}}\end{bmatrix}$
(52d) $\displaystyle\partialderivative{f^{2}(u_{2})}{\bm{u}}=$
$\displaystyle\begin{bmatrix}0\\\ \frac{u_{2}}{2r_{0}^{2}}\end{bmatrix}$ (52e)
$\displaystyle\partialderivative{h(\bm{u})}{\bm{u}}=$
$\displaystyle\begin{bmatrix}-\frac{2}{(r_{0}+u_{1})}h(\bm{u})\\\
\frac{1}{(r_{0}+u_{1})^{2}}\partialderivative{f^{2}(u_{2})}{u_{2}}\end{bmatrix}$
(52f)

and the derivatives of $\bm{\mathcal{L}}^{(1)}$, $\bm{\mathcal{L}}^{(2)}$ (Eqs
(34)) are

$\displaystyle\partialderivative{\mathcal{L}^{(1)}_{a}}{u_{b}}=$
$\displaystyle\frac{1}{\sqrt{m}}\Upsilon_{\text{nl},ab}$ (53a)
$\displaystyle\partialderivative{\mathcal{L}^{(2)}_{a}}{u_{b}}=$
$\displaystyle\frac{1}{m}\left(-2\Upsilon_{\text{nl},ab}\sum_{i=1}^{2}\Upsilon_{\text{nl},ai}u_{i}\right)$
(53b)

The derivative of the BOES with respect to the auxiliary displacements
$\bm{u}$ is

$\partialderivative{V^{\text{(BO)}}(\bm{R})}{u_{j}}=-\sum_{i=1}^{2}J^{i}_{j}f^{\text{(BO)}}_{i}$
(54)

where $\bm{J}$ is defined in Eq. (26).

We minimize with respect to $\sqrt{\bm{\mathcal{M}}}$ so that we can enforce
the positive definiteness of
$\bm{\mathcal{M}}=\sqrt{\bm{\mathcal{M}}}^{T}\cdot\sqrt{\bm{\mathcal{M}}}$

$\displaystyle\partialderivative{F_{\text{nl}}}{\sqrt{\bm{\mathcal{M}}}}=\left\langle\partialderivative{\mathcal{K}}{\sqrt{\bm{\mathcal{M}}}}\right\rangle_{\text{nl}}$
(55)
$\displaystyle+\frac{1}{2}\sum_{ijk=1}^{2}\partialderivative{\overset{-1}{\Upsilon}_{\text{nl},ij}}{\sqrt{\bm{\mathcal{M}}}}\Upsilon_{\text{nl},ik}\left\langle
u_{k}\left(\partialderivative{\mathcal{K}}{u_{j}}+\partialderivative{V^{\text{(BO)}}}{u_{j}}\right)\right\rangle_{\text{nl}}$
$\displaystyle-T\partialderivative{S_{\text{nl}}}{\sqrt{\bm{\mathcal{M}}}}$

This equation is a straightforward generalization of Eq. (47). The derivative
$\overset{-1}{\bm{\Upsilon}}_{\text{nl}}$ is computed in Ref. [20]. The
derivative of the kinetic kernel is

$\partialderivative{\mathcal{K}}{\sqrt{\bm{\mathcal{M}}}}=\sum_{a=1}^{2}\left(\mathcal{K}^{(2)}_{a}\partialderivative{\mathcal{L}^{(2)}_{a}}{\sqrt{\bm{\mathcal{M}}}}+\mathcal{K}^{(1)}_{a}\partialderivative{\mathcal{L}^{(1)}_{a}}{\sqrt{\bm{\mathcal{M}}}}\right)+\partialderivative{\mathcal{K}^{(0)}}{\sqrt{\bm{\mathcal{M}}}}$
(56)

Then, to compute

$\displaystyle\partialderivative{\mathcal{L}^{(2)}_{a}}{\sqrt{\bm{\mathcal{M}}}}$
$\displaystyle=\partialderivative{\Upsilon_{\text{nl},aa}}{\sqrt{\bm{\mathcal{M}}}}-2\sum_{ij=1}^{2}\partialderivative{\Upsilon_{\text{nl},ai}}{\sqrt{\bm{\mathcal{M}}}}\Upsilon_{\text{nl},aj}u{}_{i}u{}_{j}$
(57a)
$\displaystyle\partialderivative{\mathcal{L}^{(1)}_{a}}{\sqrt{\bm{\mathcal{M}}}}$
$\displaystyle=\sum_{i=1}^{2}\partialderivative{\Upsilon_{\text{nl},ai}}{\sqrt{\bm{\mathcal{M}}}}u{}_{i}$
(57b)

and

$\partialderivative{\mathcal{K}^{(0)}}{\sqrt{\bm{\mathcal{M}}}}=\frac{\hbar^{2}}{2m}\Tr\left[\bm{g}\cdot\left(\frac{1}{4}\partialderivative{\bm{\Upsilon}_{\text{nl}}}{\sqrt{\bm{\mathcal{M}}}}+\partialderivative{\bm{A}_{\text{nl}}}{\sqrt{\bm{\mathcal{M}}}}\right)\right]$
(58)

we employ the formulas derived in Ref. [20]. The derivative of the entropy is

$\displaystyle\partialderivative{S_{\text{nl}}}{\sqrt{\mathcal{M}}_{ab}}=\sum_{ij=1}^{2}\partialderivative{S_{\text{nl}}}{D_{\text{nl}}{}_{,ij}}\partialderivative{D_{\text{nl}}{}_{,ij}}{\sqrt{\mathcal{M}}_{ab}}$
(59)

where $\partialderivative{\bm{D}_{\text{nl}}}{\sqrt{\bm{\mathcal{M}}}}$ and
$\partialderivative{S_{\text{nl}}}{\bm{D}_{\text{nl}}}$ are given by Ref.
[20].

##  Appendix F Nonlinear SCHA simulations

The NLSCHA simulations were performed on a uniform square grid in $\bm{u}$
space of size $1000$ between $\pm 3$ Bohr. The conjugate gradient (CG)
minimization of the NLSCHA free energy was performed with the scipy [31]
function scipy.optimize.minimize setting gtol to $10^{-8}$ and maxiter $4000$.
We report the free parameters of NLSCHA both at zero (appendix F.1) and finite
temperature (appendix F.2).

###  F.1 Parameters at zero temperature

In Figs 17 18 we report the nonlinear SCHA free parameters for $E=0-E=0.06$
Ha/Bohr and $E=0-E=0.005$ Ha/Bohr. We show the $x$-component of
$\bm{\mathcal{R}}_{\text{C}}$, denoted by $x_{\text{C}}$, the radius of the
curvature, $r_{0}$, the rotational and vibrational frequencies,
$\omega_{\text{rot}}$ and $\omega_{\text{vib}}$, with the corresponding
effective mass eigenvalues,
$\mathcal{M}_{\text{rot}},\mathcal{M}_{\text{vib}}$, divided by the physical
value $\mu_{\ch{H2}}=m_{\ch{H}}/2$. The symmetry of the potential, see Eq.
(1), fixes the values of $y_{\text{C}}=0$ Bohr and $\phi_{0}=\pi$ are fixed.
Note that at $T=0$ K we find that $\bm{\mathcal{M}}$ coincides with the
physical one.

Figure 17: The free parameters of NLSCHA
($X_{C},r_{0},\omega_{\text{rot}},\omega_{\text{vib}},\mathcal{M}_{\text{rot}}/\mu_{\ch{H2}},\mathcal{M}_{\text{vib}}/\mu_{\ch{H2}}$)
at $T=0$ K from $E=0$ Ha/Bohr to $E=0.06$ Ha/Bohr.

Figure 18: The free parameters of NLSCHA
($X_{C},r_{0},\omega_{\text{rot}},\omega_{\text{vib}},\mathcal{M}_{\text{rot}}/\mu_{\ch{H2}},\mathcal{M}_{\text{vib}}/\mu_{\ch{H2}}$)
at $T=0$ K and low values of the crystal field $E$ from $E=0$ Ha/Bohr to
$E=0.005$ Ha/Bohr.

###  F.2 Parameters at finite temperature

Figs 19 20 report the NLSCHA free parameters at finite temperature from $T=0$
K to $T=1000$ K for $E=0.01$ Ha/Bohr and $E=0.06$ Ha/Bohr. For $E=0.01$
Ha/Bohr, the eigenvalue of $\bm{\mathcal{M}}$ corresponding to the rotational
mode, $\mathcal{M}_{\text{rot}}$, is lower than the physical one so that
rotational fluctuations are enhanced. Note that for the vibron mode, the
eigenvalue of $\bm{\mathcal{M}}$ is
$\mathcal{M}_{\text{vib}}\simeq\mu_{\ch{H2}}$. For $E=0.06$ Ha/Bohr,
$\mathcal{M}_{\text{vib}}<\mu_{\ch{H2}}$ at $200$ K compensates for the
increase of $\omega_{\text{rot}}$.

Figure 19: The free parameters of NLSCHA
($X_{C},r_{0},\omega_{\text{rot}},\omega_{\text{vib}},\mathcal{M}_{\text{rot}}/\mu_{\ch{H2}},\mathcal{M}_{\text{vib}}/\mu_{\ch{H2}}$)
from $T=0$ K to $T=1000$ K with $E=0.01$ Ha/Bohr.

Figure 20: The free parameters of NLSCHA
($X_{C},r_{0},\omega_{\text{rot}},\omega_{\text{vib}},\mathcal{M}_{\text{rot}}/\mu_{\ch{H2}},\mathcal{M}_{\text{vib}}/\mu_{\ch{H2}}$)
from $T=0$ K to $T=1000$ K with $E=0.06$ Ha/Bohr.

In Fig. 21 we compare the NLSCHA free energies ($E=0.01$ Ha/Bohr) obtained
optimizing $\bm{\mathcal{M}}$ and keeping it fixed
$\bm{\mathcal{M}}=\mu_{\ch{H2}}\bm{1}$. Fig. 21 b-l report the corresponding
probability distributions compared with the exact result. Note that at
$T=1000$ K optimizing $\bm{\mathcal{M}}$ gives an error of $8.5$ meV while
keeping it fixed yields $9.3$ meV.

Figure 21: Fig. a shows the exact, and NLSCHA free energies as a function of
temperature for $E=0.01$ Ha/Bohr. We report the NLSCHA solutions with
$\bm{\mathcal{M}}=\mu_{\text{\ch{H2}}}\bm{1}$ and
$\bm{\mathcal{M}}\neq\mu_{\text{\ch{H2}}}\bm{1}$. Figs b-l report the exact
(Figs b-e-h), and NLSCHA (Figs c-f-i with
$\bm{\mathcal{M}}\neq\mu_{\text{\ch{H2}}}\bm{1}$ and Figs d-g-l with
$\bm{\mathcal{M}}=\mu_{\text{\ch{H2}}}\bm{1}$) probability distribution
$\rho(R)$ (Bohr-2) at $0-500-1000$ K.

##  Appendix G Phase transition model

In Fig. 22-23 we report the exact, harmonic, SCHA and NLSCHA free energies for
$E=0.01-0.06$ Ha/Bohr form $T=0$ K to $1000-600$ K along with the probability
distributions. Note that when rotations are locked ($E=0.06$ Ha/Bohr) the
NLSCHA solution coincides with the SCHA.

Figure 22: Fig. a shows the exact, harmonic, SCHA and NLSCHA free energies as
a function of temperature for $E=0.01$ Ha/Bohr. Figs b-o report the exact
(Figs b-f-l), harmonic (Figs d-h-n), SCHA (Figs c-g-m) and NLSCHA (Figs e-i-o)
probability distribution $\rho(\bm{R})$ (Bohr-2) at $0-500-1000$ K.

Figure 23: Fig. a shows the exact, harmonic, SCHA and NLSCHA free energies as
a function of temperature for $E=0.06$ Ha/Bohr. Figs b-o report the exact
(Figs b-f-l), harmonic (Figs d-h-n), SCHA (Figs c-g-m) and NLSCHA (Figs e-i-o)
probability distribution $\rho(\bm{R})$ (Bohr-2) at $0-300-600$ K.

## References

  * Monacelli  _et al._ [2021] L. Monacelli, R. Bianco, M. Cherubini, M. Calandra, I. Errea, and F. Mauri, The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects, Journal of Physics: Condensed Matter 33, 363001 (2021).
  * Oba  _et al._ [2019] Y. Oba, T. Tadano, R. Akashi, and S. Tsuneyuki, First-principles study of phonon anharmonicity and negative thermal expansion in ${\mathrm{scf}}_{3}$, Phys. Rev. Materials 3, 033601 (2019).
  * Hellman  _et al._ [2011] O. Hellman, I. A. Abrikosov, and S. I. Simak, Lattice dynamics of anharmonic solids from first principles, Phys. Rev. B 84, 180301 (2011).
  * Souvatzis  _et al._ [2008] P. Souvatzis, O. Eriksson, M. I. Katsnelson, and S. P. Rudin, Entropy driven stabilization of energetically unstable crystal structures explained from first principles theory, Phys. Rev. Lett. 100, 095901 (2008).
  * Monserrat  _et al._ [2013] B. Monserrat, N. D. Drummond, and R. J. Needs, Anharmonic vibrational properties in periodic systems: energy, electron-phonon coupling, and stress, Phys. Rev. B 87, 144302 (2013).
  * Gambhir  _et al._ [1999] M. Gambhir, M. T. Dove, and V. Heine, Rigid unit modes and dynamic disorder: Sio2 cristobalite and quartz, Physics and Chemistry of Minerals 26, 484 (1999).
  * Giddy  _et al._ [1993] A. P. Giddy, M. T. Dove, G. S. Pawley, and V. Heine, The determination of rigid-unit modes as potential soft modes for displacive phase transitions in framework crystal structures, Acta Crystallographica Section A 49, 697 (1993), https://onlinelibrary.wiley.com/doi/pdf/10.1107/S0108767393002545 .
  * Kapil  _et al._ [2019] V. Kapil, E. Engel, M. Rossi, and M. Ceriotti, Assessment of approximate methods for anharmonic free energies, Journal of Chemical Theory and Computation 15, 5845 (2019), pMID: 31532997, https://doi.org/10.1021/acs.jctc.9b00596 .
  * Rossi  _et al._ [2016] M. Rossi, P. Gasparotto, and M. Ceriotti, Anharmonic and quantum fluctuations in molecular crystals: A first-principles study of the stability of paracetamol, Phys. Rev. Lett. 117, 115702 (2016).
  * Kitamura  _et al._ [2000] H. Kitamura, S. Tsuneyuki, T. Ogitsu, and T. Miyake, Quantum distribution of protons in solid molecular hydrogen at megabar pressures, Nature 404, 259 (2000).
  * Kohanoff  _et al._ [1997] J. Kohanoff, S. Scandolo, G. L. Chiarotti, and E. Tosatti, Solid molecular hydrogen: The broken symmetry phase, Phys. Rev. Lett. 78, 2783 (1997).
  * Morresi  _et al._ [2022] T. Morresi, R. Vuilleumier, and M. Casula, Hydrogen phase-iv characterization by full account of quantum anharmonicity, Phys. Rev. B 106, 054109 (2022).
  * Georgescu and Mandelshtam [2012] I. Georgescu and V. A. Mandelshtam, Self-consistent phonons revisited. i. the role of thermal versus quantum fluctuations on structural transitions in large lennard-jones clusters, The Journal of Chemical Physics 137, 144106 (2012), https://doi.org/10.1063/1.4754819 .
  * Brown  _et al._ [2013a] S. E. Brown, I. Georgescu, and V. A. Mandelshtam, Self-consistent phonons revisited. ii. a general and efficient method for computing free energies and vibrational spectra of molecules and clusters, The Journal of Chemical Physics 138, 044317 (2013a), https://doi.org/10.1063/1.4788977 .
  * Tadano and Tsuneyuki [2018] T. Tadano and S. Tsuneyuki, First-principles lattice dynamics method for strongly anharmonic crystals, Journal of the Physical Society of Japan 87, 041015 (2018), https://doi.org/10.7566/JPSJ.87.041015 .
  * Brown  _et al._ [2013b] S. E. Brown, I. Georgescu, and V. A. Mandelshtam, Self-consistent phonons revisited. II. A general and efficient method for computing free energies and vibrational spectra of molecules and clusters, The Journal of Chemical Physics 138, 044317 (2013b), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4788977/15461101/044317_1_online.pdf .
  * Monacelli and Mauri [2021] L. Monacelli and F. Mauri, Time-dependent self-consistent harmonic approximation: Anharmonic nuclear quantum dynamics and time correlation functions, Phys. Rev. B 103, 104305 (2021).
  * Siciliano  _et al._ [2023] A. Siciliano, L. Monacelli, G. Caldarelli, and F. Mauri, Wigner gaussian dynamics: Simulating the anharmonic and quantum ionic motion, Phys. Rev. B 107, 174307 (2023).
  * Lihm and Park [2021] J.-M. Lihm and C.-H. Park, Gaussian time-dependent variational principle for the finite-temperature anharmonic lattice dynamics, Phys. Rev. Research 3, L032017 (2021).
  * Siciliano  _et al._ [2024] A. Siciliano, L. Monacelli, and F. Mauri, Beyond gaussian fluctuations of quantum anharmonic nuclei (2024), arXiv:2407.03802 [cond-mat.mtrl-sci] .
  * Hooton [1955] D. Hooton, Li. a new treatment of anharmonicity in lattice thermodynamics: I, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 46, 422 (1955), https://doi.org/10.1080/14786440408520575 .
  * Met Office [2015] Met Office,  _Cartopy: a cartographic python library with a Matplotlib interface_, Exeter, Devon (2010 - 2015).
  * Monacelli  _et al._ [2023] L. Monacelli, M. Casula, K. Nakano, S. Sorella, and F. Mauri, Quantum phase diagram of high-pressure hydrogen, Nature Physics 10.1038/s41567-023-01960-5 (2023).
  * Howie  _et al._ [2012a] R. T. Howie, C. L. Guillaume, T. Scheler, A. F. Goncharov, and E. Gregoryanz, Mixed molecular and atomic phase of dense hydrogen, Phys. Rev. Lett. 108, 125501 (2012a).
  * Howie  _et al._ [2012b] R. T. Howie, T. Scheler, C. L. Guillaume, and E. Gregoryanz, Proton tunneling in phase iv of hydrogen and deuterium, Phys. Rev. B 86, 214104 (2012b).
  * Ramachandran and Varoquaux [2011] P. Ramachandran and G. Varoquaux, Mayavi: 3D Visualization of Scientific Data, Computing in Science & Engineering 13, 40 (2011).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with quantum ESPRESSO, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Miehlich  _et al._ [1989] B. Miehlich, A. Savin, H. Stoll, and H. Preuss, Results obtained with the correlation energy density functionals of becke and lee, yang and parr, Chemical Physics Letters 157, 200 (1989).
  * Lehoucq  _et al._ [1998] R. B. Lehoucq, D. C. Sorensen, and C. Yang,  _ARPACK Users’ Guide_ (Society for Industrial and Applied Mathematics, 1998) https://epubs.siam.org/doi/pdf/10.1137/1.9780898719628 .
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261 (2020).
  * Bianco  _et al._ [2017] R. Bianco, I. Errea, L. Paulatto, M. Calandra, and F. Mauri, Second-order structural phase transitions, free energy curvature, and temperature-dependent anharmonic phonons in the self-consistent harmonic approximation: Theory and stochastic implementation, Phys. Rev. B 96, 014111 (2017).