文/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图像的等值面面积在某些情况下可以作为解释或预测相互作用能的很好的描述符，这一点很值得在未来进一步探索。

为了考察温度对碳环-富勒烯结合的影响，以及探索结合的临界温度，文中利用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程序做波函数分析对弱相互作用研究的关键性价值！

文/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原文，给别人代算时也必须明确告知对方这一点。

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

0 前言

键级是考察化学键特征的一类重要方法，非常强大的波函数分析程序Multiwfn可以方便地计算很多类型的键级，在《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）的键级部分有具体的说明，Multiwfn手册4.9节有一些实际计算例子，在“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/workshop/WFN_content.html）里做了特别全面详细的讲解和演示。另外，Multiwfn还支持NAdO方法，可以将原子间的模糊键级以轨道形式图形化展现，对于了解成键本质极其有用，介绍和实例看《使用键级密度(BOD)和自然适应性轨道(NAdO)图形化研究化学键》（//www.umsyar.com/535）。

以上博文都是举的孤立体系的例子，而如今Multiwfn已经将上述方法扩展到了对周期性体系的分析上。将Multiwfn与CP2K做周期性DFT计算产生的波函数相结合，可以很容易地用上述方法考察周期性体系的成键情况，在本文将通过一个层状的共价有机框架化合物（COF）体系，以及一个Pd表面吸附苯的体系，对具体做法进行演示，希望读者举一反三将本文介绍的方法运用到实际研究中。使用本文的方法计算的结果若用于发表，记得需要按照Multiwfn启动时的提示恰当引用Multiwfn的原文。

读者务必使用2024-Jun-24及以后更新的Multiwfn版本，否则情况与本文所述不符。Multiwfn可以在主页//www.umsyar.com/multiwfn免费下载。不了解Multiwfn者建议阅读《Multiwfn入门tips》（//www.umsyar.com/167）和《Multiwfn FAQ》（//www.umsyar.com/452）。本文要用到非常流行、高效且免费的第一性原理程序CP2K，笔者假定读者已具备了相关常识，推荐不熟悉者通过“北京科音CP2K第一性原理计算培训班”（http://www.keinsci.com/workshop/KFP_content.html）系统性学习。本文要使用Multiwfn创建CP2K输入文件，这在《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）里有简要介绍。本文用的CP2K是2024.1版。

下面的例子涉及到的所有文件都可以在//www.umsyar.com/attach/719/file.rar里得到。

1 单层COF体系

这部分示例的COF体系的cif文件是本文文件包里的COF_16371N2.cif，用GaussView显示的结构如下，可见晶胞里有两层。

1.1 准备输入文件

此体系两层之间的相互作用主要是pi-pi堆积，对电子结构影响微乎其微，因此只计算一层就够了，可以节约时间，尤其是能显著节约NAdO分析过程中计算原子重叠矩阵（AOM）的时间。于是在GaussView里删除一层，然后保存为COF.gjf。由于当前体系不是导体且晶胞边长约15埃已经不小了，所以可以不用扩胞就直接做gamma点计算。注意CP2K导出molden文件只支持gamma点的情况，这在《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）里已明确说了。

启动Multiwfn，然后输入
COF.gjf
cp2k  //进入产生CP2K输入文件的界面
[回车]  //输入文件名用默认的COF.inp
-2  //要求产生molden文件
-7  //设置周期性
XY  //仅在平行于COF的方向考虑周期性，而垂直于COF的方向用非周期性
-9  //其它设置
16  //要求将Kohn-Sham矩阵导出到.csr文件里。此文件在Multiwfn计算NAdO轨道能量时会用到
0  //返回
0  //产生输入文件

现在当前目录下就有了COF.inp，对应PBE/DZVP-MOLOPT-SR-GTH级别的单点计算。用CP2K运行之，当前目录下会产生记录波函数信息的COF-MOS-1_0.molden和记录KS矩阵的COF-KS_SPIN_1-1_0.csr。按照《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》所述，在molden文件开头插入以下内容定义晶胞信息和当前赝势基组描述的各元素的价电子数
 [Cell]
 14.94670000     0.00000000     0.00000000
 -7.47335000    12.94422190     0.00000000
  0.00000000     0.00000000     8.00000000
 [Nval]
 C 4
 N 5
 H 1

之后就可以基于这个molden格式的波函数文件做各种周期性体系的分析了。

1.2 计算Mayer键级

首先计算非常常用的Mayer键级，它衡量原子间等效的共享电子的对数。启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
9  //键级分析
1  //Mayer键级
Multiwfn首先计算重叠矩阵，然后立刻输出了Mayer键级，只有大于阈值的（由settings.ini里的bndordthres参数控制）才直接显示了出来。要想看所有原子间的Mayer键级可以再选择y导出键级矩阵然后查看相应的矩阵元。为了便于对照，下图显示了当前体系各个原子的序号，本文主要关注其中绿色虚线圈住的那些原子。

Multiwfn输出的一些有代表性的键的Mayer键级值如下
11(C )   19(C )    1.12454138
11(C )   22(N )    1.48268711
19(C )   48(C )    1.27625346
20(C )   50(C )    1.39369508
18(C )   48(C )    1.24722702
可见这些键的Mayer键级都不同程度地明显大于1.0，因此可以推测它们都不仅仅形成了sigma键，还有一定pi成份。在后文会计算pi键级和使用NAdO分析更进一步确认这一点。上面列出的C-C键键级可以体现出体系中哪种C-C键的强度更强。可见连接萘单元和C3N3单元的C11-C19的键是相对最弱的。

1.3 计算模糊键级（fuzzy bond order）

Mayer键级怕弥散函数，而模糊键级则没有这个问题，用于任意基组都可以。Mayer键级用于MOLOPT系列基组都是没问题的，不过这里作为例子，也计算一下模糊键级。由于计算模糊键级需要先构造AOM，对较大的体系很耗时，所以对于Mayer键级适用的情况建议优先用Mayer键级。注意对孤立体系Multiwfn计算模糊键级默认是基于Becke划分的，而对于周期性体系计算模糊键级默认是基于Hirshfeld划分的，因此默认设置下周期性体系的计算结果和孤立体系的没有严格的可比性。

在键级分析主功能的菜单里选择7，然后Multiwfn就会开始计算AOM，之后立刻给出模糊键级的计算结果，前述的那些键的模糊键级如下
11(C )   19(C )    1.05202861
11(C )   22(N )    1.50866214
19(C )   48(C )    1.22313834
20(C )   50(C )    1.38494183
18(C )   48(C )    1.17984661
虽然模糊键级与Mayer键级在定量上有一些差别，但不同的键的键级的大小顺序是完全一致的，没有结论上的差异。

1.4 计算pi电子贡献的Mayer键级

我在《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）中专门介绍了pi键级的计算方法，也即pi占据轨道对Mayer键级的贡献，没读过的话务必先阅读。这一节对单层COF这个周期性体系也做这个计算，看看各个键的pi作用程度的差异。由于当前体系是精确平行于XY平面的，所以让Multiwfn自动指认pi轨道前不需要先做轨道定域化。

启动Multiwfn并载入COF-MOS-1_0.molden后，依次输入
100  //其它功能（Part 1）
22  //检测pi轨道
0  //当前的轨道都是离域形式（如正则分子轨道）
2  //设置其它轨道占据数为0
0  //返回
9  //键级分析
1  //Mayer键级
前面提到的那些键的键级计算结果如下，这对应的就是pi键级。可见C3N3六元环中的C-N键的pi共享电子作用很强，比这里列出的C-C键的还要更显著。
11(C )   19(C )    0.14583409
11(C )   22(N )    0.41775145
19(C )   48(C )    0.28202057
20(C )   50(C )    0.38270736
18(C )   48(C )    0.27414145

使用此做法发表文章时，除了引用Multiwfn原文外还建议同时引用介绍Multiwfn做pi电子结构分析的文章Theor. Chem. Acc., 139, 25 (2020)。

1.5 NAdO分析

利用前述的《使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键》中介绍的NAdO分析，可以将sigma和pi作用分别以轨道形式呈现出来，非常有价值。例如在我的研究18氮环的ChemPhysChem (2024) https://doi.org/10.1002/cphc.202400377文中就用了NAdO分析非常清楚、严格地展示了较短的N-N键的sigma+pi作用特征。本节我们用NAdO分析展示一下C11-C19的sigma和pi轨道相互作用情况。

启动Multiwfn并载入COF-MOS-1_0.molden后，依次输入
15  //模糊空间分析
3  //计算AOM。对周期性体系默认是基于Hirshfeld原子空间计算的
[回车]  //用默认的0.2 Bohr格点间距。此设置对当前体系耗时不高，但对于明显更大体系为了节约耗时可以适当用更大的格点间距如0.35 Bohr，但会多多少少牺牲精度

计算完成后当前目录下出现了记录AOM的AOM.txt。在Multiwfn里接着输入
0  //返回
200  //其它功能（Part 2）
20  //BOD和NAdO分析
-1  //要求计算NAdO的能量
2  //从外部文件读取Fock矩阵
COF-KS_SPIN_1-1_0.csr  //输入此文件的实际路径
1  //基于AOM计算原子间的相互作用
[回车]  //读取当前目录下的AOM.txt
11,19  //分析11和19号原子间的相互作用
马上NAdO分析就做完了，屏幕上显示以下内容。这说明所有NAdO的本征值加和为1.052，和上一节得到的模糊键级一致。并且当前只有前两个NAdO的本征值明显大于0，它们是最值得考察的，其它的基本可以忽略

 Eigenvalues of NAdOs: (sum=   1.05204 )
   0.78955   0.19252   0.06784   0.02133   0.00807   0.00651   0.00127
   0.00096   0.00031   0.00011   0.00010   0.00005   0.00004   0.00000
