LAMMPS粗粒化建模 — 从全原子映射到粗粒化力场拟合的实战流程

LAMMPS粗粒化建模是将复杂全原子体系简化为粗粒化模型的过程。与直接使用MARTINI等现成力场不同，LAMMPS粗粒化建模通常需要针对特定体系自定义珠子映射和力场参数，以平衡精度和效率。粗粒化建模的核心挑战是：如何用更少的自由度尽可能保留原体系的关键物理化学性质。

一、粗粒化建模的基本流程

LAMMPS粗粒化建模的标准流程包含四个步骤：

第一步：珠子映射设计。确定每个粗粒化珠子代表哪些原子。映射方案取决于体系特征——蛋白质常用Cα珠子（每个残基1个骨架珠子+1-2个侧链珠子），磷脂用4:1映射（4个重原子为1个珠子），聚合物可以每2-4个重复单元为1个珠子。

第二步：非键参数拟合。确定珠子间的Lennard-Jones参数（ε和σ），使粗粒化模型再现全原子模拟的径向分布函数、密度、表面张力等宏观性质。

第三步：键参数拟合。确定珠子间的键长、键角、二面角参数，使粗粒化模型的构象分布与全原子一致。

第四步：验证。将粗粒化模型用于计算已知的实验或全原子参考数据，检查关键性质是否再现。

二、珠子映射策略

映射方案决定了粗粒化模型的精度上限。一个好的映射应该满足三个原则：保留构型自由度的关键信息、保留化学环境的区分度、计算量显著降低。

蛋白质映射：MARTINI方案是每残基1个骨架珠子+1-2个侧链珠子。更简化的方案是Cα-only模型——只用Cα珠子，侧链信息编码到Cα-Cα的相互作用中。Cα-only模型将100残基蛋白质从约1500个原子降到100个珠子，计算量减少93%，但丢失了侧链细节。

聚合物映射：聚乙烯的粗粒化常采用2:1映射——每2个CH2单元为1个珠子。聚合度1000的PE从2000个原子降到500个珠子。更激进的4:1映射可降到250个珠子，但玻璃化转变温度的精度会下降约20 K。

核苷酸映射：DNA的粗粒化常用3Sites-per-Nucleotide（3SPN）模型——每个核苷酸用3个珠子（糖、磷酸、碱基）。双链DNA从约600原子/碱基对降到6珠子/碱基对，计算量减少99%。

三、非键参数拟合方法

非键参数是粗粒化模型最关键也最难拟合的部分。常用方法有两种：

迭代Boltzmann反演（IBI）：从全原子参考RDF出发，反推有效的粗粒化势能函数。公式：U_{n+1}(r) = U_n(r) + kT ln[g_n(r)/g_ref(r)]。迭代直到粗粒化RDF与全原子参考一致。IBI的缺点是得到的势能是数值表格（tabulated potential），不能用于不同温度或组成的体系。

力匹配（Force Matching, FM）：从全原子轨迹中提取每个珠子受到的合力，拟合一个解析形式的力场（如LJ势）使其重现这些力。FM得到的参数更有转移性，但需要大量参考数据。

在LAMMPS中使用tabulated势能：

pair_style table linear 1000
pair_coeff 1 1 cg_potential.table ENTRY1
pair_coeff 1 2 cg_potential.table ENTRY2

其中cg_potential.table是IBI拟合的势能表格文件，包含r和U(r)两列。

四、键参数校准

键参数的拟合相对简单——从全原子轨迹中提取珠子距离/角度分布，拟合到谐调势或高斯函数上。

键长分布：全原子轨迹中两个珠子（珠子质心）间距离的分布通常近似高斯分布。拟合得到平均距离r0和标准差σ，对应谐调势参数k = kT/σ²。

以PE的2:1映射为例：相邻珠子（代表2个CH2）的距离分布在2.53 Å处峰值，σ=0.08 Å。对应的谐调势参数k = kT/σ² = 0.596/(0.08²) = 93.1 kcal/(mol·Å²)。r0 = 2.53 Å。

