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

0 前言

D-pi-A型分子普遍具有较大的第一超极化率，其中D（donor）是电子给体，A（acceptor）是电子受体，pi是指具有全局pi共轭特征的片段，可以称为linker或者bridge。2019年首次通过STM观测到的由18个碳原子构成的环状体系，18碳环（cyclo[18]carbon），同时具有平面内和平面外两套全局pi共轭特征，这点在《18碳环及衍生物的十分全面系统的研究综述已在Acc. Mater. Res.期刊发表！》（//www.umsyar.com/749）介绍的综述中有充分说明，如果你对18碳环缺乏了解十分推荐阅读此文。显然，18碳环的这种特征使得它也很有希望作为pi-linker构建具有显著第一超极化率的独特的D-pi-A型分子。基于这个idea，北京科音自然科学研究中心的卢天和江苏科技大学的刘泽玉等人对18碳环与不同donor、acceptor相连接构成的D-pi-A型分子的特征，尤其是光学性质，做了充分全面的理论研究。此工作近期已经发表在Phys. Chem. Chem. Phys.期刊上，欢迎阅读和引用：

Jingbo Xu, Jiaojiao Wang, Xiaohui Chen, Wenwen Zhao, Xiufen Yan, Zeyu Liu,* Tian Lu,* Aihua Yuan,* Design of donor–π–acceptor type cyclo[18]carbon derivatives for infrared nonlinear optical materials: a theoretical perspective, Phys. Chem. Chem. Phys., 27, 11993 (2025) https://doi.org/10.1039/d5cp00736d

下文将会对这篇研究文章的主要内容进行介绍并添加一些附加信息，使读者更好地理解文章的研究思想和主要结果，图片、表格来自于文章正文或补充材料。同作者之前发表过大量其它的和18碳环及其衍生物、复合物有关的研究工作，十分推荐阅览//www.umsyar.com/carbon_ring.html里的汇总，其中有很多篇研究工作都与18碳环及相关体系的光学性质有关。

1 18碳环衍生物的几何结构

此文研究的将18碳环作为pi-linker构建的示意图如下所示，共考虑四种情况：(1)两侧都连接氢 (2)一侧连氢一侧连硝基 (3)一侧连氨基一侧连氢 (4)一侧连氨基一侧连硝基。其中(4)对应最典型的D-pi-A形式的体系。而连氢的话，在电子激发时氢不会起到明显的给/吸电子效应，它的作用主要是让相连的碳从原本在碳环中的sp杂化状态变成sp2杂化。

此文用Carbon, 165, 468 (2020)论证过描述碳环十分可靠的ωB97XD/def2-TZVP级别对上述体系做了几何优化，得到的无虚频的笛卡尔坐标在文中的补充材料里都提供了。各体系的碳环部分的长、宽（上图的d1和d2）在下表给出了，可见无论碳环连的是什么基团，环的形状都没什么区别，都是椭圆形。

各个体系的结构图如下（补充材料图S2），键长以埃为单位标注了，并根据键长对键进行了着色。由图可见，这些18碳环衍生物的C-C键长度是交替变化的，这个特征和孤立状态的碳环一致。

2 18碳环衍生物的电子结构

前面的表格里QM和Qn分别是体系的各个部分的Mulliken和NPA电荷，原子电荷的相关知识见《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！》（//www.umsyar.com/714）介绍的笔者写的综述。可见无论哪种原子电荷计算方法，H-C18-H和NH2-C18-H中，18碳环都能获得电子，故带有明显的负电。而在H-C18-NO2和NH2-C18-NO2中，由于吸电子能力更强的-NO2从18碳环上吸了一些电子，导致18碳环部分的净电荷恰好接近0。

前面表格里的μ0是体系的基态的永久偶极矩。有趣的是H-C18-H的偶极矩并不为0，因为此体系并不是精确中心对称的（哪怕用很严收敛限优化也是如此），仔细看前面图中标注的键长可以看出这点来。H-C18-NO2和NH2-C18-H的偶极矩都较大，而D-pi-A特征最鲜明的NH2-C18-NO2具有特别大的偶极矩，高达15.5 Debye。

按照《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图》（//www.umsyar.com/443）绘制的4个18碳环衍生物的范德华表面静电势填色图如下所示，由此可以清楚直观地看出它们的电荷分布特征的差异。颜色越红、越蓝分别体现出其附近的原子带正电、负电越显著。并且表面静电势最大点和最小点的位置和数值也都标注了出来。可见虽然NH2-C18-NO2的碳环中部并不怎么带净电荷，但左右两端由于连接的基团的给/吸电子特征存在巨大差异，因而形成了显著的正负电荷分离现象，这也是为什么它有最大的偶极矩。

下图是基于《使用Multiwfn绘制态密度(DOS)图考察电子结构》（//www.umsyar.com/482）介绍的方法绘制的18碳环衍生物的TDOS和PDOS图，HOMO、LUMO位置以及相应的gap都标注了。在相同计算级别下，单独的18碳环的HOMO和LUMO能级分别为-8.45和-1.70 eV，对应HOMO-LUMO gap为6.75 eV，而当前研究的18碳环衍生物的gap在4.62-4.9 eV范围。可见接上基团，令两端的碳从sp杂化变成sp2杂化从而一定程度破坏平面内的芳香性会使得gap有所缩窄。PDOS曲线体现出HOMO、LUMO轨道几乎都是C18部分贡献的，但氨基，特别是硝基，对于其它前线轨道也有明显贡献，这暗示出一些低阶电子激发态会牵扯到它们。

