求助关于使用CGMD对施加了部分位置限制的蛋白分子进行NPT模拟的mdp设置问题

DoubleSheep

各位gromacs前辈好，本人近期在跑一个病毒衣壳的粗粒度模拟，目的是为了研究病毒表面spike蛋白的稳定性。因为病毒太大所以只选择了部分的衣壳结构进行模拟，并对衣壳底座backbone施加了position restraints以保证其结构稳定，而病毒spike则不做限制。模型如图，黄色是spike，灰色是capsid base




在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），但是这样设置同样会导致体系不稳定。

非常感谢！

Frozen-Penguin

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

DoubleSheep

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

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