Predicting structure of molecular crystals from first principles

Massively Parallel Computational Chemistry with the Super Instruction Architecture and ACES4

The task of developing high-performing parallel software must be made easier and more cost-effective in order to fully exploit existing and emerging large-scale computer systems for the advancement of science. The Super Instruction Architecture (SIA) is a parallel programming platform geared toward applications that need to manage large amounts of data stored in potentially sparse multidimensional arrays during calculations. The SIA platform was originally designed for the quantum chemistry software package ACESIII. More recently, the SIA was reimplemented to overcome the limitations in the original ACESIII program. It has now been successfully employed in the new ACES4 quantum chemistry software package. This paper describes the SIA and ACES4 and illustrates their capabilities with some difficult quantum chemistry open-shell coupled-cluster benchmark calculations.

Benchmarking Swaths of Intermolecular Interaction Components with Symmetry-Adapted Perturbation Theory

A benchmark database for interaction energy components of various noncovalent interactions (NCIs) along their dissociation curve is one of the essential needs in theoretical chemistry, especially for the development of force fields and machine-learning methods. We utilize DFT-SAPT or SAPT(DFT) as one of the most accurate methods to generate an extensive stock of the energy components, including dispersion energies extrapolated to the complete basis set limit (CBS). Precise analyses of the created data, and benchmarking the total interaction energies against the best available CCSD(T)/CBS values, reveal different aspects of the methodology and the nature of NCIs. For example, error cancellation effects between the S2 approximation and nonexact xc-potentials occur, and large charge transfer energies in some systems, including heavy atoms, can explain the lower accuracy of DFT-SAPT. This method is perfect for neutral complexes containing light nonmetals, while other systems with heavier atoms should be treated carefully. In the last part, a representative data set for all NCIs is extracted from the original data.

Frontiers of molecular crystal structure prediction for pharmaceuticals and functional organic materials

The reliability of organic molecular crystal structure prediction has improved tremendously in recent years. Crystal structure predictions for small, mostly rigid molecules are quickly becoming routine. Structure predictions for larger, highly flexible molecules are more challenging, but their crystal structures can also now be predicted with increasing rates of success. These advances are ushering in a new era where crystal structure prediction drives the experimental discovery of new solid forms. After briefly discussing the computational methods that enable successful crystal structure prediction, this perspective presents case studies from the literature that demonstrate how state-of-the-art crystal structure prediction can transform how scientists approach problems involving the organic solid state. Applications to pharmaceuticals, porous organic materials, photomechanical crystals, organic semi-conductors, and nuclear magnetic resonance crystallography are included. Finally, efforts to improve our understanding of which predicted crystal structures can actually be produced experimentally and other outstanding challenges are discussed.

Many-Body Effects in Aqueous Systems: Synergies Between Interaction Analysis Techniques and Force Field Development

Interaction analysis techniques, including the many-body expansion (MBE), symmetry-adapted perturbation theory, and energy decomposition analysis, allow for an intuitive understanding of complex molecular interactions. We review these methods by first providing a historical context for the study of many-body interactions and discussing how nonadditivities emerge from Hamiltonians containing strictly pairwise-additive interactions. We then elaborate on the synergy between these interaction analysis techniques and the development of advanced force fields aimed at accurately reproducing the Born-Oppenheimer potential energy surface. In particular, we focus on ab initio-based force fields that aim to explicitly reproduce many-body terms and are fitted to high-level electronic structure results. These force fields generally incorporate many-body effects through (a) parameterization of distributed multipoles, (b) explicit fitting of the MBE, (c) inclusion of many-atom features in a neural network, and (d) coarse-graining of many-body terms into an effective two-body term. We also discuss the emerging use of the MBE to improve the accuracy and speed of ab initio molecular dynamics.

Three-body interaction effects in heterolytic hydrogen splitting by frustrated Lewis pairs

