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

IUCrJ
Volume 3| Part 4| July 2016| Pages 237-246
ISSN: 2052-2525

Quantum crystallographic charge density of urea

CROSSMARK_Color_square_no_text.svg

aComputer, Computational, and Statistical Sciences Division, Los Alamos National Laboratory, Mail Stop B256, Los Alamos, New Mexico 87545, USA
*Correspondence e-mail: mewall@lanl.gov

Edited by A. Fitch, ESRF, France (Received 24 November 2015; accepted 13 April 2016; online 8 June 2016)

Standard X-ray crystallography methods use free-atom models to calculate mean unit-cell charge densities. Real molecules, however, have shared charge that is not captured accurately using free-atom models. To address this limitation, a charge density model of crystalline urea was calculated using high-level quantum theory and was refined against publicly available ultra-high-resolution experimental Bragg data, including the effects of atomic displacement parameters. The resulting quantum crystallographic model was compared with models obtained using spherical atom or multipole methods. Despite using only the same number of free parameters as the spherical atom model, the agreement of the quantum model with the data is comparable to the multipole model. The static, theoretical crystalline charge density of the quantum model is distinct from the multipole model, indicating the quantum model provides substantially new information. Hydrogen thermal ellipsoids in the quantum model were very similar to those obtained using neutron crystallography, indicating that quantum crystallography can increase the accuracy of the X-ray crystallographic atomic displacement parameters. The results demonstrate the feasibility and benefits of integrating fully periodic quantum charge density calculations into ultra-high-resolution X-ray crystallographic model building and refinement.

1. Introduction

Efforts to increase the accuracy of charge density models from X-ray crystallography have mainly focused on fitting the Bragg data using functions that are more expressive than the usual free-atom spherical distributions. Stewart (1969[Stewart, R. F. (1969). J. Chem. Phys. 51, 4569-4577.]) proposed using general scattering factors that are the products of atom-centered orbital wavefunctions, and restrictions to better match the number of free parameters to the number of reflections in fitting (Stewart, 1970[Stewart, R. F. (1970). J. Chem. Phys. 53, 205-213.]). Coppens et al. (1971[Coppens, P., Pautler, D. & Griffin, J. F. (1971). J. Am. Chem. Soc. 93, 1051-1058.]) separated the free atom charge density into core and valence components, and allowed them to be centered on different positions. Dawson decomposed the charge into symmetric and antisymmetric components centered on each atom (Dawson, 1967a[Dawson, B. (1967a). Proc. R. Soc. A: Math. Phys. Eng. Sci. 298, 255-263.]), and expanded each atom-centered charge density in spherical harmonics (Dawson, 1967b[Dawson, B. (1967b). Proc. R. Soc. A: Math. Phys. Eng. Sci. 298, 264-288.]). Hirshfeld developed a least-squares method that models aspherical atomic charge densities using basis functions related to spherical harmonics, but with alternative symmetry properties (Hirshfeld, 1971[Hirshfeld, F. L. (1971). Acta Cryst. B27, 769-781.]). Spherical harmonic-related methods were integrated into multipole refinement computer programs that are used when charge density models are desired (Hansen & Coppens, 1978[Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909-921.]; Hirshfeld, 1977a[Hirshfeld, F. L. (1977a). Isr. J. Chem. 16, 226-229.]; Craven & Weber, 1977[Craven, B. M. & Weber, H.-P. (1977). Technical Report, Department of Crystallography, University of Pittsburgh.]; Stewart & Spackman, 1983[Stewart, R. F. & Spackman, M. A. (1983). Valray User's Manual. Carnegie Mellon University, USA, and University of Copenhagen, Denmark.]; Jelsch et al., 2005[Jelsch, C., Guillot, B., Lagoutte, A. & Lecomte, C. (2005). J. Appl. Cryst. 38, 38-54.]; Volkov et al., 2006[Volkov, A., Macchi, P., Farrugia, L. J., Gatti, C., Mallinson, P., Richter, T. & Koritsanszky, T. (2006). XD2006 - A Computer Program Package for Multipole Refinement, Topological Analysis of Charge Densities and Evaluation of Intermolecular Energies from Experimental and Theoretical Structure Factors. University of New York at Buffalo, USA.]).

Although less well exploited than multipole methods, the potential for combining quantum theory and X-ray diffraction to obtain accurate charge density models of molecular crystals has been long appreciated (Lipscomb, 1972[Lipscomb, W. (1972). Trans. Am. Crystallogr. Assoc. 8, 79-92.]). This combination has been termed quantum crystallography (Massa et al., 1995[Massa, L., Huang, L. & Karle, J. (1995). Int. J. Quantum Chem. 56, 371-384.]). The high computational cost of quantum electronic structure calculations has been a major barrier to exploiting the theory for crystallography; however, recent linear scaling methods have made calculations possible for large systems (Bowler & Miyazaki, 2010[Bowler, D. R. & Miyazaki, T. (2010). J. Phys. Condens. Matter, 22, 074207.], 2012[Bowler, D. R. & Miyazaki, T. (2012). Rep. Prog. Phys. 75, 036503.]; Goedecker, 1999[Goedecker, S. (1999). Rev. Mod. Phys. 71, 1085-1123.]; VandeVondele et al., 2012[VandeVondele, J., Borštnik, U. & Hutter, J. (2012). J. Chem. Theory Comput. 8, 3565-3573.]), and fast quantum molecular dynamics simulations for systems approaching 104 atoms with 105 time steps are now possible (Mniszewski et al., 2015[Mniszewski, S. M., Cawkwell, M. J., Wall, M. E., Mohd-Yusof, J., Bock, N., Germann, T. C. & Niklasson, A. M. N. (2015). J. Chem. Theory Comput. 11, 4644-4654.]). Methods using quantum theory to calculate crystallographic charge density models for all but the largest systems therefore might soon be within reach, not only for small-molecule crystallography (Capelli et al., 2014[Capelli, S. C., Bürgi, H.-B., Dittrich, B., Grabowsky, S. & Jayatilaka, D. (2014). IUCrJ, 1, 361-379.]) but also for macromolecular crystallography.

Several methods have been proposed for quantum crystallography, including the method of kernel projector matrices (Massa et al., 1995[Massa, L., Huang, L. & Karle, J. (1995). Int. J. Quantum Chem. 56, 371-384.]) and fitting of wavefunctions to diffraction data (Jayatilaka, 1998[Jayatilaka, D. (1998). Phys. Rev. Lett. 80, 798-801.]). One method that is showing promise in practical applications is Hirshfeld Atom Refinement (HAR) (Bruning & Feil, 1992[Bruning, H. & Feil, D. (1992). Acta Cryst. A48, 865-872.]; Capelli et al., 2014[Capelli, S. C., Bürgi, H.-B., Dittrich, B., Grabowsky, S. & Jayatilaka, D. (2014). IUCrJ, 1, 361-379.]; Jayatilaka & Dittrich, 2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.]). In HAR, the static charge density of a molecule is calculated using quantum theory and is partitioned into individual atom contributions using Hirshfeld's stockholder method (Hirshfeld, 1977b[Hirshfeld, F. L. (1977b). Theor. Chim. Acta, 44, 129-138.]). The partitioned charge is used to calculate aspherical atomic structure factors that are substituted for the usual structure factors in crystallographic refinement, considering both the atomic positions and displacement parameters (Bruning & Feil, 1992[Bruning, H. & Feil, D. (1992). Acta Cryst. A48, 865-872.]). Whereas Bruning & Feil (1992[Bruning, H. & Feil, D. (1992). Acta Cryst. A48, 865-872.]) originally decomposed the charge density into individual atom contributions using a multipole expansion; the more recent implementation of Jayatilaka & Dittrich and coworkers (Capelli et al., 2014[Capelli, S. C., Bürgi, H.-B., Dittrich, B., Grabowsky, S. & Jayatilaka, D. (2014). IUCrJ, 1, 361-379.]; Jayatilaka & Dittrich, 2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.]) directly makes use of a Becke grid for individual atom charge densities. The HAR method has been automated to apply iterative updates of the quantum electronic structure calculation during refinement of atomic positions (Capelli et al., 2014[Capelli, S. C., Bürgi, H.-B., Dittrich, B., Grabowsky, S. & Jayatilaka, D. (2014). IUCrJ, 1, 361-379.]). So far HAR has been limited to gas-phase electronic structure calculations, with cluster charges placed at symmetry-related positions to approximate the crystal environment.

Whether HAR or other quantum crystallography methods will be adopted widely depends critically on whether they will substantially increase the accuracy of X-ray crystallography models. To date the main focus on the accuracy of HAR has been whether it yields molecular geometry and atomic displacement parameters that are consistent with neutron crystallography. The results here have been promising: applications to X-ray diffraction from crystalline benzene and urea (Jayatilaka & Dittrich, 2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.]), L-phenyl­alaninium hydrogen maleate (Woińska et al., 2014[Woińska, M., Jayatilaka, D., Spackman, M. A., Edwards, A. J., Dominiak, P. M., Woźniak, K., Nishibori, E., Sugimoto, K. & Grabowsky, S. (2014). Acta Cryst. A70, 483-498.]), and a Gly L-Ala dipeptide (Capelli et al., 2014[Capelli, S. C., Bürgi, H.-B., Dittrich, B., Grabowsky, S. & Jayatilaka, D. (2014). IUCrJ, 1, 361-379.]) found that HAR bond distances agreed very well with neutron crystal structures, overcoming known deficiencies in spherical-atom charge density models (Lipscomb, 1972[Lipscomb, W. (1972). Trans. Am. Crystallogr. Assoc. 8, 79-92.]). Atomic displacement parameters from HAR similarly agreed reasonably well with the neutron crystal structures.

