使VMD讀入Gromacs產生的trr軌跡中速度信息的方法

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方向上被拉長。如預期地，往右運動的水呈紅色，往左運動的水呈藍色。由于盒子尺寸有限，往右運動的水撞到右邊相鄰盒子的水球中往左運動的水。之后，水球逐漸變成了一個水柱。在水柱剛形成時，若仔細觀察，通過顏色能看到不同區域的水的運動方向不是完全均勻的，與相鄰水球碰撞產生的振蕩造成了水的運動在局部有一定取向性。而在水柱穩定后，由于水的熱運動方向是雜亂無章的，所以紅、綠、藍在各處都均勻混合在一起。