# The effects of experimentally obtained electron correlation and polarization on electron densities and exchange-correlation potentials

Emanuel Hupf, Florian Kleemiss, Tobias Borrmann, Rumpa Pal, Joanna M Krzeszczakowska, Magdalena Woińska, Dylan Jayatilaka, Alessandro Genoni, Simon Grabowsky

► To cite this version:

Emanuel Hupf, Florian Kleemiss, Tobias Borrmann, Rumpa Pal, Joanna M Krzeszczakowska, et al.. The effects of experimentally obtained electron correlation and polarization on electron densities and exchange-correlation potentials. The Journal of Chemical Physics, 2023, 158 (12), pp.124103. 10.1063/5.0138312. hal-04046695

HAL Id: hal-04046695

https://hal.univ-lorraine.fr/hal-04046695

Submitted on 26 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Effects of Experimentally Obtained Electron Correlation and Polarization on Electron Densities and Exchange-Correlation Potentials

Emanuel Hupf,ᵃ,† Florian Kleemiss,ᵃ,ᵇ,ᶜ,† Tobias Borrmann,ª Rumpa Pal,ª Joanna M.

Krzeszczakowska,ᵈ,ᵉ Magdalena Woińska,ᵉ Dylan Jayatilaka,ᵈ Alessandro Genoni,ᵉ Simon

Grabowskyᵃ,ᵇ,ᵈ,*

ᵃ Institut für Anorganische Chemie und Kristallographie, Fachbereich 2 – Biologie/Chemie,
Universität Bremen, Leobener Straße 7, 28359 Bremen, Germany

ᵇ Departement für Chemie, Biochemie und Pharmazie, Universität Bern, Freiestrasse 3, 3012 Bern,
Switzerland

ᶜ Fakultät für Chemie und Pharmazie, Universität Regensburg, Universitätsstraße 31, 93053
Regensburg, Germany

ᵈ School of Molecular Sciences, The University of Western Australia, 35 Stirling Highway, Crawley
WA 6009, Australia

ᵉ Biological and Chemical Research Centre, Chemistry Department, University of Warsaw, Żwirki i
Wigury 101, 02-089 Warszawa, Poland

f Université de Lorraine & CNRS, Laboratoire de Physique et Chimie Théoriques (LPCT), UMR CNRS
7019, 1 Boulevard Arago, 57078 Metz, France

## KEYWORDS

X-ray constrained wavefunction fitting, electron correlation, polarisation, exchange-correlation potential, quantum crystallography

<footnote>*Correspondence to Simon Grabowsky (simon.grabowsky@unibe.ch), †These authors contributed equally.</footnote>
## ABSTRACT

In X-ray constrained wavefunction (XCW) fitting, external information such as electron correlation and polarization are 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.

## 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.<sup>1,2</sup> 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”<sup>3,4</sup> which has been successfully used to calculate molecular exchange-correlation potentials.<sup>5</sup> This was the starting point for the development of the TH<sup>6,7,8</sup> and HCTH<sup>9,10,11</sup> functionals by Tozer, Handy 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,13,14,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,18,19 The XCW method consists in the minimization of the functional

$$L = E_{QM} + \lambda \cdot GooF^2 \quad (1)$$

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 ($|F_{model}|$) and measured ($|F_{exp}|$) X-ray structure factor amplitudes, which also takes into account standard uncertainties ($\sigma$) of the data:

$$GooF^2 = \frac{1}{N_r - N_p} \sum_{\vec{h}} \frac{(t|F_{model}(\vec{h})| - |F_{exp}(\vec{h})|)^2}{\sigma_{\vec{h}}^2}. \quad (2)$$

In the previous relation, $\vec{h}$ 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 put 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-centre functions using Hirshfeld's stockholder partitioning,20 whereby approximations to two-centre functions are also available.21,22,23

The mixing factor $\lambda$ 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 $\lambda$ is to terminate the fit, known as the “halting problem” of XCW fitting.20,24,25 The relationship between the ultimate $\lambda$ and the desired $GooF^2$ value is also in the centre 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 $\lambda$ 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,30,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.

I. *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.<sup>36,37,38,39,40,41,42,43</sup> 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.<sup>44,45,46,47</sup> Later, Chandler et al. investigated whether electron correlation can be extracted from simulated theoretical structure factors via XCW fitting.<sup>48</sup> 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, Macchi et al. found that "the effect of electron correlation on electron density can be partially captured by the XC-WF method".<sup>49</sup> 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.<sup>50</sup> 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.

## II. 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 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,<sup>51,52,53</sup> or using QM/MM methods.<sup>54,55</sup> 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.<sup>56,57</sup> This study was then extended to experimental data, still utilizing the multipole and the related invariance model.<sup>58</sup> 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.⁵⁹ This has been shown once more recently.⁶⁰

III. *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.⁶¹,⁶²,⁶³,⁶⁴

IV. *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.⁶⁵,⁶⁶ 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.⁶⁷ 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.⁵⁴ 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.⁶⁸ Subsequently, in a study using 14 different data sets 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 a XCW fitting.²⁵ 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. ²⁶.

Based on these findings and problems, we have chosen to use X-ray diffraction data sets
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.⁶⁹ and used, e.g., in refs. ⁷⁰,⁷¹,⁷², as well as
L-alanine measured by Destro et al. ⁷³ and used, e.g., in refs. ⁷⁴,⁷⁵,⁷⁶,⁷⁷,⁷⁸. 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 data sets 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 the 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.

## EXPERIMENTAL AND COMPUTATIONAL SECTION

*Clarification of concepts to facilitate separation of electron correlation and polarization*

i) 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 centred on the coordinates of the corresponding nuclei within the unit cell of a crystal.⁷⁹,⁸⁰,⁸¹
ii) 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.<sup>82,83,84</sup>

iii) 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.<sup>85</sup> Although there are more advanced techniques available (such as the use of ELMOs for the simulation of the environment),<sup>86</sup> these have not been tested in the XCW framework yet.

iv) The *reconstructed crystal* density is not a quantum mechanical approximation, as in i) – iii), 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<sup>87,88</sup>) 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.⁸⁹,⁹⁰

*Hirshfeld atom refinement of urea and L-alanine*

The crystallographic information files (CIFs) and structure factor lists (HKLs) of urea⁶⁹ and L-alanine⁷³ as available from the literature were used as input files for a Hirshfeld atom refinement (HAR)⁷⁰,⁸⁵ as implemented in the software *TONTO*⁹¹ using the restricted HF method and both the *def2-TZVP* and *pob-TZVP* (with the exact coefficients as implemented in *CRYSTAL14*⁹²) 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.⁹³,⁹⁴ Figure 1 shows the experimental structures after HAR at the RHF/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. ⁴⁹.

<figure><img src="image_2.png" /></figure>

**Figure 1.** (a) HAR-refined structure of urea at RHF/pob-TZVP level of theory. Crystallographic details as of ref. ⁶⁹: space group = *P* ¯42₁m (tetragonal), *a* = 5.5780(6) Å, *b* = 5.5780(6) Å, *c* = 4.6860(7) Å, *V* = 145.80(3) Å³. Measurement details as of ref. ⁶⁹: *T* = 123(2) K, sinθmax/λ = 1.4417 Å⁻¹. Refinement details from HAR: GooF² = 7.8748, *R* = 0.0172, wR = 0.0117, min/max residual density = -0.101/0.079 eÅ⁻³.

(b) HAR-refined structure of L-alanine at RHF/pob-TZVP level of theory. Crystallographic details as of ref. ⁷³: space group = *P*2₁2₁2₁ (orthorhombic), *a* = 5.9279(10) Å, *b* = 12.2597(17) Å, *c* = 5.7939(9) Å, *V* = 421.1(1) Å³. Measurement details as of ref. ⁷³: *T* = 23(1) K, sinθmax/λ = 1.0778 Å⁻¹. Refinement details from HAR: GooF² = 1.5029, *R* = 0.0198, wR = 0.0169, min/max residual density = -0.146/0.153 eÅ⁻³. Anisotropic displacement parameters of all atoms are drawn at 50% probability level.
<figure><img src="image_3.png" /></figure>

The software Gaussian09⁹⁵ 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 CCSD (no frozen core) and HF single-point calculations, respectively. The basis sets def2-TZVP and pob-TZVP (with the exact coefficients as implemented in CRYSTAL14⁹²) 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 CRYSTAL14 and was explicitly defined by us for the software programs TONTO and Gaussian09.

