文/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-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，如果用于给别人代算也要明确告知对方这一点。

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

1 前言

笔者之前写过《使用Multiwfn做电荷分解分析(CDA)、绘制轨道相互作用图》（//www.umsyar.com/166）介绍了怎么用Multiwfn程序基于Gaussian等 程序产生的波函数文件对分子体系做电荷分解分析（CDA）。这种方法可以从片段轨道间相互作用的角度深入了解片段间是如何转移电子的。由于CDA的很高价值，Multiwfn做CDA分析已经得到非常广泛的使用。如果读者没看过此//www.umsyar.com/166的话一定要先仔细看一下，否则无法理解下文的内容。

从2024-Jun-5更新的Multiwfn开始，CDA模块已经支持了基于CP2K产生的molden文件做周期性体系的CDA分析，下面将通过一个简单的例子进行演示。Multiwfn可以在主页//www.umsyar.com/multiwfn免费下载，不了解此程序者参看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。我假定读者已经具备了CP2K的基本使用知识，不了解者推荐通过北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）完整系统地学习。

本文的例子是计算NaCl板吸附CO分子的体系，用CDA方法考察CO与NaCl板之间电子转移情况。这个体系优化后的结构如下，C和O分别是1和2号原子，3-110号原子是NaCl板。对这个体系，之前我在《使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线》（//www.umsyar.com/638）中通过密度差图的方式定性考察了电子转移，在《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）里通过片段电荷的方式定量考察了电子转移。而下面用CDA分析，可以以轨道角度在明显更深层次理解电子转移细节。

本文例子涉及的所有文件都在//www.umsyar.com/attach/716/file.7z里。CP2K用的版本是2024.1，Multiwfn是2024-Jun-5更新的版本，读者绝对不要用Multiwfn的更老版本。

2 产生输入文件

首先我们要用CP2K计算NaCl+CO、NaCl、CO各自的molden格式的波函数文件。如果你不知道怎么用CP2K产生molden文件的话，先仔细阅读《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651），我假定读者已经看过本文了。这里有几个要点：
(1)整体和片段的molden文件里的原子坐标必须对应
(2)整体和片段的计算必须使用严格相同的级别和设置，用的盒子也应当一致
(3)必须要求CP2K计算所有空轨道
(4)整体的molden文件里需要加入[Cell]字段定义晶胞信息，而片段的molden文件里加不加入不影响CDA结果。[Nval]字段可加入可不加入，也不影响CDA结果

本文文件包里的NaCl_CO.cif是以前我在PBE-D3(BJ)/DZVP-MOLOPT-SR-GTH下优化的NaCl板吸附CO的结构文件。用Multiwfn载入它，然后输入
cp2k  //创建CP2K输入文件
NaCl_CO.inp  //产生的输入文件名
-2  //要求产生molden文件
-9  //其它设置
12  //设置计算的空轨道数
-1  //计算所有空轨道
0  //返回
0  //产生输入文件

用CP2K运行NaCl_CO.inp，得到NaCl_CO-MOS-1_0.molden。然后将以下字段插入到molden文件的开头
[Cell]
16.92168000     0.00000000     0.00000000
 0.00000000    16.92168000     0.00000000
 0.00000000     0.00000000    25.00000000

将NaCl_CO.inp复制为NaCl.inp，再将NaCl.inp中&COORD里的CO部分删掉、PROJECT NaCl_CO改名为PROJECT NaCl，然后用CP2K运行之，得到NaCl-MOS-1_0.molden。将上面的[Cell]加入其中（这和CDA分析无关，加入这个的目的是为了之后在Multiwfn中能正确观看NaCl板的周期性的晶体轨道）。

将NaCl_CO.inp复制为CO.inp，再将CO.inp中&COORD里的NaCl部分删掉、PROJECT NaCl_CO改名为PROJECT CO，然后用CP2K运行之，得到CO-MOS-1_0.molden。这个可以不用刻意加上[Cell]，因为当前体系中CO离盒子边界很远，显示轨道时是否考虑周期性对看到的图像没影响。

接下来就可以进行CDA分析了。

3 进行CDA分析

启动Multiwfn，载入NaCl_CO-MOS-1_0.molden，然后依次输入
16  //CDA分析
2  //两个片段
CO-MOS-1_0.molden  //CO的波函数文件。注意要按整个体系里原子序号顺序载入片段
NaCl-MOS-1_0.molden  //NaCl板的波函数文件

现在屏幕上输出了一大堆复合物轨道的CDA结果，同时还输出了ECDA结果：
The net electrons obtained by frag. 2 = CT( 1-> 2) - CT( 2-> 1) =    0.1580
这个ECDA结果告诉你CO（片段1）往NaCl板（片段2）转移了0.158个电子，这和《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）里通过Mulliken方法计算的CO的片段电荷严格相同。

当前CDA分析输出的复合物轨道太多，一共437个（对应整个体系的占据轨道数），而绝大多数对d、b、r项的贡献微乎其微或完全为0，没有考察的意义。为了简化输出，在当前界面里输入
-3  //设置CDA输出时d、b、r项的阈值，它们之一的绝对值大于这个阈值的复合物轨道才会输出
0.005
0  //显示CDA和ECDA结果

当前结果如下
   Orb.      Occ.          d           b        d - b          r
    272    2.000000    0.006220   -0.000012    0.006232    0.000488
    275    2.000000    0.125862   -0.001874    0.127736    0.051155
    297    2.000000    0.001098    0.000305    0.000792   -0.005373
    318    2.000000    0.001053    0.000326    0.000727   -0.005720
    341    2.000000    0.000781    0.000200    0.000581   -0.006708
    352    2.000000    0.002291    0.000731    0.001560   -0.015705
    360    2.000000    0.001101    0.000298    0.000803   -0.009297
    414    2.000000    0.000635    0.000214    0.000421   -0.005229
    427    2.000000    0.000739    0.000219    0.000520   -0.005731
-------------------------------------------------------------------
Sum:     874.000000    0.153805    0.017732    0.136074   -0.037831
Note: The "Sum" includes all terms including those not printed above

CDA的总结果是CO向NaCl转移了d=0.154个电子，NaCl向CO反馈了b=0.018个电子，NaCl净获得d-b=0.136个电子。b项基本可以忽略不计，而d的主要贡献是由275号复合物轨道造成的。这个复合物轨道的贡献显然值得进一步探究，看看是由CO和NaCl板的哪些片段轨道混合所造成的。为此，输入
5  //将特定复合物轨道对CDA项的贡献进行分解
275  //要考察的复合物轨道序号
0.005  //输出阈值
现在看到如下结果

