Computational Tools for Calculating log β Values of Geochemically Relevant Uranium Organometallic Complexes

: Uranium (U VI ) interacts with organic ligands, subsequently controlling its aqueous chemistry. It is therefore imperative to assess the binding ability of natural organic molecules. We evidence that density functional theory (DFT) can be used as a practical protocol for predicting the stability of U VI organic ligand complexes, allowing for the development of a relative stability series for organic complexes with limited experimental data. Solvation methods and DFT settings were benchmarked to suggest a suitable o ﬀ -the-shelf solution. The results indicate that the IEFPCM solvation method should be employed. A mixed solvation approach improves the accuracy of the calculated stability constant (log β ); however, the calculated log β are approximately ﬁ ve times more favorable than experimental data. Di ﬀ erent basis sets, functionals, and e ﬀ ective core potentials were tested to check that there were no major changes in molecular geometries and Δ r G . The recommended method employed is the B3LYP functional, aug-cc-pVDZ basis set for ligands, MDF60 ECP and basis set for U VI , and the IEFPCM solvation model. Using the ﬁ tting approach employed in the literature with these updated DFT settings allows ﬁ tting of 1:1 U VI complexes with root-mean-square deviation of 1.38 log β units. Fitting multiple bound carboxylate ligands indicates a second, separate ﬁ tting for 1:2 and 1:3 complexes.


