First-principles molecular dynamics simulation of the Ca 2 UO 2 ( CO 3 ) 3 complex in water †

Recent experiments have shown that the neutral Ca2UO2(CO3)3 complex is the dominant species of uranium in many uranyl-containing streams. However, the structure and solvation of such a species in water has not been investigated from first principles. Herein we present a first principles molecular dynamics perspective of the Ca2UO2(CO3)3 complex in water based on density functional theory and Born–Oppenheimer approximation. We find that the Ca2UO2(CO3)3 complex is very stable in our simulation timeframe for three different concentrations considered and that the key distances from our simulation are in good agreement with the experimental data from extended X-ray absorption fine structure (EXAFS) spectroscopy. More important, we find that the two Ca ions bind differently in the complex, as a result of the hydrogen-bonding network around the whole complex. This finding invites confirmation from time-resolved EXAFS and has implications in understanding the dissociative equilibrium of the Ca2UO2(CO3)3 complex in water.


Introduction
The concentration of uranium in seawater is very minute at about 3.3 ppb, but the large mass of seawater contains about 4.5 billion metric tons of uranium. 1 Mining uranium from seawater can provide an ideal sustainable alternative to current uranium mining methods (open pit, underground, and leaching) that create much environmental concern. In addition, extraction from seawater can provide a price ceiling for uranium that is of strategic importance in evaluating the economics of long-term supply of uranium for nuclear energy. 2,3 Much research has been focused on the development of selective methods for mining uranium in seawater. Currently, the most viable method is utilizing an amidoxime-functionalized polymer sorbent known as poly(acrylamidoxime) fibers. [1][2][3][4][5][6] However, to make the uranium extraction from seawater economically viable, sorbent performance in terms of uranium uptake and uranium/vanadium selectivity needs to be further improved. 1,7 To design a better ligand and sorbent for uranium extraction, it is necessary to understand the fundamentals about uranium speciation in seawater.
In aqueous solution, uranium exists as the stable oxocation with an oxidation state of U(VI), called uranyl -UO 2 2+ . Early work has focused on prying into the prominent equilibrium species bound to the uranyl complex in seawater 8 and suggested the anionic [UO 2 (CO 3 ) 3 ] 4− complex to be the dominant species in seawater. However, over the past two decades experimental data has shifted the consensus to cationbalanced complexes. [9][10][11][12][13] Concentrations of magnesium (Mg 2+ ) and calcium (Ca 2+ ) in seawater are overwhelmingly larger than the concentration of U(VI), so the ternary Ca-UO 2 -CO 3 or Mg-UO 2 -CO 3 exists predominately in seawater. The complexation of Ca 2+ with [UO 2 (CO 3 ) 3 ] 4− has been validated experimentally by Bernhard et al. 10 and Kelly et al. 11 with extended X-ray absorption fine structure (EXAFS) spectroscopy. Most recently, Rao et al. 13 examined the thermodynamics of uranium in seawater and the complexation between Ca/Mg and [UO 2 (CO 3 ) 3 ] 4− . They concluded that in seawater pH (8.2) Ca 2 UO 2 (CO 3 ) 3 accounts for 58% of the total uranium in the solution while CaUO 2 (CO 3 ) 3 2− and MgUO 2 (CO 3 ) 3 2− account for 18% each and [UO 2 (CO 3 ) 3 ] 4− accounts for only 6%. 13 In addition, Rao et al. studied the binding of U(VI) with various types of ligands and the subsequent leaching process. [14][15][16][17][18][19] In addition, the stability constant for the speciation of calcium is larger than magnesium. 13 On the theoretical and computational side, work has been done on the binding of UO 2 2+ with ligands using a cluster model 18,20- The present work seeks to describe the structure and solvation of Ca 2 UO 2 (CO 3 ) 3 in water using first principles MD based on density functional theory (DFT-MD for short) for the first time. Our goal is to provide a fundamental baseline understanding of the structure and solvation of Ca 2 UO 2 (CO 3 ) 3 in water in terms of Ca-UO 2 (CO 3 ) 3 and Ca-water interactions. The other goal is to compare DFT-MD simulations with the previous EXAFS data and classical MD simulations.

