文/Sobereva@北京科音  2024-Feb-29

0 前言

筆者在《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）中介紹過，磁感生電流是一種重要、直觀的考察分子芳香性的方式。SYSMOIC是一個基于continuous transformation of the origin of the current density (CTOCD)方法考察化學體系的磁感生電流相關問題的程序。本文將結合苯分子為例，介紹這個程序的最關鍵用法。SYSMOIC程序和相關原理的詳細介紹可以看原文J. Chem. Inf. Model., 61, 270 (2021)。

AICD和GIMIC也是兩個重要的考察磁感生電流的程序，筆者在下文有專門的介紹。這些文章里介紹的背景知識對理解本文很有幫助
? 使用AICD 2.0繪制磁感應電流圖（http://www.shanxitv.org/294）
? 使用AICD程序研究電子離域性和磁感應電流密度（http://www.shanxitv.org/147）
? 考察分子磁感生電流的程序GIMIC 2.0的使用（http://www.shanxitv.org/491）
? Utilizing AICD and GIMIC programs to study magnetically induced current density（http://www.shanxitv.org/148）
之所以筆者寫此文介紹SYSMOIC，在于它對以上兩個程序有明顯互補性。AICD程序是將感生電流以箭頭方式繪制在AICD的等值面上，而等值面以外區域的感生電流就無法在圖上體現出來了，而且AICD沒法給出感生電流的鍵截面積分值（鍵電流強度）來定量考察電流大小。SYSMOIC則能夠給出三維空間或二維平面上的感生電流向量場圖，從而能展現出AICD圖表現不出來的信息，而且SYSMOIC還能非常方便地給出鍵電流強度。GIMIC程序雖然能做感生電流積分、結合Paraview程序能給出漂亮的向量場圖，但是GIMIC卻無法將感生電流分解成不同分子軌道的貢獻，而SYSMOIC則能夠方便地指定計算感生電流時只考慮或不考慮哪些分子軌道。可見SYSMOIC有獨特的價值。

本文介紹的情況對應SYSMOIC 2024.1版，對其它版本不一定適用，請結合相應版本程序提示和手冊隨機應變。

1 SYSMOIC的基本特征

SYSMOIC是免費而不開源的程序，Windows、Linux、Mac版都有，可以在http://SYSMOIC.chem.unisa.it/DISTRIB下載。諸如STABIN_LNX-2024.1.tar.gz就是Linux的2024.1版可執行文件包。Linux版需要用較新的系統，筆者用Rocky Linux 9可以正常使用，CentOS 7.4用不了。

SYSMOIC是由一堆子程序組成的，對應壓縮包里不同可執行文件，Windows版有.exe后綴，其它平臺的沒有。例如JBMAP子程序用于計算感生電流格點數據，BOCUST子程序用于計算鍵電流強度，等等。v3d子程序用來對其它子程序計算出來的.3d文件進行可視化，如繪制感生電流矢量圖。各個程序的用法在SYSMOIC的在線手冊http://sysmoic.chem.unisa.it/MANUAL/里都有說明。

SYSMOIC是解壓即用的程序。解壓后只需要把SYSMOIC目錄加入到操作系統的PATH環境變量里（不會的話自行google），然后就可以在任何目錄下輸入子程序名直接使用了。

SYSMOIC有不少官方例子文件，是http://sysmoic.chem.unisa.it/WORKED_EXAMPLES/里的WORKEXA.tar.gz。例子按照體系/理論方法/基組對目錄層級進行命名，里面RUN_TESTS文件是對例子體系做各種分析的shell腳本，可以從中看到都是以什么命令調用程序的，以及傳入的命令文件是什么。REFERENCES目錄下的文件是運行RUN_TESTS腳本后產生的.3d文件和程序輸出信息，對前者可以直接用v3d作圖。

SYSMOIC必須結合最主流的量子化學程序Gaussian使用，目前只支持閉殼層形式的HF和DFT。只需要用Gaussian做NMR計算的同時產生wfx文件，再用SYSMOIC自帶的unpackwfx將之轉化成fort.3、fort.11、fort.23、fort.28這四個二進制文件，之后就可以用JBMAP、BOCUST等子程序做各種計算了，算完了可以用v3d子程序直接作圖。

2 產生苯的SYSMOIC輸入文件

下面以最典型的芳香性分子苯為例，講解SYSMOIC繪制感生電流矢量圖和計算鍵電流強度的方法。其它的用法我覺得在一般研究中用到的機會很少，所以就不在本文介紹了，感興趣者請看程序手冊、原文并結合程序自帶的例子文件摸索。下文涉及的各種輸入輸出文件都可以在http://www.shanxitv.org/attach/702/file.rar中獲得。

首先寫一個在B3LYP/6-31G*級別下做NMR計算的Gaussian輸入文件benzene.gjf，內容如下，在本文的文件包里已經提供了。其中的坐標是在當前級別下已經優化過的。NMR=CSGT要求以CSGT方法做NMR計算，out(wfx,CSGT)要求產生.wfx波函數文件，并且把CSGT計算過程中產生的擾動的軌道展開系數寫入到wfx文件中，輸出位置為末尾指定的D:\benzene.wfx。

# b3lyp/6-31g(d) NMR=CSGT out(wfx,CSGT)
[空行]
benzene
[空行]
0 1
 C                  0.00000000    1.39621600    0.00000000
 C                  1.20915900    0.69810800    0.00000000
 C                  1.20915900   -0.69810800    0.00000000
 C                  0.00000000   -1.39621600    0.00000000
 C                 -1.20915900   -0.69810800    0.00000000
 C                 -1.20915900    0.69810800    0.00000000
 H                  0.00000000    2.48297800    0.00000000
 H                  2.15032200    1.24148900    0.00000000
 H                  2.15032200   -1.24148900    0.00000000
 H                  0.00000000   -2.48297800    0.00000000
 H                 -2.15032200   -1.24148900    0.00000000
 H                 -2.15032200    1.24148900    0.00000000
[空行]
D:\benzene.wfx
[空行]
[空行]

算完后就有了benzene.wfx，確保它在當前目錄下，并且已經把SYSMOIC目錄加入到了PATH環境變量中，然后運行unpackwfx benzene命令，當前目錄下立刻就出現了fort.3、fort.11、fort.23、fort.28，它們都已提供在了本文的文件包里。后文說的所有操作都是對于這四個fort文件都處于當前目錄下的情況而言的。

注意當前的苯分子是處在Z=0的XY平面上，用Multiwfn載入Gaussian的out文件然后進主功能0，或者用GaussView載入out文件并顯示笛卡爾軸，都可以確認這一點。后面我們要考察的都是外磁場方向垂直于苯環所導致的感生電流的情況。

3 繪制磁感生電流向量場圖的例子

3.1 默認方式繪制

先來看怎么產生一個最簡單的磁感生電流量場圖，之后再解釋細節。直接在命令行窗口輸入JBMAP運行之，屏幕上會出現好多默認用的參數，比如下面這行，就是告訴你默認的外磁場是(0 0 1)矢量，顯然正垂直于當前苯環，和我們期望的一致
magnetic field...........B    0.0000    0.0000    1.0000

程序問你當前顯示出來的參數是否ok，這里就直接按回車代表用默認的y（yes）。然后再依次問你三個問題：
? external grid points y/n? [n]：即是否讀取外部格點，這里直接按回車代表用默認的n（如果要計算的點不是均勻分布在一個矩形區域，而是比如分布在一個球層表面，就需要輸入y并且提供格點定義文件）
? reference arrow y/n? [n]：含義不明。直接按回車用默認的n
? B direction y/n? [n]：即是否把外磁場箭頭畫在圖上，直接按回車用默認的n。如果你想把外磁場箭頭繪制在圖上就輸入y，并且輸入箭頭的位置和長度

馬上當前目錄下就有了JBM.3d文件。輸入v3d JBM.3d通過v3d程序對這個文件作圖，就蹦出了圖形窗口，如下所示，可以看到在一個矩形區域內均勻分布的每個格點位置都有一個箭頭指示此位置的磁感生電流的方向和大小。

想退出的話在圖形窗口按ESC鍵，或者在文本窗口按ctrl+C。

上面的圖非常丑陋，箭頭特別擁擠，分子結構也看不清楚，完全沒法用。下一節講怎么對JBMAP的計算進行自定義。

3.2 自定義JBMAP的運行參數

JBMAP有兩個方式指定運行參數，比如要把外磁場矢量改為0 1 1（程序自動會做歸一化成為0.0000 0.7071 0.7071），既可以運行JBMAP -B 0 1 1命令，也可以先啟動JBMAP然后直接輸入B 0 1 1后按回車，在屏幕上顯示的參數中都會看到外磁場矢量已經變更成自設的值了。

JBMAP命令后面可以直接接的選項可以輸入JBMAP -h查看，其中幾個關鍵的在這里提一下：
-o：輸出的.3d文件名，默認是-o JBM
-f：原子球大小的倍數，默認為-f 0.2。想讓作出的圖里原子球大一些就把這個改大
-B：外磁場矢量，默認為-B 0 0 1
-j：感生電流類型，默認是-j TOT，即diatropic和paratropic的總和。寫-j PARA代表只考慮paratropic部分，寫-j DIA代表只考慮diatropic部分
-q：考慮的軌道。例如只想考慮17 20 21三個分子軌道對磁感生電流的貢獻，就寫-q +3 17 20 21，只去掉15 16號分子軌道的貢獻就寫-q -2 15 16。默認是考慮所有分子軌道的貢獻（只有占據軌道才有貢獻）
-m：設置SYSMOIC計算感生電流用的CTOCD方法的具體形式，參看http://sysmoic.chem.unisa.it/MANUAL/S4.SS12.html的說明。當基組是完備的時候這些方法結果等價。有些選項是僅對于某些方法才需要設，比如用common gauge-origin的時候需要用-g自定義原點坐標（默認為0 0 0）。默認用的CTOCD-DZ2基本上對于所有情況都是合適的選擇，也不用設額外參數。關于CTOCD不同形式的介紹和差異，可以看前述的SYSMOIC原文以及感生電流綜述文章WIREs Comput Mol Sci, 6, 639 (2016)的Continuous Transformations of the Origin of the Current Density一節。

JBMAP運行后可以修改的參數當中較重要的如下：
? B：外磁場矢量，默認B 0 0 1
? JTERM：感生電流類型，默認是JTERM 0，即diatropic和paratropic的總和。1和2分別代表paratropic和diatropic
? RI和RF：分別指定格點所均勻分布的矩形盒子的X、Y、Z最小笛卡爾坐標和最大笛卡爾坐標。默認是根據體系坐標自動確定以覆蓋整個體系。對上面的苯分子的例子，從屏幕上默認的參數可見，當前對應的是RI -6.0635 -6.6921 -2.0000和RF 6.0635 6.6921 2.0000
? FATT：給感生電流乘的系數。比如大部分箭頭在圖上都顯得很不顯著，就可以設大這個值以更便于觀看
? FMM：磁感生電流過濾范圍。例如FMM 0.01 1就代表忽略感生電流的模小于0.01和大于1的格點（是對乘上FATT之前而言的），免得干擾觀看。默認0 0.1
? STEP：格點間距，數值越小箭頭越密、計算越耗時。默認0.4 Bohr
? LUNA：0和1分別代表箭頭的長度和面積正比于感生電流的模，默認是后者
? METH：計算方法，默認的2對應CTOCD-DZ2

為了方便，我推薦以命令行方式運行JBMAP。交互式程序怎么以命令行方式運行，我在《詳談Multiwfn的命令行方式運行和批量運行的方法》（http://www.shanxitv.org/612）里有很詳細說明。下面以這種方式結合自定義參數用JBMAP再次對苯進行計算。創建一個文本文件叫CDMAP.txt，內容如下，每一行對應一條傳給JBMAP程序的命令。此例RI和RF的Z值都是1，X和Y范圍都是-6到6 Bohr，因此格點會分布在苯分子上方1 Bohr的XY平面上的正方形區域內。此文件后面的y、n、n、n對應交互式界面里程序向你提問時依次輸入的四個指令。

B   0 0 1
RI -6 -6 1
RF  6  6 1
FMM 0.01 1
FATT 3
STEP 0.4
LUNA 1
JTERM 0
y
n
n
n

運行JBMAP -f 0.5 < CDMAP.txt，然后再運行v3d JBM.3d，看到下圖，可見清楚地把苯分子上方的磁感生電流描繪了出來。

再來看只考慮pi軌道時的磁感生電流圖的繪制。把benzene.wfx文件載入Multiwfn，進入主功能0，按照《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）說的，從HOMO開始往下挨個看軌道，就能找到三個pi軌道，序號為17、20、21。原子數多的時候也可以用《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）介紹的做法讓Multiwfn自動找出來所有pi軌道序號。運行JBMAP -q +3 17 20 21 -f 0.5 < CDMAP.txt，然后再用v3d繪圖，得到下圖，這便是苯的pi電子的感生電流圖了。相對于之前的所有電子貢獻的圖，此圖少了sigma電子造成的出現在內環的逆時針（paratropic電流）特征，而且也沒有了氫原子區域的感生電流，因為pi電子不在氫上。

3.3 v3d的使用細節

下面專門細談一下v3d界面的操作。v3d圖形窗口出現后，可以用鼠標拖動體系旋轉視角。還可以在圖形窗口處于激活狀態下按h鍵，此時文本窗口會出現可以用的所有按鍵的功能說明。其中較重要的按鍵在這里說一下。

先說視角調節按鍵：
a、b、g：視角由alpha、beta、gamma值定義，按a、b、g鍵會分別增加它們1.0。如果按住shift鍵再按，也即按A、B、G，會分別增加它們10
? f和F：調節縮放因子。按一下f圖像就會放大一點，按一下F就會放大很多
? d和D：調節視角距離。按一下d就會增加0.1，按一下D就會增加1。視角距離越遠，近大遠小的透視效果越弱、越接近正交視角的效果（據我所知v3d沒法完全用正交視角）
? n和N：設置霧化距離閾值。超過視角一定距離閾值開始就會用霧化效果使得遠處物體被屏蔽。按n增加閾值0.1，按N增加1.0
? C：令圖像居中

每次按以上按鍵調節設置的時候在文本窗口都會看到當前值。如果你想把這些鍵的操作效果反過來，比如按一次f就令圖像縮小一點，那就先按一下c鍵進入反轉操作狀態（文本窗口里會看到c變成了c-1），然后再按f就行了。再次按c則還原到普通操作狀態。

按p鍵可以把當前圖像保存成當前目錄下的.ps文件。這是矢量圖像格式，在Windows下可以用裝了ghostscript插件的irfanview或photoshop等程序打開，也可以用Acrobat distiller轉成pdf再打開。在Rocky Linux、CentOS等系統里也自帶了ps文件觀看程序。

默認情況下產生的ps文件畫出來的感生電流的箭頭會有較厚的白邊，當箭頭多的時候白邊會導致后面的箭頭被遮擋，比較難看。解決辦法是按一下c鍵先進入反轉操作狀態，然后按住k鍵一會兒使得白邊量減小（文本窗口不會提示減小到了多少），然后再導出ps文件，然后就會發現白邊非常細了。

如果想保存當前的作圖設置，按S，當前目錄下就出現了SpecialValues.txt文件，里面記錄了作圖設置，可以看到f、d等數值都記錄在內。當再次啟動v3d時，程序若發現當前目錄下有這個文件，就會自動載入里面的設置。如果你希望每次產生的ps文件里箭頭的白邊都非常細，可以啟動v3d后先直接按S保存一個記錄默認狀態的SpecialValues.txt文件，然后把里面對應白邊粗細的ridfr前面的值改為0.1。這樣每次啟動v3d就會默認用這個設置了。

4 計算鍵電流強度

對鍵截面做感生電流積分的概念我在《考察分子磁感生電流的程序GIMIC 2.0的使用》（http://www.shanxitv.org/491）中已經專門講過了，這可以定量考察穿越特定化學鍵的電流強度從而便于對比，這在SYSMOIC里叫鍵電流強度。SYSMOIC算鍵電流強度遠比用GIMIC更方便，都省得自己定義積分范圍，結果還能直接可視化。下面我們對苯的一個C-C鍵計算鍵電流強度。

直接在命令行窗口中輸入BOCUST命令運行此程序，程序首先會提示
| plotted    12 atoms
| plotted    12 bonds
 ====> Pair C1    -  C2     distance=    1.396 Ang
這說明程序發現有12個原子對兒的距離都處在默認的距離下限（0.9埃）和距離上限（1.6埃）之內。其中C1-C2是第一個原子對，現在問你是否對這個C1-C2鍵計算鍵電流強度。這里輸入y然后按回車，很快就會出現結果
CONTRIBUTION  1 -0.6000526E+00 AU -0.1690912E-07 SI
CONTRIBUTION  3  0.1784109E+00 AU  0.5027513E-08 SI
TOTAL CS ==>    -0.4216417E+00 AU -0.1188161E-07 SI
(2) CS/CS0 RELBEN FICS     0.99    3.00  0
上面顯示的兩個CONTRIBUTION分別是diatropic和paratropic對鍵電流強度的貢獻，由于方向相反，顯然符號相反。二者之和取絕對值0.4216417 a.u.就是鍵電流強度大小，換算成常用的nA/T單位就乘以28.179409，即11.88 nA/T。這個值與《考察分子磁感生電流的程序GIMIC 2.0的使用》（http://www.shanxitv.org/491）在B3LYP/6-311G*級別通過GIMIC程序得到的12.086 nA/T非常接近。上面信息中的CS/CS0 RELBEN FICS     0.99代表當前值與參考值的比值是0.99。這里參考值是程序內置的對苯在某個級別計算的C-C鍵電流強度12 nA/T。

接下來程序又問你是否對C1-C6計算鍵電流強度。由于當前所有C-C鍵都是等價的，只需要計算一個就行了，因此輸入N回車，程序就結束運行了。如果你想把所有鍵都進行計算，就接著輸入11次y回車把剩下11個鍵都算完。當前體系一共6個C-C鍵、6個C-H鍵。

當前目錄下出現了BCS.3d文件。運行v3d BCS.3d，就可以看到下圖

此圖里黑色箭頭體現了C1-C2鍵電流的大小和方向。積分截面處在鍵的正中央并與鍵垂直，圖中的曲線是這個截面上的感生電流0.001 a.u.等值線，曲線內部區域是對感生電流積分的區域，紅色和藍色分別對應paratropic和diatropic部分。可見diatropic分布在環的外側，這與之前從向量場圖中看到的順時針環電流的主要分布是一致的。圖中的紅色小點對應截面上感生電流的極值點。鍵上面標的文字99對應于前述的參考比值0.99乘以100。文字出現的位置在鍵的正中央，不容易分辨，可以在圖形窗口激活狀態下按W鍵關閉文字顯示，然后再自行ps上。

如果你想一次性把所有鍵的電流強度都計算，又不想一次次麻煩地輸入y，那么可以寫一個文本文件比如叫ally.txt，里面有12行，每一行都是y。運行BOCUST < ally.txt，這樣所有鍵就都會被計算了。此時對得到的BCS.3d繪制會看到下圖

下面介紹一下BOCUST運行時可以接的選項中比較重要的：
-o、-f、-B、-m、-q等：和前面介紹的JBMAP程序的情況完全相同。即BOCUST也可以用-B指定外磁場方向、用-q指定考慮哪些軌道，等等。
-d、-D：確定哪些原子對兒被考慮時用的距離下限和上限（埃），默認-d 0.9 -D 1.6
-e：確定鍵截面上做積分的區域用的感生電流等值線的數值（a.u.），默認-e 0.001
-CS0：參考電流強度（nA/T），默認-CS0 12
-nocd：不顯示積分區域和極值點。如果想讓圖上只有箭頭而沒有曲線和小點，應該加上這個
-l [比值]：在截面上繪制感生電流的不同數值的等值線，后面跟的是相鄰等值線數值的比值
-C [粗細]：繪制截面的邊框