Occupation number of orbital   275 of the complex:  2.00000000
FragA Orb(Occ.)  FragB Orb(Occ.)      d           b        d - b          r
   5( 2.0000)     271( 2.0000)    0.000000    0.000000    0.000000    0.005183
   5( 2.0000)     289( 2.0000)    0.000000    0.000000    0.000000    0.005642
   5( 2.0000)     336( 2.0000)    0.000000    0.000000    0.000000    0.005127
   5( 2.0000)     341( 2.0000)    0.000000    0.000000    0.000000    0.008960
   5( 2.0000)     435( 0.0000)    0.009358    0.000000    0.009358    0.000000
   5( 2.0000)     439( 0.0000)    0.011934    0.000000    0.011934    0.000000
   5( 2.0000)     443( 0.0000)    0.010929    0.000000    0.010929    0.000000
   5( 2.0000)     447( 0.0000)    0.008032    0.000000    0.008032    0.000000
   5( 2.0000)     453( 0.0000)    0.009177    0.000000    0.009177    0.000000
   5( 2.0000)     465( 0.0000)    0.010755    0.000000    0.010755    0.000000
   5( 2.0000)     469( 0.0000)    0.013735    0.000000    0.013735    0.000000
   5( 2.0000)     473( 0.0000)    0.011282    0.000000    0.011282    0.000000
Sum of above terms:               0.085201    0.000000    0.085201    0.024911

可以看到CO向NaCl转移电子靠的都是CO的5号占据轨道和NaCl空轨道相互作用。NaCl板接收CO转移来的电子分散在其大量空轨道上，比如NaCl板的443号空轨道接受了0.011个电子。在当前的输出阈值下，上面输出的这些项的所有的d项加和为0.085，明显和275号复合物轨道的总d值0.153有很大差距，说明CO还向上面没输出的很多其它的NaCl板的空轨道转移了电子。

下图把275号复合物轨道、CO的5号占据轨道，以及NaCl接收CO贡献来的电子最多的两个轨道（439、469）展示在了一起，便于大家直观了解情况，等值面数值用的是0.02。复合物和片段的轨道可以用Multiwfn载入相应的波函数文件后用主功能0按照《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）说的观看。从此图可以清晰直观地认识到CO的孤对电子往吸附CO的那个Na位点的空轨道上转移了电子。还有大量也具有这种特征的NaCl板的空轨道也接收了CO的孤对电子，这里没全都画出来。

下面再看一下上图中275号复合物轨道是怎么由片段轨道构成的。按0返回到CDA菜单，然后进入选项2，输入275，看到以下信息

 Note: Only the fragment orbitals with contribution >  1.0 % will be shown below, the threshold can be changed by "compthresCDA" in settings.ini

 Occupation number of orbital   275 of the complex:  2.00000000
 Orbital     5 of fragment  1, Occ: 2.00000    Contribution:   87.42 %
 Sum of values shown above:     87.42 %

这告诉你275号复合物轨道由CO的5号占据轨道贡献了87.4%。默认情况下只有贡献大于1%的片段轨道才会输出，虽然NaCl板的一大堆非占据轨道都通过与CO的5号轨道混合而参与了这个复合物轨道的构成，但贡献都小于输出阈值，所以上面没看到。从当前信息可以认识到275号复合物轨道主要体现CO的5号轨道特征，这是为什么上面图中275号复合物轨道和CO的5号轨道看起来十分相似。

4 总结

此文通过一个简单的表面吸附的例子演示了怎么用Multiwfn的CDA功能分析周期性体系的电荷转移本质。讨论电子转移常见的套路就是绘制电子密度差、算算片段电荷，而从本文的例子看到Multiwfn可以从片段轨道相互作用对电子转移的本质进行深入的考察，掌握了本文介绍的方法明显可以给分析增光添彩。

使用Multiwfn做本文介绍的分析请记得在发表文章时引用Multiwfn程序启动时提示的Multiwfn原文，以及引用介绍CDA/GCDA方法的文章（Multiwfn用的实质上是我对CDA广义化后的GCDA形式），这在进入Multiwfn的CDA功能时屏幕上明确提示了。

文/Sobereva@北京科音  2024-Jun-13

1 前言

笔者经常看到有人问“请问如何使用高斯计算电子云密度”之类的问题。之前我在《Multiwfn FAQ》（//www.umsyar.com/452）的Q45里专门回答过这种问题，由于这个话题较大，现在觉得值得再多说说，所以写成个单独的短文。

首先要明确，电子密度（electron density）或者所谓的电子云密度，是一个三维空间函数，没有唯一的考察和分析方法。有人问“怎么计算原子的电子云密度”这种问题显得莫名其妙，波函数分析领域有不同方法如Hirshfeld、Hirshfeld-I、MBIS、Voronoi、AIM等都可以把化学体系整个三维空间划分为独立的原子空间，每个原子的空间里每个位置都有相应的电子密度值，你到底想计算哪里的？再比如“请问如何使用高斯计算电子云密度”这种提问也十分不明不白，本来Gaussian就不是专门干这个的，而且你想以什么形式展现？计算目的何在？明显基本情况都没说清楚，显然只有对电子密度的基本概念完全稀里糊涂的人才会这么问。

2 获得电子密度的方法

在说怎么考察电子密度之前先简谈一下电子密度数据能通过什么方式获得：
(1)实验方法：通过电子衍射、X光衍射、中子衍射，都可以测定出电子密度在三维空间中的分布，从而给出格点数据，描述三维空间里均匀分布的各个格点位置的电子密度值。在测定的电子密度解析度足够高的情况下还能进一步确定原子坐标，看《实验测定分子结构的方法以及将实验结构用于 计算需要注意的问题》（//www.umsyar.com/569）。
(2)理论计算方法：对孤立体系用 程序（如Gaussian），以及对周期性体系用第一性原理程序（如CP2K），都可以计算电子的波函数，并进而得到电子密度。理论方法、基组等因素都直接影响得到的电子密度质量。理论方法对基态的电子密度分布描述的精度对比可参考J. Chem. Theory Comput., 13, 4753 (2017)。在常用的基组范畴中对于价层电子的描述质量，def2-SVP或6-31G**算是可以接受的底限，6-311G(2d,p)或def-TZVP算是不错，def2-TZVPP就算很好了。对于CP2K常见的基组来说，DZVP-MOLOPT-SR-GTH是底限，TZVP-MOLOPT-GTH算是不错，TZV2PX-MOLOPT-GTH就算非常好了。

除上述外，还有特别简单的构造电子密度分布的方法，称为准分子密度（promolecular density）。这就是把各个原子在孤立状态下的电子密度分布根据当前体系的核坐标进行叠加得到的，这样的电子密度相当于没有考虑由于原子间相互作用而造成的电子转移和极化情况的“零阶”的电子密度。这种电子密度不在下文讨论范畴中。

