Molecular Simulations of Water Adsorbed on ... - ACS Publications

Jan 10, 2013 - Division of Environmental Studies, Graduate School of Frontier ... Manuel I. Velasco , M. Belén Franzoni , Esteban A. Franceschini ...
0 downloads 0 Views 3MB Size
Subscriber access provided by FORDHAM UNIVERSITY

Article

Molecular Simulations of Water Adsorbed on Mesoporous Silica Thin Films Kyohei Yamashita, and Hirofumi Daiguji J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/jp312804c • Publication Date (Web): 10 Jan 2013 Downloaded from http://pubs.acs.org on January 13, 2013

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular Simulations of Water Adsorbed on Mesoporous Silica Thin Films Kyohei Yamashita†‡ and Hirofumi Daiguji†‡*

†Division of Environmental Studies, Graduate School of Frontier Sciences, The University of Tokyo 5-1-5 Kashiwanoha, Kashiwa 277-8563, Japan ‡CREST, Japan Science and Technology Agency (JST) 7 Gobancho, Chiyoda-ku, Tokyo 102-0076, Japan

*To whom correspondence should be addressed. E-mail: [email protected], Phone: +81 4 7136 4658, Fax: +81 4 7136 4659.

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 40

ABSTRACT

Grand canonical Monte Carlo and canonical ensemble molecular dynamics simulations were used to investigate water adsorption properties in mesoporous silica thin films. The effect of pore radius on the adsorption properties was assessed using two models of mesoporous silica thin films having different pore radius and film thickness (1.38 and 5.66 nm in Model 1, respectively, and 1.81 and 7.30 nm in Model 2, respectively). In the simulations, a water adsorption layer or water menisci were formed in a mesopore accompanying the growth or shrinkage of stable adsorption layers on the upper and lower surfaces. The thickness of a water adsorption layer immediately before the onset of capillary condensation and the curvature radius of a water meniscus prior to capillary evaporation were identified. In the initial stage of layer adsorption, the surface state and radius of the pore affected the coordination structure of water around silanol groups. In the pore-filling state, the hydrogen bond network of water in the pore was similar but not the same as that in the bulk liquid; the difference was because the pore wall affected the position of water molecules. This difference may confer unique dynamic properties upon the water in mesoporous silica.

Keywords: Grand canonical Monte Carlo simulations, Molecular dynamics simulations, Adsorption isotherms, Adsorption layers, Menisci, Structural and dynamic properties

ACS Paragon Plus Environment

2

Page 3 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1. INTRODUCTION

It is known well that physical and thermodynamic properties of a liquid confined in nanopores are different to those of the bulk phase. These differences are due to the fact that the dimensions of the system are comparable to the range of intermolecular interactions, and thus, the liquid–liquid and liquid–pore wall interactions have to be considered explicitly in order to understand these properties. Water is one of the most essential molecules in nature, and it exhibits several unique properties such as large surface tension and large heat of vaporization, due to its hydrogen bond network. One of the great challenges in the study of water is to understand its various properties under nanoscale confinement and to exploit its confined properties to develop novel engineering systems. In this study, we focused on mesoporous silica as an adsorbent of water. Mesoporous silica can adsorb and desorb water vapor over a specific range of relative humidity, thus it may have potential applications in energy-efficient desiccant and humidity control systems. Mesoporous silica is an important class of porous silica materials having pore sizes of 2 to 50 nm.1,2 These materials have well-uniformed pore size distributions and highly ordered pore directions, which can be attributed to the template mechanism involving micelles formed by self-assembly of surfactant or block copolymer into the silica source solvent. In addition, pore size and porous structure can be controlled by modifying synthesis conditions. Consequently, mesoporous silica materials have received much attention for use as catalysts, filters, and adsorbents. Adsorption/desorption of water in mesoporous silica exhibits phase-change-like behavior in adsorption/desorption processes with a pronounced hysteresis loop.3-5 These phenomena consist of the following processes: (1) formation/deformation of a water adsorption layer on the pore surface, (2) abrupt increase/decrease in the amount of water by capillary condensation/evaporation, and (3) adsorption /desorption in the water filled pore. Water adsorption/desorption of mesoporous silica has been investigated experimentally. Inagaki et al.3 thrice recorded the water adsorption isotherm at 298 K of mesoporous silica FSM-16, having a pore radius from 1.4 to 2.6 nm. The results revealed that the ACS Paragon Plus Environment

3

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 40

hysteresis width between capillary condensation and evaporation in the adsorption/desorption isotherm decreased with decreasing pore radius and that hysteresis disappeared at a pore radius of 1.4 nm. Moreover, the adsorption isotherms for the second and subsequent measurements were shifted to a region of lower relative pressure of water vapor, p/p0, where p is the actual pressure of water vapor and p0 is the vapor pressure at saturation, because the surface of the silica wall was rehydroxylated following the first adsorption-desorption measurement. Russo et al.4 measured water adsorption isotherms at 298 K on Al grafted MCM-41 having different Si/Al molar ratios. As the Si/Al molar ratios decreased, capillary condensation commenced at a lower relative pressure, which was associated with an increase in the strength of the adsorbate–adsorbent interaction and a decrease in the pore size. Endo et al.5 measured water adsorption isotherms on Zr-doped mesoporous silica (Zr-MPS) between 263 and 298 K. Zr-MPS features cylindrical mesopores that are uniformly arranged in a 2D-hexagonal structure. The width of hysteresis loop slightly increased with decreasing temperature. The results of these studies demonstrated that pore radius, surface properties, and temperature have an affect the adsorption properties. The adsorption/desorption behavior of water in mesoporous silica is characterized by capillary condensation and evaporation. p/p0 at which capillary condensation and evaporation occur can be described as a function of curvature radius, ρ, surface tension, γ, molar volume, Vm, of adsorbates, and temperature, T. This relationship is known as the Kelvin equation: RT ln( p p 0 ) = − 2γVm ρ , where R is the gas constant. The Kelvin equation is usually used to explain adsorption phenomena of simple molecules in large pores since it is assumed that the pressure reduction within the pore can be described by the Young-Laplace effect and the chemical potential of adsorbates in the pore is identical to that in the bulk liquid.6 However, in water-mesoporous silica systems, the discreteness of molecules cannot be neglected, and thus, the Kelvin equation cannot be used without modification. In order to predict water adsorption in mesoporous silica more accurately and to design new adsorbents of water, it is important to understand the adsorption phenomena at a molecular level. Molecular simulations are an effective tool for elucidating the discreteness of molecules. ACS Paragon Plus Environment

4