The software denprop⁹⁶ was employed to calculate the hkl indices and static structure factor magnitudes that include the effect of electron correlation. For this purpose, the Gaussian09 wavefunction output (.wfn file) of the CCSD calculation was used as input for a model pseudo-periodic system with cell constants of 10x10x10 Å³ (cubic cell) and a resolution of 0.7 Å⁻¹, 1.4417 Å⁻¹ and 2.0 Å⁻¹ for urea, and 20x20x20 Å³ (cubic cell) with a resolution of 0.7 Å⁻¹, 1.0778 Å⁻¹ and 2.0 Å⁻¹ 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⁹¹ was used for the XCW fitting against the theoretical structure factors. A .cif file was manually written with the coordinates from the Gaussian09 .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 Å⁻¹ for urea as well as 1.0778 Å⁻¹ and 2.0 Å⁻¹ for L-alanine) or 0.2 (0.7 Å⁻¹ 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 $\lambda$-values were obtained by monitoring until which $\lambda$-value the fitting converged within a given number of cycles, or when $\lambda = 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 $GooF^2 \sim \frac{1}{\sigma^2}$ (see Introduction), a 100 times bigger multiplier $\lambda$ is necessary to realize the same perturbation when $\sigma$ of 1.0 is used (as was done, for example, in Genoni, Macchi et al. ⁴⁹) instead of $\sigma = 0.1$ (as we do here). This means that with $\sigma = 0.1$, $\lambda = 10$ represents a high perturbation. It would equal $\lambda = 1000$ with $\sigma = 1.0$ as applied in ref. ⁴⁹. Other models of error distribution are possible and were recently discussed and tested in refs. 62,97

The determination and visualisation of electron correlation effects was done in terms of calculating the differences between electron density grid files of different $\lambda$-values and of the reference *Gaussian09* calculations in the Gaussian cube format. The reference electron correlation effect was calculated as 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 *Gaussian09* 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 $\lambda=x$ (TONTO_λ=x.cube) and the unconstrained wavefunction with $\lambda=0$ (TONTO_λ=0.cube, which is identical to Gaussian09_HF.cube). In these cases, the order of subtraction is always $\lambda=x$ minus $\lambda=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 towards the cube borders must be exactly the same in both cubes. With a second script, the absolute values of each grid point (rᵢ) of the difference density grids were summed up and then divided by two to obtain the total electron count of the difference, hence the total magnitude of the electron correlation effect as it impacts on the electron density (Nₑ = (Σₙₚᵢ=₁|ρₘethod₁(rᵢ) - ρₘethod₂(rᵢ)|)/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ᵢ) and the unconstrained values ρₜ=₀(rᵢ) in the TONTO_λ=₀.cube was divided by the reference electron density effect (CCSD minus HF). To obtain the real-space R-value (RSR)⁹⁸ (eq. 4), the XCW-fitted electron density values ρₜ=x(rᵢ) were directly referred to the CCSD calculation. Any difference grid files were visualized with the software VMD⁹⁹ as isosurface plots.

$$PR(\rho_{\lambda=x}, \rho_{\lambda=0}) = 100 \cdot \frac{\sum_{i=1}^{n_p} |\rho_{\lambda=x}(r_i) - \rho_{\lambda=0}(r_i)|}{\sum_{i=1}^{n_p} |\rho_{CCSD}(r_i) - \rho_{HF}(r_i)|} \quad (3)$$

$$RSR(\rho_{CCSD}, \rho_{\lambda=x}) = \frac{\sum_{i=1}^{n_p} |\rho_{CCSD}(r_i) - \rho_{\lambda=x}(r_i)|}{\sum_{i=1}^{n_p} |\rho_{CCSD}(r_i) + \rho_{\lambda=x}(r_i)|} \quad (4)$$

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

$$PR(\lambda) = c + \sum_{i=1}^{3} a_i \exp(-b_i \lambda) \quad (5)$$

was fitted to the observed course of PR versus λ 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 behaviour 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 minus λ=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.

## Polarization

The software *CRYSTAL14*¹⁰⁰ 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 3D 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 Å⁻¹ for urea as well as 0.7 Å⁻¹, 1.0778 Å⁻¹ and 2.0 Å⁻¹ 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 3D model with the *runprop14* utility as implemented in *CRYSTAL14*. 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 *CRYSTAL14*-obtained structure factors of the 3D 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 Å⁻¹ for urea), 0.2 (for resolutions of 0.7 Å⁻¹ for urea as well as 1.0778 Å⁻¹ and 2.0 Å⁻¹ for L-alanine) or 0.02 (0.7 Å⁻¹ for L-alanine).
applying a uniform structure factor error of 0.1. Again, fitting was stopped at the maximum value of $\lambda$ for which convergence was achieved or at $\lambda$=10. The same XCW fitting at HF/pob-TZVP level was also performed upon introducing surrounding cluster charges mimicking the crystal environment. By using this self-consistent cluster of point charges and dipoles, polarization effects were already simulated in the wavefunction ansatz at $\lambda$ = 0.

The reference polarization effect on the electron density (i.e., the interaction density) was calculated as difference electron density in 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 3D minus MOLSPLIT. CRYSTAL14 grid files were obtained with the runprop14 utility and the ECH3 keyword in the input (.d3) file. The XCW-fitted polarization was calculated as the difference electron density between the fitted wavefunction at the value $\lambda$=x (TONTO_λ=x.cube) and the unconstrained wavefunction with $\lambda$=0 (TONTO_λ=0.cube), for both cases with and without cluster charges in the ansatz. The order of subtraction is always $\lambda$=x minus $\lambda$=0. The distance between grid points obtained from CRYSTAL14 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:

$$PR(\rho_{\lambda=x}, \rho_{\lambda=0}) = 100 \cdot \frac{\sum_{i=1}^{np} |\rho_{\lambda=x}(r_i) - \rho_{\lambda=0}(r_i)|}{\sum_{i=1}^{np} |\rho_{3D}(r_i) - \rho_{MOLSPLIT}(r_i)|} \quad (6).$$

For visualisation of the CRYSTAL14 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⁵⁴ for this study. Subsequent visualization of isosurfaces was carried out with VMD.

### 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 Å⁻¹ for urea and 1.0778 Å⁻¹ for L-alanine as well as each pruned to the resolution of 0.7 Å⁻¹), four fittings were carried out: 1) HF with cluster charges (cc), 2) HF without cc, 3) BLYP with cc, 4) BLYP without cc. The effects represented by the four fittings always refer to difference electron densities from calculations at λ=x and λ=0.

## RESULTS AND DISCUSSION

### I. 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 isovalues are approximately around 0.03 eÅ⁻³, 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.³⁸,⁴³,⁴⁹,⁶⁷
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Plot of the difference between CCSD and HF electron densities used as reference for the reconstruction of electron correlation effects via XCW fitting for a) urea and b) L-alanine. Isosurface values in a.u. (solid/wireframe): +0.005/+0.0025 (blue) and -0.005/-0.0025 (red). 0.005 a.u. = 0.005 e bohr⁻³ = 0.034 e Å⁻³. Total number of electrons shifted: 0.27 (urea), 0.40 (L-alanine).</figcaption></figure>

Figures 3-5 and Table 1 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 the introduction, in refs. 48 and 49 the authors doubt that XCW fitting can recover it, at least not fully. Therefore, we first observe in Figure 3 a) to c) and Figure 4 a) to 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 Figure 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.
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0138312

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Plot of the difference in electron density for different $\lambda$ values during the XCW fitting of urea using CCSD structure factors of a resolution below 1.4415 Å$^{-1}$. Individual details are given in subcaptions. Isosurface values in a.u. (solid/wireframe): +0.005/+0.0025 (blue) and -0.005/-0.0025 (red). Representations for the other resolutions are given in the supplementary material.</figcaption></figure>

Figure 5 a) shows the curves of GooF$^2$, RSR and PR as a function of the increasing $\lambda$ values. GooF$^2$, a measure in reciprocal space, and RSR, a measure in real space, both approach zero with increasing $\lambda$, 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 versus $\lambda$ graph (eq. 5, discussed in the Experimental and Computational Section, gray line in Figure 5 a)) and checked the limiting behavior $\lim_{\lambda \to \infty} PR(\lambda) = c$ (see also Table 1). 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 $\lambda_{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.
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0138312

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Plot of the difference in electron density for different λ values during the XCW fitting on L-alanine using CCSD structure factors of a resolution below 1.0778 Å⁻¹. Individual details are given in subcaptions. Isosurface values in a.u. (solid/wireframe): +0.005/+0.0025 (blue) and -0.005/-0.0025 (red). Representations for the other resolutions are given in the supplementary material.</figcaption></figure>

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.¹⁰² 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, Figures 3 d) and 4 d) show that already at λ=0, before the fitting procedure starts, the DFT electron correlation map overestimates the
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0138312

reference CCSD map significantly. All features are more intense (comparison of Figures 3 d) with 2 a) and 4 d) with 2 b)) and the PR values are 156.5% (urea) and 159.4% (L-alanine), respectively.

