gromacs蛋白质-配体复合物（tutorial 5）-初学-模拟流程记录

建议在做此教程前，先完成tutorial 1比较好

获取蛋白PDB文件：

直接百度pdb，进入pdb数据库，搜索3HTB

在页面右上角download files下选择PDB Format（我是win10）

去除结晶水和其他配体：

删除所有的CONECT，HOH，PO4，BME，将文件另存为3htb_clean.pdb

提取pdb文件中的配体：

在相同的工作目录下新建文件jz4.pdb，将3htb_clean.pdb文件中的JZ4所在行剪切到jz4.pdb中

将3htb_clean.pdb保存（ctrl+s），此时拉到文件最下端，按理说应该如下显示：

下载力场文件：

点击此处mackerell Lab-gromacs进行跳转，下载charmm36立场文件（charmm36-jul2022.ff）

拉到页面下方下载其配套的脚本文件（我这里的python版本是3.x），所以我下载第三个，至于如何查看自己的python看接下来的步骤

python版本查询：

由于苯人python是很早之前安装的，忘记版本是多少了，于是win+R，打开cmd，输入python -V，回车，得到，版本为3.8.8

获取蛋白拓扑文件：

解压刚刚下载的压缩包

此时工作目录下应该有这些文件才对，这个时候可以把下载压缩包剪切放在其他文件夹，或者删掉（当然不放也行，我只是强迫症）

在工作目录的地址栏输入cmd，输入gmx回车，出现名人名言表示可以直接输入命令了（这个是需要先配置好环境变量才可以的）

输入命令：gmx pdb2gmx -f 3htb_clean.pdb -o 3htb_processed.gro，选择力场1（就是我们刚刚下载的），回车

之后选择水模型TIP3P，这个是适用charmm力场的水模型，注意不同立场适用的水模型不一样，要查了之后再使用，选择1，回车

但是此时出现报错：atom C1 not found in buiding block 1MET while combining tdb and rtp

百度之后查看了这一篇文章：关于charmm36-jul2022.ff报错的问题（Fatal error:atom C1 not found in buid - 哔哩哔哩 (bilibili.com)

原来是2022这个文件里面力场对端基的处理不一样

解决方法一：再去官网下载2021版本的：charmm36-jul2021.ff 

解决方法二：加上-ter命令，gmx pdb2gmx -f 3htb_clean.pdb -o 3htb_processed.gro -ter（-ter选项可以交互的选择蛋白质N 端和 C 端的质子化状态，一般默认NH3+和COO-，这里可能2021和2022默认的处理端基的方式不一样？可能是这里氨基端那个MET？有问题？不是很懂这个里面MET的报错，如果有大佬路过可以解释一下吗）

重新输入命令：

gmx pdb2gmx -f 3htb_clean.pdb -o 3htb_processed.gro -ter

选择力场，水模型（和上面一样的操作）

选择开始端（氮端）的处理：

接着选择结束端（碳端）的处理：

等待输出名人名言之后，就算成功了，此时工作目录下会出现三个文件

获取配体拓扑文件（这一步比较长）：

官网推荐采用：

Avogadro这个来处理配体加氢（为什么要加氢呢，因为charmm是全原子力场，所以H也要显示出来，一般晶体里面是没有H的，所以我们得自己加），进入官网点击下载

安装好之后打开，将jz4.pdb拖入页面，build→add hydrogens加H

加H之后的分子如下：

将文件另存为mol2格式，File→save as...→jz4.mol2

修改mol2文件，将文件顶部的*****改为JZ4

修改mol2文件的残基名称和编号，因为mol2文件里面只有一个分子，所以残基名称和编号都要改为统一的，将JZ164和UNL1改为JZ4，167全部改为1，ctrl+s保存文件

   →→→→→→→

调整mol2文件的@<TRIPOS>BOND成键顺序，使用官网下载的脚本sort_mol2_bonds.pl ，将处理后的文件改为jz4_fix.mol2

ps：这里说明一下怎么看成键顺序是否有问题（参考）：@<TRIPOS>BOND字段第一列为键的序号，第二列为成键第一个原子的序号，第三列为成键第二个原子的序号，第四列为键的类型。正确的需要应该是递增的，例如第二列最开始的原子需要应该是该列中最小的，例如查看jz4.mol2文件可以发现，第二列中最小的序号为1，对应的第二个原子序号分别是（1 9）（1 13）（1 12）（1 11），第三列要在序号为1的基础上再进行排列，所以第三列应该是9，11，12，13

我们可以将修正前后的文件进行对比：

→→→→→→

