Charge densities in actinide compounds: strategies for data reduction and model building

A three part paper on charge density analysis of actinide containing compounds is presented that covers experimental protocols, augmented multipole model building and a comparison of topological analysis of experimental and theoretical electron densities.


Introduction
The nature of f-element bonding is poorly understood and presents a challenge from both experimental and theoretical perspectives (Kü chle et al., 1994;Prodan et al., 2007;Kaltsoyannis, 2013;Kaltsoyannis & Kerridge, 2014), especially for the early to middle actinides where relativistic effects cannot be neglected (Pyykkö, 1988;Onoe et al., 1993). Participation of the 5f or 6d orbitals (or both) in covalent bonding has been considered for many years, although the amount of covalent character in the highly polar actinide-element bond is still a subject of debate. In contrast to the chemically familiar orbital overlap concept of covalency in light elements, covalency in actinides may be enhanced by small energy differences between metal and ligand orbitals (Su et al., 2018;Ingram et al., 2008;Tanti et al., 2018;Gregson et al., 2016;Kerridge, 2017;Kaltsoyannis, 2013Kaltsoyannis, , 2016Kaltsoyannis, , 2018Jones et al., 2013;Kirker & Kaltsoyannis, 2011;Prodan et al., 2007;Neidig et al., 2013;Kelley et al., 2017). This observation has provoked further debate as to whether 'genuine' covalency ought to increase or decrease across the actinide series. Further complicating matters, experiment and theory indicate that the filled 6s and 6p orbitals also respond to bonding in actinide-element compounds, and are commonly described as pseudo-core orbitals (Boring et al., 1974;Onoe et al., 1993;Pyykkö, 1988;Neidig et al., 2013;Denning, 2007). Experimental insight into f-element bonding has predominately been obtained spectroscopically, for example, by means of UV-Vis (absorption and fluorescence) (Jorgensen et al., 1963), photoelectron (Dau et al., 2012), Mö ssbauer (Kalvius, 1986) and, more recently, X-ray absorption spectroscopy (Minasian et al., 2012;Jollet et al., 1997), although such measurements involve a convolution of the ground-and excited-states. Charge-density studies and topological analysis of the total electron density in the ground state provide a unique opportunity to study different bonding interactions simultaneously and result in quantitative characterization of the bonds present in a crystal, and the combined experimental and theoretical approach provides a more comprehensive understanding of such kinds of interactions. Accurate charge-density experiments are often in remarkable agreement with theoretical predictions in the case of typical light-atom structures and are well represented in the literature (Gianopoulos et al., 2016;Chua et al., 2017b,a;Tidey et al., 2017). Such studies were regarded as 'almost impossible on heavy-atom systems' (Lander et al., 1986) only 30 years ago, but since then the area of experimental X-ray charge-density investigations has greatly expanded and improved, benefitting from considerable advances in computing, diffraction instrumentation and experimental methodologies. In our opinion, these improvements have made charge-density studies on heavy-atom systems an attractive opportunity for further development of the field, although such experiments are still extremely challenging. Thus, recent work in our laboratory concerning Cs 2 UO 2 Cl 4 (Zhurov et al., 2011a,b) and [PPh 4 ] [UF 6 ]  has demonstrated that accurate charge densities of uranium complexes can be obtained from an in-house X-ray diffraction system at low temperatures (20 K). The adequacy of these results can be confirmed by topological analysis of the resulting total ground-state electron density and by comparison with results derived from theoretical charge densities in the context of Bader's Quantum Theory of Atoms in Molecules (QTAIM) (Bader, 1994). Experimentally determined bond descriptors are in fair agreement with those obtained from high-level theoretical calculations (Tanti et al., 2018;Gianopoulos et al., 2017b;Wellington et al., 2016), although theoretical results may show a dependence on the methods applied (Schreckenbach & Shamov, 2010).
Our goal herein is to demonstrate that experimental chargedensity studies of heavy-atom systems are worthwhile challenges and provide opportunities to refine our understanding of actinide-element bonding and inform theoretical developments in the area. The following covers results from three areas of our approach presented as one cohesive story. In short, we will describe (i) improvements in our data collection and reduction strategy, (ii) modifications to the traditional Hansen-Coppens multipolar formalism to increase its flexibility when applied to actinide elements and (iii) the outcome of these efforts for [PPh 4 ] [UF 6 ]. In this context, the following will relate our observations and opinions developed over the last several years, as well as difficulties arising well before the stage of model building. Some of these concerns, while well known, still merit discussion and require careful correction to the primary data. Such effects include application of a flood field adjustment, corrections for strong absorption, harmonic contamination of monochromated X-rays Kirschbaum et al., 1997), oblique incidence correction, multiple scattering and resolution-dependent radiation damage to the sample. Although heavy-atom systems are significantly more complicated when compared with lightatom structures, we believe that combined theoretical and experimental charge-density studies provide attractive opportunities for understanding the nature of actinideelement bonding.

