Insights from Nonequilibrium Molecular Dynamics - ACS Publications

Subscriber access provided by ORTA DOGU TEKNIK UNIVERSITESI KUTUPHANESI

Article

Optimising water transport through graphene-based membranes: Insights from non-equilibrium molecular dynamics Jordan Muscatello, Frederike Jaeger, Omar K. Matar, and Erich A Muller ACS Appl. Mater. Interfaces, Just Accepted Manuscript • DOI: 10.1021/acsami.5b12112 • Publication Date (Web): 28 Apr 2016 Downloaded from http://pubs.acs.org on April 29, 2016

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

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

Page 1 of 24

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

ACS Applied Materials & Interfaces

Optimising water transport through graphene-based membranes: Insights from non-equilibrium molecular dynamics Jordan Muscatello,∗,† Frederike Jaeger,‡ Omar K. Matar,† and Erich A. Müller† †Department of Chemical Engineering, Imperial College London, SW7 2AZ ‡Department of Physics, Imperial College London, SW7 2AZ E-mail: [email protected]

Abstract Recent experimental results suggest that stacked layers of graphene oxide exhibit strong selective permeability to water. To construe this observation the transport mechanism of water permeating through a membrane consisting of layered graphene sheets is investigated via non-equilibrium and equilibrium molecular dynamics simulations. The effect of sheet geometry is studied by changing the offset between the entrance and exit slits of the membrane. The simulation results reveal that the permeability is not solely dominated by entrance effects; the path traversed by water molecules has a considerable impact on the permeability. We show that contrary to speculation in the literature, water molecules do not pass through the membrane as a hydrogen-bonded chain; instead, they form well-mixed fluid regions confined between the graphene sheets. The results of the present work are used to provide guidelines for the development of graphene and graphene oxide membranes for desalination and solvent separation.

1

ACS Paragon Plus Environment


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

Keywords graphene, membranes, separation, permeation, molecular dynamics, water, confined fluids

1

Introduction

As the recent discovery of lower dimensional allotropes of carbon (fullerenes and nanotubes) generated an abundance of research around their applications, the isolation of graphene has sparked an interest in using these materials as the basis for thin membranes capable of remarkable separations. Most relevant to this manuscript are the recent papers by Geim et al. 1 that claim that thin film graphene oxide (GO) sheets are capable of a high water permeance while presenting an extremely restricted porosity. Subsequent papers have confirmed these experimental observations 2 and the field is now open for new technology-driven applications of this phenomenon, particularly in the area of membrane separations. Here one can envisage membranes with an active layer made entirely from GO 3,4 or otherwise the incorporation of GO flakes within a membrane material 5,6 with the aim of water desalination and/or solvent purification. A proposed explanation of the apparent paradox raised by the above-mentioned heliumtight but water-permeable membrane is that water forms a linear “conga line” consisting of a hydrogen bonded network which stretches from the entrance to the exit of the membrane, 7 resulting in a collective motion of a single monolayer of water molecules. This view, presented in a schematic way in the seminal papers, 1 has been reproduced in further publications, 8 and is becoming the de facto most plausible explanation. However, other works question this mechanism. 9 The present work revisits some aspects of this view and employs molecular dynamics (MD) simulations of water permeation through GO-like sheets in order to gain further insight into the mechanisms underlying the observed high water permeance. Previous studies of water transport in confined nanogeometries have presented surprising results. Most notable is the observation that water molecules can flow within carbon nan2


Page 2 of 24

Page 3 of 24

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


otubes at rates which surpass by orders of magnitude predictions based on Hagen-Poiseuille flow. 10,11 Such transport behaviour can be attributed to the occurance of slip flow at the hydrophobic carbon walls as well as confinement effects on the density distribution of water molecules. This phenomenon has since been observed experimentally. However, a microscopic understanding has been predominantly provided by simulation studies. 12 A further pertinent element of discussion is the quantification of the energetic barriers required for a water molecule to remove itself from the bulk state and enter into a nanoconfined system. For round entrances to nanotubes or perforated graphene sheets, these energy barriers can be significant and have to be considered when accounting for the use of nanoporous carbons as membrane materials. 13 Modification of these barriers is fundamental to maximizing the efficiency of transport and this mechanism is key to the ion selectivity of nanochannels in biological membranes. 14 In contrast to the observations made of water permeation, the MD studies of gas diffusion through GO layers do not seem to show any particularly remarkable results, 15 and it seems understandable that the low porosity of typical GO membranes allows only for Knudsen-type diffusion, which is extremely limited. We present here a study of the flow of water through a GO-like membrane with an appropriate geometry, i.e. a net water flow normal to the plane of graphene flakes. The results of the permeance for different geometries, calculations of the potential of mean force and of the average mean number of hydrogen bonds paint a picture of the transport which is at odds with the current explanations available. Specifically, we show that energy barriers exist which coincide with the the breaking and forming of hydrogen bonds during the passage of water through the membrane and that significant mixing of water molecules occurs within the membrane. The rest of this paper is organised as follows. The methodology is outlined in Section 2, a discussion of the results is presented in Section 3, and concluding remarks are provided in Section 4.

