research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Journal logoBIOLOGICAL
CRYSTALLOGRAPHY
ISSN: 1399-0047

Polarizable atomic multipole X-ray refinement: application to peptide crystals

CROSSMARK_Color_square_no_text.svg

aDepartment of Chemistry, Stanford, CA 94305, USA, bDepartment of Molecular and Cellular Physiology, Stanford, CA 94305, USA, and cHoward Hughes Medical Institute, USA
*Correspondence e-mail: pande@stanford.edu, brunger@stanford.edu

(Received 24 April 2009; accepted 13 June 2009; online 14 August 2009)

Recent advances in computational chemistry have produced force fields based on a polarizable atomic multipole description of biomolecular electrostatics. In this work, the Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) force field is applied to restrained refinement of molecular models against X-ray diffraction data from peptide crystals. A new formalism is also developed to compute anisotropic and aspherical structure factors using fast Fourier transformation (FFT) of Cartesian Gaussian multipoles. Relative to direct summation, the FFT approach can give a speedup of more than an order of magnitude for aspherical refinement of ultrahigh-resolution data sets. Use of a sublattice formalism makes the method highly parallelizable. Application of the Cartesian Gaussian multipole scattering model to a series of four peptide crystals using multipole coefficients from the AMOEBA force field demonstrates that AMOEBA systematically underestimates electron density at bond centers. For the trigonal and tetrahedral bonding geometries common in organic chemistry, an atomic multipole expansion through hexadecapole order is required to explain bond electron density. Alternatively, the addition of inter­atomic scattering (IAS) sites to the AMOEBA-based density captured bonding effects with fewer parameters. For a series of four peptide crystals, the AMOEBA–IAS model lowered Rfree by 20–40% relative to the original spherically symmetric scattering model.

1. Introduction

The number of X-ray crystal structures in the Protein Data Bank (PDB) with a resolution of higher than 1.0 Å continues to increase rapidly (Berman et al., 2000[Berman, H. M., Westbrook, J., Feng, Z., Gilliland, G., Bhat, T. N., Weissig, H., Shindyalov, I. N. & Bourne, P. E. (2000). Nucleic Acids Res. 28, 235-242.]). In late 2002, there were already over 100 structures available at subatomic resolution (Afonine & Urzhumtsev, 2004[Afonine, P. V. & Urzhumtsev, A. (2004). Acta Cryst. A60, 19-32.]), while as of early 2009 the number had more than tripled to well over 300. Examples include the proteins lysozyme at 0.65 Å (Wang et al., 2007[Wang, J., Dauter, M., Alkire, R., Joachimiak, A. & Dauter, Z. (2007). Acta Cryst. D63, 1254-1268.]), aldose reductase at 0.66 Å (Howard et al., 2004[Howard, E. I., Sanishvili, R., Cachau, R. E., Mitschler, A., Chevrier, B., Barth, P., Lamour, V., Van Zandt, M., Sibley, E., Bon, C., Moras, D., Schneider, T. R., Joachimiak, A. & Podjarny, A. (2004). Proteins, 55, 792-804.]) and serine protease at 0.78 Å (Kuhn et al., 1998[Kuhn, P., Knapp, M., Soltis, S. M., Ganshaw, G., Thoene, M. & Bott, R. (1998). Biochemistry, 37, 13446-13452.]), as well as nucleic acid structures such as B-DNA at 0.74 Å (Kielkopf et al., 2000[Kielkopf, C. L., Ding, S., Kuhn, P. & Rees, D. C. (2000). J. Mol. Biol. 296, 787-801.]), Z-DNA at 0.60 Å (Tereshko et al., 2001[Tereshko, V., Wilds, C. J., Minasov, G., Prakash, T. P., Maier, M. A., Howard, A., Wawrzak, Z., Manoharan, M. & Egli, M. (2001). Nucleic Acids Res. 29, 1208-1215.]) and an RNA tetraplex at 0.61 Å (Deng et al., 2001[Deng, J. P., Xiong, Y. & Sundaralingam, M. (2001). Proc. Natl Acad. Sci. USA, 98, 13665-13670.]). Crystals that diffract to high resolution are ideal for studying valence-electron distributions (Jelsch et al., 2000[Jelsch, C., Teeter, M. M., Lamzin, V., Pichon-Pesme, V., Blessing, R. H. & Lecomte, C. (2000). Proc. Natl Acad. Sci. USA, 97, 3171-3176.]; Muzet et al., 2003[Muzet, N., Guillot, B., Jelsch, C., Howard, E. & Lecomte, C. (2003). Proc. Natl Acad. Sci. USA, 100, 8742-8747.]; Zarychta et al., 2007[Zarychta, B., Pichon-Pesme, V., Guillot, B., Lecomte, C. & Jelsch, C. (2007). Acta Cryst. A63, 108-125.]; Volkov et al., 2007[Volkov, A., Messerschmidt, M. & Coppens, P. (2007). Acta Cryst. D63, 160-170.]; Coppens & Volkov, 2004[Coppens, P. & Volkov, A. (2004). Acta Cryst. A60, 357-364.]) that dictate the electrostatic properties of macromolecules. Electrostatics, in turn, is one of the driving forces in protein and nucleic acid folding, which should be understood in detail in order to predict biomolecular thermodynamics and kinetics (Snow et al., 2002[Snow, C. D., Nguyen, N., Pande, V. S. & Gruebele, M. (2002). Nature (London), 420, 102-106.], 2005[Snow, C. D., Sorin, E. J., Rhee, Y. M. & Pande, V. S. (2005). Annu. Rev. Biophys. Biomol. Struct. 34, 43-69.]; Sorin & Pande, 2005[Sorin, E. J. & Pande, V. S. (2005). J. Comput. Chem. 26, 682-690.]; Pande et al., 2003[Pande, V. S., Baker, I., Chapman, J., Elmer, S. P., Khaliq, S., Larson, S. M., Rhee, Y. M., Shirts, M. R., Snow, C. D., Sorin, E. J. & Zagrovic, B. (2003). Biopolymers, 68, 91-109.]). In this work, we contribute an improved theory and algorithm for computing the anisotropic and aspherical valence-electron density of molecules for X-ray crystal structure refinement.

Calculation of structure factors is generally based on scattering factors defined by the isolated-atom model (IAM), which assumes that the electron density around each atom is spherically symmetric. However, subatomic resolution diffraction data capture aspherical features of the electron density that result from bonding and the local chemical environment. The difference between the IAM and the true electron density is defined as the deformation density. For example, aspherical electron-density models of diamond, silicon and germanium developed by DeMarco and Weiss and later by Dawson explained the peaks of deformation density at bond midpoints observed in the experimental data (Dawson, 1967a[Dawson, B. (1967a). Proc. R. Soc. London Ser. A, 298, 264-288.],b[Dawson, B. (1967b). Proc. R. Soc. London Ser. A, 298, 379-394.],c[Dawson, B. (1967c). Proc. R. Soc. London Ser. A, 298, 395-401.]; DeMarco & Weiss, 1965[DeMarco, J. J. & Weiss, R. J. (1965). Phys. Rev. 137, A1869.]). In these works, the IAM was augmented by atom-centered spherical harmonic expansions, whose physical consequence was to redistribute electron density from nonbonding lobes into the tetragonal arrangement of bond centers.

A variety of radial functions have been used in combination with atom-centered spherical harmonic expansions. Modified Gaussians were promoted by Dawson (1967a[Dawson, B. (1967a). Proc. R. Soc. London Ser. A, 298, 264-288.]), a set of harmonic oscillator wavefunctions by Kurki-Suonio (1968[Kurki-Suonio, K. (1968). Acta Cryst. A24, 379-390.]) and more recently a formalism based on Slater-type orbitals (STO) was described by Stewart and coworkers (Epstein et al., 1977[Epstein, J., Bentley, J. & Stewart, R. F. (1977). J. Chem. Phys. 66, 5564-5567.]; Cromer et al., 1976[Cromer, D. T., Larson, A. C. & Stewart, R. F. (1976). J. Chem. Phys. 65, 336-349.]; Stewart, 1979[Stewart, R. F. (1979). Chem. Phys. Lett. 65, 335-342.], 1977[Stewart, R. F. (1977). Chem. Phys. Lett. 49, 281-284.]) and by Hansen & Coppens (1978[Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909-921.]), which represents the current standard (Jelsch et al., 2005[Jelsch, C., Guillot, B., Lagoutte, A. & Lecomte, C. (2005). J. Appl. Cryst. 38, 38-54.]; Zarychta et al., 2007[Zarychta, B., Pichon-Pesme, V., Guillot, B., Lecomte, C. & Jelsch, C. (2007). Acta Cryst. A63, 108-125.]; Volkov et al., 2007[Volkov, A., Messerschmidt, M. & Coppens, P. (2007). Acta Cryst. D63, 160-170.]; Coppens, 2005[Coppens, P. (2005). Angew. Chem. Int. Ed. 44, 6810-6811.]). However, spherical harmonics are not the only basis set available to describe the angular dependence of the deformation density.

We first present a formulation of anisotropic and aspherical atomic densities based on Cartesian Gaussian multipoles, which leads to much simpler formulae for the calculation of structure factors via direct summation in reciprocal space than the STO-based theory of Hansen & Coppens (1978[Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909-921.]). We also demonstrate that Cartesian Gaussian multipoles allow the computation of structure factors via fast Fourier transformation (FFT) of the real-space electron density (Cooley & Tukey, 1965[Cooley, J. W. & Tukey, J. W. (1965). Math. Comput. 19, 297-301.]). The latter approach, originally proposed by Ten Eyck (1973[Ten Eyck, L. F. (1973). Acta Cryst. A29, 183-191.], 1977[Ten Eyck, L. F. (1977). Acta Cryst. A33, 486-492.]), is the basis of the efficient macromolecular refinement algorithms (Brünger, 1989[Brünger, A. T. (1989). Acta Cryst. A45, 42-50.]; Afonine & Urzhum­tsev, 2004[Afonine, P. V. & Urzhumtsev, A. (2004). Acta Cryst. A60, 19-32.]; Afonine et al., 2007[Afonine, P. V., Grosse-Kunstleve, R. W., Adams, P. D., Lunin, V. Y. & Urzhumtsev, A. (2007). Acta Cryst. D63, 1194-1197.]; Agarwal, 1978[Agarwal, R. C. (1978). Acta Cryst. A34, 791-809.]) implemented in programs such as CNS (Brünger et al., 1998[Brünger, A. T., Adams, P. D., Clore, G. M., DeLano, W. L., Gros, P., Grosse-Kunstleve, R. W., Jiang, J.-S., Kuszewski, J., Nilges, M., Pannu, N. S., Read, R. J., Rice, L. M., Simonson, T. & Warren, G. L. (1998). Acta Cryst. D54, 905-921.]; Brunger, 2007[Brunger, A. T. (2007). Nature Protoc. 2, 2728-2733.]) and PHENIX (Adams et al., 2002[Adams, P. D., Grosse-Kunstleve, R. W., Hung, L.-W., Ioerger, T. R., McCoy, A. J., Moriarty, N. W., Read, R. J., Sacchettini, J. C., Sauter, N. K. & Terwilliger, T. C. (2002). Acta Cryst. D58, 1948-1954.]). The sublattice method implemented in CNS lends itself to efficient parallelization (Brünger, 1989[Brünger, A. T. (1989). Acta Cryst. A45, 42-50.]).

Boys originally proposed Cartesian Gaussian functions as basis functions to solve the many-electron Schrödinger equation (Boys, 1950[Boys, S. F. (1950). Proc. R. Soc. London Ser. A, 200, 542-554.]). The advantage of Gaussians over STOs in this context is that two-electron integrals have analytic forms, which has led to the adoption of Gaussian basis sets for many ab initio calculations (Hehre et al., 1969[Hehre, W. J., Stewart, R. F. & Pople, J. A. (1969). J. Chem. Phys. 51, 2657-2664.], 1970[Hehre, W. J., Ditchfie, R., Stewart, R. F. & Pople, J. A. (1970). J. Chem. Phys. 52, 2769-2773.]). We note that the equivalence of spherical harmonics and Cartesian tensors is well known, with key relationships having been presented by Stone (1996[Stone, A. J. (1996). The Theory of Intermolecular Forces. Oxford: Clarendon Press.]) and Applequist (1989[Applequist, J. (1989). J. Phys. A, 22, 4303-4330.], 2002[Applequist, J. (2002). Theor. Chem. Acc. 107, 103-115.]).

