Electronic structure theory to decipher the chemical bonding in actinide systems

The chemical bonding in actinide compounds is usually analysed by inspecting the shape and the occupation of the orbitals or by calculating bond orders which are based on orbital overlap and occupation numbers. However, this may not give a definite answer because the choice of the partitioning method may strongly influence the result possibly leading to qualitatively different answers. In this review, we summarized the state-of-the-art of methods dedicated to the theoretical characterisation of bonding including charge, orbital, quantum chemical topology and energy decomposition analyses. This review is not exhaustive but aims to highlight some of the ways opened up by recent methodological developments. Various examples have been chosen to illustrate this progress.


Introduction
Historically, there is an industrial motivation for a better understanding of bonding in f-element complexes [1] due to the aim to separate trivalent actinides from trivalent lanthanides in advanced nuclear fuel cycles (see J. Veliscek-Carolan [2] for the experimental methods currently used to extract radioactive actinide elements from solutions of spent nuclear fuel). The areas of study have expanded in recent years with applications in the fields of nuclear toxicology [3][4][5][6], cancer therapy (e.g., the radioisotope 225 Ac of actinium is an alpha emitter useful for targeted radiotherapy [7]), analytical chemistry [8], and in basic research in structural chemistry and reaction mechanisms [9]. All these applications have in common the need to develop pre-organised chelating agents. The problem is complicated by our relatively limited understanding (in comparison to transition metals) of the coordination behaviour of the actinides. In this context, theoretical chemistry is useful to bring additional information to experiments and even by providing in certain cases experimentally inaccessible data. There has been an ongoing debate regarding the nature of the interaction between an actinide element and an organic molecule (ligand) and more specifically on the extent of covalency in actinide-ligand bonds. Most specialists in the field agree that covalency contributes to early actinide-ligand bonding in compounds where the formal oxidation state of the metal is +3 or higher. There is also general consensus that actinide-ligand interactions are predominantly ionic for late actinide compounds in which the formal oxidation state of the metal is +3 [1].
Understanding bonding and interaction in actinide compounds is not straightforward because of their complex electronic structure, the significant electron correlation effects, and the large relativistic effect contributions. The f-block elements are probably the most challenging group within the periodic table for electronic structure theory [10]. Due to this complexity, all the methods used to analyse the nature of chemical bond in organic, inorganic or organometallic compounds are rarely applied to actinides. Common approaches to studying bonding in actinide compounds are based on the shape and the occupation of the orbitals or by calculating bond orders which are based on orbital overlap and occupation numbers. However, this may not give a definite answer because the choice of the partitioning method may strongly influence the result possibly leading to qualitatively different answers. An accurate description of bonding and interactions in actinides compounds cannot be obtained with only a single method. Various complementary tools are mandatory, namely molecular orbital, population, electron density and energy analyses.
In the following sections, the discussion is focused on the tools that theoretical chemists have in hands to investigate the nature of chemical bond and reviews several applications 2 in actinide chemistry. Some of these methods have already been described in the reference [11] and will only be summarized. A number of promising computational approaches are currently under development and have been already successfully applied to non-f-element systems. However, these methods can deepen our understanding of the chemical bond. They will be briefly noted while evaluating their limitations in the context of an application to actinide species. Relativistic effects strongly influence the chemical and physical properties of heavy elements and their compounds. Here, this aspect is outlined in section 5. This review is not exhaustive but aims to highlight some of the ways opened up by recent methodological developments. The examples are voluntarily few in number, but they have been chosen to illustrate this progress.

Population analysis
Partial atomic charges and bond order analyses have been used since the early years of computational chemistry for analysing chemical bonding. They can be used to understand charge transfer and charge flow during chemical processes. Available population analysis or partitioning schemes are based on (a) the basis functions, (b) the electrostatic potential, (c) the wave function or electron density [12].

Basis functions
To this set belong the Mulliken and the Löwdin population analyses [13][14][15][16][17][18], Löwdin partitioning [19,20] and natural population analysis (NPA) [21,22]. In the Mulliken population analysis, the electrons are distributed into the atomic basis functions. The partial charges assigned to atoms vary significantly for a given system when basis sets of different sizes are used. Consequently, results of computations using different basis sets cannot be directly compared. The Mulliken scheme does not require orthogonal basis functions. As a consequence, orbital populations may have negative values, which is physically meaningless [23], or values greater than two in violation of the Pauli exclusion principle. The Löwdin method tries to overcome the limitations of the Mulliken scheme by correcting the instability with increasing the basis set size using orthogonal basis functions. Some examples proposed by Jensen [24] clearly evidence that the Mulliken and Löwdin methods do not converge as the size of the basis set is increased and the values in general behave unpredictable or leads to absurd behaviours especially if diffuse functions are used. Despite these well-known limitations, Mulliken and Löwdin methods were widely successfully used for organic molecules with minimal or small basis sets. For transition metal, lanthanide or actinide complexes, larger basis sets are needed and these analyses lead frequently to quantities that are not chemically meaningful [25].
The LoProp (local properties of quantum chemical systems) approach [12] is designed to avoid this basis set sensitivity and to provide localised quantities like charges, multipole moments, and polarisabilities. The method requires a subdivision of the atomic orbitals into occupied and virtual basis functions for each atom in the molecular system. Several examples can be found on organic molecules [12], metal complexes [26] or bonds containing metals [27]. This method is quite promising for actinide compounds especially since it can be used at the CASSCF level of theory [28].

Electrostatic potential
To this set belong partial atomic charges derived from the electrostatic potential that depends directly on the electron density and therefore on the wave function (ρ(r) = where N is the number of electrons). They are extensively used for force field based classical molecular dynamics and Monte Carlo simulations. There are many possible methods that could be used to determine reasonable atomic charges. It is critical to model the electrostatic interactions at an acceptable level of accuracy since the resulting charges on interior atoms of a molecular system can be unstable and may have unphysical values [29][30][31]. This approach will not be discussed here. Informations can be found in reference publications, reviews and references therein [32][33][34][35][36].

