使用Multiwfn和VMD計算分子表面積和片段表面積
使用Multiwfn和VMD計算分子表面積和片段表面積
文/Sobereva@北京科音 2019-May-27
時常有人問怎么計算分子表面積、怎么計算基團的表面積,其實這用Multiwfn或VMD來做非常簡單,這里就結合例子簡單說一下。
1 關于分子表面的定義
首先要清楚定義表面積的方式非常多,在《談談分子體積的計算》(http://www.shanxitv.org/102)一文中就已經提到了好幾個。其中有一個叫做SASA(溶劑可及表面積),這個面積定義依賴于溶劑探針半徑,通過VMD等程序可以快速計算。
范德華表面是分子表面的一種,但也有不同具體定義。一種廣為接受的以量子化學方式的定義是Bader在JACS, 109, 7968 (1987)中提出的,也就是對于氣相情況,將電子密度為0.001 a.u.(1 a.u.=e/Bohr^3)的等值面作為范德華表面。這種方式定義的范德華表面平滑,并且可以反映電子結構特征(如孤對電子、pi電子等),其范圍比通過范德華球疊加得到的范德華表面略大一些。對于凝聚相環境,考慮到分子間相互作用導致范德華表面穿透,分子體積會減小,故此時建議用電子密度0.002 a.u.的等值面作為范德華表面。計算以電子密度等值面方式定義的分子表面可以用Multiwfn程序非常簡單快速地實現。
下面,筆者以多巴胺分子為例演示計算的方式。此體系結構如下,我們既打算計算其整體表面積,也想計算其氨基部分的表面積。
Multiwfn程序可以在官網http://www.shanxitv.org/multiwfn免費下載,不熟悉者強烈建議參看《Multiwfn入門tips》(http://www.shanxitv.org/167)和《Multiwfn FAQ》(http://www.shanxitv.org/452)。本文使用的是Multiwfn 3.6版。VMD是一個絕佳的化學體系可視化程序,可以在http://www.ks.uiuc.edu/Research/vmd/免費下載,本文使用的是1.9.3。
本文用到的所有文件可以在此下載:http://www.shanxitv.org/attach/487/file.rar。
2 通過Multiwfn計算范德華表面積
2.1 基于量子化學方式產生的密度計算范德華表面積
要基于電子密度計算范德華表面積,自然得先獲得電子密度,一般通過做量子化學計算來實現。使用B3LYP/6-31G*這種比較便宜的計算級別,對于當前目的就已經夠用了。Multiwfn支持大量量子化學程序產生的含有波函數信息的文件來做當前的分析,見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》(http://www.shanxitv.org/379)。這里我們使用常用的Gaussian來產生.wfn文件。
本文文件包里的dopamine.xyz是之前恰當優化后得到的多巴胺結構。將之載入Multiwfn,然后依次輸入
100 //其它功能(Part 1)
2 //轉換文件格式、產生量化程序輸入文件
10 //產生Gaussian輸入文件
dopamine.gjf //輸出的文件名
現在當前目錄下就有了dopamine.gjf,對應B3LYP/6-31G*計算級別。將out=wfn關鍵詞加進去,末尾空一行寫上.wfn文件的輸出路徑,比如C:\dopamine.wfn。然后用Gaussian計算之,得到的.wfn文件已經提供在了本文的文件包里。
Multiwfn的定量分子表面分析功能是專門考察各種函數在分子表面上分布特征的,在http://www.shanxitv.org/159和http://www.shanxitv.org/196文中都有介紹。我們目前的目的不是分析函數在分子表面上的分布特征,而是僅僅要獲得表面積,因此我們下面將把映射到分子表面上的函數切換為一個沒有任何耗時的函數來節約不必要的計算時間。
啟動Multiwfn,載入dopamine.wfn,然后依次輸入
12 //定量分子表面分析
2 //選擇被映射到分子表面的函數
-1 //用戶自定義函數。默認設置下,這個函數值處處為1,所以不會有額外計算耗時
Multiwfn的這個功能默認分析的是電子密度為0.001 a.u.的等值面。假設此處我們想改成0.002等值面來對應凝聚相情況下范德華表面的定義,我們就輸入
1 //選擇用于定義表面的函數
1 //用電子密度等值面定義表面
0.002 //等值面數值(isovalue)
然后選0開始分析。Multiwfn轉眼就分析完了。其它輸出不用管,就看下面這一行輸出
Overall surface area: 648.82739 Bohr^2 ( 181.69020 Angstrom^2)
即整個多巴胺用0.002 a.u.電子密度等值面定義的范德華表面的面積為181.7 Angstrom^2。
2.2 計算片段對應的范德華表面積
接著上一節,這里我們來算一下分子表面上對應氨基區域的局部表面積,這要利用Multiwfn的局部分子表面分析功能。在后處理菜單中依次輸入
12 //輸出某個片段的屬性
3,19,20 //氨基的原子序號
從屏幕上可見下面的輸出
Overall surface area: 99.58145 Bohr^2 ( 27.88565 Angstrom^2)
即曰這個電子密度為0.002 a.u.的等值面上由氨基貢獻的部分為27.9 Angstrom^2,占總面積比例為27.9/181.7=15.3%。
此時程序還問你是否導出locsurf.pdb。如果我們想圖形化觀察分子表面上對應氨基的區域,我們就選擇y來導出它。然后啟動VMD,將這個文件拖到VMD main窗口載入,然后選Graphics - Representation,把Drawing Method設為Points,把Coloring Method設為Beta。然后再把dopamine.xyz也載入進VMD來把分子結構顯示出來,把Drawing Method設為Licorice,背景改為白色,當前看到的圖像如下所示
圖中每個小點對應電子密度等值面上的一個頂點,藍色的區域對應氨基。由此圖可見,Multiwfn給出的對應氨基的局部表面定義得很合理,因此上面給出的氨基的面積是可靠的。
2.3 基于準分子密度計算范德華表面積
如果你沒Gaussian的版權,也不會用任何免費的量子化學程序,也照樣可以用Multiwfn來計算分子和片段的表面積,只要提供一個含有合理的分子結構的文件即可(可以用比如xyz、mol、mol2、pdb等)。此時電子密度就不是直接基于波函數算出來的了,而用的是promolecular density(準分子密度),它是將每個原子在孤立狀態下的電子密度按照分子中原子的坐標簡單疊加來產生的,整個周期表幾乎所有元素在孤立狀態下的高精度電子密度在Multiwfn里直接內置了。這種方式產生的分子的電子密度顯然比較糙,因為沒考慮原子間相互作用導致的電子的轉移、極化。但如果你要求不高,有個定性結果就夠,那還是完全可以接受的。
啟動Multiwfn,載入dopamine.xyz,然后依次輸入
12 //定量分子表面分析
1 //選擇用于定義表面的函數
2 //用一個特殊函數的等值面定義表面
1 //準分子密度
0.002 //等值面數值
0 //開始分析
值得一提的是當前被映射到表面的函數默認被切換為了user-defined function,默認設置下這個函數是個常數,計算它沒有任何額外的耗時。當前面積計算結果為
Overall surface area: 697.18104 Bohr^2 ( 195.23060 Angstrom^2)
這里給出的結果195.2 Angstrom^2和之前基于B3LYP/6-31G*密度給出的181.7 Angstrom^2定性一致。
和2.2節一樣,我們用相同的步驟再計算一下氨基對應的表面積,輸出為
Overall surface area: 112.37312 Bohr^2 ( 31.46768 Angstrom^2)
此處31.5和之前給出的27.9也差不太多,占總面積的比例31.5/195.2=16.1%和之前算的15.3%也相仿佛。這說明,基于直接通過幾何結構構建的準分子密度來考察表面積是完全靠譜的。如果這次也通過locsurf.pdb去看一下分子表面頂點和劃分情況,會發現從肉眼上看不出圖像與上一節的圖有什么差別。
值得一提的是,如果你考察的是很大體系,如果做上述分析發現太耗時(無論是基于量化計算的密度還是準分子密度),可以在主功能12的界面里選擇3 Spacing of grid points for generating molecular surface,然后輸入一個比默認更大的格點間距,比如0.4。間距越大,表面頂點就越稀疏,得到的表面積數值精度就越差,但分析耗時也越低。
3 通過VMD計算SASA面積
啟動VMD,載入dopamine.xyz,然后在命令行窗口輸入下面這一行,就可以計算整體的SASA
measure sasa 1.4 [atomselect top all]
返回的結果為343.08,單位是Angstrom^2。上面的命令中1.4代表探針半徑,探針半徑越小,得到的SASA也會越小。1.4對應于水分子的探針半徑。
將下面的語句復制到VMD命令行窗口里運行,就可以計算氨基部分暴露在SAS(溶劑可及表面)上的面積,即氨基對整體的SASA的貢獻
measure sasa 1.4 [atomselect top all] -restrict [atomselect top "serial 3 19 20"]
返回的結果為59.80 Angstrom^2,占總SASA比例為59.80/343.08=17.4%,比例值和前面我們基于電子密度范德華表面的方式算的比較接近。