或第一性原理程序在计算后，往往可以导出电子密度的格点数据。例如CP2K可导出电子密度的cube文件，记录了均匀分布的格点上的电子密度，cube文件的介绍见《Gaussian型cube文件简介及读、写方法和简单应用》（//www.umsyar.com/125），它可以直接用于VMD、VESTA等程序观看电子密度等值面等目的，也被一些程序支持用于定量分析。很多程序也可以给出记录波函数信息的文件，如Gaussian可以给出fch/wfn/wfx文件，ORCA和CP2K可以给出molden文件，都可以被Multiwfn等波函数分析程序用来做电子密度的相关分析，参考《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）和《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）。

3 考察电子密度的方法

波函数分析领域博大精深，电子密度的分析仅仅属于波函数分析的一小部分而已，下面会按类别把分析、表征电子密度的方法做简要的提及，这些分析都是Multiwfn完美支持的，不仅可以用于孤立体系也可以用于周期性体系的分析。Multiwfn（//www.umsyar.com/multiwfn）是波函数分析领域功能最丰富、最易用、最流行的程序之一，不了解者建议阅读《Multiwfn入门tips》（//www.umsyar.com/167）和《Multiwfn FAQ》（//www.umsyar.com/452）。

(1)计算某个点的电子密度。用Multiwfn载入波函数文件后在主功能1中输入要考察的坐标即可得到这个位置上电子密度在内的一大堆属性。这有一些实际用处，比如考察一些芳环体系的各个原子在分子平面上方1埃位置的电子密度可以用来考察亲电反应的优先位点，看比如《亲电取代反应中活性位点预测方法的比较》（http://www.whxb.pku.edu.cn/CN/abstract/abstract28694.shtml）中的应用。

(2)图形化考察电子密度在三维空间的分布。这主要用于定性考察电子密度分布情况，包括绘制电子密度在某条路径上的曲线图、特定平面上的图（填色图、等值线图、地形图、梯度线图、向量场图等）、等值面图（适合展现立体分布情况），可以分别用Multiwfn的主功能3、4、5实现，分别把Multiwfn手册4.3、4.4、4.5节的例子看过一遍就懂了。

注意一般来说直接对电子密度作图没什么意义，因为电子密度分布是很“单调”和“乏味”的，即一般在原子核位置是个电子密度极大的峰，向四周以指数型下降，在这样的图上很难观察到化学上感兴趣的信息。但如果你考察的是价层电子密度，则可以讨论很多化学问题，十分建议参看笔者发表的《Revealing Molecular Electronic Structure via Analysis of Valence Electron Density》（http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB201709252），里面有大量生动的例子。价层密度的绘制参考Multiwfn手册里4.6.2节的例子。

值得一提的是，电子密度的0.001 a.u.等值面有特殊的意义，通常被视为气相下分子的范德华表面，里面包围了体系大约98%的电子。电子密度的0.0015 a.u.等值面还可以用来定义动力学直径，见《使用Multiwfn计算分子的动力学直径》（//www.umsyar.com/503）。

Multiwfn还能只绘制特定轨道贡献的电子密度，看《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）。

(3)计算原子电荷（atomic charge）。原子电荷相当于原子的核电荷数减去原子带的电子数，这可以非常清晰、定量地展现实际化学环境中各个原子带的净电荷量。这种分析在文献中被普遍使用，非常有实用价值，比如《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（//www.umsyar.com/696）一文通过原子电荷考察了B9N9和B6N6C6的原子带电情况。我强烈建议阅读《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！》（//www.umsyar.com/714）和《原子电荷计算方法的对比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）以全面了解原子电荷的各种相关知识和计算方法的特点及差异。可以通过Multiwfn的主功能7的各种选项以ADCH、Hirshfeld-I、RESP、EEM、AIM等许多不同方法来计算原子电荷，把手册4.7节的例子好好看一遍就能很快掌握了，对于周期性体系的原子电荷的计算我还专门写了文章《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）。

将某个片段中的所有原子的电荷加和就得到了片段电荷（fragment charge），可以确切了解片段的带电情况、分析电子转移量。在Multiwfn中计算片段电荷很方便，在主功能7里选择选项-1，输入片段里原子的序号来定义片段，之后再照常选择一种方法计算原子电荷，则原子电荷输出后就会输出片段电荷。之前还录过演示视频《使用Multiwfn计算分子的某个片段的电荷》（https://www.bilibili.com/video/av26312703/）。

还有个概念叫布居分析（population analysis），所谓的布居就是指某个对象带的电子数。除了计算原子带的电子数外，在Multiwfn的主功能7中基于波函数信息还可以用Mulliken、SCPA等方法给出各个基函数、壳层、原子轨道上带的电子数，参考比如《利用布居分析判断基函数与原子轨道的对应关系》（//www.umsyar.com/418）。非常灵活的Multiwfn还可以给出各个原子上带的特定类型电子（如sigma、pi电子）的量，参考《使用Multiwfn基于Hirshfeld-I划分计算特定类型电子在各个原子上的分布量》（//www.umsyar.com/697）。

(4)电子密度差（有个很难听、我很讨厌的别称叫差分电荷密度）。这是在不同状态或不同体系间将电子密度相互求差得到的，由此可以展现出单个体系电子密度分布无法展现出的信息，在定量分析电子转移、体系对外场的响应、考察共价键的形成情况等方面有重要价值。电子密度差有多种类型，见《使用Multiwfn作电子密度差图》（//www.umsyar.com/113）中的说明，里面还有在Multiwfn中的绘制过程的介绍。电子密度差在周期性体系的研究中也很重要，而且对一些体系将密度差转化为电荷位移曲线后在分析讨论时可以更为清晰准确，看《使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线》（//www.umsyar.com/638）。

预测反应位点常用的福井函数和双描述符本质上也是电子密度差，看《概念密度泛函综述和重要文献合集》（http://bbs.keinsci.com/thread-384-1-1.html）和《使用Multiwfn计算双描述符势预测反应位点》（//www.umsyar.com/708）里面提到的相关资料。分析（超）极化率本质常用的（超）极化率密度也是基于密度差所定义的，看《使用Multiwfn计算（超）极化率密度》（//www.umsyar.com/305）。

(5)一些与电子密度分布密切相关的量，例如：
• 电偶极矩、电多极矩：体系的电荷分布可以做电多极展开成为单极矩+偶极矩+四极矩+八极矩...的描述形式，其中每个部分都可以分为核电荷与电子密度各自贡献的部分，Multiwfn手册3.18.3节有完整的计算公式和相关说明。偶极矩可以体现正电中心和负电中心间的分离程度，四极矩可以体现电荷分布偏离球对称的程度。在Multiwfn中载入波函数文件后，选择主功能300的子功能5即可计算它们。Multiwfn还可以计算各个原子和片段的偶极矩和多极矩，看《使用Multiwfn计算分子片段的偶极矩和复合物中单体的偶极矩》（//www.umsyar.com/558）。

