Predictive Modeling of Paracetamol Polymorphs: Linking Basis Set and Hamiltonian Effects to an Optimized Quasi-Harmonic Lattice Dynamics Workflow

Abstract

Thermodynamic properties of molecular crystals are the cornerstone of a broad range of applications, spanning pharmaceutics, organic semiconductors, and many more. The reliable modeling of these properties is usually hindered by collective influences of delicate dispersion interactions and insufficient prior knowledge. Providing reasonable accuracy and minimal reliance on experimental data, lattice dynamics (LD) based on the quasi-harmonic approximation (QHA) and density functional theory (DFT) has been widely and successfully adopted for crystal structure prediction (CSP). However, for polymorphic molecular crystals, it remains challenging to accurately establish the relative stability among multiple polymorphs that are distinguished by subtle energy differences. Despite recent intense efforts to achieve calculations that can reproduce observations in well-characterized systems, a comprehensive study elucidating how predictions are governed by errors within the theoretical framework and details of implementation is lacking. This undermines confidence in the transferability and efficiency of QHA LD workflows proposed in the literature. The need for reliable and predictive calculations necessitates systematic investigation and disentangling of the influences from key parameters such as the choice of basis set and Hamiltonian, so the cost and precision of DFT-based QHA LD calculations can be rationally balanced in practical workflows. In this study, paracetamol polymorphs are adopted as the prototype system, as they are well-characterized and display representative chemical interactions, where covalent bonds, hydrogen bonds, and London dispersions are present with sufficient complexity. The intermixed sources of error underlying free-energy predictions are resolved by comparing the structural, vibrational, and thermochemical properties of paracetamol polymorphs with high-quality references. A hierarchy of atomic and plane-wave basis sets are adopted to highlight the significance of basis set completeness, and results with and without the inclusion of Fock exchange identify the particular importance of Hamiltonian for thermochemical properties. Furthermore, an optimal QHA LD workflow is proposed and tested utilizing these insights, which combines results from plane-wave basis set with semilocal functional and local, well-converged def2-TZVP basis set with hybrid functional to efficiently converge the numerical precision for QHA LD. While demonstrated in detail for paracetamol, conclusions of this study are also expected to be of relevance to a wide range of molecular crystals. These findings rigorously set the stage for affordable and reliable QHA LD workflow for molecular crystals, where fortuitous error cancellations are clearly distinguished from converged levels of theory, thereby clarifying the numerical precision attainable before explicit anharmonic and configurational effects are considered for complex organic systems.