如何从 Metadynamics 自由能面准确计算平衡常数pKi？

ZQing

各位计算化学的老师、同仁好，
我最近在研究如何通过 Metadynamics (增强采样) 得到的自由能面 (FES) 来计算金属络合物的稳定常数 (pKi)。主要参考了近期 JCTC 上的一篇文章 (Sagresti et al., 2025, DOI: 10.1021/acs.jctc.5c01079)。
文章计算稳定常数 (pKi)的核心逻辑是通过玻尔兹曼分布结合质量守恒来求解：
首先对于反应
ML_{i-1}+L=ML_i
平衡常数有如下定义

两个物种的平衡浓度比写作

其中，n_i 是与 [ML_i] 配位物种相关的equilibrium population；ΔG_{ij} 是第 i 个和第 j 个配位状态之间的自由能差
所以K_i又可以写作

其中未知的 free ligand concentration [L] 可以根据质量守恒条件确定：

其中 [L0] 是初始配体浓度，Ns 是可能的配位态总数。

我的复现过程与疑惑：
我尝试复现文中 Cd(II)-en 体系的数据（10个金属，30个配体，盒子 330 nm³，即初始 [M_0]=0.05M, [L_0]=0.15M）。也就是文章中的图2c所对应的第一个过程
Cd(II)+en=Cd(II)en

对应的ΔG则是-17.2kJ/mol 也正如如Table S5所示

将这个值带入到上面的平衡浓度比值的式子里，也就是
[Cd(II)en]/[Cd(II)]=exp(17.2/k_B·T)= 988.3
接下来要用这个数除以游离配体浓度[L]得到K_i，
因为盒子里有10个金属，30个配体，我的理解就是[L]相当于是游离的20个配体的浓度也就是
20/(6.022×10^23)/330×10^(-24)=0.05mol/L
给988.3/0.05就得到K_i=19766
算出来logk1也就是pK1是4.3
但是文章中算出来的(Figure 2e)是5.1，差了蛮多

所以我想请教做过类似自由能算平衡常数的大佬指点迷津，看看我做的哪里有问题。非常感谢！
这篇文献以及Supplementary Information附上

ZQing

KazusaT 发表于 2026-4-29 21:42
我没有具体算这个数，注意17.2的单位是kJ/mol，你应该除的是R·T且应该是17200/(RT)

是的，我是除的RT，这里手误写错了，算得就是exp(-ΔG/RT) = 988.3，ΔG是转化成 J/mol算的。
好的，我再用这种方法算一下试试。
作者确实是做了很多比例，但是并没有求平均。

谢谢老师解答！

KazusaT

我没有看文章，也没有具体尝试算，但指出几个要注意的点，你可以看看

[Cd(II)en]/[Cd(II)]=exp(17.2/k_B·T)= 988.3

注意17.2的单位是kJ/mol，你应该除的是R·T且应该是17200/(RT)

因为盒子里有10个金属，30个配体，我的理解就是[L]相当于是游离的20个配体的浓度也就是

并非，注意到所有的[MLn]均可以用[M]表示，以Cd为例，求解[Cd]+[Cden]+[Cden2]+[Cden3]=0.05可得各物种浓度，进而得到[L]

算出来logk1也就是pK1是4.3

小bug，logK1是lgK1，pK1其实是-lgK1，不过这里文章都用反了，对你计算暂时没什么影响。而且他这个表格和图片的ΔG符号是相反的也让我挺迷惑的。

另外注意到你贴出来的文章的lgK结果图是有error bar的，并且表里面作者尝试了不同的比例去计算自由能，文章是否是对几种情况求平均得到的？