martini粗粒化聚合物模型在lammps中使用fix bond/create指令聚合成键总是失败怎么办？

anananante

求助各位大佬，第一次做粗粒化，选用了martini力场，需要异氰酸酯和醇聚合。
我的键结参数都是martini官方教程中给出的bartender拟合生成的，pair coeffs遵循了martini3的粗粒化规则，目前就是单体盒子事先跑npt和nvt没有问题，一旦开始用fix bond/create指令成键，就会在运行一段时间后报错bond atom missing。我的体系按照数量最多能成2000个键，我调整过各种参数最多跑到成三四百个键左右就会报错。
请问是我的模型构建的有问题还是参数设置不合理？还是说martini有什么特殊的设置要求？下面贴上我的映射和in文件
目前是仅仅在令IPDI和PPG反应，OPPEA准备二阶反应，暂时还没用到，盒子中有这3个单体
# ----------------- Init Section -----------------
units           real
dimension 3
boundary        p p p
neighbor 3 bin
neigh_modify every 1 delay 0 check yes

#-----style-----------
atom_style      full
pair_style      lj/cut  20.0
pair_modify             tail no
pair_modify     mix arithmetic
bond_style      harmonic
angle_style     harmonic
special_bonds lj/coul 0.0 0.0 1.0 angle yes dihedral yes

# ----------------- Atom Definition Section -----------------
read_data stage1.data
# ============= coeffes Section ===================
pair_coeff     1        1                     0.88191    4.7
pair_coeff     2        2                     0.88191    4.7
pair_coeff     3        3                     0.81021    4.7
pair_coeff     4        4                     0.56165    4.1
pair_coeff     5        5                     0.67876    4.1
pair_coeff     6        6                0.6214     4.1
pair_coeff     7        7                     0.88191    4.7
pair_coeff     8        8                     0.36089    3.4
pair_coeff     9        9                     0.88191    4.7
pair_coeff     10        10                     0.67876    4.1

# -- harmonic bondType parameter  --
bond_coeff      1                 7.3439    3.45
bond_coeff      2                 7.3439    3.45
bond_coeff      3                 7.3439    3.45
bond_coeff      4                 7.3439    3.45
bond_coeff      5                 7.3439    3.45
bond_coeff      6                 7.3439    3.45
bond_coeff      7                 7.3439    3.45
bond_coeff      8                 7.3439    3.45
bond_coeff      9                 7.3439    3.45

# -- cosine/squared angleType parameter  --
angle_coeff     1                5.6804  180
angle_coeff     2                5.6804  180
angle_coeff     3                5.6804  180
angle_coeff     4                5.6804  180
angle_coeff     5                5.6804  180
angle_coeff     6                5.6804  180
angle_coeff     7                5.6804  180
angle_coeff     8                5.6804  180

# -- harmonic bondType parameter  --
bond_coeff      10                 62.9010      3.24
bond_coeff      11                 111.7513     2.55
bond_coeff      12                 4.6309       3.08
bond_coeff      13                 85.6537      2.77
bond_coeff      14                 125.9143     2.81

# -- cosine/squared angleType parameter  --
angle_coeff     9                3.2886   125.4
angle_coeff     10                7.9372   115.28
angle_coeff     11                4.1634   119.19
angle_coeff     12                8.3339   109.88


# -- harmonic bondType parameter  --
bond_coeff      15                 6.6422    3.41
bond_coeff      16                 4.5462    3.07
bond_coeff      17                 282.3161  2.08
bond_coeff      18                 294.4805  2.09
bond_coeff      19                 33.2672   2.36
bond_coeff      20                 296.0863  2.07
bond_coeff      21                 269.0198  2.1
bond_coeff      22                 273.3664  2.1
bond_coeff      23                 292.1409  2.07
bond_coeff      24                 7.0162    3.2
bond_coeff      25                 7.0162    3.04
# -- cosine/squared angleType parameter  --
angle_coeff     13                5.5137   125.26
angle_coeff     14                12.8128  103.29
angle_coeff     15                23.0826   92.45
angle_coeff     16                27.9845  145.31
angle_coeff     17                34.0288  106.72
angle_coeff     18                35.1402  163.08
angle_coeff     19                36.0842  161.65
angle_coeff     20                233.5891  104.79
angle_coeff     21                13.7640    180.0
angle_coeff     22                13.7640    180.0

#-----NVT(298K)--------------
#velocity all create 298 161013 rot yes dist gaussian
min_style sd
minimize 1.0e-4 1.0e-6  10000 1000000
timestep 1

fix stage1h all  bond/create 200  1 5  4.2  24  prob 0.05 16101 iparam 1 9  jparam  1 10 #atype 21
dump    mydump all  custom 10000 stage1.lammpstrj id type x y z mol
fix 1 all nvt temp 298 298 100  
restart 10000 tmp.restart

thermo_style custom step elapsed temp press pe ke etotal vol density f_stage1h[1]  f_stage1h[2]
thermo 1000

run  10000000
write_data stage1out.data

yuzc

1.感觉你的参数有问题，martini大概是几千的bond k，换算成kcal/mol也不应该这么小。
2. 如果 fix bond/create 让两个分子突然沾在一起，它们原来不共享 angle ， 现在共享 angle，结果这个 angle 的 eq θ0 和当前几何角度差了几十度，乘上几百的 kθ，系统同样会瞬间爆炸。这也是官方文档里提及的注意的。
3. 我建议你跑粗粒化使用专用的粗粒化软件，从构型生成到其他功能都要更专业一些。比如hoomd-blue以及我们课题组开发的pygamd。

那年冬天风在吹

yuzc 发表于 2025-10-28 17:35
1.感觉你的参数有问题，martini大概是几千的bond k，换算成kcal/mol也不应该这么小。
2. 如果 fix bond/c ...

请问，pygamd能做蛋白的粗粒化模拟吗？

jrfjrf123

不知道答主还在做类似的体系没，我说一下我的理解：
1. 就像二楼的答主回复的那样，部分bond参数是肯定有问题的，根据Alessandri和Marrink 2019年发表在JCTC上的文章"Pitfalls of the Martini Model"里描述的，过低的K_bond是不利于模拟相分离过程的，会产生很大的偏差，至少得在500kJ/mol(也就是120kcal/mol左右)以上才能保证能有一个正确的结果。


2. 可以考虑使用bond/react功能来创建新的键，因为bond/react 有stability 功能，也就是会将设定好的反应中的原子使用nve/limit对其运动进行限制，其余未参加反应的bead会保持npt/nvt系综，能有效的避免创建的键过长导致的局部原子受力过大产生的错误。

yuzc

那年冬天风在吹 发表于 2025-12-30 11:06
请问，pygamd能做蛋白的粗粒化模拟吗？

当然可以，本身的分子动力学模拟引擎没什么不能做的。