计算化学公社

标题: 【更新1.1版本】二维材料的二维能带图绘制以及有效质量计算软件band2d.py [打印本页]

作者
Author: ggdh 时间: 2017-3-10 16:28
标题: 【更新1.1版本】二维材料的二维能带图绘制以及有效质量计算软件band2d.py
本帖最后由 ggdh 于 2021-8-22 20:24 编辑

听说2D材料最近要火？我博后期间做这个写了个程序。关键是能够产生2D能带图，1维材料只能产生1维能带图，2维材料可以产生2维能带和1维能带。看2维能带能够直观掌握整个布里渊区的信息。而3维材料的3维能带图估计只有用VR跑到布里渊区里面才能看个contour版的3D能带图。。。现在把这个脚本分享出来，亲生的脚本希望能产生更多价值。

需求：
vasp
python 3.8.3
mayavi 4.7.3

1.1版新特性
1.1版主要增加了一些风格的控制
1. --shift 选项位移能带，--shift homo/fermi/0.04 分别是把vbm位移到0，fermi位移到0，以及位移0.04eV
2. --style3d 选项可以设置2维能带图的一些风格，包括是否显示边框，坐标轴，填色风格，尺寸大小等等（具体用法见帮助）
3. --rangeC 选项设置contour图的色彩范围（具体用法见帮助）
4. --contour 选项设置二维contour填色图等高线的数量和位置（具体用法见帮助）
5. -j pc 配合 --range3d 可以显示某个能带的contour填色图

