计算化学公社

标题: 求助模拟有机分子聚集过程参数如何设置 [打印本页]

作者
Author: tiandikuoyuan 时间: 2021-6-15 15:18
标题: 求助模拟有机分子聚集过程参数如何设置
各位老师，我想模拟有机分子在不良溶剂堆积的过程，咨询一下使用gromacs模拟计算的参数如何设置以及模型构建是否合理。有机分子元素有CHN，分子不带电荷；已经通过Multiwfn计算出RESP电荷，并通过acpype生成了拓扑文件；通过packmol生成了5*5*50nm的盒子，插入了300个分子，随后通过gmx solvate在盒子中填充水分子。分别做了能量最小化、NVT平衡、NPT平衡和MD；由于教程是蛋白质模拟，参数设置可能不是很合理；由于刚接触gromacs，很多参数不太清楚，请老师帮忙看下哪里需要改进。

计算过程能不能把盒子当成孤立体系，不考虑周期性；或者不考虑XY周期，只在Z方向延长；除了手册上提到的参数，设置pbc = no, nstlist = 0, ns-type = simple，其余还有哪里需要考虑的吗?

(, 下载次数 Times of downloads: 67)

由于体系中分子不带电荷，coulombtype是否需要设置为PME，能否直接使用cut-off
参考文献Science, 2020, 368, 1472 (10.1126/science.aba0980)

(, 下载次数 Times of downloads: 63)

当前能量最小化使用的参数minim.mdp：

define = -DFLEXIBLE
integrator = steep
emtol = 250.0
nsteps = 5000
nstenergy = 1
energygrps = System
nstlist = 1
ns_type = grid
coulombtype = PME
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
constraints = none
pbc = xyz

复制代码

NVT平衡使用的参数nvt.mdp：

; 预处理选项
define = -DPOSRES ; 告诉GROMACS运行位置限制性模拟
; 运行控制参数
integrator = md
dt = 0.002 ; 时间步长(单位为ps, 我们使用了2 fs). 只用于动力学积分器(如md), 能量最小化时不需要
nsteps = 50000 ; 模拟步数(总模拟时间为nsteps*dt)
; 输出控制参数
nstxout = 500 ; 输出模拟坐标的频率(nstxout=500且dt=0.002, 所以每1 ps输出一次)
nstvout = 500 ; 速度保存频率
nstenergy = 500 ; 能量保存频率
nstlog = 500 ; log文件输出频率
; 近邻列表参数
nstlist = 5
ns_type = grid
pbc = xyz
rlist = 1.0
; 静电和VDW参数
coulombtype = PME ; 长程静电相互作用的计算方法
pme_order = 4 ; 三次插值
fourierspacing = 0.16 ; FFT间隔
rcoulomb = 1.0 ; 计算静电作用的截断值(单位nm)
vdw-type = Cut-off
rvdw = 1.0
; 温度耦合部分非常重要, 必须正确填写.
tcoupl = v-rescale ; 随机重新调整速度
tc-grps = System ; 与控温器耦合的组(模型中的每个原子或残基都用一定的索引组表示)
tau_t = 0.1 ; 温度耦合的时间常数(单位ps). 必须每个tc_grps指定一个, 且顺序对应
ref_t = 300 ; 代表耦合的参考温度(即动力学模拟的温度, 单位K). 每个tc_grp对应一个ref_t
; 色散校正
DispCorr = EnerPres ; 校正VDW截断
; 不使用压力耦合
pcoupl = no ; NVT中不能使用压力耦合
; 初始速度选项
gen_vel = yes ; 根据Maxwell分布随机产生速度
gen_temp = 300 ; 当你改变温度时, 别忘了改变gen_temp变量以生成速度
gen_seed = -1 ; 随机数生成器的种子
; 键约束选项
constraints = all-bonds ; 使用LINCS算法约束所有键
continuation = no ; 第一次运行
constraint_algorithm = lincs ; 约束算法
lincs_iter = 1 ; LINCS精度
lincs_order = 4 ; LINCS阶数, 与精度有关

复制代码

NPT平衡使用的参数npt.mdp：

define = -DPOSRES
integrator = md
dt = 0.002
nsteps = 50000
nstxout = 500
nstvout = 500
nstfout = 500
nstenergy = 500
nstlog = 500
nstlist = 5
ns-type = Grid
pbc = xyz
rlist = 1.0
coulombtype = PME
pme_order = 4
fourierspacing = 0.16
rcoulomb = 1.0
vdw-type = Cut-off
rvdw = 1.0
Tcoupl = V-rescale
tc-grps = System
tau_t = 0.1
ref_t = 300
DispCorr = EnerPres
; 压力耦合
Pcoupl = Parrinello-Rahman ; Parrinello-Rahman控压器.
Pcoupltype = Isotropic ; isotropic 指盒子可以平均地向各个方向(x, y,z)膨胀或压缩以维持一定的压力. 进行膜模拟时需要用semiisotropic.
tau_p = 2.0 ; 压力耦合的时间常数(单位ps).
compressibility = 4.5e-5 ; 溶剂的压缩系数(4.5e-5为水在300 K和标准大气压下的压缩系数).
ref_p = 1.0 ; 压力耦合的参考压力(单位bar, 1大气压约为0.983 bar).
refcoord_scaling = com
gen_vel = no ; 不产生速度
constraints = all-bonds
continuation = yes
constraint_algorithm = lincs
lincs_iter = 1
lincs_order = 4

