The Role of Hydrogen-Bonded Bridges for the Diffusion of Water Confined Between Graphene Oxide Layers

Among the carbon-based two-dimensional materials, graphene oxide (GO) has been attracting a growing interest due to its capability to be utilized in the field of water remediation. Therefore, an atomistic understanding of the transport properties of water in layered GO is pivotal for the development of novel GO membranes. Surprisingly, the very issue of the two-dimensional self-diffusion of water confined between two GO sheets appears to be controversial and simulations showing either a slow-down or no effect have been reported. In any case the formation of Hydrogen bonds, i.e. among the confined water and between water and the GO functional groups, was identified to control diffusion. However, results of molecular dynamics simulations heavily depend

on the used forces. Density functional theory and empirical force fields are on opposite when it comes to accuracy and numerical costs. As a compromise in the present study we performed molecular dynamics simulations using a density functional theory-based tight method (xTB) to investigate the diffusion of water confined between GO sheets.
Specifically, we considered six GO/water models, differing in the ordering of epoxide and hydroxyl groups as well as in the thickness of the water layer. For these models, having GO inter-layer distances between 8 and 12 Å we find a reduction of the diffusion coefficient by a factor in between two and three as compared with bulk water. One possible origin of this effect is the temporary trapping of water within Hydrogen-bonded water bridges between the GO sheets.

Introduction
There is a tremendous interest in graphene and graphene-based materials triggered by their exceptional electronic, optical, and mechanical properties. 1 Graphene oxide (GO), an oxidized form of graphene, is one of these materials. The fact that the properties of GO can be tuned by modifying its structure, in terms of the types, ratio, and spatial arrangement of functional groups, opens the door for rational design for a wide spectrum of potential applications, e.g., water remediation, energy-harvesting and -storing devices, electrocatalysis, gas separation, and proton-conducting membranes in fuel cells. [2][3][4][5][6][7][8][9] The growing number of studies on the structure of and processes at the GO/water interface underlines the importance of understanding the interactions between GO and interfacial water. For example, in contrast to the non-oxidized carbon materials, GO is known to be soluble in water and some other organic solvents. Nevertheless, the dominant interactions in GO solutions and their correlation with the two-dimensional (2D) nature of the GO flakes are not fully understood yet. For instance, Dimiev and co-workers 2,10 have argued that associating the GO solutions nature and stability with the formal macroscopic parameters is not the appropriate approach to explain the properties of such solutions. Instead, they proposed that GO solutions should be treated as true solutions where the nature of the solution is governed only by the chemical interactions between the solute functional groups and the solvent molecules. In other words, the physicochemical properties of GO solutions can be understood by focusing mainly on the GO/liquid interface. 2,10 This emphasizes the importance of modeling aqueous GO solutions at the atomic level. [11][12][13][14][15] Especially, the water transport properties in GO membranes have been of great interest, but the results of the relevant studies, as well as their interpretations, are controversial. Here an important parameter is the inter-layer spacing. Using neutron scattering, Buchsteiner et al. 16 have demonstrated that the inter-layer spacing in multi-layer GO ranges between 7 and 11 Å depending on the humidity level. Talyzin et al. 17 using X-ray diffraction have found that immersion of GO in liquid water increases the inter-layer distance from 8 to 12 Å upon changing the temperature.
An extraordinarily rapid permeation of water through thick GO membranes has been reported by Nair et al. 18,19 and explained in the light of the capillary pressure effect due to the relatively short inter-layer spacing in their experiments (from 6 to 10 Å). These findings were substantiated using molecular dynamics (MD) simulations. In general, there is a number of MD studies of water transport in GO. 20-23 Devanathan et al. 24 have calculated the diffusion coefficient of the water intercalated between GO layers. Their work has shown that the 2D channels between GO sheets (with inter-layer spacing ≈ 15 Å) slow down the diffusion of water by an order of magnitude. They attributed this to the strong Hydrogen bonds (HBs) formed between water and the OH groups of GO. In stark contrast with this study, Mouhat et al. 25 have performed ab initio MD simulations using Density Functional Theory (DFT) of hydrated GO (with inter-layer distance around 16 Å) and found that the diffusion of water at the GO/water interface is as fast as in bulk water.
Since a clear understanding of the effect of GO confinement on water diffusion is of fundamental importance, e.g., to rationalize permeation of molecular species through GO membranes, there is a need for a systematic investigation in particular as far as the amount of confined water is concerned. The present study is addressing this issue by presenting MD simulations for several hydrated GO models with different inter-layer distances and water content. In order to properly describe hydrogen bond formation and dynamics we have chosen DFT based tight-binding theory 26 as a computationally efficient compromise between empirical force fields and DFT (see also Refs. 27,28). The results have been thoroughly analyzed and discussed, in terms of the structural properties, H-bonding, and the diffusion of the confined water. Inspired by the suggestion of a H-bonding network between functional groups of GO by Dreyer et al., 3,29,30 we focused on the analysis of potential water bridges that link the GO layers through chains of H-bonds. In fact it is found that these H-bonded bridges provide a key to the atomistic understanding of the slow-down of self-diffusion of water confined between GO layers.