■ INTRODUCTION
The utilization of nuclear energy has resulted in widespread contamination of the environment with uranium. It has generated uranium-containing waste, which, depending on national policy, will likely be disposed of in a near-surface or deep geological repository. Uranium (U VI ) can be mobilized from these sites by groundwater as part of an alkaline (pH >9) or acidic plume (pH 3.5−6). 1−4 Understanding how the chemical composition of groundwater can influence aqueous U VI chemistry is key to constrain the subsequent mobility of U VI . In particular, it is important to characterize the effect of naturally occurring organic molecules because they are dominant in these environments.
One class of naturally occurring organic molecules is siderophores. These are released by plants, bacteria, and fungi primarily to solubilize Fe III . 5,6 However, siderophores are known to complex with uranium, 7−9 leading to uranium leaching 10−12 and desorption from mineral surfaces. 13,14 Determining the stability of possible U VI siderophore complexes is not practical because there are 500 known siderophores. Fortunately, there are very few functional groups in siderophores. These are the catecholate (1), hydroxamate (2), α-hydroxycarboxylate (3), α-aminocarboxylate (4), hydroxyphenyloxazolone (5), and α-hydroxyimidazole (6) functional groups, as seen in Figure 1. 5 Characterizing the stability of each functional group across pH 3.5−6 will allow us to begin to understand which class of siderophores form the most stable complexes with U VI in the acidic pH range of interest.
The relative stability of U VI −siderophore functional group complexes can be compared by calculating the stability constant, log β, of the reaction. The experimental determination of stability constants often requires multiple experimental techniques and is time-consuming and labor-intensive. It can also be difficult to use certain experimental techniques such as potentiometric titrations in high-pH and high-ionicstrength solutions due to the sodium error when using glass electrodes. 15 Density functional theory (DFT) provides a costeffective and radiation-free way to estimate the log β of U VI complexes because it computes the Gibbs energy (G) directly for each species in metal−ligand complexation reactions, allowing for the subsequent calculation of log β values.
The log β values for U VI ligand reactions were calculated before using DFT. The absolute values, however, were off by 10 log β units (>60 kJ/mol) 16,17 or >30 log β units. 18 Although absolute values differed from experiment, the relative values estimated chemical trends correctly. The accuracy and computational cost of log β calculations under aqueous conditions depend on four different model chemistry aspects in the DFT protocol suggested here. These are the solvation method employed, the effective core potential (ECP), the basis set, and the functional. This is represented by eq 1.
The functional employed, ECP (as long as it is a small-core ECP), and the basis set used do not significantly affect the solvation energies for closed-shell, f elements (e.g., U VI ) in solution. 19−21 However, previous work suggested that the errors in log β are sensitive to the solvation method employed. 19,20,22 Computational methods generally model water molecules that are not directly coordinated with U VI implicitly to save on computational cost. Calculated stability constants, however, are improved using a mixed implicit−explicit solvation method because it is more successful at incorporating solute−solvent cavity and dispersion terms. 23 Nevertheless this approach has to date only been explored for uranyl in detail for water exchange reactions. 19,22,24 A mixed approach was employed for 1:1 U VI acetate, oxalate, and catecholate complexes. 25 To counter the absolute errors in log β, a fitting approach is frequently employed. In this case, the computed log β of different complexes is plotted against experimental data, and the fitting curve is used to calculate fitted log β values. This was successfully applied for a variety of metal centers with ligands including siderophores. 26−28 Recently this approach was also applied to 1:1 U VI complexes for oxygen donor ligands, allowing for the prediction of log β values to within <1 (rootmean-square deviation). 18 The protocol developed in ref 18 was used in combination with the SMD implicit solvation model 29 in recent uranyl studies. 30,31 However, the protocol used a model chemistry with prohibitive computational costs when investigating ligands containing tens of atoms, and establishing a method for large ligands is key for modeling siderophore complexes. Furthermore, the solvation cavity used in the SMD model has not been parametrized for actinides. 29 The aim of the present study was to establish a computationally cost-effective protocol that allows us to estimate the relative stability of functional groups that compromise the reactive components in siderophores. The protocol needs to be sufficiently cost-effective so that it can be applied to full siderophore structures. Initially, we tested different solvation methods because the greatest error in calculated stability is expected to be due to solvation. This was tested for 1:1, 1:2, and 1:3 U VI /acetate complexes. Two implicit solvation methods were subsequently tested to examine whether a method that better calculates nonelectrostatic contributions improved the calculated complex stability. A mixed solvation approach was employed where each aqua ligand was represented as one U VI bound aqua ligand whose hydrogen atoms form one hydrogen bond to an explicit solvent molecule, [H 2 O···(H 2 O) 2 ]. This was done to investigate if the approach proposed by Gutowski and Dixon 22 improves calculated log β and whether this outweighs the difficulties associated with applying this method. After this, the effect of changing ECP, basis set size, and functional was tested to establish the magnitude of differences in log β for the three U VI acetate complexes (1:1, 1:2, 1:3). A fitting equation was then established using the seven different ligands presented in Figure 2. These seven ligands were selected because they are similar to or the same as the functional groups that comprise the reactive components of siderophores and have structural and stability data against which the protocol can be validated. Finally, the protocol was used to estimate the log β of 1:1, 1:2, and 1:3 uranyl carboxylate complexes, and potential fitting trends were analyzed. ■ METHODS Computational Details. Electronic structure calculations were carried out with Gaussian 16, rev. A.03 32 using GaussView version 6 33 for pre-and postprocessing. All structures were geometrically optimized, and the frequencies were calculated to confirm local minima and determine thermodynamic properties including Gibbs free energy. The PW91, 34 TPSS, 35 and M06 36 functionals were utilized in addition to B3LYP 37 to provide comparisons with a generalized gradient (GGA), meta-GGA, and recent hybrid  The Journal of Physical Chemistry A Article functional, respectively. Pople's polarized double-ζ basis set, 6-31G(d,p), 38,39 and Dunning's correlation-consistent, polarized valence n zeta basis sets augmented with polarization functions (aug-cc-pVnZ) were employed, where n = double (D), triple (T), and quadruple (Q) 40−42 was used for ligand atoms. Dolg's  MWB60-ECP 43 and MDF60-ECP 44 along with Peterson's cc-pVnZ basis sets for pseudo potentials (cc-pVnZ-PP) 45 were employed for U VI . Implicit solvation was accounted for with the integral equation formalism polarizable continuum model (IEFPCM) 46 and the universal solvation model built upon the integral equation formalism (SMD). 29 Combined implicit− explicit models were tested by using the [UO 2 [H 2 O··· (H 2 O) 2 ] 5 ] 2+ structure from Gutowski and Dixon. 22 This models aqua ligands as two water molecules hydrogen bonded to the aqua ligand interacting directly with the U VI metal center. Symmetry was not constrained, and all calculations utilized a pruned grid with 99 radial shells and 590 angular points per shell.
Reactions Used for Selection of DFT Components. The formation of 1:1, 1:2, and 1:3 U VI acetate complexes was used to investigate the different model chemistry components (ECP, basis set, functional, solvation method). Acetate is similar to the α-hydroxycarboxylate siderophore functional group, and the geometries of the species along with Δ r G have been comprehensively experimentally studied. 47−54 The generalized reactions for uranyl with one to three acetate ligands are given in eqs 2−4, respectively, where L is the ligand.
Reactants and products were optimized, and vibrational frequencies were calculated to confirm that structures were at energy minima. All calculations were conducted in aqueous phase, and free energies were calculated using ΔG gas + ΔG solvent = ΔG solvation . The Gibbs free energy of explicit water molecules in the bulk solution was corrected to standard-state concentration (55.5 M) by adding −17.88 or +17.88 kJ/mol (4.27 kcal/mol) depending on whether water is a reactant or product side, respectively. 55,56 Gibbs free energy (G) for each product and reactant was computed using the harmonic oscillator and rigid rotor approximations without accounting for hindered rotation. Theory for Calculating Stability Constants. Typically, the relative stability of different complexes is compared via the stability constant (β) rather than Δ r G. The simplest metal− ligand complexation reaction is of one metal ion (M) with one ligand (L) to form a 1:1 ML complex, represented in eq 5. In solution, this reaction is a ligand displacement reaction where each ligand liberates solvato ligands, Sol, from the coordination sphere of the metal ion, M(Sol) m . The values of m and n in eq 5 reflect the number of solvato ligands coordinated around the metal center and the denticity of the ligand, respectively. The stability of the complex is encapsulated by the equilibrium constant (K). Equilibrium constants are commonly reported as −log K (pK) because individual K values may span many orders of magnitude. Because metal ions have multiple coordination sites, complexation reactions often involve more than one ligand that bind in a series of reactions. The resulting stepwise equilibrium constants are multiplied together and reported as the stability constant, β. Stability constants are often reported with the number of metal ions and ligands consumed in the reaction written as subscripts. The reaction of one metal ion with one ligand is denoted β 11 . The stability constant of the complex is determined by the corresponding Gibbs free energy of reaction, Δ r G. This connection is shown in eq 6, where α M(Sol) m , α L , and α ML(Sol) m−n are activities for molecules in the reaction. The Gibbs free energy of reaction is calculated by subtracting the sum of the Gibbs free energy (G) of the reactant species from that of the product species, as shown in eq 7.
Organic Molecule Training Set Used for Fitting Calculated Stability Constants. Seven organic molecules were used for the 1:1 fitting and are shown in Figure 2. They were chosen because they are, or are similar to, siderophore chelating groups. The reaction used to calculate Δ r G is given in eq 7, except for L-glycine, which is given in eq 8. Equation 8 was used based on the reaction used in the experimental literature. 57 The log β values were calculated using eq 6. For testing the fitting approach for multiple bound ligands, the formation of different 1:1, 1:2, and 1:3 U VI carboxylate complexes was calculated. The ligands were glycolate, acetate, and oxalate. These were chosen to include changes in denticity and difference in charge.
The experimental stability constants (and resulting Δ r G) used for all comparisons were adjusted to zero ionic strength by using the specific ion theory and the Davies equation 58,59 depending on the ionic strength at which the experiments were performed. Further details can be found in the Supporting Information.