Collection and treatment of primary data
Although the experimental details and methods applied for the processing of the primary data have been previously described, in the following, we have reprocessed the original diffraction data obtained for [PPh 4 ][UF 6 ] at 20 K , and tested different multipole refinement protocols using both XD2006 (Volkov et al., 2006) and MoPro (Jelsch et al., 2005). We have thus obtained a significantly improved result, as well as informing the charge density community on potential protocols to follow and pitfalls to avoid. Thus we will expand on what we have found to be important considerations for obtaining charge-density quality data sets for heavy-atom containing systems, keeping in mind that the quality of those data must be significantly higher than for a light-atom charge density analysis.
To begin with, we prefer to use a 0.5 mm-wide collimator which provides no less than 95% of the maximum beam intensity for the 0.3 mm area at the crystal position rather than a sharply focused one (Zhurov et al., 2008). This removes the need for a correction for inhomogeneous illumination of the sample. Carefully optimized data integration, which takes into account 1 / 2 splitting, is essential. Developed for this purpose, the program VIIPP (Zhurov et al., 2008(Zhurov et al., , 1999 with an image plate floodfield correction, ability to linearize detector response, background and reflection profiles averaged over the whole data set, and rejection of partial or overlapped reflections provides superior data. Careful alignment of the detector system eliminates a significant number of possible errors. For example, the Rapid II detector utilizes separate photomultiplier tubes (PMTs) for strong (PMT2) and weak data (PMT1). We have found that the signal-to-noise ratio can be increased by optimizing the tube voltage for PMT1 (usually by decreasing it), followed by appropriate adjustment of the tube voltage for PMT2. Reduced measured intensities in the case of decreasing the PMT1 voltage can be compensated for by increasing exposure time. For data accuracy, the gain in signal-to-noise for weak data is more important than the decrease of measured intensities. We also have found that non-linearity of the PMTs response should be taken into account. In our case, we have measured the intensity of the direct beam with attenuators at several exposure times providing about 1000 counts s À1 at the beam center. A slight non-linearity was observed for PMT1 ( Fig. 1); whereas the PMT2 response appeared to be linear within the accuracy of the measurements. We have also observed that over time the response of the PMTs is subject to slight drifting. Periodic recalibration of our PMTs is helping to ensure that the highest quality data are obtained.
In addition to the PMT recalibration, we have also found that an oblique incidence effect (Wu et al., 2002;Zaleski et al., 1998), although small for the image plate used in RAPID detectors, still has to be accounted for. We have measured the research papers 896 Gianopoulos et al. Strategies for data reduction and model building value of this correction by equivalents comparison using a spherical ruby crystal. The magnitude of the deviation varies with source wavelength and its behavior is shown in Fig. 2 for Mo K and Ag K radiation. The correction has been incorporated into our integration software VIIPP. The effect is smaller for the less penetrating Mo K radiation (deviation at IP edges is up to $6% on average), R merge still improves from $2.0 to $1.5% for the spherical ruby test crystal. For Ag K, the improvement in R merge is doubled, from $2.5 to $1.5%.
There are several sources which contribute to over-estimation of weak data, such as multiple diffraction, contamination with cosmic radiation and electronic noise spikes. Application of I/(I)-based cutoff prior to averaging also elevates weak data. For this reason, we use all measured equivalents including negative intensities and only apply an I/(I) cutoff after merging. High data redundancy lowers statistical counting errors but it is often difficult to obtain redundancy as high as desired for the highest-order observations, which are inherently weaker. Plots of |F obs |/|F calc | for high-angle reflections demonstrate that the errors associated with weak high-order data essentially average out; however, we find that inclusion of data with I < 3(I) introduces undesirable noise in the experimental electron density, as judged on the basis of Fourier difference maps. Including data with I < 3(I) would be justifiable in the case of zero 'average' noise, unfortunately, in our system -due to the aforementioned reasons -this average is positive and produces bias in the refined scale factor, where slight errors in the scale factor are known to introduce bias in the charge-density model (Becker, 1977;Rees, 1978;Stevens & Coppens, 1975). Despite rejecting some unique data due to applying the sigma cut, the data-to-parameter ratio is still very high owing to the high resolution of the measurement Àðsin =Þ max = 1.3 Å À1 (Mo) or 1.7 Å À1 (Ag).
The use of helium cooling (20 K) provides significant improvement to data quality by removing thermal diffuse scattering, increasing intensities and accuracy of high-order data, and reducing the amplitude of both harmonic and anharmonic atomic displacements more prevalent at higher temperatures (Larsen, 1995). Proper separation of thermal motion from the aspherical electron density is necessary to obtain a meaningful multipole model. While the boost in scattering power is considerable at helium temperatures relative to nitrogen cooling, the added benefit to the reduction of thermal motion due to lower temperature cannot be ignored. Moreover, the increased scattering power, particularly for high-angle data, allows for shorter frame times and experiments, which enables the measurement of a highly redundant data set fairly quickly (1-2 days) using a large area detector (Zhurov et al., 2008).
Crystal size and shape are also extremely important. A small ($100 mm) isotropically shaped crystal would be ideal in order to minimize the effects of anisotropic absorption and extinction, while still having enough scattering power to keep frame lengths reasonable. In the case of [PPh 4 ][UF 6 ], the crystal was a prism of dimensions 0.17 Â 0.14 Â 0.09 mm.
A careful absorption correction is necessary for heavy-atom systems and to that end we have obtained our best results using a numerical absorption correction as previously described and implemented in the program CCDABS (Zhurov & Tanaka, 2003). Briefly, this method involves photographing the crystal with a high-resolution digital camera. Photographs are taken at several degree increments about the ' axis (typically 5 ), providing 72 or more still  Oblique incidence effect for RAXIS RAPID detector with (a) Mo K and (b) Ag K radiaton, determined from a spherical ruby crystal and plotted as a function of vertical distance from the equatorial plane in pixels, based on intensity ratios of equivalents to the one within AE100 pixels from the equatorial plane.

Figure 1
Calibration curve for the response of PMT1 and PMT2. The full calibration curve is shown in the inset.
images of the crystal covering a complete 360 rotation. The direct space outer contour of the crystal is then determined from pairs of mirrored images at the positions ' and ' + 180 . The outlining contours at these positions will therefore be identical, other than mirrored. Both images are used simultaneously for the boundary determination at angles of ' and ' + 180 to improve accuracy as it may be difficult to identify from only one image of the pair due to shadows and glare at certain crystal positions. From the obtained set of contours, we then reconstruct a series of limiting planes which determine the crystal shape. The recovered three-dimensional shape is then used for accurate numerical absorption correction. This method has proven to be extremely useful for crystals of arbitrary shape or when faces cannot be clearly seen, making it troublesome to apply absorption corrections based on faceindexing. For the example of [PPh 4 ][UF 6 ] and prior to other corrections described below, R merge improved from 3.23 to 2.87% when comparing an empirical method, as implemented in the program SORTAV (Blessing, 1995(Blessing, , 1997, with numerical absorption correction performed with CCDABS. Following scaling and merging of the absorption corrected data, several additional corrections may be necessary including corrections for harmonic contamination of the 'monochromatic' X-ray beam Kirschbaum et al., 1997) and radiation damage to the sample. The presence of harmonic contamination is a machine-specific property. All crystals with small absorption will be affected identically. In the case of strongly absorbing crystals, a secondary effect will also be observed due to the difference of the absorption and anomalous scattering coefficients at and different harmonics of . Although the correction for harmonic contamination is generally small, it is especially important for weak low-angle reflections, e.g. see the case of triaminotrinitrobenzene, wherein inclusion of the over-estimated, uncorrected 001 and 003 reflections in the original data set were responsible for unreasonable polarization of the electron density perpendicular to the molecular plane .
Based on our experiences to date, many crystals of uranium compounds are prone to radiation damage even at 20 K, and the [PPh 4 ][UF 6 ] system is no exception. In some cases where the radiation damage is severe, it is worthwhile to re-measure with another crystal, as we have noticed that different crystals from the same sample tend to decay at different rates depending on crystal quality. Radiation damage causes weakening of intensities which increases with resolution, the effect being similar to an increase of thermal motion. Diagnosis of severe radiation damage may be obvious, for example, color change (see Fig. 3), but can also be diagnosed by comparison of the scaling of low-resolution and high-resolution data separately. If radiation damage is present, the scaling factors of the high-resolution data will increase considerably with frame number, but for the low resolution they will only increase slightly. A traditional way of dealing with radiation damage is introducing an overall B-factor for each frame in addition to the scale factor; this has proven useful in macromolecular crystallography where radiation damage is a common problem (Lomb et al., 2011). These values can be determined by refinement during the scaling and merging process or by calculating the ratios of the scaling factors for the high-and low-resolution data, assuming that the distribution of reflections participating in scaling is similar for each frame.
The accuracy of high-order data influences the magnitude of high-frequency noise in charge density maps and therefore how well the features near nuclear positions are resolved, especially for very heavy atoms, owing to their high charge density compared with light atoms. The decay may affect intensities of high-resolution data by an order of magnitude. In this case, the reconstruction of the undisturbed electron density is very complicated. Strongly decaying systems are thus likely to be poor candidates for charge-density studies with heavy atoms, and indeed we have abandoned several compounds/datasets as radiation damage was too severe. Typically, we try to work with crystals that exhibit narrow, well formed and symmetrical diffraction peaks. This is generally observed for high-quality crystals which we have found to be less prone to decay than those of lower quality. Unfortunately, the magnitude of decay in many cases can only be estimated after data are collected.
Another significant source of error in measured intensities is the presence of multiple diffraction (MD). In the case of area detectors with limited degrees of goniostat freedom, the MD effect is unavoidable. While such effects were first described between 1920 and 1940 by Berg (1926), Wagner (1923) and Renninger (1937), contamination of primary data by MD is rarely accounted for (Zhurova et al., 1999;Tanaka et al., 1994;Sakakura et al., 2014). However, as noted by Speakman and others, MD ought not to be ignored when extremely accurate data are required (Speakman, 1965), i.e. for charge-density studies and especially for heavy-atom containing systems (Tanaka et al., 1994). Multiple diffraction events can attenuate or amplify the intensity of a measured reflection, known as aufhellung and umweganregung, respectively. In general, MD contaminates strong and weak reflections in different ways, namely, the important MD events attenuate the intensity (aufhellung) of strong data while increasing the intensity (umweganregung) of weak reflections (Tanaka et al., 1994). In the case of a four-circle diffractometer equipped with a point detector, a protocol to avoid the   collection of strongly contaminated data can be developed (Tanaka et al., 1994). However, implementation of such a routine for imaging plate data is a considerable challenge.
In order to analyze the contamination of measured intensities, we carefully examined images and equivalent reflections measured with sufficient redundancy, as previously suggested (Macchi et al., 2015;Otwinowski & Minor, 1997). A simple graphical tool visualizing the distribution of equivalents for each independent reflection, sorted by intensity (Fig. 4), can be very useful for identifying outliers. Manual inspection of the data, though time consuming, allows for more confident outlier rejection from the dataset. In general, equivalents rejected in such a manner amount to about 1% of the total measured reflections. Even though the improvement of R merge is slight (a few hundredths of a percent), the decrease in the noise level of residual density maps is noticeable. In the case of [PPh 4 ][UF 6 ], decay correction improved R merge by $1% after absorption correction, followed by manual rejection of 417 out of 75 976 (0.55%) measured intensities, which reduced R merge by $0.08% to yield a final R merge of 1.68% for all data.
As a rule of thumb, we try to work with a dataset that has R merge < 1.75%. Small differences in R merge may not appear significant, but datasets with R merge > 2% are noisier and more challenging to model to the point that we would likely abandon charge-density modeling on a dataset with R merge > 3%. While the incremental improvements described above have proven to be necessary for our heavy elements program, we note that our small-molecule datasets have also benefitted from the careful attention to improving our methodology. For example, we recently obtained R 1 = 0.68% for a multipole refinement on the energetic small-molecule salt TKX-50 at 20 K [monoclinic P2 1 /c; (sin/) max = 1.33 Å À1 ; data/parameter = 12.34; Tidey et al., 2017].

Modifications of the traditional Hansen-Coppens formalism
In the Hansen-Coppens multipolar formalism, the electron density of a crystal is modeled as the sum of pseudo-atoms, which traditionally employ three parts [equation (1)], namely (i) a spherical frozen core, (ii) a spherical valence part with P v and s parameters describing the valence charge and expansion/contraction of this term, and (iii) the aspherical portion of the valence density represented as a sum of multipolar terms with the appropriate expansion/contraction coefficients ( l ).
Generally, the difference between such a representation and a spherical atom model reflects a small but chemically meaningful perturbation. In the case of light-atom structures, the traditional multipolar formalism has been extremely successful and can accurately recover total electron densities based on the excellent agreement between experimentally and theoretically derived properties. Nevertheless, while for light elements such as carbon, a single mixed radial term (R l ) is sufficient to describe the hybridized 2s2p orbital; for heavy elements with very different orbitals involved in the valence shell, a single mixed radial term representing them all simultaneously is inappropriate. For example, the ground-state configuration of uranium is [Rn]5f 3 6d 1 7s 2 , spanning three principal quantum numbers which have very different radial behavior (Fig. 5). Moreover, there is significant overlap between the radial distribution of the valence terms and some of the core density. It is therefore unsurprising that (i) neglecting perturbation of the relevant core electron distribution when using a large frozen core and (ii) the use of a single-weighted average radial function to model the aspherical portion of the electron density do not provide a satisfactory description of the experimental electron density in the case of heavy-element compounds. In our previous work (Zhurov et al., 2011a,b Non-normalized, one-electron single-radial functions for the U atom for the valence levels 7s, 6d and 5f, depicting each separately as well as the weighted average corresponding to the 5f 3 + 6d 1 + 7s 2 configuration. The more contracted 5f radial distribution is plotted according to the right axis, whereas the other more diffuse radial functions are plotted on the left axis. employed an augmented multipole model, following the example described by Farrugia & Senn (2012) by including additional aspherical terms based on four to five radial functions representing 5f, 6s, 6p (or averaged 6s + 6p), 6d and 7s orbitals, but keeping the frozen core. We have now redone these refinements for [PPh 4 ][UF 6 ] using the data reduction protocol described above, as well as using a smaller frozen core and a larger set of radial functions as suggested by theoretical studies (Cao et al., 2003;Odoh & Schreckenbach, 2010;Kü chle et al., 1994). The modified Hansen-Coppens scheme is described below in equation (2) and is similar to those previously employed for heavy-element studies (Batke & Eickerling, 2013;Gianopoulos et al., 2017b;Zhurov et al., 2011a,b). The single weighted-average 'valence' term has been replaced by a summation over several split pseudo-atoms. The parts of the complete atom are one frozen-core + 'valence' component and several no-core 'valence' only pseudo-atoms. This results in a smaller frozen core and the sum of N pseudoatoms which accordingly model the spherical valence and aspherical portions of the heavy-atom electron density. We also note that the multipolar expansion has been expanded to the sixth-order, which is necessary to accurately describe f orbitals.
l¼0 3 l R l ð l rÞ X l m¼0 P lmAE y lmAE ðr=rÞ : Note that the index j = 1 . . . N has been omitted from the part of equation (2) in brackets for clarity, but is actually applied to each term. This problem is also well known among theoreticians, as allelectron calculations are exceedingly costly for very heavy elements. A popular approximate method, which decreases computational requirements, employs relativistic effective core potentials (RECPs). In the RECP method, a transferable pseudopotential for the frozen core is defined up to a certain level of n, and calculations proceed with the rest of the electrons (Kahn et al., 1978). Results from such calculations demonstrate that outer-core electrons play a role in the chemistry of actinides and that they may be considered as 'pseudo-valence' or 'semi-core'. In the case of uranium, it has been shown that, at the minimum, the closed 6s and 6p orbitals ought to be considered in the 'valence' space; however, significant errors may still result from excluding electrons belonging to the n = 5 shell. Among the most preferred RECPs for uranium is the small-core Stuttgart RECP with 60 frozen-core electrons. This method includes all electrons with n ! 5 in the 'valence' space, and has been shown to minimize the significant 'frozen core errors' obtained with the large-core Stuttgart RECP (78 frozen-core electrons) (Schreckenbach & Shamov, 2010). Therefore it is not surprising that the best agreement between the obtained augmented Hansen-Coppens experimental models and theory is observed when the uranium partitioning is defined in a manner similar to the small-core Stuttgart RECP employed in theoretical calculations.
We have thus explored multiple ways to construct the uranium pseudo-atom with the expectation that outer-core orbitals will indeed need to be included in the aspherical portion of the multipolar formalism. Unfortunately, we have not found a perfect solution despite considerable testing, but we have been able to make incremental progress towards improving crystallographic agreement factors as well as flattening the residual density around the uranium nucleus. In the following, we will describe two new models which both employ six pseudo-atom terms for uranium partitioned in two different ways (28 and 60 frozen-core electrons) to explore the possibility of observing polarization in deeper electron shells. These will be compared with our previous communication wherein we chose to model the uranium atom with five pseudo-atoms and a frozen core including up to the 5d level (78 frozen-core electrons). By comparison of these partitioning schemes, we hope to identify which is most preferable for understanding the electronic structure and properties of actinide containing compounds.

Geometry and symmetry
The compound [PPh 4 ][UF 6 ] crystallizes in the tetragonal space group I4, with the uranium atom located on the 4 symmetry element (Fig. 6). The equatorial fluorine(s) (F1) occupy general positions and the axial F2 ligands lie on the twofold axis which is coincident with the 4 element. The phosphorus atom of the cation also lies on a position of 4 symmetry with the independent C 6 H 5 fragment occupying general positions. The high-symmetry sites of the uranium and phosphorus atoms conveniently constrain the crystallographically allowed multipole parameters. ORTEP plot of [PPh 4 ] [UF 6 ] with ellipsoids at the 99% probability level, H atoms have been omitted for clarity. Symmetry operators: a = y, 1 À x, 1 À z; b = 1 À x, 1 À y, z; c = 1 À y, x, À z; d = À0.5 + y, 0.5 À x, 0.5 À z; e = 0.5 À y, 0.5 + x, À0.5 À z; f = À x, 1 À y, z.

Refinement strategy
Refinement of H, C and F atoms followed the traditional Hansen-Coppens scheme. Modeling of carbon atoms utilized a single s and l for all six phenyl carbons and multipole expansion up to and including l = 4 with no symmetry constraints imposed. All five hydrogen atoms were treated with a single s and l (constrained to 1.2), the multipole expansion being carried out using dipolar terms only and bond distances constrained to the average neutron values. The fluorine atoms were modeled with separate s and l parameters with a similarity restraint applied. All symmetryallowed multipole parameters up to and including l = 4 were refined, namely all multipoles for the F1 atoms occupying a general position and those allowed by twofold symmetry for the axial F2 atoms. The modeling strategy for the phosphorus atom required some care to ultimately obtain a negative Laplacian at the critical point for the covalent P-C bond. We have chosen to split the phosphorus atom with an additional component such that the P atom is described with a frozen 1s 2 core and refined P v and s on the 2s 2 2p 6 shell and a separate core-less component for the valence 3s + 3p electrons. With this approach, the expected negative Laplacian at the P-C bond critical point was obtained. It is worth noting that the 2s2p portion only deviates slightly from the expected spherical frozen core (<2.5% for s and P v refines to $7.98 e), while multipole parameters refine to $0 and were finally set to zero. For the 3s3p valence portion of the P atom, s and l were refined as well as all multipole parameters up to and including l = 4, allowed by the 4 symmetry. Employment of anharmonic thermal parameters with a non-split model provided the same result as the split-core approach but utilized significantly more variables.
Following a traditional multipole refinement, the observation of high residual peaks significantly closer to the uranium nucleus than the maxima of the radial distributions of the valence electrons suggests that perturbation of the core electron density should also be accounted for (see Table 3). Therefore the inclusion of 'semi-core' pseudo-atoms is necessary to obtain the best models in the case of heavy-atom systems and is in agreement with the recommendation to avoid large-core pseudopotentials for theoretical calculations (Schreckenbach & Shamov, 2010). In our preferred method, the uranium atom is split into the sum of six pseudo-atoms; (i) a frozen core which contains all electrons up to the 4f shell and a 7s valence term, separate valence terms for (ii) 6d and (iii) 5f, (iv) a weighted average for 6s and 6p pseudo-valence, (v) 5d semi-core and (vi) a weighted average for 5s and 5p. In order to maintain the relative shell structure of the uranium atom in the collective pseudo-atoms, we have chosen to apply two separate kappa restraints (or constraints), i.e. one for the spherical kappa sets and one for the aspherical kappas, for all uranium pseudo-atoms such that the kappa parameters have a tight similarity restraint (or to be identical).
Refinements are started using a block refinement technique, the P v and s terms being refined during early block refinement stages. Generally, multipolar terms are slowly added to the uranium pseudo-atoms, as informed by inspection of the residual density, paying particular attention to the peak heights and their distance from the nucleus, and comparing with the available radial functions (Fig. 5). Avoiding correlations between the diffuse 6d and 7s terms is also a challenge and requires extremely accurate low-order data as can be seen from the scattering curves in Fig. 7. Moreover, significant correlations may appear between similar multipolar density terms, for example, in the present case, between the z-directed (parallel to 4) multipolar functions P lm with m = 0 for different l. Thus, while gradually adding uranium multipolar terms, we initially restrict these to a small set. For example, we first add angular terms appropriate to the orbitals anticipated to be relevant, such as l = 2 angular density functions for p pseudoatoms, l = 4 for d orbitals and l = 6 for f orbitals. The complete multipolar expansion over the entire basis set may be added later as necessary. In the case of [PPh 4 ][UF 6 ], we restrict the multipolar expansion to symmetry-allowed multipoles for 4/m symmetry at the U pseudo-atoms, despite the crystallographic 4 symmetry of uranium. We justify this choice on the basis that deviation from 4/m symmetry is only slight and inclusion of l = odd multipolar terms do not significantly improve the model fit, quality of difference or residual maps, or crystallographic agreement factors. That said, with multipole expansion up to and including l = 6 for uranium, a maximum of seven multipolar parameters can be refined per U pseudo-atom using the approximate 4/m symmetry compared with a maximum of eleven multipole parameters at 4 symmetry.
Despite the overlap of the tails of n = 4 radial density functions with the valence 5f orbital (Fig. S1 of the supporting information), inclusion of split pseudo-atoms for the n = 4 electrons does not dramatically improve refinement, although in some cases it appears to stabilize refinement of the l parameters. We also note that refinement of anharmonic thermal parameters (included in model 1c, see discussion below) stabilizes the refinement of multipolar l parameters for the U pseudo-atoms which, in some cases, tend to drive towards being unrealistically contracted ( l ) 1.5). Correlated  Scattering factors (normalized to one electron) for the uranium atom and selected subshells plotted against resolution up to sin()/ = 1 Å À1 . spherical valence charges (P v ) for inner shells were restrained to be equal to the number of electrons in that shell. The spherical valence charges for 5f and 6d were freely refined in both 1b and 1c. The charge on the 6s + 6p pseudo-atom was freely refined in the case of 1b and refined to 7.40. Free refinement of the 6s + 6p charge was also attempted for model 1c and provided similar results during block refinement, but refined to a population greater than eight during full refinement. Ultimately the population was restrained to its neutral value for model 1c.
Treatment of the density due to the extremely diffuse 7s orbital is not straightforward. We have experimented with completely removing the 7s portion, using an averaged term for the diffuse 6d and 7s density functions, and have also tried to keep them separate but employ a restraint on the 7s valence charge. As depicted in Fig. 7, the major portion of information concerning the 7s density of uranium is expected to be contained in the first 16 low-order reflections (sin/ < $0.15 Å À1 , assuming s and l = 1), and the first 36 reflections in the case of 6d (sin/ < $0.20 Å À1 , assuming s and l = 1). Moreover, spatial overlap of the 7s radial function with the nuclear position(s) of the surrounding ligands results in strong correlation between the spherical 7s valence charge and the ligand P v parameters. We concluded that the best treatment of the population of the 7s pseudo-atom was to restrain it to a small value which is consistent with theoretical results (Boring & Wood, 1979;Boring et al., 1974;De Jong & Nieuwpoort, 1996;Hay et al., 1979;Straka et al., 2003).
When evaluating the quality of the model during the refinement, aside from lower R-factors, it is important to monitor the residual density, and normal probability and scale factor plots, as implemented in the program DRKplot (Stash, 2007), and the reasonableness of model parameters. In our experience, following the response observed in the scale factor plots and residual density provides suggestions for the values of well fitted model parameters. At times these may need to be frozen during certain stages of block refinement in order to stabilize the refinement of correlated parameters. Ultimately, the fact that these parameters can safely be refined at later stages suggests the adequacy of model parameters obtained by such a method. Nevertheless, the least-squares surface is typically quite flat and it is difficult to obtain stable refinements without the use of a damp parameter or chemically sensible restraints. Although this is not ideal, we are generally unable to obtain stable models without damping at this time.

Comparison of models
With this general strategy we have obtained several well fitting models, see Table 1. The models presented below (1a, 1b and 1c) were designed in order to compare the effect of the uranium atom frozen core size analogous to theoretical studies, which have shown that large-core RECPs may introduce significant 'frozen core' (FC) errors in comparison with the small-core and all electron treatments (Odoh & Schreckenbach, 2010). Specifically, model 1a includes 78 electrons in the frozen core, model 1b includes 60 FC electrons and 1c includes 28 electrons, additional details about partitioning of the aspherical portions are given in Table 1. Variations between R-factors, residual density maxima, scale factor plots and residual density (Table 1 and Figs. S3-S6) were taken into account when deciding which model is most preferred. For this reason we have also compared the results of topological analysis of the total electron density for each model with those from theory  [B3LYP (Becke, 1993a,b); F: all-electron 6-31 g(d 0 , p 0 ); U: small-core Stuttgart RSC 1997 RECP (Dolg et al., 1993)].

Topology of the total electron density
The reconstructed electron densities have been evaluated and compared based on the properties of the electron density and its derivatives at the (3, À1) bond critical points (bcp) following Bader's Quantum Theory of Atoms in Molecules (Bader, 1994). Further qualitative comparisons between the models were also made by comparing maps of the deformation densities and the Laplacian of the total electron density; these will be discussed in turn below. There is fair agreement between the models' properties and those predicted by theory (Table 2), whereas the corresponding properties obtained from a spherical atom model are quite different.
The topological properties at the U-F bcps for all three models (Table 2) demonstrate slight differences of and r 2 . Augmented Hansen-Coppens models 1b and 1c, which include aspherical modeling of deeper lying core electrons, are in closer agreement with our theoretical results compared with model 1a. Moreover, previous theoretical studies have also indicated that outer-core electrons do indeed respond to perturbations resulting from metal-ligand interactions, especially in the case of 6p (Fryer-Kanssen & Kerridge, 2018;Denning, 2007;Bartleet et al., 1992;O'Grady & Kaltsoyannis, 2002). The 'best' agreement between the experimental and theoretical (  I > 2(I) R 1 = 0.0082, wR 2 all data = 0.0181, 1.291 and À1.300 e Å À3 Multipole model 1a † 78 e FC: 1s to 5d 1b 60 e FC: 1s to 4f 1c 28 e FC: 1s to 3d Valence pseudo-atoms 6s, 6p, 5f, 6d, 7s 5s + 5p, 5d, 6s + 6p, 5f, 6d, 7s 4s to 4f, 5s to 5d, 6s + 6p, 5f, 6d, 7s Final R indices: R 1 I > 3(I) wR2 all data R 1 = 0.0066, wR 2 = 0.0095 R 1 = 0.0066, wR 2 = 0.0093 R 1 = 0.0065, wR 2 = 0.0084 Á (max and min) sin / < 1.0 Å À1 0.295 and À0.561 e Å À3 0.319 and À0.588 e Å À3 0.298 and À0.519 e Å À3 with model 1b, which has a frozen core analogous to the smallcore RECP used in the theoretical model. In this case, the electron density at the bcp and its Laplacian are within 10% of the theoretical values. On the other hand, model 1c provided the most stable refinement allowing looser restraints. In all cases, the experimental models suggest more covalent character in the U-F bonds than is predicted by theory on the basis of their QTAIM properties. Thus, when comparing bonding descriptors (Espinosa et al., 2002;Gatti, 2005) of the experimental and theoretical models, for the experiment the density at the bcp is higher, the Laplacian is less positive, the total electronic energy density, h, is more negative, the ratio of electron potential and kinetic energy densities |v|/g is larger, and the 'covalence degree' h/ is more negative. In the case of model 1b, the covalence degree descriptor indicates 7% (U-F2) and 15% (U-F1) greater covalent character than theory. For all models, both experimental and theoretical, topological results place the bonding regime in the so called transit region, where bonding results from the interplay of electrostatic and covalent contributions (Espinosa et al., 2002;Gatti, 2005). Comparison of the deformation densities (DD) obtained from theory and experiment is also informative. As can be seen from Fig. 8, the deformation density around the F atoms is consistent across the models described herein. The axial F2 atoms appear to have more structure, and their DD is characterized by positive regions in the direction of the metal and also directly behind the F atom, opposite to the direction of the U-F bond. The deformation density in the direction of the metal has two maxima (see Fig. 8). A charge depletion is observed near the F2 nuclear position and extends in the direction perpendicular to the U-F bond, suggesting polarization of the electron density consistent with -donor character in the ligand-metal interaction. The equatorial F1 atoms appear to have less structure in the deformation density maps in as much as the maximum values of the DD are not large and any maxima essentially smear together at a low isodensity level. Nevertheless, as in the case of the axial fluorine atoms, there are two maxima directed towards the metal for the equatorial F1 atoms as well.
However, the modeling strategy does affect the structure of the DD around the uranium atom. Comparing models 1a through 1c it can be seen that including asphericity of the outer-core electrons seems to 'resolve' some of the features observed in the DD. For example, the DD around the U atom Deformation density isosurface map for the UF 6 À anion obtained from models 1a-1c. Charge concentrations are depicted in blue and depletions in red. The isodensity surfaces are plotted from the AE0.15 (mesh), 0.30 and 0.45 (most opaque) eÅ À3 levels.

Table 2
Properties of the electron density at the U-F bond critical points.
(eÅ À3 ) is the electron density at the bcp; r 2 (eÅ À5 ) is the Laplacian at the bcp; d (U-cp) (Å ) is the distance from uranium to the bcp in Å ; g (a.u.) is the electron kinetic energy density at the bcp; v (a.u.) is the electron potential energy density at the bcp; h (a.u.) is the total electron energy density at the bcp; |v|/g is the adimensional bonding regime descriptor; h/ (a.u.) is the covalence degree. The electron kinetic energy density g was obtained by the Abramov approximation (Abramov, 1997) and the electron potential energy density v according to the local Virial relationship (Bader, 1994). (for all models) shows eight maxima at a distance of $0.6 Å from the U nuclear position, as expected for electron density belonging to the 5f and/or 6s-6p levels, and resembles what would be expected for the electron density due to a singly occupied 5f orbital (Kaltsoyannis & Bursten, 1995;Boring et al., 1974) (FWHM for the 5f radial function spans from $0.33 to 0.74 Å whereas the FWHM for the averaged 6s-6p radial function spans from $0.49 to 1.02 Å , see Table 3). The 'contrast' between these eight concentrations and charges of other regions improves by modeling asphericity of deeper shells in the outer-core. Maximum regions of charge depletion in the DD around U are oriented along the metal-ligand bonds, as expected from ligand-field theory arguments (Fig. 8).
It is surprising and merits further consideration that there is such a noticeable difference in the deformation density between the axial and equatorial fluorine ligands despite the fact that the tetragonal distortion is very slight (0.011 Å ), although statistically significant. Under closer scrutiny, it is apparent that the crystal surroundings could be partially responsible for these differences (Fig. 9). The axial fluorine atoms participate in five intermolecular interactions for which bond paths were characterized. Of these, four involve F2Á Á ÁH interactions and the fifth involves F2Á Á ÁF2. All of the F2Á Á ÁH interactions are between F2 and H5 on four different symmetry-related cations. There are four intermolecular interactions involving F1 bonding to three different cations. These involve F1Á Á ÁH2 and F1Á Á ÁH6 on the same PPh 4 cation, but on separate phenyl rings, as well as F1Á Á ÁH4 and F1Á Á ÁC3 on different cations. Assuming the validity of the Espinosa relationship (Espinosa et al., 1999(Espinosa et al., , 1998 for all such interactions, the cumulative dissociation energies of these interactions amounts to $23 and $21 kJ mol À1 for F2 and F1 respectively, with each FÁ Á ÁH interaction contributing around 4.5-6.5 kJ mol À1 . The dissociation energies for the FÁ Á ÁC and FÁ Á ÁF interactions were estimated to be less than 4 kJ mol À1 . The QTAIM properties of the intermolecular interactions were consistent across both of the improved models 1b and 1c (Table S1 of the supporting information).

Topology of the Laplacian of the electron density
Topological analysis of the Laplacian of the electron density, r 2 , has been shown to reveal additional information regarding the spatial redistribution of the electron density in molecules. While the qualitative pictures obtained often resemble the charge concentrations and depletions observed in deformation density maps, the Laplacian distribution does not suffer from the requirement of a promolecule reference density. The critical points of the Laplacian distribution correspond to charge concentrations in the valence and core shells, VSCC and CSCC, respectively, and valence and core shell charge depletions, VSCD and CSCD. In the case of lightatom structures, the VSCCs correspond to bonding and lone pair regions and are in remarkable agreement with the basic concepts of the valence shell electron pair repulsion model (VSEPR) of bonding. This correspondence is striking enough that the Laplacian distribution has been proposed to provide a physical basis for the VSEPR model (Popelier, 2000;Gillespie & Robinson, 1996). Moreover, the differentiation between valence and core shell charge redistributions is straightforward based simply on the distances of the maxima of charge concentration from the nuclear position. Analysis of experimental and theoretical Laplacian distributions has also been applied to metal containing molecules and structures. For example, there are eight VSCCs for the d 6 transition metal compounds Cr(CO) 6 (Farrugia & Evans, 2005) and iron disulfide FeS 2 (Schmøkel et al., 2014), which are arranged at the vertices of a cube where the metal-ligand vectors pass through the center of each face. Such features are referred to as ligand opposed charge concentrations and can be rationa- Molecular graph depicting intermolecular interactions around the independent fluorine atoms, only (3, À1) critical points are shown for clarity. F: light blue; U: blue; C: black; P: orange; H: gray; (3, À1) critical points: purple. Symmetry operators: a = y, 1 À x, 1 À z; b = y, 1 À x, Àz; c = 1 À y, x, 1 À z; d = 1 À z, 1 À y, z; e = 0.5 À x, 1.5 À y, À0.5 + z; f = 1.5 À y, 0.5 + x, 0.5 À z. Table 3 Maxima and FWHM for selected single-radial functions, R l , for valence and outer-core electronic levels. lized in the context of a simple ligand-field theory approach wherein the concentration of the metal d electrons reflects avoidance of the ligand charge concentrations. In addition, analysis of Laplacian distributions derived from theoretical methods suggest significant polarization of outer-core orbitals in transition metal complexes (Bader et al., 1998;Batke & Eickerling, 2013).
In the present case, the valence space of uranium spans three principal quantum numbers and is further complicated by the presence of outercore orbitals which may also deform in response to bonding. The charge concentrations observed in the preferred experimental models 1b and 1c are self-consistent, suggesting that asphericity of the n = 4 electrons is not a major contribution to the bonding picture (Fig. 10). In the case of UF 6 À , there are three 'sets' of charge concentrations surrounding the uranium atom all at a distance of $0.38 Å , corresponding to 'valence shell' charge concentrations near the radial maxima of the n = 5 electron density level. The first 'set' contains eight critical points arranged in a cube and resembles the ligand opposed charge concentrations as described above in the example of Cr(CO) 6 (Farrugia & Evans, 2005). However, in the case of UF 6 À the uranium-fluorine vectors do not pass perfectly through the faces of such a cube, but are rotated towards the metal-ligand vectors. In the second 'set' there are four critical points corresponding to charge concentrations which are arranged in a square in the equatorial plane. As with the first 'set', the metalligand vectors are tilted with respect to passing through the midpoint of each edge of the square. The implication of the tilting from the idealized positions as observed in this system is not immediately obvious, although we note that d-orbital tilting has previously been predicted to be physically meaningful in some cases and experimentally confirmed (Deutsch et al., 2011). A similar set of charge concentrations was observed in a theoretical study of VF 5 (i.e. three charge concentrations between equatorial ligands in this trigonal bipyramidal example; . The final set consists of two charge concentrations directed along the U-F2 bond vector. We also note the presence of another (3, +3) critical point (corresponding to charge concentration, although the Laplacian is positive at this critical point) along the bond nearer to the F2 atom and in the vicinity of the U-F2 bcp; interestingly, we do not find the analogous critical points along the U-F1 bond, although the research papers IUCrJ (2019). 6, 895-908 Gianopoulos et al. Strategies for data reduction and model building 905 Envelope diagrams [(c) and (d)] of the negative Laplacian distribution (charge concentration) around the U atom in the UF 6 À anion plotted at the À280 e Å À5 level. The corresponding (3, +3) critical points (maxima of charge concentration) of the Laplacian distribution [(a) and (b)] are depicted as red spheres. The directions of the fluorine ligands are indicated by labels, except for the axial F2 atoms in (b) and (d), where the axial ligands are above and below the plane of the image. When depicted, thermal ellipsoids are drawn at the 99% probability level. Although these distributions were obtained from model 1c, nearly identical distributions are obtained from model 1b. U-cp distances are in the range 0.37-0.38 Å for the 14 nearest cps and 1.075 Å for the distant cp. The F2-cp distance is 1.001 Å .

Figure 11
Envelope diagrams of the negative Laplacian distribution (charge concentration) around the fluorine atoms in the UF 6 À . The isosurfaces are plotted at (a) À160 and (b) À190 e Å À5 levels. Thermal ellipsoids are drawn at the 99% probability level. While these distributions were obtained from model 1c, nearly identical distributions are obtained from model 1b.
Laplacian distribution along the U-F bond paths are similar (Fig. S2).
The Laplacian distribution around the fluorine atoms also merits brief description. Neither fluorine atom demonstrates well defined features in the Laplacian distribution. Thus, until close to the maxima of charge concentration, their Laplacian distributions are nearly spherical, with weak maxima of charge concentration at a distance of $0.3 Å from the nucleus. These results suggest a strong ionic component and an expansion of the charge cloud with respect to a neutral F atom (valence shell radial distribution maximum $0.21 Å ). Nevertheless, when visualized as an isosurface near the maxima of charge concentration the Laplacian distributions around both F1 and F2 resolve into three maxima (Fig. 11). In the case of F2, two maxima are directed towards uranium while the third maximum is directly behind the F2 atom opposite the U-F bond. The maxima of charge concentration around F1 are similar in that two maxima are in the direction of the uranium atom, although with some distortion compared with F2, whereas the third maximum is directed towards H4 on a nearby cation with a maximum -F1-U angle of $120 .
In contrast to the experimental result, there is only one analogous 'set' of critical points corresponding to charge concentrations in the Laplacian distribution around uranium for the theoretical density (Fig. 12). The theoretical CCs are arranged at the vertices of a cube at a distance of $0.85 Å from the uranium atom, corresponding to the radial maxima of the n = 6 level. In this case, the uranium-fluorine bond vectors do pass through the center of each face of the cube formed by the charge concentrations, which is expected from the point of view of ligand-field theory. Nevertheless, it is from this vantage point that disagreement between experiment and theory is most striking. It is of interest to consider whether these differences are due to the choice of RECP employed in the theoretical calculations. The small-core RECPs tend to be favored as the frozen-core 'smooths' features in the 'core' levels. Theoretical results suggest that these 'core-wiggles' may be important especially when trying to understand subtle features and (electronic and derived QTAIM) properties in actinide containing complexes (Odoh & Schreckenbach, 2010). As we have suggested above, it would be of interest to compare the Laplacian distribution and topological properties of simple actinide compounds at various levels of theory and with respect to the core size for RECPs utilized.

Conclusions
Over the past decade we have sought to improve our methodology for the collection of extremely accurate, high-resolution single-crystal diffraction data at cryogenic temperatures as well as improving our data reduction techniques. These efforts have resulted in strategies for data collection, reduction and processing, and modeling, providing the opportunity to study and characterize the electron density in systems containing very heavy elements such as actinides, which were considered 'nearly impossible' only 30 years ago. Moreover, we have shown that experimental results are in very good agreement with those predicted by theory in the case where an augmented Hansen-Coppens scheme was used to model the aspherical electron density of the uranium atom in [PPh 4 ][UF 6 ]. Comparison of QTAIM properties between experimentally and theoretically derived electron densities indicates that accurate modeling requires inclusion of terms to describe deformation of the outer-core electron levels. Analogous results have previously been shown in the case of theoretical calculations. Experiment and theory both indicate that the U-F bond is of mixed character, belonging to the socalled transit region, wherein the interaction can be considered to possess both ionic and covalent character to varying degrees. While the agreement between experiment and theory is strong, there are a few noteworthy differences: first, experimentally derived electron densities suggest $10% more covalent character than is indicated by theory and secondly, the structure of the Laplacian distribution around the heavy atom is markedly different. In the case of experiment, the maxima of charge concentration around uranium (14 total) are in the vicinity expected for the n = 5 level, whereas the maxima of charge concentration for the theoretical result (8 total) correspond to the n = 6 level. It is of interest to consider whether or not this effect results from the RECP treatment of the core electrons for the theoretical model. Finally, differences between the axial and equatorial fluorine atoms are Graph depicting the (3, +3) critical points (green spheres) in the Laplacian as determined from the theoretical model of the UF 6 À ion, shown in three orientations. much more pronounced in the experimental result than suggested by theory. It is likely that crystal packing effects play some role in the differences observed in the axial and equatorial fluorine ligands. Nevertheless, we believe that these differences may provoke discussion between experimental and theoretical chemists and that such interplay will help to refine both experimental and theoretical methodologies. Moreover, these results demonstrate that meaningful results can now be obtained for experimental charge-density studies of systems containing heavy elements such as actinides.