Diffusion Monte Carlo Study of the O$_2$ Adsorption on a Single Layer Graphene

TL;DR

This study uses diffusion Monte Carlo calculations to accurately analyze O₂ adsorption on single-layer graphene, revealing preferred orientations, adsorption sites, and diffusion barriers, with results aligning well with experimental data.

Contribution

The paper provides the first DMC-based detailed characterization of O₂ adsorption sites, orientations, and energies on graphene, challenging previous DFT predictions.

Findings

O₂ prefers a horizontal orientation at equilibrium distance.

Midpoint of C-C bond (bridge site) is energetically favored for adsorption.

O₂ diffusion barrier on graphene is approximately 11 meV.

Abstract

Diffusion Monte Carlo (DMC) calculations were performed for an accurate description of the nature of the O$_2$ adsorption on a single layer graphene. We investigated the stable orientation of O$_2$ at a specific adsorption site as well as its equilibrium adsorption energy. At equilibrium adsorption distances, an O$_2$ molecule was found to prefer a horizontal orientation, where the O-O bond is parallel to the graphene surface, to the vertical orientation. However, the vertical orientation is favored at the O$_2$-graphene distances shorter than the equilibrium distance, which could be understood by the steric repulsion between O and C atoms. Contrary to previous DFT calculations, our DMC calculations show that the midpoint of a C-C bond (a bridge site) is energetically preferred for the O$_2$ adsorption to a center of a hexagonal ring (a hollow site). The lowest DMC adsorption energy was found at an intermediate point between a hollow and a bridge site, where the O$_2$ adsorption energy was estimated to be -0.142(4) eV that was in very good agreement with the recently-reported experimental value. Finally, we have found that O$_2$ is very diffusive on the surface of graphene with the diffusion barrier along a bridge-hollow-bridge path being as small as ~ 11 meV.

Figures

Tables

Table 1. Table 1: Extrapolated DMC adsorption energies E a ​ d subscript 𝐸 𝑎 𝑑 E_{ad} (eV) of O 2 at the equilibrium distances d e ​ q subscript 𝑑 𝑒 𝑞 d_{eq} (Å) from a hollow and a bridge site, along with the corresponding DFT results based on several different exchange-correlation functionals.

			PBE	DFT-D2	vdW-DF	vdW-DF2	rVV10	DMC
hollow	$A$ mode	$d_{eq}$	3.63	3.17	3.28	3.14	3.10	3.28(5)
		$E_{ad}$	-0.016	-0.148	-0.162	-0.145	-0.141	-0.130(4)
	$B$ mode	$d_{eq}$	3.63	3.17	3.28	3.14	3.10	3.30(5)
		$E_{ad}$	-0.016	-0.148	-0.164	-0.146	-0.142	-0.126(4)
	$V$ mode	$d_{eq}$	3.46	3.05	3.12	2.97	2.92	3.26(9)
		$E_{ad}$	-0.014	-0.121	-0.130	-0.116	-0.110	-0.111(7)
bridge	$A$ mode	$d_{eq}$	3.66	3.22	3.32	3.19	3.15	3.18(3)
		$E_{ad}$	-0.017	-0.141	-0.158	-0.139	-0.132	-0.117(7)
	$B$ mode	$d_{eq}$	3.71	3.23	3.33	3.21	3.16	3.21(2)
		$E_{ad}$	-0.016	-0.135	-0.156	-0.133	-0.130	-0.135(4)

Equations

$$E(\mathrm{graphene},{\mathrm{O}}_2)-E(\mathrm{graphene})-E({\mathrm{O}}_2)$$

$${\rho }_{inter}(x,y)=\rho (x,y,z_{gr}+r_{\mathrm{C}}^c<z<z_{{\mathrm{O}}_2}-r_{\mathrm{O}}^c),$$

Full text

quantum Monte Carlo method, density functional theory, oxygen, graphene, van der Waals interaction

I Introduction