Computational methods
First-principles molecular dynamics simulations based on density functional theory (DFT-MD) and Born-Oppenheimer approximation were carried out using the Vienna ab initio simulation package with plane wave basis and periodic boundary conditions. 28,29 The Kohn-Sham equations are solved with the all-electron projected augmented wave (PAW) method. 30,31 We have chosen the Perdew-Burke-Ernzerhof (PBE) functional of the generalized-gradient approximation (GGA) for electron exchange and correlation. 32 PBE is one of the most popular GGA functionals, providing a balanced description for diverse molecules and materials, instead of being designed for a special class of molecules or interactions. In the case of liquid water, it has been shown that the PBE functional can describe well the peak positions in the radial distribution functions of g OO and g OH for the liquid structure of water but it overestimates the peak heights, in comparison with the experiment, leading to over-structuring. 33,34 Using hybrid functionals together with van der Waals interactions can soften the water structure, giving a better agreement with the experiment. But hybrid functionals are usually about two orders of magnitude more expensive than a pure DFT method such as GGA-PBE. We think that PBE is a reasonable choice in balancing accuracy and efficiency. In the ESI, † we provide an orbital-resolved local density of states for the U atom (Fig. S1 †) from a snapshot of the Ca 2 UO 2 (CO 3 ) 3 complex in water to show that the bonding picture from the PBE functional is consistent with previous theoretical studies. 35,36 The MD calculations were carried out at 298 K in a canonical NVT ensemble for a periodic cubic box that contains one Ca 2 UO 2 (CO 3 ) 3 complex in a fixed number of water molecules. Three concentrations were examined: 0.53 M, 0.42 M, and 0.36 M, corresponding to one Ca 2 UO 2 (CO 3 ) 3 complex in a periodic box containing 100, 125, and 150 water molecules, respectively; the corresponding simulation box sizes and densities are also compared in Table 1. Here we note that since there is only one complex in the simulation box for the three concentrations, this approach cannot probe the correlations between complexes but serves more to test the potential presence of size artifacts. We determined the densities from constant-pressure classical MD simulations using force-field parameters from a previous study. 26 The temperature was kept constant via a Nose-Hoover thermostat. A Verlet algorithm was used to integrate Newton's equation of motion with a time step of 1 femtosecond. After equilibration at 298 K for 15 ps, another 15 ps of production run was followed. Graphical visualization and analysis of the liquid structure packing of the uranium complex was examined with VMD. 37

Interaction between calcium and carbonate
The most important structural feature of the Ca 2 UO 2 (CO 3 ) 3 species is the binding between the two Ca ions and the [UO 2 (CO 3 ) 3 ] 4− ion. This interaction is mediated by the carbonate groups. As shown in Fig. 1, the three carbonate groups bind to the uranyl group on the equatorial plane in a bidentate mode; this structural model has been established from the crystal structure of the naturally occurring mineral Liebigite [Ca 2 UO 2 )(CO 3 ) 3 ·11H 2 O], 38 fitting of the EXAFS data, 10,11 and quantum mechanical modeling. 20 There are three different oxygens in the uranium complex: the two axial oxygens (O ax ) triple-bonded to U in the uranyl structure, 35 six equatorial carbonate oxygens (O eq ) that are bonded to U, and three distal carbonate oxygens (O dis ) not directly interacting with U. The two Ca 2+ ions bind to the carbonate groups on the same plane; each Ca 2+ ion binds to two equatorial oxygen atoms from two neighboring carbonate groups. In our DFT-MD simulations, interaction of the Ca ions to the [UO 2 (CO 3 ) 3 ] 4− complex was monitored by the four Ca-O eq distances (dashed lines in Fig. 1): Ca1-O1, Ca1-O2, Ca2-O3, and Ca2-O4.
We placed an initial structure of the uranium complex as shown in Fig. 1 into a periodic water box at a concentration of 0.53 M. After equilibration at 298 K, a production run of 15 ps was used for statistical analysis. Fig. 2 shows the four Ca-O eq distances during the 15 ps production trajectory. One can see that Ca1-O1 and Ca1-O2 distances exhibit fluctuations around 2.45-2.50 Å (Fig. 2a), while Ca2-O3 and Ca2-O4 around 2.35-2.40 Å (Fig. 2b). So in our simulation timeframe, the Ca 2 UO 2 (CO 3 ) 3 complex is very stable and maintains steady Ca-O eq distances about 2.45 Å with a standard deviation of  Fig. 2a). Here we note that initially, we placed the Ca 2 UO 2 (CO 3 ) 3 complex randomly inside a water box. To test the robustness of the asymmetric structure, we tried several different initial configurations of water solvation around the Ca 2 UO 2 (CO 3 ) 3 complex and found that they always equilibrated to the asymmetric configuration after about 5 ps.
To further examine the difference between the two Ca ions, we plot the radial distribution function (RDF) of carbonate O eq atoms around each of the two Ca ions in Fig. 3. One can see that the stronger binding Ca 2+ has a narrower and higher O eq distribution ( Fig. 3b), while the weaker Ca 2+ has a broader and lower O eq distribution (Fig. 3a). In addition, there is a slight difference between the two O eq atoms binding to each Ca 2+ . For Ca1, Ca1-O1 is slightly shorter than Ca1-O2; for Ca2, Ca2-O4 is slightly shorter than Ca2-O3.   To confirm the stability of the Ca 2 UO 2 (CO 3 ) 3 complex and the asymmetry of the two Ca 2+ ions, we further simulated two lower concentrations (0.42 M and 0.36 M) and arrived at the same conclusions. The Ca 2 UO 2 (CO 3 ) 3 complex in the two lower concentrations is also stable in our simulation time frame, as shown by the steady maintaining of the binding of the two Ca 2+ ions with the [UO 2 (CO 3 ) 3 ] 4− complex. More interestingly, we found that the asymmetry between the two Ca 2+ ions also persists in the two lower concentrations, indicating that this is likely an intrinsic feature of the Ca 2 UO 2 (CO 3 ) 3 complex in water. Fig. 4 displays the four Ca-O eq distances as a function of the U concentration. Both the asymmetry between the two Ca 2+ ions and the small difference between the two O eq atoms for each Ca 2+ ion are evident.

