Exploring the Interaction Natures in Plutonyl (VI) Complexes with Topological Analyses of Electron Density

The interaction natures between Pu and different ligands in several plutonyl (VI) complexes are investigated by performing topological analyses of electron density. The geometrical structures in both gaseous and aqueous phases are obtained with B3LYP functional, and are generally in agreement with available theoretical and experimental results when combined with all-electron segmented all-electron relativistic contracted (SARC) basis set. The Pu–Oyl bond orders show significant linear dependence on bond length and the charge of oxygen atoms in plutonyl moiety. The closed-shell interactions were identified for Pu-Ligand bonds in most complexes with quantum theory of atoms in molecules (QTAIM) analyses. Meanwhile, we found that some Pu–Ligand bonds, like Pu–OH−, show weak covalent. The interactive nature of Pu–ligand bonds were revealed based on the interaction quantum atom (IQA) energy decomposition approach, and our results indicate that all Pu–Ligand interactions is dominated by the electrostatic attraction interaction as expected. Meanwhile it is also important to note that the quantum mechanical exchange-correlation contributions can not be ignored. By means of the non-covalent interaction (NCI) approach it has been found that some weak and repulsion interactions existed in plutonyl(VI) complexes, which can not be distinguished by QTAIM, can be successfully identified.


Introduction
Plutonium is the only element in the periodic table that can have appreciable amounts of four different oxidation states existing in aqueous acidic solutions simultaneously. Therefore, the structures, redox, and interactive features of plutonium complexes are needed to be investigated by using experimental or theoretical methods. Plutonium complexes were significantly less studied in experiments relative to its thorium and uranium counterparts for reasons of toxicity and radioactivity. Raman spectra of plutonyl (VI) ions in solution were measured by Madic et al. [1] to identify the vibrational transition.Pu-L 3 X-ray near edge absorption spectra for Pu (0, VII) were measured by Conradson et al. [2] for a large number of compounds. A number of novel and unexpected behaviors were observed in their experiment. The chloride plutonyl (VI) complexes were investigated in NaCl solutions using conventional absorption spectrophotometry [3]. The experiment suggests that mixed chloro-hydroxo or chloro-carbonate may form at high pH due to strong affinity of plutonyl (VI) for chloride. Recently, the crystal structure of the plutonyl (VI) dinitrate complex was characterized with spectroscopic methods by Gaunt et al. [4]. They found the dominant formation of mononitrate and negligible formation of the dinitrate in aqueous solutions. In contrast to experiments, tremendous theoretical works with quantum mechanical methods have been performed to probe the structural and electronic properties and spectra of plutonyl complexes. There are a number of theoretical data on plutonyl aquo [5][6][7], carbonate [8], hydroxo [9] and fluoride [10,11] complexes. Odoh et al. [12] have recently utilized DFT methods to investigate the structural properties of various Pu (V) and Pu (VI) complexes. This and other previous reports would suggest that Kohn-Sham DFT can be successfully used in predicting the structures of Pu (IV) and Pu (VI) complexes. The chemical bonding natures of Uranyl (VI) complexes have been recently identified by Vallet et al. [13] with electron density topological analyses methods, namely the quantum theory of atoms in molecules (QTAIM) proposed by Bader [14]. Kaltsoyannisaŕ group [15] also performed topological analyses for actinides compounds. These reports significantly indicate that QTAIM is an appropriate and successful tool to access the bonding characters of actinides complexes.
In this paper, we first determined the structures of various plutonyl (VI) complexes with OH − , Cl − , (F − /Br − ), O 2− 2 , NO − 3 , CO 2− 3 , H 2 O ligands in gaseous and aqueous phases within DFT framework. Based on the gaseous electron density, the QTAIM and interaction quantum atom (IQA) [16] analyses were performed to identify the interaction nature between Pu and coordinating ligands. In addition, another density-based analytical tool, the electron localization function (ELF) [17], was utilized to study the chemical bondings. As an excellent extension of QTAIM theory, the non-covalent interaction (NCI) index proposed by Yang et al. [18] was used to probe into the weak interactions in the complexes studied.

