Comparison of Static and Fluctuating Charge Models for Force-Field

Crystal Growth & Design , 2005, 5 (3), pp 925–933. DOI: 10.1021/ ... Publication Date (Web): February 9, 2005 ... Cite this:Crystal Growth & Design ...
3 downloads 0 Views 152KB Size
Download PDF
Comparison of Static and Fluctuating Charge Models for Force-Field Methods Applied to Organic Crystals Sven

Brodersen,*,†

Steffen

Wilke,†

Frank J. J.

Leusen,‡

and Gerhard E.

Engel†

Accelrys Ltd., 334 Cambridge Science Park, Cambridge, CB4 0WN, United Kingdom, and Institute of Pharmaceutical Innovation, University of Bradford, Bradford, BD7 1DP, United Kingdom Received September 27, 2004;

CRYSTAL GROWTH & DESIGN 2005 VOL. 5, NO. 3 925-933

Revised Manuscript Received December 15, 2004

ABSTRACT: We investigated different charge models in a force-field approach for calculating structural data of organic crystals. Charges were calculated according to the Gasteiger and charge equilibration (QEq) methods or were derived from the quantum-mechanical electrostatic potential (ESP). A fluctuating charge model, based on the QEq approach, for describing polarization effects was developed. Different variants of the fluctuating charge model as well as static charges were applied to a test set of 48 organic crystals using the Dreiding force field. On average, the inclusion of crystal field effects with the fluctuating charge model lead to slight improvements over static ESP charges, both in terms of the average errorswith respect to experimentsand in terms of the standard deviation. The empirical charge models, Gasteiger and QEq, tend to be less accurate if compared with quantum-mechanically derived charges. Introduction Many atomistic simulation tasks are too demanding to be tackled with accurate quantum mechanical (QM) methods. This is the case for large systems involving many atoms, for long simulation times in molecular dynamics, or when large numbers of calculations need to be carried out on similar systems, e.g., in polymorph prediction. Therefore, there is a need for fast force-field based methods, which of course have to meet required accuracy criteria. In the case of polymorph prediction, crystal structures and energies for typical drug-sized molecules have to be calculated very accurately.1,2 Intermolecular interactions play a crucial role in these calculations. For rigid molecules, they directly determine the relative positions and orientations of the molecules. For flexible molecules, inter- and intramolecular energies together determine the molecular conformations that strongly influence crystal packing. In contrast to van der Waals interactions, where reliable and sufficiently accurate descriptions can be found, in classical simulations it is notoriously difficult to model electrostatic interactions with high accuracy. Since the electrostatic contribution is usually the most influential part of the intramolecular interaction, an accurate representation of the charge density, and thus the electrostatic interactions, is essential for the accurate simulation of crystal structures.3 As has been shown in earlier investigations, a good reproduction of the molecular electrostatic potential (ESP) in a vacuum is a solid basis for calculating the electrostatic interactions in the crystal environment.4 This can be achieved, for instance, with multipoles centered at the atomic sites.5,6 In comparison with ESP-derived charges, multipoles lead to a considerable improvement in calculated lattice parameters for rigid molecules. It is not yet clear * To whom correspondence [email protected]. † Accelrys Ltd. ‡ University of Bradford.

should

be

addressed.

E-mail:

how to treat molecular flexibility with multipoles, as witnessed by our previous attempt which led to large errors in calculated lattice parameters for some organic crystals.7 A method for describing flexible molecules in the gas phase with multipoles including polarization energy is given in ref 8. In this investigation, the electrostatic interactions are calculated with point charges which are less affected by molecular flexibility. We compare different charge methods by calculating lattice parameters for a set of 48 organic crystals. Besides static charges calculated with the Gasteiger method,9 the charge equilibration (QEq) approach,10 and from the ESP, we investigate the performance of a newly developed fluctuating charge model. This model aims to capture polarization effects due to the crystal field and conformational changes. Various fluctuating charge models have been shown to be useful in the description of intermolecular interactions11-14 and molecular conformation.15 There are several approaches to describe polarization with atomic dipoles.16-19 In comparison with these models, fluctuating charge models are computationally less demanding while still capable of capturing a substantial part of the polarization energy.20 In the next two sections, the fluctuating charge model and the optimization of its parameters are described. Then, we describe the procedure and the results of the crystal structure optimizations. Consistent QEq Atomic charges derived from the ESP of a molecule are known to be well suited for calculating the Coulomb part of the force-field energy. As they cannot be derived in the crystalline environment, a method for incorporating polarization due to the crystal field is desirable. The aim of our fluctuating charge model is to describe induction effects due to a crystal field, conformational changes, or an external field. The total energy expression describes not only the energy change due to changes

10.1021/cg049673+ CCC: $30.25 © 2005 American Chemical Society Published on Web 02/09/2005

926

Crystal Growth & Design, Vol. 5, No. 3, 2005