键角分布：三连珠子的角度分布通常不是简单高斯，可能有多个峰（对应trans和gauche构象）。如果只用谐调势，无法重现多峰分布。解决方法：用tabulated角度势，或用多项式势。

二面角：粗粒化模型的二面角参数通常不重要——因为珠子代表的原子团已经消除了大部分二面角自由度。对于2:1映射的PE，粗粒化珠子间没有明确的二面角，可以忽略。

五、实战复盘：聚苯乙烯粗粒化模型

项目目标：建立PS的粗粒化模型，用于模拟聚合度500的PS链在薄膜中的自组装行为。

映射方案：每个苯乙烯单元用2个珠子——1个骨架珠子（代表CH-CH₂）+1个苯环珠子（代表C₆H₅）。聚合度500的PS从约10000个原子降到1000个珠子。

第一步：全原子参考模拟。聚合度20的PS链（小链做参考），OPLS-AA力场，NPT 500 K，50 ns。提取珠子质心坐标，计算RDF和键分布。

第二步：非键参数拟合。用IBI方法迭代拟合骨架-骨架、骨架-苯环、苯环-苯环三组非键势能。迭代20轮后，粗粒化RDF与全原子参考的RMS偏差<0.03（收敛标准<0.05）。

第三步：键参数拟合。骨架-骨架键长分布峰值2.52 Å，σ=0.07 Å。骨架-苯环距离分布峰值2.48 Å，σ=0.09 Å。键角分布有两个峰（约65°和115°），用tabulated角度势。

第四步：LAMMPS模型构建。

pair_style hybrid table linear 1000
bond_style harmonic
angle_style table 500

# 非键
pair_coeff 1 1 cg_bb.table BB
pair_coeff 1 2 cg_bs.table BS
pair_coeff 2 2 cg_ss.table SS

# 键
bond_coeff 1 93.1 2.52

# 角度
angle_coeff 1 cg_angle.table ANGLE_BB

第五步：验证。用粗粒化模型跑NPT 500 K，计算密度。结果0.951 g/cm³，全原子参考0.948 g/cm³，偏差0.3%。Tg计算（降温模拟）：粗粒化模型测得Tg=358 K，全原子参考Tg=365 K，偏差7 K（2%）。考虑到粗粒化模型的计算速度是全原子的45倍，这个精度完全可接受。

第六步：大规模应用。用粗粒化模型模拟聚合度500的PS薄膜（500条链，50万珠子），在GPU上跑1 μs。观察到薄膜表面出现垂直于衬底的链取向偏好，与实验GIWAXS结果一致。

六、常见问题与解决方法

IBI参数不可转移：在500 K拟合的势能表用于400 K或600 K时可能不准确——因为粗粒化自由能包含熵贡献，温度敏感。解决方法：在目标温度附近拟合，或用多温度联合拟合。

粗粒化模型过度聚集：自定义粗粒化力场中珠子间吸引力容易偏强，导致不合理的聚集。检查方法：计算体系的第二维里系数B₂，如果B₂过负说明吸引力过强。解决方法：增大排斥项参数或减小吸引项ε。

tabulated势能的数值噪声：IBI拟合的势能表格在端点附近可能有数值噪声，导致能量不守恒。解决方法：对势能表做平滑处理（如Savitzky-Golay滤波），并在LAMMPS中设置足够大的table查找点数（≥1000）。

七、专业计算服务推荐

需要LAMMPS粗粒化建模服务？

科研学术网提供专业的LAMMPS粗粒化建模服务，涵盖珠子映射设计、IBI/FM参数拟合、力场验证、大规模模拟全流程。

• ✅ 博士级工程师团队，一对一技术支持
• ✅ 计算结果可靠，可提供详细的技术报告
• ✅ 周期灵活，加急项目最快3天交付
• ✅ 价格透明，无隐形费用

立即咨询报价 →