Model Systems
Constructing realistic models for GO is crucial for understanding its behavior and properties. Despite the large number of experimental (e.g., NMR, FTIR, X-ray photoelectron spectroscopy) and theoretical studies that aimed at revealing the possible GO structures, [30][31][32][33][34][35] the issue is still controversial, with many GO models being used currently, see for example Sheka et al. 36 The state of affairs can be summarized as follows: (1) Different types of oxygencontaining groups are distributed on the graphene surface and a relatively small number of these groups could be attached to the edges. (2) There is no correlation between the oxygen content and the types of these functional groups. For the present study the construction of model systems followed the general features of the Lerf-Klinowski model 37 that fulfills the above-mentioned criteria and, hence, is the most accepted model so far.

Anhydrous GO Models
We considered two different GO models with identical chemical composition but differences in the spatial arrangement of the functional groups. Starting point is an orthorhombic periodic supercell of graphene consisting of 72 carbon atoms. In each model, 12 hydroxyl and 6 epoxide groups were grafted on the basal plane of graphene. This means that there are 24, out of 72, sp 3 C atoms that are bonded to oxygen functional groups resulting in a functionalization rate of 25%, a typical value for GO. 25,38,39 The OH groups have been systematically placed in a way that two adjacent C atoms bear a pair of OH groups coming out of the opposite sides of the graphene sheet (cf. Figure 1c). It is noteworthy that all the above-mentioned features of our models align with the GO models of Mouhat et al. 25 Two configurations have been considered in our study as shown in Figure 1. In the first GO model, OGO, there is a fully oxidized domain surrounded by a graphene-like region.
In the fully oxidized island the functional groups have been arranged in an ordered manner (similar to the chain-like structure described by Yan and Chou 40 ). In the second GO model, RGO, the positions of the functional groups have been chosen at random, where the oxidized and graphene-like regions intersect each other. Note that these two models present limiting situations and in principle a more elaborate model could include averaging with respect to variations in the specific arrangement of the functional groups (see, e.g., Ref. 25) or simultaneously occurring random and ordered regions.

Hydrated GO Models
In order to study the proporties of confined water, we designed hydrated models using double layers of the GO structures introduced above. For each GO structure (i.e., OGO and RGO),  we constructed three hydrated models that vary in the number of water layers, from one to three, and, hence, in the inter-layer distances. A water layer contains 20 molecules, a number chosen based on geometrical consideration to cover the GO surface area.
The hydrated GO systems were built using the PACKMOL program. 41 The shortest possible distance between the water molecules and the GO surfaces has been imposed to be at least 2 Å to avoid steric overload because of the functional groups on the surface of GO. The distance between the two GO sheets has been determined based on the number of water layers. We assumed a thickness of the water layer of ≈ 3 Å. Therefore, the interlayer distance (i.e., the distance between the centers of mass of the two GO sheets in the z direction) in the models with 1-3 water layers has been set to be 8-14 Å. The models will be labeled OGOn and RGOn with n = 1, 2, 3.
Since we will apply periodic boundary conditions (PBC) in all directions (see section 2.2), we added water above and below the double GO layer to ensure that the taking into account the images of the main cell there is always the same amount of confined water. The number of layers that we refer to when we discuss our models is the number of water layers that is confined between the two GO sheets, while the actual amount of water in the simulation cell equals twice that number. The different models are shown in Figure 2. Figure 2: Hydrated GO Models after equilibration: (a) GO with 1 layer of confined water (i.e., OGO1 and RGO1 with average inter-layer distances 8.0 Å and 7.5 Å, respectively), (b) GO with 2 layers of confined water (i.e., OGO2 and RGO2 with average inter-layer distances 9.3 Å and 9.2 Å, respectively), (c) GO with 3 layers of confined water (i.e., OGO3 and RGO3 with average inter-layer distances 11.7 Å and 11.8 Å, respectively).