in the atomic coordinates but also the polarization contribution due to changes in the atomic partial charges. Once such an energy expression is derived, both the optimal geometry of the atoms and the induced charges for each such geometry can be consistently obtained by minimizing the same total energy expression with respect to both coordinates and induced charges. It is worth noting that within this method different molecular conformations will generally have different partial charge distributions. Because both geometries and partial charge distributions are derived from the same energy expression, the energy differences between these conformations are still meaningful. Therefore, within this approach it becomes possible for polymorph predictions to abandon the need to use the same partial charges for all conformations of a molecule and for all crystal packings. We represented the atomic charges as a sum of reference charges q0i, usually ESP-derived charges for a molecule, and induced charges δqi:

qi ) q0i + δqi The induced charges were determined by minimizing an energy expression similar to the QEq model by Rappe´ and Goddard,10,21 which essentially is a second-order expansion of the energy with respect to the partial charges:

ECoul )

1

1

qiJijqj + ∑ Jiiqi3 + ∑φiqi ∑i χiqi + 2∑ 2ξi∈H i,j i

with χi being atomic electronegativities, Jij a shielded Coulomb interaction between atom i and atom j, ξ a parameter for hydrogen atoms introduced by Rappe´ and Goddard, and φi an external electric field. In the original QEq approach, χi and Jii are taken from experimental data, and Jij is calculated from the overlap of s-type Slater functions u(r) between atoms centered at ri and rj:

Jij )

∫∫

2

2

|ui(r bi - b r 1)| |uj(r bj - b r 2)| 3 d r1 d3r2 |r b1 - b r 2| u(r) ) Nrn-1 exp(-Rr)

Here, N denotes the normalization constant, n is the principal quantum number of the highest occupied orbital, and R is a damping factor related to the atomic radius R by the expression

R)

2n + 1 4R

Nishimoto and Mataga22 gave another, simpler, approximation of the shielded Coulomb interaction:

Jij )

Jii + Jjj (Jii + Jjj)|r bi - b r j| + 2

which had just one parameter, the on-site term Jii. A comparison of these functions for the CO molecule is given in Figure 1. For large distances r, both the Slater and the Nishimoto-Mataga functions follow the un-

Brodersen et al.

Figure 1. Comparison of the shielded Coulomb interaction according to Slater and Nishimoto-Mataga, and the unshielded Coulomb interaction 1/r (“Coulomb”) for the CO molecule.

shielded Coulomb energy 1/r, but they have a finite value at r ) 0. The original QEq theory determines the total value of partial atomic charges by minimizing the Coulomb energy expression above using electronegativities related to free atom values. In our approach, we started from the assumption that we have a fairly accurate force-field energy expression for the system in question. This force field uses an arbitrary set of fixed atomic partial charges q0i, for instance, ESP charges. These partial charges have been optimized for a given reference configuration of the system. Instead of using the Coulomb energy expression above to determine the total value of the partial charges, we required that the energy ECoul has a stationary point for the reference configuration (which enters via the distance-dependent Coulomb interaction matrix) at the reference charges. This requirement fixed the electronegativities to be

∑j

χi ) -

3Jiiq0i2 Jij q0j δiH 2ξ 0

with δiH equals one if atom i is a hydrogen atom and zero otherwise, and Jij0 being the Coulomb interaction matrix for the reference configuration. Thus, the induced charges were zero for the reference configuration, which was usually the free molecule for which ESP (i.e., reference) charges were calculated. Essentially, in the current approach the QEq expression is used to determine partial charge differences with respect to a reference configuration of the system. At each step of the force-field minimization, the Coulomb energy ECoul was minimized with respect to δqi, subject to charge conservation, and added to the forcefield energy. Therefore, the total energy not only contained the usual Coulomb energy, but also the energy required to change the charges. So, the method is dubbed “consistent” in contrast to methods that either use static charges or neglect the polarization energy. The quadratic expansion of the polarization energy was not without problems. If we neglect for the moment the cubic dependence of ECoul on the hydrogen charges, the energy is quadratic to the charges. It has a minimum if (and only if) the matrix J is positive definite.

Comparison of Static and Fluctuating Charge Models

Depending on the chosen approximation to Jij and the atomic positions, this requirement may not hold. In this case the analytical solution, obtained by inverting J, yields the saddle point while a numerical minimization diverges. Changing J in a way that leads to positive eigenvalues, e.g., by increasing the diagonal elements, can circumvent this problem. Another possibility is to diagonalize J, replace negative eigenvalues with a small positive value, and transform back the matrix. This procedure is time-consuming and changes every entry of J, although usually by a very small amount. If we allow for cubic contributions from hydrogen charges, the energy has no lower bound but can have a local minimum. A numerical minimization may, therefore, diverge or converge depending on the starting point, even if J is positive definite. Convergence can be forced by adding a function to ECoul that increases sufficiently quickly with the induced charges. We have chosen to add a quartic function in δqi, which we call a penalty function:

Epenalty )

( ) δqi

∑i δq

4

pen

