計算不同z位置水能形成氫鍵數的VMD Tcl腳本

有人請我寫一個計算垂直于溶液界面的即不同z值處的水能形成氫鍵數的腳本，下面是我的VMD Tcl腳本，或許有人也用得著。

氫鍵判據用的是常用的35度3.5埃的幾何判據，當然也可以直接在腳本里改。計算方法是，只要有水分子有一個原子在某一層里，則這個水分子就認為屬于這一層的水。對于每一幀，計算屬于每一層的水selin與其它物質selbig之間的氫鍵數，氫鍵包括了這一層中的水作為氫鍵受體和供體兩種情況，其數目分別為代碼中的變量a和b。并且加上這一層水內部之間的氫鍵數（變量c）的2倍。a+b+2c除以這一層的水數，作為這一幀這一層的每個水的平均氫鍵數。腳本中循環軌跡中的每一幀，最終得到這一層平均氫鍵數。nonum變量記錄有多少幀在所設定的范圍里沒有水，這些幀不計算。#后面那行用于調試目的，要考察每幀結果就去掉開頭的#。

首先運行下面的腳本，來加載實現這個功能的子程序
proc numhbavg {sel fps1 fps2} {
set selin [atomselect top $sel]
set selbig [atomselect top "same resid as exwithin 3.5 of $sel"]
set k 0.0
set nonum 0
for {set i $fps1} {$i<=$fps2} {incr i} {
$selin frame $i
$selin update
$selbig frame $i
$selbig update
if {[$selin num]!=0} {
set a [llength [lindex [measure hbonds 3.5 35 $selbig $selin] 0]]
set b [llength [lindex [measure hbonds 3.5 35 $selin $selbig] 0]]
set c [llength [lindex [measure hbonds 3.5 35 $selin] 0]]
set k [expr $k+($a+$b+2*$c)*3.0/[$selin num]]
#puts "fps:$i $a+$b+[expr 2*$c] num_water:[expr [$selin num]/3.0] avg:[expr $k+($a+$b+2*$c)*3.0/[$selin num]]"
} else {incr nonum}
}
if {[expr $fps2-$fps1+1]==$nonum} {return "no result"}
return [expr $k/[expr $fps2-$fps1+1-$nonum]]
}

然后下面的循環會調用這個子程序來輸出每一層的平均氫鍵數，這里假設要計算z=4.0~5.6埃的數據，間隔為0.1埃，且限定20<x<40，20<y<40，每計算完一層輸出一次。計算結果是幀數范圍為100~150內的平均值。注意x/y/z的上下限范圍不要恰頂著體系的邊界，最好至少留出5埃的距離。選擇范圍中resname SOL代表了水。

for {set i 40} {$i<=55} {incr i} {
set k [expr $i*0.1]
set now [numhbavg "same resid as resname SOL and x<40 and x>20 and y<40 and y>20 and z<[expr $k+0.1] and z>=$k" 100 150]
puts [format "%4.2f %4.2f %5.3f" $k [expr $k+0.1] $now]
}

我這里隨便算一個主要由水構成的普通的體系，水形成的氫鍵數大概在3.1左右，標準放寬到4.0埃，40度，則可形成氫鍵數約為3.6。在冰中由于結構十分有序，可形成4個氫鍵，在液態情況下分子的動能必然造成氫鍵的破壞，所以結果是很合理的。
幀數范圍越大、xy平面越大計算越慢，這個腳本計算速度比較慢，不要一下將范圍設得太大。

輸出結果如下，前兩列代表統計的z值范圍，第三列是水的平均氫鍵數

4.00 4.10 3.015
4.10 4.20 3.161
4.20 4.30 3.201
4.30 4.40 3.159
4.40 4.50 3.237
4.50 4.60 3.130
4.60 4.70 3.201
4.70 4.80 3.338
4.80 4.90 3.201
4.90 5.00 3.041
5.00 5.10 3.122
5.10 5.20 3.182
5.20 5.30 3.160
5.30 5.40 3.309
5.40 5.50 3.189
5.50 5.60 3.327