MD Simulations Details
As a first step, we performed geometry optimization and a cell relaxation for each anhydrous GO structure (i.e., OGO and RGO) to adjust the lattice parameters a and b after adding the oxygen functional groups that induce a mechanical strain on the graphene surface. Once the GO structures were fully relaxed, 30 ps MD production trajectories have been produced at T = 300 K. The simulations were performed with 0.5 fs time step in the NVT ensemble with canonical sampling through a velocity rescaling (CSVR) 42 thermostat with a time constant of 100 fs.
After preparing the preliminary hydrated GO models (as described above), the procedure is quite similar to that of the anhydrous GO structures. To be specific, geometry and cell relaxation simulations were performed to relax the water molecules and to ensure their proper chemical configuration relative to the GO functional groups. In contrast to the anhydrous GO case, in the hydrated GO cell relaxation, we applied constraints on the a and b cell vectors to keep the GO structures preserved in the 2D plane and only the c vector was allowed to change. Afterwards, two stages of water equilibration were carried out. Namely, a 1 ps trajectory was generated using massive CSVR thermostats (with a time constant of 10 fs) followed by a 5 ps trajectory using global equilibration (with a time constant of 100 fs). Finally, a 30 ps production trajectory has been produced for each hydrated GO system using the same conditions that described above for the anhydrous systems.
All the simulations were carried out using the extended tight binding (xTB) method 26 implemented in the CP2K package. 43,44 The dispersion correction D3 has been utilized. 45 PBC in all directions have been applied in all calculations. The molecular structures and MD trajectories were visualized using the Visual Molecular Dynamics (VMD) software. 46 Finally, all the results discussed below belong to the 30 ps production trajectories generated for each of the anhydrous and hydrated GO models.

Structural Properties
Before investigating the interaction between GO and water we inspected the structure of the GO models in water as well as in vacuum to see whether the presence of water affects the structural stability of GO. To this end four geometrical parameters, namely the C-C bond lengths (d CC ), C-C-C angles (θ CCC ), C-O-C (epoxide) angles (θ COC ) and C-O-H (hydroxide) angles (θ COH ), have been chosen. The distributions of d CC and θ CCC show how the pristine graphene surface is perturbed by the presence of the epoxide and hydroxide groups, while θ COC and θ COH tell about the functional groups themselves and whether the presence of water molecules affects their structure.
We observed that the distributions of these characteristic structural parameters do not show a strong dependence on the surrounding medium. In fact there is no significant change in the distances and angles of GO in vacuum versus in water, no matter how many water layers are included. Therefore, in Figure 3, we present the results for the ordered and random GO models in vacuum and for one water case only.
The distribution of d CC in Figure 3a shows an elongation of these bonds, relative to the reference distance in pristine graphene (1.42 Å), as a result of the presence of the functional groups. Obviously, for the given degree of functionalization the d CC distribution depends on the spatial arrangement of the functional groups. The pronounced peak in the OGO models around the reference value corresponds to the graphene-like domain that is completely distinct from the oxidized region (small shoulder around 1.5 Å). The situation is reversed in the RGO models. Here, the main peak is located around 1.5 Å which approximately corresponds to the distance between sp 3 C atoms. A recent study has reported that, in contrast to our findings here, the distribution of d CC in random GO models is more diffuse. 25 This discrepancy arises from the difference between the random GO configurations. Particularly, in the present random model we avoided adding functional groups to the edges of the periodic box, while the model in Ref. 25 did not have such a limitation. As a consequence our random models have relatively large graphene-like areas of sp 2 C atoms (cf. Figure 1).
Unlike the distances, d CC , the distributions of the remaining characteristic parameters