With δqpen chosen to be 0.3, this function effectively suppresses large δqi and ensures the existence of a minimum without changing the Coulomb interaction matrix J. The penalty function describes the fact that the total energy will increase considerably if induced charges are large and the charge distribution is disturbed significantly. The minimization of ECoul, including nonquadratic contributions, was carried out with a conjugate gradient algorithm. Optimization of Parameters The parameters entering the Coulomb interaction matrix Jij with Slater-type interactions had to be optimized as the original parameters given in ref 10 did not lead to an improvement in lattice parameters. To this end, induced ESP charges were calculated by perturbing the molecule in a vacuum with an external field. The difference in ESP charges relative to the zerofield calculation defined the induced ESP charges, δqESP. The external fields were produced by randomly placing a dipole with charge 1.9675 or 3.148 e and length 0.5 aB at the molecular surface (at least 2.5 aB away from any atom and at most 3.0 aB away from the nearest atom). The high charge of 3.148 e was necessary in some cases to produce a measurable effect on the ESP charges. The results for 20-30 dipole fields served as “target” values for the modeled induced charges calculated by minimizing ECoul. The deviation of the QMderived induced charges, δqESP, from the modeled induced charges, δq, was measured with the fit function

χ)

ESP (δqφ,i - δqφ,i)2 ∑ φ,i

with φ denoting the external fields and i the atoms. This function was minimized with respect to the atomic radii Ri and the Coulomb self-energy terms Jii. As derivatives of χ with respect to the parameters Ri and Jii cannot be evaluated, the minimizations were carried out with Powells method,23 which uses function values only.

Crystal Growth & Design, Vol. 5, No. 3, 2005 927

Parameters for atoms with the same atom types were forced to be identical, thereby greatly reducing the number of parameters to be fitted. Still, the function χ showed insensitivity to variation of some parameters or their combination, making the fitting procedure rather unstable. The resulting excessive wandering in parameter space was suppressed by adding a confining term

χconf )

∑t ct(Vt - Vt0)2

to χ. Here, t indexes the variables V (either Ri or Jii), with V0 being the original parameters given in ref 10. The coefficients c of these parabolas were chosen to be very small, typically 0.001, not to disturb the minimization procedure too much. In this way, parameters that had no influence on the fitting quality kept their original value. Parameters that had a pronounced effect on χ were hardly affected by this confining term. In general, induced charges calculated with the original parameters were too large. The fitting procedure led to decreased atomic radii and increased on-site terms. Both these effects resulted in smaller induced charges. For some of the investigated molecules, the fitting led to small improvements only, even when all atoms (instead of atom types) were allowed to have different parameters. Together with the unstable minimization of the fit function due to redundant parameters, this indicated that the model itself did not adequately describe the physical reality for these molecules. This could be due to a number of approximations: the pointcharge approximation of a continuous charge density, the quadratic (cubic for hydrogen) dependence of ECoul on the charges, and the approximation of the Coulomb interaction. The last approximation could be overcome by directly fitting the matrix elements Jij24,25 instead of fitting parameters in the model expression. We did not try this as it increased the complexity of the fitting procedure and provided less insight into the underlying physics. Crystal Structure Calculations The 48 validation examples are the same as in our previous investigation7 (see Figure 2 for their chemical structures). They comprise a variety of pharmaceutically relevant molecules, such as steroids, naphthalenes, substituted benzenes, salts, and zwitterions, containing the elements carbon, hydrogen, oxygen, nitrogen, sulfur, phosphorus, fluorine, chlorine, bromine, and iodine. The test set contains polar and nonpolar, rigid and flexible, large and small compounds, and a wide range of functional groups. The set can be roughly divided into five classes as indicated in Table 1. The first class contains more or less rigid molecules consisting of hydrogen, carbon, oxygen, and nitrogen only. The more flexible molecules of the second class additionally comprise halogens, sulfur, and phosphorus. The third class allows inclusion of water molecules in the crystal structure. The fourth class contains salts with small but very flexible molecules. In two cases, water molecules are included. The fifth class consists of zwitterions. The outline of the investigation is as follows: for each molecule contained in the crystal, a QM geometry

928

Crystal Growth & Design, Vol. 5, No. 3, 2005

Brodersen et al.

Figure 2. Molecular structures of the 48 crystals investigated, labeled with their CSD reference code.

optimization in a vacuum was done. The crystal was rebuilt with the resulting molecules and optimized with a molecular mechanics method, using different methods for the electrostatic energy. The calculated crystal parameters and the molecular position(s) were compared with experiment. The QM calculations were performed in a vacuum with DMol3,26-28 using density functional theory with a nonlocal exchange-correlation functional. These calcula-

tions, of which details have been given elsewhere,7 served two purposes: first, the resulting molecular geometries were used as the starting points in the forcefield calculations. Second, atomic charges were fitted to the QM ESP using the scheme proposed by Hinsen and Roux,29 which proved to be very stable. Molecular mechanics calculations on the crystals were performed with the Dreiding30 force field as implemented in the Materials Studio28 modeling software. We

Comparison of Static and Fluctuating Charge Models

Crystal Growth & Design, Vol. 5, No. 3, 2005 929

Table 1. Results of Crystal Structure Optimizations for 48 Substances Using Nine Different Electrostatic Modelsa ref code

class

DMol3

No Charges

Gasteiger

QEq (Mol)

QEq (Cry)

ESP

QEq (Mol)

ESP

ESP

ESP