■ RESULTS AND DISCUSSION
When assessing protocol modifications, one must keep in mind that a change of one log β unit requires a −5.709 kJ/mol change in Δ r G when T = 298.15 K.
Implicit Solvation Models. For all calculations, the base model chemistry employed the B3LYP functional, aug-cc-pVDZ basis set for ligands, and MDF60 ECP with the corresponding basis set for U VI . The only exception to this rule applied when the U VI ECP and basis set was explored. In this case, the calculations used the cc-pVDZ ligand basis set to be consistent with the cc-pVDZ-PP basis set for U VI . The effects The Journal of Physical Chemistry A Article of using two different implicit solvation methods, IEFPCM and SMD, on Δ r G were studied, and the results are shown in Table  1 and Figure 3a.
The SMD model should provide more accurate Δ r G than IEFPCM. This is because the SMD model is a parametrized version of IEFPCM to better incorporate nonelectrostatic contributions such as cavitation, dispersion, and changes in local solvent structure. 29 Changing the implicit solvation model from IEFPCM to SMD changes Δ r G by +59.4, +125.0, and 155.6 kJ/mol to −59.7, −94.4, and −114.9 kJ/ mol for U VI complexation with one to three acetate ligands, respectively, as shown in Table 1 and Figure 3a. SMD is significantly more accurate for calculating Δ r G when compared with the experimental data (−19.4, −29.2, and −42.5 kJ/ mol) 48 highlighting the importance of nonelectrostatic contributions.
A downside of using the SMD method with metalcontaining systems is that the solvation cavity overexposes metal density to the bulk solvent, as seen in Figure 4a. This effect is the direct result of limited experimental data to calibrate SMD metallic radii for inorganic systems. This is an artifact of the way in which the nonelectrostatic energy contributions were included in SMD. The solvation method was parametrized by adjusting the radii distance used to create the solvation cavity. This has only been applied to nonmetallic elements typically found in organic ligands (e.g., H, C, N, O), and, to date, the solvation method has not been parametrized for transition metals, lanthanides, or actinides. 29 Therefore, the relative size of the cavity surrounding the U VI metal center to the size of the cavity around the water molecules is much larger than it should be. This is not the case for the IEFPCM model in Figure 4b.
The modified method employed recently in refs 30 and 31 should likely be used cautiously due to the solvation cavity shortcomings of the SMD model, and ideally the IEFPCM model originally employed 18 should be used to predict the relative stabilities of U VI complexes instead.
Mixed Solvation with IEFPCM. It has been observed that mixed solvation methods improve the accuracy of solvation energies. 22 This was explored for water exchange reactions between bulk water and U VI bound aqua ligands. 22 The U VI structures employed in eqs 9 to 11 are shown in Figure 5. They were chosen based on maintaining isodesmic reactions. If reactions are isodesmic, then the number of bonds is preserved. If reactions are nonisodesmic and the resulting U VI −ligand complex contains more or fewer bonds, then the stability of the complex will be higher or lower than it should be. 60 This can significantly affect Δ r G.
Using the mixed solvation method should reduce the errors compared with using the implicit method. This is because the explicit solvent molecules position the electron density of U VI further into the solvation cavity. Therefore, electrostatic interactions with the bulk solvent are reduced along with Δ r G dependence on the implicit model employed. 19,22 Because IEFPCM overpredicts solvation energies, 19 applying the mixed solvation approach should shift calculated Δ r G to less negative values and closer to the experimental data. When using the [UO 2 [H 2 O···(H 2 O) 2 ] 5 ] 2+ mixed solvation method, Δ r G changes by +30.2, +34.7, and +38.7 kJ/mol to −88.9, −184.7, and −231.8 kJ/mol for one to three acetate ligands, respectively. However, this is still 69.5, 155.5, and 189.3 kJ/ mol more negative than the experimental data for one, two, and three acetate ligands, respectively. This is shown in Figure  3a and Table 2. The Journal of Physical Chemistry A