复制代码

MD过程md.mdp：

integrator = md
dt = 0.002
nsteps = 500000 ; 1 ns
nstxout = 500
nstvout = 500
nstfout = 500
nstenergy = 500
nstlog = 500
nstlist = 5
ns-type = Grid
pbc = xyz
rlist = 1.0
coulombtype = PME
pme_order = 4
fourierspacing = 0.16
rcoulomb = 1.0
vdw-type = Cut-off
rvdw = 1.0
Tcoupl = v-rescale
tc-grps = System
tau_t = 0.1
ref_t = 300
DispCorr = EnerPres
Pcoupl = Parrinello-Rahman
Pcoupltype = Isotropic
tau_p = 2.0
compressibility = 4.5e-5
ref_p = 1.0
gen_vel = no
constraints = all-bonds
continuation = yes
constraint_algorithm = lincs
lincs_iter = 1
lincs_order = 4

复制代码

作者
Author: sobereva 时间: 2021-6-15 16:33
描述得不清楚。直接把体系结构截图发出来

直接用NPT就完了，根本不需要先做NVT。平衡之前不要用PR压浴

没事不要自己改fourierspacing、nstlist，直接用默认的

绝对不要用constraints = all-bonds，一些教程严重误人子弟

用什么方式算静电作用和分子带不带净电荷完全没有直接关系。周期性体系若无特殊情况一律用PME

作者
Author: tiandikuoyuan 时间: 2021-6-15 19:57

sobereva 发表于 2021-6-15 16:33
描述得不清楚。直接把体系结构截图发出来

直接用NPT就完了，根本不需要先做NVT。平衡之前不要用PR压浴

非常感谢老师指导！我重新调整参数计算一下。

结构为芳香环+烷基链，只是烷基链长度不一样。

(, 下载次数 Times of downloads: 70)

之前试过直接使用NPT，但是体系崩溃了，可能是由于使用了Parrinello-Rahman压力耦合，我换Berendsen试一下。

(, 下载次数 Times of downloads: 77)

整个计算过程中不需要设置键约束吗？
之前使用周期性模拟，分子不能向中间聚集，与文献差距较大，有必要取消周期性吗？也可能是之前参数设置错误导致结构不合理

(, 下载次数 Times of downloads: 79)

作者
Author: sobereva 时间: 2021-6-16 00:49

tiandikuoyuan 发表于 2021-6-15 19:57
非常感谢老师指导！我重新调整参数计算一下。
结构为芳香环+烷基链，只是烷基链长度不一样。

constraint设hbonds

想往中间聚集显然得去掉周期性

作者
Author: 牧生 时间: 2021-6-16 14:13

sobereva 发表于 2021-6-16 00:49
constraint设hbonds

想往中间聚集显然得去掉周期性

我也被相似的问题困扰了很久了。。

做表面活性剂的自组装，文献中我看到最低有不到100 mmol/L的浓度即可形成棒状胶束，而我要做到接近500 mmol/L才能看到棒状胶束。。

如果我建立一个长方体的盒子，去掉了周期性，那么，分子就会在长方体内形成棒状胶束吗？

去掉周期性，就是取消pbc = xyz这一行就行了吗？教程里面的辛醇和水的分离，也形成了棒状的形状，mdp里面就有pbc = xyz，显示XYZ周期以后，看起来就很好。

作者
Author: tiandikuoyuan 时间: 2021-6-16 14:39

sobereva 发表于 2021-6-16 00:49
constraint设hbonds

想往中间聚集显然得去掉周期性

老师，我现在在原先盒子(5*5*80)的基础上增加了5nm的真空层，先做能量最小化，然后做了NVT位置限制性模拟，看能量和势能是达到了平衡；但是跑MD过程中分子还是跨过了真空层，这是因为分子太长了吗（整个分子长度将近5nm）？ (, 下载次数 Times of downloads: 85)
另外，老师能帮忙看下当前设置的参数有哪些需要改进或者错误的地方吗？十分感谢！
能量最小化使用的参数minim.mdp：

define = -DPOSRES
integrator = steep
emtol = 10.0
nsteps = 5000
nstenergy = 1
energygrps = System
nstlist = 10
ns-type = grid
coulombtype = PME
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
constraints = none
pbc = xyz