...略

所有的NAdO轨道现在都已经导出到了当前目录下的NAdOs.mwfn中。现在输入y让Multiwfn载入此文件，之后退回到Multiwfn主菜单，进入主功能0查看NAdO轨道。不熟悉Multiwfn看轨道功能的话参考《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）。在图形界面里分别选择1号和2号轨道，可分别看到下面的图像，同时在文本窗口显示的NAdO本征值和轨道能量也标上了。这里等值面数值用的是0.08。

由上可见，C11-C19同时具有sigma和pi作用特征，前者远比后者显著得多。而且由于C19处在pi共轭区域，所以NAdO方法产生的主要描述C11-C19的pi作用的轨道并不完全定域在C11和C19上，也同时一定程度离域到了周围原子上。通过NAdO轨道能量还可以看到，C11和C19之间共享的sigma电子的能量明显比pi电子要低，这十分符合一般化学常识。

2 Pd(100)表面吸附苯

本节通过Pd(100)晶面吸附苯的体系，演示一下苯与Pd(100)基底这两个片段间Mayer键级的计算，以及对于像这样很大的体系怎么尽可能节约整个NAdO分析的时间。此体系的CP2K做结构优化的输入和输出文件，以及优化任务跑完后产生的restart文件都在本文文件包里提供了。优化完的结构如下，可见吸附作用很强，以至于Pd-C的作用都令苯产生弯曲了。

由于此例子的molden文件、FOM.txt和NAdOs.mwfn文件太大，本文的文件包里就没提供了。

2.1 准备输入文件

用Multiwfn载入本文文件包里的Pd+ben_opt-1.restart读取优化完的坐标和晶胞信息，然后输入以下命令创建一个做单点且产生molden文件的任务
cp2k  //产生CP2K输入文件
Pd+ben_SP.inp  //产生的输入文件名
-7  //设置周期性
XY
-2  //要求产生molden文件
6  //由于Pd基底是导体，开启smearing
0  //产生输入文件

为了帮助SCF收敛，把Pd+ben_SP.inp里的ALPHA设为0.15，NBROYDEN设为12。用CP2K计算Pd+ben_SP.inp，得到Pd+ben_SP-MOS-1_0.molden。然后在里面开头部分加入以下内容。

 [Cell]
 23.33880000     0.00000000     0.00000000
 0.00000000    23.33880000     0.00000000
 0.00000000     0.00000000    24.79750184
 [Nval]
 Pd 18
 C 4
 H 1

2.2 片段间Mayer键级

这一节将苯和Pd(100)基底分别定义为一个片段计算它们之间的Mayer键级。启动Multiwfn，载入Pd+ben_SP-MOS-1_0.molden，然后输入
9  //键级分析
-1  //定义片段
289-300  //苯的部分
1-288  //Pd(100)基底部分
1  //计算键级

Multiwfn开始产生重叠矩阵，之后给出了原子间键级，以下是其中一部分，可见苯中的C和有的Pd之间的Mayer键级较大，即共价作用很显著。

...略
 # 1412:       246(Pd)  270(Pd)    0.21887832
 # 1413:       246(Pd)  271(Pd)    0.21771061
 # 1414:       246(Pd)  289(C )    0.08497491
 # 1415:       246(Pd)  293(C )    0.09391700
 # 1416:       246(Pd)  294(C )    0.55236636
 # 1417:       247(Pd)  271(Pd)    0.23611928
 # 1418:       247(Pd)  272(Pd)    0.23841126
...略

老有人问怎么判断固体表面对分子是物理吸附还是化学吸附，判断方法很多，如结合能的数量级、原子间距离和原子半径的关系、电子密度差、IGMH（//www.umsyar.com/621）等等，而用Mayer或模糊键级是最能说明问题的。以上Mayer键级体现的C和Pd之间的明显的共价作用意味着二者之间形成了明显的化学键作用，显然Pd(100)对苯是化学吸附。如果键级只有比如零点零几的程度，那一般可以说是物理吸附，但也不排除形成典型离子键的可能，这通过原子电荷可以判断，计算方式可参考《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）。

最后Multiwfn还给出了片段间键级：
The bond order between fragment 1 and 2:    4.309784
即苯分子和Pd(100)之间等效共享电子对数多达4.3，显然太共价了，如果这都不叫化学吸附...

2.3 NAdO分析

这一节我们将Pd(100)与苯之间显著的共价作用通过NAdO方法以轨道形式直观展现。本节的做法和1.5节示例的常规的NAdO分析流程有两个明显的不同：
(1)对当前这样原子数又多、轨道数又多的大体系计算所有原子的AOM是相当耗时间的事情，而且导出、载入AOM.txt的耗时巨高而且巨占硬盘（正比于占据轨道数的平方乘以原子数）。为了极大程度节约时间，可以使用Multiwfn提供的一个特殊的策略，也就是在主功能15里用选项33计算感兴趣的两个片段间的片段重叠矩阵（FOM）并导出为FOM.txt，然后对两个片段之间做NAdO分析时直接从中读取要用的FOM，这样不仅省得导出和载入巨大的AOM.txt，而且产生FOM的过程中不需要计算片段里没有涉及的原子的AOM，当片段里的原子数占整个体系的原子数比例较少的时候这可以大幅节约时间。对于当前的例子，苯可以作为片段1，与苯相互作用最密切的那部分Pd原子适合作为片段2，这比起把所有Pd都定义为片段2能节约巨量耗时。但用户很难判断哪些Pd原子与苯作用密切，此时可以利用Multiwfn特意提供的一种定义片段2的方法，即片段1以外的原子中，与片段1任何一个原子之间Mayer键级大于指定阈值的原子都被定义为片段2。利用这个做法可以使得与苯的轨道相互作用最显著的Pd原子被定义为片段2，其数目远少于总的Pd数目。
(2)在计算FOM和做NAdO分析之前，需要先用主功能6的子功能38令电子根据轨道能量由低到高重新排列，使得轨道占据数都为整数。这是因为当前CP2K计算用了smearing，导致费米能级附近的轨道处于分数占据状态，这不是单行列式波函数应具有的特征，而Multiwfn的NAdO分析目前只支持单行列式波函数。重新排列后，当前的波函数就成了标准的单行列式波函数形式，从而可以做NAdO分析。另一方面，这么重排占据数后在计算AOM和FOM时可以只对占据轨道之间计算，也只有这部分是NAdO分析所实际需要的，只计算占据轨道之间的AOM/FOM的耗时、对内存的占用远远低于考虑所有轨道的情况，

启动Multiwfn，载入Pd+ben_SP-MOS-1_0.molden，然后输入
6  //修改波函数
38  //按照轨道能量由低到高重排占据数
-1  //返回
15  //模糊空间分析
33  //计算FOM
3  //定义两个片段并计算FOM，不属于片段1的原子若与片段1的任意原子的Mayer键级大于指定值则被定义为片段2
289-300  //苯的部分
0.001  //Mayer键级的阈值

之后看到以下信息，告诉了你哪些原子被作为了片段2，并且可见片段1和2之间的Mayer键级为4.114，和上一节看到的苯与整个Pd(100)表面的片段间Mayer键级4.310相差很小，这说明当前自动定义的片段2是合理的。如果你把片段2的原子在Multiwfn主功能0里用菜单栏的Other settings - Set atom highlighting功能高亮显示，会看到它们都是离苯不远的原子。

 Atoms in fragment 2:
 76,77,80-83,86,87,97,101-103,106-108,112,195-197,200-202,205-207,220,221,224-22
 6,230,231,241,245-247,250-252,256,265,266,269-272,274-277,280,281
 Mayer bond order between fragments 1 and 2:   4.114

之后输入0.35，Multiwfn就开始用0.35 Bohr间距的立方格点计算AOM，然后构造出FOM并导出到当前目录下的FOM.txt。在双路7R32机子上96核并行计算，整个过程花了26分钟（耗时与格点间距的三次方呈反比，耗时太高的话可以用比如0.5 Bohr格点间距，误差也还可以接受，此时只需要不到10分钟）。然后输入
0  //返回
200  //其它功能（Part 2）
20  //BOD/NAdO分析
4  //直接从当前目录下的FOM.txt中读取NAdO分析要用的片段1和2的FOM
之后开始了NAdO的计算，结果如下

 Eigenvalues of NAdOs: (sum=   6.95512 )
   0.79476   0.77393   0.71745   0.63686   0.62698   0.35965   0.22325
   0.16243   0.14986   0.14655   0.14545   0.13818   0.12618   0.12058
   0.11388   0.11078   0.10617   0.09335   0.09007   0.08979   0.08423
   0.08224   0.07942   0.07434   0.07391   0.07171   0.06475   0.05751
   0.05507   0.05478   0.05246   0.04964   0.04669   0.04030   0.03771
...略

当前显示的所有NAdO本征的加和6.955对应于苯和Pd(100)之间在Hirshfeld原子空间划分下算的模糊键级，和前面看到的片段间Mayer键级4.31有一定差异，这很正常，毕竟对原子空间的定义截然不同。由以上数据可见有不少NAdO轨道对模糊键级的贡献都显著，尤其是前六个。现在输入y载入新产生的NAdOs.mwfn，在主功能0里观看NAdO轨道，其中本征值最大的6个NAdO轨道如下所示，等值面数值用的0.03