Graphene has drawn much interest from various areas of chemistry, physics and industry because of its unique electronic and optical properties.Novoselov et al. (2005); Zhang et al. (2005); Geim and Novoselov (2007); Geim (2009) Among others, graphene-based materials have been identified as promising candidates for future enhanced gas sensor, beyond carbon nanotubesKong et al. (2000) or nanowires,Cui et al. (2001) because of their high reactivity to adsorbates. An adsorbed gas molecule on a graphene surface was shown to play a fundamental role as a donor/acceptor and to alter significantly its transport properties.Schedin et al. (2007); Lu et al. (2009); Varghese et al. (2015) Additionally, its atomically-thin two-dimensional hexagonal structure allows all carbon atoms to work as surface atoms for molecular adsorbates, yielding a large sensing area per volume. A number of theoretical and experimental studies have recently attempted to quantify and describe the performance of a graphene-based material, including a pristine graphene sheet, as a gas detector.Wehling et al. (2008); Lee and Kim (2013); Chen et al. (2014); Ganji et al. (2015) For example, an exfoliated graphene on SiO2 was reported to be capable of detecting NO2, a well-known critical air pollutant, at very low concentrations of a few parts per billion.Schedin et al. (2007) Other graphene-based materials, such as carbon nanotubes, graphene oxides and doped graphene, have also been considered as sensing substrates to detect not only NO or NO2 but other gases such as NH3, CO2 and O2.Ratinac et al. (2010); Yuan and Shi (2013); Jaaniso et al. (2014); Novikov et al. (2016) The adsorption energy of a molecule on a substrate is a key metric in identifying a gas-sensing material. In the past decades, first-principle density functional theory (DFT) calculations have been widely used to investigate energetics of different gas adsorptions on the surface of a graphene-based material.Huang et al. (2008); Leenaerts et al. (2009); Zhou et al. (2011); Choudhuri et al. (2015); Liang et al. (2017) However, the DFT results were often strongly dependent on the choice of exchange-correlation functional, which can be attributed to some limitation of its Kohn-Sham mean-field scheme in describing subtle competition between the covalent nature of the adsorbate-substrate interaction and its van der Waals (vdW) nature.

Oxygen dimer (O2) is the second most abundant molecule in the atmosphere. Among other usages, O2 can be used to control the rate of combustion or to destruct hazardous and waste materials in incinerators. It needs to be removed from the atmosphere to avoid corrosion. Because of the crucial importance in its industrial application, the adsorption nature of O2 has been intensively studied for various sensing materials including carbon allotropes.Kauffman et al. (2009); Rajavel et al. (2015); Tang et al. (2017) Recently, the adsorption energy of O2 on graphene was experimentally measured using temperature-programmed terahertz emission microscopy.Bagsican et al. (2017) However, detailed information about the adsorption site and equilibrium distance along with its preferred molecular orientation is yet to be accurately determined.

On the theoretical side, several DFT calculations were performed to study the adsorption of O2 on graphene. While a geometry relaxation of Yan et al. Yan et al. (2012) based on a GGA functional showed the adsorption of O2 at a hollow site, the center of a hexagonal ring on the graphene surface, with the adsorption energy of -0.13 eV, Guang et al. used a GGA functional with the addition of empirical vdW dispersion of Grimme to find an O2 adsorption at a bridge site, the midpoint of a C-C bond, with the adsorption energy of -0.15 eV Guang et al. (2013). On the other hand, two independent LDA calculations predicted the adsorption of O2 at a hollow site Zou et al. (2011); Kaur et al. (2018). This discrepancy in previous DFT calculations asks for a more systematic theoretical study for the O2 adsorption on graphene with a computational method that fully incorporates electron-electron correlations.

Quantum Monte Carlo (QMC) can provide accurate estimations of the adsorption energy by solving a many-body Schrödinger equation stochastically. QMC was used successfully to study various vdW-dominated systems and provided accurate description of vdW interactions in molecular systems, bulk solidsBenali et al. (2014), and 2D materials including bilayer graphene and other low-dimensional carbon allotropes.Shin et al. (2014); Mostaani et al. (2015); Shin et al. (2017a) In this study, we report diffusion Monte Carlo (DMC) results for the energetics, the adsorption distances, and preferred orientations of an O2 molecule adsorbed at different sites of a single layer graphene. We then estimate the O2 diffusion barrier on the graphene surface using both DMC and vdW-corrected DFT calculations.

II Computational details

