The need for configurational sampling dramatically increases the cost of combined quantum mechanics/molecular mechanics (QM/MM) simulations of chemical processes in solution. We developed an averaged condensed phase environment (ACPE) model that constructs an effective polarizable environment directly from explicitly sampled molecular dynamics configurations via the K-means++ algorithm and a mathematically rigorous translation of the molecular mechanics parameters. The model captures detailed heterogeneous features in the environment that may be difficult to describe using a conventional polarizable continuum model. Instead of performing repeated QM/MM calculations for each new configuration of the environment, the ACPE approach allows one to perform a single QM calculation on an averaged configuration. Here, we demonstrate the model by computing electronic excitation energies for several small molecules in solution. The ACPE model predicts the excitation energies in excellent agreement with conventional configurational averaging yet with orders of magnitude of reduction in the computational cost.

Fragment-based methods predict nuclear magnetic resonance (NMR) chemical shielding tensors in molecular crystals with high accuracy and computational efficiency. Such methods typically employ electrostatic embedding to mimic the crystalline environment, and the quality of the results can be sensitive to the embedding treatment. To improve the quality of this embedding environment for fragment-based molecular crystal property calculations, we borrow ideas from the embedded ion method to incorporate self-consistently polarized Madelung field effects. The self-consistent reproduction of the Madelung potential (SCRMP) model developed here constructs an array of point charges that incorporates self-consistent lattice polarization and which reproduces the Madelung potential at all atomic sites involved in the quantum mechanical region of the system. The performance of fragment- and cluster-based 1H, 13C, 14N, and 17O chemical shift predictions using SCRMP and density functionals like PBE and PBE0 are assessed. The improved embedding model results in substantial improvements in the predicted 17O chemical shifts and modest improvements in the 15N ones. Finally, the performance of the model is demonstrated by examining the assignment of the two oxygen chemical shifts in the challenging γ-polymorph of glycine. Overall, the SCRMP-embedded NMR chemical shift predictions are on par with or more accurate than those obtained with the widely used gauge-including projector augmented wave (GIPAW) model.

Solid carbon dioxide exhibits a rich phase diagram at high pressures. Metastable phase III is formed by compressing dry ice above ∼10–12 GPa. Phase VII occurs at similar pressures but higher temperatures, and its stability region is disconnected from III on the phase diagram. Comparison of large-basis-set quasi-harmonic second-order Møller–Plesset perturbation theory calculations and experiment suggests that the long-accepted structure of phase III is problematic. The experimental phase III and VII structures both relax to the same phase VII structure. Furthermore, Raman spectra predicted for phase VII are in good agreement with those observed experimentally for both phase III and VII, while those for the purported phase III structure agree poorly with experimental observations. Crystal structure prediction is employed to search for other potential structures which might account for phase III, but none are found. Together, these results suggest that phases III and VII are likely identical.

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.

Molecular crystals occur widely in pharmaceuticals, foods, explosives, organic semiconductors, and many other applications. Thanks to substantial progress in electronic structure modeling of molecular crystals, attention is now shifting from basic crystal structure prediction and lattice energy modeling toward the accurate prediction of experimentally observable properties at finite temperatures and pressures. This Account discusses how fragment-based electronic structure methods can be used to model a variety of experimentally relevant molecular crystal properties. First, it describes the coupling of fragment electronic structure models with quasi-harmonic techniques for modeling the thermal expansion of molecular crystals, and what effects this expansion has on thermochemical and mechanical properties.Excellent agreement with experiment is demonstrated for the molar volume, sublimation enthalpy, entropy, and free energy, and the bulk modulus of phase I carbon dioxide when large basis second-order Møller–Plesset perturbation theory (MP2) or coupled cluster theories (CCSD(T)) are used. In addition, physical insight is offered into how neglect of thermal expansion affects these properties. Zero-point vibrational motion leads to an appreciable expansion in the molar volume; in carbon dioxide, it accounts for around 30% of the overall volume expansion between the electronic structure energy minimum and the molar volume at the sublimation point. In addition, because thermal expansion typically weakens the intermolecular interactions, neglecting thermal expansion artificially stabilizes the solid and causes the sublimation enthalpy to be too large at higher temperatures. Thermal expansion also frequently weakens the lower-frequency lattice phonon modes; neglecting thermal expansion causes the entropy of sublimation to be overestimated. Interestingly, the sublimation free energy is less significantly affected by neglecting thermal expansion because the systematic errors in the enthalpy and entropy cancel somewhat.Second, because solid state nuclear magnetic resonance (NMR) plays an increasingly important role in molecular crystal studies, this Account discusses how fragment methods can be used to achieve higher-accuracy chemical shifts in molecular crystals. Whereas widely used plane wave density functional theory models are largely restricted to generalized gradient approximation (GGA) functionals like PBE in practice, fragment methods allow the routine use of hybrid density functionals with only modest increases in computational cost. In extensive molecular crystal benchmarks, hybrid functionals like PBE0 predict chemical shifts with 20–30% higher accuracy than GGAs, particularly for 1H, 13C, and 15N nuclei. Due to their higher sensitivity to polarization effects, 17O chemical shifts prove slightly harder to predict with fragment methods. Nevertheless, the fragment model results are still competitive with those from GIPAW.The improved accuracy achievable with fragment approaches and hybrid density functionals increases discrimination between different potential assignments of individual shifts or crystal structures, which is critical in NMR crystallography applications. This higher accuracy and greater discrimination are highlighted in application to the solid state NMR of different acetaminophen and testosterone crystal forms.

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 13C and 15N 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 performance of fragment-based ab initio1H, 13C, 15N and 17O 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. 1H, 13C, 15N, and 17O 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 1H, 13C, 15N, and 17O 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.

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.

