Predicting Molecular Crystal Properties from First Principles: Finite

This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes.

Article pubs.acs.org/accounts

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

1. INTRODUCTION

form screening in the search for potential formulations or problematic polymorphs.3,4 A large proportion of molecular crystal modeling studies in recent years have focused on lattice energy prediction without consideration of temperature. In 1995, Gavezzotti and Filippini argued that enthalpies dominate over entropies in determining polymorph stability.5 This is often true, but a recent survey of 508 polymorphic molecules found that simple harmonic free energy rankings would invert the relative room-temperature

The past decade has witnessed considerable progress in modeling molecular crystals with electronic structure theory. Advances in van der Waals dispersion, including density functional theory (DFT) approximations, periodic wave function methods, and fragment-based electronic structure methods, have improved the accuracy with which one can predict crystal structures and lattice energies.1 Success rates in the blind tests of crystal structure prediction have also improved.2 Crystal structure prediction is now being used by pharmaceutical companies to complement experimental solid © 2016 American Chemical Society

Received: August 3, 2016 Published: October 18, 2016 2501

DOI: 10.1021/acs.accounts.6b00404 Acc. Chem. Res. 2016, 49, 2501−2508

Article

Accounts of Chemical Research polymorph stabilities about 10% of the time.6 Since the free energy contributions frequently oppose contributions from the lattice energy,6 finite-temperature free energy effects will matter more often when considering broader regions of the temperature/pressure phase diagram. The next frontier in molecular crystal modeling involves moving beyond modeling structures and lattice energies at 0 K to predicting a range of properties at realistic temperatures and pressures. Temperature effects in molecular crystals extend beyond standard thermal contributions to the statistical mechanical partition functions. Heating a molecular crystal to room temperature expands its unit cell volume by up to 15%. This expansion weakens the noncovalent interactions between the molecules, and it frequently softens the lattice phonon modes.6,7 These changes alter the temperature-dependence of heat capacities, enthalpies, entropies, and free energies. They can soften the crystal, reducing its bulk modulus. Changing the intermolecular separations can substantially alter the transfer integrals in organic semiconductor materials.8 As theory is increasingly used to predict such properties in molecular crystals, it becomes more important to address these finitetemperature effects. Nuclear magnetic resonance (NMR) crystallography9−12 represents another important frontier in molecular crystal modeling. X-ray diffraction (XRD) remains the standard technique for solving crystal structures, but obtaining large single crystals or solving the powder XRD structure is not always feasible. Solid-state NMR provides a complementary alternative to XRD techniques which probes local crystal packing details. NMR experiments can also be performed on impure samples or in situ. However, translating observed chemical shifts into a three-dimensional crystal structure is challenging. In practice, NMR crystallography often relies on computationally predicted chemical shifts to identify a structure that is consistent with the experimental spectrum. This requires predicting chemical shifts with sufficient accuracy to discriminate the often subtle variations that occur upon changes in crystal packing or unit cell volume. To address these sorts of challenges, we have developed fragment-based electronic structure models1,13 which enable accurate electronic structure methods to be employed in molecular crystal problems. Here, we first focus on modeling thermal expansion and evaluate its impact on thermochemistry and mechanical properties. We have recently demonstrated that coupling fragment second-order Møller−Plesset perturbation theory (MP2) or coupled cluster singles, doubles, and perturbative triples (CCSD(T)) models to the quasi-harmonic approximation for thermal expansion allows one to predict finite temperature crystal structures, thermochemistry, and bulk moduli in excellent agreement with experiment. This model has enabled us to investigate the importance of thermal expansion in the prediction molecular crystal structures and properties for several small molecule crystals. Second, we examine the performance of fragment-based NMR chemical shift calculation in the context of NMR crystallography. The gauge-including projector augmented wave (GIPAW) planewave DFT approach is widely used for chemical shift prediction in the solid state.9,10 However, fragment-based calculations can offer performance on par with or better than GIPAW, particularly through the use of hybrid density functionals. We will review examples of chemical shift prediction in biologically and pharmaceutically relevant

crystals that demonstrate the increased discrimination among candidate structures one can obtain from fragment methods.

2. CRYSTAL PROPERTIES AT FINITE TEMPERATURES 2.1. Quasi-Harmonic Hybrid Many-Body Interaction Model

To model crystals at finite temperatures and pressures, one must focus on the Gibbs free energy G at the current temperature (T) and pressure (P), G(T , P) = E + PV + Fvib(T )

(1)