Our QMC calculations were carried out with the fixed-node DMC methodReynolds et al. (1982); Foulkes et al. (2001) as implemented in the QMCPACK code.Kim et al. (2018) Slater-Jastrow trial wave functions were used in the DMC algorithm with up to three-body Jastrow functions in order to include ion-electron (one-body), electron-electron (two-body), and ion-electron-electron (three-body) correlations. All DMC calculations used a time step of 0.005 Ha*-1* and size-consistent T-moves for the variational evaluation of the non-local pseudopotentials.Casula et al. (2010) In order to reduce finite-size effects, DMC total energies of four different supercells averaged over 36 twist angles (twist-averaged boundary conditionsLin et al. (2001)) were extrapolated to the bulk limit. For atomic calculations of both carbon and oxygen, we used norm-conserving scalar-relativistic energy-consistent pseudopotentials developed by Burkatzki, Filippi, and DolgBurkatzki et al. (2007, 2008), whose accuracies were confirmed in various DMC studies on similar materials.Shin et al. (2014); Benali et al. (2016); Luo et al. (2016); Shin et al. (2017a, b). The single Slater determinant in the wave function was constructed with single-particle orbitals obtained through spin-polarized DFT calculations based on the Perdew-Burke-Ernzerhof (PBE) parameterizationPerdew et al. (1996), which were done with the QUANTUM ESPRESSO codeGiannozzi et al. (2009). To analyze effects of the dispersion force on the O2 adsorption, DFT calculations were also performed using several vdW-corrected functionals; Grimme (DFT-D2) correction to PBE exchange correlation (XC) functionalGrimme (2004, 2006); Grimme et al. (2010), self-consistent vdW-corrected XC functionals and more recent rVV10 correction.Dion et al. (2004); Lee et al. (2010); Sabatini et al. (2013) In order to minimize the interactions between O2 and its periodic images, we used a $3\times 3\times 1$ single layer graphene with one adsorbed O2 molecule in its triplet state. The distance between O2 and its first periodic image is around 7 Å, which corresponds to the dissociation limit of O2-O2 dimer.Aquilanti et al. (1999); Hernández-Lamoneda et al. (2005) Self-consistent spin-polarized DFT calculations were done using 250 Ry plane-wave basis-set cut-off and $6\times 6\times 1$ Monkhorst-Pack $k$-point grid mesh.Monkhorst and Pack (1976) The O-O bond length in O2 was chosen to be 1.208 Å from the experimental result for the equilibrium O-O bond lengthLide (2003) and the experimental value of 1.421 Å for the C-C bond length of graphite was used for that of a single layer graphene.Lynch and Drickamer (1966)

III Results

We first investigate the adsorption of a single O2 molecule at a hollow site, which was reported in previous DFT studies Zou et al. (2011); Yan et al. (2012); Kaur et al. (2018) to be a stable adsorption site for O2. To determine the preferred orientation of O2 adsorbed at a hollow site, we here consider its three different orientation modes as shown in Figure 1; the O2 molecular axis is oriented in the directions parallel to the graphene surface in the first two modes called $A$ and $B$ modes but it is perpendicular to the graphene surface in the $V$ mode. Ideally, the most stable mode could be determined through a geometry optimization of the O2-graphene complex. A robust DMC algorithm to compute forces, however, has yet to be available. Therefore, we reduce the degrees of freedom of the system only to a vertical translation of O2 with respect to the graphene surface and then compute DMC energies as a function of the distance from the surface. This allows us to construct a binding curve for each orientation mode. Because a localized puckering of the graphene surface is expected to occur near the adsorption site of O2, we performed a full DFT optimization of the adsorbant and the graphene surface using rVV10 vdW-corrected exchange-correlation functional to assess the level of deformation and its associated energy. It was found that geometries and total energies varied only a little, less than 0.02 Å and 10 meV, respectively, with the structural deformation. From this we conclude that our approach of fixing the geometry of the graphene sheet and the O-O bond length of O2 introduces a rather small bias that should not affect the quality of the calculations or the conclusions of this study.

Following the scheme described above, we performed DMC calculations to compute the adsorption energy of a single O2 molecule, defined by

[TABLE]

