LAMMPS粗粒化建模：粗粒化分子动力学方法

粗粒化（Coarse-Grained, CG）分子动力学是将多个原子合并为一个”珠子”（bead）来表示的简化模拟方法。通过降低体系自由度，粗粒化模拟可以在保持关键物理性质的前提下，将模拟的时间和空间尺度扩展数个数量级，从而研究传统全原子模拟无法触及的介观过程。

LAMMPS作为通用分子动力学软件，支持多种粗粒化力场和建模方法，是进行粗粒化模拟的常用工具。

粗粒化基本原理

粗粒化映射

将全原子结构映射到粗粒化珠子：

映射原则：

• 珠子质量集中在质心
• 保留关键拓扑结构
• 保持体系的对称性
• 平衡精度与效率

尺度扩展效应

常用粗粒化力场

1. MARTINI力场

MARTINI是最广泛使用的生物分子粗粒化力场：

特点：

• 4个重原子映射为1个珠子（4:1映射）
• 珠子类型：Q(带电)、P(极性)、N(非极性)、C(疏水)
• 每类再细分为子类（如P1-P5极性递减）
• 水分子4个映射为1个珠子

LAMMPS中的MARTINI实现：

# data文件中粒子类型
4 atom types
# 类型1: P4 (水)
# 类型2: Qa (带负电)
# 类型3: Na (疏水)
# 类型4: Qd (带正电)

# 非键相互作用
pair_style lj/cut 1.2
pair_coeff 1 1 5.0 0.47 1.2   # 水-水
pair_coeff 2 3 4.0 0.47 1.2   # Qa-Na
pair_coeff 2 4 0.0 0.47 1.2   # Qa-Qd（特殊处理）

MARTINI参数表（部分）：

2. 自下而上粗粒化

通过全原子模拟数据推导粗粒化力场参数：

力匹配方法（Force-Matching / IBM）：

1. 运行全原子MD模拟获取轨迹
2. 计算每帧中每个珠子的受力
3. 拟合粗粒化势函数使其重现这些力
4. 迭代优化直到收敛

相对熵最小化：

1. 定义全原子和粗粒化体系的概率分布
2. 最小化两个分布之间的相对熵
3. 获得最优粗粒化势函数

Boltzmann反演：

1. 从全原子轨迹获取径向分布函数g(r)
2. 通过Boltzmann反演获得势能： U(r) = -kT·ln(g(r))
3. 对势能进行迭代修正

3. DPD（Dissipative Particle Dynamics）

DPD是特殊的粗粒化方法，珠子之间有三种力：

F = F_C + F_D + F_R

• F_C：保守力（软排斥）
• F_D：耗散力（阻力）
• F_R：随机力（热浴）

LAMMPS DPD设置：

pair_style dpd 1.0 1.0 12345
pair_coeff 1 1 25.0 4.5 1.0
# 参数：A(排斥), γ(耗散), σ(随机=√(2γkT))

LAMMPS粗粒化建模流程

1. 全原子到粗粒化的结构映射

Python映射脚本示例：

# 将全原子PDB映射为粗粒化data文件
# 以MARTINI 4:1映射为例

residues = {}
for atom in all_atoms:
    res_id = atom['residue_id']
    if res_id not in residues:
        residues[res_id] = []
    residues[res_id].append(atom)

# 计算每个残基的质心作为粗粒化珠子位置
cg_beads = []
for res_id, atoms in residues.items():
    com = calculate_com(atoms)  # 质心
    bead_type = determine_bead_type(atoms)  # MARTINI类型
    cg_beads.append({'pos': com, 'type': bead_type})

2. 拓扑结构生成

键合相互作用：

LAMMPS拓扑设置：

# 键
bond_style harmonic
bond_coeff 1 1250.0 0.35  # k=1250 kJ/mol/nm², r_0=0.35nm

# 键角
angle_style harmonic
angle_coeff 1 520.0 96.0  # k=520 kJ/mol/rad², θ_0=96°

