Van Der Waals Energy Functionals

by Dmitrii Nabok & Stefan Kontur for exciting boron

Purpose: In this tutorial we show how to set up and execute van-der-Waals density-functional (vdW-DF) calculations using the exciting code. In particular, we consider the case of graphite, where the vdW forces are known to play an important role. Structurally, graphite can be viewed as a stack of graphene (a plane of hexagonally packed carbon atoms) layers which are bound by weak vdW forces. Our aim is to find the equilibrium interlayer separation and the corresponding binding energy.

0. Define relevant shell variables and download scripts

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.

  • EXECUTE-graphite.sh: Bash script to perform self-consistent groundstate calculations for a given set of interlayer distances in graphite.
  • EXECUTE-noloco.sh: Bash script to calculate the non-local (vdW) contribution to the correlation energy.
  • PLOT-graphite.sh: Bash script to visualize the results using xmgrace.
  • graphite_xmgr.par: Parameter file for visualization with xmgrace.

From now on the symbol $ will indicate the shell prompt.

Requirements: Bash shell.

Important note: all input parameters that appear in this tutorial are given in atomic units!

1. Compiling "noloco"

To calculate the non-local contribution to the correlation energy, we use the external code noloco. All details concerning implementation and usage of noloco can be found here. Please note that this program is not compiled in exciting by default. To do it, one needs to execute

$ cd $EXCITINGROOT/src/src_vdwdf
$ make

2. Graphite: Preparation of the input file

In this tutorial, we examine the equilibrium interlayer distance in graphite and determine the corresponding binding energy.

The first step is to move back to your general working directory. Here, one should create a directory where this tutorial is executed, e.g., graphite.

$ mkdir graphite
$ cd graphite

Inside this directory, we create (e.g., just copy/paste from below) the file input_graphite.xml corresponding to a template of an exciting input file. This file contains all important parameters to perform self-consistent calculations. The only parameters which are not specified (denoted by 'XXXXX' and 'YYYYY', respectively) are the doubled interlayer distance (the length of the $c$ axis) and the number of grid points along the $c$ direction where the electronic density will be calculated.

<input>
 
   <title>Graphite</title>
 
   <structure speciespath="$EXCITINGROOT/species">
      <crystal scale="1.88973" stretch="2.461 2.461 XXXXX">
         <basevect>0.5000000000 0.8660254038 0.0000000000</basevect>
         <basevect>1.0000000000 0.0000000000 0.0000000000</basevect>
         <basevect>0.0000000000 0.0000000000 1.0000000000</basevect>
      </crystal>
      <species speciesfile="C.xml" rmt="1.3">
         <atom coord="0.00000000 0.00000000 0.00000000"></atom>
         <atom coord="0.00000000 0.00000000 0.50000000"></atom>
         <atom coord="0.66666667 0.66666667 0.00000000"></atom>
         <atom coord="0.33333333 0.33333333 0.50000000"></atom>
      </species>
   </structure>
 
   <groundstate
      do="fromscratch"
      rgkmax="5.0"
      gmaxvr="12"
      ngridk="4 4 2"
      xctype="GGA_PBE"
      outputlevel="high">
   </groundstate>
 
   <properties>
      <chargedensityplot>
         <plot3d>
            <box grid="20 20 YYYYY">
               <origin coord="0 0 0" />
               <point coord="1 0 0"/>
               <point coord="0 1 0"/>
               <point coord="0 0 1"/>
            </box>
         </plot3d>
      </chargedensityplot>
   </properties>
 
</input>

At variance with examples from other tutorials, you do not need explicitly to set the variable $EXCITINGROOT inside this file. This action will be accomplished by the scripts which are used in this tutorial.

The actual input file for the exciting calculation (which is always called input.xml) will be created by the script EXECUTE-graphite.sh. For this tutorial, the parameters, such as rgkmax and ngridk in groundstate or grid in box, have been chosen in order to speed up the calculations and thus do not correspond to completely converged calculations. Note that choosing parameters in order to achieve convergence might lead to a much longer runtime of the calculations, depending on the computers where you execute it.

3. Graphite: Execute the calculations

The total energy and charge-density distribution

In this tutorial, the vdW-DF formalism is applied as a post-scf procedure (following the original paper by Dion and co-workers (PRL, 2004)):

(1)
\begin{align} E_{\rm vdW-DF} = E^{\rm PBE-GGA} - E_{\rm xc}^{\rm PBE-GGA} + \left( E_{\rm x}^{\rm revPBE-GGA} + E_{\rm c}^{\rm LDA}\right) + E_{\rm c}^{\rm NL}. \end{align}

