Hydrogen bond analysis of confined water in mesoporous silica using the reactive force field

TL;DR

This study uses molecular dynamics simulations with a reactive force field to analyze the structural and dynamical properties of water confined in mesoporous silica, focusing on hydrogen bonding behavior.

Contribution

It provides detailed insights into hydrogen bond dynamics and structural properties of confined water, validated against experimental data, using large-scale reactive force field simulations.

Findings

Hydrogen bonds involving silanol groups have longer lifetimes.

The diffusion coefficient of confined water is reduced compared to bulk water.

Hydrogen bond strength varies with radial position in the pore.

Abstract

The structural and dynamical properties of water confined in nanoporous silica with a pore diameter of 2.7 nm were investigated by performing large-scale molecular dynamics simulations using the reactive force field. The radial distribution function and diffusion coefficient of water were calculated, and the values at the center of the pore agreed well with experimental values for real water. In addition, the pore was divided into thin coaxial layers, and the average number of hydrogen bonds, hydrogen bond lifetime, and hydrogen bond strength were calculated as a function of the radial distance from the pore central axis. The analysis showed that hydrogen bonds involving silanol (Si-OH) have a longer lifetime, although the average number of hydrogen bonds per atom does not change from that at the pore center. The longer lifetime, as well as smaller diffusion coefficient, of these hydrogen bonds is attributed to their greater strength.

Equations

$$g_{12}(r)=\frac{⟨n_2(r)⟩}{4\pi r^2{\rho }_2\Delta r},$$

$$D=\underset{t\to \infty }{\lim }\frac{⟨|\vec{r}(t)-\vec{r}(0){|}^2⟩}{2dt},$$

$$r_{{\mathrm{O}}_1{\mathrm{O}}_2}\le 3.35\mathring{\mathrm{A}},$$

$$1.35\mathring{\mathrm{A}}\le r_{{\mathrm{O}}_1{\mathrm{H}}_2}\le 2.45\mathring{\mathrm{A}},$$

$$r_{{\mathrm{O}}_2{\mathrm{H}}_2}\le 1.35\mathring{\mathrm{A}},$$

$$\angle {\mathrm{H}}_2{\mathrm{O}}_2{\mathrm{O}}_1\le 30^{\circ },$$

$$C_x(t)={⟨\frac{{\sum }_{ij}{s}_{ij}(t+t_0)s_{ij}(t_0)}{{\sum }_{ij}{s}_{ij}^2(t_0)}⟩}_{t_0},$$

$${\tau }_x={\int }_0^{\infty }{C}_x(t)dt.$$

$$C_x(t)=A_1\exp \left(-\frac{t}{{\tau }_1}\right)+A_2\exp \left(-\frac{t}{{\tau }_2}\right)+(1-A_1-A_2)\exp \left(-\frac{t}{{\tau }_3}\right).$$

$${\tau }_{\mathrm{on}}={\int }_0^{\infty }t{P}_{\mathrm{on}}(t)dt.$$

$$vu=1.55-1.06r_{\mathrm{avg}}+0.186r_{\mathrm{avg}}^2,$$

$$r_{avg}=\frac{{\int }_{r=1.35}^{2.45}r{g}_{\mathrm{OH}}(r)r^{2}dr}{{\int }_{r=1.35}^{2.45}{g}_{\mathrm{OH}}(r)r^{2}dr}.$$

Full text

††articletype: ARTICLE TEMPLATE

1 Introduction

Understanding the structural and dynamical properties of water in confinement is essential to many scientific fields and technological applications, such as permeation in the ion channels of a biological membrane, capacitance of an electrical double layer in fuel cells, and controlled drug release [1, 2, 3, 4]. Mesoporous materials have attracted attention in recent years because their high specific surface area, large specific pore volume, and narrow pore size distribution allow a wide range of possible industrial applications. In particular, mesoporous silicate MCM-41, which has an amorphous structure, is one of the most widely used mesoporous material [5] and has potential application in drug delivery [6, 7, 8], in which water plays an important role [4]. In confined systems, surface properties have a significant impact on the physical properties of the confined substance. However, water behavior near an amorphous surface is much more difficult to investigate than that near a crystalline surface because the topology of an amorphous surface is ill-defined at the molecular level.

