用于水电解质溶液的短程机器学习潜力

时间:2026年5月16日
来源:ChemPhysChem

编辑推荐:

摘要 机器学习势(MLPs)扩展了原子模拟的时间和长度尺度,使得能够研究复杂的系统,例如电解质溶液。然而,大多数模型在准确性、计算成本和捕捉长程相互作用的能力之间面临着权衡。大型基础模型虽然具有普遍性,但往往伴随着巨大的开销和高能量需求。相比之下,紧凑的、特定于系统的模型可能为

广告
   X   

摘要

机器学习势(MLPs)扩展了原子模拟的时间和长度尺度,使得能够研究复杂的系统,例如电解质溶液。然而,大多数模型在准确性、计算成本和捕捉长程相互作用的能力之间面临着权衡。大型基础模型虽然具有普遍性,但往往伴随着巨大的开销和高能量需求。相比之下,紧凑的、特定于系统的模型可能为大规模模拟提供更可持续的路径。在这里,我们在水氯化钠(NaCl)溶液上对MACE架构进行了基准测试,系统地改变了模型大小和等变性的水平,以评估它们对准确性、稳定性和效率的影响。我们发现,所研究的MLPs的预测准确性对这里考虑的关键物理观测量影响很小,但对稳定性至关重要,这突显了最小化、专用模型在高效模拟电解质溶液方面的潜力。

图解摘要

机器学习势使得模拟能够超越从头算方法的范围。我们在水氯化钠(NaCl)上对小型MACE模型进行了基准测试,改变其大小和等变性,以探测准确性、稳定性和效率。模型准确性对观测量影响不大,但对稳定性至关重要,这突显了最小化、专用模型在大规模模拟方面的潜力。

1 引言

自从2007年Behler和Parrinello引入高维神经网络势(HDNNPs)以来[1],已经出现了许多不同的架构,包括高斯近似势[2]、神经网络[1, 3]和消息传递神经网络(MPNNs)[4-6]。这导致了模拟时间尺度和长度尺度的突破,现在可以高精度地研究许多复杂系统,例如界面和/或含有电解质的溶液(参见参考文献[7-14])。大多数MLPs的一个核心假设是原子环境的局部性:每个原子都由基于截断的描述符来描述,其能量是根据其局部环境预测的。虽然这保持了可扩展性,但它限制了长程相互作用或电荷转移过程的描述。为了解决这个问题,开发了一些明确包含长程静电学的模型,例如通过电荷平衡[15]、基于Wannier的描述符[16]、时空谱图神经网络[17]或潜在的Ewald求和方案[18]。这些扩展在长程相互作用至关重要的情况下提高了准确性,例如气液界面或水的蒸气相[19, 20],但不可避免地增加了计算成本。MPNNs,如MACE[21]或NequIP[4],提供了一个自然的折中方案。通过在多个层之间传递消息,它们隐式地扩展了有效相互作用范围,同时保持了计算效率[5, 22]。这些模型的典型截断半径约为6 Å,并且有一到三个消息传递层[23],已经能够捕捉液体中大部分相关相互作用,其中分子间距离大约为3 Å[24]。特别是在大量水电解质溶液中,长程静电作用被溶剂分子强烈屏蔽[25],这使得短程MPNNs成为建模这类系统的有吸引力的候选者。然而,除了准确性之外,成本仍然是MLPs的一个核心挑战。一方面,随着模型复杂性的增加,评估能量和力的成本变得更加昂贵。另一方面,生成参考训练数据通常是一个主要瓶颈,因为它需要大量的电子结构计算。为了缓解特定应用领域中的训练数据瓶颈,最近开发了大型基础模型。这些模型旨在提供一个通用力场,覆盖广泛的化学和结构空间,通常涵盖大部分元素周期表。最近的例子包括Meta的UMA[26]模型、CHGNet[27]或MACE-MP-0[28]基础模型等[29-34]。然后可以针对特定系统对这些通用模型进行微调[35]。然而,由于这些模型通常涵盖非常多样化的数据集,它们需要大量的自由参数。因此,它们也增加了计算成本和内存需求,这可能会阻碍使用MLP的主要动机,即能够快速运行大规模模拟。特别是在图形处理单元(GPU)资源有限的情况下,这可能会限制可访问的系统规模和模拟时间尺度。虽然也有人在努力从大型基础模型中提取更小的、专用的模型[36, 37],但我们在这里采取了另一种策略,即训练显式的、特定于系统的MACE模型。我们通过分析数值误差指标来研究小型、专用模型描述水电解质的能力,并检查拟合准确性如何转化为计算出的静态和动态性质的可靠预测。因此,这项研究将为未来应用如何有效选择模型提供指导。我们选择了MACE[21]作为最近的一个MLP,它包括等变的高阶消息传递,显示出出色的数据效率和准确性[23]。MACE使用原子簇展开(ACE)[38]来构建消息。在ACE中,采用了能量的多体展开。在应用所谓的密度技巧后,可以利用张量积,使得评估与体阶和系统大小成线性关系。因此,更高阶的描述符变得可行,与传统的MLP描述符(如原子位置的平滑重叠(SOAP)[39]或以原子为中心的对称函数[1])不同,后者只支持最高到三阶的体阶。这导致了更高的准确性和更高效的数据拟合,在原始ACE公式中实现了线性回归[38]。MACE采用上述方法来构建更高阶的消息。这些消息被嵌入到类似于NequIP[4]的等变消息传递方案中。引入等变性已被证明可以提高拟合的准确性和数据效率[4, 5]。此外,一个重要的超参数是消息传递中未耦合特征通道的数量()。在MACE中,这些通道对应于在体阶耦合之前的独立ACE基函数,因此直接决定了潜在特征空间的大小和模型的表示能力。在这项研究中,我们系统地研究了紧凑的、特定于系统的MACE模型在准确性和计算效率之间的平衡。通过改变未耦合通道的数量(定义模型大小)和等变性水平(0表示不变模型,1表示等变模型),我们探讨了架构选择如何影响预测准确性、外推行为和稳定性。使用氯化钠(NaCl)溶液作为特征明确的基准系统,我们评估了没有显式长程校正的小型MACE模型是否可以以相对较低的计算成本可靠地捕捉静态和动态物理性质。

