VASP+Wannier90 计算非线性光响应

（记录一下以备后用）

帮组里同学用Wannier90算二次谐波SHG。自己之前都是QE+Wannier90，但因为他用VASP的直接给我了VASP的数据，不得不搞VASP+Wannier90. 按理说VASP应该更好用更方便的。

1. 编译带Wannier90接口的VASP，此处略过（组里集群上有现成）。另安装可计算二次谐波SHG和位移电流Shift current的wannier90，版本在https://github.com/peio-gg/wannier90/tree/eps_shg_mb ，编译无坑可放心食用。

2. 做VASP的基态计算，此步同学直接给我数据。

3. 用VASP输出Wannier90的输入文件（接口）。此时按照官网教程似乎只用在INCAR加入LWANNIER90=.true.，这样跑是能跑，但不会输出wannier90.amn文件，从而导致后面运行wannier90.x报错。原因是以lib模式调用wannier90时没有projection信息。所以如果想快速demo简单体系，最极速懒人的做法是在INCAR里加入如下内容：

   LWANNIER90 = .TRUE.
   LWRITE_MMN_AMN = .TRUE.
   WANNIER_SETUP = .FALSE.
   WANNIER90_WIN = "Begin Projections
 Random
End Projections

guiding_centres=.true.
use_bloch_phases=.false."

这个WANNIER90_WIN参数里的东西就是最后会写到wannier90输入文件的额外参数。另注意设置NCORE=1或不额外设置NPAR等其他并行参数否则会报错。因此可以先开并行做VASP基态自洽计算（ICHARG=2，LWAVE=T，LCHARG=T），保存CHGCAR后再做接口计算时用非自洽（ICHARG=11，LWAVE=F，LCHARG=F）可以稍微快点。

运行VASP后就能得到粗略跑一个wannier90所需的全部输出，都是以wannier90为seedname，如输入文件wannier90.win，轨道的overlap文件wannier90.amn、wannier90.mmn等。

4. 用wannier90做局域化。运行mpirun -n 6 wannier90.x wannier90得到wannier90.chk，查看wannier.wout里的wannier轨道展宽spreading应小于晶胞参数，一般在3埃以内。如果是自旋极化的材料但是不开SOC，那VASP做wannier后会自动输出自旋上、下的两套数据，seedname分别为wannier90.1和wannier90.2，此时跑wannier90.x需要分别运行mpirun -n 6 wannier90.1.x wannier90和mpirun -n 6 wannier90.2.x wannier90，相应的查看wannier.1.wout, wannier.2.wout 。此时如果是稍微复杂点的体系，3中所用的懒人Random projection就不太够用了会导致展宽太大，此时需要把projection换成原子轨道并添加更大的局域化迭代次数num_iter，如保险起见再调大解纠缠迭代次数dis_num_iter，比如Cu2V2P6S12需要

Begin Projections
 Cu: s;p;d
 V: s;p;d
 P: s;p
 S: s;p
End Projections

dis_num_iter=100
num_iter=10000

我测试的56个核跑100个轨道的wannier，在该参数下3-5分钟跑完还是很快的。这样拟合的wannier函数才能很好局域至可用水平（平均1.6埃），实测相比Random，该材料的spreading减小了2/3. 另一个可选步骤是画Wannier能带对比VASP能带来验证瓦的好不好。此时在wannier90.win里加入

Begin kpoint_path
 G 0.0 0.0 0.0  X 0.5 0.0 0.0
 X 0.5 0.0 0.0  L 0.5 0.5 0.0
 L 0.5 0.5 0.0  G 0.0 0.0 0.0
End kpoint_path
bands_num_points = 30
bands_plot=.True."

然后再运行wannier90.x即可得到_band结尾的能带数据。

5. 后处理计算非线性光学性质。在前一步跑mpirun -n 6 wannier90.x wannier90得到wannier90.chk的基础上， 在wannier90.win中添加

berry = true
berry_task=sc #sc：位移电流shift current，shg：二次谐波second harmonic generation
berry_kmesh = 15 15 1 #需要测试收敛性，一般是scf密度的数倍
fermi_energy=-3.65 #从OUTCAR读取

kubo_freq_min= 0.0 #光频
kubo_freq_max= 2.0
kubo_freq_step= 0.02 
sc_phase_conv=2
kubo_adpt_smr=false
smr_fixed_en_width=0.025 #delta函数展宽

然后运行mpirun -n 6 postw90.x wannier90得到光响应，sc每个文件2列，分别是光频、sc conductivity，单位分别是eV、A/V^2；shg每个文件5列，分别是光频、real shg susceptibility、imag shg susceptibility、real shg conductivity、imag shg conductivity，前3列单位分别为eV、pm/V。

关于对称性：wannier90对二阶光计算的点群对称性保持的并不特别理想（虽然据说已经比其他同行好很多了；主要来源于二阶对数值噪声的敏感），比如晶格对称性要求光响应某个分量\(\chi_{abc}, (a,b,c=x,y,z)\) 为零，但实际由于数值误差算得不为零且量级并不小。为此可以手动计算对称性变换后的结构的光响应，然后求平均来消除误差。例如我算的C2对称性的材料C2轴在y，直接计算\(\chi_{xxx}≠0\)且量级与其他相比并无很大差别，此时把POSCAR里的每个原子的x、z坐标加负号做C2轴的变换，得到新的POSCAR，然后再做计算得到\(\chi_{abc}'\)，然后\(\chi''_{abc} = [\chi_{abc} + (-1)^{i_a+i_b+i_c}\chi'_{abc}]/2\)其中\(i_x,i_y,i_z=1,2,3\)，最后得到的光响应\(\chi''_{xxx}\)几乎为0符合对称性要求。更粗暴一点可以不必重新计算，直接按照对称性做组合。

算wannier90光响应对INCAR的总修改汇总：

LWANNIER90 = .TRUE.
LWRITE_MMN_AMN=.TRUE.
#NUM_WANN=112
wannier_setup=.FALSE.
WANNIER90_WIN="
# BERRY
###############
berry = true
berry_task=sc
berry_kmesh = 15 15 1
fermi_energy=-3.65

kubo_freq_min= 0.0
kubo_freq_max= 2.0
kubo_freq_step= 0.02
#kubo_adpt_smr=true
#sc_eta=0.040
sc_phase_conv=2
kubo_adpt_smr=false
smr_fixed_en_width=0.025
###############

#Begin Projections
# Random
#End Projections

Begin Projections
 Cu:s;p;d
 V:s;p;d
 P:s;p
 S:s;p
End Projections

dis_num_iter=100
num_iter=10000

guiding_centres=true

Begin kpoint_path
 G 0.0 0.0 0.0  X 0.5 0.0 0.0
 X 0.5 0.0 0.0  L 0.5 0.5 0.0
 L 0.5 0.5 0.0  G 0.0 0.0 0.0
End kpoint_path
bands_num_points = 30
bands_plot=.True."