分子动力学模拟与结合自由能计算

可视化演示

展示分子动力学模拟与结合自由能计算三大核心：MD模拟流程（系统构建→最小化→平衡→产出100ns→轨迹分析五阶段）、力场与采样（AMBER/CHARMM/OPLS三大力场+MD/REMD/MetaD/Umbrella四种增强采样）、结合自由能（MM/PBSA/FEP/TI/ABF四种方法+范德华/静电/极性/非极性能量分解总结合能-25kJ/mol）,三模式切换理解MD全貌。

舞池边的摄像机——MD 基本概念

编舞师的第一台摄像机要记录什么？不是一张静态合影，而是一段连续流动的舞蹈。分子动力学模拟做的事情完全一样：它不再把蛋白和配体看成一张冷冻照片（那是分子对接和晶体学的工作），而是把每个原子当成一个有质量、有速度的小舞者，用牛顿力学驱动它们在时间中前进，一帧一帧“拍”出整段舞。

录像的每一帧，也就是系统里所有原子在某一瞬间的坐标集合，叫作一个构象；连接所有帧的那条路径，叫作轨迹。模拟开始时会给出初始坐标和速度，然后循环做两件事：根据当前位置计算每个原子受力，再更新速度与位置。这就是“拍摄”的基本循环。由于分子世界原子数量庞大，MD 采用经典力学近似来模拟，配合控温手段（热浴）让体系维持生理温度，避免“越跳越热”。

$$\displaystyle m_i \frac{d^2 \mathbf{r}_i}{d t^2} = \mathbf{F}_i = -\frac{\partial U}{\partial \mathbf{r}_i}$$

1. 质量乘加速度：$m_i$ 是第 $i$ 个原子的质量，$\frac{d^2 \mathbf{r}_i}{d t^2}$ 是它位置对时间的二阶导，也就是加速度。质量越大越难被推动，像体重大的舞者起步更慢。

$$\displaystyle m_i \frac{d^2 \mathbf{r}_i}{d t^2}$$

2. 合外力：$\mathbf{F}_i$ 表示其他所有原子和环境施加在第 $i$ 个原子上的合力，相当于舞伴之间互相拉扯的全部力量加起来。

$$\displaystyle \mathbf{F}_i$$

3. 势能负梯度：$-\frac{\partial U}{\partial \mathbf{r}_i}$ 说明力总是指向让体系势能降低的方向，分子会顺着“更舒服”的姿态变化，正如舞者自然滑向更顺的节拍。

$$\displaystyle -\frac{\partial U}{\partial \mathbf{r}_i}$$

你按下录制键，舞池里的分子世界开始流动：第一帧是牵手，第一百帧是旋转。一段真实可信的录像，正是后面所有统计与分析的原片——但前提是，摄像机得遵循正确的“舞蹈规则”。

舞蹈规则手册——力场设置

摄像机再高级，也拍不出没有规则的舞蹈。力场就是分子世界的舞蹈规则手册：它规定化学键像弹簧一样伸缩时用多大力，键角弯曲多少会被拉回来，两个原子靠太近会产生多大排斥，带电原子之间又有怎样的静电牵拉。

操作层面，力场是一套带参数的数学表达式，主流的有 AMBER、CHARMM、GROMACS 内置力场等。选择力场就是选择一本舞厅手册：跳华尔兹和跳探戈，脚底摩擦与身体倾斜的规则大有不同。此外还要把体系泡进“水盒子”里——选择合适的水模型（如 TIP3P）和周期性边界盒子，并且完成加氢、离子中和等系统初始化步骤，舞池才算布置完毕。

$$\displaystyle U = \sum_{\text{bonds}} k_b (r - r_0)^2 + \sum_{\text{angles}} k_\theta (\theta - \theta_0)^2 + \sum_{\text{torsions}} V_n \bigl[1 + \cos(n\phi - \gamma)\bigr] + \sum_{i<j} 4\varepsilon_{ij}\!\left(\frac{\sigma_{ij}^{12}}{r_{ij}^{12}} - \frac{\sigma_{ij}^{6}}{r_{ij}^{6}}\right) + \sum_{i<j} \frac{q_i q_j}{4\pi\varepsilon_0 r_{ij}}$$

