Quantum Monte Carlo investigations of adsorption energetics on graphene

For an N-electron system, the diffusion quantum Monte Carlo (DMC) method is, in principle, an exact approach for solving the imaginary-time Schrödinger equation

$$\begin{eqnarray} -\frac{\partial \psi (\mathbf{R},\tau )}{\partial \tau }&=&(\hat {H}-{E}_{T})\psi (\mathbf{R},\tau )\nonumber\\ &=&-\frac{1}{2}{\nabla }_{\mathbf{R}}^{2}\psi (\mathbf{R},\tau )+(V(\mathbf{R})-{E}_{T})\psi (\mathbf{R},\tau ),\nonumber\\ \end{eqnarray} \tag{ 1 }$$

where τ is the imaginary time, $\mathbf{R}=({\vec{r}}_{1},{\vec{r}}_{2},\ldots ,{\vec{r}}_{N})$, constant E_T is an arbitrary energy offset, and V(R) is the Coulomb interaction potential. Equation (1) is written in atomic units, where e = m_e = ħ = 4π₀ = 1. The ground state of the Hamiltonian $\hat {H}$ in equation (1) can be obtained by taking the limit τ → ∞ [19]. Equation (1) can also be regarded as a 3N-dimensional diffusion equation with a drift term, and ψ(R,τ) can be considered as the density of walkers. However, two issues have to be addressed in practical calculations. First, the drift term (V(R) − E_T), corresponding to the Coulomb interaction, may diverge, which will lead to large fluctuations in the population process and will slow down or hinder the simulation. Second, the wavefunction ψ(R,τ) could be either positive or negative, while the density of walkers has to be positive everywhere. This is the well-known fermion sign problem.

The introduction of importance sampling can help smooth the fluctuations and, at the same time, solve the fermion sign problem [20]. Specifically, a guiding wavefunction ψ_G(R) is introduced to guide the diffusion walkers and rewrite equation (1) in terms of the new distribution function f(R,τ) = ψ_G(R)ψ(R,τ):

$$\begin{eqnarray} \fl -\frac{\partial f(\mathbf{R},\tau )}{\partial \tau }&=&-\frac{1}{2}{\nabla }_{\mathbf{R}}^{2}f(\mathbf{R},\tau )+\nabla \boldsymbol{\cdot }[F(\mathbf{R})f(\mathbf{R},\tau )]\nonumber\\ &&+~({E}_{\mathrm{ L}}(\mathbf{R})-{E}_{T})f(\mathbf{R},\tau ), \end{eqnarray} \tag{ 2 }$$

where E_L(R) is the local energy calculated with respect to ψ_G(R), ${E}_{\mathrm{L}}(\mathbf{R})=\hat {H}{\psi }_{\mathrm{G}}(\mathbf{R})/{\psi }_{\mathrm{G}}(\mathbf{R})$, and the vector F(R) = ∇ψ_G(R)/ψ_G(R) is the drift term controlling the velocity of the walkers. If one can choose a reasonably accurate ψ_G(R), E_L(R) will be a well-behaved function without large fluctuations.

The guiding wavefunction ψ_G(R) also plays an important role in the so-called fixed-node approximation [21, 22]. The new distribution function f(R,τ) in equation (2) can be treated as a Boson ground-state wavefunction (without the sign problem) if ψ_G(R) has the nodes of the exact fermion ground state that we are looking for. The fixed-node approximation is designed to limit the walkers to moving within one nodal region, and moves across the nodal surface will be rejected [23]. The introduction of the guiding wavefunction ψ_G(R) sets up the nodal surface in the simulation. For complicated materials systems, the guiding wavefunction ψ_G thus determines the accuracy of DMC results.

The guiding wavefunction ψ_G(R) can be obtained from optimized variational quantum Monte Carlo (VMC) calculations. The VMC method provides an upper bound for the exact ground-state energy. Given the Hamiltonian $\hat {H}$, the variational energy with respect to a trial wavefunction ψ_G is given by

$$\begin{eqnarray} {E}_{G}&=&\frac{\langle {\psi }_{\mathrm{G}}\vert \hat {H}\vert {\psi }_{\mathrm{G}}\rangle }{\langle {\psi }_{\mathrm{G}}\vert {\psi }_{\mathrm{G}}\rangle }\nonumber\\ &=&\int \hspace{0.167em} \frac{\vert {\psi }_{\mathrm{G}}(\mathbf{R})\vert ^{2}}{\int \vert {\psi }_{\mathrm{G}}({\mathbf{R}}^{\prime})\vert ^{2}\hspace{0.167em} \mathrm{d}{\mathbf{R}}^{\prime}}\hspace{0.167em} [\hat {H}{\psi }_{\mathrm{G}}(\mathbf{R})/{\psi }_{\mathrm{G}}(\mathbf{R})]\hspace{0.167em} \mathrm{d}\mathbf{R}. \end{eqnarray} \tag{ 3 }$$