• <r^2>：这衡量电子分布的延展程度，见《电子空间范围<r^2>和电子径向分布函数的含义以及在Multiwfn中的计算》（//www.umsyar.com/616）。

• 电子密度的径向分布函数：这可以考察相对于某个点不同距离的单位厚度壳层内的电子数目的相对差异，在//www.umsyar.com/616里也有介绍。

(6)电子密度的拓扑分析。任何三维函数都可以在Multiwfn的主功能2中做拓扑分析，从而得到它们的各种临界点（梯度为0的点）的坐标，以及相应位置的各种实空间函数的值。做电子密度的拓扑分析是Bader提出的非常知名的Atom-in-molecules (AIM)分析中用到的关键分析手段，此分析得到的键临界点（BCP）通常是原子间相互作用路径上最有代表性的一个点，通过这个位置的各种函数值可以讨论很多问题。例如《透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍》（//www.umsyar.com/513）指出利用BCP处的电子密度可以很好地预测氢键强度，Multiwfn手册4.2.1节和《计算分子内氢键键能的几种方法》（//www.umsyar.com/522）都有具体例子。再比如，AIM分析得到的环临界点上的电子密度和这个环的芳香性有密切联系，看《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）中的相关部分。关于更多的AIM和临界点的相关知识参考《AIM学习资料和重要文献合集》（http://bbs.keinsci.com/thread-362-1-1.html）和《使用Multiwfn做拓扑分析以及计算孤对电子角度》（//www.umsyar.com/108）。

(7)电子密度的盆分析。盆分析是指用特定实空间函数定义的零通量面将三维空间划分为一个个独立的子区域，这被成为盆，然后在盆中对特定函数进行积分，就可以考察很多问题。对电子密度就可以进行盆分析，得到的叫原子盆，在里面积分电子密度对应的就是AIM方法定义的原子带的电子数，用Multiwfn的主功能17可以实现，看《使用Multiwfn做电子密度、ELF、静电势、密度差等函数的盆分析》（//www.umsyar.com/179）和手册4.17.1节。对价层电子密度或ELF函数进行盆分析，然后在盆里积分电子密度，能够获得成键区域电子数、孤对电子数等信息，看《Revealing Molecular Electronic Structure via Analysis of Valence Electron Density》（http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB201709252）中的讨论。Multiwfn甚至还可以对密度差做盆分析然后在盆内积分密度差，这可以定量考察特征区域电子密度的变化量，//www.umsyar.com/179和手册4.17.4节都有例子。孤对电子的ELF盆、电子密度=0.001 a.u.等值面内和ELF=0.5等值面内交集区域内的电子数称为high ELF localization domain population (HELP)，是孤对电子数的一种特殊定义，与电离能、前线轨道能量等不少体系的属性有联系，参看手册4.17.8节。

(8)电子密度的衍生函数。有很多函数衍生自电子密度，它们的分析也是波函数分析的重要组成部分、有巨大实际意义。比如《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用
//www.umsyar.com/621》、《使用IRI方法图形化考察化学体系中的化学键和弱相互作用》（//www.umsyar.com/598）、《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（//www.umsyar.com/667）里介绍的IGMH和IRI分析就是基于电子密度及其导数做的，在分析化学体系中的相互作用方面已经非常流行。再比如电子密度的拉普拉斯函数及其椭率对于考察化学键的特征很有意义，参考《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）和《AIM键临界点处电子密度拉普拉斯值符号判断相互作用类型失败原因的图形分析》（//www.umsyar.com/161）。静电势也是基于电子密度定义的，重要性极高，我写过大量相关文章，汇总看《静电势与平均局部离子化能相关资料合集》（http://bbs.keinsci.com/thread-219-1-1.html）。信息论在化学体系中的应用已被广泛探索，其中有很多函数和量是基于电子密度定义，参看《使用Multiwfn计算各种与信息论相关的量（information-theoretic quantities）》（//www.umsyar.com/537）。

4 总结

上文对电子密度相关分析做了一个较为全面的说明。从本文可见，电子密度相关的分析超级丰富，所以提问时切勿简简单单问别人“怎么分析/计算/考察电子密度？”，必须说清楚你想怎么分析讨论。如果你对方法学不懂，就要告诉别人你的分析目的是什么，这样别人才能告诉你具体怎么做，而且如果对于你的目的不适合光靠分析电子密度来说事，别人也能及早指出。从上文也看到，电子密度相关分析主要都是Multiwfn等波函数分析程序的事，所以提问时也别问诸如“怎么用Gaussian计算电子密度？”这种问题。

限于本文篇幅，上面的各种分析方法只能是三言两语简单提及，所有这些方法在“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/workshop/WFN_content.html）里有最全面系统的讲解并给了大量例子，可以一次性完整学透，欢迎参加。

文/Sobereva@北京科音   2024-Jun-8

1 前言

原子电荷（atomic charge）即原子所带的净电荷，对于讨论化学体系中原子的带电状态有重要意义，作为原子的一种描述符它也和原子的很多属性有密切联系。强烈建议阅读《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！》（//www.umsyar.com/714）里提到的笔者写的原子电荷的综述以全面了解原子电荷。

强大的波函数分析程序Multiwfn（//www.umsyar.com/multiwfn）支持非常丰富的原子电荷计算方法，其中多数已经支持周期性体系。Multiwfn支持的方法中有一类是基于模糊式原子空间划分，也就是每个原子对应一个平滑的权重函数，权重函数与体系的电子密度的乘积的全空间积分就是这个原子带的电子数（原子的布居数），原子的核电荷数减去它就是原子电荷。Multiwfn支持的这类方法包括Hirshfeld、Hirshfeld-I、MBIS、Becke，以及对Hirshfeld电荷做后校正得到的ADCH和CM5电荷。从2024-May-25更新的Multiwfn版本开始，Hirshfeld、CM5、Hirshfeld-I、MBIS都已经支持了周期性体系，本文将具体介绍怎么用Multiwfn与非常流行、高效且免费的第一性原理程序CP2K相结合非常方便快速地计算这些原子电荷。Multiwfn计算这些电荷的功能是普适的，并不限于结合CP2K，也可以基于其它第一性原理程序产生的体系的电子密度的cube文件来算，还支持VASP产生的记录电子密度的CHGCAR文件。
注：Multiwfn更老的一些版本也支持基于周期性波函数计算Hirshfeld和CM5电荷，但那时候用的算法对于周期性体系效率极低，本文介绍的做法与之有天壤之别。

下文第2节简要介绍一下Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷的基本特点，第3节介绍用Multiwfn+CP2K计算它们的具体方法，第4节给出具体例子，同时还会介绍怎么计算片段电荷。使用Multiwfn计算原子电荷在发表文章时记得需要按Multiwfn启动时的提示恰当引用Multiwfn。