从上图中可以看到这些轨道在苯与Pd(100)之间都是相位相同方式叠加的，必然都是起到明显成键作用的，所以本征值都不小。从轨道图形上可以看出这些轨道来自于苯上C原子的垂直于苯环的p轨道与Pd的原子轨道混合，而且有的图里直接就能清楚看出Pd用的是d原子轨道。例如NAdO(2)中和苯接触的Pd明显用的是dz2轨道，从Pd区域的轨道等值面形状就能看出这一点。感兴趣的话可以进一步按照《谈谈轨道成份的计算方法》（//www.umsyar.com/131）介绍的用SCPA方法做一下轨道成分分析。在Multiwfn主菜单里输入
8  //轨道成分分析
3  //SCPA方法
2  //2号NAdO轨道

Multiwfn马上就输出了轨道成份，下面列出的是Multiwfn返回的各个原子上各个角动量基函数产生的贡献。明显可以看到Pd主要都是用D基函数，而C用的都是P基函数。因此对NAdO(2)的轨道成份的分析体现了Pd的d原子轨道与C的p原子轨道的显著混合对Pd吸附苯有关键性贡献。可见以这种方式讨论成键特征能分析得巨清楚透彻！

Composition of each shell
Shell  1370 Type: D    in atom  196(Pd) :     0.50345 %
Shell  1405 Type: D    in atom  201(Pd) :     0.88530 %
Shell  1440 Type: D    in atom  206(Pd) :     0.50251 %
Shell  1538 Type: D    in atom  220(Pd) :     0.80192 %
Shell  1545 Type: D    in atom  221(Pd) :     0.76004 %
Shell  1570 Type: S    in atom  225(Pd) :     0.53572 %
Shell  1573 Type: D    in atom  225(Pd) :    10.40186 %
Shell  1577 Type: S    in atom  226(Pd) :     0.56745 %
Shell  1580 Type: D    in atom  226(Pd) :    10.21554 %
Shell  1608 Type: D    in atom  230(Pd) :     0.80255 %
Shell  1615 Type: D    in atom  231(Pd) :     0.76058 %
Shell  1713 Type: D    in atom  245(Pd) :     0.53415 %
Shell  1717 Type: S    in atom  246(Pd) :     0.96339 %
Shell  1720 Type: D    in atom  246(Pd) :    21.79157 %
Shell  1748 Type: D    in atom  250(Pd) :     0.53493 %
Shell  1752 Type: S    in atom  251(Pd) :     0.96056 %
Shell  1755 Type: D    in atom  251(Pd) :    21.80729 %
Shell  2019 Type: P    in atom  289(C ) :     1.48815 %
Shell  2024 Type: P    in atom  290(C ) :     1.47584 %
Shell  2029 Type: P    in atom  291(C ) :     6.18539 %
Shell  2034 Type: P    in atom  292(C ) :     1.44475 %
Shell  2039 Type: P    in atom  293(C ) :     1.41651 %
Shell  2044 Type: P    in atom  294(C ) :     6.20361 %

感兴趣的读者还可以用CDA方法从苯的分子轨道与Pd表面的晶体轨道的混合角度来考察Pd对苯吸附造成的电子转移和轨道相互作用的细节。详见《使用Multiwfn结合CP2K对周期性体系做电荷分解分析（CDA）》（//www.umsyar.com/716）。

3 总结

本文通过COF二维层状体系和Pd(100)表面吸附苯作为例子，演示了如何用Multiwfn程序计算原子间和片段间的Mayer键级和模糊键级，还演示了如何计算pi键级单独考察pi电子的贡献，以及讲解了如何做NAdO分析以轨道形式清楚直观地考察原子间或片段间的共价作用的内在特征。可见使用Multiwfn可以把成键特征情况展现得超级清楚透彻。本文的做法对其它情况的周期性体系，诸如三维原子晶体、过渡态结构等等，也都是完全适用的。更多细节和相关知识请参看本文提到的相关博文。

Multiwfn对周期性体系能算的键级不止本文示例的这些。对周期性体系Multiwfn还能算Mulliken键级和Wiberg键级，还能把Mulliken和Mayer键级分解成不同轨道的贡献，后者称为轨道占据数扰动的Mayer键级。具体介绍见《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471），做法见Multiwfn手册4.8节的相关例子，本文就不特意演示了。除了本文介绍的这些外Multiwfn还有很多其它考察成键的方法对周期性体系都可以用，比如AIM拓扑分析，见《使用Multiwfn结合CP2K做周期性体系的atom-in-molecules (AIM)拓扑分析》（//www.umsyar.com/717）。

使用Multiwfn计算的结果若用于发表，记得需要按照Multiwfn启动时的提示恰当引用Multiwfn的原文。给别人代算的话也必须明确告知对方这一点。

文/Sobereva @北京科音  2024-Jul-5
不得转载

1 前言

此文对非常适合理论、计算化学投稿的期刊（也有少部分仅限于约稿），以及化学相关的综合向以及面向不同子领域的重点期刊进行罗列，共约300种期刊，并附上2024年6月公布的2023年影响因子(IF)及中科院分区。本文对于读者了解有哪些期刊值得投稿和平时关注有明显意义。这是第15次笔者写此系列的文章，往年的文章见本文文末的链接。

本文的期刊名称一般用标准缩写，有的给出了常用首字母缩写或者全名，之后是ISSN号（若只有网络版则是EISSN号）。斜杠前是2023年的两年期影响因子（即一般所谓的影响因子。有些新刊尚无影响因子数据），斜杠后是撰写此文时的最新的中科院大类分区（2023年12月公布），今年刚有IF的期刊目前没有分区信息。

附IF计算公式：2017年某期刊的影响因子=（2015+2016年此期刊的文章在2017年被引用的次数）/（2015+2016年此期刊的文章数）

本文涉及的刊物主要分这么几类：
1 主要理论、计算类刊物
2 综合
3 中国的期刊
4 无机、有机、化学信息学及其它
5 偏材料
6 偏物理
7 偏生物
第一类是最适合理论、计算类文章投的。根据文章具体的研究方向，也可选择投其它几类中的刊物。

本文中给出的期刊列表一般是按照IF来排序的，有IF的都是SCI期刊。国内核心非SCI期刊没有IF，这些期刊的影响因子的获取方式为：进入http://navi.cnki.net/knavi/Journal.html，搜索期刊名，取综合影响因子。本文里期刊普遍用缩写，如果想得到全称，去https://www.ablesci.com/journal把缩写或ISSN输进去就可以查到。

像往年本系列文章一样，我再次强调，IF这东西只是个翔，它的存在给科学发展带来的坏处绝对大于好处。要看重文章自身的价值，垃圾文章侥幸发到高IF期刊上依然是垃圾文章，好文章发到低IF期刊上依然是好文章，会被广为引用。本文中说的期刊档次/质量仅是指的平均文章档次/质量。而如今国内大学和研究机构评判科研业绩所广为采用的中科院分区更是极度荒诞的，其存在的弊远远大于利。

2 今年形势简评

今年不管是哪个方向的期刊，跌的都明显比涨的多。能不跌就等效于轻微上涨了。

这是和理论、计算化学关系最密切的一些期刊的IF的变化情况：
大涨：WCMS、CHEM PHYS IMPACT
小涨：JCC、CTC、ELECT STRUCT、STRUCT CHEM、J STRUCT CHEM+
基本没变：JCIM、JCTC、JPCA、JMGM、CPL、IJQC、JMM、J COMPUT BIOPHYS CHEM、MOL SIMULAT、JPOC、THEOR CHEM ACC、MOL PHYS
小跌：JPCB、JPCC、PCCP、ATS、CHEM PHYS
大跌：JPCL、JCP、CPC

计算化学综述期刊WCMS从11.4猛升到16.8，这令我有点意外，近些年WCMS上的文章的平均水准也就那么回事。JCP从4.4跌倒3.1，跌幅惊人，特别是对于JCP这样刊登文章量非常大的期刊来说。原因一部分可能是2020、2021年那会儿集中刊登的计算化学程序原文对IF的贡献已过，因此又打回原形了。JPCL继续下跌，从5.7跌到4.8，如今比JCTC和JCIM都低不少了。本身从现在JPCL上文章的题材和整体水准来看也不像是IF能很高的情况。PCCP挺惨，居然跌破3了，显得比之前低了一个level，虽然学术口碑还在。CPC从原来的3左右跌到2.3真是巨惨，也显得变低了一个level。CPC上的文章整体水准还是挺不错的，和JPCA半斤八两，这个IF实在容易让人低估CPC的档次。在很多老牌计算化学相关期刊IF都下跌的大环境下，JCC能从3升到3.4，十分难得，现在反超了PCCP和JCP。CTC的IF居然进一步上升，到了3，这个给人传统印象是三线的期刊如今能超过PCCP是挺难以置信的事，也进一步说明现在的IF越来越不能反映期刊整体水准和真实的业界认可度了。ELECT STRUCT小升到2.9，看来这个新刊的地位已经比较稳固了。Chemical Physics衍生出的开放访问的Chemical Physics Impact竟然从2.2升到3.8，很离谱，达到了Chemical Physics的近两倍，不知道怎么运作的。本来这个期刊很大程度上是Chemical Physics拒稿的回收站。开放访问的吸金期刊ACS Physical Chemistry Au第一次有了IF，3.7，一般般，才和回收站ACS Omega正好一样。

