Molecular Dynamics Study of Water Flow across Multiple Layers of

Molecular Dynamics Study of Water Flow across Multiple Layers of Pristine, Oxidized, and Mixed Regions of Graphene Oxide Jon A. L. Willcox† and Hyung J. Kim*,†,‡,§ †

Department of Chemistry, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, United States School of Computational Sciences, Korea Institute for Advanced Study, Seoul 02455, Korea

‡

ABSTRACT: A molecular dynamics graphene oxide model is used to shed light on commonly overlooked features of graphene oxide membranes. The model features both perpendicular and parallel water flow across multiple sheets of pristine and/or oxidized graphene to simulate “brickand-mortar” microstructures. Additionally, regions of pristine/oxidized graphene overlap that have thus far been overlooked in the literature are explored. Differences in orientational and hydrogen-bonding features between adjacent layers of water in this mixed region are found to be even more prominent than differences between pristine and oxidized channels. This region also shows lateral water flow in equilibrium simulations and orthogonal flow in non-equilibrium simulations significantly greater than those in the oxidized region, suggesting it may play a non-negligible role in the mechanism of water flow across graphene oxide membranes. KEYWORDS: graphene oxide, membranes, molecular dynamics, water flow, nanoconfined water nomalous rapid water flow across graphene oxide (GO) membranes as noted in the work by Nair et al.1 has caught the attention of experimentalists looking to improve filtration and sieving techniques for liquids and gases2−5 and has inspired a series of theoretical studies to elucidate the mechanism of the unusual transport.6−10 Understanding and optimizing the system’s features for selective permeation may be key to enabling its practical application in areas like water filtration, solvent purification, and gas separation.1−3,11−15 The difficulty in doing so with computational methods lies in the material’s nonuniform nature. Micro- and nanoscale heterogeneity like that found in graphene oxide can severely complicate molecular dynamics investigationsa technique typically used for nanosecond calculations in nanometer-scaled worlds. A variety of perspectives are needed to approach a comprehensive understanding of the roles different structural features play in water flow across the GO membrane. In their initial study, Nair et al. proposed that regions of overlapping unoxidized or “pristine” graphene facilitate rapid water permeation by nearly frictionless flow analogous to the capillary motion in carbon nanotubes, while regions of overlapping oxidized graphene contribute to the membrane structure.1 Some follow-up studies maintain these designations,2,6 whereas others have expressed doubt that capillary action is the soleor even primarycontributor to

A

© 2017 American Chemical Society

the water transport mechanism. Given the nonuniform distribution and small area of many unfunctionalized regions, the likelihood of contiguous cross-membrane graphitic channels is dramatically reduced as GO layers are added.16−18 Furthermore, drag generated in pristine regions by water interactions with adjacent oxidized regions could significantly reduce the flow rate.7,9 Alternative (or supplemental, as the case may be) theories emphasizing the role of defects, such as pores, slits, and wrinkles, have recently been gaining attention.12,13 The current literature has also placed an emphasis on comparisons between pristine (P) channels and oxidized (O) channels. Here, we contend that even with the tendency for polar−polar and nonpolar−nonpolar alignment, the heterogeneous nature of GO likely leads to regions of misalignment (i.e., mixed regions). In this study, we introduce a multilayer “brick-and-mortar” GO model to compare pristine, oxidized, and mixed membranes, incorporating in-plane and cross-plane water flow made possible by staggered slits in the membrane layers, as shown in Figure 1. The nomenclature used throughout this work is detailed in Table 1. Received: December 21, 2016 Accepted: January 20, 2017 Published: January 20, 2017 2187

DOI: 10.1021/acsnano.6b08538 ACS Nano 2017, 11, 2187−2193

Article

www.acsnano.org

Article

ACS Nano

Figure 1. PPOO membrane side view with slitswater is shown in white, carbon in purple, and hydroxyl groups in blue. Periodic boundary conditions extend in all directions.

Figure 2. Density of water (calculated from oxygen distribution) in the GO membrane with respect to bulk water density vs z-axis. Graphene layer types and positions are indicated to the right.

