Structure and dynamics of water in TiO2 nano slits:The influence of interfacial interactions and pore sizes

Mingjie Wei*,Yong Wang

State Key Laboratory of Materials-Oriented Chemical Engineering,and College of Chemical Engineering,Nanjing Tech University,Nanjing 211816,China

Keywords:TiO2 Water Molecular dynamics ReaxFF Nanopore Transport

ABSTRACT The interaction of water with TiO2 surfaces is of enduring interest because of wide applications of the TiO2 materials in aqueous environments.The structure and dynamic properties of water molecules in TiO2 nanopores are crucial as increasingly TiO2 materials are synthesized into nanoporous structures.In this work,the structural and dynamic properties of water molecules in nanoscale slit pores of TiO2 are investigated,by using three sets of force field models for the water-TiO2 interaction,as well as four TiO2 slit pore widths.It is concluded that the water-TiO2 interaction dominates the interfacial structure of water molecules,while the dynamic properties of water molecules are primarily influenced by the slit width in both interfacial and central regions.These findings indicate that both of the fluid properties and the interactions of fluids with pore wall will determine the transport properties of fluid in nanopores.If the pore size is large enough,e.g.1.0 nm or larger in this work,the transport properties will be determined most by the fluids themselves.For the cases of pores whose sizes are in the range of interfacial region,the influences of pore size and interfacial interaction will interfere each other.

1.Introduction

Titanium dioxide(TiO2)has been the subject of heightened scientific investigations since Fujishima and Honda published their discovery on UV-induced redox of water on TiO2[1].TiO2is typically applied in solar cell and hydrogen generation as a heterogeneous photocatalyst due to its semiconductor property,so that reactions at the interface are the primary focus of investigations[2].Furthermore,since water molecules are involved with almost all of the applications of TiO2materials,the interaction of TiO2with water molecules plays a critical role in most applications [3].A number of studies on the surface reaction of water molecules at various TiO2surfaces have been carried out by experimental or computational methods [4–12].For example,the dissociation of water at terminal Ti(Ti5c)sites of rutile(110)surfaces under ultraviolet condition was directly observed by scanning tunneling microscope (STM) [13].The dissociation of water at TiO2surfaces was also investigated by the quantum chemistry calculations[14,15].

While the dissociation of water at TiO2surface is currently an issue of some controversy,the transport of the reactant and product in the vicinity of reaction sites has received less attention,despite their influence on the equilibrium and kinetics of the surface reactions.Since larger surface area provides more sites for reaction,nanoporous TiO2materials have been synthesized with large specific surface area [16].Compared to the open surfaces,it is difficult for experimental methods to characterize the motion of fluid molecules in nanopores [17–19];however,molecular dynamics (MD) simulations can provide insight into in-pore structure and dynamics[20],much as it has for TiO2/water and TiO2/aqueous solution interfaces[21].The success of MD simulations on water in carbon nanotubes indicates that the structural,dynamic and phase equilibrium properties of water could be studied effectively by MD simulations [22,23].Moreover,the influence of pore wall chemistry on the transport properties of fluid can be predicted by MD simulations [24–26],and those predictions would be verified by experimental findings [27].Hence,in this work,MD simulations are applied to investigate the structure and dynamics of water in TiO2nanopores.

The water/TiO2surface interaction clearly plays an important role in the structure and dynamics of water in TiO2nanopores[28].In our previous work [29,30],the effect of surface chemistry on the fluid properties of fluid was considered by comparing two different types of surfaces (hydrophilic– TiO2– and hydrophobic–carbon-covered TiO2).It was found that the interfacial interaction influenced the motion of the water molecules through the hydrogen bond network and that there was strong interaction of water molecules with the TiO2surfaces.The water/TiO2surface interaction in that work was the model derived by Bandura and Kubicki(which we will refer to as the BK model) [31].Alimohammadi and Fichthorn [32]claimed that the Bandura/Kubicki water/ TiO2surface interaction overestimated the true interaction,and so introduced a weakened version of the BK interaction (which will be referred to as the AF model).In this paper,we focus on the BK and AF models and their impacts on the structure and dynamic properties of water molecules in TiO2slit nanopores.As a third model for the water/TiO2surface interaction,we will consider the ReaxFF reactive force field [33]as well,which is developed to describe the atom interactions while there is chemical reaction occurs in MD simulations.In this work,the surface of TiO2would react with water molecules based on ReaxFF,and the interfacial layer of water might have distinct properties.