化学综合方面：
NATURE、SCIENCE、CHEM REV都大跌超过10。CHEM SOC REV也跌不少，而NAT REV CHEM微涨，二者IF相仿佛了。ACS CENTRAL SCI巨降，从JACS明显上头跌到明显下头了。JACS Au现在还是和JACS差得很远。跌幅最狠的是SCI DATA从9.8降到5.8，估计是什么极高被引的数据集文章对IF的贡献过期了。Frontiers in Chemistry从5.5大降到3.8。CHEM COMMUN继续降，只有4.3了，当年六点几的时代已经回不去了，而且比后起之秀的开放访问的COMMUN CHEM还落后了。Scientific Reports已降到3.8，ACS OMEGA跌到3.7，RSC Advances 3.9未变，这三个大容量开放访问期刊的IF现在惊人地接近。SPECTROCHIM ACTA A只小降0.1到4.3，在大环境都在跌的时代能保住这个IF难得了。口碑一直挺好的CHEM-EUR J小跌到3.9，终于到了4以下，可惜了。而其姊妹刊CHEM-ASIAN J跌到3.5，继续保持着和它的差距。Chemistry-Methods首次有了IF，6.1，很不错，我的IRI原文（Chemistry—Methods, 1, 231 (2021)）给此刊贡献了不少引用，是此刊目前发表的文章里引用第一名。值得一提的是，一些大出版社的创刊不久的开放访问的期刊现在也纷纷成了SCI，如mdpi的Computation和Chemistry，Elsevier的Results in Chemistry、Springer的SN Applied Sciences。另外伊朗化学会的Physical Chemistry Research和Sami Publishing Company的Chemical Methodologies也都成了SCI。

中国期刊方面：
《中国科学:化学（英文版）》小升而CCS Chem小降，前者超过了后者。《结构化学》竟然从2.2猛升到5.7，涨了一倍多，很离谱，莫非也要走《物理化学学报》巨幅提升IF的道路？《物理化学学报》基本没变，今年10.8，能守住10以上不容易。《化学学报》从2.5跌到了1.7，这是为数不多的只接受中文文章的化学期刊，在此呼吁中国研究者们多多支持。《高等学校化学学报》跌到只有0.7了，没能像许多其它中国创办的化学期刊那样从当年的1左右开始在近年来逐渐攀升。

无机、有机、化学信息学及其它：
INORG CHEM FRONT跌到6.1。J MOL LIQ跌到5.3。J MOL STRUCT居然小升到4了，此刊上很多文章巨水，JMS给我送审的有的纯计算的文章写得特别业余，我觉得发JMM都悬，现在IF却能升到4，又是个典型的IF不能反映期刊实际水准的例子。小涨的INORG CHEM COMMUN已超过小跌的INORG CHEM。J ORG CHEM连年下降，继续跌，只有3.3了，惨。不错的ORGANOMETALLICS现在跌到只有2.5了，这世道...知名的ACTA CRYSTALLOGR B居然跌到只剩1.3，太惨了。

材料方面：IF升的寥寥无几，而且都只是小涨，90%以上都下跌，而且有一大批都跌得很整齐，跌幅都在0.7-1左右。新刊ACS Materials Au现在有了IF，5.7，中规中矩。

生物方面：同样几乎都在跌，而NUCLEIC ACIDS RES却涨了1.7，难得。

3 期刊列表

(1) 主要理论、计算类刊物

WIREs Comput Mol Sci (WCMS) 1759-0876 16.8/2
J CHEM THEORY COMPUT (JCTC) 1549-9618 5.7/1
J CHEM INF MODEL (JCIM) 1549-9596 5.6/2
J PHYS CHEM LETT (JPCL) 1948-7185 4.8/2
J PHYS CHEM C (JPCC) 1932-7447 3.3/3
J PHYS CHEM B (JPCB) 1520-6106 2.8/2
J PHYS CHEM A (JPCA) 1089-5639 2.7/2
Chemical Physics Impact 2667-0224 3.8/无分区 2020创刊，免费阅览，发表收费
ACS Physical Chemistry Au 2694-2445 3.7 2021创刊，免费阅览，发表收费
J COMPUT CHEM (JCC) 0192-8651 3.4/3
J CHEM PHYS (JCP) 0021-9606 3.1/2
COMPUT THEOR CHEM (CTC) 2210-271X 3/3
PHYS CHEM CHEM PHYS (PCCP) 1463-9076 2.9/3
ADV THEORY SIMUL (ATS) 2513-0390 2.9/4
Electronic Structure 2516-1075 2.9/无分区 2019创刊
CHEM PHYS LETT (CPL) 0009-2614 2.8/3
J MOL GRAPH MODEL (JMGM) 1093-3263 2.7/4
CHEM PHYS CHEM (CPC) 1439-4235 2.3/3
INT J QUANTUM CHEM (IJQC) 0020-7608 2.3/3
J MOL MODEL (JMM) 1610-2940 2.1/4
STRUCT CHEM 1040-0400 2.1/4
CHEM PHYS 0301-0104 2/3
J COMPUT BIOPHYS CHEM 2737-4165（以前叫J THEOR COMPUT CHEM (JTCC)） 2/4
MOL SIMULAT 0892-7022 1.9/4
J PHYS ORG CHEM (JPOC) 0894-3230 1.9/4
THEOR CHEM ACC (TCA) 1432-881X 1.6/4
MOL PHYS 0026-8976 1.6/4
J STRUCT CHEM+ (JSC, Journal of Structural Chemistry) 0022-4766 1.2/4
ADV QUANTUM CHEM 0065-3276 ?/无分区
    以下期刊非SCI
Computational Chemistry 2332-5968 2013创刊。免费阅览，发表收费 http://www.scirp.org/journal/cc/
International Journal of Computational and Theoretical Chemistry (IJCTC) 2376-7286 2013创刊。免费阅览，发表收费 http://www.sciencepublishinggroup.com/j/ijctc
Communications in Computational Chemistry (CiCC) 2305-7076 2013创刊（此刊已偃旗息鼓，2018年之后未更新） http://www.global-sci.org/cicc/
SDRP Journal of Computational Chemistry & Molecular Modelling (JCCMM) 2473-6260 2015创刊。免费阅览，发表收费（网站自称IF是0.827） http://www.siftdesk.org/journal-details/SDRP-Journal-of-Computational-Chemistry-&-Molecular-Modelling-/33
Living Journal of Computational Molecular Science (LiveCoMS) 2575-6524 2017创刊。免费阅览 https://www.livecomsjournal.org/
Turkish Computational and Theoretical Chemistry (TC&TC) 2017创刊。免费阅览，免费发表（黑白图片时） https://dergipark.org.tr/en/pub/tcandtc
Journal of Molecular Physics 2017创刊。免费阅览 https://scholars.direct/journal.php?jid=molecular-physics
Chemical Physics Reviews 2688-4070 2021创刊。AIP创办
Journal of Physical Chemistry & Biophysics 2161-0398 2011创刊。免费阅览，免费发表 https://www.longdom.org/physical-chemistry-biophysics.html

(2) 综合