The reaction of heterolytic dihydrogen splitting by frustrated Lewis pairs P(R)₃ and B(C₆F₅)₃ (where R = t-butyl and 1-adamantene) is driven by strong three-body contributions which originate from the induction and charge transfer effects. The three-body effect increases dramatically as a function of inter-hydrogen distance. As predicted by the symmetry adapted perturbation theory, the "frustration" of Lewis pairs originates from the dual role of the exchange effects. First, the exchange manifests itself in the first-order Pauli repulsion by keeping the pairs away. Second, and equally important, the second-order exchange-induction almost completely cancels the effects of the second-order induction. This suppression of induction effects eases up upon the interaction of the frustrated pairs with H₂. The activation of induction in this instance constitutes the three-body effect.

Assessment of three-body dispersion models against coupled-cluster benchmarks for crystalline benzene, carbon dioxide, and triazine

To study the contribution of three-body dispersion to crystal lattice energies, we compute the three-body contributions to the lattice energies for crystalline benzene, carbon dioxide, and triazine using various computational methods. We show that these contributions converge quickly as the intermolecular distances between the monomers grow. In particular, the smallest value among the three pairwise intermonomer closest-contact distances, R_min, shows a strong correlation with the three-body contribution to the lattice energy, and, here, the largest of the closest-contact distances, R_max, serves as a cutoff criterion to limit the number of trimers to be considered. We considered all trimers up to R_max=15Å. The trimers with R_min<4Å contribute 90.4%, 90.6%, and 93.9% of the total three-body contributions for crystalline benzene, carbon dioxide, and triazine, respectively, for the coupled-cluster singles, doubles, and perturbative triples [CCSD(T)] method. For trimers with R_min>4Å, the second-order Møller-Plesset perturbation theory (MP2) supplemented with the Axilrod-Teller-Muto (ATM) three-body dispersion correction reproduces the CCSD(T) values for the cumulative three-body contributions with errors of less than 0.1 kJ mol^-1. Moreover, three-body contributions are converged within 0.15 kJ mol^-1 by R_max=10Å. From these results, it appears that in molecular crystals where dispersion dominates the three-body contribution to the lattice energy, the trimers with R_min>4Å can be computed with the MP2+ATM method to reduce the computational cost, and those with R_max>10Å appear to be basically negligible.

Development and Implementation of Atomically Anisotropic First-Principles Force Fields: A Benzene Case Study

π-interactions are an important motif in chemical and biochemical systems. However, due to their anisotropic electron densities and complex balance of intermolecular interactions, aromatic molecules represent an ongoing challenge for accurate and transferable force field development. Historically, ab initio force fields for aromatics have not exhibited good accuracy with respect to bulk properties or have only been used to study gas-phase dimers. Using benzene as a proof of concept, herein we show how our own ab initio MASTIFF force field incorporates an atomically anisotropic description of intermolecular interactions to yield an accurate and robust model for aromatic interactions irrespective of phase. Compared to existing models, the MASTIFF benzene force field not only is accurate for liquid phase properties but also offers transferability to the gas and solid phases. Additionally, we introduce a computationally efficient OpenMM plugin which enables customizable anisotropic intermolecular functional forms and which can be generically used in any MD simulation where a model for nonspherical atomic features is required. Overall, our results demonstrate the importance of atomic-level anisotropy in enabling next-generation ab initio force field development.

Physical mechanisms of intermolecular interactions from symmetry-adapted perturbation theory

Symmetry-adapted perturbation theory (SAPT) is a method for computational studies of noncovalent interactions between molecules. This method will be discussed here from the perspective of establishing the paradigm for understanding mechanisms of intermolecular interactions. SAPT interaction energies are obtained as sums of several contributions. Each contribution possesses a clear physical interpretation as it results from some specific physical process. It also exhibits a specific dependence on the intermolecular separation R. The four major contributions are the electrostatic, induction, dispersion, and exchange energies, each due to a different mechanism, valid at any R. In addition, at large R, SAPT interaction energies are seamlessly connected with the corresponding terms in the asymptotic multipole expansion of interaction energy in inverse powers of R. Since such expansion explicitly depends on monomers' multipole moments and polarizabilities, this connection provides additional insights by rigorously relating interaction energies to monomers' properties.

Reliable crystal structure predictions from first principles