Multiwfn可以在官网//www.umsyar.com/multiwfn免费下载，如果对Multiwfn不了解的话建议看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。如果你用CP2K创建给Multiwfn算原子电荷用的输入文件的话，笔者假定你对CP2K的使用已经有基本且正确的了解，不具备这些知识的话非常建议通过北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）完整、系统地学习一遍。本文的例子利用到Multiwfn创建CP2K的输入文件，相关介绍见《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）。

2 Hirshfeld、Hirshfeld-I、CM5和MBIS等原子电荷简介

Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷在Multiwfn手册3.9节的相应小节里都有很详细的介绍，因此这里不介绍细节，只是简单说一下它们的特点。

Theor. Chim. Acta (Berl.), 44, 129 (1977)中提出的Hirshfeld电荷是被广为使用的原子电荷，计算方式简单，物理意义清楚，在 研究孤立体系和第一性原理研究周期性体系的领域中都用得很多。Hirshfeld电荷为人诟病的地方是原子电荷数值整体严重偏小，对偶极矩、静电势重现性很差，这在我写的《原子电荷计算方法的对比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）里的对比测试中有充分体现。

JCP, 126, 144111 (2007)提出的Hirshfeld-I是在Hirshfeld基础上的改进，通过迭代方式不断更新原子权重函数直到收敛，这使得原子权重函数所描述的原子空间可以响应实际化学环境。其原子电荷数值显著大于Hirshfeld，对静电势的重现性也有明显改进。代价是Hirshfeld-I需要做迭代，往往几十轮，计算耗时明显高于Hirshfeld，而且需要事先提供体系中各种元素的各个氧化态的原子径向密度信息，实现起来较麻烦（为了实现Hirshfeld-I，我计算了360多个原子的径向密度文件并在Multiwfn中自带）。

一定要注意，虽然CP2K程序自己也支持Hirshfeld-I电荷的计算，但我测试发现其结果完全不靠谱，我检查了其源代码也确认了其实现方式明显不符合Hirshfeld-I的标准定义，因此切勿用CP2K自身的Hirshfeld-I电荷计算功能。CP2K算的Hirshfeld电荷在结合SHAPE_FUNCTION DENSITY选项时的结果不那么离谱，但结果和Multiwfn算的也存在不可忽视的差异，应当以Multiwfn的结果为准。

JCTC, 8, 527 (2012)中提出的CM5电荷和我在J. Theor. Comput. Chem., 11, 163 (2012)中提出的ADCH电荷一样都是对Hirshfeld电荷的后校正，对大部分情况都能使得它的原子电荷变得更大、更符合化学直觉，对静电势的重现性也明显变得更好。ADCH没有CM5那么强的经验性、不牵扯一堆经验参数，因而原理更为理想、更推荐使用。CM5最初是用于分子体系的，如今在应用于晶体方面也已经有了一定探索，例如JCTC, 16, 5884 (2020)将CM5与其它一些原子电荷对比考察了它们对一些固体的电荷分布的描述。CM5还被用于一些计算化学领域的方法，比如《比SMD算溶解自由能更好的溶剂模型uESE的使用》（//www.umsyar.com/593）里介绍的uESE溶剂模型依赖于CM5电荷，CM5乘以1.2后适合结合OPLS-AA力场做动力学模拟，见《计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒人脚本》（//www.umsyar.com/585）。

ADCH电荷在分子体系的研究方面已被广为使用，之前有人问我是否会把ADCH电荷扩展到周期性体系。我目前没这个打算，因为ADCH电荷的思想对周期性体系往往不适用。ADCH方法将原子偶极矩展开为周围原子的校正电荷，但对于诸如NaCl这样的体系，每个原子所处的环境是中心对称的，因此原子偶极矩为0、校正电荷都为0，故ADCH电荷会与Hirshfeld电荷完全相同。

MBIS (Minimal Basis Iterative Stockholder)电荷于JCTC, 12, 3894 (2016)中提出。类似于Hirshfeld-I，此方法也是通过迭代方式（因此耗时较高）优化原子空间直到收敛。MBIS比Hirshfeld-I的一大好处是不需要利用各个元素各个氧化态的原子径向密度信息，因此实现起来简单得多，也更优雅。原文里的测试体现出MBIS比Hirshfeld-I整体更具有优势。但MBIS经常收敛较慢，对于某些体系甚至需要200多轮才能充分收敛。目前Multiwfn里MBIS最高支持到Rn元素。

根据我的测试，对于大多数固体体系，原子电荷大小的关系是MBIS和Hirshfeld-I最大、谁大不一定，CM5大小介中，Hirshfeld电荷最小。如果你对数值的绝对大小不那么在乎，只是想体现相对大小关系的话，最“朴素”且很便宜的Hirshfeld也可以用。但如果对定量数值关心的话，建议用其它的。CM5、Hirshfeld-I和MBIS用于晶体哪个最理想没有定论，都可以算一算看看，如果倾向于得到较大数值的结果可以优先考虑Hirshfeld-I和MBIS。如果其中有的原子电荷明显不合理那就弃掉换别的方法，比如O的电荷如果算出来-2.34那就不应该用，因为显然O最多只能带-2电荷。

特别要强调的是，绝对不要把原子电荷与氧化态搞混，也切勿用原子电荷判断氧化态！氧化态怎么正确地判断看《使用Multiwfn结合CP2K计算晶体中原子的氧化态》（//www.umsyar.com/711）和《使用Multiwfn通过LOBA方法计算氧化态》（//www.umsyar.com/362）。

还有很多其它原子电荷计算方法，诸如Mulliken电荷、Lowdin电荷、SCPA电荷、AIM电荷（被个别文章称为Bader电荷）等，Multiwfn也都支持将它们用于周期性体系，这不属于本文的范畴。Mulliken电荷结合CP2K里常用的MOLOPT系列和pob系列基组用于周期性体系也说得过去。虽然它是基于希尔伯特空间划分的原子电荷中最原始、思想最朴素的一种，但实际结果一般还算定性正确，而且计算耗时极低。被用得很多的AIM电荷实际上很糟糕，往往出现违背化学常识的结果，这在前述的《原子电荷计算方法的对比》里已经充分体现了，不建议用。

3 Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷在Multiwfn中的计算方法

3.1 CP2K用户的情况

