Ab initio molecular crystal structures, spectra, and phase diagrams

First-Principles Models of Polymorphism of Pharmaceuticals: Maximizing the Accuracy-to-Cost Ratio

Accuracy and sophistication of in silico models of structure, internal dynamics, and cohesion of molecular materials at finite temperatures increase over time. Applicability limits of ab initio polymorph ranking that would be feasible at reasonable costs currently represent crystals of moderately sized molecules (less than 20 nonhydrogen atoms) and simple unit cells (containing rather only one symmetry-irreducible molecule). Extending the applicability range of the underlying first-principles methods to larger systems with a real-life significance, and enabling to perform such computations in a high-throughput regime represent additional challenges to be tackled in computational chemistry. This work presents a novel composite method that combines the computational efficiency of density-functional tight-binding (DFTB) methods with the accuracy of density-functional theory (DFT). Being rooted in the quasi-harmonic approximation, it uses a cheap method to perform all of the costly scans of how static and dynamic characteristics of the crystal vary with respect to its volume. Such data are subsequently corrected to agree with a higher-level model, which must be evaluated only at a single volume of the crystal. It thus enables predictions of structural, cohesive, and thermodynamic properties of complex molecular materials, such as pharmaceuticals or organic semiconductors, at a fraction of the original computational cost. As the composite model retains the solid physical background, it suffers from a minimum accuracy deterioration compared to the full treatment with the costly approach. The novel methodology is demonstrated to provide consistent results for the structural and thermodynamic properties of real-life molecular crystals and their polymorph ranking.

Improved Description of Intra- and Intermolecular Interactions through Dispersion-Corrected Second-Order Møller-Plesset Perturbation Theory

ConspectusThe quantum chemical modeling of organic crystals and other molecular condensed-phase problems requires computationally affordable electronic structure methods which can simultaneously describe intramolecular conformational energies and intermolecular interactions accurately. To achieve this, we have developed a spin-component-scaled, dispersion-corrected second-order Møller-Plesset perturbation theory (SCS-MP2D) model. SCS-MP2D augments canonical MP2 with a dispersion correction which removes the uncoupled Hartree-Fock dispersion energy present in canonical MP2 and replaces it with a more reliable coupled Kohn-Sham treatment, all evaluated within the framework of Grimme's D3 dispersion model. The spin-component scaling is then used to improve the description of the residual (nondispersion) portion of the correlation energy.The SCS-MP2D model improves upon earlier corrected MP2 models in a few ways. Compared to the highly successful dispersion-corrected MP2C model, which is based solely on intermolecular perturbation theory, the SCS-MP2D dispersion correction improves the description of both inter- and intramolecular interactions. The dispersion correction can also be evaluated with trivial computational cost, and nuclear analytic gradients are computed readily to enable geometry optimizations. In contrast to earlier spin-component scaling MP2 models, the optimal spin-component scaling coefficients are only mildly sensitive to the choice of training data, and a single global parametrization of the model can describe both thermochemistry and noncovalent interactions.The resulting dispersion-corrected, spin-component-scaled MP2 (SCS-MP2D) model predicts conformational energies and intermolecular interactions with accuracy comparable to or better than those of many range-separated and double-hybrid density functionals, as is demonstrated on a variety of benchmark tests. Among the functionals considered here, only the revDSD-PBEP86-D3(BJ) functional gives consistently smaller errors in benchmark tests. The results presented also hint that further improvements of SCS-MP2D may be possible through a more robust fitting procedure for the seven empirical parameters.To demonstrate the performance of SCS-MP2D further, several applications to molecular crystal problems are presented. The three chosen examples all represent cases where density-driven delocalization error causes GGA or hybrid density functionals to artificially stabilize crystals exhibiting more extended π-conjugation. Our pragmatic strategy addresses the delocalization error by combining a periodic density functional theory (DFT) treatment of the infinite lattice with intramolecular/conformational energy corrections computed with SCS-MP2D. For the anticancer drug axitinib, applying the SCS-MP2D conformational energy correction produces crystal polymorph stabilities that are consistent with experiment, in contrast to earlier studies. For the crystal structure prediction of the ROY molecule, so named for its colorful red, orange, and yellow crystals, this approach leads to the first plausible crystal energy landscape, and it reveals that the lowest-energy polymorphs have already been found experimentally. Finally, in the context of photomechanical crystals, which transform light into mechanical work, these techniques are used to predict the structural transformations and extract design principles for maximizing the work performed.

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.

Competition between Solvent···Solvent and Solvent···Solute Interactions in the Microhydration of the Tetrafluoroborate Anion, BF4-(H2O)n=1,2,3,4

