跳转至

分子动力学模拟流程

分子动力学模拟备忘录

分子对接与动力学模拟

对小分子和蛋白质进行分子对接。生成复合物pdb文件；或直接采用已经分子对接完成的complex文件

  1. 对拼接文件进行拆分，并修改小分子名称，拆分出蛋白质和小分子两个文件

  2. 把小分子加氢，保存为mol2文件

  3. 准备文件 可在GROMACS官方教程网站下载

charmm36-jul2022.ff
ions.mdp
em.mdp
nvt.mdp
md.mdp
npt.mdp
sort_mol2_bonds.pl
cgenff_charmm2gmx_py_3_nx2,py
  1. 制备蛋白拓扑文件
gmx pdb2gmx -f protein.pdb -o protein_processed.gro -ignh -ter
#-ignh: 这是命令的一个选项。GROMACS在处理PDB文件时，忽略任何可能存在的氢原子，然后它会根据你的力场和拓扑信息自动添加所有缺少的氢原子。官网无此选项。
#此教程删除金属原子铁，因为所选力场文件暂时不识别铁离子
#依次选择1，1，0，0
  1. 制备配体拓扑文件
小分加氢，保存为mol2图格式，运行以下脚本
perl sort_mol2_bonds.pl unk.mol2 unk_fix.mol2
在网站https://cgenff.com/  处理unk_fix.mol2，获得后缀str的文件
python cgenff_charmm2gmx_py3_nx2.py UNK unk_fix.mol2 unk_fix.str charmm36-jul2022.ff 
gmx editconf -f unk_ini.pdb -o unk.gro
  1. 新建复合物拓扑文件
新建一个complex.gro文件，将蛋白质gro文件保留全部信息和unk.gro（仅保留原子信息）的内容复制到其中，并修改原子总数（加上小分子)
在toop文件中插入以下内容 插入位置
——————————————————————————
; Include forcefield parameters
#include "./charmm36-jul2022.ff/forcefield.itp"

; Include ligand parameters
#include "unk.prm"

[ moleculetype ]
——————————————————————————
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

; Include ligand topology
#include "unk.itp"`

; Include topology for ions
#include "./charmm36-jul2022.ff/ions.itp"
———————————————————————————
在文件最后
加入UNK       1
  1. 溶剂化
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 
#检查topol文件，查看文件末尾SOL是否和UNK在同一行，如果在，需要在SOL之前按回车，到下一行
  1. 添加离子
gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname SOD -nname CLA -neutral
选择15
  1. 能量最小化
gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr #若报错 添加-maxwarn 1
gmx mdrun -v -deffnm em #可以修改为 gmx mdrun -v -deffnm em -ntmpi 1 -ntomp 0 #cpu调用错误

  1. 约束配体

    gmx make_ndx -f unk.gro -o index_unk.ndx
 > 0 & ! a H*
 > q

gmx genrestr -f unk.gro -n index_unk.ndx -o posre_unk.itp -fc 1000 1000 1000
3

#在topol文件中插入以下内容，位于配体分子下面

; Ligand position restraints
#ifdef POSRES
#include "posre_unk.itp"
#endif

gmx make_ndx -f em.gro -o index.ndx
1 | 13
14 | 15
q

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr
#该步骤需要修改tc-grps行

  2. NVT平衡

    gmx mdrun -v -deffnm nvt

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

  3. NPT平衡

    gmx mdrun -v -deffnm npt

  4. 预模拟仿真，检查文件

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

  5. 模拟分析

    gmx mdrun -v -deffnm md100

结果分析绘图

#对轨迹进行周期性矫正
gmx trjconv -s md100.tpr -f md100.xtc -o md100_nojump.xtc -pbc nojump -ur compact
0

#导出rmsd数据
gmx rms -s md100.tpr -f md100_nojump.xtc -o rmsd.xvg -tu ns
1
1 
建议选择4和4

#导出ligandrmsd数据
gmx rms -s md100.tpr -f md100_nojump.xtc -o rmsdligand.xvg -tu ns
13
13
建议选择4和13

导出蛋白配体复合物rmsd
#创建新的索引文件
gmx make_ndx -f md100.tpr -o index_complex.ndx
1｜13
name 18 Protein_LIG
q
#计算rmsd
gmx rms -s md100.tpr -f md100_nojump.xtc -n index_complex.ndx -o rmsd_complex.xvg -tu ns
4
18

#导出轨迹文件-特定帧
gmx trjconv -s md100.tpr -f md100_nojump.xtc -o start.pdb -dump 0 -pbc mol -ur compact

********************************
#快速导出全程轨迹文件xtc
gmx make_ndx -f md100.tpr -o index.ndx
1|13
name 18 Protein_UNK

gmx trjconv -s md100.tpr -f md100.xtc -n index.ndx -o protein_ligand_structure.pdb -dump 0

gmx trjconv -s md100.tpr -f md100.xtc -n index.ndx -o md_protein_ligand_skip10.xtc -pbc whole -skip 10

进入目录 cd 文件夹路径

#清空状态
reinitialize

#加载文件
load protein_ligand_structure.pdb, complex
load_traj md_protein_ligand_skip10.xtc, complex

#设置样式
# --- 隐藏所有，从头开始设置 ---
hide everything, all

# --- 将蛋白质显示为灰色的卡通模式 ---
show cartoon, polymer
color gray70, polymer

# --- 将配体显示为彩色的“棍状”模型 ---
show sticks, resn UNK
util.cbam resn UNK

# --- 将视角聚焦于配体，方便观察 ---
zoom resn UNK
center resn UNK


#调节播放速度
set movie_fps, 5

#导出gyrate
gmx gyrate -s md100.tpr -f md100_nojump.xtc -o gyrate.xvg
1
1

#导出氢键
gmx hbond-legacy -f md100_nojump.xtc -s md100.tpr -num hbond_num.xvg -life hbond_life.xvg -dist hb  
##修改为  
1
13

#导出蛋白rmsf
gmx rmsf -s md100.tpr -f md100_nojump.xtc -o fws-rmsf.xvg -ox fws-avg.pdb -res -oq fws-bfac.pdb

#导出结构快照，末尾的0代表几ns
gmx trjconv -s md100.tpr -f md100_nojump.xtc -o start.pdb -dump 0


#绘制自由能形貌图
#sham.xvg为回旋半径和RMSD的数据文件的组合，第一列为时间，第二列第三列为相应时间下的回旋半径和RMSD数据
输入rmsd和回旋半径文件，利于脚本生成sham文件

gmx sham -f sham.xvg -ls FEL_sham.xpm -b 8500 -tsham 310 -nlevels 100（绘制3d图像时 为了适配转换脚本 需要设置为25）

#利用DuIvyTools绘图
dit xpm_show -f FEL_sham.xpm -3d -o FEL_sham.pdf
dit xpm_show -f FEL_sham.xpm -o FEL.pdf -cmap jet
dit xpm_show -f FEL_sham.xpm -m 3d -o FEL_3D.pdf -cmap jet

#利用B站UP主的脚本绘图 更美观
脚本文件在压缩包内