《Journal of Chemical Theory and Computation》：Beyond Partitioning: Using Force Field Science to Evaluate Electrostatics Models
编辑推荐：
本文通过力场科学（force field science）和机器学习（machine learning）方法系统评估了多种静电模型，提出从对称性适应微扰理论（SAPT）的二聚体相互作用能出发直接训练模型，可获得比传统单体静电势（ESP）拟合更精确的静电和诱导能预测，为精准分子模拟（如药物发现、生物分子研究）提供了新路径。
引言：生物系统中的静电与诱导相互作用
单原子离子或水分子与生物大分子的相互作用，对于调控细胞生理功能、机械传感、液-液相分离、感觉受体功能以及蛋白质-核酸复合物形成等诸多生物过程至关重要。理解这些相互作用的本质具有基础研究与应用生物医学的双重意义。计算化学研究，特别是基于量子化学的计算，能够深入揭示这些复杂相互作用的细节。然而，为了对大型生物系统（尤其是其动力学）进行定量建模，需要足够简洁的模型。因此，在生物体系中，静电相互作用最常使用原子点电荷（partial (point) charges）来建模。历史上，已有多种方法用于确定这些电荷，早期方法涉及对量子力学计算得到的电子密度进行原子布居划分，但这些方法所得电荷严重依赖于划分方法和所用量子理论级别。后来，考虑共价键信息的划分方法得到了更“合理”的电荷。而应用最广泛的方法之一，则是通过拟合分子产生的静电势来获得点电荷。
尽管该领域普遍认识到静电势拟合电荷存在可转移性问题，但用于拟合的静电势信息本身可能并不完整这一点或许被忽视了。泊松方程揭示了静电势与其产生电荷分布之间的局部空间关系。如果VΦ是静电势网格点采样的体积，那么电荷密度ρ(r ∈ VΦ)可以从Φ(r ∈ VΦ)确定，但ρ(r ? VΦ)则成为一个外推，其质量取决于所使用的ρ数学模型。这尤其棘手，因为两个电荷分布之间的静电相互作用Eel取决于整个空间的精确电荷分布。在有限体积VΦ内拟合电荷模型，将导致区域外电荷分布的精度未知，进而以不可预测的方式影响静电相互作用能的计算。众所周知，静电势拟合电荷强烈依赖于用于拟合的特定体积VΦ>。本文认为，与静电势拟合相比，通过使用机器学习工具直接基于二聚体相互作用能数据训练模型，可以更直接、更准确地构建物理基础的力场，从而增强分子模拟在许多科学应用中的预测能力。
理论：从静电势拟合到二聚体能量训练
静电势ΦQM(r)由量子力学电荷分布ρQM(rj)产生。分子静电模型可以通过拟合电荷qj，使得分子力学计算的静电势ΦMM(r)（基于点电荷、高斯分布电荷或电多极子）来近似ΦQM(r)。然而，如公式所示，静电相互作用能的计算涉及对整个空间的积分。传统的静电势拟合方法在有限的采样体积内进行，无法保证得到的电荷分布能在所有距离上准确再现静电相互作用，尤其是在短程相互作用区域，而这对定量确定诸如配体结合能等过程至关重要。
本文的核心论点是，与其通过单体性质（如静电势）这一间接路径推导模型，不如直接基于对称性适应微扰理论计算的二聚体相互作用能（静电和诱导分量）来训练力场参数。这种方法利用Alexandria化学工具包这一机器学习软件，从零开始生成力场模型。本文评估了多种物理模型，其复杂度介于简单点电荷模型和基于多极子的复杂模型（如Amoeba、MASTIFF或CLIFF）之间。为了模拟电荷和/或交换各向异性，模型依赖于原子上的虚拟位点、水分子二等分线上的虚拟位点或孤对电子/σ空穴。
方法：数据库构建与模型训练
本文使用的数据库包含94个带电荷的氨基酸侧链类似物、水和无机离子的二聚体。首先使用GAFF力场对二聚体结构进行最小化，然后沿定义最短原子距离的向量进行扫描，并在接近范德华半径和的距离处生成随机取向的二聚体。使用SAPT2+(CCD)δMP2方法和aug-cc-pVTZ基组，通过Psi4软件计算了所有二聚体的静电和诱导能量分量。从生成的所有结构中，根据距离和能量标准筛选出14371个二聚体结构，其中9379个用于训练，4992个作为测试集。
此外，还以“传统”方式计算了侧链类似物、水和无机离子的静电势，即在距离原子1.4、1.6、1.8和2.0倍范德华半径处生成四层网格点，并使用MP2/aug-cc-pVTZ理论水平计算静电势。这些计算用于研究在模型中包含屏蔽电荷（通过将点电荷与高斯分布或Slater分布电荷结合，并拟合到静电势）的效果。为了比较机器学习模型与现有方案，使用Gaussian软件生成了Mulliken、Hirshfeld、ESP和CM5模型的原子电荷，并使用Antechamber生成了BCC电荷和RESP电荷。还从Alexandria库中获取了约5100个化合物的ESP电荷。
使用Alexandria化学工具包来训练分裂电荷平衡算法的参数（原子硬度η、电负性χ、键硬度Δη和键合原子间电负性差Δχ）。此外，当使用例如高斯分布电荷时，还训练分布宽度ζ；对于四点位水模型，还训练二等分线上虚拟位点的位置。在适当情况下，壳层和虚拟位点上的电荷大小以及原子极化率α也是训练的一部分。本文生成的模型特征汇总于表中。
对于可极化模型PC+GS4，包含了一个额外的吸引势以表示高阶诱导效应，其形式为指数衰减函数。该模型使用混合遗传/蒙特卡洛算法进行训练。
结果：模型评估与比较
无机离子的静电势分析：为了探究设计能准确再现量子化学静电势的电荷分布模型的可行性，本文以可解析处理的单原子离子为例。由于原子由带正电的原子核和围绕它的电子云组成，将代表原子核的点电荷核心与一个或两个高斯分布或Slater分布电荷（代表电子云）结合的模型似乎是合理的。对Li+、Na+、K+、F–、Cl–和Br–在0到4.5 ?距离上的静电势计算表明，单个点电荷不足以准确再现卤素离子的静电势（对于碱金属离子，差异较小但存在）。最复杂的模型（点电荷结合1S+2S Slater分布电荷）对量子化学静电势的拟合均方根偏差最低。然而，基于单体静电势拟合的模型，其预测的离子对静电相互作用能与SAPT参考值相比，在能量极小值距离处仍存在显著偏差，特别是当使用更宽距离范围（ESP0）拟合时，模型准确性甚至回落到与简单点电荷模型相近的水平。
静电势拟合电荷模型的准确性：使用Alexandria库中约5100个化合物的数据分析表明，基于点电荷的静电势与Hartree-Fock水平计算值的平均均方根偏差为5.7 kJ/mol e。将点电荷与高斯虚拟位点或Slater虚拟位点结合的模型，其静电势拟合误差略有降低。然而，当评估这些模型对二聚体相互作用能的预测时，即使结合了分布电荷，基于静电势拟合的模型对静电相互作用的再现仍然不佳，偏差在短距离处最大。
基于二聚体能量分量的机器学习模型：鉴于基于单体性质的预测准确性有限，本文使用机器学习来训练物理基础的静电和诱导模型。通过ACT训练分裂电荷平衡算法参数，以再现SAPT计算得到的静电能量（或静电加诱导能量）。结果显示，仅使用点电荷的模型，其静电相互作用均方根偏差约为17 kJ/mol，与ESP、CM5、BCC、RESP和MBIS等方法相当。使用高斯分布电荷的模型性能稍好（RMSD为14 kJ/mol）。当在模型中添加带有高斯或Slater分布电荷的虚拟位点后，静电能的均方根偏差比点电荷模型降低了约80%，表明非极化模型有巨大的改进空间。重要的是，基于单体性质的MBIS-S模型和基于静电势训练的双电荷模型的误差，是基于二聚体能量训练的相同计算复杂度模型的四倍。这证明了基于二聚体能量训练模型的有效性和优越性。
诱导与极化：为了模拟凝聚相，非极化模型需要隐式地纳入诱导效应。当在训练中同时考虑静电和诱导能量时，所有现有模型的静电加诱导能总和均方根偏差都超过35 kJ/mol，并且系统性地高估了该能量总和。将分裂电荷平衡参数训练为同时拟合静电和诱导能，会导致静电能的误差增大，因为此时电荷绝对值显著变大。为了构建有效的非极化模型，需要同时训练包括交换和色散在内的完整模型。最后，将PC+GV4模型中的虚拟位点设置为可极化，生成了PC+GS4模型。该模型在再现SAPT静电和诱导能总和方面表现出色，均方根偏差显著降低。
与现有模型的比较：将新模型与现有的TIP4P-Ew水模型结合离子点电荷、MBIS-S模型以及结合了CHARMM Drude可极化离子参数的SWM4-NDP水模型进行比较。对于离子-水静电相互作用，新模型（PC+GV4和PC+GS4）的预测与SAPT参考值更为接近，均方根偏差显著低于已发表模型。其中，非极化的TIP4P-Ew/点电荷模型预测的阳离子-水相互作用过强，阴离子-水相互作用过弱；而SWM4-NDP/CHARMM Drude模型预测的相互作用则普遍偏弱。对于离子-离子相互作用，新模型也优于点电荷模型、MBIS-S模型以及经验性的可极化模型。在诱导能方面，PC+GS4模型的预测也比CHARMM Drude模型更接近SAPT值，尽管仍有改进空间。
氨基酸侧链类似物的相互作用：对于蛋白质侧链类似物与离子或水的相互作用，新模型（PC+GV4和PC+GS4）在预测静电能方面比广泛使用的BCC、RESP和MBIS-S模型准确得多，具有更低的均方根偏差和接近零的平均符号误差。
讨论与展望
设计分子间势能的目标是准确预测相互作用，例如在药物设计或材料设计中。这需要能够忠实再现参考数据（如来自SAPT、量子化学计算数据库乃至实验数据）的物理模型。本文的研究表明，广泛使用的通过拟合静电势来确定原子电荷的方法受到信息缺失的制约，导致其预测的相互作用能在定量甚至定性上不正确。事实上，大多数基于单体化合物确定点电荷的方法在预测静电相互作用能方面都不够精确。
研究结果重申了从点电荷模型过渡到高斯分布电荷模型，乃至更好的PC+GV4模型，能够以较低至中等的计算成本显著提高非极化经验力场的准确性。然而，基于点电荷和简单范德华势的力场之所以能合理再现物质的凝聚相性质，很大程度上依赖于不同能量项之间的误差补偿，这限制了这类力场的可转移性。
得益于SAPT等技术能够提供能量分量，多个研究小组已着手开发基于SAPT分量训练的力场。例如，MASTIFF模型部分基于SAPT能量训练，Amoeba力场也据此进行了评估。此外，CLIFF力场结合了分子间相互作用能的物理方程与机器学习模型以实现自动参数化，并利用基于SAPT的BioFragment Database进行验证。该模型采用单体静电描述，但部分阻尼参数基于SAPT二聚体数据训练。与这些工作相比，本文的结果表明，可以纯粹基于带电氨基酸类似物、无机离子和水的二聚体能量，推导出精确的静电相互作用模型。诱导参数也可以用同样的方式训练，尽管相关模型仍需改进。
静电势拟合方法因在拟合时通常忽略短距离静电势而受限于信息不足。因此，对于本文使用的数据集，静电势拟合电荷不能很好地再现静电相互作用。尽管存在更详细的基于单体的模型，但基于单体推导静电模型的方法总体上是间接的。本文证明，直接根据SAPT计算的二聚体静电和诱导能来推导模型是一种可行的替代方案，能够产生化学精度（RMSD < 1 kcal/mol）的静电能。通过引入显式极化，也可以训练模型以再现SAPT的诱导能。虽然二聚体训练也受到可用数据和所用物理模型的限制，但SAPT直接提供了关于相互作用能的信息，而这正是力场计算的主要目的。
为了将力场中的静电和诱导分量提升到超越本文描述的水平，可能需要纳入更精细的物理描述，包括电荷转移和涨落，以及精确的诱导模型。开发一个完整的力场超出了本文范围，但可以指出，在开发色散和交换模型时，需要考虑各向异性，或许还需要考虑多体效应。ACT中系统化力场开发的工具，可以通过基于不同物理模型从零开始生成新力场，来助力这一事业。总之，结合力场科学进行模型间的直接比较，以及利用机器学习设计物理基础的、与SAPT计算一致的力场，这些工具将共同推动力场开发的快速进展，并提升分子模拟在众多科学领域应用中的预测能力。
