3.5 能量最小化实操

3.5 能量最小化实操

本节摘要：能量最小化（EM）把体系构建期留下的高能接触、键角畸变逐一卸掉，为动力学模拟提供一个受力合理的起点。本节讲清最速下降法的流程与收敛判据（最大力低于 1000 kJ/mol/nm），给出从 ions 前到 EM 完成的完整会话，演示收敛曲线的三种典型形态与各自的处置，并配一段解析 log 提取收敛历史的检查脚本。EM 是一切崩溃排查的分诊台——体系有问题，几乎都在这里提前暴露。

参数（3.1 至 3.4）都定了，走台第一幕是把结构"顺平"。为什么必须先做 EM？回想 1.2 的演算：一根被拉伸 0.01 nm 的 C–C 键携带约 28 kJ/mol 的形变能，solvate 时被挤进缝隙的水分子与溶质的近距离接触更是携带成百上千 kJ/mol 的排斥能。带着这种初始力直接积分，第一步加速度就能把原子踢出盒子——EM 的目的不是找全局极小，只是把局部的大梯度削平，让第一步动力学不至于起飞。

mdp 与会话：最速下降一遍过

EM 用最速下降（integrator = steep）：沿当前合力方向走一步，梯度为零的方向不动，直到最大原子受力低于阈值。它不追求收敛到精密极小值（共轭梯度更擅长那个，但开局反而慢），只求把粗暴的力快速压下去：

; em.mdp —— 走台第一幕的全部参数
integrator   = steep        ; 最速下降
emtol        = 1000.0       ; 最大力阈值 kJ/mol/nm：低于它就算收敛
emstep       = 0.01         ; 初始步长 nm，程序自适应增减
nsteps       = 50000        ; 步数上限，防止病态体系无限磨
nstlist      = 10
cutoff-scheme = Verlet
coulombtype  = PME
rcoulomb     = 1.0
rvdw         = 1.0
constraints  = none         ; EM 阶段放开约束，让键角也参与弛豫
pbc          = xyz
​

constraints = none 这行值得单独说明：动力学阶段全约束含氢键是为了换步长，EM 阶段没有步长压力，放开约束反而让被实验结构"冻结"的键长键角一起弛豫，收敛质量更好。

$ gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr
$ gmx mdrun -deffnm em
Steepest Descents:
   Tolerance (Fmax)   =  1.00000e+03
   Number of steps    =        50000
Step  142, time  2.840e-01
...
Potential Energy  = -2.64561e+05
Maximum force     =  8.72112e+02 on atom 3175
Norm of force     =  6.88431e+01

$ printf "Potential\nPressure\n" | gmx energy -f em.edr -o em_energy.xvg
​

三行输出各有含义：Potential Energy 为大负值（水化的体系通常每原子负几百 kJ/mol 量级）说明没有悬空的高能接触；Maximum force 低于 1000 即达 emtol 判据；报告里的 atom 3175 是最大力原子的编号——若 EM 迟迟不收敛，去结构里查它是什么，十有八九是构建期埋的雷（水钻进了环空腔、离子贴着骨架、脂膜头基翻转卡位）。

收敛曲线的三种形态与处置

把 em_energy.xvg 画出来（或用 3.6 的脚本扫一遍），典型形态三种：

1. 教科书式收敛：势能陡降几个量级后平掉，最大力指数下降到阈值以下，几百到几千步完成。处置：直接放行，进 NVT。
2. 平台后跳崖：势能先卡在某个平台，几千步后突然再降一段——通常是某个水分子或离子终于"挤"过了窄缝。处置：正常现象，若反复出现，检查盒子是不是太小（2.3 的 -d 不足）。
3. 迟迟不达标：最大力卡在 2000 至 5000 徘徊不降。处置：按顺序查——最大力原子是谁；体系有没有原子重叠（gmx check 与可视化）；PME 网格有没有报异常；实在不行降 emstep 到 0.001 重新起步，或局部手工把雷挪走再跑。

图 3-4 能量最小化收敛曲线判读

收敛检查脚本

会话输出只在最后报一次最大力，要画曲线、要留档，得自己解析。下面的脚本读 mdrun 的 log 文件，抓出每百步的势能与最大力，输出达标步数与衰减率，可直接进课题仓库：

import re, sys

def parse_em_log(path="em.log"):
    steps, pot, fmax = [], [], []
    pat_step = re.compile(r"^\s*Step\s+(\d+)")
    pat_pot  = re.compile(r"Potential Energy\s*=\s*([-\d.e+]+)")
    pat_fmax = re.compile(r"Maximum force\s*=\s*([-\d.e+]+)")
    step = pot_v = fmax_v = None
    for line in open(path, encoding="utf-8", errors="ignore"):
        m = pat_step.match(line)
        if m:
            step = int(m.group(1))
        m = pat_pot.search(line)
        if m:
            pot_v = float(m.group(1))
        m = pat_fmax.search(line)
        if m:
            fmax_v = float(m.group(1))
            if step is not None and pot_v is not None:
                steps.append(step); pot.append(pot_v); fmax.append(fmax_v)
    return steps, pot, fmax

steps, pot, fmax = parse_em_log(sys.argv[1] if len(sys.argv) > 1 else "em.log")
target = 1000.0
ok = next((i for i, f in enumerate(fmax) if f < target), None)
print(f"共 {len(steps)} 个记录点")
if ok is not None:
    print(f"第 {steps[ok]} 步达标：max force {fmax[ok]:.1f} < {target}")
    print(f"势能从 {pot[0]:.3e} 降到 {pot[ok]:.3e}（降幅 {pot[0]-pot[ok]:.3e}）")
else:
    print("未达标：", f"最后 max force {fmax[-1]:.1f}" if fmax else "log 里没有力记录")
    # 半衰期式诊断：力每千步降多少
    if len(fmax) > 2 and fmax[0] > 0:
        rate = (fmax[-1] / fmax[0]) ** (1000 / max(1, steps[-1] - steps[0]))
        print(f"每千步力残余比约 {rate:.3f}；高于 0.9 视为近乎停滞，优先查最大力原子")
​

最后那句诊断是给形态三准备的：残余比接近 1 说明最速下降在狭谷里打转，此时换 emstep 或手工拆雷都比死等 50000 步体面。

速记：EM 只削梯度不找极小；steep 加 emtol 1000 加 nst 上限是全套判据；最大力原子编号是排查的路标；收敛曲线三种形态对应三种处置；EM 不过关绝不开动力学——它是后面所有环节的分诊台。起点合格了，下一幕升温升压：NVT 与 NPT。