Quantum Crystallography in the Last Decade: Developments and Outlooks

: In this review article, we report on the recent progresses in the field of quantum crystallography that has witnessed a massive increase of production coupled with a broadening of the scope in the last decade. It is shown that the early thoughts about extracting quantum mechanical information from crystallographic experiments are becoming reality, although a century after prediction. While in the past the focus was mainly on electron density and related quantities, the attention is now shifting toward determination of wavefunction from experiments, which enables an exhaustive determination of the quantum mechanical functions and properties of a system. Nonetheless, methods based on electron density modelling have evolved and are nowadays able to reconstruct tiny polarizations of core electrons, coupling charge and spin models, or determining the quantum behaviour at extreme conditions. Far from being routine, these experimental and computational results should be regarded with special attention by scientists for the wealth of information on a system that they actually contain.


Introduction
The interplay between crystallography and quantum mechanics is very tight and long-standing. It started in the early 20 th century, when both the modern (X-ray based) crystallography and Bohr's first quantum physics began. In fact, the first X-ray diffraction experiments did not escape the attention of scientists who understood the potential of this technique and its implications for the on-going scientific revolution. Prophetic was a statement by Peter Debye stated in 1915: "it seems to me that the experimental study of the scattered radiation, in particular from light atoms, should get more attention, since along this way it should be possible to determine the arrangement of the electrons in the atoms" [1]. The suggestion was intended to foster the development of a suitable model for electron movement in atoms by means of experimental evidence, replacing the physically inconsistent Bohr model, at that time in vogue [2]. However, the second (correct) quantum physics, better known as quantum mechanics, took only partial advantage of the experimental observations made through X-rays. Instead, it is more appropriate to say that crystallography took advantage of the newly proposed atomic quantum model. In fact, a full quantum treatment was necessary to replace the original atomic scattering factors based on the classical Thomson theory adapted to the Bohr atomic model [3]. Using those form factors, the predicted diffracted intensities differed quite significantly from the experimental observations by Bragg, James and Bosanquet [4]. Only when the more appropriate atomic quantum electron densities (proposed by Waller and Hartree [5]) were adopted, a more satisfactory validation of crystal structure models was possible. The subsequent development of crystallography and the implications for theoretical chemistry, especially those concerning the development of chemical bonding theories, are well known. This strengthened the correlation between quantum chemistry and crystallography, although without using experimental crystallography (and particularly, diffraction techniques) for a direct determination of quantum mechanical properties or functions.
The purpose of extracting quantum mechanical information from crystallographic experiments was resumed much later by Richard Weiss, who tried to determine the electronic configurations of metals from X-ray diffraction on single crystal samples [6] and later pioneered the investigations of electron density in momentum space [7], as well as those of the spin density distribution [8]. The main goal of Weiss' research was of retrieving the wavefunction directly from X-ray experiments. In his book on X-ray determination of electron distribution [9], he suggested that experimentally measured scattered intensities could be employed to correct the Hartree-Fock molecular wavefunction, which neglects electron correlation. This implies a merge of the experimental and theoretical frameworks, as well as the combination of the direct and Fourier transformed position space. Since crystals are periodical objects and when they are illuminated by a coherent radiation a discrete Fourier transformation is obtained, Weiss' intuition basically originated the field of quantum crystallography, although the term itself was coined only many years later and for a narrower scope by Karle, Massa and Huang [10,11]. They defined quantum crystallography as the "combining of crystallographic data with quantum mechanical techniques in such a way that it should be possible to obtain information of enhanced value. The enhancement could be increased accuracy or information more readily obtained, or both".
The quantum crystallographic research at the end of 1960s and early 1970s split into two main streams, one based on atomic functions and the other one on molecular functions. We may view these approaches as representations of the atoms in crystals or of the molecules in crystals, respectively.
The basis of the former was set by Stewart who stressed a concept earlier introduced by Dawson, namely the generalized atomic scattering factor [12]. Instead of using atomic form factors computed quantum mechanically for isolated noninteracting atoms, Stewart's idea was that of projecting onto atomic bases the calculated molecular electron densities [13] and from them derive the atomic form factors to be employed for the refinement procedures [14]. The calculation of molecular wavefunctions was anyway too complicated and the model was eventually refurbished with a set of atomic parameters to be directly refined against experimental observations [15]. The parameters are nonetheless reminiscent of theoretical atomic wave functions. In the refinement procedure, the electron density of the crystal is projected onto these basis functions and parameters are refined in order to minimize the difference with respect to experimental intensities. Thus, the so-called pseudoatom model proposed by Stewart [15] was practically an atom-centred multipolar expansion, extended up to a reasonable finite level. In research papers, the model was called also aspherical atom model or multipolar model, as it is mostly labelled today. The most adopted formalism is actually the one proposed by Hansen and Coppens [16], who devised a useful adaptation with a locally defined orientation system, having in mind the possible exportability of the atomic models. While the Stewart and Coppens models are the most well-known, in the same period, a number of other models were proposed, all of them based on very similar ideas, although lacking of a quantum mechanical basis unlike Stewart's pseudoatom model. Among them, one has survived quite long and actually has become extremely popular not only in crystallography: the Hirshfeld atom [17]. Here the atomic projection uses as basis the ground state (spherical) electron density of the isolated atom, in practice the same atomic electron density used for conventional crystal structure refinements. In this way, the electron density is partitioned on the basis of atomic weights that simply correspond to the atomic fraction of the total (molecular/crystal) electron density calculated as sum of atomic spherical electron densities. The Hirshfeld atom has been therefore used as a very-easy-to-implement partitioning scheme for the electron density and works regardless of the way in which the density itself is obtained (either through a multipolar model or a wavefunction). In quantum chemistry and crystallography, there are three main applications of the Hirshfeld atom: a) the atomic population analysis of experimental and theoretical electron densities, especially within the framework of density functional theory (DFT) [18]; b) the definition of molecular Hirshfeld surfaces [19]; c) the Hirshfeld atom refinement (HAR) [20].
The other stream of research was instead dedicated to attempts of directly obtaining molecular wavefunctions from experimental diffraction. Coppens and coworkers, for example, tested the possibility of refining directly the coefficients of a properly designed molecular wavefunction, but realized how complicated this would be even for medium-size molecules [21]. One of the major problems to solve was the need of refining coefficients for the two-centre products, i.e. products of orbital functions centred on two different atoms. These complications pushed Coppens as well toward projected density functions onto atomic basis functions.
On the contrary, Clinton and Massa [22][23][24][25][26][27][28][29][30] proposed the calculation of density matrices constrained to the representation of a given quantum mechanical operator. As the scattering operator is an operator belonging to this category, the method could be used to obtain what, in principle, can be considered as experimental density matrices. The fulfilment of the N-representability conditions for the (one-electron) reduced density matrices also enabled the automatic and consequent determination of "experimental" wavefunctions. As mentioned above, it is mainly within this stream that the term quantum crystallography was coined, although the name quantum crystal was actually introduced much earlier by Nosanow [31] to define nonclassical behaviour of atoms in crystalline matrices at very low temperatures. However, this definition has not found much application in the literature and it is no longer in use.
In this review article, we will present some of the major research breakthroughs in quantum crystallography occurred in the past decade. The research was very lively and produced quite a number of new or improved methods. For this reason, we will mainly focus on recent advancements connected to the two research paths described above. Namely, we will consider only those studies devoted to directly extract electron densities, density matrices or wave functions from experimental data rather than on other techniques that exploit electron densities resulting from standard quantum mechanical calculations to refine crystal structures. A prominent example of these latter techniques is HAR, which, however, in the last years emerged as a very promising quantum crystallographic Xray refinement method able to provide hydrogen bond lengths in excellent agreement with those resulting from neutron diffraction measurements [32][33][34][35][36].

