py4vasp声子投影与单位

本篇接续声子色散与DOS实操。已有VASP输出后，py4vasp负责读取和作图，不会替您补算缺失的声子路径或DOS。

下载实际数值、完整输入与脚本。使用 bash run.sh 重画参考图；使用 bash extract.sh 读取您完成的 work/ 计算并输出到 analysis-my-run/。

1 读取真正的声子路径#

在资料包根目录执行：

module load py4vasp/0.11.3
python - <<'PYCODE'
import py4vasp
c = py4vasp.Calculation.from_path('work/forceconstants')
b = c.phonon.band.read()
print(b.keys())
print(b['bands'].shape)
print(b['modes'].shape)
print(b['bands'][0])
PYCODE

本例实际形状为：

bands: (80, 6)
modes: (80, 6, 2, 3)

分别对应80个q点、6条原胞分支；本征矢另外有2个原胞原子和3个笛卡尔方向。模式数组可为复数，处理时保留相位，不要直接丢掉虚部。

本版本这些频率字段的单位是 THz，不是电子能带常用的eV。原生OUTCAR同时打印THz、2πTHz、cm⁻¹、meV，比较时必须选同一单位。这里480个频率与原生THz输出的最大差约6.2×10⁻⁷ THz。

这份HDF5的路径标签字段为空，资料包绘图根据实际QPOINTS的段边界标Γ、X、W、Γ、L；不是软件凭坐标自动识别出的标签。用户修改路径后，也要同步修改标记。

2 导出表格与交互图#

python - <<'PYCODE'
import py4vasp
c = py4vasp.Calculation.from_path('work/forceconstants')
c.phonon.band.to_frame().to_csv('phonon-band.csv', index=False)
fig = c.phonon.band.to_graph('Si(x,y,z)').to_plotly()
fig.write_json('phonon-band.json')
PYCODE

selections() 可以列出这份数据实际支持的原子和方向选择：

python - <<'PYCODE'
import py4vasp
c = py4vasp.Calculation.from_path('work/forceconstants')
print(c.phonon.band.selections())
PYCODE

本例原子选择是 Si、1、2，方向是 x、y、z。投影色散图的线宽是软件的可视化权重，不宜把不同图形接口的权重都直接称为概率。

3 读取方向和原子投影DOS#

python - <<'PYCODE'
import py4vasp
c = py4vasp.Calculation.from_path('work/dos-mesh12')
d = c.phonon.dos.read('Si(x,y,z)')
a = c.phonon.dos.read('1,2')
print(d.keys())
print(a.keys())
c.phonon.dos.to_frame('Si(x,y,z)').to_csv('phonon-dos.csv', index=False)
c.phonon.dos.to_plotly('Si(x,y,z)').write_json('phonon-dos.json')
PYCODE

这里两个原胞Si位置对称等价，两条原子投影曲线重合；总DOS等于两原子投影之和，也等于x、y、z三个方向贡献之和。本例逐点求和误差小于10⁻¹⁵，不需要人为把投影再归一化才能凑出总谱。

4 用积分检查是否读对了模式数#

python - <<'PYCODE'
import numpy as np
import py4vasp
c = py4vasp.Calculation.from_path('work/dos-mesh12')
d = c.phonon.dos.read('Si(x,y,z)')
print(np.trapezoid(d['total'], d['energies']))
print(np.trapezoid(d['Si_x'], d['energies']))
PYCODE

实际结果：

两个原胞原子对应6个模式。积分一致说明谱的模式数和单位读取合理，并不证明网格或超胞已经收敛；两档谱形仍应直接比较。

5 展宽和负频区怎么看#

VASP的 PHON_SIGMA 使用eV；本例0.001 eV约为0.241799 THz。py4vasp的频率横轴是THz，相应DOS是每THz的模式密度。不能把eV展宽数值不换单位就标在THz图上。官方PHON_SIGMA说明明确了输入单位。

高斯展宽在谱端外延伸，本例频率横轴因此从约−1.209 THz开始。这个负频横轴范围本身不是虚频证据。是否存在不稳定模，应查看未展宽的声子分支或模式频率，并区分数值接近零的声学模。

最后保留HDF5和OUTCAR原件、CSV及输入；若只修改画图样式，重新运行绘图即可，不要重算力常数。