The behavior of confined water has been widely studied by both theoretical [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] and experimental [26, 27, 28, 29, 30, 31] techniques in the past few decades. These works have indicated that the properties of confined water may strongly differ from that of bulk water depending on surface properties and size of confinement space. It is thus important to properly evaluate the influence of confinement on water.

In recent years, the dynamical structure factor of water in mesoporous silica has been measured by quasi-elastic neutron scattering (QENS), and the self-diffusion coefficients ($D$) of water in fast and slow modes were evaluated by dividing the space into two regions. Yamada et al. reported that $D$ significantly differs between the vicinity of the interface and pore center [32]. However, the method of dividing the space is empirical, and its relationship with the microscopic atomic structure and dynamics is not sufficiently confirmed. This question is expected to be resolved through the analysis of trajectories from molecular dynamics (MD) simulations.

In MD simulations, the choice of force field is very important in reasonably estimating the interfacial effect. Recently, Bourg and Steefel reproduced the experimental diffusion constant obtained by QENS through MD simulations using the Clay force field (CLAYFF) and extended simple point charge (SPC/E) water model [17]. They noted that the use of an ab initio method  [33, 34, 35] or the reactive force field (ReaxFF) in MD simulations allows a good description of proton transfer reactions. van Duin and co-workers developed the ReaxFF package [36], which has been applied in various research areas [37]. In recent years, MD simulations of the silica/water interface have been performed using ReaxFF [38, 39, 40, 41, 42]. This method has a lower computational cost than quantum-mechanics-level calculations such as ab initio MD simulations. A drawback of MD simulations using classical force fields is that an appropriate chemical structure must be created before starting the simulation, and this does not change during the simulation because the chemical bond is fixed. This problem is particularly conspicuous when dealing with interfacial systems. Using ReaxFF, the physical properties of confined water can be obtained with high accuracy while avoiding the problem of interaction between water and silica.

In this study, we performed MD simulations of water in mesoporous silica using ReaxFF and analyzed water dynamics to clarify the influence of the surface on the physical properties of water at the atomic and molecular levels. Specifically, the hydrogen bond (H-bond) dynamics in the inhomogeneous system composed of amorphous silica (a-SiO2) and water was predicted. To obtain sufficient statistics, we used a large system since the pore was analyzed by dividing it into thin layers.

The remainder of this paper is organized as follows. The next section provides an overview of ReaxFF potentials and details of the preparation of the silica-water system. Section 3 presents the structural properties and diffusion of water in nanoporous silica and discusses the results of the analyses of H-bond configurations and dynamics. The final section summarizes the main findings and consequent conclusions.

2 Methods

2.1 ReaxFF potentials

To account for bond breakage and formation, we performed MD simulations using the ReaxFF potential [36]. In ReaxFF, all bond orders are calculated directly from interatomic distances every time step. The calculation of bond-order-based terms in the potential allows for bond breakage and formation. Note that H-bonds are explicitly defined in the ReaxFF potential. Additionally, the charge of each atom is dynamically assigned by a charge equilibration method [43]. Detailed explanation of ReaxFF can be found in previous articles [36, 44, 45, 37].

We used the ReaxFF parameter set for silica-water system developed by Yeon and van Duin in 2015 to simulate hydrolysis reactions at the SiO2/water interface [40]. It has been used to study interactions between water and nanoporous silica and has yielded results that agreed well with ab initio MD data for reaction mechanisms, hydroxylation rates, defect concentrations, and activation energies [41]. Thus, we chose this parameter set to study water dynamics in mesoporous silica. However, this parameter set was also reported to yield an unrealistically high fraction of pentacoordinated Si atoms in glassy silica [46]. Yu et al. reported that both reactivity and hydrophilicity are influenced by the atomic topology of the surface of the glassy silica [46]. Therefore, we also used the ReaxFF parameter set developed by Pitman et al. [47, 38], which provided an atomic structure of glassy silica that agreed well with the experimental structure [48, 49]. ReaxFF also allows O–H bond extension/contraction and changes in the H–O–H angle of water during MD simulations, unlike traditional force fields for water that treat the water molecule as rigid.

2.2 Preparation of the silica-water system

