使用Multiwfn定量化和圖形化考察分子的平面性(planarity)
使用Multiwfn定量化和圖形化考察分子的平面性(planarity)
文/Sobereva@北京科音
First release: 2021-Sep-16 Last update: 2021-Oct-15
0 前言
分子的平面性(planarity)是分子的一個很重要的結構特征,有很多問題都和平面性密切相關。對pi共軛體系,由于其平面性明顯影響pi電子共軛能力,因此和芳香性、導電性、光學性質等都有密切關系。而且平面性也影響分子間pi-pi相互作用、形成堆積結構的能力。
之前有一次網上有人問我怎么衡量分子的平面性,我仔細搜索發現以前并沒有什么很理想、很普適的方法專門考察這個問題。雖然有很多文章都是通過特定的二面角來體現平面性的,但是這種做法局限性很大,如果是比如聯苯這種體系還好,靠兩個苯環之間的二面角就能體現整體的平面程度,而碰上平面性同時由很多二面角共同決定的情況,比如碗烯,就很難靠二面角來表征平面性了。為了讓分子平面性這個重要特征有可靠、簡單又普適的方法去考察,筆者提出了兩個定量表征平面性的參數MPP和SDP,以及一種圖形化展現平面性的方法,近期已發表在了Journal of Molecular Modeling, 27, 263 (2021) https://doi.org/10.1007/s00894-021-04884-0。這種分析已經在免費、流行的Multiwfn程序中實現。
下面,將對筆者提出的考察平面性的方法的思想先做簡要介紹,更多細節請看原文。之后給出幾個在Multiwfn中考察平面性的具體例子,使讀者既了解這種方法的價值也了解怎么在Multiwfn中方便地實現。讀者請務必使用2021-Sep-16及之后更新的Multiwfn,可在http://www.shanxitv.org/multiwfn免費下載。不了解此程序者參見《Multiwfn入門tips》(http://www.shanxitv.org/167)和《Multiwfn FAQ》(http://www.shanxitv.org/452)。使用Multiwfn按照本文介紹的方法研究平面性時請同時引用Multiwfn程序啟動時顯示的程序原文和上述的Journal of Molecular Modeling文章。
1 定量化和圖形化考察平面性的方法
筆者提出的考察平面性的方法依賴于擬合平面。首先要對整個分子或被考察的局部區域根據其中原子的坐標利用最小二乘法擬合出一個平面,并且要求這個平面穿越這些原子的幾何中心。這樣得到的平面是對被考慮的原子所處的平面的最佳的描述,可以用平面方程Ax+By+Cz+D=0表示。
之后可以計算分子平面性參數(molecular planarity parameter, MPP),它是各個原子距離擬合平面的方均根偏差,如下所示
其中N_atom是被考慮的原子總數。d_i代表i原子距離擬合平面的距離,可以由下式得到
距離d并不區分原子是在擬合平面的哪一側。為了能夠予以區分,筆者定義了距離擬合平面的含符號距離d_s,也即把上式的絕對值符號給去掉。d_s數值為正和為負對應原子在擬合平面的不同兩側
偏離平面跨度(span of deviation from plane, SDP)參數定義如下,其中d_s_max和d_s_min分別是所有被考慮的原子中d_s最正和最負的值
MPP和SDP這兩個參數都可以定量衡量平面性,都是數值越小平面性越強、結構越接近平面。完全平面的體系這兩個指標同時精確為0。但是它們的特征不同,展現出的信息有明顯互補性。MPP體現的被考慮的區域的總體平面性,而SDP衡量的是被考慮的區域中偏離擬合平面的最大程度,即體現了垂直于擬合平面的最大跨距。二者的共性和差異在后文的例子中將會看到。
d_s這個量除了用于計算SDP外還另有用處。筆者發現如果對原子根據其d_s值以不同顏色進行著色,可以非常鮮明直觀地體現出各個原子相對于擬合平面的位置和程度,這對于展現體系結構特征非常有益,可以管這叫做“d_s著色表示法”。
Multiwfn計算上面這些量的功能是主功能100的子功能21的一部分,也可以在主菜單里直接輸入MPP快速進入。進入后輸入要考察的原子序號范圍,立刻就能得到MPP和SDP值。之后還可以產生pqr文件,里面的原子charge屬性對應的就是原子的d_s值,可以放到VMD程序里對原子根據charge屬性進行著色來直觀展示各個原子的d_s值。這個功能不僅可以處理單個結構,還可以提供軌跡文件,讓Multiwfn一次性對里面的一批結構進行計算。
使用這個功能可以使用Multiwfn支持的任意包含體系結構信息的文件作為輸入文件,支持的格式非常豐富,比如pdb、xyz、mol、mol2、gro、cif、Gaussian和ORCA的輸入輸出文件、POSCAR、fch、mwfn等等,詳見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》(http://www.shanxitv.org/379)。輸入的軌跡文件目前只支持xyz格式,格式介紹見《談談記錄化學體系結構的xyz文件》(http://www.shanxitv.org/477)。
下面就演示用Multiwfn來對不同體系計算上面提到的這些量。
2 實例1:考察輪烯(annulene)分子的平面性
輪烯是一類環狀烴類體系,當碳數n為偶數時通式是CnHn。根據n的不同,由于位阻效應和芳香性的差異,輪烯有的是純平的有的是扭曲的。Multiwfn程序包的examples目錄下的[14]annulene.xyz是ωB97XD/def2-TZVP級別優化的[14]annulene輪烯分子結構,如下所示,可見整體有一定平面性但又不是完全平面。此例我們就用上文介紹的方法定量考察它的平面性
啟動Multiwfn,然后輸入
examples\[14]annulene.xyz
MPP //進入考察平面性的功能
h //代表考慮所有非氫原子,因為氫通常不是我們感興趣的。這里也可以輸入具體原子序號
本文介紹的方法耗時特別低,對于很大體系都不花什么時間,對于當前這種小體系更是一瞬間就得到了結果,屏幕上顯示的信息如下所示
Plane equation: A= 1.00000 B= -0.00000 C= 0.00233 D= 0.00000
Deviation of atom 1(C ) to the plane: -0.21768 Angstrom
Deviation of atom 2(C ) to the plane: 0.21768 Angstrom
Deviation of atom 3(C ) to the plane: 0.21430 Angstrom
Deviation of atom 4(C ) to the plane: -0.21430 Angstrom
Deviation of atom 5(C ) to the plane: -0.07768 Angstrom
Deviation of atom 6(C ) to the plane: 0.07768 Angstrom
Deviation of atom 7(C ) to the plane: -0.04457 Angstrom
Deviation of atom 8(C ) to the plane: 0.04457 Angstrom
Deviation of atom 9(C ) to the plane: -0.02456 Angstrom
Deviation of atom 10(C ) to the plane: 0.02456 Angstrom
Deviation of atom 11(C ) to the plane: 0.09971 Angstrom
Deviation of atom 12(C ) to the plane: -0.09971 Angstrom
Deviation of atom 13(C ) to the plane: -0.03585 Angstrom
Deviation of atom 14(C ) to the plane: 0.03585 Angstrom
Maximal positive deviation to the fitted plane is 2(C ): 0.21768 Angstrom
Maximal negative deviation to the fitted plane is 1(C ): -0.21768 Angstrom
Molecular planarity parameter (MPP) is 0.127145 Angstrom
Span of deviation from plane (SDP) is 0.435354 Angstrom
由上可見,一開始先顯示了擬合的平面方程的參數,然后顯示了各個原子偏離擬合平面的含符號距離,即d_s值。然后輸出了d_s最正和最負的兩個原子,前者減去后者就是后面給出的SDP值。上面信息最后顯示出[14]annulene上的碳原子對應的MPP參數值是0.127埃,SDP是0.435埃。這樣的值是大是小?一會兒我們和其它尺寸的輪烯做個對比,結合圖像考察就能充分明白了。
當前程序問你是否導出pqr文件,這里輸入y,然后再輸入保存的路徑,就得到了pqr文件。現在我們在VMD程序中利用此文件里的d_s值對原子進行著色。VMD可以在http://www.ks.uiuc.edu/Research/vmd/免費下載,筆者用的是VMD 1.9.3。
啟動VMD后,將pqr文件拖入VMD main窗口里載入,進入Graphics - Representation界面,將Drawing Method改為CPK,然后把Coloring Method改為Charge。在Trajectory標簽頁里把Color Scale Data Range下面文本框里的色彩刻度下限和上限分別改為-0.5和0.5(輸入后按回車生效),當前單位是埃。著色的上下限完全由自己決定,圖像效果好即可,但不同體系間橫向對比時必須統一。默認色彩變化是紅-白-藍,而筆者習慣用藍-白-紅,故進入Graphics - Colors,將Color Scale標簽頁里的Method改為BWR。最后把背景改為白色(可以按照《VMD初始化文件(vmd.rc)我的推薦設置》http://www.shanxitv.org/545里的做法定義改背景顏色的快捷鍵),此時就看到了下圖左側。為了便于讀者更好地看清楚相對位置關系,用licorice風格顯示的側視圖在下圖右側也給出了(色彩刻度條是自行ps上去的,可以把Graphics - Colors界面里顯示的色彩刻度條直接挪上去再手動加上上下限標簽)
由于當前計算沒有考慮氫原子,所以氫原子都是白色的,不用管。那些基本處在擬合平面上的碳原子的顏色都是白色的,說明d_s值非常小。明顯處于擬合平面上方和下方的碳原子分別以紅色和藍色清晰地展現了出來,顏色越深說明偏離擬合平面越遠。
以相同做法,在前述的J. Mol. Model.文章中筆者對n=4到n=22的所有n為偶數的[n]annulene都計算了MPP、SDP以及d_s著色圖,匯總如下
由上圖可看出利用本文介紹的方法,可以將不同尺寸的輪烯的平面性特征非常清晰地展現,既定量化成了具體數值,通過原子著色圖也能很容易地了解MPP和SDP數值大小的成因。上圖體現出[4]annulene和[6]annulene,也相當于環丁二烯和苯,是完全平面的。[18]annulene和[22]annulene也幾乎是平面的,MPP很小而SDP也不大,原子基本都是白色。[8]annulene是所有體系中MPP最大的,確實整體來看它彎得最厲害,最不平面。[10]annulene的MPP不如[8]annulene大,說明整體平面性還比[8]annulene好點,但是它的SDP卻是所有體系里最大的,這是因為圖中特別藍和特別紅的那兩個原子偏離整體的擬合平面特別顯著,因此偏離平面的跨度特別大,這兩個原子的d_s值也分別對應計算SDP參數時的d_s_min和d_s_max。通過對比可看出MPP和SDP在衡量平面性時是互補的,有各自的價值,沒有必然的正相關性。
以后大家研究自己的體系時,可以將實際算出來的MPP、SDP以及原子著色的情況與上面展示的不同輪烯的情況做對比,來評估自己的體系的平面性大致算是什么程度,是像比如[18]annulene那樣近乎平面,還是類似[14]annulene那樣有一定平面性,還是像[8]annulene那樣平面性極差。
3 實例2:考察BDHN-TTF晶體結構中分子的平面性
使用Multiwfn也可以對整體結構中的某個部分考察其平面性。作為例子,這里我們考察BDHN-TTF分子晶體中的一個分子的平面性。晶體的cif文件可以在http://www.shanxitv.org/attach/618/BDHN-TTF.cif下載。
將這個cif載入Multiwfn,在主功能0里直接看到的晶胞的結構如下所示(要在上方菜單中選Other settings - Toggle showing cell frame把晶胞邊框也顯示出來)
可見分子被晶胞截斷了,沒法直接考察完整的分子。為了能讓一個完整的分子在體系中央,可以用Multiwfn向晶胞的第2、3個方向各復制延展一次構成超胞,并且讓Multiwfn把被截斷的分子保留完整,做法在《Multiwfn中非常實用的幾何操作和坐標變換功能介紹》(http://www.shanxitv.org/610)里有介紹。現在關閉主功能0的圖形窗口,然后在Multiwfn主菜單里輸入
300 //其它功能(Part 3)
7 //對當前體系做幾何操作
19 //構建超胞
1 //晶胞的第1個方向保留原樣
2 //晶胞的第2個方向平移復制成原先的兩倍
2 //晶胞的第3個方向平移復制成原先的兩倍
20 //令分子保持完整
0 //觀看當前的體系
當前體系結構如下所示
可見晶胞已經夠大了,而且如上面紅圈所標注的,中間有一個完整的分子,我們就對它來考察平面性,這需要先得到這個分子里的非氫原子的原子序號。做法是先選界面上方的Other settings - Toggle showing hydrogens讓氫不顯示,然后選界面右側的Show label復選框顯示出原子序號,適當放大體系和旋轉視角,記錄中間分子中任意一個原子的序號,比如113。然后點Tools - Select fragment,輸入113,整個分子就被高亮了,而且其中的所有非氫原子的序號都顯示在文本框里了,如下所示
把文本框里的所有序號復制出來備用:24,64,75,77-78,81,84,107,109-111,113,115,233,235-237,239,241,274,285,287-288,291,294,318。之后點OK,再點RETURN按鈕關閉圖形窗口。接下來就可以效仿前例考察平面性了。依次輸入
-10 //返回到主功能300
0 //返回到Multiwfn主菜單
MPP //進入考察平面性的功能
24,64,75,77-78,81,84,107,109-111,113,115,233,235-237,239,241,274,285,287-288,291,294,318 //把中央分子的原子序號粘貼進來(如果是用的Windows,在Multiwfn窗口標題欄點右鍵選“編輯”-“粘貼”)
然后在屏幕上看到以下輸出
Molecular planarity parameter (MPP) is 0.239126 Angstrom
Span of deviation from plane (SDP) is 0.796772 Angstrom
由于MPP不算小,所以此分子的整體的平面性并不好,而且SDP頗大,說明肯定有某些原子極度偏離整體分子對應的擬合平面。
然后選y導出pqr文件,并把它載入到VMD里。在VMD里設成白背景,文本窗口輸入pbc box命令顯示出盒子邊框。然后進入Graphics - Representation,把Drawing method改為Licorice,Bond Radius改為0.1。然后點擊Create Rep按鈕新建一個表示,在Selected Atoms文本框里輸入same fragment as serial 113(113是中間分子中的一個原子的序號,當前語句選擇這個原子所在的整個片段,見《VMD里原子選擇語句的語法和例子》http://www.shanxitv.org/504了解更多信息)。然后把Drawing Method改為CPK,Bond Radius改為0.5,Coloring Method改為Charge,Material改為Edgy Shiny(并確保VMD main窗口里選了Display - Rendermode - GLSL,否則材質的效果無法如實展現)。然后像上一節例子一樣把色彩刻度范圍設為-0.5~0.5,用藍-白-紅色彩過渡方式著色。之后可以在Graphics - Materials對EdgyShiny材質的定義做一些微調以讓圖像效果更好。此時看到的圖像如下所示
可見大多數非氫原子的顏色都不淡,明顯說明整體的平面性較差。顏色最深的原子是體系中sp3雜化的碳,它們偏離擬合平面距離最遠,也正是因為它們的存在使得分子結構發生了嚴重扭曲。
4 實例3:計算18碳環在分子動力學軌跡中的平面性變化
筆者對18碳環(cyclo[18]carbon)做過大量研究并發表了系列文章,匯總見http://www.shanxitv.org/carbon_ring.html。其中《揭示各種新奇的碳環體系的振動特征》(http://www.shanxitv.org/578)一文介紹的Chem. Asian J., 16, 56 (2021) DOI: 10.1002/asia.202001228論文中筆者做了18碳環的從頭算動力學模擬,并且指出18碳環在平面內和平面外都有高度的柔性。這一節通過前文介紹考察平面性的方法研究一下整個動力學模擬過程中的18碳環的平面性的變化。298.15 K下動力學模擬的前500 fs軌跡文件是examples目錄下的C18_MD_500.xyz,每1 fs保存一次結構,因此共有501幀。初始幀(第0幀)是幾何優化后的18碳環極小點結構,是精確平面的結構。
啟動Multiwfn然后輸入
examples\C18_MD_500.xyz
MPP
a //代表選擇所有原子,a意為all
a //代表考慮軌跡文件里的所有幀。也可以輸入特定的幀號范圍
由于本文介紹的考察平面性的算法簡單,而且Multiwfn計算效率高,所以即便有501幀也在一瞬間就算完了。Multiwfn在當前目錄下輸出了MPP_SDP.txt,其中一行對應一幀,第1、2、3列分別是幀號、MPP和SDP值。可以將之通過Origin之類程序并繪制成曲線圖,如下所示(下圖是對完整的2000 fs的模擬軌跡繪制的,軌跡文件可以在http://www.shanxitv.org/multiwfn/extrafiles/C18-MD.xyz下載)
可見上圖非常清晰地展現出模擬過程中18碳環平面性的變化,其平面性是周期性波動的,大約220 fs一個周期。由于18碳環的平面外變形比較容易,因此可以看到MPP和SDP的最大值都不小。1015 fs那一刻的d_s著色圖也附在上圖了,可見這種圖把這個時刻的各個原子位置關系清晰地展現了出來,原子越紅說明處于整體平面的越上方,越藍處于整體平面的越下方。
Multiwfn還在當前目錄下輸出了ds.pqr,包含了我們考慮的總共501幀的結構,每一幀的原子的Charge屬性對應的是相應幀結構下原子的d_s值。利用這個文件,我們可以對在VMD播放軌跡時實時根據原子的d_s值進行著色,生動地展現各個時刻原子的平面性特征。不過,VMD 1.9.3版不支持從pqr文件里加載多幀,需要通過筆者寫的VMD腳本實現實時著色的目的,做法是:啟動VMD并將C18_MD_500.xyz載入,然后把Multiwfn目錄下的examples\scripts\ds.tcl腳本文件以及ds.pqr文件都拷到VMD目錄下,最后在VMD的命令行窗口輸入source ds.tcl執行腳本。之后當你播放軌跡時,原子就會根據ds.pqr里記錄的相應幀的d_s值進行著色了。注意這個腳本默認的色彩刻度范圍是-0.4~0.4,可擇情修改。為了避免有初學者對這段操作犯難,我還特意錄了操作視頻:http://www.shanxitv.org/multiwfn/res/ds_color.mp4,其中還使用了VMD的RMSD trajectory tool工具消除了碳環在模擬過程中的平動和轉動以便觀看。當前在VMD窗口里播放軌跡時看到的動畫如下,可見頗為鮮活生動地展現了各時刻原子的位置。
5 實例4:考察分子中局部區域的平面性
此例演示一下怎么考察分子中局部區域的平面性。這個體系稱為富勒烯捕手,它有兩個碗烯片段可以通過pi-pi堆積相互作用牢牢抓住富勒烯。結構文件可以在http://www.shanxitv.org/attach/618/catcher.pdb下載。現在我們對其中一個碗烯片段里的碳來考察平面性。顯然,這需要獲得這些碳的原子序號才能在Multiwfn里輸入。你可以肉眼一個一個去看記錄序號,而最方便的方法莫過于用GaussView打開,然后按住r鍵拖動鼠標通過劃框方式把要被考察的原子全都選中成為黃色(可以選多次,被選中的會累加),然后在Tools - Atom selection界面里把選中的原子序號從文本框里復制出來,如下所示,可知被選中的原子的序號是9-10,13-30
啟動Multiwfn,然后輸入
catcher.pdb
MPP //進入考察平面性的功能
9-10,13-30
馬上看到結果
Molecular planarity parameter (MPP) is 0.352135 Angstrom
Span of deviation from plane (SDP) is 0.882635 Angstrom
之后導出pqr文件,按照第2節的做法在VMD里根據原子的d_s參數進行著色展示。下圖是把色彩刻度下限和上限分別設為-0.8和0.8的情況,繪制風格用的是Licorice。此圖中粉色和藍色分別突出顯示了碗烯整體平面上方和下方的原子。
6 總結
本文介紹了筆者在J. Mol. Model., 27, 263 (2021)提出的定量考察整個分子或體系中某個局部區域的指標MPP和SDP,其定義簡單、嚴格而且普適,既可以用于比較同一個體系不同構型/構象的平面性差異,也可以橫向對比不同體系的平面性。在討論芳香性、pi共軛、分子堆積等問題都可以結合這兩個參數討論。MPP和SDP也可以視為新的分子描述符,可用于機器學習預測各種性質的目的。本文還介紹了直觀展現不同原子相對于擬合平面位置的“d_s著色表示法”,它可以令研究者根據一張圖上的原子顏色立刻了解到各個原子在體系平面的哪一側、偏離平面的程度有多大,完全不需要擔心由于視角問題造成的視覺上的誤判。這些方法都已經實現在了免費的Multiwfn程序中,并且如上文的例子可見,操作非常簡單,而且計算特別快速,對單個分子和晶體都可以考察,不僅可以考察單一結構還可以直接對整條軌跡進行分析,而且結合VMD程序在軌跡播放時還能對原子動態著色來生動地展現其偏離體系平面的情況。本文介紹的平面性的考察方法既方便又實用,預計在以后會被越來越多的化學研究者使用,也非常歡迎大家向同行推廣。