Ab initio molecular dynamics studies of hydroxide coordination of alkaline earth metals and uranyl

Abstract

Ab initio molecular dynamics (AIMD) simulations of the Mg²⁺, Ca²⁺, Sr²⁺ and UO₂²⁺ ions in either a pure aqueous environment or an environment containing two hydroxide ions have been carried out at the density functional level of theory, employing the generalised gradient approximation via the PBE exchange–correlation functional. Calculated mean M–O bond lengths in the first solvation shell of the aquo systems compared very well to existing experimental and computational literature, with bond lengths well within values measured previously and coordination numbers in line with previously calculated values. When applied to systems containing additional hydroxide ions, the methodology revealed increased bond lengths in all systems. Proton transfer events (PTEs) were recorded and were found to be most prevalent in the strontium hydroxide systems, likely due to the low charge density of the ion and the consequent lack of hydroxide coordination. For all alkaline earths, intrashell PTEs which occurred outside of the first solvation shell were most prevalent. Only three PTEs were identified in the entire simulation data of the uranium dihydroxide system, indicating the clear impact of the increased charge density of the hexavalent uranium ion on the strength of metal–oxygen bonds in aqueous solution. Broadly, systems containing more charge dense ions were found to exhibit fewer PTEs than those containing ions of lower charge density.

---

1. Introduction

The coordination environment of alkaline earth metals in aqueous environments has been the subject of both experimental and theoretical investigations for decades.^1–7 Both magnesium and calcium are of significant biological importance^8,9 and are key components in natural groundwater,¹⁰ whereas strontium is often associated with nuclear waste. Radioactive ⁹⁰Sr is produced via nuclear fission of ²³⁸U and is a potential health hazard due to chemical similarity to calcium and a consequent analogous biological absorption profile. Immobilisation of ⁹⁰Sr present in high and intermediate level nuclear waste (HLW and ILW) is an ongoing concern of the nuclear power industry.

The first generation of British nuclear reactors used solid uranium fuel rods clad with a magnesium alloy, known as Magnox. Prior to the reprocessing or final disposal of the spent nuclear fuel the fuel is stored under water in ponds for cooling and radiation shielding resulting in the formation of uranyl (UO₂²⁺). To inhibit corrosion of the Magnox, the ponds are typically held at high pH but over their operational lifetime, since the 1950s, this has not always been the case, resulting in some fuel being exposed to a pH range between 10 and 12. The pond conditions combined with the decades of fuel storage has caused a particulate sludge of Magnox corrosion product (brucite and hydromagnesite) to form at the base of the legacy ponds which readily interacts with radionuclides and other ions present in the liquor.^11,12 Many of these storage ponds are reaching the end of their viable lifetime and thus require emptying and decommissioning. The primary radionuclides in the storage ponds are ²³⁸U and two fission products, ¹³⁷Cs (t_1/2 = 31.2 y) and the aforementioned ⁹⁰Sr (t_1/2 = 28.8 y), existing primarily as carbonate and hydrated hydroxide complexes. Detailed understanding of the solubility, and ability of the group 1 and 2 metals, in particular Cs, Sr and Mg, to interact with the minerals present is of particular importance to the decommissioning process in terms of sludge removal and the removal of radioactivity from the effluent that is generated.

To fully investigate the interaction between these ions and a mineral surface it is essential to have a detailed understanding of the microsolvation of the ions in both the absence and presence of hydroxide. While the former has been well documented in literature there is severely limited research on the later with literature focusing on gas phase^13–15 investigations, while the dynamics of hydrated hydroxide complexes have received limited attention.^16–19 This contribution aims to develop understanding of ion interactions with hydroxide in aqueous environments.

Recently, molecular dynamics simulations (MD) have gained popularity as method for investigating such dynamical processes.^20–24Ab Initio Molecular Dynamics (AIMD), in the form of Born–Oppenheimer Molecular Dynamics (BOMD) or Car–Parrinello Molecular Dynamics (CPMD)^25–27 uses quantum chemical techniques to determine interatomic interactions. One of the main advantages of AIMD is that it can be used to study chemical bond-breaking and formation,²⁸ as such it has been used to study proton transport through water,²⁹ hydrogen bonding³⁰ and the microsolvation of ions.³¹ It has also been employed to predict experimental observables linked to electronic structure, including NMR,²⁵ IR³² and Raman spectra.³³

In this contribution we employ BOMD to investigate the structure and dynamics of hydrated alkaline earth metals and uranyl, UO₂²⁺, in the presence of hydroxide. Whilst of the alkaline earths considered here, Mg, Ca and Sr, only the latter is a direct product of the nuclear fission process, the others are present in the aforementioned storage ponds and in the ground water around nuclear sites. The study of Mg and Ca also serves to highlight the variation in the nature of the ion–hydroxide interaction through the alkaline earth group, whilst the inclusion of uranyl allows consideration of a more complex ion with significantly different electrostatic properties. Consequently, improving our microscopic understanding of their solvation properties is of value and, furthermore, allows for the investigation of the applicability of BOMD in discriminating between properties of closely related species. The coordination environment of aqueous Mg²⁺, Ca²⁺, Sr²⁺ and UO₂²⁺ complexes have been well studied with both computational and experimental methods and inspection of the existing literature on alkaline earth hydration show significant variation.

The literature on the hydration structure of Mg²⁺ is relatively unambiguous. X-Ray Diffraction (XRD),^2,34 Neutron Diffraction (ND),³⁵ Raman spectroscopy,^36,37 Density Functional Theory (DFT),^38,39 MD^34,36,40 and AIMD^41–44 investigations all report a coordination number (CN) of 6 with the first shell Mg–O distance is found experimentally to be between 2.10 and 2.12 Å.^{35–37,45–48} One gas phase DFT study investigating the successive binding energies of water to Mg²⁺ by Pavlov et al.⁴⁹ found stable structures up to a CN of 7. A CN of 6.8 was reported in an XRD study by Albright⁴⁸ but the author warned that the poor resolution of the peaks was the likely cause of the inconsistent result.

The hydration structure of Ca²⁺ is less well-defined, with experimental investigations including extended X-ray absorption fine structure (EXAFS),⁵⁰ XRD^5,48,51,52 and ND^35,50,53 reporting CNs between 6^49,52 and 7.9,⁴⁸ and Ca–O distances of 2.4⁵³ to 2.45 Å,⁵¹ while combined experimental computational investigations^54–56 reporting CNs of 6⁵⁵ to 9.5,⁵⁴ and Ca–O distances of 2.45⁵⁶ to 2.51 Å.⁵⁴ Katz et al.⁵⁷ used ab initio molecular orbital calculations to examine structures of Ca²⁺ complexes and similar energies were found for CNs 6, 7, and 8 while all other CNs had considerably larger energies. MD calculations^58–60 report similar variation with CNs ranging from 5.95⁵⁹ to 8⁶⁰ and Ca–O lengths of 2.35⁵⁹ to 2.45 Å.⁶¹ The quantum chemical statistical mechanical (QMSTAT) study by Tofteberg et al.⁵⁸ reported a CN of 6.9 with a Ca–O length of 2.5 Å. Various CPMD^41,42,62,63 studies of the hydration structure of Ca²⁺ calculated CN in the range 5.9⁴¹ to 7.2⁶² and Ca–O distances of 2.36⁴¹ to 2.45 Å.⁶³ A BOMD study by Mehandzhiyski et al.⁶⁴ found a CN of 7 with a Ca–O distance of 2.45 Å, in agreement with earlier CPMD studies.

The hydration structure of Sr²⁺ is similarly varied in the literature. EXAFS studies identify structures with CNs from 6⁶⁵ to 10.3⁶⁶ and first shell Sr–O distances of 2.57⁶⁵ to 2.63 Å,⁶⁶ X-ray Absorption Fine Structure (XAFS) report a smaller CN range of 6.22⁶⁷ to 7.3⁶⁸ and first shell Sr–O distances of 2.6⁶⁷ to 2.62 Å,⁶⁸ while XRD studies report a CN of around 8^{37,48,69–71} and first shell distances of 2.6⁴⁸ to 2.64 Å.³⁷ A ND study by Neilson et al.⁷² of Sr(ClO₄)₂ reported a much higher CN of 15 with Sr–O first shell distance of 2.65 Å, however this data was later re-examined in conjunction with an Anomalous X-ray Diffraction (AXD) study to find a lower CN of 9.⁷¹