<figure><img src="image_7.png" /></figure>

**Figure 5.** Graphs of structure factor and electron density agreement throughout XCW fittings a) for urea at a maximum resolution of 1.4415 Å⁻¹ and (b) L-alanine at 1.0778 Å⁻¹. GooF² (black, squares), RSR (blue, circles) and PR (red, triangles) plotted against the value of λ during the HF-XCW(solid) or BLYP-XCW(dashed) fitting. Additionally, exponential functions were fitted to the values to estimate the convergence limit of the PR (gray) (eq. 5). Graphs for the other resolutions are in the supplementary material.

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 \to \infty} PR(\lambda) = c$ using eq. 5 (gray lines in Figure 5 b)) leads to limiting values of 108.2 and 109.8% (Table 1), showing again that the vast majority of the electron correlation effect has already been included up to $\lambda_{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 towards a better description of electron correlation. How this manifests itself in the exchange-correlation (XC) potential will be shown below.

<table><caption>**Table 1.** Leading parameters $a_1$ and $b_1$ and limiting parameter $c (= PR_\infty)$ of the exponential function (see eq. 5) fitting the percentage reconstruction (PR) of the electron correlation density as a function of $\lambda$ (see Figure 5). Maximum $\lambda$ value for a converged XCW fitting, or alternatively $\lambda=10$, and the related $PR_{max}$ (in %), $RSR_{max}$ and $GooF^2_{max}$ values. First part urea, second part L-alanine, for all resolutions. The remaining values are in Tables S2 and S6 in the supplementary material.</caption><thead><tr><th>Resolution /Å<sup>-1</sup></th><th>Ansatz</th><th>a<sub>1</sub></th><th>b<sub>1</sub></th><th>c (= PR<sub>∞</sub>)</th><th>PR<sub>max</sub></th><th>RSR<sub>max</sub></th><th>GooF<sup>2</sup><sub>max</sub></th><th>λ<sub>max</sub></th></tr></thead><tbody><tr><td>0.7</td><td rowspan="3">HF<br>(urea)</td><td>-12.51</td><td>4.47</td><td>100.5</td><td>94.3</td><td>0.001240</td><td>0.000302</td><td>3.2</td></tr><tr><td>1.4415</td><td>-34.14</td><td>5.49</td><td>93.9</td><td>88.4</td><td>0.001460</td><td>0.000179</td><td>10.0</td></tr><tr><td>2.0</td><td>-49.35</td><td>8.07</td><td>88.6</td><td>74.1</td><td>0.002443</td><td>0.000302</td><td>10.0</td></tr><tr><td>0.7</td><td rowspan="3">BLYP<br>(urea)</td><td>15.79</td><td>0.43</td><td>109.1</td><td>109.6</td><td>0.001632</td><td>0.000686</td><td>2.0</td></tr><tr><td>1.4415</td><td>15.00</td><td>4.74</td><td>108.2</td><td>110.2</td><td>0.001691</td><td>0.000135</td><td>9.5</td></tr><tr><td>2.0</td><td>20.78</td><td>6.80</td><td>110.3</td><td>115.1</td><td>0.002141</td><td>0.000201</td><td>10.0</td></tr><tr><td>0.7</td><td rowspan="3">HF<br>(L-alanine)</td><td>-35.33</td><td>0.26</td><td>94.9</td><td>92.8</td><td>0.001257</td><td>0.000699</td><td>2.6</td></tr><tr><td>1.0778</td><td>-38.47</td><td>0.09</td><td>96.0</td><td>94.1</td><td>0.001127</td><td>0.000186</td><td>10.0</td></tr><tr><td>2.0</td><td>-32.56</td><td>1.77</td><td>86.4</td><td>74.1</td><td>0.002500</td><td>0.000448</td><td>10.0</td></tr><tr><td>0.7</td><td rowspan="3">BLYP<br>(L-alanine)</td><td>16.77</td><td>0.43</td><td>109.1</td><td>109.9</td><td>0.001773</td><td>0.001142</td><td>1.8</td></tr><tr><td>1.0778</td><td>16.53</td><td>0.80</td><td>109.8</td><td>110.0</td><td>0.001736</td><td>0.000392</td><td>4.5</td></tr><tr><td>2.0</td><td>15.99</td><td>4.15</td><td>113.7</td><td>116.5</td><td>0.002301</td><td>0.000313</td><td>10.0</td></tr></tbody></table>

Before we turn to XC potentials, we will discuss the resolution dependence of the results with the help of Table 1. Intuitively, one might think that reconstructions ought to be more
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0138312

successful if higher resolution data are used because truncation of resolution means information loss in general. And indeed, the $\lambda_{max}$ values in Table 1 indicate that with less data XCW fitting becomes more problematic, as at the lowest resolution of 0.7 Å⁻¹ the fitting never converges above $\lambda_{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 \to \infty} PR(\lambda) = PR_{\infty} = 100\%$ can be reached, and, in contrast, for the high resolution of 2.0 Å⁻¹, $PR_{\infty}$ and $PR_{max}$ are only 88.6 and 74.1%, respectively. The same trend is present in the RSR values. And 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 $\lambda$ at higher resolution. This means that much higher values of $\lambda$ 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 towards 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. ⁴⁹). It is worth noting that the GooF² values do not describe the same observation; they are always significantly larger at low resolution (Table 1).

## I.1 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.⁸⁷,⁸⁸ Since the XC potential is not confined to have a
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0138312

fixed integral over a certain volume of space, the calculation of the RSR or PR values would be
on an arbitrary and ambiguous scale. Further, 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.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Plot of the difference in the exchange-correlation potential in form of the BLYP functional for different $\lambda$ values during the XCW fitting of urea using CCSD structure factors of a maximum resolution of 1.4415 Å$^{-1}$. Individual details are given in subcaptions. Isosurface values in a.u. (solid/wireframe): +0.007/+0.005 (blue) and -0.007/-0.005 (red). Plots at the other two resolutions are given in the supplementary material.</figcaption></figure>

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.<sup>103,104</sup>
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 $\lambda$ should be physically meaningful. Table S1 shows that at $\lambda$=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 Figure 6 a) from $\lambda$=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 $\lambda$ (Figures 6 b) and 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 (Figures 6 d) to 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 Figures 3 d) to 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 3 d) to f)) and too high electron density in the valence regions (extension and intensity of negative (red) isosurfaces reduce from 3 d) to 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 Figures 6 d) to 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 Figure 6 to Figures S30 and S32 in the supplementary material): Features in the XC potential are much more pronounced already at small $\lambda$ values when only low-
resolution data are included in the fitting. In other words, to reproduce the same features qualitatively and quantitatively much higher $\lambda$ values are needed for higher-resolution structure factor data sets.

## II. 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 Figure 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 Figure 7 a), left functional group in Figure 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.

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0138312

**Figure 7.** Plot of the difference between HF electron densities under periodic-boundary and isolated-molecule conditions (interaction density) used as reference for the reconstruction of polarization effects via XCW fitting for a) urea and b) L-alanine. Isosurface values in a.u. (solid/wireframe): +0.0025/+0.00125 (blue) and -0.0025/-0.00125 (red). 0.0025 a.u. = 0.0025 e bohr⁻³ = 0.017 e Å⁻³. Total number of shifted electrons: 0.25 (urea), 0.44 (L-alanine). Plots of the difference between BLYP electron densities under periodic-boundary and isolated-molecule conditions look very similar (Figures S26 and S59).

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.⁵⁴,⁵⁸,⁵⁹,⁶⁰,¹⁰⁵

<figure><img src="image_11.png" /></figure>
**Figure 8.** Plots of the difference in electron density for different $\lambda$ values during the XCW fitting of urea using structure factors from periodic-boundary calculations of a resolution below 1.4415 Å$^{-1}$. Individual details are given in subcaptions. Isosurface values in a.u. (solid/wireframe): +0.0025/+0.00125 (blue) and -0.0025/-0.00125 (red). Representations for the other resolutions are given in the supplementary material. CC = cluster charges.

Figures 8-10 and Table 2 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,⁶⁰ but electron correlation only partially.⁴⁹ By comparing the build-up of the qualitative features in Figures 8 a) to c) and 9 a) to c) with the features in the corresponding reference maps in Figure 7, it becomes clear that XCW fitting indeed reconstructs the crystal effect accurately. A comparison of the polarization with the electron correlation effect (Figures 8 a) to c) with 3 a) to c) and 9 a) to c) with 4 a) to c)) shows that polarization is reconstructed more quickly, i.e. already at low $\lambda$-values a significant portion of the polarization effect has been recovered, more than for electron correlation (at $\lambda$=1.0: 62.3% vs. 44.1% for urea and 61.2% (at $\lambda$=0.2) vs 48.9% (at $\lambda$=0.5) for L-alanine).

