Structure of water clusters on graphene: A classical molecular dynamics approach

The microscopic structure of surface water adsorbed on graphene is elucidated theoretically by classical molecular dynamics simulation. At a low temperature (100 K), the main polygon consisting of hydrogen bonds in single-layered water on graphene is tetragonal, whereas the dominant polygons in double-layered water are tetragonal, pentagonal, and hexagonal. On the other hand, at room temperature, the tetragonal, pentagonal, and hexagonal water clusters are the main structures in both single- and double-layered water.


Introduction
The microscopic structure of surface water on solid surfaces is tremendously important in understanding numerous physical and chemical characteristics of surfaces. As one example, the structure of surface water influences the hydrophilicity, and hence the wettability, of the surface. 1,2) In addition, understanding the properties of surface water plays an essential role in the development and improvement of various devices, such as field-effect transistors 3) and gas sensors, 4,5) in the fields of electrical engineering and biology. Many studies on surface water have been conducted to date. [6][7][8][9][10][11][12] Some of these studies focused on the structure of water on various oxide substrates in experiments 6,7) and simulations, 8,9) revealing the layered structure of water in the normal direction to the surface, while others focused on metal substrates, [10][11][12] revealing a similarly layered structure. Although the structure of surface water along the surface normal has been systematically and rigorously investigated, the microscopic structure of surface water parallel to a solid surface has not yet been fully characterized despite its importance.
Depending on the environmental conditions, water molecules form a variety of global hydrogen bond networks consisting of various cluster structures, [13][14][15] which contribute to the unique characteristics of water. Water clusters have mainly been investigated by numerical simulations, and a number of stable clustering structures have been proposed. In the case of bulk (H 2 O) n complexation, numerical calculations based on ab initio simulations have suggested the formation of trigonal (n = 3), tetragonal (n = 4), pentagonal (n = 5), and hexagonal (n = 6) clusters. 13) Water clusters are expected to occur even in the case of surface water. Investigating cluster structures in surface water is essential for understanding the physical and chemical properties of solid=water interfaces.
In this study, we focus on graphene, which is used in the fields of semiconductor engineering, 16,17) optics, 18) and biotechnology. 19) Because graphene is an atomically flat and nonpolar material consisting of only carbon atoms, it is an ideal system for exploring the microscopic structure of surface water parallel to a solid surface. In addition, it has long been believed that water molecules cannot adsorb on graphene because it is a hydrophobic material. 20) However, the existence of water layers on a graphene surface was theoretically confirmed by classical molecular dynamics (MD) simulations. 21,22) These previous works suggested that water molecules spread out on the graphene surface. Before these observations and simulations, Homma et al. discovered that similar water molecular layers can exist around a carbon nanotube (CNT) by both photoluminescence measurements and classical MD simulations. 23) In addition, the structure of water clusters on graphene has been investigated by ab initio calculations on the basis of density functional theory (DFT), 24,25) and some stable water cluster structures have been proposed. However, as only one water cluster was considered in these DFT simulations, 24,25) the interactions among water clusters remain poorly understood, even though they constitute essential information regarding the global hydrogen network. Thus, it is clearly important to investigate the formation of surface water clusters on a graphene substrate.
In this study, we theoretically investigated the behavior of water molecules on graphene substrates by classical MD simulations, and analyzed the formation of surface water clusters, focusing on the two-dimensional water clusters at the interface. We discuss three models, in which the number of water molecules was varied. The statistical analysis of water clusters is difficult because of their short lifetime at room temperature; 26,27) however, here we performed longterm simulations of each model and succeeded in statistically analyzing the occurrence rate of various types of surface water clusters on graphene using the above-mentioned three models at both room temperature and a low temperature.

