Ab initio Molecular Dynamics Study of Adsorption of Hydroxyl Groups on Graphene Surface†

Reduced graphene oxide is the precursor to produce graphene in a large scale; however, to date, there has been no consensus on the electronic structure of reduced graphene oxide. In this study, we carried out an ab initio molecular dynamics simulation to investigate the adsorption process of hydroxyl groups on graphene surface. During the adsorption process, the OH group needs to firstly pass through a physical adsorption complex with the OH above the bridge site of two carbon atoms, next to surmount a transition state, then to be adsorbed at the atop site of a carbon atom. With a 5×5 graphene surface, up to 6 hydroxyl groups can be adsorbed on the graphene surface, indicating the concentration coverage of the hydroxyl groups on graphene surface is about 12%. The simulation results show that the negative adsorption energy increases linearly as the number of adsorbed hydroxyl groups increases, and the band gap also increases linearly with the number of adsorbed hydroxyl groups.


I. INTRODUCTION
Graphene has some supreme properties, such as strong mechanical stiffness, strength, good light transmittance, high charge carrier mobility etc. [1][2][3][4][5][6][7][8] and has been used in many fields. The most promising method of large-scale production of graphene material is chemical oxidation of graphite to produce graphene oxide (GO), then followed by reduction [9][10][11]. However, this production chemical method results in defects and oxygen-containing groups such as hydroxyl and epoxide. † Part of Special Issue "John Z.H. Zhang Festschrift for celebrating his 60th birthday". * Author to whom correspondence should be addressed. E-mail: dywang@sdnu.edu.cn Graphene produced using this method is called reduced graphene oxide (rGO).
The distribution of functional groups determines the properties of rGO. The structure of rGO depends to a large extent on the purification procedure [12], the distribution of oxide functional groups largely depends on the degree of oxidation. Nonetheless, till now, there has been no common agreement on the specific electronic structure of rGO. Scientists have proposed a variety of models of the configurations for rGO [13][14][15][16][17][18][19][20][21][22][23]. Among them, the Lerf-Klinowski model is the mostly recognized one [18,22,23]. In this model, the main functional groups adsorbed on the graphene surface are the hydroxyl groups on the atop sites of the C atoms, and the epoxy groups locate on the bridge sites. Later, this model is updated with the carbonyl and carboxyl groups attached at the graphene edges [24].
There are also theoretical studies to predict the rGO structures [25][26][27]. Boukhalov and Katsnelson [25] studied the rGO structure with various coverages of hydroxyl and oxygen groups based on density function theory (DFT). They found that rGO with less than 25% coverage only contains OH groups, as the OH adsorbed at atop sites of the carbon atoms. The 100% coverage of the surface with both OH and O groups is less energetically favorable than 75%. Ghaderi and Peressi [27] performed pseudopotential DFT to study the adsorption of hydroxyl groups on perfect graphene and defective graphene, also with epoxy groups on the surface. Their study shows that defected graphene is much stronger than that at pristine surface. In addition, water formation from adsorbed hydroxyl groups can occur over pristine surface.
In this work, without assuming any adsorption position and processes, based on ab initio molecular dynamics (AIMD), our purpose is (i) to investigate the dynamics process of the OH group adsorption on graphene, (ii) to find the maximums number of hydroxyl groups adsorbed on graphene, (iii) to determine the relationships between the number of adsorbed hydroxyl groups with the adsorption energies and band gaps.

II. METHOD
We used the Vienna ab initio Simulation Package (VASP) [28,29] package to simulate hydroxyl adsorption processes on graphene, with the AIMD method based on density functional theory with the DFT-D3 [30] correction. For the DFT calculation, we used the Projector-Augment-Wave (PAW) [31] pseudopotential, and the Perdew-Burke-Ernzerhof (PBE) [32] functional under generalized gradient approximation was used for the exchange correlation term.
The Brillouin zone of graphene surface was set up as a 5×5 periodic cell. There is a 15Å vacuum gap between the periodic graphene layers. A 5×5×1 k-points with a 400 eV cutoff energy was chosen to sample the Brillouin zone for the AIMD simulation; but in order to get more accurate structures and energies for the stationary points along the adsorption, we used a 520 eV cutoff energy to optimize the structures and energies of the stationary points. The conditions for the electronic degrees of freedom relaxation between the two steps is less than 1×10 −5 eV and the convergence conditions for ionic relaxation of all the forces are less than 0.02 eV/Å. For the AIMD calculation, we do not assume any adsorption paths and adsorption positions. We put the OH hydroxyl group 7.5Å above the graphene surface with a random orientation and speed assignment according to a relative translational energy of 0.01 eV. Then we let the OH hydroxyl group evolve on the graphene surface based on the potential energy from the ab initio calculation. After the first OH group is adsorbed on the graphene surface, we let the second OH group attack the graphene, and so on, until there are no more OH groups which can be adsorbed on the graphene surface. Then we can obtain the maximum number and configuration of hydroxyl groups adsorbed on the graphene surface. For every OH adsorption case, the adsorption energy is calculated with the following formula, Where N represents the number of OH groups adsorbed on the graphene surface, and E G+NOH , E G and E OH are the energies of the entire system with adsorbed hydroxyl groups, the energy of graphene, and the energy of a single hydroxyl radicals in the gas phase. The band gap is also analyzed at a much higher accuracy of 15×15×1 k-points to get the accurate band [33] structure. In total, 196 trajectories were run to investigate the simulation of the OH hydroxyl group adsorption on graphene surface. The propagation time for each trajectory is about 2 ps with a time step of 1 fs. For every adsorption trajectory, the stationary points (physical adsorption complex, transition state and product complex) were searched and identified. The dimer [34][35][36][37] method was employed to search the transition state which is confirmed with one imaginary frequency.

