水文一篇,介绍了如何用 VASP 添加偶极校正参数、计算 Slab 体系的功函数,以及如何找真空能级。

概念解释

功函数

将一个固体内部的电子移动到真空所需的能量。(类似于光电效应方程中的逸出功)。

真空能级

固体表面外真空中自由电子所具有的能量。换句话说,电子跑出固体表面并达到这个能级后即可认为它自由 免费 了。

偶极校正

因 VASP 所适用的体系是周期性体系,使用它来模拟实验中的 Slab 模型时会取一个相当大的真空层来隔绝相信两个周期中 Slab 的相互作用。理想情况下,真空层中的功函数应当是一条水平的直线(函数值为定值)。但如果表面的两侧并非对称,即其中一侧吸附了分子时, 这两侧的功函数存在差异,此时如果不进行偶极校正,真空中的功函数会是一条斜线;而经过偶极校正后,功函数会出现一个阶梯,阶梯两侧附近的曲接近水平。下图1是一个例子。

Figure 1: DFT 曲线为未经过偶极校正的功函数, DFT-DC 曲线是经过偶极校正后的功函数

Figure 1: DFT 曲线为未经过偶极校正的功函数, DFT-DC 曲线是经过偶极校正后的功函数

如何计算

偶极校正

VASP 中直接使用 LDIPOLIDPOL 即可开启它的偶极校正功能。

  • LDIPOL = .TRUE. 表示打开偶极校正;
  • IDIPOL = 3 表示偶极校正所修正的方向为第 3 个晶格矢量方向对应的方向,一般来说就是 \(z\) 轴;
  • DIPOL = <3 float values> 表示体系的中心,以分数坐标表示;

上面几个 Flag 中一般来说 DIPOL 不用填,因为 VASP 手册中2有一句

If the flag is not set, VASP determines, where the charge density averaged over one plane drops to a minimum and calculates the center of the charge distribution by adding half of the lattice vector perpendicular to the plane where the charge density has a minimum (this is a rather reliable approach for orthorhombic cells).

不过有时不填它会导致体系 非常 难以收敛,如果遇到这各情况,最好还是手动算一下所有原子 \(z\) 坐标的平均值,然后 DIPOL = 0.5 0.5 <averaged z> ,此时修正效果可能并不好,读者需要:

  1. 打开偶极校正和 DIPOL = 0.5 0.5 <averaged z> ,并弛豫结构到稳定;
  2. 关闭 DIPOL ,并弛豫到结构稳定;
  3. 关闭 DIPOL ,打开 LVHAR = .TRUE. ,并做静态计算,得到功函数;

偶极校正的标准是功函数出现台阶,并且台阶两边为水平,如果画出的功函数图不是这样, 就要考虑调整参数了。

功函数

在 VASP 中 LVHAR 3参数可以使 VASP 输出体系的功函数文件 LOCPOT 。LOCPOT 文件本身是 Volumetric data ,它的格式与 CHGCAR 一样4。一般而言,用户关心的功函数是垂直于表面方向上的数据,因此在得到 LOCPOT 后需要对它做一点工作,即对 \(xy\) 平面内的数据做平均,然后乘以晶胞的体积,就得到我们需要的功函数信息。这里给出一个脚本5来完成这项工作:

 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
33
34
35
36
37
38
39
40
41
42
43
44
45
def locpot_mean(fname="LOCPOT", axis='z', savefile='locpot.dat', outcar="OUTCAR"):
    '''
    Reads the LOCPOT file and calculate the average potential along `axis`.
     @in: See function argument.
    @out:
          - xvals: grid data along selected axis;
          - mean: averaged potential corresponding to `xvals`.
    '''
    def get_efermi(outcar="OUTCAR"):
        if not os.path.isfile(outcar):
            logger.warning("OUTCAR file not found. E-fermi set to 0.0eV")
            return None
        txt = open(outcar).read()
        efermi = re.search(r'E-fermi :\s*([-+]?[0-9]+[.]?[0-9]*([eE][-+]?[0-9]+)?)', txt).groups()[0]
        logger.info("Found E-fermi = {}".format(efermi))
        return float(efermi)

    logger.info("Loading LOCPOT file {}".format(fname))
    locd = VaspChargeDensity(fname)
    cell = locd.atoms[0].cell
    latlens = np.linalg.norm(cell, axis=1)
    vol = np.linalg.det(cell)

    iaxis = ['x', 'y', 'z'].index(axis.lower())
    axes = [0, 1, 2]
    axes.remove(iaxis)
    axes = tuple(axes)

    locpot = locd.chg[0]
    # must multiply with cell volume, similar to CHGCAR
    logger.info("Calculating workfunction along {} axis".format(axis))
    mean = np.mean(locpot, axes) * vol

    xvals = np.linspace(0, latlens[iaxis], locpot.shape[iaxis])

    # save to 'locpot.dat'
    efermi = get_efermi(outcar)
    logger.info("Saving raw data to {}".format(savefile))
    if efermi is None:
        np.savetxt(savefile, np.c_[xvals, mean],
                   fmt='%13.5f', header='Distance(A) Potential(eV) # E-fermi not corrected')
    else:
        mean -= efermi
        np.savetxt(savefile, np.c_[xvals, mean],
                   fmt='%13.5f', header='Distance(A) Potential(eV) # E-fermi shifted to 0.0eV')

完整的脚本文件已经放 Gist5 上,当然你也可以直接点击来下载。运行这个脚本后得到的 Workfunction.pdf 和 locpot.dat 就是 \(z\) 方向上的功函数信息。

Figure 2: Workfunction.png 示例

Figure 2: Workfunction.png 示例

真空能级

前面已经提到,真空能级可以读取 locpot.dat 真空部分的数据得到,其实当你打开 LVHAR 时,它也可以通过读取 OUTCAR 得到,比如

1
2
$ grep vacuum OUTCAR
 vacuum level on the upper side and lower side of the slab         2.807         3.188

这里的 upper side vacuum level 是指 Slab 上表面的真空能级(图 2 中 30A 处的平台), lower side vacuum level 自然就是下表面的真空能级了 (图 2 中 35A 处的平台)。

需要注意的是,从 OUTCAR 中 grep 出的真空能级没有经过费米能级修正,它需要减去 OUTCAR 中的 E-fermi 才是实验中测得的功函数的值。所幸的是 plot-workfunc.py 已经做了这个工作,用 locpot.dat 画出来的图就对应实验所测结果。


  1. 图源: http://exciting-code.org/nitrogen-dipole-correction ↩︎

  2. https://www.vasp.at/wiki/index.php/DIPOL ↩︎

  3. 在 VASP 5.2.12 之前的版本中 LVTOT 输出的文件是静电势的贡献;但在 VASP 5.2.12 及之后的版本中,交换关联势的贡献也会被写入 LOCPOT 中, LVHAR 输出的部 分才是真正静电势的贡献,也就是我们想要的部分; ↩︎

  4. CHGCAR 的格式请见 VASP 手册 https://www.vasp.at/wiki/index.php/CHGCAR↩︎

  5. 完整脚本见 https://gist.github.com/Ionizing/1ac92f98e8b00a1cf6f16bd57694ff03↩︎