Is the Induction Energy Important for Modeling Organic Crystals

522

J. Chem. Theory Comput. 2008, 4, 522-532

Is the Induction Energy Important for Modeling Organic Crystals? Gareth W. A. Welch,† Panagiotis G. Karamertzanis,† Alston J. Misquitta,†,‡ Anthony J. Stone,‡ and Sarah L. Price*,† Christopher Ingold Laboratory, Department of Chemistry, UniVersity College London, 20 Gordon Street, London WC1H 0AJ, U.K., and UniVersity Chemical Laboratory, Lensfield Road, Cambridge CB2 1EW, U.K. Received October 16, 2007

Abstract: We compare two methods for estimating the induction energy in organic molecular crystals by approximating the charge density polarization in the crystalline state. The first is a distributed atomic polarizability model combined with distributed multipole moments, derived from ab initio monomer properties. The second uses an ab initio calculation of the molecular charge density in a point-charge field. Various parameters of the models, such as the rank of polarizability model, effect of self-consistent iterations, and damping, are investigated. The methods are applied to a range of observed and predicted crystal structures of three particularly challenging molecules, namely oxalyl dihydrazide, 3-azabicyclo[3,3,1]nonane-2,4-dione, and carbamazepine, as well as demonstrating the importance of induction in the naphthalene crystal. The two models agree well considering the different approximations made, and it is shown that the induction energy can be an important discriminator in the relative lattice energies of structures with substantially different hydrogen-bonding motifs.

1. Introduction The importance of the induction energy for modeling crystal structures has been the subject of much debate.1-5 In ionic crystals, the empirical shell model has long been used to model the induced dipole from the strong electrostatic fields.6,7 The enhanced dipole moment of water in the liquid state has led to a plethora of model intermolecular potentials that include a simple polarizability model.1 Indeed, the development of the tinker force field for biological modeling now includes a dipole polarizability term,8,9 with the atomic polarizabilities derived by Thole.10 Reliable methods of estimating the polarizability models are demanding for two reasons: first they require a large basis set and high quality wavefunctions to obtain converged polarizability tensors, and second, for all but the smallest of molecules, an accurate description of the molecular polarizability can be obtained only with a distributed polarizability * Corresponding author e-mail: [email protected]. † University College London. ‡ University Chemical Laboratory.

10.1021/ct700270d CCC: $40.75

model. Until recently, these requirements have posed an almost insurmountable problem for modeling the induced moments in organic molecules. Furthermore, there has been uncertainty as to whether an approximate polarizability model would lead to greater errors than complete neglect, and almost all organic crystal structure modeling has been performed with model intermolecular potentials that do not include an explicit polarizability term.11-14 Instead, polarizability effects are to some degree approximately absorbed in the empirically fitted repulsion-dispersion parameters. Notable exceptions are the ab initio potentials developed for crystal structure prediction studies of simple alcohols and alkanes,5,15 glycol and glycerine,16 and some self-consistent molecular mechanics work on peptides.8,17 There is mounting evidence that induction effects are important within crystal structures of even nonpolar molecules. The PIXEL method, evaluating lattice energies by integrating over semiempirical functions of the in vacuo electron densities placed in the crystal lattice, shows that the induction energies are significant.4,18-20 Furthermore, experimental analysis of the naphthalene crystal shows © 2008 American Chemical Society

Published on Web 02/19/2008

Induction Energy for Modeling Organic Crystals

evidence of significant polarization on the charge density of the naphthalene molecule.21 However, the quantification of the induced moments from the difference between the molecular charge density in the crystal and in isolation is hampered by the uncertainty in how crystal charge density should be partitioned between the constituent molecules. This highlights a major difficulty with the evaluation of induction energies in molecular crystals: the charge distributions of molecules in van der Waals contact overlap so much that the long-range models of polarizability may not be valid. The induction energy is best understood from the theory of intermolecular forces. In particular, the recently developed symmetry-adapted perturbation theory using density function theory (SAPT(DFT))22-24 provides us with a computationally efficient and accurate method for calculating the induction energy of dimers of small polyatomics. The induction energies from SAPT(DFT) include penetration and chargetransfer effects and therefore provide us with an important benchmark against which approximations can be tested.25,26 However, being a two-body theory, SAPT(DFT) does not allow us to estimate the induction energy of the condensed phase. One way of doing this is to use polarizability models. Recently, some of us have developed a method for obtaining distributed polarizability models that is well suited for small polyatomic molecules of around 30 atoms. The WilliamsStone-Misquitta (WSM) method25,26 allows us to obtain distributed polarizabilities from the ab initio properties of isolated molecules that are optimal at a given rank. From comparisons with SAPT(DFT) induction energies of a variety of dimers, ranging from HF to benzene,26 we know that the damped WSM models are able to describe not only the longrange induction energy but also an induction energy at shortrange, even in the most testing area of hydrogen-bonding contacts. These models result in errors of 2-7% of the dimer interaction energy at typical contact distances. The error would be larger if we included hyperpolarizability effects that are not included in the WSM models. However, for condensed phases the errors are smaller than for van der Waals dimers because of the large number of longer range interactions, for which the WSM models are extremely accurate. Therefore, these polarizability models give us a very powerful tool for computing the induction energy of an organic crystal. Yet another way of approximating the induction energy of the crystalline phase relies on the ab initio evaluation of the molecular charge density, with the field of the surrounding molecules represented by point charges. When done selfconsistently, we obtain an electronic response to point charge field model (SCERP) which does include some of the effects of electron penetration, because the point charges are fitted to the electrostatic potential close to the van der Waals27 surface. But this model is limited by the accuracy that can be attained by the point charge model and the lack of chargetransfer effects. Yet, once again, these are short-range effects and, for reasons explained above, are not expected to make a significant contribution to the interaction energy of the crystal. The WSM polarizability model has been validated for dimers against SAPT(DFT) energies.26 SCERP provides an

J. Chem. Theory Comput., Vol. 4, No. 3, 2008 523