3



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

2

Page 4 of 24

Methodology

2.1

Boundary-driven non-equilibrium molecular dynamics

In this work we use a boundary-driven non-equilibrium molecular dynamics (NEMD) method to establish a constant pressure difference, and hence maintain a steady state mass flux, across a graphene membrane. For this purpose a force along the z-axis is applied to all molecules located in a thin slab of width dext = 10.0˚ A approximately half the system cell length away from the membrane. In doing so a density gradient is set up in this layer, resulting in a net mass flux in the simulation cell and a corresponding pressure drop across the membrane. In these simulations the force is applied only to the oxygen sites of each SPC/E water molecule. Increasing the applied force results in a proportional increase in pressure drop, thus the linearity between the pressure drop and the resultant mass flux can be investigated. This method has been thoroughly described in previous studies of mass transport across porous media. 16–18 The externally applied force is specified as a multiple of the ratio between the and σ Lennard-Jones parameters for the oxygen atom of an SPC/E water molecule 19 and will be denoted in such units, i.e.

Fext = A/σ,

(1)

where we have chosen A to vary from 0.2 to 1.0. In the application of the externally applied force, the system undergoes a constant energy input. So as to perturb the dynamics of the water molecules as little as possible, a Berendsen thermostat 20 was applied to the graphene sheet only, such that the membrane dissipates excess energy. This is preferred to the case of applying a Nosé-Hoover or Langevin thermostat to the whole of the bulk fluid, which may artifically affect the dynamics, especially if the streaming velocity of the particles is not taken into account. 12,21 In order to promote energy dissipation in the membrane the carbon atoms of the graphene sheet are tethered to their

4


Page 5 of 24

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


2 nearest neighbours via harmonic bonds with spring constant kC-C = 478.4kcal/mol/˚ A equal

to the bond strength for aromatic carbon. 22 The thermostat was set to 300K. In this work the usual velocities-Verlet algorithm is used to integrate the equations of motion of the particles. The simulations were performed with the LAMMPS parallel molecular dynamics package. 23

2.2

Computational details and system geometry

In all simulations an orthogonal system cell with extents (Lx , Ly , Lz ) = (18.4, 23.0, 105.0)˚ A was used with periodic boundary conditions applied in all directions. The membrane was formed of graphene sheets with a slit-like opening placed in the centre of the cell in the xy plane. Membranes consisting of both a single and double layer of graphene were simulated. As indicated in Figure 1, the centre-to-centre separation of carbon atoms between adjacent sheets was fixed at 6˚ A. This value is commensurate with other computational studies of hydrated graphene oxide using reactive force fields. 24,25 The carbon-carbon opening of the slit, dslit , was initally fixed at 6˚ A. The free edges of the graphene sheet were functionalised with hydrogen atoms, using the parameters for aromatic groups from Mooney et al. 22 Water molecules were represented using the SPC/E model. 19 The interactions between water and carbon were modelled using oxygen-carbon Lennard-Jones parameters fit to give the macroscopically observed wetting behaviour on graphitic surfaces. 26 A timestep δt = 2fs was used for all simulations. The SHAKE algorithm was used to constrain the bonds and angle of the SPC/E molecule. 27 Electrostatic interactions were evaluated using the modified Wolf method proposed by Fennel and Gezelter whereby the k -space part of the Ewald summation is not computed but instead replaced by a damping function. 28 This has been observed to yield equation of state behaviour in fluid phases comparable to the full Ewald sum. 29 Results obtained from non-equilibrium simulations were produced from averages over runs of 10-16ns with a 2ns transient period where the system was allowed to reach a steady-state mass flux. 5



a)

dext 6Å

Fext y b)

doff 6Å

z

800 Fext = 0.2 Fext = 0.6 Fext = 1.0

600

P / MPa

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

Page 6 of 24

400 200 0 -200 -400 0

20

40

60

80

100

z /Å

