canvas标签做统计图：VASP计算并绘制三维能带图教程

火客 2023-01-31 01:18:57 596

canvas标签做统计图：VASP计算并绘制三维能带图教程则KPOINTS可以写成下面的样子：（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)简单来说，只要在布里渊区扫点算能带就可以了，在二维能带里，我们只想看沿着高对称路径的点上的能带色散关系。在三维能带图中，我们需要计算我们感兴趣的布里渊区域的所有点。比如，我们可以把一个小的区间里布点201*201个点。问题是，VASP不支持计算太多的点，解决方法就是把这些K点分开算。这里举个Dirac材料的例子，因为Dirac cone在K（1/3 1/3 0）点，我们想看三维狄拉克锥，需要把这个点包含在我们计算的能带中，首先我们需要把布里渊区需要的点确定好。比如我们需要的空间为

有时候大家需要就算某一能量区间的三维能带图，来了解某一区域的能带色散情况，比如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.0

0.01 0.0 0.00.01 0.5 0.0

0.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.closeKf.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