Methodology
The behavior of water molecules on graphene was examined by MS ForcitePlus in Materials Studio 8.0. 28) COMPASS II was employed in all calculations. 29,30) COMPASS is one of the universal force fields based on ab initio simulations. The force field for water in COMPASS II is a flexible model incorporating the inter-and intramolecular motion of a water molecule, which is appropriate for investigating the formation of water molecules microscopically. The cutoff radius of short-range interactions was set to 1.0 nm, and the long-range Coulomb interactions were calculated by the Ewald method. 31) The charge distribution in a water molecule was set to −0.41 eC for hydrogen and 0.82 eC for oxygen, which were determined using COMPASS II. The stable structure of water molecules as determined by COMPASS II is displayed in Fig. 1(a), where the bond angle and bond distance are estimated by averaging the results of the above bulk water simulations. The carbon atoms constituting graphene are charge-neutral. The classical equation of motion was solved numerically by the velocity Verlet algorithm with a time step of 1.0 fs. The temperature was kept at 250 K for room temperature and 100 K for the low temperature using a Nosé thermostat. 32) It is known that COMPASS II overestimates the diffusion constant and density of water molecules. Indeed, we confirmed that water at 250 K had similar diffusion constant and density to that at actual room temperature (Appendix A). The simulations at 100 K were performed for comparison with the formation of water clusters at 250 K. Figure 1(b) shows an overall view of the simulation cell on which the periodic boundary conditions have been imposed. Graphene is set at 0.34 nm from the bottom of the cell. N water molecules (N = 60, 100, or 160) were randomly placed immediately above the graphene in the initial state. In addition, a large simulation cell is prepared to calculate the domain size of the water cluster structure. Each large simulation cell contains M water molecules (M = 150, 250, or 400, respectively). Simulations were performed for 100 ns for the small cells to analyze the hysteresis of the water cluster structure in order to ensure statistical reliability. On the other hand, simulations were performed for 20 ns for the large cells to visualize the detailed water cluster structure. All calculations were executed after 5.0 ns of geometric optimization. The detailed properties of the simulation cells are listed in Table I. 3. Results and discussion 3.1 Vertical distribution of surface water on graphene Figure 2 shows the distributions of the oxygen atoms of water molecules along the normal to the graphene surface (z-axis). The graphene is located at z = 0. As shown in Fig. 2(a), a single-layered structure is formed at both 100 and 250 K for a simulation cell including 60 H 2 O molecules (called the 60-H 2 O cell). The water peak at 100 K is sharper than that at 250 K. There is more water evaporation on the graphene surface at 250 K than there is at 100 K. On the other hand, as shown in Figs    formed on graphene 22) and around a CNT. 23) For the 160-H 2 O cells at 250 K, the water above the second water layer is regarded as bulk water, where water molecules move relatively freely, whereas at 100 K, the water molecules above the second layer exhibit distributional peaks. The surface layered structure influences regions further from the surface when forming ice at a low temperature.  22,23) Again, we note that the water molecules can adsorb on the graphene surface as a result of the van der Waals interaction. Here, we define the water molecules located between the graphene surface and 0.5 nm from the surface as the first layer, and those between 0.5 and 0.7 nm from the surface as the second layer. The water densities of the first layer are 12.6 nm −2 for the 1L model at 100 K, 11.3 nm −2 for the 2L and IB models at 100 K, and 10.7 nm −2 for all boxes at 250 K. As explained later in detail, these differences in density arise from the formation of water clusters. Figure 3 demonstrates the angular distribution of OH bonds of water molecules in the water layers on graphene. The cosine (cos θ) values are calculated from the angle between the normal vector to the graphene surface and the OH vector from the oxygen atom to the hydrogen atom in the water molecules. When the OH bonds are oriented parallel to the graphene surface (θ = π=2), cos θ is equal to zero. A positive cos θ value means that the OH vector points away from the graphene surface (hereafter, H-up orientation), while a negative cos θ value indicates that the OH vector points toward the graphene surface (H-down orientation). In the first layer of the 1L model, most OH bonds are oriented parallel to the graphene surface (cos θ = 0) as seen in Fig. 3(b). The H-up and H-down orientations are rarely observed. This means that almost all water molecules contribute to the intralayer hydrogen bonds, i.e., the hydrogen bonds are "closed" in the first layer. This feature is common to both the 100 K (blue curves) and 250 K (orange curves) cases, although the probability distribution at 250 K is broader than that at 100 K. By contrast, in the first layer of the 2L model [ Fig. 3(c)] and that of the IB model [ Fig. 3(d)], a peak near cos θ = 0 still exists, but the frequency around cos θ = 1 increases, which means that the number of OH bonds with H-up orientation increases. This is because the water molecules in the first layer connect to those in the second layer via hydrogen bonds. When there are more than two layers on the graphene, both intra-and interlayer hydrogen bonds are formed. As a result, the peak at cos θ ≈ −1 increases for the second layer, which again indicates that the water molecules in the second layer connect to those in the first layer via hydrogen bonds. The peaks near cos θ = 0 for the first layer in the 2L and IB models shift in the negative direction relative to the peak for the 1L model. The bond angle of a water molecule (the H-O-H angle) is not a right angle (90°) but approximately 109°in COMPASS II, as calculated in the water bulk. 33) Consequently, when one OH bond in a given water molecule is oriented upward, the other OH bond in that water molecule is oriented slightly downward. The same mechanism can explain why the peaks near cos θ = 0 for the second layer in the 2L and IB models shift in the positive direction relative to the peak in the 1L model. These orientation trends are observed at both 100 and 250 K, except that all peaks are broader at 205 K owing to thermal fluctuations. the hydrogen bond network (shown with color coding on the right). The black lines in the middle and right panels indicate the hydrogen bonds; if the O-O distance between two water molecules is less than 0.36 nm, the O-H distance is less than 0.245 nm and the H-O£O angle is less than 30°. 34,35) As seen from the snapshot in the left panel in Fig. 4(a), the water molecules fill the layer in an orderly arrangement. Each water molecule strongly interacts with other water molecules via hydrogen bonds as shown in the middle panel in Fig. 4(a). We can see from the right panel in Fig. 4(a) that the water cluster structure mainly consists of tetragons (shown in blue). In the 1L model, tetragonal clusters are the most efficient and stable to cover the surface because all OH bonds can participate in the hydrogen bonds in two dimensions. The other polygonal clusters, such as pentagonal (green) and hexagonal (pink) clusters, are less efficient. If the water layer is filled with these nontetragonal clusters, at least one OH bond must be H-up or H-down way to cover the surface, and this OH bond cannot connect with other water molecules by hydrogen bonds. The tetragonal clusters are linked to other tetragonal clusters and form a kind of "water crystal" domain. Pseudocrystal domains consisting of tetragonal clusters are interspersed on graphene, consequently; the mismatch between neighboring domains appears as a grain boundary. Along this boundary, water molecules form nontetragonal clusters, such as trigonal, tetragonal, and hexagonal clusters. The formation of water clusters along the boundary is found to be consistent with the formation of domain discordance. Figure 4(b) illustrates snapshots of the water molecules in the first layer for the 2L model at 100 K (left panel) and of the corresponding hydrogen bonds (middle panel). As seen in the left panel, the water molecules fill the layer, but they seem less densely packed than in the 1L model. As in the case of the 1L model, trigonal, tetragonal, pentagonal, and hexagonal water clusters are formed [see the right panel in Fig. 4(b)]. However, the 2L model exhibits more pentagonal and hexagonal clusters and fewer tetragonal clusters than the 1L model. These differences are caused by the presence of the second layer as explained in the following. Unlike the simulations of the 1L model, the H-up-oriented OH bonds derived from pentagonal and hexagonal clusters in the first layer can connect to the oxygen atoms in water in the second layer by hydrogen bonds. The interaction between the upper and lower layers increases the stability of the pentagonal and hexagonal clusters. The same trends are observed for the IB model at 100 K [ Fig. 4(c)].

