In X-ray constrained wavefunction (XCW) fitting, external information, such as electron correlation and polarization, is included into a single-determinantal isolated-molecule wavefunction. In a first step, we show that the extraction of these two physical effects by XCW fitting is complete and accurate by comparing to theoretical reference calculations. In a second step, we show that fitting to data from single-crystal x-ray diffraction measurements provides the same results qualitatively and how the physical effects can be separated, although always inherently convolved in the experiment. We further demonstrate that exchange–correlation potentials are systematically affected by XCW fitting in a physically meaningful way, which could be exploited for method development in quantum chemistry, subject to some remaining challenges that we also outline.

## I. INTRODUCTION

In theoretical chemistry, the idea of fitting to electron densities calculated at a higher or different level of theory has been used in various ways, e.g., to speed-up computations by approximating integrals leading to auxiliary basis sets.^{1,2} For the development of functionals in density functional theory (DFT), one of the proposed techniques is a “constrained-search method to determine electronic wave functions from electronic densities,”^{3,4} which has been successfully used to calculate molecular exchange–correlation potentials.^{5} This was the starting point for the development of the TH^{6–8} and HCTH^{9–11} functionals by Tozer and co-workers. In this paper, we ask whether this is also feasible, useful, and meaningful if the densities that we fit come from experimental x-ray diffraction studies. Experimental data would obviously be the ultimate reference since they contain all physical effects at once. However, they are also naturally error-prone due to the presence of experimental uncertainties, which could bias the results.

One of the aspects of quantum crystallography is the enhancement of the information contents of approximate model wavefunctions by fitting to experimental x-ray diffraction structure factors, the Fourier transforms of electron densities.^{12–15} The method of x-ray wavefunction fitting was pioneered by Clinton and Massa^{16} and developed into its most widespread and advanced form by Jayatilaka [X-ray constrained wavefunction (XCW) fitting].^{17–19} The XCW method consists in the minimization of the functional

where *E* is a quantum mechanical expression for the energy [usually at Hartree–Fock (HF) or DFT level] of an isolated molecule and GooF^{2} is an agreement statistic between calculated ($Fmodel$) and measured ($Fexp$) x-ray structure factor amplitudes, which also takes into account standard uncertainties (σ) of the data,

In the previous equation, $h\u20d7$ is the reciprocal lattice vector, *t* is a scale factor, *N*_{r} is the number of measured structure factors, and *N*_{p} is the number of parameters. If known from a previous crystallographic least-squares refinement, *N*_{p} should be the number of coordinates and displacement parameters of the atoms. However, normally *N*_{p} is simply set equal to 1. In most cases, the wavefunction to determine has the form of a single Slater determinant and the molecular orbital coefficients are adjusted during the fitting. To accommodate atomic vibrations via displacement parameters in the calculated structure factors, the probability distribution functions of the atoms are normally described as one-center functions using Hirshfeld’s stockholder partitioning,^{20} whereby approximations to two-center functions are also available.^{21–23}

The mixing factor *λ* in the functional *L* between the quantum-mechanical energy and the x-ray structure factors is usually increased stepwise to avoid convergence problems. However, it is still debated what the most meaningful ultimate value of *λ* is to terminate the fit, known as the “halting problem” of XCW fitting.^{20,24,25} The relationship between the ultimate *λ* and the desired GooF^{2} value is also in the center of the debate about the most meaningful descriptor for the method, namely, whether the terms “constrained” in the sense of the above-mentioned constrained search methods, “restrained” in a crystallographic sense or “regularized” are more appropriate.^{20,24,26} In this paper, we take the pragmatic approach to both problems discussed in this paragraph in that we always fit to the biggest value of *λ* before convergence ceases, and in that, we adopt the more traditional term XCW fitting for the method. For pertinent details, we refer to the most recent review about XCW fitting.^{20}

XCW fitting has been used to experimentally reconstruct properties such as refractive indices,^{27} molecular polarizabilities and hyper-polarizabilities,^{28} or bond indices.^{29–31} XCW fitting has also been modified by coupling to extremely localized molecular orbitals (XC-ELMOs)^{32,33} and extended to multi-determinant wavefunctions to obtain weights of Lewis resonance structures from x-ray diffraction measurements.^{34,35} Despite all these developments over the past two decades, there are fewer studies exploring what the “enhancement of the information contents” achieved by the fitting process physically means and to what extent physical effects can be extracted from the information-rich but error-prone experimental structure factors. These studies are discussed below.

If one uses a HF wavefunction for an isolated molecule as ansatz for the fitting, there are four main effects present and manifest in the experimental x-ray diffraction structure factors that are absent in the HF wavefunction before the XCW fitting procedure starts.

*Electron correlation.*Here, we refer to dynamic electron correlation, i.e., the instantaneous distribution of electrons within a molecule, and understand the effect of electron correlation as the difference between the true interaction of the electrons that is, in principle, present in the measured data and the inadequate treatment of electron–electron interactions of the Hartree–Fock Hamiltonian because of its use of an average electronic potential instead of the instantaneous potential. The effect of electron correlation on the electron-density distributions of molecules has been studied widely.^{36–43}In the context of this study, Bartell and Gavin were the first to discuss the role of electron correlation in x-ray diffraction data specifically.^{44–47}Later, Chandler*et al.*investigated whether electron correlation can be extracted from simulated theoretical structure factors via XCW fitting.^{48}They concluded “that it is unlikely that x-ray experiments performed under current limitations would be capable of detecting correlation effects in small molecules.” Using a similar approach based on theoretical structure factors, Genoni*et al.*found that “the effect of electron correlation on electron density can be partially captured by the XC-WF method.”^{49}It has also been studied how regions of electron density localization and concentration in the urea crystal are affected by electron correlation introduced via XCW fitting.^{50}However, it has not been investigated so far what the effect of the constrained fit would be on the exchange–correlation potential directly, not only on the electron density. In this study, we aim to fill this gap. In future studies, this could be exploited as the starting point for the design or benchmarking of DFT functionals.*Polarization.*The electron density distribution of a molecule is polarized by its neighbors in a crystal via the composite electric field imposed by them, which is called the crystal effect. The impact of this effect on the electron density, i.e., the difference between the isolated molecular electron density and that in the electric field resulting from the surrounding molecules, has been referred to as interaction density and was investigated theoretically by comparing results from HF periodic-boundary and isolated-molecule calculations^{51–53}or using QM/MM methods.^{54,55}Interaction-density investigations based on the reconstruction of the electron density from theoretical structure factors using the multipole model qualitatively retrieved the features of polarization of the electron density.^{56,57}This study was then extended to experimental data, still utilizing the multipole and the related invariom model.^{58}In the context of XCW fitting, Ernst, Genoni, and Macchi stated that the method is “able to account for the electron density and molecular orbital deformations caused by intra-crystal electric fields,” but the effect becomes less clear when experimental data are used.^{59}This has been shown once more recently.^{60}*Relativistics.*Relativistic effects become important for heavy elements; hence, they are not part of the current study that uses urea and L-alanine as target compounds. However, the possibility of extracting relativistic effects from structure factors using the XCW fitting approach was discussed in theoretical investigations of Bi-, Hg-, and Au-containing compounds and was put into perspective of electron correlation and polarization.^{61–64}*Experimental error and noise.*Beside the discussed effects, experimental structure factors always carry systematic and statistic errors within them and are also modulated by atomic and lattice harmonic and anharmonic vibrations. All of these will also be fitted to the wavefunction and will show in the final electron-density distributions. Hence, a particular importance lies in the data integration and reduction procedure as well as in the analysis of the data and the refinement strategy.^{65,66}Some of us have shown previously that if data of high quality are used (example of xylitol), the difference density due to the XCW fitting agrees well with theoretical reference calculations.^{67}On the contrary, if data of lower quality are used (radiation damage problems in an epoxysuccinyl amide), the XCW fitted effect is significantly different than that of theoretical reference calculations.^{54}Landeros-Rivera, Contreras-García, and Dominiak recently argued that this is due to the sensitivity of the XCW fitting procedure to the data treatment, which is particularly reflected in the experimental standard uncertainties.^{68}Subsequently, in a study using 14 different datasets of oxalic acid dihydrate (none of them of the usual high quality used for experimental charge-density determinations), some of us have tested the reproducibility and reliability of the information extracted from an XCW fitting.^{25}The result was that there were consistent systematic effects preserved in every dataset, some of which could be rationalized to be electron correlation and polarization effects, but others must have been caused by systematic errors in the data. We note that XCW fitting always deals with information that are just above the noise level as discussed in Ref. 26.