Figure 1. Molecules used in this investigation: (a) naphthalene, (b) oxalyl dihydrazide, (c) 3-azabicyclo[3,3,1]nonane2,4-dione, and (d) carbamazepine. Double arrows indicate that two atoms independently have a torsion angle defined along the same bond.

independent test of the polarizability models in the condensed phase. This paper uses these models to estimate the induction effects in a range of molecular crystals. Our purpose is to establish whether the contributions are sufficiently important that we should implement polarizability models in the crystal structure modeling program DMAREL.28,29 The overall aim of the paper is to determine the importance of the induction energy in organic crystal structures, particularly its relevance to the field of organic crystal structure and polymorph prediction,30 whose promise for aiding the design of new materials and the selection of solid form for pharmaceutical development31 is severely compromised by uncertainties in the estimation of relative lattice energies. Four contrasting examples are considered for which the molecular structures are shown in Figure 1. The naphthalene (1a) crystal is investigated as a nonpolar system. Charge density studies21 have shown a change in the electron distribution in the region of the C-H bond involved in a C-H‚‚‚π interaction in the crystal structure. The other examples are all tests of the differences in induction energy corresponding to different types of hydrogen bonding, as the electrostatic fields involved in hydrogen bonding are among the strongest in crystal structures of neutral organic molecules. The relative induction energies of the 5 different polymorphs of oxalyl dihydrazide32 (1b) are examined because of the plurality of hydrogen-bonding geometries sampled, including one with a significant intramolecular component. The relative induction energies for sets of experimentally observed and hypothetical crystal structures of 3-azabicyclo[3,3,1]nonane-2,4-dione (1c) and carbamazepine (1d) are computed, to investigate whether modeling induction could improve the prediction of relative lattice energies of crystal structures based on doubly hydrogenbonded dimers or chain motifs. In both cases, the predictions that the two types of crystal structure were energetically competitive inspired extensive polymorph screening studies to search for the alternative motif.33,34 For carbamazepine, all known polymorphs are based on a doubly hydrogenbonded amide dimer (although it does adopt a catemer in a solid solution35), whereas the catemer is marginally more

524 J. Chem. Theory Comput., Vol. 4, No. 3, 2008

Welch et al.

stable according to current modeling.36-38 On the other hand, 3-azabicyclo[3,3,1]nonane-2,4-dione adopts an imide catemer in all its solid forms,34 although many of the participants in the 2001 international blind test of crystal structure prediction39 predicted a dimer structure as more stable. We first define the new computational methods for estimating the induced distributed moments and induction energy contribution to the lattice energy, before applying the two models to this range of issues in organic solid-state chemistry.

2. Method Using a finite cluster of molecules, sufficiently large to obtain converged electrostatic energies, we estimate the effect of the crystalline environment on the molecular multipole moments calculated using the distributed polarizability and the SCERP models. As will be described below, the modified molecular multipole moments are then used to estimate the induction contribution to the crystal lattice energy. We first describe the cluster model, then the two methods of estimating the induced moments, and finally the method of evaluating the induction energy in the lattice. 2.1. Choice of Crystal Structures and Cluster and Molecular Models. For our calculations we use centrosymmetric crystal structures, from which the clusters are built. Numerical experimentation has shown that a cluster in which a central molecule is surrounded to at least 15 Å in all directions is large enough to converge the electrostatic energy of a molecule in the center, to that of an infinite lattice calculation using DMAREL. This typically means using a cluster of 5 × 5 × 5 unit cells. The crystal structure used for naphthalene was the 100 K X-ray structure.21 The molecular structure was optimized in vacuo at the MP2 6-31G** level and then pasted into the experimental structure by minimizing the rms overlap of the carbon atoms. Finally the crystal structure was relaxed using DMAREL with distributed multipoles derived using the same charge density as for the distributed polarizability model. For oxalyl dihydrazide, the five experimental crystal structures32 were refined to account for the X-ray determination of the proton positions that are important in the plurality of the hydrogen-bonding in these crystals.40 This DMAFLEX refinement36 optimized the lattice energy, including the MP2 6-31G** intramolecular conformational penalty, with respect to the nine torsions shown in Figure 1, and the crystallographic cell parameters and molecular positions. All covalent bond lengths and angles apart from the explicit torsion angles were reoptimized in the ab initio intramolecular calculation at each step. The key difference between the inter- and intramolecular bonding in the polymorphs (Figure 2) has been preserved, though the model for the polymorph is more dense that the experimental structure, resulting in one short N‚‚‚N distance of 2.73 Å. The rms difference41 between these refined structures and the experimental crystal structures was about 0.2 Å for the R, γ, and polymorphs and less than 0.6 Å for β and δ, for all nonhydrogen atoms in a 15-molecule cluster (see Table S1). The cluster sizes were varied (9 × 7 × 5 for R and , 7 × 5 × 7 for γ and δ, and 9 × 5 × 9 for β polymorphs) to give suitable supercells containing between 490 and 980 mol-

Figure 2. The two major intramolecular conformations of oxalyl dihydrazide. The β, γ, δ, and polymorphs contain stretched intramolecular hydrogen bonds, indicated by a dashed line. The torsion angles for all five polymorphs are given in Table S1, Supporting Information.

