Topology of the Electron Density and of Its Laplacian from Periodic LCAO Calculations on f-Electron Materials: The Case of Cesium Uranyl Chloride

The chemistry of f-electrons in lanthanide and actinide materials is yet to be fully rationalized. Quantum-mechanical simulations can provide useful complementary insight to that obtained from experiments. The quantum theory of atoms in molecules and crystals (QTAIMAC), through thorough topological analysis of the electron density (often complemented by that of its Laplacian) constitutes a general and robust theoretical framework to analyze chemical bonding features from a computed wave function. Here, we present the extension of the Topond module (previously limited to work in terms of s-, p- and d-type basis functions only) of the Crystal program to f- and g-type basis functions within the linear combination of atomic orbitals (LCAO) approach. This allows for an effective QTAIMAC analysis of chemical bonding of lanthanide and actinide materials. The new implemented algorithms are applied to the analysis of the spatial distribution of the electron density and its Laplacian of the cesium uranyl chloride, Cs2UO2Cl4, crystal. Discrepancies between the present theoretical description of chemical bonding and that obtained from a previously reconstructed electron density by experimental X-ray diffraction are illustrated and discussed.


Introduction
The degree of participation and covalency of 4 f and 5 f electrons in the chemical bonding of lanthanide and actinide complexes and materials (including oxides) is an intriguing topic in chemistry [1][2][3][4][5][6][7], with both fundamental and technological implications. Indeed, the renewed interest in the nuclear power industry and nuclear waste management has put the separation process of uranium from lanthanides and other minor actinides under the spotlight. The efficiency and selectivity of the process depends on the respective bond strength [8,9]. In particular, several factors (including strong relativistic effects, weak crystal fields and strong electron correlation) contribute to the occurrence of a composite valence manifold (comprising of 5 f , 6p, 6d and 7s orbital shells) in actinide-ligand bonds [5,[10][11][12].
The quantum theory of atoms in molecules and crystals (QTAIMAC) provides a rigorous formal framework within which multiple aspects of chemical bonding can be analyzed [13,14]. The theory builds on a topological analysis of the electron density ρ(r), which can be either determined by theoretical calculations or measured by experimental X-ray diffraction. In both cases, several technical challenges need to be overcome for effective treatment of f -electrons, to be briefly addressed below. Such an approach is often complemented by an analysis of the topology of the Laplacian of the density ∇ 2 ρ(r), which provides additional information on the spatial distribution of the electrons and, in particular, on the asphericity of (bonded) atoms [15].
So far, the electron density of very few actinide materials could be successfully reconstructed from X-ray diffraction measurements: Th(S 2 PMe 2 ) 4 , cesium uranyl chloride Cs 2 UO 2 Cl 4 and tetraphenyl phosphate uranium hexafluoride [PPh 4 ][UF 6 ] [16][17][18][19][20]. In particular, significant advances have been recently reported by Pinkerton and co-workers on data collection and analysis, which finally allow for the accurate experimental determination of the electron density of actinide materials [21]. Similarly, from a theoretical perspective, the many challenges related to a proper and simultaneous description of relativistic effects of core electrons and electron correlation of valence electrons (including f -type ones) made it possible to perform a QTAIM analysis on top of computed electron densities of actinide molecular complexes only recently [22][23][24][25][26][27][28][29].
In this paper, we present recent formal and software advances in the quantummechanical calculation of the electron density and its Laplacian as well as in their topological analyses, for f -electron materials, within periodic boundary conditions for an effective study of chemical bonding in extended systems. In our methodology, crystalline orbitals are expressed as linear combinations of atomic orbitals (LCAO). Quantum-mechanical calculations are performed with a developmental version of the CRYSTAL program [30,31], where the LCAO approach has recently been extended to g-type basis functions [32,33]. The topological analysis of the electron density ρ(r) and its Laplacian ∇ 2 ρ(r) is performed with a developmental version of the TOPOND module [14,34,35], which can be run in parallel with a high efficiency [36]. Here, we discuss the strategy that we have followed to extend this module (formerly limited to work in terms of s-, pand d-type basis functions only) to f -and g-type basis functions as well. The new implementation has already been successfully applied to the study of chemical bonding in the [PPh 4 ][UF 6 ] crystal, where the topology of the computed density and its Laplacian was found to match that from the experiment very satisfactorily [37]. In this paper, we study the topology of ρ(r) and ∇ 2 ρ(r) of the cesium uranyl chloride Cs 2 UO 2 Cl 4 crystal. The electron density of this system has been reconstructed experimentally and its topology analyzed with the QTAIMAC [18,19]. A detailed theoretical investigation of the electron distribution and chemical bonding in the [UO 2 Cl 4 ] 2− molecular fragment has been performed, which described the weakening of the U-O uranyl bond due to the equatorial ligands [38]. The computed values of ρ(r) and ∇ 2 ρ(r) at the U-Cl bond critical points were in excellent agreement with those from the experiments while a somewhat looser agreement was reported for the U-O bond. A subsequent study investigated whether the observed discrepancy could be due to missing environmental effects in the simulations and concluded that these are very small and therefore cannot explain the observed differences [39].

Computational Details
Calculations on the molecular fragment [UO 2 Cl 4 ] 2− and extended periodic crystal Cs 2 UO 2 Cl 4 are performed with the CRYSTAL program by use of the global hybrid B3LYP exchange-correlation functional [40] of the density functional theory (DFT). In the periodic calculations, reciprocal space is sampled on a regular 6×6×6 Monkhorst-Pack grid, corresponding to 68 k-points in the symmetry-irreducible Brillouin zone. Scalar relativistic effects of U must be accounted for [27][28][29]41] and, here, are described by use of smallcore effective pseudopotentials, ECP60MDF (with 60 electrons in the core for U) [42,43]. The valence of U is described by a (10s9p7d5 f 1g)/[10s9p7d5 f 1g] basis set: we indicate within round brackets the number of Gaussian primitive functions used for the various angular quantum numbers and within square brackets the number of shells in which they are contracted. In this case, we use a fully uncontracted basis. With respect to the original basis set optimized for molecular calculations, some very diffuse exponents have been removed (crucially the most diffuse p-type exponent) that were causing linear dependencies in the periodic calculations. While the program has recently been extended to the treatment of spin-orbit coupling [44][45][46][47], this relativistic effect is disregarded here. This is because, while making the calculations significantly more demanding, it has been previously shown to induce minor changes to chemical bonding for such systems [48]. Oxygen and chlorine are described by molecular def2-TZVP basis sets of (11s6p2d1 f )/[5s3p2d1 f ] and (14s9p3d1 f )/[5s5p2d1 f ] type, respectively [49]. For Cs, Hay-Wadt small-core pseudopotentials are used [50] in combination with a (4s4p1d)/[2sp1d] valence.