Various computational methods have been used to evaluate the first shell solvation structure of Sr²⁺ including DFT,⁷³ Quantum Mechanical/Molecular Mechanics (QMMM),⁷⁴ QMSTAT⁵⁸ as well as MD and AIMD.^41,64,75 A Sr CN of around 8^{40,56,64,76–79} is typically identified within a range of CNs between 6.7⁴¹ and 9.8⁴⁷ and Sr–O distances of 2.58⁴⁰ to 2.69.⁵⁸ A recent paper by D’Angelo et al.⁷⁹ combined experimental and computational techniques to investigate the coordination shell of Sr²⁺ using X-ray absorption near-edge spectroscopy (XANES) of [Sr(H₂O)₈](OH)₂, MD and CPMD. The authors concluded that a CN of 8 with a first shell Sr–O distance of 2.6 Å was the most accurate description, reporting that the 2.72 Å bond length calculated with CPMD was inaccurate compared to their classical simulations and other AIMD literature.^41,75 A BOMD study by Mehandzhiyski et al.⁶⁴ included a DFT-D2 dispersion correction to better model the interactions of the water molecules, something absent from earlier investigations and found a CN of 7.6 with 2.60 Å Sr–O first shell distance.⁶⁴

EXAFS,^80,81 High Energy X-ray Scattering (HEXS),⁸² High Field NMR,⁸³ XAFS,^84,85 X-Ray Scattering⁸⁶ and XANES⁸⁷ have been used to probe the solvation structure of aquo complexes of UO₂²⁺. They report equatorial coordination numbers (CN) of 4.5⁸⁰ to 5.3,⁸⁴ axial uranyl (U–O_yl) bond distances of 1.70⁸⁸ to 1.77 Å⁸⁶ and equatorial (U–O) bond distances of 2.41^80,84,89 to 2.45 Å.⁹⁰ Both the HEXS study of Soderholm et al.⁸² and the X-ray scattering investigation by Neuefeind et al.⁸⁶ found a dynamic equilibrium between a four and five coordinated uranyl, [UO₂(H₂O)₄]³⁺ and [UO₂(H₂O)₅]²⁺ respectively, in which the five coordinated species was favoured.

Various DFT investigations of the aquo solvation structure of uranyl found a equatorial uranyl CN of 5,^91–98 with equatorial U–O bond lengths of 2.4^93,95 to 2.53 Å.^94,97 In general both AIMD^99–101 and MD^87,102–105 simulations of the uranyl solvation environment indicated a coordination number of 5 and the U–O bond distance is identified as 2.36¹⁰² to 2.48 Å.^104,105 The MD study by Rodríguez-Jeangros et al.¹⁰⁶ identified an average CN of 4.39 as uranyl is equatorially coordinated by either 4 or 5 waters.

There is limited literature examining the interaction of Mg, Ca and Sr with hydroxide ions in an aqueous environment. Kluge et al.¹⁰⁷ used gas phase DFT with the B3LYP exchange–correlation functional and found that the introduction of a hydroxide ligand reduces the CN of the Mg²⁺ ion from 6 to 5. In repeat calculations it was found that upon the inclusion of a hydroxide ion, one water molecule migrated to the second solvation sphere causing a change in coordination geometry from octahedral to bipyramidal. The gas phase DFT studies of Felmy et al.¹⁰⁸ explored the hydrolysis of both Ca²⁺ and Sr²⁺ aquo complexes. In hydrolysing calcium aquo species from Ca(H₂O)₆ to [Ca(H₂O)₅OH]⁺, there was no change in first solvation shell CN, and no migration of the hydroxyl group from the Ca ion. However, the removal of a proton from [Sr(H₂O)₆]²⁺ and [Sr(H₂O)₈]²⁺ had a qualitatively different impact. In the [Sr(H₂O)₅OH]⁺ structure the OH⁻ ligand was found to directly bind the Sr ion whereas for [Sr(H₂O)₇OH]⁺ the OH⁻ dissociated from the central ion into the second solvation shell and formed hydrogen bonds with three first shell H₂O molecules. Various XRD,^66,109 ND¹¹⁰ and XANES⁷⁹ investigations examining the structure of Sr(OH)₂·8H₂O all indicated that 8 water molecules coordinate to the ion, and the hydroxide oxygen forms chains of acceptor and donor bonds with the first coordination shell.

The impact of hydroxide on the first solvation shell of uranyl has been documented in both experimental^80,111,112 and computational literature alike.^{80,97,98,113–122} The EXAFS and XRD analysis by Clark et al.¹¹¹ of UO₂(OH)_n²⁻ⁿ (n = 4, 5) found U–O_yl distances of 1.80 to 1.82 Å and U–O_OH distances of 2.21 to 2.26 Å. Multiple DFT studies of [UO₂(OH)₄]²⁻ found longer distances of 1.84¹¹⁴ to 1.88 Ź¹⁶ for U–O_yl and 2.29¹¹⁶ to 2.31¹¹⁵ U–O_OH. The gas phase DFT investigation by Ingram et al.¹¹⁷ into the relative energies and ground state structures of ([UO₂(H₂O)_m(OH)_n]⁽²⁻ⁿ⁾) (n + m = 5) using PBE found that as successive hydroxides are added to uranyl's first solvation shell the U–O_yl distance lengthened from 1.77 to 1.88 Å, the U–O_w distance increased from 2.49 to 2.80 Å, while the U–O_OH distance increased from 2.11 to 2.46 Å.

There is little dynamic data on the presence of hydroxides in the first solvation shell of uranyl. The computational investigation of Austin et al.¹¹⁹ into [UO₂(OH)₅]³⁻ used MD simulations to obtain solvated uranyl hydroxide structures which were then optimised using DFT with the BP86 and B3LYP functionals and a continuum solvation model. This investigation found a U–O_yl distance of 1.88 Å and U–O_OH distance of 2.42 Å. Bühl and Schreckenbach¹²³ used CPMD with an explicit 55 water molecule solvent and an NH₄⁺ counter ion and the BLYP functional to examine the exchange of the axial and equatorial oxygen atoms in [UO₂(OH)₄]²⁻. They found that the structure can be deprotonated to form [UO₃(OH)₃]³⁻ which then undergoes proton transfer via cis-[UO₂(OH)₄]²⁻ complex. The rate limiting step in the transformation is the proton transfer which is assisted by a water molecule from the solvent.

While there is broad consensus in the literature regarding solvation structures of the ions in an aqueous environment, the significant range of reported values for both Ca and Sr hydrates suggests that the systematic study presented here will be a valuable addition to the literature as well as provide an accurate basis for the novel investigation of the hydrated ions in the presence of hydroxide. The impact on the solvation environment of alkaline earth metals due to the presence of hydroxide is investigated here for the first time using AIMD methods, with particular attention paid to identification and characterisation of proton transfer events, along with their timescales and frequency.

2. Computational details

Born–Oppenheimer molecular dynamics simulations were performed using the QUICKSTEP module of CP2K version 3.0 on cells with periodic boundary conditions containing a single cation and up to 64 water molecules, depending on the number of hydroxide species present.^124,125 Temperature and pressure were kept constant using a NPT_I ensemble, where the simulation cell is isotropic with a 0.5 fs time step. Initial cubic cell parameters were set to a = b = c = 11.99 Å, an average temperature T = 400 K was maintained using a Nosè–Hoover thermostat and a barostat maintained pressure of 1 atm.¹²⁶

The Gaussian Augmented Plane Wave method (GAPW) was used for the calculation of forces and energies, which uses a Gaussian basis set with augmented plane wave pseudopotentials.¹²⁷ The Perdew–Burke–Ernzerhof (PBE) generalised gradient approximation¹²⁸ was used to calculate the exchange correlation energy, in keeping with previous studies, including the DFT-D2 dispersion correction as proposed by Grimme.¹²⁹ Whilst the DFT-D3¹³⁰ correction by the same authors provides a superior description of dispersion, this was not available in CP2K in the earlier stages of this study. The calculations used a double-ζ polarization quality Gaussian basis sets (DZVP-MOLOPT-SR-GTH) and a planewave cutoff of 500 Ry.¹³¹ The DFT+U approach was taken for all calculations involving uranium with an effective U_eff = U − J value of 3.96 eV applied to the f orbitals. The U value of 4.5 eV and J value of 0.54 eV is in keeping with many previous studies.^132–137 Charge neutrality was achieved through the use of a uniform neutralising background charge where required.