HBs between GO and Water
One of the main features of GO is its ability to form HBs with the surrounding molecules.
Several studies 10,24,25 have shown that many of the chemical and physical properties of the GO/water interface as well as the dynamical behavior of the confined water can be explained in terms of the nature of these HBs. Therefore, we performed a thorough analysis of HBs for the present hydrated GO models.
To this end we have used geometrical criteria 49  and inter-molecular HBs (i.e., HBs between functional groups of GO and the surrounding water or with functional groups of the other GO layer in the case of very short inter-layer distance). However, due to the tight confinement of water in-between the GO layers, the inter-molecular GO-water HBs are assumed to be the dominant kind of interactions at the GO/water interface, in particular when it comes to water diffusion. Therefore, we focus our analysis on these HBs. Figure 4 summarizes the results, giving in panels (a-c) the number of HBs of type O ep , HO acc and HO don , as the percentage of the maximum possible number.
Inspection of these data reveals that: 1. For all types of HBs, the random models form a larger number of HBs with the H 2 O molecules than the ordered models do. This is not surprising since the functional groups of the ordered GO are quite close to each other and, hence, can form more intra-molecular HBs that compete with the inter-molecular HBs between functional groups of the GO sheet.
2. Regardless of the arrangement of the functional groups on the graphene surface, the number of HBs in the cases of one layer of confined water is significantly smaller than in the cases of two and three water layers. This is the result of the competition for forming HBs between the water and the two GO layers, on the one hand side, and between water molecules, on the other hand side.

