Molecular Simulation of the Phase Behavior of Water Confined in

Publication Date (Web): May 16, 2007. Copyright ..... Sophisticated and Spontaneous Template-Free Organization of Silica Nanoparticles During Storage...
0 downloads 0 Views 659KB Size
Download PDF
7938

J. Phys. Chem. C 2007, 111, 7938-7946

Molecular Simulation of the Phase Behavior of Water Confined in Silica Nanopores Katsuhiro Shirono† and Hirofumi Daiguji* Institute of EnVironmental Studies, Graduate School of Frontier Sciences, The UniVersity of Tokyo, 5-1-5 Kashiwanoha, Kashiwa-shi, Chiba 277-8563, Japan ReceiVed: NoVember 8, 2006; In Final Form: March 20, 2007

Molecular dynamics simulations and grand canonical Monte Carlo simulations of water molecules in 1.04, 1.96, and 2.88 nm diameter silica pores (pores S, M, and L, respectively) were conducted at 300 K to reveal possible phase states of water in each pore and clarify the effect of pore diameter and hydration level on the structural and dynamical properties of water molecules confined in nanopores. Three types of phases appear in nanopores. In the first phase, which appears in pores S, M, and L, gas-phase water adsorption at the pore surface is below monolayer coverage. In the second phase, which appears in pores S and L, there is a condensed water monolayer at the pore surface. In the third phase, which appears in pores M and L, the pore is completely filled with water. Water molecules interact with silanol groups mainly via hydrogen bonds at low hydration when the number of water molecules is much smaller than that of silanol groups. With increasing number of water molecules, the number ratio of water molecules that adsorb on silanol groups via non-hydrogen-bonding interactions increases. Near full hydration, the translational mobility of water in the first adsorption layer is much smaller than that of bulk liquid water in all three pores, while those in the pore center are about 30, 80, and 100% of its bulk liquid value in pores S, M, and L, respectively.

1. Introduction Microporous and mesoporous silica are among the most promising materials as adsorbents of water. Techniques to synthesize ordered mesoporous silica have been developed since the synthesis of MCM-41 was first reported in 1992,1 and recent techniques make it possible to synthesize porous materials that have various desired properties.2 Determining the possible phase states of water in equilibrium (i.e., number of coexisting phases (one, two, or three in the case of a triple point)) is crucial for understanding the adsorption and diffusion properties inside nanopores and is essential for the design of adsorbents of water. Because water in confined geometries is ubiquitous in nature and widely used in various industrial processes, understanding the structural and dynamical properties of confined water, which differ significantly from the properties of the bulk liquid water, is of great fundamental and practical importance. Clarification of the possible phase states of water is also crucial for interpretation of the structural and dynamical properties of water. To clarify the possible phase states of water in nanopores, Brovchenko et al.3 reported the coexistence curves of water in nanopores of different shapes, sizes, and strengths of watersubstrate interactions by using Gibbs ensemble Monte Carlo simulations. The results showed five different two-phase coexisting regions: two different “bulklike” liquid-vapor phase transitions (with the pore wall covered or not covered with two water layers) and three quasi two-dimensional liquid-vapor surface transitions (first layering transition, second layering transition, and prewetting). Because the hydrophilicity of pore surfaces was controlled by changing the Lennard-Jones (LJ) potential parameters, the effect of the local structure of water molecules on the pore surface (i.e., the local structure of water * To whom correspondence should be addressed. Email: daiguji@ k.u-tokyo.ac.jp. Tel: +81-4-7136-4658. Fax: +81-4-7136-4659. † Current address: National Institute of Advanced Industrial Science and Technology (AIST), 1-2-1 Namiki, Tsukuba 305-8564, Japan.

molecules around silanol groups on the silica surface) on the possible phase states of water was not clarified in that study. To clarify the adsorption properties of water in silica pores, Puibasset et al.4 calculated the adsorption/desorption isotherms of water in Vycor-like mesoporous silica glass (average pore diameter of about 3.6 nm) by using GCMC simulations. Their results showed a hysteresis loop in the adsorption/desorption isotherm at 300 K, but also showed a reversible curve at 650 K. Inagaki et al.5 determined experimentally the water adsorption isotherms of 1.4-2.7 nm diameter FSM-16 at 25 °C and showed that the hysteresis loop in water adsorption isotherms becomes narrow with decreasing pore size and disappears at a pore diameter of 1.4 nm. Although similar behavior of disappearing hysteresis has also been observed for other adsorbates such as argon and oxygen,6 the pore size at which the hysteresis disappears for water is smaller than that for other adsorbates.5 The water adsorption isotherms are sensitive to the pore diameter, especially in the range of 1-3 nm; however, few theoretical studies have been done on water vapor adsorption into nanopores within this range of diameter. To understand the diffusion of water confined in silica pores, both the local structure and the dynamics of water around silanol groups as well as the phase states of water must be clarified. The structural and dynamical properties of water molecules confined in silica pores have been investigated both theoretically and experimentally. Gallo and Rovere7 reviewed molecular dynamics (MD) of water confined in silica pores as a function of both temperature and hydration level. They reported that the interaction of water in the film close to the substrate with the silica atoms induces a strong distortion of the hydrogen bond network and slows down water dynamics. Smirnov et al.8 conducted X-ray diffraction measurements of hydrated MCM41 at 223-298 K to investigate the structural properties of absorbed water on monolayer and capillary-condensed conditions. Their results showed that the structural properties of water