The Implementation
In this section, we present some formal aspects related to the implemented algorithms for the extension to f -and g-type basis functions of the topological analysis of the electron density and its Laplacian in the TOPOND program module.

Basis Functions
In the periodic LCAO framework within the CRYSTAL and TOPOND programs, the oneelectron crystalline orbitals (COs) are expressed as a linear combination of Bloch functions φ µ (r, k), depending on the spatial coordinates r of an electron and its wave vector k. A Bloch function is written, in turn, as a Fourier transform of the associated atomic orbital (AO) χ µ (r): where the sum extends to the infinite set of lattice vectors g and A µ is the position of the atom on which χ µ is centered in the reference cell. In CRYSTAL, the AOs are expressed as a linear combination of normalized real-solid spherical harmonic Gaussian-type functions (RSSH-GTF) R n,l,m (r; α λ j ): where: and . . . is the Euclidean norm. In Equation (2), d λ j and α λ j are fixed linear coefficients and exponents associated to a shell λ. In the context of the shell strategy, AOs with similar quantum numbers n and l but different m (e.g., the 2p x , 2p y and 2p z AOs) share a common set of d λ j and α λ j . Thus, formally, we may write λ = λ(n, l), although the dependence of λ on the set of quantum numbers has been suppressed in Equations (2) and (3) to simplify the notation. N λ , N l α λ j and N l,m are normalization constants. The shell normalization constant N λ reads: The m-dependent normalization constant N l,m reads: Finally, in Equation (3) X n,l,m (r) is a RSSH, which may be expressed as a homogeneous polynomial of the Cartesian components r x , r y and r z of the electron coordinates: in which D l,m (t, u, v) are linear coefficients and the sum in Equation (7) extends to all integer triplets t, u, v whose values are constrained by the equality t + u + v = 2n + l.
In the CRYSTAL and TOPOND programs, the AOs in the basis set are defined using only RSSH-GTFs with null principal quantum number, n = 0. However, RSSH-GTFs with n = 0 are still required in the CRYSTAL program because they are used as auxiliary functions for the evaluation of electronic kinetic energy integrals.