Each calculated trajectory was 20 ps long and was comprised of 40 000 steps, each of length 0.5 fs. The first 5 ps of each trajectory was treated as an equilibration period, in keeping with previous studies,^{41,63,100,101,138–140} and was not considered in subsequent analysis. It is worth noting that an alternative approach of performing initial classical MD simulations in the NPT ensemble using classical MD may lead to better equilibration.

These trajectories were analysed by considering only every 10th simulation step. Previous literature has analysed the change in coordination number according to the “direct method” proposed by Hofer and co-workers,¹⁴¹ whereby the change in coordination number is only recognised if it lasts longer than 0.5 ps. Our trajectories were analysed to consider changes in coordination number which lasted longer than 0.1 ps in line with our later analysis of proton dynamics. A comparison between analyses conducted using 0.1 ps and 0.5 ps revealed negligible differences.

In order to properly identify the first solvation shell for each ion, a cut-off defining the radial extent of the shell was chosen. Initial radial distribution functions (RDFs) were calculated to obtain the peaks for the first shell M–O bond distances, and the minima after this peak used to define the edge of the shell. This corresponded to a cutoff of 2.7, 3.0 and 3.2 Å for Mg, Ca and Sr respectively, approximately equal to first peak position for each ion plus 0.6 Å. In the case of uranyl, the equatorial U–O peak was used, giving a cutoff of 3.0 Å. These cutoffs were employed in all subsequent calculations.

Evidence for the use of the parameters selected for these simulations, including the RDFs for the aquo complexes are given in Fig. S1–S5 of the ESI.†

3. Results and discussion

3.1. Characterisation of aquo complexes

For each ion, five 20 ps trajectories were simulated using randomly selected DFT optimised initial geometries generated from a 20 ps NPT_I 64 water molecule AIMD run, each ion was added to the centre of the box and the snapshot optimised before being used as the starting structure.

3.1.1 Hydration of Mg, Ca and Sr dications. In order to assess the quality of the model, initial AIMD simulations focussed on the hydration of Mg²⁺, Ca²⁺ and Sr²⁺. Previous similar AIMD studies by Mehandzhyski⁶⁴ and De Leeuw⁴¹ found a CN of 6 for Mg with Mg–O distance of 2.10 Å and 2.08 Å, CN of 7 and 5.9 with Ca–O distances of 2.45 Å and 2.36 Å, and CN of 8 and 6.7 with Sr–O distances of 2.65 Å and 2.6 Å. Their methodology differed with Mehandayski using BOMD with the BLYP functional, a DFT-D2 dispersion correction, and De Leeuw using PBE with no dispersion correction.

Synthesizing RDFs for the entire 75 ps of our simulation data, as shown in Fig. S5 of the ESI,† yields M–O peak values of 2.10, 2.44 and 2.63 Å for Mg, Ca and Sr, respectively. Table 1 summarises the calculated M–O bond lengths and metal coordination numbers for each simulation. As expected, we find an increase in M–O separations from Mg to Sr, along with an increase in coordination number. Our calculated M–O bond lengths are slightly longer than experimentally reported values, however the latter are defined by the M–O RDF peaks.

Table 1 Calculated M–O separations (r_M–O) and mean coordination numbers (CN) for each AIMD trajectory considered. Experimental values are detailed in Table S1 of ESI

Trajectory	r M–O (Å)			CN
	Mg	Ca	Sr	Mg	Ca	Sr
1	2.133	2.499	2.701	6.00	7.45	8.12
2	2.133	2.521	2.702	6.00	7.71	8.19
3	2.136	2.526	2.677	5.99	7.61	7.91
4	2.130	2.487	2.692	5.99	7.11	8.06
5	2.137	2.510	2.688	6.00	7.64	7.83
Mean	2.134 (0.023)	2.508 (0.017)	2.692 (0.009)	6.00 (0.01)	7.50 (0.24)	8.02 (0.14)
Expt.	2.07–2.13	2.39–2.5	2.52–2.67	5.7–6	6–10	6–10

---

The literature values, both experimental and computational, for the Mg–O bond distance cluster around 2.12 Å,^{2,37,44,45,58,142} in excellent agreement with the calculated value of 2.13 Å. Typical experimental values for the Ca–O distance cluster around 2.45 Å,^{35,51,54,56,143} in excellent agreement with our calculated value of 2.44 Å. Our calculated value also compares well to other simulation data, which range from 2.35⁵⁷ to 2.68 Å.⁶² Similarly, experimentally reported values of the Sr–O distance cluster around 2.63 Å,^{37,66,68,70,144} again in excellent agreement with the value of 2.63 Å obtained from our simulations.

Turning our attention to coordination numbers, our calculated values are, in all cases, within the range of values reported experimentally. We see an increase in CN of 2 between Mg and Sr and, interestingly, the largest variation in the value of Ca, as evidenced by a standard deviation of 0.24 in our calculated value. Analysis of the simulated trajectories allows for the residence time associated with each coordination number to be determined. Table 2 shows that Mg²⁺ (which has the highest charge density of the ions considered here) has, almost exclusively, a coordination number of six. The small percentage of time at a CN of 5 likely due to considering changes in coordination number of longer than 0.1 ps (compared to the 0.5 ps of previous studies), although it may be an artefact of the simulation approach taken. Ca²⁺ exists as the hepta- and octa-aquo complex for approximately equal periods of time, indicating significant lability of the eighth coordinating water molecule and explaining the large standard deviation in the calculated coordination number. This eighth water molecule is more easily accommodated by the larger Sr²⁺ which, nonetheless, exists for significant periods of time with coordination numbers of both seven and nine.

Table 2 Percentage of time in which different coordination environments are present, averaged over a total of 75 ps of simulation time for each alkaline earth ion in a pure aqueous environment

Cation	CN
	5	6	7	8	9	10
Mg^2+	0.35	99.65	0.00	0.00	0.00	0.00
Ca^2+	0.00	3.38	45.69	48.29	2.63	0.00
Sr^2+	0.00	0.23	14.09	69.57	15.28	0.83

---

3.1.2 Hydration of uranyl. The calculated average U–O_yl and U–O bond lengths and equatorial coordination number are summarised in Table 3 with the experimental and computational ranges given for reference. The experimental literature data has a range of 0.06 Å for the axial and 0.04 Å for the equatorial U–O bonds, while the computational literature data differs by 0.15 Å for the axial and 0.17 Å for the equatorial U–O bonds.

Table 3 Calculated axial (U–O_yl) and equatorial (U–O_eq) U–O bond lengths, along with mean equatorial coordination numbers (CN) for each AIMD trajectory run. Standard deviations are given in parentheses

Trajectory	U–Oyl	U–O
	r U−Oyl (Å)	r U−Oeq (Å)	CN
1	1.809	2.417	5
2	1.806	2.420	5
3	1.806	2.416	5
4	1.804	2.422	5
5	1.806	2.420	5
Mean	1.806 (0.002)	2.419 (0.003)	5 (0.00)
Expt.	1.70–1.76	2.41–2.45	4.5–5.3
Comp.	1.70–1.85	2.36–2.53	4–5

---

The calculated mean U–O_yl bond distances are longer than those reported experimentally but fall in the middle of the range given by computational literature, comparing excellently with the CPMD⁹⁹ calculated value of 1.81 Å. The mean U–O bond length of 2.42 falls well within the computational literature range of 2.36¹⁰² to 2.53 Å.^94,97 RDFs were generated for the entire 75 ps simulation time, as shown in Fig. 1. The RDF yielded peak positions of 1.80 Å and 2.39 Å, clearly defining the axial and equatorial U–O bond distances, respectively. These values are in generally good agreement with the experimental data, with most literature predominantly reporting U–O_yl distance of 1.76 Å^81,82,84,86 and U–O distance of 2.42 Å.^82,86,88

Fig. 1 Calculated U–O radial distribution function of uranyl in a pure aqueous environment, averaged over 75 ps of simulation time.