RESULTS AND DISCUSSION Water Layer Structure. The PPOO density plot in Figure 2 displays water bilayer structures for all regions. Semi- to highly ordered structures for water systems confined by hydrophobic walls that result in nearly discrete layers have been predicted theoretically9 and indicated experimentally.19,20 There is no definitive repeating structure for water confined between pristine graphene in our systemhexagonal, pentagonal, and rhombic clusters all coexist as is shown in Figure 3. This appears to be in agreement with the simulation results presented by Algara-Siller et al. for confined water under ambient conditions.20 There is a significant difference in density between oxidized and unoxidized regionsan imbalance that becomes even more prominent between the adjacent hydrolayers in the mixed region (MPWL = 22.91 nm−3 and MOWL = 18.36 nm−3)! This is believed to result from the disruption of intra-hydrolayer order in the presence of hydroxyl groups. Hydrogen-Bonding Structure. In order to measure the extent of hydrogen-bonded interactions in the system, a geometric hydrogen-bonding definition is introduced. To qualify as a hydrogen bond, the H···O distance was required to be less than 0.24 nm (distance was chosen to encompass the first solvation shell in the radial distribution functions and is the same for water−water and water−hydroxyl interactions; plots of radial distributions are not included), and the H−O···O angle was required to be less than 30°. A similar definition has been used for water systems in previous studies.10,21 The average numbers of water−water hydrogen bonds within a hydrolayer and between adjacent hydrolayers are displayed in

Figure 3. Snapshot of water in the pristine water layer (PWL) showing coexisting rhombic, pentagonal, and hexagonal clusters. Hydrogen-bonding interactions are shown in black.

Table 2. The numbers of hydrogen bonds to basal hydroxyl groups in the mixed and oxidized regions are listed in Table 3. As is expected due to competing interactions, water layers located adjacent to an oxidized graphene layer (O-graphene) show reduced intralayer hydrogen bonding. There is, however, an increase in the number of interlayer hydrogen bonds when one or both of the hydrolayers is adjacent to O-graphene. This increase is attributed in part to the reorientation of water molecules upon interacting with the hydroxyl groups. In the water layers bordering O-graphene (i.e., MOWL and OWL), ∼50% of water molecules form hydrogen bonds to hydroxyl

Table 1. Terminology Used Throughout This Work term

definition

pristine (P) graphene oxidized (O) graphene membrane PPPP OOOO POPO PPOO channel region hydrolayer PWL OWL MPWL MOWL

unfunctionalized graphene hydroxyl functionalized graphene four periodically repeating P/O-graphene sheets pristine membrane oxidized membrane mixed membrane membrane with alternating pristine and oxidized channels space between sheets of graphene type of graphene overlap (pristine, oxidized, or mixed) a monolayer of water hydrolayer in the pristine channel hydrolayer in the oxidized channel hydrolayer in the mixed channel adjacent to P-graphene hydrolayer in the mixed channel adjacent to O-graphene 2188


Article

ACS Nano

symmetry axis in Figure 4a. The other two axes are symmetric around θ = 90°. Some similarities appear between water orientations in the pristine and mixed pristine water layers (PWL and MPWL) as well as between the oxidized and mixed-oxidized water layers (OWL and MOWL). For example, peaks indicating that the water plane is oriented parallel or nearly parallel to the P/Ographene plane (i.e., 60° ≲ θ ≲ 120° in panels a and c and θ ≈ 0 or 180° in panel b) are more prominent for both of the cases in which the hydrolayer borders P-graphene. However, other aspects of the mixed region are unique. In Figure 4b, there is a central peak (∼60−120°) for the PWL representing the water planes oriented perpendicular to the graphene sheet; this peak is greatly reduced in the MPWL. For the OWL, Figure 4b shows a relatively uniform distribution, but for the MOWL, the central peak is more prominent than in all other regions. In other words, just as in the hydrogen-bonding case discussed above, there is a greater contrast between the layers in the mixed region than between the PWL and the OWL. To clarify the above findings, we correlate the water orientation to interlayer hydrogen bonding in Figure 5 by only plotting the orientation of water molecules that form interlayer hydrogen bonds. The plots distinguish between the molecules that act as hydrogen-bond donors (pale colors) and acceptors (dark colors). As is expected, the molecular plane of donor water molecules tends to be oriented perpendicular to the P/O-graphene plane. Surprisingly, however, the molecular plane of acceptor water has a considerable tendency to be oriented parallel to the P/O-graphene. The orientation of hydrogen-bonded pairs is most likely the cause of the central peak in Figure 4b. Forty-one percent of MPWL water molecules and 52% of MOWL water molecules act as hydrogen-bond donors, which could account for the difference in the respective peak heights. Forty-four percent of PWL water molecules and 49% of OWL water molecules act as hydrogenbond donors, but because these layers are adjacent to equivalent layers, there must be a 1:1 ratio of donors to acceptors. Therefore, orientation peaks corresponding primarily to hydrogen-bond acceptors (e.g., the peripheral peaks in Figure 4b) and donors (e.g., the central peak in Figure 4b) are approximately balanced for the PWL and OWL, whereas they are imbalanced for the MPWL and MOWL. This imbalance leads to more prominent peak differences between the layers of the mixed region when compared to the peak differences between the PWL and the OWL. Finally, we considered the orientation of water molecules in the regions above and below the slits (figures not shown). We found that the peaks that appear in Figure 4 remain present but

