使用Python脚本获取VASP能带计算数据

Posted on Edited on In Views: Symbols count in article: 4.2k

由于在宿舍蹲，自己的Coding能力与日俱增。最近刚写了一份提取VASP能带数据的Python脚本，代码已经提交到了我的Github仓库中，有需要的朋友自提。这篇文章主要是介绍一下代码的结构和代码的用法，可以作为手册使用。

1.准备工作

数据的提取主要是基于Pymatgen库来实现的，脚本语言是python3，以下简单列出脚本运行需要使用的环境：

  • Python 3.9 or later
  • pymatgen
  • numpy
  • pandas

2.使用脚本

使用脚本时需要注意以下几点：

  • 需要设置费米能，一般是提取scf计算中的费米能(因为scf计算得到的费米能比能带计算得到的费米能精度更高)，提取费米能也很容易：

grep E-fermi OUTCAR

  • 需要在脚本的同级目录下放入能带计算的K点路径KPOINTS文件和产生的vasprun.xml文件，作为初始输入文件
  • 根据计算得到的能带是否分为自旋向上和自旋向下两支决定使用哪个脚本，注意不能根据是否打开了自旋来区分，因为对于有的体系即使打开了自旋最后算出来的能带还是只有一支，这种情况仍然使用no spin版本

3.简要介绍一下代码的结构

从vasp原始输出文件中读入数据

原始数据的读入主要使用的是pymatgen.io.vasp.outputs中的Vasprun类，它是pymatgen库中专门用来读取vasprun.xml文件的类。vasprun.xml是vasp运行时产生的文件，有点类似日志文件，记录vasp计算过程的所有数据，包括态密度和能带。读取之后通过get_band_structure函数获得一个BandStructure对象，这里需要传入能带计算中的KPOINTS文件作为输入参数，来得到相应的k点路径。

获取关键变量

得到BandStructure对象之后，需要获得进一步提取所需的变量，BandStructure对象本身就是字典类型，以键值对的方式存储，提取非常方便，后续计算需要提取的几个主要变量是：

  • 倒晶格基矢，相应的键是"lattice_rec" -> “matrix”
  • 倒空间中的k点坐标（包含所有），是分数坐标的形式。相应的键是"kpoints"
  • 每条能带分支的信息，包含当前分支下start_index, end_index, 高对称点的名字。相应的键是"branches"
  • 每条能带上每个k点的能量，相应的键是"bands" -> “1"或”-1"(如果是能带分了自旋向上和自旋向下两支)

将对应的K点坐标投影到X轴上

要画能带图，需要解决的一个关键问题就是：要把三维倒空间中的K点投影到一维X轴上。投影的方法是根据它们在K空间中的坐标计算相应的距离，然后对每个点逐步累加，得出的总距离即对应该点在X轴上的距离。这里有一个技巧，即vasp能带计算过程中，对于每一个分支上的能带，k点应该是均匀分布的，所以我们只需要计算出这一分支的总距离，然后平分就能得到该分支中相邻k点之间的距离，而不需要对于每两个相邻的点都计算一下，这也是上面为什么要获得每条能带分支的信息。以下代码完成的是这一工作：

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
def GetKptDistance(ls_kpt1, ls_kpt2, rec_lat):
    '''
    根据两个k点的坐标得到它们在倒空间中的距离
    '''
    mat_rec_lat = np.array(rec_lat)
    vec_k1 = ls_kpt1[0]*mat_rec_lat[0]+ls_kpt1[1]*mat_rec_lat[1]+ls_kpt1[2]*mat_rec_lat[2]
    vec_k2 = ls_kpt2[0]*mat_rec_lat[0]+ls_kpt2[1]*mat_rec_lat[1]+ls_kpt2[2]*mat_rec_lat[2]
    k_distance = np.linalg.norm(vec_k2-vec_k1)
    return k_distance