We apply Cartesian Gaussian multipoles to restrained crystallographic refinements based on the Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) force-field electrostatic model (Ponder & Case, 2003[Ponder, J. W. & Case, D. A. (2003). Adv. Protein Chem. 66, 27-85.]; Ren & Ponder, 2002[Ren, P. Y. & Ponder, J. W. (2002). J. Comput. Chem. 23, 1497-1506.], 2003[Ren, P. Y. & Ponder, J. W. (2003). J. Phys. Chem. B, 107, 5933-5947.], 2004[Ren, P. Y. & Ponder, J. W. (2004). J. Phys. Chem. B, 108, 13427-13437.]; Schnieders et al., 2007[Schnieders, M. J., Baker, N. A., Ren, P. Y. & Ponder, J. W. (2007). J. Chem. Phys. 126, 124114.]; Schnieders & Ponder, 2007[Schnieders, M. J. & Ponder, J. W. (2007). J. Chem. Theory Comput. 3, 2083-2097.]). The AMOEBA electrostatic model is based on the superposition of permanent atomic multipoles truncated at quadrupoles and induced dipoles. Permanent electrostatics represents the electron density of a group of atoms in the absence of interactions with the environment, which may include other parts of the molecule or solvent. Groups are chosen to be relatively rigid in order to avoid conformational variability in the permanent multipole moments. Conversely, the induced dipoles of AMOEBA represent polarization, the response of the electron density to the local electric field.

Force fields are widely used to restrain macromolecular refinement by contributing forces to local optimizations and molecular dynamics (Brünger et al., 1987[Brünger, A. T., Kuriyan, J. & Karplus, M. (1987). Science, 235, 458-460.]), with the latter used within simulated-annealing algorithms to promote global optimization (Brünger, 1988[Brünger, A. T. (1988). J. Mol. Biol. 203, 803-816.], 1991[Brünger, A. T. (1991). Annu. Rev. Phys. Chem. 42, 197-223.]; Brünger et al., 1989[Brünger, A. T., Karplus, M. & Petsko, G. A. (1989). Acta Cryst. A45, 50-61.], 1990[Brünger, A. T., Krukowski, A. & Erickson, J. W. (1990). Acta Cryst. A46, 585-593.], 1997[Brünger, A. T., Adams, P. D. & Rice, L. M. (1997). Structure, 5, 325-336.]; Kuriyan et al., 1989[Kuriyan, J., Brünger, A. T., Karplus, M. & Hendrickson, W. A. (1989). Acta Cryst. A45, 396-409.]; Adams et al., 1997[Adams, P. D., Pannu, N. S., Read, R. J. & Brünger, A. T. (1997). Proc. Natl Acad. Sci. USA, 94, 5018-5023.]; Brünger & Rice, 1997[Brünger, A. T. & Rice, L. M. (1997). Methods Enzymol. 277, 243-269.]). Up to now, force fields in crystallography have been largely limited to the geometric and repulsive terms and have had no influence on the atomic scattering factors. Therefore, refinement using a scattering model based on AMOEBA electrostatics is novel and lends insight into the progress being made in the development of precise, transferable force fields. Another limitation of the use of force fields for restraining X-­ray refinement has been the lack of proper treatment of long-range electrostatic interactions, which is overcome in this work via use of particle-mesh Ewald summation (PME; Darden et al., 1993[Darden, T., York, D. & Pedersen, L. (1993). J. Chem. Phys. 98, 10089-10092.]; Essmann et al., 1995[Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H. & Pedersen, L. G. (1995). J. Chem. Phys. 103, 8577-8593.]; Sagui et al., 2004[Sagui, C., Pedersen, L. G. & Darden, T. A. (2004). J. Chem. Phys. 120, 73-87.]).

In addition to AMOEBA, polarizable force fields are being studied by a number of other groups. Maple and coworkers have pursued a model similar to AMOEBA, but with the permanent moments truncated at dipole order, which has shown promising results for protein–ligand complexes (Friesner et al., 2005[Friesner, R. A., Baldwin, R. L. & Baker, D. (2005). Adv. Protein Chem. 72, 79-104.]; Maple et al., 2005[Maple, J. R., Cao, Y. X., Damm, W. G., Halgren, T. A., Kaminski, G. A., Zhang, L. Y. & Friesner, R. A. (2005). J. Chem. Theory Comput. 1, 694-715.]). As an alternative to induced dipoles, Patel and Brooks employed a fluctuating-charge model of polarization (Patel & Brooks, 2006[Patel, S. & Brooks, C. L. (2006). Mol. Simul. 32, 231-249.]), while Lamoureux and Roux have demonstrated success using classical Drude oscillators (Lamoureux et al., 2006[Lamoureux, G., Harder, E., Vorobyov, I. V., Roux, B. & MacKerell, A. D. (2006). Chem. Phys. Lett. 418, 245-249.]; Lamoureux & Roux, 2003[Lamoureux, G. & Roux, B. (2003). J. Chem. Phys. 119, 3025-3039.]). In addition to polarization, Gresh and coworkers have developed a methodology to include nonclassical effects such as electrostatic penetration and charge transfer (Gresh, 2006[Gresh, N. (2006). Curr. Pharm. Des. 12, 2121-2158.]; Gresh et al., 2007[Gresh, N., Cisneros, G. A., Darden, T. A. & Piquemal, J. P. (2007). J. Chem. Theory Comput. 3, 1960-1986.]; Piquemal et al., 2006[Piquemal, J. P., Cisneros, G. A., Reinhardt, P., Gresh, N. & Darden, T. A. (2006). J. Chem. Phys. 124, 104101.], 2007[Piquemal, J. P., Chelli, R., Procacci, P. & Gresh, N. (2007). J. Phys. Chem. A, 111, 8170-8176.]).

Although classical potentials can be validated against a range of experimental observables, for example small-molecule solvation energies (Shirts et al., 2003[Shirts, M. R., Pitera, J. W., Swope, W. C. & Pande, V. S. (2003). J. Chem. Phys. 119, 5740-5761.]; Shirts & Pande, 2005[Shirts, M. R. & Pande, V. S. (2005). J. Chem. Phys. 122, 134508.]), high-resolution diffraction data can pinpoint deficiencies in an electrostatics model with high precision. For example, we show that truncation of permanent atomic multipoles at quadrupole order limits the ability of the AMOEBA model to place charge density at bond midpoints. We use an efficient solution to this limitation by refining partial charges at bond centers as originally proposed by Afonine et al. (2007[Afonine, P. V., Grosse-Kunstleve, R. W., Adams, P. D., Lunin, V. Y. & Urzhumtsev, A. (2007). Acta Cryst. D63, 1194-1197.]).

2. Theory

2.1. Subgrid fast Fourier transform

The starting point for this work is the subgrid fast Fourier transform algorithm (SGFFT), which will be briefly summarized (Brünger, 1989[Brünger, A. T. (1989). Acta Cryst. A45, 42-50.]). In FFT-based methods, the electron density is computed over a lattice chosen to be fine enough to avoid aliasing effects at a given resolution. This computation can be made more efficient by an artificial increase in the atomic displacement parameters (ADPs) of all atoms. The optimum choice in CNS v.1.2 (Brunger, 2007[Brunger, A. T. (2007). Nature Protoc. 2, 2728-2733.]) for the ADP offset and grid size follows the work of Bricogne (2001[Bricogne, G. (2001). International Tables for Crystallography, Vol. B, edited by U. Schmueli, pp. 25-98. Dordrecht: Kluwer Academic Publishers.]). An important point is that the electron density is only computed within a cutoff radius around each atom. As the resolution increases, the cutoff is increased based on an empirical scheme to maintain agreement between direct-summation structure factors and derivatives and the SGFFT calculation (Brunger, 2007[Brunger, A. T. (2007). Nature Protoc. 2, 2728-2733.]).

