机器学习势能能否可靠计算水中酸碱解离常数？

在单个溶质的盒子里计算pKa

这篇论文要回答的是：酸碱解离常数能否自下而上地从原子尺度模拟得到，而不是靠对已知化合物做数据拟合。目标是水溶液中某个物种的绝对pKa或pKb，它在稀溶液极限下用浓度和标准浓度c° = 1 mol/L定义，pKa = −log10 Ka。对水中的共轭酸碱对，这些常数与溶剂自电离常数满足pKw = pKa + pKb，作者后来把这一关系当作不依赖实验的一致性检验。

难点在于该反应会断裂和形成共价键，因此驱动动力学的势能必须是反应性的。密度泛函理论在所需时间尺度上过于昂贵，而解离自由能可达100 kJ/mol；解离又是稀有事件，需要增强采样；质子的量子行为也不能忽略。作者的方案是用机器学习原子间势能（MLIP）替代电子结构计算这一瓶颈——这是一种经过训练、能以很小一部分成本复现参考量子化学方法的模型——再用有针对性的近似处理其余成本，而不是偷工减料。

同一自由能曲线可给出不同pK值

论文的核心方法学主张，是如何从模拟中读出平衡常数。作者使用可逆功定理：把两个解离产物之间的平均力势W(r)对缔合盆地的玻尔兹曼权重积分到截断距离rc，就得到解离常数。由于这个积分是一个概率，它衡量离子对相对于非相互作用参考态处于缔合盆地的可能性，因此任何能识别同一盆地的解离距离都会给出同一个常数。

文献中常用的另外两种估计方式不具备这一性质。取缔合态与解离态自由能极小值之差，结果依赖集合变量的雅可比，也不对应任何明确定义的标准态；对水的自电离，作者报告该估计方式与在同一曲线上用可逆功定理得到的值相差4.2个pK单位。把玻尔兹曼权重分别在缔合区和解离区积分，思路更接近，但强烈依赖所选解离区间，由此产生的盒尺寸依赖也难以校正。可逆功路线则把所有任意性集中到一个加性常数C上。

一个常数承担有限尺寸修正

加性常数C确定对分布函数的归一化，而且很关键：C上kBT ln10 ≈ 5.7 kJ/mol的误差就对应整整一个pK单位。作者确定C的办法是，在接触离子对之外的间距上，把模拟得到的平均力势匹配到两个点电荷在与模拟盒同尺寸、同形状的周期盒中的屏蔽库仑相互作用，后者用Ewald型求和计算。由于参考的是周期盒内的相互作用而非无限稀释下的相互作用，周期静电和盒几何的影响被一并吸收进C，不额外增加成本。

有两个细节让这一做法可行。间距为r时，可及立体角在超过L/2后被盒面截断，因此参考曲线必须包含与雅可比相关的修正项；作者指出，对他们的十二面体盒子，这一项通过把点电荷放在网格上数值处理。当产物带净电荷时，中和背景电荷会移动Ewald结果，作者减去相应的Wigner自能项。他们还报告，即使MLIP不含显式长程静电，只要模拟曲线在匹配窗口内符合屏蔽库仑形状，校准仍然有效。把rc从4 Å变到8 Å的敏感性测试中，没有任何物种的变化超过一个pK单位。

带符号距离把酸性与碱性分支放在同一轴

解离产物并不绑定在固定原子上：多余质子和氢氧根中的质子空穴都通过类Grotthuss机制在氢键网络中迁移。因此作者用平滑权重把每个氢分配给附近的氧位点，得到每个位点的形式电荷：水合氢离子接近+1，氢氧根接近−1，其余接近零，并在质子转移过程中平滑过渡。这种平滑性是必需的，因为该集合变量要用于加偏置。