The system composed of water in nanoporous silica was created as follows. First, bulk a-SiO2 was prepared at 3000 K and 1 atm in the NPzT ensemble (i.e., constant number of atoms, pressure along the z-direction with fixed xy-area, and temperature) using the Morse potential [50]. The z-direction is parallel to the axis of a cylindrical pore that is subsequently created. After switching to ReaxFF MD simulation, the system was relaxed for 200 ps at 3000 K and 1 atm in the NPzT ensemble, then quenched to 300 K at a cooling rate of 15 K/ps in the NVT ensemble. After a 300-ps relaxation at 300 K and 1 atm in the NPT ensemble, 4 cylindrical pores with a diameter of 2.7 nm were created by removing atoms. Furthermore, a few oxygen atoms at the pore surface were removed to obtain a Si/O ratio of 1:2. Water was then placed inside the pore, and a 20-ps MD simulation was performed only for water using the TIP3P (transferable intermolecular potential with 3 points) model [51] while silica was immobilized. The system consists of an 8.9 $\times$ 7.7 $\times$ 45.1 nm rectangular box with periodic boundary conditions employed in all directions. Finally, ReaxFF MD simulation of the silica-water system was conducted. Silanols (Si-OH) rapidly formed at the interface through chemical reactions between silica and water. After a 500-ps equilibration run at 300 K and 1 atm, a 50-ps production run was performed for analyses. All simulations were conducted using the Nosé-Hoover thermostat and barostat using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software [52]. The ReaxFF MD simulations were performed using a time step of 0.25 fs through the USER-REAXC package of LAMMPS.

3 Results and Discussions

3.1 Structural properties

Figure 1 shows the average mass density of the silica-water system at 300 K, derived from MD simulations using Yeon’s [40] and Pitman’s [47] ReaxFF parameter sets, as a function of the radial distance from the pore central axis (R). It was calculated by dividing the cylindrical pore space into hollow cylinders with a 0.5-Å thickness, determining the density of each layer, and taking the time average over the 20-ps simulation. The approximate pore radius is also shown in Fig. 1. It was calculated by slicing the cylindrical pore space along the $z$-axis into thin disks with a 0.2-Å thickness, determining the $R$ of the innermost Si atom ($R_{\rm Si-in}$) in each disk, and taking the average of all values over the 20-ps time period. The $R_{\rm Si-in}$ values, calculated using Yeon’s and Pitman’s parameter sets, are 13.24 and 12.85 Å, respectively. The density distribution of Si, which represents the roughness of the amorphous surface, extends to $R=$ 11 (Yeon) and 10 Å (Pitman). Muroyama et al measured the surface roughness of MCM-41 by using powder X-ray diffraction and concluded that the surface roughness on the silica wall is 2 Å for the as-synthesized hexagonal pore and 1 Å for the calcined circular pore [53]. They defined the surface roughness as the distance from the boundary between the pores and the wall, which is determined using a threshold density ratio within the pore to that of the wall, to a point of a half value of the boundary density between the pores and the wall. We determined the surface roughness according to their definition. The boundary in this work is defined as the point where the density of silicon reach the bulk value and is located at $R=$13.8 Å and 13.6 Å for Yeon’s and Pitman’s parameter sets. The points of a half value of the boundary density are $R=$12.7 Å and 12.4 Å respectively, and thus the surface roughnesses are 1.1 Å and 1.2 Å for Yeon’s and Pitman’s parameter sets. These values are coincident with the experimental value of 1 Å for the calcined circular pore.

The densities of bulk a-SiO2 (averaged over $R=$ 17–20 Å), calculated using Yeon’s and Pitman’s parameter sets, are 2.08 and 2.21 g/cm3, respectively. These values coincide with those obtained in previous studies: 2.10 [42] and 2.18 [46] g/cm3 at 300 K using Yeon’s and Pitman’s parameter sets, respectively. The density of a-SiO2 calculated using Pitman’s parameter set is higher than that calculated using Yeon’s parameter set and closer to the experimental value of 2.2 g/cm3 [54]. Given the same number of water molecules confined in the pores, the two models yield different densities for water; the density in Yeon’s model is lower than that in Pitman’s model. Thus, in Yeon’s and Pitman’s models, about 8364 and 7724 water molecules are confined in each pore for a total of 33457 and 30896 water molecules (since there are four pores in each system), respectively, to make the water densities at the pore center almost equal. The average densities of water at $R<$5 Å are 0.993 and 1.006 g/cm3 in Yeon’s and Pitman’s models, respectively. Note that the density of water at 300 K and 1 atm is 0.997 g/cm3.