2 结果和讨论

2.1 MLP的构建

构建MLP的第一步是生成高质量参考数据,充分覆盖相关的配置空间。因此,对基础电子结构方法进行基准测试是必不可少的。然而,这会产生一个矛盾,因为许多感兴趣的物理属性无法通过从头算方法直接获得,因为需要长时间的采样。因此,MLPs也可以作为基准测试电子结构方法的工具(参见参考文献[13])。对于水NaCl的参考方法来说,主要挑战是准确捕捉离子和水相互作用,这已经在几项研究中得到了解决[13, 40-42]。这样的基准测试不是我们研究的重点。因此,我们选择了revPBE0-D3泛函,它已被证明可以可靠地再现大量水系统的实验性质[43, 44]。我们得到的数据集大约包含40,000个结构,相当于20 ps的采样时间。我们选择了800 Ry的平面波截断来进行采样,以使模拟成为可能。然而,我们强调这远未收敛。实际上,即使在几千Ry的情况下,我们也未能得到收敛的钠力,这在计算上也是不可行的(见支持信息(SI))。因此,我们改用了CP2K中的高斯增强平面波(GAPW)[45, 46]方法[47]。因此,使用收敛的平面波截断重新计算了初始轨迹的一个子集,以确保离子力的准确性(见方法部分)。收敛平面波截断对于实现稳定的基于MLP的模拟至关重要。

2.2 数值验证

