文/Sobereva@北京科音   2025-Feb-21

1 前言

《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）里介绍的笔者提出的图形化展现相互作用的方法IGMH已被广为使用，之前我写的《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》（//www.umsyar.com/727）和《谈谈pi-pi相互作用》（//www.umsyar.com/737）里涉及到了IGMH的等值面面积，至少对pi-pi作用来说它和相互作用强度有密切的正相关性。有不少读者都问我怎么得到面积，本文就专门说一下。计算面积的同时还会顺带得到等值面内包围的体积。此文说的IGMH等值面具体是指IGMH方法里定义的δg_inter函数的等值面。如果读者不熟悉IGMH，应先把//www.umsyar.com/621看了。读者应使用Multiwfn官网//www.umsyar.com/multiwfn上的最新版本以免和本文的情况不符。不了解Multiwfn者参看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。

下面介绍两种做法，第一种方法是使用Multiwfn的定量分子表面分析（主功能12）计算IGMH等值面面积和体积，这种情况只适合存在一个等值面，且这个等值面就是你要研究的等值面的情况。另一种方法更为普适，需要借助免费的ChimeraX程序载入Multiwfn产生的cub文件显示IGMH等值面，可以测量图中任意一个等值面的面积和体积，对于《使用IRI方法图形化考察化学体系中的化学键和弱相互作用》（//www.umsyar.com/598）介绍的IRI函数、《谈谈范德华势以及在Multiwfn中的计算、分析和绘制》（//www.umsyar.com/551）介绍的范德华势等各种函数也都可以这么测量。

2 使用Multiwfn的定量分子表面分析功能计算IGMH等值面的面积和体积

这一节以18碳环与C60富勒烯的复合物为例演示怎么直接用Multiwfn得到IGMH等值面的面积和体积。前述的《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》介绍的笔者的Chem. Eur. J., 30, e202402227 (2024)研究中得到的此体系的波函数文件C60-C18.wfn在//www.umsyar.com/attach/738/file.rar里提供了。此体系的IGMH图如下所示（δg_inter函数等值面为0.002 a.u.）

首先将Multiwfn目录下的settings.ini里的iuserfunc设为91，这代表把用户自定义函数（user-defined function）设为IGMH方法的δg_inter函数。在Multiwfn手册2.7节里有可用的用户自定义函数的完整列表。之后启动Multiwfn，载入C60-C18.wfn，之后依次输入
1000   //隐藏功能
16  //定义片段
2  //定义两个片段
1-18  //18碳环里的原子序号范围
19-78   //富勒烯里的原子序号范围
12  //定量分子表面分析功能
1  //设置用于定义表面的函数
2  //某个实空间函数
100  //用户自定义函数
0.002  //定义表面用的等值面数值
6  //开始分析，不考虑被映射的函数
接下来程序开始计算δg_inter的格点数据，过一会儿，在屏幕上看到以下信息
 Volume:    69.13431 Bohr^3  (  10.24465 Angstrom^3)
 Estimated density according to mass and volume (M/V):  151.8505 g/cm^3
 Overall surface area:         323.36643 Bohr^2  (  90.55182 Angstrom^2)

在后处理菜单选-3可以观看当前考察的δg_inter函数的0.002 a.u.等值面，如下所示，确实就是前面给出的笔者的论文Chem. Eur. J., 30, e202402227 (2024)里的那个等值面。其体积是上面显示的10.24 Å^3，面积是90.55 Å^2，和文中报道的一致。

顺带一提，在后处理菜单中还可以选-2将当前算出来的δg_inter的格点数据导出为surf.cub，之后可以被第三方程序可视化和分析。

还值得一提的是进入主功能12的时候可以看到选项3 Spacing of grid points for generating molecular surface用来设置格点间距，默认是0.25 Bohr，数值设得越小计算耗时越高而统计精度越高。

Multiwfn对各种实空间函数（包括从外部的.cub等格点数据文件读入的）都可以利用定量分子表面分析功能计算其等值面的面积和体积，另一个使用例子见《使用Multiwfn计算轨道的体积》（//www.umsyar.com/734）。

3 使用Multiwfn结合ChimeraX获得IGMH等值面的面积和体积

这一节以2-pyridoxine和2-aminopyridine的二聚体为例演示利用ChimeraX程序得到特定的IGMH等值面的面积和体积。examples\2-pyridoxine_2-aminopyridine.wfn是Multiwfn程序包里自带的这个体系的波函数文件，在0.01 a.u.的δg_inter等值面数值下分子间的IGMH等值面图如下所示，可见有两个等值面，此例分别获得它们的面积和体积

首先照常对这个二聚体做IGMH分析。启动Multiwfn，然后依次输入
examples\2-pyridoxine_2-aminopyridine.wfn
20  //弱相互作用可视化分析
11  //IGMH分析
2  //两个片段
1-12  //第1个片段
c  //其它原子作为第2个片段
4  //设置格点间距（格点分布覆盖整个体系）
0.2  //格点间距为0.2 Bohr
3  //导出格点数据
当前目录下有了dg_inter.cub，是δg_inter的cub文件。

之后退回到Multiwfn主菜单，输入xyz后按回车，再输入2-pyridoxine_2-aminopyridine.xyz，当前目录下就得到了记录当前结构的2-pyridoxine_2-aminopyridine.xyz文件。

去https://www.rbvi.ucsf.edu/chimerax/download.html下载ChimeraX并安装。本文用的是ChimeraX 1.9。

启动ChimeraX，将2-pyridoxine_2-aminopyridine.xyz和dg_inter.cub依次拖入程序界面载入，然后左右随便拖动一下下图红箭头所示的竖杠，激活这个等值面的设置，然后再在下图蓝箭头所示的文本框里输入要用的等值面数值0.01，之后看到的等值面就是下图这样

之后选择窗口上方的Tools - Volume Data - Measure and Color Blobs，之后按住Alt键并点击你要考察面积和体积的那个等值面，那个等值面就被自动着色了，并且在界面右下角显示了其体积和面积，如下所示，分别为0.31 Å^3和3.28 Å^2。

值得一提的是，由于ChimeraX和Multiwfn的定量分子表面分析功能产生等值面的算法不同，因此得到的等值面的面积和体积会有轻微差异。例如用ChimeraX对上一节的18碳环和富勒烯之间的δg_inter=0.002 a.u.的等值面进行测量，得到的面积是89.90 Å^2，体积是10.65 Å^3，面积和Multiwfn给出的90.55 Å^2相差0.7%。

笔者之前还录过一个视频《使用Multiwfn和ChimeraX绘制自定义着色的电子定域化函数(ELF)等值面图》（https://youtu.be/vC48iEB8PwI、https://www.bilibili.com/video/av85684420）演示ChimeraX里的和等值面显示、测量相关的操作，感兴趣的读者可以看看。

文/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》（//www.umsyar.com/452）。

我有不少研究文章都和pi-pi作用有密切联系，是pi-pi作用研究的典型范例，非常推荐读者阅读并欢迎引用：
• Theoretical Insight into Complexation Between Cyclocarbons and C60 Fullerene, Chem. Eur. J., 30, e202402227 (2024)。在《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》（//www.umsyar.com/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堆积特征》（//www.umsyar.com/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碳环的独特吸附行为的 、波函数分析与分子动力学研究》（//www.umsyar.com/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）与石墨烯的相互作用：基于簇模型的研究一例》（//www.umsyar.com/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电子结构特征》（//www.umsyar.com/432）的做法绘制出的几何结构结合pi电子等值面图如下所示，两个分子的pi电子的分布非常清楚直观，可以看到彼此相对。注意并非必须两个分子的所有pi电子都彼此对着才算pi-pi作用，像下图这样哪怕只要一部分对着也同样算pi-pi作用。而如果完全没有对着、完全错开了，那就不算pi-pi作用了。更多的展现pi电子分布的图见笔者的Theor. Chem. Acc., 139, 25 (2020)一文。

2 pi-pi作用的物理本质

我在《谈谈“计算时是否需要加DFT-D3色散校正？” 》（//www.umsyar.com/413）提到过，原子间相互作用从物理本质上可以基本分为静电、色散、交换、Pauli互斥、轨道相互作用；我在《谈谈“计算时是否需要加DFT-D3色散校正？”》（//www.umsyar.com/413）里也说过物理本质和相互作用表象之间的关系。pi-pi作用类似于氢键、卤键等，是属于表象层面的，用来描述一种具有特定特征的相互作用，其内在本质就是前面说的来自于相距较近的pi电子之间的色散作用，这一点在Grimme的一篇经典的讨论pi-pi作用本质的文章Angew. Chem. Int. Ed., 47, 3430 (2008)中通过严格的理论计算和分析对比已经论证得十分严格充分，而且也早已是 领域内行人的共识。

有的文章试图从轨道相互作用解释pi-pi堆积出现的本质，这是极其有误导性的，却还广泛以讹传讹。2007年的时候wiki上的词条说pi-pi作用来自于pi共轭体系间p轨道的重叠，前述Angew文章的脚注中专门对此进行了批判，指出没有任何 计算能支持这种说法。这也不难理解，本来两个分子的pi轨道之间重叠程度极低（可以用《分子间轨道重叠的图形显示和计算》//www.umsyar.com/163介绍的方法考察），远低于形成共价键时原子轨道间的重叠程度，因此有效重叠不足。而且两个分子的pi占据轨道之间混合会产生全占据的成键和反键轨道，二者效果充分抵消了，对成键没贡献；而若一个分子的pi占据轨道和另一个分子的pi非占据轨道混合，意味着要出现分子间pi->pi*电子转移，这通常又没有明显能使体系能量变得更低的内在驱动力（不像常见的较近距离的n->pi*超共轭往往可以明显降低体系能量）。此外，按照《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）说的，用Multiwfn对pi-pi堆积的两个分子间去计算衡量两个部分之间有效共享电子对数的Mayer键级、模糊键级或离域化指数（DI），会发现数值很接近0，比如《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》（//www.umsyar.com/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堆积二聚体，还有很多其它体系都是色散吸引作用驱动分子间结合，而静电作用进一步影响几何结构。《静电效应主导了氢气、氮气二聚体的构型》（//www.umsyar.com/209）介绍的笔者的研究就是很好的例子，H2、N2的各种二聚体的相互作用能通过能量分解分析可以发现都是色散作用为主，但不同构型下的静电作用的差异则直接影响不同构型的稳定性乃至存在与否，并且笔者发现通过静电势互补原理可以给出很直观的解释。这种静电势互补的分析方法也在《全面探究18碳环独特的分子间相互作用与pi-pi堆积特征》（//www.umsyar.com/572）介绍的笔者的论文中也用来解释为什么18碳环的pi-pi堆积二聚体是错位的。平面环状结构的18碳环具有独特的平面内和平面外两套pi电子，此文中确认了18碳环能够形成分子间彼此平行的pi-pi堆积二聚体，下图是按照《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图》（//www.umsyar.com/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方法做非常准确、快速、方便、普适的能量分解分析》（//www.umsyar.com/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分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）以及《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（//www.umsyar.com/667）中的综述详细介绍的笔者提出的IGMH方法用于图形化展现片段间相互作用效果极佳，已被广为使用。IGMH用于展现自定义片段间的pi-pi作用、判断pi-pi作用是否存在尤为好用，因此在此给出一些简单例子予以体现。其方法的原理、细节、计算操作见上面的文章。如果你不知道pi-pi作用可能出现在哪，不方便定义片段，则应当用《使用IRI方法图形化考察化学体系中的化学键和弱相互作用》（//www.umsyar.com/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碳环的独特吸附行为的 、波函数分析与分子动力学研究》（//www.umsyar.com/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作用形式。如果读者对碳环类体系感兴趣，非常建议进入//www.umsyar.com/carbon_ring.html观看笔者发表的所有碳环相关的理论研究和对应的深入浅出的评述文章。

上图中还对利用Multiwfn的基于力场的能量分解（EDA-FF）功能计算了各个原子对分子间相互作用的贡献并对原子进行了着色，展现出了更丰富的信息。这种分析见《使用Multiwfn做基于分子力场的能量分解分析》（//www.umsyar.com/442）。

在《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》（//www.umsyar.com/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的原子电荷就能体现，原子电荷的相关知识见《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！》//www.umsyar.com/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堆积二聚体的相互作用能，作用能越负说明作用越强。但如果要讨论二聚体的热力学稳定性，则需要计算结合自由能，溶剂环境中还得考虑溶剂模型。相关知识参见《谈谈分子间结合能的构成以及分解分析思想》（//www.umsyar.com/733）。复合物AB在特定结构下，A和B的相互作用能的最常规的计算方法就是用E(AB)-E(A)-E(B)方式手动计算，做sobEDAw能量分解（//www.umsyar.com/685）时也会顺带给你相互作用能。

