3.2 积分算法与步长选择

3.2 积分算法与步长选择

本节摘要：分子动力学的积分器只做一件事：用当前位置的力把坐标与速度推进 dt。本节从"为什么欧拉法会失败"讲起，推导 Velocity Verlet 与蛙跳法的同构关系，说明 MD 界偏爱辛积分器的能量守恒理由；再用一段可运行的 Python 演算对比三种方法的能量漂移；最后推导步长上限——为什么 O–H 振动把无约束步长压到 1 fs，而 LINCS 约束把它解放回 2 fs。

3.1 的 mdp 里 integrator = md 与 dt = 0.002 这两行是全文件的物理起点，但"为什么是这两行"还没有兑现。本节兑现它。承上：1.1 说每步只做"算力、推坐标、更新速度"；启下：步长的推导直接决定 3.1 里 constraints = h-bonds 这一行是不是可有可无（它不是）。

从欧拉法的失败说起

最朴素的数值积分是显式欧拉：v(t+dt) = v(t) + a(t)·dt，r(t+dt) = r(t) + v(t)·dt。它在绝大多数物理课上第一次出现，在 MD 里第一次出场就该退场——因为它是非辛的：每一步都不守相空间面积，系统总能量单调漂移。对要跑百万步的模拟，漂移不发散也会把温度"漂"到目标的几倍开尔文之外。

MD 的标准答案是 Verlet 家族，核心技巧只有一条：用"旧速度 + 新位置的力"的折中来更新速度。Velocity Verlet 的一步是：

• v(t+dt/2) = v(t) + a(t)·dt/2（半步速度）
• r(t+dt) = r(t) + v(t+dt/2)·dt（推坐标）
• 用新位置算新力 a(t+dt)
• v(t+dt) = v(t+dt/2) + a(t+dt)·dt/2（补全速度）

蛙跳（leap-frog）与它是同一枚硬币的两面：蛙跳干脆只维护半整数时刻的速度表，GROMACS 的 md 积分器用的就是蛙跳记账、Velocity Verlet 语义，输出时再换算回整数时刻速度——所以你在 log 里看到的温度永远与全速度对应，不必自己换算。

辛性带来的实际好处有一条朴素判据：长时能量误差有界且不漂移。真能量守恒与辛误差的区别在下面的演算里看得清清楚楚。

演算：三种积分器的能量漂移实测

拿一维谐振子当小白鼠（U = ½kx²，任何多体体系里最硬的键长项都是它），分别用显式欧拉、改进欧拉（半步速度）与 Velocity Verlet 积分 5000 步，看总能量漂移：

import numpy as np

k, m = 100.0, 1.0          # 力常数与质量（无量纲单位制）
dt = 0.05                  # 相当于真实 O-H 振动的"激进"步长
x0, v0 = 1.0, 0.0          # 初始位移 1，初速 0

def energy(x, v):
    return 0.5 * m * v ** 2 + 0.5 * k * x ** 2

def run(method, nsteps=5000):
    x, v = x0, v0
    E0 = energy(x, v)
    traj = []
    for i in range(nsteps):
        a = -k * x / m
        if method == "euler":                    # 显式欧拉：注定漂移
            v += a * dt
            x += v * dt
        elif method == "verlet":                 # Velocity Verlet：辛积分器
            v_half = v + 0.5 * a * dt
            x += v_half * dt
            a_new = -k * x / m
            v = v_half + 0.5 * a_new * dt
        traj.append(energy(x, v) / E0)
    traj = np.array(traj)
    return traj

for name in ("euler", "verlet"):
    t = run(name)
    print(f"{name:>7}: 末态能量比 {t[-1]:10.3f}   峰值 {t.max():8.3f}   谷值 {t.min():8.3f}")
​

典型输出：欧拉法的能量比冲到几十甚至发散；Verlet 法的能量比始终在 1 ± 0.02 内小幅振荡、不漂移。把 dt 放大到 0.12（约等于振动周期的五分之一）再跑，Verlet 的振荡幅度明显变大——这条"幅度随 dt 增大"的边界就是稳定性区，把不同 dt 下的振幅画出来就是下图：

