来源论文: https://arxiv.org/abs/2603.01234v1 生成时间: Mar 03, 2026 08:35

执行摘要

传统的全原子分子动力学（AA MD）模拟在时间和长度尺度上存在根本性限制，难以探究复杂的分子现象。粗粒化（CG）模型通过减少系统自由度来扩展模拟尺度，但现有的基于机器学习（MLP）的CG模型通常受限于训练数据中的噪声，导致精度不足、泛化能力差、超参数选择困难以及缺乏跨热力学状态的迁移性。本文提出了一种创新的解决方案，通过神经演化势能（NEP）框架下的NEP-CG和NEP-AACG方法，克服了这些挑战。

该研究的核心创新在于生成低噪声训练数据。通过在AA MD模拟中约束CG珠子位置，并累积时间平均力与维里，可以直接获得对应平均力势（PMF）的平滑且真实的训练数据。这种方法与传统拟合瞬时力的力匹配方法形成鲜明对比，后者容易受到噪声干扰。此外，研究引入了维里校正，补偿了粗粒化过程中损失的自由度，对于精确重现状态方程至关重要。

在NEP框架下实现的NEP-CG模型，针对液态水，实现了与基于密度泛函理论（DFT）数据训练的AA模型相当的精度。它能准确重现从1巴到1吉帕的密度，并在0.5吉帕的训练限制外成功外推。对于各向异性C60单分子层，通过区分晶体学上不同的珠子类型，应力误差降低了一个数量级，并成功捕获了方向性热导率。进一步，多尺度NEP-AACG模型无缝整合了AA和CG自由度，并以金纳米线断裂模拟为例，证明了其在实验相关应变速率下进行多尺度模拟的能力。

在计算性能方面，NEP-CG模型的速度达到每天数百至数千纳秒，仅需单个消费级GPU。这比相应的AA模型实现了显著的加速（液态水提升约50倍，C60单分子层提升约1000倍），使得大规模、长时间的模拟成为可能。

这项工作提供了一个强大而通用的框架，用于构建跨多样化系统的高精度、可迁移且高效的CG模型，为化学、生物学和材料科学领域的分子模拟开辟了新途径。

核心科学问题，理论基础，技术难点，方法细节

1. 核心科学问题

分子动力学（MD）模拟是理解复杂系统微观行为的强大工具，但其应用受到计算能力和时间/长度尺度的限制。许多重要的分子现象，如蛋白质折叠、材料断裂、相变等，发生在微秒甚至更长时间尺度和纳米甚至微米空间尺度上，远超传统全原子（AA）MD模拟可达的范围。为了解决这一根本性挑战，粗粒化（CG）建模应运而生，它通过将一组原子抽象为少数几个相互作用的“珠子”，有效降低了系统的自由度，从而能够模拟更大系统和更长时间尺度。

然而，现有CG模型，特别是近年来兴起的基于机器学习（MLP）的CG模型，在实际应用中仍面临诸多瓶颈：

  1. 训练数据噪声：MLP-CG模型通常依赖于原子轨迹中瞬时力进行训练，这些瞬时力包含大量热力学噪声。这种噪声导致训练误差显著，使得模型收敛性难以评估，超参数选择变得复杂。
  2. 精度与泛化能力不足：高噪声训练数据限制了模型的内在精度，导致其在预测结构、热力学和动力学性质时不如AA模型准确。此外，模型往往仅适用于单一热力学状态（如特定密度或压力），缺乏跨不同压强或温度条件的迁移性。
  3. 数据效率低下：为了克服噪声影响，传统MLP-CG方法需要大量的训练结构，即使对于单一状态点也可能需要数千个快照，这增加了数据采集和计算成本。
  4. 多尺度模拟的复杂性：在需要同时处理原子级细节和粗粒化区域的系统中（如裂纹扩展、界面现象），如何无缝地集成AA和CG自由度，并确保界面区域的物理一致性，是一个巨大的挑战。

本文的工作旨在通过创新的训练数据生成策略，显著提高MLP-CG模型的精度、可迁移性和计算效率，并提供一个统一的多尺度模拟框架。

2. 理论基础

2.1 分子动力学与粗粒化

分子动力学（MD）通过求解牛顿运动方程，模拟原子或分子随时间的演化。在AA MD中，每个原子都作为一个独立的粒子处理，其相互作用由势能函数描述。粗粒化（CG）建模将一组AA原子映射到一个CG珠子，其位置通常由这些原子的质心决定。CG模型的势能函数，也称为有效势能，旨在重现原始AA系统的关键结构和热力学性质。

2.2 平均力势（Potential of Mean Force, PMF）

在统计力学中，粗粒化系统的有效势能可以通过平均力势（PMF）来定义。PMF UCG(R) 描述了在固定CG构型R下，所有被积分掉的原子自由度上的自由能。其数学表达式为：

$UCG(R) = -kBT \ln [ \int dr \delta(R - C \cdot r) e^{-\beta U(r)} ]$

其中，k_B是玻尔兹曼常数，T是温度，U(r)是AA系统的势能，C是将AA原子坐标r映射到CG珠子坐标R的线性映射。PMF是温度和压力的函数，因为它反映了被积分掉的原子自由度的热力学平均效应。

2.3 力匹配方法

力匹配方法（Force-Matching Method, FMM）是优化CG模型参数的常用策略。其核心思想是使CG珠子上的力尽可能接近相应AA原子组上的系综平均力。具体而言，CG珠子I上的平均力 $F_I^{CG}(R)$ 定义为：

$F_I^{CG}(R) = -\nabla_{R_I} UCG(R) = \langle \sum_{i \in I} f_i \rangle_R$

其中，$f_i$是原子i上的力，$\langle \dots \rangle_R$表示在CG构型R受限条件下的系综平均。传统力匹配方法通常通过最小化CG模型预测力与瞬时AA力之间的均方误差（RMSE）来训练CG势能：

$\chi^2 = \frac{1}{3 N_{CG}} \sum_{I=1}^{N_{CG}} \langle |f_I^{CG}(r) - F_I^{CG}(C \cdot r)|^2 \rangle$

然而，这种基于瞬时力的训练方式，尤其是对灵活的MLP，容易导致模型过拟合训练数据中的噪声，从而影响模型的准确性和泛化能力。

