# Structure and chemistry of graphene oxide in liquid water from first principles

Félix Mouhat, François-Xavier Coudert, Marie-Laure Bocquet

<figure><img src="image_1.png" /></figure>

## To cite this version:

Félix Mouhat, François-Xavier Coudert, Marie-Laure Bocquet. Structure and chemistry of graphene oxide in liquid water from first principles. Nature Communications, 2020, 11, pp.1566. 10.1038/s41467-020-15381-y. hal-02520022

HAL Id: hal-02520022

https://hal.science/hal-02520022

Submitted on 26 Mar 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ARTICLE

https://doi.org/10.1038/s41467-020-15381-y

OPEN

<figure><img src="image_3.png" /></figure>

Check for updates

# Structure and chemistry of graphene oxide in liquid water from first principles

Félix Mouhat<fnref n="1" />, François-Xavier Coudert<fnref n="1" />,<fnref n="2" /> & Marie-Laure Bocquet<fnref n="1" />,³

Graphene oxide is a rising star among 2D materials, yet its interaction with liquid water remains a fundamentally open question: experimental characterization at the atomic scale is difficult, and modeling by classical approaches cannot properly describe chemical reactivity. Here, we bridge the gap between simple computational models and complex experimental systems, by realistic first-principles molecular simulations of graphene oxide (GO) in liquid water. We construct chemically accurate GO models and study their behavior in water, showing that oxygen-bearing functional groups (hydroxyl and epoxides) are preferentially clustered on the graphene oxide layer. We demonstrated the specific properties of GO in water, an unusual combination of both hydrophilicity and fast water dynamics. Finally, we evidence that GO is chemically active in water, acquiring an average negative charge of the order of 10 mC m⁻<fnref n="2" />. The ab initio modeling highlights the uniqueness of GO structures for applications as innovative membranes for desalination and water purification.

<footnote n="1">¹PASTEUR, Département de Chimie, École Normale Supérieure, PSL University, Sorbonne Université, CNRS, 24 Rue Lhomond 75005 Paris, France.</footnote>

<footnote n="2">²Chimie ParisTech, PSL University, CNRS, Institut de Recherche de Chimie Paris, 75005 Paris, France. ☉email: marie-laure.bocquet@ens.fr</footnote>
**G**raphene oxide (GO) is a graphene-based material that has gained significant interest in the last two decades¹,² due to its straightforward, scalable, and low-cost synthesis. It has been proposed for numerous applications: for instance, GO is a promising material as an efficient sieve for water remediation³⁻⁷ and for sustainable energy production via fuel cells⁸,⁹. Such properties strongly depend on the chemical groups present at the surface and versatile synthesis conditions using different reagents and oxidation durations have been proposed to allow a slight variation of the chemical composition of the material¹⁰⁻¹³. Consequently, there is a large variety of GO surfaces with a typical O/C ratio varying from 28 to 36%¹⁴. This ratio is known to decrease down to O/C = 5–10% when a thermal exfoliation is applied to the GO to produce single sheets of functionalized graphene, as some chemical functions are removed¹⁵,¹⁶. Using different characterization techniques, such as solid NMR, XPS, or Raman spectroscopy measurements, many theoretical models attempting to describe the GO surface have been proposed so far. The most widely used model is the one proposed by Lefr and Klinowski¹⁷,¹⁸ originally for graphite oxide, and which has since then been largely used in the literature to describe monolayers of GO. Within this model, the layer is depicted as a random distribution of flat aromatic regions with unoxidized benzene rings, wrinkled regions containing hydroxyl or 1,2-ethers (epoxide) groups, and carboxylic acids grafted on the edges of the sheet. Despite some further refinements, such as the anti position of the hydroxyl pairs in the basal plane¹⁹,²⁰, there is to the best of our knowledge, no clear statement about the spatial arrangement of the oxidized functions along the graphene layer. Are they randomly distributed along the surface or not? What are the consequences on the stability of the material and its chemical properties, especially when it interacts with water? How does that impact filtration properties across GO membranes?

Indeed, the interaction of GO with water is particularly important due to the observed superpermeability of GO membranes to water molecules and its understanding requests modeling at the atomic scale. To date there are mostly classical MD simulations, hence treating the oxygen chemical groups as passive in water²¹⁻²⁵.

In this manuscript, we propose a realistic model of the GO basal plane surface, describing both water and the chemical groups at the electronic level. Within this approach, we generate manifold GO replicas at a given O/C ratio, with functions randomly distributed or not over the surface. Our as generated GO layers are first studied without solvent, and among the stable ones (see below), six of them are placed in liquid water at room temperature by means of long ab initio molecular dynamics (AIMD) simulations. Next, we perform a thorough statistical analysis by averaging over time and over replicas in order to quantitatively measures some of the physical and chemical properties of valuable interest for the community working on 2D materials and water/solid interfaces.

## Results

**Realistic graphene oxide models.** We constructed realistic GO models initially in anhydrous conditions, based on an orthohomobic periodic supercell structure of graphene comprising 72 carbon atoms. Starting from this pristine basal plane of graphene, we grafted 18 chemical functions to carbon atoms: 12 hydroxyl groups and 6 epoxide functions, which corresponds to a 2:1 ratio. This choice is consistent with suggested chemical formulas²⁶ and recent theoretical studies²⁷. By construction with periodic systems, we discard the edge functionalization with carboxylic functions, which are experimentally known to be in smaller proportion, as compared with hydroxyls and epoxide ones.