使用CGenFF生成拓扑文件，导入iz4_fix.mol2，网站将会生成str文件（这里我不能登录CGenFF网站，注册账号不成功，如果有大佬知道怎么办麻烦指教！）我直接用官网处理好的str文件

一般来说采用CGenFF生成的str文件需要先查看惩罚分数，我们这里两个都小于1说明是可以的，太大了就不行，需要自己调整，但是我不知道怎么调，恳请大佬指教

接着采用python脚本进行转换

先检查pip库中的networkx是什么版本，输入命令pip list，在下方出现的Package找到networkx

我这里是2.5版本，最好将其降级为2.3版本，不然会出现如下提示：

输入命令pip install networkx==2.3，等待安装

出现如下界面显示安装成功

建议python3.x采用3.5.2 和 3.7.3.，不然也会出现NOTE，如下（但我也没想管了）：

输入命令：python cgenff_charmm2gmx.py JZ4 jz4_fix.mol2 jz4.str charmm36-jul2022.ff

会出现如下报错

此时只需要用vs code打开py文件，然后再151行的位置加上,encoding="utf-8"

重新输入命令，出现如下界面表示成功生成拓扑文件

此时会得到这样几个文件

创建复合物坐标文件：

输入命令：gmx editconf -f jz4_ini.pdb -o jz4.gro

采用editconf模块将pdb文件转换为gro文件，要将JZ4 (添加H后) 的坐标追加到溶菌酶的坐标上，因此要转换为gro文件，我们可以得到以下文件

将3htb_processed.gro另存为complex.gro，将jz4.gro内容粘贴到complex.gro中（注意是粘贴到蛋白向量后，盒向量前），记得最上面的原子数改为2636，因为新增了22个原子

创建复合物拓扑文件：

因为拓扑文件是分别生成的（gmx不能处理配体所以分别生成），所以我们这里还要合起来

在topol.top文件里直接修改

→→→→→→

在topol.top文件里增加配体topology

→→→→→→

在topol.top文件里增加配体参数

→→→→→→

定义盒子：

gmx editconf -f complex.gro -o newbox.gro -bt dodecahedron -d 1.0

填充溶剂：

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

我们将会得到两个文件

使用vmd查看solv.gro文件（pbc box -on），会发现溶剂并没有完全在盒子内，这个是正常的，可以用这个继续跑

添加离子：

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr

此过程会生成ions.tpr文件

gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral

此时会新生成两个文件

记得检查top文件里面有没有更新离子信息。

ps：注意，使用charmm36-jul2022.ff添加离子信息时，后面不会报错，但使用charmm36-jul2021.ff时会显示报错，no such moleculetype CL，没有这个原子类型。

这是因为在2021的力场文件中有个ions.itp文件，里面没有NA和CL的原子类型，后者说有只是定义的不一样（下图为2022的文件内容），所以如果下载的是2021的力场文件话，那就要换一种阳离子和阴离子，后面的操作都是一样的（ps：后续发现2021的CL定义为CLA，NA定义为SOD）

能量最小化：

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

gmx mdrun -v -deffnm em

在我的电脑上需要400多步才收敛，我看其他都是100多步就收敛了

限制配体：

先生成配体（非H）的索引文件

gmx make_ndx -f jz4.gro -o index_jz4.ndx

输入：2 & ! a H*

接着输入q，保存并结束，可以得到索引文件

输入命令：gmx genrestr -f jz4.gro -o posre_jz4.itp -n index_jz4.ndx -fc 1000 1000 1000

选择4，回车

此时会生成配体的位置限制文件：

在topol.top文件里增加位置限制

→→→→→→

平衡前处理（对控温重新分组）：

tc-grps不要对体系中的单一物质耦合，尽量进行分组，但是考虑到配体和蛋白结合紧密，所以分为protein non-protein不合理，所以我们使用命令：gmx make_ndx -f em.gro -o index.ndx，将蛋白和jz4分为一组，剩下一组为water_and_inors

输入1|13，回车

接着输入q保存并结束，可以得到新的索引文件（原来双相系统的这个索引是这样建立的，我还在疑惑呢，总算是学会了），这时候记得将nvt.mbp文件里的tc-grps改一下

NVT平衡：

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr

gmx mdrun -deffnm nvt

这次平衡要跑20min，我这小笔记本太辛苦了（；´д｀）ゞ

终于跑完了

NPT平衡：

gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr

gmx mdrun -deffnm npt -v

MD平衡：

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_10.tpr

gmx mdrun -deffnm md_0_10

最后等待跑完即可

后续处理再写一篇吧，md平衡实在是太长时间了