Deconstructing Hydrogen-Bond Networks in Confined Nanoporous

Aug 7, 2014 - Essential topological indices of the hydrogen-bond networks of water, methanol, ethanol, and their binary mixtures adsorbed in microporo...
0 downloads 11 Views 2MB Size
Article pubs.acs.org/JPCC

Deconstructing Hydrogen-Bond Networks in Confined Nanoporous Materials: Implications for Alcohol−Water Separation Chun-Hung Wang,† Peng Bai,‡ J. Ilja Siepmann,*,‡,§ and Aurora E. Clark*,† †

Department of Chemistry and the Materials Science and Engineering Program, Washington State University, Pullman, Washington 99164, United States ‡ Department of Chemical Engineering and Materials Science, University of Minnesota, 421 Washington Avenue SE, Minneapolis, Minnesota 55455, United States § Department of Chemistry and Chemical Theory Center, University of Minnesota, 207 Pleasant Street SE, Minneapolis, Minnesota 55455, United States S Supporting Information *

ABSTRACT: Essential topological indices of the hydrogen-bond networks of water, methanol, ethanol, and their binary mixtures adsorbed in microporous silicalite-1 (a hydrophobic zeolite with potential application for biofuel processing) are analyzed and compared to their bulk liquid counterparts. These include the geodesic distribution (the shortest H-bond pathways between molecular vertices), the average length, the geodesic index, the orientation and distance of the adsorbate to the interior of the zeolite, and the sorbate−sorbate and sorbate−sorbent distributions of Hbonds. In combination, they describe how the H-bond networks are altered when going from the bulk to the confined silicalite-1 environment. The speciation of the adsorbed compounds is quantified in terms of their network connectivity, revealing that pure water has a high probability of forming long, contiguous H-bonded chains in silicalite-1 at high loading, while alcohols form small dimeric/trimeric clusters. The extent to which the H-bond network of binary water−alcohol systems is altered relative to either unary system is quantified, demonstrating an enhanced interconnectivity that is reflected in the tendency of individual H2O molecules to become co-adsorbed with alcohol clusters in the zeolite framework. Selectivity for the alcohol over water diminishes with increasing alcohol loading as the H-bonded clusters serve as favorable adsorption sites for H2O.



INTRODUCTION Graph theoretical algorithms have recently been developed with the aim of characterizing chemical networks defined by intermolecular interactions.1 Hydrogen-bond networks are perhaps the most famous class of chemical networks, and numerous studies have been dedicated to the structural and dynamic features of bulk water.2−4 Somewhat less attention has traditionally been paid to the network properties of other Hbonding solvents5 or solvent mixtures,6 and characterization of the H-bond network properties at interfaces1,7 and in confined environments8−10 is an emerging topic of high interest. These issues are relevant to the industrial and energy sectors, where microporous materials (i.e., frameworks containing pores with free diameters of less than 2 nm)11 are becoming increasingly common adsorbents for separation of complex chemical mixtures, often involving H-bonding solvents. Biofuel purification, for example, relies on separation of bioalcohol from water, an inherently difficult task for distillation because of the relatively similar saturated vapor pressures for water and ethanol and the presence of an azeotrope. The separation of alcohol from water may also be achieved by selective adsorption or pervaporation using a microporous materiala zeolite, for © 2014 American Chemical Society

instancewhose specificity relies upon the physical interactions and organization of the sorbate molecules inside the zeolite framework. The specificity is governed by the hydrophobic interior, pore size distribution, pore volume, and potentially the alteration of the H-bond network formed by the sorbates. Much prior work has focused upon the amount of adsorbate in the zeolite as a function of the solution-phase concentration or partial vapor pressure of the adsorbates,12−26 the mixture selectivity and partition coefficient,25,26 and occasionally analysis of the adsorbate clustering and H-bond properties inside the framework.23,25,27,28 Generally, only a basic understanding has been obtained regarding the H-bond network inside the confined environment of the zeolite; thus, there has been little discussion regarding how it is perturbed during the adsorption process. Toward this end, the present work utilizes a series of topological descriptors for the quantification of adsorbed water and simple alcohols (methanol and ethanol) in terms of their Received: March 22, 2014 Revised: August 5, 2014 Published: August 7, 2014 19723

