####################

### 概述 ###

####################

-----------------------------------------------------------------

-----------------------------------------------------------------

该例程来自Gromacs程序/share/tutor/目录下。整个例程大概只需要十分钟就可以完成，非常适合初学者学习。该例程是一个完整的分子动力学模拟过程，涵盖了Gromacs程序基本的使用方法。模拟内容是一个水环境下的小肽链。模拟唯一要求是该小肽链的PDB文件。

相关内容请参阅Gromacs文档，或者给Gromacs开发组询问。

-----------------------------------------------------------------

-----------------------------------------------------------------

#########################

### 环境变量设置 ###

#########################

-----------------------------------------------------------------

-----------------------------------------------------------------

在以下的例程中，所有命令都直接运行，没有添加绝对路径。所以，必须将Gromacs安装路径的bin文件夹加入到系统PATH变量中。如果不加入PATH变量，那么运行时要加入命令的绝对路径。

-----------------------------------------------------------------

-----------------------------------------------------------------

################

### PDB2GMX ###

################

-----------------------------------------------------------------

-----------------------------------------------------------------

在进行如何分子动力学模拟之前，必须建立分子的拓扑文件。Gromacs的分子拓扑文件是用pdb2gmx命令生成，文件后缀名为 .top。攑db2gmx的输入唯一文件是分子的PDB文件，可以从www.pdb.org寻找，文件后缀名 .pdb。”

不是所有的PDB文件都含有氢原子，pdb2gmx将添加所有缺失的氢原子。在pdb2gmx输出的gromacs结构文件中，包含了蛋白质结构的每一个原子，并定义了分子结构的尺寸大小。文件后缀名为 .gro。”

-----------------------------------------------------------------

-----------------------------------------------------------------

假设分子PDB文件名为MOL.pdb，命令运行格式：

pdb2gmx -f MOL.pdb -o MOL.gro -p MOL.top > output.pdb2gmx

命令得到三个文件：MOL.gro是gromacs结构文件；MOL.top分子拓扑文件；output.pdb2gmx是pdb2gmx命令输出文件。

###############

### GENBOX ###

###############

-----------------------------------------------------------------

-----------------------------------------------------------------

真空条件下进行分子模拟输出结果误差较大，所以模拟之前，必须给分子添加水环境。Gromacs中使用genbox命令完成。

Genbox命令读入gromacs的结构文件，并读入水盒子尺寸的大小，输出文件包含分子文件和水盒子。同时genbox更改原来的拓扑文件，使其包含水分子。在使用genbox之前，要使用editconf命令定义水盒子大小。

-----------------------------------------------------------------

-----------------------------------------------------------------

命令格式：

editconf -f MOL.gro -o MOL.gro -d 0.5 > output.genbox

genbox -cp MOL.gro -cs -o MOL_b4em.gro -p MOL.top >> output.genbox

##############

### GROMPP ###

##############

-----------------------------------------------------------------

-----------------------------------------------------------------

原则上讲，到此为止已经可以开始进行模拟计算了。但是，添加在准备分子系统过程中，系统有很多原子距离太近，局部能量太高。这些相互距离太近的原子多是由于genbox程序产生的，溶剂分子与蛋白分子之间存在不稳定的高能量。如果现在开始模拟计算，而不进行能量最优化，系统将可能很不稳定。

去除这些局部高能量的办法是对系统进行能量最优化。能量最优化过程是改变系统中局部高能量的原子的位置，降低这些点的能量。

在进行能量最优化之前，我们先用GROMACS的预处理程序grompp处理所有输入文件。grompp预处理拓扑文件(.top)，坐标文件(.gro)和一个参数文件(.mdp)，然后输出一个二进制拓扑文件(.tpr)。这个二进制文件包含所有模拟需要的信息，利用这个文件即可以进行能量最优化和动力学模拟。

-----------------------------------------------------------------

em.mdp文件范例，详细请参考Gromacs手册

-----------------------------------------------------------------

title = MOL

cpp = /usr/bin/cpp

define = -DFLEX_SPC

constraints = none

integrator = steep

dt = 0.002 ; ps !

nsteps = 100

nstlist = 10

ns_type = grid

rlist = 1.0

rcoulomb = 1.0

rvdw = 1.0

;

; Energy minimizing stuff

;

emtol = 1000.0

emstep = 0.01

---------------------------------------------------------------------------------

Grompp命令格式：

gompp -f em -c MOL_b4em -p MOL -o MOL_em >&! output.grompp_em

################

### MDRUN EM ###

################

-----------------------------------------------------------------

-----------------------------------------------------------------

有了二进制拓扑文件，现在可以开始进行能量最优化了。进行能量最优化的程序是mdrun，在

Gromacs程序中，所有模拟都是用mdrun程序进行。

在进行能量最优化过程中，注意查看mdrun程序的输出文件。在输出文件在中，从左到右第一个数字是模拟步数，第二个数字是计算步长，第三个数字是系统能量。如果模拟顺利利，可以看到系统能量从一个很高的值迅速降低，最后稳定在一个大负值。

-----------------------------------------------------------------