The Electron Density and Its Derivatives
The QTAIMAC is based on a topological analysis of the electron density ρ(r) (a scalar field), which requires the evaluation of its first-and second-order derivatives w.r.t. spatial coordinates in the definition of the gradient vector and Hessian matrix: The trace of the Hessian matrix defines the Laplacian of the density (another scalar field): The description of chemical bonding features passes through the identification of critical points (CPs) of the density, defined as those points in space r CP where the gradient vector of the density vanishes ∇ρ(r CP ) = 0. CPs can be classified into different types according to the signs of the corresponding eigenvalues of the Hessian matrix. Further information on the spatial distribution of the electrons can be obtained from a topological analysis of the Laplacian, whose positive and negative values correspond to regions of (relative) charge concentration and depletion, respectively. Valence shell charge concentrations (VSCCs) are particularly relevant to the rationalization of chemical bonding, and can be analyzed in terms of CPs of the Laplacian of type (3, +3), i.e., minima. It is often preferable to flip the sign of the Laplacian and work in terms of the following function: so that positive values and maxima of this function correspond to regions of charge concentration. Critical points of L(r) can be found and characterized by analyzing its firstand second-order derivatives (involving third-and fourth-order derivatives of the electron density, respectively). Thus, at the heart of the TOPOND program lies the evaluation of the electron density ρ(r) and of its derivatives w.r.t. Cartesian components of the spatial coordinate. In the following, i, j, k, l ∈ [x, y, z] indices will be used to refer to Cartesian components. The density is defined as follows in the basis of AOs, using the shorthand notation r µ = r − A µ and r g ν = r − A ν − g: where P g µν is the direct-space one-electron density matrix. For the first derivatives: For the second derivatives: For the third derivatives: and for the fourth derivatives: From inspection of the equations above, we see that the evaluation of the electron density and/or of its derivatives boils down to the efficient evaluation of pair products of two AOs and/or pair products of their derivatives. Bearing in mind the expansion in Equations (2) and (3), the algorithms implemented in the TOPOND module rely on the efficient evaluation of pair products of RSSH-GTFs and/or of their derivatives.

