A Computational Investigation into the Suitability of Purely Siliceous

ARTICLE pubs.acs.org/JPCC

A Computational Investigation into the Suitability of Purely Siliceous Zeolites as Reverse Osmosis Membranes Zak E. Hughes,*,† Louise A. Carrington,‡ Paolo Raiteri,† and Julian D. Gale† † ‡

Nanochemistry Research Institute, Department of Chemistry, Curtin University, P.O. Box U1987, Perth, WA 6845, Australia School of Chemistry, University of Southampton, Southampton, SO17 1BJ, U.K. ABSTRACT: Desalination by reverse osmosis is an increasingly important source of potable water in many countries. The interest in developing new, more effective membranes is, therefore, great. One set of materials that have been suggested as a possible new type of desalination membrane are nanoporous materials. In this work computational methods are used to investigate the behavior of water within five different zeolitic systems. Quantum mechanical calculations are used to construct a set of force-field parameters for two atomistic models. Molecular dynamics simulations of the zeolites show that water will diffuse through zeolites at a rate faster than that obtained with the composite membranes currently used in commercial desalination. In addition, the thermodynamics of salt rejection have been investigated using the free energy perturbation method. The results of these calculations show that the chloride ion finds the zeolitic environment strongly unfavorable compared to the bulk solution. In the case of the sodium ion, the energetic difference between the zeolite environment and solution is less significant, but charge separation prevents sodium from permeating the membrane.

’ INTRODUCTION Currently, traditional water sources are under greater pressure than ever before.1,2 Furthermore, as the population of the world continues to grow, and the effects of climate change become increasingly felt, this pressure will only increase. Desalination of seawater or brackish water is one obvious method of meeting the increasing demand for potable water. Desalination technologies have already been in use for many years now1,2 but the high energy cost associated with such processes has meant that the technology was largely limited to areas that had access to cheap and abundant energy.1 This situation has now changed, and desalination plants are been built throughout the world, even in countries that historically receive plenty of rainfall, such as the U.K.1 Outside of the Middle East, the method of desalination used in most modern desalination plants is reverse osmosis (RO).1 In RO, water is forced through a semipermeable membrane (originally cellulose acetate membranes, although now usually a polyamide-polysulfone (PA-PS) composite1,3) at high pressure, thereby removing the impurities, especially salt, from the water. The need for high pressure means that, although more efficient than previous thermal desalination methods, RO is still an expensive process in terms of energy. Thus, there is considerable interest in developing new, more efficient membranes that will reduce the energy cost by either having greater selectivity for salt rejection or a faster flow rate than current membranes. Carbon nanotubes (CNTs) are one alternative material that has been proposed for use as an RO membrane due to the fact that the flow rates of water through CNTs are high.4,5 There has been both simulation and experimental studies of CNTs as RO r 2011 American Chemical Society

membranes,6,7 and while the simulation results are favorable, the practicalities of making a membrane consisting of CNTs with such well-defined lengths in the quantities needed for commercial use is still problematic. However, this research has established the concept that a material with uniform hydrophobic channels may offer promise as an RO membrane. Another type of material that has been suggested as a possible alternative to the current polymer-based membranes are zeolites. There has already been some work carried out with regard to using zeolites as desalination membranes8,9 or incorporating them within traditional membranes as a composite material to improve efficiency.10 Zeolites are nanoporous crystalline materials, consisting of an aluminosilicate framework of corner-sharing SiO4 or AlO4 tetrahedra with cations present in the channels in order to charge compensate for the presence of aluminum when necessary.11 Normally the incorporation of aluminum within the zeolite framework leads to such structures being hydrophilic. However, there has been considerable recent interest in preparing purely siliceous materials, either directly or by dealumination.12-14 In this form, the zeolitic silica becomes hydrophobic, and thus more closely resembles CNTs. Despite the very different channel structures, it is possible that these siliceous materials could be tailored for the purpose of desalination in the same manner, by selecting pore dimensions that are too small for diffusion of solvated sodium or chloride ions. Received: October 6, 2010 Revised: January 18, 2011 Published: February 23, 2011 4063

dx.doi.org/10.1021/jp109591f | J. Phys. Chem. C 2011, 115, 4063–4075

The Journal of Physical Chemistry C The behavior of water within siliceous zeolites has been investigated previously by both experiment15-17 and simulation,15,18-24 with most of this work concentrating on silicalite or dealuminated zeolite Y (DAY). Very little work has looked at zeolites with other pore topologies, especially those with onedimensional channel structures that may lead to more directed diffusion of water, suitable for desalination. Molecular dynamics (MD) simulations18,20 have shown that in DAY, which has larger pores, the water within the channels behaves in a manner similar to that of bulk water, whereas in silicalite a very different type of behavior is seen, with more vapor-like characteristics being observed18,20 (at 300K and 1 atm). The different pore dimensions of the two systems thus cause significant changes in the behavior and structure of the water within the zeolite. Grand Canonical Monte Carlo simulations23 also show that the water within silicalite differs significantly from bulk water, with the dipole moment being roughly 10% smaller than that of bulk water for the confined water molecules, indicating the importance of ensuring that the water-zeolite interactions are accurate. Experimental studies of siliceous zeolites have attempted to determine the self-diffusion coefficient of water within zeolites with some success. However, factors such as defects within the zeolite or a small amount of Al contamination of the framework give rise to a number of uncertainties in these results.15,18 In this study, the bulk diffusion of water through several different zeolites has been investigated using computer simulation, based on two different force-fields: one accounting for the polarizability of ions, the second a rigid ion partial charge model. All the zeolites investigated were completely siliceous, but a variety of pore shapes and dimensions have been simulated. Initial calculations were performed for silicalite, which possesses a two-dimensional channel system, being one of the best characterized zeolites, especially in siliceous form. When selecting potential candidate materials for use in desalination, it is useful to narrow the field from the several hundred known framework topologies, to a few representative examples that may have desirable characteristics. First, zeolites that have a one-dimensional channel structure most closely resemble a bundle of aligned nanotubes and therefore might be expected to performed better than higher dimensionalities, since the diffusion is focused on a single desired direction. Second, it is likely that there will be an optimal ring size, which can be characterized by the number of tetrahedra that compose the limiting dimension. If the ring size is too small, then diffusion of water will be suppressed, whereas in the converse situation a large aperture will lead to loss of selectivity. Thus four different zeolites have been selected for study that possess one-dimensional channels with 8-, 10-, and 12-rings as the limiting pore size. In the present work, the bulk diffusion characteristics of water within the zeolites are examined. Additionally, the thermodynamics of Naþ and Clions binding to the zeolites are also investigated, while the interfacial selectivity of their surfaces is deferred to future work.

’ METHODS Simulation of molecular adsorption within zeolites has advanced considerably in the past two decades, with the use of firstprinciples techniques now being routinely possible based on density functional theory.25-28 While this is a valuable approach for studying adsorption and reactivity at acidic sites within aluminosilicate zeolites, especially where questions of proton transfer arise,25 the use of quantum mechanical (QM) techniques

ARTICLE