Page 5 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular simulations of water adsorption/desorption of mesoporous silica have been performed by several researchers. For example, Puibasset et al.7 performed the grand canonical Monte Carlo (GCMC) simulation of water adsorption/desorption in a Vycor-like mesoporous silica with an average pore radius of about 1.8 nm at 300 and 650 K. The adsorption/desorption isotherm exhibited hysteresis at 300 K, while no hysteresis was observed at 650 K. Bonnaud et al.8 performed GCMC simulation of water adsorption/desorption in an amorphous silica slit pore having widths of 1.0 and 2.0 nm at 300 K. Results revealed that the adsorption/desorption isotherm exhibited hysteresis for the slit pore with a diameter of 2.0 nm, while hysteresis was not seen for the slit pore with a diameter of 1.0 nm. Furukawa et al.9 carried out orientational-bias GCMC simulations for calculating the adsorption isotherms of pure acetone and water at 298 K on non-silylated (hydrophilic) and silylated (hydrophobic) mesoporous silica with a pore radius of about 1.5 nm. A good agreement was obtained between simulations and experiments. Shirono et al.10 conducted canonical ensemble molecular dynamics (NVT-MD) and GCMC simulations of water molecules in silica pores having a radius of 0.52, 0.98, and 1.44 nm at 300 K. The adsorption/desorption isotherms were calculated using both the NVT-MD simulation plus cavity insertion Widom (CIW) method and also the GCMC simulation. Both results showed that the phase behavior of water confined in the nanopores was dependent upon the pore radius. de la Llave et al.11 investigated the phase behavior of adsorbed water in cylindrical hydrophilic pores with diameter of 1.5 and 3 nm. They performed NVT-MD simulations of confined water as a function of the amount of adsorbed water. They showed that the hysteresis of adsorption-desorption isotherms develops above ~1.5 nm in diameter and capillary condensation appears at 27 and 34% filling for 1.5 and 3 nm diameter pores, respectively. These results were in quantitative agreement with experimental ones. In order to realize a more accurate prediction and control of adsorption/desorption isotherms and adsorption and transport processes, the objectives of this study were to elucidate the thickness of a water adsorption layer immediately before the onset of capillary condensation and the curvature radius of a water meniscus prior to capillary evaporation. We also sought to clarify structural and dynamic properties of confined water in different adsorption states. In the past, the structural and dynamic ACS Paragon Plus Environment

5

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 40

properties of confined water have been investigated by many researchers.12-20 In the majority of previous studies, the calculation systems employed were either nanopores or nanoslits with periodic boundary conditions. When the adsorption/desorption isotherm exhibits hysteresis, adsorbed water immediately preceding capillary condensation or evaporation is in the metastable state, and thus stable water adsorption layers or bubbles cannot be formed. As a result, it is difficult to identify the thickness of a water adsorption layer and the curvature radius of a water meniscus. However, if the calculation system uses mesoporous silica thin films, i.e., a finite cylindrical-like pore and the upper and lower surfaces are included in the system, a water adsorption layer or water menisci can be formed in a mesopore accompanying the growth or shrinkage of stable adsorption layers on the upper and lower surfaces. As a result, it should be possible to identify the thickness of a water adsorption layer and the curvature radius of a water meniscus.21 In this study, we performed GCMC and NVT-MD simulations of water adsorbed on mesoporous silica thin films at 300 K.

2. METHODS

2.1 Modeling of Mesoporous Silica Thin Films A mesoporous silica thin film was modeled as a thin film of α-quartz crystal with a finite cylindrical pore normal to the film surfaces. The crystal structure was sliced by two parallel planes (top and bottom planes) and a finite cylindrical pore was aligned normal to the two parallel planes. The pore surface and the top and bottom planes were fully hydroxylated with different kinds of surface functional groups: vicinal silanol groups (one hydroxyl group is bonded to a Si atom) and geminal silanol groups (two hydroxyl groups are bonded to a Si atom). In this study, two mesoporous silica thin film models having a different pore radius and film thickness were employed. The pore radius r0 and film thickness L0 were 1.44 and 5.40 nm for Model 1 and 1.87 and 7.00 nm for Model 2, respectively. The mesoporous silica thin film model was located in a rectangular parallelepiped simulation cell. The dimensions of a simulation cell and a thin film model were identical in the film tangential directions, whereas, the ACS Paragon Plus Environment

6

Page 7 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

dimension of a simulation cell was about twice larger than that of a thin film model in the film normal direction. The upper and lower surfaces of the thin film were in contact with the vapor phase. The modeling entailed the following 8-step procedure. In the first step, l, m, and n trigonal unit cells of α-quartz, whose structure has been described in detail by Wright and Lehmann,22 was arrayed in the [100], [010], and [001] directions, respectively. The values of l, m, and n were 9, 10, and 20 for Model 1 and 11, 12, and 26 for Model 2, respectively. In the second step, this α-quartz crystal was transformed from a rhombohedron to a rectangular parallelepiped.10 The dimensions of this rectangular parallelepiped were identical to those of a simulation cell of which volume was 4.42 × 4.26 × 10.81 nm3 for Model 1 and 5.40 × 5.11 × 14.05 nm3 for Model 2, respectively. In the third step, an arbitrary pore axis (z-axis) was selected parallel to the [001] direction. Si atoms within a radius r0 from the z-axis were removed. In the fourth step, O atoms that were covalently bonded to two adjacent removed Si atoms were also removed. On the other hand, O atoms bonded to only one Si atom were saturated with H atoms to prevent dangling bonds on the silica surface. In the fifth step, if all silanol groups on the pore surface produced in the fourth step were not vicinal silanol groups, then the third step was revisited and a different z-axis was selected. In the sixth step, arbitrary parallel planes with an interplanar distance of L0 were selected normal to z-axis. Si and O atoms were removed if they were located outside the thin film with a thickness of L0. As described in the fourth step, H atoms were added to dangling O atoms. In the seventh step, Si atoms bonded to more than three OH groups, were removed. By removing Si atoms, new dangling O atoms were produced. Such dangling O atoms were saturated with H atoms. In the eighth step, to minimize stress in the atom model thin film, NVT-MD simulations were conducted with all atoms mobile at 300 K. Figure 1 shows the final structure of the model mesoporous silica thin film. After NVT-MD simulations, the pore shape in top view became hexagonal, thus the pore radius rp was defined as the average distance between the z-axis and O atoms of vicinal silanol groups. The obtained values of rp were 1.38 and 1.81 nm, respectively, for Models 1 and 2. Furthermore, the film thickness was enlarged by NVT-MD simulations. The obtained film thicknesses Lp were 5.66 and 7.30 nm, respectively, for Models 1 and 2. The OH surface number density of the mesoporous silica thin film ACS Paragon Plus Environment

7

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 40

models was approximately 7.2 OH nm-2 for the inner surface of a silica pore and 9.6 OH nm-2 for the upper and lower surfaces of a thin film. The value for the inner surface of a silica pore was in good agreement with that in other experimental and computational studies.8, 21, 23

