Predicting Organic Crystal Lattice Energies with Chemical Accuracy

pubs.acs.org/JPCL

Predicting Organic Crystal Lattice Energies with Chemical Accuracy Gregory J. O. Beran* and Kaushik Nanda Department of Chemistry, University of California, Riverside, California 92521, United States

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

M

olecular crystal structure impacts the properties of pharmaceuticals, organic photovoltaics, and energetic materials. It also determines solid-state reactivity, which can be exploited for highly controlled, solvent-free chemistry. Unfortunately, a priori crystal structure prediction remains difficult.1 Ab initio quantum chemistry techniques promise to overcome the limitations of classical force-field methods, but accuracy and computational cost issues remain. Periodic density functional theory (DFT) is the most obvious technique to use, but it requires corrections for dispersion interactions. Dispersion corrections used in molecular crystals are usually empirical;2-6 although more expensive, nonempirical dispersion-corrected models have also been applied.7-9 Alternatively, second-order Møller-Plesset perturbation theory (MP2) captures dispersion naturally, and its local periodic implementations have become practical in recent years.10-14 However, periodic MP2 remains limited to unit cells with around 50 atoms in a medium-sized basis set, and post-MP2 correlation effects can be important. To overcome these difficulties, we have extended the fragment-based quantum/classical hybrid many-body interaction (HMBI) model15,16 to infinite systems. HMBI-predicted lattice energies for several different molecular crystals lie within a few kJ/mol of the experimental values, demonstrating that highly accurate lattice energy predictions are practical in small-molecule organic crystals. Fragment-based methods, which decompose a molecular cluster or crystal into a series of interacting molecules, greatly reduce the computational cost and allow essentially MP2- or coupled-cluster-quality (CCSD(T)) results even for very large systems. Other fragment methods have been adapted to treat

r 2010 American Chemical Society

infinite liquids or crystals. Examples include the fragment molecular orbital (FMO) method,17 Hirata's binary interaction model,18,19 and the incremental methods from the Stoll20,21 and Manby22 groups. One such approach predicts the lattice energy of ice very accurately by combining a periodic HartreeFock (HF) calculation with a real-space treatment of one- and two-body electron-electron correlation effects.23 This same approach has also been applied to rare gas crystals.24 A similar approach that uses periodic DFT instead of periodic HF has been applied to benzene25 and other crystals.26 Hirata and co-workers applied the binary interaction method to examine the structures and phonon dispersion curves of formic acid and hydrogen fluoride.27,28 DFT-based symmetry-adapted perturbation theory (SAPT(DFT)) represents a somewhat different fragment approach that has been used very successfully to parametrize force fields for crystal structure prediction and to predict crystal lattice energies.29-31 Our HMBI approach for modeling systems of interacting molecules uses standard quantum mechanical methods to treat the individual molecules (one-body terms) and their pairwise interactions (two-body terms). The less important but non-negligible three-body and higher “many-body” terms are approximated with a classical molecular mechanics (MM) polarizable force field, QM MM EHMBI ¼ EQM one-body þ Etwo-body þ Emany-body

ð1Þ

Received Date: October 7, 2010 Accepted Date: November 24, 2010 Published on Web Date: December 03, 2010

3480

DOI: 10.1021/jz101383z |J. Phys. Chem. Lett. 2010, 1, 3480–3487

pubs.acs.org/JPCL

The polarizable force field captures many-body induction selfconsistently, which is particularly important when polar molecules or hydrogen bonding are involved. The use of a force field to capture nonadditive interactions distinguishes15,16 this approach from other fragment methods that use electrostatic embedding.19,32,33 One can either use an existing force field15 or predict the necessary force field parameters in ab initio fashion on the fly.16 The force field used here neglects many-body dispersion effects. Three-body dispersion is usually small, but it can be important in systems with nonpolar molecules, such as in benzene or rare-gas crystals. Several features make a periodic-boundary-conditions (PBC) implementation of HMBI appealing. First, the slow 1/r decay of electrostatic interactions makes traditional periodic electronic structure methods computationally expensive. In periodic HMBI, an inexpensive classical approximation replaces the expensive quantum mechanical treatment of longrange electrostatics. Periodic HMBI is analogous to the incremental approach used in refs 23-26, except the periodic HF/DFT calculation is replaced with an inexpensive polarizable force field one. Second, because the long-range lattice sums are performed classically, it is trivial to apply any standard nonperiodic electronic structure method (MP2, CCSD(T), etc.) to capture the more difficult shorter-range interactions. Third, many fragment methods use embedding charges/densities, which means that analytical derivatives for each monomer and dimer actually depend on the coordinates of all other monomers. HMBI requires no such embedding, so analytical derivatives can be evaluated much less expensively, especially if geometry-independent force-field parameters are used, as they are here. Before presenting any results, we outline the implementation of periodic HMBI. Periodic HMBI treats all central unit-cell monomers and short-range two-body interactions (with either one or both monomers inside the central unit cell) quantum mechanically, while the long-range two-body electrostatics/induction and all many-body induction effects are handled classically (see Figure 1): X QM X damp MM ðEi - EMM fij ðΔ2 EQM EHMBI ij PBC ¼ EPBC þ i Þþ i

- Δ2 EMM ij Þ þ

1 2

Figure 1. For each monomer in the central unit cell, periodic HMBI treats the short-range interactions quantum mechanically (QM) and the long-range ones classically (MM). A damping region (in purple) is used to transition smoothly between the two regimes.