Hydroxyl groups are systematically placed in pairs, respecting an anticonfiguration (with the two groups on opposite sides of the sheet). The total number of functional groups is chosen so that there are 24 sp³ carbon atoms linked to an oxygen atom (out of 72), and the functionalization rate corresponds to O/C = 18/72 = 25%, a value in the low range of typical GO flakes²⁸,²⁹. We note that the choice of O/C ratio results from a compromise: taking values closer to the higher experimental limit would have decreased the number of possible independent replicas (for the semiordered case because the number of C atoms bonded to an oxygen atoms acquire the majority status); while a lower ratio would require more simulation replicas in order to measure statistically significant changes in behavior.

With this procedure, we generated ten GO structures, with the exact same chemical composition. Furthermore, for half of the structures, herein referred to as “random models”, the positions of the hydroxyl and epoxide groups were chosen randomly among all C atoms of the graphene sheet. For the other half structures, the groups were kept concentrated in one half of the graphene sheet. We call these the “semiordered models”—although the exact positions of the functional groups are stochastic, these models contain some correlation. By looking at the systems generated in each case (see Fig. 1), the semiordered models present by construction a nanoscale patterning, where the GO shows aromatic regions, percolating to form graphene-like wires, as was suggested in some past works³⁰,³¹.

For the ten GO structures displayed in Fig. 1, we first checked their structural stability and energetics in vacuum by carrying out a series of structure optimizations (see details in the “Methods”). We report in Table 1 the relative total energies, including both electronic and atomic contributions, of the various models after geometry optimization (atomic positions and cell parameters). The lowest-energy configuration was chosen as reference point, with $E_{\text{ref}} = -266.270$ eV per C atom. All models were sorted from most energetically stable to least stable, and an average energy was calculated for each type (random and semiordered). We first note that there is significant dispersion of the computed energies, and that it is more important for the random structures than for the semiordered ones—possibly due to the constraints on the chemical groups. We also remark that the semiordered structures are on average 222 kJ mol⁻¹ or equivalently 32 meV per C atom more stable than the random ones, hinting that the former type could be the dominant form of GO synthesized experimentally. This energetic stability is indeed substantial and, interestingly compares well with the 40 meV per C atom stabilization offered by the adsorption of a graphene layer on a strongly coupling metal-like Ru³²,³³.

Two different hypotheses, which do not exclude one another, can be formulated to account for this result. First, adding a chemical function on a carbon atom induces some strain around the carbon site. In the random structures, a significant part of the graphene lattice is stressed by the epoxide and alcohol functions, whereas in the semiordered case, some regions remain unperturbed, and are therefore lower in energy. There is an energetic gain to concentrate the defects in localized regions, instead of applying an homogeneous deformation. Secondly, in semiordered structures, the closest proximity of the oxygen-bearing functional groups leads to a higher number of internal hydrogen bonds stabilizing the edifice. Such findings have been suggested previously by theoretical work using different methodologies, such as cluster models studied at the Density Functional Theory (DFT) level³⁴, or larger systems modeled with reactive force fields³⁵.

For the remaining part of this paper, we will focus on the three most stable structures obtained for each type of model (marked with a in Table 1). Their chemical structure is sketched in Fig. 1,
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Models of graphene oxide studied in this work. Semiordered (SO, top) and random (R, bottom) models of graphene oxide generated and studied as part of this work. Each material is designed according to a 25% O/C functionalization rate, and the oxygen-bearing groups are hydroxyls and epoxides with a 2:1 ratio. Gray hexagons are a guide for the eye, indicating the areas of pristine graphene, where no chemical function is grafted. Double C=C bonds are not drawn as double bonds, for the sake of clarity. All the configurations are represented as a 2 × 2 supercell of the graphene oxide model (the cell size is indicated in dotted lines on the top-left panel). Numbering of the models corresponds to Table 1. Intramolecular hydrogen bonds in GO are indicated by light-blue dotted lines, and the special reactive motif of three consecutive hydroxyl groups in the semiordered 4 model is colored in light green.</figcaption></figure>

<table><caption>**Table 1** Relative total energy *E* − *E*<sub>ref</sub> (in meV per carbon atom) of the ten generated semiordered or random GO structures with respect to the most stable system, labeled "ref".</caption><thead><tr><th>GO configuration</th><th>Electronic energy (meV per C atom)</th></tr></thead><tbody><tr><td>Semiordered 3<sup>a</sup> (ref)</td><td>0</td></tr><tr><td>Semiordered 4<sup>a</sup></td><td>3</td></tr><tr><td>Semiordered 1<sup>a</sup></td><td>13</td></tr><tr><td>Semiordered 5</td><td>13</td></tr><tr><td>Semiordered 2</td><td>33</td></tr><tr><td>Average for semiordered</td><td>16</td></tr><tr><td>Random 2<sup>a</sup></td><td>12</td></tr><tr><td>Random 1<sup>a</sup></td><td>27</td></tr><tr><td>Random 4<sup>a</sup></td><td>40</td></tr><tr><td>Random 5'</td><td>72</td></tr><tr><td>Random 3'</td><td>81</td></tr><tr><td>Average for random</td><td>48</td></tr></tbody></table>

All structures are relaxed using the CP2K code (see "Methods").

<sup>a</sup>The three most stable structures of each GO type that are further studied in this paper and displayed in Fig. 1.

where gray hexagons are used as guide for the eye to indicate the remaining regions of pristine graphene in GO models. It becomes clear that such regions are scarce and small in the random models (bottom panel), while in the semiordered models there are large pristine sectors that percolate to form one-dimensional wires.

