使VMD实时显示gromacs轨迹中原子的受力

文/Sobereva@北京科音2009-Jan-31

显示原子受力有时很重要，但实现起来又是个麻烦事。由于这两天要用这个功能，所以写了一个VMD的tcl脚本，可以用箭头显示指定原子的受力，并且随进度条拖动实时更新，在这里做一下介绍。

VMD当前版本并不读取轨迹中的力，不能在内部直接调用，所以需要手动把轨迹中的力提取出来并读入VMD。在此例中，我们要分析序号为977的原子的受力，序号由0开始算。


这里用到一个技巧，结合gmxdump和grep可以容易地从.trr文件中提取977号原子受力，虽然比较慢，但是很方便：
gmxdump -f all.trr |grep "f. 977" > force.txt
包含f[ 977]的行就被提取出来了。force.txt内容如下：
f[ 977]={ 1.02676e+03, -1.37174e+03, 1.29441e+03}
f[ 977]={ 5.54416e+01, 9.32820e+02, -7.47177e+02}
f[ 977]={ 7.87461e+02, 7.43057e+02, -9.03874e+02}
f[ 977]={ 1.52381e+02, -5.37775e+02, 1.38345e+02}
f[ 977]={-6.90852e+02, -1.33539e+02, -9.30110e+02}
......
应当注意，为了使每一行都对应一帧结构，.mdp的nstxout和nstfout应当相同。（最后一帧未必，因为最后一步的结构不管怎么设nstxout都会被输出，而这一帧未必会正好是输出力的那步，但这个问题无关紧要）

把force.txt放到tcl启动后的默认文件夹（比如d:\study\vmd186，在控制台输入pwd可显示），启动vmd，载入轨迹，在控制台运行：
set i 0;set fx 0;set fy 0;set fz 0
set force [open force.txt r]
while {[gets $force line] >= 0} { //检测何时停止循环读取
scan [string range $line 16 55] "%f,%f,%f" fx fy fz //把force.txt的每行读入并只保留a,b,c的形式，赋值到fx fy fz
set f(x$i) $fx;set f(y$i) $fy;set f(z$i) $fz //赋值到数组
incr i
}
close $force
我们已将全部受力读入一个数组f，例如f(x43)代表第43帧时977号原子受力的x分量。

然后运行下面的内容，定义一个名为showforce的过程，用于绘制箭头
proc showforce {args} {
global vmd_frame
global f
global atom //引用全局变量
graphics 0 delete all //清空已绘制的物体
set sel [atomselect top "index $atom" frame $vmd_frame(0)]
set x [$sel get x]
set y [$sel get y]
set z [$sel get z] //得到$atom原子当前帧的坐标
set init "$x $y $z" //箭头尾巴坐标
set middle "[expr $x+$f(x$vmd_frame(0))/400] [expr $y+$f(y$vmd_frame(0))/400] [expr $z+$f(z$vmd_frame(0))/400]" //箭头锥形和柱形衔接处坐标
set after "[expr $x+$f(x$vmd_frame(0))/300] [expr $y+$f(y$vmd_frame(0))/300] [expr $z+$f(z$vmd_frame(0))/300]" //箭头头部坐标
graphics 0 color red
graphics 0 cylinder $init $middle radius 0.3 filled yes resolution 20 //绘制圆柱
graphics 0 cone $middle $after radius 1.3 resolution 20 //绘制锥形
$sel delete //释放atomselect所占用的内存
}

运行：
trace variable vmd_frame(0) w showforce
说明跟踪vmd_frame(0)变量，这个变量每当ID为0的分子的进度条移动后都会更新，内容是当前帧号。此时便触发showforce过程根据此帧的受力来重绘箭头。

最后定义要显示哪个原子受力，在此例中运行set atom 977
拖动进度条，就可以看到效果了，如图所示的红箭头，箭头的长度与受力大小成正比。


如果不想显示了，只需输入
trace vdelete vmd_frame(0) w showforce
graphics 0 delete all
如果又想显示，仍然是运行trace variable vmd_frame(0) w showforce


为了方便，我将上述内容写进了showforce.tcl脚本。简要来说使用方法如下：
1. 将force.txt和showforce.tcl放入VMD默认文件夹
2. 启动VMD，载入轨迹
3. 在控制台运行source showforce.tcl
4. set atom 977
即可。



showforce.tcl文件内容如下

set i 0;set fx 0;set fy 0;set fz 0
set force [open force.txt r]
while {[gets $force line] >= 0} {
scan [string range $line 16 55] "%f,%f,%f" fx fy fz
set f(x$i) $fx;set f(y$i) $fy;set f(z$i) $fz
incr i
}
close $force


proc showforce {args} {
global vmd_frame
global f
global atom
graphics 0 delete all
set sel [atomselect top "index $atom" frame $vmd_frame(0)]
set x [$sel get x]
set y [$sel get y]
set z [$sel get z]
set init "$x $y $z"
set middle "[expr $x+$f(x$vmd_frame(0))/400] [expr $y+$f(y$vmd_frame(0))/400] [expr $z+$f(z$vmd_frame(0))/400]"
set after "[expr $x+$f(x$vmd_frame(0))/300] [expr $y+$f(y$vmd_frame(0))/300] [expr $z+$f(z$vmd_frame(0))/300]"
graphics 0 color red
graphics 0 cylinder $init $middle radius 0.3 filled yes resolution 20
graphics 0 cone $middle $after radius 1.0 resolution 20
$sel delete
}

trace variable vmd_frame(0) w showforce


PS：还可将此脚本改写，显示每一帧某原子的速度。