自旋轨道耦合计算探索过程

自旋轨道耦合计算过程探索
1.经验总结
1）对于Bi2Se3家族材料，QL内是强的共价结合作用，QL之间是范德瓦尔斯作用力。

所以，在优化结构的时候，需要考虑范德瓦尔斯相互作用。

一般，对于一种没有算过的新材料，可以尝试以上五种方法，哪一种最合理就用哪个。

Bi2Se3家族材料，经测试最合适的是optPBE-vdW方法。

3）测试发现，对于1QL和块体，范德瓦尔斯作用的影响不是很影响；对于多个QL厚度的薄膜，QL之间范德瓦尔斯作用的影响比较明显。

4）文献上，很多人直接不优化结构，用实验上的参数，这样算，得到的结果也比较合理。

5）算soc加入LSORBIT=.TRUE.和LORBMOM=.TRUE.，
比LSORBIT=.TRUE.和GGA_COMPAT = .FALSE.得到的结果更合理。

6）薄膜优化的时候，可以用ISIF=2。

7）计算静态的时候输出CHARG，能带的时候ISTART可以等于0，ICHARG等于11。

7）薄膜的结构需要中心对称，切得时候需要注意。

8）计算vdW，需要vasp5.2.12以上的版本，并且将vdw_kernel.bindat文件放到计算的文件夹中。

9）vdW相互作用对结构的影响比较大，对后面的静态计算和能带计算电子态的影响比较小。

10）取合适的K点，可以得到较为合理的结构，对后面电子态的计算影响也不是很大。

2. 结构优化
赝势：PAW_GGA_PBE E cut=340 eV Kpoints=10×10×10
ISMER取-5，计算能带时，取0，对应SIGMA=0.05
在MS中可以在build-Symmetry -中把Bi2Se3 rhombohedral representation（菱形表示）和hexagonal representation（六角表示）相互转换
图中黑色t1、t2、t3基矢围成菱形原胞，用于计算块体，红色方框包含一个五元层计算能带的布里渊区高对称点：
块体:文献中倒空间高对称点坐标Г(0 0 0)-Z(π π π)-F(π π 0)-Г(0 0 0)-L(π 0 0)，
根据正空间和倒空间坐标的转换关系，得到正空间中高对称点的坐标：
Г(0 0 0)-Z(0.5 0.5 0.5)-F(0.5 0.5 0)-Г(0 0 0)-L(0 0 -0.5)
KPOINTS
20
Line-mode
Rec
0.0 0.0 0.0 !Г
0.5 0.5 0.5 ! Z
0.5 0.5 0.5 ! Z
0.5 0.5 0.0 ! F
0.5 0.5 0.0 ! F
0.0 0.0 0.0 !Г
0.0 0.0 0.0 !Г
0.0 0.0 -0.5 ! L
[通过比较结构，发现Ecut=580，KPOINTS=151515，得到的结构比较靠谱]
3.块体soc的计算
文献能带结构图：
块体（Bi2Se3-VASP-GGA-PAW-PBE）Sb2Te3Bi2Te3Bi2Se3
晶格参数六角a (Å) 4.250 4.383 4.138 六角c (Å) 30.35 30.487 28.64 菱形T(Å) 10.41 10.473 9.841
原子位置μ (2Bi)0.400 0.400 0.399 v (2Se2) 0.211 0.212 0.206 0 (Se1) 0 0 0
我们的结果（未考虑vdW+静态和能带都加soc计算结果与文献基本符合）：
4.薄膜的计算
薄膜：Kpoints=10×10×1
计算能带的K点和石墨烯（六角晶胞的）的K点一样:
KPOINTS
20
Lone-mode
Rec
0.66666667 0.33333333 0.0 !K
0.0 0.0 0.0 !Г
0.0 0.0 0.0 !Г
0.5 0.0 0.0 !M
考虑薄膜的对称性
由MS六角结构，沿（001）方向切割，可以得到两种以Se原子作为表面原子的薄膜，如下图，分别为1QL和3QL的两种切法，右图比左图对称性要更好一些，这一区别在计算过程中会导致巨大的区别，我们通过比较，发现，只有右图的结果，才可以得到合理的结果，尤其是在多个QL的情况。

