制作動畫分析電子結構特征

制作動畫分析電子結構特征

文/Sobereva @北京科音
First release: 2011-May-1    Last update: 2014-Apr-15

在《自寫Link生成Gaussian的IRC任務中每個點的波函數文件》（http://www.shanxitv.org/85）一文中筆者已經談了怎么獲得IRC、Scan等任務中每一步的波函數文件，本文將通過實例介紹如何利用Multiwfn對它們批量作圖并生成動畫。動畫是圖形化研究電子結構特點的最高級形式，比靜態圖像能展現更豐富的信息，即便是外行人對此也是喜聞樂見的，插在幻燈片里會使得講述妙趣橫生。注：如果限于水平搞不定自寫link的話，可以參看另一種方法，更容易實現，見《通過鍵級曲線和ELF/LOL/RDG等值面動畫研究化學反應過程》（http://www.shanxitv.org/200）。

本文使用的是Multiwfn 2.0.2版，可以在http://www.shanxitv.org/multiwfn免費下載。本文涉及的Gaussian輸入文件、波函數文件、動畫文件等等都可以在這里下載：/usr/uploads/file/20150610/20150610200942_56003.rar。筆者windows系統為WinXP，Linux為RHEL6 server，ImageMagick版本為6.5.4-7。

 

1 三重態H2拉伸過程的LOL函數填色圖動畫

這是最簡單的一個例子。定域化軌道指示函數(LOL)類似電子定域化函數(ELF)，數值在[0,1]之間，數值越大的地方說明電子定域性越大，即電子的運動被越大程度地限制在相應區域內。當兩個擁有相同自旋電子的氫原子間距逐漸縮小，由于Pauli互斥，必然導致兩個電子運動區域相互躲著，高定域性區域分別跑到氫分子的兩端去。此例就來繪制這樣一個過程的動畫，看看是不是這樣。

在Gaussian中通過scan關鍵詞，讓原子間距從0.2開始拉，共10步，每步0.1埃。輸入文件如下
%subst l555 .
#P ub3lyp/cc-pvdz scan geom=nocrowd 6d 10f extralinks=l555

Title Card Required

0 3
 H             
 H                  1            B1

   B1             0.20000000 10 0.1

/sob/H2tri/H2.wfn
其中l555是自編的link，在《自寫Link生成Gaussian的IRC任務中每個點的波函數文件》里有介紹，extralinks=l555使得l555被插入在SCF模塊l502之后，即每一步都輸出波函數文件，將依次得到H2_0000.wfn, H2_0001.wfn...H2_0010.wfn。geom=nocrowd是為了避免距離檢測，否則原子距離小于0.5埃時Gaussian會報錯。

運行完畢后將裝著波函數文件的H2tri目錄從linux中拷到windows，放在Multiwfn所在目錄下。將Multiwfn的Settings.ini里的isilent參數設為0（后文皆如此）。在Multiwfn目錄下編寫一個文本文件batchcolor.txt，內容是（//后是注釋）：
4   //繪制平面圖
10  //LOL函數
1   //填色圖
200,200   //兩方向格點數
3   //YZ平面。波函數文件中分子軸順著Z
0   //X=0
4   //顯示原子標簽
3   //用藍色字體
0   //保存圖像到當前目錄
[空行]
[空行]
若不知道上述過程是什么目的，請看Multiwfn手冊附錄部分關于silent模式運行的介紹。

然后再編輯一個文本文件也放在Multiwfn錄下，叫batchrun.bat，內容只有一行：
for /f %%i in ('dir H2tri\*.wfn /b') do Multiwfn H2tri\%%i < batchcolor.txt
這是一個DOS的批處理腳本。其中dir H2tri\*.wfn /b代表輸出H2tri目錄下所有.wfn文件名。這個腳本會使每個.wfn文件都被Multiwfn安靜地執行一遍，相當于依次執行
Multiwfn H2tri\H2_0000.wfn < batchcolor.txt
Multiwfn H2tri\H2_0001.wfn < batchcolor.txt
...
Multiwfn H2tri\H2_0010.wfn < batchcolor.txt
由于將isilent設為了1，圖像不會自動彈出，而且每一步要輸入的命令都通過batchcolor.txt重定向進了Multiwfn，所以每一步Multiwfn都會無聲無息地執行，用戶不需要干預。執行batchrun.bat后，由于體系小，幾秒鐘后當前目錄下就出現了DISLIN.PNG, DISLIN_1.PNG, DISLIN_2.PNG...DISLI_10.PNG，依次對應核間距從近到遠。

