6.2 伞形采样与增强采样路线

6.2 伞形采样与增强采样路线

本节摘要：伞形采样沿已知反应坐标铺一串谐振子窗口，把体系"押送"到能量面上它自己不肯去的地方，再用 WHAM 把各窗口的偏置概率拼回无偏 PMF；metadynamics 则不预设坐标，边采样边填坑。本节讲窗口设计与重叠验收的操作标准，从零实现一个迷你 WHAM 验证算法，并划清两条路线与 QM/MM 的适用边界。选型判据先行：反应坐标知道且要整条 PMF，选伞形；坐标难给或要自由探索，选 meta；要断键成键，直接 QM/MM。

6.1 的炼金术擅长算"差值"，但有一条曲线它给不了：沿某个真实坐标的势量平均力（PMF）——肽键水解的解离曲线、聚合物链穿过膜的牵引曲线、离子跨膜的运输剖面，都要问"沿这条坐标每走一步能量怎么变"。这类问题里势垒是真实存在的物理（不是可以绕开的路径依赖），绕道哲学失效，只能正面攻：把体系押到垒顶去采样。伞形采样与 metadynamics 是两种押送术。

伞形采样：窗口、力常数、重叠

伞形采样的配方三步。第一步定反应坐标与范围：坐标必须是"能概括慢过程"的低维量（两基团距离、二面角、穿过膜的深度），范围从束缚态扫到自由态。第二步铺窗口：沿坐标每隔 s nm 放一个谐振子中心，力常数 k 控制窗口宽度——经验配对是窗口间距 0.1 到 0.15 nm 配 k = 1000 到 2500 kJ/(mol·nm²)，保证相邻直方图在腰部重叠三成以上。第三步每窗口平衡加采样，流程与普通模拟无二（pull 代码拉到窗口中心，或直接从上一窗口末构象起步）。

mdp 关键段比 6.1 还短：

; 每窗口：pull 代码把反应坐标锁在窗口中心
pull                     = yes
pull-ngroups             = 2
pull-group1-name         = CHAIN_A
pull-group2-name         = CHAIN_B
pull-coord1-type         = umbrella      ; 谐振子窗口
pull-coord1-geometry     = distance
pull-coord1-k            = 1500          ; kJ/mol/nm2
pull-coord1-init         = 0.60          ; 本窗口中心，逐窗改
pull-coord1-rate         = 0
​

采样完 gmx wham 一条命令出 PMF；但验收不能只看曲线——直方图重叠图必须先看：相邻窗口直方图在交叠区的高度应同量级，出现"断腰"（两窗各说各话）说明间距太大或力常数太小，PMF 必然有台阶伪影。

迷你 WHAM：从零把直方图拼回去

WHAM 的思想一句话：每个窗口的偏置概率 P_i(ξ) ∝ P_0(ξ)·exp[−β·(½k(ξ−ξ_i)² − F_i)]，反解无偏分布 P_0 并自洽确定各窗口权重 F_i。下面用合成数据走一遍全流程——三个高斯窗口拼出一条带垒的一维 PMF，代码可直接跑：

import numpy as np

rng = np.random.default_rng(3)
KB = 0.0083144621
T = 350.0
beta = 1 / (KB * T)

# 真实 PMF（隐藏答案）：两洼一垒，垒顶在 0.5 nm
xi = np.linspace(0.2, 0.8, 200)
true_pmf = 6.0 * (np.exp(-((xi - 0.28) ** 2) / 0.004)
                  + np.exp(-((xi - 0.72) ** 2) / 0.004)) \
           + 14.0 * np.exp(-((xi - 0.50) ** 2) / 0.02)   # kJ/mol
P0 = np.exp(-beta * true_pmf); P0 /= P0.sum()

# 三窗口采样：中心 0.30/0.50/0.70，k=1500 kJ/mol/nm2，每窗 4000 样本
k, centers, nsamp = 1500.0, [0.30, 0.50, 0.70], 4000
edges = np.linspace(0.2, 0.8, 61)
data = []
for c in centers:
    bias = 0.5 * k * (xi - c) ** 2
    Pb = P0 * np.exp(-beta * bias); Pb /= Pb.sum()
    data.append(rng.choice(xi, size=nsamp, p=Pb))

# 迷你 WHAM：迭代自洽
histograms = [np.histogram(d, bins=edges)[0] / len(d) for d in data]
mids = 0.5 * (edges[:-1] + edges[1:])
f = np.zeros(len(centers))
for it in range(500):
    numer = np.zeros(len(mids))
    for h, c, fi in zip(histograms, centers, f):
        bias = 0.5 * k * (mids - c) ** 2
        numer += h / (np.exp(-beta * (fi - bias)) * h.sum() + 1e-300)
    P0w = numer / numer.sum()
    f = -KB * T * np.log(np.array([
        np.sum(h * np.exp(-beta * 0.5 * k * (mids - c) ** 2) / (
            P0w + 1e-300)) for h, c in zip(histograms, centers)]))
pmf = -KB * T * np.log(P0w + 1e-300)
pmf -= pmf.min()

peak = pmf[np.argmin(np.abs(mids - 0.5))]
well = pmf[np.argmin(np.abs(mids - 0.28))]
print(f"重建垒高 {peak - well:.1f} kJ/mol（设计值 14 附近，采样噪声内有偏差）")
print(f"垒顶位置 {mids[np.argmax(pmf)]:.2f} nm（设计值 0.50）")
​

跑通后你会看到重建的垒高与垒位都在合理偏差内——WHAM 就是一个"把偏置减掉、把权重配平"的自洽平均，没有魔法。对照实验：把窗口减到两个（0.30 与 0.70），垒顶没人采样，重建的垒高会系统性跑偏——这就是重叠验收必须先做的原因。

图 6-2 伞形采样：窗口直方图与拼回的 PMF

metadynamics 与 QM/MM：边界各在哪

metadynamics 的机制是"边走边填坑"：体系沿（若干个）集合变量每积累一段驻留时间，就往那里倒一小坨高斯偏置，洼地被逐渐垫平，体系被迫游历全空间；游历完后累计的偏置量本身收敛到负的 PMF。它的适用边界恰好与伞形互补——反应坐标难给（自由能面是二维以上、路径未知）或怕先入为主（预设坐标可能漏掉真正的慢变量）时用它；代价是收敛难判（要靠多walker、well-tempered 变体与多次独立重跑的一致性），且 GROMACS 生态里通常要挂 PLUMED 这类外挂一起跑。

QM/MM 是另一维度的升级：把发生化学反应的小区域用量子化学（CP2K/ORCA 接口）处理，其余环境照旧经典。它能回答经典力场原理上答不了的问题（断键、电荷转移、光激发），账单也直白：每步算力贵两到三个量级，采样时间被压到皮秒级——所以它只用于"局域化学事件"，并且几乎总是与本章的增强采样联用（QM/MM 讲的伞形采样能垒）才有统计意义。判断题就一句：你的问题里有没有化学键的生成断裂？没有就别碰 QM/MM，先把经典方法的增强采样用足。

速记：伞形三步——定坐标、铺窗口（0.1 到 0.15 nm 配 k 1000 到 2500）、平衡加采样；wham 之前先验收直方图咬合；迷你 WHAM 的思想是减偏置加配权重，两窗不重叠必翻车；meta 适合坐标难给的自由探索，收敛要多次独立重跑背书；QM/MM 只为化学事件买单，且必与增强采样联用。采样工具齐了，下一节解决"跑不跑得起"。