Pyramid(pyraimd2)背后的核心方法：energetic力误差、方向性响应与符号残差功完全解析

Pyramid（Python 包名 pyraimd2）是一个把第一性原理计算与机器学习势（MLIP）结合起来的自适应分子动力学框架。它的核心是一套 energetic 力误差方法：先定义锚定后的力残差，再用方向性响应预测力误差沿原子运动方向的增长，最后用符号残差功量化误差对轨迹的真实能量影响。这三个量共同决定了"哪一步可以用廉价势、哪一步必须请求参考计算"，是理解整个 Pyramid 自适应循环的钥匙 🔑。

为什么力误差值得"energetic"地对待？

自适应第一性原理分子动力学（AIMD）的基本分工是：

• 参考计算（如 DFT，Quantum ESPRESSO / PySCF）——慢而准，提供"标准答案"；
• 快速势（如 MACE 机器学习势）——快而近似，驱动绝大多数原子运动。

问题在于：快速势给出的力与参考力不一致时，误差到底"伤"了轨迹多少？只看误差数值大小是不够的——两个力误差相同的方向，一个可能让误差迅速放大，另一个几乎无感。Pyramid 的答案是同时跟踪三件事：

量	回答的问题	关键表达式	实现位置
力误差 e(X)	力与参考力差多远？	e(X) = max_i |R_i(X)|	response.py
方向性响应 q / 曲率 kappa	沿运动方向，误差涨多快？	kappa = u·q	response.py
符号残差功 W	误差实际注入了多少能量？	W = ΔU_ref − ΔU_base + c·ΔX	work.py

完整方法说明见项目自带的指南 docs/energetic_force_error.md，本文按"锚定 → 响应 → 功"的顺序逐步拆解。

第一步：参考锚定——力误差是怎么定义的 🎯

Pyramid 在参考构型 X0 处计算一个常数修正：

c = F_ref(X0) − F_base(X0)
F_anchor(X) = F_base(X) + c
R(X) = F_anchor(X) − F_ref(X)      （残差力）
e(X) = max_i |R_i(X)|              （力误差度量）

• 在锚点处残差被"清零"（R(X0) = 0），整个片段（segment）内 c 保持不变；
• 这里的模长是每个原子三维向量的欧氏模长再取最大，即"最不准的那个原子"说了算；
• 若启用 FixAtoms 约束，默认度量 active_dofs_max_atom 只在自由坐标上统计，固定原子不贡献位移功。

可以把它类比为用标准砝码校准天平：先在已知位置"调零"，之后讨论的偏差才有物理意义。如果基势的力是某能量的梯度，则存在对应的锚定势能 U_anchor(X) = U_base(X) − c·(X − X0)——这正是后面能量差的来源。

第二步：方向性响应——提前预测误差的增长 🧭

锚定只修正了当前点，原子一旦运动，残差就会增长。增长快慢取决于朝哪个方向运动，这就是"方向性"的含义：

1. 选一个单位方向 u（默认沿原子速度方向，也可用回调提供多个方向）；
2. 沿 u 的正负位移做成对中心探针，得到残差的方向导数：

q = [R(X0 + h·u) − R(X0 − h·u)] / (2h)

3. 用两个探针尺度（如 h1 < h2）估计响应及其对尺度的敏感度；较大的尺度给出最终响应 q，两尺度之差进入包络系数 eta。

由此导出两个关键量：

• 方向曲率 kappa = u·q：符号保留了残差沿该运动是"注入能量"还是"抽走能量"，是符号信息的起点；
• 方向性残差功系数 C_rw = kappa / g²（单位 Ų/eV，g = max|q_i|）：把主导方向功联系到力误差的平方。

对沿 u 的位移 a·u，预报为：

e_linear   = |a| · g          （线性力误差）
W_predicted = a² · kappa / 2   （预测符号功）

运行时在请求昂贵的参考计算之前，就用这些量预测"这一段运动会积累多少力误差、带来多少能量扰动"。若预报满足预设的力误差预算（force_budget，单位 eV/Å）与域约束（时限 time_cap_fs、横向偏离 transverse_cap），就继续用锚定快速力；否则请求新的参考标签并重新锚定。数值入口是 src/pyraimd2/energetics/response.py 中的 estimate_responses()，每个方向封装为 DirectionalResponse，其 forecast() 方法返回 linear_error、envelope、predicted_work、admitted 等字段，供运行时做准入决策。

