文/Sobereva@北京科音  2025-Jun-26

原子、分子間的相互作用在化學體系中無處不在。化學體系中的相互作用可以按強度大體分為化學鍵相互作用和弱相互作用，也可以按作用特征大體分為共價和非共價相互作用。這些相互作用普遍在數值層面上進行研究，如研究相互作用能及其物理成份（能量分解）、幾何參數、鍵級、AIM拓撲分析考察鍵臨界點屬性，等等。而如今，這些相互作用的可視化分析（visual analysis）日趨流行，在計算化學研究文章中被應用得越來越普遍、在文獻中越來越常見。這些方法可以令化學家們十分直觀生動地認識相互作用出現的位置、強度、本質，不僅應用十分便利，也對傳統的基于數值的研究相互作用的方式提供了關鍵性的互補。這些方法在近年來還得到了快速發展，涌現了許多新的分析手段。顯然，對這些可視化分析方法進行全面、系統的綜述是極為重要且迫切的。

近期，北京科音自然科學研究中心（http://www.keinsci.com）的盧天應知名的Angew. Chem. Int. Ed.期刊的邀請，發表了名為Visualization Analysis of Covalent and Noncovalent Interactions in Real Space的綜述，可在此訪問：https://onlinelibrary.wiley.com/doi/10.1002/anie.202504895。若未訂閱，也可以免費閱讀此文在ChemRxiv上的預印版https://doi.org/10.26434/chemrxiv-2025-9t442（內容和正式版基本一樣，寫文章時請引用正式版）。

推薦使用Multiwfn程序做文中介紹的各種可視化分析時引用此綜述！

此綜述是到目前為止最完整和系統介紹各種化學體系中相互作用可視化分析方法的文章，面向廣大化學家，特別是計算化學家。文章的寫作注重令內容清晰易懂、深入淺出，且引發廣泛讀者們的興趣。通過此文，化學研究者們可以快速且充分了解各種可視化分析方法的基本思想、適用范疇，并且結合文中提供的豐富的例子，讀者能夠充分體會到它們的重要實用價值。在看過本文有了基本的知識框架后，若想進一步深究物理、算法細節等方面的知識以及了解更多實際應用，推薦再閱讀此綜述中引用的各種方法的原文和相關文章。

此綜述介紹的可視化分析方法十分全面，并且囊括了最新進展。主要介紹了NCI（亦稱RDG）、IRI、IRI-pi、IGM、IGMH、mIGM、aNCI、amIGM、Hirshfeld表面分析、靜電勢分析、范德華勢分析、ELF、LOL、LOL-pi、電子密度拉普拉斯函數、變形密度、價層電子密度。值得一提的是，其中IRI、IRI-pi、IGMH、mIGM、amIGM、范德華勢分析、價層電子密度分析都是此文的作者所提出的。

綜述中還簡要提及了文章作者開發的十分流行的波函數分析程序Multiwfn（主頁：http://www.shanxitv.org/multiwfn），這是唯一可以實現此文介紹的所有方法的程序，有很多方法還是此程序獨家支持的。Multiwfn對于大量相互作用可視化分析方法在計算學領域的普及和推廣起到了關鍵性的作用，此綜述里幾乎所有圖都是基于Multiwfn計算的數據繪制的。Multiwfn的相關信息見《Multiwfn FAQ》（http://www.shanxitv.org/452）和介紹文章J. Chem. Phys., 161, 082503 (2024) DOI: 10.1063/5.0216272。

此外，http://www.shanxitv.org/667介紹的盧天在2024年發表的Visualization Analysis of Weak Interactions in Chemical Systems（Comprehensive Computational Chemistry Vol.2中的一章）和前述的綜述有極為密切的聯系，此文介紹的方法專注于弱相互作用，涉及的方法更少，而由于有明顯更多的篇幅，因而對它們的介紹明顯更深入細致，故和上面的綜述有極大的互補性，強烈建議閱讀！

還十分推薦閱讀以下文章及其中引用的文章獲取更多相關知識以及計算操作細節：
使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用
http://www.shanxitv.org/621（http://bbs.keinsci.com/thread-28147-1-1.html）
使用IRI方法圖形化考察化學體系中的化學鍵和弱相互作用
http://www.shanxitv.org/598（http://bbs.keinsci.com/thread-23457-1-1.html）
使用Multiwfn做Hirshfeld surface分析直觀展現分子晶體和復合物中的相互作用
http://www.shanxitv.org/701（http://bbs.keinsci.com/thread-43178-1-1.html）
談談范德華勢以及在Multiwfn中的計算、分析和繪制
http://www.shanxitv.org/551（http://bbs.keinsci.com/thread-17271-1-1.html）
在Multiwfn中單獨考察pi電子結構特征
http://www.shanxitv.org/432（http://bbs.keinsci.com/thread-10610-1-1.html）
ELF綜述和重要文獻小合集
http://bbs.keinsci.com/thread-2100-1-1.html
靜電勢與平均局部離子化能相關資料合集
http://bbs.keinsci.com/thread-219-1-1.html
使用Multiwfn研究分子動力學中的弱相互作用
http://www.shanxitv.org/186
使用Multiwfn做aNCI分析圖形化考察動態過程中的蛋白-配體間的相互作用
http://www.shanxitv.org/591（http://bbs.keinsci.com/thread-21826-1-1.html）
通過電子定域化函數(ELF)、價層電子密度分析討論親核進攻的方向
http://www.shanxitv.org/606（http://bbs.keinsci.com/thread-23982-1-1.html）
通過價層電子密度分析展現分子電子結構
http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB201709252

此外，筆者的18碳環及相關體系的研究文章中充分利用了文中介紹的許多方法，是很好的應用實例，研究匯總見http://www.shanxitv.org/carbon_ring.html。

量子化學波函數分析與Multiwfn程序培訓班（http://www.keinsci.com/WFN）對前述的綜述中幾乎所有方法都有十分全面、透徹、系統的介紹，并對各種分析方法在Multiwfn中的實際操作做了細致的演示并傳授大量技巧，信息量遠比綜述里的大得多得多，因而很推薦對這些分析方法感興趣者關注和參加。

下面是Visualization Analysis of Covalent and Noncovalent Interactions in Real Space綜述中的圖，供預覽。

文/Sobereva@北京科音   2025-Jun-11

0 前言

NICS和磁感生電流是非常常用的磁響應性質一類的考察芳香性的方法，在《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）里有簡要介紹并里面附了大量我寫過的相關博文，而在“量子化學波函數分析與Multiwfn程序培訓班”（http://www.keinsci.com/WFN）的“芳香性與離域性分析”一節里有極為全面深入的介紹。很多人對這兩個方法在衡量芳香性方面的認識局限在：
(1)NICS(1)ZZ約為0說明環是非芳香性，明顯為負說明環是芳香性且越負芳香性越強，明顯為正說明環是反芳香性且越正反芳香性越強。注意NICS的形式很多，這里僅以最常用的NICS(1)ZZ舉例
(2)感生電流的方向如果滿足左手定則（diatropic電流），電流越大芳香性越強；如果與左手定則方向相反（paratropic電流），電流越大反芳香性越強。電流大小一般通過對鍵截面上的電流密度做積分來衡量，稱為bond current strength (BCS)
但實際上，僅光憑以上知識討論芳香性或反芳香性，在一些情況下會得到不正確甚至嚴重錯誤的結論。在本文筆者介紹一些大多數人都普遍沒有注意到的用NICS和感生電流討論芳香性方面的重要問題。把這些問題都徹底搞清楚了，基于磁屬性研究芳香性時就不懼審稿人comment了。

本文中的計算用的幾何結構，若無專門提及，皆是B3LYP/6-31G*下優化的無虛頻結構。

1 NICS、感生電流與HOMO-LUMO gap的關系

感生電流和NICS計算公式都涉及到對占據軌道和空軌道的循環，并且空軌道與占據軌道能量之差出現在分母項中，這一點看《深入理解分子軌道對磁感生電流的貢獻》（http://www.shanxitv.org/703）里的公式1、2、3便知。總的來說，HOMO-LUMO gap越大，空軌道與占據軌道整體能量差越大，導致NICS和感生電流越接近0。因此，如果兩個體系的gap相差較大，不應直接憑NICS或BCS判斷二者的芳香性強弱，還應該把gap差異對結果的影響考慮進去。在《深度揭示互為等電子體的苯、無機苯和carborazine的芳香性的顯著差異》（http://www.shanxitv.org/731）里介紹的筆者的Chem. Eur. J., 30, e202403369 (2024)一文中就注意到了這一點，文中研究的carborazine的HOMO-LUMO gap（7.30 eV）明顯小于苯（11.20 eV）和無機苯（12.18 eV），因此文中提到：It is noteworthy that the aromaticity of carborazine compared to benzene and borazine is more or less quantitatively overestimated by the value of BCS and NICS。

4n電子的Huckel反芳香性體系的HOMO-LUMO gap比4n+2的Huckel芳香性體系的整體更小，Chem. Soc. Rev., 44, 6597 (2015)認為這是前者的BCS比后者更大的原因。例如反芳香性的環丁二烯和芳香性的苯在B3LYP/6-31G*級別算的HOMO-LUMO gap分別為3.692和6.800 eV，使用《使用SYSMOIC程序繪制磁感生電流圖和計算鍵電流強度》（http://www.shanxitv.org/702）介紹的SYSMOIC算出來的環丁二烯的BCS為-0.7941（長C-C鍵）和-0.7499 a.u.（短C-C鍵），而苯的C-C鍵的BCS大小則小得多，只有0.4216 a.u.。

如《正確地認識分子的能隙(gap)、HOMO和LUMO》（http://www.shanxitv.org/543）所說，HF成份越低的泛函算出來的HOMO-LUMO gap越低，因此若發現HF成份低的泛函算出來的NICS和BCS的絕對值較大，大概率也是上述原因（但并非HF成份越低的泛函算出來的總是越大，因為泛函的差異還影響到其它方面）。

2 NICS的環尺寸的依賴性

原理上，感生電流比NICS判斷芳香性更合理，因為NICS對環尺寸有依賴性。雖然有的文章說NICS的環尺寸依賴性弱，但那往往只是對較小的環之間對比做的結論，再加上本身不同體系的芳香性就有所不同，令那種說法更不確切。

電流對特定位置產生的磁場可以用Biot–Savart方程得到。如果是環形電流，則在環中心產生的磁場為μ0*I/(2*R)，其中μ0是真空磁導率，I是電流大小，R是環的半徑，更多信息見https://en.wikipedia.org/wiki/Biot%E2%80%93Savart_law。因此，對芳香性的情況，當感生電流相同時，對應芳香性的環形路徑整體半徑越大、原子距離環中心整體越遠，在環中心產生的磁屏蔽效果就越弱、NICS的絕對值越小。如果認為感生電流I的大小是芳香性強弱的準確體現，那么顯然用NICS對比環尺寸明顯不同的體系的芳香性的大小就是不甚嚴格的，因為對越大的環會越低估芳香性。所以，如果要基于磁屬性對比明顯不同大小的環的芳香性，建議用感生電流而非NICS說事。對于反芳香性，也有和上面完全類似的情況。

3 感生電流對積分截面的依賴性

對應芳香性的環上，如果只考慮對芳香性有貢獻的電子的感生電流，原理上來說在垂直于鍵（更確切來說是垂直于電流方向）的各個截面上積分電流密度得到的BCS都應該是相同的。但實際得到的BCS沒有這么理想化，環上的各個鍵的BCS往往不完全相同，在于：
(1)積分的數值誤差、有限的積分截面尺寸
(2)積分截面可能并不嚴格垂直于電流方向
(3)在一些鍵上存在方向彼此相反的感生電流的抵消作用，此效應對于多環芳烴等并環的情況尤為明顯
(4)環周圍的原子的感生電流產生的干擾，以及環上原子的定域化的電子（內核電子、sigma電子、孤對電子等）產生的局部感生電流添亂，這導致計算的BCS對鍵上的積分截面位置的選取可能存在不可忽視的依賴性
(5)基組不夠完備，導致感生電流計算不夠準確。對D2h點群的環丁二烯（反芳香性）極小點結構，用SYSMOIC在B3LYP/6-31G*下計算的短和長C-C鍵的BCS分別為0.7499和0.7941 a.u.，而基組用完備得多的def2-QZVP計算結果分別為0.7116和0.7111 a.u.，可見后者明顯更等價

若發現要考察芳香性的環上的各個鍵的BCS不等價性不可忽視，可以取感生電流較明確、被上述的(4)干擾小的鍵的BCS作為環電流值衡量芳香性，或者對環上的鍵的BCS取平均作為環電流值（適合碳環、環丁二烯等情況）。此外，也可以嘗試計算BCS時只考慮對芳香性有貢獻的電子所占據的分子軌道，例如用感生電流研究pi芳香性時就只考慮pi軌道。

4 感生電流判斷芳香性和反芳香性方面相對于NICS的優勢

NICS很有可能給出關于芳香性錯誤的結論，不如觀看和積分感生電流，以及觀看《通過Multiwfn繪制等化學屏蔽表面(ICSS)研究芳香性》（http://www.shanxitv.org/216）里介紹的ICSS圖來得靠譜，這一節給出幾個例子。

JCTC, 6, 1131 (2010)發現NICS用于考察氟化氫三聚體(HF)3的芳香性有問題。在wB97XD/def-TZVP下優化出的此體系的結構如下所示。憑基本常識就知道這個體系是非芳香性的，但是B3LYP/def2-TZVP算出來的NICS(0)ZZ為13.52 ppm，仿佛有顯著反芳香性似的。不過NICS(1)ZZ為-0.04 ppm，這倒是正確體現出此體系并沒有反芳香性。

為什么(HF)3的NICS(0)ZZ那么正？這可以從此體系所在平面上的感生電流圖上理解，如下圖所示（磁場垂直于屏幕向外）。可見此體系并沒有形成環繞整體的感生電流，但由于存在環繞各個HF的sigma鍵的dia電流，導致在環中央區域呈現出等效的para電流的特征（紅色箭頭所示），使得環中央是去屏蔽的，因而NICS(0)ZZ為正。而由于這種效應在垂直于體系平面方向上衰減很快，因此NICS(1)ZZ幾乎為0。這也體現出，對于小環體系，由于sigma鍵的感生電流產生的磁場離環中心太近，NICS(0)ZZ用來衡量芳香性是極差的做法。但對于18碳環及其等電子體B6N6C6這樣的大環倒不是什么問題，因為局域的電子對應的局部的感生電流造成的磁場在環中心處都基本衰減沒了，這也是為什么《18碳環等電子體B6N6C6獨特的芳香性：揭示碳原子橋聯硼-氮對電子離域的關鍵影響》（http://www.shanxitv.org/696）介紹的文章中用了NICS(0)ZZ討論。

實際上對于苯分子也有類似情況。下圖是苯分子的pi以外的電子產生的感生電流圖，可見由于局域的sigma電子產生的局部dia電流的存在，在環中心也有等效的para電流（紅色箭頭所示）并使得環中心產生磁去屏蔽效應。不過由于苯的pi電子造成的全局的dia電流在環中心產生的磁屏蔽效應遠強于para電流的磁去屏蔽效應，因此苯的NICS(0)ZZ還是負值。

我極其推薦閱讀《使用Multiwfn繪制一維NICS曲線并通過積分衡量芳香性》（http://www.shanxitv.org/681）中的討論，了解sigma和pi電子對NICS_ZZ的貢獻在垂直于環方向上的變化特征。

《18個氮原子組成的環狀分子長什么樣？一篇文章全面揭示18氮環的特征！》（http://www.shanxitv.org/725）里介紹的我的ChemPhysChem, 25, e202400377 (2024)文章中，對18氮環的最穩定構型計算了NICS(0)ZZ，為3.42 ppm，也是表面上體現出了輕微的反芳香性。然而通過ICSS等值面圖和AICD感生電流圖，都充分確認了此體系沒有芳香性或反芳香性特征，詳見博文里的圖。這也是個很好地體現NICS(0)ZZ展現（反）芳香性不靠譜的例子。

NICS亦有可能把非芳香性體系誤判成芳香性的。Theor. Chem. Acc., 134, 8 (2015)指出NICS衡量很多過渡金屬團簇的芳香性也是不合理的。其中一個例子是如下所示的Ti3(CO)3。基于PCCP, 13, 4576 (2011)的SI里給出的BP86/6-311+G(d)級別下優化的結構，我用PBE0/def2-TZVP算出的NICS(0)ZZ和NICS(1)ZZ分別是-141.1和-11.7 ppm，乍看起來明顯芳香性。

下圖左側和右側分別是Ti3(CO)3的分子平面上和分子平面上方1埃處的感生電流，數值過大的部分截掉了，粉色大圓球是Ti，磁場垂直屏幕向外。由圖可見，無論是在分子平面上還是在其上方1埃處，都有顯著的環繞Ti原子的para電流，很大程度等效帶來了在環中心及其上方的dia電流的效果（紅色箭頭所示），即令環中心及其上方一定區域都是磁屏蔽的，這是NICS(0)ZZ和NICS(1)ZZ都為負值的原因。然而從圖上并沒有看到環繞整體的連貫的感生電流，因此Ti3(CO)3就是非芳香性的。可見此例不僅NICS(0)ZZ，連NICS(1)ZZ都誤判了芳香性。

5 顯著且連貫的dia感生電流未必能證明有明顯的芳香性

從上一節看到，感生電流判斷芳香性比個別位置的NICS更為嚴格且直觀，但即便是感生電流，也未必總能正確衡量芳香性，這里以[N6H6]2+體系為例進行說明，更多討論見ChemistrySelect, 2, 863 (2017)。下圖展示了此體系的平面構型

此構型并不是極小點，而是有很多破壞平面的虛頻，并且其能量甚至比其開環后的線型結構的還要高，因此明顯沒有芳香穩定化作用。用Multiwfn計算出它的多中心鍵級幾乎為0（0.00035），也進一步體現出沒有pi共軛。如果用《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）的做法繪制它的LOL-pi圖，也會看到各個氮都有定域的孤對電子，根本沒有像苯的pi電子一樣在環上充分離域。然而，下圖所示的[N6H6]2+平面構型的分子上方1埃處的感生電流（磁場垂直于屏幕向外）則是完全順時針的dia電流，和具有顯著芳香性的苯的情況看起來完全一樣！

為什么[N6H6]2+的平面構型不是極小點？這是因為如下圖的軌道占據情況所示，盡管此體系表面上有10個pi電子，看似滿足Huckel的4n+2芳香性規則因而芳香性會使此構型被穩定化，但e2u對稱性的軌道的反pi特征顯著（對于苯來說是沒有占據的），因此pi共軛作用很弱，再加上如此眾多的孤對電子間的強烈互斥作用，顯然令其極不穩定。

為什么[N6H6]2+的平面構型像普通芳香性體系一樣有顯著的貫穿整個環的dia電流？這需要利用《深入理解分子軌道對磁感生電流的貢獻》（http://www.shanxitv.org/703）里介紹的知識，沒看過者強烈建議認真閱讀。此文說了，感生電流可視為由各種占據軌道向各種空軌道躍遷所貢獻，上圖藍色箭頭標注的躍遷是此結構下對感生電流貢獻最大的躍遷（因為這種躍遷能最小），其中非占據的b2g軌道比占據的e2u軌道軌道多一個節面，這種躍遷是平動允許的，因此對dia電流有貢獻。這和造成苯分子的dia電流的機制是完全相同的，只不過苯分子的dia電流主要是e1g到e2u躍遷所貢獻的。

由此例可見，光靠是否存在全局的dia感生電流判斷芳香性并不總是可靠！因為這種電流是否存在，和體系是否真的有芳香性特征并沒有內在的必然聯系，而與它關系真正最直接、最緊密的是軌道的對稱性特征！

可能大多數人覺得，有芳香性的體系有滿足左手定則的dia電流的原因在于：某個環有芳香性意味著環上有明顯離域的電子，因此在外磁場的作用下會像經典電磁學描述的導體一樣在環上產生dia電流。但實際上真正本質并非如此，不能忽視關鍵性的量子力學效應！一般的芳香性體系有dia電流實際上是一個表象，在于它們的軌道對稱特征和躍遷正好能滿足產生顯著dia電流的條件，而能產生顯著dia電流的并非只有芳香性體系。因此，要這樣認為：有顯著的貫穿整個環的dia電流是這個環具有芳香性的必要非充分條件。

由于磁屏蔽特征是感生電流的進一步衍生屬性，[N6H6]2+的平面構型既然都和苯一樣具有相同感生電流特征了，其自然NICS(0)ZZ和NICS(1)ZZ也都為負，分別為-178.1和-24.6 ppm。