dx.doi.org/10.1021/jp502867v | J. Phys. Chem. C 2014, 118, 19723−19732

The Journal of Physical Chemistry C

Article

molecules per unit cell reflect the different molecular volumes of these three compounds. Similarly, the solution-phase isotherms are shifted to much lower concentrations for the water−ethanol mixture than for the water−methanol mixtures, but both systems exhibit at a maximum in water loading at intermediate concentrations. The three solution-phase concentrations analyzed in this article reflect conditions (i) where relatively few alcohol molecules are adsorbed and the adsorption selectivities are very high, (ii) near the maximum in water loading, and (iii) when most water molecules are displaced by alcohol molecules but the adsorption selectivities are much lower.23 Molecular Dynamics. The solution-phase concentrations that describe the unary and binary systems investigated using molecular dynamics simulations are presented in Table 1, using the same force fields described above. Simulation boxes were constructed using Packmol.36 Each system was initially brought to a temperature of 350 K and allowed to evolve using the NVT ensemble for 2 ns using a 1 fs time step and the Nosé−Hoover thermostat.37−39 This was followed by cooling in 10 K increments until a temperature of 300 K was achieved. Equilibration using the NPT ensemble was performed at 300 K for 5 ns using the Nosé− Hoover thermostat and barostat,40,41 followed by a production run in the NVE ensemble for 2 ns. An Ewald summation with a relative accuracy of 10−6 was used for all simulations with a 10 Å cutoff. All the simulations were performed with the DL_Poly_3 software.42 Graph Theory Analyses. The water, methanol, and ethanol solutions have been characterized by their H-bond network topologies using the ChemNetworks software.1 For each system and frame of data, the H-bond network was represented as a graph where the entire molecule becomes a vertex (at the location of the oxygen site) and edges between vertices represent a H-bond that satisfies a set of geometric criteria. Geometric criteria are often employed to define Hbonds, and prior work has shown that an angle criterion between 140 and 150° and distance cutoff of 2.5−2.7 Å yield very similar topological properties.1 Further, these geometric criteria have been demonstrated to yield an interpretation equivalent to that obtained from energetic criteria, provided the geometric parameters are not too restrictive.43 For the sake of simplicity, the same criteria were used to define the interactions H2O···H2O, H2O···MeOH, H2O···EtOH, MeOH···MeOH, and EtOH···EtOH. As shown in Figure S1, an O···H bond distance cutoff of 2.5 Å and a bond angle cutoff of 150° were employed. Topological analyses of the graph focused primarily upon the H-bond distribution per molecule for each type of H-bond that can be formed and the geodesics (gd’s, the shortest contiguous paths between two vertices) of the H-bond networks obtained using the Floyd−Warshall algorithm,44,45 both of which can be used to identify clustering or chain formation that may occur within the confined environment of the zeolite. The N×N geodesic distance matrix has entries that are the lengths of gd’s that connect every pair of vertices in the system. For each gd that is calculated, ChemNetworks maintains a connectivity list such that all vertices that are part of the geodesic pathway are known. Inspection of the connectivity lists reveals that many of the gd’s in the geodesic distance matrix are redundant (e.g., a geodesic that connects vertices 1-2-3 also contains the 1-2 geodesic), and thus these “sub-geodesics” have been removed to create what is termed the isolated geodesic distance matrix. This sparse matrix is a better reflection of the extension of the

H-bond network structure within the zeolite known as silicalite1. Particular emphasis lies in the comparison and contrast of the analogous descriptors from the bulk solutions. These descriptors include the geodesic distribution (the shortest Hbond pathways between molecular vertices), average length, and range, as well as the geodesic index, ηgd, which is a measure of the interconnectivity within the graph, the orientation and distance of the adsorbate to the interior of the zeolite, and finally the sorbate−sorbate and sorbate−sorbent distribution of H-bonds per molecule.