where $E({\rm graphene},{\rm O}_{2})$, $E({\rm graphene})$, and $E({\rm O}_{2})$ are, respectively, the total energies of graphene with the O2 adsorbate, the pristine graphene, and an isolated O2 molecule. Figure 2 shows DMC adsorption curves of O2 at a hollow site that were simulated with a $2\times 2\times 1$ supercell involving 4 O2 molecules and 72 C atoms. As can be seen, the equilibrium binding distances are $3.2$ to $3.4$ Å and are quite similar among the three orientation modes. While equilibrium binding energies of $A$ and $B$ modes are indistinguishable from each other within error bars, the $V$ mode turns out to be the least probable one among the three orientations. By choosing the mode with the lowest adsorption energy at a specific distance from the binding energy curves of Fig. 2(b), we determine the preferred orientation profiles of an O2 molecule at a hollow site; O2 is found to prefer an orientation parallel to the graphene surface at the distances longer than 3.0 Å from the graphene surface while the vertical $V$ mode is favored at shorter distances. It is not surprising to find that the $V$ mode is most stable at short distances because O2 in this mode will be less affected by the steric repulsion between C and O atoms when approaching the graphene surface. This mode is expected to be the preferred orientation when O2 crosses through the graphene layer.

The above results suggest a rotation of O2 adsorbed at a hollow site from its parallel ($A$ or $B$) mode to the vertical $V$ mode when getting closer to the graphene surface. Figure 3 presents the energy difference between the $A$ mode and the $V$ mode as a function of the O2-graphene distance $d$. The DMC results show that a transition from the $A$ mode to the $V$ mode takes place around 2.8 Å, which seems to be reproduced, at least qualitatively, by DFT calculations based on various vdW-corrected functionals. In the case of a DFT-PBE calculation without any vdW correction, however, such a transition is hardly noticeable, suggesting that it is mostly driven by the vdW interaction.

For another possible adsorption site for O2, we now consider a bridge site that was predicted in the DFT geometry optimization study of Guang Guang et al. (2013) to be a stable adsorption site. For this only two parallel modes of $A$ and $B$ are investigated because the vertical $V$ mode was found to be less stable at the equilibrium distance than the parallel modes. Figure 4 presents the DMC adsorption energies of O2 at a bridge site, which were computed for the $2\times 2\times 1$ supercell, as functions of the distance from the graphene surface, along with the corresponding results at a hollow site. Unlike an O2 molecule adsorbed at a hollow site that showed no clear preference between $A$ and $B$ modes, O2 at a bridge site is observed to favor the $B$ mode which shows a significantly lower equilibrium adsorption energy than the $A$ mode. This suggests that an O2 molecule at a bridge site is not allowed to rotate freely in the plane parallel to the graphene surface, which can be understood by stronger repulsive interaction due to shorter O-C interatomic distance in the $A$ mode than in the $B$ mode. Furthermore, our DMC calculations for the $2\times 2\times 1$ supercell reveal that the $B$-mode equilibrium adsorption energy at a bridge site is slightly lower than, or nearly identical to, the equilibrium adsorption energies at a hollow site.

To understand the reason why the $B$ mode is preferred at a bridge site over the $A$ mode, we computed two-dimensional charge density distribution in the intermediate region between O2 and the graphene surface:

[TABLE]

where $z_{gr}$ and ${z_{\rm O_{2}}}$ denote $z$ coordinates of graphene and the oxygen dimer, respectively, and $r^{c}_{\rm C}$ ($r^{c}_{\rm O}$) represents the covalent radius of a carbon (oxygen) atom. Figure 5 shows contour plots of $\rho_{inter}$ that were computed using DMC (left), PBE-DFT (middle) and rVV10-DFT (right) methods. One can see clear difference in the DMC density peaks between $A$-bridge and $B$-bridge modes. Smaller amount of charges distributed in the intermediate region for the $A$ mode than for the $B$ mode suggests that electrons in the $A$-mode configuration tend to be pushed away from the intermediate region to reduce the repulsion due to the overlap between electron clouds of oxygen atoms and those of carbon atoms. Note that oxygen atoms are closer to their nearest carbon atoms in the $A$-bridge mode than in the $B$-bridge mode. This difference in the charge distributions between $A$ and $B$ modes is understood to account for the lower DMC adsorption energy of the $B$ mode at a bridge site. On the other hand, both PBE and rVV10 calculations do not show much difference in the amount of the intermediate charges between the $A$-bridge and the $B$-bridge modes, reflecting that the adsorption energies of these two orientation modes are similar to each other in both DFT calculations. From this we conclude that the inability of DFT to distinguish the adsorption energy between the two parallel orientation modes comes from its limitation in describing the electron-electron correlation, especially among the electrons in the intermediate region. It can be also seen in Fig. 5 that the DMC charge density distribution of the $B$-bridge mode do not show much difference from the corresponding result for the $B$-hollow mode, accounting for similar DMC adsorption energies between these two adsorption modes.