ecules that conformed to our requirement of 15 Å of material surrounding the polarizable molecule. The intermolecular electrostatic energy for all of these clusters is within 0.5 kJ mol-1 of the infinite lattice value (Supporting Information Table S2). The bicyclic structure of 3-azabicyclo[3,3,1]nonane-2,4dione42 makes it essentially rigid; therefore, we used the in vacuo MP2 6-31G** optimized molecular structure. A set of 8 low-energy crystal structures34 generated using this molecular conformation was considered to represent a range of packing arrangements within 3 kJ mol-1 of the global minimum lattice energy. We also examined the minimum obtained with the same computational model, starting from the 297 K experimental crystal structure.42 The set of structures include both the observed catemer and doubly hydrogenbonded dimer motifs in a range of space groups. The 5 × 5 × 5 unit cell clusters contained 250, 500, or 1000 molecules. For carbamazepine,33 we used DMAFLEX to relax the positions of the amide protons and torsion angles identified in Figure 1, for fifteen low-energy crystal structures36 obtained from a previous search.33 These structures covered a wide range of packings including those corresponding to known forms II, III, and IV. The carbamazepine clusters used in the polarizability calculations consisted of 5 × 5 × 5 unit cells and contained 250 to 1000 molecules. 2.2. Distributed Polarizability Model for Induced Moments. For large molecules, the single-center multipole expansion may not converge, and it is necessary to distribute the polarizabilities43-45 in an analogous way to the multipole moments. In general, the distributed polarizability R is dependent on the response of a moment at site a to a field at another site a′ in the same molecule. These nonlocal distributed polarizabilities can be obtained very accurately, and in a computationally efficient way, using the methods developed by Misquitta and Stone.46 However, the summation over two-site terms makes polarizability calculations for large molecules and clusters expensive computationally. The nonlocal polarizability description may be simplified using a mathematical transformation to remove the nonlocal terms. Such a transformation, put forward by Le Sueur and



Stone,47 effectively removes any explicit intramolecular polarization terms, resulting in a local distributed polarizability description that depends only on the electrostatic field exerted directly by the surrounding molecules. This is done at the expense of a loss in accuracy46,47 that increases with molecular size and rank of the polarizabilities. However, the localized polarizability description is still reasonably accurate, and the terms themselves are in good correspondence with what we might expect from physical arguments. Using these polarizabilities as anchors, the method described by Williams and Stone48 can now be used to refine the model. The refinement is done by tuning the distributed polarizabilities to reproduce the responses to point charges placed on a grid around the molecule. This final step in the WSM method25,26 for obtaining a local, distributed polarizability model partially absorbs the effects of higher rank polarizabilities in the lower rank terms. For this work, the distributed polarizabilities are localized to either rank 1 or rank 2 and refined by fitting to responses computed using linear-response DFT at randomly generated points26 between surfaces at 2 and 4 times van der Waals Bondi radii.49 We denote the WSM models as being L1 (induced-dipole), L2 (L1 plus induced-quadrupole), and L2/L1 (L1 on hydrogen sites, L2 for all other sites). The fitting procedure includes a penalty function to discourage drifting from the localized values,26 which has been found to result in accurate models and reduce the occurrence of unphysical values. All distributed properties, both polarizabilities and multipoles, have been derived using CamCASP,50 from charge densities calculated with DALTON51 using the asymptotically corrected52 PBE053 functional in combination with the Sadlej basis set.54 Having obtained the distributed polarizabilities, the change to the in vacuo multipole moments by polarization effects due to the crystal environment can now be estimated. In our work we obtain distributed multipole moments up to rank 4 using the revised version of Distributed Multipole Analysis.55 In the condensed phase the charge density distorts due to polarization effects, requiring an alteration in the multipole moments. The change in the multipole moments at a polarizable site a in molecule A, due to the static field of all other sites b in molecules B in the surrounding environment, is given by56 ∆Qat ) -

a ab RtV fuV(βRab)TVu (Qbu + ∆Qbu) ∑ ∑ ∑ B*A b∈B uV

(1)

a where RtV is the polarizability tensor for site a, quantifying the susceptibility of the multipole moment t on site a to be induced by the field arising from the static multipole moments Qbu and induced multipole moments ∆Qbu on all sites of other molecules. The subscripts t, u, and V, refer to the component of the multipole moments and run as ab 00,10,11c,11s... TVu is the interaction tensor containing the distance and orientational relation between sites a and b and their multipole components V and u, and fuV(βRab) is a TangToennies damping function that is assumed to depend only on distance, a damping parameter β, and the rank of multipoles represented by V and u. A damping function is used in an attempt to compensate for the divergence of the multipole expansion at small

intersite distances. Little is known about damping functions for induction;57 see ref 25 for a recent discussion. We use Tang-Toennies damping, which has been used to damp multipolar expansions of the dispersion energy. Examples show that it does not correct fully for the limitations of the multipolar model,26 but no better form has been proposed. The Tang-Toennies damping function has the form58 fuV(βRab) ) 1 -

(∑ ) n

(βRab)k

k)0

k!

exp(-βRab)

(2)

where n is the sum of the ranks of multipoles u and V, and has been effective in reducing the singular behavior of the induction energy when intersite distances are particularly short.26 The damping expression is used in the calculation of the induced moments, using β ) 2 x2IX, where IX is the first vertical ionization potential in atomic units.25 The values of β are as follows: oxalyl dihydrazide R 1.625; β 1.667; δ 1.650; 1.649; γ 1.657; naphthalene 1.547; carbamazepine 1.510; and 3-azabicyclo[3,3,1]nonane-2,4-dione 1.674. Thus, the distributed polarizability model estimates the induced moments in the crystal using eq 1, implemented in ORIENT,59 for a molecule, A, at the center of the cluster. At the first iteration with zero induced moments, the induction (2) and corresponds, within the approxienergy is Eind,d-class mations implicit in the truncation of the multipole and polarization expansions, to the second-order energy in the Rayleigh-Schro¨dinger theory60 for a 2-body system at large intermolecular separations (2) (X) Eind,pol

)

∑ r*0

|〈ΦX0 |Vˆ |ΦXr 〉|2 EX0 - EXr

(3)

where ΦXr and EXr are the eigenstates and eigenvalues of the monomer Hamiltonian of molecule X, and Vˆ is the intermolecular electrostatic potential operator arising from the rest of the system. The suffixes pol and d-class indicate, respectively, the root of the induction energy term in perturbation theory and a damped classical polarizability model. The coupled eqs 1 for the ∆Q are usually solved by (2-3) , iteration. After one iteration the energy becomes Eind,d-class where the change in energy corresponds to a third-order term in perturbation theory. After iteration to self-consistency, the induced multipoles correspond to all orders of the induction energy in the linear-response approximation. The total damped classical induction energy ∞

(2-∞) ) Eind,d-class

(k) Eind,d-class ∑ k)2

(4)