instead of the electronic energy E. This involves computing a pressure−volume term (PV) and the Helmholtz vibrational free energy Fvib. The electronic energy is computed via the fragment hybrid many-body interaction (HMBI) model,14−17 in which intramolecular and short-range pairwise intermolecular interactions are treated quantum mechanically, while long-range and manybody interactions are approximated using a polarizable force field, QM MM MM E = E1QM −body + ESR2−body + E LR2−body + Emany−body

(2)

This model allows high-level electronic structure methods like MP2 or CCSD(T) to be employed on molecular crystals with reasonable computational cost. The HMBI model predicts crystal lattice energies15,16 and structures17 in good agreement with data inferred from experiments. It has been used to investigate polymorphism in crystals such as aspirin,18 oxalyl dihydrazide,19 and ice XV.20 Performing structure optimizations and lattice dynamical phonon calculations at each new crystal volume to obtain Fvib and G(T,P) would be very computationally demanding. Instead, the quasi-harmonic approximation (QHA) provides an effective tool for evaluating the vibrational contributions in small-molecule crystals. The QHA estimates the phonon frequencies ωi at the current unit cell volume V based on a set of reference frequencies ωref i computed at a reference cell volume Vref. ⎛ V ⎞−γi ωi = ωiref ⎜ ref ⎟ ⎝V ⎠

(3)

The volume dependence of the phonon frequencies is approximated using mode-specific Grüneisen parameters γi that are computed via finite difference, ⎛ ∂ln ωi ⎞ ⎟ γi = −⎜ ⎝ ∂ln V ⎠

(4)

Using these estimated phonon frequencies, harmonic vibrational partition functions and vibrational contributions to the free energy (Fvib) can be computed readily. The QHA captures some anharmonicity through the volume dependence of the harmonic phonon modes. The quasi-harmonic approximation allows optimization of the crystal structure and computation of other crystal properties as a function of temperature and pressure. Even with the HMBI fragment approach, computing the reference phonons and Grüneisen parameters requires significant computational effort initially, especially at the MP2 level. However, subsequent quasi-harmonic free-energy structure optimizations incur computational costs that are comparable to a standard electronic energy structure optimization. See ref 7 for details. Exploitation of space-group symmetry can significantly reduce 2502


Article

Accounts of Chemical Research

Figure 1. Predicted (a) volume, (b) sublimation enthalpy, (c) sublimation entropy, and (d) errors in the free energies for carbon dioxide (phase I) as a function of temperature. Solid lines include quasi-harmonic thermal expansion, while dotted lines neglect it. No QHA refers to standard electronic energy calculations, with no quasi-harmonic contributions. Adapted with permission from ref 7. Copyright 2016 International Union of Crystallography.

Third, thermochemical properties can also be predicted accurately when thermal expansion is treated. Large basis MP2 and CCSD(T) calculations predict the sublimation enthalpy (Figure 1b) and entropy (Figure 1c) results in very good agreement with experiment.23 Favorable results were obtained for the other crystals we have examined, too.7 For the enthalpy of solid carbon dioxide in particular, the quasi-harmonic model captures the temperature dependence very well. However, the lattice energy is consistently slightly too large, causing the sublimation enthalpies to be systematically ∼1.5 kJ/mol too large. From the quasi-harmonic free energy of sublimation, one predicts a sublimation temperature of 198.0 K, in excellent agreement with the experimental value of 194.7 K.7 Neglecting thermal expansion causes overestimation of the sublimation enthalpy/entropy at higher temperatures, as seen in Figure 1b and c.23 For the enthalpy, thermal expansion typically reduces the lattice energy of the crystal (greater intermolecular spacing weakens the intermolecular interactions), and neglecting this expansion causes the lattice energy to be spuriously large at higher temperatures.7 The overestimation of the sublimation entropy at higher temperatures arises from the changes in the lattice phonon modes, which frequently soften with thermal expansion. Neglecting thermal expansion therefore typically leads to artificially high lattice phonon frequencies, which causes the isobaric heat capacity of the crystal to be underestimated at higher temperatures. The underestimated heat capacity in turn reduces the predicted entropy of the crystal Scrystal and increases sublimation entropy, ΔSsub = Sgas − Scrystal.

the computational effort required to perform the optimization and lattice dynamics phonon calculations,21 and fragment-based supercell Hessian calculations add only marginal overhead compared to conventional Γ-point only Hessian calculation,22 enabling facile treatment of phonon dispersion. 2.2. Importance of Thermal Expansion