This study addresses a major factor that so far has been lacking in evaluating quantum crystallographic methods: the accuracy of the charge density model. Here, charge density models for crystalline urea are obtained using spherical atom, atomic multipole or quantum methods. For the quantum method, the HAR method (Jayatilaka & Dittrich, 2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.]) is adapted for crystalline phase electronic structure calculations performed using VASP (Kresse & Furthmüller, 1996[Kresse, G. & Furthmüller, J. (1996). Phys. Rev. B, 54, 11169-11186.]). Electronic structure calculations using VASP previously were performed on hexachlorobenzene for comparisons to the X-ray crystallographic multipole charge density (Aubert et al., 2011[Aubert, E., Lebègue, S., Marsman, M., Bui, T. T. T., Jelsch, C., Dahaoui, S., Espinosa, E. & Ángyán, J. G. (2011). J. Phys. Chem. A, 115, 14484-14494.]), but without allowing for individual ADPs. The novel aspect of the present method therefore is the combination of a crystalline phase density-functional-theory-based electronic structure calculation with an atomic displacement model from HAR. The results indicate that HAR can yield not only molecular geometries and ADPs that are similar to the neutron crystal structure, but also both 2FoFc maps and static charge densities that are distinct from the multipole model, but that nevertheless agree comparably with the experimental data. Quantum crystallography therefore can yield accurate charge densities that are consistent simultaneously with theory and experiment.

2. Methods

2.1. Diffraction data and initial crystal structure

Ultra-high-resolution urea synchrotron diffraction data were obtained from Birkedal et al. (2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]) at //journals.iucr.org/a/issues/2004/05/00/xc5013/xc5013Isup7.hkl. These data were collected at a temperature of 123 K using a wavelength of 0.5996 (1) Å, and were merged into 1045 unique reflections (992 positively valued) extending to 0.347 Å resolution. The data were consistent with a [P\bar 42_1m] unit cell (space group 113), with a = b = 5.5780 (6), c = 4.6860 (7) Å, α = β = γ = 90°. Other data collection details are published in Birkedal et al. (2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]), Table 1[link]. The multipole refined urea crystal structure was obtained from Birkedal et al. (2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]) at //journals.iucr.org/a/issues/2004/05/00/xc5013/xc5013sup1.cif. The hydrogen parameters of this model were copied from a 123 K neutron crystal structure (Swaminathan et al., 1984[Swaminathan, S., Craven, B. M. & McMullan, R. K. (1984). Acta Cryst. B40, 300-306.]).

Table 1
Values of crystallographic agreement factors for: spherical atom (sphere); quantum starting with the spherical atom (Quant-S), Birkedal et al. (2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]) multipole (Quant-M), and neutron (Quant-N) model; and multipole (Birkedal et al., 2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]) model

  wR2F wR2I RF RSR GooF
Sphere 0.034 0.068 0.038 0.042 6.5
Quant-S 0.009 0.018 0.023 0.021 1.7
Quant-M 0.009 0.018 0.023 0.021 1.7
Quant-N 0.009 0.018 0.023 0.021 1.7
Multipole 0.011 0.021 0.026 0.023 2.0

2.2. Spherical atom and multipole crystallographic models

The program SHELXL (Sheldrick, 2008[Sheldrick, G. M. (2008). Acta Cryst. A64, 112-122.], 2015[Sheldrick, G. M. (2015). Acta Cryst. C71, 3-8.]), Version 2014/1, was used to refine a spherical atom model of urea. Atomic coordinates and anisotropic atomic displacement parameters (ADPs) were refined for all atoms, in addition to an overall scale factor (27 parameters in all). The experimental temperature of 123 K was selected for geometry restraints. SHELXL reported agreement factors for the final model are: R1 = 0.0370, wR2 = 0.0796, and SHELX goodness of fit = 0.639 for all reflections. Mean unit cell charge-density maps Fo, 2FoFc, and FoFc were calculated using the program shelx2map provided in the SHELX software distribution, using the refined .fcf file as the input, with default weighting, yielding a map of dimensions 56 × 15 × 42 for the asymmetric unit. The maps were expanded to P1 using CCP4 (Stein et al., 1994[Stein, P. E., Boodhoo, A., Armstrong, G. D., Cockle, S. A., Klein, M. H. & Read, R. J. (1994). Structure, 2, 45-57.]) mapmask and were interpolated to a 64 × 64 × 64 grid using CCP4 maprot.

The program MoPro, Version 14.06 (Guillot et al., 2001[Guillot, B., Viry, L., Guillot, R., Lecomte, C. & Jelsch, C. (2001). J. Appl. Cryst. 34, 214-223.]; Jelsch et al., 2005[Jelsch, C., Guillot, B., Lagoutte, A. & Lecomte, C. (2005). J. Appl. Cryst. 38, 38-54.]), was used to refine a multipole charge-density model of urea. Refinement was based on the ultra-high-resolution data and structure from Birkedal et al. (2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]). 20 cycles of automated density refinement were performed using the REFI DENS method. The total number of free parameters was 37: one scale factor; five valence (VAL); five κ1 (K1); five κ2 (K2); 21 Plm (PLM) multipole parameters. The atom coordinates and ADPs were kept constant. MoPro reported agreement factors for the final model are: RF = 0.0242, wR2F = 0.0107, RI = 0.0212, wR2I = 0.0212, and GooF = 2.293 for 992 nonzero reflections. Charge density maps were calculated using the MoPro supplied program VMoPro. The Fo, 2FoFc, and FoFc maps were computed by Fourier reconstruction using the FOUR method using the refined .par file and .FOUR file as inputs, with default resolution limits and the FFT method, yielding maps on a 92 × 92 × 80 grid; these maps were interpolated to a 64 × 64 × 64 grid using CCP4 maprot (Stein et al., 1994[Stein, P. E., Boodhoo, A., Armstrong, G. D., Cockle, S. A., Klein, M. H. & Read, R. J. (1994). Structure, 2, 45-57.]). The total static crystalline charge density was computed using the VMoPro STAT method, with a 10 Å selection for grid limits, grid-cube dimensions in fractional coordinates, the origin at (0,0,0), a maximum coordinate value of 0.9844 in each dimension, a 10 Å margin around the grid for contributing atoms, and 64 × 64 × 64 grid points. MoPro charge densities were scaled to yield a total charge of 64 electrons in the unit cell.

2.3. Quantum crystallographic models

