  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 thermodynamics
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 \textit{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.
