使VMD读入Gromacs产生的trr轨迹中速度信息的方法

使VMD读入Gromacs产生的trr轨迹中速度信息的方法

Making VMD load the velocity information in the trr trajectory generated by Gromacs

文/Sobereva@北京科音    First release: 2012-Jan-6


1. 前言

Gromacs的trr轨迹文件中不仅包含坐标,还可以包含速度、受力等信息。只要让mdp文件中nstvout等于nstxout,trr里每帧结构信息都将伴随对应的速度信息。gro文件中也记录着速度信息,也就是第45列后面的三列。

VMD里面每帧都有vx,vy,vz三个字段用于记录速度,如同每帧的x,y,z字段记录原子坐标。但是,一直到目前最新的1.9.1 beta版本中,VMD仍然在读取trr或gro文件时只读取坐标不读取速度,因此vx,vy,vz三个字段在读取trr/gro后仍然是空白。竟然丢弃速度这么重要的信息,VMD的这一点十分令人匪夷所思,何况解决这个问题对开发者来说只是举手之劳。而VMD在读取lammps的轨迹时则已经可以读取速度、受力等诸多信息了。


2. 思路

让VMD能读取trr文件中的速度,我想出来三个办法:
1 修改VMD的molfile插件。各种格式的文件在VMD中都是通过这个插件来读取,修改其中读取trr文件的代码,然后重新编译VMD就可以达到目的。但是C语言看起来太累,而且重新编译毕竟是麻烦事,尤其是对于图形程序,总得跟乱七八糟的图形库文件相斗争。尽管这个方法是最好的,但经过初步研究后我不打算这么做,还是留待以后由官方修改molfile插件吧。

2 修改Gromacs的trjconv代码,使转换出的trr或xtc文件中的坐标部分记录速度信息,然后在VMD载入这文件,再将其“坐标”通过tcl脚本挪到真实体系轨迹的速度字段里。但是这种做法显得古怪,改trjconv的代码也同样是麻烦事。

3 第三个办法就这是本文将介绍的方法。首先将trr文件的每一帧都转换成相应的gro文件,然后在VMD里照常读取轨迹文件,通过tcl脚本将这些gro文件里记录的速度信息写入相应帧的vx,vy,vz字段里。这种方法的好处是不需要修改和编译VMD或Gromacs的代码,原理很透明,还可以通过修改脚本实现一些特殊目的。不足之处是使用tcl脚本批量读取gro文件中的信息相对于VMD内部利用molfile插件载入轨迹文件慢不少,但仍是可以接受的。另一个缺点是转换出的gro文件总大小是相应trr文件的6倍,不过一般绝不至于没有地方存它们(何况本文的脚本很灵活,可以将整条轨迹的速度信息分次逐步载入完,对硬盘临时空间其实没有绝对要求)。


3. 方法

假设我们不仅想在VMD中观察MIO.trr中的坐标,还想将其中原子的速度信息也载入,以供后续分析(比如通过脚本分析局部的温度)。那么首先应当确保生成MIO.trr所用的mdp文件中的nstvout等于nstxout。此轨迹生成后,建立一个文件夹,比如vel。然后执行比如
trjconv -f MIO.trr -s MIO.tpr -o miku/t.gro -sep
选择system。假设轨迹有101帧,就会在miku目录下生成t0.gro, t1.gro, t2.gro ... t100.gro。

在VMD里载入一个与MIO.trr相对应的gro文件(比如生成相应tpr文件所用的gro文件),此时VMD中会有1帧,通过选择delete frame将这帧清掉,然后将MIO.trr载入到里面。此时VMD中帧的编号是从0到100,这与.gro文件名里的数字相对应。

启动TkConsole,将以下脚本复制进去执行。这里假设我们把生成的.gro文件都放到了/sob/water/miku/里面,并且文件名都是以t开头后面跟着帧号。若非如此,请根据实际情况修改rootname变量。fpsinit和fpsend对应于帧号的起始和终止。第四行中的all代表读取所有原子的速度信息。假设我们只想读取体系中水的速度,那么在执行trjconv的时候应当选择Water或SOL那项,并且将脚本第四行all替换为waters。适当修改脚本前四行并且适当利用trjconv,就可以读取任意时间范围内体系中任意组成部分的速度。