10.1021/jp067380g CCC: $37.00 © 2007 American Chemical Society Published on Web 05/16/2007

Phase Behavior of Water Confined in Silica Nanopores

J. Phys. Chem. C, Vol. 111, No. 22, 2007 7939

Figure 1. Molecular model of 1.04 nm diameter (pore S) (left), 1.96 nm diameter (pore M) (center), and 2.88 nm diameter (pore L) (right) silica pores after an equilibration calculation time of 100 ps. Cross-section of each model pore along the radial axis (r-axis) and pore axis (z-axis) is 0.6 nm thick.

were independent of temperature for their monolayer sample at relative humidity of 0.3. They concluded that this independence was due to strong hydrogen bonds between water molecules and silanol groups on the surface of MCM-41. Crupi et al.9 reported a spectroscopic study of water dynamics in bulk and in sol-gel porous glass with 2.6 nm diameter interconnected cylindrical pores at 20 °C and various hydration levels using Rayleigh wing, Raman, and quasi-elastic and inelastic neutron scattering techniques. The Rayleigh wing and Raman spectroscopy showed that the tetrahedral networks of water molecules were broken in the pores. The inelastic neutron scattering showed that the hindered translational mode was flattened and had attenuated with respect to the bulk water. Crupi et al.9 concluded that surface water loses its unique bulk properties and originates a new structural environment due to surface interactions that are driven by the hydrogen bonds. Takahara et al.10,11 reported the quasi-elastic neutron scattering study on water dynamics in MCM-41 at 238-294 K to investigate the effect of pore size on water dynamics. Their measured translational diffusion coefficient of confined water was smaller than that of bulk water and decreased with decreasing pore diameter over the entire range of temperature. Several other studies12-14 reported a smaller diffusion coefficient or longer relaxation time of translational motion of water molecules in silica pores due to strong interactions between water molecules and silanol groups, but only a few studies have reported on the molecular interactions in detail. Because the local structure and dynamics of water around silanol groups affect the diffusion as well as adsorption of water in silica pores, the modeling of porous silica is crucial for clarifying these properties in molecular simulations. In this study, canonical ensemble (NVT) MD simulations and grand canonical Monte Carlo (GCMC) simulations of water molecules in 1.04, 1.96, and 2.88 nm diameter silica pores at 300 K were conducted to reveal the effect of pore size on the possible phase states and on the structural and dynamical properties of water. In this report, Section 2 provides details of the MD and GCMC simulations, including the modeling of silica pores and the choice of potential functions and parameters. Section 3 describes the calculation results of the phase transitions and the structural and dynamical properties of water molecules. Section 4 provides a summary and states the conclusions.

2. Methodology 2.1. Modeling of Silica Pores. Realistic modeling of silica pores with surface roughness and pore size distribution has recently been reported.15,16 In our current study, to clarify the effect of pore size we modeled silica pores for a crystal structure of R-quartz. The following five-step procedure was used in this modeling. In the first step, 9, 10, and 8 trigonal unit cells of R-quartz (Si3O6), whose structure has been clarified in detail by Wright and Lehmann,17 were arrayed in the [100], [010], and [001] directions, respectively. This R-quartz crystal was transformed from a rhombohedron to a rectangular parallelepiped with a volume of 4.422 × 4.255 × 4.324 nm3 by relocating the positions of unit cells. The pore axis (z-axis) was selected as the direction parallel to the [001] direction. In step 2, Si atoms within a radius r0 from the z-axis were removed. In step 3, O atoms were also removed that were covalently bonded with the two adjacent removed Si atoms. In step 4, to prevent dangling bonds H atoms were added to O atoms bonded with only one Si atom. In step 5, to minimize stress in the model pores NVT-MD simulations were conducted with all atoms mobile at 300 K. In our calculation, we used the Tsuneyuki potential model18 in which SiOH groups are assumed to be rigid bodies. The potential functions and parameters for O and Si atoms of the SiOH groups were set to be the same as those for SiO2. To ensure that the overall charge of the entire system was neutral, the charge of H atoms of SiOH groups, QH(SiOH), was set to +0.6e, and only the Coulomb interactions were considered in the interactions between H atoms of SiOH groups and all other atoms. Figure 1 shows the final structure of silica pores obtained after an equilibration calculation time of 100 ps. The composition of the model pores depends on both the pore radius r0 and on the position of the z-axis. In this study, the position of the z-axis was chosen as the center of a six-membered cyclosilicate ring, and r0 was assumed to be 0.52, 0.98, or 1.44 nm. The silica pore models in which r0 ) 0.52, 0.98, and 1.44 nm were labeled as pores S, M, and L, respectively, and the composition of pores S, M, and L was Si2064O4176H96, Si1824O3744H192, and Si3024O1440H288, respectively. All silanol groups on the surface

7940 J. Phys. Chem. C, Vol. 111, No. 22, 2007

Shirono and Daiguji

