使用Multiwfn計算分子的長寬高以及顯示分子的主軸

使用Multiwfn計算分子的長寬高以及顯示分子的主軸

文/Sobereva @北京科音

First release: 2018-Jun-10  Last update: 2021-Oct-23

0 前言

經常有人問怎么計算分子的長寬高，之前筆者在《談談分子半徑的計算和分子形狀的描述》（http://www.shanxitv.org/190）一文中也曾討論分子的長寬高怎么算，但那個博文用的蒽的例子過于簡單，而實際分子形狀往往比較復雜，而且結構文件里的分子的朝向可能歪斜著，沒法直接用那篇博文里的做法。本文介紹通過Multiwfn程序非常方便、準確地計算長寬高的做法。讀者請務必使用最新版本Multiwfn，可以在官網http://www.shanxitv.org/multiwfn下載。如果讀者對Multiwfn不了解，看《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn FAQ》（http://www.shanxitv.org/452）。沒看過《談談分子半徑的計算和分子形狀的描述》的讀者應先看一下了解相關基本知識。

1 計算分子長寬高的原理

計算長寬高的第一步，是先把分子“擺正”，讓分子最長的軸沖著某個笛卡爾軸。手動在可視化程序里旋轉不準確，而且比較麻煩。一個自動化的做法是計算分子的慣性矩張量，如下所示：
 

式中x,y,z是原子坐標，m是原子質量。這個3*3矩陣的三個本征矢對應分子的三個主軸，彼此正交，相應的本征值是繞著主軸的轉動慣量。最小的本征值對應的主軸是分子的最長軸，最大的本征值對應的是分子的最短軸。因此，如果我們旋轉分子，讓分子的三個主軸朝向正好對應三個笛卡爾軸，分子就被擺正了。代碼實現其實很容易，就是利用慣性矩張量的本征矢矩陣來對各個原子的坐標做個線性變換而已。

擺正之后，尋找x坐標最小的原子，減去其范德華半徑，記為xmin；尋找x坐標最大的原子，加上其范德華半徑，記為xmax，則x方向尺寸就是xmax-xmin。類似地，對y、z方向也都這么計算。這樣分子的長、寬、高就都有了。

當然，能用長寬高來描述一個分子，前提是這個分子必須能夠通過矩形來比較好地描述，諸如一個分子是月亮形，或者比如碳納米管側面接出來一個長鏈，或者是一個四面體形的團簇，此時長寬高這種描述就沒什么意義，且有誤導性。而且，注意很多大分子有很多不同構象，不同構象的長寬高算出來可能差異很大，這種情況也不適合用長寬高這種概念說事。另外，還有的分子用球形、柱形、橢球形等形式描述更恰當，此時也不要局限于用長寬高來描述。

2 計算長寬高實例

例如當前有一個pdb文件，用可視化程序打開后，看到其朝向和位置如下圖所示。我們想計算其長寬高

.mol、.xyz、.pdb、.fch等含有結構信息的Multiwfn可以識別的格式都可以直接載入到Multiwfn里面計算長寬高。什么格式含有什么信息看《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。啟動Multiwfn后，載入此pdb文件，然后依次輸入
100  //主功能100
21   //輸出各種描述體系結構的信息
size  //計算分子尺寸
此時程序自動把體系平移和旋轉，各個主軸都分別沖著各個笛卡爾軸了。屏幕上顯示以下信息
Farthest distance:   21(H )  ---   45(H ):    10.250 Angstrom
vdW radius of   21(H ): 1.200 Angstrom
vdW radius of   45(H ): 1.200 Angstrom
Diameter of the system:    12.650 Angstrom
Radius of the system:     6.325 Angstrom
Length of the three sides:    12.650     7.513     7.619 Angstrom

從中可以看到，體系中相距最遠的兩個原子是H21和H45，二者的Bondi范德華半徑也輸出了，將倆原子的范德華半徑加到它們之間的距離上，就是屏幕上輸出的體系的直徑，屏幕上輸出的半徑是其一半。可見利用Multiwfn計算分子半徑/直徑很方便。Length of the three sides后面的三個值就是體系的三個邊長，可以被視為長寬高。

如果想把長寬高圖形化地展現一下，選擇屏幕上的選項1 Visualize the new orientation and molecular box，會看到下圖，圖中藍色方框的三個邊長正對應于分子的長寬高。下圖把Ratio of atomic size設為了4.0，此時圖上原子球的半徑正對應于Bondi范德華半徑，可見藍色方框正好卡著范德華表面。

如果想把坐標已經被Multiwfn旋轉后的當前體系導出成pdb文件，從而能在VMD等程序中顯示，可以選擇2 Export the geometry in new orientation as new.pdb in current folder，此時當前目錄下就會產生new.pdb。此文件中有個CRYST1字段，記錄了此體系的長寬高。如果用VMD程序（http://www.ks.uiuc.edu/Research/vmd/）載入此pdb文件，就會讀入盒子信息，并可以繪制出來。

把new.pdb拖到VMD里載入，然后在命令行窗口輸入pbc box，盒子就被顯示出來了。將背景改為白色，選Display-Orthographic用正交視角，在Graphics-Representation里把體系設置以VDW風格顯示，然后把長寬高手動ps到圖上去，就看到了下圖的效果。

順帶一提，利用VMD的繪圖命令，還可以把主軸繪制出來。把new.pdb載入到VMD后，把以下內容復制到VMD的命令行窗口運行（最后要按一次回車確保最后一條命令運行了）
set lena [molinfo top get a]
set lenb [molinfo top get b]
set lenc [molinfo top get c]
set cenx [expr $lena/2]
set ceny [expr $lenb/2]
set cenz [expr $lenc/2]
draw color red
draw cylinder "$cenx $ceny 0" "$cenx $ceny $lenc" radius 0.15 resolution 20
draw color green
draw cylinder "0 $ceny $cenz" "$lena $ceny $cenz" radius 0.15 resolution 20
draw color blue
draw cylinder "$cenx 0 $cenz" "$cenx $lenb $cenz" radius 0.15 resolution 20

調節顯示方式為CPK，然后就會看到下面的圖像，紅綠藍三個圓柱對應于三個主軸

下面是兩篇已發表的文章中使用Multiwfn的這個功能的實例

下圖中對分子不同片段表面分別著色突出顯示的做法見《使用Multiwfn和VMD計算分子表面積和片段表面積》（http://www.shanxitv.org/487）。