以上述方法算相互作用能的一个问题是，如果两个分子之间不仅仅有pi-pi作用，还有其它作用（如氢键），那么得到的只是总相互作用能。如果想只得到pi-pi作用能，有几个办法可以用：
(1)用《使用Multiwfn做基于分子力场的能量分解分析》（//www.umsyar.com/442）介绍的EDA-FF能量分解方法，将两个分子的pi作用区域定义为两个片段，让Multiwfn给出基于分子力场的这两个部分的相互作用能，取其中的范德华作用能（即色散作用和交换-互斥作用之和）。如果你只需要色散部分，还有另一种做法，见《使用Multiwfn图形化展现原子对色散能的贡献以及色散密度》（//www.umsyar.com/705）。这两种做法还都可以给出具体原子产生的贡献，并可以对原子进行着色以便通过图像直观展现和考察
(2)如果分子间同时有pi-pi作用和氢键，且其它作用可忽略不计，可以按《透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍》（//www.umsyar.com/513）介绍的方法使用Multiwfn做拓扑分析估计出氢键作用能，从总相互作用能中扣掉之作为pi-pi作用能
(3)对体系进行恰当改造，基本保留每个分子参与pi-pi作用的部分，然后将总相互作用能近似当成pi-pi作用能。

如果是分子内的pi-pi作用能的估计，还可以参考《计算分子内氢键键能的几种方法》（//www.umsyar.com/522）里的说明举一反三处理。

在特定条件下，pi-pi作用强度和IGMH的等值面面积有密切的正相关性。《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》（//www.umsyar.com/727）介绍的文章中给出了下图，是不同原子数的碳环与富勒烯之间的相互作用能和描述pi-pi作用的IGMH等值面面积的关系，可见有挺理想的线性关系。通过这种关系，可以直接根据pi-pi作用区域的等值面面积估计对应的pi-pi相互作用能。这种面积的计算方法见《计算IGMH等值面的面积和体积的方法》（//www.umsyar.com/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色散校正？” 》（//www.umsyar.com/413）、《DFT-D色散校正的使用》（//www.umsyar.com/210）、《谈谈 研究中什么时候用B3LYP泛函优化几何结构是适当的》（//www.umsyar.com/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相互作用！》（//www.umsyar.com/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电子结构特征》（//www.umsyar.com/432）介绍的方法直接把pi电子密度等值面图画出来，一目了然
(2)相互作用的pi电子之间离得不太远。比如明显超过5埃的直接就可以忽略了
(3)用Multiwfn做IGMH分析（如果得不到波函数文件或体系太大难以算得动的话，可以改用极便宜且只依赖于原子坐标的mIGM），在等值面数值调到诸如0.003 a.u.这样较小的值时，在预期出现pi-pi作用的区域能明显看到等值面
(4)两部分之间的Mayer键级或模糊键级或离域化指数非常小（远小于0.1），故而轨道相互作用可基本忽略。注：如《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）所述，Multiwfn在主功能9里计算键级之前可以用选项-1定义两个片段，之后做键级计算时会给出两个部分之间的总键级，即片段间每一对原子的键级的总和
(5)使用sobEDAw等能量分解方法分析，色散项应当占所有吸引项总和的大部分，轨道相互作用项应当占比微小

文/Sobereva@北京科音  2025-Jan-1

1 前言

有人在思想家公社QQ群（//www.umsyar.com/QQrule.html）里问怎么计算分子轨道的体积，并且贴了一张轨道等值面图。只要计算出这个等值面里包围的体积就可以达到这个目的（虽然也可以有其它方式的定义，但没这么简单直观）。这里我介绍一下怎么利用Multiwfn程序的定量分子表面分析功能来实现这个目的。定量分子表面分析功能之前在《使用Multiwfn的定量分子表面分析功能预测反应位点、分析分子间相互作用》（//www.umsyar.com/159）有简要介绍，结合Multiwfn的极度灵活的设计，这个功能还能做很多其它的分析，正如本文所展示的。

如果读者不熟悉Multiwfn，应当阅读《Multiwfn FAQ》（//www.umsyar.com/452）、《Multiwfn入门tips》（//www.umsyar.com/167）、《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）。Multiwfn可以在//www.umsyar.com/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的命令行方式运行和批量运行的方法》（//www.umsyar.com/612）里专门讲了怎么用Multiwfn批处理计算。这里提供一个脚本，可以便捷地用前文的方法批量计算特定序号范围的轨道的体积。//www.umsyar.com/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)衡量轨道的空间离域程度》（//www.umsyar.com/525），里面用我提出的ODI指数考察了与本文相同的D-pi-A体系的各个轨道的离域程度，并且绘制了ODI与轨道序号之间的柱形图。将上图与ODI的柱形图相对比，可以看到两个图有明显的互补性，即ODI越低（体现轨道越离域、越能同时覆盖很多原子），轨道体积倾向于越大。如果你的目的就是想衡量轨道的离域程度，ODI相对更严格，毕竟不依赖于轨道等值面数值的选取（轨道等值面数值设得很大的话，轨道等值面就完全看不到了，体积为0）。而用轨道体积来说事的好处是可以直接与等值面图相对应，结合图像讨论非常清楚直观。

4 其它

本文的做法不限于计算考察分子轨道的体积，对任意轨道都可以考察，如定域化轨道、NTO轨道、双正交化轨道、NAdO轨道、AdNDP轨道等等，载入含有相应轨道的波函数文件即可，产生方法在下文说了。
Multiwfn的轨道定域化功能的使用以及与NBO、AdNDP分析的对比
//www.umsyar.com/380
用于非限制性开壳层波函数的双正交化方法的原理与应用
//www.umsyar.com/448
使用键级密度(BOD)和自然适应性轨道(NAdO)图形化研究化学键
//www.umsyar.com/535
使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键
//www.umsyar.com/138

文/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程序（//www.umsyar.com/multiwfn）做波函数分析讨论芳香性，是波函数分析方法和Multiwfn程序的很好且十分典型的应用范例。Multiwfn能做的芳香性分析在《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）里有简略介绍，在北京科音 波函数分析与Multiwfn程序培训班（http://www.keinsci.com/WFN）里有多达160页以上的幻灯片极为完整、全面、系统地讲授芳香性的各种概念和全部分析手段，非常推荐参加！

与本文的研究有密切相关的是《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（//www.umsyar.com/696）介绍的Inorg. Chem., 62, 19986 (2023)一文，十分建议看完下文后认真阅读。

2 苯、无机苯和carborazine的结构和稳定性

文中首先用ωB97XD/def2-TZVP对上述三个体系做了几何优化，这是相当稳的级别，也被用于//www.umsyar.com/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氮环的特征！》（//www.umsyar.com/725）介绍的文章里通过AIMD直接证明了18氮环非常不稳定。本文对carborazine和borazine用ORCA程序在ωB97XD/6-311G*级别下做了50 ps AIMD，由于AIMD的昂贵耗时使得模拟时间很有限，文中故意用了1000 K的较高温度令势能面的采样尽量充分。文中还确认了当前温度下的AIMD过程的平均总能量是高于极小点结构下的ZPE的，因此当前的模拟不会低估实际的振动基态能量。这种模拟在北京科音高级 培训班（http://www.keinsci.com/KAQC）里有超详细的讲解，在《使用ORCA做从头算动力学(AIMD)的简单例子》（//www.umsyar.com/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的电荷分布

原子电荷的概念在《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！》（//www.umsyar.com/714）介绍的笔者的文章里有充分介绍，这是定量衡量体系电荷分布最常用的指标。本文通过流行的ADCH方法以及NPA、Hirshfeld和Hirshfeld-I、CHELPG方法都对前述三个体系计算了原子电荷，结果如下所示。可见虽然不同方法计算的原子电荷存在定量差异，但结论都是共通的，即硼和氮由于电负性一个很小一个极大，导致分别带显著的正电荷和负电荷，而在carborazine中它们的差异远小于在无机苯中，这体现出在硼和氮之间插入的碳可以极大地均衡体系中的电荷分布。另外，carborazine中的碳的原子电荷的大小很小，体现出旁边的氮和硼并没有对碳的净电荷带来明显影响。

《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/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可以非常便利地计算轨道成份，见《谈谈轨道成份的计算方法》（//www.umsyar.com/131）。使用Mulliken方法计算的上面图中展示的占据的pi轨道的成份如下所示。可见无机苯的所有pi轨道全都是由氮原子贡献所主导的，这也必定导致无机苯的pi电子几乎都集中且定域在氮原子上，因此在这点上它就不可能有显著的芳香性。而虽然carborazine的HOMO-3和HOMO-5也是以氮原子的贡献为主（因为氮的pz原子轨道的能量比碳和硼的都要低，故偏向于对能量较低的分子轨道贡献较多），但还不至于起到绝对的主导性，而其HOMO的主要贡献则来自硼和碳。由于carborazine的硼、碳、氮都同时明显参与了pi占据轨道，它势必有比无机苯强得多的六中心pi电子离域。

5 苯、无机苯和carborazine的成键特征

键级的概念在《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）中做了简要介绍。文中用Multiwfn对上述三个体系计算了Mayer键级，它体现两个原子间等效共享的电子对儿数，结果如下表左侧所示，并且还按照《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/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键临界点处电子密度拉普拉斯值符号判断相互作用类型失败原因的图形分析》（//www.umsyar.com/161）中专门说过类似情况。对B-N键，如前所示Mayer键级显著大于1，而且BCP处的能量密度明显为负，就已经充分体现出这属于共价键了。文中更进一步在补充材料里还用Multiwfn绘制了电子定域化函数（ELF）和变形密度图，分别如下图左侧和右侧所示，可见B-N之间有显著的电子高定域性区域，并且从孤立状态变成分子过程中在B-N键上有显著的电子密度增加（图中红色实线为正值部分），这都进一步证明B-N键本质上是共价键，只不过极性很强而已。关于变形密度的知识见《使用Multiwfn作电子密度差图》（//www.umsyar.com/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电子结构特征》（//www.umsyar.com/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、多中心键级研究多中心键》（//www.umsyar.com/138）和《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）。结果如下所示，可见六中心电子离域性是苯 > carborazine > 无机苯，差异非常鲜明，这和上述其它分析给出的结论非常一致。从MCBO的数值大小来说，carborazine具有苯的大约一半的芳香性，而无机苯则连carborazine芳香性的一半都不到。

7 苯、无机苯和carborazine的磁屏蔽特征

有大量的方法从体系对外磁场的屏蔽角度衡量芳香性，参考笔者写过的《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）、《使用Multiwfn绘制一维NICS曲线并通过积分衡量芳香性》（//www.umsyar.com/681）、《使用Multiwfn计算FiPC-NICS芳香性指数》（//www.umsyar.com/724）、《通过Multiwfn绘制等化学屏蔽表面(ICSS)研究芳香性》（//www.umsyar.com/216）、《使用Multiwfn巨方便地绘制二维NICS平面图考察芳香性》（//www.umsyar.com/682），在《 波函数分析与Multiwfn程序培训班》（http://www.keinsci.com/WFN）里有更全面系统的讲授。本文从磁屏蔽角度对苯、无机苯和carborazine的芳香性特征进行了定量的考察以及直观的展现。

NICS及其变体在定量衡量芳香性强弱方面极为常用，文中对苯、无机苯和carborazine计算了NICS最理想的形式之一NICS(1)ZZ，并通过《基于Gaussian的NMR=CSGT任务得到各个轨道对NICS贡献的方法》（//www.umsyar.com/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绘制磁感应电流图》（//www.umsyar.com/294）。本文对苯、无机苯和carborazine分别绘制了总的、sigma电子的和pi电子的AICD图，如下所示，从AICD(pi)图中箭头标注的感生电流上可以清楚地看出carborazine形成了和苯一样的贯穿整个六元环的显著的感生电流，而无机苯则完全不具备这样的特征。AICD图进一步有力地体现出，carborazine和苯一样都属于芳香性物质，而无机苯是近乎非芳香性的，这和ICSS分析的结论一致，彼此相互印证。

此文还使用《使用SYSMOIC程序绘制磁感生电流图和计算键电流强度》（//www.umsyar.com/702）中介绍的做法用CTOCD-DZ2方法绘制了分子平面上方1 Bohr处的感生电流图，如下所示，灰、红、蓝分别对应碳、硼、氮原子。由图可见CTOCD-DZ2方法展现的感生电流图和pi电子的AICD图的现象很类似，而这种平面图看得更清楚一些。此图体现carborazine分子平面上方pi区域的感生电流的强度虽不像苯那么均匀，在硼上方相对弱一些，但终究还是形成了明显的全局感生电流。而无机苯的感生电流则主要在氮原子周围打转，体现出pi电子主要都是定域在各个氮的局部原子空间里。

为了便于定量对比，此文还计算了六元环上各种键的电流强度（BCS）值，即在键中心位置垂直于键的平面上的感生电流的积分，结果如上图标注的数值所示，单位为nA/T。可见carborazine的感生电流强度与苯相差并不悬殊，而无机苯则相当微弱。这从定量角度上进一步验证了前面基于肉眼观看感生电流图所做出的结论。

值得一提的是，在《深入理解分子轨道对磁感生电流的贡献》（//www.umsyar.com/703）中我给了计算感生电流的具体公式，假定只有HOMO-LUMO跃迁对感生电流有显著贡献，那么感生电流大小与HOMO-LUMO gap之间会有反比关系。在前面说过，carborazine的HOMO-LUMO gap明显小于苯和无机苯，因此基于磁属性和感生电流得到的关于carborazine的芳香性的强度会有一定程度的高估。但即便如此也不影响它具有明显芳香性的结论，而且前面从多中心键级、LOL-pi等分析角度也全都体现了carborazine有远强于无机苯的芳香性。

9 carborazine中的碳原子对硼和氮原子桥联的重要意义

