为什么[18]碳环不是完全的键长平均化？

作者 @Acid @杨知守

审稿人 @常思一二

Abstract

在问为什么C18没有键长平均化之前，我首先问了自己一个相反的问题：为什么我们期待它具有键长平均化？

正如题主所说，一个很自然的答案是，因为C18有4N+2个π电子（无论是垂直于平面的π键还是平面内的π键），符合Huckel芳香性的规则，因此结构上应当表现出芳香化合物键长平均化的特征。这是一个令人熟知的经验规律，推导也十分简单，也是为什么至今有机化学还在教Huckel方法的重要原因。

但这里的逻辑真的没问题吗？仔细回顾用Huckel方法推4N+2规则的过程就会发现，其实这一过程证明的是4N+2电子数是键长平均化的必要而非充分条件，即键长平均化必有4N+2电子，但4N+2电子未必有键长平均化（当然前者也并非是绝对成立的，只是在Huckel的框架内证明的结论）。

因此我首先在Huckel的框架内试图证明芳香体系表现出键长平均化这一现象，结果惊讶的发现在Huckel方法下苯环的键长平均化的结构居然不是最优结构，而是表现出了D6h->D3h的破坏对称性的形变。

另一方面，对称性破缺常常用姜泰勒效应来解释，简并基态的对称性破缺称为一阶姜泰勒效应，是为人熟知也非常好用的理论解释。而非简并基态的对称性破缺对应于二阶姜泰勒效应，虽然也常常被拿出来说事，看起来对称性破缺后HOMO确实能量降低了，但是把同样的计算放到苯环上就会发现，苯在D3h下HOMO能量也是降低的，因此二阶姜泰勒在这里只能算是对现象的描述而非解释。

那么究竟C18为什么会发生形变降低对称性呢？我与@杨知守进行了一系列的讨论和测试，选择了氢环为测试系统，认为类似的环状系统的结构是电子结构与核排斥等因素竞争的结果，单纯的轨道能量是倾向于对称性降低的，核排斥才是氢环中维持对称性的主要原因。这一结果大大出乎了我的意料，这里分享给大家，只为了提供一个新的视角，不代表认为芳香性的本质是核排斥等扩大化的结论。

Huckel方法推导4N+2规则

让我们首先回顾一下Huckel方法推导4N+2规则的过程：首先写出π体系的Hamiltonian，Huckel方法做出了若干假设，比如总能量等于轨道能量之和，忽略轨道之间的重叠，忽略除相邻轨道以外的轨道相互作用等，因此体系的Hamiltonian只包含对角元α代表原子轨道自身的能量，及β代表相邻轨道间的相互作用：

不难发现H具有很高的对称性，因此其本征函数可以直接猜出：

并由此解得H的本征值：

其中n是碳原子的个数，k代表态的编号(0~n-1)，j代表原子编号(0~n-1)。

显然，k=0能量取得最小值εmin=α+2β，k=n/2时（如果n为偶数）取得最大值εmin=α-2β（注意β<0），此外的任何一个k和n-k对应的态相同的能量。因此如果我们暂时不考虑三重态的存在，由于轨道的简并度，体系的电子数必须满足4N+2的形式，否则体系的前线轨道就会出现不全占满的情形，从而引发姜-泰勒(Jahn-Teller)效应。

可以说，4N+2的电子计数规则是Huckel方法作为一种半经验的计算方法最重要的推论之一。但是Huckel方法除了能得到4N+2规则以外，真的能推出键长平均化吗？

Huckel方法下单双键交替的苯环更稳定！

仔细分析一下上述推导的过程，注意到这个推导的前提是Hamiltonian具有前面给出的特定形式，而这一形式事实上已经蕴涵了键长平均化，否则不同相邻原子之间的β不全相等。因此这一推导的内在逻辑是：

键长平均化+4N电子数⇒JTE（不稳定）

根据反证法，如果要求体系稳定，那就不能有JTE，于是：

键长平均化+稳定⇒4N+2电子数

因此，4N+2电子计数规则是键长平均化的必要条件，但这一逻辑没有表明4N+2是否也是键长平均化的充分条件，也就是说即使知道体系有4N+2个电子，体系是否具有键长平均化的性质仅凭上述推导不得而知的。