Figure 2. Calculated potential energy surface for water molecules and for the cluster containing one silanol group (OH-Si-(O-SiH3)3), Uad [kJ/mol], along the direction of optimized hydrogen bond for the PD structure (top) and PA structure (bottom). Optimization was done by using the DFT of B3LYP/6-31G(d) and Monte Carlo simulation at 30 K using the proposed potential energy model. r [Å] is the hydrogen bond length, and OW, HW, OS*, and HS* represent oxygen and hydrogen atoms of water molecules and silanol group, respectively.

the model pores were considered neither isolated nor germinal silanol groups, but vicinal silanol groups. (See ref 19 for the nomenclature of isolated, germinal, and vicinal silanol groups.) The final pore radii were smaller than those before equilibration and the pore structures were hexagonal-like. Coasne et al.20 reported that the difference in adsorption of the LJ particles inside smooth cylindrical and hexagonal pores is marginal. In our model, the average distance of a hydrogen bond between vicinal silanol groups was 1.57 Å, which is slightly smaller than the previously reported calculated value of 1.68∼1.71 Å for the surface of SiO2 (β-cristobalite) crystal by using ab initio MD simulation.21 2.2. Potential Functions and Parameters between WaterSilica Pores. The water molecules were assumed to interact via the extended simple point charge (SPC/E) model potential.22 The potential between water-silica pores involved two contributions: Coulomb interaction between all the atoms, and BornMayer-Huggins (BMH) interaction between the oxygen atoms of water (OW) and the oxygen atoms of silica/silanol groups (OS/OS*). In this study, van der Waals interactions were neglected between the hydrogen atoms of silanol groups, silicon atoms of silica/silanol groups (HS* and Si/Si*), and oxygen and hydrogen atoms of water molecules (OW and HW), namely, HS-OW, HS-HW, Si/Si*-OW, and Si/Si*-HW. Note that the potential functions and parameters for OS* and Si* are the same as those for OS and Si, respectively. The BMH potential function can be expressed in Tsuneyuki potential18 as:

U(rij) ) f0(bi + bj)exp[(ai + aj - rij)/(bi + bj)] - cicj/rij6

(1)

where f0 is a conversion factor in which f0 ) 1 (kcal Å-1 mol-1) ) 4.3384 × 10-2 (eV Å-1), and the parameters a and b are the size and hardness of ions, respectively, and have units of Å.

Figure 3. Chemical potential of water molecules, µ [kJ mol-1], in pore S (top), pore M (middle), and pore L (bottom) as a function of average water density, Fj [g cm-3], calculated by using MD + CIW method and GCMC. Fj is given by N mw/Veff, where N is the number of water molecules, mw [g] is the molecular mass of water and Veff [cm3] is effective pore volume defined as the volume of cylindrical pores in which water molecules were inserted. The diameter of cylindrical pores S, M, and L was 0.89, 1.81, and 2.77 nm, respectively. The number at each symbol expresses N.

This combination rule of BMH potential function was applied in this study by fitting the LJ potential function of the SPC/E water model to the BMH potential function and setting the fitted parameters as aOW ) 1.694 Å, bOW ) 0.1179 Å, and cOW ) 5.21 eV1/2 Å.3 To validate the accuracy of these potential functions and parameters, the position of water molecules adsorbed on a silanol group was compared to the position calculated using density functional theory (DFT). The structure around the silanol group of silica pores was modeled as OH-Si-(O-SiH3)3. The overall charge of this system was determined so as to be neutral. Calculations were done using Gaussian 98.23 The B3LYP hybrid exchange-correlation function and the 6-31G(d) basis sets were used in the optimization of the structures of the water molecules around the silanol group. The calculation revealed two optimized structures (schematics at right in Figure 2): a structure in which water molecules are adsorbed on the silanol group as a proton donor (PD) and a structure in which water molecules are adsorbed as a proton acceptor (PA). To validate the potential model used in this study, Monte Carlo simulations were conducted for a water-OH-Si-(OSiH3)3 system at 30 K. In this simulation, the interactions involving H atoms of SiH3 were neglected except for the Coulomb interaction, and the charges of these atoms, QH(SiH3), were determined such that the overall charge of the entire system was neutral (i.e., QH(SiH3) ) -0.6e). Only the water molecule was allowed to move while OH-Si-(O-SiH3)3 was fixed to the position obtained by optimization with DFT. For the PD structure, based on these simulations with the proposed potential

Phase Behavior of Water Confined in Silica Nanopores

J. Phys. Chem. C, Vol. 111, No. 22, 2007 7941