注意：预报是有限探针推得的经验包络，不是对每一步力的认证；探针调用本身也计入参考成本。

第三步：符号残差功——能量影响的最终裁判 ⚖️

每当轨迹到达一个新的参考端点（重新锚定），即可用两端能量差计算该片段上残差实际做的有符号功：

W(X) = [U_ref(X) − U_ref(X0)] − [U_base(X) − U_base(X0)] + c·(X − X0)

它的要点：

• 符号保留能量转移的方向：正值表示残差向系统注入能量，负值表示抽出——"误差在帮倒忙还是在帮忙"一目了然；
• 若片段内有成对力标签，还可用梯形积分 Q = Σ (R_j + R_{j+1})/2 · ΔX_j 交叉验证，两者一致说明力求积的分辨率足够；
• 在保守 NVE 片段中，参考哈密顿量变化可拆分为 ΔH_ref = W + ΔH_anchor，把残差的能量效应与积分本身的漂移干净分离。

实现见 src/pyraimd2/energetics/work.py：residual_work() 基于端点能量，integrate_residual_work() 基于形状 (T, N, 3) 的成对力数组，返回逐原子累计功。

独立检查：另一条统计验证通道 🔍

预报之外，Pyramid 还支持独立的已接受力检查：先固定将要采用的力，再按固定概率 p 抽取是否请求参考来对照。发现违规时，verification.py 中的 IndependentCheckBound 会随观测积累"接受力违规比例"的统计上界（置信水平 1−delta 下对所有前缀同时成立）。它与方向性预报相互独立、可分别阅读，共同构成可复现的精度报告。

这套方法在代码库中的位置 📁

energetic 力误差方法是纯数值层（仅依赖 NumPy），不发起参考计算、不控制轨迹，因此可以脱离 MD 单独使用：

• 数值核心：src/pyraimd2/energetics/
  • response.py：方向性响应、C_rw、forecast() 预报
  • work.py：端点功与逐原子积分功
  • verification.py：独立检查上界
• 动力学与自适应循环：EnergeticCalculator / EnergeticRunner，位于 src/pyraimd2/loop/
• 方法全文：docs/energetic_force_error.md；整体架构：docs/architecture.md
• 可复现数据包：reproducibility/force_error/——含上千个轨迹片段的功数据、原子级功分解、硅体系尺寸标度等数值记录（运行 python -B check_all.py 即可自检重建）

五分钟跑通一个完整示例 🚀

仓库提供了一个秒级完成的离线演示 examples/energetic_loop.py：内置解析的谐振荡参考势与快速势，无需任何外部程序或模型文件，即可跑一段 energetic 门控的 NVE 轨迹，并在数据库里记录校准成本、预报、端点功与实际驱动力：

uv run python examples/energetic_loop.py --output results/energetic

阅读它的 write_report() 就能看清每次参考调用被分账为"锚定 / 探针 / 检查"三类成本，这正是理解自适应循环开销结构的最好入口。

适用范围与使用须知 ⚠️

• 能量口径（锚定能量与端点功）要求每个片段内力与其能量自洽（保守力）；方法指南中声明的完整支持范围为固定晶胞 NVE（可带 FixAtoms），变胞与其他约束暂不支持；
• 力误差预算、探针尺度、数值下限、时限与横向上限都是运行前设定的操作参数，属于经验包络而非逐力认证；
• 热机、约束与变胞做功需要各自独立的记账，不要混入 W 的解读；
• 单位遵循 ASE 约定：位置 Å、能量 eV、力 eV/Å、时间 fs；
• 系数 C_rw 要求方向性响应非零，退化时 estimate_responses 会抛出 DegenerateResponseError 提醒。

小结

Pyramid（pyraimd2）把"力误差"从单一数值升级为一套有方向的能量叙事：锚定定义误差、方向性响应预测误差的增长与符号、符号残差功裁定误差对轨迹的最终能量贡献。理解了这三层，你就掌握了 docs/energetic_force_error.md 中全部核心公式的物理动机，也能看懂每次自适应运行中"为什么这一步用了快速势、那一步请求了参考计算"。