1. 键伸缩弹簧：$\sum_{\text{bonds}} k_b (r - r_0)^2$ 把每个化学键看成弹簧，键长偏离平衡值 $r_0$ 越远惩罚越大，像双人舞中两人之间绷着一条弹力带。

$$\displaystyle \sum_{\text{bonds}} k_b (r - r_0)^2$$

2. 键角弯曲弹簧：$\sum_{\text{angles}} k_\theta (\theta - \theta_0)^2$ 约束三个原子形成的角度，防止骨架弯出“不可能”的形状，像舞者的膝盖弯曲角有正常范围。

$$\displaystyle \sum_{\text{angles}} k_\theta (\theta - \theta_0)^2$$

3. 二面角扭转：$\sum_{\text{torsions}} V_n \bigl[1 + \cos(n\phi - \gamma)\bigr]$ 描述绕单键旋转的势垒，决定骨架能“转得多自由”，像舞者腰部旋转的灵活度。

$$\displaystyle \sum_{\text{torsions}} V_n \bigl[1 + \cos(n\phi - \gamma)\bigr]$$

4. 范德华作用：$\sum_{i<j} 4\varepsilon_{ij}\!\left(\frac{\sigma_{ij}^{12}}{r_{ij}^{12}} - \frac{\sigma_{ij}^{6}}{r_{ij}^{6}}\right)$ 描述原子对间的吸引与近距离强排斥，像两人太近会撞到、太远又感觉不到对方。

$$\displaystyle 4\varepsilon_{ij}\!\left(\frac{\sigma_{ij}^{12}}{r_{ij}^{12}} - \frac{\sigma_{ij}^{6}}{r_{ij}^{6}}\right)$$

5. 静电作用：$\sum_{i<j} \frac{q_i q_j}{4\pi\varepsilon_0 r_{ij}}$ 是电荷间的库仑力，异号相吸同号相斥，像舞池里带着正负气场的舞者互相靠近或弹开。

$$\displaystyle \frac{q_i q_j}{4\pi\varepsilon_0 r_{ij}}$$

手册定好了，舞者才知道什么动作“合法”。但规则正确只是第一步，接下来要回答一个更实际的问题：这段舞，到底该拍多久？

该拍多久——采样方法

一段三分钟的华尔兹不能代表舞者的全部实力。判断一对舞伴是否默契，得看他们跳遍慢节奏、快节奏、探戈、伦巴之后的表现。MD 模拟面临同样的困境：我们只能从有限长度的轨迹里，去推断整个“舞步空间”（构象空间）的平均表现——这一步叫采样。

采样不足的典型症状是“困在局部”：配体被卡在一个错误的姿势里反复振动，始终没能访问重要的构象，就像舞者一直在舞池角落重复同一个舞步。时间步长的选择也很关键：普通原子运动最快的是氢原子振动，步长通常设为 1–2 fs，步长太大会导致能量爆炸。为了跳出不舒适的角落，人们发展出增强采样技术（副本交换、元动力学、伞形采样），相当于主动给舞者发“必跳清单”。最后还要用 RMSD 收敛曲线等判据确认：录像确实覆盖了该跳的舞步。

$$\displaystyle \langle A \rangle = \lim_{T \to \infty} \frac{1}{T} \int_0^T A\bigl(\mathbf{r}(t)\bigr)\, dt$$

1. 可观测量：$A\bigl(\mathbf{r}(t)\bigr)$ 表示在 $t$ 时刻的构象下我们能读出的一个量（比如配体到口袋的距离），相当于从某一帧画面里量出两人的牵手距离。

$$\displaystyle A\bigl(\mathbf{r}(t)\bigr)$$

2. 时间平均：$\frac{1}{T}\int_0^T \cdots\, dt$ 把整段录像里的读数累加起来再除以时长，等于把每帧的“亲密度”汇总成整支舞的平均水平。

$$\displaystyle \frac{1}{T}\int_0^T A\bigl(\mathbf{r}(t)\bigr)\, dt$$

3. 无穷长极限：$\langle A \rangle = \lim_{T \to \infty} \cdots$ 断言：当录像无限长且舞步覆盖完整时，时间平均等于整个构象空间的理论均值——这就是各态历经假设。