2.2 Potential Functions and Parameters In this study, the intermolecular interactions between silica/silanol groups and silica/silanol groups– water involved both Coulomb and Born-Mayer-Huggins (BMH) interactions. Specifically,

U ij (rij ) =

Zi Z je2 rij

 ai + a j − rij + f 0 (bi + b j )exp  b +b i j 

 ci c j − ,  r6 ij 

(1)

where Uij(rij) is the intermolecular potential energy between atoms i and j as a function of interatomic distance of rij. Z is the formal ionic charge and e is the elementary electric charge. The terms ai (Å), bi (Å) and ci (eV1/2 Å3) are the effective radius, softness, and van der Waals parameters of the ith atom, respectively. The term f0 is a conversion factor in which f0 = 1 kcal Å-1 mol-1 = 4.3384 × 10-2 eV Å-1. For O and Si atoms of silica/silanol groups (Os/Os*, Sis/Sis*), the parameters of the Tsuneyuki potential model24 were used. For H atoms (Hs*) of silanol groups, only coulomb interaction was considered to all other atoms with formal ionic charge of Z = 0.6 to ensure the electroneutrality of the entire system. For water, the SPC/E rigid water model25 was employed. The SPC/E rigid water model consists of three simple point charges, i.e., O atom (Ow) and two H atoms (Hw). In this model, intermolecular interaction is described by a Lennard-Jones (LJ) potential function between Ow–Ow where the OH bond length is 1.0 Å and the angle between two OH bonds is 109.47°. The potential parameters of BMH functions between Ow–Os/Os* were calculated by fitting the LJ potential function of the SPC/E model to the BMH potential function.10 In this study, BMH interactions between Ow– Sis/Sis*, Ow–Hs*, Hw–Sis/Sis*, Hw–Os/Os*, and Hw–Hs* were not considered. The potential parameters used in this study are summarized in Table 1. ACS Paragon Plus Environment

8

Page 9 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

2.3 GCMC and NVT-MD Simulation Details In this study, GCMC simulations were conducted to obtain the adsorption isotherms of water in mesoporous silica at 300 K. The adsorption process started from an empty mesoporous silica thin film model and the chemical potential increased step-by-step. The equilibrium configuration of i th state was used as the initial configuration of (i + 1) th state. After reaching the saturation state, the desorption process started with decreasing the chemical potential. The GCMC consisted of the following four types of trial moves: translation, rotation, insertion, and removal. For each state, more than 4.0 × 105 trials were attempted for these four types of moves per water molecule for equilibration. Subsequently, more than 6.0 × 104 trials were attempt per water molecule to evaluate the number of adsorbed water molecules. The translation and rotation moves were modulated to obtain the acceptance ratio of 0.5, while the insertion and removal moves had much lower acceptance ratio. The acceptance ratio of insertion and removal moves after equilibration was about 0.005 – 0.01 and the ratio depended on the adsorption states. Entire area in a simulation cell except for the area where silica molecules were occupied was selected as target area for the insertion of water molecules. In GCMC simulations, to minimize calculation time, Si and O atoms of silica were fixed at the initial positions obtained by the procedure described in Section 2.1. The geometry of two kinds of silanol groups were fixed with bond lengths of Sis*–Os* and Os*–Hs* of 1.61 and 1.0 Å, respectively. The Sis*–Os*–Hs* angle was 109.47°. Only rotational moves around Sis* atoms were allowed for silanol groups. The chemical potential, µ, was calculated as the summation of the ideal chemical potential, µid, and the excess chemical potential, µex. The ideal chemical potential involves rotational contribution of SPC/E water molecule and given by the following equation.

 π 1 2 µ = kT ln nΛ3 − ln   σ  id

( )

12

 8π 2 I A kT    ∏ h2 A   x, y ,z

   ,  

ACS Paragon Plus Environment

(2)

9

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 40

where k is Boltzmann’s constant, n is number density of water, Λ is thermal de Broglie wavelength, σ is the symmetry number of a molecule (σ = 2 for SPC/E model), Ix, Iy and Iz are the principal moments of inertia of SPC/E model, h is Planck’s constant. NVT-MD simulations were performed to investigate the structural and dynamic properties of water adsorbed on the mesoporous silica thin films. Initially, N water molecules were placed in the system for NVT-MD simulations. The velocity Verlet method26 was implemented to integrate the equation of motion for all atoms in the calculation system with a time step of 1.0 fs. The temperature of system was maintained at 300 K using the Nose thermostat.27 The geometry of two kinds of silanol groups and SPC/E water molecules were constrained by RATTLE algorithm.26 The bond lengths of Sis*–Os* and Os*–Hs* and the Sis*–Os*–Hs* angle were constrained to be 1.61 and 1.0 Å and 109.47°, respectively. Thus, the numbers of degree-of-freedom of vicinal and geminal silanol groups were 6 and 9, respectively. The long-range coulomb interactions were calculated using Ewald method26 and periodic boundary conditions were applied in all directions. The cut off radius of intermolecular interactions and real space term of Coulomb interactions was 12.5 Å.

3. RESULTS AND DISCUSSION

3.1 Adsorption Properties Figure 2 shows the snapshots of NVT-MD simulations of water adsorbed on the mesoporous silica thin film of Model 1. Layer adsorption and pore filling were observed at N = 500 and 1200, respectively. Two different states of layer adsorption and pore filling could be obtained depending on the number of water molecules, N in both Models 1 and 2. Adsorbed water molecules were categorized into two groups depending on their position: (1) on the upper or lower surface ( z > z 0 ) and (2) in the pore ( z ≤ z 0 ), where z 0 = 2.525 and 3.35 nm in Models 1 and 2, respectively. The number of water molecules adsorbed on the upper or lower surface, N1, and in the pore, N2, were counted at each adsorption state. Figure 3a shows the N1–N2 curves calculated from ACS Paragon Plus Environment

10

Page 11 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

