ARTICLE pubs.acs.org/JPCA
Intermolecular Interaction Energies in Molecular Crystals: Comparison and Agreement of Localized MøllerPlesset 2, Dispersion-Corrected Density Functional, and Classical Empirical Two-Body Calculations Lorenzo Maschio,† Bartolomeo Civalleri,† Piero Ugliengo,†,* and Angelo Gavezzotti‡ † ‡
Dipartimento di Chimica IFM and NIS Centre of Excellence, Universita di Torino, via P. Giuria 7, 10125 Torino, Italy Dipartimento di Chimica Strutturale e Stereochimica Inorganica, Universita di Milano, via Venezian 21, 20133 Milano, Italy
bS Supporting Information ABSTRACT: A comparative analysis of the intermolecular energy for a data set including 60 molecular crystals with a large variety of functional groups has been carried out using three different computational approaches: (i) a method based on a physically meaningful empirical partition of the interaction energy (PIXEL), (ii) density functional methods with a posteriori empirical correction for the dispersion interactions (DFT-D), and (iii) a full periodic ab initio quantum mechanical method based on MøllerPlesset perturbation theory for the electron correlation using localized crystal orbitals (LMP2). Due to the large computational cost, LMP2 calculations have been restricted to a subset of seven molecular crystal comprising benzene, formic acid, formamide, succinic anhydride, urea, oxalic acid, and nitroguanidine, and the results compared with PIXEL and DFT-D data as well as with the experimental data show excellent agreement among all adopted methods. This shows that both DFT-D and PIXEL approaches are robust predictive tools for studying molecular crystals. A detailed analysis shows a very similar dispersion contribution of the two methods across the 60 considered molecular crystals. The study also confirms that pure DFT shows serious deficiencies in properly handling molecular crystals in which the dispersive contribution is large. Due to the negligible requested computational resources, PIXEL is the method of choice in screening of a large number of molecular crystals, an essential step to predict crystal polymorphism or to study crystal growth processes. DFT-D can then be used to refine the ranking emerged from PIXEL calculations due to its general applicability and robustness in properly handling short-range interactions.
’ INTRODUCTION A proper quantitative evaluation of intermolecular interaction energies is a prerequisite for an understanding of molecular aggregation modes in condensed phases, and hence for prediction and control of structural, thermodynamic, and physical properties of materials. Recently, quantum chemical methods have advanced to the point that cohesion energies of elementary molecular aggregates are more easily calculated than measured,17 and, in general, molecular modeling is nowadays accepted and even widely encouraged as a complement, and in some cases even a substitute, of expensive experimental investigation. However, the range of available computational methods, and of corresponding economical effort either in hardware or in man-time and machine software, is very wide, and researchers of moderate experience in chemical theory are often perplexed when it comes to a proper selection. The optimization of cost-effectiveness ratios for a widespread use of such methods is of particular concern. The present contribution compares the performance of several methods to deal with molecular crystals, focusing on molecular sizes from 30 to 200 Da or from 10 to 100 valence electrons, a typical range of practical relevance for materials science or biological and pharmaceutical chemistry. The methods chosen for the comparisons are (a) LMP2, the latest implementation of MP2 (second order MøllerPlesset perturbation theory, for the first time applied to molecular H2 by r 2011 American Chemical Society
Schulman and Kaufman8 in 1970 and subsequently engineered in Gaussian70 by Pople and co-workers9) recently extended to periodic systems by exploiting the locality in the electron correlation;10,11 (b) uncorrected density functional theory, DFT-B3LYP; (c) DFT-D,12 or dispersion-corrected DFT, in its B3LYP-D* form13,14 (for the latter two methods, the program package CRYSTAL0915,16 has been used); (d) canonical MP2/6-31G(d,p) molecular orbital; and (e) PIXEL, a hybrid method based on numerical integration of classical formulas over quantum chemical electron densities, in its latest formulation.17 At variance with many benchmark studies focusing exclusively on isolated molecules or small aggregates, the chosen model systems are organic crystals, considering both full lattice energies compared with experimental sublimation enthalpies, and separate interactions between molecular pairs in the crystal coordination sphere. These pairwise energies convey important information, because they reveal the hierarchy of forces at work in a crystal structure. Moreover, the PIXEL partitioning of cohesion energies into Coulomb-polarization and dispersion terms reveals the nature of Special Issue: Pavel Hobza Festschrift Received: April 5, 2011 Revised: August 23, 2011 Published: September 06, 2011 11179
dx.doi.org/10.1021/jp203132k | J. Phys. Chem. A 2011, 115, 11179–11186
The Journal of Physical Chemistry A the interaction. Because of the availability of unquestionable structural data from X-ray diffraction, crystals have recently been a primary tool and target in the optimization of quantum chemical methods for intermolecular interactions.1825 LMP2 calculations are, for the first time, presented for crystals of sizable molecules, and together with (and perhaps even better than) experimental enthalpies of sublimation, they provide a valuable cornerstone matching point. Our present results show that the two semiempirical methods, B3LYP-D* and PIXEL, give reliable information on intermolecular interaction energies at moderate implementation and running costs. Moreover, the excellent agreement between the separate estimation of dispersion energies by the two methods shows that this crucial contribution is consistently defined and well accounted for.
ARTICLE
Table 1. Calculated Lattice Energies (kJ mol1)a
benzene
LMP2
B3LYP-D*
B97-D
PIXEL
exp. ΔH(s)b
56.6
48.1
56.7
49.9
45c
formic acid
63.2
70.6
65.1
56.2
6468
formamide
75.9
84.1
79.8
73.0
72
succinic anhydride
87.0
89.6
86.7
77.5
81
urea
108.6
119.5
115.2
107.8
8899
oxalic acid
119.7
128.5
119.7
106.3
9398
nitroguanidine
156.7
150.2
143.7
157.8
143
a
Method/basis sets: B3LYP-D*/6-31G(d,p), B97-D/6-31G(d,p); LMP2/p-aug-631G(d,p). b See ref 27; see also the NIST Webbook at (http://webbook/nist/gov/chemistry) and ref 48, pp 190ff; for oxalic acid, see also ref 60. c Up to 53.9 at 153 K.
2. COMPUTATIONAL DETAILS 2.1. Crystal Structure Retrieval. For crystals, atomic posi-
tions for non-hydrogen atoms are retrieved from the singlecrystal X-ray data in the Cambridge Structural Database (CSD).26 Unreliable X-ray hydrogen-atom positions are renormalized by geometrical procedures as usual17 to CH, NH, and OH distances of 1.08, 1.00, and 1.00 Å, respectively, with bond and torsion angles derived by educated guesses starting from the X-ray data. This procedure has been shown by long experience to be preferable to reoptimization of the whole molecular structure, which leads to undesirable modifications of the whole molecular geometry including non-hydrogen atoms. The choice of structures included in the database represents a wide range of chemical functionalities, including aliphatic and aromatic hydrocarbons, esters and quinones, nitro compounds, aromatic nitrogencompounds and nitriles, chlorinated aromatics, carboxylic acids, amines and amides, and azoles (fluorine compounds represent a special problem and will be the subject of future work). The choice was dictated primarily by the simultaneous availability of an experimental sublimation enthalpy,27 and, second, by the requirement of as rigid molecular conformation as possible (see below). Our database is of unprecedented size, containing 60 crystal structures; the full list with CSD refcodes and experimental data is available as Supporting Information. Pairs of close neighbor molecules were extracted at fixed geometry from the crystal structures, and their interaction energies were calculated as for rigid gas-phase dimers by B3LYP-D*, PIXEL, and canonical MP2 methods. 2.2. Periodic Ab Initio DFT-D (B3LYP-D*) and LMP2 Calculations. B3LYP-D*13 calculations were carried out with the periodic ab initio code CRYSTAL09.16 Periodic LMP2 calculations were performed with a development version of the CRYSCOR09 code.10,11 A 6-31G(d,p) basis set was employed for all B3LYP-D* calculations. The 6-31G(d,p) basis set was augmented with low-exponent polarization functions (d- for C, N and O atoms, p- for H atoms, taken from the aug-cc-pVDZ basis set) to compute the LMP2 post-HF correction through a dual basis set technique. This technique implies the evaluation of the HF reference energy with a smaller basis, while the subsequent calculation of the LMP2 correlation energy employs the virtual space represented in a larger basis. The correction to the HF energy is evaluated using the contribution of first-order singles in the larger basis. The augmented basis set were denoted as p-aug-6-31G(d,p). Parameters for the HF calculations are as described in our previous work1316,28 with CRYSTAL tolerances ITOL13 set to 7 and ITOL45 set to 15 and 50, respectively (see ref 16 for
further details). DFT and MP2 energies are corrected for basis set superposition error by the counterpoise method.29 For the case of isolated dimers canonical MP2 calculations have been carried out using GAUSSIAN0330 computer code adopting the frozen core approximation. 2.2.1. DFT-D Correction in CRYSTAL. The DFT-D method12 uses a damped pairwise empirical potential to include long-range dispersion contributions to the computed DFT total energy and gradients (i.e., the B3LYP functional in the present work). The total energy is given by EðDFT-DÞ ¼ EðDFTÞ þ EðDÞ EðDÞ ¼ s6
∑g ∑ij
0
ij
fdmp ðRij, g Þ
C6 Rij,6 g
in which E(DFT) is the total energy of the crystal unit cell derived by solving the KohnSham equations within the chosen functional (in the present case B3LYP), and E(D) is the pure dispersive contribution deriving from the classical London formula3133 as redefined by Grimme to handle short-range interactions already accounted for by the E(DFT) term.12 Empirical parameters (i.e., van der Waals Rvdw radii and atomic C6 coefficients) were taken from Table 1 of ref 12 Rij,g is the interatomic distance between atoms i in the reference cell and j in the neighboring cells at distance |g|. A cutoff distance of 25.0 Å was used to truncate the summation over lattice vectors. The D correction is then parametric in internuclear distances, and is weighted against the DFT component of energy by a suitable cutoff function fdmp ðRij, g Þ ¼
1 1 þ expð dððRij, g =RvdW Þ 1ÞÞ
to avoid double counting of the energy contributions for shortrange interactions already accounted for by the DFT calculation. Here, Rvdw is the sum of the atomic van der Waals for the considered atomic pair ij. The s6 parameter is a global scaling factor whose value is functional dependent, as different functionals account for the short-range interactions in a different way. The original Grimme model has been modified here by rescaling van der Waals radii in the damping function (scaling factors: 1.3 for H and 1.05 for all other atoms), whereas the value of the global scaling factor s6 was set to 1. This variation of the original D term, named D*, was shown to be particularly effective for modeling molecular crystals. The proposed scaling factors were determined in previous work13 from a manual fitting procedure, 11180
dx.doi.org/10.1021/jp203132k |J. Phys. Chem. A 2011, 115, 11179–11186
The Journal of Physical Chemistry A searching for the best agreement between computed and experimental cohesive energies on a set of 14 molecular crystals. The B97-D method originally adopted by Grimme12 to include dispersion interactions in intermolecular complexes has also been used to have a further independent check of both B3LYP-D* and LMP2. It is based on the same functional form as the Becke’s GGA functional34 introduced in 1997 and supplemented by the London type formula described above; it has been shown to provide excellent results for a large number of intermolecular complexes characterized by H-bond and dispersive interactions similar to those found in the molecular crystals considered here. 2.2.2. The LMP2 Method As Implemented in CRYSCOR. The “local correlation” approach, first proposed by Peter Pulay in 198335 allows one to exploit the fundamentally local character of dynamic electron correlation, leading to quantum chemical methods that scale linearly with the size of the system. The CRYSCOR program implements this kind of approach for correlated calculations (MP2) in periodic systems employing the periodic HartreeFock (HF) solution in the basis of Gaussiantype orbitals (GTOs) provided by the CRYSTAL code.16 Wannier functions (WFs) are used to represent the occupied HF manifold: these are well-localized, symmetry-adapted, mutually orthogonal, translationally repeated functions. Conventionally, only correlations between the valence electrons are considered according to the commonly used “frozen core” approximation. The virtual manifold is described by the nonorthogonal set of projected atomic orbitals (PAOs). These orbitals are obtained by projecting individual AOs out from the occupied subspace, and they are as well localized as WFs. For a detailed description of the CRYSCOR periodic LMP2 implementation and theory and capabilities of CRYSCOR, we refer to previous work.10,11 Here we just recall its main features. The total LMP2 energy ELMP2 can be written as a sum of pair energies EijLMP2, each corresponding to two-electron excitations from the given ij-pair of WFs. Translational symmetry allows us to impose the first WF (i) to belong to the reference unit cell. For a given ij pair, the local-correlation ansatz restricts the virtual space, allowing excitations only to the PAOs in the pair-domain, i.e., those, spatially close to at least one of the WFs in the pair. This leads to the linear scaling of the LMP2 method with respect to cell size. In the present work, the local excitation domains have been restricted to the molecular units, to which the corresponding Wannier functions belong, a choice that has been shown to reduce BSSE effects.36 Approximate integral evaluation techniques, such as density fitting,37 make correlation calculations feasible for systems presently up to more than 100 atoms in the unit cell and relatively large basis sets. Translational and point symmetry are fully exploited by the code at all levels, which gives extra saving in both computer time and memory. The largest calculations in this work have been possible thanks to a newly developed parallel implementation38 of CRYSCOR. MP2 calculations are, in this framework, an essential reference in order to validate the reliability of the other methods. 2.3. PIXEL Calculations. In the PIXEL method,17,39 a molecule is described by an ab initio MP2/6-31G(d,p) molecular electron density, which is then repeated in space according to the desired aggregation geometry, either for molecular dimers, or for full crystal structures, represented by arrays of 100150 molecules obeying space-group symmetry. Electron density points from the MO calculation are contracted into charge pixels, or elementary cubes including n n n points of the original
ARTICLE
Scheme 1. PIXEL Definition of the Coulomb Interaction between Molecules A and Ba
a
Atomic nuclei as large blue dots, electron pixels as small red dots. See text for details.
density, with n = 4 or 3. A molecule is thus represented by its nuclear positions and typically some 10 000 to 20 000 charge density pixels. The charge density clouds are juxtaposed in a rigid fashion without modification due to molecular approach. PIXEL yields a nonempirical Coulombic energy integral plus separate semiempirical polarization, overlap repulsion, and dispersion terms. The Coulombic interaction energy is calculated by numerical integration over nucleinuclei, nucleipixel and pixel pixel terms. For a molecular pair AB of the crystal, the Coulomb interaction is (see also Scheme 1)
∑ ∑ qi qk =Rik þ ∑ ∑ Zjqi =Rij i, B k, A i, B j, A þ ∑ ∑ Zm qk =Rkm þ ∑ ∑ Zm Zj =Rjm m, B k, A m, B j, A
EEL, AB ¼
EEL ¼ 1=ð4πε°Þ
∑B ∑A EEL, AB
in which the nuclei of molecule A and B are labeled as Zm and Zj with the associated pixel charges qi and qk, respectively, the distances between these elements are labeled as R, and ε° is the vacuum dielectric constant (see Scheme 1). The sums to calculate EEL are computed by direct summation in PIXEL until a cutoff distance between the AB pair that ensures a convergence of the energy within fractions of kJ/mol. This approach runs into difficulties for highly polar compounds, when molecular dipoles align along a polar crystal direction. In such cases, the central cell of the lattice summation is not representative of the entire crystal and termination effects become non-negligible. This problem can be dealt with by EwaldBertaut reciprocal space methods40 or, much more conveniently, by the van EijckKroon dipole integration method,41 which gives identical results at a much lower computational and programming cost. In PIXEL, the molecular dipole is calculated from the contracted charge density map, and the total cell dipole moment p is evaluated by a vector sum of all molecular dipoles in the unit cell. The cell dipole energy correction to the direct EEL is then41 EEL ðtotalÞ ¼ EEL þ Eðcell dipoleÞ ¼ EEL 2πp2 =ð3V cell Þ in which Vcell is the cell volume. For crystals of noncharged molecules this correction seldom exceeds a few kilojoules per mole, but the correction becomes quite substantial for example in zwitterionic natural aminoacids, which are likely to pack in polar space groups. Polarization is evaluated as a sum of terms resulting from the electric field at each pixel in the linear polarization approach. Dispersion in estimated by a London-type formula involving distributed polarizabilities αi at each pixel, including a damping 11181
dx.doi.org/10.1021/jp203132k |J. Phys. Chem. A 2011, 115, 11179–11186
The Journal of Physical Chemistry A
ARTICLE
function; the total energy is a sum over pixelpixel terms: EDisp, AB ¼ ð 3=4Þ
EOS f ðRÞαi αj =½4πðε°Þ2 ðRij Þ6 ∑ ∑ i, A j, B
f ðRÞ ¼ exp½ ðD=Rij 1Þ2 ðfor Rij < DÞ D is a damping length parameter, set here to 3.0 Å or to 2.4 Å for nitro compunds. EOS is the “oscillator strength”, taken as the average local ionization potential. Repulsion is calculated as being proportional to the overlap integral between electron densities, corrected by a term in the difference of electronegativity. The scheme has now reached a decade-long maturity, and has been shown to be reliable in the description of dimerization energies and of lattice energies of organic compounds.19,39,4247
3. RESULTS 3.1. Comparison of Calculated Lattice Energies with Experimental Sublimation Enthalpies. We compare directly
experimental sublimation enthalpies, ΔH(s), to calculated lattice energies, Ulatt,; this refers to the process of transferring one mole of molecules in the solid state to one mole of rigid, noninteracting molecules in the gaseous state. This comparison neglects lattice-vibrational energies and differences in internal vibrational energy (the 2RT correction, sometimes applied, introduces perhaps more uncertainty than accuracy; see ref 48, p 210). Neglect of molecular relaxation on going from solid to vapor introduces an overestimation of ΔH(s); nearly all the compounds in our data set are, however, conformationally and tautomerically undeformable molecules. The temperature dependence of ΔH(s) is also neglected, as the experimental data are not always unequivocal, and the calculated numbers refer to an unspecified temperature; variations are of the order of 5 kJ mol1 every 100 K of temperature difference (see ref 48, p 191). Possible molecular association in the gaseous state is also neglected; this could be relevant in strongly H-bonding compounds. Large oscillations in experimental data from different sources are also frequent (see ref 48, Chapter 7). On the basis of experience gained while gathering the thermochemical data, we assume an expected variance of (10% over the experimental values. For these reasons, the discussion of discrepancies for single cases is hardly warranted, and the comparison offers mostly an overall guideline or a trend over a large database. 3.2. LMP2 versus Semiempirical Methods. Among our data set of 60 crystal structures, lattice energies were calculated by the full periodic LMP2 method for seven crystals (see Table 1), chosen to represent very weak dispersive interactions, polar forces, and strong hydrogen bonding, and compared with DFT-D and PIXEL results. DFT-D predicted lattice energies are not too sensitive to the basis set, while for LMP2 the choice is more critical;49,50 low-exponent polarization functions must be added to the 6-31G(d,p) basis set. Periodic-orbital LMP2 is still restricted by the computational cost associated with its present implementation and by the amount of man work associated with each calculation: the definition of local subspaces in local postHF methods adopted in CRYSCOR cannot be considered as a black-box procedure. Remarkably, such calculations are nevertheless feasible for crystals such as benzene, succinic anhydride, and nitroguanidine. DFT-D and periodic LMP2 are in good agreement, further validating the empirical D* formulation. Lattice energies from the semiempirical method generally exceed those by LMP2 by
Figure 1. Calculated lattice energies and experimental sublimation enthalpies for 60 organic crystals. The solid lines mark the (10% error range in ΔH(s). Single refcodes marked above the straight lines are for dichlorobenzene, trichlorobenzene, hexachlorobenzene, dichlorodibenzodioxin, and tetrachlorodibenzodioxin (FELSEU). Refcodes marked below the lines are acetic acid, formamide, oxalic and malonic acid. Pairs of points in a square frame (MLEICA, TNBENZ, ANTQUO, TNOXYL, PYMDAN, UREAXX, BENZDC, MELAMI, and CTMTNA; see text for the chemical names) are consistent indications by DFT-D and PIXEL.
510%; this observation could be a valuable suggestion for further fine-tuning the semiempirical part of DFT-D, within the known limits of the MP2 level of theory. Comparison with experiment is ridden with difficulties, as discussed above. In particular, the too large value of the lattice energy of benzene by LMP2 is expected as perturbative methods such as MP2 always overestimate the interaction energy of highly polarizable molecules. For instance, it is well-known from molecular studies51,52 that the MP2 method overestimates the interaction energy between dimers of benzene, while higher order methods such as CCSD(T) provide the correct answer. Estimates on periodic systems based on ab initio two-body interactions confirm this behavior also for the benzene crystal, with MP2 overestimating the cohesive energy by ∼10 kJ/mol with respect to CCSD(T) results, at basis set convergence.19,23,53 Our fully periodic LMP2 result is coherent with this picture. In the case of formic acid crystal, the LMP2 result is almost coincident with the experimental value and in better agreement than the other methods. The predicted value for formamide is quite good, and closer to experiment than the estimate obtained by Beran et al. using the two-body MP2 approach and comparable basis set.53 The ΔH(s) of urea and oxalic acid are peculiar cases as discussed below. PIXEL does slightly better against experiment, probably due to its larger number of parameters and to the fact that the sublimation enthalpies were among the training set in parameter optimization. For nitroguanidine, computed results are slightly overestimated with respect to experiment, but the agreement is within the expected variance of (10%. In order to provide an independent check of the LMP2 results, the B97-D method has been used to compute the lattice energies for the same subset of the seven molecular crystals subject to LMP2 calculations. Results of Table 1 show a general good agreement 11182
dx.doi.org/10.1021/jp203132k |J. Phys. Chem. A 2011, 115, 11179–11186
The Journal of Physical Chemistry A
ARTICLE
Table 2. Results of Calculations for the Crystals of Zwitterionic L-Alanine (LALNI24) and DL-Alanine (DLALNI04) PIXEL crystal L-alanine DL-alanine
a
Figure 2. PIXEL and D* dispersion energies for the structures of Figure 1. The least-squares straight line is y = 0.918x, R2 = 0.904. kJ mol1 units. The largest outliers are hexachlorobenzene, trichlorobenzene, nitroguanidine, and pyromellitic dianhydride.
between LMP2, B97-D and experiment giving reassuring evidence of the effectiveness of the LMP2 method. Mean deviation between B3LYP-D* and B97-D data is about 6 kJ/mol, with larger differences (up to 9 kJ/mol) for benzene, oxalic acid, and nitroguanidine crystals. 3.3. Lattice Energies: Comparison between PIXEL, DFT-D, and ΔH(s). The crucial issue for organic compounds is the treatment of dispersion energy,. The two semiempirical methods used here (PIXEL and DFT-D, the latter is the B3LYP-D* method) include a dispersion energy term of different analytical form, but essentially stemming from the same premises, i.e., the London approach. The comparison among DFT-D, PIXEL, and experiment was then extended by computing the lattice energies for crystals of 60 multifunctional organic materials. Under the proviso of the discussion in section 3.1, the agreement with experimental sublimation enthalpies is quite satisfactory for both kinds of calculation (Figure 1; data points for Figures 1 and 2 are available as Supporting Information, Table S1). PIXEL does slightly better than DFT-D in this respect, as already noted and for the reasons exposed in the preceding section. When points outside the brackets of the expected uncertainty in Figure 1 are considered, a few trends appear: (a) In a few cases (MLEICA: maleic anhydride; TNBENZ: trinitrobenzene; ANTQUO: anthraquinone; TNOXYL: trinitroxylene), the consistent response of DFT-D and PIXEL suggests that the experimental values may be in error by excess. (b) The DFT-D lattice energy of chloro-aromatic compounds is smaller than the observed ΔH(s), pointing to a possible correction in the dispersion contribution for chlorine; in fact, the PIXEL dispersion energies are consistently larger than the DFT-D dispersion energies for these crystals;
SG
ref
55
DFT-D
a
coul
P212121 277 269 267.3
cel dip 0.0
pol
disp
rep
tot
88.1 81.6 184.1 252.9
Pna21 274 272 109.4 161.4 61.5 83.7 185.6 230.4
B3LYP-D*/6-31G(d,p).
(c) The DFT-D lattice energy of carboxylic acids, amides, and urea is systematically larger, sometimes much larger, than the corresponding ΔH(s). The indication is consistent also in PIXEL for urea, 1,3-phthalic acid (BENZDC), and melamine. While a possible insufficient parametrization of DFT-D cannot be ruled out, more likely explanations invoke (i) some critical reconsideration of the experimental data: either an experimental value in error by defect, or some aggregation phenomenon in the gaseous state; (ii) lattice energy is computed with respect to the total energy of the molecule as cut out from the bulk, so that no relaxation effects are taken into account, which would influence more strongly the molecule than the crystal, thus reducing the computed lattice energy. Both explanations are actually very likely considering that these are all strongly hydrogen-bonding compounds. (d) The consistent indication of a lattice energy much larger than ΔH(s) for PYMDAN (pyromellitic dianhydride) and CTMTNA (trinitro-triazacyclohexane) both nonhydrogen bonding compounds, strongly suggests that the experimental value is wrong by defect. A glance at the literature for PYMDAN4 reveals a very old experimental determination, published in a non peer-reviewed journal. All considered, the comparison between calculation and experiment shows that there are indeed very few cases of large discrepancies between experiment and calculation, for which no convincing explanation can be found. The picture in Figure 1 shows that theoretical methods are able to reproduce thermochemical data, but also that theory can be sometimes more accurate than experiment and can give a considerable help in the critical evaluation of experimental data. We suggest that calculations by PIXEL or DFT-D can provide reliable estimates of sublimation enthalpies, possibly more accurate, and certainly much less expensive than experiment. 3.4. Dispersion Energies by DFT-D and by PIXEL. Figure 2 shows the excellent agreement between PIXEL and D* empirical dispersion energies, lending mutual support to both methods. The least-squares fitting indicates that D* dispersion energies are on the average about 8% smaller than PIXEL ones, with large deviations for the small D* dispersion contributions in chloro compounds, which is the cause of the underestimation of the lattice energy, as already noted. The agreement between two separate sets of dispersion energies and (with a few exceptions) of total lattice energies also implies that the sum of PIXEL Coulombic, polarization, and repulsion terms is consistent with total, nondispersive DFT energies. We take these results as a robust confirmation of the definition of dispersion energy, and as a confirmation that the dispersion energy term can be consistently and reliably evaluated by a London -type approach. 3.5. Lattice Sums in Polar Space Groups: Cell Dipole Energy Corrections. For the crystals of noncharged molecules 11183
dx.doi.org/10.1021/jp203132k |J. Phys. Chem. A 2011, 115, 11179–11186
The Journal of Physical Chemistry A
Figure 3. PIXEL, MP2/6-31G(d,p), DFT-D (B3LYP-D*), and uncorrected DFT-B3LYP dimer interaction energies in the crystal of nitroguanidine, O2NNdC(NH2)2. Data points and structure drawings are available as Supporting Information, Table S4.
ARTICLE
Figure 6. As in Figures 3 and 5, for molecular pairs extracted from the crystal structure of naphthoquinone. Data points in Table S6.
Figure 4. The E dimer (left) and the L dimer (right) in the nitroguanidine crystal (see Figure 3; H white, C gray, N blue, O red; drawing by QuteMol).61 Figure 7. The A, B and C molecular pairs from the crystal structure of naphtoquinone (see Figure 6, drawing by QuteMol61).
Figure 5. As in Figure 3, for molecular pairs extracted from the crystal structure of urea. See data points in Table S5.
listed in Figure 1, the cell dipole energy is negligible, the only exceptions being the small contributions for coumarin (4.2 kJ/ mol) and 1,3-dinitrobenzene (1.2 kJ/mol). On the other hand, a test calculation for zwitterionic alanine shows that the correction for the crystal structure in the polar space group (which unexpectedly is for the racemic crystal) is substantial. Indeed the DL-alanine crystal is a paradigm of polarity, and its different
terminations along the c-axis have been long exploited for testing the effects of tailor-made impurities on crystal growth.54 The agreement between PIXEL and DFT-D is very good for the nonpolar L-alanine crystal (Table 2), and the van Eijck Kroon41 cell dipole correction brings the PIXEL Coulombic energy of DL-alanine in line with that of the crystal of L-alanine, as it should be since the two crystals have very nearly the same lattice energy.55 The total energy of DL-alanine is still underestimated because, at the moment, PIXEL does not include an analogous correction to be applied for the polarization energy. 3.6. Molecular Pairs from Crystal Structures. The cohesive energies of molecular dimers rigidly extracted out of their crystal were also calculated by PIXEL, DFT-D, and the canonical counterpoise-corrected MP2/6-31G(d,p) method. The following examples were chosen: (i) nitroguanidine,56 with strongly stabilizing and moderately to-strongly destabilizing molecule molecule terms,45 (ii) urea,57 a benchmark of the hydrogen-bond description; and (iii) naphthoquinone,58 a non-hydrogen bonding system including a mix of dispersive and strongly polar interactions. Data points and structure drawings for the dimers in Figures 3, 5 and 6 are available as Supporting Information, Tables S2S4. Figure 3 shows the almost perfect agreement between PIXEL, DFT-D, and MP2 for dimers out of the nitroguanidine crystal. 11184
dx.doi.org/10.1021/jp203132k |J. Phys. Chem. A 2011, 115, 11179–11186
The Journal of Physical Chemistry A Uncorrected B3LYP energies are acceptable for dimers where Coulombic contributions dominate, but are in serious error for pair E (too little stabilization) and for pair L (too much destabilization) where parallel stacking (Figure 4) is dispersion-driven. The result conforms to the different physical nature of crystal forces and once again support the appropriateness of the treatment of dispersion energy. For urea, Figure 5 shows a remarkable coincidence between DFT-D and PIXEL. For short-distance dimers, B3LYP and MP2 results are less cohesive, while for other dimers there is in good agreement even with uncorrected DFT. For naphthoquinone, Figure 6 again shows an excellent agreement between DFT-D and PIXEL. Dimers A and C (Figure 7) are dispersion-dominated stacked pairs, and accordingly, uncorrected B3LYP energies are largely in error. In dimer B the molecules are coplanar with a substantial Coulombic term from approach of CdO and CH groups; accordingly, the B3LYP energy is acceptable. We underscore the valuable information on the performance of the methods, vis-a-vis the physical nature of the interaction, that can be obtained by an analysis of structure-energy relationships in carefully selected dimers.
4. CONCLUSION AND PERSPECTIVE One of the main goals of the present work is to present a description of the comparative feasibility of intermolecular energy calculations in large assemblies of organic compounds, by comparing results from increasingly accurate and costly methods, i.e., (i) PIXEL, which is based on physically meaningful empirical expressions for the interactions; (ii) DFT-D, density functional methods with a posteriori empirical correction for the dispersion interactions; and (iii) LMP2, a full ab initio quantum mechanical method which corrects the HartreeFock energy of the crystal for the electron correlation using the MøllerPlesset perturbation theory in its local formulation. As the LMP2 approach represents the state of the art in solid state post-HF methods, it is of cornerstone importance for a first-principle validation of empirical dispersion corrections. Indeed, comparisons among different methods are an essential aid in the selection of proper means for widespread application. The overall agreement among several computational schemes, as demonstrated here on a wide range of lattice energies, and the detailed agreement on particular aspects of the physics of the interaction, including the crucial dispersive effects, as demonstrated by energy-structure correlations in molecular dimers, add confidence and clarify the available range of tools and costs. PIXEL is the method of choice in the screening of a large number of molecular crystals, an essential step to predict crystal polymorphism or to study crystal growth processes. The present version does not allow for consistent optimization of the molecular structures, so that, the geometry of the unit to be replicated by the symmetry operators of a given space groups during crystal prediction should come from either experiment or quantum mechanical calculations and will be kept rigid. Furthermore, molecular crystals in which solid state reactions are taking place, as in the case of incipient proton transfer, are outside the PIXEL applicability. In all these cases, DFT-D, due to the modest computational resources, is the method of choice: it can adopt PIXEL crystal structures as initial guess for full geometry optimization (both unit cell size and internal coordinates) and can rather easily follow solid-state reactions with the ability to characterize transition state structures using well-established search algorithms. If a bright future for the
ARTICLE
characterization, prediction and eventual control of complex organic materials is in view, quantitative energetics are indispensable. Theoretical methods different from those adopted in the present work, to provide even more accurate lattice energies, can in principle be envisaged, but are not affordable for application in a mature and efficient way to a large and complex data set of molecular crystals like the one considered in this work; in the meantime, widespread use of DFT-D or PIXEL methods will foster significant advances in molecular simulation for crystal structure prediction, in molecular engineering, and in general in all molecular-level nanosciences.59 MP2 calculations are, in this framework, an essential reference in order to validate the reliability of the other methods.
’ ASSOCIATED CONTENT
bS
Supporting Information. Detailed lists of data points in Figures 13, 5, and 6 (Tables S1S4), and structural drawings for dimers extracted from crystal structures. This material is available free of charge via the Internet at http://pubs.acs.org.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT For encouragement and advice in this project and through the years we are greatly indebted to Carla Roetti and Cesare Pisani, both prematurely passed away in recent times. Carla played an essential role in the development and maintenance of the CRYSTAL code. Cesare gave an outstanding contribution to the periodic-LMP2 code CRYSCOR. Carla and Cesare were great scientists and beloved friends and will be sorely missed by all of us. ’ REFERENCES (1) Goerigk, L.; Grimme, S. J. Chem. Theor. Comput. 2010, 6, 107. (2) Hobza, P.; Sponer, J. Chem. Rev. 1999, 99, 3247. (3) Jurecka, P.; Sponer, J.; Cerny, J.; Hobza, P. Phys. Chem. Chem. Phys. 2006, 8, 1985. (4) Kannemann, F. O.; Becke, A. D. J. Chem. Theor. Comput. 2010, 6, 1081–1088. (5) Podeszwa, R.; Szalewicz, K. Chem. Phys. Lett. 2005, 412, 488–493. (6) Riley, K. E.; Pitonak, M.; Cerny, J.; Hobza, P. J. Chem. Theor. Comput. 2010, 6, 66–80. (7) Zhao, Y.; Truhlar, D. G. J. Chem. Theor. Comput. 2005, 1, 415–432. (8) Schulman, J. M.; Kaufman, D. N. J. Chem. Phys. 1970, 53, 477. (9) Binkley, J. S.; Pople, J. A. Int. J. Quantum Chem. 1975, 9, 229. (10) Pisani, C.; Busso, M.; Capecchi, G.; Casassa, S.; Dovesi, R.; Maschio, L.; Zicovich-Wilson, C. M.; Sch€utz, M. J. Chem. Phys. 2005, 122, 094113. (11) Pisani, C.; Maschio, L.; Casassa, S.; Halo, M.; Sch€utz, M.; Usvyat, D. J. Comput. Chem. 2008, 29, 2113. (12) Grimme, S. J. Comput. Chem. 2006, 27, 1787. (13) Civalleri, B.; Zicovich-Wilson, C. M.; Valenzano, L.; Ugliengo, P. CrystEngComm 2008, 10, 405. (14) Ugliengo, P.; Zicovich-Wilson, C. M.; Tosoni, S.; Civalleri, B. J. Mater. Chem. 2009, 19, 2564. (15) Dovesi, R.; Orlando, R.; Civalleri, B.; Roetti, C.; Saunders, V. R.; Zicovich-Wilson, C. M. Z. Kristallogr. 2005, 220, 571. (16) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, 11185
dx.doi.org/10.1021/jp203132k |J. Phys. Chem. A 2011, 115, 11179–11186
The Journal of Physical Chemistry A I. J.; D’Arco, P.; Llunell, M. CRYSTAL09; Universita di Torino: Torino, Italy, 2009. (17) Gavezzotti, A. Mol. Phys. 2008, 106, 1473. (18) Bucko, T.; Hafner, J.; Lebegue, S.; Angyan, J. G. J. Phys. Chem. A 2010, 114, 11814. (19) Dunitz, J. D.; Schweizer, W. B. Chem.—Eur. J. 2006, 12, 6804. (20) Dyakonenko, V. V.; Maleev, A. V.; Zbruyev, A. I.; Chebanov, V. A.; Desenko, S. M.; Shishkin, O. V. CrystEngComm 2010, 12, 1816. (21) Liu, Y.; Goddard, W. J. Phys. Chem. Lett. 2010, 1, 12550. (22) Moellmann, J.; Grimme, S. Phys. Chem. Chem. Phys. 2010, 12, 8500. (23) Ringer, A. L.; Sherrill, C. D. Chem.—Eur. J. 2008, 14, 2542. (24) Tsuzuki, S.; Orita, H.; Honda, K.; Mikami, M. J. Phys. Chem. B 2010, 114, 6799. (25) van der Streek, J.; Neumann, M. A. Acta Crystallogr. 2010, B66, 544. (26) Allen, F. H. Acta Crystallogr. 2002, B58, 380. (27) Gavezzotti, A.; Filippini, G. J. Phys. Chem. 1994, 98, 4831–4837. (28) Zicovich-Wilson, C. M.; Kirtman, B.; Civalleri, B.; RamirezSolis, A. Phys. Chem. Chem. Phys. 2010, 12, 1985. (29) Boys, S. F.; Bernardi, M. Mol. Phys. 1970, 19, 553. (30) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheesman, J. R.; Montgomery, J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmanu, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. GAUSSIAN03; Gaussian Inc.: Wallingford, CT, 2004. (31) Eisenschitz, R.; London, F. Z. Phys. 1930, 60, 491. (32) London, F. Z. Phys. 1930, 63, 245. (33) London, F. Z. Phys. Chem. 1930, B11, 222. (34) Becke, A. D. J. Chem. Phys. 1997, 107, 8554. (35) Pulay, P. Chem. Phys. Lett. 1983, 100, 151–154. (36) Sch€utz, M.; Rauhut, G.; Werner, H.-J. J. Phys. Chem. A 1998, 102, 5997. (37) Maschio, L.; Usvyat, D.; Manby, F.; Casassa, S.; Pisani, C.; Sch€utz, M. Phys. Rev. B 2007, 76, 75101. (38) Maschio, L. J. Chem. Theor. Comput. 2011, DOI: dx.doi.org/ 10.1021/ct200352g. (39) Gavezzotti, A. J. Phys. Chem. B 2002, 106, 4145–4154. (40) Williams, D. E. Acta Crystallogr. 1971, A27, 452–455. (41) van Eijck, B. P.; Kroon, J. J. Phys. Chem. 1997, 101, 1096–1100. (42) Aldridge, S.; Downs, A. J.; Tang, C. Y.; Parsons, S.; Clarke, M. C.; Johnstone, R. D. L.; Robertson, H. E.; Rankin, D. W. H.; Wann, D. A. J. Am. Chem. Soc. 2009, 131, 2231–2243. (43) Forgan, R. S.; Wood, P. A.; Campbell, J.; Henderson, D. K.; McAllister, F. E.; Parsons, S.; Pidcock, E.; Swart, R. M.; Tasker, P. A. ChemComm 2007, 4940–4942. (44) Gavezzotti, A. J. Chem. Theor. Comput. 2005, 1, 834–840. (45) Gavezzotti, A. Acta Crystallogr. 2010, B66, 396. (46) Muniappan, S.; Lipstman, S.; Goldberg, I. Chem. Commun. 2008, 1777–1779. (47) Wood, P. A.; Haynes, D. A.; Alistair, R.; Motherwell, W. D. S.; Parsons, S.; Pidcock, E.; Warren, J. E. Cryst. Growth Des. 2008, 8, 549– 558. (48) Gavezzotti, A. Molecular Aggregation; Oxford University Press: Oxford, 2007. (49) Slough, W.; Perger, W. F. Chem. Phys. Lett. 2010, 498, 97.
ARTICLE
(50) Maschio, L.; Usvyat, D.; Civalleri, B. CrystEngComm 2010, 12, 2429. (51) Sinnokrot, M.; Sherrill, C. D. J. Phys. Chem. A 2006, 110, 10656–10668. (52) Park, Y. C.; Lee, J. S. J. Phys. Chem. A 2006, 110, 5091–5095. (53) Beran, G. J. O.; Nanda, K. J. Phys. Chem. Lett. 2010, 1, 3480–3487. (54) Weissbuch, I.; Lahav, M.; Leiserowitz, L. Cryst. Growth Des. 2003, 3, 125–150. (55) Destro, R.; Soave, R.; Barzaghi, M. J. Phys. Chem. B 2008, 112, 5163–5174. (56) Murmann, R. K.; Glaser, R.; Barnes, C. L. J. Chem. Cryst. 2005, 35, 317. (57) Swaminathan, S.; Craven, B. M.; McMullan, R. K. Acta Crystallogr. 1984, B40, 300. (58) Gaultier, J.; Hauw, C. Acta Crystallogr. 1965, 18, 179. (59) Gubbins, K. E.; Moore, J. D. Ind. Eng. Chem. Res. 2010, 49, 3026. (60) Raabe, G. Z. Naturforsch. 2002, 57a, 961. (61) Tarini, M.; Cignoni, P.; Montani, C. IEEE Trans. Visualization Comput. Graphics 2006, 12, 1237.
11186
dx.doi.org/10.1021/jp203132k |J. Phys. Chem. A 2011, 115, 11179–11186