Many-body noncovalent interactions are increasingly important in large and/or condensed-phase systems, but the current understanding of how well various models predict these interactions is limited. Here, benchmark complete-basis set coupled cluster singles, doubles, and perturbative triples (CCSD(T)) calculations have been performed to generate a new test set for three-body intermolecular interactions. This “3B-69” benchmark set includes three-body interaction energies for 69 total trimer structures, consisting of three structures from each of 23 different molecular crystals. By including structures that exhibit a variety of intermolecular interactions and packing arrangements, this set provides a stringent test for the ability of electronic structure methods to describe the correct physics involved in the interactions. Both MP2.5 (the average of second- and third-order Møller–Plesset perturbation theory) and spin-component-scaled CCSD for noncovalent interactions (SCS-MI-CCSD) perform well. MP2 handles the polarization aspects reasonably well, but it omits three-body dispersion. In contrast, many widely used density functionals corrected with three-body D3 dispersion correction perform comparatively poorly. The primary difficulty stems from the treatment of exchange and polarization in the functionals rather than from the dispersion correction, though the three-body dispersion may also be moderately underestimated by the D3 correction.

Controlling charge transfer at a molecular scale is critical for efficient light harvesting, energy conversion, and nanoelectronics. Dipole-polarization electrets, the electrostatic analogue of magnets, provide a means for “steering” electron transduction via the local electric fields generated by their permanent electric dipoles. Here, we describe the first demonstration of the utility of anthranilamides, moieties with ordered dipoles, for controlling intramolecular charge transfer. Donor–acceptor dyads, each containing a single anthranilamide moiety, distinctly rectify both the forward photoinduced electron transfer and the subsequent charge recombination. Changes in the observed charge-transfer kinetics as a function of media polarity were consistent with the anticipated effects of the anthranilamide molecular dipoles on the rectification. The regioselectivity of electron transfer and the molecular dynamics of the dyads further modulated the observed kinetics, particularly for charge recombination. These findings reveal the underlying complexity of dipole-induced effects on electron transfer and demonstrate unexplored paradigms for molecular rectifiers.

First-principles chemical shielding tensor predictions play a critical role in studying molecular crystal structures using nuclear magnetic resonance. Fragment-based electronic structure methods have dramatically improved the ability to model molecular crystal structures and energetics using high-level electronic structure methods. Here, a many-body expansion fragment approach is applied to the calculation of chemical shielding tensors in molecular crystals. First, the impact of truncating the many-body expansion at different orders and the role of electrostatic embedding are examined on a series of molecular clusters extracted from molecular crystals. Second, the ability of these techniques to assign three polymorphic forms of the drug sulfanilamide to the corresponding experimental 13C spectra is assessed. This challenging example requires discriminating among spectra whose 13C chemical shifts differ by only a few parts per million (ppm) across the different polymorphs. Fragment-based PBE0/6-311+G(2d,p) level chemical shielding predictions correctly assign these three polymorphs and reproduce the sulfanilamide experimental 13C chemical shifts with 1 ppm accuracy. The results demonstrate that fragment approaches are competitive with the widely used gauge-invariant projector augmented wave (GIPAW) periodic density functional theory calculations.

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.

