Gromacs操作心得

文章来源：“分子动力学”公众号

链接：https://mp.weixin.qq.com/s/nQ0MdLmIXHPB8VLSe7CTZQ

一、进行模拟的基本步骤

以下是进行分子动力学模拟的最基本步骤，来自Gromacs维基。具体的步骤可能会因为模拟目标的不同而存在差异，这里只做基本的概述。

十分清楚要什么样的模拟系统性质和现象，明白模拟的目标。

选择适当的工具以实现模拟目标。这需要熟知其他研究人员对相似模拟系统的模拟方法等资料，要好好读别人的论文，主要包括：

模拟使用的软件，当然还要考虑使用的力场，有时因为力场关系，需要更改使用的软件。
力场描述了系统中原子的性质和相互作用，选择一个适合自己研究目标的力场至关重要，还是需要多看论文。根据sensenbobo血泪史和别人的经验，模拟核酸，可以采用amber力场，水分子可以使用TIP3P；模拟蛋白质，可以采用OPLS-AA，水分子使用TIP4P（重申一下，我本人没有跑过测试）。
获取模拟分子的坐标文件（PDB），或者自己生成一个。自己生成很好啊，强烈鼓励，但是要十分明白有一定的危险性，要明白每一步模拟系统的现象是为什么。
根据模拟软件，生成模拟系统初步坐标文件。
获取或者生成系统的拓扑文件，比如Gromacs的pdb2gmx命令，或者使用网络服务器PRODRG（google一下这几个字母），amber的leap命令和antechamber，NAMD的psfgen或者VMD。请参考相关手册。
给模拟系统添加一个盒子（盒子的翻译方法我十分不喜欢，箱子显得太大，边界显得太难理解，但是盒子听起来想到蛋糕，我们是严肃的科学家。娃哈哈。）在盒子中添加水分子，添加离子如NA和CL之类（添油加醋，生活就是这样进行的）。添加任何东西到你的系统中，就要相应的修改一下拓扑文件，以保持与系统对应，不然会出错哦。
进行能量最优化。这样做是消去系统中原子之间距离太近造成的高能量。必要的时候，分子可以现在真空中能量最优话，然后再做溶剂环境中能量最优话。
选择恰当的模拟参数模拟平衡模拟系统。需要慎重考虑力场的作用，必要的时候保持系统体积温度不变（NVT系综）进行位置限制性平衡系统，然后在保持系统温度压强不变（NPT系综）修正系统的密度。然后再选择正确的系综（NPT，NVT，NET，BT？？）进行数据收集。（真烦！）
跑足够长时间的模拟，使系统达到理想的平衡。使用模拟软件的工具或者其他软件观察系统的性质，如温度，能量瞪瞪。使用适当的参数进行数据采集模拟，注意模拟衔接过程中各个系统性质保持不变，如果原子速度之类。这里还要考虑力场的作用，同时要考虑怎么测量描述模拟目标。
模拟足够长的时间使模拟目标明显出现（可能也出现跟你期望相反的结果，出现率应该百分之五十，加上个人感情因素，可能会达到百分之八十。保重！）。
分析模拟结果，解释得到结果。要好好利用视图工具，同时要十分注意漂亮的三维图说明不了什么问题。图表，比较，统计，一个都不能少。就这样吧。

二、Gromacs使用－新手入门篇

其实gromacs的思路比较简单，大体就是以下的流程,在此先对em进行一下总结：

2.1 convert the pdb-file to a gromacs structure file(.gro / .pdb) and a gromacs topology file(.top)

这一步需要gromacs的预处理程序pdb2gmx 标准的输入命令应该是：

pdb2gmx -ignh -f speptide.pdb -p speptide.top -o speptide.gro
这个命令之后有一个force field的选择 系统给出了11种力场，具体以后实际计算中如何选择力场还是应该多参阅文献和gromacs说明书的。这步之后可以在屏幕看到：

--------- PLEASE NOTE ------------
You have succesfully generated a topology from: speptide.pdb.
The G43a1 force field and the spc water model are used.
Note that the default mechanism for selecting a force fields has
changed, starting from GROMACS version 3.2.0
--------- ETON ESAELP ------------
至此，完成了从pdb到gro和top的转变。另外，我也看过有的作这一步命令时候并不进行top文件的转换，因为在genbox（对盒子加水溶剂化时还要对top文件进行修改）

2.2 solvate the peptide in water

这一步需要用到editconf和genbox两个程序，命令的标准书写应该是：

editconf -bt cubic -f speptide.gro -o speptide_box.gro -c -d 0.5
可以在屏幕上看到：

Read 191 atoms
Volume: 8.17315 nm^3, corresponds to roughly 3600 electrons
No velocities found

system size :
1.774
3.372
1.367 (nm)