前面从各个角度极为充分、全面、严格地论证了carborazine具有无机苯不具备的显著的芳香性，而二者的差异仅在于前者给后者的氮和硼之间引入了碳原子起到桥联作用。为什么碳原子的桥联对于提升芳香性有如此关键性的影响？此文将Multiwfn、NBO、VMD相结合绘制了前述三个体系的环上的各原子的对应2pz原子轨道的预正交自然原子轨道（PNAO）等值面的叠加图，过程可参照《使用Multiwfn绘制NBO及相关轨道》（//www.umsyar.com/134）举一反三，如下所示。可见2pz的空间延展程度是硼 > 碳 > 氮，这和电负性增大的顺序相一致。氮和硼的轨道重叠颇差，电子明显难以离域过去。而给二者之间插入的碳原子，增加离域路径上的原子轨道重叠程度、降低相连原子的电负性不平衡性，明显可以起到了帮助电子离域过去的“中介”，从而使得carborazine能够具有无机苯不具备的显著芳香性。之前笔者在《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（//www.umsyar.com/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@北京科音   2024-Oct-28

0 前言

18碳环以其独特的几何和电子结构，自从2019年在凝聚相观测到后引发了化学界的巨大关注。与之相关的分子间相互作用文章已有不少，例如《全面探究18碳环独特的分子间相互作用与pi-pi堆积特征》（//www.umsyar.com/572）介绍的Carbon, 171, 514-523 (2021)、《8字形双环分子对18碳环的独特吸附行为的 、波函数分析与分子动力学研究》（//www.umsyar.com/674）介绍的Phys. Chem. Chem. Phys., 25, 16707 (2023)、《理论设计新颖的基于18碳环构成的双马达超分子体系》（//www.umsyar.com/684）介绍的Chem. Commun., 59, 9770 (2023)，等等。笔者迄今对18碳环及其衍生物的研究论文汇总和各种相关博文见//www.umsyar.com/carbon_ring.html（不断更新）。

碳单环体系和富勒烯体系是碳元素的两种关键的同素异形体，形状截然不同，一个环形一个球形，而且由于二者都有pi共轭特征，理应存在pi-pi相互作用。这显著的差异性和明显的共性无疑使得碳环与富勒烯的相互作用非常值得进行细致、充分的探索。近期，北京科音自然科学研究中心（http://www.keinsci.com）的卢天和江苏科技大学的刘泽玉在知名的Chem. Eur. J.（欧洲化学）期刊上发表了极为全面、系统、透彻的研究不同尺寸碳环（C18到C36）与不同数目C60富勒烯相互作用的理论研究文章，充分揭示了它们之间独特的相互作用，包括强度、结构和本质，内容很新颖有趣，非常欢迎阅览和引用：

Zeyu Liu, Tian Lu*, Theoretical Insight into Complexation Between Cyclocarbons and C60 Fullerene, Chem. Eur. J., 30, e202402227 (2024) DOI: 10.1002/chem.202402227
可以通过此链接免费在线阅览：https://onlinelibrary.wiley.com/share/author/PJDSJZN8IMBAVDCCBXQS?target=10.1002/chem.202402227

此文不仅研究碳环与富勒烯相互作用本身，文章在分析思想和方法学方面对于其它分子间相互作用的研究也颇有借鉴意义，十分推荐对这类研究问题感兴趣的读者阅读。此文充分运用了强大的Multiwfn（主页//www.umsyar.com/multiwfn）波函数分析程序提供的丰富的弱相互作用分析功能，包括IGMH、sobEDA、范德华势等，是Multiwfn研究弱相互作用问题的很好的范例。Multiwfn支持的弱相互作用分析功能在《Multiwfn支持的弱相互作用的分析方法概览》（//www.umsyar.com/252）有概述，在 波函数分析与Multiwfn程序培训班（http://www.keinsci.com/workshop/WFN_content.html）的“弱相互作用的分析”部分有最全面透彻的讲解和演示。

下面将对上述Chem. Eur. J.文章的主要内容进行深入浅出的介绍，便于读者更容易理解文章的研究结果，同时额外附上许多分析和计算细节以帮助读者能够将文中的研究手段举一反三运用到自己的研究中。原文里还有很多细节信息和讨论，故请阅读下文后阅读原文。

1 碳环-富勒烯二聚体

文章首先考察了碳环与富勒烯形成的二聚体结构。文中用Gaussian 16在ωB97XD/6-311G*级别下优化了碳环Cn (n=18, 20, 22, 24, 26, 30, 32, 34, 36)与富勒烯形成的二聚体C60@Cn，并确认了无虚频。ωB97XD泛函曾用于//www.umsyar.com/carbon_ring.html里罗列的笔者的各种18碳环的研究中，不仅可以正确描述18碳环的几何结构（见Carbon, 165, 468 (2020)中的方法对比测试），也能合理描述弱相互作用，故很适合用于优化涉及碳环的复合物。得到的二聚体中比较有代表性的5个如下所示。可以看到随着碳环尺寸的增大，富勒烯逐渐陷入碳环中，到了C60@C34的时候富勒烯正好不偏不倚精确嵌入在碳环的正中央，富勒烯的中心和碳环的中心正好重合，这是一个完美的纳米土星（nano-saturn）结构。而当碳环尺寸进一步增大到C36，为了最大化分子间相互作用，富勒烯的中心自发偏离了碳环的中心。

笔者提出的IGMH是目前非常流行、效果理想的可视化片段间相互作用的分析方法，见《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）和《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（//www.umsyar.com/667）介绍的综述。上图的等值面是IGMH方法定义的sign(λ2)ρ着色的δg_inter等值面，描述富勒烯和碳环这两个片段间的相互作用。图中非常清楚直观地展现出了碳环-富勒烯之间的相互作用区域。等值面几乎都是绿色，说明相互作用区域的电子密度数值非常低，pi-pi相互作用区域普遍具有这种特点，再加上碳环和富勒烯的接触面也正是pi电子分布的区域相互接触，无疑富勒烯和碳环之间的结合是典型的pi-pi相互作用所驱动的。

为了严格考察碳环和富勒烯的结合强度，基于Gaussian优化的几何结构，文中用ORCA程序精确计算了不同碳环与富勒烯形成二聚体对应的结合能，即E_bind = E(C60@Cn) - E(C60) - E(Cn)，其中E是电子能量，每个能量都在各自分别优化的结构下得到。计算级别用的是ωB97M-V/def2-QZVPP，ωB97M-V不仅像ωB97XD一样是长程极限HF成份为100%的范围分离泛函因而能够合理描述碳环的电子结构，而且大量测试都体现了其计算弱相互作用能相当优秀，在《简谈 计算中DFT泛函的选择》（//www.umsyar.com/272）中以及北京科音中级 培训班（http://www.keinsci.com/workshop/KBQC_content.html）里的“弱相互作用的计算与相关问题”专题部分也专门说过这一点。def2-QZVPP是相当高质量的基组，算分子间相互作用能的BSSE问题小到可以忽略，在ORCA里利用RIJCOSX加速技术计算像C60@C36这样约100个非氢原子体系的单点能完全算得动。此外，为了从热力学角度严格考察C60@Cn二聚体形成的可能性，文中还计算了标况下的结合自由能，其中对电子能量的自由能热校正量通过《使用Shermo结合 程序方便地计算分子的各种热力学数据》（//www.umsyar.com/552）介绍的Shermo程序基于ωB97XD/6-311G*振动分析的输出文件得到。如552博文所述，像这种含有大量很低频的体系的自由能的计算一定要用Shermo支持的quasi-RRHO模型得到的自由能热校正量，而切勿直接用Gaussian等程序基于RRHO模型给出的，否则低频模式对熵的贡献可能被高估得离谱，导致结合自由能误差很大。另外，为了考察溶剂效应对结合的影响，文章还使用Gaussian的SMD溶剂模型按照《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》（//www.umsyar.com/327）说的做法获得了水环境下的自由能，并求差得到了水中的结合自由能。计算结果汇总如下

由上图可见，从C18到C32，随着碳环尺寸的增大，它与富勒烯的结合强度近乎严格地线性增强，这是由于二者间的接触面积逐渐增大，这点从上文的IGMH图的等值面面积的变化就可以清晰地反映出来。而从C34变成C36，结合强度反倒减弱了，这正如IGMH图所示的，C36有的地方与富勒烯离得相对较远、作用较弱，导致δg_inter等值面也没有在相应区域出现。

上图体还现出标况气相的结合自由能和基于电子能量算的结合能的变化趋势基本一致，但前者没有后者那么负，这是因为结合过程中熵减小对结合造成了严重不利的影响。对比上图中的红线和蓝线可见水溶剂环境下的结合自由能比气相下的更负几kcal/mol，说明极性溶剂会促使碳环与富勒烯的结合，这本质上就是疏水效应，众所周知疏水效应会驱动非极性物质间在水中的结合。

很值得一提的是，文中发现碳环-富勒烯之间的δg_inter=0.002 a.u.的等值面面积与它们的结合能有颇好的线性关系，如下所示。因此IGMH图像的等值面面积在某些情况下可以作为解释或预测相互作用能的很好的描述符，这一点很值得在未来进一步探索。这种面积的计算方法见《计算IGMH等值面的面积和体积的方法》（//www.umsyar.com/738）。

为了考察温度对碳环-富勒烯结合的影响，以及探索结合的临界温度，文中利用Shermo程序做了自由能随温度的扫描，并进而得到了结合自由能与温度的关系，如下所示。可见温度越高越不利于结合，在气相常压下，C18碳环与富勒烯之间从366 K开始就无法结合了，而C34和富勒烯之间的二聚体在高达605 K以下都是可以形成的。C36的情况处于二者之间。

碳环与富勒烯的复合势必造成单体的结构变化，由于单体结构相对于极小点结构改变导致单体的电子能量在复合过程中的升高称为变形能。文中补充材料里的表S1给出了不同碳环与富勒烯结合时的碳环的变形能，发现C18到C34的变形能都<=0.3 kcal/mol，可以忽略不计，而C36的变形能则达到了1.0 kcal/mol，这一方面是因为这样较大的碳环有足够显著的柔性，另一方面在于它与富勒烯之间的不对称的相互作用（如全面IGMH图所展示的），这造成了这种碳环在复合物中相对于圆形的极小点结构的不可忽略的形变。下图是按照《在VMD中计算RMSD衡量两个结构间的差异以及叠合两个结构》（//www.umsyar.com/290）的方法绘制的C36在孤立状态（红线）和与富勒烯结合状态（蓝线）的叠加图，可见富勒烯的存在诱使C36碳环被拉长了。

碳环与富勒烯的作用虽然形式上属于“弱相互作用”类别，但结合强度真不是一般的强，而是相当的强！如《透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍》（//www.umsyar.com/513）所列举的，像是水二聚体这样的一般强度氢键的结合能才-5 kcal/mol左右，而本文研究的与富勒烯结合最弱的18碳环，它与富勒烯的结合能都达到了-14.4 kcal/mol，将近三个普通氢键的强度。而且这比《全面探究18碳环独特的分子间相互作用与pi-pi堆积特征》（//www.umsyar.com/572）提到的曾经笔者研究的两个18碳环间的结合能-9.2 kcal/mol强得多，也比J. Chem. Theory Comput., 13, 274 (2017)高精度计算的两个富勒烯之间的结合能-8.3 kcal/mol强得多！

2023年笔者提出的能量分解方法sobEDA和sobEDAw的详细介绍见《使用sobEDA和sobEDAw方法做非常准确、快速、方便、普适的能量分解分析》（//www.umsyar.com/685），其中sobEDAw用于考察碳环+富勒烯的相互作用成份极为合适，速度足够快而且结果也很准确。sobEDAw只对部分泛函拟合了参数，由于碳环体系需要HF成份较高的泛函才能正确描述，因此文中使用了HF成份达到50%的BHandHLYP泛函结合DFT-D3(BJ)色散校正和6-311+G(2d,p)基组做碳环-富勒烯之间的sobEDAw能量分解，并且考虑了counterpoise校正，结果是总相互作用能为-12.0 kcal/mol，其中静电贡献-8.4、交换互斥贡献25.6、轨道相互作用贡献-2.4、色散作用贡献-26.8 kcal/mol。可见色散作用对碳环-富勒烯之间的吸引作用起绝对主导效果，静电作用相对次要，而轨道相互作用可忽略不计。这完全符合pi-pi堆积的典型特征，在《全面探究18碳环独特的分子间相互作用与pi-pi堆积特征》（//www.umsyar.com/572）介绍的Carbon, 171, 514 (2021)中使用sSAPT0/jun-cc-pVDZ级别对18碳环的pi-pi作用形成的二聚体做能量分解基本上也是这个情况。

此文还使用Multiwfn通过笔者提出的ADCH原子电荷计算方法，以及很常用的Mulliken方法，计算了C60@C18中18碳环部分的片段电荷。如果对原子电荷缺乏了解的话建议阅读《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！》（//www.umsyar.com/714）里提到的笔者的原子电荷综述。这两种方法给出的18碳环的电荷分别为0.015和-0.027，都十分接近于0，体现出碳环与富勒烯之间的基态的电荷转移效应可以忽略不计。文中也用Multiwfn的主功能9的选项-1将18碳环和富勒烯分别定义成片段1和片段2，然后计算了它们之间的总Mayer键级，数值仅为可忽略的0.041，即基本没有共享电子作用，再次体现出两个分子间完全是非共价相互作用。如果读者对键级不了解，建议看《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）中键级的相关部分。