was used to examine the effect of polarization on the relative stability of the crystal structures. Figure 3 demonstrates that this iteration procedure makes a significant difference to the calculated induction energy, stabilizing the crystal. The rank 1 model converges rapidly, but higher-ranking polarizabilities do require damping. In practice, Figure 3 shows that the infinite summation in (3) can be truncated, depending on the model, to 5-8 iterations, which is sufficient to achieve convergence of 0.5 kJ mol-1.


Welch et al.

Figure 3. Convergence of Eind,d-class for polymorphs of oxalyl dihydrazide for several induction energy models. The plot shows (2-n+1) (2-∞) the error, Eind ,d-class - Eind,d-class, in the induction for different truncations of the infinite sum.

Figure 4. Fragments of carbamazepine used to calculate its atomic polarizabilities. The transferability of the polarizabilities calculated for these molecules to carbamazepine is described in the Supporting Information. The circled atoms are added in place of the 6-membered ring.

The carbamazepine molecule was too large for the WSM polarizability analysis due to computational limitations, and we adopted a different scheme for this molecule. The distributed multipoles were calculated in Gaussian03 using the nonasymptotically corrected PBE0 functional with the Sadlej basis set. The difference between the corrected and uncorrected functionals is insignificant for the calculation of electrostatic energies using distributed multipole moments, but the correction is essential for accurate polarizabilities. The polarizabilities were constructed from two molecular fragments, indicated in Figure 4. The structure of the fragments were held rigid at the MP2 6-31** in vacuo optimized geometry of carbamazepine, except the positions of the hydrogen atoms added in place of the 6-membered ring, which were optimized at the same level. Although polarizability is a molecular property, influenced by all sites,

it has been necessary to make the approximation of transferability (see the Supporting Information) for the polarizabilities calculated for these smaller molecules to the larger DMAFLEX minimized structures. 2.3. Self-Consistent Electronic Response to Point Charge Field Model (SCERP). We also present an alternative method of evaluating the effect of induction on the charge distribution directly using the Gaussian03 ab initio package.61 The CHELPG potential derived charges,27 which are fitted to a grid of points between the van der Waals atomic radii and 2.8 Å from the nuclei, were obtained for the isolated molecule from an aug-cc-pVTZ charge density with the PBE0 functional. These charges were placed on all the atomic sites of the same clusters as described in paragraph 2.1, except the central molecule which is described using augcc-pVTZ basis functions. A DFT calculation using the PBE0 functional is conducted for this molecule within the cluster of charges. The polarized charge density was analyzed by GDMA2.255 to obtain Qbu + ∆Qbu, and hence the induced multipole moments (up to hexadecapole) were obtained by subtraction of the multipoles obtained from the in vacuo calculation. The potential derived charges of the polarized charge distribution were then used in a further cluster calculation, and the process was repeated until the calculated induction energy had converged as for the distributed polarizability model. Figure 3 shows that, as with the WSM model, iteration is required to capture a significant part of the energy and that around half a dozen iterations are sufficient for convergence within 0.5 kJ mol-1. This method is more computationally expensive than using the multipole expansion and cannot be used for lattice energy minimization but can be used for testing aspects of the polarizability model. The resources required are almost independent of the number of charges used, and so very large



Eind,d-class(A) )

1

b ∆Qat f(tu)(βRab)Tab ∑ ∑ ∑∑ tu Qu tu

2 a∈A B*A b∈B

(5)

where the omission of the superscript implies that ∆Qat are converged induced moments. If the damping function is set to unity, then this equation is almost identical to the expression for the electrostatic energy Eelectrostatic(A) )

Figure 5. Electrostatic field difference around R oxalyl dihydrazide. The plot shows the norm of the difference between the electrostatic field vectors calculated from distributed multipole moments and from point charges. The surface is the van der Waals surface scaled by 1.8, which is accessible by the hydrogen-bonding protons. The maximum field difference displayed is 0.226 V/Å.

clusters could be used to check convergence with cluster size. The electron density is simply calculated in the fields of the surrounding background charges. Although the use of point charges to model electrostatic field is relatively crude, they are used here to induce a second-order response of the molecular charge density and within the self-consistent nature of the process. Penetration effects from the overlap of the charge distributions in the cluster are absent, except insofar as they are included in the fitting of the potential derived charges to points so close to the molecule. If the polar hydrogen sites are considered to have a van der Waals radius of zero,62 the region of interaction with surrounding nuclei in hydrogen bonding arrangement may be approximated by the van der Waals surface scaled by 1.8. In Figure 5 we present a comparison of the electrostatic field at this surface when calculated using multipole moments or point charges, in terms of the norm of the difference vectors at 19 814 points on the surface for R oxalyl dihydrazide. The mean difference is 0.08 V/Å (standard deviation 0.04 V/Å) which is less than 9% of the largest field, 0.92 V/Å, with the multipole moments. The highly localized nature of the error in the electrostatic field can be plainly seen in Figure 5, as dark blemishes around the hydrogen sites. For both the R and the polymorphs these regions coincide with the shortest hydrogen bonds seen in any of the crystal structures in this work, hence we anticipate the largest errors in our calculations to be for these crystals. 2.4. Calculation of the Induction Contribution to the Lattice Energy. We can evaluate the induction energy for a given crystal structure by the following method using the induced multipole moments that we have calculated by the methods described above. The following easily implemented method is not suitable for optimizing a crystal structure but allows a quick assessment of the importance of induction energy in organic molecular crystals. The classical polarization model for the induction energy is25,56

1

b Qat Tab ∑ ∑ ∑ ∑ tu Qu 2 a∈A B*A b∈B tu

(6)

and this can be exploited to estimate the induction energy of the crystal using the routines already implemented in DMAREL29 that evaluate this function and perform the lattice summations. Equation 5 has only one molecule bearing just the induced moments interacting with the electrostatic field of the rest of the crystal and hence cannot be directly calculated by DMAREL, assumes that all symmetry related sites bear equal (or inverted) multipole moments. However, substituting (Qt + ∆Qt/2) into eq 6 gives 1

∑ ∑ ∑∑ tu

2 a∈A B*A b∈B )

1

