Molecular dynamics study of aqueous uranyl in hydrophilic mesoporous confinement: the case of slit-like pore in amorphous silica

Molecular dynamics simulations are used to study structural and dynamic properties of water and aqueous uranyl ion adsorbed in a slit-like pore of amorphous silica. Calculations are performed for the flexible SPC/E water model in the atomistically detailed pores with sizes in the range of 2.0–5.0 nm. The hydroxyl groups on the pore surfaces lead to a strong adsorption and strongly affect the mobility of water molecules. The uranyl ion and its aqueous environment adsorbed in the pores are studied at the room temperature for different amounts of water. The effect of hydroxylated silica pores on the formation of uranyl hydrate complexes is discussed within the present study.


Introduction
Besides the fundamental interest to the radionuclides adsorbed in pores of various materials based on SiO 2 its study turns out to be of great practical importance, especially from the ecological point of view.The most applied problems considered are aimed at preventing contamination of large areas caused by radioactive compounds that are capable of migrating through minerals composing the soils and subsurface sediments from the sites where the radionuclides are concentrated: uranium deposits and uranium ore milling sites, nuclear wastes repositories and sites of nuclear accidents.Uranium is a major contaminant due to its predominantly heavy metal content of the spent nuclear fuel.Dissolved uranium dioxide is mostly hexavalent and it can be very mobile in the water environment depending on the material in which it is confined.The dioxide of hexavalent uranium is known as uranyl (UO 2+  2 ) and it is positively charged.Thus, a coulombic interaction is very important in such systems.We present the study of some dynamic and structural properties of uranyl ion and its water environment that are affected by a pore of silica in which they are adsorbed.While one can find numerous experimental results for uranyl ion in various environments, there are only a few investigations where the uranyl ion interacting with some substrate is considered using the methods of computer simulations.A molecular dynamics study of aqueous uranyl near α-quartz surface was presented by Greathouse et al. [1], where a number of radial distribution functions were obtained for different configurations of uranyl and its hydrate and carbonato complexes.The self-diffusion coefficients of water molecules in different layers relatively to the quartz surface were calculated.Later, the computer simulations were performed [2,3] in order to study a sorption of the aqueous uranyl on the clay mineral of montmorillonite.In contrast to the previous studies where the substrate surface had an ordered structure, we focus on an amorphous structure since we are interested in the problems concerned with uranyl ion transport and the properties of its aqueous environment in the so-called lava-like fuel-containing materials formed in the reactor building of Chornobyl Nuclear Power Plant as a consequence of well known heavy nuclear accident which happened in 1986.The lava-like fuel-containing materials have a complicated structure and contain various compounds, but in general they consist of high-radioactive silica glasses.The pores in these glasses form a disordered network of meso-scaled channels, where different kinds of dissolved radionuclides and other light and heavy metals can migrate.Recently, we presented a study of water dynamics in the model disordered porous material [4].However, in the current study we concentrate on the partial case of a confinement, namely on the single-pore model.In this case we neglect the randomness of porous medium for a better comprehension of structural and dynamic properties of water and aqueous uranyl in a separate pore.In this way we can try to better understand different mechanisms of the effect of porous silica on uranyl ion and its aqueous environment.A slit-like pore model of amorphous silica with atomistically detailed structure is chosen.The surfaces of pore are saturated by hydroxy groups which appear due to the chemical reactions of water with non-bonded oxygens on the silica surface.Therefore, the pore is hydrophilic and the adsorption effects are anticipated.We also consider the cases where hydroxyls are partially or completely removed from the surfaces.Several widths of meso-pore are studied to demonstrate a size effect of confinement on the water density profiles and to see how this affects uranyl sorption dynamics.The average values of self-diffusion coefficients are calculated in the hydrophobic and hydrophilic pore of 2 nm size.