Fitting the Wavefunction
As already mentioned above, the determination of wavefunctions and density matrices from experimental scattering data has been a tantalizing perspective since the early days of quantum crystallography. In fact, according to the postulates of quantum mechanics, wavefunctions and density matrices intrinsically contain all the information about a system and the possibility of obtaining those objects from experiments would automatically lead to evaluating all the properties of interest for the system under examination.
In this context, different attempts have been proposed over the years. Just to cite a few of them, we mention the pioneering works of Clinton, Massa and their coworkers since the early 1960s [22][23][24][25][26][27][28][29][30] or the more recent molecular orbital occupation number (MOON) approach [37,38] and Tanaka's X-ray atomic orbital (XAO) [39] and X-ray molecular orbital (XMO) methods [40]. Nevertheless, within this family of techniques, the strategies that emerged the most or witnessed a significant evolution in the last decade were essentially two: i) the X-ray constrained/restrained wavefunction (XCW/XRW) method originally devised by Jayatilaka and ii) the approaches proposed by Gillet and collaborators to refine (spin-resolved) one-particle reduced density matrices (1-RDMs). The recent applications and developments of these two strategies in the last years will be reviewed in the following paragraphs.
The technique introduced by Jayatilaka [41,42] stems from the observation that, according to the Gilbert corollary [43] to Coleman's theorem [44], an infinite number of density matrices (and, therefore, wavefunctions) may actually correspond to one given electron density, with the implicit consequence that, unlike experimental electron densities (see section 3 of this review), "experimental" wavefunctions cannot be obtained through a simple fitting of parameters (e.g., molecular orbitals coefficients) against experimental data, as also pointed out by Coppens in the 1970s [21]. Prompted by this fact and by a possible solution proposed by Henderson and Zimmermann [45] to circumvent the problem, Jayatilaka adopted a different strategy that allows the extraction of plausible wavefunctions that not only reproduce experimental X-ray diffraction data within the limit imposed by the unavoidable experimental uncertainties, but also minimize the energy of the system under investigation [41,42]. In other words, the variational principle, at the heart of many methods of quantum chemistry, remains a guide to obtain the wavefunction, but with the addition of some restraints having the form of the statistical agreement between experimental and theoretical structure factors amplitudes. In a nutshell, the method consists in determining the wavefunction parameters (usually the coefficients of the molecular orbitals) that minimize the following functional: where the first term represents the energy of the system, while the second one accounts for the fitting of the experimental X-ray diffraction data; is an external multiplier that is adjusted during the calculations and determines the weight of the fitting in the procedure, Δ is the desired average agreement between theoretical and observed structure factor amplitudes, which is given by the reduced 2 statistical agreement given by the following expression: with as the number of experimental data, as the number of adjustable parameters, as the triad of Miller indices for the reflection under exam, as the uncertainty corresponding to the experimental structure factor amplitude � � and as an overall multiplicative factor, which puts the computed structure factor amplitudes � � on the same scale of � �. The method is mainly known in the literature as the X-ray constrained wavefunction (XCW) fitting approach. However, in the absence of a defined target for the agreement with the experimental data (note that Δ in equation (1) is a desired average agreement), it is more appropriate to see the "perturbation" of the experimental data as a restraint rather than a real constraint. Therefore, in the rest of this manuscript we will refer to Jayatilaka's strategy as X-ray restrained wavefunction (XRW) fitting technique.
It is worth mentioning that, in the ideal case, the XRW calculations should stop when the reduced 2 is equal to 1.0 (i.e., when the calculated structure factor amplitudes are on average within one standard deviation of the experimental data). However, also depending on the quality of the experimental measurements used as restraints, reaching the ideal agreement with the experimental values is not always possible. Therefore, alternative convergence criteria to halt the search of X-ray restrained wavefunctions have been proposed by different research groups working in the field, although this is still an open question that will have to be addressed in the near future.
Although initially proposed in 1998 [41] and consolidated in the early 2000s [42,[46][47][48][49][50], the XRW technique has seen a vivid revival and improvement in the last ten years, both in terms of applications and methodological development.
The first applications to materials science were proposed by Spackman, Jayatilaka and coworkers, who exploited XRWs to extract effective molecular (hyper)polarizabilities and crystal refractive indices from X-ray diffraction measurements on some simple molecular crystals [51,52]. Their studies were motivated by the fact that, unlike what was generally and erroneously believed, those kinds of properties cannot be correctly determined only on the basis of electron density distributions (e.g., from electron density distributions obtained through traditional multipole model refinements) because fundamental many-electron contributions would be missing. A wavefunction would not suffice either, because the polarizability requires a sum over all states, especially those closer in energy to the ground state. However, thanks to a proper approximation, the ground state wavefunction may be sufficient to estimate the molecular polarizability. Spackman and Jayatilaka demonstrated that also an X-ray restrained wavefunction may serve this purpose because it correctly accounts for effects of the crystal environment and can be treated as if it was a true ground state wavefunction. The refractive indices calculated for a few molecular crystals were in very good agreement with those from direct measurements and extrapolated to the limit of an infinite wavelength, as the calculations do not include the electric field frequencies. Later on, Cole et al. [53,54] further exploited Jayatilaka's method for a series of compounds exhibiting interesting nonlinear optical properties. The studies confirmed the capability of the X-ray restrained wavefunction fitting method in capturing solid-state effects in the condensed phase and the superiority of the XRW approach compared to the multipole model techniques when properties depending on two-electron terms are concerned.
The capability of the Jayatilaka approach in accounting for crystal field effects was also the topic of a very recent and focused work [55], in conjunction with a complementary investigation aimed at evaluating the amount of electron correlation effects [56]. The two effects are inherently contained in converged X-ray restrained wavefunctions because the procedure can be seen as a correction of the Hamiltonian for both effects. Both studies showed that Jayatilaka's method can definitely retrieve the deformations of electron densities and molecular orbitals due to intra-crystal electric fields and electron correlation, although it was also pointed out that the resolution of the diffraction data used as external restraints plays a crucial role. In fact, the recovery of the two effects surprisingly decreases with the increase of the resolution. This is ascribed to the down-weighting of the low-angle reflections, which are extremely important for a good modelling of the valence electrons (namely, the electrons most affected by the above-mentioned effects), but which become a lower fraction of the collected data when one increases the diffraction resolution.
In this respect, the wavefunctions obtained with the XRW method were adopted for accurate chemical bonding analyses. The goal was to unequivocally quantify the influence of the environment on the system under investigation and to recover traditional chemical pictures from "experimental" wave functions in order to rationalize interesting and controversial chemical problems. The first example in this direction is the investigation conducted by Jayatilaka and Grimwood, who computed electron localization functions (ELFs) for X-ray restrained Hartree-Fock wavefunctions corresponding to different molecular crystals and clearly pointed out the differences compared to the corresponding gas phase results [57].
Following this direction, Grabowsky and coworkers afterwards applied other chemical bonding analyses in conjunction with the Jayatilaka approach. They determined XRW-based electron localizability indicator domains to shed light on substituent and crystal effects for a series of acceptorsubstituted epoxide derivatives [58] and to rationalize the reactivity differences between α,βunsaturated carbonyl and hydrazine compounds [59]. Grabowsky et al. also combined the XRW approach with the delocalization index δ, the Roby bond index τ, the electron localizability indicator ELI and the quantum theory of atoms in molecules (QTAIM) approach to tackle the controversial problem of hypervalency in sulphur dioxide (SO2). The consensus bond orders of about 1.5 resulting from their analyses showed that the S-O bonds are mainly characterized by a ionic and multicentre character, thus excluding the possibility of hypervalent resonance structures and explaining the shortening of the S-O bonds in terms of electrostatic forces associated with the intrinsic ionicity of the compound [60]. In a follow-up study, the concept of hypervalency was also re-examined for the phosphate, sulphate and perchlorate anions, always on the basis of chemical bonding analyses carried out on X-ray restrained wavefunctions. While hypervalency of phosphorus and sulphur atoms was excluded for the description of P-O and S-O bonds, respectively, the analysis suggested a hypervalent character for the chlorine atom, probably due to the hyperconjugation of the p-type oxygen lone pairs with the σ* molecular orbitals associable with the Cl-O bonds [61]. Other notable examples of a combined use of the Jayatilaka approach with traditional methods of the chemical bonding analysis are the investigations conducted by Thomas et al., as the one that aimed at characterizing the nature of the intramolecular S···O chalcogen bond in acetazolamide [62].
So far, we have mainly focused on the most recent applications of the XRW approach to chemical and physical problems. However, as already anticipated, in the last ten years, the technique also went through methodological developments. In fact, since the Jayatilaka strategy was originally proposed only in the framework of the basic restricted Hartree-Fock and restricted DFT formalisms, several efforts have been recently done to couple the X-ray restrained wavefunction philosophy with other methods of quantum chemistry.
Bučinský and collaborators developed the unrestricted and the relativistic versions of the X-ray restrained wavefunction approach [63][64][65][66]. This clearly opened new possibilities for the electronic structure investigation of solid-state systems. On the one hand, the use of the unrestricted formalism gave access to experimental spin densities (see Figure 1) obtained at X-ray restrained level, comparable to those resulting from multipole-model or 1-RDM joint refinements (see section 3 and the last part of this section for more details). On the other hand, the new relativistic XRW methods paved the way to use Jayatilaka's strategy for chemical/physical properties of crystal compounds containing heavy elements and, therefore, strongly affected by important relativistic effects (see again Figure 1). Another stream of methodological development aimed at directly introducing the traditional chemical perception into the X-ray constrained wavefunction formalism without the need of resorting to chemical bonding analyses. Most of these developments were based on the concept of extremely localized molecular orbitals (ELMOs) [67][68][69][70][71][72][73], namely molecular orbitals strictly localized on small molecular fragments, which can be obtained by defining a priori (and according to the chemical intuition) a localization scheme that subdivides the molecule under exam into subunits. Therefore, by introducing this fragmentation into the machinery of Jayatilaka's approach, it has been possible to extract X-ray restrained molecular orbitals strictly corresponding to atoms, bonds, lonepairs and functional groups [74][75][76][77]. In other words, it was possible to obtain X-ray restrained molecular orbitals very close to traditional concepts for chemists. This had the non-negligible consequence of getting back for the XRW technique one of the main features of the multipole model strategy, namely the possibility of looking at the global electron density as the sum of contributions coming from elementary units. We could even state that X-ray restrained ELMOs are for the Jayatilaka approach what pseudoatoms are for the multipole models. To further pursue the parallelism, it is worth noting that XR-ELMOs are as transferable entities as multipole model pseudoatoms. For this reason, in the future, they could be possibly exploited to extend the recently constructed libraries of theoretical extremely localized molecular orbitals [78][79][80], which have been already used to develop multiscale quantum chemistry embedding techniques [81,82], quickly detect noncovalent interactions in large systems [83], and refine crystal structures of polypeptides and small proteins in the framework of the Hirshfeld atom refinement [84].
Extremely localized molecular orbitals were also used in combination with the Jayatilaka philosophy to develop the so-called X-ray restrained ELMO valence bond (XR-ELMO-VB) method [85,86], which can be considered the first prototype multideterminant X-ray restrained wavefunction technique. The strategy was introduced to determine the weights of molecular resonance structures from X-ray diffraction data. In this method, the global wavefunction is written as a linear combination of predetermined unrestrained ELMO wavefunctions associated with the different resonance structures of the investigated system: By keeping fixed the precomputed ELMO wavefunctions �Ψ , �, the strategy consists in determining the coefficients � � that simultaneously minimize the electronic energy and the statistical disagreement between experimental and computed structure factors magnitude (namely, minimization of the usual Jayatilaka functional given by equation (1) with respect to the coefficients � �). It is important to note that these coefficients do not directly provide the weights of the resonance structures because the ELMO wavefunctions �Ψ , � in expansion (3) are nonorthogonal. In order to get the weights, one needs to resort to the corresponding Chirgwin-Coulson coefficients: With �Ψ , �Ψ , � as the overlap integral between the ELMO wave functions Ψ , and Ψ , . The technique was interestingly used to further support the conclusions of a recent charge density investigation conducted on the syn-1,6:8-13-biscabonyl [14]annulene (BCA), which indicated a partial suppression of the compound aromaticity when pressure increases [87]. The XR-ELMO-VB computations confirmed the trend showing that, while at ambient pressure the two resonance structures of BCA are basically equivalent, at high pressure one of the two becomes clearly predominant [86] (see Figure 2). Always with the goal of extracting useful chemical information from experiment X-ray diffraction data, the more recent X-ray restrained spin-coupled (XRSC) approach has been also introduced [88,89]. It results from the coupling of the Jayatilaka method with the spin-coupled technique of the valence bond theory [90][91][92] and it can be considered as a step forward compared to the previous XR-ELMO-VB strategy. In fact, in the novel XRSC method, it is possible to extract both orbitals and resonance structure weights from X-ray data without providing a priori any preliminary information, namely without specifying any localization scheme or using any pre-computed wavefunction in the calculations. Test calculations have shown that spin-coupled orbitals and resonance structure weights obtained from XRSC calculations present non-negligible differences compared to those obtained through traditional gas-phase spin-coupled computations. This further indicates the intrinsic capability of the Jayatilaka approach of accounting for the effects of the crystal field on the molecular electronic structure.
Finally, it is also worth noting that recently the X-ray restrained wavefunction (XWR) fitting method has been coupled with the Hirshfeld atom refinement (HAR), giving rise to the so-called Xray wavefunction refinement (XRW) technique [60,93]. In this strategy, HAR and XWR fitting are alternated until convergence. The former obviously refines structural parameters (i.e., atomic positions and thermal parameters), while the latter accounts for the electronic ones (e.g., molecular orbital coefficients in case of a single Slater determinant wavefunction). Test refinements conducted on amino acids and polypeptides have indicated that the new XRW approach seems to provide crystal structures and electron densities in much better agreement with the experimental diffraction data than those resulting from traditional multipole model refinements [93].
As already anticipated at the beginning of this section, other than the Jayatilaka approach, another group of techniques that recently showed promising advances in the context of extracting wavefunctions/density matrices from experimental data is the one represented by the methods proposed by Gillet and coworkers, who particularly aimed at reconstructing (spin-resolved) 1-RDMs by simultaneously considering experimental data obtained by means of different kinds of scattering experiments.
This research field stems from the observation that the diagonal part of the one-particle reduced density matrix (namely, the electron density) can be obtained only from elastic (Bragg) X-ray diffraction, whereas the off-diagonal parts are associated with inelastic (Compton) scattering. At the same time, if we focus on the one-particle spin density matrix, it is worth noting that its diagonal part (namely, the usual spin-density) can be reconstructed from polarized neutron diffraction (PND) measurements or, in principle, from magnetic X-ray diffraction, whereas the off-diagonal terms can be potentially obtained only from magnetic Compton scattering experiments. This is the context of Gillet and collaborators' most recent studies, who also took inspiration from the first pioneering investigations in this field, such as those of Schmider, Smith and Weyrich [94][95][96], and those of Becker, Gillet himself and Cortona [97,98].
In particular, two main research directions were explored: i) the reconstruction of 1-RDMs by simultaneous refinement of X-ray diffraction and directional Compton scattering data; ii) the determination of spin-resolved 1-RDMs from polarized neutron diffraction and magnetic Compton scattering measurements.
Pertaining to the first aspect, in 2007 Gillet devised a method that extended the well-known Hansen and Coppens multipole model for electron density distributions to one-particle density matrices [99]. Although it was tested only on two diatomic systems (HF and CO), the strategy showed that data from complementary experiments are useful to retrieve important chemical bonding features. More recently, a completely different strategy has been adopted. In fact, instead of following the multipole-model formalism, De Bruyne and Gillet preferred to express the one-particle density matrix in terms of orthogonalized atom-centred basis functions [100]. Consequently, they determined the corresponding population-matrix elements by means of a constrained least-squares fitting scheme that took into account the needed N-representability conditions for 1-RDMs through the use of semidefinite programming. Preliminary test calculations on the crystal of dry ice provided a very good agreement with the results of periodic ab initio computations (see Figure 3). Concerning the second research direction, after preliminary investigations where the spin density of YTiO3 was reconstructed both in position and momentum space from polarized neutron diffraction data and theoretical/experimental incoherent x-ray magnetic Compton scattering profiles [101,102], a more advanced model to reconstruct spin-resolved one-electron reduced density matrices has been proposed [103]. As in the recent paper by De Bruyne and Gillet [100], the spin-resolved 1-RDM is expanded in terms of atom-centred Gaussian basis functions, with the only difference that, for this method, the exponents of the Gaussians are not fixed, but can be rescaled during the refinement procedure. Therefore, once the optimal scale-factor for the Gaussian exponents is defined, the elements of the spin population-matrix are determined by minimizing an objective function that accounts for the agreement with the available experimental data (in this case, magnetic structure factors and magnetic Compton profiles). Again, the fulfilment of N-representability conditions is imposed in order to obtain a quantum mechanically rigorous spin-resolved one-particle reduced density matrix. The preliminary test of the model was conducted on an artificial magnetic crystal of urea and showed that the new joint-refinement method enables to get more accurate results compared to those strategies that take into account only polarized neutron diffraction data [103]. In fact, it was clearly shown that accounting for magnetic Compton scattering profiles does not affect only the off-diagonal terms of the spin-resolved 1-RDM, but also the diagonal ones, thus leading to a much better global description of the spin density (see Figure 4). The technique was afterwards used to determine the spin-resolved one-electron density matrix of YTiO3, confirming the results obtained from the preliminary test refinements and enabling the study of the magnetic properties of the crystal along the Ti-O-Ti bonding pattern [104].

