【科学代码】能带反折叠VaspBandUnfolding 程序包

蔻享学术代码下载页面：

https://code.koushare.com/#/code/codeDetail?codeId=212

（请复制链接至浏览器打开，或点击左下方“阅读原文”跳转）

VASP（Vienna Ab-inito Simulation Package）是基于贋势平面波基组的第一性原理密度泛函计算程序，是目前材料模拟和计算物质科学研究中十分流行的商用软件之一。在VASP众多的输出文件中，WAVECAR中写入了Kohn-Sham赝波函数的平面波展开系数等信息。为了减少文件大小，该文件采用了二进制方式读写，增加了读取该文件的难度。比如，在VASP中要查看某个轨道，只能通过写PARCHG的方式来进行，且PARCHG中的轨道其实是波函数模平方，而不是真正的波函数，缺少波函数相位的信息。

VaspBandUnfolding程序包利用PYTHON来读取WAVECAR中的相关信息，可以做到以下几个功能：

1. 读取WAVECAR中赝波函数的平面波展开系数，通过快速傅里叶变化（FFT）将生成实空间的贋波函数并写入文件，可以直接用VESTA打开查看。

2. 读取超胞WAVECAR中的信息并结合超胞-单胞转换矩阵，可以做到能带反折叠。

3. 简单介绍一下MoS₂中掺杂Ce原子的能带反折叠计算。

先用如下代码生成KPOINTS

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import numpy as np
from unfold import make_kpath, removeDuplicateKpoints, find_K_from_k, save2VaspKPOINTS

# high-symmetry point of a Hexagonal BZ in fractional coordinate
kpts = [[0.0, 0.5, 0.0],            # M
        [0.0, 0.0, 0.0],            # G
        [1./3, 1./3, 0.0],          # K
        [0.0, 0.5, 0.0]]            # M

# create band path from the high-symmetry points, 30 points inbetween each pair
# of high-symmetry points
kpath = make_kpath(kpts, nseg=30)
K_in_sup = []
for kk in kpath:
    kg, g = find_K_from_k(kk, M)
    K_in_sup.append(kg)

# remove the duplicate K-points
reducedK, kmap = removeDuplicateKpoints(K_in_sup, return_map=True)

if not os.path.isfile('KPOINTS'):
    # save to VASP KPOINTS
    save2VaspKPOINTS(reducedK)

第九行，kpts为能带路径。

然后计算非自洽，生成WAVECAR。

计算完成以后，使用如下脚本后处理。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import numpy as np
from procar import procar
from unfold import unfold, EBS_scatter
from unfold import make_kpath, removeDuplicateKpoints, find_K_from_k, save2VaspKPOINTS


# The tranformation matrix between supercell and primitive cell.
M = [[3.0, 0.0, 0.0],
     [0.0, 3.0, 0.0],
     [0.0, 0.0, 1.0]]

# high-symmetry point of a Hexagonal BZ in fractional coordinate
kpts = [[0.0, 0.5, 0.0],            # M
        [0.0, 0.0, 0.0],            # G
        [1./3, 1./3, 0.0],          # K
        [0.0, 0.5, 0.0]]            # M

# basis vector of the primitive cell
cell = [[ 3.1903160000000002,    0.0000000000000000,    0.0000000000000000],
        [-1.5951580000000001,    2.7628940000000002,    0.0000000000000000],
        [ 0.0000000000000000,    0.0000000000000000,   30.5692707333920026]]

# create band path from the high-symmetry points, 30 points inbetween each pair
# of high-symmetry points
kpath = make_kpath(kpts, nseg=30)
K_in_sup = []
for kk in kpath:
    kg, g = find_K_from_k(kk, M)
    K_in_sup.append(kg)

# remove the duplicate K-points
reducedK, kmap = removeDuplicateKpoints(K_in_sup, return_map=True)

if not os.path.isfile('KPOINTS'):
    # save to VASP KPOINTS
    save2VaspKPOINTS(reducedK)

if os.path.isfile('WAVECAR'):
    if os.path.isfile('awht.npy'):
        atomic_whts = np.load('awht.npy')
    else:
        p = procar()
        # The atomic contribution of Ce, Mo, S to each KS states
        # index starting from 0
        atomic_whts = [
            p.get_pw(0)[:,kmap,:],
            p.get_pw("1:18")[:,kmap,:],
            p.get_pw("18:54")[:,kmap,:]
        ]
        np.save('awht.npy', atomic_whts)

    if os.path.isfile('sw.npy'):
        sw = np.load('sw.npy')
    else:
        WaveSuper = unfold(M=M, wavecar='WAVECAR')
        sw = WaveSuper.spectral_weight(kpath)
        np.save('sw.npy', sw)

    EBS_scatter(kpath, cell, sw,
                atomic_whts,
                atomic_colors=['blue', "red", 'green'],
                nseg=30, eref=-1.0671,
                ylim=(-3, 3), 
                kpath_label = ['M', 'G', "K", "M"],
                factor=20)

通过调节参数可以知道导出不同原子的权重。

该脚本为二硫化钼中掺杂Ce原子，第一个原子是Ce，2-18个为Mo原子，19-54是S原子，如果不用分原子的导出，只想导出总的能带，只要把43-54，以及64-65行注释掉即可。

参考：

1.https://github.com/QijingZheng/VaspBandUnfolding

2.http://staff.ustc.edu.cn/~zqj/posts/Band-unfolding-tutorial/

编辑：黄琦 

为满足更多科研工作者的需求，蔻享平台开通了各科研领域的微信交流群。进群请添加微信18019902656（备注您的科研方向）小编拉您入群哟！

蔻享网站www.koushare.com已开通up主自主上传功能，期待您的分享！

---

欢迎大家提供各类学术会议或学术报告信息，以便广大科研人员参与交流学习。

联系人：李盼 18005575053（微信同号）