FourPhonon：与TDEP比较声子频率及Γ点
用原生路径和27点网格说明频率单位、q点匹配及精确Γ点与方向性LO–TO极限的区别。
阅读本篇目录
下载ZnS极性TDEP与FourPhonon示例包，保存为example.zip：
python3 -m zipfile -e example.zip .
cd polar-tdep-fourphonon-examples
bash results.sh
bash read.sh
bash plot.sh
以上命令读取本例的实算结果并画图，不会重算VASP。需要自己重新拟合时，运行sbatch run.slurm，然后用bash read.sh work读取新的原生输出。CPU队列按整节点申请56核，本例小程序实际使用1个MPI进程和1个线程；申请核数不能用来缩减整节点计费。绘图使用ase-3.29.0的Python，HDF5读取使用phonopy-4.8.3的Python，脚本已写好解释器路径。
本例使用TDEP25.03（d38f435）与FourPhonon1.3-fix20260924（5b814b3）的CPU程序，Intel oneAPI2024.2、Intel MPI2021.13、MKL顺序库。四个位移结构的力由VASP6.6.1实算获得：NVHPC25.3、CUDA12.8、HPC-X2.22.1、MKL2023.2，P100/V10016G单卡。
1. 先画原生TDEP能带#
提交脚本里已包含：
module load tdep/25.03-8173m
cd work/TDEP-reference
phonon_dispersion_relations --dumpgrid --qpoint_grid 3 3 3 --meshtype 2 --nq_on_path 5
head outfile.dispersion_relations
路径由原生TDEP确定，为Γ–X–U|K–Γ–L。图和CSV保留原生路径坐标，未把它改成另一套距离定义。--dumpgrid同时输出可与FourPhonon网格匹配的HDF5数据。
2. 比较前先统一单位和q点#
| 原生输出 | 频率定义 | 转换到rad/ps |
|---|---|---|
FourPhonon BTE.omega |
rad/ps | 直接使用 |
TDEP网格HDF5 frequencies |
属性为Hz (angular) |
乘1E−12，不能再乘2π |
本例TDEP路径HDF5 frequencies |
THz普通频率 | 乘2π |
最容易出错的是拿路径频率THz直接与网格角频率比较。公开读取脚本会检查文件属性，再用BTE.qpoints_full逐点匹配，最后才比较分支。
bash read.sh
head data/matched-frequencies.csv
head data/TDEP-reference-path.csv
当前TDEP版本这份网格文件的q点属性写1/A，但实际数值按逆Bohr的周波坐标才能与全部27个网格点一一对应。脚本采用q_frac = q_native @ cell.T / 0.529177210903再按倒格矢周期归一；这只是对当前文件和版本核对后的约定，不能不检查就推广到别的版本。原始HDF5不修改。
3. 为什么Γ点最后一支不同#
本例保留以下原生数值，单位rad/ps：
| 结果 | 光学1 | 光学2 | 光学3 |
|---|---|---|---|
| FourPhonon精确q=0 | 49.288014 | 49.288014 | 49.288014 |
| TDEP方向性Γ极限 | 49.285120 | 49.285120 | 61.808766 |
两条横光学支匹配，方向性纵光学支因LO–TO分裂而不同。这两种Γ点定义不能逐支强行相等；本例非Γ对比单独列出，Γ点也完整保存。不能靠删去最后一支、取平均值或改写原生输出来“对齐”。
4. 换成自己的极性材料#
依次检查：原胞和原子顺序；质量；Born电荷和介电张量；短程/长程分离的参数；q点实际坐标；角频率与普通频率；Γ点取极限方向。确认这些一致后，再分析真正来自超胞、IFC截断和数值积分的差异。
完整逐点CSV包含Γ标记、两程序频率及只用短程IFC的频率；路径CSV可直接重新绘图。资料包同时提供read.py与plot.py，可以查看全部换算过程。