Gromacs命令行参数解析

Gromacs命令行参数解析

命令参数解析

gmx pdb2gmx

命令行功能
gromacs在处理数据之前要将pdb文件转化为gromacs能够处理的文件格式, 该步骤相当于是一个准备工作,准备之后进行能量最小化等操作所需要的结构文件。

gmx pdb2gmx  -f input.pdb -o input.gro -water spce

• -f 参数用来指定输入的pdb文件（可以是其他格式的结构文件），该pdb文件中的原子类型要与力场中的原子类型要一致，假如存在力场中不存在的原子，那么命令执行会失败。比如我利用Discovery Studio构建的抗体的pdb模型，加了H原子之后，利用gmx pdb2gmx命令就失败了，显示某些氢原子的命名方式与gmx中给定的命名方式不同，但是没有经过DS加氢原子操作的pdb文件就可以利用gmx进行分析。
• -o 指定输出的gro文件的名称，假如没有指定名称那么会有一个默认的输出名称（conf.gro）
• -water 指定水分子的模型，可选模型有none, spc, spce, tip3p, tip4p, tip5p

该步骤的任务
gmx pdb2gmx读取一个pdb文件（或者gro文件，或者一些数据库文件），向体系中添加氢原子（本来只有重原子），然后以gromacs的格式生成坐标文件（可以选择是否生成pdb文件），生成的文件可以作为后续的分析的输入文件。

• 选择力场
  gmx pdb2gmx将会在当前文件夹以及的子文件夹或者gromacs的文库中搜索力场文件（itp文件），默认条件下力场的选择是交互的（运行命令后会让选择），自己可以通过-ff命令来关闭这种交互指定一个力场，比如可以用-ff charmm27(小写)来指定力场，这样用一个命令行就可以实现，便于批量的分析。
  选择力场之后，所有的文件都将会从力场文件夹中读取，假如你想调整或者添加一个残基类型（注意：力场中不存在的pdb中的残基或导致无法执行命令），那么可以从力场文件夹中拷贝出来到当前工作目录，然后进行修改。比如，你想添加一个残基类型，那么你需要修改residuetype.dat文件，详细情况看第5章。

• 残基的数据库文件
  程序可以读取很多数据库文件，这些数据库文件使得程序能够产生特定的键，这个过程可以手动完成。这个程序可以加快用户选择LYS,ASP,GLU,CYS,HIS残基的类型。对于LYS，可以选择中性或者质子化的状态（在NZ上有两个质子）或者三个质子。对于HIS质子可以是ND1/NE2或者两者皆有。默认情况下这些选择都是自动完成的。对于HIS，这是基于最优的氢键构象完成的。氢键是根据简单的几何学标准来定义的（通过最大的氢键供受体的角度以及供受体的距离），这些参数可以通过-angle和-dist来设置。
  N端和C端的质子化状态可以通过交互式的方式来选择（利用-ter 标签），默认情况下末端是质子化的（NH3+,COO-）。一些力场支持两性离子的形式（对于一个条链只有一个残基的条件下），但是对于多肽这些选项不应该被选择。AMBER力场有位移的末端残基形式，这些与-ter参数是不兼容。你需要利用N/C前缀来对你的N/C端的残基进行处理从而使用这些形式。同时也可以利用末端残基的形式，比如ACE,NME。

注意事项
想要利用抗体进行动力学分析,先执行第一步gmx pdb2gmx
gmx pdb2gmx -f RIB.pdb -o RIB_merge.pdb -water spce -merge interactive -chainsep id_and_ter -ff charmm27

注意：

• 我使用的pdb文件是同源建模得到的抗体的文件，其中有H/L链之分。
• 这里没有用交互模式选择力场，而是用-ff指定。
• 命令执行后，会为每一条链生成对应的位置文件和力场文件
• 在执行命令之后会让选择是否将H/L链融合在一起，选择No即可，否则会把H链直接接在L链之后
• 观察输出的top文件，发现文件中已经对原子的位置文件名称（posre_Protein_chain_H.itp，posre_Protein_chain_L.itp）和原子的力场详细信息文件（topol_Protein_chain_H.itp，topol_Protein_chain_H.itp）进行了记录.
• 此时不产生gro文件，将会产生一个pdb文件,同时含有轻链和重链的信息（相当于是一个加H原子的过程）。