主要功能：
对于二维材料:
a.计算1维能带结构，并绘制如下所示的图(左边是能带图，右边是K路径和布里渊区示意图）：
(, 下载次数 Times of downloads: 223)

b.计算2维能带结构，并绘制如下所示的图（这里只绘制了top of valence band）：
(, 下载次数 Times of downloads: 190)

c.计算有效质量，并在能带平面图上标出方向（分别是空穴和电子的有效质量）：
(, 下载次数 Times of downloads: 194)

安装说明:
先安装anaconda作为python环境
https://www.continuum.io/downloads
这里装python3.8.3版本的，可以装windows或者linux版本
安装好之后再linux Termial 或者是windows cmd 里面输入：
pip install mayavi
pip install PyQt5

使用说明：
1. 准备工作
用vasp优化好你的二维结构，或者是拿优化好的结构算个单点，例子输入文件如下，注意你的slab必须在xy平面，而z方向是真空层，同时晶体的ab轴在xy平面内，c轴和z平行：

INCAR文件如下：
PREC=Normal
ISMEAR=0
SIGMA=0.05
LREAL=AUTO
NSW=100
IBRION=2
ISIF=4

KPOINT文件如下：
0
Gamma
4 4 1
计算完成后产生的OUTCAR IBZKPT文件会被本程序读取。

2.检查当前情况以及设置条件
将band2d.py 拷贝到刚才vasp计算的文件夹下，运行
python band2d.py -j c
或者在linux 下把band2d.py 放到PATH路径下，直接运行
band2d.py -j c
对当前目录进行检查，并产生输出如下（红色是我的说明文字，不是输出的部分）
10 k points found in OUTCAR  #在OUTCAR中找到10个K点
no mesh_band data found #没有找到网格数据，网格数据是用来产生2D能带的
no mass_points data found #没有找到质量数据，质量数据是用来计算有效质量的
no points detected for peripheral band #没有找到外周能带数据点，外周能带就是布里渊区的边
no points detected for radical band #没有找到径向能带数据点，径向能带就是从Gamma点到布里渊区边缘的能带
0 k points found in 2D data #2维数据里面找到了0个K点，2维数据就是上面的mesh,mass,peripheral,radical点之和，因为这一步还没有计算，所以没有找到任何点。
Using eigen values from OUTCAR to determine the following infor #从OUTCAR里面的K点数据来计算下表中的数据
Gap    HOMO    LUMO    H_width L_width H_kpt             L_kpt
1.1188 -4.9184 -3.7996 0.1497 0.2720 0.0,0.0          0.0,0.0 #分别是能带，HOMO，LUMO，价带宽度，导带宽度，价带顶点坐标，导带底坐标（这两个坐标程序会自动选取价带顶和导带底，也可以手动选。有可能2d能带算了以后会找到一个更高的价带顶，那么就需要重新算有效质量）
step size for h, e are 0.060, 0.035 with delte E 0.020000 #计算有效质量的步长，因为价带宽度小一些（平滑一些），所以改变相同的能量（0.02）所需要的步长就大一些，用-s选项设置能量变化
angle between a b is 59.9979
length of a, b is 0.4273, 0.4273
hex reciprocal lattice detected with precision 1 #根据ab和夹角检测到是一个正六变形的晶胞, 这里对称性高的话，下面所需要算的k点数量也少，可以用-p选项提高precision 来强制降低对称性
number of kpts along a is 8 with step size (-G) of 0.050000 #根据间隔是0.05计算出沿着a轴要算8个k点，这部分k数量点包括了peripheral 和 radical，决定了1维能带是否平滑，用-G选项改变间隔
number of grid pt in x is 4 with step size (-g) of 0.050000 #根据间隔是0.05 x和y方向需要各算4个k点，这部分k数量点包括了mesh，决定了2维能带是否平滑，用-g选项改变间隔
number of grid pt in y is 4 with step size (-g) of 0.050000

另外还可以产生一个图片如下，其中的小黑点表明了mesh k点的密集度
， (, 下载次数 Times of downloads: 192)

3.产生输入文件
python band2d.py -j i
然后程序会创建新的文件夹，复制WAVECAR和CHGCAR，产生INCAR，KPOINTS
完成后会发现当面目录下多了一串 2D开头的文件夹
(, 下载次数 Times of downloads: 159)

4.计算
cd到每个2D开头的目录下，运行vasp计算。这里我是用脚本批量提交到服务器上算。输入vasps 2D* 就行。但是每个服务器的作业调度系统不一样。所以这个脚本要根据具体的服务器来编写。

5.产生结果文件
算完以后将所有的文件拷贝到本地机器上，一般我是把CHGCAR和WAVECAR 删了再拷。
for i in 2D*
do
rm $i/WAVECAR
rm $i/CHGCAR
done
然后在当前目录（当前目录包括了所有2D开头的文件夹以及之前单点/优化计算的结果）下运行
python band2d.py -j p
在终端会产生下面的输出：
10 k points found in OUTCAR
Warning! Number of points (12) in mesh line 0 not equal to 0.
Change the grid_np_x to 12 now
12 X 4 points detected for mesh grid
25 X 2 points detected for eff. mass grid
7 points detected for peripheral band
11 13 points detected for radical band
129 k points found in 2D data #可以看到这里找到了很多新算的k点
Using 2D data to determine the following infor #采用2Dk点来决定下面的属性，可以看到导带和价带都变宽了，这是布里渊区取样更全导致。
Gap    HOMO    LUMO    H_width L_width H_kpt             L_kpt
1.1188 -4.9262 -3.8074 0.2324 0.4061 0.0,0.0          0.0,0.0
step size for h, e are 0.040, 0.025 with delte E 0.020000
angle between a b is 59.9979
length of a, b is 0.4273, 0.4273
hex reciprocal lattice detected with precision 1
number of kpts along a is 8 with step size (-G) of 0.050000
number of grid pt in x is 12 with step size (-g) of 0.050000
number of grid pt in y is 4 with step size (-g) of 0.050000
12 X 4 points detected for mesh grid
25 X 2 points detected for eff. mass grid
7 points detected for peripheral band
11 13 points detected for radical band
然后会依次产生开头出现的那四个图。调整大小后保存就可以了。另外会产生如下的数据文件
DATA_BAND1D.out  绘制1维能带的数据文件（包含所有的能带）
DATA_BAND2D.out  绘制2维能带的数据文件（只包含当前绘制的能带，在--range3d里面设置）
DATA_MASS_K.out  用来计算有效质量的K点能量数据
DATA_MESH_K.out  用来产生2维能带的K点能量数据
DATA_PERI_K.out  用来产生1维外周能带的K点能量数据（外周能带是指那些不经过Γ点的K路径）
DATA_RADI_K.out 用来产生1维径向能带的K点能量数据（径向能带是指那些经过Γ点的K路径）
注意：
使用band2d.py -h 查看可用的选项和帮助
下面举一个使用例子
band2d.py -j p -g 0.03,0.03 -G 0.04 --range2d h5,l2 --range3d h2,h=-5,0  -p  0120340
说明：-j p 运行处理数据的任务，p for process
-g 0.03,0.03 对于2维能带，在x,y方向上的取样间隔是0.03
-G 0.04 对于1维能带，取样间隔是 0.04 （单位都是angstrom^-1)
--range2d h5,l2对于1维能带，绘制HOMO-5 到 LUMO+5 的能带
--range3d h2,h=-5,0 对于2维能带，绘制HOMO-2 和HOMO的能带，然后把z坐标的范围改成-5到0eV
-m h1,h,e 计算HOMO-1，HOMO，LUMO 三条能带的有效质量
-p 0120340 对于1维能带，k点的路径是从Γ点（0）出发，到外周的第一个点（1），然后逆时针转到外周第二个点（2），然后回到Γ，然后到外周的第三个点（3）。。。

