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