GCMC and NVT-MD simulations. These results are in good agreement with each other. Assuming Henry’s law of adsorption, i.e., N1 is proportional to the pressure of the calculation system, the N1–N2 curves will become the adsorption/desorption isotherms of water in Models 1 and 2. Because the number of water molecules in the gas phase was almost zero under all calculation conditions, Henry’s law of adsorption is a good approximation. Even if Henry’s law of adsorption is not applicable, N1 should still be a function of pressure, thus the N1–N2 curves and the adsorption/desorption isotherms should have similar shapes. Figure 3b shows the chemical potential of water vs. the number of water molecules (µ–N) curves calculated from GCMC simulations and Figures 3c and 3d show the decomposed µ–N curves: the chemical potential of water vs. the area number density of water molecules on the upper or lower surface (µ–n1) curves (n1 was defined as N1 per unit upper or lower surface area) and the chemical potential of water vs. the number density of water molecules in the pore (µ–n2) curves (n2 was defined as N2 per unit pore volume), respectively. For the µ–n1 curves, n1 remained almost constant at low µ regions and increased rapidly above a threshold value of µ. The µ–n1 curves of Models 1 and 2 were approximately same. These curves did not exhibit hysteresis, suggesting that the chemical potential of water adsorbed on the upper or lower surface depended on the area number density of water molecules n1 and did not depend on the adsorption/desorption process. On the other hand, the µ–n2 curves are a function of the pore models and exhibited hysteresis. These curves expressed the adsorption/desorption isotherms. The capillary condensation and evaporation began at higher µ and wider hysteresis could be observed in the larger pore radius of Model 2. In order to investigate the validity of GCMC calculation, water adsorption-desorption isotherms were also predicted by the Derjaguin-Broekhoff-de Boer (DBdB) theory and compared with the results of GCMC calculation. Figure 3d also shows the theoretical prediction of µ–n2 curves based on DBdB theory (see supporting information S1).4,28 The solid and dashed lines correspond to the µ–n2 curves of pore radius of 1.38 and 1.81 nm, respectively. The µ–n2 curves obtained by GCMC calculation and a theoretical model based on DBdB theory agree with each other. Furthermore, to evaluate the internal energy of the system, the

ACS Paragon Plus Environment

11

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 40

differential heat of adsorption was calculated (see supporting information S2). The calculated differential heat of adsorption was close to the heat of vaporization of water. Figure 4 shows the contour plots of water density as a function of N in Models 1 and 2. The density contour plots were calculated by assuming that the spatial and time average of water density is symmetrical around the z-axis (see supporting information S3). As for the adsorption process, stable water adsorption layers were formed discreetly and extra water that was not present in the water adsorption layer flowed out to the upper or lower surfaces. Figures 5a and 5b show the profiles of the vapor-liquid interface in the layer adsorption process for Models 1 and 2, respectively. The profiles did not change monotonically with respect to N, and the first and the second water adsorption layers appeared discreetly. These results suggest the stable formation of the first and second water adsorption layers; however, the onset of capillary condensation occurs prior to the formation of a third layer. Previous simulation studies of water confined in silica pore12-14 and cylindrical pore with hydrophilic smooth surface15,16 reports the confined water exhibit specific two water layers were formed on the pore surface. These reports agree with stable adsorbed water layers were formed up to second layer on the surface of silica pore. Figure 5c shows the plot of average thickness of adsorption layer versus the number of water molecules in the whole system and in the pore (t–N and t–N2) curves. The average thickness of adsorption layer, t, increased with an increase in the number of water molecules, and t can be expressed as a function of N2 by assuming the cylindrical adsorption layer onto the pore surface with uniform density:

t (N 2 ) = rp − rp2 − vN 2 2 z 0π ,

(3)

where v is the average volume of one water molecule. The calculated t–N2 curves were fitted by Eq. 3 within the range of 397 ≤ N2 ≤ 531 (500 ≤ N ≤ 800) in Model 1 and 834 ≤ N2 ≤ 1345 (1000 ≤ N ≤ 1900) in Model 2 using v as a parameter. Figure 5c also shows the fitted curves. The values for v were 41.26 12 ACS Paragon Plus Environment

Page 13 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

and 37.57 Å3 for Models 1 and 2, respectively, suggesting that the average water density in the adsorption layer of Model 2 was slightly higher than that of Model 1. The maximum values obtained for t were 6.24 and 8.65 Å in Models 1 and 2, respectively. These values were estimated to be the thickness

of water adsorption layers at the onset of the capillary condensation. According to the DBdB theory, the thicknesses of water adsorption layers at the onset of capillary condensation were predicted to be 6.04 and 7.46 Å for the pore radius of 1.38 and 1.81 nm, respectively, suggesting that the predicted thickness of the DBdB theory was slightly smaller than the calculated one of molecular simulations. In the desorption process, stable water menisci were formed on the top and bottom ends of the pore. Figures 6a and 6b show the vapor-liquid interface profiles for the desorption process in Model 1 and 2, respectively. The water meniscus changed smoothly from flat to concave with decreasing N. The curvature radius of the meniscus, ρ, was calculated by fitting a spheroid to the meniscus of water. Specifically, a local curvature radius of a spheroid, ρ1, was calculated by fitting an elliptical arc to the meniscus curve in the vapor-liquid interface profile in the r–z plane as shown in Figs. 6a and 6b within the range 0 ≤ r ≤ r0, where r0 is the inflection point of the meniscus curve. Another local curvature radius of the spheroid, ρ2, was calculated in the direction normal to both the tangential and normal directions of the fitted elliptical arc. The mean curvature radius of a surface element of a spheroid was defined as follows:

ρ ′ −1 = 2(ρ1 + ρ 2 )−1 .

(4)

The mean curvature radius of a spheroid was defined as follows:

ρ=

1 s0



s0

0

ρ ′ds ,

(5)

ACS Paragon Plus Environment

13

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 40

where ds is the line element of the fitted elliptical arc and s0 is the total length of the fitted elliptical arc. Figure 6c shows the curvature radius of menisci versus number of water molecules (ρ–N) curves. It can be seen that ρ decreases with decreasing N and approaches a constant value in both Models 1 and 2. These curves were well fitted by a power function ρ = ρ 0 + BN C , where ρ0, B, and C were parameters. The values obtained for ρ0 were 11.31 and 13.46 Å for Models 1 and 2, respectively. These values were estimated to be the curvature radii of a meniscus at the onset of the capillary evaporation.

3.2 Structural Properties 3.2.1 Water Distribution around Silanol Groups

Figure 7 shows the radial distribution functions of oxygen in water around oxygen in silanol groups (gOs*-Ow), hydrogen in water around oxygen in silanol groups (gOs*-Hw), and oxygen in water around hydrogen in silanol groups (gHs*-Ow) in Models 1 and 2. In gOs*-Ow, the first peak is located at r = 2.8 Å in both models, suggesting that this peak expresses the position of oxygen in water interacting with a silanol group due to hydrogen bonding. The second and third peaks can be seen at r = 3.7 and 4.5 Å in Model 1 and r = 3.6 and 4.8 Å in Model 2, respectively. In the bulk water, the second peak is located around 4.5 Å because of the tetrahedral ice-like structure of bulk water.29 The second peak position is about 2 2 3 times greater than that of the first peak. On the other hand, for water confined in Models 1 and 2, the second peak can be attributed to water interacting with adjacent silanol groups and the third peak is equivalent to the second peak of bulk water. Furthermore, the third peak position of Model 2 is greater than that of Model 1 and the oscillation in radial distribution function of Model 2 is much more distinct than that of Model 1, suggesting that silanol groups in Model 2 can affect the coordination structure of water over a wider range. The tetrahedrality order parameter of water was also calculated along the radial direction (see supporting information S4).30 The result shows that the water structure in the first adsorption layer in Model 2 is more highly ordered than that in Model 1. In gOs*-Hw, the first peak is located at 1.90 Å, whereas, for gHs*-Ow, the first peak can be seen at 1.75 Å. These values