遇到bug，或者有使用问题，或者有新功能建议在本帖留言或者加本人QQ：32589927。

作者
Author: 不取俗名 时间: 2017-3-11 16:51
很感兴趣，楼主介不介意放一两个例子，有相关的2D DOS的文章那就更好了，多谢

作者
Author: helpme 时间: 2017-3-11 19:38
这个估计会用到，先备份一个。

作者
Author: lindlar 时间: 2017-3-12 19:56
谢谢，学习啦

作者
Author: luckstar18 时间: 2017-3-27 04:36
good, vary good!
Thanks for sharing!

作者
Author: lindlar 时间: 2017-4-7 21:14
先收藏了，非常好的应用，谢谢啦！

作者
Author: solarman 时间: 2017-4-9 10:18
很好的东西，赞一个，不过应该称之为二维布里渊区的三维能带图更合适些？

作者
Author: 追枫少年 时间: 2017-5-10 09:55
您好，我是编程的小白，在尝试您的脚本的时候提示如下错误：
File "band2d.py", line 731
for i in range(len(self.radi_frac))}.copy()
^
SyntaxError: invalid syntax
望解答是什么原因？先拜谢了

作者
Author: ggdh 时间: 2017-5-10 17:52

追枫少年发表于 2017-5-10 09:55
您好，我是编程的小白，在尝试您的脚本的时候提示如下错误：
File "band2d.py", line 731
for i in ...

你装了anaconda没有
具体哪一步运行出错的？用 python band2d.py 的方式运行试试看

作者
Author: 追枫少年 时间: 2017-5-10 19:44

ggdh 发表于 2017-5-10 17:52
你装了anaconda没有
具体哪一步运行出错的？用 python band2d.py 的方式运行试试看

我安装的是anaconda2 版本对这个脚本有影响吗？
我尝试了 Python band2d.py的命令还是提示相同的错误

作者
Author: ggdh 时间: 2017-5-11 23:23
本帖最后由 ggdh 于 2017-5-11 23:27 编辑

追枫少年发表于 2017-5-10 19:44
我安装的是anaconda2 版本对这个脚本有影响吗？
我尝试了 Python band2d.py的命令还是提示相同的错 ...

anaconda2 是没有问题的
你试试输入
python --version 看看版本是不是2.7.12
另外使用dos2unix band2d.py
另外重新下载band2d.py
最后可以加我qq 32589927 ，有问题可以实时跟我反馈

作者
Author: 追枫少年 时间: 2017-5-13 09:08

ggdh 发表于 2017-5-11 23:23
anaconda2 是没有问题的
你试试输入
python --version 看看版本是不是2.7.12

非常感谢得到您的指点

作者
Author: shuihuntian 时间: 2017-5-29 09:09
谢谢楼主分享，请问能处理hse06的结果（计算有效质量）吗？