例如這里計算C1-C2的sigma電子的鍵電流強度（即扣除三個pi軌道貢獻），并在截面上顯示等值線、顯示邊框，就運行BOCUST -q -3 17 20 21 -l 3 -C 5，然后輸入y回車，再N回車。之后用v3d對BCS.3d作圖看到下圖。可見整個截面，以及里面的感生電流等值線，都顯示得很清楚。由于這個截面上diatropic和paratropic的大小都差不多（畢竟苯分子沒有sigma芳香性特征），所以鍵電流強度僅為0.007 a.u.。

5 總結

此文對分析磁感生電流的程序SYSMOIC做了介紹并給出了具體例子。由例子可見，通過SYSMOIC結合Gaussian繪制感生電流向量場圖很方便，不需要借助第三方程序。SYSMOIC做鍵電流強度計算也很自動化，不需要像GIMIC那樣謹慎地定義積分區域，而且鍵電流大小和方向還能直接通過箭頭直觀顯示出來。SYSMOIC的計算速度也比較快。

本文只涉及了SYSMOIC最常用的功能，還另外有一些其它功能，比如對感生電流密度做拓撲分析，計算各種感生電流密度相關的場（ACID及其各向異性、磁屏蔽密度、vorticity張量的跡及其各向異性等）并繪制成等值面和填色平面圖、計算任意格點上的感生電流張量和導數等等。

文/Sobereva@北京科音

First release: 2024-Feb-16   Last update: 2024-Feb-17

0 前言

Hirshfeld surface分析（以下簡稱HS分析）是展現分子晶體中一個或多個分子與周圍分子間相互作用的方法，它也同樣可以用于展現孤立體系中特定片段與周圍原子間的相互作用。HS分析的原理容易理解，圖像比較直觀，在分子晶體的研究領域已經用得非常普遍。Multiwfn的主功能12中很早以前就已經實現了HS分析，如今已經被不少文章使用。本文將結合許多例子，專門詳細具體講解一下Multiwfn做HS分析的各方面細節、操作和技巧。

本文的讀者請務必使用2024-Feb-16及以后更新的Multiwfn版本，否則與本文所述情況會有很多不同。Multiwfn可以在其主頁http://www.shanxitv.org/multiwfn免費下載。不了解Multiwfn者請參看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。使用Multiwfn做HS分析在寫文章時請按照Multiwfn啟動時的提示對程序進行恰當的引用。

很值得一提的是，筆者提出的IGMH方法也非常適合展現分子晶體中的分子間相互作用，HS分析與IGMH分析的展現形式有明顯區別且有一定程度的互補性，二者亦可以同時使用以提供更多視角。請閱讀《使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用》（http://www.shanxitv.org/621）和《一篇最全面介紹各種弱相互作用可視化分析方法的文章已發表！》（http://www.shanxitv.org/667）提到的綜述了解IGMH的相關知識，里面也專門有IGMH用于分子晶體的實際例子。

1 Hirshfeld surface分析的基本思想

在筆者講授的量子化學波函數分析與Multiwfn程序培訓班（http://www.keinsci.com/workshop/WFN_content.html）和Multiwfn手冊3.15.5節對HS分析的原理有很具體、詳細的講解，下文只是把原理的最關鍵的部分簡要介紹一下，對于正確做HS分析基本夠了。

HS分析最早由Spackman等人于Chem. Phys. Lett., 267, 215 (1997)提出，后來又得到了發展，Acta Cryst., B60, 627 (2004)和CrystEngComm, 11, 19 (2009)是其兩篇綜述文章。HS分析關鍵思想是對分子晶體中的特定分子構造出Hirshfeld surface，它相當于這個分子在分子環境中的表面，然后再將一些有特殊意義的實空間函數映射到這個表面上，由此可以對分子晶體中的分子間相互作用特征進行考察。

簡要說一下Hirshfeld surface是怎么定義的。Hirshfeld在Theoret. Chim. Acta (Berl.), 44, 129 (1977)最早提出了一種定義化學體系中原子空間的方式，它給每個原子定義了Hirshfeld權重函數來描述這個原子在三維空間中各個位置所占權重，數值從0到1平滑變化，0和1分別對應于此位置完全不屬于和完全屬于這個原子。每個位置所有原子的權重函數加和為1。這是一種典型的模糊式原子空間定義方式。具體定義細節見《原子電荷計算方法的對比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）的2.5節和Multiwfn手冊3.9.1節，看完了就會知道產生Hirshfeld權重需要有各個原子的坐標，以及體系中的各種元素的原子在孤立狀態下的球對稱化的電子密度。將一個分子中所有原子的Hirshfeld權重函數加和就定義了這個分子的權重，因此分子晶體中各個分子都有各自的權重函數。某個分子的Hirshfeld surface就對應于它的權重函數數值為0.5的等值面。例如下圖的曲線展現了某個平面上各個分子的Hirshfeld surface對應的輪廓。可見，Hirshfeld surface算是分子環境中各個分子與其它分子接觸面的一種定義方式。下文所說的“表面”都是指Hirshfeld surface。

現實當中構造Hirshfeld surface是利用我在J. Mol. Graph. Model., 38, 314 (2012)介紹的Marching Tetrahedron或類似的算法實現的。這個表面被描述為大量小三角形的集合，每個三角形由三個表面頂點構成。由于構造Hirshfeld surface對應的等值面用的算法和《使用Multiwfn的定量分子表面分析功能預測反應位點、分析分子間相互作用》（http://www.shanxitv.org/159）涉及的構造電子密度等值面對應的分子范德華表面的算法在本質上相同，所以Multiwfn中HS分析也是在主功能12（定量分子表面分析）中實現的。

實際中做HS分析的時候是選取一個（也允許是多個）感興趣的分子并對它構造Hirshfeld surface進行分析，這個分子在下文管它叫“中心分子”，它周圍的分子也可以稱為“環境分子”。

HS分析定義了一些映射到Hirshfeld surface上的三維實空間函數，常見的有：
? d_i：表面內部（interior）的原子（也即中心分子的原子）到當前點的最近距離
? d_e：表面外部（exterior）的原子（也即周圍分子的原子）到當前點的最近距離
? d_norm：歸一化的（normalized）分子間接觸距離。由d_i、d_e和原子范德華半徑計算出來。數值越小體現表面上此處附近的內、外原子有越近的接觸，暗示此處的相互作用越強
? shape index（形狀指數）：數值在[-1,1]范圍，越負說明此處表面越凹，反之越凸
? curvedness（曲度）：在[-4,0.4]范圍，-4對應表面此處完全平坦，越正越凸，0對應單位球面的曲度

d_norm經常用來對Hirshfeld surface進行著色來直觀展現分子環境中的分子間的相互作用，d_norm數值較小的區域對應于較強烈的分子間相互作用。筆者發現用電子密度來著色也很有意義，Hirshfeld surface上電子密度越大的地方體現相互作用越強，效果比用d_norm著色時明顯更好，色彩過渡更為平滑，物理意義也更強，也沒有范德華半徑選取的任意性。由于對大體系做量子化學計算產生比較精確的電子密度比較耗時，因此只需要用準分子近似的電子密度（promolecular density）對Hirshfeld surface著色就夠了，它直接由各個原子孤立狀態的電子密度疊加得到，耗時極低，我實測和使用量子化學計算的電子密度著色的效果差不多。

由于各種元素的孤立狀態的電子密度以及范德華半徑在Multiwfn中是內置的，故產生Hirshfeld surface以及計算以上提及的各種函數只需要原子坐標和元素信息就夠了。由于HS分析依賴的信息非常簡單，不牽扯基于波函數的計算，因此耗時極低，用起來也很方便。

HS分析還經常繪制指紋圖（fingerprint map），是把d_i和d_e作為散點圖的橫軸和縱軸，然后把Hirshfeld surface上的各個頂點根據其d_i和d_e的數值繪制在指紋圖上作為一個個小點。根據指紋圖上散點分布的位置可以對中心和周圍分子間的相互作用特征進行討論，后文會結合具體例子來講。

HS分析中還經常做局部接觸（local contact）分析。完整的HS分析描述了中心分子的所有原子和周圍分子的所有原子間的相互作用，而局部接觸分析可以指定在HS分析中只考慮中心分子的哪些原子和周圍分子的哪些原子的相互作用。比如可以了解在整個Hirshfeld surface中體現中心分子的氧和周圍分子的H之間的相互作用的區域的位置和面積。

Multiwfn還支持Becke surface分析，是我自己提出的概念。它和HS分析的唯一差別是用Becke權重函數而非Hirshfeld權重函數的0.5等值面來定義表面。二者實際效果差別不大，一般沒必要用Becke surface，但它的一個特殊好處是允許表面出現在電子密度為0的區域，此處沒法定義Hirshfeld權重并構造Hirshfeld surface，詳見本文的第5節。Becke權重函數的具體定義方式參見《密度泛函計算中的格點積分方法》（http://www.shanxitv.org/69），它基于原子坐標和原子共價半徑得到。出現在相互作用的原子間的Becke surface會離半徑較小的原子較近、離半徑較大的原子較遠，Hirshfeld surface也有這樣的特點，這是由于定義它的準分子密度分布特征所自然而然帶來的。

2 Multiwfn的HS分析的功能

Multiwfn中做HS分析需要提供含有原子信息的文件作為輸入文件，如.pdb、.xyz、.mwfn、.cif、.fch、.mol2、.gjf等等等等，詳見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。

HS分析大多研究的是分子晶體，cif是最常用的記錄晶體結構的格式。通常不能載入cif文件后上來就做HS分析，因為HS分析一般需要提供一個中心分子+環境分子的簇模型，這樣才能靠HS分析考察中心分子與環境分子的相互作用，而cif文件記錄的是晶胞里各個原子的坐標，分子往往是被截斷的，中心分子或環境分子一般都不完整。因此首先需要用《Multiwfn中非常實用的幾何操作和坐標變換功能介紹》（http://www.shanxitv.org/610）中介紹的自動挖團簇的功能構造簇模型。如果你研究的不是分子晶體的情況，就是比如分子二聚體中兩個分子間的相互作用，那么就不牽扯挖簇的過程了，直接提供含有二聚體坐標信息的文件當輸入文件就行了。

HS分析在Multiwfn主功能12實現。進入這個功能后，首先選1把定義表面的方式切換為Hirshfeld surface。此時被映射到表面的函數會自動改為電子密度（對于輸入文件沒有波函數信息的情況具體是指準分子電子密度），你也可以選2把被映射的函數改成其它的。之后選0就開始計算了，Hirshfeld surface會被構造出來，組成它的所有表面頂點上的被映射的函數值會被計算出來，并顯示出Hirshfeld surface的面積和包圍的體積。然后會看到后處理菜單，里面有豐富的選項，提示得都很清楚，Multiwfn手冊3.15.5節都有解釋。利用后處理菜單的選項，可以導出用于在VMD程序中繪制著色的Hirshfeld surface圖要用的.cub文件，還可以繪制指紋圖、做局部接觸分析，這些在后文的例子里都有體現。

后文的例子涉及到的各種文件可以在http://www.shanxitv.org/attach/701/file.zip中獲得。

3 用Multiwfn對NAOB晶體做Hirshfeld surface分析實例

這一節以(Z)-4-((2-nitrophenyl)amino)-4-oxobut-2-enoic acid (NAOB)晶體為例做HS分析。NAOB的分子結構如下，其晶體在DOI: 10.1007/s13738-023-02904-9里進行了研究，文章的補充材料里給了它的cif文件的信息，是本文文件包里的NAOB.cif。

3.0 準備工作：構造簇模型

NAOB.cif對應的晶體結構如下，可見連一個完整的分子都沒有，因此在進行HS分析之前，我們顯然得先構造出中心分子+周圍分子的團簇結構才行。

為了構造簇模型，啟動Multiwfn，然后載入NAOB.cif，之后輸入
300  //主功能300
7   //幾何相關操作
25  //構造“中心分子+臨近一圈分子”的團簇
1  //當前晶胞里1號原子所在的分子將被作為中心分子（由于NAOB晶體里所有分子都是等價的，所以當前隨便輸入一個原子序號即可）
[回車]  //代表若一個周圍分子與中心分子間最近原子對距離小于這倆原子的范德華半徑和的1.2倍，則這個周圍分子就被整個納入團簇

Multiwfn瞬間就構造出了團簇，從屏幕上的提示可看到這個簇有375個原子，并且屏幕上還巨貼心地把中心分子中的原子序號給了出來，此例為1-10,17-19,22-25,205-210,214,215。把這個序號記下來，之后HS分析時要用到。

現在可以選當前菜單中的選項0看一眼新構造出的團簇是什么樣，如下所示，可見非常理想，確實是中心分子被周圍一層分子所圍繞（為了中心分子看得清楚，在Multiwfn圖形界面的菜單欄里選了Other settings - Set atom highlighting，然后輸入了前述的中心分子里的原子序號，使中心分子用青色高亮了）。

點圖形界面右上角的RETURN按鈕關閉圖形窗口，然后選擇-2 Output system to .pdb file并輸入NAOB_cluster.pdb，以將當前簇結構導出為當前目錄下的這個文件。這個pdb文件在本文的文件包里也提供了。

注：如《實驗測定分子結構的方法以及將實驗結構用于量子化學計算需要注意的問題》（http://www.shanxitv.org/569）所強調的，由于X光衍射實驗一般難以確定氫原子的準確位置，因此原理上做HS分析之前最好先固定重原子而優化一下所有氫原子的位置。用免費高效的CP2K程序對晶體結構優化氫原子位置然后再用Multiwfn摳團簇，或是先摳團簇再用Gaussian等量子化學程序優化氫，都是可以的。

3.1 繪制電子密度著色的Hirshfeld surface圖

這一節演示繪制基于準分子近似的電子密度著色的Hirshfeld surface圖，這是HS分析最重要的一種圖。除了Multiwfn外還會用到非常流行的VMD可視化程序，可以在http://www.ks.uiuc.edu/Research/vmd/下載，使用筆者現在用的VMD 1.9.3版肯定沒問題，用其它版本不保證能按照本文的例子正常作圖。

啟動Multiwfn，載入上一節產生的NAOB_cluster.pdb，然后輸入
12  //定量分子表面分析
1  //選擇定義表面的方式
5  //Hirshfeld surface
1-10,17-19,22-25,205-210,214,215   //中心分子的原子序號
0  //開始分析

從屏幕上可以看到許多信息，以下兩條是值得注意的，第一個是Hirshfeld surface包圍的體積，相當于分子晶體中屬于這個中心分子的體積，第二個是Hirshfeld surface的面積
Volume:  1644.06256 Bohr^3  ( 243.62494 Angstrom^3)
Overall surface area:         859.01840 Bohr^2  ( 240.54965 Angstrom^2)

現在看到了后處理菜單。選擇-2在當前目錄下導出surf.cub，它是中心分子的Hirshfeld權重的格點數據。再選13，Multiwfn會計算被映射的函數的格點數據并導出為當前目錄下的mapfunc.cub。在這個界面里如果選擇-3也可以直接在Multiwfn里預覽Hirshfeld surface，但沒有著色效果。

將剛剛產生的surf.cub和mapfunc.cub以及Multiwfn自帶的examples\scripts\目錄下的hirsh_rho.vmd作圖腳本一起拷到VMD目錄下（即啟動VMD后在其文本窗口里運行pwd命令看到的目錄）。啟動VMD，在文本窗口里輸入source hirsh_rho.vmd來執行作圖腳本，然后就看到了下圖。

上圖顯示出了中心分子的Hirshfeld surface，并根據電子密度進行了著色。通過hirsh_rho.vmd腳本里的mol scaleminmax top 1 0.0 0.015這條命令可知，默認用的色彩刻度范圍的下限和上限分別為0和0.15。根據腳本里color scale method BWR這行命令可知，當前用的色彩變化是藍-白-紅，色彩變化示意圖在Graphics - Colors - Color Scale里可以看到。因此上圖中越紅的地方就是Hirshfeld surface上電子密度越大的地方，無疑對應于越強的相互作用。上圖中有一塊非常大的紅色，這對應于中心分子羧基的O-H作為氫鍵給體、周圍一個分子的O作為氫鍵受體形成的顯著的氫鍵作用區域。它的旁邊還有一小塊淡紅色區域，對應于中心分子的羧基氧作為氫鍵受體與周圍分子的C-H形成的弱氫鍵。上圖中還有一些發白的區域，在VMD里旋轉圖像仔細觀察的話可以看到對應的是比較遠距離的C-H...O-N很弱氫鍵，以及pi-pi堆積和普通色散吸引作用顯著的區域。這些特征區域都可以自行用powerpoint之類畫個箭頭標注在圖上便于讀者看清楚。上圖還有很多偏藍色的區域，這些地方電子密度接近0，因此不牽扯任何值得一提的分子間相互作用。特別藍的地方也往往對應于分子晶體中的孔洞區域，這種地方電子密度自然特別低，非常建議感興趣的讀者按照《使用Multiwfn計算晶體結構中自由區域的體積、圖形化展現自由區域》（http://www.shanxitv.org/617）介紹的方法作圖考察。

要注意中心分子與各個方向的周圍分子都有相互作用，光是靠一張圖的話很難展現完整，因此文章中可以多給幾張圖展現不同視角的Hirshfeld surface圖。

下面再說一下怎么改進作圖效果。Hirshfeld surface圖的效果受到多方面影響：
(1)光源。可以通過VMD的Display菜單里的Lighting選項設置打開哪些光源。如果選Mouse - Move light，然后在圖形窗口中拖動，還可以移動特點光源的位置。
(2)材質。hirsh_rho.vmd默認對等值面用Translucent材質，可以自行在Graphics - Materials界面里對這個材質的具體定義進行調節。
(3)Graphics - Representation界面里的作圖設置。里面可以創建更多Rep，每個Rep都可以獨立設置顏色和材質，并且通過選擇語句可以定義各個rep顯示的原子，不懂選擇語句怎么寫的話參考《VMD里原子選擇語句的語法和例子》（http://www.shanxitv.org/504）。特別值得一提的是，當前的圖中每個分子都有一個獨立的residue編號并被分子中所有原子所共享，因此可以利用residue選擇特定分子。若想查詢某個分子的residue號，就選Mouse - Query，然后點擊這個分子上任意一個原子的正中央，在VMD的文本窗口中就能看到它的residue號了。

為了讓上面例子的圖像效果更好，我在VMD的Graphics - Representation里做了些修改：把第1個rep設為了residue 0專門用于顯示中心分子，用CPK方式顯示并把Bond Radius設為了0.5。點擊Create Rep按鈕新增一個Rep，選擇語句用residue 7 11 5使得三個與中心分子作用顯著的分子顯示出來，Material設Edgy，Drawing Method用Licorice，把Bond Radius設成0.1。點擊用來顯示Isosurface的那個Rep，再選擇Trajectory標簽頁，在里面把色彩刻度上限設為0.012（注意文本框只顯示兩位小數，輸入0.012后按回車會如實設為0.012），使得著色的色彩顯得更鮮明。在Display菜單里把所有四個光源都打開。最后選File - Render，選擇Tachyon(Internal)進行渲染。之后在圖像上展現特征作用的地方通過powerpoint進行標注，并且把Graphics - Color - Color Scale里顯示的色彩刻度條挪到圖上并適當拉伸、標記上色彩刻度上下限。最終得到下圖，可見對相互作用展現得非常直觀清楚。

如果要繪制d_norm著色的Hirshfeld surface圖，與上面的例子只有兩個差別
(1)進入主功能12并選擇用Hirshfeld surface方式定義分子表面后，選擇2 Select mapped function，再選d_norm。然后再開始分析
(2)作圖使用Multiwfn目錄下examples\scripts\里的hirsh_dnorm.vmd代替hirsh_rho.vmd
由于d_norm著色的圖明顯不如電子密度著色的圖色彩變化平滑，在很多地方有顏色的突越不好看，故這里就不多說了。

Multiwfn還可以給出Hirshfeld surface上被映射的函數（當前為電子密度）的極大點的位置和數值，便于定量對比討論。在Multiwfn后處理菜單中選8 Export all surface vertices and surface extrema as vtx.pqr and extrema.pqr，此時當前目錄下就出現了vtx.pqr和extrema.pqr，分別記錄了所有表面頂點和表面極值點的坐標和被映射的函數數值。如屏幕上的提示所示，extrema.pqr里碳和氧原子分別用于記錄表面極大點和極小點，此文件在本文的文件包里也提供了。按前文用hirsh_rho.vmd作圖后，將extrema.pqr載入VMD，將其顯示方式設為VDW并把Sphere Scale設為0.2，顏色用黃色，此時看到下圖，每個黃色小球都對應Hirshfeld surface上電子密度極大點位置。

