Molecular dynamics of Protein-Ligand by Gromacs
一、关于GROMACS 安装及运行
Gromacs适合在NVIDIA GPU运行? AMD的APU（5700g）装不上显卡版Gromacs
安装目录避免空格/中文
具体参照“基于Ubuntu22.04的Linux系统安装GROMACS”
https://zhuanlan.zhihu.com/p/628218852
二、动力学模拟部分
（以下参照“蛋白质配体复合物-分子动力学模拟Gromacs”https://blog.csdn.net/chaojFF/article/details/123112216）
1. 准备蛋白质结构文件、分子对接后的配体结构文件
1.1首先得到分子对接后的蛋白质结构文件与配体结构文件（通常得分最高的构象）
1.2把蛋白质和配体分子的坐标文件分开保存。配体的坐标可以通过如下命令保存：
grep Protein Ligand.pdb > Ligand.pdb
然后从Protein Ligand.pdb文件中删除与Ligand有关的几行.
#或者通过pymol 手动分离，选中配体后选择摘出，然后分别保存蛋白和配体（如果直接选中配体保存，可能后续盒子中出现蛋白配体分离现象）
2. 准备拓扑、选择力场
2.1蛋白拓扑
gmx pdb2gmx -f Protein.pdb -o Protein_processed.gro -water spc
gmx会提示选择一个力场，输入6，选择Amber99sb-ildn力场
# 如果出现类似报错：“Fatal error:Atom HD1 in residue HIS 76 was not found in rtp entry HISB with 12 atomswhile sorting atoms.”，说明pdb文件中对氢原子的命名与力场rtp文件中命名规则与顺序不同。解决方法：①如果需要保留初始氢原子，需要找到该力场的rtp文件查看文件中的命名规则对氢原子重命名。②在上面的命令行中末尾加上 -ignh 命令，用此选项忽略文件中的氢原子, 统一使用氢原子命名规则, 以免出现氢原子名称不一致的情况。
选择完力场得到三个文件: Protein_processed.gro, topol.top, posre.itp
2.2小分子拓扑
使用ACPYPE在线版处理配体，进入网站Bio2Byte home page, 上传小分子pdb文件，全部选项默认，提交后等待十分钟以内，获取压缩文件
提取压缩文件中Ligand_GMX.itp, Ligand_GMX.igro 放到工作目录。
# 小分子文件用文本编辑器打开，如果文件内的文件名等信息有错，需要手动修改
# 如果生成失败，更改小分子配体摘取方法
3. 创建复合物拓扑
3.1将之前处理得到的蛋白质文件（Protein_processed.gro）复制并改名为complex.gro：cp Protein_processed.gro complex.gro
3.2接下来复制Ligand_GMX.gro的坐标部分粘贴到complex.gro文件中蛋白质原子最后一行的下面(也就是文件末尾一行的上方)。
由于添加了配体原子到.gro文件中, 将complex.gro文件中第二行的数字增加配体原子数量。
3.3在topol.top文件中插入一行，内容为
; Include ligand topology
#include "Ligang_GMX.itp"
如图：
3.4最后在文档的末尾修改[ molecules ]部分，添加PNP信息
（PNP 改为Ligand）
#后面如果运行出错，注意这里是否写错了
4. 添加盒子、溶剂
4.1定义单元盒子gmx editconf -f complex.gro -o newbox.gro -bt dodecahedron -d 1.0
得到newbox.gro
# -d决定了盒子的尺寸，即盒子边缘距离分子边缘 1.0nm （10Å）。理论上在绝大多数系统中，-d 都不能小于0.85nm。
# editconf 也可以用来进行gromacs文件（*.gro）和pdb 文件（*.pdb）的相互转化。例如：editconf –f newbox gro –o file.pdb 则将newbox. gro 转换为 file.pdb
# 这里检查一下蛋白在不在盒子里，蛋白配体位置是否与对接结果原位置一致。
4.2填充溶剂水 gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro
得到solv.gro文件
5. 添加离子
5.1 点击下载em.mdp文件，存放到同工程目录下(下载蛋白配体的那个文件夹里的)。
然后输入命令
gmx grompp -f em.mdp -c solv.gro -p topol.top -o next.tpr
如果出现平衡电荷数警告，在命令后端加上 -maxwarn 1
gmx grompp -f em.mdp -c solv.gro -p topol.top -o next.tpr -maxwarn 1
5.2翻看上面的输出信息，可以看到整个体系携带的电荷数。
由于生命体系中不存在净电荷，必须在体系中添加离子，电荷为负添加钠离子，电荷为正添加氯离子。（-pname 阳离子的名称 -nname 阴离子的名称 -np 阳离子个数 -nn 阴离子个数）
gmx genion -s next.tpr -o solv_ions.gro -p topol.top -pname NA -np 19
(19就是添加的离子数)
运行上面的命令后会提示选择溶剂，选择15 SOL；得到solv_ions.gro文件。
#根据前面选择的力场，有的力场文件里氯离子/钠离子表示方式可能不同，需要去力场文件里查看。（一般为NA / CL）
6. 能量最小化
gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr
得到em.tpr文件，运行EM：gmx mdrun -v -deffnm em
得到四个文件.edr .trr .log .gro
# em.mdp也是从上述网站下载
7. NVT、NPT平衡
7.1下载nvt.mdp npt.mdp文件，存放到同工程目录下。
平衡蛋白质配体复合物与平衡任何其他蛋白质水溶液体系一样. 但需要有些特殊考虑:
a. 对配体施加限制gmx genrestr -f Ligand_GMX.gro -o posre_Ligand.itp -fc 1000 1000 1000
在topol.top文件中添加如下几行
; Ligand position restraints
#ifdef POSRES
#include "posre_Ligand.itp"
#endif
b.处理温度耦合组 热浴
对只有几个原子的组在控制其动能的涨落时, 温度耦合算法不够稳定，不要对体系中的每一单个物种使用独立的耦合。创建一个特殊的组, 其中包含蛋白质和配体。
gmx make_ndx -f em.gro -o index.ndx
通过下面的命令合并"Protein"和"PNP"，其中">"是make_ndx的提示符：
> 1 | 13 回车
> q 回车退出
7.2开始nvt平衡
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr -n index.ndx
运行之后得到nvt.tpr文件：gmx mdrun -deffnm nvt
# nvt画图：gmx energy -f nvt.edr -o temperature.xvg
7.3接下来开始npt平衡
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr -n index.ndx
运行之后得到nvt.tpr文件gmx mdrun -deffnm npt
# npt画图：gmx energy -f npt.edr -o pressure.xvg
#两步平衡后再进行模拟会节约时间
8. 开始模拟; 依次运行两条命令
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr
gmx mdrun -deffnm md_0_1 -v (-v 可显示计算结束时间)
#如果整个过程有bug, 查看topol文件中是否蛋白/配体信息需要修改？有些时候小分子文件需要写成UNK
#中断后运行：MD作业中断后可以继续用 mdrun 接着续跑，只需要加上 -cpi md.cpt 和 -s md.tpr 的参数。gmx mdrun -s md.tpr -cpi md.cpt -deffnm md
mdrun 默认会将新产生的轨迹添加到原始文件末尾，最终文件会包括中断前与续跑后的所有内容。
#续跑：续跑10ns且续写入源文件
gmx convert-tpr -s md.tpr -extend 10000 -o md.tpr
gmx mdrun -V -deffnm md -cpi md1.cpt
三、分析结果
参考教程：
https://www.bilibili.com/video/BV1Xh411p7ay/?spm_id_from=333.999.0.0
https://www.bilibili.com/video/BV1M8411z7Vm/?spm_id_from=333.999.0.0
https://blog.csdn.net/m0_55294055/article/details/119169553
# xpm转eps:
gmx xpm2ps -f FEL_sham.xpm -o FEL_sham.eps -rainbow blue
gmx xpm2ps -f enthalpy.xpm -o enthalpy.eps -rainbow blue
gmx xpm2ps -f entropy.xpm -o entropy.eps -rainbow blue
gmx xpm2ps -f prob.xpm -o prob.eps -rainbow blue
#报错参考文档