Structure factors are computed by FFT of the electron density of an asymmetric unit of atoms (Agarwal, 1978[Agarwal, R. C. (1978). Acta Cryst. A34, 791-809.]). The SGFFT is based on factorizing this computation into smaller FFTs that are computed separately on sublattices, which allows efficient parallelization since these tasks are independent (Brünger, 1989[Brünger, A. T. (1989). Acta Cryst. A45, 42-50.]; Kay Diederichs, private communication). CNS v.1.21 has implemented this approach via an OpenMP environment (courtesy Kay Diederichs, University of Konstanz; available at https://cns-online.org ). Crystallographic symmetry is then applied to the structure factors, and the target function and its derivatives with respect to structure factors are evaluated. Symmetry operators are applied to the derivatives of the target function with respect to the structure factors followed by inverse Fourier transform. Using the chain rule, derivatives of the target function with respect to atomic parameters are then computed by multiplication and summation over the local neighborhood around each atom of the derivatives of the electron density with respect to atomic parameters.

Although the original SGFFT method was developed with an isolated-atom description of electron density and isotropic ADPs, it is generalizable to aspherical Cartesian Gaussian multipoles and anisotropic ADPs. All that is needed are formulae for the electron density and the derivatives of the electron density with respect to atomic parameters, which then can be inserted into equations (29) and (40) of Brünger (1989[Brünger, A. T. (1989). Acta Cryst. A45, 42-50.]). In the following sections, we develop these necessary formulae.

2.2. Isolated-atom Gaussian density

The key mathematical property of Gaussians with respect to efficient calculation of structure factors is that they are an eigenfunction of the Fourier transform (FT). In other words, a Gaussian in real space transforms to a Gaussian in reciprocal space and vice versa. Consider the canonical spherically symmetric Gaussian atomic scattering factor (Agarwal, 1978[Agarwal, R. C. (1978). Acta Cryst. A34, 791-809.]),

[f^{(n,\kappa)}({\bf r}) = \kappa^3 (4\pi)^{3/2}{\textstyle \sum\limits_{i = 1}^n} {{a_i} \over {b_i^{3/2} }} \exp\left (- {{ 4\pi^2 \kappa ^2|{\bf r}|^2 } \over {b_i}}\right), \eqno (1)]

where ai and bi are constant parameters fitted to ab initio calculations on isolated atoms (this work is based on a sum of six Gaussians; n = 6; Su & Coppens, 1998[Su, Z. & Coppens, P. (1998). Acta Cryst. A54, 357.]), κ is an expansion/contraction parameter used to adjust the width of the density and r is a position vector relative to the center of the atom. Its FT is given by

[\hat f^{(n,\kappa)} ({\bf s}) = {\textstyle \sum\limits_{i = 1}^n} a_i \exp \left (- {{b_i|{\bf s}|^2} \over {4\kappa^2}}\right), \eqno (2)]

where s is the reciprocal-lattice vector and we have used the FT definition given in Appendix A[link]. The reciprocal-lattice vector is s = htA−1 = (A−1)th, where h is a column vector with the Miller indices of a Bragg reflection and A is the fractional­ization matrix that transforms coordinates r with respect to a Cartesian basis to fractional coordinates rfrac as defined in a crystallographic basis set. The Debye–Waller factor (Waller, 1923[Waller, I. (1923). Z. Phys. A, 17, 398-408.]) is given by

[\hat t({\bf s}) = \exp (-2\pi^2{\bf s}^t{\bf Us}) \eqno (3)]

in reciprocal space, where each element of the symmetric positive-definite matrix U is defined via a Cartesian basis consistent with PDB ANISOU records (Trueblood et al., 1996[Trueblood, K. N., Bürgi, H.-B., Burzlaff, H., Dunitz, J. D., Gramaccioli, C. M., Schulz, H. H., Shmueli, U. & Abrahams, S. C. (1996). Acta Cryst. A52, 770-781.]; Grosse-Kunstleve & Adams, 2002[Grosse-Kunstleve, R. W. & Adams, P. D. (2002). J. Appl. Cryst. 35, 477-480.]). Multiplication of (3)[link] by the atomic form factor from (2)[link] gives the scattering factor

[\hat \rho^{(n,\kappa)}({\bf s}) = \hat f^{(n,\kappa)}({\bf s})\hat t({\bf s}) = \textstyle \sum\limits_{i = 1}^n a_i \exp(- 2\pi^2 {\bf s}^t {\bf U}_i {\bf s}) \eqno (4)]

based on Ui that are defined by

[{\bf U}_i = \left(\matrix{U_{11} & U_{12} & U_{13} \cr U_{21} & U_{22} & U_{23} \cr U_{31} & U_{32} & U_{33}} \right) + I_3 \left({{b_i} \over {8\pi^2\kappa ^2}} + U_{\rm add} \right), \eqno (5)]

where Uadd is the artificial isotropic increase or decrease in the ADP discussed above and I3 is a 3 × 3 identity matrix. Removal of Uadd analytically from each structure factor after the FT is straightforward. The only difference, therefore, between each Ui is the isolated-atom scattering parameter bi.

Application of the inverse FT to (4)[link] gives the real-space anisotropic electron density

[\rho^{(n,\kappa)}({\bf r}) = (2\pi)^{-3/2} \textstyle \sum\limits_{i = 1}^n a_i|{\bf U}_i|^{-1/2} \exp(-{1 \over 2}{\bf r}^t {\bf U}_i^{-1} {\bf r}), \eqno (6)]

where |Ui| is the determinant of matrix Ui and Ui−1 is its inverse. This expression can also be viewed as the convolution of the Gaussian form factor of (1)[link] with the inverse Fourier transform of the Debye–Waller factor of (3)[link]. Although the underlying isolated-atom scattering factor is spherically symmetric, convolution with anisotropic ADPs can lead to an angular dependence in ρ(n,κ)(r). Using the relationship that B = 8π2U, one can show that (6)[link] reduces to the isotropic density expression reported by Brünger in equation (16) of Brünger (1989[Brünger, A. T. (1989). Acta Cryst. A45, 42-50.]) if all diagonal elements of Ui are equal to Uiso + bi/8π2 + Uadd with zero off-diagonal components.

2.3. Polarizable atomic multipole electron density

For the derivation of an atomic multipole expansion from a collection of point charges we begin with the Taylor expansion of the electric potential V(r) at r arising from n partial point charges that represent the electron density of an atom,

[\eqalignno {V({\bf r}) & = {\textstyle \sum\limits_{i = 1}^n} {{c_i} \over {|{\bf r}-\Delta_i|}} \cr & = {\textstyle \sum\limits_{i = 1}^n} c_i \biggr [{1 \over r} + \Delta_{i,\alpha}\left({\partial \over {\partial \Delta_{i,\alpha}}} {1 \over {|{\bf r} - \Delta _i|}}\right)_{\Delta_{i} = 0} \cr &\ \quad +\ {1 \over 2}\Delta _{i,\alpha} \Delta _{i,\beta} \left({\partial \over {\partial \Delta_{i,\alpha}}}{\partial \over {\partial \Delta _{i,\beta}}}{1 \over {|{\bf r} - \Delta _i|}} \right)_{\Delta_i = 0} + \ldots \biggr] \cr & = {\textstyle \sum\limits_{i = 1}^n}c_i \left[{1 \over r} - \Delta_{i,\alpha} {\partial \over {\partial r_\alpha}}{1 \over r} + {1 \over 2}\Delta_{i,\alpha} \Delta_{i,\beta} {\partial \over {\partial r_\alpha }}{\partial \over {\partial r_\beta }}{1 \over r} - \ldots \right] \cr & = {\textstyle \sum\limits_{i = 1}^n} \left(1 - \Delta_{i,\alpha} \nabla_\alpha + {1 \over 2}\Delta_{i,\alpha} \Delta_{i,\beta} \nabla_\alpha \nabla _\beta - \ldots \right) {{c_i} \over r}, & (7)}]

where Δi is the position of partial charge ci, ∇α = ∂/∂rα is one component of the del operator, α ∈ {xyz} and the Greek subscripts {α, β} represent the use of the Einstein summation convention for summing over tensor elements (Stone, 1996[Stone, A. J. (1996). The Theory of Intermolecular Forces. Oxford: Clarendon Press.]). We omit the constant factor of 1/4π0 throughout for com­pactness. Let the monopole, dipole and traceless quadrupole moments be defined as

[\eqalignno {q & = \textstyle \sum\limits_{i = 1}^n c_i, \cr d_\alpha & = \textstyle \sum\limits_{i = 1}^n c_i \Delta _{i,\alpha}, \cr \Theta_{\alpha \beta} & = \textstyle{3 \over 2} \textstyle\sum\limits_{i = 1}^n c_i (\Delta _{i,\alpha} \Delta _{i,\beta} - {1 \over 3}\Delta _i^2 \delta _{\alpha \beta}), & (8)}]

where removal of the trace in the definition of the quadrupole moment is allowed because the potential satisfies the Laplace equation (i.e.2V = 0). Substitution of the relationships in (8)[link] into the final expression of (7)[link] gives the electric potential in terms of a Cartesian multipole expansion, which we truncate at quadrupole order

[V({\bf r}) = (q - d_\alpha \nabla_\alpha + {\textstyle {1 \over 3}}\Theta_{\alpha \beta} \nabla_\alpha \nabla_\beta ){1 \over r}. \eqno (9)]

We now replace the Coulomb potential of (9)[link] with the potential from the sum of Gaussians from (1)[link], which is given by

[\varphi^{(n,\kappa)} (r) = {\textstyle \sum\limits_{i = 1}^n} a_i {{ {\rm erf}(2\pi \kappa r/ b_i^{1/2}) } \over r}\eqno (10)]

and find

[\Phi({\bf r}) = (q - d_\alpha \nabla_\alpha + \textstyle {1 \over 3}\Theta _{\alpha \beta} \nabla_\alpha \nabla_\beta)\varphi^{(n,\kappa)}(r). \eqno (11)]

We now introduce unique superscripts on the charge, dipole and quadrupole Gaussian basis sets, denoted by {nqndnΘ} and {κq, κd, κΘ}, to allow them to differ in number and width.

[\Phi ({\bf r}) = q\varphi^{(n_q, \kappa _q)}(r) - d_\alpha \nabla_\alpha \varphi ^{(n_d, \kappa _d)}(r) + \textstyle {1 \over 3}\Theta_{\alpha \beta } \nabla_\alpha \nabla_\beta \varphi^{(n_\Theta, \kappa _\Theta)} (r). \eqno (12)]

The potential of the charge density of (12)[link] quickly approaches the Coulomb potential as r increases since the error function goes to unity such that at large r this potential satisfies the Laplace equation and the use of a traceless quadrupole tensor is still justified. Application of the Laplace operator to both sides of (12)[link] gives the negative of a continuous charge density based on Cartesian Gaussian multipoles,

[\rho({\bf r}) = - qf^{(n_q, \kappa _q)}({\bf r}) + d_\alpha \nabla_\alpha f^{(n_d, \kappa _d)} ({\bf r}) - \textstyle {1 \over 3}\Theta_{\alpha \beta} \nabla_\alpha \nabla_\beta f^{(n_\Theta, \kappa _\Theta)}({\bf r}). \eqno (13)]

In crystallography the convention is that electron density is positive, so we will keep the negative sign. Therefore, a negative partial charge equates to positive scattering density.

Inclusion of ADPs is described by convolution of (13)[link] with the real-space temperature factor,

[\rho_{\rm ADP} ({\bf r}) = \rho({\bf r}) * {\rm t}({\bf r}).\eqno (14)]

Based on the convolution differentiation rule

[[\nabla_\alpha \rho ({\bf r})] * t({\bf r}) \equiv \nabla_\alpha [\rho ({\bf r}) * t({\bf r})] \eqno (15)]

the solution to (14)[link] is given by substituting for f(r) in (13)[link] with the corresponding ρ(r) from (6)[link] to give

[\eqalignno {\rho _{\rm ADP} ({\bf r}) &= -q\rho^{(n_q, \kappa _q)}({\bf r}) + d_\alpha \nabla_\alpha \rho^{(n_d, \kappa _d)}({\bf r}) \cr &\ \quad -\ \textstyle {1 \over 3}\Theta_{\alpha \beta} \nabla_\alpha \nabla_\beta \rho^{(n_\Theta, \kappa _\Theta)}({\bf r}). & (16)}]

However, since q only represents partial atomic charges, the contributions from valence and core electrons need to be added. Additionally, the AMOEBA force field divides each atomic dipole moment into permanent (d) and induced (u) contributions to account for polarization. Therefore, we construct the total atomic electron density at a location r relative to the center of atom j by adding the contribution of core and valence electron density to (16)[link] and splitting the dipole into permanent and induced components to give

[\eqalignno {\rho_{{\rm ADP},j}({\bf r}) & = P_j^{(c)} \rho_j^{(6,1)}({\bf r}) + [P_j^{(v)} - q_j] \rho_j^{(6,\kappa _v)}({\bf r}) \cr &\ \quad +\ (d_{j,\alpha } + u_{j,\alpha })\nabla_\alpha \rho_j^{(1,\kappa _d)}({\bf r}) - \textstyle {1 \over 3}\Theta_{j,\alpha \beta} \nabla_\alpha \nabla_\beta \rho_j^{(1,\kappa _\Theta)}({\bf r}), \cr &&(17)}]

where Pj(c) is the integer number of core electrons (carbon has two) and Pj(v) is the integer number of valence electrons (carbon has four). The superscripts on the anisotropic Gaussian form factors ρj(n,κ)(r) have been made explicit for our model. We make the reasonable choice of using the isolated-atom scattering parameters for both core and valence electron densities. The width of the core electron density is frozen at the isolated-atom description (κ = 1) based on the observation that chemical bonding does not perturb it significantly (Hansen & Coppens, 1978[Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909-921.]). On the other hand, the width of the valence electron density expands or contracts relative to the isolated-atom model owing to a gain or reduction, respectively, of electron density from or to covalently bonded atoms. This effect is modeled by the width parameter of the valence density κv. In this work, the dipole and quadrupole densities are described by a single Gaussian (nd = nΘ = 1) based on a and b parameters set to unity. The widths of the dipole and quadrupole densities are controlled by the κd and κΘ parameters. In this work, the width parameters {κv, κd, κΘ} are optimized against the diffraction data for each AMOEBA multipole type. The multipole moments are fixed by the AMOEBA force field and are not refined against the data.

The partial derivatives through second order of the anisotropic and aspherical density defined in (6)[link], which are required for the real-space multipolar density given in (17)[link], are

[\eqalignno {\nabla_\alpha \rho^{(n,\kappa)}({\bf r}) & = -(2\pi)^{- 3/2} \textstyle \sum\limits_{i = 1}^n a_i |{\bf U}_i|^{- 1/2} \cr &\ \quad {\times} \exp(-\textstyle{1 \over 2}{\bf r}^t {\bf U}_i^{- 1} {\bf r}) ({\bf r}^t {\bf U}_i^{-1} {\bf u}_\alpha) \cr \nabla_\alpha \nabla_\beta \rho^{(n,\kappa)}({\bf r}) & = (2\pi)^{- 3/2}\textstyle \sum\limits_{i = 1}^n a_i|{\bf U}_i|^{ - 1/2} \cr &\ \quad {\times} \exp (- \textstyle {1 \over 2}{\bf r}^t {\bf U}_i^{-1} {\bf r})[({\bf r}^t{\bf U}_i^{-1}{\bf u}_\alpha)({\bf r}^t {\bf U}_i^{-1} {\bf u}_\beta) - U_{i,\alpha \beta }^{ - 1}], \cr &&(18)}]

where uα is a unit vector in the α direction with α ∈ {xyz}. In addition, the third-, fourth- and fifth-order terms of the expansion are presented as supplementary information along with a Mathematica notebook.1

To the best of our knowledge, (17)[link] is the first expression reported in the literature for a real-space form factor that is the convolution of an atomic multipolar electron density with anisotropic ADPs. This equation opens the door to exploring precise polarizable atomic multipole refinements in tandem with efficient computation of structure factors via FFT.

Given a molecular conformation, the AMOEBA permanent multipole moments for each atom in the global coordinate frame (q, d, Θ) are converted via rotation from a local frame. For example, as shown in Fig. 1,[link] the z axis of the local frame for the carbonyl O atom of the peptide bond is in the direction of the bond to the carbonyl C atom. Its positive x axis is located in the O=C—Cα plane in the direction of the Cα atom and the y axis is chosen to give a right-handed coordinate system (Ren & Ponder, 2002[Ren, P. Y. & Ponder, J. W. (2002). J. Comput. Chem. 23, 1497-1506.]). The induced dipole (u) on each atom is determined via a self-consistent field (SCF) calculation, where the field is a sum of contributions from the permanent atomic multipoles and induced dipoles. The AMOEBA polarization model is described in greater detail in work by Ren & Ponder (2002[Ren, P. Y. & Ponder, J. W. (2002). J. Comput. Chem. 23, 1497-1506.]).

[Figure 1]
Figure 1
The local multipole frame of the carbonyl O atom of the peptide backbone is shown. The positive z axis is along the C=O bond and the x axis is chosen in the O=C—Cα plane in the direction of the Cα atom. The y axis is directed into the page in order to achieve a right-handed coordinate system. Also shown are the nonzero multipole moments of the O atom and a qualitative representation of their shape. The dz Cartesian Gaussian dipole (in Debye units) places electron density along the C=O bond, while the trace of the Cartesian Gaussian quadrupole (in Buckingham units) positions electron density approximately at lone-pair positions.

2.4. Derivatives of the electron density

2.4.1. Atomic coordinates

As a simplification, the derivation up to this point has assumed that the atomic center was the origin of the coordinate system. However, for this section on the derivatives with respect to atomic coordinates we place atom j at rj in the global frame. In order to keep the derivation manageable, we split the total electron density into that produced by permanent charges ρperm and that of induced charges ρind,

[\rho_{\rm total} ({\bf r}) = \textstyle \sum\limits_{j = 1}^n \rho_{{\rm perm}, j} ({\bf r} - {\bf r}_j) + \rho _{{\rm ind},j} ({\bf r} - {\bf r}_j). \eqno (19)]

The derivative of the permanent multipole electron density of atom j with respect to the α coordinate of atom j is given by

[\eqalignno {{{\partial \rho_{{\rm perm}, j} ({\bf r} -{\bf r}_j)} \over {\partial r_{j,\alpha}}} &= P_j^{(c)}{{\partial \rho^{(6,1)}_j ({\bf r} - {\bf r}_j)} \over {\partial r_{j, \alpha }}} \cr + \big[P_j^{(v)} & - q_j\big] {{\partial \rho_j^{(6,\kappa_v)}({\bf r} - {\bf r}_j)} \over {\partial r_{j, \alpha} }} + d_{j,\beta} {{\partial [\nabla_\beta \rho_j^{(1,\kappa _d)} ({\bf r} - {\bf r}_j)]} \over {\partial r_{j, \alpha} }} \cr + {{\partial d_{j,\beta}} \over {\partial r_{j, \alpha} }} &\nabla_\beta \rho_j^{(1,\kappa _d)}({\bf r} - {\bf r}_j) - {1 \over 3} \Theta_{j,\beta \gamma} {{\partial [\nabla_\beta \nabla_\gamma \rho_j^{(1,\kappa _\Theta)}({\bf r} - {\bf r}_j)]} \over {\partial r_{j, \alpha} }} \cr &- {1 \over 3}{{\partial \Theta_{j,\beta \gamma}} \over {\partial r_{j, \alpha}}}\nabla_\beta \nabla_\gamma \rho_j^{(1,\kappa _\Theta)}({\bf r} - {\bf r}_j), & (20)}]

where the derivative of the dipole and quadrupole densities are each composed of two terms owing to the chain rule. As described above, the dipole and quadrupole moments of each atom are implicitly a function of its coordinates and the coordinates of a few of its bonded neighbors (atoms k) that define the local frame of the multipole. Therefore, the derivative of the permanent multipole electron density of atom j with respect to the α coordinate of atoms k must also be considered,

[\eqalignno {{{\partial \rho_{{\rm perm}, j}({\bf r}-{\bf r}_j)} \over {\partial r_{k,\alpha} }} &= {{\partial d_{j,\beta} } \over {\partial {r}_{k,\alpha} }} \nabla_\beta \rho_j^{(1,\kappa_d)}({\bf r} - {\bf r}_j) \cr &\ \quad -\ {1 \over 3} {{\partial \Theta_{j,\beta \gamma} } \over {\partial r_{k,\alpha} }}\nabla_\beta \nabla_\gamma \rho_j^{(1,\kappa _\Theta)}({\bf r} - {\bf r}_j), & (21)}]

where the derivatives of spherically symmetric terms are zero with respect to the coordinates of atom k because they have no dependence on the orientation of the local frame. Note that the partial derivative of an anisotropic and aspherical density tensor with respect to an atomic coordinate is the negative of the partial derivatives given in (18)[link], simply due to the negative sign on rj. The derivatives of the polarizable density with respect to atomic coordinates are very specific to the AMOEBA electrostatic model and are discussed in Appendix B[link]. However, we note that computing the derivatives of a polarizable density with respect to atomic coordinates is O(n2logn) using PME, which quickly becomes the most expensive part of the overall calculation.

2.4.2. ADPs

The derivative of the anisotropic electron density of atom j with respect to an anisotropic displacement parameter Uj,τυ is given by

[\eqalignno {{{\partial \rho _j ({\bf r})} \over {\partial U_{j,\tau \upsilon} }} & = P_j^{(c)} {{\partial \rho_j^{(6,1)}({\bf r})} \over {\partial U_{j,\tau \upsilon} }} + [P_j^{(v)} - q_j] {{\partial \rho_j^{(6,\kappa _v)}({\bf r})} \over {\partial U_{j,\tau \upsilon} }} \cr &\ \quad +\ (d_{j,\alpha} + u_{j,\alpha}){{\partial [\nabla_\alpha \rho _j^{(1,\kappa _d)}({\bf r})]} \over {\partial U_{j,\tau \upsilon } }} \cr &\ \quad -\ {1 \over 3}\Theta _{j,\alpha \beta } {{\partial [\nabla_\alpha \nabla_\beta \rho_j^{(1,\kappa _\Theta)}({\bf r})]} \over {\partial U_{j,\tau \upsilon }}}& (22)}]

and requires the partial derivatives of the Cartesian Gaussian tensors with respect to ADP components. Introducing a few relationships facilitates their presentation. Firstly, based on the equality

[{{\partial |{\bf U}|} \over {\partial U_{\tau \upsilon} }} = |{\bf U}|U_{\tau \upsilon}^{-1} (2 - \delta_{\tau \upsilon }) \eqno (23)]

we have

[{{\partial |{\bf U}|^{-1/2} } \over {\partial U_{\tau \upsilon}}} = \textstyle -{1 \over 2}| {\bf U}|^{-1/2} U_{\tau \upsilon}^{-1} (2 - \delta_{\tau \upsilon }), \eqno (24)]

where the Kronecker delta δτυ is unity for diagonal elements of U and zero otherwise. Differentiating an identity from matrix algebra U−1U = I gives the following relationship

[{{\partial {\bf U}^{-1}} \over {\partial U_{\tau \upsilon}}} = -{\bf U}^{-1} {{\partial {\bf U}} \over {\partial U_{\tau \upsilon}}}{\bf U}^{-1}, \eqno (25)]

which makes it possible to differentiate U instead of its inverse. This is preferred since only one or two elements of ∂U/∂Uτυ are equal to unity and the rest are zero. Specifically, a single element is equal to unity if τ equals υ, while two elements are equal to unity otherwise, since Uτυ and Uυτ represent the same variable in this case. For convenience, we define a 3 × 3 matrix J(τυ),

[{\bf J}^{(\tau \upsilon)} = -{\bf U}^{-1} {{\partial {\bf U}} \over {\partial U_{\tau \upsilon}}} {\bf U}^{-1}, \eqno (26)]

and based on the chain rule we have

[{\partial \over {\partial U_{\tau \upsilon} }}\exp(- {\textstyle {1 \over 2}}{\bf r}^t {\bf U}^{-1} {\bf r}) = - {\textstyle {1 \over 2}} {\bf r}^t {\bf J}^{(\tau \upsilon)} {\bf r}\exp(-{\textstyle {1 \over 2}}{\bf r}^t {\bf U}^{-1} {\bf r}). \eqno (27)]

Differentiating (6)[link] with respect to Uτυ and using (24)[link], (27)[link] and the product rule gives

[\eqalignno {{{\partial \rho^{(n,\kappa)}({\bf r})} \over {\partial U_{\tau \upsilon}}} =\ &(2\pi)^{-3/2} \textstyle \sum\limits_{i = 1}^n a_i | {\bf U}_i|^{-1/2} \cr &\times \exp(- {\textstyle {1 \over 2}}{\bf r}^t{\bf U}_i^{-1}{\bf r}){\textstyle {1 \over 2}}[-{\bf r}^t {\bf J}_i^{(\tau \upsilon)}{\bf r} - U_{i,\tau \upsilon}^{-1} (2 - \delta _{\tau \upsilon})] \cr {{\partial [\nabla_\alpha \rho^{(n,\kappa)}({\bf r})]} \over {\partial U_{\tau \upsilon} }} & = -(2\pi)^{- 3/2} \textstyle \sum\limits_{i = 1}^n a_i|{\bf U}_i|^{-1/2}\exp(- {1 \over 2}{\bf r}^t {\bf U}_i^{-1} {\bf r}) \cr \times\ \{ \textstyle {1 \over 2}[- &{\bf r}^t {\bf J}_i^{(\tau \upsilon)} {\bf r} - U_{i,\tau \upsilon }^{-1}(2 - \delta _{\tau \upsilon})]({\bf r}^t{\bf U}_i^{-1} {\bf u}_\alpha) + ({\bf r}^t {\bf J}_i^{(\tau \upsilon)} {\bf u}_\alpha)\} \cr {{\partial[\nabla_\alpha \nabla_\beta \rho^{(n,\kappa)}({\bf r})]} \over {\partial U_{\tau \upsilon} }} & = (2\pi)^{-3/2}\textstyle \sum\limits_{i = 1}^n a_i|{\bf U}_i|^{-1/2} \exp(-{1 \over 2}{\bf r}^t {\bf U}_i^{-1} {\bf r}) \cr \times\ \{ \textstyle {1 \over 2}[-{\bf r}^t {\bf J}_i^{(\tau \upsilon)} &{\bf r} - U_{i,\tau \upsilon}^{-1}(2 - \delta_{\tau \upsilon})][({\bf r}^t{\bf U}_i^{-1}{\bf u}_\alpha)({\bf r}^t {\bf U}_i^{-1}{\bf u}_\beta) - U_{i,\alpha \beta}^{-1}] \cr +\ [{\bf r}^t{\bf J}_i^{(\tau \upsilon)}{\bf u}_\alpha]&({\bf r}^t {\bf U}_i^{-1} {\bf u}_\beta) + ({\bf r}^t{\bf U}_i^{-1} {\bf u}_\alpha)[{\bf r}^t {\bf J}_i^{(\tau \upsilon)}{\bf u}_\beta] - J_{i,\alpha \beta}^{(\tau \upsilon)}\}.& (28)}]

2.4.3. Gaussian width

The Gaussian width parameter κ controls radial expansion and contraction of the Cartesian Gaussian multipoles. Analogous parameters are used to optimize the STOs within the Hansen and Coppens scattering model (Hansen & Coppens, 1978[Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909-921.]). The derivative of the electron density with respect to this parameter is similar to the gradient for the ADP parameters. Two chain-rule terms are necessary. Firstly, the gradient of the normalizing term

[{\partial \over {\partial \kappa }}(|{\bf U}_i|^{-{1/2}}) = - {1 \over 2}|{\bf U}_i|^{-3/2} {{\partial |{\bf U}_i|} \over {\partial \kappa}}, \eqno (29)]

where

[\eqalignno {{{\partial |{\bf U}_i|} \over {\partial \kappa }} = &- {{3b_i^3} \over {256\pi ^6 \kappa ^7 }} - {{b_i^2 (U_{11} + U_{22} + U_{33} + 3U_{\rm add})} \over {16\pi ^4 \kappa ^5 }} \cr &+ \{b_i [U_{12}^2 + U_{13}^2 + U_{23}^2 - U_{11} U_{22} - U_{11} U_{33} - U_{22} U_{33} \cr & - 2U_{\rm add} (U_{11} + U_{22} + U_{33}) - 3U_{\rm add}^2]\}/4\pi ^2 \kappa ^3. & (30)}]

Secondly, the gradient of the inverse ADP matrix is most conveniently expressed using the gradient of the original ADP matrix,

[{{\partial {\bf U}_i^{-1}} \over {\partial \kappa }} = - {\bf U}_i^{-1} {{\partial {\bf U}_i } \over {\partial \kappa }}{\bf U}_i^{-1}, \eqno (31)]

where

[{{\partial {\bf U}_i } \over {\partial \kappa }} = - {{b_i } \over {4\pi ^2 \kappa ^3 }}{\bf I}_3. \eqno (32)]

For convenience the matrix Ji(κ) is defined to more compactly represent this result,

[{\bf J}_i^{(\kappa)} = {{b_i} \over {4\pi ^2 \kappa ^3 }} {\bf U}_i^{-1} {\bf U}_i^{-1}. \eqno (33)]

Differentiating (6)[link] with respect to κ and using (29)[link], (33)[link] and the product rule gives

[\eqalignno {{{\partial \rho^{(n,\kappa)}({\bf r})} \over {\partial \kappa }} & = (2\pi)^{- 3/2} {\textstyle \sum\limits_{i = 1}^n} a_i |{\bf U}_i|^{-1/2} \exp(-{\textstyle {1 \over 2}}{\bf r}^t {\bf U}_i^{-1} {\bf r}) \cr &\ \quad {\times}\ {1 \over 2}\left[-{\bf r}^t{\bf J}_i^{(\kappa)} {\bf r} -|{\bf U}_i|^{-1} {{\partial |{\bf U}_i|} \over {\partial \kappa }}\right], \cr {{\partial [\nabla_\alpha \rho^{(n,\kappa)}({\bf r})]} \over {\partial \kappa }} & = - (2\pi)^{-3/2} {\textstyle \sum\limits_{i = 1}^n} a_i|{\bf U}_i|^{-1/2} \exp(-{\textstyle {1 \over 2}}{\bf r}^t {\bf U}_i^{-1} {\bf r}) \cr &\ \quad {\times}\ \biggr\{{1 \over 2}\left[-{\bf r}^t {\bf J}_i^{(\kappa)} {\bf r} - |{\bf U}_i|^{-1} {{\partial |{\bf U}_i|} \over {\partial \kappa }} \right] \cr &\ \quad {\times}\ ({\bf r}^t {\bf U}_i^{-1} {\bf u}_\alpha) + [{\bf r}^t {\bf J}_i^{(\kappa)}{\bf u}_\alpha]\biggr\}, \cr {{\partial[\nabla_\alpha \nabla_\beta \rho^{(n,\kappa)}({\bf r})]} \over {\partial \kappa }} & = (2\pi)^{-3/2} {\textstyle \sum\limits_{i = 1}^n} a_i |{\bf U}_i|^{-1/2} \exp(-{\textstyle {1 \over 2}}{\bf r}^t {\bf U}_i^{-1} {\bf r}) \cr &\ \quad {\times}\ \biggr\{ {1\over 2}\left[- {\bf r}^t {\bf J}_i^{(\kappa)} {\bf r} - |{\bf U}_i|^{-1} {{\partial |{\bf U}_i|} \over {\partial \kappa}} \right]\cr &\ \quad {\times}\ [({\bf r}^t {\bf U}_i^{-1} {\bf u}_\alpha)({\bf r}^t {\bf U}_i^{-1} {\bf u}_\beta) - U_{i,\alpha \beta}^{-1}] \cr &\ \quad +\ [{\bf r}^t {\bf J}_i^{(\kappa)}{\bf u}_\alpha]({\bf r}^t {\bf U}_i^{-1}{\bf u}_\beta) \cr &\ \quad +\ ({\bf r}^t {\bf U}_i^{-1} {\bf u}_\alpha)[{\bf r}^t {\bf J}_i^{(\kappa)}{\bf u}_\beta] - J_{i,\alpha \beta}^{(\tau \upsilon)} \biggr\}, & (34)}]

together with the third- and fourth-order terms available as supplementary information1.

2.5. Fourier transform of the polarizable atomic multipole electron density

Remarkably, the FT of the anisotropic and aspherical density given in (17)[link] is simply

[\eqalignno {{\hat \rho}_{{\rm ADP},j}({\bf s}) =&\ \{P_j^{(c)}\hat f^{(6,1)}_j ({\bf s}) + [P_j^{(v)} - q_j] \hat f_j^{(6,\kappa _v)}({\bf s}) \cr &\ -\ (d_{j,\alpha} + u_{j,\alpha})2\pi is_\alpha \hat f_j^{(1,\kappa _d)}({\bf s}) \cr &\ +\ \textstyle {1 \over 3}\Theta_{j,\alpha \beta} 4\pi ^2 s_\alpha s_\beta \hat f_j^{(1,\kappa _\Theta)}({\bf s}) \}\hat t_j ({\bf s}), & (35)}]

where the dipole and quadrupole terms in (35)[link] depend on the FT of the partial derivatives defined in (18)[link]. Through fifth order the reciprocal-space tensors are

[\eqalignno {F[\nabla_\alpha \rho^{(n,\kappa)}({\bf r})]({\bf s}) & = - 2\pi is_\alpha \hat f^{(n,\kappa)}({\bf s})\hat t({\bf s}) \cr F[\nabla _\alpha \nabla _\beta \rho^{(n,\kappa)} ({\bf r})]({\bf s}) & = - 4\pi ^2 s_\alpha s_\beta \hat f^{(n,\kappa)}({\bf s})\hat t({\bf s}) \cr F[\nabla _\alpha \nabla _\beta \nabla _\gamma \rho^{(n,\kappa)} ({\bf r})]({\bf s}) & = 8\pi ^3 is_\alpha s_\beta s_\gamma \hat f^{(n,\kappa)}({\bf s})\hat t({\bf s}) \cr F[\nabla _\alpha \nabla _\beta \nabla _\gamma \nabla _\delta \rho^{(n,\kappa)} ({\bf r})]({\bf s}) & = 16\pi ^4 s_\alpha s_\beta s_\gamma s_\delta \hat f^{(n,\kappa)}({\bf s})\hat t({\bf s}) \cr F[\nabla _\alpha \nabla _\beta \nabla _\gamma \nabla _\delta \nabla _\varepsilon \rho^{(n,\kappa)} ({\bf r})]({\bf s}) & = - 32\pi ^5 is_\alpha s_\beta s_\gamma s_\delta s_\varepsilon \hat f^{(n,\kappa)}({\bf s})\hat t({\bf s}) & (36)}]

and in compressed tensor notation the general expression for order u + v + w is

[F[\nabla _x^u \nabla _y^v \nabla _z^w \rho^{(n,\kappa)} ({\bf r})]({\bf s}) = (-2\pi i)^{u + v + w} s_a^u s_b^v s_c^w \hat f^{(n,\kappa)}({\bf s})\hat t({\bf s}). \eqno (37)]

This expression is considerably more compact than any reported previously for an aspherical scattering factor in reciprocal space, particularly the formulation based on STOs and spherical harmonics (Hansen & Coppens, 1978[Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909-921.]). Notably, our formulation has no dependence on cumbersome Fourier Bessel transforms of Slater-type functions (Dawson, 1967a[Dawson, B. (1967a). Proc. R. Soc. London Ser. A, 298, 264-288.]; Hansen & Coppens, 1978[Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909-921.]; Su & Coppens, 1990[Su, Z. & Coppens, P. (1990). J. Appl. Cryst. 23, 71-73.]). Our equation (35)[link] has been implemented by `direct summation' for com­parison to the performance of the FFT algorithm.

3. Scattering models

Four scattering models were implemented by modifying and combining the CNS (Brünger et al., 1998[Brünger, A. T., Adams, P. D., Clore, G. M., DeLano, W. L., Gros, P., Grosse-Kunstleve, R. W., Jiang, J.-S., Kuszewski, J., Nilges, M., Pannu, N. S., Read, R. J., Rice, L. M., Simonson, T. & Warren, G. L. (1998). Acta Cryst. D54, 905-921.]) and TINKER (Ponder, 2004[Ponder, J. W. (2004). TINKER: Software Tools for Molecular Design, v.4.2.]) code bases. The scattering models were added to the CNS code base, while TINKER was used to compute AMOEBA chemical forces and to supply CNS with polarizable multipoles in the global frame.

3.1. Isolated atom

The first scattering model (`IAM') is the conventional IAM based on the relativistic elastic scattering factors described by Su & Coppens (1998[Su, Z. & Coppens, P. (1998). Acta Cryst. A54, 357.]).

3.2. Isolated atom with inter-atomic scattering

The second scattering model (`IAM–IAS') augments the IAM with inter-atomic scattering sites at bond centers (Afonine et al., 2007[Afonine, P. V., Grosse-Kunstleve, R. W., Adams, P. D., Lunin, V. Y. & Urzhumtsev, A. (2007). Acta Cryst. D63, 1194-1197.]). Unlike the model of Afonine and coworkers, our implementation does not include IAS sites at lone pairs or at the center of aromatic rings. We have neglected these sites based on the rationale that the AMOEBA electrostatic model is sufficient to capture these details of the electron density, which we provide further evidence for below when discussing the refinement of a Tyr-Gly-Gly tripeptide.

In our approach, chemically equivalent bonds are constrained to use the same IAS parameters. Charge density that is added to or removed from bond centers is exactly balanced by changing the net charge of the bond-defining atoms. For example, a bond charge of −0.2 e requires atomic charge increments that sum to 0.2 e. In this way, all molecules retain their original net charge. Each bond type requires three parameters: the charge increments of both atoms and the Gaussian width of the scattering site. Bond types are defined based on the concatenation of the AMOEBA force-field atom types.

3.3. AMOEBA

The third scattering model (`AMOEBA') is based on the polarizable atomic multipoles of the AMOEBA force field. Each chemically unique multipole type requires three Gaussian width parameters as described in §[link]2. The induced dipoles were iterated to self-consistency using PME whenever any atomic coordinates were changed during refinement (Darden et al., 1993[Darden, T., York, D. & Pedersen, L. (1993). J. Chem. Phys. 98, 10089-10092.]; Sagui et al., 2004[Sagui, C., Pedersen, L. G. & Darden, T. A. (2004). J. Chem. Phys. 120, 73-87.]; Essmann et al., 1995[Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H. & Pedersen, L. G. (1995). J. Chem. Phys. 103, 8577-8593.]).

3.4. AMOEBA with inter-atomic scattering

The final scattering model (`AMOEBA–IAS') augments AMOEBA electrostatics with inter-atomic scattering sites. It became clear during the course of this study that an atomic multipole expansion truncated at quadrupole order is insufficient to capture bond charge density for most molecular geometries. This is consistent with theoretical observations by Stone and coworkers that the convergence of a distributed multipole analysis (DMA) may be improved by using both atoms and bond centers as expansion sites (Stone & Alderton, 1985[Stone, A. J. & Alderton, M. (1985). Mol. Phys. 56, 1047-1064.]; Stone, 2005[Stone, A. J. (2005). J. Chem. Theory Comput. 1, 1128-1132.]). Furthermore, experimental data from the X-ray scattering of diamond and silicon, simple examples of tetrahedral bonding geometry, are explained by the superposition of one atomic octopole moment and one atomic hexadecapole moment (Dawson, 1967a[Dawson, B. (1967a). Proc. R. Soc. London Ser. A, 298, 264-288.],b[Dawson, B. (1967b). Proc. R. Soc. London Ser. A, 298, 379-394.]). The characteristics of the four scattering models are further clarified below with respect to four peptide test cases.

The following computational details were constant across all of the refinements. The isotropic ADP offset Uadd was set to 1/(4π2), which is equivalent to Badd = 8π2Uadd = 2, the FFT grid factor to 0.33 (as appropriate for crystal structures at sub­atomic resolution), and the electron-density cutoff around each atom was 18 (specified by the Elim parameter in CNS). These conservative parameters led to close agreement between direct summation and FFT computation of structure factors. The CNS parameter wA that controls the weighting of X-ray target function relative to the force-field energy was set to 1.0, although we also tested 0.2.

[E_{\rm total} = w_{\rm A} E_{{\rm X}{\hbox {-}} {\rm ray}} + E_{\rm force\,\,field}. \eqno (38)]

This raised Rfree values by less than 0.1% and lowered the AMOEBA potential energy differences between refinements presented below, but did not alter any trends or our conclusions. It should be noted that force-field restraints are not necessarily required for refinement at subatomic high resolution. However, their use in this study gives an insight into the relative energetic cost of the structural changes arising from differences in the four scattering models. A modified version of the refine.inp CNS task file was used for all refinements using the MLI target function.

4. Applications

To demonstrate the behavior of X-ray refinements based on Cartesian Gaussian multipoles, we present two sets of applications. The first set is simply to illustrate the performance of direct summation versus FFT and SGFFT computation of structure factors as a function of system size. The second set describes refinements on a series of four peptide crystals that diffract to 0.59 Å resolution or better. All examples use the AMOEBA force field for chemical forces, instead of the default CNS force field based on Engh & Huber parameters (Engh & Huber, 1991[Engh, R. A. & Huber, R. (1991). Acta Cryst. A47, 392-400.]). Although the refinements were performed in the native space group of each crystal, AMOEBA energies and gradients as computed using the TINKER code base required expanding to P1. This did not increase the number of refined variables, but suggests the need for an AMOEBA code that takes advantage of crystal symmetry.

4.1. Runtime scaling on protein data sets

Evaluation of the target function and its derivatives by direct-summation calculation of structure factors via (35)[link] and (36)[link] is O(Natoms × Nreflections × Nsymm). Alternatively, the FFT algorithm based on (17)[link] and (18)[link] is O(Ngrid × logNgrid), where the number of grid points Ngrid depends on the resolution of the diffraction data. Aspherical refinements based on the Hansen–Coppens formalism are currently limited to direct summation, since the real-space form of the electron density convolved with ADPs is unknown. Therefore, the performance of X-ray refinements based on Cartesian Gaussian multipoles and FFT is of particular interest. The results are summarized in Table 1[link] and are plotted in Fig. 2[link]. Although the performance difference is only about a factor of two for the small protein crambin, over an order of magnitude improvement is achieved for both ribonuclease A and aldose reductase. Parallelization with the SGFFT method results in a further significant speedup (a speedup of a factor of nearly four relative to a single processor on a four-processor machine).

Table 1
Comparison of computational efficiency of direct-summation, FFT and SGFFT methods for the computation of the Cartesian Gaussian multipole scattering factors

The permanent multipole expansion was truncated at atomic quadrupoles and polarization was included via induced dipoles. The FFT method shows a speedup factor of 1.8–14.5 relative to direct summation. Parallelization by SGFFT provided an additional factor of 3.7–3.9 using four processors. All calculations were performed on a MacPro workstation with 2 × 2.66 GHz Dual Core Intel Xeon processors.

PDB code Atoms Reflections Nsymm Atoms × reflections × Nsymm × 10−6 Direct (s) FFT (s) Direct/FFT SGFFT (s) Direct/SGFFT
1ejg 642 112209 2 144.1 49.9 28.1 1.8 7.3 6.8
2vb1 2544 187165 1 476.1 301.8 91.5 3.3 23.6 12.8
1fn8 4294 158550 1 680.8 245.1 45.8 5.4 12.4 19.8
1dy5 4835 159422 2 1541.6 505.6 37.0 13.7 9.7 52.1
1us0 6865 511265 2 7019.7 2346.2 162.3 14.5 42.3 55.5
[Figure 2]
Figure 2
The scaling of the Cartesian Gaussian multipole model, truncated at quadrupole order, is plotted on a log–log scale for computation of the intensity-based maximum-likelihood target function (MLI) for direct summation, FFT and SGFFT. Direct summation scales linearly with the product of the number of atoms, the number of reflections and the number of symmetry operators. Computation of the crystallographic target function by FFT of the Cartesian Gaussian multipole electron density shows a speedup of a factor of between 1.8 and 14.5 compared with direct summation. A further speedup factor of nearly four is achieved using the SGFFT method on a four-processor machine.

4.2. Refinement of peptide crystals

In principle, a more precise scattering model based on Cartesian Gaussian multipoles with coefficients from the AMOEBA electrostatics model should improve the quality of refinements relative to the IAM as judged by both Rfree and the potential energy of the asymmetric unit. Furthermore, the quality of the AMOEBA potential energy function can also be assayed, since it is reasonable to expect that potential energy and Rfree should be correlated.

The peptide crystals studied include YG2 (Pichon-Pesme et al., 2000[Pichon-Pesme, V., Lachekar, H., Souhassou, M. & Lecomte, C. (2000). Acta Cryst. B56, 728-737.]), cyclic P2A4 (Dittrich et al., 2002[Dittrich, B., Koritsánszky, T., Grosche, M., Scherer, W., Flaig, R., Wagner, A., Krane, H. G., Kessler, H., Riemer, C., Schreurs, A. M. M. & Luger, P. (2002). Acta Cryst. B58, 721-727.]) and AYA with three waters or with an ethanol molecule (Chęcińska, Forster et al., 2006[Chęcińska, L., Förster, D., Morgenroth, W. & Luger, P. (2006). Acta Cryst. C62, o454-o457.]; Chęcińska, Mebs et al., 2006[Chęcińska, L., Mebs, S., Hubschle, C. B., Forster, D., Morgenroth, W. & Luger, P. (2006). Org. Biomol. Chem. 4, 3242-3251.]). Detailed descriptions of the unit-cell parameters, number of atoms, resolution and measured reflections are given in Table 2[link]. The refinement results are summarized in Table 3[link] and compared with previous work below.

Table 2
Refinement systems

Molecule Space group and unit-cell parameters (Å, °) Non-H atoms H atoms Bonds dmin (Å) Reflections
YG2 P212121, a = 7.98, b = 9.54, c = 18.32 22 19 40 0.43 4766
P2A4 P212121, a = 10.13, b = 12.50, c = 19.50 35 36 72 0.37 24878
AYA + 3 waters P21, a = 8.12, b = 9.30, c = 12.53, β = 91.21 26 27 50 0.59 5019
AYA + ethanol P21, a = 8.85, b = 9.06, c = 12.36, β = 94.56 26 27 52 0.59 5258

Table 3
Refinement statistics and the relative AMOEBA potential energy per asymmetric unit are given for four small peptide crystals using the IAM, IAM–IAS, AMOEBA and AMOEBA–IAS scattering models

In all cases, the lowest Rfree was found using the AMOEBA–IAS scattering model. Furthermore, the structure with the lowest AMOEBA potential energy per asymmetric unit also corresponded to AMOEBA–IAS refinement.

        Rwork/Rfree (%)  
Molecule Scattering model Nparam Ndata/Nparam Iobs/σ(Iobs) > 0 Iobs/σ(Iobs) > 3 Energy (kcal mol−1)
YGG IAM 274 17.4 4.73/4.74 4.41/4.60 36.5
  IAM–IAS 349 13.7 3.93/4.01 3.59/3.86 7.2
  AMOEBA 355 13.4 4.50/4.56 4.16/4.37 6.8
  AMOEBA–IAS 430 11.1 3.54/3.72 3.17/3.50 0.0
PPAAAA IAM 339 73.4 4.25/4.22 3.65/3.73 32.2
  IAM–IAS 417 59.7 3.56/3.48 3.00/3.01 18.3
  AMOEBA 417 59.7 4.24/4.23 3.69/3.77 12.9
  AMOEBA–IAS 495 50.3 3.42/3.42 2.86/2.94 0.0
AYA + 3 waters IAM 342 14.7 2.75/2.79 2.67/2.71 17.5
  IAM–IAS 411 12.2 2.24/2.47 2.16/2.39 4.1
  AMOEBA 423 11.9 2.40/2.55 2.31/2.47 4.7
  AMOEBA–IAS 492 10.2 1.72/2.03 1.64/1.95 0.0
AYA + ethanol IAM 342 15.4 3.30/3.50 3.20/3.33 23.1
  IAM–IAS 423 12.4 2.32/2.66 2.21/2.49 14.8
  AMOEBA 435 12.1 3.42/3.75 3.32/3.58 3.7
  AMOEBA–IAS 516 10.2 1.90/2.25 1.79/2.08 0.0
†1 kcal mol−1 = 4.186 kJ mol−1.
4.2.1. YG2

The Rfree values of the IAM and IAM–IAS refinements of YG2 (4.60 and 3.86%, respectively) are slightly lower than those reported by Afonine and coworkers (4.72 and 4.06%, respectively; Afonine et al., 2007[Afonine, P. V., Grosse-Kunstleve, R. W., Adams, P. D., Lunin, V. Y. & Urzhumtsev, A. (2007). Acta Cryst. D63, 1194-1197.]). The Rfree value of the AMOEBA–IAS refinement (3.50%) is a significant improvement. The Rwork value (3.17%) of the AMOEBA–IAS refinement is also lower and is comparable to multipolar refinements reported by Volkov and coworkers using transferred or refined multipole coefficients (3.66% and 3.42%, respectively; Volkov et al., 2007[Volkov, A., Messerschmidt, M. & Coppens, P. (2007). Acta Cryst. D63, 160-170.]). Cross-validation-based comparisons are unavailable in this case. We note that the AMOEBA–IAS refinement used a reflections-to-parameters ratio of 11.1, which is slightly higher than the value of 10.6 reported by Volkov and coworkers using refined multipole coefficients. This is computed based on the number of reflections reported in Table 2[link] and the number of parameters given in Table 3[link].

Electron-density maps of the tyrosine ring for the four scattering models are shown in Fig. 3[link], which lend visual insight into their properties. The non-H atom positions are apparent in the 2FoFc contours for each refinement. The standard IAM scattering model underestimates the electron density at bond centers and at the oxygen lone-pair sites, as shown by the FoFc con­tours. Our IAM–IAS scattering model explains the electron density at bond centers, but does not capture lone-pair electron density. Conversely, the AMOEBA model places electron density approximately at the lone-pair positions but not at bond centers. Finally, the AMOEBA–IAS model explains much of the lone-pair and bonding electron densities.

[Figure 3]
Figure 3
(a) IAM, (b) IAM–IAS, (c) AMOEBA and (d) AMOEBA–IAM refinements, respectively, for GY2. The FoFc and 2FoFc σA-weighted electron-density maps are contoured at 3.5σ and shown in green and gray, respectively. Both the IAM and AMOEBA models fail to explain the electron density at bond centers seen in the data. In addition, the IAM model does not account for lone-pair density on the O atom.
4.2.2. P2A4

The Rfree values of our IAM and IAM–IAS refinements of P2A4 (3.73 and 3.01%, respectively) agree closely with the values of Afonine and coworkers (3.63 and 3.23%, respectively; Afonine et al., 2007[Afonine, P. V., Grosse-Kunstleve, R. W., Adams, P. D., Lunin, V. Y. & Urzhumtsev, A. (2007). Acta Cryst. D63, 1194-1197.]). The Rfree value of the AMOEBA–IAS refinement (2.94%) is lower by 0.07%, which is the least amount of improvement seen for AMOEBA–IAS relative to IAM–IAS in this study. The Rwork value (2.86%) of the AMOEBA–IAS refinement is slightly higher, but com­parable to those reported by Volkov and coworkers using transferred or refined multipole coefficients (2.60% and 2.53%), although this work uses a higher reflections-to-parameters ratio (50.3 compared with 43.6; Volkov et al., 2007[Volkov, A., Messerschmidt, M. & Coppens, P. (2007). Acta Cryst. D63, 160-170.]). As for YG2, cross-validation was not performed. The similarity of the R values for YG2 and P2A4 between the AMOEBA–IAS refinements and the multipolar refinements of Volkov and coworkers is consistent with the principle that bond scattering sites capture density that is represented by higher order atomic moments missing in the AMOEBA model (octopole and hexadecapole).

In Fig. 4[link] the precision of the Rwork and Rfree values computed using discrete FTs are compared with analytic direct summation for P2A4 under the AMOEBA scattering model. Agreement to four decimal places is seen for Badd values between 0 and 3 Å2, which serves as validation of the correctness of (17)[link] and (35)[link]. These results support the conclusion that FFT-based computation of structure factors is appropriate for anisotropic and aspherical scattering models.

[Figure 4]
Figure 4
The precision of numerical computation of the Rwork and Rfree values via FFT is compared with analytic direct summation as a function of the isotropic increase Badd in ADP parameters for P2A4 under the AMOEBA scattering model. Note that Badd = 8π2Uadd.
4.2.3. AYA

The AYA data sets were chosen because of the extremely low temperature achieved during the measurement of structure factors (9 K for AYA + three waters and 20 K for AYA + ethanol). For AYA + water, Chęcińska and coworkers (Chęcińska, Forster et al., 2006[Chęcińska, L., Förster, D., Morgenroth, W. & Luger, P. (2006). Acta Cryst. C62, o454-o457.]; Chęcińska, Mebs et al., 2006[Chęcińska, L., Mebs, S., Hubschle, C. B., Forster, D., Morgenroth, W. & Luger, P. (2006). Org. Biomol. Chem. 4, 3242-3251.]) originally reported an R value of 2.4%, which is in agreement with the R value of our IAM refinement (2.67%). Addition of IAS lowered the Rfree statistic from 2.71% to 2.39%, while addition of polarizable atomic multipole electron density showed a further improvement to an Rfree of 1.95%. For AYA + ethanol the Rwork value of the IAM (3.20%) is comparable to that reported originally by Chęcińska and coworkers (2.9%). IAM–IAS lowered Rfree from 3.33 to 2.49%, while AMOEBA–IAS achieved 2.08%.

4.3. Refinement summary

The results for all four peptide refinements are summarized in Fig. 5[link]. In every case, use of the AMOEBA–IAS scattering model relative to the IAM scattering model lowered both Rfree and the potential energy of the crystal. When the IAM scattering model is used, molecular conformations are highly strained to compensate. For example, H—C atom bonds are too short because the IAM model centers electron density at the hydrogen nucleus. In the crystal structures, this electron density is shifted towards the C atom. As the description of the electron density is improved, the molecular conformation relaxes by approximately 16 kJ mol−1 per residue. The precise amount of relaxation depends on the weighting between the crystallographic target and the force field. Unrestrained refinements with an IAM scattering model could adopt even more unphysical conformations. This suggests that accurate chemical restraints are necessary even for ultrahigh-resolution refinements unless an anisotropic and aspherical scattering model is used.

[Figure 5]
Figure 5
The improvement arising from the AMOEBA–IAS scattering model, relative to the IAM model, is plotted as a function of relative percentage improvement in Rfree value and the relative AMOEBA potential energy per residue. For all data sets, the best Rfree value and lowest potential energy per residue were achieved using the AMOEBA–IAS scattering model. 1 kcal mol−1 = 4.186 kJ mol−1.

In Fig. 6[link], we present plots of the IAS sites that were refined for each peptide system. Their Gaussian full-width at half-maximum (FWHM) is plotted against charge magnitude for both the IAM–IAS and the AMOEBA–IAS models. The majority of the charges under the IAM–IAS model and all of the charges under the AMOEBA–IAS model refined to negative partial charge values (or positive scattering density), which is consistent with the physical concentration of charge density at chemical bonds. The similarity of the refined charges between the IAM–IAS and the AMOEBA–IAS models suggests that an atomic multipole description of electron density truncated at quadrupole order underestimates density at trigonal and tetrahedral bond centers.

[Figure 6]
Figure 6
For the inter-atomic scattering sites of the IAM–IAS (a) and AMOEBA–IAS (b) scattering models, the refined Gaussian full-width at half-maximum (FWHM) is plotted versus partial charge magnitude. The majority of charges for the IAM–IAS model and all charges for the AMOEBA–IAS are negative. The sub-angstrom FWHM values are consistent with very localized bond densities.

5. Conclusions

Cartesian Gaussian multipoles offer an efficient alternative to the Hansen and Coppens formulation of aspherical scattering. They eliminate the use of Slater-type functions and allow structure factors to be computed by FFT. Numerical tests show that that FFT and direct-summation implementations of Cartesian Gaussian multipoles agree to high precision. For subatomic resolution biomolecular data sets such as ribo­nuclease A and aldose reductase, parallelized computation of structure factors using the SGFFT method results in a speedup of one to two orders of magnitude compared with direct summation.

Ideally, a force-field electrostatics model should be accurate enough to explain the electron density observed in X-ray diffraction experiments. Although the AMOEBA polarizable multipole force field energetic model shows promise, truncation of the permanent moments at quadrupole order systematically underestimates electron density at bond centers. Our results suggest that the added computational expense of including hexadecapole moments in the atomic scattering factor computation is justified. As supplementary information we have provided a Mathematica notebook and formulae that allow computation of Cartesian Gaussian multipoles up to the fourth order in anticipation of further improvements to force fields.

In the near future, we will present the results of applying our polarizable atomic multipole refinement methodology to macromolecules. For ultrahigh-resolution macromolecular data sets, such as HEWL at 0.65 Å (Wang et al., 2007[Wang, J., Dauter, M., Alkire, R., Joachimiak, A. & Dauter, Z. (2007). Acta Cryst. D63, 1254-1268.]), our scattering model significantly improves refinement statistics, as it does for the simpler peptide cases presented here. Equally exciting will be the use of the AMOEBA force field and in particular the electrostatic forces to orient water molecules in the absence of clear H-atom electron density. We anticipate that refinement of hydrogen-bonding networks will enhance the usefulness of X-ray crystallography experiments with respect to explaining pKa shifts, ligand-binding affinities and enzymatic mechanisms.

APPENDIX A

Fourier transform definition

The definition and notation for the Fourier transform as used in this work is given by

[\eqalignno {\hat f(k) & = \textstyle \int\limits_{- \infty}^\infty f(t)\exp (2\pi ikt) \,{\rm d}t\cr & = F[f(t)](k) & (39)}]

and the corresponding inverse Fourier transform by

[\eqalignno {f(t) & = \textstyle \int\limits_{-\infty}^\infty \hat f(k)\exp (-2\pi ikt) \,{\rm d}k \cr & = F^{-1} [\hat f(k)](t). &(40)}]

APPENDIX B

Derivative of the polarizable electron density with respect to atomic coordinates

The total polarizable electron density arising from the induced dipole of all atoms is given by

[\rho_{\rm ind}({\bf r}) = \textstyle \sum\limits_{i = 1}^n {\bf u}_{i,\beta } \nabla_{i,\beta } \rho _i^{(1,\kappa _d)}({\bf r}). \eqno (41)]

The gradient of this density with respect to the α component of atom j is

[\eqalignno {{{\partial \rho _{\rm ind} ({\bf r})} \over {\partial r_{j,\alpha} }} & = {\partial \over {\partial r_{j,\alpha } }}\textstyle \sum\limits_{i = 1}^n u_{i,\beta} \nabla_{i,\beta } \rho _i^{(1,\kappa _d)} ({\bf r}) \cr & = {\textstyle \sum\limits_{i = 1}^n} \left[{{\partial u_{i,\beta}} \over {\partial r_{j,\alpha} }}\nabla_{i,\beta} \rho _i^{(1,\kappa _d)}({\bf r}) \right] + \delta_{ij} u_{i,\beta} \nabla_\alpha \nabla_\beta \rho_i^{(1,\kappa _d)}({\bf r}). & (42)}]

The second term is nonzero only for i = j and is simple to calculate. The first term, however, depends on ∂ui,β/∂rj,α which is the derivative of a component of the induced dipole of atom i with respect to the α component of atom j. In other words, perturbing the position of atom j affects not only its own scattering but that of all polarizable atoms. The induced dipole on atom i arises from the self-consistent crystal field multiplied by the polarizability,

[\eqalignno {{\bf u}_i & = \alpha_i {\bf E}_i \cr & = \alpha_i \left[\textstyle\sum\limits_{k \ne i} {\bf T}_{ik}^{(1)} {\bf M}_k + \sum\limits_{k \ne i} {\bf T}_{ik}^{(11)} {\bf u}_k \right], & (43)}]

where αi is the atomic polarizability of atom i, Tik(1) is a matrix of tensors that produces the field at site i

[{\bf T}_{ik}^{(1)} = \left(\matrix{ \displaystyle{\partial \over {\partial x_i }} & \displaystyle{{\partial ^2 } \over {\partial x_i \partial x_k }} & \displaystyle{{\partial ^2 } \over {\partial x_i \partial y_k }} & \displaystyle{{\partial ^2 } \over {\partial x_i \partial z_k }} & \ldots \cr\displaystyle {\partial \over {\partial y_i }} & \displaystyle{{\partial ^2 } \over {\partial y_i \partial x_k }} & \displaystyle{{\partial ^2 } \over {\partial y_i \partial y_k }} & \displaystyle{{\partial ^2 } \over {\partial y_i \partial z_k }} & \ldots \cr \displaystyle{\partial \over {\partial z_i }} & \displaystyle{{\partial ^2 } \over {\partial z_i \partial x_k }} & \displaystyle{{\partial ^2 } \over {\partial z_i \partial y_k }} & \displaystyle{{\partial ^2 } \over {\partial z_i \partial z_k }} & \ldots} \right){1 \over {r_{ik}}} \eqno (44)]

owing to the multipole Mk at site k

[{\bf M}_k = (q_k, d_{k,x}, d_{k,y}, d_{k,z}, \Theta _{k,xx}, \Theta _{k,xy}, \Theta _{k,xz}, \ldots,\Theta _{k,zz})^t \eqno (45)]

and Tik(11) is the matrix of tensors that produces the field at site i

[{\bf T}_{ik}^{(11)} = \left(\matrix{ \displaystyle{{\partial ^2 } \over {\partial x_i \partial x_k }} & \displaystyle{{\partial ^2 } \over {\partial x_i \partial y_k }} & \displaystyle{{\partial ^2 } \over {\partial x_i \partial z_k }} \cr \displaystyle{{\partial ^2 } \over {\partial y_i \partial x_k }} & \displaystyle{{\partial ^2 } \over {\partial y_i \partial y_k }} & \displaystyle{{\partial ^2 } \over {\partial y_i \partial z_k }} \cr \displaystyle{{\partial ^2 } \over {\partial z_i \partial x_k }} & \displaystyle{{\partial ^2 } \over {\partial z_i \partial y_k }} & \displaystyle{{\partial ^2 } \over {\partial z_i \partial z_k }}} \right){1 \over {r_{ik}}} \eqno (46)]

owing to the induced dipole uk at site k. For simplicity, we have not formulated (43)[link] using PME electrostatics. Therefore, the sum over k includes all atoms in the crystal except atom i. The derivative of (43)[link] with respect to coordinate rj,α is given by

[\eqalignno {{{\partial u_i } \over {\partial r_{j,\alpha} }} & = \alpha_i \left\{ {\textstyle \sum\limits_{k \ne i}^n} \biggr[{{ \partial {\bf T}_{ik}^{(1)} } \over {\partial r_{j,\alpha} }} {\bf M}_k + {\bf T}_{ik}^{(1)} {{\partial {\bf M}_k } \over {\partial r_{j,\alpha} }}\right] \cr &\ \quad +\ {\textstyle \sum\limits_{k \ne i}^n} \left[{{\partial {\bf T}_{ik}^{(11)} }\over {\partial r_{j,\alpha} }} {\bf u}_k + {\bf T}_{ik}^{(11)} {{\partial {\bf u}_k } \over {\partial r_{j,\alpha} }} \right] \biggr\} \cr & = \alpha_i \biggr[\delta_{jk} {{\partial {\bf T}_{ik}^{(1)}} \over {\partial r_{j,\alpha} }} {\bf M}_k + {\textstyle \sum\limits_{k = \{j\}}} {\bf T}_{ik}^{(1)} {{\partial {\bf M}_k } \over {\partial r_{j,\alpha} }} \cr &\ \quad +\ \delta_{jk} {{\partial {\bf T}_{ik}^{(11)} } \over {\partial r_{j,\alpha } }} {\bf u}_k + {\textstyle \sum\limits_{k \ne i}^n} {\bf T}_{ik}^{(11)} {{\partial {\bf u}_k } \over {\partial r_{j,\alpha } }} \biggr].&(47)}]

The first three terms on the right-hand side are not difficult to compute. However, the fourth term shows that the gradients of the polarizable scattering are O(n3) without use of PME. Specifically, there are 3n × 3n induced dipole density derivatives, each of which is the sum of 3n terms. In this work, we have computed these derivatives by finite differences using PME, which is O(n2logn).

Supporting information


Footnotes

1Supplementary material has been deposited in the IUCr electronic archive (Reference: DZ5164 ). Services for accessing this material are described at the back of the journal.

Acknowledgements

The authors wish to thank Jay W. Ponder and Chuanjie Wu for carefully editing the manuscript. We also thank Pengyu Ren, Paul D. Adams and Thomas A. Darden for helpful discussions. This work was supported by an award from the NSF to Vijay S. Pande, Jay W. Ponder, Teresa Head-Gordon and Martin Head-Gordon for `Collaborative Research: Cyberinfrastructure for Next Generation Biomolecular Modeling' (Award No. CHE-0535675) and by the Howard Hughes Medical Institute.

References

First citationAdams, P. D., Grosse-Kunstleve, R. W., Hung, L.-W., Ioerger, T. R., McCoy, A. J., Moriarty, N. W., Read, R. J., Sacchettini, J. C., Sauter, N. K. & Terwilliger, T. C. (2002). Acta Cryst. D58, 1948–1954.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationAdams, P. D., Pannu, N. S., Read, R. J. & Brünger, A. T. (1997). Proc. Natl Acad. Sci. USA, 94, 5018–5023.  CrossRef CAS PubMed Web of Science Google Scholar
First citationAfonine, P. V., Grosse-Kunstleve, R. W., Adams, P. D., Lunin, V. Y. & Urzhumtsev, A. (2007). Acta Cryst. D63, 1194–1197.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationAfonine, P. V. & Urzhumtsev, A. (2004). Acta Cryst. A60, 19–32.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationAgarwal, R. C. (1978). Acta Cryst. A34, 791–809.  CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationApplequist, J. (1989). J. Phys. A, 22, 4303–4330.  CrossRef Web of Science Google Scholar
First citationApplequist, J. (2002). Theor. Chem. Acc. 107, 103–115.  CrossRef CAS Google Scholar
First citationBerman, H. M., Westbrook, J., Feng, Z., Gilliland, G., Bhat, T. N., Weissig, H., Shindyalov, I. N. & Bourne, P. E. (2000). Nucleic Acids Res. 28, 235–242.  Web of Science CrossRef PubMed CAS Google Scholar
First citationBoys, S. F. (1950). Proc. R. Soc. London Ser. A, 200, 542–554.  CrossRef CAS Google Scholar
First citationBricogne, G. (2001). International Tables for Crystallography, Vol. B, edited by U. Schmueli, pp. 25–98. Dordrecht: Kluwer Academic Publishers.  Google Scholar
First citationBrünger, A. T. (1988). J. Mol. Biol. 203, 803–816.  CAS PubMed Web of Science Google Scholar
First citationBrünger, A. T. (1989). Acta Cryst. A45, 42–50.  CrossRef Web of Science IUCr Journals Google Scholar
First citationBrünger, A. T. (1991). Annu. Rev. Phys. Chem. 42, 197–223.  Google Scholar
First citationBrunger, A. T. (2007). Nature Protoc. 2, 2728–2733.  Web of Science CrossRef CAS Google Scholar
First citationBrünger, A. T., Adams, P. D., Clore, G. M., DeLano, W. L., Gros, P., Grosse-Kunstleve, R. W., Jiang, J.-S., Kuszewski, J., Nilges, M., Pannu, N. S., Read, R. J., Rice, L. M., Simonson, T. & Warren, G. L. (1998). Acta Cryst. D54, 905–921.  Web of Science CrossRef IUCr Journals Google Scholar
First citationBrünger, A. T., Adams, P. D. & Rice, L. M. (1997). Structure, 5, 325–336.  CrossRef CAS PubMed Web of Science Google Scholar
First citationBrünger, A. T., Karplus, M. & Petsko, G. A. (1989). Acta Cryst. A45, 50–61.  CrossRef IUCr Journals Google Scholar
First citationBrünger, A. T., Krukowski, A. & Erickson, J. W. (1990). Acta Cryst. A46, 585–593.  CrossRef Web of Science IUCr Journals Google Scholar
First citationBrünger, A. T., Kuriyan, J. & Karplus, M. (1987). Science, 235, 458–460.  PubMed Web of Science Google Scholar
First citationBrünger, A. T. & Rice, L. M. (1997). Methods Enzymol. 277, 243–269.  PubMed CAS Web of Science Google Scholar
First citationChęcińska, L., Förster, D., Morgenroth, W. & Luger, P. (2006). Acta Cryst. C62, o454–o457.  Web of Science CSD CrossRef IUCr Journals Google Scholar
First citationChęcińska, L., Mebs, S., Hubschle, C. B., Forster, D., Morgenroth, W. & Luger, P. (2006). Org. Biomol. Chem. 4, 3242–3251.  Web of Science PubMed Google Scholar
First citationCooley, J. W. & Tukey, J. W. (1965). Math. Comput. 19, 297–301.  CrossRef Web of Science Google Scholar
First citationCoppens, P. (2005). Angew. Chem. Int. Ed. 44, 6810–6811.  Web of Science CrossRef CAS Google Scholar
First citationCoppens, P. & Volkov, A. (2004). Acta Cryst. A60, 357–364.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationCromer, D. T., Larson, A. C. & Stewart, R. F. (1976). J. Chem. Phys. 65, 336–349.  CSD CrossRef CAS Web of Science Google Scholar
First citationDarden, T., York, D. & Pedersen, L. (1993). J. Chem. Phys. 98, 10089–10092.  CrossRef CAS Web of Science Google Scholar
First citationDawson, B. (1967a). Proc. R. Soc. London Ser. A, 298, 264–288.  CrossRef CAS Google Scholar
First citationDawson, B. (1967b). Proc. R. Soc. London Ser. A, 298, 379–394.  CrossRef CAS Google Scholar
First citationDawson, B. (1967c). Proc. R. Soc. London Ser. A, 298, 395–401.  CrossRef CAS Google Scholar
First citationDeMarco, J. J. & Weiss, R. J. (1965). Phys. Rev. 137, A1869.  CrossRef Google Scholar
First citationDeng, J. P., Xiong, Y. & Sundaralingam, M. (2001). Proc. Natl Acad. Sci. USA, 98, 13665–13670.  Web of Science CrossRef PubMed CAS Google Scholar
First citationDittrich, B., Koritsánszky, T., Grosche, M., Scherer, W., Flaig, R., Wagner, A., Krane, H. G., Kessler, H., Riemer, C., Schreurs, A. M. M. & Luger, P. (2002). Acta Cryst. B58, 721–727.  Web of Science CSD CrossRef CAS IUCr Journals Google Scholar
First citationEngh, R. A. & Huber, R. (1991). Acta Cryst. A47, 392–400.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationEpstein, J., Bentley, J. & Stewart, R. F. (1977). J. Chem. Phys. 66, 5564–5567.  CrossRef CAS Web of Science Google Scholar
First citationEssmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H. & Pedersen, L. G. (1995). J. Chem. Phys. 103, 8577–8593.  CrossRef CAS Web of Science Google Scholar
First citationFriesner, R. A., Baldwin, R. L. & Baker, D. (2005). Adv. Protein Chem. 72, 79–104.  Web of Science CrossRef PubMed Google Scholar
First citationGresh, N. (2006). Curr. Pharm. Des. 12, 2121–2158.  Web of Science CrossRef PubMed CAS Google Scholar
First citationGresh, N., Cisneros, G. A., Darden, T. A. & Piquemal, J. P. (2007). J. Chem. Theory Comput. 3, 1960–1986.  Web of Science CrossRef PubMed CAS Google Scholar
First citationGrosse-Kunstleve, R. W. & Adams, P. D. (2002). J. Appl. Cryst. 35, 477–480.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationHansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909–921.  CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationHehre, W. J., Ditchfie, R., Stewart, R. F. & Pople, J. A. (1970). J. Chem. Phys. 52, 2769–2773.  CrossRef CAS Web of Science Google Scholar
First citationHehre, W. J., Stewart, R. F. & Pople, J. A. (1969). J. Chem. Phys. 51, 2657–2664.  CrossRef CAS Web of Science Google Scholar
First citationHoward, E. I., Sanishvili, R., Cachau, R. E., Mitschler, A., Chevrier, B., Barth, P., Lamour, V., Van Zandt, M., Sibley, E., Bon, C., Moras, D., Schneider, T. R., Joachimiak, A. & Podjarny, A. (2004). Proteins, 55, 792–804.  Web of Science CrossRef PubMed CAS Google Scholar
First citationJelsch, C., Guillot, B., Lagoutte, A. & Lecomte, C. (2005). J. Appl. Cryst. 38, 38–54.  Web of Science CrossRef IUCr Journals Google Scholar
First citationJelsch, C., Teeter, M. M., Lamzin, V., Pichon-Pesme, V., Blessing, R. H. & Lecomte, C. (2000). Proc. Natl Acad. Sci. USA, 97, 3171–3176.  Web of Science CrossRef PubMed CAS Google Scholar
First citationKielkopf, C. L., Ding, S., Kuhn, P. & Rees, D. C. (2000). J. Mol. Biol. 296, 787–801.  Web of Science CrossRef PubMed CAS Google Scholar
First citationKuhn, P., Knapp, M., Soltis, S. M., Ganshaw, G., Thoene, M. & Bott, R. (1998). Biochemistry, 37, 13446–13452.  Web of Science CrossRef CAS PubMed Google Scholar
First citationKuriyan, J., Brünger, A. T., Karplus, M. & Hendrickson, W. A. (1989). Acta Cryst. A45, 396–409.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationKurki-Suonio, K. (1968). Acta Cryst. A24, 379–390.  CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationLamoureux, G., Harder, E., Vorobyov, I. V., Roux, B. & MacKerell, A. D. (2006). Chem. Phys. Lett. 418, 245–249.  Web of Science CrossRef CAS Google Scholar
First citationLamoureux, G. & Roux, B. (2003). J. Chem. Phys. 119, 3025–3039.  Web of Science CrossRef CAS Google Scholar
First citationMaple, J. R., Cao, Y. X., Damm, W. G., Halgren, T. A., Kaminski, G. A., Zhang, L. Y. & Friesner, R. A. (2005). J. Chem. Theory Comput. 1, 694–715.  Web of Science CrossRef CAS Google Scholar
First citationMuzet, N., Guillot, B., Jelsch, C., Howard, E. & Lecomte, C. (2003). Proc. Natl Acad. Sci. USA, 100, 8742–8747.  Web of Science CrossRef PubMed CAS Google Scholar
First citationPande, V. S., Baker, I., Chapman, J., Elmer, S. P., Khaliq, S., Larson, S. M., Rhee, Y. M., Shirts, M. R., Snow, C. D., Sorin, E. J. & Zagrovic, B. (2003). Biopolymers, 68, 91–109.  Web of Science CrossRef PubMed CAS Google Scholar
First citationPatel, S. & Brooks, C. L. (2006). Mol. Simul. 32, 231–249.  Web of Science CrossRef CAS Google Scholar
First citationPichon-Pesme, V., Lachekar, H., Souhassou, M. & Lecomte, C. (2000). Acta Cryst. B56, 728–737.  Web of Science CSD CrossRef CAS IUCr Journals Google Scholar
First citationPiquemal, J. P., Chelli, R., Procacci, P. & Gresh, N. (2007). J. Phys. Chem. A, 111, 8170–8176.  Web of Science CrossRef PubMed CAS Google Scholar
First citationPiquemal, J. P., Cisneros, G. A., Reinhardt, P., Gresh, N. & Darden, T. A. (2006). J. Chem. Phys. 124, 104101.  Web of Science CrossRef PubMed Google Scholar
First citationPonder, J. W. (2004). TINKER: Software Tools for Molecular Design, v.4.2.  Google Scholar
First citationPonder, J. W. & Case, D. A. (2003). Adv. Protein Chem. 66, 27–85.  CrossRef PubMed CAS Google Scholar
First citationRen, P. Y. & Ponder, J. W. (2002). J. Comput. Chem. 23, 1497–1506.  Web of Science CrossRef PubMed CAS Google Scholar
First citationRen, P. Y. & Ponder, J. W. (2003). J. Phys. Chem. B, 107, 5933–5947.  Web of Science CrossRef CAS Google Scholar
First citationRen, P. Y. & Ponder, J. W. (2004). J. Phys. Chem. B, 108, 13427–13437.  Web of Science CrossRef CAS Google Scholar
First citationSagui, C., Pedersen, L. G. & Darden, T. A. (2004). J. Chem. Phys. 120, 73–87.  Web of Science CrossRef PubMed CAS Google Scholar
First citationSchnieders, M. J., Baker, N. A., Ren, P. Y. & Ponder, J. W. (2007). J. Chem. Phys. 126, 124114.  Web of Science CrossRef PubMed Google Scholar
First citationSchnieders, M. J. & Ponder, J. W. (2007). J. Chem. Theory Comput. 3, 2083–2097.  Web of Science CrossRef CAS Google Scholar
First citationShirts, M. R. & Pande, V. S. (2005). J. Chem. Phys. 122, 134508.  Web of Science CrossRef PubMed Google Scholar
First citationShirts, M. R., Pitera, J. W., Swope, W. C. & Pande, V. S. (2003). J. Chem. Phys. 119, 5740–5761.  Web of Science CrossRef CAS Google Scholar
First citationSnow, C. D., Nguyen, N., Pande, V. S. & Gruebele, M. (2002). Nature (London), 420, 102–106.  Web of Science CrossRef PubMed CAS Google Scholar
First citationSnow, C. D., Sorin, E. J., Rhee, Y. M. & Pande, V. S. (2005). Annu. Rev. Biophys. Biomol. Struct. 34, 43–69.  Web of Science CrossRef PubMed CAS Google Scholar
First citationSorin, E. J. & Pande, V. S. (2005). J. Comput. Chem. 26, 682–690.  Web of Science CrossRef PubMed CAS Google Scholar
First citationStewart, R. F. (1977). Chem. Phys. Lett. 49, 281–284.  CrossRef CAS Web of Science Google Scholar
First citationStewart, R. F. (1979). Chem. Phys. Lett. 65, 335–342.  CrossRef CAS Web of Science Google Scholar
First citationStone, A. J. (1996). The Theory of Intermolecular Forces. Oxford: Clarendon Press.  Google Scholar
First citationStone, A. J. (2005). J. Chem. Theory Comput. 1, 1128–1132.  Web of Science CrossRef CAS Google Scholar
First citationStone, A. J. & Alderton, M. (1985). Mol. Phys. 56, 1047–1064.  CrossRef CAS Web of Science Google Scholar
First citationSu, Z. & Coppens, P. (1990). J. Appl. Cryst. 23, 71–73.  CrossRef Web of Science IUCr Journals Google Scholar
First citationSu, Z. & Coppens, P. (1998). Acta Cryst. A54, 357.  Web of Science CrossRef IUCr Journals Google Scholar
First citationTen Eyck, L. F. (1973). Acta Cryst. A29, 183–191.  CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationTen Eyck, L. F. (1977). Acta Cryst. A33, 486–492.  CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationTereshko, V., Wilds, C. J., Minasov, G., Prakash, T. P., Maier, M. A., Howard, A., Wawrzak, Z., Manoharan, M. & Egli, M. (2001). Nucleic Acids Res. 29, 1208–1215.  Web of Science CrossRef PubMed CAS Google Scholar
First citationTrueblood, K. N., Bürgi, H.-B., Burzlaff, H., Dunitz, J. D., Gramaccioli, C. M., Schulz, H. H., Shmueli, U. & Abrahams, S. C. (1996). Acta Cryst. A52, 770–781.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationVolkov, A., Messerschmidt, M. & Coppens, P. (2007). Acta Cryst. D63, 160–170.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationWaller, I. (1923). Z. Phys. A, 17, 398–408.  CrossRef CAS Google Scholar
First citationWang, J., Dauter, M., Alkire, R., Joachimiak, A. & Dauter, Z. (2007). Acta Cryst. D63, 1254–1268.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationZarychta, B., Pichon-Pesme, V., Guillot, B., Lecomte, C. & Jelsch, C. (2007). Acta Cryst. A63, 108–125.  Web of Science CrossRef CAS IUCr Journals Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.

Journal logoBIOLOGICAL
CRYSTALLOGRAPHY
ISSN: 1399-0047
Follow Acta Cryst. D
Sign up for e-alerts
Follow Acta Cryst. on Twitter
Follow us on facebook
Sign up for RSS feeds