因此接下来我们依旧在Huckel方法的框架内，尝试证明4N+2的充分性：已知体系包含4N+2个电子，能推出体系是否呈现键长平均化吗？

这一推导也不难，我们只需把相邻原子之间的作用稍作区分，单键连接的两个碳原子之间的作用记为β1，双键记为β2，即：

此时矩阵依然具有周期性的对称性，但周期变成了之前的两倍，不难解得

当β1=β2=β时，体系约化到之前的情形。如果我们从这一情形出发进行微扰，即有β1=β+δ，β2=β-δ，则

很容易发现，δ的存在是净增加能级的分裂的，能量降低的一支（成键轨道）随着δ的增加进一步降低，升高的一支（反键轨道）则进一步升高。ε(k)对δ画图如下：

可以看到，无论n和k是多少，δ的增加会使得成键轨道能量降低。因此按照Huckel理论的算法，无论是苯环还是C18，Huckel的预测都是单双键交替的。

同样的证明也可以用来说明聚乙炔的单双键交替结构。聚乙炔是聚合物，周期性边界条件下紧束缚模型(tight binding, TB)本质上跟Huckel方法是一样的，只不过k从离散的变成了连续的，因此这一证明可以清晰的体现出聚乙炔单双键交替结构的稳定性，在周期性体系中这一对称性破缺的现象称为Peierls distortion。但考虑到Huckel预测苯也是单双键交替的，因此这里用TB来说明聚乙炔的单双键交替并不具有说服力。这一计算更重要的结论是聚乙炔在单双键交替的结构下会打开一个能隙，金属变成了绝缘体。

二阶姜泰勒效应（second order Jahn-Teller effect, SOJTE）

文献中用来解释C18采取长短键交替的polyynic结构的理由大都是二阶姜泰勒效应，可以参见本问题的其它回答以及相关人员的专栏文章。[1][2]

由于JTE的定义有一定的争议（物理上有时特指vibronic coupling带来的symmetry breaking，而化学上通常指的还是电子效应带来的symmetry breaking），我们暂且抛开这个效应的适用性、与vibronic coupling的关系以及对计算水平提出的要求不谈，只单独分析用二阶JTE解释C18形变的逻辑是否自洽。

二阶JTE对分子从高对称性结构形变到低对称性结构的解释是，占据轨道（通常为HOMO）与非占据轨道（通常为LUMO）之间发生混合，形变时HOMO能量降低，LUMO能量升高，从而使分子变得更稳定。这个逻辑看似没有问题，但结合上面Huckel模型的计算就会发现，苯环形变时也应当表现出同样的HOMO降低而LUMO升高的趋势，而计算验证的确如此。

我在PBE0/def2TZVP下对苯分子进行柔性扫描，扫描的变量是C原子偏离平衡位置的角度，C原子到中心的距离以及H原子的位置放开优化，整体结构固定在D3h对称性，在DFT下能量显然是随着扫描的进行逐渐升高的。在PBE0扫描出的结构的基础上用不同算法做单点，得到的势能曲线如下图所示：

Gaussian中的Extended Huckel method (EHM)有三套不同的参数实现，我都跑了一遍，顺便测试了几种半经验方法（EHM也是半经验方法，但是与CNDO等几种方法还是有着本质的不同，这里的半经验方法指CNDO系的几种算法，下同）也跑了一下。本来还准备跑个CCSD不过耗时较长我懒得跑了，对于苯环我相信DFT还是可信的。从图中可见不同DFT的结果之间是很接近的，EHM不出所料随着D6h到D3h的形变能量的确是降低的。其它几种半经验方法略微高估了形变能，但定性还是正确的，这几种办法本质上是对HF方法中涉及到的积分做的近似，但包含原子核-原子核及电子-电子排斥，比起Huckel系方法简单的对轨道能量求和显然还是准确不少。

顺便提一句，Gaussian软件支持EHM的计算，平常根本想不起来用。这次我拿EHM对苯环做优化，发现无论初始结构如何，优化都是一步结束的。考虑到对DFT scan的结构做单点的结果是EHM倾向于D3h的结构，这大概就是为什么Gaussian的Huckel方法没有进行实质性的优化的原因吧。

