DFT calculation (VASP) in surface science

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定义为：

$$\rho(r,E)=\sum|\psi_\mu(r，E)|^2\delta(E_\mu-E)$$

即能量为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

  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

#!/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区间态密度的积分。

$$\int_{Vi}^{Vf}\rho(r,E)dE =\int_{Vi}^{Vf}|\psi_\mu(r，E_\mu)|^2 dE$$

这其实是个空间坐标的函数，在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

#!/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)$可以通过下面的公式求得：

$$\rho(r,E_0)=\frac{\int_{E_0-\varepsilon}^{E_0+\varepsilon}\rho(r,E)dE}{2\varepsilon}$$

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

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

#!/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。

$$\int_{V_F}^{V}|\psi(r，E)|^2 dE$$