COMPUTATIONAL DETAILS For this work, trajectories obtained from Monte Carlo and molecular dynamics simulations are analyzed. Configurational-Bias Monte Carlo in the Gibbs Ensemble (CBMC-GE). The solution-phase adsorption of water, methanol (MeOH), ethanol (EtOH), and their binary mixtures was calculated using configurational-bias Monte Carlo simulations29,30 in the NpT-version of the Gibbs ensemble,31,32 in which both the adsorbed and the bulk solution phases were modeled explicitly using separate simulation boxes. These simulations were performed using a zeolite simulation box containing 12 unit cells and a total of 1100 molecules for the entire system with overall compositions ranging from 64:1 to 1:64 water:alcohol, at 298 K and an external pressure of 1 atm. The interactions of alcohols, water, and silicalite-1 are described by the united-atom version of the TraPPE (transferable potentials for phase equilibria) force field,33 the TIP4P (transferable intermolecular potential−four point) model for water,34 and the TraPPE-zeo force field.35 For complete details regarding the Monte Carlo simulations, the reader is referred to a previous paper.23 Figure 1 shows the adsorption isotherms for

Figure 1. Adsorption isotherms at T = 298 K.23 (Left) Unary adsorption for water, methanol, and ethanol. (Right) Binary adsorption from aqueous solutions of methanol and ethanol. The filled symbols denote the state points analyzed in this work.

neat water, methanol, and ethanol and for solution phases containing water + (methanol or ethanol). Due to the hydrophobic nature of silicalite-1 and the different strengths of dispersive interactions of the framework with these three compounds, very different external pressures are required to achieve saturation loading, and the numbers of adsorbed 19724

dx.doi.org/10.1021/jp502867v | J. Phys. Chem. C 2014, 118, 19723−19732

The Journal of Physical Chemistry C

Article

Table 1. Mole Fraction (χH2O) and Molecular Composition of Unary and Binary Systems Unary no. of molecules ρ (g/cm3)

H2O

MeOH

EtOH

1460 0.997

651 0.791

452 0.789

Binary H2O/MeOH

H2O/EtOH

χH2O

0.25

0.50

0.75

χH2O

0.25

0.50

0.75

no. of H2O no. of MeOH ρ (g/cm3)

864 288 0.940

469 469 0.890

194 582 0.839

no. of H2O no. of EtOH ρ (g/cm3)

744 248 0.941

366 366 0.889

142 426 0.841

unbroken H-bond pathways within the graph (see Supporting Information for a complete discussion). Only the isolated geodesic distance matrix has been used for the analyses presented here; all subsequent references to “geodesics” imply the isolated gd’s and their respective connectivity lists which are used to determine the geodesic index (vide inf ra). Note that descriptors based upon the gd’s are sensitive to the system size; however, the size dependence is smaller than the differences in descriptor values emphasized in this work, and thus does not impact the conclusions herein (see Supporting Information).



RESULTS AND DISCUSSION Hydrogen-Bond Networks of Water in the Bulk and Absorbed within Silicalite-1. The average distribution of Hbonds for neat water is shown in Figure 2A. In the liquid phase, 25%, 39%, and 26% of water molecules are involved in four, three, and two H-bonds, respectively, with only small fractions either being isolated (not having any H-bonds) or possessing five H-bonds. The current definition of a H-bond leads to a slightly smaller average number of H-bonds per water molecule than would be obtained from the simple integration of the nonbonded radial distribution function (see Figure S2 in Supporting Information); however, trends in the data will be consistent across the conditions studied. Geodesic analysis is the study of the shortest edge pathways (geodesics, or gd’s) between verticesin this case, the shortest contiguous H-bond path between water molecules. The distribution of isolated gd’s (see Supporting Information) is indicative of the extension of the H-bond network from any reference water molecule in the system. When combined with the actual Cartesian distance between vertices, the length scale of interconnectivity within the network can be elucidated. The distribution of gd’s in the bulk liquid is essentially normal (or Gaussian) in shape, with a peak at 13, an average path length of 12 (see Figure 2B and Table 2), and an average Euclidean distance of ∼20 Å. Thus, on average, each water molecule is connected to other water molecules that are ∼20 Å away through 12 contiguous Hbonds. The upper bound in the geodesic distribution is 25, yielding a range in the distribution of 24 (see Table 2). As illustrated in Figure 2C, there is a large spread in the distances that the different gd’s adopt. The extension over Euclidean space increases with the number of contiguous H-bonds and is reflective of the complexity of the H-bond network, where the shortest H-bond pathway between waters can be circular, an extended chain, or a variety of three-dimensional shapes in between. Another metric to assess the interconnectivity within the network is the geodesic index:46

