Interest in molecular crystals has grown thanks to their relevance to pharmaceuticals, organic semiconductor materials, foods, and many other applications. Electronic structure methods have become an increasingly important tool for modeling molecular crystals and polymorphism. This article reviews electronic structure techniques used to model molecular crystals, including periodic density functional theory, periodic second-order Møller-Plesset perturbation theory, fragment-based electronic structure methods, and diffusion Monte Carlo. It also discusses the use of these models for predicting a variety of crystal properties that are relevant to the study of polymorphism, including lattice energies, structures, crystal structure prediction, polymorphism, phase diagrams, vibrational spectroscopies, and nuclear magnetic resonance spectroscopy. Finally, tools for analyzing crystal structures and intermolecular interactions are briefly discussed.

Organic crystals frequently adopt multiple distinct polymorphs exhibiting different properties. The ability to predict not only what crystal forms might occur, but under what experimental thermodynamic conditions those polymorphs are stable would be immensely valuable to the pharmaceutical industry and others. Starting only from knowledge of the experimental crystal structures, this study successfully predicts the methanol crystal polymorph phase diagram from first-principles quantum chemistry, mapping out the thermodynamic regions of stability for three polymorphs over the range 0-400 K and 0-6 GPa. The agreement between the predicted and experimental phase diagrams corresponds to predicting the relative polymorph free energies to within ∼0.5 kJ mol-1 accuracy, which is achieved by employing fragment-based second-order Møller-Plesset perturbation theory and coupled cluster theory plus a quasi-harmonic treatment of the phonons.

Molecular crystals expand appreciably upon heating due to both zero-point and thermal vibrational motion, yet this expansion is often neglected in molecular crystal modeling studies. Here, a quasi-harmonic approximation is coupled with fragment-based hybrid many-body interaction calculations to predict thermal expansion and finite-temperature thermochemical properties in crystalline carbon dioxide, ice Ih, acetic acid and imidazole. Fragment-based second-order Möller-Plesset perturbation theory (MP2) and coupled cluster theory with singles, doubles and perturbative triples [CCSD(T)] predict the thermal expansion and the temperature dependence of the enthalpies, entropies and Gibbs free energies of sublimation in good agreement with experiment. The errors introduced by neglecting thermal expansion in the enthalpy and entropy cancel somewhat in the Gibbs free energy. The resulting ∼ 1-2 kJ mol(-1) errors in the free energy near room temperature are comparable to or smaller than the errors expected from the electronic structure treatment, but they may be sufficiently large to affect free-energy rankings among energetically close polymorphs.

Ab initio predictions of chemical shifts and electric field gradient (EFG) tensor components are frequently used to help interpret solid-state nuclear magnetic resonance (NMR) experiments. Typically, these predictions employ density functional theory (DFT) with generalized gradient approximation (GGA) functionals, though hybrid functionals have been shown to improve accuracy relative to experiment. Here, the performance of a dozen models beyond the GGA approximation are examined for the prediction of solid-state NMR observables, including meta-GGA, hybrid, and double-hybrid density functionals and second-order Møller-Plesset perturbation theory (MP2). These models are tested on organic molecular crystal data sets containing 169 experimental ¹³C and ¹⁵N chemical shifts and 114 ¹⁷O and ¹⁴N EFG tensor components. To make these calculations affordable, gauge-including projector augmented wave (GIPAW) Perdew-Burke-Ernzerhof (PBE) calculations with periodic boundary conditions are combined with a local intramolecular correction computed at the higher level of theory. Within the context of typical NMR property calculations performed on a static, DFT-optimized crystal structure, the benchmarking finds that the double-hybrid DFT functionals produce errors versus experiment that are no smaller than those of hybrid functionals in the best cases, and they can be larger. MP2 errors versus experiment are even bigger. Overall, no practical advantages are found for using any of the tested double-hybrid functionals or MP2 to predict experimental solid-state NMR chemical shifts and EFG tensor components for routine organic crystals, especially given the higher computational cost of those methods. This finding likely reflects error cancellation benefiting the hybrid functionals. Improving the accuracy of the predicted chemical shifts and EFG tensors relative to experiment would probably require more robust treatments of the crystal structures, their dynamics, and other factors.

Molecular crystal structures, thermodynamics, and mechanical properties can vary substantially with temperature, and predicting these temperature-dependencies correctly is important for many practical applications in the pharmaceutical industry and other fields. However, most electronic structure predictions of molecular crystal properties neglect temperature and/or thermal expansion, leading to potentially erroneous results. Here, we demonstrate that by combining large basis set second-order Møller-Plesset (MP2) or even coupled cluster singles, doubles, and perturbative triples (CCSD(T)) electronic structure calculations with a quasiharmonic treatment of thermal expansion, experimentally observable properties such as the unit cell volume, heat capacity, enthalpy, entropy, sublimation point and bulk modulus of phase I crystalline carbon dioxide can be predicted in excellent agreement with experiment over a broad range of temperatures. These results point toward a promising future for ab initio prediction of molecular crystal properties at real-world temperatures and pressures.