作者
Author: ggdh 时间: 2017-5-31 15:13

shuihuntian 发表于 2017-5-29 09:09
谢谢楼主分享，请问能处理hse06的结果（计算有效质量）吗？

可以

作者
Author: luobo123456 时间: 2017-6-6 13:01
谢谢楼主分享请问可以把图片的数据导出来吗？

作者
Author: ggdh 时间: 2017-6-7 10:38

luobo123456 发表于 2017-6-6 13:01
谢谢楼主分享请问可以把图片的数据导出来吗？

目前不行。需要的话，可以加这个功能。。

作者
Author: hakuna 时间: 2017-6-7 17:45
感谢原创分享，学习了

作者
Author: biao.biao 时间: 2017-7-12 11:23
you did a very nice work

作者
Author: biao.biao 时间: 2017-7-12 17:30
三维体系的材料可以用吗？，我的计算大多是三维的~

作者
Author: ggdh 时间: 2017-7-14 18:13

biao.biao 发表于 2017-7-12 17:30
三维体系的材料可以用吗？，我的计算大多是三维的~

不好用，因为这里不支持三维布里渊区

作者
Author: 重新再来668 时间: 2017-9-1 22:39
急需，很及时，但是软件安装老出错，不知道数据能不能导出来？

作者
Author: ggdh 时间: 2017-9-4 18:48

重新再来668 发表于 2017-9-1 22:39
急需，很及时，但是软件安装老出错，不知道数据能不能导出来？

已经更新，现在K点数据还有绘图的坐标数据都能自动产生了

作者
Author: ggdh 时间: 2017-9-4 18:49

biao.biao 发表于 2017-7-12 17:30
三维体系的材料可以用吗？，我的计算大多是三维的~

如果你只关心三维材料的2个维度，那还是可以用的。调整晶体结构。把你关心的两个维度放到xy平面上即可

作者
Author: wuhuawujiu 时间: 2017-12-20 09:01
非常好的软件！2维能带能否改变其中心，现在是gamma

作者
Author: ggdh 时间: 2017-12-20 14:41

wuhuawujiu 发表于 2017-12-20 09:01
非常好的软件！2维能带能否改变其中心，现在是gamma

改变中心的目的是？

作者
Author: wuhuawujiu 时间: 2017-12-21 08:07

ggdh 发表于 2017-12-20 14:41
改变中心的目的是？

有一些二维体系的价带顶和导带底不是在gamma点处，比如K点、M点等或者一些非高对称点(石墨烯，或一些node-line体系等)。想画这些体系的3D band，中心可以移动会更好一些。

作者
Author: zhuceid02 时间: 2018-3-23 09:38
你好，请问这个计算单点指的是价带顶和导带底必须是在Gamma的情况吗？其他情况不能算吗？还是所有的情况都可以？

作者
Author: ggdh 时间: 2018-3-24 22:55

zhuceid02 发表于 2018-3-23 09:38
你好，请问这个计算单点指的是价带顶和导带底必须是在Gamma的情况吗？其他情况不能算吗？还是所有的情况都 ...

所有情况都可以

作者
Author: maxwong2018 时间: 2018-4-25 10:06
好东西，感谢分享。刚开始做2D材料。

作者
Author: maxwong2018 时间: 2018-4-28 22:37
楼主，有个问题，想请教一下，为什么第二步自动产生的文件无法进行vasp计算？
我查看了一下，自动产生的K点文件实际上就是IBZKPT里的K点，但是数量又不一致，这样导致进行非自洽计算时，不能计算。这是怎么回事？

作者
Author: chengfang 时间: 2018-6-11 08:01
好漂亮的工作！

作者
Author: baibing111 时间: 2019-1-24 05:17
感谢无私奉献的楼主！！！

作者
Author: hakuna 时间: 2019-1-24 10:50
快两年了，又被翻出来了，那个水母图看起来还是那么性感

如楼主所言，2D的确火了

作者
Author: 张丽 时间: 2020-10-29 23:22
好贴