An inexpensive and reliable method for molecular crystal structure predictions (CSPs) has been developed. The new CSP protocol starts from a two-dimensional graph of crystal's monomer(s) and utilizes no experimental information. Using results of quantum mechanical calculations for molecular dimers, an accurate two-body, rigid-monomer ab initio-based force field (aiFF) for the crystal is developed. Since CSPs with aiFFs are essentially as expensive as with empirical FFs, tens of thousands of plausible polymorphs generated by the crystal packing procedures can be optimized. Here we show the robustness of this protocol which found the experimental crystal within the 20 most stable predicted polymorphs for each of the 15 investigated molecules. The ranking was further refined by performing periodic density-functional theory (DFT) plus dispersion correction (pDFT+D) calculations for these 20 top-ranked polymorphs, resulting in the experimental crystal ranked as number one for all the systems studied (and the second polymorph, if known, ranked in the top few). Alternatively, the polymorphs generated can be used to improve aiFFs, which also leads to rank one predictions. The proposed CSP protocol should result in aiFFs replacing empirical FFs in CSP research.

FFCASP: A Massively Parallel Crystal Structure Prediction Algorithm

A new algorithm called Fast and Flexible CrystAl Structure Predictor (FFCASP) was developed to predict the structure of covalent and molecular crystals. FFCASP is massively parallel and able to handle more than 200 atoms in the unit cell (in other terms, it allows global optimization around 100 individual parameters). It uses a global optimizer specialized for Crystal Structure Prediction (CSP) which combines particle swarm and simulated annealing optimizers. Three different molecular crystals, including diverse intermolecular interactions, namely, cytosine, coumarin, and pyrazinamide, have been selected to evaluate the performance of FFCASP. While cytosine polymorphs have been searched by employing two different force fields (a DFT-SAPT based intermolecular potential and generalized amber force field (GAFF)) up to Z = 16, only GAFF has been used both in coumarin and pyrazinamide polymorph searches up to Z = 4. For these three molecular crystals, FFCASP generated more than 20 000 crystal structures, and the unique ones have been further treated by DFT-D3. A combination of data mining and a machine learning approach was introduced to determine the unique structures and their distribution into different clusters, which ultimately gives an opportunity to retrieve the common features and relations between the resulting structures. There are two known experimental crystal structures of cytosine, and both were successfully located with FFCASP. Two of the reported crystal structures of coumarin have been reproduced. Similarly, in pyrazinamide, three known experimental structures have been rediscovered. In addition to finding the experimentally known structures, FFCASP also located other low-energy structures for each considered molecular crystals. These successes of FFCASP offer the possibility to discover the polymorphic nature of other important molecular crystals (e.g., drugs) as well.

Data-efficient machine learning for molecular crystal structure prediction

The combination of modern machine learning (ML) approaches with high-quality data from quantum mechanical (QM) calculations can yield models with an unrivalled accuracy/cost ratio. However, such methods are ultimately limited by the computational effort required to produce the reference data. In particular, reference calculations for periodic systems with many atoms can become prohibitively expensive for higher levels of theory. This trade-off is critical in the context of organic crystal structure prediction (CSP). Here, a data-efficient ML approach would be highly desirable, since screening a huge space of possible polymorphs in a narrow energy range requires the assessment of a large number of trial structures with high accuracy. In this contribution, we present tailored Δ-ML models that allow screening a wide range of crystal candidates while adequately describing the subtle interplay between intermolecular interactions such as H-bonding and many-body dispersion effects. This is achieved by enhancing a physics-based description of long-range interactions at the density functional tight binding (DFTB) level-for which an efficient implementation is available-with a short-range ML model trained on high-quality first-principles reference data. The presented workflow is broadly applicable to different molecular materials, without the need for a single periodic calculation at the reference level of theory. We show that this even allows the use of wavefunction methods in CSP.

Structures and Spectroscopic Properties of Large Molecules and Condensed-Phase Systems Predicted by Generalized Energy-Based Fragmentation Approach