-----------------------------------------------------------------

命令格式：

mdrun -nice 4 -s MOL_em -o MOL_em -c MOL_b4pr -v >& ! output.mdrun_em

#################

### GROMPP PR ###

#################

-----------------------------------------------------------------

能量最优化之后，首先保持蛋白质不动，对蛋白质周围水环境进行动力学模拟，该过程称为位置限制性分子动力学(position restrained MD)

位置限制分子动力学保持蛋白质位置不变，对溶剂分子进行平衡计算，可以使溶剂分子填补空间空洞。这个可能存在的空洞是genbox程序产生的。

首先对输入文件进行预处理得到二进制拓扑文件，这些输入文件就是能量最优化得到的输出文件。然后再写一个参数配置文件和索引文件。

默认情况下，程序对模拟系统是分部分的。我们利用两个部分进行位置限制性模拟：Protein和

SOL部分，分别表示蛋白质和溶剂。这这个过程中，即保持protein位置不变。

参数文件(.mdp)包含了位置限制性模拟的所有参数，包括步长，步数，温度等等。同时参数文件告诉GROMACS模拟的类型，如能量最优化、位置限制性或者是分子动力学模拟。

-----------------------------------------------------------------

pr参数文件

-----------------------------------------------------------------

title = MOL position restraining

cpp = /usr/bin/cpp

define = -DPOSRES

constraints = all-bonds

integrator = md

dt = 0.002 ; ps !

nsteps = 500 ; total 1.0 ps.

nstcomm = 1

nstxout = 10

nstvout = 1000

nstfout = 0

nstlog = 10

nstenergy = 10

nstlist = 10

ns_type = grid

rlist = 1.0

rcoulomb = 1.0

rvdw = 1.0

; Berendsen temperature coupling is on in two groups

Tcoupl = berendsen

tau_t = 0.1 0.1

tc-grps = protein sol

ref_t = 300 300

; Pressure coupling is not on

Pcoupl = no

tau_p = 0.5

compressibility = 4.5e-5

ref_p = 1.0

; Generate velocites is on at 300 K.

gen_vel = yes

gen_temp = 300.0

gen_seed = 173529

----------------------------------------------------------------------------------------

命令格式

grompp -f pr -c MOL_b4pr -r MOL_b4pr -p MOL -o MOL_pr >& ! output.grompp_pr

################

### MDRUN PR ###

################

-----------------------------------------------------------------

-----------------------------------------------------------------

现在进行1ps的位置限制性分子动力学模拟(the Position restrained Molecular Dynamics simulation)在这里进行1ps只是为了节省例程演示实践时间，在实际研究中，1ps模拟将太短，不符合实际情况。

-----------------------------------------------------------------

-----------------------------------------------------------------

命令格式

mdrun -nice 4 -s MOL_pr -o MOL_pr -c MOL_b4md -v >& ! output.mdrun_pr

#################

### GROMPP MD ###

#################

-----------------------------------------------------------------

-----------------------------------------------------------------

到此为止，分子系统已经可以进行正式的分子动力学模拟了。再一次利用grompp程序将输入文件转化成二进制拓扑文件(.tpb/.tpr文件后缀)。

分子动力学模拟配置文件(参数详细意义请参考Gromacs文档)

-----------------------------------------------------------------

title = MOL MD

cpp = /usr/bin/cpp

constraints = all-bonds

integrator = md

dt = 0.002 ; ps !

nsteps = 5000 ; total 5 ps.

nstcomm = 1

nstxout = 50

nstvout = 0

nstfout = 0

nstlist = 10

ns_type = grid

rlist = 1.0

rcoulomb = 1.0

rvdw = 1.0

; Berendsen temperature coupling is on in two groups

Tcoupl = berendsen

tau_t = 0.1 0.1

tc-grps = protein sol

ref_t = 300 300

; Pressure coupling is not on

Pcoupl = no

tau_p = 0.5

compressibility = 4.5e-5

ref_p = 1.0

; Generate velocites is on at 300 K.

gen_vel = yes

gen_temp = 300.0

gen_seed = 173529

----------------------------------------------------------------

命令格式：

grompp -f md -c MOL_b4md -p MOL -o MOL_md >& ! output.grompp_md

################

### MDRUN MD ###

################

-----------------------------------------------------------------

-----------------------------------------------------------------

现在进行分子动力学模拟计算。详细观察模拟步数的增加。(整体模拟步数为2500，时间长度为5ps。再一次强调，模拟时间太短，仅适合教学演示。)

-----------------------------------------------------------------

-----------------------------------------------------------------

mdrun -nice 4 -s MOL_md -o MOL_md -c MOL_after_md -v >& ! output.mdrun_md

############

### NGMX ###

############

-----------------------------------------------------------------

-----------------------------------------------------------------

最后，可以使用Gromacs的ngmx程序观看动力学模拟轨迹(文件后缀.trj)。轨迹文件包含了分子动力学过程的坐标，原子速度和受力情况等等。可供分析使用。Ngmx的详细使用方法，请参考Gromacs文档。由于分子动力学模拟分析方法各异，这里不再讨论。