**Structural properties.** In the following, we analyze the structural properties of the selected semiordered and random GO configurations, first in vacuum, and then in presence of liquid water. These properties were obtained through ab AIMD at room temperature (see Supplementary Fig. 2 and the "Methods" for details). The size of the simulation box along the z-axis, perpendicular to the GO sheet, imposes the scale of confinement of the water. We chose a value (c ≃ 14.5 Å) in good agreement with the typical interlayer spacing for GO laminates that swell in water (~14 Å)<sup>36,37</sup>. Extending first the stability analysis in the presence of water, the conclusions obtained in vacuum hold true: on average the stability of semiordered GO layers compared with random GO ones is increased, from 2.4 eV in vacuum to 3.3 eV in water (see details in Supplementary Table 1). Next, we looked at four characteristic distances and angles, displayed in Fig. 2. In a nutshell, the histograms of C–C distances and the CCC angles indicate how the original graphene layer is perturbed by the presence of chemical functions (alcohols and epoxides). Moreover, the distribution of the COC and the COH angles allows us to look at the structure of the functional groups themselves, and to see how the GO interacts with the surrounding H<sub>2</sub>O molecules.

We first point out that the *d*<sub>C–C</sub>, CCC and COC distributions for either the semiordered and random GO models do not significantly change with the environment, i.e., in vacuum or in the presence of water. Panels a–c of Fig. 2 are therefore only shown in water. As expected the C–C distance (Fig. 2a) is elongated by the presence of oxidizing groups with respect to the reference distance in graphene, *d*<sub>C–C,graphene</sub> = 1.42 Å<sup>38</sup>. Moreover, the semiordered and random models of GO give markedly
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Structural properties of graphene oxide. Histograms of the C–C distance *d*<sub>C–C</sub> (a), the CCC angle (b), and the COC angle (c) for semiaromatic (black) and random (red) graphene oxide models. Blue dashed lines correspond to the values in reference systems: *d*<sub>C–C</sub> = 1.42 Å in graphene, CCC<sub>sp<sup>3</sup></sub> = 109.28°, CCC<sub>sp<sup>2</sup></sub> = 120°, and COC = 60° for epoxide. d Histogram of the COH angle for the GO in vacuum (dashed lines) and GO solvated in water (solid lines). e Snapshots of the different H bonds (visualized as blue dashed lines) types classified in Table 3. The atoms of each type of chemical function (hydroxyl group or epoxide) involved in the H bonds between the surface and H<sub>2</sub>O are highlighted in orange.</figcaption></figure>

different distributions. There is a more pronounced peak around the graphene reference value for the semiaromatic configuration (black curve), corresponding to the regions that are not functionalized (i.e., pristine or graphene-like), whereas the distribution is more diffuse in the random case (red curve), as the chemical functions are spread over the entire GO sheet.

In contrast the CCC angles distributions (Fig. 2b) coincide for the semiaromatic and random models, exhibiting two peaks corresponding to sp<sup>2</sup> and sp<sup>3</sup> carbon atoms, with different peak heights. We tried to fit these distributions by two Gaussians:

$$f(x) = \frac{\lambda}{\sqrt{2\pi\sigma_{sp^3}^2}} e^{-\frac{(x-\mu_{sp^3})^2}{2\sigma_{sp^3}^2}} + \frac{1-\lambda}{\sqrt{2\pi\sigma_{sp^2}^2}} e^{-\frac{(x-\mu_{sp^2})^2}{2\sigma_{sp^2}^2}}, \quad (1)$$

where symbols (black circles and red diamonds) in Fig. 2 correspond to the fit. The main fitting parameter λ is an estimation of the functionalization rate O/C, which we expect to be around λ₀ = 0.25 (25% sp³ C atoms) if the model is appropriate. We report in Table 2 the parameters obtained from this analysis. A very good agreement is found between λ<sub>fit</sub> = 0.22–0.23 and λ₀. The values of μ<sub>sp³</sub> and μ<sub>sp²</sub> are in good agreement with the typical values tabulated for sp³ and sp² angles, confirming that the CCC angles distribution can be a useful tool to determine accurately the functionalization rate of GO. While such angle distributions may not be directly accessible experimentally, other quantities can be correlated to these angle distributions, such as bending vibration frequencies.

The distributions of epoxide COC angles, displayed in Fig. 2c, are sharp and well-defined, yet at the same time slightly shifted to values > 60°. Despite their stability, these functions are thus likely to interact with neighboring –OH groups or H₂O molecules. We also note the presence of small, broad peak at unexpectedly large

<table><caption>**Table 2** Fitting parameters of the CCC angle distributions for the semiaromatic and the random graphene oxide models, using the two Gaussian functions of Eq. (1).</caption><thead><tr><td>Fitting parameter</td><td>Semiaromatic</td><td>Random</td></tr></thead><tbody><tr><td>λ</td><td>0.23</td><td>0.22</td></tr><tr><td>μ<sub>sp³</sub></td><td>112.2</td><td>111.4</td></tr><tr><td>σ<sub>sp³</sub></td><td>4.3</td><td>3.6</td></tr><tr><td>μ<sub>sp²</sub></td><td>119.7</td><td>119.6</td></tr><tr><td>σ<sub>sp²</sub></td><td>2.6</td><td>2.8</td></tr></tbody></table>

values the COC angle. These correspond to specific situations, discussed in detail in the Supplementary Information, such as the formation in rare cases of a 1,3-epoxide (Supplementary Fig. 1) after rearrangement of the structure at *T* = 0 K (see Supplementary Fig. 2). Large COC angles are also seen in cases of strongly strained structures where two epoxides are near-neighbors on a same hexagonal pattern, inducing large local stress that stretches the involved C–C distances and therefore flattens the corresponding epoxides (see Supplementary Fig. 3).

Finally the COH angle distribution (Fig. 2d) is narrow and symmetric around COH = 105° in vacuum, relatively unaffected by the order of the functions. In sharp contrast, in the presence of water the peak shifts to higher values and displays a broad tail— with values up to COH = 160° being observed. This can be understood as a competition between intra- and intermolecular hydrogen bonds: in vacuum, to minimize its internal energy, some H bonds are formed between alcohols and epoxides functions of the GO. In the presence of water, in addition to the already present intramolecular H bonds, the GO can form new
intermolecular H bonds, between its oxygen-containing groups and the H₂O molecules of the liquid. This competition is a first indicator of the presence numerous H bonds between the GO sheet and the liquid water, suggesting that the GO might be reactive.

