Application of Van Der Waals Density Functionals to Two Dimensional Systems Based on a Mixed Basis Approach
Chung-Yuan Ren Yia-Chung Chang Chen-Shiung Hsue
physics.comp-ph cond-mat.mes-hall



 ^a Department of Physics, National Kaohsiung Normal University,
Kaohsiung 824, Taiwan  	

^b Research Center for Applied Sciences, Academia Sinica, Taipei 115, Taiwan 

^c Department of Physics, National Cheng-Kung University, Tainan 701, Taiwan

^d Department of Physics, National Tsing Hua University,
Hsinchu 300, Taiwan 

† E-mail address: cyren@nknu.edu.tw 



A van der Waals (vdW) density functional was implemented in the mixed basis 
approach previously developed for studying two dimensional systems, in which 
the vdW interaction plays an important role. 
The basis functions here are taken to be the localized B-splines 
for the finite non-periodic dimension and 
plane waves for the two periodic directions. 
This approach will 
significantly reduce the size of the basis set, especially for large systems, 
and therefore is computationally 
efficient for the diagonalization of the Kohn-Sham Hamiltonian.  
We applied the present algorithm to calculate the binding energy 
for the two-layer graphene case and the results are consistent with data 
reported earlier.    
We also found that, due to the relatively weak vdW interaction, 
the charge density obtained self-consistently for the whole bi-layer 
graphene system is not significantly different from the simple 
addition of those for the two individual one-layer system, except when the
interlayer separation is close enough that the  
strong electron-repulsion dominates. This finding suggests an efficient way to 
calculate the vdW interaction for large complex systems involving the 
Moiré pattern configurations.  

PACS: 71.15.Mb, 73.20.-r	


 Application of Van Der Waals Density Functionals to Two Dimensional Systems 
  Based on a Mixed Basis Approach 
    Chung-Yuan Ren^a,†, 
 Yia-Chung Chang^b,c, and Chen-Shiung Hsue^d  
    December 30, 2023
================================================================================================================




§ INTRODUCTION

The electronic properties of two-dimensional (2D) systems are fundamentally
different from those in higher dimensions due to their unusual collective
excitations. Among these 2D materials, graphite is the most well known. 
Graphite has a layered planar structure and 
is electrically conductive along the planes, whereas diamond, another allotrope
of carbon, is an insulator.  
Graphene is an isolated sheet of graphite and can be stacked via 
the weak van der Waals (vdW) interaction to form the graphite structure. 
There has been increasing interest in vdW graphene-based composite systems,
e.g., alkali metal/graphite adsorption systems <cit.> or 
MoS2/graphene heterostructures <cit.>.   
Expectations concerning the creation  
of improved functional electrodevices with 
better performance characteristics are rising from 
the intensive exploration of such graphene-based materials.






First-principles methods based on the density functional theory 
within local density approximation (LDA) and 
generalized gradient approximation(GGA) have proven to
be powerful and successful in investigating static and dynamic 
properties of materials with strong ionic, covalent and metallic interactions. 
Unfortunately, these methods fail to describe the weak vdW 
dispersion interaction properly. 
For example, GGA calculations show no relevant binding 
between graphite sheets <cit.>.
While the LDA approach predicts an underestimated minimum for 
graphite <cit.>-<cit.>, it cannot 
capture the vdW physics <cit.>. 
Neither of these traditional functionals has basis to 
address issues of transferability for soft-matter problems involving the weak 
vdW bonding. To remedy this situation, a new approach using a van der Waals
density functional (vdW-DF) with a nonlocal correlation energy 
has been developed by Dion et al <cit.>. 
This formalism accounts for the dominant dispersion energy, 
which is not correctly treated in standard DFT functionals. 

In this work, we implement the vdW-DF functionals 
using the mixed basis approach developed previously <cit.>. 
The basis functions here are taken to be the localized B-splines 
for the finite non-periodic dimension and 2D
plane waves for the two periodic directions. 
B-splines are highly localized and piecewise polynomials 
<cit.>, which  
have proven to be an excellent tool for the description of 
wavefunctions in a real-space approach <cit.>-<cit.>. 
Such a mixed-basis method <cit.> avoids 
the use of  artificial vacuum layers of large thickness 
introduced by the supercell modeling,  
reducing  significantly the number of the basis functions, and   
therefore easing the computational burden for the 
diagonalization of the Kohn-Sham Hamiltonian. Another advantage of the 
present mixed basis method is that, 
for charged systems, the spurious Coulomb interaction between the defect, 
its images and the compensating background charge in the supercell approach can
be automatically avoided. No further modification needs to be made in the 
total-energy calculation <cit.>.
 