文中发现18碳环和富勒烯实际上有两种极小点结构（见原文图S1），次低的比最低的能量高0.9 kcal/mol，它们之间互变的过渡态的虚频为9.8i cm-1，振动模式如下所示，可见明显对应于富勒烯相对于18碳环的旋转。在ωB97M-V/def2-QZVPP级别下计算出的正向和逆向势垒分别为0.90和0.04 kcal/mol，极低的势垒体现出复合物中富勒烯可以非常自由地旋转。实际上pi-pi堆积作用的体系之间普遍容易发生滑移（如18碳环二聚体，在//www.umsyar.com/572里做了动力学模拟直接体现了这一点），因为滑移造成的能量变化很小。

2 碳环-富勒烯 2:1三聚体

文章系统考察了一个富勒烯结合两个碳环（C18至C34）形成的三聚体的结构，优化得到的结果如下图所示。这种体系中既有碳环与富勒烯之间的相互作用，也有两个碳环之间的相互作用，为了能直观区分，这两类相互作用的δg_inter等值面分别用绿色和青色着色，这样的IGMH图真是巨直观！可以看到，随着碳环的增大，两个碳环间的夹角逐渐减小，富勒烯越来越多地被两个碳环“吃”了进去，富勒烯始终与碳环保持着紧密的接触和充分的pi-pi作用。当碳环大到C34时，两个碳环之间完全平行并形成了全面的pi-pi相互作用，而且和//www.umsyar.com/572提到的18碳环二聚体一样是较长的C-C键对着较短的C-C键，此时富勒烯已精确嵌入到了两个碳环的正中间。

ωB97M-V/def2-QZVPP计算远超100个非氢原子的体系时过于昂贵，因此文中在计算三聚体和更多聚体时在ORCA里使用的是ωB97M-V/def2-TZVP结合gCP方式的经验性的BSSE校正，精度比ωB97M-V/def2-QZVPP差不了多少而耗时低得多得多。使用ωB97M-V/def2-TZVP+gCP，文中计算了富勒烯结合第一个碳环的结合能ΔE_1_bind，以及在此基础上再结合第二个碳环的结合能ΔE_2_bind（即C60@Cn与Cn形成C60@2Cn的结合能），总的结合能（三聚体结合能）为ΔE_1_bind与ΔE_2_bind之和。结果如下图(a)所示。可见随着碳环增大，结合能越来越负、结合作用越来越强，而且ΔE_2_bind比ΔE_1_bind明显更负，这体现出多个碳环与富勒烯结合的非常明显的协同作用。如上面的IGMH图所示，已有的一个碳环可以和另一个碳环形成pi-pi作用，而且碳环越大时pi-pi作用越显著（青色等值面越来越大），这是为什么碳环越大时协同效应越强。

为了更进一步考察三聚体结构中不同片段之间的相互作用强度，此文还基于优化后的三聚体结构计算了片段间的相互作用能，如上图(b)所示，两个碳环间的相互作用能对应红线，两个碳环与富勒烯之间的相互作用能对应蓝线。可见随着碳环增大，碳环间的相互作用逐渐增强，特别是两个C34碳环之间的相互作用格外强，这是因为如前面的IGMH所示，它们能形成完整的遍及整个碳环的pi-pi作用，而不再只是局部区域的pi-pi作用。还可见从C18到C30之间，碳环越大，由于可以和富勒烯相互作用的原子越多，碳环-富勒烯的相互作用越强，但是从C30到C34这种相互作用反倒稍微变弱了，这是因为在C60@2C34中C34与富勒烯的接触不像C60@2C30中C30与富勒烯的接触那么充分，这一点从IGMH图上就可以明显看出来，即C60@2C34中绿色的两条环状等值面明显比C60@2C30的窄很多。这体现出恰当利用IGMH方法可视化分子间相互作用，甚至可以把不同体系间轻微的差异性也给明确揭示出来。《直观解释分子间相互作用如何影响不对称催化：Nature Chemistry上一个很好的IGMH分析范例》（//www.umsyar.com/700）里的例子也同样体现了这一点。

各个片段在结合成C60@2Cn三聚体过程中的总变形能如上图的黑线所示，由于富勒烯的刚性很强，因此变形能基本都来自于碳环的变形。可见随着碳环的增大，变形能也逐渐上升，这在于越大的碳环柔性越强。从前面的结构图甚至肉眼都能明显看出来C30碳环在三聚体结构中有明显的弯曲。而C34的变形能则小于C30，这是因为C60@2C34复合物中的两个C34基本是圆形、平面的结构，和孤立状态结构特征相差较小，而且此结构中C34与富勒烯的相互作用没有C30的那么强。

3 碳环-富勒烯 1:2三聚体

下面再来看一个碳环与两个富勒烯结合成2C60@Cn（n=18到30）三聚体的情况。优化后的复合物结构和IGMH图如下所示，绿色等值面展现富勒烯与碳环之间的相互作用，黄色等值面展现两个富勒烯之间的相互作用。标注的d_min对应两个富勒烯之间最近原子距离（可以将结构文件载入Multiwfn，进入主功能100的主功能21，输入dist，然后依次输入两个片段里的原子序号得到此值）。由图可见，从C18到C26碳环，碳环越大，由于对两个富勒烯之间的接触阻碍越小，富勒烯之间的距离越近、富勒烯之间的相互作用越显著。碳环从C22变成C26时，富勒烯与碳环之间的IGMH等值面明显变窄，体现相互作用明显变弱，和富勒烯之间的相互作用强度的变化形成了此消彼长的关系。对于2C60@C30的情况，碳环与两个富勒烯的相互作用不再对称，上面的富勒烯歪到一侧去了，与碳环的相互作用显著小于下面的富勒烯，因此可以预料2C60@C30的稳定性明显不及2C60@C26。

下图(a)给出了2C60@Cn的三聚体总结合能（黑线），以及结合第一个和第二个富勒烯时候的结合能（红线和蓝线）。可见碳环越大，对第一个富勒烯的结合越强，这在前文的碳环-富勒烯二聚体的研究中已经说过。由于已有的富勒烯会对第二个结合的富勒烯产生吸引作用，对于2C60@C18、2C60@C22、2C60@C26，第二个富勒烯的结合能比第一个富勒烯的结合能更负，即两个富勒烯与C18、C22、C26碳环的结合有明显的协同效应。而对于2C60@C30，由于第一个富勒烯严重妨碍第二个富勒烯的结合，导致第二个富勒烯只能歪斜着与C30碳环发生较弱的相互作用，同时第二个富勒烯结合时还稍微削弱了第一个富勒烯与碳环的相互作用，因此第二个富勒烯的结合能只有第一个富勒烯结合能的一半。

上图(b)更进一步考察了2C60@Cn的三聚体中各个片段间的相互作用能。可见由于碳环越大、越不给富勒烯之间的接触碍事，图中红线对应的富勒烯之间的相互作用能是随着碳环的增大而越来越负的。图中蓝线体现C22和C26与两个富勒烯的相互作用最强，而C30由于只能松散地与其中一个富勒烯相互作用，因此它与两个富勒烯的总相互作用更弱一些。

根据以上信息，文中提出了新颖的“分子胶水”的概念。大小恰到好处的碳环，如C22，可以把两个富勒烯特别牢固地粘在一起。C22对第二个富勒烯的结合能，也等价于C60@C22与一个富勒烯的结合能，达到了约-42 kcal/mol，这比起两个富勒烯之间的结合能-8.3 kcal/mol强太多了，甚至都轻微超过了C-C单键键能的一半！！！

为了更充分地探究碳环对富勒烯结合起到的胶水作用，此文运用了《谈谈范德华势以及在Multiwfn中的计算、分析和绘制》（//www.umsyar.com/551）中介绍的方法，基于GAFF力场参数绘制了2C60@C22和2C60@C30三聚体结构中碳环产生的范德华势，以碳原子为探针原子，如下所示。图(a)的黄色等值面对应范德华势=-0.8 kcal/mol，富勒烯的原子按照原子所在位置的范德华势着色，颜色越蓝体现富勒烯的相应原子与碳环间的色散吸引作用越强。图(b)是利用VESTA程序基于Multiwfn产生的范德华势cube文件绘制的填色平面图。由图可见C22碳环的范德华势最负的区域从碳环中心向上下两侧延伸出来，没过了每个富勒烯大约1/3的区域，显然C22碳环的这种范德华势分布能恰到好处地把两个富勒烯粘起来，是富勒烯之间最简单且最完美的胶水分子！而碳环也不宜过大，如图中C30的情况，虽然其范德华势很负的区域能覆盖下方的富勒烯将近一半的原子，但对上方的富勒烯的色散吸引则较为有限。尽管如此，C30对第二个富勒烯的结合能-17.1 kcal/mol仍然比两个富勒烯之间的结合能-8.3 kcal/mol负得多。因此哪怕碳环偏大，照样能促进两个富勒烯的结合。

4 一个富勒烯最多能结合多少个18碳环？

从前面给出的两个18碳环与一个富勒烯形成的复合物来看，富勒烯表面实际上还有空间可以结合更多的18碳环。那么一个富勒烯最多能结合多少个18碳环？多少个碳环才能令富勒烯的表面饱和？为了探究这个有趣的问题，文中根据碳环和富勒烯的结构，搭建了6个碳环均匀分布在富勒烯上下左右前后的结构并做了几何优化。由于此体系多达168个碳原子，用ωB97XD/6-311G*优化和振动分析实在过于昂贵，所以这部分的计算对碳环用6-311G，而对富勒烯用6-31G。这么做的合理性一方面在于碳原子对于极化函数的要求远低于杂原子，另一方面是以C60@C18进行测试，6-311G&6-31G优化的结果与6-311G*的结果的差异确实可忽略不计，分别如下图的红线和蓝线所示（原文图S4），可见几乎完全重合。注意18碳环的基组不能再减小到6-31G，这点笔者在《我对一篇存在大量错误的J.Mol.Model.期刊上的18碳环研究文章的comment》（//www.umsyar.com/584）中专门指出过。

优化得到的富勒烯结合6个18碳环的结构如下图(a)所示。为了清楚起见，在VMD程序里通过绘图命令用黄色圆球把18碳环中心位置显示了出来，并且把相邻的黄球用蓝线连了起来，由此可以清楚地看出6个碳环构成了一个包围富勒烯的近乎理想的正八面体空间。下图(b)是IGMH图，依然是富勒烯与碳环之间的作用用绿色展示，碳环之间的作用用青色展示，可见这个C60@6C18结构中每个碳环都与富勒烯产生了非常显著、充分的pi-pi作用，而且每一对相邻碳环之间也产生了特别显著的pi-pi作用，这是何等完美的结构啊！6个18碳环与富勒烯的结构匹配程度完美到令人惊叹！

在ωB97M-V/def2-TZVP+gCP级别下计算的C60@6C18的总结合能高达-101.3 kcal/mol，已经达到了常规化学键的数量级！此体系中平均每个18碳环的结合能是-101.3/6=-16.9 kcal/mol，这比起C60@2C18中18碳环平均结合能-15.0 kcal/mol更负，体现出C60@6C18的形成过程的协同作用真是巨强。这在于此体系中每个18碳环都能同时与周围四个18碳环充分地吸引。C60@6C18在气相标况下的结合自由能为-13.4 kcal/mol，体现出此结构在热力学上很容易自发形成。

基于C60@6C18的结构，文章还大胆地设想了富勒烯与18碳环形成的共晶的结构，示意图如下所示（原文图S4）。此结构里每个碳环、富勒烯都最大程度利用了自己的色散吸引能力与对方结合，因此应当是一个颇为稳定、大概率在未来能被实验合成出来的晶体。

5 总结

本文介绍的Chem. Eur. J., 30, e202402227 (2024)一文通过严谨的 计算并充分运用Multiwfn实现波函数分析，首次预测了各种尺寸碳环与不同数目C60富勒烯形成的复合物的结构，并深入探讨了相互作用强度和本质。此文体现出富勒烯与碳环这两种碳的同素异形体之间通过pi-pi作用表现出极强的亲合性。碳环体系以其特殊的环形几何结构以及独特的两套全局共轭的pi电子，还可以作为分子胶水将两个富勒烯牢牢粘在一起。文章还证明了一个富勒烯最多能结合6个18碳环，而且结合过程中有显著的协同性，结合越多越容易。因此靠富勒烯吸附碳环，或许在未来能成为一种富集富勒烯的手段。文章还预测了富勒烯与18碳环形成共晶的可能，在未来有可能能以这种形态将不稳定的18碳环稳定地储存起来。

此文是通过理论计算研究新颖的分子间复合物的很好的例子，兼具重要的理论意义和实际意义。同时此文也是利用Multiwfn做波函数分析探究弱相互作用的很好的范例，把相互作用特征研究得十分通透，尤其是IGMH方法把富勒烯与碳环的复合物中的弱相互作用展现得超级生动直观、一目了然，此文充分体现出掌握Multiwfn程序做波函数分析对弱相互作用研究的关键性价值！