3 18碳环衍生物的电子激发特征和电子光谱

本研究使用TDDFT在ωB97XD/def2-TZVP级别下做了电子激发计算，算了最低50个态，并模拟了电子光谱，如下所示。可见当前研究的18碳环衍生物的吸收分为三个部分：500多nm的可见光区的弱吸收、300多nm的近紫外区的弱吸收、200nm左右的更远紫外区的强吸收。如Carbon, 165, 461 (2020)的计算研究所示，独立的18碳环是没有可见光吸收的，而当前研究表面给它恰当接上H、氨基、硝基后能够令体系显色。基团的种类一定程度影响最大峰位置和峰高度，但影响程度很有限。

为了更好地理解这些体系的吸收光谱的内在本质，文中按照《使用Multiwfn绘制电荷转移光谱(CTS)直观分析电子光谱内在特征》（//www.umsyar.com/628）介绍的方法绘制了电荷转移光谱，如下所示。可见不管接上什么基团，>300 nm区域的吸收几乎完全都是18碳环片段上的局域激发，而氨基和硝基则对于200nm左右的吸收有少量贡献，使得对应的电子激发有一些电荷转移特征，但仍以18碳环区域的激发为主。

能量最低同时振子强度又显著的电子激发态被称为crucial state（关键态），对于当前研究的体系来说对应于可见光区的吸收。各个体系的关键态的信息如下所示，包括激发能、振子强度、激发态偶极矩相对于基态偶极矩的变化的大小

为了更好地了解这些激发的本质特征，文中按照《使用Multiwfn做空穴-电子分析全面考察电子激发特征》（//www.umsyar.com/434）介绍的方法对各个体系的关键态绘制了空穴（蓝色等值面）、电子（绿色等值面）分布图和热图，如下所示，图中还标注了两种衡量电子激发特征的定量指标，在434博文里也都详细介绍了。可见尽管这些激发主要都来自于碳环部分的pi电子，而结合定量指标D和Sr，还是能看出差别的。四个体系中，NH2-C18-NO2的关键态具有最显著的电荷转移特征，而且从蓝色和绿色等值面上可以看到是NH2向NO2方向整体一定程度转移电子。而H-C18-H的电荷转移特征几乎没有，空穴和电子几乎都是相对于体系中心对称分布的。H-C18-NO2的情况较接近H-C18-H，而NH2-C18-H的情况较接近NH2-C18-NO2。

电子激发态特征和第一超极化率之间有密切联系，通过SOS公式关联起来，在《使用Multiwfn基于完全态求和(SOS)方法计算极化率和超极化率》（//www.umsyar.com/232）里有详细介绍。SOS简化后可以得到双、三能级公式，见《使用Multiwfn对第一超极化率做双能级和三能级模型分析》（//www.umsyar.com/512）和《谈谈计算第一超极化率的双能级公式》（//www.umsyar.com/361）。双能级公式用来讨论影响第一超极化率的本质因素十分有用且流行。双能级公式指出静态第一超极化率与关键态的振子强度呈正比、与激发态相对于基态的偶极矩变化的大小呈正比、与激发能的三次方呈反比。根据前面表格里的关键态信息可见，相对于其它三个18碳环衍生物，NH2-C18-NO2的关键态具有最大的振子强度、最大的偶极矩变化，而激发能和它们则差不多，根据双能级公式可以预期NH2-C18-NO2有最大的第一超极化率。但双能级公式只适合用来解释和定性预测第一超极化率，更严格的结论、准确的数值还是需要依赖于用严格的方法定量计算，见下一节。

4 18碳环衍生物的非线型光学性质

此文使用Gaussian 16在ωB97XD/aug-cc-pVTZ(-f,-d)级别下以CPKS方法解析地计算了前述四种18碳环衍生物的二次谐波生成（SHG）形式的静态和动态的第一超极化率，β(-2ω;ω,ω)。aug-cc-pVTZ(-f,-d)基组是常用的aug-cc-pVTZ基组去掉非氢原子的f和氢原子的d极化函数，根据笔者以前的充分的测试，这种简化可以节约很多计算时间而并不会令计算精度下降太多。文中具体考察的量包括第一超极化率的总大小（βtot）、平行于偶极矩的分量（βvec）、各笛卡尔分量（βX、βY、βZ），这些量都可以通过Multiwfn基于Gaussian输出文件直接产生，见《使用Multiwfn分析Gaussian的极化率、超极化率的输出》（//www.umsyar.com/231）。

文中计算出的静态（∞ nm）和1460、1907 nm下的动态第一超极化率的总大小如下图所示。可见H-C18-H的β相当小，光是把一个H替换成硝基也没什么变化，但替换成氨基则能令β显著增大，而同时连上氨基和硝基构成典型D-pi-A体系后β变得颇大，这体现出氨基和硝基在提升电场响应性质方面的显著的协同作用。从下图还可以看到随着外场波长的减小，即外场频率的增大，对应的β有大幅提升，即这些体系具有很显著的频率-色散效应。

NH2-C18-NO2的静态的βtot为4462 a.u.，算是颇大了，同样是D-pi-A型分子的p-nitroaniline (NH2-C6H4-NO2)的实验气相值只有1400 a.u.。可见18碳环作为pi-linker有独特价值，是个颇好的pi-linker。

