使VMD根據pdb文件中的CONECT字段設定原子連接關系

雖然pdb文件當中定義了CONECT字段用于明確定義原子間連接關系，然而VMD并不會利用之，無論有沒有CONECT字段，VMD在讀取pdb后都是根據距離范圍自動猜測原子間連接關系。這種處理方式往往不能滿足我們的特殊要求，利用DynamicBonds的顯示方式也并不能完全解決問題。雖然也有conect2psf程序（http://www.ks.uiuc.edu/Development/MDTools/conect2psf/）可以將CONECT的連接關系轉換為psf文件，載入psf文件后就可以向VMD提供自定義的拓撲信息，然而conect2psf程序只有IRIX 5.x下預編譯好的版本，一般用戶沒法用。本文提供一個Tcl腳本，可以讓VMD載入pdb文件時能夠根據CONECT字段的信息設定連接關系。

首先簡要介紹一下pdb文件中的CONECT字段。CONECT字段一般出現在原子坐標定義之后，諸如這樣：
CONECT 1787 1786 1788 1789                                                      
CONECT 1788 1787 1793  
它代表了1787號原子和1786、1788、1789號原子相連，1788原子和1787、1793號原子相連。普通的殘基、水分子不需要CONECT字段描述原子連接關系，因為根據原子類型和鍵長可以很明確地確定哪些原子是相連的。而對于涉及到非標準鍵長、特殊殘基（比如血紅素）、特殊殘基與標準殘基相連的情況，程序在判斷哪些原子間應當相連時存在不確定性，CONECT字段提供的連接信息就用來解決這個問題。另外，二硫鍵也是靠CONECT來描述的。

在VMD中，每個原子都有連接關系列表，用戶可以查看和修改。在Tcl窗口中輸入下面的命令選擇3至5號原子
set sel [atomselect top "serial 3 to 5"]
然后輸入$sel getbonds，返回比如：
{20 8} {10 21} {5 19 18}
這代表當前選擇范圍中的1號原子（對應 serial 3）與index 20、index 8原子相連；當前選擇范圍中的2號原子（對應 serial 4）與index 10、index 21原子相連；當前選擇范圍中的3號原子（對應 serial 5）與index 5、index 19、index 18原子相連。

我們可以修改連接關系，比如讓選擇范圍中的1號原子只與index 4原子相連，讓2號不與任何原子相連，讓3號在原有基礎上還與index 8原子相連，那么就輸入：$sel setbonds {4 {} {5 19 18 8}}
此時會看到圖形窗口中連接關系立刻更新了。你可能會覺得圖上顯示得很奇怪，有很多化學鍵只顯示了一半。這是因為兩個原子的連接關系列表里都有對方時才肯定會顯示完整的鍵，如果A的連接列表里有B，而B里的沒有A，那么往往只會顯示出從A延伸向B的半個鍵。

下面這個Tcl腳本就是根據上述原理通過讀取CONECT字段來設定VMD中原子間連接關系。拷貝到Tcl窗口下即可運行。
natm是體系中總原子數，nskipline代表在CONECT段落之前有多少行，這些行會被空過去。itype如果為0，CONECT字段定義的連接關系就是最終原子間連接關系；如果為1，代表CONECT字段定義的是附加的連接關系，而不是將原有連接關系給覆蓋掉。
此腳本允許CONECT字段中每個原子最多與12個其它原子相連（實際上這也是VMD自身的上限）。CONECT字段的第一列，即參考原子的序號，要求從小到大排序。不用給所有原子都用CONECT設定連接關系，需要設哪些原子連接關系就寫哪些原子即可。

#Programmed by Sobereva, 2012-Jan-25
set natm 8312
set nskipline 8313
#itype=0 means reset connectivity by user-defined list, =1 means add self-defined connectivity list to the original one
set itype 0
set rdpdbcon [open "D://CM//my program//Multiwfn//vtx.pdb" r]

#Skip other lines
for {set i 1} {$i<=$nskipline} {incr i} {
gets $rdpdbcon line
}

#Cycle each atom
set ird 1
for {set iatm 1} {$iatm<=$natm} {incr iatm} {

if {$ird==1} {
 for {set i 1} {$i<=12} {incr i} {set cn($i) 0}
 gets $rdpdbcon line
 scan [string range $line 6 84] "%d %d %d %d %d %d %d %d %d %d %d %d %d" self cn(1) cn(2) cn(3) cn(4) cn(5) cn(6) cn(7) cn(8) cn(9) cn(10) cn(11) cn(12)

 set tmplist {}
 #Formation of connectivity list
 for {set i 1} {$i<=12} {incr i} {
  if {$cn($i)==0} {break}
  lappend tmplist [expr $cn($i)-1]
 }
}

if {$self==$iatm} {
 #puts Atom\ serial:\ $iatm\ \ User-connectivity:\ $tmplist
 set sel [atomselect top "serial $iatm"]
 if {$itype==0} {
  $sel setbonds "{$tmplist}"
 } else {
  set orglist [$sel getbonds]
  $sel setbonds "{[concat [lindex $orglist 0] $tmplist]}"
 }
 $sel delete
 set ird 1
} else {
 set ird 0
}

}
close $rdpdbcon


我們來看一個使用例子。我們想把多巴胺分子的氮原子和對位的羥基氧和羥基氫原子連上。這三個原子分別對應于pdb文件中的3、2、22號原子（對應于VMD里的index 2、index 1、index 22）。那么首先在多巴胺分子的pdb的末尾從
...
HETATM   21  H   MOL A   1      -1.829   2.693  -0.089  1.00  0.00           H
HETATM   22  H   MOL A   1      -3.743  -1.368   0.461  1.00  0.00           H
END
改為
...
HETATM   21  H   MOL A   1      -1.829   2.693  -0.089  1.00  0.00           H
HETATM   22  H   MOL A   1      -3.743  -1.368   0.461  1.00  0.00           H   <--注：此為第23行
CONECT 2 3   <--注：不要寫在CONECT 3 2 22的后頭，如前所述，必須先定義原子編號小的再定義大的。
CONECT 3 2 22
CONECT 22 3
END
之后將此pdb文件載入進VMD。根據pdb文件里的信息將腳本中的natm設為22，nskipline設為23，將腳本中的文件路徑設對，確認itype為1，然后拷貝進tcl窗口中，圖像就會立刻變成這樣，和預期的一致：



實際上，之所以寡人要寫這個腳本，是與寡人近期研究工作有關。寡人正在為Multiwfn程序添加一個新功能，也就是利用Marching tetrahedron算法做分子表面分析，這個算法可以將等值面細分為一堆三角形。為了檢驗程序代碼是否正確，寡人將這個算法產生的大量表面頂點作為原子，通過pdb文件載入VMD以顯示之，下圖的例子是一個水分子的電子密度為0.001的表面，共9654個頂點：



以line方式繪制出的VMD自動判斷的連接關系完全是一團糟：



寡人寫的程序中全部頂點的連接關系會輸出到pdb文件的CONECT字段，利用上面的腳本重設原子連接關系后（即itype應為0），得到下面的圖，正確顯示出分子表面已經被化成一堆三角來描述了。很多分子可視化軟件中在顯示格點文件等值面圖時都可以用mesh方式來顯示，會看到網狀圖案，實際上原理與此處是一致的。