model the hydrogen bond length between the water-silanol group was 1.91 Å and the adsorption heat was -22.56 kJ/mol, compared with 1.93 Å and -25.12 kJ/mol calculated using DFT. For the PA structure, the calculated values were 1.81 Å and -45.99 kJ/mol, compared with 1.83 Å and -39.14 kJ/mol calculated using DFT. Figure 2 shows the potential surfaces along the direction of the hydrogen bond in the optimized positions calculated using DFT (B3LYP/6-31G(d)) and using Monte Carlo simulation at 30 K. The potential surface calculated using the proposed potential model agrees well with that calculated using DFT and the molecular orbital (MO) method (MP2/6-31G(d)). The calculations for the optimization and zero point energy in MO were conducted using DFT/6-31G(d). 2.3. Simulation Details. MD simulations with NVT ensemble and GCMC simulations were performed at 300 K. NVT MD simulations were done for 0, 12, 18, 24, 36, 48, and 54 water molecules in pore S, for 0, 24, 48, 96, 144, 192, 240, and 288 water molecules in pore M, and for 0, 24, 48, 96, 192, 288, 336, 384, 480, 576 and 672 water molecules in pore L. Initially, water molecules were placed randomly within the pores to a desired concentration. The temperature was controlled at 300 K by using a Nose´ thermostat.24,25 Periodic boundary conditions were applied to the pore model described in Section 2.1, and the electrostatic interactions were calculated using the Ewald method.26 Both the cutoff distance of van der Walls interaction and the real space term of the Coulomb interaction were 12.5Å. To minimize calculation time, all atoms except for water molecules and silanol groups were fixed at the initial position prepared by the procedure described in Section 2.1. Water molecules and silanol groups were allowed to equilibrate at 300 K for 200 ps. In all the simulations, the time step was 1 fs and the total run time was 1200 ps. Furthermore, the chemical potential of water molecules was calculated by the cavity insertion Widom (CIW) method27,28 using the data obtained by MD simulation. The chemical potential, µ, of each system was calculated as the summation of the ideal part of the chemical potential involving the density, µid, and the excess part of the chemical potential involving the interactions, µex

[

µid ) kT ln ln

N 3 Λ Veff

[ (

)(

2 π1/2 8π IAkT σ h2

1/2

h2 ∆U N+1 µex ) - kT ln exp kT

[ (〈 (

)(

8π2IBkT

1/2

)〉)

) ]]

8π2ICkT h2

]

Dz ) lim

tf∞

〈∆z2〉 2t

(4)

where is MSD along the z-axis and t is time. In this study, Dz was thus estimated by using the data of from t ) 200 to 800 ps. 3. Results

1/2

- ln〈Pcav〉

angles. For each hydration state, chemical potential calculations were conducted by using 1000 different atomic positions obtained every 1 ps. The calculation code was calibrated based on the chemical potential of bulk liquid water. The calculated chemical potential of SPC/E water at a density of 1.000 g/cm3 and 300 K was -52.7 kJ/mol. Whitley et al.29 reported that the calculated chemical potential of SPC/E water at the density of 0.977 g/cm3 and 298K was -43.5 kJ/mol when the ideal chemical potential was given by µid ) -kT ln(FΛ3). If the ideal chemical potential was given by eq 2, the corresponding value would be -53.2 kJ/mol, which agreed well with our calculated value of -52.7 kJ/mol. However the chemical potentials calculated by the CIW method normally include large calculation error. Jedlovszky and Mezei28 reported that the error ranges from kT to 2kT for bulk water. In this study, GCMC simulations were also conducted to verify the possible phase states of water in silica pores. A GCMC move consisted of one of the following: translation, rotation, insertion, and removal. Equal numbers of trials for these four moves were attempted and the number of trials was more than 1.5 × 104 for each hydration state after equilibration. For a translation move, the SMART algorithm26 was employed. In the SMART algorithm, the probability density function to chose a trial displacement of molecule i from state m to state n, δrnmi, was given by a Gaussian function: Rnmi ) (4Aπ)-3/2 exp(-(δrnmi - Afmi/kT)2/4A), where fmi is the force exerted on molecule i at state m and A is the adjustable parameter. The parameter A was modulated to obtain the acceptance ratio of 0.1. For a rotation move, the maximum rotational displacement in the Euler angle was set 0.1π. The local diffusion coefficient of water, D(r), was calculated by using the data of mean square displacements (MSD) of water from 2 to 10 ps. This method was used in ref 30. The diffusion coefficient along the z-axis, Dz, was calculated as

(2) (3)

where k is Boltzmann’s constant, T is temperature, N is number of water molecules, Veff is the effective volume for the adsorption of water molecules, Λ is the de Broglie wavelength of molecules, σ is the symmetry number of a molecule (σ ) 2 for SPC/E water molecule), IA, IB, and IC are the principal moments of inertia, h is Planck’s constant, ∆UN+1 is the energy of the inserted molecule at N, and Pcav is the probability of finding cavities from random insertion. For pores S, M, and L, water molecules were inserted into 0.89, 1.81, and 2.77 nm diameter cylindrical pores, respectively. The effective volume of each pore, Veff, was regarded as the volume of these cylindrical pores. When the distance between the oxygen atoms of inserted water and the nearest oxygen atom of water or silanol group was longer than 2.3 Å, the insertion location was regarded as a cavity. In the calculations, a water molecule was inserted into a cylindrical pore 200 000 times and the inserted water molecule was rotated around the oxygen atom at 10 different

3.1. Possible Phase States of Water in Silica Pores. Figure 3 shows the chemical potential versus average water density (µ-Fj) curves for pores S, M, and L calculated by the following two methods: (i) a combination of NVT MD simulations and chemical potential calculations by the CIW method and (ii) GCMC simulations. Fj is given by N‚mW/Veff, where N is the number of water molecules, mw is the molecular mass of water, and Veff is effective pore volume defined as the volume of cylindrical pores in which water molecules were inserted. In NVT MD simulations and chemical potential calculations by the CIW method, the calculated µ at the highest hydration levels (i.e., N ) 54, 288, and 672 for pores S, M, and L, respectively) is -55.9, -54.1, and -56.5 kJ/mol, respectively. These values are close to the calculated µ of bulk liquid water at 300 K (-52.7 kJ/mol). In Figure 3, the calculated µ does not increase monotonically with respect to N. If ∂µ/∂N < 0, the calculation system becomes unstable thermodynamically and the phase transition occurs over this regime. In pores S and M, there is only one regime where ∂µ/∂N < 0 over the entire range of N, whereas in pore L there are two regimes where ∂µ/∂N < 0. The possible phase states of water appear on both sides of the regime where ∂µ/∂N < 0. In GCMC simulations, the calculated systems

