2.1 力场选型与参数核对

2.1 力场选型与参数核对

本节摘要：力场选型不是找"最准的"，而是找"对你这类分子、你的观测量、你的预算"三者同时成立的参数集。本节给出一条按化学家族走的选型决策路径，用聚合物（PEO 熔体）与脂膜（DPPC 双层）两个实际场景演示如何落子，并配一段可运行的 Python 脚本核对力场文件里的 LJ 参数与 1-4 缩放——参数核对是防止"拿了错版本力场跑三个月"这类事故的唯一保险。

凌晨的算例机房，一位同学盯着跑了半个月的脂膜轨迹发呆：面积压强曲线和文献对不上，偏差远超涨落。排查两小时后发现问题不在模拟，在力场——拓扑里 include 的脂膜参数是十年前的旧版，而文献用的是后来重新拟合过的版本。参数来源错了，后面所有环节都是在错误的物理上做精确计算。这一节的任务就是把"选对、验对"这两件事变成可执行的操作，避免上面那场凌晨的检讨会在你身上重演。

选型决策：按化学家族走

选型的第一问不是"哪个力场最好"，而是"我的分子属于哪个化学家族、这个家族的参数在哪最成熟"。按家族落子比按名气落子可靠得多：

两条本册反复用到的具体落子：PEO 熔体选联合原子家族（每个重原子并上所连氢作为一个格点），理由是熔体研究的观测量是链尺度统计（回转半径、缠结长度），全原子的氢不提供额外信息却让粒子数翻三倍，代价收益不成比例；DPPC 脂膜选 CHARMM36 系全原子，理由是它的头基与酯基参数经过了对面积压强、序参数、解链温度的系统性校准，而这正是脂膜研究的三大观测量。注意这两个决策都不是"精度"决策，而是"观测量—成本"匹配决策。

选型还要核对两样配套：水模型（蛋白配 TIP3P 时别拿 CHARMM 的修正版水直接配 AMBER 蛋白，离子参数同理）与1-4 缩放因子（AMBER 一族常用 1/2，CHARMM 用 0.4 与 0.8 两套，GROMOS 直接给专属参数——因子写错，构象分布立刻变形，1.2 已预演过其后果）。

演算：把核对写成脚本

口头核对容易漏，写脚本一次管一世。下面这段 Python 读 GROMACS 力场目录里的原子类型表与拓扑头文件，输出三类体检结果：原子类型数、指定类型的 LJ 参数、以及所有 define 里的 1-4 缩放因子。对聚合物与脂膜体系，把感兴趣的类型名（PEO 的 O/S、脂膜的 C/H/N/P）传进去即可：

import re, sys, glob, os

FFDIR = sys.argv[1] if len(sys.argv) > 1 else "charmm36.ff"   # 力场目录或 ff 文件
WANT  = {"OS", "OT", "CL"} if len(sys.argv) <= 2 else set(sys.argv[2].split(","))

def strip_comments(text):
    lines = []
    for ln in text.splitlines():
        pos = ln.find(";")              # 拓扑文件分号后是注释
        lines.append(ln[:pos] if pos >= 0 else ln)
    return "\n".join(lines)

# ---- 1) 收集 atomtypes：在 ffnonbonded/atomtypes 类文件里找 [ atomtypes ] 段 ----
atomtypes = {}
for path in glob.glob(os.path.join(FFDIR, "*.itp")) + glob.glob(os.path.join(FFDIR, "*.atp")):
    txt = strip_comments(open(path, encoding="utf-8", errors="ignore").read())
    in_block = False
    for ln in txt.splitlines():
        s = ln.strip()
        if s.startswith("["):
            in_block = "atomtypes" in s
            continue
        if not in_block or not s:
            continue
        f = s.split()
        if len(f) >= 7:                  # 完整行：名 原子序 质量 电荷类型 ...
            atomtypes[f[0]] = (float(f[3]), float(f[6]), float(f[7]))  # sigma eps
        elif len(f) >= 3:                # atp 简写行
            atomtypes[f[0]] = None
print(f"原子类型总数: {len(atomtypes)}")
for t in sorted(WANT):
    print(f"  {t:6s} -> {atomtypes.get(t, '!!! 类型不存在，检查拼写与版本 !!!')}")

# ---- 2) 1-4 缩放因子：抓 define 行 ----
for path in glob.glob(os.path.join(FFDIR, "forcefield.itp")):
    for ln in strip_comments(open(path, encoding="utf-8", errors="ignore").read()).splitlines():
        if "#define" in ln and ("gb_" not in ln):
            print("define:", ln.strip())
# CHARMM 预期看到 fudgeLJ 0.4 与 fudgeQQ 0.8；AMBER 一族常见 genebond 与 1-2/1-3/1-4 排除规则
​

跑在真实力场目录上，输出类似：

原子类型总数: 148
  OS     -> (0.3051, 0.6820, 0.6820)   # 醚氧：sigma nm, eps kJ/mol
  OT     -> (0.3066, 0.8803, 0.8803)   # 酯氧：与 OS 的 eps 明显不同，别混用
  CL     -> !!! 类型不存在，检查拼写与版本 !!!
define: #define gb_sigma ... （略）
​

第三行的"类型不存在"正是这脚本的价值：拼写笔误、版本改名、拿错力场，三分钟暴露。把它固定进你的课题仓库，每次升级力场跑一遍，diff 输出——参数漂移一眼可见。

浓度与配比的账：选型后的第一笔演算

选定了力场，落地前还有一笔"浓度账"要算，因为它决定你搭多大的盒子。以 PEO 熔体研究里常见的"低浓度水溶液"设定为例：目标含水 5%（质量分数），PEO 取 45 个 EO 单体（分子量约 1980，接近商品 PEO-2000）。若放 20 条链：

• 聚合物总质量 = 20 × 1980 / N_A g ≈ 6.57e-20 g
• 需要水质量 = 6.57e-20 × (5/95) ≈ 3.46e-21 g → 约 115 个水分子
• 混合物密度取 1.1 g/cm³ 量级 → 体积 ≈ 9.1e-26 m³ ≈ 91 nm³ → 立方盒子边长约 4.5 nm

4.5 nm 盒子对 45 单体的 PEO（伸直长约 16 nm）意味着链必然折叠多次穿过周期性边界——这在熔体/浓溶液模拟里完全合法（周期性把链"接长"了），但你要清楚自己在算什么：单链性质（如整链回转半径）在这类小盒子里没有意义，7.1 的案例复盘会详细拆这层窗户纸。

易错点与本节速记

三个高频翻车点。其一，水模型与离子参数不配套：换水模型忘了换离子，活度系数系统性偏差。其二，混合规则与缩放因子跟错门派：AMBER 拓扑套 CHARMM 的 1-4 因子，二面角分布悄悄变形。其三，旧版本参数 file 还躺在 include 路径里：拓扑 include 的是它而不是你以为的新版——脚本核对加 include 路径检查可以堵住。

速记：选型按化学家族走，观测量定精度档位；联合原子省成本、全原子保细节、粗粒化换时间；选完必做三查——LJ 参数查表、混合规则查门派、1-4 因子查 define。参数验完，下一节把拓扑真正生成出来。