计算化学公社

标题: 求助关于使用CGMD对施加了部分位置限制的蛋白分子进行NPT模拟的mdp设置问题 [打印本页]

作者
Author: DoubleSheep 时间: 2024-3-8 13:21
标题: 求助关于使用CGMD对施加了部分位置限制的蛋白分子进行NPT模拟的mdp设置问题
各位gromacs前辈好，本人近期在跑一个病毒衣壳的粗粒度模拟，目的是为了研究病毒表面spike蛋白的稳定性。因为病毒太大所以只选择了部分的衣壳结构进行模拟，并对衣壳底座backbone施加了position restraints以保证其结构稳定，而病毒spike则不做限制。模型如图，黄色是spike，灰色是capsid base

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

在production前期的预平衡阶段的步骤如下：

1. Energy minimization with softcore potentials: 100 steps, 0.02ps time step, flexible water, steep integrator, V-rescale temperature coupling, 310K temperature, no pressure coupling.
2. Energy minimization with positional restraints: 10,000 steps, 0.02ps time step, flexible water, steep integrator, no temperature and pressure coupling.
3. Equilibration with positional restraints: 10,000 steps, 0.001ps time step, rigid water, sd integrator, V-rescale temperature coupling, 310K temperature, no pressure coupling.
4. Solvation, ionization and Equilibration with positional restraints: 100,000 steps, other parameters are identical to 2)
5. NVT equilibration with positional restraints: 500,000 steps, 0.002ps time step, rigid water, sd integrator, V-rescale temperature coupling, 310K temperature, no pressure coupling.
6. Short NPT run with positional restraints: 750,000 steps, 0.02ps time step, rigid water, md integrator, V-rescale temperature coupling, 310K temperature, Berendsen barostat, isotropic pressure.

然后进行production，mdp参数如下：

define = -DPOSRES
integrator = md
dt = 0.02
nsteps = 50000000
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 10000
nstxout-compressed = 10000
compressed-x-precision = 1000
comm-mode = linear
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005
coulombtype = reaction-field
rcoulomb = 1.1
epsilon_r = 15 ; 2.5 (with polarizable water)
epsilon_rf = 0
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1
tcoupl = v-rescale
tc-grps = system
tau-t = 1.0
ref-t = 310
continuation = yes
gen-vel = no
Pcoupl = berendsen
Pcoupltype = isotropic
tau-p = 4.0
compressibility = 3e-4
ref-p = 1.0
refcoord-scaling = com

想请教大家，目前的mdp对于模拟该体系是合理的吗。本人目前主要对【comm-mode】【refcoord-scaling】及【Pcoupl】的参数选择不太确定。pcoupl 选择 berendsen 而非 parrinello-rahman 主要是由于后者会造成体系崩溃。按网上一些资料的说法和gromacs的建议，使用position restraints时应当不移除COM的移动（comm-mode设置为none），但是这样设置同样会导致体系不稳定。

非常感谢！

作者
Author: Frozen-Penguin 时间: 2024-3-8 14:14
comm-mode是用来消除分子整体平动转动的，你的体系用了位置限制，所以不会出现上述情况，这个选项没有影响。
refcoord-scaling用来修改位置约束时参考位置的变化，由于压强平衡会改变盒子大小，位置约束也要随之改变，如果你希望被约束的分子始终位于盒子中央，就应该选择com。
berendsen得到的系综不是严格的系综，所以一般只用于平衡阶段，正式模拟还是要用PR，或者如果版本比较新可以用C-rescale，另外控压算法是否稳定还与tau-p的数值有关，粗粒化模拟的步长通常比全原子大，相应的tau-p和tau-t也要改大。

作者
Author: DoubleSheep 时间: 2024-3-8 14:42

Frozen-Penguin 发表于 2024-3-8 14:14
comm-mode是用来消除分子整体平动转动的，你的体系用了位置限制，所以不会出现上述情况，这个选项没有影响 ...

好的，谢谢老师的建议！我再多试一些参数

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