Gromacs动力学模拟

Gromacs动力学模拟

实验要求：

实验对象：选取目标体系，可自行选择一个蛋白质体系，也可以用Modeller 建模最佳模型结构.pdb进行操作。
软件： Gromacs-5.1.2 www.gromacs.org (manual/ tutorial)

实验步骤：

1. 加力场。
   gmx pdb2gmx –h 打开帮助菜单。 选力场的时候选择 Amber99sb…,溶剂类型选Tip3p。
2. 加模拟盒子。
   gmx editconf -bt ( boxtype: 做三个盒子的对比cubic/triclinic/dodecahedron ), -d 1.0 , 比较三种类型的盒子水分子数目差别。
3. 加溶剂。
   Gmx solvate –h
4. 做能量优化。 参数文件：em.mdp（可从官网上下载） etol: 500，达到收敛，实验完成。
5. 平衡体系，将体系升温。 从0K升温到300K，在30ps内完成。
6. 动力学模拟采样。模拟时间：1ns，步长：2fs。坐标保存的频率为每100ps保存一帧结果，整个轨迹共200个frame.
7. 结果分析:
   7.1. 全体系的alpha-C原子的均方根偏差(RMSD)结果获取及分析，g_rmsd。
   7.2. 全体系的alpha-C原子的均方根涨落(RMSF)结果获取及分析, g_rmsf。
   7.3. 体系的总势能变化曲线分析，采用g_energy命令。
   7.4. 分析蛋白质的回旋半径变化，采用g_gyrate命令
   7.5. 分析溶剂的可及表面积（SASA），采用g_sas命令
   7.6. 将采样最后的构象与初始构象进行叠加比较，分析构象的变化情况。
   7.7. 从模拟的轨迹中将体系中的蛋白质单独取出来，另存为一个轨迹文件protein.xtc, 用VMD的插件”movie maker”做成一个小电影。

二、操作过程记录及结果

总体流程图

图表 1流程

1、加力场

输入如下命令，在命令行的交互式操作中，选择力场(输入5)AMBER99SB protein，选择溶剂模型(输入1)：TIP3P

gmx pdb2gmx -f 5mgh_pro.pdb -o conf.gro -p topol.top

图表 2 加力场命令

2、加盒子

输入如下命令：分别加三个模拟盒子，比较三种类型的盒子水分子数目差别。

gmx editconf -f conf.gro -bt cubic -d 1.0 -o cubic_out.gro
gmx editconf -f conf.gro -bt triclinic -d 1.0 -o triclinic_out.gro
gmx editconf -f conf.gro -bt dodecahedron -d 1.0 -o dodecahedron_out.gro

3、加溶剂

输入如下命令：

gmx solvate -cp cubic_out.gro -o cubic.pdb -p topol.top
gmx solvate -cp triclinic_out.gro -o triclinic.pdb -p topol.top
gmx solvate -cp dodecahedron_out.gro -o dodecahedron.pdb -p topol.top

立方体盒子里的水分子最多。

4、做能量优化

由于是示例而已，从官网上下载一个em.mdp的例子文件，下载地址如下
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/complex/Files/em.mdp
检测离子平衡
输入如下命令，检测体系的离子平衡：

gmx grompp -f em.mdp -c cubic.pdb -p topol.top -o cubic_em

图表 3 检测结果

添加离子
发现总共是6个负电荷，需要加6个正电荷
添加离子命令：

gmx genion -s cubic_em.tpr -o ION.gro -p topol.top -np 6 -pname NA

再次检测离子平衡
再次使用新的文件ION.gro进行检测：

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

图表 4 再次检测

能量最小化

gmx mdrun -s cubic_em.tpr -deffnm cold –v

能量查看(potential)

gmx energy -f cold.edr -o cold.xvg

图表 5 低温能量
将cold.xvg导入WPS绘制表格，确实能量在降低

5、体系升温

根据生成的mdout.mdp文件，创建upgrade.mdp，将体系升温至300k，步长为2fs
在这里插入图片描述

图表 6 升温参数文件

体系升温

gmx grompp -f upgrade.mdp -c cold.gro -p topol.top -o cubic_em_hot.tpr
gmx mdrun -s cubic_em_hot.tpr -v -deffnm hot

能量查看(temperature)

gmx energy -f hot.edr -o hot.xvg 

图表 7 高温能量
将hot.xvg导入WPS绘图，发现确实升温到了300k

6、采样

体系采样

gmx grompp -f sample.mdp -c hot.gro -p topol.top -o cubic_em_hot_sample.tpr
gmx mdrun -s cubic_em_hot_sample.tpr -v -deffnm hot_sample

能量查看(temperature)

gmx energy -f hot_sample.edr -o hot_sample.xvg 

总能量确实达到了300K

图表 8 采样能量

结果分析与讨论

全体系的alpha-C原子的均方根偏差(RMSD)结果获取及分析

gmx rms -s cubic_em_hot_sample.tpr -f hot_sample.trr -o rmsd.xvg

图表 10均方根偏差