此例彰顯出判斷芳香性最好不要光基于磁屬性，即便這種判斷芳香性的做法極為流行。同時結合多中心鍵級等基于電子結構的方法，以及結合芳香穩定化能或離域化能這種能量層面的方法，結論會更可靠、更有說服力。本節的[N6H6]2+的例子也體現出，形式上滿足4n+2規則的體系未必真的有芳香性，若電子都沒法在環上充分離域、不構成共軛體系，滿足4n+2也仍是非芳香性的。

6 顯著且連貫的para感生電流未必能證明有明顯的反芳香性

J. Org. Chem., 88, 14831 (2023)通過一系列分子指出，即便環上存在顯著且連貫的para電流，也未必能可靠地證明這個環就有明顯的反芳香性。這一節就用其中最簡單的體系，苯+1電荷的陽離子自由基（具有5個pi電子）和-1電荷的陰離子自由基（具有7個pi電子）來說明這一點，這種體系的（反）芳香性沒法基于經典的Huckel規則判斷。

如下所示，JOC文中給出了苯陽離子、中性苯、苯陰離子的分子平面上方1 Bohr處的感生電流，磁場垂直于屏幕向外，計算用的是對各自狀態優化的幾何結構，電場強度越大箭頭越藍、越小越紅。可見，苯陽離子和苯陰離子在碳環上方的電流都是逆時針的para電流，與中性苯的情況相反，乍看起來苯的陰、陽離子都是顯著反芳香性的。

然而，JOC這篇文章計算了苯的陰、陽離子的芳香穩定化能（用到的反應見文章的SI），發現這兩個體系是有芳香穩定化作用的，因此從能量角度來說是有芳香性的。我還通過Multiwfn在文中給的苯離子的極小點結構下在B3LYP/6-31G*級別計算了六中心鍵級，苯的陽離子和陰離子的結果分別為0.044和0.045，只有中性苯的0.086的約一半。具體來說，苯陽離子和陰離子的六中心鍵級分別幾乎來自于alpha和beta電子的貢獻（且和中性苯的相應自旋電子的貢獻差不多），而另一個自旋幾乎沒貢獻。也就是說，相對于苯，離子狀態下電子數改變1的那個自旋的電子對芳香性沒明顯貢獻，而電子數未變的自旋的電子對芳香性的貢獻基本保持不變。

由于從能量和電子結構角度都確認了苯的陰、陽離子是芳香性的，那么為什么它們的感生電流是para的，顯得像普通反芳香性體系一樣？JOC文中給了下面的圖對感生電流進行了解釋。由于Jahn-Teller效應，陽離子（5π）和陰離子（7π）體系的對應于苯的二重簡并的pi軌道發生了輕微分裂。圖中藍色大箭頭標注了對感生電流產生主導的躍遷。對苯離子體系的情況，由于能級分裂造成的能量差較小，明顯小于1→2的躍遷，因此對陽離子和陰離子分別是1→1和2→2躍遷對感生電流產生了主導性的貢獻。由于這種躍遷涉及的兩個軌道的節面數相同，是轉動允許的，因此如《深入理解分子軌道對磁感生電流的貢獻》所述，貢獻的是para電流。

由本節和上一節的例子可見，雖然感生電流對大多數體系判斷芳香性是合理的，而且比起只看個別位置的NICS_ZZ可靠不少，但在極個別體系上還是能暴露出其局限性。因此若有可能，除了基于磁屬性分析外，還是建議再結合一兩種靠譜的手段，特別是可靠又容易計算的多中心鍵級，來衡量芳香性，減少誤判芳香性、被審稿人懟的可能（本來芳香性這方面的“說法”和“觀點”就相當多）。

但也沒必要盲目用過多的方法衡量芳香性，否則一方面太冗余，另一方面本身有的芳香性指標就不那么靠譜、普適，納入討論反倒無益于揭示芳香性的真相。有文章認為鍵的均衡化、出現dia和para的感生電流，屬于芳香性的二階（次級）特征，而能量的穩定化和電子的離域，才算是芳香性的一階（即最關乎本質的）特征，這個觀點我是認可的。

文/Sobereva@北京科音   2025-May-31

0 前言

18碳環（cyclo[18]carbon）自從于2019年在凝聚相實驗中被觀測到后，引發了化學家們非常廣泛的興趣，目前已有巨量量子化學研究文章發表，筆者對18碳環及相關體系的理論研究工作匯總見http://www.shanxitv.org/carbon_ring.html。完全由氮、硼原子交替構成的環B9N9，以及18碳環部分由硼、氮取代的結構，也已被研究過，例如《18碳環等電子體B6N6C6獨特的芳香性：揭示碳原子橋聯硼-氮對電子離域的關鍵影響》（http://www.shanxitv.org/696）介紹的筆者的Inorg. Chem., 62, 19986 (2023)文章。而硼氮替換如何影響碳環的光學性質，之前尚未有系統的研究。

近期，江蘇科技大學的劉澤玉等人和北京科音自然科學研究中心的盧天，共同在知名的Chemical Engineering Journal (CEJ)刊物上發表文章，全面揭示了硼氮摻雜對碳環的光學性質的影響，并發現恰到好處的摻雜能夠獲得具有出色非線型光學性質的物質。此文的研究也對硼氮摻雜改性大共軛碳材料提供了十分有益的啟示。非常歡迎閱讀以及引用此文：

Xiaohui Chen, Zhibo Xie, Jiaojiao Wang, Wenwen Zhao, Xiufen Yan, Zeyu Liu,* Cai Ning,* Tian Lu,* Obtaining (BN)4C10 with excellent optoelectronic properties by screening boron-nitrogen substituted cyclo[18]carbon, (BN)nC(18-2n) (n = 1–9), Chem. Eng. J., 515, 163236 (2025) https://doi.org/10.1016/j.cej.2025.163236

在2025年6月29日之前可以通過https://authors.elsevier.com/c/1l44R4x7R2o3f3免費閱讀此文。

下面，本文就對這篇Chem. Eng. J.文章的研究工作的主要內容進行深入淺出的介紹，并對不少細節和要點做一些附加說明，以便讀者更順利地了解文章的最關鍵的內容，以及能在其它研究中借鑒此文的研究方式。更具體的研究結果和更多的討論請閱讀原文。下文的圖片皆來自原文的正文或補充材料。

1 (BN)nC(18-2n)的幾何結構

此文研究的對像是不同數目硼氮（BN）單元取代的18碳環，通式為(BN)nC(18-2n)，其中n=1-9，各個BN單元連續排布。因此，此文研究了從氮化硼環B9N9到18碳環的完整過渡過程。在ωB97XD/6-311+G(2d)下優化出的體系的極小點結構如下所示，紅色、藍色、灰色分別對應硼、氮、碳，結構均為純平面。ωB97XD在http://www.shanxitv.org/carbon_ring.html列舉的各種碳環及衍生物的研究中都被證實可以很合理描述18碳環及其硼、氮摻雜體系的電子和幾何結構。

2 (BN)nC(18-2n)的電子結構

下圖展示出了不同的硼氮取代的18碳環的HOMO、LUMO能級和HOMO-LUMO gap，通過Multiwfn結合VMD按照《用VMD繪制藝術級軌道等值面圖的方法）（http://www.shanxitv.org/449）方法繪制的等值面圖也給出了。可見隨著BN單元數n的增加，LUMO先輕微減小然后再明顯增大，而HOMO先輕微增加再明顯減小，這使得HOMO-LUMO gap在n=3時達到了最小值。

如《使用Multiwfn基于完全態求和(SOS)方法計算極化率和超極化率》（http://www.shanxitv.org/232）里的公式所示，激發能出現在計算極化率的SOS公式的分母部分，而HOMO-LUMO gap越小很大程度暗示激發能整體越小，因此由此可以預期HOMO-LUMO gap最小的(BN)3C12或相鄰的(BN)2C14、(BN)4C10具有最大的極化率，這與后文的實際計算結果一致。

當n不很接近9時，從上圖可以看到HOMO、LUMO幾乎都出現在碳的部分，這和下圖（原文中圖S1）所示的基于《使用Multiwfn繪制態密度(DOS)圖考察電子結構》（http://www.shanxitv.org/482）介紹的方法繪制的PDOS圖相吻合，碳構成的片段的PDOS總是主導最高占據能級和最低非占據能級。

以上信息直接暗示出電子激發最容易的區域、反應活性最高的區域、對（超）極化率貢獻最大的部分，都是(BN)nC(18-2n)中的共軛碳鏈部分（至少對于n不很接近9時）。

從18碳環開始，隨著取代的BN單元數逐漸增加，HOMO-LUMO gap先輕微下降然后逐漸顯著上升的本質原因何在？文中認為原因在于以下幾點：
? 有顯著芳香性的體系的gap傾向于比芳香性更弱的相似物的更大，例如《深度揭示互為等電子體的苯、無機苯和carborazine的芳香性的顯著差異》（http://www.shanxitv.org/731）中提到的，具有很強芳香性的苯的gap比芳香性更差的carborazine的gap明顯更大就是這種情況。18碳環有明顯的芳香性，這在Carbon, 165, 468 (2020)中筆者已經充分論證過了。當18碳環的C-C逐漸被BN單元所取代后，芳香性被逐漸破壞，導致gap減小。
? 當BN單元數變得略多時，上面的效應弱化，而導致gap增大的一個因素變成了主導，使得(BN)3C12具有最小的gap。這個因素是隨著BN單元數增加、環上的碳鏈區域變短，導致pi共軛區域的長度逐漸減小，從而令gap變大。這類似于一維無限深勢阱，當勢阱越窄，相鄰能級差會越大。在《正確地認識分子的能隙(gap)、HOMO和LUMO》（http://www.shanxitv.org/543）里給出過共軛聚合物的例子，也是共軛空間范圍越窄，gap越大。
? 另一方面，當BN單元數已經很多并進一步增加時，體系也越來越接近B9N9，自然gap也越來越接近B9N9。缺乏pi共軛的B9N9的gap顯著大于18碳環，因此gap在BN單元數增加的后期提升得很快。

再來看永久偶極矩（μ0）、垂直電離能（VIP）、垂直電子親和能（VEA）、電子硬度（η）隨BN單元數的變化，如下所示。可見VIP和VEA的變化趨勢分別與HOMO能量和LUMO能量正相反，因為按照Koopmans定理，VIP和VEA分別近似等于HOMO和LUMO能量的負值（雖然這個近似的精度往往很爛，但可以解釋趨勢）。硬度定義為VIP-VEA，數值越小電子軟度越大，一定程度暗示反應活性越高，在Koopmans近似下等于HOMO-LUMO gap。可見硬度的變化趨勢確實和HOMO-LUMO gap一樣，也在(BN)3C12處達到了最小值。

由于18碳環和B9N9都是中心對稱的，所以它們的永久偶極矩都為0。如上所示，在BN單元增加的過程中偶極矩先增大后減小。有趣的是，HOMO-LUMO gap最小的(BN)3C12正好擁有最大的偶極矩。

3 (BN)nC(18-2n)的（超）極化率

此文的（超）極化率使用Gaussian 16在ωB97XD/aug-cc-pVTZ級別下計算，數據的提取和相關物理量通過《使用Multiwfn分析Gaussian的極化率、超極化率的輸出》（http://www.shanxitv.org/231）介紹的方式實現。此級別下計算的各個(BN)nC(18-2n)體系的靜態極化率如下，分子在XY平面上。曾經筆者在Carbon, 165, 461 (2020)中全面深入研究了18碳環的光學性質，并指出了18碳環在平面內具有巨大的極化率。當前研究進一步體現出，隨著不具備pi共軛特征的BN單元越多、pi共軛的碳鏈區域越短、電子分布在外場驅動下能夠輕易轉移的范圍越狹窄，平面內的總極化率（αxx+αyy）逐漸顯著減小，最終的B9N9的平面內總極化率遠小于18碳環的。另外由圖可見，垂直于分子平面方向的極化率αzz基本不依賴于BN單元數，且B9N9和18碳環的基本沒區別，分別為98和99 a.u.，體現出這兩種等電子體雖然元素不同，但電子結構層面的差異并不會明顯造成在垂直于環平面方向上電子被外電場極化能力的差異。下圖的αani是極化率的各向異性，可見隨著從18碳環逐漸變向B9N9，αani逐漸減小，這是因為B9N9不具備18碳環那樣的在環上的全局pi共軛，因此在不同方向上的極化率差異明顯更不顯著。

文中使用《使用Multiwfn極為方便地繪制(超)極化率密度和三維空間對(超)極化率的貢獻》（http://www.shanxitv.org/683）的方法繪制了極化率密度，其中下面這張圖展現了(BN)4C10的垂直于環平面的極化率分量的極化率密度的±0.12 a.u.等值面，等值面包裹的區域是對αzz有主要貢獻的區域。可見主要貢獻的部分是氮的孤對原子，以及碳鏈部分的pi電子很豐富的較短的C-C鍵。注意BN取代的18碳環和原本的18碳環一樣都具有C-C鍵鍵長交替的特征。

(BN)nC(18-2n)的第一超極化率的總值（βtot）、在偶極矩方向的投影的大小|βvec|，以及《使用Multiwfn計算與超瑞利散射(HRS)實驗相關的量》（http://www.shanxitv.org/499）介紹的超瑞麗散射（βHRS）實驗測量的βHRS在下圖給出了。由于18碳環和B9N9是中心對稱的，二者的β皆精確為0。從18碳環開始，隨著BN單元數的增加，各種β都是先增大后減小。這體現出BN單元數可以對體系的非線型光學（NLO）性質起到充分的調控作用，這個發現給設計具有特定性能的NLO材料提供了一種新思路。其中，(BN)4C10具有最大的βtot和βHRS，并且數值相當大，這是本文理論篩選出的具有關鍵價值的分子！

文章還使用《使用Multiwfn通過單位球面表示法圖形化考察（超）極化率張量》（http://www.shanxitv.org/547）介紹的方法用Multiwfn結合VMD對上述最具特殊性的(BN)4C10繪制了單位球面表示法圖像。其中的小箭頭展現出無論兩個電磁場同時向什么方向打向此體系，耦合作用產生的誘導偶極基本上都是順著綠色箭頭出現的，這也正是順著碳鏈pi共軛的主體方向。

此外，文章還考察了外場頻率對極化率和第一超極化率的影響，證實了外場頻率越大，各種(BN)nC(18-2n)的第一超極化率都會明顯越大，即具有正的頻率-色散效應。

4 (BN)nC(18-2n)的電子吸收光譜

文章全面考察了不同的(BN)nC(18-2n)的電子吸收光譜，如下圖左側所示，可見它們在可見光范圍內都沒有明顯的吸收，說明是無色物質。隨著BN單元數的增多，光譜明顯平滑地從18碳環向B9N9逐漸過渡，變化細節見原文的討論。對于前述很有特點的(BN)4C10，此文在補充材料里給出了按照《使用Multiwfn繪制電荷轉移光譜(CTS)直觀分析電子光譜內在特征》（http://www.shanxitv.org/628）介紹的CTS方法將總光譜特征進行分解的圖，如下圖右側所示。可見在300多nm部分的峰主要來自于碳鏈區域內的電子重排的貢獻，而200nm左右的吸收的成份更為復雜，同時牽扯到碳鏈和BN單元部分。

上面這個300多nm的峰對應于激發波長為323.8 nm的S0-S7激發，由于它沒有主導性的單一分子軌道躍遷，因此文中按照《使用Multiwfn做空穴-電子分析全面考察電子激發特征》（http://www.shanxitv.org/434）介紹的極為普適的方法對其繪制了空穴和電子分布圖，如下所示，紫色和黃色分別對應于空穴和電子分布。可見此激發的空穴和電子確實主要分布在碳鏈部分，僅有極少部分分布在氮和硼上，因此可以基本指認為碳鏈部分的局域激發。值得一提的是，按照《談談計算第一超極化率的雙能級公式》（http://www.shanxitv.org/361）和《使用Multiwfn對第一超極化率做雙能級和三能級模型分析》（http://www.shanxitv.org/512）里說的雙能級分析，S7是關鍵態，對(BN)4C10的第一超極化率有最大貢獻，他是激發能最低的振子強度明顯的態（f=0.603），而且這種激發會帶來很顯著的偶極矩的變化（Δμ=2.27 Debye）。

文中還研究了溶劑效應如何影響（超）極化率和電子光譜，結果見原文。

5 總結

本文介紹的近期發表的Chem. Eng. J., 515, 163236 (2025)一文通過嚴謹的量子化學計算并結合Multiwfn做波函數分析，非常全面、系統地探究了理論設計的硼氮取代的18碳環，即(BN)nC(18-2n) n=1-9。此文展現了硼氮單元數目是如何影響體系的電子結構、（超）極化率和電子光譜特征的。本文的工作還篩選出了具有優秀非線型光學性質的(BN)4C10，并對它的電子激發和光譜特點進行了具體的分析。論文還專門有一節4. Prospects for molecular crystal of (BN)4C10，其中展望了由(BN)4C10構成的分子晶體的各方面特征以及可能的制備方式。本文的研究對于硼氮摻雜的碳環材料的潛在應用提供了重要的理論指導，同時本文的研究方式對于其它同類的研究也是很有價值的參考。歡迎讀者仔細閱讀原文了解更多信息。

文/Sobereva@北京科音   2025-Apr-30

0 前言

18碳環（cyclo[18]carbon, C18）是由18個sp雜化的碳原子相連構成的環狀物質，它以其獨特的幾何和電子結構，自2019年在凝聚相中觀測到后引發了化學界的巨大關注。在2023年的時候16碳環（下文簡稱C16）也被以類似方式制備了出來，見Nature, 623, 977 (2023)。C16和C18一樣都有in-plane和out-of-plane兩套pi電子，即pi(in)和pi(out)，詳見Carbon, 165, 468 (2020)里的介紹。C18每套pi電子有18個電子，滿足Huckel的4n+2規則因此基態具有雙芳香性，筆者已在Carbon, 165, 468 (2020)文中通過詳盡的波函數分析嚴格、全面證實了這一點。C16每套pi電子有16個電子，滿足Huckel 4n規則因此理應基態具有雙反芳香性。然而，激發態芳香性和基態芳香性往往是具有極大差別的。根據Baird規則，筆者意識到若讓C16的pi(in)電子和pi(out)電子都被激發的話，極有可能處于五重態激發態的C16將和C18基態一樣具有顯著的雙芳香性，而C16的三重態激發態則可能同時內在兼具芳香性和反芳香性！顯然，研究非同尋常的C16的基態和激發態芳香性特征是一個極為新穎、有趣的課題。

基于以上思路，北京科音自然科學研究中心的盧天和江蘇科技大學的劉澤玉等人共同深入全面研究了16碳環的單重態基態（S0）、最低三個三重態激發態（T1、T2、T3）和最低五重態激發態（Q1）的芳香性，成果近期發表在了知名的Chemistry - A European Journal期刊上，非常歡迎閱讀：
Yang Wu, Jingbo Xu, Xiufen Yan, Mengdi Zhao, Zeyu Liu*, Tian Lu*, Interesting Aromaticity of Franck-Condon Excited States of Cyclo[16]carbon, Chem. Eur. J., 31, e202404138 (2025) DOI: 10.1002/chem.202404138
可以通過此鏈接免費在線閱覽：https://onlinelibrary.wiley.com/share/author/SB6H76ATMQF2RNABQAM2?target=10.1002/chem.202404138

此文的Table of Contents如下，體現了此文的關鍵性結果

筆者此前還對碳單環體系及其衍生物的各方面特征做過十分廣泛的理論研究，發表的大量相關論文和深入淺出的介紹博文匯總見http://www.shanxitv.org/carbon_ring.html（不斷更新）。

下面將對這篇Chem. Eur. J.文章的主要內容進行淺顯易懂的介紹，便于讀者更容易理解文章的研究結果，同時額外附上不少分析和計算細節以幫助讀者能夠將文中的研究手段舉一反三運用到自己的研究中。原文里還有很多細節信息和討論，故請閱讀下文后閱讀原文。此文不僅研究C16芳香性本身，文章在分析思想和方法學方面對于其它分子的芳香性研究也很有借鑒意義。

1 一些研究和計算細節

介紹結果之前，有一些研究和計算細節值得說明。

C16的基態和激發態勢能面的極小點結構是不一樣的，因此激發態極小點結構下的激發態芳香性與基態極小點結構下激發態的芳香性也是不同的，而且可以相差很大。此文專注于研究C16基態極小點結構下的基態、三重態和五重態芳香性，這避免了引入激發態勢能面上結構的弛豫造成的芳香性改變。而若把結構弛豫效應考慮進去會導致分析討論變得復雜，也無法“干凈”地體現電子激發的改變本身對電子結構的直接影響。激發態勢能面上的基態結構對應的那個位置被稱為Franck-Condon（FC）點，因此此文研究的C16的激發態芳香性被稱為Franck-Condon激發態芳香性。

此文對C16做的計算主要使用Gaussian 16程序在ωB97XD/def2-TZVP級別下做，這在http://www.shanxitv.org/carbon_ring.html中匯總的筆者的巨量和碳環有關的文章中都被證實可以很穩健地描述碳環類體系的幾何和電子結構。文中對C16還做了衡量靜態相關常用的T1診斷，對S0、T1、Q1態結果分別為0.015、0.017、0.013，都是較小的值，體現出使用更復雜昂貴的多參考方法研究C16不是必須的。

此文計算C16的T1/T2/T3和Q1態時用的是ΔSCF方法結合UKS計算，而不是算激發態常用的TDDFT方法，這是因為在TDDFT下沒法算磁性質來考察芳香性，而且也沒法要求算五重態；而用ΔSCF不僅沒有這個局限性，而且還能考慮軌道弛豫效應，因此原理上比TDDFT更理想。

ΔSCF方式算Q1時直接把自旋多重度設成5做DFT計算即可，通過用Multiwfn程序查看Gaussian計算得到的軌道的等值面圖，可以發現這就是想要研究的pi(in)→pi(in)*結合pi(out)→pi(out)*的激發，星號代表空軌道，電子躍遷時自旋都發生了beta→alpha的翻轉（后同）。算三重態情況略復雜。根據能量順序，最低三個三重態對應的電子激發的情況是這樣的：
T1：pi(in)→pi(out)*，簡稱T1(in-out)
T2：pi(out)→pi(out)*，簡稱T2(out-out)
T3：pi(in)→pi(in)*，簡稱T3(in-in)
? 直接把自旋多重度設為3做計算得到的是T1(in-out)。
? 想得到T3(in-in)，除了自旋多重度設3外，還需要用Gaussian的guess(read,alter)關鍵詞，讀取S0計算的chk文件當初猜，并在坐標末尾空一行寫49 50來交換初猜軌道的順序，以令默認處于占據的49號軌道pi(out)和默認處于非占據的50號軌道pi(in)軌道交換。
? 想得到T2(out-out)更為折騰，我實際發現若直接通過guess(read,alter)調換軌道順序做UKS計算試圖得到這個態，總收斂到T3(in-in)。最后解決辦法是先用ROKS結合軌道調換設置收斂到期望的T2(out-out)波函數，然后把fch文件用Multiwfn的主功能6的選項37將ROKS軌道split成為UKS波函數形式，之后導出fch并轉成chk。最后UKS計算時將自旋多重度設成3并結合guess=read讀取這個chk文件里的波函數當初猜即可。

T4態應當對應pi(out)→pi(in)*。實際發現由于變分坍塌問題難以得到這個態，并且可以預計這個態的芳香性和T1(in-out)并不會有顯著區別，因此文中就沒有研究T4。

2 電子結構和激發能

ΔSCF方式得到的上述激發態的激發能為
T1(in-out)：2.082 eV
T2(out-out)：2.106 eV
T3(in-in)：2.114 eV
Q1(in-in;out-out)：4.295 eV
可見三個最低三重態的能量差很小，本來也是因為pi(in)和pi(out)能級分布特征極為接近，這從筆者的Carbon, 165, 461 (2020)的圖2展示的MO-PDOS圖就可以清楚地看出來。Q1態由于涉及兩個電子的激發且分別來自pi(in)和pi(out)、彼此正交，因此激發能是T態的大約兩倍。

此文研究的不同電子態的軌道的能級，以及對應的S0極小點所屬的D8h點群下的不可約表示，如下所示。紅色和藍色分別對應pi(in)和pi(out)軌道

由上圖可以看出T和Q態存在顯著的自旋極化，使得alpha和beta電子的前線軌道能級明顯不匹配。并且根據不可約表示的匹配關系，可以看出激發特征是
T1(in-out): b1g→b1u
T2(out-out): b2u-b1u
T3(in-in): b1g→b2g
Q1(in-in;out-out)：pi(in) b1g→b2g、pi(out) b2u→b1u

3 不同電子態的原子受力

優化出來的C16基態結構的長、短C-C鍵的鍵長分別為1.363和1.213 ?，而C18為1.334和1.221 ?，前者的鍵長交替（bond length alternation, BLA）值0.150 ?明顯大于后者0.113 ?。芳香性傾向于令共軛環上的鍵長平均化，因此從幾何結構上雖然無法判斷C16的基態是否有明顯的反芳香性，但至少足矣看出C16的芳香性遠弱于C18。倘若C16的激發態的芳香性顯著強于其基態，那么在相應的Franck-Condon結構下，激發態的原子受力理應傾向于讓BLA減小。基于這個考慮，文中計算了C16的不同激發態的Franck-Condon結構下的受力，并借助VMD的作圖命令畫成了圖，如下所示，箭頭長度正比于原子受力大小

可見被激發到Q1態后原子受力顯著，極為傾向于讓較長的C-C鍵變得更短、讓較短的C-C鍵變得更長，從而減小BLA，這直接暗示出Q1具有很顯著的芳香性。T1/2/3態的原子受力情況極其接近、難以區分，故在上圖沒有區分給出，由圖可見雖然此時的受力也有令BLA減小的特征，但遠不如Q1態顯著，故從原子受力角度直接體現出T1態有比基態顯著得多的芳香性但又遠不如Q1態。此文的Franck-Condon結構的受力分析給激發態芳香性的分析提供了全新的視角。

4 不同電子態的NICS

NICS是最常用，而且很穩健、普適的衡量芳香性和反芳香性的定量指標，對基態和激發態都能用，因此本文也使用了NICS考察了不同激發態的芳香性。關于NICS及其各種變體的介紹詳見《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）、《使用Multiwfn計算FiPC-NICS芳香性指數》（http://www.shanxitv.org/724）的相關部分，其中通常來說最可靠的是NICS(1)ZZ。NICS(1)ZZ的計算結果如以下表格左半邊所示，比0負得越多說明芳香性越強、比0正得越多說明反芳香性越強。可見用GIAO和CSGT方式得到的總值差不多，而CSGT還可以分解成不同軌道的貢獻便于更好地了解本質，見《基于Gaussian的NMR=CSGT任務得到各個軌道對NICS貢獻的方法》（http://www.shanxitv.org/670）。本文將CSGT的NICS(1)ZZ分解成了pi(in)、pi(out)和其它部分。

從以上數據可見，C16的基態確實是顯著的反芳香性，因為S0的NICS(1)ZZ明顯為正，并且pi(in)和pi(out)對反芳香性的貢獻基本一樣。而Q1具有極強的芳香性，總NICS(1)ZZ達到了驚人的-94.4 ppm，基本上均等地來自pi(in)和pi(out)的貢獻。相比之下，Carbon, 165, 468 (2020)中報道的具有基態雙芳香性的C18的NICS(1)ZZ才只不過-23.7 ppm，被視為芳香性代表的基態苯分子的NICS(1)ZZ也只有-29.9 ppm。基態是反芳香性的C16在特定的激發態下竟有如此強的芳香性，實在令人吃驚！

C16的T1、T2、T3的芳香性相差不大，從NICS(1)ZZ的數值來看都算得上是很顯著芳香性。基于CSGT的NICS(1)ZZ分解，可見它們的芳香性的本質十分不一樣。T1(in-out)的芳香性同時來自于pi(in)和pi(out)電子，貢獻程度相仿佛。T2(out-out)的芳香性完全由pi(out)系統的電子貢獻，而未被激發的pi(in)電子則和S0態一樣產生反芳香性貢獻，并且由于pi(out)對芳香性的貢獻遠大于pi(in)貢獻的反芳香性，因此T2(out-out)整體是顯著芳香性的。T3(in-in)和T2(out-out)的情況類似，只不過pi(in)和pi(out)產生的作用幾乎正好反過來。

Baird規則是經典的預測處于S1和T1激發態的環狀體系的芳香性的規則，例如S0態是芳香性的苯在S1和T1態是反芳香性的。當前研究的C16的Q1態可以稱為雙重Baird芳香性，這是極為獨特的特征。T2(out-out)和T3(in-in)態都算是Baird芳香性與Huckel芳香性的混合，前者占主導。具有雙芳香性的T1(in-out)則不屬于Baird芳香性，因為沒有牽扯到單套pi系統內的電子激發，它的狀態類似于兩套pi系統都處于單自由基狀態。

為了更進一步確認UKS方式算的C16不同電子態的NICS(1)ZZ的結果的合理性，文中還用Dalton程序在CASSCF級別下做了NICS(1)ZZ計算，這算是對于可能存在較顯著靜態相關的體系算NICS(1)ZZ能用的最穩健的方法了（盡管未充分考慮動態相關）。計算使用CAS(12,12)活性空間，利用《利用Multiwfn令Dalton計算時使用其它程序產生的軌道作為初猜》（http://www.shanxitv.org/740）介紹的方法，將Gaussian做對稱破缺UHF計算產生的波函數經由Multiwfn轉化出的UNO軌道作為CASSCF的初猜軌道。計算出的結果是S0、T1、Q1態的NICS(1)ZZ分別為33.6、-29.1、-95.6 ppm，可見結論和UKS方式算的沒有什么區別，進一步確認了UKS計算的C16不同自旋態的NICS(1)ZZ已經很可靠了，沒什么可質疑的。

NICS_ZZ指數的數值依賴于計算位置的選取，這是它的一個不足，只考慮距離環中心1埃位置的NICS(1)ZZ對芳香性的展現有可能不夠全面。《使用Multiwfn繪制一維NICS曲線并通過積分衡量芳香性》（http://www.shanxitv.org/681）介紹了對NICS在垂直于環方向進行掃描的方法，下圖是C16的不同電子態的這種曲線圖，從這個圖上可以看出無論NICS_ZZ的計算位置取在哪，前面基于NICS(1)ZZ的分析結論都是不變的，更進一步鞏固了結論的可靠性。

對上面這種曲線進行積分還能得到∫NICS指數，能夠比NICS(1)ZZ更嚴格地衡量芳香性。此文對前述的電子態計算了∫NICS，對芳香性的判斷結論與NICS(1)ZZ完全一致，并且發現NICS(1)ZZ與∫NICS之間有非常理想的線性相關性。

5 AICD和ICSS分析

《使用AICD 2.0繪制磁感應電流圖》（http://www.shanxitv.org/294）介紹的AICD圖是通過磁場作用下的感生電流角度直觀展現芳香性特征的很好、很常用的方法，在筆者的大量的與18碳環相關的芳香性和衍生體系的研究文章中也都使用了此方法，如Carbon, 165, 468 (2020)、Inorg. Chem., 62, 19986 (2023)、Chem. Eur. J., 29, e202300348 (2023)、Chem. Eur. J., 29, e202300348 (2023)、ChemPhysChem, 25, e202400377 (2024)。此文對C16的前述電子態都繪制了AICD圖，并分解成了pi(in)和pi(out)電子貢獻，結果如下所示，外磁場與環平面垂直并且指向圖的外側。為了看得清楚，AICD等值面上的小箭頭的整體效果用手畫的大箭頭表示了出來。

可見相對于NICS(1)ZZ，上圖更直觀地展現了不同激發態的芳香性本質。例如Q1態的圖上確實可以清晰看出這個電子態下體系同時在環平面內和環平面外產生了很顯著的滿足左手定則的（順時針）感生電流，體現出了很鮮明的雙芳香性特征。而比如T2(out-out)，則可明確看出其pi(out)電子具有和Q1態一樣的很強的順時針感生電流、貢獻顯著的芳香性，而其pi(in)電子則具有和S0態一樣的一般強度的逆時針感生電流、貢獻一定反芳香性，從AICD(tot)圖整體來看順時針環電流更明顯，故是芳香性的。

《通過Multiwfn繪制等化學屏蔽表面(ICSS)研究芳香性》（http://www.shanxitv.org/216）介紹的ICSS_ZZ圖對芳香性的展現遠比NICS_ZZ及其一維掃描更全面直觀，但計算耗時也高得多。文中繪制了不同電子態的ICSS_ZZ等值面圖（等值面數值為± 12 ppm）和截面填色圖，如下所示，橙色和藍色分別對應于磁屏蔽和去屏蔽區域。可見S0態在環內是去屏蔽的，而環的外圍是磁屏蔽的，這是典型的反芳香性體系具有的特征。而Q1則在環內和環外分別展現了極強的磁屏蔽和去屏蔽，再次體現出其超強的芳香性。相比之下，T1、T2、T3態雖然也和Q1的特征一樣，能看出整體是具有明顯芳香性的，但在環外的去屏蔽區域的等值面比Q1窄得多，而且環內屏蔽區域的等值面在垂直于環方向的延展程度不及Q1、填色圖中環內紅色區域也遠不及Q1的鮮艷，因此這些三重態的芳香性都顯著弱于Q1。有趣的是，由于三個T態對芳香性、反芳香性有所貢獻的pi電子不一樣，導致它們的ICSS_ZZ等值面的具體形狀有可查覺的不同。例如T2(out-out)由分布在環的上、下方的pi(out)電子貢獻其芳香性部分，這使得其環內屏蔽區域的等值面在垂直于環方向上的延伸幅度明顯大于T3(in-in)。

此文還對ICSS_ZZ在環內的柱形區域進行了積分以便定量對比，積分區域如上圖(f)部分的灰色區域所示，數值也羅列在圖上了。可見不同態之間的ICSS_ZZ積分結果的差異和NICS(1)ZZ、∫NICS的非常類似，例如Q1的芳香性都是T態的三倍左右。

6 總結

本文介紹的Chem. Eur. J. (2025) DOI: 10.1002/chem.202404138一文通過嚴謹的量子化學計算、原子受力，以及各種基于磁屬性芳香性判斷方法，充分、嚴格、全面地揭示了C16的基態、最低三個三重態激發態、最低五重態激發態的芳香性特征，不同芳香性衡量方法的結論高度吻合，并且文中還對造成芳香性差異的電子結構本質進行了具體討論。由于C16具有兩套獨立的pi電子體系，導致其激發態芳香性特征與一般芳香性或反芳香性分子具有顯著差異。表面上看起來芳香性差不多的T1、T2、T3態，其芳香性的內在構成截然不同。而具有雙Baird芳香性的Q1態的芳香性則強得令人印象十分深刻，NICS(1)ZZ接近-100 ppm，遠遠強于苯分子的芳香性。希望這篇內容新穎的文章能夠令讀者對碳環體系和激發態芳香性產生興趣，并在研究芳香性方面獲得啟發。

文/Sobereva@北京科音  2025-Apr-5

1 前言

Dalton這個量子化學程序雖然用起來復雜，但有很多獨特的亮點，例如可以在SOPPA(CCSD)和CASSCF級別下計算各種磁屬性（筆者最近在Chem. Eur. J., e202404138 (2025)文中就用Dalton在CASSCF下算了NICS），可以用響應方式較準確地計算磷光壽命，可以算雙光子吸收（北京科音高級量子化學培訓班http://www.keinsci.com/KAQC里專門有一節講了)，可以用高級別耦合簇方法算激發態（見《使用Dalton通過CC3方法極高精度計算激發態》http://www.shanxitv.org/463），等等。可惜Dalton的SCF收斂方面做得一般，比不上Gaussian，而且單靠Dalton無法產生一些重要的軌道，比如UHF計算產生的UNO軌道經常被用于CASSCF的初猜，但Dalton甚至連UHF都做不了，因為Dalton里開殼層只支持RO形式。

為解決以上問題，從2025-Apr-5更新的Multiwfn開始，Multiwfn產生Dalton輸入文件的功能支持了將Multiwfn中當前的軌道波函數寫入到Dalton的輸入文件.dal里作為初猜波函數。因此借助Multiwfn，Dalton用戶可以用Gaussian、ORCA、GAMESS-US等諸多主流的量子化學程序產生的波函數當初猜，從而解決Dalton遇到SCF不收斂的問題。而且，Dalton用戶還可以直接用Multiwfn產生的UNO、定域化分子軌道（LMO，見http://www.shanxitv.org/380）等當CASSCF初猜。

下面就通過例子簡單演示一下具體實現方法，本文中涉及的文件都可以在http://www.shanxitv.org/attach/740/file.zip下載。Multiwfn可以在http://www.shanxitv.org/multiwfn免費下載，相關常識看《Multiwfn FAQ》（http://www.shanxitv.org/452）。如果讀者還不知道Dalton怎么安裝、運行、產生輸入文件的話，看《量子化學程序Dalton的編譯方法和運行方式簡介》（http://www.shanxitv.org/463）。本文用的Multiwfn是2025-Apr-5更新的版本，Gaussian是G16 B.01，Dalton是2022版。

如果本文介紹的Multiwfn的功能給你的研究帶來了幫助，請在論文中引用Multiwfn剛啟動時提示的Multiwfn程序的原文，這是對Multiwfn開發和維護最好的支持！

2 把Gaussian收斂的HF波函數當Dalton做CCSD(T)的初猜

這一節以CH3NH2為例演示把Gaussian收斂的HF/cc-pVTZ波函數當Dalton做CCSD(T)/cc-pVTZ的初猜。首先用Gaussian做一個普通的單點計算，輸入文件是本文文件包里的CH3NH2.gjf，內容如下所示。算完后把得到的CH3NH2.chk用formchk轉成CH3NH2.fch。

%chk=D:\CH3NH2.chk
# HF/cc-pVTZ int=NoBasisTransform
[空行]
test
[空行]
0 1
 C                  0.05159500    0.70381800    0.00000000
 H                  0.59439800    1.06209400    0.88166000
 H                  0.59439800    1.06209400   -0.88166000
 H                 -0.94293100    1.18451300    0.00000000
 N                  0.05159500   -0.76078400    0.00000000
 H                 -0.45830300   -1.10306000    0.81258800
 H                 -0.45830300   -1.10306000   -0.81258800
[空行]
[空行]

之后啟動Multiwfn，載入CH3NH2.fch，依次輸入
100  //其它功能（Part 1）
2  //產生輸入文件
19  //產生Dalton輸入文件
CCSDpT.dal  //要產生的輸入文件名
[回車]  //要產生的.mol文件名用默認的CH3NH2.mol
y  //由于當前內存里有基函數信息，Multiwfn問你是否把軌道展開系數寫入到.dal文件里
1  //使用Dalton標準格式（4F18.14）方式寫入。即四個值換一行，每個值占18列、小數位數14個

此時當前目錄下就有了CCSDpT.dal和CH3NH2.mol。把CH3NH2.mol里的基組名都替換為cc-pVTZ_emsl。目前CCSDpT.dal里的計算設置對應B3LYP算單點，自行把關鍵詞改為CCSD(T)計算用的，此時CCSDpT.dal的內容如下
**DALTON INPUT
.RUN WAVE FUNCTIONS
**WAVE FUNCTIONS
.CC
*CC INPUT
.CC(T)
*ORBITAL
.MOSTART     //代表從**MOLORB部分讀取軌道展開系數作為初猜
 FORM18
**MOLORB (Punched by Multiwfn)
 -0.00001750284370  0.00060444251600  0.00013849165900 -0.00001024308040
  0.00002442752580 -0.00001726827680 -0.00014237439200  0.00004058047020
...略
**END OF INPUT

用Dalton進行計算，可以看到SCF一輪就收斂了，達到了我們的目的：
@ *** DIIS converged in   1 iterations !
@     Converged SCF energy, gradient:    -95.252488149534    1.44D-07
    - total time used in SIRFCK :              0.00 seconds
之后CCSD(T)計算也順利跑完了。

一些注意事項和相關信息：
? 對于廣義或部分廣義收縮基組，都應當用int=NoBasisTransform關鍵詞，拿不準就始終帶上。詳見《將Gaussian等程序收斂的波函數作為ORCA的初猜波函數的方法》（http://www.shanxitv.org/517）里對此關鍵詞的解釋。
? 把Gaussian的cc-pVDZ、cc-pVTZ的波函數當Dalton的初猜時，Dalton要用cc-pVDZ_emsl、cc-pVTZ_emsl基組，而不要用不帶_emsl后綴的，否則還是需要迭代很多次才能收斂（甚至可能不收斂），但帶不帶_emsl的版本收斂后的結果是相同的。而如果是def、def2系列基組，沒什么特別要注意的，比如Gaussian里def-TZVP和def2-TZVP做的計算，在Dalton里基組名分別寫成Turbomole-TZVP、def2_tzvp即可，注意對大小寫有要求，必須對應Dalton目錄下的basis子目錄里的基組文件名。
? Dalton默認是基于球諧型基函數做的計算，因此如果在Gaussian里使用的是比如6-31G*等D殼層默認是笛卡爾型基函數的基組，應當寫上5d關鍵詞強行要求用球諧型，詳見《談談5d、6d型d殼層基函數與它們在Gaussian中的標識》（http://www.shanxitv.org/51）。雖然Dalton也可以用笛卡爾型基函數，但Multiwfn往Dalton輸入文件里寫入初猜軌道系數的功能不支持此情況。
? 只要一個量子化學程序能輸出Multiwfn可以識別的帶有基函數信息的波函數文件，并且用的是球諧型基函數，則波函數都可以通過以上方式作為Dalton的初猜。支持的文件詳見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。例如ORCA做計算產生的molden文件也可以載入Multiwfn后像上面那樣產生帶初猜波函數的Dalton輸入文件。

3 把Gaussian收斂的DFT波函數當Dalton做DFT的初猜

此例以順鉑體系為例，把Gaussian收斂的PBE0波函數當Dalton做PBE0計算的初猜，對Pt用Stuttgart小核贗勢，對配體用6-311G*。Gaussian輸入文件如下（本文文件包里的Ptcoord.gjf）：

%chk=D:\Ptcoord.chk
#P PBE1PBE/genecp int=NoBasisTransform
[空行]
b3lyp/def2TZVP opted    //眾所周知這是純給自己看的標題行，別問為什么寫成這個
[空行]
0 1
[坐標略...]
[空行]
N H Cl
6-311G*
****
Pt
SDD
****
[空行]
Pt
SDD
[空行]
[空行]

按照和上一節例子相同的過程，把chk轉成fch，再用Multiwfn產生.dal和.mol文件。.dal里的泛函名改成PBE0，然后設置.mol里的基組和贗勢成為下面這樣，當前對Pt用的贗勢和贗勢基組和Gaussian里的情況一樣（若不確定的話，可以Gaussian里用gfinput關鍵詞輸出定義，并和Dalton的基組、贗勢目錄下的相應文件對照）。

ATOMBASIS
test molecule
Generated by Multiwfn
Atomtypes=4 Angstrom Nosymmetry charge=0
Charge=1.0 Atoms=6 Basis=6-311G*
H1   -0.82596600    1.64390200    2.14978600
H2    0.00000000    2.40772900    0.93575300
H3    0.82596600    1.64390200    2.14978600
H4   -0.82596600   -1.64390200    2.14978600
H5   -0.00000000   -2.40772900    0.93575300
H6    0.82596600   -1.64390200    2.14978600
Charge=7.0 Atoms=2 Basis=6-311G*
N1    0.00000000    1.59755500    1.56108400
N2   -0.00000000   -1.59755500    1.56108400
Charge=17.0 Atoms=2 Basis=6-311G*
Cl1    0.00000000    1.70827400   -1.36819100
Cl2   -0.00000000   -1.70827400   -1.36819100
Charge=78.0 Atoms=1 Basis=stuttgart_rsc_1997_ecp ECP=stuttgart_rsc_1997_ecp
Pt1   -0.00000000   -0.00000000    0.18195700

用Dalton計算后，經過4輪SCF就輕易收斂了，結果為-1152.597056，和Gaussian給出的-1152.596983非常接近。而如果此例Gaussian和Dalton都用HF的話，會和上例一樣，結果精確相同且SCF一輪可收斂。

4 把UNO軌道當Dalton做CASSCF的初猜

之前在《CASSCF計算雙自由基以及雙自由基特征的計算》（http://www.shanxitv.org/264）里示例過用Gaussian基于UNO初猜軌道做CASSCF(2,2)計算的方法，不了解這種計算和UNO概念的話建議先看看。這一節演示使用Gaussian產生UHF波函數，讀入Multiwfn并產生UNO軌道，然后寫入.dal文件用于Dalton做CASSCF(2,2)計算的初猜軌道的流程。示例體系還是264號博文里的C4H8雙自由基，基組還是此文里的6-31G*。

首先用以下Gaussian輸入文件（本文文件包里的C4H8.gjf）對C4H8做UHF計算得到對稱破缺波函數。不了解為什么用guess=mix的話看《談談片段組合波函數與自旋極化單重態》（http://www.shanxitv.org/82）。

%chk=D:\C4H8.chk
# UHF/6-31G* guess=mix 5d
[空行]
ub3lyp/6-31g(d) opted
[空行]
0 1
C     -0.74742092    1.77656753    0.00000000
H     -0.62438907    2.32965189    0.92333358
H     -0.62438907    2.32965189   -0.92333358
C     -0.74742092    0.30933780    0.00000000
H     -1.24978876   -0.09122321    0.88294562
H     -1.24978876   -0.09122321   -0.88294562
C      0.74742092   -0.30933780    0.00000000
H      1.24978876    0.09122321   -0.88294562
H      1.24978876    0.09122321    0.88294562
C      0.74742092   -1.77656753    0.00000000
H      0.62438907   -2.32965189   -0.92333358
H      0.62438907   -2.32965189    0.92333358
[空行]
[空行]

計算完畢后把C4H8.chk轉成C4H8.fch，載入Multiwfn，然后依次輸入
200  //其它功能（Part 2）
16  //基于fch文件里的密度矩陣產生自然軌道。這個功能在《在Multiwfn中基于fch產生自然軌道的方法與激發態波函數、自旋自然軌道分析實例》（http://www.shanxitv.org/403）里有專門的介紹
SCF  //讀取fch里的SCF類型的密度矩陣
1  //產生空間自然軌道
此時從屏幕上看到如下信息，這便是UNO軌道的占據數。其中最高占據自然軌道（HONO）和最低非占據自然軌道（LUNO）的占據數分別是1.123225和0.876775，都偏離整數顯著，很滿足期望。

 Occupation numbers:
    2.000000    2.000000    1.999999    1.999999    1.999904    1.999898
    1.999881    1.999835    1.999309    1.999120    1.999004    1.998977
    1.996829    1.996536    1.996053    1.123225    0.876775    0.003947
    0.003464    0.003171    0.001023    0.000996    0.000880    0.000691
...略

接著輸入
n  //不把自然軌道導出成new.mwfn文件并重新載入，因為當前不需要如此（如果你想看一下當前的UNO軌道圖像的話，應選y，之后用Multiwfn主功能0觀看軌道）
0  //返回主菜單
100  //其它功能（Part 1）
2  //產生輸入文件
19  //產生Dalton輸入文件
CASSCF.dal  //要產生的輸入文件名
[回車]  //要產生的.mol文件名用默認的C4H8.mol
y  //把當前內存里的軌道（即UNO）的展開系數寫入到.dal文件里
1  //使用Dalton標準格式（4F18.14）方式寫入

當前C4H8.mol里的基組默認就是6-31G*，不用改。修改CASSCF.dal，把以下兩行
.DFT
 B3LYPg
替換為CASSCF(2,2)計算的設置：
.MCSCF
*CONFIGURATION INPUT
.CAS SPACE
2
.ELECTRONS
2
.INACTIVE
15

Dalton計算完畢后，活性空間內自然軌道占據數為1.2679和0.7320，能量為-156.024938 Hartree，和Gaussian在同級別做CASSCF精確一致。你若想看Dalton做CASSCF產生的自然軌道，把計算產生的.tar.gz文件里的molden.inp解壓出來載入Multiwfn，并效仿《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）的方式觀看即可。

5 關于線性依賴基函數的問題

較大基組（尤其是帶彌散函數的），計算原子分布較致密的體系很容易出現線性依賴問題，Gaussian和Dalton此時都會自動去除線性依賴基函數，導致最終求解出來的軌道數會少于基函數數目。這種情況下，Gaussian輸出信息中的NBsUse（實際用的基函數數目，等同于實際求解出來的軌道數）小于NBasis，Dalton輸出文件中的Total number of orbitals小于Number of basis functions。這給上文的傳軌道帶來了麻煩。有兩個解決方法：
(1)不去除線性依賴基函數：Gaussian計算時帶著IOp(3/32=2)關鍵詞、Dalton計算時在**WAVE FUNCTIONS里面加上以下內容，從而讓兩個程序都不去除線性依賴基函數
*ORBITAL
.AO DELETE
1E-100
這種做法的缺點是可能導致數值不穩定性。
(2)自動去除線性依賴基函數，但讓Dalton少讀取軌道：比如Gaussian計算顯示NBsUse=50、NBasis=54，因此自動去除了4個線性依賴基函數。因此當Multiwfn載入fch文件后，最后4個軌道的展開系數都會為0，導出到.dal文件里也是如此。此時Dalton若照常讀取全部軌道，計算一開始會報錯*** ORTHO-FATAL ERROR ***。解決方法是在*ORBITAL里寫上
.DELETE
4
這樣Dalton就不會讀取最后4個軌道了，計算可以正常進行。

6 .dal里的軌道展開系數存在星號而無法讀取的問題

有時候Multiwfn寫入到.dal里的軌道展開系數會存在一堆星號，導致Dalton讀取軌道時提示input conversion error錯誤。這不是Multiwfn的問題，而是Dalton自身設計不周。Dalton定義的**MOLORB字段中記錄展開系數的格式在前面說了，每個數值都是F18.14格式，因此當有的軌道展開系數特別大，超過了其整數位記錄上限，就會被記錄成一堆星號。這種問題在基組較彌散而又不讓自動去除線性依賴基函數的情況下出現概率大。這是為什么Multiwfn導出Dalton輸入文件的時候特意讓你選擇記錄的格式，前面的例子都是選1用4F18.14格式，還可以選2用4E20.12格式，這是使用科學計數法輸出，因此就不會出現數值過大記錄成星號的問題了。相應地，必須自己修改Dalton的讀取軌道展開系數的源代碼。對于Dalton 2022版，需要修改Dalton源代碼目錄下的DALTON/sirius/sirgp.F文件，把PFMT = '(4F18.14)'替換為PFMT = '(4E20.12)'，之后重新編譯即可。

7 其它

Multiwfn就相當一個大熔爐，能讀取諸多量子化學程序產生的波函數，又能導出波函數作為諸多量子化學程序的初猜，這在《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）里專門說了。對于傳軌道到ORCA的情況我之前還專門寫過《將Gaussian等程序收斂的波函數作為ORCA的初猜波函數的方法》（http://www.shanxitv.org/517）。

MOKIT也是很好的可以實現不同程序間傳軌道的工具，fch到.dal和.mol的轉換可以用里面的fch2dal子程序。相對于MOKIT，本文介紹的用Multiwfn傳軌道有幾個額外的好處：
(1)Multiwfn不需要安裝MOKIT所需的Python環境。尤其是對于Windows用戶，Multiwfn解壓即用明顯更為方便
(2)Multiwfn特意考慮了.dal里記錄軌道展開系數的格式的問題并給了解決辦法，上一節說了
(3)Multiwfn自身可以產生諸多類型的軌道（定域化分子軌道、空間/自旋自然軌道、NAdO、AdNDP、NTO、雙正交化軌道等等）并輸出作為Dalton的初猜
而fch2dal也有額外的優點，它會把fch里的基組、贗勢的具體定義直接寫到輸出的.mol里，因此如果Gaussian計算時用的基組、贗勢在Dalton里沒有自帶也能直接計算而不需要自定義。并且fch2dal支持用笛卡爾型基函數的情況。

文/Sobereva@北京科音

First release: 2025-Feb-18   Last update: 2025-Feb-19

0 前言

pi-pi相互作用（pi-pi interaction）是化學體系中很常見而且很重要的一種弱相互作用，本文全面談談這種相互作用的各方面特征，以令讀者對其有充分的認識、能在自己的研究中正確分析討論。我在互聯網上廣泛答疑時也時不時看到有關于pi-pi作用的非常錯誤的理解，比如有人被一些文章誤導，竟以為pi-pi作用來自于軌道相互作用或靜電相互作用，而且這種說法還廣泛以訛傳訛，筆者希望通過此文以正視聽。此外，氫鍵、范德華作用等概念在結構化學書里普遍都有，但pi-pi作用鮮有提及，我也推薦相關書籍作者和結構化學教師參考此文的內容將pi-pi作用納入書籍和教學。本文中大量涉及非常流行的波函數分析程序Multiwfn，相關信息見《Multiwfn FAQ》（http://www.shanxitv.org/452）。

我有不少研究文章都和pi-pi作用有密切聯系，是pi-pi作用研究的典型范例，非常推薦讀者閱讀并歡迎引用：
? Theoretical Insight into Complexation Between Cyclocarbons and C60 Fullerene, Chem. Eur. J., 30, e202402227 (2024)。在《全面揭示各種碳環與富勒烯之間獨特的pi-pi相互作用！》（http://www.shanxitv.org/727）中有詳細介紹，專門研究了不同尺寸的碳環和富勒烯之間的pi-pi作用
? Intermolecular interaction characteristics of the all-carboatomic ring, cyclo[18]carbon: Focusing on molecular adsorption and stacking, Carbon, 171, 514-523 (2021)。在《全面探究18碳環獨特的分子間相互作用與pi-pi堆積特征》（http://www.shanxitv.org/572）中有詳細介紹，其中涉及兩個18碳環之間的pi-pi作用
? Molecular assembly with a figure-of-eight nanohoop as a strategy for the collection and stabilization of cyclo[18]carbon, Phys. Chem. Chem. Phys., 25, 16707 (2023)。在《8字形雙環分子對18碳環的獨特吸附行為的量子化學、波函數分析與分子動力學研究》（http://www.shanxitv.org/674）中有詳細介紹，其中涉及18碳環與OPP雙環分子的pi-pi作用
? Comment on “18 and 12 – Member carbon rings (cyclo[n]carbons) – A density functional study”, Mat. Sci. Eng.: B, 273, 115425 (2021)。在《18碳環（cyclo[18]carbon）與石墨烯的相互作用：基于簇模型的研究一例》（http://www.shanxitv.org/615）中有詳細介紹，其中涉及18碳環與石墨烯的pi-pi作用

1 pi-pi作用的基本特征

pi-pi作用的定義和說法比較亂，這里給出我認為最嚴格準確的定義：“pi-pi作用是相距較近的兩個片段上彼此朝向相對的pi電子之間的獨特的色散吸引作用”。這里做幾個注釋，結合后文的實例讀者會了解得更充分：

(1)pi-pi作用可以是分子間的也可以是分子內的，滿足以上定義即可
(2)諸如苯分子里面的pi電子之間的作用不叫pi-pi作用，那屬于共享電子作用
(3)兩套pi電子的分布必須近乎彼此相對才可能算pi-pi作用。而諸如兩套pi電子近乎肩并肩挨著就不能算pi-pi作用，像是不能說環丁二烯里面兩套近乎定域的pi電子之間是分子內pi-pi作用
(4)“相距較近”一般是在4.0-4.5埃以內。也不是說更遠距離就完全沒有pi-pi作用了，只不過由于色散作用隨作用距離r呈-1/r^6快速衰減行為，因此距離稍微一遠pi-pi作用就非常弱了，就不太值得一提了。但相互作用的片段間距離也不能太近，若顯著小于相接觸的原子的范德華半徑之和，則顯著的位阻互斥作用會遠大于pi-pi吸引作用，使得pi-pi作用也相對不值得一提
(5)pi-pi作用最常出現在碳原子間，因為碳最容易帶顯著的pi電子。碳的Bondi和CSD范德華半徑分別為1.70和1.77埃。如果沒有特殊的因素影響相互作用距離的話，在平衡結構（勢能面極小點結構）下，C-C間的pi-pi作用出現在3.4-3.6埃左右是最常見的
(6)pi-pi作用屬于弱相互作用和非共價相互作用范疇。雖然名義上算作弱相互作用，但實際強度可大可小，直接取決于作用面積，詳見后文
(7)有一個詞叫pi-pi堆積（pi-pi stacking），這個詞往往和pi-pi作用混用。但在我來看，這個詞更適合用來形容由于pi-pi作用的吸引效果使得相互作用的兩個部分發生緊密結合從而產生的彼此堆積的結構特征

下圖就是一個再典型不過的存在顯著pi-pi作用的體系，是暈苯的二聚體。高精度理論計算的極小點結構下平面間距大約在3.45埃

苯二聚體有不同構型，彼此垂直的T型二聚體是能量最低構型，還有一種能量與之非常接近的局部極小點構型是平行位錯構型，屬于pi-pi堆積結構。按照《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）的做法繪制出的幾何結構結合pi電子等值面圖如下所示，兩個分子的pi電子的分布非常清楚直觀，可以看到彼此相對。注意并非必須兩個分子的所有pi電子都彼此對著才算pi-pi作用，像下圖這樣哪怕只要一部分對著也同樣算pi-pi作用。而如果完全沒有對著、完全錯開了，那就不算pi-pi作用了。更多的展現pi電子分布的圖見筆者的Theor. Chem. Acc., 139, 25 (2020)一文。

2 pi-pi作用的物理本質

我在《談談“計算時是否需要加DFT-D3色散校正？” 》（http://www.shanxitv.org/413）提到過，原子間相互作用從物理本質上可以基本分為靜電、色散、交換、Pauli互斥、軌道相互作用；我在《談談“計算時是否需要加DFT-D3色散校正？”》（http://www.shanxitv.org/413）里也說過物理本質和相互作用表象之間的關系。pi-pi作用類似于氫鍵、鹵鍵等，是屬于表象層面的，用來描述一種具有特定特征的相互作用，其內在本質就是前面說的來自于相距較近的pi電子之間的色散作用，這一點在Grimme的一篇經典的討論pi-pi作用本質的文章Angew. Chem. Int. Ed., 47, 3430 (2008)中通過嚴格的理論計算和分析對比已經論證得十分嚴格充分，而且也早已是量子化學領域內行人的共識。

有的文章試圖從軌道相互作用解釋pi-pi堆積出現的本質，這是極其有誤導性的，卻還廣泛以訛傳訛。2007年的時候wiki上的詞條說pi-pi作用來自于pi共軛體系間p軌道的重疊，前述Angew文章的腳注中專門對此進行了批判，指出沒有任何量子化學計算能支持這種說法。這也不難理解，本來兩個分子的pi軌道之間重疊程度極低（可以用《分子間軌道重疊的圖形顯示和計算》http://www.shanxitv.org/163介紹的方法考察），遠低于形成共價鍵時原子軌道間的重疊程度，因此有效重疊不足。而且兩個分子的pi占據軌道之間混合會產生全占據的成鍵和反鍵軌道，二者效果充分抵消了，對成鍵沒貢獻；而若一個分子的pi占據軌道和另一個分子的pi非占據軌道混合，意味著要出現分子間pi->pi*電子轉移，這通常又沒有明顯能使體系能量變得更低的內在驅動力（不像常見的較近距離的n->pi*超共軛往往可以明顯降低體系能量）。此外，按照《Multiwfn支持的分析化學鍵的方法一覽》（http://www.shanxitv.org/471）說的，用Multiwfn對pi-pi堆積的兩個分子間去計算衡量兩個部分之間有效共享電子對數的Mayer鍵級、模糊鍵級或離域化指數（DI），會發現數值很接近0，比如《全面揭示各種碳環與富勒烯之間獨特的pi-pi相互作用！》（http://www.shanxitv.org/727）介紹的文章里對18碳環與富勒烯之間算的Mayer鍵級僅為可忽略不計的0.04，也充分體現出pi-pi作用中的軌道相互作用成份基本可忽略不計。所以無論從哪個角度，pi-pi堆積的出現都不是軌道相互作用驅動的，頂多是在色散作用驅動下，pi-pi堆積結構出現后順帶出現了點可忽略的軌道相互作用而已。所以pi-pi作用在主流學術界被歸為非共價相互作用范疇，這是完全正確的。然而有一些論文試圖從軌道相互作用分析pi-pi作用的特征，還擺出一堆軌道圖，比如Phys. Chem. Chem. Phys., 15, 9397 (2013)就是典型文章，切勿把這類文章太當回事，很多都是瞎討論、尬討論。