The radial distribution function (RDF) between two atoms is defined as

[TABLE]

where $r$ is the distance between a pair of type 1 and 2 atoms, $n_{2}(r)$ is the number of type 2 atoms in the shell between $r$ and $r+\Delta r$ around type 1 atoms and $\rho_{2}$ is the density of type 2 atoms. For example, in the O–H RDF, the type 1 and 2 atoms are O and H, respectively. The term in angular brackets represents the average over all atoms of interest and different starting times of the RDF calculation. To calculate the RDFs of water in the nanopore, atoms located at $R<$4 Å were taken as the type 1 atoms in Eq. 1, and the cut-off radius for the $g(r)$ calculation was 8.0 Å. Therefore, the RDF calculations were performed for atoms located at $R<$12 Å, excluding those near the interface. The RDFs calculated using Yeon’s and Pitman’s parameter sets are displayed in Fig. 2. Note that the first peak of the O–H RDF is the intramolecular O–H distance, while the second one represents the nearest intermolecular O–H distance. The locations of the peaks and minima in RDFs are listed in Table 3.1 with experimental values for comparison. The results in bulk water using classical water models are also listed as reference. The RDFs calculated using Yeon’s parameter set agree well with experimental values, while the water structure obtained using Pitman’s parameter set deviates from that of real water. The locations of the first minimum in the O–O and O–H RDFs and second minimum in the O–H RDF were used for H-bond analysis, as described in Section 3.3.

3.2 Water diffusion in nanoporous silica

The D of atoms can be calculated from the long-time behavior of the mean square displacement (MSD) using the Einstein relation:

[TABLE]

where d is the dimension and the term in angular brackets represents the average over all atoms of interest and different starting times of the $D$ calculation. Applying this equation to our system, however, is nontrivial because $D$ varies with the distance from the interface. On the basis of previous studies on water confined in nanospace [59, 17], we divided the pore space into hollow cylinders with a 0.2-Å thickness and calculated the D in each layer using the slope of MSD from t = 5 ps to t = 10 ps in the z-direction. For this layer thickness, this time scale was found to be sufficiently long for the molecular motion of water to be in the diffusive regime, but sufficiently short for individual water molecules not to probe regions with very different D [17]. The MSDs in the z-direction were calculated for oxygen atoms in each region over the 20-ps time period. The D along the z-direction ($D\rm_{z}$) was calculated from the MSD data using Eq. 2with $d=$1 and is shown in Fig. 3 as a function of R. The D of bulk water is 2.3 $\times$ 1011 Å2/s at 298 K based on several experiments [60]. Bourg and Steefel evaluated the D of real water at 300 K as 2.41 $\times$ 1011 Å2/s by interpolation from the D at 298 and 309 K using the Arrhenius relation [17]. Using Yeon’s parameter set, the $D\rm_{z}$ averaged over values at $R<$6 Å is 2.41 $\times$ 1011 Å2/s, which agrees with the experimental value for bulk water. On the other hand, using Pitman’s parameter set, the $D\rm_{z}$ is about 0.85 $\times$ 1011 Å2/s in the inner area of the pore. This value is much lower than that of bulk water, but corresponds well with the D of water in the clay-zeolite composite [47] used to develop the parameter set.

Yamada et al. examined the water dynamics in mesoporous silica with a 2.7-nm pore using QENS [32]. The water dynamics was divided into two modes, slow and fast, and the D of the latter corresponds to that of bulk water. Because the fast-mode water is presumably located near the pore center, water in the center of the 2.7-nm pore of mesoporous silica behaves like bulk water. This agreement with the results obtained using Yeon’s parameter set indicates that the parameter set is suitable for the investigation of the dynamics of water confined in mesoporous silica. In addition, the RDFs obtained using Yeon’s parameter set agrees better with the experimental values than those obtained using Pitman’s parameter set. Thus, the following sections only discuss the results calculated using Yeon’s parameter set. Note that Yu et al. demonstrated that both reactivity and hydrophilicity are determined by the atomic topology of the surface and Pitman’s parameter set provides a better description of the glassy silica structure [46]. Thus, it would be interesting to perform simulations using Pitman’s parameter set for SiO2 and other ReaxFF parameters for water.

