gromacs-4.0.3 笔记(2)

二 第一次模拟      利用GROMACS自带的例子进行第一次分子动力学模拟实验。步骤如下:建立project目录,将需要的数据文件复制到此文件夹下:

$ mkdir -p ~/project
$ cp ~/gromacs/gromacs-4.0.3/share/gromacs/tutor/speptide/*.pdb ~/project

$ cd ~/project

这是一个S-多肽的PDB结构文件,大家可以用任何一款支持PDB格式的分子可视化软件(如MoldenRasmolJmol等)进行查看,如图所示:

对此分子进行溶液相中的分子动力学模拟大致分为以下步骤:

  • 将PDB文件转换成GROMACS的结构文件(*.gro)和拓扑文件(*.top
  • 在水箱中对分子进行溶剂化
  • 在水箱中对分子进行能量最小化
  • 如果需要的话,向溶剂中添加离子(此处忽略)
  • 在有束缚的情况下,对分子进行一次短暂的分子动力学模拟
  • 在没有束缚的情况下进行完全的分子动力学模拟
  • 分析结果

下面是具体进行的操作:

1. 利用pdb2gmx程序将PDB文件转化为GROMACS格式。

执行以下命令:

$ pdb2gmx -f speptide.pdb -p speptide.top -o speptide.gro

运行后会要求用户选择合适的分子力场。这一步将会由原始的PDB文件得到结构文件speptide.gro以及相应的拓扑文件speptide.top,二者都是ASCII文件,可以使用文本编辑器进行查看。此外,还有一个posre.itp文件,被称做”包含拓扑文件”。对于pdb2gmx程序的其他用法请参阅GROMAS的在线参考手册或运行:

$ pdb2gmx -h

2. 在水箱中对分子进行溶剂化。
这一步通过使用editconfgenbox两个程序实现。前者用来生成一个包裹整个分子的长方体空箱,后者则会读取结构文件并向其中加入水。执行以下命令:

$ editconf -f speptide -o out.gro -d 0.5
$ genbox -cp out.gro -cs -p speptide.top -o b4em.gro

genbox还会修改相应的.top文件,执行tail *.top即可看到genbox加入的溶剂分子的量。

3. 溶剂相的能量最小化。

这一步骤的目的在于消除加入溶剂后肽分子内部的张力,并调整间距小于反的话半径的过于接近的分子。
这一步骤将会用到mdrun程序。在运行之前,需要对已有的.gro文件和.top文件进行预处理,并且编写mdrun需要的参数文件energy_minimization.mdp。
首先来编写参数文件energy_minimization.mdp,其内容如下:

cpp                 =  /usr/bin/cpp
define              =  -DFLEX_SPC
constraints         =  none
integrator          =  steep
nsteps              =  100
emtol               =  2000emstep              =  0.1
nstcomm             =  1
ns_type             =  grid
rlist               =  1
rcoulomb            =  1.0
rvdw                =  1.0
Tcoupl              =  no
Pcoupl              =  nogen_vel             =  no

接下来进行.gro文件,.top文件和.mdp文件的预处理。按照mdrun程序的要求,必须首先将.gro文件,.top文件以及.mdp文件预处理成为.tpr文件(run input file),然后输入mdrun进行处理。我们用到的预处理程序名为grompp。运行以下命令:

$ grompp -v -f energy_minimization.mdp -c b4em.gro -p speptide.top -o energy_minimization.tpr

运行完成后即可得到mdrun的输入文件energy_minimization.tpr。接下来我们可以通过mdrun来进行能量最小化了:

$  mdrun -v -s energy_minimization.tpr -o after_energy_minimization.trr -c after_energy_minimization.gro

得到三个文件:after_energy_minimization.gro,after_energy_minimization.trr和ener.edr,第一个是优化后得到的结构文件,第二个是优化过程的轨迹文件,最后一个则是系统能量变化。

4. 束缚的分子动力学优化。
这一步骤的目的在于在限制原子运动空间范围的情况下,对体系进行短暂的优化。需要对相应的.top文件进行修改以加入束缚条件。在speptide.top中有这样一段:

#ifdef POSRES
#include "posres.itp"
#endif

事实上这是一个条件包含语句,仅仅当变量POSRES的值被设置时,才会包含posres.itp中的内容。这样我们就不需要修改.top文件,使用相同的.top文件也可进行自由状况下的分子动力学模拟。因此,我们只要在相应的.mdp文件中对POSRES变量的值进行设置即可实现束缚状态下的动力学优化。建立参数文件position_restraints.mdp,内容如下:

title               =  Yo
cpp                 =  /usr/bin/cpp
define              =  -DPOSRES
constraints         =  all-bonds
integrator          =  md
dt                  =  0.002    ; ps !
nsteps              =  5000     ; total 10 ps.
nstcomm             =  1
nstxout             =  50
nstvout             =  1000
nstfout             =  0
nstlog              =  10
nstenergy           =  10
nstlist             =  10
ns_type             =  grid
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
Tcoupl              =  berendsen
tc-grps             =  Protein  SOL
tau_t               =  0.1      0.1
ref_t               =  300      300
energygrps          =  Protein  SOL
Pcoupl              =  no
tau_p               =  0.5
compressibility     =  4.5e-5
ref_p               =  1.0
gen_vel             =  yes
gen_temp            =  300.0
gen_seed            =  173529
      然后执行:           $ grompp -f position_restraints.mdp -c after_energy_minimization.gro -r after_energy_minimization.gro -p speptide.top -o position_restraints.tpr
$ mdrun -v -s pr -e position_restraints.tpr -o after_position_restraints.trr -c after_position_restraints.gro >& pr.job &通过$ tail -f pr.job查看运行进度。

5. 无束缚条件下的完全动力学模拟
首先建立对应的参数文件full.mdp,内容如下:

title               =  Yo
cpp                 =  /usr/bin/cpp
constraints         =  all-bonds
integrator          =  md
dt                  =  0.002    ; ps !
nsteps              =  5000     ; total 10 ps.
nstcomm             =  1
nstxout             =  250
nstvout             =  1000
nstfout             =  0
nstlog              =  100
nstenergy           =  100
nstlist             =  10
ns_type             =  grid
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
Tcoupl              =  berendsen
tc-grps             =  Protein  SOL
tau_t               =  0.1      0.1
ref_t               =  300      300
energygrps          =  Protein  SOL
Pcoupl              =  berendsen
Pcoupltype          = isotropic
tau_p               =  0.5
compressibility     =  4.5e-5
ref_p               =  1.0
gen_vel             =  no
gen_temp            =  300.0
gen_seed            =  173529

然后从经过优化的结构文件生成对应的.tpr文件:

$ grompp -v -f full.mdp -o full.tpr -c after_position_restraints.gro -p speptide.top

最后进行完全的分子动力学模拟:

$ mdrun -v -s full.tpr -e full.edr -c after_full.gro -o full.trr >& full.job &

6. 结果分析
在有X-server的机器上运行:

$ ngmx -s position_restraints.tpr -f full.trr

即可查看优化过程中的结构变化,如图所示。

参考文献:http://www.gromacs.org/documentation/reference/online/speptide.html

y2pOO0_Hc2wel9L66mt5gtwdZRebMyqPffLXK7qpaQnBAdEgQISl5TLI4_CDgCLnrN_PvDH7dt3g4Kp0Kdb1DphxQ y2pm-5LTL5_O85qPhnOqFMWk3LkZMLMWfotOi9Z1ti36ZPu-5HJFZidUm4vgli8T9zkaPAvvjvZ89l5QRzTKct4DA

Author: armadillo

傻傻的笨蛋,什么都不懂的Small Kids,总是在幻想,轻轻地走来,静静地站在那里,默默地看着一切,细细地思考,然后悄悄地离开……永远都不愿意留在这里……You mustn't allow yourself to be chained to fate, to be ruled by your genes. Human beings can choose the kind of life that they want to live. What's important is that you choose life... and then live.

2 thoughts on “gromacs-4.0.3 笔记(2)”

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.