通過鍵級曲線和ELF/LOL/RDG等值面動畫研究化學反應過程
通過鍵級曲線和ELF/LOL/RDG等值面動畫研究化學反應過程
文/Sobereva @北京科音 2013-Sep-11
很久很久以前，寡人寫過一篇帖子名為《制作動畫分析電子結構特征》（http://www.shanxitv.org/86）。在此帖中筆者介紹了怎么制作動畫研究IRC、SCAN過程中電子結構的變化。本文算是對那篇文章的進一步擴展。這里將介紹怎么利用Multiwfn (http://www.shanxitv.org/multiwfn)繪制出IRC過程中鍵級變化的曲線圖，如何做出ELF/LOL/RDG等值面變化的動畫。
在《制作動畫分析電子結構特征》當中筆者獲取IRC/SCAN過程中的波函數文件是通過自寫Gaussian的link實現的。那個方法雖然不用每一步重新做一次單點計算，但是對于初學者來說實現起來比較困難。本文中獲得IRC/SCAN過程中每個點的.wfn/.fch文件的方法是使用IRCsplit和SCANplit程序實現的，這對于初學者來說明顯更容易實現（盡管需要對每個點重新做一次單點任務因而需要額外的耗時），見《產生Gaussian的IRC和SCAN任務每個點的波函數文件的工具》（http://www.shanxitv.org/199），操作方法在此帖已有介紹故在本文不再詳述。
本文使用Multiwfn 3.2.1(dev)版計算（不要用老版本），VMD1.9生成等值面圖，G09生成波函數。Linux系統為RHEL6-U1 64bit。若未注明，操作皆在Win7 64bit下進行。本文用到Multiwfn的silent運行模式，見Multiwfn手冊5.2節的介紹。也用到DOS下的批處理腳本，可以參考《從高斯windows下的批量執行談dos批處理文件》（http://www.shanxitv.org/6）一文的對DOS批處理文件的簡介。
本文使用乙炔三聚化過程的IRC作為例子，對應的IRC動畫如下所示。
PS：去過Multiwfn 2013年8月的培訓班的人應該記得筆者留了個作業，就是通過各種方法研究這個三聚化過程中電子結構的變化，本文一定程度上也算是對這個作業的部分解答。
1 產生三聚化過程IRC路徑中每個點的fch文件
IRCsplit程序包中的examples\trimerization.out就是這個三聚化過程的IRC任務的輸出文件，計算由G09在B3LYP/6-31G**下完成，IRC向前和向后的路徑分別有34和25個點，加上過渡態的點整個路徑有60個點。IRC路徑兩端皆已延伸至反應物和產物極小點，分別對應三個乙炔和苯分子。利用IRCsplit程序結合標準單點任務文件examples\trimerization_SP.gjf，就可以將trimerization.out中每個點的坐標提取并生成一批新的Gaussian單點輸入文件，批量運行之并經過formchk批量轉換之后就可以得到每個點的.fch文件。假設這些.fch文件都放在了c:\IRC目錄下，文件名為IRC0001.fch、IRC0002.fch...IRC0060.fch。
注1：由于計算Mayer鍵級、多中心鍵級都需要基函數信息，而.wfn文件不包含這樣的信息，所以這里生成.fch文件而不生成.wfn文件
注2：用IRCsplit處理之前，應當把trimerization_SP.gjf中補上nosymm關鍵詞。這樣新生成的Gaussian輸入文件也就都帶著nosymm關鍵詞了。如果沒有nosymm關鍵詞，Gaussian做單點計算時由于會自動調整分子朝向到標準朝向，會導致某些點產生的.chk/.wfn文件中體系的平面和其它點不對應，隨后做等值面變化動畫的時候體系就會中途發生奇怪的翻轉。如果不打算做等值面變化的動畫的話則無需nosymm，寫了它之后做單點時由于不啟用對稱性，計算還會多耗時很多。
2 繪制Mayer鍵級變化曲線
我們先繪制Mayer鍵級變化曲線。由于Multiwfn默認只輸出>0.05的Mayer鍵級，所以需要修改settings.ini文件，將bndordthres參數改為0.0，之后所有鍵級就都會輸出了。Multiwfn中計算Mayer鍵級的過程是輸入9然后輸入1，為了批量計算，我們寫個文本文件（假設叫batch.txt）放在Multiwfn所在目錄下，文件內容只有兩行，即為
9
1
然后寫個文本文件（假設叫batchrun.bat，即批處理文件）放在Multiwfn目錄下，內容為
for /f %%i in ('dir c:\IRC\*.fch /b') do Multiwfn c:\IRC\%%i < batch.txt > c:\IRC\%%~ni.txt
之后，雙擊batchrun.bat圖標，Multiwfn就會依次使用那些.fch文件計算Mayer鍵級，輸出信息依次產生在c:\IRC目錄下的IRC0001.txt、IRC0002.txt...IRC0060.txt里面。
將這些輸出文件拷貝到Linux系統某文件夾下，并在這個文件夾下運行
grep "1(C ) 6(C )" * > out.txt
所有IRC點上的C1-C6鍵級就都存到out.txt了。用Ultraedit的列模式（或者諸如Linux的awk命令）只把out.txt的最后一列提取出來，然后導入到諸如origin之類作圖程序里。然后我們從trimerization.out的末尾可以看到這樣的信息
Summary of reaction path following
--------------------------------------------------------------------------
Energy Rx Coord
1 -0.05598 -2.49599
2 -0.05520 -2.39872
3 -0.05432 -2.29904
...
59 -0.33639 3.13057
60 -0.33639 3.13555
其中Energy就是相對于過渡態能量的每個點的能量，單位為a.u.。Rx Coord是反應坐標。用Ultraedit的列模式把Rx Coord那一列的數據也提取出來導入到origin里。之后，在origin里繪制曲線圖，X軸就是Rx Coord值，Y軸就是鍵級數據，就得到了下圖
橫坐標零點對應于過渡態。可見，隨著反應的進行，鍵級由3.0左右不斷降低到1.5左右，體現出乙炔的三重鍵逐漸變成苯上的C-C間的1.5重鍵。而且反應過程一開始（乙炔相互接近并略微變形）和最后（氫原子位置調整）過程中鍵級變化不大，在反應中途明顯涉及到新C-C鍵形成舊C-C鍵削弱時鍵級變化才迅速。
通過執行grep "1(C ) 2(C )" * > out.txt以及grep "1(C ) 7(H )" * > out.txt，我們以類似方法把C1-C2和C1-H7鍵級變化也提取出來并一起繪制到圖上。為了便于分析，能量變化曲線也一起畫上去。如下所示
可見原本沒成鍵的C1-C2和原先的三重鍵C1-C6最后鍵級都變得一致了，約1.5。由于C-H鍵在反應過程中沒有直接受到影響，所以鍵級一直基本約為1.0。
我們還可以通過比如grep "Atom 1(C )" * > out.txt將C1的原子價導出并用于作圖。對于這樣的閉殼層體系，原子價就是相應原子涉及的所有鍵的鍵級加和。我們可以看到盡管每個C-C鍵的鍵級都變化了很大，但碳的原子價一直接近于4.0，沒有明顯波動，表現了反應過程中總鍵級守恒性，也就是鍵級會此消彼長，一個鍵的增強伴隨著另一個鍵的等量削弱。
有興趣的讀者建議也自行嘗試繪制下拉普拉斯鍵級(LBO)、離域化指數、各種原子電荷、AIM鍵臨界性質等屬性隨著反應進行的變化曲線。
3 繪制多中心鍵級變化曲線
我們還可以研究乙炔三聚化過程的碳原子六中心鍵級的變化。把前述的batch.txt文件內容改為
9
2
1,2,3,4,5,6
然后還是運行batchrun.bat，并且把生成的.txt文件都弄到Linux某文件夾里，然后進入此文件夾并運行
grep "The bond order is" * > out.txt
提取out.txt中的最后一列六中心鍵級數據，并結合能量變化一起作圖，如下所示
顯然乙炔三聚體肯定沒有6中心鍵，所以一開始六中心鍵級為0。形成苯之后六中心鍵級達到了最大值，因此苯環上的電子很容易整體離域，這也是為什么苯具有極強的芳香性。在三聚化的過程中六中心鍵級的變化比較奇特，并不是單調增加的，而是在過渡態稍靠后一些的位置上有個局部極大點，其原因何在？對于平面體系，多中心鍵級可以精確分離為sigma和pi多中心鍵級，我們不妨把它們分別繪制出來，看看sigma電子和pi電子在反應過程中是如何影響總多中心鍵級的。
我們先計算sigma多中心鍵級。這個過程需要利用Multiwfn主功能100里的子功能22，它會自動將平面體系的pi軌道找出來，選擇把pi軌道占據數都設成0，之后再算多中心鍵級，得到的就是sigma多中心鍵級了。具體來說還是修改batch.txt文件，將內容改為下面這樣
100
22
0
1
0
9
2
1,2,3,4,5,6
然后再次按照如上過程運行batchrun.bat并提取數據，就得到了sigma多中心鍵級。
接下來算pi多中心鍵級，操作和上面完全一樣，只不過用到的batch.txt文件內容改為下面這樣，計算多中心鍵級前就會把所有pi軌道以外的軌道的占據數都設成0來扣除它們的貢獻
100
22
0
2
0
9
2
1,2,3,4,5,6
把sigma和pi多中心鍵級和之前的圖作到一起，如下所示
可以看到，之所以在過渡態稍靠后一些的位置六中心鍵級有個局部極大點，是因為sigma電子造成的。然后隨著反應的繼續進行，sigma電子的離域性迅速衰退，在形成苯之后sigma離域性幾乎可忽略不計，因此苯沒有sigma芳香性。pi電子的離域性則是隨著反應的進行單調增強的，但是增強基本是在過渡態以后才比較明顯，到了苯的基本結構快要形成階段增強得極快，顯示出苯的pi芳香性的迅速形成。
本節研究結論和ELF、ELF-pi、ELF-sigma對此問題的研究結論定性一致，如下所示，有興趣者可以參看Theoretical Aspects of Chemical Reactivity的第五章。用多中心鍵級來研究芳香性比起用ELF方便、省時得多，而且比較上圖和下圖也可以看出多中心鍵級的結果更為清楚，ELF-sigma的曲線會讓不熟悉ELF的人會誤以為即便形成了苯之后仍然有很強的sigma芳香性。多中心鍵級是筆者最為推薦的研究芳香性的方法。
另外，有興趣的讀者也建議通過Multiwfn來計算其它各種芳香性指標在這個三聚化過程中的變化，Multiwfn能算的芳香性指標多達十余種，詳見《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）。
4 制作ELF等值面圖變化的動畫
本節將要繪制乙炔三聚化過程中電子定域化函數(ELF)等值面圖變化的動畫。制作ELF平面圖的動畫在《制作動畫分析電子結構特征》中已經詳細介紹了，這里就不再談怎么做了。制作等值面圖的動畫很重要，可以讓函數在三維空間中的變化看得很清楚。對于當前的例子，既涉及到sigma電子又涉及到pi電子的變化，如果作平面圖就得做至少兩個平面上的動畫才行，既麻煩也不直觀。
首先需要生成每個點對應的ELF的cube文件。把前述的batch.txt改為以下內容
5
9
2
2
再將前述的batchrun.bat批處理文件改為如下內容，
for /f %%i in ('dir c:\IRC\*.fch /b') do (
Multiwfn c:\IRC\%%i < batch.txt
rename ELF.cub %%~ni.cub
)
然后雙擊運行此文件，將會在當前目錄下產生IRC0001.cub、IRC0002.cub...IRC0060.cub。
這里用VMD程序來產生每個cube文件對應的等值面圖像文件。把剛得到的cube文件都移動到VMD所在目錄下面。在VMD目錄下寫個文本文件isoall.tcl，內容如下
set isoval 0.8
axes location Off
for {set i 1} {$i<=60} {incr i} {
set name IRC[format %04d $i]
puts "Processing $name.cub..."
mol default style CPK
mol new $name.cub
scale to 0.3
rotate x by -30
translate by 0.000000 0.20000 0.000000
mol addrep top
mol modstyle 1 top Isosurface $isoval 0 0 0 1 1
mol modcolor 1 top ColorID 8
render snapshot $name.bmp
mol delete top
}
腳本開頭isoval用來設定等值面的數值，這里對應于繪制ELF=0.8的等值面。for {set i 1} {$i<=60} {incr i}代表這里要處理編號為0001到0060的cube文件。腳本中scale to、rotate x by、translate by都是用來調節視角的。對于其它體系應根據實際情況對參數進行適當調整。
啟動VMD后在VMD的控制臺里面運行source isoall.tcl，在VMD目錄下就會輸出IRC0001.bmp、IRC0002.bmp...IRC0060.bmp。圖像尺寸和當前VMD窗口大小是一致的。
然后隨便找個視頻制作軟件將這些圖片合成動畫就行了，諸如會聲會影。為了便于放在網頁上，這里通過Linux下的convert工具制作成.gif動畫。也就是把所有生成的bmp文件拷到Linux下某文件夾，然后進入此文件夾并運行
convert -resize 300 -delay 12 -colors 30 -monitor *.bmp ELF_IRC.gif
很快就得到了動畫文件ELF_IRC.gif，如下所示。convert命令的使用詳見《制作動畫分析電子結構特征》里的介紹。
從動畫中可以看到一開始乙炔上的三重鍵對應于血小板形狀的高定域性區域，隨著反應的進行，血小板形狀逐漸收縮成啞鈴狀（體現出sigma鍵+部分pi鍵特征），同時相鄰的乙炔間的未成鍵C-C之間也逐漸出現了同樣的啞鈴狀等值面。這個動畫清楚地展現出反應過程中化學鍵是如何形成、變化的。
5 制作LOL等值面圖變化的動畫
制作LOL等值面變化的動畫和上一節的過程幾乎沒有差異，把所用的腳本稍微變一下就行了。
batch.txt改為以下內容
5
10
2
2
batchrun.bat批處理文件改為如下內容，
for /f %%i in ('dir c:\IRC\*.fch /b') do (
Multiwfn c:\IRC\%%i < batch.txt
rename LOL.cub %%~ni.cub
)
isoall.tcl改為如下內容
color Display Background white
set isoval 0.6
axes location Off
for {set i 1} {$i<=60} {incr i} {
set name IRC[format %04d $i]
puts "Processing $name.cub..."
mol default style CPK
mol new $name.cub
scale to 0.4
rotate x by -30
translate by 0.000000 0.20000 0.000000
mol modstyle 0 top CPK 0.500000 0.300000 12.000000 10.000000
mol addrep top
mol modstyle 1 top Isosurface $isoval 0 0 0 1 1
mol modcolor 1 top ColorID 3
render snapshot $name.bmp
mol delete top
}
最終得到的動畫如下所示。表現的電子結構的變化和ELF是基本一致的。
作為練習，建議大家嘗試繪制電子密度拉普拉斯函數、ELF-pi、變形密度等實空間函數在此IRC過程中的變化，方法是大同小異的，稍微改下腳本就能實現。只要能作出一個波函數文件對應的圖，就肯定能直接作出整個IRC過程的變化圖。
6 繪制RDG填色等值面變化
約化密度梯度(RDG)填色等值面圖對于研究弱相互作用非常有用而且直觀，在《使用Multiwfn圖形化研究弱相互作用》（http://www.shanxitv.org/68）和《使用Multiwfn研究分子動力學中的弱相互作用》（http://www.shanxitv.org/186）當中有十分詳細的討論。雖然化學反應主要是強相互作用范疇，但是研究這個過程中的弱相互作用也是有益的。
制作RDG填色等值面變化的動畫和前面兩節的流程也是完全一致，區別僅在于所用腳本。batch.txt的內容應改為
100
1
15,13
-10
2
2
3
batchrun.bat批處理文件改為如下內容
for /f %%i in ('dir c:\IRC\*.fch /b') do (
Multiwfn c:\IRC\%%i < batch.txt
rename func1.cub f1_%%~ni.cub
rename func2.cub f2_%%~ni.cub
)
isoall.tcl改為如下內容
color Display Background white
color scale method BGR
set isoval 0.5
axes location Off
for {set i 1} {$i<=60} {incr i} {
set name IRC[format %04d $i]
puts "Processing f1_$name.cub and f2_$name.cub..."
mol default style CPK
mol new f1_$name.cub
mol addfile f2_$name.cub
scale to 0.4
rotate x by -30
translate by 0.000000 0.20000 0.000000
mol modstyle 0 top CPK 0.500000 0.300000 12.000000 10.000000
mol addrep top
mol modstyle 1 top Isosurface $isoval 1 0 0 1 1
mol modcolor 1 top Volume 0
mol scaleminmax top 1 -0.04 0.02
render snapshot $name.bmp
mol delete top
}
做出的動畫如下
一開始相鄰乙炔之間以及三個乙炔中央有綠色RDG等值面，體現的是它們之間的范德華相互作用。當反應開始，三者相互接近，等值面由綠變藍，說明分子間開始呈現較強吸引作用。并且中間的小等值面由綠變棕，說明三個乙炔間開始出現了一定位阻，這對反應的進行會有所阻礙。由于默認情況下，Multiwfn計算RDG的時候會屏蔽掉密度較高的區域，所以當相鄰的乙炔間開始形成C-C鍵故而電子密度驟增后，相應的RDG等值面在動畫中就消失了。但是在體系中央的那個RDG等值面依然存在，且越變越紅，說明形成苯環過程中環的位阻是愈發增大的。
7 繪制RDG填色等值面+ELF等值面的變化
一般認為RDG主要用處是表現弱相互作用，而ELF又只能表現強相互作用，在JCTC,8,3993(2012)中作者提出可以將這兩種函數的等值面同時繪制在一起，這樣圖中就把強、弱相互作用一起表現了，就顯得更完整了。這里我們也繪制這樣的RDG+ELF等值面的動畫。
把本文第4節產生的ELF格點數據IRC0001.cub、IRC0002.cub...IRC0060.cub，以及我們剛才產生的名字為f1_xxxx.cub、f2_xxxx.cub這些格點數據一起放到VMD目錄下。然后將isoall.tcl改為如下內容
color Display Background white
color scale method BGR
set isovalRDG 0.5
set isovalELF 0.8
axes location Off
for {set i 1} {$i<=60} {incr i} {
set name IRC[format %04d $i]
puts "Processing $name.cub, f1_$name.cub and f2_$name.cub..."
mol default style CPK
mol new f1_$name.cub
mol addfile f2_$name.cub
mol modstyle 0 top CPK 0.500000 0.300000 12.000000 10.000000
mol addrep top
mol modstyle 1 top Isosurface $isovalRDG 1 0 0 1 1
mol modcolor 1 top Volume 0
mol scaleminmax top 1 -0.04 0.02
mol new $name.cub
mol modstyle 0 top Isosurface $isovalELF 0 0 0 1 1
mol modcolor 0 top ColorID 8
scale to 0.3
rotate x by -30
translate by 0.000000 0.20000 0.000000
render snapshot $name.bmp
mol delete top
mol delete top
}
做出的動畫如下
可見RDG和ELF等值面有機地融合到了同一個畫面上。當描述弱相互作用的RDG等值面消失后，馬上ELF等值面就出現了，很好地表現出弱相互作用到強相互作用的無縫轉換。這樣的動畫特別是對于研究化學反應過程明顯涉及弱相互作用的情況十分有幫助，比如反應物或產物結構是復合物、反應過程涉及到氫鍵形成與斷裂。
8 總結
本文以乙炔三聚化過程為例通過繪制鍵級變化曲線、制作ELF/LOL/RDG等值面動畫展現了怎么研究IRC過程中電子結構變化。SCAN過程、幾何優化過程、分子振動過程、從頭算動力學過程等等動態過程中電子結構、分子屬性或分子間相互作用等方面的變化也都可以對相應的性質繪圖來表現。無論是把曲線圖放在正文中還是把動畫放在補充材料里，都可以使研究內容更充實、細致、深入、直觀易懂。對于量子化學的外行、本科生們，通過這樣的圖形描述，也可以讓他們迅速、清楚地了解到究竟化學過程中發生了什么，比起一堆抽象的數據和語言描述也有趣得多。
作這些圖其實在技術上都沒有什么難度，根據以上例子和《制作動畫分析電子結構特征》里的例子，在理解每一步操作的基礎上舉一反三，多查閱資料即可。如果遇到實在解決不了的技術問題也可以直接寫郵件給筆者(sobereva@sina.com)。另外需要提醒的是，對于實際研究的體系作圖時切勿直接照搬本文提供的VMD腳本，很多作圖設定需要自己修改才能達到滿意效果，所以應該弄懂腳本的含義。
Multiwfn能計算的性質、繪制的圖數不勝數，諸如電荷、原子多極矩、自旋布居、軌道成份、AIM盆或模糊空間的定域/離域化指數、分子表面定量性質、臨界點位置和性質、鍵徑、電子密度圖、能量密度圖、分子/原子體積、片段差值圖等等等等，他們在化學過程中的變化都很值得分析，必定能由此得到很有意義的發現并帶來新觀點。希望本文的例子能給讀者帶來一些啟發。研究化學反應不應局限在算個渡態結構和勢壘，然后頂多再算個IRC就了事的程度，實際上還有無數有用的信息有待進一步挖掘、還有豐富多彩的展現化學反應過程的角度。