DFT-1/2

by Ronaldo Rodrigues Pela for exciting carbon

Purpose: In this tutorial, you will learn how to run DFT-1/2 calculations using exciting. An explicit example is given for Si and GaN.


0. Define relevant environment variables

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, you can find a list of the scripts which are relevant for this tutorial, with a short description.

  • SETUP-dft0.5.py: Python script for generating a set of input files varying the attribute cut.
  • EXECUTE-single.sh: (Bash) shell script for running a single exciting calculation.
  • EXECUTE-elastic-strain.sh: (Bash) shell script for running a series of exciting calculations.
  • PLOT-plot.py: Python visualization tool which generates plots for files with two columns.

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 are in atomic units!


1. Introduction

The DFT-1/2 method, in principle, is suitable only for correcting the Kohn-Sham KS eigenvalues around the bottom of conduction band and the top of valence band. This means that the method is appropriate for calculating single-particle band gaps, but not for calculating properties that depend, e.g., on the total energy. For compounds with not too large band gaps (typically not larger than ~ 10 eV), it can reach good accuracy with the same computational cost of the standart (semi)local functionals. Within this approach, the KS potential VKS(r) of a standard DFT calculation is replaced by a modified KS potential Vmod-KS(r)=VKS(r)- VS(r). The "correction" VS(r), named as self-energy potential (because its similarity to a classical eletrostatic self-energy), is calculated for the atoms which form the crystal, and is transferred to the crystal after a trimming process defined by the function:

(1)
\begin{align} \Theta (r) = \left\{ \begin{array}{ll} A\left[1-\left(\frac{r}{r_{\rm cut}} \right)^{\!\!n} \right]^3 & r \le r_{\rm cut} \\ 0 & r>r_{\rm cut} \end{array} \right. \end{align}

the $r_{\rm cut}$ is to be determined variationally (see Section 3), $n$ is usually taken as 8, and $A$, the amplitude of the correction, is equal to 1. Note that setting $r_{\rm cut}=0$ corresponds to a standard DFT calculation. Depending on the exchange and correlation employed, one can have, e.g., LDA-1/2 or GGA-1/2. Here, we will employ LDA-1/2. The suffix "-1/2" comes from the fact that the corrective potential mimics the half-occupation scheme of Slater. The self-energy potential is defined for each atom that comprises the crystal as a difference of the KS potential for the neutral atom and the KS potential for the half-ionized one (here, the shell to be ionized is that one that gives rise to the valence band).


2. Running a simple LDA-1/2 calculation

Create a working directory tutorial-lda05 and move inside it.

$ mkdir tutorial-lda05
$ cd tutorial-lda05

Here is one input file (input.xml) for bulk Si.

<input>
 
   <title>Bulk Silicon: LDA-1/2 example</title>
 
   <structure speciespath="$EXCITINGROOT/species">
      <crystal scale="10.26">
         <basevect>0.0   0.5   0.5</basevect>
         <basevect>0.5   0.0   0.5</basevect>
         <basevect>0.5   0.5   0.0</basevect>
      </crystal>
      <species speciesfile="Si.xml" rmt="2.1">
         <atom coord="0.00  0.00  0.00"></atom>
         <atom coord="0.25  0.25  0.25"></atom>
         <dfthalfparam 
            cut="3.90" 
            ampl="1" 
            exponent="8">
            <shell number="0" ionization="0.25" />
         </dfthalfparam>
      </species>
   </structure>
 
   <groundstate
      do="fromscratch"
      rgkmax="7.0"
      gmaxvr="14"
      ngridk="6 6 6"
      outputlevel="high"
      xctype="LDA_PW">
      <dfthalf printVSfile="false"/>
   </groundstate>
 
</input>

Let us take a closer look at this input file. Please, note that in production calculations one should always check the appropriate value of the attribute cut as well as the values of the attributes of the element shell which identify the shell that should be ionized, namely number and ionization). The element dfthalf triggers the DFT-1/2 calculation: Without it, a standard DFT calculation is performed. The description of the attributes and elements concerning the DFT-1/2 are given in the table below:

Parameter Description
dfthalfparam/cut Specifies the value of $r_{\rm cut}$ (in Bohrs) for the self-energy potential, which determines the range of the correction. Note that this can/should be specified for each species.
dfthalfparam/ampl Specifies the amplitude of the self-energy potential.
dfthalfparam/exponent The exponent $n$ that appears in the function that trims the correction.
shell/number, shell/ionization The value of these attributes determines the number of the shell that is to be ionized and the ionizaction that must be applied to this shell. The attribute number follows the same convention used in the species file, where 0 is the shortcut that defines the last atomic shell. According to the convention, the 1s shell is always identified by number="1" (in our example for Si, a number equal to 7 would have the same result as 0). The attribute ionization specifies the amount of charge that is to be removed from this shell when generating the corresponding VS potential.
dfthalf/printVSfile Optional attribute. When set to "true", the self-energy correction potential VS (as defined in the DFT-1/2 method) is calculated for each constituent atomic species and written into the files VS_S*.OUT, where * ranges from 1 to the number of atomic species. The exciting run quits after the printing. In this case, a serial calculation is suggested. It is useful to visualize the self-energy potential, or for debugging purposes.

You can see that for Si, because there are two identical atoms in the unit cell, we need to remove 1/4 electron of the last shell in order to have the ionization scheme of 1/2 in the valence band. The shell that is being ionized is the one that corresponds to 3p orbital, which will further form the top of valence band in the bulk. When performing DFT-1/2 calculations, we always suggest to use

   <groundstate
      ...
      outputlevel="high"
      ...>
     ...
   </groundstate>

within the groundstate element: This will print in the output INFO.OUT the contribution of the VS potential to the total energy (basically, the integral of this potential times the electronic density evaluated in the entire unit cell).

To perform the actual calculation, copy and paste the text from above into a file called input.xml within a directory of your choice. Make sure to set $EXCITINGROOT to the correct exciting root directory in the speciespath attribute using the command

$ SETUP-excitingroot.sh

Now, you can start the calculation by invoking the exciting executable (serial version)

$ EXECUTE-single.sh lda05-test-01

The output files of this run are stored inside the directory lda05-test-01. In order to analyze the results, move to the directory lda05-test-01 when the execution has been completed.

$ cd lda05-test-01

After running this calculation, one can inspect the INFO.OUT and look for the band gap. In this case, it is obtained a gap of about 0.0435 Ha, or 1.18 eV, which is quite close to the experimental band gap of Si (1.17 eV).


3. Optimizing the choice of the parameter "cut"

3.1 Bulk Si

The parameter $r_{\rm cut}$, which defines the range of the corrective potential within the DFT-1/2 method, is specified by the value of the attribute cut. The "right" value of $r_{\rm cut}$, to be used in production calculations, must be determined variationally, by making the band gap extreme. This is achieved in the example of Si using the following procedure. First, $r_{\rm cut}$ is varied between 0 and 5 Bohr (a step of 0.5 Bohr is a good choice, but should be refined). Then, run the calculations and search for the one that gives a maximum to the band gap of Si. You can do this manually, but it would be a tedious task. We will provide a script that takes care of this job. Now, go back to the tutorial-lda05 directory and create a second folder (cut-Si) to perform this kind of calculation. Then enter the new directory.

$ cd ../
$ mkdir cut-Si
$ cd cut-Si/

Now, create an input.xml file as the one given in Section 2. Make sure to set $EXCITINGROOT to the correct exciting root directory in the speciespath attribute using the command

$ SETUP-excitingroot.sh

To execute the series of calculations varying the $r_{\rm cut}$, you need first to run the script SETUP-dft0.5.py. This script always generates a working directory containing input files for different values of $r_{\rm cut}$. The results of the current calculations will be also stored in this working directory, whose name can be specified by adding an argument when running SETUP-dft0.5.py:

$ SETUP-dft0.5.py DIRECTORYNAME

When no name is given, the script uses as default the directory workdir. Very important: The workdir directory is overwritten each time you execute the script SETUP-dft05.py. Therefore, it is a wise choice to change the name when performing a different calculation. Choosing coarse as working directory, you should run SETUP-dft0.5.py as follows:

$ SETUP-dft0.5.py coarse

Enter minimum value of r_cut_min [Bohr]            >>>> 0

Enter maximum value of r_cut_max [Bohr]            >>>> 5.0

Enter the number of steps in [r_cut_min,r_cut_max] >>>> 11

In this example, (on screen) input entries are preceded by the symbol ">>>>". You must type the entry values when requested. The first entry (here 0 in the example) determines the minimum value of $r_{\rm cut}$ in the series of calculations. The second entry (5.0) is, on the other hand, the maximum value of $r_{\rm cut}$. The third (last) entry (11) determines the number of steps to be considered between the minimum and maximum values of $r_{\rm cut}$ (here the number 11 provides a step of 0.5, as we recommended before).