pi-pi堆積的體系有一個普遍特點是在極小點結構下，pi-pi作用區域的原子間通常不是正好對著，而是相互錯位。下圖是暈苯二聚體的俯視圖，其中一個暈苯用紅色顯示，可以很清楚看到錯位特征，即一層的碳對著另一層的六元環中心。一層層堆疊的石墨中每一層之間也同樣是這樣錯位的。出現自發錯位有不同的解釋并且存在一定爭議，有人認為是因為錯位的結構下比嚴格面對面的結構下的原子間的Pauli互斥更低，例如Chem. Sci., 11, 6758 (2020)持這種觀點。也有人認為是因為錯位的結構下的靜電相互作用能更負（靜電吸引作用更強），在Phys. Chem. Chem. Phys., 24, 8979 (2022)中作者通過能量分解論證了這一點并反駁了Chem. Sci.那篇文章。靜電作用引起位錯的原因從邏輯上容易理解：形成pi-pi堆積的兩部分都有豐富的pi電子，面對面構型下pi電子間靜電互斥會較大，錯開后互斥自然就沒那么強了。雖然筆者更同意靜電效應是位錯出現的主導，但我也不否認降低Pauli互斥也可能是產生位錯的另一個驅動力，盡管相對次要。

順帶一提，不止是pi-pi堆積二聚體，還有很多其它體系都是色散吸引作用驅動分子間結合，而靜電作用進一步影響幾何結構。《靜電效應主導了氫氣、氮氣二聚體的構型》（http://www.shanxitv.org/209）介紹的筆者的研究就是很好的例子，H2、N2的各種二聚體的相互作用能通過能量分解分析可以發現都是色散作用為主，但不同構型下的靜電作用的差異則直接影響不同構型的穩定性乃至存在與否，并且筆者發現通過靜電勢互補原理可以給出很直觀的解釋。這種靜電勢互補的分析方法也在《全面探究18碳環獨特的分子間相互作用與pi-pi堆積特征》（http://www.shanxitv.org/572）介紹的筆者的論文中也用來解釋為什么18碳環的pi-pi堆積二聚體是錯位的。平面環狀結構的18碳環具有獨特的平面內和平面外兩套pi電子，此文中確認了18碳環能夠形成分子間彼此平行的pi-pi堆積二聚體，下圖是按照《使用Multiwfn+VMD快速地繪制靜電勢著色的分子范德華表面圖和分子間穿透圖》（http://www.shanxitv.org/443）繪制的二聚體極小點結構下的單體的靜電勢著色的分子范德華表面的疊加圖，可以看到兩個碳環表面靜電勢為正和為負的區域在錯位的結構下是以互補的方式重疊的，很清晰地解釋了錯位的產生原因，即這種構型下的靜電吸引作用比一個個碳原子彼此精確對著的時候更強。

要注意的是，有些文章作者極度夸大了靜電作用對于pi-pi堆積形成的作用，不僅認為平行位錯構型的出現是因為此時靜電吸引作用最有利，還鼓吹pi-pi堆積結構的形成的本質就是靜電吸引作用，這是大錯特錯！最典型的這種文章就是Chem. Sci.,3 , 2191 (2012)，這篇文章的發表嚴重誤導了很多讀者。此文雖然是在前述的Grimme的pi-pi作用的研究之后寫的，而且文中還引了那篇文章，但給我的感覺是此文的作者完全沒好好看那篇文章，似乎也完全不懂什么叫電子相關作用，而依然基于古老且過于簡單（忽略了pi電子結構和穿透效應）的Hunter和Sanders的四極矩觀點盲目、主觀、武斷、信誓旦旦地解釋pi-pi堆積的位錯結構形成的本質，而完全沒有可信、嚴格的理論和計算數據作為依據，各種自說自話，并且還錯誤地解釋一些文獻里的理論計算數據，甚至還不承認pi-pi作用的概念。作者還執拗地認為非得是存在精確面對面堆積才能說明pi-pi作用的存在，并因為實際pi-pi堆積體系的極小點結構不是面對面堆積就被他用來否認存在pi-pi作用。此文文中給出了下圖，令一些初學者不明覺厲，誤以為不僅清晰解釋了為什么平行位錯結構比精確面對面結構更有利，還誤以為這種二聚體的形成主要靠的就是靜電吸引作用，完全沒分清楚pi-pi堆積二聚體的形成中色散作用和靜電作用的主次關系。還有不少其它文章如J. Chem. Soc., Dalton Trans., 2000, 3885也是類似地對pi-pi作用存在嚴重誤解、給出誤導性的示意圖，讀者看到這種文章時切勿當回事！

能量分解分析在pi-pi作用的研究中應用得非常普遍，對于認識其本質非常有幫助。這類方法將總的相互作用能分解成不同物理成份。例如將《使用sobEDA和sobEDAw方法做非常準確、快速、方便、普適的能量分解分析》（http://www.shanxitv.org/685）介紹的筆者提出的流行的sobEDAw能量分解方法用于苯的兩種二聚體構型可得到如下結果（計算級別：B3LYP-D3(BJ)/6-311+G(2d,p)并考慮counterpoise，單位為kcal/mol），可見色散項ΔE_disp遠比靜電作用項ΔE_els更負，軌道相互作用項ΔE_orb更是微乎其微，故色散作用對苯的分子間結合起到主導作用。特別是在平行位錯的結構下的ΔE_disp比T構型下的明顯負得多（同時ΔE_els的大小還變小），充分體現出平行位錯結構下才有的pi-pi作用的本質是色散作用，這和上文對pi-pi作用本質的討論相一致。平行位錯結構下，色散作用對總吸引作用項的貢獻百分比是：-7.6/(-7.6-0.8-1.6)*100=76%。

順帶一提，pi-pi堆積的結構在現實環境中一般很容易發生分子間相對滑移（除非有額外的位阻效應等阻礙），這是由于滑移導致能量變化很小。Carbon, 171, 514-523 (2021)文中做18碳環二聚體的勢能面掃描、Phys. Chem. Chem. Phys., 24, 8979 (2022)中做的苯分子平行二聚體的勢能面掃描都體現了這一點。并且由于滑移方向的勢能面很緩，哪怕室溫下也非常容易出現滑移。在Carbon, 171, 514-523 (2021)文中我給出了18碳環二聚體在100K下的4 ps的從頭算動力學的軌跡疊加圖（消除了下方的18碳環的整體運動以著重表現相對運動），如下所示，分子位置根據模擬時間按照藍-白-紅變化，可見在模擬過程中的確出現了顯著的滑移

3 pi-pi作用的強度

pi-pi堆積的極小點結構下的每一對近距離接觸的原子的pi-pi作用都是非常弱的，就是普通范德華作用的程度，但是當涉及pi-pi作用的原子很多時，就可以達到很高的強度，甚至能達到近乎化學鍵的程度。最強的pi-pi作用就是從微觀來看無限延展的石墨中的層間pi-pi作用了。下面看一些實例，以令讀者更好地理解pi-pi作用的強度特征。

Angew. Chem. Int. Ed., 47, 3430 (2008)文中計算了T型苯二聚體（下圖a）、平行位錯苯二聚體（下圖b）、環己烷二聚體（下圖c和d是兩個視角）的相互作用能，并且還將它們的環數n從1逐漸加到4（分別對應下圖e、f、g的結構），以考察了三種作用類型的相互作用能隨環數的變化，分別對應下圖的aromatic T-shaped、aromatic stacked和saturated stacked三條曲線。由圖可見，當pi-pi作用涉及的原子很少時，即苯二聚體，由于pi-pi作用還不夠強，因此平行位錯二聚體的能量比靜電吸引作用更強的苯T型二聚體略微更高，相互作用強度比起無pi-pi作用特征的環己烷二聚體沒優勢。但隨著環數增加，pi-pi作用強度不斷增大，使得具有pi-pi堆積結構的二聚體的相互作用變得顯著強于不具備這種特征的二聚體。四并苯的二聚體的相互作用能已達到-16 kcal/mol（67 kJ/mol），這已經是水二聚體這種普通強度氫鍵作用能的3倍左右了。

可見，前面提到的Chem. Sci.,3 , 2191 (2012)那篇文章拿兩個苯及其衍生物在液體和晶體中普遍缺少平行位錯構型的出現，以及平行位錯的苯二聚體不是能量最低構型，以此鼓吹pi-pi作用不是一種客觀存在的作用，這明顯是極其狹隘的。

對本文第1節的兩個暈苯的二聚體，在DLPNO-CCSD(T)級別下算出來的相互作用能達到-20 kcal/mol，也是相當強了。

前面給出的18碳環二聚體的相互作用能，在Carbon, 171, 514-523 (2021)中我用ωB97X-V/def2-QZVPP結合counterpoise校正的計算結果為-9.2 kcal/mol（-38.5 kJ/mol），達到平行位錯的苯二聚體的相互作用能的大約三倍，無疑也是非常顯著的pi-pi作用了。

不是只有平面區域之間才可能有pi-pi作用，例如下面三個體系的pi-pi作用區域都是有顯著曲率的。下圖左邊的結構是Org. Lett., 17, 5292 (2015)中合成出的主-客體復合物晶體的一部分，由于客體分子與富勒烯之間的pi-pi作用面積頗大，B97D/QZVP級別計算的相互作用能達到了-50 kcal/mol，已經是弱化學鍵的強度了。下圖右側是碗烯的pi-pi堆積二聚體，相互作用的研究見J. Phys. Chem. A, 116, 11920 (2012)。

4 圖形化展現pi-pi作用區域

《使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用》（http://www.shanxitv.org/621）以及《一篇最全面介紹各種弱相互作用可視化分析方法的文章已發表！》（http://www.shanxitv.org/667）中的綜述詳細介紹的筆者提出的IGMH方法用于圖形化展現片段間相互作用效果極佳，已被廣為使用。IGMH用于展現自定義片段間的pi-pi作用、判斷pi-pi作用是否存在尤為好用，因此在此給出一些簡單例子予以體現。其方法的原理、細節、計算操作見上面的文章。如果你不知道pi-pi作用可能出現在哪，不方便定義片段，則應當用《使用IRI方法圖形化考察化學體系中的化學鍵和弱相互作用》（http://www.shanxitv.org/598）介紹的筆者的IRI方法，可以把體系中所有相互作用區域全都展現出來，也包括化學鍵作用區域。

不是只有純碳的pi區域之間才可能有pi-pi作用，pi-pi作用也可以涉及到其它原子。眾所周知，DNA結構中相鄰層的堿基與堿基之間是平行堆疊的，再加上堿基區域有大量pi電子，所以這也屬于典型的pi-pi作用。下圖是從DNA中截取的堆疊的GCGC堿基對，兩層分子各定義為一個做IGMH分析的片段，兩層之間的綠色的IGMH方法定義的delta_g_inter函數的0.004 a.u.等值面非常清晰、優雅地展現出了pi-pi作用的主體區域。圖中還有個面積較小的等值面出現在O...O之間，雖然C=O鍵的O有pi電子，但這倆O的pi電子不是彼此對著的，所以這就只能算是普通的范德華作用區域了。

下圖是《8字形雙環分子對18碳環的獨特吸附行為的量子化學、波函數分析與分子動力學研究》（http://www.shanxitv.org/674）介紹的筆者的論文中的一張圖，研究的是OPP雙環分子結合兩個18碳環后的復合物。圖中綠色等值面直觀地展現出在什么區域18碳環與OPP主體分子間有顯著的pi-pi作用。可見pi-pi作用主要出現在OPP和18碳環近距離接觸的體系的兩端，而OPP靠中間區域離18碳環較遠因此缺乏pi-pi作用。順帶一提，18碳環是極少有的同時具有平面內（in-plane）和平面外（out-of-plane）兩套全局離域的pi電子的體系，當前體系中18碳環主要憑借的是平行于環方向分布的pi電子與OPP大環上的垂直于苯環的pi電子產生pi-pi作用，這是極為新穎、獨特的pi-pi作用形式。如果讀者對碳環類體系感興趣，非常建議進入http://www.shanxitv.org/carbon_ring.html觀看筆者發表的所有碳環相關的理論研究和對應的深入淺出的評述文章。

上圖中還對利用Multiwfn的基于力場的能量分解（EDA-FF）功能計算了各個原子對分子間相互作用的貢獻并對原子進行了著色，展現出了更豐富的信息。這種分析見《使用Multiwfn做基于分子力場的能量分解分析》（http://www.shanxitv.org/442）。

在《全面揭示各種碳環與富勒烯之間獨特的pi-pi相互作用！》（http://www.shanxitv.org/727）介紹的文章中，筆者全面研究了碳環與富勒烯的相互作用，其中給出了下圖，展示了不同尺寸的兩個碳環與富勒烯形成的三聚體。并且為了直觀區分富勒烯-碳環和碳環-碳環兩種pi-pi作用，分別將對應的等值面用綠色和青色著色。可見此圖極好地將體系中兩種pi-pi作用出現的主要位置展現了出來，比起在文中只是給出個幾何結構圖信息豐富得多。

標準的IGMH圖是通過sign(λ2)ρ函數對等值面著色的，ρ是電子密度。按照IGMH方法的標準的色彩刻度著色的話，ρ接近0的區域的等值面都是綠色。由于pi-pi作用區域、普通范德華作用、極弱的氫鍵（如C-H作為氫鍵給體時）等情況的作用區域的電子密度都很低（0.01 a.u.以內），因此相應的標準方式著色的IGMH等值面都會是綠色或很接近綠色。因此解釋IGMH圖的等值面時必須結合不同特征的弱相互作用的基本定義，不能單純根據顏色去亂解釋。比如網上有人問我“下面這個圖里紅圈部分是pi-pi作用嗎？”，這實在是匪夷所思的問題。跨越那個等值面和右下方苯環的pi電子區域直接作用的是一個氫原子，氫原子又沒pi電子，怎么可能是pi-pi作用！？而且右下方的苯環也明顯不可能和左上方的苯環有pi-pi作用，一方面二者離得老遠，另一方面二者的pi電子區域也根本不相互對著，顯然無論如何也不可能算是pi-pi作用。

那么上圖中紅圈里的等值面算什么作用？明顯這是pi-氫鍵，C-H作為氫鍵給體，其中氫帶一定正電荷（這一點用H的原子電荷就能體現，原子電荷的相關知識見《一篇深入淺出、完整全面介紹原子電荷的綜述文章已發表！》http://www.shanxitv.org/714這篇綜述），右下角的苯環上的pi電子令苯環上方顯負電性而作為氫鍵受體（用靜電勢可直接體現，參見http://bbs.keinsci.com/thread-219-1-1.html里匯總的靜電勢相關資料和我的相關博文）。

下面再舉一個IRI圖展現pi-pi作用區域的例子。Multiwfn的原文之一J. Chem. Phys., 161, 082503 (2024)里給出了144-輪烯的LOL-pi函數的等值面圖，如下圖左側所示，此圖清晰直觀地展現出了這個全局大共軛體系的pi電子的主要離域路徑。根據此圖展現的pi電子的分布情況，憑直覺就知道這個體系里肯定存在pi-pi作用。下圖右側是對這個體系的其中一個局部區域繪制的IRI等值面圖，藍色的等值面展現出了化學鍵作用區域，綠色的等值面清楚直觀地展現出了pi-pi作用區域。可見這個體系非常有趣，pi-pi作用區域綿延不斷貫通整個體系！也正是有分子內pi-pi作用的存在，此體系才能形成螺旋狀結構，要不然就散了（如同用不能描述色散作用的理論方法優化DNA結構時結構會散掉）。

5 衡量pi-pi作用強度的方法

這一節說一下如何衡量pi-pi作用強度。最簡單的方法就是計算pi-pi堆積二聚體的相互作用能，作用能越負說明作用越強。但如果要討論二聚體的熱力學穩定性，則需要計算結合自由能，溶劑環境中還得考慮溶劑模型。相關知識參見《談談分子間結合能的構成以及分解分析思想》（http://www.shanxitv.org/733）。復合物AB在特定結構下，A和B的相互作用能的最常規的計算方法就是用E(AB)-E(A)-E(B)方式手動計算，做sobEDAw能量分解（http://www.shanxitv.org/685）時也會順帶給你相互作用能。

以上述方法算相互作用能的一個問題是，如果兩個分子之間不僅僅有pi-pi作用，還有其它作用（如氫鍵），那么得到的只是總相互作用能。如果想只得到pi-pi作用能，有幾個辦法可以用：
(1)用《使用Multiwfn做基于分子力場的能量分解分析》（http://www.shanxitv.org/442）介紹的EDA-FF能量分解方法，將兩個分子的pi作用區域定義為兩個片段，讓Multiwfn給出基于分子力場的這兩個部分的相互作用能，取其中的范德華作用能（即色散作用和交換-互斥作用之和）。如果你只需要色散部分，還有另一種做法，見《使用Multiwfn圖形化展現原子對色散能的貢獻以及色散密度》（http://www.shanxitv.org/705）。這兩種做法還都可以給出具體原子產生的貢獻，并可以對原子進行著色以便通過圖像直觀展現和考察
(2)如果分子間同時有pi-pi作用和氫鍵，且其它作用可忽略不計，可以按《透徹認識氫鍵本質、簡單可靠地估計氫鍵強度：一篇2019年JCC上的重要研究文章介紹》（http://www.shanxitv.org/513）介紹的方法使用Multiwfn做拓撲分析估計出氫鍵作用能，從總相互作用能中扣掉之作為pi-pi作用能
(3)對體系進行恰當改造，基本保留每個分子參與pi-pi作用的部分，然后將總相互作用能近似當成pi-pi作用能。

如果是分子內的pi-pi作用能的估計，還可以參考《計算分子內氫鍵鍵能的幾種方法》（http://www.shanxitv.org/522）里的說明舉一反三處理。

在特定條件下，pi-pi作用強度和IGMH的等值面面積有密切的正相關性。《全面揭示各種碳環與富勒烯之間獨特的pi-pi相互作用！》（http://www.shanxitv.org/727）介紹的文章中給出了下圖，是不同原子數的碳環與富勒烯之間的相互作用能和描述pi-pi作用的IGMH等值面面積的關系，可見有挺理想的線性關系。通過這種關系，可以直接根據pi-pi作用區域的等值面面積估計對應的pi-pi相互作用能。這種面積的計算方法見《計算IGMH等值面的面積和體積的方法》（http://www.shanxitv.org/738）。

Chem. Commun., 48, 9239 (2012)提出的LOLIPOP方法值得一提。基于單個分子的波函數文件，就可以用Multiwfn非常容易地計算體系中的各個六元環的LOLIPOP指數，此值越小說明這個環發生pi-pi堆積的能力越強，原文對不同體系以苯分子作為探針分子進行了測試證實了這一點。Multiwfn手冊3.100.14節有LOLIPOP的詳細介紹，4.100.14節有計算實例。

6 不同理論方法描述pi-pi作用的能力

色散作用的更深層本質是電子的庫侖相關作用，因此做量子化學或第一性原理計算時，若想描述好pi-pi作用，用的方法必須能描述好庫侖相關，顯然高精度后HF方法都沒問題，比如CCSD(T)在計算pi-pi作用上可以算是金標準。至于低級別的后HF方法MP2則傾向于明顯高估pi-pi作用，絕對不要用。

對于特別常用的DFT，描述pi-pi作用的能力基本等同于描述普通色散作用的能力，如果泛函原本描述色散作用爛（如PBE、PBE0）或者完全不能描述（如B3LYP），就必須帶色散校正，參考《談談“計算時是否需要加DFT-D3色散校正？” 》（http://www.shanxitv.org/413）、《DFT-D色散校正的使用》（http://www.shanxitv.org/210）、《談談量子化學研究中什么時候用B3LYP泛函優化幾何結構是適當的》（http://www.shanxitv.org/557）。本來就帶色散校正的wB97M-V、wB97X-D3等直接就能很好描述pi-pi作用。M06-2X描述pi-pi作用雖然定性正確但不算太好，加了零阻尼DFT-D3色散校正后對pi-pi作用的描述有明顯改進，但還是不如B3LYP-D3(BJ)，對比測試見考慮了很多pi-pi作用體系的L7測試集的原文（J. Chem. Theory Comput., 9, 3364 (2013)）。雙雜化泛函由于帶有MP2項，所以都有描述pi-pi作用的能力，但一般在考慮了色散校正后才會變得足夠好，如revDSD-PBEP86-D3(BJ)。

實際上DFT-D那種形式的色散校正在原理上對于描述pi-pi作用并非很理想，因為pi電子不是繞著原子核球對稱分布的，而DFT-D校正能的公式依賴的只是原子間距離（這里不考慮三體校正項），沒體現出pi電子在原子核周圍的具體分布特征。不過這倒也不是明顯問題，常用的DFT-D3、DFT-D4色散校正對于描述pi-pi作用從實際效果上來看并沒有什么問題。