若想查詢表面極值點處的電子密度數值，就在VMD里選Mouse - Query，然后點擊要考察的黃球的正中心，文本窗口就出現了它的index號，從0開始記。上圖箭頭所指的那個對應O-H...O作用的極值點的index為9，相當于從1開始記的編號為10。打開extrema.pqr，找到對應10號碳的下面這一行，倒數第3列的數值0.04228468就是此處的電子密度了，單位為a.u.。
HETATM   10  C   MOL A   1      -1.214   3.606   5.014   0.04228468   1.0000 C
以類似的方式查詢旁邊那個C-H...O氫鍵對應的Hirshfeld surface上的電子密度極大點，數值為0.01093053，可見作用顯著弱于O-H...O。

3.2 繪制指紋圖

這一節繪制HS分析中很常見的指紋圖。先按上一節的過程做HS分析并進入到主功能12的后處理菜單，然后輸入
20  //指紋圖與局部接觸分析
0  //開始分析
計算很快就完成了，然后看到一個新的后處理菜單，里面的選項不言自明。直接選0就會在屏幕上顯示完整的指紋圖，選1就會把指紋圖保存為當前目錄下的pdf文件（程序用pdf格式是因為它可以無損縮放、線條平滑），此pdf文件已經提供在了本文的文件包里（NAOB_HS.pdf）。當前的圖像如下

這個指紋圖里每一個小點都是一個Hirshfeld surface上的頂點，并且根據點在指紋圖上的分布密度進行了著色，越密的地方顏色越黃，越稀疏的地方顏色越紫。Multiwfn自動根據指紋圖上的最大分布密度設置色彩變化范圍的上限。Multiwfn畫指紋圖用的這種色彩過度方式稱為viridis，用Google搜圖功能一搜viridis就能找到相應的色條。

上圖在左下角有兩個顯著的尖兒（spike。兒化音明確體現這是個名詞），是圖像的特征區域，相當于“指紋”。靠左的那個尖兒的頂端大約是d_i≈0.7、d_e≈1.1位置，此處d_e顯著大于d_i，這是此體系具有氫鍵給體特征的體現（如上一節所示此體系的羧基氫確實是氫鍵給體）。為什么這個尖兒的存在能體現此分子存在氫鍵給體？因為這說明在Hirshfeld surface的這個區域，中心分子的原子到表面的最近距離（d_i）顯著小于周圍分子的原子到表面的最近距離（d_e），只有當氫鍵跨越這個區域，半徑很小的氫在表面內、半徑較大的重原子在表面外的時候才會出現這種情況，如下圖所示。

顯然不難理解前面的指紋圖中靠下的那個在d_e≈0.7、d_i≈1.1的尖兒體現的是中心分子作為氫鍵受體與周圍分子形成的氫鍵跨越了Hirshfeld surface。所以指紋圖說明當前研究的NAOB分子同時作為氫鍵給體和氫鍵受體。

指紋圖還蘊含了更多信息，但當前的指紋圖的中間區域密密麻麻一片，很難直接討論，進一步考察就需要利用下一節的局部接觸分析了。

順帶一提，如果你想改變指紋圖上的點的密度，可以在主功能12里做分析之前先選擇3 Spacing of grid points for generating molecular surface修改格點間距，格點間距越小則產生的Hirshfeld surface上的點就越多，指紋圖也會越密。

3.3 局部接觸分析與局部指紋圖

這一節演示一下怎么對中心分子特定的原子與周圍原子特定的原子進行局部接觸分析并繪制與之對應的指紋圖。作為例子，這里考察中心分子的氧與周圍分子的氫的局部接觸情況。

按照上一節說的進入到主功能12的后處理菜單中的20 Fingerprint plot and local contact analyses選項后，輸入以下內容
1  //設置內部原子范圍。默認是所有中心分子的原子都包括
[回車]  //代表對原子序號不做限制
O  //必須是氧元素
現在中心分子的所有氧原子都納入到了要考慮的內部原子范圍了。然后再輸入
2  //設置外部原子范圍。默認是所有周圍分子的原子都包括
[回車]  //代表對原子序號不做限制
H  //必須是氫元素
現在周圍分子的所有氫原子都納入到了外部原子要考慮的范圍了。

選擇0開始分析，算完后屏幕上顯示以下信息，告訴你當前考察的這種接觸面積是61.5 Angstrom^2，占Hirshfeld surface總面積的25.57%。
The area of the local contact surface is    61.504 Angstrom^2
The area of the total contact surface is   240.550 Angstrom^2
The local surface occupies   25.57% of the total surface

然后在后處理菜單選擇繪制指紋圖，看到下圖（對應本文文件包里的NAOB_HS_O-H.pdf）。此圖中只有對應于當前考察的局部接觸面上的頂點才被繪制為彩色，可見指紋圖中d_e≈0.7、d_i≈1.1的尖兒確實對應于內O與外H的接觸。

這種局部接觸還可以在立體結構圖上體現。在當前界面里選擇4 Export surface points to .pqr file in current folder，然后當前目錄下就產生了finger.pqr和finger_all.pqr，它們都提供在了本文的文件包里。前者記錄的是當前指定的局部接觸表面上的頂點，后者記錄的是完整的Hirshfeld surface上的頂點。它們可以用文本編輯器打開，可以看到每個表面頂點在pqr文件里用一個碳原子表示，倒數第三列記錄的是Charge屬性，當前被用來記錄表面頂點上被映射的函數數值，對當前來說就是準分子電子密度，單位為a.u.。

為了得到同時展現體系結構和局部接觸表面的圖，現在將NAOB_cluster.pdb載入VMD并恰當設置顯示方式，再把finger.pqr載入VMD，在Graphics - Representation里將它的Drawing Method設為Points并恰當設置Size，Coloring Method設為Charge（即根據Charge屬性著色），在Trajectory標簽頁里把色彩刻度設成與之前繪制等值面圖用的相同的0到0.012。確保Display - Rendermode已經設為了GLSL。在Graphics - Colors - Color Scale里把色彩刻度設為與之前相同的BWR，現在看到的圖如下，確實這些局部表面只對應中心分子的氧和周圍分子的氫的接觸。

也可以先用hirsh_rho.vmd照常繪制完整的Hirshfeld surface圖，然后再載入finger.pqr，把Drawing Method設為VDW并把Sphere Scale設為0.1，并且對Charge屬性用0到0.012色彩刻度范圍著色，此時看到的圖像如下，局部接觸部分以小球形式著重展現了。

下面再做個演示，對中心分子的芳環和周圍分子的芳環之間考察局部接觸，展現它們之間的pi-pi堆積。在GaussView里，或者在VMD里使用我在《在VMD中顯示原子序號的方法》（http://www.shanxitv.org/197）中提供的腳本，都可以在它們載入NAOB_cluster.pdb后顯示出原子序號，可以看到中心分子的芳環的原子序號是7-9,205,207,209。與中心分子芳環較近因而有pi-pi堆積作用的周圍分子的芳環有兩個，原子序號為11,13,15,112,114,116,201-203,356-358。因此，進入前述的Multiwfn的20 Fingerprint plot and local contact analyses選項后，輸入
1  //設置內部原子范圍
7-9,205,207,209  //原子序號范圍
[回車]  //對元素不做限制
2  //設置外部原子范圍
11,13,15,112,114,116,201-203,356-358  //原子序號范圍
[回車]  //對元素不做限制
0  //開始分析
4   //導出finger.pqr和finger_all.pqr

然后按照上一節的方法，用hirsh_rho.vmd繪制出Hirshfeld surface，再將finger.pqr載入VMD并恰當設置顯示方式。顯示體系結構的rep的選擇語句寫residue 0 1 11，其中residue 0是中心分子，residue 1和11對應與它有pi-pi堆積作用的上、下兩個分子。然后在Display - Material里把Translucent材質的Opacity改為0.5以降低透明度，并打開Angle-Modulated Transparency選項使得等值面立體感更強一些。之后看到的圖如下，可見小圓球把中心分子芳環上下兩側的pi-pi堆積對應的接觸區域都清晰展示了出來。

在Multiwfn里顯示出相應的局部接觸表面的指紋圖，如下所示。由于碳的原子半徑不小，因此pi-pi堆積對應的局部表面與表面內和表面外的碳原子都有一定距離，而碳原子間距離若太遠也不會有pi-pi堆積效應，因此這些表面頂點在指紋圖中的位置是d_i和d_e都不大不小的區域。

3.4 統計不同元素間的接觸面積

如果想得到中心分子與周圍分子的每一對元素之間的局部接觸面積，雖然通過前面演示的局部接觸分析可以手動實現，但有多少種元素組合就得操作多少次，很麻煩。因此Multiwfn提供了一次性完成的功能。進入前述的Multiwfn的20 Fingerprint plot and local contact analyses選項后，不需要定義內部外部原子，直接選擇選項3 Calculate contact area between different elements，就可以馬上得到以下信息

Inside element, outside element, their contact area (Angstrom^2) and percentage (%)
 H-H        47.330      19.676
 H-C        10.423       4.333
 H-N         0.521       0.216
 H-O        56.618      23.537
 C-H        12.821       5.330
 C-C        22.304       9.272
 C-N         2.614       1.087
 C-O         4.552       1.892
 N-H         0.876       0.364
 N-C         2.546       1.059
 N-N         0.464       0.193
 N-O         1.910       0.794
 O-H        61.504      25.568
 O-C         4.280       1.779
 O-N         1.816       0.755
 O-O         9.970       4.145

The same as above, but do not distinguish inside and outside elements
 H-H              47.330      19.676
 H-C/C-H          23.244       9.663
 H-N/N-H           1.396       0.581
 H-O/O-H         118.121      49.105
 C-C              22.304       9.272
 C-N/N-C           5.160       2.145
 C-O/O-C           8.831       3.671
 N-N               0.464       0.193
 N-O/O-N           3.727       1.549
 O-O               9.970       4.145

Area of total contact surface is   240.550 Angstrom^2

可見以上信息包含了每一對元素的分析結果。例如O-H        61.504      25.568這一行代表中心分子的氧元素與周圍分子的氫元素之間的局部接觸面積是61.504 Angstrom^2，占總面積的25.568%，這和3.4節我們專門算的結果完全一致。而后面輸出的比如H-O/O-H         118.121      49.105代表不管O和H誰在表面內、誰在表面外，兩種元素之間的接觸面積總共是118.121 Angstrom^2，占總面積的49.105%。

為了看起來更直觀，可以把上面的三列形式的數據拷到txt文件里，然后拖到Origin中導入，再繪制餅形圖，如下所示，一目了然。圖中只有占比相對較大的幾種接觸直接標注了數值。

3.5 Hirshfeld surface等值面的截斷問題

Multiwfn在構造Hirshfeld surface之前，需要先對一個矩形區域（稱為“盒子”）中均勻分布的格點計算你設定的體系片段的Hirshfeld權重。盒子范圍是對你定義的片段往各個方向延展一定距離來確定的，延展距離是位于最邊界的原子的范德華半徑乘以一個系數得到的。默認的系數值是1.7，一般來說夠大了，但如果Hirshfeld surface延伸到距離當前片段較遠的地方，導致超過了盒子范圍，則Hirshfeld surface就會在相應地方被截斷，等值面在那個地方看起來就會有窟窿。此時如果你想要讓等值面完整、完全封閉，顯然就需要增大系數。這一節拿富勒烯晶體舉個例子。

本文文件包里的C60.cif是C60富勒烯晶體結構文件。按照3.0和3.1節的做法摳團簇、繪制電子密度著色的Hirshfeld surface，會得到下圖。這里我在Graphics - Representation里選擇顯示等值面的那個Rep后在界面右下角Show旁邊的下拉框里選擇了Box+Isosurface，這樣除了等值面外，盒子范圍還會用細線同時顯示出來。由下圖可見，等值面上有個難看的窟窿，因為6個富勒烯之間有孔洞區域，Hirshfeld surface實際上會延伸到這里，但當前被尺寸有限的盒子截斷了。

當你發現存在這樣的窟窿，想通過增大延展距離來避免，就應當從主功能12的后處理菜單返回之前的菜單（即剛進入主功能12看到的菜單），然后輸入
4  //高級選項
1  //設置確定延展距離用的范德華半徑的倍數
2.3   //設一個比原本更大的值，數值可以反復嘗試。設得越大，盒子就越大，格點數就越多，計算耗時就越高、產生的cub文件尺寸就越大。由于之前的等值面只有一個小窟窿，所以盒子再大一點就夠了，因此當前嘗試的2.3比之前默認的1.7沒大太多
0  //返回
0  //開始計算

之后按照之前的步驟繪制Hirshfeld surface圖，如下所示，窟窿已經沒了。當前的Hirshfeld surface是個正十二面體，每個面都有塊白色且中間偏紅的區域，清楚地展現出中間的富勒烯和周圍富勒烯在此處有顯著的pi-pi堆積作用，由于這些區域離碳原子不遠因此電子密度不是很小。而表面上藍色部分都是晶體中的縫隙、孔洞區域，電子密度極低。

順帶一提，此體系的指紋圖如下所示。繪圖前用Multiwfn菜單里的選項3 Set range of axes適當加大了上限。可見散點的分布區域很窄，各處d_i和d_e數值都差不多大，涵蓋從中等到較大范圍。這體現出表面頂點基本都是位于中間分子和周圍分子之間的正中央區域，并且中間分子和周圍分子有的地方挨得近（d_i和d_e為中等大小），有的地方離得遠（d_i和d_e較大）

4 將Hirshfeld surface分析用于分子復合物

這一節示例HS分析用于展現孤立體系中的片段間相互作用。片段是指一批原子的集合，可以根據需要自由定義，比如既可以是分子復合物中的一個或多個分子，也可以是一個分子中的一個基團，也可以是一個過渡金屬配合物中的一個或多個配體，等等。對于孤立體系，由于感興趣的片段通常不是像晶體環境一樣被其它原子所完整包圍的，所以相應的Hirshfeld surface通常是開放的而不是完全封閉的。

筆者之前做了大量的和18碳環及其衍生物有關的研究工作，匯總見http://www.shanxitv.org/carbon_ring.html。其中研究了具有雙環特征的OPP分子與18碳環的結合問題，成果介紹見《8字形雙環分子對18碳環的獨特吸附行為的量子化學、波函數分析與分子動力學研究》（http://www.shanxitv.org/674）和《理論設計新穎的基于18碳環構成的雙馬達超分子體系》（http://www.shanxitv.org/684）。研究中使用了IGMH方法直觀展現了OPP與18碳環之間的弱相互作用，此例通過繪制電子密度著色的Hirshfeld surface圖也來展現一下。

上述研究中使用DFT優化出來的OPP結合一個18碳環的結構文件是本文文件包里的C18-OPP.pdb。將之載入Multiwfn，然后進入主功能0，會發現此結構文件中18碳環是斜著的，而由于HS分析用的盒子的邊框總是平行于笛卡爾軸的，因此直接對它產生Hirshfeld surface的話表面開口的地方也將是斜著的，看著很惡心。因此此例做HS分析之前還需要先用《Multiwfn中非常實用的幾何操作和坐標變換功能介紹》（http://www.shanxitv.org/610）里介紹的功能令18碳環平行于XY平面，如下所示。

啟動Multiwfn，載入C18-OPP.pdb，然后輸入
300  //其它功能（Part 3）
7  //對當前體系做幾何操作
11  //令一批原子擬合的平面平行于某個笛卡爾平面
225-242  //18碳環的原子序號
1  //平行于XY平面
現在可以選0在圖形界面觀看當前結構，會發現確實18碳環已經在XY平面上了。點RETURN關閉圖形窗口，然后接著輸入
-10  //返回
0   //返回到主菜單
12   //定量分子表面分析
1  //設置表面的定義
5  //Hirshfeld surface
225-242  //18碳環的原子序號
4  //高級選項
1  //設置定義盒子延展距離用的范德華半徑的倍數
2.3  //比默認值更大，實測不這樣的話會導致不好看的截斷
0  //返回
0  //開始分析
-2  //導出surf.cub
13  //導出mapfunc.cub

現在基于surf.cub和mapfunc.cub用hirsh_rho.vmd照常繪制Hirshfeld surface圖。這里我把VMD的四個光源都打開了，18碳環用CPK方式顯示，雙環OPP分子用Edgy材質。為了著色效果更鮮明，色彩刻度上限設為了比默認更小的0.01。此時圖像如下所示，可見白色和偏紅的條狀區域把18碳環與OPP的pi-pi作用最強烈的區域展現得挺清楚，此處電子密度比周圍相對更大。由于OPP的大環不是像18碳環一樣理想的圓形，形狀沒有完美匹配，因此它們之間的相互作用明顯不是處處均勻的。通過Hirshfeld surface的顏色可見，在大環末端，特別是兩個大環連接區域，大環與18碳環的pi-pi作用都很弱、電子密度很低。當前的Hirshfeld surface在18碳環上、下方都是開口的，因為在這里被盒子邊界截斷了，顯然這樣的截斷是理應有的，因為上下沒有其它原子了，原理上就不可能是封閉的。

下面這張圖是Comput. Biol. Chem., 101, 107786 (2022)中用Multiwfn+VMD畫的Hirshfeld surface圖，展現了配體和附近的蛋白質的氨基殘基間的相互作用情況，效果不錯，作用位點能看得比較清楚。

值得一提的是HS分析不僅限于用于考察弱相互作用，讓Hirshfeld surface跨越化學鍵也完全可以。在《一篇最全面介紹各種弱相互作用可視化分析方法的文章已發表！》（http://www.shanxitv.org/667）介紹的我的綜述文章中給出了一個九并苯共價二聚體的例子，這里我也繪制了其電子密度著色的Hirshfeld surface，其中一個九并苯被定義為了片段，色彩刻度用0-0.01，圖像如下所示。可見Hirshfeld surface正好嚴格平行于兩個九并苯、精確處于二者之間。在兩個九并苯之間形成共價鍵的區域及附近電子密度很大，顏色都為紅色，而其它白色區域體現了九并苯之間的明顯的pi-pi堆積作用。

5 Becke surface分析一例：DNA

距離原子非常遠的地方電子密度精確為0，這樣的地方是沒法計算Hirshfeld權重的，Hirshfeld surface也沒法在這個區域出現，這是HS分析的一個缺點。如果極個別情況下被研究的片段與其它部分的接觸面就是會牽扯到距離所有原子都很遠的區域，那么這個時候就必須用我提出的Becke surface代替Hirshfeld surface了，因為在電子密度為0的地方也是可以計算片段的Becke權重的。除了這種極個別情況外沒有用Becke surface的必要，圖像效果不會更好，而且對大體系耗時更高。

這一節就以DNA雙螺旋結構為例繪制做Becke surface分析。Multiwfn自帶的examples目錄下的DNA.pdb是一個DNA片段的結構文件，其中一條鏈的原子序號是1-319，它將被定義為片段。

啟動Multiwfn，載入DNA.pdb，然后輸入
12  //定量分子表面分析
1  //選擇表面的定義
6  //Becke surface
1-319  //第一個片段的原子序號
4  //高級選項
1  //設置定義盒子延展距離用的范德華半徑的倍數
0  //倍數設為0，使得盒子緊貼著邊界原子。因為感興趣的相互作用區域是在兩條鏈間，這樣減小盒子尺寸可以避免Becke surface過大、延伸到不感興趣的區域去
0  //返回
0  //開始計算
在筆者的i9-13980HX筆記本上，用16核并行，花了15分鐘。嫌慢的話可以用核很多的服務器。對于預覽目的，也可以把格點間距設大到比如0.5，耗時只有原本的約八分之一，但此時圖像會明顯更為粗糙。