In contrast to the data reported for the alkaline earth metals reported in Table 2, no variation in coordination was found in any of the uranyl simulation trajectories. The calculated coordination number of 5 is an agreement with both experimental and computational literature. While some studies^82,86,106 have reported a variation in the coordination number between 4 and 5 with a dominance of the latter, this was not found here. The invariance of the coordination environment can be attributed to the fact that, while uranyl is a dication, the uranium centre is in the +6 oxidation state and so electrostatic interactions with equatorial ligands is expected to be significantly stronger than in the alkaline earths.

3.2 Characterisation of dihydroxide complexes

The agreement of the simulations reported in Section 3.1 with both experiment and previous computational studies gave us confidence that our approach could be applied to the coordination of alkaline earth ions by hydroxide. We have previously considered this in the context of static DFT simulations¹⁴ but, to our knowledge, this is the first dynamical quantum chemical study of hydroxide complexation. Beyond the broad goal of characterising the nature of hydrated hydroxide complexes in bulk environments, our motivation here was to further investigate the recent finding of a low barrier to proton transfer between the first and second solvation shells of strontium dihydroxide complexes, along with the apparent similar stabilities of complexes in which hydroxide resides in the first and second solvation shells.¹³

To minimise bias with respect to initial conditions, we trebled the number of trajectories calculated for the alkaline earth aquo complexes discussed in Section 3.1.1. For each ion, we considered five DFT-optimised starting structures in which either zero, one or two hydroxide species were present in the first solvation shell, giving a total of 15 starting structures. Each of these was then used as the starting point for a BOMD simulation and so, assuming a 5 ps equilibration time for each simulation, a total of 225 ps of analysable simulation time was generated for each ion in the study.

3.2.1 Dihydroxide complexes of the alkaline earths. Here we consider how the coordination of hydrated hydroxide complexes of alkaline earth metals compares to those of aquo complexes. We have previously investigated this for strontium systems containing a partially and fully occupied second solvation shell,^13,14 finding a decrease in total coordination number as hydroxide coordination increases, but these simulations did not allow for migration of the hydroxide species into the bulk.

M–O RDFs for each ion in a dihydroxide environment are shown in Fig. 2. Peaks are found at 2.12 Å, 2.41 Å and 2.59 Å for the Mg-, Ca- and Sr-containing systems, respectively. These values represent a slight increase in the size of the first solvation shell for Mg, which manifests the strongest interaction with the hydroxide species according to coordination numbers, and a slight reduction for the other ions, where the interaction is weaker. Table 4 summarises the calculated average M–O bond lengths, the total and hydroxide CNs for each of the 15 trajectories for each ion. The overall CN, hydroxide CN and bond lengths averaged over 225 ps are also given.

Fig. 2 Calculated M–O RDFs for rare earths in a dihydroxide environment, generated using a total of 225 ps of simulation time ion.

Table 4 Calculated M–O separations (r_M–O), total coordination numbers (CN) and hydroxide coordination numbers (CN_OH) and accompanying standard deviation (SD) for each AIMD trajectory (Traj.) considered in this study

Traj.	r M–O (Å)			CN			CNOH
	Mg^2+	Ca^2+	Sr^2+	Mg^2+	Ca^2+	Sr^2+	Mg^2+	Ca^2+	Sr^2+
1	2.150	2.471	2.705	5.97	6.84	7.99	0.51	0.61	0.41
2	2.151	2.449	2.681	5.97	6.43	7.56	0.88	1.08	0.20
3	2.149	2.470	2.708	5.98	6.61	7.84	0.80	0.75	0.26
4	2.159	2.436	2.676	6.00	6.15	7.44	1.38	0.98	0.25
5	2.155	2.473	2.668	5.98	6.69	7.36	1.39	0.52	0.60
6	2.149	2.488	2.688	5.99	6.97	7.49	0.72	0.46	0.52
7	2.148	2.435	2.667	6.00	6.15	7.42	0.67	1.28	0.11
8	2.127	2.458	2.689	5.57	6.48	7.64	1.23	0.92	0.06
9	2.151	2.470	2.689	5.95	6.63	7.73	0.98	0.78	0.14
10	2.153	2.490	2.686	5.98	6.96	7.51	0.99	0.51	0.77
11	2.141	2.490	2.699	5.67	6.82	7.97	1.56	0.71	0.05
12	2.152	2.472	2.681	6.00	6.55	7.53	0.98	1.03	0.13
13	2.154	2.481	2.679	5.99	6.83	7.61	0.97	0.60	0.15
14	2.155	2.471	2.682	5.85	6.72	7.63	1.29	0.35	0.08
15	2.147	2.479	2.683	5.99	6.82	7.47	0.38	0.32	0.57
Mean (SD)	2.150 (0.008)	2.469 (0.018)	2.685 (0.012)	5.93 (0.13)	6.64 (0.26)	7.61 (0.19)	0.98 (0.34)	0.73 (0.28)	0.29 (0.23)

---

The introduction of hydroxide ions appears to have had an effect on the bonding in the first solvation shell. Compared to the analysis of the aquo environments given in Section 3.1.1 the average M–O bond length for Mg²⁺ increased by ∼0.02 Å, while the average bond length for Ca²⁺ and Sr²⁺ decreased by ∼0.04 Å and 0.01 Å, respectively. Comparing Tables 2 and 5, a moderate reduction in mean coordination number of 0.07, 0.86 and 0.39 is found for the Mg, Ca and Sr ions, respectively. Again, the Ca complex is found to exhibit the greatest lability, with comparable time spent in 6-fold and 7-fold coordination environments. The Sr complex, while showing a propensity for 8-fold coordination, also spent significant time in a 7-fold environment.

Table 5 Percentage of time in which different coordination environments are present, along with mean coordination number, averaged over a total of 225 ps of simulation for each alkaline earth ion. Standard deviations are given in parentheses

Cation	CN						〈CN〉
	4	5	6	7	8	9
Mg^2+	0.13	6.54	93.33	0.00	0.00	0.00	5.93 (0.13)
Ca^2+	0.00	0.21	41.29	52.46	6.03	0.00	6.64 (0.26)
Sr^2+	0.00	0.21	2.67	37.97	54.16	4.99	7.61 (0.19)

---

Since the reduction in coordination numbers exhibited no obvious trend we investigated the contribution to the coordination number from the hydroxide species themselves. The hydroxide ions in the system could be identified by an increased negative oxygen charge and this was used to track the species through the simulation trajectories. Any substantial change in the charge of the oxygen associated with the hydroxide was noted as being potentially indicative of a proton transfer. Then, any change in hydroxide oxygen charge which lasted for less than 0.1 ps was discarded, in order to eliminate “proton rattling” which occurs on the timescale of the Eigen/Zundel interconversion (<0.1 ps).¹⁴⁵ This value is roughly in line with the short bursts of activity seen as part of the Grotthuss mechanism for proton transfer through water.^{17,28,146–149} Although further studies would be useful in order to eliminate all possibility of proton rattling, we also note an absence of the proton returning to the original oxygen within two successive PTEs.

Table 6 reports the hydroxide coordination numbers (CN_OH) and residence times for each ion. Interestingly, while the mean value for Mg is very close to unity, the residence times demonstrate that a simple view of a static complex with a single hydroxide ion in the first solvation shell is incorrect: in fact, this coordination environment was only present for approximately half of the simulation time. A similarly dynamic picture was found for Ca. Here, a lower mean hydroxide coordination number of 0.73 was found due to an increased simulation time in which no hydroxide species coordinated the ion. A very low hydroxide coordination number of 0.29 was found for Sr due to the fact that this ion existed uncoordinated by hydroxide for almost 75% of the simulation time, while simultaneous coordination by both hydroxides was only found during 3% of the simulation.

Table 6 Percentage of time in which different hydroxide coordination environments are present, along with mean coordination number, averaged over a total of 225 ps of simulation for each alkaline earth ion. Standard deviations are given in parentheses

Cation	CNOH			〈CNOH〉
	0	1	2
Mg^2+	23.45	54.95	21.60	0.98 (0.34)
Ca^2+	42.99	41.36	15.65	0.73 (0.28)
Sr^2+	74.22	23.00	2.78	0.29 (0.23)

---