**Hydrogen bonds and dynamics.** We therefore performed a systematic analysis of the hydrogen bonds present, differentiating between donor (D) and acceptor (A) atom types using Chemfiles/cfiles³⁹,⁴⁰. To do so, we adopted selection criteria typical for H bonds of moderate strength, as they mostly are in proteins⁴¹: from the trajectory, we consider a potential H bond as formed if the D-A distance $d_{DA} \le 3.5$ Å, and $\overline{DHA} \le 30^\circ$. Figure 2e represents a zoom-in view of an MD snapshot focusing on the H bonds (marked in dotted blue lines) involving hydroxyl groups (left) and epoxide functions (right) of GO and interfacial water molecules.

In our analysis, we assume that two H bonds per epoxide (two lone pairs) and three per alcohol (two lone pairs + one donor site) can be formed with surrounding H₂O molecules. Therefore, the theoretical maximum number of water-GO hydrogen bonds in our systems is $2 \times 6 + 3 \times 12 = 48$. Table 3 shows that a large fraction of the functional groups of the GO (roughly 65%) are directly involved in a H bond, either with a neighboring function or with a H₂O molecule. As expected, the fraction of intramolecular H bonds is more important in the semiordered models than in the random structures since the functional groups are closer to each other in the highly functionalized regions. Because of the numerous lone pairs on oxygen atoms, the GO sheet accepts more H bonds that it donates. However, regarding the occupancy rates of the donor/acceptor sites of the surface in the right side of Table 3, epoxides seem to be discarded from the counting and hence less cooperative than hydroxyl groups. Almost all the protons at the GO surface are involved in a H bonds, at variance with the oxygens sites which are partially participating. Indeed, the oxygen position is more constrained, leaving less flexibility for the spontaneous formation/breaking of H bonds with surrounding H₂O molecules.

Finally, we turn to the impact of the GO on the transport properties of the liquid water. We estimated the water diffusion coefficients for semiordered GO models, and compared them with an ab initio MD simulation of bulk water in the same conditions (see Supplementary Fig. 9 and Supplementary Table 2). As displayed in Supplementary Fig. 9, we find that the average lateral diffusion value for water near the semiordered GO is around $0.67 \times 10^{-6}$ cm² s⁻¹ and a similar value for bulk water. This is an important and surprising conclusion, in stark contrast with the previous literature on confined water: contrary to other confined systems at 1–2 nm scale⁴²,⁴³, the water dynamics near the GO sheets is not slowed down, but its transport is as fast as in the bulk liquid. This occurs despite strong hydrogen bonding

<table><caption>**Table 3** Average percentage of intra- and intermolecular hydrogen bonds between graphene oxide and the surrounding H₂O molecules, and donor/acceptor sites occupancy analysis of the intermolecular H bonds.</caption><thead><tr><th></th><th></th><th>Semiordered</th><th>Random</th></tr></thead><tbody><tr><td rowspan="3">% H bonds</td><td>Intermolecular</td><td>49</td><td>57</td></tr><tr><td>Intramolecular</td><td>11</td><td>7</td></tr><tr><td>Total</td><td>60</td><td>64</td></tr><tr><td rowspan="2">% occupancy</td><td>Donor</td><td>70</td><td>79</td></tr><tr><td>Acceptor</td><td>36</td><td>37</td></tr></tbody></table>

Each type of H bond is illustrated in Fig. 2e.

between the water and the GO, hinting at the very dynamic nature of such H bonds.

**Reactive processes at the graphene oxide/water interface.** As shown above, there are numerous hydrogen bonds between GO and water: they are natural sites at which chemical reactions may be initiated. We analyzed the ab initio MD trajectories of semiordered and random GO models by dynamically checking the coordination of each oxygen atom of the GO sheet. Hence, we could identify three types of reactive events that we detail now.

Figure 3 illustrates two such processes: the ring opening of an epoxide (a) and the deprotonation of a hydroxyl group (b). In the first case, the C-O bond breaking leads to a zwitterionic form of the opened epoxide, with a negatively charged alcoholate group accompanied by a positive charge on an adjacent carbon atom. This process does not create net charge on the GO, unlike the second process observed, where deprotonation yields a negative charge on the GO sheet, balanced by an excess proton in the liquid water (hydronium cation). These two types of chemical events are the common ones, occurring in all our MD simulations, regardless of the GO model. We highlight, however, that we did not observe disruption of the C-C bond of an epoxide leading to a 1,2-ether function, as proposed in a recent study⁴⁴.

We also observe a third type of chemical reaction, much rarer, depicted in Fig. 4. In one of the semiordered models, the distribution of OH groups results in a close proximity of three hydroxyls, with two intramolecular hydrogen bonds between them. This chain of H bonds is suitable for proton transfer as in the Grotthuss mechanism⁴⁵, and indeed, we observed the proton of one OH group jumping to the next, forming an adsorbed water molecule and an alcoholate. The desorption of this adsorbed water molecule into the liquid is a dehydration reaction and is concerted with the transformation of the alcoholate into an epoxide. This highlights the possible lability over time of the oxygen-bearing functional groups in hydrated GO, affecting both the C/O ratio and the distribution of hydroxyl and epoxide groups. The structural motif that results in dehydration (highlighted in light green in Fig. 1) is only present once in our models, in the semiordered model 4. In other models, three consecutive hydroxyls pointing to the same direction can sometimes be identified but, unlike in model 4, they are participating in other intramolecular H bonds with neighboring epoxide groups (see Fig. 1).