diameter
:
3.517
(nm)

center
:
2.650
1.453
2.417 (nm)

box vectors :
1.774
3.372
1.366 (nm)

box angles
:
90.00
90.00
90.00 (degrees)

box volume
:
8.17
(nm^3)

shift

-0.391
0.806 -0.158 (nm)
new center
:
2.259
2.259
2.259 (nm)
new box vectors :
4.517
4.517
4.517 (nm)
new box angles
:
90.00
90.00
90.00 (degrees)
new box volume
:
92.18
(nm^3)
因为程序的默认，所以这一步命令可以写作：

editconf –f speptide –o –d 0.5
这里－f之后的sp文件是默认gro文件或者pdb文件，－o后面生成的output file是放入盒子的体系的gro文件 按照这样的格式，输出的文件被默认命名为－out.gro

体系放入盒子之后在盒子中注入水

genbox -cp speptide_box.gro -cs -p speptide.top -o speptide_water.gro
也可简写为 ：

genbox –cp out –cs –p sprptide –o b4em

Screen–
Output configuration contains 9035 atoms in 2967 residues
Volume
:
92.1813 (nm^3)
Density
:
994.449 (g/l)
Number of SOL molecules:
2948
之后看加水之后的top文件

用命令

tail speptide.top
可以看到：

[ system ]
; Name
Protein in water

[ molecules ]
; Compound
#mols
Protein
1
SOL
2948
其中SOL后面的数字2948就是genbox加入盒子的水分子的数目。

Editconf程序的另一个用途是讲gro文件转化回pdb这时可以讲speptide_water.gro转化回pdb观察

editconf –f speptide_water.gro –o speptide.pdb
拖回本机 用spbdv或者vmd观察 。

下一步是生成索引文件，需使用 make _ ndx 来生成一个。

make_ndx –f b4em

2.3 perform an enery minimization of the peptide in solvent

现在模拟系统已经差不多准备好了。在我们开始动力学之前，我们必须进行能量最小化，以减轻系统中可能存在的任何不良接触(原子重叠导致明显的排斥，从而在模拟中引起数值问题)。

grompp –v –f em.mdp –c speptide_water.gro –p speptide.top –o speptide_em.tpr
或者简单地使用这个命令:

grompp –v –f em –c b4em –p speptide –o em
在这个命令之后，坏的接触已经被移除。所以我们现在可以尽量减少能量消耗。

mdrun –v –s speptide_em.tpr –o speptide_em.trr –c after_em.gro –g emlog.log
或者

mdrun –v –s em –o em –c after_em –g emlog
屏幕可见：

Steepest Descents:

Tolerance (Fmax)

2.00000e+03
Number of steps