我们使用上述数据集训练了几个MACE势,并进行了数值验证。我们系统地改变了等变性水平()和未耦合通道的数量(),因为这两个超参数显著影响了训练和模拟的计算成本[23]。我们使用参考文献[48]中提出的方法计算了均方根误差(RMSE)来进行力向量的评估。图1显示了不同模型的结果误差指标。我们在SI中包括了常用的分量力RMSE,以便更好地比较。以下,我们使用与参考文献[23]中相同的术语:例如,32–0表示具有未耦合通道的不变模型。图1显示了训练有素的MLPs在不同和组合下的数值验证结果,我们使用与参考文献[23]中相同的术语:。上面板显示了能量的RMSE,下面板显示了所研究的水NaCl系统的力的元素级和平均RMSE。一般来说,增加复杂性会导致能量和力预测的一致性改进。引入等变性可以将误差减少大约30%–50%,这与Batzner等人的发现一致[4]。相比之下,单独增加只会带来较小的收益。能量误差趋于平稳,而力误差甚至在增加后还会增加,这表明模型的容量已经被超出,进一步增加其复杂性会导致过拟合。我们还注意到,最精确的不变模型64–0仍然无法胜过最小的等变模型8–1。此外,报告的误差与当前的最先进水平相符;参见例如参考文献[23]。正如引言中已经概述的,选择尽可能小的模型的动机是为了降低成本,从而能够运行大规模模拟。本节的最后一段讨论了训练和评估MLP的伴随计算成本。

2.3 MLP的稳定性和外推分析

我们使用获得的MLPs来运行模拟并分析轨迹的稳定性。因此,对于每个MLP,我们在最初的100 ps平衡后运行了五个1 ns的轨迹。图2显示了从不同初始配置开始的平均成功模拟时间比例。作为不稳定性的标准,我们选择了一个原子在单个积分步骤中移动超过一个盒长度的非物理情况。图2显示了不同架构复杂度的水NaCl模型的模拟稳定性。热图显示了每个模型成功模拟的五个1 ns轨迹的平均比例及其标准差,模型按照未耦合通道的数量()和最大等变性水平()进行排列。颜色强度表示稳定性,绿色表示完整运行。数值表示五次运行的平均值和标准差。当原子在一步中移动超过一个盒长度时,模拟将终止,这将导致非物理结构。我们观察到所有等变MLP都为这个模拟长度产生了稳定的模拟。相比之下,具有的不变模型表现出相当大的不稳定性,尽管最小的8–0模型平均仍能运行超过3 ns。随着未耦合通道数量的减少,提供不稳定轨迹的倾向增加,这可能是由于力预测不准确造成的,首先导致配置具有非物理上过短的原子间距离,最终由于陡峭的排斥势导致极大的力。或者,这些不稳定性可能反映了模型在外推行为上的差异或训练数据中的限制。为了更好地理解其原因,我们研究了不变模型中观察到的失败是由于在训练期间相关配置空间的覆盖不足,还是由于力预测中的外推误差引起的,后者引发了一系列非物理配置。理解这种区别很重要,因为迭代重新训练以修复不稳定性可能非常昂贵,而简单地使用更大的模型也同样昂贵。因此,知道失败是由于外推还是其他效应造成的对于选择更强大的模型至关重要。我们使用了参考文献[49]中定义的不同配置之间的SOAP向量核作为相似性度量,

(1)其中是新配置的SOAP向量,是训练数据集中的第个配置,是一个正整数,我们使用了与参考文献[49]中相同的值(详见计算方法)。我们计算了一个外推分数,定义为

(2)

以便具有最高相似性的训练配置决定分数。此外,我们还使用密度泛函理论(DFT)参考方法计算了一个1 ns轨迹的能量和力以及MACE模型的结果误差。结果展示在图3中,我们还报告了与在训练数据相同水平上计算的参考点相比的误差。图3可以在图形查看器中打开。

我们对每个模型的代表性轨迹进行了外推分析,所有轨迹都从相同的初始配置开始。每个模型都是独立平衡的,导致轨迹不相关。展示的是外推得分(方程(2),顶部面板)、能量差异(第二面板)、原子力误差的中位数大小(第三面板),以及相对于DFT的原子力误差的最大值(底部面板),适用于所有模型架构。不变模型的结果()显示在左列,等变模型的结果显示在右列。我们观察到,对于所有模型,外推得分随时间增加,这与其他轨迹中的观察结果一致(见补充信息)。相比之下,不同模型架构的预测准确性有所不同。不变的32–0和128–0模型表现出的能量偏差明显超过了它们的测试集误差(分别为0.15和0.13 meV/原子;见图1),达到了超过10 meV/原子的值。尽管中位数力误差最初保持稳定,但最大原子力误差急剧增加,表明局部力预测可能变得不物理。这些错误的力通过局部环境传播,最终导致中位数力误差的增加。重要的是,仅凭外推得分无法捕捉到这种失败模式。一些模型的能量误差基本保持不变,这表明训练数据集已经足够收敛,观察到的偏差主要反映了架构差异。特别是,不变模型似乎更容易出现局部力误差,这解释了它们的不稳定性。相比之下,等变模型不仅通常显示出更低的误差,与它们的测试集性能一致,而且也显示出更大的稳定性;除了小的16–1模型的能量有所偏差外,能量误差保持在3 meV/原子以下,并没有进一步增加。这表明等变模型可能从力预测的扰动中恢复得更好。在接下来的部分,我们将评估这些问题如何在物理可观测量中表现出来,因为这最终决定了模型的可靠性。