Moreover, in contrast with simulations in liquid water, reactive events are scarce for GO replicas simulated in vacuum: no surface hydroxyl deprotonation, less frequent epoxide opening. For the semiordered model 4 exhibiting the reactive motif, we observe that dehydration also occurs in vacuum, but at a longer time (see Supplementary Section II.1 and Supplementary Fig. 3 for details).

**Charge of graphene oxide in water.** We now aim at quantifying these chemical events over time, and measure the net charge carried by the GO in presence of water. Figure 5 displays the time evolution of the presence of neutral and charged species, both in the GO sheet and the interfacial water, for the three semiordered GO models (data for three random models are displayed in Supplementary Fig. 1). Neutral species are hydroxyl and epoxide groups, while charged species are alcoholate groups and hydrogen cations. We first observed that, over the course of our MD simulations, the three different semiordered replicas behave very differently over time—this shows the impact of functional group distribution on the properties of the GO, and justifies the need for replicas in modeling GO. It is noticeable that all replicas start with five epoxides and not six showing that the initial geometry optimization process has resulted into one epoxide opening,
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Snapshots taken from molecular dynamics simulations. These snapshots are taken from trajectories of semiaromatic graphene oxide models in liquid water. **a** Opening of an epoxide function, leading to a zwitterionic form of the graphene oxide, but creating no net charge of the GO sheet. **b** Deprotonation of a surface hydroxyl group, leading to a surface alcoholate (blue shaded circles) and an excess proton (orange) in the liquid water. The labile hydrogen atom is highlighted in green. Schematics are presented on the right side, for clarity.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4** Spontaneous dehydration of the graphene oxide. We highlight the proton transfer induced by the presence of a strong hydrogen bond network. **a**, **b**, and **c** are snapshots along the MD trajectory of the semiaromatic 4 model displayed in blue in Fig. 5. The leaving water molecule is colored in green. Schematics are presented on the right side, for clarity. See Fig. 3 for color scheme. This reactive motif of three consecutive hydroxyl groups linked by two intramolecular hydrogen bonds is only present in semiaromatic model 4 and colored in light green in Fig. 1.</figcaption></figure>

reducing the strain in the material. The semiaromatic model 1 (black) remains moderately reactive with a few closing and opening of epoxide functions. The semiaromatic model 4 (blue curve) appears to be the most reactive one, and displays the dehydration mechanism described above. Finally the semiaromatic model 3 (red curve) is also reactive, displaying a few proton transfers and producing transient charged species. This latter model is the most energetically stable model (see Table 1).

We now turn to the impact of these chemical events on the charge of the GO layer. To do so, looking at the presence of alcoholate groups is not sufficient, as they can arise from both deprotonation and epoxide opening (which does not induce a net charge on the GO). Looking at the number of hydronium cations, which is by design the opposite of the net GO charge, the results are shown in Table 4. The average charge per replica for the semiaromatic models is negative, and its absolute value spreads from 2.5 to 13 mC m⁻²—with an average value of 6 mC m⁻². These values, obtained from our ab initio simulations, compare favorably with recent measurements stating 16 mC m⁻²²⁴⁶ or with values extracted from zeta potential measurements at neutral pH being 11 mC m⁻² for GO and 3 mC m⁻² for reduced GO (rGO)⁴⁷. It should be stressed that our GO replicas resemble more their rGO samples, as they contain no carboxylic functions. Moreover, random models yield different negative charges but remaining in the same order of magnitude (see Table S1).

It is an important feature that a GO sheet should no longer be modeled neutral when dipped into water, but instead carrying a modest negative charge, ranging from a few to tens of mC m⁻² and due to by partial deprotonation of its hydroxyl groups. It confirms that the pKₐ of GO is lower than the value of alcohols (16–19), and presumably lower than 14. This surface charge may explain electrostatic repulsion between flakes that induces their perfect layering at confined distances.
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Dynamics of graphene oxide in liquid water. Time evolution of the number of neutral functional groups of the GO surface (left) and charged species present at the GO surface or in the water bulk (right). **a** -OH groups. **b** Epoxides. **c** Surface alcoholate. **d** $H_3O^+$. The black, red, and blue curves are associated to the semiordered 1, 3, and 4 GO models, respectively. Some reactive events occurring in the GO dynamics in water are also identified and marked with capital letters.</figcaption></figure>

<table><caption>**Table 4** Excess charge/GO surface (in $mC m^{-2}$) on the GO interfacing bulk water (equivalent with opposite sign to the excess charge distributed over the $H_2O$ molecules of the bulk).</caption><thead><tr><td>GO configuration</td><td>Excess surface charge (mC m<sup>-2</sup>)</td></tr></thead><tbody><tr><td>Semiordered 1</td><td>-2.46</td></tr><tr><td>Semiordered 3</td><td>-4.27</td></tr><tr><td>Semiordered 4</td><td>-12.96</td></tr><tr><td>Average</td><td>-6.56</td></tr></tbody></table>

## Discussion

In summary, by means of extensive ab AIMD simulations, we have built several realistic models of GO at a fixed O/C ratio, comprising hydroxyls and epoxide groups, and analyzed their respective stability, as well as physical and chemical behavior in liquid water. Our computational study emphasizes the importance of partly ordered models of GO in order to tackle its physical and chemical properties. Although computationally expensive, it is the combination of such models with the use of extensive ab initio MD that allows us to describe the chemical processes occurring at the GO/liquid water interface. In particular, we demonstrate the impact of functional group distribution on the GO sheet, and show that semiordered models—with correlated functional groups and some regions of pristine graphene—are the most stable structures in vacuum as well in liquid water. We also demonstrate the formation of a net negative charge on the GO induced by the reactivity with water, with a moderate charge in the same range of silica materials in water, but two orders of magnitude less than hexagonal boron nitride layers$^{48,49}$. Altogether this combination may help understand the

