Quantum crystallographic charge density of urea
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.
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) 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). 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 atom-centered charge density in spherical harmonics (Dawson, 1967b). 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). Spherical harmonic-related methods were integrated into multipole refinement 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 X-ray 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 molecular dynamics simulations for systems approaching 104 atoms with 105 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 small-molecule 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 Refinement (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 refinement, 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 refinement of atomic positions (Capelli et al., 2014). 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), L-phenylalaninium hydrogen maleate (Woińska et al., 2014), and a Gly L-Ala dipeptide (Capelli et al., 2014) found that HAR bond distances agreed very well with neutron crystal structures, overcoming known deficiencies in spherical-atom 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 X-ray 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 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 2Fo − Fc 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.
Ultra-high-resolution 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 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), Table 1. The multipole refined urea crystal structure 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 crystal structure (Swaminathan et al., 1984).
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: R1 = 0.0370, wR2 = 0.0796, and SHELX goodness of fit = 0.639 for all reflections. Mean unit cell charge-density maps Fo, 2Fo − Fc, and Fo − Fc 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) 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 charge-density model of urea. Refinement was based on the ultra-high-resolution data and structure from Birkedal et al. (2004). 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, 2Fo − Fc, and Fo − Fc 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, 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.
A custom implementation of the original Hirshfeld atom refinement method (Jayatilaka & Dittrich, 2008) 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); and the multipole model of Birkedal et al. (2004). An expanded unit cell with 16 atoms was generated using the Computational Crystallography Toolbox (cctbx) (Grosse-Kunstleve et al., 2002) by applying 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). Instead of pseudopotentials, the PAW method was used, with PAW_PBE parameters (Kresse & Joubert, 1999). 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), which was originally designed for analysis and modeling of diffuse X-ray 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
Similar to Bruning & Feil (1992), 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), 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), 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).
In the original Hirshfeld refinement method (Jayatilaka & Dittrich, 2008), the value of Ai(hkl ) after a coordinate shift 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). The structure factors are obtained by transforming to the grid coordinates , decomposing these coordinates into integer (u0v0w0) and fractional (u1v1w1) parts such that , , and , and using the following equation to calculate
The unit-cell structure factor Fc(hkl ) was then calculated as
where 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.
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 refinement was performed by minimizing the goodness-of-fit (GooF) statistic
where Ic(hkl ) = | Fc(hkl ) |2, Io(hkl ) and 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 . The experimental Fo, 2Fo − Fc and Fo − Fc 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, 2Fo − Fc, and Fo − Fc maps obtained using the multipole model as an input structure are provided in the supporting information.
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)] was used as the refinement target. The weighted R-squared factor for amplitudes, wR2F, was calculated as
the R factor for amplitudes, RF, was calculated as
the real-space R-factor, RSR, was calculated as
where is the experimental Fo 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.
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). 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).
Three-dimensional visualizations of the 2Fo − Fc and Fo − Fc maps for the spherical atom, quantum-M, and multipole models are shown in Fig. 1. The spherical atom and multipole 2Fo − Fc 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 2Fo − Fc real-space 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 two-dimensional section including the C=O bond reveals that the quantum and multipole 2Fo − Fc 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 quantum-M model is smoother (Fig. 2a). 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. 2a), but these are lower in magnitude compared with the multipole model. The quantum-M and multipole model Fo − Fc 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 2Fo − Fc 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). There are visible differences (Figs. 3a and b); however, the differences are much smaller than in the 2Fo − Fc maps (Figs. 2a and b), and they coincide with atoms and bonds. The comparison suggests that the ripples in the 2Fo − Fc 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 Fo − Fc map for the multipole model (Fig. 2d); by comparison, the Fo − Fc map of the quantum-M 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 quantum-M 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 quantum-M, quantum-N, neutron and multipole models are all very similar (Table 3). 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): 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): 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(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 three-dimensional 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 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) 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/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).
To assess the convergence of the quantum refinement, as was done in previous HAR implementations (Capelli 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 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(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 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).
The agreement of the quantum crystallographic models of urea with ultra-high-resolution data compares favorably to the multipole model. Both the 2Fo − Fc map and the total static charge density are substantially different between the quantum and multipole models, however. The differences in 2Fo − Fc 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.5-electron 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 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).
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) and qualitatively (Fig. 4) 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 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 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). The following ADP values were obtained for H atoms (Jayatilaka & Dittrich, 2008) (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); 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), 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), leading to corresponding deviations in geometry (Tables 4 and 5). 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) 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 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; 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.
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 small-molecule crystallography workflows, and would provide a complementary quantum crystallographic alternative to multipole refinement (Jelsch et al., 2000) for obtaining high-resolution charge density models of molecular crystals, including macromolecular crystals.
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) 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 step-wise 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 u0 plus a positive fractional component 0 ≤ u1 < 1 yields
which is only nonzero for the terms and . Performing the integrals for these terms yields
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.
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., Benali-Cherif, 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
Grosse-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
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., 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
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., 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
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., 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
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 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.