  Molecular crystals often exist in multiple competing polymorphs, showing
significantly different physico-chemical properties. Computational crystal
structure prediction is key to interpret and guide the search for the most
stable or useful form: A real challenge due to the combinatorial search space,
and the complex interplay of subtle effects that work together to determine the
relative stability of different structures. Here we take a comprehensive
approach based on different flavors of thermodynamic integration in order to
estimate all contributions to the free energies of these systems with
density-functional theory, including the oft-neglected anharmonic contributions
and nuclear quantum effects. We take the two main stable forms of paracetamol
as a paradigmatic example. We find that anharmonic contributions, different
descriptions of van der Waals interactions, and nuclear quantum effects all
matter to quantitatively determine the stability of different phases. Our
analysis highlights the many challenges inherent in the development of a
quantitative and predictive framework to model molecular crystals. However, it
also indicates which of the components of the free energy can benefit from a
cancellation of errors that can redeem the predictive power of approximate
models, and suggests simple steps that could be taken to improve the
reliability of ab initio crystal structure prediction.