Nevertheless, the initially steep slope for PR vs. $\lambda$ flattens so that PR converges for polarization towards a value lower than 100% (Figure 10), similar to electron correlation (Figure 5). In fact, Table 2 shows that the limiting PR value from a fit of PR vs. $\lambda$ 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 1. This means that both effects can be reconstructed by XCW fitting as accurately and efficiently, only polarization builds up at low values of $\lambda$ more quickly. The consequence of this is, in turn, that electron density maps generated from XCW fitting to
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0138312

experimental data, as shown in the next section, should be dominated by polarization at low values of $\lambda$ and electron correlation at high values of $\lambda$.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Plots of the difference in electron density for different $\lambda$ values during the XCW fitting of L-alanine using structure factors from periodic-boundary calculations of a resolution below 1.0778 Å⁻¹. Individual details are given in subcaptions. Isosurface values in a.u. (solid/wireframe): +0.0025/+0.00125 (blue) and -0.0025/-0.00125 (red). Representations for the other resolutions are given in the supplementary material. CC = cluster charges.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Graphs of structure factor and electron density agreement throughout XCWs at a resolution below 1.4415 Å⁻¹ for urea (a) and 1.0778 Å⁻¹ for L-alanine (b). GooF² (black) and the PR (red) plotted against the value of λ during the HF-XCW (solid) or HF-CC-XCW (dashed). Additionally exponential functions were fitted to the values to estimate the convergence limit of the PR (gray) according to eq. 5, except for HF-CC in a) as the function is too horizontal. CC = cluster charges.</figcaption></figure>

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 towards the reference values.
The analogue of this theoretical experiment for polarization is the inclusion of a cluster of point
charges and dipoles at symmetry-generated positions around the central atom to simulate the
crystal effect (abbreviated CC here for “cluster charges”).⁸⁵,⁸⁶ As for electron correlation, the effect is clearly overestimated in the cluster charges, see Figures 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 Figures 8 e) and f) as wells as 9 e) and 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 section III. 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.⁸⁶

Interestingly, this behavior cannot be seen in the GooF² values that simply follow the λ-values. The higher λ, the smaller GooF² (Figures 8 d) to f) and 9 d) to f), Tables S4 and S7). This was also seen in the discussion of Table 1, and it again shows in Table 2 where overall higher λmax values correspond to lower GooF²max values. This means that a discussion of the halting problem must not only be based on the GooF² vs. λ relationship,²⁰,²⁴,²⁵ but it should involve derived chemical properties. Here, one could draw a parallel to the variational principle of quantum chemistry: just minimising 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⁵⁹,⁶⁰ agrees with the trends observed for electron correlation, but is less pronounced (see Table 2). 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 $\lambda$. This means, on the other hand, that convergence ceases much earlier for low-resolution data sets. The resolution dependence of the maximum achievable PR value on resolution is clearly less pronounced than for electron correlation (compare Table 2 with Table 1).

<table><caption>**Table 2.** Leading parameters $a_1$ and $b_1$ and limiting parameter $c (= PR_\infty)$ of the exponential function (see eq. 5) fitting the percentage reconstruction (PR) of the interaction density as a function of $\lambda$ (see Figure 10). Maximum $\lambda$ value for a converged XCW fitting, or alternatively $\lambda=10$, and the related $PR_{max}$ (in %) and $GooF^2_{max}$ values. First part urea, second part L-alanine, for all resolutions. N/A...fit failed since the slope was nearly constant.</caption><thead><tr><th>Resolution /Å<sup>-1</sup></th><th>Ansatz</th><th>a<sub>1</sub></th><th>b<sub>1</sub></th><th>c (= PR<sub>∞</sub>)</th><th>PR<sub>max</sub></th><th>GooF<sub>max</sub><sup>2</sup></th><th>λ<sub>max</sub></th></tr></thead><tbody><tr><td>0.7</td><td rowspan="3">HF<br>(urea)</td><td>N/A</td><td>N/A</td><td>N/A</td><td>94.0</td><td>0.000771</td><td>1.6</td></tr><tr><td>1.4415</td><td>-32.23</td><td>0.97</td><td>96.4</td><td>92.1</td><td>0.000150</td><td>9.5</td></tr><tr><td>2.0</td><td>-33.99</td><td>1.24</td><td>90.7</td><td>82.8</td><td>0.000177</td><td>10.0</td></tr><tr><td>0.7</td><td rowspan="3">HF_CC<br>(urea)</td><td>N/A</td><td>N/A</td><td>N/A</td><td>116.2</td><td>0.000294</td><td>1.6</td></tr><tr><td>1.4415</td><td>N/A</td><td>N/A</td><td>N/A</td><td>116.2</td><td>0.000062</td><td>9.5</td></tr><tr><td>2.0</td><td>N/A</td><td>N/A</td><td>N/A</td><td>117.1</td><td>0.000052</td><td>10.0</td></tr><tr><td>0.7</td><td rowspan="3">HF<br>(L-alanine)</td><td>-34.74</td><td>0.05</td><td>96.6</td><td>89.9</td><td>0.007940</td><td>0.4</td></tr><tr><td>1.0778</td><td>-29.52</td><td>0.26</td><td>88.6</td><td>86.7</td><td>0.003566</td><td>1.0</td></tr><tr><td>2.0</td><td>-35.17</td><td>0.93</td><td>95.9</td><td>88.9</td><td>0.000498</td><td>7.4</td></tr><tr><td>0.7</td><td rowspan="3">HF_CC<br>(L-alanine)</td><td>N/A</td><td>N/A</td><td>N/A</td><td>120.2</td><td>0.003277</td><td>0.46</td></tr><tr><td>1.0778</td><td>3.02</td><td>0.01</td><td>119.0</td><td>120.1</td><td>0.001114</td><td>1.6</td></tr><tr><td>2.0</td><td>N/A</td><td>N/A</td><td>N/A</td><td>120.9</td><td>0.000271</td><td>6.8</td></tr></tbody></table>

### III. 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 $\lambda_0$ and $\lambda_{max}$ in the following way: $\Delta\rho = \rho_{\lambda_{max,HF}} - \rho_{\lambda_{0,HF}} = \Delta\rho_{EC} + \Delta\rho_{pol} + \Delta\rho_{defect}$. As discussed in the section “Clarification of concepts to facilitate separation of electron correlation and polarization” in detail, it is unclear whether these effects can be separated from each other since
they are convolved in the experimental data. Also, sections I. and II. have shown that there will be residual effects left in the fitted wavefunction. Therefore, we will make the following idealised 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: $\Delta\rho_{pol}$, **polarization effects only** (correlation is included at $\lambda = 0$, but not polarization);

- HF with cluster charges: $\Delta\rho_{EC}$, **correlation effects only** (polarization is included at $\lambda = 0$, but not correlation);

- BLYP with cluster charges: $\Delta\rho_{defect}$, **defect density** (correlation and polarization are included at $\lambda = 0$);

- HF without cluster charges: $\Delta\rho$, **all effects** (nothing is included at $\lambda = 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 Figure 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 fitting. This is also reflected in the total number of electrons shifted during the fitting ($N_e$), which is 0.32/0.52e (urea/L-alanine), and it was 0.25/0.44e in the reference calculations (Figure 7). We therefore conclude that the chosen approach indeed qualitatively extracts $\Delta\rho_{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 Figure 2. Again, there is qualitative
similarity, and it is noteworthy that the observations in Figures 11 and 12 versus 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<sub>e</sub>) is 0.28/0.42e (urea/L-
alanine), and it was 0.27/0.40e in the reference calculations (Figure 2). Again, we conclude that
the chosen approach indeed qualitatively extracts Δρ<sub>EC</sub> 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 persists for the exchange-correlation potential, when comparing Figure 13 a) with Figure
6 a) to c) and Figure 13 b) with Figure 6 d) to f).
<figure><img src="image_14.png" /></figure>

a) $\lambda_{0.022,BLYP} - \lambda_{0.0,BLYP}$, represents **polarization**; $N_e=0.32e$

b) $\lambda_{0.046,HF\_CC} - \lambda_{0.0,HF\_CC}$, represents **correlation**; $N_e=0.28e$

<figure><img src="image_15.png" /></figure>

c) $\lambda_{0.046,HF} - \lambda_{0.0,HF}$, represents **all effects**; $N_e=0.35e$

d) $\lambda_{0.022,BLYP\_CC} - \lambda_{0.0,BLYP\_CC}$, represents **defect density**; $N_e=0.25e$