用左边结构得到的结果（Bi2Se3）：
用左边结构得到的结果（Bi2Te3）：
用右边结构得到的结果（Bi2Se3）：
1QL—根据块体的数据得到薄膜，分以下两种情况计算：
1.不优化结构，scf不加soc，bands加soc
2.用块体的参数，加vdW优化结构，scf不加soc，bands加soc
1QL—在静态中也加入soc
1.不优化结构，scf和bands加入LORBMOM=.TRUE.，LSORBIT=.TRUE.
2.不优化结构，scf和bands加入LORBMOM=.TRUE.，GGA_COMPA T=.FALSE.
3.优化结构，scf和bands加入LORBMOM=.TRUE.，LSORBIT=.TRUE.
4.优化结构，scf和bands加入LORBMOM=.TRUE.，GGA_COMPA T=.FALSE.
文献结果：
上图是没有进行离子弛豫的1QL ~6QL的Bi2Se3薄膜能带结构
上图采用optPBE-vdW 泛函进行离子弛豫1QL~6QL的Bi2Se3薄膜能带结构
上图是实验观测的1QL ~6QL（12356）的Bi2Se3薄膜能带结构.
5.调试过程错误总结
错误1：VERY BAD NEWS! Internal error in subroutine IBZKPT: Reciprocal lattice and k-lattice belong todifferent class of lattices. Often results are still useful (48)
Internal内部 subroutine子程序 Reciprocal倒数的
非常严重的错误！子程序IBZKPT中内部错误：倒格子和k点网格属于不同类型的格子。

通常结果还是有用的。

解决方案：根据所用集群，修改INCAR中NPAR。

将NPAR=4变成NPAR=1，已解决！
错误2：
internal ERRORRSPHER:running out of buffer 0 0
13 1 0
nonlr.F:Out of bufferRSPHER
解决方案：根据CPU的数量，修改INCAR中NPAR，将NPAR=1修改成4（或者2），问题得以解决。

错误3：WARNING: Sub-Space-Matrix is not hermitian in DAV 4 -4.681828688433112E-002 Sub-Space-Matrix 子空间矩阵、亚空间矩阵 Hermitian厄米共轭
警告：戴维森方法（DAV）中的子空间矩阵不是厄米共轭的。

解决方案：只需调整 AMIX, BMIX的值，把他们设置小一些。

一般采用其默认值，除非在电子迭代难以收敛的情况，才手动设置AMIX和BMIX等参数值。

经对Mixing方法的调试，通过将默认AMIX=0.4，修改成AMIX=0.2（或0.3），问题得以解决。

Mixing方法：
IMIX=type ofmixing混合、混频， AMIX=linear mixing parameter， AMIN=minimal mixing parameter,
BMIX=cutoffwave vector for Kerker mixing scheme, AMIX_MAG=linear mixing parameter for magnetization,
BMIX_MAG=cutoffwave vector for Kerker mixing scheme for mag, WC=weight factor for each step inBroyden mixing scheme,
INIMIX=type ofinitial for each step in Broyden mixing scheme, MIXPRE=type of preconditioning in Broyden mixing scheme,
MAXMIX=maximumnumber steps stored in Broyden mixer.
错误4：WARNING: Sub-Space-Matrix is not hermitian in DAV 1 -7.626640664998020E-003 解决方案：在INCAR中加上IALG=Fast 已解决！（1QL、2QL已解决，3QL以上未解决）IALG=Fast (两种方法混用)
IALGO=38 IALG=Normal 电子优化采用blockedDavidson方法
IALGO=48 IALG=Very_Fast 电子优化采用RMM-DIIS算法
错误5：ADVICE TO THIS USER RUNNING 'VASP/VAMP' (HEAR YOUR MASTER'S
VOICE ...): You have a(more or less) 'small supercell' and for smaller cells it is recommended to usethe reciprocal-space projection scheme! The real space optimization is notefficient for small cells and it is also less accurate ... Therefore setLREAL=.FALSE. in the INCAR file
解决方案：对于较小的晶胞（原子数小于20），设置LREAL=.FALSE.，计算结果比较精确。

