1. Gromacs是用来干什么的？

分子动力学模拟是分子模拟中最接近实验条件的模拟方法, 能够从原子层面给出体系的微观演化过程, 直观地展示实验现象发生的机理与规律，gromacs就是一个用于分子动力学模拟的程序。

2. Gromacs的步骤

        以T4溶菌酶L99A/M102Q为例

2.1 初始化结构

        通过实验数据或者某些工具得到体系内的每一个分子的初始结构坐标文件。之后, 我们需要按我们的想法把这些分子按照一定的规则或是随机的放在一起, 从而得到整个体系的初始结构, 这也是我们模拟的输入文件。(重点关注结构文件是否有问题）

• gro: GROMACS格式 
• g96: GROMOS-96格式
• pdb:  蛋白质数据库格式 

通用结构格式: gro, g96, pdb, tpr, tpb或tpa

结构+质量(db): tpr, tpb, tpa, gro, g96或pdb. 用作分析工具的输入文件. 当使用gro或pdb时, 从质量数据库中读取近似质量.

2.2 初始化参数

有了结构输入文件, 还需要力场参数输入文件。

• mdp ：运行参数文件, 用作gmx grompp与gmx convert-tpr的输入文件

• m2p ：用作gmx xpm2ps的输入文件

2.3 生成拓扑文件

分两步来准备体系的拓扑:

1. 使用gmx pdb2gmx准备蛋白质的拓扑

2. 使用外部工具准备配体的拓扑

上面的命令将从3HTB_clean.pdb文件中抽取出含JZ4的行, 保存到JZ4.pdb文件中. 

 gmx pdb2gmx将会处理结构文件, 并提示选择力场如下图所示。（-f是打开，-o为生成）pdb2gmx 的目的是生成三个文件： 分子的拓扑（3HTB_processed.gro)、 位置约束文件(.top)、 后处理的结构文件(drg.itp)。

• top:  体系拓扑(文本) 
• itp: 包含拓扑(文本) 
• rtp: 残基拓扑(文本) 
• ndx: 索引文件

 

 生成拓扑完成。

2.4 创建复合物

根据gmx pdb2gmx的输出, 我们已经有了一个名称为3HTB_processed.gro的文件, 它包含了经处理后与力场相容的蛋白质结构. 我们也有一个来自PRODRG的文件jz4.gro, 其中包含了所有需要的H原子. 将3HTB_processed.gro文件另存为一个待处理的新文件, 如complex.gro, 因为将JZ4添加到蛋白质就可以得到蛋白配体的复合物. 接下来, 复制jz4.gro的坐标部分粘贴到complex.gro文件中。

2.5 建立盒子、添加溶剂

• -f: 输入复合物的结构文件

• -o: 输出复合物的结构文件, 包含新盒子的信息

• -bt: 定义盒子形状: triclinic(三斜体), cubic(立方体), dodecahedron(十二面体), octahedron(八面体)

• -d: 定义盒子边缘距离分子边缘的距离(单位: nm), 通常不能小于0.85 nm

  gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro

• -cp: 输入含盒子的结构文件

• -cs: 指定spc水盒子

• -p: 输出含溶剂的复合物拓扑文件

• -o: 输出含溶剂的复合物结构文件

• -cp: 输入含盒子的结构文件

• -cs: 指定spc水盒子

• -p: 输出含溶剂的复合物拓扑文件

• -o: 输出含溶剂的复合物结构文件

2.6 添加离子

现在有了一个已经溶剂化的体系, 其中包含一个带电的蛋白质。由于生命体系中不存在净电荷, 我们必须在体系中添加离子。

使用gmx grompp整合.tpr文件, 可利用任何.mdp文件，此处使用了ions.mdp。

 将所得的.tpr文件传给gmx genion

• -s ions.tpr: 输入预处理得到的tpr文件

• -o solv_ions.gro: 输出添加了抗衡离子后的结构文件

• -p topol.top: 输出拓扑文件

• -pname NA: 添加的阳离子类型

• -nname CL: 添加的阴离子类型

• -nn 6: 添加的阴离子个数

 2.7 能量最小化

        由于初始构象可能会存在两个原子靠得太近的情况(称之为bad contact), 所以需要在正式模拟开始的第一步对体系进行能量最小化。比较常用的能量最小化方法有两种, 最速下降法和共轭梯度法。

使用gmx grompp并利用输入参数文件来创建一个二进制输出文件 

 运行能量最小化

2.8 平衡模拟

需要设置适当的模拟参数, 并且保证这些参数的设置与力场的构造过程相一致。常见的模拟思路是, 先在NVT下限制住你的溶质(剂)做限制性模拟, 这是一个升温的过程, 当温度达到你设定的温度后, 接着做NPT模拟, 此过程将调整体系的压强进而使体系密度收敛。

2.8.1 NVT模拟

平衡我们的蛋白质配体复合物需要有些特殊考虑:

1. 对配体施加限制

2. 处理温度耦合组

 限制配体，PRODRG并不会为我们的配体生成一个类似posre.itp的文件, 但GROMACS提供的gmx genrestr模块可帮助我们完成这项工作. 利用PRODRG给出的jz4.gro文件运行gmx genrestr。

• -f jz4.gro: 输入配体的结构文件

• -o posre_jz4.itp: 输出配体施加位置限制的itp文件

• -fc 1000 1000 1000: 位置限制的力常数(单位: kJ/mol-nm2)

 在拓扑文件中包含这些信息

接下来进行热浴。防止体系爆开（ 因为对只有几个原子的组(即, JZ4和CL), 在控制其动能的涨落时, 温度耦合算法不够稳定.）在温度耦合时将JZ4当作蛋白质的一部分。 为此, 我们需要一个特殊的组, 其中包含蛋白质和JZ4。可以使用gmx make_ndx模块来完成这个工作, 只要提供完整体系的任意坐标文件即可。

进行NVT平衡：

        gmx grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr

运行：

        gmx mdrun -deffnm nvt

2.8.2 NPT平衡

2.9 分子模拟

 

 

3.其他参考

GROMACS文件类型|Jerkwin

科学网-GROMACS中文教程之: 蛋白质配体复合物-李继存的博文

 GROMACS程序文档|Jerkwin