Cluster structure of water molecules
The schematics in Figs. 4(d) and 4(e) pertain to the second layers in the 2L and IB models at 100 K, respectively. The snapshots show the dense and sparse areas of water molecules. The influence of graphene on the water molecules decreases with increasing distance from the graphene surface; consequently, the density of the second layer (9.07 nm −2 ) is lower than that of the first layer (11.3 nm −2 ), as estimated from Fig. 2(b). This results in the distribution of dense and sparse areas. In addition, it can be seen that there are few clusters in the second layers, and these tend to have tetragonal, pentagonal, and hexagonal structures. Figure 5 shows some of the ordinary two-dimensional water cluster structures observed in the present simulations. The trigonal cluster is illustrated in Fig. 5(a). It is the only case of a (H 2 O) 3 cluster. Figures 5(b) and 5(c) show the tetragonal clusters. In the cluster in Fig. 5(b), all oxygen and hydrogen atoms are located in the same plane. This structure usually appears in the simulation of the 1L model at 100 K, because the OH bonds form part of the hydrogen bonds with neighboring water clusters. In the cluster in Fig. 5(c), one hydrogen atom is oriented H-up or H-down. This structure frequently occurs in the 250 K simulations, because the water cluster becomes isolated owing to the absence of a crystal structure. Figure 5   six clusters were reported in a previous theoretical study that investigated bulk (H 2 O) n clusters. 13) In this paper, we suggest that such two-dimensional clusters are formed and take on the main structure immediately above the graphene. This diversity of (H 2 O) 6 clusters explains why there are only a small number of hexagonal clusters observed under all simulation conditions. (As shown in Sect. 3.4, the proportion of hexagonal clusters is less than 20% under all conditions.) Some illustrations of hydrogen bonds for the first layer of the 1L model at 250 K are shown in Fig. B·1 of Appendix B. The lifetimes of hydrogen bonds at room temperature corresponding to an MD temperature of 250 K are very much shorter than those at lower temperatures. In addition, the water molecules at 250 K move toward and away from the liquid= vapor interface more frequently than those at 100 K. In order to determine the difference between the water cluster structures at 100 and 250 K, we counted the number of water clusters in the 500 snapshots extracted from each simulation and integrated the ratio of each snapshot. Snapshots were taken once every 200 ps during each simulation. Figure 6 shows the numbers of various water clusters that appeared during the calculations and the proportions of the water clusters formed, e.g., trigonal, tetragonal, pentagonal, and hexagonal clusters in the first layer, for each set of simulation conditions, viz., the 1L and 2L models at 100 K [Figs. 6(a) and 6(b), respectively] and the 1L and 2L models at 250 K [Figs. 6(c) and 6(d), respectively].

