Vaspkit对格点文件进行平面平均

Jincheng Liu VASPKIT 格点文件静电势功函

Created At : 2019-05-21 00:00

Views 👀 :

Example, Work-function of Au(111) slab with five atomic layers.

本文章为原创，版权归作者刘锦程所有，文章转载请先取得作者的同意，非常欢迎转发文章链接！严禁以任何方式挪用本文内容，用于以盈利为目的各种活动。

VASP计算结果的电荷密度(CHGCAR，PARCHG)，静电势（LOCPOT），ELF（ELFCAR）都是格点文件，可以用VESTA直接打开。

VASP输出的所有格点文件Fortran语言的书写顺序如下：

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

用VASPKIT做平面平均可以得到比较直观的结果。原理就是选定一个方向（比如z方向），在其他两个方向（比如x和y）做平均。

For one selected direction, such as z direction. VASPKIT do an average in XY plane on each NZ,
$$
\sum_{i j} \Delta x_{i} \Delta y_{j} \rho_{i, j} / \sum_{i j} \Delta x_{i} \Delta y_{j}
$$
unfold band

The output POTPAVG.dat contains two columns. First column is z coordinate (in Å), 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:

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

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

Example, 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:

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

(4) Draw the Plane-averaged electrostatic potential figure.

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

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