International
Tables for Crystallography Volume B Reciprocal space Edited by U. Shmueli © International Union of Crystallography 2010 |
International Tables for Crystallography (2010). Vol. B, ch. 3.5, pp. 458-460
Section 3.5.1. Introduction^{a}Laboratory of Structural Biology, National Institute of Environmental Health Sciences, 111 T. W. Alexander Drive, Research Triangle Park, NC 27709, USA |
High-precision single-crystal X-ray structural analysis of small organic molecules, yielding the space group, the unit-cell parameters and the fractional coordinates of the atoms making up the molecule(s) in the asymmetric unit, has become a routine matter as long as crystals of sufficient quality can be obtained. The thermodynamic stability of the crystal, as described by the enthalpy of sublimation ΔH_{sub}, can also be determined experimentally (although not always to high precision). Theoretical models for calculating intermolecular interaction energies can be used to connect the crystal structure to the molar enthalpy of sublimation using the relationshipwhere the lattice or packing energy E_{lattice} is the total (molar) intermolecular interaction energy between all the molecules in the crystal, which are treated as rigid entities with zero-point energies of intra- and intermolecular vibrations neglected. Connection to experimentally accessible heats of sublimation at higher temperatures involves thermodynamic corrections. Methods for calculating thermodynamic quantities of solids are discussed in Gavezzotti (2002a) and (in more detail) in Frenkel & Smit (2002).
Thus, given a parameterized intermolecular potential-energy function, or if computationally affordable a first-principles approach such as density-functional theory (or preferably, when it becomes feasible for crystals, a good-quality post-Hartree–Fock potential-energy surface that describes dispersion interactions), one can sum the intermolecular energies to obtain the lattice energy as a function of the above parameters defining the crystal structure. If such an energy function is used together with a method for systematic search of the crystal-structure parameters, one could in principle predict the minimum-lattice-energy crystal structure for a rigid organic molecule. To extend this approach to flexible molecules one would need to minimize the sum of the intramolecular energy plus the lattice energy. If the experimental crystal structure corresponds to the thermodynamic minimum-energy structure (i.e. it is not a metastable state determined by crystal-growth kinetics), one could in principle predict the experimental crystal structure of an organic compound through this minimization protocol. Moreover, one could ideally predict the additional metastable forms of the crystal.
Prediction of the structure of crystals of an organic molecule from its molecular structure is a difficult problem that has been compared to the protein folding problem (Dunitz, 2003; Dunitz & Scheraga, 2004). Like the protein folding problem, a solution of the crystal prediction problem has significant practical ramifications. A compound is often polymorphic, that is it has more than one crystal structure, and it may be difficult to characterize the conditions under which a particular crystal structure is formed (Dunitz & Bernstein, 1995). Polymorphs may have very different physical properties. An obvious example is diamond versus graphite, but other commercially important examples include food additives, various solid forms of explosives and the bioavailability of various forms of a drug such as ritonavir (Chemburkar et al., 2000). A method for predicting the possible crystal structures of the compound, and ideally for predicting the dominant crystal structure given the experimental conditions, would thus be very valuable. Note that due to the subtleties of crystallization, the lowest-free-energy polymorph at given temperature and pressure may not be the likeliest to form. The kinetics of growth of microcrystals may largely determine which low-energy polymorph appears (Dunitz, 2003). However, it is generally agreed that accurate calculation of the relative free energy of polymorphs is a prerequisite for predicting crystal structures.
To assess progress towards solving this latter problem, a series of blind tests of crystal-structure prediction has been undertaken (Day, Motherwell, Ammon et al., 2005). The results of these tests have highlighted the need for continued improvements in sampling methods and intermolecular energy potentials. Since extensive sampling of the crystal-structure parameters is necessary [between 10^{4} and 10^{5} starting structures, each followed by parameter minimizations (Price & Price, 2005)], there is a trade-off in the computational cost versus accuracy of the intermolecular energy functions used. Calculating the work of transforming between polymorphs is yet more ambitious in terms of sampling. Consequently empirical force fields are likely to be needed for the near term at least.
In the remainder of the introduction we outline some of the approaches to empirical potentials used in the calculation of the lattice energy, and then, motivated by these developments, discuss techniques for efficient summation of the electrostatic and other slow-decaying interaction terms that occur in these potential functions.
Methodological developments in the intermolecular force fields used in crystal-structure prediction from early times to the present state of the art have been reviewed (Price & Price, 2005). Until recently, these force fields were made up of atom–atom interactions. The earliest involved only repulsion and dispersion, usually in the `exp-6' formwhere U^{MN} denotes the intermolecular potential energy between molecules M and N and r_{ij} is the distance between atoms and . Sometimes the exponential form in the above equation is replaced by a simpler power law, as in the Lennard–Jones potential. As was pointed out by Dunitz (2003), in comparison with more sophisticated force fields, this repulsion–dispersion form readily allows analysis of the significance of particular atom–atom interactions, since the interactions are short-ranged and thus can be localized. That is, the r^{−6} form of the attractive dispersion energy means that interaction energies are halved for every 12% increase in distance. In contrast, introduction of long-range Coulombic interactions not only entails subtleties in lattice summation (the subject of this contribution), but greatly complicates the assignment of `key' atom–atom interactions. Gavezzotti and Fillipini systematically explored the use of the exp-6 potential in fitting organic crystal structures with and without hydrogen-bond interactions (Gavezzotti & Fillipini, 1994). They were surprisingly successful in accounting for weak hydrogen bonding in this way, but selective use of point charges improved the directionality of the potential. Earlier, Williams derived exp-6 parameters for the atoms C, H, N, O, Cl, F and polar H for use in organic crystal structures, but found it necessary (Williams & Cox, 1984) to supplement these with selected point charges, both atomic and at off-atom sites. Price and co-workers have used a more sophisticated electrostatic model based on distributed multipole analysis (DMA) with multipoles up to hexadecapole order positioned at atom centres (Stone, 1996). A similar approach has been adopted by Mooji et al., with the addition of dipole polarizabilities to account for intermolecular induction energy (Mooji et al., 1999). The AMOEBA force field has similar capabilities (Ren & Ponder, 2003).
Recently, a systematic study (Day et al., 2004; Day, Motherwell & Jones, 2005) of the effect of the electrostatic model on crystal-structure prediction for rigid organic molecules was carried out. The exp-6 parameters from the W99 force field (Williams, 2001) were used in conjunction with point charges derived either from a bond-increment model or electrostatic potential (ESP) fit, or with a DMA model as used by Price and co-workers. Importantly, they found a systematic and significant improvement in prediction success going from the simple bond-increment model to the ESP charges to point multipoles. Similar results were found by Broderson et al. (2003). This is a positive indicator for the use of theory in attacking this difficult prediction problem. On the other hand, Broderson et al. found that more sophisticated electrostatic models did not help in predicting the packing of flexible molecules. As discussed by Price & Price (2005), crystal-packing energies of flexible molecules are more problematic for at least two reasons. First the intramolecular energy must be fully compatible with the intermolecular energy in order that energy deformation from the gas-phase minimum be on the same scale as the concomitant gain in packing energy. Secondly, the conformational dependence of the electrostatic model must be accurately accounted for, or any gains from a more sophisticated model are immediately negated. Existing force fields are unable to meet these challenges, which may explain the lower success rate for predicting crystal structures of flexible molecules. Price has developed a methodology for flexible molecules wherein intramolecular conformations are systematically explored. For each conformer studied, the gas-phase intramolecular energy and DMA electrostatic model are calculated using quantum-chemical methods. The optimal packing of that conformer is then obtained by treating it as a rigid molecule, and the resulting minimum lattice energy is added to its intramolecular energy to give its total energy. The total energy for that conformer is then compared with that of the other conformers to predict the optimal conformation and packing of the flexible molecule. This approach has proved successful in a number of cases, for example with the successful prediction of a novel second polymorph of 1-hydroxy-7-azabenzotriazole (Nowell et al., 2006).
The above electrostatic models are all derived from gas-phase quantum-chemical calculations. A different line of development has been the parameterization of electrostatic models from condensed-phase experimental or theoretical electron densities, which should better represent the polarization due to the crystal environment. This is accomplished using aspherical multipole refinement of experimental X-ray charge densities (Destro et al., 2000) or of theoretical charge densities from quantum-chemical calculations carried out under the space-group symmetry (Dovesi et al., 2005). Crystal structure analysis is usually carried out assuming the independent-atom model (IAM), where the individual atomic densities are modelled as those of the spherically symmetric isolated atoms, with deviations in the experimental density being fitted to anisotropic temperature factors modelling the thermal motion of the nuclei. However, the IAM breaks down for data of sufficiently high quality, and the residual or deformation density is then fitted by a linear combination of atom-centred density basis elements resembling the well known Slater atomic basis sets from quantum chemistry (i.e. solid harmonics times a decaying exponential in the distance r from the atom). Technical details are provided in Coppens (1997) and Chapter 1.2 of this volume. The theoretical density can be similarly approached as a sum of pro-molecule density (superposition of isolated atomic densities in the molecular geometry) plus deformation density modelled by a sum of atomic Slater-like density elements. Until recently, the theoretical densities were thought to be of lower quality than the experimental densities, but continued improvements in the level of theory that can be applied have resulted in much closer agreement (Coppens & Volkov, 2004). Coppens and co-workers as well as others (Spackman et al., 1988) have applied these modelled densities in studies of intermolecular electrostatic energies. The Slater-like functions are typically difficult to work with in this context, so multipole approximations are used except at close range, where numerical integration is used instead (Volkov et al., 2004).
A similar density-based approach using gas-phase theoretical density was developed previously (Gavezzotti, 2002b). Gavezzotti splits the density into numerical `pixels' (small cubes of volume) within which the density is considered constant, allowing in principle exact numerical calculation of intermolecular electrostatic interactions. He sums the pixel–pixel interactions discretely, using a hierarchical approach to increase efficiency. The electrostatic interactions are supplemented by induction and repulsion terms, also modelled by pixel-based interactions. These developments, in particular the large size of the positive and negative components of the electrostatic interactions treated in terms of actual electron density, have caused a reassessment of the idea that intermolecular interactions can be accurately decomposed even in principle into localized atom–atom interactions (Dunitz & Gavezzotti, 2005).
Cisneros et al. have also begun implementing an intermolecular potential based on the explicit modelling of gas-phase theoretical electron density, using the methodology of auxiliary basis set fitting (Cisneros et al., 2005). The density is expanded into atom and off-atom sites using Gaussian density functions similar to those used as atomic basis sets in quantum chemical codes. More recently, other components of the potential such as induction, charge transfer and non-isotropic site–site repulsion based on density overlap have been implemented (Piquemal et al., 2006) and efficient calculation of terms in periodic boundary conditions has also been implemented (Cisneros et al., 2006), generalizing previous developments for point multipoles (Sagui et al., 2004).
Summarizing the above developments, it appears that accurate evaluation of crystal lattice energies necessitates the use of complex electrostatic interactions beyond the spherical-atom approximation, and ultimately including an account of penetration effects through explicit modelling of electron density. The electrostatic lattice energies are differences between large positive and negative contributions that are long-ranged, necessitating an accurate treatment of their lattice summation. Dispersion interactions, modelled by higher-order inverse powers of distance, are short-ranged but their calculation can be greatly accelerated through lattice-summation techniques.
Ewald summation in ideal infinite crystals, for interactions that depend on inverse powers of site–site distance, is covered in Chapter 3.4 of this volume, where there is also a numerical demonstration of the advantages of the Ewald approach over direct numerical summation. In this contribution we focus on the lattice-summation problem for large but finite crystals, highlighting the role of summation order and leading to a derivation of the surface term and polarization response. We also consider extensions to Coulombic interactions between Gaussian density elements, both spherical and Hermite (multipolar), and methods to further accelerate the lattice summation to (almost) linear scaling.
These developments are organized into three sections. The first section treats lattice summation of point charges. We begin with definitions and an explanation of the origin of the conditional convergence problem for these sums, anticipating the nature of the results to follow. Following that we provide for the benefit of less mathematically oriented readers the standard physics-based derivation of the Ewald sum for point charges, minus any accounting for the order of convergence or resulting surface term. This is followed by a more careful derivation focused on the order of summation, including results for higher-order inverse powers as well as the Coulomb interactions, finishing with a discussion of the surface term, the polarization response to an external dielectric and the pressure calculation. Although pitched at a higher mathematical level, we have tried to provide many details in the derivations unavailable elsewhere for the benefit of the reader. The alert reader will note that the key identity, equation (3.5.2.19), that separates 1/r^{n} into short-ranged and long-ranged components, is the same as that used in Chapter 3.4 of this volume, highlighting the essential similarity of our approach. In the second section we generalize these point-charge results to lattice sums of interacting spherical Gaussian and Hermite Gaussian-based densities. The third and final section provides a discussion of two methods that greatly increase the numerical efficiency of these lattice sums in the case of large unit cells by utilizing the fast Fourier transform. Again we have tried to provide mathematical details not available elsewhere.
References
Broderson, S., Wilke, S., Leusen, F. J. L. & Engel, G. (2003). A study of different approaches to the electrostatic interaction in force field methods for organic crystals. Phys. Chem. Chem. Phys. 5, 4923–4931.Chemburkar, S. R., Bauer, J., Deming, K., Spiwek, H., Patel, K., Morris, J., Henry, R., Spanton, S., Dziki, W., Porter, W., Quick, J., Bauer, P., Donabauer, J., Narayanan, B. A., Soldani, M., Riley, D. & McFarland, K. (2000). Dealing with the impact of ritonavir polymorphs on the late stages of bulk drug process development. Org. Process Res. Dev. 4, 413–417.
Cisneros, G. A., Piquemal, J. P. & Darden, T. A. (2005). Intermolecular electrostatics using density fitting. J. Chem. Phys. 123, 044109.
Cisneros, G. A., Piquemal, J. P. & Darden, T. A. (2006). Generalization of the Gaussian electrostatic model: Extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods. J. Chem. Phys. 125, 184101.
Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. New York: Oxford University Press.
Coppens, P. & Volkov, A. (2004). The interplay between experiment and theory in charge-density analysis. Acta Cryst. A60, 357–364.
Day, G. M., Chisholm, J., Shan, N., Motherwell, W. D. S. & Jones, W. (2004). An assessment of lattice energy minimization for the prediction of molecular organic crystal structures. Cryst. Growth Des. 4, 1327–1340.
Day, G. M., Motherwell, W. D. S., Ammon, H. L., Boerrigter, S. X. M., Della Valle, R. G., Venuti, E., Dzyabchenko, A., Dunitz, J. D., Schweizer, B., van Eijck, B. P., Erk, P., Facelli, J. C., Bazterra, V. E., Ferrarro, M. B., Hofmann, D. W. M., Leusen, F. J. J., Liang, C., Pantelides, C. C., Karamertzanis, P. G., Price, S. L., Lewis, T. C., Nowell, H., Torrisi, A., Scheraga, H. A., Arnautova, Y. A., Schmidt, M. U. & Verwer, P. (2005). A third blind test of crystal structure prediction. Acta Cryst. B61, 511–527.
Day, G. M., Motherwell, W. D. S. & Jones, W. (2005). Beyond the isotropic atom model in crystal structure prediction of rigid molecules: Atomic multipoles versus point charges. Cryst. Growth Des. 5, 1023–1033.
Destro, R., Roversi, P., Barzaghi, M. & Marsh, R. E. (2000). Experimental charge density of α-glycine at 23 K. J. Phys. Chem. A, 104, 1047–1054.
Dovesi, R., Orlando, R., Civalleri, B., Roetti, C., Saunders, V. R. & Zicovich-Wilson, C. M. (2005). Crystal: a computational tool for the ab initio study of the electronic properties of crystals. Z. Kristallogr. 220, 571–573.
Dunitz, J. D. (2003). Are crystal structures predictable? Chem. Commun. 5, 545–548.
Dunitz, J. D. & Bernstein, J. (1995). Disappearing polymorphs. Acc. Chem. Res. 28, 193–200.
Dunitz, J. D. & Gavezzotti, A. (2005). Molecular recognition in organic crystals: Directed intermolecular bonds or nonlocalized bonding? Angew. Chem. 44, 1766–1787.
Dunitz, J. D. & Scheraga, H. A. (2004). Exercises in prognistication: Crystal structures and protein folding. Proc. Natl Acad. Sci. USA, 101, 14309–14311.
Frenkel, D. & Smit, B. (2002). Understanding Molecular Simulation. San Diego: Academic Press.
Gavezzotti, A. (2002a). Structure and intermolecular potentials in molecular crystals. Modelling Simul. Mater. Sci. Eng. 10, R1–R29.
Gavezzotti, A. (2002b). Calculation of intermolecular interaction energies by direct numerical interaction energies I. Electrostatic and polarization energies in molecular crystals. J. Phys. Chem. B, 106, 4145–4154.
Gavezzotti, A. & Fillipini, G. (1994). The geometry of the intermolecular x—Hy (x, y = N, O) hydrogen bond and the calibration of empirical hydrogen bond potentials. J. Phys. Chem. 98, 4831–4837.
Mooji, W. T. M., van Duijneveldt, F. B., van Duijneveldt-van de Rijdt, J. G. C. M. & van Eijck, B. P. (1999). Transferable ab initio intermolecular potentials. 1. Derivation from methanol dimer and trimer calculations. J. Phys. Chem. A, 103, 9872–9882.
Nowell, H., Frampton, C. S., Waite, J. & Price, S. L. (2006). Blind crystal structure prediction of a novel second polymorph of 1-hydroxy-7-azabenzotriazole. Acta Cryst. B62, 642–650.
Piquemal, J. P., Cisneros, G. A. & Darden, T. A. (2006). Towards a force field based on density fitting. J. Chem. Phys. 124, 104101.
Price, S. L. & Price, L. S. (2005). Modelling intermolecular forces for organic crystal structure prediction. Struct. Bond. 115, 81–123.
Ren, P. Y. & Ponder, J. W. (2003). Polarizable atomic multipole water model for molecular mechanics simulation. J. Phys. Chem. B, 107, 5933–5947.
Sagui, C., Pedersen, L. G. & Darden, T. A. (2004). Towards an accurate representation of electrostatics in classical force fields: Efficient implementation of multipolar interactions in biomolecular simulations. J. Chem. Phys. 120, 73–87.
Spackman, M. A., Weber, H. P. & Craven, B. M. (1988). Energies of molecular-interactions from Bragg diffraction data. J. Am. Chem. Soc. 110, 775–782.
Stone, A. J. (1996). The Theory of Intermolecular Forces. New York: Oxford University Press.
Volkov, A., Koritsansky, T. & Coppens, P. (2004). Combination of the exact potential and multipole methods (ep/mm) for evaluation of intermolecular interaction energies with pseudoatom implementation of electron densities. Chem. Phys. Lett. 391, 170–175.
Williams, D. E. (2001). Improved intermolecular force field for molecules containing H, C, N and O atoms, with application to nucleoside and peptide crystals. J. Comput. Chem. 22, 1154–1166.
Williams, D. E. & Cox, S. R. (1984). Nonbonded potentials for azahydrocarbons: the importance of the Coulombic interactions. Acta Cryst. B40, 404–417.