Model and simulation details
The pore is modelled as a space between two slabs of amorphous silica.The space between these slabs is filled with water molecules of some density.An uranyl ion is involved in the equillibrated aqueous environment.During the simulation processes the values of (x, y, z)-coordinates of uranyl ion are stored to analyze its trajectories, while the averaged radial distribution functions, such as water density profiles and mean-square displacement functions, were produced during simulation runs.The whole system is atomistically detailed.Therefore, all the pair interactions in the system can be described in a realistic way, although in this case the simulations are quite time consuming and cannot be very long-time scaled.
The model of amorphous silica is prepared using molecular dynamics simulations.The βcristobalite of SiO 2 is chosen as the base atomic structure of silica slab which consists of 5400 atoms and the slab size is about 50 × 43 × 32 Å.The pair interaction between atoms is described by the Morse potential with the parameters given in [5].All simulations are performed in the constant (NVT) ensemble with periodic boundary conditions.Amorphous silica is obtained from crystalline one following the regime 1-XIII in [5], which includes the heating of crystalline structure of SiO 2 to the temperature of 8000 K and a sequential cooling to the room temperature.After equilibration of the system at 300 K during 40 ps, the resulting sample is kept at the temperature of 300 K and at a constant of pressure of 1 atm.The total time of the heating/cooling processes takes 310 ps.At the final stage, in order to make a silica surface relax, the simulation run of the prepared structure is performed during 50 ps without periodic boundary conditions in z-direction.
The next step concerns a preparation of silica surface.We remove from the silica surface all the silicium atoms that are bonded with less than four oxygens like it was proposed in [6,7].To saturate the free bonds of non-bridging oxygen atoms (O h ), hydrogens are attached to each O h at the distance of 1.0 Å and at the average Si-O-H angle 116 • .In our study the obtained number density of hydroxy groups on the surface is around 4.9 OH/nm 2 .This value is in a good agreement with theoretical and experimental data for amorphous silica [8][9][10].Thus, a silica slab is cut off from the obtained sample.Thereby, two identical silica slabs with thickness of ∼ 15 Å and the surface area of 48.29×41.82Å2 are used in order to construct two parallel walls of slit-like pore with the width of d = 2 nm and boundary conditions in x and y directions.The prepared pore can be filled with water molecules of corresponding densities.The different densities of water in the pores are considered.The uranyl ion UO 2+ 2 is involved in the pore as well.To describe an interaction between the atoms in the system, the pair potentials and the partial charges are used in accordance with [1], where the interaction parameters between uranyl ion and water molecules were chosen following Guilbaud and Wipff [11].However, in order to provide an electroneutrality in the whole system the original values of silicon and oxygen atom charges were slightly adjusted.During the simulation of water and aqueous uranyl ion in the pore, the silica structure remains stationary, e.g. the coordinates of silicium and oxygen atoms are fixed, except for hydrogen atoms in hydroxyl groups -OH connected with the corresponding non-bridging oxygen atoms by flexible bonds.In contrast to [1], where the flexible version of SPC model was used as a water model, we prefer a flexible SPC/E one [12].We noticed that flexible SPC model gives highly overestimated values for a self-diffusion coefficient of water molecules, while the self-diffusion coefficient of SPC/E water in bulk is much closer to the known experimental data (2.3•10 −9 m 2 /s), although being a little bit lower (2.0•10 −9 m 2 /s).The bonds between uranium and oxygen atoms in UO 2+ 2 are considered flexible as well [1].The snapshot of simulation box with the model system of aqueous uranyl in the hydrophilic silica pore is presented in figure 1.The potentials and the corresponding parameters for interatomic interactions are presented in table 1 in [1].To study the density profile of water and its effect on transport phenomena in the pore, different pore loadings of water molecules are considered.The water density varies in the range of 0.21 − 1.00 g/cm 3 .Depending on the pore size, the completely loaded pore (1.00 g/cm 3 ) contains 1350, 2024 and 3370 water molecules for the pore sizes of 2, 3 and 5 nm, respectively.The computer simulations of water and aqueous uranyl in the slit-like pore of silica are performed using the method of molecular dynamics at a constant temperature T = 300 K in (NVT)-ensemble with twodimensional boundary conditions.The equations of motion are integrated with a Verlet leapfrog algorithm.The interaction cutoff distance 10 Å with long range correction is used.Also the twodimensional version of Ewald summation algorithm for electrostatic interactions is applied.The average values of the desired temperatures were maintained using the algorithm according to the Nosé-Hoover thermostat with a relaxation time constant t ∼ 0.1 ps while the simulation time step, ∆t, was set to 0.0005 ps.The time of the system equilibration was 100 ps.During the production time, 600 ps, the averaged radial distribution functions and the mean-square displacement functions are calculated.The functions obtained allow us to analyse the structural and dynamic properties of the components of the system.A random removal of certain number of hydroxyl groups from the silica surface changes the hydrophilicity of the pore walls.A complete dehydroxylization leads to the hydrophobicity of the pore due to the absence of polar groups on the surface.