2.4 物理验证

为了评估拟合精度如何转化为物理上有意义的属性,我们首先分析了静态结构可观测量——具体来说是成对的径向分布函数(RDFs)。图4展示了每个离子周围的溶剂化壳层的RDFs。作为参考,我们还包括了来自完整轨迹的DFT数据。虽然这是使用未收敛的平面波截止计算的,但O’Neill等人的先前工作[13]已经表明,结构可观测量(如RDFs)不受影响。相比之下,所有MLPs都是在使用完全收敛的截止计算的数据上训练的,以获得稳定的势能。图4可以在图形查看器中打开。

图中显示了使用不同架构的MLPs获得的钠和氯溶剂化壳层的成对RDFs。为了比较,还包括了来自底层DFT训练数据集(灰色)和实验参考[50]的第一个峰位的RDFs(虚线垂直线)。我们观察到,尽管拟合精度不同,但从不同模型架构获得的RDFs彼此非常吻合。此外,它们与原始数据集的DFT参考数据也非常匹配。这表明拟合误差不会显著影响静态属性(如RDFs)的准确性。此外,RDF峰值位置与实验数据[50]非常接近,支持了底层DFT泛函在捕捉基本溶剂化结构方面的可靠性。此外,我们在图5中展示了基于成对阳离子-阴离子RDF的平均力(PMF),为系统的长程行为提供了洞察。从不同模型架构获得的PMF也非常接近,偏差通常低于0.2 kcal/mol。一个值得注意的例外是8–0模型,它预测的接触离子对(CIP)最小值较浅,过渡到溶剂分离的离子对(SSIP)的障碍较高。然而,由于上述不稳定性,这种模型只能生成目标的5 ns轨迹中的3.15 ns,因此观察到的偏差可能部分反映了采样不足而不是模型本身的不准确性。此外,我们获得了与O’Neill等人使用相同DFT参考方法计算的PMF的定性一致。其余差异可能源于使用了不同的盐浓度。PMF还突出了离子-离子相互作用的长程性质:NaCl势能在大约8–10 Å处收敛到零,表明在水溶液中有效屏蔽。为了明确评估长程静电的作用,我们还使用在包含更大模拟盒子的增强数据集上训练的32–1模型进行了模拟,这些模拟包括了和使用或不使用显式长程校正(基于潜在的Ewald求和[18])。得到的PMFs与使用相应短程模型获得的PMFs几乎相同,无论是在接触离子对最小值的深度还是在更大分离时的渐进衰减方面。在离子-离子RDFs(Na–Na和Cl–Cl)中也观察到了一致的行为,这些在支持信息中报告(见图S30),并且显示了短程和长程模拟的相同长程衰减。这些发现与经典分子动力学(MD)和从头算分子动力学(AIMD)研究[13, 51, 52]一致,支持了短程MPNNs用于模拟稀薄电解质溶液的适用性。

2.4 物理验证

为了评估拟合精度如何转化为物理上有意义的属性,我们首先分析了静态结构可观测量——具体来说是成对的径向分布函数(RDFs)。图4展示了每个离子周围的溶剂化壳层的RDFs。作为参考,我们还包括了来自完整轨迹的DFT数据。虽然这是使用未收敛的平面波截止计算的,但O’Neill等人的先前工作[13]已经表明,结构可观测量(如RDFs)不受影响。相比之下,所有MLPs都是在使用完全收敛的截止计算的数据上训练的,以获得稳定的势能。图4可以在图形查看器中打开。