而对于较大的晶胞，设置LREAL=Auto，这样计算速度比较快。

对于1QL 2QL 3QL原子数分别为5、10、15，LREAL=.False.
对于4QL 5QL 6QL原子数分别为20、25、30，LREAL=Auto
错误6：自旋轨道耦合计算时，静态和能带计算中出现的错误：
ERROR: non collinear calculations require that VASP is compiled withoutthe flag -DNGXhalf and
-DNGZhalf
错误：非线性计算需要编译过的VASP，VASP中不包含-DNGXhalf 和-DNGZhalf
解决方案：重新编译VASP。

don't forget that you may have to re-compile vasp without any of theprecompiler (CPP) flags set:
-DNGXhalf, -DNGZhalf, -DwNGXhalf, -DwNGZhalf , asnecessary for non-collinear runs in general for non-collinear magnetism
不要忘记如果你用的vasp不包含任何预编译程序命令 -DNGXhalf, -DNGZhalf, -DwNGXhalf,
-DwNGZhalf ，你必须重新编译vasp只有编译过，因为这些参数通常对于非线性磁性计算是必要的.
错误7: SB-3QL计算总是出现这个错误 VERY BAD NEWS! internal error in subroutine SGRCON: Found some non-integer element in rotation matrix 3
解决办法：首先检查POSCAR是否有问题。

MS重新获取POSCAR，参与计算，已解决。

错误8: WARNING: random wavefunctions but no delay for mixing, default for NELMDL
建议解决办法：一般不是很影响，可以继续算。

可以设置NELMDL=-5或者其它的数,是前5步保持电荷密度不变,只弛豫波函数。

错误9: 计算结构得到一个空的，然后算静态也没有数据，能带也没有数据，但都算完了。

输出文件中提示：
mpirun has exited due to process rank 21 with PID 8455 on node node14 exiting improperly. There are two reasons this could occur:1. this process did not call "init" before exiting, but others in the job did. This can cause a job to hang indefinitely while it waits for all processes to call "init". By rule, if one process calls "init", then ALL processes must call "init" prior to termination.2. this process called "init", but exited without calling "finalize". By rule, all processes that call "init" MUST call "finalize" prior to exiting or it will be considered an "abnormal termination"This may have caused other processes in the application to be terminated by signals sent by mpirun (as reported here).
解决方案：这是VASP版本的问题，考虑vdW作用，有的版本的vasp算不了，一般5.2.12以上的版本可以，换一个版本就可以计算。

错误10：计算总是卡住，输出的log文件中只有前面的描述，没有后面的数据，比如：
running on 32 total cores[2016-05-19 18:42:45] distrk: each k-point on 32 cores, 1 groups[2016-05-19 18:42:45] distr: one band on 8 cores, 4 groups[2016-05-19 18:42:45] using from now: INCAR[2016-05-19 18:42:45] vasp.5.3.3 18Dez12 (build May 19 2015 15:36:57) complex[2016-05-19 18:42:45] [2016-05-19 18:42:45] POSCAR found type information on POSCAR Se Bi[2016-05-19 18:42:45] POSCAR found : 2 types and 15 ions[2016-05-19 18:42:45] LDA part: xc-table for Pade appr. of Perdew[2016-05-19 18:42:45] WARNING: stress and forces are not correct[2016-05-19 18:42:46] POSCAR, INCAR and KPOINTS ok, starting setup[2016-05-19 18:42:46] FFT: planning ...[2016-05-19 18:42:47] WAVECAR not read[2016-05-19 18:42:48] charge-density read from file: Bi2Se3scf[2016-05-19 18:42:49] magnetization density of overlapping atoms calculated
解决方案：这种情况一般是由于计算量过大，算不动导致的，通过减少精度，就可以解决。

比如，减少K点数量（主要的影响因素）。

错误11：VERY BAD NEWS! internal error in subroutine SGRCON:
Found some non-integer element in rotation matrix 2
解决问题：一般是结构的问题。

重新得到POSCAR，可以解决。