Based on these findings and problems, we have chosen to use x-ray diffraction datasets that are well-known to be of superb quality and have been used for similar proof-of-concept studies in the past: urea measured by Birkedal *et al.*^{69} and used, e.g., in Refs. 70–72, as well as L-alanine measured by Destro *et al.*^{73} and used, e.g., in Refs. 74–78. Here—unlike in any other XCW fitting study that discussed the effects I to III—after a theoretical evaluation of the effects I and II, we proceed to the experimental datasets of urea and L-alanine with real experimental errors: In the first and second part of this study, we will investigate to which extent and percentage XCW fitting is able to retrieve electron correlation and polarization from theoretical x-ray structure factors that include either of these two effects exclusively. We will also analyze quantitatively and qualitatively how pronounced the two effects are in comparison to each other. In the third part of this study, we will use the findings of the first and second part as references for the effects extracted from experimental data. Ultimately, this should answer the questions if and to which extent we can separately obtain electron correlation and polarization from experiment.

## II. EXPERIMENTAL AND COMPUTATIONAL SECTION

### A. Clarification of concepts to facilitate separation of electron correlation and polarization

The

*atom-based procrystal*density is an approximation of the electron density of a crystal built up by summation over spherically averaged atomic electron densities from*ab initio*calculations, each centered on the coordinates of the corresponding nuclei within the unit cell of a crystal.^{79–81}The

*molecule-based procrystal*density is an improved approximation of the electron density of a molecular crystal built up by summation over the total non-spherical electron densities of the molecules in the unit cell previously obtained from*ab initio*calculations on the isolated species. This concept has been used in theoretical chemistry under different names and approximations.^{82–84}The

*embedded molecule-based procrystal*density is an even better approximation of the electron density of a molecular crystal, since the constituting molecular electron densities are obtained from previous*ab initio*calculations with an embedding given by self-consistent point charges and dipoles placed at symmetry-generated positions of the surrounding molecules to simulate the crystal effect.^{85}Although there are more advanced techniques available (such as the use of ELMOs for the simulation of the environment),^{86}these have not been tested in the XCW framework yet.The

*reconstructed crystal*density is not a quantum mechanical approximation but a model of the crystal electron density that describes the measured diffraction pattern. In this study, we use the reconstructed crystal density after XCW fitting.

In (iii), the effect of *polarization* is already approximately included in the theoretical HF ansatz before fitting. In this study, we will test whether this approximation is sufficient or whether fitting to theoretical structure factors from periodic-boundary calculations will add additional information to the fitted wavefunction beyond the ansatz. If the approximation is sufficient, then the difference between the reconstructed crystal density (iv) and the embedded molecule-based procrystal density (iii) is an approximate experimental measure of *electron correlation* (plus experimental error and noise).

The method for the XCW fitting can be changed from HF to DFT (here we use the functional BLYP^{87,88}) so that *electron correlation* is already approximately included in the wavefunction before the fitting. We will test in this study whether electron correlation as approximately included in BLYP is sufficient or whether fitting to theoretical structure factors from highly correlated coupled-cluster calculations will add additional information to the fitted wavefunction beyond the ansatz. If the approximation in BLYP is sufficient, then the difference between the reconstructed crystal density (iv) and the BLYP molecule-based procrystal density (ii) is an approximate experimental measure of *polarization* (plus experimental error and noise).

Any difference between the reconstructed crystal density (iv) and the BLYP embedded molecule-based procrystal density (iii) will include only experimental error and noise plus any other unmodeled effect, e.g., if the cluster charges and BLYP approximations are not sufficient. We term this remaining density the *defect* density, which cannot be clearly interpreted anymore physically. However, it is known that, e.g., for highly correlated inorganic materials, a DFT approximation of electron correlation is insufficient, so the information in the defect density might still be valuable to both experimentalists and theoreticians.^{89,90}

### B. Hirshfeld atom refinement of urea and L-alanine

The crystallographic information files (CIFs) and structure factor lists (HKLs) of urea^{69} and L-alanine^{73} as available from the literature were used as input files for a Hirshfeld atom refinement (HAR)^{70,85} as implemented in the software *TONTO*^{91} using the restricted HF method and both the *def*2*-TZVP* and *pob-TZVP* (with the exact coefficients as implemented in *CRYSTAL*14)^{92} basis sets for both compounds. A self-consistent cluster of point charges and dipoles around the central molecule for every symmetry-generated molecule within a radius of 8 Å was used in the crystallographic refinement to simulate the crystal field. This procedure ensures the best possible experimental structures available today from x-ray diffraction data, including A–H bond lengths similar to those derived from neutron-diffraction studies and anisotropic displacement parameters for hydrogen atoms.^{93,94} Figure 1 shows the experimental structures after HAR at the HF/pob-TZVP level of theory, with refinement details summarized in the caption. The resulting HAR CIFs at this level of theory are deposited with the Cambridge Structural Database under the deposition numbers CCDC-2182609 (urea) and 2182613 (L-alanine) and can be obtained free of charge via https://www.ccdc.cam.ac.uk/structures. They were used as starting points for the XCW fittings against the experimental structure factors. For the calculation of theoretical structure factors, the atomic coordinates obtained from HAR were fixed. Note that urea is planar in its crystal structure, whereas it is twisted when optimized quantum-mechanically as an isolated molecule (compare Ref. 49).

### C. Electron correlation

The software *Gaussian*09^{95} was used to calculate the reference electron correlation effects for urea and L-alanine. The coordinates of the HAR-refined final geometries were used as input geometries for coupled-cluster singles and doubles (CCSD) (no frozen core) and HF single-point calculations, respectively. The basis sets *def*2*-TZVP* and *pob-TZVP* (with the exact coefficients as implemented in *CRYSTAL*14)^{92} were used to test basis-set dependencies of the electron correlation effects. However, the basis set comparison is only discussed in the supplementary material. In the following, all results are described based on the *pob-TZVP* basis set only, which is implemented in *CRYSTAL*14 and was explicitly defined by us for the software programs *TONTO* and *Gaussian*09.