PS：有些人希望Multiwfn變成全GUI的，我之所以不這樣做，而是將之開發成現在這樣GUI+文本模式，一方面是因為那樣并不會使操作簡化，完全是多余的，另外屆時也沒法像上文這樣，通過腳本來方便地批量執行任務。

接下來就是將這些圖像文件合成為動畫，我傾向于用gif格式，可以方便地嵌入到網頁中，圖像清楚，盡管色彩數上限為256，但足夠了。合成gif的軟件非常多，在Windows下比如可以用Atani，但是最方便的還是用Linux下的convert命令，它是ImageMagick圖形軟件包中的一員，通常在Linux發行版本中都自帶也默認安裝。先將剛才得到的所有.PNG文件拷貝到/sob/H2tri目錄下，然后進入/sob目錄，在shell下執行（這是一整行命令）：

convert -delay 20 -resize 500 -colors 128 -monitor -reverse -fill red -annotate +5+395 "Created by Sobereva" H2tripic/DISLIN.PNG H2tripic/DISLIN_?.PNG H2tripic/DISLI_??.PNG H2tripic/H2tri_all.gif

這條命令表明H2tripic下的DISLIN.PNG、DISLIN_?.PNG、DISLI_??.PNG（問號代表通配符）按順序依次被convert命令讀入，然后組合成H2tri_all.gif。不能為了省事直接寫DISLI*.PNG，這樣的話輸入文件順序是錯的。文件名前面那一串是對圖像的操作符或者設定參數的命令，delay指定每一幀停留的時間，越大則播放越慢。resize是將圖像寬度調整為500像素而橫縱比不變，因為Multiwfn默認輸出的圖像尺寸太大，不這樣做則生成的.gif體積會較大，不適合放到網頁，如果用于插入幻燈片倒是可以設大些以保證圖像清晰。colors設定此gif動畫的色彩數，將默認的256色降為128色不影響效果，文件體積也會減小。monitor代表監控convert命令內部執行過程。reverse代表將讀入的圖像順序反轉，因為上面的scan過程是由近及遠，而此例想讓過程由遠及近。annotate命令用來在動畫上寫上文字，便于保護版權或進行內容說明，+5+395代表字符串的左下角相對于圖像左上角的位置。前面的-fill red說明文字用紅色。

不一會兒，H2tri_all.gif就生成了。如圖：

可見，兩原子離得越近三重態H2兩端定域化區域越顯著，和預期的一致。然而，仔細分析會發現這個過程還給出了新的信息，也就是當兩個原子核離得很近時，二者中央很小區域內定域性也明顯增加了！這看似與Pauli互斥有些違背，但實際上這并不奇怪。因為原子核離得越近，中間區域的電子感受到的核吸引勢越強，自然定域性會升高。如果進一步離近，則兩個氫核就融合成He了，與三重態He的LOL圖相對比，會發現動畫的最后一幀（距離0.2埃）與之十分相似。這樣來看，這個動畫實際上蘊含了兩個截然不同的過程：動畫前一半，兩個氫核越來越近，兩個自旋相同的電子因Pauli互斥相互避開，導致H2兩端的高定域化區域膨脹。動畫后一半，三重態H2逐漸顯示三重態He特征，中間緊湊的高定域性區域對應于s軌道的電子，而兩端高定域性區域實際上是占據p軌道的電子。這兩個過程雖本質不同，但在動畫中看上去是平滑過渡的。

下面的動畫是上面這個scan過程的電子密度差圖，從中可以得到與上述一致的結論，而且顯得更清晰。這個圖的做法在看過下一節的介紹后就會明白。圖中虛線代表相對于兩原子都在孤立狀態時電子密度減小的區域，實線代表增加的區域。可見，一開始原子離得遠，相互作用很弱，對原始電子密度幾乎沒有影響，幾乎沒出現等值線。隨著它們離近，由于Pauli互斥導致中間電子密度大幅減少，電子紛紛轉移向兩端了。離近到一定距離后，中間一小部分電子密度開始增加，說明氦核開始形成了，s電子開始聚集了。僅僅是一個H2拉伸過程的動畫就包含了很多信息，可見通過動畫展現電子結構特征的威力之大！

