As the solid-state nanochannels are hourglass-shaped as opposed to having a straight structure, little effort has been made on the transport phenomena across aquaporin-mimicking ultra-fine pores. Moreover, to the best of the authors’ knowledge, the existing studies on the water transportation through hourglass-shaped pore structures are limited to the symmetric geometries. In order to benefit from the advantage of combining the properties of nanoporous graphene and the unique structure of the aquaporin water channels, this chapter aims to study designing a highly effective permeable membrane for water transportation through nanopores with hourglass shape in multilayer structure using molecular dynamics (MD) simulation. Keeping the layer separation constant at 6 Å, we intend to examine the geometry effect on the pressure-induced transport of water molecules through hourglass-shaped nanopores. Three models of hourglass-shaped nanopores are selected: a symmetric structure α = 0 and asymmetric structure of α = 1/3 and α = −1/3. Our MD findings indicate that the permeability of water molecules highly increases across hourglass-shaped pores with an asymmetric structure of α = 1/3 because of the increased number of hydrogen bonds resulting from the length effect. Thus, we can conclude that an hourglass-shaped pore with a proper design can remarkably improve the rate of water transportation.
- water transport
- aquaporin water channel
- hourglass-shaped pore
- molecular dynamics simulation
The growing demand for water in everyday human activities makes producing fresh water a great challenge. Using membrane separation technology in the areas of water treatment is widely implemented. Producing portable water from existing resources such as surface water, seawater, brackish water, or groundwater can be done by using membrane processes. The membrane acts as a barrier with which two phases are separated, and transport of different chemical species are restricted selectively. It is basically referred to as a thin layer which separates materials depending on their chemical/physical properties when an external force (e.g., gradient pressure, electrical force, etc.) is applied across the membrane. Nanoporous graphene has gained great attention as promising membrane material for RO desalination because of its exceptional properties such as being mechanically strong, chemically robust, and extremely thin. Furthermore, using graphene—with such an extremely thin thickness—in desalination plants can remarkably yield higher water fluxes with significant reduction in the operating costs, for example, cost in terms of energy and cost in terms of capital and infrastructure requirements . Considering the graphene-based membrane as highly selective and permeable, the generation of pore with various pore shape and size can be achieved in graphene by using electron beam irradiation . Recent studies in both theory and experiment revealed that graphene [3, 4, 5, 6, 7, 8] and graphene oxide [9, 10] were highly selective and permeable and could be used as filtration membranes. Though single-layer graphene has demonstrated to be a promising candidate as a reverse osmosis membrane, synthesizing multilayer graphene is economically affordable [11, 12, 13, 14]. Considering parameters including pore size and chemistry, a wide range of studies have considered transport phenomena through single-layer graphene. However, less effort has been made on transport properties through nanoporous multilayer graphene. Experimentally, the studies on graphene oxide have revealed that constructing carbon-based material in molecular-scale porous structure could be achieved with tuned pore sizes and interlayer spacing [15, 16, 17, 18]. A better selectivity and higher permeability for artificial membranes could be achieved mimicking bio-inspired membranes such as aquaporins. Aquaporins are pore-forming proteins and ubiquitous in living cells, which provide fast water transportation in a selective manner resulting from their unique structure of hourglass shape . Recent study  showed that the desalination performance could be significantly improved incorporating the properties of aquaporin water channel into desalination membranes. Revealing an hourglass-shaped structure for the solid-state pores , transport phenomena pores in hourglass-shaped structure is not fully understood. Although the existing studies on water transportation properties across hourglass-shaped pore systems are limited to the symmetrical structures, to the best of our knowledge [22, 23, 24], further molecular-level studies on the transport phenomena across asymmetrical geometries are still needed [25, 26, 27]. The water transportation phenomena in confined environments (e.g., nanochannels), where the critical dimension can be a few molecular diameters, differ from macroscopic scales which explains the transport phenomena based on continuum theories (e.g., Navier-Stokes equations). Accordingly, a much-needed molecular-level understanding of the water transportation phenomena at the nanoscale is necessary. Molecular dynamics (MD) has become a powerful tool for investigating the effects of confinement on fluid properties and calculating phenomenological transport coefficients for use in multi-scale fluid flow models. A significant effort has been devoted, often with the help of MD simulation, to understand the water transportation properties in carbon nanotubes (CNTs) and graphene [28, 29, 30, 31]. Several researchers have shown a freezing state for water confined in CNTs at low temperature and the appearance of ice phases stabilizing by confinement [32, 33]. Using MD simulations, Koga and colleagues  showed that water molecules inside CNTs with diameters varying from 11 to 14 Å indicate phase transitions with applied pressure of 50–550 MPa. Experimentally, Kolesnikov and colleagues  observed a combination of water-chain and ice-shell structure for water molecules confined in CNTs having a liquid phase state even at low temperature. The appearance of single-layer and bi-layer ice structures were also disclosed for water molecules confining in between graphene slits under applied pressure [34, 35]. Another study revealed that water confined in graphene slits freezes at ambient temperature while no pressure is applied . Chakraborty and colleagues  have shown a diversity of results on the structure, dynamics, and thermodynamics of confined water molecules in CNTs and graphene oxide sheets. They found the appearance of the single-file ordered structure for water molecules confined in narrow CNTs and the appearance of the n-gonal ice structure in larger CNTs. Their findings also revealed the appearance of two-dimensional ordered structures for water molecules confined in slit-pore geometries. Considering the graphene properties combining with aquaporin water channel structure, this study aims to design an optimal and effective permeable membrane for water transportation across hourglass-shaped pores in multilayer nanoporous graphene using MD technique.
2. Simulation setup
Molecular dynamics (MD) is used to investigate the transport of water molecules across hourglass-shaped nanopores in nanoporous multilayer graphene. The schematic of the simulation models are given in Figure 1. Three pore models of hourglass shape of asymmetric (α = 1/3 and α = −1/3) and symmetric (α = 0) are chosen. Seven layers of graphene (31.9 × 34.3 Å2) with a constant interlayer separation of d = 6 Å with total length of 36 Å is represented as the multilayer structure. Considering the single-file water arrangement across the narrowest part and a bulk structure of water molecules in the largest segments, we considered being the narrowest (a) and the largest diameters (D) of hourglass-shaped structure in 4 and 7 Å, respectively, embedding in a periodic box of over 2000 SPC/E water molecules. To generate the pore in graphene layers, we removed adjacent carbon atoms in the center and considered the pores with no sharp edges (for more information, see Ref. ). To generate water flow in the models, we applied the pressure difference along –z direction across α = 1/3, α = −1/3, and α = 0 systems. A pressure-induced flow, ΔP = nf/A (denoting n and A as the number of water molecules and the membrane surface area, respectively), is achieved by applying a constant force, f, on water molecules inside the reservoir . All the simulations are performed using large-scale atomic/molecular massively parallel simulator (LAMMPS) simulation package , implementing the constant number (N), volume (V), and energy (E) (NVE) ensemble with dissipative particle dynamic (DPD) thermostat  aiming at maintaining the temperature constant at 300 K. We have chosen a time step of 0.5 fs in order for ensuring that water molecules between the bulk region and layers are accurately simulated. Considering a periodic boundary condition in three directions, a cutoff radius of 10 Å is selected for the Lennard-Jones (LJ) potentials. To reduce the computational cost, we have considered the carbon atoms of the graphene to be frozen. We have modeled only the interaction between water molecules (oxygen) and the graphene (C) by using the Lennard-Jones parameters as σO─O = 3.169 Å, ƐO─O = 0.1515 kcal/mol, σO─C = 3.190 Å, and ƐO─C = 0.0956 kcal/mol. Using the Ewald summation method  for the coulomb interactions, the SHAKE algorithm is applied to hold the water molecules rigid . Conducting more than 15 ns, MD runs for every three cases (α = 1/3, α = −1/3, and α = 0), the last 12 ns trajectory was considered for data analysis.
3. Results and discussion
In our simulations, we applied the pressure difference, ΔP, ranging from 365 to 635 MPa across the multilayer nanopore systems. Three pore models of α = 1/3, α = −1/3, and α = 0 with the long, short, and equal entrance lengths are considered. Figure 2a shows the average number of water occupancy with the applied pressures inside each of the three models. The water occupancy is calculated by counting the number of water molecules, which could diffuse through the 36-Å-long multilayer systems during the simulations. It has been found that the average water occupancy in the case of α = 1/3 is higher than that of α = −1/3, and α = 0 systems. Since the pore volume increases, at α = 1/3, for example, the influx of water molecules could be encouraged stemming from an increased cone length. Though a significant variation of the number of water occupancy in relevance to applied pressure is found in α = −1/3 model, generally, the value of water occupancy inside the pores is more affected by the pore structure compared to the applied pressure (e.g., α = 1/3 > α = 0 > α = −1/3). The number of water flux is calculated by counting the number of water molecules that cross the pore thoroughly, as shown in Figure 2b. Qualitatively, the results indicate that the water flux increases linearly with applied pressure, which is consistent with the recent study . Moreover, the MD results show that the number of water flux is higher in the case of α = 1/3 compared to α = 0 and α = −1/3 systems. Interestingly, the symmetric hourglass-shaped pore (α = 0) shows a lower flux of water molecule compared to the asymmetric pore structure (e.g., α = 1/3), which attributes to the entrance effects. Therefore, it is understood that the entrance effects have a significant contribution to flux variation inside pores with the structure of hourglass shapes. The permeability of water molecules through pores in hourglass shape is significantly restricted due to the energy barriers at the entrances. In case of α = 1/3 for example, a large amount of water molecules entering the long cone can contribute to reduce the energy barriers caused by the narrowest section of pore (of which accommodates water molecules only in single-file configuration), resulting from the increase in the number of hydrogen bonds. It is also inferred that the decrease in the cone angle would enhance the permeability of water molecules across hourglass-shaped pores in a multilayer structure (indeed, this implication very recently discussed for the permeability of water molecules in solid-state nanochannels of hourglass shape ). To assess the transport capacity of a membrane, one could calculate an important parameter named osmotic permeability. By drawing a best-fit slope through flux versus applied pressures, the osmotic permeability of a membrane can be achieved as shown in Figure 2b. In agreement with the results reported for flux, the corresponding magnitudes are (4 ± 1) × 10−14, (7.52 ± 9.76) × 10−15, and (2.337 ± 4.33) × 10−15 cm3 Pa−1 s−1 for the pore systems of α = 1/3, α = 0, and α = −1/3, respectively. It is also inferred that achieving high permeability of water molecules across the pore of hourglass shape would strongly depend on the structure of the pore. Therefore, a design of hourglass-shaped pore with which adsorbs bulk water molecules can be a promising geometry aimed at obtaining high permeability. Within this context, the water molecules in bulk-like manner contribute to remarkably reduce the energy barriers inside the pores because of the size effects. Next, we investigate the viscosity dissipation inside asymmetric pore systems at the smallest section of the pores. Very recently, a numerical study  has shown that the entrance effects (defined as viscous dissipation) would dramatically restrict the water permeability inside solid-state pores of hourglass shape. Basically, the molecules of water confined in graphene nanopores show different viscosity compared to the carbon nanotubes (CNTs). Considering a similar characteristic of pore size, a recent study has shown that the molecules of water confining in a graphene pore pose higher viscosity compared to the water molecules confining in CNTs. Increasing the pore diameter leads to the increase in the viscosity of confined water molecules in CNTs approaching the bulk viscosity [29, 43]. In graphene nanopores, the viscosity of water molecules increases, while decreasing the pore diameter in accordance with the structure of water molecules . The structure of water molecules along the flow direction inside graphene nanopores shows different behavior compared to CNTs’ nanochannels. Recent MD study  has shown that the viscosity of water molecules across graphene nanopore is linearly increased in accordance to the layering effect. In fact, the viscosity of water molecules through graphene nanopores is in relation to the layering effect of water molecules. Through the density distribution, one can obtain the water layering effect along the flow direction in the graphene system. Here, we measured the water layering factor as (ρmax − ρmin)/ρbulk, to quantify the water layering effect inside only the asymmetric pore structures, noting that ρmax and ρmin are the maximum and minimum density numbers in the vicinity of the layer with the narrowest pore diameter, respectively, as shown in Figure 3. The ρbulk is estimated to be in the range of bulk fluid densities reported, 0.8σO─O−3 ≤ ρbulk ≤ 1.1σO─O−3 . When α = −1/3, the density profile near the layer with the smallest pore diameter shows greater distribution compared to α = 1/3, suggesting a significant effect of water layering inside α = −1/3. Within this context, the water layering effect increases the viscosity of water molecules while crossing the smallest diameter in single-file arrangement. Accordingly, a significant reduction in the local permeability inside α = −1/3 is expected.
Furthermore, in case of α = 1/3, large amount of water molecules (bulk-like) approaching the narrowest diameter of the pore should contribute to diminish highly the water viscosity of single-file configuration because of the increase in hydrogen bonds. Although the dipole orientation of water molecules inside CNTs maintains oversimulation time, in graphene nanopores, the dipole orientation of water molecules is frequently flipped . Here, we intend to show how dipole orientation of water molecules may behave across pores in multilayer system. Taking into account all the water molecules inside pores, the averaged dipole orientation was estimated by computing the angle between water dipole and the axial direction (perpendicular to the surface of the multilayer). Figure 4a and b shows the profiles of the averaged dipole orientation of water molecules inside asymmetric pores α = 1/3 and α = −1/3, respectively. Over 4000 calculated values are used to plot the average dipole orientation of water molecules inside asymmetric structures (α = 1/3 and α = −1/3). In case of α = 1/3, the average dipole orientation fell into two range groups of 149–153° and 158–162°, while in the case of α = −1/3, the value of average dipole orientation fell only in the range of 44–47°. This implies a degree of water molecule ordering inside the pore in the case of α = 1/3, which enhances the formation of hydrogen bonds between water molecules. Moreover, the variation of the dipole orientation tends to strengthen the water-water interaction resulting in higher hydrogen bonds. This is responsible for fast transport of water molecules in monolayer graphene . We believe that this implication could be also valid in multilayer graphene systems.
Although, the average dipole orientation inside α = 1/3 and α = −1/3 systems is only flipped during the course of few picoseconds, it maintains for whole the simulation time showing a higher magnitude of dipole inside α = 1/3, as given in Figure 4c. At constant layer spacing (d = 6 Å), the obtained results show that the dipole orientation of water molecules inside the multilayer graphene membrane is different from that across single-layer graphene, in which the dipole orientation has frequently flipped. This implies that every single layer in multilayer system provides the energy barrier which allows the water molecules to freely flip. It would be expected that by increasing the layer spacing up to several angstroms, this phenomenon changes. The fluctuation of free energy occupancy in relation to the number of water molecules across asymmetric pores is given in Figure 5a and b. We computed the free energy using the equation, G(N) = −Ln[P(N)]/b, where P(N) is the probability of finding the exact N molecules of water in the system, b = (kBT)−1, KB and T stand for the Boltzmann constant and temperature, respectively. Ranging from 896 to 926 water molecules, the most probable of water molecules is being with 904 molecules inside α = 1/3. However, considering 722–756 water molecules inside α = −1/3, the 728 is the most probable. In comparison to α = −1/3, the water molecules are increasingly accommodated inside α = 1/3 because of the strong interplay between molecules, which results in forming the stabilized chain of water molecules. In comparison to α = −1/3, a wider distribution of probability with higher number of water molecules are found inside α = 1/3. As a result, a high rate of water conductivity can be found inside fully filled asymmetric model of α = 1/3. In addition, the probability of existence of higher numbers of water molecules is scarcely observed inside α = −1/3 due to the strong water-wall interplay.
The free energy of occupancy in relation to the hydrogen bonds number is given in Figure 5c. The hydrogen bond was calculated in accordance to the method reported in Ref. . As the figure shows, the average number of hydrogen bonds varies from 1.48 to 1.61 inside α = −1/3, whereas it ranges from 1.56 to 1.73 for the case of α = 1/3. Then, a wider distribution for probability of formation of higher numbers of hydrogen bonds inside α = 1/3 suggests an enhancement in the water density distribution. In other word, the establishment of the stable chains of water molecules inside α = 1/3 can support formation of the long-lasting network of hydrogen bonds which is significantly contributed to achieving higher flux. Radial distribution function (RDF) is an important parameter to analyze the structure of molecules in confined environments. Figure 6a depicts the profile of the RDFs between the water molecules (oxygen atoms) inside asymmetric pore models. A rapid decay for the RDFs can be seen inside both pore models of α = 1/3 and α = −1/3. According to the high magnitudes of the RDFs at low interatomic distances, the formation of solid clustering between water molecules inside α = 1/3 is predicted, which triggers potential reduction in the energy barrier caused by the vacant spaces between molecules, as shown in Figure 6b. Furthermore, these results of the RDF are in agreement with the results of the free energy fluctuations. The transport mechanism of molecules in confined spaces could be monitored by using the potential mean force (PMF). Using the RDF which illustrates how density varies as a function of distance from a reference particle, one can calculate the PMF since the density would be a key index to explain probability . Figure 6b illustrates the profiles of the potential mean forces (PMFs) for water molecules inside α = 1/3 and α = −1/3 pore systems. Because there is a weak interaction between the molecules of water inside α = −1/3 pore model, a large energy barrier caused by vacant spaces between the molecules of water is understood. We believe that the large energy barrier could be attributed to the increased viscosity caused by the pore size at the smallest section, which dramatically reduces the passage of molecules of water. Encouraging the water molecules in bulk structure to get through the pore of α = 1/3, they can significantly contribute overcoming the energy barrier. Therefore, water molecules with a tight hydrogen bonding network can pass the smallest section of the pore, which contributes strongly improving the entire passage of molecules of water across the pore. To further understand the mechanism of water transportation in confined geometries, an important parameter known as diffusion coefficient can be used. We computed the axial diffusion coefficient by calculating the mean square displacement (MSD), as given in Ref. . The diffusion rate is dependent on density because it stems from the collision between molecules and the interplay between molecules and the pore surface . The axial diffusion coefficient of water molecules in the −z coordinate of the symmetric and asymmetric pore systems is given in Figure 7a. The figure shows that the axial diffusion coefficient decreases with applied pressure inside α = 1/3; however, it is increased inside α = −1/3 and α = 0 pore systems indicating a larger slope for diffusion inside α = −1/3. It can be found that the increase in the density with applied pressure could be the responsible for the increase in the diffusion coefficient of water molecules inside α = −1/3 and α = 0 pore models. Looking at the low density states (in relation to low applied pressures), we believe that the water layering would be extremely small and the water molecules confined in between the layers would provide very limited spaces for molecules to move due to the strong interaction between water layer surfaces. Considering states with higher densities, the large number of hydrogen bonds would contribute to form the larger clusters yielding larger movement. In case of α = 1/3 pore, by increasing the applied pressure, the diffusion coefficient decreases. We believe that the water density does not have a significant effect on the mechanism of diffusion across α = 1/3 pore, as can be realized through the result reported for the water occupancy (see Figure 2a). The decreased diffusion coefficient with applied pressure could be ascribed to the formation of stabilized hydrogen bonding networks resulting from weak interaction between water layer surfaces. In comparison with α = 0 and α = 1/3 pore systems, a sharp slope for diffusion coefficient is observed inside α = −1/3 which suggests a density-dependent effect. Since the applied pressure increases, a collective diffusion mechanism can be understood inside α = 1/3 because of the formation of the stabilized hydrogen bonding networks. A discreet mechanism for filling behavior of water molecules with applied pressure can be realized across α = −1/3 and α = 0 systems due to existence of the large energy barrier. At the overlap pressure (444 MPa), the similar magnitudes for the diffusion coefficient of water molecules inside each of the three pore systems are obtained. The average residence time of water molecules inside α = 1/3 system as a function of the pressure difference is given in Figure 7b. Considering a molecule entering from one side and leaving the pore from the other side, one can calculate the residence time tres = tout − tin. In Figure 7b, we can see that the residence time inside α = 1/3 shows a negligible variation when the applied pressure increases from 365 to 444 MPa, whereas it highly decreases while the applied pressure increases from 444 to 635 MPa which suggests the complete passage of a molecule across the pore in a short course of simulation time. Accordingly, it is inferred that the time residence is dependent on the applied pressure. Since the applied pressure increases, water molecules in bulk structure are encouraged to the cone and exit the pore in a fast behavior. This suggests that the structure of the pore has less contribution to the movement of water molecules. The profiles of the axial and total velocities of water molecules inside pores are shown in Figure 8a and b as the function of applied pressure. By increasing the applied pressure, the axial velocity of water molecules is strongly increased inside α = 1/3, whereas a slight increase is found for the axial velocity inside α = 0 system. In case of α = −1/3, the increase in the applied pressure results in the reduction in the axial velocity inside the pore. It can be found that by increasing the applied pressure, the occupancy number of water molecules inside α = −1/3 is increased. Due to the strong interaction between water layer surfaces it should be difficult for the large displacement of water molecules across the pore model of α = −1/3, which results in producing low axial velocity. In addition, the figure describes high speed of transport of water molecules with applied pressure across α = 1/3 and α = 0 systems because of the weak interaction of water layer surfaces and the strong water–water interaction. Accordingly, the formation of tight hydrogen bonding networks is being expected inside α = 1/3 and α = 0 pores. It is realized that the pore structure plays a dramatic role on the transport phenomena across pores in hourglass structure. Figure 8b indicates the profile of the total velocities of water molecules in relation to applied pressure inside each of the three pore models. Generally, the total velocity (Vx + Vy + Vz) is higher inside α = −1/3 compared to α = 1/3 and α = 0 systems. At applied pressure of 365 MPa, our simulation results show the same value for the total velocity inside symmetric and asymmetric pore systems. It is understood that at given applied pressure, the velocity in the x and y directions has a significant impact on total velocity inside α = −1/3 model compared to α = 1/3 and α = 0 structures. The results reported in Figure 8a and b together suggest that the contribution of velocity in the z direction (Vz) should indicate the permeation ability in our simulation models, though a negligible effect of the velocities in x and y directions on the transport mechanism could be understood.
4. Concluding remarks
Graphene with one-atom thickness of sp2 carbon atoms organizing a honeycomb structure possess an excellent mechanical strength, supreme electrical and thermal conductivity, and other amazing properties. Lately, the application of nanoporous graphene in desalination and water purification processes has been greatly increased. Using experimental techniques [51, 52, 53, 54], one can generate nano-/micro-scale pores of various diameters and shapes on graphene sheets. Considering the properties of nanoporous graphene and taking advantage of aquaporin water channel structure, this computer simulation study aims to design an optimal and effective permeable membrane for water transportation across hourglass-shaped pores in multilayer nanoporous graphene.
In this study, by conducting MD simulations, we investigated the mechanism of water transportation across hourglass-shaped nanopores in symmetric and asymmetric structures constructed by multilayer nanoporous graphene. Considering a constant interlayer spacing of 6 Å, we particularly focused on the applied pressure effects on the transport properties through symmetric (α = 0) and asymmetric (α = 1/3 and α = −1/3) pores with hourglass shape in multilayer graphene. The results show a higher occupancy of water molecules in relation to applied pressure inside α = 1/3 compared to α = −0 and α = −1/3. In case of α = 1/3, because of the increase in the volume of the pore (the increase in the cone length), the flow of molecules of water is fostered. Then, the occupancy of water molecules inside the pore with hourglass shape is geometry dependent compared to pressure effect. The number of flux of water increases linearly as a function of applied pressure inside the pores, suggesting higher flux inside α = 1/3 pore. In comparison with α = −1/3, the water molecules crossing the smallest pore diameter inside α = −1/3 system have shown a larger layering effect which strongly reduce the water transport rate. The findings revealed that the asymmetric model of α = 1/3 is able to transfer a larger number of fluxes in comparison to the symmetric pore of hourglass shape (α = 0). This could be ascribed to the small viscosity in relation to small water layering which comes from the length effects. It is also understood a proportionality between the permeability and the pore structure with a higher value for the case of α = 1/3. In comparison with α = −1/3, the MD trajectory indicated that the dipole orientation is widely distributed inside α = 1/3, falling in the ranges of 149–153° and 158–162°. Basically, a large distribution of dipole orientation triggers the fast flow of water inside the pore. It can be seen that the diffusion coefficient decreases with the pressure difference inside α = 1/3 pore because of the strong clustering of water molecules. A collective and discreet diffusion mode is being inferred for transport mechanisms inside α = 1/3 and α = −1/3 systems, respectively. Finally, our MD results described that the contribution of the velocity perpendicular to the layers (z-axis) exposes the influx ability of molecules of water inside the pores.
This work was supported by the National Research Foundation of Korea (NRF) Grant funded by the Korean government (MEST) (no. NRF-2014R1A2A2A01003618) and the Pioneer Research Center Program (Grant no. NRF-2012-0009578).