Dipole Correction for Surface Calculations

by Qiang Fu & Dmitrii Nabok for exciting carbon

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 dipole 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 fluoride 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.

• planeAverage.py: Python script for extracting plane-averaged electrostatic potential in given direction.
• tutorial-dc.par: The parameter file to modify the format of the graph plotted by using xmgrace.

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 a 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 DipCor and move inside it.

$ mkdir DipCor
$ cd DipCor

For the calculation to be performed without applying 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 word, it is a default 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 fluoride atoms locate 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.

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 excitingser

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.

$ planeAverage.py VCL3D.xml z

Then you will get a new file named VCL3D_Z.dat. In order to visualize the results, just type

$ xmgrace -nxy VCL3D_Z.dat -param $EXCITINGSCRIPTS/tutorial-dc.par >/dev/null 2>&1 &

The resulting electrostatic potential will look like this:

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 dipole 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

Let us go back to the parent directory DipCor, create a new one named dft-dc, and then move inside it.

$ cd ..
$ 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 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 excitingser

For the visualization of the electrostatic potential, you can use the same procedure described above for the uncorrected case. The resulting electrostatic potential will look now like this:

In this figure, you can see that the electrostatic potential has been corrected. It is flat on both sides and has a jump. The platform 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).