ACS Paragon Plus Environment

14

Page 15 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

describe the hydrogen bond lengths between hydrogen in proton donor water – oxygen in silanol groups and oxygen in proton accepter water – hydrogen in silanol groups, respectively, suggesting that the hydrogen bond length between water – silanol groups are different depending on the types of water: namely proton donor water and proton accepter water.

3.2.2 Hydrogen Bond Network of Water Confined in Mesoporous Silica

Several methods to define ‘hydrogen bond’ have been proposed using data obtained from molecular simulations.31 In this study, we employed geometric criteria similar to those used by Luzar and Chandler.32 Specifically, the geometric criteria were as follows: (1) the distance between oxygen in a proton donor molecule (Od)–oxygen in a proton acceptor molecule (Oa), ROd-Oa, was smaller than the critical value of ROd-Oa(c), and (2) the distance between hydrogen in a proton donor molecule (Hd)–Oa, RHd-Oa, was smaller than the critical value of RHd-Oa(c). For hydrogen bonds between water molecules, ROd-Oa(c) and RHd-Oa(c) were determined by the position of the first local minimum of gOw-Ow and gHw-Ow

calculated from MD simulation of SPC/E water in the bulk liquid. Values of ROd-Oa(c) = 3.35 Å and RHdOa

(c)

= 2.45 Å were obtained which are close to results reported for a similar system.33 Hydrogen bonds

between water–silanol groups, ROd-Oa(c) and RHd-Oa(c) were similarly determined by gOs*-Ow and gOs*-Hw. Values of ROd-Oa(c) = 3.2 Å and RHd-Oa(c) = 2.6 Å were obtained. Figure 8 shows the number of hydrogen bonds per one water molecule, nHB, as a function of the number of water molecules inside the pore, N2. Subsequently, nHB was decomposed into the number of hydrogen bonds between water molecules and water-silanol groups. In the initial stage of layer adsorption (N2 = 231 (N = 300) in Model 1 and N2 = 532 (N = 600) in Model 2), nHB between water molecules and water-silanol groups was about 2.2 and 0.8, respectively, thus the total nHB was about 3.0. The value of nHB between water molecules increased monotonically with increasing N2, whereas nHB between water-silanol groups decreased monotonically with increasing N2. In the pore-filing state (N2 = 834 (N = 1600) in Model 1 and N2 = 2041 (N = 3400) in Model 2), nHB between water molecules and water-silanol groups was about 3.5 and 0.3, respectively,

ACS Paragon Plus Environment

15

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 40

giving a total value of nHB of about 3.8. Since nHB of bulk water was 3.8, nHB in the pore-filing state was close to nHB in the bulk liquid. Figures 7, S2 and S3 indicate that not only the local coordination structure of water around silanol groups but also the coordination structure and density of water near the center of the pore depends on the surface state and radius of the pore. Figure 8 shows that nHB in the initial stage of layer adsorption deviates significantly from nHB in the bulk liquid, but nHB approaches the bulk value with increasing N2. These results suggest that in the initial stage of layer adsorption, the surface state and radius of the pore have a direct affect on the coordination structure of water around silanol groups. In contrast, in the porefilling state, the hydrogen bond network of water in the pore is similar to that of the bulk liquid; however, the relative position of water molecules is not exactly the same as that of the bulk liquid due to the pore wall. This subtle difference may help confer unique properties upon water confined in mesoporous silica.

3.3 Dynamic Properties 3.3.1 Translational Motion

Figure 9 shows the spectral densities of the translations of water molecules in mesoporous silica (Model 2) from MD simulations. These spectral densities were calculated by Fourier transformation from normalized velocity autocorrelation functions of translational velocity.34 Figure 9a illustrates the effect of adsorption region (surface, middle and center) at N = 3000 when the pore is filled with water. Three concentric cylindrical shells are defined as surface, middle, and center. The radii of the boundaries between center–middle and middle-surface are r = 5.65 and 9.55 Å in Model 1 and r = 8.95 and 13.7 Å in Model 2, respectively. The first peak at 50 cm-1 can be assigned to hydrogen bond bending (OW···OW···OW) motions. The broad peak around 200 cm-1 is attributed to cluster vibrations involving combinations of hydrogen bond stretching (OW···OW) and bending (OW–HW···OW) motions.29,34,35 For surface water, the first peak is shifted by ~10 cm-1 to higher frequencies, the curve for the total shows an increase within the range of 60 to 200 cm-1, and is in good agreement with that for bulk water above 200 ACS Paragon Plus Environment

16

Page 17 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

cm-1. On the other hand, for middle and center water, the curves correspond closely to that for bulk water. These results suggest that surface water molecules are strongly bonded to silanol groups and they form a less flexible hydrogen bond network, giving rise to the observed shift of the first peak to higher frequencies. The hydrogen bond stretching (OW···OW) and bending (OW–HW···OW) motions do not undergo a significant change. Figure 9b shows the effect of the number of adsorbed water molecules: N = 800, 1500, and 3000. The first and second adsorption layers are formed at N = 800 and 1500, respectively. The N = 800 curve exhibits similar behavior to that of surface water, but the peak around 200 cm-1 is lower than that of bulk water due to the interfacial area between vapor and adsorption layers. As N increases, the curve becomes more aligned to that of bulk water. However, the peak at 50 cm-1 is shifted to a higher frequency due to surface water. Similar trends can also be observed in Model 1. Le Caër et al.36 measured the infrared spectra of water confined in well controlled pore glasses which have a diameter ranging from 8 to 320 nm. They showed that the peak of hindered translational band of confined water was shifted to higher frequencies as decreasing pore diameter. The trend of simulated spectral densities agrees with this result. However, the difference in peak frequency could not be observed between Models 1 and 2. The local self-diffusion coefficient, D, was obtained by fitting a line to the mean square displacement vs. time curve in the range of 2–5 ps (see supporting information S5). The results are summarized in Table 2. The calculated self-diffusion coefficient in the bulk water was 2.73 × 10-9 m2 s-1 which closely corresponds to results obtained for a similar system: 2.7–2.8 × 10-9 m2 s-1.33 In the pore-filling state (N = 1400 in Model 1 and N = 3000 in Model 2), D in the middle and surface is 74–75 and 30–32% of D in the center, respectively, and D of Model 2 is smaller than the corresponding value of D of Model 1. This result is supported by the radial distribution functions and density profiles of water shown in Figs. 7 and S2, which suggests that the pore wall constrained water more strongly nearer the wall and more widely in Model 2.

3.3.2 Librational Motion ACS Paragon Plus Environment