图中显示了使用不同架构的MLPs获得的钠和氯溶剂化壳层的成对RDFs。为了比较,还包括了来自底层DFT训练数据集(灰色)和实验参考[50]的第一个峰位的RDFs(虚线垂直线)。我们观察到,尽管拟合精度不同,但从不同模型架构获得的RDFs彼此非常吻合。此外,它们与原始数据集的DFT参考数据也非常匹配。这表明拟合误差不会显著影响静态属性(如RDFs)的准确性。此外,RDF峰值位置与实验数据[50]非常接近,支持了底层DFT泛函在捕捉基本溶剂化结构方面的可靠性。此外,我们在图5中展示了平均力(PMF),它基于成对阳离子-阴离子RDF(有关PMF计算的详细信息见补充信息),并提供了对系统长程行为的洞察。从不同模型架构获得的PMFs也非常接近,偏差通常低于0.2 kcal/mol。一个值得注意的例外是8–0模型,它预测的接触离子对(CIP)最小值较浅,过渡到溶剂分离的离子对(SSIP)的障碍较高。然而,由于上述不稳定性,这种模型只能生成5 ns轨迹中的3.15 ns,因此观察到的偏差可能部分反映了采样不足而不是模型本身的不准确性。此外,我们获得了与O’Neill等人报告的PMF的定性一致。剩余的差异可能源于使用了不同的盐浓度。PMF还突出了离子-离子相互作用的长程性质:NaCl势能在大约8–10 Å处收敛到零,表明在水溶液中有效屏蔽。为了明确评估长程静电的作用,我们还使用在包含更大模拟盒子的增强数据集上训练的32–1模型进行了模拟,这些模拟包括了使用潜在的Ewald求和的显式长程校正。得到的PMFs与使用相应短程模型获得的PMFs几乎相同,无论是在接触离子对最小值的深度还是在更大分离时的渐进衰减方面。在离子-离子RDFs(Na–Na和Cl–Cl)中也观察到了一致的行为,这些在支持信息中报告(见图S30),并显示了短程和长程模拟的相同长程衰减。这些发现与经典分子动力学(MD)和从头算分子动力学(AIMD)研究[13, 51, 52]一致,支持了短程MPNNs用于模拟稀薄电解质溶液的适用性。

2.4 物理验证

为了评估拟合精度如何转化为物理上有意义的属性,我们首先分析了静态结构可观测量——特别是成对的径向分布函数(RDFs)。图4展示了每个离子周围的溶剂化壳层的RDFs。作为参考,我们还包括了来自完整轨迹的DFT数据。虽然这是使用未收敛的平面波截止计算的,但O’Neill等人的先前工作[13]已经表明,结构可观测量(如RDFs)不受影响。相比之下,所有MLPs都是在使用完全收敛的截止计算的数据上训练的,以获得稳定的势能。图4可以在图形查看器中打开。