作者
Author: The_only_oneZLY 时间: 2021-5-31 22:25
大佬你好，我运行的时候报了如下的错，我用的是Python3.6，报错提示说指标超出边界，我修改了一下代码，但也没解决，想请教一下楼主是哪里的问题
Import  error
Spin polarized system detected. HOMO orbital is 30
120 k points found in OUTCAR
no mesh_band data found
no mass_points data found
no points detected for peripheral band
no points detected for radical band
0 k points found in 2D data
Using eigen values from OUTCAR to determine the following infor
Traceback (most recent call last):
  File "band2d.py", line 1433, in <module>
br.read_cell_param('OUTCAR')
  File "band2d.py", line 158, in read_cell_param
HOMO_band = np.array(eigval).T[self.HOMO-1]
IndexError: index 29 is out of bounds for axis 0 with size 24

作者
Author: The_only_oneZLY 时间: 2021-6-1 22:34

The_only_oneZLY 发表于 2021-5-31 22:25
大佬你好，我运行的时候报了如下的错，我用的是Python3.6，报错提示说指标超出边界，我修改了一下代码，但 ...

已解决，感谢楼主

作者
Author: 824475291 时间: 2021-7-19 22:27
请问，在哪下载啊？

作者
Author: zjxitcc 时间: 2021-7-20 08:26

824475291 发表于 2021-7-19 22:27
请问，在哪下载啊？

暴露了没仔细看帖子，就在帖子里。

作者
Author: 824475291 时间: 2021-7-21 08:45

zjxitcc 发表于 2021-7-20 08:26
暴露了没仔细看帖子，就在帖子里。

谢谢！

作者
Author: The_only_oneZLY 时间: 2021-10-12 13:11
本帖最后由 The_only_oneZLY 于 2021-10-12 13:13 编辑

(, 下载次数 Times of downloads: 86) (, 下载次数 Times of downloads: 94) trans_up.zip
楼主您好，最近在弄一个布里渊区的透射谱图，自己得到的数据画出来是图1的样子，横纵坐标单位都是2pi/a，想通过转换数据画成图2六边形的样子，参考了楼主的脚本学到了很多，但自己折腾了很久也没弄出来图2的效果，所以想请教一下楼主，想请楼主指点一下或者提供一点思路，附件是从MATLAB里面保存下来的数据，行和列都是-0.5到0.5范围的数据，万分感谢楼主[抱拳]。

作者
Author: ggdh 时间: 2021-10-14 15:53

The_only_oneZLY 发表于 2021-10-12 13:11
trans_up.zip
楼主您好，最近在弄一个布里渊区的透射谱图，自己得到的数据画出来是图1的样子，横纵坐标单 ...

我的做法是作图的时候加一个遮罩（mask）把6边形外的数据遮住。
python是肯定有这个功能的（matplotlib）
matlab也应该有这个功能

作者
Author: The_only_oneZLY 时间: 2021-10-18 13:08

ggdh 发表于 2021-10-14 15:53
我的做法是作图的时候加一个遮罩（mask）把6边形外的数据遮住。
python是肯定有这个功能的（matplotlib ...

感谢楼主，大概明白这个思路了，我再去试试。

作者
Author: The_only_oneZLY 时间: 2021-10-18 20:34
本帖最后由 The_only_oneZLY 于 2021-10-18 21:26 编辑

ggdh 发表于 2021-10-14 15:53
我的做法是作图的时候加一个遮罩（mask）把6边形外的数据遮住。
python是肯定有这个功能的（matplotlib ...

感谢楼主~

作者
Author: 尾巴上 时间: 2021-11-15 04:31
请问2D_mesh_* 中的KPOINTS需要改吗？我执行2D_mesh_*中的vasp计算都中止了，错误信息是Error reading KPOINTS file

作者
Author: 尾巴上 时间: 2021-11-15 04:32

maxwong2018 发表于 2018-4-28 22:37
**** 作者被禁止或删除内容自动屏蔽 ****

请问你的问题解决了吗？我也遇到了这个问题

作者
Author: ggdh 时间: 2021-11-15 11:52