ConspectusThe structures and spectroscopic properties of molecules and condensed-phase systems are usually experimentally characterized by X-ray, infrared (IR), Raman, nuclear magnetic resonance (NMR), and electronic absorption/emission spectra. Quantum mechanics (QM) calculations are critical in quantitatively understanding the relationship between the structure and physicochemical properties of various chemical systems. However, it is very challenging to apply traditional QM methods to large molecules and condensed-phase systems with large unit cells due to their steep computational scaling with the system size. To overcome this difficulty, theoretical chemists have developed various linear (or low) scaling QM methods, among which energy-based fragmentation methods have achieved great success for large molecules or clusters. One of the most popular energy-based fragmentation methods is the generalized energy-based fragmentation (GEBF) approach developed by us.In this approach, the ground-state energy of a large molecule can be evaluated from the ground-state energies of a series of embedded subsystems. In this Account, we focus on the recent developments and applicability of the GEBF approach for the structures and spectroscopic properties of complicated large molecules and condensed-phase systems. With new fragmentation schemes, the GEBF approach can now describe ionic liquid clusters and metal-containing supramolecular systems accurately and can provide accurate binding energies for host-guest complexes. In addition, the GEBF approach is now available for describing the localized excited states of large systems including a chromophore. More importantly, the GEBF approach under periodic boundary conditions (PBC-GEBF) has been developed to deal with periodic molecular crystals and liquids. Then, the ground-state energy (or property) per unit cell of a periodic condensed phase system can be predicted with QM calculations on nonperiodic embedded subsystems. This feature enables accurate electron correlation calculations on molecular crystals and liquids to be feasible on ordinary workstations. The PBC-GEBF approach has been applied to predict the crystal structures, lattice energies, and spectroscopic properties of some typical molecular crystals and solutions. By combining the GEBF method and machine learning (ML) method, a GEBF-ML force field has been developed for long normal alkanes, and the IR spectra of long alkanes can be obtained from the GEBF-ML molecular dynamics (MD) simulations. The GEBF and its periodic variant are expected to play increasingly important roles in investigating real-life chemical systems of broad interests at the ab initio levels.

Many-Body Dispersion

A broad range of approaches to many-body dispersion are discussed, including empirical approaches with multiple fitted parameters, augmented density functional-based approaches, symmetry adapted perturbation theory, and a supermolecule approach based on coupled cluster theory. Differing definitions of "body" are considered, specifically atom-based vs molecule-based approaches.

XO-PBC: An Accurate and Efficient Method for Molecular Crystals

In this work, we propose the XO-PBC method, which combines the eXtended ONIOM method (XO) with the periodic boundary condition (PBC) for the description of molecular crystals. XO-PBC tries to embed a finite cluster cut out from the solid into the periodic environment, making it feasible to employ advanced molecular quantum chemistry methods, which are usually prohibitively expensive for direct PBC calculations. In particular, XO-PBC utilizes the results from force calculations to design the scheme to fragment the molecule when crystals are made of large molecules and to select cluster model systems automatically consisting of dimer up to tetramer interactions for embedding. By applying an appropriate theory to each model, a satisfactory accuracy for the system under study is ensured, while a high efficiency is achieved with massively parallel computing by distributing model systems onto different processors. A comparison of the XO-PBC calculations with the conventional direct PBC calculations at the B3LYP level demonstrates its accuracy at substantially low cost for the description of molecular crystals. The usefulness of the XO-PBC method is further exemplified, showing that XO-PBC is able to predict the lattice energies of various types of molecular crystals within chemical accuracy (<4 kJ/mol) when the doubly hybrid density functional XYG3 is used as the target high level and the periodic PBE as the basic low level. The XO-PBC method provides a general protocol that brings the great predictive power of advanced electronic structure methods from molecular systems to the extended solids.

SAPT codes for calculations of intermolecular interaction energies

Symmetry-adapted perturbation theory (SAPT) is a method for calculations of intermolecular (noncovalent) interaction energies. The set of SAPT codes that is described here, the current version named SAPT2020, includes virtually all variants of SAPT developed so far, among them two-body SAPT based on perturbative, coupled cluster, and density functional theory descriptions of monomers, three-body SAPT, and two-body SAPT for some classes of open-shell monomers. The properties of systems governed by noncovalent interactions can be predicted only if potential energy surfaces (force fields) are available. SAPT is the preferred approach for generating such surfaces since it is seamlessly connected to the asymptotic expansion of interaction energy. SAPT2020 includes codes for automatic development of such surfaces, enabling generation of complete dimer surfaces with a rigid monomer approximation for dimers containing about one hundred atoms. These codes can also be used to obtain surfaces including internal degrees of freedom of monomers.

