Fourthorder重构与位移量检查
原生reap重构四阶力常数，比较三个位移量并用独立构型检查受力预测。
阅读本篇目录
下载Si真实结构、受力、力常数、8张图和终端脚本，保存为 si-fourphonon.zip。所有本系列文章使用同一个包：
python3 -m zipfile -e si-fourphonon.zip .
cd fourthorder-fourphonon-silicon-wave288
bash preview.sh
ls preview/figures
包中 reference/ 保存本次实算参考；work/ 是您重新生成的文件，results/data/ 是提取后的CSV。受力来自真实VASP计算；下载包不含POTCAR、波函数、完整电子输出或第三方程序。
1 从真实受力重建#
bash prepare.sh
bash reconstruct.sh
head -n 12 work/fourthorder/H040/FORCE_CONSTANTS_4TH
head -n 5 work/thirdorder/H010/FORCE_CONSTANTS_3RD
已经执行过prepare/reconstruct时，直接查看结果，不必再次创建目录。原生VASP接口的核心命令为：
module load fourthorder/0.1-fix20260618
cd work/fourthorder/H040
find "$PWD" -maxdepth 1 -name 'H040-*.forces.xml' | sort > ordered-forces.txt
python ../../programs/H040/Fourthorder_vasp.py reap 2 2 2 -1 < ordered-forces.txt
上述 ordered-forces.txt 是按.01至.48列出的XML路径；配套重建脚本使用相同顺序经标准输入传入，您自己收集的新DFT列表则是前篇生成的 fourthorder-xml.list。请先准备正确的列表再运行原生命令。
本例每组产生58个四原子组合块，每块含3⁴=81个笛卡尔分量。文件中的晶格平移使用Å，四阶常数使用eV/Å⁴。它是势能的四阶导数；用力做差分时必须考虑负号，不能直接把它当成力的四阶导数。
2 为什么比较H#
四阶重构含三次力差分，受力误差会被H³放大。H太小会放大电子噪声，太大又可能混入更高阶响应。因此不能简单认为“越小越准确”。
| H（Å） | 最大平移和规则残差（eV/Å⁴） | 用于本例输运 |
|---|---|---|
| 0.01 | 8.51×10⁻⁶ | 否，保留作噪声对照 |
| 0.02 | 8.60×10⁻⁸ | 比较用 |
| 0.04 | 1.23×10⁻⁸ | 本例选用 |
本例采用1E-6 eV/Å⁴的平移和检查尺度。H=0.01 Å在这一尺度上偏差较大，不将其用于热导；程序输出保留原值，没有另施加人为归零。H=0.02与0.04 Å的完整张量相对差约0.307%，只能说明这两组接近，不能由此宣布所有位移量、超胞与截断都已收敛。
3 用没参与所选重构的结构检查#
H040四阶与H010三阶重构之外，另选36个实际H020构型，以二、三、四阶展开预测受力，再对比DFT值。力分量RMSE为 0.00018555 eV/Å，相对RMSE约0.197%。仅用二、三阶时RMSE约0.00021338 eV/Å。
这一步检查截断模型在本次小位移附近是否能描述实际受力，不等于证明高温大振幅、远邻作用或材料热导已经准确。
4 二阶和三阶从哪里来#
三阶用 thirdorder/1.1.4 对同一原胞、2×2×2超胞、最近邻截断生成8个逻辑位移并原生重构，得到26块三阶常数。二阶复用同一Si16、同PAW/450 eV/3³电子k点的真实谐波DFT，显式转换到本例原子顺序。VASP原始 force_constants 是dF/dR，接入能量二阶导数时取负号；质量加权Hessian也不能直接当作eV/Ų的常数。
二阶文件名为 FORCE_CONSTANTS_2ND，三阶为 FORCE_CONSTANTS_3RD，四阶为 FORCE_CONSTANTS_4TH。下一步按FourPhonon入门组合使用。