GROMACS-常规MD教程
1） 准备同源建模得到的nsp14结构文件（nsp14_model_ions.pdb）
2） H++网站预测nsp14蛋白质子化情况：
提交nsp14_model_ions.pdb
查看result，并下载result.pdb文件
pymol中查看result.pdb文件中His质子化结果:
select HISs, resn his
Selector: selection "HISs" defined with 342 atoms.
HISs=cmd.get_model('HISs').atom
print len(set(map(lambda i: i.resi,HISs)))
20
print sum(map(lambda i: i.partial_charge,HISs))
0.0
结果表明该pdb文件中共包含20个His，这20个His均不带电荷。
3） 在pymol中删去nsp14_model_ions.pdb中的SAM配体小分子，得到nsp14_model_ions_noligand.pdb文件
4） 在服务器中新建文件夹nsp14，并将nsp14_model_ions_noligand.pdb、以及所需参数文件em.mdp和npt.mdp放入文件夹nsp14
5） echo 6 1|gmx_mpi pdb2gmx -f nsp14_model_ions_noligand.pdb -ignh -o conf.pdb
echo 1|gmx_mpi editconf -f conf.pdb -o box.gro -c -d 1.0 –princ
gmx_mpi solvate -cp box.gro -cs -o box_water.pdb -p topol.top
gmx_mpi grompp -f em.mdp -p topol.top -c box_water.pdb -o ions.tpr –maxwarn 1
echo 16| gmx_mpi genion -s ions.tpr -p topol.top -o box_ions.pdb -neutral -conc 0.1
#这一步中在“select a continuous group of solvent molecules”时选择“SOL”
Q:为什么要选择SOL，溶剂不是water吗？
因为在结构文件中对溶剂的标识为SOL
Q:离子浓度为什么设定为0.1？
#在topol.top文件中查看，添加了48个钠离子和58个氯离子
gmx_mpi grompp -f em.mdp -p topol.top -c box_ions.pdb -o em.tpr -maxwarn 1
gmx_mpi mdrun -ntomp 32 -v -deffnm em
gmx_mpi grompp -f npt.mdp -c em.gro -p topol.top -o npt.tpr -maxwarn 1
gmx_mpi mdrun -ntomp 32 -v -deffnm npt
6） 进行常规MD模拟
cp /home/jawang/chpeng/numd2/pocket/1c85/0043/md.mdp ./
#复制md.mdp参数文件到nsp14文件夹中，更改模拟时长参数如下：
dt = 0.002 (ps)
nsteps = 50000000; 3 micro-second
gmx_mpi grompp -f md.mdp -c npt.gro -p topol.top -o md.tpr -maxwarn 1
#准备md.tpr文件
nohup gmx_mpi mdrun -ntomp 32 -v -deffnm md &
#开始常规MD模拟
#nohup：与&配合，表示后台运行
MD模拟结果的后续处理学习
模拟体系的继续运行与停止：
1）可以在100 ns的MD模拟基础上，继续运行该体系的MD模拟
gmx_mpi convert-tpr -s md.tpr -f md.cpt -o mdnew.tpr -extend 40000
# convert-tpr：用于准备tpr文件，具体指延长MD模拟时间（通常该任务由于某些原因已经停止）时所需要的输入文件，即mdnew.tpr
#-s：指定已存在的tpr文件
#-f：指定轨迹文件
#.cpt：模拟断点文件（check point），该文件在模拟过程中每隔固定的时间间隔产生，保存了模拟系统的所有信息。若模拟因某些原因中断，可以使用该文件重新在断点处开始模拟，也可以从该断点文件开始，延长模拟时间
#-extend：指定延长的时间，单位ps，比如40000指40ns
nohup gmx_mpi mdrun -s mdnew.tpr -noappend -ntomp 32 -v &
#nohup：与&配合，表示后台运行
#-noappend：模拟数据按照已有顺序继续产生，不会覆盖掉已有的模拟数据
2）若要终止运行，在top界面，点击K（kill），之后输入任务的名字，即可终止该任务。
轨迹文件的处理与分析-常用命令
gmx_mpi trjconv -f md.trr -s md.tpr -pbc mol -ur compact -o look.pdb -skip 10
#trjconv：压缩trr轨迹文件，并转化为look.pdb文件
#-pbc mol -ur compact：去除周期性盒子，保持结构完整性
#-skip 10：每10个构象输出1个
#可在VMD中打开look.pdb文件，观察蛋白质的运动情况
计算完整蛋白质结构间的RMSD值
gmx_mpi rms -f md.trr -s md.tpr -o rms.xvg
#rms：用于计算RMSD值
#-f：指定轨迹文件
#-s：指定参考结构的文件
#-o：指定输出文件
#.xvg：二维图标文件，二维画图工具xmgrace的默认文件，可以使用xmgrace打开
#Select group for least squares fit：选择Protein-H，不比较H的结构
#Select group for RMSD calculation：Protein-H
计算指定残基的RMSD值：
gmx_mpi make_ndx -f npt.gro -o index.ndx
#make_index：该程序产生gromacs的索引文件，即.ndx文件。索引文件在gromacs中非常重要，可以对模拟系统进行分割，从而实现不同的实验目的
#-f：指定结构文件，同一个蛋白质的任何结构文件都可以（pdb、gro、tpr等）
#-o：指定输出文件
#index.ndx：产生的索引文件
#输入q退出该界面
vi index.ndx
#在npt.gro文件里面寻找对应残基的原子序号，然后添加到index.ndx文件中，相当于手动添加一个自定义的索引，比如residue
gmx_mpi rms -f md.trr -s md.tpr -o rms_activate.xvg -n index.ndx
#-n：指定索引文件
# Select group for least squares fit：这里选择protein和protein-H都可以
# Select group for RMSD calculation：这里选择在上一步中手动创建的索引名称,比如residue
#计算完成后得到的xvg文件中第一列为构象序号，第二列为对应的RMSD值
sort -rn -k2 rms_activate.xvg|less
#利用sort命令对第2列进行排序，查看最大值和最小值