gmx editconf

命令行功能
该命令的主要功能是将pdb或者gro文件中获取得到的分子放在一个盒子中（该盒子作为模拟的单位）

参数解释

• -f （指定输入的结构文件）

• -o(指定输出结果文件)

• -bt (盒子类型)
  相当于box type的简写，用来控制要添加的盒子的类型，可选参数：

  • triclinic:三斜方体
  • cubic：正方体
  • dodeahedron:十二面体
  • octahedron:八面体

• -d （溶液与盒子的距离）
  由于盒子是等边长的立方体，因此，我们只用-d参数设定一个距离即可。单位为纳米。

• -c（放在盒子中心）
  表明将会把分子放在盒子的中心位置。假如没有指定改参数，那么就是单纯的把分子放在盒子中。

• -box
  指定box的边长，直接-box a b c，a，b,c表示三个数字。指定该参数之后，再生成的文件里的最后一行能看到指定的a,b,c。这里注意，当a,b,c不同时，系统会直接取最小值，作为盒子的边长（因为已经用-bt指定了盒子是一个立方体，边长必然相同）

• -center（指定盒子中心的位置）
  默认情况下盒子的中心为0,0,0,此时我们可以设置其他的坐标来作为中心，这样分子也会处在我们指定的这个中心上。
  比如
  我用以下命令加的盒子

gmx editconf -f RIB.gro -o newbox_RIB.gro -d 1.0 -c -bt cubic 

那么它生成的gro文件最后的位置信息如下

(Abs_Met_Oxydation) [fanxuezhe@fanxuezhe test_gro_pdb]$ tail newbox_RIB.gro 
  227ALA     CA 3373   4.453   2.372   1.658
  227ALA     HA 3374   4.354   2.387   1.665
  227ALA     CB 3375   4.499   2.399   1.514
  227ALA    HB1 3376   4.457   2.333   1.453
  227ALA    HB2 3377   4.473   2.492   1.487
  227ALA    HB3 3378   4.599   2.390   1.509
  227ALA      C 3379   4.501   2.233   1.697
  227ALA    OT1 3380   4.454   2.135   1.632
  227ALA    OT2 3381   4.584   2.222   1.791
   8.02702   8.02702   8.02702

此外我利用-center命令指定中心之后

gmx editconf -f RIB.gro -o newbox_RIB.gro -d 1.0 -bt cubic -center 2 2 2 

文件最后的信息如下

(Abs_Met_Oxydation) [fanxuezhe@fanxuezhe test_gro_pdb]$ tail newbox_RIB.gro 
  227ALA     CA 3373   2.439   0.358  -0.355
  227ALA     HA 3374   2.340   0.373  -0.348
  227ALA     CB 3375   2.485   0.385  -0.499
  227ALA    HB1 3376   2.443   0.319  -0.560
  227ALA    HB2 3377   2.459   0.478  -0.526
  227ALA    HB3 3378   2.585   0.376  -0.504
  227ALA      C 3379   2.487   0.219  -0.316
  227ALA    OT1 3380   2.440   0.121  -0.381
  227ALA    OT2 3381   2.570   0.208  -0.222
   8.02702   8.02702   8.02702

这里可以看出位置信息发生了巨大的变化，但是盒子的大小并未发生改变。

• angles
  坐标空间中的坐标之间的角度，默认为90,90,90。

gmx genion

命令功能
利用gmx grompp命令生成的tpr文件为系统添加离子。在执行gmx grompp命令之后，会对整个体系的电荷情况进行检查，正常情况下，体系的电荷应该是整数的，这里可以看到我们的系统呈现出一个小数，此时会出现一个warning。之后加离子的时候也是根据系统的正负电荷的状态来添加离子。比如这里我们是正6，那么我们可以添加6个氯离子在系统中从而形成电中性的系统。

NOTE 2 [file topol.top, line 48]:
  System has non-zero total charge: 6.000002
  Total charge should normally be an integer. See
  http://www.gromacs.org/Documentation/Floating_Point_Arithmetic
  for discussion on how close it should be to an integer.