Results and discussion
Considering various densities of water in the hydrophilic pore, the adsorption effect on the pore surfaces is distinguished.This effect can be illustrated by the water density profiles in figure 2. The sharp peaks of the water density profiles indicate that the layering of the water molecules appears near the pore walls for the water densities of 0.21 g/cm 3 and 0.73 g/cm 3 .Due to hydrophilic nature of the surface, all water molecules are attracted towards the pore surface.In all cases one can observe a lower density at the center of pore (z ∼ 0.0) than for bulk water.In the case of a small amount of water (for example, ρ=0.21 g/cm 3 ) the water molecules can be practically missing in the pore center.In figures 3 and 4 one can see the water density profiles depending on the pore size.If the pore size is quite large (more than 5 nm), the water molecules in the center of the pore are not affected by the pore walls and the density profile is uniform in the case of a completely loaded pore.Therefore, the behavior of the water molecules in this region can be treated in the same way as in the bulk.When the pore size decreases, the water density is very inhomogeneous and the distance of the water molecules and the uranyl ion to the wall becomes very important.In this study, two locations of uranyl ion are studied -in the center of the pore and near the pore surface.When uranyl ion is far from the pore surface, it is not directly affected by hydroxyl groups.However, when it is near the surface, the binding of uranyl ion and the corresponding hydroxyl can occur.This binding is accompanied by changes in the hydration shell of the ion.The hydroxyl group is involved into the hydration shell during the binding process, while one of the five molecules of hydration shell leaves the ion.The radial distribution functions that illustrate the mentioned effects are presented in figure 5.The density profile of water molecules has a significant effect on the dynamics of the uranyl ion in the pore.It is noticed that the uranyl ion, placed initially in the pore center, remains in the central area of the pore during the simulation run and moves together with its hydrate complex in the XY surface, e.g.along the pore.The hydrate complex has a pentagonal symmetry, the same symmetry as in the case of uranyl ion in bulk water (figure 5a).The hydrate shell of the uranyl ion consists of five molecules of water.It was noticed that the binding between the uranyl ion and the pore wall can appear when the ion is initially placed at the distance of ∼ 2 Å to the wall.In this case the hydroxyl of silica surface replaces one of the five molecules of water in the hydrate shell of the uranyl.This phenomenon is illustrated by the corresponding averaged radial distribution functions and the running coordination numbers in figure 5b.Despite the attraction between uranyl ion and the pore walls the uranyl ion does not bind with the pore surface by itself when it is far enough ( 4 Å).This can be explained by the adsorption phenomena that lead to the formation of a dense layer of water, which resists the penetration of the uranyl ion towards the wall.In the case of the hydrophobic pore, where the water density is distributed more or less uniformly in the pore volume (figure 2, dashed line), the uranyl ion easily migrates to the pore surface and the close vicinity of the uranyl ion to the pore surface reduces its mobility.In all cases a pentagonal hydrate shell formed at the distance of ∼ 2.51 − 2.52 Å to the center of uranyl ion is observed.Such a configuration has already been observed in the previous studies by molecular dynamics simulation [1] as well as by quantum chemistry calculations [13].
When the water density in the pore increases to 1 g/cm 3 , the uranyl ion initially placed in the pore center rapidly reaches the pore surface area (figure 6) like it was observed in the hydrophobic case before.This can be explained by the density profile uniformity along z-direction (figure 4).However, similar to the previous case, the uranyl ion, being near the pore surface, cannot bind with hydroxyl group by itself due to its hydration shell that needs some charge redistribution, which should appear in such complexes.The method of classical molecular dynamics cannot provide any charge redistribution, because this is a kind of quantum mechanics problem.Thus, some extensions should be applied to take this effect into account.In our study all the charges are fixed.Therefore, the whole adsorption process could not be observed in dynamics.However, when the uranyl ion is placed at the distance ∼ 2 − 3 Å to the silica surface, it binds easily with a corresponding hydroxyl group and keeps in its place for a long time during simulation runs ( 1000 ps).
In figure 6 one can see a number of z-trajectories of uranyl ion, which is adsorbed in the pores of different sizes (2-5 nm).It is shown that in the narrow pore (2 nm), the uranyl ion moves fast towards the pore wall.When the pore increases the averaged velocity of uranyl ion in z-direction decreases.Moreover, a time delay appears before the ion begins to move to the wall.When the pore size is of 5 nm, the uranyl ion located in the middle of the pore becomes indifferent to the pore walls.Moving around, the uranyl ion can approach one of the walls, but the time needed to do that is much longer -more than a few orders of magnitude larger than in the case of 2 nm and 3 nm.A tendency to the motion of the ion towards the wall was not observed even after 1000 ps.Thus, we should determine the effect of the pore hydrophilicity on the water dynamics, the averaged values of the total self-diffusion coefficient of the water molecules in the pore of 2 nm and the different surface densities of hydroxyl groups.For this purpose the functions of mean-square displacement (MSD) are calculated (figure 7).The three values of the surface densities of hydroxyls are chosen, that correspond to the full hydroxylated surface (100% or 4.9 OH/nm 2 for the system considered), partial hydroxylated (2.5 OH/nm 2 ) and completely dehydroxylated surface (0%).The self-diffusion increase is observed when the surface density of hydroxyls decreases.The dashed lines denote MSD functions in z-direction.This particular form of the MSD functions is appropriate to the case of the molecular dynamics in confined space, where the motion in some directions is restricted by the pore walls and the maximum value of the MSD function -which is reached in this case -defined by a pore size.The calculated values of the self-diffusion coefficient are close to that predicted in [14,15] for water adsorbed in Vycor-like glass at the similar conditions.