首先说对于CP2K用户怎么计算。Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷都是完全基于电子密度计算的，不直接牵扯到波函数，在Multiwfn中计算它们有两种可以用的输入文件：
(1)用CP2K产生的体系的电子密度cube文件作为输入文件。这种情况在CP2K计算时可以考虑k点，是我最推荐的，耗时非常低。如果不了解cube文件的格式的话，参看《Gaussian型cube文件简介及读、写方法和简单应用》（//www.umsyar.com/125）。注意CP2K有bug（起码直到我撰文时最新的CP2K 2024.1为止），产生的cube文件里的原子有效电荷都为0，这会导致Multiwfn无法正确计算原子电荷。这里说的原子的有效核电荷是指当前基组描述的原子的电子数，对于赝势基组它对应价电子数。为此，要么用文本编辑器手动改cube文件中每个原子的有效核电荷，要么把cube文件第一行写成各个元素的有效核电荷，例如N 5 B 3 Ti 12，这代表当前体系中的N、B、Ti的有效核电荷数分别为5、3、12
(2)用CP2K产生的记录了体系波函数的molden文件作为输入文件。详见《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）。这个方法的缺点在于没法考虑k点，对于小晶胞体系必须先扩胞到足够大，从而在计算时只考虑gamma点，这显然使得在CP2K中和Multiwfn中的计算耗时都会比用(1)的方式高很多

用Multiwfn载入上述输入文件其中一种后，进入主功能7，然后选择相应选项即可计算Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷。对于molden文件当输入文件时，Multiwfn还会让你输入计算用的格点间距，格点间距越小耗时越高而计算精度越高，通常用0.2 Bohr就可以，是精度和耗时的较好权衡。而使用cube文件做为输入文件时，计算原子电荷用的格点间距与此文件的格点间距直接对应，CP2K的CUTOFF设得越大则产生的cube文件里的格点间距会越小。

由于Multiwfn计算上述电荷用的是均匀分布的格点，不可能准确积分较重原子的内核区域电子密度，因此上述输入文件应当是CP2K做GPW计算得到的，也需要用赝势基组，此时cube文件记录的电子密度或者molden文件里记录的波函数只对应于价电子。

顺带一提，如果你要计算Mulliken电荷，用Multiwfn载入上述的molden文件，然后进主功能7，依次选择5、1即可得迅速到结果。Multiwfn中计算AIM电荷的方法笔者以后会另文说明。

3.2 其它程序用户的情况

如果你是VASP用户，可以用记录晶胞中价电子密度的CHGCAR文件给Multiwfn用于计算前述原子电荷，等同于3.1节说的第2类输入文件。只要文件名里包含CHGCAR字样就会被Multiwfn当做是CHGCAR文件来载入。注意为了让Multiwfn得知各类原子的有效核电荷数，CHGCAR文件的第一行需要以Nval为开头，后面依次写上各类元素的名字和有效核电荷，并以空格分隔。比如写为Nval Na 1 Cl 7，这代表Na和Cl的有效核电荷数分别为1和7。

如果你是其它第一性原理程序的用户，若程序可以产生记录晶胞中价电子密度的cube文件（即等同于3.1节说的第2类输入文件），也都可以用于给Multiwfn计算前述的原子电荷。

之后在Multiwfn中的计算过程和CP2K用户没任何区别。

4 实例

下面给出Multiwfn结合CP2K计算Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷的具体例子。相关的文件都可以在//www.umsyar.com/attach/712/file.zip中得到（4.2、4.3和4.5节的cube/molden文件除外，因为文件过大）。

4.1 TiO2晶体

此例对TiO2的金红石晶型计算前述的原子电荷。它的实验晶体结构的cif文件是本文文件包里的TiO2-Rutile.cif，结构图如下所示

先用CP2K计算出这个晶胞的电子密度的cube文件，用于之后Multiwfn做原子电荷计算。首先用Multiwfn创建CP2K输入文件，启动Multiwfn并载入TiO2-Rutile.cif，然后输入
cp2k  //创建CP2K的输入文件
[回车]  //产生的文件名用默认的TiO2-Rutile.inp
-3  //产生cube文件
1  //电子密度
8  //设置k点
4,4,6  //与各个晶胞边长相乘皆为18埃左右，够用了
0  //产生输入文件，计算级别为默认的PBE/DZVP-MOLOPT-SR-GTH

用CP2K运行刚得到的TiO2-Rutile.inp，马上就算完了，得到了TiO2-Rutile-ELECTRON_DENSITY-1_0.cube文件，在本文的文件包里已经提供了。用文本编辑器打开此文件，把第一行改为Ti 12 O 6并保存。这代表Ti和O的有效核电荷分别为12和6，直接对应于输入文件里它们分别用的DZVP-MOLOPT-SR-GTH-q12和DZVP-MOLOPT-SR-GTH-q6基组名里面的q后面的值。

启动Multiwfn，输入
TiO2-Rutile-ELECTRON_DENSITY-1_0.cube  //写此文件的实际路径
y  //从第一行载入各个元素有效核电荷
7  //布居分析与原子电荷计算
1  //Hirshfeld电荷

立刻看到如下Hirshfeld电荷计算结果。由于格点间距不是无穷小，因此等价的原子算出来的原子电荷也会有轻微的数值差异，可以自己手动取平均。如果用更大的CUTOFF，则cube文件的格点间距会更小，不等价性会更低。
Final atomic charges:
Atom    1(Ti):     0.59349459
Atom    2(O ):    -0.29461903
Atom    3(O ):    -0.29472601
Atom    4(Ti):     0.58363165
Atom    5(O ):    -0.29385236
Atom    6(O ):    -0.29385236
现在Multiwfn还问你是否把原子电荷导出为chg文件，这里选n不导出。

接下来计算CM5电荷。在Multiwfn主功能7里选择16，马上就得到了结果
                      ---------- CM5 charges ----------
Atom:    1Ti  CM5 charge:    1.348393  Hirshfeld charge:    0.593495
Atom:    2O   CM5 charge:   -0.672093  Hirshfeld charge:   -0.294619
Atom:    3O   CM5 charge:   -0.672191  Hirshfeld charge:   -0.294726
Atom:    4Ti  CM5 charge:    1.338580  Hirshfeld charge:    0.583632
Atom:    5O   CM5 charge:   -0.671306  Hirshfeld charge:   -0.293852
Atom:    6O   CM5 charge:   -0.671306  Hirshfeld charge:   -0.293852
Summing up all CM5 charges:     0.00007649
可见CM5电荷的大小明显大于Hirshfeld电荷，看起来更合理，而Hirshfeld电荷值明显偏小。

下面再演示计算Hirshfeld-I电荷。把Multiwfn自带的examples目录下的atmrad目录挪到当前目录下，这样Multiwfn计算Hirshfeld-I电荷时就会自动用此目录下的各种元素各个价态的径向电子密度。之后在主功能7里选15，再选1用默认的设置进行Hirshfeld-I电荷计算，经过16轮迭代得到了结果：
Final atomic charges:
Atom    1(Ti):     3.01789464
Atom    2(O ):    -1.50888824
Atom    3(O ):    -1.50886131
Atom    4(Ti):     3.01787620
Atom    5(O ):    -1.50897240
Atom    6(O ):    -1.50897240
可见Hirshfeld-I电荷的大小不仅显著大于Hirshfeld电荷，比CM5电荷还要大很多。虽然数值很大，但完全在合理范围内，毕竟都没有大过原子的氧化态。