**Figure 11.** Plots of the difference in electron density during XCW fittings of urea using different effects in the *ansatz* based on experimental structure factors of a resolution below 1.4415 Å⁻¹. Individual details are given in subcaptions. Isosurface values in a.u. (solid/wireframe): +0.0025/+0.00125 (blue) and -0.0025/-0.00125 (red). Representations for the resolution of 0.7 Å⁻¹ are given in the supplementary material. CC = cluster charges.
<figure><img src="image_16.png" /></figure>

a) $\lambda_{0.34,BLYP} - \lambda_{0.0,BLYP}$, represents **polarization**; N<sub>e</sub>=0.52e

b) $\lambda_{0.53,HF\_CC} - \lambda_{0.0,HF\_CC}$, represents **correlation**; N<sub>e</sub>=0.42e

<figure><img src="image_17.png" /></figure>

c) $\lambda_{0.53,HF} - \lambda_{0.0,HF}$, represents **all effects**; N<sub>e</sub>=0.52

d) $\lambda_{0.31,BLYP\_CC} - \lambda_{0.0,BLYP\_CC}$, represents **defect density**; N<sub>e</sub>=0.43

**Figure 12.** Plots of the difference in electron density during XCW fittings of L-alanine using different effects in the *ansatz* based on experimental structure factors of a resolution below 1.0778 Å⁻¹. Individual details are given in subcaptions. Isosurface values in a.u. (solid/wireframe): +0.0025/+0.00125 (blue) and -0.0025/-0.00125 (red). Representations for the resolution of 0.7 Å⁻¹ are given in the supplementary material. CC = cluster charges.
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0138312

<figure><img src="image_18.png" /><figcaption>**Figure 13.** Plot of the difference in the exchange-correlation potential in form of the BLYP functional during the XCW fitting of urea using experimental structure factors of a maximum resolution of 1.4415 Å⁻¹ and using cluster charges. Individual details are given in subcaptions. Isosurface values in a.u. (solid/wireframe): +0.007/+0.0055 (blue) and -0.007/-0.0055 (red). Plots at other values of λ and at the resolution of 0.7 Å⁻¹ are given in the supplementary material.</figcaption></figure>

Figures 11 c) and 12 c) show the merger of all effects at $\lambda_{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 (Figures 11 b) and 12 b)) much more closely than those of polarization (Figures 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 $\lambda_{max}$ electron correlation dominates over polarization; the combined effect is not a combination of the features of both as one might have expected. Comparing Figures 11 b) with c) and 12 b) with c) shows, though, that the effects are not identical, there are more pronounced and extended features in the maps 11/12 c) representing all effects. The N<sub>e</sub> values are also higher (0.35 vs. 0.28e for urea; 0.52 vs. 0.42e 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 section II. 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 λ<sub>max</sub> presented here in Figure 11 and
12 c). Therefore, in the supplementary material (Figures 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 section II., that electron correlation dominates the plots in all ranges of
λ, even at very small values.

Finally, we discuss the defect densities shown in Figures 11 d) and 12 d). Δρ<sub>defect</sub> 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.<sup>25</sup> 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. The previous two sections 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 I. has also shown that XCW fitting mostly corrects this overestimation in BLYP, whereas section II. 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 Figures 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 section II, 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, $\Delta\rho_{defect}$ really represents the defect in the procedure and not a measure of
the experimental random and systematic error.

CONCLUSIONS

*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 $\lambda$, than electron correlation. Nevertheless, if fitted to $\lambda_{max}$, both polarization and electron correlation effects can be recovered to the same extent and with the same accuracy and reliability. An extrapolation to $\lambda_{\infty}$ 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
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0138312

maximum resolution of the structure factors set; a medium-resolution data sets allows faster
and more complete reconstruction than an ultra-high resolution data set.

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.

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. And 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.
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0138312

The results in this paper are based on two data sets of high quality, extensively used as test sets for benchmarking new methods in the field of quantum crystallography for decades. In a different study,²⁵ 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.05e 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 towards general use of experimental diffraction data in method development for quantum chemistry, a reproducibility study using very high-quality data sets 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 GR 4451/1-1 (S.G.), for a postdoctoral fellowship HU 2512/1-1 (E.H.), and for a Walter-Benjamin fellowship KL 3500/1-1 (F.K.). M Woińska acknowledges financial support from the Polish National Science Centre (NCN) under Opus grant DEC-2018/31/B/ST4/02142.

## AUTHOR DECLARATIONS

**Conflict of Interest.** The authors have no conflicts to disclose.

**Author Contributions.** Emanuel Hupf: Conceptualization (supporting), Formal Analysis (equal), Funding Acquisition (supporting), Investigation (equal), Methodology (equal), Software (equal), Visualization (equal), Writing/Original Draft Preparation (supporting); Florian Kleemiss: Conceptualization (supporting), Formal Analysis (equal), Funding Acquisition (supporting), Investigation (equal), Methodology (equal), Software (equal), Visualization (equal), Writing/Original Draft Preparation (supporting); Tobias Borrmann:
<figure><img src="image_19.png" /></figure>

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 Preparation (lead), Writing/Review & Editing (supporting).

DATA AVAILABILITY

The data that supports 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

<footnote n="1">¹ K. Eichkorn, O. Treutler, H. Öhm, M. Häser, and R. Ahlrichs, “Auxiliary basis sets to approximate Coulomb potentials,” *Chem. Phys. Lett.* **240**, 283–290 (1995).</footnote>

<footnote n="2">² M. Franchini, P. H. T. Philipsen, E. van Lenthe, and L. Visscher, “Accurate Coulomb potentials for periodic and molecular systems through density fitting,” *J. Chem. Theory Comput.* **10**, 1994–2004 (2014).</footnote>

<footnote n="3">³ Q. Zhao and R. G. Parr, “Constrained-search method to determine electronic wave functions from electronic densities,” *J. Chem. Phys.* **98**, 543–548 (1993).</footnote>
<figure><img src="image_20.png" /></figure>

<footnote n="4">⁴ Q. Zhao, R. C. Morrison, and R. G. Parr, “From electron densities to Kohn-Sham kinetic energies, orbital energies, exchange-correlation potentials, and exchange-correlation energies,” *Phys. Rev. A* **50**, 2138–2142 (1994).</footnote>

<footnote n="5">⁵ D. J. Tozer, V. E. Ingamells, and N. C. Handy, “Exchange-correlation potentials,” *J. Chem. Phys.* **105**, 9200–9213 (1996).</footnote>

<footnote n="6">⁶ D. J. Tozer and N. C. Handy, “The development of new exchange-correlation functionals,” *J. Chem. Phys.* **108**, 2545–2555 (1998).</footnote>

<footnote n="7">⁷ D. J. Tozer and N. C. Handy, “Development of new exchange-correlation functionals. 2,” *J. Phys. Chem. A* **102**, 3162–3168 (1998).</footnote>

<footnote n="8">⁸ D. J. Tozer and N. C. Handy, “The development of new exchange-correlation functionals: 3,” *Mol. Phys.* **94**, 707–715 (1998).</footnote>

<footnote n="9">⁹ F. A. Hamprecht, A. J. Cohen, D. J. Tozer, and N. C. Handy, “Development and assessment of new exchange-correlation functionals,” *J. Chem. Phys.* **109**, 6264–6271 (1998).</footnote>

<footnote n="10">¹⁰ A. D. Boese, N. L. Doltsinis, N. C. Handy, and M. Sprik, “New generalized gradient approximation functionals,” *J. Chem. Phys. B* **112**, 1670–1678 (2000).</footnote>

<footnote n="11">¹¹ A. D. Boese and N. C. Handy, “New exchange-correlation density functionals: The role of the kinetic energy density,” *J. Chem. Phys.* **116**, 9559–9569 (2002).</footnote>

<footnote n="12">¹² S. Grabowsky, A. Genoni, and H.-B. Bürgi, “Quantum crystallography,” *Chem. Sci.* **8**, 4159–4176 (2017).</footnote>

<footnote n="13">¹³ A. Genoni, L. Bučinský, N. Claiser, J. Contreras-García, B. Dittrich, P. M. Dominiak, E. Espinosa, C. Gatti, P. Giannozzi, J.-M. Gillet, D. Jayatilaka, P. Macchi, A. Ø. Madsen, L. Massa, C. F. Matta, K. M. Merz, P. N. H. Nakashima, H. Ott, U. Ryde, K. Schwarz, M. Sierka, and S. Grabowsky, “Quantum crystallography: Current developments and future perspectives,” *Chem. Eur. J.* **24**, 10881–10905 (2018).</footnote>

<footnote n="14">¹⁴ A. Genoni and P. Macchi, “Quantum crystallography in the last decade: Developments and outlooks,” *Crystals* **10**, 473: 1–20 (2020).</footnote>