$$\displaystyle \langle A \rangle = \lim_{T \to \infty} \frac{1}{T} \int_0^T A\bigl(\mathbf{r}(t)\bigr)\, dt$$

当录像足够长、舞步足够全，时间平均才等于真实水平——这是后面所有“打分”的前提。那么，分数本身到底是什么？

关系牢固度怎么打分——自由能基础

录像拍好了，该回答最初的问题：这对舞伴的关系有多牢固？舞会上，有人一见钟情，有人跳两步就散——分子世界用“结合自由能”给这种牢固程度打分。$\Delta G_{\text{bind}}$ 越负，说明从“各自独舞”到“牵手共舞”越容易发生，配体越愿意待在蛋白身边。

严格定义上，结合自由能是复合物状态与分离状态（蛋白单独 + 配体单独）的自由能之差。它还可以拆成两笔账：焓变（能量上的“投缘”，包括氢键、范德华、静电的直接贡献）与熵变（结合后自由度损失的代价，如舞步被约束）。实验上，它与解离常数 $K_d$ 直接关联：$\Delta G_{\text{bind}} = RT \ln K_d$，因此自由能越负，解离常数越小，活性越强——这就是“亲和力度量”的由来。

$$\displaystyle \Delta G_{\text{bind}} = G_{\text{complex}} - G_{\text{protein}} - G_{\text{ligand}} = \Delta H - T\Delta S$$

1. 三个状态的自有分数：$G_{\text{complex}}$、$G_{\text{protein}}$、$G_{\text{ligand}}$ 分别是复合物、单独蛋白、单独配体的自由能，相当于独舞和共舞各自的“能量水平”。

$$\displaystyle G_{\text{complex}} - G_{\text{protein}} - G_{\text{ligand}}$$

2. 结合自由能差值：$\Delta G_{\text{bind}}$ 就是“牵手后的状态”减去“各自独舞的状态”，差值越负，说明牵手状态越占优、越难被拆散。

$$\displaystyle \Delta G_{\text{bind}} = G_{\text{complex}} - G_{\text{protein}} - G_{\text{ligand}}$$

3. 焓熵拆账：$\Delta H - T\Delta S$ 把差值拆成焓项（直接相互作用的净变化，像两人磁场投不投缘）与熵项（结合后自由度损失，像舞步被约束得不自由）。

$$\displaystyle \Delta G_{\text{bind}} = \Delta H - T\Delta S$$

明白了“分数”的含义，剩下的问题是怎么从录像里把它算出来。这正是后面两节的分岔路口：统计平均，还是换人实验。

逐帧回放算平均——MM/PBSA 计算

MM/PBSA 的思路像评委逐帧回放录像：每一帧都算一次“瞬间打分”，然后全部加起来取平均。瞬时打分包含三部分：配体与蛋白直接接触的分子力学能（键合、范德华、静电），以及“泡在水里”的溶剂化成本——极性部分用泊松-玻尔兹曼（PB）方程估算，非极性部分用溶剂可及表面积（SA）模型估算，最后再扣掉熵的损失。

实际操作中，我们从 MD 轨迹里等间隔抽帧（比如每 100 ps 一帧），对每帧重新计算能量与溶剂化修正。因为熵项的计算极其昂贵且粗糙，实践中常省略或用正则模（NM）近似。正因如此，MM/PBSA 的价值在于“相对排名”：用同一套打分标准给一批化合物排队，谁的分数靠前，谁更值得进入下一轮实验。

$$\displaystyle \Delta G_{\text{bind}} \approx \bigl\langle E_{\text{MM}} \bigr\rangle + \bigl\langle G_{\text{PB}} \bigr\rangle + \bigl\langle G_{\text{SA}} \bigr\rangle - T\bigl\langle S \bigr\rangle$$

1. 分子力学能量平均：$\bigl\langle E_{\text{MM}} \bigr\rangle$ 是轨迹各帧的气相力学能平均值，包含键合、范德华与静电三部分，相当于逐帧统计两人直接接触的“瞬时互动强度”。

$$\displaystyle \bigl\langle E_{\text{MM}} \bigr\rangle = \bigl\langle E_{\text{bonded}} + E_{\text{vdW}} + E_{\text{elec}} \bigr\rangle$$