7942 J. Phys. Chem. C, Vol. 111, No. 22, 2007

Shirono and Daiguji

Figure 4. Possible phase states of water in pores S, M, and L. In pore S, (S-I) adsorption on silanol groups at N ) 12; (S-II) monolayer adsorption at N ) 54. In pore M, (M-I) adsorption on silanol groups at N ) 48; (M-II) capillary condensation at N ) 288. In pore L, (L-I) adsorption on silanol groups at N ) 48; (L-II) monolayer adsorption at N ) 288; (L-III) capillary condensation at N ) 672.

are always in the stable or metastable states, and the unstable hydration states do not appear. In Figure 3, the µ-Fj curves for pores S and M calculated by GCMC simulations showed only one phase transition, whereas those for pore L showed two phase transitions. The results of GCMC simulations verified the possible phase states of water in silica pores estimated by NVT MD simulations and chemical potential calculations by the CIW method. In more detail, the µ-Fj curve of GCMC did not agree well with that of MD + CIW method at low hydration in pore S and at high hydration in pores M and L. But these deviations did not affect the judgment of the possible phase states of water in pores S, M, and L. In the following sections, the structural and the dynamical properties were analyzed at N ) 12, 48, and 54 in pore S, N ) 24, 48, 240, and 288 in pore M, and N ) 24, 48, 288, 336, 624, and 672 in pore L. Figure 4 shows snapshots of the MD simulations for the possible phase states of water in pores S, M, and L. In pore S, there are two phase states of water. In the first phase (S-I), which appears at N ) 12, water molecules are adsorbed around silanol groups. In the second phase (S-II), which appears at N ) 48 and 54, water molecules form a monolayer over the entire pore surface. The first phase of pore M (M-I) corresponds to S-I, namely, water molecules are adsorbed around silanol groups, whereas in the second phase of pore M (M-II), the pore is filled with liquid water. In pore M, capillary condensation starts before the entire pore surface is covered with a water film, that is, no layer adsorption occurs. The M-I and M-II phases appear at N ) 24 and 48 and at N ) 240 and 288, respectively. In pore L, the first phase (L-I) corresponds to S-I and M-I. In the second phase (L-II), a water

monolayer is formed on the pore surface, and in the third phase (L-III), pores are filled with water. The L-I, L-II, and L-III phases appear at N ) 24 and 48, at N ) 288 and 336, and at N ) 624 and 672, respectively. The phase transition between L-I and L-II corresponds to the first layering phase transition. Brovchenko et al.3 concluded from their simulation study that this type of phase transition is observed only in hydrophilic pores. Figure 5 shows the µex profiles along the r-axis in pores S, M, and L. In pore L, the µex profile at dehydration (N ) 0) shows that the minimum µex is -63 kJ/mol around r ) 12 Å and that µex f0 when rf0. With increasing N, µex near the pore surface increases because the area fraction of the pore surface where no water molecules are present decreases. With increasing N, µex near the pore center decreases because water molecules located near the pore center interact with the adsorbed water molecules and thus are stabilized in terms of energy. If the µex profile at N ) 48 (L-I) is dislocated 2 Å toward the center of the pore, it coincides with the µex profile at N ) 288 (L-II). At N ) 672 (L-III), µex is constant along the r-axis except for the region near the pore surface. In pore M, the µex profile at dehydration (N ) 0) is almost the same as that for pore L at 5 Å < r < 14 Å. In contrast to µex ) 0 kJ/mol around r ) 0 in pore L, µex is about -6.5 kJ/mol in pore M. These results suggest that the potential field near the pore surface is independent of the curvature of pores, although the local potential field generated by each part of the pore surface overlaps near the center of pore M. Although the µex profile at N ) 48 of pore M is almost the same as that of pore L at 5 Å

Phase Behavior of Water Confined in Silica Nanopores

J. Phys. Chem. C, Vol. 111, No. 22, 2007 7943

Figure 6. Typical adsorption structure of a water molecule around silanol groups without hydrogen-bonding interaction. Blue atom shows the oxygen atom of the water molecule, red atoms show oxygen atoms and brown atoms show silicon atoms of silanol groups, and white atoms show hydrogen atoms.

Figure 5. Profiles of excess chemical potential, µex [kJ/mol], along the r-axis in pore S (top), pore M (middle), and pore L (bottom). Curves were obtained by fitting the calculated µex shown as cross marks.