Compressed intramolecular dispersion interactions

The feasibility of the compression of localized virtual orbitals is explored in the context of intramolecular long-range dispersion interactions. Singular value decomposition (SVD) of coupled cluster doubles amplitudes associated with the dispersion interactions is analyzed for a number of long-chain systems, including saturated and unsaturated hydrocarbons and a silane chain. Further decomposition of the most important amplitudes obtained from these SVDs allows for the analysis of the dispersion-specific virtual orbitals that are naturally localized. Consistent with previous work on intermolecular dispersion interactions in dimers, it is found that three important geminals arise and account for the majority of dispersion interactions at the long range, even in the many body intramolecular case. Furthermore, it is shown that as few as three localized virtual orbitals per occupied orbital can be enough to capture all pairwise long-range dispersion interactions within a molecule.

Fantasy versus reality in fragment-based quantum chemistry

Since the introduction of the fragment molecular orbital method 20 years ago, fragment-based approaches have occupied a small but growing niche in quantum chemistry. These methods decompose a large molecular system into subsystems small enough to be amenable to electronic structure calculations, following which the subsystem information is reassembled in order to approximate an otherwise intractable supersystem calculation. Fragmentation sidesteps the steep rise (with respect to system size) in the cost of ab initio calculations, replacing it with a distributed cost across numerous computer processors. Such methods are attractive, in part, because they are easily parallelizable and therefore readily amenable to exascale computing. As such, there has been hope that distributed computing might offer the proverbial "free lunch" in quantum chemistry, with the entrée being high-level calculations on very large systems. While fragment-based quantum chemistry can count many success stories, there also exists a seedy underbelly of rarely acknowledged problems. As these methods begin to mature, it is time to have a serious conversation about what they can and cannot be expected to accomplish in the near future. Both successes and challenges are highlighted in this Perspective.

CrystaLattE: Automated computation of lattice energies of organic crystals exploiting the many-body expansion to achieve dual-level parallelism

We present an algorithm to compute the lattice energies of molecular crystals based on the many-body cluster expansion. The required computations on dimers, trimers, etc., within the crystal are independent of each other, leading to a naturally parallel approach. The algorithm exploits the long-range three-dimensional periodic order of crystals to automatically detect and avoid redundant or unnecessary computations. For this purpose, Coulomb-matrix descriptors from machine learning applications are found to be efficient in determining whether two N-mers are identical. The algorithm is implemented as an open-source Python program, CrystaLattE, that uses some of the features of the Quantum Chemistry Common Driver and Databases library. CrystaLattE is initially interfaced with the quantum chemistry package Psi4. With CrystaLattE, we have applied the fast, dispersion-corrected Hartree-Fock method HF-3c to the lattice energy of crystalline benzene. Including all 73 symmetry-unique dimers and 7130 symmetry-unique trimers that can be formed from molecules within a 15 Å cutoff from a central reference monomer, HF-3c plus an Axilrod-Teller-Muto estimate of three-body dispersion exhibits an error of only -1.0 kJ mol^-1 vs the estimated 0 K experimental lattice energy of -55.3 ± 2.2 kJ mol^-1. The convergence of the HF-3c two- and three-body contributions to the lattice energy as a function of intermonomer distance is examined.

Assessing Electronic Structure Methods for Long-Range Three-Body Dispersion Interactions: Analysis and Calculations on Well-Separated Metal Atom Trimers