Wave function or electron density
Beyond these relatively simple approximations reported above, much more stable strategies which are more relevant to the study of actinide compounds have been proposed in recent years based on the electron density. They are fully consistent with the philosophy of DFT. They provide a clear partitioning of the electron density and they are less basis set dependent.
The Voronoi deformation density (VDD) [37] and the Hirshfeld [38] population analyses determine atomic charges by integration of the wave function over a region of space. They are quite similar in the sense that charges are obtained from density differences as seen in equation (1a): The VDD is based on a partitioning of space into non-overlapping atomic volumes (1b), the Voronoi cell of each atom, defined as the volume which is closest to the nucleus of that atom. Guerra et al. [39] showed that an analysis of the deformation density in Voronoi cells -bonded minus free atoms -provides chemically intuitive partial atomic charges.The VDD of an atom monitors the flow of charge into or out of the atomic Voronoi cell as a result of the chemical interactions between atoms. Hirshfeld charges are derived by partitioning the electron density into the contribution of each atom in proportion to the electron density of free neutral atoms [38]. The idea behind this approach is that the electron density at every point is distributed between all atoms based on a weighting function for each atom (1c). This allows for atomic electron densities to overlap. The Voronoi deformation density and Hirshfeld atomic charge of an atom A can both be written as the integral over all space of the deformation density times a weight function w M ethod A (Method = VDD, HPA) [40]. Integration of the atomic deformation densities defines net atomic charges.

