VASP Hands-on: Bulk systems

Part 4: Van der Waals forces


10 Interlayer binding in graphite using the Tkatchenko-Scheffler method

By the end of this tutorial, you will be able to:

  • include the effect of Van der Waals interaction in density-functional-theory (DFT) calculations
  • set LWAVE and LCHARG
  • compute the difference between the total energy of two structures

10.1 Task

Determine the interlayer binding energy of graphite using the Tkatchenko-Scheffler method to account for Van der Waals interactions.

Semilocal DFT at the level of the generalized gradient approximation (GGA) underestimates long-range interactions. In the case of graphite, PBE predicts the interlayer binding energy of ~1 meV/atom which is too small compared to the RPA reference of 0.048 eV/atom S. Lebègue et al., Phys. Rev. Lett. 105, 196401 (2010). This calls for a different density functional of the GGA type for Van der Waals interactions. The Tkatchenko-Scheffler method is based on Becke's power‐series ansatz from 1997 and is explicitly parameterized by including damped atom‐pairwise dispersion corrections of the form $C_6 \cdot R^{−6}$.

For details of the implementation of the Tkatchenko-Scheffler method in VASP and a number of tests see Bucko et al., Phys. Rev. B 87, 064110 (2013).

10.2 Input

The input files to run this example are prepared at $TUTORIALS/bulk/e10_VdW-TS-binding. Check the subdirectories!

graphene/POSCAR


graphene
1.0
1.22800000 -2.12695839  0.00000000
1.22800000  2.12695839  0.00000000
0.00000000  0.00000000  20.
2
direct
   0.00000000  0.00000000  0.25000000
   0.33333333  0.66666667  0.25000000

graphite/POSCAR


graphite
1.0
1.22800000 -2.12695839  0.00000000
1.22800000  2.12695839  0.00000000
0.00000000  0.00000000  6.71
4
direct
   0.00000000  0.00000000  0.25000000
   0.00000000  0.00000000  0.75000000
   0.33333333  0.66666667  0.25000000
   0.66666667  0.33333333  0.75000000

graphene/INCAR and graphite/INCAR


LWAVE  = False
LCHARG = False

ISMEAR = -5
SIGMA  = 0.01

ALGO   = Fast
PREC   = Accurate
EDIFF  = 1e-6

# Van der Waal setting
IVDW       = 20      ! Tkatchenko-Scheffler method
LVDW_EWALD = True    ! (optional) Ewald summation for pairwise interactions

graphene/KPOINTS


K-Points
 0
Monkhorst Pack
 15 15 1
 0  0  0

graphite/KPOINTS


K-Points
 0
Monkhorst Pack
 15 15 7
 0  0  0

graphene/POTCAR and graphite/POTCAR


Pseudopotential of C.


There is no interaction of layers in $z$ direction for graphene, so we need only one $\mathbf{k}$ point in this direction.

Make comments in the INCAR file for each tag and read about the specific setting on the VASP Wiki! LWAVE and LCHARG are set to false, so that the wavefunction and charge density is not written to the WAVECAR and CHGCAR, respectively. As we only use the total energy and do not want to perform a band-structure or DOS calculation on top, these are not needed and we can save some space by not writing them.

10.3 Calculation

Go ahead and run VASP for both, graphene and graphite, in the corresponding subdirectories of this example's directory with

cd $TUTORIALS/bulk/e10_*/graphene
mpirun -np 2 vasp_std
cd $TUTORIALS/bulk/e10_*/graphite
mpirun -np 2 vasp_std

You can visualize a $3\times 3 \times 3$ grid of the structure using py4vasp!

In [2]:
import py4vasp

graphite = py4vasp.Calculation.from_path("./e10_VdW-TS-binding/graphite")
graphite.structure.plot(3)

You can see that the interlayer distance is vastly different. In order to obtain the interlayer binding energy per atom, get the total energy from the OUTCAR files and divide it by the number of atoms per unit cell! The difference between these energies is the interlayer binding energy per atom!

What is the interlayer binding energy per atom of graphite using the Tkatchenko-Scheffler method?

Click to see the answer!

Use the following lines to get the total energy from the OUTCAR files and print the calculation explicitly.

cd $TUTORIALS/bulk/e10_*/
    en1=$(grep "free  ene" graphene/OUTCAR |tail -1|awk '{print $5}')
    en2=$(grep "free  ene" graphite/OUTCAR | tail -1 | awk '{print $5}')
    echo $en2/4 - $en1/2

Finally, this yields -37.42908542/4 - -18.54564785/2 and hence the binding energy is -0.08444743 eV/atom.

Even though the Tkatchenko-Scheffler method predicts a reasonable geometry, see Example 11, it overestimates the energetics strongly: the computed binding energy is too large compared to the RPA reference of 0.048 eV/atom. This overestimation is, at least in part, due to neglecting the many-body interactions. Let us try to rectify the issue by using the many-body-dispersion-energy (MBD) method proposed in Tkatchenko et al., Phys. Rev. Lett. 108, 236402 (2012).

Change the Van der Waals settings for both, graphene and graphite, in the corresponding INCAR file to include many-body corrections and compute the interlayer binding energy per atom! In particular, replace

IVDW       = 20      ! Tkatchenko-Scheffler method
LVDW_EWALD = True    ! (optional) Ewald summation for pairwise interactions

with

IVDW       = 202     ! MBD variant of Tkatchenko-Scheffler method
LVDWEXPANSION = True ! (optional) display 2- to 6- body contributions to dispersion energy
Click to see the answer!

Use the following lines to get the total energy from the OUTCAR files and print the calculation explicitly.

