使用Multiwfn計算軌道的體積
使用Multiwfn計算軌道的體積
Using Multiwfn to calculate volume of orbitals
文/Sobereva@北京科音 2025-Jan-1
1 前言
有人在思想家公社QQ群（http://www.shanxitv.org/QQrule.html）里問怎么計算分子軌道的體積，并且貼了一張軌道等值面圖。只要計算出這個等值面里包圍的體積就可以達到這個目的（雖然也可以有其它方式的定義，但沒這么簡單直觀）。這里我介紹一下怎么利用Multiwfn程序的定量分子表面分析功能來實現這個目的。定量分子表面分析功能之前在《使用Multiwfn的定量分子表面分析功能預測反應位點、分析分子間相互作用》（http://www.shanxitv.org/159）有簡要介紹，結合Multiwfn的極度靈活的設計，這個功能還能做很多其它的分析，正如本文所展示的。
如果讀者不熟悉Multiwfn，應當閱讀《Multiwfn FAQ》（http://www.shanxitv.org/452）、《Multiwfn入門tips》（http://www.shanxitv.org/167）、《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。Multiwfn可以在http://www.shanxitv.org/multiwfn免費下載，不要用上古版本。
2 實例
這里以Multiwfn計算下面這個D-pi-A體系的HOMO軌道的體積為例，圖中的軌道等值面數值為0.05，此體系的波函數文件是Multiwfn目錄下的examples\excit\D-pi-A.fchk。注意等值面包圍的體積直接依賴于等值面數值的選取，橫向對比時必須統一。
啟動Multiwfn，然后輸入
examples\excit\D-pi-A.fchk
5 //計算格點數據
4 //軌道波函數
h //HOMO
4 //設置格點間距
0.25 //0.25 Bohr是精度和耗時的較好權衡。體系大的話也可以用大一些的格點間距來節約耗時但會損失一些精度
0 //返回主菜單
13 //處理格點數據
11 //格點數據運算
13 //取絕對值。軌道波函數有正有負，這一步使之都變成正值
-1 //返回主菜單
12 //定量分子表面分析
1 //選擇定義表面的方式
11 //用內存中的格點數據定義等值面
0.05 //等值面數值
6 //開始表面分析且不考慮被映射的函數
從屏幕上可以看到以下結果，即體積為14.2 Angstrom^3，順帶著軌道等值面的面積69.0 Angstrom^2也顯示出來了。
Volume: 95.86111 Bohr^3 ( 14.20515 Angstrom^3)
Estimated density according to mass and volume (M/V): 25.0417 g/cm^3
Overall surface area: 246.50998 Bohr^2 ( 69.02983 Angstrom^2)
感興趣的話，還可以選擇選項-3 Visualize the surface看一下當前等值面對應的圖像，體積就是這個等值面里面包圍的部分
對比一下，看看下面這個第18號分子軌道的體積如何
退回到主菜單，把下面這一串命令復制到Multiwfn窗口里即可，結果是7.3 Angstrom^3，可見只有HOMO的一半，和肉眼看上去的軌道體積差異相一致。
5
4
18
4
0.25
0
13
11
13
-1
12
1
11
0.05
6
3 批量計算
在《詳談Multiwfn的命令行方式運行和批量運行的方法》（http://www.shanxitv.org/612）里專門講了怎么用Multiwfn批處理計算。這里提供一個腳本，可以便捷地用前文的方法批量計算特定序號范圍的軌道的體積。http://www.shanxitv.org/attach/734/orbvol.sh是個Bash shell腳本，用來計算examples\excit\D-pi-A.fchk這個體系所有占據軌道（1到56號）的體積，內容如下所示。其中istart和iend是起始和終止的軌道序號，這倆值和被處理的波函數文件路徑都可以根據實際需要修改。運行之前，Multiwfn可執行文件所在目錄需要添加到PATH環境變量下使得能夠通過Multiwfn命令啟動之。
#!/bin/bash
istart=1
iend=56
for ((i=$istart;i<=$iend;i=i+1))
do
cat << EOF >> orbvol.txt
5
4
$i
4
0.25
0
13
11
13
-1
12
1
11
0.05
6
-1
-1
EOF
done
echo "q" >> orbvol.txt
echo "Running Multiwfn..."
Multiwfn examples/excit/D-pi-A.fchk < orbvol.txt >> result.txt
grep "Volume:" result.txt | nl -v$istart
rm ./orbvol.txt result.txt
運行后，屏幕上看到各個軌道的體積：
1 Volume: 2.88332 Bohr^3 ( 0.42726 Angstrom^3)
2 Volume: 2.89404 Bohr^3 ( 0.42885 Angstrom^3)
3 Volume: 2.42980 Bohr^3 ( 0.36006 Angstrom^3)
4 Volume: 2.40622 Bohr^3 ( 0.35656 Angstrom^3)
...略
54 Volume: 95.10824 Bohr^3 ( 14.09359 Angstrom^3)
55 Volume: 95.86240 Bohr^3 ( 14.20534 Angstrom^3)
56 Volume: 95.86111 Bohr^3 ( 14.20515 Angstrom^3)
將軌道體積相對于軌道序號作柱形圖，如下所示。可見前16號軌道的體積非常小，這是因為它們都是內核軌道。
筆者之前寫過《通過軌道離域指數(ODI)衡量軌道的空間離域程度》（http://www.shanxitv.org/525），里面用我提出的ODI指數考察了與本文相同的D-pi-A體系的各個軌道的離域程度，并且繪制了ODI與軌道序號之間的柱形圖。將上圖與ODI的柱形圖相對比，可以看到兩個圖有明顯的互補性，即ODI越低（體現軌道越離域、越能同時覆蓋很多原子），軌道體積傾向于越大。如果你的目的就是想衡量軌道的離域程度，ODI相對更嚴格，畢竟不依賴于軌道等值面數值的選取（軌道等值面數值設得很大的話，軌道等值面就完全看不到了，體積為0）。而用軌道體積來說事的好處是可以直接與等值面圖相對應，結合圖像討論非常清楚直觀。
4 其它
本文的做法不限于計算考察分子軌道的體積，對任意軌道都可以考察，如定域化軌道、NTO軌道、雙正交化軌道、NAdO軌道、AdNDP軌道等等，載入含有相應軌道的波函數文件即可，產生方法在下文說了。
Multiwfn的軌道定域化功能的使用以及與NBO、AdNDP分析的對比
http://www.shanxitv.org/380
用于非限制性開殼層波函數的雙正交化方法的原理與應用
http://www.shanxitv.org/448
使用鍵級密度(BOD)和自然適應性軌道(NAdO)圖形化研究化學鍵
http://www.shanxitv.org/535
使用AdNDP方法以及ELF/LOL、多中心鍵級研究多中心鍵
http://www.shanxitv.org/138
效仿本文的做法還可以計算其它任意實空間函數的等值面內包圍的體積。比如計算自旋密度的等值面內的體積，只需要在主功能5里選擇被計算的函數的那一步選擇自旋密度即可。如果是像ELF那樣處處為正的函數，計算其等值面內的體積前無需像前文那樣需要先進入主功能13將其取絕對值；反之如果函數處處為負，則也需要先取絕對值。