∑ ∑ ∑∑ tu

2 a∈A B*A b∈B

(

Qat

(

Qat

Tab tu

+

Qbu

)(

∆Qat 2

+

Tab tu

2

b ∆Qat Tab tu Qu

)

∆Qbu

+

Qbu

+

2 a ∆Qbu Tab tu Qt

2

+

)

b ∆Qat Tab tu ∆Qu

4

) Eelectrostatic(A) + Eind,d-class(A) + ∆Eerror(A)

(7)

Thus the induction energy can be calculated from three evaluations of the “electrostatic” contribution to the lattice energy, one where all molecules have the distributed multipole moments (Qt + ∆Qt/2) to get Eelectrostatic + Eind,d-class + ∆Eerror, a second with distributed multipole moments ∆Qt/2 to give ∆Eerror, and a third using only Qt to give Eelectrostatic. All three evaluations use Ewald summation for the charge-charge, charge-dipole, and dipole-dipole terms and sum all the other contributions in direct space for all molecules whose center of mass is within 15 Å. Since there is no facility to include damping of the electrostatic interactions in DMAREL, the necessary damping of the induction energy (5) is included for each iteration of the interaction of induced and static multipole moments in the cluster but is not applied in the final lattice energy calculation.

3. Results 3.1. Oxalyl Dihydrazide: The Effects of Rank, Refinement, and Damping. First we compare the energies calculated using the self-consistent electronic response to potential derived charges (SCERP) with WSM models (Figure 6). We find using the SCERP model for Eind,d-class that the R structure is stabilized the most, followed by . The β, γ, and δ structures are stabilized less than the form but by similar magnitudes to one another. Each of the WSM models follow the SCERP results except the L2 models,


Welch et al.

Figure 6. The induction energy of oxalyl dihydrazide for various polarizability models.

which find the relative polarization of the structures to be >R.β, γ, δ. This deviation can be explained by examining the crystal structures. In these five polymorphs there are very short intermolecular contacts. The shortest is an N‚‚‚H-N contact in which is shorter than the experimental value and indeed unphysically short. At such short distances the induction energy is very sensitive to distance, and in this case is almost certainly too large. The sensitivity can be reduced by using only rank 1 polarizabilities on hydrogen atoms, which does not lead to significant loss of accuracy overall. Experimentally, it has been difficult to fully characterize the relative stability of these polymorphs of oxalyl dihydrazide, due to a self-reaction that takes place prior to melting.32 However, lattice-energy methods that only model the intermolecular repulsion, dispersion, and electrostatic forces, including the conformational energy differences from ab initio gas-phase calculations, predict that the lattice energy of the R form is approximately -110 kJ mol-1, whereas the other four forms range from -130 to -138 kJ mol-1 (Supporting Information Table S3). Such a large energy difference is considered to be outside the range of possible polymorphism.63 By including a correction for the induction energy of the lattice, the predicted lattice energy of the R form becomes comparable with that of the β, γ, and δ polymorphs. It seems apparent that it is important to model charge density polarization for polymorphs that exhibit different intra- and intermolecular hydrogen bonding. This issue is being explored further using electronic structure calculations on oxalyl dihydrazide and other polymorphic systems.40 Our results relating to oxalyl dihydrazide strongly suggest that an iterated, damped polarizability model, based on the L1 or mixed L2/L1 models, agrees reasonably well with the self-consistent electronic response to point charges method. 3.2. Naphthalene. Induction is important not only for hydrogen-bonded systems. The crystal structure of naphthalene has been previously analyzed for experimental evidence of induced changes in the charge density.21 Our SCERP point charge model predicts an induction energy of -1.9 kJ mol-1 for the 100 K experimental crystal structure, using the molecular geometry optimized in vacuo. Although small in absolute terms, this is 31% of the electrostatic energy. A

Figure 7. Induced electrostatic energy surface for naphthalene. The energy is calculated from the SCERP model, for the van der Waals + 1.1 Å surface that is accessible by shortcontact nuclei. The atom numbering system reflects the symmetry of contacts within the crystal structure, not of the isolated molecule. The energy is calculated using a unit charge probe and ranges from -5.23 kJ mol-1 to +6.82 kJ mol-1.

damped WSM2/1 polarizability model estimates the induction energy to be 25% of the electrostatic energy. In comparison, the SCERP induction energy for oxalyl dihydrazide polymorphs is 18-38% of the electrostatic energy. Thus, in relative terms, even the charge density of naphthalene is significantly affected by the surrounding molecules in the lattice. By analyzing the change in electrostatic energy due to the induced moments interacting with a unit charge probe, we may indirectly observe the change in charge distribution caused by the crystalline environment. Figure 7 plots the change in the electrostatic energy using the SCERP induced moments, on the van der Waals plus 1.1 Å surface that is sampled by the atomic sites of the surrounding molecules. The anisotropic nature of the induction is clear. The increased electrostatic potential around the C(4)-H bond, in contrast to the C(2)-H bond, shows that the close contact with the π-electrons of the surrounding molecules in the crystal has significantly polarized this bond, as observed in the experimental charge density.21 3.3. 3-Azabicyclo[3,3,1]nonane-2,4-dione. 3-Azabicyclo[3,3,1]nonane-2,4-dione presents several challenges in terms of our polarizability calculations: the size of the molecule, in terms of basis functions required and associated computational limits, as well as the volume of space to be sampled for the point-to-point polarizabilities and the Cs symmetry in the molecule. Despite this, and the fact that symmetry of the molecules is not explicitly enforced by CamCASP at any stage, after refinement and localization the resulting polarizabilities are reassuringly symmetric. We find the induction energy for 3-azabicyclo[3,3,1]nonane-2,4-dione to be 33-36% of the electrostatic energy,


Figure 8. Induction energies for 3-azabicyclo[3,3,1]nonane2,4-dione. The crystal structures are ordered left-to-right by decreasing lattice stability, as calculated from the distributed static multipole + empirical repulsion-dispersion potential.

