7.1 聚合物熔体案例复盘

7.1 聚合物熔体案例复盘

本节摘要：本节把一个 PEO 熔体课题从头到尾走一遍：50 条 20 单体链的联合原子模型，退火平衡解决慢松弛，链尺度统计（回转半径、末端距、RDF）避开周期性假象。每个决策点都标注"为什么这么选、不这么选会怎样"。聚合物是检验前六章功底的最好考题——它没有现成 pdb2gmx 流程、平衡慢到以十纳秒计、链与自己的镜像纠缠，每一步都在考基本功。

先立题。研究问题设为："中等分子量 PEO 熔体的链尺度统计与局域堆积结构"。这个题目小到两周能跑完、大到能牵出前六章的每一章——是理想的复盘载体。体系规格沿用 2.1 与 2.2 定下的方案：20 单体链（约 880 g/mol），50 条链共 3000 个联合原子格点，4.0 nm 立方盒子（2.1 演算过的密度账），TraPPE 风格联合原子参数。

决策点一：装配——没有 pdb2gmx 的世界

聚合物走的是 2.2 的手写 itp 路线：单体参数化一次、脚本生成 50 条链、拓扑 [ molecules ] 区一行 "PEO20 50"。复盘时值得强调三个当时容易犹豫的决策。其一，初始构象用随机行走生成而非伸直链：50 条伸直链（每条伸直约 5.3 nm）塞进 4 nm 盒子必然互相穿插过密，EM 会把重叠推开却留下大量人为张力；自避随机行走生成的坍缩线团初态温和得多。其二，盒子不用留 1.2 nm 的溶剂化边距——熔体没有溶剂，-d 的意义只剩周期性假象防护，链的回转半径约 1 nm，盒子 4 nm 意味着链与其镜像的最小距离约 2 nm，超过两倍 Rg，够用。其三，EM 后必须做退火（下一节），因为装配态离平衡隔着好几个"链松弛时间"。

决策点二：平衡——退火是熔体的正道

聚合物熔体的平衡是全册最慢的平衡：链构象的松弛时间随分子量增长极快（蛇行理论给指数律），20 单体链的端到端重排就要若干纳秒，装配初态的"随机行走记忆"要靠真实动力学洗掉。纯等温平衡效率太低，标准做法是升温退火：把体系短暂加热到远高于使用温度（如 500 到 600 K），链运动加快、拓扑纠缠重排，再降温回到目标温度（如 413 K，PEO 熔体的常用研究温度），最后等温平衡到 3.6 的判据全过：

; anneal.mdp 关键段：500 ps 内从 550 K 线性降到 413 K
annealing        = single
annealing-npoints = 3
annealing-time   = 0 100 500
annealing-temp   = 550 550 413
ref-t            = 413
dt               = 0.002
nsteps           = 250000          ; 500 ps 退火
​

退火有一个必须自查的副作用：高温段可能"蒸发"掉构象空间的其他 basin（比如把某些局部有序结构熔掉）——自查方法是把退火后的密度与文献熔体密度对表（PEO 熔体约 1.06 到 1.13 g/cm³ 随温度变化），密度对上了，退火就没跑偏。等温平衡段再接 5 到 10 ns，用 3.6 的脚本判稳，用 5.2 的 RMSD 确认平台（熔体的 RMSD 平台高度天然较高——链构象离初始态很远，别拿蛋白的 0.2 nm 直觉来读）。

决策点三：统计——链尺度量与周期性假象

熔体的观测量分两层。局域结构（5.3 三件套直接适用）：O-O 与 CH2-O 的 RDF 给出堆积壳层；特征峰位与文献联合原子数据对表是力场校验的硬指标。链尺度统计才是聚合物课题的主角：回转半径 Rg、末端距分布、链的内取向关联。Rg 的定义与逐链计算口径：

Rg² = (1/N) Σ_i | r_i − r_cm |²

Python 按帧逐链计算（不重复造轮子的部分：读轨迹用 MDAnalysis 一类库，此处给核心算法演示全质心与回转半径的口径）：

import numpy as np

def rg_of_chain(coords):
    """coords: N×3，一条链的联合原子坐标。返回 Rg（nm）。"""
    cm = coords.mean(axis=0)
    return np.sqrt(((coords - cm) ** 2).sum(axis=1).mean())

# 自测：高斯链的标度律 Rg ∝ N^0.588（良溶剂）到 N^0.5（理想链）
rng = np.random.default_rng(11)
for n_beads in (10, 20, 40):
    steps = rng.normal(size=(n_beads, 3))
    chain = np.cumsum(steps, axis=0) * 0.15
    print(f"N={n_beads:3d}  Rg={rg_of_chain(chain):.3f} nm")

