衡量芳香性的方法以及在Multiwfn中的计算.doc

衡量芳香性的方法以及在Multiwfn中的计算补充：在此文撰文之后Multiwfn又支持了ICSS方法图形化分析芳香性，见此贴通过Multiwfn绘制等化学屏蔽表面研究芳香性（/sobereva/item/9e672259b67d4d9408be1768）衡量芳香性的方法以及在Multiwfn中的计算文/SoberevaFirst release: 2013-Jan-23 Last update: 2014-Nov-20图文摘要本文首先先介绍芳香性的基本概念，然后对比较重要、目前比较流行的和近年来新提出的衡量芳香性的方法依次进行介绍，并谈谈笔者的看法。其中可以在Multiwfn程序中实现的方法会介绍操作过程。最后对实际应用中芳香性指标的选择进行讨论。1什么是芳香性芳香性是一个十分古老，重要，又含糊的化学概念。苯是最具典型的芳香性分子，也是芳香性的原型分子。与芳香性有关的文章近十几年来增速十分迅猛，如今可以说已经过热、被过分炒作。新的衡量芳香性的指标也不断被提出，同样在近年来增速迅猛，目前总计已有数十种。这些指标体现了芳香性的不同侧面，其中绝大部分依赖于量子化学计算。芳香性难以给出一个确切、唯一的定义，实际上芳香性这个词包含的内容被不断地广义化，以至于越来越不可能给出一个既简单、精确又能被所有研究者接受的定义。芳香性分子在能量、结构、反应性、磁性质、电子结构等方面都展现出独特的特征。对于大量芳香性体系，利用统计分析，可以发现它们总是同时具有很多特征，如：键长均衡化、电子呈高度整体离域性、外磁场下能形成整体诱导环电流、有较大离域化能、结构稳定等等，因此这样来看只要考察某分子的某个方面特征就能衡量其芳香性高低。但是对于不少分子，它们具有的各种特征之间的相关性并不满足上述大规律，这被称为芳香性的多维性，因此为了合理地描述芳香性就不得不同时考虑数种基于分子不同性质的芳香性指标。下面的表格列举了芳香性、非芳香性和反芳香性体系一部分常见特征，有的内容并不很准确，仅供参考，取自Chem.Rev.,105,3716(2005)。虽然芳香性这个话题越来越复杂化，不同研究者的观点不同。一个经典、古老、被化学工作者所普遍认知的关于芳香性的说法是Huckel规则：环上的分子轨道若具有4n+2电子占据则环具有芳香性。但是这个说法比较狭义，从理论化学角度也没触及物理本质。如果要对芳香性的本质做一下解释，那么寡人的观点是：对于一个环，若电子能够在整个环上高度离域，则环具有芳香性。虽然如上所述，芳香性具有多维特征，但很多都是外在表象，真正的源头是电子的离域。若电子没法在环上充分离域，那么芳香性就无从谈起。这里说的环是指电子的整体离域路径构成明显的环路，环的几何结构可以偏离圆形，可以是扭曲的不在一个平面的。环形离域路径不一定非得沿着分子的拓扑结构，比如某个局部结构是-C1-C2-C3-，若C1和C3离得近且轨道匹配，电子可以直接从C1跨过C2向C3离域；甚至还可以在分子片段间，乃至不同分子间有很强离域而构成整体环形离域路径，这被称为跨键或跨空间芳香性。另外还有很多文章研究球形芳香性，但这不属于本文涉及范围，相关综述见Chem.Rev.,105,3613(2005)。通常能够明显整体离域的是pi电子（绝大多数有机分子中sigma键定域性太强），这样的芳香性就是pi芳香性。如果整体靠的是sigma电子，就叫sigma芳香性。如果又有pi又有sigma的贡献，叫二重芳香性。虽然研究的最多的是有机分子的芳香性，但小的金属团簇（或金属和非金属混合团簇）的芳香性，即全金属（半金属）芳香性也受到广泛关注，对于这些体系sigma键的电子和delta键（d轨道构成）的电子和金属芳香性密切相关，有的团簇有人指出同时具备sigma,pi,delta三重芳香性，有的有pi,delta二重芳香性，等等，种类繁多。对于过渡金属团簇还可以利用f电子产生芳香性。本文并非芳香性的综述，不会去回溯芳香性的历史，不会去列举五花八门稀奇古怪的芳香性体系，不会去讨论一些新奇的概念。本文的目的在于归纳、简要介绍有历史意义的、目前流行的和新提出的衡量芳香性的方法，并对之发表一些寡人的观点。此文并非一览无余地囊括所有已经提出的衡量芳香性的方法，因为数目太多，而且有很多实际意义不大，或者纯属跟风。对于那些能在Multiwfn程序中计算的方法本文会介绍操作过程。用Multiwfn目前的版本已经可以算大量芳香性指数，包括：AdNDP、FLU/FLU-pi、PDI、ATI、PLR、HOMA、Bird指数、多中心键级、Shannon芳香性、LOL/ELF-sigma/pi、RCP处的电子密度及垂直于环平面的密度曲率、EL指数。本文所用Multiwfn为2.6.1版，可在此处免费下载：http:/M。以后版本可能操作步骤略有不同，以相应版本手册或屏幕上的提示为准。不要用老版本，本文许多和Multiwfn相关的内容是从2.6.1版才开始支持的。由于基于Huckel规则判断芳香性的方法属于经典、有重大意义的方法，所以先在第二节对之进行介绍。之后，会按照基于磁性质、基于几何结构、基于能量、基于电子离域性、基于电子分布的顺序来介绍各种衡量芳香性的方法。最后，将对本文涉及的方法进行比较和总结，谈谈各种场合下笔者推荐的方法。2基于Huckel规则的方法Huckel规则源自他在1931年的研究，但是直到1951年，Huckel规则在Doering的研究文章JACS,73,876中才被明确化为众所周知的4n+2规则。Huckel规则原本是基于MO轨道来分析，但是碰到含有多元环情况通常就瞎了，此时需要基于AdNDP轨道才能分析局部的环的芳香性。2.1MO结合休克尔规则判断芳香性这种判断芳香性的方法就是考察分子轨道(MO)图形，找出环上的MO，如果这些轨道上总电子数满足4n+2，则这个环具有芳香性。如果满足的是4n，则具有反芳香性。比如苯，以下是其双占据的三个pi轨道图，它们一起构成六中心的共轭环路，共6个电子，满足4n+2，故具有典型的芳香性。再看个复杂点的例子，中性状态的Li7团簇。图取自文献。HOMO是构成pi芳香性的轨道。如果是两个电子占据，即n=0，则将是典型的pi芳香性。但目前它只有一个电子占据（准确来讲应当叫SOMO），因此比典型pi芳香性要弱，这可以算作是部分pi芳香性。HOMO-1、HOMO-1和HOMO-2是构成sigma芳香性的轨道（轨道相位类似于苯的pi轨道），它们都被双电子占据，n=1，故是完美的sigma芳香性。由于此体系sigma和pi都对芳香性有贡献，所以算双芳香性。对于一个体系，若多种类型轨道都对芳香性有贡献，则整体的芳香性会因累加而变得更强。如果有的轨道呈芳香性而有的呈反芳香性，则它们会彼此抵消，整体是芳香性还是反芳香性就看是谁起主导作用了。不同电子态对应不同电子排布，因此也影响芳香性。例如某闭壳层体系单重态是4n反芳香性，如果把原先参与反芳香性的轨道上的一个电子电离掉，那么4n就被破坏了，从而反芳香性被大大削弱。如果这个电子被激发到相同环上的空轨道上，构成三重态，那么当前电子结构可以视为是本应构成4n+2的一批双占据轨道中有目前两个是单占据的（或者说，如果把这两个单占据轨道都再补上一个电子就满足了4n+2），故是不完美的4n+2，此时会有部分芳香性，即电子激发实现了反芳香性和芳香性的转化。正如Li7例子所示的，并非只有平面结构才能根据考察轨道来判断芳香性。但是，用于结构复杂，特别是非平面的情况，这种根据MO分析芳香性的方式往往很不好用、很含糊、很麻烦。因为MO图形往往较为复杂，有很多节点，涉及到体系诸多原子，时常没法判断哪些轨道是参与芳香性的；而且很多MO涵盖的空间范围广（尤其是大体系），它们只是部分地参与构成某些环的芳香性，计算4n/4n+2时是否把它们算进去根本说不清楚，或者说这么判断根本在原理上就行不通。而且，这种看轨道图形数电子数的方法的判断只是定性的，没法定量地判断芳香性的强弱。如今经常有文章利用这种方法讨论芳香性，讨论的结果令人匪夷所思。一方面是这种分析方法局限性太大，另一方面也是作者自身有问题，轨道选得莫名其妙，令人不禁觉得：凭什么把这些轨道算进4n/4n+2里？。我不很鼓励使用这种方法分析芳香性，或者至少也得结合使用一种可靠的能定量计算芳香性的方法。Huckel规则是对一般体系来说的。对于Mobius型分子，即共轭环的链扭转了(2M+1)*piM=0,1,2.弧度的分子，满足的是Mobius芳香性。相当于把Huckel规则反过来，即4n对应芳香性，而4n+2对应反芳香性。2.2AdNDP轨道结合休克尔规则判断芳香性上一节介绍的方法的主要局限性之一是MO离域性太强，即往往同时涵盖很多原子。当体系原子较多，有很多环可能同时具备芳香性/反芳香性，就很难靠MO来判断。例如菲就是这样，没法直接靠高度离域的MO来讨论它的三个环的芳香性。这时需要使用AdNDP方法，原文见PCCP,10,5207(2008)。这种方法获得的AdNDP轨道是“半定域化轨道”。AdNDP轨道既不像MO离域程度那么高，也不像传统的定域化轨道定域程度那么强（但中心或双中心轨道），而是介于两者之间，既在一定程度上定域化，又能表现体系局部区域的离域性。关于AdNDP的原理和在Multiwfn中的操作另有一帖专门介绍，就不在这里详谈了，见使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键（/sobereva/item/02b67c3ca3aa46c7382ffac6）菲的在边缘芳香环上的pi型AdNDP轨道图形如下：可见这三个轨道和苯的三个pi轨道很相似，因此它们构成了六中心共轭。这三个轨道占据数并非像苯一样都是完美的2.0，而是在1.82.0之间，基本符合4n+2规则。这说明这个环虽然有明显的芳香性，但是和苯相比要弱一些。假设基本符合的是4n规则，那么就是反芳香性。如菲的例子所见，AdNDP方法很有用，不仅能考察局部的芳香性或反芳香性，还能根据占据数大致说明其程度的强弱。但是AdNDP的威力被某些文章作者（包括一些大牛）过分鼓吹、胡作非为乱发文章。AdNDP的最大软肋就是搜索AdNDP轨道的过程对一些复杂体系任意性往往太强，人为干预太大，不同研究者可以得出不同结果，以至于，在JOC上某篇文章，coronene中心的根本算不上是芳香性（按照NICS还应该算反芳香性）的环居然被作者胡乱分析成芳香性。反正AdNDP的结果怎么解释都有道理，别人也不好质疑，但是结果可靠性有多高（也即AdNDP搜索过程含糊性有多大），只有作者自己心里有数。寡人绝非不鼓励用AdNDP，只是说，对于某些体系，如果用的过程中自己觉得心里很没底，干脆就别用了，勿在文章中强词夺理。目前乱用AdNDP分析芳香性的文章实在多得很，对于复杂的体系的分析结论，大家阅读时不要不加怀疑地相信。如果某个体系的芳香性比较重要，笔者鼓励大家自行重新找一遍AdNDP轨道以检验文章结论。3基于磁性质的方法本节介绍的方法基于体系对外磁场的响应行为。3.1核独立化学位移（NICS）及相关的方法可以说NICS是目前用得最多、接受度最高的衡量芳香性的方法。原文是JACS,118,6317(1996)，在Chem.Rev.,105,3842(2005)中有详细的综述，在Org.Lett.,8,863(2006)对几种定义进行了对比分析。NICS原始的定义是：环上重原子的几何中心处的各向同性化学屏蔽值的负值，以ppm为单位。因为考察的这个点不是原子核位置，所以叫“核独立”。对于芳香性体系，一般这个屏蔽值是正值，即NICS为负，本质是因为外磁场导致的共轭环上的环电流产生的感生磁场会一定程度对外磁场有抵消（屏蔽）作用。对于反芳香性体系，环电流的感生磁场和外磁场方向相同，故会加强外磁场，所以这个屏蔽值为负，NICS即为正。实际上，在实验上也有类似于NICS的做法来测定芳香性，但用的是比较轻的原子核，诸如He_3、Li+离子、桥氢原子来作为探针，让它们与芳香环构成复合物稳定在芳香环的面上，通过测定它们的化学位移就可以衡量环的芳香性。虽说这用的是实实在在的原子核来测定，但和NICS的思想本质是完全一致的。多数成熟的量子化学程序都能很容易地计算NICS，这是NICS之所以这么流行的主要原因之一。这里以Gaussian计算菲的边缘的环和中间的环为例来说明。获得几何中心的X、Y、Z坐标就是把环上重原子的X、Y、Z坐标分别取平均。假设边缘和中央的环的几何中心为(0.000000,-2.155707,-0.338053)和(0.000000,0.000000,0.859657)，就在Gaussian输入文件原子坐标后面写上：Bq0.000000-2.155707-0.338053Bq0.0000000.0000000.859657其中Bq代表虚原子，这使得Gaussian输出这些位置的NMR信息。然后RouteSection部分写上NMR关键词，以及所用理论方法、基组，然后开始运算即可。默认使用的是GIAO方法计算NMR，绝大部分文献报道的NICS值都是这个方法计算的。如果用其它方法，如NMR=CSGT，当基组不很完备时结果会与GIAO的稍有不同。建议大家的研究文章中要注明计算NMR所用的方法，以便读者重现。计算NICS用不着很高级别的方法和基组就能得到合理结果，通常用B3LYP/6-31+G*就可以。从输出文件中可找到这两个环中心的NMR信息25BqIsotropic=8.7039Anisotropy=5.0411XX=12.0647YX=0.0000ZX=0.0000XY=0.0000YY=6.3843ZY=-0.7977XZ=0.0000YZ=-0.8627ZZ=7.6628Eigenvalues:5.97588.071312.064726BqIsotropic=5.4377Anisotropy=3.7470XX=1.3628YX=0.0000ZX=0.0000XY=0.0000YY=7.0146ZY=0.0000XZ=0.0000YZ=0.0000ZZ=7.9357Eigenvalues:1.36287.01467.9357其中3*3矩阵是磁屏蔽张量。之所以描述为张量的形式是因为电子对不同方向射来的磁场屏蔽强度是不同的。这个张量矩阵乘以外磁场矢量得到的矢量的负矢量就是感生磁场矢量。磁屏蔽张的对角元的平均就是Isotropic后面的值，也就是各项同性屏蔽值，这也是NMR实验所对应的结果。边缘环的NICS即是-8.7039，中央环的NICS即是-5.4377，由于都是负值表明都呈现芳香性；由于前者更负，所以表明边缘的环的芳香性更强。最初的NICS定义，也被称为NICS(0)，被指出在不少体系都不很适用。主要毛病是它取的是各项同性值，而非垂直于环平面的屏蔽张量的分量。然而只有这个分量值，才和环上离域的电子对应的环电流所产生的感生磁场直接相关，才真正能清楚地展现芳香性特征，这个问题在Org.Lett.,8,863(2006)有详细讨论。这篇文章推荐用NICS(1)_ZZ来衡量芳香性，笔者也确实发现这个指标远比NICS(0)好得多得多。到现在仍有很多文献没完没了地批判NICS，实际上批的都是NICS(0)，大部分NICS(0)不适用的体系（诸如(HF)3环状三聚体被NICS(0)误判为芳香性）在NICS(1)_ZZ上都没问题。奇怪的是，这些文献总是无视NICS(1)_ZZ的存在却紧咬着NICS(0)不放，简直是莫名其妙。NICS(1)_ZZ就是指的在环平面上方1埃处（具体来说，是以环中心为起点向垂直于环平面方向挪1埃），垂直于环平面方向的屏蔽张量分量值的负值。假设分子平面是XY、YZ、XZ平面，那么这个分量值就是分别指的ZZ=、XX=、YY=后面的值。由于在环上方1埃处sigma轨道效应影响很小，所以NICS(1)_ZZ衡量的主要是pi芳香性。如果想在NICS框架内把sigma和pi对芳香性的贡献单独分析，可参见将核独立化学位移(NICS)分解为sigma和pi轨道的贡献（/sobereva/item/9d7a96ab8f6adea828ce9dc6）一文的讨论。对于被研究的环不处在笛卡尔平面上，或者环是扭曲的情况，计算NICS(1)_ZZ是很麻烦的，不过Multiwfn有专门的功能来解决这个问题，见利用Multiwfn计算倾斜、扭曲环的NICS_ZZ（/forum.php?mod=viewthread&tid=297）。计算NICS用的环中心怎么定义有含糊性，虽然原文用的是几何中心，但这明显不是理论上最合适的。使用质心，或者AIM理论的电子密度的环临界点(RCP)都是可以的。有文章（如doi:10.1002/wcms.1115）表明使用RCP比使用几何中心在结果上更合理。RCP可以通过Multiwfn的拓扑分析功能很方便地得到，见使用Multiwfn做拓扑分析以及计算孤对电子角度（/sobereva/item/c84fa5089156806dd55a11c3）。对于计算几何中心和质心，在Multiwfn也提供了一个小功能来实现，这比手算要方便。在启动Multiwfn并载入含有分子结构信息的文件后（pdb/xyz/fch/wfn等文件类型皆可），选100，再选21，然后输入环上的重原子序号，如3,4,8,9,10,7，则它们的几何中心和质心就会输出出来。在NICS的基础上，有人还做了进一步扩展，比如不考虑某个特定的点的屏蔽值，而是对垂直于环平面的整个面的屏蔽值进行扫描来分析，这等于将核屏蔽值作为一个实空间函数来考虑；再比如NICS-rate方法（CPL,493,376(2010)），它以NRR(NICS-RatesRatio)来衡量芳香性，NRR=|rate(Max)/rate(Min)|，这里rate(Max)和rate(Min)分别代表穿过环中心垂直于环平面的直线上NICS导数的最大值和最小值。不过这些方法计算起来麻烦，比起NICS(1)_ZZ没什么显著优点，所以一般没必要考虑。有人认为NICS不适合衡量尺寸不同的环的芳香性，但也有证据表明NICS对环的尺寸依赖性虽然有但是不太大。3.2磁感应电流密度与AICD相关问题和操作细节笔者已经在使用AICD程序研究电子离域性和磁感应电流密度（/sobereva/item/07d17c23d959de8e9d63d1c6）和使用GIMIC计算和分析磁感应电流密度（/sobereva/item/2fcd6479efa95d346cc37ccc）中有过详细探讨，这里就不细说了，只是十分粗略地涉及一下。磁感应环电流在上一节已经提到了。当环上电子离域性很强，能在环上自由移动，类似于导体，外加磁场时就会像环形线圈一样产感应电流，并且会在一定程度上抵消外磁场而表现抗磁性。感应电流方向如下所示电流越强芳香性越强。如果感应电流的方向和左手定则相反，则会加强外磁场，表现出顺磁性，说明环具有反芳香性，而且电流越强反芳香性越强。如果没有产生明显的贯穿整个环的感应电流，说明这个环是非芳香性的。感应电流不同位置的强弱和方向可以通过作图来直观考察，也可以对某些键的截面进行积分来获得电流大小的定量数据，这用GIMIC程序可以做到。AICD方法则是根据屏蔽张量矩阵元定义了一种实空间函数，它的数值越大的地方对应电流密度越大处。由于它是标量函数，所以可以很方便地通过绘制等值面图来考察，并可以认为数值较大的等值面包含的区域内电子有较强离域性，其功能某种程度上类似于后面要提到的ELF。而且，这个实空间函数的二分点（详见后文的ELF）大小表明了两个区域间电子相互离域的程度。如果某处的AICD等值面只有将数值调得很大的时候才分离为两个等值面，那么说明电子容易在这两个等值所含区域间离域，环电流也容易通过，这对于证实是否有跨空间芳香性很有用，也是一种判断芳香性强弱的定量指标。实现AICD方法有同名的程序，这程序可以将电流密度矢量直接标注AICD等值面图上，考察感应电流分布十分方便，例如：绘制感应电流密度图是十分有用的方法，能让电子离域路径、芳香性/反芳香性/非芳香性区域很客观地展现出来，但是对于较大的多环芳烃这样有多个共轭环连续出现的体系分析起来容易造成误导，内部电流密度的相互抵消导致一些环的芳香性特征在感应电流密度图上被掩盖掉了。如下图所示的coronene，环电流在方向相反的地方抵消掉了（右图只是简化的示意，其实这么表示不严格），因此导致表面看起来仿佛6个外围的六元环不是芳香环，而内部的环表现出反芳香性的假象。当芳香环数目更多，误导性往往更强。当外磁场方向与环垂直时，NICS_ZZ和磁感应电流密度在物理本质上有很直接的对应关系。如果看到某个环的环电流与左手定则相同/相反，那么表现在NICS_ZZ数值上也会是芳香性/反芳香性。这里NICS_ZZ未必是NICS(1)_ZZ，例如环电流是sigma电子贡献的，那么考察NICS(0)_ZZ就比考察NICS(1)_ZZ更妥当，因为平面上方1埃处sigma电子的效应就比较弱了。由于感应电流密度图分析方法用于分析较大的多环芳烃的芳香性上不很适合，这也导致NICS_ZZ不十分适合用于这类体系。GIMIC相对于AICD的优点是：(1)不仅支持HF/DFT，还支持后HF。特别是开壳层体系相关作用对电流密度影响强，尤为适合用CCSD。(2)电流密度的分析能达到定量级别，比如对截面进行积分相对缺点是(1)可视化不方便(2)没法将电流密度分离为正则轨道的贡献，因此不适合分析多重芳香性，特别是金属芳香性体系的芳香性本质。而AICD则允许只考察指定的分子轨道产生的电流密度。论初学者入门门槛，其实AICD和GIMIC都稍高。GIMIC程序就不支持最常用的Gaussian，得用turbomole、CFOUR/ACESII这样的程序；而使用AICD还得自己有Gaussian的源码并且能够成功编译。除了前面已提及的AICD和GIMIC，还有一种CTOCD-DZ(continuoustransformationoforiginofcurrentdensity-diamagneticzero)方法也被不少研究者用于绘制感应电流密度，见JPCA,105,9553(2001)。此方法支持HF/DFT，并可以分解为轨道的贡献。然而只有修改版的SYSMO程序能实现此方法，而此程序又不公开，所以这里就不多提它了。MRE是基于图论和Huckel理论的计算感应电流的方法，这在后文有详细介绍。3.3Aromaticringchemicalshieldings(ARCS)按道理来说，先得已知感应环电流，才能得到感生磁场，才得到了磁屏蔽值。ARCS（PCCP,1,3429(1999)）将这个过程反过来，先用量化程序算出穿过环中心且垂直于环的直线上各点的屏蔽值，以此为目标数据，通过拟合来得到感应环电流，并根据环电流大小和方向来判断芳香性。ARCS的拟合模型用的是电磁学上最常见的无限细的圆环导体，其半径也是被拟合的参数之一，通常比较接近于环的几何半径。实际上ARCS方法和按照前面所讨论的先计算电流密度再对某个键的截面进行积分来得到环电流数值在物理本质上并无差异，只是计算过程不同罢了。ARCS方法计算环电流有很大局限性，目前只能用于单个环，而且环必须接近圆形才行。如果将拟合的模型改为椭圆，倒也有可能用在诸如萘等链状分子上得到分子整体的环电流。目前ARCS是对每个环用一个环形导线模型描述，但如果用多个，或许描述得更为细致，比如对于sigma+pi双芳香性体系有可能同时拟合出sigma环电流和pi环电流。原则上只要自己用的量化程序能算NICS就能靠ARCS方法得到环电流，某种程度上比分析感应电流密度方便点，不过拟合ARCS环电流的程序得自己写，笔者认为其可靠性也不如直接分析感应电流密度。3.4磁化率提高(Magneticsusceptibilityexaltation)和抗磁化率各向异性磁化率是顺磁化率（正贡献）和抗磁化率（负贡献）之和，它们都用二阶张量（矩阵）描述。磁化率提高（JACS,90,811(1968)）主要用于实验测定芳香性。假设化合物的各向同性磁化率为A，而对应的电子定域化的异构体的各向同性磁化率为B，则A减B就是，体现了电子离域效应造成的额外的对磁屏蔽的效果。这里所指的磁化率是各项同性值，即磁化率张量对角元的平均值。越负代表电子离域效应使实际体系磁化率降低得越多，芳香性越强；接近0代表非芳香性；越正代表反芳香性越强。和NICS有很大雷同之处，物理本质也是相似的。从计算化学角度上看没什么价值，算起来也麻烦，算NICS就行了。（注意此方法名字中的“提高”是指抗磁化率大小因电子离域导致的增加，增量越大也就是指抗磁化率数值变得更负，总磁化率更小。有些文中用到时特意将名字中的Magnetic写为diamagnetic以免混淆。另外注意，此方法的JACS原文及其它一些文献中的的数值的正负号与如今习俗是相反的）抗磁化率的各向异性值自古就被用来衡量芳香性，因为其各向异性本质是由于环上电子离域性所产生的，和芳香性直接相关，这个概念的提出可以追溯到L.C.Pauling1936年的文章JCP,4,673。=_11-(_22+_33)/2。这里_11、_22、_33指的是矩阵的绝对值由大到小的三个本征值，数值越负表明芳香性越强。也可以通过实验或理论计算获得。但是由于实验上测定各向异性值需要先对化合物进行结晶，所以不如获得那么方便。顺磁化率/抗磁化率可以在Gaussian程序中通过NMR=Susceptibility计算。所得输出文件中会有诸如这样的抗磁化率信息Diamagneticsusceptibilitytensor(au)XX=-263.8815YX=0.0000ZX=0.0000XY=0.0000YY=-145.1938ZY=0.0000XZ=0.0000YZ=0.0000ZZ=-136.8484Isotropicdiamagneticsusceptibility=-181.9745由于当前体系是C2v对称性，且计算时开启了对称性，使得输出的抗磁化率张量直接就是对角矩阵，本征值就是对角元，故=-263.8815-(-145.1938-136.8484)/2=-122.8604a.u.。如果矩阵不是对角化的，需要自行计算矩阵的本征值再算。4基于几何结构的方法本节介绍的方法只依赖于几何结构信息，从键长的变化角度衡量芳香性。4.1Harmonicoscillatormeasureofaromaticity(HOMA)HOMA是最为常用也是历史十分悠久的基于几何结构的衡量芳香性的方法。最初是在TetrahedronLett.,36,3839(1972)中定义的，后来在Chem.Inf.Comput.Sci.,33,70(1992)中被广义化，表达式为其中N是环上的原子数，i循环环上的每个原子，j是它的下一个原子，R是键长。和R_ref都是事先计算好的常数，和键的类型有关，在原文中已提供，后者是理想的芳香性体系的平衡键长。HOMA的函数形式基于谐振势，这是名字中Harmonicoscillator的由来。如果环上的所有键都和相应的R_ref相同，显然体系就是理想的芳香性体系（例如苯），所以HOMA越接近1说明芳香性越强；若HOMA接近0，说明是非芳香性（的确定方式就是使kekule结构的苯HOMA值精确为0）；如果是负值，表现出键长极度不均衡，反映出反芳香性特征。还有一些研究者在HOMA上动一些手脚，例如JPCA,107,7496(2003)中就把HOMA的键长替换为离域化指数，定义为指数；再比如Bultinck等人在J.Phys.Org.Chem.,18,706(2005)中提出的BOIA(Bondorderindexofaromaticity)定义为BOIA=1-(1/6)(B_i-B_ref)，和HOMA主要差异也就是把键长替换为了键级。这类在HOMA上鬼鬼祟祟做一些小修改来得到新的芳香性指数的研究很无聊，非常没创意，和HOMA定性结果往往也大同小异，笔者也就不屑于将之纳入Multiwfn了。由于键级和模糊空间下的离域化指数在Multiwfn中都能算，因此和BOIA想算的话自己手算即可，也不麻烦。在Multiwfn中提供了计算HOMA的工具，并且涉及到的参数都可以自定义。还是以菲为例，启动Multiwfn后首先输入包含菲的结构的文件，fch/wfn/wfx/pdb/xyz等等类型文件都包含分子结构信息所以都可以。然后选100，再选13进入HOMA计算界面。如果想自定义参数，就选1。按0可以开始计算，用户需要输入环上的原子序号，例如3,4,8,9,10,7（对应菲的中央的环）。这些原子必须依次相邻，顺时针还是逆时针顺序都可以。之后结果马上输出出来，如AtompairContributionBondlength3(C)-4(C):-0.0656981.4271114(C)-8(C):-0.2104551.4580008(C)-9(C):-0.0656981.4271119(C)-10(C):-0.0960161.43528210(C)-7(C):-0.0336731.3600007(C)-3(C):-0.0960161.435282HOMAvalueis0.432442每个键的键长以及它对HOMA的贡献值都会输出。由于C4-C8偏离理想芳香键键长很大，导致HOMA值下降很多。接下来再计算边缘的环，结果是0.855126，由于比中央的环更接近于1，所以明显比中央的环芳香性更强。HOMA的一个缺点是预先提供的参数比较有限，包括C-C、C-N、C-O、C-P、C-S、N-N、N-O，如果有其它元素参与就没法用了。而且根本没法用于研究诸如三个乙炔聚合形成苯环过程、Diels-Alder加成反应过程中芳香性的动态变化。HOMA的一个主要优点是只要有分子结构就行，不依赖于电子结构信息，计算量可以忽略不计。4.2Bird芳香性指数Bird芳香性指数在Tetrahedron,41,1409(1985)当中被提出，表达式为其中加和符号循环环上的每个键，R是键长，N是Gordy关系式计算的键级。N上带一个横杠是环上的键的N值的算术平均值。n是键的数目，V_K是参考值，对于五元环和六元环分别为35和33.2。a和b是Gordy关系式预设的参数，在Tetrahedron,57,5715(2001)当中给出了新的a和b参数。Bird芳香性指数越接近100，说明键的均衡性越好，芳香性越强。Bird芳香性流行程度远不及HOMA，且比起HOMA适用的元素类型更少，且只适合五、六元环。在Multiwfn中计算Bird的步骤和计算HOMA的步骤几乎一样，唯一不同的是，进入主功能100里的功能13后，不是选择0，而是选择2。利用选项3可以自定义各个元素之间的键的a、b值，利用选项4可以自定义V_k值。5基于能量的方法本节介绍的方法基于能量标准来衡量芳香性。5.1离域化能及HOSE(Harmonicoscillatorstabilizationenergy)离域化能是指体系从定域化的电子结构（由lewis式描述）转化到实际的离域化的电子结构过程中的能量的变化。是否将几何结构的变化对能量影响考虑进去，又可以分为垂直和绝热离域化能。苯的离域化能就是从kekule式描述的三个双键+三个单键转变为实际的6个单键+1个六中心六电子大pi键的能量变化。芳香性越强，电子离域性越强，离域化能也就越大。离域化能这个概念很好，但是计算时很困难。一方面是以什么结构lewis结构作为参考标准的问题；另一方面更棘手，是用什么程序怎么实际去算的问题，下面谈几种可行的手段：(1)离域化能在VB框架中称为共振能(Resonanceenergy,RE)。价键理论(VB)用在计算离域化能很合用，因为VB计算时可以明确指定成键方式，因此能直接得到电子定域化状态的能量。然后把所有可能的lewis结构在计算时一起考虑，对它们各自的系数变分，得到的就是在这些定域化状态间相互共振的状态的能量，这等同于电子离域状态。将这个能量减去定域化状态的能量即得到共振能。不幸的是VB在目前是比较小众的计算方法，门槛略高。(2)NBO的二阶稳定化能(E2)分析。NBO是定域化的轨道，占据的NBO轨道与空的NBO轨道间的E2对应了电子从前者向后者的离域化造成的能量降低。把体系中各种与环上电子离域相关的E2都考虑进去，就能估计环上电子的离域化能。但是这个方法很不推荐，或者只能很定性地研究。因为E2不具备加和性，而且E2只是低阶的微扰理论，而芳香环上电子的离域远非“微”的级别。而且还有其它一些原理上的缺陷，这需要对NBO有一定认识才容易理解，这里就不多讨论了。另外还有人用类似Fock矩阵元消除法来获得稳定化能。(3)BLW(block-localizedwavefunction)方法。此方法原文见JCP,109,1687(1998)。这个方法是基于MO理论的，可以用在HF或DFT上。BLW方法是将体系的基函数划分为位于特定片段上的数个子集，在SCF过程中只允许MO在这些子空间里分别展开。于是每个子空间里的MO是正交的，而每个子空间内的MO与其它子空间的MO是非正交的。利用BLW方法可以得到与lewis式对应的定域化的波函数，并得到对应的能量。令体系实际能量减去BLW波函数的能量就是离域化能。BLW可以结合GAMESS-US和XMVB程序计算，见http:/wiki.lct.jussieu.fr/workshop/index.php/BLW。不过虞忠衡指出使用BLW计算离域化能有严重缺陷，见/blog-94786-423823.html。虞忠衡也提出了计算稳定化能的方法，见其blog上一系列文章，这里就不多提了。以上三种方法是基于电子结构的方法来计算垂直或绝热离域化能。而HOSE(Harmonicoscillatorstabilizationenergy)是基于几何结构的变化按照谐振势来计算绝热离域化能，其表达式如下式中R_r和R_r代表实际分子中的pi键键长，n1和n2分别代表kekule式对应结构的单键和双键的数目，k是力常数。因此HOSE描述了分子从实际几何结构变化到Kekule式描述的电子定域化的几何结构过程中所需的能量的负值，或者说是源自电子离域效应产生的实际结构相对于kekule结构的稳定化能。关于计算离域化能其它的一些细节讨论和历史上的方法可参见Chem.Rev.,105,3773(2005)。5.2芳香稳定化能(Aromaticstabilizationenergy,ASE)ASE和离域化能都是基于能量来分析芳香性的标准，体现电子离域导致的能量变化，但ASE的计算方法与之截然不同。ASE计算过程是构建一个开环的等键化学反应。反应能如果算出来是正值，即吸热，说明体系是芳香性；如果是负，即放热，说明是反芳香性。例如苯的ASE可按如下方式计算，能量是B3LYP/6-31G*下得到的一般认为ASE5kcal/mol就可以算是有芳香性了。计算ASE时构建的等键反应有任意性，不同的构建方式结果会有所不同，如果构建得不合适，甚至结果会不靠谱，详见ComputationalChemistry-IntroductiontotheTheoryandApplicationsofMolecularandQuantumMechanics2ed(ErrolLewars)一书307页的讨论。ASE最大的弱点就是计算极其麻烦。对于一个复杂的大体系，或者诸如金属芳香性体系，构建一个合理、有意义的等键反应几乎不可能。构建方式的不唯一性也导致不容易对不同文献的结果进行对比。所以我不认为如今使用ASE讨论芳香性有多大价值，除了对于简单小体系。5.3磁共振能(Magneticresonanceenergies,MRE)MRE是一种计算ASE和环电流的方法，详见JACS,128,2873(2006)、JPCA,111,8873(2007)。MSE基于图论和Huckel理论，只要知道分子拓扑关系即可进行计算，计算时只需要一个参数即Huckel共振积分。此方法只能用于pi芳香性研究。这个方法计算时要考虑体系全部可能的环路。比如菲不是只包含三个六元环路，而是包括三个六元环路、两个10元环路和一个14元环路。对于每个环路，比如第i个环路，可以计算相应的A_i值，它为正值和负值说明这个环路的感应环电流起抗磁（电流方向符合左手定则）和顺磁（电流方向与左手定则相反）作用。A_i也被称为环路共振能(CRE)，将当前体系的所有A_i加起来，得到的就是当前体系的MRE，可视为是离域化能。根据I_i=4.5*I_0*A_i*S_i/S_0可以计算第i个环路的环电流，其中I_0是苯的环电流，S_i是第i个环的面积，S_0是苯环的面积。将所有环路的环电流累加到一起，就能绘制出整体感应电流图。这和AICD那些方法得到的不一样，那些方法可以得到空间中每个具体的点感应电流密度，而MRE的电流图给出的只是每个键上电流大小和方向信息，如下所示。MRE方法的优点是不依赖于量化计算，不依赖于具体几何结构，只根据拓扑关系就能迅速得到离域化能和感应电流分布。而且不像大多数芳香性判据那样只能算每个小环的芳香性（比如对菲只能算三个六元环），而是可以研究体系中任意环路的稳定化能和环电流。此方法的缺点也很多，即只能用于平面pi体系、精度仅是HMO级别、不能用于非平衡构型或非基态或开壳层、体现不了键的拉伸产生的影响、不能用于诸如金属簇等体系等等。MRE最适合的体系就是多环芳烃（可以含杂原子），特别是尺寸很大的情况。实现MRE的程序据笔者所知目前没有公开的，不过估计也不会很难写。6基于电子离域性的方法本节介绍的方法都比较直接地与电子离域性相联系。6.1多中心键级多中心键级定义为其中P是单电子密度矩阵，S是重叠矩阵，a,b,c.代表基函数序号，A,B,C.代表环上的原子序号，并且这些原子是按照拓扑关系依次相邻的。PatrickBultinck等人的文章鼓吹应当把加和号内的项的各种置换组合都考虑进去，实际上根本没这个必要，分析结果和不考虑置换时一样，而计算量还增大了很多。多中心键级计算速度很快，理论依据强，是寡人强烈推荐使用的衡量芳香性的方法。它的数值越大说明芳香性越强。它也可以用在开壳层，只要把密度矩阵换成alpha和beta电子的密度矩阵分别计算，然后加和即可。由于对称性的原因，sigma和pi电子的多中心键级是可加和的，因此多中心键级可以分别研究sigma和pi芳香性。多中心键级的一个遗憾是没法用来分析反芳香性，反芳香性的体系只是多中心键级很小而已，故和芳香性很弱的体系没法区分开。另外，虽然多中心键级计算量不大，对于六元环瞬间就能出结果，但是计算量却随着环内原子数以阶乘方式增加，当环很大的时候，如9中心时，计算量就挺大了。需要注意的是计算多中心键级时切不可使用弥散函数，否则结果会缺乏意义。虽然也有诸如多中心离域化指数（Theor.Chem.Acc.,105,292(2001)）等类似的指标能够用来衡量芳香性，可以减小多中心键级的基组依赖性，尤其是能够解决不适合弥散函数的问题，但计算明显更费时。多中心键级可以在Multiwfn当中计算，目前最多支持六中心，下个版本支持到10中心。启动Multiwfn后输入Gaussian的fch文件的路径（不能用其它类型文件作为输入），然后选9，再选2，输入环上的原子序号（必须在环上依次相邻），如3,7,10,9,8,4，然后六中心键级数值会立刻输出（输入几个原子序号就计算几中心键级）。注意顺时针和逆时针输入环上的原子序号结果会有微小差异，原理上都对，取谁都可以，也可以取平均值。如果想计算仅由pi或仅由sigma电子贡献的多中心键级，则应先进入主功能100里的功能22。分子在哪个平面上就选哪个平面，然后程序会自动探测出哪些MO是pi轨道。如果接下来想计算sigma部分，就选1将pi轨道的占据数都设为0；如果想计算pi部分，就选2把除了pi轨道以外的MO占据数都设为0。之后再照常计算多中心键级，所用的密度矩阵就是扣掉了sigma或pi贡献的密度矩阵，也即得到了pi或sigma的多中心键级。（PS：设定轨道占据数也可以用主功能6的选项26，由此可以手动将指定的MO的占据数设为0来排除它的贡献。但是想排除sigma或pi轨道的贡献，用户就得自己去观察轨道图形来找出pi轨道，这比较费事，远不如直接用主功能100里的功能22那么方便。）6.2ELF-sigma/pi和LOL-sigma/pi电子定域化函数(Electronlocalizationfunction,ELF)是Becke提出的衡量电子定域性的实空间函数，目前有极为广泛的应用。此函数的在电子定域性的图形分析（/sobereva/item/4a8aeae2291d8d0f8c3ea8d8）中从图形的角度进行了介绍，在物理化学学报,27,2786(2011)中作者对ELF的函数形式和物理