To understand different nature of the O2 adsorption at a bridge site from the hollow-site adsorption, we now investigate the extent of charge redistribution induced by the graphene-O2 interaction. Figure 6 presents the difference between one-dimensional charge density distribution of the O2-graphene complex and that of isolated graphene and O2, which is projected to the vertical axis. It can be seen that the amounts of charges redistributed through the O2 binding processes at both adsorption sites are larger in DMC calculations than in PBE and rVV10 calculations. For both hollow- and bridge-site adsorptions, no noticeable charge accumulation is observed at the intermediate region between graphene and O2, confirming the van der Waals nature of the O2 binding to graphene. According to our DMC results, charges of graphene tend to be pushed to the opposite side of O2 in both hollow- and bridge-site adsorptions. A significant qualitative difference between DMC and DFT results is observed in the bridge-site adsorption, where the DMC result shows charges in O2 being pushed away from the graphene side while both PBE and rVV10 calculations produce minimal charge redistribution on the O2 side. We understand that lack of charge redistribution around the bridge-site O2 at the DFT level accounts for no significant difference in the DFT adsorption energies between the $A$-bridge and the $B$-bridge modes.

Finally we make finite size analysis of our DMC calculations to determine the most favorable adsorption site for O2 along with its preferred orientation. To reach the thermodynamic limit for each adsorption mode, we extrapolate DMC adsorption energies computed at equilibrium distances for three different supercell sizes to the bulk limit, more details of which can be found in the Supplemental Material. 111See Supplemental Material at [URL will be inserted by publisher] for more information on finite-size analysis for the adsorption energy of O2 on graphene. Our extrapolated DMC results for equilibrium adsorption energies ($E_{ad}$) and the equilibrium distances ($d_{eq}$) are summarized in Table 1, which also presents the corresponding DFT results computed with several different exchange-correlation functionals. While DFT calculations with vdW-corrected functionals show clear binding of O2 at both adsorption sites, the PBE calculation without a dispersion force predicts that an O2 molecule is barely bound to either a hollow or a bridge site with its estimated binding energy less than 0.02 eV. This reflects the van der Waals nature of the O2 binding to the graphene surface. We note that vdW-DFT calculations tend to overestimate the O2 binding energies but to underestimate its equilibrium distances when compared to the DMC results. Among exchange-correlation functionals considered here, the O2 adsorption energies are found to be best described, except for the $A$-bridge mode, with the rVV10 functional while vdW-DF (vdW-DF2) functional gives rise to the O2-graphene equilibrium distances at a hollow (bridge) site that are most consistent with our DMC results. No density functional is found to be able to recover simultaneously the DMC adsorption energies and equilibrium distances of O2 on graphene. This limitation of the DFT calculations was also observed in recent theoretical studies for the interlayer binding of phosphorus Shulenburger et al. (2015); Ahn et al. (2018) and carbon layers. Shin et al. (2017a)

According to the results presented in Table 1, regardless of the computational method used, two parallel modes of $A$ and $B$ have similar adsorption energies and equilibrium distances at a hollow site. This suggests a free planar rotation of O2 at a hollow site, which can be triggered by temperature or any external perturbation. While all DFT calculations predict that these two parallel modes are nearly degenerate in energy even at a bridge site, our DMC results show that the $B$ mode is clearly favored over the $A$ mode for O2 adsorbed at a bridge site. This indicates that free planar rotation between $A$ and $B$ modes is no longer allowed for O2 adsorbed at a bridge site, which is understood by stronger influence of repulsive interaction due to shorter O-C interatomic distance than for O2 at a hollow site. Furthermore, all DFT calculations presented in Table 1 favor the adsorption of O2 at a hollow site over its adsorption at a bridge site, regardless of the orientation of O2. This is in contrast to our DMC result that the adsorption energy in the $B$-bridge mode is slightly lower than its hollow-site adsorption energies.

