Van-Der-Waals Energy Functionals

for lithium version of Exciting

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 shell variables are already defined and that the excitingscripts directory has already been downloaded, as specified in Tutorial scripts and environment variables. 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 will appear will be given in atomic units!

1. Compiling noloco

To calculate the non-local contribution to the correlation energy, we use an external program 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 create a directory for each system that you want to investigate. Therefore, one should create a directory, e.g., 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 c axis) and the number of grid points along c direction where the electronic density will be calculated.

<input>
  <title>Graphite</title>
  <structure speciespath="$EXCITINGROOT/species" autormt="true">
    <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">
      <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"
    mixer="msec"
    maxscl="200"
    xctype="GGAPerdew-Burke-Ernzerhof"
    rgkmax="5.0"
    gmaxvr="12"
    ngridk="4 4 2" 
    vkloff="0.5 0.5 0.5"
    >
  </groundstate>
  <properties>
    <chargedensityplot>
      <plot3d>
          <box grid="10 10 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.

Please also note that the actual input file for the exciting calculation (which is always called input.xml) will be created by the script EXECUTE-graphite.sh. Furthermore, for some parameters (rgkmax and ngridk in groundstate as well as grid in box) rather low values are used as it would be required by convergence tests and used here only to speed up the calculations.

Warning! The execution time of the script EXECUTE-graphite.sh can take a while 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)). In order to do this, we perform first a self-consistent run using the PBE-GGA xc functional and generate the charge density, which is then used to evaluate the vdW-DF total energy given by expression:

(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}

where $E^{\rm PBE-GGA}$ and $E_{\rm xc}^{\rm PBE-GGA}$ are correspondingly the total and exchange-correlation energies taken from the previous run. $E_{\rm x}^{\rm revPBE-GGA}$ is the exchange energy calculated using the revPBE-GGA exchange-correlation potential and $E_{\rm c}^{\rm LDA}$ is the correlation energy obtained within LDA.

To perform actual calculations, one needs to execute 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) in order to obtain the first four terms appearing in the right side of Eq. (1). It also generates a real-space 3D charge-density and the density-gradient distribution in a format suitable for visualization with the XCrySDen program. 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, EcNL, of Eq. (1). The non-local correlation correction EcNL 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. To execute the program, the user must provide an input file, as shown below. For the purpose of this tutorial, this file is generated automatically in the script EXECUTE-noloco.sh. An example of the input is shown below.

# Integration parameters
  ${vdwdf_path}          ! path to the tabulated vdW-DF kernel (\phi(r,r'))
  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-03                 ! 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)]
  ${dens}                ! density

Note that the following parameters are system-dependent and must be chosen based on convergence tests:

  1. three integer numbers which determine the integration super cell (obtained by the unit cell translations along (+/-) a, b, and c directions). Another interpretation of these numbers is a cutoff radius of vdW interactions;
  2. the relative and (or) absolute accuracy of the MC integration estimator (the specifics of MC integration).

To modify them in this tutorial, one should edit the script EXECUTE-noloco.sh.

To obtain EcNL, one just needs to execute the following command 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.dat

Now, you can visualize the result with xmgrace.

$ xmgrace -nxy plot.dat

To obtain nice plots one can load the parameters file $EXCITINGSCRIPTS/graphite_xmgr.par in xmgrace -> Plot -> Load parameters… The resulting plot should look like

Ebind-graphite.jpg

From the resulting figure one can clearly see that GGA performs very badly when calculating the equilibrium interlayer distance and binding energy in graphite. In contrast, vdW-DF predicts the both quantities to be close to experimental ones (see for comparison, e.g., Adsorption of benzene and naphthalene on graphite, PRL 2006).

Additional exercises
  • Perform convergence tests to find optimal values for rgkmax, ngridk.
  • Improve the 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) and checking the convergence of the non-local energy with respect to the integration volume and the absolute integration accuracy $\epsilon$.
  • Test the performance of the different exchange-correlation potentials (LDA, GGAs).
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License