Figure 1: The geometry of the double-sheet system (a) and resultant pressures profiles for various applied forces (b). The bottom panel was generated from the c = 0.0 simulation. The inter-sheet distance and carbon-carbon slit opening is fixed at 6˚ A. The dashed red lines in the bottom figure indicate the positions of the atoms comprising the graphene sheets.

6


Page 7 of 24

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


For the case of a double-sheet membrane the geometry of the system can be specified by changing the offset doff between the openings of the graphene sheets as a fraction of the cell length along the y-axis. In this work we consider three offsets doff = cLy , where c = {0.0, 0.25, 0.5}, to study the effect of increasing the relative distance between the openings in the sheet. Each carbon atom in the membrane is tethered to a fixed point in space by 2

harmonic bonds with spring constant ktether = 100kcal/mol/˚ A such that the membrane does not change overall position in the cell when a flux is applied to the water molecules. This emulates a membrane held in place by an external support, in that we do not consider the mechanical properties of the membrane under applied pressure in this work. The doublesheet system is illustrated in Figure 1 along with the resultant hydrostatic pressure profiles for varying applied forces in a system with zero offset, c = 0.0, between the graphene sheets. The pressure profile was calculated by taking the average of the diagonal components of the Irving-Kirkwood expression for the pressure tensor. 30 The transmembrane pressure drop was calculated using the difference in hydrostatic pressure of the bulk fluid on either side of the membrane.

3 3.1

Results Permeability

Figure 2 shows the molar flux of water molecules versus the transmembrane pressure drop for three offset geometries for the double-sheet system as well as for a single sheet with dslit fixed at 6˚ A. We assume that a linear relationship holds between the molar flux and the pressure drop for this range of applied pressures. In making this assumption a relationship resembling Darcy’s law can be applied, 31,32

Jz = k

7

∆P lmem


(2)


30

3

c = 0.0 c = 0.25 c = 0.5 Single Sheet

20

2

Jz / mol/m s × 10

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

Page 8 of 24

10

0 0

20

40

60 80 ∆P / MPa

100

120

Figure 2: Molar flux of water through membranes with increasing offset values against transmembrane pressure drop. Symbols are simulations results. Dashed lines are a linear fitting. Errors bars represent the standard error.

8


Page 9 of 24

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


where Jz is the molar flux through the membrane, ∆P is the transmembrane pressure drop, k is the coefficient of permeability and lmem is considered to be an effective membrane width. Here we define lmem to be the region over which the pressure drop occurs, taken from the density profile as the distance from one bulk region to another. These peaks are clearly visible in Figure 1. This is in accordance with the definition of the width of a single graphene sheet used by other authors. 33 In doing so we can extract values for the permeability of each system geometry from the slope of the linear fitting to the data. It should be noted that in this case the linear expression is purely empirical, as hydrodynamic arguments break down at pore sizes of molecular scale. 34 The permeability values obtained are shown in Table 1. It can be seen that the permeability of two sheets decreases with increasing sheet offset. This dependence shows that the permeability depends not only on the porosity of the membrane (which is fixed in this case) but also on the length of the path taken by the water molecules. The permeability for a single sheet was found to be approximately three times that of the double-sheet with inline slits. Care should be taken in comparing the permeabilities for the single and double slit membranes. Permeability is expressed as a material property of homogeneous membranes (normalised by the width of the membrane). However due to the discrete number of sheets, the single and double sheets cannot be considered to have equivalent material properties in this case. An additional energetic penalty is present in the double sheet due to the presence of the monolayer of water (see discussion below). In the limit of a large number of sheets a graphene based membrane will appear homogeneous at the scale of the total membrane and the permeability should reach a constant value. In this study we do not attempt to find this limit.

9



Table 1: Table of effective membrane widths (lmem ), permeabilities (k), mean residence times (τres ), mean hydrogen bond lifetimes (τHB ), and diffusion coefficients (Dmono ), for all double slit simulations and if applicable also for a single sheet. To obtain permeability values in units of m2 s−1 Pa−2 , the molar flux was converted to a volumetric flux using the density of water at standard temperature and pressure.

PMF(kJ/mol)

c = doff /Ly 0.0 0.25 0.5 Single Sheet

lmem (˚ A) 34.0 34.0 34.0 28.0

k (×10−17 m2 s−1 Pa−1 ) 0.35 ± 0.05 0.22 ± 0.03 0.08 ± 0.02 1.02 ± 0.05

3 2

nHB

τres (ps) 13.7 ± 1.6 12.0 ± 0.7 17.3 ± 1.5 -