在半經驗方法層面，專門考慮了對色散作用描述的GFN2-xTB、PM6-D3H4X'、PM6-ML在表現pi-pi作用方面優秀，見J. Chem. Theory Comput., 21, 678 (2025)里3圖基于L7測試集的測試，PM6-D3和PM7只能算是定性正確。至于沒專門考慮對色散作用描述的諸如PM6、AM1等方法則完全失敗。

主流的分子力場，如GAFF、AMBER、CHARMM、OPLS-AA、MMFF94等，對pi-pi作用的描述雖然跟像樣的量子化學方法比算不上出色，但至少也算定性正確，因為它們都有描述色散作用的能力。J. Chem. Inf. Model., 49, 944 (2009)的測試專門體現了這點（不過這篇文章也有不少漏洞和槽點）。

7 疏水作用與pi-pi作用的關系

眾所周知，疏水效應使得水環境下非極性物質傾向于發生聚集，本質是溶劑的熵效應。兩個石墨烯片段在水中會自發堆積在一起，這算疏水作用還是pi-pi作用？實際上二者都有，對于堆積結構的形成起到協同作用。疏水效應更為普遍，無論兩個溶質的接觸區域是否有pi電子，只要溶質是基本無極性的，在水中都有疏水作用促使它們發生結合。而對于pi電子區域暴露的兩個溶質，疏水效應則在pi-pi作用的基礎上進一步促進了它們的pi-pi堆積結構的出現。如《全面揭示各種碳環與富勒烯之間獨特的pi-pi相互作用！》（http://www.shanxitv.org/727）介紹的文章的理論計算所示，在水環境下碳環與富勒烯之間的結合自由能的大小顯著大于在真空下，充分體現了這一點。

順帶一提，前述的Chem. Sci.,3 , 2191 (2012)一文居然誤以為溶劑環境下pi共軛體系間出現堆積結構僅僅是因為溶劑效應，并似乎試圖靠這個否認真空環境下也存在pi-pi相互作用，甚至說the terms "pi-stacking" or "pi–pi interactions" do not describe any physically meaningful interaction，實在是難以理喻！！！PS：這樣的文章能通過Chem. Sci.的審稿真是離譜。

8 總結&判斷pi-pi作用的標準

本文系統地對pi-pi作用的各個方面進行了介紹，包括其基本特征、物理本質、強度范圍、圖形化展現方法、考察強度的方法、理論計算方法的精度、疏水作用與它的關系。通過本文，讀者應該已經對pi-pi作用有了較全面的了解，能夠進行正確的分析討論，并且能認識到哪些文章或書籍里的說法是有誤導性的。我也推薦讀者接下來閱讀本文一開始提到的筆者的一系列和pi-pi作用有關的研究文章和對應的介紹博文。

最后再總結一下pi-pi作用的常規判斷標準，便于讀者能確切判斷哪些作用算是pi-pi作用。通過下面第1、2條就可以進行粗略判斷，3、4、5可以作為進一步檢驗。一般意義的pi-pi作用應當能同時滿足所有條件
(1)相互作用的兩部分都有pi電子且其分布彼此相互對著。如果拿不準有沒有pi電子分布、分布朝向如何，可以按照《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）介紹的方法直接把pi電子密度等值面圖畫出來，一目了然
(2)相互作用的pi電子之間離得不太遠。比如明顯超過5埃的直接就可以忽略了
(3)用Multiwfn做IGMH分析（如果得不到波函數文件或體系太大難以算得動的話，可以改用極便宜且只依賴于原子坐標的mIGM），在等值面數值調到諸如0.003 a.u.這樣較小的值時，在預期出現pi-pi作用的區域能明顯看到等值面
(4)兩部分之間的Mayer鍵級或模糊鍵級或離域化指數非常小（遠小于0.1），故而軌道相互作用可基本忽略。注：如《Multiwfn支持的分析化學鍵的方法一覽》（http://www.shanxitv.org/471）所述，Multiwfn在主功能9里計算鍵級之前可以用選項-1定義兩個片段，之后做鍵級計算時會給出兩個部分之間的總鍵級，即片段間每一對原子的鍵級的總和
(5)使用sobEDAw等能量分解方法分析，色散項應當占所有吸引項總和的大部分，軌道相互作用項應當占比微小

文/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

文/Sobereva@北京科音

 First release: 2024-Dec-31   Last update: 2025-Mar-28

在我網上回答大量量子化學的問題的過程中，經常看到有些人搞不清楚對分子間相互作用有貢獻的量之間的關系、不知道怎么將結合能分解成不同部分考察影響結果的內在因素。這使我決定作一張圖，把分子間結合自由能的構成一次性梳理清楚，同時展現我的對分子間結合能的分解分析的思想，如下所示：

下面我會把圖中各項進行簡要的解釋，便于讀者理解每項是怎么回事、知道怎么計算，務必結合上圖看。上圖的源頭是“溶劑中的結合自由能”，因為這個是包含的成份最多，內涵最復雜的量。本文僅考慮兩個單體分子A、B間結合成為復合物AB的情況。

? 溶劑中的結合自由能ΔG_bind(sol)：這是指在溶劑（sol）環境下，兩個分子從各自孤立狀態相互接近并形成二聚體過程中的自由能的總變化。即ΔG_bind(sol) = G_AB(sol) - G_A(sol) - G_B(sol)。如《談談隱式溶劑模型下溶解自由能和體系自由能的計算》（http://www.shanxitv.org/327）所述，溶劑下的體系的自由能G(sol)可以寫為氣相自由能G(gas)和溶解自由能ΔG_solv的加和，因此可以分解成兩個部分：

? 溶解自由能的貢獻Δ(ΔG_solv)：這體現溶劑效應對分子間結合產生的影響，即復合物的溶解自由能減去每個分子的溶解自由能。這部分可正可負，為正代表溶劑效應不利于結合，通常出現在兩個單體間以靜電作用為主的情況，即溶劑效應對單體的總穩定化程度比對復合物的穩定化程度更顯著；為負代表溶劑效應有利于分子間結合，例如《全面揭示各種碳環與富勒烯之間獨特的pi-pi相互作用！》（http://www.shanxitv.org/727）所講的情況，溶劑的疏水效應能顯著促進富勒烯和碳環之間的結合。

溶解自由能可分為極性和非極性兩部分，因此討論溶劑效應對分子間結合產生的影響時可以具體討論這兩部分各自的影響是多大。隱式溶劑模型下溶解自由能的標準計算方式在前述的《談談隱式溶劑模型下溶解自由能和體系自由能的計算》里說了，是溶劑模型下的電子能量減去氣相下的電子能量。Gaussian做考慮SMD溶劑模型的單點計算時可以直接給出非極性部分對溶解自由能的貢獻，把溶解自由能減去這部分就是極性部分對溶解自由能的貢獻。復合物的溶解自由能的極性（非極性）貢獻減去每個單體的這部分，就是溶劑效應的極性（非極性）部分對結合自由能的貢獻。

? 氣相結合自由能ΔG_bind(gas)：通過G_AB(gas) - G_A(gas) - G_B(gas)得到。氣相下G的計算方式見《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》（http://www.shanxitv.org/552）。注意，如果你本身就是研究溶劑環境下的結合，而且溶劑效應對幾何結構的影響不可忽視，那么復合物和單體的opt和freq任務都應當在溶劑模型下做，即計算氣相自由能用的幾何結構是在溶劑模型下優化的，自由能熱校正量也是基于溶劑模型下振動分析信息得到的。

由于G=H-TS，故ΔG_bind = ΔH_bind - TΔS_bind，遂有了下面兩部分。

? ΔH_bind：ΔG_bind(gas)中的焓效應部分，即復合物的焓減去各個單體的焓

? -TΔS_bind：ΔG_bind(gas)中的熵效應部分，即復合物的熵減去各個單體的熵然后乘以-T。由于分子間結合會令熵大幅降低，ΔS_bind為負，因此這部分總是不利于分子間結合，稱為熵罰，而且溫度（T）越大明顯越不利結合。如果分子間無法形成足夠強烈的相互作用令ΔH_bind足夠負，那么它很可能被-TΔS_bind抵消掉導致分子間無法結合。我在北京科音初級量子化學培訓班（http://www.keinsci.com/KEQC）里給了相關計算例子，體現出在氣相標況下，水二聚體在熱力學上是無法自發形成的，但是能在分子間形成兩個顯著的氫鍵的甲酰胺二聚體是能形成的。

溫度對ΔH_bind的影響較小，而對-TΔS_bind影響巨大，因此分子間能否結合極大程度取決于溫度。在《8字形雙環分子對18碳環的獨特吸附行為的量子化學、波函數分析與分子動力學研究》（http://www.shanxitv.org/674）介紹的筆者的文章Phys. Chem. Chem. Phys., 25, 16707 (2023)里通過Shermo程序做了碳環、8字形分子OPP以及二者的復合物的熱力學量相對于溫度的掃描，并進而給出了ΔG_bind、ΔH_bind、-TΔS_bind隨溫度的變化，準確指出了分子間能結合的臨界溫度。前述的《全面揭示各種碳環與富勒烯之間獨特的pi-pi相互作用！》介紹的筆者的Chem. Eur. J., 30, e202402227 (2024)也做了這樣的分析。

如Shermo程序手冊附錄部分所講的熱力學數據計算常識所示，熵是由體系的平動、轉動、振動以及電子（躍遷或簡并）共同貢獻的，Shermo可以分別輸出，而且Shermo還能給出各個振動模式對熵的貢獻。通過這些信息可以深入了解熵的構成，并對各項在復合物和單體之間求差以進一步了解ΔS_bind的內在本質。

由于焓等于電子能量與焓的熱校正量之和，即H = E + H_corr，故ΔH_bind又可以分為以下兩部分：

? 氣相結合能ΔE_bind：即E(AB)-E(A)-E(B)，其中E(AB)、E(A)、E(B)分別是AB、A、B的氣相的（即不帶溶劑模型的）電子能量，并且計算用的幾何結構是對它們各自分別優化后的。ΔE_bind體現出在氣相下，A、B從孤立狀態下結合為AB過程中的電子能量的變化。如果你不清楚什么叫電子能量，看《談談該從Gaussian輸出文件中的什么地方讀電子能量》（http://www.shanxitv.org/488）。注意，平時說結合能的時候可以指ΔG_bind也可以指ΔE_bind，只有前者才能體現熱力學上結合的可能性，后者只能用來討論分子間內在的作用強度、結合的驅動力。沒有前提的情況下說結合能的時候習俗上指的是更容易算的ΔE_bind，本文也是這個習俗。

? ΔH_corr：即復合物的H_corr減去各個單體的。這體現出分子間結合過程中熱力學效應對焓變的貢獻。理想氣體近似下H=U+RT，故雙分子結合的情況有ΔH_corr = ΔU_corr - RT，其中ΔU_corr是結合過程中內能的熱校正量的變化。類似于熵，U_corr也是能夠通過Shermo輸出平動、轉動、振動（包括ZPE在內）、電子躍遷的貢獻從而了解ΔU_corr的內在本質的。

ΔE_bind = ΔE_int + ΔE_deform，即氣相結合能可以分解為氣相相互作用能和變形能，下面專門說一下這兩項。

? 氣相相互作用能ΔE_int：即復合物的電子能量減去各個單體的電子能量，計算用的幾何結構都是優化后的復合物結構（單體的結構直接從優化后的復合物結構中摳出來，不再做優化）。

ΔE_int可以通過筆者提出的高效、普適、很流行的sobEDAw方法分解為分子間的靜電相互作用、交換互斥作用、軌道相互作用、色散作用，細節這里就不說了，在《使用sobEDA和sobEDAw方法做非常準確、快速、方便、普適的能量分解分析》（http://www.shanxitv.org/685）里有極其詳細的說明。

其中，軌道相互作用部分還可以通過《使用Multiwfn通過ETS-NOCV方法深入分析片段間的軌道相互作用》（http://www.shanxitv.org/609）介紹的ETS-NOCV方法進一步將其本質討論得很清楚。

《使用Multiwfn做基于分子力場的能量分解分析》（http://www.shanxitv.org/442）介紹的筆者提出的EDA-FF能量分解方法可以將靜電、交換互斥、色散作用都分解為原子對之間的貢獻和原子的貢獻，從而能對它們的本質和來源有確切的了解。但注意EDA-FF得到的這三項和sobEDAw給出的有所不同，畢竟計算原理相差很大，前者基于力場計算，而后者基于色散校正的DFT理論計算。另外，色散部分還可以按《使用Multiwfn圖形化展現原子對色散能的貢獻以及色散密度》（http://www.shanxitv.org/705）介紹的方式深入地分析。

? 變形能ΔE_deform：在分子形成復合物的整個過程中，單體分子會從孤立狀態的結構（單獨做優化得到的勢能面極小點結構）變成它在復合物中的結構，這會造成分子的能量增加，稱為分子的變形能，可以很容易地手動自行計算。結合能中的變形能ΔE_deform是各個分子的變形能的總和，顯然可以討論各個分子各自的貢獻。分子的柔性越強、分子間相互作用越強，其變形能傾向于越大。分子變形過程對應各個內坐標的變化，因此變形能可以分解為不同內坐標的貢獻來具體討論。還可以用類似《使用Dushin分解重組能和計算Huang-Rhys因子》（http://www.shanxitv.org/330）中提到的分解重組能的方式，把變形能分解為極小點結構下的不同振動模式的貢獻來討論。

以上算是把分子間結合自由能里的所有構成做了詳細的拆分介紹。把這些原理搞清楚了，在討論分子間結合強弱的時候就容易分析得比較深入透徹。如果你由于缺乏基本常識，看了上文之后還是不會計算里面的項，推薦通過北京科音初級量子化學培訓班（http://www.keinsci.com/KEQC）系統性學習一遍，自然就明白了。

順帶一提，有時有人問什么叫吸附能（adsorption energy）、和結合能的區別是什么，這里說一下。吸附能本質就是結合能，只不過前者是后者在特定語境下的稱呼，適用于描述一個相對較小的東西結合到一個相對較大的東西上的情況。比如《8字形雙環分子對18碳環的獨特吸附行為的量子化學、波函數分析與分子動力學研究》（http://www.shanxitv.org/674）一文的情況，18碳環比OPP雙環分子小得多，前者是客體，后者是主體，所以討論二者的結合時可以用吸附能這個詞，叫結合能也完全正確。而一個小分子吸附到一個固體表面上或者多孔材料中的能量變化，這時候就必須用吸附能這個詞描述了。雖然此時說結合能也沒本質錯誤，但由于固體表面/多孔材料比小分子大太多了，所以叫結合能太別扭、不符合習俗用詞。和結合能完全一樣，吸附能也可以用電子能量求差計算也可以用熱力學量求差計算，如果用自由能求差計算，建議明確叫吸附自由能以避免有誤解。

還有個詞叫復合能（complexation energy），和結合能完全是一碼事。注意如《談談BSSE校正與Gaussian對它的處理》（http://www.shanxitv.org/46）所說，Gaussian程序做counterpoise任務的時候會輸出complexation energy，但這是沒有考慮單體的變形能的，因此實質上是當前結構下的相互作用能。讀者應當根據語境恰當區分、變通這些詞。

文/Sobereva@北京科音  2024-Dec-12

1 前言

眾所周知，苯具有極其顯著的芳香性特征。它的重要的等電子體borazine（硼吖嗪）也被稱為無機苯，具有B3N3H6的化學組成，是把苯的六元碳環上的碳替換為交替相連的氮和硼。無機苯的芳香性雖然早就被很多文章研究過，都指出其芳香性明顯不及苯，但不同文章里關于它的芳香性強度說法不一，甚至出入很大，還有的文章使用不嚴謹的方法對其芳香性進行表征得到了誤導性的結論。New J. Chem., 39, 2483 (2015)提出了carborazine（我將它譯為碳硼吖嗪）分子，是苯的另一種等電子體，它的結構介于苯與無機苯之間，如下所示，雖然此文指出了它具有芳香性，但對芳香性研究得并不深入和嚴格。

為了真正全面、透徹闡釋苯及其上述兩種等電子體的芳香性特征、揭示芳香性差異的內在本質，以及了解carborazine的特點，北京科音自然科學研究中心的盧天和江蘇科技大學的劉澤玉等人近期在知名的Chem. Eur. J.（歐洲化學）期刊上發表了專門的研究文章充分研究了上述體系的電子結構和芳香性，非常歡迎大家閱讀和引用：

Yang Wu, Xiufen Yan, Zeyu Liu,* Tian Lu,* Mengdi Zhao, Jingbo Xu,
Jiaojiao Wang, Aromaticity in Isoelectronic Analogues of Benzene, Carborazine and Borazine, from Electronic Structure and Magnetic Property, Chem. Eur. J., 30, e202403369 (2024) DOI: 10.1002/chem.202403369

此文可通過以下鏈接免費在線閱讀：
https://onlinelibrary.wiley.com/share/author/PYSFYG86FNTHBBWXITBI?target=10.1002/chem.202403369

這篇文章還被選為了Chem. Eur. J.的封面文章（DOI: 10.1002/chem.202486603）：

下文將對這篇文章的主要內容做深入淺出的介紹，便于讀者快速順利地了解文章的研究結論，并且同時附加一些額外信息，使得讀者可以更充分了解研究思想和計算細節，從而能夠參考本文的內容更好地做自己的研究。

本文充分利用了強大的Multiwfn程序（http://www.shanxitv.org/multiwfn）做波函數分析討論芳香性，是波函數分析方法和Multiwfn程序的很好且十分典型的應用范例。Multiwfn能做的芳香性分析在《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）里有簡略介紹，在北京科音量子化學波函數分析與Multiwfn程序培訓班（http://www.keinsci.com/WFN）里有多達160頁以上的幻燈片極為完整、全面、系統地講授芳香性的各種概念和全部分析手段，非常推薦參加！

與本文的研究有密切相關的是《18碳環等電子體B6N6C6獨特的芳香性：揭示碳原子橋聯硼-氮對電子離域的關鍵影響》（http://www.shanxitv.org/696）介紹的Inorg. Chem., 62, 19986 (2023)一文，十分建議看完下文后認真閱讀。

2 苯、無機苯和carborazine的結構和穩定性

文中首先用ωB97XD/def2-TZVP對上述三個體系做了幾何優化，這是相當穩的級別，也被用于http://www.shanxitv.org/carbon_ring.html列舉的18碳環及其各種衍生物的研究中，后文的各種分析也都是在這個級別下做的。得到的苯的C-C鍵長1.387埃和電子衍射實驗測定的1.399埃相符極好，得到的無機苯的B-N鍵長1.426埃也完全在X光衍射測定的1.422-1.429埃范圍內。所有三個結構都是純平面的，苯、無機苯和carborazine分別屬于D6h、D3h和C2h點群。

尚未合成的carborazine的穩定性并不得而知，跑從頭算動力學（AIMD）是一種有效的檢驗手段，例如在《18個氮原子組成的環狀分子長什么樣？一篇文章全面揭示18氮環的特征！》（http://www.shanxitv.org/725）介紹的文章里通過AIMD直接證明了18氮環非常不穩定。本文對carborazine和borazine用ORCA程序在ωB97XD/6-311G*級別下做了50 ps AIMD，由于AIMD的昂貴耗時使得模擬時間很有限，文中故意用了1000 K的較高溫度令勢能面的采樣盡量充分。文中還確認了當前溫度下的AIMD過程的平均總能量是高于極小點結構下的ZPE的，因此當前的模擬不會低估實際的振動基態能量。這種模擬在北京科音高級量子化學培訓班（http://www.keinsci.com/KAQC）里有超詳細的講解，在《使用ORCA做從頭算動力學(AIMD)的簡單例子》（http://www.shanxitv.org/576）里講了一些粗略的信息。模擬軌跡多幀疊加圖如下所示，每1 ps繪制一次，按照模擬時間顏色按藍-白-紅變化。可見模擬過程中carborazine和無機苯一樣始終保持結構穩定，并未解離或異構化，而且結構波動范圍相仿佛，這很大程度體現了carborazine像無機苯一樣是穩定的物質，在未來很有可能被成功合成出來。

文中還在ωB97XD/def2-TZVP下計算了原子化能，carborazine的原子化能-1200.9 kcal/mol和無機苯的-1225.6 kcal/mol非常接近，進一步體現出carborazine具有顯著的穩定性。

值得一提的是carborazine的結構并非是其B2C2N2H6的化學組成的所有的可能異構體中能量最低的，這從文中在SI里給出的各種可能的異構體相對于carborazine的能量可以看出，如下所示，可見讓兩個碳挨在一起的能量更低。但由于原子的遷移要破壞原先的C-N和C-B鍵，決速步的能壘必定巨高，因此這樣的自發異構化并不會在現實中發生。