2.4 机器学习势能（MLP）与神经演化势能（NEP）

机器学习势能（MLP）利用神经网络或其他机器学习算法来拟合原子相互作用，能够捕获复杂的多体效应。与传统的经验势能不同，MLP可以从量子力学（QM）计算中学习，达到高精度。

神经演化势能（NEP）是GPUMD包中实现的一种高效MLP架构。NEP通过前馈神经网络表示每个原子的局部势能，势能是基于原子局部环境描述符向量q的函数。描述符包括径向和角向分量，它们由球谐函数、切比雪夫多项式和余弦截止函数构建。NEP的参数通过可分离自然演化策略（SNES）优化，损失函数包括能量、力、维里的RMSE以及L1/L2正则化项。

3. 技术难点

现有的MLP-CG方法面临的主要技术挑战在于：

  1. 训练数据噪声过大：瞬时力匹配导致的高RMSE（例如，水模型达到3.49到11.8 kcal/mol/Å，远高于AA模型的~0.05 eV/Å），这使得模型难以准确学习底层的PMF。
  2. 收敛性评估困难：高噪声训练数据掩盖了模型收敛的真实情况，难以判断模型何时达到最佳性能。
  3. 超参数选择复杂：噪声使得超参数（如神经网络大小、截止半径、正则化强度）的优化变得繁琐，因为训练误差不能清晰地反映模型性能。
  4. 数据效率低：为了克服噪声，模型需要大量训练数据，增加了计算成本和时间。
  5. 模型迁移性差：传统MLP-CG模型通常只在单一热力学状态下训练，在其他压力或温度条件下性能下降，缺乏迁移性。
  6. 维里（压力）的准确性：CG模型在 coarse-graining 过程中会丢失部分自由度，导致维里（及压力）的计算不准确，影响状态方程的重现。
  7. 各向异性系统处理：对于具有方向性键合的各向异性系统，如C60单分子层，简单的各向同性CG珠子难以捕获其复杂的力学和热学性质。
  8. 多尺度模拟的统一框架：在AA和CG区域之间进行无缝过渡，同时保持物理一致性，以及处理固定分辨率限制等问题。

4. 方法细节

本文提出的NEP-CG和NEP-AACG方法通过以下创新性细节解决了上述挑战：

4.1 低噪声训练数据生成：基于PMF的约束MD模拟

这是本文的核心创新。为了克服瞬时力噪声，研究团队提出了一种生成低噪声训练数据的方法，该数据直接对应于PMF定义的系综平均力：

  1. AA系统平衡：首先，在目标温度和压力/密度下，使用AA NEP模型对原子系统进行热平衡。
  2. CG珠子约束：在平衡后，对CG珠子的位置施加约束（例如，将水分子中三个原子的质心固定），然后继续进行NVE系综模拟。
  3. 累积时间平均力与维里：在NVE模拟过程中，累积每个CG珠子上的瞬时力 $f_i$ 和维里 $W_i$。通过长时间（例如1-10纳秒）的模拟，这些瞬时量的时间平均值将收敛到真实的系综平均值（对应于PMF），噪声水平远低于瞬时力，从而易于MLP学习。

这种方法生成的训练数据本质上更平滑、更准确，解决了传统力匹配方法中噪声过大的问题。

4.2 维里校正（Virial Correction）

为了确保CG模型能准确重现状态方程和跨压力条件的迁移性，研究引入了维里校正。在粗粒化过程中，由于被积分掉的原子自由度（动能部分）被移除，CG模型会系统性地低估压力或高估密度。为此，瞬时维里被校正为：

$W \rightarrow W + (N_{AA} – N_{CG}) k_B T I$

其中，$N_{AA}$是AA原子数量，$N_{CG}$是CG珠子数量，$k_B$是玻尔兹曼常数，$T$是温度，$I$是单位张量。这个校正项补偿了损失的自由度，确保了正确的压力响应。

4.3 NEP-CG模型

NEP-CG模型旨在纯CG模拟中重现系综平均量。NEP架构保持不变，但站点能量 $U_i$ 现在表示CG珠子的自由能贡献。训练数据包括珠子力（经过时间平均）和维里张量（经过时间平均和校正）。最小化NEP损失函数（包含力、维里和能量的RMSE）使得模型能够以与AA模型相当的精度进行训练，同时具有更高的数据效率。

4.4 NEP-AACG模型

NEP-AACG（全原子-粗粒化）模型旨在统一多尺度模拟，将AA和CG自由度集成在一个模型中：

  1. 混合分辨率处理：NEP架构可以自然处理混合分辨率系统，将AA原子和CG珠子视为不同的“物种”。
  2. 统一描述符框架：描述符从局部环境中构建，其中可能包含AA原子和CG珠子，使用相同的径向和角向框架。不同粒子对可能需要不同的截止半径，以优化局部环境的捕获。
  3. 训练数据生成：在生成训练数据时，对于每个构型，所有站点（AA原子和CG珠子）都被约束，并累积每个站点的系综平均量。NEP-AACG模型通过同时预测所有量，学习跨两种分辨率的一致自由能曲面。
  4. 无缝过渡：AA和CG区域之间的耦合自然地从训练数据中产生，无需额外的ad hoc方案。统一的NEP框架支持区域之间的平滑过渡和空间分辨率的变化，从而实现无界面伪影的多尺度模拟。

4.5 NEP架构细节

  • 神经网络结构：使用前馈神经网络，一个隐藏层，tanh激活函数。站点能量 $U_i$ 是描述符向量q的函数。
  • 描述符：径向描述符 $G_n(r_{ij}) = \sum_{j\ne i} g_n(r_{ij})$ 和角向描述符 $A_n^{lm} = \sum_{j\ne i} g_n(r_{ij}) Y_{lm}(\theta_{ij}, \phi_{ij})$。其中 $g_n(r_{ij})$ 使用切比雪夫多项式和余弦截止函数 $f_c(r_{ij})$ 进行展开。
  • 参数优化：采用可分离自然演化策略（SNES）。损失函数是能量、力、维里的均方误差（RMSE）与L1和L2正则化项的加权组合，所有权重均为可调超参数。

这些方法细节共同构成了一个强大、灵活且高效的框架，能够为各种复杂系统构建高精度CG和多尺度模型。

关键 benchmark 体系，计算所得数据，性能数据

本研究通过三个具有代表性的基准体系，全面展示了NEP-CG和NEP-AACG方法的强大能力：液态水、C60单分子层网络和金。每个体系都旨在突出方法在精度、迁移性、各向异性处理和多尺度集成方面的关键优势。

1. 液态水 NEP-CG 模型

液态水是传统CG和MLP-CG方法中广泛研究的系统，是验证新方法的理想基准。

1.1 训练数据生成与模型超参数

  • AA参考模型：使用基于CCSD(T)级别MB-pol参考数据训练的NEP-AA模型作为基准，该模型能够准确预测水的结构、热力学和输运性质。
  • 训练数据生成：
    1. 在300 K下，AA NEP模型在NPT系综下平衡原子系统，压力范围从1巴到0.5吉帕（分五个离散点）。
    2. 对每个压力点，进行1纳秒NPT平衡。
    3. 关闭恒温器和恒压器，对CG珠子（每个CG珠子对应一个水分子的质心）施加约束。
    4. 进行10纳秒NVE模拟，累积时间平均力与维里数据。原子系统包含1536个原子，映射为512个CG珠子。
    5. 总训练数据集仅包含5个结构，总计2560个水珠，验证了方法的高数据效率。
  • 模型超参数：
    • 隐藏层神经元数量 $N_{neu}$：NEP-AA模型为60，NEP-CG模型经过广泛测试后确定为10。这证实了CG模型的势能面更加平滑，需要的参数更少。
    • 截止半径：径向截止6 Å，角向截止4 Å，与NEP-AA模型相同，因为水分子质心与氧原子大小相当。

1.2 模型精度与性能

  • 力与应力奇偶性图：

    • 传统瞬时力方法（在1巴压力下训练）：力与应力奇偶性图显著偏离y=x线，斜率被低估，存在明显散射。RMSE较高（力0.15 eV/Å，应力0.14 GPa）。这表明瞬时力不是正确的训练目标，真实平均力被噪声淹没。
    • 本文系综平均力方法（跨1巴到0.5吉帕压力采样数据）：奇偶性图紧密跟随y=x线。RMSE显著降低（力0.080 eV/Å，应力0.0084 GPa）。这证明了系综平均力是正确的训练目标，且模型在不同压力条件下具有出色的迁移性。

  • 结构性质与状态方程：

    • 径向分布函数（RDF）：NEP-CG模型在1巴和0.5吉帕压力下计算的RDFs与NEP-AA模型高度一致，再次验证了模型在不同压力条件下的结构迁移性。
    • 密度-压力曲线：NEP-CG模型准确重现了NEP-AA模型从1巴到1吉帕的密度（密度从约1 g/cm³增加到1.2 g/cm³）。尽管训练数据仅扩展到0.5吉帕，模型仍能成功外推到1吉帕，凸显了系综平均方法的鲁棒性。
    • 维里校正的重要性：未经维里校正的NEP-CG模型会系统性地高估密度（如图2b所示），证实了维里校正对于准确重现状态方程的不可或缺性。

2. 富勒烯单分子层网络 NEP-CG 模型

第二个例子应用于共价键合C60分子的准六方相（QHP）单分子层，该体系具有各向异性的力学和热学性质，是验证模型处理复杂键合和各向异性能力的关键。

2.1 体系特性与珠子类型区分

  • 体系特性：QHP-C60单分子层具有各向异性键合：每个C60分子与六个邻居共价连接，键合类型因方向而异（例如，[2+2]环加成键沿y轴，C-C单键沿[110]和[110]方向）。
  • 珠子类型区分：为了捕获这种各向异性，研究引入了两种不同的CG珠子类型，对应于晶胞中两种C60分子。单一珠子类型无法区分x和y方向，无法捕获各向异性。

2.2 训练数据生成与模型超参数

  • AA参考模型：使用最近开发的用于研究QHP-C60单分子层热输运的NEP-AA模型作为基准。
  • 训练数据生成：
    1. 使用NEP-AA模型对7200原子系统（映射为120个CG珠子）进行NPT模拟，温度300 K，面内压力范围从1巴到1吉帕（两个方向）。
    2. 1纳秒平衡后，约束C60质心位置，进行10纳秒生产运行，累积力与维里数据。
    3. 总共生成了23个训练结构，捕获了不同面内应变状态下的各向异性响应。
  • 模型超参数：
    • 截止半径：径向和角向描述符的最佳截止半径为25 Å。由于C60分子尺寸大且分子间相互作用范围广，需要较大的截止半径来准确捕获相互作用和各向异性键合环境。
    • 隐藏层神经元数量 $N_{neu}$：CG模型使用5个神经元即可实现高精度训练，因为粗粒化过程平滑了势能面，所需的神经元数量显著减少。

2.3 模型精度与热输运

  • 力与应力奇偶性图：

    • 单一珠子类型模型：表现出系统性误差（力RMSE 0.11 eV/Å，应力RMSE 0.083 GPa），未能捕获共价键合的方向依赖性。
    • 两种珠子类型模型：显著改善了模型与参考数据的一致性。力RMSE降至0.090 eV/Å，应力RMSE更是降低了一个数量级，达到0.0025 GPa。这强调了在粗粒化各向异性系统时，基于晶体学等价性引入不同珠子类型的重要性。

  • 热输运：

    • 使用均质非平衡分子动力学（HNEMD）方法计算热导率。单分子层的有效厚度为8.785 Å。
    • NEP-CG模型成功重现了QHP-C60单分子层的各向异性热传导特性（y方向热导率大于x方向）。
    • 通过60倍的比例因子（由于自由度减少），NEP-CG的热导率（κx = 64 ± 12 W/mK，κy = 95 ± 9 W/mK）与NEP-AA参考值（κx = 102 ± 3 W/mK，κy = 137 ± 7 W/mK）在数量级上一致。这表明CG模型在显著降低自由度的同时，仍然捕获了定性各向异性和近似的热导率。

3. 金 NEP-AACG 模型

第三个案例研究展示了NEP-AACG框架在金系统中的应用，通过构建混合分辨率模型来研究金纳米线在拉伸下的变形和断裂，体现了多尺度模拟的优势。

3.1 训练数据生成与模型超参数

  • AA参考模型：使用UENP-v1势能（涵盖16种金属及其合金，性能优于传统嵌入原子势）作为AA参考模型。
  • 训练数据生成：
    1. 在300 K下，使用NEP-AA模型进行NVT模拟，涵盖多种应变状态：各向同性应变(-10%至5%)，双轴应变(-10%至5%)，单轴应变(-10%至10%)，剪切应变(-5%至5%)。
    2. 模拟使用2048个金原子，5飞秒时间步长。1纳秒NVT平衡后，进行10纳秒约束模拟以累积参考数据。
    3. 对于每个应变构型，设计两种映射方案：
      • 纯CG映射：每个CG珠子代表4个金原子（一个FCC晶胞），将2048个原子减少到512个珠子。
      • 混合分辨率AACG映射：一半原子进行CG映射，另一半保持AA分辨率。
  • 模型超参数：
    • 粒子类型：定义两种粒子类型：AA金原子和CG珠子，它们具有不同的有效尺寸。
    • 截止半径：针对不同粒子类型使用不同的截止半径。
      • UENP-v1模型：径向6 Å，角向5 Å。
      • CG珠子：径向8 Å，角向7 Å。
      • 混合AA-CG相互作用：径向7 Å，角向6 Å（取纯类型截止的算术平均值）。
    • 隐藏层神经元数量 $N_{neu}$：选择中间值30，以平衡模型容量和计算效率。

3.2 模型精度与力学验证

  • 力与应力奇偶性图：

    • 统一的NEP-AACG模型在所有三个数据集（纯AA、纯CG、混合AACG）上都表现出出色的精度。
    • 纯AA数据集：力RMSE 0.10 eV/Å，应力RMSE 0.50 GPa，与原始UENP-v1模型相当，证实添加CG和混合分辨率数据不会损害原子级精度。
    • 纯CG和混合AACG数据集：力RMSE 0.077 eV/Å，应力RMSE 0.086 GPa。较低的应力误差反映了更平滑的CG势能面和CG相关数据集中更窄的力与应力范围。
    • 关键结论：所有预测均由同一个统一的NEP-AACG模型生成，具有两种粒子类型和相互作用特异性截止半径。模型同时描述原子级、粗粒化和混合分辨率构型的能力，表明它学习到了跨多分辨率级别的一致自由能曲面，实现了无缝多尺度模拟。

  • 本体力学验证：

    • 体系：沿x方向进行单轴拉伸变形的体心立方金，工程应变 $\epsilon_x$ 高达5%。
    • 模拟：纯AA模拟（16384个原子）和纯CG模拟（4096个珠子）。
    • 结果：轴向应力 $\sigma_x$ 在弹性区与 $\epsilon_x$ 呈线性关系。横向应变 $\epsilon_y$ 表现出预期的泊松收缩（泊松比约为0.4）。NEP-AACG模型的结果在整个应变范围内与NEP-AA参考模型定量一致。这证实了CG模型（通过应变构型的系综平均数据训练）能忠实重现本体力学响应。

3.3 多尺度纳米线断裂模拟

  • 体系：构建了一个总长约80纳米的金纳米线拉伸测试模型，包括几个区域：
    • 中心原子级纳米线区域：约37纳米长，3纳米厚，预计发生断裂。
    • 两个末端CG区域：延伸至边界，作为施加真实机械约束的AA区域的储层。
    • 两个AA过渡区域：连接AA和CG区域。
  • 模拟：以 $10^7 s^{-1}$ 的恒定工程应变速率施加单轴拉伸，该速率在高速加载实验可达到的范围内。
  • 结果：
    • 应力-应变响应：在0.01到0.05的应变范围内，最大应力波动约30兆帕，纳米线保持结构完整性。应变超过0.07后，纳米线开始在中间变细，应力降至10兆帕以下。应变超过0.1时，纳米线完全断裂。
    • 结论：NEP-AACG模型在纳米线断裂模拟中成功展示了多尺度模拟的能力，能够在断裂区域保持原子级分辨率，同时在较远区域使用CG模型，有效避免了边界条件伪影，并实现了对材料变形复杂过程的研究。

4. 计算性能提升

粗粒化模型的关键优势在于其计算效率，能够访问更长的时间和空间尺度。研究通过两个指标对NEP-CG模型与NEP-AA模型进行了基准测试：粒子吞吐量（粒子-步长/秒）和等效空间尺寸下的有效模拟速度（纳秒/天）。

  • 液态水：

    • 吞吐量：NEP-CG模型具有显著更高的粒子吞吐量，因为它采用更小的神经网络（$N_{neu}$=10 vs. 60）和更少的邻居数量。
    • 时间步长：NEP-AA模型使用0.5 fs，而NEP-CG模型可以安全地使用2 fs，提高了4倍。
    • 粒子数量减少：从原子到珠子3:1的减少。
    • 总加速：综合这些因素，NEP-CG模型在等效空间尺寸下实现了约50倍的模拟速度提升（纳秒/天）。

  • C60单分子层：

    • 吞吐量：NEP-CG模型（$N_{neu}$=5 vs. 50）具有极高的吞吐量。尽管截止半径较大（25 Å），但由于相互作用站点减少，每个珠子的邻居数量仍然较少。
    • 时间步长：NEP-AA模型使用1 fs，而NEP-CG模型允许使用20 fs，提高了20倍。
    • 粒子数量减少：从原子到珠子60:1的减少。
    • 总加速：综合这些因素，NEP-CG模型在等效空间尺寸下实现了约1000倍的惊人模拟速度提升（纳秒/天）。

  • 整体性能：NEP-CG模型的计算速度在单个消费级GPU（NVIDIA RTX 5090）上达到每天数百至数千纳秒。这种显著的性能提升使得以前无法实现的大规模、长时间模拟成为可能。

代码实现细节，复现指南，所用的软件包及开源 repo link

本研究提出的NEP-CG和NEP-AACG方法是在GPUMD分子动力学模拟包中实现的，利用其高效的神经演化势能（NEP）框架。以下是代码实现的关键细节、复现指南以及相关开源资源。

1. 所用软件包与开源资源

  • GPUMD: GPUMD (GPU-accelerated Molecular Dynamics) 是一个高性能分子动力学模拟软件包，特别适用于使用机器学习势能进行模拟。NEP框架，特别是第四代NEP (NEP4)，已集成在GPUMD中，能够高效处理多组分系统。本研究中所有计算均使用GPUMD-v4.9.1及其未来版本完成。
    • GPUMD GitHub 仓库：https://github.com/brucefan1983/GPUMD
  • NEP 数据仓库: 本工作生成的所有训练数据集和训练好的机器学习势能模型均在nep-data仓库中免费提供。
    • NEP Data GitLab 仓库：https://gitlab.com/brucefan1983/nep-data

2. NEP框架实现细节

2.1 势能表示

NEP模型使用一个前馈神经网络来表示原子的局部势能 $U_i$。该网络包含一个单隐藏层，使用tanh激活函数。

  • 描述符：每个原子 $i$ 的局部环境通过一个描述符向量 $\mathbf{q}$ 来表示，该向量由径向和角向分量构成。
    • 径向分量：由径向函数 $g_n(r_{ij})$ 的和构成，用于描述原子间距离。
    • 角向分量：包括三体、四体和五体项，由径向函数 $g_n(r_{ij})$ 和球谐函数 $Y_{lm}(\theta_{ij}, \phi_{ij})$ 构成，用于描述键角。
  • 径向函数展开：径向函数 $g_n(r_{ij})$ 本身是通过一系列切比雪夫多项式 $T_k$ 和余弦截止函数 $f_c(r_{ij})$ 进行展开的，确保势能在截止距离处平滑衰减到零。
  • 原子类型编码：原子类型信息直接编码在展开系数 $c_{nk}$ 中，从而实现了对多组分系统的高效处理，无需独立的描述符分支。

2.2 参数优化

NEP的参数（连接权重和偏置）通过**可分离自然演化策略（SNES）**进行优化。SNES是一种有效的黑盒优化算法，适用于高维非凸优化问题。

  • 损失函数：优化目标是最小化一个综合损失函数，该函数包含以下几项的均方误差（RMSE）：
    • 能量RMSE：原子总能量的预测误差。
    • 力RMSE：原子力的预测误差。
    • 维里RMSE：系统维里张量的预测误差，对于准确捕捉压力响应至关重要。
  • 正则化：为了防止过拟合，损失函数中还包含了L1和L2正则化项，并由可调超参数 $\lambda_e, \lambda_f, \lambda_v, \lambda_1, \lambda_2$ 加权。

3. 复现指南

复现本研究结果主要涉及以下步骤：GPUMD软件的安装、训练数据的准备、NEP模型的训练以及MD模拟的执行和结果分析。

3.1 准备工作：安装GPUMD

  1. 访问GitHub仓库：首先，从GPUMD的官方GitHub仓库 (https://github.com/brucefan1983/GPUMD) 获取最新版本的源代码。
  2. 编译与安装：遵循GPUMD仓库中的编译指南进行安装。通常，这涉及使用Cmake构建系统，并确保CUDA工具包已正确配置，以便利用GPU加速。

3.2 训练数据准备

本研究的核心在于生成低噪声的训练数据。有两种主要方式获取训练数据：

  1. 下载预生成数据：从nep-data GitLab仓库 (https://gitlab.com/brucefan1983/nep-data) 下载本研究中使用的所有预生成训练数据集。这些数据集已经包含了经过时间平均和维里校正的力、维里和能量信息。
  2. 自行生成训练数据：如果需要针对新体系或新条件训练模型，可以按照论文第二部分“II. METHODS”中描述的步骤自行生成数据：
    • AA平衡模拟：使用一个已有的高精度AA NEP模型（如针对水、C60或金的参考模型）在目标温度和压力/密度下进行NPT系综MD模拟，使系统达到平衡。
    • CG珠子映射与约束：定义CG珠子的映射规则（例如，水分子质心、C60质心或金的FCC晶胞质心），然后对这些CG珠子的位置施加几何约束。
    • NVE生产模拟与数据累积：切换到NVE系综，并进行足够长时间（例如1-10纳秒）的模拟，在此期间，累积每个CG珠子上的瞬时力、维里，并计算它们的时间平均值。
    • 维里校正：根据公式 $(N_{AA} – N_{CG}) k_B T I$ 对累积的维里进行校正。

3.3 NEP模型训练

GPUMD提供了训练NEP模型的工具。训练过程主要涉及配置训练参数和指定训练数据。

  1. 配置文件：准备一个训练配置文件（通常为 .in 或 .json 格式），指定神经网络的架构（隐藏层神经元数量 $N_{neu}$）、描述符参数（截止半径、径向/角向分量数量）、优化器设置（SNES参数）、损失函数权重（能量、力、维里、L1/L2正则化项的权重）以及训练数据的路径。
  2. 多组分系统：对于NEP-AACG模型，需要在配置文件中定义不同的粒子类型（例如AA金和CG金珠），并为不同粒子对指定不同的截止半径，以优化混合环境下的相互作用。
  3. 执行训练：使用GPUMD的训练模块执行训练命令。训练过程将根据配置文件中的设置，迭代优化NEP模型的参数，直至收敛。

3.4 分子动力学模拟

训练好的NEP模型可用于执行大规模的MD模拟。

  1. 模型部署：将训练好的NEP模型文件（通常为 .nep 格式）放置在GPUMD模拟可访问的路径。
  2. 模拟配置文件：准备一个模拟配置文件，指定系统初始构型、温度、压力、系综类型（NVT, NPT, NVE）、时间步长、以及加载的NEP模型文件。
  3. 时间步长调整：CG模型由于其势能面更平滑，通常可以采用比AA模型更大的时间步长，从而显著提高模拟效率。
  4. 多尺度模拟：对于NEP-AACG模型，初始构型将包含AA和CG区域。GPUMD将根据模型中定义的粒子类型，自动应用AA和CG相互作用。

3.5 结果分析

模拟结束后，可以对输出的轨迹和数据进行分析，以验证模型的性能：

  • 结构性质：计算径向分布函数（RDFs）等，与AA参考结果进行比较。
  • 热力学性质：计算密度-压力曲线、热导率（通过HNEMD等方法）等，与参考数据进行对比。
  • 力学性质：进行应力-应变分析，评估模型的力学响应。
  • 性能评估：通过粒子吞吐量和有效模拟速度（纳秒/天）来量化计算性能提升。

通过遵循这些指南，研究人员可以复现本研究的结果，并利用NEP-CG和NEP-AACG框架开发新的粗粒化模型。

关键引用文献，以及你对这项工作局限性的评论

1. 关键引用文献

本文的工作建立在分子动力学、粗粒化方法和机器学习势能领域的深厚基础上，并在此基础上进行了创新。以下是与本文密切相关的关键引用文献类别及其代表性工作：

1.1 力匹配方法（Force-Matching Method, FMM）

  • Izvekov and Voth (2005)²：首次提出将力匹配方法应用于粗粒化模型，为从原子轨迹中提取有效相互作用提供了途径。
  • Noid et al. (2008)³：为力匹配方法奠定了严格的统计力学基础，将CG模型的有效势能与平均力势（PMF）联系起来，强调了系综平均力的重要性。

1.2 机器学习粗粒化（MLP-CG）

  • John and Csanyi (2017)¹³：基于高斯近似势（GAP）框架开发CG模型，展示了其在无溶剂生物分子系统中比对势模型更高的精度和更快的模拟速度。
  • Zhang et al. (2018)¹⁵：使用力匹配方法和深度势能（DP）方法构建DP-CG水模型，精确重现了氧原子的结构特征。
  • Wang et al. (2019)¹⁷：引入CGnets，表明MLP可以仅用少量CG珠子捕获显式溶剂自由能面。
  • Loose et al. (2023)²³：讨论了使用等变神经网络进行粗粒化，并报告了传统MLP-CG方法在训练数据中存在的较高RMSE值，这与本文所指出的噪声问题相符。

1.3 神经演化势能（NEP）框架

  • Fan et al. (2021)²⁴：提出了NEP框架，一种结合高精度和低计算成本的机器学习势能，并应用于原子模拟和热输运。
  • Xu et al. (2025)²⁵：GPUMD 4.0的介绍，其中NEP框架得到了高效实现，并广泛应用于各种材料模拟。
  • Ying et al. (2025)²⁶：总结了NEP在复杂材料建模中的进展，突出了其在解决材料科学挑战方面的潜力。
  • Fan et al. (2022)³¹：GPUMD软件包的详细介绍，包括NEP的描述符和网络架构细节，为本文的实现提供了基础。

1.4 各向异性粗粒化模型

  • Argun and Statt (2025)²⁷：研究了用于高效模拟各向异性胶体的机器学习势能，探讨了各向异性粒子在CG模型中的重要性。
  • Nguyen and Huang (2022)³⁸ 和 Wilson and Huang (2023)³⁹：探讨了通过力匹配结合扭矩匹配，使用各向异性粒子进行系统性底向上分子粗粒化，以捕获定向自由度。
  • Campos-Villalobos et al. (2024)⁴³：机器学习的各向异性形状和相互作用粒子粗粒化势能，进一步强调了处理各向异性在CG建模中的挑战。

2. 对这项工作局限性的评论

本文提出的NEP-CG和NEP-AACG方法在提升粗粒化模型精度和效率方面取得了显著进展，但论文作者也坦诚地指出了当前工作的一些局限性，这些局限性也为未来的研究指明了方向。以下是对这些局限性以及其他潜在限制的评论：

2.1 论文明确指出的局限性

  1. 温度依赖性：

    • 局限性：所有训练数据均在单一温度（300 K）下生成。平均力势本质上是温度依赖的，这意味着模型可能无法在显著不同温度下保持其高精度和迁移性。
    • 评论：这是一个普遍存在于许多CG模型中的问题。虽然本文模型在不同压力下表现出良好的迁移性，但跨温度的迁移性更为复杂。未来的工作需要开发能够将多个温度条件整合到训练数据中的策略，例如结合温度传输性CG（temperature-transferable CG）模型或使用多温度系综平均。这将使得模型能够用于模拟相变或其他涉及温度变化的现象。

  2. 各向同性珠子：

    • 局限性：当前模型仅使用各向同性珠子。虽然对于C60单分子层，通过区分两种晶体学上不同的珠子类型成功捕获了各向异性，但这种策略无法代表需要定向自由度的复杂系统，如聚合物、液晶或更通用的各向异性分子。
    • 评论：解决这一问题需要扩展NEP框架以适应各向异性粒子表示，例如引入方向描述符或旋转自由度。这在理论和实现上都更具挑战性，但对于拓宽NEP-CG/AACG在软物质和生物分子系统中的应用至关重要。

  3. 固定分辨率划分：

    • 局限性：NEP-AACG模型中AA和CG区域的划分是预定义且固定的，缺乏动态分辨率适应能力。这意味着在模拟过程中，如果最初粗粒化的区域发生了需要原子级细节的关键事件（例如化学反应、断裂开始），模型无法动态调整其分辨率。
    • 评论：动态分辨率适应（dynamic resolution adaptation）是一个前沿且困难的研究方向，通常涉及主动学习或基于局部物理量（如密度、能量梯度）的自动分辨率切换机制。如果能实现，将极大扩展多尺度模拟的应用范围，使其能够处理更加复杂的、演化中的系统。

2.2 其他潜在局限性和未来研究方向

  1. 映射方案的普适性：

    • 评论：本文主要采用质心映射（center of mass mapping）。虽然这对于水和C60等相对刚性的分子是有效的，但对于高度柔性的分子或聚合链，质心映射可能不足以捕获所有重要的构象信息。未来可以探索更复杂的映射方案，如虚拟站点（virtual sites）或模糊映射（fuzzy mapping），以提高柔性体系的粗粒化精度。

  2. AA-CG界面处理：

    • 评论：虽然NEP-AACG通过统一框架实现了AA和CG的无缝集成，但在AA-CG界面区域，仍然可能存在一些由于分辨率差异导致的微小伪影。虽然本文在纳米线断裂模拟中没有明显提及，但在某些高度敏感的界面现象中，这可能需要更精细的界面力校正或区域划分策略。

  3. 多组分系统的复杂性：

    • 评论：NEP4声称对多组分系统高效，本文示例包括水（一种CG珠子）、C60（两种CG珠子）和金（AA/CG金珠）。未来更复杂的、化学性质迥异的多组分系统（如聚合物混合物、生物大分子与溶剂相互作用）将是更严格的考验，以验证NEP在处理复杂化学环境方面的普适性。

  4. 长程相互作用：

    • 评论：机器学习势能通常是短程的，依赖于截断半径内的局部环境。对于涉及长程静电或范德华力的系统，如果不能很好地在截断半径内捕获这些效应，可能会引入误差。虽然本文的截断半径对于C60较大，但在更广泛的生物分子体系中，长程相互作用的处理仍然是一个挑战，可能需要结合长程校正或专门的MLP设计。

  5. 不确定性量化：

    • 评论：当前MLP模型通常提供确定性预测。对于粗粒化模型，量化预测的不确定性对于评估模型在未见过区域的可靠性至关重要。未来可以将不确定性量化方法（如高斯过程或贝叶斯神经网络）集成到NEP框架中，以提供更鲁棒的预测和指导主动学习。

尽管存在这些局限性，本文的工作通过引入低噪声训练数据和维里校正，为高精度、可迁移的MLP-CG模型设定了新的标准，是分子模拟领域的重要一步。

其他你认为必要的补充

1. 低噪声训练数据：粗粒化建模的范式转变

本文最核心且最具影响力的创新，无疑是通过约束分子动力学模拟生成低噪声训练数据。这一方法解决了传统机器学习粗粒化（MLP-CG）模型长期存在的根本性难题，为CG建模带来了一场范式转变。其重要性体现在以下几个方面：

  • 克服瞬时力噪声：传统的力匹配方法直接拟合原子轨迹中的瞬时力，而瞬时力中包含大量的热力学噪声。这种噪声使得MLP难以学习底层的、真实的平均力势（PMF），导致训练误差居高不下，模型精度受限。本文通过对CG珠子施加约束并进行长时间（1-10纳秒）的时间平均，使得力信号从瞬时噪声中解耦出来，收敛到PMF所定义的系综平均力。这些平均力更加平滑、真实，极大提升了训练数据的质量。
  • 实现高精度：低噪声训练数据使得NEP模型能够更准确地学习势能面，达到与AA模型（基于DFT或高精度QM数据训练）相当的精度。这在液态水案例中得到了充分验证，力与应力RMSE显著降低，并且模型能够准确重现水在宽泛压力范围内的结构和状态方程。
  • 简化超参数选择：在噪声数据中，训练误差的波动会掩盖模型收敛的真实情况，使超参数（如神经网络层数、神经元数量、截止半径、正则化强度）的优化变得复杂。低噪声数据使得训练曲线更加清晰，收敛性评估更直观，从而更容易选择最佳超参数组合。
  • 提高数据效率：由于数据质量高，模型所需的训练结构数量显著减少。液态水模型仅用5个结构（2560个水珠）就能实现出色的性能，而传统方法可能需要数千个快照。这种数据效率的提升对于昂贵的AA计算来说至关重要。
  • 增强模型迁移性：通过在不同压力下采集的低噪声数据进行训练，模型能够学习到在各种热力学条件下的正确响应。液态水模型在0.5 GPa训练限制外，能成功外推至1 GPa，这证明了模型在跨压强条件下的良好迁移性。

总而言之，低噪声训练数据是本文方法成功的基石，它不仅提高了MLP-CG模型的精度和效率，更使得粗粒化模型从拟合噪声走向学习真实物理，推动了整个领域的发展。

2. 维里校正的深层考量

维里校正（Virial Correction）是确保粗粒化模型准确重现宏观热力学性质（如压力和密度）的关键步骤。其背后涉及粗粒化过程中自由能和动能的精妙处理。

  • 粗粒化与自由度损失：在从全原子系统到粗粒化系统的映射过程中，每个CG珠子代表一组AA原子。这意味着CG模型中包含了更少的有效自由度（例如，水分子的三个原子映射为一个珠子，丢失了两个旋转和三个振动自由度）。根据能量均分定理，这些被积分掉的原子自由度的动能贡献被忽略了。在全原子系统中，系统的维里（与压力直接相关）包含两部分：动能贡献（与温度和粒子数相关）和构型贡献（与粒子间相互作用力相关）。当粗粒化时，CG珠子的动能只代表CG珠子自身的运动，而其内部原子运动的动能（“理想气体”部分）在形式上被移除了。
  • 压力低估与密度高估：如果不对这种自由度损失进行补偿，粗粒化模型将倾向于低估系统的总动能贡献，从而导致计算出的压力低于真实值。为了在特定温度下达到目标压力，系统将通过压缩自身（即增加密度）来“补偿”这种压力不足，最终表现为密度过高（如图2b所示）。
  • 校正机制：本文提出的维里校正 $(N_{AA} – N_{CG}) k_B T I$ 正是用于弥补这部分损失的动能贡献。其中，$N_{AA} – N_{CG}$ 代表了粗粒化过程中被有效积分掉的自由度数量，$k_B T$ 是每个自由度对应的动能，$I$ 是单位张量。通过将这部分能量加回维里，CG模型能够更准确地反映原始AA系统的总动能贡献，从而修正压力计算，确保模型在宽泛的压力范围内准确重现密度-压力状态方程。

这个校正项的引入，不仅是一个技术细节，更是对粗粒化物理深刻理解的体现。它使得粗粒化模型不仅能在微观结构层面与AA模型一致，也能在宏观热力学性质层面提供准确预测，极大地增强了模型的鲁棒性和迁移性。

3. 各向异性处理的多层次策略

各向异性是许多复杂材料和生物分子体系的固有特性。本文在粗粒化模型中处理各向异性采取了多层次的策略，并指出了未来发展的方向：

  • 简单各向同性珠子（初始尝试）：在C60单分子层案例中，最初尝试使用单一各向同性珠子来代表所有C60分子。结果表明，这种简单方法无法捕获系统固有的各向异性键合（例如不同的键合类型沿不同方向），导致应力预测出现系统性误差和精度显著下降。这证实了在处理各向异性系统时，简单各向同性模型存在局限性。
  • 通过区分珠子类型编码各向异性（有效策略）：针对C60单分子层的具体情况，研究团队引入了两种不同类型的各向同性CG珠子，分别对应于晶胞中晶体学上不等价的C60分子。这种策略成功地将应力误差降低了一个数量级，并准确捕获了各向异性热导率。这表明，通过巧妙地将晶体学信息编码到不同的珠子类型中，即使珠子本身是各向同性的，也能在一定程度上反映系统的各向异性。
  • 真正各向异性珠子（未来方向）：论文明确指出，当前模型仍使用各向同性珠子，对于需要直接捕获定向自由度（如聚合物、液晶）的系统，现有策略可能不足。未来的发展方向是扩展NEP框架，使其能够容纳具有旋转自由度和各向异性形状的珠子表示。这意味着需要开发新的描述符来捕获珠子的取向，并在势能函数中考虑这些额外的自由度，从而实现更普适的各向异性CG建模。

这种从尝试到改进再到展望的策略，展示了在粗粒化建模中处理复杂物理现象的渐进过程，以及MLP框架在适应这些需求方面的潜力。

4. 多尺度模拟的革命性潜力：以金纳米线断裂为例

NEP-AACG模型通过无缝整合全原子（AA）和粗粒化（CG）自由度，为多尺度模拟提供了一个统一且高效的框架。金纳米线断裂的案例生动地展示了这一方法的革命性潜力：

  • 桥接不同尺度：材料变形和断裂通常是典型的多尺度现象。裂纹萌生和局部断裂发生在原子尺度，而应力波传播、能量耗散和宏观形变则发生在更大的尺度上。纯AA模拟由于计算成本过高，无法模拟足够大的系统和足够长的时间来观察宏观断裂过程；纯CG模拟则无法提供断裂核心区域的原子级细节。
  • 克服边界条件伪影：在纳米线断裂模拟中，传统上为了模拟无限大系统，通常需要复杂的周期性边界条件或外部势能。NEP-AACG通过在断裂核心区域保持AA分辨率，同时在远端使用CG区域作为“储层”，有效地将原子级区域与宏观环境连接起来，避免了传统AA模拟中由于有限尺寸效应引起的边界条件伪影。
  • 模拟实验相关过程：该模型使得金纳米线断裂模拟能够在实验相关的 $10^7 s^{-1}$ 应变速率下进行。这在纯AA模拟中几乎是无法企及的，因为达到如此高的应变速率需要极长的模拟时间，或在原子尺度上进行非常小的系统模拟，这会引入严重的尺寸效应。
  • 统一的自由能曲面：NEP-AACG的关键优势在于它学习了一个跨AA和CG分辨率的一致自由能曲面。这意味着AA和CG区域之间的过渡是平滑的，没有界面伪影或能量不连续性，这对于准确模拟跨尺度耦合现象至关重要。

金纳米线断裂模拟是NEP-AACG框架的有力证明，它展示了多尺度方法如何使研究人员能够以前所未有的细节和规模探索复杂材料行为，从而加速材料设计和优化。

5. 软件生态系统与前沿工具

本研究的成功离不开强大的软件生态系统支持，特别是**GPUMD（GPU-accelerated Molecular Dynamics）**软件包。GPUMD作为一个高性能MD模拟平台，为NEP框架的开发和部署提供了坚实的基础。

  • GPU加速：GPUMD充分利用GPU的并行计算能力，极大地加速了MD模拟，这对于训练MLP和执行大规模AA或CG模拟至关重要。研究中提到的数百到数千纳秒/天的模拟速度，正是GPU加速的直接体现。
  • NEP框架集成：GPUMD集成了NEP框架，使得研究人员可以方便地定义、训练和使用NEP势能。NEP的描述符构建、神经网络评估和参数优化（通过SNES）都高度优化，确保了效率和精度。
  • 多功能性：GPUMD不仅支持NEP，还支持多种MD系综和分析工具，使其成为一个多功能的计算平台，适用于广泛的材料模拟问题。
  • 开源与社区：GPUMD作为开源软件（GitHub链接已提供），促进了科学的可重复性和协作。开放的源代码使得其他研究人员可以审查、修改和扩展其功能，共同推动领域发展。

NEP-CG和NEP-AACG方法与GPUMD的结合，代表了高性能计算与机器学习在分子模拟领域融合的典范。这种工具的进步不仅加速了当前的研究，也为未来更复杂、更精细的分子模拟提供了强大的平台。

6. 展望与更广泛的意义

本文提出的方法为分子模拟领域开辟了广阔的前景，超越了当前的工作局限，有望在多个科学领域产生深远影响。

  • 生物分子模拟：蛋白质折叠、药物-靶点相互作用、膜蛋白动力学等生物分子过程涉及广泛的时间尺度。NEP-AACG可以提供一个统一的框架，在结合区域保持原子级精度，而在溶剂或远离相互作用位点的地方使用CG模型，从而实现对复杂生物系统的更长时间、更大规模的模拟。
  • 聚合物和软物质：聚合物的自组装、液晶的相变、胶体颗粒的行为等都强烈依赖于其柔性和各向异性。虽然目前模型使用各向同性珠子，但未来各向异性珠子表示的开发将极大地扩展NEP在这些系统中的应用，帮助理解和设计新型软物质材料。
  • 材料科学与工程：从金属的断裂韧性到电池材料的离子输运，再到催化剂的表面反应，许多材料现象都是多尺度的。NEP-AACG能够为材料设计提供精确的预测工具，加速新材料的发现和优化。
  • 相变和临界现象：通过开发温度可迁移的CG模型，NEP可以用于模拟涉及温度变化的相变过程，如液体的冻结、熔化或液-液相变，这对于理解材料的热力学性质至关重要。
  • 与量子化学的桥梁：作为一名量子化学研究人员，我深知从高精度量子力学计算（如DFT、CCSD(T)）到大规模经典MD模拟的鸿沟。NEP-AA模型能够以QM精度复现原子相互作用，而NEP-CG/AACG则在保持这种物理保真度的前提下，极大地扩展了可模拟的时空尺度。这提供了一条从第一性原理数据出发，构建高度准确且高效的多尺度模型，真正实现了“从头算到宏观”的计算材料科学愿景。

综上所述，NEP-CG和NEP-AACG代表了机器学习势能在分子模拟领域的一个重要里程碑。通过解决核心的噪声问题，并提供一个灵活的多尺度框架，该工作不仅提升了现有模拟的精度和效率，更为未来探索更复杂、更具挑战性的科学问题奠定了坚实基础。