4. 聚合物体系能量最小化与预平衡

用分子动力学研究聚合物时，初始构型几乎不可能直接可用。无论是通过自回避随机行走生成的单链、用 Materials Studio 等工具搭建的无定形晶胞，还是从实验密度反推的填充结构，原子之间往往存在严重的重叠、过短的键长、异常的键角，甚至出现两个原子核几乎重合的“硬接触”。这些缺陷会让势能面出现极大的正斥能，一旦直接开始 MD，力会瞬间达到 10³～10⁵ kcal·mol⁻¹·Å⁻¹ 量级，积分器立刻发散，体系“爆炸”。因此，能量最小化与预平衡不是可选的装饰步骤，而是把一条粗糙的初始构型逐步“驯化”为可稳定积分的平衡态的关键工程环节。

为什么不能直接跑 MD

聚合物体系的初始内应力主要来自三处：一是链段生成算法本身不满足真实键长键角分布；二是周期性边界下多条链被随机塞入同一晶胞，链间距离远小于范德华接触半径；三是晶胞密度若直接取实验值，而构型尚未弛豫，会产生巨大的瞬时压力。以常见的聚乙烯无定形晶胞为例，若初始密度设为 0.95 g·cm⁻³ 而构型由随机填充得到，初始势能可能高达 10⁴ kcal·mol⁻¹ 以上，其中 90% 以上是原子对之间的排斥项。此时若用 1 fs 步长直接积分，温度会在几十步内飙升到数千开尔文，键长被拉断，结构彻底破坏。所以必须先“削峰”——把高能接触点通过最小化消除，再逐步让体系在受控条件下弛豫到目标温度和密度。

最速下降法与共轭梯度法

能量最小化的目标是寻找势能面上离初始点最近的局部极小值。最速下降法（Steepest Descent, SD）沿负梯度方向移动，每一步都重新计算梯度，方向简单、鲁棒性强，对初始高能构型特别有效，但接近极小值时会在谷底来回振荡，收敛慢。共轭梯度法（Conjugate Gradient, CG）则利用前一步的方向信息构造共轭方向，在二次型势能面附近收敛速度远快于 SD，适合在 SD 把高能接触消除后做精细收敛。

工程上通常采用“先 SD 后 CG”的两段式策略。以 LAMMPS 为例，典型命令为：

# 第一阶段：最速下降，快速消除原子重叠
min_style sd
minimize 1.0e-4 1.0e-6 1000 10000

# 第二阶段：共轭梯度，精细收敛到局部极小
min_style cg
minimize 1.0e-6 1.0e-8 1000 10000

其中 minimize etol ftol maxiter maxeval 分别控制能量收敛容差、力收敛容差、最大迭代步数和最大力评估次数。判断最小化是否成功，不能只看“是否跑完”，而要看能量曲线是否趋于平台、最大力是否降到 1～10 kcal·mol⁻¹·Å⁻¹ 以下、以及是否出现原子被“冻”在非物理位置。若最小化后仍有大量原子受力超过 100 kcal·mol⁻¹·Å⁻¹，说明初始构型存在严重重叠，需要回到建模阶段检查链填充或密度设置。

退火策略：跨越局部极小

单纯的最小化只能到达离初始点最近的局部极小，聚合物链的缠结、局部取向和链间堆积可能被困在亚稳态。退火（annealing）通过“升温—保温—降温”循环，给体系足够的热能跨越势垒，再在降温中落入更深的极小值。典型退火协议为：从 300 K 升温到 600～800 K（高于聚合物玻璃化转变温度但低于分解温度），保温 100～500 ps，再以 1～5 K·ps⁻¹ 的速率缓慢降温回 300 K，循环 3～5 次。升温阶段让链段获得活动能力，降温阶段让体系在低能构型附近重新组织。对于聚乙烯、聚丙烯等半晶聚合物，退火还能帮助消除初始构型中的虚假取向，使密度和能量更接近平衡值。

逐步升温与密度弛豫流程

