GROMACS使用教程

GROMACS使用教程 GROMACS教程 一 Gromacs基本模拟流程 .......................................................................................... 3 1 下载pdb文件..................................................................................................... 3 2 用pdb2gmx 处理 pdb 文件 ............................................................................ 3 3 建立盒子............................................................................................................. 3 5 设置能量最小化................................................................................................. 4 6 用grompp程序进行文件处理 .......................................................................... 6 7 使用 genion 和tpr文件添加离子 ................................................................... 6 8 用fws_ion.pdb来产生能量最小化的输入文件............................................... 6 7 9 在后台运行能量最小化(在命令后加&) ..................................................... 二 设置位置限制性动力学模拟.................................................................................. 7 三 设置非限制性动力学模拟...................................................................................... 9 1 如何重启一个计算........................................................................................... 11 2 如何延长一个计算........................................................................................... 11 3 如何设置并行计算........................................................................................... 11 五 模拟结果分析........................................................................................................ 12 1 如何将特定帧的轨迹保存成*.pdb文件......................................................... 12 2 用ngmx观察轨迹文件(也可以用VMD观察轨迹文件) ........................... 12 3 比较常用的分析工具....................................................................................... 14 3.3 g_covar 计算斜方差............................................................................... 16 3.4 g_energy 能量数据作图，如压力、体积、密度等 ........................... 16 3.5 g_gyrate 测量回旋半径........................................................................ 17 3.6 g_rms 与 g_rmsdist 计算结构的RMSD 值 ................................... 17 3.7 g_rmsf 计算原子位置的根均方波动( rmsf ) ................................. 18 3.8 do_dssp 计算模型的二级结构 .............................................................. 20 3.9 g_hbond 计算模拟过程中分子间的氢键的数目、距离或角度........ 21 3.10 g_saltbr 分析模拟中残基间的盐桥..................................................... 21 1 GROMACS 是一个使用经典分子动力学理论研究蛋白质动力学的高端的高效的工具。GROMACS是遵守GNU许可的免费软件，可以从以下站点下载:，并且可以在linux和 Windows上使用。 在本教程中，将研究一个从漏斗形蜘蛛的毒液中分离的毒素。我们将使用显性溶剂动力学的方法来进行研究。首先比较真空中和溶解的模型。我们将把毒素肽溶在水盒子里，紧接着用牛顿运动定律加以平衡。我们还将比较偿离子在显性溶剂动力学中的影响。 更全面的用法指导请参考官网的GROMACS 用户手册 注意:在本教程中，将要生成的gromacs(*.gro)结构文件，可以用VMD(下)查看。 2 一 Gromacs基本模拟流程 1 下载pdb文件 1OMB.pdb () 2 用pdb2gmx 处理 pdb 文件 pdb2gmx –ignh –ff G43a1 –f 1OMB.pdb –o fws.pdb –p fws.top –water spce pdb2gmx 此命令将pdb文件转换成gromacs文件并产生拓扑文件。 -ignh 因为本pdb文件是由 NMR产生的，含有氢原子，因此用-ignh选项 忽略文件中的氢原子。 -ff 指定力场(G43a1是Gromos96力场，一个通用原子力场)。 -f 读入pdb文件， -o 指定一个新产生的pdb文件(也可以是其它多种类型文件)的文件名。 指定新产生的拓扑文件名。拓扑文件包含了所有力场参数(基于一开-p 始选择的力场)，因此非常重要。 -water 来指定水模型研究表明SPC/E 水模型在水盒子模拟中表现最好。用 SPC/E 水模型研究长程静电相互作用较好。 #注:对于下面将要用到的任何命令，都可以使用“-h”查看该命令的使用方法，比如，对于命令pdb2gmx 可以使用: pdb2gmx –h 3 建立盒子 editconf -bt cubic –f fws.pdb –o fws.pdb –d 0.9 用上面的命令建立了一个简单的立方体盒子. -d 决定了盒子的尺寸，即盒子边缘距离分子边缘 0.9nm (9Å)。理论上在绝大多数系统中，-d 都不能小于0.85nm。 注:editconf 也可以用来进行gromacs文件(*.gro)和pdb 文件(*.pdb)的相互转化。 例如:editconf –f file.gro –o file.pdb 则将file.gro 转换为 file.pdb 现在就可以用产生的文件进行真空模拟了。真空模拟就是先能量最小化，然后进行动态模拟。 4 在盒子中放入溶剂 3 genbox –cp fws.pdb –cs spc216.gro –o fws_b4em.pdb –p fws.top genbox命令在editconf产生的盒子基础上生成水盒子。上面的命令行指定了SPC水盒子。 genbox命令可以在给定尺寸的盒子中加入正确数目的水分子。 5 设置能量最小化 em.mdp文件:Gromacs用*.mdp 文件指定所有计算的参数。 它用最速下降法消除原子位置碰撞。编辑文件，将 nsteps 变成400。如果最小化不能收敛，就用nsteps=500 再做一次。(最小化在400步内一般是能收敛的，但不同的平台可能结果会不一样。)要重做的话，必须重新运行grompp(注意:预处理器的位置在你的机器上可能不同，用which命令来定位，即 which cpp) em.mdp文件内容: 4 title –标题随便取(最长64 个字，简单点好) cpp –指定预处理器的位置 define –传递给预处理器的一些定义。–DFLEXIBLE 告诉 grompp将spc水 模型而非刚性 SPC包含进拓扑文件，以便用最陡下降法进一步最 小化能量。 constraints –设置模型约束 integrator – steep，告诉 gompp 使用最速下降法进行能量最小化。cg则代表使 用共轭梯度法。 dt –能量最小化用不用。只在动力学模拟中用(如 md)。 nsteps –在能量最小化中，指定最大运行步数。 nstlist –更新邻居列表的频率。 nstlist = 10表示每10 步更新一次。 rlist –短程邻居列表的阈值。 coulombtype –告诉gromacs如何计算静电。PME为particle mesh ewald 法(参 见 Gromacs用户手册) rcoulomb –指定库仑力阈值 vdwtype –告诉Gromacs如何计算范德华作用(cut-off, Shift 等) rvdw –指定 LJ 或Buckingham势能距离阈值 EM Stuff emtol –最大的力如果小于此值则能量最小化收敛(结束)(单位kJ mol–1 nm–1) emstep –初始步长(nm) 5 6 用grompp程序进行文件处理 grompp是预处理程序(即the gromacs pre-processor 的缩写) grompp –f em.mdp –c fws_b4em.pdb –p fws.top –o fws_em.tpr -f 标签指定输入参数文件(*.mdp )。 -c 输入结构文件(pdb文件，*.pdb); -p 输入拓扑文件 -o 输出mdrun的输入文件(*.tpr )。 7 使用 genion 和tpr文件添加离子 对生成的 tpr文件加入补偿离子以中和系统中的净电荷。我们的模型中有+ 2.00静电，因此加入两个氯离子。将fws_em.tpr 文件拷贝到“ionwet ”子目录，并且将fws.top 和posre.itp拷贝到这个目录。用genion 命令添加氯离子: o fws_ion.pdb –nname CL- –nn 2 –g fws_ion.log genion –s fws_em.tpr – -nname 指定阴离子名称(在Gromos G43a1力场中，用CL-表示氯离子。参 见ions.itp 查看wrt力场中离子详细信息) -nn 是指定加入的阴离子数目。 -g 输出genion 的log文件。 运行这个命令时，提示提供一个连续的溶剂组，应该是组12(SOL )。输入12 ，回车。程序会告知你有两个溶剂分子被氯离子代替。现在你必须修改fws.top 文件: 添加 #include “ions.itp” (注意:3.2及以后版本会自动添加) 经过包含声明后，力场在最后减掉两分子SOL ，加入两分子Cl 。 8 用fws_ion.pdb来产生能量最小化的输入文件 你还需要修改pr_md.mdp和md.mdp两个文件中的温度耦合参数。 加氯离子后的 pr_md.mdp 和md.mdp文件的温度耦合参数 ; Berendsen temperature coupling using vel rescaling is on Tcoupl = v-rescale tau_t = 0.1 0.1 tc_grps = protein non-protein ref_t = 300 300 6 记住:如果要加入氯离子，需要重新运行第6步的grompp。首先删除旧的 fws_em.tpr 文件，然后运行下面的 grompp命令: grompp –f em.mdp –c fws_ion.pdb –p fws.top –o fws_em.tpr 9 在后台运行能量最小化 nohup mdrun –v –s fws_em.tpr –o fws_em.trr –c fws_b4pr.pdb –e em.edr –g em.log & nohup...& 使任务后台运行 用tail命令检查最小化的进程 tail –15 em.log 当能量最小化结束，你将看到log文件中有如下总结文字，表明最速下降收敛了。 用tail -50 em.log : 二 设置位置限制性动力学模拟 什么是位置限制性模拟,你限制(或部分冻结)大分子中的原子位置，而允许溶剂分子运动。这样做像是将水分子浸入大分子。水分子松弛时间约为 10ps。因此我们要进行超过 10ps 的位置限制性模拟。本实例中用 20ps，大的模型(大蛋白或脂)可能需要更长的平衡时间，50ps 或100ps或更长。下面的设置在这个gromacs力场中运行良好。其他力场请参考用户手册(例如在GROMOS96力场中，建议 nstlist = 10 and rvdw = 1.4)。在coulombtype, PME 代表 “Particle Mesh Ewald”静电势。PME是计算长程静电势的最优算法(给出最可信的能量评估，尤其在用 Na+, Cl-, Ca2+等作为补偿离子的体系)。由于这个蛋白具有暴露的带电残基，使系统带有+2静电荷，所以适用 PME算法，更为有益的是用补偿离子使系统处于电中性。 constraints中的all-bonds 选项可以应用线形限制算法确定系统中的所有键长(当 dt > 0.001 ps时尤为重要)。学习一下下面的 mdp 文件。 pr.mdp: 7 define 声明中的–DPOSRE 告诉Gromacs运行位置限制动力学模拟。 constraints 声明如前所述。 all-bonds 设定LINCS 算法限制所有键。 integrator 告诉 gromacs进行何种动态算法(另外的选项“ sd ”代表stochastic dynamics) 8 dt 是每步的时间(我们选择了2fs;但此处的单位一定是 ps ～) nsteps 是运行的步数(总模拟时间 = nsteps * dt)。 nstxout 告诉gromacs轨迹文件收集模拟快照(坐标)的频率( nstxout = 250 且dt = 0.002，所以每 0.5ps收集一张快照) coulombtype 选择gromacs计算原子静电相互作用方法( PME 代表 particle mesh ewald;另外还可以用cut-off )。 Rcoulomb和rvdw 是计算静电和范德华作用的阈值(单位nm，1.0nm=10.0埃) 温度耦合部分非常重要，必须正确填写: Tcoupl = v-rescale [8, 9] (用随机条件重新调解速度的温度耦合类型。) tau_t 温度耦合的时间常数(单位 ps)。必须每个 tc_grps 指定一个，且 顺序对应。 tc_grps 与调温器耦合的组(模型中的每个原子或残基都用一定的索引组 表示) ref_t 代表耦合的参照温度(即动力学模拟的温度，单位 K)。每个 tc_grp 对应一个 ref_t. 当你改变温度时，别忘了改变gen_temp 变量以生成速度。 pcoupl –Parrinello-Rahman 恒压器。 pcoupltype –isotropic指“box”可以平均地向各个方向(x, y,z)膨胀或压缩，来 维持一定的压力。注意 :进行膜模拟时用 semiisotropic。 压力耦合的时间常数(单位 ps )。 tau_p – compressibility –溶剂在每bar的可压缩性(上面的设置是水在300K和1 大气 压下的可压缩性)。 ref_p –压力耦合的参照压力(单位 bar，1大气压~ 0.983 bar )。 grompp –f pr.mdp –c fws_b4pr.pdb –p fws.top –o pr.tpr –maxwarn 3 nohup mdrun –s pr.tpr –o pr.trr –c fws_b4md.pdb –e pr.edr –g pr.log & 用tail 命令检查 pr.log文件。 三 设置非限制性动力学模拟 md.mdp文件和pr.mdp文件相仿。有几处不同，define声明不再需要，因我们不再做位置限制模拟。用于explicit solvation 的md.mdp文件内容 (特殊注释:做真空模拟，去掉温度耦合中的“sol”部分。在有补偿离子的模拟中，为离子加入相应的温度耦合参数) 9 grompp –f md.mdp –c fws_b4md.pdb –p fws.top –o md.tpr –maxwarn 3 nohup mdrun –s md.tpr –o md.trr –c fws_md.pdb –e md.edr –g md.log & 用tail 命令查看md.log 文件。 可以用trjconv命令压缩轨迹文件以节省硬盘空间。 trjconv –f md.trr –o md.xtc 得到*.xtc 文件后就可以删除 *.trr 文件了。 10 四 并行计算、重启、延长模拟计算 1 如何重启一个计算 tpbconv -s prev.tpr -f prev.trr -e prev.edr -o restart.tpr mdrun -s restart.tpr -deffnm myrestart (–deffnm 将mdrun中的所有文件名设成默认名字) 2 如何延长一个计算 tpbconv -f traj.trr -s topol.tpr -e ener.edr -o tpxout.tpr –time $VALUE –until $VALUE 其中$VALUE 单位是ps (例如你要将2ns模拟延长到5ns，则$VALUE取5000) 3 如何设置并行计算 grompp–np N–f md.mdp–c pr.gro–p fws.top–o md.tpr -np标签设定并行计算的节点数，N代表节点数，自行设定。 然后用mdrun_mpi进行并行任务: mpirun –np # /products/gromacs/bin/ mdrun_mpi –deffnm md 11 五 模拟结果分析 1 如何将特定帧的轨迹保存成*.pdb文件 用-dump 选项，如 trjconv -f traj.xtc -s file.tpr -o time_3000ps.pdb -dump 3000 2 用ngmx观察轨迹文件(也可以用VMD观察轨迹文件) ngmx –f md.trr (or md.xtc) –s md.tpr 当观察器启动后，将看到一个多选项的对话框。选择标“protein ”的多选框，点击OK。 12 选择“protein ”可以只看蛋白分子，而不受盒子中另外约3000水分子的影响。 用X-Rotate 上下旋转盒子(鼠标左键向上，右键向下)。用Y-Rotate 左右旋转盒子(左键向左，右键向右)。最下面的 Scale用来放大或缩小视图(左键放大，右键缩小)。 观察模型中的其它组，点击Display > Filter … ，初始对话框就会出现，允许选择观察另外的索引组(如backbone ) 要观察模拟轨迹动画，点击Display > Animate 。动画播放控制在窗口的底部。点击中间的箭头按钮逐帧观看。点向前的双箭头观看整个轨迹动画，点暂停按钮停止动画。点向左的双箭头按钮重置动画。保存并观察*.pdb文件最好的方法是用visual molecular dynamics (VMD) (下载地址:，它是学术免费，且在linux和windows下运行)结果分析Gromacs的一个主要优点是有一系列分析轨迹文件的小程序。 13 3 比较常用的分析工具 3.1 组 make_ndx 程序make_ndx用来生成组(你想分析的某些特定原子或残基的 ID标签)。Gromacs缺省已经定义了一些组，普通分析可能够用了。但如果你想深入分析，则要用make_ndx程序标注模型中的特定项。 如何使用make_ndx建立索引文件(ndx )。为了固定某些特定组，或获得一些能量信息，可以用make_ndx指定这些组。我们来看一个固定蛋白 N端和 C端的例子。通常用 make_ndx建立索引组供grompp程序调用。 在本例中，我们有一个三螺旋的胶原蛋白结构文件，进行位置限制动力学模拟，我们想固定N 端和 C端来进行模拟。首先，确定结构文件(clg_b4md.pdb )的 N端和 C端残基号。用如下简单命令: make_ndx –f clg_b4md.pdb –o clg_ter.ndx 你将看到如下输出信息(我们省略了开头的一些描述性信息)，后面是命令提示符(>) 14 用“ r”命令输入代表三螺旋 N端和 C端的残基号。 注意:你也可以用连接符指定残基范围(如确定残基1到36 ，用> r 1-36) 新建索引组的缺省名字( r_1_36_37_72_73_108 ) 很繁琐，可以用name命令修改。我们在 命令中用索引号# (15 )。 用“v ”命令查看名字是否改成功了。用“ q”保存并退出。 现在怎么固定组呢,简单，在md.mdp文件中加入下面几行: 记住当用新的mdp 文件时，首先用grompp将新索引文件加入 tpr文件。用 grompp的-n 标签，例如: grompp –f md.mdp –c pr.gro –p clg.top –n clg_ter.ndx –o md.tpr 15 3.2 特性研究 g_confrms 要比较最后结构和初始PDB 文件的差异，用g_confrms ( 用g_confrms –h 查看详细信息 )。此程序计算两个结构的最小二乘拟合。 g_confrms –f1 1OMB.pdb –f2 md.gro –o fit.pdb 你将被提示选择一个组(两次都选(组4 ))。程序将RMSD 值，并产生一个输出文件( fit.pdb )。输出文件中包含两个位置重叠的结构。 3.3 g_covar 计算斜方差 也可用于从动态轨迹计算平均结构。如计算1ns动态模拟的后200ps的平均结构: g_covar –f traj.xtc –s topol.tpr –b 801 –e 1000 –av traj_avg.pdb 警告- 平均结构往往较粗糙，需进一步执行能量最小化。 3.4 g_energy 能量数据作图，如压力、体积、密度等 g_energy –f md.edr –o fws_pe.xvg 首先要选择输出(*.xvg)的数据。输出文件是一个电子数据表，可以用Xmgr或Grace打 开。它是一个文本文件，在进行一些小的改动后可以用 Microsoft Excel打开。 如用上面的命令，你将看到如下结果(你的可能不同): 如计算势能，输入“Potential”，回车 再按一次回车 我们得到一个平均势能和RMSD 的总结(单位kJ/mol ) 16 输入如下命令用Grace打开 *.xvg文件: xmgrace -nxy fws_pe.xvg 可以到以下地址下载Grace 只能在linux和unix 上运行。如果没有Grace或Xmgr，可以作为空格分隔文件导入 MS Excel 。 3.5 g_gyrate 测量回旋半径 这个指标用于度量结构的紧密度。此程序计算某(些)原子质量与分子重心的关系。 g_gyrate –f md.trr –s md.tpr –o fws_gyrate.xvg 3.6 g_rms 与 g_rmsdist 计算结构的RMSD 值 用g_rms 计算动态模拟过程中的结构与初始结构的结构偏差。(-dt 10 选项告诉程序每10 帧计算一次) g_rms –s *.tpr –f *.xtc –dt 10 计算与NMR结构的rmsd值用如下命令: g_rms –s em.tpr –f md.trr –o fws_rmsd.xvg 选择4 组( Backbone)计算最小二乘拟合。程序生成一个 rmsd 随时间变化的图( rmsd.xvg )。以空格分隔文件的形式导入 Excel。 选择一个 RMSD 图上平衡的范围(用g_rms计算)。例如: 17 上面的实例是一个1ns模拟(你的结果可能不同)。这个模拟需要延长到完全平衡。 3.7 g_rmsf 计算原子位置的根均方波动( rmsf ) 与g_covar相似，此程序也可以计算平均结构。例如，计算一个 2ns(2000ps )模拟的后500ps的平均结构，用如下命令: g_rmsf –f traj.xtc –s topol.tpr –b 1501 –e 2000 –o traj_rmsf.xvg –ox traj_avg.pdb 上面的例子我们用 200-500ps范围计算平均结构是因为我们看到这段比较稳定(与原始结构比)。命令如下: g_rmsf –s md.tpr –f md.trr –b 200 –e 500 –ox fws_avg.pdb 提示时选择组1“Protein”。 推荐的真空能量最小化的em.mdp文件。先用最陡下降法，再用共轭梯度法。警告～需要用pdb2gmx 重新生成拓扑文件，尤其是当你选择特定组(而非整体系统)计算平均结构时。 18 程序g_rmsf 也可以用来计算温度因子。计算的温度因子可以和 X光晶体结构的温度因子比较。 f md.xtc –o rmsf.xvg –oq fws_bfac.pdb g_rmsf –s md.tpr – 19 仍选择“Backbone” 组。 3.8 do_dssp 计算模型的二级结构 前提是你必须在电脑中(/usr/local/bin)安装了dssp 程序( )。 do_dssp –s md.tpr –f md.trr –o fws_ss.xpm (Protein )。用xpm2ps将xpm 文件转成 eps格式。然后用选择计算组1 ImageMagick 转化程序将eps文件转成 png文件或其他格式文件。 xpm2ps –f fws_ss.xpm –o fws_ss.eps convert fws_ss.eps fws_ss.png 残基数在y 轴，时间(ps )在x轴。看下面的NMR结构: 从上面的 dssp 图上，我们看到3个红色区域代表3 个beta 片层(下图中的黄色部分)。中间的较短区是最不稳定的。下边的图是用pymol 程序( )做的。 20 3.9 g_hbond 计算模拟过程中分子间的氢键的数目、距离或角度 g_hbond –f md.trr –s md.tpr –num fws_hnum.xvg Gromacs 4.0的缺省值为: r?30.35nm a?30º 用-r 和-a 标签设定其它限值。默认地，g_hbond 计算共体受体距离(r DA)。你也可以用-da no 标签改变这种计算，设置 da为no 则计算rHA距离。 3.10 g_saltbr 分析模拟中残基间的盐桥 程序会输出一系列*.xvg文件。给出-/- ，+/-(盐键)和+/+残基间的距离。 g_saltbr –f md.trr –s md.tpr 21