Computational Details
Calculations were performed with two different packages, ORCA-3.0.0 [19] and G09 [20]. All calculations in gaseous and aqueous phases were done at DFT level of theory with functional forms of B3LYP [21]. The geometry structures of titled complexes in gaseous and aqueous phases were optimized without symmetry constrain, and were identified to be one minimum in potential energy surface through vibrational frequency calculations. All solution phase calculations were carried out in water.
In ORCA calculations, the scalar relativistic effect was evaluated with zero-order regular approximation (ZORA) [22]. Segmented all-electron relativistic contracted (SARC) Gaussian-type basis sets [24] are assigned to Pu atom for optimization, and other atoms in the studied complexes were considered with all-electron TZV-ZORA basis set [23]. The conductor-like screening model (COSMO) [25] is utilized to evaluate the solvent effect. In G09 calculations, the B3LYP functional in conjunction with relativistic effective core pseudopotential (RECP) was utilized for optimization, the Stuttgart small-core relativistic effective core potential (RECP) with 60 core electrons frozen and assisted basis set [26,27] were employed to describe Pu atom, and aug-cc-PVDZ basis set [28] for light elements. The ultrafine grid (99, 590) was used in the integral calculations to accurately describe heavy element Pu. The solution phase calculations in G09 were carried out with the polarizable continuum model (PCM) [29]. The calculational scheme with G09 is defined as B3LYP/SDD in the following.
The QTAIM, IQA, energy decomposition scheme and NCI analyses were performed to deeply understand the interaction natures between Pu and different ligands in studied plutonyl (VI) complexes. The wave-function used in these bonding analyses were obtained at the B3LYP level of theory. The QTAIM and IQA calculations were carried out with AIMAll [30], and NCI with Multiwfn [31] programs.