HACCIX HACNUU LABHAX LABHAX01 VOXJEX WASROX WERWUL DETBAA03 DIWWEL DIWWEL01 DIWWEL02 ACSALA AMSALA01 IBPRAC JEKNOC ACMDOB ACNPHB BANGOM01 BANGOM02 BETPIZ CICVOZ JECYIZ CENRIW DUPTOX DUPTOX01 EABMEZ EABSIJ EACBOZ HAXBUD HECNAE JAWJIA JIGFUA KUHPAE CHAPEP01 CHAPEP02 GEBMEF LEKRIC LEKROI PEAMAN01 XIMBZB XIMBZB01 HECFIE JILTUT BEYZIO BEYZIO01 CAXMOD CAXMOD01 DIPTOX10

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5

0.106 0.099 0.064 0.333 0.077 0.353 0.069 0.091 0.115 0.208 0.097 0.057 0.040 0.274 0.224 0.230 0.063 0.077 0.075 0.321 0.058 0.346 0.112 0.146 0.126 0.174 0.131 0.341 0.149 0.132 0.256 0.140 0.243 0.397 0.097 0.174 0.237 0.288 0.117 0.125 0.100 0.434 0.409 0.037 0.047 0.052 0.155 0.062

Static 47 30 42 324 28 51 16 15 53 24 367 77 177 520 128 36 724 29 285 30 135 36 32 274 22 27 33 548 543 22 414 60 82 275 188 69 102 211 431 158 41 46 74 38 65 77 2491 406

Static 50 27 36 295 28 51 14 12 56 78 379 42 179 528 125 28 367 23 125 105 374 283 41 169 18 23 225 194 233 14 110 51 73 108 262 1482 46 538 23 468 178 123 31 57 65 115 822 206

Static 129 22 74 243 35 40 10 16 70 20 367 14 253 556 135 18 1959 385 1266 96 199 62 44 21 15 15 43 284 168 10 47 41 106 313 263 337 85 207 76 25 26 72 767 70 59 172 1511 1766

Static 116 20 37 192 26 32 6 19 120 14 327 11 191 521 129 14 1369 375 858 23 53 50 30 22 9 13 155 739 174 8 36 30 86 581 181 34 35 181 70 142 16 79 76 46 50 103 677 172

Static 56 20 39 253 32 44 9 18 57 26 424 11 328 539 140 46 805 20 236 59 173 78 134 18 17 53 60 110 12 14 186 32 60 128 224 87 30 242 40 261 48 83 33 27 62 60 461 125

Slater 119 20 68 250 27 68 10 15 96 15 353 11 58 513 120 58 645 532 997 19 123 59 47 20 13 12 162 812 178 11 20 38 91 284 202 48 58 178 53 112 20 61 141 52 57 154 2908 2098

Nish-Mat 48 16 33 116 29 30 7 18 61 16 409 11 434 526 126 49 835 18 203 64 34 60 138 168 20 59 35 77 145 10 93 32 79 182 186 49 40 364 144 347 156 63 94 29 48 72 445 116

Slater 48 17 35 278 28 47 10 15 62 15 427 12 414 522 121 57 855 17 219 30 37 58 133 20 18 52 38 63 134 12 223 39 82 177 176 18 31 343 118 311 130 92 84 30 55 76 2728 100

Slater, fit 47 17 34 277 25 62 11 13 52 16 415 7 58 521 120 60 853 18 243 14 124 56 111 17 15 40 33 49 163 14 213 41 77 182 180 24 20 165 31 103 79 76 53 29 58 88 3178 98

a CSD reference codes (first column), class numbers (second column), rms deviations in Å of QM-optimized molecules relative to experiment (third column), and structural drift factors F (fourth to 12th column) for the nine electrostatic models investigated.

used default parameters throughout with a LennardJones type potential for the van der Waals energy and explicit hydrogen bonding with a 12-10 Lennard-Jones term. All nonbonded interactions were summed with the Ewald technique. Full optimization of crystal parameters and atomic positions were carried out with a series of minimization algorithms. We investigated nine variants for calculating the electrostatic energy: (1) No charges (2) Static Gasteiger charges (3) Static QEq charges, calculated for the free molecule (“QEq(mol)-static”) (4) Static QEq charges, calculated for the crystal (“QEq(cry)-static”) (5) Static ESP charges (“ESP-static”) (6) Fluctuating QEq charges, calculated for the molecule, with Slater-type interaction (“QEq(mol)-Slater”)

(7) Fluctuating ESP charges with Nishimoto-Matagatype interaction matrix (“ESP-Nish-Mat”) (8) Fluctuating ESP charges with Slater-type interaction matrix (“ESP-Slater”) (9) Fluctuating ESP charges with Slater-type interaction matrix and optimized parameters (“ESP-Slater-fit”) In the first case (“No Charges”), no electrostatic interaction was calculated. In the next four methods, with the attribute “static”, the charges were fixed at their starting values. The last four methods used the fluctuating charge model presented in the former section but started with different reference charges (ESP or QEq) or used different approximations to the interaction matrix Jij (Nishimoto-Mataga or Slater-type approximation with original or fitted parameters). For the fluctuating charge models, Jij was calculated for 1-2 and 1-3 neighbors only. Electrostatic interactions of 1-4 neighbors and beyond were treated with the unshielded Coulomb form.

930

Crystal Growth & Design, Vol. 5, No. 3, 2005

