使用Multiwfn繪制構象權重平均的光譜

1 前言

眾所周知，柔性分子有大量構象，并且其中許多低能的構象在當前研究的溫度下都是有不可忽略的分布比例的，熱平衡狀態時它們的比例通常使用Boltzmann公式基于高精度計算的各個構象的相對自由能來計算，見《根據Boltzmann分布計算分子不同構象所占比例》（http://www.shanxitv.org/165）。

不同構象的光譜往往是不同的，甚至有定性差異，特別是圓二色譜的構象敏感性頗高。為了得到能夠和實際較相符的理論光譜，我們計算柔性體系光譜時應當對分布>10%（更嚴格時，分布>5%）的所有構象分別算出的光譜曲線，并按照構象分布比例進行加權線性組合。下面管這種方式得到的理論光譜叫做權重平均光譜或構象平均光譜。

Multiwfn（官網http://www.shanxitv.org/multiwfn）雖然主要用處是波函數分析，但經常關注此程序的人一定知道Multiwfn也具有非常強大的繪制光譜的功能，在《使用Multiwfn繪制紅外、拉曼、UV-Vis、ECD、VCD和ROA光譜圖》（http://www.shanxitv.org/224）中做了介紹，沒看過的話一定要先仔細看看。從Multiwfn 3.4版開始，Multiwfn支持直接繪制構象權重平均光譜，這大大簡化了繪制這種光譜的操作步驟，本文就通過實例介紹怎么實現，比起使用任何其它程序都要方便得多。

本文的量化計算使用Gaussian09 D.01完成，最新版Multiwfn可在其官網免費下載。不熟悉Multiwfn的建議參看《Multiwfn入門tips》（http://www.shanxitv.org/167）。本文的計算除了上述博文外還涉及以下博文的知識，缺乏基礎者應當先看看：
《談談隱式溶劑模型下溶解自由能和體系自由能的計算》（http://www.shanxitv.org/327）
《Gaussian中用TDDFT計算激發態和吸收、熒光、磷光光譜的方法》（http://www.shanxitv.org/314）
《談談諧振頻率校正因子》（http://www.shanxitv.org/221）

2 實例：繪制plumericin的構象平均電子光譜

這里通過plumericin（雞蛋花素）作為例子，演示怎么繪制其構象平均的UV-Vis和ECD光譜，為了簡單起見不考慮溶劑效應。plumericin有四種值得考慮的構象，分別用a,b,c,d表示，如下所示



a與b、c與d的差異在于圖右端酯基的方向，a與c、b與d的差異在于O11是在內側還是外側。

2.1 準備

對這四個構象，依次做以下計算
(1)在B3LYP/6-31G*下進行優化，得到以下計算用的結構
(2)在B3LYP/6-31G*下做振動分析，計算時帶著此級別下的ZPE校正因子0.9806，得到自由能校正量
(3)在M06-2X/def2-TZVP下做單點計算
(4)將(2)的自由能校正量與(3)的單點能相加，得到準確度較好的自由能。然后根據玻爾茲曼公式計算常溫下的分布比例
(5)用PBE1PBE/TZVP TD(nstates=20)關鍵詞在PBE0/def-TZVP級別下以TDDFT方法計算最低的20個激發態

第(5)步計算的輸出文件都已經提供在了Multiwfn壓縮包里的examples\spectra\weighted文件夾中。

計算構象分布比例需要把自由能算得非常準。當前計算級別，只能說能給出個定性正確的分布比例，但要求定量結果較好的話，應當用更高級別計算單點能，比如DSD-PBEP86-D3/def2-QZVP。如果考察溶劑下的情況，一定要帶著溶劑模型才行，否則分布比例可能嚴重錯誤。當前計算得到的a,b,c,d構象分布比率分別為60.46%、19.50%、16.86%、3.17%。因此，構象平均的UV-Vis或ECD光譜應當計算如下：
構象平均光譜=0.6046*a的光譜+0.1950*b的光譜+0.1686*c的光譜+0.0317*d的光譜