def GetLabelPosiInAxisX(ls_brchs_info, rec_lat, ls_kpnts):
    '''
    获取每个k点在X轴上的位置，即是将三维k点坐标转化到一维空间坐标系
    '''
    brchs_count = len(ls_brchs_info)
    sum_distance = 0.0
    df_klabels = pd.DataFrame(data=ls_brchs_info)
    ls_labels_distance = []
    ls_kpnts_distance = []
    for brchs_index in range(brchs_count):
        dict_brch = ls_brchs_info[brchs_index]
        st_index = dict_brch['start_index']
        ed_index = dict_brch['end_index']
        brch_distance = GetKptDistance(ls_kpt1=ls_kpnts[st_index], ls_kpt2=ls_kpnts[ed_index], rec_lat=rec_lat)
        ls_labels_distance.append(sum_distance+brch_distance)
        for i in range(st_index, ed_index+1):
            kpt_distance = sum_distance+(i-st_index)*brch_distance/(ed_index-st_index)
            ls_kpnts_distance.append(kpt_distance)
        sum_distance = sum_distance+brch_distance
    df_klabels = df_klabels.assign(label_sum_distance=ls_labels_distance)
    df_klabels.to_csv("klable.csv")           ### klabel.csv文件存储了每个k点的符号和在X轴上的位置
    return ls_kpnts_distance

这一过程中还生成了klabel.csv文件，这一文件的作用是标出每条分支中的高对称点和相应的在X轴上的坐标。

将能带数据存入文件

开始已经获取了存储相应每条能带上每个K点的能量，同时上一步也获取了每个K点在X轴上对应的坐标，因此就可以将完整的能带数据存入文件中了，以下函数完成这一过程：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def WriteBandsEnergyToTxt(filename, dt_bands, ls_kpts, fermi_energy):
    '''
    把能带的能量数据写入到文本文件中,其中filename是文件名,dt_bands是读取到的能带文件，fermi_energy是费米能
    '''
    bds_count = len(dt_bands)
    kps_count = len(ls_kpts)
    with open(file=filename, mode='w') as bd_file:
        bd_file.write("kpos\tenergy(eV)\tenergy_with_zero_fermi_energy\n")
        for i in range(bds_count):
            bd_energy = dt_bands[i]
            for j in range(kps_count):
                bd_energy_kj = bd_energy[j]
                bd_energy_kj_zero_efermi = bd_energy_kj-fermi_energy
                distance_kj = ls_kpts[j]
                str_oneline = '{:.5f}\t{:.5f}\t{:.5f}\n'.format(distance_kj, bd_energy_kj, bd_energy_kj_zero_efermi)
                bd_file.write(str_oneline)
            bd_file.write("\n")

保存的文件结构类似以下表格：

band no. kpos(k点对应的位置) energy(能量) energy_with_zero_fermi_energy(费米能平移为0后每条能带对应的能量)
band no.1 kpoint no.1 e_kp_1 e_kp_1-efermi
... ... ...
kpoint no.n e_kp_n e_kp_n-efermi
band no... kpoint no.1 e_kp_1 e_kp_1-efermi
... ... ...
kpoint no.n e_kp_n e_kp_n-efermi

spin和no spin版本的区别

相较于no spin版本，spin版本的改动主要有以下几处：

  • vasprun.xml文件读入时设置separate_spins=True
  • 读取能带能量时键为"1"的数据对应自旋向上，键为"-1"的数据对应自旋向下
  • 处理能带数据分为自旋向上和向下分别处理，新增两列。

4.数据后处理

程序运行完成之后会生成两个文件：klabel.csv和bd.txt，其中klabel.csv文件存储的是每条分支能带的信息，包括每个高对称点处的坐标，可以用于后续画能带图。而bd.txt存储的就是每条能带上每个点的在X轴上的坐标和相应的能量大小。如下为一个示例：

  • klabel.csv
    klabel
  • bd.txt
    bands
    如下图是用origin通过以上数据绘制的能带图(将费米面平移到纵轴零值处)：
    originband
    为了验证结果的正确性，这里同样调用pymatgen.electronic_structure.plotter库的相应函数自动绘制了能带图，结果如下，是一致的，说明了提取数据方法的正确性。
    pymatgenband