文中还使用《使用Multiwfn通过单位球面表示法图形化考察（超）极化率张量》（//www.umsyar.com/547）介绍的方法直观考察了各个体系的第一超极化率的各向异性，并基于《使用Multiwfn计算（超）极化率密度》（//www.umsyar.com/305）介绍的方法考察了不同区域、片段对第一超极化率的贡献，请感兴趣的读者阅读原文中的讨论。

为了考察溶剂效应对光学性质的影响，文中还在IEFPCM隐式溶剂模型描述的环己烷环境下也做了计算，发现环己烷会导致光谱略微红移、振子强度略微下降，但整体特征不变。环己烷还会造成绝大多数情况β的显著增加，但不同体系的β大小的趋势并未改变。

5 总结

本文介绍的Phys. Chem. Chem. Phys., 27, 11993 (2025)文章首次对18碳环与H、氨基、硝基构成的分子的几何结构、电子结构和光学性质做了充分的理论研究，研究证实18碳环可以作为很理想的pi-linker构造典型的D-pi-A体系以获得具有显著非线性光学性质的分子，文中构建的NH2-C18-NO2的第一超极化率已达到了很大的值。此文的研究工作显著拓宽了化学家们对碳环的衍生物的认识，对于未来利用碳环作为骨架构造非线性光学材料提供了重要的参考。值得一提的是，目前从小到大很多碳环皆已实验合成，以不同尺寸的碳环作为pi-linker，有望获得具有不同非线性光学性能的分子。

The way of plotting electrostatic potential for local region of molecular surface using Multiwfn in combination with VMD

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

1 前言

18碳环自从2019年在NaCl覆盖的Cu表面上被合成并观测到后，由于其独特的几何和电子结构，引发了理论化学家们的巨大关注并做了大量研究，研究也已逐渐扩展到其它尺寸的碳环以及碳环的衍生物、复合物。独特的18碳环与不同碱金属离子形成的结构是什么样，具有什么性质，是明显非常有意思、很值得探究的问题。

近期，北京科音自然科学研究中心的卢天和江苏科技大学的刘泽玉等人对18碳环与Li+、Na+、K+、Rb+、Cs+形成的复合物的各方面特征做了极为全面的研究，包括几何结构、结合强度、电荷分布、相互作用物理本质、光学吸收、非线性光学性质等。此研究对于将碳环类物质用于碱金属离子储存材料、构建光电材料等方面都有显著的指导意义。非常欢迎阅读此论文以及引用：

Yang Xiao, Xia Wang, Xiufen Yan, Zeyu Liu,* Mengdi Zhao,* Tian Lu*, Structure and Optical Properties of Alkali-Metal Ion (Li+, Na+, K+, Rb+, and Cs+) Endohedral Cyclo[18]carbon, ChemPhysChem, 26, e202500009 (2025) https://doi.org/10.1002/cphc.202500009

此文章可以通过此链接免费阅读：https://onlinelibrary.wiley.com/share/author/YZFCSBHKXAUFBD5NNJBU?target=10.1002/cphc.202500009

此文还被选为期刊的当期封面文章：

同作者之前对碳环类体系已经做过十分广泛、全面的理论研究并发表过大量研究论文，完整汇总见//www.umsyar.com/carbon_ring.html，里面还包含许多论文的评述、解读和附加讨论，欢迎查看。

下面，本文将对前述研究工作的内容的关键部分做简要介绍，使得读者快速、容易地理解文章主要内容，并对一些研究细节做附加说明。更具体的分析讨论请读者看完本文后阅读原文。此文研究的体系用M+@C18来表示，其中M+为Li+、Na+、K+、Rb+、Cs+。同作者还另有两篇论文与本文关系很大，《理论设计由18碳环与锂原子构成的电场可控的光学开关》（//www.umsyar.com/630）介绍的论文中考察了Li@C18，研究的侧重点是外场诱导Li原子在环内外的切换带来的效应，而《一篇文章深入揭示外电场对18碳环的超强调控作用》（//www.umsyar.com/570）介绍的论文中考察了Na+和Mg2+与18碳环的复合物，讨论了这两种离子对18碳环产生的外电场效应。

2 M+@C18的几何结构

此文通过Gaussian 16使用在ωB97XD/ma-TZVP级别下优化了各种M+@C18的无虚频结构（坐标在文章的补充材料里提供了），如下所示。俯视图和侧视图都给出了，同时考虑了M+在环内和环外两种结合构型。可见除了半径最大的Cs+结合在环内的情况时体系不是纯平面外，所有结构都是纯平面的。只有正电荷密度最高且半径最小的Li+、Na+结合在环内时令碳环的结构的变形肉眼可见，其它情况碳环依然保持像孤立状态一样的圆形。

文章使用ORCA程序在明显更高级的ωB97X-V/def2-QZVPP级别下对这些结构又计算了电子能量，并且结合上述级别的振动分析得到的自由能热校正量，计算了M+结合在碳环内和碳环外的自由能之差ΔG，并根据ΔG=-RTlnK计算了常温下的平衡常数K，如下所示。可见随着碱金属原子序数的增大，碱金属离子明显越来越倾向于结合在环内。这是因为周期越靠后的碱金属离子结合在环内时与碳环的色散作用越强，而结合在环外时这种效应虽然也有但没那么显著，毕竟结合在环内时碱金属离子能和所有碳原子较近接触，而结合在环外时只能接触一部分，且色散作用随作用距离衰减又很快。由于碱金属离子结合在碳环内的比率占绝对主导，因此后文不再考虑结合在环外的构型。