The software *denprop*^{96} was employed to calculate the hkl indices and static structure factor magnitudes that include the effect of electron correlation. For this purpose, the *Gaussian*09 wavefunction output (.wfn file) of the CCSD calculation was used as input for a model pseudo-periodic system with cell constants of 10 × 10 × 10 Å^{3} (cubic cell) and a resolution of 0.7, 1.4417, and 2.0 Å^{−1} for urea and 20 × 20 × 20 Å^{3} (cubic cell) with a resolution of 0.7, 1.0778, and 2.0 Å^{−1} for L-alanine. In each case, the medium resolution agrees with the resolution of the experimental .hkl file. Effects of atomic displacement were not simulated, and the artificially large unit cell constants prevent intermolecular interactions to play any role.

The software *TONTO*^{91} was used for the XCW fitting against the theoretical structure factors. A .cif file was manually written with the coordinates from the *Gaussian*09 .wfn file and the cell as specified in *denprop.* It is worth noting that the “thermal smearing model” in the *TONTO* input file needs to be set to *none*. The XCW fitting was performed from λ = 0 to the indicated maximum λ-values in steps of 0.5 (for resolutions of 1.4417 and 2.0 Å^{−1} for urea as well as 1.0778 and 2.0 Å^{−1} for L-alanine) or 0.2 (0.7 Å^{−1} for urea and L-alanine) applying a uniform structure factor error of 0.1. The fitting was done with the HF and BLYP methods as well as both basis sets, but only the pob-TZVP results are discussed in the main text. The maximum λ-values were obtained by monitoring until which λ-value the fitting converged within a given number of cycles or when λ = 10 was reached. Since the relationship between the statistical agreement GooF^{2} (against which the XCW fitting procedure is carried out) and the standard uncertainty is $GooF2\u223c1\sigma 2$ (see Sec. I), a 100 times bigger multiplier λ is necessary to realize the same perturbation when σ of 1.0 is used (as was done, for example, in the work of Genoni *et al.*^{49}) instead of σ = 0.1 (as we do here). This means that with σ = 0.1, λ = 10 represents a high perturbation. It would equal λ = 1000 with σ = 1.0 as applied in Ref. 49. Other models of error distribution are possible and were recently discussed and tested in Refs. 62 and 97.

The determination and visualization of electron correlation effects were done in terms of calculating the differences between electron density grid files of different λ-values and of the reference *Gaussian*09 calculations in the Gaussian cube format. The reference electron correlation effect was calculated as the difference between the CCSD electron density grid (Gaussian09_CCSD.cube) and the HF electron density grid (Gaussian09_HF.cube). The order of subtraction is always CCSD minus HF. The *Gaussian*09 grid files were calculated with the *cubegen* utility, and the *ntps* option was set to −1, which allowed for specifying the grid to be exactly the same as in the *TONTO*-specified grids. The distance between grid points is 0.0472 a.u. in all grid files. The XCW-fitted electron correlation was calculated as the difference electron density between the fitted wavefunction at the value λ = x (TONTO_λ=x.cube) and the unconstrained wavefunction with λ = 0 (TONTO_λ=0.cube, which is identical to Gaussian09_HF.cube). In these cases, the order of subtraction is always λ = x − λ = 0.

For the difference electron density grids, each grid point (r_{i}) of one of the grid files was subtracted from the same grid point (r_{i}) of the second grid file. The obtained grid points containing the difference electron density were then written to a file in the same Gaussian cube file format. Restrictions to this script are that the geometry of the molecule and the orientation of that molecule toward the cube borders must be exactly the same in both cubes. With a second script, the absolute values of each grid point (r_{i}) of the difference density grids were summed up and then divided by two to obtain the total electron count of the difference, namely, the total magnitude of the electron correlation effect as it impacts on the electron density ($Ne=(\u2211i=1np\rho method1(ri)\u2212\rho method2(ri))/2$). To obtain the percentage reconstruction (PR) of the XCW fitting procedure [Eq. (3)], the difference between the TONTO_λ=x.cube with electron density values *ρ*_{λ=x}(*r*_{i}) and the unconstrained values *ρ*_{λ=0}(*r*_{i}) in the TONTO_λ=0.cube was divided by the reference electron density effect (CCSD minus HF). To obtain the real-space R-value (RSR)^{98} [Eq. (4)], the XCW-fitted electron density values *ρ*_{λ=x}(*r*_{i}) were directly referred to the CCSD calculation. Any difference grid files were visualized with the software *VMD*^{99} as isosurface plots,

To judge the trend of the PR as a function of λ, a series of exponential functions of the form

was fitted to the observed course of PR vs λ, giving rise to an estimate of the expected convergence limit at *λ* → ∞. This form of the function was chosen as it fitted the curves with the best R-values of all tested functions, and the limiting behavior can easily be interpreted.

For the visualization of the effect of XCW fitting on the exchange–correlation (XC) potential, the XC potential was calculated with Tonto in the BLYP approximation and saved in the form of a Gaussian cube file. The difference XC potential was calculated as the difference between the fitted wavefunction at the value λ = x and the unconstrained wavefunction with λ = 0 with the order of subtraction always being λ = x − λ = 0. Each grid point of one of the grid files was subtracted from the same grid point of the second grid file. The obtained grid points containing the difference XC potential were then written to a new file in the same Gaussian cube format. These files were used for visualization. No further numerical analysis was undertaken.

### D. Polarization