Intramolecular force-field parameters do not contribute. It is possible, for example, to determine the electrostatic forcefield parameters (distributed multipole moments and polarizabilities) from QM,16 but here we simply use the existing Amoeba force field.34 The transition between short-range quantum and longrange classical treatments requires some attention. To ensure smooth and continuous potential energy surfaces, the interactions in this transition region are damped using a function that smoothly decays from 1 at radius r1 to 0 at radius r0,35 1 damp fij ðRÞ ¼ ð3Þ 2jr r j=ðr RÞ - jr1 - r0 j=ðR - r0 Þ 1 0 1 1þe where R is the shortest intermolecular distance between any two atoms in the pair of molecules i and j. For any dimer where the shortest intermolecular separation is less than r1, the two-body interaction is treated quantum mechanically. If the shortest intermolecular distance is greater than r0, it is approximated classically. In the damping region between r1 and r0, the dimer interaction energy is a linear combination of quantum and classical interactions. To demonstrate the usefulness of the periodic HMBI approach, we have predicted the lattice energies of five different molecular crystals: hexagonal ice, formamide, acetamide, imidazole, and benzene. The first three crystals are held together by hydrogen bonds, while π-stacking-type dispersion interactions are important in benzene. The imidazole crystal exhibits both π-stacking and hydrogen-bonding interactions. The source publication or Cambridge Crystallographic Data Center (CCDC) reference code for the initial structures and the lattice parameters are given in Table 1. All atoms were relaxed with fixed lattice constants. These particular crystals were selected because (1) they exhibit a representative range of packing forces, (2) their experimental heats of sublimation/ lattice energies are known, (3) and Amoeba force-field parameters already exist for the monomers. Experimental and predicted lattice energies of the optimized crystal structures are presented in Table 2. Experimental lattice energies are obtained from the measured heats

ij

X images X i

k~

damp

fi~k

ðΔ2 EQM - Δ2 EMM Þ i~k i~k

ð2Þ

Here, i and j run over molecules in the central unit cell, while ~k runs over all periodic image molecules within some cutoff distance of molecule i. Ei corresponds to the energy of monomer i, and Δ2Eij is the interaction energy between monomers i and j. Both can be calculated either quantum mechanically (QM) or with a force field (MM). The full crystal lattice energy, damp EMM PBC, is computed classically via Ewald summation. The fij term is discussed below. Analytical nuclear gradients of the HMBI energy expression are straightforward and have been implemented for geometry-independent force-field parameters and fixed lattice parameters. The force field is used here to capture many-body induction and long-range two-body electrostatics and dispersion.


3481


pubs.acs.org/JPCL

Table 1. Crystal Calculation Properties and HMBI Timings name molecular formula

ice H2O

formamide HCONH2

acetamide CH3CONH2

imidazole C3N2H4

benzene C6H6

unit cell/CCDC code

ref 36.

FORMAM02

ACEMID05

IMAZOL06

BENZEN01

# of molecules in the unit cell a

16

4

18

4

4

QM cutoff r0b

7Å

10 Å

10 Å

10 Å

10 Å

# of dimer calcs for sp energy a

471

264

1008

214

160

DB-RI-MP2/aug-cc-pVTZ timing c

8 min

1.8 h

28 h

11 h

23 h

DB-RI-MP2/aug-cc-pVQZ timing c

40 min

14 h

245 h

84 h

177 h

a

These numbers could be reduced using space group symmetry. b r0 = r1 - 1.0 Å. c Using six 2.5 GHz Intel Xeon E5420 quad cores (24 cores total).

Table 2. HMBI-Predicted Crystal Lattice Energies (kJ/mol) QM levela

ice

formamide

b

force-field contribution

13.9

-5.4

DB-RI-MP2/aug-cc-pVDZ

53.2

67.5

DB-RI-MP2/aug-cc-pVTZ DB-RI-MP2/aug-cc-pVQZ

57.1 58.7

72.3 74.1

DB-RI-MP2/CBS

60.3 1

error vs expt.

acetamide

imidazole

benzene

3.0

4.9

-2.0

74.1

95.3

62.2

78.5 80.3

99.1 99.7

62.1 62.3

76.0

81.7

101.7

63.1

-6

-4

15

11

ΔCCSD(T)c

0.4

1.8

-0.1

-14.2

-10.4

DB-RI-MP2/CBS þΔCCSD(T)

60.6

77.8

81.6

87.5

52.7

2

-4

-4

-4

1

error vs expt. est. lattice param. relax., ΔErelax lattice

0.0

0.1

n/a

0.3

3.7

best estimated

60.6

77.9

81.6

87.8

56.4

error vs expt.

2

-4

-4

-3

4

59