Fitting the Density
In principle, the scattered wave is an observable and, as such, is associated with a quantum mechanical operator. This implies that, according to the postulates of quantum mechanics, if we have the wavefunction describing the state of the system under exam, the average value of the scattered wave can be obtained as expectation value of the above-mentioned operator, namely: The scattering operator for a system of electrons is the one-electron operator ̂= ∑

• =1
, where is a vector that is the difference between the incident and diffracted wave vectors, and is the position of electron j. In practice, this operator sums the waves emitted by each electron under the perturbation of the radiation electric field along a given direction and the wavefunction defines the probability of finding each electron at each position in space.
If we integrate over all electrons but one and we take into account that electrons are indistinguishable, equation (5) becomes where ( ) is the one-electron density and ⟨ ⟩ is the measurable scattered wave, which obviously depends on the scattering vector . As is well known, we cannot measure the whole scattered wave, but only a quantity that is proportional to the square of its amplitude.
Because of equation (6), it is quite usual to think that the scattered radiation is causally related to the one-electron density more than to the wavefunction, despite Equation (6) is at the heart of all electron density maps used by crystallographers to solve crystal structures or inspect refined models through the so-called deformation density maps, where the electron density of a model (for example the sum of spherical atoms) is taken as a benchmark. Equation (6) easily explains the reason why the refinement of an electron density function from X-ray diffraction data is much easier than that of a wavefunction. This roadmap was pursued by Stewart, Coppens and others in the early 1970s, after they realized that a direct refinement of wavefunction coefficients was too complicated and the alternative methods to obtain an experimental wavefunction (see section 2) were not consolidated yet.
This led to the so-called multipolar model, i.e., a technique where the electron density is refined as a parameterized model, under the assumption that it can be partitioned into atomic contributions. A schematic formula is: Where , ℎ− ( ) and , ℎ− ( ) are the spherical core and spherical valence contributions to the atomic electron density, whereas , ℎ ( ) is the aspherical or deformation part, which is fitted by a series of atom-centred spherical harmonics extended up to a given level (generally, up to = 4). The aspherical part is typically considered only at the level of the valence shell, as the core is often taken as unperturbed with respect to the atomic ground state (meaning that the electron density is spherically distributed and it strictly corresponds to the number of core electrons).
Each term in equation (8) depends on the positions of the atoms, as well as on many coefficients of the density functions, as, for example, the number of electrons associated with the core or valence density or the fraction of electrons shifted around the atom in a dipolar, quadrupolar, octupolar, etc., shape. Moreover, some additional parameters describe how the atomic density expands or contracts with respect to the isolated ground state density.
The most adopted version of multipolar model is the one proposed by Hansen and Coppens [16] although other modifications have been used as well. Over the years, the model has been improved and adapted to different situations. Studies in the late 1990s and early 2000s were mainly dedicated to improving the quality of the radial functions adopted for , ℎ− ( ) and especially for , ℎ ( ) [105]. Connected with that was also the development of recipes for refining the contraction or expansion of the atomic densities [106], which are difficult to determine from diffraction intensities.
Instead, in the past decade, the improvement mainly concerned the possibility of modelling core deformations or modelling simultaneously the electron charge and spin density distributions. Last but not least, a number of databank schemes have been proposed, each of them intended to store transferable atomic electron density parameters for modelling large or complicated molecules (especially biological macromolecules) that are more problematic to investigate accurately, as required for a proper charge density analysis.
All these improvements are possible thanks to the inherent flexibility of a model like the Hansen and Coppens's one. A special technical feature is the possible definition of a local coordinate system for each atom in the crystal, instead of referring to the crystal orientation matrix. This enables the transferability of a set of multipoles from one atom to another of the same element, provided that it is also embedded in a sufficiently similar chemical environment. In fact, this was the purpose of the original model, although implemented only much later [107] when sufficient amount of data allowed the construction of databases of experimentally determined atomic multipoles for special classes of molecules [107,108], to be used to improve structural refinement of poorly diffracting species and to enable a sensible reconstruction and mapping of the electrostatic potential and electric moments. This approach was later generalized by means of theoretically calculated atomic multipoles, with the advantage of allowing the study of (macro) molecules for which experimental models cannot be easily devised [109][110][111]. While these approaches are now consolidated and their potential is wellknown, new applications emerge, such as those aiming at using these databases for modelling crystal structures measured with an emerging technique such as electron diffraction [112].
GThe flexibility of the multipole model was especially adopted to detect special features, such as the tiny deformations occurring to core electrons in atoms bonded in molecules or covalent solids [113]. For this reason, an extended Hansen and Coppens' model was proposed, where the aspherical atomic density refers not only to the valence, but also to the core electrons: Noteworthy, the contraction/expansion of the core-shell is also refined, which has a fundamental importance, not only for the chemical implications, but also for the correct determination of the atomic displacement parameters that otherwise would be biased. A small but significant corecontraction was for example observed on carbon in diamond, which is interpreted as a kind of reaction of the core density against the localized accumulation of electron charge in the chemical bonds.
Another important improvement of the multipolar model was proposed for a simultaneous description of electron charge and spin densities, which was made possible thanks to the combination of a technique which is sensitive to the electron charge density (like Bragg X-ray diffraction) and one which is subject to the electron magnetization (like the elastic neutron scattering). The latter is complicated by the (typically larger) cross section of atomic nuclei for neutrons that make the diffraction pattern a superposition of nuclear and magnetic scattering. In addition, the magnetic scattering is not only due to the atomic spin moments, but also to the atomic orbital angular moments. For this reason, for a proper determination of the spin density, a rather sophisticated technique is necessary, which is the diffraction of polarized neutrons from (quite large) single crystals. Once the two kinds of data are obtained, a joint refinement is possible, if the model is made more flexible by a further separation into a charge density for spin α electrons and a charge density for spin β electrons: The larger number of parameters exacerbates the correlation among them, in view of the limited number of data obtainable from polarized neutron diffraction experiments and the inherent linear dependency of the parameters for the X-ray diffraction part. One should also consider that even if Xray diffraction data are nowadays obtainable up to very high resolutions and, therefore, in very large number, only few of them (namely the lower resolution ones) contain information on the most deformable part of the electron density (the valence shell). This physical limit cannot be overcome, and it imposes limits to the flexibility of the multipolar model. However, thanks to a sensible combination of constraints, a joint refinement of electron charge and spin densities in paramagnetic metal complexes or magnetically active inorganic solids was possible [114][115][116], which paves the way toward applications in spintronics.
The model flexibility has often been adopted also for applications to measurements of species in unconventional conditions, such as under radiation excitation, under electric field, under high pressure or high temperature, or on powder samples instead of single crystals [117]. In all these circumstances, the classical prescriptions for accurate electron density modelling are not fulfilled. For example, in high-temperature studies, the atomic motion is clearly enhanced, which severely affects the deconvolution of the electron density from the thermal motion. Nevertheless, sometimes it is the environmental condition that makes a material interesting, and a quantum crystallographic study cannot be conducted in the ideal setting.
Among these applications, the study of crystals at high pressure has been proposed [86,87,118]. At variance from high temperature, pressure itself is not a drawback for an accurate mapping of quantum mechanical functions in crystals. However, the equipment typically adopted to measure crystals in those conditions (the diamond anvil cell) has a negative impact on the diffraction quality. In particular, it reduces the portion of measurable reciprocal space and it produces an additional and strong background, thus affecting a proper measurement of the diffracted intensities. Another pitfall is that pressure cannot be transmitted hydrostatically to a sample above a given pressure, which implies a frequent damage of the crystal quality, hence of the diffraction. It is instead easy to demonstrate that pressure reduces the atomic motion mimicking the effect of the temperature. This is particularly cogent for molecular crystals, especially if lacking of strong supramolecular interactions. From data available so far [86,87,119], at ca. 5-10 GPa (achievable hydrostatically) the atomic displacement parameters are reduced as much as with a cooling from ambient temperature to "liquid nitrogen temperature" (ca. 100K), which is the standard in typical electron density determinations. A further contraction to atomic displacements at "helium temperature" (ca. 10K) is instead much less likely since it would require a pressure in the so-called megabar regime (> 100 GPa), where hydrostaticity is impossible. Despite all the above-mentioned problems and drawbacks, some multipole refinements of molecular crystals under high pressure have been successfully reported [87], although the flexibility of the model could not be fully exploited.

Conclusions and Outlook
In this short review article, we have summarized the main last-decade-achievements and improvements concerning the "electronic structure methods" of quantum crystallography, which has witnessed a lively debate on new approaches, new techniques and new applications. In particular, the growing success of the wavefunction-based methods has enormously broadened the potential of the performed studies and the interconnections with other branches of crystallography, physics, materials science and biochemistry.
There is no doubt that the improved accuracy of quantum crystallographic methods, when coupled with simplified approaches, is able to attract the attention of chemists for the inherently richer information obtainable compared to a routine crystal structure determination. In fact, the challenge for the next decade is facilitating the technical applicability of sophisticated methods and improving the interpretation of the information therein contained.
On the other hand, the sophistication of quantum crystallographic methods makes them excellent candidates for another challenge: exploiting emerging techniques which are not based on X-ray diffraction, but on electron diffraction and other different forms of microscopy. Since all these methods explore crystals at a quantum mechanical level, the methodologies unavoidably need quantum crystallographic methods, otherwise the enormous potential of these investigations might easily vanish.