to examine diffusivity is greatly limited by the short time scales accessible with current computers. For strongly activated diffusion, transition state theory is appropriate, and the use of such QM-based methods is possible. However, for the diffusion of water within siliceous zeolites, under conditions relevant to their application as membranes, the potential energy surface necessarily consists of many distinct minima connected by low energy barriers; thus MD is the technique of choice for configurational sampling. Furthermore, the weak nature of the interaction between water and the siliceous framework represents a difficult case for first-principles methods of the form that can be readily utilized in the solid-state. As a result, the use of a force-field description of the energy surface is not only the most efficient, and readily extended to the complex situation of the zeolite surface-aqueous interface, but is also likely to be equally as accurate as any feasible QM approach, provided it is well-parametrized. Due to their many interesting fundamental properties and industrial applications, there has been a wide range of force-fields developed for study of the behavior of water within zeolites.15,18-20,23,23,24,29-31 In choosing a force-field, a number of considerations arise. One such choice is whether the use of a rigid zeolite framework is appropriate or not. A rigid zeolite will be less computationally expensive, but potentially lead to a loss of accuracy. A number of previous works have investigated the effect of approximating the system as being rigid.32 Of particular interest in this case is the study of Demontis et al.,18 where the results of flexible and rigid silicalite are compared against each other. The overall conclusion is that the lattice vibrations of the zeolite play an important role in the diffusion of species within the zeolite, thus a flexible zeolite model is to be preferred. The large number of different water models used makes it difficult to determine the importance of flexibility in the water molecules. The accuracy of the model in its description of the water-water interactions is likely to be more important. Whether to use a force-field that accounts for the polarizability of atoms in some manner (and if so by what method) is another factor in the choice of model. Including polarizability in the forcefield allows species to adapt to their local environment to a greater degree, but again is accompanied by an increase in the computational cost of the simulations. Methods that have been used to treat polarization in zeolites include using an electric field dependent potential,18,23 the electronegativity equalization method,19 and use of a shell model.29 Rigid ion models have difficulty capturing the low-temperature monoclinic distortion of silicalite33 and are also less likely to predict the order of phase stability of silica correctly.34 We note that a recent report24 of a modified partial charge rigid-ion force-field reproducing the monoclinic distortion of silicalite was as a result of discontinuities in the energy surface. Accurate summation of the C6 terms in this model (see eq 3), through the use of an Ewald-like technique35 causes the monoclinic form to optimize to the orthorhombic structure. To explore the influence of differing force-field types, it was decided that two models would be used;one incorporating polarizability, and one neglecting it;while in both models the zeolite framework would be flexible. Comparing the results of the two models against each other will determine whether it is possible to use just a simple partial charge description (thus speeding up the calculations) while not significantly affecting the results that are relevant to this study. Shell Model. For the first model, a polarizable ion description of both the zeolite and water has been used, based on the shell 4064

dx.doi.org/10.1021/jp109591f |J. Phys. Chem. C 2011, 115, 4063–4075

The Journal of Physical Chemistry C

ARTICLE

Table 1. Interatomic Potential Parameters Used for the Shell Modela ions

core charge/e

Ion Parameters shell charge/e

Ob

0.84819

-2.84819

74.92038

OWc

1.71636

-1.71636

43.3275

Hc

0.55733

Dc

-1.11466

Sib

4.00000

ion pair Si-Ob Si-OWd

Buckingham Potential Parameters A/eV F/Å 1283.90734

0.32052

10.66158

983.55600

0.32052

10.66158

22764.3

0.14900

27.88

O-OWd

22764.3

0.14900

28.92

0.25000

0.00

396.27

ion pair

Lennard-Jones (12-6) Potential Parameters ε/eV

σ/Å

OW-OWc

0.0091532166

3.138395

ions

Three Body Potential Parameters k/eV rad-2

O-Si-Ob

ion pair

model of Dick and Overhauser.36 Here the polarizability of certain ions (specifically oxygen in the present case) is modeled by having a massive core attached to a massless shell, from which it is Coulombically screened, by a harmonic potential. The potentials used for the zeolite are those derived by Sanders et al.,29 which have been extensively used in previous simulations37-40 and have been shown to accurately describe the stability of tetrahedral phases of silica.41 De Leeuw and Parker have developed a water model designed for use in combination with the Sanders et al. silica potential.42 Unfortunately, it was later found that this water model would freeze after about 1 ns of simulation.43 As a result of these findings, the model was modified by Kerisit and Parker38 in an attempt to avoid this incorrect behavior. However, while the modified parametrization shows considerable improvement, there is evidence that water still freezes at longer time scales.44 Additionally, the density of bulk water is overestimated by 27%, and the oxygen-oxygen radial distribution function calculated from this model differs appreciably from that obtained experimentally. Due to these problems, it was decided to explore alternative force-fields for water that are known to describe the liquid phase accurately. Figure 1a shows the oxygen-oxygen radial distribution functions of three water models;the Kerisit-Parker water model, the SWM4-NDP model developed by Lamoureux et al.,45 and the TIP4P-Ew46 model;all at 300 K and 1 atm, as well as

C6/eV Å6

O-Ob O-Hd

Figure 1. Radial distribution functions of bulk water, from experiment at 298 K and 1 bar and for the SWM4-NDP,45 TIP4P-Ew46 and KerisitParker38 models at 300 K and 1 atm: (a) oxygen-oxygen and (b) oxygen-hydrogen.

kcs/eV Å-2

Si-OW O-OW

2.09274

109.47

Modified Buckingham Potential Parameters A/eV F/Å 983.55600 22764.3

θ0/o

C6/eV Å6

0.32052

0.00

0.14900

17.14

a

Here O and OW represent oxygens in silica and water, respectively, and D represents the Drude particle on the SWM4-NDP water molecule. A cut-off of 10 Å was used with a tapering function to ensure the potential decayed to zero applied from 9 Å. b Taken from ref 29. c Taken from ref 45. d Taken from refs 39 and 40.

the experimental disribution function, at 298 K and 1 bar.47 The gOO(r) of SWM4-NDP is nearly identical to that obtained from experiment,45 while the TIP4P-Ew model is almost as accurate. Ultimately SWM4-NDP was chosen, as this is a rigid four-site model, similar to TIP4P,48 but with a shell added to the oxygen. In addition to its accuracy in describing the radial distribution functions, another advantage of this model is the fact that SWM4NDP is rigid, therefore, allowing the use of a longer MD time step, provided an adiabatic algorithm is employed to handle the shell motions.49 For the zeolite-water interaction potentials, it was decided to start from the parameters developed by de Leeuw and Parker for the interaction of water with R-quartz.39,50 However, given the change of water model and the different application context, it is important to validate the choice of parameters and refit if necessary. The interaction parameters for the shell model are detailed in Table 1. The core-shell interactions of the oxygen anions are described by Ucore - shell ¼ 4065

1 kcs rcs2 2

ð1Þ



ARTICLE

Figure 2. Structure of the silicate cluster used for QM calculations of the adsorption energy of water. Dangling Si-O bonds are terminated by hydrogens. The silicon, oxygen, and hydrogen atoms are colored yellow, red, and white, respectively. The water molecule is shown in the configuration with the lowest energy.

where kcs is the force constant for the core-shell interaction, and rcs is the distance between the core and the shell. The form of the Lennard-Jones potential used in this work was " 6 # σ 12 σ ULJ ¼ 4ε ð2Þ r r where ULJ is the energy, ε is the depth of the potential well, σ is the distance at which the energy is zero (other than dissociation), and r is the distance between the two atoms. The Buckingham potential is given by -r C6 ð3Þ UBuck ¼ A exp F r where r is the distance between the two atoms, UBuck, the energy of the interaction and A, F, and C6 constants. For the three-body interactions, the energy, Uthree, is given by Uthree ¼

1 kðθ - θ0 Þ2 2

ð4Þ

where k is the force constant, θ0 is the equilibrium angle, and θ is the current angle. In the absence of direct experimental information to calibrate the heat of adsorption of water in purely siliceous zeolites, we have used QM methods to provide relevant data. While most methods are capable of describing the repulsive interaction between water and silica, the challenge is to obtain a reasonable estimate of the attractive contributions, especially those due to van der Waals forces. For this reason, we have chosen to utilize a cluster model for a siliceous framework, rather than a periodic treatment, since this readily allows the comparison of wave function-based theories instead of density functional theory alone. Given that the objective of the QM calculations is purely to estimate the water-silica interaction strength, rather than to accurately model the results for a specific material, a small 3T site cluster has been utilized, as shown in Figure 2. All terminal bonds to silicon are saturated with hydroxyl groups. The cluster alone, with the water molecule, was optimized at the PBE/cc-pvtz51-53 level within C2V symmetry. Noting that water adsorption on siliceous zeolites will only weakly perturb the framework geometry, the cluster was subsequently held fixed at the optimized configuration during the adsorption of water. Three different orientations for water adsorption were then examined, comprising the dipole moment of water being directed either toward the nearest silicon atom or away from it, and with the plane of the