We tested the present algorithm by studying the binding energy between 
two graphene sheets stacked in both AA and AB types, as depicted in  
Fig. <ref>. In addition, the charge density 
around each graphene sheet was found not to be significantly affected by the 
existence of another one, except when the two sheets are so close 
that the electron distributions of individual sheets overlap with each other and 
the strong electron repulsion dominates. 
This revelation would allow us to calculate the binding energy by the 
rigid-density model, i.e., the whole charge density for the system is simply 
the sum of those self-consistently calculated 
for the individual layers, instead of using the very 
time-consuming self-consistent density calculation for the whole system.  
The justified rigid-density model enables a simpler yet accurate evaluation of  
the vdW interaction for large complex systems including the 
Moiré pattern configurations.  
The results will be presented and discussed in details. 






§ METHOD OF CALCULATION



 §.§ B-splines

For the sake of completeness, we
first briefly summarize the B-spline formalism.
More details can be found in Refs. <cit.> and <cit.>.
In general, B-spline of order κ consists
of positive polynomials of degree κ-1, over κ
adjacent intervals. These
polynomials vanish everywhere outside the subintervals τ_i < z <
τ_i+κ.
The B-spline basis set of order κ with the knot sequence {τ_i}
is generated by the following relation :

    B_i,κ(z)=z-τ_i/τ_i+κ-1-τ_iB_i,κ-1(z)+τ_i+κ-x/τ_i+κ
    -τ_i+1B_i+1,κ-1(z),