Structures in Gaseous and Solution Phases
We first optimized the structures of bare plutonyl moieties (PuO 2/1+ 2 ) and pentaaquo complexes to test the reliability of DFT functionals and basis set employed in our present work. The calculated structural parameters are listed in Table 1. Our results indicate that that the Pu-O yl bond lengths calculated with B3LYP/SARC-ZORA are closer to experimental values [1,32] relative to the B3LYP/SDD results. Our B3LYP results are in agreement with other theoretical values obtained at ab initio [33,34] and DFT levels [12]. The consideration of solution effect results in the elongation of Pu-O yl bonds and the shortening of Pu-OH 2 bonds for all studied aquo complexes. Similar effects have been observed in recent DFT calculations [12]. It is worth noting that the Pu-O yl and Pu-OH 2 bond lengths are closer to the experimental values when the solution effect is considered. The geometrical structures of the plutonyl (VI) complexes with ligands, F − , Cl − , Br − , OH − , NO − 3 , O 2− 2 and CO 2− 3 obtained at the B3LYP/SARC-ZORA level of theory are depicted in Figure 1. The typical bond lengths calculated with B3LYP/SARC-ZORA and B3LYP/SDD, respectively, are collected in Table 2 for the sake of comparison. On the whole, the Pu-O yl bond lengths calculated with solution effect considered are slightly larger than those obtained in gas-phase calculation for all studied complexes. On the other hand, the consideration of solution results in contraction of the Pu-OH 2 bond lengths and elongation of of Pu-Ligand bonds. In PuO 2 F 2− 4 , PuO 2 Cl 2− 4 , PuO 2 Br 2− 4 , PuO 2 (OH) 2− 4 complexes, the Pu-L bond lengths show the order of Pu-F<Pu-OH<Pu-Cl<Pu-Br, and the Pu-F, Pu-Cl, Pu-Br and Pu-OH bond lengths in aqueous phase were calculated as 2.201, 2.706, 2.853 and 2.243 Å, respectively, when the B3LYP/SARC-ZORA was employed. We also investigated the geometry properties of mono-and bis-chloro complexes which have been identified in NaCl solutions by Runde et al. [3]. Energetically, two isomers of PuO 2 Cl 2 (H 2 O) 2 are very close to each other, with the difference of 4.82, 0.19 kcal/mol, respectively, in gas-phase and aqueous phase at B3LYP/SARC-ZORA level of theory. Both isomers show the similar structural features as shown in Table 2. The trans-isomer of PuO 2 Cl 2 (H 2 O) 3 complexes is slightly more stable than cis-isomer by 8.32 kcal/mol (gas), 2.11 kcal/mol (aquo), respectively. There also exists no great difference in geometry parameters for both isomers. It should be noted that the structural parameters of PuO 2 Cl 2 (H 2 O) 3 in aqueous phase with B3LYP/SARC-ZORA are excellent in agreement with available experimental values [3]. The mono-chloro complex (PuO 2 Cl(H 2 O) + 4 ) possesses shorter Pu-O yl and Pu-OH 2 bonds relative to bis-chloro complexes. This result agrees well with previous experimental [3] and theoretical investigations [12]. The distance between Pu and Cl atoms in PuO 2 Cl(H 2 O) + 4 is 2.722 Å at B3LYP/SARC-ZORA level, larger than previous DFT calculations [12] by 0.08 Å, is however closer to the experimental value of 2.75 Å. In both gaseous and aqueous phases, the bond lengths of Pu-O yl and Pu-OH − in hydroxide complex, PuO 2 (OH) 2 (H 2 O) 3 , are significantly shorter than that in PuO 2 (OH) 2− 4 . For PuO 2 (O) 2 (H 2 O) 3 complex, the distance between Pu and coordinating oxygen atom of O 2− 2 ligand is longer than Pu-OH − bond length in PuO 2 (OH) 2 (H 2 O) 3 complex. It should be mentioned that the U-OH − bond is shorter than U-O 2− 2 bond in respective triaquo U(VI) analogies. This indicates that plutonyl (VI) shows weak affinity to O 2− 2 ligand relative to plutonyl (VI). In experiments, Gaunt et al. [4] recently characterized the structure of Plutonyl (VI) diaquo and dinitrate. The Pu-O yl , Pu-O L (L = NO − 3 ) and Pu-OH 2 bonds lengths were calculated with B3LYP functional as 1.767, 2.483 and 2.466 Å, respectively, and are excellent in agreement with experimental measures [4] shown in Table 2. The typical bonds in trinitrate complex are not significantly different from that in dinitrate complex. For plutonyl (VI) carbonate complexes, the Pu-O yl , Pu-O L (L = CO 2− 3 ) and Pu-OH 2 bond lengths are all larger than the counterparts in plutonyl (VI) nitrate complexes.

Natural Bonding Orbital Analyses
The charge (q) based on natural population analysis (NPA) [33] and electronic configurations (EC) of Plutonium, Oxygen in plutonyl (O yl ) and coordinating atoms in ligands (L) of studied Plutonyl (VI) complexes are summarized in Table 3. The Pu atoms with about 1e lost act as electron donor in all complexes as expected. One can find from Table 3 that the Pu atoms possess 5 f 5 6d 1.5 7p 0.5 -like electron configuration. Comparing with the population of free Pu atom, 5 f 6 7s 2 , there are significant electron promotions from 7s orbital to 6d and 7p orbitals when Pu atoms are involved in chemical bonds. On the other hand, the oxygen atoms in plutonyl moiety and coordinating atoms in studied Lewis acid ligands accept electrons from Pu atoms. The Wiberg bond index (WBI) in Table 3 indicates that the Pu-O yl bonds in studied complexes have double bond features. Other chemical bonds-Pu-F, Pu-Cl, Pu-Br and Pu-OH-show single bond nature. In the peroxide, nitrate and carbontrate complexes, the bond orders between Pu and coordinating oxygen are significantly smaller than other Pu-L bonds. We also found that the interaction between Pu and H 2 O ligand is very weak from WBI values.
The WBI values of Pu-O yl bonds vs. corresponding bond lengths and charge of oxygen atoms in plutonyl moiety are depicted in Figure 2, from which one can find that there is good linear correlation with R 2 = 0.94 between Pu-O yl bond lengths and corresponding WBI values. This can be easily understood that the bigger the bond orders, the shorter the Pu-O yl bond lengths. It is very interesting that the Pu-O yl bond orders also show perfect linear dependence on the charge of oxygen atoms in plutonyl moiety (R 2 = 0.98). The oxygen atoms in Pu-O yl chemical bonds which accept fewer electrons generally possess the larger bond orders. Therefore, the Pu-O yl chemical bonds are not mainly derived from charge transfer mechanism. This implies that the interactive nature between Pu and O yl is covalent predominantly, not ionic features. Table 3. The Natural population analyses (NPA) charge (q) and electronic configurations (EC) of Plutonium, Oxygen in plutonyl (O yl ) and coordinating atoms in ligands (L) of studied Plutonyl (VI) complexes, the Wiberg bond index (WBI) of Pu=O and Pu-L chemical bonds.