Here $E^{\rm PBE-GGA}$ is the total energy obtained from PBE, $E_{\rm xc}^{\rm PBE-GGA}$ is the exchange-correlation energy from PBE, $E_{\rm x}^{\rm revPBE-GGA}$ is the exchange energy from revPBE, and $E_{\rm c}^{\rm LDA}$ is the correlation energy from LDA. $E_{\rm c}^{\rm NL}$ is the non-local correlation correction (see below). In order to obtain the first three contributions for the different interlayer distances, we perform for each of them self consistent calculations for LDA, PBE, and revPBE.

To perform these calculations,you need to execute the script EXECUTE-graphite.sh in your work directory:

$ EXECUTE-graphite.sh

This script performs a series of runs for a set of graphite interlayer distances (specified inside the script). It also generates a real-space 3D charge density in a format suitable for visualization with XCrySDen. All important files will be stored a directory called output.

Calculation of the non-local correlation energy correction

The charge-density distribution is the only input for the program noloco (located in $EXCITINGROOT/src/src_vdwdf), which evaluate the last term, $E_{\rm c}^{\rm NL}$, of Eq. (1). The non-local correlation correction $E_{\rm c}^{\rm NL}$ is given by the six-dimensional integral:

(2)
\begin{align} E_{\rm c}^{\rm NL} = \frac{1}{2} \int \int n(\textbf{r})\, \phi(\textbf{r},\textbf{r}')\, n(\textbf{r}')\, d\textbf{r}\, d\textbf{r}', \end{align}

where the integration is over the crystal volume. noloco performs the integration by using the Monte-Carlo (MC) method implemented in the external CUBA library. In this tutorial, you do not need to run explicitly the noloco program. Indeed, this is done automatically by using the script $EXCITINGSCRIPTS/EXECUTE-noloco.sh. If you give a look to this script, you can identify the input file for noloco in the file input.in, which looks like the following

# Integration parameters
  $vdwdf_path            ! path to vdW-DF kernel 
  vdW-DF                 ! vdW-DF type (vdW-DF, vdW-DF2, VV09)
  4 4 4                  ! number of unit cells in x/y/z-direction
  1.d-16                 ! relative integration accuracy
  1.d-04                 ! absolute integration accuracy in Hartree
  10000                  ! number of Monte-Carlo sampling points
  1000000000             ! maximum number of function evaluations

# Density data file      
  bohr                   ! bohr/angs: charge density units [e/(bohr^3) or e/(angs^3)]
  ${file}                ! density

Some of the parameters defined inside this input file are system-dependent and must be changed when performing convergence tests, e.g.:

  • the "number of unit cells". These are three integer numbers which determine the size (in both the + and - directions along the axis $a$, $b$, and $c$, respectively) of the supercell considered as integration volume in Eq. (2). These numbers can also be interpreted as "cutoff radii" for the vdW interaction;
  • the "relative accuracy" and "absolute accuracy" of the MC integration.

In order to modify the parameters, the script EXECUTE-noloco.sh should be modified, too. Therefore, is a good practice, before proceeding further, to make a copy (which can be later modified) of the script $EXCITINGSCRIPTS/EXECUTE-noloco.sh inside the actual directory.

$ cp $EXCITINGSCRIPTS/EXECUTE-noloco.sh ./EXECUTE-noloco.sh

To obtain $E_{\rm c}^{\rm NL}$, one just needs to execute the script ./EXECUTE-noloco.sh after the previous calculations have been finished:

$ ./EXECUTE-noloco.sh
Visualization of results

As a final step, we visualize the result of the calculation. The script PLOT-graphite.sh processes all output files.

$ PLOT-graphite.sh > plot-data

Now, you can visualize the result with xmgrace.

$ xmgrace -nxy plot-data &

To obtain nice plots, you can copy the parameter file $EXCITINGSCRIPTS/graphite_xmgr.par to your actual directory

$ cp $EXCITINGSCRIPTS/graphite_xmgr.par ./

and load it in xmgrace using the options Plot -> Load parameters. The resulting plot should look like the following one.


EB.png

Notice that, due to the stochastic character of the MC integration, results obtained performing different runs with the same initial conditions may not be the same. Indeed, this effect is related to statistical errors that can be minimized only for "converged" MC calculation (see below the Exercises section).

Additional exercises
  • Perform convergence tests to find the optimal values for rgkmax and ngridk.
  • Improve vdW-DF results by increasing the number of grid points where the charge density is calculated (specified by the attribute grid of the element box).
  • Check the convergence of the non-local energy with respect to the integration volume, the number of MC sampling points, and the absolute integration accuracy. Change these parameters by editing directly the script ./EXECUTE-noloco.sh. Note, that to obtain small absolute errors, you might have to increase the maximum number of function evaluations.
  • Compare your results to theoretical and experimetal values available in the literature (see, e.g., Adsorption of benzene and naphthalene on graphite, PRL 2006) .
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License