此文还计算了M+的原子电荷，确认了M+@C18中M几乎完全带一个正电荷，即18碳环上的电子向M+的转移/离域微乎其微，18碳环几乎完全处于电中性。而且还发现M+与碳环中碳原子的Wiberg键级数值甚微（尤其是半径大、极化作用弱的K+、Rb+、Cs+的情况），体现出M+与碳环之间的轨道相互作用甚微。如果你不了解键级的话，可参考《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）中的相应部分。

3 M+与18碳环的结合作用

利用Multiwfn程序，此文使用了我提出的非常流行的可视化展现片段间相互作用的IGMH方法以展现M+与18碳环的相互作用，IGMH方法定义的δg_inter=0.001 a.u.函数的等值面见下图，等值面的着色变化使用IGMH方法的标准色彩刻度。IGMH方法的介绍见《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）。由此图可以十分直观地看出不同碱金属离子与碳环的最主要作用区域。Li+和K+由于半径较小，因此在环内吸附时只与一部分碳原子有明显的相互作用，而K+、Rb+、Cs+与碳环的作用区域基本没区别，都是能同时与所有碳原子等同地作用。这些图的作用区域的等值面颜色都为绿色，体现出这些区域的电子密度都很低，这主要在于M+与碳原子之间共价作用甚微。

文章表S3给出了不同M+与18碳环的基于电子能量计算的相互作用能（ΔE_int）、结合能（ΔE_bind），以及标况下的结合自由能（ΔG_bind），如下所示，单位是kcal/mol。它们之间的关系介绍和相关讨论见《谈谈分子间结合能的构成以及分解分析思想》（//www.umsyar.com/733）。从相互作用能来看，M+与碳环结合强度是Li+ >> Na+ > K+≈Rb+≈Cs+。后三者与碳环作用强度相仿佛这一点和上面的IGMH图展现的情况一致。ΔE_bind和ΔE_int之间相差片段的变形能，可见只有Li+和Na+的情况造成18碳环在结合M+时有不可忽视的变形能，这和之前给出的结构图相对应，即只有Li+@C18和K+@C18中碳环偏离圆形是肉眼能明显察觉的。由于结合时的熵罚效应，标况下的ΔG_bind没有ΔE_bind那么负，但它们对于不同M+@C18的变化趋势是一致的。ΔG_bind都为明显负值，体现出标况气相下碱金属离子与碳环的结合在热力学上是充分自发的。

为了进一步剖析为什么不同的碱金属离子与18碳环结合作用存在差异，文中做了能量分解，结果如下所示。能量分解的相关常识见《使用sobEDA和sobEDAw方法做非常准确、快速、方便、普适的能量分解分析》（//www.umsyar.com/685）及里面提及的文章。可见Li+、Na+与碳环相互作用之所以显著强于更重的碱金属离子，关键在于它们的诱导项明显更大，即对碳环的极化作用明显更强。而K+、Rb+、Cs+由于原子半径明显更大、正电荷密度更小，因此极化能力弱得多，诱导项没那么负。还可以看到，周期越靠后的碱金属离子与碳环的色散作用越强，这在于同一族里周期表越靠后的元素极化率整体越大。然而，由于交换-互斥项也是随着碱金属离子周期越往后越大，和色散项产生了很大程度抵消，这导致K+、Rb+、Cs+和碳环的相互作用强度恰好都相仿佛。至于静电作用项，无论是哪种M+@C18，其起到的作用都微乎其微，这主要在于M+@C18体系中18碳环不仅净电荷基本为0，而且18碳环也不具备明显的极性。18碳环的极性问题在笔者的Carbon, 171, 514 (2021)论文里有专门的讨论。

4 M+@C18的光学吸收

此文对不同的M+@C18通过TDDFT在ωB97XD/ma-TZVP级别下计算了UV-Vis光谱，如下图所示。可见每个体系只有两个有明显光学活性的电子激发，碱金属离子的种类对于吸收光谱并没什么明显影响，都是仅在紫外区域有显著的吸收，这一点和孤立的18碳环很像，孤立的18碳环的电子光谱在Carbon, 165, 461 (2020)中笔者专门做了详细的研究。

文中还使用我提出的电荷转移光谱（CTS）方法考察了电子光谱的本质特征，在补充材料的图S1给出了。CTS方法的介绍见《使用Multiwfn绘制电荷转移光谱(CTS)直观分析电子光谱内在特征》（//www.umsyar.com/628）。由CTS图可以看出M+@C18的电子光谱几乎完全来自于碳环部分的局域激发，而碳环与M+之间的电荷转移激发只对光谱起到了微量的贡献。因此，碱金属离子的存在并不会对碳环的电子光谱的本质产生显著影响，只不过其带来的外势会多多少少对碳环的光谱造成一些定量层面的影响。

为了更清楚展现M+@C18有光学活性的电子激发的本质特征，此文按照《使用Multiwfn做空穴-电子分析全面考察电子激发特征》（//www.umsyar.com/434）介绍的方法绘制了此类体系的振子强度明显大于0的两个激发态的空穴和电子等值面，在下图里分别对应绿色和蓝色，这种分析远比看分子轨道图像来考察充分全面得多（当前这些激发每个都有4对MO跃迁有显著贡献，因此观看MO很不方便考察）。可见空穴和电子都分布在碳环部分，而且空穴和电子的等值面都同时环绕C-C键轴，因此电子激发同时具有平面内和平面外pi-pi激发特征。更仔细看的话还会发现这些激发的空穴主要分布在较短的C-C键上，而电子主要分布在较长C-C键上，即这些激发具有较短C-C键的pi电子向较长C-C键的非占据pi轨道跃迁的特征。

