plane-average potential and work function

########## citation ##########

V. Wang, N. Xu, Jin-Cheng Liu VASPKIT: A Pre- and Post-Processing Program for the VASP Code. http://vaspkit.sourceforge.net.

########## citation ##########

Plane-average for grid files is one of the basic function of vaspkit (426). It supports all grid files of vasp output, such as CHGCAR, LOCPOT, ELFCAR. All such grid files write by following commands in Fortran:

WRITE(IU,FORM) (((C(NX,NY,NZ),NX=1,NGXC),NY=1,NGYZ),NZ=1,NGZC)

For one selected direction, such as z direction. Vaspkit do an average in XY plane on each NZ, and output POTPAVG.dat with two columns. First column is z coordination (in Angstrom), Second column is Planar Average-Potential (in eV) or Densitiy (in e).

Here, take LOCPOT as an example. It contains total local potential (in eV) when LVTOT = .TRUE. exists in the INCAR file, or contains electrostatic potential (in eV) when LVHAR = .TRUE. exists in the INCAR file. Electrostatic potential is desirable for the evaluation of the work-function, because the electrostatic potential converges more rapidly to the vacuum level than the total potential. To get the work-function, potential along z direction (i.e. slab normal direction) must be averaged in the slab planes.

Hamiltonian of Kohn-Sham is:

$$hat{H}=-frac{1}{2} nabla^{2}-sum_{A} frac{Z_{A}}{left|boldsymbol{r}-boldsymbol{R}_{A}right|}+int_{infty} frac{rholeft(boldsymbol{r}^{prime}right)}{left|boldsymbol{r}-boldsymbol{r}^{prime}right|} d boldsymbol{r}^{prime}+v_{x c}[rho(boldsymbol{r})]$$

Electrostatic potential is the summation of second (ionic) and third (hartree) parts of Hamiltonian .

For example, calculate the work-function of Au(111) slab with five atomic layers.

(1) Do an optimization and a single-point calculation to get the LOCPOT file, single-point INCAR as following:

##### initial parameters I/O #####
SYSTEM = Au
NCORE = 5 
ISTART = 1 
ICHARG = 1 
LWAVE = .TRUE. 
LCHARG = .TRUE. 
LVTOT = .FALSE. 
LVHAR = .TRUE. 
LELF = .FALSE. 

#### Electronic Relaxation ####
ENCUT = 400 
ISMEAR = 1 
SIGMA = 0.2 
EDIFF = 1E-6 
NELMIN = 5 
NELM = 300 
GGA = PE 
LREAL = Auto 
IDIPOL = 3

(2) Run vaspkit 426 function, and select the surface normal direction:

426
 +-------------------------- Warm Tips --------------------------+
                  You MUST Know What You are Doing  
   Check Convergence of Planar Average Value on Vacuum-Thickness!
 +---------------------------------------------------------------+
 ===================== Specified Direction =======================
 Which Direction is for the Planar Average?           
 1) x Direction                                       
 2) y Direction                                       
 3) z Direction                                       

 0) Quit                                              
 9) Back                                              
 ------------>>
3
 -->> (1) Reading Structural Parameters from LOCPOT File...
 -->> (2) Reading Local Potential From LOCPOT File...
 -->> (3) Written POTPAVG.dat File!
 +-------------------------- Warm Tips --------------------------+
  Check the Convergence of Vacuum-Level Versus Vacuum Thickness!      
 +---------------------------------------------------------------+
 +-------------------------- Summary ----------------------------+
 Vacuum-Level  (eV):    6.695
 +---------------------------------------------------------------+

So, the Vacuum-Level is 6.695, POTPAVG.dat as following:

 # Planar Distance (in Angstrom), Planar Average-Potential (in eV) or -Densitiy (in e)
  0.0000    0.27377302E+01
  0.1025    0.94832051E+00
  0.2049   -0.13842911E+01
  0.3074   -0.40872693E+01
  0.4099   -0.69356914E+01
  0.5123   -0.96882699E+01
  0.6148   -0.12150765E+02
  0.7172   -0.14164741E+02
  ...

(3) get the fermi level by:

grep E-fermi OUTCAR 
E-fermi :   1.5199     XC(G=0):  -6.7056     alpha+bet : -6.4642

(4) Draw figure with origin or matplotlib:

python3 ./planeavg.py 1.5199

import matplotlib as mpl
mpl.use('Agg')

import numpy as np
import matplotlib as mpl
mpl.rcParams['font.size'] = 42.
#mpl.rc('font',family='Times New Roman')
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import sys

fig = plt.figure(figsize=(35,21))    # set figure size
ax = fig.add_subplot(111)            # describing the position of the subplot
plane=np.loadtxt("./POTPAVG.dat")    # read
ax.plot(plane[1:,0],plane[1:,1],lw=4,c='blue')    # remove the first line
xmax, xmin = max(plane[1:,0]), min(plane[1:,0])
ymax, ymin = max(plane[1:,1]), min(plane[1:,1])
plt.xlim(xmin, xmax)
plt.ylim(ymin-1, ymax+1)
plt.xticks(np.arange(0, xmax, 5))
plt.yticks(np.arange(int(ymin), int(ymax), 5))

plt.ylabel("Potential (eV)")
plt.xlabel(r"Distance ($mathrm{AA}$)",labelpad=40)
ax.yaxis.set_minor_locator(ticker.MultipleLocator(1.5))
ax.legend(('Planar average of potential'),
          loc=0,shadow=True,labelspacing=0.1,fontsize=35)
plt.tight_layout()
fig.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9,hspace=0.55)
plt.hlines(float(sys.argv[1]), xmin, xmax, colors = "c", linestyles = "dashed")
plt.savefig('plane_average.png',img_format='eps',dpi=100)

The obtained Plane-averaged electrostatic potential figure:

So, work function can be calculated as:

$$Phi=E_{mathrm{vac}}-E_{mathrm{F}}$$

$Phi$ = 6.695 - 1.5199 = 5.1751 eV

转载请注明来源，欢迎对文章中的引用来源进行考证，欢迎指出任何有错误或不够清晰的表达。可以可以邮件至 [email protected]