Van-der-Waals Corrections

by Sven Lubeck & Pasquale Pavone for exciting carbon

Purpose: In this tutorial, you will learn how to perform exciting calculations for materials that are significantly affected by van der Waals (vdW) forces. Since most of the common exchange-correlation functionals fail to describe these long-range dispersion interactions adequately, it is necessary to introduce further corrections. Two rather empirical, but therefore efficient, correction methods are available in exciting. There is, on the one hand, DFT-D2 by Stefan Grimme , and on the other hand, TS-vdW by Alexandre Tkatchenko and Matthias Scheffler . We will use these two methods to investigate the binding between the layers in graphite.

Read the following paragraphs before starting with the rest of this tutorial!

Before starting, be sure that relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts. Here is a list of the scripts which are relevant for this tutorial with a short description.

• SETUP-interlayer-distance.py: Python script for generating structures with different interlayer distances.
• EXECUTE-elastic-strain.sh: (Bash) shell script for running a series of exciting calculations.
• PLOT-energy.py: Python visualization tool for energy-vs-strain curves.
• PLOT-compare-vdW.py: Python visualization tool for multiple energy-vs-distance curves.

From now on the symbol will indicate the shell prompt. Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries. Important note: all input parameters that appear in this tutorial are given in atomic units! ### 1. Introduction In order to take vdW interaction into account, we may add an empirical energy correction to the standard DFT total energy (1) \begin{align} E_{\text{DFT}+\text{disp}}=E_{\text{DFT}} + E_{\text{disp}} . \end{align} This energy correction arises from a sum over pairwiseC_6/R^6-potentials (2) \begin{align} E_{\text{disp}}=-\frac{1}{2}s_6 \sum_{A \neq B}\frac{C_6^{AB}}{\left(R^{\!AB}\right)^6}\,f_{\text{dmp}}\!\left( R^{AB}, R^A_0, R^B_0\right) . \end{align} Here,s_6$is a global scaling factor for all dispersion coefficients$C_6^{AB}$. The sum goes over all atoms$A$and$B$that are in the system. In practice, this sum is of course truncated (see cutoff parameter in the input reference).$R^{AB}$denotes the distance between atom$A$and$B$.$f_{\text{dmp}}\!\left( R^{AB}, R^A_0, R^B_0\right)$is a damping function that ensures that the correction is only applied in the long-range domain. This function should turn to zero in case the atomic separation$R^{AB}$is much smaller than the sum of the atomic van der Waals radii$R^A_0$and$R^B_0$. In this region, the$C_6/R^6$-potential in is not a valid description and should therefore be suppressed. In the case that$R^{AB}is much larger than the sum of the individual van der Waals radii, the damping function should turn to one, since the correction should be fully accounted for in this region. A possible choice for a damping function is (3) \begin{align} f_{\text{dmp}}\!\left( R^{AB}, R^A_0, R_B^0\right) = \left\lbrace 1+\exp\left[-d\left( \frac{R^{AB}}{s_{_{\!R6}}\left(R^A_0 + R^B_0 \right)} - 1 \right)\right]\right\rbrace^{-1} \end{align} which is similar to a Fermi-Dirac distribution. It includes two parametersd$and$s_{_{\!R6}}$. The first parameter$d$determines the steepness of the damping function. The second parameter$s_{_{\!R6}}$is a global scaling factor for all atomic van der Waals radii$R^A_0$. The main difference between the two methods DFT-D2 and TS-vdW is that the$C_6$and$R^A_0$coefficients in the TS-vdW scheme depend on the electron density, whereas in DFT-D2 they are static parameters which do not depend on the chemical environment of the considered atom. ### 2. Graphite: Set up the calculations ##### i) Preparation of the input file First of all, create a directory for each system that you want to investigate. $ mkdir graphite-vdW
$cd graphite-vdW  Secondly, create an exciting input file called input.xml for the system under investigation. As an example, you can find an input file for graphite below. <input> <title>Graphite</title> <structure speciespath="$EXCITINGROOT/species">
<crystal scale="1.88973" stretch="2.46 2.46 6.60">
<basevect>0.5 0.8660254038 0.0</basevect>
<basevect>1.0 0.0000000000 0.0</basevect>
<basevect>0.0 0.0000000000 1.0</basevect>
</crystal>
<species speciesfile="C.xml" rmt="1.2">
<atom coord="0.00000000 0.00000000 0.0"></atom>
<atom coord="0.00000000 0.00000000 0.5"></atom>
<atom coord="0.66666667 0.66666667 0.0"></atom>
<atom coord="0.33333333 0.33333333 0.5"></atom>
</species>
</structure>

<groundstate
rgkmax="6.0"
gmaxvr="20"
ngridk="10 10 4"
xctype="GGA_PBE"
vdWcorrection="DFTD2">
</groundstate>

</input>


Do not forget to insert the correct path for the species files. For that, simply use the following command.