17

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 40

Figure 10 shows the spectral densities of the librations of water molecules in mesoporous silica (Model 2) from MD simulations. These spectral densities were calculated by Fourier transformation from normalized velocity autocorrelation functions of rotational velocity.34 Figure 10a shows the effect of adsorption region (surface, middle, and center) at N = 3000. The infrared spectrum of liquid water is far more complex than that of water vapor due to vibrational overtones and combinations of librations, which are due to the restrictions imposed by hydrogen bonding. The librational bands of L1 and L2 are at 395.5 and 686.3 cm-1, respectively, for liquid water at 0°C. The absorbance of L1 increases with increasing temperature, whereas the absorbance of L2 decreases but broadens with reduced wavenumber with increasing temperature.37 The spectral densities of the surface, middle, and center are similar to that of bulk water. In detail, the spectral density of the surface below 400 cm-1is slightly greater than that of bulk water, whereas the spectral density of the surface between 400–600 and 800–900 cm-1 is less than that of bulk water. This result suggests that the spectral density of the surface water has similar characteristics to bulk water at higher temperatures. Figure 10b shows the effect of the number of adsorbed water molecules: N = 800, 1500, and 3000. The N = 800 curve exhibits a similar trend to that of surface water. As N increases, the curve approaches that of bulk water. The same is also true of Model 1. The local rotational diffusion coefficient, DR, was calculated by angular velocity autocorrelation functions again using the Green-Kubo relationship. The results are summarized in Table 3. In the porefilling state (N = 1400 in Model 1 and N = 3000 in Model 2), DR in the surface is 82–83% of DR in the center, suggesting that the librational motion is not constrained to the same degree in comparison with translational motion. Furthermore, the calculated DR in bulk water (2.75 × 10-9 m2 s-1) is similar to the value of DR in the center, implying that librational motion of water near the center of the pore is similar to that of bulk water. However, DR of Model 2 is less than DR of Model 1 over all regions, which confirms the notion that the surface state and radius of pores should affect the librational motions of water inside the pore. Furthermore, DR in the surface increases with decreasing N, suggesting that librational motion is enhanced at the interface between vapor and adsorption layers. A similar ACS Paragon Plus Environment

18

Page 19 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

justification can be advanced for the observation that DR in the middle is greater than the bulk value when the first adsorption layer is formed (at N = 400 in Model 1 and N = 800 in Model 2) and DR in the center is greater than the bulk value when the second adsorption layer forms (at N = 750 in Model 1 and N = 1500 in Model 2).

The calculated dynamic properties of water show that the confining effect on librational motion was smaller than that on translational motion. This result agrees with the experimental observation reported by Faraone et al.38 They investigated the translational and rotational dynamics of water molecules in mesoporous silica using the incoherent quasielastic neutron scattering (QENS) and nuclear magnetic resonance (NMR). They reported the rotational relaxation time is less affected by confinement as compared with the translational dynamic and this result supports the idea that the rotational process is not directly influenced by constraints effective on a length scale of 20 Å.

4. CONCLUSION

GCMC and NVT-MD simulations were performed for water molecules adsorbed on two models of mesoporous silica thin films having a radius and thickness of rp = 1.38 and Lp = 5.66 nm, respectively, in Model 1 and rp = 1.81 and Lp = 7.30 nm , respectively, in Model 2. The following conclusions were drawn from this study.

Water adsorption layers or water menisci were formed in mesopores accompanying the growth or shrinkage of stable adsorption layers on the upper and lower surfaces. In both models, stable first and the second water adsorption layers were formed. However, capillary condensation commenced before the formation of a third layer. The curvature radius of a water meniscus decreased monotonically and approached a constant value, followed by the onset of capillary evaporation. The thickness of a water adsorption layer immediately prior to capillary condensation and the curvature radius of a water meniscus just before the onset of capillary evaporation were identified. ACS Paragon Plus Environment

19

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 40

In the initial stage of layer adsorption, the surface state and radius of the pore directly affected the coordination structure of water around silanol groups, whereas in the pore-filling state, the hydrogen bond network of water in the pore was similar to that in the bulk liquid. However, the pore wall affected the position of water molecules, leading to a slight discrepancy between the two states. This subtle difference may give rise to different dynamic properties of water in mesoporous silica.

In the pore-filling state, the local self-diffusion coefficient, D, in the surface was 30–32% of D in the center, whereas the local rotational diffusion coefficient, DR, in the surface was 82–83% of DR in the center, suggesting that the librational motion was not constrained to the same degree in comparison with translational motion. On the other hand, DR in the surface increased with decreasing N, suggesting that librational motion was enhanced at the interface between vapor and adsorption layers. Finally, the values of D and DR of Model 2 were smaller than those of Model 1; this leads to the conclusion that Model 2 should impose a greater constraint on both the translational and librational motions of water inside the pore despite the pore radius of Model 2 being greater than that of Model 1. The dynamic properties of confined water may be sensitive to the surface state and radius of the pores.

ACKNOWLEDGMENTS

The authors are grateful to CREST of JST for funding this work.

SUPPORTING INFORMATION Details of theoretical model based on DBdB theory, differential heat of adsorption, water density profiles along radial direction, profiles of tetrahedrality order parameter along radial direction and

ACS Paragon Plus Environment

20

Page 21 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

calculation of local self-diffusion coefficient form mean square displacement. This information is available free of charge via the Internet at http://pubs.acs.org.

REFERENCES (1) Yanagisawa, T.; Shimizu, T.; Kuroda, K.; Kato, C. The Preparation of AlkyltrimethylammoniumKanemite Complexes and Their Conversion to Microporous Materials. Bull. Chem. Soc. Jpn. 1990, 63, 988–992. (2) Kresge, C.; Leonowicz, M.; Roth, W.; Vartuli, J.; Beck, J. Ordered Mesoporous Molecular Sieves Synthesized by a Liquid-Crystal Template Mechanism. Nature 1992, 359, 710–712. (3) Inagaki, S.; Fukushima, Y. Adsorption of Water Vapor and Hydrophobicity of Ordered Mesoporous Silica, FSM-16. Microporous Mesoporous Mater. 1998, 21, 667–672. (4) Russo, P. A.; Carrott, M. M. L. R.; Padre-Eterno, A.; Carrott, P. J. M.; Ravikovich, P. I.; Neimark, A. V. Interaction of Water Vapour at 298 K with Al-MCM-41 Materials Synthesised at Room Temperature. Microporous Mesoporous Mater. 2007, 103, 82–93. (5) Endo, A.; Yamashita, K.; Yamaura, T.; Matsuoka, F.; Hihara, E.; Daiguji, H. Water Adsorption– Desorption Isotherms of Two-Dimensional Hexagonal Mesoporous Silica around Freezing Point. J. Colloid Interface Sci. 2012, 367, 409–414.

