MDAnalysis聚合物链构象与持久长度

分析聚合物时，原子编号顺序未必就是沿链顺序。先从真实键拓扑得到链，再进行周期性边界处理，最后才计算键方向相关及拟合。本文使用已有LAMMPS珠簧模拟，不把短链示例的拟合值当作材料参数。

1 下载和运行#

下载水、泛素与珠簧聚合物的实际轨迹、分析脚本和完整数据，保存为mdanalysis-water-polymer-wave330-v2.zip：

module load pymatgen/2026.9.24
python -m zipfile -e mdanalysis-water-polymer-wave330-v2.zip .
cd mdanalysis-water-polymer
bash run.sh
source env.sh

run.sh只对已完成的真实轨迹做后处理，生成results/和六张plots/图片，不启动新的分子动力学。请在新的解压目录运行，避免覆盖自己的输出。MDAnalysis为2.9.0，环境通过现有qvasp 2.26的Python选择，没有调用qvasp或VASP；其他机器可通过ANALYSIS_PYTHON指定兼容环境运行run.sh。

水体系来自已有GROMACS水模拟，884个TIP3P水、2652原子、1001帧，时间0–100 ps。水分析每20帧取一帧，共51帧。轨迹为NPT，实际盒长约29.87494–30.29525 Å，不能当作固定体积轨迹交给只支持固定正交盒的LinearDensity。本例使用能明确处理所需空间网格或平均体积的工具；短轨迹统计仍要检查采样收敛。

2 本例的单位和数据#

聚合物来自真实珠簧聚合物示例：16条链，每条12个珠子，共192粒子、176条FENE键，采用WCA排斥和LAMMPS units lj。101帧覆盖0–40τ，每200个积分步保存一次，原始积分步长为0.002τ。

因此这里长度是σ、时间是τ，不是Å和ps。MDAnalysis默认单位标签不能替代实际单位换算。LAMMPSDUMP读取器用原始积分步数乘dt生成时间，所以这里传入dt=0.002；不能把保存间隔0.4直接作为dt，否则时间轴会放大200倍。

3 按真实连接关系排序并展开每条链#

python - <<'PYCODE'
import MDAnalysis as mda
from MDAnalysis.analysis.polymer import sort_backbone
from MDAnalysis.transformations import unwrap
u = mda.Universe('inputs/polymer/polymer.data',
    'inputs/polymer/trajectory.lammpstrj', topology_format='DATA',
    format='LAMMPSDUMP', atom_style='id resid type x y z', dt=0.002)
u.trajectory.add_transformations(unwrap(u.atoms))
chains = [sort_backbone(g) for g in u.atoms.fragments]
print(len(chains), [len(c) for c in chains])
u.trajectory[-1]
print('Final reduced time / tau:', u.trajectory.ts.time)
PYCODE

得到16条、每条12粒子的连通链，最后时间为40τ。sort_backbone依赖键图，适合线性主链；交联、支化或环状分子不能只按粒子ID排序后直接套用。自己的DATA格式也可能不是本例id resid type x y z，须对应实际atom_style。

unwrap按每条链的键连接处理跨盒问题，不把不同分子拼成一条链。图中每条链保持完整；不同分子的整体位置没有被当成持久长度信息。

4 PersistenceLength：先看相关，再看拟合#

python - <<'PYCODE'
import MDAnalysis as mda
from MDAnalysis.analysis.polymer import sort_backbone, PersistenceLength
from MDAnalysis.transformations import unwrap
u = mda.Universe('inputs/polymer/polymer.data',
    'inputs/polymer/trajectory.lammpstrj', topology_format='DATA',
    format='LAMMPSDUMP', atom_style='id resid type x y z', dt=0.002)
u.trajectory.add_transformations(unwrap(u.atoms))
chains = [sort_backbone(g) for g in u.atoms.fragments]
a = PersistenceLength(chains).run()
print(a.results.bond_autocorrelation)
print('native bond scale / sigma:', a.results.lb)
print('native fitted lp / sigma:', a.results.lp)
PYCODE

工具先统计间隔n根键的单位键向量点积，再对相关曲线拟合exp(-n*lb/lp)。官方聚合物手册给出了定义。原生程序本例输出键长尺度约 0.9700000σ、拟合lp约 1.2694003σ。

需要特别理解这个版本的尺度取法：分析迭代结束后轨迹回到初始帧，lb在此时读取，因此本例拟合横轴采用初始帧的平均键长，不是101帧全部键长的平均。完整逐帧键长在配套图中展示，不能把这个数值的含义自行改成时间平均。

这组原生拟合的相对L2残差约 24%，长间隔相关尾部明显没有被单一指数描述好。链只有12个珠子，轨迹又从有序初态出发，因此这个lp只用于展示工具输出，不作为已平衡材料的持久长度。实际研究应增加链长与平衡/采样，查看原始相关尾部、拟合范围和误差，再决定指数模型是否适用。

5 保存自己需要的曲线#

python - <<'PYCODE'
import numpy as np
z = np.load('results/FULL-FIELDS.npz')
np.savetxt('my-chain-correlation.csv', np.column_stack([
    z['polymer_fit_x_sigma'],z['polymer_bond_C'],z['polymer_fit']]),
    delimiter=',', header='native_scale_sigma,C_bond,native_fit', comments='')
print(z['polymer_positions_sigma'].shape)
PYCODE

完整坐标为(101,16,12,3)。资料包同时保留键相关、拟合值和实际构象，便于检查输入顺序、跨盒处理和拟合适用性。若要换成真实化学聚合物，应使用相应力场和主链原子选择，并统一真实单位；珠簧模型中的σ不能直接标成Å。