1. 项目概述当机器学习势函数遇上局部压力计算在分子动力学模拟的世界里压力或应力张量是连接微观原子运动与宏观材料力学性能的桥梁。无论是研究金属的塑性变形、聚合物的粘弹性还是分析血液在微血管中的流动我们最终都需要从原子轨迹中提取出可靠的应力信息。传统上维里应力张量因其计算简便成为了绝大多数模拟软件如LAMMPS、GROMACS的默认选择。然而当系统变得不均匀——例如在固液界面、裂纹尖端、剪切流动前沿或是任何远离热力学平衡的状态时——维里应力的物理基础就开始动摇其计算结果可能严重偏离真实情况。这就像用一把尺子去测量一个扭曲表面的曲率工具本身就不对。对于非平衡分子动力学模拟和流体动力学研究而言一个能在局部严格满足动量守恒的应力定义至关重要。这时“平面方法”就登场了。它的核心思想极其物理直观应力就是力除以面积。MoP通过统计穿过空间中一个假想平面的所有原子间作用力来定义该平面上的应力。这个定义直接源于Irving和Kirkwood的统计力学奠基性工作天生就与动量守恒方程紧密耦合。近年来机器学习势函数的爆发式发展特别是像MACE这样的模型让我们能够以前所未有的精度和速度进行分子模拟。MACE等模型通过图神经网络和消息传递架构巧妙地融合了量子力学的精度与经典分子动力学的效率。但是如何为这些复杂的、非两体的、蕴含量子精度的相互作用力定义一个正确且可计算的局部应力这正是本文要解决的核心问题。我们将深入探讨如何将平面方法拓展到MACE势函数并验证其在非平衡体系中的正确性为复杂界面和流动问题的精确模拟铺平道路。2. 核心原理从统计力学基础到平面方法2.1 传统维里应力的局限与平面方法的物理图像要理解平面方法的优势首先得看清维里应力在非均匀体系中的问题所在。维里应力的常见形式如IK1近似将原子间相互作用对压力的贡献简单地按“一人一半”的方式分配给相互作用的两个原子。在均匀流体中这种平均化处理没有问题。但在界面附近一个原子在液相另一个在固相它们之间的相互作用本质上是跨越界面的“表面力”。将这份力对半劈开分别计入液相和固相的压力物理上是不清晰的会导致界面处的应力分布出现人为的模糊甚至错误。平面方法则回归了应力的力学本源单位面积上的力。想象在流体中放置一个无限薄的平面例如zz0平面。任何试图改变该平面一侧流体动量的企图都必须通过这个平面传递力。因此统计所有在时间Δ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_{i1}^{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-平面所夹的平板区域时侧面对动量的贡献在统计平均下为零假设系统在xy方向均匀或周期性边界动量变化就简化为两个主平面上的应力差[ \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{i1}^{N} \frac{\mathbf{p}i p{i,z}}{m_i} \delta(z_i - z_0) \right\rangle \frac{1}{4A} \left\langle \sum_{i1}^{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等机器学习势函数结合为我们打开了一扇高精度研究非均匀、非平衡复杂分子体系的大门。它不仅仅是一个新的计算工具更是确保我们从这些高度复杂的“黑箱”势函数中提取出物理上可靠、与守恒律自洽的宏观量的关键保障。随着机器学习势函数在电化学、生物物理、极端条件材料等领域日益广泛的应用掌握这种严格的本征应力分析方法将成为计算研究者的一项必备技能。