The trial wavefunction has the usual Slater–Jastrow form

$$\begin{eqnarray} &&{\psi }_{\mathrm{G}}(\mathbf{R})={D}^{\uparrow }{D}^{\downarrow }\exp (J), \end{eqnarray} \tag{ 4 }$$

where D^↑ and D^↓ are the up- and down-spin Slater determinants of single-particle orbitals, respectively. These single-particle wavefunctions are usually generated from density functional theory (DFT) or Hartree–Fock calculations. J is the Jastrow factor that contains the electron–electron (u), nucleus–electron (χ), and nucleus–electron–electron (f) terms [24]:

$$\begin{eqnarray} J&=&\mathop{\sum }\limits _{i}\mathop{\sum }\limits _{j=i+1}u({r}_{i j})+\mathop{\sum }\limits _{i}\mathop{\sum }\limits _{I}{\chi }_{I}({r}_{i I})\nonumber\\ &&+~\mathop{\sum }\limits _{I}\mathop{\sum }\limits _{i}\mathop{\sum }\limits _{j=i+1}{f}_{I}({r}_{i I},{r}_{j I},{r}_{i,j}), \end{eqnarray} \tag{ 5 }$$

where r_ij is the distance between electrons i and j, and r_iI denotes the distance between electron i and ion I. The u, χ and f terms in equation (5) are power expansions containing variational parameters and dependent on the spins of electrons. Since the electron–electron Coulomb potential exhibits discontinuous first derivatives with respect to r_ij when two electrons get close, the parameters in u are constrained in such a way that the Kato cusp condition [25] is satisfied. The Slater–Jastrow trial wavefunction in equation (4) can be optimized by minimizing either the variance of the energy [26] or the energy alone [27].

For DMC calculations of an extended system, one needs to consider a finite number of electrons in one large supercell and construct the many-body wavefunction from one-particle orbitals associated with a single k-point in the new Brillouin zone (BZ) of this supercell [28]. To make these calculations feasible, we have used a 4 × 4 hexagonal unit cell including a total of 32 C atoms for both graphene and graphane studies. The one-particle orbitals are obtained from spin-polarized DFT calculations with norm-conserving pseudopotentials. For these systems, the quality of these one-particle orbitals turns out to be weakly dependent on the exchange–correlation functionals employed, therefore, most of our DMC calculations have started with LDA wavefunctions. We have also tested the effect of switching to the PBE wavefunctions in the graphane calculations, but found little difference.

For the single-atom adsorption on graphene, three possible high-symmetry adsorption sites are considered: the bridge (B) site between two C atoms, the top (T) sites on top of a single C atom, and the hollow (H) site at the center of the hexagon. The geometry is shown in figure 1, together with the 4 × 4 unit cell used in the calculations. The structure of graphane is assumed to have fully loaded hydrogen on both sides of the surface.

Zoom In Zoom Out Reset image size

We carry out spin-polarized DFT calculations using norm-conserving pseudopotentials with the plane-wave basis as implemented in the CASTEP package [29]. We use the OPIUM code [30] to generate norm-conserving pseudopotentials and Casula's approach [31] is used in DMC calculations to minimize the nonlocal pseudopotential error. The plane-wave cutoff energy is 2000 eV, and the optimized force on each atom is less than 0.03 eV  Å⁻¹. In the single-atom adsorption study, a 4 × 4 supercell with a lattice constant of 4 × 2.464 Å is used, and the Brillouin zone is sampled by a 3 × 3 k-point grid to obtain the optimized adsorption structures. We then use a single special k point at this relaxed geometry to do the DMC calculations. For graphane, we have used a 1 × 1 unit cell and a 4 × 4 k-point grid in DFT calculations, and its accuracy is checked by denser grids. We then use this initial wavefunction (equivalent to one k point in a 4 × 4 unit cell) in DMC calculations. A layer separation of 10 Å is sufficient to effectively eliminate the interaction between neighboring layers. Calculations with the projector augmented wave (PAW) method as implemented in VASP [32, 33] are also carried out, and the results are very close to those using norm-conserving pseudopotentials. For isolated single atoms and diatomic molecules, we have repeated the calculations using the Quantum ESPRESSO package [34] and obtained consistent results. Different exchange–correlation functionals are compared, including the LDA [15], PBE [16], PW91 [17], and revPBE [18] functionals.

