Gabriele
Tocci
*a,
Maria
Bilichenko
a,
Laurent
Joly
bc and
Marcella
Iannuzzi
a
aDepartment of Chemistry, Universität Zürich, 8057 Zürich, Switzerland. E-mail: gabriele.tocci@chem.uzh.ch
bUniv Lyon, Univ Claude Bernard Lyon 1, CNRS, Institut Lumière Matière, F-69622, Villeurbanne, France
cInstitut Universitaire de France (IUF), France
First published on 5th May 2020
Despite relevance to water purification and renewable energy conversion membranes, the molecular mechanisms underlying water slip are poorly understood. We disentangle the static and dynamical origin of water slippage on graphene, hBN and MoS2 by means of large-scale ab initio molecular dynamics. Accounting for the role of the electronic structure of the interface is essential to determine that water slips five and eleven times faster on graphene compared to hBN and to MoS2, respectively. Intricate changes in the water energy landscape as well as in the density correlations of the fluid provide, respectively, the main static and dynamical origin of water slippage. Surprisingly, the timescales of the density correlations are the same on graphene and hBN, whereas they are longer on MoS2 and yield a 100% slowdown in the flow of water on this material. Our results pave the way for an in silico first principles design of materials with enhanced water slip, through the modification of properties connected not only to the structure, but also to the dynamics of the interface.
Here, we tackle this challenge in full strength by using AIMD to investigate the slippage of water confined between two-dimensional materials. We concentrate our efforts on connecting the atomistic features and electronic structure of the interface to the friction of water on graphene, hBN and MoS2 sheets, motivated by a number of experimental works on nanofluidic transport at such interfaces.8,23–25 We find that the friction on graphene is about five and eleven times smaller than on hBN and MoS2, respectively. The structure of water in the direction perpendicular to the sheets is very similar across the three materials. Yet, slight changes in the corrugation of the water energy landscape of at most a few kBT‡ in the direction parallel to the sheets impact the friction dramatically. Interestingly, slippage is also controlled by the density correlations of the contact layer, and their time relaxation provides the main dynamical origin of the friction. These results illustrate how a delicate interplay of interfacial forces, which are governed by the electronic structure of materials, can have a tremendous impact on nanofluidic transport. Additionally, favourable comparison with experimental slip lengths on carbon nanotubes with small curvature and with BN nanotubes illustrates the relevance of the electronic structure on nanofluidic transport.
We performed a series of extensive AIMD simulations of water confined between graphene, hBN and MoS2 sheets for different confinement widths, as well as with water films supported on the three types of sheets for systems containing up to 2000 atoms. AIMD simulations were carried out in the Born–Oppenheimer approximation with the electronic structure computed at the DFT level using the optB88-vdW exchange–correlation functional26,27 and the dynamics was propagated in the canonical ensemble for about 100 ps at 330 K. Snapshots of the simulated systems for MoS2 with different confinements are shown in Fig. 1.
Fig. 1 Water density profiles for different confinement heights. Black, red and green graphs are density profiles for graphene (GRA), hBN and MoS2, respectively. The snapshots in (a), (e), (i) and (m) illustrate the configurations of water confined between MoS2 sheets for different confinements. The density profiles for the supported water films are shown in the ESI.† |
We start our discussion on the structural properties of water under confinement. A key property of confined fluids is their density profile as a function of the distance from the sheets,28 which is shown in Fig. 1. The water density profiles were computed with respect to the instantaneous height between each fluid molecule and the nearest sheet's atom, in the spirit of ref. 29 (ESI†). There is a strong similarity in the layering of water between the three materials for each confinement height. The properties of the first water layer, so-called contact layer, are most relevant to understand nanofluidic water slip.30 Therefore, particular emphasis is placed on the understanding of its structure. Except for the smallest confinement width, the contact layer density peaks at approximately 3.5 g cm−3, regardless of the type of material. For confinement heights between 1.2 and 1.5 nm (Fig. 1f–h) two additional smaller peaks appear beyond the contact layer. Gradually less pronounced density oscillations are also visible beyond the first two layers for confinement widths between 2.0 and 3.0 nm (Fig. 1j–p). A bulk liquid region is almost recovered in systems whose confinement width reaches about 3.0 nm (Fig. 1n–p). The differences observed in the first density peaks seen in Fig. 1(b)–(d) and in the second peaks in Fig. 1(f)–(h) result from commensurability between the size of the water molecules and the confinement width of each system,31 along with differences in the magnitude of the water/sheet interactions.§
We now turn to water slip, as the most relevant contribution of this work. At liquid–solid interfaces where the major contribution to the water/surface interaction comes from dispersion forces – such as those considered here – the tangential velocity of the fluid in contact with the solid, so-called slip velocity vslip, is non-zero and the friction force parallel to the flow is proportional to the slip velocity via the so-called friction coefficient λ, following the relation ,32,33 with the surface area. Whereas the relation of proportionality between the total force acting on the liquid and the slip velocity can be used to calculate λ from nonequilibrium molecular dynamics simulations,34 here we apply the following Green–Kubo relation and compute the friction coefficient from the equilibrium fluctuations of the force acting on the liquid extracted from our equilibrium AIMD trajectories:33,35
(1) |
It can be seen that water flow is fastest on graphene, as λ on graphene is about 5 times smaller compared to hBN and about 11 times smaller than MoS2. Additionally, there is no apparent effect of confinement on friction, as λ is approximately constant as a function of the fluid thickness for all materials.
To provide mechanistic insights into the observed trends of the friction coefficient, eqn (1) is formulated in terms of the second moment of the force 〈F2〉 and the force correlation time τF:
(2) |
(3) |
From eqn (2) it is possible to disentangle the dependence of the friction coefficient on static and dynamical properties of the interface, where 〈F2〉 and τF entail a static and dynamical contribution, respectively. We first illustrate the dependence of λ on 〈F2〉 as shown in Fig. 3 for the systems considered. A clear correlation between the friction and the mean squared force is observed: upon going from graphene to hBN to MoS2 increases, and so too does λ, for all confinement widths. Based on dimensional grounds, the dependence of λ on can be further expressed as , where 〈ΔV2〉 is a measure of the corrugation of the potential energy of the fluid in the contact layer and a0 is the lattice constant of the sheets.11,30 The corrugation of the water energy landscape can be estimated from the calculation of the free energy profile defined as ΔG(x,y) = −kBTlog[PO(x,y)/PO(xmax,ymax)], where PO(x,y) is the two-dimensional probability distribution of an oxygen atom in a water molecule in the contact layer and the subscript “max” indicates the point where the probability is at the maximum. The free energy profile thus computed is shown in Fig. 3(b)–(d).
Fig. 3 Friction coefficient λ as a function of 〈F2〉/A (a), and free energy landscape of water in the contact layer for (b) graphene, (c) hBN and (d) MoS2, where snapshots of the sheets are superimposed on part of the graphs. In (a) 〈F2〉/A is normalized by its value on graphene at a confinement width of about 0.9 nm, and the different symbols refer to systems with different fluid thicknesses, as per Fig. 2. |
The level of corrugation can be extracted from the free energy barrier, obtained as the maximum in ΔG, which is of about 0.01 eV, 0.02 eV and 0.05 eV for graphene, hBN and MoS2, respectively. Whereas the free energy barriers in graphene and hBN are smaller than the thermal energy at room temperature (about 0.025 eV), the free energy barrier on MoS2 is about 0.05 eV, and therefore larger than kBT, indicating definitely a more corrugated overlayer. The increased level in the corrugation of the free energy explains the observed relation between λ and and it provides a static and structural origin for the observed behaviour of the friction coefficient. Looking closely at Fig. 3, some changes in the dynamics can be anticipated by noting that, while λ and 〈F2〉 appear to be proportional for graphene and hBN, the friction of water on MoS2 is too large or it to be explained simply based on its dependence on 〈F2〉.
Indeed, we provide evidence of an intriguing mechanism underlying friction, inherent to the dynamics of the water contact layer. To this end we plot λ as a function of the force correlation time τF for the systems discussed so far, see Fig. 4. The most apparent feature is that τF is approximately 50 fs for graphene and hBN, considering all fluid thicknesses, whereas it is about 100 fs for MoS2. The larger value of τF indicates a slower relaxation of the lateral force acting on the fluid at the water–MoS2 interface compared to the other two interfaces. The dependence of friction on τF can be tuned by changing the amplitude of the bond length oscillations of the sheets,15 and by modifying the magnitude of the electrostatic and dispersion interactions (see ESI†). Because depositing two-dimensional materials on different substrates may lead to a change in the electrostatic and/or van der Waals interactions at the interface, it is possible to test experimentally this dependence.25,37,38
Fig. 4 Dynamical origin of water slip on two-dimensional materials, as shown by the dependence (a) of the friction coefficient λ on the force correlation time τF, and (b) of the density correlation time τρ on τF. The same symbols are used as in Fig. 2. |
In order to provide further insights into the molecular origin of this phenomenon, eqn (1) is approximated in terms of the structure and dynamics of the liquid contact layer (ESI†):30,39
(4) |
(5) |
Altogether, eqn (2) and (4) imply a correlation between the force correlation time and the density correlation time, i.e. τF ∼ τρ. As it can be seen in Fig. 4(b), both τF and τρ are approximately twice as small in graphene and hBN compared to MoS2. The apparent correlation between τF and τρ provides the main explanation for the different dynamical contribution to the friction. The density correlation time on graphene and hBN is smaller than on MoS2 because (i) the density modes on the former materials are associated with a smaller lattice constant a0 and therefore decay more rapidly40 and (ii) the collective diffusion of the fluid probed at a specific wave-vector (Dq∥) is larger on graphene and hBN compared to MoS2. This can be seen by expressing τρ in terms of the wave-vector and of Dq∥ as τρ = 1/(Dq∥q2).30 By taking τρ for graphene (or hBN) and MoS2 to be about 400 fs and 1000 fs, respectively, we can extract Dq∥ to be about 0.38 Å2 ps−1 for graphene and hBN and 0.25 Å2 ps−1 for MoS2. As a result, we find that the differences in the lattice constant (i.e. the q−2 term), and in the collective diffusion coefficient equally contribute to the observed behaviour of τρ. In the ESI (see Fig. S13†) we also show that stretching and compressing the underlying lattice has a strong effect on τρ, but the larger MoS2 lattice constant is not sufficient to explain alone the larger value of τρ on MoS2. Instead, there is a remaining contribution to τρ on MoS2 beyond the obvious q−2 dependence, which arises from the slower collective diffusion of water on MoS2 compared to the other two materials, as it can be seen by compressing MoS2 to approximately the same lattice constant of hBN and graphene. Despite the evident correlation between τF and τρ, the latter is smaller by almost a factor of 10. A possible reason for the discrepancy is that the density correlation time results from an average over the density modes for the whole contact layer, whereas the force correlation time is more sensitive to those modes closer to the contact layer.30 Note also that the derivation of eqn (4) is based on truncating the Fourier expansion of the potential energy landscape of the wall, but higher order modes could contribute.
Our results can be used in practice to test experimentally a possible dependence of the viscosity on confinement, which remains still an open question.31,41–44 Indeed, liquid–solid slip is most often characterized experimentally through the slip length b, defined as the ratio between the slip velocity vslip and the shear rate at the interface: b = vslip/. Because at the interface the viscous shear stress η (with η the shear viscosity) is equal to the wall shear stress λvslip, the slip length can be expressed as the ratio between the viscosity and the friction coefficient:45b = η/λ. Using the experimental viscosity of bulk water (η = 10−3 Pa s), we obtain on average a slip length of 19.6 ± 1.9 nm, 4.0 ± 0.4 nm and 1.8 ± 0.2 nm for graphene, hBN and MoS2, respectively. As we report little to no effect of friction on confinement in all three materials, any dramatic change in the slip length measured in experiments can be ascribed to a change in the viscosity, via the expression b = η/λ. The surface force apparatus, atomic force microscopy as well as scanning tunneling microscopy are some of the most promising techniques to investigate changes in the slip length and in the viscosity under confinement.25,46
Additionally, measurements of the slip length in nanotubes using optical microscopy have shown that water essentially does not slip inside BN nanotubes, whereas it does slip inside carbon nanotubes, the slip length being larger in nanotubes with smaller radii.8 In the carbon nanotubes with the lowest curvature (100 nm in diameter), the reported slip length of about 23 nm compares favorably with our predicted value for graphene, 19.6 ± 1.9 nm. For BN nanotubes, no slip is detected above the detection limit of 5 nm, consistent with our prediction of 4.0 ± 0.4 nm. This comparison indicates that the difference in the slip length between carbon and BN nanomaterials may originate to a large extent from the different underlying corrugation of the energy landscapes of these materials. Nevertheless, the proclivity for adsorption of hydroxide ions may provide an additional reason for the smaller slip length on BN nanotubes.8,22,47–49 Carrying out experiments at different levels of pH would likely provide a definitive answer to this question.
It has been pointed out before that there is a large uncertainty over the simulated slip length of water on graphene or graphite, with values between 1 and 80 nm.50 In this work we have contributed to reduce this uncertainty, also given the good agreement with recent experiments for large carbon nanotubes, stressing the relevance of ab initio methods for studies of liquid/solid friction. Instead, using previously reported water/sheet interaction potentials, the slip length on the three materials obtained with FFMD may be up to five times larger compared to obtained AIMD (ESI†). Moreover, the mechanisms involved in water slippage may critically depend on the particular choice of the force field, as seen in the strong dependence of τF and λ on the values of the Lennard-Jones interaction and in the values of the partial charges (ESI†). This indicates that empirical interaction potentials, usually developed to describe bulk systems and sometimes combined with ad hoc mixing rules to describe interfaces, may be unreliable for the comparison of slippage mechanisms across different systems.
Overall, in this work we have unraveled the dependence of water slippage on the energy landscape corrugation and on the correlation time of the density fluctuations of water on three prototypical two-dimensional materials. Obtaining such fundamental mechanistic insights from 100 ps-long AIMDs through a systematic investigation of water slippage under different confinement regimes on different materials was not possible just a few years ago, when a first AIMD study of liquid/solid friction appeared,21 and has been made possible thanks to recent advances in computational power and algorithm developments.51 With the aim of designing materials with faster water slip one can envision different mechanisms to further reduce the friction, not only by decreasing the corrugation of the energy landscape but also by reducing the force correlation time. The recent discovery of a few thousands “exfoliatable” two-dimensional materials and the rise of van der Waals layered heterostructures may present a number of additional routes to further improve water slip.52–55 So far graphene is the two-dimensional material which exhibits the lowest friction to water. It remains to be seen whether water slip on newly proposed two-dimensional materials or van der Waals heterostructures can be faster than on graphene. Furthermore, inclusion of non-additive many-body effects may reveal new routes to enhance slp, as a result of repulsive dispersion forces arising from confinement.56,57 It will be interesting to extend AIMD studies to other domains of nanofluidics including problems in osmotic energy conversion, in a time where also ionic transport is accessible experimentally under extreme confinement regimes.
AIMD simulations within the Born–Oppenheimer approximation were performed in the canonical ensemble for about 100 ps at 330 K to partly compensate for the overstructuring of water with the optB88-vdW functional.60,61 The initial atomic configurations for each system were obtained from pre-equilibrated FFMD trajectories. Additional AIMD and FFMD simulations ensured that our results were not affected by the finite size and time-scales accessible to AIMD and by the type of ensemble chosen (ESI†).
Footnotes |
† Electronic supplementary information (ESI) available: Further computational and theoretical details, as well as tests on the structure and dynamics of water at the interface with graphene, hBN and MoS2 from force field and ab initio simulations. See DOI: 10.1039/D0NR02511A |
‡ kB being Boltzmann constant and T room temperature. |
§ The magnitude of the peaks for the smallest confinement widths (Fig. 1b–d) and the presence of one or two additional peaks beyond the contact layer (Fig. 1f–h) may be marginally affected by the number of water molecules present in our simulations, and as such care must be taken while interpreting these differences in the density profiles. |
This journal is © The Royal Society of Chemistry 2020 |