而不同算法对形变后苯的HOMO/LUMO能量如下图，其中Y轴画出的是以D6h结构下的HOMO/LUMO能量为0/1进行平移和放缩后的轨道能量，<0/1表示该算法下HOMO/LUMO能量降低，反之则升高。由此可见，无论是Huckel、半经验还是DFT，形变后苯的LUMO都是升高的，而HOMO在所有DFT下都是降低的，只有CNDO/INDO两个半经验方法给出了升高的HOMO。

结合前文的Huckel模型下的计算，我得出的结论是苯在形变时HOMO能量降低、LUMO能量升高，但分子总能量不等于轨道能量，苯形变时毫无疑问总能量是升高的。

而把同样的计算在C18上重复一遍得到了更有意思的结果，先看总能量：

从图中可以得到许多有意思的结论：

1. EHM对于C18的结果与前往预测的相同，都是倾向于偏离高对称性的结构

2. DFT的结果开始出现差别，虽然整体形状依然有相似之处，但部分DFT的最优结构并非出现在x=0（高对称性结构），表现出了实验上STM给出的polyynic结构。

3. 半经验算法依然表现的较为夸张，无论是倾向于cumulenic的还是polyynic的结构，能量变化都比DFT给出的幅度大不少，不过CNDO/INDO/MNDO这三种方法得到的势能曲线与HF非常接近了，可见其对HF的近似还是比较准确的。

由于我试了十多种不同的泛函，上图中挤在一起看不清楚，我挑了一个结构出来把每一个泛函的能量标记了出来：

这个图就有意思了，撇开EHM和半经验方法不谈，DFT的表现是很有规律的，HF成分越高的泛函越倾向于polyynic结构，而纯泛函则清一色的倾向于键长平均化的cumulenic结构。这样的结果也在情理之中，纯泛函过于倾向电子离域而HF倾向电子定域是DFT理论界熟知的事实，range-separated functional则表现出的比较中庸，只略微倾向于polyynic结构。

轨道能量则与原始计算文献一致，所有算法都表现出了HOMO降低LUMO升高的趋势：

因此，我认为用二阶JTE解释C18的形变并不具有说服力，因为苯环随着形变的发生HOMO能量也是下降的。虽然我们常常把总能量的变化归因于某些轨道的能量变化，但倘若只看轨道能量，那其实就是Huckel模型采用的近似，而我们现在已经知道了Huckel模型不能解释苯环的键长平均化，那么其对于C18的polyynic结构的预测只能当作是一个巧合。

DFT的形式比较复杂，不便于理解，但HF也正确的预测了C18的polyynic结构，而DFT相比于HF则是削弱了polyynic结构的稳定性，说明电子相关应当是倾向于cumulenic的（当然做出这个结论如果用post-HF的计算结果来说明会更有说服力）。

但是既然HF能预测正确的结构而Huckel不能，于是我得出猜测，HF与EHM的区别应当就是维持结构高对称性的主要因素，这也是为什么我还测试了CNDO等几种半经验方法，因为CNDO等方法是对HF做出的近似，而他们的结论也十分相近，但模型上CNDO则更接近EHM，唯一增加的项就是核-核排斥以及（部分）电子-电子相互作用。这不禁让我怀疑，难道核排斥或者电子之间的排斥才是维持高对称性的因素，而平常讨论的轨道相互作用反而是破坏对称性的因素吗？

氢环模型

在与@杨知守的讨论中，最初我们的关注点仍然在电子结构上，Huckel忽略了次近邻相互作用、忽略了重叠积分、忽略了电子相互作用，他提出可以利用氢环代替π体系模拟这些因素对能量的影响，熟悉凝聚态物理的人应该能看出来有各种凝聚态模型的影子。在他的帮助下，我计算了Hn（n = 6, 10, 14, …）在不同键长和不同形变角度时的能量。但是在对比了考虑次近邻相互作用的结果跟完整的HF结果之后，我们并没有得到有意义的结论。而且EHM计算与最粗糙的Huckel模型不同，已经考虑了非近邻相互作用以及轨道重叠积分的贡献，但EHM计算得到的苯形变依然是有利的。只有核-核排斥和电子-电子排斥是EHM模型相比于其它方法完全缺失掉的部分，因此我利用氢环模型，尝试对总能量按照核-核排斥、核-电子吸引、电子-电子排斥进行了分解。

