计算化学公社

标题: martini粗粒化聚合物模型在lammps中使用fix bond/create指令聚合成键总是失败怎么办？ [打印本页]

作者
Author: anananante 时间: 2025-6-25 23:00
标题: martini粗粒化聚合物模型在lammps中使用fix bond/create指令聚合成键总是失败怎么办？
本帖最后由 anananante 于 2025-6-25 23:08 编辑

求助各位大佬，第一次做粗粒化，选用了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

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

作者
Author: 那年冬天风在吹 时间: 2025-12-30 11:06

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

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

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

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

作者
Author: yuzc 时间: 2026-1-4 15:05

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

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

欢迎光临计算化学公社 (http://bbs.keinsci.com/)