Three-body dispersion interactions are much weaker than their two-body counterpart. However, their importance grows quickly as the number of interacting monomers rises. To explore the numerical performance of correlation methods for long-range three-body dispersion, we performed calculations on eight very simple dispersion-dominated model metal trimers: Na₃, Mg₃, Zn₃, Cd₃, Hg₃, Cu₃, Ag₃, and Au₃. One encouraging aspect is that relatively small basis sets of augmented triple-ζ size appear to be adequate for three-body dispersion in the long-range. Coupled cluster calculations were performed at high levels to assess MP3, CCSD, CCSD(T), empirical density functional theory dispersion (D3), and the many-body dispersion (MBD) approach. We found that the accuracy of CCSD(T) was generally significantly lower than for two-body interactions, with errors sometimes reaching 20% in the investigated systems, while CCSD and particularly MP3 were generally more erratic. MBD is found to perform better than D3 at large distances, whereas the opposite is true at shorter distances. When computing reference numbers for three-body dispersion, care should be taken to appropriately represent the effect of the connected triple excitations, which are significant in most cases and incompletely approximated by CCSD(T).

Understanding non-covalent interactions in larger molecular complexes from first principles

Non-covalent interactions pervade all matter and play a fundamental role in layered materials, biological systems, and large molecular complexes. Despite this, our accumulated understanding of non-covalent interactions to date has been mainly developed in the tens-of-atoms molecular regime. This falls considerably short of the scales at which we would like to understand energy trends, structural properties, and temperature dependencies in materials where non-covalent interactions have an appreciable role. However, as more reference information is obtained beyond moderately sized molecular systems, our understanding is improving and we stand to gain pertinent insights by tackling more complex systems, such as supramolecular complexes, molecular crystals, and other soft materials. In addition, accurate reference information is needed to provide the drive for extending the predictive power of more efficient workhorse methods, such as density functional approximations that also approximate van der Waals dispersion interactions. In this perspective, we discuss the first-principles approaches that have been used to obtain reference interaction energies for beyond modestly sized molecular complexes. The methods include quantum Monte Carlo, symmetry-adapted perturbation theory, non-canonical coupled cluster theory, and approaches based on the random-phase approximation. By considering the approximations that underpin each method, the most accurate theoretical references for supramolecular complexes and molecular crystals to date are ascertained. With these, we also assess a handful of widely used exchange-correlation functionals in density functional theory. The discussion culminates in a framework for putting into perspective the accuracy of high-level wavefunction-based methods and identifying future challenges.

Theory and practice of modeling van der Waals interactions in electronic-structure calculations

The accurate description of long-range electron correlation, most prominently including van der Waals (vdW) dispersion interactions, represents a particularly challenging task in the modeling of molecules and materials. vdW forces arise from the interaction of quantum-mechanical fluctuations in the electronic charge density. Within (semi-)local density functional approximations or Hartree-Fock theory such interactions are neglected altogether. Non-covalent vdW interactions, however, are ubiquitous in nature and play a key role for the understanding and accurate description of the stability, dynamics, structure, and response properties in a plethora of systems. During the last decade, many promising methods have been developed for modeling vdW interactions in electronic-structure calculations. These methods include vdW-inclusive Density Functional Theory and correlated post-Hartree-Fock approaches. Here, we focus on the methods within the framework of Density Functional Theory, including non-local van der Waals density functionals, interatomic dispersion models within many-body and pairwise formulation, and random phase approximation-based approaches. This review aims to guide the reader through the theoretical foundations of these methods in a tutorial-style manner and, in particular, highlight practical aspects such as the applicability and the advantages and shortcomings of current vdW-inclusive approaches. In addition, we give an overview of complementary experimental approaches, and discuss tools for the qualitative understanding of non-covalent interactions as well as energy decomposition techniques. Besides representing a reference for the current state-of-the-art, this work is thus also designed as a concise and detailed introduction to vdW-inclusive electronic structure calculations for a general and broad audience.

Nature of halogen bonding involving π-systems, nitroxide radicals and carbenes: a highlight of the importance of charge transfer