3.3 Hydrogen-bond configurations and dynamics

The locations of the first minimum in the O–O RDF and second minimum in the O–H RDF were generally used for H-bond analysis. Note that the first and second peaks in the O–H RDF in Fig. 2 represent the nearest intramolecular and intermolecular O–H distances, respectively. The location of the first minimum in the O–H RDF was additionally used as a criterion for H-bonding since ReaxFF is composed of atom-based potentials. Thus, the H-bond O${}_{1}\cdots$H2–O2 is defined using the following criteria:

[TABLE]

where O1 is the H-bond acceptor and H2–O2 is the H-bond donor. The angular cut-off is widely accepted in the literature [61, 62, 13, 63, 64].

First, the average number of H-bonds per oxygen atom ($\langle n_{\rm HB}\rangle$) was calculated and is shown in Fig. 4 as a function of R. The calculations were performed by dividing the pore space into hollow cylinders with a 0.2-Å thickness. In bulk water, the average number of H-bonds per water molecule is between 2 and 4 depending on the experimental technique [65]. $\langle n_{\rm HB}\rangle$ in the center of the pore is about 3.5 at 300 K, which is a reasonable value for bulk-like water. It is almost constant at $R<$11 Å and decreases at $R>$11 Å. It indicates that the H-bond network becomes less structured at $R>$11 Å.

The H-bond lifetime can be evaluated through the following autocorrelation function [66]

[TABLE]

where the subscript $x$ represents the two different definitions of lifetime: intermittent (I) and continuous (C). In an intermittent lifetime, $s_{ij}$ is equal to 1 if a pair of atoms, i and j, are bonded at time t and 0 if the bond is absent. Note that only pairs of atoms having $s_{ij}$ = 1 at t = 0 were included in the calculation. In a continuous lifetime, $s_{ij}$ is allowed a single transition from 1 to 0 when the bond is first broken, but is not allowed to return to 1 if the same bond is reformed. Namely, $C_{\rm I}$(t) and $C_{\rm C}$(t) focus on the elapsed time until the final and first breakage of the bond, respectively. The term in angular brackets refers to the average over different starting times $t_{0}$ during the production run. The calculations were performed by dividing the pore space into hollow cylinders with a 1.0-Å thickness and only for acceptor oxygens remaining in each layer during the simulation.

The average H-bond lifetime $\tau_{\rm HB}$ can be evaluated by integrating $C_{x}$(t) using the following equation [63, 67],

[TABLE]

To obtain the H-bond lifetime, the $C_{x}$(t) calculated directly from the MD trajectories was fitted to the sum of three exponential functions and integrated analytically.

[TABLE]

The calculated $C_{x}$(t) for several layers are shown in Fig. 5 with the fitting functions. For the center of the pore, both $C_{\rm{I}}$(t) and $C_{\rm{C}}$(t) can be fitted to the sum of two exponential functions. However, for the interface, the sum of three exponential functions show better fitting than the sum of two exponentials to both $C_{\rm I}$(t) and $C_{\rm C}$(t). Thus, $C_{x}$(t) was fitted to three exponentials for all layers.