A. Adsorption of the first OH group
Eighteen AIMD trajectories were run for the first hydroxyl group coming from the asymptotic region to attack the graphene and each of trajectories showed that the OH is adsorbed at the atop site of a C atom on the graphene surface.
The AIMD trajectories show that as the reactant OH attacks the graphene surface, it first passes through a potential minimum which is a physical adsorption complex; then it surmounts a potential maximum, con-  FIG. 1(a) is the reactants with the OH is 7.5Å above the graphene with the OH bond distance of 0.98Å. For the physical adsorption complex, the OH groups locate at the bridge site above two C atoms, with one O−C distance of 2.89Å and the other one of 3.02Å. The angle ∠HOC is 156.0 • . As the reaction evolves into the transition state at FIG. 1(c), which is confirmed with an imaginary frequency at 232 cm −1 ; moreover, the O−C distance is decreased to 2.08Å, and the angle ∠HOC is decreased to 99.8 • . At the product complex, it shows that the OH group is adsorbed on the atop site of one C atom, which is consistent with a previous study [27]. The O−C distance is further decreased to 1.51Å, whereas the ∠HOC, 106.5 • , is a little larger than the one of the transition state. Throughout the whole reaction process, the OH distance is always kept at 0.98Å, which means the OH actually serves as a spectator for the whole reaction process.
The schematic potential reaction profile for the first OH adsorption along the reaction path is shown in FIG. 2. The physical adsorption complex is about 0.47 eV lower than the reactants, while the transition state is about 0.03 eV higher than the physical adsorption complex, and the transition state is lower than the reactants with a clean graphene surface, every C atom has the same chemical property that OH can be equivalently adsorbed on any atop site of C atoms.

B. Adsorption of the second OH group
With one OH already adsorbed on the graphene surface, the second OH group is placed at the asymptotic region to collide with the surface. We ran 28 AIMD trajectories for this case, and there are 26 of them with the OH group adsorbed on the graphene, and two of them form into water molecule OH+OH * →H 2 O+O * , and the produced O atom is adsorbed at the bridge site to form epoxy group with the C atom. Since our focus of this study is the OH adsorption on graphene, the H 2 O reaction route will not be discussed in this work.
The stationary points have also been identified and optimized from the second OH adsorption trajectories. The reactants, physical adsorption complex, transition state and product state are displayed in FIG. 3. As the second OH comes in, it first forms the physical adsorption complex with the OH at the bridge above two C atoms. One O−C distance is 2.64Å and the other O−C distance is 2.75Å. Then one of the O−C distance decreases to 2.04Å at the transition state, finally it reaches 1.52Å at product state. After the second OH is adsorbed on the atop site of a C atom on the surface, the two adsorbed OH hydroxyl groups almost have the same geometries with the O−C distance around 1.52Å and ∠HOC is around 106.0 • with OH at 0.98Å. It is worth mentioning that the orientation of the first adsorbed hydroxyl group is also changing during the reaction process, from facing against the second hydroxyl at the reactant state to facing towards the second hydroxyl at the product state.
The transition state for the second OH adsorption is confirmed to have one imaginary frequency of 306 cm −1 . The schematic plot of the second OH adsorption on graphene along the adsorption reaction path is also shown in FIG. 2. We find in this case the transition state barrier is the same as the first OH adsorption case and lower than its corresponding reactants of −0.44 eV, and the physical adsorption complex is about 0.51 eV lower than the reactants. The product state is about 1.11 eV lower than the reactants.

C. Adsorption of the third OH group
There are 32 trajectories run for the third OH group attacking the graphene surface. Among them, there are 21 with the third OH group adsorbed on the surface, the rest with the third OH reacts with the OH * absorbed on surface to form water molecules. In addition, these wa- ter molecules are not adsorbed on the graphene surface in this case.
Again, the stationary points along the reaction trajectories were searched and identified. As the third OH coming from the asymptotic region, it also passes through a physical adsorption complex (FIG. 4(b)) with the OH locating at the bridge site of two C atoms. One of the O−C distance is 2.56Å and the other is 2.70Å, with ∠HOC of 93.7 • . At the transition state (FIG. 4(c)), the distance between the O atom of the hydroxyl group and the C atom is 2.06Å, and the ∠HOC is 99.1 • . At the product state (FIG. 4(d)), the third hydroxyl group is also adsorbed at the atop site, having the identical geometries with the other two adsorbed OH hydroxyl groups except for their spatial orientations. The three adsorbed OH hydroxyl groups all have the O−C bond distances of 1.50Å, the ∠HOC angle of ∼106.4 • .
The transition state is confirmed to have one imaginary frequency of ∼226 cm −1 . The reaction profile along the adsorption reaction path is also plotted in