These data can be understood in terms of the reduced charge density as the group 2 elements are descended and the strength of the ionic interaction with both the water and hydroxide species reduces correspondingly. For the charge-dense Mg²⁺ ion, one hydroxide is preferentially bound in the first coordination sphere, replacing a water to approximately maintain the coordination number of the aquo complex. The Ca²⁺ ion also preferentially binds a single hydroxide, but its charge density is insufficient to maintain the coordination number of 7.5 found for the aquo complex, reducing by nearly one. Finally, the charge density of the Sr²⁺ ion is sufficiently low that it is no longer strongly energetically favourable for hydroxide to bind to the ion which therefore has a coordination environment very similar to that of the aquo complex and a correspondingly similar coordination number.

3.3 Uranyl dihydroxide complexes

It was originally intended that an analogous study of hydroxide complexation would be carried out for uranyl, however attempts to optimise initial structures in which the hydroxide species resided outside of the first solvation shell of uranyl were unsuccessful, resulting in proton migration from the first solvation shell to generate two hydroxides directly coordinating the uranyl ion. Six of these optimised structures were nonetheless used as starting points for our AIMD simulations: each simulation generated a 20 ps trajectory, resulting in a total of 90 ps of analysable trajectory time.

The total U–O RDF for the 90 ps simulation time is shown in Fig. 3. In contrast to those of the alkaline earths, this RDF is qualitatively different to that found for the aquo complex, with three well defined peaks at a U–O separation of less than 3 Å. The peak at 1.83 Å corresponds to the axial (U–O_yl) bond whereas the peaks at 2.21 Å and 2.45 Å correspond to interactions with hydroxide (U–O_OH) and water (U–O_w) oxygens, respectively.

Fig. 3 Calculated U–O RDF of uranyl in a dihydroxide environment, averaged over 90 ps of simulation time.

The calculated mean U–O_yl, U–O_w and U–O_OH bond lengths and uranium equatorial coordination number are summarised in Table 7. There is an increase of 0.24 Å and 0.11 Å in the U–O_yl and U–O_w bond lengths, respectively, reflecting the trend reported by both Ingram et al.¹¹⁷ and Cao et al.¹²² The calculated mean U–O_w and U–O_OH bond length are 0.08 Å and 0.07 Å longer, respectively, than those reported by Cao et al.,¹²² and both are 0.07 Å longer than those reported by Ingram et al.¹¹⁷ However, neither of these previous studies included explicit water molecules beyond the first solvation shell, and relied on continuum solvent models to model the long range interactions of the water, which has been shown to impact the accuracy of the first solvation shell in studies with uranyl in water.^91,93,150

Table 7 Mean U–O bond lengths and equatorial coordination numbers for each AIMD simulation. Standard deviations are given in parentheses

Trajectory	r U−Oyl (Å)	r U−Ow (Å)	r U−OOH (Å)	〈CN〉	〈CNOH〉
1	1.827	2.575	2.258	5.00	1.99
2	1.826	2.564	2.246	5.00	2.00
3	1.838	2.549	2.229	5.00	2.00
4	1.827	2.502	2.227	4.49	2.00
5	1.830	2.483	2.213	4.29	2.00
6	1.832	2.485	2.212	4.37	2.00
Mean	1.830 (0.005)	2.526 (0.041)	2.231 (0.018)	4.69 (0.34)	2.00 (0.01)

---

In three of these AIMD trajectories, there is migration of a water molecule to the 2nd solvation shell, reducing the equatorial coordination number to 4. Considering the 90 ps of simulation data in its entirety, uranyl spent 69.18% of the time as a five-coordinated species, and the remaining time four-coordinated. The ion was coordinated by two hydroxide species for 99.77% of the simulation time. Whilst the variation in total coordination number is broadly consistent with the alkaline earth data, the invariance of the hydroxide coordination is again a manifestation of the stronger electrostatic interactions with the hexavalent uranium centre.

3.4 Dynamics of hydroxide coordination

The large standard deviations in hydroxide coordination number of the alkaline earths provide further evidence of the dynamic nature of hydroxide coordination in these complexes. To investigate this further, we explicitly evaluated the total and hydroxide coordination numbers throughout each trajectory in order to understand and characterise any fluctuations found. Representative examples of trajectories analysed in this manner for each ion are presented in Fig. 4 and the remainder can be found in Fig. S6–S14 of the ESI.†

Fig. 4 Total (dotted red) and hydroxide (solid blue) coordination numbers and proton transfer events (△: 1st shell → 2nd shell, ▽: 2nd shell → 1st shell, ○: intrashell) of representative AIMD simulations.

The solid blue and dotted red lines in Fig. 4 indicate the hydroxide and total coordination numbers, respectively, showing significant variation in both during the course of the simulation, particularly for the larger Ca²⁺ and Sr²⁺ ions. The variation in total coordination bears a degree of anticorrelation with the hydroxide coordination, with higher values tending to be found when hydroxide coordination is low, as might be expected based on the comparison of coordination numbers above. There are numerous points in each simulation, however, where hydroxide coordination changes without any change in total coordination, which is indicative of a proton transfer event (PTE). Mulliken charges were evaluated at each timestep in the simulation. This allowed to us to identify when a proton transfer event occurred and, furthermore, to determine whether this event corresponded to proton transfer from the first to second solvation shell (△), second to first solvation shell (▽), or within one shell (○). This data is also presented in Fig. 4, with the position of each symbol indicating the distance between the ion and the oxygen in the water molecule donating the proton. In almost all cases, the change in hydroxide coordination is accompanied by a PTE indicating that the vast majority of the dynamics of the coordination environment is due to PTEs trajectories of each alkaline earth ion. The dashed black line indicates the first solvation shell cutoff distance.

The total number of PTEs found during the simulations is summarised in Table 8. There is a significant increase in the number of PTEs as the group 2 elements are descended, however this increase is not reflected in the number of PTE involving a 1st shell species. In fact, the number of intrashell PTEs in the first solvation shell significantly decreases, although this should be interpreted in the context of the lower probability of finding direct coordination by hydroxide to the larger ions (see Table 6).

Table 8 Number and character of proton transfer events, obtained from 225 ps of simulation data for each alkaline earth ion considered in this study. Percentage values are given in parentheses

Cation	PTE
	1st–1st (○)	1st–2nd (△)	2nd–1st (▽)	2nd–2nd (○)	Total
Mg^2+	15 (4)	77 (24)	79 (24)	155 (48)	326
Ca^2+	9 (2)	80 (22)	82 (22)	200 (54)	371
Sr^2+	6 (1)	59 (14)	55 (13)	306 (72)	426

---

For the alkaline earths considered here the number of PTEs involving transfer into the first solvation shell is approximately equal to the number involving transfer out of the first shell, as would be expected if there was no bias in the starting configurations. The total number of PTEs with this character also drops with increasing ion size and this decrease should again be interpreted in the context of the data presented in Table 6.

Only intrashell PTEs outside of the first solvation shell follow the overall increase reported in Table 8. This can, in part, be attributed to the increased number of hydroxide species outside of the first solvation shell as the group 2 elements are descended but this fails to explain the corresponding increase in the total number of PTEs. Hellström and Behler have demonstrated that proton transfer from water molecules directly coordinating a Na⁺ ion is less likely than that from a bulk water molecule due to an increased energetic barrier to donation in the former and, more generally, that proton transfer is affected by modulation of the hydrogen-bonding environment of the proton donor.¹⁵¹ We suggest that a similar argument may hold here. As shown in Fig. 4, the Mg²⁺ ion structures water into well-defined first and second solvation shells, in which the number of hydrogen bonds is maximised, more effectively than either Ca²⁺ or Sr²⁺. Analogously to the argument above, it would appear that this ordering suppresses proton transfer events, which therefore become more likely for the less charge dense ions. Summarising, we find that, on average, PTEs occur every 0.7, 0.6 and 0.5 ps, in the Mg-, Ca- and Sr-containing systems, respectively.

The results of the uranyl dihydoxide simulations provide a stark contrast to those discussed above: in the entire 90 ps simulation time, only a pair of intershell PTEs were identified: a proton migrated from the 2nd shell into the 1st shell, reducing the hydroxide coordination number to 1. This proton then transferred back out before one from an adjacent water molecule migrated in. The entire process took ∼0.2 ps. The trajectory analysis plots for all 6 trajectories can be found in Fig. S15 of the ESI.†