$SETUP-excitingroot.sh  Please note that the values of the attributes rgkmax and ngridk do not correspond to accurate calculations. These values are chosen to reduce the computing time. In order to obtain meaningful results, these parameters must be converged with respect to the quantity of interest. In the input file above, we use DFT-D2 as a vdW correction. This is controlled by the attribute vdWcorrection. There are three possible values for this attribute, namely "DFTD2", "TSvdW", and "none", where the last one is the default value. ##### ii) Generation of input files for different interlayer spacings In order to generate input files for a series of interlayer distances, you have to run the script SETUP-interlayer-distance.py. Notice that the script SETUP-interlayer-distance.py always generates a working directory containing input files for different interlayer distances. Results of the current calculations will also be stored in the working directory. The directory name can be specified by adding the name in the command line. $ SETUP-interlayer-distance.py DIRECTORYNAME


If no name is given, the script uses the default name workdir. Very important: The working directory is overwritten each time you execute the script SETUP-interlayer-distance.py. Therefore, choose different names for different calculations.

The script SETUP-interlayer-distance.py produces the following output on the screen.

$SETUP-interlayer-distance.py DFTD2 Enter minimum interlayer distance dmin [Bohr] >>>> 5 Enter maximum interlayer distance dmax [Bohr] >>>> 8 Enter the number of distances in [dmin,dmax] >>>> 10 Enter interlayer distance at infinity [Bohr] >>>> 20  In this example, the working directory is named DFTD2, the (on screen) input entries are preceded by the symbol ">>>>". The entry values must be typed on the screen when requested. Here, we generate input files for 11 different structures. 10 of these have interlayer separations between 5 Bohr and 8 Bohr. The last structure (the 11th one) has an interlayer separation of 20 Bohr. This last calculation provides a reference energy for all other calculations. Indeed, the interlayer separation for the reference energy would be the same as for an isolated monolayer and therefore infinite. However, the corresponding calculation would not be feasible due to the enormous computational costs. Therefore, convergence tests are needed in order to select the appropriate interlayer separation of the reference structure. The minimum total energies of all structures with respect to the (converged) reference energy can be interpreted as binding energies. In order to obtain meaningful binding energies, then, you should converge the reference energies with respect to the interlayer distance. ### 3. Graphite: Perform the calculation To execute the series of calculation with input files created by SETUP-interlayer-distance.py you have to run the script EXECUTE-elastic-strain.sh. If a name for the working directory has been specified, then you must give it here, too. $ EXECUTE-elastic-strain.sh DFTD2

===> Output directory is "DFTD2" <===

Running exciting for file input-01.xml ----------------------------------

...

Run completed for file input-oo.xml -------------------------------------

$  After the complete run, move to the working directory DFTD2. $ cd DFTD2


Inside this directory, results of the calculation for the input file input-i.xml are contained in the subdirectory rundir-i where i runs from 01 to the total number of structures. For the reference calculation at "infinite" interlayer separation i takes the value oo. The data for binding energy curves are contained in the file energy-vs-strain.

### 4. Graphite: Binding energy curve

Inside the directory where the file energy-vs-strain is located, you can run the script PLOT-energy.py. This will produce two files named PLOT.ps and PLOT.png which both contain the following figure. If the the reference calculation in directory rundir-oo is already performed, the script will also produce the file normalized-energy which contains the effective binding energies with respect to the reference energy from rundir-oo.

### 4. Graphite: Compare correction methods

In order to obtain a binding energy curve for the TS-vdW correction, simply replace the following line in the input file input.xml

 vdWcorrection="DFTD2">


by

 vdWcorrection="TSvdW">


and repeat the procedure from above.

$SETUP-interlayer-distance.py TSvdW Enter minimum interlayer distance dmin [Bohr] >>>> 5 Enter maximum interlayer distance dmax [Bohr] >>>> 8 Enter the number of distances in [dmin,dmax] >>>> 10 Enter interlayer distance at infinity [Bohr] >>>> 20$ EXECUTE-elastic-strain.sh TSvdW

===> Output directory is "TSvdW" <===

Running exciting for file input-01.xml ----------------------------------

...

Run completed for file input-oo.xml -------------------------------------

$cd TSvdW$ PLOT-energy.py

You will obtain the following figure. You can also repeat this procedure for the uncorrected PBE case, i.e., replace the following line in the input file
 vdWcorrection="TSvdW">


by

 vdWcorrection="none">


or just delete this attribute.

If you perform this series of calculations in a directory named PBE-uncorrected, you can create a figure with all three binding energy curves by running the script PLOT-compare-vdW.py as follows.

\$ PLOT-compare-vdW.py normalized-energy DFTD2 TSvdW PBE-uncorrected


The output file PLOT.png will contain the following figure. Bibliography
1. Grimme, Stefan. "Semiempirical GGA-type density functional constructed with a long-range dispersion correction." Journal of Computational Chemistry 27.15 (2006): 1787-1799.
2. Alexandre Tkatchenko and Matthias Scheffler. "Accurate molecular van-der-Waals interactions from ground-state electron density and free-atom reference data." Phys. Rev. Lett. 102, 073005 (2009).