< r < 14 Å, when N is increased further µex becomes constant along the r-axis at N ) 288. In pore M, monolayer adsorption apparently does not occur. If the µex profile at N ) 48 (M-I) is dislocated 2 Å toward the center of the pore as it is in pore L, µex around r ) 0 is about -15 kJ/mol, that is, if a water monolayer is formed, then µex around r ) 0 is close to that of bulk liquid water (-22 kJ/mol). This result suggests that rather than adsorb on the surface and form a water monolayer, water molecules condense inside the pore. In pore S, the µex profile at dehydration (N ) 0) shows that the minimum µex is -49 kJ/mol around r ) 3 Å and that µex at r ) 0 is -22 kJ/mol. The minimum µex in pore S is larger than that in pores M and L. Besides, µex at r ) 0 is smaller than that in pore M or L because the local potential field generated by each part of the pore surface overlaps near the center of pore S. At the highest hydration (N ) 54), µex at r ) 0 is about -10 kJ/mol, suggesting that water molecules are so closely packed that there is smaller space available for adsorption at the center of the pore. 3.2. Structural Properties. On the basis of the simulation results, if the distance between OW and OS*, namely, rOW-OS*, is less than 3.2 Å and the angle between OS*, HW, and OW, namely, ∠OS*-HW-OW, is larger than 120°, the interaction between HW and OS* can be defined as a hydrogen bond and the absorption of water molecules as a PD, as shown in Figure 2 (top). In contrast, if rOW-OS* < 3.2 Å and ∠OW-HS*-OS* > 120°, then the interaction between OW and HS* can be defined as a hydrogen bond and the absorption of water molecules as a PA, as shown in Figure 2 (bottom). In these definitions, both rHW-OS* and rOW-HS* are less than 2.58 Å. At the highest hydration in this study (N ) 54, 288, and 672 for pores S, M, and L, respectively), the ratio of the silanol groups on which water molecules are adsorbed as PD and PA simultaneously is 0.2, 0.4 and 0.4% in pores S, M, and L, respectively. The ratio of silanol groups on which water molecules are adsorbed as PD is 16.3, 22.3, and 21.7% in pores

Figure 7. Number density of oxygen atoms of water molecules (OW) around oxygen atoms of silanol groups (OS*) reduced by the averaged number density of water inside the pore (N/Veff [cm-3]), n*OW-OS*, as a function of the distance between OS* and OW, rOS*-OW [Å]. N is the number of water molecules and Veff is effective pore volume defined as the volume of cylindrical pores. The diameter of cylindrical pores S, M, and L was 0.89, 1.81, and 2.77 nm, respectively.

S, M, and L, respectively, and that as PA is 5.1, 13.1, and 12.0%, respectively. Other silanol groups show no hydrogen-bonding interaction with water molecules in this definition. Figure 6 shows a typical adsorption structure of water molecules around silanol groups without hydrogen-bonding interaction. In this structure, a water molecule interacts with two adjacent silanol groups via Coulomb force between HS*-OW (rOS*-OW ∼ 3.6 Å and θ ) ∠OS*-HS*-OW ∼ 100°). Figure 7 shows n*OW-OS* as a function of rOW-OS* in pores S, M, and L, where n*OW-OS* is defined as the number density of OW around OS* normalized by the averaged number density of water inside the

7944 J. Phys. Chem. C, Vol. 111, No. 22, 2007

Figure 8. Profiles of calculated water density, F, along the r-axis in pore S (top), pore M (middle), and pore L (bottom). Dashed lines represent local minima of the density profile to show the boundary of the nth adsorption layer at full hydration (at N ) 288 in pore M and N ) 672 in pore L).

pore (N/Veff). In pores M and L, at low hydration (N ) 48), the peaks are located at rOW-OS* ≈ 2.9 Å, which corresponds to the length of a hydrogen bond between OW and OS*. This result suggests that water molecules interact with silanol groups mainly via hydrogen bonds at low hydration. At high hydration (N ) 288 in pore M, and N ) 288 and 672 in pore L), the maximum peaks are located at rOW-OS* ≈ 3.6 Å, which corresponds to the non-hydrogen-bonded structure as shown in Figure 6. This adsorption structure was not observed at low hydration in pores M and L. In pore S, the ratio of silanol groups that interact with water molecules via hydrogen-bonding interactions is less than that of pores M and L. The non-hydrogen-bonded structure was observed at either N ) 12 or 54. Smirnov et al.8 reported the X-ray radial distribution functions characterized by the nonhydrogen-bonded OW-OW interactions at ∼3.3 Å and OWSi interactions at ∼3.8 Å in the hydrated cylindrical pores of MCM-41, suggesting that the OW-OS* interactions are not negligible in X-ray radial distribution functions. Figure 8 shows the profiles of water density, F, along the r-axis at various hydration levels N for pores S, M, and L. In pore L, water molecules adsorb on the surface at N ) 48 (L-I). With increasing N, the first layering phase transition occurs and the first adsorption layer is formed at N ) 288 (L-II). Because the number of silanol groups in pore L is 288, the number of water molecules to form the monolayer is almost the same as the number of silanol groups. At N ) 336, the water monolayer still remains and neither multilayer adsorption nor condensation occurs. With further increase in N, liquid-vapor transition occurs and the pore is filled with water at N ) 624 (L-III), although F near the center of the pore at N ) 624 is smaller than that of bulk liquid water. On the basis of these F profiles and the µ-Fj curves shown in Figure 3, the distribution of water at N ) 336 and 624 is apparently not in the stable state but in

Shirono and Daiguji