<footnote n="15">¹⁵ S. Grabowsky, A. Genoni, S. P. Thomas, and D. Jayatilaka, “The advent of quantum crystallography: form and structure factors from quantum mechanics for advanced structure refinement and wavefunction fitting,” *Structure and Bonding*, *21st Century Challenges in Chemical Crystallography II*, pp. 65–144 (2020).</footnote>

<footnote n="16">¹⁶ W. L. Clinton and L. J. Massa, “Determination of the electron density matrix from X-ray diffraction data,” *Phys. Rev. Lett.* **29**, 1363–1366 (1972).</footnote>

<footnote n="17">¹⁷ D. Jayatilaka, “Wave function for beryllium from X-ray diffraction data,” *Phys. Rev. Lett.* **80**, 798–801 (1998).</footnote>
¹⁸ D. Jayatilaka and D. J. Grimwood, “Wavefunctions derived from experiment. I. Motivation and theory,” *Acta Cryst.* **A** *57*, 76–86 (2001).

¹⁹ D. J. Grimwood and D. Jayatilaka, “Wavefunctions derived from experiment. II. A wavefunction for oxalic acid dihydrate,” *Acta Cryst.* **A** *57*, 87–100 (2001).

²⁰ M. L. Davidson, S. Grabowsky, and D. Jayatilaka, “X-ray constrained wavefunctions based on Hirshfeld atoms. I. Method and review,” *Acta Cryst.* **B** *78*, 312–332 (2022).

²¹ R. F. Stewart, “Generalized X-ray scattering factors,” *J. Chem. Phys.* **51**, 4569–4577 (1969).

²² P. Coppens, T. V. Willoughby, and L. N. Csonka, “Electron population analysis of accurate diffraction data. I. Formalisms and restrictions,” *Acta Cryst.* **A** *27*, 248–256 (1971).

²³ K. Tanaka, “X-ray analysis of wavefunctions by the least-squares method incorporating orthonormality. I. General formalism,” *Acta Cryst.* **A** *44*, 1002–1008 (1988).

²⁴ A. Genoni, “On the termination of the X-ray constrained wavefunction procedure: reformulation of the method for an unequivocal determination of $\lambda$,” *Acta Cryst.* **A** *78*, 302–308 (2022).

²⁵ M. L. Davidson, S. Grabowsky and D. Jayatilaka, “X-ray constrained wavefunctions based on Hirshfeld atoms. II. Reproducibility of electron densities in crystals of $\alpha$-oxalic acid dihydrate,” *Acta Cryst.* **B** *78*, 397–415 (2022).

²⁶ H.-B. Bürgi and A. Genoni, “Remarks on X-ray constrained/restrained wavefunction fitting,” *Acta Cryst.* **B** *78*, 298–304 (2022).

²⁷ A. E. Whitten, D. Jayatilaka, and M. A. Spackman, “Effective molecular polarizabilities and crystal refractive indices estimated from x-ray diffraction data,” *J. Chem. Phys.* **125**, 174505:1–14 (2006).

²⁸ D. D. Hickstein, J. M. Cole, M. J. Turner, and D. Jayatilaka, “Modeling electron density distributions from X-ray diffraction to derive optical properties: Constrained wavefunction versus multipole refinement,” *J. Chem. Phys.* **139**, 064108:1–14 (2013).

²⁹ S. Grabowsky, D. Jayatilaka, S. Mebs, and P. Luger, “The electron localizability indicator from X-ray diffraction data—A first application to a series of epoxide derivatives,” *Chem. Eur. J.* **16**, 12818–12821 (2010).

³⁰ S. Grabowsky, M. Weber, D. Jayatilaka, Y.-S. Chen, M. T. Grabowski, R. Brehme, M. Hesse, T. Schirmeister, and P. Luger, “Reactivity differences between $\alpha,\beta$-unsaturated carbonyls and hydrazones investigated by experimental and theoretical electron density and electron localizability analyses,” *J. Phys. Chem. A* **115**, 12715–12732 (2011).
<footnote n="31">³¹ S. Grabowsky, P. Luger, J. Buschmann, T. Schneider, T. Schirmeister, A. N. Sobolev, and D. Jayatilaka, “The significance of ionic bonding in sulfur dioxide: Bond orders from X-ray diffraction data. *Angew. Chem. Int. Ed.* **51**, 6776–6779 (2012).</footnote>

<footnote n="32">³² A. Genoni, “Molecular orbitals strictly localized on small molecular fragments from X-ray diffraction data,” *J. Phys. Chem. Lett.* **4**, 1093–1099 (2013).</footnote>

<footnote n="33">³³ A. Genoni, “X-ray constrained extremely localized molecular orbitals: theory and critical assessment of the new technique,” *J. Chem. Theory Comput.* **9**, 3004–3019 (2013).</footnote>

<footnote n="34">³⁴ A. Genoni, “A first-prototype multi-determinant X-ray constrained wavefunction approach: the X-ray constrained extremely localized molecular orbital–valence bond method,” *Acta Cryst.* **A** **73**, 312–316 (2017).</footnote>

<footnote n="35">³⁵ A. Genoni, D. Franchini, S. Pieraccini, and M. Sironi, “X-ray constrained spin-coupled wavefunction: a new tool to extract chemical information from X-ray diffraction data,” *Chem. Eur. J.* **24**, 15507–15511 (2018).</footnote>

<footnote n="36">³⁶ R. F. W. Bader, and A. K. Chandra, “A view of bond formation in terms of molecular charge distributions,” *Can. J. Chem.* **46**, 953–966 (1968).</footnote>

<footnote n="37">³⁷ V. H. Smith, “Theoretical determination and analysis of electronic charge distributions,” *Phys. Scr.* **15**, 147–162 (1977).</footnote>

<footnote n="38">³⁸ M. E. Stephens and P. J. Becker, “Virial partitioning analysis of electron correlation and nuclear motion in diatomic molecules,” *Mol. Phys.* **49**, 65–89 (1983).</footnote>

<footnote n="39">³⁹ J. P. Ritchie, H. F. King, and W. S. Young, “Structures and energies for C₄,” *J. Chem. Phys.* **85**, 5175–5182 (1986).</footnote>

<footnote n="40">⁴⁰ R. Moszynski, and K. Szalewicz, “Accurate electron densities of the hydrogen molecule,” *J. Phys. B* **20**, 4347 (1987).</footnote>

<footnote n="41">⁴¹ C. Gatti, P. J. McDougall, and R. F. W. Bader, “Effect of electron correlation on the topological properties of molecular charge distributions,” *J. Chem. Phys.* **88**, 3792–3804 (1988).</footnote>

<footnote n="42">⁴² R. J. Boyd and L. C. Wang, “The effect of electron correlation on the topological and atomic properties of the electron density distributions of molecules,” *J. Comput. Chem.* **10**, 367–375 (1989).</footnote>

<footnote n="43">⁴³ K. B. Wiberg, C. M. Hadad, T. J. LePage, C. M. Breneman, and M. J. Frisch, “Analysis of the effect of electron correlation on charge density distributions,” *J. Phys. Chem.* **96**, 671–679 (1992).</footnote>

<footnote n="44">⁴⁴ L. S. Bartell and R. M. Gavin, “Effects of electron correlation in X-ray and electron diffraction,” *J. Am. Chem. Soc.* **86**, 3493–3498 (1964).</footnote>

<footnote n="45">⁴⁵ L. S. Bartell and R. M. Gavin, “Effects of electron correlation in X-ray and electron diffraction. II. Influence of</footnote>
nuclear charge in two-electron systems,” *J. Chem. Phys.* **43**, 856–861 (1965).

46 R. M. Gavin and L. S. Bartell, “Effects of electron correlation in X-ray and electron diffraction. III. Experimental electron–electron distribution functions,” *J. Chem. Phys.* **44**, 3687–3691 (1966).

47 R. M. Gavin and L. S. Bartell, “Effects of electron correlation in X-ray and electron diffraction. IV. Approximate treatment for many-electron atoms,” *J. Chem. Phys.* **45**, 4700–4703. (1966)

48 I. Bytheway, G. Chandler, B. Figgis, and D. Jayatilaka, “Are intramolecular dynamic electron correlation effects detectable in X-ray diffraction experiments on molecular crystals?” *Acta Cryst.* **A** **63**, 135–145 (2007).

49 A. Genoni, L. H. R. dos Santos, B. Meyer and P. Macchi, “Can X-ray constrained Hartree–Fock wavefunctions retrieve electron correlation?” *IUCR J* **4**, 136–146 (2017).

50 D. Ramírez-Palma, B. Landeros-Rivera, A. Genoni, F. Cortés-Guzmán, and J. Contreras-García, Visualizing correlation regions: The case of the ammonia crystal. *Chem. Methods* **2**, e202100045:1–8 (2022).