Interaction between calcium and water
The interaction between calcium and carbonate in the Ca 2 UO 2 (CO 3 ) 3 complex is the most important information that we obtained from our DFT-MD simulations. The asymmetry between the two Ca 2+ ions must be closely related to the water molecules around the Ca 2 UO 2 (CO 3 ) 3 complex. We now analyze the interaction between the two Ca 2+ ions and the water molecules. Fig. 5 shows radial distribution functions (RDFs) of oxygen atoms from the water molecules around the two Ca 2+ ions both separately and together. One can see that the solvation shell around Ca1 has an average Ca-O water distance of 2.45 Å (with a standard deviation of 0.12 Å) and the integrated RDF (with a cutoff at 3.0 Å) gives the coordination number of five; in other words, there are five molecules around Ca1 in addition to the two O eq atoms from two carbonate groups. On the other hand, Ca2 has four water molecules in the solvation shell with an average Ca-O water distance of 2.35 Å (with a standard deviation of 0.09 Å). So together, the average coordination number of the two Ca 2+ ions is 4.5 in terms of water mole-cules. We further examined the RDF of water oxygens around Ca 2+ ions for the two lower concentrations and found the same trend of five water molecules around Ca1 and four water molecules around Ca2. In comparison, previous classical MD simulations predicted that both calcium ions have five water molecules in the first hydration shell, 26 similar to the case of Ca1 in our simulation.