quite unique properties of GO in terms of water and ion transport, which are at the root of its use in water filtration and remediation. Our results highlight fast water dynamics, impeded neither by the confinement nor by the presence of hydrophilic oxygen groups. Notably, the favorable clustering of oxygen functions may open fast dynamic pathways for water transport on the remaining pristine graphene regions. This could explain the large hydrodynamic slip-page of water across GO, as reported experimentally$^{37}$.

Because GO is experimentally very diverse in its chemical composition, and in the interlayer spacing of hydrated GO, a lot of perspectives emerge in ab initio modeling of GO in different conditions. In particular, the impact of the confinement thickness on the water dynamics is of interest for the future. The present work also opens the way for further studies on the impact of oxygen concentration (number of functional groups) on the GO/water interactions, by varying the O/C ratio and introducing carboxyl groups. The dynamics of charged ions at this charged interface, including hydronium and hydroxide ions, is also a wide open question, that quantum modeling at the electronic scale should be able to tackle.

## Methods

**Density functional theory (DFT) calculations.** All the simulations discussed in the paper and in this section have been performed at the ab initio level, using the DFT approach. DFT has therefore been used to both optimize the materials' geometry at zero temperature and to evaluate the ionic forces acting on each atom along the AIMD trajectories. Most of the calculations have been carried out with the CP2K software$^{50}$ and some preliminary electronic structure relaxations have been performed with VASP$^{51}$. With C2PK, the DFT method is straightforwardly employed using the Quickstep module$^{52,53}$ of the package. As in ref. $^{54}$, DZVP-MOLOPT-SR-GTH basis sets were used$^{55}$ along with planewaves expanded to a 800 Ry absolute energy cutoff and a 40 Ry relative cutoff. Electronic cores were represented by Geodecker–Teter–Hutter pseudopotentials$^{56–58}$. The Perdew, Burke, and Ernzerhof functional was used$^{59}$ with the DFT-D3 dispersion
correction scheme$^{60,61}$. Technical details about the convergence of the calculation settings are provided in the Supplementary Information Section III.1.

**Construction of graphene oxide models.** In detail, we first built each empty semiarid or random GO replica by adding chemical functions according to the constraints detailed in the first part of the paper. Then, these structures have been placed into a 15 × 13 × 16 Å³ orthorhombic cell, corresponding to a 6 × 6 graphene-like lattice with an interlayer spacing of 16 Å. Then, we performed a preliminary geometry optimization and a cell relaxation using VASP. Indeed, adding chemical functions on the GO surface induces a mechanical strain and the lattice parameters *a* and *b* need to be adjusted. This procedure is then repeated with the C2PK code to ensure that the minimum energy configuration has been reached and to check the consistency of the results with two different codes. Once the empty surface is fully relaxed (i.e., the pressure of the system is about 10–100 MPa, a standard value for such materials), we generated a single 30 ps AIMD trajectory for each replica at *T* = 300 K, using periodic boundary conditions. The deuterium mass is substituted for all protons to reduce the timestep size needed for energy conservation in our Born Oppenheimer dynamics AIMD and to limit nuclear quantum effects of the proton. All the AIMD simulations were carried out with a 0.5 fs timestep in the NVT ensemble, using a stochastic velocity rescaling (SVR) thermostat$^{62}$ with a time constant of 1 ps.

**Construction of hydrated graphene oxide models.** After having ensured the stability and studied the structural properties of the anhydrous GO structures, we embedded them into explicit water. This is achieved by adding 80 H₂O molecules using the PACKMOL program$^{63}$. We imposed the water molecules to be distant of at least 2 Å from the GO surface to avoid steric overload due to the presence of hydroxyl/epoxide functions on the surface. The procedure is then very similar to the anhydrous case. First, a geometry and a simulation cell relaxation is necessary to make the H₂O molecules organized in a relevant chemical configuration. Let us notice that only the *c* cell parameter along the *z*-axis is allowed to change to preserve the GO structure in the *xy* plane. Later, a short AIMD trajectory of 1 ps at *T* = 300 K is with a small time constant of 100 fs for the SVR thermostat is generated to make the water liquid around the GO surface. Afterwards, a new cell relaxation is applied and a *Δc* = 1.5–1.8 Å variation is typically observed during this step. Finally, AIMD trajectories are generated after a 5 ps equilibration period: 60 ps long for each semiarid replica while we stopped the calculation at 30 ps for the random structures. A 60-ps-long trajectory of bulk liquid water at the experimental density *ρ* = 0.99657 g cm⁻³ is also generated to discuss the diffusion coefficients on similar trajectories (see Supplementary Information Section II.9 for further details).

Finally, the observables presented in the paper are obtained using block averages of the AIMD trajectories of each material, using Chemfiles/cfiles$^{39,40}$.

**Data availability**

Data supporting this study is available online on the repository at https://github.com/fxcoudert/citable-data. They can also be made available upon request to the corresponding author.

Received: 22 November 2019; Accepted: 4 March 2020;
Published online: 26 March 2020

**References**

1. Dreyer, D. R., Park, S., Bielawski, C. W. & Ruoff, R. S. The chemistry of graphene oxide. *Chem. Soc. Rev.* **39**, 228–240 (2010).

2. Zhu, Y. et al. Graphene and graphene oxide: synthesis, properties, and applications. *Adv. Mater.* **22**, 3906–3924 (2010).

3. Joshi, R. K. et al. Precise and ultrafast molecular sieving through graphene oxide membranes. *Science* **343**, 752–754 (2014).

4. Su, Y. et al. Impermeable barrier films and protective coatings based on reduced graphene oxide. *Nat. Commun.* **5**, 3906 (2014).