Brodersen et al.

Table 2. Average Results and Standard Deviations for the Full Set of Structures and a Reduced Set as Described in the Text ref charges

No Charges

Gasteiger

QEq (Cry)

ESP

ESP

ESP

ESP

interaction full set (48) average 〈F〉 standard deviation σF outliers 〈F〉 (without outliers) σF (without outliers) reduced set (36) average 〈F〉 standard deviation σF outliers 〈F〉 (without outliers) σF (without outliers)

Static

Static

QEq (Mol) Static

Static

Static

Slater

Nish Mat

Slater

Slater, fit

206 380 CAX1 158 178

185 256 GEB 157 172

261 450 ACN,DIP 191 303

172 269 ACN 146 205

125 159 ACN 111 125

250 527 CAX1,DIP 152 217

132 165 ACN 117 130

179 409 CAX1 125 163

170 468 CAX1 106 152

102 135 EAC 90 113

129 247 GEB 90 86

153 313 DIP 107 149

81 134 EAC 62 73

64 62 AMS 56 44

146 368 DIP 91 155

68 76 AMS 58 43

64 75 AMS 54 45

49 41 JAW 45 30

The quality of the calculated crystal structures were assessed with the so-called structural drift factor F:31

F ) (∆φ/2)2 + (∆x/0.1 Å)2 + (100∆a/a)2 + (100∆b/b)2 + (100∆c/c)2 + ∆R2 + ∆β2 + ∆γ2 with ∆φ denoting the net rotation of the molecule with respect to the experimental orientation, ∆x is the translation (in Å) of the molecule with respect to the experimental position, and ∆a, ∆b, ∆c, ∆R, ∆β, ∆γ are the differences between the calculated and experimental lattice parameters (in Å and degrees). Typical values for F were between 10 (very close to experiment) and 1000 (far off experiment). As in all our calculations temperature effects were neglected, and absolute agreement of lattice parameters cannot be expected. Thus, values of F below, say, 20 should be considered as reliable. Results and Discussion The results of the crystal structure calculations with the nine electrostatic methods for all 48 substances are presented in Table 1. The first two columns show the CSD (Cambridge Structural Database)32 reference codes of the crystals and the class numbers as mentioned in the preceding section. The third column contains the rms deviations in atomic positions of the QM-optimized molecules relative to experiment. To take into account poor accuracy in the experimental hydrogen positions, distances were weighted with the atomic mass. If there is more than one molecule per asymmetric unit, the highest rms value is given. The deviations from experiment are not (only) due to inaccuracies of the QM calculations but also reveal the influences of the crystal fields on the molecular conformations. As such, they are an indication of the flexibility of the molecules as well as the strengths of the crystal fields. The following nine columns in Table 1 show the structural drift factors F for the nine electrostatic models. In Table 2, we show the average values 〈F〉 and standard deviations σF for the 48 samples. Because in a few cases F has very large values, we also show average and standard deviations without these outliers. They were defined as having F-values higher than 〈F〉 + 3σF and are indicated in the row labeled “Outliers”. In the following, results are discussed in terms of the values without outliers unless noted otherwise. There are a few structures whose lattice parameters were not well reproduced with any of the charge models.

QEq (Mol)

This indicates that the Dreiding force field is inadequate for these structures. Thus, in Table 2 we also show average values for all structures that have F-values below 100 for at least two electrostatic methods. This selection is arbitrary, but the definition of which structures can be described adequately by Dreiding does not have a major effect on the conclusions. According to this definition, the following structures can be disregarded: LABHAX01, DIWWEL02, IBPRAC, JEKNOC, ACNPHB, BANGOM02, HAXBUD, CHAPEP01, CHAPEP02, LEKROI, XIMBZB, and CAXMOD01. The results for the reduced set of 36 substances are given in the lower five rows of Table 2. The Dreiding force field was originally developed with either no charges or Gasteiger charges.30 Our results show that the two approaches led to almost the same average values (1578 and 157, respectively) of the structural drift factor. However, individual results, see Table 1, were similar for about half of the test set. Considerably different results were obtained for the 10 salts, four out of five zwitterions and about 10 others, mainly from class 2. Surprisingly, the crystal structures of some salts were better described without charges than with Gasteiger charges, and in two cases (XIMBZB01, HECFIE) the F-values are even very low (41 and 46, respectively). The very poor result for GEBMEF with Gasteiger charges is due to a large rotation of the anion with respect to its experimental orientation, which has a strong effect on the lattice parameters. Like the Gasteiger method, the QEq method offers a fast, empirical charge model. Unlike Gasteiger charges, QEq charges account for the conformations and environments of molecules. This allows us to estimate the effect of the crystal field by comparing the results obtained with QEq charges calculated for the molecules in a vacuum to the results obtained with QEq charges calculated for the molecules in their crystal environments. The lattice parameters obtained with these two methods show a significant improvement when the effect of the crystal field is included: the average F-value drops from 191 (“QEq(mol)-static”) to 146 (“QEq(cry)-static”). The consistency of the improvement is noteworthy: only five structures (DIWWEL, EABSIJ, EACBOZ, CHAPEP01, and XIMBZB) were described considerably worse when the crystal environment was taken into account and over 20 structures (e.g., the zwitterions) improved noticeably. The usage of QM-derived ESP charges causes the average F-value to drop further to 111 even though these charges were derived for the molecules in a