图 3-1 步长与稳定区：能量误差如何随 dt 爬坡

步长上限的推导：从振动周期到 2 fs

经验红线"dt ≈ 周期的十分之一"可以推出来。谐振子的稳定条件（对 Verlet）是 ω·dt < 2，实际使用留足余量取 ω·dt ≈ 0.2 到 0.3。用真实数字走一遍：

• O–H 伸缩波数约 3000 至 3400 cm⁻¹，取 3300 cm⁻¹；周期 T = 1/(c·ν̃) = 1/(3e10 cm/s × 3300 cm⁻¹) ≈ 1.0e-14 s = 10 fs
• 无约束时的安全步长：0.2 × 10 fs = 2 fs？——不行，注意约束的真正作用：ω·dt 的"ω"由最硬振动给出，无约束时就是这根 O–H 键；理论极限虽是 2/ω ≈ 3.2 fs，但误差振荡在接近极限时已不可接受，行业实践取 1 fs
• 全约束含氢键后，最快振动变成 C–O–C 类角振动（约 700 至 900 cm⁻¹，周期 37 至 48 fs）→ 安全步长 0.2 × 37 fs ≈ 7 fs，取 2 fs 仍有约三倍余量

所以"constraints = h-bonds 加 dt = 0.002"不是两行独立参数，而是一笔绑定交易：用约束算法冻住最硬的振动，换取步长翻倍。约束算法本身（LINCS 迭代解约束方程、SETTLE 解析解算刚性水）每步多花一点算力，换来产能净翻倍，稳赚。LINCS 的阶数（lincs-order）默认 4，约束力方向修正到四阶精度；若 4.3 的排查手册里见到 LINCS 警告，第一反应不是调阶数，而是回头查 3.5 的最小化做没做够。

时间可逆性：为什么这点值得在意

Verlet 家族还有一个容易被忽略的性质：时间可逆——把所有速度反号，轨迹会原路退回。配上可逆的恒温恒压算法（3.3 的 Nosé-Hoover 一族），整个模拟在统计力学意义上严格对应一个可逆系综，自由能微分式的误差项才可控。这也是 3.3 表格里"Berendsen 不可逆、只配平衡期"的根源：不可逆的耦合破坏系综的严格性，平衡期没人追究，生产期必须换掉。

积分器选型的完整菜单

Velocity Verlet 之外，mdp 里还能遇到的积分器各有一席之地，选型按问题性质走：

粗粒化体系（如 MARTINI）有个例外值得知道：粗粒化把最快振动抹掉了，dt 可以放大到 0.02 到 0.04 ps，这也是粗粒化"换时间"的红利来源——但 sd 积分器的摩擦系数与随机力要配套调，动力学量（扩散系数）的解读要乘回因子，不是"白送的时间"。

两个高频疑问

步长 2 fs 上下，产能差一倍，能不能偷偷用 2.5 fs？ 偷偷不行，明着也难：稳定边界由最硬的未约束振动决定，2.5 fs 已越过多数全约束体系的保险区，温度会缓慢漂高、约束告警会在几十纳秒后随机爆发——你买到的那 25% 步长会在排错里加倍还回去。真要省机时，去 6.3 谈并行与参数。

温度耦合与积分器的顺序有讲究吗？ 有但程序替你管了：GROMACS 每步先算力、积分、再施加约束、最后做耦合修正，协变量的顺序固定且经过长期验证。你唯一要做的决定是 3.3 的"选哪个耦合算法"，别试图手搓"每几步手动重标一次速度"这类私货——那是在重新发明一个不受支持的恒温器。

速记：欧拉法能量漂移，Verlet 家族有界振荡；蛙跳与 Velocity Verlet 同构，GROMACS 的 md 就是它；dt 上限由最快振动周期定，O–H 约 10 fs 压出无约束 1 fs；LINCS/SETTLE 冻住含氢键，换来 2 fs 与产能翻倍；可逆性是自由能计算的隐性前提。步长定了，下一节看恒温恒压怎么选算法。