# 轨迹统计的正确姿势：逐帧逐链算 Rg，再对帧取均值与误差棒
# 错误姿势：把 50 条链的原子混在一起算一个"总 Rg"——那不是任何物理量
​

周期性假象在这个环节最会咬人：若某条链在轨迹里因 PBC 被切开（5.1 的病），它的 Rg 会被算成两个碎片的跨盒距离——数值巨大且随帧跳变。验收手段：逐链 Rg 的时间序列里若有离群链（比中位数大好几倍），先回去补 -pbc whole 再重算。另一个假象是"自相互作用"：链长超过半盒宽时链会与其镜像"自己也算非键"，Rg 被人为压扁——盒子够大的设计（决策点一）就是防它。

决策点四：产能与收尾

按 6.3 的口径收个尾：3000 原子的联合原子体系在单节点（64 核加一卡）典型产能 100 到 300 ns/天——这就是聚合物圈偏爱联合原子的现实理由。全案预算：退火 0.5 ns + 平衡 10 ns + 生产 200 ns（链尺度统计的误差棒要够窄），单节点一周内可交卷；若要研究缠结（需链长上一个量级），产能账立刻翻两番，6.2 的粗粒化路线就该上桌了。

演算：末端距分布与特征比

除了 Rg，聚合物课题的第二个标配量是末端距 Re 与特征比 C∞ = ⟨R²⟩ / (N·b²)（N 为格点数、b 为键长），它把"链有多伸展"归一成可与理论对比的数。演算口径：

import numpy as np

rng = np.random.default_rng(23)
# 三种链的 Re/Rg 对比：理想链 Re/Rg = sqrt(6)
def chain_stats(coords):
    cm = coords.mean(axis=0)
    rg = np.sqrt(((coords - cm) ** 2).sum(axis=1).mean())
    re = np.linalg.norm(coords[-1] - coords[0])
    return rg, re

# 自由行走参考
walk = np.cumsum(rng.normal(size=(59, 3)), axis=0) * 0.15
rg, re = chain_stats(walk)
print(f"自由行走: Rg={rg:.2f} Re={re:.2f} Re/Rg={re/rg:.2f}（理想值 {np.sqrt(6):.2f}）")

# 轨迹统计口径：逐帧逐链算 Re/Rg，报告分布而非单值
# 报告格式：均值 ± 块平均标准误（6.3 的产能换算决定了样本够不够）
​

判读基准：熔体是理想链环境（排除体积被链间屏蔽），熔体内链的 Re/Rg 应接近 2.45；明显偏大说明链被人为拉伸（初态记忆未洗掉——回炉退火），明显偏小说明链坍缩（温度或密度设错）。一个数就能给平衡质量打折或盖章。

全流程时间线复盘

把本案例的执行顺序与耗时立成表，方便你替换成自己的体系：

这张表真正的教学点是平衡与生产的时长比：熔体课题里平衡段常常吃掉总预算的一成到三成——"平衡很贵"是聚合物圈与蛋白圈最大的文化差异，蛋白圈 10 分之一量级的直觉在熔体上不够用。

两个高频疑问

联合原子模型算 Tg 靠谱吗？ 粗略趋势可以、绝对值要小心：联合原子抹掉了氢键与部分扭转势垒，玻璃化温度通常系统性偏移几十开尔文。要报 Tg 建议同口径下做相对比较（不同链长、不同端基），并与全原子短模拟或实验锚一个位移量——绝对直报是这类模型最常见的滥用。

链间缠结怎么看？ 定性看拓扑约束对链运动的限制（蛇行），定量算缠结长度需要本体涨落或纯几何（Z1 类算法）分析，属进阶课题。入门阶段更实际的替代问题是"我们的链长够不够 Entangled"：先估缠结分子量（PEO 约 2000 到 3000 g/mol 量级），再决定 20 单体链（880 g/mol）够不够格谈缠结——本案例刻意选在缠结线以下，统计更干净。

速记：聚合物没有 pdb2gmx，手写 itp 加脚本复用；初始态随机行走、盒子按 Rg 留、溶剂化边距概念不适用；平衡靠退火，密度对表验退火没跑偏；Rg 逐帧逐链统计，混链总 Rg 不是物理量；离群 Rg 回查 PBC，链长过半盒宽自查自相互作用。熔体考完了耐心，下一节换个考细心的：脂膜。