5. Kumar, S. et al. Graphene oxide–MnFe₂O₄ magnetic nanohybrids for efficient removal of lead and arsenic from water. *ACS Appl. Mater. Interfaces* **6**, 17426–17436 (2014).

6. Favaro, M. et al. Multiple doping of graphene oxide foams and quantum dots: new switchable systems for oxygen reduction and water remediation. *J. Mater. Chem. A* **3**, 14334–14347 (2015).

7. Dutta, D., Thiyagarajan, S. & Bahadur, D. SnO₂ quantum dots decorated with graphene oxide nanocomposites for efficient water remediation. *Chem. Eng. J.* **297**, 55–65 (2016).

8. Yadav, R., Subhash, A., Chemmenchery, N. & Kandasubramanian, B. Graphene and graphene oxide for fuel cell technology. *Ind. Eng. Chem. Res.* **57**, 9333–9350 (2018).

9. Farooqui, U., Ahmad, A. & Hamid, N. Graphene oxide: a promising membrane material for fuel cells. *Renew. Sustain. Energy Rev.* **82**, 714–733 (2018).

10. Nakajima, T., Mabuchi, A. & Hagiwara, R. A new structure model of graphite oxide. *Carbon* **26**, 357–361 (1988).

11. Szabó, T. et al. Evolution of surface functional groups in a series of progressively oxidized graphite oxides. *Chem. Mater.* **18**, 2740–2749 (2006).

12. Stankovich, S., Piner, R. D., Nguyen, S. T. & Ruoff, R. S. Synthesis and exfoliation of isocyanate-treated graphene oxide nanoplatelets. *Carbon* **44**, 3342–3347 (2006).

13. Stankovich, S. et al. Synthesis of graphene-based nanosheets via chemical reduction of exfoliated graphite oxide. *Carbon* **45**, 1558–1565 (2007).

14. Yang, D. et al. Chemical analysis of graphene oxide films after heat and chemical treatments by X-ray photoelectron and micro-raman spectroscopy. *Carbon* **47**, 145–152 (2009).

15. Schniepp, H. C. et al. Functionalized single graphene sheets derived from splitting graphite oxide. *J. Phys. Chem. B* **110**, 8535–8539 (2006).

16. McAllister, M. J. et al. Single sheet functionalized graphene by oxidation and thermal expansion of graphite. *Chem. Mater.* **19**, 4396–4404 (2007).

17. He, H., Klinowski, J., Forster, M. & Lerf, A. A new structural model for graphite oxide. *Chem. Phys. Lett.* **287**, 53–56 (1998).

18. Lerf, A., He, H., Forster, M. & Klinowski, J. Structure of graphite oxide revisited. *J. Phys. Chem. B* **102**, 4477–4482 (1998).

19. Yan, J.-A., Xian, L. & Chou, M. Y. Structural and electronic properties of oxidized graphene. *Phys. Rev. Lett.* **103**, 086802 (2009).

20. Lee, D. W. et al. The structure of graphite oxide: Investigation of its surface chemical groups. *J. Phys. Chem. B* **114**, 5723–5728 (2010).

21. Bagri, A. et al. Structural evolution during the reduction of chemically derived graphene oxide. *Nat. Chem.* **2**, 581–587 (2010).

22. Wei, N., Peng, X. & Xu, Z. Understanding water permeation in graphene oxide membranes. *ACS Appl. Mater. Interfaces* **6**, 5877–5883 (2014).

23. Devanathan, R., Chase-Woods, D., Shin, Y. & Gotthold, D. W. Molecular dynamics simulations reveal that water diffusion between graphene oxide layers is slow. *Sci. Rep.* **6**, 457 (2016).

24. Cohen-Tanugi, D., Lin, L.-C. & Grossman, J. C. Multilayer nanoporous graphene membranes for water desalination. *Nano Lett.* **16**, 1027–1033 (2016).

25. Hou, D. & Yang, T. A reactive molecular dynamics study of graphene oxide sheets in different saturated states: structure, reactivity and mechanical properties. *Phys. Chem. Chem. Phys.* **20**, 11053–11066 (2018).

26. Cassagneau, T., Guérin, F. & Fendler, J. H. Preparation and characterization of ultrathin films layer-by-layer self-assembled from graphite oxide nanoplatelets and polymers. *Langmuir* **16**, 7318–7324 (2000).

27. Boukhvalov, D. W. & Katsnelson, M. I. Modeling of graphite oxide. *J. Am. Chem. Soc.* **130**, 10697–10701 (2008).

28. Yang, Q. et al. Ultrathin graphene-based membrane with precise molecular sieving and ultrafast solvent permeation. *Nat. Mater.* **16**, 1198–1202 (2017).

29. Zhou, K.-G. et al. Electrically controlled water permeation through graphene oxide membranes. *Nature* **559**, 236–240 (2018).

30. Kudin, K. N. et al. Raman spectra of graphite oxide and functionalized graphene sheets. *Nano Lett.* **8**, 36–41 (2008).

31. Yan, J.-A. & Chou, M. Y. Oxidation functional groups on graphene: structural and electronic properties. *Phys. Rev. B* **82**, 17 (2010).

32. Wang, B., Bocquet, M.-L., Marchini, S., Günther, S. & Wintterlin, J. Chemical origin of a graphene Moiré overlayer on Ru(0001). *Phys. Chem. Chem. Phys.* **10**, 3530–3534 (2008).

33. Wang, B., Bocquet, M.-L., Günther, S. & Wintterlin, J. Comment on “Periodically rippled graphene: growth and spatially resolved electronic structure”. *Phys. Rev. Lett.* **101**, 099703 (2008).

34. Zhou, S. & Bongiorno, A. Origin of the chemical and kinetic stability of graphene oxide. *Sci. Rep.* **3**, 457 (2013).