(6) Melrose J. C. Thermodynamic Aspects of Capillarity. Ind. Eng. Chem. 1968, 60, 53–70. (7) Puibasset, J.; Pellenq, R. J.-M. Water Adsorption in Disordered Mesoporous Silica (Vycor) at 300 K and 650 K: A Grand Canonical Monte Carlo Simulation Study of Hysteresis. J. Chem. Phys. 2005, 122, 094704.

ACS Paragon Plus Environment

21

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 40

(8) Bonnaud, P. A.; Coanse, B.; Pellenq, R. J.-M. Molecular Simulation of Water Confined in Nanoporous Silica. J. Phys.: Condens. Matter 2010, 22, 284110. (9) Furukawa, S.; Nishiumi, T.; Aoyama, N.; Nitta, T.; Nakano, M. Molecular Simulation of Water Confined in Nanoporous Silica Modified by Pore Surface Silylation. J. Chem. Eng. Jpn. 2005, 38, 999– 1007. (10) Shirono, K.; Daiguji, H. Molecular Simulation of the Phase Behavior of Water Confined in Silica Nanopores. J. Phys. Chem. C 2007, 111, 7938–7946. (11) de la Llave, E.; Molinero, V.; Scherlis, D. A. Water Filling of Hydrophilic Nanopores. J. Chem. Phys. 2010, 133, 034513.

(12) Gallo, P.; Ricci, M. A.; Rovere, M. Layer Analysis of the Structure of Water Confined in Vycor Glass. J. Chem. Phys. 2002, 116, 342–346. (13) Gallo, P.; Rovere, M.; Spohr, E. Glass Transition and Layering Effects in Confined Water: A Computer Simulation Study. J. Chem. Phys. 2000, 113, 11324–11335. (14) Gallo, P.; Rovere, M.; Chen, S-H. Anomalous Dynamics of Water Confined in MCM-41 at Different Hydrations. J. Phys.: Condens. Matter 2010, 22, 284102. (15) Brovchenko, I.; Geiger, A.; Oleinikova, A.; Paschek, D. Phase Coexistence and Dynamic Properties of Water in Nanopores. Eur. Phys. J. E 2003, 12, 69–76. (16) Brovchenko, I.; Geiger, A.; Oleinikova, A. Water in Nanopores. I. Coexistence Curves from Gibbs Ensemble Monte Carlo Simulations. J. Chem. Phys. 2004, 120, 1958–1972. (17) Giovambattista, N.; Rossky P. J.; Debenedetti, P. G. Effect of Temperature on the Structure and Phase Behavior of Water Confined by Hydrophobic, Hydrophilic, and Heterogeneous Surfaces. J. Phys. Chem. B 2009, 113, 13723–13734.

ACS Paragon Plus Environment

22

Page 23 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(18) Rovere, M.; Ricci, M. A.; Vellati, D.; Bruni, F. A Molecular Dynamics Simulation of Water Confined in a Cylindrical SiO2 Pore. J. Chem. Phys. 1998, 108, 9859–9867. (19) Zhang, Q.; Chan, K.; Quirke, N. Molecular Dynamics Simulation of Water Confined in a Nanopore of Amorphous Silica. Mol. Simul. 2009, 35, 1215–1223. (20) Milischuk, A. A.; Krewald, V.; Ladanyi, B. M. Water Dynamics in Silica Nanopores: the SelfIntermediate Scattering Functions. J. Chem. Phys. 2012, 136, 224704. (21) Coasne, B.; Galarneau, A.; Di Renzo, F.; Pellenq, R. J. -M. Molecular Simulation of Nitrogen Adsorption in Nanoporous Silica. Langmuir 2010, 26, 10872–10881. (22) Wright, A. F.; Lehmann, M. S. The Structure of Quartz at 25 and 590°C Determined by Neutron Diffraction. J. Solid State Chem. 1981, 36, 371–380. (23) Siboulet, B.; Coanse, B.; Dufreche, J. F.; Turq, P. Hydrophobic Transition in Porous Amorphous Silica. J. Phys. Chem. B 2011, 115, 7881–7886. (24) Tsuneyuki, S.; Tsukada, M.; Aoki, H. First-Principles Interatomic Potential of Silica Applied to Molecular Dynamics. Phys. Rev. Lett. 1988, 61, 869–872. (25) Brendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6259–6271.

(26) Allen, M. P.; Tildesley, D. J. Computer simulation of Liquids; Clarendon Press: Oxford, 1987, 78–82. (27) Nose, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 51, 511–519.

(28) Gor, G. Y.; Neimark, A. V. Adsorption-Induced Deformation of Mesoporous Solids. Langmuir

2010, 26, 13021–13027. ACS Paragon Plus Environment

23

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 40

(29) Rahman, A.; Stillinger, F. H. Molecular Dynamics Study of Liquid Water. J. Chem. Phys. 1971, 55, 3336–3359.

(30) Alabarse, F.; Haines, J.; Cambon, O.; Levelut, C.; Bourgogne, D.; Haidoux, A.; Granier, D.; Coasne, B. Freezing of Water Confined at the Nanoscale. Phys. Rev. Lett. 2012, 109, 035701. (31) Marti, J. Analysis of the Hydrogen Bonding and Vibrational Spectra of Supercritical Model Water by Molecular Dynamics Simulations. J. Chem. Phys. 1999, 110, 6876–6886. (32) Luzar, A.; Chandler, D. Structure and Hydrogen Bond Dynamics of Water–Dimethyl Sulfoxide Mixtures by Computer Simulations. J. Chem. Phys. 1993, 98, 8160–8173. (33) Mark, P.; Nilsson, L. Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. J. Phys. Chem. A 2001, 105, 9954–9960. (34) Heinzinger, K. Computer Simulations of Aqueous Electrolyte Solutions. Physica 1985, 131B, 196–216. (35) Marti, J.; Padro, J. A.; Guardia, E. Molecular Dynamics Simulation of Liquid Water along the Coexistence Curve: Hydrogen Bonds and Vibrational Spectra. J. Chem. Phys. 1996, 105, 639–649. (36) Le Caër S.; Pin, S.; Esnouf, S.; Raffy, Q.;Renault, J. Ph.; Brubach, J.-B.; Creff, G.; Roy, P. A Trapped Water Network in Nanoporous Material: the Role of Interfaces. Phys. Chem. Chem. Phys. 2011, 13, 17658–17666.

(37) Zelsmann, H. R. Temperature Dependence of the Optical Constants for Liquid H2O and D2O in the Far IR Region. J. Mol. Struct. 1995, 350, 95–114. (38) Faraone, A.; Liu, L.; Mou, C.-Y.; Shih, P.-C.; Copley, J. R. D. Chen, S.-H. Translational and Rotational Dynamics of Water in Mesoporous Silica Materials: MCM-41-S and MCM-48-S. J. Chem. Phys. 2003, 119, 3963–3971.

