- 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有所差异,所以如果使用不同的软件进行氢键分析,会获得不同的结果。