research papers
Quantum crystallographic charge density of urea
^{a}Computer, Computational, and Statistical Sciences Division, Los Alamos National Laboratory, Mail Stop B256, Los Alamos, New Mexico 87545, USA
^{*}Correspondence email: mewall@lanl.gov
Standard Xray crystallography methods use freeatom models to calculate mean unitcell charge densities. Real molecules, however, have shared charge that is not captured accurately using freeatom models. To address this limitation, a charge density model of crystalline urea was calculated using highlevel quantum theory and was refined against publicly available ultrahighresolution 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 Xray crystallographic atomic displacement parameters. The results demonstrate the feasibility and benefits of integrating fully periodic quantum charge density calculations into ultrahighresolution Xray crystallographic model building and refinement.
Keywords: charge density; quantum theory; spherical atom model.
1. Introduction
Efforts to increase the accuracy of charge density models from Xray crystallography have mainly focused on fitting the Bragg data using functions that are more expressive than the usual freeatom spherical distributions. Stewart (1969) proposed using general scattering factors that are the products of atomcentered orbital wavefunctions, and restrictions to better match the number of free parameters to the number of reflections in fitting (Stewart, 1970). Coppens et al. (1971) 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), and expanded each atomcentered charge density in spherical harmonics (Dawson, 1967b). Hirshfeld developed a leastsquares method that models aspherical atomic charge densities using basis functions related to spherical harmonics, but with alternative symmetry properties (Hirshfeld, 1971). Spherical harmonicrelated methods were integrated into multipole computer programs that are used when charge density models are desired (Hansen & Coppens, 1978; Hirshfeld, 1977a; Craven & Weber, 1977; Stewart & Spackman, 1983; Jelsch et al., 2005; Volkov et al., 2006).
Although less well exploited than multipole methods, the potential for combining quantum theory and Xray diffraction to obtain accurate charge density models of molecular crystals has been long appreciated (Lipscomb, 1972). This combination has been termed quantum crystallography (Massa et al., 1995). 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, 2012; Goedecker, 1999; VandeVondele et al., 2012), and fast quantum simulations for systems approaching 10^{4} atoms with 10^{5} time steps are now possible (Mniszewski et al., 2015). 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 smallmolecule crystallography (Capelli et al., 2014) but also for macromolecular crystallography.
Several methods have been proposed for quantum crystallography, including the method of kernel projector matrices (Massa et al., 1995) and fitting of wavefunctions to diffraction data (Jayatilaka, 1998). One method that is showing promise in practical applications is Hirshfeld Atom (HAR) (Bruning & Feil, 1992; Capelli et al., 2014; Jayatilaka & Dittrich, 2008). 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). The partitioned charge is used to calculate aspherical atomic structure factors that are substituted for the usual structure factors in crystallographic considering both the atomic positions and displacement parameters (Bruning & Feil, 1992). Whereas Bruning & Feil (1992) 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; Jayatilaka & Dittrich, 2008) 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 of atomic positions (Capelli et al., 2014). So far HAR has been limited to gasphase electronic structure calculations, with cluster charges placed at symmetryrelated 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 Xray 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 Xray diffraction from crystalline benzene and urea (Jayatilaka & Dittrich, 2008), Lphenylalaninium hydrogen maleate (Woińska et al., 2014), and a Gly LAla dipeptide (Capelli et al., 2014) found that HAR bond distances agreed very well with neutron crystal structures, overcoming known deficiencies in sphericalatom charge density models (Lipscomb, 1972). 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) is adapted for crystalline phase electronic structure calculations performed using VASP (Kresse & Furthmüller, 1996). Electronic structure calculations using VASP previously were performed on hexachlorobenzene for comparisons to the Xray crystallographic multipole charge density (Aubert et al., 2011), but without allowing for individual ADPs. The novel aspect of the present method therefore is the combination of a crystalline phase densityfunctionaltheorybased 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 but also both 2F_{o} − F_{c} 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
Ultrahighresolution urea synchrotron diffraction data were obtained from Birkedal et al. (2004) 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 (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), Table 1. The multipole refined urea was obtained from Birkedal et al. (2004) 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 (Swaminathan et al., 1984).