QTAIM Topological Analyses
The wavefunction obtained at B3LYP level of theory is used to analyze the electron densities and further topological properties within QTAIM framework for Plutonyl (VI) species studied. The gas-phase electron densities corresponding to (3, −1) bonding critical points (BCPs) in the equatorial plane of the studied complexes are depicted in Figure 3. In bare PuO 2+ 2 moiety, there is significant electron accumulation between Pu and O atoms, consistent with WBI results mentioned above. However, there is practically no built-up of electron density between Pu and H 2 O/CO 2− 3 /O 2− 2 /NO − 3 , implying the closed-shell interaction characteristics. It should be mentioned that previous works found significant electron density accumulation between Uranium and coordinating oxygen atoms in O 2− 2 and CO 2− 3 ligands. However, Plutonyl peroxide and carbonate complexes do not show similar interaction feature to Uranyl analogies. [13]. For PuO 2 F 2− 4 , PuO 2 Cl 2− 4 , PuO 2 Br 2− 4 , PuO 2 (OH) 2− 4 complexes, the minor electron accumulations is presented between Plutonium and F − /Cl − /Br − /OH − .  To deeply understand the interactive natures between Pu and coordinating atoms in complexes studied, the electron density (ρ) at (3,−1) BCPs, Laplacian of electron density (∇ 2 ρ), energy density H(r), delocalization indices δ(A, B) and −V(r)/G(r) ratio are collected in Table 4. According to QTAIM theory, the covalent interaction has a negative ∇ 2 ρ at the critical point (CP). However, this criterion has been proved to be not sufficient to describe the bond natures of heavy atoms in previous work [36]. Another property, the total energy density H(r) (defined as the sum of local kinetic energy density G(r) and the local potential energy density V(r)) proposed by Cremer [37] was proved to be very appropriate to characterize the degree of covalency of a bond, the more negative the H(r) value, the more stabilizing the interaction. One can also employ the −V(r)/G(r) ratio as another useful description, the −V(r)/G(r)< 1 is characteristic of a typical ionic bond; and −V(r)/G(r)> 2 is diagnostic of a classical covalent interaction. As shown in Table 4, the ∇ 2 ρ values are positive for all studied Pu-L BCPs, and show the sign of closed-shell interaction nature. The strong Pu-O yl covalent interaction is clearly revealed by great negative H(r) value, large δ(Pu, L) and −V(r)/G(r) ratio being close to 2. The positive ∇ 2 ρ and H(r) values of Pu-OH 2 BCPs in all complexes studied strongly suggest closed-shell interactions, which is also supported by −V(r)/G(r)< 1 and low WBI and δ values. In the Fluoride, Chloride, Bromide and Hydroxide complexes, the value of H(r) is negative with 1.05 ≤−V(r)/G(r)≤ 1.31, indicating weak covalent natures of Pu-F − /Cl − /Br − / OH − bonds. It is in particular important to note that there are more electron accumulations located at Pu-OH − BCPs in PuO 2 (OH) 2− 4 complexes relative to other Pu-ligand bonds, and values of H(r) are slightly larger in absolute values than for the other ligands. This indicates that the affinity between plutonyl (VI) and OH − group is strongest among studied ligands. It is also worth noting that the interaction between Pu and Cl atoms in Plutonyl (VI) chloro complexes are significantly enhanced when the H 2 O ligands are coordinated, as illustrated by the H(r) and −V(r)/G(r) values in Table 4