计算选取的氢环模型结构如下：记H原子个数为n，H原子间的平均键长为d，形变程度为α，则第k个H原子的坐标为：

这里其它参数的选取应当是显然的，偏移量

选为偏移角度乘以n是为了方便不同H原子个数之间进行比较。

对于不同的平均键长d、H原子个数n及偏移量α，在HF/STO-3G水平下计算Hn的能量如下图所示：

可以看出，键长较短（相互作用较强）时，体系倾向于键长平均化的结构；当键长很长时，体系倾向于长短键交替的结构（也就是以H2分子的形式存在）。最神奇的是对于中间键长（d = 0.7~0.8），H原子较少时倾向于键长平均化，H原子较多时则倾向于长短键交替，这与碳环中观测到的苯环键长平均化、C18长短键交替的结果一致。

为了确认一下这个模型的能量是可信的，对几个关键结构做了FCI/STO-3G和HF/cc-pVTZ的计算，总能量曲线与HF/STO-3G几乎完全一致。

然后根据上面的猜想对d=0.8的情形（因为此时不同n的曲线不同）按照核-电子吸引、电子-电子排斥、核-核排斥进行分解：

可以发现：

1. 核-电子吸引Ene是单调递减的，也就是说这一项倾向于长短键交替的结构，

2. 核-核排斥Enn是单调增的，表示核间排斥倾向于高对称性的键长平均化结构

3. 电子-电子排斥Eee很小，放大可见其也是单调增的

4. 不同n对应的曲线中Enn、Ene差不太多，Eee也都很小，加起来抵消后才出现了local minimum

电子能Ee=Ene+Eee主要由Ene控制，这与Huckel模型的预测一致，说明Huckel对电子能的估计定性是正确的，而总能量是Enn和Ee相互制约的结果，在这个体系中Enn是避免分子发生对称性破缺的主力，而Huckel模型和EHM计算中没有Enn的存在，因此给出了定性错误的结论。联系前往碳环中CNDO等半经验方法包含核-核排斥，其势能曲线与Huckel的结果有本质上的不同，也在一定程度上说明了核排斥项的重要性。

当然，这一分解只是能量分解的一种手段而非唯一手段，不同的分解方式可能给出不同的理解思路。尤其是这一分解方式并不能直接套到碳环上去，因为Enn和内层电子之间的排斥应当相互抵消，否则出现非常大的数减非常大的数得到一个很小的数的情形，并没有讨论价值（总不能说碳环的核电荷数高，所以有更强的对称化倾向吧）。如果把结论进一步扩大，说芳香性的根源是核排斥，显然也是没有意义的。

总结

芳香性的概念已经被玩坏了，为了描述分子的芳香性，各种各样的描述符被提出来，从4N+2的电子计数规则，到键长平均化等几何特征，再到化学反应性和现在文献常用的磁特征（NICS），这些描述符之间甚至只有弱相关性，到底什么是芳香性已经没人说得清楚了。不过我原以为至少对于碳环来说，芳香性还是well-defined，没想到简单的Huckel计算把我这一点仅存的幻想也摧毁了。芳香性绝不仅仅停留在4N+2电子规则，而是各种因素复杂的交织在了一起

二阶JTE也是一个常用的解释理论，可以用来解释H2O的弯曲构型、Sn2H4的非平面构型等，但是轨道能量和总能量只有相关性，难说有因果性，有些时候甚至可以是负相关的。这种时候把总能量的降低归咎于二阶JTE上，某种程度说何尝不也是一种cherry-picking

化学上的解释太过直觉和经验，一旦把背后的逻辑理清楚常常蕴含着错误的逻辑，但偏偏化学理论就建立在这些极度经验的规则上。而计算化学往往都是依赖高精度的计算，很少像物理系那样用toy model来说明问题（当然很大程度上也是因为物理系关心的体系常常无法负担从头算的解法）。因此本文试图用一个偏建模的方式去理解C18发生形变的原因，而不是通过“统计”的方式（也就是画一条总能量和轨道能量的曲线，把相关性当作因果性），希望可以提供一些新的idea