Dmono (×10−5 cm2 s−1 ) 0.65 ± 0.21 0.72 ± 0.21 0.73 ± 0.21 -

τHB (ps) 1.86 ± 0.10 2.28 ± 0.10 2.19 ± 0.10 -

(a)

(b)

(c)

(d)

(e)

(f )

1 0 −1 3.6 3.4 3.2 3.0 2.8

ρ (g/cm3 )

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

Page 10 of 24

3.0 2.5 2.0 1.5 1.0 0.5 0.0 35

40

45

50 ˚ z(A)

55

60

65

35

40

45

50

55

60

65

70

˚ z(A)

Figure 3: PMF (panels (a) and (b)), hydrogen bonds per molecule (panels (c) and (d )), and density (panels (e) and (f )) profiles for single sheet (left) and double-sheet (right, c = 0.0) systems with dslit = 6˚ A. The positions of the graphene sheets are indicated by dashed red lines. Free energy barriers coincide with density peaks due to molecular ordering at the sheet entrances and the breaking of hydrogen bonds as molecules pass through the membrane. Error bars on the PMFs were obtained via a Monte Carlo bootstrap algorithm 35

3.2

Potentials of Mean Force (PMF)

To investigate the energetics of water molecules entering the membrane, the potential of mean force (PMF) was calculated along a reaction coordinate which crosses perpendicularly

10


Page 11 of 24

˚ Distance from center (A)

8 7

PMFbarrier (kJ/mol)

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


6 5 4 3 2

6

˚ dslit = 5 A ˚ dslit = 6 A

4

˚ dslit = 9 A ˚ dslit = 15 A

2 0 −2 −4 −6 0

1

2 3 4 3 ρ (g/cm )

10

12

5

6

1 0 4

6

8

˚ dslit (A)

14

16

Figure 4: Value of potential of mean force at the center of a slit for a range of slit widths in a single-sheet system. The inset plot shows the density profile in the y-direction as a function of distance from the slit center. the entrance to the graphene sheet. An umbrella sampling method in the canonical ensemble was used for this purpose. 36 The Nosé-Hoover thermostat was used at T = 300K with a damping parameter µdamp = 100fs. A molecule is chosen at random and attached to a harmonic bias potential, with equilibrium position at a point along the z-axis, aligned with the centre of the slit opening. The harmonic potential is described by Vspring = 12 k(z − z0 )2 , where k is the spring constant and z0 the equilibrium position. After an equilibration period, the position distributions were collected along the reaction coordinate in steps of 0.2˚ A over a period of 6ns using a spring constant ranging from 5.0-12.6 kJ/mol-˚ A2 . The data was then analysed using the weighted histogram analysis method (WHAM). 35 Figure 3 (panels (a) and (b)) shows the PMFs for the single sheet and a zero offset doublesheet membrane with reference to the bulk with slit width of 6˚ A. In both cases free energy

11



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

barriers are shown to emerge at the entrance to the membrane. When compared to the density profile (Figure 3, panels (e) and (f )) it can be seen that these barriers coincide with molecular adsorption on the graphene sheets. For the case of the double-sheet the energetic barrier at the slit is significantly higher than for the single sheet. This can be interpreted as being due to the presence of a monolayer of water in the intermembrane region exhibiting a different structure and higher density from the bulk. Such high density configurations have been frequently reported in confined systems. 37 PMFs were also calculated for single sheets with varying slit widths. Figure 4 shows the height of the energy barrier at the slit for dslit = {5, 6, 9, 15}˚ A. In general, the barrier increases with decreasing slit width. For larger slit widths confinement effects become less prevalent indicating a return to bulk-like behaviour, as shown in the inset plot illustrating the density distribution within the slit. As the slit width is increased, a transition from a single monolayer to a layer at each edge is observed. This is analogous to the ring distributions of water molecules observed in carbon nanotubes. 38 Density peaks appear in quantitatively similar positions.

3.3

Hydrogen bonding profile

The degree of hydrogen bonding was calculated along the z-axis. A hydrogen bond is considered to exist between two molecules if the following distance and orientation criteria are met: ROO < 3.4˚ A, ROH < 2.425˚ A, and the H-O..O angle is less than 30◦ , where ROO and ROH are the radial oxygen-oxygen and oxygen-hydrogen distances respectively on neighbouring molecules. 39,40 The centre panels in figure 3 ((c) and (d)) show the average number of hydrogen bonds per molecule calculated in slabs along the simulation cell for both a single and a double-sheet membrane. The average number of hydrogen bonds in the bulk is ≈ 3.5 in both cases. The degree of hydrogen bonding drops rapidly at the entrance to the membrane, indicating that there is an energetic penalty incurred as a molecule enters the membrane; hydrogen bonds 12