The almost complete absence of PTEs, with uranyl coordinated by both hydroxides species for the vast majority of the trajectory time again indicates the increased strength of the U–O_OH bond in comparison to the M–O_OH bonds of the alkaline earth. In the latter, we found that the number of PTEs identified inside the first solvation shell was less than the number identified outside. As the uranyl ion is only coordinated by less than 2 hydroxides for a very brief period, the opportunity for intrashell PTEs is almost entirely eliminated.

3.5 Relative energetics of strontium hydroxides

Recently, we reported static DFT simulations of strontium hydroxides using the meta-GGA TPSS exchange–correlation functional, finding that the most stable complex in which both hydroxide species resided in the first solvation shell lay just 3.0 kJ mol⁻¹ higher in energy than that in which only one hydroxide directly coordinated the Sr²⁺ ion. These simulations employed 22 water molecules to explicitly represent the first and second hydration shells of the ion, with the effects of the bulk solvent approximated using the COSMO continuum solvation model. We were interested in establishing if this energy difference could be replicated here with our more realistic description of solvation and so optimised analogous complexes, obtained from snapshots from our dynamic simulations, with appropriate coordination environments. In order to ensure that calculated energies could be directly compared, we fixed the size of the simulation cell to be 11.99 Å in both. We found the complex in which both hydroxides resided in the first solvation shell to be 3.1 kJ mol⁻¹ less stable than that in which only one hydroxide coordinated the ion, in excellent agreement with our previous study.

4. Summary and conclusions

We have reported the results of ab initio molecular dynamics (AIMD) simulations of the Mg²⁺, Ca²⁺, Sr²⁺ and UO₂²⁺ ions in either a pure aquo environment or an environment containing two hydroxide ions.

Mean coordination numbers of the first solvation shell coordinated numbers were calculated, and for systems containing hydroxide species, mean hydroxide coordination was also evaluated. The M–O bond lengths in the first solvation shell of the aquo systems compared very well to existing experimental and computational literature, with bond lengths well within values measured previously and coordination numbers in line with previously calculated values. The expected increase in both average bond length and coordination number was identified when descending the alkaline earths from Mg²⁺ to Sr²⁺.

The accuracy of the results obtained for the aquo complexes gave us confidence in applying the methodology to systems additionally containing hydroxide ions. The addition of hydroxide increased bond lengths in all cases, as previously identified through DFT and AIMD calculations. A reduction in coordination number was also found when two hydroxides were present in the system.

A robust analysis for the hydroxide dynamics in the simulation box over the timescale of an AIMD trajectory was presented and this analysis allowed for the identification of proton transfer events (PTEs) in the dihydroxide systems. PTEs were found to be most prevalent in the strontium hydroxide systems, likely due to the low charge density of the ion and the consequent lack of hydroxide coordination. For all alkaline earths, intrashell PTEs which occurred outside of the first solvation shell were most prevalent, with the numbers of transfers from the first to second shell and vice versa approximately equal. This corroborates the Hellström and Behler¹⁵¹ study which, although using a different methodology, also identified PTEs as most likely to occur away from the ion. Only three PTEs events were identified in the entire simulation data of uranium dihydroxide systems, indicating the clear impact of the increased charge density of the hexavalent uranium ion on the strength of metal–oxygen bonds in aqueous solution. This can be seen as an extreme of the trend found in the alkaline earths, where systems containing the more charge dense Mg²⁺ ion were found to have significantly less PTEs than those containing the less charge dense Sr²⁺.

The results presented here demonstrate that AIMD simulations are able to produce quantitatively accurate structural data for solvated ions and that, with careful analysis, detailed dynamical information can also be extracted, giving highly resolved data regarding the nature of chemical interactions in these systems. Future research efforts will focus on the effects of more complicated ligand environments, along with the energetics and dynamics of ion–surface interactions. In particular, we aim to investigate if and how the frequency of PTEs relates to the ion–surface interactions strength and whether this correlates with experimental data such as the mineral saturation index or dissolution/precipitation kinetics.

Conflicts of interest

There are no conflicts of interest to declare.

References