ACS Paragon Plus Environment

24

Page 25 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 1 Parameters of potential functions.10, 24 Atom

Z

ai (Å)

bi (Å)

ci (eV1/2 Å3)

Os, Os*

-1.2

2.0474

0.17566

14.657

Sis,Sis*

+2.4

0.8688

0.03285

4.828

Hs*

+0.6

0.0

0.0

0.0

Ow

-0.8476

1.694

0.1179

5.21

Hw

+0.4238

ACS Paragon Plus Environment

25

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 40

Table 2 Calculated local self-diffusion coefficients of water adsorbed on mesoporous silica thin films. The unit is × 10-9 m2 s-1. Model 1

Model 2

N = 400

N = 750

N = 1400

N = 800

N = 1400

N = 3000

Center

-

-

2.50

-

-

2.28

Middle

-

2.13

1.86

-

1.83

1.71

Surface

0.74

0.99

0.81

0.64

0.76

0.68

ACS Paragon Plus Environment

26

Page 27 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 3 Calculated local rotational diffusion coefficients of water adsorbed on mesoporous silica thin films. The unit is × 10-9 m2 s-1. Model 1

Model 2

N = 400

N = 750

N = 1400

N = 800

N = 1500

N = 3000

Total

3.02

3.02

2.68

2.57

2.70

2.43

Center

-

8.49

2.99

-

5.45

2.62

Middle

5.32

3.75

2.75

5.37

3.12

2.54

Surface

2.87

2.50

2.48

2.28

2.21

2.16

ACS Paragon Plus Environment

27

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 40

Figure Captions Figure 1. Top view (left) and cross-sectional view (right) of the atomic model of a mesoporous silica thin film. The pore radius of rp = 1.38 nm and the film thickness of Lp = 5.66 nm (Model 1). Figure 2. Snapshots of NVT-MD simulations of water adsorbed on a model mesoporous silica thin film. The number of water molecules of N = 500 (left) and 1200 (right) (Model 1). Figure 3. (a) N1–N2 curves, (b) chemical potential of water vs. the number of water molecules (µ–N) curves, (c) chemical potential of water vs. the area number density of water molecules on the upper or lower surface (µ–n1) curves, (d) chemical potential of water vs. the number density of water molecules in the pore (µ–n2) curves and the theoretical model base on the DBdB theory. Figure 4. Contour plots of water density as a function of N in Models 1 and 2. Figure 5. Profiles of the vapor-liquid interface in the layer adsorption process in (a) Model 1 and (b) Model 2, (c) average thickness of adsorption layer vs. number of water molecules (t–N) curves. Figure 6. Profiles of the vapor-liquid interface in the desorption process in (a) Model 1 and (b) Model 2, (c) curvature radius of meniscus vs. number of water molecules (ρ–N) curves. Figure 7. Radial distribution functions of oxygen in water (Ow) around oxygen in silanol group (Os*), and hydrogen in water (Hw) around oxygen in silanol group (Os*), and oxygen in water (Ow) around hydrogen in silanol group (Hs*) in (a) Model 1 and (b) Model 2. Figure 8. The number of hydrogen bonds per one water molecule, nHB, as a function of the number of water molecules inside the pore, N2. nHB was decomposed into the number of hydrogen bonds between water molecules and water-silanol group’s water molecule.

ACS Paragon Plus Environment

28

Page 29 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 9. Spectral densities of the translations of water molecules in mesoporous silica (Model 2) from MD simulations. (a) Effect of adsorption region (surface, middle and center) at N = 3000 and (b) effect of the number of adsorbed water molecules (N = 800, 1500, and 3000). Figure 10. Spectral densities of the librations of water molecules in mesoporous silica (Model 2) from MD simulations. (a) Effect of adsorption region (surface, middle, and center) at N = 3000 and (b) effect of the number of adsorbed water molecules (N = 800, 1500, and 3000).

ACS Paragon Plus Environment

29

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 40

Figure 1. Top view (left) and cross-sectional view (right) of the atomic model of a mesoporous silica thin film. Pore radius of rp = 1.38 nm and film thickness of Lp = 5.66 nm (Model 1).

ACS Paragon Plus Environment

30

Page 31 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 2. Snapshots of NVT-MD simulations of water adsorbed on a model mesoporous silica thin film. Number of water molecules of N = 500 (left) and 1200 (right) (Model 1).

ACS Paragon Plus Environment

31

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 40

Figure 3. (a) N1–N2 curves, (b) chemical potential of water vs. the number of water molecules (µ–N) curves, (c) chemical potential of water vs. the area number density of water molecules on the upper or lower surface (µ–n1) curves, and (d) chemical potential of water vs. the number density of water molecules in the pore (µ–n2) curves and the theoretical model base on the DBdB theory.

ACS Paragon Plus Environment

32

Page 33 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 4. Contour plots of water density as a function of N in Models 1 and 2.

ACS Paragon Plus Environment

33

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 40

Figure 5. Profiles of the vapor-liquid interface in the layer adsorption process in (a) Model 1 and (b) Model 2, (c) average thickness of adsorption layer vs. number of water molecules (t–N) curves.

ACS Paragon Plus Environment

34

Page 35 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 6. Profiles of the vapor-liquid interface in the desorption process in (a) Model 1 and (b) Model 2, (c) curvature radius of meniscus vs. number of water molecules (ρ–N) curves.

ACS Paragon Plus Environment

35

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 40

Figure 7. Radial distribution functions of oxygen in water (Ow) around oxygen in silanol group (Os*), and hydrogen in water (Hw) around oxygen in silanol group (Os*), and oxygen in water (Ow) around hydrogen in silanol group (Hs*) in (a) Model 1 and (b) Model 2.

ACS Paragon Plus Environment

36

Page 37 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 8. The number of hydrogen bonds per one water molecule, nHB, as a function of the number of water molecules inside the pore, N2. nHB was decomposed into the number of hydrogen bonds between water molecules and water-silanol group’s water molecules.

ACS Paragon Plus Environment

37

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 40

Figure 9. Spectral densities of the translations of water molecules in mesoporous silica (Model 2) from MD simulations. (a) Effect of adsorption region (surface, middle and center) at N = 3000 and (b) effect of the number of adsorbed water molecules (N = 800, 1500, and 3000).

ACS Paragon Plus Environment

38

Page 39 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 10. Spectral densities of the librations of water molecules in mesoporous silica (Model 2) from MD simulations. (a) Effect of adsorption region (surface, middle and center) at N = 3000 and (b) effect of the number of adsorbed water molecules (N = 800, 1500, and 3000).

ACS Paragon Plus Environment

39

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 40

TABLE OF CONTENTS

ACS Paragon Plus Environment

40