Table 2. Average Number of Water−Water Hydrogen Bonds within a Hydrolayer (Intra) and between Hydrolayers (Inter) per Water Moleculea intra inter intra inter

(channel) (channel) (slit) (slit)

PWL

MPWL

MOWL

OWL

2.49 0.89 2.64 0.71

2.52 0.91 2.64 0.74

1.77 0.91 2.22 0.74

1.86 0.98 2.23 0.77

The first two rows exclude water directly above and below the slit, and the second two rows only pertain to water directly above and below the slits.

a

Table 3. Average Number of Water−Hydroxyl Hydrogen Bonds Per Hydroxyl Group hydroxyl role

OWL

MWL

donor acceptor

0.53 0.69

0.56 0.71

groups, and 97% of the water molecules that are hydrogen bonded to hydroxyl groups were also found to form an interlayer hydrogen bond. Water molecules directly above and below the slits form fewer interlayer hydrogen bonds than those bordering the basal plane. The number of intralayer hydrogen bonds, however, is found to increase near the slits; water is able to reorient for more favorable interactions in the less constrained environment near the slits, disrupting the layer structure. The hydroxyl groups lining the slits were found to have minimal hydrogenbonded interactions with water. On average, there were ∼0.1 hydrogen bonds per slit hydroxyl group, primarily formed via water acting as a donor. Approximately 60% of slit hydroxyl groups formed hydrogen bonds with other slit hydroxyl groups. By contrast, a hydroxyl group away from the slit region forms 1.2−1.3 hydrogen bonds with water (Table 3). Water Orientation. The results for the probability distribution P(θ) of water orientations in the channels between slits are shown in Figure 4, where θ is the angle between the specified molecular axis and the z-axis, and P(θ) is normalized as

∫ P(θ)sin θ dθ = 1

(1)

The water orientations are defined such that at θ = 0°, the molecular axis vector points toward the adjacent water layer. Note that this distinction is only consequential for the C2-

Figure 4. P(θ)sin θ vs θ for the PWL (blue), MPWL (green), MOWL (purple), and OWL (red), where θ is the angle between the depicted water molecular axis and the z-axis. 2189


Article

ACS Nano

Figure 5. P(θ) sin θ vs θ for water involved in interlayer hydrogen bonds, where θ is the angle between the depicted water molecular axis and the z-axis. Water molecules in the PWL (blue), MPWL (green), MOWL (purple), and OWL (red) are displayed in pale colors if they act as hydrogen-bond donors and dark colors if they act as hydrogen-bond acceptors.

Figure 6. Layer lifetime correlation functions for water in the OWL (a), MOWL (b, purple), MPWL (b, green), and PWL (c).

found in the MPWL and a ∼45% chance it will be found in the MOWL due to the difference in layer populations. The secondary relaxations of C(t) correspond to characteristic times of ∼8−10 ns and can be attributed to water crossing through slits in P/O-graphene sheets. The average time a water molecule remained in a hydrolayerits layer lifetime, τlwas calculated via integrating S(t), where