另外，如果你对pi-pi研究感兴趣，强烈建议阅读《谈谈pi-pi相互作用》（//www.umsyar.com/737）。

文/Sobereva@北京科音  2024-Aug-31

0 前言

2019年首次在凝聚相中发现的18个碳原子相连构成的环状体系18碳环（cyclo[18]carbon）已经广为知晓，笔者对此体系及衍生物陆续做过大量理论研究，成果汇总见//www.umsyar.com/carbon_ring.html。在研究18碳环的过程中有一个问题引发了笔者的好奇心：18个氮原子形成的18氮环（cyclo[18]nitrogen）会是什么结构？能否存在？具有什么性质？之前无论是实验还是理论研究，都完全没有18氮环或其它的较大的纯氮环体系的报道。无疑通过 计算和电子波函数分析回答这些问题非常有理论和实际意义。应ChemPhysChem期刊的邀请，笔者近期在此期刊上发表了18氮环的专题研究，欢迎阅读和引用：

Tian Lu, Theoretical Prediction and Comprehensive Characterization of an all-Nitrogenatomic Ring, Cyclo[18]Nitrogen (N18), ChemPhysChem, 25, e202400377 (2024) DOI: 10.1002/cphc.202400377

此文可以在此免费在线阅读：https://onlinelibrary.wiley.com/share/author/YPJDJ5XPMVT8SD7VDDYQ?target=10.1002/cphc.202400377

下面，笔者将对这篇文章的关键内容进行深入浅出的介绍，同时对研究思想和细节做一些补充说明，以帮助读者更好、更容易地理解这篇文章的工作。此文的研究内容和手段对于理论探究其它特征新奇的物质也很有借鉴意义。

1 18氮环的构型

理论研究一个完全未知的物质，最先要研究的是它的几何结构，然后再说其它的。有的体系有可能存在多个有意义的极小点构型，就都得优化出来并对比能量，以弄清楚哪个是热力学上最稳定的、最需要关注的，以及不同构型的分布比例如何。ωB97XD/def2-TZVP是优化大部分体系很靠谱的级别，笔者之前做过的18碳环及各种衍生物的研究也都是用这个级别优化，因此对18氮环也首先用Gaussian 16在这个级别下做了优化，总共得到了以下结构，三个极小点分别对应C1、D3、D9点群。图中为了令其结构特征看得尽可能清楚，给了俯视图和侧视图，并且让分子最大程度平行于XY平面，并在VMD里根据Z坐标按照色彩刻度条进行了着色。可见，18氮环的极小点结构并不是18碳环那样严格纯平面的，而是弯折的，这和18氮环具有孤对电子而缺乏全局离域的pi电子有关，详见后文。

此文用ORCA在非常精确的DLPNO-CCSD(T1)/cc-pVQZ级别下对上述优化完的结构计算了电子能量，不同构型间相对电子能量如图中的ΔE所示（kcal/mol）。文章还利用Shermo程序（//www.umsyar.com/552）计算了标况下的自由能热校正量并与电子能量相加得到自由能，标况下的相对自由能如图中ΔG所示（kcal/mol）。可见，能量关系是D9>>D3>C1。因此看似结构理想、对称性特别高、像皇冠一样的D9结构实际上无法在现实中出现。根据《根据Boltzmann分布计算分子不同构象所占比例》（//www.umsyar.com/165）介绍的方法，可以算出标况下D3的出现比率也可以忽略不计，C1和它出现比率是158:1。所以本文后面的研究基本都只基于C1结构来做。

为了确认是否有比上述C1能量更低的18氮环的结构，此文还借助molclus程序（http://www.keinsci.com/research/molclus.html）做了构型搜索。具体来说，此文在相对便宜的ωB97XD/def2-SVP级别下控温在较高温度跑了5 ps从头算动力学模拟对势能面进行采样，每隔0.4 ps提取一帧用molclus做批量优化，最后发现所有结构都收敛到了C1，因此可以确信不存在能量更低的结构。注：实测ωB97XD/def2-SVP级别会严重高估18氮环的解离势垒，所以用较高温度跑动力学时也不至于出现解离。

可能有读者想问上述三种极小点结构是怎么获得的。由于18氮环实际长什么样事先完全无法估计，笔者的做法是把18碳环里面每个碳都替换为氮，然后反复进行优化和做消虚频的操作，最终找出无虚频的结构。这个过程中会遇到同时存在许多虚频的结构，虚频大多都会破坏局部对称性，这种情况消虚频常用的做法是按照虚频模式调结构并重新优化，反复如此直到没有任何虚频，这在《Gaussian中几何优化收敛后Freq时出现NO或虚频的原因和解决方法》（//www.umsyar.com/278）中专门说过。取决于以这种方式消虚频的消除顺序，最终得到的无虚频结构可能不同。前述的三种极小点结构是以不同方式调结构消虚频得到的。这么搞不排除遗漏某些极小点的可能，但由于做前述的构型搜索过程中并没有得到其它能量较低构型，因此可以认为至少不存在值得关注的其它能量很低的极小点结构。

前面图中最右边的D3h点群的结构是消除了所有破坏环平面的虚频后的结构，它依然有平面内的虚频。它的能量明显也很高。18氮环不仅能量最低结构不是纯平面的，纯平面结构就连对应的极小点都没有。

毕竟18氮环是个新颖的体系，为了100%确保ωB97XD优化的结构合理，此文还在常用的B3LYP-D3(BJ)、PBE0-D3(BJ)、M06-2X泛函下，以及可靠且昂贵的CCSD级别下也都做了优化（对比见补充材料），结果和ωB97XD的高度一致，而且T1诊断体现出18氮环的极小点结构的多参考特征不强而没必要用多参考方法，因此可以认为本文给出的结构是非常可靠的，笛卡尔坐标在补充材料里提供了。

此文利用《使用Multiwfn计算Bond length/order alternation (BLA/BOA)和考察键长、键级、键角、二面角随键序号的变化》（//www.umsyar.com/501）介绍的方法，利用Multiwfn很便利地得到了整个环上的键长、键角的变化图，如下所示，这从定量层面更进一步展现了18氮环的结构特征。由左图可见，18氮环上的N-N键键长是很明显长-短交替变化的，这也体现出N-N键的强弱是显著交替变化的，这从后文对成键本质的分析上可以了解原因。值得一提的是18碳环也具有类似的明显长-短交替的C-C键。D9结构的较长N-N键比C1结构中的明显更长，一定程度体现出D9的稳定性更弱、更易解离。从键角变化来看，C1和D3极小点结构的键角是在一定范围内波动的，而高对称性的D9中所有键角等同，而且比所有的C1和D3的键角都明显要小。这过小的键角无疑是为了满足其高对称性所致，也必定会因此带来明显的键角张力，这是D9具有很高能量的关键原因。

2 18氮环的动力学稳定性

为了研究18氮环的动力学稳定性，以及它的分解机理，此文对C1结构的一个较长的N-N键按照《详谈使用Gaussian做势能面扫描》（//www.umsyar.com/474）说的方法进行了逐渐拉长的柔性扫描，发现在扫描路径上有个极大点，并且根据它和后面一个点的结构可以确认它适合作为搜索解离过程的过渡态的初猜，果然基于它进一步优化过渡态后得到了对应于18氮环解离的过渡态，并由此按照《在Gaussian中计算IRC的方法和常见问题》（//www.umsyar.com/400）中的做法进一步跑了IRC，而且要求IRC尽可能跑得完整、理想。下图是在ωB97XD/def2-TZVP级别下跑的IRC曲线，以及过渡态和IRC最后一帧的结构。过渡态标记为C1-TS，它具有的唯一虚频对应的振动模式按照《在VMD中绘制Gaussian计算的分子振动矢量的方法》（//www.umsyar.com/567）的做法用黄色箭头绘制了出来。可以看到18氮环的解离不是一次只断一个N-N键，而是三个N-N键同时断开（虚频模式对应它们仨同时显著的伸缩运动），直接解离产物是两个氮气分子和一个14个氮组成的链状结构。

上面的IRC图左端对应C1结构，以它为能量零点，可见这个反应经历了18.9 kcal/mol的势垒。实际上这个势垒数据并不准确。笔者在高精度的DLPNO-CCSD(T1)/cc-pVQZ级别下对过渡态和C1极小点结构计算了电子能量并求差，得到的高精度势垒是8.8 kcal/mol，明显ωB97XD/def2-TZVP严重高估了势垒。再进一步考虑自由能热校正量后，得到的高精度的标况自由能垒是4.4 kcal/mol。将之代入《基于过渡态理论计算反应速率常数的Excel表格》（//www.umsyar.com/310）提供的表格里计算反应速率常数，可知常温下18氮环解离的反应速率常数是3.7E9 /s，对应于半衰期为180 ps。因此，在常温下18氮环的寿命极短，转瞬即解离。然而，在很低温下它则具有一定稳定性，例如70 K的时候自由能垒为5.5 kcal/mol，对应于半衰期为14.8小时（如果再考虑隧道效应会更短）。所以在极低温下产生并检测到18氮环还是有希望的。

3 18氮环的分子动力学行为

前面都是从静态角度研究18氮环。为了了解其动力学行为，以及从动态角度考察解离过程，此文按照《使用ORCA做从头算动力学(AIMD)的简单例子》（//www.umsyar.com/576）介绍的做法做了从头算动力学模拟。如前所述，ωB97XD/def2-TZVP严重高估了解离势垒，势必会导致在动力学过程中严重高估结构的稳定性，因此必须选择一个又足够便宜又能较好符合DLPNO-CCSD(T1)/cc-pVQZ高精度计算的势垒的级别。测试发现B3LYP-D3(BJ)/def2-SVP可以满足这个要求，因此AIMD在此级别下做。

下图左侧是在300K下模拟20 ps过程中18氮环的结构变化，用VMD绘制，每1 ps绘制一次并叠加显示，根据模拟时间按照蓝-白-红着色以区分。可见在算得动的很有限的模拟时间尺度内，18氮环的骨架结构保持得较好，始终都在C1极小点结构附近波动。

上图右侧是分别在较高的500K、750K、1000K下模拟得到的轨迹相对于初始结构的RMSD曲线图。可见三个温度下在前3 ps内都发生了解离，而且温度越高解离时间越早，解离的时刻对应于RMSD突然飙升的位置。从500K曲线上标注的结构可见，在刚解离时，其结构正好和前面优化出的解离过渡态结构C1-TS高度一致，解离产物也和之前跑的IRC的产物端一致，这体现出研究一个反应可以同时在静态和动态两个视角下进行。

温度不光影响解离发生的快慢，还同时影响解离的直接产物。如下图所示，在750K下18氮环不是先解离出两个氮气分子，而是一下子彻底解离成9个氮气分子。这充分体现出18氮环完全解离成氮气在高温下是一瞬间的事，且不经历中间体。

4 18氮环的能量相关属性

此文从能量属性角度对18氮环的特征做了一系列考察。首先此文在DLPNO-CCSD(T1)/cc-pVQZ//ωB97XD/def2-TZVP级别下用Shermo计算了18氮环的标况下的生成自由能和生成焓，分别为685.1和598.7 kcal/mol，前者很大体现出从氮气分子直接合成18氮环在热力学上是极度不利的，后者很正体现出当18氮环分解为氮气分子时会释放巨大热量，因此是极为高能的物质。

值得一提的是18氮环并不是N18最稳定的异构体。有前人研究过N3(N5)3分子，和18氮环一样化学组成都是N18，而N3(N5)3的电子能量和标准自由能分别比18氮环低62.3和56.2 kcal/mol。尽管如此，这并不意味着18氮环就一定不能在实验上得到。一方面如前所述，本身18氮环在很低温下就有一定稳定性，另一方面，利用特殊合成手段，亚稳的环状物质本来也可能得到。例如20个碳连成环状的20碳环如今在实验上已经观测到了，然而高精度理论计算指出碗状的C20比环状的C20的能量低得多。

18氮环可以与氧气反应发生燃烧变成NO2气体，此反应计算出的标准反应焓是-426.6 kcal/mol。这体现出若18氮环固体能制备出来，它可以作为可燃材料，而且由于氮气是空气中含量最多的分子，它在原理上还可以无限产生。

环张力是环状分子重要的特征，在《谈谈如何计算环张力能：以CPP和碳单环体系为例》（//www.umsyar.com/698）里有详细的计算方式的介绍。为了考察18氮环的环张力的大小，此文通过同联反应法计算了其最稳定的C1结构的环张力能，结果为6.3 kcal/mol。相比之下，18碳环的环张力能61.7 kcal/mol以及含有18个碳的[18]环聚乙炔的环张力能70.0 kcal/mol都远大于之，这体现出18氮环环张力能极小，也因此环张力并不是18氮环具有很高能量的原因。之所以它的环张力能相对非常小，在于它由于缺乏整体的pi共轭，使得它的骨架柔性较大，可以自发避免显著形成环张力。