2. 极性溶剂化项：$\bigl\langle G_{\text{PB}} \bigr\rangle$ 用泊松-玻尔兹曼方程估算水等极性溶剂的“溶解成本”，相当于舞池里的空气阻力对每个舞步收了多少“场地费”。

$$\displaystyle \bigl\langle G_{\text{PB}} \bigr\rangle$$

3. 非极性溶剂化项：$\bigl\langle G_{\text{SA}} \bigr\rangle$ 用溶剂可及表面积模型估算疏水相关的非极性成本，相当于“沾水就湿鞋”的隐性消耗。

$$\displaystyle \bigl\langle G_{\text{SA}} \bigr\rangle$$

4. 熵项扣减：$-T\bigl\langle S \bigr\rangle$ 统计结合后舞步自由度的损失，通常用正则模或简谐近似估算，是整个公式中最粗糙、最昂贵的一项。

$$\displaystyle -T\bigl\langle S \bigr\rangle$$

“平均打分法”擅长给一批化合物排队，但绝对数值仍与实验有差距。如果你想知道“某个微小化学改造到底能不能增强结合”，就需要更精密的工具——换人实验。

精准换人实验——FEP 自由能微扰

与其把整晚录像反复统计，不如直接做实验：把男舞伴的鞋底从硬底换成软底，其他什么都不变，观察两人的关系分数变化多少。自由能微扰（FEP）就是这个思路——它只计算“改变一点点”（比如把配体的一个氢换成氟）带来的相对自由能变化 $\Delta\Delta G$，而不是费力去算绝对结合自由能。

为了避免一步换人太猛导致画面崩溃，FEP 引入 λ 参数把 A 状态平滑插值到 B 状态，切成几十个中间窗，每个窗都做独立的 MD 平衡与采样。热力学循环保证：只要分别算得出“在水里变的成本”和“在口袋里变的成本”，两者之差就是结合的相对变化。每个窗采样充分、前后状态重叠良好，结果是目前计算化学中最接近实验精度的预测之一。

$$\displaystyle \Delta G_{A \to B} = -k_B T \ln \left\langle \exp\!\left( -\frac{U_B - U_A}{k_B T} \right) \right\rangle_A$$

1. 两个状态的能量差：$U_B - U_A$ 表示在同一个构象上分别用 A、B 两套势能函数算出的能量差，相当于“同一段舞蹈在不同鞋底下的费力程度差”。

$$\displaystyle U_B - U_A$$

2. 玻尔兹曼权重：$\exp\!\left( -\frac{U_B - U_A}{k_B T} \right)$ 把能量差折算成统计权重：差值越大，这一帧的贡献越小，像“鞋底差异太大时这段舞步几乎不能代表新状态”。

$$\displaystyle \exp\!\left( -\frac{U_B - U_A}{k_B T} \right)$$

3. A 态平均：$\left\langle \cdots \right\rangle_A$ 表示在 A 状态的大量采样构象上求平均，强调“换人实验”的基准是原状态拍到的全部画面。

$$\displaystyle \left\langle \exp\!\left( -\frac{U_B - U_A}{k_B T} \right) \right\rangle_A$$

4. 换算自由能：$-k_B T \ln(\cdots)$ 把平均权重倒过来换算成自由能差，这就是 Zwanzig 方程的全部秘密。

$$\displaystyle \Delta G_{A \to B} = -k_B T \ln \left\langle \exp\!\left( -\frac{U_B - U_A}{k_B T} \right) \right\rangle_A$$

换人实验精度高，成本也高——它要求每个中间态的采样都足够充分，相当于录像里的每一帧都经得起慢放。而无论是统计平均还是换人实验，背后都有一笔逃不掉的计算账单。

机位、内存与卡顿排查——性能与调试

高清录像一场舞会并不便宜：机位越多（并行计算），每秒能拍的帧数越多；但每当需要把各机位的画面合成一条时间轴，总有一段“串行”环节慢下来。这就是 Amdahl 定律描述的现实。此外还有各种“卡顿”要排查：模拟中途崩溃、体系温度暴涨、配体飞出盒子……