Page 12 of 24

Page 13 of 24

must be broken as a molecule passes through the slit. This provides a physical explanation of the free energy barriers shown in the top panel. For the case of the double-sheet the degree of hydrogen bonding in the region between the sheets approaches that found in the bulk. This indicates that even though the width of the spacing can only support a monolayer of water and is therefore too narrow to support a fully tetrahedrally ordered network hydrogen bonding is still present. The relative increase in the height of the free energy barrier at the slit compared to the single sheet case can be interpreted as due to an increased disruption and subsequent rearrangement of the hydrogen bonded network in order to form the monolayer.

3.4

Residence times and water molecule dynamics

0.6 0.5

p(tres)

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


c = 0.0 c = 0.25 c = 0.5

1 pHB

0.4

0.1

0.3

0.01

0.2

0.001 1

10

100

0.1 0 0

10

t res / ps

20

30

Figure 5: Normalised residence time distributions for water molecules in the double-sheet membrane with different offset values. The inset plot shows the same data on log-log axes, as well as the hydrogen bond lifetime distribution for c = 0.0. The residence time distribution of water molecules in the membrane for the double-sheet system was calculated for each sheet offset value. The times were logged only for molecules contributing to the net mass flux, i.e. the residence time tres = tout − tin is logged only if a 13



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

molecule enters the membrane from the left sheet and subsequently exits through the right sheet (with respect to the direction of flux). The normalised residence time distributions are shown in figure 5. In all cases the distribution decays approximately as a power law, as indicated by the inset log-log plot. This shows that although a water molecule is more likely to spend a short period of time in the membrane, some molecules can remain diffusing in the monolayer. The mean residence time (see table 1) for the c = 0.5 case is approximately 25% longer than the c = 0.25 case. This may be surprising, as the average path length for the offset cases, assuming an equal probability of a molecule taking the most direct route out of the membrane, is the same, since in the c = 0.25 offset there is also a longer path the molecules can take (0.75 × Ly ). Analogous to the case of fluid elements in a chemical reactor, the residence time distribution reveals information regarding the degree of mixing of water molecules in the membrane. 41 For the hypothesised case of a single hydrogen bonded sheet of water molecules moving collectively through the membrane the flow would resemble a ‘plug-flow’-like scenario in which the residence time distribution is single-valued, i.e. a delta function. Whilst there is a strong peak for shorter residence times, the distributions all feature a slowly decaying tail. This indicates that a significant contribution is made by molecules remaining in the membrane for longer times. An animated visualisation of this diffusive behaviour is higlighted in the Supporting Information. The inset log-log plot in figure 5 shows the residence time distributions compared to the hydrogen bond lifetime distribution for the c = 0.0 case. At short timescales the distributions appear similar. However, the hydrogen bond lifetime does not feature the long decay of the residence time distribitions. This becomes apparent when the mean values are compared (see table 1), since the mean residence times are greater than the hydrogen bond lifetime by a factor of at least five, despite the similarity of the distributions at short times. The hydrogen bond lifetime and associated probability distribution function were calculated for each offset for NEMD simulations with A = 1.0. The hydrogen bond probability

14


Page 14 of 24

Page 15 of 24

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


distribution was computed for water molecules located in the membrane in accordance with Han et al., 42 using the same criteria for the existence of a hydrogen bond as mentioned in section 3.3. The hydrogen bond lifetime is given by 42 ∞

Z

tpHB (t)dt

τHB =

(3)

0

where pHB (t) is the probability distribution. Figure 6 shows pHB for c = 0.25. The values for τHB are given in table 1. These values are small compared to the average residence times, suggesting that during the time that a water molecule remains in the membrane, the hydrogen bonded network is undergoing continuous breaking and reforming. To further investigate the dynamic behaviour of water molecules within the monolayer the velocity autocorrelation function (VACF) was calculated for NEMD runs with A = 1.0 in the xy plane such that only the diffusive motion of water molecules in the monolayer was taken into account, i.e. only the vx and vy components of the velocity were used. Figure 6 shows the results for the c = 0.25 case and that of bulk water, where the normalised VACF is given by Ψ(t) =

hv(0) · v(t)i h|v(0)|2 i

(4)

and vx v= vy for molecules located between the carbon sheets. Integrating the non-normalised VACF yields the diffusion coefficient for molecules diffusing in the xy plane, given by

Dmono

1 = d

Z 0

∞

hv(0) · v(t)idt