After running SETUP-dft0.5.py coarse, a directory called coarse is created, containing input files for different values of $r_{\rm cut}$. To execute the series of calculation with input files created by SETUP-dft0.5.py coarse, 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. As a final remark, it is worth mentioning that, despite its name, the script EXECUTE-elastic-strain.sh, in DFT-1/2 calculations, is not supposed to execute strain calculations, as done in Energy vs. Strain Calculations - here the same script is reused with a different purpose.

$ EXECUTE-elastic-strain.sh coarse

===> Output directory is "coarse/" <===

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

...

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

This set of calculations should last some minutes (typically 1 or 2). After the complete run, move to the working directory coarse.

$ cd coarse

Inside this directory, results of the calculation for the input file input-i.xml are contained in the subdirectory rundir-i where i is running from 01 to the total number of steps (in our case, 11). During the calculation, results for the value of the band gap (in eV) as a function of the $r_{\rm cut}$ parameter (in Bohr) are stored in the file bandgap-vs-cut. The content of this file should look like this

 0.00000000           0.6059471787
 0.50000000           0.6190981743
 1.00000000           0.6225785119
 1.50000000           0.6392533833
 2.00000000           0.7316289982
 2.50000000           0.9076366622
 3.00000000           1.1026218260
 3.50000000           1.1982429441
 4.00000000           1.1690799467
 4.50000000           1.0394505699
 5.00000000           0.8522402462
The content of the bandgap-vs-cut file can be also visualized opening with the appropriate tool the files PLOT.png and PLOT.ps. These files are automatically generated by the script EXECUTE-elastic-strain.sh during the run. In the example presented here, you should get something similar to the following plot.
bandgap-Si-coarse.png

Now, you can refine the search for the $r_{\rm cut}$ that maximizes the band gap. Go back to the previous directory.

$ cd ../

and setup a new refined calculation (using the working directory fine).

$ SETUP-dft0.5.py fine

Enter minimum value of r_cut_min [Bohr]            >>>> 3.6

Enter maximum value of r_cut_max [Bohr]            >>>> 3.9

Enter the number of steps in [r_cut_min,r_cut_max] >>>> 4

After this, you can run this new set of calculations by entering

$ EXECUTE-elastic-strain.sh fine

Inside the directory fine, the file bandgap-vs-cut will contain these data

 3.60000000           1.2020190696
 3.70000000           1.2007953731
 3.80000000           1.1947438307
 3.90000000           1.1840883922

In order to comprise all the results that you have so far, you can concatenate the file bandgap-vs-cut inside the directory fine, with that one inside coarse. Supposing that your current directory is cut-Si, you can easily do this using the following command.

$ cat coarse/bandgap-vs-rcut fine/bandgap-vs-rcut | sort > ./bandgap-with-all

The file bandgap-with-all, located in the current directory, will have all the data from the previous two series of calculations, sorted by the value of $r_{\rm cut}$. Its content is:

 0.00000000           0.6059471787
 0.50000000           0.6190981743
 1.00000000           0.6225785119
 1.50000000           0.6392533833
 2.00000000           0.7316289982
 2.50000000           0.9076366622
 3.00000000           1.1026218260
 3.50000000           1.1982429441
 3.60000000           1.2020190696
 3.70000000           1.2007953731
 3.80000000           1.1947438307
 3.90000000           1.1840883922
 4.00000000           1.1690799467
 4.50000000           1.0394505699
 5.00000000           0.8522402462

If you want to visualize this new file bandgap-with-all, you can use the script PLOT-plot.py:

$ PLOT-plot.py bandgap-with-all "o-" "rcut [Bohr]" "Bandgap [eV]"
This generates in the current folder two files PLOT.png and PLOT.ps with the graph of band gap vs. $r_{ \rm cut}$. The corresponding plot should appear like the following.
bandgap-Si-with-all.png

Now, you see that the $r_{ \rm cut}$ that maximizes the band gap lies between 3.6 and 3.7 Bohr (if you want, you can refine it further), and the band gap obtained with LDA-1/2 is 1.20 eV. Note that, for a more precise result, you should have changed since the beginning to a better ngridk (for Si, this does matter, because it has an indirect band gap).

3.2. Bulk GaN

Once you have already applied the LDA-1/2 method to bulk Si, you can try now to apply the method to GaN. The ground-state structure of this material is the hexagonal wurtzite. However, we will employ here the zinc-blende phase, because it is easier and faster to be calculated. Go back to the previous directory where you were saving the calculations done in this tutorial, create a new folder and enter inside it.

$ cd ../
$ mkdir cut-GaN
$ cd cut-GaN/