The software *CRYSTAL*14^{100} was used to calculate the reference polarization effects of urea and L-alanine. The coordinates of the HAR were used as input geometry for a fully periodic single-point calculation at the experimental crystal symmetry and crystal lattice with the HF/*pob-TZVP* level of theory, termed the 3*D* model. In a separate calculation, the keyword *MOLSPLIT* was used to obtain the molecule-based procrystal electron density that does not include polarization effects.

The software *denprop* was used only to calculate the list of hkl indices for a certain resolution, given the experimental lattice parameters. The considered resolutions are 0.7, 1.4417, and 2.0 Å^{−1} for urea as well as 0.7, 1.0778, and 2.0 Å^{−1} for L-alanine. The predicted hkl indices were then used as input (.d3 file) for the calculation of the respective structure factor magnitudes for the 3*D* model with the *runprop*14 utility as implemented in *CRYSTAL*14. The effects of atomic displacement were not simulated; the resulting structure factor magnitudes are static.

The software *TONTO* was used for the XCW fitting, and the procedure is similar to the one described for obtaining the correlation effects. The *CYRSTAL*14-obtained structure factors of the 3*D* model served as input to be fitted against with a molecular HF/pob-TZVP wavefunction as ansatz. The XCW fitting was performed from λ = 0 to the indicated maximum λ-values in steps of 0.5 (for resolutions of 1.4417 and 2.0 Å^{−1} for urea), 0.2 (for resolutions of 0.7 Å^{−1} for urea as well as 1.0778 and 2.0 Å^{−1} for L-alanine), or 0.02 (0.7 Å^{−1} for L-alanine) applying a uniform structure factor error of 0.1. Fitting was stopped at the maximum value of λ for which convergence was achieved or at λ = 10. The same XCW fitting at the HF/pob-TZVP level was also performed upon introducing surrounding cluster charges mimicking the crystal environment. By using such a self-consistent cluster of point charges and dipoles, polarization effects were already simulated in the wavefunction ansatz at λ = 0.

The reference polarization effect on the electron density (i.e., the *interaction density*) was calculated as the difference electron density in the form of a grid file between the fully periodic calculation (CRYSTAL14_3D.cube) and the molecule-based procrystal electron density (CRYSTAL14_MOLSPLIT.cube). The order of subtraction is always 3*D* minus *MOLSPLIT*. *CRYSTAL*14 grid files were obtained with the *runprop*14 utility and the *ECH*3 keyword in the input (.d3) file. The XCW-fitted polarization was calculated as the difference electron density between the fitted wavefunction at the value λ = x (TONTO_λ=x.cube) and the unconstrained wavefunction with λ = 0 (TONTO_λ=0.cube) for both cases with and without cluster charges in the ansatz. The order of subtraction is always λ = x minus λ = 0. The distance between grid points obtained from *CRYSTAL*14 is 0.0527 a.u. (urea) and 0.0563 a.u. (L-alanine), and from *TONTO*, it is 0.0472 a.u. for both molecules. RSR values could not be calculated because of the lack of the definition of a molecule inside the periodic crystal. The calculation of PR values was similar to the case of electron correlation,

For visualization of the *CRYSTAL*14 interaction density, symmetry-generated molecules and their electron densities had to be removed from the supercells that constitute the grid files; otherwise, the view onto the central molecule would be obstructed. Here, we programmed and followed a procedure based on the quantum theory of atoms in molecules (QTAIM) partitioning of space originally published in Ref. 101 to keep only the desired molecule and its electron density in the supercell. The required routines have been implemented in the software cuQCT^{54} for this study. Subsequent visualization of isosurfaces was carried out with *VMD*.

### E. Fitting against experimental structure factors

To evaluate electron correlation and polarization effects from the experimental structure factors of urea and L-alanine, XCW fittings were performed using the pob-TZVP basis set starting with the HAR-obtained .cif-files and the experimental .hkl-files in λ-steps of 0.02 (for urea) and 0.01 (for L-alanine). For each set of structure factors (maximum experimental resolutions 1.4417 Å^{−1} for urea and 1.0778 Å^{−1} for L-alanine as well as each pruned to the resolution of 0.7 Å^{−1}), four fittings were carried out: (1) HF with cluster charges (cc), (2) HF without cc, (3) BLYP with cc, and (4) BLYP without cc. The effects represented by the four fittings always refer to difference electron densities from calculations at λ = x and λ = 0.

## III. RESULTS AND DISCUSSION

### A. Electron correlation

Figure 2 presents the theoretical reference maps for the XCW fitting. Blue regions indicate electron gain due to the inclusion of electron correlation into the method. These are around the core regions of the atoms. The bonds and lone pairs lose electron density upon the inclusion of electron correlation, so overall it is a redistribution of electron density from valence to core regions. The effect is very systematic, but not large for organic molecules: the isovalue is 0.034 *e* Å^{−3}, and the total number of electrons shifted is 0.27 (urea) and 0.40 (L-alanine), which corresponds to exactly 0.067 electrons per non-H atom (C, N, O) in both cases. The direction, distribution, and magnitude of the effect of electron correlation on the electron density is consistent with earlier theoretical studies.^{38,43,49,67}

Figures 3–5 and Table I refer to the XCW fitting results for both compounds against the theoretical CCSD structure factors at the intermediate resolution to test if and to which extent the electron correlation effect can be recovered by the fitting technique. As mentioned in Sec. I, in Refs. 48 and 49, the authors doubt that XCW fitting can recover it, at least not fully. Therefore, we first observe in Figs. 3(a)–3(c) and 4(a)–4(c) how the effect of electron correlation is built up in the single determinantal HF wavefunction upon increasing the λ value stepwise. Certainly, the features match those of the reference calculations in Fig. 2, and they become more intense with increasing λ value. At the maximum λ value, λ_{max} = 10.0, 88.4% and 94.1% of the effect has been recovered for both compounds, respectively.

Resolution (Å^{−1})
. | Ansatz
. | a_{1}
. | b_{1}
. | c (=PR_{∞})
. | PR_{max}
. | RSR_{max}
. | $GooFmax2$ . | λ_{max}
. |
---|---|---|---|---|---|---|---|---|

0.7 | HF (urea) | −12.51 | 4.47 | 100.5 | 94.3 | 0.001240 | 0.000302 | 3.2 |

1.4415 | −34.14 | 5.49 | 93.9 | 88.4 | 0.001460 | 0.000179 | 10.0 | |

2.0 | −49.35 | 8.07 | 88.6 | 74.1 | 0.002443 | 0.000302 | 10.0 | |

0.7 | BLYP (urea) | 15.79 | 0.43 | 109.1 | 109.6 | 0.001632 | 0.000686 | 2.0 |

1.4415 | 15.00 | 4.74 | 108.2 | 110.2 | 0.001691 | 0.000135 | 9.5 | |

2.0 | 20.78 | 6.80 | 110.3 | 115.1 | 0.002141 | 0.000201 | 10.0 | |

0.7 | HF (L-alanine) | −35.33 | 0.26 | 94.9 | 92.8 | 0.001257 | 0.000699 | 2.6 |

1.0778 | −38.47 | 0.09 | 96.0 | 94.1 | 0.001127 | 0.000186 | 10.0 | |

2.0 | −32.56 | 1.77 | 86.4 | 74.1 | 0.002500 | 0.000448 | 10.0 | |

0.7 | BLYP (L-alanine) | 16.77 | 0.43 | 109.1 | 109.9 | 0.001773 | 0.001142 | 1.8 |

1.0778 | 16.53 | 0.80 | 109.8 | 110.0 | 0.001736 | 0.000392 | 4.5 | |

2.0 | 15.99 | 4.15 | 113.7 | 116.5 | 0.002301 | 0.000313 | 10.0 |

Resolution (Å^{−1})
. | Ansatz
. | a_{1}
. | b_{1}
. | c (=PR_{∞})
. | PR_{max}
. | RSR_{max}
. | $GooFmax2$ . | λ_{max}
. |
---|---|---|---|---|---|---|---|---|

0.7 | HF (urea) | −12.51 | 4.47 | 100.5 | 94.3 | 0.001240 | 0.000302 | 3.2 |

1.4415 | −34.14 | 5.49 | 93.9 | 88.4 | 0.001460 | 0.000179 | 10.0 | |

2.0 | −49.35 | 8.07 | 88.6 | 74.1 | 0.002443 | 0.000302 | 10.0 | |

0.7 | BLYP (urea) | 15.79 | 0.43 | 109.1 | 109.6 | 0.001632 | 0.000686 | 2.0 |

1.4415 | 15.00 | 4.74 | 108.2 | 110.2 | 0.001691 | 0.000135 | 9.5 | |

2.0 | 20.78 | 6.80 | 110.3 | 115.1 | 0.002141 | 0.000201 | 10.0 | |

0.7 | HF (L-alanine) | −35.33 | 0.26 | 94.9 | 92.8 | 0.001257 | 0.000699 | 2.6 |

1.0778 | −38.47 | 0.09 | 96.0 | 94.1 | 0.001127 | 0.000186 | 10.0 | |

2.0 | −32.56 | 1.77 | 86.4 | 74.1 | 0.002500 | 0.000448 | 10.0 | |

0.7 | BLYP (L-alanine) | 16.77 | 0.43 | 109.1 | 109.9 | 0.001773 | 0.001142 | 1.8 |

1.0778 | 16.53 | 0.80 | 109.8 | 110.0 | 0.001736 | 0.000392 | 4.5 | |

2.0 | 15.99 | 4.15 | 113.7 | 116.5 | 0.002301 | 0.000313 | 10.0 |

Figure 5(a) shows the curves of GooF^{2}, RSR, and PR as a function of the increasing λ values. GooF^{2}, a measure in reciprocal space, and RSR, a measure in real space, both approach zero with increasing λ, as expected, but do not reach it. Simultaneously, PR approaches 100% but does not reach it. To find out whether, in principle, a 100% reconstruction of the calculated CCSD effect would be possible with the HF ansatz, we fitted a curve through the PR vs λ graph [Eq. (5), discussed in Sec. II, gray line in Fig. 5(a)] and checked the limiting behavior $lim\lambda \u2192\u221ePR\lambda =c$ (see also Table I). For urea, 93.9% can be reached as a limit, whereas for L-alanine, 96.0% can be reached. Hence, the values of 88.4% and 94.1% at λ_{max} = 10.0 are already close to the optimum reconstruction. Although 100% cannot be reached, most of the effect of electron correlation can be recovered, and qualitatively, all the chemical features in lone pair, bond, and core regions are reproduced upon XCW fitting.

There can be several reasons why the reference effect was largely but not fully recovered in this theoretical experiment: (i) Limited resolution when the wavefunction is projected into structure factors, (ii) insufficiently flexible method, (iii) insufficiently flexible basis set. (i) will be investigated further below. For (iii), it has been shown before that triple-zeta basis sets are sufficiently flexible for XCW fitting.^{102} We also investigate the effect of different triple-zeta basis sets in the supplementary material (def2-TZVP vs pob-TZVP) and find that these differences are negligible. Hence, we test here whether the introduction of electron correlation into the method by means of the DFT functional BLYP leads to the same or a similar result when fitted to the CCSD structure factors. In fact, Figs. 3(d) and 4(d) show that already at λ = 0, before the fitting procedure starts, the DFT electron correlation map overestimates the reference CCSD map significantly. All features are more intense [comparison of Fig. 3(d) with Fig. 2(a) and Fig. 4(d) with Fig. 2(b)] and the PR values are 156.5% (urea) and 159.4% (L-alanine), respectively.

With these BLYP wavefunctions that overestimate electron correlation at the starting point already, the XCW fitting counterbalances this overestimation and corrects the BLYP wavefunctions, approaching PR of 100% as in the case of HF. Figures 3(f) and 4(f) show that the final maps are qualitatively very similar to the CCSD reference and to the HF-fitted maps, and the final PR values are 110.2% (urea) and 110.0% (L-alanine). The same type of extrapolation for $lim\lambda \u2192\u221ePR\lambda =c$ using Eq. (5) [gray lines in Fig. 5(b)] leads to limiting values of 108.2% and 109.8% (Table I), showing again that the vast majority of the electron correlation effect has already been included up to *λ*_{max} but also that 100% cannot be reached in principle. Nevertheless, this result is promising as we have shown here for the first time that XCW fitting can systematically and significantly correct DFT functionals toward a better description of electron correlation. How this manifests itself in the exchange–correlation (XC) potential will be shown below.

Before we turn to XC potentials, we will discuss the resolution dependence of the results with the help of Table I. Intuitively, one might think that reconstructions ought to be more successful if higher resolution data are used because truncation of resolution means information loss in general. Indeed, the λ_{max} values in Table I indicate that with less data, XCW fitting becomes more problematic, as at the lowest resolution of 0.7 Å^{−1}, the fitting never converges above λ_{max} = 1.8–3.2. In contrast, the reconstruction of electron density in real space is significantly more successful if only low-resolution data are used. For urea, $lim\lambda \u2192\u221ePR\lambda =PR\u221e=100%$ can be reached, and in contrast, for the high resolution of 2.0 Å^{−1}, PR_{∞} and PR_{max} are only 88.6% and 74.1%, respectively. The same trend is present in the RSR values. This resolution dependence also holds for the BLYP results, although here 108.2% is the optimum reconstruction, found in urea at the intermediate resolution. The leading terms of the exponential fit a_{1} and b_{1} support the trend, with the magnitude of the values being (mostly) the largest at high resolution and the smallest at low resolution. This corresponds to a slower convergence with respect to increasing λ at higher resolution. This means that much higher values of λ would be needed to observe the same effect as at lower resolution. This has to do with the non-linear increase of the number of reflections with the expansion of the Ewald sphere toward higher resolution so that weak-intensity reflections start to dominate the fitting and attenuate the information present in the high-intensity reflections in the low order (see also Ref. 49). It is worth noting that the GooF^{2} values do not describe the same observation; they are always significantly larger at low resolution (Table I).

### B. Exchange–correlation potential

The XC potential is a scalar field in real space, so it can be evaluated on a grid in the same way as the electron density. However, unlike the electron density, a universal form of the XC potential is inexistent, so it needs to be taken in the definition of a certain functional. Here, for consistency, we use the BLYP functional.^{87,88} Since the XC potential is not confined to have a fixed integral over a certain volume of space, the calculation of the RSR or PR values would be on an arbitrary and ambiguous scale. Furthermore, the XC potential is directly related to the electron density, whose changes upon fitting have been quantified above. Therefore, only the qualitative effects on the XC potential will be discussed in this subsection and only for the urea molecule.

Starting with HF densities to improve DFT functionals is not an entirely new idea. In density-corrected DFT (DC-DFT) or HF-DFT, a HF density is simply given to DFT functionals without further modification, and the results are convincing if used in the right context.^{103,104} From our study, we already know (see above) that upon XCW fitting, the HF wavefunction absorbs nearly all of the physically correct electron correlation effects (compared to the CCSD reference calculation), which we have validated above in terms of the electron-density distribution. Therefore, calculating the BLYP XC potential from a fitted HF wavefunction at high λ should be physically meaningful. Table S1 shows that at λ = 2.5, the HF wavefunction has already reconstructed 66% of the electron density features assigned to the electron correlation effect. Therefore, we start to observe the XC potential differences upon fitting in the HF framework in Fig. 6(a) from λ = 2.5. The main features are in the valence region of the atoms, specifically in the oxygen lone pairs and around the carbon atom. With increasing λ [Figs. 6(b) and 6(c)], features appear in the C–O and C–N bonds as well as becoming more diffuse overall.

When we use the BLYP functional for XCW fitting [Figs. 6(d)–6(f)], the picture is entirely different. The correction that the XCW fitting introduces into the BLYP Kohn-Sham-type wavefunction manifests itself only in the core region of the atoms. We know from Figs. 3(d)–3(f) that the general overestimation of the electron correlation effect by the BLYP functional specifically means too low electron density in the core regions [extension and intensity of positive (blue) isosurfaces reduce from Figs. 3(d)–3(f)] and too high electron density in the valence regions [extension and intensity of negative (red) isosurfaces reduce from Figs. 3(d)–3(f)]. This means that the electron density is increased upon fitting in the atomic cores. The consequence of this is an increasingly more negative XC potential in the core region as shown by an increase in the extension and intensity of negative (red) isosurfaces from Figs. 6(d)–6(f). Although, in general, the relationship between the XC potential and the underlying electron density is not easy to interpret, it seems that here the technique of XCW fitting allows a systematic improvement of the BLYP functional via the XC potential.

With increasing resolution of included data, the same trend as for the electron density is observed (compare Fig. 6 to Figs. S30 and S32 in the supplementary material): Features in the XC potential are much more pronounced already at small λ values when only low-resolution data are included in the fitting. In other words, to reproduce the same features qualitatively and quantitatively, much higher λ values are needed for higher-resolution structure factor datasets.

### C. Polarization

Figure 7 presents the theoretical electron density reference maps for the XCW fitting in the case of polarization as the difference between HF wavefunctions under periodic-boundary and isolated-molecule conditions, also called interaction density. Blue regions indicate electron gain due to interaction with neighboring molecules, whereas red regions indicate electron loss. In contrast to the electron correlation effect (see Fig. 2), here the core regions are less affected but bond and lone pair regions are both involved in electron gain and loss, featuring dipolar character. The two C=O carbonyl groups in urea and L-alanine, both involved in N–H⋯O hydrogen bonding [top functional group in Fig. 7(a), left functional group in Fig. 7(b)], are polarized recognizably similarly across the two different compounds: at the oxygen atom, there is a negative region close to the atomic core separated from the negative region in the bond closer to the carbon atom by a positive surface spanning from the oxygen lone pairs through the bond in a crescent shape.

Overall, this systematic effect of polarization of the electron density leads in total only to a small redistribution of electrons, with the total number of shifted electrons being 0.25 (urea) and 0.44 (L-alanine). These numbers are nearly the same as for electron correlation, which means that the magnitude of the impact of these two effects on the electron density is nearly the same for small organic molecules, but not the direction and distribution. The direction, distribution, and magnitude of the interaction density observed here are consistent with earlier quantum crystallographic studies.^{54,58–60,105}

Figures 8–10 and Table II refer to the XCW fitting results for both compounds against the theoretical periodic-boundary structure factors at the intermediate resolution to test if and to which extent the interaction density originating from the crystal effect can be recovered by the fitting technique. As discussed in the Introduction, it was concluded in previous studies that the crystal effect can be incorporated efficiently into isolated-molecule wavefunctions by XCW fitting^{60} but electron correlation only partially.^{49} By comparing the build-up of the qualitative features in Figs. 8(a)–8(c) and 9(a)–9(c) with the features in the corresponding reference maps in Fig. 7, it becomes clear that XCW fitting indeed reconstructs the crystal effect accurately. A comparison of the polarization with the electron correlation effect [Figs. 8(a)–8(c) with Figs. 3(a)–3(c) and Figs. 9(a)–9(c) with Figs. 4(a)–4(c)] shows that polarization is reconstructed more quickly, i.e., already at low λ-values, a significant portion of the polarization effect has been recovered, more than for electron correlation [at λ = 1.0: 62.3% vs 44.1% for urea and 61.2% (at λ = 0.2) vs 48.9% (at λ = 0.5) for L-alanine].

Resolution (Å^{−1})
. | Ansatz
. | a_{1}
. | b_{1}
. | c (=PR_{∞})
. | PR_{max}
. | $GooFmax2$ . | λ_{max}
. |
---|---|---|---|---|---|---|---|

0.7 | HF (urea) | N/A | N/A | N/A | 94.0 | 0.000771 | 1.6 |

1.4415 | −32.23 | 0.97 | 96.4 | 92.1 | 0.000150 | 9.5 | |

2.0 | −33.99 | 1.24 | 90.7 | 82.8 | 0.000177 | 10.0 | |

0.7 | HF_CC (urea) | N/A | N/A | N/A | 116.2 | 0.000294 | 1.6 |

1.4415 | N/A | N/A | N/A | 116.2 | 0.000062 | 9.5 | |

2.0 | N/A | N/A | N/A | 117.1 | 0.000052 | 10.0 | |

0.7 | HF (L-alanine) | −34.74 | 0.05 | 96.6 | 89.9 | 0.007940 | 0.4 |

1.0778 | −29.52 | 0.26 | 88.6 | 86.7 | 0.003566 | 1.0 | |

2.0 | −35.17 | 0.93 | 95.9 | 88.9 | 0.000498 | 7.4 | |

0.7 | HF_CC (L-alanine) | N/A | N/A | N/A | 120.2 | 0.003277 | 0.46 |

1.0778 | 3.02 | 0.01 | 119.0 | 120.1 | 0.001114 | 1.6 | |

2.0 | N/A | N/A | N/A | 120.9 | 0.000271 | 6.8 |

Resolution (Å^{−1})
. | Ansatz
. | a_{1}
. | b_{1}
. | c (=PR_{∞})
. | PR_{max}
. | $GooFmax2$ . | λ_{max}
. |
---|---|---|---|---|---|---|---|

0.7 | HF (urea) | N/A | N/A | N/A | 94.0 | 0.000771 | 1.6 |

1.4415 | −32.23 | 0.97 | 96.4 | 92.1 | 0.000150 | 9.5 | |

2.0 | −33.99 | 1.24 | 90.7 | 82.8 | 0.000177 | 10.0 | |

0.7 | HF_CC (urea) | N/A | N/A | N/A | 116.2 | 0.000294 | 1.6 |

1.4415 | N/A | N/A | N/A | 116.2 | 0.000062 | 9.5 | |

2.0 | N/A | N/A | N/A | 117.1 | 0.000052 | 10.0 | |

0.7 | HF (L-alanine) | −34.74 | 0.05 | 96.6 | 89.9 | 0.007940 | 0.4 |

1.0778 | −29.52 | 0.26 | 88.6 | 86.7 | 0.003566 | 1.0 | |

2.0 | −35.17 | 0.93 | 95.9 | 88.9 | 0.000498 | 7.4 | |

0.7 | HF_CC (L-alanine) | N/A | N/A | N/A | 120.2 | 0.003277 | 0.46 |

1.0778 | 3.02 | 0.01 | 119.0 | 120.1 | 0.001114 | 1.6 | |

2.0 | N/A | N/A | N/A | 120.9 | 0.000271 | 6.8 |

Nevertheless, the initially steep slope for PR vs λ flattens so that PR converges for polarization toward a value lower than 100% (Fig. 10), similar to electron correlation (Fig. 5). In fact, Table II shows that the limiting PR value from a fit of PR vs λ according to Eq. (5) is PR_{∞} = 96.4% for urea and 88.6% for L-alanine. For electron correlation, these values are 93.9% and 96.0% (see Table I). This means that both effects can be reconstructed by XCW fitting as accurately and efficiently, only polarization builds up at low values of λ more quickly. The consequence of this is, in turn, that electron density maps generated from XCW fitting to experimental data, as shown in Sec. III D, should be dominated by polarization at low values of λ and electron correlation at high values of λ.

For electron correlation, we have investigated how the inclusion of the effect into the wavefunction ansatz by means of the DFT BLYP functional impacts on the fitting. Electron correlation is overestimated in BLYP and is corrected upon fitting toward the reference values. The analog of this theoretical experiment for polarization is the inclusion of a cluster of point charges and dipoles at symmetry-generated positions around the central molecule to simulate the crystal effect (abbreviated CC here for “cluster charges”).^{85,86} As for electron correlation, the effect is clearly overestimated in the cluster charges [see Figs. 8(d) and 9(d)] with PR values of 115.1% and 131.0%, respectively. In contrast to electron correlation, the effect is not corrected by XCW fitting. The features in the difference density maps of Figs. 8(e) and 8(f) as wells as Figs. 9(e) and 9(f) stay as intense as without fitting. Figure 10 shows that the slope of the PR vs λ curve is essentially horizontal. From these findings, we hypothesize that during XCW fitting against experimental structure factors, cluster charges should not be used if polarization is to be extracted from the experimental data as they bias the result and do not allow to recover the reference effect from the data. This contrasts with electron correlation, where a DFT functional can be used in fitting as it does not hamper to recover the reference electron correlation effect from the data. These hypotheses are tested in Sec. III D. In the future, it should be tested whether a quantum-mechanical model of the crystal effect instead of a model based on point charges behaves differently upon XCW fitting.^{86}

Interestingly, this behavior cannot be seen in the GooF^{2} values that simply follow the λ-values. The higher λ, the smaller GooF^{2} [Figs. 8(d)–8(f) and 9(d)–9(f), Tables S4 and S7]. This was also seen in the discussion of Table I, and it again shows in Table II where overall higher λ_{max} values correspond to lower GooF^{2}_{max} values. This means that a discussion of the halting problem must not only be based on the GooF^{2} vs λ relationship,^{20,24,25} but it should involve derived chemical properties. Here, one could draw a parallel to the variational principle of quantum chemistry: just minimizing the energy does not assure that the determined wavefunction is the best one for other properties. In addition, the role of the distribution of standard uncertainties upon fitting must be clarified in future studies.

The resolution dependence of the fitting of the polarization effect^{59,60} agrees with the trends observed for electron correlation but is less pronounced (see Table II). With low-resolution data, the physical effect is incorporated into the wavefunction at lower λ-values more quickly. From the higher exponential coefficients *b*_{1} of the higher resolution datasets, it is quite clear that these datasets contain less information about polarization per value of λ. This means, on the other hand, that convergence ceases much earlier for low-resolution datasets. The resolution dependence of the maximum achievable PR value on resolution is clearly less pronounced than for electron correlation (compare Table II with Table I).

### D. Experimental XCW fitting

As discussed in the Introduction, XCW fitting using a HF wavefunction should incorporate a superposition of electron correlation (EC), polarization (pol), and experimental errors into the model wavefunction. Hence, one could formulate the difference in density between λ_{max} and λ_{0} in the following way:

As discussed in Sec. II A in detail, it is unclear whether these effects can be separated from each other since they are convolved in the experimental data. In addition, Secs. III A and III C have shown that there will be residual effects left in the fitted wavefunction. Therefore, we will make the following idealized assumption as a working hypothesis and test it further below: If the XCW fittings are performed at the specific levels given in the list, the effects given in bold font should become visible.

BLYP without cluster charges: Δρ

_{pol},**polarization effects only**(correlation is included at λ = 0, but not polarization);HF with cluster charges: Δρ

_{EC},**correlation effects only**(polarization is included at λ = 0, but not correlation);BLYP with cluster charges: Δρ

_{defect},**defect density**(correlation and polarization are included at λ = 0);HF without cluster charges: Δρ,

**all effects**(nothing is included at λ = 0).

Figures 11(a) and 12(a) show the effect of polarization extracted from the experimental structure factors according to the above hypothesis. Indeed, a comparison with Fig. 7 that represents the purely theoretical polarization effect shows some degree of qualitative similarity. The distribution of positive (blue, in the bonds and lone pairs) and negative regions (red, atomic cores) is similar. Even the crescent-shaped positive area at the carbonyl oxygen atom is present in both compounds urea and L-alanine. Only the intensity and extension of both positive and negative effects is more pronounced in the experimental fittings. This is also reflected in the total number of electrons shifted during the fitting (N_{e}), which is 0.32/0.52 *e* (urea/L-alanine), and it was 0.25/0.44 *e* in the reference calculations (Fig. 7). We therefore conclude that the chosen approach indeed qualitatively extracts Δρ_{pol} separately from the experimental structure factors and overestimates it only slightly compared to the purely theoretical reference calculations. In fact, it could conversely mean that the effect is underestimated in the reference calculations and expressed correctly here in the experimental fitting. Evidence for this interpretation is presented at the end of this section.

Figures 11(b) and 12(b) show the effect of electron correlation extracted from the experimental structure factors according to the above hypothesis. Here, comparison goes with the purely theoretical electron correlation effects in Fig. 2. Again, there is qualitative similarity, and it is noteworthy that the observations in Figs. 11 and 12 vs Fig. 2 are completely independent. The typical features of positive (blue) difference electron density around the atomic cores and negative (red) features in the lone pairs and bonds, including the ring-shaped lone-pair region at the oxygen atom in urea, are all reproduced in the experimental maps. The total number of electrons shifted during the fitting (N_{e}) is 0.28/0.42 *e* (urea/L-alanine), and it was 0.27/0.40 *e* in the reference calculations (Fig. 2). Again, we conclude that the chosen approach indeed qualitatively extracts Δρ_{EC} separately from the experimental structure factors, and quantitatively, the effect is only slightly higher in the experimental fitting. It remains to be seen whether this small quantitative difference is significant enough to develop and benchmark DFT functionals against the experimentally fitted electron densities. In any case, some qualitative similarities but different extent of intensity and extension of the features also persist for the exchange–correlation potential, when comparing Fig. 13(a) with Figs. 6(a)–6(c) and Fig. 13(b) with Figs. 6(d)–6(f).

Figures 11(c) and 12(c) show the merger of all effects at *λ*_{max} according to the above hypothesis. The effect is clearly still systematic, not random, so it is not dominated by experimental *random* noise. In fact, all features for both urea and L-alanine resemble the features of electron correlation [Figs. 11(b) and 12(b)] much more closely than those of polarization [Figs. 11(a) and 12(a)] and are consequently not influenced by *systematic* experimental errors either. Clearly, the positive regions of the core electron density and the negative regions in the oxygen lone pairs are signs of residual electron correlation. This means that at *λ*_{max}, electron correlation dominates over polarization; the combined effect is not a combination of the features of both as one might have expected. Comparing Fig. 11(b) with Fig. 11(c) and Fig. 12(b) with Fig. 12(c) shows, though, that the effects are not identical; there are more pronounced and extended features in the maps Figs. 11(c) and 12(c) representing all effects. The N_{e} values are also higher (0.35 vs 0.28 *e* for urea and 0.52 vs 0.42 *e* for L-alanine). This means that there is some added level of noise and the effect of polarization does modulate the electron correlation features, but polarization remains hidden, its features being not discernible.

It was discussed in Sec. III C that polarization builds up more quickly in the difference electron density maps than electron correlation at small values of λ, whereas electron correlation builds up at larger values of λ. Therefore, it could be expected that features of polarization are visible in maps at much smaller values than the maps at *λ*_{max} presented here in Figs. 11(c) and 12(c). Therefore, in the supplementary material (Figs. S29 and S62), we show the evolution of the fitting effect in small steps of λ. These figures show, importantly, and in contrast to the theoretical prediction in Sec. III C, that electron correlation dominates the plots in all ranges of λ, even at very small values.

Finally, we discuss the defect densities shown in Figs. 11(d) and 12(d). Δρ_{defect} will contain information about untreated experimental errors and untreated physical effects, not limited to untreated absorption, extinction, relativistic effects, measurement and machine errors, thermal diffuse scattering, untreated anharmonic atomic motion, etc.^{25} This does explicitly also include those parts of electron correlation and polarization that are treated incorrectly or insufficiently in the model wavefunction ansatz represented by BLYP and cluster charges. Sections III A and III C have shown that both methods, in fact, overestimate the electron correlation and polarization effects relative to more sophisticated theoretical estimates of them (coupled-cluster calculations and periodic-boundary conditions, respectively). Section III A has also shown that XCW fitting mostly corrects this overestimation in BLYP, whereas Sec. III C has shown that this is not true for the overestimation in cluster charges.

Figures 11(d) and 12(d) show clearly systematic features, not just the remaining noise from the experiment, and are very similar to Figs. 11(a) and 12(a), in which the polarization effect was extracted from the experimental structure factors successfully. This means that when cluster charges are used in the model wavefunction, they do not filter out the effect of polarization, but polarization is nevertheless fitted from the experimental data. In turn, this means that the cluster-charge approximation is insufficient, and the experiment provides the real polarization effect, which is absorbed by the fitted wavefunction beyond and despite the cluster charges being present in the ansatz. This agrees with the findings from Sec. III C, where XCW fitting has also shown to be insensitive to the cluster charges. Consequences of this observation are (i) that the cluster-charge approach should be replaced in XCW fitting by a better, more quantum mechanical, approximation—as suggested in Ref. 86 for Hirshfeld Atom Refinement—if the polarization effect is to be filtered out by the ansatz and (ii) that the importance of systematically distributed experimental standard uncertainties affecting the XCW fitting procedure has not been understood correctly yet. Moreover, we cannot exclude that the non-optimal form of the BLYP functional also prevents the defect density to be free from systematic effects. In summary, Δρ_{defect} really represents the defect in the procedure and not a measure of the experimental random and systematic error.

## IV. CONCLUSIONS

### A. Conclusions from XCW fitting against theoretical data

The effect of polarization builds up in the fitted wavefunction more quickly, i.e., at lower values of λ, than electron correlation. Nevertheless, if fitted to λ_{max}, both polarization and electron correlation effects can be recovered to the same extent and with the same accuracy and reliability. An extrapolation to λ_{∞} shows that inherently XCW fitting is not more sensitive to the one or the other effect. The success of the reconstruction of both effects depends on the maximum resolution of the structure factors set; a medium-resolution dataset allows faster and more complete reconstruction than an ultra-high resolution dataset.

We show for the first time how the flow-on effect of the fitting to correlated structure factors on the exchange–correlation potential looks like and that it is physically meaningful, even if a HF wavefunction is used for fitting. Importantly, XCW fitting corrects the overestimation of the electron correlation effect in the used BLYP functional relative to a CCSD reference in both electron density and XC potential. In contrast, an overestimation of the polarization effect by a cluster of point charges and dipoles is not corrected by XCW fitting relative to a fully periodic calculation.

### B. Conclusions from XCW fitting against experimental data

Polarization and electron correlation can be separated and extracted separately from the experimental structure factors by a judicious choice of the wavefunction ansatz upon fitting. The results of the experimental fitting are qualitatively and quantitatively similar to the independent high-level reference calculations used. This is remarkable since it has not been shown before that the effects of polarization and electron correlation extracted from experiment are phenomenologically identical to theoretical assumptions.

The prospect of these results is that XCW fitting can indeed be used as a method to train wavefunctions for different physical effects separately. We have shown that the use of experimental data for this purpose is in view. In order to use this idea directly for the development of DFT functionals, though, a generalized, functional-independent way of extracting the XC potential from experiment is needed. Some of us are working on two different ways of achieving this goal at present.

### C. Conditions and uncertainties

The results in this paper are based on two datasets of high quality, extensively used as test sets for benchmarking new methods in the field of quantum crystallography for decades. In a different study,^{25} some of us have recently demonstrated that with data of not the same high quality, the reproducibility of the fitted effects is limited. This has to do with the fact that the discussed effects are very small. For urea and L-alanine, they amount to roughly 0.05 *e* per non-hydrogen atom per physical effect. Moreover, as a ballpark figure, 95% (corresponding to a crystallographically refined R-value of 5%) of the measured physical effects are reconstructed by simply summing up calculated spherical atomic electron densities, another few percent by accounting for atomic non-sphericity in the chemical and crystal field of neighbors (chemical bonding effect), so only a very small percentage of the reconstructed electron density is due to the discussed effects of polarization and electron correlation. Therefore, as the next step toward the general use of experimental diffraction data in method development for quantum chemistry, a reproducibility study using very high-quality datasets is needed and envisaged. Although the effects are small on an electron-density scale relative to experimental uncertainties, one should not forget that such small permanent, fluctuating or induced electron-density differences give rise to very significant intermolecular forces (electrostatic, van-der-Waals, London dispersion etc.). Roughly, already a change of 1% in the molecular energy of urea or L-alanine amounts to several thousand kJ/mol.

Finally, we have shown that the approximation of the crystal effect with cluster charges is insufficient for XCW fitting so that better quantum-mechanical models must be tested in the future. The appearance of systematic effects in the defect density also shows that we still lack understanding of how the experimental standard uncertainties influence the fitting procedure as they systematically vary in the experiment and can be seen as a second set of independent experimental information. An investigation of the information contents of the uncertainties is therefore also crucial for application of XCW fitting for method development in quantum chemistry.

## SUPPLEMENTARY MATERIAL

The supplementary material for this article is available as a pdf document and includes numerical details of the XCW fittings, difference density plots for each compound at each resolution, as well as a basis-set dependency study.

## ACKNOWLEDGMENTS

S. Grabowsky, E. Hupf, and F. Kleemiss thank the German Research Foundation (Deutsche Forschungsgemeinschaft DFG) for funding within the Emmy-Noether scheme (Grant No. GR 4451/1-1) (S.G.), for a postdoctoral fellowship (Grant No. HU 2512/1-1) (E.H.), and for a Walter–Benjamin fellowship (Grant No. KL 3500/1-1) (F.K.). M. Woińska acknowledges financial support from the Polish National Science Centre (NCN) under Opus Grant No. DEC-2018/31/B/ST4/02142.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

E.H. and F.K. contributed equally to this work.

**Emanuel Hupf:** Conceptualization (supporting), Formal Analysis (equal), Funding Acquisition (supporting), Investigation (equal), Methodology (equal), Software (equal), Visualization (equal), Writing – Original Draft (supporting); **Florian Kleemiss:** Conceptualization (supporting), Formal Analysis (equal), Funding Acquisition (supporting), Investigation (equal), Methodology (equal), Software (equal), Visualization (equal), Writing – Original Draft (supporting). **Tobias Borrmann:** Methodology (supporting), Software (supporting). **Rumpa Pal:** Methodology (supporting). **Joanna M. Krzeszczakowska:** Conceptualization (supporting). **Magdalena Woińska:** Investigation (supporting), Conceptualization (supporting). **Dylan Jayatilaka:** Conceptualization (supporting), Software (equal). **Alessandro Genoni:** Methodology (supporting), Validation (lead), Writing – Review & Editing (lead). **Simon Grabowsky:** Conceptualization (lead), Funding Acquisition (lead), Methodology (equal), Project Administration (lead), Resources (lead), Supervision (lead), Validation (supporting), Writing – Original Draft (lead), Writing – Review & Editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material. In addition, the crystallographic refinement results have been deposited with the Cambridge Structural Database and can be obtained under the deposition numbers CCDC-2182609 (urea) and 2182613 (L-alanine) free of charge via https://www.ccdc.cam.ac.uk/structures.

## REFERENCES

_{4}

_{2}meets relativistic quantum crystallography. How to teach relativity to a non-relativistic wavefunction

_{2}and Ce

_{2}O

_{3}