2.2. Spherical atom and multipole crystallographic models
The program SHELXL (Sheldrick, 2008, 2015), 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: R_{1} = 0.0370, wR_{2} = 0.0796, and SHELX goodness of fit = 0.639 for all reflections. Mean chargedensity maps F_{o}, 2F_{o} − F_{c}, and F_{o} − F_{c} 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 The maps were expanded to P1 using CCP4 (Stein et al., 1994) mapmask and were interpolated to a 64 × 64 × 64 grid using CCP4 maprot.
The program MoPro, Version 14.06 (Guillot et al., 2001; Jelsch et al., 2005), was used to refine a multipole chargedensity model of urea. was based on the ultrahighresolution data and structure from Birkedal et al. (2004). 20 cycles of automated density 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 P_{lm} (PLM) multipole parameters. The atom coordinates and ADPs were kept constant. MoPro reported agreement factors for the final model are: R_{F} = 0.0242, wR^{2}F = 0.0107, R_{I} = 0.0212, wR^{2}I = 0.0212, and GooF = 2.293 for 992 nonzero reflections. Charge density maps were calculated using the MoPro supplied program VMoPro. The F_{o}, 2F_{o} − F_{c}, and F_{o} − F_{c} 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). The total static crystalline charge density was computed using the VMoPro STAT method, with a 10 Å selection for grid limits, gridcube 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 ) 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 of Swaminathan et al. (1984); and the multipole model of Birkedal et al. (2004). An expanded with 16 atoms was generated using the Computational Crystallography Toolbox (cctbx) (GrosseKunstleve et al., 2002) by applying symmetry to the fiveatom Ab initio density functional theory calculations were performed using VASP, Version 5.3.3 (Kresse & Furthmüller, 1996). Instead of pseudopotentials, the PAW method was used, with PAW_PBE parameters (Kresse & Joubert, 1999). The electronic structure was computed using 4^{3} = 64 MonkhorstPack 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 LREAL = .FALSE. was used to evaluate projection operators in as recommended in the VASP documentation. The valence charge density v(x) was calculated for the expanded P1 In addition, using the same VASP PAW method as for the molecular calculation, 16 crystalline core charge densities c_{i}(x) and 16 crystalline freeatom (`promolecule') charge densities f_{i}(x ) were obtained for each individual atom i.
method (Jayatilaka & Dittrich, 2008To achieve the desired model accuracy, all VASP charge densities were calculated on a 128 × 128 × 128 grid spanning the For crystallographic 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.
Xray structure factors were calculated using both new and existing tools in Lunus software (Wall, 2009), which was originally designed for analysis and modeling of diffuse Xray scattering data (Wall et al., 1997a,b, 2014). The effect of ADPs was modeled using a Stockholder method (Bruning & Feil, 1992; Hirshfeld, 1977b). The total valence density v(x) was partitioned into atomic contributions using the equation
Hirshfeld partitioning (Hirshfeld, 1977b) was used with weights q_{i}(x ) defined using the freeatom charge density f_{i}(x )
Similar to Bruning & Feil (1992), ADPs were modeled by treating each partitioned atom charge density a_{i}(x ) = c_{i}(x ) + v_{i}(x ) as a rigid distribution, displaced along with the atom. However, in contrast to Bruning & Feil (1992), instead of using a multipole expansion, the charge density a_{i}(uvw ) was sampled on a rectilinear grid spanning the indexed by uvw. This method is similar to that of Jayatilaka & Dittrich (2008), who used a radialangular 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 A_{i}(hkl ) = DFT[a_{i}(uvw ) ], where DFT denotes a discrete Fourier transform. The DFT was computed using a fast Fourier transform (FFT) algorithm (Press et al., 1999).
In the original Hirshfeld ), the value of A_{i}(hkl ) after a coordinate shift was obtained by multiplying A_{i}(hkl ) by a phase factor. Although multiplication by a phase factor is appropriate for arbitrary translations of a continuous distribution or an atomcentered 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 a_{i}(uvw ) on the original grid (Appendix). The structure factors are obtained by transforming to the grid coordinates , decomposing these coordinates into integer (u_{0}v_{0}w_{0}) and fractional (u_{1}v_{1}w_{1}) parts such that , , and , and using the following equation to calculate
method (Jayatilaka & Dittrich, 2008The unitcell F_{c}(hkl ) was then calculated as
where is the Debye–Waller factor for the matrix U_{i} of ADPs for atom i, and s_{hkl} is the scattering vector corresponding to hkl.
2.4. Quantum model refinement
Quantum refinements were performed starting with the spherical atom (S), multipole (M) (Birkedal et al., 2004), and neutron crystallography (N) (Swaminathan et al., 1984) atomic coordinates and ADPs. Model was performed by minimizing the goodnessoffit (GooF) statistic
where I_{c}(hkl ) =  F_{c}(hkl ) ^{2}, I_{o}(hkl ) and are the values and errors of the observed intensities, and the number of NDF = 965 is the number of data points (= 992 nonnegative 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 U_{i} 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 U^{23} = 0 for C, O atoms; Y = X+0.5 for all atoms; and U^{22} = U^{11} and U^{13} = U^{23} 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 including the scale factor between the data and model (the same number as for spherical atom but without geometry restraints).
The mean ρ were calculated using Fourier reconstruction as . The experimental F_{o}, 2F_{o} − F_{c} and F_{o} − F_{c} maps were calculated by applying the model phases to the observations. Values of  F_{c}(hkl )  were used in place of missing values of  F_{o}(hkl ) . Only complete grids were used for FFT calculations on the quantum models; reflections were not truncated using a resolution cutoff. The F_{o}, F_{c}, 2F_{o} − F_{c}, and F_{o} − F_{c} maps obtained using the multipole model as an input structure are provided in the supporting information.
charge densities2.5. Agreement factors
The agreement of all models with the diffraction data was assessed using several standard statistics: GooF, wR^{2}F, wR^{2}I, RF, and RSR. The GooF [equation (5)] was used as the target. The weighted Rsquared factor for amplitudes, wR^{2}F, was calculated as
where  F_{o}(hkl )  and are the experimental reflection amplitudes and errors, and F_{c}(hkl ) is calculated using equation (4). The weighted Rsquared factor for intensities, wR^{2}I, was calculated as
the R factor for amplitudes, RF, was calculated as
the realspace Rfactor, RSR, was calculated as
where is the experimental F_{o} map. Calculated and observed values were scaled to minimize the RMSD prior to using equations (5)–(9), and both and were offset to have zero mean prior to using equation (9). 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): wR^{2}F (target) = 0.9, wR^{2}I = 1.8, RF = 2.3, RSR = 2.1 and GooF = 1.7 (Table 1). 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).
Threedimensional visualizations of the 2F_{o} − F_{c} and F_{o} − F_{c} maps for the spherical atom, quantumM, and multipole models are shown in Fig. 1. The spherical atom and multipole 2F_{o} − F_{c} 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 2F_{o} − F_{c} realspace maps, using an appropriately modified equation (9). 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).
Visualization of contours in a twodimensional section including the C=O bond reveals that the quantum and multipole 2F_{o} − F_{c} maps are very different (Figs. 2a and b). (Much of this difference might be an artifact in the multipole map calculation, as mentioned below.) Compared with the multipole model (Fig. 2b), the quantumM model is smoother (Fig. 2a). The multipole model shows ripples surrounding core atoms and peaks away from atoms, including between bonded heavy atoms. The quantumM model has some peaks away from the atom cores (Fig. 2a), but these are lower in magnitude compared with the multipole model. The quantumM and multipole model F_{o} − F_{c} difference maps are broadly similar (Figs. 2c 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).
To investigate further the differences between the 2F_{o} − F_{c} maps of the quantumM 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). There are visible differences (Figs. 3a and b); however, the differences are much smaller than in the 2F_{o} − F_{c} maps (Figs. 2a and b), and they coincide with atoms and bonds. The comparison suggests that the ripples in the 2F_{o} − F_{c} 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).
Subtracting the static charge densities from VASP and MoPro reveals substantial differences in the charge distribution along the C=O bond (Fig. 3c). These differences show a similar pattern of peaks and troughs as in the F_{o} − F_{c} map for the multipole model (Fig. 2d); by comparison, the F_{o} − F_{c} map of the quantumM model shows smaller differences along the C=O bond (Fig. 2c). Combined, Figs. 2 and 3 indicate that the multipole static charge density contains deviations from the data in the C=O bond that are decreased in the quantumM model.
The VASP and MoPro calculations were further compared using a Bader analysis of the net atom charges (Tang et al., 2009; Table 2). 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.

The atomic coordinates of the quantumM, quantumN, neutron and multipole models are all very similar (Table 3). The differences between these models and either the spherical atom or quantumS 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 quantumS structure compared with the neutron structure (Table 4): the bond length is 1.006 Å in the neutron structure compared with 0.911 Å in the spherical atom and 0.810 Å in the quantumS structure. The differences also lead to decreases in the C—N—H1 and C—N—H2 bond angles for both the spherical atom and quantumS structures compared with the neutron structure (Table 5): the angles in the spherical atom model are about 2° smaller, and the angles in the quantumS 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 quantumS structure. The angle deviations for the quantumS model are visible in the stick diagram in Fig. 3(a); the effect is smaller but still perceptible for the spherical atom model (not shown).



The similarity of ADPs was assessed using the S statistic, which describes the deviation of the threedimensional positional distribution of the atoms defined by the ADPs (Whitten & Spackman, 2006). The ADPs for the heavy atoms in the quantum models are very similar to the neutron models (Table 6; Fig. 4): 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 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 highlevel quantum theoretical calculation of vibrations of urea to obtain ADPs (Madsen et al., 2013) 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) 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/631G(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).
To assess the convergence of the quantum et al., 2014; Jayatilaka & Dittrich, 2008), the electronic structure calculation was iteratively applied to each of the quantum models. In the iteration, the refined atomic coordinates were used to recompute all charge densities using VASP, and the model was rerefined against the data using the new densities. The quantumM and quantumN 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 quantumM model in Fig. 3(d); the agreement factors remained the same as in Table 1; 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 quantumS 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 yposition 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).
as was done in previous HAR implementations (Capelli4. Discussion
The agreement of the quantum crystallographic models of urea with ultrahighresolution data compares favorably to the multipole model. Both the 2F_{o} − F_{c} map and the total static charge density are substantially different between the quantum and multipole models, however. The differences in 2F_{o} − F_{c} 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. 2a 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. 3c) and in the 0.5electron higher negative charge associated with the C atom for the multipole model (Table 2). The difference in the C atom charge is consistent with the multipole charge density study of Birkedal et al. (2004), 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) 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) which found that the charge density of urea could not be determined uniquely using multipole 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).
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
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) and qualitatively (Fig. 4) similar to the neutron This was even the case for the quantumS model, despite the lack of convergence seen in a second iteration of and deviations in the geometry with respect to the neutron model (Tables 4 and 5). This finding is consistent with studies in which the neutron crystallographic temperature factors of urea (Jayatilaka & Dittrich, 2008) and other systems (Capelli et al., 2014; Woińska et al., 2014) were found to be reproduced reasonably well using HAR. In particular, the previous urea study (Jayatilaka & Dittrich, 2008) used BLYP density functional theory with a ccpVTZ basis and surrounding charge clusters to mimic periodic boundary conditions, and used the same starting structure as the present quantumM (Birkedal et al., 2004). The following ADP values were obtained for H atoms (Jayatilaka & Dittrich, 2008) (U^{11}, U^{33}, U^{12}, and U^{13} 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); 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 quantumM 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), 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 quantumM and quantumN models indicate the potential advantages of quantum crystallography for accurate charge density studies. The lack of convergence and geometry deviations of the quantumS model, however, indicate that challenges remain for the general applicability of these methods. The deviations of the quantumS 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), leading to corresponding deviations in geometry (Tables 4 and 5). These deviations are increased rather than decreased in the quantumS model, which prevents the quantum 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 Xray crystallographic LAla using spherical atoms models like those used here as input structures (Capelli et al., 2014) 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) converges using the present spherical atom model of urea as an input (as mentioned above, the study of Jayatilaka & Dittrich, 2008, made use of the same model as the quantumM 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; Bendeif & Jelsch, 2007; Dadda et al., 2012; Dittrich et al., 2005, 2009), which can place H atoms to within O(10^{−2}) Å of the positions in neutron crystal structures.
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–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) to quantum crystallography. Increasing the speed of calculations would enable the wider adaptation of HAR by adapting existing smallmolecule crystallography workflows, and would provide a complementary quantum crystallographic alternative to multipole (Jelsch et al., 2000) for obtaining highresolution 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 onedimensional case; the extension to three dimensions in equation (3) is straightforward. The structure factors are given by the discrete Fourier transform (DFT) of a periodic charge density sampled at N fixed grid points u, with . Let the charge density correspond to the following continuous stepwise distribution
where the Heaviside distribution for and otherwise. Translation of by a shift yields a new charge density
Discrete sampling of this shifted charge density on the original fixed grid at integer u yields
Decomposing the shift into an integer component u_{0} plus a positive fractional component 0 ≤ u_{1} < 1 yields
which is only nonzero for the terms and . Performing the integrals for these terms yields
Define the structure factors and . is obtained from A(h ) by applying the DFT shift theorem to equation (14)
which is the onedimensional version of equation (3). Equation (15) is exact, and demonstrates that if includes a nonzero fractional shift u_{1}.
Supporting information
Compressed archive containing charge density maps in CCP4 format. DOI: 10.1107/S2052252516006242/fc5014sup1.gz
Acknowledgements
This study was supported by the US Department of Energy under Contract DEAC5206NA25396 through the LaboratoryDirected Research and Development Program at Los Alamos National Laboratory.
References
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. Web of Science CrossRef CAS PubMed Google Scholar
Bą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
Bendeif, E. & Jelsch, C. (2007). Acta Cryst. C63, o361–o364. Web of Science CSD CrossRef CAS IUCr Journals Google Scholar
Birkedal, 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
Bowler, D. R. & Miyazaki, T. (2010). J. Phys. Condens. Matter, 22, 074207. Web of Science CrossRef PubMed Google Scholar
Bowler, D. R. & Miyazaki, T. (2012). Rep. Prog. Phys. 75, 036503. Web of Science CrossRef PubMed Google Scholar
Bruning, H. & Feil, D. (1992). Acta Cryst. A48, 865–872. CrossRef CAS Web of Science IUCr Journals Google Scholar
Capelli, 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
Collaborative Computational Project, Number 4 (1994). Acta Cryst. D50, 760–763. CrossRef IUCr Journals Google Scholar
Coppens, P., Pautler, D. & Griffin, J. F. (1971). J. Am. Chem. Soc. 93, 1051–1058. CrossRef CAS Web of Science Google Scholar
Craven, B. M. & Weber, H.P. (1977). Technical Report, Department of Crystallography, University of Pittsburgh. Google Scholar
Dadda, N., Nassour, A., Guillot, B., BenaliCherif, N. & Jelsch, C. (2012). Acta Cryst. A68, 452–463. Web of Science CSD CrossRef CAS IUCr Journals Google Scholar
Dawson, B. (1967a). Proc. R. Soc. A: Math. Phys. Eng. Sci. 298, 255–263. Google Scholar
Dawson, B. (1967b). Proc. R. Soc. A: Math. Phys. Eng. Sci. 298, 264–288. Google Scholar
Dittrich, 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
Dittrich, 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
Emsley, P., Lohkamp, B., Scott, W. G. & Cowtan, K. (2010). Acta Cryst. D66, 486–501. Web of Science CrossRef CAS IUCr Journals Google Scholar
Goedecker, S. (1999). Rev. Mod. Phys. 71, 1085–1123. Web of Science CrossRef CAS Google Scholar
GrosseKunstleve, 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
Guillot, 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
Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909–921. CrossRef CAS IUCr Journals Web of Science Google Scholar
Hirshfeld, F. L. (1971). Acta Cryst. B27, 769–781. CSD CrossRef CAS IUCr Journals Web of Science Google Scholar
Hirshfeld, F. L. (1977a). Isr. J. Chem. 16, 226–229. CrossRef CAS Google Scholar
Hirshfeld, F. L. (1977b). Theor. Chim. Acta, 44, 129–138. CrossRef CAS Web of Science Google Scholar
Jayatilaka, D. (1998). Phys. Rev. Lett. 80, 798–801. Web of Science CrossRef CAS Google Scholar
Jayatilaka, D. & Dittrich, B. (2008). Acta Cryst. A64, 383–393. Web of Science CrossRef CAS IUCr Journals Google Scholar
Jelsch, C., Guillot, B., Lagoutte, A. & Lecomte, C. (2005). J. Appl. Cryst. 38, 38–54. Web of Science CrossRef IUCr Journals Google Scholar
Jelsch, C., Teeter, M. M., Lamzin, V., PichonPesme, V., Blessing, R. H. & Lecomte, C. (2000). Proc. Natl Acad. Sci. USA, 97, 3171–3176. Web of Science CrossRef PubMed CAS Google Scholar
Kresse, G. & Furthmüller, J. (1996). Phys. Rev. B, 54, 11169–11186. CrossRef CAS Web of Science Google Scholar
Kresse, G. & Joubert, D. (1999). Phys. Rev. B, 59, 1758–1775. Web of Science CrossRef CAS Google Scholar
Lipscomb, W. (1972). Trans. Am. Crystallogr. Assoc. 8, 79–92. CAS Google Scholar
Macrae, C. F., Bruno, I. J., Chisholm, J. A., Edgington, P. R., McCabe, P., Pidcock, E., RodriguezMonge, 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
Madsen, A. Ø., Civalleri, B., Ferrabone, M., Pascale, F. & Erba, A. (2013). Acta Cryst. A69, 309–321. Web of Science CrossRef CAS IUCr Journals Google Scholar
Massa, L., Huang, L. & Karle, J. (1995). Int. J. Quantum Chem. 56, 371–384. CrossRef Google Scholar
Mniszewski, S. M., Cawkwell, M. J., Wall, M. E., MohdYusof, 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
Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (1999). Numerical Recipes in C, 2nd ed. Cambridge University Press. Google Scholar
Sheldrick, G. M. (2008). Acta Cryst. A64, 112–122. Web of Science CrossRef CAS IUCr Journals Google Scholar
Sheldrick, G. M. (2015). Acta Cryst. C71, 3–8. Web of Science CrossRef IUCr Journals Google Scholar
Stein, 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
Stewart, R. F. (1969). J. Chem. Phys. 51, 4569–4577. CrossRef CAS Web of Science Google Scholar
Stewart, R. F. (1970). J. Chem. Phys. 53, 205–213. CrossRef CAS Web of Science Google Scholar
Stewart, R. F. & Spackman, M. A. (1983). Valray User's Manual. Carnegie Mellon University, USA, and University of Copenhagen, Denmark. Google Scholar
Swaminathan, S., Craven, B. M. & McMullan, R. K. (1984). Acta Cryst. B40, 300–306. CSD CrossRef CAS Web of Science IUCr Journals Google Scholar
Tang, W., Sanville, E. & Henkelman, G. (2009). J. Phys. Condens. Matter, 21, 084204. Web of Science CrossRef PubMed Google Scholar
VandeVondele, J., Borštnik, U. & Hutter, J. (2012). J. Chem. Theory Comput. 8, 3565–3573. Web of Science CrossRef CAS PubMed Google Scholar
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. Google Scholar
Vries, R. Y. de, Feil, D. & Tsirelson, V. G. (2000). Acta Cryst. B56, 118–123. Web of Science CrossRef IUCr Journals Google Scholar
Wall, M. E. (2009). Methods Mol. Biol. 544, 269–279. CrossRef PubMed CAS Google Scholar
Wall, M. E., Clarage, J. B. & Phillips, G. N. (1997a). Structure, 5, 1599–1612. Web of Science CrossRef CAS PubMed Google Scholar
Wall, 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
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. Web of Science CrossRef CAS PubMed Google Scholar
Whitten, A. E. & Spackman, M. A. (2006). Acta Cryst. B62, 875–888. Web of Science CrossRef CAS IUCr Journals Google Scholar
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. Web of Science CSD CrossRef IUCr Journals Google Scholar
This is an openaccess article distributed under the terms of the Creative Commons Attribution (CCBY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.