100
Step=
0, Dmax= 1.0e-02 nm, Epot= -8.02562e+04 Fmax= 2.70468e+04, atom= 6792
Step=
1, Dmax= 1.0e-02 nm, Epot= -8.69641e+04 Fmax= 1.24531e+04, atom= 6657
Step=
2, Dmax= 1.2e-02 nm, Epot= -9.48067e+04 Fmax= 5.54072e+03, atom= 9009
Step=
3, Dmax= 1.4e-02 nm, Epot= -1.02743e+05 Fmax= 2.86957e+03, atom= 8766
Step=
4, Dmax= 1.7e-02 nm, Epot= -1.08941e+05 Fmax= 1.53075e+04, atom= 150
Step=
5, Dmax= 2.1e-02 nm, Epot= -1.10083e+05 Fmax= 1.31653e+04, atom= 150
Step=
6, Dmax= 2.5e-02 nm, Epot= -1.10564e+05 Fmax= 2.16454e+04, atom= 150
Step=
7, Dmax= 3.0e-02 nm, Epot= -1.11528e+05 Fmax= 1.82398e+04, atom= 150
Step=
9, Dmax= 1.8e-02 nm, Epot= -1.12767e+05 Fmax= 3.87548e+03, atom= 120
Step=
11, Dmax= 1.1e-02 nm, Epot= -1.13598e+05 Fmax= 1.23013e+04, atom= 120
Step=
12, Dmax= 1.3e-02 nm, Epot= -1.14464e+05 Fmax= 5.37656e+03, atom= 120
Step=
13, Dmax= 1.5e-02 nm, Epot= -1.14787e+05 Fmax= 1.74494e+04, atom= 120
Step=
14, Dmax= 1.9e-02 nm, Epot= -1.15874e+05 Fmax= 8.35820e+03, atom= 120
Step=
16, Dmax= 1.1e-02 nm, Epot= -1.16420e+05 Fmax= 6.65421e+03, atom= 120
Step=
17, Dmax= 1.3e-02 nm, Epot= -1.16807e+05 Fmax= 1.08106e+04, atom= 120
Step=
18, Dmax= 1.6e-02 nm, Epot= -1.17248e+05 Fmax= 1.12834e+04, atom= 120
Step=
19, Dmax= 1.9e-02 nm, Epot= -1.17425e+05 Fmax= 1.45145e+04, atom= 120
Step=
20, Dmax= 2.3e-02 nm, Epot= -1.17590e+05 Fmax= 1.67089e+04, atom= 120
Step=
22, Dmax= 1.4e-02 nm, Epot= -1.18845e+05 Fmax= 2.65542e+03, atom= 121
Step=
23, Dmax= 1.7e-02 nm, Epot= -1.19198e+05 Fmax= 2.39072e+04, atom= 120
Step=
24, Dmax= 2.0e-02 nm, Epot= -1.20588e+05 Fmax= 6.17310e+03, atom= 120
Step=
26, Dmax= 1.2e-02 nm, Epot= -1.20901e+05 Fmax= 1.04224e+04, atom= 120
Step=
27, Dmax= 1.4e-02 nm, Epot= -1.21229e+05 Fmax= 8.89909e+03, atom= 120
Step=
28, Dmax= 1.7e-02 nm, Epot= -1.21346e+05 Fmax= 1.51125e+04, atom= 120
Step=
29, Dmax= 2.1e-02 nm, Epot= -1.21686e+05 Fmax= 1.29390e+04, atom= 120
Step=
31, Dmax= 1.2e-02 nm, Epot= -1.22216e+05 Fmax= 3.21419e+03, atom= 120
Step=
33, Dmax= 7.5e-03 nm, Epot= -1.22523e+05 Fmax= 6.52631e+03, atom= 120
Step=
34, Dmax= 8.9e-03 nm, Epot= -1.22796e+05 Fmax= 6.31916e+03, atom= 120
Step=
35, Dmax= 1.1e-02 nm, Epot= -1.22998e+05 Fmax= 8.04922e+03, atom= 120
Step=
36, Dmax= 1.3e-02 nm, Epot= -1.23167e+05 Fmax= 9.93619e+03, atom= 120
Step=
37, Dmax= 1.5e-02 nm, Epot= -1.23317e+05 Fmax= 1.10029e+04, atom= 120
Step=
39, Dmax= 9.3e-03 nm, Epot= -1.23857e+05 Fmax= 1.88770e+03, atom= 810
writing lowest energy coordinates.
Steepest Descents converged to Fmax < 2000 in 40 steps
Potential Energy
= -1.2385705e+05
Maximum force

1.8876959e+03 on atom 81
Norm of force

1.3663420e+04
NOTICE!!!
If the potential enery after minimization is lower than -1.1e+05kJ/mol, it is acceptable and the structure can be used for MD caculations

三、gromac的一些笔记

用了一个星期，基本操作差不多了。自顶向下大概是这样：1．核心步骤当然是mdrun了，当然后续的数据分析才是对理论水平的要求，先不讨论。这里需要一个tpr输入文件——由grompp生成的计算体系和计算控制详细信息的二进制文件。计算生成trr轨迹文件、xtc压缩的轨道文件、edr能量等状态参量记录文件、log计算记录等。后续处理靠这些输出文件和 gmx提供的命令。2.grompp预处理文件，将三种文件整合到一个tpr文件作为mdrun的输入。分别是：

mdp文件：计算控制文件，里面要给md模拟时的计算参数加以描述，控制时间、精度、温度什么的。
gro文件：分子坐标文件，记录要模拟的体系中全部（要计算的）原子的坐标、名称等。
top文件：力场参数文件，记录gro文件中的原子之间相互作用立场参数和原子分子类别信息的文件。
3．mdp文件一般可以看手册根据需要自己填写，或者看模版。4．gro、top文件，用editconf加周期边界元胞，genbox加溶剂。那么最初的gro、top文件如何得到呢？

如果是从网站上直接当下来的pdb文件，用pdb2gmx转化。（有时需要先去掉里面的水，有时需要加-ignh不考虑多余的氢，有时需要加一些离子（genion）来平衡蛋白质上的电荷）
要是还要加入自己的分子，比如在MS什么的软件中先画了一个，存成pdb文件，这时不能直接转，因为残基gmx不认识，立场库里面没有描述，需要自己写top文件。不过这很累，还要自己计算，又不会，怎么办呢，网上推荐用ProDrg软件（只找到网络版的网址google一下就有），把pdb文件 copy进去，run一下就出来gro，和itp文件了（还有其他的一坨，用不到）。把自己gro文件里的内容copy到转出来的gro里，总原子数和序号改一下，再把itp文件include到老top文件中，[ molecules ] 标签下面添上加入的分子名和数，就可以grompp了。