thirdorder 收集原子力与三阶常数
用真实VASP XML按位移顺序运行reap，理解重复结构、原子顺序、导数符号和输出单位。
阅读本篇目录
reap读取每一个位移结构的原子力，按对称性和有限差分关系重构三阶力常数。最关键的是结构、原子顺序、位移编号和电子设置匹配。
下载完整Si实例包。在新目录解压后操作；包内有真实VASP输入和原始力输出、位移映射、三阶/二阶力常数、ShengBTE结果、CSV和高清图，不含POTCAR。
1. 先完成sow#
bash prepare.sh
ls workflow/sow/si16-c0p30nm/3RD.POSCAR.*
完整包已经保存本集群实际VASP计算的原始XML与OUTCAR，可先用它们练习力收集：
bash reap.sh
脚本按 DFT-MAP.json查找每份位移对应的真实DFT输出，逐份核对POSCAR和电子设置，再生成有序XML清单。它不把人工排好的数字伪装成VASP输出。
2. 底层reap命令#
cd workflow/sow/si16-c0p30nm
cat vasprun-files.txt
thirdorder_vasp.py reap 2 2 2 0.3 < vasprun-files.txt
sow和reap的重复数、截断必须一致。清单每行一份XML，顺序对应 3RD.POSCAR.01、.02等原生编号。不要按随机生成的计算目录名排序后直接交给reap；程序并不会替你判断文件列表中的结构是否放错。
3. 为什么有的XML会重复出现#
原生位移清单中可能存在抵消后相同的结构。本例112个位移槽位包含98种不同的POSCAR字节；加上匹配谐波位移，共102个不同DFT条件。相同POSCAR、INCAR、KPOINTS和赝势才可以共用一次求力。
XML在清单中重复出现，表示两个位移编号确实对应同一份结构与设置。不能因为原子数相同或总能接近，就把不相同的结构合并。
4. 读取输出#
head FORCE_CONSTANTS_3RD
首行为相互作用块数。每块包含平移向量、原胞原子编号和27个Cartesian张量分量。本例原生块数为26、182和266，分别对应Si16/3 Å、Si16/4 Å和Si54/4 Å。
| 内容 | 单位或约定 |
|---|---|
| 输出平移向量 | Å |
| 三阶力常数 | eV/ų |
| 分量编号 | 1、2、3对应x、y、z |
| 输入cutoff正数 | nm，不能与输出Å混用 |
原子力是能量梯度的负值，因此从力的二阶变化得到能量三阶导数时包含负号。配套检查直接比较原始四组VASP原子力差分与三阶张量的预测。
三个有限模型的全体混合力曲率相对RMSE约2.26%、1.25%和2.09%。这是有限差分、截断及对称性重构后的数值对照，不是对材料热导率误差的估计。
5. 二阶常数不要随便拼接#
本包同时用相同电子设置的正负谐波位移求二阶常数，保存Phonopy原生顺序的 FORCE_CONSTANTS_2ND。三阶块按原胞与平移向量描述相互作用，不能简单把三阶数组的原子行号当成Phonopy超胞行号。
脚本通过实际周期几何和格式读取器核对映射；更换结构、超胞或原子排列后，应重新核对整个对应关系。ShengBTE接入继续使用成套二、三阶常数。