Old Strategy Based on the Expansion in Hermite Gaussian Type Functions
The previously existing strategy implemented in the TOPOND module by Victor R. Saunders (limited to s, p and d AOs) was based on the expansion of any RSSH-GTF pair product into an auxiliary basis [51]: (16) where the E[. . . ] are linear coefficients and the Λs are Hermite Gaussian-type functions (HGTF): In Equation (16), the values of γ and P are obtained from the Gaussian product theorem (i.e., γ = α + β and P = (αA + βB)/γ and the sum over the indices t, u, v is extended to all integer triplets that satisfy the inequality t + u + v ≤ 2n + 2ñ + l +l. The E[. . . ] coefficients can be conveniently generated from recurrence relations [51]. Moreover, once these coefficients are generated for a product pair of RSSH-GTFs with a given set of quantum numbers, they can be re-used to generate the expressions for those product pairs of RSSH-GTFs with lower quantum numbers but involving derivatives, which results in a general and efficient algorithm. However, let us note that the recurrence relations were never coded as such and instead the explicit expressions for all the relevant E[. . . ] coefficients were pre-determined and hard-coded. While making the implementation very efficient, this strategy produced a code that could not be easily extended beyond d-type AOs. For this reason, here we have followed an alternative strategy for the generalization of the code to f -and g-type AOs (i.e., with l = 3, 4).

New Strategy Based on a Direct Evaluation of the Electron Density in the RSSH-GTF Basis
Instead of evaluating the electron density and its derivatives through an expansion of pair products of RSSH-GTFs in HGTFs as in Equation (16), we have devised a new strategy that works for l = 0, 1, 2, 3, 4 (i.e., for s-, p-, d-, f -, and g-type AOs). The expressions for the AO pairs and their derivatives (up to order 4), required for evaluating the electron density and its derivatives in Equation (11), are directly evaluated in the RSSH-GTF basis, following Equation (2).
In a first step, we pre-calculated the symbolic expressions for the RSSH-GTF pairs from the computer algebra system for symbolic computation available in MATLAB by exploiting the following two recurrence relations [52]: X n,l+1,±(l+1) (r) = (2l + 1) r x X n,l,±l (r) ∓ r y X n,l,∓l (r) , and: Equations (18) and (19) allow us to generate the RSSHs for the entire set of l and m values for a given n (here, we limit our considerations to n = 0). The starting point for the recurrence is with the s-type RSSH for which X 0,0,0 = 1 and, in connection with Equation (18), the convention X 0,0,−0 = 0 is understood. The algorithm then proceeds with the calculation of X 0,l+1,l+1 and X 0,l+1,−(l+1) through Equation (18) for l = 1, 2, 3, 4. The RSSH X 0,l,m with m ∈ [−l + 1, l − 1] are then generated for successively higher l. The steps of the used algorithm are summarized in the following, where quantum numbers in bold are being increased: (1)X 0,0,0 = 1 (2)X 0,l,l (3)X 0,l,−l (4)X 0,1,m (5)X 0,l,m Once the full set of RSSH-GTF pairs were generated, the expressions for the derivatives of the RSSH-GTF pairs were determined by symbolic differentiation. Finally, we have generated FORTRAN routines containing the symbolic expressions for the (derivatives of the) RSSH-GTF pairs.
The new code has the obvious advantage of being generally applicable to basis sets with AOs of s-, p-, d-, f -, and g-type compared to the old one that was limited to work with s-, pand d-type AOs. At the same time, it is found to be less efficient than the original one. For this reason, the old code is active by default for those systems where f -, and g-type AOs are not included in the basis set, with the new code coming into play only when needed (i.e., in the presence of f -, and g-type AOs in the basis set).

The Electron Localization Function
Besides ρ(r) and L(r), another scalar field which is commonly analyzed to extract useful information on the electron distribution is the electron localization function η(r) (ELF). Roughly speaking, the ELF function provides a (local) measure of electron localization in a system relative to that of an electron gas having the same electron density ρ(r) of the system [14]. The various, though largely compatible, definitions of η(r) share the general formula: where c(r) is the relevant kernel of the ELF and where η is scaled relative to c to bound its values between 0 and 1. ELF was first introduced by Becke and Edgecombe (BE) using as a kernel, c BE (r), the ratio of the curvatures of the spherically averaged samespin conditional pair densities for the system under study and for the corresponding (i.e., same density) uniform electron gas [53]. Later on, Savin et al. proposed a more general and simpler definition, c Savin (r) = t p (r)/t ph (r), where t p (r) is the Pauli kinetic energy density of the system and t ph (r) that of the uniform electron gas with the same density [54]. The Pauli kinetic energy density represents the local increase of the kinetic energy due to the redistribution of the electrons caused by the Pauli principle; that is, the local kinetic energy excess relative to the kinetic energy density of a bosonic system. Use of either c BE (r) or c Savin (r) leads to formally equivalent ELF expressions, within the Hartree-Fock approximation originally adopted by BE. The interpretation due to Savin et al. has, however, a great conceptual advantage: regardless of the kind of adopted wave function, the same expression for the ELF is retained. In addition, being based on the kinetic energy density rather than on the electron pair density, the ELF can also easily be evaluated within the DFT, as in our case.

Results and Discussion
Cesium uranyl chloride, Cs 2 UO 2 Cl 4 , crystallizes in a monoclinic lattice with space group C2/m, whose atomic structure is shown in Figure 1. Each U atom forms two symmetry-equivalent bonds with O atoms (bond length of 1.776 Å in the experimental geometry) as well as four symmetry-equivalent bonds with Cl equatorial ligands (bond length of 2.670 Å). Each Cs atom is connected with eight Cl atoms (four symmetry-independent pairs with bond lengths in the range 3.502-3.624 Å) and one O atom (bond length of 3.259 Å). For a better comparison with the experimental electron density [18,19], we have performed our quantum-mechanical simulations on the experimental geometry. Calculations are performed both on the actual periodic model (here referred to as cry-Cs 2 UO 2 Cl 4 ) and on an extracted molecular fragment [UO 2 Cl 4 ] 2− to quantify the effect of the crystalline environment on the chemical bonding features of the U atom.  Table 1 reports computed atomic charges for the different species in the two systems (the extended crystal and the molecular fragment). Atomic charges are computed according to two schemes: a simple Mulliken orbital partitioning (M) and Bader's integration of the electron density over atomic basins (B). The two approaches provide similar values for Cl and Cs while the QTAIMAC systematically results in larger charges for U and O, both in the molecular fragment and in the crystal. The computed atomic charges of the molecular fragment are in excellent agreement with those previously reported by Vallet et al. [38]. From inspection of the computed QTAIMAC charges on the extended crystal, we find that the predicted electronic structure of the system is very close to the limit ionic [Cs 2 ] 2+ [UO 2 Cl 4 ] 2− one. When passing from the isolated molecular fragment to the actual periodic crystal, the most affected atomic charges are those of the U and O atoms, which both increase in absolute value. A closer inspection of the computed atomic charges already suggests a stronger polar character of the U-Cl bonds than the U-O ones. Indeed, each Cl is found to host an extra 0.74 electrons with respect to its neutral state for a total of 0.74 × 4 = 2.96 extra electrons hosted by the four Cl atoms surrounding the U atom. This value almost perfectly matches the atomic charge of U (+2.94), thus suggesting that the U atom donates about three electrons to the four Cl atoms while it participates in a much more covalent interaction with O. The agreement between our computed charges and those derived by the experiment is remarkable and corroborates this picture. Chemical bonding features of cesium uranyl chloride have been analyzed by Pinkerton and co-workers by performing a QTAIMAC topological analysis of the experimentallyreconstructed electron density [18,19]. Table 2 reports several descriptors evaluated at the various bond critical points of the crystal from the experiments, along with the cor-responding values we obtained from our quantum-mechanical calculations both on the actual periodic extended lattice (Cry) and on the extracted molecular fragment [UO 2 Cl 4 ] 2− (Mol). Calculations are performed at the experimental geometry to ensure a more direct comparison. Overall, the two chemical bonds involving uranium (U-O and U-Cl) can both be classified [55] as of "incipient covalent" character (having 1 < |V|/G < 2, a negative total energy density H < 0, and positive Laplacian ∇ 2 ρ > 0) but show very different chemical bonding features. In particular, the much higher covalent character of the U-O bond is clearly seen in terms of (i) the much larger absolute value of the bond degree |H|/ρ of 0.9 a.u. compared to 0.2 a.u. of U-Cl and (ii) the significantly larger value of the |V|/G ratio of about 1.8 compared to 1.3 of U-Cl. All chemical interactions involving Cs are very weak and of "closed-shell" type with |V|/G < 1, a positive total energy density H > 0, and a positive Laplacian ∇ 2 ρ > 0. The overall agreement between experimental and computed descriptors at the bond critical points is excellent for all the weak interactions (including some subtle aspects such as the Cs-O interaction being characterized by the lowest electron density at the bond CP among all interactions involving Cs), good for the U-Cl interaction and definitely less satisfactory for the U-O interaction. Indeed, as per the U-O bond, while some descriptors (such as the bond degree H/ρ) show a rather good agreement, others (such as the electron density and its Laplacian at the bond CP) show large discrepancies. This was previously noted [38], tentatively attributed to missing environmental effects, and finally shown not to depend on crystal field effects [39]. This latter conclusion is further corroborated by our present results in Table 2 where differences between the molecular fragment [UO 2 Cl 4 ] 2− and the actual extended solid are seen to be very small if not negligible on both U-O and U-Cl interactions. Table 2. Descriptors of chemical bonding of cesium uranyl chloride from the QTAIMAC: bond length l, distance between first atom and bond critical point d CP , value of several local quantities at the bond critical point such as the electron density ρ, the Laplacian of the density ∇ 2 ρ, the ratio between the potential energy density and kinetic energy density |V|/G, and the bond degree H/ρ (i.e., ratio between total energy density and electron density). Calculations are performed at the experimental geometry. Experimental data are taken from Ref. [19]. Calculations are performed both on the actual periodic extended lattice (Cry) and on the extracted molecular fragment [UO 2 Cl 4 ] 2− (Mol). In order to analyze further the origin of the discrepancy in the description of the U-O interaction between theory and experiments, we present deformation density (DD) maps in Figure 2. Deformation density (relative to a neutral atomic reference), ∆ρ(r), contour maps of the cesium uranyl chloride crystal around the U atom are shown in three different planes: (left) through the O-U-O axis and the b crystal lattice vector, (center) the equatorial plane of the four Cl atoms, (right) through the O-U-O axis and a pair of Cl atoms. Upper panels report the experimental deformation density from Ref. [19] while bottom panels report the computed one in this study. The same iso-values have been used for the computed and experimental contour maps. In the panels (A 1 -A 3 ), DD contour lines are overlaid over the ∇ρ trajectories of the crystal electron density, so as to show the intersection of the atomic boundaries with the three considered planes. Both quantitative and qualitative differences can be seen in the description of the electron distribution in the bonding region around the U atom. Panels (A 1 ) and (B 1 ) both show the nearly axial symmetry of the U-O interactions. In particular, the DD of present quantum-mechanical calculations in panel (B 1 ) corroborates [11,56,57] the previously suggested "triple bond" nature of the U-O interaction with a sp hybridization of the oxygen and the formation of a σ bond along the U-O axis and, supposedly, two π bonds with a maximum of charge deformation at about 0.71 Å off the axis. The charge accumulation close to the U atom along the U-O axis in the σ bond in panel (B 1 )-suggestive of a covalent character-is instead missing to a large extent in the experimental DD in panel (A 1 ). Other qualitative differences are observed in (i) a higher degree of localized electron build-up in the region of π bonds in the theoretical DD and (ii) a more pronounced charge deformation behind the oxygen atoms in the experimental DD. Large differences are also observed in the equatorial plane of the four Cl atoms in panels (A 2 ) and (B 2 ). In particular, the expected nearly 4-fold symmetry of this plane seems to be lost in the experimental DD of panel (A 2 ) while it is still largely there in the computed DD of panel (B 2 ). The large departure from the expected symmetry was acknowledged in Ref. [19] and tentatively attributed to the different crystal environment at the second and third nearest neighbor level. While present calculations do show some asymmetry in the equatorial plane (for instance, inspection of the DD along the two green lines in panel (B 2 ) reveals that those two directions are not exactly symmetry equivalent), they predict it to be very subtle while preserving the overall symmetric distribution of the density around the U atom in this plane. However, both theory and experiments describe a charge depletion close to U and a charge accumulation close to Cl along the U-Cl bonds, indicative of a higher ionic character of this interaction relative to the U-O one. Panels  The analysis of the distribution of the electron density ρ(r) is often complemented by the analysis of its Laplacian ∇ 2 ρ(r) or, equivalently, of L(r) = −∇ 2 ρ(r), which provides additional information on the spatial distribution of the electrons and in particular on the asphericity of (bonded) atoms [15]. In particular, some critical points of the Laplacian correspond to charge concentrations and depletions in the core and valence shells. Valence shell charge concentrations (VSCCs) are particularly relevant to the rationalization of chemical bonding, and can be identified in terms of critical points of L(r) of (3, −3) type (i.e., maxima). Figure 3B shows the position in space of the VSCCs critical points (yellow spheres) as obtained from our quantum-mechanical calculations on the [UO 2 Cl 4 ] 2− molecular fragment. The number and position in space of the VSCCs around the U atom are very different with respect to those of the [UF 6 ] − molecule that we studied previously [37]. In that case, 14 VSCCs were found around U while here only 6, all of which lying along the axes of the bonds with Cl and O, thus indicating that charge concentration around U only occurs in the direction of the bonds. The radial distance from U of these 6 critical points is of 0.30 Å and corresponds to the maxima of the VSCC of the principal quantum number n = 6. Along the direction of the U-Cl bond, a critical point is found at 0.64 Å from Cl, which corresponds to the chlorine atom VSCC of n = 3. No other (3,-3) L(r) critical points are found behind the Cl atoms in the molecular fragment, thus indicating that their otherwise spherical density distribution is only altered along the Cl-U internuclear axis. Such an almost spherical electron distribution is broken in the crystal because of the weak interactions with the four Cs atoms surrounding each Cl, with the emergence of further critical points (see below). Finally, along the U-O axis, two critical points corresponding to the O valence charge concentrations are observed: one at 0.35 Å from O toward U and one at 0.34 Å behind O, which suggests a sp-like hybridization of O along the U-O direction. Such a feature could not be clearly seen from our computed DD maps in Figure 2 while it can be further corroborated by inspection of Figure 3C 1 ,C 2 , where maps of the electron localization function (ELF) [58,59] in the region of the O-U-O bonds are shown. Panel (C 1 ) reports the computed ELF from a superposition of non-interacting atomic densities while panel (C 2 ) the ELF from the actual computed density of the interacting system. A η(r) value equal to 0.5 indicates electron localization similar to that of the uniform electron gas, while η(r) values greater or smaller than 0.5 (red and blue iso-contours in the figure) denote higher or lower electron localization compared to the uniform electron gas. The comparison of the two maps clearly shows the higher localization of the electrons both in the σ and π region of the U≡O triple bond and behind the O atoms.

O-U Cl-U O-Cs Cl-Cs Cl-Cs Cl-Cs Cl-Cs
Panels (A 1 ) and (A 2 ) of Figure 3 show the evolution of L(r) (upper panel), ∆L(r) (middle panel) and ∆ρ(r) (bottom panel) along the U-Cl and U-O axes, respectively, for a more quantitative analysis of the features of chemical bonding discussed above. Analogous panels, with a highly magnified interval range on the y-axis and reduced interval range on the x-axis, are reported in Figure 4 to show details of the same functions in the U atom VSCC region. These panels clearly show the difference between the highly polar U-Cl interaction and the definitely more electron shared U-O bond. The L(r) profile in the VSCC and VSCD (valence shell charge depletion) regions of U indicate that the former region is much larger and the latter much smaller when U is bonded to O rather than to Cl, and so are the magnitudes of the corresponding maxima and minima along the profile. In fact, the value of L(r) at the (3,-3) CP in the U VSCC is 177.2 and 210.8 e/Å 5 for U-Cl and U-O bonding interactions, respectively, and the values of ρ(r) at this CP follow a similar increasing trend being 6.0 and 6.5 a.u., respectively. The different nature of the two interactions is confirmed by their ∆L(r) and ∆ρ(r) profiles in the region around the bond critical point (green vertical line in the profiles in Figure 3, with U-O showing a definitely larger charge concentration (larger ∆L(r)) and charge build-up (larger ∆ρ(r)) relative to the U-Cl case. Moreover, such distinct behaviors in the bonding region are observed in a much wider r interval for the U-O rather than for the U-Cl interaction. Finally, it is interesting to note that the two (3,-3) CPs in the VSCC of O are of similar nature (compatible with sp-like hybridization) but not perfectly equivalent. The one facing the U atom and directly involved in the U-O interaction is slightly more displaced from the O nucleus (0.35 Å versus 0.34 Å) and has a larger ρ(r) value (0.877 versus 0.862 a.u.) as a clear sign of (partially) covalent bonding. Moreover, a smaller L(r) value is observed at this (3,-3) CP (3.8 versus 3.9 e/Å 5 ) since the associated VSCC region has expanded toward the facing U atom and become larger in size for promoting the (partially) covalent bonding interaction.  As anticipated above, the electron distribution around the Cl atoms in the crystal gets perturbed by the weak interactions with the four inequivalent Cs atoms surrounding it. This breaking of an otherwise almost spherical distribution is well captured by the appearance of two extra CPs of (3,+3) type of the Laplacian of the electron density, as shown in Figure 5. These two CPs, corresponding to regions of valence charge concentration, are found along directions in between neighboring Cs atoms and not along Cl-Cs axes, denoting the presence of an ionic rather than a Cl-Cs dative interaction.

Conclusions and Perspectives
The TOPOND module (previously limited to work in terms of s-, pand d-type basis functions only) of the CRYSTAL program has been generalized to f -and g-type basis functions, within the linear combination of atomic orbitals (LCAO) approach. The main formal and computational aspects behind this generalization have been presented. The new algorithms now allow for an effective characterization of chemical bonding of materials containing lanthanides and actinides through a topological analysis of the electron density and its Laplacian.
The program has been applied to the analysis of chemical bonding of the cesium uranyl chloride, Cs 2 UO 2 Cl 4 , crystal. Discrepancies between the present quantum-mechanical description of the electron density distribution and that previously reconstructed by experimental X-ray diffraction are illustrated and the implications on the description of chemical bonding discussed.
As a next step, a collaborative study is planned with the researchers who ran the X-ray diffraction experiments and re-constructed the corresponding experimental electron density about 10 years ago to further analyze the observed discrepancies. Indeed, their multipolar model refinement strategy has significantly improved in recent years, which may lead to a closer agreement with present quantum-mechanical predictions.