(5)

where d = 2 for the case of diffusion in a plane. Comparing the VACFs of confined and bulk water, it can be seen that there is a stronger negative correlation for the confined molecules than those of the bulk. This is analogous 15



to the VACF observed in amorphous solids or liquids below the glass transition. 43 The diffusion coefficient was calculated from the VACFs for all offsets, yielding values of ≈ 0.7 × 10−5 cm2 s−1 . This value is lower than that of bulk water by a factor of five, rather than orders of magnitude lower as would be expected for glassy system, indicating that liquid-like diffusion is still prevalent. Furthermore, the timescale associated with diffusion in the sheet is much shorter than the average residence time. Therefore significant diffusive mixing is to be expected for a molecule remaining in the sheet. 1.5

1 c = 0.25 Bulk water

1 pHB(t)

0.5 Ψ(t)

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

Page 16 of 24

0.5 0

0

100

200

300

400

0 0

500

2

t / fs

(a)

4

t / ps

6

8

10

(b)

Figure 6: a) Normalised velocity autocorrelation function for water molecules located between carbon sheets for c = 0.25 and for bulk water. b) Hydrogen bond lifetime distribution for water molecules located between carbon sheets for c = 0.25.

4

Conclusions

In this study we have observed the permeation of water through stacked graphene layers via non-equilibrium molecular dynamics simulation. In particular, we investigated the effect of slit widths and sheet geometries on the transport of water molecules. For narrow entrances there are energetic barriers on either side of the membrane due to molecular layering of the adsorbed water as well as a barrier associated with the passage through the opening, all of which serve to hinder transport. As the slit width is increased bulk behaviour is recovered, where no barriers are present. In addition, there is a further energetic penalty associated with the formation of a confined phase between multiple sheets. This is due to molecular 16


Page 17 of 24

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


rearrangement that must occur when a molecule passes from the bulk to the monolayer. Hydrogen bonds are broken and then reformed as a molecule enters the membrane. Results readily reveal that the sheet offset also has a considerable effect on the resultant permeability of the membrane, whereby it is favourable for the slits to be as close as possible in adjacent stacking layers. The previously proposed mechanism for confined water transport through graphene-based membranes involves the continuous plug-like flow of a hydrogen bonded sheet of water. Other authors have shown alternative mechanisms for fast flow in confined geometries that support more than one layer of water. 9 In addition we have shown that for highly confined single monolayer water the previously proposed description is not accurate. The residence time analysis instead reveals a high degree of mixing as molecules permeate through the membrane, supported by analysis of the dynamical properties of the hydrogen bonded monolayer. Analysis performed on an additional system permitting only a single route through the membrane indicates that the above conclusions are not an effect of the periodicity of the simulations. These results are presented in the supporting information. The above conclusions suggest clear guiding principles for the design of graphene-based membranes in terms of maximising the transport of water, whilst still maintaining the potential for use in separation processes. This involves minimising the energy barriers upon entrance to the membrane. We have shown that unimpeded entry to the sheets can be achieved by using slit widths greater than ≈ 3 molecular diameters. It has been suggested that tunable sheet spacings may be achieved by covalently bonded “linkers”. 44,45 Rejection properties may also be manipulated in the same way by ensuring that the spacing between sheets only admits a single monolayer of water. In summary, by optimising the configuration at the molecular scale, the viability of graphene oxide as a potential material for separations (e.g. desalination, alcohol dehydration etc.) can be maximised.

17



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

5

Supporting Information

Additional simulation results describing flow through a membrane only permitting a single path. A video highlighting a single molecule diffusing through the membrane with c = 0.5.

Acknowledgement The authors acknowledge funding from the Thomas Young Centre under grant number TYC101 and EPSRC grants EP/E016340 and EP/J014958. FJ was supported through a studentship in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London funded by the Engineering and Physical Sciences Research Council under grant number EP/G036888/1. FJ thanks Dr. Carlos Braga for advice regarding the calculation of PMF profiles. JM thanks Prof. Adrian Sutton and Dr. Arash Mostofi for comments and suggestions.

References (1) Nair, R. R.; Wu, H. A.; Jayaram, P. N.; Grigorieva, I. V.; Geim, A. K. Unimpeded Permeation of Water Through Helium-Leak-Tight Graphene-Based Membranes. Science 2012, 335, 442–444. (2) Kim, H. W.; Yoon, H. W.; Yoon, S.-M.; Yoo, B. M.; Ahn, B. K.; Cho, Y. H.; Shin, H. J.; Yang, H.; Paik, U.; Kwon, S.; Choi, J.-Y.; Park, H. B. Selective Gas Transport Through Few-Layered Graphene and Graphene Oxide Membranes. Science 2013, 342, 91–95. (3) Cohen-Tanugi, D.; Grossman, J. C. Nanoporous Graphene as a Reverse Osmosis Membrane: Recent Insights from Theory and Simulation. Desalination 2015, 366, 59 – 70.