This study systematically examines the interactions of the tetrafluoroborate anion (BF4-) with up to four water molecules (BF4-(H2O)n=1,2,3,4). Full geometry optimizations and subsequent harmonic vibrational frequency computations are performed using a variety of density functional theory (DFT) methods (B3LYP, B3LYP-D3BJ, and M06-2X) and the MP2 ab initio method with a triple-ζ correlation consistent basis set augmented with diffuse functions on all non-hydrogen atoms (cc-pVTZ for H and aug-cc-pVTZ for B, O, and F; denoted as haTZ). Optimized structures and harmonic vibrational frequencies were also obtained with the CCSD(T) ab initio method and the haTZ basis set for the mono- and dihydrate (n = 1, 2) structures. The 2-body:Many-body (2b:Mb) technique, in which CCSD(T) computations capture the 1- and 2-body contributions to the interactions and MP2 computations recover all higher-order contributions, was used to extend these demanding computations to the tri- and tetrahydrate (n = 3, 4) systems. Four, five, and eight new stationary points have been identified for the di-, tri-, and tetrahydrate systems, respectively. The global minimum of the monohydrate adopts a symmetric double ionic hydrogen bond motif with C2v symmetry and an electronic dissociation energy of 13.17 kcal mol-1 at the CCSD(T)/haTZ level of theory. This strong solvent···solute interaction, however, competes with solute···solute interactions in the lowest-energy BF4-(H2O)n=2,3,4 minima that are not seen in the other di-, tri-, or tetrahydrate minima. The latter interactions help increase the 2b:Mb dissociation energies to more than 26, 41, and 51 kcal mol-1 for n = 2, 3, and 4, respectively. Structures that form hydrogen bonds between the solvating water molecules also exhibit the largest shifts in the harmonic OH stretching frequencies for the waters of hydration. These shifts can exceed -280 cm-1 relative to an isolated H2O molecule at the 2b:Mb/haTZ level of theory.

Discriminating and understanding molecular crystal polymorphism

Polymorph discrimination for a molecular crystal has long been a challenging task, which, nonetheless, is a major concern in the pharmaceutical industry. In this work, we have investigated polymorph discrimination on three different molecular crystals, tetrolic acid, oxalic acid, and oxalyl dihydrazide, covering both packing polymorphism and conformational polymorphism. To gain more understanding, we have performed energy decomposition analysis based on many-body expansion, and have compared the results from the XO-PBC method, that is, the eXtended ONIOM method (XO) with the periodic boundary condition (PBC), with those from some commonly used dispersion corrected density functional theory (DFT-D) methods. It is shown here that, with the XYG3 doubly hybrid functional chosen as the target high level to capture the intra- and short-range intermolecular interactions, and the periodic PBE as the basic low level to take long range interactions into account, the XO-PBC(XYG3:PBE) method not only obtains the correct experimental stability orderings, but also predicts reasonable polymorph energy ranges for all three cases. Our results have demonstrated the usefulness of the present theoretical methods, in particular XO-PBC, while highlighted the importance of a better treatment of different kinds of interactions to be beneficial to polymorph control.

Benchmarking two-body contributions to crystal lattice energies and a range-dependent assessment of approximate methods

Using the many-body expansion to predict crystal lattice energies (CLEs), a pleasantly parallel process, allows for flexibility in the choice of theoretical methods. Benchmark-level two-body contributions to CLEs of 23 molecular crystals have been computed using interaction energies of dimers with minimum inter-monomer separations (i.e., closest contact distances) up to 30 Å. In a search for ways to reduce the computational expense of calculating accurate CLEs, we have computed these two-body contributions with 15 different quantum chemical levels of theory and compared these energies to those computed with coupled-cluster in the complete basis set (CBS) limit. Interaction energies of the more distant dimers are easier to compute accurately and several of the methods tested are suitable as replacements for coupled-cluster through perturbative triples for all but the closest dimers. For our dataset, sub-kJ mol^-1 accuracy can be obtained when calculating two-body interaction energies of dimers with separations shorter than 4 Å with coupled-cluster with single, double, and perturbative triple excitations/CBS and dimers with separations longer than 4 Å with MP2.5/aug-cc-pVDZ, among other schemes, reducing the number of dimers to be computed with coupled-cluster by as much as 98%.

Predicting the structures and vibrational spectra of molecular crystals containing large molecules with the generalized energy-based fragmentation approach