The DMC adsorption energies presented in Table 1 are a little higher than the experimentally-reported value of -0.15 eV for the O2 adsorption energy on graphene. Bagsican et al. (2017) This leads us to explore the graphene surface for more favorable adsorption sites for O2 than either a bridge or a hollow site and to investigate in the process diffusion mobility of O2 on the graphene surface. For this, the adsorption energy of O2 is computed as it moves ”point by point” in the direction parallel to the graphene surface from a hollow site to a bridge site (see the top panel of Fig. 7) while the O2-graphene distance is fixed to its hollow-site equilibrium distance. The bottom panel of Fig. 7 presents the O2 adsorption energies along the diffusion path relative to the corresponding values at a hollow site. This reflects the above finding that the $B$ mode has a much lower DMC adsorption energy than the $A$ mode at a bridge site (note that the $A$ mode is nearly degenerate to the $B$ mode at a hollow site). And it clearly shows the $B$-mode DMC adsorption energy is lower at a bridge site than at a hollow site. Furthermore, our DMC results show that the $A$ orientation mode of O2 has the lowest adsorption energy at an intermediate point between a hollow and a bridge site, where the adsorption energy is lower than, at least as low as, the adsorption energy of the $B$-bridge mode. The DMC $A$-mode adsorption energy at an intermediate site is estimated to be $-0.142(4)$ eV, which is quantitatively consistent with the recently-reported experimental value. Bagsican et al. (2017) According to our DMC results, either an intermediate point ($A$ mode) or a bridge site ($B$ mode), depending on the orientation mode, is energetically the most stable adsorption site for O2. On the other hand, all DFT results presented in Fig. 7 show that the adsorption energy of O2 increases steadily as it moves from a hollow site to a bridge site, indicating that a bridge site is the peak point in the diffusion path for O2 on the graphene surface. This contradicts our DMC prediction that a hollow site is the peak point in the O2 diffusion path. As far as the diffusion barrier of O2 along a bridge-hollow-bridge diffusion path is concerned, our DMC estimate is 11(4) meV, similar to its rVV10 and vdW-DF2 estimates.

IV Conclusions

We have investigated the relative stability of three orientation modes of O2 adsorbed on a single layer graphene using DMC and vdW-corrected DFT calculations. While the two orientation modes parallel to the graphene surface are found to be more favorable than the vertical mode at the equilibrium adsorption distances, the vertical mode is favored at shorter distances where strong repulsive force between graphene and O2 is more dominant. While vdW-corrected functionals were able to describe, at least qualitatively, the orientation profile of O2 adsorbed on graphene, no functional was able to reproduce simultaneously the equilibrium adsorption distance and adsorption energy estimated in the DMC calculations.

Similar adsorption energies for its two parallel orientation modes of $A$ and $B$ led us to predict a nearly free planar rotation of O2 adsorbed at a hollow site. On the other hand, our DMC calculations for the O2 adsorption at a bridge site show that the $B$ mode has a significantly lower adsorption energy than the $A$ mode, breaking the rotational symmetry between these two modes. Further DMC investigation of the O2 adsorption reveals that in contrast to DFT calculations, a hollow site is not the most stable adsorption site for O2. The site with the lowest O2 adsorption energy was found to be a bridge site for the $B$ orientation mode and an intermediate site between a hollow and a bridge site for the $A$ mode. The stable $A$-mode adsorption of O2 at an intermediate site can be attributed to the interplay between repulsive interaction and vdW interaction. As an O2 molecule in the $A$-mode orientation approaches a bridge site from a hollow site, an attractive interaction due to dispersion force gets stronger but the strength of the repulsive interaction due to charge overlap between O and C atoms also increases, resulting in the lowest adsorption energy at an intermediate site.

Unlike DFT results based on several different exchange-correlation functionals, our DMC study predicts that the $B$-bridge mode is one of the stable adsorption modes. Based on low diffusion barrier of $\sim 11$ meV, we conclude that the adsorbed O2 molecule is very diffusive on the graphene surface. The most stable O2 adsorption site is concluded to be degenerate between an intermediate site and a bridge site, depending on the molecular orientation of O2. This can be hardly reproduced within DFT framework because the O2-graphene interaction in the $B$-bridge mode is not accurately described even with vdW-corrected exchange-correlation functionals. Finally we expect that our DMC density distributions can provide a guideline to find better exchange-correlation functionals for accurate theoretical estimations of gas sensing properties of graphene, possibly other 2D materials.

Acknowledgements.

H.S, Y.L and A.B were supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, as part of the Computational Materials Sciences Program and Center for Predictive Simulation of Functional Materials. Y.K was supported by the Basic Science Research Program (2018R1D1A1B07042443) through the National Research Foundation of Korea funded by the Ministry of Education. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC02-06CH11357.