基于神经网络势的、与体系结构无关的绝对溶剂化自由能计算

《The Journal of Physical Chemistry Letters》：Architecture-Independent Absolute Solvation Free Energy Calculations with Neural Network Potentials

本文提出一种基于邻居列表的4D脱耦方法，用于通过神经网络势（NNP）进行自由能模拟（FES）中的原子或分子脱耦合。该方法通过调整原子间距离实现溶剂与溶质的相互作用剥离，适用于任意NNP架构，避免了传统力场或QM/MM方法的局限性。验证实验显示该方法在循环闭合测试中误差极低（0.0±0.4 kcal/mol），且MACE-OFF23(S/M)模型计算出的溶解自由能（ASFE）与实验值吻合较好（平均误差0.4 kcal/mol），但模型复杂度增加可能导致误差。研究为NNP在FES中的应用提供了通用方法论。

　　在现代计算化学中，自由能模拟（Free Energy Simulations, FES）已成为预测分子系统热力学性质的重要工具。自由能差（ΔG）可以提供关于反应或过程自发方向的关键信息，因此，其计算在药物设计、材料科学以及生物分子相互作用研究等领域中具有广泛的应用价值。例如，通过计算溶剂化自由能差（Solvation Free Energy Differences, SFED）和分配系数，研究人员能够理解分子在不同环境中的稳定性和溶解性。在这些应用中，绑定自由能的计算尤为重要，因为它们被广泛用于药物分子的优化和筛选。

然而，传统分子力学（MM）力场在计算自由能差时存在一定的局限性。尽管MM力场在大多数FES中被使用，但由于其近似性质，所得到的ΔΔG值往往不够精确。量子化学（QM）方法理论上可以提供更高的精度，但其计算成本较高，难以在大规模或长时间的模拟中使用。近年来，神经网络势能（Neural Network Potentials, NNP）的快速发展为解决这一问题提供了新的思路。NNP通过训练大量高精度的量子力学数据，能够以较低的计算成本获得接近量子化学精度的能量和力场信息。这使得NNP在自由能模拟中的应用成为可能，尤其是在涉及多个纳米秒的模拟系统时，如溶质-溶剂或小蛋白质-配体体系。

尽管NNP在计算精度上具有优势，但在进行非物理的化学变换（如原子或分子的“消失”）时仍面临挑战。这些变换是FES的核心步骤，涉及在非物理状态下逐步改变分子间相互作用。与传统力场不同，NNP对原子间的距离更加敏感，这意味着当原子具有不合理的键数（如三价或五价碳）时，NNP无法像力场那样容忍这些错误。因此，如何在不引入非物理状态的前提下，实现这些变换成为研究的重点。

在这一背景下，研究者提出了一种基于邻居列表（Neighbor List）修改的方法，以实现对原子或分子的直接解耦（Decoupling）。该方法的核心思想是通过调整邻居列表中涉及特定原子的距离，使得这些原子与系统其余部分的相互作用被逐步关闭。邻居列表是所有NNP计算的基础，它记录了在一定截断距离（r_cut）范围内所有原子对的相互作用。因此，这种方法具有普适性，适用于各种NNP架构，而无需针对特定模型进行调整。

为了验证该方法的正确性，研究者采用了“循环闭合”（Cycle Closure）作为标准。循环闭合是一种用于检验相对自由能计算是否一致的方法，通过比较不同路径下计算出的自由能差是否一致，来判断模拟方法是否可靠。例如，在计算两种同分异构体（Tautomer）的相对溶剂化自由能差（Relative Solvation Free Energy, RSFE）时，可以通过两种不同的路径进行计算：一种是通过计算两种同分异构体的绝对溶剂化自由能差（Absolute Solvation Free Energy, ASFE）的差值；另一种是通过直接的化学变换计算相对自由能差。如果两种方法得到的结果一致，则说明该方法是正确的。

在实验中，研究者测试了四种不同的距离调整方案，分别对应不同的λ（耦合参数）函数。这些方案分为线性调整（Linear Shifting）和4D调整（4D Shifting）两类。线性调整方案通过在距离上施加线性偏移，确保在λ = 1时，所有涉及特定原子的距离超过截断距离，从而将其排除在邻居列表之外。4D调整方案则引入了一个4维向量，其前三个分量是原始的距离向量，第四个分量是耦合参数λ与截断距离的乘积。这种方法更接近于传统力场中使用的软核势（Soft-Core Potential），因为它通过改变距离的方式逐步关闭相互作用。

在模拟细节方面，研究者使用了OpenMM软件包，因为它支持NNP和混合NNP/MM系统。为了确保计算的稳定性，他们采用了单精度（float32）的数值计算，并选择了特定的λ值进行模拟。对于小分子的溶剂化自由能计算，他们使用了15个λ值，从0到1，以确保逐步解耦的准确性。在每种λ状态下，进行了1纳秒的兰格-辛格动力学（Langevin Dynamics）模拟，并结合蒙特卡洛压力统计（Monte Carlo Barostat）来维持恒定的压力和温度。模拟过程中，每0.25皮秒保存一次坐标数据，前20%的数据被排除，以确保系统达到平衡状态。最终，仅使用每四个保存的坐标数据进行自由能计算，以减少计算误差。

在自由能差的计算中，研究者最初尝试使用MBAR（Multistate Bennett Acceptance Ratio）方法，但由于数据在λ = 1时出现异常，导致MBAR收敛困难。因此，他们改用BAR（Bennett Acceptance Ratio）估计器，计算相邻λ状态之间的自由能差，并将这些结果相加。这种方法避免了因距离调整方式不同而产生的误差，使得自由能差的计算更加稳健。

在验证方法的正确性方面，研究者对两种同分异构体对进行了循环闭合测试。结果显示，两种同分异构体对的循环闭合误差分别为0.0和-0.2 kcal/mol，考虑到误差范围在0.2-0.4 kcal/mol之间，这些误差可以视为可接受的。此外，他们还比较了不同距离调整方案下的结果，发现四种方案在计算绝对溶剂化自由能时结果一致，进一步证明了该方法的有效性。

在实际应用中，研究者还计算了六种小分子的溶剂化自由能，并与实验数据进行了比较。结果显示，MACE-OFF23(S)模型在大多数情况下与实验数据吻合较好，而MACE-OFF23(M)模型则在部分分子上出现了较大的偏差。这表明，尽管NNP在计算精度上具有优势，但其在描述复杂系统（如溶剂化过程）时仍需进一步优化。此外，研究者还指出，当前大多数NNP模型主要基于单分子数据训练，因此在预测多体相互作用时表现有限。

总的来说，基于邻居列表的解耦方法为使用NNP进行自由能模拟提供了一种新的思路。该方法不仅适用于现有的NNP架构，也为未来的发展奠定了基础。尽管当前的模拟结果与实验数据之间仍存在一定差距，但随着NNP训练数据的扩展和模型的优化，这种方法有望在更广泛的化学体系中得到应用。同时，研究者也指出，该方法在某些情况下需要结合额外的约束（如平底谐波势）以确保模拟的稳定性，这表明在实际应用中还需进一步探索和优化。未来的工作将聚焦于如何更好地结合不同类型的NNP模型，以提高其在复杂系统中的预测能力，并进一步验证其在自由能计算中的鲁棒性和稳定性。