CQPES：基于GPU加速的置换不变多项式神经网络全维精确势能面构建软件包

《Chemistry》：CQPES: A GPU-Aided Software Package for Developing Full-Dimensional Accurate Potential Energy Surfaces by Permutation-Invariant-Polynomial Neural Network

本工作报道了CQPES v1.0软件包，它通过整合置换不变多项式（PIP）与神经网络（NN），利用GPU（TensorFlow）加速训练和Levenberg–Marquardt优化器，实现了全维精确势能面（PES）的高效自动化构建。该软件提供与主流动力学软件（如Gaussian 16、Polyrate）的直接接口，为分子动力学与光谱学研究提供了强大工具。

　　

摘要

精确的势能面（PESs）是分子动力学和光谱学精确研究的先决条件。置换不变多项式神经网络（PIP-NN）方法已被证明在构建气相分子系统的全维PESs方面非常成功。基于十多年的发展，我们推出了CQPES v1.0（重庆势能面），这是一个旨在自动化和加速PES构建的开源软件包。CQPES将数据准备、PIP基生成和模型训练集成到一个现代化的基于Python的工作流中，同时保留了用于处理动力学接口的高效Fortran内核。其关键特性包括通过TensorFlow进行GPU加速训练、用于高精度拟合的鲁棒Levenberg–Marquardt优化器、通过Jupyter和Tensorboard进行实时监控，以及建立在此基础上的主动学习模块。我们通过四个代表性案例研究展示了CQPES的能力：CH₄用于基准测试高对称性处理，CH₃CN用于典型的单分子异构化反应，OH + CH₃OH用于测试大型系统上的GPU训练加速，以及Ar + H₂O用于验证主动学习模块。此外，CQPES提供了与成熟动力学软件（如Gaussian 16、Polyrate 2017-C、VENUS96C、RPMDRate v2.0和Caracal v1.1）的直接接口，使其能够立即应用于化学动力学和动力学模拟。

1. 引言

势能面（PESs）架起了电子结构和核运动之间的桥梁。它们植根于Born–Oppenheimer近似，将电子能量表示为核坐标的函数。虽然直接从头算分子动力学（“on-the-fly”方法）很流行且被广泛使用，但对于复杂系统的长时间模拟，由于电子结构计算的高成本，往往难以实现。与其实时计算能量，不如将PES表示为一个数学函数更为高效。然而，并没有一种万全之策可以表示所有分子系统的PES。因此，PES通常通过拟合大量从头算数据点来构建。这需要计算大量构型的电子能量。获取拟合数据是一项计算密集型任务，因为计算成本随系统规模和理论水平迅速增加。此外，将这些离散数据点转换为连续的高维势能景观是另一个挑战。尽管存在这些障碍，但回报是丰厚的。一旦建立，高质量的PES就成为一个计算平台，使得能够以远低于从头算计算成本的方式广泛研究光谱、动力学和动态，将以前不可能的模拟变为常规工作。

近年来，在开发中小型分子的全维化学精确PESs方面取得了显著进展。它们的电子能量通常在耦合簇和组态相互作用理论水平上确定。Bowman等人最近综述了这些方法20年发展的全面视角。已经提出了各种方法，常见的例子包括原始使用的置换不变多项式（PIPs）、神经网络（NNs）、作为NN输入的PIPs（PIP-NN）、非冗余PIPs（NR-PIP-NN）、多体PIPs（MB-PIP-NN）以及作为NN输入的基本不变量多项式（FI-NN）。此外，基于核的回归方法（例如，KRR和GPR）是提供统计不确定性的高效非参数策略。

其中，将不变描述符与神经网络相结合的方法尤其强大。它们将基于物理的对称性与机器学习的灵活性相结合。描述符保证了关于相同原子的平移、旋转和置换的不变性，而NN则作为通用逼近器来重现数据点。该方法成功地拟合了多达约10个原子的系统，使用了高质量数据（10³~10⁶个构型，CCSD(T)或MRCI水平）。更大的系统难以处理，因为PIP的数量呈阶乘增长。Hao等人的最新工作释放了FIs的能力，发现可以为超过20个原子的系统（即阿司匹林C₉H₈O₄）生成描述符。除此之外，高效的采样策略至关重要。低水平从头算分子动力学或在原始PESs上的动力学对于高维采样非常有效。为了进一步降低成本，已经提出了几何相似性、能量截断和主动采样方法来过滤掉冗余数据点并收集信息丰富的构型。

目前有几个可用的软件包旨在自动化构建PESs，例如PES-Learn、ROBOSURFER、AUTOSURF和Asparagus。这些软件包提供了有用的工具来自动化从头算数据生成和拟合。然而，对于涉及柔性原子的高维反应系统，计算成本仍然很高。随着系统规模和数据集增长，在CPU上训练高保真神经网络成为一项繁重的任务。现有解决方案侧重于自动化，但通常缺乏用于训练的原生GPU支持。此外，提供与反应动力学软件的直接、高效接口至关重要。这定义了我们旨在填补的特定空白。