導出mapfunc.cub和surf.cub并照常用hirsh_rho.vmd繪圖，得到的圖像如下所示，可見Becke surface把兩條DNA鏈間的接觸面非常直觀地展現了出來。在每一層的兩個堿基之間都有兩塊紅色區域，體現出兩個典型氫鍵的存在導致相應地方電子密度相對較大。仔細看的話，會發現有的兩塊紅色區域旁邊還有一小塊白色區域，這對應于C-H...O很弱的氫鍵。Becke surface的其它區域都是藍色，說明兩條鏈在其它地方并沒有值得一提的相互作用（注意不能描述為“沒有相互作用”，用詞必須嚴謹。重原子間的色散作用能即便到了10埃也沒完全衰減到0，兩條鏈的骨架之間的色散作用對相互作用能的貢獻不可完全忽略，只不過強度屬于“不值得專門一提”的程度）

按照HS分析的做法也可以顯示Becke surface的指紋圖，如下所示，可見左下方存在兩個尖兒，體現出每條DNA鏈既作為氫鍵給體也作為氫鍵受體和另一條鏈形成了氫鍵，和HS分析能展現的信息基本一致。但尖兒的具體位置不可能十分一致，畢竟Hirshfeld和Becke權重函數的定義方式差異顯著。

6 總結

本文介紹了HS分析的基本思想，并通過大量例子非常詳細介紹了Multiwfn做HS分析的方法和很多要點。可見HS分析用起來靈活方便，可以較好地直觀展現化學體系的特定片段與周圍其它原子之間的相互作用。同樣的目的用Multiwfn做基于電子波函數的IGMH分析也同樣可以達到，而且IGMH原理更為嚴格，通過sign(lambda2)rho函數對等值面著色還能區分相互作用類型，還能給出各個原子的貢獻量。HS分析可以作為IGMH分析的展現形式的補充，并且由于HS分析只依賴于原子坐標而且計算量很低，因此可以快速地用于很大體系，也不需要事先做量子化學計算產生電子波函數。注：只有原子坐標信息時也可以用Multiwfn做IGM分析，耗時也極低，但圖像效果明顯不及IGMH，參見《使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用》（http://www.shanxitv.org/621）。

最后提醒一下，使用Multiwfn做HS分析或Becke surface分析時請按Multiwfn啟動時的提示恰當引用Multiwfn原文，對于給別人代算的目的也需要明確告知對方這一點。

文/Sobereva@北京科音  2024-Feb-10

在《使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用》（http://www.shanxitv.org/621）中介紹的我在J. Comput. Chem., 43, 539 (2022)中提出的可視化分析弱相互作用的IGMH方法目前已經很流行，在《一篇最全面介紹各種弱相互作用可視化分析方法的文章已發表！》（http://www.shanxitv.org/667）介紹的綜述文章里對此方法也有全面透徹的介紹。

我在去年審了一篇Nature Chemistry的文章，看到作者通過對比對應不同類型相互作用的IGMH等值面來對不對稱催化的原因做了很好的解釋，我覺得很值得作為IGMH分析例子說一下，對一些人可能有啟發。現在這篇文章已經正式發表了：Nature Chemistry, 15, 1672 (2023)。

這是一篇實驗為主理論計算為輔的文章，研究對象是下圖的過程。圖的左上角所示的1是個氨基丙二酸類物質，在加入下圖右上角所示的有手性的SPA作為催化劑后，就可以得到脫羧后的結構2，是個有手性的氨基酸類物質。實驗發現，改變底物1中的AG基團會明顯影響產率和對映體過剩率e.e.，如圖所示AG為Bz時是75% e.e.，改為Bz1時80% e.e.，再引入OR取代基后e.e.更高，且R的選擇會進一步影響e.e.，如圖中方框部分所示AG=Bz7時e.e.最高，此時對應于R為CH2-CF2-CF3。下圖右邊示意了這個結構的底物可以與SPA形成顯著的弱相互作用，既有與SPA側鏈的pi-pi堆積作用，底物上的大量的F也與SPA的骨架的C-H有顯著的靜電和色散混合的吸引作用，這些作用可以顯著穩定化SPA對底物的結合。

上面那種AG=Bz7的結構在文中稱為S6，能夠自發脫羧變成下圖的Int-1結構。如下圖左側所示，Int-1可以和S6實現自催化反應（經歷TS-5a），也可以自己發生1-4質子轉移（經歷TS-5b），但勢壘都較高，不是研究的重點。重點是圖中所示的Int-1在SPA的作用下，會經歷一系列過渡態和中間體，最終達到氨基酸類似物結構6'，產生S和R手性的這個產物經歷的勢壘有很大差別。圖中可見TS-3是決速步，TS-3S的勢壘遠低于TS-3R，這是SPA催化作用下e.e.很大的關鍵原因。

由于反應物相同，因此TS-3S和TS-3R勢壘的差異體現在這兩個結構的能量差別上，差別最主要來自于SPA與Int-1分子的弱相互作用。為了說明這一點，文中用Multiwfn對這兩個過渡態結構做了IGMH分析，給出了下圖，對兩個分子間存在的相互作用類型做了標注。其中N-H...O無疑是挺強的氫鍵，C-H...pi和C-H...O/F都是很弱的氫鍵，苯環間存在pi-pi相互作用，Bz7基團上的一堆F的孤對電子與苯環有色散主導的相互作用。圖中標注Bonding的地方是Int-1與SPA之間發生氫轉移而一定程度成鍵的地方，對TS-3S和TS-3R這個成鍵作用不會有什么明顯差別。

上面兩張IGMH圖的等值面很多，在表面上看起來不太好對比相互作用強度，但如果像本文一樣，把每種相互作用一一羅列出來做細致對比，是完全可以對總作用強度進行區分和解釋的。如上可見，通過仔細觀察等值面可發現TS-3S比TS-3R對應C-H...X作用多出來2個，并且TS-3S的pi-pi堆積作用區域的面積明顯大于TS-3R的，所以前者被標注為strong而后者被標注為weak。

為了能令弱相互作用情況看得更清楚，文中還在分子結構圖上直接用虛線對IGMH呈現的信息做了明確的標注，如下圖所示，字母對應上面的列表里的項，諸如c1、c2...對應不同的C-H...X相互作用。從這個簡化的圖上能更清楚地看到TS3-S中SPA的C-H與底物的Bz7基團上的F之間的相互作用比TS3-R中的更豐富，這一定程度解釋了為什么底物接上F原子較多的Bz7基團時被SPA催化時可以達到較高的e.e.。此外，下面這張圖上還標記了d: C-H...C-H作用的位置，顯然這是普通色散作用，這個作用在上面的IGMH等值面上肉眼不太好辨別（在VMD程序里旋轉等值面時才容易識別），而像下圖這么用箭頭標記一下就容易令讀者看清楚出現在什么位置了。

這篇文章體現出基于IGMH分析能清晰直觀地對分子間相互作用導致勢壘存在的差異和由此帶來的實驗現象予以解釋，是很好的IGMH應用范例，值得在其它研究中借鑒。此文的利用弱相互作用實現e.e.的控制也是新穎的催化設計的思路。

文/Sobereva@北京科音  2024-Feb-9

0 前言

量子化學研究中在優化完極小點結構然后做振動分析時遇到虛頻是家常便飯的煩人的問題。之前我專門寫過博文《Gaussian中幾何優化收斂后Freq時出現NO或虛頻的原因和解決方法》（http://www.shanxitv.org/278）講了怎么消虛頻。除非程序存在bug，否則原理上不存在按上文處理完全解決不掉的虛頻。一個現實問題是：雖然靠各種嘗試、折騰，最終確實總能消掉虛頻，但這個過程對于大體系來說過于浪費計算時間，也耗費研究者的精力。那么是否帶有虛頻的結構就一定不能拿來算能量相關的問題呢？是否虛頻總是要不惜一切代價消除呢？此文就談談我個人的看法，給出理論分析，并且通過一個例子對我的觀點予以驗證。

本文牽扯到熱力學數據計算的各種常識知識，沒有這些常識的話先看《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》（http://www.shanxitv.org/552），以及Shermo程序手冊的附錄，里面有計算原理的簡潔易懂的概述。在北京科音基礎（中級）量子化學培訓班（http://www.keinsci.com/workshop/KBQC_content.html）里對熱力學量的計算有非常全面深入的介紹，歡迎關注和參加。

本文說的小虛頻和過渡態必然存在的虛頻沒有絲毫聯系，初學者絕對不要把本文內容胡亂套用到過渡態的情況。

1 導致小虛頻的原因

為了避免誤解，首先回顧一下優化極小點后振動分析時出現虛頻的各種常見原因：
(1)振動分析時用的影響勢能面的設置與優化時的不同。這通常是沒有基本常識的外行人才會犯的低級錯誤。注：JCTC, 17, 1701 (2021)中提出的single-point Hessian (SPH)方法通過對勢能面添加偏移勢，從而允許振動分析用的計算級別低于優化用的級別來節約時間，但這個做法非主流、經驗性強，本文不涉及
(2)初始結構的點群高于實際極小點的點群，且如果優化后的結構維持住了初始的對稱性，振動分析后就會發現破壞對稱性的虛頻，即按照虛頻方式移動結構會導致結構的對稱性下降。有時候這種有虛頻的結構恰對應于過渡態
(3)程序bug。比如Gaussian 09 D.01的DFT-D3色散校正對解析Hessian的貢獻的計算有bug，容易導致算勢能面平緩體系出現虛頻。這點在《DFT-D色散校正的使用》（http://www.shanxitv.org/210）我專門做了說明
(4)某些方法對勢能面引入明顯數值噪音。如Gaussian中的SMD溶劑模型的這個問題就比較明顯
(5)幾何優化收斂得不精確。幾何優化收斂限過于寬松或Hessian的質量不太好會導致這種問題。
(6)積分格點質量不夠好。這個問題對于M06-2X等對積分格點敏感的泛函尤為顯著。關于量子化學程序用的原子中心的積分格點參考《密度泛函計算中的格點積分方法》（http://www.shanxitv.org/69）。如果是CP2K等利用平面波的第一性原理程序，用的均勻分布的格點依賴于平面波截斷能的設置。

以上前四種情況導致的虛頻是一定要解決掉的，而且本來解決就非常簡單。只要不犯低級錯誤自然就沒有(1)的情況出現，解決(2)按照虛頻模式調節結構后進一步優化即可，解決(3)就換成沒bug的版本或其它計算程序，解決(4)可以換成其它的能在優化+振動分析中實現基本相同目的而沒嚴重數值噪音的方法，如SMD換成IEFPCM，或者用SMD時結合SAS方式定義溶質孔洞表面。對于(5)和(6)情況帶來的虛頻，只要計算能力和精力允許，能消盡量要消，但如果由于計算資源或精力所限實在太難消掉、被迫要容忍虛頻的存在，那么后續的能量相關問題的計算就是后文討論的內容了。

對于柔性大分子、弱相互作用復合物這樣勢能面的某些維度非常平緩的體系，在做完幾何優化和振動分析后很容易發現存在大小非常小的虛頻，如波數為9.8i、27.3i cm^-1這種，普遍小于50，這一般是前述(5)或(6)情況所導致的。此類虛頻的振動模式通常對應于甲基旋轉運動、柔性骨架的變形運動、分子間相對運動。這類(5)或(6)情況所導致的虛頻以下簡稱為“小虛頻”。但要注意也不是所有波數小于50的虛頻都屬于這種情況，前述的其它四種情況也完全可能導致數值很小的虛頻，所以絕對別光拿虛頻大小來判斷情況。

2 有小虛頻時能量相關問題的計算

當遇到小虛頻時，顯然最完美的做法（具體選擇視情況而定）是用諸如更嚴格的幾何優化收斂限、更好的Hessian、更好的積分格點精度進一步做優化和振動分析，使得虛頻完全消除后再計算能量相關問題。但是對于普通泛函算幾百個原子有機體系，很可能算一次Hessian就要花一整天時間，幾何優化一步也可能花費一個小時左右，最終把虛頻完全消掉可能得花幾天甚至一個禮拜，十分痛苦，甚至不可能實現。那么，在后續算能量相關問題時怎么考慮？下面根據被計算的量專門討論。

2.1 計算電子能量

諸如Gaussian這樣的主流量子化學程序默認的幾何優化收斂限是以受力和位移作為判據的。在默認收斂限下優化正常完成后，即便結構還沒有與實際極小點特別精確一致，但通常相差很小。而且當前的最大和平均受力肯定已經很小了（因為收斂限考慮了它們），由于受力對應于能量梯度的負矢量，所以當前結構稍微再挪一點達到精確的極小點的過程中能量變化肯定非常小。即曰，在有小虛頻的情況下，計算的電子能量一般是可以用的，也即結果和消掉虛頻后再算的能量相差極小。

2.2 計算0 K下的內能U(0)

U(0)是電子能量與零點振動能（ZPE）之和。在諧振近似模型下，每個振動模式對ZPE的貢獻是1/2*h*ν，h是普朗克常數，ν是振動頻率。頻率很低的振動模式對ZPE的貢獻是非常小的，每波數的振動頻率對ZPE的貢獻僅為5.98 J/mol，因此諸如30波數的振動模式只對ZPE有0.18 kJ/mol（0.04 kcal/mol）的貢獻，數量級遠遠低于理論方法、基組、溶劑模型等因素帶來的誤差。順帶一提，Shermo計算熱力學量的時候如果把settings.ini里的prtvib設為1，每個振動模式對熱力學量的貢獻都會明確輸出。

主流量子化學程序，以及Shermo程序在默認情況下，在計算振動對熱力學量的貢獻時都是不考慮虛頻產生的貢獻的，因此小虛頻對熱力學量的貢獻會被直接忽略掉，這不會對ZPE的計算造成實際問題。因為如果把這個小虛頻消了，從而產生一個實頻，其頻率肯定也是非常低的，對ZPE的貢獻可忽略不計。也就是不管消不消這個虛頻，算出的ZPE都沒多大差別。不過，實際情況遠沒這么簡單，因為在虛頻消掉前后，由于二者所在勢能面的位置的些許不同，或改用不同檔次積分格點導致頻率計算精度的差異，其它的實頻模式也會多多少少受到影響，不過頻率值受影響相對明顯的主要是低頻模式（可能影響幾個甚至10個左右波數，而高頻模式受影響一般小于2波數），它們對ZPE的貢獻恰恰又很小，所以虛頻消掉前后算的ZPE雖然未必差異小到諸如前面說的0.2 kJ/mol左右這種程度，但至少也不會大，比如對柔性大體系頂多也就影響個1-2 kJ/mol。

如上所述，鑒于是否小虛頻對電子能量和ZPE的影響都不太明顯，因此算電子能量和U(0)時一般可以容忍一個或多個小虛頻的存在。

2.3 計算從0 K升溫到當前溫度對內能的影響

研究U(0)沒太大實際意義，畢竟實際化學過程也不可能在0 K發生。計算有限溫度，也就是高于0 K時候的熱力學量，就需要考慮從0 K升溫到當前溫度對內能的改變量，其中振動產生的貢獻下面用Uvib(0→T)表示。首先要知道一點，低頻對Uvib(0→T)的貢獻是明顯高于高頻模式的。這點從下面的J. Phys. Chem., 100, 16502 (1996)的圖1可以明顯看出來，圖中ΔHvib(T)和這里說的Uvib(0→T)是一碼事，圖是對T=298.15K的情況繪制的。

雖然低頻模式的貢獻大，但是由圖中放大的0-100波數的區間可見，在小幾十波數范圍以內，頻率值對Uvib(0→T)的影響雖然不可忽視，但不算多大，都在2.0-2.4 kJ/mol范圍內。這也就是說，由于存在一個小虛頻導致少考慮一個實低頻模式，會給Uvib(0→298.15K)帶來-2.2 kJ/mol左右的誤差，這雖然不算小，但還算可以湊合接受。如前所述消虛頻后還可能會對其它低頻模式產生不可忽視的影響，這對Uvib(0→T)的改變一般也不至于太大。

綜上所述，若帶著小虛頻計算U(T)，會造成一定低估。對于有一個小虛頻的情況，低估程度在0.5 kcal/mol左右。如果有多個小虛頻，那還是別計算U(T)了，除非T遠低于常溫，或者結合后文說的把小虛頻當成實頻的做法。

2.4 計算熵、自由能

低頻對熵的貢獻遠大于高頻，而且在很低頻區間內，頻率越低貢獻明顯越大。在包括Gaussian在內，大多數量子化學程序（ORCA除外）計算熱力學量用的是傳統的rigid-rotor harmonic oscillator (RRHO)模型，振動熵的計算完全基于諧振近似模型。此時在低頻區間內，隨著頻率趨近于0，振動模式對熵的貢獻猛增，如下面的Chem. Eur. J., 18, 9955 (2012)文中圖2的虛線所示。Grimme在這篇文章中還提出了一種quasi-RRHO模型（QRRHO，后來也叫modified RRHO、mRRHO），以下稱為QRRHO(Grimme)，它將自由轉子和RRHO用的諧振近似模型算的振動熵做插值，頻率越接近0自由轉子的權重越大、頻率越高諧振近似的權重越大，如下圖實線所示。QRRHO(Grimme)對高頻模式算的熵貢獻和RRHO完全一樣，而對于低頻模式（特別是波數在100以下的）考慮得比RRHO明顯更合理。由下圖還可見QRRHO(Grimme)算的低頻模式的熵比RRHO整體小得多，沒有RRHO的嚴重高估的問題。此外，QRRHO(Grimme)下低頻模式的振動熵對頻率值的敏感性低于RRHO很多。

直接計算熵的用處不大，算熵的主要目的是用來算自由能，因為自由能里有-T*S項。顯然低頻對熵的貢獻計算的準確度顯著影響到自由能的計算準確度。假設體系有一個小虛頻，并且在消掉之后會多一個20 cm^-1的實頻，并且用的是QRRHO(Grimme)算298.15 K的熵，那么在有虛頻的結構下由于缺失了這個實頻，根據上圖可知會造成熵的誤差約為-18 J/mol/K，對自由能造成的誤差約為-298.15*(-18)/1000=5.4 kJ/mol，略高于1 kcal/mol，這就不小了，都能趕上高質量ωB97M-V泛函結合大基組算普通有機反應的誤差了。不過實際中小虛頻帶來的誤差未必那么大，因為有和沒有小虛頻的結構下實低頻模式的頻率也往往有不可忽視的不同，例如下面第3節的例子在有虛頻的結構下實低頻模式的頻率比精確極小點結構下的整體偏小（勢能面形狀原因所致），導致高估了不少實低頻模式的熵，藉由-T*S項進而對自由能造成低估，于是和少考慮一個低實頻對自由能的高估產生了很大程度的相互抵消，下面管這稱為“熵的誤差抵消”。

根據以上討論，有小虛頻的情況下計算熵和自由能都有一定風險。如果就有一個小虛頻，算常溫的情況往往還說得過去，而如果有多個小虛頻，或計算高溫的情況（注意T越大-T*S越大，熵的誤差造成自由能的誤差越大），我就不建議冒險帶著虛頻了。而且就算非要帶著小虛頻算，也至少應當用QRRHO模型，而切勿用RRHO模型，因為RRHO算的振動熵對低頻頻率值的敏感程度顯著高于QRRHO，因而用有虛頻的結構算的熵的誤差上限比QRRHO的明顯要大。

要注意QRRHO模型不止Grimme提出的這一種。如Shermo手冊附錄部分所介紹的，Truhlar的QRRHO是把所有小于一定閾值（通常為100 cm^-1，可以由Shermo的ravib參數控制）的低頻統一提升為這個閾值再算它們對Uvib(0→T)和熵的貢獻，這個做法的物理意義不甚明確，沒有QRRHO(Grimme)的插值做法優雅。目前還有一種QRRHO是Minenkov等人在J. Comput. Chem., 44, 1807 (2023)中提出的，它對振動熵部分的考慮和QRRHO(Grimme)一樣，但還同時將這種自由轉子與諧振近似做插值的思想同時應用到了計算振動對內能的熱校正量Ucorr上（此時沒法單獨計算ZPE），原理更理想，實際結果還略好一點點。這些QRRHO模型在最新版Shermo中都是支持的，通過ilowfreq參數控制。

2.5 將小虛頻視為實頻的做法