Chemical shift prediction plays an important role in the determination or validation of crystal structures with solid-state nuclear magnetic resonance (NMR) spectroscopy. One of the fundamental theoretical challenges lies in discriminating variations in chemical shifts resulting from different crystallographic environments. Fragment-based electronic structure methods provide an alternative to the widely used plane wave gauge-including projector augmented wave (GIPAW) density functional technique for chemical shift prediction. Fragment methods allow hybrid density functionals to be employed routinely in chemical shift prediction, and we have recently demonstrated appreciable improvements in the accuracy of the predicted shifts when using the hybrid PBE0 functional instead of generalized gradient approximation (GGA) functionals like PBE. Here, we investigate the solid-state ¹³C and ¹⁵N NMR spectra for multiple crystal forms of acetaminophen, phenobarbital, and testosterone. We demonstrate that the use of the hybrid density functional instead of a GGA provides both higher accuracy in the chemical shifts and increased discrimination among the different crystallographic environments. Finally, these results also provide compelling evidence for the transferability of the linear regression parameters mapping predicted chemical shieldings to chemical shifts that were derived in an earlier study.

The dispersion-corrected second-order Møller-Plesset perturbation theory (MP2C) approach accurately describes intermolecular interactions in many systems. MP2C, however, expends much computational effort to compute the long-range correlation with MP2, only to discard and replace those contributions with a simpler long-range dispersion correction based on intermolecular perturbation theory. Here, we demonstrate that one can avoid calculating the long-range MP2 correlation by attenuating the Coulomb operator, allowing the dispersion correction to handle the long-range interactions inexpensively. With relatively modest Coulomb attenuation, one obtains results that are very similar to those from conventional MP2C. With more aggressive attenuation, one can remove just enough short-range repulsive exchange-dispersion interactions to compensate for finite basis set errors. Doing so makes it possible to approach complete basis set limit quality results with only an aug-cc-pVTZ basis, resulting in substantial computational savings. Further computational savings could be achieved by reformulating the MP2C algorithm to exploit the increased sparsity of the two-electron integrals.

The ability to predict not only what organic crystal structures might occur but also the thermodynamic conditions under which they are the most stable would be extremely useful for discovering and designing new organic materials. The present study takes a step in that direction by predicting the temperature- and pressure-dependent phase boundary between the α and β polymorphs of resorcinol using density functional theory (DFT) and the quasi-harmonic approximation. To circumvent the major computational bottleneck associated with computing a well-converged phonon density of states via the supercell approach, a recently developed approximation is employed, which combines a supercell phonon density of states from dispersion-corrected third-order density functional tight binding [DFTB3-D3(BJ)] with frequency corrections derived from a smaller B86bPBE-XDM functional DFT phonon calculation on the crystallographic unit cell. This mixed DFT/DFTB quasi-harmonic approach predicts the lattice constants and unit cell volumes to within 1%-2% at lower pressures. It predicts the thermodynamic phase boundary in almost perfect agreement with the experiment, although this excellent agreement does reflect fortuitous cancellation of errors between the enthalpy and entropy of transition.

The performance of fragment-based ab initio(1)H, (13)C, (15)N and (17)O chemical shift predictions is assessed against experimental NMR chemical shift data in four benchmark sets of molecular crystals. Employing a variety of commonly used density functionals (PBE0, B3LYP, TPSSh, OPBE, PBE, TPSS), we explore the relative performance of cluster, two-body fragment, and combined cluster/fragment models. The hybrid density functionals (PBE0, B3LYP and TPSSh) generally out-perform their generalized gradient approximation (GGA)-based counterparts. (1)H, (13)C, (15)N, and (17)O isotropic chemical shifts can be predicted with root-mean-square errors of 0.3, 1.5, 4.2, and 9.8 ppm, respectively, using a computationally inexpensive electrostatically embedded two-body PBE0 fragment model. Oxygen chemical shieldings prove particularly sensitive to local many-body effects, and using a combined cluster/fragment model instead of the simple two-body fragment model decreases the root-mean-square errors to 7.6 ppm. These fragment-based model errors compare favorably with GIPAW PBE ones of 0.4, 2.2, 5.4, and 7.2 ppm for the same (1)H, (13)C, (15)N, and (17)O test sets. Using these benchmark calculations, a set of recommended linear regression parameters for mapping between calculated chemical shieldings and observed chemical shifts are provided and their robustness assessed using statistical cross-validation. We demonstrate the utility of these approaches and the reported scaling parameters on applications to 9-tert-butyl anthracene, several histidine co-crystals, benzoic acid and the C-nitrosoarene SnCl2(CH3)2(NODMA)2.