随后，分离度被定义为从酸碱中心参考点出发、按电荷加权的最小镜像距离之和。酸性解离给出正的带符号距离，碱性解离给出负值，于是像H2PO4−这样的两性离子的两个反应出现在同一条一维曲线的原点两侧，未解离物种位于零附近。对水的自电离，两个产物都是流动性的，同一思路写成类偶极表达式，参考点用周期性的圆均值确定。实际好处是只需收敛一维自由能曲线，而不必像此前的集合变量那样收敛二维曲面。

元动力学、副本交换与质量积分修正

采样分两步：元动力学先构建偏置，推动体系沿带符号距离移动；随后在固定偏置下做副本交换，得到收敛的玻尔兹曼权重。每个物种从不同初始构型出发做十次独立运行，直接给出统计不确定度估计。对水，这些运行之间的离散度对应解离区几kJ/mol的差异，以及经典pKw上±0.37个pK单位的不确定度。作者还指出，从缔合态到解离态密度增加超过1%，因此采样在恒压下进行。

核量子效应单独处理，而不是沿整条路径跑完整的路径积分模拟。作者通过对核质量做热力学积分，计算自由能差的量子修正，这只需收敛两个端点态的量子动能。这一步对主要数字影响最大：水的修正为−3.54 ± 0.07个pK单位，磷酸各步解离的修正从−0.11到−3.05个pK单位不等。由于同一种酸的各级解离之间修正差异如此之大，作者认为不能用常数平移替代，即使在同一个势能下也不行。

量子修正让磷酸各级解离进入化学精度

论文报告了300 K下水以及六个磷酸解离步骤的经典值和量子修正值。对水，经典pKw为17.86 ± 0.37，经−3.54修正后为14.31 ± 0.37，实验值为13.93。对磷酸，第一级pKa从3.00 ± 0.06变为2.10 ± 0.09，实验值2.16；第二级从9.00 ± 0.22变为6.87 ± 0.23，实验值7.19；第三级从15.06 ± 0.28变为12.02 ± 0.29，实验值12.33。整组数据相对实验的平均绝对误差从2.19降到0.45个pK单位，作者认为这落在传统1 kcal/mol化学精度阈值所对应的约0.7个pK单位之内。

最明显的离群值是H2PO4−的碱性解离：经典值15.22 ± 0.13，经−2.38修正后为12.85 ± 0.14，实验值为11.77。作者指出，这是整组中唯一一个产物（中性H3PO4和OH−）之间不存在电荷-电荷相互作用的反应，因此无法对它应用C的库仑校准。他们表示，这一偏差可能来自参考方法的残余误差、模型的误差，或小盒子中反离子的存在，不做进一步的参考计算就无法区分。

共轭对求和提供不依赖实验的检验

由于六个磷酸常数并不独立，把每种酸的pKa与其共轭碱的pKb相加，可得到三个独立的pKw估计。量子值为：H3PO4/H2PO4−对14.95，H2PO4−/HPO4²−对14.05，HPO4²−/PO4³−对14.20，与直接计算的14.31比较。后两者在合并不确定度内一致；第一个偏离0.6个pK单位。经典求和为18.22、17.60和17.35，对应直接值17.86，偏差相当且没有明显离群者，因此第一对的偏差来自经典曲线约0.4个单位的偏差加上量子修正约0.3个单位的差异。

作者建议把这一检验作为任何多元酸或共轭对的常规诊断，理由是：与实验吻合可能来自误差相互抵消，而违反pKw = pKa + pKb则不可能。论文的实际限制在势能上：基础模型让新溶质的反应性MLIP只需几千次参考计算即可获得，但微调所用的理论水平决定最终精度，作者认为GGA泛函不太可能给出定量结果。他们使用revPBE0-D3，该泛函在水体系上有良好记录，同时指出即便这一选择也可能受益于误差抵消。可迁移的部分是流程本身：从水转到六个磷酸解离步骤，序参量定义和分析都无需改动，只需按解离自由能的大小调整元动力学参数。

本文引用的文献