5 M+@C18的（超）极化率

极化率和超极化率是描述化学体系对外电场响应的关键的量，超极化率特征直接决定了体系用于非线性光学材料的可能性。此文用ωB97XD基于CPKS方法解析计算了M+@C18的极化率（α）和第一超极化率（β），并计算了半数值的第二超极化率（γ），所考察的具体指标都可以用《使用Multiwfn分析Gaussian的极化率、超极化率的输出》（//www.umsyar.com/231）介绍的方法基于Gaussian的polar任务的输出文件用Multiwfn获得。由于超极化率的准确计算对基组的弥散函数要求很高，为了得到尽可能准确结果，此文使用了很大的对每个角动量带两层弥散的d-aug-cc-pVTZ基组，对于此基组没定义的元素则使用d-aug-cc-pVTZ-PP赝势基组。

静态外场下计算的（超）极化率的结果如下图所示

由上图可见，碱金属离子的类型对体系的极化率和第二超极化率并没什么影响，但是第一超极化率的差异特别明显。Li+的βvec很显著，Na+略逊一些但也很大，这是因为如之前的结构图所示，只有这两个体系在环平面内偏离中心对称性明显。K+@C18和Rb+@C18的βvec都精确为0，这是因为它们都是精确的中心对称体系。Cs+@C18虽然βvec不为0但数值也很小，这是因为它仅仅在垂直于碳环的方向上轻微偏离中心对称性。

动态外场下的（超）极化率的计算结果如下图所示，可见外场波长越短、频率越高，对应的（超）极化率越大，即当前体系有明显的频率-色散效应。

为了对M+@C18的（超）极化率的本质特征了解得更充分，此文对Li+@C18还按照《使用Multiwfn极为方便地绘制(超)极化率密度和三维空间对(超)极化率的贡献》（//www.umsyar.com/683）和《使用Multiwfn计算（超）极化率密度》（//www.umsyar.com/305）介绍的方法绘制了（超）极化率密度图，如下所示。可见在所有图中Li上面几乎都没有等值面出现，直接体现了Li+对于（超）极化率没有任何贡献，这主要也是在于它没有价电子，只剩电子被束缚得很紧的内核部分，自然对外场的响应微乎其微。

文中还对Li+@C18绘制了单位球面表示法的图像，这对于考察（超）极化率的各向异性很有价值，绘制方法和分析方式见《使用Multiwfn通过单位球面表示法图形化考察（超）极化率张量》（//www.umsyar.com/547）里的介绍。如下图可见，对于极化率和第二超极化率来说，分子平面内的分量远大于垂直于分子平面的分量，这来自于环平面内碳环具有大范围全局pi共轭特征；另外，在分子平面内Li+@C18对外场的响应几乎是精确各向同性的。第一超极化率在垂直于环平面方向上精确为0，而在环平面内的分量则十分显著，完全平行于碳环被拉长的方向，这也正是体系偶极矩的方向。

6 其它

为了回应审稿人要求考察计算结果受泛函影响的要求，文中还给出了下图，全面对比了ωB97XD的结果和计算（超）极化率常用的BHandHLYP和CAM-B3LYP的计算结果。可见泛函的选择虽然会产生一些定量影响，但不同泛函算出来的不同的量之间的差异趋势是完全相同的，无疑当前文中的结论是可靠的。

极化率同时由电子和核振动共同贡献，通常只讨论前者，因为绝大多数情况下它占主导，而后者的贡献虽然通常较小但也值得了解。此文对不同的M+@C18也计算了振动极化率，如下所示，振动对极化率的贡献的比例确实相当小，仅有百分之几，极化率几乎都来自于电子分布对外场的响应。

7 总结

从以上信息可以看到，通过 计算与波函数分析手段，ChemPhysChem, 26, e202500009 (2025)这篇文章对18碳环与不同碱金属离子的复合物做了十分全面的理论研究，既拓展了化学家们对18碳环的认识，也对18碳环及衍生物在未来的潜在应用提供了很有价值的参考，文中的研究方式也值得在其它类似的研究中借鉴。如果你对碳环相关的非线性光学特征感兴趣，建议阅读《全面揭示各种尺寸的碳单环体系的独特的光学性质》（//www.umsyar.com/608）、《从18碳环的硼氮取代物中理论筛选出具有优异光学性质的分子：一篇CEJ期刊文章介绍》（//www.umsyar.com/742）、《深入揭示18碳环的重要衍生物C18-(CO)n的电子结构和光学特性》（//www.umsyar.com/640）了解更多信息。

自2019年18碳环（cyclo[18]carbon）首次在凝聚相中的观测被报道后，目前已经有约200篇和碳氮环有关的研究文章诞生。北京科音的卢天和江苏科技大学的刘泽玉等人到目前为止已经共同发表了20多篇和碳环有关的纯理论研究文章，内容充分涵盖了碳环的各个方面，包括化学键、电子离域、基态与激发态的芳香性、分子振动特征、环张力能、光学吸收，非线型光学性质、弱相互作用、外电场的影响、硼-氮掺杂效应、引入外部基团效应，等等，研究汇总见//www.umsyar.com/carbon_ring.html，此网页中还包含了笔者对很多研究论文写的深入浅出的介绍和重要的附加说明、讨论。

