3. 聚合物体系初始构型生成工程

分子动力学模拟的可靠性，从第一帧构型就已经被决定。如果初始结构里存在原子重叠、链穿越或密度严重偏离平衡值，后续即使使用再小的积分步长，也会在能量最小化阶段出现巨大受力、体系爆炸或虚假交联。本章围绕单链构建、无规线团生成、熔体平衡密度填充、周期性边界条件设置，以及避免原子重叠与链穿越的实操方法展开，目标是让读者能够独立生成一个“物理上合理、数值上稳定”的聚合物初始构型。

一、单链构建：从拓扑到坐标

单链构建的第一步不是写坐标，而是确定链的拓扑与连接关系。以聚乙烯（PE）为例，若聚合度 N=100，则主链碳原子数为 100，氢原子数为 202，总原子数 302。需要明确每个原子的类型、成键关系、键长与键角参考值。常见做法是先用分子编辑工具或脚本生成拓扑文件（如 LAMMPS 的 data 文件、GROMACS 的 .top/.gro），再通过几何算法生成初始坐标。

单链坐标生成最直接的方法是“自回避随机行走”（Self-Avoiding Walk, SAW）。以键长 b=1.54 Å 为步长，从原点出发，每一步在球面上随机取方向，同时检查新原子与已有原子的距离是否小于排斥阈值（通常取 0.8b～0.9b）。若发生冲突则重新采样方向。对于 N=100 的链，SAW 可以快速得到无重叠的线团。但 SAW 生成的链局部键角分布较宽，需要后续用键角约束或短程松弛修正。

更工程化的做法是“模板生长法”：先按反式或旁式构象放置前几个原子，固定键长和键角（如 C-C-C 键角约 109.5°），再通过旋转二面角逐步延伸。这样得到的初始链在局部几何上已经接近真实聚合物，能量最小化时收敛更快。

二、无规线团生成与熔体填充

单链生成后，需要把多条链放入模拟盒子，形成无规线团或熔体。这里有两种典型策略：

• 稀溶液线团：盒子尺寸远大于链的回转半径 Rg，链之间几乎不接触。适合研究单链性质，但需要大量溶剂分子，计算成本高。
• 熔体填充：盒子尺寸由目标密度决定。以 PE 为例，实验熔体密度约 0.85 g/cm³。若体系含 20 条 N=100 的链，总碳原子数 2000，氢原子数 4040，总质量可由原子量求和得到。根据 ρ = m/V，可反推盒子边长 L。例如总质量约 2.8×10⁻²⁰ g，则 V ≈ 3.3×10⁻²⁰ cm³，L ≈ 32 Å。这个尺寸与链的 Rg 相当，链之间会相互穿插，形成真正的熔体。

熔体填充的关键是“先放链、再压缩”或“先放链、再退火”。直接按目标密度随机放置多条链，几乎必然产生原子重叠。工程上常用“分步填充法”：先在较大盒子中放入所有链，使体系密度仅为目标值的 0.3～0.5 倍，然后通过逐步缩小盒子尺寸（或施加外部压力）压缩到目标密度。每一步压缩后都进行能量最小化和短程弛豫，避免局部应力积累。

三、周期性边界条件设置

周期性边界条件（PBC）是消除表面效应的标准手段。设置时需注意三点：第一，盒子必须是长方体或三斜晶胞，且三个方向都满足最小镜像约定；第二，原子坐标应折回主盒子（wrap），但成键原子对不能因为折回而断裂；第三，非键相互作用的截断半径 r_cut 必须小于盒子最短边长的一半，即 r_cut < L_min/2。例如 L=32 Å 时，r_cut 不应超过 16 Å，实际常取 10～12 Å。

在 LAMMPS 中，PBC 通过 boundary p p p 声明，坐标折回用 pbc wrap 或 reset_atoms image 处理。GROMACS 中则通过 pbc = xyz 和 periodic_molecules = yes 控制。需要特别注意的是，如果链本身跨越边界，必须保证整条链的成键拓扑在逻辑上连续，否则分析 Rg、末端距时会出现错误。

四、避免原子重叠与链穿越的实操方法

原子重叠和链穿越是初始构型生成中最常见的两类致命缺陷。原子重叠表现为两个非键原子距离远小于范德华半径之和，导致能量最小化时受力发散；链穿越则表现为两条链的非键原子间距小于 0.9 Å，形成物理上不可能的“穿透”，后续模拟中会表现为虚假的拓扑约束或异常慢的弛豫。

工程上推荐以下组合策略：

1. 低密度预填充：先在 0.3～0.5 倍目标密度下放置所有链，给每条链留出足够的空间。此时即使有局部重叠，也可以通过简单的软核势或逐步压缩消除。
2. 软核势预弛豫：在正式力场之前，先用软核势（soft-core potential）或纯排斥势进行 1000～5000 步的能量最小化，让重叠原子对缓慢推开，而不是瞬间弹开。
3. 逐步压缩：每压缩 5% 体积，进行一次能量最小化和 1～2 ps 的 NVT 弛豫，观察势能是否单调下降。若势能突增，说明压缩过快，应回退并减小压缩步长。
4. 链穿越检测：在压缩完成后，用脚本遍历所有非键原子对，统计最小距离。若存在 d < 0.9 Å 的原子对，且它们不属于同一链的相邻原子，则判定为链穿越，需要重新生成或局部退火。
5. 退火消除残余应力：将体系从 300 K 升温至 600 K，再缓慢降温至 300 K，循环 2～3 次。高温段允许链局部重排，低温段恢复局部有序，可有效消除残余重叠和链穿越。

五、结果判读与常见问题

初始构型是否合格，可以通过三个指标快速判断：第一，最小非键原子间距是否大于 2.5 Å；第二，体系势能是否在能量最小化后收敛到合理范围（PE 熔体约 -4～-6 kcal/mol 每原子）；第三，Rg 是否与理论标度关系 Rg ∝ N^ν（ν≈0.5 对于理想链，ν≈0.6 对于良溶剂）大致吻合。若势能异常高或 Rg 明显偏小，通常意味着存在残余重叠或链穿越，需要回到压缩或退火步骤重新处理。

常见问题包括：压缩过快导致局部密度过高，表现为势能曲线出现尖峰；PBC 设置不当导致原子在边界处“断裂”，表现为成键原子距离突然增大；链穿越未被检测到，表现为模拟初期温度失控或压力异常。解决这些问题的核心思路是“慢压缩、勤检查、多退火”，让体系在每一步都有足够的时间弛豫到局部平衡。

完成初始构型后，建议保存一份“干净”的 data 文件，并记录盒子尺寸、原子数、密度、最小非键距离和 Rg 等关键参数。这份构型将作为后续平衡模拟的起点，其质量直接决定整个模拟项目的成败。