Relative proportions of various polygonal cluster structures
As seen in Fig. 6(a), the tetragonal cluster is the main structure of the 1L model at 100 K with 68.72% of the total number of clusters. The second most frequent is the pentagonal cluster (17.92%). On the other hand, as seen in Fig. 6(b), the percentages of the tetragonal (31.14%) and pentagonal clusters (29.83%) are comparable for the 2L model at 100 K. The hexagonal structure is also more common than in the 1L model. As mentioned above, the interactions between the water molecules in the first and second layers give the pentagonal and hexagonal clusters stability. The variety of water clusters increases the  domain discordance, which in turn increases the occurrence of other clusters with heptagonal, octagonal, and higher-order polygonal structures. The heptagonal and octagonal clusters each constitute roughly 7% of the total (included in the clusters labeled "others"). Most of the heptagonal and octagonal clusters appear to be combinations of defective trigonal, pentagonal, and hexagonal clusters that do not form complete loops (Fig. C·1 in Appendix C). This was expected from the DFT calculation results, which suggested that water clusters involving at least seven H 2 O molecules form three-dimensional cluster structures. 10,36,37) Tetragonal clusters are more common than pentagonal and hexagonal clusters; thus, such differences in the microscopic structure of the water clusters leading to the difference in the water density of the first layer have already been mentioned in Sect. 3.1.
As shown in Fig. 6(c), for the 1L model at 250 K, the tetragonal cluster is the most frequent, as in the case of 100 K; however, it accounts for no more than 39.78%. Interestingly, the relative proportions of the water clusters in the 1L model are very similar to those in the first layer of the 2L model at 250 K. The reasons for this are as follows. First, the watercrystal-like structure of the 1L model is destroyed because the frequency with which water molecules move toward and away from the liquid=vapor interface increases with temperature. Thus, the advantage of the tetragonal formation, in which all OH bonds form part of the hydrogen bonds, decreases as the temperature increases. Second, the proportions of H-up and H-down orientations in the 1L model at 250 K are higher than those at 100 K (Fig. 3) because of the presence of liquid= vapor interfaces. Indeed, some OH bonds in the liquid=vapor interfaces are oriented H-up or H-down according to a previous study. 38) Therefore, at 250 K, the relative proportions of the water clusters remain unchanged regardless of the number of water molecules. The present results suggest that the water clusters at the interface typically have trigonal, tetragonal, pentagonal, and hexagonal structures.