这些研究工作无疑非常值得进行系统性的总结。近期，卢天和刘泽玉在美国化学会的Accounts of Materials Research期刊上共同发表了名为Cyclo[18]carbon and Beyond-New Materials, New Properties, and New Opportunities的文章，访问地址：https://doi.org/10.1021/accountsmr.5c00131，目前可以免费阅览，欢迎阅读和引用！此综述对于希望一次性较为完整、系统地了解碳单环各方面主要特征的研究者来说极具价值，也非常适合作为踏入碳环研究领域的人第一篇阅读的文章。此文除了对作者的理论研究工作做系统总结外，对于碳环研究的历史、实验方面的背景信息也都做了介绍，并且还对未来碳环体系的进一步研究指出了方向。

由于此期刊对篇幅的严格限制，因此此文着重综述作者的包含18碳环在内的碳单环体系以及其衍生物方面的工作，而碳环与其它物质（如小分子、富勒烯、碱金属离子、超碱原子、8字形双环主体分子）形成的复合物的大量研究工作将在未来另文综述。

文/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、https://pan.baidu.com/s/17x4Pm0LINmt9iHDBLXWhhw?pwd=86rv。若未订阅，也可以免费阅读此文在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（主页：//www.umsyar.com/multiwfn），这是唯一可以实现此文介绍的所有方法的程序，有很多方法还是此程序独家支持的。Multiwfn对于大量相互作用可视化分析方法在计算学领域的普及和推广起到了关键性的作用，此综述里几乎所有图都是基于Multiwfn计算的数据绘制的。Multiwfn的相关信息见《Multiwfn FAQ》（//www.umsyar.com/452）和介绍文章J. Chem. Phys., 161, 082503 (2024) DOI: 10.1063/5.0216272。

此外，//www.umsyar.com/667介绍的卢天在2024年发表的Visualization Analysis of Weak Interactions in Chemical Systems（Comprehensive Computational Chemistry Vol.2中的一章）和前述的综述有极为密切的联系，此文介绍的方法专注于弱相互作用，涉及的方法更少，而由于有明显更多的篇幅，因而对它们的介绍明显更深入细致，故和上面的综述有极大的互补性，强烈建议阅读！

还十分推荐阅读以下文章及其中引用的文章获取更多相关知识以及计算操作细节：
使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用
//www.umsyar.com/621（http://bbs.keinsci.com/thread-28147-1-1.html）
使用IRI方法图形化考察化学体系中的化学键和弱相互作用
//www.umsyar.com/598（http://bbs.keinsci.com/thread-23457-1-1.html）
使用Multiwfn做Hirshfeld surface分析直观展现分子晶体和复合物中的相互作用
//www.umsyar.com/701（http://bbs.keinsci.com/thread-43178-1-1.html）
谈谈范德华势以及在Multiwfn中的计算、分析和绘制
//www.umsyar.com/551（http://bbs.keinsci.com/thread-17271-1-1.html）
在Multiwfn中单独考察pi电子结构特征
//www.umsyar.com/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研究分子动力学中的弱相互作用
//www.umsyar.com/186
使用Multiwfn做aNCI分析图形化考察动态过程中的蛋白-配体间的相互作用
//www.umsyar.com/591（http://bbs.keinsci.com/thread-21826-1-1.html）
通过电子定域化函数(ELF)、价层电子密度分析讨论亲核进攻的方向
//www.umsyar.com/606（http://bbs.keinsci.com/thread-23982-1-1.html）
通过价层电子密度分析展现分子电子结构
http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB201709252

此外，笔者的18碳环及相关体系的研究文章中充分利用了文中介绍的许多方法，是很好的应用实例，研究汇总见//www.umsyar.com/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中的计算》（//www.umsyar.com/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计算公式都涉及到对占据轨道和空轨道的循环，并且空轨道与占据轨道能量之差出现在分母项中，这一点看《深入理解分子轨道对磁感生电流的贡献》（//www.umsyar.com/703）里的公式1、2、3便知。总的来说，HOMO-LUMO gap越大，空轨道与占据轨道整体能量差越大，导致NICS和感生电流越接近0。因此，如果两个体系的gap相差较大，不应直接凭NICS或BCS判断二者的芳香性强弱，还应该把gap差异对结果的影响考虑进去。在《深度揭示互为等电子体的苯、无机苯和carborazine的芳香性的显著差异》（//www.umsyar.com/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程序绘制磁感生电流图和计算键电流强度》（//www.umsyar.com/702）介绍的SYSMOIC算出来的环丁二烯的BCS为-0.7941（长C-C键）和-0.7499 a.u.（短C-C键），而苯的C-C键的BCS大小则小得多，只有0.4216 a.u.。

如《正确地认识分子的能隙(gap)、HOMO和LUMO》（//www.umsyar.com/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)研究芳香性》（//www.umsyar.com/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独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（//www.umsyar.com/696）介绍的文章中用了NICS(0)ZZ讨论。

实际上对于苯分子也有类似情况。下图是苯分子的pi以外的电子产生的感生电流图，可见由于局域的sigma电子产生的局部dia电流的存在，在环中心也有等效的para电流（红色箭头所示）并使得环中心产生磁去屏蔽效应。不过由于苯的pi电子造成的全局的dia电流在环中心产生的磁屏蔽效应远强于para电流的磁去屏蔽效应，因此苯的NICS(0)ZZ还是负值。