尾巴上发表于 2021-11-15 04:31
请问2D_mesh_* 中的KPOINTS需要改吗？我执行2D_mesh_*中的vasp计算都中止了，错误信息是Error reading KP ...

你看看你是不是用了4面体的k点，不支持使用4面体k点

作者
Author: 尾巴上 时间: 2021-11-17 07:41

ggdh 发表于 2021-11-15 11:52
你看看你是不是用了4面体的k点，不支持使用4面体k点

谢谢回复。

作者
Author: 尾巴上 时间: 2021-11-17 07:42

ggdh 发表于 2021-11-15 11:52
你看看你是不是用了4面体的k点，不支持使用4面体k点

嗯嗯，不换了K点，计算就没问题了，大写的赞👍！

作者
Author: oucheng 时间: 2023-6-12 16:39
请问一下，我的导带极小值是在K点的，但是有效质量m1和m2的交点为什么不在K点呢？尝试用-m 修改位置但是还是同样的结果

作者
Author: ggdh 时间: 2023-6-13 13:24

oucheng 发表于 2023-6-12 16:39
请问一下，我的导带极小值是在K点的，但是有效质量m1和m2的交点为什么不在K点呢？尝试用-m 修改位置但是还 ...

估计是因为数值误差，导致这个地方的能量比k点还要低。
你可以手动改一下输出OUTCAR，把这个点的能量改高一点。从而使得band2d找到正确的CBM点

作者
Author: oucheng 时间: 2023-6-14 10:50

ggdh 发表于 2023-6-13 13:24
估计是因为数值误差，导致这个地方的能量比k点还要低。
你可以手动改一下输出OUTCAR，把这个点的能量改 ...

谢谢您的回复，请问这个点的坐标应该是-j p过程中的L_kpt坐标吗？我这里并不是，我的L_kpt是-K点的坐标（0.6667，-0.3333），然后我现在大概从图上判断这个点的坐标（0.333~0.5，<0.3333）

作者
Author: oucheng 时间: 2023-6-14 11:41

ggdh 发表于 2023-6-13 13:24
估计是因为数值误差，导致这个地方的能量比k点还要低。
你可以手动改一下输出OUTCAR，把这个点的能量改 ...

好像明白了，一开始SCF的k点文件不包含K点，算完2D能带才找到这个最低点，电子的有效质量是按一开始的SCF结果来的，所以并不包含K点，有效质量要重算

作者
Author: 00gi 时间: 2024-5-31 00:31
本帖最后由 00gi 于 2024-5-31 00:49 编辑

楼主您好我在计算最后一步的时候出现了下面的问题没有出来最后的彩图麻烦楼主帮我看看是什么原因十分感谢
(, 下载次数 Times of downloads: 58)

作者
Author: 00gi 时间: 2024-5-31 00:53

00gi 发表于 2024-5-31 00:31
楼主您好我在计算最后一步的时候出现了下面的问题没有出来最后的彩图麻烦楼主帮我看看是什么原因十分 ...

感谢已经解决啦环境重新配置了一下~

作者
Author: oucheng 时间: 2024-11-1 01:03
楼主请问这个提示如何应对呢，我用-s尝试了不同的值还是会有

Warning! The difference between ELECTRON mass calculated by h and 2h is larger than 20% (m1:60.60%, m2:60.04%). You may want to try different delta E(current value is 0.0400 )

作者
Author: PBEsol 时间: 2024-12-3 11:59
感觉会用到，楼主牛皮

作者
Author: ggdh 时间: 2024-12-3 21:37

oucheng 发表于 2024-11-1 01:03
楼主请问这个提示如何应对呢，我用-s尝试了不同的值还是会有

Warning! The difference between ELECTRON ...

不用管。这是个警告而已
如果不同的-s值差别不大，随便取一个就行，保持你研究的多个体系-s值统一就行

作者
Author: @yang 时间: 2024-12-25 16:55
请问老师，我进行第一步的时候出现这样的情况应该如何解决

欢迎光临计算化学公社 (http://bbs.keinsci.com/)