In order to increase the efficiency of the DMC calculations, the one-particle wavefunctions are subsequently fitted by localized B-spline functions [35]. The scale of the B-spline grid spacing is given by the parameter b = π/G_max, where G_max is the length of the largest plane-wave component. A grid spacing of b/1.5 is used in the calculations for adatoms on graphene, while a grid spacing of b/2 is used for the study of graphane. The number of walkers in each DMC calculation is 1000 with a time step of 0.01 au for both graphene and graphane calculations. The trial wavefunctions with a Jastrow form in equation (4) are optimized by the energy minimization method [27]. All the DMC calculations described above are performed using the CASINO code [36].

We first present energetics results for the adsorption of single O, F and H atoms on graphene. The adsorption energy of these atoms on graphene is evaluated with respect to that of the isolated atom: E_ad = E(graphene + atom) − E(graphene) − E(atom). We will also discuss the results with respect to the corresponding molecules later. The calculated energies for the isolated single atoms and molecules are discussed in the appendix.

An important issue in DMC calculations is the choice of the special k-point, which effectively determines the imposed boundary condition for the many-body wavefunctions [28]. A careful selection of this special k point is desirable in order to reduce the finite-size error and to obtain a reasonably accurate energy using a finite supercell. The adequacy of a particular k-point choice can be assessed by DFT calculations that are much less computationally demanding. For the 4 × 4 supercell used in the DMC calculations in this work, we find the M point (1/2, 0) at the zone boundary of the supercell BZ to be a good choice for the special k-point. We use the DFT-LDA calculations to compare the energy results obtained from this special k-point with those using denser k-grids in the same BZ, including a 2 × 2 grid without offset, a 2 × 2 grid with an offset of (1/4, 1/4), and a 3 × 3 grid without offset. These results for a single H atom at three different adsorption sites calculated using the CASTEP package with norm-conserving pseudopotentials are shown in figure 2. It can be seen that the adsorption energies for the top and hollow sites obtained using the special k-point are quite close to those obtained using a 3 × 3 k-point mesh. As shown in figure 2, the LDA energy difference between the results using the special k-point and the 3 × 3 k-point set is less than 0.05 eV for H on the top site, the lowest-energy site. For the high-energy bridge site, the relatively larger energy difference is probably due to the reduction of the graphene C_6v symmetry to C_2v for the bridge site. This choice of the zone-boundary special k-point yields real-valued determinants, gives an excellent energy result for the lowest-energy H adsorption site, and results in a calculated DMC cohesive energy of 7.385 eV/atom for graphene, which is very close to the experimental value of the cohesive energy of graphite (7.374 eV/atom) [37].

Zoom In Zoom Out Reset image size

For the adsorption of a single O atom at the bridge site, the LDA adsorption energy difference between the special k-point and the 3 × 3 k-point set is also small, only about 0.16 eV. For this case we have further carried out more DMC calculations for each of the point on a 3 × 3 k-grid for the supercell and obtained a final average that is 0.29 eV higher than the energy value using the single special k-point. On the other hand, the LDA overbinding amount for this geometry is around 1.7–2.0 eV, much larger than the energy uncertainty resulting from the k-point choice.

To check the quality of the calculations with norm-conserving pseudopotentials, we have also performed a calculation using the VASP package with the PAW method. As shown in figure 2, the adsorption energy results for a single H in a 4 × 4 supercell with the 3 × 3 k-mesh agree with those using the CASTEP package with norm-conserving pseudopotentials to within 4–17 meV for the three H adsorption sites. The calculations are spin-polarized, resulting in an unpaired electron centered at the H site and a net magnetic moment of 1 μ_B per supercell.

