5、你可以试着用其他热力学数值对Time Step作图,熟悉使用NAMD plot。然后关闭NAMD plot。
6、首先,打开VMD,选择File→New Molecule菜单项,在弹出的文件浏览中找到 common目录下的文件usq_ws.psf,点击按钮Load。这样我们首先载入了泛素蛋白的结构信息。(注意图形窗口中是空白的,并没有显示出蛋白结构,因为psf文件不存储各个原子的坐标)
7、此时Molecule File Browser窗口上端一栏应该显示Load files for: 0: ubq_ws.psf (图)说明如果再次载入新文件,所载入的文件将与usq_wb.psf结合起来一同显示。我们将载入原子运动轨迹文件,和psf结构文件配合即可显示出蛋白质原子的运动轨迹。因此点击Browse,找到1-2-sphere文件夹下的ubq_wb_eq.dcd,点击Load,这时关闭Molecule File Browser窗口。可以在图形窗口中看到已经载入的泛素和球状水体。
8、我们将要使用一个脚本(script)计算RMSD值。用写字板打开1-2-sphere目录下的脚本rmsd.tcl,可以看到脚本的内容为:
set outfile [open rmsd.dat w]; set nf [molinfo top get numframes]
set frame0 [atomselect top \set sel [atomselect top \# rmsd calculation loop
for {set i 1 } {$i < $nf } { incr i } { $sel frame $i
$sel move [measure fit $sel $frame0] puts $outfile \}
close $outfile 脚本的每一句含义解释如下:
? ?
set outfile [open rmsd.dat w]: 新建文件rmsd.dat作为输出文件
set nf [molinfo top get numframes]:声明变量nf,并把轨迹文件中的帧数(frame)
赋值给变量nf ?
set frame0 [atomselect top \声明变量frame0,frame0储存第一帧中蛋白质骨架非氢原子的位置。 \and backbone and noh\表示我们选择的是蛋白质骨架上的所有原子,并且不包含氢。定义这个变量是为了以第一帧作为基准,其他各帧都和第一帧比较计算RMSD值(注意我们在定义的时候,
RMSD值是以原子运动的平均位置为基准的。这里为了计算方便我们以第一帧为基准) ?
set sel [atomselect top \:声明变量sel,变量暂时没有确定的值。以后每个循环依次将各帧的蛋白质骨架中的非氢原子放入这一变量中。
? ?
# rmsd calculation loop:注释:RMSD值循环计算开始
for {set i 1 } {$i < $nf } { incr i } {:和C语言中的for语句用法一样,这一语句表示:定义变量i=1,每次循环i的增量为1,i > nf 时跳出循环。相当于C语言中的语句:for(int i=1;i < nf;i++) 此后由{ }括起来的内容就是循环体。每循环执行一次
? ? ?
$sel frame $i : 将变量sel设定为所要比较的那一帧(第i帧)。
$sel move [measure fit $sel $frame0] :将第i帧和第1帧的蛋白骨架联配(align)。 puts $outfile \$sel $frame0] :计算第i帧和第1帧之间比较的RMSD值,并输出到我们在第一条语句中新建的文件:rmsd.dat。
9、如果刚才关闭了VMD,请再次打开,重新载入common/ubq_ws.psf 和 1-2-sphere/ubq_ws_eq.dcd。下面,在VMD主窗口中选择Extensions →TK Console打开控制台,在控制台中以cd命令转换目录到namd-tutorial/1-2-sphere, 然后输入source rmsd.tcl 运行脚本。脚本会在1-2-sphere目录下生成文件rmsd.dat。这一文件记录了蛋白骨架的RMSD值随时间的变化。
得到rmsd.dat数据文件后,我们可以用多种绘图软件进一步作图。这里以最常见的Excel为例,做出蛋白质骨架RMSD随时间的变化。
10、关闭VMD,打开Excel,选择 “文件”→“打开” ,在弹出的文件浏览中找到namd-tutorial/1-2-sphere目录,文件类型选择所有文件,找到rmsd.dat打开它(如果选择不显示已知文件类型的扩展名,可能看不到后缀“.dat”)。
11、这时会弹出“文本导入向导”对话框(图)。直接单击“完成”导入文本到Excel。
导入后,数据位于A1至A9,就是蛋白质骨架的RMSD值。我们还需要输入时间以进行作图。前面在配置文件中,我们的动力学模拟时间设定是2500 timesteps,1 timestep = 2fs(飞秒),所以共计5000fs = 5ps(皮秒),又因为我们设定dcd文件的输出频率是250 timesteps输出一次(参见2.3.1配置文件),所以每次输出对应0.5ps。因此在B1-B9中依次输入1, 1.5,…,4.5,5。注意:0.5ps是第一帧,被用作计算的基准帧。
12、为了作图方便,我们将A列的RMSD值数据剪切→粘贴到C列。A列空出即可(图)
13、选择“插入”→“图表”开始作图。在弹出的对话框中选择“XY 散点图”中的右下角子图表“折线散点图”。点击“完成”即可看到做出的初步图形(图)
0.90.80.70.60.50.40.30.20.100123456系列1 14、进一步进行修饰并加上图题等内容可以得到最终结果(图)。
RMSD0.90.80.70.60.50.40.30.20.1000.511.5球形水体中泛素RMSD(5ps)T/ps22.533.544.555.5 其实上面的这张RMSD图中的信息量很少,因为本次分子动力学模拟时,能量平衡进行的时间太短,写入dcd文件的频率又低,2500 timesteps 中,一共输出了10帧,因此作图呈离散点状。下面是一张典型的RMSD值图:
可以看到,RMSD值在最后逐渐趋于恒定,在一定值上下波动。这说明体系已经达到能量平衡。虽然从我们自己做的RMSD图中得不出明显的结论,但是基本分析方法我们已经掌握了。