在Angew. Chem. Int. Ed. 2022, e202205735中的3.2節Grimme等人認為存在小虛頻的情況下，只要使用QRRHO模型并且將小虛頻視為實頻處理，比如18i cm^-1直接當做18 cm^-1，則得到的熵就是可以用的、小虛頻的存在就不再是個問題。這種做法表面看起來似乎有點莫名其妙、太任性了，但確實有一定道理，可以這么理解：在消除小虛頻后，必定會得到相同數目的波數也很小的實頻模式，它們都對熵有明顯貢獻。相對于完全忽略掉這些小虛頻的貢獻而導致明顯低估熵，干脆把他們當做實頻處理，從而“湊”出來相應數目的實低頻是更好的做法，而且本身在QRRHO模型下，數值在十幾、幾十波數區間的頻率對熵的貢獻也都相差不懸殊，因此也用不著對這么人為湊的實頻頻率太較真。由于如前所述，Uvib(0→T)也對低頻數值不敏感而低頻的貢獻又沒小到可以隨便忽略的程度，所以計算Uvib(0→298.15K)時也使用這種處理是有道理的。從Shermo 2.5版開始，可以通過settings.ini里的imagfreq選項設置視為實頻的虛頻閾值，比如設為50，就代表大小在50 cm^-1以內的虛頻都當做實頻處理來計算振動對各種熱力學量的貢獻，此時從Shermo程序輸出的頻率信息中也會看到已經沒有虛頻了。

然而，這種將小虛頻視為實頻的做法目前很非主流，其思想在我來看過于主觀和樂觀、把問題想得太簡單了，目前也缺乏系統性的檢驗，我個人不建議輕易使用。這種做法有一個明顯問題：如前所述，有小虛頻和消掉虛頻的兩種情況相差的不僅僅是一個小虛頻轉變成一個小實頻，許多其它低頻模式的頻率也會有不小變化，而且有時是整體明顯增加的。不做這個虛頻到實頻的轉換的情況下能夠有上一節說的熵的誤差抵消效應，而這么做了之后，忽略虛頻模式導致對熵的低估問題確實很大程度解決了，但在有小虛頻結構下其它頻率往往偏低而導致熵的高估帶來的誤差就沒處抵消了。此外，如果當前體系的低頻還有很多都是10 cm^-1 左右這種特別低的，這個范圍附近即便是QRRHO(Grimme)或QRRHO(Minenkov)模型算的熵也對頻率很敏感，問題會更大。

所以，將小虛頻視為實頻的做法雖然確實對一些情況更好地計算熵和自由能是有幫助的，但也很可能起到負面效果。所以這絕不是普適、穩健、理想的做法，更千萬不能一看到存在小虛頻就總是想靠這個做法來完全無視虛頻。

3 有虛頻時算熱力學量的實例

筆者最近通過Gaussian 16在ωB97XD/6-311G*級別研究的一個以pi-pi堆積方式形成的共128個原子的三分子復合物，在優化和振動分析后發現存在6.44i cm^-1的小虛頻，之后進一步優化時用opt=calcfc關鍵詞提供精確的初始Hessian，優化完再做振動分析就發現沒虛頻了，此時最低的頻率成了10.35 cm^-1。正好這個例子可以拿來考察消虛頻前后計算的熱力學量的差異，由此檢驗帶著那個小虛頻時以不同方式算的熱力學量有多大誤差。要強調的是，僅僅這么一個測試，顯然不足以充分證明某種做法理想，因為這必須通過大量體系做統計分析，但至少誤差較大的情況可以用來體現相應做法不夠穩健和普適。

先看一下有虛頻相對于無虛頻時各個實頻的頻率變化情況。先把實頻頻率按照由低到高排列，用有虛頻時候它們的頻率減去無虛頻時候的頻率作為縱坐標值，而兩種情況頻率的平均值作為橫坐標，作的圖如下所示。可見，至少對于此例，有虛頻時數值較低的實頻是普遍明顯偏低不少的，而>200 cm^-1的相對而言的中、高頻則受到的影響很小。這種情形雖然不能說是普遍情況，但至少也是挺常見的一種。

此例有虛頻相對于無虛頻的情況計算的電子能量僅僅相差0.04 kcal/mol，小到完全可以忽略不計。用Shermo程序基于Gaussian振動分析輸出文件計算的有虛頻相對于無虛頻時的ZPE和氣態標況下的熵、內能、自由能熱校正量的差異如以下表格所示。RRHO以及Truhlar、Grimme、Minenkov三種QRRHO分別對應Shermo的ilowfreq=0/1/2/3的情況，imagreal=0和50相當于忽略那個虛頻以及把那個虛頻當大小相同的實頻考慮的情況。

首先看常規的忽略虛頻貢獻的情況。由表格可見：
? ZPE：有虛頻時和無虛頻時ZPE相差僅-0.3 kcal/mol，屬于通常可以接受的差異
? 熵：由于誤差抵消的巧合，恰好RRHO下有小虛頻時的熵和無虛頻時的差異很小。QRRHO(Truhlar)的情況差異挺大，在于它把實低頻都統一當成了100 cm^-1，這部分在有和沒有小虛頻情況之間幾乎都抵消了，而有小虛頻時相當于少考慮了一個100 cm^-1振動模式對熵的貢獻，因此造成熵的嚴重低估。QRRHO(Grimme)和QRRHO(Minenkov)對熵的部分處理相同，由于前述的熵的誤差抵消效應導致有和沒有小虛頻時熵的差異不算太大
? 內能的熱校正量Ucorr：等于ZPE+Uvib(0→T)。由于Uvib(0→T)對低頻的具體數值不太敏感，有無小虛頻情況下實低頻的貢獻變化很小，而有小虛頻時由于缺失了一個低實頻對Uvib(0→T)的貢獻，導致有小虛頻時Uvib(0→T)偏低不少，Ucorr因此比ZPE偏低得明顯更多。QRRHO(Minenkov)用的計算Ucorr的方式和其它模型不同，對于此例使得有小虛頻時Ucorr偏低程度比其它模型更小，僅為-0.5 kcal/mol。因此在計算有小虛頻情況下的U(T)時，用QRRHO(Minenkov)可能會比其它模型更有優勢。
? 自由能的熱校正量Gcorr：相對于其它模型，有小虛頻相對于無虛頻時Gcorr的差異在RRHO下是最大的，這在于有無小虛頻時RRHO算的S的差異對此例恰好頗小，導致-T*S項的差異才0.03 kcal/mol，因此沒能充分對Ucorr的差異產生抵消效應。而各種QRRHO模型下Gcorr的差異明顯較小，這在于這種情況下S的差異不小，而且-T*S項的差異的符號和Ucorr的差異相反，因此相互抵消了很多。其中Gcorr差異最小的是QRRHO(Minenkov)，才-0.11 kcal/mol，體現出它可能適合作為有小虛頻情況下算自由能的優先選擇。

再來看將小虛頻當實頻處理的情況。根據以上表格可見，這種處理使得Ucorr在有和沒有小虛頻情況下計算的結果差異減小很多，因此對于算內能的目的，這種做法有一定意義。但是，除了QRRHO(Truhlar)外，這個做法卻造成了有小虛頻時算的S和Gcorr與無虛頻時的差異變得顯著更大，對RRHO最為嚴重，對QRRHO(Grimme)和QRRHO(Minenkov)問題也挺嚴重，這來自于有小虛頻時實頻模式偏低造成的對振動熵的高估無法與因為少考慮一個實低頻導致的低估一定程度相抵消。把小虛頻當實頻處理時結合QRRHO(Truhlar)對此例倒是很不錯，有小虛頻相對于無虛頻時的熵只相差0.27 cal/mol/K，這在于此種QRRHO將低頻全都上拉到100 cm^-1，此時又把小虛頻轉化成了小低頻，因此無虛頻和有虛頻時相當于有同等數目的100 cm^-1的實頻，自然兩種情況下算的熵很接近。

4 總結

根據以上討論和實際體系的測試可知，對于因(5)和(6)情況造成有小虛頻的情況，算電子能量和U(0)沒大問題，不是非得消虛頻。

對于計算Ucorr或與之相差RT的焓的熱校正量Hcorr，用QRRHO(Minenkov)或QRRHO(Grimme)并同時把小虛頻當實頻處理是較好選擇，不會因為存在小虛頻的問題而帶來太大誤差。

有小虛頻時最難搞的是熵，以及溫度不很低時候Gcorr的計算，難點在于小虛頻消掉前后實低頻模式頻率也可能受到不小影響且影響情況不易估計，然而即便用QRRHO(Grimme或Minenkov)時低頻的具體數值對熵的影響也不算太小，所以明顯不是把小虛頻當成實頻處理以避免少考慮一個實低頻模式就能較好解決的，而且如第3節的測試可見，這種做法還可能造成熵和Gcorr在有小虛頻時計算結果更差。所以當你對熵、Gcorr要求精度高時，小虛頻要盡量消。實在搞不定的話，那就姑且用原理上最好的QRRHO(Minenkov)，并且建議不用小虛頻當實頻的處理，畢竟此做法非主流而且起到正面效果的可能性很有限。雖然把小虛頻當實頻并結合QRRHO(Truhlar)的情況下小虛頻和無虛頻時熵的差異很小，但在J, Comput. Chem., 44, 1807 (2023)的測試中QRRHO(Truhlar)算團簇反應自由能表現得遠不如QRRHO(Minenkov)，甚至比RRHO沒明顯改進，所以我也不推薦。

文/Sobereva@北京科音   2024-Jan-29

0 前言

環狀化學體系具有不同程度的環張力，環張力和環狀體系的諸多問題緊密相關，諸如生成焓、穩定性、合成難易程度等。環張力能（strain energy, SE）是衡量環張力大小的最直接的指標，它對應當前環狀結構轉變為假想的沒有環張力的結構過程中能量變化的負值。

北京科音自然科學研究中心（www.keinsci.com）的盧天等人之前受邀發表了一篇專門精確考察碳環（cyclocarbon）體系環張力能的文章，文中對環張力的計算思想有詳細介紹并給出了諸多討論：
Accurate theoretical evaluation of strain energy of all-carboatomic ring (cyclo[2n]carbon), boron nitride ring, and cyclic polyacetylene, Chin. Phys. B, 31, 126101 (2022) DOI: 10.1088/1674-1056/ac873a（不方便獲取者也可以看預印版https://doi.org/10.26434/chemrxiv-2022-v8w9h，內容基本一樣）

下面，將基于此論文中的內容簡要說明環張力能的計算原理并給出具體例子。感興趣的讀者請務必仔細閱讀論文原文以了解更多信息，也推薦以類似方式計算其它體系的環張力能時引用此論文。另外，讀者如果對碳環體系感興趣，很建議查看同作者之前發表的此類體系的各種研究文章和相關博文：http://www.shanxitv.org/carbon_ring.html。

1 同聯反應法

設計同聯（homodesmotic）反應是計算環張力能的很常用方法。這個方法設計一個化學反應，既使得環被打開使得環張力被完全釋放，同時又盡可能不影響原本的成鍵情況，顯然這個反應能的負值就對應了環張力能。

例如，要計算氧雜環丁烷的環張力能，就可以設計以下同聯反應。紅色部分對應具有顯著環張力的氧雜環丁烷，將它插入到甲乙醚所示位置就變成了右邊的產物，顯然此時環張力就被完全釋放掉了。與此同時，反應物和產物的成鍵特征并未發生變化，例如下圖紅字的第一個碳在反應前后始終是sp3雜化狀態、始終連著一個氧原子和一個CH2，因此化學環境變化導致的成鍵特征的極輕微改變對這個反應能的影響可以忽略不計，也即這個反應能幾乎只體現出環張力的釋放。

同聯反應法計算環張力能的缺點是反應的設計并不唯一，比如以上反應的甲乙醚也可以是甲醚、乙醚等等，結果或多或少會有點差異。顯然反應應當設計得盡可能令各個鍵的特征在反應過程中不發生改變，但這并非總能很好實現，后文提到的碳環體系就是一例。

2 同聯反應法計算實例：[6]CPP

這一節示例使用同聯反應法計算[6]CPP的環張力能。此體系的幾何結構，以及設計的同聯反應如下所示。可見是把[6]CPP插入到了聯苯中間變成了直線的八聯苯，從而完全釋放了環張力。

上圖涉及的三個分子都在ωB97XD/6-31G*級別下用Gaussian 16程序進行了幾何優化，并做振動分析確認了無虛頻，輸入輸出文件都在http://www.shanxitv.org/attach/698/file.rar里。基于輸出文件里最后一次輸出的電子能量，用反應物能量之和減去產物能量，就得到了環張力能：627.51*(-463.1444911-1385.7232492+1849.0260518) = 99.3 kcal/mol。這個值比C-C單鍵的鍵能還要大一些，可見其環張力是相當強的。

如果想得到更精確結果，可以基于優化完的結構在更好的級別下算高精度單點再求差，比如基組可以改用明顯更好的def2-TZVP。這對結果可能造成幾kcal/mol程度的影響。

3 外推法

如果被考察的環狀體系是由重復單元構成的，比如前面的[6]CPP是苯環作為重復單元，就也可以用這一節介紹的外推法計算其環張力能。這個方法計算過程如下

(1)計算從小到大不同重復單元數（n）的環的每單元的能量E/n，這里E是整體能量。顯然E/n是依賴于n的，n較小時環尺寸較小，自然環張力大，平均到每個單元的環張力能也大
(2)用Origin之類程序通過恰當的函數（一般是二次函數）擬合E/n與n之間的解析關系，并外推得到n無窮大時的E/n，這對應于沒有環張力時候的每單元能量
(3)將當前考察的環狀體系的E/n減去n無窮大時的E/n，就得到了當前體系的每單元環張力能（SE/n）
(4)將當前體系的SE/n乘以當前的體系的n即是當前體系的環張力能

相較于同聯反應法，外推法的優點：
?原理非常嚴格
?沒有同聯反應法設計反應的任意性，并可以用于同聯反應法不適合的體系
?可以順帶得到環張力能與n的解析關系，由此可以直接預測出任意重復單元數的環張力能
缺點：
?需要計算從小到大不同尺寸的環狀體系，總耗時明顯高于同聯反應法
?計算步驟略繁瑣，得自己做擬合

如果你只是要考察當前一個環狀體系的環張力能，而且可以設計出很合適的同聯反應，那么就沒必要用昂貴且略麻煩的外推法。

4 外推法計算實例：18碳環

2019年首次在凝聚相實驗中觀測到的18碳環分子由18個sp雜化的碳依次相連構成環狀。由于sp雜化的碳明顯傾向于形成180度鍵角，因此18碳環必然有顯著的環張力，具體環張力能是多少非常值得研究，這正是前述Chin. Phys. B, 31, 126101 (2022)文章研究的主要內容。

此文中指出碳環類體系不適合同聯反應法計算環張力能。碳環類體系是長、短C-C鍵交替構成，兩個碳是一個重復單元，若用同聯反應法計算化學式為C_2n的碳環的環張力能，可以設計如下反應式，相當于把碳環的每個重復單元插入到了一定長度的碳炔里。碳炔是sp雜化的碳形成的鏈狀分子，兩端由氫封端。PS：對碳炔感興趣者建議閱讀《氫封端碳鏈H-(C≡C)n-H (n = 3-9, 15)的電子光譜的尺寸依賴性：性質分析及對碳炔的預測》（http://www.shanxitv.org/679）

然而，如下圖可見，將18碳環的重復單元插入到不同長度（通過m體現）的碳炔里，按照同聯反應算出來的環張力能有明顯區別，因此靠同聯反應法難以得到嚴格的18碳環的環張力能。

碳炔和18碳環都具有C-C鍵長、短交替的特征。文中還發現，不管碳炔的m取多少，處于它中央的較長的C-C鍵和18碳環中的較長的C-C鍵之間的鍵長差異都不可忽視，對較短的C-C鍵也是如此。因此原理上沒法設計一個同聯反應使得18碳環的C-C鍵特征在轉移到產物時幾乎完全不發生改變，這也必然導致靠同聯反應法不可能對此體系得到完全精確的環張力能，也即成鍵特征差異造成的能量改變會不可避免地摻入到計算出的環張力能中。

由于發現了以上問題，文中使用外推法計算碳環的環張力能。為了追求盡可能好的精度，文中用ORCA程序在精度很好的DLPNO-CCSD(T)/cc-pVTZ級別下對n=7到n=24的碳環都計算了單點能（基于ωB97XD/def2-TZVP優化的結構），ORCA輸出文件在http://www.shanxitv.org/attach/699/file.rar里都提供了。每重復單元能量（E/n）與n的關系如下所示。由于n較小的碳環的電子結構特征有特殊性、和較大的碳環存在一定差異，因此文中只對n>=13的數據點做了二次函數的擬合，由下圖的紅色曲線可見擬合質量超極好，R^2近乎精確為1。

令擬合公式E/n = 0.8844/n^2 - 76.00997中的n為無窮大，就得到了無環張力對應的無窮大的碳環的每單元能量-76.00997 Hartree。每單元環張力能SE/n = E/n - 76.00997，于是有SE/n = 0.8844/n^2，也即碳環類體系的環張力能SE = 0.8844/n Hartree。18碳環的n=9，代進去就得到了其環張力能0.8844/9 Hartree = 61.7 kcal/mol。

值得一提的是，如果令前述同聯反應法中碳炔的m設得很大，從而使得環張力能隨m收斂，最終算出來的結果是59 kcal/mol，可見和嚴格的外推法的結果相比誤差并不可忽略。

前述論文中也測試了用ωB97XD/def2-TZVP的能量結合外推法算碳環的環張力能，結果和上面用的昂貴的DLPNO-CCSD(T)/cc-pVTZ能量的結果沒多大差別。這在一定程度上體現出ωB97XD/def2-TZVP級別算碳環類體系很不錯。

文中還基于不同尺寸碳環的環張力能和鍵角之間的關系，推導出了適合碳環體系用的C-C-C鍵角力常數56.23 kcal/mol/rad^2。

此外，文中還對以下所示的18碳環的等電子體硼-氮環，以及具有單套pi共軛特征的環聚乙炔的環張力做了計算，發現硼-氮環類體系的環張力明顯小于碳環類，而環聚乙炔類體系的環張力則略大于碳環類。文中利用Multiwfn程序（http://www.shanxitv.org/multiwfn）做波函數分析從電子結構角度清晰解釋了它們的環張力與碳環之間存在差異的原因。文中還發現用常用的B3LYP泛函算全局pi共軛的環狀聚合物的合理性較差，E/n隨n的變化明顯不合理。這些方面這里就不多提了，請讀者務必閱讀論文原文。

5 總結

環張力是環狀分子的重要特性。本文簡明扼要地介紹了兩種最常用的計算環狀分子的環張力能的方法，并以[6]CPP和18碳環作為具體例子進行了演示。同聯反應法比較方便省事，但對于一些全局pi共軛的環狀體系，諸如18碳環，則更適合外推法，明顯嚴格得多，同時外推法還很有助于探究環尺寸與環張力之間的關系。還應當注意并不是所有環狀體系都能用這兩種方法之一計算環張力能，比如苯分子，它存在顯著的和其環狀結構緊密相關的六中心鍵，因而環張力本來就是難以被定義的。

文/Sobereva@北京科音  2024-Jan-26

在《18碳環等電子體B6N6C6獨特的芳香性：揭示碳原子橋聯硼-氮對電子離域的關鍵影響》（http://www.shanxitv.org/696）中提到，筆者在Inorg. Chem., 62, 19986 (2023)一文考察B6C6N6分子的電荷分布時，專門計算了此分子的平面內π軌道、平面外π軌道、σ軌道和內核軌道上的電子是怎么分布在各個原子上的。原子空間定義的方法不唯一，此文用的是流行的Hirshfeld-I原子空間劃分，這種方式劃分的原子空間物理意義較強，可以較合理體現外環境導致的原子空間收縮和膨脹。這種分析方法對于讀者研究很多其它體系也很有益處。本文就演示一下怎么用Multiwfn實現這種計算，以計算B6C6N6的平面內π電子（π-in）的分布為例。由于這個計算需要充分利用Multiwfn的靈活性，牽扯一些細節，這是為什么我專門寫個文章來具體說明。