set rootname /sob/water/miku/t
set fpsinit 0
set fpsend 100
set sel [atomselect top all]
set natom [$sel num]

for {set fps $fpsinit} {$fps<=$fpsend} {incr fps} {
set realname $rootname$fps.gro
puts Loading\ $realname
set rdgro [open $realname r]
gets $rdgro line
gets $rdgro line
$sel frame $fps
set tmpvx {};set tmpvy {};set tmpvz {}

for {set iatm 0} {$iatm<=[expr $natom-1]} {incr iatm} {
gets $rdgro line
scan [string range $line 45 69] "%f %f %f" grovx grovy grovz
lappend tmpvx $grovx
lappend tmpvy $grovy
lappend tmpvz $grovz
}

$sel set vx $tmpvx
$sel set vy $tmpvy
$sel set vz $tmpvz
close $rdgro
}

这个脚本的原理是:循环帧号,每次循环中打开对应帧的.gro文件。空过去.gro中前两行,然后一行行地读取后三列(x/y/z方向的速度)。每读一行时都把刚读到的速度添加到tmpvx、tmpvy、tmpvz这个临时速度列表的末尾。读完一帧gro后将临时速度列表挪到所选原子范围的vx、vy、vz字段中。

在读取时还可以顺便做一些运算。速度有正负,若感兴趣的只是速度大小,那么对于x方向,就可以把lappend tmpvx $grovx改为lappend tmpvx [expr abs($grovx)]


4 一个使用速度信息的例子

在较新版本的VMD的绘图方式中,Coloring method里面有Trajectory-Velocity,选择它以后每个原子都会根据速度的大小用不同颜色显示,考察一堆原子的运动行为的时候往往挺有用。如果我们不将速度信息载入vx,vy,vz字段的话,那么这个选项等于毫无意义。Velocity这种着色方式只能根据速度的总大小进行着色,却不能根据速度的分量进行着色。然而在Trajectory-User里面可以看到VMD能根据User,User2,User3和User4字段的数值进行着色,因此为了达到根据速度分量进行着色的目的,我们可以修改上面的脚本,将x,y,z方向速度分别载入进User,User2,User3里,即将相应部分改为
$sel set user $tmpvx
$sel set user2 $tmpvy
$sel set user3 $tmpvz

这里,我们利用速度分量着色的方法,直观地考察下真空中的由11417个水构成的纳米水球在x方向上施加了2V/nm的均匀电场下,在100ps内原子在x方向上的运动速度的动态变化。在模拟中用了周期边界条件(无控温空压),因此水球被极化后会与相邻盒子的水球撞在一起。

按照前述方法载入轨迹,并将x,y,z方向速度信息分别载入User,User2,User3字段。然后将Graphics-Color...-Color Scale-Method设成BGR(我比较喜欢这种色彩刻度)。建立一个显示方式,所选范围设为name OW,Coloring method设为Trajectory-User-User,Drawing Method设为Points,Size设为3,Material设为Transparent。Trajectory标签页中Color Scale Data Range输入-0.3和0.3然后点Set,这样设可以让速度的差异通过色彩充分地展现。在TkConsole输入pbc box -on使盒子显示出来。这样设好后,蓝点和红点将分别代表朝着x负方向(向左)和正方向(向右)有较大速度的氧原子,绿点代表x方向速度接近0的氧原子。播放轨迹,如下面的视频所示:



可以看到,在电场作用下,水球迅速在x方向上被拉长。如预期地,往右运动的水呈红色,往左运动的水呈蓝色。由于盒子尺寸有限,往右运动的水撞到右边相邻盒子的水球中往左运动的水。之后,水球逐渐变成了一个水柱。在水柱刚形成时,若仔细观察,通过颜色能看到不同区域的水的运动方向不是完全均匀的,与相邻水球碰撞产生的振荡造成了水的运动在局部有一定取向性。而在水柱稳定后,由于水的热运动方向是杂乱无章的,所以红、绿、蓝在各处都均匀混合在一起。