2.1 电子性质1
自洽计算
VASP做静态自洽的目的主要是为了得到好的电子密度和波函数文件，为了使后续的性质计算可以读取电子密度和波函数，进而增加收敛速度。
首先进行 高精度 的结构优化。然后复制文件， cp CONTCAR POSCAR 做一次自洽计算。
自洽计算要求有密的k网格点，在计算条件允许的情况下要求自洽的k网格点大致为优化时的2倍左右。
修改标签
1
2
3
4IBRION=-1
NSW=0
LWAVE = .T. (Write WAVECAR )
LCHARG = .T. (Write CHGCAR )
这一步完成后的输出文件CHGCAR即是电荷密度，导入VESTA直接查看。
bader电荷
自洽时INCAR再加一个参数：
1
2
3IBRION=-1
NSW=0
LAECHG=.T.
输出文件：
AECCAR0 全电子中的芯电子部分（POTCAR）
AECCAR1 自洽计算前的交叠原子电荷（POTCAR，bader 分析不用）
AECCAR2 全电子中的价电子部分（自洽计算后）
CHGCAR 自洽计算后赝电子
输入：
1 | chgsum.pl AECCAR0 AECCAR2 |
得到两个文件：
- BCF.dat 电荷极大值序号和坐标
体积内的 Bader 电荷（积分）
距离极大值最近的原子和距离 - ACF.dat 原子序号和坐标
原子电荷（极大值的电荷之后）
到零通量面最小距离
原子体积
李强博士《learn vasp the hard way》https://www.bigbrosci.com/2011/12/24/A08/ Bader电荷与原子电荷着色：
https://www.bigbrosci.com/2011/12/23/A07/ https://mp.weixin.qq.com/s/BNzjhz8SI_HXkaEnyL29Bg
成键反键
自洽时添加：
1
2
3
4
5NBANDS=96
ISYM=-1
IBRION=-1
NSW=0
LWAVE = T
vi lobsterin
原子编号改为自己的，例如
lobsterin：
1 | COHPstartEnergy -14 |
前两行代表提取的能量范围；第三行代表选择的基组，如果不设置(空行)就是默认基组，如果要覆盖默认基集，可使用此关键字进行覆盖。它目前支持bunge、koga和pbeVaspFit2015；第四行代表为每个元素指定所使用的基函数；最后一行代表选择一对两个原子进行键合分析。
lobster < lobsterin
1、删除存在的WAVECAR文件，做单点计算，计算需要保存波函数文件；
2、做完单点计算在当前目录下运行lobster软件；
3、用wxDragon软件打开COHPCAR.lobster文件。
4、COHPCAR.lobster 文件包含COHP结果，第二列为pCOHP，将第二列乘-1可得-pCOHP
https://github.com/JaGeo/LobsterPy
https://zhuanlan.zhihu.com/p/470592188
注：
- orbitalWise 这个参数一般加在cohpBetween或者cohpGenerator之后，用以把原子之间的相互作用具体分到s p d 这些轨道上，比如计算s-s, p-p, d-d之间的pCOOP, pCOHP；
比如：cohpBetween atom 1 atom 2 orbitalWise
cohpGenerator from 1.4 to 1.5 type Ga type Sr orbitalWise
不要使用US-PP超软赝势，使用PAW赝势；请使用vasp_std版本，不要使用vasp_gam
；不支持处理对称:ISYM=0或ISYM=-1由于需要把基于平面波的波函数投影到原子轨道的基函数上，所以计算的能带的数量要多于投影的基函数，VASP计算的时候，默认的NBANDS达不到要求，需要加大NBANDS，只要查看OUTCAR里面的NBANDS多少，
grep NBANDS OUTCAR；在INCAR里面，NBANDS关键词设置成1.5倍NBANDS，足以满足计算的投影要求，在计算过程中一定要记住设置足够的NBANDS，否则会出现报错，一直不能正常算出COHPcohpGenerator 这里通过距离或者元素类别告诉想要计算什么原子之间的COHP；
cohpGenerator from 1.5 to 2.3 type Fe type N orbitalwise#输出Fe-N原子对距离在1.5到2.3的所有数据之和。gaussianSmearingWidth 高斯展宽，如果VASP中用的ISMEAR=0，那么gaussianSmearingWidth的数值和SIGMA一样，默认值是0.2 eV，但是过大了，所以需要减小。设置格式：
gaussianSmearingWidth 0.05
ELF电子局域函数计算
电子局域函数（Electron Localization Function—ELF）是研究电子结构的手段之一。用来描述以某个位置处的电子为参考，在其附近找到与他同自旋的电子的概率，可以表征这个作为参考的电子的局域化程度，也是一种描述在多电子体系中的电子对概率的方法。
ELF的取值范围是[0，1]。
ELF=1 对应完全局域化；
ELF=1/2，对应类电子气型的成对概率；
ELF=0，对应电子完全离域化。
计算流程：做自洽的时候在INCAR中加入：
1
2
3
4PREC=A
LELF=T
IBRION=-1
NSW=0
运行vasp得到ELFCAR文件。
将ELFCAR拖入Vesta http://www.jp-minerals.org/vesta/en/download.html
点击 Utiltties - 2D Data Display。
点击Slice，确定切面和fractional coordinates。
导出图片：之后可以用PS工具将原子模型加入。
费米面
费米面计算，这里我们以Cu为例，见vaspkit/examples/Cu_fermi_surface。
第一步，准备好POSCAR，注意一定是原胞（primitive cell ）和用于静态计算的INCAR文件；
第二步，运行VASPKIT，输入261命令，产生用于计算计算费米面的KPOINTS和POTCAR文件，
第三步，提交VASP作业；第四步，待VASP计算完成后，再次运行VASPKIT，输入262命令，生成包含费米面数据的FERMISURFACE.bxsf文件；第五步，运行XcrySDen软件。
发现第五条能带穿过费米面，因此我们选择Band Number 5，然后点击Selected，然后得到Cu的费米面。推荐使用FermiSurfer可视化费米面。
对于半导体不存在费米面，如果想可视化导带或价带中某个能级处的等能面，新建FERMI_ENERGY.in文件然后在该文件第二行（第一行为注释行）中人为指定能级大小，注意该值一定要保证穿过价带或导带。接下来调用vaspkit-252或263命令提取可视化。
STM模拟
参考：https://zhuanlan.zhihu.com/p/538950799
https://zhuanlan.zhihu.com/p/667558468
https://www.vasp.at/wiki/index.php/STM_of_graphite
首先做一个自洽计算，输出 WAVECAR 和 CHGCAR 文件。
加入一些特殊参数计算，输出我们想要的 PARCHG 文件。
1
2
3
4
5LPARD = .TRUE.
LSEPK = .FALSE.
LSEPB = .FALSE.
NBMOD = -3
EINT = -0.1 0.1用p4vasp软件打开PARCHG文件或者vaspkit 后处理。由于 VASP 在对于远离原子表面时波函数的衰减指数描述不明确，所以模拟探针高度的时候不应该太高，否则会出现非物理的现象，比如 7 埃，就算比较远了。所以一般高度的选择要接近材料表面，这里选择了 0.5 埃。而原胞沿着 ab 平面重复的标准是 20 埃，也就是说原胞晶格 * 重复单元 = 20 埃就可以。vaspkit 里选择的是 Tersoff-Hamman 方法，现在大多数模拟 STM 的都是选择这个方法。
vaspkit处理后的输出文件包括X.grd, Y.grd, STM_*.grd三个文件，可以用以下脚本画图。
Usage: python plot.pyimport numpy as npimport matplotlib as mpl#mpl.use(‘Agg’) #silent modefrom matplotlib import pyplot as plt
x_mesh=np.loadtxt(‘X.grd’) y_mesh=np.loadtxt(‘Y.grd’) v_mesh=np.loadtxt(‘STM.grd’)
fig = plt.figure()cmap = plt.cm.get_cmap(“Greys_r”) plt.contourf(x_mesh,y_mesh,v_mesh,cmap=cmap)
plt.axis(‘equal’)plt.axis(‘off’)#plt.show()plt.savefig(“STM.jpg”,dpi=300)
或者用以下脚本处理后，得到data.txt文件，将该文件导入origin画图。Usage: sh get_data.shawk ‘{for(i=1;i<=NF;i++) print $i}’ X.grd >xawk ‘{for(i=1;i<=NF;i++) print $i}’ Y.grd >yawk ‘{for(i=1;i<=NF;i++) print $i}’ STM_0.50.grd >spaste x y s >data.txt
态密度
自洽计算之后，读取产生的CHGCAR、WAVECAR
上面例如bader电荷、ELF标签都关闭。
INCAR中修改：
1
2
3
4
5ISTART=1
ICHARG=11
ISMEAR=-5
NEDOS=2001
LORBIT=11
其他参数和自洽时保持一致。
运行后vaspkit的111/112/113/114命令处理即可得到。
分波态密度和局域态密度
态密度很多分析和能带的分析结果可以一一对应，很多术语也和能带分析相通。但是因为它更直观，因此在结果讨论中用得比能带分析更广泛一些。在电子能级为准连续分布的情况下，单位能量间隔内的电子态数目，即为态密度。能态密度与能带结构密切相关，是一个重要的基本函数。固体的许多特性，如电子比热、光和X射线的吸收和发射等，都与能态密度有关。在VASP中，LORBIT=10计算的就是LDOS，也就是每个原子的局域态密度 (local DOS)，是分解到每个原子上面的s,p,d；LORBIT=11，计算的就是PDOS,投影态密度(projected DOS)或分波态密度(partial DOS)，不仅分解到每个原子的s,p,d,而且还进行px,py,pz分解。
在 /vaspkit.0.72/examples/LDOS_PDOS/Partial_DOS_of_CO_on_Ni_111_surface 目录下启动 vaspkit ，输入命令 115 选择子功能 The Sum of Projected DOS for Selected Atoms and orbitals 。我们需要提取的是O的s，p轨道，C的s，p轨道以及表面的dxy和d3z2-r2(dz2)轨道。首先提示你选择元素（累加），第一次输入元素O，回车后提示你输入提取的轨道名（累加），第一次输入轨道s。回车后重复以上两个操作。如果想结束输入，在元素选择行直接按回车键结束输入。
元素行接受自由格式输入，你可以输入以下内容元素行接受自由格式输入，你可以输入以下内容1-3 4 Ni表示选择元素1,2,3,4和Ni元素的PDOS进行累加。轨道行只支持标准输入，如果使用了LORBIT=10不投影轨道，你只能选择 s p d f中的一个或多个，如果使用了LORBIT=11投影了轨道，你可以从s p d f和s py pz px dxy dyz dz2 dxz dx2 f-3 f-2 f-1 f0 f1 f2 f3中选择一个或多个输入。轨道行还支持输入all计算所有轨道的DOS之和。
比如元素行输入C O ,对应的轨道选择s px，那么提取的就是所有C和O元素的s轨道和px轨道的PDOS之和。 PDOS_USER.dat中的轨道名 为#Energy C&O_s&px，显而易懂。
电荷密度差
Vaspkit 可以方便的处理 VASP 输出的电荷密度文件 (CHGCAR)，做电荷密度差分处理，VTST 脚本或者 VESTA 也可以做，但是不如 vaspkit方便灵活。
注：此文中讲述的功能同样适用于任何的实空间函数形式的文件处理，比如静电势差值图 (LOCPOT)，电子定域化函数差值图 (ELFCAR)，能带分解的电荷密度插值图 (PARCHG) 等等。
电荷密度差分（charge density difference）是研究电子结构的重要手段之一。可以直观的得到两个片段相互作用后的电子流向，原子形成分子过程中电子密度的变化、探究化学键的本质。电荷密度差分有以下几种形式：
（1）体系的电荷密度减去组成它的两个或几个片段的密度，也是最常见的电荷密度差分形式：
（2）自洽计算收敛以后体系的电荷密度减去该原子构型下每个原子的球对称的电荷密度(即初猜电荷密度)，也称为变形电荷密度(Deformation charge density)
（3）在某个状态的密度减去这个体系在另外一个状态的密度。比如：外加电场作用下的电荷密度减去没有外势场的电荷密度。再比如：激发态的密度减去它在基态时的密度。
以上三种电荷密度差，无论哪一种计算都需要保证原子坐标坐标必须一致！并且保证晶胞参数和 FFT-mesh 的格点密度 (NGXF, NGYF, NGZF) 一致！
(1) 的计算只需要优化 AB 的结构，A、B片段保持 AB 结构中原子坐标，不能再对片段进行单独优化。
(2) 只要优化自洽计算时候的几何结构，生成原子的球对称的电荷密度时候 (ICHARG = 12) 保持坐标不变。
(3) 也需要共用其中一个状态的几何构型（比如激发态或基态的几何构型）。
CHGCAR文件格式
CHGCAR 是包含电子密度信息的格点文件，对于自旋非极化体系（ISPIN = 1）计算只包含电荷密度，对于自旋极化体系（ISPIN = 2）计算还包含自旋电子密度。可以用 VESTA，Jmol 等程序打开。CHGCAR 的第一部分和 POSCAR, CONTCAR 的格式是完全一样的，它包含了最终结构的晶格矢量，原子核坐标等信息。紧接着是三个数字，这三个数是实空间函数的网格密度，对应于 NGXF，NGYF，NGZF 三个变量。然后是电荷密度信息 ρ(r) * Vcell。一共有 NGXF * NGYF * NGZF 个数值。比如：下面是一个 Al2O3 晶胞的 CHGCAR 文件
1 | Al2O3 Cell opt |
Vaspkit计算几个片段的电荷密度差
以 CO 分子吸附在 Ni(100) 表面为例：
步骤一：优化 CO/Ni(100) 的结构；
步骤二：分别计算 Ni(100) 和 CO 的单点能 CO 和 Ni(100) 片段的坐标从 CO/Ni(100) 的 CONTCAR 里直接截取，不要再结构优化！！计算时也要保证三次自洽计算所采用的 FFT mesh 一致（NGXF,NGYF,NGZF）
步骤三：用 vaspkit 314 功能做电荷差分。依次输入三个片段的 CHGCAR 路径:
1 | 314 |
生成 CHGDIFF.vasp 包含电荷密度差的信息，可以直接导入到 VESTA 里作图。
图片
Vaspkit计算变形电荷密度差
变形电荷密度差：自洽计算收敛以后体系的电荷密度减去该原子构型下每个原子的球对称的电荷密度(即初猜电荷密度)。
以 CO 分子为例：
步骤一：先自洽计算优化 CO 分子。
步骤二：新键文件夹，在优化好的结构基础上用 ICHARG = 12 做非自洽计算。
步骤三：用 vaspkit 314 功能做电荷差分。依次输入两个 CHGCAR 路径。
1 | 314 |
步骤四：得到 CHGDIFF.vasp 文件，导入到 VESTA 里做图。
Vaspkit计算外加电场下的电荷密度差
外加电场作用下的电荷密度减去没有外势场的电荷密度。
以 InSe 二维单层材料为例：
步骤一：先优化没有外加电场的结构。
步骤二：在同样结构下计算外加电场下做单点自洽计算。添加关键词，EFIELD 控制电场力的大小（eV/Angstrom）F=qE， 对应电场的单位是 E=F/q， 场强单位是 V/Angstrom
1 | EFIELD = 0.05 |
步骤三：运行 Vaspkit 314 功能做电荷差分。依次输入二个 CHGCAR 路径。
1 | 314 |
电荷密度差作图
得到的电荷密度差文件必须作成图才便于被直观地考察，常用的图分为三类：
（1）3D 等值面图，isosurface.
（2）2D 切面图.
（3）1D，平面平均的数值，planar averaged charge.
通常这几种图可以配合使用讨论。等值面图和切面图用 VESTA 画最方便，平面平均图用 vaspkit 生成方便。得到的 CHGDIFF.vasp 用 VESTA 打开，isosurface 青色部分电荷密度减小，黄色密度电荷密度增加。注意 VESTA 对电荷密度的默认单位是 e/Bohr3。
vaspkit 319 功能可以把VASP的格点文件转化成 .cube的方法，.cube 在 VMD 里作图可以得到更好看的效果：
能带计算
参考链接：
https://yh-phys.github.io/2019/10/04/vasp-band-structure-dos/
自洽计算之后，读取产生的CHGCAR、WAVECAR
1
2ISTART=1
ICHARG=11
其他参数和自洽时保持一致。
能带计算的KPOINTS文件不同于前面，称为高对称点。
学术之友微信公众号历史记录搜索文章《五种方法生成能带结构计算的高对称点》
1、http://cryst.ehu.es/
2、https://www.materialscloud.org/work/tools/seekpath
3、http://www.aflowlib.org/aflow-online (已更新)
4、http://pymatgen.org/ (pymatgen)
5、http://www.xcrysden.org/ (XCrySDen)
Accurate DOS and Band-structure calculations：
http://cms.mpi.univie.ac.at/vasp/vasp/Accurate_DOS_Band_structure_calculations.html
可以通过vaspkit的303命令产生KPATH.in文件，cp KPATH.in KPOINTS即可。
计算完成后通过vaspkit的211命令提取能带。
注意能带路径生成操作需在自洽计算前完成。
1.准备POSCAR，调用vaspkit- 303（体相材料）或-302（二维材料）得到KPATH.in和PRIMCELL.in文件；对于二维体系，需要检查POSCAR文件的真空层是否沿z方向，如果没有，可调用vaspkit-923或vaspkit-407强制真空层沿z方向。另外，如果需要改变识别结构对称性精度（symprec默认值为1E-5），可通过vaspkit -task 302 -symprec 1E-6实现或者在~/.vaspkit中修改SYMPREC默认值。
2.执行cp PRIMCELL.vasp POSCAR命名。（如果PRIMCELL.vasp 和POSCAR不一样，并且要读取电荷密度和波函数，要从SCF计算保持结构一致。）KPATH.in推荐能带路径只针对于PRIMITIVE CELL，缺少这一步，你可能得到错误的结果。如果有必要，比较KPATH.in文件中的能带路径是否与在线能带路径产生工具SeeK-Path产生的一致，包括比较PRIMCELL.vasp和HIGH_SYMMETRY_POINTS文件。需要指出的是SeeK-Path只用于体相结构能带路径的产生。自1.3.0版本，对于二维体系不需要执行cp PRIMCELL.vasp POSCAR。
投影电荷密度
复制能带结构文件，
vi INCAR
1
2
3增加标签 LPARD=.T.
增加标签 IBAND=16 17
增加标签 LSEPB=.T.
- LPARD：它的取值为.TRUE.或.FALSE.，它的默认值是.FALSE.，当为.TRUE.时，表示读入自
洽收敛的 CHGCAR 和 WAVECAR 并进行 partial charge density计算。 - IBAND：它的取值就是你想要计算的第几条能带或哪几条能带(比如要计算第 4、 5、 6 条能带，那么就设置 IBAND = 4 5 6)，此时 NBMOD的值就是所要计算的能带的条数。
- EINT：另一种方式指定所要计算的能带，它是指定计算某能量范围的partial charge density，一般是设置为两个实数，比如EINT= 4.00 5.00，此时NBOMD的值设置为-2。如果只设置了一个数，那么表示计算从 EINT 到费米能级这个范围内的，此时NBOMD 的值为-3。
- KPUSE：指定所要计算的 k 点(哪个或哪几个)。
- LSEPB：指定是不是要把计算的partial charge density按每个带分别写到各自对应的文件
- PARCHG.nb (设置为.TRUE.)中，还是把它们合并写到一个文件中(相当于把各个带对应的partial charge density 加起来，设置为.FALSE.)，默认值为.FALSE.。
- LSEPK：指定是不是把要计算的partial charge density按每个k分别写到各自对应的文件PARCHG.nk (设置为.TRUE.).
三个方法：
1、读取自恰计算的波函数，KPOINTS文件采用特定的k点
1 | #partial charge density |
1 | CsPbI3-cubic |
http://muchong.com/html/200906/1404880.html
2、读取自恰计算的波函数，KPOINTS文件采用自恰计算的k点
1
2
3
4
5
6
7
8#partial charge density
LPARD = .TRUE.
IBAND = 17 18
KPUSE = 35
LSEPB = .TRUE.
LSEPK = .TRUE.
NELM = 200
ICHARG = 1
在scf自恰计算OUTCAR中
找到对应导带底和价带顶
对应的k-point编号
1
2
3
4
5CsPbI3-cubic
0
Gamma
8 8 8
0.0 0.0 0.0
3、读取能带计算的波函数，KPOINTS文件采用能带计算的k点
1 | #partial charge density |
KPOINTS和能带计算一样。在band能带计算OUTCAR中找到对应导带底和价带顶对应的k-point编号。
轨道投影态密度和能带。
1 | LORBIT=11 |
提取DOS：vaspkit -task 111/113
提取band：211/212/213/214
如果3D能带图不够平滑，可通过能带插值进行改善。
VASPKIT升级到1.3.0或更新版本，在~/.vaspkit文件中设置以下三个参数：
1
2
3GET_INTERPOLATED_DATA .TRUE. 是否插值
INTERPOLATION_SPACING 0.04 插值间距，决定平滑程度
INTERPOLATION_METHOD 'cubic' 插值方法
HSE06计算电子结构：
结构优化（PBE）
静态自洽scf计算（PBE）KPOINTS和POTCAR文件不变，写波函数和电荷密度
(可选）静态自洽scf计算（HSE06）KPOINTS，POSCAR和POTCAR文件不变
DOS/band计算（HSE06）vaspkit -task 303
vaspkit→251→2→0.03→0.03生成KPOINTS，运行VASP
vaspkit-252
http://vaspkit.cn/index.php/29.htmlPBE优化结构
HSE自洽计算(得到波函数)
HSE能带计算/态密度计算
{.line-numbers} 1
2
3
4
5
6ALGO = Damped
HFSCREEN = 0.2
LHFCALC = .True.
LMAXFOCK = 4
PRECFOCK = Fast
TIME = 0.4
+SOC计算电子结构：
- 结构优化（PBE）
- 静态自洽scf计算（PBE+SOC）KPOINTS和POTCAR文件不变，写波函数和电荷密度
- 非自洽计算（PBE+SOC）
Spin-Orbit Coupling Calculation{.line-numbers} 1
2
3
4
5
6
7
8
9ISPIN = 2
NELMIN = 6 (Min electronic SCF steps)
LSORBIT = .TRUE. (Activate SOC)
GGA_COMPAT = .FALSE. (Apply spherical cutoff on gradient field)
#VOSKOWN = 1 (Enhances the magnetic moments and the magnetic energies)
LMAXMIX = 2 (For d elements increase LMAXMIX to 4, f: LMAXMIX = 6)
SAXIS = 0 0 1 (Direction of the magnetic field)
#MAGMOM = 0 0 3 (Set this parameters manually, Local magnetic moment parallel to SAXIS, 3*NIONS*1.0 for non-collinear magnetic systems)
NBANDS = 32 (Set this parameters manually, 2 * number of bands of collinear-run)
调节SOC强度
- 调节自旋轨道耦合作用的大小
要调优VASP中的自旋轨道耦合强度，需要修改源文件vasp_source_code_path/src/relativistic.F 并重新编译。如果要将强度减少一半，可以在relativistic.F的129行L●S项乘以0.5d0，如下：
DO I=0,1DO J=0,1DO M =1,2LL+1DO MP=1,2LL+1 DLLMM(LMP+MP-1,LM+M-1,J+2I+1)=DLLMM(LMP+MP-1,LM+M-1,J+2I+1)+ & 0.5d0SUMLS(M,MP,I+2*J+1,LL) !!!line 129 relativistic.F fileEND DOEND DOEND DOEND DO
代码片段：可切换语言，无法单独设置文字格式
注意，需要使用重新编译的vasp_ncl来运行计算。另外，最好比较0d0L●S，1d0L●S和vasp_std，vasp_ncl计算结果。
- 如果只想给一些元素加SOC，而其他元素不加该怎么办呢？拓扑材料计算互助群的Eurus-清华给出了VASP中给特定元素加自旋轨道耦合相互作用的方法（感谢Eurus-清华）。在原来VASP，src文件夹下的relativistic.F文件中注释掉46行的
! REAL(q), PARAMETER :: INVMC2=7.45596E-6
然后53行位置加入如下6-17行的内容，然后重新编译。
INTEGER, PARAMETER :: LMAX=3, MMAX=LMAX*2+1 COMPLEX(q) DUMMY(MMAX,MMAX,3,LMAX) COMPLEX(q) LS(MMAX,MMAX,4,LMAX) REAL(q) SUM, SCALE
! add parameter REAL(q), PARAMETER :: INVMC2_ORIG=7.45596E-6 REAL(q) INVMC2! add end ! add control IF(ABS(Z-83).LT.1.0E-5)THEN INVMC2=INVMC2_ORIG ELSE INVMC2=INVMC2_ORIG*0.0001 END IF! add end
LS=(0._q,0._q) CALL SETUP_LS(1,THETA,PHI,DUMMY(1:3,1:3,1:3,1),LS(1:3,1:3,1:4,1)) CALL SETUP_LS(2,THETA,PHI,DUMMY(1:5,1:5,1:3,2),LS(1:5,1:5,1:4,2)) CALL SETUP_LS(3,THETA,PHI,DUMMY(1:7,1:7,1:3,3),LS(1:7,1:7,1:4,3))
原理：这样做的原理就是将除Z-83的元素外的其他元素缩小一万倍。注意：这个脚本的第12行, IF(ABS(Z-83))，表明除了83号元素打开SOC之外，其他所有的元素SOC都是相当于关闭的。如果想打开某个特定元素的SOC，就把83改成你想要的元素序号就行。
另外，其他元素SOC关闭方式就是将其缩小一万倍。所以使用这种方法，任然是非共线计算，如果要进行wannier90计算，只要按照正常的wannier90设置就可以。
HSE06+SOC band structure calculation
Step-I: HSE06 SCF calculation with SOC (INCAR):
{.line-numbers} 1
2
3
4
5
6
7
8
9
10#SOC-related:
LSORBIT=.TRUE.
SAXIS= 0 0 1
ISYM=0
#HSE06-related:
LHFCALC = .TRUE.
HFSCREEN = 0.2
ALGO = Damped #Damped/Normal/All
TIME = 0.4Step-II: HSE06 NSCF calculation with SOC:
{.line-numbers} 1
2
3
4
5
6
7
8
9
10#SOC-related:
LSORBIT=.TRUE.
SAXIS= 0 0 1
ISYM=0
#HSE06 related:
LHFCALC = .TRUE.
HFSCREEN = 0.2
ALGO = Normal !use Normal
TIME = 0.4
能带反折叠
vaspkit:
- 准备POSCAR
- 运行vaspkit -task 302/303生成KPATH.in文件
- cp PRIMCELL.vasp POSCAR
- 运行vasp优化结构，cp CONTCAR POSCAR
- 运行vaspkit -task 400/400生成超胞结构POSCAR，并执行cp TRANSMAT TRANSMAT.in
- 运行vasp优化含有缺陷的超胞结构，cp CONTCAR POSCAR
- 运行vaspkit -task 281生成KPOINTS文件
- 准备INCAR（静态计算，设置LWAVE = T，适当增加NBANDS值）并执行VASP计算
- 运行vaspkit -task 282提取effective band structure(画法类似于投影能带)
- 如果不会计算变换矩阵（TRANSMAT.in文件不存在时），可以执行cp PRIMCELL.vasp PRIMCELL.in，281 ，这时vaspkit分别读入超胞基矢（POSCAR）和原胞基矢（PRIMCELL.in）计算变换矩阵并生成TRANSMAT.in文件。
- vaspkit/examples/band_unfolding/ebs.py画图。生成的ENERGY.grd，WEIGHT.grd和MOMENTUM.grd三个文件，可使用vaspkit/examples/band_unfolding/ebs_k_resolved_dos_plot_matlalb.m脚本画图。
bandup:
第一步 得到超胞CHGCAR:输入文件参考BANDUP的example的第一步，提交任务之后得到CHGCAR。
第二步 主要是为了得到超胞非自洽计算的K点路径，分别给出单胞和超胞的晶格坐标文件prim_cell_lattice.in与supercell_lattice.in。还有单胞的高对称K点路径。
下一步命令是bandup kpts-sc-get -no_symm_avg,这一步会生成KPOINTS_supercell.out文件后面会用到。-no_symm_avg这个标签是为了减少需要计算的k点个数，减小计算量和WAVECAR大小。
第三步 利用前面第一步生成 CHGCAR; POSCAR POTCAR与第一步的一样，K点从第二步那个复制过来,然后提交任务得到波函数WAVECAR
第四步就要开始做能带反折叠了，一句命令，bandup unfold -emin -13 -emax 6 -dE 0.050 -no_symm_avg
会生成一个unfolded_EBS_not-symmetry_averaged.dat文件，画出来就行了。
https://github.com/band-unfolding/bandup
http://staff.ustc.edu.cn/~zqj/posts/Band-unfolding-tutorial/
可以处理能带反折叠的软件汇总：
https://mp.weixin.qq.com/s/-5lN2lwxsALA-rhNQxK1dw
https://vaspkit.com/tutorials.html#band-unfolding
https://github.com/QuantumMaterialsModelling/bands4vasp
https://github.com/QuantumMaterialsModelling/UnfoldingPatch4vasp
3D能带
__第一步__，准备好石墨烯的POSCAR和用于静态计算的INCAR文件；
__第二步__：运行VASPKIT并选择231生成用于3D能带计算 的KPOINTS文件。为了获得平滑的3D能带面，用于产生KPOINTS文件的分辨率值应该设置在0.005左右 ；
__第三步__，提交VASP作业；
__第四步__，待VASP计算完成后，再次运行VASPKIT，输入232或233命令。233命令可一次性输入包含价带顶的BAND.HOMO.grd和导带底的BAND.LUMO.grd文件的能带；232可以得到其它任意能带的计算数据，这里我们选择233。
__第五步__：运行python how_to_visual.py，绘制3D能带图(python2.7环境)。
如果3D能带图不够平滑，可通过能带插值进行改善。
有效质量：
除了通过能带，也可以通过调用VASPKIT-911命令来得到带边位置，接下来我们计算K点处沿着K -> Γ和K -> M方向的电子及空穴载流子的有效质量；
第一步：准备POSCAR文件以及VPKIT.in文件，VPKIT.in文件包含以下内容；
1 | 1 设置为1将会生成计算有效质量的KPOINTS文件，2则计算有效质量 |
第二步：运行vaspkit并选择912或者913产生KPOINTS，POTCAR文件及静态计算INCAR文件；
第三步：提交VASP作业；
第四步：计算完成后把VPKIT.in文件中第一行的1修改为2，然后再次运行vaspkit并选择913，得到结果：
为了得到可靠的有效质量计算值，括号里的精度至少要在1E-7这个量级，如果误差比较大，可以适当调整k-cutoff值，一般减少该值但不能太小。
手动vaspkit处理的能带：$m^∗=3.815/B$
GW方法
参考：https://www.vasp.at/wiki/index.php/Bandgap_of_Si_in_GW
步骤：
- PBE基态计算
- 获得虚拟轨道（50-100/atom）
INCAR：
1 | ALGO = Exact |
完成后：
1 | cp WAVECAR WAVECAR.DIAG |
- GW计算
读取上一步
INCAR：在OUTCAR中quasi-particle (QP) energies.1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20# Frequency dependent dielectric tensor including
# local field effects within the RPA (default) or
# including changes in the DFT xc-potential (LRPA=.FALSE.).
# N.B.: beware one first has to have done a
# calculation with ALGO=Exact, LOPTICS=.TRUE.
# and a reasonable number of virtual states (see above)
ALGO = GW0 ; LSPECTRAL = .TRUE. ; NOMEGA = 50
# be sure to take the same number of bands as for
# the LOPTICS=.TRUE. calculation, otherwise the
# WAVEDER file is not read correctly
NBANDS = 64
# Add this to update the quasiparticle energies
# in the Green's function (GW0)
#NELM = 4
ISMEAR = 0
SIGMA = 0.05
EDIFF = 1E-81
2
3
4
5
6
7QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 1
for sc-GW calculations column KS-energies equals QP-energies in previous step
and V_xc(KS)= KS-energies - (<T + V_ion + V_H > + <T+V_H+V_ion>^1 + <V_x>^1)
k-point 1 : 0.0000 0.0000 0.0000
band No. KS-energies QP-energies sigma(KS) V_xc(KS) V^pw_x(r,r') Z occupation
.........
LWANNIER90计算能带结构
重新开始第三步
INCAR：
1 | ALGO = GW0 ; LSPECTRAL = .TRUE. ; NOMEGA = 50 |
wannier90.win:
1 | num_wann=8 |
Compute Wannier functions
run wannier90:
wannier90.x wannier90
This run generates the wannier90 standard output (wannier90.wout) and the file wannier90.chk needed for the wannier interpolation (next step)
Obtain bandstructure (Wannier interpolation)
Uncomment the bandstructure plot flags in wannier90.win and rerun (restart) wannier90:
wannier90.x wannier90
This run generates the following bandstructure files which can be visualized using xmgrace or gnuplot:
wannier90_band.agr
wannier90_band.dat
wannier90_band.gnu
to plot the band structure using gnuplot: gnuplot -persist wannier90_band.gnu
Beyond the random-phase-approximation
To include local field effects beyond the random-phase-approximation in the description of the frequency dependent dielectric response function (local field effects in DFT) add the following line to your INCAR file:
LRPA = .FALSE.
and again restart from the WAVECAR and WAVEDER files from step 2.
Beyond G0W0: GW0
The most usual step beyond single-shot GW (G0W0) is to iterate the quasi-particle energies in the Greens functions. This is the socalled GW0 approximation. To have VASP do, for instance, 4 iterations of the QP-energies in G, add the following line to the INCAR file:
NELM = 4
and again restart from the WAVECAR and WAVEDER files from step 2.
To quickly find the QP-energy of the highest lying occupied state after 4 iterations of the QP energies in G, type:
./gap_GW.sh OUTCAR