An ab initio quantum-mechanical theoretical framework is presented to compute the thermal properties of molecular crystals. The present strategy combines dispersion-corrected density-functional-theory (DFT-D), harmonic phonon dispersion, quasi-harmonic approximation to the lattice dynamics for thermal expansion and thermodynamic functions, and quasi-static approximation for anisotropic thermoelasticity. The proposed scheme is shown to reliably describe thermal properties of the urea molecular crystal by a thorough comparison with experimental data.Molecular crystals have increasingly attracted great attention due to their peculiar chemical and physical properties, which make them suitable as high energy-density materials, 1-3 active pharmaceutical ingredients (APIs), 4-6 constituents of opto-electronic devices for their linear and non-linear optical properties, 7-9 etc.Nonetheless, from a theoretical view-point, they still represent a major challenge to the state-of-the-art quantum-chemical methods as many kinds of chemical interactions (covalent intra-molecular, electrostatic, hydrogen-bond, long-range dispersive) need to be accurately described simultaneously. Only in recent years, have different theoretical approaches been devised in order to predict their structural and energetic properties (with the main goal of discriminating between competing polymorphs): from force-field to high-level molecular fragment-based schemes, from periodic dispersion-corrected density functional theory (DFT-D) to periodic many-body wave-function techniques. [10][11][12][13][14][15][16][17][18] However, once a reliable and balanced description of the various chemical interactions has been achieved by means of any of the above-mentioned quantum-chemical methods, the extension of their applicability to more complex properties of technological and industrial relevance, which would greatly increase their predictiveness, such as mechanical, elastic, optical and thermodynamic responses, [19][20][21][22] has to be performed. Apart from the intrinsic high degree of complexity of the required theoretical techniques and algorithms, the main difficulty here is represented by the fact that most of those properties are largely affected by thermal effects, 23-25 even at room temperature, such as zero-point energy, harmonic and anharmonic thermal nuclear motion, thermal lattice expansion, etc.Most quantum-chemical ab initio methods describe the groundstate of a system at zero pressure and temperature. If the inclusion of pressure in the computed structural and elastic properties is a relatively easy task, 16,[26][27][28] this is definitely not the case when temperature has to be accurately accounted for. Indeed, we are still far from having effective schemes formally developed and efficiently implemented in a solid state context, particularly so when anharmonic terms to the lattice potential have to be included into the formalism. When the harmonic approximation (HA) to the lattice potential is used, the vibrational contribution to the free ene...