性能层面，长程静电是 MD 最贵的单项计算：直接两两求和是 $O(N^2)$。粒子网格 Ewald（PME）用傅里叶网格把它降到 $O(N\log N)$，GPU 加速则让千万原子体系成为日常。调试层面，常见错误包括：力场参数缺失或类型不匹配、水盒子太小导致周期镜像干扰、时间步长过大导致原子“穿模”（能量爆炸）、以及跳过能量最小化就启动模拟。用 profiler 定位瓶颈、按清单逐项排查，是每个计算化学团队的必修课。

$$\displaystyle S(N) = \frac{1}{(1-P) + P/N}$$

1. 并行比例：$P$ 表示整个计算时间中可以被并行处理的部分占比，相当于整场晚会里能同时开多个机位拍摄的分量。

$$\displaystyle P$$

2. 串行比例：$1-P$ 是必须排队完成的部分，比如把多个机位的画面合成一个时间轴，这段只能一个人慢慢做。

$$\displaystyle (1-P)$$

3. 加速比上限：$S(N)$ 表示用 $N$ 个处理单元能获得的最大加速，公式说明即使 $N$ 无穷大，加速也逼近 $1/(1-P)$——串行部分永远卡住天花板。

$$\displaystyle S(N) = \frac{1}{(1-P) + P/N}$$

当机位合理、流程规范，整套“录像—统计—实验”管线才能高吞吐地跑起来。这不仅是学术论文里的方法，更是工业界药物发现流水线上每天发生的日常。

AI 应用场景

Python 演示

import mathimport random  def run_langevin(steps, dt, mass, k_spring, temp, friction=0.5, seed=42):    """一维朗之万动力学: 配体在谐振势阱中采样"""    rng = random.Random(seed)    kb = 0.001987  # 玻尔兹曼常数, kcal/(mol*K)    x = 0.3        # 初始位置    v = 0.0        # 初始速度    sigma = math.sqrt(2 * friction * kb * temp * dt / mass)    xs = []    for _ in range(steps):        force = -k_spring * x        v += (force / mass) * dt - friction * v * dt        v += rng.gauss(0.0, sigma)  # 水的随机碰撞        x += v * dt        xs.append(x)    return xs  def zwanzig_delta_g(xs_a, k_a, k_b, temp):    """用 Zwanzig 方程估算 k_a 到 k_b 的自由能差"""    kb = 0.001987    total = 0.0    n = 0    for x in xs_a:        u_a = 0.5 * k_a * x * x        u_b = 0.5 * k_b * x * x        total += math.exp(-(u_b - u_a) / (kb * temp))        n += 1    return -kb * temp * math.log(total / n)  def analytical_delta_g(k_a, k_b, temp):    """谐振子自由能差解析解: 0.5*kb*T*ln(k_b/k_a)"""    kb = 0.001987    return 0.5 * kb * temp * math.log(k_b / k_a)  def main():    temp = 300.0    xs = run_langevin(20000, 0.002, 1.0, 2.0, temp)    est = zwanzig_delta_g(xs, 2.0, 2.5, temp)    exact = analytical_delta_g(2.0, 2.5, temp)    print("采样帧数    : %d" % len(xs))    print("FEP 估算    : %.4f kcal/mol" % est)    print("解析解      : %.4f kcal/mol" % exact)    print("绝对误差    : %.4f kcal/mol" % abs(est - exact))  if __name__ == "__main__":    main()

1. sigma = math.sqrt(2 * friction * kb * temp * dt / mass)涨落耗散定理：噪声大小与摩擦力、温度绑定，随机碰撞不是乱拍的。
2. v += rng.gauss(0.0, sigma)每步给速度加一个高斯随机增量，模拟水分子对配体的撞击。
3. total += math.exp(-(u_b - u_a) / (kb * temp))对每个采样帧累加玻尔兹曼因子，能量差小的帧贡献大。
4. return -kb * temp * math.log(total / n)把平均权重取对数并乘上 $-k_B T$，得到 A→B 的自由能差。
5. xs = run_langevin(20000, 0.002, 1.0, 2.0, temp)主流程：采样 2 万个帧、弹簧常数 2.0、温度 300 K。
6. exact = analytical_delta_g(2.0, 2.5, temp)解析解是 $0.5 k_B T \ln(k_B/k_A)$，用来验证估算是否靠谱。

常见误区

面试题