Conclusions
The model of a slit like pore of silica is considered within this study.The water density profiles are obtained to study an effect of the pore size and its loadings on the uranyl ion dynamics.If the pore is narrow (2-3 nm) and the loading is close to 1 g/cm 3 of the water density, the uranyl ion adsorbs very fast on the pore surface.In the hydrophobic pore the uranyl ion reaches the surface of a partly loaded pore faster than in the hydrophilic one.However, in the hydrophobic pore, where hydroxyl groups are absent, the ion cannot associate with the surface.
The significant effect of the arrangement of the water molecules in the pore on the uranyl ion behavior and its hydration shell inspired us to consider them depending on the pore parameters.The layering structure of the density of water molecules near the pore surface is observed in the hydrophilic pore.The first peak in the water density profile becomes lower when the surface density of hydroxyls decreases.Moreover, the water density distribution through the pore volume is more homogeneous when there is a small amount of the hydroxyl groups on the surface.In the case of partially loaded pore, the water molecules are practically missing in the center of the pore .When the pore surface is hydrophobic (small amount of hydroxyl groups), the difference between the densities of water molecules in the center and near the surface of the pore is much smaller.In the pore that is partially loaded with water molecules, one can see a high inhomogeneity of water density which is caused by adsorption of water molecules on the hydrophilic pore surface.Due to a small amount of water that is incapable of filling the whole pore, one can notice an asymmetry in the water density profile.An asymmetric profile indicates that the most part of water molecules prefers one of the two pore walls in the time interval observed.In the case of a pore completely loaded with water molecules (ρ = 1 g/cm 3 ), the average self-diffusion coefficient of water decreases when the surface density of hydroxyl groups increases.
The aqueous uranyl ions confined in the silica pore are attracted to the pore walls and bind with the silica surface.Due to a strong association, the uranyl ion remains on the silica surface for a long time.A reverse process (e.g.disconnection of uranyl ion) was not observed during simulation runs (> 600 ps).Near the silica surface, the hydrate shell of uranyl ion consists of four water molecules and one hydroxyl group bearing on the surface.A complex of four water molecules and one hydroxyl group was obtained only when uranyl ion was placed at the distance of ∼ 2 Å from the silica surface.
The report presented here is only preliminary and should be extended and perhaps studied more in detail in some directions.However, the obtained results give us a flavor of the structure as well as show the dynamics of the system considered.The problems that emerged during this study suggest us to make some modifications related to the model of the system as well as to the very simulation process.The first methodical problem concerns the charge redistributions during the binding processes of uranyl ion and water/hydroxyl that should be taken into account in order to correctly describe the adsorption process.Another problem consists in the large time consumption while simulating the dynamics in such systems.Being a very heavy compound compared to water molecules, the uranyl ion has an extremely slow dynamics that cannot be observed properly during relatively short-time scale simulations (∼ 1000 ps) within the framework of the model considered in this study.Therefore, a more simple model of silica structure is needed, where the number of sites of interaction are much less.

Acknowledgements
TP thanks the World Federation of Scientists and ICSC "World Laboratory" for the partial financial support within the framework of the National Scholarship Program.The computational facilities were provided by the computing cluster of the Institute for Condensed Matter Physics.

Figure 1 .
Figure 1.A snapshot of the simulated system.The uranyl ion (in the center) and water in the slit-like pore of amorphous silica.

Figure 2 .
Figure 2. The water density profiles for different pore loadings.The dashed line corresponds to the case of hydrophobic pore at water density of ρ = 0.73 g/cm 3 .

Figure 4 .
Figure 4.The water density profiles for different pore sizes (2-5 nm) and different surface density of hydroxyl groups.Dotted lines correspond to the partly hydrated surface of the pores.

Figure 5 .
Figure 5.The radial distribution functions U-Ow for the aqueous uranyl and its hydration shell in a slit-like pore of amorphous silica at the temperature 300 K: (a) uranyl ion in the center of the pore; (b) uranyl ion near the surface of pore wall.

Figure 6 .
Figure 6.z-trajectories of uranyl ion in the pores of different sizes and loadings.

Figure 7 .
Figure7.The mean-square displacement functions and the total self-diffusion coefficients of water molecules in the slit-like pore of amorphous silica for different levels of the pore hydrophilicity.The pore size is equal to 2 nm.