The recently developed adiabatic absolutely localized molecular orbital energy decomposition analysis (ALMO-EDA) has proven to be useful in determining the effects of different energy components on the geometries of complexes bound by intermolecular interactions. The authors have applied it to systems such as the water dimer, water-ion complexes, metallocenes and lone-pair type halogen-bonded (XB) dimers. In this study, we have extended the second-generation ALMO-EDA method to 40 different XB complexes by benchmarking against its classical counterpart and symmetry-adapted perturbation theory (SAPT). In addition, we have examined the nature of halogen bonding involving less studied XB acceptors, namely π-systems, radicals and carbenes, using the adiabatic ALMO-EDA analyses, particularly to shed light on how each energy component affects the geometries of the XB complexes. Our results show that the second-generation ALMO-EDA predicts a higher electrostatic energy contribution in all XB complexes compared to SAPT and classical ALMO-EDA schemes. On the other hand, when comparing across different XB acceptors, all three partition schemes produced the same qualitative finding. The adiabatic ALMO-EDA analyses indicate that while the inclusion of a charge transfer contribution is important in achieving accurate XB bond lengths and interaction energies, as well as recovering the binding site specificity of XB involving benzene and naphthalene acceptors, it is sufficient to obtain the linearity of the XB complexes in the frozen approximation.

Fast and accurate quantum Monte Carlo for molecular crystals

Computational approaches based on the fundamental laws of quantum mechanics are now integral to almost all materials design initiatives in academia and industry. If computational materials science is genuinely going to deliver on its promises, then an electronic structure method with consistently high accuracy is urgently needed. We show that, thanks to recent algorithmic advances and the strategy developed in our manuscript, quantum Monte Carlo yields extremely accurate predictions for the lattice energies of materials at a surprisingly modest computational cost. It is thus no longer a technique that requires a world-leading computational facility to obtain meaningful results. While we focus on molecular crystals, the significance of our findings extends to all classes of materials.

Computer simulation plays a central role in modern-day materials science. The utility of a given computational approach depends largely on the balance it provides between accuracy and computational cost. Molecular crystals are a class of materials of great technological importance which are challenging for even the most sophisticated ab initio electronic structure theories to accurately describe. This is partly because they are held together by a balance of weak intermolecular forces but also because the primitive cells of molecular crystals are often substantially larger than those of atomic solids. Here, we demonstrate that diffusion quantum Monte Carlo (DMC) delivers subchemical accuracy for a diverse set of molecular crystals at a surprisingly moderate computational cost. As such, we anticipate that DMC can play an important role in understanding and predicting the properties of a large number of molecular crystals, including those built from relatively large molecules which are far beyond reach of other high-accuracy methods.

The nature of three-body interactions in DFT: Exchange and polarization effects

We propose a physically motivated decomposition of density functional theory (DFT) 3-body nonadditive interaction energies into the exchange and density-deformation (polarization) components. The exchange component represents the effect of the Pauli exclusion in the wave function of the trimer and is found to be challenging for density functional approximations (DFAs). The remaining density-deformation nonadditivity is less dependent upon the DFAs. Numerical demonstration is carried out for rare gas atom trimers, Ar₂-HX (X = F, Cl) complexes, and small hydrogen-bonded and van der Waals molecular systems. None of the tested semilocal, hybrid, and range-separated DFAs properly accounts for the nonadditive exchange in dispersion-bonded trimers. By contrast, for hydrogen-bonded systems, range-separated DFAs achieve a qualitative agreement to within 20% of the reference exchange energy. A reliable performance for all systems is obtained only when the monomers interact through the Hartree-Fock potential in the dispersion-free Pauli blockade scheme. Additionally, we identify the nonadditive second-order exchange-dispersion energy as an important but overlooked contribution in force-field-like dispersion corrections. Our results suggest that range-separated functionals do not include this component, although semilocal and global hybrid DFAs appear to imitate it in the short range.

First-Principles Models for van der Waals Interactions in Molecules and Materials: Concepts, Theory, and Applications