water molecule parallel to the mirror plane of the cluster or orthogonal to it. The preferred configuration was identified at the PBE/cc-pvtz level and then a subsequent optimization was performed at the MP2/cc-pvtz level. A further calculation at the MP2/cc-pvqz level was performed using the MP2/cc-pvtz geometry. Binding energies were corrected for basis-set superposition error using the counterpoise correction.54 All calculations were performed using NWChem.55 The binding energies calculated using the above QM methods were then compared against those of the de Leeuw and Parker force-field model for the same system, under equivalent structural constraints, using the GULP program.56,57 Note that in all calculations the shells are allowed to relax, even if the cluster is held fixed, since these degrees of freedom are equivalent to the electronic polarization that occurs during the self-consistent field process in the QM methods. Partial Charge Rigid Ion Model. For the water model, a shellless analogue of SWM4-NDP, namely TIP4P-Ew46 was used. For the zeolite, a new model was developed from that of Sanders et al. First, the charges were reset to the average on-diagonal Born effective charge determined from the Sanders et al. model (Si = þ3.3, O = -1.65), which is known to agree with first-principles calculations values.58 The form of the potentials was kept the same as that for the model of Sanders et al., but refitted to match the mechanical properties (elastic constants), static dielectric constants, and the extremes of the phonon spectra of the Rquartz structure (Table 2 shows the parameters for the rigid ion model). As found for other rigid ion partial charge models of silica, the order of phase stability is incorrectly predicted with Rcristobalite being the most stable structure, followed by tridymite and then R-quartz. This result does not change if lower charges (such as those of BKS model59) are used. Given that we are not interested in the phase transformations of the zeolite framework, the failings with regard to the silica phase diagram are inconsequential in the present context, but should be noted for broader application. The interaction parameters between the zeolite and water employed the same values as per the shell model, but with a small adjustment in order to give a binding energy very close to those obtained from the MP2 calculation and the shell model for the silica cluster. Molecular Dynamics. Once the parameters of the force-fields had been modified to obtain better agreement with the QM results MD simulations of water diffusion within the zeolite systems were performed. Here the Velocity Verlet algorithm was used to integrate the equations of motion, as implemented in the program DL_POLY, version 2.19.60 The simulations were carried out in the NPT ensemble with the Nose-Hoover thermostat and Hoover barostat,61,62 with relaxation parameters of 1.0 and 10.0 ps, respectively. The Coulombic interactions were handled using a smooth particle mesh Ewald sum63 (with a relative accuracy of 1.0 10-6). The systems were first equilibrated during a short run of 500 ps, using a time step of 1 fs, before the 4 ns production runs were carried out, all at standard temperature and pressure. The zeolite systems chosen for study were silicalite (where previous work has been carried out using a variety of different models, thus allowing comparisons to be made) and four different zeolites with one-dimensional channel structures: MTF (8), SFF (10), VET (12), and GON (12), where the number in parentheses indicates the smallest ring size of the channel of each framework. In all cases, simulations were carried out at different 4066



ARTICLE

Table 2. Interatomic Potential Parameters for the Rigid Ion Modela Ion parameters Ions

Charge/e

O

-1.65000 b

OW H

b

method

ΔUads /kJ mol-1

PBE/cc-pvtz

-18.3

0.00000

MP2/cc-pvtz

-21.5

0.52422

MP2/cc-pvqz//MP2/cc-pvtz

-21.9

force-field shell model: initial parameters

-25.5

force-field shell model: final parameters force-field partial charge model: initial parameters

-21.2 -22.7

Force-field partial charge model: final parameters

-21.3

D

-1.04844

Si

3.30000

b

Table 3. Adsorption Energies of Water on the Hydroxyl Terminated 3T Silica Cluster as Determined from QM and Force-Field Calculations

Buckingham potential parameters Ion pair

A/eV

F/Å

C6/eV Å6

Si-O

332.54584

0.408412

0.00

Si-OW

983.56

0.32052

0.00

O-O

22764.3

0.14900

32.00

O-OW

22764.30

0.149

17.14

O-Hc

396.27

0.25

0.00

Lennard-Jones (12-6) potential parameters Ion pair

ε/eV

σ/Å

OW-OWb

0.0070520704

3.16435

Three body potential parameters Ions

k/eV rad-2

θ0/o

O-Si-O

1.10730

109.47

Modified Buckingham potential parameters Ion pair

A/eV

F/Å

C6/eV Å6

O-OW

22764.3

0.14900

15.0

a

Here O and OW represent oxygens in silica and water, respectively, and D represents the Drude particle on the TIP4P-Ew water molecule. Potentials had a tapering function applied from 4 Å to a cut-off of 12 Å. b Taken from ref 46. c Taken from refs 39 and 40.

levels of water loading. The diffusion coefficients of water within the different zeolites were calculated from the gradient of a plot of mean square displacement (MSD) versus time.

’ RESULTS AND DISCUSSION Quantum Mechanical Calculations. The binding energies of water to the zeolite system calculated from the QM results and those obtained using the initial force-field model are shown in Table 3. The lowest energy geometry was found to be that with the water molecule in the same plane as the cluster and with the hydrogens directed toward to the cluster. The binding energy calculated using MP2 theory is approximately 3 kJ mol-1 more exothermic than that calculated using the PBE method. This difference is to be expected, as MP2 more accurately describes the van der Waals interactions present in the system. Increasing the basis set from cc-pvtz to cc-pvqz only results in a small increase in binding. The original parameters used for the zeolitewater interaction in the de Leeuw-Parker force-field result in too strong an attraction. Aside from the electrostatic interaction, the dominant attractive force arises from the van der Waals contribution. As it is not desirable to alter charges of the species, the C6 terms in the Buckingham potential of the Si-OW and OOW interactions are the parameters that should be modified to

reduce the strength of the attraction of the water to the zeolite. After these modifications were made to the force-field, the binding energies are in good agreement with the results of the QM calculations. In the case of the TIP4P-Ew water and rigid ion zeolite model, the binding energy was initially calculated using the final interaction parameters for the SWM4-NDP system. This was slightly too large, and so the parameters were refitted again to give a binding energy closer to that of the SWM4-NDP model and the MP2 calculations. Silicalite. Silicalite is an all-silica zeolite having two different interconnecting channel systems, with both channels consisting of 10-membered rings. Parallel to the y-axis are elliptical straight channels approximately 5.7 5.2 Å in diameter, while parallel to the x-axis are sinusoidal channels, circular in cross-section with a diameter of ∼5.4 Å. At low temperature, silicalite possesses monoclinic symmetry (space group P21/n11), while on heating there is a reversible phase transition to orthorhombic symmetry (space group Pnma) somewhere between 300 and 350 K.65-67 As the cell parameters for the two systems are very similar and there is only a very slight change in one angle from 90.67 to 90.0, it was decided to simply use the orthorhombic cell for all simulations. This also makes comparisons with previous simulations easier as most previous work has been carried out in the orthorhombic phase. Table 4 shows that the cell parameters for the orthorhombic phase of silicalite obtained from both of the force-fields used in this study only differ by a maximum of 1.2% when compared against the values obtained by experiment. There has been some previous work determining the adsorption energy of water in siliceous zeolites, with quite a wide range of values. For silicalite, previous simulations have given values of -17,64 -20,22 -32.5,18 -34.7/-38.319 and -4323 kJ mol-1. One recent study24 identified two different types of adsorption sites in silicalite at -21.5 and -39.8 kJ mol-1. The authors of this last model used lower point charges for the zeolite framework than many other models, including the one outlined within this study, as they were concerned that a highly charged framework could result in species adsorbing too strongly to the zeolite. Using the two force-fields described in this study, the adsorption energy of water in silicalite is found to be -29.3 and -28.9 kJ mol-1 for the shell model and rigid-ion model, respectively, placing both models within the range of literature values. There have been a number of previous MD simulations of water diffusing within silicalite,15,18-20,24,68,69 usually carried out with eight water molecules per unit cell (as this corresponds to the maximum loading of the samples used for experimental measurements of the diffusion coefficient20,70,71). Table 5 summarizes the results of the present study alongside those results 4067