Solvation environments of the two calcium ions
From the above discussion of the Ca 2 UO 2 (CO 3 ) 3 complex in water, we can clearly see that the difference between the two Ca ions is reflected in both the Ca-carbonate and the Cawater interactions. The two interactions are in fact correlated: Ca1 has weaker binding with the [UO 2 (CO 3 ) 3 ] 4− complex, five molecules in the solvation shell, and a total of seven coordination bonds; Ca2 has stronger binding with the [UO 2 (CO 3 ) 3 ] 4− complex, four molecules in the solvation shell, and a total of six coordination bonds. Ca2 has a tighter solvation shell, so both average Ca2-O water and Ca2-O carbonate distances are shorter than Ca1-O water and Ca1-O carbonate distances, respectively.
What causes the asymmetry of binding and solvation between the two Ca ions in the Ca 2 UO 2 (CO 3 ) 3 complex? To answer this question, we analyzed the solvation environment of the complex from the views of Ca-carbonate, Ca-water, and carbonate-water interactions together, as shown in Fig. 6. One can see that there are three water molecules in the equatorial plane coordinating to Ca1, instead of two in the case of Ca2. The reason why Ca1 can have one more water molecule in the equatorial plane is that two of the three water molecules (water1 and water2 in Fig. 6) are interacting with both Ca1 and the carbonates. From Fig. 6, one can see that both water1 and water2 form hydrogen bonding (hb1 and hb2) with the two distal oxygen atoms of the two carbonate groups. These two  hydrogen bonds pull water1 and water2 closer to the Ca 2 UO 2 (CO 3 ) 3 complex, thereby leaving space for a third water molecule to enter the equatorial plane. In other words, it is the hydrogen-bonding network around the Ca 2 UO 2 (CO 3 ) 3 complex that leads to the difference in solvation and binding between the two Ca ions. We further examined the two lower concentrations and found the same solvation environment around the Ca 2 UO 2 (CO 3 ) 3 complex that confirmed the role of the hydrogen-bonding network in differentiating the two Ca ions in the Ca 2 UO 2 (CO 3 ) 3 complex.
To further explain the asymmetry between the two Ca ions, we show a schematic drawing (Fig. 7) of the equatorial plane around U. One can see that Ca1 is coordinated by both water1 (W1) and water2 (W2), while W1 is hydrogen bonded to O5 of carbonate1 (C1) and W2 is hydrogen bonded to O6 of carbon-ate2 (C2). As a result, the hydrogen bonding pulls the two carbonate groups closer (indicated by the two black arrows), so Ca1 is "squeezed" a little further away from O1 and O2. On the other hand, the O3-U-O4 angle becomes wider (indicated by the red double arrow), thereby allowing Ca2 to come closer to O3 and O4. Another way to think about this is via carbonate1. If Ca1 and Ca2 were symmetric in binding, the hydrogen bonding around carbonate1 (C1) would be symmetric. But as shown in Fig. 7, the hydrogen bonding around carbonate1 is asymmetric that eventually leads to the asymmetry in binding between Ca1 and Ca2.

Solvation environments of the whole complex
The discussion above shows the importance of the hydrogenbonding network in dictating the complex geometry. To further analyze this network, we examined the first solvation shell of the whole complex, namely, the water molecules in direct interaction with the complex. Since we have analyzed the water solvation around the two Ca ions, here we focus our discussion on the carbonate and uranyl oxygens. One can see from Fig. 8 that the top uranyl oxygen (O A1 ) has two water Fig. 6 A snapshot of the Ca 2 UO 2 (CO 3 ) 3 complex in water at 0.53 M, showing only the water molecules directly interacting with the two Ca ions; hb1 and hb2 denote hydrogen bonding between the two water molecules (water1 and water2) and the two distal oxygen atoms of the two carbonate groups around Ca1; atom-atom distances are labeled in Å.  molecules hydrogen-bonded to it, while the bottom uranyl oxygen (O A2 ) has one. Moreover, one can see strong solvation of the carbonate distal oxygens by water: O D1 is hydrogenbonded by three water molecules, O D2 by four, O D3 by two. In addition, the two carbonate equatorial oxygens not interacting with the Ca ions are also solvated by water. Together with the water molecules around the two Ca ions, we found that there are 21 molecules in the first solvation shell. This large solvation shell indicates the necessity of using the explicit solvation model to address the structure, thermodynamics, and chemistry of the aqueous Ca 2 UO 2 (CO 3 ) 3 complex.

Ca-U distances
Besides speciation studies based on thermodynamics, 10,12,13 the most direct characterization of the Ca 2 UO 2 (CO 3 ) 3 complex in water has been EXAFS studies of the coordination shells around the central U atom. 10,11 Since the Ca-U distance is a key piece of information available from fitting the EXAFS spectra, we examined in detail the Ca-U distances for the Ca 2 UO 2 (CO 3 ) 3 complex in water. Fig. 9 shows the RDF of Ca ions around the U atom at three different concentrations. One can see that the asymmetry between the two Ca ions is also reflected in the Ca-U distances: the stronger-binding Ca2 is about 4.05 Å away from U and has a narrower distribution of the Ca2-U distance (standard deviation: 0.10 Å), while the weaker-binding Ca1 is about 4.15 Å away from U and has a broader distribution of the Ca1-U distance (standard deviation: 0.12 Å). In addition, the three concentrations show very consistent distributions of Ca-U distances (Fig. 9).

