前言
看了几篇文献后,先对MD模拟进行了个简单的实操,对水中溶菌酶进行分子动力学模拟,pdb
文件已经通过Pymol
进行了处理,删除了水分子,对这个实操做了笔记。
操作过程
产生拓扑文件
使用pdb2gmx
模块产生拓扑文件。 在WSL中打开windows系统中的文件系统:
cd /mnt
输入命令:
gmx pdb2gmx -f 1aki_clean.pdb -o 1aki_processed.gro -water spce -ignh
出现立场选择提示:
选择15(现在不懂)。 输出结果:工作目录产生三个新的文件: 1aki_processed.gro
:Gromacs的结构文件格式,包含力场中定义的所有原子 topol.top
:体系的拓扑文件 posre.itp
:包含了用于限制重原子位置的信息
定义盒子并溶剂化
模拟一个简单的蛋白水溶液体系,可在适当的情况下模拟其他不同溶剂中的蛋白质或其他分子,只要有合适的力场参数即可。 首先使用editconf
模块定义盒子的大小,输入以下命令:
gmx editconf -f 1aki_processed.gro -o 1aki_newbox.gro -c -d 1.0 -bt cubic
使用solvate
模块给盒子中填充溶剂(水),输入以下命令:
gmx solvate -cp 1aki_newbox.gro -cs spc216.gro -o 1aki_solv.gro -p topol.top
其中spc216.gro
为溶剂构型文件,为Gromacs自带的。 可视化软件Pymol查看溶剂化体系:
添加抗衡离子
由pdb2gmx
模块的输出结果,可以看到蛋白质带有+8e的净电荷。分定动力学模拟通常是在电中性下进行的,因此需要加入抗衡离子。
ions.mdp
,内容如下:(文件内容含义待学习)
; ions.mdp - used as input into grompp to generate ions.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = cutoff ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; Short-range electrostatic cut-off
rvdw = 1.0 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
使用grompp
模块生成.tpr
文件:
gmx grompp -f ions.mdp -c 1aki_solv.gro -p topol.top -o ions.tpr
产生了.tpr
类型的未见,提供了体系的原子级别的描述,包含了模拟参数说明。 之后使用genion
模块添加抗衡离子,将一些水分子替换为离子,输入命令:
gmx genion -s ions.tpr -o 1aki_solv.gro -p topol.top -pname NA -nname CL -neutral
选择溶液13
能量最小化
分子动力学模拟之前,要保证体系结构正常,几何构型合理,原子无冲撞,此时需要先将溶剂化体系弛豫至能量较低的状态,即能量最小化。 需要一个参数文件minim.mdp
,内容如下:
; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; Short-range electrostatic cut-off
rvdw = 1.0 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
使用grompp
模块将结构、拓扑和模拟参数整合到.tpr
文件中,输入命令:
gmx grompp -f minim.mdp -c 1aki_solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em
得到四个文件: em.log
:EM过程log文件 em.edr
:二进制能量文件 em.trr
:二进制全精度坐标文件 em.gro
:能量最小化后的解构 有两个因素评价EM是否成功: 第一个是势能:Epot为负,丙炔对于水中的蛋白质在10^5~10^6量级,具体取决于系统大小和水分子数量 第二个是Fmax,其目标在minim.mdp
中设置,emtol=1000.0
,意思是目标Famx的值要低于1000KJ/(mol nm)。 分析em.edr
文件中包含了Gromacs在EM期间收集的所有能量项。采用energy
模块:
gmx energy -f em.edr -o potential.xvg
生成图如下: 这是系统处于能量最小值,开始模拟。
平衡
EM确保我们在构象和溶剂取向方面具有合理的起始结构。在开始正式的动力学模拟之前,必须平衡蛋白质周围的溶剂和离子。如果这时候直接进行都力学模拟,体系可能会崩溃。原因是溶剂一般都在自身内部优化结构,而不一定需要跟溶质一起优化。需要将溶剂调整到模拟所需的温度然后在溶质周围形成温度的构象。在达到准确温度之后(基于动能),我们会给体系施加压力指定体系密度达到合理值。
平衡过程锋味两个阶段,第一阶段实在NVT系综(正则系综,等温等容系综)下进行。其中体系的粒子数(N),体积(V)和平均温度(T)保持不变。如果温度尚未稳定,则需要额外的时间,一般50-100ps.进行100ps的NVT平衡。需要一个nvt.mdp
参数文件。(Index of /gmx/lysozyme/Files (mdtutorials.com)) 输入以下命令执行NVT平衡:
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt
第一条命令运行结果如下:
第二条命令执行需要时间,结果如下:分析温度变化过程,使用energy
模块:
gmx energy -f nvt.edr -o temperature.xvg
输入16 0
已选择体系温度并退出。
npt.mdp
参数文件,其中增加了压力耦合选项,使用的是Parrinello-Rahman
压浴,其他一些更改:
- continuation=yes:我们正在从NVT平衡阶段继续仿真
- gen_vel=no:不生成速度,而是从轨迹文件中读取 跟NVT平衡模拟一样,还是使用
grompp
和mdrun
。
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tp
gmx mdrun -deffnm npt
分析压力变化过程:
gmx energy -f npt.edr -o pressure.xvg
选择18 0
选择体系压力并退出,作图如下:
energy
模块输入24 0
:
gmx energy -f npt.edr -o density.xvg
作图如下:
上述两图说明压力和体系达到了很多的平衡。MD模拟
完成两步平衡后,体系处于温度和压力下的平衡状态,现在我们将位置限制取消,进行MD模拟。 运行一个1ns的MD模拟,参数文件为md.mdp
:
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr
执行mdrun
:
gmx mdrun -deffnm md_0_1
分析
对体系进行分析。介绍一些工具: 第一个是trjconv
,该后处理工具可以提取坐标、纠正周期性或者手动调整轨迹(时间单位、帧频等)。在本例中,将使用trjconv
消除体系的周期性。在模拟中,蛋白质在扩散过程中可能会穿过盒子边界,这可能会导致蛋白质的一部分出现在盒子的另一侧,考虑到这种行为,使用下面的命令:
gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center
居中的组为1(Protein)
,输出的组为0(System)
。对修正的轨迹进行所有的分析。
Gromacs
有计算rmsd
计算的内置使用程序,成为rms
。输入命令:
gmx rms -s md_0_1.tpr -f md_0_1_noPBC.xtc -o rmsd.xvg -tu ns
计算结晶结构的RMSD
,输入以下命令:
gmx rms -s em.tpr -f md_0_1_noPBC.xtc -o rmsd_xtal.xvg -tu ns
两个结果放在一起:
表明结构非常稳定,图之间的细微差异表明,t=0ns时的结构与该晶体结构略有不同。 蛋白质的均方回转半径(RG)可以描述其紧密程度。如果蛋白质处于稳固的折叠状态,Rg会处于有一个相对稳定值。如果蛋白没有折叠,Rg会一直变化。分析Rg:gmx gyrate -s md_0_1.tpr -f md_0_1_noPBC.xtc -o gyrate.xvg
选择1(Protein)
进行分析:
总结
经过实操,对MD模拟有了基本的认识,也深感自己的不足,需要继续看文献以及阅读Gromacs
手册,学习软件的使用方法以及丰富自己的理论知识。 其次,想要利用Pymol
做一个模拟后的轨迹动画,还在学习,对软件的熟练使用仍需加强。
PS:Pymol是世界上最烂的软件
参考
[1] GROMACS Tutorial