VASP与Wannier90计算非线性位移电流
用GaAs实算展示波函数、Wannier插值、响应张量与数值对照的完整流程。
阅读本篇目录
非线性位移电流是光激发产生直流响应的一种机制，通常在缺少反演对称性的材料中研究。它需要带间矩阵元和相位信息，不能从普通能带图直接读出；VASP的LOPTICS线性光学计算也不能直接替代它。
本例采用鸢算科技集群实际计算的GaAs波函数，经过Wannier90得到频率相关位移电流张量。重点是让你知道先准备什么、怎样运行、输出怎样检查。
1. 示例设置#
两原子闪锌矿GaAs，晶格常数固定为5.65 Å，Ga位于原点、As位于原胞分数坐标(0.25,0.25,0.25)，空间群216、点群−43m。采用PBE PAW64 Ga/As，分别3和5个价电子；Ga的3d在芯内。本例不含SOC，也没有激子或剪刀修正。
VASP设置为ENCUT=550 eV、EDIFF=1E-6、6³Gamma网格、16条能带、固定原子与晶胞。该网格上的PBE带隙为0.5932 eV，不能把它当成实验光学带隙。
每个原子投影四个sp³轨道，共8个Wannier轨道。冻结窗口上限位于费米参照以上5 eV；程序同时检查冻结态数不超过8。与同一电荷密度下直接计算的高对称路径比较，费米参照−8至+3 eV范围内的能带误差为：RMSE 0.0189 eV，最大0.0640 eV。
VASP使用6.6.1-v100-wannier90（NVHPC 25.3/CUDA 12.8）；独立Wannier90/postw90使用4.0.3（Intel oneAPI 2024.2 ifx/AVX2）。后处理为CPU计算。
2. 工作流：从波函数到响应曲线#
在附件的两个DFT目录中准备合法POTCAR，顺序为Ga、As。回到附件根目录运行：
sbatch run.sh
脚本按顺序完成：
- gaas-sp3：VASP自洽，保存电荷密度，写出Wannier接口所需的投影、相邻k点重叠和本征值。
- 补齐Wannier输入，建立8轨道模型。
- gaas-band-check：使用自己的CHGCAR，在同一路径直接计算DFT能带，对照插值。
- 分别运行20³、40³积分网格，再改变中间态正则化参数及相位约定作比较。
- 提取18个独立张量分量，补齐后两个指标的交换对称，并按晶体点群作张量平均。
这一步的输入不能只剩HR文件；原始重叠矩阵等接口文件也要保留。附件中的脚本会在步骤之间传递本次产生的文件，不包含预先打包的波函数。
响应计算的核心设置：
berry = true
berry_task = sc
berry_kmesh = 40 40 40
kubo_freq_min = 0.0
kubo_freq_max = 3.0
kubo_freq_step = 0.05
kubo_adpt_smr = false
kubo_smr_fixed_en_width = 0.1
sc_eta = 0.04
sc_phase_conv = 1
这里的0.1 eV控制光谱展宽；sc_eta=0.04 eV用于处理中间态分母，二者用途不同。6³是生成DFT波函数的网格，40³是后续响应积分的网格，也不能混为一谈。
3. 得到哪些结果#
postw90输出wannier90-sc_xyz.dat等文件。第一列为光子能量，单位eV；第二列为张量分量，单位A/V²。表中为了便于阅读换成μA/V²。
| 积分网格 | sc_eta / eV | 相位约定 | 区间内最大绝对值所在能量 / eV | 对应σxyz / μA/V² |
|---|---|---|---|---|
| 20³ | 0.04 | 1 | 2.60 | -9.152 |
| 40³ | 0.04 | 1 | 3.00 | -8.264 |
| 40³ | 0.08 | 1 | 3.00 | -8.412 |
| 40³ | 0.04 | 2 | 3.00 | -8.264 |
表格和主响应曲线使用24个点群操作平均后的张量。GaAs在这一坐标约定下只有一类独立允许分量xyz及其指标排列；xxx应为零。张量平均使结果满足已知晶体对称性，但不会自动消除有限网格、模型和矩阵元误差。
因此图中同时保留原始xxx作为误差提示：40³组的最大绝对值约1.572 μA/V²，并不小，说明还需要继续改进和收敛模型。不能因为平均后的禁阻分量为零，就认为原始结果已经充分收敛。
20³与40³允许分量的最大差异约3.022 μA/V²，仍明显；把sc_eta从0.04改到0.08 eV时最大差异约0.251 μA/V²。两种相位约定给出的总响应一致。40³组表中的最大值落在3 eV扫描上限，不能据此认定谱峰位置。
4. 提取自己的数据#
全部任务结束后运行：
module load ase/3.29.0
python extract_shift.py
生成GAAS-SHIFT-DATA.json，包含原始完整张量、点群平均后的σxyz、能带插值对照和四组数值比较。参考数据在RESULTS.json，图在figures。
用于新材料时，先确认空间群与坐标方向，再选择投影、冻结窗口和响应能区；更换结构后不要沿用本例的GaAs点群结论。若材料带隙与实验差异较大，也要先明确采用哪种电子结构近似，再讨论光谱峰位。
5. 这些数值能说明什么#
本例给出一条实际可运行的位移电流工作流，尚不是GaAs的收敛材料定值。正式研究需要进一步增加DFT网格和响应积分网格、检查Wannier矩阵元与对称性、测试展宽，并考虑SOC、泛函和激子等因素。
σxyz是响应系数，还没有代入实际光场。它不等于器件短路电流、量子效率或光伏转换效率；光强、偏振、接触和载流子过程需要另外建模。原子朝向与坐标约定还会影响张量符号，比较不同资料时要先统一这些约定。
参考与下载#
Wannier90位移电流方法示例、postw90参数与输出定义。本文图表来自本集群VASP/Wannier90计算。
跟着输入文件操作#
下载本页附件后，先解压到单独目录。下面的目录名与附件一致。
unzip vasp-properties-batch30-*.zip -d vasp-example-30
cd vasp-example-30
cat README.txt
以下展开的是 calculations/gaas-sp3 这一步的真实示例输入。多结构扫描、位移法和多阶段响应仍按正文顺序准备全部任务；改变材料时，结构、元素顺序和方法参数要一起修改。
cd calculations/gaas-sp3
ls INCAR POSCAR KPOINTS
nano INCAR
INCAR：展开查看完整示例
SYSTEM = GaAs nonlinear response reference
ENCUT = 550
PREC = Accurate
EDIFF = 1E-6
IBRION = -1
NSW = 0
ISMEAR = 0
SIGMA = 0.05
LREAL = .FALSE.
LASPH = .TRUE.
ISYM = -1
NELM = 160
ALGO = Normal
NBANDS = 16
NUM_WANN = 8
LWANNIER90 = .TRUE.
LWAVE = .TRUE.
LCHARG = .TRUE.
POSCAR：展开查看完整示例
Ga As
1.0000000000000000
0.0000000000000000 2.8250000000000002 2.8250000000000002
2.8250000000000002 0.0000000000000000 2.8250000000000002
2.8250000000000002 2.8250000000000002 0.0000000000000000
Ga As
1 1
Direct
0.0000000000000000 0.0000000000000000 -0.0000000000000000
0.2500000000000000 0.2500000000000000 0.2500000000000000
KPOINTS：展开查看完整示例
Gamma mesh
0
Gamma
6 6 6
0 0 0
关键参数对照：
| 参数 | 本例值 | 使用时注意 |
|---|---|---|
IBRION |
-1 |
决定离子处理方式：−1为静态，0为分子动力学，2常用于优化；5/6等有限位移任务见正文。 |
NBANDS |
16 |
控制计算能带数；光谱、响应和Wannier窗口通常需要额外空带。 |
LWANNIER90 |
.TRUE. |
设为.TRUE.时输出Wannier接口；后续仍需独立运行Wannier90并检查插值。 |
LCHARG |
.TRUE. |
控制是否写出CHGCAR；后续固定密度计算需要前一步保存该文件。 |
EDIFF |
1E-6 |
电子自洽收敛阈值，单位eV；本组通常为1E−6。 |
POSCAR的元素顺序决定POTCAR拼接顺序。Direct表示分数坐标，Cartesian表示笛卡尔坐标；有Selective dynamics时，F/T还限制相应原子的移动。KPOINTS第二行为0通常使用自动网格；显式列表、Line-mode路径和均匀网格用途不同，不能互换来做同一种积分。
这个输入由附件根目录的总脚本组织。按正文为各阶段准备POTCAR后，返回附件根目录提交整个流程：
cd ../..
sbatch run.sh
squeue -u "$USER"
run.sh：整个流程的执行顺序
#!/bin/bash
#SBATCH --job-name=shift-current
#SBATCH --partition=gpu-8v100-16g-hs
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH --gres=gpu:v100-16g:1
#SBATCH --time=04:00:00
#SBATCH --output=shift-%j.log
set -eo pipefail
source /etc/profile
study_dir="$SLURM_SUBMIT_DIR"
export OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 OMPI_MCA_coll_hcoll_enable=0 NO_STOP_MESSAGE=1
ulimit -s unlimited
module purge
module load vasp/6.6.1-v100-wannier90
cd "$study_dir/calculations/gaas-sp3"
mpirun --bind-to none -np 1 vasp_std > stdout.log 2>&1
module purge
module load ase/3.29.0
python "$study_dir/prepare_wannier.py"
module purge
module load wannier90/4.0.3-serial
wannier90.x wannier90 > wannier.stdout 2>&1
module purge
module load ase/3.29.0
python "$study_dir/prepare_bandcheck.py"
cd "$study_dir/calculations/gaas-band-check"
module purge
module load vasp/6.6.1-v100
mpirun --bind-to none -np 1 vasp_std > stdout.log 2>&1
module purge
module load ase/3.29.0
python "$study_dir/check_bands.py"
python "$study_dir/prepare_response.py"
module purge
module load wannier90/4.0.3-serial
for name in mesh20 mesh40 eta08 phase2; do
cd "$study_dir/calculations/$name"
postw90.x wannier90 > post.stdout 2>&1
grep -q 'All done: postw90 exiting' wannier90.wpout
done
作业结束后，在实际VASP输出目录查看最后几步：
tail -n 12 OSZICAR
grep "EDIFF is reached" OUTCAR | tail -n 3
grep "reached required accuracy" OUTCAR
ls -lh CONTCAR
静态计算不要求出现离子优化完成提示；优化任务还要检查最终自由原子力。性质专用输出和后处理仍按正文执行。附件中的参考图与RESULTS文件不会随着你的新任务自动更新。
需要选择后处理软件时，可对照实算教程中的辅助工具表。