复制代码

NVT平衡使用的参数nvt.mdp：

define = -DPOSRES
integrator = md
dt = 0.002
nsteps = 50000
nstxout = 500
nstvout = 500
nstenergy = 500
nstlog = 500
nstlist = 10
ns-type = grid
pbc = xyz
rlist = 1.0
coulombtype = PME
pme-order = 4
fourierspacing = 0.12
rcoulomb = 1.0
vdw-type = Cut-off
rvdw = 1.0
tcoupl = v-rescale
tc-grps = System
tau-t = 0.1
ref-t = 300
DispCorr = EnerPres
pcoupl = no
gen-vel = yes
gen-temp = 300
gen-seed = -1
constraints = h-bonds
continuation = no
constraint-algorithm = lincs
lincs-iter = 1
lincs-order = 4

复制代码

MD过程md.mdp：

integrator = md
dt = 0.002
nsteps = 500000 ; 1 ns
nstxout = 5000
nstvout = 5000
nstfout = 5000
nstenergy = 5000
nstlog = 5000
nstlist = 10
ns-type = Grid
pbc = xyz
rlist = 1.0
coulombtype = PME
pme-order = 4
fourierspacing = 0.12
rcoulomb = 1.0
vdw-type = Cut-off
rvdw = 1.0
Tcoupl = v-rescale
tc-grps = System
tau-t = 2.0
ref-t = 300
DispCorr = EnerPres
Pcoupl = no
gen-vel = no
constraints = h-bonds
continuation = yes
constraint-algorithm = lincs
lincs-iter = 1
lincs-order = 4

复制代码

作者
Author: sobereva 时间: 2021-6-17 11:44

tiandikuoyuan 发表于 2021-6-16 14:39
老师，我现在在原先盒子(5*5*80)的基础上增加了5nm的真空层，先做能量最小化，然后做了NVT位置限制性模拟 ...

有真空区的体系明显不应当用DispCorr = EnerPres

都说了要去掉周期性，你还写pbc = xyz干嘛，明显要设no

另外，跑你这种体系，做位置限制的动力学这一步完全是多余的，直接做常规动力学就完了

作者
Author: tiandikuoyuan 时间: 2021-6-17 13:50

sobereva 发表于 2021-6-17 11:44
有真空区的体系明显不应当用DispCorr = EnerPres

都说了要去掉周期性，你还写pbc ...

谢谢老师；我将pbc设为no之后，运行时报错了，说是cutoff-scheme = Verlet不支持pbc = no，但是当前版本cutoff-scheme只有Verlet一个选项，这应该怎么解决？
minim.mdp：

integrator = steep
emtol = 10.0
nsteps = 5000
nstenergy = 1
energygrps = System
nstlist = 10
coulombtype = PME
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
constraints = none
pbc = no

复制代码

错误信息：

ERROR 1 [file minim.mdp]:
With Verlet lists only full pbc or pbc=xy with walls is supported
ERROR 2 [file minim.mdp]:
Can not have Ewald with pbc=no

复制代码

作者
Author: 牧生 时间: 2021-6-18 09:53
我也试了下，cutoff-scheme = Verlet不支持pbc = no

但是cutoff-scheme只有Verlet一个选项才能支持GPU加速。

作者
Author: sobereva 时间: 2021-6-19 00:22

tiandikuoyuan 发表于 2021-6-17 13:50
谢谢老师；我将pbc设为no之后，运行时报错了，说是cutoff-scheme = Verlet不支持pbc = no，但是当前版本[ ...

没PBC了，当然不能用Ewald及其变体
改成cut-off算静电作用
并且改用group方式构建邻居列表代替verlet方式（但版本太新的话不行）

作者
Author: tiandikuoyuan 时间: 2021-6-19 16:53

sobereva 发表于 2021-6-19 00:22
没PBC了，当然不能用Ewald及其变体
改成cut-off算静电作用
并且改用group方式构建邻居列表代替verlet方 ...

好的，老师，我已经安装2019.6版本计算了一下。
请问我想统计MD过程中某一帧中所有分子内某两个原子间的距离分布，比如说计算烷基链末端碳原子和咔唑上氮原子的距离范围和概率，gromacs或者vmd中有现成的命令吗，还是需要自己去编代码计算原子距离？主要是想探讨分子聚集过程中烷基链末端与咔唑之间距离变化范围。

作者
Author: sobereva 时间: 2021-6-20 17:45

tiandikuoyuan 发表于 2021-6-19 16:53
好的，老师，我已经安装2019.6版本计算了一下。
请问我想统计MD过程中某一帧中所有分子内某两个原子间的 ...

得自己写脚本

作者
Author: 小维 时间: 2024-9-13 08:42
本帖最后由小维于 2024-9-13 10:31 编辑

哥，在线吗？关于模拟有机分子聚集过程参数如何设置的一些问题想请教下您。

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