Comparison of Static and Fluctuating Charge Models

vacuum and ignore crystal field effects. This is not so much due to an improvement of the majority of structures in comparison with QEq(cry)-charges (less than half of them have lower F-values with ESP charges). Instead, some structures improved a lot and only very few structures (such as IBPRAC and ACNPHP) have high F-values. This led to low standard deviations (159 for all results, 125 without outliers), which can be regarded as a measure for the reliability of this method. The goal of our fluctuating charge model is to combine the accuracy of ESP charges with inclusion of the crystal field. However, our approach proved to be only partially successful. Taking molecular QEq charges as a reference point, the fluctuating charge model improved the results and yielded an average F-value (152) close to the value obtained with static crystalline QEq charges (146). We take this as a sign that the crystal field effects were by and large captured correctly. The salts especially (except XIMBZB) improved by adding polarization charges. Considering the results obtained with ESP charges, the fluctuating charge model did not significantly improve over static charges regardless of which model was used to shield the Coulomb interactions. With the Nishimoto-Mataga approximation, averages and standard deviations (with and without outliers) are close to the results obtained with static ESP charges. A slight improvement can be observed for all steroids (except KUHPAE, which includes a solvent molecule), while 7 out of 10 salts were described considerably worse with the fluctuating charge model. The more sophisticated Slater approximation to the shielded Coulomb interaction yielded slightly worse results (125) than the Nishimoto-Mataga approximation. Both approaches gave similar figures for most structures with only few exceptions, such as LABHAX01, BETPIZ, DUPTOX, JAWJIA, GEBMEF, and, most strikingly, CAXMOD01. The results obtained with the fluctuating charge model using the Slater approximation with optimized parameters (106) were better than the Slater results without optimized parameters (125), although they were only marginally better than the results with static ESP charges (111). Comparison to the results obtained with original parameters (i.e., “ESP-Slater”) shows that for many structures the lattice parameters hardly changed. In one case, CICVOZ, it even worsened the agreement with experiment. Noticeable improvement was achieved for AMSALA01 and most of the salts. Considering the reduced set of structures (Table 2, lower part) confirms that the fluctuating charge model with original parameters yielded approximately the same results as static ESP charges. With fitted parameters, the improvement over static charges was more convincing: 〈F〉 drops from 56 to 45. In Tables 4 and 5, we show the value of atomic charges for the two molecules DIWWEL and JECYIZ, respectively. The first column indicates the atoms as labeled in Figure 3. Hydrogen atoms that are not shown in the chemical structure inherit their name from the carbon atom they are bonded to, and multiple hydrogen atoms get an additional index a, b, c, etc. In the second and third column, the ESP charges before and after the crystal structure optimization (with Slater-type interaction) are shown. The difference in the charges are very

Crystal Growth & Design, Vol. 5, No. 3, 2005 931 Table 3. Energy Differences of Polymorphs in kcal/mol Per Asymmetric Unit Calculated with Static and Fluctuating ESP Chargesa name

ESP static

ESP Slater

LAB-LAB1 IBP-JEK BAN1-BAN2 DUP-DUP1 CHA1-CHA2 XIM-XIM1 BEY-BEY1 CAX-CAX1 DIW-DIW1 DIW-DIW2

2.156 -1.625 1.082 -1.547 21.529 -1.943 -0.437 -3.617 0.181 -3.043

4.708 -2.217 3.688 -1.275 18.096 4.121 -0.849 -5.537 -0.004 -4.405

a In the first column, an abbreviation of the CSD code is given. For example, LAB-LAB1 means energy (LABHAX) - energy (LAHAX01).

Table 4. Comparison of Charges for DIWWEL C1 H1 C2 H2 C3 H3 C4 C5 H5 C6 H6 C7 H7 C8 C9 C10 C11 H11a H11b H11c C12 H12a H12b H12c N13 N14 O14a O14b F

ESP static

ESP Slater

-0.312 0.174 -0.053 0.151 -0.452 0.201 0.468 -0.387 0.185 -0.041 0.148 -0.335 0.222 0.135 -0.255 0.264 -0.518 0.176 0.145 0.183 -0.414 0.151 0.121 0.171 0.051 0.582 -0.368 -0.392 57

-0.321 0.193 -0.069 0.167 -0.464 0.220 0.462 -0.372 0.208 -0.034 0.159 -0.340 0.236 0.124 -0.267 0.256 -0.509 0.165 0.165 0.181 -0.412 0.161 0.129 0.161 0.052 0.565 -0.397 -0.418 62