Figure 9. Profiles of calculated local diffusion coefficient normalized by the diffusion coefficient of bulk liquid water (2.49 × 10-9 m2/s), D*, along the r-axis in pore S (top), pore M (middle), and pore L (bottom). Dashed lines represent local minima of the density profile to show the boundary of the nth adsorption layer at full hydration (at N ) 288 in pore M and N ) 672 in pore L).

the metastable state. At the highest hydration (N ) 672), four adsorption layers could be observed. The density oscillation in the third and fourth adsorption layers is much smaller than that in the first and second adsorption layers. In pore M, water molecules adsorb on the surface at N ) 24 and 48 (M-I). Before water molecules completely cover the surface, capillary condensation starts and large voids disappear at N ) 240. At full hydration (N ) 288), three adsorption layers are formed. The fluctuation in F in the first and second adsorption layers is much larger than that in the third adsorption layer, suggesting that water molecules in the first and second adsorption layers have different structure from that of bulk liquid water. In pore S, because the radius of pore S is almost the same as the thickness of the first adsorption layer, all water molecules are located inside the first adsorption layer. Water molecules form the adsorption layer on the pore surface, thus F at r ) 0 is almost zero even at high hydration (N ) 48 and 54). 3.4. Diffusion Coefficient. Figure 9 shows the profiles of calculated D* along the r-axis for pores S, M, and L, where D* is defined as the calculated local diffusion coefficient normalized by the diffusion coefficient of bulk liquid water (Db ) 2.49 × 10-9 m2/s). For all three pores, D* increases from the pore surface to the pore center. In pore S, D* < 1 except for the center of the pore at N ) 12 and D* < 1 over entire range of r at N ) 54, suggesting that water molecules are constrained on the surface due to strong interaction between water molecules and silanol groups. In M-I at N ) 48, D* < 1 in the first adsorption layer but increases rapidly with decreasing r in the second adsorption layer, and the water molecules around the center of the pore behave like water vapor. In M-II at N ) 288, D* increases with decreasing r (in the first and second adsorption layers), and D*at the center of the pore is about 0.8, suggesting

Phase Behavior of Water Confined in Silica Nanopores

Figure 10. Calculated diffusion coefficients along the z-axis normalized by the diffusion coefficient of the SPC/E water model in the bulk liquid (2.49 × 10-9 m2/s), DZ*, as a function of average water density, Fj [g cm-3], in pores S, M, and L. Fj is given by N mw/Veff, where N is the number of water molecules, mw [g] is the mass of water molecules and Veff [cm3] is effective pore volume defined as the volume of cylindrical pores in which water molecules were inserted. The diameter of cylindrical pores S, M, and L was 0.89, 1.81, and 2.77 nm, respectively. The number at each symbol expresses N.

that water molecules near the center of the pore are still constrained on the surface due to strong interaction between water molecules and silanol. The D* profiles in L-I at N ) 48 and L-II at N ) 288 are similar to that in M-I at N ) 48, whereas the D* profile in L-III at N ) 672 is similar to that in M-II at N ) 288. In L-III, D* near the center of the pore is about 1, suggesting that water molecules near the center of the pore are not constrained on the surface and behave like liquid water. In conclusion, D* increases with decreasing r in the first and second adsorption layers and then reaches 1 in pore L, whereas in pore M, D* is smaller than 0.8 even at the center of the pore. Figure 10 shows Dz* as a function of Fj in pores S, M, and L, where Dz* is defined as the calculated diffusion coefficient along the z-axis normalized by the diffusion coefficient of the SPC/E water model in the bulk liquid. Results show that Dz* increases with increasing r0 at a fixed Fj, and that Dz* increases with increasing N at a fixed r0 except for high hydration (N ) 48 and 54 in pore S, N ) 240 and 288 in pore M, and N ) 624 and 672 in pore L). At high hydration, F near the center of the pore increases with increasing N, thus the diffusivity near the center of the pore deceases and approaches that in the bulk liquid water. Takahara et al.11 conducted neutron diffraction measurements to clarify the diffusion coefficient of water molecules in hydrated cylindrical pores of 2.14 and 2.84 nm diameter MCM41 at 298 K. Their measured diffusion coefficient was 1.07 × 10-9 m2/s for the 2.14 nm diameter pore (i.e., 0.43Db) and 1.45 × 10-9 m2/s for the 2.84 nm diameter pore (i.e., 0.58Db). The calculated Dz* values were 0.80 × 10-9 m2/s for pore M (i.e., 0.32Db) and 1.42 × 10-9 m2/s for pore L (i.e., 0.57Db) at the highest hydration (N ) 288 and N ) 672 for pore M and pore L, respectively) agree well with these measured values. 4. Conclusions MD and GCMC simulations of 1.04, 1.96, and 2.88 nm diameter hydrated silica pores (pores S, M, and L, respectively) were conducted at 300 K. The following conclusions could be drawn from this study.

