py4vasp声子色散与DOS实操
Si小超胞真实力常数、完整声子路径和均匀q网格DOS，展示VASP计算、复用力常数及py4vasp读取绘图。
阅读本篇目录
已有Γ点模式不等于声子色散。要画随q点变化的分支，先从超胞得到力常数，再做傅里叶插值；声子DOS使用均匀q网格，不能用一条高对称路径代替。
本篇在鸢算科技集群实际计算Si16原子超胞，得到80个路径点、每点6条原胞声子分支，并用同一力常数计算6³和12³均匀网格DOS。下载完整输入、数值和绘图程序，保存为 py4vasp-phonon.zip。
1 先看结构与结果#
采用晶格常数5.43 Å的Si金刚石原胞，重复2×2×2得到16原子超胞。这里固定晶胞，使用对称参考结构计算有限位移力；参考结构的最大原子力为0。它是有限超胞教学例，没有完成超胞尺寸、k网格和截断能的材料级收敛研究。
VASP 6.6.1 HDF5 GPU构建使用NVHPC 25.3、CUDA 12.8和串行HDF5 1.14.6；P100/V100 16G使用对应模块和配套MPI。py4vasp固定0.11.3，绘图用ASE 3.29.0环境中的Matplotlib。
先重画已有数值，不提交新计算：
python3 -m zipfile -e py4vasp-phonon.zip .
cd py4vasp-phonon-band-dos-wave282
bash run.sh
ls analysis/figures
2 准备力常数计算#
bash prepare.sh
选择 1（力常数与色散），GPU选择 1 为P100、2 为V100 16G，再按提示提交：
cd work/forceconstants
sbatch run.sh
程序会复制完整输入和本地授权Si势文件。实际任务使用1个MPI进程与1张GPU，核心命令是：
module load vasp/6.6.1-p100-hdf5
mpirun --bind-to none -np 1 vasp_std > vasp.log 2>&1
完整POSCAR：
Si
1.0000000000000000
0.0000000000000000 5.4299999999999997 5.4299999999999997
5.4299999999999997 0.0000000000000000 5.4299999999999997
5.4299999999999997 5.4299999999999997 0.0000000000000000
Si
16
Direct
0.0000000000000000 0.0000000000000000 -0.0000000000000000
0.1250000000000000 0.1250000000000000 0.1250000000000000
0.0000000000000000 0.0000000000000000 0.5000000000000000
0.1250000000000000 0.1250000000000000 0.6250000000000000
0.0000000000000000 0.5000000000000000 -0.0000000000000000
0.1250000000000000 0.6250000000000000 0.1250000000000000
0.0000000000000000 0.5000000000000000 0.5000000000000000
0.1250000000000000 0.6250000000000000 0.6250000000000000
0.5000000000000000 0.0000000000000000 -0.0000000000000000
0.6250000000000000 0.1250000000000000 0.1250000000000000
0.5000000000000000 0.0000000000000000 0.5000000000000000
0.6250000000000000 0.1250000000000000 0.6250000000000000
0.5000000000000000 0.5000000000000000 -0.0000000000000000
0.6250000000000000 0.6250000000000000 0.1250000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.6250000000000000 0.6250000000000000 0.6250000000000000
完整INCAR：
SYSTEM=Si16-native-phonon-band
ENCUT=450
PREC=Accurate
EDIFF=1E-6
EDIFFG=-0.02
IBRION=6
NFREE=2
POTIM=0.015
NSW=1
ISIF=2
ISYM=2
NBANDS=48
KPAR=1
ISMEAR=0
SIGMA=0.05
ALGO=Normal
NELM=160
LREAL=.FALSE.
LASPH=.TRUE.
LWAVE=.FALSE.
LCHARG=.FALSE.
LPHON_DISPERSION=.TRUE.
PHON_NWRITE=1
IBRION=6 通过有限差分计算力常数，并利用对称性减少独立位移；NFREE=2 使用正负位移，POTIM=0.015 是位移幅度，单位Å。本例实际有1个独立位移自由度，原始参考加正负位移共3次力评估。
EDIFF=1E-6 控制电子自洽。计算声子前先确认参考结构的原子力足够小，通常可按0.02 eV/Å的结构优化阈值准备；不能要求每个故意位移的结构也满足零力。
电子KPOINTS：
Si16 electronic mesh
0
Gamma
3 3 3
0 0 0
声子路径单独写在QPOINTS中，电子k点和声子q点不要混用：
Si primitive phonon reference path
20
Line-mode
Reciprocal
0 0 0 ! G
0.5 0 0.5 ! X
0.5 0 0.5 ! X
0.5 0.25 0.75 ! W
0.5 0.25 0.75 ! W
0 0 0 ! G
0 0 0 ! G
0.5 0.5 0.5 ! L
每段20点，四段共80点，段边界有重复端点。LPHON_DISPERSION=.TRUE. 会利用刚得到的力常数计算这条原胞路径；PHON_NWRITE=1 把各q点频率写入OUTCAR。
3 检查输出并保留力常数#
tail vasp.log
ls -lh vaspout.h5
grep 'q-point No.' OUTCAR | head
保留 vaspout.h5、OUTCAR和输入。后面的DOS可以直接复用力常数，不需要再做有限位移与SCF。
4 生成均匀网格DOS#
返回资料包根目录：
cd ../..
bash prepare.sh
选择 2 准备6³网格DOS，GPU类型按实际资源选择，再提交：
cd work/dos-mesh6
sbatch run.sh
返回根目录选择 3，准备并提交12³网格DOS。脚本把 work/forceconstants/vaspout.h5 复制为各目录的 vaspin.h5，并使用以下INCAR：
SYSTEM=Si16-reused-FC-DOS6
LPHON_READ_FORCE_CONSTANTS=.TRUE.
LPHON_DISPERSION=.TRUE.
PHON_DOS=1
PHON_SIGMA=0.001
PHON_NEDOS=301
PHON_NWRITE=1
ENCUT=450
EDIFF=1E-6
IBRION=-1
NSW=0
ISYM=2
ISMEAR=0
SIGMA=0.05
LREAL=.FALSE.
LWAVE=.FALSE.
LCHARG=.FALSE.
DOS的QPOINTS：
Uniform phonon DOS mesh
0
Gamma
6 6 6
0 0 0
LPHON_READ_FORCE_CONSTANTS 读取已有力常数，VASP完成声子后处理后直接退出；本例两档DOS实际没有新的SCF迭代。PHON_DOS=1 是高斯展宽，PHON_SIGMA=0.001 的单位是eV，约等于0.241799 THz。
DOS必须用均匀网格。即使程序接受路径点并给出谱形，也不能把路径采样当作代表整个布里渊区的态密度。VASP官方流程对两类QPOINTS及力常数复用有相同区分。
5 提取并画自己的结果#
两档DOS都结束后，回到资料包根目录：
bash extract.sh
ls analysis-my-run/tables
ls analysis-my-run/figures
脚本调用真实py4vasp接口，导出CSV与Plotly数据，再绘制静态图；没有重新调用VASP。它还会把80×6个频率与原生OUTCAR逐项比较，并检查DOS的模式数与投影求和。
Γ点三个光学模式均约15.429901 THz，三个声学模式的数值在约10⁻⁷ THz量级，属于这里的数值零。原胞有两个原子，因此每个q点应有6条分支；不要把16原子超胞的48个位移自由度误认为插值后的原胞分支数。
接下来见声子投影、单位与积分。换成其他材料时，需重新建超胞、生成相应声子路径，并检查极性材料的长程修正与超胞收敛；不能只替换元素名后照抄Si的结果。