LAMMPS POD描述符与导数矩阵

本组实际完成一个小型POD机器学习势流程：为32个Cu108近FCC结构计算ASE EMT参考能量和力，24个用于训练、8个独立测试；调用LAMMPS原生 fitpod，再做CPU/P100推断、描述符和导数检查。

参考标签是半经验ASE EMT，不是DFT。 这个例子用于学习输入格式和工具流程，不是可直接用于任意铜体系的第一性原理势。原生fitpod日志沿用 DFT energy/force 列名，本包这些列实际对应EMT参考值。

下载全部输入、参考标签、模型、原始结果与图片，在个人目录执行：

mkdir pod-demo
cd pod-demo
curl -fL https://yuansuan.top/assets/downloads/full-function/lammps-pod-wave81-13e715e5efcc.zip -o examples.zip
unzip examples.zip
sbatch submit.slurm

脚本申请一个CPU核、一个P100，自动生成标签、训练、分析已有模型、CPU与Kokkos推断并画图。模块为 ase/3.29.0、lammps/30Sep2026-p100-gpu 和 lammps/30Sep2026-p100-kokkos；LAMMPS基础构建使用GCC13.3、OpenMPI5.0.10、CUDA12.8，Kokkos推断为双精度。

完成后在整包目录运行：

module purge
module load ase/3.29.0
python analyze.py work
python export_descriptors.py work
python plot_results.py work
head -10 work/analysis/heldout-errors.csv

reference/ 保存已有实算，可把后处理命令的 work 换成 reference 先看结果。新计算使用新解压目录。整包一直采用单MPI进程；本版本四种POD描述符compute也要求单MPI，不要直接把它们改成多MPI提交。

1 用小模型看清矩阵#

这里另用4个原子的周期性FCC胞，并把一个原子位移 (0.03,−0.02,0.01) Å，避免完美晶体的零力掩盖符号错误。仍使用前面训练的16系数POD模型；这是描述符代数检查，不是对4原子材料模型的独立精度结论。

原生输入的完整描述符输出部分如下，后面的六个位移单点保存在下载包中：

units metal
atom_style atomic
atom_modify map array
boundary p p p
read_data structure.data
pair_style pod
pair_coeff * * params.pod coefficients.pod Cu
neighbor 1.0 bin
compute ea all pe/atom
compute dp all pod/atom params.pod coefficients.pod Cu
compute dd all podd/atom params.pod coefficients.pod Cu
compute dg all pod/global params.pod coefficients.pod Cu
compute dl all pod/local params.pod coefficients.pod Cu
fix globalout all ave/time 1 1 1 c_dg[*] file global.dat mode vector format " %.16g"
fix localout all ave/time 1 1 1 c_dl[*] file local.dat mode vector format " %.16g"
dump descriptors all custom 1 descriptors.dump id c_ea c_dp[*] c_dd[*]
dump_modify descriptors sort id format float %.16g
thermo_modify norm no format float %.16g
run 0
print "CENTER $(pe:%.16g)" file center.dat screen no
write_dump all custom forces.dump id type x y z fx fy fz c_ea modify sort id format float %.16g
undump descriptors
unfix globalout
unfix localout

四种compute在本版本均使用单MPI进程；pod/local 名字中的local指按中心原子分解描述符，其返回值是全局数组，本例用 fix ave/time ... mode vector 输出，不用 dump local。

2 四种实际输出的形状#

一体常数系数另计，因此全局模型有16个系数。其他元素数、环境簇或描述符参数会改变维度；这里的数组读取与核对专门对应这个单元素、单环境簇模型。

podd/atom 的内存随原子数平方和描述符数增长，大体系不要为了“多输出”全部开启。实际只需要能量/力预测时，保留 pair_style pod 就够了。

head -6 work/cases/descriptor-control/global.dat
head -6 work/cases/descriptor-control/local.dat
head -11 work/cases/descriptor-control/descriptors.dump

3 原生导数采用力的符号#

这份30Sep2026实际输出中，原生导数矩阵对应负坐标梯度。pod/global 第一行乘系数给总能量，其余行乘系数直接给力，不再额外加一个负号。我们用六个真实位移单点，逐描述符核对了这一符号，并用独立能量差分核对实际力。

图中力来自原生POD计算、设计矩阵乘系数和独立能量差分，三者一致。导数逐项差分的最大差约 1.37e-8（本例原生描述符单位/Å）。

4 导出可以直接读的CSV#

module load ase/3.29.0
python export_descriptors.py work
head -5 work/analysis/pod-global-design.csv
head -5 work/analysis/pod-force-signed-derivatives.csv

工具得到全局设计矩阵、720条带中心/目标原子/方向/描述符索引的力符号导数，以及16个模型系数。下面直接用CSV重构能量和力：

python - <<'PY'
import numpy as np
c=np.loadtxt('work/analysis/pod-coefficients.csv',delimiter=',',skiprows=1,usecols=1)
G=np.loadtxt('work/analysis/pod-global-design.csv',delimiter=',',skiprows=1,usecols=range(1,17))
values=G@c
print('Energy / eV:',values[0])
print('Forces / (eV/Angstrom):')
print(values[1:].reshape(4,3))
PY

逐原子能量还可由 c0 + D_i·c_rest 重构。完整原生dump、模型参数和位移输入都保留在包内，便于对照索引，不要跨参数文件混用描述符与系数。

5 检查自己的修改#

python analyze.py work
python check_results.py work
python plot_results.py work

默认脚本检查数据来源、独立测试、CPU/P100一致性、四类矩阵维度、力符号及能量/描述符差分。更换模型或体系后，需要同步修改对应假设，不能把本例16列和4个原子的读取方式当作所有POD文件的通用格式。

继续：POD训练、模型预测。原生接口与模型说明见官方POD页面。