The generalized energy-based fragmentation (GEBF) approach under periodic boundary conditions (PBCs) has been developed to facilitate calculations of molecular crystals containing large molecules. The PBC-GEBF approach can help predict structures and properties of molecular crystals at different theory levels by performing molecular quantum chemistry calculations on a series of non-periodic subsystems constructed from the studied systems. A more rigorous formula of the forces on translational vectors of molecular crystals was proposed and implemented, enabling more reliable predictions of crystal structures. Our benchmark results on several typical molecular crystals show that the PBC-GEBF approach could reproduce the forces on atoms and the translational vectors and the optimized crystal structures from the corresponding conventional periodic methods. The improved PBC-GEBF approach is then applied to predict the crystal structures and vibrational spectra of two molecular crystals containing large molecules. The PBC-GEBF approach can provide a satisfactory description on the crystal structure of a molecular crystal containing 312 atoms in a unit cell at density-fitting second-order Møller-Plesset perturbation theory and density functional theory (DFT) levels and the infrared vibrational spectra of another molecular crystal containing 864 atoms in a unit cell at the DFT level. The PBC-GEBF approach is expected to be a promising theoretical tool for electronic structure calculations on molecular crystals containing large molecules.

Range-dependence of two-body intermolecular interactions and their energy components in molecular crystals

Routinely assessing the stability of molecular crystals with high accuracy remains an open challenge in the computational sciences. The many-body expansion decomposes computation of the crystal lattice energy into an embarrassingly parallel collection of computations over molecular dimers, trimers, and so forth, making quantum chemistry techniques tractable for many crystals of small organic molecules. By examining the range-dependence of different types of energetic contributions to the crystal lattice energy, we can glean qualitative understanding of solid-state intermolecular interactions as well as practical, exploitable reductions in the number of computations required for accurate energies. Here, we assess the range-dependent character of two-body interactions of 24 small organic molecular crystals using the physically interpretable components from symmetry-adapted perturbation theory (electrostatics, exchange repulsion, induction/polarization, and London dispersion). We also examine correlations between the convergence rates of electrostatics and London dispersion terms with molecular dipole moments and polarizabilities, to provide guidance for estimating convergence rates in other molecular crystals.

Structures and properties of ionic crystals and condensed phase ionic liquids predicted with the generalized energy-based fragmentation method

The generalized energy-based fragmentation (GEBF) approach is extended to facilitate ab initio investigations of structures, lattice energies, vibrational spectra and ¹ H NMR chemical shifts of ionic crystals and condensed-phase ionic liquids (ILs) with the periodic boundary conditions (PBC). For selected periodic systems, our results demonstrate that the so-called PBC-GEBF approach can provide satisfactory descriptions on ground-state energies, structures, and vibrational spectra of ionic crystals and IL crystals. The PBC-GEBF approach is then applied to three realistic condensed phase systems. For three ionic crystals (LiCl, NaCl, and KCl), we apply the PBC-GEBF approach with MP2 theory as well as some popular DFT methods to investigate their crystal structures and lattice energies. Our calculations indicate that the crystal structures obtained with PBC-GEBF-MP2/6-311 + G** are very close to the corresponding X-ray structures, while PBC-GEBF-ωB97X-D/6-311 + G** provides satisfactory prediction for crystal structures and lattice energies. For two polymorphs of [n-C₄ mim][Cl] crystals, we find that the PBC-GEBF approach at the M06-2X/6-311 + G** level can give a satisfactory descriptions on structures and Raman spectra of these two crystals. Furthermore, for [C₂ mim][BF₄ ] ILs, we demonstrate that their ¹ H NMR chemical shifts can be estimated from averaging over 5 typical snapshots (extracted from MD simulations) with the PBC-GEBF approach at the B97-2/pcSseg-2 level. The calculated results account for the observed experimental data quite well. Therefore, we expect that the PBC-GEBF approach, combined with various quantum chemistry methods, will become an effective tool in predicting structures and properties of ionic crystals and condensed-phase ILs.

In Situ Determination of Speciation and Local Structure of NaCl-SrCl2 and LiF-ZrF4 Molten Salts

Understanding the local environment of the metal atoms in salt melts is important for modeling the properties of melts and predicting their behavior and thus helping enable the development of technologies such as molten salt reactors and solar-thermal power systems and new approaches to recycling rare-earth metals. Toward that end, we have developed an in situ approach for measuring the coordination of metals in molten salt coupling X-ray absorption spectroscopy (XAS) and Raman spectroscopy. Our approach was demonstrated for two salt mixtures (1.9 and 5 mol % SrCl₂ in NaCl, 0.8 and 5 mol % ZrF₄ in LiF) at up to 1100 °C. Near-edge (X-ray absorption near-edge structure, XANES) and extended X-ray absorption fine structure (EXAFS) spectra were measured. The EXAFS response was modeled using ab initio FEFF calculations. Strontium's first shell is observed to be coordinated with chlorine (Sr²⁺-Cl^-) and zirconium's first shell is coordinated by fluorine (Zr⁴⁺-F^-), both having coordination numbers that decrease with increasing temperature. Multiple zirconium complexes are believed to be present in the melt, which may interfere and distort the EXAFS spectra and result in an anomalously low zirconium first shell coordination number. The use of boron nitride (BN) powder as a salt diluent for XAFS measurements was found to not interfere with measurements and thus can be used for investigations of such systems.

