Computer simulation study of the diffusion of water molecules confined in silica gel

The molecular dynamics simulations are used to study the dynamic properties of SPC/E water molecules confined in silica gel. The model proposed in this study covers three main features that characterize systems of that kind: the excluded volume due to the substrate presence, strong attraction between pore surface and water molecules, randomness of confinement structure. The gel-like structure is modelled as an array of silica composite spheres that are randomly dispersed in volume. Therefore, the spherical SiO2 composite bearing silanol groups on the surface is designed as an element of the porous medium of silica gel. The diffusion coefficients of water molecules are calculated at different temperatures and water densities. To illustrate the model proposed and the behavior of molecules near the pore surface, several snapshots are presented.


Introduction
A knowledge of the properties of the sorbate in specific porous solids is highly valuable in various practical applications.Therefore, a study of liquids confined in porous media is of great interest to both fundamental and applied science.Porous materials being extensively used as catalysts and adsorbents in the chemical industry appear to be a powerful stimulus for statistical-mechanics studies.The computer simulations are widely used not only to make a prediction of the results that can be employed in some further research, but also can be used as the final results in practice, especially when a real experiment cannot be physically performed or is too expensive.
Considering a model of porous media we should focus on the three major features of the porous materials that affect the dynamics and other properties of the adsorbed fluids.These features are related to the volume excluded by confinement, to the type c T.Patsahan, M.Holovko of confinement and to the degree of interaction between the sorbent molecules and the surface of pore surfaces, i.e. they are related to the wetting phenomena.The first one is obvious when we deal with the confinement.The second one can be recognized as the geometrical aspect of the effect of the porous environment on the sorbate while the third one can be associated with the physical and/or chemical aspects of the adsorption processes.A number of simulation studies have been accomplished during the recent years to model and to investigate the properties of the water adsorbed into the porous media.The research primarily concerns the idealized single pore spaces like slit [1,2], cylindrical [2,3] and spherical [4,5].However, more sophisticated models with an atomistically described interaction between water molecules and substrate have been studied as well [6][7][8].
The model of Vycor glass confining water molecules proposed by Rovere et al. [7] and utilized in a number of recent studies [9][10][11][12][13] is a single pore model.It is quite acceptable due to the low porosity found in this kind of glasses where the interaction between the fluid molecules from different pores can be neglected while the details of pore structure can be taken into account.In the particular case, when the hydrophilic groups are located on the surface of the porous silica-like Vycor glass, Rovere et al. [7] used the model of the system considered on the site-site interactions between the adsorbed molecules and the atoms of pore material.This model is based on the one proposed previously by Brodka and Zerda [14] for the study of liquid acetone in silica pores and it is adopted for the case of TIP4P water guest fluid [7] and later for SPC/E water [11].Computer simulation studies of the model proposed by Rovere et al. have proved to be effective in describing the structural and dynamic properties of the fluid adsorbed into a single pore of silica glass.However, in the case of silica gel we deal with a random porous structure built of numerous interconnected cavities of various shapes and different sizes.Obviously, we need a lot of atoms per simulation box to build a silica gel like confinement, causing the simulation study of that system to become very time consuming.Another problem appears when one needs to create pore surfaces of tortuous structure of confinement.
A relatively simple but still nontrivial model for simple fluids adsorbed into a heterogeneous porous solid was formulated by Kaminsky and Monson (KM) [15].Although KM model is simpler than the one used by Rovere et al. [7], it essentially differs from the previously studied models [16][17][18].The last ones describe the adsorption of fluid into non-attractive porous media modelled as hard-sphere matrix, while the KM model includes the attractive fluid-matrix and fluid-fluid interactions.Hence, the effect of temperature on the properties of an adsorbed fluid is present.The fluid-matrix attraction is much stronger compared to the fluid-fluid interaction.Moreover, this model is also characterized by quite a large asymmetry of the diameters of matrix obstacles and of the adsorbed fluid molecules.The asymmetry of energies and diameters may result in a strong accumulation of the fluid molecules in the vicinity of matrix particles.This has already been confirmed by Monte Carlo (MC) computer simulations [15,17].The main attention in the studies [17][18][19][20][21] based on the KM model has been paid so far to the structure and thermodynamic properties.However, there are some studies of dynamic properties as well [22,23].
The main drawback of the KM model is inability to describe the interaction between the silica surface and the polar molecules.The effective potential suggested to be used in KM is obtained as the result of integration through the uniform continuum of interaction sites of silica composite.Hence, the electrostatic contribution vanishes and the hydrophilic surface cannot be described.As far as we consider water molecules in this study, we develop a new model which is based on the two above mentioned models: the first is proposed by Rovere et al. [7] and the second one is KM model [15].We combine these two models to benefit from their advantages and to be able to describe the system of water molecules confined in a porous heterogeneous structure which represents silica gel material.Within the framework of our model, the porous medium consists of a number of spherical silica composites bearing silanol groups on the surface.The silica composite (hereafter called silica particle) represents sui generis a "brick" of our model porous structure.Using a number of such "bricks" we can build any configuration of confinement.In this study an arbitrary configuration of randomly dispersed silica particles is considered.
In this paper we present the description of the proposed model and perform molecular dynamic simulations to study diffusion processes of water molecules confined in silica gel.The mean-square displacement functions are obtained and several values of diffusion coefficients are estimated at the two temperatures.Two levels of hydration are considered in this study.The radial distribution function which describes the molecule layout in the vicinity of silica particle is calculated.Several snapshots of the modelled system are presented.