图中显示了使用不同架构的MLPs获得的钠和氯溶剂化壳层的成对RDFs。为了比较,还包括了来自底层DFT训练数据集(灰色)和实验参考[50]的第一个峰位的RDFs(虚线垂直线)。我们观察到,尽管拟合精度不同,但从不同模型架构获得的RDFs彼此非常吻合。此外,它们与原始数据集的DFT参考数据也非常匹配。这表明拟合误差不会显著影响静态属性(如RDFs)的准确性。此外,RDF峰值位置与实验数据[50]非常接近,支持了底层DFT泛函在捕捉基本溶剂化结构方面的可靠性。此外,我们在图5中展示了平均力(PMF),它基于成对阳离子-阴离子RDF(有关PMF计算的详细信息见补充信息),并提供了对系统长程行为的洞察。从不同模型架构获得的PMF也非常接近,偏差通常低于0.2 kcal/mol。一个值得注意的例外是8–0模型,它预测的接触离子对(CIP)最小值较浅,过渡到溶剂分离的离子对(SSIP)的障碍较高。然而,由于上述不稳定性,这种模型只能生成目标的5 ns轨迹中的3.15 ns,因此观察到的偏差可能部分反映了采样不足而不是模型本身的不准确性。此外,我们获得了与O’Neill等人使用的相同DFT参考方法计算的PMF的定性一致。其余差异可能源于使用了不同的盐浓度。PMF还突出了离子-离子相互作用的长程性质:NaCl势能在大约8–10 Å处收敛到零,表明在水溶液中有效屏蔽。为了明确评估长程静电的作用,我们还使用在包含更大模拟盒子的增强数据集上训练的32–1模型进行了模拟,这些模拟包括了和使用或不使用潜在的Ewald求和的显式长程校正。结果得到的PMFs与使用相应短程模型获得的PMFs几乎相同,无论是在接触离子对最小值的深度还是在更大分离时的渐进衰减方面。在离子-离子RDFs(Na–Na和Cl–Cl)中也观察到了一致的行为,这些在支持信息中报告(见图S30),并显示了短程和长程模拟的相同长程衰减。这些发现与经典分子动力学(MD)和从头算分子动力学(AIMD)研究[13, 51, 52]一致,支持了短程MPNNs用于模拟稀薄电解质溶液的适用性。

2.4 物理验证

为了评估拟合精度如何转化为物理上有意义的属性,我们首先分析了静态结构可观测量——具体来说是成对的径向分布函数(RDFs)。图4展示了每个离子周围的溶剂化壳层的RDFs。作为参考,我们还包括了来自完整轨迹的DFT数据。虽然这是使用未收敛的平面波截止计算的,但O’Neill等人的先前工作[13]已经表明,结构可观测量(如RDFs)不受影响。相比之下,所有MLPs都是在使用完全收敛的截止计算的数据上训练的,以获得稳定的势能。图4可以在图形查看器中打开。

图中显示了使用不同架构的MLPs获得的钠和氯溶剂化壳层的成对RDFs。为了比较,还包括了来自底层DFT训练数据集(灰色)和实验参考[50]的第一个峰位的RDFs(虚线垂直线)。我们观察到,尽管拟合精度不同,但从不同模型架构获得的RDFs彼此非常吻合。此外,它们与原始数据集的DFT参考数据也非常匹配。这表明拟合误差不会显著影响静态属性(如RDFs)的准确性。此外,RDF峰值位置与实验数据[50]非常接近,支持了底层DFT泛函在捕捉基本溶剂化结构方面的可靠性。此外,我们在图5中展示了平均力(PMF),它基于成对阳离子-阴离子RDF(有关PMF计算的详细信息见补充信息),并提供了对系统长程行为的洞察。从不同模型架构获得的PMF也非常接近,偏差通常低于0.2 kcal/mol。一个值得注意的例外是8–0模型,它预测的接触离子对(CIP)最小值较浅,过渡到溶剂分离的离子对(SSIP)的障碍较高。然而,由于上述不稳定性,这种模型只能生成目标的5 ns轨迹中的3.15 ns,因此观察到的偏差可能部分反映了采样不足而不是模型本身的不准确性。此外,我们获得了与O’Neill等人使用的相同DFT参考方法计算的PMF的定性一致。其余差异可能源于使用了不同的盐浓度。PMF还突出了离子-离子相互作用的长程性质:NaCl势能在大约8–10 Å处收敛到零,表明在水溶液中有效屏蔽。为了明确评估长程静电的作用,我们还使用在包含更大模拟盒子的增强数据集上训练的32–1模型进行了模拟,这些模拟包括了和使用或不使用显式长程校正(基于潜在的Ewald求和[18])。结果得到的PMFs与使用相应短程模型获得的PMFs几乎相同,无论是在接触离子对最小值的深度还是在更大分离时的渐进衰减方面。在离子-离子RDFs(Na–Na和Cl–Cl)中也观察到了一致的行为,这些在支持信息中报告(见图S30),并显示了短程和长程模拟的相同长程衰减。这些发现与经典分子动力学(MD)和从头算分子动力学(AIMD)研究[13, 51, 52]一致,支持了短程MPNNs用于模拟稀薄电解质溶液的适用性。

2.4 物理验证