注意彩圖gif動畫比較占空間，上面的填色圖動畫有11幀，占了535KB，而等值線動畫有25幀，才238KB。

2 單重態H2拉伸過程的電子密度差等值線動畫

做電子密度差圖必須先獲得涉及到的各個元素原子的波函數文件。每次執行Multiwfn時程序會從當前目錄下atomwfn目錄里尋找是否已有原子的波函數文件，沒有則需要調用Gaussian重算。為了方便、省時，這里先事先算好氫原子的波函數文件，放在atomwfn目錄下，并命名為H .wfn，注意空格必須有，否則Multiwfn不認。基組建議與計算分子時一致（若基組級別已較高，不一致時問題也不大）。其實計算密度差圖還有個既重要又麻煩的球對稱化原子密度問題，由于氫的密度本來就是球對稱的，這里就不提了。

執行壓縮包內的H2sin.gjf，這是單重態H2的scan任務，從0.2埃開始24步，每步0.1埃。注意盡管是單重態，但是拉遠后對稱破缺態是基態，所以要用非限制性方法，且需要加上guess(always,mix)打破對稱性。計算完畢后將儲存波函數文件的目錄H2sin拷到Multiwfn所在目錄下。

編寫一個文本文件batchcontourfix.txt放在Multiwfn目錄下，內容是：
4
-2
1
2
200,200
6  //自行輸入格點設定
0.00000  -4.50000  -6.95664  //原點
0.00000   0.04500   0.00000  //平移向量1
0.00000   0.00000   0.06957  //平移向量2
0
[空行]
[空行]
其中-2代表要做變形屬性，由于函數選了1，計算的就是變形密度（密度差）。第一節的例子中坐標軸尺寸和刻度在不斷變動，這是為了使當前幀的結構能充滿圖像，然而有些人不喜歡這樣，因此此例中自行設定原點和平移向量，坐標軸就固定住了。顯然很難直觀判斷數值應該設多少，我建議用體系涉及空間尺度最廣的一幀作參考獲取這些信息。比如此例H2_0024.wfn對應拉得最遠的一幀，就先手動對這個波函數的任意一個函數做圖（X=0的YZ平面，延展距離用默認的4.5 bohr），Multiwfn會輸出作圖時所用的格點信息：
Origin is   0.00000  -4.50000  -6.95664 bohr
The x/y/z of transition vector 1 is   0.00000   0.04500   0.00000 bohr
The x/y/z of transition vector 2 is   0.00000   0.00000   0.06957 bohr
令所有幀作圖時都使用這個格點設定，就能保證每一幀感興趣的區域都不出界。

然后，將前面batchrun.bat里的H2tri替換為H2sin，batchcolor.txt替換為batchcontourfix.txt，執行之，就得到了每幀的.PNG圖像。將圖像復制到/sob/H2sinpic目錄下，并進入/sob目錄，執行:

convert -delay 20 -resize 500 -shave '100x0' -colors 128 -monitor -reverse -fill red -annotate +5+395 "Created by Sobereva" H2sinpic/DISLIN.PNG H2sinpic/DISLIN_?.PNG H2sinpic/DISLI_??.PNG H2sinpic/H2sin_all.gif

和第一節的命令沒什么區別，只是目錄換了，并且加入了一個新的設定-shave '100x0'，這是因為這次的圖像的坐標軸橫軸固定了，左右的白邊長度也固定了，將每幀在resize為寬度為500后，左右都截掉100像素可以讓動畫的畫布更緊湊。最終所得動畫H2sin_all.gif如圖所示：

隨距離拉近后，電子密度迅速在原子間聚集，表明共價鍵正在加強。聚集的區域一開始逐漸變大并逐漸變得扁長，但后來卻開始收縮并變得接近于球型，最后一幀電子密度減少區域和增加區域都幾乎變成了球對稱狀，前者包著后者。這是因為原子序數越大，核吸引勢越強，電子密度分布范圍收縮得越厲害，基態He的1s電子密度延展范圍小于基態H的，所以當兩個氫核距離接近到氦核特征明顯出現時內部電子密度就增加了，相應地外側減少。