small, at most 0.03 au. For DIWWEL, most charges become larger in value during the optimization, whereas for JECYIZ the modulus of the charges kept similar on average. Thus, it seems that either the polarization model is not sensitive enough to the crystal field or the ESP charges already include some polarization effect. As the value of the ESP charges depend on the underlying QM method and its parameters, e.g., the basis set, it can be that the ESP charges lead to too high dipole moments. This fact is exploited in the AMBER force field.33 However, we did not investigate the dependence of the crystal structure optimization on the QM method. The very bad results for the zwitterions CAXMOD01 (with all methods) and DIPTOX10 (with molecular QEq charges) are due to the flexibility of the molecules which led the lattice optimization to change the hydrogen bonding patterns. For CAXMOD, the hydrogen bonding pattern was described correctly in most cases. The poor results for ACNPHB calculated with every electrostatic model are probably due to inadequate forcefield parameters for the iodine atom. The similar molecule ACMDOB (with chlorine instead of iodine but

932

Crystal Growth & Design, Vol. 5, No. 3, 2005

Figure 3. Two sample molecules for which charges are given in Tables 4 and 5. Table 5. Comparison of Charges for JECYIZ C1 H1a H1b C2 H2a H2b C3 H3a H3b C4 H4 C5 H5 C6 H6a H6b C7 H7a H7b H7c C8 H8a H8b H8c C9 H9 C10 H10a H10b C11 H11a H11b H11c P12 O13 N14 N15 N16 N17 N18 F

ESP static

ESP Slater

-0.131 0.063 0.063 -0.085 0.052 0.057 -0.177 0.083 0.085 0.170 0.063 -0.048 0.139 -0.140 0.083 0.095 -0.469 0.169 0.157 0.203 -0.287 0.141 0.139 0.134 0.050 0.131 -0.098 0.105 0.104 -0.441 0.127 0.136 0.120 0.785 -0.558 -0.473 -0.303 -0.538 0.683 -0.389 78

-0.120 0.087 0.057 -0.086 0.045 0.075 -0.188 0.084 0.075 0.180 0.079 -0.041 0.130 -0.127 0.093 0.118 -0.458 0.164 0.176 0.189 -0.301 0.131 0.153 0.115 0.073 0.162 -0.106 0.099 0.083 -0.451 0.122 0.127 0.124 0.766 -0.579 -0.461 -0.298 -0.535 0.663 -0.419 58

similar size and force field types) is described well. As ACNPHB is a planar molecule another possible source of error is out-of-plane polarization which cannot be described with charges bound to lie in the molecular plane. This might also hold for AMSALA01 which is planar too. For IBPRAC, the high F-values are due to a rotation of the molecule with respect to its experimental orientation. Molecular conformation, the hydrogen bonding pattern, and even the lattice parameters are still correct. The poor results for BANGOM02 with QEq charges are due to a distortion of the torsion angles in the bridge connecting the two aromatic rings, which leads to a distortion of the unit cell. For static molecular QEq charges, this is caused by an intermolecular hydrogen bond that is not present in the experimental structure. A possible source of errors, which is difficult to assess, is the dependence of the force-field parameters on the

Brodersen et al.

electrostatic model. The Dreiding force field was developed with Gasteiger charges. This investigation as our former work7 shows that ESP-derived electrostatics (charges or multipoles) perform significantly better than Gasteiger charges, especially for the more complicated molecules (other than class 1). Since the electrostatic energy is the most influential part of the nonbonded interaction, we believe that it is justified to neglect the dependence of the van der Waals parameters on the electrostatic model. Energies of Polymorphs In nine cases, the test set comprises two (LABHAX, IBPRAC, BANGOM, DUPTOX, CHAPEP, XIMBZB, BEYZIO, and CAXMOD) or three (DIWWEL) polymorphs. As described above, the polarization model allows changes in partial charges due to conformational and/or crystal packing effects. Although we do not have experimental data on the relative stability of these polymorphs, it is interesting to study how different charge models affect predicted energy differences between different polymorphs. In Table 3, energy differences ∆E in kcal/mol per asymmetric unit of these polymorphs are reported as calculated for the electrostatic models “ESP-static” and “ESP-Slater”. To ensure common zero points in energy, the same ESP charges (i.e., reference charges for the fluctuating charge model) were used for the polymorphs. These ESP charges were derived by fitting to all conformations of each molecule found in the polymorphic crystal structures at once. In all but one case (CHAPEP) the energy differences are below 6 kcal/mol, which is a reasonable value as it is of the order of the values usually found in experiments. For the fluctuating charge model (“ESP-Slater”) this finding shows that the inclusion of the polarization energy leads to comparable energies even though the charge distributions of different polymorphs are not the same after the crystal optimization. This is a necessary requirement for polymorph prediction where the decision on the most stable form is usually based on the energy ranking. For CHAPEP, the large differences of 21.5 kcal/mol for “ESP-static” and 18.1 kcal/mol for “ESP-Slater” indicate a failure of the force field. Examination of the results showed distortions in the PO3OH group of the anion. This may be an intrinsic Dreiding problem as in the original Dreiding publication a molecule with a PO3OH group (CSD code AFCDYP) also led to poor results when charges (Gasteiger charges with dielectric constant equal to unity) were used. However, a particular problem with phosphates was not mentioned in the original Dreiding publication. Both methods, “ESP-static” and “ESP-Slater”, agree in the sign of ∆E in seven out of nine cases (not considering CHAPEP) but not in the magnitude. As we are not aware of the experimental order of stability of the investigated polymorphs, no statement can be made regarding the correctness of the calculated energy differences. Conclusions Whereas in our former investigation7 we showed that the performance of multipoles for describing crystal