3 苯、無機苯和carborazine的電荷分布

原子電荷的概念在《一篇深入淺出、完整全面介紹原子電荷的綜述文章已發表！》（http://www.shanxitv.org/714）介紹的筆者的文章里有充分介紹，這是定量衡量體系電荷分布最常用的指標。本文通過流行的ADCH方法以及NPA、Hirshfeld和Hirshfeld-I、CHELPG方法都對前述三個體系計算了原子電荷，結果如下所示。可見雖然不同方法計算的原子電荷存在定量差異，但結論都是共通的，即硼和氮由于電負性一個很小一個極大，導致分別帶顯著的正電荷和負電荷，而在carborazine中它們的差異遠小于在無機苯中，這體現出在硼和氮之間插入的碳可以極大地均衡體系中的電荷分布。另外，carborazine中的碳的原子電荷的大小很小，體現出旁邊的氮和硼并沒有對碳的凈電荷帶來明顯影響。

《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）介紹了Multiwfn對pi電子的分析功能。本文使用Multiwfn通過Mulliken方法計算了前述體系的非氫原子的pi布居數以考察它們pi電子的分布情況。結果是苯(C=0.995)、無機苯(B=0.429, N=1.558)、carborazine (B=0.586, C=0.985, N=1.411)。此結果清楚地體現了相對于無機苯，carborazine中在硼和氮之間引入的碳原子可以平衡它們的pi電子數的差異。而carborazine中的碳的pi布居數和它在苯中幾乎沒區別。

為什么要討論環上的pi電子分布？因為它的分布量和分布均衡性與芳香性有關鍵的聯系。較強的pi芳香性必須同時滿足環上的pi電子豐富、pi電子分布較均衡，以及pi電子在環上的離域性較強這三個條件。雖然無機苯的pi電子數和苯一樣，都是6個電子，但環上的原子的pi布居數的最小值，即硼的pi布居數，僅有0.429，因此姑且不論其電子的離域性與苯的差距有多大，光憑這點，其芳香性至多也只有苯的約40%。反之，carborazine六元環上硼的pi布居數達到0.586，因此可以認為其芳香性的上限約有苯的60%。

4 苯、無機苯和carborazine的pi軌道

pi分子軌道的特征與pi電子的芳香性有密切聯系，所以此文也考察了pi軌道的特征。前述三個體系的HOMO、LUMO軌道都是pi軌道，它們的能級和等值面圖如下所示，其它兩個占據的pi軌道也一起顯示了。可見carborazine的HOMO-LUMO gap比苯和無機苯明顯更小，一定程度暗示可能其穩定性相對更弱一些。

Multiwfn可以非常便利地計算軌道成份，見《談談軌道成份的計算方法》（http://www.shanxitv.org/131）。使用Mulliken方法計算的上面圖中展示的占據的pi軌道的成份如下所示。可見無機苯的所有pi軌道全都是由氮原子貢獻所主導的，這也必定導致無機苯的pi電子幾乎都集中且定域在氮原子上，因此在這點上它就不可能有顯著的芳香性。而雖然carborazine的HOMO-3和HOMO-5也是以氮原子的貢獻為主（因為氮的pz原子軌道的能量比碳和硼的都要低，故偏向于對能量較低的分子軌道貢獻較多），但還不至于起到絕對的主導性，而其HOMO的主要貢獻則來自硼和碳。由于carborazine的硼、碳、氮都同時明顯參與了pi占據軌道，它勢必有比無機苯強得多的六中心pi電子離域。

5 苯、無機苯和carborazine的成鍵特征

鍵級的概念在《Multiwfn支持的分析化學鍵的方法一覽》（http://www.shanxitv.org/471）中做了簡要介紹。文中用Multiwfn對上述三個體系計算了Mayer鍵級，它體現兩個原子間等效共享的電子對兒數，結果如下表左側所示，并且還按照《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）所述的方法將sigma和pi電子的貢獻分別給出了。可見雖然B-N、C-N、C-B的pi鍵級都比C-C鍵的更小，但數值也都不低，所以它們的pi作用都是不容忽視的，也因此切勿以為B-N鍵只能算作純sigma鍵。carborazine六元環上的各個鍵的pi鍵級都大于無機苯的B-N鍵，這點已在一定程度上能體現出carborazine的整體pi共軛是強于無機苯的。

文中還對前述體系用Multiwfn做了流行的atom-in-molecules (AIM)拓撲分析，以從鍵臨界點（BCP）處的實空間函數角度考察體系的成鍵特征，結果如上表右側所示。AIM的相關知識和Multiwfn里的操作見AIM學習資料和重要文獻合集（http://bbs.keinsci.com/thread-362-1-1.html）里提到的內容。共價鍵在BCP處的能量密度（H）應當為負，并且eta指數（η）大于1是最典型共價鍵的特征。當前的結果體現出苯的C-C鍵和carborazine的C-N鍵都具有最典型的共價鍵特征，而carborazine和無機苯的N-B的共價性相對最低、離子性相對最強，因此H沒那么負且η顯著小于1，這無疑是硼和氮的巨大電負性差異所帶來的。

要注意切勿把B-N當做離子鍵。有些文章盲目使用BCP位置的電子密度拉普拉斯函數（▽2ρ）的正負作為共價鍵的判據，如Inorg. Chim. Acta., 360, 619 (2007)以這個為理由說B-N是離子鍵，這是很誤導的。之前我在《AIM鍵臨界點處電子密度拉普拉斯值符號判斷相互作用類型失敗原因的圖形分析》（http://www.shanxitv.org/161）中專門說過類似情況。對B-N鍵，如前所示Mayer鍵級顯著大于1，而且BCP處的能量密度明顯為負，就已經充分體現出這屬于共價鍵了。文中更進一步在補充材料里還用Multiwfn繪制了電子定域化函數（ELF）和變形密度圖，分別如下圖左側和右側所示，可見B-N之間有顯著的電子高定域性區域，并且從孤立狀態變成分子過程中在B-N鍵上有顯著的電子密度增加（圖中紅色實線為正值部分），這都進一步證明B-N鍵本質上是共價鍵，只不過極性很強而已。關于變形密度的知識見《使用Multiwfn作電子密度差圖》（http://www.shanxitv.org/113）和《通過價層電子密度分析展現分子電子結構》（http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB201709252），關于ELF的知識見《ELF綜述和重要文獻小合集》（http://bbs.keinsci.com/thread-2100-1-1.html）和《電子定域化函數的含義與函數形式》（http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB20112786），繪制方法參考Multiwfn手冊4.4節，在量子化學波函數分析與Multiwfn程序培訓班（http://www.keinsci.com/WFN）里有更豐富的例子、更細致的講解和各種技巧的傳授。

6 苯、無機苯和carborazine的電子離域特征

如《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）和筆者的Theor. Chem. Acc., 139, 25 (2020)所提到的，pi電子的局域軌道指示函數（LOL），即LOL-pi，對于展現pi電子的離域性極為直觀。對前述三個體系用Multiwfn計算并用VMD繪制的LOL-pi等值面圖如下所示，可見苯的六中心共軛作用最強，等值面均勻、理想地貫穿六元環。雖然carborazine和無機苯的LOL-pi等值面都沒有分布得很均勻，特別是在B-N鍵處受到了阻礙，算是電子共軛的“瓶頸”，但是carborazine的每個B-C-N單元上的等值面較為均勻、連貫，這是無機苯不具有的，因此從LOL-pi展現的電子離域性圖上可以期望carborazine有明顯更強的芳香性。

文中還繪制了EDDB_P-π等值面圖，詳見原文。LOL-pi展現的是不同區域電子的離域程度，而EDDB_P-π展現的是不同區域離域的pi電子量，二者體現的信息存在互補性。EDDB_P-π等值面圖及其布居分析展現出無論是carborazine還是無機苯，硼上的離域的pi電子都比環上的其它原子更少，原因之一是硼上的pi電子數本來就非常少，這在前面的pi布居分析中已經體現了。相比之下，carborazine的硼的離域的pi電子為0.360，明顯多于無機苯的0.184，這進一步體現出carborazine具有更強的pi芳香性。并且0.360這個值是苯的碳原子的離域pi電子數0.885的41%，這說明carborazine至少有苯41%的芳香性。

進一步，文章還使用Multiwfn計算了多中心鍵級（MCBO），以及sigma+core和pi電子各自的貢獻量，相關知識參看《使用AdNDP方法以及ELF/LOL、多中心鍵級研究多中心鍵》（http://www.shanxitv.org/138）和《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）。結果如下所示，可見六中心電子離域性是苯 > carborazine > 無機苯，差異非常鮮明，這和上述其它分析給出的結論非常一致。從MCBO的數值大小來說，carborazine具有苯的大約一半的芳香性，而無機苯則連carborazine芳香性的一半都不到。

7 苯、無機苯和carborazine的磁屏蔽特征

有大量的方法從體系對外磁場的屏蔽角度衡量芳香性，參考筆者寫過的《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）、《使用Multiwfn繪制一維NICS曲線并通過積分衡量芳香性》（http://www.shanxitv.org/681）、《使用Multiwfn計算FiPC-NICS芳香性指數》（http://www.shanxitv.org/724）、《通過Multiwfn繪制等化學屏蔽表面(ICSS)研究芳香性》（http://www.shanxitv.org/216）、《使用Multiwfn巨方便地繪制二維NICS平面圖考察芳香性》（http://www.shanxitv.org/682），在《量子化學波函數分析與Multiwfn程序培訓班》（http://www.keinsci.com/WFN）里有更全面系統的講授。本文從磁屏蔽角度對苯、無機苯和carborazine的芳香性特征進行了定量的考察以及直觀的展現。

NICS及其變體在定量衡量芳香性強弱方面極為常用，文中對苯、無機苯和carborazine計算了NICS最理想的形式之一NICS(1)ZZ，并通過《基于Gaussian的NMR=CSGT任務得到各個軌道對NICS貢獻的方法》（http://www.shanxitv.org/670）介紹的做法分解成了sigma+core和pi電子各自的貢獻，數值越負說明芳香性越強。由數據可見，carborazine的pi芳香性雖然不及苯，但也相當顯著，而無機苯的芳香性極為微弱，接近非芳香性。sigma電子對芳香性的貢獻可忽略不計。

如《使用Multiwfn繪制一維NICS曲線并通過積分衡量芳香性》所述，對NICS_ZZ在垂直于環的方向上積分，即∫NICS_ZZ或稱integrated NICS (INICS)，是比NICS(1)ZZ還更為嚴格的基于磁屏蔽特征衡量芳香性的方法，因為它的結果不依賴于計算磁屏蔽張量的位置的選取。從上面表格的數據可見∫NICS_ZZ的結論與NICS(1)ZZ相一致，更進一步體現了之前結論的可靠性。∫NICS_ZZ計算過程中利用到的穿過環中心且垂直于環方向上的NICS_ZZ曲線可以繪制出來并進行直觀對比，如下所示，橫坐標是相對于環中心（零點）的位置。根據pi電子的NICS_ZZ曲線，可以清楚地看出carborazine的曲線與苯比較接近，這將carborazine顯著的pi芳香性特征體現得非常確切，而無機苯的pi芳香性相比之下顯得微乎其微，其曲線在全范圍都很接近0。三個體系的sigma電子的NICS_ZZ曲線幾乎完全重合，很嚴格地體現出它們的sigma芳香性基本沒有任何區別。

下圖左列是Multiwfn結合VMD繪制的ICSS_ZZ等值面圖，等值面數值為3 ppm，橙色和深藍色分別是對外磁場屏蔽和去屏蔽區域。從這種圖上能非常完整地了解體系的磁屏蔽特征。可見carborazine和苯一樣都是在分子外圍形成了環狀的去屏蔽區域，又一次體現出了carborazine類高度似于苯的顯著的芳香性，而無機苯的去屏蔽區域等值面則明顯顯得不夠連續，體現出其芳香性遠比carborazine弱得多。

上圖還給出了Multiwfn直接繪制的填色等值線圖，中列和右列分別是分子平面以及分子的側截面，顏色越紅和越藍分別對應相應位置的磁屏蔽和去屏蔽效應越強。可以看到carborazine在分子中心是磁屏蔽的，顏色和苯一樣是偏紅的，而borazine在分子中心卻是去屏蔽的，顏色偏藍，這充分反映了carborazine和苯在芳香性上是類似物，而無機苯可近似歸為非芳香性分子一類。

8 苯、無機苯和carborazine的感生電流

在外磁場作用下在環狀區域內離域的電子可以形成感生環電流，這是衡量電子離域性和芳香性的非常直觀的方法，其中AICD展現方式尤為常見，詳見《使用AICD 2.0繪制磁感應電流圖》（http://www.shanxitv.org/294）。本文對苯、無機苯和carborazine分別繪制了總的、sigma電子的和pi電子的AICD圖，如下所示，從AICD(pi)圖中箭頭標注的感生電流上可以清楚地看出carborazine形成了和苯一樣的貫穿整個六元環的顯著的感生電流，而無機苯則完全不具備這樣的特征。AICD圖進一步有力地體現出，carborazine和苯一樣都屬于芳香性物質，而無機苯是近乎非芳香性的，這和ICSS分析的結論一致，彼此相互印證。

此文還使用《使用SYSMOIC程序繪制磁感生電流圖和計算鍵電流強度》（http://www.shanxitv.org/702）中介紹的做法用CTOCD-DZ2方法繪制了分子平面上方1 Bohr處的感生電流圖，如下所示，灰、紅、藍分別對應碳、硼、氮原子。由圖可見CTOCD-DZ2方法展現的感生電流圖和pi電子的AICD圖的現象很類似，而這種平面圖看得更清楚一些。此圖體現carborazine分子平面上方pi區域的感生電流的強度雖不像苯那么均勻，在硼上方相對弱一些，但終究還是形成了明顯的全局感生電流。而無機苯的感生電流則主要在氮原子周圍打轉，體現出pi電子主要都是定域在各個氮的局部原子空間里。

為了便于定量對比，此文還計算了六元環上各種鍵的電流強度（BCS）值，即在鍵中心位置垂直于鍵的平面上的感生電流的積分，結果如上圖標注的數值所示，單位為nA/T。可見carborazine的感生電流強度與苯相差并不懸殊，而無機苯則相當微弱。這從定量角度上進一步驗證了前面基于肉眼觀看感生電流圖所做出的結論。

值得一提的是，在《深入理解分子軌道對磁感生電流的貢獻》（http://www.shanxitv.org/703）中我給了計算感生電流的具體公式，假定只有HOMO-LUMO躍遷對感生電流有顯著貢獻，那么感生電流大小與HOMO-LUMO gap之間會有反比關系。在前面說過，carborazine的HOMO-LUMO gap明顯小于苯和無機苯，因此基于磁屬性和感生電流得到的關于carborazine的芳香性的強度會有一定程度的高估。但即便如此也不影響它具有明顯芳香性的結論，而且前面從多中心鍵級、LOL-pi等分析角度也全都體現了carborazine有遠強于無機苯的芳香性。

9 carborazine中的碳原子對硼和氮原子橋聯的重要意義

前面從各個角度極為充分、全面、嚴格地論證了carborazine具有無機苯不具備的顯著的芳香性，而二者的差異僅在于前者給后者的氮和硼之間引入了碳原子起到橋聯作用。為什么碳原子的橋聯對于提升芳香性有如此關鍵性的影響？此文將Multiwfn、NBO、VMD相結合繪制了前述三個體系的環上的各原子的對應2pz原子軌道的預正交自然原子軌道（PNAO）等值面的疊加圖，過程可參照《使用Multiwfn繪制NBO及相關軌道》（http://www.shanxitv.org/134）舉一反三，如下所示。可見2pz的空間延展程度是硼 > 碳 > 氮，這和電負性增大的順序相一致。氮和硼的軌道重疊頗差，電子明顯難以離域過去。而給二者之間插入的碳原子，增加離域路徑上的原子軌道重疊程度、降低相連原子的電負性不平衡性，明顯可以起到了幫助電子離域過去的“中介”，從而使得carborazine能夠具有無機苯不具備的顯著芳香性。之前筆者在《18碳環等電子體B6N6C6獨特的芳香性：揭示碳原子橋聯硼-氮對電子離域的關鍵影響》（http://www.shanxitv.org/696）介紹的文章中也使用了類似的角度進行分析。

文章里還研究了其它4種苯的等電子體，以考察碳的數目和位置對芳香性的影響。它們的結構如下所示，NICS(1)ZZ也標在了圖上，注意carborazine的值如前所述是-23.1 ppm。通過對比可見，環上的碳原子數多不意味著芳香性就強，例如下圖第1、3個分子雖然有4個碳原子，但單從NICS(1)ZZ來看其芳香性還沒只有兩個碳的carborazine強。下圖第4個分子和carborazine的碳原子數一樣，但碳沒有插入在硼和氮之間，芳香性也遠不如carborazine強。而只有下圖第2個分子，碳橋聯了氮和硼，同時碳原子數比carborazine還更多，這才帶來了稍微更強的芳香性。這個對比體現出對只有硼、氮的環上引入碳，而且能起到對它們的橋聯作用時，才能顯著提升芳香性。

10 總結

本文介紹的Chem. Eur. J. 2024, e202403369一文十分全面透徹探討了互為等電子體的苯、無機苯（borazine，B3N3H6）和尚未實驗合成的carborazine（B2N2C2H6）的電子結構和芳香性的差異。此文的重點是綜合運用各種分析手段，充分論證了含有碳原子橋聯硼、氮原子特征的carborazine具有明顯的芳香性（盡管強度一定程度弱于苯），同時還嚴格地證明了無機苯的芳香性極弱，甚至很大程度上體現出非芳香性分子的特征。文章對此現象的本質從電子結構層面進行了詳細的分析和解釋。此文的研究工作對于化學家們更好地了解由碳、氮、硼構成的物質的電子離域性和芳香性有顯著的意義。

值得一提的是，此文的研究也是使用Multiwfn等程序做波函數分析深入剖析化學物質特征的非常典型的文章，可以作為波函數分析的很好的范例作為參考。此研究涉及的研究方法在前面提及的各種博文里有一些粗略、初步的介紹，通過量子化學波函數分析與Multiwfn程序培訓班（http://www.keinsci.com/WFN）可以一次性全面深入透徹掌握文中用到的一切波函數分析方法的背景知識和使用具體程序計算的方法，從而能順利地復現此文的研究并舉一反三探索更多的體系。

文/Sobereva@北京科音

First release: 2024-Dec-8   Last update: 2025-Feb-12

筆者開發的免費的分子構象和團簇構型搜索程序Molclus（http://www.keinsci.com/research/molclus.html）如今已經非常流行，在網上經常有人提問此程序的使用問題。這里我把一些比較常見的問題進行統一回答。日后還會對本文內容有進一步追加。

Q1：Molclus是干什么？
A：如果你沒有任何量子化學背景知識，先把《辨析計算化學中的任務類型和理論方法》（http://www.shanxitv.org/680）看了。多原子分子、原子/分子團簇體系的勢能面上一般有不止一個極小點，直接基于自己猜的初始結構去用Gaussian等程序做幾何優化，極可能會優化到不是能量最低的極小點，這會導致之后做的分析討論沒實際意義、結果有誤導性。Molclus可以幫助你獲得能量最低的極小點，即全局最小點的結構，這是體系最有意義的結構。

此外，柔性分子、團簇的勢能面上通常有一批能量都很低、能量相差不大的極小點結構，在實際研究的溫度下都有一定出現比例（參見《根據Boltzmann分布計算分子不同構象所占比例》http://www.shanxitv.org/165），計算構象權重的光譜（參見《使用Multiwfn繪制構象權重平均的光譜）http://www.shanxitv.org/383和《使用Multiwfn繪制NMR譜》http://www.shanxitv.org/565）以及計算構象權重的熱力學量（參見《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》http://www.shanxitv.org/552）需要得到這些低能量極小點的結構和出現百分比，Molclus也可以很好地實現這個目的。

Q2：怎么引用Molclus程序？
A：計算化學程序的引用在《寫計算化學文章時引用理論方法、基組、程序時應注意的問題》（http://www.shanxitv.org/370）里明確強調了。使用Molclus發文章（包括給別人代算）必須按照主頁http://www.keinsci.com/research/molclus.html上明確寫明的方式引用，引用必須在正文里體現。Molclus完全免費，按要求正確引用Molclus是對此程序開發和維護的最好的支持。