Interaction Quantum Atom (IQA) Analyses
To interpret the physical nature of the Pu-L interatomic interaction and its local contribution to the overall energy of one Plutonyl (VI) complex, the computationally expensive but highly useful IQA energy decomposition scheme is considered. IQA defines the interaction E AB int between two atoms as a competing contribution made by classical components (interaction energy between electrons and nuclei as well as coulombic interaction between electrons of atom A and B) conveniently combined as V AB cl , and quantum-mechanical contribution, as V AB xc . The calculated decomposition of Pu-L interaction energy in the studied Plutonyl (VI) complexes is collected in Table V, from which one can see that all Pu-L bonds except for Pu-N(NO − 3 ) and Pu-C(CO 2− 3 ) shows negative E AB int values, suggesting the stabilizing contribution. Among all Pu-L chemical bonds,the strong Pu-OH bond in PuO 2 (OH) 2− 4 is characterized by the largest |E AB int | value. This result is in line with QTAIM parameters and WBI values mentioned above. When the relative values of V AB cl were considered, we found that all Pu-L bonds except for Pu-N L (L = NO − 3 ) and Pu-C L (L = CO 2− 3 ) are predominately derived from the electrostatic attractive interaction. However, it is important to note that the V AB xc values from quantum-exchange correlation are in the range of −38.91 to −81.58 kcal/mol. These quantities can not be ignored and also make a significant contribution to total interaction energy. Therefore, the quantum-exchange term is also important for the stabilizing of interaction. On the other hand, the Pu-N and Pu-C interactions are characterized by extremely large electrostatic repulsion with positive values of V Pu−N cl = 247.87, V Pu−C cl = 596.13 kcal/mol, respectively, which is responsible for ring critical point (RCP) in QTAIM. We also calculated the V AB xc /E AB int ratio, conveniently called the exchange-interaction ratio (EIR). The V AB xc term is always negative, therefore, when EIR> 0, it can identify a locally stabilizing interatomic interaction, while locally destabilizing is characterized with EIR< 0. From Table 5

Electron Localization Function (ELF) Analyses
We now consider the behavior of another density-based analytical tool, the electron localization function (ELF, denoted by η(r)). η(r) is a scalar field, and its topology can be analyzed in a manner similar to that of ρ(r) [38]. According to the ELF criterion, for covalent bonding, a typical maximum value of the ELF connecting two atoms is closer to 1.0, the higher the ELF value is, the stronger the bond is. The values of ELF associated with (3,−1) BCP between Pu and O yl and ligand are collected in Table 6, one can see that the ELF values corresponding to Pu-O yl bonds are significantly larger than those of Pu-ligand, suggesting the strong covalent nature of Pu-O yl bonds. These results correlate excellently with conclusions from QTAIM. Moreover, one can find the variation of ELF in these Pu-O yl bifurcation points is very small, indicating that different ligand fields do not significantly affect the Pu-O yl interaction-for the ELF values associated with Pu and ligands suggest the weak interaction between them. Comparing ELF values, one can obtain that Pu-OH − , Pu-F − ,Pu-Cl − ,Pu-Br − bonds are slightly stronger than Pu-L bonds in rest complexes studied.

