JiaoYuan's Blog
近期做了一次分子动力学模拟实验，命令比较多，整理一下步骤，参考了生信实战的教程，但他们的教程存在很多问题，这里完善一下。
以江苏大学一位师兄提供的蛋白文件为基础，前提是要去除小分子，否则 gromacs 会报错提示有非标准残基。
环境准备
使用 conda 来安装软件：
conda create -n gmx python=3.9 #创建虚拟环境
conda activate fmx
conda install -c bioconda gromacs
安装完成后会 WARING，需要安装 pocl：
conda install -c conda-forge pocl
数据准备
准备好蛋白文件，我这里的蛋白文件是YdcW_C284S.pdb。
首先去除本实验中不需要的水分子，将去除水分子的蛋白文件命名为YdcW_C284S.clean.pdb ：
grep -v HOH YdcW_C284S.pdb > YdcW_C284S.clean.pdb
生成拓扑文件，并且将拓扑文件命名为YdcW_C284S_processed.gro：
gmx pdb2gmx -f YdcW_C284S.clean.pdb -o YdcW_C284S_processed.gro -water spce
这时会提示我们进行选择，输入15选择OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)力场。
执行后记录下体系的电荷情况，例如本例：Total charge 9.000 e。
接下来需要下载一些参数文件，为后面的分析做准备：
wget http://www.mdtutorials.com/gmx/lysozyme/Files/ions.mdp
wget http://www.mdtutorials.com/gmx/lysozyme/Files/minim.mdp
wget http://www.mdtutorials.com/gmx/lysozyme/Files/nvt.mdp
wget http://www.mdtutorials.com/gmx/lysozyme/Files/npt.mdp
wget http://www.mdtutorials.com/gmx/lysozyme/Files/md.mdp
定义模拟体系的尺寸并加水
先使用 editconf 模块定义晶胞容器的尺寸，并且将处理后的文件命名为YdcW_C284S_newbox.gro：
gmx editconf -f YdcW_C284S_processed.gro -o YdcW_C284S_newbox.gro -c -d 1.0 -bt cubic
使用 solvate 调用溶剂模块向晶胞容器中填充水：
gmx solvate -cp YdcW_C284S_newbox.gro -cs spc216.gro -o YdcW_C284S_solv.gro -p topol.top
为模拟体系添加离子
生成ions.tpr文件：
gmx grompp -f ions.mdp -c YdcW_C284S_solv.gro -p topol.top -o ions.tpr
加入中和离子：
gmx genion -s ions.tpr -o YdcW_C284S_solv_ions.gro -p topol.top -pname NA -nname CL -neutral
选择 13 后回车。
体系能量最小化
生成能量最小化的拓扑文件：
gmx grompp -f minim.mdp -c YdcW_C284S_solv_ions.gro -p topol.top -o em.tpr
执行能量最小化：
gmx mdrun -v -deffnm em
这一步需要很长时间，所以建议使用下列命令挂在后台运行：
nohup gmx mdrun -v -deffnm em &
gmx energy 模块来分析能量最小化的结果：
gmx energy -f em.edr -o potential.xvg
输入 10 0 来选择势能 Potential(10), 并用零 (0) 来结束输入。
xvg 文件可以使用 DuIvyTools 来进行查看，DuIvyTools 需要用 pip 来安装：
pip config set global.index-url https://pypi.tuna.tsinghua.edu.cn/simple #设置为清华源
pip install DuIvyTools
使用下列命令查看
dit xvg_show -f potential.xvg
NVT 温度平衡
执行：
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
再执行下列命令，所需时间比较长，依旧挂在后台：
nohup gmx mdrun -deffnm nvt &
平衡结果分析：
gmx energy -f nvt.edr -o temperature.xvg
输入 16 0，选择温度代号 16，选择 0 退出。
查看温度平衡情况：
dit xvg_show -f temperature.xvg
NPT 压力平衡
执行：
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
nohup gmx mdrun -deffnm npt &
压力分析：
gmx energy -f npt.edr -o pressure.xvg
根据提示输入 18 0 获得压力的数据结果。
查看压力分析结果：
dit xvg_show -f pressure.xvg
密度分析：
gmx energy -f npt.edr -o density.xvg
根据提示 输入 24 0 获得结果。
查看压力分析结果：
dit xvg_show -f pressure.xvg
正式的动力学模拟
执行 grompp 生成模拟的 tpr 文件：
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr
执行模拟：
nohup gmx mdrun -deffnm md_0_1 & #挂在后台
分子动力学模拟结果分析
使用下列命令进行修复：
gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center
选择 1 (Select 1 (“Protein”) as the group to be centered and 0 (“System”) for output) 完成后生成 md_0_1_noPBC.xtc，我们将基于这个修复的数据进行后续的数据分析。
计算 RMSD：
gmx rms -s md_0_1.tpr -f md_0_1_noPBC.xtc -o rmsd.xvg -tu ns
选择 4 (“Backbone”) for both the least-squares fit and the group for RMSD calculation。
查看 RMSD 结果：
dit xvg_show -f rmsd.xvg
计算相对于模拟之前晶体的结构差异，可以使用下面的命令：
gmx rms -s em.tpr -f md_0_1_noPBC.xtc -o rmsd_xtal.xvg -tu ns
也是选择 4 (“Backbone”)。
计算回旋半径：
gmx gyrate -s md_0_1.tpr -f md_0_1_noPBC.xtc -o gyrate.xvg
选择 group 1 (Protein)，执行完成后新生成了 gyrate.xvg。
查看结果：
dit xvg_show -f gyrate.xvg
这一套流程都是按照默认参数和步骤来的，实际上我也不太懂，但是大多数情况下使用默认参数是足以满足实验需求。