此文还在ωB97XD/aug-cc-pVTZ级别下计算了18氮环的第一垂直电离能、第一垂直电子亲和能、fundamental gap（等于电子硬度）、电子软度和电负性，并与18碳环和氮气分子做了对比，如下表所示。可见18氮环的电离能比氮气分子小得多，主要在于18氮环的较长的N-N键远比氮气分子的N-N键弱得多，因此成键轨道能量较高，自然其电子更容易电离。氮气分子的电子亲和能为负，说明没法再结合额外的电子，这在于它的LUMO轨道是反pi特征，被电子占据后自然会由于削弱成键作用而令能量变得更高。而18氮环的电子亲和能则为正，说明可以再结合电子形成(N18)-阴离子。之所以它还有再额外结合电子的能力，在于它的最低空轨道并不完全对应反pi特征，而是对部分N-N键来说还具有成键轨道特征，因此被电子占据后还能令体系能量变得更低。18氮环的电子软度比氮气分子大很多，说明18氮环的电子整体更容易变形、被极化。与18碳环相比，由VIP和VEA可见18氮环更难失电子而更容易得电子，这也正对应于表中它具有明显更大的电负性，相对来说是更好的电子受体。

5 18氮环的分子光谱

为了令实验化学家在未来可以通过光谱技术检测18氮环，本文理论预测了它的振动光谱和电子光谱。按照《使用Multiwfn绘制红外、拉曼、UV-Vis、ECD、VCD和ROA光谱图》（//www.umsyar.com/224）的做法用Multiwfn模拟的18氮环的基于谐振近似的红外光谱如下所示。由于C1结构的18氮环不像18碳环那样具有高对称性，因此也没有对称禁阻现象，使得它的光谱特征比碳环类体系的丰富的多。笔者发表的不同尺寸碳环的振动光谱和振动行为的专题研究介绍见《揭示各种新奇的碳环体系的振动特征》（//www.umsyar.com/578），里面给出了18碳环的红外光谱。相比之下，18氮环不具备18碳环那样在某些振动模式上具有特别强的红外吸收强度的特征。同样在ωB97XD/def2-TZVP级别下做谐振近似的振动分析，18氮环红外强度最强的振动模式的强度是23.7 km/mol，而18碳环则高达224.4 km/mol。

Multiwfn具有绘制分振动态密度图（partial vibrational density-of-state maps, PVDOS)的功能，由此可以了解不同频率范围的振动模式主要对应的是什么区域或什么特征的运动。下图绘制了总的振动态密度图，并且键伸缩、键角弯曲、二面角扭转三类运动模式的PVDOS也分别给出了。可见在整个波数范围内，三类运动模式之间都有显著的耦合，但导致18氮环骨架大幅变化的二面角扭转模式主要贡献的是中、低频部分，而高频部分更多来自于刚性的键伸缩运动模式。

文中还用Multiwfn模拟了18氮环的UV-Vis光谱图，如下所示，使用的是TDDFT结合PBE0/def2-TZVP。对于不牵扯里德堡激发、电荷转移激发的单重态价层激发的情况，PBE0通常是较好选择，没有明显高估和低谷激发能的趋势，这在《乱谈激发态的计算方法》（//www.umsyar.com/265）里说过。从模拟的谱图上可见18碳环在可见光区的吸收可忽略不计，因此基本上可认为是无色的，起码是对于当前研究的孤立状态来说。PS：想更严格预测可以用《通过 计算和Multiwfn程序预测化学物质的颜色》（//www.umsyar.com/662》里说的方法。

对于18碳环这样基态是单重态的分子，通常最关键的电子激发是S0到S1的激发。为了考察其电子激发本质，文中按照《使用Multiwfn做自然跃迁轨道(NTO)分析》（//www.umsyar.com/377）所述的方法用Multiwfn对S0-S1激发做了NTO分析，并发现此电子激发基本可以由上图中的hole-NTO到electron-NTO的跃迁来描述。从等值面图可以看出，hole-NTO同时具有孤对电子(n)和sigma轨道的特征，而electron-NTO则只有反pi（pi*）轨道特征，因此S0-S1激发可以认为是n-pi*和sigma-pi*的杂化激发。这种情况对于普通有机体系是极其罕见的，因为其sigma占据轨道能量一般都较低，很难在S0-S1激发中牵扯到。

6 18氮环的电子结构

18氮环具有非同寻常的电子结构，这是本文研究的关键重点，必须充分运用波函数分析的手段才能考察。

LOL是一个重要的考察定域化电子出现区域的实空间函数，参见《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）里的简介。文中使用Multiwfn绘制了18氮环的LOL函数的0.55的等值面图，如下所示，上面标注的数字是Multiwfn算的Mayer键级值。由此图可以明显看出18氮环有显著的N-N共价键，同时每个氮上还都有显著的孤对电子。

很值得一提的是，前面第1节给出的D3和D9点群的18氮环的几何结构像极了下图所示的[18]annulene和[18]trannulene，它们都有长-短键交替的特征，关键差别在于每个C-H单元在18氮环里对应一个氮原子，也可以视为每个C-H共享的电子对对应于氮原子上的孤对电子。显然，18氮环中的氮原子的杂化状态可以视为与[18]annulene和[18]trannulene中的碳原子一样是sp2杂化。

由于18氮环中的孤对电子之间距离较近，因此孤对电子之间可能存在位阻-互斥效应，文中补充材料中S1节对此进行了专门的分析讨论并通过计算证实了这一点。D9结构下孤对电子之间的位阻-互斥效应尤为明显，这是其结构能量很高的另一个原因。

Mayer键级反映了原子间等效的共享电子对数。从前面的LOL等值面图上标注的Mayer键级可见，18氮环中N-N键近似可以视为具有单-双键交替特征。也由此，可以把18氮环的形成用下图来示意，即它是由氮气分子聚合而成的环状产物，聚合过程中每个氮气分子的N-N三键变成双键，并与每一侧的另一个氮气分子形成一个单键。这非常类似于乙炔形成聚乙炔的过程，只不过聚合产生的氮链远没那么稳定。值得一提的是，氮链类物质已经在高压合成的一些晶体中观测到了，只不过目前观测到的物质中氮链都是与过渡金属配位的状态。

为了能更充分地理解18氮环的成键，文中还使用Multiwfn按照《使用键级密度(BOD)和自然适应性轨道(NAdO)图形化研究化学键》（//www.umsyar.com/535）和《使用Multiwfn对周期性体系做键级分析和NAdO分析考察成键特征》（//www.umsyar.com/719）中介绍的NAdO方法对18氮环的N-N键进行了考察。每个化学键的模糊键级可以分解为一系列NAdO轨道的贡献，通过观看NAdO轨道的特征，就可以确切了解键级内在对应了什么样的相互作用。对一个较短和一个较长N-N键进行的NAdO分析的结果如下，只有贡献相对显著（大于0.2）的NAdO轨道被列了出来，括号外的是NAdO轨道对键级的贡献，括号内的是轨道的能量。可见较短的N-N键基本上可以视为一个sigma和一个pi键构成，因为相应的两个NAdO轨道对模糊键级的贡献0.926和0.858都很大且接近1，而其余的NAdO的贡献远远小于它们。对较长N-N键的模糊键级贡献接近1的仅仅有一个sigma特征的NAdO轨道。由此可以看到，18氮环中的N-N键其实在成键特征上并没有什么特别之处，可以近似视为是较典型的sigma+pi双键与sigma单键的交替出现。值得注意的是，如标注的NAdO轨道能量所体现的，较长N-N键的sigma轨道的能量是明显高于较短N-N键的，而且前者对键级的贡献更小，这说明较长N-N键的sigma轨道相对更弱。

文中还用Multiwfn对18氮环做了AIM拓扑分析，这是讨论化学键非常常用的方法，见《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）以及《AIM学习资料和重要文献合集（http://bbs.keinsci.com/thread-362-1-1.html）。下图左侧图中橙色小圆球是键临界点（BCP），对具有代表性的一个较长和较短的N-N键的BCP计算的结果如下图右侧所示。BCP位置的电子密度（ρ）及每电子能量密度（E/ρ）体现出较短的N-N键相对更强。电子密度拉的普拉斯函数值▽2ρ和能量密度都为负体现出两类N-N键都属于典型的共价键。BCP位置很大的ELF值体现出成键的电子具有很强的定域性，这也是共价键的典型特征。较短的N-N键的BCP处的电子密度的椭率ε明显大于0，体现出它具有显著的单套pi作用特征，而较长的N-N键的这个值则基本为0，说明成键区域电子密度几乎完全轴对称，故基本上是纯sigma的作用。这些结论和前述的其它分析相一致。

分子的静电势对于讨论分子间相互作用有极其重要的意义，相关信息参看《静电势与平均局部离子化能相关资料合集》（http://bbs.keinsci.com/thread-219-1-1.html）。基于《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图》（//www.umsyar.com/443）和《使用Multiwfn结合VMD分析和绘制分子表面静电势分布》（//www.umsyar.com/196）中介绍的方法，本文绘制了18氮环范德华表面的静电势填色图并做了表面静电势分布的定量统计，如下图所示。从左图可见，18氮环表面静电势分布很不均匀，既有很正的地方，最高达到30 kcal/mol，也有较负的地方，最负为-12.1 kcal/mol。这体现出取决于具体位置，18氮环既可以表现出明显的局部Lewis酸性特征，也可以表现出一定Lewis碱性特征（主要来自于孤对电子）。

上图右侧的不同静电势范围的面积统计图中同时包含18氮环和18碳环的情况。通过对比可明显看出18氮环的静电势分布范围较广，而18碳环仅分布在数值接近0的较窄的范围内，这在《全面探究18碳环独特的分子间相互作用与pi-pi堆积特征》（//www.umsyar.com/572）介绍的笔者的研究文章里也专门讨论过。这体现出18氮环有形成静电主导的分子间相互作用的能力，而18碳环则只能形成范德华作用主导的相互作用。笔者曾基于表面静电势提出过MPI指数衡量分子的等效极性，详见《谈谈如何衡量分子的极性》（//www.umsyar.com/518）。18氮环和18碳环的MPI分别为8.0和2.6 kcal/mol，对比可见18氮环的极性显著大于18碳环。值得一提的是靠偶极矩是无法如实区分它们的极性的，C1结构下的18氮环的偶极矩仅为0.161 Debye，和偶极矩精确为0的中心对称的18碳环几乎没什么区别。

本文还计算了18氮环的原子电荷以考察其电荷分布特征。原子电荷概念的介绍见《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！》（//www.umsyar.com/714）里提到的笔者的综述文章。无论是用ADCH原子电荷还是CHELPG原子电荷，得到的18氮环的原子电荷都分布在很窄的范围内（ADCH电荷在-0.023至0.029之间），因此18氮环里每个原子所带的净电荷甚微、感受到的化学环境高度相似。

18碳环具有明显的芳香性和整体的pi电子的离域性，这在笔者的论文Carbon, 165, 468-475 (2020)中层做过全面、深入的讨论。18氮环是否也有这样的特征？为了揭晓答案，文中首先使用Multiwfn计算了多中心键级，这是衡量电子多中心离域强度的非常流行的指标，见《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）和《使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键》（//www.umsyar.com/138）中的介绍。18氮环的多中心键级的计算结果近乎精确为0，也远远小于18碳环，因此从这一点上已经证明18氮环不具备像18碳环一样明显的整体电子离域特征，因而也没有芳香性。18氮环尽管有很多pi电子，但这些pi电子完全定域在一个个化学键上，并没有有效联通为整体。

《通过Multiwfn绘制等化学屏蔽表面(ICSS)研究芳香性》（//www.umsyar.com/216）介绍的通过绘制ICSS_ZZ等值面图考察芳香性的方法虽然昂贵，但是非常严格而且直观，也比流行的NICS更有说服力。为了进一步确认18氮环的芳香性，文中使用Multiwfn绘制了18氮环的ICSS_ZZ等值面图，如下所示，C1和D9构型的图都给出了，等值面数值也标出来了。图中绿色和蓝色区域分别是对垂直于环方向上施加的磁场产生屏蔽和去屏蔽的区域。无论哪种构型，其等值面特征都和具有典型芳香性分子的情况截然不同（即环中心区域完全是磁屏蔽，而环外侧是一圈连贯的去屏蔽），进一步证明了无论哪种构型，18氮环都不具备芳香性。

最后，文章还根据《使用AICD 2.0绘制磁感应电流图》（//www.umsyar.com/294）的做法绘制了18氮环两种构型的感生电流图。由下图可见，无论哪种构型都没形成像Carbon, 165, 468-475 (2020)中的18碳环那样环绕整体的磁感生电流，再次确认了18氮环的非芳香性特征。

7 总结

本文浅显易懂地介绍了近期发表的专门研究新颖的18氮环特征的文章ChemPhysChem, 25, e202400377 (2024)的主要内容，更多细节请读者阅读原文。此文不仅首次研究了18个氮原子形成的大环分子，也是首次系统性考察孤立状态下纯氮构成的长链状物质的特征，它可以视为是氮气分子作为单体产生的聚合物。相信本文可以明显拓宽大多数读者对纯氮物质的认识。本文的很多研究思想和分析方式也可以作为范例，在理论预测和分析其它新颖的化学物质时予以借鉴。同时本文也充分体现出使用Multiwfn程序做波函数分析考察新颖物质的电子结构的重要意义。

文/Sobereva@北京科音  2024-Aug-21

0 前言