接下来再计算MBIS电荷。在主功能7里选择20，再选1用默认的设置进行MBIS电荷计算，经过19轮迭代得到了结果：
Final atomic charges:
Atom    1(Ti):     1.73473636
Atom    2(O ):    -0.86733169
Atom    3(O ):    -0.86731253
Atom    4(Ti):     1.73467041
Atom    5(O ):    -0.86734303
Atom    6(O ):    -0.86734303
对于当前例子，MBIS的数值介于CM5和Hirshfeld-I之间。但也有很多反例，MBIS也经常比Hirshfeld-I电荷还大不少，后文有例子。

以上原子电荷的计算顺序是随意的，不要误以为必须按上述顺序计算，你想算哪个就算哪个即可。

4.2 AlN晶体

这一节计算AlN晶体的原子电荷。完全可以按照和上一节相同的做法基于电子密度的cube文件来计算，但本节演示一下基于超胞的molden文件的计算过程。AlN原胞的cif文件是本文文件包里的AlN.cif，将之载入Multiwfn，然后输入以下命令创建CP2K输入文件
cp2k  //创建CP2K输入文件
AlN_443.inp  //产生的输入文件名
-11  //进入几何操作界面
19  //构造超胞
4  //a方向扩胞成原本的4倍
4  //b方向扩胞成原本的4倍
3  //c方向扩胞成原本的3倍
-10  //返回
-2  //要求产生molden文件
0  //产生输入文件

用CP2K运行刚产生的AlN_443.inp，得到AlN_443-MOS-1_0.molden。在此文件开头插入以下内容定义晶胞和各个元素的价电子数
[Cell]
12.44400000     0.00000000     0.00000000
-6.22200000    10.77682012     0.00000000
 0.00000000     0.00000000    14.93400000
[Nval]
Al 3
N 5

启动Multiwfn，载入AlN_443-MOS-1_0.molden，然后输入
7  //布居分析与原子电荷计算
16  //CM5电荷
[回车]  //使用均匀分布的格点计算
[回车]  //使用默认的0.2 Bohr格点间距。格点间距越小结果越精确但计算越耗时
在i9-13980HX CPU上24核并行，50秒就算完了，结果为

                      ---------- CM5 charges ----------
Atom:    1Al  CM5 charge:    0.467945  Hirshfeld charge:    0.302813
Atom:    2N   CM5 charge:   -0.467945  Hirshfeld charge:   -0.302812
Atom:    3Al  CM5 charge:    0.467945  Hirshfeld charge:    0.302812
Atom:    4N   CM5 charge:   -0.467946  Hirshfeld charge:   -0.302814
...略

Al的电负性很小，N的电负性很大，因此二者的原子电荷理应相差悬殊。当前的Hirshfeld电荷明显偏小。CM5电荷的数值虽然更大，但还是显得有些偏低。

下面再计算Hirshfeld-I电荷。在主功能7里输入
15  //Hirshfeld-I
1  //开始计算
[回车]  //使用默认的0.2 Bohr格点间距

经过37轮收敛，结果如下。可见Hirshfeld-I电荷的数量级比较合理，体现出AlN中成键的极强离子性
 Atom    1(Al):     1.80317384
 Atom    2(N ):    -1.80317397
 Atom    3(Al):     1.80317257
 Atom    4(N ):    -1.80317539
...略

最后再计算下MBIS电荷。在主功能7里输入
20  //MBIS
1  //开始计算
[回车]  //使用默认的0.2 Bohr格点间距

对此体系MBIS收敛很慢，经过190轮才达到收敛，在i9-13980HX上24核并行算了5分钟，在双路7R32 96核并行下算了不到一分钟。结果如下。可见对此体系，MBIS电荷的数值很大，和氧化态几乎正好是一样的。所以，如前所述，若你就是想算出来比较大的原子电荷，可以先考虑MBIS，但也有一些例外，后面会提到。
Atom    1(Al):     2.99599978
Atom    2(N ):    -2.99601488
Atom    3(Al):     2.99600129
Atom    4(N ):    -2.99599935

4.3 MOF-5晶体

MOF-5晶体含有Zn、O、C、H。其cif文件是本文文件包里的MOF-5.cif。由于此体系的晶胞边长约26埃，因此计算电子密度格点数据时不用考虑k点。按照3.1节的做法，基于PBE/DZVP-MOLOPT-SR-GTH级别的电子密度cube文件（第一行改为Zn 12 O 6 C 4 H 1），用Multiwfn计算各种原子电荷，结果如下。此体系中O有两类，分别给出，C则处于较多不同的化学环境，其电荷就不给出了。

Hirshfeld：
Zn=0.372 O=-0.332/-0.200 H=0.056

CM5：
Zn=0.644 O=-0.584/-0.305 H=0.116

Hirshfeld-I（36轮收敛，2*7R32机子上96核并行4分多钟算完）
Zn=1.568 O=-1.917/-0.746 H=0.139

MBIS（39轮收敛，2*7R32机子上96核并行3分多钟算完）
Zn=1.546 O=-1.662/-0.782 H=0.190

对这个体系，MBIS和Hirshfeld-I相符得很好，而且Zn的电荷大小很合理。Hirshfeld电荷再次数值最小、CM5比之大一些。

4.4 其它：MoS2、BaTiO3、MgF2、NaCl

注意MBIS电荷并非总是较大。例如《使用Multiwfn结合CP2K计算晶体中原子的氧化态》（//www.umsyar.com/711）里的单层MoS2的例子，在PBE/DZVP-MOLOPT-SR-GTH级别的波函数下，Mo的电荷的计算结果是Hirshfeld=0.330，CM5=0.827，Hirshfeld-I=1.987，MBIS=-0.059，Mulliken=-0.771。对这个体系来说CM5和Hirshfeld-I显得靠谱，MBIS电荷接近0有点说不通，Mulliken电荷明显有误导性。

BaTiO3在PBE/DZVP-MOLOPT-SR-GTH级别下的原子电荷计算结果如下。Hirshfeld还是严重偏小。Hirshfeld-I大得离谱，Ba的原子电荷都超过氧化态+2了。相较之下CM5和MBIS的结果比较靠谱，正好二者算的O的原子电荷差不多，差异主要在Ba和Ti的电荷的相对大小上。
Atom:    1Ba  CM5 charge:    1.103253  Hirshfeld charge:    0.404357
Atom:    2Ti  CM5 charge:    1.200467  Hirshfeld charge:    0.537212
Atom:    3O   CM5 charge:   -0.767963  Hirshfeld charge:   -0.310211
Atom:    4O   CM5 charge:   -0.767892  Hirshfeld charge:   -0.315693
Atom:    5O   CM5 charge:   -0.767892  Hirshfeld charge:   -0.315693