上述文章考察的B6N6C6的波函數文件是http://www.shanxitv.org/attach/697/B6C6N6_OS.rar，解壓后是B6C6N6_OS.fchk，是由Gaussian 16在wB97XD/def2-TZVP級別下以對稱破缺方式計算得到的。本文例子用的Multiwfn是2024-Jan-21更新的版本。Multiwfn可以在http://www.shanxitv.org/multiwfn免費下載，不了解者看《Multiwfn FAQ》（http://www.shanxitv.org/452）。

啟動Multiwfn，載入B6C6N6_OS.fchk。首先要做的是構造Hirshfeld-I原子空間，把Multiwfn的examples目錄下的atmrad子目錄挪到當前目錄下，這是因為此目錄下有各種元素不同價態的徑向電子密度信息，在構造Hirshfeld-I原子空間時要用到（若缺乏計算機常識不了解什么叫“當前目錄”，看http://www.shanxitv.org/237）。然后在Multiwfn里依次輸入

15  //模糊空間分析
-1   //選擇原子空間
1  //開始構造Hirshfeld-I原子空間
很快就構造完了。Hirshfeld-I原子空間權重數據現已被儲存在了內存里，輸入0返回主菜單。

因為我們要考察π-in分子軌道上的電子分布，故接下來需要把這類軌道以外的分子軌道的占據數都清零。可以按照《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）說的，在Multiwfn主功能0里一個一個觀看占據軌道的圖形判斷，把π-in軌道的序號都記錄下來，最終找到的序號如下。注意當前是對稱破缺計算，因此兩種自旋要分別考察。
Alpha自旋的π-in占據軌道：38,41,42,45,46,48,50,53,54
Beta自旋的π-in占據軌道：596,599,600,603,604,606,608,611,612（注意beta軌道序號在Multiwfn中的記錄規則，在http://www.shanxitv.org/269專門說了）

在Multiwfn里接著輸入
6   //修改波函數
26  //修改軌道占據數
0   //選擇所有軌道
0  //把所有軌道占據數清零
38,41,42,45,46,48,50,53,54,596,599,600,603,604,606,608,611,612   //所有π-in軌道序號
1  //占據數還原為原本的1.0
q   //返回
-1   //返回主菜單

現在就可以開始正式計算了。在Multiwfn里接著輸入
15  //模糊空間分析。進入后從選項-1的文字上可以看到當前的原子空間仍是Hirshfeld-I
1   //對特定實空間函數在各個原子空間中積分
1  //被積函數是電子密度。顯然當前對應的是π-in電子密度

馬上看到如下結果

  Atomic space        Value                % of sum            % of sum abs
    1(B )            0.64940465             3.607804             3.607804
    2(C )            1.00713759             5.595209             5.595209
    3(N )            1.34345778             7.463655             7.463655
    4(B )            0.64926116             3.607007             3.607007
    5(C )            1.00733834             5.596325             5.596325
    6(N )            1.34338216             7.463235             7.463235
    7(B )            0.64944181             3.608010             3.608010
    8(C )            1.00710182             5.595011             5.595011
    9(N )            1.34348431             7.463802             7.463802
   10(B )            0.64924542             3.606919             3.606919
   11(C )            1.00734233             5.596347             5.596347
   12(N )            1.34339152             7.463287             7.463287
   13(B )            0.64942036             3.607891             3.607891
   14(C )            1.00713369             5.595188             5.595188
   15(N )            1.34344824             7.463602             7.463602
   16(B )            0.64928283             3.607127             3.607127
   17(C )            1.00730639             5.596147             5.596147
   18(N )            1.34341808             7.463434             7.463434
Summing up above values:         17.99999847
Summing up absolute value of above values:         17.99999847

如上可見，B、C、N原子的π-in電子數分別是0.649、1.007、1.343，和如下所示的Inorg. Chem., 62, 19986 (2023)中的表2中的數據完全一致。

最后再提醒一下，必須按以上說明先產生Hirshfeld-I原子空間、修改軌道占據數，最后再在Hirshfeld-I原子空間里積分，而不能先修改軌道占據數然后再進入子功能15產生Hirshfeld-I原子空間并做積分。因為修改軌道占據數之后，此時的電子密度函數就不再是總電子密度了，而Hirshfeld-I原子空間在構造時用到的電子密度函數必須對應總電子密度。

文/Sobereva@北京科音   2024-Jan-26

0 前言

由18個碳原子連接形成的環狀體系18碳環（cyclo[18]carbon）具有十分特殊的幾何和電子結構，自從它于2019年首次在凝聚相中被觀測到后，18碳環及其衍生物就得到了化學界的廣泛關注，并陸續有大量理論研究工作發表。特別是在Carbon, 165, 468-475 (2020)中，北京科音自然科學研究中心（http://www.keinsci.com）的盧天和江蘇科技大學的劉澤玉等人證實18碳環具有明顯的雙pi芳香性特征，這是18碳環非常與眾不同的一點。

18碳環的一種等電子體是B9N9，由B和N原子交替相連構成環狀。此體系早已被研究過，被證實幾乎不具備像18碳環那樣的明顯的芳香性。最近，盧天、劉澤玉和西班牙多諾斯蒂亞國際物理中心的Mesías Orozco-Ic發現，由6個(硼-碳-氮)重復單元依次相連構成的18碳環的另一種等電子體B6C6N6具有和18碳環很接近的芳香性。為什么在無芳香性的B9N9中摻雜進去碳原子后，或者說讓碳原子橋聯硼和氮后，B9N9的芳香性就能得到巨大提升？這是非常有意思的問題。于是以上研究者對B6C6N6的特征做了全面深入的理論研究，詳細分析討論了其芳香性特征和本質，成果近期在美國化學會的知名期刊Inorganic Chemistry上發表，非常推薦閱讀：

Exploring the Aromaticity Differences of Isoelectronic Species of Cyclo[18]carbon (C18), B6C6N6, and B9N9: The Role of Carbon Atoms as Connecting Bridges, Inorg. Chem., 62, 19986 (2023) https://doi.org/10.1021/acs.inorgchem.3c02675

此研究的研究意義不僅在于討論了B6C6N6本身，還在于揭示了碳原子橋聯硼和氮原子對電子離域特征的關鍵性影響，這對于通過碳摻雜的方式對純硼、氮構成的體系（如硼-氮納米管、二維層狀h-BN體系等）進行改性提供了重要啟示。此研究通過強大的Multiwfn程序（http://www.shanxitv.org/multiwfn）利用了各種波函數分析手段對B6C6N6進行了充分的考察并和18碳環、B9N9進行了對比，這也是充分運用波函數分析探究新穎體系的很好的范例。

下面本文就對Inorg. Chem., 62, 19986 (2023)這篇文章的主要研究思想、結論進行深入淺出的介紹，并對一些研究細節進行附加說明，以便于讀者更好地理解此文，以及將此文的研究手段用于自己的研究。同作者之前還對18碳環及其衍生物的各方面特征做過充分的理論研究并得到了同行的廣泛關注，成果匯總見http://www.shanxitv.org/carbon_ring.html。

1 B6N6C6的幾何結構

做B6N6C6的各方面研究之前先得得到其可靠的幾何結構。18碳環及衍生物用ωB97XD/def2-TZVP級別做幾何優化是很穩妥的，在Carbon, 165, 468-475 (2020)以及http://www.shanxitv.org/carbon_ring.html里列舉的此類體系的大量研究工作中都已經充分證明了這點。因此此文用Gaussian 16在此級別下優化了B6N6C6的結構。特別需要注意的一點是，B6N6C6基態是單重態，但是其閉殼層單重態波函數存在不穩定性，用DFT計算時需要作為對稱破缺單重態來處理方可得到真正的基態，相關討論見《談談片段組合波函數與自旋極化單重態》（http://www.shanxitv.org/82）。因此，此文對B6N6C6基態的各種研究都是基于對稱破缺單重態的結構和波函數做的。如果誤當成了閉殼層單重態來算，作者發現得到的幾何結構將明顯不合理，而且電子結構也明顯錯誤，比如芳香性嚴重偏高。實際上之前有一篇論文也算過B6N6C6，但由于那些作者沒注意這一點而誤當成了閉殼層計算，導致結果是沒有任何意義的。此例體現出DFT計算B6N6C6這種電子結構不同尋常的體系一定要注意做波函數穩定性測試。而ωB97XD/def2-TZVP級別下的B9N9、18碳環的穩定波函數都是單重態閉殼層的。

ωB97XD/def2-TZVP級別下優化得到的無虛頻的B6N6C6、B9N9和18碳環的結構如下所示，紅/灰/藍分別對應硼/碳/氮，笛卡爾坐標在文章的SI里提供了。可見B6N6C6屬于C6h點群，而B9N9屬于D9h點群。

畢竟B6N6C6是一個新穎的結構，不排除靜態相關很顯著導致ωB97XD描述不理想的可能，故為了絕對穩妥，此文還用ORCA在CASSCF(12,12)/def2-TZVP級別下也做了幾何優化。此活性空間設定下，最終會有12個占據數明顯偏離整數的pi型自然軌道，而若活性空間進一步擴大到CAS(14,14)，則會再多出兩個占據數接近整數的pi型自然軌道，這體現CAS(12,12)恰好足夠充分描述B6N6C6的pi電子的靜態相關效應。CASSCF(12,12)優化出的結構和ωB97XD的非常相近，體現出ωB97XD得到的此體系的幾何結構靠譜。為了進一步體現ωB97XD計算出的波函數也合理，文中也順帶對比了一下這兩個級別得到的單電子分布情況。文章按照《談談自旋密度、自旋布居以及在Multiwfn中的繪制和計算》（http://www.shanxitv.org/353）計算了ωB97XD波函數的自旋密度，而由于CASSCF沒法得到自旋密度，故文章使用《使用Multiwfn計算odd electron density考察激發態單電子分布》（http://www.shanxitv.org/583）介紹的方法對CASSCF(12,12)波函數產生了單電子密度（ODD），等值面圖對比如下所示。在不考慮自旋密度符號（體現在等值面顏色上，不同顏色對應單電子不同自旋方向）的情況下，可見二者展現出的單電子分布是基本吻合的，即單電子主要都分布在碳上。這在很大程度上體現出ωB97XD對稱破缺計算產生的波函數也是合理的，適合用于之后的波函數分析。

2 B6N6C6的電荷分布

原子電荷是定量考察化學體系中電荷分布最常用的指標之一。18碳環由于所有的碳都是等價的，顯然原子電荷都為0，而對B6N6C6和B9N9計算原子電荷，可以直接考察這倆體系中不同原子帶的凈電量、了解二者電荷分布的差異。原子電荷計算方法并不唯一，不同方法基于不同的思想、有不同的特點，在《原子電荷計算方法的對比》（DOI: 10.3866/PKU.WHXB2012281）中有大量介紹和對比。考慮到B6N6C6電子結構的特殊性，為了得到穩妥的結論，文中用Multiwfn計算了很常用的ADCH、Hirshfeld、Hirshfeld-I、CHELPG原子電荷（它們的計算例子看Multiwfn手冊4.7節）、用NBO程序計算了常見的NPA原子電荷，結果如下所示。可見雖然不同原子電荷給出的具體數值不同，但可以得到共同的結論，就是氮明顯帶負電，硼明顯帶正電，這體現出了它們顯著的電負性差異。而在B6N6C6中，硼和氮的凈電荷差異明顯小于它們在B9N9中，這體現出碳元素在硼、氮之間的插入明顯促進了體系整體電荷分布的均衡性。而B6N6C6中碳自己的原子電荷則接近于0，即沒怎么受周圍化學環境的影響。

18碳環具有四類分子軌道：平面內和平面外兩套π軌道，分別稱為π-in和π-out，此外還有σ軌道，以及碳的內核原子軌道組成的內核分子軌道。B6C6N6和B9N9也有這四類軌道。為了深入考察各個原子上的電子是怎么分布在這些軌道上的，文中用Multiwfn基于Hirshfeld-I原子空間劃分對三個體系計算了電子在這些軌道上的布居數，如下所示。具體計算方法我專門寫了篇文章說明：《使用Multiwfn基于Hirshfeld-I劃分計算特定類型電子在各個原子上的分布量》（http://www.shanxitv.org/697）。

從以上數據可以看到，在每個原子上，π-in和π-out電子分布量都是基本一致的，沒有什么偏向性，體現出這些體系里π-in和π-out電子特征的近似等價性。B6C6N6和18碳環中的碳的電子結構特征非常相近，π-in和π-out電子數在兩種體系里都基本為1.0，明顯都處于sp雜化的狀態。

3 B6N6C6的成鍵情況

鍵級是考察化學鍵特征非常常見的方式，參見《Multiwfn支持的分析化學鍵的方法一覽》（http://www.shanxitv.org/471）。文中用Multiwfn對B6C6N6、B9N9和18碳環計算了流行的Mayer鍵級，并且分解為了不同類型軌道的貢獻，如下所示。Mayer鍵級體現了原子間等效的共享電子對數，由數據可見，三個體系里所有鍵都包含典型的單個σ鍵特征（即σ+core電子貢獻的Mayer鍵級接近非常接近于1），同時π-in和π-out電子都構成了一定程度的π鍵特征，且二者貢獻相近。B6C6N6中的π作用程度是N-B ≈ C-N > B-C。B9N9中所有B-N等價。18碳環由長、短兩種C-C鍵構成，后者的π作用顯著強于前者。文中還用Multiwfn計算了模糊鍵級以進一步確認了這些結論，數據見補充材料。

《使用IRI方法圖形化考察化學體系中的化學鍵和弱相互作用》（http://www.shanxitv.org/598）中介紹了盧天提出的IRI-π函數，可以直觀地展現化學體系中π電子相互作用情況。此文對B6C6N6和B9N9繪制了IRI-π=1.2的等值面圖并用sign(λ2)ρ函數著色，如下所示。可見每個鍵上都出現了環狀的等值面，非常直觀地體現出這些鍵都存在明顯的π作用，這和Mayer鍵級展現出的信息相吻合。

4 B6N6C6的芳香性

B6N6C6的芳香性是當前研究的重點。和18碳環一樣，B6N6C6和B9N9的π-in和π-out電子皆為18個，都滿足Hückel的4n+2芳香性規則，因此都有可能存在雙π芳香性。衡量芳香性的方法非常多，在《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）里有充分的介紹。其中基于電子結構和基于磁性質的芳香性衡量方法的可靠性和普適性普遍較好，可以用于考察B6N6C6的芳香性并與18碳環和B9N9的進行對比。

AV1245*1000在《使用Multiwfn計算AV1245指數研究大環的芳香性》（http://www.shanxitv.org/519）里專門介紹過，是基于波函數定義的適合用于定量考察共軛大環體系芳香性的指數。文中首先使用Multiwfn計算了AV1245*1000，結果如下。除了總值外，還分別給出了π-out電子和其它電子的貢獻（注：π-in、σ、core電子對AV1245的貢獻無法精確分離故沒有單獨給出。但顯然σ和core的離域性極低，不會對AV1245有什么貢獻，故“其它電子貢獻”實質上等于π-in電子的貢獻）。由數據可見，B6C6N6的芳香性顯著強于B9N9，并且π-out的芳香性比π-in還略強一點，這可能在于體系的π-in電子的共軛程度相對略低，源自于環的曲率使得平行于環的p原子軌道彼此交疊程度低于垂直于環平面的p原子軌道。單從AV1245*1000的數值來看，B6C6N6的芳香性甚至比18碳環還要更強一點。

NICS是非常流行和穩健的基于磁屬性衡量芳香性的指標，它的不同形式在http://www.shanxitv.org/176里有專門介紹。其中NICS(0)zz對本文研究的這些環狀共軛體系最為適合，結果如下所示，數值越負說明芳香性越強。為了更充分了解芳香性的來源，文中還按照《基于Gaussian的NMR=CSGT任務得到各個軌道對NICS貢獻的方法》（http://www.shanxitv.org/670）介紹的方法分別計算了π-in、π-out和σ+core電子的貢獻。NICS(0)zz體現出B6C6N6的芳香性接近18碳環，并且π-in和π-out的貢獻相仿佛。而B9N9雖然形式上也滿足4n+2規則，但由于其NICS(0)zz基本為0，因此可斷定是非芳香性的，這必然是由于其π電子缺乏整體離域能力所致。雖然B9N9的各個B-N鍵如前述Mayer鍵級和IRI-π所示有明顯π作用，但顯然不意味著其π電子能夠容易地離域。對這些體系，σ和core電子都由于其高度定域性而對芳香性基本沒有影響。

《使用AICD程序研究電子離域性和磁感應電流密度》（http://www.shanxitv.org/147）和《使用AICD 2.0繪制磁感應電流圖》（http://www.shanxitv.org/294）介紹的AICD圖是使用得較普遍的圖形化直觀展現化學體系磁感生電流的方法。為了進一步從磁感生電流的角度展現B6C6N6的芳香性以及和B9N9的差異，文中對這些體系的π-in、π-out以及全部π電子分別繪制了AICD圖，如下所示。外磁場方向垂直于環由下朝上施加。由于等值面上的箭頭太小，為了看得清楚，圖中還用大的橙色箭頭做了清晰的標注。此圖直觀體現出，B6C6N6確實存在顯著的π芳香性，因為感生電流方向符合左手定則并且連貫地環繞整個體系，這正是典型的π芳香性體系具備的特征。而缺乏芳香性的B9N9則明顯不具備這樣的特征，感生電流幾乎只是繞著原子核或原子間區域轉。

此文還按照《考察分子磁感生電流的程序GIMIC 2.0的使用（含24分鐘演示視頻）》（http://www.shanxitv.org/491）介紹的方法繪制了GIMIC磁感生電流的動畫，如下所示，更進一步地體現出芳香性明顯的B6C6N6能夠形成整體感生環電流而非芳香性的B9N9不具備這樣的特征。
B6C6N6：http://www.shanxitv.org/attach/696/B6C6N6_GIMIC.mp4
B9N9：http://www.shanxitv.org/attach/696/B9N9_GIMIC.mp4

文中還對穿越B6C6N6、B9N9和18碳環的化學鍵的感生電流進行了積分，數值分別為21.64、0.70、22.50 nA/T，這定量體現出B6C6N6和18碳環的芳香性相仿佛，而B9N9完全不具備芳香性。

在《通過Multiwfn繪制等化學屏蔽表面(ICSS)研究芳香性》（http://www.shanxitv.org/216）和《使用Multiwfn巨方便地繪制二維NICS平面圖考察芳香性》（http://www.shanxitv.org/682）中介紹了ICSS函數，它直觀地體現了化學體系的磁屏蔽和去屏蔽區域的位置和強度，是NICS的重要擴展，已被非常廣泛應用于討論芳香性問題。此研究使用Multiwfn對B6C6N6和B9N9都繪制了ICSS_ZZ = 10 ppm的等值面圖以及垂直于環的截面的填色等值線圖，如下所示。等值面圖中橙色和藍色分別對應屏蔽和去屏蔽區域。由圖可見，B6C6N6的環內部完全是屏蔽區域，環外面完全是去屏蔽區域，這種特征和苯這樣的典型芳香性體系完全相同，更進一步嚴格地證明了B6C6N6有顯著的芳香性。而且B6C6N6的ICSS的這種特征和Carbon, 165, 468-475 (2020)報道的18碳環的高度一致，體現出二者芳香性的高度共性。而B9N9的ICSS等值面圖的屏蔽和去屏蔽區域彼此交錯，沒有形成遍及整體的連貫的等值面，體現出此體系既沒有芳香性也沒有反芳香性，而是非芳香性。

5 B6N6C6的電子離域性

芳香性的本質來自于電子的離域，因此文中還特意對B6N6C6的電子離域性從多方面進行了考察。在《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）中介紹了Multiwfn可以計算和繪制的LOL-π函數，這對于考察π電子的離域情況極其有用，如今已經非常流行。下圖對比了B6N6C6和B9N9的LOL-π=0.5等值面，明顯可見B6N6C6的等值面比B9N9顯得連貫得多，明顯體現出B6N6C6的π電子有強得多的在18個原子上全局離域的能力，也一定程度解釋了前者的芳香性比后者強得多的原因。