A custom implementation of the original Hirshfeld atom refinement method (Jayatilaka & Dittrich, 2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.]) was used to obtain quantum crystallographic models. Quantum charge density calculations were performed using atomic coordinates from each of three different models: the SHELX refined spherical atom structure; the neutron crystal structure of Swaminathan et al. (1984[Swaminathan, S., Craven, B. M. & McMullan, R. K. (1984). Acta Cryst. B40, 300-306.]); and the multipole model of Birkedal et al. (2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]). An expanded unit cell with 16 atoms was generated using the Computational Crystallography Toolbox (cctbx) (Grosse-Kunstleve et al., 2002[Grosse-Kunstleve, R. W., Sauter, N. K., Moriarty, N. W. & Adams, P. D. (2002). J. Appl. Cryst. 35, 126-136.]) by applying [P\bar 42_1m] symmetry to the five-atom asymmetric unit. Ab initio density functional theory calculations were performed using VASP, Version 5.3.3 (Kresse & Furthmüller, 1996[Kresse, G. & Furthmüller, J. (1996). Phys. Rev. B, 54, 11169-11186.]). Instead of pseudopotentials, the PAW method was used, with PAW_PBE parameters (Kresse & Joubert, 1999[Kresse, G. & Joubert, D. (1999). Phys. Rev. B, 59, 1758-1775.]). The electronic structure was computed using 43 = 64 Monkhorst-Pack k points. Partial occupancies were calculated using Fermi smearing at the experimental temperature of 123 K. As there are fewer than 20 atoms in the expanded urea unit cell, LREAL = .FALSE. was used to evaluate projection operators in reciprocal space, as recommended in the VASP documentation. The valence charge density v(x) was calculated for the expanded P1 unit cell. In addition, using the same VASP PAW method as for the molecular calculation, 16 crystalline core charge densities ci(x) and 16 crystalline free-atom (`promolecule') charge densities fi(x ) were obtained for each individual atom i.

To achieve the desired model accuracy, all VASP charge densities were calculated on a 128 × 128 × 128 grid spanning the unit cell. For crystallographic refinement, the densities were decimated to a 64 × 64 × 64 grid. The decimated total static charge density calculated in VASP is provided in the supporting information, along with the difference between the VASP and MoPro multipole static charge density.

X-ray structure factors were calculated using both new and existing tools in Lunus software (Wall, 2009[Wall, M. E. (2009). Methods Mol. Biol. 544, 269-279.]), which was originally designed for analysis and modeling of diffuse X-ray scattering data (Wall et al., 1997a[Wall, M. E., Clarage, J. B. & Phillips, G. N. (1997a). Structure, 5, 1599-1612.],b[Wall, M. E., Ealick, S. E. & Gruner, S. M. (1997b). Proc. Natl Acad. Sci. USA, 94, 6180-6184.], 2014[Wall, M. E., Van Benschoten, A. H., Sauter, N. K., Adams, P. D., Fraser, J. S. & Terwilliger, T. C. (2014). Proc. Natl Acad. Sci. USA, 111, 17887-17892.]). The effect of ADPs was modeled using a Stockholder method (Bruning & Feil, 1992[Bruning, H. & Feil, D. (1992). Acta Cryst. A48, 865-872.]; Hirshfeld, 1977b[Hirshfeld, F. L. (1977b). Theor. Chim. Acta, 44, 129-138.]). The total valence density v(x) was partitioned into atomic contributions using the equation

[{v_i}\left(x \right) = {q_i}\left(x \right)v\left(x \right). \eqno (1)]

Hirshfeld partitioning (Hirshfeld, 1977b[Hirshfeld, F. L. (1977b). Theor. Chim. Acta, 44, 129-138.]) was used with weights qi(x ) defined using the free-atom charge density fi(x )

[{q_i}\left(x \right) = {{{f_i}\left(x \right)} \over {\mathop \sum \nolimits_i {f_i}\left(x \right)}}. \eqno (2)]

Similar to Bruning & Feil (1992[Bruning, H. & Feil, D. (1992). Acta Cryst. A48, 865-872.]), ADPs were modeled by treating each partitioned atom charge density ai(x ) = ci(x ) + vi(x ) as a rigid distribution, displaced along with the atom. However, in contrast to Bruning & Feil (1992[Bruning, H. & Feil, D. (1992). Acta Cryst. A48, 865-872.]), instead of using a multipole expansion, the charge density ai(uvw ) was sampled on a rectilinear grid spanning the unit cell, indexed by uvw. This method is similar to that of Jayatilaka & Dittrich (2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.]), who used a radial-angular Becke grid for sampling. Here a rectilinear grid is chosen, as it corresponds precisely both to the VASP results and to the discrete sampling by the Bragg peaks in the crystallographic experiment. The partitioned atom structure factors were defined as Ai(hkl ) = DFT[ai(uvw ) ], where DFT denotes a discrete Fourier transform. The DFT was computed using a fast Fourier transform (FFT) algorithm (Press et al., 1999[Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (1999). Numerical Recipes in C, 2nd ed. Cambridge University Press.]).

In the original Hirshfeld refinement method (Jayatilaka & Dittrich, 2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.]), the value [A^{\prime}_{i}\left({hkl} \right)] of Ai(hkl ) after a coordinate shift [x^{\prime}] was obtained by multiplying Ai(hkl ) by a phase factor. Although multiplication by a phase factor is appropriate for arbitrary translations of a continuous distribution or an atom-centered grid, it is not appropriate for translations by fractional grid points on a fixed rectilinear grid such as is used here. The correct transformation instead requires a resampling of the shifted distribution ai(uvw ) on the original grid (Appendix[link]). The structure factors are obtained by transforming [x^{\prime}] to the grid coordinates [u^{\prime}v^{\prime}w^{\prime}], decomposing these coordinates into integer (u0v0w0) and fractional (u1v1w1) parts such that [0 \le {u_1} \lt 1], [0 \le {v_1} \lt 1], and [0 \le {w_1} \lt 1], and using the following equation to calculate [A^{\prime}_{i}\left({hkl} \right)]

[\eqalignno { A{^\prime _i}\left({hkl} \right) &= {A_i}\left({hkl} \right){e^{{{ - 2\pi ih{u_0}} \over {{N_1}}}}}{e^{{{ - 2\pi ik{v_0}} \over {{N_2}}}}}{e^{{{ - 2\pi il{w_0}} \over {{N_3}}}}} \cr & \times \biggl [(1 - {u_1})(1 - {v_1})(1 - w_1) \,+ \, e^{{{- 2\pi ih}\over{N_1}}} u_1 (1 - v_1)(1 - w_1) \cr & + e^{{{- 2\pi ik}\over{N_2}}}(1 - u_1)v_1(1 - w_1) \,+ \, e^{{{- 2\pi il}\over{N_3}}} (1 - u_1)(1 - v_1){w_1} \cr & + e^{{{ - 2\pi ih}\over{N_1}}} e^{{{- 2\pi ik}\over{N_2}}} {u_1}{v_1}\left({1 - {w_1}} \right) + e^{{{- 2\pi ih}\over{N_1}}} e^{{{- 2\pi il}\over{N_3}}}u_1\left(1 - {v_1}\right){w_1} \cr & + e^{{{- 2\pi ik}\over{N_2}}} e^{{{- 2\pi il}\over{N_3}}} (1 - u_1){v_1}{w_1} + e^{{{-2\pi ih}\over{N_1}}} e^{{{- 2\pi ik}\over{N_2}}} e^{{{- 2\pi il}\over{N_3}}} u_1 v_1 w_1\biggr].\cr & &(3)}]

The unit-cell structure factor Fc(hkl ) was then calculated as

[{F}_{\rm c}\left(hkl\right) = {\sum }_{i}{A^\prime}_{i}\left(hkl\right){e}^{-{2{\pi }^{2}s}_{hkl}\cdot {U}_{i}\cdot {s}_{hkl}}, \eqno (4)]

where [{e}^{-{2{\pi }^{2}s}_{hkl}\cdot {U}_{i}\cdot {s}_{hkl}}] is the Debye–Waller factor for the matrix Ui of ADPs for atom i, and shkl is the scattering vector corresponding to Miller indices hkl.

2.4. Quantum model refinement

Quantum refinements were performed starting with the spherical atom (S), multipole (M) (Birkedal et al., 2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]), and neutron crystallography (N) (Swaminathan et al., 1984[Swaminathan, S., Craven, B. M. & McMullan, R. K. (1984). Acta Cryst. B40, 300-306.]) atomic coordinates and ADPs. Model refinement was performed by minimizing the goodness-of-fit (GooF) statistic

[GooF = {\left [{{{ \sum \nolimits_{hkl} {{\left({{I_{\rm o}}\left({hkl} \right) - {I_{\rm c}}\left({hkl} \right)} \right)}^2}/\sigma _I^2\left({hkl} \right)} \over {NDF}}} \right]^{1/2}}, \eqno (5)]

where Ic(hkl ) = | Fc(hkl ) |2, Io(hkl ) and [\sigma _I\left({hkl} \right)] are the values and errors of the observed intensities, and the number of degrees of freedom NDF = 965 is the number of data points (= 992 non-negative intensity values in the merged data set), minus the number of free parameters in the fit (= 27, see below). A value of the GooF for each set of coordinates and ADPs was obtained by minimizing it with respect to an arbitrary scale factor between the calculated and observed reflection amplitudes. Each matrix Ui was decomposed into eigenvalues and eigenvectors, and three Euler angles were computed from the eigenvectors, to obtain a set of independent parameters for efficient optimization. Optimization with respect to atom positions and eigenvalues and Euler angles from U matrices was performed in python using the scipy.optimize.minimize Powell method, using default settings. Due to the use of the Powell method, error bars were not obtained for the fitted parameter values. Eigenvalues were constrained to be positive. Symmetry of the atomic coordinates and ADPs was enforced explicitly using the following equations: X = 0 and U23 = 0 for C, O atoms; Y = X+0.5 for all atoms; and U22 = U11 and U13 = U23 for all atoms. Enforcement of symmetry reduced the number of free parameters from 9 to 4 for the C, O atoms and to 6 for the N, H1, and H2 atoms. There were a total of 27 free parameters in the refinement, including the scale factor between the data and model (the same number as for spherical atom refinement, but without geometry restraints).

The mean unit cell charge densities ρ were calculated using Fourier reconstruction as [\rho \left({uvw} \right) = DFT\left [{F\left({hkl} \right)} \right]]. The experimental Fo, 2FoFc and FoFc maps were calculated by applying the model phases to the observations. Values of | Fc(hkl ) | were used in place of missing values of | Fo(hkl ) |. Only complete grids were used for FFT calculations on the quantum models; reflections were not truncated using a resolution cutoff. The Fo, Fc, 2FoFc, and FoFc maps obtained using the multipole model as an input structure are provided in the supporting information.

2.5. Agreement factors

The agreement of all models with the diffraction data was assessed using several standard statistics: GooF, wR2F, wR2I, RF, and RSR. The GooF [equation (5)[link]] was used as the refinement target. The weighted R-squared factor for amplitudes, wR2F, was calculated as

[wR^2F = {\left [{{{ \sum \nolimits_{hkl} {{\left({\left| {{F_{\rm o}}\left({hkl} \right)} \right| - \left| {{F_{\rm c}}\left({hkl} \right)} \right|} \right)}^2}/\sigma _F^2\left({hkl} \right)} \over { \sum \nolimits_{hkl} {{\left| {{F_{\rm o}}\left({hkl} \right)} \right|}^2}/\sigma _F^2\left({hkl} \right)}}} \right]^{1/2}}, \eqno (6)]

where | Fo(hkl ) | and [\sigma _F\left({hkl} \right)] are the experimental reflection amplitudes and errors, and Fc(hkl ) is calculated using equation (4)[link]. The weighted R-squared factor for intensities, wR2I, was calculated as

[wR^2I = {\left [{{{ \sum \nolimits_{hkl} {{\left({{I_{\rm o}}\left({hkl} \right) - {I_{\rm c}}\left({hkl} \right)} \right)}^2}/\sigma _I^2\left({hkl} \right)} \over { \sum \nolimits_{hkl} {I_{\rm o}}\left({hkl} \right)/\sigma _I^2\left({hkl} \right)}}} \right]^{1/2}}\semi \eqno (7)]

the R factor for amplitudes, RF, was calculated as

[RF = {{\sum \nolimits_{hkl} \left| {\left| {{F_{\rm o}}\left({hkl} \right)} \right| - \left| {{F_{\rm c}}\left({hkl} \right)} \right|} \right|} \over {\sum \nolimits_{hkl} \left| {{F_{\rm o}}\left({hkl} \right)} \right|}}; \eqno (8)]

the real-space R-factor, RSR, was calculated as

[RSR = {{ \sum \nolimits_{uvw} \left| {{\rho _{\rm o}}\left({uvw} \right) - {\rho _{\rm c}}\left({uvw} \right)} \right|} \over {\sum \nolimits_{uvw} \left| {{\rho _{\rm o}}\left({uvw} \right) + {\rho _{\rm c}}\left({uvw} \right)} \right|}}, \eqno (9)]

where [{\rho _{\rm o}}] is the experimental Fo map. Calculated and observed values were scaled to minimize the RMSD prior to using equations (5)–(9)[link][link][link][link][link], and both [{\rho _{\rm o}}] and [{\rho _{\rm c}}] were offset to have zero mean prior to using equation (9)[link]. To enable fair comparison, all agreement factors were calculated using Lunus software tools. Values reported in primary references were very similar to those computed using Lunus.

3. Results

The agreement factors for all quantum crystallographic models are the same (in % units): wR2F (target) = 0.9, wR2I = 1.8, RF = 2.3, RSR = 2.1 and GooF = 1.7 (Table 1[link]). These are slightly better than the multipole model, which has values 0.2–0.3% higher for each. The quantum and multipole models agree much better with the data than the spherical atom model (Table 1[link]).

Three-dimensional visualizations of the 2FoFc and FoFc maps for the spherical atom, quantum-M, and multipole models are shown in Fig. 1[link]. The spherical atom and multipole 2FoFc maps appear to be more similar to each other than they are to the quantum model. This appearance is supported quantitatively using a RSR statistic calculated between each pair of 2FoFc real-space maps, using an appropriately modified equation (9)[link]. A value of 4.9% was obtained between the spherical atom and multipole models. By comparison, the RSR values between the quantum model and either the spherical atom (8.2%) or the multipole model (7.4%) were much greater. These values are all higher than the RSR of any of the models with the data (Table 1[link]).

[Figure 1]
Figure 1
Comparison of 2FoFc and FoFc maps for (a) spherical atom, (b) Quantum-M, and (c) multipole models. Level charge density surfaces in 2FoFc maps are rendered using a blue wireframe at a level of 1-sigma. Level surfaces in FoFc maps are shown in green (positive electron density, negative charge) and red (negative electron density, positive charge) wireframes at 3-sigma. The figure was created using COOT (Emsley et al., 2010[Emsley, P., Lohkamp, B., Scott, W. G. & Cowtan, K. (2010). Acta Cryst. D66, 486-501.]).

Visualization of contours in a two-dimensional section including the C=O bond reveals that the quantum and multipole 2FoFc maps are very different (Figs. 2[link]a and b). (Much of this difference might be an artifact in the multipole map calculation, as mentioned below.) Compared with the multipole model (Fig. 2[link]b), the quantum-M model is smoother (Fig. 2[link]a). The multipole model shows ripples surrounding core atoms and peaks away from atoms, including between bonded heavy atoms. The quantum-M model has some peaks away from the atom cores (Fig. 2[link]a), but these are lower in magnitude compared with the multipole model. The quantum-M and multipole model FoFc difference maps are broadly similar (Figs. 2[link]c and d), with the larger deviations from the data in the multipole model along the C=O axis, consistent with the slightly higher values of agreement factors for this model (Table 1[link]).

[Figure 2]
Figure 2
Comparison of two-dimensional contours in the reconstructed mean unit cell charge density for the quantum-M (left), and multipole (right) models in the y = 0 section. (a)–(b) 2FoFc maps in 0.05 e Å−1 contours, to a maximum of 3 e Å−1. (c)–(d) FoFc difference maps in 0.05 e Å−1 contours. Negative electron density contours in each panel are colored red. The view is along the same direction as that in Fig. 4[link], in the plane of the C=O bond. The orientation is such that x increases along the horizontal, and z increases along the vertical. The image was created using mapslicer in the CCP4 suite (CCP4, 1994[Collaborative Computational Project, Number 4 (1994). Acta Cryst. D50, 760-763.]).

To investigate further the differences between the 2FoFc maps of the quantum-M and multipole models, we compared the static total charge densities calculated using either VASP or MoPro. Both charge densities correspond to the multipole geometry (Birkedal et al., 2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]). There are visible differences (Figs. 3[link]a and b); however, the differences are much smaller than in the 2FoFc maps (Figs. 2a and b), and they coincide with atoms and bonds. The comparison suggests that the ripples in the 2FoFc map from the multipole model are an artifact of the FOUR method implementation in VMoPro (e.g. a truncation of reflections beyond 0.347 Å resolution).

[Figure 3]
Figure 3
Two-dimensional contours in the static total charge densities derived from the multipole model atom coordinates (Birkedal et al., 2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]) (the section and orientation is the same as in Fig. 2[link]). (a) Theoretical density computed using VASP. (b) Multipole charge density refinement in MoPro. (c) Density in (b) subtracted from density in (a). (d) Difference density computed by subtracting the theoretical density using quantum-M refined atom coordinates from the density in (a) [the view of the total density using the quantum-M structure is indistinguishable from (a)]. Contours in all panels are in 0.05 e Å−1 intervals, to a maximum of 3 e Å−1. Negative electron density contours in panels (c) and (d) are colored red. The image was created using mapslicer in the CCP4 suite (CCP4, 1994[Collaborative Computational Project, Number 4 (1994). Acta Cryst. D50, 760-763.]).

Subtracting the static charge densities from VASP and MoPro reveals substantial differences in the charge distribution along the C=O bond (Fig. 3[link]c). These differences show a similar pattern of peaks and troughs as in the FoFc map for the multipole model (Fig. 2[link]d); by comparison, the FoFc map of the quantum-M model shows smaller differences along the C=O bond (Fig. 2[link]c). Combined, Figs. 2[link] and 3[link] indicate that the multipole static charge density contains deviations from the data in the C=O bond that are decreased in the quantum-M model.

The VASP and MoPro calculations were further compared using a Bader analysis of the net atom charges (Tang et al., 2009[Tang, W., Sanville, E. & Henkelman, G. (2009). J. Phys. Condens. Matter, 21, 084204.]; Table 2[link]). The theoretical VASP charges are similar for the spherical atom, multipole and neutron structures. The main difference between these and the multipole model charges is for the C atom, which has a value of 4.06–4.08 electrons from the theoretical density, and 4.57 electrons from the multipole model density. This substantial 0.5 electron difference is compensated by smaller differences in the charges on the other atoms, which are between 0.06 and 0.09 electrons smaller in the multipole model.

Table 2
Atom charges based on Bader analysis of total static charge densities

VASP calculations correspond to the spherical atom (S), multipole (M), and neutron (N) structures. The MoPro calculation corresponds to the multipole (M) structure. Units are negative charge in electrons.

  C O N H1 H2
VASP (S) 4.08 9.26 8.33 0.51 0.50
VASP (M) 4.06 9.26 8.31 0.51 0.53
VASP (N) 4.06 9.27 8.31 0.51 0.53
MoPro (M) 4.57 9.20 8.22 0.46 0.44

The atomic coordinates of the quantum-M, quantum-N, neutron and multipole models are all very similar (Table 3[link]). The differences between these models and either the spherical atom or quantum-S model are small for the heavy atoms, but are larger for the H atoms. The differences lead to a substantial deviation in the N—H1 bond for the spherical atom and quantum-S structure compared with the neutron structure (Table 4[link]): the bond length is 1.006 Å in the neutron structure compared with 0.911 Å in the spherical atom and 0.810 Å in the quantum-S structure. The differences also lead to decreases in the C—N—H1 and C—N—H2 bond angles for both the spherical atom and quantum-S structures compared with the neutron structure (Table 5[link]): the angles in the spherical atom model are about 2° smaller, and the angles in the quantum-S model are about 4° smaller than in the neutron structure. There is a corresponding increase in the H1—N—H2 angle for each compared with the neutron structure: 4° for the spherical atom and 8° for the quantum-S structure. The angle deviations for the quantum-S model are visible in the stick diagram in Fig. 3[link](a); the effect is smaller but still perceptible for the spherical atom model (not shown).

Table 3
Coordinates for atoms in the asymmetric unit

Model labels in column 2 are as in Table 1[link], with the addition of the neutron model (Swaminathan et al., 1984[Swaminathan, S., Craven, B. M. & McMullan, R. K. (1984). Acta Cryst. B40, 300-306.]). Hydrogen model parameters for the multipole model were copied from the neutron model (Birkedal et al., 2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]). Units are fractions of unit-cell dimensions.

    X Y Z
C Sphere 0 0.5 0.3281
Quant-S 0 0.5 0.3279
Quant-M 0 0.5 0.3281
Quant-N 0 0.5 0.3277
Neutron 0 0.5 0.3280
Multipole 0 0.5 0.3282
O Sphere 0 0.5 0.5964
Quant-S 0 0.5 0.5967
Quant-M 0 0.5 0.5966
Quant-N 0 0.5 0.5963
Neutron 0 0.5 0.5962
Multipole 0 0.5 0.5963
N Sphere 0.1450 0.6450 0.1783
Quant-S 0.1452 0.6452 0.1782
Quant-M 0.1446 0.6446 0.1796
Quant-N 0.1447 0.6447 0.1785
Neutron 0.1447 0.6447 0.1785
Multipole 0.1447 0.6447 0.1790
H1 Sphere 0.2438 0.7438 0.2793
Quant-S 0.2316 0.7316 0.2722
Quant-M 0.2571 0.7571 0.2837
Quant-N 0.2571 0.7571 0.2838
Neutron 0.2557 0.7557 0.2841
Multipole 0.2557 0.7557 0.2841
H2 Sphere 0.1382 0.6382 −0.0363
Quant-S 0.1324 0.6324 −0.0371
Quant-M 0.1432 0.6432 −0.0343
Quant-N 0.1432 0.6432 −0.0344
Neutron 0.1431 0.6431 −0.0348
Multipole 0.1431 0.6431 −0.0348

Table 4
Bond lengths for alternative models

Model labels are as in Table 1[link]. Units are Å.

  C=O N—O N—H1 N—H2
Sphere 1.257 2.268 0.911 1.006
Quant-S 1.259 2.271 0.810 1.013
Quant-M 1.258 2.263 1.011 1.001
Quant-N 1.258 2.266 1.012 0.996
Neutron 1.257 2.266 1.006 1.000
Multipole 1.257 2.265 1.005 1.002

Table 5
Bond angles for alternative models

Model labels are as in Table 1[link]. Units are degrees.

  O—C—N N—C—N C—N—H1 C—N—H2 H1—N—H2
Sphere 121.54 116.91 117.28 118.42 124.30
Quant-S 121.49 117.04 115.68 115.82 128.50
Quant-M 121.40 117.21 119.87 120.81 119.32
Quant-N 121.49 117.01 119.30 121.02 119.68
Neutron 121.54 116.92 118.99 120.82 120.20
Multipole 121.50 117.01 119.15 120.78 120.07

The similarity of ADPs was assessed using the S statistic, which describes the deviation of the three-dimensional positional distribution of the atoms defined by the ADPs (Whitten & Spackman, 2006[Whitten, A. E. & Spackman, M. A. (2006). Acta Cryst. B62, 875-888.]). The ADPs for the heavy atoms in the quantum models are very similar to the neutron models (Table 6[link]; Fig. 4[link]): the value of S for the C atom ranges from 0.04 to 0.06%; the value for the O atom from 0.13 to 0.14%; and the value for the N atom from 0.27 to 0.32%. The similarities are comparable for the multipole model. The ADPs for the H atoms in the quantum models are also similar to the neutron crystal structure, but to a lesser degree than the heavy atoms: the value of S for H1 varies from 1.61 to 2.25%; and the value for H2 varies from 2.34 to 3.46%. A high-level quantum theoretical calculation of vibrations of urea to obtain ADPs (Madsen et al., 2013[Madsen, A. Ø., Civalleri, B., Ferrabone, M., Pascale, F. & Erba, A. (2013). Acta Cryst. A69, 309-321.]) yielded a comparable similarity for the heavy atoms (S = 0.12, 0.16, and 1.1% for C, O, and N, respectively) and a higher similarity for the H atoms (S = 0.13 and 0.05% for H1 and H2, respectively). However, the similarities in Madsen et al. (2013[Madsen, A. Ø., Civalleri, B., Ferrabone, M., Pascale, F. & Erba, A. (2013). Acta Cryst. A69, 309-321.]) were computed after applying an overall scale factor with respect to the experimental ADPs; the similarities are considerably lower without applying the scale factor (S = 1.36, 0.9, 1.7, 0.5, and 0.65% for C, O, N, H1, and H2, respectively, using the B3LYP/6-31G(d,p) method). The similarity of the spherical atom heavy atom ADPs to the neutron structure is high (S = 0.04, 0.13, and 0.21% for C, O, and N), and the similarity for the H atoms is low, as expected for a spherical atom model (S = 25.65 and 5.91%). The multipole model hydrogen parameters were copied from the neutron structure and therefore are identical (Birkedal et al., 2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]).

Table 6
Values of ADPs for atoms in the asymmetric unit

Labels in column 2 are as in Table 2[link]. Units are Å2. U11 = U22 and U13 = U23 by symmetry. The similarity statistic (S) with respect to the ADPs of the neutron model is computed following Whitten & Spackman (2006[Whitten, A. E. & Spackman, M. A. (2006). Acta Cryst. B62, 875-888.]), in % units. Hydrogen model parameters for the multipole model were copied from the neutron model (Birkedal et al., 2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]).

    U11 = U22 U33 U12 U13 = U23 S
C Sphere 0.0150 0.0070 0.0000 0.0000 0.04
Quant-S 0.0141 0.0061 0.0001 0.0000 0.05
Quant-M 0.0141 0.0061 0.0001 0.0000 0.04
Quant-N 0.0141 0.0060 0.0001 0.0000 0.06
Neutron 0.0147 0.0065 0.0001 0.0000 0.00
Multipole 0.0152 0.0068 −0.0004 0.0000 0.03
O Sphere 0.0199 0.0066 0.0020 0.0000 0.13
Quant-S 0.0194 0.0059 0.0019 0.0000 0.14
Quant-M 0.0194 0.0060 0.0020 0.0000 0.14
Quant-N 0.0194 0.0061 0.0020 0.0000 0.13
Neutron 0.0197 0.0063 0.0001 0.0000 0.00
Multipole 0.0196 0.0067 0.0016 0.0000 0.10
N Sphere 0.0293 0.0096 −0.0155 0.0001 0.21
Quant-S 0.0285 0.0087 −0.0158 0.0000 0.32
Quant-M 0.0285 0.0085 −0.0156 0.0001 0.29
Quant-N 0.0286 0.0087 −0.0156 0.0001 0.27
Neutron 0.0286 0.0095 −0.0147 0.0002 0.00
Multipole 0.0293 0.0096 −0.0157 0.0000 0.23
H1 Sphere 0.0295 0.0468 −0.0127 −0.0184 25.65
Quant-S 0.0550 0.0259 −0.0392 −0.0019 1.61
Quant-M 0.0495 0.0168 −0.0295 0.0019 2.13
Quant-N 0.0490 0.0172 −0.0296 0.0022 2.25
Neutron 0.0440 0.0216 −0.0223 −0.0031 0.00
Multipole 0.0440 0.0216 −0.0223 −0.0031 0.00
H2 Sphere 0.0415 0.0206 0.0099 −0.0024 5.91
Quant-S 0.0380 0.0227 −0.0191 0.0015 2.34
Quant-M 0.0409 0.0270 −0.0187 −0.0013 3.43
Quant-N 0.0410 0.0272 −0.0186 −0.0012 3.46
Neutron 0.0430 0.0141 −0.0159 0.0020 0.00
Multipole 0.0430 0.0141 −0.0159 0.0020 0.00
[Figure 4]
Figure 4
Displacement ellipsoids at 50% probability for (a) Quantum-S, (b) Quantum-M, and (c) neutron diffraction models. In each structure, the O atom is red, the C dark grey, the N blue, and the H1 and H2 light grey, with the H2 at the bottom of the molecule. The Quantum-N model is not shown as it is indistinguishable from the Quantum-M model; similarly, the multipole model is not shown as it is indistinguishable from the neutron. The latter is due to the fact that the hydrogen parameters of the multipole model were copied from the neutron model (Birkedal et al., 2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]). The image was created using Mercury (Macrae et al., 2008[Macrae, C. F., Bruno, I. J., Chisholm, J. A., Edgington, P. R., McCabe, P., Pidcock, E., Rodriguez-Monge, L., Taylor, R., van de Streek, J. & Wood, P. A. (2008). J. Appl. Cryst. 41, 466-470.]).

To assess the convergence of the quantum refinement, as was done in previous HAR implementations (Capelli et al., 2014[Capelli, S. C., Bürgi, H.-B., Dittrich, B., Grabowsky, S. & Jayatilaka, D. (2014). IUCrJ, 1, 361-379.]; Jayatilaka & Dittrich, 2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.]), the electronic structure calculation was iteratively applied to each of the quantum models. In the iteration, the refined atomic coordinates were used to re-compute all charge densities using VASP, and the model was re-refined against the data using the new densities. The quantum-M and quantum-N models were essentially unchanged in the second iteration: the initial static charge densities were very similar to those for the first iteration, as shown for the quantum-M model in Fig. 3[link](d); the agreement factors remained the same as in Table 1[link]; all of the fractional atomic coordinates changed by less than 5 × 10−3, with maximal changes of 1 × 10-3 for heavy atom coordinates; and the similarity statistic for the ADPs was 0.06% or lower for all atoms between the first and second iterations. In contrast, the quantum-S model showed divergent behavior: the agreement factors were slightly larger (by 0.1–0.2%) for the second iteration; hydrogen fractional coordinates changed by as much as 0.05 (a 0.2 Å shift of the x- and y-position of the H1 atom); and the similarity statistic for the ADPs was as high as 2.6% (H1 atom), which is comparable to the value computed between the quantum models and the neutron model (Table 6[link]).

4. Discussion

The agreement of the quantum crystallographic models of urea with ultra-high-resolution data compares favorably to the multipole model. Both the 2FoFc map and the total static charge density are substantially different between the quantum and multipole models, however. The differences in 2FoFc appear largely to be due to an artifact in the multipole map, as they contain ripples that do not coincide with atom positions or bonds (Figs. 2[link]a and b). The differences in the static charge density, however, appear to be real, with notable differences both in the electronic structure of the C=O bond (Fig. 3[link]c) and in the 0.5-electron higher negative charge associated with the C atom for the multipole model (Table 2[link]). The difference in the C atom charge is consistent with the multipole charge density study of Birkedal et al. (2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]), which reported a 0.7–0.8 electron larger negative charge for the C atom in the multipole model compared with theoretical charge density calculations.

Whereas Birkedal et al. (2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]) concluded the difference between their multipole model of urea and the theoretical charge density was due to inaccuracies in the quantum electronic structure calculation, this study suggests that the difference might instead be due to inaccuracies in the multipole model. The quantum charge densities were obtained using quantum theory and are consistent with the experimental data. The multipole model, although also consistent with the experimental data, is more weakly tied to the underlying theory, and relies on the fitting of many parameters. The possibility of inaccuracies in the multipole model is supported by a controlled study using synthetic data (De Vries et al., 2000[Vries, R. Y. de, Feil, D. & Tsirelson, V. G. (2000). Acta Cryst. B56, 118-123.]) which found that the charge density of urea could not be determined uniquely using multipole refinement; however, this support is tempered by the fact that the synthetic data did not extend to a resolution as high as the data in Birkedal et al. (2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]).

The present results indicate that it would be worthwhile investigating whether HAR might produce more reliable interaction density models than are currently obtained using multipole methods. Compared with multipole refinement, HAR uses fewer parameters and relies on quantum theory for the increased expressiveness needed to model the aspherical component of the charge density. In addition, in HAR, the same quantum electronic structure method used for the crystal phase calculation can be used for the gas phase. Thus, whereas calculating the multipole interaction density involves subtracting two densities that were obtained using substantially different methods, the HAR interaction density can be obtained by subtracting densities that are more comparable.

The thermal ellipsoids in the quantum models are both quantitatively (Table 6[link]) and qualitatively (Fig. 4[link]) similar to the neutron crystal structure. This was even the case for the quantum-S model, despite the lack of convergence seen in a second iteration of refinement and deviations in the geometry with respect to the neutron model (Tables 4[link] and 5[link]). This finding is consistent with studies in which the neutron crystallographic temperature factors of urea (Jayatilaka & Dittrich, 2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.]) and other systems (Capelli et al., 2014[Capelli, S. C., Bürgi, H.-B., Dittrich, B., Grabowsky, S. & Jayatilaka, D. (2014). IUCrJ, 1, 361-379.]; Woińska et al., 2014[Woińska, M., Jayatilaka, D., Spackman, M. A., Edwards, A. J., Dominiak, P. M., Woźniak, K., Nishibori, E., Sugimoto, K. & Grabowsky, S. (2014). Acta Cryst. A70, 483-498.]) were found to be reproduced reasonably well using HAR. In particular, the previous urea study (Jayatilaka & Dittrich, 2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.]) used BLYP density functional theory with a cc-pVTZ basis and surrounding charge clusters to mimic periodic boundary conditions, and used the same starting structure as the present quantum-M refinement (Birkedal et al., 2004[Birkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371-381.]). The following ADP values were obtained for H atoms (Jayatilaka & Dittrich, 2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.]) (U11, U33, U12, and U13 in Å2 units): (0.0550, 0.0170, −0.0350, 0) for H1, and (0.0450, 0.0260, −0.0190, −0.0020) for H2. These values are similar to those found here (Table 6[link]); in addition, values for heavy atoms differed by less than 0.001 Å2 compared with those found here. The similarity statistics S computed with respect to the ADPs from the quantum-M model are (in % units): 0.01, 0.03, 0.06, 0.26, and 0.08 for the C, N, O, H1, and H2 atoms, respectively. The similarity of the ADPs here with those in Jayatilaka & Dittrich (2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.]), in addition to the consistency of the quantum refined ADPs in this study using different starting structures, indicates that HAR can yield estimates of ADPs that are robust to differences in starting structures and DFT methods.

The results for the quantum-M and quantum-N models indicate the potential advantages of quantum crystallography for accurate charge density studies. The lack of convergence and geometry deviations of the quantum-S model, however, indicate that challenges remain for the general applicability of these methods. The deviations of the quantum-S model can be traced back to deviations in the spherical atom model: although the heavy atoms are consistent with the neutron structure, the deviation in the hydrogen positions is more substantial (Table 3[link]), leading to corresponding deviations in geometry (Tables 4[link] and 5[link]). These deviations are increased rather than decreased in the quantum-S model, which prevents the quantum refinement from converging on successive iterations.

Because it is not currently feasible to obtain neutron crystal structures for all systems of interest, the generalization of the quantum crystallography methods developed here to routine X-ray crystallographic structure determination will require improved modeling of hydrogen positions in the starting structure. The successful application of iterative electronic structure calculations in HAR applications to ammonia and Gly–L-Ala using spherical atoms models like those used here as input structures (Capelli et al., 2014[Capelli, S. C., Bürgi, H.-B., Dittrich, B., Grabowsky, S. & Jayatilaka, D. (2014). IUCrJ, 1, 361-379.]) indicates that a spherical atom model is adequate for at least some molecular crystals. It would be interesting to determine whether iterative HAR using the implementation of Capelli et al. (2014[Capelli, S. C., Bürgi, H.-B., Dittrich, B., Grabowsky, S. & Jayatilaka, D. (2014). IUCrJ, 1, 361-379.]) converges using the present spherical atom model of urea as an input (as mentioned above, the study of Jayatilaka & Dittrich, 2008[Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383-393.], made use of the same model as the quantum-M refinement here). It is possible that hydrogen positions in spherical atom models would be sufficiently improved using methods that leverage information in structure databases (Bąk et al., 2011[Bąk, J. M., Domagała, S., Hübschle, C., Jelsch, C., Dittrich, B. & Dominiak, P. M. (2011). Acta Cryst. A67, 141-153.]; Bendeif & Jelsch, 2007[Bendeif, E. & Jelsch, C. (2007). Acta Cryst. C63, o361-o364.]; Dadda et al., 2012[Dadda, N., Nassour, A., Guillot, B., Benali-Cherif, N. & Jelsch, C. (2012). Acta Cryst. A68, 452-463.]; Dittrich et al., 2005[Dittrich, B., Hübschle, C. B., Messerschmidt, M., Kalinowski, R., Girnt, D. & Luger, P. (2005). Acta Cryst. A61, 314-320.], 2009[Dittrich, B., Weber, M., Kalinowski, R., Grabowsky, S., Hübschle, C. B. & Luger, P. (2009). Acta Cryst. B65, 749-756.]), which can place H atoms to within O(10−2) Å of the positions in neutron crystal structures.

There are many ways HAR may be extended, targeting, e.g., more accurate models of structure variation than ADPs, and larger systems. For larger systems, it will be especially important to assess the applicability of fast, approximate quantum electronic structure calculations (Mniszewski et al., 2015[Mniszewski, S. M., Cawkwell, M. J., Wall, M. E., Mohd-Yusof, J., Bock, N., Germann, T. C. & Niklasson, A. M. N. (2015). J. Chem. Theory Comput. 11, 4644-4654.]) to quantum crystallography. Increasing the speed of calculations would enable the wider adaptation of HAR by adapting existing small-molecule crystallography workflows, and would provide a complementary quantum crystallographic alternative to multipole refinement (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.]) for obtaining high-resolution charge density models of molecular crystals, including macromolecular crystals.

APPENDIX A

Transformation of the structure factors upon translation of the charge density

Here the transformation is illustrated using the one-dimensional case; the extension to three dimensions in equation (3)[link] is straightforward. The structure factors are given by the discrete Fourier transform (DFT) of a periodic charge density [\rho \left(u \right)] sampled at N fixed grid points u, with [\rho \left({u + N} \right) = \rho \left(u \right)]. Let the charge density [\rho \left(u \right)] correspond to the following continuous step-wise distribution

[\rho \left(x \right) = \sum \limits_{u = 0}^{N - 1} \rho \left(u \right)\left [{\theta \left({x - u} \right) - \theta \left({x - u - 1} \right)} \right], \eqno (10)]

where the Heaviside distribution [\theta \left(x \right) = 1] for [x\ge 0] and [\theta \left(x \right) = 0] otherwise. Translation of [\rho \left(x \right)] by a shift [\Delta x] yields a new charge density

[\eqalignno{{\rho _{\Delta x}}\left(x \right) =\, & \rho \left({x - \Delta x} \right)\cr =\, & \sum \limits_u^{} \rho \left(u \right)\left [{\theta \left({x - \Delta x - u} \right) - \theta \left({x - \Delta x - u - 1} \right)} \right]. \cr & & (11)}]

Discrete sampling of this shifted charge density on the original fixed grid at integer u yields

[\eqalignno { \rho _{\Delta x}(u) =\, & \int_u^{u + 1} {\rm d}x\rho (x - \Delta x) \cr =\, & \int_u^{u + 1} {\rm d}x \sum_{{u^\prime}}\rho (u^\prime) [\theta (x - \Delta x - u^\prime)\cr & - \theta (x - \Delta x - u^\prime - 1)].\cr & & (12)}]

Decomposing the shift [\Delta x = {u_0} + {u_1}] into an integer component u0 plus a positive fractional component 0 ≤ u1 < 1 yields

[\eqalignno{{\rho _{\Delta x}}\left(u \right) =\, & \sum \limits_{{u^\prime}}^{} \rho \left({{u^\prime}} \right) \int \limits_u^{u + 1} dx \big[\theta \left({x - {u_0} - {u^\prime} - {u_1}} \right) \cr & - \theta \left({x - {u_0} - {u^\prime} - {u_1} - 1} \right) \big], &(13)}]

which is only nonzero for the terms [u^{\prime} = u - {u_0}] and [u^{\prime} = u - {u_0} - 1]. Performing the integrals for these terms yields

[{\rho _{\Delta x}}\left(u \right) = \rho \left({u - {u_0}} \right)\left({1 - {u_1}} \right) + \rho \left({u - {u_0} - 1} \right){u_1}. \eqno (14)]

Define the structure factors [A\left(h \right) = DFT\left [{\rho \left(u \right)} \right]] and [{A_{\Delta x}}\left(h \right) = DFT\left [{{\rho _{\Delta x}}\left(u \right)} \right]]. [{A_{\Delta x}}\left(h \right)] is obtained from A(h ) by applying the DFT shift theorem to equation (14)[link]

[{A_{\Delta x}}\left(h \right) = {e^{{{ - 2\pi ih{u_0}} \over N}}}A\left(h \right)\left [{\left({1 - {u_1}} \right) + {e^{{{ - 2\pi ih} \over N}}}{u_1}} \right], \eqno (15)]

which is the one-dimensional version of equation (3)[link]. Equation (15)[link] is exact, and demonstrates that [{A_{\Delta x}}\left(h \right) \ne {e^{{{ - 2\pi ih\Delta x} \over N}}}A\left(h \right)] if [\Delta x] includes a nonzero fractional shift u1.

Supporting information


Acknowledgements

This study was supported by the US Department of Energy under Contract DE-AC52-06NA25396 through the Laboratory-Directed Research and Development Program at Los Alamos National Laboratory.

References

First citationAubert, E., Lebègue, S., Marsman, M., Bui, T. T. T., Jelsch, C., Dahaoui, S., Espinosa, E. & Ángyán, J. G. (2011). J. Phys. Chem. A, 115, 14484–14494.  Web of Science CrossRef CAS PubMed Google Scholar
First citationBąk, J. M., Domagała, S., Hübschle, C., Jelsch, C., Dittrich, B. & Dominiak, P. M. (2011). Acta Cryst. A67, 141–153.  Web of Science CrossRef IUCr Journals Google Scholar
First citationBendeif, E. & Jelsch, C. (2007). Acta Cryst. C63, o361–o364.  Web of Science CSD CrossRef CAS IUCr Journals Google Scholar
First citationBirkedal, H., Madsen, D., Mathiesen, R. H., Knudsen, K., Weber, H.-P., Pattison, P. & Schwarzenbach, D. (2004). Acta Cryst. A60, 371–381.  Web of Science CSD CrossRef CAS IUCr Journals Google Scholar
First citationBowler, D. R. & Miyazaki, T. (2010). J. Phys. Condens. Matter, 22, 074207.  Web of Science CrossRef PubMed Google Scholar
First citationBowler, D. R. & Miyazaki, T. (2012). Rep. Prog. Phys. 75, 036503.  Web of Science CrossRef PubMed Google Scholar
First citationBruning, H. & Feil, D. (1992). Acta Cryst. A48, 865–872.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationCapelli, S. C., Bürgi, H.-B., Dittrich, B., Grabowsky, S. & Jayatilaka, D. (2014). IUCrJ, 1, 361–379.  Web of Science CSD CrossRef CAS PubMed IUCr Journals Google Scholar
First citationCollaborative Computational Project, Number 4 (1994). Acta Cryst. D50, 760–763.  CrossRef IUCr Journals Google Scholar
First citationCoppens, P., Pautler, D. & Griffin, J. F. (1971). J. Am. Chem. Soc. 93, 1051–1058.  CrossRef CAS Web of Science Google Scholar
First citationCraven, B. M. & Weber, H.-P. (1977). Technical Report, Department of Crystallography, University of Pittsburgh.  Google Scholar
First citationDadda, N., Nassour, A., Guillot, B., Benali-Cherif, N. & Jelsch, C. (2012). Acta Cryst. A68, 452–463.  Web of Science CSD CrossRef CAS IUCr Journals Google Scholar
First citationDawson, B. (1967a). Proc. R. Soc. A: Math. Phys. Eng. Sci. 298, 255–263.  Google Scholar
First citationDawson, B. (1967b). Proc. R. Soc. A: Math. Phys. Eng. Sci. 298, 264–288.  Google Scholar
First citationDittrich, B., Hübschle, C. B., Messerschmidt, M., Kalinowski, R., Girnt, D. & Luger, P. (2005). Acta Cryst. A61, 314–320.  Web of Science CSD CrossRef CAS IUCr Journals Google Scholar
First citationDittrich, B., Weber, M., Kalinowski, R., Grabowsky, S., Hübschle, C. B. & Luger, P. (2009). Acta Cryst. B65, 749–756.  Web of Science CSD CrossRef CAS IUCr Journals Google Scholar
First citationEmsley, P., Lohkamp, B., Scott, W. G. & Cowtan, K. (2010). Acta Cryst. D66, 486–501.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationGoedecker, S. (1999). Rev. Mod. Phys. 71, 1085–1123.  Web of Science CrossRef CAS Google Scholar
First citationGrosse-Kunstleve, R. W., Sauter, N. K., Moriarty, N. W. & Adams, P. D. (2002). J. Appl. Cryst. 35, 126–136.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationGuillot, B., Viry, L., Guillot, R., Lecomte, C. & Jelsch, C. (2001). J. Appl. Cryst. 34, 214–223.  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 citationHirshfeld, F. L. (1971). Acta Cryst. B27, 769–781.  CSD CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationHirshfeld, F. L. (1977a). Isr. J. Chem. 16, 226–229.  CrossRef CAS Google Scholar
First citationHirshfeld, F. L. (1977b). Theor. Chim. Acta, 44, 129–138.  CrossRef CAS Web of Science Google Scholar
First citationJayatilaka, D. (1998). Phys. Rev. Lett. 80, 798–801.  Web of Science CrossRef CAS Google Scholar
First citationJayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383–393.  Web of Science CrossRef CAS IUCr Journals 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 citationKresse, G. & Furthmüller, J. (1996). Phys. Rev. B, 54, 11169–11186.  CrossRef CAS Web of Science Google Scholar
First citationKresse, G. & Joubert, D. (1999). Phys. Rev. B, 59, 1758–1775.  Web of Science CrossRef CAS Google Scholar
First citationLipscomb, W. (1972). Trans. Am. Crystallogr. Assoc. 8, 79–92.  CAS Google Scholar
First citationMacrae, C. F., Bruno, I. J., Chisholm, J. A., Edgington, P. R., McCabe, P., Pidcock, E., Rodriguez-Monge, L., Taylor, R., van de Streek, J. & Wood, P. A. (2008). J. Appl. Cryst. 41, 466–470.  Web of Science CSD CrossRef CAS IUCr Journals Google Scholar
First citationMadsen, A. Ø., Civalleri, B., Ferrabone, M., Pascale, F. & Erba, A. (2013). Acta Cryst. A69, 309–321.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationMassa, L., Huang, L. & Karle, J. (1995). Int. J. Quantum Chem. 56, 371–384.  CrossRef Google Scholar
First citationMniszewski, S. M., Cawkwell, M. J., Wall, M. E., Mohd-Yusof, J., Bock, N., Germann, T. C. & Niklasson, A. M. N. (2015). J. Chem. Theory Comput. 11, 4644–4654.  Web of Science CrossRef CAS PubMed Google Scholar
First citationPress, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (1999). Numerical Recipes in C, 2nd ed. Cambridge University Press.  Google Scholar
First citationSheldrick, G. M. (2008). Acta Cryst. A64, 112–122.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationSheldrick, G. M. (2015). Acta Cryst. C71, 3–8.  Web of Science CrossRef IUCr Journals Google Scholar
First citationStein, P. E., Boodhoo, A., Armstrong, G. D., Cockle, S. A., Klein, M. H. & Read, R. J. (1994). Structure, 2, 45–57.  CrossRef CAS PubMed Web of Science Google Scholar
First citationStewart, R. F. (1969). J. Chem. Phys. 51, 4569–4577.  CrossRef CAS Web of Science Google Scholar
First citationStewart, R. F. (1970). J. Chem. Phys. 53, 205–213.  CrossRef CAS Web of Science Google Scholar
First citationStewart, R. F. & Spackman, M. A. (1983). Valray User's Manual. Carnegie Mellon University, USA, and University of Copenhagen, Denmark.  Google Scholar
First citationSwaminathan, S., Craven, B. M. & McMullan, R. K. (1984). Acta Cryst. B40, 300–306.  CSD CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationTang, W., Sanville, E. & Henkelman, G. (2009). J. Phys. Condens. Matter, 21, 084204.  Web of Science CrossRef PubMed Google Scholar
First citationVandeVondele, J., Borštnik, U. & Hutter, J. (2012). J. Chem. Theory Comput. 8, 3565–3573.  Web of Science CrossRef CAS PubMed Google Scholar
First citationVolkov, A., Macchi, P., Farrugia, L. J., Gatti, C., Mallinson, P., Richter, T. & Koritsanszky, T. (2006). XD2006 – A Computer Program Package for Multipole Refinement, Topological Analysis of Charge Densities and Evaluation of Intermolecular Energies from Experimental and Theoretical Structure Factors. University of New York at Buffalo, USA.  Google Scholar
First citationVries, R. Y. de, Feil, D. & Tsirelson, V. G. (2000). Acta Cryst. B56, 118–123.  Web of Science CrossRef IUCr Journals Google Scholar
First citationWall, M. E. (2009). Methods Mol. Biol. 544, 269–279.  CrossRef PubMed CAS Google Scholar
First citationWall, M. E., Clarage, J. B. & Phillips, G. N. (1997a). Structure, 5, 1599–1612.  Web of Science CrossRef CAS PubMed Google Scholar
First citationWall, M. E., Ealick, S. E. & Gruner, S. M. (1997b). Proc. Natl Acad. Sci. USA, 94, 6180–6184.  CrossRef CAS PubMed Web of Science Google Scholar
First citationWall, M. E., Van Benschoten, A. H., Sauter, N. K., Adams, P. D., Fraser, J. S. & Terwilliger, T. C. (2014). Proc. Natl Acad. Sci. USA, 111, 17887–17892.  Web of Science CrossRef CAS PubMed Google Scholar
First citationWhitten, A. E. & Spackman, M. A. (2006). Acta Cryst. B62, 875–888.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationWoińska, M., Jayatilaka, D., Spackman, M. A., Edwards, A. J., Dominiak, P. M., Woźniak, K., Nishibori, E., Sugimoto, K. & Grabowsky, S. (2014). Acta Cryst. A70, 483–498.  Web of Science CSD CrossRef 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.

IUCrJ
Volume 3| Part 4| July 2016| Pages 237-246
ISSN: 2052-2525