GROMACS教程【02】分子动力学模拟结果分析常用命令(二)

  • A+
所属分类:杂谈天下

1.导出固定的帧

很多时候我们需要导出分子动力学轨迹中某一时刻的结构构象,比如比较动力学前后的构象差异,就需要导出轨迹中的第一帧和最后一帧。在Gromacs中,导出固定帧十分容易,常用命令如下:

gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o first.pdb -dump 0

以上命令可输入轨迹中的第一帧,即0 ps时的结构,通过“-dump”可指定时间,单位为ps。2.间隔固定时间导出帧

除了单帧导出外,可能还需要间隔帧导出,比如制作分子动力学轨迹动画时,直接使用原轨迹制作,轨迹文件较大,容易造成卡顿,就可以间隔一定的时间导出帧,比如每隔100 ps导出一帧到PDB文件,可输入以下命令:

gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o traj.pdb -dt 100

以上命令导出的全部帧都保存为一个PDB文件,如果需要将每一帧都输出为独立的PDB文件,可以通过“-sep”命令来指定。3.计算平均结构

gmx covar”主要用来计算并对角化(质量加权的)结构协方差矩阵,也可以用来计算平均结构

gmx covar -s md_0_10.tpr -f md_0_10_center.xtc -b 500 -e 800 -av traj_aveg.pdb

以上命令可以得到分子动力学模拟轨迹中500到800 ps的分子平均结构,结果输出为traj_aveg.pdb

此外,上期内容为大家介绍了使用“gmx rmsf”进行RMSF分析,计算分子结构的均方根涨落,也可使用该命令获取某一时间段的平均结构,如输入以下命令获得分子动力学模拟轨迹中500到800 ps的分子平均结构:

gmx rmsf -s md_0_10.tpr -f md_0_10_center.xtc -b 500 -e 800 -ox traj_aveg.pdb

4.叠合两结构并输出RMSD

gmx confrms”可将两个分子结构进行叠合,并且计算RMSD值,值得一提的是,叠合的两个结构的原子数可以不相同,只要用于叠合的两个索引组一样即可。常用命令为:

gmx confrms -f1 structure_1.gro -f2 structure_2.gro -o fit.pdb -n1 index_1.ndx -n2 index_2.ndx

通过以上命令,叠合的两个结构会写入至“fit.pdb”文件中。5.分子动力学轨迹聚类聚类(Clustering) 是按照某个特定标准(如距离)把一个数据集分割成不同的类或簇,使得同一个簇内的数据对象的相似性尽可能大,同时不在同一个簇中的数据对象的差异性也尽可能地大。

在分子动力学中,可以依据RMSD对轨迹进行聚类,把相似的构象聚集在一起,聚类的常用命令为:

gmx cluster -f md_0_10_center.xtc -s md_0_10.tpr -cutoff 0.8 -b 1000 -wcl 10 -method gromos -sz

其中:-b:指定从多少ps开始分析,一般情况下轨迹前期结构变化较大,可以从稳定后如1000 ps开始分析;-method:指定算法,默认为linkage,可以选jarvis-patrick、monte-carlo、diagonalization、gromos;-cutoff:截断值,可以根据实际情况需要反复尝试,cutoff设置的越大,簇里的帧之间的RMSD允许差异就越大,聚类的数目越少。通过以上命令,每个簇当中最有代表性的结构一起被输出为多帧文件cluster.pdb,可以通过各类可视化软件如VMD、PyMOL等进行观看。6.氢键分析氢键在分子动力学中相互作用的分析中十分常用,在Gromacs中可通过“gmx hond”命令对氢键的数量、稳定性等进行分析。在氢键分析中需要指定两个不同的组,它们必须完全相同或者彼此之间无任何重叠。

分析氢键常用的命令为:

gmx hbond -s md_0_10.tpr -f md_0_10_center.xtc -n index.ndx -num hbond_num.xvg -dt 10 -life hbond_life.xvg -ac hbond_ac.xvg

通过以上命令可输出hbond_ac.xvg、hbond_life.xvg、hbond_num.xvg文件。其中:hbond_ac.xvg输出了氢键的平均存在周期(forward lifetime),可以作为衡量氢键稳定性的一个指标;hbond_life.xvg主要需要关注p(t);p(t)是t的概率密度函数,就是统计不同的t出现的频率得到的结果,即氢键从0能维持到t时刻的概率;此外,t*p(t)对时间积分,用统计学语言说这个值就是时间的期望值,也就是平均寿命;hbond_num.xvg:第二列表示的是氢键数目,而第三列显示的是Pairs within 0.35nm,表示距离在0.35nm不考虑氢键键角判据的氢键数,实际意义不大。需要注意的是,hbond程序使用几何准则来对氢键加以判定,当氢键供受体距离小于3.5埃且氢-供体-受体所成角度小于30度时,即认为其为一个氢键。而不同软件如VMD、PyMOL在默认情况下对于氢键的判定与Gromacs有所差异,所以如果使用不同的软件进行氢键分析,会获得不同的结果。

weinxin
我的微信公共号
我的微信公招扫一扫

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: