DeePMD-kit v2：深度势能模型的软件包

翻译自 J. Chem. Phys., 2023, 159, 54801.

**作者：**Jinzhe Zeng, Duo Zhang, Denghui Lu, Pinghui Mo, Zeyu Li, Yixiao Chen, Marián Rynik, Li'ang Huang, Ziyao Li, Shaochen Shi, Yingze Wang, Haotian Ye, Ping Tuo, Jiabin Yang, Ye Ding, Yifan Li, Davide Tisi, Qiyu Zeng, Han Bao, Yu Xia, Jiameng Huang, Koki Muraoka, Yibo Wang, Junhan Chang, Fengbo Yuan, Sigbjørn Løland Bore, Chun Cai, Yinnian Lin, Bo Wang, Jiayan Xu, Jia-Xin Zhu, Chenxing Luo, Yuzhi Zhang, Rhys E. A. Goodall, Wenshuo Liang, Anurag Kumar Singh, Sikai Yao, Jingchao Zhang, Renata Wentzcovitch, Jiequn Han, Jie Liu, Weile Jia, Darrin M. York, Weinan E, Roberto Car, Linfeng Zhang, and Han Wang

Jinzhe Zeng 译

© 2023 Jinzhe Zeng. 译文翻译自原文（原文作者版权所有，以CC BY 4.0协议授权），并以CC BY-SA 4.0协议授权。

DeePMD-kit是一个强大的开源软件包，它使用被称为深度势能模型（DeepPotential，简称DP）的机器学习势能（MLP）来促进分子动力学模拟。该软件包于2017年发布，并广泛应用于物理学、化学、生物学和材料科学领域，用于研究原子尺度体系。DeePMD-kit的当前版本提供了许多先进功能，如DeepPot-SE模型、基于注意力机制的描述符、混合描述符、拟合张量性质、类型嵌入、模型偏差、深度势能·范围校正模型（DPRc）、深度势能·长程模型（DPLR）、支持自定义操作符的GPU加速、模型压缩、非冯·诺依曼分子动力学（NVNMD）和改进的易用性，包括文档、编译的二进制软件包、图形用户界面（GUI）和应用程序编程接口（API）。本文概述了DeePMD-kit软件包的当前主要版本，突出其特点和技术细节。此外，本文还介绍了作为代表应用的分子动力学的通用流程，评估了不同模型的准确性和效率，并讨论了正在进行的开发。

1.引言

近年来，机器学习势能（MLP）模型的日益流行，彻底改变了分子动力学（MD）模拟在多个领域的应用，包括神经网络势能（NNP）^2-19，消息传递模型^7,21-24，以及其它机器学习模型^26-28。已开发出许多软件包来支持MLP的使用。^13,30-40MLP被广泛采用的主要原因之一是其出色的速度和准确性，优于传统的分子力学（MM）和从头算量子力学（QM）方法。^41,42因此，基于MLP的分子动力学（MD）模拟在该领域已经无处不在，并越来越被认为是研究原子体系的有价值工具。^44-49

DeePMD-kit是一个开源软件包，用于利用神经网络势能进行分子动力学（MD）模拟。该软件包于2017年首次发布²⁹，自那以后经历了快速发展，并得到了许多开发者的贡献。DeePMD-kit实现了一系列被称为深度势能（DP）模型的MLP模型^9,10,51-54，已被广泛应用于物理学、化学、生物学和材料科学等领域，用于研究各种原子体系。这些体系包括金属材料⁵⁵、非金属无机材料^57-60、水^62-71、有机体系^10,72、溶液^52,74-76、气相体系^78-80、大分子体系^81,82和界面^84-87等。此外，DeePMD-kit能够模拟包含元素周期表中几乎所有元素的体系⁵¹，在广泛的温度和压力范围内运行⁸⁸，并能够处理类似药物的分子^72,89、离子^73,76、过渡态^75,77和激发态⁹⁰。因此，DeePMD-kit是一个强大而多功能的工具，可用于模拟各种原子体系。在这里，我们展示三个示例，以突显其多样化的应用。

理论研究水相图是一个巨大的挑战，因为需要高精度的水分子间相互作用模型^91,92。因此，它对模型的准确性提出了极高的要求，并提供了验证分子动力学模拟中所需的软件实现用于相图计算的手段⁹²。Zhang等人⁸⁸使用DeePMD-kit构建了一个水体系的深度势能模型，涵盖了从0到2400K和0到50GPa的一系列热力学状态。该模型采用了SCAN近似交换相关泛函的密度泛函理论数据进行训练，在相关热力学范围内具有一致的准确性（均方根误差小于1meV/H${}_2$O）。此外，它准确预测了范围内的液态、分子态和离子态以及除III相和XV相以外所有稳定的冰多晶型物。该研究广泛调查了从冰VII向VII"的两个一级相变以及从VII"向离子态流体的原子机制，利用了该模型在预测水分子电离方面的能力和高准确性。

另一个具有挑战性的领域是凝聚相分子动力学模拟，长程相互作用对于模拟凝聚态异质体系至关重要。电静相互作用是最长的，但也是最容易理解的，存在着用于点电荷⁹³、多极^94,95和量子力学^96,97级别高效计算的线性缩放方法。可以开发快速半经验量子力学方法^98,99，可以准确高效地模拟长程电荷密度和多体效应，但在中程范围（通常小于8Å）可能缺乏定量准确性。这限制了凝聚相模拟方法的预测能力。Zeng等人⁵²提出了一种名为深度势能·范围校正（DPRc）的新$\Delta$-MLP方法，用于与量子力学/分子力学（QM/MM）的组合势能模型相结合，该方法校正了由快速线性缩放的低级半经验QM/MM理论得到的势能，使其成为高级从头算QM/MM理论。与许多正在发展中的$\Delta$-MLP模型不同，该DPRc模型校正了QM-QM和QM-MM相互作用，以一种在MM原子进入（或离开）QM区域接近时保持能量守恒的方式。这使得该模型能够作为分子模拟软件中的中程势能校正轻松集成，该软件使用了非键列表，即对于每个原子，有固定截断距离（通常为8-12 Å）内的其他原子的列表。训练过的具有6 Å范围校正的DPRc模型被应用于长时间尺度下模拟溶液中的RNA2${}^{\prime }$-O转磷酸酯化反应⁷⁵，并通过广义加权热力学扰动（gwTP）方法的推广获得了更好的自由能估计¹⁰⁰。最近，Zeng等人⁷²使用一种名为量子深度势能相互作用（QD$\pi$）的$\Delta$-MLP校正模型训练了包括互变型和质子化状态在内的类药物分子，发现其优于其他半经验方法和纯MLP模型⁸⁹。

第三个重要的应用是在纳秒时间尺度上进行大规模反应分子动力学模拟，并由此建立复杂反应体系的交织反应网络¹⁰¹，而不是仅关注研究单个反应。这些模拟需要精确且计算效率高的势能模型，覆盖可能发生的化学反应的化学空间。Zeng等人⁷⁷开发了一个深度势能模型用于模拟1纳秒的甲烷燃烧反应，并使用了ReacNetGenerator软件包¹⁰²在空间和时间上识别出了798种不同的化学反应。在复杂反应过程中，采用了同步学习流程¹⁰³，并证明其在探索已知和未知的化学空间中的作用至关重要。随后，进一步的研究将这些模拟扩展到了更复杂的反应体系，包括线性烷烃热解⁷⁸、爆炸物分解^79,104和多环芳烃生长⁸⁰。

与最初发布的版本²⁹相比，DeepMD-kit发展迅速，当前版本（v2.2.1）提供了广泛的功能。这些功能包括：DeepPot-SE、基于注意力机制的描述符、混合描述符^10,50,51,53，拟合张量性质^105,106，类型嵌入，模型偏差^103,107，深度势能·范围校正（DPRc）模型^52,75，深度势能·长程（DPLR）模型⁵³，针对自定义操作符的图形处理单元（GPU）支持¹⁰⁸，模型压缩¹⁰⁹，非冯·诺伊曼分子动力学（NVNMD）¹¹⁰，以及各种可用性的改进，如文档、编译的二进制软件包、图形用户界面（GUI）和应用程序编程接口（API）。本文概述当前主要版本的DeePMD-kit的特点和技术细节，重点介绍其功能，并讨论正在进行的开发。

2.功能

在这一部分，我们从组件的角度介绍DeePMD-kit中的功能（如图 [1]所示）。一个组件表示计算单元，并组织为DeePMD-kit软件包中的一个Python类，相应的TensorFlow静态图将在运行时创建。

图1 DeePMD-kit软件包的各组件。箭头的方向表示组件之间的依赖关系。蓝色方框表示一个可选组件。

2.1 模型

深度势能（DeepPotential，简称DP）模型，记作$\mathcal{M}$，可以一般地表示为

$${\mathbf{y}}_i=\mathcal{M}({\mathbf{x}}_i,\{{\mathbf{x}}_j{\}}_{j\in n(i)};\mathbf{\theta })=\mathcal{F}(\mathcal{D}({\mathbf{x}}_i,\{{\mathbf{x}}_j{\}}_{j\in n(i)};{\mathbf{\theta }}_d);{\mathbf{\theta }}_f),\text{\hspace{1em}(1)}$$ 其中${\mathbf{y}}_i$为拟合属性，$\mathcal{F}$为拟合网络（在第2.1.3节介绍），$\mathcal{D}$为描述符（在第2.1.2节介绍）。$\mathbf{x}=({\mathbf{r}}_i,{\alpha }_i)$，其中${\mathbf{r}}_i$为笛卡尔坐标，${\alpha }_i$为化学种类，表示原子$i$的自由度。原子$i$的邻近原子的索引（即在一定截断半径内的原子）由符号$n(i)$给出。注意，笛卡尔坐标可以是周期边界条件（PBC）下的，也可以是真空中的（开放边界条件下）。网络参数由$\mathbf{\theta }=\{{\mathbf{\theta }}_d,{\mathbf{\theta }}_f\}$表示，其中${\mathbf{\theta }}_d$和${\mathbf{\theta }}_f$分别给出了描述符的网络参数（如果有）和拟合网络的网络参数。根据公式(1)，可以通过如下方式计算体系的全局性质：

$$\mathbf{y}=\sum _{i=1}^N{\mathbf{y}}_i,\text{\hspace{1em}(2)}$$ 其中$N$为一帧中原子的数量。例如，如果$y_i$表示原子$i$的势能贡献，则$y$给出了该帧的总势能。在以下文本中，$N_c$为预期的最大邻近原子数，对于所有帧中的所有原子来说，该常数相同。如果邻近原子的数量小于$N_c$，则会将矩阵的维度填充为$N_c$。

2.1.1 神经网络

神经网络（NN）函数$\mathcal{N}$是多个层${\mathcal{L}}^{(i)}$的组合：

$$\mathcal{N}={\mathcal{L}}^{(n)}\circ {\mathcal{L}}^{(n-1)}\circ \cdots \circ {\mathcal{L}}^{(1)}.\text{\hspace{1em}(3)}$$ 在DeePMD-kit软件包中，层$\mathcal{L}$可以是以下形式之一，具体取决于是否使用ResNet¹¹¹以及节点的数量：

$$\mathbf{y}=\mathcal{L}(\mathbf{x};\mathbf{w},\mathbf{b})=\{\begin{array}{ll}\hat{\mathbf{w}}\odot \mathbf{\varphi }({\mathbf{x}}^T\mathbf{w}+\mathbf{b})+\mathbf{x},& \text{ResNet and }{N}_2=N_1,\\ \hat{\mathbf{w}}\odot \mathbf{\varphi }({\mathbf{x}}^T\mathbf{w}+\mathbf{b})+\{\mathbf{x},\mathbf{x}\},& \text{ResNet and }{N}_2=2N_1,\\ \hat{\mathbf{w}}\odot \mathbf{\varphi }({\mathbf{x}}^T\mathbf{w}+\mathbf{b}),& \text{otherwise},\end{array}\text{\hspace{1em}(4)}$$ 其中$\mathbf{x}\in {\mathbb{R}}^{N_1}$是输入向量，$\mathbf{y}\in {\mathbb{R}}^{N_2}$是输出向量。$\mathbf{w}\in {\mathbb{R}}^{N_1\times N_2}$和$\mathbf{b}\in {\mathbb{R}}^{N_2}$分别表示权重和偏置，两者都是可训练的。$\hat{\mathbf{w}}\in {\mathbb{R}}^{N_2}$可以是可训练向量，表示跳跃连接中的"时间步"；或者是由所有元素为1的向量$1=\{1,1,\dots ,1\}$，表示禁用时间步。$\mathbf{\varphi }$是激活函数。理论上，激活函数可以是任何形式，DeePMD-kit软件包中提供了以下函数：双曲正切（tanh）、修正线性单元（ReLU）¹¹²、ReLU6、softplus¹¹³、sigmoid、高斯误差线性单元（GELU）¹¹⁴和恒等函数。在这些激活函数中，ReLU和ReLU6的一阶导数不连续，其它函数的一阶导数和二阶导数均连续。

2.1.2 描述符

DeePMD-kit支持多种原子描述符，包括局部框架描述符，二体和三体嵌入DeepPot-SE描述符，基于注意力机制的描述符和混合描述符，后者定义为多个描述符的组合。在以下文本中，我们使用${\mathcal{D}}^i=\mathcal{D}({\mathbf{x}}_i,\{{\mathbf{x}}_j{\}}_{j\in n(i)};{\mathbf{\theta }}_d)$表示原子$i$的描述符。

2.1.2.1 局部框架。

局部框架描述符${\mathcal{D}}^i\in {\mathbb{R}}^{N_c\times \{1,4\}}$（有时简称为DPMD模型），是DP描述符的第一个版本⁹，通过使用完整信息或仅径向信息构建：

$$({\mathcal{D}}^i{)}_j=\{\begin{array}{ll}\{\begin{array}{cccc}\frac{1}{r_{ij}}& \frac{x_{ij}}{r_{ij}}& \frac{y_{ij}}{r_{ij}}& \frac{z_{ij}}{r_{ij}}\end{array}\},& \text{完整},\\ \{\begin{array}{c}\frac{1}{r_{ij}}\end{array}\},& \text{仅径向},\end{array}\text{\hspace{1em}(5)}$$ 其中$(x_{ij},y_{ij},z_{ij})$是原子$i$和原子$j$之间相对位置的三个笛卡尔坐标，即局部框架中${\mathbf{r}}_{ij}={\mathbf{r}}_i-{\mathbf{r}}_j=(x_{ij},y_{ij},z_{ij})$。$r_{ij}=|{\mathbf{r}}_{ij}|$是其范数。在公式(5)中，邻近原子$j$的顺序按照离原子$i$的距离升序排列。通过公式(6)，将${\mathbf{r}}_{ij}$从全局参考坐标${\mathbf{r}}_{ij}^0$转换为

$${\mathbf{r}}_{ij}={\mathbf{r}}_{ij}^0\cdot R_i,\text{\hspace{1em}(6)}$$其中

$$R_i=\{{\mathbf{e}}_{i1},{\mathbf{e}}_{i2},{\mathbf{e}}_{i3}\}\text{\hspace{1em}(7)}$$ 是由以下方式构造的旋转矩阵： $${\mathbf{e}}_{i1}=\mathbf{e}({\mathbf{r}}_{i,a(i)}),\text{\hspace{1em}(8)}$$ $${\mathbf{e}}_{i2}=\mathbf{e}({\mathbf{r}}_{i,b(i)}-({\mathbf{r}}_{i,b(i)}\cdot {\mathbf{e}}_{i1}){\mathbf{e}}_{i1}),\text{\hspace{1em}(9)}$$ $${\mathbf{e}}_{i3}={\mathbf{e}}_{i1}\times {\mathbf{e}}_{i2},\text{\hspace{1em}(10)}$$ 其中$\mathbf{e}({\mathbf{r}}_{ij})={\mathbf{r}}_{ij}/r_{ij}$表示对向量进行归一化的操作。$a(i)\in n(i)$和$b(i)\in n(i)$是用于定义原子$i$的局部参考系轴的两个轴原子，一般来说，它们是两个最近的原子，不考虑它们的种类，以及中心原子$i$。

局部框架描述符的局限性在于在截断半径处及最近邻之间的顺序交换（即$a(i)$和$b(i)$的交换）不平滑，因此其使用受到限制。我们注意到，在所有DP描述符中，局部框架描述符是唯一一个非平滑的描述符，因此我们建议在通常的体系中使用其它描述符。

2.1.2.2 两体嵌入DeepPot-SE。

DP描述符的两体嵌入平滑版本${\mathcal{D}}^i\in {\mathbb{R}}^{M\times M_{<}}$通常被称为DeepPot-SE描述符¹⁰。注意，描述符是原子$i$的局部环境的多体表示；我们称之为"两体嵌入"，是因为嵌入网络仅考虑原子$i$和$j$之间的距离（见下文），但并不意味着描述符只包含$i$和其邻居之间的一对信息。描述符可以使用完整信息或仅径向信息，表示为不同的形式：

$${\mathcal{D}}^i=\{\begin{array}{ll}\frac{1}{N_c^2}({\mathcal{G}}^i{)}^T{\mathcal{R}}^i({\mathcal{R}}^i{)}^T{\mathcal{G}}_{<}^i,& \text{完整},\\ \frac{1}{N_c}{\sum }_j({\mathcal{G}}^i{)}_{jk},& \text{仅径向},\end{array}\text{\hspace{1em}(11)}$$ 其中${\mathcal{R}}^i\in {\mathbb{R}}^{N_c\times \{1,4\}}$是坐标矩阵，${\mathcal{R}}^i$的每一行可以构造为：

$$({\mathcal{R}}^i{)}_j=\{\begin{array}{ll}\{\begin{array}{cccc}s(r_{ij})& \frac{s(r_{ij})x_{ij}}{r_{ij}}& \frac{s(r_{ij})y_{ij}}{r_{ij}}& \frac{s(r_{ij})z_{ij}}{r_{ij}}\end{array}\},& \text{完整},\\ \{\begin{array}{c}s(r_{ij})\end{array}\},& \text{仅径向},\end{array}\text{\hspace{1em}(12)}$$ 其中${\mathbf{r}}_{ij}={\mathbf{r}}_j-{\mathbf{r}}_i=(x_{ij},y_{ij},z_{ij})$为相对坐标，$r_{ij}=\Vert {\mathbf{r}}_{ij}\Vert$为其范数。切换函数$s(r)$定义为

$$s(r)=\{\begin{array}{ll}\frac{1}{r},& r<r_s,\\ \frac{1}{r}[x^3(-6x^2+15x-10)+1],& r_s\le r<r_c,\\ 0,& r\ge r_c,\end{array}\text{\hspace{1em}(13)}$$其中$x=\frac{r-r_s}{r_c-r_s}$在$r_s$处从0切换到$r_c$处的1，$[x^3(-6x^2+15x-10)+1]$在$r_s$处从1切换到$r_c$处的0。切换函数$s(r)$在二阶导数连续的意义下是光滑的。关于五次多项式$[x^3(-6x^2+15x-10)+1]$的推导过程可以在附录1中找到。

嵌入矩阵${\mathcal{G}}^i\in {\mathbb{R}}^{N_c\times M}$的每一行由$s(r_{ij})$的NN函数${\mathcal{N}}_{e,2}$的输出层的$M$个节点组成：

$$({\mathcal{G}}^i{)}_j={\mathcal{N}}_{e,2}(s(r_{ij})),\text{\hspace{1em}(14)}$$ 其中NN函数$\mathcal{N}$在式(4)中给出，下标"$e,2$"用于区分DP模型中其他NN函数。在式(14)中，网络参数未显式写出。${\mathcal{G}}_{<}^i\in {\mathbb{R}}^{N_c\times M_{<}}$仅取${\mathcal{G}}^i$的前$M_{<}$列以减小${\mathcal{D}}^i$的大小。$r_s$、$r_c$、$M$和$M_{<}$是用户提供的超参数。与局部框架描述符相比，DeepPot-SE在其定义域上在一阶和二阶导数连续。

2.1.2.3 三体嵌入DeepPot-SE。

三体嵌入DeepPot-SE描述符结合了键角信息，使模型更准确⁵⁰。描述符${\mathcal{D}}^i$可以表示为

$${\mathcal{D}}^i=\frac{1}{N_c^2}({\mathcal{R}}^i({\mathcal{R}}^i{)}^T):{\mathcal{G}}^i,\text{\hspace{1em}(15)}$$ 其中${\mathcal{R}}^i$由式(12)定义。目前，三体嵌入只支持${\mathcal{R}}^i$的全信息情况。类似于公式(14)，${\mathcal{G}}^i\in {\mathbb{R}}^{N_c\times N_c\times M}$中的每个元素都来自于NN${\mathcal{N}}_{e,3}$函数输出层的$M$个节点： $$({\mathcal{G}}^i{)}_{jk}={\mathcal{N}}_{e,3}(({\theta }_i{)}_{jk}),\text{\hspace{1em}(16)}$$ 其中$({\theta }_i{)}_{jk}=({\mathcal{R}}^i{)}_j\cdot ({\mathcal{R}}^i{)}_k$考虑到了两个邻居（$j$和$k$）之间的角度形式。公式(15)中的符号":"表示矩阵${\mathcal{R}}^i({\mathcal{R}}^i{)}^T$与张量${\mathcal{G}}^i$的前两个维度之间的收缩。网络参数在公式(16)中也没有明确写出。

2.1.2.4 处理由多个化学元素组成的体系。

对于一个包含多个化学元素（$|\{{\alpha }_i\}|>1$）的体系，在公式(14)和(16)中，嵌入网络${\mathcal{N}}_{e,\{2,3\}}$的参数如下：

$$({\mathcal{G}}^i{)}_j={\mathcal{N}}_{e,2}^{{\alpha }_i,{\alpha }_j}(s(r_{ij}))\mathrm{or}({\mathcal{G}}^i{)}_j={\mathcal{N}}_{e,2}^{{\alpha }_j}(s(r_{ij})),\text{\hspace{1em}(17)}$$ $$({\mathcal{G}}^i{)}_{jk}={\mathcal{N}}_{e,3}^{{\alpha }_j,{\alpha }_k}(({\theta }_i{)}_{jk}).\text{\hspace{1em}(18)}$$ 因此，将会有$N_t^2$个或者$N_t$个嵌入网络，其中$N_t$是化学元素的数量。为了提高矩阵运算的性能，$n(i)$被划分为具有不同化学元素的块。每个维度为$N_c$的矩阵被划分为相应的块，并且每个块都单独填充到$N_c^{{\alpha }_j}$。这种方法的限制是当有大量化学物种时，例如OC2M数据集中的57种元素^115,116，嵌入网络的数量将会变成3249或57，需要大内存，降低了计算效率。

2.1.2.5 类型嵌入。

为了在化学元素数量庞大时减少NN参数的数量并提高计算效率，引入了类型嵌入$\mathcal{A}$，表示为原子类型$\alpha$的NN函数${\mathcal{N}}_t$：

$${\mathcal{A}}^i={\mathcal{N}}_t(\text{one hot}({\alpha }_i)),\text{\hspace{1em}(19)}$$ 其中${\alpha }_i$被转换为表示化学物种的独热向量，然后作为输入提供给NN（已在公式 (4)介绍）。基于公式(14)和(16)，中心原子和邻近原子的类型嵌入${\mathcal{A}}^i$和${\mathcal{A}}^j$被添加为嵌入网络${\mathcal{N}}_{e,\{2,3\}}$的额外输入：

$$({\mathcal{G}}^i{)}_j={\mathcal{N}}_{e,2}(\{s(r_{ij}),{\mathcal{A}}^i,{\mathcal{A}}^j\})\mathrm{or}({\mathcal{G}}^i{)}_j={\mathcal{N}}_{e,2}(\{s(r_{ij}),{\mathcal{A}}^j\}),\text{\hspace{1em}(20)}$$ $$({\mathcal{G}}^i{)}_{jk}={\mathcal{N}}_{e,3}(\{({\theta }_i{)}_{jk},{\mathcal{A}}^j,{\mathcal{A}}^k\}).\text{\hspace{1em}(21)}$$ 通过类型嵌入方法，所有化学元素共享相同的网络参数。

2.1.2.6 基于注意力的描述符。

基于注意力的描述符${\mathcal{D}}^i\in {\mathbb{R}}^{M\times M_{<}}$在可预训练的 DPA-1⁵¹ 模型中被提出，其定义如下：

$${\mathcal{D}}^i=\frac{1}{N_c^2}({\hat{\mathcal{G}}}^i{)}^T{\mathcal{R}}^i({\mathcal{R}}^i{)}^T{\hat{\mathcal{G}}}_{<}^i,\text{\hspace{1em}(22)}$$ 其中 ${\hat{\mathcal{G}}}^i$是经过额外的自注意力机制¹¹⁷处理后的嵌入矩阵${\mathcal{G}}^i$，${\mathcal{R}}^i$在公式 (12) 中定义。需要注意的是，在此描述符中，默认使用类型嵌入方法，从公式 (20) 中获得 ${\mathcal{G}}^i$。

为了执行自注意力机制，首先得到查询${\mathcal{Q}}^{i,l}\in {\mathbb{R}}^{N_c\times d_k}$，键${\mathcal{K}}^{i,l}\in {\mathbb{R}}^{N_c\times d_k}$，和值${\mathcal{V}}^{i,l}\in {\mathbb{R}}^{N_c\times d_v}$： $${\left({\mathcal{Q}}^{i,l}\right)}_j=Q_l\left({\left({\mathcal{G}}^{i,l-1}\right)}_j\right),\text{\hspace{1em}(23)}$$ $${\left({\mathcal{K}}^{i,l}\right)}_j=K_l\left({\left({\mathcal{G}}^{i,l-1}\right)}_j\right),\text{\hspace{1em}(24)}$$ $${\left({\mathcal{V}}^{i,l}\right)}_j=V_l\left({\left({\mathcal{G}}^{i,l-1}\right)}_j\right),\text{\hspace{1em}(25)}$$ 其中 $Q_l$，$K_l$，$V_l$表示三个可训练的线性变换，其输出的查询和维度为 $d_k$ 的键，以及维度为$d_v$ 的值，$l$ 为注意力层的索引。 注意力层的输入嵌入矩阵，记为${\mathcal{G}}^{i,0}$，选择为两体嵌入矩阵 (14)。

然后使用缩放点积注意力方法^117,118：

$$A({\mathcal{Q}}^{i,l},{\mathcal{K}}^{i,l},{\mathcal{V}}^{i,l},{\mathcal{R}}^{i,l})=\phi \left({\mathcal{Q}}^{i,l},{\mathcal{K}}^{i,l},{\mathcal{R}}^{i,l}\right){\mathcal{V}}^{i,l},\text{\hspace{1em}(26)}$$ 其中$\phi \left({\mathcal{Q}}^{i,l},{\mathcal{K}}^{i,l},{\mathcal{R}}^{i,l}\right)\in {\mathbb{R}}^{N_c\times N_c}$是注意力权重。在原始的注意力方法中，通常有$\phi \left({\mathcal{Q}}^{i,l},{\mathcal{K}}^{i,l}\right)=\mathrm{softmax}\left(\frac{{\mathcal{Q}}^{i,l}({\mathcal{K}}^{i,l}{)}^T}{\sqrt{d_k}}\right)$，其中$\sqrt{d_k}$是归一化温度。为了加入角度信息，稍作修改：

$$\phi \left({\mathcal{Q}}^{i,l},{\mathcal{K}}^{i,l},{\mathcal{R}}^{i,l}\right)=\mathrm{softmax}\left(\frac{{\mathcal{Q}}^{i,l}({\mathcal{K}}^{i,l}{)}^T}{\sqrt{d_k}}\right)\odot {\hat{\mathcal{R}}}^i({\hat{\mathcal{R}}}^i{)}^T,\text{\hspace{1em}(27)}$$ 其中${\hat{\mathcal{R}}}^i\in {\mathbb{R}}^{N_c\times 3}$表示归一化的相对坐标，${\hat{\mathcal{R}}}_j^i=\frac{{\mathbf{r}}_{ij}}{\Vert {\mathbf{r}}_{ij}\Vert }$，$\odot$表示逐元素相乘。

然后通过层归一化以残差的方式添加，最终得到自注意力局部嵌入矩阵${\hat{\mathcal{G}}}^i={\mathcal{G}}^{i,L_a}$，其中共有$L_a$个注意力层：

$${\mathcal{G}}^{i,l}={\mathcal{G}}^{i,l-1}+\mathrm{LayerNorm}(A({\mathcal{Q}}^{i,l},{\mathcal{K}}^{i,l},{\mathcal{V}}^{i,l},{\mathcal{R}}^{i,l})).\text{\hspace{1em}(28)}$$

2.1.2.7 混合描述符。

混合描述符${\mathcal{D}}_{\text{hyb}}^i$将多种类型的描述符连接成一个描述符：⁵³

$${\mathcal{D}}_{\text{hyb}}^i=\{\begin{array}{cccc}{\mathcal{D}}_1^i& {\mathcal{D}}_2^i& \cdots & {\mathcal{D}}_n^i\end{array}\}.\text{\hspace{1em}(29)}$$ 描述符的列表可以是不同类型的，也可以是具有不同参数的相同描述符。通过这种方式，可以为不同的描述符设置不同的截断半径。

2.1.2.8 压缩。

DP模型的压缩使用三种技术进行，即制表推理、运算符合并和精确的邻居索引，以提高在模型参数经过适当训练时的模型训练和推断性能 ¹⁰⁹。

为了获得更好的性能，对于输入维度为一的情况，可以使用制表函数求值来替代神经网络推断。由(14)定义的嵌入网络${\mathcal{N}}_{e,2}$和由(16)定义的嵌入网络${\mathcal{N}}_{e,3}$属于此类型。其思想是通过分段多项式拟合来近似神经网络的输出。将输入域（在$\mathbb{R}$上的一个紧致域）均匀地分为$L_c$个间隔，在每个间隔中应用一个五阶多项式$g_m^l(x)$来逼近神经网络函数的第$m$个输出分量：

$$g_m^l(x)=a_m^l{x}^5+b_m^l{x}^4+c_m^l{x}^3+d_m^l{x}^2+e_m^{l}x+f_m^l,x\in [x_l,x_{l+1}),\text{\hspace{1em}(30)}$$ 其中$l=1,2,\dots ,L_c$是区间的索引，$x_1,\dots ,x_{L_c},x_{L_c+1}$是区间的端点，$a_m^l$, $b_m^l$, $c_m^l$, $d_m^l$, $e_m^l$, 和 $f_m^l$是拟合参数。

拟合参数可以通过以下方程计算： $$a_m^l=\frac{1}{2\Delta x_l^5}[12h_{m,l}-6(y_{m,l+1}^{\prime }+y_{m,l}^{\prime })\Delta x_l+(y_{m,l+1}^{\prime \prime }-y_{m,l}^{\prime \prime })\Delta x_l^2],\text{\hspace{1em}(31)}$$ $$b_m^l=\frac{1}{2\Delta x_l^4}[-30h_{m,l}+(14y_{m,l+1}^{\prime }+16y_{m,l}^{\prime })\Delta x_l+(-2y_{m,l+1}^{\prime \prime }+3y_{m,l}^{\prime \prime })\Delta x_l^2],\text{\hspace{1em}(32)}$$ $$c_m^l=\frac{1}{2\Delta x_l^3}[20h_{m,l}-(8y_{m,l+1}^{\prime }+12y_{m,l}^{\prime })\Delta x_l+(y_{m,l+1}^{\prime \prime }-3y_{m,l}^{\prime \prime })\Delta x_l^2],\text{\hspace{1em}(33)}$$ $$d_m^l=\frac{1}{2}{y}_{m,l}^{\prime \prime },\text{\hspace{1em}(34)}$$ $$e_m^l=y_{m,l}^{\prime },\text{\hspace{1em}(35)}$$ $$f_m^l=y_{m,l},\text{\hspace{1em}(36)}$$ 其中$\Delta x_l=x_{l+1}-x_l$表示区间的大小。$h_{m,l}=y_{m,l+1}-y_{m,l}$。$y_{m,l}=y_m(x_l)$,$y_{m,l}^{\prime }=y_m^{\prime }(x_l)$, 和$y_{m,l}^{\prime \prime }=y_m^{\prime \prime }(x_l)$分别表示目标NN函数的第$m$个分量在区间点$x_l$处的值、一阶导数和二阶导数。一阶和二阶导数可以通过NN函数的反向传播轻松计算。

在标准DP模型推理中，以双体嵌入描述符为例，矩阵乘积$({\mathcal{G}}^i{)}^T\mathcal{R}$需要在寄存器和主机/设备内存之间传输张量${\mathcal{G}}^i$，由于GPU的相对较小的内存带宽，这通常成为计算的瓶颈。压缩DP模型将矩阵乘法$({\mathcal{G}}^i{)}^T\mathcal{R}$与制表推理步骤合并。具体来说，一旦计算出$({\mathcal{G}}^i{)}^T$的一列，就立即将其与寄存器中的一个环境矩阵的一行进行乘法运算，并将外积存放到$({\mathcal{G}}^i{)}^T\mathcal{R}$的结果中。通过操作符合并技术，避免了${\mathcal{G}}^i$的分配和在寄存器和主机/设备内存之间的内存移动。三体嵌入的操作符合并可以类似地推导出来。

环境（${\mathcal{R}}^i$）和嵌入（${\mathcal{G}}^i$）矩阵的第一个维度，$N_c$，是预期的最大邻居数。如果某个原子的邻居数小于$N_c$，则矩阵的相应位置填充零。实际上，如果真实的邻居数明显小于$N_c$，则在填充零的乘法运算中会花费大量的操作。在压缩DP模型中，在表格推理阶段精确地索引邻居数，进一步节省了计算成本。

2.1.3 拟合网络

拟合网络可以拟合体系的势能、力和维里，以及张量性质如偶极矩和电极化率。

2.1.3.1 拟合势能。

在DP模型 (1)中，我们使用拟合网络${\mathcal{F}}_0$将描述子${\mathcal{D}}^i$映射到一个数量，其中下标"0"表示输出为零阶张量（即标量）。然后可以使用该模型通过以下方式预测体系的总势能：

$$E=\sum _i{E}_i=\sum _i{\mathcal{F}}_0({\mathcal{D}}^i),\text{\hspace{1em}(37)}$$ 其中拟合网络的输出被视为原子势能贡献，即$E_i$。输出的标量也可以被视为原子上定义的其他标量性质，例如原子$i$的局部电荷。

在某些情况下，可以将原子特定或帧特定的参数，例如电子温度¹¹⁹，作为拟合网络的额外输入。我们用${\mathbf{P}}^i\in {\mathbb{R}}^{N_p}$（其中$N_p$是维度）和$\mathbf{Q}\in {\mathbb{R}}^{N_q}$（其中$N_q$是维度）分别来表示原子和帧特定的参数。

$$E_i={\mathcal{F}}_0(\{{\mathcal{D}}^i,{\mathbf{P}}^i,\mathbf{Q}\}).\text{\hspace{1em}(38)}$$ 原子力${\mathbf{F}}_i$和维里张量$\mathbf{\Xi }=({\Xi }_{\alpha \beta })$（如果应用了周期性边界条件）可以由势能$E$导出：

$$F_{i,\alpha }=-\frac{\partial E}{\partial r_{i,\alpha }},\text{\hspace{1em}(39)}$$ $${\Xi }_{\alpha \beta }=-\sum _{\gamma }\frac{\partial E}{\partial h_{\gamma \alpha }}{h}_{\gamma \beta },\text{\hspace{1em}(40)}$$ 其中，$r_{i,\alpha }$和$F_{i,\alpha }$分别表示原子$i$坐标和力的第$\alpha$个分量。$h_{\alpha \beta }$是模拟区域的第$\alpha$个基向量的第$\beta$个分量。

2.1.3.2 拟合张量性质。

为了表示一阶张量性质（即矢量性质），我们让拟合网络${\mathcal{F}}_1$输出一个$M$维的矢量，然后表示如下:

$$(T_i^{(1)}{)}_{\alpha }=\frac{1}{N_c}\sum _{j=1}^{N_c}\sum _{m=1}^M({\mathcal{G}}^i{)}_{jm}({\mathcal{R}}^i{)}_{j,\alpha +1}({\mathcal{F}}_1({\mathcal{D}}^i){)}_m,\alpha =1,2,3.\text{\hspace{1em}(41)}$$ 我们让拟合网络${\mathcal{F}}_2$输出一个$M$维的矢量，然后二阶张量性质（矩阵性质）可以表达为:

$$(T_i^{(2)}{)}_{\alpha \beta }=\frac{1}{N_c^2}\sum _{j=1}^{N_c}\sum _{k=1}^{N_c}\sum _{m=1}^M({\mathcal{G}}^i{)}_{jm}({\mathcal{R}}^i{)}_{j,\alpha +1}({\mathcal{R}}^i{)}_{k,\beta +1}({\mathcal{G}}^i{)}_{km}({\mathcal{F}}_2({\mathcal{D}}^i){)}_m,\alpha ,\beta =1,2,3,\text{\hspace{1em}(42)}$$ 其中，${\mathcal{G}}^i$和${\mathcal{R}}^i$可以在方程 (14)和(12)（完整情况）中找到。因此，张量拟合网络要求描述符与DeepPot-SE描述符具有相同或相似的形式。NN函数${\mathcal{F}}_1$和${\mathcal{F}}_2$在方程 (4)中给出。总张量$\mathbf{T}$（总偶极子${\mathbf{T}}^{(1)}$或总极化率${\mathbf{T}}^{(2)}$）是原子张量的总和：

$$\mathbf{T}=\sum _i{\mathbf{T}}_i.\text{\hspace{1em}(43)}$$ 张量模型可用于计算红外光谱¹⁰⁵和拉曼光谱¹⁰⁶。

2.1.3.3 处理由多种元素成分组成的体系。

与嵌入网络类似，如果不使用类型嵌入方法，则拟合网络参数是按化学元素分别的，并且有$N_t$组拟合网络参数。为了性能，原子按其化学成分${\alpha }_i$提前排序。以一个例子来说明，原子能量$E_i$根据方程 (38)表示如下：

$$E_i={\mathcal{F}}_0^{{\alpha }_i}({\mathcal{D}}^i).\text{\hspace{1em}(44)}$$ 当使用类型嵌入时，所有化学元素共享相同的网络参数，并且类型嵌入被插入到方程 (38)中的拟合网络的输入中：

$$E_i={\mathcal{F}}_0(\{{\mathcal{D}}^i,{\mathcal{A}}^i\}).\text{\hspace{1em}(45)}$$

2.1.4 深度势能·范围校正（DPRc）

深度势能·范围校正（DPRc）^52,75最初设计用于将来自快速的、线性扩展的低级半经验量子力学/分子力学（QM/MM）理论的势能修正为高级从头算QM/MM理论，在范围修正方式下量化修正短程和中程的非键相互作用，利用分子力学力场，如AMBER常用的非键列表进行分子动力学模拟的定性校正¹²⁰。通过这种方式，长程静电相互作用可以高效地使用粒子网格Ewald方法¹²⁰及其在多极^94,95和QM/MM^96,97势能上的推广来建模。在DPRc模型中，方程 (13)中的切换函数被修改为禁用MM-MM相互作用：

$$s_{\text{DPRc}}(r_{ij})=\{\begin{array}{ll}0,& \text{如果i∈MM∧j∈MM},\\ s(r_{ij}),& \text{其它情况},\end{array}\text{\hspace{1em}(46)}$$ 其中$s_{\text{DPRc}}(r_{ij})$是新的开关函数，$s(r_{ij})$是公式 (13)中的旧开关函数。 这确保了MM原子之间的力为零，即

$${\mathbf{F}}_{ij}=-\frac{\partial E}{\partial {\mathbf{r}}_{ij}}=0,i\in \text{MM}\wedge j\in \text{MM}.\text{\hspace{1em}(47)}$$ 公式 (38)中的拟合网络被修正以消除MM原子的能量偏差：

$$E_i=\{\begin{array}{ll}{\mathcal{F}}_0({\mathcal{D}}^i),& \text{if i∈QM},\\ {\mathcal{F}}_0({\mathcal{D}}^i)-{\mathcal{F}}_0(0),& \text{if i∈MM},\end{array}\text{\hspace{1em}(48)}$$其中 $0$ 是一个零矩阵。值得一提的是，DPRc的使用不仅限于其最初的设计用于QM/MM校正，还可以扩展到任何类似的相互作用¹²¹。

2.1.5 深度势能·长程（DPLR）

深度势能·长程（DPLR）模型将静电能添加到总能量中⁵³：

$$E=E_{\text{DP}}+E_{\text{ele}},\text{\hspace{1em}(49)}$$ 其中$E_{\text{DP}}$是短程成分，构建为公式 (37)中对$(E^{\ast }-E_{\text{ele}})$进行拟合的标准能量模型。$E_{\text{ele}}$是由一组高斯分布引入的静电能，在傅里叶空间中计算如下：

$$E_{\text{ele}}=\frac{1}{2\pi V}\sum _{m\ne 0,\Vert m\Vert \le L}\frac{\exp (-{\pi }^2{m}^2/{\beta }^2)}{m^2}{S}^2(m),\text{\hspace{1em}(50)}$$ 其中$\beta$是一个可自由调节的参数，控制高斯分布的扩展程度。 $L$是Fourier空间中的截断参数，结构因子$S(m)$可表示为

$$S(m)=\sum _i{q}_i{e}^{-2\pi ım{\mathbf{r}}_i}+\sum _n{q}_n{e}^{-2\pi ım{\mathbf{W}}_n},\text{\hspace{1em}(51)}$$ 其中$ı=\sqrt{-1}$ 表示虚数单位，${\mathbf{r}}_i$ 表示离子坐标，$q_i$是离子$i$的电荷，$W_n$ 是第$n$个Wannier重心(WC)，可从方程 (42)中的分离二极子模型中得到。根据证明，高斯近似引入的静电能误差主要来自二极-四极相互作用的求和，其衰减形式为$r^{-4}$，其中$r$为二极与四极之间的距离 ⁵³。

2.1.6 使用成对势函数进行插值

在辐射损伤模拟等应用中，原子间的距离可能变得过近，使得DFT计算失败。在这种情况下，通常会用经验势函数来代替DP模型，如辐射损伤模拟中的Ziegler-Biersack-Littmark(ZBL) ¹²²屏蔽核斥力势能 ¹²³。DeePMD-kit软件包支持DP和经验成对势函数之间的插值关系，用下式表示：

$$E_i=(1-w_i)E_i^{\mathrm{DP}}+w_i{E}_i^{\mathrm{pair}},\text{\hspace{1em}(52)}$$ 其中$w_i$为插值权重，$E_i^{\mathrm{pair}}$为原子对势$u^{\mathrm{pair}}(r)$的贡献，即

$$E_i^{\mathrm{pair}}=\sum _{j\in n(i)}{u}^{\mathrm{pair}}(r_{ij}).\text{\hspace{1em}(53)}$$ 插值权重$w_i$的定义如下： $$w_i=\{\begin{array}{ll}1,& {\sigma }_i<r_a,\\ u_i^3(-6u_i^2+15u_i-10)+1,& r_a\le {\sigma }_i<r_b,\\ 0,& {\sigma }_i\ge r_b,\end{array}\text{\hspace{1em}(54)}$$ 其中，$u_i=({\sigma }_i-r_a)/(r_b-r_a)$。附录2中可以找到等式 (54)的推导过程。在范围$[r_a,r_b]$内，DP模型平滑关闭，并且成对势函数从$r_b$平滑过渡到$r_a$。${\sigma }_i$是原子$i$与其邻居之间距离的软最小值，

$${\sigma }_i=\frac{\sum _{j\in n(i)}{r}_{ij}{e}^{-r_{ij}/{\alpha }_s}}{\sum _{j\in n(i)}{e}^{-r_{ij}/{\alpha }_s}},\text{\hspace{1em}(55)}$$ 其中，${\alpha }_s$是可调整的原子间距离的尺度。成对势能$u^{\text{pair}}(r)$由用户定义的表格定义，该表格提供了在从0到截断距离的等间隔离散化网格上的$u^{\text{pair}}$的值。

2.2 训练器

基于等式 (1)中定义的DP模型$\mathcal{M}$，还需要定义一个训练器来训练模型中的参数，包括等式 (4)中的权重和偏置项。训练器应该给出学习率$\gamma$，损失函数$L$和训练过程。

2.2.1 学习率

学习率$\gamma$按指数衰减：

$$\gamma (\tau )={\gamma }^0{r}^{\lfloor \tau /s\rfloor },\text{\hspace{1em}(56)}$$ 其中，$\tau \in \mathbb{N}$是训练步骤的索引，${\gamma }^0\in \mathbb{R}$是第一步的学习率，衰减率$r$由以下公式给出：

$$r={\left(\frac{{\gamma }^{\text{stop}}}{{\gamma }^0}\right)}^{\frac{s}{{\tau }^{\text{stop}}}},\text{\hspace{1em}(57)}$$ 其中，${\tau }^{\text{stop}}\in \mathbb{N}$，${\gamma }^{\text{stop}}\in \mathbb{R}$和$s\in \mathbb{N}$分别是停止步骤、停止学习率和衰减步骤，所有这些都是事先提供的超参数。

2.2.2 损失函数

损失函数$L$是由不同拟合属性损失$L_p$的加权和给出：

$$L(\mathbf{x};\mathbf{\theta },\tau )=\frac{1}{\mathcal{B}}\sum _{k\in \mathcal{B}}\sum _{\eta }{p}_{\eta }(\tau )L_{\eta }({\mathbf{x}}^k;\mathbf{\theta }),\text{\hspace{1em}(58)}$$ 其中$\mathcal{B}$是数据的小批量。$\mathbf{x}=\{{\mathbf{x}}^k\}$是数据集。${\mathbf{x}}^k=({\mathbf{x}}_1^k,\dots ,{\mathbf{x}}_N^k)$是从该集合中抽取的单个数据帧，由原子的所有自由度组成。$\eta$表示要拟合的属性。 对于每个属性，$p_{\eta }$是一个前置因子，给定为

$$p_{\eta }(\tau )=p_{\eta }^{\text{limit}}(1-\frac{\gamma (\tau )}{{\gamma }^0})+p_{\eta }^{\text{start}}\frac{\gamma (\tau )}{{\gamma }^0},\text{\hspace{1em}(59)}$$ 其中$p_{\eta }^{\text{start}}$和$p_{\eta }^{\text{limit}}$是超参数，分别给出第一步训练时和无限步骤训练时的前置因子。$\gamma (\tau )$是由式 (56)定义的学习速率。

特定拟合属性$L_{\eta }$的损失函数由数据帧的均方误差(MSE)定义，并且如果$\eta$是一个由原子属性线性组合的帧属性，则经过原子数$N$进行归一化。以能量模型的拟合为例，如式 (37)所示，这些属性$\eta$可以是能量$E$，力$\mathbf{F}$，势能矩阵$\mathbf{\Xi }$，相对能量$\Delta E$⁷²，或它们之间的任意组合，它们的损失函数分别为：

$$L_E(\mathbf{x};\mathbf{\theta })=\frac{1}{N}(E(\mathbf{x};\mathbf{\theta })-E^{\ast }{)}^2,\text{\hspace{1em}(60)}$$ $$L_F(\mathbf{x};\mathbf{\theta })=\frac{1}{3N}\sum _{k=1}^N\sum _{\alpha =1}^3(F_{k,\alpha }(\mathbf{x};\mathbf{\theta })-F_{k,\alpha }^{\ast }{)}^2,\text{\hspace{1em}(61)}$$ $$L_{\Xi }(\mathbf{x};\mathbf{\theta })=\frac{1}{9N}\sum _{\alpha ,\beta =1}^3({\Xi }_{\alpha \beta }(\mathbf{x};\mathbf{\theta })-{\Xi }_{\alpha \beta }^{\ast }{)}^2,\text{\hspace{1em}(62)}$$ $$L_{\Delta E}(\mathbf{x};\mathbf{\theta })=\frac{1}{N}(\Delta E(\mathbf{x};\mathbf{\theta })-{\Delta E}^{\ast }{)}^2,\text{\hspace{1em}(63)}$$ 其中$F_{k,\alpha }$是第$k$个原子上的第$\alpha$个分量的力，上标"$\ast$"表示属性的标签需要提前给出。使用$N$确保每个拟合属性的损失在在其对总损失的贡献之前通过权重对原子贡献进行平均。

如果一部分原子比其他原子更重要，例如，在计算自由能曲线或动力学同位素效应时，某些原子起着重要作用^52,75，则可以使用带有前置因子$q_k$的原子力的均方误差作为损失函数：

$$L_F^p(\mathbf{x};\mathbf{\theta })=\frac{1}{3N}\sum _{k=1}^N\sum _{\alpha }{q}_k(F_{k,\alpha }(\mathbf{x};\mathbf{\theta })-F_{k,\alpha }^{\ast }{)}^2.\text{\hspace{1em}(64)}$$ 具有较大前置因子的原子力将比其他原子中的原子力拟合得更准确。

如果某些力非常大，例如，在高温反应模拟中，力可以大于60 eV/Å^77,78，这时也可以选择让力的损失函数相对于力的大小，而非公式 (61)：

$$L_F^r(\mathbf{x};\mathbf{\theta })=\frac{1}{3N}\sum _{k=1}^N\sum _{\alpha }{\left(\frac{F_{k,\alpha }(\mathbf{x};\mathbf{\theta })-F_{k,\alpha }^{\ast }}{|{\mathbf{F}}_k^{\ast }|+\nu }\right)}^2.\text{\hspace{1em}(65)}$$ 其中 $\nu$ 是一个小常数，用于保护 一个原子，其中 ${\mathbf{F}}_k^{\ast }$的大小很小，从而让$L_F^r$很大。得益于相对力损失函数，大小较小的力可以更精确地拟合。

2.2.3 训练过程

在训练过程中，通过随机梯度下降算法 Adam ¹²⁴最小化损失函数。理想情况下，得到的参数是损失函数最小化时对应的参数，

$${\mathbf{\theta }}^{\ast }=\underset{\mathbf{\theta }}{\mathrm{argmin}}\underset{\tau \to +\infty }{\lim }L(\mathbf{x};\mathbf{\theta },\tau ).\text{\hspace{1em}(66)}$$ 在实践中，Adam优化器在第${\tau }_{\text{stop}}$步停止，并且学习率根据公式 (56) 变化。${\tau }_{\text{stop}}$是超参数，常设置为几百万。

2.2.4 多任务训练

多任务训练过程能够同时处理具有不同特性的数据集，这些特性无法用同一个网络拟合（例如，不同交换相关泛函或不同基组的DFT计算结果的特性）。这些数据集用${\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(n_t)}$表示。对于每个数据集，定义一个训练任务如下：

$$\underset{\mathbf{\theta }}{\min }{L}^{(t)}({\mathbf{x}}^{(t)};{\mathbf{\theta }}^{(t)},\tau ),t=1,\dots ,n_t.\text{\hspace{1em}(67)}$$ 在多任务训练过程中，所有任务共享一个带有可训练参数${\mathbf{\theta }}_d$的描述符，而每个任务都有自己的拟合网络，其可训练参数为${\mathbf{\theta }}_f^{(t)}$，因此${\mathbf{\theta }}^{(t)}=\{{\mathbf{\theta }}_d,{\mathbf{\theta }}_f^{(t)}\}$。在每个训练步骤中，从$1,\dots ,n_t$中随机选择一个任务，然后执行Adam优化器来最小化$L^{(t)}$一步以更新参数${\mathbf{\theta }}^{(t)}$。如果不同拟合网络具有相同的架构，可以共享某些层的参数以提高训练效率。

2.3 模型偏差

模型偏差${ϵ}_y$是由独立训练使用相同数据集初始化的模型集合${\mathcal{M}}_1,\dots ,{\mathcal{M}}_{n_m}$所推断出的属性$\mathbf{y}$的标准偏差。DeePMD-Kit支持将$\mathbf{y}$视为原子力${\mathbf{F}}_i$和虚拟张量$\mathbf{\Xi }$。模型偏差被用来估计模型在特定数据帧上的误差，记为$\mathbf{x}$，包含所有原子的坐标和化学物种。我们介绍原子力和维里张量的模型偏差 $${ϵ}_{\mathbf{F},i}(\mathbf{x})=\sqrt{⟨\Vert {\mathbf{F}}_i(\mathbf{x};{\mathbf{\theta }}_k)-⟨{\mathbf{F}}_i(\mathbf{x};{\mathbf{\theta }}_k)⟩{\Vert }^2⟩},\text{\hspace{1em}(68)}$$ $${ϵ}_{\mathbf{\Xi },\alpha \beta }(\mathbf{x})=\frac{1}{N}\sqrt{⟨({\Xi }_{\alpha \beta }(\mathbf{x};{\mathbf{\theta }}_k)-⟨{\Xi }_{\alpha \beta }(\mathbf{x};{\mathbf{\theta }}_k)⟩{)}^2⟩},\text{\hspace{1em}(69)}$$ 其中${\mathbf{\theta }}_k$是模型${\mathcal{M}}_k$的参数，集合平均$⟨\cdot ⟩$通过以下方式估计：

$$⟨\mathbf{y}(\mathbf{x};{\mathbf{\theta }}_k)⟩=\frac{1}{n_m}\sum _{k=1}^{n_m}\mathbf{y}(\mathbf{x};{\mathbf{\theta }}_k).\text{\hspace{1em}(70)}$$ 小的${ϵ}_{\mathbf{F},i}$意味着模型已经学习了给定的数据；否则，它还没有覆盖，需要扩展训练数据。如果${\mathbf{F}}_i$或$\mathbf{\Xi }$的大小相当大，可以使用相对模型偏差${ϵ}_{\mathbf{F},i,\text{rel}}$或${ϵ}_{\mathbf{\Xi },\alpha \beta ,\text{rel}}$代替绝对模型偏差：⁷⁸

$${ϵ}_{\mathbf{F},i,\text{rel}}(\mathbf{x})=\frac{|{ϵ}_{\mathbf{F},i}(\mathbf{x})|}{|⟨{\mathbf{F}}_i(\mathbf{x};{\mathbf{\theta }}_k)⟩|+\nu },\text{\hspace{1em}(71)}$$ $${ϵ}_{\mathbf{\Xi },\alpha \beta ,\text{rel}}(\mathbf{x})=\frac{{ϵ}_{\mathbf{\Xi },\alpha \beta }(\mathbf{x})}{|⟨\mathbf{\Xi }(\mathbf{x};{\mathbf{\theta }}_k)⟩|+\nu },\text{\hspace{1em}(72)}$$ 其中$\nu$是一个小常数，用于保护由于${\mathbf{F}}_i$或$\mathbf{\Xi }$的大小较小而产生大型模型偏差的原子。

${ϵ}_{\mathbf{F},i}$和${ϵ}_{\mathbf{\Xi },\alpha \beta }$的统计数据可以提供，包括在原子索引$i$和分量索引$\alpha ,\beta$上的最大、平均和最小模型偏差。在原子索引$i$和分量索引$\alpha ,\beta$上分别描述了最大的力模型偏差${ϵ}_{\mathbf{F},\text{max}}$。将力的最大模型偏差${ϵ}_{\mathbf{F},\text{max}}$作为同步学习算法或主动学习算法中最佳的错误指标发现在一个框架中。^103,107

3.技术实现

除了引入了新的强大功能外，DeePMD-kit还设计了以下目标：高性能、高可用性、高扩展性和社区参与度。这些目标对于DeePMD-kit成为在各种计算领域广泛使用的平台至关重要。在本节中，我们将介绍为实现这些目标而采取的几项技术实现。

3.1 代码架构

DeePMD-kit利用TensorFlow的计算图架构构建其DP模型¹²⁵，该模型由多种操作符组成，包括用C++实现的自定义操作符，如环境矩阵、Ewald求和、压缩操作符及其反向传播。TensorFlow提供的自动微分机制用于计算DP模型相对于输入原子坐标和模拟单元张量的导数。为了优化性能，一些关键的自定义操作符使用CUDA或ROCm工具库实现以在GPU上执行。TensorFlow提供的自动微分机制用于计算DP模型相对于输入原子坐标和仿真单元张量的导数。为了优化性能，一些关键的自定义操作符使用CUDA或ROCm工具库实现以在GPU上执行。如图2所示，DeePMD-kit的代码包括以下模块:

图2 DeePMD-kit代码的架构。红色框表示DeePMD-kit包内的模块（绿色框），橙色框表示计算图，蓝色框表示DeePMD-kit的依赖项，黄色框表示与DeePMD-kit集成功能的第三方包（包括LAMMPS，i-PI，GROMACS，AMBER，OpenMM，ABACUS，ASE，MAGUS，DP-Data，DP-GEN和MLatom）。自定义操作符是指TensorFlow不提供的操作符，包括原子环境矩阵、基于成对势能的插值和嵌入矩阵的制表推理。黑色箭头$\mathrm{A}\to \mathrm{B}$表示模块A依赖于模块B。红色和紫色箭头分别表示"定义"和"使用"。

• 核心C++库提供自定义操作符的实现，如原子环境矩阵、邻居列表和压缩神经网络。需要注意的是，核心C++库是在没有TensorFlow的C++接口的情况下独立构建和测试的。

• GPU库（CUDA¹²⁶或ROCm¹²⁷）是核心C++库的可选部分，用于在GPU设备上计算自定义操作符。该库依赖于GPU工具库（NVIDIACUDA Toolkit或AMD ROCm Toolkit），也是独立构建和测试的。

• DP操作符库包含一些TensorFlow不支持的自定义操作符¹²⁵。TensorFlow提供了Python和C++接口来实现一些自定义操作符，其中TensorFlow的C++库被打包在其Python包内。

• "模型定义"模块用Python编写，用于生成由TensorFlow操作符、DP自定义操作符和模型参数组成的计算图。该图可以保存为可以用于推断的文件。它依赖于TensorFlowPythonAPI（版本1，tf.compat.v1）和其他Python依赖项，如NumPy¹²⁸和H5Py¹²⁹包。

• Python应用程序编程接口（API）用于推断，可以从文件中读取计算图，并使用TensorFlowPython API执行该图。

• C++ API建立在TensorFlow C++接口之上，与PythonAPI在推断方面具有相同的功能。

• C API是C++ API的包装，提供与C++ API相同的功能。与C++ API相比，CAPI具有更稳定的应用程序二进制接口（ABI）并确保向后兼容性。

• 仅包含头文件的 C++ API 是 C API 的封装，提供与 C++ API相同的接口。它具有与 C API 相同的稳定ABI，但仍然利用了 C++的灵活性。

• 命令行界面（CLI）提供给普通用户和开发者，并用于训练和推理。它依赖于模型定义模块和Python API。

CMake 构建系统¹³⁰ 管理所有模块，而 pip 和scikit-build¹³¹ 包用于将 DeePMD-kit 分发为 Python 包。标准的 Python 单元测试框架¹³² 用于所有 Python代码的单元测试，而 GoogleTest 软件¹³³ 用于所有 C++代码的测试。 GitHub Actions 自动化构建、测试和部署流程。

3.2 性能

3.2.1 硬件加速

在TensorFlow框架中，静态图将多个运算符与输入和输出组合在一起。在训练或推断过程中，有两种耗时的运算符。第一种是TensorFlow的原生神经网络运算符（见第 2.1.1节）和矩阵运算符，它们已经被TensorFlow框架本身充分优化，适用于CPU和GPU架构¹²⁵。第二，DeePMD-kit的自定义运算符用于计算原子环境（Eq. (6)和(12)），用于与成对势函数的插值，以及用于嵌入矩阵的表格推断（Eq. (30)）。这些运算符不受TensorFlow框架的支持，但可以使用OpenMP¹³⁴、CUDA¹²⁶和ROCm¹²⁷在CPU和GPU上并行加速（附录2列出的不支持GPU的功能除外）。

环境矩阵运算符包含两个步骤¹⁰⁸：邻居列表的格式化和计算矩阵$\mathcal{R}$的元素。在格式化步骤中，原子$i$的邻居根据它们的类型${\alpha }_j$、到原子$i$的距离$r_{ij}$和最后的索引$j$进行排序。为了提高在GPU上的排序性能，原子类型、距离和索引被压缩为一个64位整数$S\in \mathbb{N}$用于排序：

$$S={\alpha }_j\times {10}^{15}+\lfloor r_{ij}\times {10}^8\rfloor \times {10}^5+j.\text{\hspace{1em}(73)}$$ 排序后的邻居索引从排序后的 $S$ 中解压缩，然后用于格式化邻居列表。

3.2.2 MPI 实现多设备训练和分子动力学模拟

用户可以选择利用多个 CPU 核心、GPU或多个节点上的硬件，在训练或分子动力学（MD）模拟期间实现更快的性能和更大的内存。为了方便此操作，DeePMD-kit增加了面向消息传递接口（MPI）的实现^135,136，用于多设备训练和MD 模拟，具体描述如下。

多设备训练使用分布式训练框架Horovod¹³⁷。 Horovod通过将数据批次 $\mathcal{B}$沿批次大小的轴均匀分配给工作器，在数据并行模式下工作¹³⁸。

在训练过程中，每个工作器在不同的偏移量上消耗切片的输入记录，并且只对可训练参数梯度与同行进行平均。该设计避免了批次大小和张量形状的冲突，并减少了进程间需要通信的字节数。mpi4py 包¹³⁹ 用于去除冗余日志。

多设备 MD 模拟通过利用第三方 MD 包的并行功能来实现。 例如，LAMMPS通过优化分区、通信和邻居列表实现 CPU的并行¹⁴⁰。 AMBER 在与DeePMD-kit接口中构建了类似的邻居列表^52,54,141。DeePMD-kit 支持本地原子环境计算，并接受从其他软件中获取的邻居列表$n(i)$，以替代原生的邻居列表计算¹⁰⁸。在一个设备中，来自其他设备的邻居被视为不对该设备的总势能 $E$贡献原子能量 $E_i$ 的"幽灵"原子。

3.2.3 非冯·诺依曼分子动力学（NVNMD）

在 CPU 和 GPU上进行分子动力学（MD）模拟时，很大一部分时间和能量（例如超过 95%）被 DP模型推理所消耗。这个推理过程受到冯·诺依曼（vN）架构的"内存墙"和"功耗墙"瓶颈的限制，这意味着大量时间和能量（例如超过90%）被浪费在处理器和内存之间的数据传输上。因此，提高计算效率变得困难。

为了解决这些挑战，非冯·诺依曼分子动力学（NVNMD）使用非冯·诺依曼（NvN）架构芯片加速推理过程。NvN芯片包含可用于实现 DP 算法的处理和存储单元。在 NvN芯片中，硬件算法完全流水线运行。模型参数在初始化过程中从离芯片存储器加载后存储在芯片上的内存中。因此，避免了两个数据传输的组成部分：（1）从/向离芯片存储器读取/写入中间结果，和（2）在计算过程中从离芯片存储器加载模型参数。结果，DP模型确保了与 NVNMD 相关的高准确性，而 NvN芯片确保了高计算效率。更多详情请参见文献¹¹⁰。

3.3 可用性

3.3.1 文档

随着越来越多的开发工作，DeePMD-kit的功能和参数也在迅速增长。为了解决这个问题，我们引入了Sphinx¹⁴²和Doxygen¹⁴³来管理和生成开发者文档，这些文档是从代码中的文档字符串生成的。我们使用DArgs包（见第 节）自动生成用户输入参数的Sphinx文档。目前，这些文档托管在Read theDocs上（https://docs.deepmodeling.org/projects/deepmd/）。此外，我们也努力使DeePMD-kit引发的错误消息对用户更加清晰。此外，GitHubDiscussion论坛允许用户提问和获取答案。最近，我们还发表了一些教程^49,54，以帮助新用户快速学习DeePMD-kit。

3.3.2 便捷安装

如图2所示，DeePMD-kit依赖于TensorFlow的Python和C++库，这可能会使新用户从源代码构建TensorFlow和DeePMD-kit变得困难且耗时。因此，我们提供了通过pip、Conda（DeepModeling和conda-forge¹⁴⁴渠道）、Docker和离线包在Linux、macOS和Windows平台上分发的预编译二进制包。在这些预编译的二进制包的帮助下，用户可以在短短几分钟内安装DeePMD-kit。这些二进制包包括DeePMD-kit的LAMMPS插件、i-PI驱动程序和GROMACS补丁。由于LAMMPS在其最新版本中提供了插件模式，不必重新编译LAMMPS即可编译DeePMD-kit的LAMMPS插件¹⁴⁰。我们提供了一个编译的二进制包，其中包括C API和仅包含头文件的C++API，使其更容易与诸如AMBER^52,54,141的复杂软件集成。

3.3.3 用户界面

DeePMD-kit提供了命令行界面（CLI）用于训练、冻结和测试模型。除了CLI参数外，用户还必须提供一个完成了第2节中列出的各组件的参数的JSON¹⁴⁵或YAML¹⁴⁶文件。DArgs包（见第3.5节）解析这些参数以检查用户输入是否正确。如文献 所示，我们提供了一个使用用户界面的示例。用户还可以使用DP-GUI（见第3.5节）在交互式网页中填写参数并将其保存到一个JSON¹⁴⁵文件中。

DeePMD-kit提供了一种自动算法，帮助新用户决定一些参数。例如，自动批处理大小$\mathcal{B}$确定训练或推断过程中的最大批处理大小，以充分利用GPU上的内存。自动邻域大小$N_c$通过对训练数据进行描述来确定邻居的最大数量，以减少模型的内存使用量。自动概率确定了在训练过程中使用各体系的概率。这些自动参数减少了学习和使用DeePMD-kit的难度。

3.3.4 输入数据

为了训练和测试模型，用户需要按照指定的格式提供拟合数据。DeePMD-kit支持两种数据输入文件格式：NumPy二进制文件¹²⁸和HDF5文件¹⁴⁷。与文本文件相比，这些格式被设计为在程序使用并行算法读取时具有更高的性能。HDF5文件的优势在于可以在单个文件中存储多个数组，使得在不同机器之间传输更加容易。Python软件包"DP-Data"（参见第 节）可以从电子结构计算软件的输出生成这些文件。

3.3.5 模型可视化

DeePMD-kit支持大部分TensorBoard可视化功能¹²⁵，如跟踪和可视化指标、查看模型图、张量直方图、可训练变量摘要和调试配置文件。

3.4 可扩展性

3.4.1 应用编程接口和第三方软件

DeePMD-kit提供多种API，包括Python、C++、C和仅头文件的C++API，以及命令行接口（CLI），如图2所示。这些API主要用于开发人员和高级用户在不同情境下进行推断。Sphinx¹⁴²生成了文档中的API详细信息。

这些API可以方便地被各种第三方软件使用。 例如，PythonAPI被第三方Python软件包如ASE¹⁴⁸、MAGUS¹⁴⁹和DP-Data（见第3.5节）所使用。 C++, C或者仅包含头文件的C++API还被整合到了几个第三方分子动力学包中，例如LAMMPS^140,150，i-PI¹⁵¹，GROMACS¹⁵²，AMBER^52,54,141，OpenMM^153,154和ABACUS¹⁵⁵。此外，命令行界面(CLI)还被各种第三方工作流程包调用，例如DP-GEN¹⁰⁷和MLatom³⁴。虽然ASE计算器、LAMMPS插件、i-PI驱动和GROMACS补丁是在DeePMD-kit代码内部开发的，但其他部分是独立分发的。通过将这些API整合到他们的程序中，研究人员可以进行模拟和最小化运算，而不受DeePMD-kit软件特性的限制^72,75,83,156。此外，如果需要的话，他们还可以将DP模型与DeePMD-kit包外的其他势能结合起来^52,72,157。

DP模型的主要应用之一是分子动力学模拟，它可以与DeePMD-kit包接口的几个第三方包协同使用，提供广泛的支持功能:

• LAMMPS¹⁴⁰通过DeePMD-kit项目内开发的专用插件与DeePMD-kit无缝集成。此插件支持MPI，如第 3.2.2节中所讨论的，并提供了力计算等基本功能。此外，它还可以在同步学习期间实时计算模型偏差(如[[eq:68,eq:69,eq:70,eq:71,eq:72]]所示)。该插件可以从多种来源获取(38)中的原子和帧特定参数，包括常数、LAMMPS计算得出的电子温度，或来自LAMMPS的任何compute样式。LAMMPS还支持使用Ewald求和和快速算法粒子-粒子粒子-网格Ewald(PPPM)来计算经典点电荷的长程（库仑）相互作用。这些方法中的k空间部分，涉及将高斯电荷分布在傅里叶空间中计算电荷间相互作用¹⁵⁸(如(50)所示)，被DPLR方法用于处理长程相互作用。

• i-PI¹⁵¹通过DeePMD-kit项目内提供的专用驱动程序与DeePMD-kit集成。该驱动程序支持由i-PI引擎驱动的路径积分分子动力学(PIMD)，并与MolSSIDriver Interface(MDI)包¹⁵⁹和ASE包中的类似接口¹⁴⁸兼容。然而，i-PI驱动程序和引擎之间的通信依赖于UNIX域套接字或网络接口，这可能限制性能。为了克服这个限制，开发人员将PIMD功能整合到LAMMPS软件包中，实现与DeePMD-kit的无缝集成。

• 通过定制的源代码，将AMBER¹⁴¹与DeePMD-kit软件包集成在一起⁵⁴。AMBER/DeePMD-kit接口允许使用DPRc模型进行有效的QM/MM+DPRc模拟⁵²。该接口扩展了QM/QM相互作用，包括对QM/MM相互作用的范围修正。DeePMD-kit软件包仅推断所选的QM区域（由AMBER掩码分配）及其在QM区域半径截断内的MM缓冲区。与LAMMPS集成类似，该接口支持MPI，如第 3.2.2节所讨论的，并允许在同步学习期间实时计算模型偏差。AMBER/DeePMD-kit接口还允许执行变能的自由能模拟，利用AMBER的GPU加速自由能引擎¹²⁰以及用于MM转化的新特性^161-162，并使用间接的MM$\to$QM/$\Delta$-MLP方法¹⁶³将末态校正到更高层次。

• 被广泛采用的分子动力学引擎OpenMM¹⁵³通过OpenMM插件与DeePMD-kit集成。该插件允许使用DP模型进行标准分子动力学模拟，并支持混合DP/MM类型的模拟。在混合模拟中，体系可以使用固定的DP区域或在模拟过程中自适应地改变区域¹⁵⁴。

• GROMACS¹⁵²通过对GROMACS进行补丁，与DeePMD-kit集成。该补丁通过分配由DeePMD-kit推断的原子类型，实现了DP/MM模拟。

• ABACUS¹⁵⁵支持DeePMD-kit提供的C和C++接口。此外，ABACUS还支持基于不同方法的各种分子动力学，例如使用LJ势能的经典分子动力学和基于KSDFT、SDFT和OFDFT等方法的第一性原理分子动力学。DeePMD-kit与这些方法结合的可能性需要进一步探索。

这些与现有软件包的整合和接口为研究人员提供了灵活性，使他们能够将DeePMD-kit与其他强大的工具相结合，增强分子动力学模拟的能力。

3.4.2 自定义插件

DeePMD-kit基于面向对象的设计，第[2]节中讨论的每个组件对应一个Python类。这种设计的一个优点是为这些组件提供了插件系统的可用性。有了这个插件系统，开发人员可以创建和整合自定义组件，而无需修改DeePMD-kit包。这种方法加速了他们的想法的实现。此外，插件系统还便于在DeePMD-kit包中添加新的组件。

3.5 DeepModeling社区

DeePMD-kit是一个免费且开源的软件，采用LGPL-3.0许可证，使开发人员能够修改和整合DeePMD-kit到他们自己的软件包中。作为核心，DeePMD-kit于2021年创立了一个名为DeepModeling的开源社区，专门管理科学计算的开源软件包。从那时起，许多与科学计算相关的开源软件包要么被创建，要么加入了DeepModeling社区，如DP-GEN¹⁰⁷，DeePKS-kit¹⁶⁴，DMFF¹⁶⁵,ABACUS¹⁵⁵，DeePH¹⁶⁶，和DeepFlame¹⁶⁷，等等，无论是直接还是间接与DeePMD-kit相关。与DeePMD-kit相关的DeepModeling软件包如下所示。

1. Deep Potential GENerator(DP-GEN)¹⁰⁷是一个实现同步学习过程¹⁰³、能够以最小的人为干预和计算成本生成一致准确的DP模型的软件包。DP-GEN2是该软件包的下一代版本，构建在工作流平台Dflow上。

2. 深度势能·热力学积分（DP-Ti）是一个Python软件包，它使用户能够计算自由能、执行热力学积分，并确定具有DP模型的材料的压力-温度相图。

3. DP-Data是一个Python软件包，帮助用户在不同格式之间转换原子数据，并通过电子计算和MLP软件包计算原子数据。它可以用于生成DeePMD-kit的训练数据文件，并通过3Dmol.js¹⁶⁸可视化结构。该软件包支持插件系统，并与ASE¹⁴⁸兼容，使其能够支持任何数据格式，不受软件包代码的限制。

4. DP-Dispatcher是一个Python软件包，用于生成高性能计算（HPC）调度程序的输入脚本，将其提交到HPC系统，并在完成之前监视其进度。它最初是作为DP-GEN软件包的一部分开发的¹⁰⁷，但现在已成为一个独立的软件包，可以为其他软件包提供服务。

5. DArgs是一个Python软件包，用于管理和过滤用户输入参数。它提供了一个Sphinx¹⁴²扩展，用于生成参数的文档。

6. DP-GUI是一个基于Web的图形用户界面（GUI），采用Vue.js框架¹⁶⁹构建。它允许用户在网页上交互式填写参数，并将其保存到JSON¹⁴⁵文件中。DArgs用于在GUI中提供参数的详细信息和文档。

4. 示例应用：分子动力学

图3 以同步学习方式执行深度势能分子动力学的一般工作流程。

本节介绍了使用同步学习¹⁷⁰从头开始执行深度势能分子动力学的一般工作流程，如图3所示。目标模拟可以包括各种条件，例如温度、压力、经典动力学或路径积分动力学，带有或不带有增强采样方法，在平衡或非平衡态下，以及不同的尺度和时间尺度。需要注意的是，本节不是用户手册或教程，也不深入研究特定的体系。

初始步骤涉及准备初始数据集。该数据集通常通过在与目标模拟相同条件下的小规模、短时间的MD模拟中进行采样生成。模拟级别可以不同，从从头算¹⁷⁰到半经验⁵²或力场⁷⁷，取决于计算成本。随后，使用高精度的从头算方法重新标记这些构型。

准备好初始数据后，下一步是执行同步学习循环，这对于提高目标模拟的准确性至关重要。每个循环包括三个步骤：训练、探索和标记。在训练步骤中，DeePMD-kit使用短训练步骤对现有的目标数据集训练多个模型（通常为四个模型）。这些模型可以从不同的随机种子或前一次迭代训练的模型初始化。在探索步骤中，使用其中一个模型执行目标模拟并对构型空间进行采样。如果目标模拟涉及非平衡过程，则随着同步学习循环的进行，模拟时间可以逐渐增加。从满足以下条件的构型（或构型中的一部分原子以减少计算成本⁷⁷）中随机选择：

$$\{R_n|n\in I_{\text{cand}},I_{\text{cand}}=\{n|{\theta }_{\text{low}}\le {ϵ}_{\mathbf{F},\text{max}}<{\theta }_{\text{high}}\}\}\text{\hspace{1em}(74)}$$ 其中${ϵ}_{\mathbf{F},\text{max}}$在第 [2.3]节中给出，${\theta }_{\text{low}}$应设置为高于目标数据集中的大多数${ϵ}_{\mathbf{F},\text{max}}$的值，${\theta }_{\text{high}}$通常设置为比${\theta }_{\text{low}}$高约0.15 eV/Å的值。这些阈值确保仅选择尚未添加到目标数据集中的构型。使用一致的从头算方法对选定的构型进行标记，并在标记步骤中将其添加到目标数据集中，然后继续下一次迭代。

如果一个模拟中准确构型的比例(${ϵ}_{\mathbf{F},\text{max}}<{\theta }_{\text{low}}$)在同步学习的后续循环中保持不变，则可以认为该模拟为目标模拟，并停止迭代。这样的模拟轨迹可以进一步进行分析。

上述工作流程可以手动执行，也可以使用DP-GEN软件包¹⁰⁷进行自动执行。

5. 基准测试

我们在多个数据集上对使用DeePMD-kit软件包内开发的不同描述符的各种势能模型进行了基准测试，以展示其精度和性能。数据集、模型、硬件和结果将在以下部分描述和讨论。

5.1 数据集

表1 用于基准测试的数据集。

数据集	帧数	元素	DFT水平	参考文献
水	140000	H, O	PBE0+TS/PAW ($E_{\text{cutoff}}$ = 1565 eV)	9,61
铜	15366	Cu	PBE/PAW ($E_{\text{cutoff}}$ = 650 eV)	107
HEA	8160	Ta, Nb, W, Mo, V, Al	PBE/PAW($E_{\text{cutoff}}$ = 1200 eV)	51,171
OC2M	2000000	Ag, Al, As, Au, B, Bi, C, Ca, Cd, Cl, Co, Cr, Cs, Cu, Fe, Ga, Ge, H, Hf, Hg, In, Ir, K, Mg, Mn, Mo, N, Na, Nb, Ni, O, Os, P, Pb, Pd, Pt, Rb, Re, Rh, Ru, S, Sb, Sc, Se, Si, Sn, Sr, Ta, Tc, Te, Ti, Tl, V, W, Y, Zn, Zr	RPBE/PAW ($E_{\text{cutoff}}$ = 350 eV)	115,116
SPICE	1132808	H, Li, C, N, O, F, Na, Mg, P, S, Cl, K, Ca, Br, I	$\omega$B97M-D3BJ/def2-TZVPPD	172
二肽	33850	H, C, N, O, S	$\omega$B97M-D3BJ/def2-TZVPPD	172

我们使用的数据集包括水^9,61，铜（Cu）¹⁰⁷，高熵合金（HEA）^51,171，OpenCatalyst2020（OC20）中的OC2M子集^115,116，小分子/蛋白质相互作用化学能量（SPICE）¹⁷²，以及SPICE中的二肽子集¹⁷²，如表1所示，并列举如下：

• 水数据集包含了从路径积分从头算MD模拟和经典从头算MD模拟中收集的个水的构型。构型使用了混合版本的Perdew-Burke-Ernzerhof（PBE0）¹⁷³+Tkatchenko-Scheffler（TS）函数和投影增强波（PAW）方法进行标记¹⁷⁴。能量截断设置为115 Ry（1565 eV）。

• 铜数据集包含了FC C、HCP 和BCC晶体中的个构型。MD模拟在温度范围为50至 K和压力范围为1至5$\times 10^4$巴内采样了这些构型。采用同步学习方案¹⁷⁰选择了关键构型，提高了用于估计模型预测误差的模型集合的准确性。使用了Perdew-Burke-Ernzerhof（PBE）函数¹⁷⁵和PAW方法，能量截断为650 eV。

• 高熵合金（HEA）数据集包括6种元素：钽、铌、钨、钼、钒和铝^51,171。这些元素占据一个$2\times 2\times 2$的体心立方（BCC）晶格，包含16个以随机方式排列的原子。钽、铌、钨、钼和钒的浓度涵盖了整个组成空间，而铝被视为一种添加剂，其最大数量小于六个。MD模拟在温度范围为50至388.1K和压力范围为1至5$\times 10^4$巴内采样了这些构型。采用同步学习方案¹⁷⁰选择了关键构型，提高了用于估计模型预测误差的模型集合的准确性。该数据集包含由密度泛函理论用PBE近似计算¹⁷⁵的标记构型个。使用了PAW方法，能量截断为1200 eV，k空间采样网格大小为0.12 Å${}^{-1}$。

• OC20数据集中的OC2M子集¹¹⁶将OC20数据集¹¹⁵中的200万个构型取出，包含57个元素。OC20数据集包含了1281040个构型，涵盖了广泛的材料、表面和吸附物，并通过修订后的PBE函数¹⁷⁵在周期性边界条件下进行标记。使用了PAW方法，能量截断为350 eV。

• 小分子/蛋白质相互作用化学能量（SPICE）数据集¹⁷²是一个类似药物的数据集，包括各种子集：二肽、溶剂化氨基酸、PubChem分子、DES370K二聚体、DES370K单体和离子对。数据集由标记为$\omega$B97M-D3BJ/def2-TZVPPD水平的个非周期配置组成^176,177。它包含15种元素并包含带电荷的配置。我们采用了与参考文献中描述的方法相同的方法，将每种唯一的元素与形式电荷的组合视为不同的原子类型。

• SPICE中的二肽子集合¹⁷²包括由20种天然氨基酸及其常见质子化变体形成的所有可能的二肽。该子集包含具有H、C、N、O和S元素的个配置，对应于氨基酸。

以上数据集是具有代表性的，因为它们包含液体、固体和气体，周期性和非周期性边界条件下的配置，跨越一系列温度和压力的配置，以及处于不同质子化状态的离子和药物样分子。在化学物理领域中研究所有这些体系是必不可少的。

我们将所有数据集分成一个包含95%数据的训练集和一个包含剩余5%数据的验证集。

5.2 模型与硬件

我们比较了不同的描述符，包括本地参考框架（loc_frame）、两体嵌入全信息DeepPot-SE（se_e2_a）、具有两体嵌入全信息和径向信息DeepPot-SE的混合描述符（se_e2_a+se_e2_r）、具有两体嵌入全信息和三体嵌入DeepPot-SE的混合描述符（se_e2_a+se_e3），以及基于注意力的描述符（se_atten）。在所有模型中，如果适用，我们将$r_s$设置为0.5 Å，$M_{<}$设置为16，$L_a$设置为2。两体嵌入网络${\mathcal{N}}_{e,2}$使用了(25,50,100)个神经元，三体嵌入网络${\mathcal{N}}_{e,3}$使用了(2,4,8)个神经元，拟合网络${\mathcal{F}}_0$使用了(240,240,240,1)个神经元。对于具有两体嵌入全信息和半径信息的混合描述符DeepPot-SE（se_e2_a+se_e2_r）的全信息部分（se_e2_a）和具有两体全信息和三体DeepPot-SE的混合描述符（se_e2_a+se_e3）的三体部分，我们将$r_c$设置为4 Å；对于OC2M体系，我们将$r_c$设置为9 Å，而在其他情况下，我们将$r_c$设置为6 Å。

我们使用双精度（FP64）和单精度（FP32）的神经网络对每个模型进行了固定步数的训练（步用于水、铜和双肽，步用于HEA，步用于OC2M和SPICE）。我们使用LAMMPS软件包¹⁴⁰进行了水、Cu和HEA的分子动力学模拟，尽可能使用了更多的原子。我们将压缩模型的性能与原始模型进行了比较¹⁰⁸。用于基准测试性能的平台包括128核AMD EPYC 7742、NVIDIA GeForce RTX 3080 Ti(12 GB)、NVIDIA Tesla V100 (40 GB)、NVIDIA Tesla A100 (80 GB)、AMDInstinct MI250，以及Xilinx Virtex Ultrascale+ VU9PFPGA用于NVNMD¹¹⁰。注意，目前，模型压缩功能仅支持se_e2_a、se_e2_r和se_e3描述符，NVNMD仅支持FP64精度中具有不超过4种化学元素的体系的常规se_e2_a。关于se_atten的模型压缩功能正在开发中。

值得注意的是，这些模型旨在在相同条件下比较软件包支持的不同描述符和浮点数精度，以推荐最佳模型使用。然而，需要强调的是，训练步数有限，并且超参数（如神经网络中的神经元数量）未针对任何特定体系进行调整。因此，不建议将这些模型用于生产目的，并且与其他参考文献中报道的成熟模型进行比较是没有意义的^9,51,107,172。

此外，不推荐将这些模型与其他软件包生成的模型进行比较，因为确保不同模型的超参数完全相同，或者使所有模型在不同软件包中消耗相同的计算资源并不容易。

5.3 结果与讨论

表2 水、铜、HEA、OC2M、二肽和SPICE 验证集的能量每原子以及力的平均绝对误差（MAE）和均方根误差（RMSE）。

		loc_frame		se_e2_a		se_e2_a+se_e2_r		se_e2_a+se_e3		se_atten
体系	指标	FP64	FP32	FP64	FP32	FP64	FP32	FP64	FP32	FP64	FP32
Water	E RMSE	0.7	0.7	1.0	1.0	0.9	1.0	1.0	1.0	1.5	1.2
	F RMSE	40.0	39.2	49.0	48.4	48.6	50.0	46.5	45.9	44.4	42.3
Cu	E RMSE	12.7	19.2	3.0	2.8	4.8	4.9	2.5	2.6	3.2	3.6
	F RMSE	84.7	105	17.7	17.9	21.4	22.0	16.8	16.6	16.9	16.9
HEA	E RMSE	⋯	⋯	15.4	15.3	13.5	14.5	12.1	17.2	5.5	6.4
	F RMSE	⋯	⋯	134	137	163	158	136	180	90.7	98.3
OC2M	E RMSE	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	15.1	14.3
	F RMSE	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	155	148
Dipeptides	E RMSE	⋯	⋯	9.5	9.6	12.5	11.7	14.9	16.8	12.8	12.7
	F RMSE	⋯	⋯	97.9	97.7	98.6	101	160	222	99.7	96.7
SPICE	E RMSE	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	80.9	78.3
	F RMSE	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	233	234

表3 水、Cu、HEA、OC2M、二肽和SPICE体系的训练性能（每步毫秒）。“FP64”表示双精度浮点数，“FP32”表示单精度浮点数，“FP64c”和“FP32c”分别表示双精度和单精度压缩训练 。“EPYC”在128个AMD EPYC 7742内核上执行，“3080 Ti”在NVIDIA GeForce RTX 3080Ti显卡上执行，“V100”在NVIDIA Tesla V100显卡上执行，“A100”在NVIDIA Tesla A100显卡上执行，“MI250”在AMD Instinct MI250图形计算单元（GCD）上执行。

		loc_frame		se_e2_a				se_e2_a+se_e2_r				se_e2_a+se_e3				se_atten
体系	硬件	FP64	FP32	FP64	FP32	FP64c	FP32c	FP64	FP32	FP64c	FP32c	FP64	FP32	FP64c	FP32c	FP64	FP32
Water	EPYC	14.7	9.20	97.3	45.0	28.4	16.2	63.7	32.5	29.9	15.4	141	85.2	34.0	20.6	1210	383
	3080 Ti	7.00	4.80	24.6	10.3	9.70	6.40	26.3	11.6	12.0	8.20	52.8	17.2	16.3	6.80	199	26.9
	V100	7.90	8.50	11.1	8.20	5.90	4.80	13.6	10.9	6.90	6.40	23.5	14.0	8.60	7.30	69.6	31.7
	A100	10.7	10.0	8.20	9.30	4.90	5.70	14.5	10.8	7.80	6.30	24.5	12.0	7.50	7.20	30.8	21.2
	MI250	11.7	10.9	20.3	13.1	7.70	7.00	27.3	19.7	11.5	10.9	278	27.7	12.8	11.2	125	31.7
Cu	EPYC	4.90	3.30	33.7	12.8	8.00	5.40	19.9	10.0	10.5	5.30	45.5	24.2	9.10	6.50	226	89.1
	3080 Ti	3.20	2.20	6.50	5.10	4.60	3.90	8.70	6.30	5.90	3.40	11.8	4.80	7.20	5.70	36.8	8.80
	V100	3.20	3.80	4.20	4.80	3.20	3.70	6.50	5.30	5.50	4.10	7.90	5.60	6.00	5.80	15.6	11.9
	A100	4.00	3.90	3.80	3.70	3.10	3.00	5.40	5.30	4.10	4.10	8.00	5.60	4.80	4.60	11.6	11.2
	MI250	4.80	4.90	6.90	6.40	5.10	5.00	9.10	9.40	7.40	7.00	49.9	10.1	8.00	7.30	23.6	18.6
HEA	EPYC	⋯	⋯	53.4	30.5	19.4	12.2	52.3	29.3	27.7	16.7	83.7	51.1	26.6	15.7	159	60.1
	3080 Ti	⋯	⋯	38.4	25.2	11.2	9.10	71.4	41.8	16.3	12.7	93.6	41.0	19.7	15.0	35.9	9.10
	V100	⋯	⋯	33.2	29.8	11.8	11.1	63.2	47.4	17.5	16.5	65.5	49.6	27.4	18.7	15.6	11.9
	A100	⋯	⋯	30.5	28.6	10.9	10.4	51.6	67.4	16.9	21.2	61.7	52.9	18.6	18.8	11.7	11.5
	MI250	⋯	⋯	48.8	42.7	18.5	18.0	72.3	69.3	28.7	27.3	134	88.4	32.7	32.3	21.6	19.5
OC2M	EPYC	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	2070	625
	3080 Ti	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	352	46.0
	V100	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	120	52.8
	A100	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	51.4	30.9
	MI250	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	171	55.7
Dipeptides	EPYC	⋯	⋯	49.7	30.5	21.2	19.4	52.0	35.3	30.1	21.2	89.5	61.1	35.0	21.2	214	91.5
	3080 Ti	⋯	⋯	54.8	39.5	17.3	11.3	90.0	64.3	19.0	15.3	131	67.7	25.4	19.2	26.1	12.0
	V100	⋯	⋯	54.1	52.6	14.8	14.8	88.0	84.3	20.5	21.7	96.2	103	30.1	30.8	14.3	10.6
	A100	⋯	⋯	50.2	50.8	14.3	14.3	89.0	75.9	20.7	19.9	91.1	82.7	26.6	26.7	13.2	11.1
	MI250	⋯	⋯	66.2	67.8	23.1	22.9	117	112	35.0	32.4	155	129	45.9	44.9	19.6	16.8
SPICE	EPYC	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	244	98.0
	3080 Ti	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	35.4	15.3
	V100	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	17.3	15.9
	A100	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	11.9	12.2
	MI250	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	29.0	24.1

表4 水、Cu和HEA体系的MD性能（μs/step/atom）。“FP64”表示双精度浮点，“FP32”表示单精度浮点，而“FP64c” 和“FP32c” 分别表示双精度和单精度浮点的压缩模型 。“EPYC”表示在128 核的AMD EPYC 7742 测试，“3080 Ti” 表示在NVIDIA GeForce RTX 3080 Ti 卡测试，“V100” 表示在NVIDIA Tesla V100 卡测试，“A100” 表示在NVIDIA Tesla A100 卡测试，“MI250” 表示在AMD Instinct MI250 Graphics Compute Die (GCD)测试，而“VU9P”表示在Xilinx Virtex Ultrascale+ VU9P FPGA 运行NVNMD 。

		loc_frame		se_e2_a				se_e2_a+se_e2_r				se_e2_a+se_e3				se_atten
体系	硬件	FP64	FP32	FP64	FP32	FP64c	FP32c	FP64	FP32	FP64c	FP32c	FP64	FP32	FP64c	FP32c	FP64	FP32
Water	EPYC	1.25	0.699	19.3	8.73	3.89	2.61	8.33	3.43	3.78	1.86	37.2	15.1	5.04	3.63	221	83.8
	3080 Ti	12.9	8.63	29.0	4.21	9.71	1.73	20.8	3.43	9.06	1.99	69.5	10.5	18.5	2.89	294	32.3
	V100	16.1	16.8	8.25	4.59	1.94	1.51	6.21	3.53	2.22	1.62	22.2	11.3	3.31	2.41	91.2	37.2
	A100	35.7	33.9	4.37	3.01	1.56	1.42	4.11	2.44	2.07	1.53	12.5	7.17	2.64	2.25	35.6	22.4
	MI250	40.2	39.6	7.74	3.96	1.74	1.41	6.03	3.20	2.00	1.54	30.5	18.8	3.51	2.64	55.0	30.2
	VU9P	⋯	⋯	0.306	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
Cu	EPYC	1.14	0.702	22.2	9.38	3.43	2.04	11.9	5.28	3.09	1.56	47.9	19.5	4.20	2.73	200	62.1
	3080 Ti	14.9	8.98	30.5	4.18	8.52	1.51	18.8	3.15	7.98	1.81	74.6	11.2	14.7	2.32	294	33.0
	V100	15.7	15.7	8.73	4.81	1.56	1.27	5.71	3.18	1.84	1.38	24.3	12.2	2.60	1.83	91.1	37.3
	A100	36.9	36.9	4.41	2.65	1.36	1.15	3.35	2.15	1.63	1.42	13.5	7.49	2.15	1.78	36.2	21.0
	MI250	39.0	39.1	8.27	4.13	1.37	1.21	5.62	2.98	1.59	1.35	26.9	12.6	2.56	2.00	55.4	29.5
	VU9P	⋯	⋯	0.310	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
HEA	EPYC	⋯	⋯	32.8	13.0	7.04	4.58	15.3	7.64	6.83	3.80	81.0	33.4	8.56	5.68	156	45.9
	3080 Ti	⋯	⋯	65.3	9.72	10.5	2.51	36.1	6.83	11.9	3.24	171	24.9	29.6	5.37	290	32.8
	V100	⋯	⋯	20.1	10.9	2.88	2.39	12.3	6.86	12.3	2.85	55.2	28.4	9.42	5.47	91.2	37.4
	A100	⋯	⋯	10.4	6.09	2.13	1.83	7.25	5.48	2.98	2.83	30.1	17.1	4.21	4.22	35.0	20.0
	MI250	⋯	⋯	20.1	11.6	4.57	4.22	16.2	12.0	7.01	6.44	76.0	44.9	9.09	7.61	55.7	30.5

表5 GPU卡对于水、铜和高熵合金体系可以模拟的最大原子数量（$10^3$）。“FP64”表示双精度浮点精度，“FP32”表示单精度浮点精度，“FP64c”和“FP32c”分别表示压缩模型 中的双精度和单精度浮点精度。“3080 Ti”表示在NVIDIA GeForce RTX 3080 Ti卡（12GB）上进行测试，“V100”表示在NVIDIA Tesla V100卡（40GB）上进行测试，“A100”表示在NVIDIA Tesla A100卡（80GB）上进行测试。

		se_e2_a				se_e2_a+se_e2_r				se_e2_a+se_e3				se_atten
体系	硬件	FP64	FP32	FP64c	FP32c	FP64	FP32	FP64c	FP32c	FP64	FP32	FP64c	FP32c	FP64	FP32
Water	3080 Ti	27	51	127	141	44	74	94	128	10	21	86	95	3	7
	V100	73	135	415	493	114	196	274	430	28	55	250	265	9	19
	A100	189	332	987	1128	288	488	651	900	76	147	618	736	22	49
Cu	3080 Ti	18	35	214	202	40	77	125	151	7	14	83	144	3	7
	V100	54	106	606	635	122	183	337	461	19	39	271	330	9	19
	A100b	141	244	1534	1615	286	453	706	1074	50	99	697	867	22	49
HEA	3080 Ti	11	19	60	69	21	33	47	59	5	9	42	52	3	7
	V100	30	53	175	184	57	87	131	166	13	25	117	142	9	19
	A100	76	132	447	468	140	218	323	408	35	66	292	365	22	48

^a MI250 卡上的结果没有报告，因为当原子数增加时，ROCm Toolkit 报告了一个段错误。

^b 使用了两个MPI进程，因为当操作符中的元素数量超过 $2^{31}$ 时，TensorFlow发生了整数溢出。

我们在表2中展示了不同模型的验证误差，以及表3和表4中的各个平台上的训练和分子动力学性能，同时在表5中展示了平台可以模拟的最大原子数。在所有数据集中，没有任何一个模型在准确性方面胜过其他模型。对于水体系而言，非平滑的局部坐标描述符实现了最佳的准确性，能量RMSE为0.689meV/原子，力RMSE为39.2meV/Å。此外，该模型在CPU上的计算性能是所有模型中最快的，尽管尚未实现在GPU上（见附录2）。局部坐标描述符在某些情况下准确性较高，但存在限制，限制了其广泛应用性。其中之一是它不是平滑的。此外，该描述符对铜体系的性能表现不佳，而铜体系在温度和压力范围内收集了大量数据¹⁰⁷。另一个限制是它要求所有体系具有相似的化学物种来构建局部坐标，这使得在HEA、OC2M、二肽和SPICE等数据集中应用变得困难。

另一方面，DeepPot-SE描述符在准确性和性能方面具有更好的泛化能力。压缩模型在训练和推理方面比原模型快1倍到10倍，并且NVNMD比常规分子动力学快50倍到100倍，这两者都展示了令人印象深刻的计算性能。预计未压缩的水se_e2_aFP64模型（单个V100显卡上的每步/原子耗时为8.25 $\mu$s）的分子动力学性能接近于文献 中报告的分子动力学性能（每个V100显卡上的每步/原子耗时为8.19 $\mu$s），而压缩模型（每步/原子耗时为1.94 $\mu$s）在该情况下约快3倍。此外，与未压缩模型相比，压缩模型可以在单个显卡上模拟6倍的原子数。三体嵌入描述符在理论上比两体嵌入描述符包含更多信息，更准确但速度较慢。尽管对于水和铜体系而言是正确的，但在HEA和二肽数据集中观察不到预期的准确性顺序。需要进一步研究来确定这种差异的原因，但很可能是因为更多化学种类导致更多可以训练的参数时，损失不能在相同的训练步骤内收敛。此外，这两个数据集上的性能会随着神经网络的增多而降低。

具有类型嵌入的注意力模型在HEA体系中具有更好的准确性，在二肽体系中具有相同的准确性。这些模型在GPU上训练速度更快，并通过减少神经网络数量在这两种体系上达到相同的准确性。但是，在CPU或分子动力学模拟中并未观察到此优势，因为注意力层的计算开销很大，这需要进一步改进。此外，当化学物种很多时，注意力描述符所需的CPU或GPU内存比其他模型少，因为它具有更少的神经网络。这个特性使得它能够应用于具有超过60个化学元素的OC2M数据集和大约20个化学元素的SPICE数据集。

值得注意的是，在几乎所有体系中，FP32的速度比FP64快0.5倍到2倍，并且展示了类似的验证误差。在这种情况下，由于我们只将FP32应用于神经网络而在其他组件中保持精度，模型精度也被称为"混合精度"。这个结果与其他软件包中广泛采用混合精度的事实相一致^38,140。因此，FP32应该被广泛应用于大多数应用程序。此外，FP32可以在性能差的硬件上实现高性能，例如消费级GPU或CPU。

6.总结

DeePMD-kit是一个功能强大且多用途的社区开发的开源软件包，用于使用机器学习势能（MLP）进行分子动力学（MD）模拟。其出色的性能、易用性和可扩展性使其成为各个领域研究人员的热门选择。DeePMD-kit采用LGPL-3.0许可证，任何人都可以自由使用、修改和扩展该软件。

由于其良好设计的代码架构，DeePMD-kit具有高度可定制性，可以在各个方面轻松扩展。模型以面向对象的设计方式组织为Python模块，并保存为计算图，使添加新模型变得更容易。计算图由TensorFlow和自定义操作符组成，使得针对特定硬件架构和特定操作符优化软件包变得更容易。该软件包还具有丰富灵活的API，使其更容易与其他分子模拟软件包集成。DeePMD-kit对计算科学研究人员的贡献是开放的，我们希望社区在未来继续开发和增强其功能。

代码可用性

DeePMD-kit托管在GitHub库中，地址为https://github.com/deepmodeling/deepmd-kit。本研究使用的数据集、模型、模拟体系和基准测试脚本可以从GitHub库https://github.com/deepmodeling-activity/deepmd-kit-v2-paper下载。本研究支持结果的其他数据可在合理请求的情况下由通讯作者提供。

致谢

作者感谢Yihao Liu、Xinzijian Liu、Haidi Wang、HailinYang和GitHub用户ZhengdQin对DeePMD-kit的代码贡献，以及Zhi X. Chen、JincaiYang和Tong Zhu对本文的建议，Ming Li对突出显示图片的设计。D.T.感谢Stefano Baroni、Riccardo Bertossa、Federico Grasselli和PaoloPegolo对本工作的启发性讨论。ChatGPT被用于以"润色英语"指令润色语言，其输出结果经过人工审核。J.Z.和D.M.Y.的工作得到了美国国立卫生研究院的支持（向D.M.Y.资助，基金编号为GM107485）和美国国家科学基金会的支持（基金编号为2209718）。J.Z.对美国新泽西州立大学化学与化学生物学系Van Dyke奖项表示感谢。Y.C.、YifanLi和R.C.的工作得到由美国能源部资助的"溶液和界面化学"（CSI）中心资助（基金编号为DE-SC0019394）。M.R.的工作得到VEGA项目资助（基金编号为1/0640/20），并得到斯洛伐克研究和发展署根据合同编号APVV-19-0371的支持。Q.Z.的工作得到湖南省科技创新计划的支持（基金编号为2021RC4026）。S.L.B.的工作得到挪威研究委员会的支持，通过卓越中心Hylleraas量子分子科学中心（基金编号为262695）。C.L.和R.W.的工作得到美国能源部（DOE）的支持（基金编号为DE-SC0019759）。H.W.的工作得到国家重点研发计划的支持（基金编号为2022YFA1004300），以及国家自然科学基金的支持（基金编号为12122103）。计算资源由以下提供： 深势科技的Bohrium云平台； 新泽西州立大学Office ofAdvanced Research Computing（OARC）；高级网络基础设施协调生态系统：服务和支持 (ACCESS)程序，该程序由美国国家科学基金会提供资助，编号为#2138259、#2138286、#2138307、#2137603和#2138296(通过CHE190067分配给SDSC的超级计算机Expanse)；德克萨斯大学奥斯汀分校的德克萨斯先进计算中心(TACC)，网址：http://www.tacc.utexas.edu(通过CHE20002分配给Frontera超级计算机)； AMD公司的 AMD 云平台；普林斯顿大学的普林斯顿研究计算资源，该资源由普林斯顿计算科学与工程研究所(PICSciE) 和信息技术办公室的研究计算领导的一组联合体组成。

附录1：五阶多项式

定义一个分段定义的函数 $f(x)$，其中 $f(x)$ 在区间 $[0,1)$内是一个五次多项式，在其他区间是一个常数: $$f(x)=\{\begin{array}{ll}1,& x<0\\ a{x}^5+b{x}^4+c{x}^3+d{x}^2+ex+f,& 0\le x<1,\\ 0,& x\ge 1.\end{array}\text{\hspace{1em}(75)}$$令 $f(x)$ 以及它的一阶导数 $f^{\prime }(x)$ 和二阶导数 $f^{\prime \prime }(x)$ 在$x=0$ 和 $x=1$ 处连续， $$\{\begin{array}{l}{\lim }_{x\to 0^{-}}f(x)={\lim }_{x\to 0^{+}}f(x),\\ {\lim }_{x\to 0^{-}}{f}^{\prime }(x)={\lim }_{x\to 0^{+}}{f}^{\prime }(x),\\ {\lim }_{x\to 0^{-}}{f}^{\prime \prime }(x)={\lim }_{x\to 0^{+}}{f}^{\prime \prime }(x),\\ {\lim }_{x\to 1^{-}}f(x)={\lim }_{x\to 1^{+}}f(x),\\ {\lim }_{x\to 1^{-}}{f}^{\prime }(x)={\lim }_{x\to 1^{+}}{f}^{\prime }(x),\\ {\lim }_{x\to 1^{-}}{f}^{\prime \prime }(x)={\lim }_{x\to 1^{+}}{f}^{\prime \prime }(x).\end{array}\text{\hspace{1em}(76)}$$ 求解方程组 (76)，得到解：

$$a=-6,b=15,c=-10,d=0,e=0,f=1.\text{\hspace{1em}(77)}$$$f(x)$ 的最终形式为： $$f(x)=\{\begin{array}{ll}1,& x<0\\ -6x^5+15x^4-10x^3+1,& 0\le x<1,\\ 0,& x\ge 1.\end{array}\text{\hspace{1em}(78)}$$

附录2：GPU不支持的功能

目前，以下功能不支持 GPU：

• 公式 (5)、(6)、(7)、(8)、(9) 和 (10) 中的局部框架描述符。

• 公式 (52)、(53)、(54) 和(55) 中的成对势。

• 从给定数据中计算截断半径内的最大邻居数 $N_c$。

• 公式 (68)、(69)、(70)、(71) 和 (72) 中的模型偏差。

• 所有 NVNMD 相关的功能。

• LAMMPS 插件中 DPLR 的 KSpace 求解器。