CHEM REV 0009-2665 51.4/1
NATURE 0028-0836 50.5/1
SCIENCE 0036-8075 44.7/1
CHEM SOC REV 0306-0012 40.4/1
NAT REV CHEM (Nature Reviews Chemistry) 2397-3358 38.1/1
Nature Chemistry 1755-4330 19.2/1
CHEM (Elsevier旗下) 2451-9294 19.1/1
ACCOUNTS CHEM RES (ACR) 0001-4842 16.4/1
NATL SCI REV (National Science Review) 2095-5138 16.3/1
ANGEW CHEM INT EDIT 1433-7851 16.1/1
NAT COMMUN 2041-1723 14.7/1 免费阅览，发表收费
J AM CHEM SOC (JACS) 0002-7863 14.4/1
ADV SCI (Advanced Science。Wiley旗下) 2198-3844 14.3/1 免费阅览，发表收费
Trends in Chemistry 2589-5974 14/2 2019创刊 Cell出版社
ACS CENTRAL SCI 2374-7943 12.7/1 免费阅览，免费发表
ANNU REV PHYS CHEM 0066-426X 11.7/1
SCI ADV (Science Advances。Science旗下) 2375-2548 11.7/1 免费阅览，发表收费
P NATL ACAD SCI USA (PNAS) 0027-8424 9.4/1 发表收费
JACS Au 2691-3704 8.5/无分区 2020创刊，免费阅览，发表收费
Research 2639-5274 8.5/1
Cell Rep. Phys. Sci. 2666-3864 7.9/2 免费阅览，发表收费
CHEM SCI (Chemical Science) 2041-6520 7.6/1 免费阅览，免费发表
Chemistry-Methods 2628-9725 6.1/无分区 2020创刊，Wiley旗下。免费阅览，发表收费
COMMUN CHEM (Communications Chemistry) 2399-3669 5.9/2 2018创刊 Nature出版社。免费阅览，发表收费
SCI DATA (Scientific Data) 2052-4463 5.8/2 免费阅览，发表收费
Arabian J. Chem. 1878-5352 5.3/2
INT J MOL SCI (IJMS, International Journal of Molecular Sciences) 1422-0067 4.9/2 免费阅览，发表收费
iScience 2589-0042 4.6/2 Cell旗下，2018创刊。免费阅览，发表收费
CHEM COMMUN 1359-7345 4.3/2
BMC Chemistry 2661-801X 4.3/2 免费阅览，发表收费
SPECTROCHIM ACTA A 1386-1425 4.3/2
Molecules 1420-3049 4.2/2 免费阅览，发表收费
Dyes and Pigments 0143-7208 4.1/3
CHEM-EUR J (Chemistry-A European Journal) 0947-6539 3.9/2
RSC Advances 2046-2069 3.9/3 免费阅览，发表收费
SCI REP (Scientific Reports) 2045-2322 3.8/2 免费阅览，发表收费
Frontiers in Chemistry 2296-2646 3.8/3 有理论化学版块。免费阅览，发表收费
ACS OMEGA 2470-1343 3.7/3 免费阅览，发表收费
CHEM-ASIAN J (Chemistry-An Asian Journal) 1861-4728 3.5/3
Chemical Methodologies 2645-7776 3.5/无分区 2017创刊。免费阅览，发表收费。http://www.chemmethod.com
Heliyon 2405-8440 3.4/3 免费阅览，发表收费。全学科
FARADAY DISCUSS 1364-5498 3.3/3
J COMPUT SCI 1877-7503 3.1/3
Royal Society Open Science 2054-5703 2.9/3 2014创刊。免费阅览，发表收费
SN Applied Sciences 2523-3963 2.8/无分区 2019创刊，Springer旗下。免费阅览，发表收费
NEW J CHEM 1144-0546 2.7/3
INT REV PHYS CHEM 0144-235X 2.5/2
ChemistryOPEN 2191-1363 2.5/4 免费阅览，发表收费
Results in Chemistry 2211-7156 2.5/无分区 2019创刊，Elsevier旗下。免费阅览，发表收费
SoftwareX 2352-7110 2.4/4 免费阅览，发表收费，专门发表免费程序的介绍文章
Chemistry 2624-8549 2.4/无分区 2019创刊，MDPI旗下。免费阅览，发表收费 https://www.mdpi.com/journal/chemistry
ISRAEL JOURNAL OF CHEMISTRY 0021-2148 2.3/4
B KOREAN CHEM SOC 0253-2964 2.3/4
J IRAN CHEM SOC 1735-207X 2.2/4
J COMPUT ELECTRON 1569-8025 2.2/4
OPEN CHEM 2391-5420 2.1/4 免费阅览，发表收费
CHEMICAL PAPERS 0366-6352 2.1/4
ChemistrySelect 2365-6549 1.9/4
Computation 2079-3197 1.9/无分区 MDPI旗下。免费阅览，发表收费 https://www.mdpi.com/journal/computation
J CHEM SCI 0974-3626 1.7/4
Physical Chemistry Research 2322-5521 1.4/无分区 伊朗化学会创办，2013创刊。免费阅览。http://www.physchemres.org
CAN J CHEM 0008-4042 1.1/4
AUST J CHEM 0004-9425 1/4
RUSS J PHYS CHEM A+ 0036-0244 0.7/4
CROAT CHEM ACTA 0011-1643 0.7/4
INDIAN J CHEM 0019-5103 0.4/无分区
  以下期刊非SCI
Nature Computational Science 2662-8457 2021创刊，Nature旗下
Chemical Physics Reviews 2688-4070 2020创刊，AIP旗下
Electronic Materials 2673-3978 2021创刊，MDPI旗下。免费阅览，发表收费
Physical Sciences Reviews 2365-659X 2016创刊 只接受邀请稿件 https://www.degruyter.com/journal/key/psr/html?lang=en
Current Physical Chemistry 1877-9468 2011创刊，http://benthamscience.com/journal/index.php?journalID=cpc
General Chemistry 2414-3421 2015创刊，免费阅览和发表 http://www.genchemistry.org
Natural Science 2150-4091 2009创刊 免费阅览，发表收费 http://www.scirp.org/journal/ns/
Open Journal of Physical Chemistry 2162-1969 2011创刊。免费阅览，发表收费。http://www.scirp.org/journal/ojpc/
International Journal of Chemistry 1916-9698 2012创刊。免费阅览，发表收费。http://www.ccsenet.org/journal/index.php/ijc
American Journal of Physical Chemistry 2327-2430 2012创刊。免费阅览，发表收费。http://www.sciencepublishinggroup.com/j/ajpc
Science Journal of Chemistry (SJC) 2330-0981 2013创刊。免费阅览，发表收费。http://www.sciencepublishinggroup.com/j/sjc
American Journal of Chemistry and Application (AASCIT) 2375-3765 2014创刊。免费阅览，发表收费。http://www.aascit.org/journal/about?journalId=905
American Journal of Chemistry 2165-8749 2011创刊。免费阅览，发表收费。http://www.sapub.org/Journal/articles.aspx?journalid=1091
Physical Chemistry 2167-7042 2011创刊。免费阅览，发表收费。http://www.sapub.org/journal/articles.aspx?journalid=1022
Journal of Atomic and Molecular Sciences (JAMS) 2075-1303 2010创刊。http://www.global-sci.org/jams/
Chemical Review and Letters 2645-4947 2018创刊。发表和阅览都免费。http://chemrevlett.com（网站自称IF是0.92）
Journal of Open Research Software (JORS) 2049-9647 2013创刊，专门收轻量级程序介绍文章，免费阅览，可申请免费发表。https://openresearchsoftware.metajnl.com
Journal of Open Source Software (JOSS) 2475-9066 专门收开源程序的轻量级介绍文章，免费阅览和发表。https://joss.theoj.org
Chemical Reports 2591-7943 2019创刊。免费阅览，发表收费 https://www.syncsci.com/journal/CR/about
F1000 Research 2046-1402 免费阅览，发表收费 https://f1000research.com/
Journal of Chemistry: Education Research and Practice(JCERP) 2578-7365 （网站自称IF是0.94） 免费阅览，发表收费 https://www.opastpublishers.com/journal/journal-of-chemistry-education-research-and-practice
PeerJ Physical Chemistry 2689-7733 免费阅览，发表收费 https://peerj.com/physical-chemistry/

(3) 中国的期刊

中国期刊英文版：
Science Bulletin 2095-9273 18.8/1 全学科，收版面费
中国科学:化学（英文版）SCI CHINA CHEM 1674-7291 10.4/1。以前叫SCI CHINA SER B
CCS Chem 2096-5745 9.4/无分区 2019创刊，免费访问，免费发表
中国化学快报 CHINESE CHEM LETT 1001-8417 9.4/1
结构化学 CHINESE J STRUC CHEM 0254-5861 5.9/4
Fundamental Research 2667-3258 5.7/3 2021创刊。免费阅览，发表收费。国家自然科学基金委创办
中国化学 Chinese Journal of Chemistry 1001-604X 5.5/1
高等学校化学研究 CHEM RES CHINESE U 1005-9040 3.1/4
中国化学会会志（台湾） J CHIN CHEM SOC-TAIP 0009-4536 1.6/4
化学物理学报 CHINESE J CHEM PHYS 1674-0068 1.2/4

中国期刊（中文为主）：
物理化学学报 ACTA PHYS-CHIM SIN 1000-6818 10.8/2
有机化学 CHINESE J ORG CHEM 0253-2786 1.8/4
化学学报 Acta Chim Sinica 0567-7351 1.7/4
化学进展 PROG CHEM 1005-281X 1/4
物理学报 Acta Phys Sinica 1000-3290 0.8/4
高等学校化学学报 CHEM J CHINESE U 0251-0790 0.7/4
  以下期刊非SCI
中国科学:化学 Scientia Sinica Chimica 1674-7224 （以前叫 中国科学B）0.789 核心
化学通报 Chemistry 0441-3776 0.745 核心
化学研究与应用 Chemical Research and Application 1004-1656 0.741 核心
分子科学学报 Journal of Molecular Science 1000-9035 0.359 核心
原子与分子物理学报 Journal of Atomic and Molecular Physics 1000-0364 0.341 核心
计算机与应用化学 Computers and Applied Chemistry 1001-4160 非核心
物理化学进展 Journal of Advances in Physical Chemistry 2168-6122 非核心，2012创刊。免费阅览，发表收费 http://www.hanspub.org/journal/japc

(4) 无机、有机、化学信息学及其它

COORDIN CHEM REV 0010-8545 20.3/1
INORG CHEM FRONT 2052-1553 6.1/1
J MOL LIQ 0167-7322 5.3/2
ORG CHEM FRONT 2052-4129 4.6/1
ORG LETT 1523-7060 4.9/1
INORG CHEM 0020-1669 4.3/2
Dalton Transactions 1477-9226 3.5/3
LANGMUIR 0743-7463 3.7/2
J MOL STRUCT 0022-2860 4/2
INORG CHEM COMMUN 1387-7003 4.4/3
J ORG CHEM (JOC) 0022-3263 3.3/2
MOL INFORM 1868-1743 (以前叫QSAR & Combinatorial Science 1611-020X) 2.8/4
J COMPUT AID MOL DES (Journal of Computer-Aided Molecular Design) 0920-654X 3/3
J CHEM EDUC (JCE) 0021-9584 2.5/3
Inorganics 2304-6740 3.1/4 免费阅览，发表收费
ORGANOMETALLICS 0276-7333 2.5/3
EUR J ORG CHEM 1434-193X 2.5/3
POLYHEDRON 0277-5387 2.4/3
MATCH-COMMUN MATH CO (MATCH Communications in Mathematical and in Computer Chemistry) 0340-6253 2.9/2
EUR J INORG CHEM 1434-1948 2.2/4
J ORGANOMET CHEM 0022-328X 2.1/3
TETRAHEDRON 0040-4020 2.1/3
ACTA CRYSTALLOGR B 2052-5206 1.3/3
TETRAHEDRON LETT 0040-4039 1.5/4
EUR PHYS J D 1434-6060 1.5/4
J MATH CHEM 0259-9791 1.7/3
MACROMOL THEOR SIMUL 1022-1344 1.8/4
LETT ORG CHEM 1570-1786 0.7/4
TETRAHEDRON CHEM 2666-951X 2022创刊。免费阅览，发表收费。非SCI
Organics 2673-401X 2020创刊，MDPI旗下。免费阅览，发表收费。非SCI
Journal of Computer Chemistry, Japan -International Edition (JCCJIE) 2189-048X 非SCI，免费阅览，发表收费