3 HCN異構化過程的LOL等值線動畫

如何獲得IRC過程中每個點的波函數文件的方法在《自寫Link生成Gaussian的IRC任務中每個點的波函數文件》已經介紹，就不累述，此例的Gaussian輸入文件是附件里的HCN_IRC.gjf，加上過渡態最終將得到21個波函數文件。將它們都拷入Multiwfn目錄下HCN子目錄內，編寫一個文本文件batchcontourfixldset.txt放在Multiwfn目錄下，內容是：
4
10
2
200,200
6
-6.47878  -5.57534   0.00000
0.05565   0.00000   0.00000
0.00000   0.05626   0.00000
3   //進入等值線設定
7   //從外部文件載入等值線設定
ELFctrsp.txt  //記錄著適合繪制LOL、ELF函數的等值線設定的文件（假設在當前目錄），從0至1，步長0.05。
10  // 加粗等值線
1   // 被加粗的只有一條
11  // 加粗第11條等值線，對應0.5，它所包圍的區域一般是化學上感興趣的區域
1   // 保存等值線設定
0
[空行]
[空行]
將batchrun.bat里的路徑都改為HCN，重定向文件名也改為batchcontourfixldset.txt，執行此文件。將所得的DISLIN.PNG、DISLIN_1.PNG...DISLI_10.PNG放到/sob/HCNpic/b，將DISLIN_11.PNG...DISLI_20.PNG放到/sob/HCNpic/f，這兩批圖像分別對應IRC的兩個方向（DISLIN.PNG對應TS）。進入/sob目錄，然后依次執行下述命令：

convert -delay 25 -resize 500 -colors 128 -reverse -fill red -annotate +5+395 "Created by Sobereva" -monitor HCNpic/b/DISLIN.PNG HCNpic/b/DISLIN_?.PNG HCNpic/b/DISLI_??.PNG HCNpic/b.gif

convert -delay 25 -resize 500 -colors 128 -fill red -annotate +5+395 "Created by Sobereva" -monitor HCNpic/f/DISLI_??.PNG HCNpic/f.gif

convert HCNpic/b.gif HCNpic/f.gif HCNpic/HCN_all.gif
最后一步用來將前后兩個方向的動畫合并到一起，得到的HCN_all.gif就是我們想要的：

從加粗的等值線上能很容易考察電子結構特征的變化。在接近兩種穩定異構體的構型下C-N鍵的定域化域扁平而中間凹陷，顯示出明顯多重鍵特征。氫轉移開始后，隨著氫遠離氮，氮的孤對電子逐漸顯現，并且氫開始與C-N鍵的定域化域連通，說明此時形成了三中心鍵。當氫與氮的孤對電子定域化域分離后，馬上就與碳的孤對電子區域相連，隨后碳的孤對電子逐漸消失，氫也不再參與三中心鍵，最終完成了整個異構化過程。

4 參考點在CH3Br的C-Br鍵徑上移動過程中Fermi穴函數的變化

這是個特殊的例子，不是使用多個波函數，而是使用一個波函數但是通過變化參數來獲得動畫。

Fermi穴函數是個六維量，需選擇一個參考點才能用圖形來表征。它代表的含義是已知一個電子在參考點位置，由于Pauli互斥作用而在其它位置發現另一個同自旋電子的概率變化，其值處處為負，全空間積分值為-1。電子只能運動到以它的位置為參考點時Fermi穴函數有數值分布的區域（數值越大機會越大），如果范圍很小，則這個區域有較高的定域性。定域化函數與Fermi穴函數是有很大共通之處的，這些問題我會單獨撰文介紹，這里就不再多談了。

為了方便寫腳本，Multiwfn中Fermi穴函數的參考點可以在主功能1000里設定（也可以在settings.ini里設定）。為了繪制參考點移動過程的動畫，就需要生成一批重定向文件，每個文件里依次使用不同的參考點位置。此例的CH3Br中C和Br的x、y坐標皆為0，z坐標分別為-2.88和0.79 bohr，我們想讓參考點從z=-4.1移動到z=2.1，x、y都為0，過程共31幀。執行下面的shell腳本createinpctr.sh就可以在當前目錄生成31個重定向文件了。