slightly dulled for water in the hydrolayer opposite the slit, and there is convergence between the PWL and MPWL distributions as well as between those of the OWL and the MOWL. For the layer next to the slit, the water orientations become almost uniformly distributed. Hydrolayer Dynamics. In the equilibrium simulation, the water residence time in a hydrolayer was calculated using two binary time correlation functions, S(t) and C(t), analogous to those used by Thar et al. to describe hydrogen-bonding dynamics.22−24 Here, S(t) corresponds to the probability that a water molecule has continually remained in its original hydrolayer at all intermediate times between the initial time (i.e., t = 0) and the observation time t. C(t), on the other hand, only requires the water molecule to be in the same layer at the initial and the observation time. To ensure that water had truly left its initial layer, a buffer region (0.15 nm) was introduced between hydrolayers. If a water molecule in this region subsequently returned to its original hydrolayer, it was considered never to have left. However, if a water molecule in this region moved to the adjacent water layer, it was considered to have left from the time that it entered the buffer region. Both correlation functions are plotted in Figure 6. For each plot, C(t) shows two relaxation time scales. An initial relaxation occurs on the same time scale as S(t) (on the order of 10−1 ns) to 0.5 for the pristine and oxidized water layers, PWL and OWL, and ∼0.55 and ∼0.45 for the mixedoxidized and mixed-pristine water layers, MOWL and MPWL, respectively. These fast decays can be attributed to the exchange of water molecules with the adjacent hydrolayer. After this relaxation, there is a 50% chance for the PWL and the OWL that a water molecule will be found in its original layer and a 50% chance that it will be found in its adjacent layer. There is a ∼55% chance a water in the mixed region will be

⎡ t⎤ S(t ) = exp⎢ − ⎥ ⎣ τl ⎦

(2)

Values of τl are displayed in Table 4. Table 4. Layer Lifetimes in Different Hydrolayer Types layer type

PWL

MPWL

MOWL

OWL

τl (ps)

120

133

113

179

The difference in layer lifetimes between the PWL and the OWL is expected to be a result of interactions between water and hydroxyl groups in the OWL. At a glance, this rationale seems to be contradicted by the layer lifetimes of the MPWL and MOWL. However, because the system is in equilibrium, the rate at which water leaves a hydrolayer must be equal to the rate at which water enters the hydrolayer. As was discussed above, water exchange on this time scale occurs primarily with the adjacent hydrolayer. Therefore, the larger value of τl for the MPWL when compared to that of the MOWL is due to its greater density. In other words, there are more water molecules in the MPWL than in the MOWL, so it takes longer for them all to be exchanged, producing a greater average layer lifetime. Equilibrium Water Transport. The mean square displacement of water parallel to the graphene membrane (x−y MSD) 2190


Article

ACS Nano in the equilibrium simulation is plotted in Figure 7 for individual hydrolayers and for channels. For hydrolayers, the

Figure 8. Time correlation functions describing the time a water molecule remains in its original channel in the non-equilibrium simulations.

Figure 7. Water MSD parallel to the graphene membrane. Solid lines pertain to water in an individual hydrolayer, and dashed lines pertain to water in a channel.

These time constants correspond to average rates of water crossing between channels of 193 molecules ns−1 in the PPPP membrane, 146 molecules ns−1 in the POPO membrane, and 94 molecules ns−1 in the OOOO membrane. This trend is expected based on the relative rates of lateral diffusion found in the equilibrium simulations above (Figure 7).

MSD calculation was terminated when a water molecule left its original layer (see Hydrolayer Dynamics section for details), and for channels, the calculation was terminated when a water molecule left its original channel (i.e., crossed a slit). We examined the time dependence of the MSD using a power law 2

⟨|Δr(t )| ⟩ ∝ t

γ

CONCLUSIONS Channels composed of overlapping pristine graphene sheets or oxidized graphene sheets have been compared in numerous studies. Here, we argue that the heterogeneous structure of GO likely leads to the overlap of oxidized and unoxidized graphene, producing a “mixed” region type. In this region, the two hydrolayers show significantly different structural and dynamic features. In the mixed region, water density is greater in the MPWL than in the the MOWL, even though they are directly adjacent. The difference in numbers of intralayer hydrogen bonds is greater between the MPWL and the MOWL than it is between the PWL and the OWL. Similar disparities are found for water orientation in the mixed hydrolayers, which was correlated to interlayer hydrogen bonding. From our study, we conclude that this region may play a significant role in the rapid transport of water across GO membranes. In both equilibrium and non-equilibrium simulations, the dynamics in the mixed region lie between those of the oxidized and the pristine regions. For additional insight, we have performed equilibrium simulations in the presence of graphene sheets oxidized by epoxide groups. Their preliminary results indicate that the basic trend of water transport remains essentially the same as the case considered here (i.e., graphene oxidation by hydroxyl groups). While the full analysis of the epoxide case is postponed for a future study, this finding suggests that our conclusion on water transport in GO membranes is not limited just to the hydroxyl case.