Using quasi-harmonic HMBI fragment modeling, we have investigated the impact of thermal expansion on thermochemical and mechanical crystal properties in several small molecule crystals: carbon dioxide, ice, acetic acid, and imidazole. Several interesting conclusions can be drawn from these studies. First, the quasi-harmonic approximation can accurately describe important thermal expansion effects up to room temperature in these rigid molecules.7 As shown in Figure 1a, large-basis MP2 and CCSD(T) underestimate the crystal volume only slightly and predict its temperature dependence well.23 Second, zero-point vibrational energy often introduces a substantial volume expansion in molecular crystals. In the crystals we have examined, it ranges from 30−60% of the overall expansion between 0 K and 200−300 K.7 In carbon dioxide, for example, about 30% of the overall expansion between the optimal electronic energy structure and the volume at the sublimation point (195 K) stems from zero-point contributions (Figure 1a). While this phenomenon is well understood in principle, it is often ignored in practice. Notably, even when low-temperature experimental structures are available for comparison, one should not seek perfect agreement between predicted and experimental structures unless vibrational effects are treated appropriately. 2503


Article


Figure 2. Quasi-harmonic fragment MP2/CBS predicts (red) the phase I carbon dioxide (a) pressure−volume curves and (b) bulk modulus in good agreement with experiment (gray). Neglecting thermal expansion leads to underestimation of the volume and substantial overestimation of the room-temperature bulk modulus. The experimental phase change of CO2 (yellow region) is not modeled here. Adapted with permission from ref 23. Published 2016 Royal Society of Chemistry.

Fourth, based on the discussion above, the errors in ΔHsub and ΔSsub introduced by neglecting thermal expansion frequently exhibit the same sign, which means that they partially cancel in ΔGsub. Figure 1d plots the errors in these three thermodynamic quantities with and without thermal expansion. Although the errors in ΔHsub and ΔSsub when neglecting thermal expansion become appreciable near the sublimation point, the corresponding errors in the Gibbs free energy of sublimation for carbon dioxide are smaller and are comparable in magnitude whether or not thermal expansion is modeled. Indeed, the predicted sublimation temperature when neglecting thermal expansion is fortuitously accurate at 190.4 K,7 versus 194.7 K experimentally. Fifth, thermal expansion also contributes significantly to mechanical properties like the bulk modulus. Figure 2a plots the predicted and experimental pressure−volume curves for carbon dioxide at 300 K.23 Up to ∼10 GPa, where phase I carbon dioxide undergoes a transition to phase III,24 the quasiharmonic model with thermal expansion matches the experimental data very accurately. Interestingly, neglecting thermal expansion leads to visible underestimation of the volume even at high pressures. Fitting pressure−volume curves obtained at a given temperature to the Vinet equation of state allows prediction of the bulk modulus (Figure 2b). Once again, thermal expansion plays a substantial role here, with the bulk modulus decreasing several fold between 0 K and room temperature. Though the bulk moduli obtained from various experimental studies exhibit some variability, the fragment MP2/CBS predictions are in generally good agreement across all temperatures.23 Bulk modulus predictions which neglect thermal expansion (“No QHA” in Figure 2b) substantially overestimate the bulk modulus at room temperature.23,25,26 As these examples demonstrate, thermal expansion can play an appreciable role in predicting a range of crystal properties. Its effects on free energies are modest (a couple kJ/mol or less) in the examples we have studied, but such differences are sometimes enough to rerank polymorph stabilities. The impact on other properties that depend on the structure of the crystal, such as mechanical properties, can be larger. The quasiharmonic approximation provides an effective tool for modeling thermal expansion and describing these temperature-dependent

properties, at least for rigid molecules. Future work will need to investigate the performance of quasi-harmonic models for crystals involving larger, more flexible species. Doing so will also require using lower-cost electronic structure methods, particularly for computing the reference phonon frequencies and Grüneisen parameters which typically form the computational bottleneck.7

3. NMR CRYSTALLOGRAPHY 3.1. Fragment-Based NMR Chemical Shifts

Fragment-based techniques are also promising for solid-state NMR chemical shift prediction in molecular crystals. Our group is not the first to use fragment approaches for NMR. As discussed in a recent review,1 early models for predicting chemical shifts in molecular crystals used charge-embedded molecules, while various embedded fragment approaches are used in biological systems. More recently, plane wave GIPAW DFT calculations (usually employing the PBE functional) have become the tool of choice for predicting chemical shifts in molecular crystals.9,10 Interest in finite cluster and fragment models has renewed in the past few years27−34 The NMR chemical shielding tensor σ on nucleus A can be written as the second derivative of the energy with respect to the magnetic field B and the nuclear magnetic moment μA,

σA =

∂ 2E ∂B∂μ A

(5)