with good agreement between SCERP and the damped WSM models (Figure 8). For the crystal structures considered, the induction energy varies by less than 3 kJ mol-1, but this is significant relative to the difference in lattice energies of these structures calculated using a repulsion-dispersion model potential,64 which range from -95.08 to -97.64 kJ mol-1. Hence, more realistic modeling of the intermolecular interactions to include the induction energy would certainly rerank the structures. However, the observed hydrogen-bonded chain motif is not favored relative to many of the competitive dimer structures,39 and there is no clear-cut correlation with the hydrogen-bonding motif. Hence neglect of the induction energy does not appear to be the only problem in modeling the relative stability of crystal structures of 3-azabicyclo[3,3,1]nonane-2,4-dione.34 3.4. Carbamazepine. For carbamazepine, we contrast the SCERP with polarizabilities derived from fragment molecules (Figure 4). Despite the additional assumptions, there is still reasonably good agreement in the relative induction energies between SCERP and the damped L1 polarizability model, accounting for an increase in stability of 10.5-18.2 kJ mol-1 in the lattice energy. Both models find that the dimer-based structures, and particularly the experimental forms III and IV, are stabilized more by induction than the chain-based structures, and all hydrogen-bonded structures are stabilized more than the structure (ab41) with no hydrogen-bonding (Figure 9). This is significant, as the published crystal structure predictions33 for carbamazepine found that a structure with a hydrogen-bonded chain motif was more stable than the experimentally known dimer based structures. Improving the modeling of the electrostatic energies by using distributed multipoles from the better charge distribution used in the current work also alters the relative stabilities (Supporting Information), favoring the most stable observed polymorph form III. Hence, more accurate modeling of the electrostatics and adding the induction clearly gives a significant energy lowering to the most stable dimer based structures, which is in accord with experiment.


Figure 9. Induction energies for crystal structures of carbamazepine. The structures are ordered left-to-right by decreasing lattice stability, as calculated from the distributed multipoles described, plus an empirical64 repulsion-dispersion potential. The lattice-energy range for the structures shown is 16 kJ mol-1. The horizontal line indicates the average induction energy with the SCERP model to illustrate the discrimination of structural motifs by the polarizability model.

4. Discussion 4.1. How Important Is the Induction Energy for Organic Crystals? We have used two very different models for estimating the induced moments in organic crystals: an ab initio response to an applied field due to point charges representing the crystal environment, and the use of distributed polarizabilities in the field arising from a distributed multipole representation of the surrounding molecules. The induction energy contribution to the lattice energy, evaluated from these induced moments, is significant. Over this diverse range of crystal structures, the models agree that the induction energy is often between 20 and 40% of the electrostatic contribution to the lattice energy. This order of magnitude is consistent with estimates of the induction energy relative to the electrostatic energy for small polyatomic molecules,15,65-67 using equally rigorous or better models for the induction energy, although often the polarization is not iterated to self-consistency. It is also comparable with the less rigorous modeling of the induction contribution to lattice energies of neutral organic molecules derived by the pixel method,3,4,20 where experimental atomic polarizabilities are evenly distributed over the atomic charge density. More importantly, for different known and predicted crystal structures that otherwise have very similar lattice energies, the two models agree on the relative magnitude of the induction energy. In the case of oxalyl dihydrazide, inclusion of the intermolecular induction is essential for the calculated relative lattice energies to be consistent with the experimental observation of the polymorphs. This is an extreme case, as the intermolecular induction for the R polymorph compensates for the intramolecular hydrogen bonding in the other conformational polymorphs. In the case of carbamazepine, the induction energy favors the observed doubly hydrogen-bonded dimer based structures over the hypothetical catemer based structures. The differences in the induction energies of the low-energy computed structures of 3-azabicyclo[3,3,1]nonane-2,4-dione cannot be so simply ascribed to the hydrogen-bonding motif, but this reflects the


relative weakness of the hydrogen bonds for this imide, which forms a plastic phase.35 In each of these comparisons of known and hypothetical crystal structures, the differences in induction energies are small, only a few kilojoules per mole, but this is sufficient to provide a significant reordering of the relative stability of structures that are virtually equienergetic according to models which do not explicitly model the induction. To correctly model crystal lattice energy, intermolecular potentials require a reparametrization of the entire repulsiondispersion potential: adding the induction energy to lattice energies calculated using an empirically fitted potential involves a high degree of double counting. This is sufficient to lead to structures which are too dense if we attempt to minimize crystal structures with an induction term in addition to potentials which have been empirically fitted without the explicit inclusion of induction. It is also important that the model for the induction can be readily implemented in a program that minimizes lattice energies of organic crystal structure. 4.2. Practical Consideration for Using Polarizability Models in the Organic Solid State. A local polarizability model can be implemented in lattice energy minimization packages that use distributed multipole moments. It appears to be feasible to calculate WSM polarizabilities from a reasonable quality ab initio charge density for quite large molecules, with 3-azabicyclo[3,3,1]nonane-2,4-dione probably being the limit with current resources. This is an acceptable limitation, given that the transferable polarizability model calculated from fragments of carbamazepine gave reasonable results compared with the SCERP calculations that used the complete molecule. Thus, it seems that transferable polarizability models could be derived for use in modeling larger molecules. The induction energy does depend on the order of the polarizabilities included. We have noted some anomalous behavior where rank 2 polarizabilities are used on hydrogen, particularly when involved in short contacts within the crystal structure (most notably on oxalyl dihydrazide ). Given the small amount of charge density associated with polar hydrogen atoms, it seems reasonable that polarizabilities for these sites should be limited to rank 1 for applications to dense systems. The differences between L2 and L1 WSM models for the other atoms are comparable to those between them and the SCERP model. The error in modeling charge overlap effects in particularly short hydrogen-bonding geometries probably explains the larger variance with polarizability model observed in our oxalyl dihydrazide induction energies, relative to those for the 3-azabicyclo[3,3,1]nonane-2,4-dione crystals which do not have such short contacts. The WSM polarizability model does not account for any density overlap effects, and we have shown that damping is required in order to avoid unreasonable energies for the shorter intermolecular contacts found in hydrogen-bonded crystal structures. We have already shown26 the agreement between SAPT(DFT) induction energies and WSM models to be very good, and so the WSM polarizability method of modeling the induction energy has a firm foundation. This investigation has shown that damped