51 C. Gatti, V. R. Saunders, and C. Roetti, “Crystal field effects on the topological properties of the electron density in molecular crystals: The case of urea,” *J. Chem. Phys.* **101**, 10686–10696 (1994).

52 C. Gatti, B. Silvi, and F. Colonna, “Dipole moment of the water molecule in the condensed phase: a periodic Hartree-Fock estimate,” *Chem. Phys. Lett.* **247**, 135–141 (1995).

53 J. A. Platts and S. T. Howard, “Periodic Hartree–Fock calculations on crystalline HCN,” *J. Chem. Phys.* **105**, 4668–4674 (1996).

54 F. Kleemiss, E. K. Wieduwilt, E. Hupf, M. W. Shi, S. G. Stewart, D. Jayatilaka, M. J. Turner, K. Sugimoto, E. Nishibori, T. Schirmeister, T. C. Schmidt, B. Engels, and S. Grabowsky, “Similarities and differences between crystal and enzyme environmental effects on the electron density of drug molecules,” *Chem. Eur. J.* **27**, 3407–3419 (2021).

55 F. Kleemiss, P. Puylaert, D. Duvinage, M. Fugel, K. Sugimoto, J. Beckmann, and S. Grabowsky, “Ibuprofen and sila-ibuprofen: polarization effects in the crystal and enzyme environments,” *Acta Cryst.* **B** **77**, 892–905 (2021).

56 R. Y. de Vries, D. Feil, and V. G. Tsirelson, “Extracting charge density distributions from diffraction data: a model study on urea,” *Acta Cryst.* **B** **56**, 118–123 (2000).

57 M. A. Spackman, P. G. Byrom, M. Alfredsson, and K. Hermansson, “Influence of intermolecular interactions on multipole-refined electron densities,” *Acta Cryst.* **A** **55**, 30–47 (1999).

58 B. Dittrich and M. A. Spackman, “Can the interaction density be measured? The example of the non-standard amino acid sarcosine,” *Acta Cryst.* **A** **63**, 426–436 (2007).
<footnote n="59">⁵⁹ M. Ernst, A. Genoni, and P. Macchi, “Analysis of crystal field effects and interactions using X-ray restrained ELMOs,” *J. Mol. Struct.* **1209**, 127975:1–11 (2020).</footnote>

<footnote n="60">⁶⁰ G. Macetti, P. Macchi, and A. Genoni, “X-ray restrained extremely localized molecular orbitals for the embedding of quantum mechanical calculations,” *Acta Cryst.* **B 77**, 695–705 (2021).</footnote>

<footnote n="61">⁶¹ L. Bučinský, D. Jayatilaka, and S. Grabowsky, “Importance of relativistic effects and electron correlation in structure factors and electron density of diphenyl mercury and triphenyl bismuth,” *J. Phys. Chem. A* **120**, 6650–6669 (2016).</footnote>

<footnote n="62">⁶² L. Bučinský, D. Jayatilaka, and S. Grabowsky, “Relativistic quantum crystallography of diphenyl- and dicyanomercury. Theoretical structure factors and Hirshfeld atom refinement,” *Acta Cryst.* **A 75**, 705–717 (2019).</footnote>

<footnote n="63">⁶³ M. Podhorský, L. Bučinský, D. Jayatilaka, and S. Grabowsky, “HgH₂ meets relativistic quantum crystallography. How to teach relativity to a non-relativistic wavefunction,” *Acta Cryst.* **A 77**, 54–66 (2021).</footnote>

<footnote n="64">⁶⁴ S. Pawlędzio, M. Malinska, M. Woińska, J. Wojciechowski, L. A. Malaspina, F. Kleemiss, S. Grabowsky, and K. Woźniak, “Relativistic Hirshfeld atom refinement of an organo-gold (I) compound,” *IUCR J* **8**, 608–620 (2021).</footnote>

<footnote n="65">⁶⁵ L. Krause, B. Niepötter, C. J. Schürmann, D. Stalke, and R. Herbst-Irmer, “Validation of experimental charge-density refinement strategies: when do we overfit?” *IUCR J* **4**, 420–430 (2017).</footnote>

<footnote n="66">⁶⁶ V. V. Zhurov, E. A. Zhurova, and A. A. Pinkerton, “Optimization and evaluation of data quality for charge density studies,” *J. Appl. Cryst.* **41**, 340–349 (2008).</footnote>

<footnote n="67">⁶⁷ L. A. Malaspina, A. Genoni, D. Jayatilaka, M. J. Turner, K. Sugimoto, E. Nishibori, and S. Grabowsky, “The advanced treatment of hydrogen bonding in quantum crystallography,” *J. Appl. Cryst.* **54**, 718–729 (2021).</footnote>

<footnote n="68">⁶⁸ B. Landeros-Rivera, J. Contreras-García, and P. M. Dominiak, “The influence of refinement strategies on the wavefunctions derived from an experiment,” *Acta Cryst.* **B 77**, 715–727 (2021).</footnote>

<footnote n="69">⁶⁹ H. Birkedal, D. Madsen, R. H. Mathiesen, K. Knudsen, H.-P. Weber, P. Pattison, and D. Schwarzenbach, “The charge density of urea from synchrotron diffraction data,” *Acta Cryst.* **A 60**, 371–381 (2004).</footnote>

<footnote n="70">⁷⁰ D. Jayatilaka and B. Dittrich, “X-ray structure refinement using aspherical atomic density functions obtained from quantum-mechanical calculations,” *Acta Cryst.* **A 64**, 383–393 (2008).</footnote>

<footnote n="71">⁷¹ M. E. Wall, “Quantum crystallographic charge density of urea,” *IUCR J* **3**, 237–246 (2016).</footnote>

<footnote n="72">⁷² A. A. Hoser and A. Ø. Madsen, “Dynamic quantum crystallography: lattice-dynamical models refined against diffraction data. I. Theory,” *Acta Cryst.* **A 72**, 206–214 (2016).</footnote>

<footnote n="73">⁷³ R. Destro, R. E. Marsh, and R. Bianchi, “A low-temperature (23 K) study of L-alanine,” *J. Phys. Chem.* **92**, 966–973 (1988).</footnote>
<footnote n="74">74 C. Gatti, R. Bianchi, R. Destro, and F. Merati, "Experimental vs. theoretical topological properties of charge density distributions. An application to the L-alanine molecule studied by X-ray diffraction at 23 K," *J. Mol. Struct. Theochem* **255**, 409-433 (1992).</footnote>

<footnote n="75">75 P. Munshi, A. Ø. Madsen, M. A. Spackman, S. Larsen, and R. Destro, "Estimated H-atom anisotropic displacement parameters: a comparison between different methods and with neutron diffraction results," *Acta Cryst. A* **64**, 465-475 (2008).</footnote>

<footnote n="76">76 B. Dittrich and D. Jayatilaka, "Reliable measurements of dipole moments from single-crystal diffraction data and assessment of an in-crystal enhancement," *Structure and Bonding, Electron Density and Chemical Bonding II* **147**, pp. 27-46 (2012).</footnote>

<footnote n="77">77 L. A. Malaspina, E. K. Wieduwilt, J. Bergmann, F. Kleemiss, B. Meyer, M. F. Ruiz-López, R. Pal, E. Hupf, J. Beckmann, R. O. Piltz, A. J. Edwards, S. Grabowsky, and A. Genoni, "Fast and accurate quantum crystallography: from small to large, from light to heavy," *J. Phys. Chem. Lett.* **10**, 6973-6982 (2019).</footnote>

<footnote n="78">78 E. K. Wieduwilt, G. Macetti, L. A. Malaspina, D. Jayatilaka, S. Grabowsky, and A. Genoni, "Post-Hartree-Fock methods for Hirshfeld atom refinement: are they necessary? Investigation of a strongly hydrogen-bonded molecular crystal," *J. Mol. Struct.* **1209**, 127934:1-13 (2020).</footnote>

<footnote n="79">79 M. A. Spackman and E. N. Maslen, "Electron density and the chemical bond. A reappraisal of Berlin's theorem," *Acta Cryst. A* **41**, 347-353 (1985).</footnote>

<footnote n="80">80 M. A. Spackman and E. N. Maslen, "Chemical properties from the promolecule," *J. Phys. Chem.* **90**, 2020-2027 (1986).</footnote>

<footnote n="81">81 R. T. Downs, G. V. Gibbs, M. B. Boisen, and K. M. Rosso, "A comparison of procrystal and ab-initio model representations of the electron-density distributions of minerals," *Phys. Chem. Minerals* **29**, 369-385 (2002).</footnote>