退火解决的是“构型采样”问题，而逐步升温与密度弛豫解决的是“热力学状态点”问题。工程上常把两者串成一条完整的预平衡流水线：

  1. 能量最小化：SD + CG，消除硬接触，最大力降到 10 kcal·mol⁻¹·Å⁻¹ 以下。
  2. 低温短程 MD：在 50～100 K 下用 NVT 跑 10～50 ps，步长 0.5 fs，让原子在局部小幅调整，避免高温下键断裂。
  3. 逐步升温：以 0.5～1 K·ps⁻¹ 的速率从 100 K 升到目标温度（如 300 K 或 500 K），每升温 50～100 K 保温 20～50 ps，让体系逐步适应。
  4. NPT 密度弛豫：在目标温度下切换 NPT 系综，压力设为 1 atm，跑 200～500 ps，让晶胞体积自由调整，直到密度和压力趋于稳定。
  5. NVT 平衡：固定弛豫后的密度，在 NVT 下再跑 200～500 ps，监测温度、能量、回旋半径等量是否围绕平均值波动。

判断预平衡是否完成，常用三个指标：势能随时间不再有系统性漂移，只在均值附近波动；体系压力在 NPT 阶段稳定在设定值附近；聚合物链的回旋半径 Rg 和端到端距离 Re 的统计平均值趋于稳定。若 Rg 仍在持续增大或减小，说明链仍在缓慢松弛，需要延长平衡时间。

聚合物体系能量最小化与预平衡流程 ① 能量最小化 SD → CG 消除原子重叠 Fmax < 10 kcal/mol/Å ② 低温短程 MD NVT, 50–100 K 步长 0.5 fs 10–50 ps ③ 逐步升温 0.5–1 K/ps 每 50–100 K 保温 升至目标温度 ④ NPT 1 atm 密度弛豫 200–500 ps ⑤ NVT 平衡 固定密度，200–500 ps 监测 E, T, Rg, Re 退火循环（可选，3–5 次） 300 K → 600–800 K → 保温 100–500 ps → 1–5 K/ps 降温 → 300 K 跨越局部极小，消除虚假取向与缠结 预平衡完成判据 能量判据 势能无系统漂移 仅围绕均值波动 波动幅度 < 2% 压力/密度判据 NPT 压力稳定在 1 atm 密度趋于平台值 波动 < 0.01 g/cm³ 构象判据 Rg、Re 统计平均稳定 链段取向无系统变化 可进入生产 MD

工程化操作中的常见陷阱

第一，最小化容差设得过松。若 ftol 设为 1.0，最大力可能仍有几十 kcal·mol⁻¹·Å⁻¹，后续升温时局部仍会“爆”。建议 ftol 至少 1.0e-4，最好 1.0e-6。第二，升温速率过快。以 10 K·ps⁻¹ 从 100 K 升到 600 K，体系来不及重新分布能量，容易出现局部过热和键断裂。第三，NPT 弛豫时间不足。聚合物链的松弛时间可达纳秒量级，200 ps 的 NPT 可能只完成了局部密度调整，Rg 仍在缓慢变化。此时应延长 NPT 或先做退火再 NPT。第四，忽略周期性边界下的镜像原子重叠。当晶胞尺寸小于截断半径的两倍时，原子可能与自身的镜像发生非物理相互作用，需要在建模阶段保证晶胞边长大于 2×截断半径。

实用建议：对于常见的聚乙烯无定形晶胞（约 5000～10000 原子），一套可复现的预平衡参数为：SD 最小化 5000 步 → CG 最小化 5000 步 → NVT 100 K/50 ps → 逐步升温至 300 K（0.5 K/ps，每 50 K 保温 20 ps）→ NPT 300 K/1 atm/300 ps → NVT 300 K/300 ps。总耗时约 1～2 天（视硬件而定），可得到能量、密度和构象均稳定的初始平衡态。

完成预平衡后，体系才真正具备进入生产 MD 的条件。此时应保存一份平衡态构型（data 文件或 restart 文件），后续所有生产模拟都从同一平衡态出发，以保证不同批次模拟之间的可比性。能量最小化与预平衡看似只是“跑几步”，但它决定了后续所有动力学结果的物理可信度——初始构型的内应力若未消除，后面再长的模拟也只是在错误的状态点上浪费时间。

← 上一章 📖 返回目录 下一章 →