2.2 繪制構象平均的UV-Vis譜

我們寫一個叫做multiple.txt的文件（必須叫這個名字），內容如下。里面每一行是上面電子激發計算的輸出文件路徑，之后是相應構象的出現比率。
examples\spectra\weighted\a.out 0.6046
examples\spectra\weighted\b.out 0.1950
examples\spectra\weighted\c.out 0.1686
examples\spectra\weighted\d.out 0.0317

文件路徑寫相對路徑還是絕對路徑都行。這里寫的是相對路徑，是對于“當前目錄”而言的，如果不知道什么叫“當前目錄”，看《將文件快速載入Multiwfn程序的幾個技巧》（http://www.shanxitv.org/237）。注意如果你用的是Linux環境，multiple.txt里涉及的有的文件的路徑里有/符號或空格，則必須把路徑兩邊加上雙引號，否則文件沒法載入。

如手冊3.13.2節所示，Multiwfn繪制光譜絕不僅僅支持Gaussian輸出文件，比如ORCA、sTDA的輸出文件也支持，因此multiple.txt里引入的文件也可以是其它類型的，甚至可以不同來源文件混用。

啟動Multiwfn，載入multiple.txt，然后依次輸入
11  //繪制光譜的主功能
3  //UV-Vis
0  //顯示光譜

馬上看到下面的圖


圖中紅色粗曲線是權重平均后的UV-Vis光譜，可以直接和實驗對照。如圖例所示，綠、藍、粉、黑色曲線分別是a,b,c,d四種構象的光譜。從這樣的圖上，非常便于比較各個構象光譜以及權重平均光譜的差異。由于構象a的比例達到60.5%，占了絕對主導，因此權重平均的光譜整體上也和構象a的比較相似。

用Multiwfn作過單個體系的光譜的人都知道圖中黑色豎線對應的是量化程序直接算出來的激發能和強度數據，經過展寬就產生了理論光譜。當前圖上黑色豎線是所有構象的豎線的集合，但每個構象對應的豎線的高度都已經乘過了相應的權重值。因此，可以認為當前圖中紅色粗曲線是由圖上這些黑色豎線直接展寬產生的。

Multiwfn還可以把各個構象對權重平均的光譜的貢獻曲線直接繪制出來，由此便于討論實際光譜是如何構成的。關閉剛才的圖，然后選一次選項18，然后再選0重新繪制光譜，得到下圖


此圖中每個構象各自的光譜曲線高度已經乘上了相應的分布比率，因此，圖中四個構象光譜曲線高度之和就是紅色粗線；換句話說，圖中綠、藍、粉、黑線分別對應于a,b,c,d構象對權重平均光譜的貢獻。從此圖上明顯看出，a構象的貢獻是最顯著的。

此圖中豎線位置和高度和前一張圖一樣，但是此圖里不同構象對應的豎線用了不同的顏色，和曲線顏色相對應。由于此圖中無論是曲線高度還是豎線高度都已經是乘過構象權重值的，所以可以認為每個構象的光譜曲線是由相同顏色的豎線展寬產生的。

2.3 繪制構象平均的ECD譜

所有輸入文件同上一節。還是啟動Multiwfn，載入multiple.txt，然后依次輸入
11  //繪制光譜的主功能
4  //ECD
2  //載入速度表象的轉子強度
 
為了下面便于討論，我們先選16讓程序把權重平均的ECD譜的極大點和極小點位置都找出來，如下所示
Local maximum X:       181.6653      Value:        24.5974
Local maximum X:       223.9452      Value:        44.9789

Local minimum X:       158.1958      Value:        -0.0081
Local minimum X:       200.9999      Value:       -30.7966
Local minimum X:       253.6460      Value:        -5.2244

然后選0顯示光譜，看到下圖