笔者之前在《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）里介绍过大量衡量化学体系芳香性的方法，其中NICS十分流行，并且有大量变体。本文将介绍Inorg. Chem., 53, 3579 (2014)提出的一种和NICS密切关联的芳香性指数，称为free of in-plane component NICS (FiPC-NICS)，不仅比最原始的NICS(0)强得多，比认可度很高的NICS(1)ZZ在原理上也更好。笔者之前写过《使用Multiwfn绘制一维NICS曲线并通过积分衡量芳香性》（//www.umsyar.com/681），如今Multiwfn对此文介绍的功能进行了扩展，在后处理菜单中新加入了计算FiPC-NICS的功能。

下面就介绍FiPC-NICS方法的思想，并给出在Multiwfn中计算的完整实例。上面提到的博文如果还没看过的话一定要先仔细看一看，笔者假定读者已经读过了。Multiwfn可以在官网//www.umsyar.com/multiwfn免费下载。如果对Multiwfn不了解，建议阅读《Multiwfn入门tips》（//www.umsyar.com/167）和《Multiwfn FAQ》（//www.umsyar.com/452）。使用Multiwfn计算FiPC-NICS用于发表文章的话必须按照Multiwfn启动时的提示恰当引用Multiwfn的原文，给别人代算时也必须明确告知这一点。

1 原理

对于一个平面体系（至少是对于计算NICS的区域来说是局部平面的体系），任意一个点计算的NICS可以视为平行于环的in-plane分量（NICS_in）和垂直于环的out-of-plane分量（NICS_out）之和。假设环平面平行于XY平面，有以下关系
NICS = NICS_in + NICS_out = (NICS_XX + NICS_YY + NICS_ZZ)/3
NICS_in = (NICS_XX + NICS_YY)/3
NICS_out = NICS_ZZ/3
其中诸如NICS_ZZ等于相应位置的磁屏蔽张量的ZZ分量的负值。

FiPC-NICS原文认为NICS_in对环中定域化的电子（内核电子、sigma键的电子、孤对电子）敏感，而芳香性本质上来自于电子在环上的离域，那些定域化的电子对芳香性实质上是没有贡献的。因此，如果以NICS的方式来衡量芳香性，不应当纳入NICS_in部分。但也不是从环中心出发垂直于环方向上随便找个地方计算NICS_out就行，如果距离环中心太远，则NICS_out数值太小、对芳香性不敏感；而如果距离环中心太近，则定域化的电子会“污染”NICS_out，导致它不能充分反映芳香性。因此Inorg. Chem., 53, 3579 (2014)提出一种较严格地利用NICS思想衡量芳香性的指数FiPC-NICS，它定义为从环中心出发垂直于环的扫描路径上NICS_in为0的位置的NICS_out值，这可以尽可能避免定域化的电子产生的影响。

很常用的NICS(1)ZZ和FiPC-NICS有密切联系。设NICS_in为0的位置距离环中心垂直距离为x，则FiPC-NICS=[NICS(x)ZZ]/3。根据原文的测试，对于碳环的情况，x在1.2埃左右，和计算NICS(1)ZZ所用的x=1埃相近。这也体现出NICS(1)ZZ用于衡量芳香性的内在合理性，即统一用距离环中心1埃处来计算NICS_ZZ是个近似合适的位置。而对于非碳环类体系，考虑到不同元素原子大小的差异，用FiPC-NICS来对比芳香性会比用NICS(1)_ZZ更为严格。

如《使用Multiwfn绘制一维NICS曲线并通过积分衡量芳香性》（//www.umsyar.com/681）、《基于Gaussian的NMR=CSGT任务得到各个轨道对NICS贡献的方法》（//www.umsyar.com/670）和《将核独立化学位移(NICS)分解为sigma和pi轨道的贡献》（//www.umsyar.com/145）所述，计算NICS时还可以只考虑pi电子的贡献，pi电子对芳香性是有真正贡献的电子。因此至少对于碳环类体系来说，只考虑pi电子贡献的NICS(1)ZZ有可能比FiPC-NICS在衡量芳香性上还更理想，但FiPC-NICS原文没有进行对比。FiPC-NICS相对来说好处是在程序实现上更容易，不需要考虑 程序是否支持将pi和非pi贡献在NMR计算时进行分离，也不需要事先确定pi轨道序号。

FiPC-NICS原文还提出将NICS_in和NICS_out绘制成相关图的思想，如下所示，横坐标是NICS_in值，纵坐标是NICS_out值，图中每个点对应于从环中心出发垂直于环扫描路径上的一个点。扫描从环中心开始，对应下图最左端（r=0.0Å处）。图中各个体系的曲线都是逐渐收敛到(0,0)位置，因为距离环中心很远处NICS的各个分量都会收敛到0。图中黑线是最典型的芳香性分子苯，可见曲线是下凹的，绿线是最典型的反芳香性分子环丁二烯，可见曲线是上凸的，而其它都是非芳香性分子，曲线接近直线。由此可见，不同芳香性特征的分子的这种NICS_in vs. NICS_out相关性图的特征截然不同，这给判断芳香性提供了一个新的手段，值得借鉴。

2 实例

下面我们就通过Multiwfn结合Gaussian来计算一下FiPC-NICS。Gaussian用的是16 C.02版，Multiwfn是2024-Jul-13版。下面牵扯到的大多数操作和《使用Multiwfn绘制一维NICS曲线并通过积分衡量芳香性》是一致的，所以就不详细解释了。所有用到的文件都可以在//www.umsyar.com/attach/724/file.rar中得到。我们试图重现FiPC-NICS原文里苯的计算结果，本例使用和原文中相同的PBE0/6-311++G**计算级别。

启动Multiwfn，输入
opt.out  //这是Gaussian 16对苯在PBE0/6-311++G**下做几何优化得到的输出文件，在本文的文件包里。Multiwfn会读取最后一步的结构，分子平面平行于XY
25  //芳香性与离域性分析
13  //NICS-1D扫描、积分与FiPC-NICS计算
2    //通过一批原子定义环，将扫描的两个端点置于环中心上方和下方一定距离处，且连线垂直穿越环中心
1-6   //用于1-6号碳原子定义环
[回车]  //用环上的原子的几何中心作为环中心
0  //首端位于环中心
10  //另一个端点位于环上方10埃处
200  //扫描的点数。此例相当于每隔约0.05埃一个点，足够精细了。点数越多计算越耗时
1  //产生Gaussian输入文件
template_NMR_6-311++Gxx.gjf  //这是Gaussian做NMR计算的模板文件，在本文的文件包里。打开它便知计算是在PBE0/6-311++G**下进行的，和普通输入文件唯一的差异是坐标部分用[geometry]代替

当前目录下现在出现了NICS_1D.gjf。用Gaussian运行之，得到NICS_1D.out。然后在Multiwfn中接着输入
2  //读取NICS-1D扫描结果
NICS_1D.out  //在本文的文件包里
6  //计算FiPC-NICS

现在在屏幕上看到以下结果，说明FiPC-NICS为-9.58 ppm，是在苯环中央上方1.192埃处计算出来的。
FiPC-NICS is   -9.577714 ppm, at   1.192 Angstrom

FiPC-NICS原文中的结果为-9.59 ppm，位于环中心上方1.18埃，可见我们计算的与文献相符极好（文献用的是Gaussian 09）。

Multiwfn在当前目录下还输出了FiPC-NICS.txt文件，其中一部分内容如下。第1列是扫描的点的序号，第二列是距离环中心的距离，第3、4列分别是NICS_in和NICS_out。

...略
   20     0.955   -0.656100   -9.899500
   21     1.005   -0.500933   -9.902900
   22     1.055   -0.354133   -9.864733
   23     1.106   -0.215967   -9.788600
   24     1.156   -0.086633   -9.678467
   25     1.206    0.033833   -9.538367
   26     1.256    0.145400   -9.372367
   27     1.307    0.248267   -9.184367
...略

对FiPC-NICS.txt最后两列在Origin里绘制scatter+line图，如下所示，可见和上一节展示的FiPC-NICS原文里的图1中的苯的曲线完全吻合！

3 总结

本文介绍了FiPC-NICS芳香性指数，它比常用的NICS(1)ZZ更为严格，尤其是对于讨论无机的环状区域来说（例如FiPC-NICS原文里讨论的N和P交替构成的六元环）。从例子可见通过Multiwfn计算这个指数以及绘制NICS_in vs. NICS_out扫描图的过程相当容易，值得大家在实际芳香性研究中考虑使用。

注意用这个功能时扫描路径必须平行于X或Y或Z轴，如果环平面不平行于某个笛卡尔平面的话，需要利用《Multiwfn中非常实用的几何操作和坐标变换功能介绍》（//www.umsyar.com/610）中介绍的功能令要考察的环平行于某个笛卡尔平面，并且模板文件里应带上nosymm关键词避免Gaussian自动旋转朝向，参考《谈谈Gaussian中的对称性与nosymm关键词的使用》（//www.umsyar.com/297）。

使用Multiwfn计算FiPC-NICS用于发表文章的话请记得按照程序启动时的提示恰当引用Multiwfn的原文。

文/Sobereva@北京科音  2024-Jul-31

0 前言

衡量化学体系的芳香性的方法非常多，见《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）。研究芳香性的文章数目甚巨，但大多数都是对分子、原子团簇这样的孤立体系研究的，主要在于很少有程序能够支持将衡量芳香性的方法用于周期性体系。如今Multiwfn已经可以将很多方法用于周期性体系了，包括多中心键级、AV1245、AVmin、PDI、FLU、FLU-pi、PLR、HOMA、Bird、ELF二分值、Shannon芳香性指数、芳香环的环临界点属性，等等，还可以通过绘制LOL-pi函数图像直观考察共轭情况。在本文中，将对其中大部分方法在Multiwfn中的操作进行演示。作为例子的体系是一个共价有机框架化合物（COF）的一层，此体系在《使用Multiwfn对周期性体系做键级分析和NAdO分析考察成键特征》（//www.umsyar.com/719）中已经作为例子分析过，结构如下所示，本文要通过不同方法对比绿色虚线里三个环的芳香性的大小。可见1号环是C3N3，原子序号为1,2,21,22,11,12。2号环和3号环相当于萘片段中的两个六元环，2号环序号为19,20,50,49,18,48，3号环序号为14,48,18,44,46,16。这里序号是按照成键关系排的。

笔者假定读者已经阅读过《衡量芳香性的方法以及在Multiwfn中的计算》充分了解了各种衡量芳香性方法的特征，也假定读者有了Multiwfn的基本使用常识，不了解者建议阅读《Multiwfn入门tips》（//www.umsyar.com/167）和《Multiwfn FAQ》（//www.umsyar.com/452）。Multiwfn可以从官网//www.umsyar.com/multiwfn免费下载，注意务必使用2024-Jul-1及以后更新的版本（注意看Multiwfn启动时的更新日期的提示），否则情况和本文所述不符。

本文涉及的多数芳香性分析方法都是基于波函数的，上面COF这个体系的由CP2K程序在PBE/DZVP-MOLOPT-SR-GTH级别计算得到的molden格式的波函数文件和//www.umsyar.com/719里用的是一致的，请从此文的链接里下载，是其中的COF-MOS-1_0.molden。此文件的产生方法在此文里也详细说了，不会用CP2K者必须仔细看。本文涉及到的诸如HOMA这样的芳香性指数是基于几何结构的，只要给Multiwfn提供的文件里包含结构信息就可以，如cif、pdb、xyz、gjf、mol2等常见格式都可以。如果被分析的区域是跨晶胞的，则输入文件必须能给Multiwfn提供晶胞信息，哪些格式能提供晶胞信息见《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）里的说明。

2  计算多中心键级

多中心键级是非常好、很严格的考察芳香性的指标。这一节计算一下当前体系中三个六元环的六中心键级，以此考察它们六中心共轭的强弱，这正比于它们的芳香性。多中心键级的概念参看《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）和《使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键》（//www.umsyar.com/138）。

启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
25  //离域性与芳香性分析
1  //多中心键级
1,2,21,22,11,12  //计算1号环。原子序号按原子连接关系输入
此时得到多中心键级结果为0.0399。类似地，再分别输入19,20,50,49,18,48 [回车]和14,48,18,44,46,16 [回车]计算2、3号环，多中心键级分别为0.0212和0.0281。

根据六中心键级大小可见此体系中的C3N3环（1号环）的芳香性显著强于六元碳环，而不与C3N3环连接的六元碳环（3号环）的芳香性比与C3N3环连接的六元碳环（2号环）更强。

3 计算AV1245和AVmin

AV1245的概念和计算方法在《使用Multiwfn计算AV1245指数研究大环的芳香性》（//www.umsyar.com/519）中介绍过，这里就不多说了。这里用AV1245算一下前述的三个环的芳香性。

启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
25  //离域性与芳香性分析
2  //AV1245
1,2,21,22,11,12  //计算1号环。原子序号按连接关系输入
看到以下结果
AV1245 times 1000 for the selected atoms is    7.35835761
AVmin times 1000 for the selected atoms is     7.358354 (   11   12    2   21)
即AV1245乘以1000和AVmin乘以1000都为7.358。

再输入19,20,50,49,18,48计算2号环，结果为
AV1245 times 1000 for the selected atoms is    4.46380673
AVmin times 1000 for the selected atoms is     3.441219 (   20   50   18   48)

再输入14,48,18,44,46,16计算3号环，结果为
AV1245 times 1000 for the selected atoms is    6.30980750
AVmin times 1000 for the selected atoms is     4.014836 (   46   16   48   18)