(5) 偏材料

Nature Materials 1476-1122 37.2/1
ADV MATER 0935-9648 27.4/1
InfoMat 2567-3165 22.7/1 免费阅览，前三年不收发表费
MATER TODAY 1369-7021 21.1/1
ADV FUNCT MATER 1616-301X 18.5/1
ACS Nano 1936-0851 15.8/1
Small 1613-6810 13/2
Mater Horiz 2051-6347 12.2/2
J Mater Chem A 2050-7488 10.7/2
Carbon 0008-6223 10.5/2
Nano Lett 1530-6984 9.6/1
NPJ COMPUT MATER 2057-3960 9.4/1 免费阅览
J MATERIOMICS 2352-8478 8.4/1 免费阅览
Communications Materials 2662-4443 7.5/无分区
CHEM MATER 0897-4756 7.2/2
Science China Materials 2095-8226 6.8/2
Applied Surface Science 0169-4332 6.3/2
J Mater Chem B 2050-750X 6.1/3
MATER CHEM FRONT 2052-1537 6/2
Nanoscale 2040-3364 5.8/3
J Mater Chem C 2050-7526 5.7/2
ACS Materials Au 2694-246 5.7/无分区 2021创刊，ACS旗下，免费阅览，发表收费。非SCI
APL Materials 2166-532X 5.3/2 免费阅览，发表收费
Materials Advances 2633-5409 5.2/无分区 2020创刊，RSC旗下。免费阅览，发表收费。非SCI
2D Materials 2053-1583 4.5/3
Nanomaterials 2079-4991 4.4/3
MATER CHEM PHYS 0254-0584 4.3/3
MATER TODAY COMMUN 2352-4928 3.7/3
J MATER SCI 0022-2461 3.5/3
J Solid State Chem 0022-4596 3.2/3
COMP MATER SCI 0927-0256 3.1/3
ORG ELECTRON 1566-1199 2.7/4
Computational Condensed Matter 2352-2143 2.6/无分区 2014创刊，Elsevier旗下。非SCI
FRONT MATER 2296-8016 2.6/4
Physica Status Solidi (RRL) - Rapid Research Letters 1862-6254 2.5/4
J PHYS-CONDENS MAT 0953-8984 2.3/4
Surface Science 0039-6028 2.1/4
MODEL SIMUL MATER SC (MODELLING AND SIMULATION IN MATERIALS SCIENCE AND ENGINEERING) 0965-0393 1.9/4
Bulletin of Materials Science 0250-4707 1.9/4
Chemistry of Inorganic Materials 2949-7469 2023创刊，Elsevier旗下，免费阅览，发表收费。非SCI

(6) 偏物理

Reviews of Modern Physics (RMP) 0034-6861 45.9/1
Nature Physics 1745-2473 17.6/1
Physical Review X (PRX) 2160-3308 11.6/1
PHYS REV LETT (PRL) 0031-9007 8.1/1
COMPUT PHYS COMMUN 0010-4655 7.2/2
JPhys Materials 2515-7639 无/3
J COMPUT PHYS 0021-9991 3.8/2
PHYSICAL REVIEW RESEARCH 3.5/无分区 2019创刊，免费阅览，发表收费
PHYS REV B (PRB) 2469-9950 3.2/2
PHYS REV A (PRA) 2469-9926 2.6/2
PHYS REV E (PRE) 2470-0045 2.2/3
Physica B: Physics of Condensed Matter 0921-4526 2.8/3
CHINESE PHYS B 1674-1056 1.5/4
AIP Advances 2158-3226 1.4/4 免费阅览，发表收费

(7) 偏生物

NUCLEIC ACIDS RES (NAR) 0305-1048 16.6/2
EMBO J 0261-4189 9.4/1
PLOS BIOL 1544-9173 7.8/1
J MOL BIOL 0022-2836 4.7/2
COMPUT STRUCT BIOTEC 2001-0370 4.4/2 免费阅览，发表收费
METHODS 1046-2023 4.2/3
J BIOL CHEM (JBC) 0021-9258 4/2
PLoS Comput Biol 1553-734X 3.8/2
ACS CHEM BIOL 1554-8929 3.5/2
BIOPHYS CHEM 0301-4622 3.3/3
BIOPHYS J 0006-3495 3.2/3
PROTEINS (proteins: Structure, Function, and Bioinformatics) 0887-3585 3.2/4
BIOCHEMISTRY (ACS的) 0006-2960 2.9/3
PLoS One 1932-6203 2.9/3 发表收费
ORG BIOMOL CHEM 1477-0520 2.9/3
BBA Biomembranes 0005-2736 2.8/3
J BIOMOL STRUCT DYN 0739-1102 2.7/3
COMPUT BIOL CHEM (Computational Biology and Chemistry) 1476-9271 2.6/4 2002年及以前叫Computers & chemistry
BBA-PROTEINS PROTEOM 1570-9639 2.5/4
J BIOL PHYS 0092-0606 1.8/4
J COMPUT BIOL 1066-5277 1.4/4
Journal of Biophysical Chemistry (JBPC) 2153-036X 非SCI，2010创刊。免费阅览，发表收费 http://www.scirp.org/journal/jbpc
J Org Biomol Simul 2325-2170 非SCI，2013创刊。免费阅览，发表收费 http://thescipub.com/jobs.toc
Computational Molecular Bioscience 2165-3445 非SCI，2011创刊。免费阅览，发表收费 http://www.scirp.org/journal/cmb/

附：历年的影响因子和点评

适合理论、计算化学投稿的期刊及其2022年影响因子（2023年公布）
//www.umsyar.com/675
适合理论、计算化学投稿的期刊及其2021年影响因子（2022年公布）
//www.umsyar.com/646
适合理论、计算化学投稿的期刊及其2020年影响因子（2021年公布）
//www.umsyar.com/603
适合理论、计算化学投稿的期刊及其2019年影响因子（2020年公布）
//www.umsyar.com/560
适合理论、计算化学投稿的期刊及其2018年影响因子（2019年公布）
//www.umsyar.com/492
适合理论、计算化学投稿的期刊及其2017年影响因子（2018年公布）
//www.umsyar.com/427
适合理论、计算化学投稿的期刊及其2016年影响因子（2017年公布）
//www.umsyar.com/382
适合理论、计算化学投稿的期刊及其2015年影响因子（2016年公布）
//www.umsyar.com/335
适合理论、计算化学投稿的期刊及其2014年影响因子（2015年公布）
//www.umsyar.com/296
适合理论、计算化学投稿的期刊及其2013年影响因子（2014年公布）
//www.umsyar.com/248
适合理论、计算化学投稿的期刊及其2012年影响因子（2013年公布）
//www.umsyar.com/192
适合理论、计算化学投稿的期刊及其2011年影响因子（2012年公布）
//www.umsyar.com/149
适合理论、计算化学投稿的期刊及其2010年影响因子（2011年公布）
//www.umsyar.com/92
适合理论、计算化学投稿的期刊及其2009年影响因子（2010年公布）
//www.umsyar.com/64

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

0 前言

Bader提出的非常知名的atom-in-molecules (AIM)理论在分析孤立体系和周期性体系的电子结构方面有非常广泛的应用，相关资料汇总见《AIM学习资料和重要文献合集》（http://bbs.keinsci.com/thread-362-1-1.html）。AIM分析框架中最常见的分析手段之一是电子密度的拓扑分析，也常被称为AIM拓扑分析，例如根据这种拓扑分析得到的键临界点（BCP）的位置可以用来讨论成键和弱相互作用特征，参考《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）和《Multiwfn支持的弱相互作用的分析方法概览》（//www.umsyar.com/252）里的相关介绍。非常强大的波函数分析程序Multiwfn（//www.umsyar.com/multiwfn）很早以前就已经支持了结合Gaussian、ORCA等 程序产生的波函数文件对孤立体系做AIM拓扑分析，见《使用Multiwfn做拓扑分析以及计算孤对电子角度》（//www.umsyar.com/108）和Multiwfn手册4.2节的大量例子。如今Multiwfn也已支持了对周期性体系做AIM拓扑分析，本文的目的是对此进行演示，以使得读者可以轻松地举一反三充分将这种手段应用于实际研究中。使用Multiwfn做AIM拓扑分析在发表文章时记得需要按照Multiwfn启动时的提示恰当引用Multiwfn的原文。