#!/bin/bash
nstep=30
refx=0
refy=0
refz=0
refzmin=-4.1
refzmax=2.1
refz=$refzmin
rangez=`echo "$refzmax-($refzmin)"|bc`
dz=`echo "scale=5;$rangez/$nstep"|bc`
for ((i=1;i<=$(($nstep+1));i=i+1))
do
cat << EOF > `printf "%4.4i" $i`.txt
1000              #特殊功能
1                 #設定Fermi穴參考點
$refx,$refy,$refz #Fermi穴參考點數值
4
17
2
200,200
3
0                  #X=0的YZ平面，是CH3Br中的一個對稱面
-7                 #Multiwfn 2.02新增功能，對所有數據乘上某個值
-1                 #乘上的值為-1
0

EOF
refz=`echo "$refz+$dz"|bc`  #每個重定向文件內的參考點z值依次增加
done
由于shell腳本不能直接支持浮點運算，所以這個腳本利用了linux下bc命令來做運算。之所以所有數據乘上-1，是因為默認等值線設定下負值和正值分別使用虛線和實線表示，而Fermi穴函數都為負值，滿屏虛線看起來有點亂，乘上-1后就會通過實線來表示了，會好看些。將這個腳本生成的0001.txt至0031.txt文件拷到Multiwfn所在目錄下的inpctr子目錄下。

在Multiwfn目錄下寫一個批處理文件batchrunsinfil.bat，內容如下
for /f %%i in ('dir inpctr\*.txt /b') do Multiwfn CH3Br.wfn < inpctr\%%i
將CH3Br.wfn也放到Multiwfn所在目錄下。這個波函數是Hartree-Fock方法生成的，注意研究Fermi穴函數時切勿使用DFT波函數，因為DFT波函數本意用于產生基態密度和獲得較準確動能，而并不是用于逼近真實波函數，它的Fermi穴函數意義不清。

執行batchrunsinfil.bat，將得到的31幅圖像都放到/sob/CH3Br目錄下，然后進入/sob目錄，執行
convert -delay 20 -resize 600 -colors 128 -shave '50x0' -monitor -fill red -annotate +5+475 "Created by Sobereva" CH3Br/DISLIN.PNG CH3Br/DISLIN_?.PNG CH3Br/DISLI_??.PNG CH3Br/CH3Br_Fermi.gif
得到的CH3Br_Fermi.gif如下圖所示：

藍色米字符號表示的是參考點位置。有兩個氫離平面太遠所以沒有標出。從圖中可見，在一開始參考點與C有一段距離的時候Fermi穴函數分布很廣，也包圍了氫原子，這表明在參考點的電子有很大離域性，這個參考點附近電子定域性因此也不大。而到了碳原子核附近后，Fermi穴函數分布緊密集中在碳原子核附近，這說明這個位置的電子的運動被明顯束縛在有限空間了，實際上這正是內層s電子的情況。有文章指出具有近完美的定域性的區域只會出現在原子內核。當參考點移動到C-Br鍵中間時雖然電子也離域到C-Br的兩端，但中間有很廣且數值較高的區域，這是為什么共價鍵區域定域化函數比較大。當參考點與Br開始接近的時候Fermi穴函數分布范圍開始急劇縮小，和C的情況類似，當穿過Br之后離域性又開始變大。

5 結語

限于時間和精力，本文只給出了幾個基本化學過程的動畫的制作方法并對其包含的信息進行了討論，但我想仍足矣證明動畫這種形式的強大表現力，也足以闡明制作的基本流程。雖然制作過程看似步驟較多，實際上熟練之后會發現很容易，無非是生成波函數、用Multiwfn批量繪圖、用convert合并。制作不同化學過程的動畫需要對上述過程做不同細微的調整，需隨機應變。曲線的動畫和等值面圖的動畫在本文沒有涉及，前者對于研究一維上的問題很適合（比如前文的H2拉伸過程，用曲線圖在定量上更清楚），后者對研究在空間上比較復雜的過程是很有用的，制作并不困難，方法參見文章開頭提到的《通過鍵級曲線和ELF/LOL/RDG等值面動畫研究化學反應過程》。