J. Phys. Chem. C, Vol. 111, No. 22, 2007 7945 1. A molecular model of silica pores was proposed that was based on R-quartz and the pair potential functions between water molecules and silica pores to reproduce two adsorption structures, that is, a structure in which water molecules are adsorbed on the silanol group as a PD and a structure in which water molecules are adsorbed as a PA. 2. Three types of phases appear in nanopores. In the first phase, which appears in pores S, M, and L, gas-phase water adsorption at the pore surface is below monolayer coverage. In the second phase, which appears in pores S and L, there is a condensed water monolayer at the pore surface. In the third phase, which appears in pores M and L, the pore is completely filled with water. A water monolayer observed in the second phase is due to capillary condensation in pore S, while due to layering transition in pore L. 3. At low hydration, water molecules interact with silanol groups mainly via hydrogen bonds. With increasing water density, the number ratio of water molecules that adsorb onto silanol groups via non-hydrogen-bonding interactions increases. At full hydration, only one adsorption layer is present in pore S, whereas three and four adsorption layers are present in pores M and L, respectively. The large density oscillation could be observed only in the first and second adsorption layers from the surface, suggesting that water molecules in the first and second adsorption layers have different structure from that of bulk liquid water. 4. The radial dependence of the translational mobility of water was clarified at various hydration levels. Near full hydration, the translational mobility of water in the first adsorption layer is much smaller than that of bulk liquid water in all three pores, while those in the pore center are about 30, 80, and 100% of its bulk liquid value in pores S, M, and L, respectively. In the axial direction, the translational mobility of water molecules increases with increasing pore diameter. In all three pores, the translational mobility increases with increasing water density but decreases near full hydration. References and Notes (1) Beck, J. S.; Vartuli, J. C.; Roth, W. J.; Leonowicz, M. E.; Kresge, C. T.; Schmitt, K. D.; Chu, C. T-W.; Olson, D. H.; Sheppard, E. W.; McCullen, S. B.; Higgins, J. B.; Schlenkert, J. L. J. Am. Chem. Soc. 1992, 114, 10834. (2) Lin, H.-P.; Mou, C.-Y. Acc. Chem. Res. 2002, 35, 927. (3) Brovchenko, I.; Geiger, A.; Oleinikova, A. J. Chem. Phys. 2004, 120, 1958. (4) Puibasset, J.; Pellenq, R. J. M. J. Chem. Phys. 2005, 122, 094704. (5) Inagaki, S.; Fukushima, Y. Microporous and Mesoporous Mater. 1998, 21, 667. (6) Branton, P. J.; Hall, P. G.; Sing, K. S. W.; Reichert, H.; Schu¨th, F. S.; Unger, K. K. J. Chem. Soc., Faraday Trans. 1994, 90, 2965. (7) Rovere, M.; Gallo, P. Eur. Phys. J. E 2003, 12, 77. (8) Smirnov, P.; Yamaguchi, T.; Kittaka, S.; Takahara, S.; Kuroda, Y. J. Phys. Chem. B 2000, 104, 5498. (9) Crupi, V.; Dianoux, A. J.; Majolino, D.; Migliardo, P.; Venuti, V. Phys. Chem. Chem. Phys. 2002, 4, 2768. (10) Takahara, S.; Nakano, M.; Kittaka, S.; Kuroda, Y.; Mori, T.; Hamano, H.; Yamaguchi, T. J. Phys. Chem. B 1999, 103, 5814. (11) Takahara, S.; Sumiyama, N.; Kittaka, S.; Yamaguchi, T.; BellissentFunel, M-C. J. Phys. Chem. B 2005, 109, 11231. (12) Mansour, F.; Dimeo, R. M.; Peemoeller, H. Phys. ReV. E 2002, 66, 041307. (13) Takamuku, T.; Yamagami, M.; Wakita, H.; Masuda, Y.; Yamaguchi, T. J. Phys. Chem. B 1997, 101, 5730. (14) Bellissent-Funnel, M.-C. Eur. Phys. J. E 2003, 12, 83. (15) Sonwane, C. G.; Jones, C. W.; Ludovice, P. J. J. Phys. Chem. B 2005, 109, 23395. (16) Coasne, B.; Hung, F. R.; Pellenq, R. J. M.; Siperstein, F. R.; Gubbins, K. E. Langmuir 2006, 22, 194. (17) Wright, A. F.; Lehmann, M. S. J. Solid State Chem. 1981, 36, 371. (18) Tsuneyuki, S.; Tsukada, M.; Aoki, H.; Matsui, Y. Phys. ReV. Lett. 1988, 61, 869. (19) Zhuravlev, L. T. Colloids Surf., A 2000, 173, 1.

7946 J. Phys. Chem. C, Vol. 111, No. 22, 2007 (20) Coasne, B.; Pellenq, R. J. M. J. Chem. Phys. 2004, 120, 2913. (21) Iarlori, S.; Ceresoli, D.; Bernasconi, M.; Donadio, D.; Parrinello, M. J. Phys. Chem. B 2001, 105, 8007. (22) Brendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (23) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.;

Shirono and Daiguji Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, revision A.11; Gaussian, Inc.: Pittsburgh, PA, 1998. (24) Nose´, S. Mol. Phys. 1984, 52, 255. (25) Nose´, S. J. Chem. Phys. 1984, 81, 511. (26) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (27) Widom, B. J. Chem. Phys. 1963, 39, 2808. (28) Jedlovszky, P.; Mezei, M. J. Am. Chem. Soc. 2000, 122, 5125. (29) Whitley, H. D.; Smith, D. E. J. Chem. Phys. 2004, 120, 5387. (30) Brovchenko, I.; Geiger, A.; Oleinikova, A.; Paschek, D. Eur. Phys. J. E 2003, 12, 69.