Simulations of the structural properties of water in cylindrical[34]and slit [35]pores have been previously reported.Given the high degree of characterization of the TiO2(110) surface,we consider a slit pore composed of two (110) surfaces as our model,although we note that recently amorphous TiO2nanotubes have been prepared for use in applications [16,36].As MD simulations of water in carbon nanotubes have shown,the structure of water is evidently influenced by the pore size[37].In order to investigate how the confinement affects water structure and dynamics in the slit pores under flow,we performed simulations of water confined in nanoscale TiO2slit pores with various pore widths.

2.Simulation Details

Similar to our previous work[35],the slit model is used to simulate the TiO2pores filled with water molecules,but the slit size is much smaller in this work.The slit model is composed of a slab of rutile TiO2with the (110) faces exposed to the fluid.The surface size is 2.60 nm (x) by 2.37 nm (y),above which there is a slit pore of w 0.8,1.2,1.6 or 2.0 nm(described by the coordinate z).The slit pore is filled with water molecules.Using periodic boundary conditions in all three dimensions,an infinitely long and wide slit model is built (shown in Fig.1).It is noted that the slit width was measured by the distance between the planes of bridging oxygen atoms (BO) rather than the planes of five coordinated Ti (Ti5c)atoms in this work.

As noted in the introduction,three sets of force fields that describe the interaction between water molecules and the (110)surface of TiO2in the rutile structure are adopted in this work.The first one is that derived by Bandura and Kubicki [31](BK)and the second one is the minor modification by Alimohammadi and Fichthorn[32](AF)based on a reformulation of the BK parameters.The third is the reactive force field (ReaxFF) optimized for water/(110)rutile surface of TiO2[38,39].The interaction between O atom of water and Ti and O atoms of rutile for classical MD simulations(BK and AF)are listed in Table 1,and the ReaxFF potential file was obtained from the supporting information of Kim et al.[39].Since Alimohammadi and Fichthorn believed the BK potential model overbinded water on TiO2surfaces [32],they modified the geometric value which describes the interaction of titanium (Ti)and oxygen (O) atoms of rutile and oxygen atoms of water (Ow)in the AF potential model (shown in bold in Table 1).The water model used in both BK and AF simulations is SPC/E[40],consistent with the water model choice of Predota,Cummings and co-workers [7,21,35].The SPC/E is among the best of the non-polarizable models for water as it reproduces the vapor–liquid phase envelope as well as other structural and transport properties at ambient conditions [40,41].

Fig.1.Snapshot for simulation box (2.0 nm slit shown only):Pink,Ti;Cyan,O of rutile;Purple,O of water;Lime,H.Red,green and blue axes represent x,y and z direction,respectively.(a) The rutile (110) surface exposed to water.(b) The first step of 10 ns simulations with obtained water density in slits,3D periodical boundary condition used to form the slits.(c) The first step of 1 ns NPT simulation to obtain precise water density in the slit.

Since the slit widths (0.8–2.0 nm) considered here are smaller than those in most other studies [21,35],the precise density of water in the slits is crucial.For this reason,the isothermalisobaric(NPT) ensemble is used to determine the density of water in the slits.However,if the typical three-dimensional NPT ensemble algorithm is applied in all directions,the slit width cannot be preserved.To solve this problem,we applied NPT ensemble only in one direction (the y direction in this work) while placing two water reservoirs outside the slits in this direction (as shown in Fig.1c).After 1-nanosecond (ns) of simulation time,the average density of water in the slits can be calculated (densities and the numbers of water in slits are listed in Table 2).The pressure in the NPT simulations was set to 100 bars (1 bar=0.1 MPa) to avoid the appearance of nano bubbles.By this method,the system was minimized to save computational resources.

After initializing the water in the pores,10 ns simulations were performed without the water reservoirs (i.e.,restoring periodic boundary conditions in all three directions,as shown in Fig.1b).The atoms in the rutile slabs were frozen to keep the slit widths unchanged during the simulation and the thermostat was coupled with water molecules at 300 K.All simulations were performed with LAMMPS [42].

Table 1 Potential parameters used in classical MD simulations.

Table 2 Densities and numbers of water molecules in the slits.

3.Results and Discussions

3.1.Interfacial density