The model and the details of simulation
The disordered porous medium is constructed in the shape of linked cavities between randomly dispersed silica particles.The surfaces of silica particles represent the structure of a tortuous pore surface which is peculiar to silica gels.Each particle is of a nearly spherical form and consists of the atoms of silicon and oxygen.Besides that, the hydrogen atoms are placed on the surface of a silica particle creating the shell of silanol groups.
In order to make a silica particle, a β-cristobalite crystalline structure of SiO 2 is built.Since we consider a spherical particle, a sphere of ≈ 3.2 nm in diameter is cut off.Further, we remove from the surface all silicon atoms that are bonded to less than four oxygens similarly to [7,8,14].As a result, on the surface we obtain a number of oxygen atoms bonded only to silicon, the so-called non-bridging oxygens (nbO).Oxygens bonded to two silicon atoms are referred to as bridging oxygens (bO).To compensate for the free valences of non-bridging oxygen atoms, hydrogens are attached to each nbO at the equilibrium distance 0.95 Å.After the energy minimization using MM2 force field, the average Si-O-H angle takes a value of 116 • .In our study the obtained number density of silanol groups on the surface is around 6.5 OH/nm 2 and this value is close to some experimental results 5 ± 1 OH/nm 2 [24,25] and 7 OH/nm 2 [26], but it is much greater than the one in P.Gallo [11].It might concern β-cristobalite crystalline structure and convex form of our silica particle.It is known that the surface of amorphous silica has a smaller number density of silanols than the crystalline one, while in the case of β-cristobalite, the surface concentration can reach ≈ 8 OH/nm 2 of silanol groups [27].Moreover, in the study of water adsorption in a mesoporous substrate based on crystobalite silica structure, OH surface density was 7 OH/nm 2 [8].
After the above described transformations, the average diameter of the silica particle is around 2.7 nm.The total number of atoms in silica particle is 840.In figure 1 the resulting form of silica particle is shown.Now the prepared silica particle can be used to build a porous medium of silica gel-like structure by randomly placing its 32 replicas in the cubic box with the size of 9.463 nm and the periodic boundary conditions applied.Besides that, the silica particles are randomly directed to introduce more disorder and to minimize the total electrostatic moment in the system.The atoms composing the porous structure are frozen, i.e. their coordinates are fixed during the simulation run.Thereby, within the model considered we represent the atoms of porous medium of silica gel as a number of immovable interactive sites confining water molecules.In our study we use SPC/E model for water [28].According to [7,14], the site-site interactions between water molecules and silica particles are described by the Lennard-Jones (LJ) and by Coulombic potentials with the appropriate parameters.At that, silicon and hydrogen atoms are interacing with charged water sites only via Coloumbic forces.It should be noted that the use of the original charges for silica sites as was done in [14,7] leads to a large excessive positive charge of the whole system.Thus, in order to preserve electroneutrality of the system the silicon site charge is a little bit reduced .The LJ parameters and fractional charges for SiO 2 sites and SPC/E water are given in table 1.The molecular dynamics (MD) simulations are performed in NVT ensemble for low and high densities of SPC/E water molecules.In this report we consider the following two densities of water: low (ρ w = 0.0206 g/cm −3 ) and high (ρ w = 0.1900 g/cm −3 ), where the high density system corresponds to the hydration level of 22 %.The interaction cutoff distance 9 Åwith long range correction is used.Also, the Ewald summation algorithm for electrostatic interactions is applied.The dynamic properties of water molecules are studied at the temperatures T = 300 and 350 K.The average values of the desired temperatures were maintained using the velocity rescaling procedure according to the weak-coupling scheme of Berendsen [29] with a coupling time constant τ * ≈ 0.5 ps, while the simulation time step, ∆t, was set to 0.001 ps.All MD simulations were performed by the parallel version of the standard DL POLY package [30].The whole simulation box with water molecules introduced into the pores is shown in figure 2. The total number of water molecules in the simulation box is 585 in the case of the low density and 5384 for the high water density.

