使用Multiwfn圖形化展示分子動力學模擬體系中的孔洞、自由區域
使用Multiwfn圖形化展示分子動力學模擬體系中的孔洞、自由區域
文/Sobereva@北京科音
First release: 2020-Mar-5 Last update: 2021-Sep-14
1 前言
最近有人問我怎么獲得類似下面文獻中的圖
此圖像展示出分子動力學模擬體系中的孔洞區域,或者叫自由區域(free region)也可以,這體現出了盒子里哪些區域沒有被原子占據。
之前筆者寫過《使用Multiwfn可視化分子孔洞并計算孔洞體積》(http://www.shanxitv.org/408)一文,那篇文章的做法是用于展現分子體系的特定的內部孔洞的,沒法像上圖這樣對大量體系展示盒子里所有自由區域,而且也根本算不動。為了能得到上面的圖,筆者在Multiwfn中加入了相應功能,作為主功能300的子功能1出現,在此文介紹一下。此功能還可以給出自由區域的體積的具體數值便于定量分析。
Multiwfn可以在http://www.shanxitv.org/multiwfn免費下載,讀者務必使用2021-Sep-14及以后更新的版本。不了解Multiwfn者可參看《Multiwfn FAQ》(http://www.shanxitv.org/452)和《Multiwfn入門tips》(http://www.shanxitv.org/167)。本文用的VMD是1.9.3版,可在http://www.ks.uiuc.edu/Research/vmd/免費下載。
筆者后來還寫了《使用Multiwfn計算晶體結構中自由區域的體積、圖形化展現自由區域》(http://www.shanxitv.org/617),是本文的重要補充,且介紹了更多作圖技巧,建議讀者閱讀。
2 實例
我們先直接看一個具體計算實例,之后再說相關細節。本例用到的文件是Multiwfn文件包中的examples目錄下的coal.pdb,這是位網友用Lammps程序跑動力學得到的多孔介質體系的一幀,如下所示。明顯能看出有不少稀疏、空白區域。
啟動Multiwfn,然后輸入
examples\coal.pdb
300 //其它功能(Part 3)
1 //觀看自由區域和計算自由區域體積
4 //修改用于平滑邊界的展寬函數
1 //Gaussian展寬
1.8 //用原子范德華半徑的1.8倍作為Gaussian函數的半高全寬(FWHM)
1 //設置格點并開始計算
[按回車] //用默認的原點位置,即(0,0,0)處,這和當前情況一致
[按回車] //用晶胞的三個邊長作為格點數據計算范圍的三個邊長
[按回車] //用0.25埃格點間距。格點間距越小,要算的格點數就越多,耗時就越高,圖像也會越平滑,因此應根據實際情況決定。當前格點間距下的圖像質量已經不錯了
普通8核機子花大約一分鐘就能算完。算完后可以看到
Volume of entire box: 30038.649 Angstrom^3
Free volume: 14783.074 Angstrom^3, corresponding to 49.21 % of whole space
這說明整個盒子總體積,即三個邊長的乘積是30038.6 Angstrom^3,自由體積是14783.1 Angstrom^3,占總體積的49.21%。自由體積是這樣算的:對于每個格子,如果它距離任意一個原子的距離在此原子的范德華半徑以內,這個格子就被認為是占據的。所有沒有占據的格子數乘上格子體積就是自由體積。
之后看到后處理菜單,選擇3,就出現了圖形界面,將窗口右側的Show molecule和Show data range都選上,可看到下圖
上圖中的等值面直觀地展現出了盒子中未被原子占據的區域。為了得到更好的圖像效果,我們選擇4將格點數據導出為當前目錄下的free_smooth.cub文件。啟動VMD后將之載入,進入Graphics - Representation,點Create Rep,把設置改為下面的樣子來顯示出等值面。
最后在控制臺輸入pbc box和color Display Background white顯示出盒子邊框并改成白背景,此時圖像如下
可見效果非常理想。
如果你把Isovalue設小,等值面會占更多的空間區域,如果把Isovalue設大,等值面則會收縮,且較小的等值面會消失。因此可以利用這個設置控制等值面顯示的視覺效果。例如將Isovalue設為0.9后看到的圖像如下,相對于上圖,此時只有最顯著的孔洞還可以看到,而相對次要、較小的孔洞就不顯示出來了
值得順帶一提的是,以上述方式得到的孔洞的等值面與VMD的VDW(范德華球疊加)方式顯示的分子表面是互補的。將分子的Drawing method改為VDW后,VMD顯示的圖像如下所示
可見黃色的對應孔洞的等值面與對應分子的范德華表面緊緊貼合,二者合在一起正好充滿整個盒子。故這兩種顯示方式有高度互補性,應根據實際被研究的體系和問題恰當使用。
3 相關細節說明
下面把Multiwfn上述功能涉及到的各種細節說一下。
使用當前功能可以用各種Multiwfn支持的包含晶胞/盒子信息(即記錄三個晶胞矢量,或者三個邊的長度+夾角)的文件作為輸入文件。比如含有CRYST1字段的pdb文件、GROMACS程序的gro文件、cif文件、VASP的POSCAR、CP2K的restart文件、含有平移矢量(TV)的Gaussian輸入文件,等等,完整列舉見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》(http://www.shanxitv.org/587)的相應部分。此時計算的格點數據的三個邊是分別平行于晶胞的三條邊的。如果你用的輸入文件里沒有晶胞/盒子信息,在當前功能里選擇選項1設置盒子并計算時,程序也會讓你手動輸入盒子原點坐標、邊長、格點間距,但注意此時會假定盒子的三個邊是分別平行于X、Y、Z軸的(因此此時不適合非正交盒子的情況)。
在后處理菜單選擇觀看等值面的時候,只有晶胞/盒子的三條邊是正交的情況才能正確顯示。對于非正交的情況,必須導出cub文件然后用VMD才能正確顯示。
在當前功能的界面里,可以看到Toggle considering periodic boundary condition選項,默認狀態是Yes,即考慮周期邊界條件。如果選擇此選項將之狀態切換為No,則計算耗時可以顯著降低(約一個數量級),但構造格點數據時就不考慮原子的周期鏡像了,這時候會導致誤把盒子邊緣的一些本來是被占據的地方顯示成沒有占據。下圖左邊是考慮周期性的,右邊是沒考慮周期性的,可見考慮周期性是非常重要的。
當前功能里有個選項Toggle making isosuface closed at boundary,默認狀態是Yes。如果切換為No的話,等值面在盒子邊緣就不是封閉的,而是鏤空的,如下所示,部分鏤空地方我用箭頭高亮了。通常來說這樣不如封閉的美觀。這個設置不影響耗時。
此功能效率很高,用于甚至幾萬原子都可以。代碼做了充分的并行化,如果體系很大算起來很慢,有以下幾種方法降低耗時:
(1)用核數很多的機子算。別忘了需要將settings.ini里的nthreads設為實際的CPU的物理核數
(2)用更大的格點間距,從而減少要計算的點數。注意對于很大體系,格點間距太小的話不僅算不動,還會由于格點數太多導致內存存不下,致使程序崩潰
(3)不考慮周期性,前面已經說過了。這樣做的話等值面的可靠性會下降,而且算出來的自由體積通常明顯不準確
當前功能有個選項Toggle calculating smoothed grid data of free regions,默認狀態是Yes,也就是既產生“原始(raw)格點數據”,也產生“平滑后(smoothed)的格點數據”,將之設為No不產生后者的話可以節約幾分之一的時間。本文前面展示的等值面都是基于平滑后的格點數據繪制的。這里說一下兩類格點數據的產生機制。原始格點數據在構造時是循環每個格點,某個點與任意一個原子的距離小于此原子的范德華半徑的話,這個格點的數值就為0,否則為1,因此數值為1的格點就對應于自由區域,前面提到的自由體積也是相當于是所有數值為1的格點數的總和乘上格子體積。這種方式構造的格點數據并不適合展現孔洞區域,比如如果你在后處理菜單選擇1來觀看它的等值面的話,看到的是下圖
可見效果巨爛,就像是奶酪中在有原子的地方挖走半徑等于范德華半徑的窟窿,在當前圖像中甚至連孔洞的基本模樣都看不出來。相比之下,平滑后的格點數據在構造的時候優雅得多。Multiwfn支持通過三種不同的切換函數實現平滑化,即使得原子的邊界不是一刀切,而是隨原子徑向距離增加從1平滑地變化為0。支持的切換函數包括:
(1)Gaussian函數。變化特征依賴于FWHM,數值越大收斂越慢
(2)誤差函數。變化特征依賴于scale值,scale越大收斂越慢
(3)Becke函數。變化特征依賴于迭代次數(n_iter),迭代次數少收斂越慢
你可以在當前功能里通過選項4選擇用哪種函數做平滑,并且之后可以輸入所用參數。這幾個函數結合不同參數時的變化示意圖如下所示,它們等于0.5的位置被設為了碳原子范德華半徑處(3.21 Bohr)
前文的例子我們用的是FWHM為1.8倍原子范德華半徑時的Gaussian函數作為展寬,此時在原子范德華半徑0.9倍處的切換函數恰為0.5。Multiwfn是這樣計算平滑后的格點數據的:對于某個點,起初格點數據為1,令它減去所有原子在此處的切換函數值,如果最后發現為負則把數值設為0。由前面的例子可見,基于這樣平滑后的格點數據顯示出的等值面非常光滑,能很好地勾勒出孔洞的位置,而且還能通過調節等值面數值來決定是只觀看重要孔洞還是也觀看次要孔洞。如果大家對所得圖像不滿意,可以嘗試不同展寬函數結合不同參數。
使用Gaussian函數用于平滑化對大多數情況效果不錯,有一個缺點是對原子分布得很密集的體系比如C60晶胞,你會發現平滑后的格點數據處處為0。這是由于Gaussian函數收斂得比較慢,每個原子都會影響到較遠的區域,即對較遠處的平滑化的格點數據都有明顯負貢獻。故當原子很密集時,可能導致很大范圍區域的格點算出來的平滑化的格點數據都為負,最后都被設為了0。解決辦法就是改用誤差函數或Becke函數來做展寬,Multiwfn會令原子范德華半徑處正好等于這倆函數為0.5的位置,而且由于如上圖所示它倆在0.5附近衰減得較快,因此可以避免上述用Gaussian函數時的問題。如果堅持用Gaussian函數來平滑化,可將FWHM設得比平時更小(比如1.2)來解決,這使得它收斂得更快。
當前功能有個選項Set scale factor of vdW radii for calculating free volume and raw free region,默認的scale factor是1.0。這個設的是構造前述的原始格點數據時給原子乘的范德華半徑,因此影響基于原始格點數據看的等值面的效果,而且由于自由體積是對原始格點數據統計來得到的,因此這個半徑也直接決定了計算出的自由體積。比如如果將前例的刻度因子設為1.5的話,算出來的自由體積就只有16.6%了。這個選項對平滑后的格點數據無影響,因為在它產生的過程中不考慮這個刻度因子。