发布时间 : 星期日 文章Gromacs教程II-MD结果分析 - 图文更新完毕开始阅读469eaee7dd3383c4ba4cd243
Gromacs 中文教程（II）：结果分析
淮海一粟
MD结果分析
模拟结束后，就可以分析数据了。分析包括三个阶段。首先，有必要对模拟的质量进行检查，如果检查结果表明模拟良好，那么就可利用该模拟对所研究的问题作出回答；最后，不同的模拟结果可以合并。
注：文件名应该反映文件内容，这根据你的模拟系统不同而不同。这里我们假定使用默认文件名，那么就会产生以下文件：
? ? ? ? ? ?
topol.tpr: 模拟开始时包含完整系统描述的输入文件； confout.gro: 结构稳健，包含最后一步的坐标和速度； traj.trr: 全部精确轨迹，包括位置、速度和随时间变化的力
traj.xtc: 压缩的轨迹文件，只包含低精度(0.001 nm)的坐标信息； ener.edr: 随时间变化的有关能量数据 md.log: 模拟过程的日志
附注：许多分析工具都能生成.xvg文件。这些文件能用xmgr或xmgrace查看，也可用python脚本程序xvg2ascii.py在终端显示出来。
Each group writes one report. For general questions a single answer should be given in the report. Questions specific to each simulation should be given in table, indicated with ( T ). Answers to questions from one section should be combined in a single table if possible.
先检查一下结果
在开始分析前，要确定模拟是否正确地完成。有很多原因会导致模拟中断，尤其是与力场和系统平衡不充分引起的问题。为了检查模拟是否正确完成，运行gmxcheck程序：
gmxcheck -f traj.xtc
看看是否执行了10纳秒的模拟。
==Q== How many frames are in the trajectory file and what is the time resolution? ( T )
另一个重要的信息源是日志文件。在文件md.log的末尾有模拟过程的统计数据；包括内存和CPU的使用情况和模拟时间。看看日志文件的末尾，使用'less'命令时，你可以用'G' (shift-g) 命令，跳到文件末尾。
==Q== How long did the simulation run in real time (hours), what was the simulation speed (ns/day) and how many years would the simulation take to reach a second? ( T )
==Q== Which contribution to the potential energy accounts for most of the calculations?
Don't be afraid to use the Gromacs online manual, to search in the Gromacs mailing list or even to use Google to get information about the terminology used by Gromacs.
结果可视化
好玩的环节开始了。虽然很多分析都能归结为从轨迹文件中提取图像，MD当然首先要关注系统的移动。来看看轨迹文件。
首先用gromacs提供的查看器ngmx来看看。虽然该软件的完善程度和视觉效果不及其他查看器，但它能够在拓扑文件信息的基础上画出键。其他查看器可能隐含远程键，这可能导致这些键被认为太长而不画出，或者会在非常接近的原子之间画出键。这是对模拟结果分析的一个常见错误源。使用ngmx载入拓扑和轨迹文件：
ngmx -s topol.tpr -f traj.xtc
看看程序菜单，试试不同的选项。播放动画。观看过程通过右边的选项控制。右击或左击选择选项来改变查看。
==Q== What happens if the protein diffuses over the boundary of the box?
为了视觉美观，我们将从轨迹文件中提取1000帧 (-dt 10)并忽略水分子（当软件询问时，选择 Protein）。而且，我们还将忽略边界的跨越，作出连续轨迹（-pbc nojump）。为了做这些工作，我们使用瑞士军刀般的gromacs工具trjconv，该工具有1001个选项组合。我们用它写出一个多模型的pdb文件，从而能在Pymol中观看。
trjconv -s topol.tpr -f traj.xtc -o protein.pdb -pbc nojump -dt 10
在Pymol中提取轨迹文件： pymol protein.pdb
当所有的帧装入完毕，播放动画。 Mplay
动画播放时，其他控制键仍在运行，可以用鼠标旋转、放大或缩小图像，也可以改变分子外观。 Spectrum Show cell
如果没错的话，你现在能看到蛋白质扩散、翻转跳跃。但我们对内部运动更感兴趣，而不是总体行为。在Pymol中，你能使用命令intra_fit将其他所有帧与第一帧对齐。随后，你可以用定向工具设定蛋白质中心： intra_fit protein and (name c,n,ca) Orient
现在，所有的帧应该都对齐了，你可以看到蛋白质的哪一部分移动得更厉害。这些差异将在以后定量分析。
当然，在cartoon模式下，蛋白质看起来更舒服，试试这个命令： show cartoon
因为.pdb文件里面没有二级结构信息，你可能会看到碳骨架是个粗管状，而看不到正确的二级结构。Pymol可以自动计算蛋白质的二级结构，但只计算一帧，并将其映射到其它的帧。例如，下述命令可以计算第一帧的二级结构： dss
通过设定状态，用于计算的帧可以改变： dss state=1000
最后，同时查看所有帧并检查蛋白质的柔性和刚性部分。 set all_states=1
请随便练习Pymol的使用。试试放大柔性或刚性区域，并检查侧链构象。使用'ray'和 'png'制作一份图像，即使浪费点CPU时间也不要紧。但记住，如果图像太复杂，可能会导致pymol的插件ray-tracer崩溃，这种情况下，你可以直接用'png' 得到屏幕上的图像。
The following part, up to 'quality assurance', is optional and it may be best to first finish the other sections.
如果你有足够机时做此教程，或者你已经做完了前面的功课，那就来做个不错的电影吧。你可能注意到，这些轨迹噪音非常多，那是热噪音，是正常的；但是对于制作好的电影会有影响。我们可以屏蔽这些高频运动，只保留低速和平滑的总体移动。为了达到这个目的，使用g_filter程序：
g_filter -s topol.tpr -f traj.xtc -ol filtered.pdb -fit -nf 5
现在，在Pymol中倒入轨迹文件。计算二级结构 (dss)并显示(show cartoon)，隐藏碳骨架上的侧链 (hide lines, not (name c,n,o)) 并给蛋白质上色，然后orient分子设好视角。现在开始制作电影了： viewport 640,480 set ray_trace_frames,1 mpng frame_.png
现在退出Pymol (quit) 并显示文件路径 (ls)。你能看到，文件的数量多了好些，包括1000个图像。以每秒30帧的速度，将会产生30秒的电影。下载mpeg_encode程序和参数文件movie.param，用它来产生单帧图像的电影（你可能需要编辑参数文件来改变文件名）： mpeg_encode movie.param
质量确认
进行了最初的轨迹图像查看后，该对模拟的质量进行彻底检查了。这个质量检查(QA)包括热动力参数（如温度、压力、势能和动能）的收敛情况。更一般地说， QA试图评价（模拟）是否达到平衡状态。 结构上的收敛也需要检查，这个用起始结构和平均结构的均方差 (RMSD)来表示。随后，还要检查邻近的周期性图像之间没有互相作用，因为这将导致非物理学效应。最后，QA检查包括原子的均方差，这个可以与晶体学数据b-factors进行比较。
能量收敛
我们首先从能量文件中提取一些热动力学数据。研究以下性质：温度、压力、势能、动能、单位盒子体积、密度和盒子尺寸。这些大部分性质已在系统准备步骤中检查过了。用工具g_energy进行能量分析。该程序读出能量文件，也就是模拟过程中产生的扩展名为.edr的文件。g_energy程序将会问需要提取什么参数并将产生一幅图像。输入下面的命令：