DFT calculation (VASP) in surface science

VASP计算过渡态

cNEB计算,Optimizer (IOPT)的选择

cNEB的Optimizer即寻找过渡态的方法,对应与结构优化的搜索方法IBRION,只不过cNEB是同时计算MEP路径上的多个image,由于过渡态计算的受力是投影到MEP路径上的受力,能量不再像结构优化时一样与受力对应,因此过渡态的搜索方法都只能是基于受力的。
为了使用VTST提供的搜索方法,需要先关闭VASP自带的结构优化搜索方法:
IBRION=3; POTIM=0
结构受力大或者计算的受力不是十分准确时,用IOPT=3IOPT=7。 这表明,一开始优化的时候可以用这两个参数,而EDIFF不用特别小。
受力计算十分精确时,用IOPT=1或者IOPT=2,这两种方法收敛较快,optimize globally。 这意味着用IOPT=7; EDIFF=1E-4; EDIFFG=5E-1 进行初步优化后, 可以再用IOPT=1; EDIFF=1E-7; EDIFFG=5E-2 进行进一步优化。

在过渡态计算的时候,IOPT=3和IOPT=7的收敛性时最好的, 受力总的趋势是变小。 但是唯一的担心就是这两种优化的过渡态是global的吗? VTST网站上说We recommend using CG or LBFGS when accurate forces are available. This is essential for evaluating curvatures., 那是不是应该用IOPT=2和IOPT=1得到的曲线比较正确? 原则上讲,应该只要收敛了,搜索方法就不重要的,所有的方法应该都能用,只不过效率会不同。

Simulation of dI/dV spectrum and mapping

dI/dV 对应于LDOS

STM 表征中dI/dV直接与理论计算中局域电子态密度LDOS相对应,LDOS为空间坐标和能量的函数$\rho(r,E)$,dI/dV谱对应特定位置r0的LDOS分布$\rho(r_0,E)$,dI/dV mapping 对应特定能量的DOS的空间分布$\rho(r,E_0)$。
LDOS定义为:

即能量为E的电子态的空间概率分布。回想VASP中模拟STM图的过程正是将$E_F$至$E_F-eV$范围内的态密度进行叠加,我们也可以通过类似的过程进行di/dV模拟。
整个计算过程包括

  • 计算体系的DOSCAR文件,确定DOS中态密度范围
  • 计算DOS中态密度下限至不同能量区间$E_i$的积分LDOS对应得PARCHG文件
  • 对PARCHG文件进行分析,获得特定位置积分LDOS,并微分得LDOS

    计算体系的DOSCAR文件,确定DOS中态密度范围

    DOS 计算得时候参数和静态计算一致,只不过得加上LORBIT=10或11。
    INCAR
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    ISTART = 0  
    ICHARG = 2
    PREC = Accurate
    ENCUT = 480
    EDIFF = 1.0e-06
    ALGO = Fast
    ISMEAR = 0; SIGMA = 0.01
    LREAL = .FALSE.
    LCHARG = .TRUE.
    LAECHG = .TRUE.
    LORBIT = 10

计算得到的DOSCAR包含了DOS信息,其中DOS是通过积分DOS微分而来。

用p4v软件将DOSCAR作图获得DOS图。(直接打开DOSCAR获得的DOS图费米能级未移到横坐标0的位置,而打开vasprun.xml获得的DOS图费米能级位于0)

由图可知,DOS中态密度不为零的能量范围为(-20,6)。(此处费米能级位于0的位置)

计算积分LDOS

复制上述DOS计算后的文件夹,并进入该文件夹,用下面的脚本进行积分LDOS的计算。

Shell script for PARCHG calculation

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
#!/bin/bash
for i in {1..260}
do
Vi=-20
Vf=$(echo "scale=2; -20+26/260*$i"|bc)
echo $i $Vi $Vf
cat >INCAR<<!
ISTART = 1
ICHARG = 0
ENMAX = 480
ISMEAR = 0 ;SIGMA = 0.01
ALGO= Fast
PREC = Single
IVDW = 11
EDIFF = 1E-4

partial charge densities:
LPARD = .TRUE.
LSEPK = .FALSE.
LSEPB = .FALSE.
NBMOD = -3 % Vi and Vf are vs. fermi energy
EINT = $Vi $Vf

mechine:
LPLANE = .TRUE
NCORE = 8
NSIM = 4
!
mpirun -np 32 vasp_std
mv PARCHG $i-PARCHG_$Vf
done

该脚本计算得到-20 到-19.9, 19.8, 19.7…5.8, 5.9, 6.0的260个能量范围的空间态密度积分,对应得文件为1-PARCHG-19.9, 2-PARCHG-19.8 … 259-PARCHG_5.9, 260-PARCHG_6.0。
PARCHG 记录了Vi至Vf区间态密度的积分。

这其实是个空间坐标的函数,在PARCHG中以一个三维矩阵记录态密度信息。

该文件中将Unit cell 分成84x84x140的空间格子,并将各个格子内的态密度与unite cell体积的乘积在PARCHG中记录,记录的顺序为(1,0,0)(2,0,0)…(84,0,0)(0,1,0)…(84,84,0)(0,0,1)…(84,84,140)

计算LDOS谱和LDOS mapping

LDOS谱的计算需要选定空间位置r,此处我们选择距离分子中心高度为$2\overset{\circ}{A}$的位置,对应的分数坐标为:0.5508 0.6530 0.3203,对应的空间网格为46 55 46,对应PARCHG中第(46-1)x84x84+(55-1)x84 +46=322102个数据,PARCHG中密度数据记录是按照一行十列,因此位于第32210行,第2列。按照下面的脚本可以从计算的i-PARCHG-Vf 文件中读取该点的数据,获得该位置积分LDOS曲线。

dIdV_dataAnalyse

1
2
3
4
5
6
7
8
9
10
11
12
13
#!/bin/bash
rm LDOS_data
e=0.01
for i in {1..260}
do
Vi=-20
Vf=$(echo "scale=2; -20+26/260*$i"|bc)
LDOS=$(awk -v line=$(awk '/84 84 140/{print NR}' $i-PARCHG_$Vf) '{if(NR==line+32210){print $2}}' $i-PARCHG_$Vf)
echo $Vf $LDOS >>LDOS_data
echo $i $Vf $LD0S
done
gnuplot
plot LDOS_data

在origin或者matlab中对LDOS_data中的数据进行作图获得积分LDOS曲线,作图的时候数据除去了Unit cell 的体积。

进行微分可获得LDOS曲线

LDOS maping

dI/dV mapping 对应特定能量的LDOS的空间分布$\rho(r,E_0)$可以通过下面的公式求得:

通过下面的脚本可以计算上式积分部分,计算时需要给定$E_0$和极小值$\varepsilon$。

计算$E0±\varepsilon$范围电子态密度空间分布

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
32
#!/bin/bash
E0=-2
for e in 0.005 0.01 0.025 0.05 0.075
do
Vi=$(echo "scale=6; $E0-$e"|bc)
Vf=$(echo "scale=6; $E0+$e"|bc)
echo $Vi $Vf
cat >INCAR<<!
ISTART = 1
ICHARG = 0
ENMAX = 480
ISMEAR = 0 ;SIGMA = 0.01
ALGO= Fast
PREC = Single
IVDW = 11
EDIFF = 1E-4

partial charge densities:
LPARD = .TRUE.
LSEPK = .FALSE.
LSEPB = .FALSE.
NBMOD = -3
EINT = $Vi $Vf

mechine:
LPLANE = .TRUE
NCORE = 8
NSIM = 4
!
echo "EINT = $i "; mpirun -np 48 vasp_std
mv PARCHG PARCHG_$E0_$e
done

由于分母$2\varepsilon$为常数,所以PARCHG中的数据直接对应$\rho(r,E_0)$,因此直接通过PARCHG可以模拟dI/dV mapping。

Guaussian计算电子亲和能

电子亲和能是电子受体和电子给体的一个重要的参数,描述了原子和分子得电子得能力,其定义为的电子后,中性分子的能量和带电荷的分子能量之差。

按照这个定义,亲和能为正值表示倾向于得到电子,而负值倾向于失去电子。 比如典型的电子受体C60的亲和能约为+2.67,具有得电子能力,而苯环的电子亲和能为-1.15,其亲合能受到给电子基团和得电子基团的影响.
那计算方面就很简单了,计算一个中性分子的能量,再计算一个带有一个净电荷的分子的能量,二者相减就是电子亲和能。 中性分子肯定要结构优化,再计算能量,那带有一个净电荷的呢? 多加一个净电荷,体系最稳定的结构肯定会与中性的有所差别。 这时候,分两种情况: 带净电荷的分子就一中性的结构为基础直接进行能量计算,这么算出来的是vertical EA,即垂直电子亲和能,而带净电荷的分子独立地进行结构优化再进行能量计算, 此时算出来的电子亲和能为Adiabatic绝热电子亲和能。

软件或者脚本的安装

vtotav.py

vtotav.py是用来对VASP计算得到的势函数或者电荷密度进行面积分的一个python脚本。该软件的使用需要python环境和ase包。以下是在conda环境下安装该脚本所需要的软件包的命令:
conda install numpy
下载并进入解压目录安装ase软件包,并把bin路径添加到环境变量中
python3 setup.py install --user
由于python3.9不支持time.clock 命令, 将vtotav.py中的time.clock全部替换为time.perf_counter().

文章目录
  1. 1. VASP计算过渡态
    1. 1.1. cNEB计算,Optimizer (IOPT)的选择
  2. 2. Simulation of dI/dV spectrum and mapping
    1. 2.1. dI/dV 对应于LDOS
    2. 2.2. 计算体系的DOSCAR文件,确定DOS中态密度范围
    3. 2.3. 计算积分LDOS
    4. 2.4. 计算LDOS谱和LDOS mapping
    5. 2.5. LDOS maping
  3. 3. Guaussian计算电子亲和能
  4. 4. 软件或者脚本的安装
    1. 4.1. vtotav.py
,