vasp画能带图，vasp能带图怎么画

有时候大家需要就算某一能量区间的三维能带图，来了解某一区域的能带色散情况，比如Dirac cone、 nodal-line之类的奇特能带结构。

这里小编给大家分享之前自己在研究Dirac cone的性质时候写的脚本，其实当时自己写了好多个子脚本，为了大家使用方便，这里小编尽量把脚本弄的简单一些，让读者用起来方便

首先讲一下原理：

其实三维能带就是二维能带图的叠加，就像切土豆片一样，每一个土豆片是一个二维能带图，把它们合在一起就是一个三维能带图

简单来说，只要在布里渊区扫点算能带就可以了，在二维能带里，我们只想看沿着高对称路径的点上的能带色散关系。在三维能带图中，我们需要计算我们感兴趣的布里渊区域的所有点。比如，我们可以把一个小的区间里布点201*201个点。问题是，VASP不支持计算太多的点，解决方法就是把这些K点分开算。

这里举个Dirac材料的例子，因为Dirac cone在K（1/3 1/3 0）点，我们想看三维狄拉克锥，需要把这个点包含在我们计算的能带中，

首先我们需要把布里渊区需要的点确定好。比如我们需要的空间为

（0.0 0.0 0.0） ... ... （0.0 0.5 0.0）: :（0.5 0.0 0.0） ... ... （0.5 0.5 0.0）

如果K点比较少，可以写成一个文件 (比如0到0.5之间布点51个，则每段之间的间隔为0.5/50=0.01)

则KPOINTS可以写成下面的样子：

k-points along high symmetry lines51Line-moderec0.0 0.0 0.0 0.0 0.5 0.00.01 0.0 0.00.01 0.5 0.00.02 0.0 0.00.02 0.5 0.0: : : 这里省略了中间的部分: : :0.5 0.0 0.00.5 0.5 0.0

如果布点比较多，就分开写。（小编已将脚本分享在文末）

准备工作：读者准备以下变量：

1. K点区域比如四方相的布里渊区里：（0.0 0.0 0.0） ... ... （0.0 0.5 0.0） : :（0.5 0.0 0.0） ... ... （0.5 0.5 0.0）2. 布点密度n，如果301*301布点，则n=3013. 能带指标，比如看Dirac cone需要价带顶（VBM），导带底（CBM）两个能带，（如果在你的体系中VBM，CBM分别是16，17）则变量nb = 2, nb1 = 16， nb2 = 174. 分成m个文件来计算，比如K点太多，需要分成十个KPOINTS来算 m=10

脚本file-pack.py变量含义如下：

Lb=[0.0, 0.0, 0.0] #K区间左下角坐标Lt=[0.0, 0.5, 0.0] #K区间左上角坐标Rb=[0.5, 0.0, 0.0] #K区间右下角坐标Rt=[0.5, 0.5, 0.0] #K区间右上角坐标n=301 #K点密度nb=2 #绘制能带数目nb1=16 #第16条带nb2=17 #第17条带m=10 #分成10个文件Ef=0.3682 #费米能级，这里是手动指定，下面有个提取E-fermi的命令，如果不想手动给，就使用下面的命令：Ef=`grep E-fermi OUTCAR |awk '{print $3}'|cut -d- -f 2` #(注意这一行跟上一行，注释其中一行，选一种方法给处Ef值就可以了，读者需要修改费米能级，就用上面一行)

接下来开始上脚本，第一个 file-pack.py 用来产生KPOINT以及提取数据的sh脚本

n=301 nb=2 nb1=9 nb2=10 m=10Ef=0.3682totbnd=16Lb=[0.0, 0.0, 0.0] Lt=[0.0, 0.5, 0.0] Rb=[0.5, 0.0, 0.0] Rt=[0.5, 0.5, 0.0] import numpy as nphead = \"k-points along high symmetry lines\n\" + str(n) + \"\nLine-mode\nrec \n\"#wirte the K-path to the file \"KFILE\"d13=[Rb[i]-Lb[i] for i in range(3)]d24=[Rt[i]-Lt[i] for i in range(3)]n=int(n)Kf=open(\"KFILE\",'w+')for i in range(n): A=[Lb[j]+i*d13[j]/(n-1) for j in range(3)] B=[Lt[j]+i*d24[j]/(n-1) for j in range(3)] Kf.writelines(str(A[0])+' '+ str(A[1])+' '+ str(A[2])+\"\n\") Kf.writelines(str(B[0])+' '+ str(B[1])+' '+ str(B[2])+\"\n\") Kf.writelines(\"\n\")Kf.close()Kf=open(\"KFILE\",'r+')kp_num = n//mlastline_num = n % mmm = kp_num * 3lm = lastline_num * 3for kk in range(m): f=open(\"KPOINTS_\"+str(kk),'w+') f.writelines(head)  for k in range(mm): line=Kf.readline() f.writelines(line)  f.close()f=open(\"KPOINTS_\"+str(m-1),'a+')f.writelines(head) for kk in range(lm):  line=Kf.readline() f.writelines(line) f.close()Kf.close()#prepare for data-export.shf=open(\"data-export.sh\",'w+')f.writelines(\"n=\"+str(n)+\"\n\" \+\"nb=\"+str(nb)+\"\n\" \+\"nb1=\"+str(nb1)+\"\n\" \+\"nb2=\"+str(nb2)+\"\n\" \+\"m=\"+str(m)+\"\n\" \+\"Ef=\"+str(Ef)+\"\n\" \+\"totbnd=\"+str(totbnd)+\"\n\" \+\"d=$((m-1))\n \Ef=`grep E-fermi OUTCAR |awk '{print $3}'|cut -d- -f 2`\n \for z in `seq 0 $d` \n \do \n \grep ^band $z/PROCAR |awk '{if(NR % '$totbnd'=='$nb1')print $5-'$Ef'}' >>VBM \n \grep ^band $z/PROCAR |awk '{if(NR % '$totbnd'=='$nb2')print $5-'$Ef'}' >>CBM \n \done \n \\")f.close()

比如我们需要做一个301*301K点3D能带图，需要分成十次算能带。1. 首先做静态计算（IBRION=-1; NSW=0; LCHARG=.TRUE.）

2. 把脚本复制到该目录下，

运行命令 python file-pack.py 产生10个KPOINTS_{1..10}文件，以及一个data-export.sh (用来提取能带能量; 临时文件KFILE)

修改INCAR （ICHARG=11; ISMEAR=0）（能带计算的INCAR参数）

3. 运行命令sh vasp-pack.sh 将10个KPOINTS文件以及其他计算能带的文件复制到0..9目录下

4. 提交计算（0..9目录分别提交）

5. 能带计算完成后，运行命令 sh data-export.sh 产生VBM，CBM两个能带的数据

6. 将VBM，CBM，surface-3D.py 文件复制到同一目录， 打开anaconda 中的spyder运行surface-3D.py

最后产生图 3d.png，其中surface-3D.py 中几个可调参数：

num=301 ！K点密度min = -2 ！能带图能量范围最小值max = 2 ！能带图能量范围最大值：：：ax.view_init(elev=5,azim=30) ！能带图倾斜角度

脚本分享：https://pan.baidu.com/s/1x3IQChoGAqC-1U8A-Ni7Xg 提取码: pz83 脚本中小编已经给了VBM，CBM 数据例子，可以先直接在spyder中测试surface-3D.py。参考小编的一篇文章：DOI: 10.1103/PhysRevB.102.165404