Comparison of Static and Fluctuating Charge Models

structures was mainly influenced by the flexibility of the molecules, this appeared to play a minor role in the present investigation. Therefore, the comparison of the different methods was discussed in terms of the class the structures were assigned to and their chemical similarity. However, the number of substances belonging to each class is still small, which limits the validity of our conclusions. Comparison of the static charge models confirmed that the empirical methods (i.e., no charges, Gasteiger, and QEq), averaged over all samples, are inferior to QMderived ESP charges. From the empirical models, QEq charges calculated for the crystal were considerably better than the other two models. The comparison of the different fluctuating charge models yielded a less clear picture. On one hand, the model seemed to capture crystal field effects as it improved over molecular QEq charges. On the other hand, it did not improve significantly over static ESP charges. As molecular QEq charges are not well suited for calculating the electrostatic energy in crystals, it seems that any improvement, even if not very accurate, helps a lot when using these charges. For ESP charges, which already yield good results, the fluctuating charge model has to be improved to make a greater impact. With original parameters taken from the QEq approach, the current model was not able to describe polarization effects accurately enough. Optimization of the parameters, which effectively reduces the polarization, yielded a slight improvement over static charges. Still, significant further research toward a more detailed model of polarization is needed. Acknowledgment. This work was supported by the Marie Curie Industry Host Fellowship programme of the European Commission. References (1) Verwer, P.; Leusen, F. J. J. In Reviews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; Wiley-VCH: John Wiley and Sons, New York, 1998; Vol. 12. (2) Lewis, T. C.; Tocher, D. A.; Price, S. L. Cryst. Growth Des. 2004, 4, 979. (3) Stone, A. J.; Price, S. L. J. Phys. Chem. 1988, 92, 3325.

Crystal Growth & Design, Vol. 5, No. 3, 2005 933 (4) Coombes, D. S.; Price, S. L.; Willock, D. J.; Leslie, M. J. Phys. Chem. 1996, 100, 7352. (5) Stone, A. J. Chem. Phys. Lett. 1981, 83, 233. (6) Stone, A. J.; Alderton, M. Mol. Phys. 2002, 100, 221. (7) Brodersen, S.; Wilke, S.; Leusen, F. J. J.; Engel, G. Phys. Chem. Chem. Phys. 2003, 5, 4923. (8) Ren, P.; Ponder, J. W. J. Comput. Chem. 2002, 23, 1497. (9) Gasteiger, J.; Marsili, M. Tetrahedron 1980, 36, 3219. (10) Rappe´, A. K.; Goddard, W. A., III. J. Phys Chem. 1991, 95, 3358. (11) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141. (12) Rick, S. W.; Berne, B. J. J. Am. Chem. Soc. 1996, 118, 672. (13) Winn, P. J.; Ferenczy, G. G.; Reynolds, C. A. J. Comput. Chem. 1999, 20, 704. (14) Patel, S.; Brooks, C. L., 3rd. J. Comput. Chem. 2004, 25, 1. (15) Ogawa, T.; Kitao, O.; Kurita, N.; Sekino, H.; Tanaka, S. Chem-Bio Informatics J. 2003, 3, 78. (16) Applequist, J.; Carl, J. R.; Fung, J. K. J. Am. Chem. Soc. 1972, 94, 2952. (17) van Duijnen, P. T.; Swart, M. J. Phys. Chem. A 1998, 102, 2399. (18) Ma, B.; Lii, J.-H.; Allinger, N. L. J. Comput. Chem. 2000, 21, 813. (19) Kaminski, G. A.; Friesner, R. A.; Zhou, R. J. Comput. Chem. 2003, 24, 267. (20) Ferenczy, G. G.; Reynolds, C. A. J. Phys. Chem. A 2001, 105, 11470. (21) Kitao, O.; Ogawa, T. Mol. Phys. 2003, 101, 3. (22) Nishimoto, K.; Mataga, N. Z. Phys. Chem. 1957, 12, 335. (23) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in C; Cambridge University Press: 1988. (24) Banks, J. L.; Kaminski, G. A.; Zhou, R.; Mainz, D. T.; Berne, B. J.; Friesner, R. A. J. Chem. Phys. 1999, 110, 741. (25) Stern, H. A.; Kaminski, G. A.; Banks, J. L.; Zhou, R.; Berne, B. J.; Friesner, R. A. J. Phys. Chem. B 1999, 103, 4730. (26) Delley, B. J. Chem. Phys. 1990, 92, 508. (27) Delley, B. J. Chem. Phys. 2000, 113, 7756. (28) DMol3 and Materials Studio are products of Accelrys Inc., 9685 Scranton Road, San Diego, CA 92121-3752, USA. (29) Hinsen, K.; Roux, B. J. Comput. Chem. 1997, 18, 368. (30) Mayo, S. L.; Olafson, B. D.; Goddard, W. A., III. J. Phys. Chem. 1990, 94, 8897. (31) Fillippi, G.; Gavezzotti, A. Acta Crystallogr. B 1993, 49, 868. (32) Allen, F. H. Acta Crystallogr. B 2002, 58, 380. (33) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179.

CG049673+