Noncovalent van der Waals (vdW) or dispersion forces are ubiquitous in nature and influence the structure, stability, dynamics, and function of molecules and materials throughout chemistry, biology, physics, and materials science. These forces are quantum mechanical in origin and arise from electrostatic interactions between fluctuations in the electronic charge density. Here, we explore the conceptual and mathematical ingredients required for an exact treatment of vdW interactions, and present a systematic and unified framework for classifying the current first-principles vdW methods based on the adiabatic-connection fluctuation-dissipation (ACFD) theorem (namely the Rutgers-Chalmers vdW-DF, Vydrov-Van Voorhis (VV), exchange-hole dipole moment (XDM), Tkatchenko-Scheffler (TS), many-body dispersion (MBD), and random-phase approximation (RPA) approaches). Particular attention is paid to the intriguing nature of many-body vdW interactions, whose fundamental relevance has recently been highlighted in several landmark experiments. The performance of these models in predicting binding energetics as well as structural, electronic, and thermodynamic properties is connected with the theoretical concepts and provides a numerical summary of the state-of-the-art in the field. We conclude with a roadmap of the conceptual, methodological, practical, and numerical challenges that remain in obtaining a universally applicable and truly predictive vdW method for realistic molecular systems and materials.

Exploring polymorphism of benzene and naphthalene with free energy based enhanced molecular dynamics

Prediction and exploration of possible polymorphism in organic crystal compounds are of great importance for industries ranging from organic electronics to pharmaceuticals to high-energy materials. Here we apply our crystal structure prediction procedure and the enhanced molecular dynamics based sampling approach called the Crystal-Adiabatic Free Energy Dynamics (Crystal-AFED) method to benzene and naphthalene. Crystal-AFED allows the free energy landscape of structures to be explored efficiently at any desired temperature and pressure. For each system, we successfully predict the most stable crystal structures at atmospheric pressure and explore the relative Gibbs free energies of predicted polymorphs at high pressures. Using Crystal-AFED sampling, we find that mixed structures, which typically cannot be discovered by standard crystal structure prediction methods, are prevalent in the solid forms of these compounds at high pressure.

Small Atomic Orbital Basis Set First-Principles Quantum Chemical Methods for Large Molecular and Periodic Systems: A Critical Analysis of Error Sources

In quantum chemical computations the combination of Hartree-Fock or a density functional theory (DFT) approximation with relatively small atomic orbital basis sets of double-zeta quality is still widely used, for example, in the popular B3LYP/6-31G* approach. In this Review, we critically analyze the two main sources of error in such computations, that is, the basis set superposition error on the one hand and the missing London dispersion interactions on the other. We review various strategies to correct those errors and present exemplary calculations on mainly noncovalently bound systems of widely varying size. Energies and geometries of small dimers, large supramolecular complexes, and molecular crystals are covered. We conclude that it is not justified to rely on fortunate error compensation, as the main inconsistencies can be cured by modern correction schemes which clearly outperform the plain mean-field methods.

Modeling Polymorphic Molecular Crystals with Electronic Structure Theory

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.

Next-Generation Force Fields from Symmetry-Adapted Perturbation Theory

Symmetry-adapted perturbation theory (SAPT) provides a unique set of advantages for parameterizing next-generation force fields from first principles. SAPT provides a direct, basis-set superposition error free estimate of molecular interaction energies, a physically intuitive energy decomposition, and a seamless transition to an asymptotic picture of intermolecular interactions. These properties have been exploited throughout the literature to develop next-generation force fields for a variety of applications, including classical molecular dynamics simulations, crystal structure prediction, and quantum dynamics/spectroscopy. This review provides a brief overview of the formalism and theory of SAPT, along with a practical discussion of the various methodologies utilized to parameterize force fields from SAPT calculations. It also highlights a number of applications of SAPT-based force fields for chemical systems of particular interest. Finally, the review ends with a brief outlook on the future opportunities and challenges that remain for next-generation force fields based on SAPT.

Benchmark Calculations of Three-Body Intermolecular Interactions and the Performance of Low-Cost Electronic Structure Methods

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.

Free-energy calculations using classical molecular simulation: application to the determination of the melting point and chemical potential of a flexible RDX model

Several methods are used in sequence to determine the chemical potential of atomistic RDX in the solid and liquid phases, and its corresponding melting point. Results yield the thermodynamic melting point of 488.75 K at 1.0 atm.

Fragment-Based Electronic Structure Approach for Computing Nuclear Magnetic Resonance Chemical Shifts in Molecular Crystals

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 (13)C spectra is assessed. This challenging example requires discriminating among spectra whose (13)C 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 (13)C 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.