GroMacs 跑起来！
开始之前感谢一下Justin A. Lemkul, Ph.D.教授的GROMACS Tutorials。现在已经几乎成为新手必学行业标准了。强烈推荐可以去上他的课。这篇blog是我根据他的教程总结出来，并加上了自动化的部分，让一切更加方便！
我们会采用3NPO pdb file 和一个autodock vina对接后的pdbqt文件开始。你可以从下面下载他们，但也可以参考我之前的教程一步一步拿到自己的pdbqt flie：
一、准备pdb file
你会发现autodock vina dock之后的pdbqt里面有多个构象，他们都被MODEL分割。如下图：
所以我们处理的第一步就是分割，把这个含有多个model的变成单独的一个model。利用下面几行sh就可以轻松拿到第一个构象。当然，也可以直接用手复制粘贴出来。不过可能使用vina的vina_split会更可靠一些
vina_split --input dock_result.pdbqt
mv dock_result_ligand_01.pdbqt ligand.pdbqt # rename to liagnd
rm dock_result_ligand_* # remove all other pdbqt
拿到第一个构象
pdbqt是autodock专用格式，我们下面转换成ligand常用的通用格式：
obabel -ipdbqt ligand.pdbqt -omol2 -O ligand.mol2 -p 7.0转换为mol2
3npo.pdb)还有一个mol2 file!这里有一个我写好的sh文件，可以用来自动化这个流程，跟着问答走就好
二、准备拓扑结构
如果你一步一步做上来，肯定已经疑问很久了，pdb file里面只有x y z坐标，那些键是哪里来的？为啥使用pymol打开他有键的信息？没错，pymol是猜的！所以经常显示错误。
拓扑结构生成还是比较耗费算力的，下面就来生成一个
2.1 ligand topology
ligand相对复杂，因为topology并没有现成的。毕竟你想想，氨基酸就那么几个，实际蛋白在生成拓扑结构的时候就是做简单的查字典，但是小分子可不一样，小分子变化多端没有字典可查，而是要老老实实的算薛定谔方程。这个过程非常复杂，还好有一个好用的工具帮我们做到，他就是acpype.
acpype -i ligand.mol2 -c bcc没错，一行命令，10分钟等待
这里我们指定了AM1-BCC 电荷模型，它是半经验半计算的，让他比传统RESP速度大幅提升。我们小分子虽然小，但是实际上也不小，使用BCC模型非常常见。
里面含有_GMX为后缀的就是给gromacs用的文件。里面有用的是一个gro文件和一个itp。
2.2 protein topology
下面是相对简单的蛋白拓扑图。Gromacs直接有对应的工具，直接使用即可
gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top -tergmx pdb2gmx 一句话搞定
命令中也可以不加-ter而是-ignh让他直接根据默认方式加H，如果有报错的话。
三、合并拓扑和坐标
现在我们有两个分开的topology和gro的xyz文件，在真正跑gromacs的时候，我们要一个文件包含所有信息。这一步似乎基本上只能手动，不过我这里提供脚本可以变得全自动，但是我们从原理开始
gromacs最终跑需要：
- 一个gro文件，也就是xyz坐标文件，类似pdb。
- 一个top文件，也就是topology。itp也是top文件，并且可以被top include
现在我们逐个搭建：
3.1 gro文件搭建
gro文件本质上就是xyz，所以我们直接粘贴过去就好。复制整个ligand_GMX.gro里面的信息，粘贴在protein的gro的最下面，box信息之前，并且修改第二行的总原子个数。
下面是等效的sh file用于自动化操作
head -n -1 protein.gro > complex.gro
sed '1,2d;$d' ligand_GMX.gro >> complex.gro
tail -n 1 protein.gro >> complex.gro
PROT_ATOMS=$(sed -n '2p' protein.gro | awk '{print $1}')
LIG_ATOMS=$(sed -n '2p' ligand_GMX.gro | awk '{print $1}')
TOTAL_ATOMS=$((PROT_ATOMS + LIG_ATOMS))
echo "Total atoms will be: $TOTAL_ATOMS"
sed -i "2s/.*/$TOTAL_ATOMS/" complex.gro3.2 top文件搭建
top文件相对简单。只要在力场文件后面加上include ligand的itp。并且在最后一行加上1个ligand组成成分即可。下面是等效的sh文件：
echo "ligand 1" >> topol.top
sed -i '/#include "amber14sb.ff\/forcefield.itp"/a #include "ligand_GMX.itp"' topol.top四、搭建溶液体系并加ion中和
从这里开始流程就相对固定，可以高度自动化了。并且软件会处理的比较好，简单了解流程即可。
首先要定义一个有溶液的框框：
gmx editconf -f complex.gro -o complex_box.gro -c -d 1.0 -bt dodecahedron
gmx solvate -cp complex_box.gro -cs spc216.gro -o complex_solv.gro -p topol.top定义box并加入溶液
下面添加ion中和可能的离子。gromacs体系是不带电的，通常情况下如果你的protein+ligend带正电，就加Cl-，如果是负电就加Na+。
这里需要做grompp所以需要一个配置文件。autodock中配置文件是config.txt,这里的配置文件就是**.mdp。对于这个ion.mdp可以用下面的链接下载，这玩意比较通用，除非你是大神否则没必要改里面的内容。
gmx grompp -f ion.mdp -c complex_solv.gro -p topol.top -o ions.tpr
printf "SOL\n" | gmx genion -s ions.tpr -o complex_neutral.gro -p topol.top -pname NA -nname CL -neutral 五、控制ligand
不能让ligand在平衡模拟阶段乱跑，在平衡模拟阶段水分子并不均匀，很有可能和ligand重叠一开始就直接把ligand撞飞。为了防止这个情况发生，我们采用强锁定ligand的方式
首先我们选中所有ligand做成一个组，然后生成一个规则名为posre_ligand.itp，注意这是一个可以被添加进top文件的拓扑补充文件，后缀是之前就见过的itp。
printf "0 & ! a H*\nq\n" | gmx make_ndx -f ligand_GMX.gro -o index_ligand.ndx
printf "3\n" | gmx genrestr -f ligand_GMX.gro -n index_ligand.ndx -o posre_ligand.itp -fc 1000 1000 1000 然后include这个itp file，这里的意思是如果配置文件mdp里面写了POSRES_LIG就include，如果没有写就不include。因为我们只有平衡阶段要引入他，正式模拟还是希望小分子也会运动的，所以要这样写
; Ligand position restraints
#ifdef POSRES_LIG
#include "posre_ligand.itp"
#endif或者使用sh自动化完成添加的过程
sed -i '/#include "ligand_GMX.itp"/a \\n; Ligand position restraints\n#ifdef POSRES_LIG\n#include "posre_ligand.itp"\n#endif' topol.top对于从三到五也有一个自动化脚本可以一键完成所有事情, 感兴趣的可以试试看。也是跟着提示一步一步做即可。
六、开始模拟！
模拟主要分为两个阶段，平衡模拟和正式模拟。平衡模拟在可控的环境下进行（也就是我们刚刚做的控制ligand步骤），让水分子、温度、压力完全达到配置文件mdp中写出来的那样。正式模拟就让所有原子自由运动了。
在模拟之前，我有很多配置文件要给你，你可以打开看他们，在大多数情况下不需要更改，这都是非常基础和通用的模板，最可能要改的就是步数和保存间隔了。如果你希望跑的多一些，就把步数调整高一些。所有的mdp文档也都是从两位教授的教程中拿到的。
第一步是能量最小化：
gmx grompp -f inputs/em.mdp -c complex_neutral.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em第二步是Canonical Ensemble，压力P 和 能量E 变化，体积和温度不变
gmx grompp -f inputs/nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -v -deffnm nvt下一步是Isothermal-Isobaric Ensemble，压力和温度都不变，能量和体积可变
gmx grompp -f inputs/npt.mdp -c nvt.gro -r nvt.gro -p topol.top -o npt.tpr
gmx mdrun -v -deffnm npt这一共3步平衡模拟以后，就可以开始正式模拟了！开始前记得修改到自己想要跑的步数。目前mdp里面是跑200ns的，也就是发文章最低标准。
gmx grompp -f inputs/md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_noPBC.tpr
gmx mdrun -v -deffnm md然后？
然后等着吧，如果你要求的步数多的话，几天都是有可能的。没错，几天也就能跑出来蛋白质在几纳秒里面的变化，不得不感慨，大自然真的是一台超级超级计算机。
量子计算其实就是再利用大自然给我们做计算，我跑完gromacs真的有点想去学量子计算了！
另外，跑出来的文件就可以用来进行轨迹分析等等了。