Third-Order Many-Body Expansion of OSV-MP2 Wave Function for Low-Order Scaling Analytical Gradient Computation

We present a many-body expansion (MBE) formulation and implementation for efficient computation of analytical energy gradients from the orbital-specific-virtual second-order Møllet-Plesset perturbation theory (OSV-MP2) based on our earlier work (Zhou et al. J. Chem. Theory Comput. 2020, 16, 196-210). The third-order MBE(3) expansion of OSV-MP2 amplitudes and density matrices was developed to adopt the orbital-specific clustering and long-range termination schemes, which avoids term-by-term differentiations of the MBE energy bodies. We achieve better efficiency by exploiting the algorithmic sparsity that allows us to prune out insignificant fitting integrals and OSV relaxations. With these approximations, the present implementation is benchmarked on a range of molecules that show an economic scaling in the linear and quadratic regimes for computing MBE(3)-OSV-MP2 amplitude and gradient equations, respectively, and yields normal accuracy comparable to the original OSV-MP2 results. The MPI-3-based parallelism through shared memory one-sided communication is further developed for improving parallel scalability and memory accessibility by sorting the MBE(3) orbital clusters into independent tasks that are distributed on multiple processes across many nodes, supporting both global and local data locations in which selected MBE(3)-OSV-MP2 intermediates of different sizes are distinguished and accordingly placed. The accuracy and efficiency level of our MBE(3)-OSV-MP2 analytical gradient implementation is finally illustrated in two applications: we show that the subtle coordination structure differences of mechanically interlocked Cu-catenane complexes can be distinguished when tuning ligand lengths; and the porphycene molecular dynamics reveals the emergence of the vibrational signature arising from softened N-H stretching associated with hydrogen transfer, using an MP2 level of electron correlation and classical nuclei for the first time.

Monte Carlo MP2-F12 for Noncovalent Interactions: The C60 Dimer

A scalable stochastic algorithm is presented that can evaluate explicitly correlated (F12) second-order many-body perturbation (MP2) energies of weak, noncovalent, intermolecular interactions. It first transforms the formulas of the MP2 and F12 energy differences into a short sum of high-dimensional integrals of Green's functions in real space and imaginary time. These integrals are then evaluated by the Monte Carlo method augmented by parallel execution, redundant-walker convergence acceleration, direct-sampling autocorrelation elimination, and control-variate error reduction. By sharing electron-pair walkers across the supermolecule and its subsystems spanned by the joint basis set, the statistical uncertainty is reduced by one to 2 orders of magnitude in the MP2 binding energy corrected for the basis-set incompleteness and superposition errors. The method predicts the MP2-F12/aug-cc-pVDZ binding energy of 19.1 ± 4.0 kcal mol^-1 for the C₆₀ dimer at the center distance of 9.748 Å.

Advances, challenges and perspectives of quantum chemical approaches in molecular spectroscopy of the condensed phase