读者务必使用2024-Jun-16及以后更新的Multiwfn版本，否则情况与本文所述不符。Multiwfn可以在主页//www.umsyar.com/multiwfn免费下载。不了解Multiwfn者建议阅读《Multiwfn入门tips》（//www.umsyar.com/167）和《Multiwfn FAQ》（//www.umsyar.com/452）。本文要用到非常流行、高效且免费的第一性原理程序CP2K，不熟悉者推荐通过“北京科音CP2K第一性原理计算培训班”（http://www.keinsci.com/workshop/KFP_content.html）系统性学习。本文要使用Multiwfn创建CP2K输入文件，这在《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）里有简要介绍。本文用的CP2K是2024.1版。

下面的例子涉及到的所有文件都可以在//www.umsyar.com/attach/717/file.rar里得到。

1 实例1：NaCl晶体

这一节对NaCl做AIM拓扑分析，要使用CP2K产生的molden文件作为Multiwfn的输入文件。molden文件的产生方法在《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）里有详细说明。由于考虑k点时CP2K没法产生molden文件，因此必须使用足够大的超胞模型并只考虑gamma点。本文文件包里NaCl.cif是NaCl单胞的晶体结构文件，基于它我们创建CP2K算单点任务的超胞的输入文件，同时得到molden文件。

用Multiwfn载入NaCl.cif，然后输入
cp2k  //创建CP2K输入文件的界面
[回车]  //产生的输入文件名用默认的NaCl.inp
-2  //要求产生molden文件
-11  //几何操作界面
19  //构造超胞
2  //第1个方向平移复制为原先的2倍
2  //第2个方向平移复制为原先的2倍
2  //第3个方向平移复制为原先的2倍
-10  //返回
0  //产生输入文件

现在当前目录下就有了NaCl.inp，对应的是PBE/DZVP-MOLOPT-SR-GTH计算级别，对此体系是合适的。当前晶胞边长为11.28埃，对于绝缘体来说只考虑gamma点基本够用。用CP2K运行NaCl.inp，得到NaCl-MOS-1_0.molden。用文本编辑器打开NaCl.inp，将以下内容插入文件开头来定义晶胞信息从而使得Multiwfn在分析时能考虑周期性，并且定义Na和Cl的价电子数分别为9和7（不了解为什么这么设置的人看前述的《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》）。
 [Cell]
11.28112000     0.00000000     0.00000000
 0.00000000    11.28112000     0.00000000
 0.00000000     0.00000000    11.28112000
 [Nval]
 Na 9
 Cl 7

下面开始做AIM拓扑分析。启动Multiwfn，载入NaCl-MOS-1_0.molden，然后输入
2  //拓扑分析（默认是对电子密度做拓扑分析）
2  //以每个原子核位置为初猜搜索临界点

当前从屏幕上会看到找到了64个(3,-3)临界点，分别对应64个原子核位置的电子密度极大点。之所以虽然当前用的是赝势基组，在原子核位置还能找到电子密度极大点，是因为Multiwfn自动对这些原子指认了合适的EDF信息来描述内核电子，详见《在赝势下做波函数分析的一些说明》（//www.umsyar.com/156）。

然后再输入
3  //以每两个原子间连线中点作为搜索临界点的初猜位置
4  //以每三个原子的中心位置作为搜索临界点的初猜位置
0  //观看结果

此时文本窗口显示了各种临界点的数目，如下所示，并且提示Poincare-Hopf关系已经满足了，因此大概率所有临界点都已经找全了。
 The number of critical points of each type:
 (3,-3):    64,   (3,-1):   384,   (3,+1):   384,   (3,+3):    64
 Poincare-Hopf verification:    64  -   384  +   384  -    64  =     0
 Fine, Poincare-Hopf relationship is satisfied, all CPs may have been found

在Multiwfn的图形窗口可以看到下图，棕红色大球是Na，绿色大球是Cl，橙色、黄色、绿色小球分别是电子密度的(3,-1)、(3,+1)和(3,+3)临界点，而(3,-3)临界点都被原子球挡住了。在Multiwfn手册3.14节对Multiwfn的拓扑分析功能有完整、非常详细的说明，建议阅读。

下面再把AIM拓扑分析里常涉及的键径产生出来。在图形界面右上角点Return关闭之，然后选8产生(3,-3)和(3,-1)之间的拓扑路径，也即键径。在笔者的i9-13980HX CPU上24核并行很快就产生完了。之后再选0进入图形界面观看，看到下图，可见所有键径（棕色细线）都产生了出来，既连接相邻的Na和Cl，也连接每一对临近的两个Cl。

默认情况下Multiwfn把处于晶胞边缘的所有镜像原子和临界点也显示了出来，如果不想显示它们的话（也即只显示晶胞中唯一的那些），在图形界面顶端依次选择Other settings - Toggle showing all boundary atoms和Toggle showing all boundary CPs and paths。在文本窗口也会看到当前设置的状态。

下面来考察一下Na和Cl之间的BCP的属性。用选项0的图形界面右边的选项把其它临界点关闭显示而只显示(3,-1)临界点，选择CP labels复选框要求显示临界点序号，关闭键径显示并要求显示原子，然后就能看到下图，这是当前图像的一个边缘区域（点击save picture按钮可在当前目录下得到高分辨率图像文件，下图是剪切出的一部分）

可见189、512、238等序号的临界点都对应Na-Cl的BCP，考察哪个都可以，都是等价的。点击Return按钮关闭图形窗口，然后选择7进入临界点属性计算界面，然后再输入临界点序号189，此时就看到了这个临界点处各种函数的值，如下所示。这些函数在Multiwfn手册2.6节都有简要介绍，且在“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/workshop/WFN_content.html）里有十分全面深入详细的讲解。

CP type: (3,-1)
 Density of all electrons:  0.1144132659E-01
 Density of Alpha electrons:  0.5720663296E-02
 Density of Beta electrons:  0.5720663296E-02
 Spin density of electrons:  0.0000000000E+00
 Lagrangian kinetic energy G(r):  0.1168592222E-01
 G(r) in X,Y,Z:  0.7296255010E-03  0.1022667133E-01  0.7296253875E-03
 Hamiltonian kinetic energy K(r): -0.3025433987E-02
 Potential energy density V(r): -0.8660488230E-02
 Energy density E(r) or H(r):  0.3025433987E-02
 Laplacian of electron density:  0.5884598270E-01
 Electron localization function (ELF):  0.1993303578E-01
 Localized orbital locator (LOL):  0.1249064363E+00
...略

如我在全面提到的《Multiwfn支持的分析化学键的方法一览》文中所述，BCP位置上电子密度拉普拉斯函数和能量密度都为正是离子键的典型特征。当前BCP这两个值分别为0.0588 a.u.和0.00302 a.u.，都为正值，完全体现出Na-Cl是离子键的实事。此外，当前BCP的电子密度仅为0.01144 a.u.、ELF仅为0.0199，BCP处很小的电子密度和ELF值也是典型的离子键的普遍共性。

我在《使用Multiwfn+VMD快速地绘制高质量AIM拓扑分析图（含视频演示）》（//www.umsyar.com/445）介绍了一种极好的将Multiwfn和VMD联用快速、理想地绘制AIM拓扑分析图的做法，比在Multiwfn里直接显示的效果更好，而且视角可以自由旋转，没看过的话务必先看一下。对于当前体系，也可以像此文一样用Multiwfn文件包里的examples\scripts\AIM.bat和AIM.txt来作图，得到的图像如下。

上图效果明显不好，显得残缺不全，这是因为此时显示的原子、临界点、键径都只有晶胞里唯一的那些，而处于棱和面上的边界原子、临界点、键径的周期镜像都没有显示出来。为了避免这个问题，对于当前体系有原子、临界点、键径处于棱和面上的情况，应当改用examples\scripts\目录下的AIM+boundary.bat和AIM+boundary.txt来绘制，其用法与AIM.bat和AIM.txt的组合完全一致。只不过前者传递给Multiwfn的指令要求Multiwfn导出mol.pdb、CPs.pdb和paths.pdb文件时把边界的镜像也都纳入进去。利用AIM+boundary.bat和AIM+boundary.txt来绘制的结果如下，可见此时的效果就很好了，和前面在Multiwfn里看到的完全相同。注意此时这三个pdb文件里只有序号靠前的那些才是晶胞中唯一的那些。

2 层状COF

下面再举个例子，用Multiwfn对一个层状的共价有机框架化合物（COF）体系做AIM拓扑分析。此体系的cif文件是本文文件包里的COF-12000N2.cif，其结构图如下所示。

由于此体系的晶胞尺寸不小，计算时可以直接对原胞做只考虑gamma点的计算来得到molden文件。如《实验测定分子结构的方法以及将实验结构用于 计算需要注意的问题》（//www.umsyar.com/569）所说，由于X光衍射测的晶体结构中氢的原子位置通常不准确，最好做分析前先冻结重原子而优化一下氢原子的坐标，但此例为了省事就跳过这一步了。如果你要用AIM拓扑分析考察氢键，由于氢的位置是重中之重，那显然是绝对不能略过这一步的。

当前体系用PBE结合CP2K的GTH赝势基组计算完全可以，而此例作为演示，这回使用PBE结合pob-TZVP-rev2全电子基组计算。启动Multiwfn，载入COF-12000N2.cif，然后输入
cp2k  //创建CP2K输入文件的界面
[回车]  //产生的输入文件名用默认的COF-12000N2.inp
-2  //要求产生molden文件
2  //选择基组
15  //pob-TZVP-rev2
4  //开OT，对当前用的基组通常比对角化明显更容易收敛
0  //产生输入文件