Two copper borates, CuB2O4 and Cu3B2O6, produce either H2 or O2 under visible light in the presence of sacrificial agents. Both copper borates exhibit intrinsic midgap states between the valence bands and conduction bands. However, structural differences in the two compounds lead to midgap states with remarkably different photocatalytic behaviors. In CuB2O4, the midgap states facilitate visible light absorption and enhance photocatalytic activity, while in Cu3B2O6 they act as trap states that reduce photocatalytic activity.

Powder neutron diffraction and Raman spectroscopy experiments for ice XV, the recently discovered proton-ordered polymorph of ice VI, suggest that the protons arrange in an antiferroelectric structure with P1̅ symmetry, contrary to several density functional theory predictions of a ferroelectric Cc structure. Here, we find that higher-level fragment-based second-order perturbation theory (MP2) and coupled cluster theory (CCSD(T)) electronic structure calculations predict that the experimentally proposed proton ordering is indeed slightly more stable than the other possible structures. These calculations reveal a close competition between the structure with the strongest local hydrogen bonding (Cc) and the one with the most favorable “delocalized” hydrogen bond cooperativity effects (P1̅), with the latter being preferred by ∼0.4 kJ/mol per molecule. The results reiterate the importance of viewing ice networks as a whole instead of focusing on pairwise hydrogen-bonding interactions.Keywords: fragment methods; ice; molecular crystal; order−disorder; polymorphism;

Crystalline oxalyl dihydrazide has five experimentally known polymorphs whose energetics are governed by subtle balances between intra- and intermolecular interactions, providing a severe challenge for theoretical crystal structure modeling. Previous work has shown that many common density functional methods that neglect van der Waals dispersion cannot correctly describe this system, but it has been argued that empirically dispersion-corrected DFT-D performs much better. Here, we examine these crystals with second-order Møller–Plesset perturbation theory (MP2) and related levels of theory using the fragment-based hybrid many-body interaction method. The energetics prove sensitive to the treatment of electron–electron correlation, the basis set, many-body induction, three-body dispersion, and zero-point contributions. Nevertheless, our best predictions for the polymorph energy ordering based on dispersion-corrected MP2C calculations agree with the available experimental data. In contrast, lower levels of theory, including the common B3LYP-D* and D-PW91 dispersion-corrected density functional approximations, fail to reproduce experimental observations and/or the high-level calculations.

A graphical abstract is available for this content

Significant advances in fragment-based electronic structure methods have created a real alternative to force-field and density functional techniques in condensed-phase problems such as molecular crystals. This perspective article highlights some of the important challenges in modeling molecular crystals and discusses techniques for addressing them. First, we survey recent developments in fragment-based methods for molecular crystals. Second, we use examples from our own recent research on a fragment-based QM/MM method, the hybrid many-body interaction (HMBI) model, to analyze the physical requirements for a practical and effective molecular crystal model chemistry. We demonstrate that it is possible to predict molecular crystal lattice energies to within a couple kJ mol−1 and lattice parameters to within a few percent in small-molecule crystals. Fragment methods provide a systematically improvable approach to making predictions in the condensed phase, which is critical to making robust predictions regarding the subtle energy differences found in molecular crystals.

We perform the first high-level ab initio calculations (MP2) on crystalline aspirin using a newly developed fragment-based QM/MM method. Contrary to earlier density functional theory predictions, the two polymorphs are virtually degenerate, which is consistent with experimentally observed intergrowth structures. This near-degeneracy arises “accidentally” from a competition between intramolecular relaxation (form I) and intermolecular hydrogen bonding (form II).

We combine quantum and classical mechanics in a fragment-based many-body interaction model to predict organic molecular crystal lattice energies. Individual molecules in the central unit cell and their short-range pairwise interactions are modeled quantum mechanically, while long-range pairwise and many-body interactions are approximated classically. The classical contributions are evaluated using an accurate ab initio force field that is constructed on-the-fly from quantum mechanical calculations on the individual molecules in the unit cell. The force field parameters include ab initio distributed multipole moments, distributed polarizabilities, and isotropic two- and three-body atomic dispersion coefficients. This QM/MM fragment model reproduces full periodic MP2 lattice energies to within a couple kJ/mol at substantially reduced cost. When high-level electronic structure methods are coupled with the ab initio force field, molecular crystal lattice energies are predicted to within 2 kJ/mol of their experimental values for six of the seven crystals examined here. Finally, Axilrod–Teller–Muto three-body dispersion energy plays a nontrivial role in several of the molecular crystals studied here.