我极其推荐阅读《使用Multiwfn绘制一维NICS曲线并通过积分衡量芳香性》（//www.umsyar.com/681）中的讨论，了解sigma和pi电子对NICS_ZZ的贡献在垂直于环方向上的变化特征。

《18个氮原子组成的环状分子长什么样？一篇文章全面揭示18氮环的特征！》（//www.umsyar.com/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电子结构特征》（//www.umsyar.com/432）的做法绘制它的LOL-pi图，也会看到各个氮都有定域的孤对电子，根本没有像苯的pi电子一样在环上充分离域。然而，下图所示的[N6H6]2+平面构型的分子上方1埃处的感生电流（磁场垂直于屏幕向外）则是完全顺时针的dia电流，和具有显著芳香性的苯的情况看起来完全一样！

为什么[N6H6]2+的平面构型不是极小点？这是因为如下图的轨道占据情况所示，尽管此体系表面上有10个pi电子，看似满足Huckel的4n+2芳香性规则因而芳香性会使此构型被稳定化，但e2u对称性的轨道的反pi特征显著（对于苯来说是没有占据的），因此pi共轭作用很弱，再加上如此众多的孤对电子间的强烈互斥作用，显然令其极不稳定。

为什么[N6H6]2+的平面构型像普通芳香性体系一样有显著的贯穿整个环的dia电流？这需要利用《深入理解分子轨道对磁感生电流的贡献》（//www.umsyar.com/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碳环及相关体系的理论研究工作汇总见//www.umsyar.com/carbon_ring.html。完全由氮、硼原子交替构成的环B9N9，以及18碳环部分由硼、氮取代的结构，也已被研究过，例如《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（//www.umsyar.com/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在//www.umsyar.com/carbon_ring.html列举的各种碳环及衍生物的研究中都被证实可以很合理描述18碳环及其硼、氮掺杂体系的电子和几何结构。

2 (BN)nC(18-2n)的电子结构

下图展示出了不同的硼氮取代的18碳环的HOMO、LUMO能级和HOMO-LUMO gap，通过Multiwfn结合VMD按照《用VMD绘制艺术级轨道等值面图的方法）（//www.umsyar.com/449）方法绘制的等值面图也给出了。可见随着BN单元数n的增加，LUMO先轻微减小然后再明显增大，而HOMO先轻微增加再明显减小，这使得HOMO-LUMO gap在n=3时达到了最小值。

如《使用Multiwfn基于完全态求和(SOS)方法计算极化率和超极化率》（//www.umsyar.com/232）里的公式所示，激发能出现在计算极化率的SOS公式的分母部分，而HOMO-LUMO gap越小很大程度暗示激发能整体越小，因此由此可以预期HOMO-LUMO gap最小的(BN)3C12或相邻的(BN)2C14、(BN)4C10具有最大的极化率，这与后文的实际计算结果一致。

当n不很接近9时，从上图可以看到HOMO、LUMO几乎都出现在碳的部分，这和下图（原文中图S1）所示的基于《使用Multiwfn绘制态密度(DOS)图考察电子结构》（//www.umsyar.com/482）介绍的方法绘制的PDOS图相吻合，碳构成的片段的PDOS总是主导最高占据能级和最低非占据能级。

以上信息直接暗示出电子激发最容易的区域、反应活性最高的区域、对（超）极化率贡献最大的部分，都是(BN)nC(18-2n)中的共轭碳链部分（至少对于n不很接近9时）。

从18碳环开始，随着取代的BN单元数逐渐增加，HOMO-LUMO gap先轻微下降然后逐渐显著上升的本质原因何在？文中认为原因在于以下几点：
• 有显著芳香性的体系的gap倾向于比芳香性更弱的相似物的更大，例如《深度揭示互为等电子体的苯、无机苯和carborazine的芳香性的显著差异》（//www.umsyar.com/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》（//www.umsyar.com/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的极化率、超极化率的输出》（//www.umsyar.com/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极为方便地绘制(超)极化率密度和三维空间对(超)极化率的贡献》（//www.umsyar.com/683）的方法绘制了极化率密度，其中下面这张图展现了(BN)4C10的垂直于环平面的极化率分量的极化率密度的±0.12 a.u.等值面，等值面包裹的区域是对αzz有主要贡献的区域。可见主要贡献的部分是氮的孤对原子，以及碳链部分的pi电子很丰富的较短的C-C键。注意BN取代的18碳环和原本的18碳环一样都具有C-C键键长交替的特征。

(BN)nC(18-2n)的第一超极化率的总值（βtot）、在偶极矩方向的投影的大小|βvec|，以及《使用Multiwfn计算与超瑞利散射(HRS)实验相关的量》（//www.umsyar.com/499）介绍的超瑞丽散射（βHRS）实验测量的βHRS在下图给出了。由于18碳环和B9N9是中心对称的，二者的β皆精确为0。从18碳环开始，随着BN单元数的增加，各种β都是先增大后减小。这体现出BN单元数可以对体系的非线型光学（NLO）性质起到充分的调控作用，这个发现给设计具有特定性能的NLO材料提供了一种新思路。其中，(BN)4C10具有最大的βtot和βHRS，并且数值相当大，这是本文理论筛选出的具有关键价值的分子！