Article
The calculated Δ r G could be improved by releasing water as a cluster rather than as individual water molecules because it would keep the reaction isodesmic. For reactions 9 and 10, the liberated water molecules are assumed to form a six-ring water structure. This was chosen rather than the three more energetically stable conformers 61 to maintain a consistent number of water molecules. In reaction 11, the final three water molecules were assumed to be released as a trimer. The standard state correction changes from −17.88 kJ/mol to −17.88/n kJ/mol for a cluster of (H 2 O) n molecules. 55,56 The reaction energies changed by +4.3, +8.5, and +4.4 kJ/mol to −84.6, −176.2, and −227.4 kJ/mol, respectively, for reactions 9−11, respectively, suggesting that there is little improvement in releasing the water molecules individually.
Mixed Solvation with SMD. One interesting idea is to combine the parametrized SMD model with the mixed [H 2 O··· (H 2 O) 2 ] method. The extra water molecules reduce the U VI atom contribution to the solvation cavity by positioning the U VI atom further into the structure, away from the bulk water. This leads to a better representation of the solvation cavity.
Using the mixed solvation method, [H 2 O···(H 2 O) 2 ], in combination with the SMD solvation model results in less accurate Δ r G than modeling solvation implicitly by −15.9, −64.1, and −70.0 kJ/mol for the 1:1, 1:2, and 1:3 uranyl acetate complexes. This is visualized in Figure 3a. This is related to the different errors in solvation. Reference 19 has explored the errors associated with adding increasing numbers of explicit water molecules. The change in free energy of solvation in a reaction, ΔG solv , is the sum of the change in the gas-phase hydration energy, ΔG hyd , and the change in the solvation correction for the reaction ΔG corr (ΔG solv = ΔG hyd + ΔG corr ). Parmar et al. 19 observed that when only aqua ligands   The Journal of Physical Chemistry A Article are complexed around U III and U IV , a large solvation correction ΔG corr is required to counter the errors associated with the gasphase hydration energy, ΔG hyd. The ΔG corr term is dominated by ΔG elec , which is sensitive to the volume and surface area of the solvation cavity. Reference 19 observed that the small cavity volume and large surface area of the SMD solvation model led to an overestimated ΔG solv . This means that the Δ r G from adding aqua ligands to a U III or U IV ion is overpredicted (Δ r G is too negative). In the U VI acetate reactions, water molecules are released rather than added; therefore, the predicted Δ r G is too positive. Addition of the explicit water molecules leads to a reduction in the magnitude of ΔG corr required as ΔG hyd converges toward ΔG solv . 19 The extra water molecules position the electron density of U VI further into the solvation cavity, reducing G corr at the edge of the solvation cavity and reducing the overestimation of solvation energies. This leads to a shift of SMD Δ r G to more negative values, as observed in Figure 3a. Conversely, ref 19 observed that the large surface area and volume of the UFF IEFPCM meant that ΔG solv is underestimated (Δ r G of adding aqua ligands is too positive); therefore, in the U VI acetate reactions where water molecules are released, Δ r G will be too negative. Adding explicit molecules in this case reduces Δ r G, and so IEFPCM   25 ] 2+ ). 9 Keeping reactions isodesmic using such a structure is challenging, and computational cost increases significantly. This would not be a practical approach for exploring a relative stability series.
The results from combining the mixed solvation method with SMD and the 6-31g(d,p) basis set for ligands are worth noting due to the high accuracy of Δ r G when comparing to the experimental data. The calculated Δ r G for one to three acetate molecules are −21.57, −29.16, and −49.98 kJ/mol compared with the experimental values of −19.4 ± 4.87, −29.2 ± 5.02, and −42.5 ± 3.61 kJ/mol, respectively. This needs to be treated with caution, however, because there can be a large truncation error for smaller basis sets. 62,63 The errors associated with the different model components could be working in favor of the reaction energetics, and the SMD solvation model is still uncalibrated in these complexes, as shown Figure 6c,d.
U VI ECP and Basis Set. Two small core ECPs have been tested, that is, the quasi-relativistic MWB60 ECP and the fully relativistic MDF60 ECP, to examine if changing U VI ECP from quasi-relativistic to fully relativistic improves Δ r G and structure geometry. In these calculations, the cc-pVDZ basis set was used for ligands while maintaining the use of B3LYP and IEFPCM. The calculated Δ r G is expected to improve when moving from MWB60 to MDF60 ECP, as MDF60 is fully relativistic, allowing for a more accurate representation of U VI valence energies. The results from these calculations are presented in Table 3 and Figure 3b. Changing from the MWB60 ECP to the MDF60 ECP resulted in a change in Δ r G of −2.0, −2.4, and +3.4 kJ/mol to −132.0, −230.2, and −329.6 kJ/mol for reactions with one to three acetate ligands, respectively. These changes are small compared with the differences in the experimental data, and they all show the same stability trends.
Changing the basis set for U VI valence electrons from the MDF60 basis set to the cc-pVDZ-PP basis set should further improve U VI valence energies and improve overall reaction energetics. The Δ r G changes by +0.6, −1.2, and −8.1 kJ/mol to −131.4, −231.4, and −337.7 kJ/mol for one to three acetate ligands, respectively, indicating no significant change. Using the MDF60 basis set is adequate in terms of Δ r G. Table 3 suggests that geometries are very similar, with average UO, U−O water , and U−O acetate bonds within 1.14, 0.49, and 0.083% of each other, respectively. The UO length is underpredicted compared with experimental data 47,52−54 by MWB60 ECP in the uranyl complex and 1:1 uranyl acetate complex by 1.4 to 0.4%, respectively. The UO bonds continue to lengthen and are over predicted in the 1:2 and 1:3 uranyl acetate complexes by 0.23 and 0.28%, respectively. The UO length is over predicted when using MDF60 ECP with and without the cc-pVDZ-PP basis set for all complexes, by up to 1.4%. The U−O acetate length increases when switching from MWB60 ECP to MDF60 ECP but remains the same when switching to the cc-pVDZ-PP basis set. These bond lengths are up to 2.8% shorter than the experimental data in the 1:1 uranyl acetate complex, 47 with lengths increasing until they are 0.6% longer than the experimental data for the 1:3 uranyl acetate complex. The U−O water lengths tend to be longer than the experimental data, increasing as acetate molecules are added, and are up to 6.7% longer in the 1:2 uranyl acetate complex. A general trend emerges where UO and U−O water bond lengths increase in the order MWB60 (shortest), MDF60, and MDF60 ECP with cc-pVDZ-PP basis set (longest). The opposite trends occur with OUO angle decreasing when switching from MWB60 ECP to MDF60 ECP. The change in OUO decreases as water is removed from the complex, with a difference in the ECP and basis set being 1.69% for [UO 2 (H 2 O) 5 ] 2+ and 0% with 180°OUO angle for [UO 2 (Acet) 3 Overall, there is little change in the geometries and Δ r G between the tested U VI ECP and basis sets. Therefore, to reduce computational cost, it is best to use either the MWB60 The Journal of Physical Chemistry A Article or the MDF60 basis set rather than the cc-pVDZ-PP basis set. The MDF60 ECP was used in subsequent calculations because it is fully relativistic.
Ligand Basis Set. Molecular geometries converge much faster than molecular energies when increasing basis set size, and geometries typically converge for double-ζ basis sets. 23 To test this hypothesis, the structures were optimized with Dunning's diffuse augmented correlation-consistent split valence basis sets, aug-cc-pVnZ, where n = D (double), T (triple), and Q (quadruple). The results presented in Table 4 suggest that there is virtually no change in the geometries of the uranyl complexes as the basis size increases. Typically, the UO and U−O water bond lengths and OUO angles vary by <0.1% between the different sized basis sets. The exception is U−O water bond distances of the quadruple-ζ basis set, which increases by 0.4% compared with the double-ζ basis set. The bond distances of U−O acetate vary by <0.2% between the three basis sets. Using the aug-cc-pVDZ basis set is adequate when wishing to predict the geometries of uranyl complexes.
To understand further the dependence of Δ r G on basis set size, the structures have been optimized with Dunning's diffuse augmented correlation-consistent split valence basis sets, aug-  The Journal of Physical Chemistry A Article cc-pVnZ, where n = D (double), T (triple), and Q (quadruple). The results are presented in Table 4 and Figure  3c. The diffuse functions were included due to the presence of acetate anions. As the number of basis functions increases, the energies of the complexes should increase, which may, in turn, affect Δ r G. Increasing basis set size from aug-cc-pVDZ to augcc-pVTZ changes Δ r G by −3.4, −4.8, and −4.9 kJ/mol to −122.5, −224.2, and −275.4 kJ/mol for the formation of 1:1, 1:2, and 1:3 uranyl acetate complexes, respectively. Changing from aug-cc-pVTZ to aug-cc-pVQZ changes Δ r G by a further −1.9, +0.3, and −3.2 kJ/mol to −124.4, −224.0, and −278.6 kJ/mol for the three reactions. In future calculations, the augcc-pVDZ basis set is used, as the geometries have converged at this basis set size, and the calculated Δ r G does not significantly change.
Functional. A comparison of Δ r G between the PW91, TPSS, B3LYP, and M06 functionals is found in Table 5 and is illustrated in Figure 3d. The calculated Δ r G typically improves with respect to experimental Δ r G in the order TPSS, PW91, B3LYP, and M06, with Δ r G predicted by TPSS between 18.1 and 26.6 kJ/mol more negative than M06. This indicates that the hybrid functionals are the most accurate at predicting Δ r G. The Δ r G predicted by the M06 functional differs from the B3LYP functional by +6.9, +5.0, and −12.6 kJ/mol for reactions 1 to 3, respectively, indicating that M06 is more accurate for reactions 1 and 2, whereas B3LYP is more accurate for reaction 3.
The average UO, U−O water , and U−O acetate bond lengths presented in Table 5 for the uranyl and uranyl acetate complexes vary by 3.06, 0.65, and 1.81%, respectively, for the four functionals. In general, the average UO bond length decreases in the order of PW91 (GGA), TPSS (meta-GGA), B3LYP (hybrid), and M06 (newer hybrid). The M06 functional underestimates experimental UO bond length by up to 2.31%, whereas PW91, TPSS, and B3LYP overestimate it by up to 2.13%. The average U−O water bond lengths increase in the order TPSS, PW91, M06, and B3LYP. The average U−O water bond lengths are overpredicted by the functionals, and this overprediction increases with the addition of acetate ligands. The difference between calculated and experimental U−O water can be up to 7.2% in the 1:2 U VI /L complex. The average U−O acetate bond lengths increase in the order TPSS, PW91, B3LYP, and M06, and tend to under predict experimental values for 1:1 and 1:2 U VI /L complex by up to 3% relative to the experimental data. 47 The exception is M06, which overpredicts the average U−O acetate bond length in the 1:2 complex by 0.6%. All of the functionals overpredict the average U−O acetate bond lengths for the 1:3 complex by up to 1.1% relative to the experimental values. The OUO bond angle is bent (167°) when employing the PW91 and TPSS functionals compared with the linear angle (178 to 179°) observed when the B3LYP and M06 functionals are employed.
Overall, a hybrid functional is selected for the calculations because exact exchange needs to be included in DFT methods for thermochemical calculations. The B3LYP functional is chosen over M06 to reduce computational cost.
Experimental Calibration. Finding isodesmic structures using mixed solvation is very challenging; the computational cost significantly increases, and there is no guarantee that the structures are at their global minima due to the flexible nature of the hydrogen-bonded water molecules. Using a fitting approach with [UO 2 (H 2 O) 5 ] 2+ as the starting structure, as proposed by ref 18, is currently the most practical method for scientists who are interested in exploring the relative stabilities of different U VI complexes. The fitting approach presented by ref 18 is modified here so that it can be employed in the future for ligands with many tens of atoms such as siderophores. This was done by employing Dunning's aug-cc-pVDZ basis set for the ligands, MDF60 ECP, and the corresponding basis set for U VI . This protocol keeps B3LYP as the functional and IEFPCM to model bulk solvation. The structures are optimized fully in the implicit solvation model rather than single point energy calculations employed in Vukovic et al. 18 Vukovic et al. 18 calculated the relative stability of 1:1 U VI complexes for monovalent and divalent ligands within 0.1:1.9 log β units with root-mean-square deviation (rmsd) < 1 log β. We extended this approach to organic molecules similar to the functional groups found in siderophores. The binding modes of the organic molecules tested have been selected based on previous experimental and computational work. 18,47,64−68 The ligands tested form bidentate complexes, with the exception of glycolate, which forms a monodentate complex. The reaction used for the stability constant is given in eq 2. This is with the exception of L-glycine, which uses eq 12, in line with experimental derivation. 57 Experimental and calculated stability constants are presented in Table 6. 48  from the data plotted in Figure 7 and is used to fit the calculated log β of 1: The log β values of U VI complex formation with formic and nicotinic acid were subsequently calculated to test the accuracy of the fitting curve for calculating the log β of ligands beyond the training set. These two ligands have been selected for the two different reactions shown in eq 2 for formic acid and eq 12 for nicotinic acid. Formic acid forms a monodentate complex 73 with log β 2.70, 74 corrected to zero ionic strength using SIT theory. The calculated log β for the complex is 15.81, and applying the fitting equations reduces the log β to 2.18. This is within 0.52 log β units of the experimental data and is less than the standard error (1.64 log β units) and rmsd (1.38 log β units) of the equation. To our knowledge, no structural data are available for the U VI nicotinic acid complex. Protonation data indicate that the amine group is protonated and the ligand binds to U VI via the carboxylate end of the molecule. 75 This could be in a bidentate fashion, as shown in eq 12, or a monodentate fashion, with one less water molecule released. The experimental log β 75 is 8.83 after correction to zero ionic strength using the Davies equation. The calculated stability constants for the monodentate and bidentate complexes are 37.39 and 38.12, respectively, and applying the fitting equation reduces the log β to 10.59 and 10.88. The difference in log β to the experimental data for the monodentate complex is 1.76 log β units, and that for the bidentate complex is 2.04 log β units. These are higher than both the rmsd (1.38 log β units) and standard error (1.64 log β units) of the fitting equation, however the experimental error has not been considered. The experimental error is ±0.17 log β units calculated from the average of four independent measurements. 75 The error increases to ±0.41 log β units after correcting to zero ionic strength. Applying the ionic strength corrected error can decrease the difference to the fitted data to 1.37 and 1.66 log β units for the monodentate and bidentate complex respectively, close to the rmsd (1.38 log β units) and standard error (1.64 log β units).
These test ligands indicate that the fitting curve significantly improves the estimated stability constant for ligands outside the ligand training set.
Sequential Complexation. To our knowledge, using a fitting approach for multiple ligands coordinated around U VI has not been explored in detail to date. The fitting for 1:2 U VI ligand complexes was studied recently; 31 however, this was conducted with the SMD solvation model rather than the IEFPCM solvation model that we employ here. To examine if a fitting can be used for 1:1 to 1:3 U VI ligand complexes on one fitting curve, the stability of U VI glycolate, acetate, and oxalate ligands was explored. All three glycolate ligands bind in a monodentate style via the carboxylate end (at pH <3.5), whereas all three acetate ligands bind in a bidentate fashion. 47,68 Oxalate binds in a bidentate fashion in the 1:1 and 1:2 U VI /L complex, while the third oxalate binds in a monodentate fashion. 64,65 Figure 8a shows that the experimental log β loosely increases exponentially with respect to calculated log β. The correlation is 0.889 for the nine complexes investigated during this study. Further inspection highlights the complexity of such a fitting.   Figure 8b, the data have been split by ligand charge, that is, monoanionic charged ligands, glycolate, and acetate (green) and dianionic charged oxalate (black). Both of these curves are linearly correlated (R 2 = 0.9131 for −1 charge glycolate and acetate and R 2 = 0.9918 for −2 charged oxalate). The trends and the level of overestimation of the calculated log β are similar between glycolate and acetate. Oxalate underpredicts the calculated stability compared with glycolate and acetate and follows a steeper trend. Vukovic et al. split the fitting of 1:1 U VI /L complexes via ligand charge and improved accuracy in the fitted log β values. The results in Figure 8b support the evidence that splitting via ligand charge may be important, particularly when exploring multiple bound ligands on one fitting curve. Figure 8c splits the complexes via the number of ligands complexed with U VI , where blue, green, and black represent one, two, and three ligands complexed around U VI , respectively. The results highlight that 1:2 and 1:3 complexes follow the same trend with similar correlations (R 2 = 0.9092 and 0.9056 for 1:2 and 1:3 complexes, respectively). This means that 1:2 and 1:3 complexes can be fitted on the same curve regardless of ligand charge, indicating that if 1:2 or 1:3 experimental data is missing, then it can be inferred from this plot. The trend for 1:1 complexes has a similar correlation (R 2 = 0.922); however, the calculated stability constants are underpredicted (shifted to the left) compared with 1:2 and 1:3 complexes. This indicates that the DFT calculations suggest it is more difficult to complex the first ligand, but once the first ligand is in, it becomes easier for the next two ligands to complex with U VI .

■ CONCLUSIONS
The first objective of the study presented here was to modify the protocol for estimating the relative stability of U VI complexes proposed by ref 18. This adjusted protocol allows us to estimate the relative stability of functional groups that comprise the reactive components in siderophores while also providing a protocol whose computational cost will not prohibit future exploration of siderophores with tens of atoms. This is important because siderophores complex with U VI , significantly affecting U VI aqueous chemistry and mobility in natural systems. The second objective was to assess whether the difficulties associated with mixed solvation methods are outweighed by improvements in Δ r G. The third objective was to explore whether the fitting approach can be expanded to multiple bound ligands, focusing on carboxylate type ligands. The suggested protocol is the B3LYP functional, aug-cc-pVDZ ligand basis set, MDF60 ECP and basis set for U VI , and the IEFPCM solvation model. The complexes were fully optimized in the solvation model. Mixed solvation models based on previous research 22 were tested for U VI acetate complexes. We demonstrate that applying the approach employed in ref 22 significantly improves calculated Δ r G; however, a fitting approach would still need to be employed to determine accurate log β values. To save on computational cost and to avert difficulties in keeping reactions isodesmic, all water molecules outside of the uranyl complex [UO 2 (H 2 O) 5 ] 2+ should be modeled implicitly with the IEFPCM solvation model. Using the protocol outlined above, in combination with a linear regression allows the user to obtained fitted log β for 1:1 U VI /L complexes with rmsd of 1.38 log β units and a standard error of 1.64 log β units. This enables us to develop a relative stability series for organic ligands and to aid future experiments by providing molecular level insight into macroscale experiments.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.8b06863.
Details of correcting ionic strength using specific ion theory and the Davies equation along with the full Gaussian 16 rev. A.03 reference, and structure energies (PDF) Structure coordinates (PDF)