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
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
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
K-Points
0
Monkhorst Pack
15 15 1
0 0 0
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!
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¶
- Are effects of Van der Waals interaction well-estimated in the framework of DFT, when using GGA?
- How does the Tkatchenko-Scheffler method account for Van der Waals interaction?
- Is it meaningful to compare the total energy of two structures per unit cell or per atom?
- Is it possible to avoid that CHGCAR and WAVECAR are written?
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!
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
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
K-Points
0
Monkhorst Pack
15 15 7
0 0 0
Pseudopotential of C.
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.
from IPython.display import Image
PATH = "./e11_interlayer-distance/"
Image(filename = PATH + "loop.png", width=800)
- How could the KPOINTS file be changed to completely neglect interactions between layers, even though this is clearly not desired here?
- What is the motivation to use Ewald summation and what is the basic idea?