相关参数

• -s
  指定gmx grompp命令生成的tpr文件

• -o
  `指定新生成的gro文件的名称，在这个文件里将会显示出新增加的氯离子。

• -nname
  表明要新增加的负电原子的种类，可以根据体系的性质进行添加，系统带5个正电，则可以加上带5个负电的离子，默认是CL离子

• -nn
  指定添加的带负电离子的个数

• -pname
  与-nname类似，用来表明要添加的正离子的类型，默认为NA离子（钠离子）

• -np
  与nn类似，表明要添加的正离子的个数

一个实例

gmx genion -s ion.tpr -o solv-RIB.gro -p topol.gro -nname CL -nn 6

命令执行后会让选择将来替换为离子的时候选择哪些分子来执行。

Reading file topol.tpr, VERSION 5.1.4 (single precision)
Reading file topol.tpr, VERSION 5.1.4 (single precision)
Will try to add 0 NA ions and 6 CL ions.
Select a continuous group of solvent molecules
Group     0 (         System) has 51192 elements
Group     1 (        Protein) has  3381 elements
Group     2 (      Protein-H) has  1724 elements
Group     3 (        C-alpha) has   227 elements
Group     4 (       Backbone) has   681 elements
Group     5 (      MainChain) has   906 elements
Group     6 (   MainChain+Cb) has  1111 elements
Group     7 (    MainChain+H) has  1125 elements
Group     8 (      SideChain) has  2256 elements
Group     9 (    SideChain-H) has   818 elements
Group    10 (    Prot-Masses) has  3381 elements
Group    11 (    non-Protein) has 47811 elements
Group    12 (          Water) has 47811 elements
Group    13 (            SOL) has 47811 elements
Group    14 (      non-Water) has  3381 elements
Select a group: 13
Selected 13: 'SOL'
Number of (3-atomic) solvent molecules: 15937

Processing topology

Back Off! I just backed up temp.topDDa7Ex to ./#temp.topDDa7Ex.1#
Replacing 6 solute molecules in topology file (topol.top)  by 0 NA and 6 CL ions.

Back Off! I just backed up topol.top to ./#topol.top.3#
Replacing solvent molecule 1600 (atom 8181) with CL
Replacing solvent molecule 4436 (atom 16689) with CL
Replacing solvent molecule 3496 (atom 13869) with CL
Replacing solvent molecule 14523 (atom 46950) with CL
Replacing solvent molecule 14838 (atom 47895) with CL
Replacing solvent molecule 2529 (atom 10968) with CL

我们这里选择水分子来替换为氯离子。

(Abs_Met_Oxydation) [fanxuezhe@fanxuezhe test_gro_pdb]$ tail solv_RIB.gro 
16158SOL     OW51172   7.876   7.872   7.968
16158SOL    HW151173   7.906   7.800   7.906
16158SOL    HW251174   7.837   7.832   8.051
16159CL      CL51175   0.940   3.560   6.490
16160CL      CL51176   1.592   5.353   7.211
16161CL      CL51177   0.083   6.844   6.608
16162CL      CL51178   2.272   0.813   4.975
16163CL      CL51179   6.136   7.644   4.609
16164CL      CL51180   7.755   0.063   4.342
   8.02702   8.02702   8.02702

通过检查文件的结尾处的信息，我们可以看出有6个水分子被替换为了CL。
下面是原来的gro文件的最后几行，可以看出编号从16158到16164的水分子被替换为氯离子。

(Abs_Met_Oxydation) [fanxuezhe@fanxuezhe test_gro_pdb]$ tail \#solv_RIB.gro.2# 
16162SOL     OW51184   7.523   7.793   7.481
16162SOL    HW151185   7.431   7.765   7.452
16162SOL    HW251186   7.554   7.870   7.425
16163SOL     OW51187   7.520   7.614   7.766
16163SOL    HW151188   7.503   7.697   7.712
16163SOL    HW251189   7.610   7.577   7.744
16164SOL     OW51190   7.876   7.872   7.968
16164SOL    HW151191   7.906   7.800   7.906
16164SOL    HW251192   7.837   7.832   8.051
   8.02702   8.02702   8.02702

gmx grompp

命令行功能
将能量最小化的mdp文件转化为可执行的tpr文件。
执行流程
首先建立一个mdp文件，该文件的作用是控制能量最小化的流程，随后利用gmx grompp命令将这个文件转化为可执行文件（tpr文件）。
一个实例（生成tpr文件）

integrator     =  steep
emtol          =  1000
emstep         =  0.01
nsteps         =  10000
nstlist        =  1
cutoff-scheme  =  Verlet
ns_type        =  grid
coulombtype    =  PME 
rcoulomb       =  1.0 
rvdw           =  1.0 
pbc            =  xyz 

接下来执行命令生成可执行的tpr文件

gmx grompp -f ion.mdp -c solv_RIB.gro  -p topol.top

然后在文件中可以看到topol.tpr的输出文件

相关参数

• -f: 指定mdp文件，这个文件需要自己编辑
• -c: 指定gro文件，gro文件指的是加完水分子的gro文件
• -p: 指定top文件，这里的top文件指的是加完水分子的top文件
• -o: 指定输出文件

gmx mdrun

使用 gmx mdrun 的 -rerun 选项指定原来的轨迹文件再跑一次模拟, 这个过程很快.

gmx genrestr

生成urea分子的restraction文件

gmx genrestr -f urea.pdb -o porse_urea.itp -fc 1000 1000 1000
## 或者
gmx genrestr -f jz4.gro -o posre_jz4.itp -fc 1000 1000 1000

gmx insert-molecules

功能：插入输入文件中-nmol指定的系统的副本-ci。插入发生在给出的溶质构象的空白空间中-f，或者发生在给出的空盒中-box。同时指定-f 和的-box行为类似于-f，但是在插入之前在溶质周围放置一个新框。存在的任何速度都将被丢弃。

指定输入文件的选项：

• -f [<.gro / .g96 /…>]（protein.gro）（可选）
  现有的配置插入到：GRO G96 PDB BRK耳鼻喉科ESP TPR
• -ci [<.gro / .g96 /…>]（insert.gro）
  配置插入：GRO G96 PDB BRK耳鼻喉科ESP TPR
• -ip [<.dat>]（positions.dat）（可选）
  预定义的插入试验位置
• -n [<.ndx>]（index.ndx）（可选）
  额外索引组
  指定输出文件的选项：
• -o [<.gro / .g96 /…>]（out.gro）
  插入后的输出配置：gro g96 pdb brk ent esp
  其他选项：
• -replace <选择>
  重叠时可以去除的原子
• -sf <文件>
  提供文件中的选择
• -selrpos <枚举>（原子）
  选择参考位置：atom，res_com，res_cog，mol_com，mol_cog，whole_res_com，whole_res_cog，whole_mol_com，whole_mol_cog，part_res_com，part_res_cog，part_mol_com，part_mol_cog，dyn_res_com，dyn_res_cog，dyndres_comg，dyn_res_cog，
• -box <向量>（0 0 0）
  包装盒尺寸（单位：nm）
• -nmol （0）
  要插入的额外分子数
• -try （10）
  尝试插入-nmol倍-try倍
• -seed （0）
  随机生成器种子（0表示生成）
• -radius <真实>（0.105）
  默认范德华距离
• -scale <真实>（0.57）
  在share / gromacs / top / vdwradii.dat中从数据库乘以范德华半径的比例因子。默认值为0.57，则水中蛋白质的密度接近1000 g / l。
• -dr <向量>（0 0 0）
  允许从-ip文件位置以x / y / z位移
• -rot <枚举>（xyz）
  随机旋转插入的分子：xyz，z，无

分析

gmx distance

功能：计算原子之间的距离

gmx distance -s md_0_10.tpr -f md_0_10_center.xtc -select 'resname "JZ4" and name OAB plus resid 102 and name OE1' -oall

gmx distance -s -f -select "com of resnr 1 to 4 plus com of resnr 7 to 9" -oav -oall

gmx check

校验轨迹信息

gmx check-f trj_40-50ns.xtc #轨迹检查