ARTICLE

Table 4. The Orthorhombic Cell Parameters of Silicalite Determined by Experiment84 at a Temperature of 298 K and from Both Static Calculations and MD Simulations at 300 K Using Both the Shell and Rigid Ion Modelsa shell model parameter

rigid-ion model

experiment

static

MD

static

MD 5246.07 ( 28.93

volume/Å

5365.24

5266.06

5348.58 ( 11.74

5221.17

a/Å

20.070

20.122

20.130 ( 0.015

19.826

20.000 ( 0.037

b/Å c/Å

19.920 13.420

19.813 13.209

19.823 ( 0.014 13.404 ( 0.010

19.773 13.319

19.695 ( 0.036 13.318 ( 0.024

phonon frequency/cm-1

21

34

The phonon frequency given in the table refers to the lowest phonon frequency at Γ point, excluding the three translational modes. This mode is real in both cases, indicating the phonon stability with respect to the monoclinic distortion.

a

Table 5. Comparison of Diffusion Coefficients, β and δ for Water in Silicalite at a Loading of Eight Water Molecules Per Unit Cell at 300 K, as Determined by Simulations and Experiment

Table 6. Diffusion Coefficients for Water in Silicalite as a Function of Loading at 300 K and 1 atm diffusion coefficients/10-9 m2 s-1 model

diffusion coefficients/10-9 m2 s-1 β

shell

δ

loading per unit cell

D

Dx

Dy

Dz

1

4.62 ( 1.51

7.61

5.74

0.86

8

1.65 ( 0.19

2.16

2.36

0.43 0.29

D

Dx

Dy

Dz

Bussai et al. 15

3.3

2.6

6.5

0.19

1.06 5.76

16

1.00 ( 0.09

1.20

1.47

Demontis et al. 18

8.6

7.9

15.9

2.0

1.19 5.95

32

0.98 ( 0.06

1.21

1.31

0.31

Smirnov et al. 20

2.3/4.1

1

3.47 ( 1.02

3.13

6.27

1.00

study/model

Fleys et al. 69 Ari et al. 24

8.83 1.94

Bordat et al.24

rigid-ion 2.33 0.69

6.17 1.24

0.67 0.00

1.14 6.34

8

1.23 ( 0.15

1.35

1.89

0.29

3.5

6.35

0.9

1.11 5.47

16 32

1.01 ( 0.03 0.91 ( 0.04

1.35 1.23

1.58 1.15

0.21 0.28

PG-NMR 15

1.7

SWM4-NDP

1.65

2.16

2.36

0.43

1.18 5.24

TIP4P-Ew

1.18

1.35

1.89

0.29

1.23 5.62

of previous MD simulations and experimental data, determined from pulsed field gradient (PFG) NMR measurements.15 In addition to the overall three-dimensional diffusion coefficient, the components along each axis are also given. Also calculated are two parameters that are often used to compare the diffusion of water in microporous materials,72 β and δ. The first of these parameters, β, measures the discrepancy from a simple random walk model: β¼

c2 =Dz a2 =Dx þ b2 =Dy

ð5Þ

where a, b, and c are the dimensions of the silicalite unit cell and Dx, Dy, and Dz are the components of the diffusion coefficient along the x, y, and z axes, respectively. The random walk model applies when β = 1. If β < 1, then the water molecules prefer to move between the different channel types, whereas if β > 1 (which is the case for water in silicalite under these conditions), then the water molecules prefer to diffuse along the same channel type after passing an intersection. The diffusion anisotropy parameter, δ, is given by δ¼

1=2ðDx þ DyÞ Dz

ð6Þ

and in silicalite should be ∼4.4 for a random walk model.72 Diffusion parameters can vary widely with small changes in force-field; indeed it is worth noting that even when the same force-field is used, different simulation parameters can result in

quite a large change in D, as the results of Ari et al.69 and Fleys et al.20 show for this very system. It is, therefore, not surprising that the different simulations of silicalite (all at 300 K and eight water molecules per unit cell) show a spread of values for the diffusion coefficient (1.94-8.83 10-9 m2 s-1). In this study, the diffusion rates obtained for silicalite at the same loading are 1.65 10-9 m2 s-1 for SWM4-NDP and 1.18 10-9 m2 s-1 for TIP4P-Ew. These values are a little lower than those determined from previous simulations, but are within the same order of magnitude, and show good agreement, especially in the case of the shell model, with the experimental value. There is also good agreement between the two force-fields and previous simulations for the values of β and δ. The effect of water loading on diffusion has also been investigated, with simulations performed with 4, 8, 16, and 32 water molecules per unit cell, with the results being shown in Table 6. As one would expect, for both models the rate of diffusion decreases as the number of water molecules increases. In addition, the higher the loading, the better the agreement between the two models. This is partly due to the fact that, with more water molecules in the system the sampling will be improved, but also indicates that at high loadings the diffusive behavior is dominated by the water-water interactions. The radial distribution functions of the two force-fields at the different loadings are shown in Figures 3 and 4. While there are some differences in gOO(r) and gOH(r) for the two force-fields, the overall pattern is broadly the same: gOO(r) displays a sharp first peak at about 2.8 Å (TIP4P-Ew) or 2.9 Å (SWM4-NDP) but the second peak, characteristic of the tetrahedral arrangement of water molecules in the bulk liquid, is broad and ill-defined, tending to become a shoulder to the first peak. In the case of 4068



Figure 3. Oxygen-oxygen radial distribution function of (a) SWM4NDP water and (b) TIP4P-Ew water within silicalite as a function of loading.

gOH(r), both force-fields show two peaks (at about 1.9 and 3.2 Å) with an extended shoulder on the second peak. The radial distribution functions indicate that while hydrogen bonds are present, the structures of the water clusters in the zeolites are less well ordered than in bulk water, in line with the results of previous simulations.18,20 As the amount of water within the zeolite is increased, the number of hydrogen bonds formed (given by the integral of the first peak in the gOH(r)) increases, as shown in Table 7. This suggests that at low loadings the water molecules do not have a liquid structure but are more vapor-like, and as the water loading is increased the structure of the water becomes more like that of the bulk liquid. Figure 5 shows how the structure of the water within silicalite changes as the loading increases. At lower loadings the tendency for water to form hydrogen bonded chains is clearly evident. The above results for the well-characterized case of silicalite confirm that the modified water-zeolite potentials that have been examined in this paper are suitable to use for modeling the behavior of water within siliceous zeolites. Zeolites with One-Dimensional Channel Systems. The four zeolites with one-dimensional channels were chosen such that a variety of channel sizes and shapes could be investigated (Figure 6). MTF, VET, and SFF all have nearly circular channels of different sizes, while GON has an elliptical channel structure. The channels of MTF are smaller than those of silicalite, consisting of an eight-ring channel (diameter 3.75 Å); those of SFF are

ARTICLE

Figure 4. Oxygen-hydrogen radial distribution function of (a) SWM4NDP water and (b) TIP4P-Ew water within silicalite as a function of loading.

Table 7. Number of Hydrogen Bonds Per Water Molecule in Silicalite Given as a Ratio to the Bulk Liquid as a Function of Loadinga model shell

rigid-ion

a

loading per unit cell

ratio of hydrogen bonds to bulk water

8

0.470

16

0.566

32

0.578

8

0.599

16

0.626

32

0.655

In bulk liquid a water molcule forms approximately 3.5 hydrogen bonds.

almost the same size (10-membered ring, diameter 5.5 Å) while in VET the channels are a little larger (12-ring structure, 5.9 Å diameter). GON also has channels larger than silicalite, consisting of 12-membered rings with dimensions of approximately 5.9 5.0 Å. Unlike silicalite, there has been little previous simulation or experimental work that has studied these zeolites. This makes it difficult to know how much water would enter the zeolite, particularly under the conditions needed for RO. Thus simulations with a number of different water loadings have been performed. The levels of loading for the one-dimensional zeolites 4069



ARTICLE

Table 8. Diffusion Coefficients of Water within the Four Zeolites with One-Dimensional Channel Structures at Different Water Loadings (all at 300 K and 1 atm) water loading zeolite MTF

SFF

Figure 5. Snapshots from the simulations of TIP4P-Ew water within silicalite at (a) eight waters per unit cell and (b) 32 waters per unit cell. Silicon, oxygen and hydrogen are colored yellow, red, and white, respectively.

VET

GON

Figure 6. Structures of four zeolites considered in the present work, as viewed along the along the one-dimensional channel: (a) MTF, (b) SFF, (c) VET, (d) GON. Silicon and oxygen atoms are colored yellow and red, respectively.

span a similar range to the loadings in silicalite, from a single water per channel (∼3 mg/g), to a loading level equivalent to eight waters per unit cell of silicalite (∼25 mg/g) and high loadings of around 100 mg/g. The diffusion coefficients for different water loadings are presented in Table 8. The quantitative agreement between the two water models is not as good for these four zeolites as it is for silicalite. However, the two models do agree with each other qualitatively in a number of aspects. For both models, the rate of diffusion decreases as the number of water molecules within the channels is increased. Additionally, at the highest loadings the diffusion coefficients of SFF, VET, and GON are very similar, while the rate of diffusion in MTF is somewhat lower. This indicates that once a channel has reached a diameter of ∼5.5 Å, the size of the channel does not strongly influence the rate of diffusion. In the case of the two zeolites with smaller channels, MTF and SFF, the diffusion coefficients of the two different water models agree well once the loading per channel is greater than a single water molecule. This is not true for VET and GON, where (except at the highest loading) the diffusion coefficients obtained using the TIP4P-Ew model are larger than those obtained using the shell model. In addition, despite the fact that the channels of the two zeolites are approximately the same size, the shell model gives a

per channel

diffusion coefficients/10-9 m2 s-1 SWM-NDP

TIP4P-Ew

1

2.26 ( 0.92

5.96 ( 1.43

3

0.62 ( 0.11

0.44 ( 0.14

5

0.37 ( 0.07

0.40 ( 0.08

17

0.24 ( 0.02

0.15 ( 0.02

1 3

1.91 ( 0.61 1.00 ( 0.40

4.21 ( 0.64 1.96 ( 0.56

5

0.73 ( 0.16

0.49 ( 0.09

10

0.53 ( 0.12

0.45 ( 0.03

32

0.49 ( 0.02

0.38 ( 0.02

1

2.39 ( 1.39

22.09 ( 5.11

3

2.06 ( 0.47

9.43 ( 2.98

6

1.49 ( 0.29

5.23 ( 1.3

32 1

0.41 ( 0.01 8.01 ( 3.68

0.34 ( 0.02 25.43 ( 6.26

3

7.16 ( 0.89

9.57 ( 3.15

6

2.36 ( 0.55

5.11 ( 1.10

32

0.51 ( 0.03

0.32 ( 0.02

diffusion coefficient for GON, which is considerably larger than that of VET, whereas the rates of diffusion calculated using the rigid-ion model are very similar. To investigate the reason for the difference between VET and GON, the energy profile of a single water molecule passing through the channels of both structures was calculated. The procedure used for these calculations was to place a single water molecule at a series of points along the z-axis, defined to be the axis parallel to the channel (the water molecules were initially placed in the center of the pore with respect to the x and y dimensions). The z-coordinate of the water oxygen was held fixed while the system was optimized. The results of these calculations (Figure 7) show that for the shell model the water molecule adsorbs more strongly to the zeolite VET than GON. For the rigid-ion model, the difference between the two zeolites is minor, and the adsorption energies show far less variation as a function of the z-axis position than for the shell model. It can be concluded, therefore, that the rigid-ion model does not manage to capture the zeolite-water interaction as well as a model incorporating polarizability, and that at low water concentrations this can lead to differences in behavior of water within the zeolite. However, as the number of water molecules present in each channel increases, the differences in the rates of diffusion for the two zeolites becomes less marked. This is most likely due to the fact that, at the higher concentrations, the water-water interactions become the dominant force in determining behavior. Consequently, the water molecules are unable to get “trapped” at points along the zeolite due to the forces exerted on them by other water molecules. In regard to their potential for desalination, all the zeolites exhibit diffusion rates that compare favorably with those found in simulations of the traditional polyamide membranes, which have a diffusion rate of ∼0.2 10-9 m2 s-1.73,74 Therefore, in terms of flow rate through the membrane, siliceous zeolites are suitable for use as RO membranes. 4070



ARTICLE

Figure 7. Adsorption energy of a water molecule within the zeolites VET and GON, as a function of fractional displacement along the z-axis of the unit cell: (a) SWM4-NDP and (b) TIP4P-Ew.

The structure of water within the channels of these zeolites is similar to the structure of water within silicalite; at low loadings the water behaves more like a vapor, but as the water concentration is increased, the number of hydrogen bonds the water molecules form also increases, indicating that the structure of the water is becoming more like that of the bulk liquid. Furthermore, as for silicalite, the rigid-ion system with the TIP4P-Ew water has a structure closer to bulk water than the SWM4-NDP water. Interaction of the Zeolites with Salt Ions. Flow rate is only one factor that will determine the suitability of a material as a RO membrane; the salt rejection rate is another important criterion that needs evaluating. While it is not possible to determine the salt rejection rate of the zeolites when only simulating the bulk zeolites, it is possible to investigate how favorable the ions find the zeolite channels compared with bulk solution. This can be computed using the free energy perturbation method (FEP),75 in which the free energy difference between two states can be calculated from ΔG ¼ -

1 ln hexp½ - βΔUi0 β

ð7Þ

where ΔU is the difference in potential energy between the two states and Æ...æ0 is the ensemble average of configurations sampled in state 0. FEP calculations are often assisted by using a process called staging, where a series of intermediate states are constructed

between the two states of interest, and the total free energy difference is the sum of all the free energy differences between neighboring states. An example where FEP has been used successfully in this manner is in the calculation of the hydration/ dehydration free energy of water76,77 or ions.44,78 In these calculations, the free energy of hydration (dehydration) is calculated by creating (destroying) the solute of interest inside the solvent by gradually increasing (decreasing) the solute-solvent interactions, from nothing to their full values. The favorability of the zeolite environment for the ions can, therefore, be determined using FEP methods. In the desalination process it is likely that if an ion does enter the zeolite it will be hydrated. Initially, however, the FEP method was used on ions in dry zeolites, and the results checked against results obtained from Mott-Littleton79 (ML) calculations in order to ensure the validity of the FEP technique for these systems. After validation, further simulations with the hydrated ions were performed. All FEP calculations described below were performed using the rigid-ion model, and the potential parameters for the Naþ and Cl- ions were taken from the work of Joung and Cheatham,78 where the ion parameters were specifically fitted for use with TIP4P-Ew water. The parameters for the interactions between the ions and the zeolite were also taken from previous work.80,81 In the case of the simulations carried out in this paper, the same procedure as described by Raiteri et al. for the calculation of the hydration free energy of Ca2þ was used,44 except where noted. The simulations used for the FEP calculations involve a single ion placed in the zeolite and consist of 10 ps of equilibration followed by a 1 ns production period. In addition to calculating the free energy needed to remove the ions from all four zeolite systems, the free energy released for the reverse process was also calculated for a number of the systems as a check. The differences between the results was within the statistical uncertainty of the free energies. The results of the FEP simulations are also compared against results of binding energies calculated using the static ML method,79 determined using an inner region of 12 Å and outer region of 60 Å. The GULP program was used for the ML calculations, while the FEP simulations were performed using a modified version of DL_POLY 2.19. The entropic and enthalpic contributions of the free energies were determined using the method described by Smith et al.82 The results of both the FEP simulations and the ML calculations are shown in Table 9. For comparison, the free energies of dehydration for Naþ and Cl- are 371.1 and 374.9 kJ mol-1.78 The agreement between the two methods is good, with the enthalpies determined from the ML calculations being within the statistical uncertainty of the enthalpies obtained from the FEP results. The difference that there is between the two sets of results is to be expected, as the ML calculations are performed at 0 K while the FEP simulations took place at 300 K and thus include a contribution from the heat capacity. As might be expected for a single ion binding to a crystalline material, the entropic contribution to the free energy is very small. Both methods show that the sodium ion finds the zeolite environment more favorable than the chloride ion, but for both ions the channels of all four zeolites are unfavorable when compared with bulk solution. The minimum energy configuration for both ions in all four of the zeolites is a configuration where the ion is close to the wall of the channel (as shown in Figure 8). Here, the sodium ion bridges two oxygens in the zeolite framework, at a distance of 2.8-3.0 Å. The larger Cl- cannot get as close to the zeolite framework, but similarly adopts a position where it is approximately equidistant from two silicon atoms, at 3.6-3.8 Å. 4071



ARTICLE

Table 9. The Free Energy Needed to Remove a Naþ or Cl- Ion from a Zeolite and the Associated Enthalpy Calculated Using the FEP Technique, Compared with the Enthalpy of the Same Process Calculated Using the ML Methoda removal of ion from zeolite ion Na

Cl

a

zeolite

-1

ΔG (FEP) /kJ mol

ΔH (FEP)/ kJ mol-1

ΔH (ML)/ kJ mol-1

-ΔGaq-zeo/ kJ mol-1

MTF

320.0 ( 1.6

322.8 ( 8.0

323.6

51.1

SFF

293.0 ( 1.1

293.0 ( 8.6

291.4

78.1

VET GON

324.7 ( 1.2 318.4 ( 1.9

326.1 ( 9.9 325.5 ( 10.4

325.1 320.8

46.4 52.7

MTF

16.2 ( 1.3

15.9 ( 4.6

16.3

358.7

SFF

24.2 ( 1.9

26.6 ( 5.6

31.5

350.7

VET

3.4 ( 1.6

3.4 ( 4.3

1.8

371.5

GON

12.7 ( 1.4

14.1 ( 4.9

12.9

362.2

The overall free energy needed to remove the ion from bulk solution and incorporate it within the zeolite, -ΔGaq-zeo, is also given.

Figure 8. Examples of the minimum energy configurations of the unsolvated ions in the MTF zeolite frameworks: (a) Naþ viewed down the z-axis (b) Cl- viewed down the z-axis, (c) Naþ viewed down the yaxis, and (d) Cl- viewed down the y-axis. Naþ, Cl-, oxygen, and silicon are colored blue, cyan, red, and yellow, respectively. The closest contacts of the zeolite framework to the ions are shown as spheres.

After the above validation of the method, the free energies needed to remove an ion from a zeolite containing water were calculated;a more relevant process for an RO membrane than the previous anhydrous case. For these simulations, the ion and its first solvation shell (six and seven water molecules for Naþ and Cl-, respectively) were placed in a channel of each of the four zeolites, and the free energy associated with the removal of the ion calculated as described above. Snapshots of the hydrated ions in the VET zeolite are shown in Figure 9, the dimensions of the channel cause some disruption of the solvation shell, especially in the case of the larger chloride anion, but it is clear that the ions preferentially interact with the water molecules rather than the zeolite framework. The free energy needed to remove the solvation shell, one water molecule at a time, was also calculated for all four zeolites. The uncorrected dehydration free energy of a TIP4P-Ew water molecule at 298 K is 29.0 kJ mol-1.76 There are a number of corrections to this value proposed by Horn et al., principally arising from the fact that TIP4P-Ew is a rigid molecule and, as such, has the same structure and properties in both the gas and aqueous phases. For this paper, any corrections to the free

Figure 9. Snapshot from the FEP simulations of the solvated ions in the VET zeolite frameworks: (a) Cl- viewed down the z-axis (b) Naþ viewed down the z-axis, (c) Cl- viewed down the y-axis, and (d) Naþ viewed down the y-axis. Naþ, Cl-, oxygen, hydrogen, and silicon are colored blue, cyan, red, white, and yellow, respectively.

energy calculations for TIP4P-Ew water have been neglected, as we are concerned with the water migrating between the liquid phase and the solid state, rather than a molecule in the gas phase. The free energies to remove the ions from the solvated zeolites are given in Table 10, while the total free energy needed to remove six TIP4P-Ew water molecules from the zeolites is 46.0, 53.6, 50.5, and 46.5 kJ mol-1 for MTF, SFF, VET, and GON, respectively. The incremental free energies are also given in Table 11. As in the case of the dehydrated zeolite, the free energies are far more favorable for Naþ than for Cl-, with the result that the chloride ion is highly unlikely to enter the zeolite. In the case of the sodium ion, the total free energy difference for the process described below þ Naþ ðaqÞ þ 6H2 OðaqÞ þ Zeolite f Zeolite 3 Na 3 6H2 O

ð8Þ

is -3.1, 12.2, -9.1, and 0.6 kJ mol-1 for MTF, SFF, VET, and GON, respectively. This result would seem to indicate that the rejection rate of sodium ions for a zeolitic desalination membrane 4072



ARTICLE

Table 10. The Free Energy Needed to Remove a Naþ or ClIon to the Gas Phase from Four Different Hydrated Zeolites ion

ΔG/kJ mol-1

zeolite

Na

Cl

MTF

502.2 ( 3.7

SFF VET

479.3 ( 2.2 503.7 ( 5.1

GON

498.0 ( 2.9

MTF

18.5 ( 1.7

SFF

25.2 ( 2.1

VET

2.8 ( 1.6

GON

13.3 ( 0.8

Table 11. The Free Energy Needed to Remove Six TIP4P-Ew Water Molecules from the Zeolites SFF and GON ΔG/kJ mol-1 water molecule

MTF

SFF

VET

GON

first

11.1 ( 2.4

13.2 ( 1.2

11.6 ( 1.7

10.9 ( 2.6

second

10.3 ( 1.1

12.2 ( 0.9

11.1 ( 1.5

10.8 ( 1.0

third

8.4 ( 1.3

10.9 ( 1.1

9.9 ( 0.9

10.1 ( 0.6

fourth fifth

8.1 ( 0.8 7.0 ( 0.9

9.4 ( 1.1 6.5 ( 0.9

9.6 ( 1.0 7.2 ( 1.1

7.9 ( 0.8 6.3 ( 1.4

sixth

1.0 ( 0.9

1.4 ( 0.9

1.4 ( 1.0

0.5 ( 1.0

Total

46.0

53.6

50.5

46.5

might be rather poor. However, there are two important factors that need to be considered. First, unless a counterion (such as Cl- or OH-) enters the zeolite simultaneously there will be a extra free energy penalty due to charge separation. As the dielectric constant of siliceous zeolites is low (for the four onedimensional zeolites in this study the dielectric constant was calculated at 4-5) when compared to water (78.4), this penalty will be significant. This means that, while Naþ ions might enter the zeolite channels, they are unlikely to diffuse more than a few angstroms from the membrane interface. Second, the above calculations cannot make any prediction regarding the free energy barrier to pass through the entrance of the zeolite channel. Therefore, while it is possible that some sodium ions will diffuse into the zeolite membranes, the overall salt rejection rate is likely to be higher than the FEP calculations indicate. When an ion and water molecules are both present in a zeolite channel, the solvation shell of the ion, whether Naþ or Cl-, is largely maintained, with the ion near the center of the zeolite channel with the water molecules surrounding it. This indicates that the ion-water interaction is the dominant interaction.

’ CONCLUSIONS The applicability of zeolites as desalination membranes has been investigated using computational methods. Two different types of atomistic force-fields, a shell model incorporating polarizability and a nonpolarizable rigid-ion model, were used in the course of this investigation. The parameter sets of both forcefields were refined based on the results of MP2 calculations of a water molecule binding to a hydroxyl-terminated silicate cluster. Using the modified force-fields MD simulations of water in five different zeolite structures were carried out, with the aim of determining their appropriateness for use as RO membranes. Pore shapes, pore size, the dimensionality of the zeolites channel

structure, and the amount of water within the pore were all varied. The zeolite silicalite was used as a test system due to fact that there have been a number of previous studies investigating the behavior of water within silicalite. As in previous studies, both new models exhibit a structure of water within the zeolite that is significantly different from that of bulk water. The diffusion coefficients obtained using the force-fields described in this paper are similar to those obtained from previous simulations and experimental work. This indicates that the potentials developed provide a reasonable description of the water-zeolite interaction. The other four zeolites investigated all possess one-dimensional channel structures, but with a variety of different pore shapes and sizes. The agreement between the two force-fields was less consistent than in the case of silicalite, particularly when the loading of water within the zeolite was low. The differences between the two models were particularly evident for the GON and VET zeolites. For the rigid ion model, there was little difference in the strength of adsorption of a water molecule to the zeolites. With the shell model, however, the binding of water to VET was significantly greater than for GON. The fact that the shell model incorporates polarizability allows it to better capture the interaction of water molecules with the zeolite framework at low loadings. As the water concentration increases, the differences between the two models are relatively minor, making the computationally cheaper rigid-ion model preferable. Regardless of the model used, the diffusion coefficients of water, even at the highest loadings, are greater for all but one zeolite (MTF) than those obtained using the typical membranes currently favored in RO desalination. The rate of diffusion of water molecules, in short, defect-free CNTs of a similar pore radius (an internal diameter of 4.7 Å), was calculated as being half that of a water molecule in the aqueous phase.83 In comparison, at the highest loadings, the diffusion rate of water in the one-dimensional zeolites simulated in this report was approximately a sixth for the shell model, or an eighth for the rigid-ion model, of the diffusion rate of bulk water. Therefore, although water diffusion in siliceous zeolites is slower than in CNTs, they would be competitive for practical application In addition to RO membranes requiring a fast flow rate, they also need a high salt selectivity. In order to gain information on this aspect, the free energy of the interaction of salt ions with the zeolites was investigated for the rigid-ion model. The results for the chloride ion are promising; the insertion of Cl-, either as single ion or with a solvation shell, into the zeolite was strongly unfavorable compared with the ion remaining in solution. In the case of the sodium ion, the ion by itself is unlikely to enter the zeolite without its solvation shell. When the first solvation shell of the Naþ ion is also taken into account, the ion does not show a significant preference for the bulk solution over a zeolitic environment. However, the need for counterions to balance the charge means that the salt rejection rate will be higher than simple consideration of the free energies suggests. Even so, a build up of sodium ions near the membrane interface may be an issue, as, for a one-dimensional zeolite, channels could be obstructed and the flow rate reduced. This blockage of pores could also occur due to foulant molecules present in the seawater. Moreover, as the one-dimensional zeolites do not show a higher rate of diffusion in comparison with silicalite, zeolites with two-dimensional channel systems may be equally as, or more, effective as RO membranes, since they will be less susceptible to flow rate reduction due to fouling. 4073


The Journal of Physical Chemistry C To further investigate the effectiveness of siliceous zeolites for RO, the interactions of Naþ, Cl-, and water at the membrane surface, and especially around the channel openings also needs to be investigated. Overall, the results outlined in this study indicate that zeolites show promise as a new type of desalination membrane, and that simulation can provide an appropriate basis for further analysis to be carried out.

’ ACKNOWLEDGMENT This research was supported by the Commonwealth Scientific and Industrial Research Organisation (CSIRO) under the “Water For A Healthy Country” Flagship. J.D.G. and P.R. thank the Australian Research Council (ARC) for fellowships. The authors thank iVEC and the NCI for the provision of computational resources. ’ REFERENCES (1) Fritzmann, C.; Loewenberg, J.; Wintgens, T.; Melin, T. Desalination 2007, 216, 1–76. (2) Charcosset, C. Desalination 2009, 245, 214–231. (3) Petersen, R. J. J. Membr. Sci. 1993, 83, 81–150. (4) Brovchenko, I.; Geiger, A.; Oleinikova, A.; Paschek, D. Eur. Phys. J. E 2003, 12, 69–76. (5) Striolo, A. Nano Lett. 2006, 6, 633–639. (6) Corry, B. J. Phys. Chem. B 2008, 112, 1427–1434. (7) Kalra, A.; Garde, S.; Hummer, G. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 10175–10180. (8) Li, L. X.; Dong, J. H.; Nenoff, T. M.; Lee, R. J. Membr. Sci. 2004, 243, 401–404. (9) Duke, M. C.; O’Brien-Abraham, J.; Milne, N.; Zhu, B.; Lin, J. Y. S.; da Costa, J. C. D. Sep. Purif. Technol. 2009, 68, 343–350. (10) Jeong, B.; Hoek, E.; Yan, Y.; Subramani, A.; Huang, X.; Hurwitz, G.; Ghosh, A.; Jawor, A. J. Membr. Sci. 2007, 294, 1–7. (11) Breck, D. W. Zeolite Molecular Sieves: Structure, Chemistry and Use; Wiley: New York, 1974. (12) Muller, M.; Harvey, G.; Prins, R. Microporous Mesoporous Mater. 2000, 34, 135–147. (13) Morales-Pacheco, P.; Alvarez-Ramirez, F.; del Angel, P.; Bucio, L.; Dominguez, J. M. J. Phys. Chem. C 2007, 111, 2368–2378. (14) Apelian, M. R.; Fung, A. S.; Kennedy, G. J.; Degnan, T. F. J. Phys. Chem. 1996, 100, 16577–16583. (15) Bussai, C.; Vasenkov, S.; Liu, H.; Bohlmann, W.; Fritzsche, S.; Hannongbua, S.; Haberlandt, R.; Karger, J. Appl. Catal. A: Gen. 2002, 232, 59–66. (16) Caro, J.; Hocevar, S.; Karger, J.; Riekert, L. Zeolites 1986, 6, 213–216. (17) Caro, J.; Bulow, M.; Richter-Mendau, J.; Karger, J.; Hunger, M.; Freude, D.; Rees, L. V. C. J. Chem. Soc., Faraday Trans. 1 1987, 83, 1843. (18) Demontis, P.; Stara, G.; Suffritti, G. B. J. Phys. Chem. B 2003, 107, 4426–4436. (19) Smirnov, K.; Bougeard, D. Chem. Phys. 2003, 292, 53–70. (20) Fleys, M.; Thompson, R. W.; MacDonald, J. C. J. Phys. Chem. B 2004, 108, 12197–12203. (21) Desbiens, N.; Boutin, A.; Demachy, I. J. Phys. Chem. B 2005, 109, 24071–24076. (22) Fleys, M.; Thompson, R. J. Chem. Theory. Comput. 2005, 1, 453–458. (23) Puibasset, J.; Pellenq, R. J. M. J. Phys. Chem. B 2008, 112, 6390–6397. (24) Bordat, P.; Cazade, P.-A.; Baraille, I.; Brown, R. J. Chem. Phys. 2010, 132, 094501. (25) Shah, R.; Payne, M. C.; Lee, M. H.; Gale, J. D. Science 1996, 271, 1395–1397. (26) Haase, F.; Sauer, J. Microporous Mesoporous Mater. 2000, 3536, 379–385.

ARTICLE

(27) Benco, L.; Demuth, T.; Hafner, J.; Hutschka, F.; Toulhoat, H. J. Chem. Phys. 2001, 114, 6327–6334. (28) Hytha, M.; Stich, I.; Gale, J. D.; Terakura, K.; Payne, M. C. Chem.—Eur. J. 2001, 7, 2521–2527. (29) Sanders, M. J.; Leslie, M.; Catlow, C. R. A. J. Chem. Soc., Chem. Commun. 1984, 1271–1273. (30) Faux, D.; Smith, W.; Forester, T. J. Phys. Chem. B 1997, 101, 1762–1768. (31) Ockwig, N. W.; Cygan, R. T.; Criscenti, L. J.; Nenoff, T. M. Phys. Chem. Chem. Phys. 2008, 10, 800–807. (32) Demontis, P.; Suffritti, G. B. Microporous Mesoporous Mater. 2009, 125, 160–168. (33) Bell, R. G.; Jackson, R. A.; Catlow, C. R. A. J. Chem. Soc., Chem. Commun. 1990, 10, 782–783. (34) Martonak, R.; Donadio, D.; Oganov, A. R.; Parrinello, M. Phys. Rev. B 2007, 76, 014120. (35) Williams, D. E. Cryst. Rev. 1989, 2, 3–23. (36) Dick, B. G.; Overhauser, A. W. Phys. Rev. 1958, 112, 90–103. (37) Du, Z.; de Leeuw, N. H. Surf. Sci. 2004, 554, 193–210. (38) Kerisit, S.; Parker, S. C.; Harding, J. H. J. Phys. Chem. B 2003, 107, 7676–7682. (39) Du, Z.; de Leeuw, N. H. Dalton Trans. 2006, 2623–2634. (40) de Leeuw, N. H.; Higgins, F. M.; Parker, S. C. J. Phys. Chem. B 1999, 103, 1270–1277. (41) Henson, N. J.; Cheetham, A. K.; Gale, J. D. Chem. Mater. 1994, 6, 1647–1650. (42) de Leeuw, N. H.; Parker, S. C. Phys. Rev. B 1998, 58, 13901–13908. (43) van Maaren, P. J.; van der Spoel, D. J. Phys. Chem. B 2001, 105, 2618–2626. (44) Raiteri, P.; Gale, J. D.; Quigley, D.; Rodger, P. M. J. Phys. Chem. C 2010, 114, 5997–6010. (45) Lamoureux, G.; Harder, E.; Vorobyov, I. V.; Roux, B.; MacKerell, A. D. Chem. Phys. Lett. 2006, 418, 245–249. (46) Horn, H.; Swope, W.; Pitera, J.; Madura, J.; Dick, T.; Hura, G.; Head-Gordon, T. J. Chem. Phys. 2004, 120, 9665–9678. (47) Soper, A. K. Chem. Phys. 2000, 258, 121–137. (48) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (49) Lindan, P. J. D.; Gillan, M. J. J. Phys.: Condens. Matter 1993, 5, 1019–1030. (50) de Leeuw, N. H.; Parker, S. C. Mol. Simul. 2000, 24, 71–86. (51) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865–3868. (52) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007–1023. (53) Woon, D. E.; Dunning, T. H. J. Chem. Phys. 1993, 99, 3730–3737. (54) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553–566. (55) Kendall, R. A.; Apra, E.; Bernholdt, D. E.; Bylaska, E. J.; Dupuis, M.; Fann, G. I.; Harrison, R. J.; Ju, J. L.; Nichols, J. A.; Nieplocha, J.; Straatsma, T. P.; Windus, T. L.; Wong, A. T. Comput. Phys. Commun. 2000, 128, 260–283. (56) Gale, J. D. J. Chem. Soc., Faraday Trans. 1997, 93, 629–637. (57) Gale, J. D.; Rohl, A. L. Mol. Simul. 2003, 29, 291–341. (58) Gonze, X.; Allan, D. C.; Teter, M. P. Phys. Rev. Lett. 1992, 68, 3603–3606. (59) van Beest, B. W. H.; Kramer, G. J.; van Santen, R. A. Phys. Rev. Lett. 1990, 64, 1955–1958. (60) Smith, W.; Forester, T. R. J. Mol. Graphics 1996, 14, 136–141. (61) Nose, S. J. Chem. Phys. 1984, 81, 511–519. (62) Hoover, W. G. Phys. Rev. A 1985, 31, 1695–1697. (63) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577–8593. (64) Cailliez, F.; Stirnemann, G.; Boutin, A.; Demachy, I.; Fuchs, A. H. J. Phys. Chem. C 2008, 112, 10435–10445. (65) Hay, D. G.; Jaeger, H. J. Chem. Soc., Chem. Commun. 1984, 1433–1433. (66) Fyfe, C. A.; Kennedy, G. J.; Kokotailo, G.; Lyerla, J. R.; Fleming, W. W. J. Chem. Soc., Chem. Commun. 1985, 740–742. 4074



ARTICLE

(67) Grau-Crespo, R.; Acuay, E.; Ruiz-Salvador, R. R. Chem. Commun. 2002, 2544–2545. (68) Bussai, C.; Hannongbua, S.; Fritzsche, S.; Haberlandt, R. Chem. Phys. Lett. 2002, 354, 310–315. (69) Ari, M. U.; Ahunbay, M. G.; Yurtsever, M.; Erdem-Senatalar, A. J. Phys. Chem. B 2009, 113, 8073–8079. (70) Giaya, A.; Thompson, R. W.; Denkewicz, R. Microporous Mesoporous Mater. 2000, 40, 205–218. (71) Giaya, A.; Thompson, R. W. Microporous Mesoporous Mater. 2002, 55, 265–274. (72) Karger, J.; Ruthven, D. M. Diffusion in Zeolites and Other Microporous Solids; Wiley: New York, 1992. (73) Kotelyanskii, M. J.; Wagner, N. J.; Paulaitis, M. E. J. Membr. Sci. 1998, 139, 1–16. (74) Hughes, Z. E.; Gale, J. D. J. Mater. Chem. 2010, 20, 7788–7799. (75) Chipot, C.; Pohorille, A. Free Energy Calculations; Springer: New York, 2007. (76) Horn, H.; Swope, W.; Pitera, J. J. Chem. Phys. 2005, 123, 194504. (77) Karamertzanis, P. G.; Raiteri, P.; Galindo, A. J. Chem. Theory Comput. 2010, 6, 1590–1607. (78) Joung, I. S.; Cheatham, T. E. J. Phys. Chem. B 2008, 112, 9020–9041. (79) Mott, M. F.; Littleton, M. J. Trans. Faraday Soc. 1938, 34, 485–499. (80) Binks, D. J.; Ph.D. Thesis, University of Surrey, U.K., 1994. (81) Cygan, R.; Liang, J.; Kalinichev, A. J. Phys. Chem. B 2004, 108, 1255–1266. (82) Smith, D. E.; Zhang, L.; Haymet, A. D. J. J. Am. Chem. Soc. 1992, 114, 5875–5876. (83) Kalra, A.; Garde, S.; Hummer, G. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 10175–10180. (84) Olson, D. H.; Kokotailo, G. T.; Lawton, S. L.; Meier, W. M. J. Phys. Chem. 1981, 85, 2238–2243.

4075


A Computational Investigation into the Suitability of Purely Siliceous

Recommend Documents