Comparison with the literature
To our knowledge, the present study is the first DFT-MD simulation of the Ca 2 UO 2 (CO 3 ) 3 complex in water. It would be very informative to compare the present DFT-MD results with previous experiments and molecular-mechanical MD simulations (MM-MD) based on empirical force fields. For comparison with the experiment, we focus mainly on the liquid-phase EXAFS analysis on the structure of the Ca 2 UO 2 (CO 3 ) 3 complex from Kelly et al. 11 and Bernhard et al. 10 In the model fitting of the EXAFS spectra, they assumed that the two U-Ca distances are the same. To directly compare with their data, we therefore computed the total RDF of Ca ions around U and obtained an average Ca-U distance of about 4.07 Å at the peak of the RDF. Table 2

Implications of the present findings
As we discussed above, a key finding from the present DFT-MD simulation is the asymmetry between the two Ca ions in the Ca 2 UO 2 (CO 3 ) 3 complex. A key issue here is whether and how often the two Ca ions can switch their bonding environments,  namely, from Ca1-weak binding/Ca2-strong binding to Ca1strong binding/Ca2-weak binding. Such switching will be closely related to water exchange in the first solvation shell of the Ca ions. In our limited simulation timeframe (∼50 ps), we did not observe such switching. This implies that our bruteforce DFT-MD is unlikely to address this issue due to its limited accessible timescale that is too short in comparison with the timescale of such switching. We are currently pursuing two lines of research to address this issue that will be published in the near future: (a) DFT-MD coupled with metadynamics to estimate the free-energy profile of such switching; (b) classical MD based on force fields to increase the timescale to about ∼100 ns. Despite the limited timescale of the present DFT-MD simulation, our finding of the asymmetry between the two Ca ions in the aqueous Ca 2 UO 2 (CO 3 ) 3 complex in the timescale of 10 to 100 ps may be confirmed by time-resolved EXAFS that can measure the variation in bond length in ps timescale. 40 Here we suggest an experiment to use time-resolved EXAFS to measure the Ca-U distances of the aqueous Ca 2 UO 2 (CO 3 ) 3 complex at ps snapshots. Another implication from our finding concerns dissociation of Ca 2 UO 2 (CO 3 ) 3 to CaUO 2 (CO 3 ) 3 2− . Rao et al. found that in seawater conditions, Ca 2 UO 2 (CO 3 ) 3 and CaUO 2 (CO 3 ) 3 2− account for 58% and 18% of total U(VI), respectively. 13 In other words, Ca 2 UO 2 (CO 3 ) 3 is in equilibrium with CaUO 2 (CO 3 ) 3 2− and free Ca 2+ in seawater.
Our finding suggests that Ca1 is much more likely to break away from Ca 2 UO 2 (CO 3 ) 3 than Ca2, to form CaUO 2 (CO 3 ) 3 2− . This information will be useful for studies of the mechanism of Ca 2 UO 2 (CO 3 ) 3 dissociation with or without an attacking ligand. We plan to also use DFT-MD coupled with metadynamics to examine the free-energy profile of the dissociation mechanism. We suspect that there may exist some intermediate states of the Ca 2 UO 2 (CO 3 ) 3 complex before it becomes [CaUO 2 (CO 3 ) 3 ] 2− and free Ca 2+ . For example, one likely configuration can have one Ca ion coordinating to one equatorial and one distal oxygen from the same carbonate group, while the other Ca ion coordinates "normally" to two equatorial oxygens of two different carbonate groups.

Conclusion
We have simulated the neutral Ca 2 UO 2 (CO 3 ) 3 complex in water using first principles molecular dynamics based on density functional theory (DFT-MD). Three concentrations (0.53, 0.42, and 0.36 M) feasible to DFT-MD simulations were examined. In the accessible timescale (∼30 ps), we found that the structure of the Ca 2 UO 2 (CO 3 ) 3 complex is very stable where the two Ca ions bind to the carbonate groups on the same equatorial plane. We found that one Ca ion binds to the center UO 2 (CO 3 ) 3 4− anion more stronger than the other Ca ion. This asymmetry of binding between the two Ca ions is reflected in several aspects: the stronger binding Ca has shorter Ca-O carbonate bonds, shorter Ca-U distance, and four coordinating water molecules, while the weaker binding Ca has longer Ca-O carbonate bonds, longer Ca-U distance, and five coordinating water molecules. This finding suggests that using timeresolved EXAFS spectra may confirm the asymmetry in binding of the two Ca ions in the aqueous Ca 2 UO 2 (CO 3 ) 3 complex, since our DFT-MD simulation shows in general good agreement in terms of key distances with the EXAFS experiments.