Results and discussion
The hydrophilic surface of the silica particles leads to the wetting effect near the pore surface.In figure 3 one can see a single silica particle with its vicinity.The water molecules are arranged around the silica particle forming a shell.However, the shell is not concentrated due to a low level of hydration 22 %.The radial distribution functions (RDF) of water oxygens relatively to the silica particle are presented in figure 4. As one can see, the water shell around the silica particle for a low density  of water is remarkably thinner than in the case of a high water density.However, no secondary shell is observed in both cases, while there is some amount of free water molecules that are not located near the pore surface (figure 2).The second hydration level might not be distinguished due to the very rough surface of the silica particle.Large dispersity of the distances between the center of the silica particle and its surface sites (1.2 − 1.4 nm) may lead to the smoothing of peaks of the first and the second hydration levels on the RDFs.From the simulation trajectories, the mean-square displacement (MSD) functions are calculated for the low and high water densities and two temperatures T = 300 K and 350 K (figure 5).The two upper lines correspond to the temperature of 350 K and the other two correspond to 300 K.The minor fits are made to correct MSD data for the low density.Diffusion coefficients are obtained from the MSD functions shown in table 2. From the values obtained, the increase of diffusion with temperature is observed.Although, for a high water density, this increase is slower (figure 6).The same effect was observed in our previous study of the methane in silica gel [23], where the temperature dependency is getting weaker with the fluid density increase.
In spite of the fact that the high density is by an order of magnitude greater than the low density, the diffusion coefficient of the low density water at 300 K is smaller than the diffusion coefficient of the high density at 350 K.It shows that the 50 K increase of the temperature greatly effects the mobility of water molecules in the present conditions.Obviously, the amount of water molecules located in the pore space (free molecules) is marked to be larger when the water density is high compared to the case when it is low.The mobility of the free molecules is much higher than those in the shell around the silica particle.Thus, the contribution of the free molecules in the total diffusion coefficient is essential.It should be noticed that the diffusion coefficients presented here are the averaged diffusion coefficients of all water molecules and they depend only on the ratio of the free water molecules to the molecules in the shell and the total amount of water in the system.The surface and bulk diffusions are not separately considered in this paper due to the complications concerning the determination of hydration shell and the need of a more precise study of structural properties in the vicinity of the silica particle.
The results obtained in this report are mainly intended as a display of the applicability of the above presented model and they appear to be mostly preliminary.As it is seen there is an expediency to continue utilizing this model in further studies.Moreover, it is necessary to consider the system with the higher level of hydration since the amount of water used here is quite low for such a large silica gel porosity we have herein (≈ 0.75).

Figure 2 .
Figure 2. SPC/E water molecules in the model silica gel.

Figure 3 .
Figure 3.The silica particle surrounded by water molecules.The water molecules density is ρ w = 0.1900 g/cm −3 and temperature T = 300 K.

Figure 4 .
Figure 4.The radial distribution function of water molecules around silica particle for low and high water densities.

Figure 5 .
Figure 5.The mean-square displacement data used to determine the diffusion coefficients for water molecules.The solid and dotted lines correspond to the high and the low densities of water molecules, correspondingly.The dashed lines present the linear fit of the MSD data for the low density.

Figure 6 .
Figure 6.The diffusion coefficient of water molecules confined in silica gel.

Table 1 .
Interaction potential parameters for water (SPC/E) and silica particle sites.

Table 2 .
The diffusion coefficients of water molecules in silica gel (m 2 /s).