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

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