Dipole Correction for Surface Calculations

by Qiang Fu & Dmitrii Nabok for exciting oxygen

Purpose: In this tutorial, you will learn how to set up and execute exciting calculations for dipole correction on the electrostatic potential of two-dimensional slabs containing a dipole moment in the vertical direction. An explicit example is the half-hydrogenated and half-fluorinated graphene monolayer (HCCF), in which the hydrogen atoms bond with carbon atoms on one side and fluorine atoms on the other side. It is also shown how to calculate and visualize the electrostatic potential, with and without the above dipole correction.


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

  • EXECUTE-planarAverage.py: Python script for extracting planar-averaged electrostatic potential in a given direction.
  • PLOT-files.py: Python visualization tool (see The python script "PLOT-files.py").

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

DFT calculations of surface systems can be performed in exciting by employing the "supercell approach", that is, embedding a slab of material in a periodic supercell. This approach is similar to the procedure used for calculating properties of molecules, as described in How to run calculations for simple molecules.

Using the supercell approach to simulate asymmetric slab containing a net surface dipole density will cause an artificial electric field across the slab. The dipole correction scheme compensates for the artificial dipole field within the context of periodic supercell calculations by introducing an additional ramp-shaped potential that cancels that field. Assuming that the slab is set to be in the xy plane, the compensating potential is given by

(1)
\begin{align} V_{\rm dip}(z) = 4 \pi m \left( \frac{z}{z_m}-\frac{1}{2} \right), \end{align}

where $m$ is a surface dipole moment and $z_m$ is the thickness of the box in the z direction. The correction introduces a jump in electrostatic potential which should, of course, be placed within the vacuum region of the supercell.


2. Running a default calculation without applying dipole correction

Create a working directory dipole-correction and move inside it.

$ mkdir dipole-correction
$ cd dipole-correction

For the calculation to be performed without applying the dipole correction, create a working directory DFT and move inside it.

$ mkdir DFT
$ cd DFT
i) Ground-state calculation for electrostatic potential

Here is the input file (input.xml) for HCCF without applying the dipole correction. In other words, it is a standard exciting calculation.

<input>
 
   <title>HCCF</title>
 
   <structure speciespath="$EXCITINGROOT/species" 
              cartesian="true" tshift="false">
      <crystal>
         <basevect>  4.92460000    0.00000000    0.00000000</basevect>
         <basevect>  2.46230000    4.26482870    0.00000000</basevect>
         <basevect>  0.00000000    0.00000000   28.00000000</basevect>
      </crystal>
      <species speciesfile="C.xml" rmt="1.25">
         <atom coord="  0.000000    0.000000   10.549058"/>
         <atom coord="  4.924600    2.843219   11.426802"/>
      </species>
      <species speciesfile="F.xml" rmt="1.20">
         <atom coord="  4.924600    2.843219   13.999961"/>
      </species>
      <species speciesfile="H.xml" rmt="0.80">
         <atom coord="  0.000000    0.000000    8.452276"/>
      </species>
   </structure>
 
   <groundstate 
      do="fromscratch"
      xctype="LDA_PW"
      ngridk="6 6 1"
      gmaxvr="14"
      rgkmax="3.5"
      outputlevel="normal">
   </groundstate>
 
   <properties>
      <exccplot>
         <plot3d>
            <box grid="24 24 125">
               <origin coord=" 0.0  0.0  0.0"/>
               <point  coord=" 1.0  0.0  0.0"/>
               <point  coord=" 0.0  1.0  0.0"/>
               <point  coord=" 0.0  0.0  1.0"/>
            </box>
         </plot3d>
      </exccplot>
   </properties>
 
</input>

Make sure to set $EXCITINGROOT to the correct exciting root directory in the speciespath attribute using the command

$ SETUP-excitingroot.sh

If the visualization program XCrySDen is set up appropriately (find here how to do this: XCrySDen Setup for exciting), you can visualize the structure of HCCF in the exciting input file by executing

$ xcrysden --exciting input.xml >/dev/null 2>&1 &
You will find the structure of the HCCF monolayer, like in the figure below. Here, all the fluorine atoms are located on one side, and the hydrogen atoms on the other side. Because of the electronegativity difference between F, H, and C atoms, there will be charge accumulation at the F side, and charge depletion at the H side. In other words, the monolayer has a dipole moment.
HCCF.jpg

A detailed description on most of the attributes used here inside the element properties can be found at How to Visualize Kohn-Sham States. Here, the attribute exccplot is a new one, from which the electrostatic potential and the exchange-correlation potential can be plotted in real space.

Then, start the ground-state calculation by executing the following command:

$ time exciting_smp

The electrostatic potential and the exchange-correlation potential are stored inside the file VCL3D.xml and VXC3D.xml, respectively.

ii) Visualization of the electrostatic potential

The file VCL3D.xml contains all information necessary for the visualization of the electrostatic potential. For the two-dimensional slab, it is always a good idea to visualize the plane-averaged electrostatic potential along the vertical (z) direction.

$ EXECUTE-planarAverage.py VCL3D.xml z

Then you will get a new file named VCL3D_Z.dat. In order to visualize the results, just move back to the parent directory dipole-correction and execute there the script PLOT-files.py (for a detailed description of the script arguments see The python script "PLOT-files.py") as follows

$ cd ../
$ PLOT-files.py -f planarAverage_z  -d DFT  -x 0 28  -y -0.4 0.3  -lx 'z [bohr]'  -ly '$V_{elec}$ [Ha]'  -t ""  -lp 4  -ll DFT

The resulting plot of the electrostatic potential Velec(z) is contained in the file PLOT.png and will look like this:

no-dipole-correction.png

From this figure, you can see that the electrostatic potential within the vacuum layer is not flat, and the vacuum energy level cannot be identified. In principle, the true electrostatic potential of a two-dimensional slab, which contains a dipole moment in the vertical direction, will be flat within the vacuum layer but has a jump. However, owing to the periodic boundary condition applied in the "supercell approach", the electrostatic potential should be continuous at the boundary of the supercell, which makes the electrostatic potential at the right side decreases, and that at the left side increases, in order to maintain the continuity. By applying the dipole correction, the electrostatic potential will be flat, and the potential jump will be recovered, as presented in the next section.


3. Running a calculation with dipole correction

Inside the directory dipole-correction, create a new one named DFT-dc, and then move inside it.

$ mkdir DFT-dc
$ cd DFT-dc

Then, copy the files input.xml and STATE.OUT, from where your last calculation was performed, to the current directory.

$ cp ../DFT/input.xml ./
$ cp ../DFT/STATE.OUT ./

The input file input.xml must be changed in order to perform the new calculation for HCCF including the dipole correction. In this case, three changes have to be made inside the groundstate element, as shown in the following.

...
   <groundstate 
      do="fromfile"
      dipolecorrection="true"
      dipoleposition="0.9"
      ...>
   </groundstate>
...
  1. Change the attribute do="fromscratch" to do="fromfile". It means that you will use the pre-converged STATE.OUT file from the previous calculation. In some cases, this will improve the performance of self-consistent convergence after the corrected potential is added on the electrostatic potential.
  2. Add the attribute dipolecorrection="true". Then, the functionality "dipole correction" is switched on.
  3. Add the attribute dipoleposition="0.9". The value of the attribute dipoleposition indicates the position of the jump in electrostatic potential, after the compensating potential (i.e., the dipole correction) is applied. The value of "0.9" is a fractional coordinate in the vertical direction. Please note that this jump position should be located within the vacuum region enough far away from the atomic layers, otherwise the compensating potential cannot be correctly applied. It is recommended to put the jump position at the middle of the vacuum layer.

Then, start a new ground-state calculation by executing the command:

$ time exciting_smp
$ EXECUTE-planarAverage.py VCL3D.xml z

To visualize the effect of the applied dipole correction, one can follow the same procedure described above

$ cd ../
$ PLOT-files.py -f planarAverage_z  -d DFT DFT-dc  -x 0 28  -y -0.4 0.3  -lx 'z [bohr]'  -ly '$V_{elec}$ [Ha]'  -t ""  -lp 4

The resulting electrostatic potentials with and without the dipole correction will look now like this:

dipole-correction.png

In this figure, you can see how the electrostatic potential has been corrected (blue line). It is flat on both sides and has a jump. The plateau in the vacuum layer refers to the vacuum energy level. By comparing its position with those of the corresponding valence band maximum (VBM) and conduction band minimum (CBM), you can calculate the ionization energy (IE) and electron affinity (EA), accordingly. Please note that here, due to the existence of a dipole, two platforms appear in the plane-averaged electrostatic potential. It means that the vacuum energy level is not the same at the two sides of HCCF, and, consequently, the IE and EA either.


Literature

  • Lennart Bengtsson, "Dipole correction for surface supercell calculations", Phys. Rev. B 59, 12301 (1999).
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License