文章还使用《使用Multiwfn通过单位球面表示法图形化考察（超）极化率张量》（//www.umsyar.com/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)直观分析电子光谱内在特征》（//www.umsyar.com/628）介绍的CTS方法将总光谱特征进行分解的图，如下图右侧所示。可见在300多nm部分的峰主要来自于碳链区域内的电子重排的贡献，而200nm左右的吸收的成份更为复杂，同时牵扯到碳链和BN单元部分。

上面这个300多nm的峰对应于激发波长为323.8 nm的S0-S7激发，由于它没有主导性的单一分子轨道跃迁，因此文中按照《使用Multiwfn做空穴-电子分析全面考察电子激发特征》（//www.umsyar.com/434）介绍的极为普适的方法对其绘制了空穴和电子分布图，如下所示，紫色和黄色分别对应于空穴和电子分布。可见此激发的空穴和电子确实主要分布在碳链部分，仅有极少部分分布在氮和硼上，因此可以基本指认为碳链部分的局域激发。值得一提的是，按照《谈谈计算第一超极化率的双能级公式》（//www.umsyar.com/361）和《使用Multiwfn对第一超极化率做双能级和三能级模型分析》（//www.umsyar.com/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如下，体现了此文的关键性结果

笔者此前还对碳单环体系及其衍生物的各方面特征做过十分广泛的理论研究，发表的大量相关论文和深入浅出的介绍博文汇总见//www.umsyar.com/carbon_ring.html（不断更新）。

下面将对这篇Chem. Eur. J.文章的主要内容进行浅显易懂的介绍，便于读者更容易理解文章的研究结果，同时额外附上不少分析和计算细节以帮助读者能够将文中的研究手段举一反三运用到自己的研究中。原文里还有很多细节信息和讨论，故请阅读下文后阅读原文。此文不仅研究C16芳香性本身，文章在分析思想和方法学方面对于其它分子的芳香性研究也很有借鉴意义。

1 一些研究和计算细节

介绍结果之前，有一些研究和计算细节值得说明。

C16的基态和激发态势能面的极小点结构是不一样的，因此激发态极小点结构下的激发态芳香性与基态极小点结构下激发态的芳香性也是不同的，而且可以相差很大。此文专注于研究C16基态极小点结构下的基态、三重态和五重态芳香性，这避免了引入激发态势能面上结构的弛豫造成的芳香性改变。而若把结构弛豫效应考虑进去会导致分析讨论变得复杂，也无法“干净”地体现电子激发的改变本身对电子结构的直接影响。激发态势能面上的基态结构对应的那个位置被称为Franck-Condon（FC）点，因此此文研究的C16的激发态芳香性被称为Franck-Condon激发态芳香性。

此文对C16做的计算主要使用Gaussian 16程序在ωB97XD/def2-TZVP级别下做，这在//www.umsyar.com/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.344和1.221 Å，前者的键长交替（bond length alternation, BLA）值0.150 Å明显大于后者0.123 Å。芳香性倾向于令共轭环上的键长平均化，因此从几何结构上虽然无法判断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中的计算》（//www.umsyar.com/176）、《使用Multiwfn计算FiPC-NICS芳香性指数》（//www.umsyar.com/724）的相关部分，其中通常来说最可靠的是NICS(1)ZZ。NICS(1)ZZ的计算结果如以下表格左半边所示，比0负得越多说明芳香性越强、比0正得越多说明反芳香性越强。可见用GIAO和CSGT方式得到的总值差不多，而CSGT还可以分解成不同轨道的贡献便于更好地了解本质，见《基于Gaussian的NMR=CSGT任务得到各个轨道对NICS贡献的方法》（//www.umsyar.com/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计算时使用其它程序产生的轨道作为初猜》（//www.umsyar.com/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曲线并通过积分衡量芳香性》（//www.umsyar.com/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绘制磁感应电流图》（//www.umsyar.com/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)研究芳香性》（//www.umsyar.com/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方法极高精度计算激发态》//www.umsyar.com/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，见//www.umsyar.com/380）等当CASSCF初猜。

下面就通过例子简单演示一下具体实现方法，本文中涉及的文件都可以在//www.umsyar.com/attach/740/file.zip下载。Multiwfn可以在//www.umsyar.com/multiwfn免费下载，相关常识看《Multiwfn FAQ》（//www.umsyar.com/452）。如果读者还不知道Dalton怎么安装、运行、产生输入文件的话，看《 程序Dalton的编译方法和运行方式简介》（//www.umsyar.com/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的初猜波函数的方法》（//www.umsyar.com/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中的标识》（//www.umsyar.com/51）。虽然Dalton也可以用笛卡尔型基函数，但Multiwfn往Dalton输入文件里写入初猜轨道系数的功能不支持此情况。
• 只要一个 程序能输出Multiwfn可以识别的带有基函数信息的波函数文件，并且用的是球谐型基函数，则波函数都可以通过以上方式作为Dalton的初猜。支持的文件详见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/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计算双自由基以及双自由基特征的计算》（//www.umsyar.com/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的话看《谈谈片段组合波函数与自旋极化单重态》（//www.umsyar.com/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产生自然轨道的方法与激发态波函数、自旋自然轨道分析实例》（//www.umsyar.com/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观看分子轨道》（//www.umsyar.com/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支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）里专门说了。对于传轨道到ORCA的情况我之前还专门写过《将Gaussian等程序收敛的波函数作为ORCA的初猜波函数的方法》（//www.umsyar.com/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》（//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等能量分解方法分析，色散项应当占所有吸引项总和的大部分，轨道相互作用项应当占比微小