为了评估拟合精度如何转化为物理上有意义的属性,我们首先分析了静态结构可观测量——具体来说是成对的径向分布函数(RDFs)。图4展示了每个离子周围的溶剂化壳层的RDFs。作为参考,我们还包括了来自完整轨迹的DFT数据。虽然这是使用未收敛的平面波截止计算的,但O’Neill等人的先前工作[13]已经表明,结构可观测量(如RDFs)不受影响。相比之下,所有MLPs都是在使用完全收敛的截止计算的数据上训练的,以获得稳定的势能。图4可以在图形查看器中打开。

图中显示了使用不同架构的MLPs获得的钠和氯溶剂化壳层的成对RDFs。为了比较,还包括了来自底层DFT训练数据集(灰色)和实验参考[50]的第一个峰位的RDFs(虚线垂直线)。我们观察到,尽管拟合精度不同,但从不同模型架构获得的RDFs彼此非常吻合。此外,它们与原始数据集的DFT参考数据也非常匹配。这表明拟合误差不会显著影响静态属性(如RDFs)的准确性。此外,RDF峰值位置与实验数据[50]非常接近,支持了底层DFT泛函在捕捉基本溶剂化结构方面的可靠性。此外,我们在图5中展示了平均力(PMF),它基于成对阳离子-阴离子RDF(有关PMF计算的详细信息见补充信息),并提供了对系统长程行为的洞察。从不同模型架构获得的PMF也非常接近,偏差通常低于0.2 kcal/mol。一个值得注意的例外是8–0模型,它预测的接触离子对(CIP)最小值较浅,过渡到溶剂分离的离子对(SSIP)的障碍较高。然而,由于上述不稳定性,这种模型只能生成目标的5 ns轨迹中的3.15 ns,因此观察到的偏差可能部分反映了采样不足而不是模型本身的不准确性。此外,我们获得了与O’Neill等人使用的相同DFT参考方法计算的PMF的定性一致。其余差异可能源于使用了不同的盐浓度。PMF还突出了离子-离子相互作用的长程性质:NaCl势能在大约8–10 Å处收敛到零,表明在水溶液中有效屏蔽。为了明确评估长程静电的作用,我们还使用在包含更大模拟盒子的增强数据集上训练的32–1模型进行了模拟,这些模拟包括了和使用或不使用显式长程校正(基于潜在的Ewald求和[18])。结果得到的PMFs与使用相应短程模型获得的PMFs几乎相同,无论是在接触离子对最小值的深度还是在更大分离时的渐进衰减方面。在离子-离子RDFs(Na–Na和Cl–Cl)中也观察到了一致的行为,这些在支持信息中报告(见图S30),并显示了短程和长程模拟的相同长程衰减。这些发现与经典分子动力学(MD)和从头算分子动力学(AIMD)研究[13, 51, 52]一致,支持了短程MPNNs用于模拟稀薄电解质溶液的适用性。

2.4 物理验证

为了评估拟合精度如何转化为物理上有意义的属性,我们首先分析了静态结构可观测量——具体来说是成对的径向分布函数(RDFs)。图4展示了每个离子周围的溶剂化壳层的RDFs。作为参考,我们还包括了来自完整轨迹的DFT数据。虽然这是使用未收敛的平面波截止计算的,但O’Neill等人的先前工作[13]已经表明,结构可观测量(如RDFs)不受影响。相比之下,所有MLPs都是在使用完全收敛的截止计算的数据上训练的,以获得稳定的势能。图4可以在图形查看器中打开。