1. P. R. Smirnov and V. N. Trostin, Russ. J. Gen. Chem., 2007, 77(12), 2101–2107.
2. P. R. Smirnov and V. N. Trostin, Russ. J. Gen. Chem., 2008, 78(9), 1643–1649.
3. P. R. Smirnov and V. N. Trostin, Russ. J. Gen. Chem., 2011, 81(2), 282–289.
4. P. R. Smirnov and V. N. Trostin, Russ. J. Gen. Chem., 2009, 79(8), 1600–1607.
5. H. Ohtaki and T. Radnai, Chem. Rev., 1993, 93(3), 1157–1204.
6. Y. Marcus, Chem. Rev., 2009, 109(3), 1346–1370.
7. B. M. Rode, C. F. Schwenk, T. S. Hofer and B. R. Randolf, Coord. Chem. Rev., 2005, 249(24), 2993–3006.
8. J. Cowan, Biometals, 2002, 15, 225–235.
9. J. J. R. F. da Silva and R. J. P. Williams, The biological chemistry of the elements: the inorganic chemistry of life, Oxford University Press, 2001.
10. E. Bellion, J. Chem. Educ., 1992, 69(12), A326.
11. J. D. Farr, M. P. Neu, R. K. Schulze and B. D. Honeyman, J. Alloys Compd., 2007, 444–445, 533–539.
12. G. R. Bochkarev and G. I. Pushkareva, J. Min. Sci., 2009, 45(3), 104–109.
13. E. Makkos, A. Kerridge and N. Kaltsoyannis, Dalton Trans., 2015, 44, 11572–11581.
14. A. Kerridge and N. Kaltsoyannis, Dalton Trans., 2011, 40(42), 11258.
15. R. Ayala, J. M. Martínez, R. R. Pappalardo and E. Sánchez Marcos, J. Phys. Chem. B, 2012, 116(51), 14903–14914.
16. S. Borah and P. Padma Kumar, Phys. Chem. Chem. Phys., 2016, 18(18), 26755–26763.
17. M. Tuckerman, K. Laasonen, M. Sprik and M. Parrinello, J. Phys. Chem., 1995, 99(16), 5749–5752.
18. I. Tunon, D. Rinaldi, M. F. Ruiz-Lopez and J. L. Rivail, J. Phys. Chem., 1995, 99(11), 3798–3805.
19. M. E. Tuckerman, D. Marx and M. Parinello, Nature, 2002, 417, 925.
20. M. Hellström and J. Behler, Phys. Chem. Chem. Phys., 2017, 19(1), 82–96.
21. O. Demerdash, L. P. Wang and T. Head-Gordon, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2017, 8, e1355.
22. L. M. Varela, T. Méndez-Morales, J. Carrete, V. Gómez-González, B. Docampo-Álvarez, L. J. Gallego, O. Cabeza and O. Russina, J. Mol. Liq., 2015, 210, 178–188.
23. D. Laage, T. Elsaesser and J. T. Hynes, Chem. Rev., 2017, 117(16), 10694–10725.
24. A. Ramĺrez-Solĺs, J. I. Amaro-Estrada, J. Hernández-Cobos and L. Maron, J. Phys. Chem. A, 2017, 121(11), 2293–2297.
25. B. Kirchner, P. J. di Dio and J. Hutter, Multiscale Molecular Methods in Applied Chemistry, 2011, pp. 109–153.
26. M. E. Tuckerman, J. Phys.: Condens. Matter, 2002, 14(50), R1297–R1355.
27. R. Car and M. Parrinello, Phys. Rev. Lett., 1985, 55(22), 2471–2472.
28. N. Agmon, Chem. Phys. Lett., 1995, 244(5–6), 456–462.
29. N. Agmon, H. J. Bakker, R. K. Campen, R. H. Henchman, P. Pohl, S. Roke, M. Thämer and A. Hassanali, Chem. Rev., 2016, 116(13), 7642–7672.
30. D. Marx, A. Chandra and M. E. Tuckerman, Chem. Rev., 2010, 110(4), 2174–2216.
31. A. A. Hassanali, J. Cuny, V. Verdolino and M. Parrinello, Philos. Trans. R. Soc., A, 2014, 372(2011), 20120482.
32. Y.-K. Choe, E. Tsuchida and T. Ikeshoji, Int. J. Quantum Chem., 2009, 109(9), 1984–1990.
33. M. Thomas, M. Brehm, R. Fligg, P. Vöhringer and B. Kirchner, Phys. Chem. Chem. Phys., 2013, 15(18), 6608–6622.
34. G. Pálinkás and T. Radnai, Z. Naturforsch., A: Phys. Sci., 1982, 37(9), 1049–1060.
35. F. Bruni, S. Imberti, R. Mancinelli and M. A. Ricci, J. Chem. Phys., 2012, 136(6), 64520.
36. K. M. Callahan, N. N. Casillas-Ituarte, M. Roeselová, H. C. Allen and D. J. Tobias, J. Phys. Chem. A, 2010, 114(15), 5141–5148.
37. M. A. Marques, M. I. Cabaco, M. I. Marques, B. de and A. M. Gaspar, J. Phys.: Condens. Matter, 2002, 14(32), 303.
38. C. W. Bock, G. D. Markham, A. K. Katz and J. P. Glusker, Theor. Chem. Acc., 2006, 115(2–3), 100–112.
39. C. W. Bock, A. Kaufman and J. P. Glusker, Inorg. Chem., 1994, 33(3), 419–427.
40. J. P. Larentzos and L. J. Criscenti, J. Phys. Chem. B, 2008, 112(45), 14243–14250.
41. D. Di Tommaso and N. H. de Leeuw, Cryst. Growth Des., 2010, 10(10), 4292–4302.
42. T. Ikeda, M. Boero and K. Terakura, J. Chem. Phys., 2007, 127(7), 074503.
43. L. Bernasconi, E. J. Baerends and M. Sprik, J. Phys. Chem. B, 2006, 110(23), 11444–11453.
44. F. C. Lightstone, E. Schwegler, R. Q. Hood, F. Gygi and G. Galli, Chem. Phys. Lett., 2001, 343(5–6), 549–555.
45. R. Caminiti, G. Licheri, G. Piccaluga and G. Pinna, J. Appl. Crystallogr., 1979, 12(1), 34–38.
46. W. Bol, G. J. a. Gerrits and C. L. van Panthaleon Eck, J. Appl. Crystallogr., 1970, 3(6), 486–492.
47. E. Spohr, G. Palinkas, K. Heinzinger, P. Bopp and M. Probst, J. Phys. Chem., 1988, 92(23), 6754–6761.
48. J. N. Albright, J. Chem. Phys., 1972, 56(8), 3783–3786.
49. M. Pavlov, P. E. M. Siegbahn and M. Sandström, J. Phys. Chem. A, 1998, 102(1), 219–228.
50. Y. S. Badyal, A. C. Barnes, G. J. Cuello and J. M. Simonson, J. Phys. Chem. A, 2004, 108(52), 11819–11827.
51. T. Yamaguchi, S. I. Hayashi and H. Ohtaki, Inorg. Chem., 1989, 28(12), 2434–2439.
52. G. Licheri, J. Chem. Phys., 1976, 64(6), 2437.
53. S. Cummings, J. E. Enderby and R. A. Howe, J. Phys. C: Solid State Phys., 1980, 13(1), 1–8.
54. D. Spångberg, K. Hermansson, P. Lindqvist-Reis, F. Jalilehvand, M. Sandström and I. Persson, J. Phys. Chem. B, 2000, 104(45), 10467–10472.
55. V. I. Chizhik, A. V. Egorov, M. S. Pavlova, M. I. Egorova and A. V. Donets, J. Mol. Liq., 2016, 224, 730–736.
56. L. X. Dang, G. K. Schenter and J. L. Fulton, J. Phys. Chem. B, 2003, 107(50), 14119–14123.
57. A. K. Katz, J. P. Glusker, S. A. Beebe and C. W. Bock, J. Am. Chem. Soc., 1996, 118(24), 5752–5763.
58. T. Tofteberg, A. Öhrn and G. Karlström, Chem. Phys. Lett., 2006, 429(4–6), 436–439.
59. M. Kohagen, P. E. Mason and P. Jungwirth, J. Phys. Chem. B, 2014, 118(28), 7902–7909.
60. S. Obst and H. Bradaczek, J. Phys. Chem., 1996, 100(39), 15677–15687.
61. S. Koneshan, J. C. Rasaiah, R. M. Lynden-Bell and S. H. Lee, J. Phys. Chem. B, 1998, 102(21), 4193–4204.
62. M. M. Naor, K. Nostrand Van and C. Dellago, Chem. Phys. Lett., 2003, 369(1–2), 159–164.
63. S. Bogatko, E. Cauët, E. Bylaska, G. Schenter, J. Fulton and J. Weare, Chem. – Eur. J., 2013, 19(9), 3047–3060.
64. A. Y. Mehandzhiyski, E. Riccardi, T. S. van Erp, T. T. Trinh and B. A. Grimes, J. Phys. Chem. B, 2015, 119(33), 10710–10719.
65. P. A. O’Day, M. Newville, P. S. Neuhoff, N. Sahai and S. A. Carroll, J. Colloid Interface Sci., 2000, 222(2), 184–197.
66. P. D’Angelo, H.-F. Nolting and N. V. Pavel, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 53(2), 798–805.
67. L. Axe, G. B. Bunker, P. R. Anderson and T. A. Tyson, J. Colloid Interface Sci., 1998, 199(1), 44–52.
68. D. M. Pfund, J. G. Darab, J. L. Fulton and Y. Ma, J. Phys. Chem., 1994, 98(50), 13102–13107.
69. I. Persson, M. Sandstrom, H. Yokoyama and M. Chaudhry, Z. Naturforsch., A: Phys. Sci., 1995, 50, 21–37.
70. M. I. Cabaço, M. I. De Barros Marques, M. A. Marques, A. M. Gaspar and M. L. De Almeida, J. Mol. Liq., 2005, 117(1–3), 69–76.
71. S. Ramos, G. W. Neilson, A. C. Barnes and M. J. Capitán, J. Chem. Phys., 2003, 118(12), 5542.
72. G. W. Neilson and R. D. Broadbent, Chem. Phys. Lett., 1990, 167(5), 429–431.
73. A. Kerridge and N. Kaltsoyannis, Chem. – Eur. J., 2011, 17(18), 5060–5067.
74. T. S. Hofer, B. R. Randolf and B. M. Rode, J. Phys. Chem. B, 2006, 110(41), 20409–20417.
75. D. J. Harris, J. P. Brodholt and D. M. Sherman, J. Phys. Chem. B, 2003, 107(34), 9056–9058.
76. B. J. Palmer, D. M. Pfund and J. L. Fulton, J. Phys. Chem., 1996, 100(32), 13393–13398.
77. A. Boda, S. De, S. M. Ali, S. Tulishetti, S. Khan and J. K. Singh, J. Mol. Liq., 2012, 172, 110–118.
78. T. M. Seward, C. M. B. Henderson, J. M. Charnock and T. Driesner, Geochim. Cosmochim. Acta, 1999, 63(16), 2409–2418.
79. P. D’Angelo, V. Migliorati, F. Sessa, G. Mancini and I. Persson, J. Phys. Chem. B, 2016, 120(17), 4114–4124.
80. U. Wahlgren, H. Moll, I. Grenthe, B. Schimmelpfennig, L. Maron, V. Vallet and O. Gropen, J. Phys. Chem. A, 1999, 103(41), 8257–8264.
81. C. Hennig, J. Tutschku, A. Rossberg, G. Bernhard and A. C. Scheinost, Inorg. Chem., 2005, 44(19), 6655–6661.
82. L. Soderholm, S. Skanthakumar and J. Neuefeind, Anal. Bioanal. Chem., 2005, 383(1), 48–55.
83. N. Bardin, P. Rubini and C. Madie, Radiochim. Acta, 1998, 83(4), 189–194.
84. P. G. Allen, J. J. Bucher, D. K. Shuh, N. M. Edelstein and T. Reich, Inorg. Chem., 1997, 36(21), 4676–4683.
85. M. A. Denecke, Coord. Chem. Rev., 2006, 250(7–8), 730–754.
86. J. Neuefeind, L. Soderholm and S. Skanthakumar, J. Phys. Chem. A, 2004, 108(14), 2733–2739.
87. C. Den Auwer, D. Guillaumont, P. Guilbaud, S. D. Conradson, J. J. Rehr, A. Ankudinov and E. Simoni, New J. Chem., 2004, 28(8), 929.
88. M. Aaberg, D. Ferri, J. Glaser and I. Grenthe, Inorg. Chem., 1983, 22(26), 3986–3989.
89. C. Hennig, K. Schmeide, V. Brendler, H. Moll, S. Tsushima and A. C. Scheinost, Inorg. Chem., 2007, 46(15), 5882–5892.
90. N. W. Alcock and S. Esperås, J. Chem. Soc., Dalton Trans., 1977, 9, 893–896.
91. B. Siboulet, C. J. Marsden and P. Vitorge, Chem. Phys., 2006, 326(2–3), 289–296.
92. K. E. Gutowski and D. A. Dixon, J. Phys. Chem. A, 2006, 110(28), 8840–8856.
93. N. Kumar and J. M. Seminario, J. Phys. Chem. A, 2015, 119(4), 689–703.
94. P. J. Hay, R. L. Martin and G. Schreckenbach, J. Phys. Chem. A, 2000, 104, 6259–6270.
95. D. Hagberg, G. Karlstro, B. O. Roos and L. Gagliardi, J. Am. Chem. Soc., 2005, 127(6), 14250–14256.
96. C. Clavaguéra-Sarrio, V. Brenner, S. Hoyau, C. J. Marsden, P. Millié and J. P. Dognon, J. Phys. Chem. B, 2003, 107(13), 3051–3060.
97. V. Vallet, U. Wahlgren, B. Schimmelpfennig, Z. Szabó and I. Grenthe, J. Am. Chem. Soc., 2001, 123(48), 11999–12008.
98. H. P. Hratchian, J. L. Sonnenberg, P. J. Hay, R. L. Martin, B. E. Bursten and H. B. Schlegel, J. Phys. Chem. A, 2005, 109(38), 8579–8586.
99. M. Bühl, R. Diss and G. Wipff, J. Am. Chem. Soc., 2005, 127(39), 13506–13507.
100. P. Nichols, E. J. Bylaska, G. K. Schenter and W. de Jong, J. Chem. Phys., 2008, 128(12), 124507.
101. R. Atta-Fynn, D. F. Johnson, E. J. Bylaska, E. S. Ilton, G. K. Schenter and W. A. de Jong, Inorg. Chem., 2012, 51(5), 3016–3024.
102. S. P. Tiwari, N. Rai and E. J. Maginn, Phys. Chem. Chem. Phys., 2014, 16(17), 8060–8069.
103. P. Guilbaud and G. Wipff, J. Phys. Chem., 1993, 97(21), 5685–5692.
104. S. Kerisit and C. Liu, J. Phys. Chem. A, 2013, 117(30), 6421–6432.
105. S. Kerisit and C. Liu, Geochim. Cosmochim. Acta, 2010, 74(17), 4937–4952.
106. N. Rodríguez-Jeangros and J. M. Seminario, J. Mol. Model., 2014, 20(3), 2150.
107. S. Kluge and J. Weston, Biochemistry, 2005, 44(12), 4877–4885.
108. A. R. Felmy, D. A. Dixon, J. R. Rustad, M. J. Mason and L. M. Onishi, J. Chem. Thermodyn., 1998, 30(9), 1103–1120.
109. H. G. Smith, Acta Crystallogr., 1953, 6(7), 604–609.
110. J. S. Ricci, R. C. Stevens, R. K. McMullan and W. T. Klooster, Acta Crystallogr., Sect. B: Struct. Sci., 2005, 61(4), 381–386.
111. D. L. Clark, S. D. Conradson, R. J. Donohoe, D. W. Keogh, D. E. Morris, P. D. Palmer, R. D. Rogers and C. D. Tait, Inorg. Chem., 1999, 38, 1456–1466.
112. H. Moll, T. Reich and Z. Szabó, Radiochim. Acta, 2000, 88(7), 411.
113. G. Schreckenbach, P. J. Hay and R. L. Martin, J. Comput. Chem., 1999, 20(1), 70–90.
114. G. Schreckenbach, P. J. Hay and R. L. Martin, Inorg. Chem., 1998, 37(17), 4442–4451.
115. J. L. Sonnenberg, P. J. Hay, R. L. Martin and B. E. Bursten, Inorg. Chem., 2005, 44(7), 2255–2262.
116. M. García-Hernández, C. Willnauer, S. Krüger, L. V. Moskaleva and N. Rösch, Inorg. Chem., 2006, 45(3), 1356–1366.
117. K. I. M. Ingram, L. J. L. Häller and N. Kaltsoyannis, Dalton Trans., 2006, 2403–2414.
118. G. A. Shamov and G. Schreckenbach, J. Am. Chem. Soc., 2008, 130(41), 13735–13744.
119. J. P. Austin, M. Sundararajan, M. A. Vincent and I. H. Hillier, Dalton Trans., 2009, 5902.
120. V. Vallet, U. Wahlgren, B. Schimmelpfennig, H. Moll, Z. Szabó and I. Grenthe, Inorg. Chem., 2001, 40(14), 3516–3525.
121. S. Tsushima and T. Reich, Chem. Phys. Lett., 2001, 347(1–3), 127–132.
122. Z. Cao and K. Balasubramanian, J. Chem. Phys., 2009, 131(16), 164504.
123. M. Bühl and G. Schreckenbach, Inorg. Chem., 2010, 49(8), 3821–3827.
124. J. Vandevondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167(2), 103–128.
125. J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4(1), 15–25.
126. S. Nosé, Mol. Phys., 1984, 52(2), 255–268.
127. J. Hutter, M. Parrinello and G. Lippert, Theor. Chem. Acc., 1999, 103, 124–140.
128. J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77(18), 3865–3868.
129. S. Grimme, J. Comput. Chem., 2004, 25(12), 1463–1473.
130. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
131. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127(11), 114105.
132. M. Krack, Phys. Scr., 2015, 90(9), 94014.
133. B. Dorado, B. Amadon, M. Freyss and M. Bertolus, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79(23), 235125.
134. J. Rabone and M. Krack, Comput. Mater. Sci., 2013, 71, 157–164.
135. B. Dorado, M. Freyss, B. Amadon, M. Bertolus, G. Jomard and P. Garcia, J. Phys.: Condens. Matter, 2013, 25(33), 333201.
136. S. L. Dudarev, D. N. Manh and A. P. Sutton, Philos. Mag. B, 1997, 75(5), 613–628.
137. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, Z. Szotek, W. M. Temmerman and A. P. Sutton, Phys. Status Solidi, 1998, 166(1), 429–443.
138. J. A. Morrone and R. Car, Phys. Rev. Lett., 2008, 101(1), 1–4.
139. M. J. McGrath, I.-F. W. Kuo and J. I. Siepmann, Phys. Chem. Chem. Phys., 2011, 13(44), 19943.
140. I. Bakó, J. Hutter and G. Pálinkás, J. Chem. Phys., 2002, 117(21), 9838–9843.
141. T. S. Hofer, H. T. Tran, C. F. Schwenk and B. M. Rode, J. Comput. Chem., 2004, 25(2), 211–217.
142. A. Tongraar and B. Michael Rode, Chem. Phys. Lett., 2001, 346(5–6), 485–491.
143. F. Jalilehvand, D. Spångberg, P. Lindqvist-Reis, K. Hermansson, I. Persson and M. Sandström, J. Am. Chem. Soc., 2001, 123(3), 431–441.
144. I. Persson, P. D’Angelo, S. De Panfilis, M. Sandström and L. Eriksson, Chem. – Eur. J., 2008, 14(10), 3056–3066.
145. A. Chandra, M. E. Tuckerman and D. Marx, Phys. Rev. Lett., 2007, 99(14), 1–4.
146. C. A. Wraight, Biochim. Biophys. Acta, Bioenerg., 2006, 1757(8), 886–912.
147. A. Hassanali, F. Giberti, J. Cuny, T. D. Kühne and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2013, 110(34), 13723–13728.
148. P. L. Geissler, C. Dellago, D. Chandler, J. Hutter and M. Parinello, Sci. Mag., 2001, 291(5511), 2121–2124.
149. G. Tocci and A. Michaelides, J. Phys. Chem. C, 2014, 5, 474.
150. M. Bühl, H. Kabrede, R. Diss and G. Wipff, J. Am. Chem. Soc., 2006, 128(19), 6357–6368.
151. M. Hellström and J. Behler, J. Phys. Chem. Lett., 2016, 7(17), 3302–3306.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9cp00142e

---