18


Page 18 of 24

Page 19 of 24

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


(4) Mahmoud, K. A.; Mansoor, B.; Mansour, A.; Khraisheh, M. Functional Graphene Nanosheets: The Next Generation Membranes for Water Desalination. Desalination 2015, 356, 208 – 225. (5) Xia, S.; Ni, M. Preparation of Poly(Vinylidene Fluoride) Membranes with Graphene Oxide Addition for Natural Organic Matter Removal. J. Membr. Sci. 2015, 473, 54 – 62. (6) Goh, K.; Setiawan, L.; Wei, L.; Si, R.; Fane, A. G.; Wang, R.; Chen, Y. Graphene Oxide as Effective Selective Barriers on a Hollow Fiber Membrane for Water Treatment Process. J. Membr. Sci. 2015, 474, 244 – 253. (7) M¨ uller, E. A. Purification of Water Through Nanoporous Carbon Membranes: A Molecular Simulation Viewpoint. Curr. Opin. Chem. Eng. 2013, 2, 223 – 228. (8) Liu, G.; Jin, W.; Xu, N. Graphene-based Membranes. Chem. Soc. Rev. 2015, 44, 5016– 5030. (9) Wei, N.; Peng, X.; Xu, Z. Understanding Water Permeation in Graphene Oxide Membranes. ACS Appl. Mater. Interfaces 2014, 6, 5877–5883. (10) Hummer, G.; Rasaiah, J. C.; Noworyta, J. P. Water Conduction through the Hydrophobic Channel of a Carbon Nanotube. Nature 2001, 414, 188–190. (11) Majumder, M.; Chopra, N.; Andrews, R.; Hinds, B. J. Nanoscale Hydrodynamics: Enhanced Flow in Carbon Nanotubes. Nature 2005, 438, 44. (12) Thomas, M.; Corry, B. Thermostat Choice Significantly Influences Water Flow Rates in Molecular Dynamics Studies of Carbon Nanotubes. Microfluid. Nanofluid. 2014, 18, 41–47. (13) Cohen-Tanugi, D.; Grossman, J. C. Water Desalination across Nanoporous Graphene. Nano Lett. 2012, 12, 3602–3608. 19



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

(14) Agre, P.; Preston, G. M.; Smith, B. L.; Jung, J. S.; Raina, S.; Moon, C.; Guggino, W. B.; Nielsen, S. Aquaporin CHIP: The Archetypal Molecular Water Channel. Am. J. Physiol.: Renal, Fluid Electrolyte Physiol. 1993, 265, F463–F476. (15) Jiao, S.; Xu, Z. Selective Gas Diffusion in Graphene Oxides Membranes: A Molecular Dynamics Simulations Study. ACS Appl. Mater. Interfaces 2015, 7, 9052–9059. (16) Zhu, F.; Tajkhorshid, E.; Schulten, K. Pressure-Induced Water Transport in Membrane Channels Studied by Molecular Dynamics. Biophys. J. 2002, 83, 154 – 160. (17) Frentrup, H.; Avenda˜ no, C.; Horsch, M.; Salih, A.; M¨ uller, E. A. Transport Diffusivities of Fluids in Nanopores by Non-Equilibrium Molecular Dynamics Simulation. Mol. Simul. 2012, 38, 540–553. (18) Huang, C.; Choi, P. Y. K.; Kostiuk, L. W. A Method For Creating a Non-Equilibrium NT(P1 - P2) Ensemble in Molecular Dynamics Simulation. Phys. Chem. Chem. Phys. 2011, 13, 20750–20759. (19) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269–6271. (20) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684–3690. (21) Bernardi, S.; Todd, B. D.; Searles, D. J. Thermostating Highly Confined Fluids. J. Chem. Phys. 2010, 132, 244706–10. (22) Mooney, D.; M¨ uller-Plathe, F.; Kremer, K. Simulation Studies for Liquid Phenol: Properties Evaluated and Tested over a Range of Temperatures. Chem. Phys. Lett. 1998, 294, 135 – 142.

20


Page 20 of 24

Page 21 of 24

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