Welch et al.

polarizability models are also suitable for modeling the induction energy in large clusters that represent crystals, with many-body effects, qualitatively different field anisotropy and short contacts. We consistently find that the ab initio SCERP model falls midway between the L1 and L2/L1,WSM models and that the relative ordering of the energies is consistent. The SCERP model has also approximated the electrostatic field around the molecules (Figure 5), which does lead to significant errors in the hydrogen-bonding region. Thus, we conclude that we cannot at present model the induction energy more accurately than the range indicated by the differences between the SCERP and the L1 and L2/L1 WSM models. It is, however, clear (Figure 3) that the induced moments will need iterating to self-consistency.

5. Conclusions We have presented two distinct methodologies for approximating the effect of the different crystal environment on the charge distribution of four organic molecules, one using a self-consistent, ab initio derived induced multipoles approach and the other an ab initio derived, distributed, and localized polarizability (WSM) model. We have shown that these models can reproduce experimentally observed changes in the charge density when comparing gas phase and crystalline molecules in the case of naphthalene. We have also shown that the induction energy contribution to the lattice energy of organic molecules is significant and that properly describing the induction energy in lattice energy calculations may improve the relative ranking of the structures to be more in line with experimental observation. The WSM polarizability model and damping scheme can be extended from small polyatomics to crystal structure modeling, on the basis of its rigorous testing for smaller systems, and agreement with an alternative model in crystals (SCERP). A self-consistent WSM polarizability model would be a worthwhile addition to our ability to model molecular crystals. The considerable programming involved is in progress, alongside efforts to improve the theoretical basis of models for all terms in the intermolecular energy. Acknowledgment. This research was funded by EPSRC EP/C539109/1. The computed crystal structures used are held on STFC Data Portal. G.W.A.W. thanks Dr. Maurice Leslie for useful discussions; A.J.M. thanks Girton College, Cambridge for a Research Fellowship. Supporting Information Available: S1, structural information for polymorphs of oxalyl dihydrazide; S2, comments of transferability of polarizabilities; and S3, sensitivity of the relative lattice energies of carbamazepine to the electrostatic model. This material is available free of charge via the Internet at http://pubs.acs.org. References (1) Halgren, T. A.; Damm, W. Curr. Opin. Struct. Biol. 2001, 11, 236-242. (2) Fowler, P. W.; Stone, A. J. J. Phys. Chem. 1987, 91, 509511.

Induction Energy for Modeling Organic Crystals (3) Gavezzotti, A. J. Phys. Chem. B 2002, 106, 4145-4154. (4) Gavezzotti, A. J. Phys. Chem. B 2003, 107, 2344-2353. (5) Mooij, W. T. M.; van Eijck, B. P.; Kroon, J. J. Phys. Chem. A 1999, 103, 9883-9890.

J. Chem. Theory Comput., Vol. 4, No. 3, 2008 531 (33) Florence, A. J.; Johnston, A.; Price, S. L.; Nowell, H.; Kennedy, A. R.; Shankland, N. J. Pharm. Sci. 2006, 95, 1918-1930.

(6) Lindan, P. J. D.; Gillan, M. J. J. Phys.: Condens. Matter 1993, 5, 1019-1030.

(34) Hulme, A. T.; Johnston, A.; Florence, A. J.; Fernandes, P.; Shankland, K.; Bedford, C. T.; Welch, G. W. A.; Sadiq, G.; Haynes, D. A.; Motherwell, W. D. S.; Tocher, D. A.; Price, S. L. J. Am. Chem. Soc. 2007, 129, 3649-3657.

(7) Catlow, C. R. A.; Norgett, M. J. J. Phys. C Solid State 1973, 6, 1325-1339.

(35) Florence, A. J.; Leech, C. K.; Shankland, N.; Shankland, K.; Johnston, A. CrystEngComm 2006, 8, 746-747.

(8) Ren, P. Y.; Ponder, J. W. J. Comput. Chem. 2002, 23, 14971506.

(36) Karamertzanis, P. G.; Price, S. L. J. Chem. Theory Comput. 2006, 2, 1184-1199.

(9) Ren, P. Y.; Ponder, J. W. J. Phys. Chem. B 2003, 107, 59335947.

(37) Cabeza, A. J. C.; Day, G. M.; Motherwell, W. D. S.; Jones, W. Cryst. Growth Des. 2006, 6, 1858-1866.

(10) Thole, B. T. Chem. Phys. 1981, 59, 341-350. (11) Cox, S. R.; Hsu, L. Y.; Williams, D. E. Acta Crystallogr., Sect A 1981, 37, 293-301. (12) Williams, D. E. Acta Crystallogr., Sect. A: Found. Crystallogr. 1984, 40, C95. (13) Price, S. L. CrystEngComm 2004, 6, 344-353. (14) Price, S. L.; Price, L. S. Modelling Intermolecular Forces for Organic Crystal Structure Prediction. In Intermolecular Forces and Clusters I; Wales, D. J., Ed.; Springer-Verlag: Berlin, Heidelberg, Germany, 2005; pp 81-123. (15) Mooij, W. T. M.; van Duijneveldt, F. B.; van Duijneveldtvan de Rijdt, J. G. C. M.; van Eijck, B. P. J. Phys. Chem. A 1999, 103, 9872-9882. (16) Mooij, W. T. M.; van Eijck, B. P.; Kroon, J. J. Am. Chem. Soc. 2000, 122, 3500-3505.