Conclusions
We investigated the surface water cluster structure on a graphene surface by classical MD simulations. Even when the number of water molecules on graphene is small, they cover the entire graphene surface. We prepared only a single-layer (1L) model, a double-layer (2L) model, and an including bulk (IB) model, and performed simulations with each model at 100 and 250 K, the latter temperature taken as room temperature.
In the case of the 1L model at 100 K, most of the OH bonds are oriented parallel to graphene and connect to the other water molecules in the first layer by hydrogen bonds. Water mainly forms tetragonal clusters, because all OH bonds can form hydrogen bonds only in the case of tetragonal clusters in the 1L model. In particular, tetragonal clusters connect with each other, constructing a crystal-like structure at 100 K. Pentagonal and hexagonal clusters are located at the boundary of this water crystal domain.
In comparison with the 1L model, in the 2L model at 100 K, there is an increased number of pentagonal and hexagonal clusters in the first layer. The OH bonds with H-up orientation in the first layer can connect with the oxygen atoms in the second layer, which stabilize the pentagonal and hexagonal clusters. By contrast, in the second layer, some OH bonds adopt a H-down orientation in order to form hydrogen bonds with the water in the first layer. Water clusters rarely occur, but tetragonal, pentagonal, and hexagonal clusters can be observed. At 250 K, the relative proportions of the water cluster structures in the 1L and 2L models are comparable. This is because at 250 K, the crystal-like structure collapses and some OH bonds are oriented H-up or H-down owing to the liquid= vapor interface. The theoretical calculations in this study demonstrate that the surface water clusters mainly consist of trigonal, tetragonal, pentagonal, and hexagonal clusters.
In this Appendix, we discuss the temperature dependences of the diffusion constant and density of bulk water molecules. There are 512 H 2 O molecules in the 1.971 nm cubic box with periodic boundary conditions. The temperature is maintained at 200, 240, 250, 260, or 300 K by a Nosé thermostat. The pressure is maintained at 1 atm by a Berendsen barostat. The other simulation conditions follow those in the main paper. The density is estimated by taking the running average of each simulation. The diffusion constant is calculated using Einstein's diffusion equation. As shown in Fig. A·1, we can see that water at 250 K in the present MD simulation has similar diffusion constant and density to that at actual room temperature. 39,40) Appendix B: Dynamics of hydrogen bonds in the first layer in the 1L model at 250 K Figure B·1 shows snapshots of the time evolutions of hydrogen bonds for the first layer in the 1L model at 250 K.
Appendix C: Heptagonal and octagonal clusters on graphene is too long to make a hydrogen bond. Figure C·1(b) shows a heptagonal cluster consisting of defective trigonal and hexagonal clusters. The 2-7 bond is too long. Figure C·1(c) shows an octagonal cluster consisting of four defective tetragonal clusters. The 2-7 and 3-6 bonds are too long. Figure C·1(d) shows an octagonal cluster consisting of two defective pentagonal clusters. The 4-8 bond is too long. Figure C·1(e) shows an octagonal cluster consisting of two defective pentagonal clusters.