Water Diffusion
The lateral (two-dimensional, 2D) diffusion coefficient, D W , of the confined water has been calculated from the mean square displacement for all the six water models (for details, see the Supplementary   bridges between the two parallel GO layers. These HB bridges seem to play a crucial role in slowing down the diffusion of the confined water. In the next section, we discuss these HB bridges and their role in detail.

H-bonded Bridges Between GO Layers
A HB bridge is defined as a chain of GO-water (and, in some cases, water-water) intermolecular HBs that co-exist and extend from one GO sheet to the other. These bridges will be classified into different types based on the number of involved H 2 O molecules. Thereby, an n-th order bridge exists if there are n water molecules bridging the two GO layers. Although longer H-bonded chains may in principle exist, we will concentrate on the shortest possible bridge, which is of order n for OGOn and RGOn.
First, we identified the aforementioned bridges using the geometric criteria given above and counted their occurrence in the corresponding trajectories. Figure 6 displays histograms for the investigated bridges in relation to the number of possible donor-acceptor groups (pairs of functional groups on different GO layers close enough to be connected by an n-th order bridge). As can be seen in Figure 6a, for the shortest inter-layer distance there is a clear difference between the OGO1 and RGO1 models. The weighted average (given in table 1) of the number of bridges/pair of functional groups in the random system RGO1 is 0.93, whereas it is 0.62 for OGO1. Increasing the amount of confined water (hence, the inter-layer distance) from one to two layers increases the probability of forming more HB bridges involving different water molecules. Indeed the number of 2nd order bridges in OGO2 increases by about 0.2 as compared to OGO1 (cf. the blue lines in Figure 6a vs. 6b). In contrast, the number of 2nd order bridges in RGO2 increases only by about 0.01. In general both distributions are broader than in case of a single water layer. Figure 6c shows that for 3rd order bridges the difference between the ordered and random systems essentially vanishes (difference is 0.01). Thus, we can conclude that: (1) the spatial arrangement of the GO functional groups has an insignificant effect on the number of the HB bridges as the distance between GO layers increases (> 12 Å), and (2) to have the maximum number of HB bridges, the optimum number of water layers is 1-2 for the random GO, and 2-3 for the ordered GO.  Table 1. From the table one notices that, in all studied models, the OH don -OH acc bridge pattern (cf. Figure 7b) has by far the highest occurrence probability (cf. the highlighted row in Table 1). The H-bonded bridges are very dynamic, they form and break up many times during the trajectory. Therefore, we did not attempt to quantify their average lifetime. Instead, in order to better understand their effect on water diffusion, we investigated the distribution of HB donor-acceptor distances along the trajectory. Figure 8 shows the result for the dominant OH don -OH acc bridge and the case of a single water layer. Here, R 1 and R 2 is the distance between the oxygen of the bridging water and the oxygen atom in OH don and OH acc , respectively. The actual time evolution of R 1 and R 2 are given in Figure ?? of the Supplementary Information. In Figure 8 the range of possible HB bridges according to our definition is shown by the rectangle.
From Figure 8 we can conclude that the HB bridge of the random GO (panel (b)) is more stable than its analogue in the ordered GO (panel (a)). Specifically, the OH acc acceptor in the case of OGO1 is weaker than its counterpart in RGO1 (tendency for larger R 2 values). This difference can be traced to the spatial arrangement of the GO functional groups in OGO1, which are closer to each other compared to RGO1. As a result neighboring functional groups on the same GO layer compete with functional groups from the other GO layer to form a HB with the water layer. In both cases the distribution of R 1 values indicates that a

Conclusion
We have investigated the properties of two-dimensional water confined in between two graphene oxide layers using extended tight-binding (xTB) based molecular dynamics. The main focus has been on the influence of confinement on the self-diffusion of water, where previously published work came to conflicting conclusions. 24,25 In these studies different structural models as well as different methods for the forces (DFT versus force field) were used. Therefore our aim has not been a direct comparison. Instead we approached the issue by designing six models capturing different layer distances and ordering of the GO functional groups. The interlayer distances were in the range of experimental values reported previously. 16,17 Using xTB whose accuracy is expected to be intermediate between DFT and force fields, we provide a computational frame that can be systematically extended to larger systems in future work. In the present work our main focus has been on the role of H-bonded bridges connecting the GO layers for water diffusion. Although suggested more than ten years ago, 3 this motif has never been explored by computational studies so far.
For the ordered and random GO models we found that there is no significant structural variation upon solvation. The oxygen functions (i.e., epoxide and hydroxide groups) form 30-70 % of the possible HBs with the surrounding water molecules, with smaller values for the respective ordered GO models. Based on geometrical criteria, most of GO-water HBs are of moderate strength.
Despite the dynamic nature of these HBs, the confined water diffusion between the GO layers has been found to be slower than in bulk water by a factor in between two and three.
This holds irrespective of the number of water layers and ordering of the GO sheets. Besides H-bonding to a single GO layer, H-bonded bridges connecting the two layers have been identified for all models. Their number is larger for random GO sheets and in that case also decreases with increasing interlayer separation. Both types of HBs, i.e. GO-water and GO-water-GO, cause a temporary trapping of water molecules. However, in particular the bridges are part of the H-bonded intra-water network. Hence there is a competition between intra-water liberating forces and GO trapping, which was identified as a mechanism behind the slow down of the diffusion. It will be interesting to explore even larger systems such as to approach the bulk-like limit and to connect to experimental data on diffusion of molecular species.

Supporting Information Available
The SI contains more details on the water self-diffusion calculations, a time evolution representation of the data of Figure 8, and two movies that visualize the trapping of water molecules in the HB bridges.

SI Water Diffusion Calculations
The data for the lateral self-diffusion coefficient, D W , shown in Figure 5 have been calculated as follows. First, the 2D mean squared displacement (MSD) (i.e., MSD in the xy plane for the water in between the two GO layers of the simulation box) have been obtained as follows (N number of MD steps): where r i are the atomic coordinates. Then, the D W values have been estimated from the slope of the linear fit of MSD with respect to time (cf. Figure S1 and Table S1) according to The error in the slope has been obtained from the sum of squared residuals of the regression SSR and the time points t i according to One might argue that behavior of the trapped H 2 O molecules is not likely to be of random walk like Brownian motion and a description in terms of sub-or superdiffusion is more appropriate. 1,2 In order to test this assumption, we have also fitted the MSD data to the equation To this end we have used python module "Diffusion Analysis MD Simulations". 3,4 The results are given in Figure S2 and Table S1. First, we notice that similar to the normal diffusion model, all GO/water models show a decrease of the diffusion coefficient with respect to bulk water. In contrast to normal diffusion, OGO1, has a different behavior as compared to the other cases, i.e. the slow-down is not that pronounced. Moreover, while all other cases show subdiffusive behavior -as one might have expected for diffusion in an inhomogenous environment 1 -OGO1 is superdiffusive. Based on the analysis of OGO1 there appears to be no reason for a qualitatively different behavior. Therefore, we tentatively assign this finding to be caused by insufficient sampling and postpone further discussion of anomalous diffusion to a future study.   Table S1: The values of the lateral diffusion coefficients, D W (×10 −5 cm 2 /s) and D α W (×10 −5 cm 2 /s α ), of the bulk water and the confined water in the six GO/water models. The table shows the results obtained from both the linear ( Figure S1) and non-linear ( Figure S2) fits of the MSD.

SII H-bonded bridges between GO layers:
SII.i Time evolution of the distances whose distribution is given in Figure 8 in the main text.  Figure S3: Time evolution of R 1 and R 2 values in OH don -OH acc HB bridges of (a) OGO1, (b) RGO1. R 1 is the donor-acceptor distance in the GO-water HB OH don , while R 2 is the donor-acceptor distance in the GO-water HB OH acc .

SII.ii Supplementary Movies
In order to provide a visualization of the mechanism by which GO inhibits the confined water self-diffusion, we have generated two short movies using the Visual Molecular Dynamics (VMD) software. 5