LOL-π函數比較適合直觀定性考察，而《在Multiwfn中單獨考察pi電子結構特征》里介紹的ELF-π函數的二分點數值則適合定量分析對比。下圖是B6N6C6和B9N9的ELF-π函數等值面圖，等值數值設為了令等值面剛好分裂時的數值，這被稱為二分值，數值越大說明電子越容易跨越二分點區域發生離域。在圖上也標注了二分點位置和二分值。由圖可見，B6N6C6的π-in和π-out的二分值都明顯大于B9N9的，這給B6N6C6的整體電子離域性更強又進一步提供了定量證據。

電子的費米穴是考察電子離域特征的一個重要的函數，它包含兩個三維位置坐標，作圖考察時通常將一個位置作為參考點，然后對另一個坐標繪圖。Atoms-in-molecules理論的提出者Bader指出an electron can go where its hole goes and, if the Fermi hole is localized, then so is the electron，也就是假設一個電子出現在某個參考點位置，那么其費米穴函數明顯分布在哪，這個電子就容易離域到哪去。因此，通過繪制費米穴圖，可以更進一步了解體系處于不同位置的電子的離域能力。文中通過Multiwfn對B6N6C6、B9N9和18碳環分別繪制了費米穴的二維圖，如下所示。參考點位置以綠色箭頭標注，設在了不同原子的價層區域距離原子核特定距離的位置，此時圖中的費米穴的分布就體現了這個位置的價電子容易離域到哪去。用Multiwfn繪制這種圖的操作在《制作動畫分析電子結構特征》（http://www.shanxitv.org/86）里有相關說明。

由上圖等值線分布可見，B6N6C6的硼和碳的價電子的費米穴分布范圍較廣，體現出電子有往較遠區域離域的能力，特別是其碳的費米穴分布特征和已知具有較顯著全局離域性的18碳環的很接近。相比之下，B9N9的硼和氮的費米穴分布廣度明顯不及B6N6C6的相應元素的，等值線只分布在距離參考點較近的原子上，體現出B9N9的價電子的離域能力相對較差。

上面的等值線圖適合定性討論，為了能把以上信息轉化成定量形式來更好地對比分析，文中定義了一個新的量叫做原子遠程離域指數（atomic remote delocalization index, ARDI）。離域性指數（DI）在《Multiwfn支持的分析化學鍵的方法一覽》（http://www.shanxitv.org/471）里專門介紹過，是衡量特定兩個原子間等效共享電子對數用的，而ARDI則是對單個原子定義的，定義為距離它兩個化學鍵以上的原子與當前原子的所有DI之和，顯然基于Multiwfn計算的所有原子間的DI可以簡單地手算出來各個原子的ARDI。ARDI體現了特定原子上的電子的遠程離域能力，數值越大說明遠程離域能力越強。環上各個原子上的電子有顯著的遠程離域能力顯然是這個環具備明顯芳香性的前提。各原子ARDI如下所示，可見B6N6C6上的原子的ARDI都較大，平均值與18碳環的碳的相仿佛，而B9N9的硼和氮的ARDI都遠小于B6N6C6相應元素的原子，這更充分體現了B9N9的電子的全局電子離域能力遠不如B6N6C6，進一步展現了B9N9芳香性遠低于B6N6C6的原因。

6 碳摻雜進B9N9使得芳香性巨幅提升的本質

根據前面的討論，文章已經充分證明了芳香性的關系是B6C6N6 ≈ C18 ? B9N9。為什么將碳元素摻雜進入無芳香性的B9N9，就能使其芳香性巨幅提升、電子全局離域能力大幅增加？實際上原因并不難理解。盡管前面的Mayer鍵級體現出B9N9中的硼-氮之間有一定π作用，但由于硼和氮的電負性差異非常大，這無疑會導致電子傾向于在氮上定域而難以往相鄰的硼上離域，至于往更遠的原子上離域就更是難上加難了。而碳的電負性則正好介于硼和氮之間，顯然硼-碳和氮-碳鍵的共價性顯著強于硼-氮鍵，因此電子更容易離域過去，相當于碳給硼和氮之間架設了能夠令電子容易離域過去的橋梁，充分調和了硼和氮之間電負性差異過大導致的電子跨鍵離域太難的矛盾。

為了更好地展現碳原子起到的價值，文中將垂直于環的各原子的2p原子軌道繪制了出來，π-out電子正是分布在這些p軌道上，π-out分子軌道也正是它們線性組合產生的。具體來說，由于NBO程序產生的預自然原子軌道（PNAO）可以較合理地描述分子環境中的原子軌道，因此文中按照《使用Multiwfn繪制NBO及相關軌道》（http://www.shanxitv.org/134）中的做法，用Multiwfn結合VMD將相應的PNAO軌道用等值面圖形式做了展示，如下所示。等值面數值選為了適合對比的0.1。圖中也把PNAO的能量（eV）標注在了上面。由圖可清楚地看出，B9N9中的氮和硼的p軌道的空間延展程度差異很大，而且軌道能量差異也很大，勢必電子難以在它們之間離域。而從B6C6N6的圖中可見，碳的p軌道延展范圍和軌道能量都恰好介于硼和氮之間，它的引入無疑能夠顯著幫助電子離域過去。

7 B6N6C6異構體的芳香性

以上研究的那種以B-C-N作為重復單元排列的B6N6C6結構是否真的是B6N6C6的能量最低結構？是否B6N6C6化學組成的體系還有芳香性更強的異構體？這是個值得考察的有趣的問題。為此，文中還另外考慮了兩種B6N6C6構型，優化后得到的無虛頻結構如下所示，可見B6N6C6'構型讓所有碳盡可能連在一起，其它部分部分保持B9N9那樣硼-氮交錯的結構，而B6N6C6''構型則是由三個C-C-B-N-B-N重復單元構成的結構。這兩個結構的基態是閉殼層的，能量分別比前面研究的B-C-N為重復單元構成的B6N6C6結構能量低280.0和250.2 kcal/mol。這體現出B6N6C6中B-C和C-N的平均鍵能比B-N和C-C是要更小的。雖然前面研究的B6N6C6結構能量不是這種化學組成的能量最低結構，但并不意味著它是不穩定的。根據《使用ORCA做從頭算動力學(AIMD)的簡單例子》（http://www.shanxitv.org/576）中介紹的方法，文中對此體系做了較高溫度（500 K）下20 ps長度的基于wB97X-D3/6-311G*級別的從頭算分子動力學模擬，用VMD每1 ps繪制一次的結構疊加圖如下所示（顏色按照藍-白-紅變化），可見環狀結構始終很好維持著，并沒有發生異構化、解離等現象，體現出B6N6C6具有一定熱力學穩定性。

文中還對B6N6C6'和B6N6C6''構型計算了NICS(0)zz，發現分別僅有1.8和1.4 ppm，可見這倆構型沒有芳香性。這說明只有讓碳插入到硼和氮中間，從而以B-C-N作為重復單元時，碳才能真正起到顯著提升純硼-氮環體系的電子離域性的作用。

8 B6N6C6最低三重態激發態的芳香性

前面考察的都是B6N6C6的單重態基態（S0）。最后，文章還研究了B6N6C6最低三重態激發態（T1態）的芳香性。計算發現S0-T1垂直和絕熱激發能分別為1.62和1.17 eV。下圖(a)是優化后的T1態B6N6C6的結構，可見對稱性相對于S0態發生了下降。下圖(b)是按照《使用Multiwfn作電子密度差圖》（http://www.shanxitv.org/113）繪制的在B6N6C6的S0結構下T1態與S0態的密度差，藍色和綠色分別為電子密度降低和增加部分。由圖可見S0到T1態的電子激發是π-in → π-out的激發。由于垂直激發導致的電子密度變化破壞了B6N6C6基態的C6h對稱性，這也是為什么在T1態勢能面上從Franck-Condon點往極小點弛豫過程中幾何結構對稱性發生了下降。

文中對B6N6C6的T1態計算了NICS(0)zz，并分解為了不同類型電子的貢獻。而且為了考察幾何結構弛豫的影響，在S0極小點和T1極小點分別做了這種計算，結果如下所示。可見無論在哪種結構下，T1態由于NICS(0)zz遠大于0，因此是顯著反芳香性的，這和Baird規則所述情況一致。Baird規則指出單重態基態滿足Hückel 4n+2規則的體系的最低三重態是反芳香性的。從NICS(0)zz還看到，T1態的反芳香性特征在T1極小點結構下比S0極小點結構下更弱，這體現出B6N6C6自發減小T1態的反芳香性可能是其S0→T1垂直激發后結構弛豫的重要驅動力。

9 總結

此文基于量子化學計算和Multiwfn等程序做波函數分析，非常全面系統研究了18碳環的等電子體B6N6C6的幾何結構和電子結構，從不同角度全面考察了其芳香性和電子離域性特征，并且和18碳環及其另外一個重要的等電子體B9N9進行了細致的對比。本文體現出以碳原子橋聯硼和氮原子，可以顯著提升純硼-氮體系的電子離域性。這不僅對電子結構產生重要影響，必然也同時會影響其它性質，諸如光學性質、電子輸運性質等。本文雖然研究的只是一維體系，可以預想電負性介于硼和氮之間的碳元素的摻雜必定也會對二維、三維純硼-氮體系產生顯著影響，而且摻雜方式的選取可能會對這些體系的性質起到一定調控作用。

如果讀者對于18碳環及其相關體系感興趣，非常推薦查看http://www.shanxitv.org/carbon_ring.html里匯總的本文作者的巨量其它研究工作，包含大量深入淺出的論文內容和研究思想介紹文章，以及對計算技術細節的附加說明。

此文對筆者完整看過的2023年度完結的79部全長版TV動畫和7部泡面番進行總評，約20萬字。劇場版、OVA之類的都沒有納入其中，因為另有文章。此文只有“2024年1月完結”部分是新寫的，2023年4月、7月、10月完結的動畫的評論在作品剛完結之時就已經寫好并發在了筆者的blog上，在此文只是匯總而已（所以有些文字現在看起來已經out了），不過很多作品的總評分為了保持評價的一致性，在此文又做了微調。以前寫的總評見（不同年份的評分標準有一些系統性不同，不適合橫向對比）：

2022年度完結動畫總評：http://www.shanxitv.org/669
2021年度完結動畫總評：http://www.shanxitv.org/633
2020年度完結動畫總評：待發布
2019年度完結動畫總評：http://www.shanxitv.org/544
2018年度完結動畫總評：http://www.shanxitv.org/470
2017年度完結動畫總評：http://www.shanxitv.org/401
2016年度完結動畫總評：http://www.shanxitv.org/358
2015年度完結動畫總評：http://www.shanxitv.org/316
2014年度完結動畫總評：http://www.shanxitv.org/273
2013年度完結動畫總評：http://www.shanxitv.org/219
2012年度完結動畫總評：http://www.shanxitv.org/174
2011年度完結動畫總評：http://www.shanxitv.org/118

2023年筆者觀看的非TV版動畫另有文章進行評論：
2022下半年和2023上半年觀看的非TV版動畫簡評（http://www.shanxitv.org/673）
2023下半年觀看的非TV版動畫簡評（http://www.shanxitv.org/694）
 
筆者正在追的2023年10月開播的《葬送的芙莉蓮》和《藥屋少女的呢喃》尚未完結，所以不納入本文。

筆者是堅定的掃番黨，所有新番無論喜好至少看第一話，不少都觀望到第3話才決定是否棄。沒有納入此文的動畫，按筆者的觀點來說大部分就屬于無聊、跟風、庸俗、粗制濫造之作（可能讀者會發現在視頻網站上播放數高的作品，反倒出現在本文的幾率低，因為主流的品味實在是不敢恭維，而真正的好作品卻總是受到冷落）。也有極少數作品雖然可能做得不錯卻沒出現在此文，是因為那些題材著實不是灑家的菜，諸如恐怖惡心之類，所以也就沒看下去。

評分主體區間是90~95分（低于90分根本不追，超過95的則是超神作水準，極度罕見），不了解此系列文章評分傳統的勿試圖吐槽這點。評分依據是作品畫面、音樂、表現力、創意、劇情、人設等因素，以及很少部分個人主觀觀感和興趣因素。

注意，本文的評論有的包含一些劇透成份，會影響觀賞樂趣。如果目的是想找一些好看的，只要看頒獎部分就夠了，等看完了再看評論為宜。

每個獎項里面的順序是隨意的，不分先后。

題外話：我只看日本動畫，因此此文里看不到其它國家的動畫。我總看到有一些人無視每年層出不窮的大量杰出的日本動畫作品，卻往往拿撐死了也就占70%比例的那些很俗套的商業化日本動畫作品（這些我也基本都不追）為由貶低和唱衰日本近年來的動畫，對于他們我表示極度遺憾，錯過了太多佳作，而且他們的言論也嚴重偏激和誤導他人，令我十分反感。若做不到像我一樣完整掃番（就算哪怕只掃1/3也行）、認真觀賞，就根本沒有做評價的資格。上述這樣的人也有些是因為自己年齡和心態發生變化導致觀賞體驗發生變化，卻歸咎于是當下日本動畫本身的問題，這顯然也是很不負責任的。還有些人對日本動畫主觀地懷有惡意，甚至會因為一些動畫里有一些其不喜歡的元素（如賣萌賣肉）就無腦地批判和否定整部作品乃至業界。比如有人說“現在日本動畫千篇一律，都是賣萌賣肉、穿越異世界，沒意思”，我就想反問這些人，你完整看過《石紀元》？完整看過《16bit的感動》？看過《藥屋少女的呢喃》？這算賣萌賣肉？這是穿越異世界？它們的出色，明明靠的是杰出的創意、劇情、人設以及畫面/音樂/配音等各方面高水準的制作質量。倘若你看到了自己不喜歡類型的作品，別張口就隨便否定業界，近年來的日本動畫里有的是其它類型的非常優秀的作品可推薦的。我從2008年開始認真掃番至今，每年完整追100部左右TV版動畫（另加許多劇場版、OVA），可以非常準確客觀地在此表示，日本動畫整體品質的平均值從未有任何下降，而且出色作品的數量也根本沒有減少過，依然是百花齊放，繁榮昌盛的狀態。我想強調的是，做出評論一定要負責，諸如我雖然不看其它國家的動畫，并且對某些國家的動畫有一些長久以來的負面印象和看法，但由于我沒認真全面觀賞過其近年的作品，無法確定是否有一些優秀作品是我不知道的、乃至我看了后說不定會追的，所以我就不會輕易做出評價。我一直致力于寫動畫評論，主要目的之一是為了讓真正優秀的日本動畫得到更多的關注，給觀賞時間有限的（而且觀賞品味不離譜的）人提供指引，讓他們更容易接觸到出色且適合他們觀賞口味的作品，我認為這是很有意義的。

文/Sobereva  2024-Jan-2

這里對我從2023年7月1日開始截止到2023年12月31日觀看的非TV版動畫進行簡單評論。由于這類動畫不好對比和打分，我也就不給出評分了。以往的此系列文章和TV版動畫評論見http://www.shanxitv.org右方的“動畫評論”分類。

黃金拼圖 Thank you!!
類型：劇場版
長度：1h21m
我作為TV版老觀眾對這個劇場版很滿意，內容很有趣，觀感很好。這個劇場版主要講述高三畢業前夕大家一起經歷的故事，包括休學旅游、日常、參加大學考試等劇情。這個劇場版百合濃度很高，有很多高橘的畫面和對話，能看得人挺開心。

令我再次感嘆可憐的人設真出彩，顏值真是相當高，賞心悅目，是此作里我最喜歡的角色。居然可憐留在日本上大學，而忍和愛麗絲一起去英國上大學了，從此可憐和愛麗絲這青梅竹馬終于分開了。

穗乃花真是喜歡可憐，老是想著在修學旅行過程中和可憐單獨拍照，但又不好意思，她還挺漂亮的。可憐的老相好愛麗絲如今和她不在一起，可憐又很友善又很豁達和open，如果穗乃花多主動一些引起可憐的注意力，我感覺穗乃花在大學期間還是挺有希望和她愛慕的可憐發展出親密關系的。

豬熊陽子和小路綾之間繼續橘里橘氣，在這劇場版里綾對陽子表現出各種心動和臉紅。這倆人不僅上了同一所大學，而且之后還租了房子一起同居了，這發展挺不得了，大橘已定！綾是那種看起來學習優秀的孩子，而躁動的陽子的學習成績應該不如她。陽子能和綾考上同一所大學，不知道是不是綾為了能和陽子在一起而遷就她選了好考一些的大學。

忍繼續表現出對愛麗絲的各種癡迷。想不到忍的媽媽當年和愛麗絲的媽媽也有過美好的回憶，可惜后來都和男性結婚了，百合緣分就只能留給她們的女兒繼承了。愛麗絲她媽真是挺漂亮，感覺比愛麗絲更有風采，要是在此作中多出鏡就好了。

看到最后，綾、陽子和可憐去英國看望在那里讀書的忍和愛麗絲的時候，我再次心想，11區的護照就是方便...

總之就是非常可愛 女子高中篇 02話
類型：OVA
長度：23m40s
這個OVA主要講述星空在女子高中當老師時遇到的女學生的事。內容包括司和星空去店里吃拉面、司換上女高中生校服、好幾個女學生來澡堂、星空自己動手在院子里制造水塘。這個OVA很有趣。

那個家境富有的女高中生的設定像極了《愛吃拉面的小泉同學》，走遍各個拉面店并在網上寫評論，感覺作者肯定是借鑒了小泉的人設。

吃拉面的時候星空還挺貼心，看到司點了個巨辣的拉面時，刻意點了個不辣的拉面，免得到時候可愛的妻子被辣到時沒的吃，這令和他們同桌的那個女學生被閃到。

星空的好幾個女學生看上去不漂亮，但到了澡堂里把頭發放下來后一下子顏值增加了許多。這些女學生還蠻有趣的，要是能再在續作中看到就好了。這些女生對星空還挺感興趣的，不過并沒有做出倒貼的行為，也知道他已經有了超級可愛的妻子。但我還是挺想在此作里看到這些女生挑逗星空令他臉紅的畫面的。雖然星空和這些女學生之間什么都沒發生，但還是令司感到吃醋，挺有趣的，她對星空的占有欲還真強。

星空居然憑自己的能力就把水塘造出來了，弄得還有模有樣，就像日本庭院里那種專門建造的水池，真是不簡單，腦力和動手能力真是相當強。

星空和穿了高中制服的司在屋后陰影處除了接吻外都做了什么，令人浮想聯翩，就像男女高中生在學校的角落做不可告人的事一樣。

碧藍航線 女王號令 OVA 01、02
長度：2*24m
我對《碧藍航線》沒有絲毫興趣，一丁點都沒有，興趣是負值。只是偶然看到出了這個OVA就隨便看了看。發現此作簡直就是世界杯，各種壯觀的又大又白。這些人設也太媚宅了。巨r角色的比例實在太高了。

OVA01劇情沒什么營養，內容就是女王伊麗莎白下象棋輸給了英勇，然后英勇當了一天女王，而傻乎乎的伊麗莎白去當了一天下人，算是體驗生活了。劇情沒什么意思，完全就是觀球了。

OVA02劇情也沒啥營養，內容就是要搞陣營聯歡會，伊麗莎白去各個陣營查看準備情況的事。這一話出現了俾斯麥，人設明明就是抄襲艦娘的俾斯麥嘛！

虹咲學園學園偶像同好會 NEXT SKY
類型：OVA
長度：29m40s
之前在TV版最后看到步夢在英國和兩個妹子在一起，令我有點好奇。這個OVA講述其中一個褐皮妹子艾拉來日本體驗偶像文化，大家一起陪她玩的事。感覺這OVA也沒講什么內容。英國的另一個白發的妹子當初令我眼前一亮，可惜在這個OVA里只是打了醬油，沒有什么描寫，感覺她比艾拉的人設更有亮點，希望在以后的作品里她能有一些正經的戲份。可能艾拉以后和她會在英國搞學園偶像活動了。艾拉給我的感覺一般般，個性不怎么鮮明。這OVA里鹽子的戲份不少，之前她一直很正經，甚至有點死板的感覺，結果在這OVA最后終于也決定放飛自我了，是個重要轉變。眾人和艾拉一起逛秋葉原的時候看到了UTX大樓，這一幕令我長嘆，當年UTX的超級學園偶像二賴子再也看不到了，A-Rise組合真的是非常出彩，當年被ll官方浪費掉了太太太太太可惜了！！！這OVA開頭的說唱部分真難聽。居然這OVA有兩個ED。

