Nature and Strength of Interlayer Binding in Graphite

We compute the interlayer bonding properties of graphite using an ab initio many-body theory. We carry out variational and diffusion quantum Monte Carlo calculations and find an equilibrium interlayer binding energy in good agreement with most recent experiments. We also analyze the behavior of the total energy as a function of interlayer separation at large distances comparing the results with the predictions of the random phase approximation.

Figure 1

Binding energy (BE) curve for graphite with $AB$ stacking obtained at the DFT-LDA level of theory, using plane waves (PW) and Gaussian basis sets. The results of PW calculations (solid black line) converged as a function of the kinetic energy cutoff (90 Ry) are in excellent agreement with those of $8s8p4d$ Gaussian basis sets (red triangles), carried out with the same cell (32 atoms and the $\Gamma$ point). Smaller Gaussian basis sets ($4s4p2d$, full circle) yield a much larger BE and a slightly smaller equilibrium distance. Fully converged PW calculations (4-atoms unit cell, 90 Ry, $20\times 20\times 8$ $\mathbf{k}$ points) yield a BE of $24\mathrm{meV}/\mathrm{\text{atom}}$ and an equilibrium interlayer position $D$ of 3.30 Å. Note that $D$ is significantly affected by convergence as a function of $\mathbf{k}$-point sampling, while the value of the BE depends weakly on it.

Figure 2

BE of graphite as a function of separation between planes (a) and of the number of layers ($n$) included in the periodic cell used in our calculations (b). Note the logarithmic scales. In (a) results obtained with VMC and LRDMC calculations are reported. In (b) only LRDMC results are shown.

Figure 3

BE of graphite as obtained at the VMC (black solid circles) and LRDMC (red squares) level of theory, using a $2\times 2\times 1$ supercell, and at the LRDMC level (blue triangles) with a $2\times 2\times 2$ supercell. All data include finite-size corrections using the scheme of Ref. 23. Dotted lines and solid lines are obtained with a fit with the function $a\exp (\delta D)+b/D^4$. The experimental interlayer distance [31] is shown by the dashed line. The large difference in BE between the $2\times 2\times 1$ and $2\times 2\times 2$ supercell calculations arises from large spurious interactions between periodic images of planes in the case of the $2\times 2\times 1$ cell. These spurious interactions are still present in the case of the $2\times 2\times 2$ cell, but they are greatly reduced, as shown by the difference ($\sim 5\mathrm{meV}/\mathrm{\text{atom}}$) of BE obtained with 4 planes in the cell and the extrapolated value.