在这项工作中，我们报告了CQPES，它建立在我们课题组过去12年积累的PIP-NN开发基础上。底层代码已在许多分子系统中得到验证，这保证了软件的可靠性。CQPES将这些经过验证的方法集成到一个现代化的软件包中。它具有Python的灵活性、GPU加速的快速训练以及通过Fortran内核确保效率的与动力学软件的直接接口。

2. PIP和PIP-NN理论

2.1. PIP

笛卡尔坐标缺乏整体平移、旋转和置换不变性。因此，改用距离矩阵坐标，其元素r_ij表示原子i和j之间的核间距离。通过取距离矩阵上三角部分的非零元素并将其排列成向量，得到距离向量r，其长度为k = N(N-1)/2，其中N是分子系统中的原子数。

在PIP方法中，核间距离首先转换为类Morse变量y_ij= exp[-r_ij/α]，这可以保证当r_ij增大时的渐近行为。对于距离向量r，它变为y = exp(-r/α)，其中α是可调范围参数。PIP的显式形式为G = ? ∏_{i< />}^Ny_ij^l_ij，其中?是包含相同原子间所有置换操作的对称化算子，l_ij是y_ij的度数，∑_{i< />}^Nl_ij是每个单项式的总度数。

在PIP拟合中，势能V_PIP表示为V_PIP= c^TG = ∑_ic_iG_i，其中向量c包含所有系数，上标T表示转置。基于充分采样的构型（几何结构和能量），可以通过线性最小二乘优化确定系数，从而获得PIP PES。

2.2. PIP-NN

在PIP-NN拟合中，从几何结构到势能的映射以NN形式描述。输入PIPs（G）和输出势能（V）都进行归一化，使其缩放到[-1, 1]的范围。通常，具有2个隐藏层的PIP-NN估计的势能V_PIP-NN由公式（5）给出，其中W^(i,j)是权重，b^(i,j)是偏置，在相邻层i和j（j = i+1, i=0,1,2）之间建立连接。在输出层之外，在层0、1和2之间使用激活函数f。常用的激活函数包括tanh函数和sigmoid函数。此外，Shang和Zhang报道使用反平方根单位（ISRU）激活函数f(t) = t / √(t²+1) 在保持高拟合精度的同时提供了略快的速度。

通常，对于回归任务，选择均方误差（MSE）作为损失函数，其通用公式为L_MSE= (1/N_data) ∑_i=1^N_data(T_i^true- T_i^pred)²。均方根误差（RMSE）用于评估训练模型的性能，定义为L_RMSE= √(L_MSE)。

对于NN优化，即确定最优参数W和b，已经开发了许多算法。对于回归任务，LM（Levenberg–Marquardt）和L-BFGS（有限内存Broyden–Fletcher–Goldfarb–Shanno算法）二阶优化方法在数学上比一阶方法更精确和鲁棒。LM和L-BFGS方法旨在使用近似的Hessian矩阵来加速NN训练中的反向传播步骤。准牛顿二阶优化器（包括LM和L-BFGS）更新NN参数的一般公式为θ_k+1= θ_k- H^-1?L。

对于LM方法，更新公式为θ_k+1= θ_k- (J^TJ + μI)^-1J^Te，其中J是包含NN误差对NN参数的一阶导数的雅可比矩阵，e表示第k次循环的NN误差，μ是自适应阻尼因子，I是单位矩阵。

在NN训练开始时，初始参数θ₀（W₀和b₀）是随机设置的。对于回归任务，由Glorot和Bengio提出的Xavier初始化方案有利于快速收敛和较小的最终损失，因为它可以防止训练过程中的梯度消失问题，其中权重w应服从该区间的均匀分布。NN的参数在训练循环（也称为epoch）中不断更新。经过几百个epoch的训练后，NN的损失不再减少或达到稳定状态，从而获得收敛且最优的NN。

如果关注某些特定区域（例如，静态点和反应路径周围的低能区域），可以通过为这些区域的构型分配较大权重来拟合PES。通过这种加权策略，高能量（或围绕某些特定区域）的数据点被分配比低能量数据（或其他区域）更低的权重。例如，对于CH₄，能量低于2.0 eV（V_THR）的数据点权重设置为1.0，而能量高于2.0 eV的数据点权重设置为0.1。V_THR的选择由目标属性指导；这里，2.0 eV覆盖了相关的振动态。通过为高能区域（如接近解离或排斥壁）分配较低权重，迫使模型专注于平衡点的精度，这对于可靠的振动频率计算至关重要。