蠟筆小新 2023年劇場版
這是我于2023年9月初在大阪陪別人看的。其實我一直以來看不上《蠟筆小新》，并且感覺低幼+低俗。但這個劇場版令我感覺挺出乎意料的，做得著實很不錯，歡樂搞笑、驚險刺激、感人催淚的元素全都有，而且這些元素融合得非常理想，毫無生硬和突兀感。整部作品畫面制作和劇情構思都很下功夫，真是一部老少咸宜的用心的佳作。此作的畫面是三渲二，一開始感覺有點奇怪，看過之后感覺此作也只能用這種方式制作，畢竟此作里面大場面相當多，手繪難度太大。而且三渲二效果還挺自然，令觀感并未打什么折扣。

此作主要內容是有兩個隕石撞到了地球上，一個撞到了小新和妹妹身上，另一個撞到了一個性格陰暗、失去理想的男青年身上。那個男青年的人生一直非常灰暗，特別缺愛，作為心靈支柱的他鐘愛的偶像妹子又突然宣布和異性交往了，導致男青年在得到超能力后決定報復社會。小新在此作里利用自己的超能力和他對抗。先是挫敗了他綁架小新所在幼兒園的行動，之后又將變得巨大化的他擊敗。最后男青年失去神志，被超能力和陰暗心理所支配，變成了異形，小新被吸入他肚子里之后接觸到了他過往的記憶和靈魂。這一段內心對話構思得真是相當好，如果沒有這部分，此作就是簡單歡樂過個眼癮的作品，而有了這個部分就給此作內涵和深度拔高了很多，挺催淚的。這個部分很好地展現了外在顯得沙雕的小新的內在品質優良的一面，小新努力和心靈孤單、飽受欺凌的男青年做朋友，友情以及他人的支持也正是那個男青年所最需要的。最終心態非常積極的小新改變了陰暗男青年，并且在小新他爸臭襪子作為武器的情況下戰勝了男青年內心負面情感幻化出的三個反派。最后男青年回歸人形，也在大家的鼓勵下決心好好努力面對人生，可喜可賀。

看了這個劇場版，真是令我再次感到11區的非常多的動畫真是非常用心和良心，哪怕是面向子供的作品也真是做得毫不含糊，反觀某國的低齡動畫那品質真是連低齡兒童都不想看。此作令我對《蠟筆小新》的印象有很大的改觀。結尾看到小新的2024年劇場版的預告，是關于恐龍的，不知道是不是小新要穿越回古代，想來《機器貓》里也有類似的劇情。如果到時候我人在日本的話打算去看。

王者天下3 真人電影（Kingdom 3）
我頗為喜歡《王者天下》的動畫，之前也聽說過有真人電影，但一直沒機會看。2023年9月初我在大阪的TOHO CINEMA意外地看到有這個真人版電影，就果斷看了，票價折合100 RMB。

這個真人電影講的劇情我印象中在動畫版里大部分都有，但是是很多很多年以前的劇情了，應該是王者天下第二季那會兒，講王騎率兵跟趙國交戰的事。王騎讓當時還是百人將的信去偷襲敵人后方，信經歷殊死拼殺后成功斬了敵軍總大將。故事中還回憶了嬴政少時被紫夏帶著逃離趙國時，路途中被追殺的那段故事。

這電影版的前1/3以宮殿里的對話為主，感覺比較沉悶，特別是對于聽不懂日語的人來說。嬴政擺脫追殺，以及后面的戰爭劇情，都做得很不錯，很有看頭。

王者天下原作里的角色普遍都很有個性，也令我這位動畫版老觀眾對此作的選角頗有看法。

王騎演的只能算是及格。那種爆棚的自信，以及骨子里的那種奇妙的幽默和妖嬈感，都演得不錯，但關鍵是欠缺威武感和壓倒性的氣場，他剛走進宮殿時我一瞬間都認不出他來。明明動畫里王騎是個形象高大的武將，結果這真人版里他只是個一般身高的角色，總大將的氣場實在是展現不出來，甚至覺得他是個有點滑稽的角色。

騰的演員基本沒說話，外形并不像動畫版那么顯得像洋人，不過總體氣質還是挺接近騰的，令我比較滿意。

嬴政演得還算可以，但眼神和智慧感和動畫版有明顯差距，80分吧。信演得說不上多出彩，但也能令我滿意，85分。

河了貂居然是橋本環奈演的。之前我覺得她演河了貂大抵只是因為橋本環奈的人氣才被選上，但看了此作之后感覺她演得不錯。橋本環奈長得可愛，而且個子又很小只，演河了貂真是很合適，討人喜歡，換成其她演員也不會明顯更好。可以打93分。可惜河了貂在此電影里主要就是打打醬油，沒什么戲份。

蒙武演得很不錯，長相和動畫版如出一轍，很有野性和霸氣的感覺。美中不足的就是演員個子還不夠高大，要不然就可以接近滿分了。

此作演員選的最糟糕的是羌瘣，完全不及格。完全沒有動畫版的氣質，演員很顯老，明顯不像是靈秀的少女，表情也挺木訥的，一點神采都沒有。不知道當初選角怎么選的，懷疑是不是當初還可以，但過了N年拍第3部的時候其容貌已經發生了很大變化，但又不方便中途換演員。

楊端和的演員還可以，雖然顏值達不到原作中她的level，但整體感覺也還可以了，不令我失望。可惜她在這一部里總共出現了也就10秒鐘，完全是怒領工資（似乎第1或第2部里她有不少細分）

紫夏演的很出色，可以給97分，顏值不錯，特別是演技很好，充分展現出她對嬴政的重視和愛護，最后她為了保護嬴政而被射死的時候相當壯烈。

呂不韋演的很爛，完全沒有動畫版里的氣場，感覺這真人版里他就是個糟老頭子，看不出有任何過人之處。完全不及格。

肆氏一出場就看得出是肆氏，發型很獨特，發際線有幾個尖角。

信的部隊里有好幾個很有個性的小兵，看動畫版的時候我就挺注意到他們，他們在這真人版里被演得活靈活現，這批又不起眼又不可忽視的小角色的選角很不錯。馬桶蓋、大齙牙、大棒槌、老漢等小兵，在此作里全都演得很鮮活。

邪神與廚二病少女 世紀末篇
類型：OVA
長度：23m40s
這是小邪神的OVA。這個OVA的內容繼續無厘頭，又挺有趣的。講述百合鈴成了城保町的獨cai者，然后小邪神到各個地方連哄帶騙呼朋喚友找幫手去找百合鈴對決，美其名曰是推翻獨cai、把百合鈴說得像拉奧似的，實際上就是想找個由頭再次搞掉自己的宿敵。這個OVA的一個巨大槽點是拳四郎居然出現在了此作中，北斗神拳的主題曲的令人印象深刻的伴奏都出現了。小邪神雖然從他那里學了暗殺拳，但底子太弱，再次不知天高地厚地挑戰百合鈴時果然又被虐。更想不到居然百合鈴還從從他面前路過的拳四郎那里學了暗殺拳，這回小邪神被她這招虐得比以往還要慘，腸子紛飛，還像烤魚一樣被木桿貫穿身體叉在公共場所示眾，好慘。小邪神真是從不長記性，如果有下一季肯定還是要沒完沒了作死。也不知道到底以什么方式才能徹底搞死小邪神，有種把她剁碎了包餃子被眾人吃了都還能復活的感覺。這OVA登場角色真夠多的，幾乎所有TV版的角色都出現了，CV成本肯定很高。

文/Sobereva@北京科音   2023-Dec-30

0 前言

耦合簇方法可以通過線性響應（linear response, LR）計算激發態。在常見范疇內，精度和耗時關系是LR-CC3>LR-CCSDR(3)≈EOM-CCSD(T)>LR-CCSD=EOM-CCSD>=CC2。LR-CC3如今被普遍視為是激發態計算的金標準（多參考特征很強以及雙電子激發特征占主導的情況除外），激發能能精確到零點零幾eV。但LR-CC3結合像樣基組只能算得動個位數原子，對更大體系則可以用明顯更便宜的LR-CCSDR(3)，精度也已經非常好了。

Dalton程序在耦合簇方面功能豐富，可以做LR-CC3、LR-CCSDR(3)、LR-CCSD、LR-CC2激發能計算，對于LR-CC2/CCSD/CC3還可以得到包括振子強度在內的諸多躍遷屬性。Dalton的輸入文件對于初學者來說比較復雜抽象，本文用盡可能簡單易懂的語言完整介紹一下用免費的Dalton做這些激發態耦合簇計算的方法，使用有C2v對稱性的非常典型的分子甲醛為例。讀者如果對Dalton不熟悉，務必先閱讀我之前寫的《量子化學程序Dalton的編譯方法和運行方式簡介》（http://www.shanxitv.org/463），我假定讀者已經認真看過此文。

值得一提的是，雖然Dalton做LR-CC3計算如今在文獻中很常見，但Dalton在這方面不是最快的，而且代碼還沒有并行化。如果追求更好的效率，可以用專門做耦合簇計算的e^T程序（https://www.etprogram.org），它的LR-CC3是所有程序中明顯最快的。另外，LR-CCSD激發能和振子強度沒任何必要非得用Dalton計算，常用的Gaussian和ORCA的EOM-CCSD也都做得很好，結果和Dalton的LR-CCSD是等同的。Dalton做耦合簇的激發能計算不支持考慮溶劑模型，而G16的EOM-CCSD則可以。還值得一提的是，Dalton的衍生程序LSDalton還額外支持RI-CC2，可以比Dalton的LR-CC2更高效率地計算較大體系，這不屬于本文的范疇。

本文涉及的輸入輸出文件都可以在http://www.shanxitv.org/attach/693/file.zip里獲得。計算用的是2023-Dec-22通過git下載的Dalton最新的開發版。

本文下面的例子首先要重復知名的TBE1激發能測試集（J. Chem. Phys., 128, 134110 (2008)）里的LR-CC3/def-TZVP算的甲醛的1-1A2、1-1B1、2-1A1三個單重態的垂直激發能。如文章表1所示，結果分別為3.95、9.18、10.45 eV。此文用的幾何結構是MP2/6-31G*級別優化的，對應的結構文件是file文件包里的H2CO.xyz，本文的計算也將基于這個結構。

1 確定不可約表示順序

為了指定各個不可約表示的激發態各算幾個態，需要先知道C2v點群的各個不可約表示在Dalton程序中的順序。最簡單的方法是做個DFT單點計算，看一開始輸出的不可約表示的順序。為此，我們用Multiwfn（http://www.shanxitv.org/multiwfn）創建這個任務的輸入文件。啟動Multiwfn，然后依次輸入
H2CO.xyz
100   //其它功能（Part 1）
2   //導出各種文件
19   //Dalton
DFT.dal  //在當前目錄下產生DFT.dal，對應B3LYP/6-31G*計算設置
H2CO.mol   //在當前目錄下產生Dalton格式的體系定義文件H2CO.mol

由于當前計算要考慮對稱性，因此把H2CO.mol里的Nosymmetry刪掉。然后基于DFT.dal和H2CO.mol做計算，從輸出文件DFT_H2CO.out里可以看到在SCF開始之前就輸出了當前自動判斷的點群，以及相應的各個不可約表示的順序

@    Full point group is: C(2v)         
@    Represented as:      C2v

@  * The irrep name for each symmetry:    1: A1     2: B1     3: B2     4: A2

2 做LR-CC3激發能計算

寫一個文本文件LRCC3.dal，內容如下（這是最簡單寫法，默認設置就已經適合的選項就沒再做額外設置）

**DALTON
.RUN WAVE FUNCTIONS
**WAVE FUNCTION
.CC
*CC INP
.CC3
*CCEXCI
.NCCEXCI
2 1 0 1
**END OF DALTON

對當前體系默認是做閉殼層HF計算得到參考態波函數。**DALTON控制任務類型，.RUN WAVE FUNCTIONS要求算單點并得到波函數。**WAVE FUNCTION設置波函數計算類型，.CC要求做耦合簇計算。*CC INP設置做什么形式耦合簇計算，.CC3要求做CC3計算。*CCEXCI要求做耦合簇的激發能計算，.NCCEXCI下面按照不可約表示的順序指定各個不可約表示的能量最低激發態各計算幾個。當前為了重復TBE1測試集的數據，要求A1算兩個（由此得到1-1A1和2-1A1），B1算1個（得到1-1B1），B2不算，A2算1個（得到1-1A2）。

然后把前述的H2CO.mol復制為H2CO_TZVP.mol，手動把里面的6-31G*替換為Turbomole-TZVP，這是Dalton內置的def-TZVP基組的寫法，對應Dalton目錄下的basis子目錄中的Turbomol-TZVP這個文件的名字。

然后基于LRCC3.dal和H2CO_TZVP.mol做計算，得到LRCC3_H2CO_TZVP.out。本節的文件都已經提供在了前述的file文件包里。打開輸出文件，從里面可以看到如下激發能信息，2-1A1為10.44842 eV，1-1B1為9.18343 eV，1-1A2為3.94711 eV，可見和TBE1原文里的10.45、9.18、3.95 eV完全一致。表中||T1||是單激發貢獻，和TBE1原文里給出的也是一致的。

 +=============================================================================+
 |  sym. | Exci.  |        CC3        Excitation energies            | ||T1||  |
 |(spin, |        +------------------------------------------------------------+
 | spat) |        |     Hartree    |       eV.      |     cm-1       |    %    |
 +=============================================================================+
 | ^1A1  |    1   |     0.3502215  |       9.53001  |     76864.734  |  91.00  |
 | ^1A1  |    2   |     0.3839723  |      10.44842  |     84272.170  |  91.32  |
 +-----------------------------------------------------------------------------+
 | ^1B1  |    1   |     0.3374848  |       9.18343  |     74069.360  |  90.93  |
 +-----------------------------------------------------------------------------+
 | ^1A2  |    1   |     0.1450537  |       3.94711  |     31835.610  |  91.16  |
 +=============================================================================+

如果想要計算三重態激發態，就在.NCCEXCI下面的第二行定義。比如下面的寫法，代表除了計算如上那些單重態激發態以外，還計算1-3A2和1-3A1三重態激發態
.NCCEXCI
2 1 0 1
1 0 0 1
如果不需要計算單重態激發態，就寫
.NCCEXCI
0 0 0 0
1 0 0 1

筆者試了幾個Dalton版本，包括2016、2018、2022開發版，算出來的三重態激發能都是錯的，和TBE1表2里的1-3A2和1-3A1三重態激發能明顯不符，而且||T1||嚴重偏低，我認為是bug。如果不利用對稱性，即.mol文件里帶著Nosymmetry，并且.NCCEXCI下面直接指定計算的激發態的總態數，則問題輕得多，可以看到||T1||都顯著高于90%。

值得一提的是LR-CC3是對激發態一個一個獨立計算的，故算的態數越多耗時明顯越高。

如果只需要計算CC3基態能量，去掉*CCEXCI及下面的部分即可，相對于LR-CC3算激發態部分的耗時來說可以忽略不計。

3 做LR-CC3激發能+振子強度的計算

如果對上面算的激發態還要計算振子強度，就創建LRCC3mom.dal文件，寫入以下內容，然后進行計算。其中**INTEGRAL指定要算基函數之間哪些類型的積分，.DIPLEN要求算長度形式的偶極積分，這是因為算振子強度要算躍遷偶極矩，會用到這類積分。*CCOPA模塊用來計算耦合簇的基態到激發態的躍遷屬性，對當前用的Dalton版本支持CCS、CC2、CCSD、CC3，其中CC3僅限單重態激發態。.DIPOLE要求計算長度形式的躍遷偶極矩及振子強度。

**DALTON
.RUN WAVE FUNCTIONS
**INTEGRAL
.DIPLEN
**WAVE FUNCTION
.CC
*CC INP
.CC3
*CCEXCI
.NCCEXCI
2 1 0 1
*CCOPA
.DIPOLE
**END OF DALTON

從輸出文件（file文件包里的LRCC3mom_H2CO_TZVP.out）可以看到各個激發態的躍遷偶極矩和振子強度，比如2-1A1的如下所示，振子強度為0.34867845。這個值很接近TBE1原文表VII里這個態的LR-CC2和LR-CCSD振子強度（這篇文章里沒給CC3的振子強度），分別為0.368和0.374。

     Transition from ground state to:
     number, multiplicity, symmetry :    2  ^1A1
     frequency :   0.3839722612 a.u.    10.44842 e.V.     84272.2 cm^-1

     +-----------+-----------------+-----------------+---------------------+
     | operator  |   left moment   |  right moment   | transition strength |
     +-----------+-----------------+-----------------+---------------------+
     | XDIPLEN   |      0.00000000 |      0.00000000 |       0.00000000    |
     | YDIPLEN   |      0.00000000 |      0.00000000 |       0.00000000    |
     | ZDIPLEN   |     -0.83921304 |     -1.62309631 |       1.36212358    |
     +-----------+-----------------+-----------------+---------------------+
       oscillator strength (length gauge)   :      0.34867845

上面表格里的ZDIPLEN（偶極矩的Z分量算符）的transition strength值1.36212358對應的是躍遷偶極矩Z分量的平方。

如果你只需要躍遷偶極矩的特定分量，比如Z分量，可以把.DIPOLE改為
.OPERATOR
ZDIPLEN
這樣耗時會比計算所有分量時稍微低一點。

4 LR-CC2/CCSD/CCSDR(3)激發態計算

LR-CC2、LR-CCSD、LR-CCSDR(3)計算只需要把前面例子里的.CC3分別改成.CC2、.CCSD、.CCR(3)即可。注意LR-CCSDR(3)只能算激發能而無法算包括振子強度在內的躍遷屬性。

這些級別的激發能計算的輸入輸出文件在file文件包里的other目錄下都給了，耗時跟LR-CC3比都不值得一提。整體來說按LR-CC2、LR-CCSD、LR-CCSDR(3)的順序，結果和金標準LR-CC3的越來越接近。比如它們的1-1B1激發能分別為9.348、9.257、9.186 eV，LR-CC3的為9.183 eV。所以LR-CC3激發能很難算得動的時候用LR-CCSDR(3)是很好的平替。

5 凍核設置

Dalton默認是不凍核的。顯然恰當凍核可以節約耗時而幾乎不影響精度。這里演示對H2CO做LR-CC3計算時凍結能量最低兩個分子軌道（對應C和O的1s軌道）的方式。首先得得知它們的不可約表示，從之前B3LYP/6-31G*的輸出文件中可以看到各個占據軌道能量的負值（Koopmans定理對應的電離能，相關知識參考《正確地認識分子的能隙(gap)、HOMO和LUMO》http://www.shanxitv.org/543）。

@ Summary of DFT KT ionization potentials [eV]:

@ Symmetry 1 (A1 ) :   521.511     279.798      28.621      17.429      12.139

@ Symmetry 2 (B1 ) :    10.711

@ Symmetry 3 (B2 ) :    13.454       7.303

@ Symmetry 4 (A2 ) :  No IPs in this symmetry

顯然能量最低兩個軌道是A1。

寫一個LR-CC3激發能計算輸入文件LRCC3_FC.dal，內容如下。其中.FROIMP下面第一行指定對各個不可約表示凍結多少個能量最低占據軌道，下面第二行指定對各個不可約表示凍結多少個能量最高空軌道。當前凍結的是能量最低的兩個A1軌道。

**DALTON
.RUN WAVE FUNCTIONS
**WAVE FUNCTION
.CC
*CC INP
.CC3
.FROIMP
2 0 0 0
0 0 0 0
*CCEXCI
.NCCEXCI
2 1 0 1
**END OF DALTON

從file文件包里的此任務的輸出文件LRCC3_FC_H2CO_TZVP.out可看到，現在凍核后耗時是1分49秒，1-1B1激發能是9.18872 eV，之前沒凍核時是3分32秒，1-1B1激發能是9.18343 eV。明顯凍核對激發能的影響可完全忽略不計，而耗時大為降低，因而有重要實際意義。