Although the interfacial structure of water near the rutile(110)surface was well investigated by experimental and computational methods,we herein try to find out whether the slit width or potential model will affect the interfacial structure of the fluid in this paper.The density profiles of water oxygen (Ow) atoms were calculated to study the interfacial structure of water in slits (shown in Fig.2).In Fig.2(a-c),the first and second peaks coincide independently of the slit width,indicating that the interfacial density profiles of water are not influenced by the slit width.But in the narrowest slit (D=0.8 nm),the third peak is slightly different,because there is no enough space for the narrowest slit to form a complete third layer.

Fig.2.Density profiles for water molecules in slits.Here we show only the left half part of the density profiles.The 0.0 at the x-axis represents the plane of the surface formed by bridging O atoms.(a) Results from the BK potential simulations (solid lines).(b) Results from AF potential simulations (dashed lines).(c) Results from the ReaxFF simulations (dotted lines).(d) Comparison of interfacial structure for various potentials in the slits of 2.0 nm.

By comparing the results from various potentials for the same slit width(shown in Fig.2d),it is observed that the interfacial density profiles are different from each potential model,which indicates that the interfacial density is affected by the potential model.However,in the central region,the densities of water for various potential models are trending to be almost the same,as if the effect of potential models is shielded by the interfacial structure of water.

If the results of interfacial density are compared with the geometric values in potential models of BK and AF,it will be found that these geometric values take effect much.The large geometric value results in the farther position of peaks in density profiles.In order to clarify this,the pair density profiles between surface titanium atoms (Ti5c) and oxygen atoms of water (Ow) were calculated,as well as those between surface bridging oxygen (BO) atoms and water oxygen atoms(shown in Fig.3,more details could be found in Fig.S1).The position of the first and the second peaks in Fig.3a stands for the distance between Ti5cand the Ow atoms in the 1st and the 2nd layer of water,respectively.Hence,the results from Fig.3a are quite similar to those in Fig.2d.Moreover,similar to the results from density profiles,the large geometric value gives the longer distance.These results reveal that the geometric values in various potentials do affect the interfacial density dramatically.

Fig.3.Pair density profiles of Ti5c -Ow (a) and BO -Ow (b) in the slits of 2.0 nm width.Solid,dashed and dotted lines denote results of the BK,AF and ReaxFF simulations,respectively.

In Fig.3b,the positions of the first peaks of BO -Ow curves are 0.260,0.300 and 0.325 nm for BK,AF and ReaxFF,respectively.From the experimental work of x-ray crystal truncation rod (CRT)[10],the first BO -Ow peak appears at (0.242 ± 0.005) nm.The results of BK models predict a little larger distance than the experimental result but shorter than the one of Ow-Ow in bulk water at room temperature (0.28 nm) [43].That means the interaction between BO atoms and water being stronger than the one between water molecules themselves in the models of BK.For AF potential models,the distance between BO and Ow atoms are larger than the one between bulk water oxygen atoms,which suggests the weaker interaction of water molecules with the surface as this potential model planned to act.More discussions about the interaction between water molecules and the surface can be found in the section of Hydrogen bond lifetime below.

The first peak for ReaxFF potential model in Fig.3b is much larger than the two other’s,which predicts the weaker water– surface interaction.However,there is a shoulder near the surface.By analyzing the dissociation of water molecules at the surfaces,it is found more than 98%of 1st layer water molecules are dissociated,so that the bridging oxygen atoms transfer into hydroxyls,which prevent the 2nd layer water getting close to the bridging oxygen atoms.The shoulder in the dotted line of Fig.3b denotes the molecular adsorption of water molecules that are no more than 2%.Since the dissociation of water is not the purpose of this paper,we will not discuss this topic further.More discussion about dissociation of water molecules on rutile(110)surface could be found in the paper of Wesolowski et al.[44],in which much more evidence about dissociation was presented.

3.2.Dipole orientation profile

In the previous section,it is obvious that various interaction potential models predict different interfacial density,which indicates that the interfacial density profiles are governed by the potential models.To examine this dependence more thoroughly,dipole orientation profiles of water molecules are calculated and compared in Fig.4.The definition of the angle θ could be found in our previous papers [29].

Similar to the results of interfacial density,the dipole orientation profiles for different pore widths almost coincide in Fig.4(ab).It is concluded that the width of the slit takes nearly no effect on the dipole profile of water molecules in the slits,except the narrowest one(D=0.8 nm).As has been already discussed above there is no enough space for water molecules to form the complete 3rd layer in the narrowest pores.Thus,the dipole orientation of the 3rd layer water in the slits of 0.8 nm differs from the others.Above the 3rd layer,the dipole orientations become random suggesting that the properties of water molecules in the central regions are close to the bulk ones.These results confirm that the effect of interfacial potential model is shielded by the interfacial structure of water molecules again.