Noncovalent Interaction (NCI) Analyses
NCI index based on the relation between the electron density and reduced density gradient (RDG) was proposed by Yang et al. [18] to identify the weak interactions. The RDG is defined as follows: NCI resorts to the sign of the second eigenvalue of the Hessian (λ 2 ) of the electron density due to its sign dependence on different situations. More specifically, the λ 2 value will be negative in bonding cases and positive in most of the non-bonding repulsive situations. Generally, NCI procedure maps the product of the sign of the second eigenvalue of the Hessian, sign (λ 2 ) times the density onto an isosurface of RDG (s), accessing a pictorial display of the most relevant attractive (blue) and repulsion interaction (red) in a system. In addition, the 2-dimension (2D) scatter plots between sign(λ 2 ) * ρ and RDG (s) can be used for NCI analyses. NCI allows for distinguishing which of the weak interactions are attractive or repulsive, and it can reveal those interactions that are not easily detected by QTAIM.
As we know that the preference for actinides to form complexes generally follows their effective ion charges, in the order of AnO 2 2+ > AnO 2 + , therefore, the PuO 2 (H 2 O) 2+ 5 (VI) complex was more favored in stability than PuO 2 (H 2 O) + 5 (V). As  2+ 5 (VI), which is responsible for its higher preference in stability. Apart from the bonding interaction between Pu and OH 2 , there are two gradient spikes in s(ρ) of Figure 4 (bottom), one is attractive interaction derived from the van der Waals interaction between neighboring H 2 O molecules (around −0.01 spikes), and the other corresponds to steric repulsion (around 0.01 spikes) between Pu and neighboring two H 2 O molecules. For other studied complexes with Cl − , OH − , O 2− 2 , CO 2− 3 and NO − 3 ligands, the NCI isosurfaces and corresponding 2D scatter plots are pictorially presented in Figure 5. Both 3D isosurfaces and 2D scatter graphs show that bonding interaction between Pu and OH − which corresponds to low gradient spike of -0.10 is significantly stronger than Pu-Cl − and Pu-O 2− 2 bonds, and the interaction strength sequence is OH − > Cl − > O 2− 2 , in line with QTAIM and IQA results. For carbonate and nitrate complexes, the low gradient spike corresponding to Pu-O(NO 3 − ) interaction in nitrate complexes is directed toward more negative sign(λ 2 ) * ρ values relative to that in carbonate complexes, indicating the slightly stronger Pu-O(NO 3 − ) bonding in nitrate complexes. We also note that there exist significant steric interactions between Pu and C/N atoms in both carbonate and nitrate complexes which have been revealed by IQA analyses, and this correlates very well with the QTAIM analyses where the bond path between Pu and C/N atoms are not observed. We also should mention that the Pu-OH 2 interaction strengths, corresponding to the peaks from −0.01 to −0.05 are very different for studied complexes.

Conclusions
The geometric structures of plutonyl(VI) complexes in gaseous and aqueous phases were obtained with B3LYP functional in conjunction with relativistic effective core pseudopotential (RECP) and SARC basis set, respectively. The predicted geometry structures by all-electron SARC basis set are more accurate than RECP results. The consideration of solution effect results in the elongation of Pu-O yl bonds and the shortening of Pu-OH 2 bonds for all studied aquo complexes. Moreover, the Pu-O yl and Pu-OH 2 bond lengths are more close to the experimental values when solution effect is considered.
NPA analyses showed that there are significant electron promotions from 7s orbital to 6d and 7p orbitals when Pu atoms are involved in chemical bonds. The Pu-O yl chemical bonds, in which oxygen atoms accept fewer electrons, correspond to the stronger bond strength. The interesting linear dependence between Pu-O yl bond orders and the charge of oxygen atoms revealed the covalent characters of Pu-O yl bonds.
By employing the wavefunction obtained at B3LYP level of theory, we have investigated the binding natures of a series of plutonyl (VI) complexes utilizing QTAIM, IQA, ELF and NCI methods. Our QTAIM analyses indicate that the significant closed-shell interactions are suggested for Pu-Ligand bonds in almost all complexes studied. It is worth noting that the Pu-Ligand interaction shows a slight degree of covalency in some complexes. Moreover, the ligand affinity of plutonyl (VI) is different from uranyl(VI) analogies. IQA energy decomposition scheme was utilized to probe into the interaction between Pu and ligands, the electrostatic attractive interaction is characterized as the dominate stabilizing contribution, but the quantum-exchange term is also non-ignorable, making obvious contribution to total interaction energy. ELF was also employed to identify the bonding natures, and was found to draw similar conclusions to QTAIM.The NCI index results indicate that NCI is a successful extension of QTAIM. Some weak interactions, like Pu-OH 2 , H 2 O-OH 2 , which can not be distinguished in QTAIM, can be successfully identified with NCI. Moreover, the steric repulsion can be found with NCI analyses. We hope our electron density topological analyses are useful for a deep understanding of Pu-Ligand chemical bondings in plutonyl (VI) complexes.