Q3：使用Molclus程序遇到問題如何求助？
A：在計算化學公社論壇（http://bbs.keinsci.com）的“量子化學”版塊發帖，發帖時選擇Molclus分類。也可以在筆者創建的思想家公社QQ群里提問，群介紹和群號見http://www.shanxitv.org/QQrule.html。我看到后都會盡量回答。提問時必須把使用細節交代得盡可能清楚，避免含糊其辭或引起歧義，注意看《在網上求助計算化學問題的時候必須把問題描述得詳細、具體、準確、清楚》（http://www.shanxitv.org/620）。

Q4：Molclus程序怎么上手？
A：Molclus主頁http://www.keinsci.com/research/molclus.html上明確給了具體教程，寫得極為詳細，認真一個字一個字閱讀、照著操作即可。

需要注意的是，Molclus雖然非常易于使用、靈活方便，但不代表這是完全黑箱程序，使用者不能一點計算化學常識都沒有。Molclus需要調用第三方程序來做單點、幾何優化、振動分析等任務。Molclus調用哪個程序，用戶就應當有相應程序的最基本常識。例如使用Molclus調用Gaussian做優化和單點能計算是最常涉及的情況，因此用戶就得懂Gaussian最基本的使用，要不然連Gaussian怎么安裝和運行、Gaussian模板文件里怎么寫關鍵詞、最常見的Gaussian報錯都不知道怎么解決可不行。而且僅當你具備基本的Gaussian計算常識的情況下，才能知道怎么針對當前的體系根據精度要求和計算資源選擇合適的級別用Molclus做構型或構象搜索。不具備相關知識的話非常建議通過北京科音初級量子化學培訓班（http://www.keinsci.com/KEQC）快速且系統地學習。

Q5：Molclus做構型/構象搜索的基本流程是什么？
Molclus的最簡單的使用過程是：
(1)用Molclus的gentor或genmer或第三方程序產生一批初始構型/構象記錄在traj.xyz文件里，這是多幀的xyz格式的文件，見《談談記錄化學體系結構的xyz文件》（http://www.shanxitv.org/477）
(2)用Molclus對當前目錄下的traj.xyz里的構型/構象依次做優化，得到的結構和能量都會儲存在isomers.xyz里
(3)用Molclus里的isostat工具對isomers.xyz中的結構去重、按能量排序，得到cluster.xyz，里面第一個就是全局能量最小點結構，后面的是能量次低的結構

更準確而且也是實際中特別常用的流程如下，這里以普通有機體系為例
(1)如上所述先產生一批初始構型/構象
(2)用Molclus結合xtb程序在很便宜但也較糙的GFN-xTB級別下對這些構型/構象依次做很快速的優化（對柔性大體系這一步可能要對幾百甚至上千個結構做計算)
(3)用isostat對優化后的結構去重、按能量排序得到cluster.xyz
(4)將cluster.xyz改名為traj.xyz，用Molclus調用Gaussian在精度夠用而且也不貴的B3LYP-D3(BJ)/6-31G*級別下對之前能量最低的幾kcal/mol的結構（通常在幾十個以內）進行更準確的優化以及做振動分析（后者是為了檢驗虛頻情況以及獲得自由能熱校正量），并且之后自動在M06-2X/def2-TZVP級別下做較高精度的單點能計算。此時isomers.xyz里存的是各個結構及其較準確的自由能
(5)用isostat對isomers.xyz里的結構去重、按自由能排序得到cluster.xyz，并且可以輸入溫度讓isostat直接算出來各個結構的Boltzmann出現比率

記住如果是做溶液情況下的構型/構象搜索，單點計算必須帶溶劑模型（建議用SMD），優化和振動分析最好也帶上以圖穩妥（建議用IEFPCM）。如果嫌振動分析貴也可以不做，最后isostat會將按照電子能量而非自由能排序，給出的Boltzmann分布的準確性可能會差不少。

沒有哪種構型/構象搜索流程對所有情況都是最理想的，精度和效率不可能同時特別理想，上面說的只是對一般有機體系的精度和效率比較均衡的搜索流程。使用者可根據實際情況隨機應變。具體操作細節在Molclus官網上的教程里都有體現。

上述提到的Gaussian是收費的，因版權問題實在不方便用Gaussian的話也可以讓Molclus調用免費的ORCA量子化學程序來代替，再加上xtb和Molclus都是免費的，因此利用Molclus做構象搜索可以全程免費，只不過ORCA在優化+振動分析方面的整體穩健性方面不如Gaussian。用ORCA的話，上述優化+振動分析部分建議用B97-3c，單點計算部分建議用wB97M-V/def2-TZVP結合gCP形式的BSSE校正。ORCA和xtb的使用在《北京科音高級量子化學培訓班》（http://www.keinsci.com/KAQC）里有極其全面的介紹。

Q6：哪種方式產生初始構型/構象最適合我？
A：不同方式產生初始構型/構象各有長處，在Molclus官網上寫明了，如下所示。有不清楚的可以在公社論壇或QQ群里向我咨詢
? Molclus自帶的genmer程序結合Molclus：用于分子團簇或原子團簇構型搜索，使用簡便。對于分子團簇搜索這通常是首選。唯一局限性是對分子的構象考慮不夠充分，因為genmer將分子當成剛性來堆積，因此不適合單分子柔性很強的情況，這種情況更適合后面說的做動力學程序結合Molclus。相關介紹和示例看《genmer：生成團簇初始構型結合molclus做團簇結構搜索的超便捷工具》（http://bbs.keinsci.com/thread-2369-1-1.html）
? Molclus自帶的gentor程序結合Molclus：用于分子構象搜索。使用簡便，但不支持環狀區域構象搜索（環狀區域是指諸如環己烷這樣有多種構象的環狀區域）。對于可旋轉的鍵不很多的分子的構象搜索這是首選，可控性最強。相關介紹和示例看《gentor：掃描方式做分子構象搜索的便捷工具》（http://bbs.keinsci.com/thread-2388-1-1.html）
? 第三方的Confab或Frog2構象生成程序結合Molclus：用于分子構象搜索。使用最為傻瓜化，但不支持環狀區域構象搜索、不支持普通有機分子以外的情況。相關介紹和示例看《將Confab或Frog2與Molclus聯用對有機體系做構象搜索》http://bbs.keinsci.com/thread-20063-1-1.html）
? xtb程序跑動力學軌跡結合Molclus：普適性極強，用于構象搜索、原子/分子團簇構型搜索皆可，使用容易。但由于模擬時間長度有限，有動力學采樣不充分導致遺漏構型/構象的風險（風險隨動力學模擬時間增長而減小）。相關介紹和示例看《使用Molclus結合xtb做的動力學模擬對瑞德西韋(Remdesivir)做構象搜索》（http://bbs.keinsci.com/thread-16255-1-1.html）
? GROMACS等經典力場分子動力學程序結合Molclus：用于分子團簇搜索和分子構象搜索。使用稍繁瑣，因為得創建拓撲文件，且被動力學模擬的體系必須有恰當的力場，好處是動力學模擬耗時極低，因此動力學采樣不充分導致遺漏構型/構象的風險很小。相關介紹和示例看《使用molclus程序做團簇構型搜索和分子構象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）。不了解GROMACS使用的話推薦通過《北京科音分子動力學與GROMACS培訓班》（http://www.keinsci.com/KGMX）學習相關使用知識。

Q7：為什么運行Molclus的時候閃退/崩潰？
只有兩個可能：
(1)被Molclus所讀取的traj.xyz的格式有問題導致程序崩潰。要按照《談談記錄化學體系結構的xyz文件》（http://www.shanxitv.org/477）說的用文本編輯器檢查其內容看是否符合標準格式。也可以試圖往免費的VMD可視化程序里載入，如果載入不了也說明格式有問題
(2)被Molclus所讀取的settings.ini文件有問題。典型情況是用戶修改里面設置的時候把格式改亂了，導致無法被讀取。也有的用戶干了多余的事，自作主張地把他覺得不需要的內容給刪掉了，導致Molclus沒法讀取要讀取的內容。還有的用戶直接把老版本Molclus的settings.ini放到新版本Molclus目錄下用，由于新、老版本此文件里面包含的選項不同，導致新版本Molclus無法正常讀取之（所以不要干這種事，要改什么設置一定要在當前版本Molclus自帶的settings.ini基礎上按照原本的格式改）
死活搞不明白的話，到計算化學公社論壇發帖，把整個Molclus目錄壓縮后作為附件上傳。

Q8：為什么運行Molclus自帶的genmer或gentor后程序閃退/崩潰？
A：要么是涉及的xyz文件格式有問題，自行檢查；要么是genmer的genmer.ini、gentor的gentor.ini里的格式沒寫對或者寫的內容有bug。死活搞不明白的話，到計算化學公社論壇發帖，把整個genmer或gentor目錄壓縮后作為附件上傳。

Q9：為什么產生的isomers.xyz內容為空？
A：這說明連一個結構都沒跑成功。一定要仔細看Molclus運行窗口中的提示，Molclus在執行什么、運行狀態如何，都提示得清清楚楚。導致此問題常見的原因如下，是哪種情況根據Molclus輸出的提示就能很容易地判斷（實在看不懂的話在提問時把完整截圖貼出來）
(1)Molclus調用Gaussian、xtb等第三方程序失敗。通常是第三方程序在機子里沒裝好，或者settings.ini里調用它們的命令沒寫對
(2)連一個任務都還沒跑完。如果你確信你的計算設置沒有離譜的地方，就繼續等著
(3)所有任務都運行失敗了。如果把settings.ini里的ibkout設為1，每個任務的被調用的第三方程序的輸出文件都會被保留，應當根據里面的信息結合相應程序的使用常識判斷當前是什么原因所致。常見的問題包括：模板文件里關鍵詞寫錯了、凈電荷或自旋多重度設錯了、SCF/幾何優化不收斂。

Q10：怎么Molclus老也算不完？
A：如果根據Molclus的提示能確信任務正在算著，這種情況無非就是以下原因
(1)CPU忒爛導致耗時太高。一些初學者對耗時沒基本常識，在個位數核心的CPU上就想讓Molclus調用Gaussian或ORCA用DFT優化上百原子體系、對其做振動分析，耗時不高才怪。牽扯到大批量上百原子體系的DFT計算都應該有個像樣的性能較好的服務器，要不然做構象搜索可能得跑N天甚至一周以上。也可能使用者的CPU并不差，但沒恰當并行，或者非要用Windows版的Gaussian，導致CPU性能無法充分發揮
(2)用的計算級別過于昂貴。諸如讓Molclus調用Gaussian對挺大的體系用def2-TZVP這么貴的基組做優化、對其用雙雜化泛函做單點計算甚至優化，耗時不高才怪
(3)遇到了難收斂。應根據被調用的程序的使用常識，根據其臨時輸出文件內容判斷當前是什么情況。Molclus調用第三方程序用的命令在屏幕上提示得清清楚楚，被調用的程序產生的臨時輸出文件是哪個也體現得明明白白，比如Gaussian就是gau.out。對于調用Gaussian做幾何優化，若一個任務老跑不完，顯然應該想到用GaussView打開gau.out，看看是遇到了震蕩難收斂還是什么其它情況，以及弄清楚當前跑多少步了、是不是每一步的耗時太高，等等。如果真遇上幾何優化難收斂，Gaussian用戶可以參考《量子化學計算中幫助幾何優化收斂的常用方法》（http://www.shanxitv.org/164）在模板文件里加上適當的幫助收斂的關鍵詞。對于Molclus調用ORCA也是類似的處理思路

Q11：Molclus對各個結構優化時有個別沒成功，結果可以用么？
A：我只能說沒優化成功的任務占總結構數目比例越小，遺漏重要結構的概率越低。如果你用ibkout=3把失敗的任務的輸出文件都備份了出來，可以看看其優化最后一幀的結構和isostat排序后能量最低一批結構的差異，如果看上去優化失敗的結構繼續跑下去也不太可能收斂到能量較低的結構去，那這樣的結構就可以舍棄。而如果你懷疑它繼續優化有可能能得到一個不可忽視的能量很低的結構，那就把他回爐結合適當的關鍵詞重新優化。值得注意的是settings.ini里可以通過ngeom指定對traj.xyz里哪些幀進行處理，iappend設1可以把新得到的結果在之前的isomers.xyz上續寫，因此可以方便地將回爐的結果累加到之前的上面。

Q12：Molclus能在Windows下用么？
A：能是能，但在Linux下用更理想。xtb程序雖然有Windows版，但有些bug，Windows下使用時容易莫名其妙中斷。因此筆者在Windows主機下用Molclus調用xtb時都是在VMware虛擬機里裝的Linux客戶機里用。而且Gaussian的Windows版的計算效率不及Linux版，尤其是CPU核數較多時相差很懸殊，這點在《Gaussian的安裝方法及運行時的相關問題》（http://www.shanxitv.org/439）里專門說了。

Q13：為什么Windows下用Molclus調用Gaussian計算時蹦出Gaussian圖形窗口后卡住不動？
A：那是因為你把settings.ini里gaussian_path的路徑設成了比如g16w.exe（對于Gaussian 16來說），帶w后綴的只是Gaussian 16圖形界面這層殼，g16.exe才是真正的Gaussian 16程序本體，gaussian_path必須指定成它的路徑。

Q14：為什么Molclus調用xtb不能運行，屏幕上提示諸如：'xtb'不是內部或外部命令，也不是可運行的程序
A：顯然是xtb程序沒在機子里裝好，倘若裝好了就直接可以用xtb命令調用xtb程序，自然不會有這種提示。仔細看《將Gaussian與Grimme的xtb程序聯用搜索過渡態、產生IRC、做振動分析》（http://www.shanxitv.org/421）了解xtb怎么安裝。

Q15：為什么Windows下用Molclus調用Gaussian計算時提示No executable for file l1.exe？
A：沒有把Gaussian的文件夾加入GAUSS_EXEDIR環境變量所致。這點在《使用molclus程序做團簇構型搜索和分子構象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）里專門說了，在文中搜GAUSS_EXEDIR。或者你加了但路徑沒寫對，也可能環境變量尚未生效（嘗試重啟機子再試）。

Q16：genmer產生構型時怎么考慮單分子的構象？是否需要先對單分子做構象搜索？
A：genmer把每個分子當做剛性看待，所以沒法考慮。如果單分子的柔性較強，因此形成復合物的時候構象可能發生極大的變化，那么建議使用分子動力學方式產生團簇初始結構的traj.xyz文件。對于柔性分子先對單分子做構象搜索意義不大，畢竟分子間相互作用可能導致復合物中的分子的構象相對于在孤立狀態下發生很大變化。

Q17：genmer產生多少個初始構型合適？
A：產生的初始構型越多，顯然經過批量優化后遺漏掉能量最低構型的可能性就越低。對于分子較小、分子數較少的團簇的構型搜索，比如水的四聚體，初始構型產生大幾十個就足夠了。分子越大、分子數越多的團簇，由于勢能面維度更高、極小點更多，原理上應該產生越多的初始構型。鑒于調用xtb程序做GFN-xTB級別的初篩耗時非常低（調用GFN-FF分子力場的話耗時更低），優化大幾百、上千個都完全不是難事，因此在計算量可以接受的情況下對較大的體系產生的初始構型原則上多多益善。也可以根據之前得到的結果憑直覺判斷是不是還有遺漏能量更低構型的可能，如果懷疑有這種可能性，可以讓genmer再產生一批構型，讓Molclus再次批量優化并續寫到之前的isomers.xyz上。

Q18：gentor該旋轉哪些鍵？
A：旋轉的鍵越多，產生的初始構象數顯著越多，導致構象搜索耗時明顯越高，因此不能太多。只有柔性的、現實環境下易于旋轉的鍵才應該在gentor里要求旋轉來產生各種初始構象，這憑借基本常識就能判斷。現實中不易旋轉的鍵有：(1)有單套pi共軛作用的鍵，比如肽鍵就有pi共軛使其C-N鍵的旋轉勢壘巨高、在常溫下不可能自發旋轉 (2)旋轉后會嚴重破壞之前存在的顯著的吸引性的分子內弱相互作用，比如導致破壞好多氫鍵 (3)旋轉后會出現嚴重的位阻作用而且無法靠旋轉其它部分回避。

也不是所有易于旋轉的鍵都應當納入考慮，這憑常識就能判斷，比如甲基雖然在常溫下極易旋轉，但旋轉前后的結構是相同的，自然旋轉它毫無意義。

如果你實在缺乏常識，索性用confab產生初始構象，但不像gentor那樣有可控性。

Q19：isostat的能量和結構去重閾值怎么設置合適？
A：若無特殊情況就用默認設置即可。閾值多大合適和計算級別的精度無關，而取決于幾何優化的收斂精度。默認的去重閾值適合幾何優化收斂限偏嚴的情況（如Gaussian默認的收斂限）。當你用了較寬松的幾何優化收斂限，去重閾值也需要相應地放寬。

如果發現得到的cluster.xyz中有肉眼看上去特別接近的結構，說明去重閾值偏小，導致本應該只保留其中一個的情況下多個結構都被保留了，此時可以稍微加大閾值（如變成默認的兩倍）重新嘗試。但閾值也不能太大，否則能量和結構比較接近的對應實際不同極小點的構象會只保留一個。

Q20：isostat判斷重復結構的原理是什么？
A：在《使用molclus程序做團簇構型搜索和分子構象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）里專門說了，仔細看相關部分。

Q21：怎么對過渡態做構象搜索？
A：不小的體系的化學反應只發生在體系的一部分區域，可以叫反應區域，整個反應路徑中只有這部分區域的原子涉及的內坐標有特別顯著變化、牽扯化學鍵的改變，基于這點可以做過渡態的構象搜索。首先照常搜索過渡態，然后用gentor對非反應區域做各種二面角的旋轉產生traj.xyz，在Molclus做批量優化時用settings.ini里的freeze設置將反應區域的原子凍住，之后照常用Molclus做剩下的構象搜索步驟即可。當然，反應區域和非反應區域會存在一定耦合，非反應區域處于不同構象時反應區域的結構也會多多少少有所不同。因此為了嚴謹起見，建議以上述構象搜索得到的能量最低一批結構作為初猜再次優化過渡態，并取其中能量最低的過渡態。

Q22：為什么用OpenBabel通過Confab方法產生一批構象時失敗？
A：和OpenBabel版本可能有關，讀者請下載《將Confab或Frog2與Molclus聯用對有機體系做構象搜索》（http://bbs.keinsci.com/thread-20063-1-1.html）里直接提供的OpenBabel版本，并且先確保對帖子里的體系能正常產生構象然后再做自己的體系。并且要參考帖子確保輸入的OpenBabel命令合理，包括選項的大小寫。OpenBabel不要裝到默認的C盤的目錄下，免得由于操作系統可能的權限限制導致OpenBabel無法成功創建新文件。也要確保提供給OpenBabel的mol2文件格式合理，參考《談談記錄化學體系結構的mol2文件》（http://www.shanxitv.org/655）。

Q23：Molclus怎么在超算上用？
A：筆者不用超算，也并未特意考慮Molclus在超算上用的情況。但Molclus是完全可以在超算上跑的，參考以下帖子
《consearch：一鍵提交slurm的molclus構象搜索腳本》（http://bbs.keinsci.com/thread-43932-1-1.html）
《分享molclus配套PBS作業提交腳本》（http://bbs.keinsci.com/thread-15598-1-1.html）
《分解合并Molclus任務的小腳本splitMC》（http://bbs.keinsci.com/thread-14361-1-1.html）

Q24：為什么Molclus做構型/構象搜索得到的能量最低結構比我一開始自己手動建模和計算得到的還更高？
A：有兩個可能原因
(1)數據沒有可比性。比如一個取的是電子能量，一個取的是自由能。或者你自己讀數據沒讀對地方。或者計算級別（理論方法和基組）、設置（溶劑模型、DFT積分格點、算熱力學量用的溫度等）不嚴格相同。或者用的計算程序版本不同。諸如此類問題必須自己反復認真檢查，別人在遠程難以跟你說清楚。
(2)molclus做構型/構象搜索采樣不夠充分，這是產生traj.xyz的方式不理想所致（如genmer產生的構型數太少、gentor忽略了一些重要的二面角的旋轉、跑動力學用的溫度太低或時間不夠長，等等），或者搜索流程不當所致（比如預篩選的時候用的計算級別不當而導致把關鍵結構給篩掉了），這會使得能量最低的結構沒被最終搜索出來，而你自己擺的結構又恰好在優化后收斂到了能量更低的結構。
molclus做構型/構象搜索的原理易于理解、全流程都是十分透明的，結合量子化學計算的常識、對勢能面的理解、對Molclus的使用邏輯的認識，用戶很容易就能搞清楚導致這種問題的原因。