参考文献

[1] https://zhuanlan.zhihu.com/p/419864383

[2] https://zhuanlan.zhihu.com/p/165381617

[3] https://goldbook.iupac.org/terms/view/J03361

[4] https://en.wikipedia.org/wiki/Jahn%E2%80%93Teller_effect

[5] https://en.wikipedia.org/wiki/Pseudo_Jahn%E2%80%93Teller_effect

真 · 参考文献

[1] Kaiser, K.; Scriven, L. M.; Schulz, F.; Gawel, P.; Gross, L.; Anderson, H. L. An Sp-Hybridized Molecular Carbon Allotrope, Cyclo[18]Carbon. Science 2019, 365 (6459), 1299–1301. https://doi.org/10.1126/science.aay1914.

[2] Pereira, Z. S.; da Silva, E. Z. Spontaneous Symmetry Breaking in Cyclo[18]Carbon. J. Phys. Chem. A 2020, 124 (6), 1152–1157. https://doi.org/10.1021/acs.jpca.9b11822.

[3] Brémond, É.; Pérez-Jiménez, Á. J.; Adamo, C.; Sancho-García, J. C. Sp-Hybridized Carbon Allotrope Molecular Structures: An Ongoing Challenge for Density-Functional Approximations. J. Chem. Phys. 2019, 151 (21), 211104. https://doi.org/10.1063/1.5133639.

利益相关：发表过碳环相关的计算化学论文。

（1）

需要注意的是，Cn环和C6H6最大的区别在于前者的碳是sp杂化，而后者是sp2杂化，剩余的p电子数分别为12个和6个。如果你做过简单的休克尔轨道计算，就知道休克尔规则的4n+2是指6个p电子而非碳原子数，而sp杂化的碳环显然不符合这个规律。

因此虽然所有的论文都会提休克尔规则，但是我认为其实这个体系不适合用休克尔规则。

（休克尔规则拓展阅读：Luyao Zou：「知乎知识库」- 休克尔规则）

（2）

目前所发现的C2n环，可以参考这几篇论文（https://doi.org/10.1063/1.2838200 https://doi.org/10.1246/bcsj.20200345 https://doi.org/10.1016/j.cplett.2015.04.006），所有的碳环都是polyynic结构的。除了一篇有些可疑的论文以外（https://doi.org/10.1016/j.cplett.2018.05.062），我所接触到的偶数碳环都是单三键交替的polyynic而不是双键首尾相连的cumulenic结构。

比较简单粗暴的解释如下图：

其他研究者则做了一些有益的尝试。这篇论文（https://doi.org/10.1021/acs.jpca.0c02577）总结了一些常见的碳环的结构尝试回答为什么sp杂化的4n+2碳环并非cumulenic结构。预测在n>5的时候，将转变为cumulenic结构。我个人的看法是，因为此文中使用的方法的可靠性问题，我并不是很相信它的结论。但是n<5的情况因为研究较多，除了上文提到的那个例子以外（而且是很小的能量差），所有已知的C2n体系都是polyynic结构的。

（3）

至于说用jahn-teller效应或者说是vibronic coupling来解释polyynic而不是cumulenic碳环，我只能说这是个思路，但是目前我还没看到相关的使用可靠方法的计算工作，所以在这里我不做结论。

PS：为什么要强调方法可靠呢？因为用DFT做碳环结构优化时有可能会得出错误结论。可以参见这两篇论文中对DFT的benchmark，常用的DFT并不可靠，因此最好使用像CCSD这样可靠的研究方法（https://doi.org/10.1063/1.5133639 https://doi.org/10.1016/j.cplett.2015.04.006）。事实上以我的经验来说，不仅仅需要至少CCSD方法，还需要至少triple-zeta水平的基组才能得出有意义的结果。我的导师曾经考虑过研究用所谓的KDC model研究碳环的polyynic是否是vibronic coupling的结果，受限于计算方法和计算资源，最终作罢。这个问题目前还属于悬而未决。我认为如果想从理论上解决这个问题，至少需要理论化学/计算化学高年级博士生的水平。