曲線和豎線的含義和上一節相同，不再累述。前面說過圓二色譜對構象比較敏感，確實，當前的ECD圖上，各個構象光譜的差異比起它們在UV-Vis中明顯更大。相對來說，a和b的譜較相似，c和d的較相似。雖然構象a的出現比例占絕對優勢，因此其ECD譜（綠曲線）和權重平均的譜（紅曲線）相似程度最高，但是，此圖也反映，只考慮這一個構象還是不充分的，因為在<190nm的區域，如果不考慮其它構象的影響，則181.7nm處的峰就難以凸顯出來。

關閉上圖。為了繪制出各個構象對權重平均ECD譜的貢獻，選一次18，再選0，看到下圖


從這個圖上可以非常直觀地看出權重平均ECD譜的構成。181.7nm的峰很明顯是粉色曲線（構象c）所貢獻的，因為它幾乎和紅粗線重合，而構象a,b,d在此處的貢獻不僅小，而且由于貢獻有正有負，彼此也很大程度抵消了。對201.0nm、253.6nm的負峰和223.9nm正峰，構象a（綠線）顯然起到了主導性的貢獻，因為其曲線很接近紅粗線，構象b對201.0nm和223.9nm處的峰的形成也起到了推波助瀾的作用，因為其貢獻和構象a的符號是相同的。構象c（粉線）在223.9nm處搗蛋，產生顯著負貢獻，要是沒有它，此處的紅粗線的峰高度能明顯更高。

3 總結

本文通過很簡單的例子說明如何結合主流量化程序和Multiwfn繪制構象權重平均的光譜。可見過程極為簡單，寫個multiple.txt文件列出各個構象的輸入文件以及相應比例即可。而且還可以把所有構象的曲線同時繪制出來，或者把各個構象產生的貢獻繪制出來，對討論權重平均光譜的特征極為有益。Multiwfn還提供了大量選項可以調節繪圖效果，可以得到能直接用于發表的圖。

另外，在繪圖界面里用選項2，可以導出spectrum_curve.txt和curveall.txt到當前目錄下。前者記錄的是權重平均的光譜曲線數據，而curveall.txt記錄的是各個構象的曲線數據，每個構象對應一列。這些文件可以導入Origin等程序來作圖，屆時可以更靈活地調節圖像效果。

在繪圖界面里還有選項21 Set status of showing weighted curve and curves of individual systems，進入其中后可以選擇是只顯示構象權重平均的曲線、只顯示每個體系各自的曲線，還是二者同時顯示。

本文的體系需要考慮的構象比較少，手工計算不麻煩，但如果體系可能有很多構象，自己懶得去一一搜索，那么最理想的做法是使用Molclus程序做構象搜索，參見《使用molclus程序做團簇構型搜索和分子構象搜索》（http://bbs.keinsci.com/thread-577-1-1.html），使用其中的gentor組件可以很方便地大批量產生初猜構象，見《gentor：掃描方式做分子構象搜索的便捷工具》（http://bbs.keinsci.com/thread-2388-1-1.html）。將gentor+molclus+量化程序結合使用，你會發現超級靈活、超級便捷，而且希望精度高一些（但計算耗時高）還是希望速度快一些（計算精度低），完全可以通過調節所調用的量化程序的關鍵詞來自由控制，用熟了之后就基本再也不想用其它構象搜索程序了。

附：同時繪制多個體系的光譜

在2017-Aug-14及之后更新的Multiwfn當中也可以同時直接繪制多個體系的光譜，比如可以把不同分子、不同級別下的光譜繪制在一起。做法很簡單，把前文的multiple.txt里的權重值改成圖例文字即可，這樣就不會產生和顯示總光譜，而且圖例可以自定義。比如multiple.txt可以寫成這樣
examples\spectra\indigo\ZINDO.out ZINDO
examples\spectra\indigo\TD-PBE0.out TD-PBE0/6-31G*
examples\spectra\indigo\TD-PBE0_TZVP.out TD-PBE0/def-TZVP
結果如下，同時把三種不同計算級別的譜圖顯示出來了，便于對比。完整的繪制的例子見Multiwfn手冊4.11.6節。