Investigation of the water exchange mechanism of the Plutonyl(VI) and Uranyl(VI) ions with quantum chemical methods

The water exchange reactions of [PuO2(OH2)5]2+ and [UO2(OH2)5]2+ were investigated with density functional theory (DFT) and wave function theory (WFT). Geometries and vibrational frequencies were calculated with DFT and CPCM hydration. The electronic energies were evaluated with general multiconfiguration quasi-degenerate second-order perturbation theory (GMC-QDPT2). Spin-orbit (SO) effects, computed with SO configuration interaction (SO–CI), are negligible. Both Actinyl(VI) ions react via an associative exchange mechanism, most likely Ia. The Gibbs activation energies (ΔG‡) at 25 °C are 33–34 and 30–37 kJ mol−1 for [PuO2(OH2)5]2+ and [UO2(OH2)5]2+, respectively. ΔG‡ for dissociative mechanisms (D, Id) is higher by more than 15 kJ mol−1. Graphical Abstract


Introduction
In the Plutonyl(VI) aqua ion, five water molecules are coordinated in the equatorial plane [1]. The exchange of these water molecules, reaction (1), was measured with 1 H NMR by Bardin et al. [2], whereby the rate constant extrapolated to 25°C was estimated to be approximately 10 4 s −1 .
For the corresponding reaction (2) of the Uranyl(VI) aqua ion, they obtained the same value at 25°C, 10 4 s −1 .
It is lower by about two orders of magnitude than the rate constant of (1.30 ± 0.05) × 10 6 s −1 measured by Farkas et al. [3] via 17 O NMR. In this study, erroneously [4], the dissociative (D) mechanism was attributed to reaction (2). Assuming a transmission coefficient (κ) of 1, they obtained a Gibbs activation energy (ΔG ‡ ) of 38.0 kJ mol −1 .
Later, Kerisit and Liu [5] studied reaction (2) with classical molecular dynamics (MD) simulations and determined ΔG ‡ = 35.1 kJ mol −1 and a transmission coefficient κ = 0.194 via the reactive flux method (using model 1 in Ref. [5]). The exchange mechanism of reaction (2) was attributed based on quantum chemical calculations. Among numerous studies, the most reliable ones using the most advanced quantum chemical methods and implicit solvation or the most elaborate treatment of hydration are presented. Bühl and coworkers [6,7] performed ab initio MD simulations of [UO 2 ] 2+ with 59 H 2 O molecules in a cubic box using the BLYP functional [8,9] and periodic boundary conditions. For the associative (A) and the less favorable dissociative mechanisms, they obtained Helmholtz activation energies of 28.0 and 45.2 kJ mol −1 , respectively. In an ab initio molecular orbital (MO) study [4], geometries and vibrational frequencies were computed at the complete active space self-consistent field (CASSCF) level, whereby hydration was treated using the polarizable continuum model (PCM) [10,11]. Total energies were calculated with the multiconfiguration quasi-degenerate second-order perturbation theory [12,13], yielding Gibbs activation energies of 35.3 and 57.1 kJ mol −1 for the A and D mechanisms, respectively. The value for A is in perfect agreement with the above-mentioned experimental value and, hence, an associative mechanism operates for reaction (2).
While the water exchange rate and mechanism of Uranyl(VI) are established, this is not the case for Plutonyl(VI). Since the exchange rate constant for reaction (2) measured by Bardin et al. [2] is too low by a factor of about 100, the corresponding value for reaction (1) might also be inaccurate. In classical MD simulations of Actinyl ions in water, Tiwari et al. [14] calculated ΔG ‡ for reactions (1) and (2) following the associative interchange (I a ) pathway as 28.7 ± 0.8 and 21.0 ± 0.1 kJ mol −1 . Since the value for (2) is underestimated by 13-14 kJ mol −1 , the value of 28.7 kJ mol −1 for Plutonyl(VI) cannot be considered reliable without reservation. To assess the Gibbs activation energy and the mechanism of reaction (1), these two water exchange reactions were studied with density functional theory (DFT) and wave function theory (WFT). was used. For O and H, modified [19] Karlsruhe def2-TZVP basis sets [20,21] were taken. Figures 1 and 2 were generated with MacMolPlt [22].
The transition states (TS) were located by maximizing the energy along the reaction coordinate (the imaginary mode) via eigenmode following. They exhibited a single imaginary frequency. Reactants, products, and intermediates were obtained via computation of the intrinsic reaction coordinate. All of their computed vibrational frequencies were real. Selected atomic coordinates of the investigated species are given in tables S1-S15 (see online supplemental material at http://dx.doi.org/10.1080/00958972.2015.1059425).
Methods for the calculation of the Gibbs activation energy (ΔG ‡ ) for associative and dissociative exchange mechanisms The water exchange reactions (1) and (2) can be formulated as unimolecular processes using water-adducts [34][35][36]. Alternatively, the reactions may involve a water molecule from the bulk solution.
Solvent molecules are not involved in the activation step of the dissociative (D) pathway (5a) [37][38][39], but a water molecule is eliminated in the formation of the respective intermediate (5b).
ΔG ‡ is calculated using equation (6): ΔE ‡ is the electronic activation energy in the gas phase, ΔG solv ‡ is the difference in solvation energies, and Δg° ‡ is the difference of the thermal corrections to the Gibbs energy (including the zero point energy) at 1 atm and 25°C. ΔE ‡ and ΔG solv ‡ are calculated at the geometries of the solvated species and Δg° ‡ is based on the vibrational frequencies of the solvated species. ΔG for reactions (3b) and (5b) is calculated analogously.
As already mentioned, bulk solvent molecules are not involved in reaction (5a). The other reactions, (3), (4) and (5b), are rewritten to take into account the participation of water from the bulk solvent: Reactions (7), (8) and (9) are not unimolecular. Therefore, equation (6) cannot be used for the calculation of ΔG ‡ and ΔG without corrections for the standard states in the gas phase (1 atm) and solution (1 M). For reactions (7) and (8), ΔG ‡ is obtained via equation (10) [40]: V°is the volume of 1 mol (ideal gas) at 1 atm and 25°C, and V* is equal to one liter. Because the concentration of the water solvent is not 1 M, the RT ln[H 2 O] term is also required. Since a water molecule is liberated in equation (9), the corrections have the opposite sign, and ΔG has to be evaluated via equation (11): Both methods exhibit strengths and limitations. With water-adducts, all of the reactions are unimolecular. Hence, no corrections for the various standard states are required. The differences in electronic energies correspond approximately to those of the Gibbs energies. As a disadvantage, the number of isomers is larger for the water-adducts. Therefore, elucidation of the global minimum is more demanding. Furthermore, the bridged structures of the water-adducts do not correspond to the structures in solution, and the reactants for the A, I a , I, and I d mechanisms differ from that for D. For the reactions with water from the bulk solution, corrections for the different standard states in the gas phase and solution are required, and the solvent concentration has to be taken into account as well. The gas phase electronic energy differences may differ strongly from the corresponding Gibbs energy differences in solution. As a benefit, all of the mechanisms are tractable starting with the same reactant structure.
In this study, reactions (1) and (2) were investigated using both methods.

Gibbs activation energies (ΔG ‡ ) for associative and dissociative water exchange mechanisms
To assess the accuracy of the above-mentioned functionals, Gibbs activation energies (table 3) for associative and dissociative water exchange mechanisms of reaction (1) were computed at the DFT geometries with DFT and WFT. Due to the presence of static correlation in Actinyl complexes, multiconfiguration methods are required [4,60,61]. Thus, ΔE ‡ and ΔG solv ‡ in equations 6, 10 and 11 (ΔE and ΔG solv for the intermediates) were computed with DFT and GMC-QDPT2 and SO-CI (Computational Details), whereby the SO energy differences were negligible (<0.01 kJ mol −1 ) because both of the f δ levels are occupied by a single electron and because of the absence of low-lying triplet and singlet (excited) states. WFT hydration was calculated with MCSCF and CPCM hydration. Δg° ‡ or Δg°was always computed with DFT (because the geometries were optimized with DFT). The relative Gibbs energy (ΔG) of the intermediates for the A and the D pathways with an increased or a reduced coordination number are also included in table 3. It should be noted that the LC-BOP-LRD and the ωB97X functionals yielded TS for the I a mechanism (figure 1), whereas with the other functionals, PBE0-D3, TPSSh, TPSSh-D3, and LC-BLYP, TSs for the A pathway (figure 2) were obtained. ΔG ‡ for A is somewhat higher than for I a . The distinction of the A from the I a , and the D from the I d mechanisms is very demanding [4,36]; this issue will be discussed later. For the D mechanism, WFT and DFT yielded considerably different ΔG ‡ values. The differences are smallest for LC-BOP-LRD and ωB97X. Virtually the same differences (within 1 kJ mol −1 ) were observed for the corresponding electronic energy differences (ΔE ‡ and ΔE) in the gas phase. Hence, these differences are not due to errors in the CPCM hydration. For the I a and the A mechanisms, the WFT-DFT differences in ΔG ‡ are smaller. For all of the investigated mechanisms, D, I a , and A, the best ΔG ‡ /ΔG and ΔE ‡ /ΔE agreement was obtained with LC-BOP-LRD, followed by ωB97X.
While the computed ΔG ‡ values based on DFT do not allow the attribution of the exchange mechanism of reaction (1), WFT clearly shows that ΔG ‡ for I a and A is lower than for D (table 3). For the most reliable functionals LC-BOP-LRD and ωB97X, I a is favored over D by > 15 kJ mol −1 . ΔG ‡ for I a and A, and ΔG for A and D were evaluated based on both methods, that involving the water-adducts, equations (3)- (5), and the other one with the participation of water from the bulk solution, equations (5a), (7)-(9). (ΔG ‡ for D does not involve water from the bulk solution). For the I a and the A mechanisms, ΔG ‡ calculated using both of these methods differs considerably for some functionals. The same holds for ΔG for D. It is noteworthy that in the water-adducts, there are two hydrogen bonds, which, in general, are difficult to treat accurately with DFT [62].
The intermediates for the A and the D mechanisms should exhibit a lower Gibbs energy than the corresponding TSs. This is, however, not always the case (table 3) because of the approximations in the calculation of ΔG ‡ and ΔG. These are for example the geometry optimizations, for which the electronic energy is minimized for the rigid species at 0 K (not G at 25°C), the equations of statistical mechanics, which were derived for the free species in the gas phase (although the vibrational frequencies of the solvated species were used for the calculation of Δg° ‡ and Δg°), the neglected variations of the thermodynamic properties of the bulk solvent as well as specific solute-solvent interactions, and the harmonic approximation. In cases with a small Gibbs energy difference between TS and intermediate, which gives a hypothetical life-time of the intermediate being smaller than the duration of the vibration leading to the reactant or the intermediate, the mechanism would be rather I a or I d than A or D [4]. Thus, for the italicized entries in table 3, for which the intermediate exhibits a higher energy than the corresponding TS, and for the entries with small ΔG ‡ -ΔG differences as well, the I a or the I d mechanism might operate instead of A or D.
As mentioned above and in previous work [4,36], the distinction of the concerted (I a , I d ) from the step-wise (A, D) mechanisms is subtle. For reaction (1), depending on the functional, TSs for the I a (figure 1) and the A (figure 2) mechanisms were obtained (table 3). Due to the lower ΔG ‡ for I a and the small or even negative ΔG ‡ -ΔG difference for A, it is preferable to attribute the I a mechanism to reaction (1). As already seen, ΔG ‡ and ΔG based on DFT are not reliable and DFT does not predict the preference for an associative mechanism of reaction (1) unequivocally. In Actinyl(VI) complexes, static electron correlation is present; its appropriate treatment [4,[59][60][61] is mandatory to obtain reliable results. For reaction (1), the most reliable data were obtained with LC-BOP-LRD or ωB97X geometries and frequencies, and GMC-QDPT2/SO-CI energies. Hence, ΔG ‡ for the I a and the unfavorable D or I d mechanisms is 33-34 and 46-49 kJ mol −1 , respectively.
To assess the accuracy of the present computations, the well-studied Uranyl(VI) water exchange reaction (Introduction) was re-investigated using the most accurate functionals LC-BOP-LRD and ωB97X, and GMC-QDPT2 (table 4). ΔG ‡ calculated with DFT is inaccurate and lower than ΔG ‡ based on WFT, which amounts to 30-37 kJ mol −1 (table 4) for the A mechanism. This value agrees with the experimental value of 34-35 kJ mol −1 (Introduction). The D mechanism, exhibiting a higher ΔG ‡ (by > 15 kJ mol −1 ), is unfavorable. Also for this reaction (2), it is preferable to attribute the I a mechanism due to the small (LC-BOP-LRD/WFT) or negative (ωB97X/WFT) ΔG ‡ -ΔG difference. This mechanism was also attributed by Kerisit and Liu based on classical MD simulations [5].

Summary
The water exchange reactions of the f 0 and f δ 2 Actinyl(VI) aqua ions of U and Pu, respectively, proceed via an associative mechanism, most likely I a . For both reactions, ΔG ‡ is virtually equal (30-37 kJ mol −1 for U and 33-34 kJ mol −1 for Pu). For the unfavorable D pathway, ΔG ‡ is higher by more than 15 kJ mol −1 for both Actinyl(VI) ions. The reactivity (mechanism and ΔG ‡ ) of these two ions is equal in spite of their disparate, but symmetric f orbital occupations (the occupation is asymmetric in f 1 and f 3 systems).
The An=O bond lengths calculated with DFT are somewhat too short for most functionals because of an imperfect treatment of static correlation. Furthermore, DFT is often inaccurate for hydrogen bonding. Therefore, the computed ΔG ‡ and ΔG are also inaccurate. In such cases, reliable results can be obtained with state-of-the-art WFT. In systems exhibiting static correlation, the adequacy of DFT and MP2 should be assessed with WFT.
Reactions involving a solvent molecule from the bulk can be treated with or without water-adducts, provided that the various standard states are taken into account.

Supplementary material
Modified ECPs for U and Pu. Atomic coordinates of the reactants, TSs, and intermediates for the I a , A, and D mechanisms. An=O, An-O(H 2 ), and An⋯O(H 2 ) bond lengths in the I a , A, and D TSs (An=Pu, U).