The purpose of this review is to demonstrate advances, challenges and perspectives of quantum chemical approaches in molecular spectroscopy of the condensed phase. Molecular spectroscopy, particularly vibrational spectroscopy and electronic spectroscopy, has been used extensively for a wide range of areas of chemical sciences and materials science as well as nano- and biosciences because it provides valuable information about structure, functions, and reactions of molecules. In the meantime, quantum chemical approaches play crucial roles in the spectral analysis. They also yield important knowledge about molecular and electronic structures as well as electronic transitions. The combination of spectroscopic approaches and quantum chemical calculations is a powerful tool for science, in general. Thus, our article, which treats various spectroscopy and quantum chemical approaches, should have strong implications in the wider scientific community. This review covers a wide area of molecular spectroscopy from far-ultraviolet (FUV, 120-200 nm) to far-infrared (FIR, 400-10 cm-1)/terahertz and Raman spectroscopy. As quantum chemical approaches, we introduce several anharmonic approaches such as vibrational self-consistent field (VSCF) and the combination of periodic harmonic calculations with anharmonic corrections based on finite models, grid-based techniques like the Numerov approach, the Cartesian coordinate tensor transfer (CCT) method, Symmetry-Adapted Cluster Configuration-Interaction (SAC-CI), and the ZINDO (Semi-empirical calculations at Zerner's Intermediate Neglect of Differential Overlap). One can use anharmonic approaches and grid-based approaches for both infrared (IR) and near-infrared (NIR) spectroscopy, while CCT methods are employed for Raman, Raman optical activity (ROA), FIR/terahertz and low-frequency Raman spectroscopy. Therefore, this review overviews cross relations between molecular spectroscopy and quantum chemical approaches, and provides various kinds of close-reality advanced spectral simulation for condensed phases.

Ab initio determination of crystal stability of di-p-tolyl disulfide

With the rapid growth of energy demand and the depletion of existing energy resources, the new materials with superior performances, low costs and environmental friendliness for energy production and storage are explored. Di-p-tolyl disulfide (p-Tol₂S₂) is a typical lubricating material, which has been applied in the field of energy storage. The conformational properties and phase transformations of p-Tol₂S₂ have been studied by pioneers, but their polymorphs and the polymorphism induced crystal structure changes require further analysis. In this study, we perform the crystal structural screening, prediction and optimization of p-Tol₂S₂ crystal with quantum mechanical calculations, i.e., density functional theory (DFT) and second-order Møller–Plesset perturbation (MP2) methods. A series of crystal structures with different molecular arrangements are generated based on the crystal structure screening. As compared to long-established lattice energy calculation, we take an advantage of using more accurate technique, which is Gibbs free energy calculation. It considers the effects of entropy and temperature to predict the crystal structures and energy landscape. By comparing the Gibbs free energies between predicted and experimental structures, we found that phase α is the most stable structure for p-Tol₂S₂ crystal at ambient temperature and standard atmospheric pressure. Furthermore, we provide an efficient method to discriminate different polymorphs that are otherwise difficult to be identified based on the Raman/IR spectra. The proposed work enable us to evaluate the quality of various crystal polymorphs rapidly.

Sublimation Properties of α,ω-Diamines Revisited from First-Principles Calculations

Sublimation enthalpies of alkane-α,ω-diamines exhibit an odd-even pattern within their homologous series. First-principles calculations coupled with the quasi-harmonic approximation for crystals and with the conformation mixing model for the ideal gas are used to explain this phenomenon from the theoretical point of view. Crystals of the odd and even alkane-α,ω-diamines distinctly differ in their packing motifs. However, first-principles calculations indicate that it is a delicate interplay of the cohesive forces, phonons, molecular vibrations and conformational equilibrium which governs the odd-even pattern of the sublimation enthalpies within the homologous series. High molecular flexibility of the alkane-α,ω-diamines predetermines higher sensitivity of the computational model to the quality of the optimized geometries and relative conformational energies. Performance of high-throughput computational methods, such as the density functional tight binding (DFTB, GFN2-xTB) and the explicitly correlated dispersion-corrected Møller-Plesset perturbative method (MP2C-F12), are benchmarked against the consistent state-of-the-art calculations of conformational energies and interaction energies, respectively.

Fragment-based quantum mechanical approach to biomolecules, molecular clusters, molecular crystals and liquids

To study large molecular systems beyond the system size that the current state-of-the-art ab initio electronic structure methods could handle, fragment-based quantum mechanical (QM) approaches have been developed over the past years, and proved to be efficient in dealing with large molecular systems at various ab initio levels. According to the fragmentation approach, a large molecular system can be divided into subsystems (fragments), and subsequently the property of the whole system can be approximately obtained by taking a proper combination of the corresponding terms of individual fragments. Therefore, the standard QM calculation of a large system could be circumvented by carrying out a series of calculations on small fragments, which significantly promotes computational efficiency. The electrostatically embedded generalized molecular fractionation with conjugate caps (EE-GMFCC) method is one of the fragment-based QM approaches which has been developed by our research group in recent years. This Perspective presents the theoretical framework of this fragmentation method and its applications in biomolecules, molecular clusters, molecular crystals and liquids, including total energy calculation, protein-ligand/protein binding affinity prediction, geometry optimization, vibrational spectrum simulation, ab initio molecular dynamics simulation, and prediction of excited-state properties.

Ab Initio Prediction of the Phase Transition for Solid Ammonia at High Pressures

Ammonia is one of the most basic components on the planet and its high-pressure characteristics play an important role in planetary science. Solid ammonia crystals frequently adopt multiple distinct polymorphs exhibiting different properties. Predicting the crystal structure of these polymorphs and under what thermodynamic conditions these polymorphs are stable would be of great value to environmental industry and other fields. Theoretical calculations based on the classical force fields and density-functional theory (DFT) are versatile methods but lack of accurate description of weak intermolecular interactions for molecular crystals. In this study, we employ an ab initio computational study on the solid ammonia at high pressures, using the second-order Møller-Plesset perturbation (MP2) theory and the coupled cluster singles, doubles, and perturbative triples (CCSD(T)) theory along with the embedded fragmentation method. The proposed algorithm is capable of performing large-scale calculations using high-level wavefunction theories, and accurately describing covalent, ionic, hydrogen bonding, and dispersion interactions within molecular crystals, and therefore can predict the crystal structures, Raman spectra and phase transition of solid ammonia phases I and IV accurately. We confirm the crystal structures of solid ammonia phases I and IV that have been controversial for a long time and predict their phase transition that occurs at 1.17 GPa and 210 K with small temperature dependence, which is in line with experiment.

Complexity Reduction in Density Functional Theory Calculations of Large Systems: System Partitioning and Fragment Embedding

With the development of low order scaling methods for performing Kohn-Sham density functional theory, it is now possible to perform fully quantum mechanical calculations of systems containing tens of thousands of atoms. However, with an increase in the size of the system treated comes an increase in complexity, making it challenging to analyze such large systems and determine the cause of emergent properties. To address this issue, in this paper, we present a systematic complexity reduction methodology which can break down large systems into their constituent fragments and quantify interfragment interactions. The methodology proposed here requires no a priori information or user interaction, allowing a single workflow to be automatically applied to any system of interest. We apply this approach to a variety of different systems and show how it allows for the derivation of new system descriptors, the design of QM/MM partitioning schemes, and the novel application of graph metrics to molecules and materials.

The many-body expansion approach to ab initio calculation of electric field gradients in molecular crystals

Accurate calculation of electric field gradients (EFGs) in molecular crystals, despite big advances in ab initio techniques, is still a challenge. Here, we present a new approach to calculate the EFGs in molecular crystals by employing the many-body expansion (MBE) technique with electrostatic embedding. This allows for (i) a reduction in the computational cost or an alternative increase in the level of theory (we use the MP2/6-311++G) and (ii) the ability to monitor EFG convergence by progressively adding more surrounding molecules and/or adding higher many-body interactions. We focus on the ¹⁴N EFG and study four (model) compounds in more detail: solid nitrogen, ethylamine, methylamine, and ammonia. Solid nitrogen is rather insensitive to neighbors; for ethylamine and methylamine, the 3-body interactions are found sufficient for a converged EFG, whereas for ammonia, even the inclusion of 5-body interactions is insufficient although convergence is anticipated. We then validate our technique by comparing the experimental and ab initio¹⁴N EFGs for 116 organic compounds utilizing their known crystal structures and published EFG. Overall, we find a very good agreement, with a small EFG rms error, which is probably due to other sources, rather than the MBE approximation.

Overcoming the difficulties of predicting conformational polymorph energetics in molecular crystals via correlated wavefunction methods

Molecular crystal structure prediction is increasingly being applied to study the solid form landscapes of larger, more flexible pharmaceutical molecules. Despite many successes in crystal structure prediction, van der Waals-inclusive density functional theory (DFT) methods exhibit serious failures predicting the polymorph stabilities for a number of systems exhibiting conformational polymorphism, where changes in intramolecular conformation lead to different intermolecular crystal packings. Here, the stabilities of the conformational polymorphs of o-acetamidobenzamide, ROY, and oxalyl dihydrazide are examined in detail. DFT functionals that have previously been very successful in crystal structure prediction perform poorly in all three systems, due primarily to the poor intramolecular conformational energies, but also due to the intermolecular description in oxalyl dihydrazide. In all three cases, a fragment-based dispersion-corrected second-order Møller-Plesset perturbation theory (MP2D) treatment of the crystals overcomes these difficulties and predicts conformational polymorph stabilities in good agreement with experiment. These results highlight the need for methods which go beyond current-generation DFT functionals to make crystal polymorph stability predictions truly reliable.

Phase Transition of Ice at High Pressures and Low Temperatures

The behavior of ice under extreme conditions undergoes the change of intermolecular binding patterns and leads to the structural phase transitions, which are needed for modeling the convection and internal structure of the giant planets and moons of the solar system as well as H2O-rich exoplanets. Such extreme conditions limit the structural explorations in laboratory but open a door for the theoretical study. The ice phases IX and XIII are located in the high pressure and low temperature region of the phase diagram. However, to the best of our knowledge, the phase transition boundary between these two phases is still not clear. In this work, based on the second-order Møller–Plesset perturbation (MP2) theory, we theoretically investigate the ice phases IX and XIII and predict their structures, vibrational spectra and Gibbs free energies at various extreme conditions, and for the first time confirm that the phase transition from ice IX to XIII can occur around 0.30 GPa and 154 K. The proposed work, taking into account the many-body electrostatic effect and the dispersion interactions from the first principles, opens up the possibility of completing the ice phase diagram and provides an efficient method to explore new phases of molecular crystals.

Ab initio-enabled phase transition prediction of solid carbon dioxide at ultra-high temperatures

Carbon dioxide is one of the fundamental chemical species on Earth, but its solid-phase behavior at high pressures is still far from well understood and some phases remain uncertain or unknown, which increases the challenge to predict its structures.

Prediction of properties from first principles with quantitative accuracy: six representative ice phases

Based on a high-level MP2 theory with the fragment approach, the crystal structure, vibration spectra and phase transitions of six representative ice phases (II, VI, VII, VIII, IX, and XV) are predicted.

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.

Cohesive properties of the crystalline phases of twenty proteinogenic α-aminoacids from first-principles calculations

Cohesive properties (lattice and cohesive energy of the crystal and corresponding sublimation enthalpy) of the complete set of twenty enantiopure anhydrous proteinogenic amino acids are investigated using first-principles calculations. In contrast to neutral amino acid molecules in the vapor phase, all amino acids form crystals in their zwitterionic form. Therefore, reliable ab initio calculations of the proton transfer energy are an indispensable step of such calculations. Simplifying procedures, designed to rationalize the computational cost of the quasi-harmonic approximation, which proves too demanding if performed fully at the given quantum level of theory, are presented and tested. For this purpose, atomic multipoles (up to the quadrupoles) for the amoeba force field are parametrized for all amino acid zwitterions. While the calculated lattice energies of the amino acids range from 235-458 kJ mol-1 in absolute value, the proton transfer energies typically amount to 100-220 kJ mol-1, which translates to sublimation enthalpies ranging from 117-202 kJ mol-1, appreciably exceeding the sublimation enthalpy values common for nonionic molecular crystals. Critically assessed experimental data on sublimation enthalpies are used as a benchmark for comparison of the data calculated in this work. Cohesive properties of most amino acids calculated in this work, combining the PBE-D3(BJ)/PAW and CCSD(T)-F12/aug-cc-pVDZ levels of theory used for predictions of the lattice energies and of the proton transfer energies, respectively, exhibit a reasonable agreement with the experiment. At the same time, this work contains the first published data on cohesive properties for several enantiopure amino acids.

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.

Low-Temperature Polymorphic Transformation of β-Lactam Antibiotics

Abstract: Polymorphic screening and transformation of molecular crystals are presently popular research areas in pharmaceutical studies. In this study, we developed an ab initio method to examine the structures, spectra, and stabilities of β-lactam (trans-13-azabicyclo[10.2.0]tetradecan-14-one), an important component of antibiotics. Based on the density functional theory (DFT) and second-order Møller-Plesset perturbation (MP2) methods, the present work demonstrated that forms I and II have isomorphic structures but can be distinguished by their Gibbs free energies and vibrational spectra. Forms I and II show a low-temperature polymorphic transformation at 308 K, where form I is stable below 308 K and form II is stable above 308 K. The proposed method suggests that the theoretical calculation can be used as a tool to effectively distinguish the isomorphic structures, and temperature-induced polymorphic transformation has far-reaching significance for drug storage and design.

Crystal Structure Optimization and Gibbs Free Energy Comparison of Five Sulfathiazole Polymorphs by the Embedded Fragment QM Method at the DFT Level

Molecular crystal plays an important role in many fields of science and technology, but it often crystallizes in different polymorphs with different physical properties. To guide the experimental synthesis of candidate materials, the atomic-scale model is frequently used to predict the most stable polymorph and its structural properties. Here, we show how an ab initio method can be used to achieve a rapid and accurate prediction of sulfathiazole crystal polymorphs (an antibiotic drug), based on the Gibbs free energy calculation and Raman spectra analysis. At the atmospheric pressure and the temperature of 300 K, we demonstrate that form III (FIII) is the most stable structure of sulfathiazole. The agreement between the predicted and experimental crystal structures corresponds to the order of stability for five sulfathiazole polymorphs as FI < FV < FIV < FII < FIII, which is achieved by employing the density functional theory (DFT) calculations.

Towards reliable ab initio sublimation pressures for organic molecular crystals - are we there yet?

Knowledge of molecular crystal sublimation equilibrium data is vital in many industrial processes, but this data can be difficult to measure experimentally for low-volatility species. Theoretical prediction of sublimation pressures could provide a useful supplement to experiment, but the exponential temperature dependence of sublimation (or any saturated vapor) pressure curve makes this challenging. An uncertainty of only a few percent in the sublimation enthalpy or entropy can propagate to an error in the sublimation pressure exceeding several orders of magnitude for a given temperature interval. Despite this fundamental difficulty, this paper performs some of the first ab initio predictions of sublimation pressure curves. Four simple molecular crystals (ethane, methanol, benzene, and imidazole) have been selected for a case study showing the currently achievable accuracy of quantum chemistry calculations. Fragment-based ab initio techniques and the quasi-harmonic approximation are used for calculations of cohesive and phonon properties of the crystals, while the vapor phase is treated by the ideal gas model. Ab initio sublimation pressure curves for model compounds are compared against their experimental counterparts. The computational uncertainties are estimated, weak points of the computational methodology are identified, and further improvements are proposed.

Phonons, phase transitions and thermal expansion in LiAlO2: an ab initio density functional study

We have used the ab initio density functional theory technique to understand the phase transitions and structural changes in various high temperature/pressure phases of LiAlO₂. The electronic band structure as well as phonon spectra is calculated for various phases as a function of pressure. The phonon entropy used for the calculations of Gibbs free energy is found to play an important role in the phase stability and phase transitions among various phases. A sudden increase in the polyhedral bond lengths (Li/Al-O) signifies the change from the tetrahedral to octahedral geometry at high-pressure phase transitions. The activation energy barrier for the high-pressure phase transitions is calculated. The phonon modes responsible for the phase transition (upon heating) from high pressure phases to ambient pressure phases are identified. Moreover, ab initio lattice dynamics calculations in the framework of quasi-harmonic approximations are used to study the anisotropic thermal expansion behavior of γ-LiAlO₂.

Theoretical predictions suggest carbon dioxide phases III and VII are identical

Accurate electronic structure calculations for the structures and simulated Raman spectra of high-pressure carbon dioxide suggest phases III and VII are identical, and the phase diagram should be revised.

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.

State-of-the-Art Calculations of Sublimation Enthalpies for Selected Molecular Crystals and Their Computational Uncertainty

A computational methodology for calculation of sublimation enthalpies of molecular crystals from first principles is developed and validated by comparison to critically evaluated literature experimental data. Temperature-dependent sublimation enthalpies for a set of selected 22 molecular crystals in their low-temperature phases are calculated. The computational methodology consists of several building blocks based on high-level electronic structure methods of quantum chemistry and statistical thermodynamics. Ab initio methods up to the coupled clusters with iterative treatment of single and double excitations and perturbative triples correction with an estimated complete basis set description [CCSD(T)/CBS] are used to calculate the cohesive energies of crystalline phases within a fragment-based additive scheme. Density functional theory (DFT) calculations with periodic boundary conditions (PBC) coupled with the quasi-harmonic approximation are used to evaluate the thermal contributions to the enthalpy of the solid phase. The properties of the vapor phase are calculated within the ideal-gas model using the rigid-rotor harmonic-oscillator model with correction for internal rotation using a one-dimensional hindered rotor approximation and a proper treatment of the molecular rotational degrees of freedom in the vicinity of 0 K. All individual terms contributing to the sublimation enthalpy as a function of temperature are discussed and their uncertainties estimated by comparison to critically evaluated experimental data.

Accurate prediction of energetic properties of ionic liquid clusters using a fragment-based quantum mechanical method

Accurate prediction of physicochemical properties of ionic liquids (ILs) is of great significance to understand and design novel ILs with unique properties.

Ab initio thermodynamic properties and their uncertainties for crystalline α-methanol

To investigate the performance of quasi-harmonic electronic structure methods for modeling molecular crystals at finite temperatures and pressures, thermodynamic properties are calculated for the low-temperature α polymorph of crystalline methanol and their computational uncertainties are analyzed.

Computational assignments of lattice vibrations of ice Ic

Herein,viathe first-principles density functional theory, CASTEP code, we investigated the 15 vibrational normal modes of ferroelectric hydrogen-ordered phase of ice Ic and the two peaks of hydrogen bond are clarified.

Predicting Molecular Crystal Properties from First Principles: Finite-Temperature Thermochemistry to NMR Crystallography

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 ¹H, ¹³C, and ¹⁵N nuclei. Due to their higher sensitivity to polarization effects, ¹⁷O 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.

Quantum mechanical calculation of electric fields and vibrational Stark shifts at active site of human aldose reductase

Recent advance in biophysics has made it possible to directly measure site-specific electric field at internal sites of proteins using molecular probes with C = O or C≡N groups in the context of vibrational Stark effect. These measurements directly probe changes of electric field at specific protein sites due to, e.g., mutation and are very useful in protein design. Computational simulation of the Stark effect based on force fields such as AMBER and OPLS, while providing good insight, shows large errors in comparison to experimental measurement due to inherent difficulties associated with point charge based representation of force fields. In this study, quantum mechanical calculation of protein's internal electrostatic properties and vibrational Stark shifts was carried out by using electrostatically embedded generalized molecular fractionation with conjugate caps method. Quantum calculated change of mutation-induced electric field and vibrational Stark shift is reported at the internal probing site of enzyme human aldose reductase. The quantum result is in much better agreement with experimental data than those predicted by force fields, underscoring the deficiency of traditional point charge models describing intra-protein electrostatic properties.

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.

Thermodynamic Properties of Molecular Crystals Calculated within the Quasi-Harmonic Approximation

A computational study of the possibilities of contemporary theoretical chemistry as regards calculated thermodynamic properties for molecular crystals from first-principles is presented. The study is performed for a testing set of 22 low-temperature crystalline phases whose properties such as densities of phonon states, isobaric heat capacities, and densities are computed as functions of temperature within the quasi-harmonic approximation. Electronic structure and lattice dynamics are treated by plane-wave based calculations with optPBE-vdW functional. Comparison of calculated results with reliable critically assessed experimental data is especially emphasized.