能带反折叠计算

最近在研究一些关于电荷密度波CDW的文献，需要研究CDW相变之前和CDW相变之后的能带，由于发生CDW相变都会伴随着晶格周期的加倍，这时候布里渊区会相应的减小，能带会产生折叠，如果直接看CDW相的能带以及未发生晶格畸变的能带，会发现差距比较大，难以进行比较，所以需要把能带反折叠。
关于能带反折叠的教程网上有很多，但是小编感觉都不太详细，于是自己重新整理了以下，具体步骤如下。

第一步

准备一个正常相的POSCAR，也就是没有扩胞的结构。

第二步

准备KPATH.in，具体内容如下：

K-Path Generated by VASPKIT.
   20
Line-Mode
Reciprocal
   0.5000000000   0.0000000000   0.0000000000     M    
   0.0000000000   0.0000000000   0.0000000000    #GAMMA          
          
   0.0000000000   0.0000000000   0.0000000000    #GAMMA             
   0.3333333333   0.3333333333   0.0000000000     K              
 
   0.3333333333   0.3333333333   0.0000000000     K              
   0.5000000000   0.0000000000   0.0000000000     M  

第三步

产生超胞的变化矩阵TRANSMAT.in,这一个矩阵可以自己手动改，嫌麻烦的话也可以用vaspkit自动生成。

Read transformation matrix from the TRANSMAT.in file if it exists. 
    2    0    0          # must be three integers
    0    2    0          # must be three integers
    0    0    1          # must be three integers

上面这个就是我的TRANSMAT.in文件，因为我的CDW相相比于普通相，相当于扩胞2*2*1，是原来基矢的整数倍，所以对角线上都是0，如果不是整数倍，那就不是一个对角矩阵，具体的矩阵需要自己求出来。
当然，这个矩阵也可以用vaspkit生成，用vaspkit的401命令产生，这就是为什么第一步我们需要原胞的POSCAR的原因。

第四步

准备VASP计算的输入文件，首先是POSCAR，注意，这个POSCAR是CDW相，也就是扩胞后发生晶格畸变的POSCAR！！！也就是说，第一步的POSCAR已经没用了，可以删了，第一步的主要作用就是为了能让你用vaspkit方便的产生TRANSMAT.in文件，如果你这个文件自己手动输入，第一步就可以省略！！！
接着准备POTCAR，输入vaspkit-281，可以利用KPATH.in产生KPOINTS，最后准备自洽计算的INCAR，我的INCAR设置如下，最后提交任务计算即可。

!control
ISTART=0
ICHARG=2
ENCUT = 450
NELM = 200
ISMEAR=0
SIGMA=0.01
ISPIN=1
EDIFF=1E-7
PREC=A
ALGO = N
GGA = PE

#####Ion steps:
IBRION=-1
NSW=0
ISIF=3
EDIFFG=-0.001

#####DOS and bands: 
LORBIT = 11

!speed
NPAR=4
LPLANE = .TRUE.
LSCALAPACK=F
LSCALU=F
NSIM   = 4

!ISYM=2
!SYMPREC=1E-7

!magnetic
!MAGMOM = 6*1

LREAL = F


!+U
LDAU = T
LDAUTYPE = 2
LMAXMIX = 4
LDAUL=  -1  2  
LDAUU= 0.0  3.95  
LDAUJ= 0.0  0.0  
LASPH=T
#MAGMOM=0 0 0 0 1 -1

第五步

计算完成之后，我们就可以进行数据处理了，输入vaspkit-282产生EBS.dat，最后利用vaspkit自带的python小程序画图，输入：python vaspkit/examples/band_unfolding/ebs_plot.py，就会产生能带的pdf文件，如果这个小程序没改，它的能窗可能比较大，如下图所示：

这时我们可以修改代码，将ax.set_ylim(-10,10)改成ax.set_ylim(-0.5,0.5)，最终效果如下：

对比下图正常的未发生畸变的结构，可以明显看出哪里的能带发生变化：