Accordingly, a many-body expansion for the shielding tensor is derived simply by differentiating the standard many-body expansion for the energy. This chemical shielding tensor expansion converges rapidly, especially when electrostatic embedding is employed.30 For 1H, 13C, and 15N, a one- and two-body fragment expansion (i.e., monomer and dimer chemical shift calculations out to 6 Å) combined with electrostatic point-charge embedding to capture long-range and many-body effects allows accurate chemical shift prediction.32,33 This model is similar to the binary interaction fragment model22 Hirata and co-workers have applied to other molecular crystal problems.35 2504


Article


PBE approach improves the fragment result considerably, to the 8.8 ppm shown in Figure 3. Second, switching from PBE to a hybrid density functional like PBE0 or B3LYP improves the statistical accuracy of the predicted 13C or 15N shifts by up to ∼30%. For 1H, the differences among functionals are negligible. For 17O, cluster/ fragment PBE0 reduces the error to 7.6 ppm, which is comparable with the 7.2 ppm GIPAW PBE result, especially since the experimental uncertainties in the isotropic 17O shifts studied range ∼0.5−5 ppm. In other word, even though fragment-based approaches behave less well for oxygen nuclei, the improved accuracy from PBE0 compensates for errors in the description of many-body effects. The higher accuracy of hybrid functionals relative to GGAs is perhaps the chief advantage of fragment methods. Hybrid density functionals are generally computationally cost-prohibitive in a plane wave basis set, but they can be used routinely in the Gaussian orbital fragment approach. Thus, fragment methods provide a practical route to higher-accuracy chemical shifts. Computationally, the two-body fragment model scales linearly with the number of symmetrically unique molecules in the crystal, Z′, and it is independent of the total number of molecules Z in the cell. Each individual fragment job can be run over multiple processors, and multiple fragment jobs can be run simultaneously. This enables the fragment approach to scale effectively to hundreds of processor cores, and often allows one to compute chemical shifts for even very large organic crystals within a few hours of wall time. Chemical shift referencing is an important issue in modeling. Theory predicts absolute chemical shieldings, but experiments measure the chemical shift, which is the difference between the chemical shielding of the compound of interest and that of a standard. Many referencing strategies exist,36,37 but a common one involves mapping between the calculated shieldings σi and observed shifts δi via a linear regression,

Many-body effects are more significant for the shielding tensors of oxygen nuclei. However, these effects can be captured using a hybrid embedded cluster/fragment method in which one performs chemical shift calculations on a cluster consisting of the molecule of interest and its nearest neighbors (∼10−15 molecules), followed by inexpensive fragment twobody corrections for longer-range interactions. In principle, the fragment chemical shielding contributions can be computed with any electronic structure method. Even with DFT, fragment approaches can offer a number of advantages, as discussed in the next section. 3.2. Performance of the Fragment Approach

To evaluate these fragment chemical shift methods, we have performed benchmarking on 63 molecular crystals with experimentally known structures and 328 isotropic 1H, 13C, 15 N, or 17O chemical shifts spanning much of the practical chemical shift range observed in organic compounds. These studies have revealed a number of important features. First, for a generalized gradient approximation (GGA) density functional like PBE, the fragment and GIPAW approaches perform similarly for 1H, 13C, and 15N (Figure 3). The most notable

Figure 3. Root-mean-square errors vs experiment in solid-state NMR isotropic chemical shifts predicted using GIPAW PBE (green), fragment PBE (blue), and fragment PBE0 (red) for four molecular crystal benchmark sets.33 The 2-body fragment method was used for 1 H, 13C, 15N, and the cluster/fragment approach was used for 17O. Note: the fragment PBE and PBE0 lines for 1H are superimposed.

δi = aσi + b

(6)

where a and b are empirical parameters. Our benchmark studies provide robust sets of a and b parameters for mapping shieldings to shifts in an explicit condensed-phase chemical environment (as distinct from gas-phase or polarizable continuum environments). These regression parameters are robust with respect to the composition of the training set, as demonstrated via statistical cross-validation. The regression parameters can then be applied to new systems with no additional empiricism. They have been used,

exception occurs for 17O, where the increased magnitude of the many-body terms30 leads to reduced performance with the 2body fragment approach (RMS error 11.6 ppm, vs 7.2 ppm for GIPAW PBE33). Switching to a combined cluster/fragment

Figure 4. Comparison of predicted and experimental 13C spectra for the three polymorphs of acetaminophen using (a) fragment PBE0 and (b) GIPAW PBE. Disagreements in assignment are highlighted in green. Adapted with permission from ref 34. Copyright 2016 American Chemical Society. 2505


Article


Figure 5. Reduced χ̃2 values for blind assignment of the monomers in (a) the three polymorphs of acetaminophen or (b) two forms of testosterone. Compared to GIPAW or fragment PBE, fragment PBE0 provides enhanced discrimination among the best (red) and other assignments. Adapted with permission from ref 34. Copyright 2016 American Chemical Society.