35. Kumar, P. V. et al. New insights into the thermal reduction of graphene oxide: impact of oxygen clustering. *Carbon* **100**, 90–98 (2016).

36. Xia, S., Ni, M., Zhu, T., Zhao, Y. & Li, N. Ultrathin graphene oxide nanosheet membranes with various *d*-spacing assembled using the pressure-assisted filtration method for removing natural organic matter. *Desalination* **371**, 78–87 (2015).

37. Abraham, J. et al. Tunable sieving of ions using graphene oxide membranes. *Nat. Nanotechnol.* **12**, 546–550 (2017).

38. Khoei, A. R. & Khorrami, M. S. Mechanical properties of graphene oxide: a molecular dynamics study. *Fuller. Nanotub. Carbon Nanostruct.* **24**, 594–603 (2016).

39. Fraux, G. Chemfiles, available online at https://zenodo.org/record/3490304#.Xab2mS3M1QI.

40. Fraux, G. cfiles, available online at https://github.com/chemfiles/cfiles.

41. Jeffrey, G. A. *An Introduction to Hydrogen Bonding* (Oxford University Press, Oxford, 1997).

42. Coudert, F.-X., Cailliez, F., Vuilleumier, R., Fuchs, A. H. & Boutin, A. Water nanodroplets confined in zeolite pores. *Faraday Discuss.* **141**, 377–398 (2009).

43. Fogarty, A. C., Coudert, F.-X., Boutin, A. & Laage, D. Reorientational dynamics of water confined in zeolites. *ChemPhysChem* **15**, 521–529 (2014).
44. Savazzi, F., Risplendi, F., Mallia, G., Harrison, N. M. & Cicero, G. Unravelling some of the structure–property relationships in graphene oxide at low degree of oxidation. *J. Phys. Chem. Lett.* **9**, 1746–1749 (2018).

45. Agmon, N. The Grotthuss mechanism. *Chem. Phys. Lett.* **244**, 456–462 (1995).

46. Nair, R. Private communication.

47. Konkena, B. & Vasudevan, S. Understanding aqueous dispersibility of graphene oxide and reduced graphene oxide through pKₐ measurements. *J. Phys. Chem. Lett.* **3**, 867–872 (2012).

48. Siria, A. et al. Giant osmotic energy conversion measured in a single transmembrane boron nitride nanotube. *Nature* **494**, 455–458 (2013).

49. Grosjean, B. et al. Chemisorption of hydroxide on 2D materials from DFT calculations: graphene versus hexagonal boron nitride. *J. Phys. Chem. Lett.* **7**, 4695–4700 (2016).

50. Hutter, J., Iannuzzi, M., Schiffmann, F. & VandeVondele, J. CP2K: atomistic simulations of condensed matter systems. *WIREs Comput. Mol. Sci.* **4**, 15–25 (2014).

51. Hafner, J. Ab-initio simulations of materials using VASP: density-functional theory and beyond. *J. Comput. Chem.* **29**, 2044–2078 (2008).

52. Vandevondele, J. et al. Quickstep: fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. *Comput. Phys. Commun.* **167**, 103–128 (2005).

53. Lippert, G., Hutter, J. & Parrinello, M. A hybrid gaussian and plane wave density functional scheme. *Mol. Phys.* **92**, 477–488 (1997).

54. Grosjean, B., Bocquet, M.-L. & Vuilleumier, R. Versatile electrification of two- dimensional nanomaterials in water. *Nat. Commun.* **10**, 2541 (2019).

55. VandeVondele, J. & Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. *J. Chem. Phys.* **127**, 114105 (2007).

56. Goedecker, S., Teter, M. & Hutter, J. Separable dual-space Gaussian pseudopotentials. *Phys. Rev. B* **54**, 1703–1710 (1996).

57. Hartwigsen, C., Goedecker, S. & Hutter, J. Relativistic separable dual-space Gaussian pseudopotentials from H to Rn. *Phys. Rev. B* **58**, 3641–3662 (1998).

58. Krack, M. Pseudopotentials for H to Kr optimized for gradient-corrected exchange-correlation functionals. *Theor. Chem. Acc.* **114**, 145–152 (2005).

59. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. *Phys. Rev. Lett.* **77**, 3865–3868 (1996).

60. Grimme, S., Antony, J., Ehrlich, S. & Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. *J. Chem. Phys.* **132**, 154104 (2010).

61. Grimme, S., Ehrlich, S. & Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. *J. Comput. Chem.* **32**, 1456–1465 (2011).

62. Bussi, G., Donadio, D. & Parrinello, M. Canonical sampling through velocity rescaling. *J. Chem. Phys.* **126**, 014101 (2007).

63. Martínez, L., Andrade, R., Birgin, E. G. & Martínez, J. M. Packmol: a package for building initial configurations for molecular dynamics simulations. *J. Comput. Chem.* **30**, 2157–2164 (2009).

## Acknowledgements

F.M. and M.-L.B. acknowledge funding from the H2020-424 FETOPEN project NANOPHLOW number 766972. This work was granted access to the HPC resources of TGCC and CINES under grants A5-A0070807364 and A0070807069 by GENCI. We thank Vasu Kalangi and Lydéric Bocquet for insightful discussions about GO membranes, and Guillaume Fraux for his help with Chemfiles/cfiles.

## Author contributions

M.-L.B. designed the research. F.M. performed the simulations. F.M., F.-X.C., and M.-L.B. analyzed the data, wrote and revised the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41467-020-15381-y.

Correspondence and requests for materials should be addressed to M.-L.B.

Peer review information *Nature Communications* thanks the anonymous reviewers for their contribution to the peer review of this work. Peer reviewer reports are available.

Reprints and permission information is available at http://www.nature.com/reprints

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020