<footnote n="82">82 M. A. Spackman and P. G. Byrom, "A novel definition of a molecule in a crystal," *Chem. Phys. Lett.* **267**, 215-220 (1997).</footnote>

<footnote n="83">83 P. D. Walker and P. G. Mezey, "Molecular electron density lego approach to molecule building," *J. Am. Chem. Soc.* **115**, 12423-12430 (1993).</footnote>

<footnote n="84">84 L. Massa, T. Keith, Y. Cheng, and C. F. Matta, "The kernel energy method applied to quantum theory of atoms in molecules-energies of interacting quantum atoms," *Chem. Phys. Lett.* **734**, 136650:1-4 (2019).</footnote>

<footnote n="85">85 S. C. Capelli, H.-B. Bürgi, B. Dittrich, S. Grabowsky, and D. Jayatilaka, Hirshfeld atom refinement. *IUCrJ* **2014**, *1*, 361-379.</footnote>
<figure><img src="image_21.png" /></figure>

<footnote n="86">86 E. K. Wieduwilt, G. Macetti, and A. Genoni, “Climbing Jacob’s Ladder of structural refinement: Introduction of a localized molecular orbital-based embedding for accurate X-ray determinations of hydrogen atom positions,” *J. Phys. Chem. Lett.* **12**, 463–471 (2020).</footnote>

<footnote n="87">87 A. D. Becke, “Density-functional exchange-energy approximation with correct asymptotic behavior,” *Phys. Rev. A* **38**, 3098–3100 (1988).</footnote>

<footnote n="88">88 C. T. Lee, W. T. Yang, and R. G. Parr, “Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density,” *Phys. Rev. B* **37**, 785–789 (1988).</footnote>

<footnote n="89">89 J. Graciani, A. M. Márquez, J. J. Plata, Y. Ortega, N. C. Hernández, A. Meyer, C. M. Zicovich-Wilson, and J. F. Sanz, “Comparative study on the performance of hybrid DFT functionals in highly correlated oxides: the case of CeO₂ and Ce₂O₃,” *J. Chem. Theory Comput.* **7**, 56–65 (2011).</footnote>

<footnote n="90">90 J. P. Malrieu, R. Caballol, C. J. Calzado, D. de Graaf, and N. Guihery, “Magnetic interactions in molecules and highly correlated materials: physical content, analytical derivation, and rigorous extraction of magnetic Hamiltonians,” *Chem. Rev.* **114**, 429–492 (2014).</footnote>

<footnote n="91">91 D. Jayatilaka and D. J. Grimwood, “Tonto: afortan based object-oriented system for quantum chemistry and crystallography,” In *International Conference on Computational Science* (pp. 142–151). Springer, Berlin, Heidelberg. Revision 3.6.0 v. 3320ef5 (2003).</footnote>

<footnote n="92">92 M. F. Peintinger, D. V. Oliveira, and T. Bredow, “Consistent Gaussian basis sets of triple-zeta valence with polarization quality for solid-state calculations,” *J. Comp. Chem.* **34**, 451–459 (2013).</footnote>

<footnote n="93">93 M. Woińska, S. Grabowsky, P. M. Dominiak, K. Woźniak, and D. Jayatilaka, “Hydrogen atoms can be located accurately and precisely by x-ray crystallography,” *Sci. Adv.* **2**, e1600192 (2016).</footnote>

<footnote n="94">94 M. Fugel, D. Jayatilaka, E. Hupf, J. Overgaard, V. R. Hathwar, P. Macchi, M. J. Turner, J. A. K. Howard, O. V. Dolomanov, H. Puschmann, B. B. Iversen, H.-B. Bürgi, and S. Grabowsky, “Probing the accuracy and precision of Hirshfeld atom refinement with HARt interfaced with Olex2,” *IUCRJ* **5**, 32–44 (2018).</footnote>

<footnote n="95">95 Gaussian 09, Revision A.02, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi,</footnote>
J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, Gaussian, Inc., Wallingford CT, 2016.

⁹⁶ A. Volkov, T. Koritsanszky, M. Chodkiewicz, and H. F. King, “On the basis-set dependence of local and integrated electron density properties: Application of a new computer program for quantum-chemical density analysis,” *J. Comput. Chem.* **30**, 1379–1391 (2009).

⁹⁷ R. Pal, D. Jayatilaka, and E. Nishibori, “Structure factors and charge density description of aluminum: A quantum crystallographic study,” *J. Phys. Chem. A* **126**, 2042–2049 (2022).

⁹⁸ T. A. Jones, J.-Y. Zou, S. T. Cowan, and M. Kjeldgaard, Improved methods for building protein models in electron density maps and the location of errors in these models. *Acta Cryst.* **A** *47*, 110–119 (1991).

⁹⁹ W. Humphrey, A. Dalke, and K. Schulten, “VMD: visual molecular dynamics,” *J. Mol. Graphics* **14**, 33–38 (1996).

¹⁰⁰ R. Dovesi, R. Orlando, A. Erba, C. M. Zicovich-Wilson, B. Civalleri, S. Casassa, L. Maschio, M. Ferrabone, M. de la Pierre, P. d’Arco, Y. Noel, M. Causa, M. Rerat, and B. Kirtman, “CRYSTAL14: A program for the ab initio investigation of crystalline solids,” *Int. J. Quantum Chem.* **114**, 1287–1317 (2014).

¹⁰¹ W. Tang, E. Sanville, and G. Henkelman, “A grid-based Bader analysis algorithm without lattice bias,” *J. Phys. Cond. Matter* **21**, 084204:1–7 (2009).

¹⁰² A. Genoni, and B. Meyer, “X-Ray constrained wave functions: Fundamentals and effects of the molecular orbitals localization,” *Adv. Quant. Chem.* **73**, 333–362 (2016).

¹⁰³ M.-C. Kim, H. Park, S. Son, E. Sim, and K. Burke, “Improved DFT potential energy surfaces via improved densities,” *J. Phys. Chem. Lett.* **6**, 3802–3807 (2015).

¹⁰⁴ S. Song, S. Vuckovic, E. Sim, and K. Burke, “Density-corrected DFT explained: Questions and answers,” *J. Chem. Theory Comput.* **18**, 817–827 (2022).

¹⁰⁵ B. Dittrich, E. Sze, J. J. Holstein, C. B. Hübschle, and D. Jayatilaka, “Crystal-field effects in l-homoserine: Multipoles versus quantum chemistry,” *Acta Cryst.* **A** *68*, 435–442 (2012).
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>
<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>
<figure><img src="image_32.png" /></figure>
<figure><img src="image_33.png" /></figure>
<figure><img src="image_34.png" /></figure>
<figure><img src="image_35.png" /></figure>
<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>
<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>
<figure><img src="image_43.png" /></figure>
<figure><img src="image_44.png" /></figure>
<figure><img src="image_45.png" /></figure>
<figure><img src="image_46.png" /></figure>
<figure><img src="image_47.png" /></figure>
<figure><img src="image_48.png" /></figure>
<figure><img src="image_49.png" /></figure>
<figure><img src="image_50.png" /></figure>
<figure><img src="image_51.png" /></figure>
<figure><img src="image_52.png" /></figure>
<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /></figure>
<figure><img src="image_55.png" /></figure>
<figure><img src="image_56.png" /></figure>
<figure><img src="image_57.png" /></figure>
<figure><img src="image_58.png" /></figure>

Accepted
10.1039/c8bp38312
<figure><img src="image_59.png" /></figure>
<figure><img src="image_60.png" /></figure>
<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /></figure>

<figure><img src="image_63.png" /></figure>
<figure><img src="image_64.png" /></figure>
<figure><img src="image_65.png" /></figure>

<figure><img src="image_66.png" /></figure>

<figure><img src="image_67.png" /></figure>
<figure><img src="image_68.png" /></figure>
<figure><img src="image_69.png" /><figcaption>red</figcaption></figure>
<figure><img src="image_70.png" /></figure>
<figure><img src="image_71.png" /></figure>
<figure><img src="image_72.png" /></figure>
<figure><img src="image_73.png" /></figure>
<figure><img src="image_74.png" /></figure>
<figure><img src="image_75.png" /></figure>
<figure><img src="image_76.png" /></figure>
<figure><img src="image_77.png" /></figure>
<figure><img src="image_78.png" /></figure>
<figure><img src="image_79.png" /></figure>
<figure><img src="image_80.png" /></figure>

<figure><img src="image_81.png" /></figure>
<figure><img src="image_82.png" /></figure>
<figure><img src="image_83.png" /><figcaption>4</figcaption></figure>
<figure><img src="image_84.png" /></figure>
<figure><img src="image_85.png" /></figure>
<figure><img src="image_86.png" /></figure>
<figure><img src="image_87.png" /></figure>
<figure><img src="image_88.png" /></figure>