Figure 6. Comparison of predicted and experimental low-frequency region 13C solid-state NMR spectra of α-testosterone using (a) Fragment PBE0 and (b) GIPAW PBE. Adapted with permission from ref 34. Copyright 2016 American Chemical Society.

ppm). Both models generally predict the same 13C assignments, except for a couple shifts which are separated by less than 1 ppm (i.e., below the expected prediction accuracy). This increased accuracy results in better discrimination among different crystallographic environments. To see this, consider a “blind” assignment of the predicted shifts for the for crystallographically unique monomers (one each from forms I and II, and two from form III). There are 24 possible ways to assign the four unique monomers to four sets of shifts, and Figure 5a compares the reduced χ̃2 values,

for example, to study molecular crystal polymorphs (section 3.3) and to solve the structure for photoreacted dimers of 9tert-butyl anthracene dimers.38 The scaling parameters can also be applied to nonperiodic systems; they helped interpret the chemical shifts for substrates bound in the β-subunit active site of the enzyme tryptophan synthase.39 3.3. Polymorph Discrimination

Solid-state chemical shift prediction is primarily used to interpret experimental NMR spectra or infer crystal structures. For example, one might compare predicted and experimental chemical shifts to help assign an NMR spectrum, to rule out candidate structures generated via crystal structure prediction, or to refine experimental crystal structures by relaxing it to lower its energy and improve agreement with the NMR shifts. The fundamental challenge in each of these examples lies in discriminating between correct and incorrect structures/ assignments based on how chemical shifts vary in distinct crystallographic environments. Increased discrimination can be achieved through the use of different magnetically active nuclei, multidimensional NMR correlation experiments, or full shielding tensor measurements (instead of only isotropic shifts). The increased accuracy offered by fragment methods and hybrid density functionals provides another, complementary route to enhance the resolving power of theoretical predictions. Here, we highlight two examples where higher accuracy fragment PBE0 calculations provide increased discrimination compared to GIPAW PBE: acetaminophen and testosterone. Acetaminophen is one of the most widely used pharmaceuticals. For the three known polymorphs, Figure 4 shows that fragment PBE0 13C shifts offer improved agreement with experiment compared to GIPAW PBE (RMS error 1.2 vs 1.7

χ̃2 =

1 N

N

∑ i

(δipred − δiexpt)2 σrms 2

(7)

for each. This expression sums the square of the error in each shift divided by the RMS error from the expected error distribution (determined from the benchmark tests in Section 3.2). All methods shown predict the same overall assignment as having the smallest χ̃2. However, the larger GIPAW PBE errors lead to smaller differences in the χ̃2 values among the different potential assignments. In contrast, fragment PBE0 increases the separation between the different χ̃2 values, thereby increasing discrimination. The hardest discrimination occurs between the two crystallographically inequivalent monomers in form III (red and blue lines in Figure 5a), whose resonances are ill-resolved experimentally. Chemical shift predictions have also proved helpful in testosterone. In the α crystal form, the asymmetric unit consists of two inequivalent monomers u and v with 19 carbons each. Assigning all 38 13C shifts was feasible only by coupling twodimensional 13C NMR experiments and GIPAW chemical shift predictions.40 Figure 6 shows the improved agreement between theory and experiment with fragment PBE0 for the dense low2506



■

frequency region of the α-testosterone 13C spectrum. Combining results for the α and β (monohydrate) forms, the RMS error with fragment PBE0 is a third smaller than that obtained with GIPAW PBE (2.1 vs 3.3 ppm). Like in acetaminophen, fragment PBE0 provides improved discrimination (Figure 5b) among different possible assignments of the three crystallographically unique monomers found in the two forms.

ACKNOWLEDGMENTS This research was made possible through funding from the National Science Foundation (CHE-1362465), and supercomputer time is from XSEDE (TG-CHE110064). Collaborations with G. M. Day, L. J. Mueller, and B. Schatschneider on the NMR studies are gratefully acknowledged.

■

REFERENCES

(1) Beran, G. J. O. Modeling polymorphic molecular crystals with electronic structure theory. Chem. Rev. 2016, 116, 5567−5613. (2) Reilly, A. M.; Cooper, R. I.; Adjiman, C. S.; Bhattacharya, S.; Boese, A. D.; Brandenburg, J. G.; Bygrave, P. J.; Bylsma, R.; Campbell, J. E.; Car, R.; Case, D. H.; Chadha, R.; Cole, J. C.; Cosburn, K.; Cuppen, H. M.; Curtis, F.; Day, G. M.; DiStasio, R. A., Jr; Dzyabchenko, A.; van Eijck, B. P.; Elking, D. M.; van den Ende, J. A.; Facelli, J. C.; Ferraro, M. B.; Fusti-Molnar, L.; Gatsiou, C.-A.; Gee, T. S.; de Gelder, R.; Ghiringhelli, L. M.; Goto, H.; Grimme, S.; Guo, R.; Hofmann, D. W. M.; Hoja, J.; Hylton, R. K.; Iuzzolino, L.; Jankiewicz, W.; de Jong, D. T.; Kendrick, J.; de Klerk, N. J. J.; Ko, H.Y.; Kuleshova, L. N.; Li, X.; Lohani, S.; Leusen, F. J. J.; Lund, A. M.; Lv, J.; Ma, Y.; Marom, N.; Masunov, A. E.; McCabe, P.; McMahon, D. P.; Meekes, H.; Metz, M. P.; Misquitta, A. J.; Mohamed, S.; Monserrat, B.; Needs, R. J.; Neumann, M. A.; Nyman, J.; Obata, S.; Oberhofer, H.; Oganov, A. R.; Orendt, A. M.; Pagola, G. I.; Pantelides, C. C.; Pickard, C. J.; Podeszwa, R.; Price, L. S.; Price, S. L.; Pulido, A.; Read, M. G.; Reuter, K.; Schneider, E.; Schober, C.; Shields, G. P.; Singh, P.; Sugden, I. J.; Szalewicz, K.; Taylor, C. R.; Tkatchenko, A.; Tuckerman, M. E.; Vacarro, F.; Vasileiadis, M.; Vazquez-Mayagoitia, A.; Vogt, L.; Wang, Y.; Watson, R. E.; de Wijs, G. A.; Yang, J.; Zhu, Q.; Groom, C. R. Report on the sixth blind test of organic crystal structure prediction methods. Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater. 2016, 72, 439−459. (3) Price, S. L.; Reutzel-Edens, S. M. The potential of computed crystal energy landscapes to aid solid-form development. Drug Discovery Today 2016, 21, 912−923. (4) Neumann, M. A.; van de Streek, J.; Fabbiani, F. P. A.; Hidber, P.; Grassmann, O. Combined crystal structure prediction and highpressure crystallization in rational pharmaceutical polymorph screening. Nat. Commun. 2015, 6, 7793. (5) Gavezzotti, A.; Filippini, G. Polymorphic Forms of Organic Crystals at Room Conditions: Thermodynamic and Structural Implications. J. Am. Chem. Soc. 1995, 117, 12299−12305. (6) Nyman, J.; Day, G. M. Static and lattice vibrational energy differences between polymorphs. CrystEngComm 2015, 17, 5154− 5165. (7) Heit, Y. N.; Beran, G. J. O. How important is thermal expansion for predicting molecular crystal structures and thermochemistry at finite temperatures? Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater. 2016, 72, 514−529. (8) Coropceanu, V.; Cornil, J.; da Silva Filho, D. A.; Olivier, Y.; Silbey, R.; Brédas, J.-L. Charge Transport in Organic Semiconductors. Chem. Rev. 2007, 107, 926−952. (9) Ashbrook, S. E.; McKay, D. Combining Solid-State NMR Spectroscopy with First-Principles Calculations - A Guide to NMR Crystallography. Chem. Commun. 2016, 52, 7186−7204. (10) Bonhomme, C.; Gervais, C.; Babonneau, F.; Coelho, C.; Pourpoint, F.; Azaïs, T.; Ashbrook, S. E.; Griffin, J. M.; Yates, J. R.; Mauri, F.; Pickard, C. J. First-Principles Calculation of NMR Parameters Using the Gauge Including Projector Augmented Wave Method: A Chemist’s Point of View. Chem. Rev. 2012, 112, 5733− 5779. (11) Salager, E.; Day, G. M.; Stein, R. S.; Pickard, C. J.; Elena, B.; Emsley, L. Powder Crystallography by Combined Crystal Structure Prediction and High-Resolution 1H Solid-State NMR Spectroscopy. J. Am. Chem. Soc. 2010, 132, 2564−2566. (12) Dudenko, D. V.; Williams, P. A.; Hughes, C. E.; Antzutkin, O. N.; Velaga, S. P.; Brown, S. P.; Harris, K. D. M. Exploiting the Synergy of Powder X-ray Diffraction and Solid-State NMR Spectroscopy in

4. CONCLUSIONS In summary, rapid advances are currently being made in the prediction of experimentally observable quantities in molecular crystals. As demonstrated here, fragment-based methods offer a number of advantages for predicting such properties accurately. More theoretical work is still needed, however. Even with the fragment approach, the need for phonons makes the quasiharmonic calculations expensive computationally compared to conventional energy and gradient calculations. Can these phonons be computed more approximately? In addition, how well does the quasi-harmonic behave in systems where anharmonicity is particularly large? On the NMR side, can the quality of the chemical shift calculations be improved via better electrostatic embedding and/or treatment of Pauli repulsion?41 The simple point-charges used here do not behave as well for 17O chemical shifts, for which many-body effects are important. Perhaps a better description of the environment could enable more uniformly accurate predictions for all nuclei.

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. Biographies Gregory J. O. Beran was born near San Francisco, California. He studied chemistry at the University of California San Diego (B.S. 2000) and went on to earn a Ph.D. with Martin Head-Gordon at the University of California Berkeley (2005). After 2 years of postdoctoral research in Chemical Engineering with William H. Green at the Massachusetts Institute of Technology, he started as a faculty member at the University of California Riverside. His research focuses on the development and application of fragment-based treatments of molecular crystals and electronic structure methods for noncovalent interactions. Joshua D. Hartman was born in Wheat Ridge, Colorado. He studied biochemistry at Colorado State University in Fort Collins, Colorado (B.S. 2007) and the University of Denver (M.S. 2009). He went on to earn a Ph.D. in chemistry with Gregory Beran at the University of California Riverside (2015). He currently holds a postdoctoral research position at the University of California Riverside. His research focuses on the development and application of computational techniques aimed a predicting magnetic properties in molecular crystals. Yonaton N. Heit graduated from Fairleigh Dickinson University with a bachelor’s in chemistry in 2010. In 2016, he completed his Ph.D. in chemistry from University of California, Riverside with Dr. Gregory Beran. His research focused on modeling molecular crystal thermal expansion, thermodynamic properties at finite temperatures, and exploitation of space group symmetry to accelerate fragment based calculations. 2507


Article

Accounts of Chemical Research Structure Determination of Organic Molecular Solids. J. Phys. Chem. C 2013, 117, 12258−12265. (13) Wen, S.; Nanda, K.; Huang, Y.; Beran, G. J. O. Practical quantum mechanics-based fragment methods for predicting molecular crystal properties. Phys. Chem. Chem. Phys. 2012, 14, 7578−7590. (14) Beran, G. J. O. Approximating quantum many-body intermolecular interactions in molecular clusters using classical polarizable force fields. J. Chem. Phys. 2009, 130, 164115. (15) Beran, G. J. O.; Nanda, K. Predicting organic crystal lattice energies with chemical accuracy. J. Phys. Chem. Lett. 2010, 1, 3480− 3487. (16) Wen, S.; Beran, G. J. O. Accurate molecular crystal lattice energies from a fragment QM/MM approach with on-the-fly ab initio force-field parameterization. J. Chem. Theory Comput. 2011, 7, 3733− 3742. (17) Nanda, K.; Beran, G. J. O. Prediction of organic molecular crystal geometries from MP2-level fragment quantum mechanical/ molecular mechanical calculations. J. Chem. Phys. 2012, 137, 174106. (18) Wen, S.; Beran, G. J. O. Accidental degeneracy in crystalline aspirin: New insights from high-level ab initio calculations. Cryst. Growth Des. 2012, 12, 2169−2172. (19) Wen, S.; Beran, G. J. O. Crystal polymorphism in oxalyl dihydrazide: Is empirical DFT-D accurate enough? J. Chem. Theory Comput. 2012, 8, 2698−2705. (20) Nanda, K.; Beran, G. J. O. What governs the proton-ordering in ice XV? J. Phys. Chem. Lett. 2013, 4, 3165−3169. (21) Heit, Y.; Beran, G. J. O. Exploting space group symmetry in fragment-based molecular crystal calculations. J. Comput. Chem. 2014, 35, 2205−2214. (22) Hirata, S. Fast electron-correlation methods for molecular crystals: an application to the α, β1, and β2 modifications of solid formic acid. J. Chem. Phys. 2008, 129, 204104. (23) Heit, Y. N.; Nanda, K. D.; Beran, G. J. O. Predicting finitetemperature properties of crystalline carbon dioxide from first principles with quantitative accuracy. Chem. Sci. 2016, 7, 246−255. (24) Aoki, K.; Yamawaki, H.; Sakashita, M.; Gotoh, Y.; Takemura, K. Crystal Structure of the High-Pressure Phase of Solid CO2. Science 1994, 263, 356−358. (25) Li, J.; Sode, O.; Hirata, S. Second-Order Many-Body Perturbation Study on Thermal Expansion of Solid Carbon Dioxide. J. Chem. Theory Comput. 2015, 11, 224−229. (26) Gohr, S.; Grimme, S.; Söhnel, T.; Paulus, B.; Schwerdtfeger, P. Pressure dependent stability and structure of carbon dioxide-A density functional study including long-range corrections. J. Chem. Phys. 2013, 139, 174501. (27) Holmes, S. T.; Iuliucci, R. J.; Mueller, K. T.; Dybowski, C. Density functional investigation of intermolecular effects on 13C NMR chemical-shielding tensors modeled with molecular clusters. J. Chem. Phys. 2014, 141, 164121. (28) Holmes, S. T.; Iuliucci, R. J.; Mueller, K. T.; Dybowski, C. Critical Analysis of Cluster Models and Exchange-Correlation Functionals for Calculating Magnetic Shielding in Molecular Solids. J. Chem. Theory Comput. 2015, 11, 5229−5241. (29) Reid, D. M.; Collins, M. A. Calculating nuclear magnetic resonance shieldings using systematic molecular fragmentation by annihilation. Phys. Chem. Chem. Phys. 2015, 17, 5314−5320. (30) Hartman, J. D.; Beran, G. J. O. Fragment-based electronic structure approach for computing nuclear magnetic resonance chemical shifts in molecular crystals. J. Chem. Theory Comput. 2014, 10, 4862−4872. (31) Hartman, J. D.; Neubauer, T. J.; Caulkins, B. G.; Mueller, L. J.; Beran, G. J. O. Converging nuclear magnetic shielding calculations with respect to basis and system size in protein systems. J. Biomol. NMR 2015, 62, 327−340. (32) Hartman, J. D.; Monaco, S.; Schatschneider, B.; Beran, G. J. O. Fragment-based 13C nuclear magnetic resonance chemical shift predictions in molecular crystals: An alternative to plane-wave methods. J. Chem. Phys. 2015, 143, 102809.

(33) Hartman, J. D.; Kudla, R. A.; Day, G. M.; Mueller, L. J.; Beran, G. J. O. Benchmark fragment-based 1H, 13C, 15N and 17O chemical shift predictions in molecular crystals. Phys. Chem. Chem. Phys. 2016, 18, 21686−21709. (34) Hartman, J. D.; Beran, G. J. O. Enhanced NMR discrimination of pharmaceutically relevant molecular crystal forms through fragmentbased ab initio chemical shift predictions. Cryst. Growth Des. 2016, DOI: 10.1021/acs.cgd.6b01157. (35) Hirata, S.; Gilliard, K.; He, X.; Li, J.; Sode, O. Ab initio molecular crystal structures, spectra, and phase diagrams. Acc. Chem. Res. 2014, 47, 2721−30. (36) Harris, R. K.; Hodgkinson, P.; Pickard, C. J.; Yates, J. R.; Zorin, V. Chemical shift computations on a crystallographic basis: Some reflections and comments. Magn. Reson. Chem. 2007, 45, S174−S186. (37) Lodewyk, M. W.; Siebert, M. R.; Tantillo, D. J. Computational prediction of 1H and 13C chemical shifts: a useful tool for natural product, mechanistic, and synthetic organic chemistry. Chem. Rev. 2012, 112, 1839−1862. (38) Yang, C.; Zhu, L.; Kudla, R. A.; Hartman, J. D.; Al-Kaysi, R. O.; Monaco, S.; Schatschneider, B.; Magalhaes, A.; Beran, G. J. O.; Bardeen, C. J.; Mueller, L. J. Crystal structure of the meta-stable intermediate in the photomechanical, crystal-to-crystal reaction of 9tertbutyl anthracene ester. CrystEngComm 2016, 18, 7319−7329. (39) Young, R. P.; Caulkins, B. G.; Borchardt, D.; Bulloch, D. N.; Larive, C. K.; Dunn, M. F.; Mueller, L. J. Solution-State 17O Quadrupole Central-Transition NMR Spectroscopy in the Active Site of Tryptophan Synthase. Angew. Chem., Int. Ed. 2016, 55, 1350− 1354. (40) Harris, R. K.; Joyce, S. A.; Pickard, C. J.; Cadars, S.; Emsley, L. Assigning carbon-13 NMR spectra to crystal structures by the INADEQUATE pulse sequence and first principles computation: a case study of two forms of testosterone. Phys. Chem. Chem. Phys. 2006, 8, 137−143. (41) Cui, Q.; Karplus, M. Molecular Properties from Combined QM/MM Methods. 2. Chemical Shifts in Large Molecules. J. Phys. Chem. B 2000, 104, 3721−3743.

2508


Predicting Molecular Crystal Properties from First Principles: Finite

Recommend Documents