with

    B_i,1(z)= {[              1, τ_i ≤ z < τ_i+1;              0,    otherwise  . ].

The first derivative of the B-spline of order κ is given by

    d/dzB_i,κ(z)=κ-1/τ_i+κ-1-τ_iB_i,κ-1(z)-κ-1/τ_i+κ-τ_i+1B_i+1,κ-1(z).


Therefore, the derivative of B-splines of order κ is simply a 
linear combination of B-splines of order κ-1, which is
also a simple polynomial and is continuous across the knot sequence.
Obviously, B-splines are flexible to accurately represent any localized 
function of z with a modest number of the basis by only
increasing the density of the knot
sequence where it varies rapidly.



 §.§ vdW-DF functional

The nonlocal energy functional proposed  by 
Dion et al. <cit.> is  

    E^vdW-DF_xc= 
    E^revPBE_x+
    E^LDA_c+
    E^nl_c.
	
The first two parts are simply revPBE exchange <cit.> and 
LDA correlation <cit.>.  
E^nl_c is a non-local correlation functional 
that was introduced to accounts for dispersion interactions, and is  
given as 

    E^nl_c= 1/2∫∫ d r_1 d r_2
    n( r_1) 
    
    ϕ(q_1,q_2,r_12) 
    n( r_2),

where r_12=|r_1-r_2|, and 
q_1,q_2 are the values of a universal function q_0 at 
r_1 and r_2.  
It turns out that the kernel ϕ depends on 
r_1 and r_2 only through two variables 
d_1=q_1r_12 and d_2=q_2r_12, 
and can be expressed as 

    ϕ(d_1,d_2)=2/π^2∫_0^∞a^2da	
    	      ∫_0^∞b^2db  W(a,b)
    		T(ν(a),ν(b),ν^'(a),ν^'(b)),

where W and T are defined as 

    W(a,b)   =   2[(3-a^2)bcos bsin a+
                            (3-b^2)acos asin b+	
           (a^2+b^2-3)sin asin b-3abcos acos b]/a^3b^3, 	
    
    
    T(w,x,y,z)   =   1/2[1/w+x+ 1/y+z] ×
    [1/(w+y)(x+z)+ 1/(w+z)(y+x)].

The quantities ν and ν^' are given by
ν(u)=u^2/2h(u/d_1) and 
ν^'(u)=u^2/2h(u/d_2) with 
h(t)=1-exp (-4π t^2/9).

The universal function q_0 reads as 

    q_0( r)=-4π/3ϵ_xc^LDA n( r) - Z_ab/9
    s^2( r)
    k_F( r).

Here, the Fermi wave vector k_F and the reduced gradient s  are  

    k^3_F ( r)=3π^2n( r),  
    s( r)=∇ n( r)/
    2k_F( r)n( r),

and 
Z_ab=-0.8491.

The nonlocal correlation energy in Eq. (<ref>) is expressed as a double
spatial integral. To alleviate the O(N^2) evaluation of such the integral, 
we adopt the algorithm by 
Román-Pérez and Soler <cit.>, 
which transforms the double real space integral to reciprocal space 
and reduces the computational effort.

First, ϕ was interpolated as 

    ϕ(q_1,q_2,r_12)=∑_α,βϕ(q_α,q_β,r_12)
    p_α(q_1) p_β(q_2),

where q_α are fixed  values, chosen to ensure a good 
interpolation of function ϕ. 
Here, we use cubic splines interpolation, in which p_α(q) 
is a succession of cubic polynomial in every interval 
[q_β,q_β+1], matching in value and the first two derivatives at 
every point q_β.

Substituting Eq. (<ref>) into (<ref>),

    E^nl_c= 1/2∑_α,β∫∫ d r_1 d r_2θ_α( r_1) 
    θ_β( r_2) 
    ϕ_αβ(r_12),

with θ_α( r)=n( r)p_α(q_0( r)) and 
ϕ_αβ(r_12)≡ϕ(q_α,q_β,r_12).
Now, with use of the convolution, just like the Coulomb energy, 
E^nl_c becomes

    E^nl_c= 1/2∑_α,β∫ d kθ_α^*( k) 
    θ_β( k) 
    ϕ_αβ(k),

where θ_α( k) and ϕ_αβ(k) are the 
corresponding Fourier transforms. In practice, 
 ϕ_αβ(k) was pre-calculated 
in spherical radial mesh of points k.
Then, ϕ_αβ(k) and its second derivative via 
cubic spline interpolation were stored for later use.
A logarithmic mesh of interpolation points q_α, 
of which the total number is 20 in the present calculation, 
was used to describe ϕ up to a cutoff q_c of 5.0 a.u..



 §.§ GGA charge density in the mixed-basis approach 

With one set of B-splines for the non-periodic z 
direction and 2D plane waves for the periodic xy plane,
the present mixed basis used to expand the wavefuction
is defined as

    <  r  |  k_∥ +  G_∥ ; j, κ > =
    1/√(A)e^i(  k_∥ +  G_∥ ) ·ρ B_j,κ(z),

where
G_∥
denotes an in-plane reciprocal lattice vector and   
k_∥ is the in-plane Bloch wave vector.
A is the surface area of the system.
Therefore, the charge density can be written in the form  







    n( r)  = ∑_ g n( g,z)
      e^ i  g·ρ ,

where g=  G_∥- G^'_∥.

Because the GGA energy functional 
depends upon |∇ n( r)|, the 
corresponding potential v_xc is a functional of not only  
|∇ n|, but also of 
∇^2 n and  ∇ n ·∇ |∇ n |. In order to 
efficiently and precisely obtain v_xc,   
we used the method by White and Bird <cit.>. 
First of all, we interpolate 
n( g,z) along the z direction by
using the Fourier interpolation technique:  

    n( g,z)  = ∑_g_z n( g,g_z)
      e^ i g_z z  .

Then,  

    n( r)=∑_ G   n (  G) 
    e^i  G· r,

where G is a compact notation for ( g,g_z).

Following the procedure in Ref. <cit.>, 

    ∇ n( r)=∑_ G  i  G  n (  G) 
    e^i  G· r=
    1/N∑_ G,R  i  G  n (  R) 
    e^i  G· ( r-R)



    =1/N∑_ G  i  G ( 
    ∑_ R n (  R) 
    e^-i  G· R) 
    e^ i  G· r,

with N real space points R of the fast Fourier-transform (FFT) grid set. 

We define f_xc such that 

    E_xc[n]=∫ f_xc( n( r), |∇ n( r)|)d r.

E_xc[n] can be approximated by 

    E_xc[n] ≈V/N∑_ R  f_xc( n( R), |∇ n( R)|).





The associated xc potential at the FFT grid point R
can be obtained efficiently through

    v_xc( R)    =   N/Vd E_xc/d n( R)
       =   ∂ f_xc/∂ n ( R)
    +∑_ R^'∂ f_xc/∂∇ n ( R^')·d ∇ n ( R^')/d n ( R)
       =   ∂ f_xc/∂ n ( R)
    + 1/N∑_ G,R^' i  G ·∂ f_xc/∂∇ n ( R^')
    e^i  G· ( R^'-R)
       =   ∂ f_xc/∂ n ( R)
    + 1/N∑_ G,R^' i  G ·∇ n ( R^')/n ( R^')∂ f_xc/∂ |∇ n ( R^')|
    e^i  G· ( R^'-R).

Given the charge density on the FFT grid points, only eight FFT's are required
to obtain v_xc. That is computationally moderate with respect to 
the derivation of the
second derivative needed to evaluate the conventional potential via  

    v_xc( r)= ∂ f_xc/∂ n ( r)-
    ∇·∂ f_xc/∂∇ n ( r).





 §.§ vdW-DF total energy
 
For the vdW-DF total energy functional, we first performed the self-consistent
total energy calculation using the GGA-PBE functional <cit.>. With the 
converged charge density obtained in the previous step,  
the revPBE exchange energy <cit.> and LDA correlation energy <cit.>
 are evaluated and substituted for 
the GGA-PBE counterparts, and the nonlocal correlation energy E^nl_c 
is added. Now, the vdW-DF energy functional is written 

    E^vdW-DF= E^PBE- E^PBE_xc+  (
    E^revPBE_x+
    E^LDA_c+
    E^nl_c  ).

The last three terms in the above equation are 
treated as a post-GGA perturbation because of their low sensitivity to 
the choice of GGA electronic density. By the present approach, E^nl_c was obtained 
via 

    E^nl_c= 1/2A
    ∑_α,β∑_ g∫ dk_zθ_α^*( k)
    θ_β( k)
    ϕ_αβ(k),

with k=( g,k_z).


The interplanar binding energy per surface atom E_b is defined as 

    E_b = (E_bilayer-2E_graphene)/N_b

where E_bilayer and E_graphene  
are respectively the total energy of the bilayer graphene system and   
that of the system containing only one graphene sheet. 
N_b is the number of atoms in one single graphene sheet. 




§ APPLICATIONS OF PRESENT METHOD
 

To test the present approach, we apply it to investigate the
vdW interaction between graphene sheets. 
The calculations were carried out with
the unit cell containing two graphene sheets. 
Here, we study both AA and AB stacking, as depicted in Fig. <ref>. 
 In the former case,
the carbon atoms of the adjacent sheets are aligned directly on top 
of each other. In the latter case, 
the energetically more stable structure, the 
graphene layers are shifted relative to each other such that 
half of the atoms are located exactly over the center of a hexagon and 
another half lie directly on top of the atoms in the second graphene sheet. 

All C atoms in the graphene sheet were kept at the ideal positions. The
in-plane lattice constant a_0 was fixed to 
the experimental value of 2.461 Å, while the inter-plane distance d 
is allowed to vary. 
A FFT mesh with the grid spacing of  
0.08 Å for the charge density are chosen for accurate total
energy calculations. A mixed basis set with 
34 B-splines distributed over a maximum range of 6.0 a_0 
and 2D plane waves with an energy cutoff of 30 Ry 
are used to expand the wavefunction.   
The 7 × 7 Monkhorst-Pack grids including Γ point
were taken to sample the surface Brillouin zone. 
We used the Vanderbilt's ultra-soft pseudopotential (USPP) <cit.>.
The C USPP was generated from
the Vanderbilt's code <cit.> and its quality was   
examined previously <cit.>. 
The potential is determined self-consistently until its
change is less than 10^-7 Ry. Finally, the vdW-DF 
total energy and the associated binding energy are calculated according to the 
procedure described in Sec. <ref>.
For comparison, we also performed calculations by using 
the standard supercell approach implemented in the 
popular VASP code with the projector-augmented-wave potential 
(PAW) <cit.>.  
A typical vacuum space of 10 Å  required in VASP was 
used in the calculation.






The binding energy of bilayer graphene in the AA stacking 
as a function of interlayer 
separation d is shown in Fig. <ref>.  
The VASP counterpart is also plotted for comparison. 
Obviously, the GGA-PBE calculations show no relevant minimum for the graphite 
binding energy, reflecting the failure to include the proper long-range 
dispersive interaction within the GGA approximation. On the other hand, 
with the vdW-DF xc functional expressed in Eq. (<ref>), 
we obtained for the graphene pair a binding energy 
of 47 meV/atom for the AA stacking. 
More importantly, the results obtained with our algorithm agree nicely with those by 
the popular VASP code. We are then convinced that the present program 
has been implemented successfully for the vdW interaction and the 
outcomes are very reliable. 

We also calculated for the AB stacking and    
the resulting binding energy is displayed in Fig. <ref>. 
Clearly, the AB stacking is  
energetically more stable than the AA stacking, in agreement with the 
experiment that natural graphite occurs mainly with AB stacking order <cit.>. 

The binding energy for the AB stacking was found to be 50.5 meV/atom at the distance 
of 3.7 Å  and that for the less favored AA stacking was 47.0 meV/atom at the distance
of 3.8 Å.  The results are consistent with data reported in the 
literature <cit.>,<cit.>-<cit.>. It can be seen from  
this figure that the E_b curve for the AB stacking merges into that for the AA
stacking at large separation d, which seemingly indicates that 
the total energy will be less sensitive to the orientation of the two graphene
sheets if the separation is not too close.    

Actually, contrary to the covalent bond, 
the charge density distribution near the individual graphene sheet would not 
be noticeably affected via 
the weak vdW interaction from other sheets unless the sheet separation  
is close enough such that it begins to 
'contact' or even overlaps with the density from the adjacent sheet. In that case, 
the shape of the density distribution around the graphene sheet
will be distorted because of the dominant strong electron repulsion.   

To justify this assertion, we use
a rigid-density model, i.e., the whole charge density of the system in the AA stacking 
is simply assumed to be the sum of those self-consistently calculated 
for the individual single layer. Then we employed such charge density to
re-calculated the total energy E^tot,AA_rigid and compared it 
to the total energy E^tot,AA_scf with self-consistent density calculations. 
The results are summarized in Table <ref>. Here, we chose seven cases with 
various interlayer separation d, as also indicated in Fig. <ref>.
Clearly, the total energy per surface atom 
by the rigid-density model is very similar to the self-consistent total energy.  
The energy difference only becomes notable (2.7 meV/atom)  
for Case 7 with d=2.95 Å, which already enters the electron-repulsive
region. With detailed analysis of energy components, we found that even 
though the energy difference from the kinetic energy part is somewhat 
sizable (for example, 32 meV/atom for Case 4 with d=3.69 Å), but it is 
largely compensated for with the Hartree energy part, leading to 
almost the same value in total. 

When the charge density of the single graphene sheet 
was shifted a bit to the second one, we obtained similar conclusion, as 
shown in Table <ref> for the AB stacking. Therefore, it is reasonable to 
expect that if the first graphene sheet is rotated with respect to the 
second to become a Moiré rotated pattern, the rigid-density model can still hold to  
efficiently predict a reliable binding energy for such a 
large complex system. This finding is promising in searching for the true
ground-state atomic configuration for vdW-dominated graphene-based materials,
like MoS2/graphene heterostructures <cit.>. Instead of using the  
very time-consuming self-consistent approach for such materials, 
the rigid-density model allows us to only focus on accurate charge density calculations
of every individual slab of different type.  





§ CONCLUSIONS
 
In conclusion, we have successfully implemented
the van der Waals (vdW) density functional proposed by Dion et al. <cit.>
in our mixed-basis approach for investigating the bi-layer graphene system. 
As compared to the conventional supercell model with
alternating slab and vacuum regions, it is a real space approach along
the non-periodic direction. Therefore, the number of the 
basis functions used to expand the wavefunction is significantly reduced, 
especially for large complex systems. 

We also found that the self-consistent total energy 
obtained for the bilayer 
system is not significantly different from that with charge density 
assumed to be the simple 
sum of those for the two individual single-layer system, except when the
distance between the two layers is close enough that the  
strong electron-repulsion dominates. Such observations help us to propose a
rigid-density model which can efficiently 
calculate the binding of vdW-dominated 2D systems with  Moiré pattern configurations.  

This work was supported by Ministry of Science and Technology under 
grant numbers MOST 106-2112-M-017 -003 and MOST 106-2112-M-001-022 
and by National Center for Theoretical Sciences of Taiwan.


* 


§ 


99
BKLW
M. Breitholtz, T. Kihlgren, S. -Å. Lindgren, and L. Walldén, 
Phys. Rev. B 66, 153401 (2002).
ZKSH
E. Ziambaras, J. Kleis, E. Schröder, P. Hyldgaard,
Phys. Rev. B 76, 155425 (2007).
PHABNP
D. Pierucci,H. Henck, J. Avila, A. Balan, C. H. Naylor, G. Patriarche,
Y. J. Dappe, M. G. Silly, F. Sirotti, A. T. C. Johnson, M. C. Asensio,
and A. Ouerghi,
Nano Lett. 16, 4054 (2016). 
RJHSLL
H. Rydberg, N. Jacobson, P. Hyldgaard, S. I. Simak, B. I. Lundqvist, and 
D. C. Langreth,
Surf. Sci. 532-535, 606 (2003). 
LM
I. -H. Lee and R. M. Martin, 
Phys. Rev. B 56, 7197 (1997).
B
J. C. Boettger, Phys. Rev. B 55, 11202 (1997).
NPA
D. Nabok, P. Puschnig, and C. Ambrosch-Draxl, 
Comp. Phys. Comm. 182, 1657 (2011).
DRSLL
M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist,
Phys. Rev. Lett. 92, 246401 (2004), 95, 109902(E) (2005). 
RHC
C. Y. Ren, C. S. Hsue and Y.-C. Chang, Comp. Phys. Comm. 188, 94 (2015).
RCH
C. Y. Ren, Y.-C. Chang, and C. S. Hsue, Comp. Phys. Comm. 202, 188 (2016).
deBoor
Carl deBoor, A practical Guide to Splines, (Springer, New York, 1987).
JBS
W. R. Johnson, S. A. Blundell, and J. Sapirstein, Phys. Rev. A 37, 307 (1988).
JH
H. T. Jeng, and C. S. Hsue, Phys. Rev. B 62, 9876 (2000).
RJH
C. Y. Ren, H. T. Jeng, and C. S. Hsue, Phys. Rev. B 66, 125105 (2002).
LC
G.-W. Li and Y.-C. Chang, Phys. Rev. B 48, 12032 (1993).
LC1
G.-W. Li and Y.-C. Chang, Phys. Rev. B 50, 8675 (1994).
ZY
Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998). 
PW
 J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992) and references therein.
