3.6 NVT 与 NPT 平衡实操

3.6 NVT 与 NPT 平衡实操

本节摘要：平衡分两幕：NVT 在定容下把动能温度带到目标值、让约束与恒温器磨合；NPT 在恒压下让盒子呼吸到实验密度、把溶剂化空腔压塌到位。本节给出两段平衡的完整 mdp 与命令会话，讲清位置限制（POSRES）怎么开怎么撤，并给出一套可量化的放行判据——温度、压强、密度三条曲线各自的平稳标准，配一段自动判稳脚本。平衡判据是"敢不敢开始采集数据"的闸门，主观感觉不算数。

EM 给了受力合格的起点，但这个起点还泡在"人为安排"里：分子被摆成了理想构象、盒子体积是随手定的、温度名义上是 0 K 的坐标加随机速度。走台第二、三幕（NVT、NPT）的任务就是把这些人为痕迹冲掉，让体系自己长成平衡态的样子。两幕的分工要说清：NVT 管温度，NPT 管密度——升温在定容下做（压强乱飞没关系），等温度稳了再放压强进来，一步到位同时控两者只会让两套耦合器互相打架。

第一幕：NVT——只管温度

在 3.1 的 nvt.mdp 基础上，会话如下。速度的初始化藏在 gen_vel 两行里：

gen_vel                 = yes        ; 生成初始速度
gen_temp                = 350        ; 按目标温度的麦克斯韦分布采样
gen_seed                = 173529     ; 随机种子：写进记录，复现靠它
​

$ gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
$ gmx mdrun -deffnm nvt
​

两个细节常被忽略。-c 与 -r 都填 em.gro：-c 是坐标起点，-r 是位置限制的参考位置（define = -DPOSRES 激活拓扑里的 [ position_restraints ]，把聚合物重原子用 1000 kJ/mol/nm² 的弹簧拴在 EM 后的位置上——弹簧常数就写在 posre.itp 里）。限位的目的：升温初段溶剂会在聚合物周围重新排布，若不加限位，链可能整体平移旋转"蹭"走溶剂化功，测得的平衡状态里混进了非物理的质心漂移。生产段必须撤掉 POSRES（3.1 说过的差异之一），否则你测的扩散系数是被拴着的假数。

第二处是压强的预期行为：NVT 段 pcoupl = no，压强会以几十到几百 bar 的幅度大幅涨落——这不是错误，是定容下的正常物理（初始盒子偏大则负压、偏小则正压）。新手在这里的恐慌占 NVT 阶段提问量的大头，提前打针。

第二幕：NPT——把密度交还物理

NPT 的 mdp 在 NVT 基础上改三行（其余不动）：define 删掉或保留看你的方案（聚合物研究常在 NPT 前半段保留限位）；pcoupl 换上 3.3 选定的 C-rescale；nsteps 按 500 ps 到 1 ns 给足。会话照抄：

pcoupl                  = C-rescale
pcoupltype              = isotropic
tau_p                   = 5.0
ref_p                   = 1.0
compressibility         = 4.5e-5
​

$ gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
$ gmx mdrun -deffnm npt
​

-t nvt.cpt 把速度与耦合器状态从 NVT 无缝接过来（比只接坐标干净得多，4.1 会展开 checkpoint 的全部用途）。NPT 段要盯的曲线是密度：2.3 用 -d 1.0 搭的盒子，初装水接近 1 g/cm³ 但溶剂化壳层未弛豫，前几十 ps 密度会明显爬升或回落，最终平台值就是你的体系在 350 K 下的真实密度——对 PEO 水溶液，落在 1.0 到 1.1 g/cm³ 区间都对得上物性表。

图 3-5 NVT 到 NPT：三条曲线的放行判据

演算：把"平稳"变成数字

凭"看着稳了"放行是目测迷信，判据要能写进脚本。三条曲线的量化标准：温度——末段（如最后一半轨迹）均值与目标的偏差小于 2 到 3 K；密度——末段均值与前段均值的漂移小于 0.5 到 1 个百分点；能量——势能末段趋势线的斜率与涨落之比足够小（漂移速度远慢于涨落幅度）。下面的脚本读 gmx energy 导出的 xvg，自动执行这三条检查：

import numpy as np, sys

def load_xvg(path):
    data = []
    for line in open(path, encoding="utf-8"):
        if line[0] in "@#":          # 跳过图例与注释
            continue
        data.append([float(x) for x in line.split()])
    return np.array(data)

def equilibrium_check(path, columns=(1, 2, 4), names=("T", "P", "rho")):
    d = load_xvg(path)
    t = d[:, 0]
    half = len(t) // 2
    print(f"总时长 {t[-1]-t[0]:.0f} ps，取后一半判定")
    verdict = []
    for col, name in zip(columns, names):
        if col >= d.shape[1]:
            continue
        y_first, y_late = d[:half, col], d[half:, col]
        drift = abs(y_late.mean() - y_first.mean()) / max(1e-12, abs(y_first.mean()))
        slope = np.polyfit(t[half:], y_late, 1)[0]
        spread = y_late.std()
        slow = abs(slope) * (t[-1] - t[half]) / max(spread, 1e-12)  # 漂移/涨落
        verdict.append((name, y_late.mean(), drift, slow))
        print(f"{name:>4}: 末段均值 {y_late.mean():10.3f}  "
              f"前后漂移 {drift*100:5.2f} %  漂移涨落比 {slow:.3f}")
    return verdict

# 用法：gmx energy 导出含 T/P/rho 的 xvg 后
verdict = equilibrium_check(sys.argv[1] if len(sys.argv) > 1 else "npt_energy.xvg")
fails = [v for v in verdict if (v[2] > 0.01 and v[0] != "P") or v[3] > 0.1]
print("放行" if not fails else f"回炉：{[f[0] for f in fails]} 漂移过大")
​

压强在判定里被豁免漂移检查（第 2 列均值百 bar 级摆动是正常物理），但生产期报告压强均值时要用长轨迹取平均并报出标准误差——这是 7.2 质控清单的条目之一。

易错点与本节速记

平衡段的三个高频失误：其一，NVT 直接上生产长度却不看温度曲线，恒温器与约束磨合期（前 50 ps 左右）的坏数据混进平衡统计；其二，NPT 忘记带 -t 续 checkpoint，速度从高斯分布重新采样，温度曲线出现一个尖峰起步；其三，posres 忘撤就跑生产，扩散系数与构象采样双双缩水。此外记住一句行话：平衡判据防的是系统性漂移，不是涨落——涨落是物理，漂移是病。

速记：先 NVT 后 NPT，两幕各司其职；-r 与 POSRES 在平衡期拴链、生产期必撤；-t 接力 checkpoint 是无缝续算的正道；放行判据三条数字化——温度均值贴目标、密度漂移小于一个百分点、能量漂移涨落比小于一成；判稳脚本进仓库，拒绝目测放行。走台到此合格，下一章正式开演。