# 二面角
dihedral_style harmonic
dihedral_coeff 1 0.0 1 0  # 根据需要设置

3. 力场参数确定

非键相互作用参数确定流程：

1. 目标性质选择：决定需要重现的实验/全原子数据
   • 密度
   • 径向分布函数
   • 界面张力
   • 分配系数
2. 初始参数估计：从文献或Boltzmann反演获取
3. 迭代优化：

   while 未收敛:
       运行粗粒化MD
       计算目标性质
       与参考数据对比
       调整力场参数
   

4. 验证：检查未用于参数化的性质是否一致

4. LAMMPS输入脚本

典型粗粒化MD脚本：

粗粒化建模的关键问题

1. 可转移性问题

粗粒化力场通常在参数化条件下表现良好，但转移到其他条件时可能失效：

2. 熵效应

粗粒化丢失了原子尺度的自由度，导致熵减少：

补偿方法：

• 在有效势中包含熵贡献
• 使用温度相关势函数
• 增加有效排斥力补偿缺失的构型熵

3. 时间尺度映射

粗粒化的时间不等于真实时间，需要时间尺度映射因子：

τ_CG / τ_AA = λ_T

λ_T通常在1-100之间，取决于粗粒化程度和体系。

确定方法：

• 对比全原子和粗粒化的扩散系数
• 对比特定动力学过程的速率
• 标定时间映射因子

常见粗粒化体系

聚合物熔体

Kremer-Grest模型：

# FENE键 + WCA排斥
bond_style fene
bond_coeff 1 30.0 1.5 1.0 1.0

pair_style lj/cut 1.12
pair_coeff 1 1 1.0 1.0 1.12
pair_modify shift yes

脂质双分子层

蛋白质

常见问题与解决方案

问题1：粗粒化力场重现性差

原因：力场参数未充分优化或体系超出参数化范围。

解决方案：

• 扩大参数化的热力学条件范围
• 使用迭代Boltzmann反演(IBI)优化
• 针对目标性质专门参数化
• 验证多个独立性质

问题2：粗粒化结构不稳定

原因：键合参数设置不当或缺乏约束。

解决方案：

• 增加键弹簧常数
• 使用约束算法（SHAKE/RATTLE）
• 检查初始结构合理性
• 减小时间步长

问题3：相行为不正确

原因：粗粒化力场可能无法准确重现相变。

解决方案：

• 在多个相态下参数化
• 使用状态点依赖的势函数
• 验证临界温度和压力
• 考虑使用多态粗粒化方法

问题4：与全原子结果不一致

• 检查映射方案的合理性
• 验证非键参数是否正确
• 检查静电相互作用处理
• 考虑有效电荷的标定

实操案例：聚合物自组装粗粒化模拟

1. 体系：两嵌段共聚物A₅B₅，1000条链
2. 映射：每个单体映射为1个珠子
3. 力场：Kremer-Grest（FENE + LJ）
4. 参数：
   • A-A: ε=1.0, σ=1.0
   • B-B: ε=1.0, σ=1.0
   • A-B: ε=0.5, σ=1.0（排斥）
5. 设置：
   • NVT, T=1.0
   • dt=0.005
   • 10⁸步
6. 结果：
   • A和B嵌段分离
   • 形成层状结构
   • 周期约10σ
   • 与RPA理论预测一致

总结

LAMMPS粗粒化建模通过将多个原子映射为珠子，大幅扩展了分子动力学模拟的时空尺度，使研究介观过程成为可能。在实际应用中，需要根据目标体系选择合适的粗粒化策略和力场参数，通过系统验证确保粗粒化模型能够准确重现关键物理性质。

我们提供专业的LAMMPS粗粒化分子动力学计算服务，涵盖聚合物、生物膜、蛋白质和自组装体系等各类介观模拟，支持从粗粒化映射到力场参数化的完整建模流程。