RS
G. Román-Pérez, and J. M. Soler, Phys. Rev. Lett. 103, 096102 (2009). 
WB
J. A. White and D. M. Bird, Phys. Rev. B 50, 4954(R) (1994).
PBE
J. P. Perdew, K. Burke, and  M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 
DV
D. Vanderbilt, Phys. Rev. B 41, 7982 (1990).
vancode
http://www.physics.rutgers.edu/ dhv/uspp/.
KJ
G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
KF
G. Kresse and J. Furthmüller, Comput. Mater. Sci.  6, 15 (1996).
THZ
P. H. Tan, W. P. Han, W. J. Zhao, Z. H. Wu, K. Chang, H. Wang, Y. F. Wang,
N. Bonini, N. Marzari, N. Pugno, G. Savini, A. Lombardo, A. C. Ferrari,
Nat. Mater. 11, 294 (2012).
CSLL
S. D. Chakarova-Käck, E. Schröder, B.I. Lundqvist, and D.C. Langreth, 
Phys. Rev. Lett. 96, 146107 (2006). 
CKHS
S. D. Chakarova-Käck, J. Kleis, P. Hyldgaard, and E. Schröder, 
New J. Phys. 12, 013017 (2010). 






FIGURE CAPTIONS    

Fig. 1: Atomic structure of the graphite in both AA and AB stacking.  
 

Fig. 2: (Color online) Binding energy for the two-layer graphenes in AA stacking 
        calculated with
        GGA-PBE and vdW-DF xc functionals. The counterparts obtained by VASP are 
	also displayed for comparison.  
 

Fig. 3: (Color online) Binding energy for the two-layer graphenes in both AA and 
        AB stacking calculated with vdW-DF xc functional.  
 










