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

  1. 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:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
##### 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:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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:

1
2
3
4
5
6
7
8
9
10
# 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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
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:

unfold band unfold band

So, work function can be calculated as:
$$
\Phi=E_{\mathrm{vac}}-E_{\mathrm{F}}
$$
$\Phi$ = 6.695 - 1.5199 = 5.1751 eV


欢迎在评论区讨论问题、补充内容、指出错误。请勿发送占楼、沙发、点赞等无意义的回复。
目录