使VMD實時顯示gromacs軌跡中原子的受力
顯示原子受力有時很重要，但實現起來又是個麻煩事。由于這兩天要用這個功能，所以寫了一個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：還可將此腳本改寫，顯示每一幀某原子的速度。