The H-bond lifetimes $\tau_{\rm I}$ and $\tau_{\rm C}$ were calculated through Eq. 5 using these fitting functions and are displayed in Fig. 6. They are almost constant at $R<$8 Å, but increase at $R>$8 Å, which corresponds to the trend in $D\rm_{z}$ (Fig. 3). The $\tau_{\rm I}$ in the center of the pore is $\sim$0.4 ps, which is shorter than those obtained in other studies [63, 61, 62, 68]. The $\tau_{\rm I}$ in bulk water is reported to be 2 to 9 ps depending on the water model. One of the reasons is that $\tau_{\rm I}$ was calculated only for atoms that remain in each thin layer. Atoms that moved to the adjacent layer were not counted in the $C_{x}$(t) calculation if they reformed the same H-bond. In fact, the $\tau_{\rm I}$ calculated over atoms located at $R<$6 Å was about 3.2 ps, which is comparable with literature values [63, 61, 62, 68]. Another reason is the short simulation time. H-bonds that reformed over the 20-ps simulation were not included in the $\tau_{\rm I}$ calculation, although they would raise $\tau_{\rm I}$. From the $\tau_{x}$ calculated in this work, however, the dynamical behavior of confined water can be captured. Figure 4 shows that $\langle n_{\rm HB}\rangle$ is almost constant up to $R$ = 11 Å, although the H-bond dynamics begin to change from $R\sim$ 8 Å (Fig. 6). This is because a few H-bond pairs change from water to silanol. The density profiles of four types of oxygen, H2O, –Si–O–H, –Si–O*-, and –Si–O–Si–, are displayed as a function of $R$ in Fig. 7. Other types of oxygen, such as OH-* or H3O*+*, have much smaller densities and thus are not shown in the figure. The density of silanol oxygen (OSi-H) starts to increase from $R\sim$ 10 Å and peaks at $R\sim$ 12 Å. Thus, oxygen in water (Ow) within 8 Å$<R<$10 Å can form H-bonds with silanols because the H-bond distance is 1.35–2.45 Å. Oxygen atoms within 8 Å$<R<$10 Å, where H-bonds include silanol, have an average of 3.5 H-bonds (Fig. 3) like those near the center of the pore, where H-bonds are formed only by water. The decrease in the $D\rm_{z}$ of oxygen with $R$ at $R>$8 Å (Fig. 3) implies that the H-bond network including silanol slows down water dynamics. The influence of interfacial silanol on the H-bond lifetime was also reported for confined glycerol in a silica nanopore [69].

$\tau_{\rm C}$ slightly increases from $R\sim$ 8 Å and peaks at around $R$ = 12–13 Å. To see the first breakage time of H-bond more clearly, we calculated another time constant, $\tau_{\rm{on}}$. It can be evaluated from $P_{\rm{on}}(t)$ [61], which is defined as the distribution of time from the formation of a H-bond to its destruction. $\tau_{\rm{on}}$ is calculated using the following equation,

[TABLE]

Figure 8 shows $\tau_{\rm on}$, which was calculated over a 50-ps time period. At $R<$8 Å, $\tau_{\rm on}$ is almost constant at $\sim$0.34 ps, which is comparable to that of bulk water in Ref. [61]. $\tau_{\rm{on}}$ also increases with $R$ in the range of 8 Å$<R<$11 Å, but decreases with $R$ at $R>$11 Å. This behavior is similar to that of $\tau_{\rm C}$, although the peak position is slightly different. The behavior of $\tau_{\rm on}$ at $R>$8 Å is elucidated by considering the valence unit per hydrogen bond (vu), which was introduced as an index of H-bond strength by Machesky et al. [70].

$vu$ is obtained from

[TABLE]

where the average H-bond distance $r_{\rm avg}$ is calculated as

[TABLE]

$g_{\rm OH}(r)$ in Eq. 9 is calculated for oxygen atoms in each layer with a width of 1.0 Å. The integration range of $r$ = 1.35–2.45 Å corresponds to the second peak of $g_{\rm OH}$, namely, the nearest intermolecular O–H distance. Figure 9 shows $r_{\rm avg}$ and $vu$ as a function of $R$. $vu$ at $R<$8 Å is about 0.19, which agrees with the value for a normal H-bond in bulk water at 298 K (0.2) [71]. The plots in Figs. 8 and 9 suggest that the behavior of $vu$ against $R$ is similar to that of $\tau_{\rm on}$. $vu$ increases between 8 and 11 Å and reflects the decrease in $r_{avg}$ (Fig. 9), which indicates that H-bonds including silanol (8 Å$<R<$11 Å) are stronger than that of bulk water. Stronger H-bonds lead to a longer H-bond lifetime (Fig. 8) and slower diffusion (Fig. 3); thus, H-bonds involving silanol are long lived.

On the other hand, $D\rm_{z}$ values at $R>$11 Å are quite small compared with those at $R<$11 Å (Fig. 3), which indicates that molecules located at $R>$11 Å are almost immobile. Moreover, at $R>$11 Å, $\tau_{\rm{I}}$ (i.e., the final breakage time) increases with $R$ (Fig. 6), while $\tau_{\rm on}$ (i.e., the first breakage time) decreases with $R$ (Fig. 8). These observations suggest that the H-bond acceptor (O) at $R>$11 Å frequently and repeatedly forms and destroys the H-bond with the same H-bond donor (OH). This is caused by the less-structured H-bond network, as shown in Fig. 4. Fig. 1(a) shows that the density distribution of Si starts to increase at $R=$ 11 Å and continues until $R=$ 14 Å. The density distribution of O in SiO2 also increases within this region, as shown in Fig. 7. Namely, the roughness of a-SiO2 surface disrupts the formation of H-bond network.

Finally, we focused on the H-bond network and surface topology of water located at $R=$ 10 Å (Fig. 10(a)) and $R=$ 13 Å (Fig. 10(b)). Water at $R=$ 10 Å forms four H-bond networks, two each with water and silanol. The H-bond lengths with silanol are 1.77 and 1.80 Å, while those with water are 1.84 and 1.87 Å. The shorter (and stronger) H-bond length with silanol attracts water near the interface and leads to a slight increase in the density of Ow around $R=$ 8–10.5 Å (Fig. 7), which also contributes to the slower diffusion. Water molecules at $R=$ 13 Å  on the other hand, are located in a small hole on the rough silica surface and surrounded by three silanols. It is separated from the water H-bond network and form H-bonds only with silanols. The configuration of water at $R>$11 Å can be estimated from Fig. 11. It shows the distribution of water atoms $P_{\rm Ow}$ and $P_{\rm Hw}$, normalized so that the sum of each probability becomes 1. At $R>$11 Å, $P_{\rm Hw}$ is slightly larger than $P_{\rm Ow}$. It indicates that water hydrogen is closer to the wall than water oxygen in this region. Water configuration at $R>$11 Å is thus presumed to be captured in a silica hall and with silanols as shown in Fig. 10(b). Because of this strong H-bond network with silanols and steric hindrance by surrounding silica atoms, they are almost immobile.

4 Summary and Conclusion

We conducted large-scale MD simulations of confined water in nanoporous silica with a pore diameter of 2.7 nm using ReaxFF. We compared Yeon’s and Pitman’s parameter sets for our system. We found that the former is suitable for the analyses of water structure and dynamics in a 2.7-nm silicate pore, although the density of a-SiO2 was better using Pitman’s parameter set. We calculated the radial distribution function and diffusion coefficient of water using Yeon’s parameter set, and the values in the center of the pore agreed well with those of real water obtained from several experimental studies. In addition, the diffusion coefficient of water essentially corresponded to previous simulation studies using the classical MD force fields, CLAYFF and SPC/E water [17]. Thus, these force fields can capture the water diffusion in silica nanopores well.

We divided the pore into hollow cylindrical layers and calculated the diffusion coefficient, average number of H-bonds, H-bond lifetime, and H-bond strength as a function of the radial distance from the pore central axis ($R$). Based on H-bond structure and dynamics, the pore can be divided into three regions: $R<$8 Å, 8 Å$<R<$11 Å, and $R>$11 Å. At $R<$8 Å, water behaves like bulk water. Within 8 Å$<R<$11 Å, H-bonds involving silanol have a longer H-bond lifetime, although the average number of H-bonds per atom does not change from that in the pore center. H-bonds involving silanol become stronger, resulting in a longer H-bond lifetime and smaller diffusion coefficient. Finally, at $R>$11 Å, the H-bond network becomes less structured owing to fewer surrounding water. This is caused by steric hindrance by interfacial Si and O in SiO2 and surrounding silanols, which arise from the amorphous surface. Thus, the elapsed time for the final H-bond breakage ($\tau_{\rm I}$) becomes longer and the diffusion coefficient becomes smaller.

Water behavior near the interface depends on the roughness of the surface and fraction of silanols present, which will be investigated in our future work. The MD simulation with ReaxFF provides highly reliable results, as evidenced by the good agreement with experimental values. Thus, this method is also expected to accurately predict water behavior near the interface, which is difficult to access in experiments. In addition, H-bonds involving silanol become stronger than that of bulk water because of the shorter O–H distance, although the detailed structure of the H-bond network is unclear. In future work, we will also examine how the network structure causes slower dynamics in detail.

Acknowledgement(s)

This research used the computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Project ID: hp170021), of the Supercomputer Center of the Institute for Solid State Physics provided by The University of Tokyo, and of the Research Center for Computational Science at Okazaki, Japan. The authors are partially supported by the Research Center for Computational Science, National Institute of Natural Sciences. This work was partially supported by the Japan Society for the Promotion of Science KAKENHI grant nos. JP15K05244, JP18H04494 and JP19K05209.