(3)

In the diffusion regime, MSD increases linearly with time, that is, γ becomes unity ⟨|Δr(t )|2 ⟩ = 4Dt

(4)

where D is the two-dimensional translational diffusion coefficient. On the relevant time scales (i.e., t < τl), water in the oxidized and mixed-oxidized water layers (OWL and MOWL) did not reach the diffusive regime, but water in the pristine and mixed-pristine water layers (PWL and MPWL) did (respective diffusion coefficients were 1.5 and 0.92 nm2 ns−1). Diffusion coefficients obtained for the pristine, mixed, and oxidized channels were found to be 1.45 , 0.80, and 0.35 nm2 ns−1, respectively. As has been found in numerous previous works, water transport becomes hindered by its interactions with hydroxyl groups.1,7 Interestingly, however, the MPWL shows considerably faster translational dynamics than the MOWL probably because the MPWL interacts only with the hindered but fluid MOWL and does not interact directly with spatially fixed GO hydroxyl groups. This result indicates that mixed regions could play a prominent role in the lateral motion of water permeation across a GO membrane via the MPWL. This was particularly unexpected given the findings by Wei et al., where a sidepinning effect drastically reduced water flow in pristine capillaries bordered by oxidized regions.7 Non-equilibrium Water Transport. In order to explore the relative rates of water percolation in the presence of an auxiliary force, a time correlation functionsimilar to that used for the hydrolayer dynamics abovewas implemented for water in a channel. The results in the presence of a force field that yields an acceleration of 0.1 nm ps−2 in the direction perpendicular to GO sheets are shown in Figure 8. The correlation function decays as water passes across a slit. By integrating over the function, we find the average time a water molecule spends in a channel. The functions could all be fit with a single exponential, and the time constants for the pristine (PPPP), mixed (POPO), and oxidized (OOOO) membranes were found to be 2.5, 3.0, and 4.3 ns, respectively.

METHODS The GO membrane model used for this study consists of four 806carbon (3.2 nm × 6.4 nm) periodically repeating layers of P/Ographene with an intersheet distance of 0.9 nm and staggered, hydroxyl-functionalized 1.0 nm wide slits. Adjacent hydroxyl group pairsone on either face of the basal planeare randomly distributed across O-graphene sheets, oxidizing 20% of the basal carbons (156 hydroxyl groups). A similar degree of oxidation has been used in analogous studies.7−9 Alternate functional groups such as carboxylic acids along the edges and epoxides along the basal plane frequently occur at high oxidation states but are less common at lower oxidation states (e.g., 20%).25 Here, as well as in similar works, functionalization has been restricted to hydroxyl groups in order to maintain focus on the comparison of region types without variation in oxidation.7−9 2191


Article

ACS Nano §

The density of water in the membrane was established by equilibrating three P/O-graphene sheets (with slits) that are in contact with bulk water at a density of 1.00 g cm−3 and averaging the interlayer densities. Results were then used for equilibrium simulations in the PPOO membrane. The water densities in each region-type after the PPOO equilibration differ slightly from those established via the above protocol (see Table 5). The densities found in the PPOO equilibration are used for non-equilibirium simulations.

Permanent address for H.J.K.: Carnegie Mellon University.

