求助 gromacs 水盒子能量最小化问题及一些初级使用问题

lqhenwunai

大家好，我想试算300多个水分子的MD。但是在能量最小化这一步似乎有点问题。

我的拓扑文件如下：
#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/tip3p.itp"

[ System ]
TIP3P in water

[ Molecules ]

参考自https://www.jianshu.com/p/3740fd5094e0

然后通过gmx solvate -cs tip3p -o conf.gro -box 2.2 2.2 2.2 -p topol.top
定义水盒子的大小。

我的mdp文件如下，参考自tutorial第一节里的mdp：
; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep         ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1         ; Frequency to update the neighbor list and long range forces
cutoff-scheme   = Verlet    ; Buffered neighbor searching
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
coulombtype     = PME       ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions

但是我得到的势能曲线确是下图这样的。
1）请问哪里出错了？
2）关于gmx的命令 比如
gmx solvate -cs tip3p -o conf.gro -box 2.2 2.2 2.2 -p topol.top
请问哪里有关于gmx各个module的list？比如，如果我不知道gmx存在solvate这么一个module，那我可能就去用packmol了。哪里能找到各个module及功能介绍的总结？
3）有什么好办法能减少gmx的命令行输入？比如存成一些别名之类的，例如alias gmx-mdrun = 'gmx mdrun $1 $2 ...' 这类的？可能用的还不多，记不住各种options。

wuzhiyi

你这不是能量曲线啊？是virial在YY轴上的分量

sobereva

2 自行把手册里各个子程序的介绍过目一遍。

PS：北京科音分子动力学与GROMACS培训班里会把所有常用子程序的基本使用都结合实际应用例子讲一遍，其它的不常用子程序也都会简要提及，不至于有现成的子程序还以为不存在。

3 这种做法毫无意义。把-h显示的解释看一遍，把常用的、对自己有用的在记事本里记下来就完了

lqhenwunai

wuzhiyi 发表于 2020-8-13 23:30
你这不是能量曲线啊？是virial在YY轴上的分量

谢谢！

lqhenwunai

wuzhiyi 发表于 2020-8-13 23:30
你这不是能量曲线啊？是virial在YY轴上的分量

您好，不好意思再请教一下，我现在是392个水分子在300K下nvt系综下跑，之后又续跑到50000ps。但是总能量还是显示没有收敛。请问这样正常么？我最后跑完的nvt8.part0007.gro如下图。总能量变化如下图。

以下是我的mdp文件：
title                   = water in a box

; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 100000     ; 2 * 100000 = 200 ps
dt                      = 0.002     ; 2 fs
; Output control
nstxout                 = 500       ; save coordinates every 1.0 ps
nstvout                 = 500       ; save velocities every 1.0 ps
nstenergy               = 500       ; save energies every 1.0 ps
nstlog                  = 500       ; update log file every 1.0 ps
; Bond parameters
continuation            = no        ; first dynamics run
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
DispCorr                = EnerPres  ; account for cut-off vdW scheme
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT

; Pressure coupling is off
pcoupl                  = no        ; no pressure coupling in NVT
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Velocity generation
gen_vel                 = yes       ; assign velocities from Maxwell distribution
gen_temp                = 300       ; temperature for Maxwell distribution
gen_seed                = -1        ; generate a random seed

wuzhiyi

lqhenwunai 发表于 2020-8-14 15:16
您好，不好意思再请教一下，我现在是392个水分子在300K下nvt系综下跑，之后又续跑到50000ps。但是总能量 ...

我觉得挺怪的50ns早该平衡了 不过我一般搞npt 搞不好nvt比较慢