ALAMODE高阶力常数：从原子力拟合三阶与四阶项
用Si64真实GPU原子力，演示位移生成、提力、NORDER与NBODY设置、三阶四阶拟合及新结果接入。
阅读本篇目录
谐波教程主要拟合二阶力常数。如果希望继续研究非谐效应，需要增加三阶、四阶等项。本页用同一个 Si64经典Stillinger–Weber（SW）势模型 演示完整操作；它不是DFT硅的实算结果，也不代表已经算出了热导率或自洽声子。
1. 先理解这次算什么#
| 设置 | 本例内容 |
|---|---|
| 结构 | 金刚石Si，2×2×2常规晶胞，共64原子，a=5.43Å |
| 参考原子力 | LAMMPS双精度Kokkos，P100与V100 16G实际运行 |
| 位移数据 | 32个随机方向组，每组3个位移幅度、正负配对，共192个快照 |
| 拟合数据 | 24个方向组；另外8个方向组独立留出 |
| 对比模型 | 仅FC2、FC2+FC3、FC2+FC3+FC4 |
另外的平衡参考与硬件对照使原始计算共有200个静态力条件；它们服务于同一个高阶力常数工作流，不是200种不同性质。
幅度0.005、0.015、0.03Å指所有原子笛卡尔分量的均方根位移，不是每个原子的最大位移。生成时去除整体平移，随机种子固定为20261105188。图中位移放大5倍便于看清，计算仍使用原始幅度。
2. 下载后先复现已完成的数据#
下载高阶力常数输入、原子力、脚本与结果（ZIP）。在新目录解压：
unzip alamode-higher-order-wave188-*.zip -d alamode-highorder
cd alamode-highorder
cat README.md
进入已分配且支持AVX-512的CPU环境后运行：
module load alamode/1.5.0
bash run.sh
脚本按顺序重建位移结构、调用原生提力工具核对单位、运行ALM拟合、提取验证结果和画图。默认使用附件里真实计算过的原子力，不自动提交GPU作业。结果分别保存在 generated-model/、extracted/、fit/、analysis/、figures/；参考数据保留在 reference/。
3. 位移与原子力怎样准备#
本例位移由ASE根据固定种子生成，可单独执行：
module load alamode/1.5.0
python make-displacements.py
ls generated-model/p100
该命令是默认 run.sh 的第一步；已有 generated-model/ 时会保留目录并停止，重复操作请用新解压目录。它会逐字节核对200个参考结构，不能在只改结构、不改位移元数据的情况下继续套用旧数据。
实际LAMMPS原子力输入为：
package kokkos neigh half
newton on
units metal
atom_style atomic/kk
boundary p p p
read_data structure.lmp
mass 1 28.0855
pair_style sw/kk
pair_coeff * * Si.sw Si
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes
dump f all custom 1 force.dump id xu yu zu fx fy fz
dump_modify f sort id format float %24.17e
thermo 1
thermo_style custom step atoms pe press fmax
run 0
run 0 只计算给定结构的力，不做分子动力学，也不优化掉人为位移。units metal 的力单位为eV/Å；dump_modify ... %24.17e 保留足够数值精度。
要重新计算，先在根目录生成P100批次工作目录：
python3 new-case.py p100
cd work/p100
sbatch run.sh
squeue -u "$USER"
回到附件根目录，对V100批次重复：
python3 new-case.py v100
cd work/v100
sbatch run.sh
两批都完成后，回到根目录接入新力结果：
bash run-new.sh
它检查新输出与配套模型输入，将全部结果写入独立 new-workflow/ 再拟合和画图。原 force-cases/ 参考结果保持不动；一个批次尚未结束时，不要先把部分数据当成完整训练集。
4. 提取位移和力，注意单位#
默认脚本调用：
python extract-forces.py
head extracted/DFSET_native
其中使用ALAMODE自己的 extract.py；单个文件的等价操作示例为：
python /public/software/alamode/1.5.0/8173m/tools/extract.py \
--LAMMPS model/Si64.lammps --unit Rydberg \
force-cases/p100/d00-a005-p/force.dump > DFSET_one
ALM读取的位移单位为Bohr，力为Ry/Bohr。脚本将原生提取结果与原始dump独立对比；extracted/DFSET_native保留原生文本。拟合的数据文件由原始dump做完整精度单位转换，避免提取器文本输出的小数截断。
这里保留正负位移各自的原始力，不像谐波示例那样只取反对称部分。 正负力的偶次分量也含有三阶能量项的信息，全部丢掉会影响高阶拟合。
5. NORDER、NBODY与截断半径#
| 目标 | NORDER | 最高能量项 |
|---|---|---|
| 谐波 | 1 | 二阶FC2 |
| 加入三阶 | 2 | 三阶FC3 |
| 加入四阶 | 3 | 四阶FC4 |
本例四阶模型使用：
&interaction
NORDER = 3
NBODY = 2 3 3
/
&cutoff
Si-Si 8.5 8.5 8.5
/
NBODY限制各阶项里不同原子的体数，不等于Taylor多项式阶数。SW势本身最多含三体相互作用，所以本例四阶项仍允许最多3个不同原子；DFT材料不能直接照抄这个限制。
截断单位是Bohr，8.5Bohr约4.50Å。本例邻近Si–Si键及邻居之间的间距落在这个范围内。换材料、超胞或相互作用模型后，应重新检查范围、数据量和留出误差。
下面是实际四阶全训练输入，完整晶胞与原子位置也包含在内：
展开查看完整fit.in
&general
PREFIX = si64_highorder
FCSYM_BASIS = C
HESSIAN = 1
MODE = optimize
NAT = 64; NKD = 1
KD = Si
/
&interaction
NORDER = 3
NBODY = 2 3 3
/
&cell
20.5224257265890415
1 0 0
0 1 0
0 0 1
/
&cutoff
Si-Si 8.5 8.5 8.5
/
&position
1 0.0000000000000000 0.0000000000000000 0.0000000000000000
1 0.1250000000000000 0.1250000000000000 0.1250000000000000
1 0.0000000000000000 0.2500000000000000 0.2500000000000000
1 0.1250000000000000 0.3750000000000000 0.3750000000000000
1 0.2500000000000000 0.0000000000000000 0.2500000000000000
1 0.3750000000000000 0.1250000000000000 0.3750000000000000
1 0.2500000000000000 0.2500000000000000 0.0000000000000000
1 0.3750000000000000 0.3750000000000000 0.1250000000000000
1 0.0000000000000000 0.0000000000000000 0.5000000000000000
1 0.1250000000000000 0.1250000000000000 0.6250000000000000
1 0.0000000000000000 0.2500000000000000 0.7500000000000000
1 0.1250000000000000 0.3750000000000000 0.8750000000000000
1 0.2500000000000000 0.0000000000000000 0.7500000000000000
1 0.3750000000000000 0.1250000000000000 0.8750000000000000
1 0.2500000000000000 0.2500000000000000 0.5000000000000000
1 0.3750000000000000 0.3750000000000000 0.6250000000000000
1 0.0000000000000000 0.5000000000000000 0.0000000000000000
1 0.1250000000000000 0.6250000000000000 0.1250000000000000
1 0.0000000000000000 0.7500000000000000 0.2500000000000000
1 0.1250000000000000 0.8750000000000000 0.3750000000000000
1 0.2500000000000000 0.5000000000000000 0.2500000000000000
1 0.3750000000000000 0.6250000000000000 0.3750000000000000
1 0.2500000000000000 0.7500000000000000 0.0000000000000000
1 0.3750000000000000 0.8750000000000000 0.1250000000000000
1 0.0000000000000000 0.5000000000000000 0.5000000000000000
1 0.1250000000000000 0.6250000000000000 0.6250000000000000
1 0.0000000000000000 0.7500000000000000 0.7500000000000000
1 0.1250000000000000 0.8750000000000000 0.8750000000000000
1 0.2500000000000000 0.5000000000000000 0.7500000000000000
1 0.3750000000000000 0.6250000000000000 0.8750000000000000
1 0.2500000000000000 0.7500000000000000 0.5000000000000000
1 0.3750000000000000 0.8750000000000000 0.6250000000000000
1 0.5000000000000000 0.0000000000000000 0.0000000000000000
1 0.6250000000000000 0.1250000000000000 0.1250000000000000
1 0.5000000000000000 0.2500000000000000 0.2500000000000000
1 0.6250000000000000 0.3750000000000000 0.3750000000000000
1 0.7500000000000000 0.0000000000000000 0.2500000000000000
1 0.8750000000000000 0.1250000000000000 0.3750000000000000
1 0.7500000000000000 0.2500000000000000 0.0000000000000000
1 0.8750000000000000 0.3750000000000000 0.1250000000000000
1 0.5000000000000000 0.0000000000000000 0.5000000000000000
1 0.6250000000000000 0.1250000000000000 0.6250000000000000
1 0.5000000000000000 0.2500000000000000 0.7500000000000000
1 0.6250000000000000 0.3750000000000000 0.8750000000000000
1 0.7500000000000000 0.0000000000000000 0.7500000000000000
1 0.8750000000000000 0.1250000000000000 0.8750000000000000
1 0.7500000000000000 0.2500000000000000 0.5000000000000000
1 0.8750000000000000 0.3750000000000000 0.6250000000000000
1 0.5000000000000000 0.5000000000000000 0.0000000000000000
1 0.6250000000000000 0.6250000000000000 0.1250000000000000
1 0.5000000000000000 0.7500000000000000 0.2500000000000000
1 0.6250000000000000 0.8750000000000000 0.3750000000000000
1 0.7500000000000000 0.5000000000000000 0.2500000000000000
1 0.8750000000000000 0.6250000000000000 0.3750000000000000
1 0.7500000000000000 0.7500000000000000 0.0000000000000000
1 0.8750000000000000 0.8750000000000000 0.1250000000000000
1 0.5000000000000000 0.5000000000000000 0.5000000000000000
1 0.6250000000000000 0.6250000000000000 0.6250000000000000
1 0.5000000000000000 0.7500000000000000 0.7500000000000000
1 0.6250000000000000 0.8750000000000000 0.8750000000000000
1 0.7500000000000000 0.5000000000000000 0.7500000000000000
1 0.8750000000000000 0.6250000000000000 0.8750000000000000
1 0.7500000000000000 0.7500000000000000 0.5000000000000000
1 0.8750000000000000 0.8750000000000000 0.6250000000000000
/
&optimize
DFSET = DFSET_train
LMODEL = ols
ICONST = 11
NDATA = 144
/
FCSYM_BASIS=C使用笛卡尔基底，ICONST=11施加平移约束。本例采用OLS，脚本一共生成三个阶数、全训练及三折组合的12个原生拟合目录。
6. 查看力常数输出#
cd fit/order3-full
ls si64_highorder.xml si64_highorder.fcs si64_highorder.hessian
grep 'Number of free' fit.log
.xml供ANPHON后续计算读取；.fcs保留力常数条目；.hessian是二阶Hessian矩阵，不是四阶张量。怎样判断模型是否学到了未参与拟合的数据，继续看分组交叉验证与留出检查。
ALAMODE 1.5.0采用Intel oneAPI 2024.2、Intel MPI和MKL，当前CPU构建包含AVX-512；参考原子力采用LAMMPS 30Sep2026、Kokkos 5.2.1、CUDA 12.8双精度路径。Python结构与分析操作使用现有ASE环境。官方高阶力常数教程、ALM输入说明