ACKNOWLEDGMENTS This work was supported in part by the National Science Foundation through NSF Grant No. CHE-1223988. 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) Joshi, R. K.; Carbone, P.; Wang, F. C.; Kravets, V. G.; Su, Y.; Grigorieva, I. V.; Wu, H. a.; Geim, a. K.; Nair, R. R. Precise and Ultrafast Molecular Sieving through Graphene Oxide Membranes. Science 2014, 343, 752−754. (3) Sun, P.; Zhu, M.; Wang, K.; Zhong, M.; Wei, J.; Wu, D.; Xu, Z.; Zhu, H. Selective Ion Penetration of Graphene Oxide Membranes. ACS Nano 2013, 7, 428−437. (4) Shen, J.; Liu, G.; Huang, K.; Chu, Z.; Jin, W.; Xu, N. Subnanometer Two-Dimensional Graphene Oxide Channels for Ultrafast Gas Sieving. ACS Nano 2016, 10, 3398−3409. (5) Shen, J.; Zhang, M.; Liu, G.; Jin, W. Facile Tailoring of the TwoDimensional Graphene Oxide Channels for Gas Separation. RSC Adv. 2016, 6, 54281−54285. (6) Boukhvalov, D. W.; Katsnelson, M. I.; Son, Y.-W. Origin of Anomalous Water Permeation through Graphene Oxide Membrane. Nano Lett. 2013, 13, 3930−3935. (7) Wei, N.; Peng, X.; Xu, Z. Understanding Water Permeation in Graphene Oxide Membranes. ACS Appl. Mater. Interfaces 2014, 6, 5877−5883. (8) Wei, N.; Lv, C.; Xu, Z. Wetting of Graphene Oxide: a Molecular Dynamics Study. Langmuir 2014, 30, 3572−3578. (9) Wei, N.; Peng, X.; Xu, Z. Breakdown of Fast Water Transport in Graphene Oxides. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2014, 89, 012113. (10) Muscatello, J.; Jaeger, F.; Matar, O. K.; Müller, E. A. Optimizing Water Transport through Graphene-Based Membranes: Insights from Nonequilibrium Molecular Dynamics. ACS Appl. Mater. Interfaces 2016, 8, 12330−12336. (11) Sun, P.; Zheng, F.; Zhu, M.; Song, Z.; Wang, K.; Zhong, M.; Wu, D.; Little, R. B.; Xu, Z.; Zhu, H. Selective Trans-Membrane Transport of Alkali and Alkaline Earth Cations through Graphene Oxide Membranes Based on Cation-π Interactions. ACS Nano 2014, 8, 850−859. (12) 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. (13) Li, H.; Song, Z.; Zhang, X.; Huang, Y.; Li, S.; Mao, Y.; Ploehn, H. J.; Bao, Y.; Yu, M. Ultrathin, Molecular-Sieving Graphene Oxide Membranes for Selective Hydrogen Separation. Science 2013, 342, 95− 98. (14) Mi, B. Materials Science. Graphene oxide Membranes for Ionic and Molecular Sieving. Science 2014, 343, 740−742. (15) Williams, C. D.; Carbone, P. Selective Removal of Technetium from Water Using Graphene Oxide Membranes. Environ. Sci. Technol. 2016, 50, 3875−3881. (16) Erickson, K.; Erni, R.; Lee, Z.; Alem, N.; Gannett, W.; Zettl, A. Determination of the Local Chemical Structure of Graphene Oxide and Reduced Graphene Oxide. Adv. Mater. (Weinheim, Ger.) 2010, 22, 4467−4472. (17) Pacilé, D.; Meyer, J. C.; Fraile Rodríguez, a.; Papagno, M.; Gómez-Navarro, C.; Sundaram, R. S.; Burghard, M.; Kern, K.; Carbone, C.; Kaiser, U. Electronic Properties and Atomic Structure of Graphene Oxide Membranes. Carbon 2011, 49, 966−972. (18) Talyzin, A. V.; Hausmaninger, T.; You, S.; Szabó, T. The Structure of Graphene Oxide Membranes in Liquid Water, Ethanol and Water-Ethanol Mixtures. Nanoscale 2014, 6, 272−281.

Table 5. Interlayer Densities Are Measured in Units of Molecules per nm3 pristine

mixed

oxidized

23.0 22.7

20.5 20.6

19.0 19.0

bulk equilibration PPOO equilibration