The procedure to apply the LDA-1/2 to GaN is similar to that one you learned for Si, but now, in principle, you should care about "ionizing" Ga and N. As the orbitals of N are the ones that most contribute to the top of valence band of the crystal, a very good approximation is to consider only the self-energy potential for N. This is what we will do here.

To determine the $r_{ \rm cut}$ of N that maximizes the band gap, you can start by inserting the following content

<input>
 
   <title>Zinc-blende GaN: LDA-1/2 example</title>
 
   <structure speciespath="$EXCITINGROOT/species">
      <crystal scale="8.5038">
         <basevect>0.0   0.5   0.5</basevect>
         <basevect>0.5   0.0   0.5</basevect>
         <basevect>0.5   0.5   0.0</basevect>
      </crystal>
      <species speciesfile="Ga.xml" rmt="1.6">
         <atom coord="0.00 0.00 0.00"/>
         <dfthalfparam 
            cut="0.00" 
            ampl="1" 
            exponent="8">
            <shell number="-2" ionization="0.5" />
         </dfthalfparam>
      </species>
      <species speciesfile="N.xml" rmt="1.6">
         <atom coord="0.25 0.25 0.25"/>
         <dfthalfparam 
            cut="2.0" 
            ampl="1" 
            exponent="8">
            <shell number="0" ionization="0.5" />
         </dfthalfparam>
      </species>
   </structure>
 
   <groundstate
      do="fromscratch"
      rgkmax="7.0"
      gmaxvr="14"
      ngridk="6 6 6"
      outputlevel="high"
      xctype="LDA_PW">
      <dfthalf printVSfile="false"/>
   </groundstate>
 
</input>

into a file input.xml inside the current directory. Then, proceed analogously as for Si. As a remark, it is worth mentioning that the script SETUP-dft0.5.py will always change the value of $r_{ \rm cut}$ of the last species that appears inside the input.xml file.

We provide some data (of $r_{ \rm cut}$ and band gap) for you to check your results.

 0.00000000           1.7579272119
 0.50000000           1.7457199075
 1.00000000           1.8212402394
 1.50000000           1.8342524569
 2.00000000           2.4185175951
 2.50000000           3.1272909254
 2.60000000           3.2256144072
 2.70000000           3.2989148334
 2.80000000           3.3474207797
 2.90000000           3.3727222080
 3.00000000           3.3756896107
 3.10000000           3.3563031237
 3.20000000           3.3163902642
 3.30000000           3.2592117738
 3.40000000           3.1877309736
 3.50000000           3.1036047656
 4.00000000           2.5624661528
 4.50000000           2.0000485003
 5.00000000           1.6091127124

Using as above the script PLOT-plot.py and visualizing the obtained PLOT.png or PLOT.ps files, you should get something like the following.

bandgap-GaN-with-all.png

This leads you to the conclusion that the $r_{ \rm cut}$ of N (that maximizes the band gap) is 3.0 (within an accuracy of 0.1) and, that the band gap of GaN calculated through the LDA-1/2 method is 3.38 eV (which exhibits nice agreement with the experimental one, of 3.3 eV).

3.3. Band structure and DOS

Of course it is possible to plot band-structures and density of states when using DFT-1/2. You just need to follow the same procedure as in Electronic Structure Calculations, adding to input.xml the elements and attributes that are specific for DFT-1/2 calculations. You can plot these graphs for Si and GaN as an exercise.

3.4. Transferability

The value of the parameter $r_{ \rm cut}$ that maximizes the band gap should be obtained for each material to be studied with the LDA-1/2 method. However, in many cases, you can assume that it is transferable (indeed, this has been verified in many calculations). Therefore, if you are interested in calculations of interfaces, you can determine the $r_{ \rm cut}^{\rm optimized}$ separately for each material that comprises the interface, and then employ these "optimal" $r_{ \rm cut}$ for the interface. Another example of transferability would appear considering our example of GaN and supposing that you want to study, for example, AlN: in principle, the same $r_{ \rm cut}^{\rm optimized}$ obtained for N in GaN can be applied to AlN.


Literature

  • More details about the DFT-1/2 formalism: L.G. Ferreira, M. Marques, L.K. Teles. "Approximation to density functional theory for the calculation of band gaps of semiconductors", Phys. Rev. B 78, 125116 (2008); and L.G. Ferreira, M. Marques, L.K. Teles. "Slater half-occupation technique revisited: the LDA-1/2 and GGA-1/2 approaches for atomic ionization energies and band gaps in semiconductors", AIP Adv. 1, 032119 (2011).
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License