Hirshfeld-I
Atom    1(Ba):     2.37169239
Atom    2(Ti):     3.02878350
Atom    3(O ):    -1.78799416
Atom    4(O ):    -1.80625420
Atom    5(O ):    -1.80625434

MBIS：
Atom    1(Ba):     0.57391985
Atom    2(Ti):     1.69416825
Atom    3(O ):    -0.74678425
Atom    4(O ):    -0.76066530
Atom    5(O ):    -0.76066537

再看两个典型离子化合物的情况，都在PBE/DZVP-MOLOPT-SR-GTH级别下计算。对NaCl计算的Na的原子电荷为：Hirshfeld=0.217，CM5=0.443，Hirshfeld-I：1.053，MBIS：1.043，Mulliken：0.596
对MgF2计算的Mg的原子电荷为：
Hirshfeld=0.414，CM5=0.828，Hirshfeld-I：2.048，MBIS：1.888，Mulliken：1.439
可见对于强离子性化合物，Hirshfeld-I和MBIS电荷都在氧化态附近，并且由于方法定义的缺陷，往往电荷的大小还轻微大于氧化态的大小（但如果用大核赝势，对Na和Mg只描述最外层的3s价电子，则不会有这个问题）。其它原子电荷的大小关系为Mulliken>CM5>Hirshfeld。

以MgF2为例，这里顺带一提CP2K算的Hirshfeld-I电荷明显不对，千万不要用。在CP2K输入文件的&DFT/&PRINT里加上以下内容要求输出Hirshfeld-I电荷后，CP2K直接给出的Mg的Hirshfeld-I电荷为0.522，和Multiwfn算的相距甚远。
&HIRSHFELD
  SHAPE_FUNCTION DENSITY
  SELF_CONSISTENT T
&END HIRSHFELD

4.5 NaCl板吸附CO的片段电荷

NaCl板吸附CO是《使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线》（//www.umsyar.com/638）里的一个例子体系，这里基于这个模型计算一下原子电荷和CO片段的电荷。CP2K优化这个体系后产生的NaCl_CO-1.restart文件在本文的文件包里提供了，将之载入Multiwfn，然后输入
cp2k  //创建CP2K输入文件
NaCl_CO.inp  //产生的输入文件名
-7  //设置周期性
XY
-3  //要求产生cube文件
1  //电子密度
-2  //产生molden文件（用于之后算Mulliken电荷的目的）
0  //产生输入文件

计算完毕后，把NaCl_CO-ELECTRON_DENSITY-1_0.cube文件第一行改为Na 9 Cl 7 C 4 O 6。先用Hirshfeld方法各个原子的电荷和CO的片段电荷。启动Multiwfn，载入此cube文件，然后输入
7  //布居分析与原子电荷计算
-1  //定义片段
109,110  //CO的原子序号
1  //Hirshfeld电荷
结果为
...略
 Atom  109(C ):     0.13622804
 Atom  110(O ):    -0.03514734

Fragment charge:    0.10108070
Fragment population:    9.89891930

即CO的Hirshfeld片段电荷为0.101。

之后再选择CM5方法，结果为
...略
 Atom:  109C   CM5 charge:    0.147306  Hirshfeld charge:    0.136228
 Atom:  110O   CM5 charge:   -0.081972  Hirshfeld charge:   -0.035147
 Summing up all CM5 charges:     0.00510030

 Fragment charge:    0.06533361

之后再选择Hirshfeld-I方法，结果为
 Atom  109(C ):     0.09790054
 Atom  110(O ):    -0.08744136

 Fragment charge:    0.01045918

之后再选择MBIS方法，结果为
...略
 Atom  109(C ):     0.19642171
 Atom  110(O ):    -0.15937966

 Fragment charge:    0.03704205

再来算Mulliken电荷。用Multiwfn载入加入了恰当的[Cell]和[Nval]信息的CP2K计算后产生的NaCl_CO-MOS-1_0.molden文件，输入
7  //布居分析与原子电荷计算
-1  //定义片段
109,110  //CO的原子序号
5  //Mulliken电荷
1  //在屏幕上输出结果
结果为
...略
Atom   109(C )    Population:  3.60671601    Net charge:  0.39328399
Atom   110(O )    Population:  6.23484007    Net charge: -0.23484007
Total net charge:   -0.00000086

Fragment charge:    0.15844392

可见不同方法计算的结果在定量上有一定差异，但都有共性，即CO整体带少量正电，即电子从CO往NaCl板上转移了一些。并且C和O分别带正电和负电，且原子电荷的绝对值都是C大于O。

5 总结

本文简要介绍了四种知名的基于模糊式原子空间划分计算原子电荷的方法的特点，包括Hirshfeld、Hirshfeld-I、CM5和MBIS，并介绍了它们在Multiwfn中的计算方法，并且给了很多具体例子，以使得读者能了解如何使用Multiwfn结合CP2K或其它第一性原理程序很容易地计算它们。同时还说明了怎么直接计算片段电荷，这是定量讨论体系中两部分间电子转移的最严格的方式。推荐将这些原子电荷应用于实际问题的研究当中以考察固体与表面体系的电荷分布情况。记得届时应恰当引用Multiwfn启动时提示的Multiwfn原文以及相应的原子电荷的原文。

通过本文的一些体系的对比可以明显看到，Hirshfeld电荷原理上最简单，但明显缺点是整体严重偏小，CM5通常比它大不少，而且电荷的稳定性不错。虽然结合化学直觉来看，CM5电荷的绝对大小仍然往往能偏小，但至少也不会明显离谱，总比Hirshfeld电荷更建议使用。Hirshfeld-I和MBIS电荷普遍比较大，谁更大不一定，至少几乎总比CM5还大不少，从绝对大小来看往往比CM5更合理，但稳健性比CM5弱一些，对极个别体系可能会表现得不符合常识。对于离子性特征很强的化合物，Hirshfeld-I和MBIS电荷往往接近氧化态，由于方法并不完美，有时候原子电荷甚至比氧化态还要大一点点。当然，没有任何原子电荷计算方法是绝对完美的，本文介绍的这些原子电荷在实际当中都可以用，可以根据它们的特点和实际结果决定选用，在横向对比时方法必须统一。

顺带一提，对于计算固体表面或者多孔物质的原子电荷以用于基于力场的动力学模拟的目的，我目前最推荐CP2K算REPEAT电荷，北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）里讲“电子结构的分析”的部分做了专门的介绍。但REPEAT电荷需要在真空区定义拟合点，因而无法用于诸如本文所示的AlN、TiO2等致密的固体，而且对于远离真空区的原子也没法拟合出质量靠谱的原子电荷。所以REPEAT电荷有特定用处，但不是像本文介绍的那些原子电荷一样属于普适性的原子电荷计算方法。