1. 项目概述：当机器学习势函数遇上局部压力计算
在分子动力学模拟的世界里，压力或应力张量是连接微观原子运动与宏观材料力学性能的桥梁。无论是研究金属的塑性变形、聚合物的粘弹性，还是分析血液在微血管中的流动，我们最终都需要从原子轨迹中提取出可靠的应力信息。传统上，维里应力张量因其计算简便，成为了绝大多数模拟软件（如LAMMPS、GROMACS）的默认选择。然而，当系统变得不均匀——例如在固液界面、裂纹尖端、剪切流动前沿，或是任何远离热力学平衡的状态时——维里应力的物理基础就开始动摇，其计算结果可能严重偏离真实情况。
这就像用一把尺子去测量一个扭曲表面的曲率，工具本身就不对。对于非平衡分子动力学模拟和流体动力学研究而言，一个能在局部严格满足动量守恒的应力定义至关重要。这时，“平面方法”就登场了。它的核心思想极其物理直观：应力就是力除以面积。MoP通过统计穿过空间中一个假想平面的所有原子间作用力，来定义该平面上的应力。这个定义直接源于Irving和Kirkwood的统计力学奠基性工作，天生就与动量守恒方程紧密耦合。
近年来，机器学习势函数的爆发式发展，特别是像MACE这样的模型，让我们能够以前所未有的精度和速度进行分子模拟。MACE等模型通过图神经网络和消息传递架构，巧妙地融合了量子力学的精度与经典分子动力学的效率。但是，如何为这些复杂的、非两体的、蕴含量子精度的相互作用力，定义一个正确且可计算的局部应力？这正是本文要解决的核心问题。我们将深入探讨如何将平面方法拓展到MACE势函数，并验证其在非平衡体系中的正确性，为复杂界面和流动问题的精确模拟铺平道路。
2. 核心原理：从统计力学基础到平面方法
2.1 传统维里应力的局限与平面方法的物理图像
要理解平面方法的优势，首先得看清维里应力在非均匀体系中的问题所在。维里应力的常见形式（如IK1近似）将原子间相互作用对压力的贡献，简单地按“一人一半”的方式分配给相互作用的两个原子。在均匀流体中，这种平均化处理没有问题。但在界面附近，一个原子在液相，另一个在固相，它们之间的相互作用本质上是跨越界面的“表面力”。将这份力对半劈开，分别计入液相和固相的压力，物理上是不清晰的，会导致界面处的应力分布出现人为的模糊甚至错误。
平面方法则回归了应力的力学本源：单位面积上的力。想象在流体中放置一个无限薄的平面（例如z=z0平面）。任何试图改变该平面一侧流体动量的企图，都必须通过这个平面传递力。因此，统计所有在时间Δt内穿过该平面的分子间作用力，再除以平面的面积和统计时间，就得到了该平面上的应力。这种方法不依赖于对相互作用路径的人为分割，只要两个原子分居平面两侧，它们之间的作用力就会完整地贡献给该平面的应力。这使得MoP在物理上尤为干净，尤其适合分析界面、剪切层等梯度剧烈的区域。
2.2 Irving-Kirkwood公式与控制体积法
平面方法的严格数学基础来自Irving和Kirkwood于1950年建立的统计力学框架。他们从相空间中的刘维尔方程出发，推导出了质量、动量和能量守恒方程的微观表达式。对于动量方程，其关键步骤是将系统总势能U对位置的导数（即力）的统计平均，表达为一个应力张量的散度。
其中，处理原子间相互作用力贡献的部分，涉及一个著名的技巧：将两个狄拉克δ函数之差（δ(r−r_i) − δ(r−r_j)）展开。一种常见的处理是采用所谓的“IK路径积分”，假设相互作用沿原子i和j的连线直线传播。然而，对于MACE这样的高阶多体势，相互作用的“路径”在物理上并非简单的直线，而是隐含在图神经网络复杂的消息传递中。强行指定一条直线路径（如IK轮廓）会引入不必要的假设。
一个更稳健、且与守恒律直接挂钩的方法是采用控制体积法。考虑空间中的一个有限体积V，其边界为闭合曲面S。根据散度定理，体积内的动量变化率等于通过表面S进入该体积的净动量通量（即表面应力）。通过对Irving-Kirkwood公式在控制体积V内积分，我们可以得到一个极其有力的关系式：
[ \frac{\partial}{\partial t} \sum_{i=1}^{N} \langle \mathbf{p}i \vartheta_i \rangle = -\sum{\beta \in \text{faces}} \iint_{S_\beta} \mathbf{P}^c \cdot d\mathbf{S}_\beta ]
这里，(\vartheta_i)是一个指示函数，当原子i在体积V内时为1，否则为0。(\mathbf{P}^c)是构型应力张量。这个等式的物理意义非常明确：控制体内总动量的变化，完全由穿过其各个表面β的应力所决定。平面方法正是这个一般性结论的一个特例：当我们的控制体积是由两个平行的无限大平面（如z+和z-平面）所夹的平板区域时，侧面对动量的贡献在统计平均下为零（假设系统在x，y方向均匀或周期性边界），动量变化就简化为两个主平面上的应力差：
[ \frac{\partial}{\partial t} \sum_{z^- < z_i < z^+} \mathbf{p}i = \mathbf{P}{z^+} - \mathbf{P}_{z^-} ]
这个公式是验证任何局部应力定义是否正确的“试金石”。一个正确的局部应力定义，必须在其对应的控制体积上，严格满足上述瞬时动量守恒（无需系综平均）。我们将在后续的验证部分看到，平面方法应力满足这一强约束，而局部化的维里应力则不满足。
2.3 MACE势函数中的“成对力”分解
MACE是一个典型的消息传递神经网络势函数。它的能量预测依赖于以每个原子为中心的局部原子环境构成的图结构。虽然能量表达式本质上是多体的（MACE-MP-0的有效体阶高达13），但其对原子位置的梯度（即原子所受的力）总可以形式上写成一个反对称的成对力和的形式：
[ \mathbf{F}i = -\frac{\partial U}{\partial \mathbf{r}i} = \sum{j \neq i} \mathbf{f}{ij}, \quad \text{其中} \quad \mathbf{f}{ij} = \frac{\partial U}{\partial \mathbf{r}{ij}} - \frac{\partial U}{\partial \mathbf{r}_{ji}} ]
这里的关键在于 (\partial U / \partial \mathbf{r}{ij} \neq - \partial U / \partial \mathbf{r}{ji})。这是因为在MACE的图网络中，从原子i看向j的“视角”与从j看向i的“视角”所依赖的局部环境特征（通过消息传递聚合的信息）是不同的。因此，(\mathbf{f}{ij}) 不能简单地解释为“原子j施加在原子i上的力”，而应理解为“由于原子j的存在而对原子i所受力的贡献”。尽管如此，由上述定义可以自动保证 (\mathbf{f}{ij} = -\mathbf{f}_{ji})，从而满足整个系统的线动量守恒（牛顿第三定律）。
注意：这种分解方式并非唯一。另一种常见的分解是 (\tilde{\mathbf{f}}{ij} = -\partial U_j / \partial \mathbf{r}i)，其中 (U_j) 是分配给原子j的那部分势能。这种分解能自然地满足系统能量守恒，但通常不满足 (\tilde{\mathbf{f}}{ij} = -\tilde{\mathbf{f}}{ji})。在平面方法的框架下，只要力是成对定义的且满足作用力与反作用力，就可以用于动量通量的计算。而能量通量的计算则需要使用后一种分解。本文主要关注动量与应力，因此采用前一种定义。
实操心得：在代码实现中，获取 (\mathbf{f}{ij}) 需要利用自动微分工具。以PyTorch为例，在计算得到系统总能量U后，我们需要对每一对相互作用矢量 (\mathbf{r}{ij}) 进行求导，而非直接对原子坐标 (\mathbf{r}i) 求导后再进行组合。这是因为自动微分库在反向传播时，沿着图结构计算的梯度会自然汇聚到发送者（sender）和接收者（receiver）节点上，我们需要显式地捕获这些边（edge）上的梯度。一个常见的做法是构建一个N×N×3的稠密或稀疏张量来存储所有 (\partial U / \partial \mathbf{r}{ij})，然后通过上述反对称操作得到 (\mathbf{f}_{ij})。虽然这会牺牲一些MACE利用图结构带来的计算效率，但由于应力计算通常不是每一步都进行（而是在采样间隔进行），这部分开销是可以接受的。
3. 平面方法应力在MACE中的实现与计算
3.1 应力张量的具体计算公式
对于一个垂直于z轴、位于 (z = z_0) 的平面，平面方法给出的应力向量 (\mathbf{P}^{MoP}_{z_0})（包含x, y, z三个方向的分量）由两部分构成：动能部分和构型（势能）部分。
[ \mathbf{P}^{MoP}{z_0} = \frac{1}{A} \left\langle \sum{i=1}^{N} \frac{\mathbf{p}i p{i,z}}{m_i} \delta(z_i - z_0) \right\rangle + \frac{1}{4A} \left\langle \sum_{i=1}^{N} \sum_{j \neq i} \mathbf{f}{ij} , d{ij}(z_0) \right\rangle ]
其中：
(A) 是平面在x-y方向的有效面积（对于周期性系统，就是模拟盒子的横截面积 (L_x \times L_y)）。
第一项是动能贡献，来源于原子自身运动携带的动量穿过平面。只有当原子i恰好位于平面 (z_0) 上时，其贡献才不为零。在实际的离散化计算中，我们通常会将空间划分为薄层（bin），将位于某个薄层内的原子的贡献平均到该层对应的平面上。
第二项是构型贡献，来源于原子间相互作用力。(d_{ij}(z_0)) 是一个符号函数： [ d_{ij}(z_0) = \text{sgn}(z_0 - z_j) - \text{sgn}(z_0 - z_i) ] 其取值只有三种可能：
- (d_{ij} = +2)：原子i在平面下方 ((z_i < z_0))，原子j在平面上方 ((z_j > z_0))。
- (d_{ij} = -2)：原子i在平面上方，原子j在平面下方。
- (d_{ij} = 0)：两个原子位于平面的同一侧。
因此，求和项 (\sum_{i,j} \mathbf{f}{ij} d{ij}) 实际上只统计那些“跨越”了目标平面的原子对。因子 (1/4) 的出现，是因为每个相互作用在遍历所有原子对时被重复计算了两次（i-j和j-i），并且 (d_{ij}) 的绝对值为2。
3.2 计算流程与代码实现要点
基于上述公式，在MACE分子动力学模拟中计算平面方法应力的流程可以概括如下：
模拟与采样：使用MACE势函数进行NEMD或平衡MD模拟。在系统达到稳态后，以远大于应力自相关时间的时间间隔采集快照（snapshots）。对于界面体系，需要确保模拟时间足够长，以获得良好的统计平均。
力分解：对于每一帧快照，利用自动微分计算系统总能量U对每一对相互作用边 (\mathbf{r}{ij}) 的梯度，得到张量 (\partial U / \partial \mathbf{r}{ij})。然后通过反对称操作 (\mathbf{f}{ij} = \partial U / \partial \mathbf{r}{ij} - (\partial U / \partial \mathbf{r}{ji})^T) 得到成对力。注意确保力满足 (\mathbf{f}{ij} = -\mathbf{f}_{ji})，这是动量守恒的保证。
空间划分与平面指定：沿着需要计算应力分布的方向（如垂直于界面的z方向），将空间划分为一系列等间距的薄层（slabs）或直接定义一系列z坐标的平面。平面的面积A由模拟盒子的x和y维度决定（需考虑周期性边界）。
贡献累加：
- 动能部分：遍历所有原子。对于原子i，确定其所在的薄层索引k（对应平面 (z_k)）。将其动量通量 (\mathbf{p}i p{i,z} / m_i) 累加到该平面k的动能应力累加器中。
- 构型部分：遍历所有非键相互作用对 (i, j)。对于每一对，检查其连线是否跨越了某个平面 (z_k)。这可以通过判断 (z_i) 和 (z_j) 是否位于 (z_k) 两侧来完成。如果跨越，则计算 (d_{ij}(z_k) = \pm 2)，并将力 (\mathbf{f}{ij}) 乘以 (d{ij}/4) 后，累加到平面k的构型应力累加器中。注意：一个原子对可能跨越多个薄层（如果薄层厚度小于原子间距离），在离散化时需要妥善处理，通常采用线性插值将贡献分配到沿途穿过的所有平面上。
平均与输出：对所有采集的快照重复步骤2-4。最后，将每个平面上的累加器除以总采样帧数、平面面积A，即得到该平面上的平均应力向量 (\mathbf{P}^{MoP}(z))。完整的应力张量需要通过分析不同方向平面的应力来构建。
代码片段示意（关键逻辑）：
import torch def calculate_mop_stress(frames, z_planes, box_xy_area): """ frames: 列表，每个元素为包含原子位置、动量、成对力f_ij的快照字典 z_planes: 一维张量，平面z坐标 box_xy_area: 平面面积 Lx * Ly """ num_planes = len(z_planes) stress_kin = torch.zeros((num_planes, 3)) # 动能应力 stress_config = torch.zeros((num_planes, 3)) # 构型应力 for snapshot in frames: pos = snapshot['positions'] # [N, 3] mom = snapshot['momenta'] # [N, 3] mass = snapshot['masses'] # [N] f_ij = snapshot['pair_forces'] # [N, N, 3] 稀疏或稠密矩阵 # 1. 动能贡献 pz = mom[:, 2] flux_kin = mom * pz[:, None] / mass[:, None] # [N, 3] # 将每个原子的贡献分配到最近的平面上（或通过插值） bin_indices = torch.bucketize(pos[:, 2], z_planes) # 简单分桶 for i in range(num_atoms): k = bin_indices[i] if 0 <= k < num_planes: stress_kin[k] += flux_kin[i] # 2. 构型贡献 # 遍历所有相互作用对 (i, j)，这里假设f_ij是稀疏存储 for i, j, f_vec in iterate_sparse_pairs(f_ij): zi, zj = pos[i, 2], pos[j, 2] # 找出所有被跨越的平面 z_min, z_max = min(zi, zj), max(zi, zj) # 找到z_planes中位于(z_min, z_max)区间内的平面索引 mask = (z_planes > z_min) & (z_planes < z_max) crossed_plane_indices = torch.where(mask)[0] sign = 2.0 if zi < zj else -2.0 # 确定d_ij的符号 for k in crossed_plane_indices: stress_config[k] += sign * 0.25 * f_vec # d_ij/4 * f_ij # 平均 num_frames = len(frames) stress_total = (stress_kin + stress_config) / (num_frames * box_xy_area) return stress_total # [num_planes, 3]注意事项：在实际计算中，处理周期性边界条件需要小心。对于跨越周期性边界的原子对，其相对位置矢量 (\mathbf{r}_{ij}) 必须使用最小镜像约定。在判断是否跨越平面时，也需要考虑原子可能因为周期性边界而出现在盒子的“另一侧”。一个稳健的做法是，将所有原子的坐标都还原到其原始映像（unwrapped coordinates）再进行判断，但计算力时使用的仍是基于最小镜像约定的相对位置。
4. 验证案例：氧化锆-水界面体系的应力分析
为了验证平面方法在MACE���函数下的正确性，我们选择一个具有挑战性的体系：氧化锆（ZrO₂）与水的固液界面。这个体系存在强烈的非均匀性、电荷转移和氢键网络，是检验局部应力定义的理想场景。
4.1 模拟体系与设置
我们构建了一个类似于文献��的三斜晶系模拟盒子，包含一个ZrO₂晶体壁面和填充的水分子。为了获得更好的统计并观察体相性质，我们将流体区域在垂直界面的方向进行了扩展。体系在x和y方向采用周期性边界条件，在z方向（垂直于界面）采用固定壁面或可调节的周期性边界。
- 势函数：使用MACE-MP-0预训练模型。该模型在广泛的材料数据库上训练，能较好地处理Zr、O、H等多种元素及其相互作用。
- 模拟细节：采用LAMMPS或ASE等支持MACE的接口进行模拟。首先对体系进行能量最小化，然后在NVT系综下进行充分平衡，使温度和密度达到稳定。对于非平衡模拟（如施加剪切），在平衡后切换到NEMD算法。
- 应力计算：我们同时计算三种应力进行对比：
- 平面方法应力：按照上述流程计算。
- 局部维里应力：将传统的维里公式应用到局部空间格点（bin）中，即 (P^{\text{virial}}(z) = \frac{1}{V_{\text{bin}}} \left\langle \sum_{i \in \text{bin}} \left( \frac{\mathbf{p}i \mathbf{p}i}{m_i} + \frac{1}{2} \sum{j \neq i} \mathbf{f}{ij} \otimes \mathbf{r}_{ij} \right) \right\rangle)。
- 全局维里应力：对整个模拟盒子计算一个平均的维里应力，作为参考。
4.2 结果分析与守恒律验证
下图示意了在平衡状态下，沿z方向（从固体内部穿过界面到液体内部）的法向应力（P_zz）分布。
| z 位置区域 | 平面方法 P_zz | 局部维里 P_zz | 现象与解释 |
|---|---|---|---|
| ZrO₂ 固体内部 | 较大负值（压应力） | 波动剧烈的正值/负值 | 固体内部应力由晶格约束决定。平面方法给出平滑合理的分布。局部维里在原子尺度剧烈震荡，物理意义不明确。 |
| 界面区域 (~1 nm) | 从固体应力平滑过渡，在界面处可能出现峰值（表面张力贡献） | 出现非物理的振荡，甚至符号变化 | 界面处原子密度和成键环境剧烈变化。平面方法正确捕捉了跨越界面的力传递。局部维里因“对半劈分”规则，在化学键跨越bin边界时产生人为噪声。 |
| 液体体相 | 恒定值（体系压力） | 接近但略偏离体系压力，且有微小波动 | 在均匀液体中，两者应趋近于同一恒值。平面方法结果平稳。局部维里结果受bin内原子数统计涨落影响，波动更大。 |
最关键验证：控制体积动量守恒
我们选择一个由两个平行平面z1和z2界定的水层作为控制体积（CV）。根据动量守恒，CV内总动量随时间的变化率，应等于作用在其两个边界平面上的净应力（平面方法应力）：
[ \frac{d}{dt} \sum_{z_1 < z_i < z_2} \mathbf{p}_i(t) = A \left[ \mathbf{P}^{MoP}(z_2, t) - \mathbf{P}^{MoP}(z_1, t) \right] ]
我们在模拟中每隔很短的时间步（如10 fs）就计算一次等式左右两边的值。对于平面方法应力，我们发现在每一个时间步，上述等式都精确成立（仅受浮点数精度限制）。这是一个瞬时守恒律，无需系综平均。这强有力地证明了平面方法应力定义与牛顿运动定律的一致性。
相反，如果我们用局部维里应力（在z1和z2处的bin内平均值）来代替上式右边的 (\mathbf{P}^{MoP})，该等式将不再成立。即使在长时间平均后左右两边可能数值接近，但瞬时值存在显著差异。这说明局部维里应力不能作为描述局部动量通量的精确力学量，在非平衡瞬态模拟中会引入误差。
4.3 在非平衡剪切流中的应用
为了进一步展示平面方法在NEMD中的价值，我们在水层中施加了一个沿x方向的剪切流（例如，通过移动上下壁面或使用Lees-Edwards边界条件）。在稳态剪切下，我们计算了剪切应力P_xz沿z方向的分布。
- 平面方法结果：会显示出一个在液体内部基本恒定、在界面附近发生变化的剪切应力剖面。这与连续介质流体力学在简单剪切流中的预期一致。我们可以从该剖面直接计算流体的表观粘度，并研究界面滑移现象。
- 局部维里结果：在剪切方向均匀的体相区域，其平均值可能与平面方法接近。但在界面附近和速度梯度大的区域，它会再次出现非物理的振荡，使得精确提取壁面剪切应力变得困难，从而影响对滑移长度等关键界面属性的预测。
实操心得：在强剪切或快速瞬变过程中，动量通量的时间相关性很强。平面方法因其固有的守恒性，能够准确追踪动量在空间中的传递，这对于耦合分子动力学与连续计算流体动力学（CFD）的多尺度模拟至关重要。基于平面方法应力提供的边界条件，可以确保微观与宏观模拟之间的动量交换是严格守恒的。
5. 常见问题、挑战与进阶技巧
5.1 计算效率与优化策略
计算所有原子对之间的跨越关系是平面方法中最耗时的部分，复杂度为O(N^2)。对于大规模体系，这可能是瓶颈。以下是一些优化策略：
- 邻居列表与区域划分：利用分子动力学中已有的邻居列表（neighbor list）。对于每个平面，可以预先构建一个“可能跨越”该平面的原子对列表。由于相互作用有截断半径 (r_c)，只有当两个原子在z方向的距离小于 (r_c) 且x、y坐标相近时，才可能跨越某个平面。可以基于此大幅减少需要检查的原子对数量。
- 稀疏矩阵存储：(\mathbf{f}_{ij}) 张量非常稀疏。使用稀疏矩阵格式（如COO、CSR）存储和计算，可以节省大量内存和计算时间。
- 采样频率：应力通常比原子坐标弛豫得更慢。无需每步都计算应力，可以以远大于应力自相关时间的时间间隔进行采样和计算，这能极大降低开销。
- GPU加速：MACE的力计算本身已在GPU上高度优化。应力计算中的循环和判断逻辑也可以编写成GPU核函数，利用其并行能力。例如，可以将“判断原子对是否跨越平面”的任务映射到大量的GPU线程上并行执行。
5.2 处理复杂相互作用与多体势的哲学思考
MACE的高体阶意味着其描述的相互作用远非简单的两体中心力。这引发了一个根本性问题：对于这样的势函数，“原子间作用力” (\mathbf{f}{ij}) 还有明确的物理意义吗？在量子力学中，只有总力 (\mathbf{F}i) 是明确的。我们采用的分解方式 (\mathbf{f}{ij} = \partial U/\partial \mathbf{r}{ij} - \partial U/\partial \mathbf{r}_{ji}) 是一种数学构造，它保证了系统的动量守恒，并且与图网络的边（edge）梯度直接对应。
尽管这种分解在数学上不唯一，但平面方法的美妙之处在于，只要使用满足 (\mathbf{f}{ij} = -\mathbf{f}{ji}) 的分解，并且只统计那些真正跨越平面的相互作用（即原子i和j分居平面两侧），那么最终计算出的平面应力就是唯一定义明确的，并且满足控制体积的动量守恒。这个应力是作用于该平面上的净力的度量，它不依赖于对相互作用路径的想象。因此，即使我们无法描绘MACE中“力”是如何在原子间“传递”的，我们也能可靠地测量它穿过某个平面时产生的宏观效应。
5.3 与其他物理量的耦合：能量流与热通量
应力关注的是动量流。在非平衡模拟中，能量流（热通量）同样重要。基于类似的平面方法思想，可以推导出热通量的表达式。关键区别在于，能量流计算需要用到另一种力分解 (\tilde{\mathbf{f}}_{ij} = -\partial U_j / \partial \mathbf{r}_i)，以及原子自身的能量输运。热通量的平面方法公式会包含动能、势能和对流项的复杂组合。验证其正确性的黄金标准是控制体积内的能量守恒方程。实现能量流计算是未来将ML势函数应用于热输运问题的重要步骤。
5.4 在真实科研工作流中的整合建议
- 先验证，后生产：在将平面方法应力用于重要科学问题之前，务必在简单的均匀体系（如体相水）和已知有解析解或可靠结果的模型体系（如Lennard-Jones流体）中进行验证。确保你的实现能复现正确的体相压力，并在平衡时给出均匀的应力分布。
- 与全局量对比：平面方法应力在均匀区域的空间平均，应与对整个系统计算的全局维里应力（对于均匀体系是有效的）在统计误差内一致。这是一个有用的交叉检查。
- 可视化与诊断：始终绘制应力剖面图。非物理的跳跃、剧烈震荡或趋势错误通常是代码bug或物理定义问题的信号。对于界面体系，检查应力在界面处的行为是否符合物理直觉（如法向应力连续）。
- 误差估计：由于应力计算涉及采样平均，需要报告统计误差。可以使用块平均法（block averaging）或直接计算标准误差。对于非平衡稳态模拟，确保采样时间远大于系统最慢模式的弛豫时间。
将平面方法与MACE等机器学习势函数结合，为我们打开了一扇高精度研究非均匀、非平衡复杂分子体系的大门。它不仅仅是一个新的计算工具，更是确保我们从这些高度复杂的“黑箱”势函数中，提取出物理上可靠、与守恒律自洽的宏观量的关键保障。随着机器学习势函数在电化学、生物物理、极端条件材料等领域日益广泛的应用，掌握这种严格的本征应力分析方法，将成为计算研究者的一项必备技能。