The DMC adsorption energies for a single O, F, and H atom at different adsorption sites on graphene are plotted in figure 3 and listed in table 1, in comparison with DFT results using different exchange–correlation functionals, including the LDA [15], PBE [16], PW91 [17], and revPBE [18] functionals. The DMC calculations are performed using structures optimized by the LDA, while the structures in DFT calculations are optimized individually. For H, O and F adatoms on the preferred adsorption sites, no significant differences in the structural parameters are found using different functionals. The C–H, C–O and C–F bond length values obtained using the LDA are 1.13 Å (top), 1.44 Å (bridge), and 1.50 Å (top), respectively, while the bond length values obtained using the PBE functional are 1.13 Å, 1.47 Å, and 1.56 Å, respectively. The difference between LDA and PBE bond lengths is therefore smaller than 0.01, 0.03 and 0.06 Å for C–H, C–O and C–F, respectively. The numbers of the spin-up and spin-down electrons per cell differ by one for F and H, giving rise to a magnetic ground state. In all three cases, the most stable sites predicted by the DMC and DFT calculations are in agreement, with O favoring the bridge site and F and H favoring the top site. However, the DMC results indicate that DFT calculations vastly overestimate the binding of O and F atoms on graphene. The PBE, PW91, and revPBE results for the adsorption energies are generally better than the LDA results, but they significantly underestimate the energy difference between various sites of F on graphene, giving rise to questionable results for the diffusion barrier. The situations are much better for H on graphene, for which PBE, PW91, and revPBE calculations give quite accurate adsorption energies for the top and hollow sites. Overall, the revPBE adsorption energies are closest to the DMC results, although the relative energy order at different adsorption sites obtained by different functionals is generally consistent. The absolute error of the LDA energies for adsorbed O and F is of the order of 2 eV and 1.5 eV, respectively. Most of the energy errors of the PBE results come from the solid calculation, since PBE is able to give quite accurate atomic energies, as shown in the appendix. The revPBE errors are noticeably smaller compared with LDA results, about 0.2 eV and 0.5 eV for O and F at their preferred sites, respectively. In all three cases, the adsorption energy is negative at their lowest-energy site compared with the atomic state, with O having the strongest energy gain by forming bonds with two C atoms at the bridge site.

Zoom In Zoom Out Reset image size

An alternative way to examine the adsorption energetics is to compare the adatom energies on the surface with that of the molecular state. The DFT results may become better due to error cancelation. These comparisons are shown in figure 4 in the same format as in figure 3. Indeed, the energy errors are reduced, but are less consistent. For example, the PBE and revPBE energy results for O at the lowest-energy bridge sites are higher than the DMC value, while the LDA value is lower. In addition, the results show that single O and H atoms on graphene have a higher energy than their molecular state by 0.6 eV and 1.5 eV, respectively. Only a single F atom on graphene is more favorable than in the molecule, with an energy gain of 0.7 eV/atom. The detailed energy values are given in table 1.

Zoom In Zoom Out Reset image size

We next investigate the graphene sheet fully loaded by H. At the saturation concentration, H atoms are adsorbed on both sides of the graphene plane. The structure of graphane, which was first named by Sofo et al [14], is illustrated in figure 5. Similar to the previous case of single-atom adsorption on graphene, we need to first determine the special k-point to be used in our DMC calculations. The unit cell of graphane contains two C atoms and two H atoms. The results of our LDA calculations using norm-conserving pseudopotentials indicate that the total energy obtained with a 4 × 4 k-point mesh with an offset of (0.5, 0.5) is very close to that obtained by using a 12 × 12 mesh. The former is then used in the DMC calculations, which is equivalent to using a 4 × 4 supercell with a special k-point of (0.5, 0.5). The total number of electrons included in the DMC calculations for graphane is 160.

Zoom In Zoom Out Reset image size

Figure 6 shows the DFT and DMC results for the formation energy of graphane defined as E_f = E(graphane) − [E(graphene) + 2E(H atom)] as a function of the in-plane lattice constant a. The equilibrium lattice constants from LDA, PBE and DMC calculations determined by the fitted curves are 2.50, 2.52 and 2.50 Å, respectively. The DMC energies have a statistical error of 50 meV, which will not affect the lattice constant result. We have also examined the effect of using different one-particle wavefunctions in constructing the Slater determinant, and the difference in the DMC energies using LDA and PBE orbitals is completely negligible. The DMC calculation gives a lattice constant of 2.50 Å, which is similar to the results from our DFT calculations and previous DFT work [14]. However, this value is larger than the reported experimental value of 2.42 Å in a TEM measurement [13]. The DFT calculations with the LDA or PBE exchange–correlation functional give similar lattice constants, with the PBE result slightly larger. Our DMC and DFT results indicate that the lattice constant of graphane is greater than that of graphene (2.46 Å). Therefore, the small lattice constant reported previously [13] remains unexplained, and could arise from the substrate interaction. It is noted that DFT calculations do not provide an accurate value for the formation energy. The errors are about 1 eV for LDA and 0.5 eV for PBE.

Zoom In Zoom Out Reset image size