By comparing the results of various potential models(Fig.4c),it is found that the dipole orientation profiles are a little different.This confirms that the interfacial structure of water molecules strongly depends on the potential models chosen to describe the interaction of water with the surface.Since the water molecules are not rigid in the ReaxFF simulations,it is hard to define the dipole orientation for single molecules.Hence,the dipole orientation profiles for ReaxFF potential model are not shown here.

3.3.Average hydrogen bond number (NHB )

Fig.4.Dipole orientation profiles of water molecules near wall.The 0.0 at x-axis denotes the plane of surface bridging O atoms.Horizontal dashed line represents the value of cosine of the angle is 0.0,where angle is composed of dipole moment of water molecules and z-axis.The vertical dotted lines are the edge of 1st and 2nd layer of water.(a)Results from BK potential simulations (square symbols).(b) Results from AF potential simulations (circle symbols).(c) Comparison of interfacial structure for various potentials in the slit of 2.0 nm.The x-axis of AF is modified to fit the 1st and 2nd layer edge.

Hydrogen bond network is one of the most important structure properties of water.To investigate the effect of potential and slit width on water structure,the NHBin interfacial and central regions of the pores is calculated by the geometric definition of hydrogen bonds [45,46],based on three criteria [47].The distance of O－O atoms is less than 0.35 nm,the distance of O－H atoms less than 0.245 nm and the angle of H···O－H less than 30°.The NHBin interfacial and central regions of the pores is calculated by the geometric definition of hydrogen bonds [47].The NHBfor the interfacial region (between the 2nd layer water and the 1st layer water or BO atoms) and the central region (water molecules outside the 2nd layer) is plotted as empty and filled symbols in Fig.5,respectively.

All NHBvalues for each potential model keep almost the same regardless of the slit width in both interfacial and central regions.By comparing the absolute values of NHBfor various potentials,it is obviously that the potential model strongly influences the amount of hydrogen bonding in the interfacial region.This could be induced by observing that the three dashed lines parallel to each other but with various values in Fig.5.These results also support the conclusion that the interfacial structure of water is strongly dependent on the potential models while the slit width could not affect the interfacial water structure.In Fig.5,all filled symbols have very near values,no matter with the potential models.This once again confirms that the effect of potential models is shielded by the interfacial water layers.More details on the NHBfor interfacial water molecules could be found in Fig.S2.

Fig.5.Average hydrogen bond numbers are plot as function of slit width.Filled symbols represents the NHB for water molecules in central region.Empty symbols are those of the interfacial region.

3.4.Diffusion coefficients

The results shown so far indicate the interfacial structure of water molecules is governed by the potentials rather than slit width,and the structure of water molecules outside the 3rd layer are closed to the bulk ones.In the coming sections,the dynamic properties of water molecules will be discussed.Since all simulations are slit modeled,the diffusion along the plane parallel to the surfaces (x &y) is more meaningful than the one along the direction normal to the surface(z).Hence,the self-diffusion coefficients of water molecules are calculated only in xy plane by

where MSD denotes the Mean square displacement(MSD).In order to better fit the MSD curve,we make a conversion of equation (1)into

From the results given by filled symbols,it is obvious that the diffusion coefficients increase with the slit width.Since the value of the diffusion coefficient is calculated by averaging the MSD of all water molecules in the slits and the 1st layer of water is immobile [35],such increasing might be due to the larger proportion of water outside 1st layer in larger slit pores.Thus,the diffusion of water molecules is calculated in the volume excluding the 1st layer(empty symbols in Fig.6).The absolute Dxyvalues are a little larger for the water molecules outside the 1st layer.However,the tendencies keep the same,indicating that the diffusivity of water in slits is mainly influenced by the slit width.The diffusion coefficients of water molecules by layers are shown in Fig.S3,which confirms the diffusivity of water in the central region changing with the slit width.

Fig.6.Diffusion coefficients (along xy plane) of water molecules in the slits.The Filled symbols denote results for the water in the whole pore volume,while the empty symbols denote results calculated for the volume with the 1st layer of water excluded.