4
The VDD and Hirshfeld charges are not free from issues. The major drawback is that the atomic charges depend upon the (arbitrary) choice of atomic reference densities (the neutral isolated atom), that may not be suitable for charged systems. Guerra et al. [39] have pointed out this aspect and have investigated various alternatives to the VDD framework. In the same way, in the original Hirshfeld approach, the weighting functions for each atom were based on the electron density of the isolated atoms in the gas phase, leading to atomic electron densities that are dependent on the choice of the isolated atoms electron configuration [41]. This dependence is mostly removed by the iterative Hirshfeld method (Hirshfeld-I) introduced by Bultinck et al. [42]. However calculations of gas-phase electron densities are still required. To solve the issue several approaches were suggested these last years (see for example references [41,43]). These methods are relatively insensitive to basis set. Unfortunately, they are rarely implemented in quantum chemistry programs. For example, Hirshfeld charges can be obtained by post-processing the wave function with the program HPAM written by Elking et al. [44] or Multiwfn written by Lu [45].
More recently, Marenich et al. proposed a novel approach to deriving partial atomic charges from population analysis, called Charge Model 5 (CM5) [46]. CM5 charges are derived by mapping from the Hirshfeld population analysis (HPA) through the following equations [46]: where the indexes i and j run over all atoms in the molecule, Z i and Z j are the corresponding atomic numbers, R Z is the atomic covalent radius based on the data from refs [47][48][49], and q i is the partial atomic charge from Hirshfeld population analysis. The quantities α and T ij (where T ij = −T ji ) are model parameters to be determined. The quantity B ij defined by (2b) was called the Pauling bond order in their work by analogy to a similar equation proposed by Pauling (see eq. 1 in ref [50]).
In the Quantum Theory of Atoms in Molecules (QTAIM), Bader [51] proposed to partition the space of a molecular system into non-overlapping mononuclear regions, the atomic basins Ω. They are separated by interatomic surfaces S(r) that satisfy the boundary condition of zero flux of the gradient vector field of the electron charge density (3): ∇ρ (r) · n (r) = 0 for all r belonging to the surface S(r) where r is the position vector and n(r) the unit vector normal to the surface S(r). An atom in a molecule is defined as the union of a nucleus and its associated basin [52]. By integrating the electronic density (4) within the atomic basin Ω A the total charge on an atom A can be calculated [44] as: One of the most widely used wave function based approach for the investigation of chemical bonding (including in f-element chemistry) is Weinholds Natural Bonding Orbital (NBO) analysis [21,22]. NBO is part of a more comprehensive framework which comprises a sequence of transformations from the input basis set (non-orthogonal atomic orbitals (AOs)) to various localized orthonormal basis sets (natural atomic orbitals (NAOs), hybrid orbitals (NHOs), bond orbitals (NBOs), and localised molecular orbitals (NLMOs)) [21,[53][54][55]: In the open-shell case, the analysis is performed in terms of "different NBOs for different spins", based on distinct density matrices for α and β spin [55]. Because the localisation is based on a maximum occupation criterion, the NBO analysis provides a description of the molecular wave function close to the chemical concept of Lewis structures. As a consequence, these methods only work well for systems that can be described in this formalism. However, for simple electron configuration, they usually give a good qualitative picture of the bonding and are quite robust when increasing the basis set. The multipole derived charges (MDC) analysis [56] uses the atomic multipoles (obtained from the fitted density), and reconstructs these multipoles exactly by distributing charges over all atoms. This is achieved by using Lagrange multipliers and a weight function to keep the multipoles local. The recommended level is to reconstruct up to quadrupole, i.e., MDC-q charges [56].
As an illustrative example, we consider the gas-phase organoactinyl complexes possessing discrete An -C bonds (An = U, Np, Pu) recently synthesized by Dau et al. [57] in a quadrupole ion trap by endothermic decarboxylation of [AnO 2 (O 2 C -R) 3 ]anion complexes. A formally AnO 2 2+ actinyl core is coordinated by three carboxylate ligands, with R= CH 3 (methyl), CH 3 CC (1-propynyl), C 6 H 5 (phenyl), C 6 F 5 (pentafluorophenyl). The chemical properties of these bonds were analysed by QTAIM [57]. Here, we calculate the partial atomic charges for the decarboxylation products, [(R)AnO 2 (O 2 C -R) 2 ]with R = CH 3 and An=U. The results are given in Table 1 for several population analysis methods and with increasing basis set size. The partial charges assigned to atoms using Mulliken population analysis vary significantly for the same system when different basis set sizes are used. This instability with increasing basis set size might be unacceptable to calculate the electronic structure and properties of actinide compounds. For all other methods, the variation with basis set size is quite small or negligible. The uranium charge (Table 1) is clearly underestimated by Hirshfeld, VDD, CM5 and NPA population analysis. As comparison, q(U)=2. 78 and q(O yl =-0.39 were obtained from multireference calculations and a multipole analysis on uranyl -water complex by Hagberg et al. [58]. Although the populations analysis schemes described above are useful for deriving chemical insights, one should remember that neither calculated atomic charges nor bond orders are observables from a quantum mechanical point of view.

Topological analysis of the electron density
The electron density is a physical observable obtained experimentally by high resolution X-ray diffraction and subsequent multipole refinement. The comparison of key topologi-cal indices of quantum chemically computed and experimentally determined charge density distributions has become an important field. In many cases the qualitative features of experimentally and theoretically derived densities agree well but the question whether a topological analysis provides measurable quantities is still open [59]. The theoretical approach should rather be seen as providing an intelligible link between experiment and theory. Bond descriptions based on the topological analysis of various molecular fields, such as the electron density [60,61], the Laplacian of the electron density [51,62], the electron localisation function [63,64] or a measure of electron localisability [65] were introduced during these last years. Once again, they have rarely been applied to the study of actinide compounds. This section contains a summary of the main concepts of recently developed methods. Some detailed reviews on the analysis of the molecular fields including the description of electron delocalisation can be found in the literature [51,[66][67][68][69][70][71] or in this issue with the contribution of Lepetit and Silvi.
Within QTAIM, the electron density ρ(r) is the central property and provides physical basis for many chemical concepts. The electron density is a scalar field and its topological properties are determined by the analysis of its associated gradient vector field. The first important feature is the location in the space of the critical points (CPs). A CP in the electron density is a point were the gradient of the density vanish, i.e. ∇ρ = 0. In a 3D scalar field such as ρ(r), there are three types of CPs: a maximum, a minimum or a saddle point [70]. In three dimensions there are two types of saddle points. Critical points can be characterized by the eigenvalues λ i (i = 1, 2, 3) of the Hessian matrix H xyz of ρ, evaluated at the CPs, i.e.: The critical points are classified according to their rank, R, the number of non-zero eigenvalues, and the signature S, that is the algebraic sum of the sign (sgn) of the eigenvalues [66,70]: For example, one type of saddle point has two strictly negative (i.e. non-zero) eigenvalues and one strictly positive one. Consequently, its rank is three, and its signature is (-1)+(-1)+1=-1 [66,70]. We designate this point as a (3, -1) CP, where the first index refers to the rank and the second to the signature. This particular type of CP is called a bond critical point (BCP) because it indicates the existence of a bond between two nuclei of a molecule in an equilibrium geometry [70]. A (3, -3) CP is a local maximum, a (3, +1) is a minimum in two directions and a maximum in the other direction, and finally, a (3, +3) is a local minimum. Each type of critical point described above is identified with an element of chemical structure: (3,-3) nuclear critical point (NCP), (3,+1)ring critical point (RCP), and (3,+3) cage critical point (CCP). The bond critical points are linked to the nuclei via the atomic interaction line. This line consists of a pair of gradient paths, each of which originates at the bond CP and terminates at a nucleus. ln a bound system, an atomic interaction line is called a bond path [72]. The set of all bond paths occurring in a molecule (or molecular complex) is called a molecular graph [70,73]. Figure 1 shows the molecular graph of the bent U(η 8 - [74] discussed in more detail later in this section. A molecular graph and the characteristics of the density at the bond critical points provide a concise summary of the bonding within a molecule. Following Bader [51], there are basically two kinds of bonding interactions: sharedelectron interaction and closed-shell (or more generally unshared-electrons) interaction. Covalent and dative interactions are classified as shared-interactions whereas ionic, electrostatic, hydrogen and van der Waals bonds are classified as closed-shell. Metal-metal bonds lie in an intermediate situation between typical ionic and typical covalent bonds [52,75,76]. For a classification of bonds, the sign of the Laplacian of the density ∇ 2 (r) (the trace of the Hessian of the density) at the BCP is checked out. A negative value indicates a local concentration of the electron density, whereas charge depletion is characterized by a positive Laplacian. This point has been amply reviewed in the literature and can be summarized in Table 2 [77]. However, for weak or highly polar bonds, the analysis of the density and its Laplacian (close to zero) does not suffice to characterise the bonding. Additional information can be obtained from the eigenvalues (e.g., η = |λ 1 | λ 3 ) or the total electronic energy density H(r) (6) as the sum of the kinetic G(r) and the potential energy V (r) [78]: AIM also allows to obtain localisation λ(A) and delocalisation indices δ(A, B) for atoms and pairs of atoms in molecules [79] in order to gain quantitative information on the extent to which electrons are localized to a given atom or shared with others. The AIM theory recovers the chemical picture of the molecule made of atoms but the partition does not explicitly reveal a sub-structure corresponding to the cores and the valence shell of the atoms. The electron localisation function (ELF) approach [63,64] belongs to the same type of methodology. It attempts to overcome the conceptual limits of the topological analysis of the sole electron density [80]. A different partition of space is performed using ELF. As already specified, in the AIM theory, the basins are defined as a region of space bounded by a zero-flux surface in the gradient vectors of the one-electron density, ρ (r), or by infinity. As ELF is a scalar function, the analysis of its gradient field can be carried out to locate its attractors (local maxima) and the corresponding basins [81]. For an N-electron single determinantal closed-shell wave function built from Hartree-Fock (HF) or Kohn-Sham (KS) orbitals, ψ i , the ELF function can be written using the following equation (7): ELF represents a continuous and differentiable scalar field η(r) as does the electron density ρ(r). Therefore, topological analysis of η(r) is technically similar to the one put forward by Bader and coworkers for ρ(r). The topological analysis of η(r) is performed via its associated gradient vector field ∇η(r). This field is characterized by critical points. They represent local maxima, minima and saddle points of η(r).
Another interesting approach was proposed by Kohout [65,82], i.e., the electron localizability indicator (ELI). It is a functional that describes the correlation of electronic motion. In fact ELI designates a family of indicators. Among them, one of the most pertinent for the study of the bonding in molecular systems is the ELI-D, symbol Υ D (8). ELI-D can be seen as being proportional to the charge that is needed to form an electron pair. Based on the pair density, it is not explicitly accessible from a Density Functional Theory (DFT) calculation [83]. Consequently, the ELI-D is strictly justified within wave function-based theories, e.g., Hartree-Fock (HF) and post-HF methods such as the complete active space (CAS) and configuration interaction (CI) methods [84][85][86]. However, at least formally, the functional can be derived from the Kohn-Sham orbitals within the framework of DFT yielding reliable results [87]. A detailed description of the formulation is available in references [82,[87][88][89]. In brief, Υ D (r) can be expressed as the product of the electron density with the pair-volume function (8): At HF level for a spin coordinate σ: For closed-shell systems at HF level, the ELI formula for the same-spin electron pairs resembles the kernel of the ELF one of Becke and Edgecombe [63], but without any reference to the homogeneous electron gas. This is one of the key feature of the ELI functional. The ELI is a real function with positive values and is not limited to an interval [0,1] as in the case of the ELF. As a scalar field, the topology of ELI-D can be analysed in term of attractors, basins of attractors and properties at critical points as the density Laplacian or the kinetic energy density. Localization functions based methods provide a meaningful representation of atomic shells, lone pairs, and covalent bonds, but usually, they do not reveal non-covalent interactions. Johnson and co-workers [90] introduced a method for studying non-covalent interactions (NCIs) on the basis of the electron density ρ (r), the reduced gradient (RDG) of the density (s = |∇ρ| /(2(3π 2 ) 1/3 ρ 4/3 )) and the Laplacian of the density ∇ 2 ρ (r). The approach allows one to identify the interactions in real-space, thus enabling a graphical visualisation of the regions where non-covalent interactions occur. NCI visualises both intra-and intermolecular weak interactions through RDG isosurfaces at low electron density values, however, it does not reveal covalent bonds. Combined NCI/ELF [91] or NCI/ELI-D [92] analysis provide a way to visualize simultaneously non-covalent and covalent interactions despite a different value range. A more elegant way is the Density Overlap Regions Indicator analysis proposed by de Silva and Corminboeuf [93]. It should be seen as a modification of the Single-Exponential Decay Detector, a density-based bonding descriptor [94,95]. This is a scalar field, which reveals both covalent and non-covalent interactions in the same value range. It is derived from the electron density and its derivatives (9): θ (r) becomes infinite in bonding regions. A transformation (10) is applied for mapping to the [0,1] range: γ (r) is close to 1 in bonding regions and close to 0 at nuclei and far from the molecule. In fact, γ does not probe electron localisation but rather geometrical features of the electron density in terms of deviations from single-exponentiality. The idea is based on the fact that atomic densities are approximately piecewise exponential [96][97][98], whereas in the interaction regions they are not, due to the overlap of two or more atomic densities [93]. γ can be identified as a detector of strong relative density overlap regions and is named Density Overlap Regions Indicator (DORI). As DORI identify regions of overlapping densities, it does not quantify the strength of the interaction and does not distinguish between attraction and repulsion.
In line with the NCI index, this limitation is resolved by combining the analysis of DORI with that of the electron density Laplacian (∇ 2 ρ(r)). In particular, the second eigenvalue λ 2 < 0 identifies bonding regions, while λ 2 > 0 indicates non-bonding interactions. Along with its sign, the magnitude of the interaction is estimated from the values of the density itself. de Silva and Corminboeuf [93] use sgn (λ 2 ) ρ (r) as a complementary scalar field. DORI simultaneously uncovers regions of covalent as well as strong and weak non-covalent interactions. Non-interacting lone pairs are, however, not visible within DORI. We should note, that the introduced scalar field depends only on the electron density; thus, it is welldefined at any level of theory. In particular, the DORI analysis is easily applicable to densities obtained from post-Hartree-Fock methods. This is of particular importance for actinide chemistry.
As an illustrative example, we consider the bent U( orbital and ELF analysis [74]. Despite the broken symmetry (rigorous D 8h symmetry for U(η 8 -C 8 H 8 ) 2 ), the gain in electrostatic interaction and a significant U -CNorbital interaction are sufficient to stabilise the bent CNcomplex with respect to the bare uranocene. As previously found for the U(η 8 -C 8 H 8 ) 2 compound [100,101], 6d, and to a less extend 5f, uranium orbitals have a significant participation in the interaction both with the aromatic rings and the cyanide ligand. The results support the hypothesis of a non-negligible covalent character in the U···CNinteraction. The QTAIM molecular graph ( Figure 1) shows a bond path and a BCP between the U···C atoms. At this BCP the Laplacian of the density is positive and the total energy density is negative indicating a covalent nature of the U···CNinteraction [72,102,103]. The interaction between the U 4+ ion and the CNligand is also characterized by the formation of a ELF V(U,C) basin with a population of about one electron. However, the visualisation of the ELF function is not sufficient to provide a clear understanding of the nature of the bonding in the whole system. A more comprehensive picture can be drawn using regions of density overlap. Figure 2 shows the DORI=0.9 isosurface for U(η 8 -C 8 H 8 ) 2 (CN)complex with color-coded values of sgn(λ 2 )ρ(r). The bonding domains are coded in red color, the repulsive atom-atom contacts in blue color and the noncovalent interaction in green color. The covalent interaction between U and CNis clearly identified in Figure 2 as the "red button". This covalent interaction allows the formation of a stable bent complex between U(η 8 -C 8 H 8 ) 2 and CN -.
These results highlight the interest of the use of several topological analyses for understanding bonding in complex systems. However, limitations exist in the case of actinide compounds. They will be analysed in section 5.

Energy decomposition analysis
For a more complete view of the chemical bond formation, insights could be provided by an energy decomposition analysis (EDA) into different terms, which can be interpreted in a physically meaningful way (e.g., electrostatics, exchange-repulsion, polarisation (or induction), and charge transfer). During the past several decades, several EDA fragment-based approaches have been developed and used in a wide range of applications in chemistry. All of them have evolved from the early EDA of Kitaura and Morokuma [104][105][106]. These developments have taken place aiming to overcome limitations of the original schemes and provide more chemical significance to the energy components, which are not uniquely defined. We can cite, for example, CSOV (Constrained Space Orbital Variations) [107], RVS (Reduced Variational Space Self-Consistent-Field) [108], SAPT (Symmetry-Adapted Perturbation Theory) [109], NEDA (Natural Energy Decomposition Analysis) [110][111][112], LMOEDA (Localized Molecular Orbital Energy Decomposition Analysis) [113], ALMO-EDA (Absolutely Localized Molecular Orbital [114]), FMO (Fragment Molecular Orbital [115][116][117]) and a Morokuma-type EDA developed by Ziegler and Rauk [118][119][120]. Certainly, this list is not exhaustive, but recalls some popular currently used EDA schemes. A more comprehensive description of the various approaches can be found, e.g., in references [121,122]. As pointed out by Phipps et al. [122] the application of a scheme to a particular system may be limited by the level of theory at which the scheme is implemented. Schemes implemented at the higher levels of theory often include dispersion and other correlation components which are not available at the HF level. Due to the complexity of the electronic structure of the actinide complexes, electron correlation and relativity effects have to be taken into account. Therefore, an usable method is not always accessible, each method having well defined application domains. As an illustration, the localised molecular orbital LMOEDA scheme of Su and Li [113] is implemented for the analysis of both covalent bonds and intermolecular interactions on the basis of single-determinant HF (restricted closed shell HF, restricted open shell HF, and unrestricted open shell HF) wave functions and their DFT equivalents. For HF methods, the total interaction energy is decomposed into electrostatic, exchange, repulsion, and polarisation terms. Dispersion energy is obtained from second-order Møller-Plesset perturbation theory and coupled cluster methods up to CCSD(T). Recently, has been proposed an energy decomposition analysis of open and closed shell f-elements mono aqua complexes with the CSOV scheme using a MCSCF wave function [123]. The procedure is derived from the Bagus and Bauschlicher one [124] and is currently valid for open shell high-spin complexes. In the various methods commonly used, a charge transfer (CT) term may be present or not. When the CT contribution is not separated it is rather included within the induction term. In this case, one usually groups all of the terms assignable to the bond (polarisation, charge transfer, etc.) into a single orbital interactions term. The physical significance of the CT term is currently an important subject of controversy related to the arbitrary definitions of the interaction energy components. For example, at short intermolecular distances and especially with large basis sets the separation between charge transfer and polarisation becomes increasingly ill-defined. At larger intermolecular distances, CT becomes more easily separable from polarisation [122]. Novel approaches are under development in order to provide more meaningful and robust definitions for the various components of the intermolecular interaction energy including CT. As an example, Rezác and de la Lande [125] propose an approach to calculate the energy associated to intermolecular charge transfer. This component of the interaction energy is calculated as the difference between a standard DFT calculation (with full relaxation of the electron density) and a calculation in which CT between the molecular fragments is prohibited by means of constrained DFT (cDFT) methodology [126]. This approach was originally formulated by Wu et al. [127]. The present implementation supports only GGA and meta-GGA functionals and will be extended to hybrid functionals. Lao and Herbert [128] proposes a composite framework with cDFT-based definition of the CT energy in conjunction with a SAPT-based definition of the remaining energy components. SAPT/cDFT when combined with the new implementation of open-shell SAPT [129], can be applied to supramolecular complexes involving molecules, ions, and/or radicals [128]. With the same aims, in the framework of orbital-based EDAs, Levine et al. [130] developed an extension of the ALMO-EDA [114] specifically for the variational energy decomposition analysis of chemical bonds relative to radical fragments. The total interaction energy is broken up into four terms: frozen interactions, spin-coupling, polarisation, and charge-transfer. This method includes a new EDA term, spin-coupling, to describe the covalency of the bond in an energetic way. Currently, this method is only implemented for single bonds corresponding to the CAS(2,2) wave function [130]. Of course, these current developments are far from an application to actinide complexes but they are very promising for the future interpretations.

Some aspects of relativistic quantum chemistry
It was P. Pyykkö who showed that relativistic effects can strongly influence the chemical properties of the heavy elements and their compounds [131,132]. The starting point is quite simple but it has great impact: for the heaviest elements, the inner electrons are approaching the speed of light leading to a mass increase and a following smaller Bohr radius. These effects have three main consequences on the orbitals of many-electron atoms (e.g., see Figure  3) • Contraction and stabilisation of all s and most p-type orbitals (direct effects) • Expansion of d-and f -type orbitals (as a consequence of the stronger screening of the nuclear attraction by the contracted s and p) • Spin-orbit splitting of shells with l > 0 Many examples are discussed in the literature, e.g. in references [133][134][135]. The RTAM searchable database [136] also includes an up-to-date list of relativistic calculations for atoms and molecules. A well-known example is the lanthanide contraction of the Ln -X bond lengths in LnX 3 molecules from Ln = La to Lu. Relativistic effects account for 9-23% of the total contraction depending on the system [137]. One can also observe relativity in everyday life, e.g., the yellow color of gold and the liquidity of mercury [135]. The lead-acid battery, which is commonly used in cars, strongly relies on the effects of relativity [138]. As put by Pyykkö "relativistic effects in chemistry are more common than you thought" [135]. In the past two decades, large progresses have been made in the development of both relativistic methodology and computer codes [139,140]. But fundamental questions still remain open in relativistic electronic structure theory. In 1928, Dirac proposed a relativistic formulation of the quantum mechanics [141]. In addition, Dirac's theory is the basis for modern quantum electrodynamics, one of the most accurate quantum theories to date. The relativistic equation corresponding to the Hartree-Fock equation is the Dirac-Fock equation, which in its time-independent form can be written as (11): where α and β are the Dirac 4 x 4 matrices. α is written in term of the three Pauli 2 x 2 spin matrices σ and β in term of a 2 x 2 identity matrix I:  [24]. The usual notation of the wave function is in the form (12): where Ψ L and Ψ S are the large and small components of the wave function and α and β the spin functions. For electrons, the large component reduces to the solutions of the Schrödinger equation when c → ∞ (the non-relativistic limit), and the small component disappears. The small component of the electronic wave function corresponds to a coupling with the positronic states (negative-energy states) [24]. Moreover, the use of the Dirac Hamiltonian requires considerably larger computational resources compared to the use of the Schrödinger Hamiltonian. In chemical applications, one is usually concerned with the electronic (or positive-energy) states only. Therefore, some reductions of the four-component wave function have been proposed leading to the two-component methods of relativistic quantum chemistry [142], e.g., the Douglas-Kroll-Hess (DKH) transformation [143,144], the zeroth-order regular approximation (ZORA) [145], the infinite-order two-component theory (IOTC) [146,147] and the exact two-component method (X2C). All of them account for spin-orbit interaction. These developments increased the size of the system that can be treated computationally. Nowadays, more extensive savings may be achieved by introducing the use of effective core potentials. For first-row systems, the savings are marginal but for heavier elements the computational cost may be significantly reduced. Compounds involving heavy elements have a large number of core electrons. These are often considered as unimportant in a chemical sense, but it is necessary to use a large number of basis functions to expand the corresponding orbitals, otherwise the valence orbitals will not be properly described (due to a poor description of the electron-electron repulsion) [24]. These problems and the role of relativity could be addressed simultaneously by modelling the core electrons using a suitable function, and treating only the valence electrons explicitly. This function is called an effective core potential (ECP). In this framework, the relativistic terms are implicitly incorporated into the potential considered by parametrization with respect to suitable relativistic allelectron reference data. A detailed description of the formalism and limitations can be found in references [148][149][150]. The vast majority of electronic structure calculations on actinide species uses this simple one-component ECP method including only scalar relativistic effects. For some recent comprehensive guides to the theoretical background and the computational applications of relativistic quantum chemistry, we quote Schwerdtfeger [139,140], Hess [151], Dyall and Faegri [152], Barysz and Ishikawa [153], Reiher and Wolf [154], or in reviews and reports from Saue [155], Pyykkö [156] or Liu [157].

Bonding analysis
Regarding the analysis of chemical bonding with the previously described methods in section 2 to section 4, we must distinguish between scalar-relativistic approaches and the two-and four component ones.
Everything available in the non-relativistic framework can be used in scalar-relativistic calculations. Population and energy decomposition analyses are currently available in several quantum chemistry package. If there are no issue with these methods, it is important to be very careful when ECPs are used for topological analyses [158][159][160]. The absence of core electrons for an atom results in the absence of a corresponding nuclear maximum of the electron density. Therefore, the qualitative change in the topology of the electron density in the neighbourhood of the atom, makes meaningful topological analyses difficult or impossible [160]. Several studies tried to recover this issue (see i.e. references [160][161][162][163]). The missing electron density contribution can be restored using an atomic core density obtained from an all-electron calculation for the free atom. When a small-core or medium-core ECP is used, representing the corresponding core electron density using a tightly localized electron density function may be sufficient to obtain a correct electron density topology and perform QTAIM analyses with qualitative meaningful results [160]. This is not true when large-core ECPs are used. In general, a quantitative analysis cannot be obtained without using an all-electron basis set.
As previously seen, scalar-relativistic effects can be introduced through DKH or ZORA Hamiltonians for example, and spin-orbit coupling via a two-or four-component approaches. While the population analysis can be performed from one-to four-component calculations, EDA and even more topology analyses are either scarce or non-existent for two-or fourcomponent calculations. It is currently possible to perform topological analyses with a two-component DKH or ZORA Hamiltonian. To the best of our knowledge, no more than the electron localisation function, tested only for closed shell systems, has been implemented in an all-electron four-component approach (DIRAC program [164]).
The question is: can one find an effect due to spin-orbit coupling so large, that one would see a difference in the topology of the electron density? In which chemical cases? At a first glance, the largest effect due to relativity is expected within the inner shells of the heavy atoms. The critical points, however, are located in the valence region of the atoms so that one might expect the deviations due to relativistic effects to be small at these points. Very few systematic studies of the effects of the relativity on the topology of the electron density exists in the literature for actinides or other heavy elements of the periodic table (e.g. see [165,166]). Amaouch et al. [167] propose a QTAIM and ELF topological analyzes of the bonding in At 2 and UO 2 2+ in the framework of two-component relativistic DFT comparing scalar-relativistic and spin-orbit coupling calculations. In UO 2 2+ , among ELF and QTAIM descriptors, the most affected one by spin-orbit effects is the Laplacian of the density. ∇ 2 ρ(r) is increased by 10% in UO 2 2+ with respect to scalar relativistic calculations, showing a more depleted electron density at the U-O BCPs [166,168]. Anderson et al. [169] perform DFT and QTAIM calculations at both the non-relativistic level and using the scalar relativistic ZORA for H 2 O, X(C 5 H 5 ) 2 (X=Fe, Ru, Os), Au 4 and UF 6 compounds. In the case of UF 6 , the difference in the electron density and its Laplacian at U-F BPCs is ca 2% and 19%, respectively between the SR-ZORA and the non-relativistic frameworks.
The main conclusion of the various studies (including non-actinides heavy elements) is that relativistic effects are not only restricted to the core electron density of heavy elements but are also clearly discernible in the outer core and the valence region within adjacent bonds [165]. In most cases scalar relativistic variants of the Hamiltonian operators as ZORA and DKH are sufficient to obtain a reasonable description of the electron density [166,168]. These effects must be taken into account for a meaningful topological analysis of bonding in actinide compounds. As the Laplacian of the electron density at the bond critical points is a very sensitive quantity to detect changes in a density distribution, one finds significantl larger relative deviations including spin-orbit coupling. However, twocomponent calculations including spin-orbit effects can provide a detailed insight into the topology of the electron density for atoms and molecules. Four-component calculations are most of the time not required.

Selected applications in actinide chemistry
In this section, we shall focus on some recent examples that illustrate applications of analysis techniques presented above. Further examples can be found in reference [170]. From the past two decades, DFT becomes the unchallenged method for describing electronic structure and bonding in medium-and large-size molecular systems. Examples are reported in reviews dedicated to inorganic or organometallic compounds [171][172][173][174]. DFT is not free from issues in metal or f-element chemistry. Kepp recently provided a detailed discussion of specific systematic effects often ignored in mainstream DFT modelling [172]. Their review concerns recent efforts to accurately and consistently describe bonds across the s-, p-, and d-blocks chemistry. Physical effects and ingredients in functionals, their systematic errors, and approaches to deal with them are discussed, in order to identify broadly applicable methods for inorganic chemistry. There are numerous examples in the literature where the bonding analysis is achieved using Mulliken or NBO methods despite the fact that they provide at best a qualitative description and can give somewhat contradictory results in assessing covalency [175].

Beyond the 18-electron rule with a central actinide atom
Molecules may become particularly stable when the electron count at a central atom reaches a magic number. Two well-known examples are Lewis' octets and Langmuir's 18electron principle, corresponding to filling s + p-and s + p + d-like shells. For the history and actual interpretation of the 18-electron principle, see Pyykkö [176]. If, in addition to singly occupied s + p + d shells, f-or f-like shell at the central atom would be occupied, one would obtain systems with the magic number 32 [177]. Obviously, the central atom would have to be an actinide, with a relatively diffuse and energetically available 5f shell. We combine charge, orbital, quantum chemical topology and energy decomposition analysis to propose An@Pb 12 [178], An@Sn 12 [179], An@C 28 [180] and An@Si 20 [181] series as 32electron systems. For a comprehensive review, see [177]. Let us consider the example of An@Pb 12 with An=U, Pu, Am, Cm. The insertion of an actinide atom highly stabilizes the cage owing to a strong Pb-An interaction (I h symmetry, see Figure 4) which originates in an attractive electrostatic effect and in a large orbital interaction (see Table 3). The Pauli repulsion remains nearly the same from U to Cm. The orbital term contributes ∼ 45% to the total attractive interactions and is determined by the 6d orbital energy for the Pu and Am clusters and by the 5f orbital energy for the Cm cluster. Further analysis from the valence molecular orbital (MO) diagram revealed the formation of sixteen MOs between the cage and the central atom (a g , g u , h g , t 1u and t 2u ). The Pu 7s, 7p, 6d and 5f orbitals are involved in the hybridization with the Pb 2− 12 cage orbitals. In the [An@Pb 12 ] n clusters, the ELI-D map exhibits localisation domains between the central atom and the cage (diagram (a) in Figure 5). They result from the interaction between Pb and Pu for example, and are a consequence of the formation of the a g , g u , h g , t 1u and t 2u orbitals. e.g., the diagram (b) in Figure 5 depicts the partial orbital contribution (pELI-D) of the 1h g orbitals in the real space. This accounts for the existence of a strong interaction between the central atom and the cage. The combination of charge, orbital, quantum chemical topology and energy decomposition analyses provides an unambiguous picture of the bonding in the [An@Pb 12 ] n clusters with notable covalency between the central actinide atom and the surrounding cage. These results have enabled us to conclude that Pu@Pb 12 and [Am@Pb 12 ] + compounds are the first examples of centered 32-electron systems.

Assessing covalency in nitrogen ligand-actinide complexes
Actinide recycling by separation and transmutation is considered worldwide as one of the most promising strategies to reduce the inventory of radioactive waste [182]. A liquidliquid extraction and treatment strategy is currently implemented [183]. One of the most challenging task is the partitioning of trivalent actinides, Am 3+ and Cm 3+ from trivalent lanthanides. The two series of ions are both hard Lewis acids, with similar charges and ionic radii. Few years ago, it was discovered that nitrogen ligands are able to exhibit a selectivity between An and Ln [184]. Therefore, the researches were focused on the chemistry of polyazine extractants [185]. The design of extractant ligands for lanthanide/actinide separation exploits the possibility of a significant covalency in the 5f series. But, as stated by Girnt et al. [186] and by Kaltsoyannis [175], "the level of understanding of BTPs selectivity on a molecular level is insufficient to target the design of new, more efficient, and selective partitioning reagents or fine-tune partitioning process conditions. Such advances are presently empirical, on a trial and error basis". Most of the theoretical chemistry works tried to estimate the covalency in the f-element nitrogen bond and to determine if there are differences between the minor actinide and lanthanide compounds that would account for the observed separation factors [175]. These works are mainly based on DFT calculations with "standard" functionals and population analysis tools. The results are often contradictory and failed to improve our knowledge of bonding in such compounds. di Pietro and Kerridge [187] considered four different methods for studying the bonding in BTP ( Figure  6a) and isoamethyrin (Figure 6b) complexes of uranyl. To avoid the ambiguity which can arise from orbital based methods, the author focused solely on properties derived from the electron density. Their analyses involve the use of QTAIM, ELF, RDG for the regions of weak interaction and the density differences upon complexation. They studied the nature of U -N bonding in these complexes and the consequent effects on the highly covalent U -O bond of uranyl. The analyses demonstrated weak, but non-negligible, covalent character in the U -N bonding region of both complexes. As might be expected, the covalent character of the bonds increases as the U -N bond length is shortened. Analyses revealed a strong effect on the uranyl U -O bonds upon complexation, namely a noticeable reduction in electron sharing in the U -O bonding region, with charge instead localising on the oxygen centres. This leads to an increase of ionic character in the U -O bond. This, of course, also corresponds to a reduction in covalency.

The actinocenes An
As mentioned earlier, the electronic structure of actinides compounds can be highly complex mainly due to the strong electron correlation within the valence orbital manifold and the effects of relativity. One of the most typical example is the actinocene series. In the development of modern organometallic and theoretical chemistry, uranocene U(η-C 8 H 8 ) 2 (or U(COT) 2 ) has played a central role in understanding the fundamental bonding, electronic structure and chemical properties of actinide complexes. The possible existence of the uranium complex U(COT) 2 was suggested in 1963 by Fischer [188] who noted that there might be "a possible additional gain in energy due to the participation of the f-orbitals of heavy central atoms". The synthesis in 1968 of this first linear sandwich compound of an f-element, by Streitwieser and Mueller-Westerhoff [189,190] who christened uranocene to highlight its similarity with ferrocene, was a milestone in the history of actinide chemistry [191][192][193]. Since the theoretical prediction [188] and subsequent experimental discovery of uranocene, a series of actinocenes An(COT) 2 (An=Th, Pa, U, Np, Pu, Am) or their anions have been synthesized and characterized. Historically, the first theoretical work on U(COT) 2 including relativistic effects was done by Pyykkö and Lohr in 1981 [100] using relativistically parametrised extended Hückel calculations. Whereas a single-determinant wave function could describe the ground states of thorocene and uranocene [194,195], a multiconfigurational description of the wave function is needed for the later actinides [196,197]. For example, Kerridge and Kaltsoyannis [196] performed all-electron spin-orbit coupling complete active space self-consistent field calculations including dynamic correlation via second order perturbation theory (SOC-CASPT2) for the ground and low-lying excited states of M(COT) 2 (M=Th, U, Pu, Cm). The multiconfigurational character was associated to the occupation of the 5f δ and π 2u orbitals [198]. In all the above mentioned calculations, the complexes were assumed to be planar and of D 8h point group symmetry. The electronic structure was mainly described in term of molecular orbital energy diagrams and Mulliken population analysis. Nevertheless, the analyses methods employed did not provide a fully clear picture of bonding. The orbital based measures of covalency still remain ambiguous. Recently, Kerridge has used a topological analysis based on QTAIM and a CASSCF wave function in the actinocenes series (An=Th-Cm) [198]. Combining analyses of ρ and its Laplacian at the An -C bond critical points (see the molecular graph in Figure 7), of the localisation and delocalisation indices λ(An) and δ(An, C) and the electron poulation at the BCPs, the author provides an unambiguous method for assessing and characterising covalency in these compounds. Increasing at the beginning of the series, the covalency reached a maximum for the elements Pa-Pu before decreasing for the later actinides. ρ at the BCP can be considered as an indicator of covalent interaction, in line with an accumulation of electronic charge density in the bonding region. While 5f contributions to covalency in these complexes are smaller in magnitude than 6d contributions, the variation in covalency is almost entirely accounted for by the variation in the 5f contribution [198].

Concluding remarks
In this review, we summarized the current state-of-the-art of methods dedicated to the theoretical characterisation of bonding. They are well established for organic, inorganic or organometallic compounds. However, some of them cannot be applied to characterise bonding in actinide compounds, but, the recent methodological developments suggest that there will be significant advances in the near future. The f-block still remains a challenge for the theoretical chemist. The very nature of the interaction is far from completely understood. Multiple analysis methods (charge, orbital, quantum chemical topology and energy decomposition analyses) are required to obtain a complete description of bonding. While several simple approaches have been mentioned already, the complexity of the electronic structure of f-element systems makes that the validity of the theoretical models employed must be checked before drawing conclusions. Understanding covalent mixing in actinide metal-ligand bonds is of particular interest in nuclear sciences and more generally for studies associated with environmental releases of radionuclides and, e.g., to ultimately allow for the design of more selective extractants.       (a) (b) Figure 6: Molecular structure of BTP (a) and the isoamethyrin dianion (b). Symmetry-distinct coordinating nitrogens are labelled (figure was reproduced from ref. [187] with permission of the copyright holders).  [198] with permission of the copyright holders).