(38) Cabeza, A. J. C.; Day, G. M.; Motherwell, W. D. S.; Jones, W. Cryst. Growth Des. 2007, 7, 100-107. (39) Motherwell, W. D. S.; Ammon, H. L.; Dunitz, J. D.; Dzyabchenko, A.; Erk, P.; Gavezzotti, A.; Hofmann, D. W. M.; Leusen, F. J. J.; Lommerse, J. P. M.; Mooij, W. T. M.; Price, S. L.; Scheraga, H.; Schweizer, B.; Schmidt, M. U.; van Eijck, B. P.; Verwer, P.; Williams, D. E. Acta Crystallogr., Sect. B: Struct. Sci. 2002, 58, 647-661. (40) Karamertzanis, P. G.; Day, G. M.; Welch, G. W. A.; Kendrick, J.; Leusen, F. J. J.; Neumann, M. A.; Price, S. L. J. Chem. Phys. 2008, submitted. (41) Chisholm, J. A.; Motherwell, S. J. Appl. Crystallogr. 2005, 38, 228-231. (42) Howie, R. A.; Skakle, J. M. S. Acta Crystallogr., Sect. E 2001, 57, o822-o824.

(17) Gascon, J. A.; Leung, S. S. F.; Batista, E. R.; Batista, V. S. J. Chem. Theory Comput. 2006, 2, 175-186.

(43) Angyan, J. G.; Jansen, G.; Loos, M.; Hattig, C.; Hess, B. A. Chem. Phys. Lett. 1994, 219, 267-273.

(18) Gavezzotti, A. J. Chem. Theory Comput. 2005, 1, 834-840.

(44) Stone, A. J. Mol. Phys. 1985, 56, 1065-1082.

(19) Gavezzotti, A. Struct. Chem. 2005, 16, 177-185.

(45) Le Sueur, C. R.; Stone, A. J. Mol. Phys. 1993, 78, 12671291.

(20) Gavezzotti, A. Z. Kristallogr. 2005, 220, 499-510. (21) Oddershede, J.; Larsen, S. J. Phys. Chem. A 2004, 108, 1057-1063. (22) Misquitta, A. J.; Jeziorski, B.; Szalewicz, K. Phys. ReV. Lett. 2003, 91, art. no.-033201. (23) Misquitta, A. J.; Szalewicz, K. J. Chem. Phys. 2005, 122, art-214103. (24) Misquitta, A. J.; Podeszwa, R.; Jeziorski, B.; Szalewicz, K. J. Chem. Phys. 2005, 123, 214103. (25) Misquitta, A. J.; Stone, A. J. J. Chem. Theory Comput. 2008, 4 (1), 7-18. (26) Misquitta, A. J.; Stone, A. J.; Price, S. L. J. Chem. Theory Comput. 2008, 4 (1), 19-32.

(46) Misquitta, A. J.; Stone, A. J. J. Chem. Phys. 2006, 124, 024111. (47) Le Sueur, C. R.; Stone, A. J. Mol. Phys. 1994, 83, 293307. (48) Williams, G. J.; Stone, A. J. J. Chem. Phys. 2003, 119, 4620-4628. (49) Bondi, A. J. Phys. Chem. 1964, 68, 441-451. (50) Misquitta, A. J.; Stone, A. J. CamCASP: A program for studying intermolecular interactions and for the calculation of molecular properties in distributed form, Version 5; University of Cambridge: 2007. See: http://www-stone.ch. cam.ac.uk/programs.html#CamCASP (accessed Feb 2008).

(27) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361-373.

(51) DALTON, a molecular structure program, release 2.0; 2005. See: http://www.kjemi.uio.no/software/dalton/dalton.html (accessed June 21, 2007).

(28) Willock, D. J.; Price, S. L.; Leslie, M.; Catlow, C. R. A. J. Comput. Chem. 1995, 16, 628-647.

(52) Allen, M. J.; Tozer, D. J. J. Chem. Phys. 2000, 113, 51855192.

(29) DMAREL, Version 4.1.1; Price, S. L.; Willock, D. J.; Leslie, M.; Day, G. M. 2004.

(53) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 61586170.

(30) Gavezzotti, A. CrystEngComm 2002, 4, 343-347.

(54) Sadlej, A. J. Collect. Czech. Chem. C 1988, 53, 1995-2016.

(31) Price, S. L. AdV. Drug DeliVery ReV. 2004, 56, 301-319.

(55) Stone, A. J. J. Chem. Theory Comput. 2005, 1, 1128-1132.

(32) Ahn, S. Y.; Guo, F.; Kariuki, B. M.; Harris, K. D. M. J. Am. Chem. Soc. 2006, 128, 8441-8452.

(56) Stone, A. J. The Theory of Intermolecular Forces, 1st ed.; Clarendon Press: Oxford, 1996.

532 J. Chem. Theory Comput., Vol. 4, No. 3, 2008 (57) Stone, A. J.; Misquitta, A. J. Int. ReV. Phys. Chem. 2007, 26, 193-222. (58) Tang, K. T.; Toennies, J. P. J. Chem. Phys. 1984, 80, 37263741. (59) Stone, A. J.; Dullweber, A.; Engkvist, O.; Fraschini, E.; Hodges, M. P.; Meredith, A. W.; Nutt, D. R.; Popelier, P. L. A.; Wales, D. J. Orient: a program for studying interactions between molecules, Version 4.6; University of Cambridge: 2006. See: http//www-stone.ch.cam.ac.uk/ programs.html#Orient (accessed Feb 2008). (60) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Chem. ReV. 1994, 94, 1887-1930. (61) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, 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.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.;

Welch et al. 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. Gaussian 03; Gaussian Inc.: Wallingford, CT, 2003. (62) Buckingham, A. D.; Fowler, P. W. Can. J. Chem. 1985, 63, 2018-2025. (63) Bernstein, J. Polymorphism in Molecular Crystals; Clarendon Press: Oxford, 2002. (64) Coombes, D. S.; Price, S. L.; Willock, D. J.; Leslie, M. J. Phys. Chem. 1996, 100, 7352-7360. (65) Chipot, C.; Luque, F. J. Chem. Phys. Lett. 2000, 332, 190198. (66) Chipot, C.; Angyan, J. G. New J. Chem. 2005, 29, 411420. (67) Jansen, G.; Hattig, C.; Hess, B. A.; Angyan, J. G. Mol. Phys. 1996, 88, 69-92. CT700270D

Is the Induction Energy Important for Modeling Organic Crystals

Recommend Documents