In Fig.6,the various potentials give very similar diffusivity.It seems that the diffusivity of water is not strongly dependent on the interaction potential model,but mainly depends on the width of the slit.As has been discussed in the previous sections of interfacial density and dipole orientation profiles,the influence of potential model cannot extend beyond the 3rd layer of water(about 0.5 nm away from the surface),so the structure properties of water in the central region of pores are close to those of bulk water.However,the diffusion coefficients observed here are much different from the one of bulk water(2.75×10-9m2·s-1at 300 K)[48].The diffusion of water in these slits is limited by the nano confinement not only in z direction but also in the xy plane.The possible reason is that large slits provide more space for water molecules (excluding the interfacial ones) to move in 3D dimensions,which increases the diffusivity.

3.5.Hydrogen bond lifetime

The results for the diffusion coefficients indicate that the slit width does influence the diffusivity of water,while it does not significantly affect the structure of water(in either interfacial or central region).As was discussed above,the space provided by different slits(with different width)might be the key factor,which governs the behavior of the diffusion coefficient.To confirm this,the hydrogen bond lifetime is calculated from the autocorrelation function

where si（t ）=｛0,1}is the existence function of the hydrogen bond i at time t.The integral of C（τ ）gives a n estimate of the average hydrogen bond lifetime [49].

Generally,in the central regions,the hydrogen bond lifetime of water molecules drops as slit width increases (Fig.7a),indicating weaker interaction between water molecules and correspondingly stronger motion in large slits.This is in agreement with the results of diffusion coefficients that water molecules in larger slits have higher diffusion coefficients.In the section of average hydrogen bond numbers,the values of NHBof the central region are very close regardless of the slit width,but the hydrogen bond lifetimes change regularly with slit width.This confirms that the dynamic properties of water are influenced by the slit width dramatically.

It is noticeable that the values for interfacial region (Fig.7b-c)have the same tendency with those in central region.This indicates that,in the interfacial regions,the dynamic properties of water molecules are also affected by the slit width.The diffusion coefficients of the 2nd layer water molecules confirm this as well(Shown in Fig.S4).

Fig.7.Hydrogen bond lifetime is plotted as a function of slit width.Solid marks with dotted lines represent the hydrogen bond lifetime for water molecules in central region.Hollow marks with dashed lines are those in the interfacial region.

From Fig.7b-c,it is clear that the hydrogen bond lifetimes of water– surface in the BK potential models are larger than those in AF and ReaxFF potential models are,indicating a stronger interaction of water with surfaces in the BK potential model.The results of BO and 2nd layer of water (Fig.7c) are in agreement with the results of pair density distribution shown in Fig.3b.

It should be noted that the values hydrogen bond lifetime obtained in this paper is much higher than those in other works are [47,50–55],which possibly originate from the different methods of calculating hydrogen bond lifetime.In this work,the values are compared only within our calculations,so as to focus more on the tendency of the hydrogen bond lifetime depending on the slit width,rather than its absolute value.

No matter of the interfacial interactions,the hydrogen bond lifetime is evidently larger for the cases of slit width being 0.8 nm,indicating the diffusion mode is ‘‘plane-whole”diffusion,similar to the ‘‘single-file”diffusion in CNTs [56].While slit width is large enough,the diffusion mode would turn to Fickian diffusion as water molecules can move freely like in bulk phase.

4.Conclusions

In this paper,we used classical and reactive force field MD simulations to study the effect of potential model and slit width on the structure and dynamic properties of water in slits of rutile TiO2.By comparing the results from three sets of potential and four various slit widths,it is found that both structure and dynamics properties of water molecules in the interfacial region(within around 0.5 nm away from the surface)are governed mainly by the potentials chosen to describe the interaction between water molecules and surface atoms.Outside the interfacial region,the structure and dynamics properties of water molecules are independent of the interfacial potential models and closed to bulk one,generally due to the shield effect of the interfacial water layer.While the structure properties of water are almost independent of the slit width,the effect of slit width only works on the dynamic properties of water molecules in both interfacial and central regions of the slits.That illustrates that the geometry size plays an important role since narrow slits limit the thermal motion of water molecules in these nano slits.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

Financial support from the National Key Research and Development Program of China(2017YFC0403902)and the Jiangsu Natural Science Foundations (BK20190085).We thank the High Performance Computing Center of Nanjing Tech University for supporting the computational resources.We are also grateful to the Program of Excellent Innovation Teams of Jiangsu Higher Education Institutions and the Project of Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).

Supplementary Material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.cjche.2020.10.028.