82 ( 0.3

86 ( 2

experimente

91 ( 4

52 ( 3

Counterpoise corrected. Classical polarization contribution. It is already included in the reported lattice energies. Post-MP2 correction, Δ = d DB-RI-MP2/CBS e - EMP2 þ ΔCCSD(T) þ ΔErelax ECCSD(T) lattice lattice, using the basis sets described in the text. Best estimate = E lattice. Reported errors are the standard deviation among the set of extrapolated 0 K lattice energies. Actual experimental errors may be larger. See Supporting Information. a

b

c

of sublimation. For ice, we used the values reported by Whalley.37,38 For each of the other crystals, we started with the set of experimental heat of sublimation measurements collected by Chickos and Acree.39 Each heat of sublimation was extrapolated to 0 K using the approximate relationship, Ulattice = ΔHsub(T) þ 2RT,40 and we computed the median and standard deviation of the 0 K values. For each crystal, we then removed the zero-point energy30,41-46 from the median 0 K value to obtain the final lattice energies. Table 2 lists the final lattice energies and the standard deviation among the various extrapolated 0 K lattice energies. Overall, these experimental lattice energies are probably accurate to within a few kilojoules per mole. See the Supporting Information for more details on these calculations. For the three hydrogen-bonded crystals (ice, formamide, and acetamide), MP2 predicts lattice energies fairly accurately, particularly if triple-ζ or larger basis sets are used. Basis-set extrapolation reduces the errors to only ∼1-6 kJ/mol relative to the experimental values. For benzene and imidazole, π-stacking dispersion interactions are important, and complete-basis-set MP2 unsurprisingly overestimates the


CCSD(T)

lattice energies by 11-15 kJ/mol. MP2 is known to overestimate the binding energy of the benzene dimer and many other π-stacking systems.47,48 The HMBI model here also neglects many-body intermolecular dispersion terms, which can be important.49 Next, we apply a standard additive correction for post-MP2 correlation, ΔCCSD(T) = ECCSD(T) - EMP2, which is calculated in a computationally affordable basis set for each system: augcc-pVTZ for ice, aug-cc-pVDZ for formamide and imidazole, and 6-31þG* for acetamide and benzene. The ΔCCSD(T) corrections are small (less than 2 kJ/mol) for ice, formamide, and acetamide. For benzene and imidazole, however, including post-MP2 correlation drastically improves the lattice energies, lowering the errors to less than 4 kJ/mol. Finally, we estimate the errors arising from the use of fixed experimental lattice parameters during the geometry optimizations. We adopted a procedure similar to the one used in ref 11: we uniformly scaled the lattice constants in increments of 1% and optimized the geometry for each fixed lattice using periodic DFT. We then performed single-point HMBI energy calculations at the dual-basis RI-MP2/aug-cc-pVTZ level along

3482


pubs.acs.org/JPCL

these one-dimensional potential energy curves. The estimated lattice energy increase due to lattice parameter relaxation, ΔErelax lattice, is given in Table 2. Combining this energy with the DB-RI-MP2/CBS þΔCCSD(T) energies gives our best estimate for the lattice energy. Further details are provided in the Supporting Information. For ice, formamide, and imidazole, these estimates indicate that relaxing the lattice parameters changes the lattice energies only slightly (0.3 kJ/mol or less). Very similar results have been found previously for ice.23 Analogous calculations on acetamide proved computationally impractical, but acetamide probably behaves similarly to the other hydrogenbonded crystals. For benzene, however, lattice-parameter relaxation increases the lattice energy by a much larger 3.7 kJ/mol. This raises the error in our best estimate for benzene to 4 kJ/mol, which is comparable to the errors for the other crystals. Overall, the HMBI predictions are within 2-4 kJ/mol of the experimental lattice energies. Three additional main sources of error affect these lattice energy predictions. First, our previous gas-phase work has shown that the Amoeba force field overestimates many-body effects in water and underestimates them in formamide.16 We see the same pattern here in the predicted lattice energies. Similar considerations probably apply to the other crystals. Improved many-body induction parameters can be predicted quantum mechanically.16 Second, the Amoeba force field neglects many-body dispersion interactions. For benzene, these contribute -6.5 kJ/mol.30 They are also probably important in imidazole. In principle, inexpensive estimates for three-body dispersion could be incorporated into the HMBI many-body force field to help capture these effects.50 Third, practical computational considerations limited the size of the basis set used in our CCSD(T) calculations. Comparing our data against aug-cc-pVDZ results in refs 25 and 51, we estimate that going from 6 to 31þG* to aug-cc-pVDZ would increase the predicted benzene lattice energy by several kilojoules per mole. In other words, the neglect of many-body dispersion in benzene is largely canceled by basis set errors in our ΔCCSD(T) correction. We now compare these HMBI results to previously published results obtained with other techniques (see Table 3). Benchmark, large-basis periodic MP2 calculations on the NH3 and CO2 crystals predict lattice energies that differ from experiment by only 1-2 kJ/mol,11 which is comparable to what HMBI predicts with MP2 for ice. Triple-ζ local MP2 calculations gave errors of 1-6 on a set of six different crystals.14 HMBI is less computationally expensive than periodic MP2. As mentioned previously, local periodic MP2 is practical in systems with up to ∼50 atoms in the unit cell in a medium-sized basis.11 For comparison, the acetamide unit cell considered here contains 162 atoms (72 heavy atoms and 90 hydrogen atoms). In addition, HMBI and other fragment methods enable the straightforward inclusion of post-MP2 correlation effects when needed (as for benzene and imidazole). The results reported here also compare well with previous high-level calculations on ice and benzene. Periodic HF/6311G** combined with one- and two-body MP2/aug-cc-pVTZ interaction energies predicts the lattice energy of ice to be 55.7 kJ/mol.23 HMBI similarly predicts 57.1 kJ/mol in the same


basis, despite replacing periodic HF with a polarizable force field estimate of electrostatic induction. Periodic HF calculations can be difficult to converge in large Gaussian basis sets, and periodic HMBI avoids such difficulties. A number of previous high-level calculations have been performed on the benzene crystal as well. Ringer and Sherrill used MP2 and CCSD(T) with a truncated many-body expansion and the experimental geometry to predict a value of 56.4 kJ/mol.51 Compared to this calculation, HMBI includes additional long-range dimer interactions and classical many-body induction contributions. A similar calculation that also used DFT to capture long-range and many-body effects and fixedcell optimization predicted a lattice energy of 50.5 kJ/mol, which is fairly close to our fixed-cell value of 52.7 kJ/mol.27 Podeszwa and co-workers used SAPT(DFT) to predict twoand three-body interactions in the experimental benzene crystal structure and obtained a lattice energy of 50.3 kJ/ mol.30 Their calculations include a three-body dispersion contribution of -6.5 kJ/mol, which is missing from the HMBI calculations here. As noted above, however, this error is largely canceled by the small-basis-set error in our ΔCCSD(T) correction. Lu and co-workers used their adiabatic-connection DFT approach and a fully relaxed structure to obtain a lattice energy of 47.0 kJ/mol.8,9 All of these predictions lie within a few kilojoules per mole of each other, and the experimental errors are probably similar in magnitude. Finally, we compare these results to those from traditional and empirically corrected DFT. Fully relaxed planewave DFT calculations predict the lattice energy of ice to be 109.9 (LDA), 67.7 (PW91), or 64.6 (PBE) kJ/mol.53 All three functionals overbind the ice crystal, although PBE and PW91 perform much more reasonably than LDA. Many generalized-gradient approximation functionals perform well for water-water interactions in the gas and liquid phases, so the decent performance of PBE and PW91 here is unsurprising. For four of the crystals examined here, Feng and Li performed periodic B3LYP/6-21G** calculations in which they optimized the structure with a fixed unit cell and included empirical dispersion corrections.52 B3LYP itself fails to capture dispersion interactions54 and underestimates the lattice energy by ∼50-75%. Adding an empirical dispersion correction drastically improves the B3LYP results in all cases. It should be noted that the 6-21G* basis set is quite small. Increasing the basis to 6-311G** raises the benzene lattice energy by 2.6 to 52.4 kJ/mol,52 which is very close to experiment and our predicted value of 52.7 kJ/mol. Of course, simple empirical DFT dispersion corrections of the sort used here also neglect three-body dispersion terms. A similar periodic B3LYP-D/6-31G* (using Grimme's empirical dispersion correction55) study also considered the benzene and formamide crystals.4 Their lattice energies without the dispersion correction and with fully relaxed structures are similar to those predicted by Feng and Li, but the corrected values differ more substantially. B3LYP-D/6-31G* gives lattice energies of -58.3 (benzene) and -84.0 (formamide) kJ/mol, both of which are more than 8 kJ/mol larger than those reported in ref 52. On one hand, periodic DFT performs well with moderate basis sets and low cost. HMBI and other fragment methods

3483


pubs.acs.org/JPCL

Table 3. Summary of Recently Published Lattice Energy Predictions for These Molecular Crystals (kJ/mol)a source

ice

formamide

acetamide

imidazole

benzene

Experimental Crystal Geometries: SAPT(DFT)

b

50.3

CCSD(T) 1 þ 2 bodyc

56.4

Relaxed Cell with Experimental Lattice Parameters: this work, DB-RI-MP2/CBS þΔ

CCSD(T)

60.6

77.8

81.6

87.5

CCSD(T) 1 þ 2 body þ PBE many-bodyd

52.7 50.5

B3LYP/6-21G* e

42.1

40.1

39.6

13.8

B3LYP/6-21G* þ dispersion e

75.9

84.8

93.0

49.8

81.6

87.8

56.4

Fully Relaxed Geometries: this work, best estimate

60.6

MP2/aug-cc-pVTZ þ HF/6-311G** f

55.7

ACDFT g LDA h

109.9

PW91 h

67.7

PBE h

64.6

77.9

47.0

B3LYP/6-31G* i

43.2

15.9

B3LYP-D/6-31G* i

84.0

58.3

experiment

59

82 ( 0.3

86 ( 2

91 ( 4

52 ( 3

a

The results are grouped based on the degree of crystal structure relaxation performed in obtaining the lattice energies. b Reference 30. c Reference 51. d Reference 25. e Reference 52. f Reference 23. g Reference 8. h Reference 53. i Reference 4.

make wave function models affordable in large systems, but large-basis-set MP2 and CCSD(T) are typically much more expensive than DFT. On the other hand, empirically corrected DFT results clearly depend on the particular details of the dispersion correction.4,52 A recent quantum Monte Carlo benchmark study also found that a widely used empirical dispersion correction scheme substantially overestimated the energy difference between two polymorphs of diiodobenzene crystal.6 While periodic DFT is cheaper than many of the HMBI calculations reported here, the HMBI computational requirements are reasonable on modern multiprocessor clusters, as shown in Table 1. The HMBI algorithm is embarrassingly parallel: each of the time-consuming quantum mechanical two-body (dimer) calculations can be performed on a separate processor (or processor group). The MP2-level ice calculations require only minutes for a single-point energy on 24 processors. The MP2 costs are significantly higher, but still affordable, for the other crystals. Unsurprisingly, even the small-basis CCSD(T) calculations here require 3-10 times more computer time than the aug-cc-pVQZ RI-MP2 calculations, depending on the specific molecules and basis sets. Two factors affect the cost: unit-cell size and the size of the individual monomers. The acetamide unit cell contains 18 molecules, which means that 1008 dimer interactions must be computed using our conservative cutoffs. Fortunately, the HMBI approach is linear-scaling with respect to the number of monomers in the unit-cell size: each central-unit cell molecule


interacts with all central-cell and periodic-image molecules within the cutoff radius r0. Doubling the number of monomers in the unit-cell simply doubles the number of central-unit cell monomers whose interactions need to be evaluated. The monomer size affects the computational cost for each dimer interaction. At the MP2 level (N5 scaling), the 1008 dimer calculations in acetamide end up requiring the most computational effort of the five crystals considered here. In contrast, the 160 dimer calculations in benzene crystal become more expensive than the 1008 acetamide dimer ones at the CCSD(T) level due to the N7 scaling of each dimer calculation. Further computational savings could be obtained by taking advantage of crystal space-group symmetry and by adjusting the cutoff at which the model transitions to purely classical electrostatics. Here, only simple translational symmetry was used, and a conservative cutoff of 10 Å was used for all crystals except ice. For ice, we examined more aggressive cutoffs and found that the lattice energy is little changed even at a 7 Å cutoff (see Supporting Information). Shrinking the cutoff radius reduced the number of required dimer calculations by a factor of 3. More aggressive cutoffs for the other species may be possible, although we have not investigated this aspect. In summary, we have extended the hybrid quantum/classical HMBI model to infinite systems using periodic boundary conditions. We demonstrated that it allows one to compute the lattice energy of small-molecule organic crystals within several kJ/mol (∼1 kcal/mol). Such accuracy requires a complete-basis-set MP2-type treatment of two-body intermolecular interactions,

3484


pubs.acs.org/JPCL

and sometimes a small-basis estimate for post-MP2 correlation. These HMBI calculations compare well with other highlevel benchmark calculations. Empirically corrected periodic DFT also predicts reasonably accurate lattice energies, although it appears sensitive to the details of the dispersion correction and may have difficulties for predicted energy differences between crystal polymorphs.6 In this context, the HMBI approach provides a robust and systematically improvable method for predicting crystal lattice energies. A few limitations hinder our current approach. First, the Amoeba force field predicts only moderately accurate induction energies, and it completely neglects many-body dispersion interactions. We are currently adapting a much more accurate ab initio electrostatic induction force field16 to periodic systems and exploring possibilities for including three-body dispersion terms. Second, although the HMBI approach reduces the problem of a molecular crystal to a series of monomer and dimer calculations, even dimer calculations remain computationally expensive with high-level electronic structure methods. Further developments in accurate, inexpensive electron-electron correlation methods would be especially valuable. The ability of our approach to predict accurate crystal structures in the absence experimental structural data also remains to be seen. In any case, as the results presented here demonstrate, fragment-based approaches are particularly promising for applications in molecular crystal structure prediction and other molecular condensed-phase systems.

To estimate the impact of relaxing the lattice parameters from their experimental values, a one-dimensional potential energy scan was performed by uniformly scaling the lattice parameters a, b, and c in increments of 1%. For each set of lattice parameters, the structure was optimized with planewave DFT using the PBE functional, ultrasoft pseudopotentials, and a 340 eV planewave cutoff, as implemented in the Dacapo software package.69 Single-point HMBI DB-RI-MP2/ aug-cc-pVTZ calculations were performed for each lattice parameter set, and cubic splines were used to identify the optimal parameters and change in the lattice energy. These calculations were performed for all crystals except acetamide, where the planewave DFT calculations proved expensive and poorly convergent.

COMPUTATIONAL METHODS

ACKNOWLEDGMENT Computer time for this research was provided by the National Science Foundation Teragrid (TG-CHE090099).

SUPPORTING INFORMATION AVAILABLE Further details regarding the experimental lattice energies, the crystal structure optimization, the calculated/extrapolated total energies for each crystal, and the quantum/classical cutoff parameters. This material is available free of charge via the Internet at http://pubs. acs.org/.

AUTHOR INFORMATION Corresponding Author: *To whom correspondence should be addressed. E-mail: [email protected].

Quantum mechanical energies were evaluated using frozen core, resolution of the identity MP2 (RI-MP2),56 dualbasis57 RI-MP2 (DB-RI-MP2) or standard frozen core CCSD(T).58,59 DB-RI-MP2 achieves nearly canonical MP2 accuracy at vastly lower cost. All quantum mechanical calculations were performed with the Q-Chem60 (MP2) or PSI361 (CCSD(T)) software packages. The classical calculations were performed using the Amoeba force field34 as implemented in the TINKER software package.62 Tinfoil boundary conditions were used in the Amoeba Ewald summation. The value of the cutoff distance r0 beyond which the interactions are treated purely classically is given in Table 1. The inner edge of the damping region is given by r1=r0 - 1.0 Å, which performs well in empirical testing. Crystal geometries were optimized at the RI-MP2/aug-ccpVDZ:Amoeba level using the DL-FIND63 open-source geometry optimizer. For ice, the geometry was refined further in the aug-cc-pVTZ basis set. In all cases, crystal geometries were relaxed with fixed experimental lattice parameters. In computing the lattice energies, counterpoise (CP)-corrected64 DB-RI-MP2 dimer energies were calculated in each of the aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets65 and the corresponding RI auxiliary66 basis sets. The triple- and quadruple-ζ basis set results were used to extrapolate the CP-corrected HF and correlation energies to the completebasis-set limit.67,68 Lattice energies were computed relative to monomer energies that were optimized at the same level of theory. Post-MP2 ΔCCSD(T) corrections were calculated as the difference between CCSD(T) and canonical MP2 (both with frozen core) in the basis set described in the main text.


REFERENCES (1)

(2)

(3)

(4)

(5)

(6)

(7)

3485

Price, S. L. Computed Crystal Energy Landscapes for Understanding and Predicting Organic Crystal Structures and Polymorphism. Acc. Chem. Res. 2009, 42, 117–126. Li, T.; Feng, S. Empirically Augmented Density Functional Theory for Predicting Lattice Energies of Aspirin, Acetaminophen Polymorphs, and Ibuprofen Homochiral and Racemic Crystals. Pharm. Res. 2006, 23, 2326–2332. Neumann, M. A.; Leusen, F. J. J.; Kendrick, J. A Major Advance in Crystal Structure Prediction. Angew. Chem., Int. Ed. 2008, 47, 2427–2430. Civalleri, B.; Zicovich-Wilson, C. M.; Valenzano, L.; Ugliengo, P. B3LYP Augmented with an Empirical Dispersion Term (B3LYP-D*) as Applied to Molecular Crystals. CrystEngComm 2008, 10, 405–410. Sorescu, D. C.; Rice, B. M. Theoretical Predictions of Energetic Molecular Crystals at Ambient and Hydrostatic Compression Conditions Using Dispersion Corrections to Conventional Density Functionals (DFT-D). J. Phys. Chem. C 2010, 114, 6734–6748. Hongo, K.; Watson, M. A.; Sanchez-Carrera, R. S.; Iitaka, T.; Aspuru-Guzik, A. Failure of Conventional Density Functionals for the Prediction of Molecular Crystal Polymorphism: A Quantum Monte Carlo Study. J. Phys. Chem. Lett. 2010, 1, 1789–1794. Kleis, J.; Lundqvist, B. I.; Langreth, D. C.; Schr€ oder, E. Towards a Working Density-Functional Theory for Polymers: FirstPrinciples Determination of the Polyethylene Crystal Structure. Phys. Rev. B 2007, 76, 1002001.


pubs.acs.org/JPCL

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

Lu, D.; Li, Y.; Rocca, D.; Galli, G. Ab Initio Calculation of van der Waals Bonded Molecular Crystals. Phys. Rev. Lett. 2009, 102, 206411. Li, Y.; Lu, D.; Nguyen, H.-V.; Galli, G. Van der Waals Interactions in Molecular Assemblies From First-Principles Calculations. J. Phys. Chem. A 2010, 114, 1944–1952. Pisani, C.; Maschio, L.; Casassa, S.; Halo, M.; Sch€ utz, M.; Usvyat, D. Periodic Local MP2 Method for the Study of Electronic Correlation in Crystals: Theory and Preliminary Applications. J. Comput. Chem. 2008, 29, 2113–2124. Maschio, L.; Usvyat, D.; Sch€ utz, M.; Civalleri, B. Periodic Local Møller-Plesset Second Order Perturbation Theory Method Applied to Molecular Crystals: Study of Solid NH3 and CO2 Using Extended Basis Sets. J. Chem. Phys. 2010, 132, 134706. Erba, A.; Pisani, C.; Casassa, S.; Maschio, L.; Sch€ utz, M.; Usvyat, D. MP2 Versus Density-Functional Theory Study of the Compton Profiles of Crystalline Urea. Phys. Rev. B 2010, 81, 165108. Marsman, M.; Grueneis, A.; Paier, J.; Kresse, G. Second-Order Møller-Plesset Perturbation Theory Applied to Extended Systems. I. Within the Projector-Augmented-Wave Formalism Using a Plane Wave Basis Set. J. Chem. Phys. 2009, 130, 184103. Maschio, L.; Usvyaat, D.; Civalleri, B. Ab Initio Study of van der Waals and Hydrogen-Bonded Molecular Crystals with a Periodic Local-MP2 Method. CrystEngComm 2010, 12, 2429– 2435. Beran, G. J. O. Approximating Quantum Many-Body Intermolecular Interactions in Molecular Clusters Using Classical Polarizable Force Fields. J. Chem. Phys. 2009, 130, 164115. Sebetci, A.; Beran, G. J. O. Spatially Homogeneous QM/MM for Systems of Interacting Molecules with On-the-Fly Ab Initio Force-Field Parameterization. J. Chem. Theory Comput. 2010, 6, 155–167. Nagayoshi, K.; Ikeda, T.; Kitaura, K.; Nagase, S. Computational Procedure of Lattice Energy Using the Ab Initio MO Method. J. Theory. Comput. Chem. 2003, 2, 233–244. Hirata, S.; Valiev, M.; Dupuis, M.; Xantheas, S. S.; Sugiki, S.; Sekino, H. Fast Electron Correlation Methods for Molecular Clusters in the Ground and Excited States. Mol. Phys. 2005, 103, 2255–2265. Kamiya, M.; Hirata, S.; Valiev, M. Fast Electron Correlation Methods for Molecular Clusters without Basis Set Superposition Errors. J. Chem. Phys. 2008, 128, 074103. Paulus, B.; Rosciszewski, K.; Gaston, N.; Schwerdtfeger, P.; Stoll, H. Convergence of the Ab Initio Many-Body Expansion for the Cohesive Energy of Solid Mercury. Phys. Rev. B 2004, 70, 165106. Stoll, H.; Paulus, B.; Fulde, P. On the Accuracy of CorrelationEnergy Expansions in Terms of Local Increments. J. Chem. Phys. 2005, 123, 144108. Manby, F. R.; Alfe, D.; Gillan, M. J. Extension of Molecular Electronic Structure Methods to the Solid State: Computation of the Cohesive Energy of Lithium Hydride. Phys. Chem. Chem. Phys. 2006, 8, 5178–5180. Hermann, A.; Schwerdtfeger, P. Ground-State Properties of Crystalline Ice from Periodic Hartree-Fock Calculations and a Coupled-Cluster-Based Many-Body Decomposition of the Correlation Energy. Phys. Rev. Lett. 2008, 101, 183005. Hermann, A.; Schwerdtfeger, P. Complete Basis Set Limit Second-Order Møller-Plesset Calculations for the fcc Lattices of Neon, Argon, Krypton, and Xenon. J. Phys. Chem. 2009, 131, 244508. Bludsky, O.; Rubes, M.; Soldan, P. Ab Initio Investigation of Intermolecular Interactions in Solid Benzene. Phys. Rev. B 2008, 77, 092103.


(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

(34)

(35)

(36)

(37) (38)

(39)

(40)

(41)

(42)

(43)

(44)

3486

Tsuzuki, S.; Orita, H.; Honda, K.; Mikami, M. First-Principles Lattice Energy Calculation of Urea and Hexamine Crystals by a Combination of Periodic DFT and MP2 Two-Body Interaction Energy Calculations. J. Phys. Chem. B 2010, 114, 6799– 6805. Hirata, S. Fast Electron-Correlation Methods for Molecular Crystals: An Application to the R, β1, and β2 Modifications of Solid Formic Acid. J. Chem. Phys. 2008, 129, 204104. Sode, O.; Keceli, M.; Hirata, S.; Yagi, K. Coupled-Cluster and Many-Body Perturbation Study of Energies, Structures, and Phonon Dispersions of Solid Hydrogen Fluoride. Int. J. Quantum Chem. 2009, 109, 1928–1939. Podeszwa, R.; Bukowski, R.; Rice, B. M.; Szalewicz, K. Potential Energy Surface for Cyclotrimethylene Trinitramine Dimer From Symmetry-Adapted Perturbation Theory. Phys. Chem. Chem. Phys. 2007, 9, 5561–5569. Podeszwa, R.; Rice, B. M.; Szalewicz, K. Predicting Structures of Molecular Crystals From First-Principles. Phys. Rev. Lett. 2008, 101, 115503. Misquitta, A. J.; Welch, G. W. A.; Stone, A. J.; Price, S. L. A First Principles Prediction of the Crystal Structure of C6Br2ClFH2. Chem. Phys. Lett. 2008, 456, 105–109. Fedorov, D. G.; Kitaura, K. Extending the Power of Quantum Chemsitry to Large Systems with the Fragment Molecular Orbital Method. J. Phys. Chem. A 2007, 111, 6904–6914. Dahlke, E. E.; Truhlar, D. G. Electrostatically Embedded Many-Body Correlation Energy, with Applications to the Calculation of Accurate Second-Order Møller-Plesset Perturbation Theory Energies for Large Water Clusters. J. Chem. Theory Comput. 2007, 3, 1342–1348. Ren, P.; Ponder, J. W. Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation. J. Phys. Chem. B 2003, 107, 5933–5947. Subotnik, J. E.; Sodt, A.; Head-Gordon, M. The Limits of Local Correlation Theory: Electronic Delocalization and Chemically Smooth Potential Energy Surfaces. J. Chem. Phys. 2008, 128, 034103. Morrison, I.; Li, J.-C.; Jenkins, S.; Xantheas, S. S.; Payne, M. C. Ab Initio Total Energy Studies of the Static and Dynamical Properties of Ice Ih. J. Phys. Chem. B 1997, 101, 6146–6150. Whalley, E. The Difference in the Intermolecular Forces of H2O and D2O. Trans. Faraday Soc. 1957, 53, 1578–1585. Whalley, E. The Hydrogen Bond in Ice. In The Hydrogen Bond; Schuster, P., Zundel, G., Sandorfy, C., Eds.; North-Holland: Amsterdam, 1976; Vol. 3, Chapter 29, pp 1425-1470. Chickos, J. S.; Acree, W. E. Enthalpies of Sublimation of Organic and Organometallic Compounds. J. Phys. Chem. Ref. Data 2002, 31, 537–698. deWit, H. G. M.; van Miltenburg, J. C.; de Kruif, C. G. Thermodynamic Properties of Molecular Organic Crystals Containing Nitrogen, Oxygen, and Sulphur. 1. Vapour Pressures and Enthalpies of Sublimation. J. Chem. Thermodyn. 1983, 15, 651–663. Tam, C. N.; Bour, P.; Eckert, J.; Trouw, F. R. Inelastic Neutron Scattering Study of Hydrogen-Bonded Solid Formamide at 15 K. J. Phys. Chem. A 1997, 101, 5877–5884. King, S. T. Infrared Study of the NH2 “Inversion” Vibration for Formamide in the Vapor Phase and in an Argon Matrix. J. Phys. Chem. 1971, 75, 405–410. Day, G. M.; Price, S. L.; Leslie, M. Atomistic Calculations of Phonon Frequencies and Thermodynamic Quantities for Crystals of Rigid Organic Molecules. J. Phys. Chem. B 2003, 107, 10919–10933. Ganeshsrinivas, E.; Sathyanarayana, D. N.; Machida, K.; Miwa, Y. Simulation of the Infrared Spectra of Acetamide


pubs.acs.org/JPCL

(45)

(46)

(47)

(48)

(49)

(50)

(51)

(52)

(53)

(54)

(55)

(56)

(57)

(58)

(59)

(60)

(61)

(62)

by an Extended Molecular Mechanics Method. J. Mol. Struct. (THEOCHEM) 1996, 361, 217–227. Uni, T.; Machida, K.; Saito, Y. Infrared Spectra of Partially Deuterated Acetamide. Bull. Chem. Soc. Jpn. 1969, 42, 897– 904. Uni, T.; Machida, K.; Saito, Y. Out-of-Plane Vibrations of Acetamide and Partially N-Deuterated Acetamide. Spectrochim. Acta A 1971, 27, 833–844. Sinnokrot, M.; Sherrill, C. D. High-Accuracy Quantum Mechanical Studies of π-π Interactions in Benzene Dimers. J. Phys. Chem. A 2006, 110, 10656–10668. y, J.; Hobza, P. Benchmark DataJurecka, P.; Sponer, J.; Cern base of Accurate (MP2 and CCSD(T) Complete Basis Set Limit) Interaction Energies of Small Model Complexes, DNA Base Pairs, and Amino Acid Pairs. Phys. Chem. Chem. Phys. 2006, 8, 1985–1993. Podeszwa, R. Comment on “Beyond the Benzene Dimer: An Investigation of the Additivity of π-π Interactions. J. Phys. Chem. A 2008, 112, 8884–8885. von Lilienfeld, O. A.; Tkatchenko, A. Two- and Three-Body Interatomic Dispersion Energy Contributions to Binding in Molecules and Solids. J. Chem. Phys. 2010, 132, 234109. Ringer, A. L.; Sherrill, C. D. First-Principles Computation of Lattice Energies of Organic Solids: The Benzene Crystal. Chem.;Eur. J. 2008, 14, 2542–2547. Feng, S.; Li, T. Predicting Lattice Energy of Organic Crystals by Density Functional Theory with Empirically Corrected Dispersion Energy. J. Chem. Theory Comput. 2006, 2, 149–156. Thierfelder, C.; Hermann, A.; Schwerdtfeger, P.; Schmidt, W. G. Strongly Bonded Water Monomers on the Ice Ih Basal Plane: Density-Functional Calculations. Phys. Rev. B 2006, 74, 045422. Kristyan, S.; Pulay, P. Can (Semi)Local Density Functional Theory Account for the London Dispersion Forces? Chem. Phys. Lett. 1994, 229, 175-180. Grimme, S. Accurate Description of van derWaals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463–1473. Feyereisen, M. W.; Fitzgerald, G.; Komornicki, A. Use of Approximate Integrals in Ab Initio Theory. An Application in MP2 Energy Calculations. Chem. Phys. Lett. 1993, 208, 359–363. Steele, R. P.; Distasio, R. A., Jr.; Shao, Y.; Kong, J.; HeadGordon, M. Dual-Basis Second-Order Møller-Plesset Perturbation Theory: A Reduced Cost Reference for Correlation Calculations. J. Chem. Phys. 2006, 125, 074108. Purvis, G. D., III; Bartlett, R. J. A Full Coupled-Cluster Singles and Doubles Model: The Inclusion of Disconnected Triples. J. Chem. Phys. 1982, 76, 1910–1918. Raghavachari, K.; Trucks, G.; Pople, J. A.; Head-Gordon, M. A Fifth-Order Perturbation Comparison of Electron Correlation Theories. Chem. Phys. Lett. 1989, 157, 479–483. Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown, S. T.; Gilbert, A. T. B.; Slipchenko, L. V.; Levchenko, S. V.; O'Neill, D. P.; et al. Advances in Methods and Algorithms in a Modern Quantum Chemistry Program Package. Phys. Chem. Chem. Phys. 2006, 8, 3172–3191. Crawford, T. D.; Sherrill, C. D.; Valeev, E. F.; Fermann, J. T.; King, R. A.; Leininger, M. L.; Brown, S. T.; Janssen, C. L.; Seidl, E. T.; Kenny, J. P.; et al. PSI3: An Open-Source Ab Initio Electronic Structure Package. J. Comput. Chem. 2007, 28, 1610–1616. Ponder, J.W. TINKER, v4.2, 2004, http://dasher.wustl.edu/ tinker/. Accessed January 23, 2008.


(63)

(64)

(65)

(66)

(67)

(68)

(69)

3487

K€ astner, J.; Carr, J. M.; Keal, T. W.; Thiel, W.; Wander, A.; Sherwood, P. DL-FIND: An Open-Source Geometry Optimizer for Atomistic Simulations. J. Phys. Chem. A 2009, 113, 11856– 11865. Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553–566. Dunning, T. H., Jr. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron Through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. Weigend, F.; K€ ohn, A.; H€ attig, C. Efficient Use of Correlation Consistent Basis Sets in the Resolution of the Identity MP2 Calculations. J. Chem. Phys. 2002, 116, 3175–3183. Karton, A.; Martin, J. M. L. Comment on: “Estimating the Hartree-Fock Limit From Finite Basis Set Calculations. Theor. Chem. Acc. 2006, 115, 330–333. Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. Basis-Set Convergence of Correlated Calculations on Water. J. Chem. Phys. 1997, 106, 9639–9646. Dacapo v2.7.7, https://wiki.fysik.dtu.dk/dacapo (accessed May 2, 2006).


Predicting Organic Crystal Lattice Energies with Chemical Accuracy

Recommend Documents