可见，无论是从AV1245还是AVmin上来看，芳香性都是1号环>3号环>2号环，这和多中心键级的结论完全一致。

4 计算PDI

这一节用PDI衡量芳香性。注意PDI只能用于六元环。Multiwfn中PDI是基于模糊空间定义的离域化指数算的，对周期性体系默认用的是Hirshfeld原子空间划分。

启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
15  //模糊空间分析
5  //PDI
[回车]  //用默认的格点间距计算原子重叠矩阵（AOM）。对较大体系可以用更大的比如0.35 Bohr格点间距明显节约这一步的计算时间，对结果影响甚微
1,2,21,22,11,12  //第一个环里的原子序号
结果如下

Delocalization index of     1(C )   --   22(N ):    0.104520
Delocalization index of     2(N )   --   11(C ):    0.104520
Delocalization index of    21(C )   --   12(N ):    0.104520
PDI value is    0.104520

再依次输入19,20,50,49,18,48 [回车]和14,48,18,44,46,16 [回车]分别计算第2、3个环的PDI，结果分别为0.083658和0.096507。PDI越大芳香性越强，可见PDI给出的芳香性大小的结论也是1号环>3号环>2号环。

5 计算FLU和FLU-pi

先计算FLU芳香性指数。启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
15  //模糊空间分析
6  //FLU
[回车]  //用默认的格点间距计算原子重叠矩阵（AOM）
现在从屏幕上可以看到做FLU计算用的参数。输入第一个环里的原子序号1,2,21,22,11,12，得到如下结果

        Atom pair         Contribution          DI
   1(C )  --    2(N ):       0.000800        1.474298
   2(N )  --   21(C ):       0.000312        1.508760
  21(C )  --   22(N ):       0.000800        1.474298
  22(N )  --   11(C ):       0.000312        1.508761
  11(C )  --   12(N ):       0.000800        1.474298
  12(N )  --    1(C ):       0.000312        1.508761
FLU value is    0.003335

再依次输入19,20,50,49,18,48 [回车]和14,48,18,44,46,16 [回车]分别计算第2、3个环的FLU，结果分别为0.017882和0.012199。由于FLU越小芳香性越强，因此FLU的芳香性大小的结论是1号环>3号环>2号环，和前述芳香性指标结论相同。

下面再来计算FLU-pi。算这个需要告诉Multiwfn当前体系的所有pi占据轨道序号。可以按照《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）肉眼一个个看来记录序号，但太麻烦。建议按照《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）介绍的做法自动指认。启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
100  //其它功能（Part 1）
22  //自动检测pi轨道
0  //当前轨道是离域的（分子/晶体轨道属于这类）
马上Multiwfn就显示了占据的pi轨道序号，为47,50,60,61,64,70-73,77,81,82,84-89,92-94。下面计算FLU-pi，接着输入
0  //对识别出的pi轨道什么都不做
0  //返回主菜单
15  //模糊原子空间分析
7  //FLU-pi
[回车]  //用默认的格点间距计算原子重叠矩阵（AOM）
47,50,60,61,64,70-73,77,81,82,84-89,92-94  //pi占据轨道的序号
之后输入第一个环里的原子序号1,2,21,22,11,12，得到如下结果

Average of DI-pi is    0.387267
        Atom pair         Contribution          DI
   1(C )  --    2(N ):       0.000191        0.375470
   2(N )  --   21(C ):       0.000191        0.399063
  21(C )  --   22(N ):       0.000191        0.375471
  22(N )  --   11(C ):       0.000191        0.399063
  11(C )  --   12(N ):       0.000191        0.375470
  12(N )  --    1(C ):       0.000191        0.399063
FLU-pi value is    0.001148

再依次输入19,20,50,49,18,48 [回车]和14,48,18,44,46,16 [回车]分别计算第2、3个环的FLU-pi，结果分别为0.028167和0.026396。由于FLU-pi越小芳香性越强，因此FLU-pi的芳香性大小的结论是1号环>3号环>2号环，和前述芳香性指标结论相同。

6 计算HOMA

HOMA是流行的基于环上的键长特征衡量芳香性的方法。给Multiwfn提供的输入文件里有结构信息就行了。由于COF-MOS-1_0.molden里也包含结构信息，所以此例还是用这个作为输入文件。启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
25  //芳香性分析
6  //HOMA和Bird芳香性指数
0  //开始计算HOMA（基于默认参数）
现在屏幕上显示了计算HOMA用的参数。输入第一个环里的原子序号1,2,21,22,11,12，得到如下结果

        Atom pair         Contribution  Bond length(Angstrom)
   1(C )  --    2(N ):      -0.003592        1.349180
   2(N )  --   21(C ):      -0.001191        1.342742
  21(C )  --   22(N ):      -0.003592        1.349181
  22(N )  --   11(C ):      -0.001191        1.342741
  11(C )  --   12(N ):      -0.003592        1.349181
  12(N )  --    1(C ):      -0.001191        1.342741
HOMA value is    0.985651

可见HOMA为0.985651，并且环上各个键的键长以及对HOMA的贡献量都给出了。再依次输入19,20,50,49,18,48 [回车]和14,48,18,44,46,16 [回车]分别计算第2、3个环的HOMA，结果分别为0.598160和0.781579。由于HOMA越接近1芳香性越强，因此HOMA的芳香性大小的结论是1号环>3号环>2号环，和前述芳香性指标结论相同。

Bird芳香性指数的计算方法和HOMA非常类似，只不过是在主功能25里的子功能6里选择2而非此例的0而已，这里就不演示了。

7 计算Shannon芳香性指数

Shannon芳香性指数是基于被考察的环上的各个键的键临界点的电子密度定义的，因此计算它之前必须先用Multiwfn搜索临界点。这方面的操作在《使用Multiwfn结合CP2K做周期性体系的atom-in-molecules (AIM)拓扑分析》（//www.umsyar.com/717）里有详细说明，这里就不对细节做具体介绍了。

启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
2  //拓扑分析
3  //用每一对原子连线的中点作为初猜点搜索临界点
0  //在图形窗口里查看结果

在图形界面里要求只显示(3,-1)临界点，也就是键临界点，并且要求把临界点序号显示出来，看到下图。可见第一个环上的临界点序号为72,76,67,59,55,61。

点击图形界面右上角的Return按钮，然后输入
20  //计算Shannon芳香性
72,76,67,59,55,61  //第1个环上的BCP序号
马上看到结果：

Electron density at CP 72:   0.3352748688  Local entropy:   0.2993061920
Electron density at CP 76:   0.3318277479  Local entropy:   0.2979424307
Electron density at CP 67:   0.3352750629  Local entropy:   0.2993062683
Electron density at CP 59:   0.3318276569  Local entropy:   0.2979423945
Electron density at CP 55:   0.3352752224  Local entropy:   0.2993063310
Electron density at CP 61:   0.3318279246  Local entropy:   0.2979425011
Total electron density:   2.0013084835
Total Shannon entropy:   1.7917461175
Expected maximum Shannon entropy:   1.7917594692

Shannon aromaticity index:    0.0000133518

即Shannon芳香性指数为0.0000133518，每个键临界点的电子密度和局部熵也都给出了，这些是计算Shannon芳香性指数的中间量。

类似地也输入第2、3个环上的键临界点的序号，分别为62,64,54,45,43,51和64,71,82,84,75,65，结果分别为0.0013759152和0.0010309007。由于Shannon芳香性指数越小芳香性越强，因此芳香性大小的结论是1号环>3号环>2号环，和前述芳香性指标结论相同。

8 基于环临界点属性判断芳香性

在某个环中央区域的环临界点位置上，若垂直于环的方向上电子密度曲率越负，说明这个环的芳香性越强。这种判断芳香性的方法需要利用AIM拓扑分析得到的环临界点的属性。为了用这种分析，和上一节一样先用主功能2做拓扑分析，然后在选项0里只显示(3,+1)临界点，即环临界点，看到的图如下。可见第1、2、3号环的环临界点序号分别为66、53、73。

关闭图形窗口后输入
21  //计算顺着某个方向的电子密度梯度和曲率
66  //第1个环的环临界点序号
1  //指定某个方向
0,0,1  //当前体系平行于XY平面，因此计算Z方向的电子密度的梯度和曲率
结果如下

Electron density is                            0.0287379670 a.u.
Electron density gradient is                   0.0000000000 a.u.
Electron density curvature is                 -0.0261250099 a.u.

可见，在第一个环的环临界点位置上垂直于这个环的电子密度曲率为-0.0261250099 a.u.。类似地再使用这个功能考察第2、3个环分别对应的第53、73号环临界点垂直于环平面的电子密度曲率，结果分别为-0.0154770974和-0.0159388407 a.u.。由于数值越负芳香性越强，因此芳香性大小的结论是1号环>3号环>2号环，和前述芳香性指标结论相同。

9 ELF-pi二分点考察芳香性

这一节通过各个环上ELF-pi函数的等值面首次发生二分的ELF-pi数值（ELF-pi二分值）来衡量芳香性的大小。这个值既可以通过反复微调等值面数值来获得，也可以通过对ELF-pi做拓扑分析来获得，前者更直观，后者更精确，相关信息见《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）。这一节演示前者的做法。

首先需要产生ELF-pi的格点数据。启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
100  //其它功能（Part 1）
22  //检测pi轨道
0  //轨道都是离域形式
2  //设置其它轨道占据数为0
0  //返回主菜单
5  //计算格点数据
9  //ELF
9  //利用晶胞平移矢量定义格点信息
[回车]  //坐标原点为(0,0,0)
[回车]  //盒子三个方向尺寸和相应晶胞边长一致
0.15  //用较精细的0.15 Bohr格点间距，使得通过观看等值面获得二分值能尽可能准确
2  //将格点数据导出为cube文件

现在当前目录下就有了ELF.cub，对应ELF-pi的格点数据。用VMD显示其等值面（不会操作的话用《在VMD里将cube文件瞬间绘制成效果极佳的等值面图的方法》//www.umsyar.com/483里提供的脚本）。选Display - Orthographic用正交视角，在Graphics - Representation里将等值面数值（isovalue）由小到大调节，会发现在等值面数值为0.761的时候正好1号环上等值面首次发生二分，位置如下图箭头所示。因此1号环的二分点数值为0.761。

由上图可见，在1号环ELF-pi等值面二分的时候，2、3号环早就发生了二分、等值面间间隔已经很大了，这说明1号环的芳香性远强于2、3号环。

再将等值面数值调节到0.623，会看到2、3号环上首次出现了ELF-pi的二分，如下图红色箭头所示。这个键是被两个六元环共享的，怎么区分哪个芳香性更强？这需要再看其它的键。在下图紫色箭头所示的2号环上的位置，只要等值面数值再增加一点点就会二分，而3号环就没有这个情况，说明3号环的芳香性比2号环更强、环上的电子共轭遇到的瓶颈更少。因此ELF-pi的分析结论和其它方法完全一致。

10 LOL-pi图形分析

作为前述分析的扩展和补充，这一节绘制LOL-pi图直观展示一下单层COF体系的pi共轭情况，这能十分清楚地让大家理解三个环上电子离域特征的差异。这种分析在前述的《在Multiwfn中单独考察pi电子结构特征》文中有大量介绍，在笔者的Theor. Chem. Acc, 139, 25 (2020)中也有很多例子，欢迎阅读和引用。

启动Multiwfn并载入COF-MOS-1_0.molden后，依次输入
100  //其它功能（Part 1）
22  //检测pi轨道
0  //轨道都是离域形式
2  //设置其它轨道占据数为0
0  //返回
4  //绘制平面图
10  //LOL
1  //填色图
[回车]  //用默认的格点数
1  //XY平面
4.25  //当前体系中的原子都在Z=3.25位置，绘制平面上方1 Bohr处的图像显然应该设Z=4.25 Bohr

利用后处理菜单对作图效果的一番调节，得到下图。如果不会调的话，仔细把Multiwfn手册4.4节的所有例子都仔细领会了，结合后处理菜单中提示得很清楚的选项名字去理解，自然就明白了。当前用的色彩刻度是0到0.65，超过上限的区域显示为白色。

由上图可见，在1号环上pi共轭很显著，LOL-pi在环上分布得比较均匀。而如图上我标注的箭头所示，2号环上有三个键的pi共轭显著低于其它地方，因此这个环上的六中心pi共轭作用明显受到了很大遏制。而3号环上的pi共轭整体相对好点，但也不是很均匀，而是在不同键上有大有小，所以芳香性只是介于1和2号环之间。

11 总结

Multiwfn支持极为丰富的考察化学体系芳香性的方法，本文对其中已支持周期性体系的大部分方法的操作过程进行了简明扼要的演示。从结果可见尽管这些方法思想差异很大，但对于当前研究的这个体系中的三个六元环，它们给出的芳香性顺序完全一致，互相印证了彼此的可靠性。当然对于一些特殊情况，由于一些方法原理上的局限性、稳健性和普适性的不足也有可能导致结果存在差异。本文示例的只是一个很标准、简单的体系，希望大家能充分举一反三，将Multiwfn应用于广泛体系的芳香性的研究中。

使用Multiwfn做任何分析在发表文章时都请务必记得按照程序启动时的提示恰当引用Multiwfn原文，给别人代算时也必须明确告知对方这一点。