LAMMPS粗粒化建模：把几万个原子缩减到几百个珠子，精度不是白送的

粗粒化（CG）的本质是一笔交易——用化学细节换时空尺度。把苯乙烯单体的16个原子映射成2个CG珠子，体系的总粒子数从12万降到1.5万，时间步长从1fs可以拉到10fs，模拟能覆盖的时空范围瞬间扩张了一个数量级。但代价呢？

代价必须一笔一笔地算。一个做聚合物自组装的项目找到我们的时候，全原子模型已经跑不动了——40条聚苯乙烯链（每条100个单体）在显式溶剂里的全原子模拟，跑100ns需要两周。而他们要看到微相分离的完整动力学，至少需要微秒量级。

粗粒化是唯一出路。

映射方案：珠子不是越少越好

第一次尝试用了”一个单体=一个珠子”的映射——苯乙烯的苯环和主链碳被揉成一个CG珠子，质量取104amu。珠子少意味着快，但跑出来的Rg（回转半径）比全原子结果小了25%。

不是随机偏差，是系统性的。苯乙烯单体的苯环侧基在主链上形成了局部刚性——CG珠子是球对称的，丢失了这一各向异性。珠子之间的有效排斥半径必须调大才能逼近真实的链刚度，但调大之后珠子之间的堆积行为又变了。

返工选了双珠映射方案：苯环和主链各成一个CG珠子。珠子的非键参数（LJ势的σ和ε）通过Iterative Boltzmann Inversion（IBI）拟合，目标分布来自全原子模拟的径向分布函数g(r)。第一次IBI迭代：初始LJ σ=4.7Å，epsilon=3.8kJ/mol，10轮迭代后Rg从全原子的2.15nm变到CG的2.08nm，偏差收窄到3%。

双珠方案比单珠慢了约30%，但Rg和g(r)都回到可接受范围了。这里的一条经验：芳环侧基不能映射成一个各向同性珠子。如果CG珠子需要同时承担主链柔性和侧基排斥，物理上不可能——珠子没有”方向”这个自由度。

键合项拟合：分布才是真标准

映射定下来之后，键长和键角的势函数参数怎么取？不是用平衡值——键长取1.53Å、键角取109.5°就可以了。CG模型的键合项要拟合的是分布，不是均值。

全原子模拟里苯环与主链之间的虚拟键长分布不是高斯型——在4.2Å处有个主峰，在5.1Å处有个长尾。这是苯环旋转异构造成的双峰结构。如果用单高斯去拟合，宽峰会同时低估近邻距离和高估远程距离。

我们用了多高斯叠加的表格势（tabulated potential），直接把全原子分布的负对数转换成势能面。代码上LAMMPS的`bond_style table`和`angle_style table`直接支持。键长偏差在全原子分布的±0.1Å范围内还原了98%的构型。

做完键合项拟合，CG模型的键长分布、键角分布与全原子参考的Kullback-Leibler散度降到了0.05以下——意味着构型统计上两个模型几乎等效。

动力学标定：加速因子不是恒定值

一个很多人忽略的问题是：CG模型少了摩擦，动力学比全原子快。快多少、快得是否均匀，必须标定。

这个项目里，CG模型的链中心均方位移（MSD）在100ns之后的斜率是全原子模型的4.2倍。但这不是一个恒定的加速因子——在50ns之前加速比约3.5，50~200ns加速比约4.8。原因是短时间尺度的局部松弛被CG势的软性放大了，而长时间的扩散又受拓扑缠结影响（缠结在CG里部分丢失）。

我们最终采用了一个分段的扩散系数标定：在数据分析和报告里，CG时间轴被除以4.2（取100-200ns的平均值）映射回等效全原子时间。不完美——缠结效应无法完全还原——但在当前CG分辨率下，这是能做到的最好方案。

一个不回避的事实：粗粒化的精度上限在映射方案定下来的那一刻就部分锁死了。单珠换来的速度、双珠换来的精度——二者不可兼得。选哪个，取决于你的科学问题到底要什么。