图中显示了使用不同架构的MLPs获得的钠和氯溶剂化壳层的成对RDFs。为了比较,还包括了来自底层DFT训练数据集(灰色)和实验参考[50]的第一个峰位的RDFs(虚线垂直线)。我们观察到,尽管拟合精度不同,但从不同模型架构获得的RDFs彼此非常吻合。此外,它们与原始数据集的DFT参考数据也非常匹配。这表明拟合误差不会显著影响静态属性(如RDFs)的准确性。此外,RDF峰值位置与实验数据[50]非常接近,支持了底层DFT泛函在捕捉基本溶剂化结构方面的可靠性。此外,我们在图5中展示了平均力(PMF),它基于成对阳离子-阴离子RDF(有关PMF计算的详细信息见补充为此,我们从经典分子动力学(MD)模拟中选取了包含一个或两个NaCl离子对的较大模拟盒子,这些模拟盒子中分别含有146个或149个水分子。随后,根据上述相同的第一性原理密度泛函(DFT)方法对这些配置进行了重新计算,以获得能量和力数据。有关结构生成协议的更多细节,请参见补充信息(SI)。

4.2 多层感知器(MLPs)的构建
我们构建了训练数据集,其中包含1500个含有单个NaCl离子对的结构、750个含有两个离子对的结构,以及64个水分子。立方模拟盒子的边长分别为12.484 Å(NaCl溶液浓度较低时)和12.553 Å(浓度较高时)。此外,还通过沿正常模式扫描生成了32个单独的水分子结构,以增强MLP的稳定性。在验证过程中使用了5%的训练数据。测试集包含750个单离子对结构和375个双离子对结构,这些结构在训练过程中未被使用。我们总共训练了10个MACE [21] 0.3.5版本模型。截断半径被设置为某个特定值(这里未提及具体值),两个层次在信息传递过程中遍历了总范围(此处未提及具体范围)。所有其他超参数在模型架构间保持不变,并在补充信息中提供。为了评估长程相互作用的影响,我们使用了从经典分子动力学模拟中得到的增强数据集训练了两个额外的模型。这个增强数据集包括90个含有一个离子对的结构,盒子边长为16.742 Å,以及84个含有两个离子对的结构,盒子边长为16.674 Å。测试数据集与用于短程模型的数据集相同。我们使用MACE 0.3.14版本训练了两个模型:一个是短程模型,另一个模型通过潜在Ewald求和法[18]加入了显式的长程电静力修正。训练程序和所有其他超参数均与短程模型使用的相同。

4.3 MLP模拟
使用经过训练的MLP进行的模拟主要通过MD引擎LAMMPS [66]进行,该引擎集成了MACE。时间步长为0.5 fs,采用Nosé–Hoover恒温器在300 K下的正则(NVT)系综中运行模拟。恒温器的时间常数设置为50 fs。模拟中使用了包含512个水分子和八个离子对的盒子,盒子尺寸为24.968 Å,并在周期性边界条件下进行。从测试数据集中随机抽取了五个初始配置,并进行了另外100 ps的平衡处理。所有的生产运行持续时间为1 ns。对于长程相互作用的研究,模拟使用了Atomic Simulation Environment(ASE)[67, 68]。这些模拟采用了相同的起始配置、时间步长和目标温度,但使用了具有摩擦系数的Langevin恒温器(摩擦系数未提及)。

4.4 外推分析
在外推分析中,我们使用了mlptrain [49]包的一个改编版本。每个配置通过周期性逐个原子的SOAP描述符进行描述,该描述符的截断距离为6 Å,径向基函数的最大阶数为某个值(未提及具体值),径向基函数的数量为某个值(未提及具体值)。根据参考文献[49],核函数是通过某个指数计算的。

4.5 硬件规格
我们使用NVIDIA A100 GPU进行模型训练,每个GPU配备了80 GB内存,单精度计算性能高达19.5 TFLOPS。

致谢
我们感谢德国研究基金会(DFG)在德国卓越战略(EXC 2089/1-390776260,e-conversion项目)下的财政支持。同时,我们也感谢德国研究基金会(DFG)CRC 1487项目“Iron, upgraded!”(项目编号443703006)的资助。作者们还对莱布尼兹超级计算中心(ww.lrz.de)提供的计算和数据资源表示感谢。L.H. 感谢“化学工业基金”通过Kekulé奖学金提供的财政支持。开放获取项目的资金支持和组织工作由Project DEAL完成。

本研究得到了德国研究基金会(Grant 2089/1-390776260和443703006)的支持。

利益冲突声明
作者声明没有利益冲突。

数据可用性声明
训练数据、测试数据、CP2K输入文件和轨迹快照可作为Zenodo仓库获取:https://doi.org/10.5281/zenodo.18108882。

生物通微信公众号
微信
新浪微博


生物通 版权所有