All-atom optimized potentials for liquid simulation (OPLS-AA) parameters were used for the graphene sheets and hydroxyl groups as previously reported,26 though the hydroxyl charges and out-of-plane distortion of functionalized carbons were taken from the adapted model for GO by DeYoung et al.27,28 Test simulations with flexible graphene showed some variations in the interlayer distance, leading to the coexistence of water monolayers and bilayers. The majority of regions in the membrane allowed for the formation of water bilayers. In this initial study, we focused on the behavior of water in the bilayer regions. We therefore maintained the interlayer distance by freezing graphitic carbons for the duration of the simulation. The extended simple-point charge (SPC/E) model was used for water.29 OPLS-AA parameters in combination with the SPC/E water model have been reported to provide results in agreement with GO aggregation, interfacial behavior, and water contact angle data.9,26 All MD simulations were performed in the canonical (NVT) ensemble at 300 K with Nosé-Hoover temperature coupling using the Gromacs simulation package.30−33 Electrostatic interactions were computed using the particle mesh Ewald method with a 1.3 nm cutoff. Equilibrium simulations consisted of 10 ns of production time with a time step of 0.5 fs, preceded by 2 ns of prerun to ensure that the system had reached equilibrium. Equilibrium simulations were carried out using the PPOO membrane, wherein all environments (pristine, oxidized, and mixed) are present. Non-equilibrium simulations were performed in the PPPP, OOOO, and POPO membranes by applying an auxiliary acceleration perpendicular to the GO layers (along the z-axis) to all water molecules in order to promote passage across the membrane. Due to the time required for MD calculations, an acceleration comparable to those in ambient settings is impractical, thus, we applied a force that yields an acceleration of 0.1 nm ps−2. The intention of the data displayed in this section is to compare relative water flow rate in different membrane types and not as a direct comparison to experimental work. Non-equilibrium simulations were performed over 30 ns with 2 ns of prerun to ensure that the system has reached a steady state. All other non-equilibrium simulation details were equivalent to those for the equilibrium simulations discussed above. In addition to the simulations already listed, simulations with an alternative distribution of hydroxyl groups on the basal plane were run in equilibrium (PPOO-alt) and non-equilibrium (POPO-alt). Results (not shown) aligned quite well with those using the original hydroxyl distribution, indicating that our analysis is robust.

AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. ORCID

Jon A. L. Willcox: 0000-0002-3782-4770 Hyung J. Kim: 0000-0003-4334-1879 Notes

The authors declare no competing financial interest. 2192


Article

ACS Nano (19) He, K. T.; Wood, J. D.; Doidge, G. P.; Pop, E.; Lyding, J. W. Scanning Tunneling Microscopy Study and Nanomanipulation of Graphene-Coated Water on Mica. Nano Lett. 2012, 12, 2665−2672. (20) Algara-Siller, G.; Lehtinen, O.; Wang, F. C.; Nair, R. R.; Kaiser, U.; Wu, H. A.; Geim, A. K.; Grigorieva, I. V. Square Ice in Graphene Nanocapillaries. Nature 2015, 519, 443−445. (21) Bursulaya, B. D.; Kim, H. J. Molecular dynamics simulation study of water near critical conditions. I. Structure and solvation free energetics. J. Chem. Phys. 1999, 110, 9646−9655 and references therein. (22) Luzar, A.; Chandler, D. Hydrogen-Bond Kinetics in Liquid Water. Nature 1996, 379, 55−57. (23) Thar, J.; Brehm, M.; Seitsonen, A. P.; Kirchner, B. Unexpected Hydrogen Bond Dynamics in Imidazolium-Based Ionic Liquids. J. Phys. Chem. B 2009, 113, 15129−15132. (24) Rapaport, D. Hydrogen Bonds in Water. Mol. Phys. 1983, 50, 1151−1162. (25) Krishnamoorthy, K.; Veerapandian, M.; Yun, K.; Kim, S.-J. The Chemical and Structural Analysis of Graphene Oxide with Different Degrees of Oxidation. Carbon 2013, 53, 38−49. (26) Shih, C.-J.; Lin, S.; Sharma, R.; Strano, M. S.; Blankschtein, D. Understanding the pH-Dependent Behavior of Graphene Oxide Aqueous Solutions: a Comparative Experimental and Molecular Dynamics Simulation Study. Langmuir 2012, 28, 235−241. (27) DeYoung, A. D.; Park, S.-W.; Dhumal, N. R.; Shim, Y.; Jung, Y.; Kim, H. J. Graphene Oxide Supercapacitors: A Computer Simulation Study. J. Phys. Chem. C 2014, 118, 18472−18480. (28) Park, S.-W.; DeYoung, A. D.; Dhumal, N. R.; Shim, Y.; Kim, H. J.; Jung, Y. Computer Simulation Study of Graphene Oxide Supercapacitors: Charge Screening Mechanism. J. Phys. Chem. Lett. 2016, 7, 1180−1186. (29) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (30) Lindahl, E.; Hess, B.; van der Spoel, D. GROMACS 3.0: a package for molecular simulation and trajectory analysis. J. Mol. Model. 2001, 7, 306−317. (31) Berendsen, H.; van der Spoel, D.; van Drunen, R. GROMACS: A message-passing parallel molecular dynamics implementation. Comput. Phys. Commun. 1995, 91, 43−56. (32) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (33) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447.

2193


Molecular Dynamics Study of Water Flow across Multiple Layers of

Recommend Documents