Figure 2. (A) H-bond distribution of bulk liquid TIP4P water, and water adsorbed in silicalite-1 at varying loadings. (B) Distribution of geodesic pathways (contiguous H-bond paths). (C) Euclidean distance average (symbol) and range (bar) between vertex end points for bulk water (red) and water inside silicalite-1 at a loading of 38.5 water molecules per unit cell (blue) for geodesics up to length 13. M ⎧ σ ̅ m) ⎫ ⎪∑ ⎪ 1 gd( ⎬ ηgd = 100 × ⎨ ⎪ N × M ⎪ ⎭ ⎩

19725

(1)

dx.doi.org/10.1021/jp502867v | J. Phys. Chem. C 2014, 118, 19723−19732

The Journal of Physical Chemistry C

Article

Table 2. Essential Topological Indices of the Unary and Binary Liquids (Subscript “liq”) and Adsorbed Phases in Silicalite-1 (Subscript “zeo”)a ⟨nads/unit⟩

system H2Oliq H2Ozeo @ pext = 252 MPa MeOHliq MeOHzeo @ pext = 12.9 kPa MeOHzeo @ pext = 2.5 kPa H2O···H2O in 3:1 H2O:MeOHliq MeOH···MeOH in 3:1 H2O:MeOHliq entire network in 3:1 H2O:MeOHliq H2O···H2O in 1:3 H2O:MeOHliq MeOH···MeOH in 1:3 H2O:MeOHliq entire network in 1:3 H2O:MeOHliq EtOHliq EtOHzeo @ pext = 7.2 kPa H2O···H2O in 3:1 H2O:EtOHliq EtOH···EtOH in 3:1 H2O:EtOHliq entire network in 3:1 H2O:EtOHliq H2O···H2O in 1:3 H2O:EtOHliq EtOH···EtOH in 1:3 H2O:EtOHliq entire network in 1:3 H2O:EtOHliq

38.5 17.9 14.7

14.1

⟨nHB⟩

ηgd

⟨gd⟩

r

2.9 1.7 1.8 0.9 0.8 2.3 0.4 2.6 0.9 1.2 2.0 1.7 0.7 2.3 0.4 2.5 0.9 1.2 1.9

30.7 − 3.6 − − 53.3 9.0 × 10−2 16.6 0.8 0.5 45.8 3.1 − 48.2 0.1 15.9 1.5 0.6 33.6

12 10 10 1.3 1.2 15 1.3 12.8 2.4 2.9 20.7 7.8 1.1 14.8 1.3 12.6 2.7 3.1 22.8

24 61 62 6 4 40 6 28 15 23 61 60 2 46 5 29 13 25 62

The average number of adsorbed molecules per unit cell is presented as ⟨nads/unit⟩, with applied pressures chosen such that the topological properties at saturation loading and similar absolute loadings can be compared. Other values given are the average number of H-bonds per molecule, ⟨nHB⟩; the geodesic index, ηgd; the average geodesic length, ⟨gd⟩; and the range, r, as a measure of the spread in the geodesic distributions.

a

zeolite is characterized by straight and zigzag channels with nearly spherical intersections, as can be seen in Figure 3. Each

which has σgd ̅ =

i ∑gd

i gdA(i)+ ∑gd

gdB(i) + ··· + N

i ∑gd

gdN (i) (2)

for an individual graph within the M snapshots analyzed. In eq 2, gdA(i) is a geodesic path i in which vertex A participates, and the sum is over all gd’s wherein A may be a linking vertex or a terminal vertex. While the maximum number of gd’s where A can be a terminal vertex (either a beginning or an ending of the geodesic) is N, vertex A can be a linking vertex in many more gd’s that connect other pairs of vertices. The individual sum, ∑igd gdA(i), is size extensive, as the isolated geodesic distance matrix is N×N. The quantity σ̅gd is also size extensive, as it has N vertex summations. Thus, the numerator in eq 1 is divided by both N and the total number of snapshots M that are analyzed. The scaling factor of 100 is introduced due to the sparsity of the isolated geodesic distance matrix that is being analyzed. For a complete discussion, see the Supporting Information, in which the convergence properties of this quantity are also illustrated. As a reference, consider that a perfect network of water found in hexagonal ice, where each H2O has exactly four H-bonds, has an ηgd value greater than 100, which converges to ∼160 with system size (>1500 H2O in the simulation box, see Figure S3 in Supporting Information). Given the deviation of the H-bond distribution of liquid water away from this idealized case (see Figure 2A), it is not surprising that one finds ηgd = 30.7 for the TIP4P water model having 1460 H2O in the simulation box (see Table 2). Comparisons of ηgd across the solution-phase composition provide insight into the deviations in the essential connectivity within the H-bond networks in the presence of a co-solvent. To understand the effects of confinement and hydrophobicity of silicalite-1 upon the H-bond networks of adsorbates, it is pertinent to briefly discuss silicalite-1’s structural features. The unit cell is orthorhombic, with dimensions of a = 20.1 Å, b = 19.7 Å, and c = 13.1 Å.47 The

Figure 3. (A) The (010) plane of silicalite-1 showing the straight channels, (B) the (001) plane showing the zigzag channels, and (C) a simplified structure indicating the straight and zigzag channels and their intersections.

channel is composed of a 10-membered oxygen ring; however, the diameter of each channel varies, with the zigzag channel along the (100) direction being 5.1 Å × 5.5 Å, while the straight channel along the (100) direction dimensions of 5.3 Å × 5.6 Å. The H-bond network properties of bulk water change dramatically upon adsorption into silicalite-1 due to the hydrophobicity of the interior and structural confinement. The hydrophobic nature of the interior is demonstrated by several features of the H-bond network. For example, the interaction of water with the interior Ozeo-atom of the −Si−O− Si− units was examined as a function of both distance between the water hydrogen, Hw and Ozeo, as well as the Ozeo−Hw−Ow angle at all loadings. As can be seen in Figure 4, the most likely location for adsorbed H2O is found for angles in the range of 130−140° at a distance of 3.0 Å, which has the water O-atom pointing slightly toward a neighboring silicon atom in the wall of the silicalite-1 pore. Using either a distance, or distance + angle definition of the H-bond interaction between Ozeo and the Hw of water leads to negligible amounts of H-bonding between the sorbate and sorbent. Indeed, if a distance + angle 19726

dx.doi.org/10.1021/jp502867v | J. Phys. Chem. C 2014, 118, 19723−19732

The Journal of Physical Chemistry C

Article

branching that is characteristic for bulk water. Necessarily, the average number of H-bonds per water molecule inside the zeolite is significantly diminished (1.7 H-bonds/H2O at a saturation loading of 38.5 molecules per unit cell versus 2.9 Hbonds/H2O in the bulk liquid phase). Hydrogen-Bond Networks within Bulk Methanol and Ethanol and Unary Adsorption into Silicalite-1. Recent graph-theoretical studies of bulk water in comparison to methanol have attempted to understand the H-bonding networks in water and methanol liquids in the context of small-world and or Erdős−Rényi networks known within the topological models in traditional network/graph theory.5 The topological analyses presented here complement the prior approach but focus upon a simpler, chemically intuitive analysis scheme. Bulk methanol and ethanol have distinct network characteristics relative to water, in part because of the decrease in H-atoms available for donation to H-bonds. As observed from the radial distribution functions in the bulk (see Figure S2 in the Supporting Information), both alcohols have on average two neighboring molecules with a maximum probability at ∼2.7 Å in the first solvation shell. Examination of the H-bond distribution in liquid methanol reveals that 72% of the methanol molecules have two H-bonds, while 22% have only one H-bond, with minor fractions having zero or three Hbonds (see Figure 5). This distribution is qualitatively the same for liquid ethanol, which has 70% and 24% of ethanol molecules with two H-bonds and a single H-bond, respectively. The difference in the molecular H-bond distributions relative to water does not necessarily imply that the H-bond network is less connected; however, analysis of the geodesic pathways of methanol and ethanol liquids reveals a very fragmented network with significantly more disorder than water. The average gd length in methanol is 10 (slightly less than water), but this number derives from a right-skewed distribution of Hbond pathways. As can be seen in Figure 5, there is a slight preference for H-bonded methanol dimers and trimers, but their concentration is rather low (