Article pubs.acs.org/JPCC
Deconstructing the Confinement Effect upon the Organization and Dynamics of Water in Hydrophobic Nanoporous Materials: Lessons Learned from Zeolites Tiecheng Zhou,*,† Peng Bai,‡ J. Ilja Siepmann,‡ and Aurora E. Clark*,† †
Department of Chemistry and Materials Science and Engineering Program, Washington State University, Pullman, Washington 99164, United States ‡ Department of Chemistry and Chemical Theory Center and Department of Chemical Engineering and Materials Science, University of Minnesota, 207 Pleasant Street SE, Minneapolis, Minnesota 55455, United States S Supporting Information *
ABSTRACT: The properties of confined water are relevant to many chemical, geological, and biological phenomena, where they underpin essential changes to molecular scale reactivity− perturbing both the energetic and configurational landscape. Though much prior literature has focused on hydrophilic confinement, the hydrophobic confinement of water is less well understood. Here, we use molecular dynamics simulations to investigate the structures and dynamics of water in hydrophobic all-silica zeolites that have sequentially smaller pore dimensions. Of special interest is the role that pure geometric restriction imparts, relative to the rugged potential energy landscape for water interacting with the atomistic pore surface. These two effects were studied via the hydrogen bond dynamics, specifically the rates and mechanisms of hydrogen bond breakage and formation. Measuring the dynamic features as a function of scaling the water:zeolite interaction energy revealed that geometric restriction is responsible for 67%−86% of the total perturbations to water upon confinement in MFI (depending on the property) while the water:surface interactions are responsible for 14%− 33%. The relative magnitude of the interaction of water with the pore surface was confirmed by second order Møller−Plesset perturbation theory. Thus, in a highly confined environment, the weak water-surface interaction should not be neglectedeven in hydrophobic adsorbents to which zeolites and other materials like carbon nanotubes belong. chemical nature of the water−surface interaction.23,37,40 In addition to changes in translational dynamics, orientational dynamics are also perturbed. Using ultrafast infrared (IR) pump−probe spectroscopy, the orientational relaxation dynamics measured from the anisotropy decay have been observed to be much slower than bulk water for confined water in reverse micelles.3,5 Using MD simulations, the reorientation dynamics (measured from the reorientation time−correlation function) have been found to be slower for confined water in hydrophilic silica pores39,45 and hydrophobic all-silica Linde type A (LTA) zeolite46 (which has a cage diameter of 13 Å). The confinement effect has two major components, consisting of geometric restriction and the water:surface interaction energy.3,29,35,47 Ultrafast IR pump−probe spectroscopy experiments have indicated that geometric restriction plays the dominant role in governing the dynamics of water in hydrophilic reverse micelles.3,29 In contrast, MD simulations have indicated that instead the chemical nature of the surface atoms (and thus the interaction energy) has a significant effect on the dynamics of the confined water in hydrophilic silica
1. INTRODUCTION Confined water is present in a variety of environments, from reverse micelles or protein cavities,1−5 to water transported through membranes,6−8 rocks,9−11 or other porous materials.12−17 Additionally, it exhibits unusual properties that are not observed in the bulk, including the observation of a nonfreezable water layer inside mesoporous silica,17 a liquidto-amorphous phase transition at low temperature for water in silica nanopores and zeolites,9,11,12,18 and a high proton conductivity for water in porous materials.10,19−21 A variety of experimental and theoretical methods have been employed to study confined water, including quasi-elastic neutron scattering (QENS)12,22−24 and inelastic neutron scattering (INS),19,25 Raman scattering,17 nuclear magnetic resonance (NMR),26,27 nuclear spin echo measurements,28 ultrafast infrared (IR) spectroscopy,3,5,6,29 optical Kerr effect spectroscopy, 3 0 − 3 2 molecular dynamics simulations (MD),9,11,33−39 and Monte Carlo simulations (MC).40−42 As anticipated, confined water undergoes significant structural changes (destruction of the 3-dimensional hydrogen bond (HB) network accompanied by loss of hydrogen bonds for each molecule)41−44 commensurate with a decrease in the diffusion coefficient,12,23,24,27 which is related to the distance between individual water molecules and the confining surface and the © XXXX American Chemical Society
Received: May 23, 2017 Revised: August 12, 2017 Published: September 19, 2017 A
DOI: 10.1021/acs.jpcc.7b04991 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C nanopores.35,47 However, deconstructing the influence of the geometric restriction versus the water:surface interaction energy is challenging in hydrophilic systems due to the strong and favorable water:surface interactions. Instead, more may be learned about this issue by studying water under hydrophobic confinement. Though less studied, water confined in hydrophobic environments does exhibit some similarities to the hydrophilic confinement system. A decreased number of hydrogen bonds (HBs) has been found both for water in a hydrophobic slit34,43 and for water in the hydrophobic all-silica zeolites LTA and MFI,41,42 a feature that also influences the heats of water intrusion.48 Shorter HB lifetimes have been found from the rate of HB breaking kinetics of the intermittent HB correlation function for water in carbon nanotubes using MD simulations.49 The HB lifetime, determined from the decay of the continuous HB population correlation function, has also been shown to be shorter than bulk water when confined in graphene sheets.44 Similarly, the HB lifetime, calculated from the HB persistence distribution, has been found to be smaller for water within hydrophobic slits using MD simulations.34,43 These findings suggest faster HB dynamics for water under hydrophobic confinement. However, the relaxation of the intermittent HB correlation function of water is slower within carbon nanotubes,49 and the reorientation time scale obtained both from the reorientation time-correlation function and from the extended jump model (EJM) is longer for water in hydrophobic LTA zeolites.38 In combination, these studies indicate that more than one component of the HB dynamics may be affected by hydrophobic confinement, and that these components may be manifested in different ways as a result of the underlying and perhaps competing forces that result from geometric restriction versus the water:surface interaction energy. This work attempts to decompose and quantify the importance of these two components, using MD simulations and intermolecular network theory (INT)50 to study the structure and dynamics of water adsorbed within the hydrophobic all-silica zeolites FAU and MFI. These adsorbents exhibit a well-defined difference in pore size with a largest included sphere diameter of 11.2 Å (FAU (Faujasite) and 6.4 Å (MFI, also known as (ZSM-5). The structure of the HB network is altered significantly within this series, such that the extended 3-dimensional HB network of bulk water is transformed into small water clusters in FAU, while 1dimensional H-bonded water chains are formed in the channels of MFI. Analysis of all hydrogen bonds reveals that confinement shortens the hydrogen bond (HB) lifetime and lengthens the non-hydrogen-bond (NHB) lifetimes compared to the bulk. In addition, the HB dynamics of confined water have been examined by a detailed analysis of two types of HB breakage/ reformation motions (see Computational Methods). Interestingly HB fluctuation motions (those that do not lead to changes in hydrogen bond partners) are observed to increase in frequency, while HB reorientation (with switching of HB partners) decreases in frequency for water adsorbed in FAU and MFI. These dynamic features are compared as a function of scaling of the water:zeolite interaction energy, making it apparent that geometric restriction in MFI plays a dominant role in determining both water HB structural features and dynamic properties (responsible for 67%−86% of the total perturbations to water upon confinement, depending on the property). Yet importantly, even in these hydrophobic systemswhere the water:zeolite interaction is weakthe
interaction energy also has a considerable influence on the properties of confined water (responsible for 14%−33% of the total effect, depending on the property). In contrast to typical computational approaches, this implies that the weak water:surface interaction energy should not be neglected when studying the structures and dynamics of water under hydrophobic confinement.
2. COMPUTATIONAL METHODS System Information. Adsorbed water confined in all-silica zeolites FAU and MFI (which only consist of silicon and oxygen atoms) was investigated. The zeolite FAU has a cagewindow structure, where its supercage can accommodate a sphere of diameter 11.2 Å and its windows are formed by a ring containing 12 oxygen atoms with a cross section of 7.4 Å. The zeolite MFI has a two-dimension channel-like structure, where its channels are formed by 10-member rings that allow for diffusion of a sphere with a maximum diameter of 4.7 Å, while the channel intersection can accommodate a sphere of diameter 6.4 Å (see Figure 1). The information on the simulated systems is summarized in Table 1. The crystal structures of the all-silica zeolites (FAU, faujasite; and MFI, ZSM-5) were obtained from the International Zeolite Association (www.iza-online.org). The initial structures of the fully hydrated zeolites for molecular dynamics simulations were obtained from Gibbs ensemble Monte Carlo simulations. Gibbs Ensemble Monte Carlo Simulations. Configurational-bias Monte Carlo simulations51 in the NpT−Gibbs ensemble52,53 (CB-GEMC) were carried out using the in-house MC simulation software package MCCCS-MN (Monte Carlo for Complex Chemical Systems - Minnesota)54 to calculate water adsorption isotherms and obtain equilibrium configurations at saturation loadings. The zeolite simulation box consisted of 2 × 2 × 2 (FAU) or 2 × 2 × 3 (MFI) unit cells, and a total of 2500 (FAU) or 1500 (MFI) water molecules were initially placed in the bulk liquid box and allowed to transfer with the zeolite box via a fictitious gas phase. The simulations were equilibrated at 298 K for about 100 000 Monte Carlo cycles, and production periods consisted of 50 000 cycles (pressures noted in Figure 3). Each cycle consisted of N randomly selected moves (N being the number of water molecules in the system). The fraction of volume moves (applied only to the bulk liquid box) and that of configurational-bias particle transfer moves were each adjusted to yield approximately one accepted move per cycle, and the remainder of the moves were equally divided between translations and rotations for water molecules, wherein the zeolite framework atoms were kept rigid. Water molecules were represented by the TIP4P55 model, and their interactions with zeolite framework atoms were described by the TraPPE-zeo56 force field (see in Table 2), where the parameters for unlike Lennard-Jones interactions were determined from the Lorentz−Berthelot combining rules.57 Coulomb interactions were handled using the Ewald summation technique.58 The TIP4P model was employed as it yields only a relatively small overprediction of the vapor pressure of neat water (by about 30%), which is a significant improvement over later models like TIP4P/2005.59 Molecular Dynamics Simulations. LAMMPS60 (version 10 August, 2015) was used to create an equilibrated ensemble of time-dependent configurations to study the structure and dynamics of the confined water in zeolites FAU and MFI at the saturation loading. The positions of the zeolite atoms were B
DOI: 10.1021/acs.jpcc.7b04991 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C
Table 1. System Information Regarding the MD Simulations system
number of waters
simulation box (Å)
bulk water in FAU water in MFI
900 1773 463
30.000 × 30.000 × 30.000 48.376 × 48.376 × 48.376 40.044 × 39.798 × 40.149
Table 2. Force Field Parameters Used in MD Simulations force field TIP4P
traPPE-zeo
atoms
q (|e|)
εii (kcal/mol)
σii (Å)
Ow−Mw (Å)
Ow Hw Mw Oz Si
− 0.520 −1.040 −0.75 1.50
0.1550 − − 0.1053 0.0437
3.1536 − − 3.3 2.3
0.15
TIP4P55 force field for water (Table 2) were employed in the MD runs. The Shake algorithm61 was used to constrain the bond length and the bond angle for the TIP4P water model. The particle−particle particle-mesh solver was used to calculate the long-range Coulomb interaction with a cutoff of 12.8 Å and a relative error in force of 0.001. The repulsive-dispersive interactions were represented by 12−6 Lennard-Jones potentials truncated at 12.8 Å. The Lorentz−Berthelot combining rules57 were used for the unlike Lennard-Jones parameters. The simulated systems were first equilibrated in the NVT ensemble at T = 298 K for 1 ns to reach the equilibrium. The production runs were then performed in the NVE ensemble for 3 ns in total with configurations saved every 10 fs. The total 3 ns trajectories were divided into three sets of data (1 ns each) to estimate the statistical error bars in the results. The same simulation protocols were applied to the bulk water system for comparison. Hydrogen Bond (HB) Network Analysis. Intermolecular network theory (INT) analyses50 were performed using the ChemNetworks software,62 and used to quantify the structural and dynamical properties of the HB networks for bulk water, and confined water in FAU and MFI. A geometric definition for hydrogen bonds was employed wherein two water molecules were considered to be H-bonded if their intermolecular O··H distance was less than 2.50 Å and also the intermolecular ∠OH··O was larger than 145°. The whole trajectories from the MD simulation were converted into the time-dependent HB networks that were then data-mined to characterize all the HB structural and dynamical information. The water structure was studied in terms of the radial distribution functions (RDFs), the hydrogen-bond number distribution, the average number of HBs per water molecule, and the distribution of geodesic path lengths. The hydrogenbond number distribution is the probability of a water molecule that have 0−5 hydrogen bonds. The average number of HBs per water is calculated as the weighted average of this distribution. Geodesic paths are defined as the shortest unbroken path of hydrogen bonds that connects any pairs of water molecules. It is thus an index that measures the extension of the HB network. The water HB dynamics were studied in terms of the frequency of the HB breakage/reformation events (see Figure 2), and the corresponding contiguous hydrogen bond (HB) and the non-hydrogen-bond (NHB) lifetimes. To ensure that these quantities are less dependent upon the exact geometric criterion associated with the definition of a hydrogen bond, we have employed an algorithm63 that examines the history of
Figure 1. Snapshots of the water structures for (a) bulk water, (b) water adsorbed in the FAU pore viewed along the [111] direction, and (c) water adsorbed in the MFI channel viewed along the [010] direction. The oxygen and hydrogen atoms of water molecules are shown as red and white spheres, and the dotted blue lines indicate hydrogen bonds between water molecules. The zeolite atoms are colored in yellow and red and one of the pore windows is highlighted with thick lines. Orthogonal views of the same structures are presented in Figure S4, with a side view for FAU and MFI (the side view of MFI clearly reveals the 1-dimensional chain structure of water in MFI).
fixed, as is customary for small, nontight-fitting sorbate molecules. The TraPPE-zeo56 force field for zeolites and the C
DOI: 10.1021/acs.jpcc.7b04991 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C
Figure 2. Two types of HB breakage/reformation motions.
molecular dynamics simulations for confined water in MFI. Figure S5 plots the total time for each rotational event (beginning with the snapshot prior to HB breakage and ending with the snapshot where the new HB is formed), where it is observed that many rotational events occur on a short time scale (50% of all events last