D. Adsorption of the fourth, fifth and sixth OH group
As subsequent OH group coming in to collide with the OH-adsorbed graphene surface, we found that the 5×5 graphene surface can adsorb up to six OH hydroxyl groups. As the seventh OH comes in from the asymptotic, the OH will not be adsorbed on the substrate and react with the six absorbed hydroxyl groups to form water molecules and epoxy groups. Therefore, the maximum adsorption concentration of OH hydroxyl on G surface is about 12 % with respect to the number of C atoms in the 5×5 super cell.
For the fourth, fifth, and sixth hydroxyl radical adsorption processes, we did not locate the transition states and the physical adsorption structures on their trajectories. FIG. 5(a−c) shows the geometries with four, five, and six OH hydroxyl groups adsorbed on the G surface. For the adsorbed hydroxyl groups of these three cases, each OH group almost has the identical adsorption structure with the O−C bond length around 1.50Å, the O−H bond length about 0.99Å, and the ∠HOC angle around 106.5 • . Table I shows the number of AIMD trajectories we ran, the number of the adsorption trajectories, and the  number of forming water trajectories from the fourth to the seventh OH hydroxyl adsorption. For instance, for the fourth OH adsorption, we ran 30 AIMD simulations, 17 paths show the OH group adsorbed on the substrate; 13 simulations of them form water molecules and epoxy groups with the already adsorbed hydroxyl groups, of which one of them has the water molecule adsorbed on the surface of the substrate. Table I indicates that as the number of hydroxyl groups adsorbed on the surface increases, the incoming OH group is more likely to react with the adsorbed hydroxyl groups on the surface to form water molecules and epoxy groups. And as the number of OH adsorption increases, the generated water molecules are more likely to be adsorbed on the surface. We think that it might be that, as the number of hydroxyl groups on the surface increases, the adsorption interactions between the OH groups and the substrate increase. This is proven by the increase of the average adsorption energy because; for the fourth, fifth and sixth OH group adsorbed on the substrate, the adsorption energies decrease from −1.03 eV for the fourth OH adsorption to −1.14 eV for the fifth OH, then to −1.23 eV for the sixth OH.
Three energy band diagrams are shown in FIG. S5, S6, and S7 (Supplementary materials) for the four, five and six adsorptions. And we find the largest lap for these three cases are 0.57 eV, 0.78 eV and 0.87 eV, respectively.

E. Relationship between the adsorption energy, band gap and the number of adsorbed OH groups
The relationship between the adsorption energy and the number of adsorbed OH groups on the graphene surface is plotted in FIG. 6(a). It has a linear relationship between them as the number of the adsorbed OH groups increases, the adsorption energy becomes more negative which means the interactions between the adsorbed hydroxyl groups and the surface increase. However, as the number of the OH on G surface reaches six, which is the maximum number we find the hydroxyl groups can be adsorbed on surface, the density of the OH is too high, which results in other coming-in OH interaction with the adsorbed OH groups to produce water molecules and epoxy groups on the surface. FIG. 6(b) shows the relationship between the largest band gap and the number of adsorbed OH groups on graphene surface. It shows that the band gap basically increases linearly with the increase of the number of hydroxyl groups. The relationship between the average adsorption energy versus the largest band gap is plotted in FIG. 6(c). There is also a linear relationship between the two, as the negative adsorption energy increases the band gap increases too. Note that the largest band gap with four OH groups is off the linear-fitted lines, in FIG. 6 (b) and (c), more than the other cases. This might be that we only run limited numbers of trajectories, that the largest band gap just does not show up in the 17 trajectories adsorbed in the four OH case.

IV. CONCLUSION
The ab initio molecular dynamics method was used to study the adsorption of hydroxyl groups on graphene surface. We found for the first, second and third hydroxyl adsorption, the OH group needs to pass a physical adsorption complex, then surmount a barrier to be adsorbed on the graphene surface. The barrier height is almost the same at about −0.44 eV with all three cases. In the physical adsorption complex, the hydroxyl group is located above the bridge site of two C atoms. The maximum number of the hydroxyl groups adsorbed on the 5×5 graphene supercell is six, indicating the concentration coverage of the hydroxyl groups on graphene is about 12%. As the number of adsorbed hydroxyl groups increases on the surface, the negative adsorption energy increases linearly, showing that the exoergic interaction increases as the number of adsorption increases. Furthermore, the band gap also increases linearly as the number of adsorbed hydroxyl groups increases on the surface. Our study here provides dynamics information of the OH adsorption dynamics process on graphene and the structure of rGO with the adsorption of hydroxyl group.
Supplementary materials: Figures of the potential energy versus time and the band gaps with corresponding structures are available.

V. ACKNOWLEDGMENTS
This work is supported by the National Natural Science Foundation of China (No.11774206) and Taishan Scholarship Fund from Shandong Province. The computation was carried out at the Shenzhen Supercomputer Center in China.