虽然这种硬截断方法通过明确分离势能阱对简单的单分子系统效果很好，但对于复杂反应，加权函数可能更合适以避免不连续性。另一个广泛使用的权重函数是W_i^V= Δ / (Δ + V_i)，其中V_i是第i个采样几何结构的势能，Δ是一个可修改的参数，用于调整NN偏好的能量范围。另一个权重函数被证明可以提高ANI型NNs（精确神经网络分子能量引擎）在静态点附近区域的精度。这些权重然后通过加权MSE损失函数影响NN训练：L = (1/N_data) ∑_i=1^N_dataW_i^V(T_i^true- T_i^pred)²。因此，加权PIP-NN模型在权重较高的区域表现出改进的性能。有趣的是，注意PIP的第一个元素总是1。因此方程（4）可以重写为V = c^Tp = c₁+ ∑_i=2^kc_ip_i= WA + b，这意味着PIP可以被视为没有隐藏层的特殊PIP-NN。

3. 探索构型空间

构建PES模型的目的是让研究人员专注于分子动力学和光谱学，而不必担心获取势能的成本。然而，采样和量子化学计算是构建PES中最基本和最耗时的步骤。在量子化学计算中，计算成本受三个因素影响：电子结构方法和相应的基组、系统规模和采样构型的规模。因此，对于具有确定数量采样结构的特定系统，选择一种平衡精度和计算成本的合适电子结构方法至关重要。全面的讨论超出了本工作的范围，但感兴趣的读者可参考文献[36]。需要注意的是，目标PES的精度和有效范围取决于PES的预期用途。PES模型可以是全局的或局部的，全维的或降维的，反应的或非反应的。以下讨论侧重于全局全维反应PESs，它们代表了最复杂的情况。其他情况则更简单。

数据采样的第一步是确定采样构型空间的范围。对于多通道反应，通常足以覆盖研究目标所规定的可及反应通道，并且可能不需要包括具有较高势垒的反应通道。此外，对于表现出长程相互作用的系统，必须特别注意这些区域。第二步是在反应构型空间内采样结构。有几种方法可用，包括从头算直接动力学、简正模式采样等。在这里，我们根据经验介绍一些实用的采样结构方法。反应PES的关键区域包括反应物、复合物、过渡态和产物——统称为静态点。对于每个静态点，通过使用自适应调整的因子随机扭曲优化的平衡结构，生成几百个新结构以涵盖所有维度。对于最小能量路径（MEP）区域，将此方法应用于沿内禀反应坐标（IRC）的每个结构以覆盖MEP区域。对于其他区域，通常采用两种方法。第一种是进行广泛的刚性扫描，这非常实用，因为它允许灵活确定网格密度和扫描范围。建议最初用稀疏网格采样相对较宽的区域。第二种方法是使用低水平电子结构方法执行从头算直接动力学。然后可以通过拟合所有采样结构及其相应的势能来构建初步PES。这个初步PES通常是粗糙的。分子动力学模拟对于探索整个动力学构型空间和识别“空洞”（PES上缺乏足够结构的区域）特别有用。如果轨迹落入“空洞”，由于该区域结构表示不足，PES将给出显著负的势能。为了改进初步PES，应将围绕这些“空洞”的额外结构纳入采样数据集。这个改进过程可能需要重复多次，直到PES良好收敛。

对于中小型系统，可以使用上述过程生成全局、全维且精确的反应PES。然而，对于相对较大的系统，使用上述过程构建高保真全维PES的计算成本变得过高。一个有前途的解决方案是基于PIP-NN的增量机器学习（?-ML）方法。

基于PIP-NN方法的?-ML方法的实现可以分为六个关键步骤：（i）选择低水平和高水平计算电子结构方法；（ii）使用低水平方法构建良好收敛的PES；（iii）选择低水平数据集的一个子集来获取低水平和高水平方法之间的能量差；（iv）将选定的数据集和相应的能量差拟合到一个校正势；（v）应用校正PES来获得剩余数据集的能量差，从而将能量差与低水平能量相加来估计高水平能量；（vi）拟合具有高水平能量的整个数据集以产生最终的高水平PES。关键要注意第三步是最关键的。为了增强这个选择过程的效率，我们采用了Behler开发的方法，该方法利用了NNs的有限外推和非线性拟合能力。该方法利用了这样一个观察结果：几个PIP-NN PESs，尽管共享相同的NN拟合参数，但在采样不良的区域表现不同。因此，应基于初始数据集拟合多个初步校正PESs。通常，使用3~5个初步校正PESs来确定低水平数据集的未选择部分。然后使用这些校正PESs之间的平均能量差来确定一个构型是否应添加到校正数据集中。较大的平均能量差表明该构型应包含在校正数据集中。这个选择过程也应迭代重复，直到校正PES良好