(23) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1 – 19. (24) Compton, O. C.; Cranford, S. W.; Putz, K. W.; An, Z.; Brinson, L. C.; Buehler, M. J.; Nguyen, S. T. Tuning the Mechanical Properties of Graphene Oxide Paper and Its Associated Polymer Nanocomposites by Controlling Cooperative Intersheet Hydrogen Bonding. ACS Nano 2012, 6, 2008–2019. (25) Medhekar, N. V.; Ramasubramaniam, A.; Ruoff, R. S.; Shenoy, V. B. Hydrogen Bond Networks in Graphene Oxide Composite Paper: Structure and Mechanical Properties. ACS Nano 2010, 4, 2300–2306. (26) Werder, T.; Walther, J. H.; Jaffe, R. L.; Halicioglu, T.; Koumoutsakos, P. On the WaterCarbon Interaction for Use in Molecular Dynamics Simulations of Graphite and Carbon Nanotubes. J. Phys. Chem. B 2003, 107, 1345–1352. (27) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23, 327 – 341. (28) Fennell, C. J.; Gezelter, J. D. Is the Ewald Summation Still Necessary? Pairwise Alternatives to the Accepted Standard for Long-range Electrostatics. J. Chem. Phys. 2006, 124 . (29) Muscatello, J.; Bresme, F. A comparison of Coulombic Interaction Methods in NonEquilibrium Studies of Heat Transfer in Water. J. Chem. Phys. 2011, 135, 234111–8. (30) Irving, J. H.; Kirkwood, J. G. The Statistical Mechanical Theory of Transport Processes. IV. The Equations of Hydrodynamics. J. Chem. Phys. 1950, 18, 817–829. (31) Roque-Malherbe, R. M. Adsorption and Diffusion in Nanoporous Materials; CRC Press, 2007. 21



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

(32) Mulder, M. Basic Principles of Membrane Technology; Kluwer Academic Publishers, 1996. (33) Suk, M. E.; Aluru, N. R. Molecular and Continuum Hydrodynamics in Graphene Nanopores. RSC Adv. 2013, 3, 9365–9372. (34) Falk, K.; Coasne, B.; Pellenq, R.; Ulm, F.-J.; Bocquet, L. Subcontinuum Mass Transport of Condensed Hydrocarbons in Nanoporous Media. Nat. Commun. 2015, 6 . (35) Grossfield, A. WHAM: An Implementation of the Weighted Histogram Analysis Method. (36) Roux, B. The Calculation of the Potential of Mean Force Using Computer Simulations. Comput. Phys. Commun. 1995, 91, 275–282. (37) Striolo, A. From Interfacial Water to Macroscopic Observables: A Review. Adsorpt. Sci. Technol. 2011, 29, 211–258. (38) Wang, L.; Dumont, R. S.; Dickson, J. M. Nonequilibrium Molecular Dynamics Simulation of Water Transport Through Carbon Nanotube Membranes at Low Pressure. J. Chem. Phys. 2012, 137, 044102. (39) Guàrdia, E.; Mart´ı, J.; Garc´ıa-Tarrés, L.; Laria, D. A Molecular Dynamics Simulation Study of Hydrogen Bonding in Aqueous Ionic Solutions. J. Mol. Liq. 2005, 117, 63 – 67. (40) Luzar, A.; Chandler, D. Structure and Hydrogen Bond Dynamics of Water–Dimethyl Sulfoxide Mixtures by Computer Simulations. J. Chem. Phys. 1993, 98, 8160–8173. (41) Danckwerts, P. Continuous Flow Systems: Distribution of Residence Times. Chem. Eng. Sci. 1953, 2, 1–13. (42) Han, S.; Kumar, P.; Stanley, H. E. Hydrogen-bond Dynamics of Water in a Quasi-twodimensional Hydrophobic Nanopore Slit. Phys. Rev. E 2009, 79, 041202. 22


Page 22 of 24

Page 23 of 24

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


(43) Rahman, A.; Mandell, M. J.; McTague, J. P. Molecular Dynamics Study of an Amorphous LennardJones System at Low Temperature. J. Chem. Phys. 1976, 64, 1564–1568. (44) Mi, B. Graphene Oxide Membranes for Ionic and Molecular Sieving. Science 2014, 343, 740–742. (45) Hu, M.; Mi, B. Enabling Graphene Oxide Nanosheets as Water Separation Membranes. Environ. Sci. Technol. 2013, 47, 3715–23.

23


ACS Applied Materials & Interfaces Page 24 of 24

1 2 3 4 5


Insights from Nonequilibrium Molecular Dynamics - ACS Publications

Recommend Documents