全体系的alpha-C原子的均方根涨落(RMSF)结果获取及分析

gmx rmsf -s cubic_em_hot_sample.tpr -f hot.xtc -o rmsf.xvg

图表 11均方根涨落

体系的总势能变化曲线分析

gmx energy -f hot_sample.edr -o potential.xvg

由于采样只选了0.1纳秒，不足1纳秒，所以总势能变化看上去有一点点像在上升，其实如果继续模拟，应该是在波动而已。

图表 12总势能变化

分析蛋白质的回旋半径变化(只绘制了半径的曲线图)

gmx gyrate -f hot_sample.trr -s cubic_em_hot_sample.tpr -o gyrate.xvg

图表 13回旋半径

分析溶剂的可及表面积（SASA）

gmx sasa -f hot_sample.trr -s cubic_em_hot_sample.tpr -o area.xvg

图表 14表面积

将采样最后的构象与初始构象进行叠加比较，分析构象的变化情况。RMSD不是很高，变化不太大。

gmx confrms -f1 ION.gro -f2 hot_sample.gro -o fit.pdb

图表 15构象叠加比较

VMD小电影制作

从模拟的轨迹中将体系中的蛋白质单独取出来，另存为一个轨迹文件protein.xtc, 用VMD的插件“movie maker”做成一个小电影。(以下命令选取其一即可)

gmx trjconv -s cubic_em_hot_sample.tpr -f hot_sample.trr -o pengpeng.pdb
gmx trjconv -s cubic_em_hot_sample.tpr -f hot_sample.trr -o pengpeng.trr

xtc文件是通过trjconv命令将轨迹文件.trr（是轨迹文件，包含坐标，速度，力等）压缩成的
百度搜索后下载VideoMach，将pdb或者trr加载进VMD中，
选择Extension->visualization->movie maker，Movie duration(seconds)设置生成视频的总时长，此处建议生成大约10秒左右，（之前nstxout指定过输出频率）

附录

Gromacs相关命令

gmx pdb2gmx -f 5mgh_pro.pdb -o conf.gro -p topol.top

gmx editconf -f conf.gro -bt cubic -d 1.0 -o cubic_out.gro
gmx editconf -f conf.gro -bt triclinic -d 1.0 -o triclinic_out.gro
gmx editconf -f conf.gro -bt dodecahedron -d 1.0 -o dodecahedron_out.gro

gmx solvate -cp cubic_out.gro -o cubic.pdb -p topol.top
gmx solvate -cp triclinic_out.gro -o triclinic.pdb -p topol.top
gmx solvate -cp dodecahedron_out.gro -o dodecahedron.pdb -p topol.top

wget http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/complex/Files/em.mdp

gmx grompp -f em.mdp -c cubic.pdb -p topol.top -o cubic_em

gmx genion -s cubic_em.tpr -o ION.gro -p topol.top -np 6 -pname NA

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

gmx mdrun -s cubic_em.tpr -deffnm cold -v

gmx energy -f cold.edr -o cold.xvg

依据mdout.mdp参数，生成upgrade.mdp
gmx grompp -f upgrade.mdp -c cold.gro -p topol.top -o cubic_em_hot.tpr

gmx mdrun -s cubic_em_hot.tpr -v -deffnm hot

gmx energy -f hot.edr -o hot.xvg

gmx grompp -f sample.mdp -c hot.gro -p topol.top -o cubic_em_hot_sample.tpr

gmx mdrun -s cubic_em_hot_sample.tpr -v -deffnm hot_sample

gmx energy -f hot_sample.edr -o hot_sample.xvg 

gmx rms -s cubic_em_hot_sample.tpr -f hot_sample.trr -o rmsd.xvg
gmx rmsf -s cubic_em_hot_sample.tpr -f hot.xtc -o rmsf.xvg（选择C-alpha）
gmx energy -f hot_sample.edr -o potential.xvg
gmx gyrate -f hot_sample.trr -s cubic_em_hot_sample.tpr -o gyrate.xvg
gmx sasa -f hot_sample.trr -s cubic_em_hot_sample.tpr -o area.xvg

gmx confrms -f1 ION.gro -f2 hot_sample.gro -o fit.pdb

gmx trjconv -s cubic_em_hot_sample.tpr -f hot_sample.trr -o pengpeng.pdb

gmx trjconv -s cubic_em_hot_sample.tpr -f hot_sample.trr -o pengpeng.trr

mdp文件重要参数

integrator 动力学模拟方法，即整合牛顿力学定理的方法
md： 使用跳蛙算法(leap-frog)整合牛顿定律。
steep：使用最速下降法进行能量优化。
(能量优化最大位置移动用emstep设定，能量最大容忍度由emtol决定)
emtol：最大容许力。默认为10.0。当最大作用力小于此值，认为最小化过程收敛。

nsteps：最大模拟步数。
dt :时间步长(fs)
nstxout:坐标保存的频率(每多少步保存一帧)

PS: nsteps× dt=模拟时间总长
注意设置nstxout采样保存