Quantum and classical mechanics are combined in a hybrid many-body interaction model to enable the computationally affordable study of systems containing many interacting molecules. This model treats intramolecular and pairwise intermolecular interactions quantum mechanically, while many-body electrostatic induction effects are approximated using a polarizable force field. In this paper, we demonstrate that parametrizing the force field with distributed multipoles and atom-centered polarizabilities obtained on-the-fly from ab initio quantum mechanical monomer calculations makes the model very accurate and eliminates nearly all empiricism. Test calculations on water, formamide, hydrogen fluoride, and glycine−water clusters, all of which exhibit strong many-body interactions, are presented. The performance of the hybrid model is competitive with related point-charge embedding models.

Using cyclic voltammetry, we examined the dependence of the reduction potentials of six quinones on the concentration of the supporting electrolyte. An increase in the electrolyte concentration, resulting in an increase in the solution polarity, caused positive shifts of the reduction potentials. We ascribed the observed changes in the potentials to the dependence of the solvation energy of the quinones and their anions on the media polarity. Analysis of the reduction potentials, using the Born solvation energy equation, yielded unfeasibly small values for the effective radii of the quinone species, that is, the experimentally obtained effective radii were up to 4-fold smaller than the radii of the solvation cavities that we calculated for the quinones. The nonspherical shapes of the quinones, along with the uneven charge density distribution in their anions, encompassed the underlying reasons for the discrepancies between the obtained experimental and theoretical values for the radii of these redox species. The generalized Born approach, which does not treat the solvated species as single spheres, provided means for addressing this discrepancy and yielded effective radii that were relatively close to the measured values.

A fast, fragment-based hybrid many-body interaction model is used to optimize the structures of five small-molecule organic crystals (with fixed experimental lattice parameters) and predict their lattice energies with accuracies of ∼2−4 kJ/mol compared to experiment. This model treats individual molecules in the central unit cell and their short-range two-body interactions quantum mechanically, while long-range electrostatics and many-body induction are treated with a classical polarizable force field. For the hydrogen bonded ice, formamide, and acetamide crystals, MP2 calculations extrapolated to the complete-basis-set limit provide good accuracy. However, MP2 exhibits difficulties for crystals such as benzene and imidazole, where π-stacking dispersion interactions are important, and post-MP2 corrections determined from small-basis-set CCSD(T) calculations are required to achieve chemical accuracy. Using these techniques, accurate crystal lattice energy predictions for small-molecule organic crystals are feasible with currently available computing power.Keywords (keywords): fragment methods; molecular crystals; structure prediction;

Significant advances in fragment-based electronic structure methods have created a real alternative to force-field and density functional techniques in condensed-phase problems such as molecular crystals. This perspective article highlights some of the important challenges in modeling molecular crystals and discusses techniques for addressing them. First, we survey recent developments in fragment-based methods for molecular crystals. Second, we use examples from our own recent research on a fragment-based QM/MM method, the hybrid many-body interaction (HMBI) model, to analyze the physical requirements for a practical and effective molecular crystal model chemistry. We demonstrate that it is possible to predict molecular crystal lattice energies to within a couple kJ mol−1 and lattice parameters to within a few percent in small-molecule crystals. Fragment methods provide a systematically improvable approach to making predictions in the condensed phase, which is critical to making robust predictions regarding the subtle energy differences found in molecular crystals.

Two copper borates, CuB2O4 and Cu3B2O6, produce either H2 or O2 under visible light in the presence of sacrificial agents. Both copper borates exhibit intrinsic midgap states between the valence bands and conduction bands. However, structural differences in the two compounds lead to midgap states with remarkably different photocatalytic behaviors. In CuB2O4, the midgap states facilitate visible light absorption and enhance photocatalytic activity, while in Cu3B2O6 they act as trap states that reduce photocatalytic activity.

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.

The performance of fragment-based ab initio1H, 13C, 15N and 17O 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. 1H, 13C, 15N, and 17O 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 1H, 13C, 15N, and 17O 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.