用CP2K计算COF-12000N2.inp，然后将以下晶胞设置信息插入到计算得到的COF-12000N2-MOS-1_0.molden文件的开头
 [Cell]
 22.55600000     0.00000000     0.00000000
-11.27800000    19.53406901     0.00000000
  0.00000000     0.00000000     6.80000000

启动Multiwfn，载入COF-12000N2-MOS-1_0.molden，然后输入
2  //拓扑分析（默认是对电子密度做拓扑分析）
2  //以每个原子核位置为初猜位置搜索临界点
3  //以每两个原子间连线中点作为搜索临界点的初猜位置
4  //以每三个原子的中心位置作为搜索临界点的初猜位置
8  //产生(3,-3)到(3,-1)的拓扑路径
0  //观看结果
现在看到下图。可见临界点和键径都成功地得到了，不仅有化学键对应的BCP和键径，COF内氢键N-H...O作用以及层之间相互作用相对应的BCP和键径都有。而且不仅有晶胞内两层COF之间相互作用的BCP和键径，当前晶胞里的COF和相邻镜像晶胞里的COF相互作用对应的BCP和键径也都有，显然周期性完全充分考虑了。

虽然Multiwfn的文本窗口提示当前并未符合Poincare-Hopf关系，但由于感兴趣的临界点和键径都有了，所以就没必要进一步搜索了。

再利用前述名为AIM+boundary.bat的Windows批处理脚本结合AIM+boundary.txt里记录的指令，用Multiwfn+VMD绘制体系结构+临界点+键径图。然后再在VMD的文本窗口里输入pbc box命令显示盒子，VMD main界面里选择Display - Orthographic用正交视角。此时VMD显示的图像如下，可见非常理想！

顺带一提，对于当前体系要想完整展现层间相互作用的话，最理想的做法是用我提出的IGMH，见《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）和《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（//www.umsyar.com/667）。

3 对Ag晶胞基于电子密度格点数据做AIM拓扑分析

Multiwfn不仅可以基于CP2K产生的波函数文件做AIM拓扑分析，还可以基于从cube等文件里载入的电子密度的格点数据做拓扑分析，此时各个位置的电子密度、梯度和Hessian矩阵基于格点数据以B-spline方式插值得到。由于这个特性，Multiwfn做周期性体系的AIM拓扑分析不限于CP2K用户，诸如从VASP产生的CHGCAR文件里读取电子密度格点数据做AIM拓扑分析也可以。Multiwfn支持哪些格点数据文件格式见手册2.5节的说明。而且CP2K用户还可以在考虑k点的情况下对原胞计算得到电子密度cube文件结合Multiwfn做AIM拓扑分析，而不必在晶胞较小时非得弄成超胞模型再算。但基于电子密度格点数据做AIM拓扑分析的关键不足是对临界点只能得到电子密度及其导数信息以及一些与之直接相关的量，而依赖于波函数才能计算的如动能密度、ELF等函数就没法计算了。另外，如果电子密度格点数据是在赝势下计算产生的，基于它插值出的电子密度函数中，原子核处将没有电子密度极大点，也因此无法产生键径。

此例演示用Multiwfn对Ag的单胞基于电子密度cube文件做AIM拓扑分析的过程。使用的级别是PBE/TZVPP-MOLOPT-PBE-GTH-q19，这个基组是CP2K 2024.1的data目录下的BASIS_MOLOPT_UZH基组文件里的。之所以不用常用的BASIS_MOLOPT里的Ag的基组，在于那里面Ag的基组都是-q11的，即只描述11个价电子。使用赝势的情况下做AIM拓扑分析时，被描述的价电子越多，即使用越小核赝势，结果越准确、越可能接近全电子计算的结果。

启动Multiwfn，载入本文文件包里的Ag单胞的晶体结构文件Ag.cif，然后输入
cp2k  //创建CP2K输入文件的界面
[回车]  //产生的输入文件名用默认的Ag.inp
-3  //要求产生cube文件
1  //电子密度
6  //开smearing
8  //设置k点
10,10,10
0  //产生输入文件

把Ag.inp里BASIS_SET_FILE_NAME设为BASIS_MOLOPT_UZH，POTENTIAL_FILE_NAME设为POTENTIAL_UZH。Ag的BASIS_SET设为TZVPP-MOLOPT-PBE-GTH-q19，POTENTIAL设为GTH-PBE-q19。然后用CP2K计算，之后会得到记录晶胞内价电子密度的Ag-ELECTRON_DENSITY-1_0.cube。如果不了解cube格式的话，看《Gaussian型cube文件简介及读、写方法和简单应用》（//www.umsyar.com/125）。用文本编辑器打开这个cube文件，在第一行里写上box2cell。这样Multiwfn载入这个cube文件时，一旦发现第一行的内容包含box2cell，就会把格点数据的盒子信息转化为晶胞信息，从而在计算的时候能够考虑周期性（也可以不修改cube文件，而是载入cube文件后用隐藏的主功能1000里的子功能18来实现盒子信息向晶胞信息的转化）。

把Multiwfn的settings.ini里的iuserfunc设为-3，使得user-defined function对应基于载入的格点数据通过B-spline插值的函数。然后启动Multiwfn，载入Ag-ELECTRON_DENSITY-1_0.cube，之后输入
2  //拓扑分析
-11  //修改被分析的函数
100  //User-defined function，当前对应于插值方式产生的价层电子密度
3  //将每两个原子之间的中点作为搜索临界点的初猜
然后进入选项0，并把原子标签显示出来、把键设细，就会看到下图。可见临近原子之间的BCP都已经找到了

下面再把更多的临界点也找出来。当前情况不适合选择选项4、5分别用3、4个原子中心作为初猜点来搜索临界点。选项3用每两个原子之间中点当做初猜点时是考虑循环相邻晶胞中的镜像原子的，而选项4、5则只会循环当前晶胞里那些唯一的原子，因此此例用选项4、5搜索后还是会遗漏许多重要的临界点。当前最适合的方式是选择选项6里的子选项-1，此时会在每个原子附近特定半径内撒一批初猜点来搜索临界点。这么选过之后，会发现找到了一大堆新的临界点，进入选项0后看到的图如下所示。下图左侧把所有临界点都显示了出来，可以看到离原子核较近的区域临界点有不少，而且各种类型的都有。这是因为当前用了赝势，由于电子密度不是从原子核开始向四周单调下降的，因此会在价层区域出现各种各样的临界点，这些不是我们感兴趣的。把原子球设大掩盖这些没什么意义的临界点，并且旋转体系使得视角与晶轴平行后，就看到了下面的右图。可见在原子间相互作用区域找到的临界点的分布整齐、对称，因此可以认为有意义的临界点都找到了。

有一个概念叫平面性指数，在“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/workshop/WFN_content.html）里有一页幻灯片做了介绍，如下所示。我们这里来对Ag晶体计算一下这个平面性指数，看看是否和文章里说的情况吻合。注意这里说的平面性指数和《使用Multiwfn定量化和图形化考察分子的平面性（planarity）》（//www.umsyar.com/618）里介绍的我提出的衡量体系几何结构平面性的参数是两码事。

首先我们得确定要考察的笼临界点也即(3,+3)的临界点编号。下图左边只显示了笼临界点，可见当前体系有两类笼临界点，一种是诸如晶胞中心的那个（15号），它处于八面体中心。另一种是处于四面体中心的，比如63号临界点。下图右边只显示了所有BCP，此体系中所有的BCP都是等价的，我们考虑其中任意一个，比如1号。利用当前拓扑分析界面里的选项7依次对这三个临界点的属性进行考察，可以得到的这三个位置的基于格点数据插值产生的电子密度，也即User-defined real space function后面的值，结果如下（注意不要读Multiwfn直接输出的Density of electrons后面的值，那个是用孤立状态原子的电子密度叠加得到的准分子密度）
CP(15): 0.1572533472E-01 a.u.
CP(63): 0.2212960297E-01 a.u.
CP(1): 0.2706526462E-01 a.u.

可见两类笼临界点里电子密度最小值是0.01572 a.u.，BCP处的电子密度为0.02706 a.u.，因此平面性指数为0.01572/0.02706=0.58，在0.5左右，符合幻灯片里说的“其它金属及合金”的情况。

4 总结&其它

本文通过三个有代表性的例子详细介绍了在Multiwfn中对周期性体系做电子密度的拓扑分析，也即一般说的AIM拓扑分析的完整过程，所有要注意的细节都做了解释。通过例子可见，Multiwfn做周期性体系的拓扑分析的过程很简单，并不比分析孤立体系更复杂。读者如果实际操作一遍，还会感受到Multiwfn做这些分析的速度相当快。AIM拓扑分析对于考察周期性体系的电子结构有重要意义，推荐读者在实际研究中运用，使文章增光添彩。

很值得一提的是，Multiwfn的拓扑分析功能极其普适，绝不仅限于电子密度的拓扑分析。对于Multiwfn自身可以直接计算的函数，以及载入的格点数据里记录的函数，Multiwfn都可以做拓扑分析，对分子体系的例子参看比如《使用Multiwfn对静电势和范德华势做拓扑分析精确得到极小点位置和数值》（//www.umsyar.com/645）对静电势、范德华势的拓扑分析，以及《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）里面对ELF-pi函数的拓扑分析。对周期性体系，Multiwfn也完全可以以相同的方式对任意函数做拓扑分析。

使用本文的做法做AIM拓扑分析在发表文章时必须按照Multiwfn启动时的提示恰当引用Multiwfn，如果用于给别人代算也要明确告知对方这一点。