cd $TUTORIALS/bulk/e10_*/graphene
    
    cp INCAR INCAR.bk
    cat INCAR.bk | grep -v VDW > INCAR
    echo "IVDW       = 202" >> INCAR 
    echo "LVDWEXPANSION = True" >> INCAR
    
    mpirun -np 2 vasp_std
    
    cd $TUTORIALS/bulk/e10_*/graphite
    
    cp INCAR INCAR.bk
    cat INCAR.bk | grep -v VDW > INCAR
    echo "IVDW       = 202" >> INCAR 
    echo "LVDWEXPANSION = True" >> INCAR   
    
    mpirun -np 2 vasp_std
    
    cd $TUTORIALS/bulk/e10_*/
    en1=$(grep "free  ene" graphene/OUTCAR |tail -1|awk '{print $5}')
    en2=$(grep "free  ene" graphite/OUTCAR | tail -1 | awk '{print $5}')
    echo $en2/4 - $en1/2

Finally, this yields -37.44318263/4 - -18.61910927/2 and hence the binding energy is -0.0512410225 eV/atom.

The computed value is now fairly close to the RPA reference of $0.048$ eV/atom Lebègue et al., PRL 105, 195401 (2010).

10.4 Questions

  1. Are effects of Van der Waals interaction well-estimated in the framework of DFT, when using GGA?
  2. How does the Tkatchenko-Scheffler method account for Van der Waals interaction?
  3. Is it meaningful to compare the total energy of two structures per unit cell or per atom?
  4. Is it possible to avoid that CHGCAR and WAVECAR are written?

11 Interlayer distance in graphite

By the end of this tutorial, you will be able to:

  • compute the interlayer distance of a layered compound
  • explain the basic concept of Ewald summation

11.1 Task

Determine the interlayer distance of graphite in the stacking direction using the Tkatchenko-Scheffler method to account for Van der Waals interactions.

The long-range Coulomb potentials that need to be integrated not just but also within the Tkatchenko-Scheffler method can optionally be treated using Ewald summation. Basically, Ewald summation treats the short-range contribution in real space, and the long-range contribution in reciprocal space. This is in order to avoid singularities.

In this example, we use the Tkatchenko-Scheffler method because GGA in the framework of DFT underestimates long-range dispersion interactions such as Van der Waals interactions. This problem causes an overestimation of the lattice constant in graphite along the stacking direction. In particular, GGA-PBE yields $8.84$ Å, while the experimental result is $6.71$ Å.

11.2 Input

The input files to run this example are prepared at $TUTORIALS/bulk/e11_interlayer-distance. Check them out!

POSCAR


graphite
1.0
1.22800000 -2.12695839  0.00000000
1.22800000  2.12695839  0.00000000
0.00000000  0.00000000  c
4
direct
   0.00000000  0.00000000  0.25000000
   0.00000000  0.00000000  0.75000000
   0.33333333  0.66666667  0.25000000
   0.66666667  0.33333333  0.75000000

INCAR


ISMEAR = 0
SIGMA  = 0.1

ALGO   = Fast
PREC   = Accurate
ENCUT  = 620
EDIFF  = 1e-6

# Van der Waal setting
IVDW       = 20      ! Tkatchenko-Scheffler method
LVDW_EWALD = True    ! (optional) Ewald summation for pairwise interactions

KPOINTS


K-Points
 0
Monkhorst Pack
 15 15 7
 0  0  0

graphite/POTCAR


Pseudopotential of C.


In the POSCAR file, the placeholder c should be replaced by trial values until you find the optimal distance of layers.

Recall the meaning of each tag in the INCAR file.

11.3 Calculation

Let's do this in steps! Change to this example's directory:

cd $TUTORIALS/bulk/e11_*

Determine the optimal interlayer distance up to 2 decimal places starting from a guess interval of $[6.6,6.8]$ Å!

Click to see the answer!

loop.sh


rm loop.dat
vasp_rm

for c in $(seq 6.6 0.01 6.8)
do
  cat>POSCAR<<!
graphite
1.0
1.22800000 -2.12695839  0.00000000
1.22800000  2.12695839  0.00000000
0.00000000  0.00000000  $c
4
direct
   0.00000000  0.00000000  0.25000000
   0.00000000  0.00000000  0.75000000
   0.33333333  0.66666667  0.25000000
   0.66666667  0.33333333  0.75000000
!

  mpirun -np 2 vasp_std

  cp OUTCAR OUTCAR.c$c
  energy=$(grep "free  ene" OUTCAR.c$c|awk '{print $5}')

  echo $c $energy >> loop.dat
done

The optimal length of the lattice vector $\mathbf{c}$ normal to the stacking direction is determined in a series of calculations with varied values of $|\mathbf{c}|$, while all other degrees of freedom are fixed at their experimental values. The computed $|\mathbf{c}|$ vs. energy dependence is written in the file loop.dat.

Plot the result with loop.gp


set term png
set output "loop.png"

set title "Graphite"
set xlabel "interlayer distance"
set ylabel "Total energy (eV)"

plot "loop.dat" using 1:2 w lp

Then, open loop.png.

Using the Tkatchenko-Scheffler method, the computed interlayer distance of 6.65 Å agrees well with the experimental value, 6.71 Å. This might be surprising as in Example 9 you have seen that the interlayer binding energy is too large compared to the RPA reference.

In [3]:
from IPython.display import Image

PATH = "./e11_interlayer-distance/"
Image(filename = PATH + "loop.png", width=800)
Out[3]:
No description has been provided for this image
  1. How could the KPOINTS file be changed to completely neglect interactions between layers, even though this is clearly not desired here?
  2. What is the motivation to use Ewald summation and what is the basic idea?

Awesome! You have finished all parts of Session 4!

Go to Top $\uparrow$