//www.umsyar.com/category/QC //www.umsyar.com/727 2024-10-28T20:02:00+08:00 全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！ 文/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程序做波函数分析对弱相互作用研究的关键性价值！ //www.umsyar.com/725 2024-08-31T23:40:00+08:00 18个氮原子组成的环状分子长什么样？一篇文章全面揭示18氮环的特征！ What does a ring molecule composed of 18 nitrogen atoms look like? A paper fully reveals the characteristics of cyclo[18]nitrogen! 文/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程序做波函数分析考察新颖物质的电子结构的重要意义。 //www.umsyar.com/714 2024-06-01T18:29:00+08:00 一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！ 文/Sobereva@北京科音   2024-Jun-1 原子电荷（atomic charges）也称分数电荷（partial charges）是波函数分析领域的极其重要的概念，不仅对于了解原子在化学环境中的带电状态有关键意义，与分子力场、反应位点、反应活性、分子间相互作用、静电势以及很多属性（如pKa和化学位移）等都有密切联系，而且很多计算化学领域的方法都利用或者建立在原子电荷基础上。 笔者应刘述斌教授的邀请，撰写了Wiley出版社2024年出版的专著Exploring Chemical Concepts Through Theory and Computation中的第6章，名为Partial Charges。文章访问链接为https://onlinelibrary.wiley.com/doi/10.1002/9783527843435.ch6 （链接：https://pan.baidu.com/s/1y-3dnd5971QyNmeCf7OnGw 提取码：98ya） 在此文中，笔者对原子电荷的概念、理论意义、实际价值、计算方法，做了最完整全面的介绍。目前所有流行的原子电荷计算方法在本文中都做了说明。本文写作力求深入浅出，在介绍原子电荷计算方法方面尽量避免用复杂的数学公式和抽象的概念，而尽可能用通俗易懂的语言讲清楚各种计算方法的基本思想、特点以及适合的使用范畴。这一章里还通过一批常见的分子，对不同原子电荷计算方法得到的结果的特点做了对比，使得读者可以直观认识到不同原子电荷计算方法结果的差异性。最后，本文还列举了一系列原子电荷的计算程序（可看到Multiwfn是目前支持原子电荷计算方法最多的程序），给出了原子电荷在实际研究中的使用建议，对研究者很有指导意义。限于篇幅和文章的定位，本文并未对原子电荷的算法和实际计算细节做展开讨论，非常推荐读者阅读超级详细的Multiwfn程序（//www.umsyar.com/multiwfn）的手册的3.9节，里面的内容在这些方面是对此综述的关键性补充。 总之，此综述是从头、完整了解原子电荷极为推荐优先阅读的文章，欢迎阅读和引用！此文可以按照以下格式引用（可根据期刊的要求修改） Tian Lu,* Qinxue Chen, Partial Charges, In Exploring Chemical Concepts Through Theory and Computation, WILEY-VCH GmbH: Weinheim (2024); pp. 161-187. DOI: 10.1002/9783527843435.ch6 另外，笔者2012年发表的《原子电荷计算方法的对比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）中对12种原子电荷进行了全面的对比，内容与上面的文章互相补充，也非常推荐阅读。此外，在笔者讲授的“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/workshop/WFN_content.html）中专门有一节极为全面讲授原子电荷的相关知识、实际计算和应用，内容远比上文广阔得多，欢迎关注和参加。 上文里对各种原子电荷的对比是对于孤立体系做的，关于原子电荷用于周期性体系的讨论，见《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）。 //www.umsyar.com/706 2024-05-05T14:31:00+08:00 2024年计算化学公社举办的计算化学程序和DFT泛函的流行程度投票结果 Results of the Computational Chemistry Commune 2024 poll on the popularity of computational chemistry programs and DFT functionals 文/Sobereva@北京科音  2024-May-5 0 前言 2024年4月4号，在北京科音建立的人气最高、专业性最强的综合性计算化学论坛“计算化学公社”（http://bbs.keinsci.com）上开展了为期一个月的五项投票： 你最常用的做 计算的DFT泛函投票（http://bbs.keinsci.com/thread-44387-1-1.html） 你最常用的做第一性原理计算的DFT泛函投票（http://bbs.keinsci.com/thread-44391-1-1.html） 你最常用的 程序投票（http://bbs.keinsci.com/thread-44388-1-1.html） 你最常用的分子动力学程序投票（http://bbs.keinsci.com/thread-44389-1-1.html） 你最常用的第一性原理程序投票（http://bbs.keinsci.com/thread-44390-1-1.html） 现对投票结果进行总结和简单评论。未来预计每三年重新开展一次投票。要强调的是，这个投票只是体现流行程度，和方法/程序的好坏并没直接关系。本投票结果主要反映中国的计算化学领域的专业、内行群体的情况，不反映业余/外行群体的情况。本次投票的结果也有助于业余/外行研究者正确认清什么才是主流，减少被他人忽悠、听信歪曲说辞误入歧途的概率。 上一次的投票于2021年举行，当时的结果和很多相关评论见下文： 2021年计算化学公社论坛“你最常用的计算化学程序和DFT泛函”投票结果统计 //www.umsyar.com/599（http://bbs.keinsci.com/thread-23482-1-1.html） 1 你最常用的做 计算的DFT泛函投票 本次可投的泛函有43种，带不带色散校正算同一种泛函。在2021年的投票条目基础上有所增减，特别是增加了双杂化泛函，明显几乎不会有人用的泛函没纳入可投范围。投票范畴仅限做 计算的情况，不包含第一性原理计算的情况。关系特别近的，比如M05-2X和M06-2X、wB97X和wB97XD和wB97X-D3、SCAN和r2SCAN、revDSD-PBEP86-D3(BJ)和DSD-PBEP86-D3(BJ)等等当做同一个泛函来计。此次投票者共713人。本投票每个人最多选6项，且所投的泛函必须占平时全部研究工作的5%以上。按照得票率（票数除以总投票人数）绘制的图如下。为了避免条目过多，只把得票前30名的列出。此图中诸如某泛函对应50%就代表有50%的人平时较多使用此泛函，后文的统计图同理。 总的来说本年度的投票结果和2021年时没太大变化，前六名的顺次没有改变，还是依次为B3LYP、(M05/M06)-2X、PBE0、wB97X-(/D/D3)、PBE、CAM-B3LYP老几样，各自有各自的用处（参看我对2021年投票结果的评论//www.umsyar.com/599，这里不再赘述）。它们的得票率的相对比例也基本没变，也就是量化领域里用处相对有限的PBE的比例有一定降低，以后肯定还得跌。2021年时候的第7名M06虽然得票率如今还是5%左右，但排名已下滑到第10，被wB97M-V、SCAN/r2-SCAN、PWPB95-D3(BJ)所超过。M06在我来看用处着实不大，虽然计算过渡金属配合物体系有一定用处，但PBE0-D3(BJ)/D4多数情况是更好的选择，而且M06还有后继者MN15可用。wB97M-V的得票率从2018年的3.1%提升到了如今的10%，已经算是增幅很快了，再过3年统计时肯定还会增高。此泛函在国内量化研究者中一定程度的流行，很大程度在于在计算化学公社论坛和思想家公社QQ群的讨论中时常被提及、在《简谈 计算中DFT泛函的选择》（//www.umsyar.com/272）博文中和我在北京科音基础（中级） 培训班（http://www.keinsci.com/workshop/KBQC_content.html）中的推荐、被免费的ORCA程序支持。提出时间不算很长的纯泛函SCAN及其改进版r2SCAN现在得票率能到6%着实不容易，2021年时得票率还不到1%，这主要在于有越来越多的程序已经支持此泛函，而且综合素质整体强于更早的经典泛函PBE，因而在纯泛函范畴内能有重要的位置。 2021年投票的时候没纳入双杂化泛函，这次得票率超过1%的双杂化泛函的排名顺序是PWPB95-D3(BJ)(5.9%) > (rev)DSD-PBEP86-D3(BJ)(3.1%) > B2PLYP-D3(BJ) (2.7%) > ωB97X-2-D3(BJ) (2.0%)。PWPB95-D3(BJ)和(rev)DSD-PBEP86-D3(BJ)能在国内用户中变得流行和上述wB97M-V的情况很类似。本身这俩泛函各有长处，有流行开来的必然性。PWPB95-D3(BJ)比较robust，算过渡金属配合物能量问题较好，能在ORCA里用；而revDSD-PBEP86-D3(BJ)算主族体系反应能、势垒以及弱相互作用能都是双杂化里顶尖的，而且在Gaussian里通过《Gaussian中非内置的理论方法和泛函的用法》（//www.umsyar.com/344）我介绍的方法能直接用。此外，ORCA中DSD-PBEP86适合算TDDFT和NMR目的也是其加分项。这俩泛函现在流行度能超过经典且最早提出的双杂化泛函B2PLYP是其应得的。 BLYP这次的排名降幅很大，从第10名已跌到第22名，本身这个泛函如今鲜有用武之地。BLYP一般也就算算主族体系，在Gaussian里用这个明显不如用B3LYP，耗时也持平，而以前在ORCA里用这个搭配def2-SVP结合RIJ加速做有机体系几何优化速度效率高是个用处，以前我在《大体系弱相互作用计算的解决之道》（//www.umsyar.com/214）里也提过，但如今也不如改用B97-3c来跑。 2 你最常用的做第一性原理计算的DFT泛函投票 可投的泛函有26种，带不带色散校正算同一种泛函。此投票在2021年没有，是本次新加的。此次投票者共442人。本投票每个人最多选6项，且所投的泛函必须占平时全部研究工作的5%以上。结果如下，零票的未显示 96年提出的PBE至今仍稳居第1的位置，如同B3LYP在 领域的地位，而且和第二名相差更悬殊。PBE能有这样的地位是必然的，PBE提出年代早、被程序支持得极为广泛，计算便宜，有丰富的专门为其搞的赝势，几何优化和分子动力学目的大多数时候够用（加色散校正后又拓宽了其普适性），而且在基态能量相关问题方面依然有使用价值而且没被已流行的其它纯泛函甩开特别多（这和B3LYP在 领域的情况截然不同，B3LYP算能量早过时了，很难再发得出去文章，见http://bbs.keinsci.com/thread-12773-1-1.html）。B3LYP在这次投票里得了第2，令我有点意外，大概率是很多人不好好看投票帖子的说明，误把 研究用的泛函也在此进行投票了。PBE0能排第3不意外，需要一个HF成分适当的杂化泛函做TDDFT/TDDFPT算激发态、算强相关体系的问题时经常用得着。HSE06流行度排得上第4主要来自于其计算带隙、能带方面公认很好，以及晶胞参数优化方面表现不错。PBEsol是优化晶体结构、晶胞参数的好把式，而且还是便宜的纯泛函，能排到第5很正常。M06-2X能排第6令我有点意外，一方面必定是有人误当成 计算的情况来投，另一方面是计算主族晶体/液体相关的化学反应、吸附的相关能量问题必定有一些人在用。SCAN/r2SCAN已经有一定流行度，由于在文献中出现频率越来越高，在未来的流行度必定也会逐渐提升，但流行度超越PBE在可预见的未来还不太可能，毕竟对于大量PBE就已经表现得够用的范畴，作为更贵但没带来显著优势的meta-GGA的SCAN/r2SCAN不可能显著侵占这方面的市场。第一性原理领域里用BLYP的人我不很理解是什么心态。revPBE和与之相似的RPBE能有一定流行度在于算化学吸附方面不错，被不少人用。第一性原理方面的泛函投票中CAM-B3LYP显得远不如 领域里来得流行，估计用这个的大部分都是CP2K用户用来算激发态和UV-Vis谱方面，只占投票的少量群体。算化学吸附很好的BEEF-vdW和算物理吸附很好的optB88-vdW能有一定票数不算意外。纯泛函中矮子里拔将军算带隙往往可以接受的HLE17在本次投票中获得了一点流行度，略意外的是算带隙整体更好些的纯泛函mBJ反倒在这次投票中显得无人问津，可能是前者在CP2K里能直接用而后者不能所致。作为PBE后继提出来的知名的TPSS流行度那么低有点令我意外，倒也确实实际用处不太大，现在又有了理论上以及实际整体表现得更好的r2SCAN。PW91虽然在文献里出现得很多，但这次得票率相当低，毕竟实际中有PBE就基本没有更老的PW91能派上用场的时候。B97M-rV能有少量票数，主要在于算热力学量方面在纯泛函里是表现得较突出的。 3 你最常用的 程序投票 可投程序有29种，投票者共679人。本投票每个人最多选三项，且所投的程序必须占平时全部研究工作的10%以上。按照得票率绘制的图如下，只显示了得票前20名的 前三位和2021年投票的结果一样，还是Gaussian > ORCA > xtb，Gaussian依然是约90%的 研究者日常必用的程序，稳稳占据主导位置。而ORCA和xtb的得票率比2021年时都有了约10%增长，这是意料之中的。实际上这三个程序也是我自己最常用的：xtb拿来快速预优化和结合molclus（http://www.keinsci.com/research/molclus.html）做构象搜索的初筛，Gaussian做基于DFT的opt、freq、扫描、IRC等涉及几何变化的任务以及算一些属性（NMR、超极化率等），ORCA算高精度单点。这三个程序的流行度远远甩开了其它程序，一方面是它们比较容易安装和使用，一方面各有各的独特优势，有被大量使用的刚性原因。而且它们把大部分 计算的应用领域都给覆盖了，对于日常应用性研究来说其它程序难以有和它们竞争的显著空间。Dmol3、ADF、Q-Chem、Turbomole这四个商业程序日子愈发不好过。在量化计算方面没有长处还巨贵的Dmol3的用户越来越少，从2021年的6.2%已经进一步萎缩到了4.3%，以后肯定还会明显进一步萎缩。ADF从2021年时候的仅仅2.2%进一步萎缩到了1.5%。Q-Chem从2021年的3.0%萎缩到了1.0%。Turbomole从2021年的1.6%萎缩到了1.0%。以GPU加速为卖点的TeraChem更不幸，2021年时候还有1人投票，今年变成了0人。 4 你最常用的第一性原理程序投票 可投程序有25种，投票者共565人。本投票每个人最多选三项，且所投的程序必须占平时全部研究工作的10%以上。按照得票率绘制的图如下（0票的没显示） 根据这次投票的结果可见，至少在计算化学公社论坛里，CP2K的流行程度已经赶超了VASP。这令我90%程度感到意外，但也有10%程度感觉是在情理之中。由于Multiwfn在2021年开始提供了极其易用的创建CP2K输入文件的功能（//www.umsyar.com/587），我后来又对此功能反复打磨并在Multiwfn中提供了对CP2K绘制DOS、能带、STM、电子激发、成键分析等许多功能，再加上2023年3月、2023年10月和2024年3月开办了三期北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）非常全面系统讲解了CP2K的使用，无疑令中国的CP2K用户猛增。但即便我已预料CP2K的得票率肯定会远高于2021年时候的33.3%，但也没预料到这次居然能达到60.9%，直接翻了将近一倍，甚至把一直霸占流行度榜首的VASP给挤下去了。由于免费且十分高效的CP2K的用户还在不断激增，而且CP2K更新很快，不断完善和发展新功能，Multiwfn在未来还会给其提供更多的相关分析处理功能，CP2K的位置在以后无疑还会更加牢固。相比之下，VASP的流行度从2021年投票时候的65.8%降到了现在的58.1%。由于现在有非常强大的竞争者CP2K，而且CP2K不具备的一些功能还有免费的Quantum ESPRESSO能用来平替VASP，售价较昂贵且算大体系速度远不如CP2K的VASP在未来的前景不乐观。以前很多人一提到第一性原理计算就言必称VASP，以后就不再是如此了。除了CP2K的流行度猛增外，其它程序的流行度都普遍出现了下降，如CASTEP和Dmol3分别从2021年的19.0%和9.3%下降到了13.8%和6.4%。Wien2k今年更是连一票都没有，而2021年时还有3票。那些程序流行度百分比减少，一方面是它们的票数大多数确实有一定减少，用户发现有更理想或免费的程序可用，另一方面原因应当是有很多通过CP2K程序开始入手第一性原理计算的人加入了投票，他们只对CP2K的得票率有贡献而间接地拉低了其它程序的得票率。 5 你最常用的分子动力学程序投票 可投程序有18种，投票者共551人。本投票每个人最多选三项，且所投的程序必须占平时全部研究工作的10%以上。按照得票率绘制的图如下 GROMACS依旧是用户数的龙头老大，而且流行度从2021年投票时的69.3%还进一步略微提升到了71.3%，得票数大约等于其它所有程序用户数之和，和曾经的情况一样。第2、3位依然分别是Lammps和AMBER。Lammps和OpenMM得票率略涨，而AMBER、Forcite和NAMD的流行度都有较多降幅，GULP、DL_POLY和CHARMM更是快跌没了。 //www.umsyar.com/703 2024-03-03T12:07:24+08:00 深入理解分子轨道对磁感生电流的贡献 Deep understanding of contribution of molecular orbitals to magnetically induced currents 文/Sobereva@北京科音   2024-Mar-3 0 前言 施加外磁场会导致化学体系产生磁感生电流，考察磁感生电流是研究分子芳香性的重要方法。我之前写过诸多博文介绍过，如《使用SYSMOIC程序绘制磁感生电流图和计算键电流强度》（//www.umsyar.com/702）、《考察分子磁感生电流的程序GIMIC 2.0的使用》（//www.umsyar.com/491）、《使用AICD 2.0绘制磁感应电流图》（//www.umsyar.com/294）、《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）等。最常见的是考察整个分子的磁感生电流，而对于苯等平面体系，也经常区分成sigma和pi电子的感生电流分别考察，以更充分了解不同类型电子的离域性、对芳香性的贡献。对于CTODZ和与之本质相同的CSGT计算感生电流方法，总的感生电流是可以精确写为所有占据轨道贡献之和的形式的，也因此原理上可以给出任意占据轨道对应的感生电流。轨道的感生电流会带来感生磁场，对外磁场有屏蔽或去屏蔽效果，因此相应地可以得到各个轨道对NICS的贡献，如《基于Gaussian的NMR=CSGT任务得到各个轨道对NICS贡献的方法》（//www.umsyar.com/670）所展示的。 然而，单一轨道对感生电流和NICS的贡献实际上是存在任意性的，过度解释可能得到十分误导性的结论。本文就专门详细谈一谈这个问题，以令读者正确理解轨道对感生电流的贡献。并且还同时介绍一下如何根据占据轨道到空轨道之间的跃迁模型来预测和解释一些体系的磁感生电流。本文许多内容和讨论取材自J. Phys. Chem. A, 105, 9553 (2001)一文，这是一篇很经典的探讨轨道与磁感生电流关系的文章。但这篇文章里很多部分写得难懂，本文则尽可能说得清楚、明白、易懂。 1 分子轨道对感生电流贡献的表达式 此文只讨论闭壳层的情况，假设体系有N个电子，因此有N/2个占据轨道。如下式(1)所示，总的感生电流可以写为所有占据轨道的加和形式，n是轨道序号。下式(2)给出了n轨道对r处感生电流矢量的表达式，B是外磁场矢量，d是测量原点（gauge origin），φ是无磁场下计算的分子轨道波函数，φ(1)是以二阶微扰理论得到的外磁场作用下的φ的一阶校正。 由上式可见，轨道的感生电流由paramagnetic（顺磁）和diamagnetic（抗磁）两部分构成，以下分别简称para和dia。dia电流是楞次定律说的给导体施加外磁场时产生的那种感生电流，正比于轨道上的电子的概率密度，它产生的感生磁场会在环形感生电流内部区域对外磁场产生削弱作用。而para部分完全来自于量子力学效应，直接牵扯到轨道波函数，方向与dia相反。 在完备基组下，体系总的感生电流计算结果是不依赖于d的选取的，但构成它的dia和para电流各自大小却依赖于d，因此dia和para的划分没有严格的物理意义、也不存在唯一的划分。而对分子轨道而言，不仅d的选取会影响此轨道的dia和para电流，对二者的总和也会产生影响。上式的dia部分里直接就有d，而para部分的φ(1)当中的φ(d)也依赖于d，具体表达式如下 ε是分子轨道能量，p是单电子动量算符，l(0)=r×p是相对于笛卡尔坐标原点的单电子角动量算符。φ(p)和φ(d)分别涉及当前轨道向各个空轨道的转动跃迁和平动跃迁。 如果计算所有位置的感生电流时使用统一的d，叫做common gauge-origin。这种做法平时基本不用，因为d的选取任意性太强而又会直接影响结果（也就某些情况有相对适合、意义较强的，比如算苯的时候可以把d设为环中央）。如果计算r处的感生电流用的d可表达为d(r)形式，即d也是像r一样连续变化的且依赖于r，则被称为CTOCD方法。SYSMOIC程序就是基于此方法计算感生电流。CTOCD有不同具体形式，其中CTOCD-DZ方法的名字中的DZ代表diamagnetic zero，它将d设为r，即算哪个位置的感生电流就把此时用的d设在哪里，使得上面(2)式中的dia项为0。虽然不能说把d设为r就有什么确切的物理意义，但起码没有common gauge-origin的d选取的表面上的任意性。 下面我把一般情况、d=r的情况、d=0的情况的感生电流表达式放在一起，便于大家弄清楚它们之间的异同、para和dia的定义。 CTOCT-DZ方法由于d=r，所以上面“一般情况”中的B×(d-r)项就没了，再把φ(1)=φ(p)+φ(d)代进去就得到了CTOCT-DZ的表达式。可见其中一部分由涉及平动跃迁的φ(d)所贡献，这部分被习俗地视为CTOCT-DZ的dia项，另一部分由涉及转动跃迁的φ(p)所贡献，这部分被习俗地视为CTOCT-DZ的para项。而对应d=0的common gauge-origin的情况中φ(d)就不复存在了，φ(1)=φ(p)，就有了上面的式子。由上式可见，CTOCT-DZ和d=0的情况的para项是相同的，但是dia项不同，因此两种情况计算的轨道感生电流是不同的。 基组够大的情况，总的感生电流近乎不依赖于d，因此可以用总的感生电流讨论芳香性。而由于轨道的感生电流明显依赖于d，所以绝对不能轻易拿一个轨道的感生电流的特征去判断它对体系芳香性、电子离域有什么样的贡献，否则会造成误导。但这不代表做感生电流的轨道层面的分解讨论就没实际意义，因为完备基组下某些轨道的感生电流的总和也是不依赖于d的。对纯平面体系，sigma电子的感生电流和pi电子的感生电流可以精确分离，且二者都是不依赖于d的。具体为什么，是因为存在以下关系 由上面最后一个式子可见，某分子轨道的CTOCD-DZ和d=0的感生电流只相差一部分，它依赖于此轨道与其它占据轨道之间的平动跃迁。当所有占据轨道的感生电流加和时，n轨道的感生电流表达式中的n→m跃迁部分的贡献和m轨道的感生电流表达式中的m→n跃迁部分的贡献精确抵消了（因为轨道能量差是φ(d,occ)的分母项，ε_n-ε_m和ε_m-ε_n的符号相反），所以CTOCD-DZ和d=0得到的体系总感生电流是相同的。为什么对于纯平面体系pi（或sigma）电子的总感生电流也不依赖于d？这是因为由于对称禁阻原因sigma与pi轨道之间的<φ_m|p|φ_n>项精确为0（即在垂直于外磁场方向即便平移sigma轨道后也仍然与pi轨道正交），φ(d,occ)因此又可以分为只循环sigma轨道和只循环pi轨道部分的加和，即sigma和pi是两个不同的block。所以所有pi轨道的感生电流加和时由于所有牵扯到的n→m和m→n（n和m皆为占据的pi轨道）的贡献精确抵消原因，它不依赖于d。 2 苯分子的不同pi轨道的感生电流 这一节看一个实例，以更好地了解上一节的讨论。下图给出了苯的最低pi轨道（MO 17）以及两个简并的HOMO pi轨道（MO 20 21)各自的Multiwfn绘制的0.05 a.u.轨道等值面图，以及SYSMOIC程序绘制的分子平面上方1 Bohr处的感生电流图。所有三个轨道加和对应的总的pi感生电流图也一起给出了。计算分别使用CTOCD-DZ形式（对应SYSMOIC里JBMAP结合-m DZ1选项），以及d=0的common-gauge情况（对应SYSMOIC里JBMAP结合-m CO选项，默认的d坐标(0,0,0)此时对应苯环中心）。几何优化和产生波函数文件都是在B3LYP/6-31G*下进行的。 由上图可见，CTOCD-DZ和d=0的总的pi感生电流基本没有可查觉的差别，但是这两种方法算出来的不同pi轨道的感生电流就相差显著了。CTOCD-DZ的MO 17的感生电流极小，在当前作图设置下甚至都看不到，而d=0时它的感生电流则很显著，和其它pi轨道的差不多。其原因在于MO 17与低阶的空轨道之间的平动和转动跃迁都是禁阻的，只可能往能量很高的空轨道跃迁，又由于轨道能量差出现在φ(d)和φ(p)的分母项，自然CTOCD-DZ方法算出来的它的感生电流特别小。而反之，d=0方式计算轨道感生电流的公式里有B×rφ^2项，φ^2对应当前轨道概率密度，因此由于MO 17在整个六元环上分布得很均匀且密集，使得MO 17的d=0方法算的感生电流很显著。也由于d=0的感生电流计算公式的这个特点，MO 20和21的感生电流显著分布区域也和轨道波函数显著分布的区域高度类似，而CTOCD-DZ的感生电流分布和轨道分布则没有视觉上的显著联系。 对MO 17，CTOCD-DZ和d=0的图像展现的信息有天壤之别。如果不懂以上原理的话，看CTOCD-DZ图就会盲目以为MO 17的电子几乎没有在六元环上的离域性、跟pi芳香性没联系。而事实上，在碳环上完全均匀分布的MO 17显然对苯的pi电子的全局离域有重要贡献。所以，特定轨道的感生电流和轨道实质上体现的电子离域特征并没有必然关联，这点一定要注意！ 那么是否轨道的感生电流一定不能用于讨论轨道对芳香性的贡献呢？在我来看，这种讨论方式虽然由于上述原因而非常不推荐，但在特殊情况下还是可行的，比如上面d=0方法算的MO 17的感生电流确实体现出MO 17对芳香性有明显贡献。这里我阐述我的个人观点：如果一个轨道对电子在某个环上的多中心离域完全没有贡献，那么无论怎么选择d也算不出它有连贯环绕整个环的dia电流，例如内核电子轨道；而如果能有特定的d可以算出来它有连贯环绕整个环的dia电流，那么这个轨道对多中心离域就是有明显贡献的，如苯的MO 17。从前面的图中还可以看到，d=0时的MO 20和21各自的感生电流都没有在整个苯环上连成整体，但在CTOCD-DZ的情况下则连成了整体，根据以上我的观点，这直接说明了这俩轨道确实对pi电子的六中心离域有直接贡献。 CTOCD-DZ算的MO 17的感生电流极小，也很大程度反映在《基于Gaussian的NMR=CSGT任务得到各个轨道对NICS贡献的方法》（//www.umsyar.com/670）文中以CSGT方式算的苯的各个分子轨道对NICS(1)zz的贡献上，MO 17的贡献（-1.58 ppm）远低于MO 20和21的贡献（-14.07 ppm）。CSGT方法是使用(2)式计算各处的感生电流密度，用Becke多中心积分方法通过Biot-Savart定律计算特定位置的磁屏蔽张量，每个积分格点用的d是此处的坐标往距离它最近的原子核一定程度移动得到的（基于各个原子的Becke权重确定）。SYSMOIC程序给出的CSGT方法算的MO 17的感生电流密度和CTOCD-DZ一样也是微乎其微的。 3 基于轨道跃迁模型的解释感生电流的例子 CTOCD-DZ算的轨道感生电流没法展现出苯的最低pi轨道对六中心离域的贡献，是否说明CTOCD-DZ没啥实际优点？答案是否定的。一方面，CTOCD-DZ计算轨道感生电流时给定了明确的d的选取方式，因此消除了人为选择的任意性；另一方面，CTOCD-DZ计算轨道感生电流的公式是可以严格分解为当前轨道向不同空轨道的跃迁贡献的，这种形式被称为完全态求和（sum-over-state, SOS）。SOS形式使得研究者可以从轨道间的平动和转动跃迁角度来定性预测不同轨道的感生电流该是什么样，以及深层次解释和理解CTOCD-DZ实际算出来的轨道感生电流为什么是那样，从而获得重要的physical insight。其实还有不少其它的分子属性也可以表达为SOS形式，例如我在《使用Multiwfn基于完全态求和(SOS)方法计算极化率和超极化率》（//www.umsyar.com/232）里介绍了（超）极化率就可以通过SOS形式来计算，并且经过简化，还可以得到两态和三态模型来分析讨论影响第一超极化率的关键因素，见《使用Multiwfn对第一超极化率做双能级和三能级模型分析》（//www.umsyar.com/512）和《谈谈计算第一超极化率的双能级公式》（//www.umsyar.com/361）。 对于CTOCD-DZ方法计算的占据轨道n的感生电流，n向某个空轨道a的跃迁若能够对它有显著的贡献，需要满足三个条件 (1)对称匹配。n→a必须要么平动跃迁允许（对dia电流有贡献），要么转动跃迁允许（对para电流有贡献），要么都允许。对于有点群对称性的体系来说，如式(3)所示，平动跃迁如果是允许的就意味着<φ(a)|p|φ(n)>不能精确为0，这要求n与a轨道不可约表示的直积对应的表示包含垂直于磁场的平移。转动跃迁如果是允许的就意味着<φ(a)|l(0)|φ(n)>不能精确为0，这要求n与a轨道不可约表示的直积对应的表示包含绕着磁场方向的旋转 (2)轨道重叠显著。如果a和n的空间分布没有明显的重叠，自然平动和转动跃迁肯定非常小、n→a跃迁不可能对n的感生电流有可查觉的贡献 (3)能量相近。由于a和n轨道的能量差出现在式(3)的分母部分，因此只有较高的占据轨道与较低的空轨道之间的跃迁才可能对感生电流有明显贡献 下图是JPCA文中在HF/6-31G*级别算的D6h点群的苯的前线轨道能级图，小圆点是占据的电子，黑色箭头是平动允许的轨道跃迁，苯处在XY平面上。例如e1g与e2u之间的直积是b1u+b2u+e1u，查特征标表可知e1u包含(x,y)，与Z方向的外加磁场相垂直，因此平动跃迁允许。下图所标注的a2u到e1g、e2u到b2g也类似地都是平动跃迁允许的。这些轨道之间没有转动跃迁是允许的，因为它们的不可约表示的直积都没有包含绕Z轴旋转的行为（Rz）。对于CTOCD-DZ方法，前面说过它只涉及占据轨道到空轨道的跃迁，因此此体系的感生电流仅来自于下图中e1g（HOMO）到e2u（LUMO）的跃迁，而且完全是dia电流。D6h点群的不可约表示里只有e1u含有x和y行为，能量最低的占据pi轨道是a2u，仅和e1g直积可以得到e1u，而空轨道里最低的e1g轨道是能量很高的，因此这个占据pi轨道向e1g轨道跃迁导致的对dia电流的贡献可忽略不计，这是为什么上一节看到CTOCD-DZ算的苯的MO 17的感生电流在图上都看不到。 下图是CTOCD-DZ计算的pyracylene的感生电流以及轨道跃迁示意图，分子在YZ平面上，外磁场由屏幕内向着屏幕外顺着X方向加（图片取自前述JPCA文章，为了和我的习俗相同，我把文中的感生电流方向反转了一下，顺时针和逆时针分别对应dia和para电流）。下图右侧示意图中实箭头对应允许的平动跃迁（因为b3u与b1g的直积b2u包含y，au与b1g的直积b1u包含z），故b3u和au轨道贡献顺时针的dia电流；空心箭头对应允许的转动跃迁（因为b2g与b1g的直积b3g包含Rx），故b2g轨道贡献逆时针的para电流。由于更低阶pi轨道对感生电流贡献甚微，所以总的pi感生电流的图像特征基本等于b3u和au贡献的dia电流以及b2g轨道贡献的para电流的总和。由于同时有4个电子贡献dia电流而只有2个电子贡献para电流，故前者整体强度更大，这是为什么在这两种电流方向相反的分子的上、下区域，总电流体现的是顺时针特征。 如果不通过群论知识，怎么判断什么样的轨道之间的转动跃迁是比较显著的呢？这可以结合图形进行判断和理解。若绕着磁场旋转一个轨道后就可以令二者不再正交，那么它们之间的转动跃迁就是允许的，这通常出现在两个轨道的节面数目是相似的情况下。例如下图所示的pyracylene的HOMO (b2g)→LUMO (b1g)的跃迁就是如此，最右边的图是把HOMO旋转90度后和LUMO的重叠图，可见相位完全匹配且高度重叠，因此这种跃迁能够对para电流有显著贡献。 如果在垂直于磁场方向平移一个轨道后能够与另一个轨道之间不再正交（一个轨道比另一个轨道多一个节面往往是这种情况），那么二者间的平动跃迁就是允许的。例如下图所示的pyracylene的HOMO-1 (au)→LUMO (b1g)的跃迁就是如此，HOMO-1有两个节面，LUMO有三个节面，最右边的图是把HOMO-1垂直于磁场平移后和LUMO的重叠图，可见能是以相位匹配的方式重叠，重叠积分明显不为0，因此这种跃迁能够对dia电流有显著的贡献。 通过以上例子，可以看出利用占据-非占据轨道跃迁模型，确实可以对CTOCD-DZ方法计算的轨道感生电流予以预测，还能对总感生电流的特征从轨道角度进行深层次的解释。前述JPCA文章里还有更多分析例子，包括苯的阴/阳离子、萘、并六苯、晕苯、碗烯，感兴趣者可以看看。文章还总结出CTOCD-DZ计算的感生电流的规律，对满足4n+2的普通碳单环体系，只有处于二重简并的HOMO轨道上的4个pi电子对pi感生电流有明显贡献，而更低的轨道的贡献可以忽略。而对于更复杂的多环分子，只有少量的前线pi占据轨道对感生电流有显著贡献，这是因为低阶pi占据轨道往平动/转动允许的pi空轨道上跃迁对应的能级差太大而必定贡献很小，而低阶pi占据轨道往前线pi占据轨道跃迁即便是平动/转动允许的也没用，因为占据-占据轨道的跃迁贡献不在CTOCD-DZ方法中出现）。 //www.umsyar.com/702 2024-02-29T10:22:00+08:00 使用SYSMOIC程序绘制磁感生电流图和计算键电流强度 Using SYSMOIC program to plot magnetically induced current map and calculate bond current strength 文/Sobereva@北京科音  2024-Feb-29 0 前言 笔者在《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）中介绍过，磁感生电流是一种重要、直观的考察分子芳香性的方式。SYSMOIC是一个基于continuous transformation of the origin of the current density (CTOCD)方法考察化学体系的磁感生电流相关问题的程序。本文将结合苯分子为例，介绍这个程序的最关键用法。SYSMOIC程序和相关原理的详细介绍可以看原文J. Chem. Inf. Model., 61, 270 (2021)。 AICD和GIMIC也是两个重要的考察磁感生电流的程序，笔者在下文有专门的介绍。这些文章里介绍的背景知识对理解本文很有帮助 • 使用AICD 2.0绘制磁感应电流图（//www.umsyar.com/294） • 使用AICD程序研究电子离域性和磁感应电流密度（//www.umsyar.com/147） • 考察分子磁感生电流的程序GIMIC 2.0的使用（//www.umsyar.com/491） • Utilizing AICD and GIMIC programs to study magnetically induced current density（//www.umsyar.com/148） 之所以笔者写此文介绍SYSMOIC，在于它对以上两个程序有明显互补性。AICD程序是将感生电流以箭头方式绘制在AICD的等值面上，而等值面以外区域的感生电流就无法在图上体现出来了，而且AICD没法给出感生电流的键截面积分值（键电流强度）来定量考察电流大小。SYSMOIC则能够给出三维空间或二维平面上的感生电流向量场图，从而能展现出AICD图表现不出来的信息，而且SYSMOIC还能非常方便地给出键电流强度。GIMIC程序虽然能做感生电流积分、结合Paraview程序能给出漂亮的向量场图，但是GIMIC却无法将感生电流分解成不同分子轨道的贡献，而SYSMOIC则能够方便地指定计算感生电流时只考虑或不考虑哪些分子轨道。可见SYSMOIC有独特的价值。 本文介绍的情况对应SYSMOIC 2024.1版，对其它版本不一定适用，请结合相应版本程序提示和手册随机应变。 1 SYSMOIC的基本特征 SYSMOIC是免费而不开源的程序，Windows、Linux、Mac版都有，可以在http://SYSMOIC.chem.unisa.it/DISTRIB下载。诸如STABIN_LNX-2024.1.tar.gz就是Linux的2024.1版可执行文件包。Linux版需要用较新的系统，笔者用Rocky Linux 9可以正常使用，CentOS 7.4用不了。 SYSMOIC是由一堆子程序组成的，对应压缩包里不同可执行文件，Windows版有.exe后缀，其它平台的没有。例如JBMAP子程序用于计算感生电流格点数据，BOCUST子程序用于计算键电流强度，等等。v3d子程序用来对其它子程序计算出来的.3d文件进行可视化，如绘制感生电流矢量图。各个程序的用法在SYSMOIC的在线手册http://sysmoic.chem.unisa.it/MANUAL/里都有说明。 SYSMOIC是解压即用的程序。解压后只需要把SYSMOIC目录加入到操作系统的PATH环境变量里（不会的话自行google），然后就可以在任何目录下输入子程序名直接使用了。 SYSMOIC有不少官方例子文件，是http://sysmoic.chem.unisa.it/WORKED_EXAMPLES/里的WORKEXA.tar.gz。例子按照体系/理论方法/基组对目录层级进行命名，里面RUN_TESTS文件是对例子体系做各种分析的shell脚本，可以从中看到都是以什么命令调用程序的，以及传入的命令文件是什么。REFERENCES目录下的文件是运行RUN_TESTS脚本后产生的.3d文件和程序输出信息，对前者可以直接用v3d作图。 SYSMOIC必须结合最主流的 程序Gaussian使用，目前只支持闭壳层形式的HF和DFT。只需要用Gaussian做NMR计算的同时产生wfx文件，再用SYSMOIC自带的unpackwfx将之转化成fort.3、fort.11、fort.23、fort.28这四个二进制文件，之后就可以用JBMAP、BOCUST等子程序做各种计算了，算完了可以用v3d子程序直接作图。 2 产生苯的SYSMOIC输入文件 下面以最典型的芳香性分子苯为例，讲解SYSMOIC绘制感生电流矢量图和计算键电流强度的方法。其它的用法我觉得在一般研究中用到的机会很少，所以就不在本文介绍了，感兴趣者请看程序手册、原文并结合程序自带的例子文件摸索。下文涉及的各种输入输出文件都可以在//www.umsyar.com/attach/702/file.rar中获得。 首先写一个在B3LYP/6-31G*级别下做NMR计算的Gaussian输入文件benzene.gjf，内容如下，在本文的文件包里已经提供了。其中的坐标是在当前级别下已经优化过的。NMR=CSGT要求以CSGT方法做NMR计算，out(wfx,CSGT)要求产生.wfx波函数文件，并且把CSGT计算过程中产生的扰动的轨道展开系数写入到wfx文件中，输出位置为末尾指定的D:\benzene.wfx。 # b3lyp/6-31g(d) NMR=CSGT out(wfx,CSGT) [空行] benzene [空行] 0 1  C                  0.00000000    1.39621600    0.00000000  C                  1.20915900    0.69810800    0.00000000  C                  1.20915900   -0.69810800    0.00000000  C                  0.00000000   -1.39621600    0.00000000  C                 -1.20915900   -0.69810800    0.00000000  C                 -1.20915900    0.69810800    0.00000000  H                  0.00000000    2.48297800    0.00000000  H                  2.15032200    1.24148900    0.00000000  H                  2.15032200   -1.24148900    0.00000000  H                  0.00000000   -2.48297800    0.00000000  H                 -2.15032200   -1.24148900    0.00000000  H                 -2.15032200    1.24148900    0.00000000 [空行] D:\benzene.wfx [空行] [空行] 算完后就有了benzene.wfx，确保它在当前目录下，并且已经把SYSMOIC目录加入到了PATH环境变量中，然后运行unpackwfx benzene命令，当前目录下立刻就出现了fort.3、fort.11、fort.23、fort.28，它们都已提供在了本文的文件包里。后文说的所有操作都是对于这四个fort文件都处于当前目录下的情况而言的。 注意当前的苯分子是处在Z=0的XY平面上，用Multiwfn载入Gaussian的out文件然后进主功能0，或者用GaussView载入out文件并显示笛卡尔轴，都可以确认这一点。后面我们要考察的都是外磁场方向垂直于苯环所导致的感生电流的情况。 3 绘制磁感生电流向量场图的例子 3.1 默认方式绘制 先来看怎么产生一个最简单的磁感生电流量场图，之后再解释细节。直接在命令行窗口输入JBMAP运行之，屏幕上会出现好多默认用的参数，比如下面这行，就是告诉你默认的外磁场是(0 0 1)矢量，显然正垂直于当前苯环，和我们期望的一致 magnetic field...........B    0.0000    0.0000    1.0000 程序问你当前显示出来的参数是否ok，这里就直接按回车代表用默认的y（yes）。然后再依次问你三个问题： • external grid points y/n? [n]：即是否读取外部格点，这里直接按回车代表用默认的n（如果要计算的点不是均匀分布在一个矩形区域，而是比如分布在一个球层表面，就需要输入y并且提供格点定义文件） • reference arrow y/n? [n]：含义不明。直接按回车用默认的n • B direction y/n? [n]：即是否把外磁场箭头画在图上，直接按回车用默认的n。如果你想把外磁场箭头绘制在图上就输入y，并且输入箭头的位置和长度 马上当前目录下就有了JBM.3d文件。输入v3d JBM.3d通过v3d程序对这个文件作图，就蹦出了图形窗口，如下所示，可以看到在一个矩形区域内均匀分布的每个格点位置都有一个箭头指示此位置的磁感生电流的方向和大小。 想退出的话在图形窗口按ESC键，或者在文本窗口按ctrl+C。 上面的图非常丑陋，箭头特别拥挤，分子结构也看不清楚，完全没法用。下一节讲怎么对JBMAP的计算进行自定义。 3.2 自定义JBMAP的运行参数 JBMAP有两个方式指定运行参数，比如要把外磁场矢量改为0 1 1（程序自动会做归一化成为0.0000 0.7071 0.7071），既可以运行JBMAP -B 0 1 1命令，也可以先启动JBMAP然后直接输入B 0 1 1后按回车，在屏幕上显示的参数中都会看到外磁场矢量已经变更成自设的值了。 JBMAP命令后面可以直接接的选项可以输入JBMAP -h查看，其中几个关键的在这里提一下： -o：输出的.3d文件名，默认是-o JBM -f：原子球大小的倍数，默认为-f 0.2。想让作出的图里原子球大一些就把这个改大 -B：外磁场矢量，默认为-B 0 0 1 -j：感生电流类型，默认是-j TOT，即diatropic和paratropic的总和。写-j PARA代表只考虑paratropic部分，写-j DIA代表只考虑diatropic部分 -q：考虑的轨道。例如只想考虑17 20 21三个分子轨道对磁感生电流的贡献，就写-q +3 17 20 21，只去掉15 16号分子轨道的贡献就写-q -2 15 16。默认是考虑所有分子轨道的贡献（只有占据轨道才有贡献） -m：设置SYSMOIC计算感生电流用的CTOCD方法的具体形式，参看http://sysmoic.chem.unisa.it/MANUAL/S4.SS12.html的说明。当基组是完备的时候这些方法结果等价。有些选项是仅对于某些方法才需要设，比如用common gauge-origin的时候需要用-g自定义原点坐标（默认为0 0 0）。默认用的CTOCD-DZ2基本上对于所有情况都是合适的选择，也不用设额外参数。关于CTOCD不同形式的介绍和差异，可以看前述的SYSMOIC原文以及感生电流综述文章WIREs Comput Mol Sci, 6, 639 (2016)的Continuous Transformations of the Origin of the Current Density一节。 JBMAP运行后可以修改的参数当中较重要的如下： • B：外磁场矢量，默认B 0 0 1 • JTERM：感生电流类型，默认是JTERM 0，即diatropic和paratropic的总和。1和2分别代表paratropic和diatropic • RI和RF：分别指定格点所均匀分布的矩形盒子的X、Y、Z最小笛卡尔坐标和最大笛卡尔坐标。默认是根据体系坐标自动确定以覆盖整个体系。对上面的苯分子的例子，从屏幕上默认的参数可见，当前对应的是RI -6.0635 -6.6921 -2.0000和RF 6.0635 6.6921 2.0000 • FATT：给感生电流乘的系数。比如大部分箭头在图上都显得很不显著，就可以设大这个值以更便于观看 • FMM：磁感生电流过滤范围。例如FMM 0.01 1就代表忽略感生电流的模小于0.01和大于1的格点（是对乘上FATT之前而言的），免得干扰观看。默认0 0.1 • STEP：格点间距，数值越小箭头越密、计算越耗时。默认0.4 Bohr • LUNA：0和1分别代表箭头的长度和面积正比于感生电流的模，默认是后者 • METH：计算方法，默认的2对应CTOCD-DZ2 为了方便，我推荐以命令行方式运行JBMAP。交互式程序怎么以命令行方式运行，我在《详谈Multiwfn的命令行方式运行和批量运行的方法》（//www.umsyar.com/612）里有很详细说明。下面以这种方式结合自定义参数用JBMAP再次对苯进行计算。创建一个文本文件叫CDMAP.txt，内容如下，每一行对应一条传给JBMAP程序的命令。此例RI和RF的Z值都是1，X和Y范围都是-6到6 Bohr，因此格点会分布在苯分子上方1 Bohr的XY平面上的正方形区域内。此文件后面的y、n、n、n对应交互式界面里程序向你提问时依次输入的四个指令。 B   0 0 1 RI -6 -6 1 RF  6  6 1 FMM 0.01 1 FATT 3 STEP 0.4 LUNA 1 JTERM 0 y n n n 运行JBMAP -f 0.5 < CDMAP.txt，然后再运行v3d JBM.3d，看到下图，可见清楚地把苯分子上方的磁感生电流描绘了出来。 再来看只考虑pi轨道时的磁感生电流图的绘制。把benzene.wfx文件载入Multiwfn，进入主功能0，按照《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）说的，从HOMO开始往下挨个看轨道，就能找到三个pi轨道，序号为17、20、21。原子数多的时候也可以用《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）介绍的做法让Multiwfn自动找出来所有pi轨道序号。运行JBMAP -q +3 17 20 21 -f 0.5 < CDMAP.txt，然后再用v3d绘图，得到下图，这便是苯的pi电子的感生电流图了。相对于之前的所有电子贡献的图，此图少了sigma电子造成的出现在内环的逆时针（paratropic电流）特征，而且也没有了氢原子区域的感生电流，因为pi电子不在氢上。 3.3 v3d的使用细节 下面专门细谈一下v3d界面的操作。v3d图形窗口出现后，可以用鼠标拖动体系旋转视角。还可以在图形窗口处于激活状态下按h键，此时文本窗口会出现可以用的所有按键的功能说明。其中较重要的按键在这里说一下。 先说视角调节按键： a、b、g：视角由alpha、beta、gamma值定义，按a、b、g键会分别增加它们1.0。如果按住shift键再按，也即按A、B、G，会分别增加它们10 • f和F：调节缩放因子。按一下f图像就会放大一点，按一下F就会放大很多 • d和D：调节视角距离。按一下d就会增加0.1，按一下D就会增加1。视角距离越远，近大远小的透视效果越弱、越接近正交视角的效果（据我所知v3d没法完全用正交视角） • n和N：设置雾化距离阈值。超过视角一定距离阈值开始就会用雾化效果使得远处物体被屏蔽。按n增加阈值0.1，按N增加1.0 • C：令图像居中 每次按以上按键调节设置的时候在文本窗口都会看到当前值。如果你想把这些键的操作效果反过来，比如按一次f就令图像缩小一点，那就先按一下c键进入反转操作状态（文本窗口里会看到c变成了c-1），然后再按f就行了。再次按c则还原到普通操作状态。 按p键可以把当前图像保存成当前目录下的.ps文件。这是矢量图像格式，在Windows下可以用装了ghostscript插件的irfanview或photoshop等程序打开，也可以用Acrobat distiller转成pdf再打开。在Rocky Linux、CentOS等系统里也自带了ps文件观看程序。 默认情况下产生的ps文件画出来的感生电流的箭头会有较厚的白边，当箭头多的时候白边会导致后面的箭头被遮挡，比较难看。解决办法是按一下c键先进入反转操作状态，然后按住k键一会儿使得白边量减小（文本窗口不会提示减小到了多少），然后再导出ps文件，然后就会发现白边非常细了。 如果想保存当前的作图设置，按S，当前目录下就出现了SpecialValues.txt文件，里面记录了作图设置，可以看到f、d等数值都记录在内。当再次启动v3d时，程序若发现当前目录下有这个文件，就会自动载入里面的设置。如果你希望每次产生的ps文件里箭头的白边都非常细，可以启动v3d后先直接按S保存一个记录默认状态的SpecialValues.txt文件，然后把里面对应白边粗细的ridfr前面的值改为0.1。这样每次启动v3d就会默认用这个设置了。 4 计算键电流强度 对键截面做感生电流积分的概念我在《考察分子磁感生电流的程序GIMIC 2.0的使用》（//www.umsyar.com/491）中已经专门讲过了，这可以定量考察穿越特定化学键的电流强度从而便于对比，这在SYSMOIC里叫键电流强度。SYSMOIC算键电流强度远比用GIMIC更方便，都省得自己定义积分范围，结果还能直接可视化。下面我们对苯的一个C-C键计算键电流强度。 直接在命令行窗口中输入BOCUST命令运行此程序，程序首先会提示 | plotted    12 atoms | plotted    12 bonds  ====> Pair C1    -  C2     distance=    1.396 Ang 这说明程序发现有12个原子对儿的距离都处在默认的距离下限（0.9埃）和距离上限（1.6埃）之内。其中C1-C2是第一个原子对，现在问你是否对这个C1-C2键计算键电流强度。这里输入y然后按回车，很快就会出现结果 CONTRIBUTION  1 -0.6000526E+00 AU -0.1690912E-07 SI CONTRIBUTION  3  0.1784109E+00 AU  0.5027513E-08 SI TOTAL CS ==>    -0.4216417E+00 AU -0.1188161E-07 SI (2) CS/CS0 RELBEN FICS     0.99    3.00  0 上面显示的两个CONTRIBUTION分别是diatropic和paratropic对键电流强度的贡献，由于方向相反，显然符号相反。二者之和取绝对值0.4216417 a.u.就是键电流强度大小，换算成常用的nA/T单位就乘以28.179409，即11.88 nA/T。这个值与《考察分子磁感生电流的程序GIMIC 2.0的使用》（//www.umsyar.com/491）在B3LYP/6-311G*级别通过GIMIC程序得到的12.086 nA/T非常接近。上面信息中的CS/CS0 RELBEN FICS     0.99代表当前值与参考值的比值是0.99。这里参考值是程序内置的对苯在某个级别计算的C-C键电流强度12 nA/T。 接下来程序又问你是否对C1-C6计算键电流强度。由于当前所有C-C键都是等价的，只需要计算一个就行了，因此输入N回车，程序就结束运行了。如果你想把所有键都进行计算，就接着输入11次y回车把剩下11个键都算完。当前体系一共6个C-C键、6个C-H键。 当前目录下出现了BCS.3d文件。运行v3d BCS.3d，就可以看到下图 此图里黑色箭头体现了C1-C2键电流的大小和方向。积分截面处在键的正中央并与键垂直，图中的曲线是这个截面上的感生电流0.001 a.u.等值线，曲线内部区域是对感生电流积分的区域，红色和蓝色分别对应paratropic和diatropic部分。可见diatropic分布在环的外侧，这与之前从向量场图中看到的顺时针环电流的主要分布是一致的。图中的红色小点对应截面上感生电流的极值点。键上面标的文字99对应于前述的参考比值0.99乘以100。文字出现的位置在键的正中央，不容易分辨，可以在图形窗口激活状态下按W键关闭文字显示，然后再自行ps上。 如果你想一次性把所有键的电流强度都计算，又不想一次次麻烦地输入y，那么可以写一个文本文件比如叫ally.txt，里面有12行，每一行都是y。运行BOCUST < ally.txt，这样所有键就都会被计算了。此时对得到的BCS.3d绘制会看到下图 下面介绍一下BOCUST运行时可以接的选项中比较重要的： -o、-f、-B、-m、-q等：和前面介绍的JBMAP程序的情况完全相同。即BOCUST也可以用-B指定外磁场方向、用-q指定考虑哪些轨道，等等。 -d、-D：确定哪些原子对儿被考虑时用的距离下限和上限（埃），默认-d 0.9 -D 1.6 -e：确定键截面上做积分的区域用的感生电流等值线的数值（a.u.），默认-e 0.001 -CS0：参考电流强度（nA/T），默认-CS0 12 -nocd：不显示积分区域和极值点。如果想让图上只有箭头而没有曲线和小点，应该加上这个 -l [比值]：在截面上绘制感生电流的不同数值的等值线，后面跟的是相邻等值线数值的比值 -C [粗细]：绘制截面的边框 例如这里计算C1-C2的sigma电子的键电流强度（即扣除三个pi轨道贡献），并在截面上显示等值线、显示边框，就运行BOCUST -q -3 17 20 21 -l 3 -C 5，然后输入y回车，再N回车。之后用v3d对BCS.3d作图看到下图。可见整个截面，以及里面的感生电流等值线，都显示得很清楚。由于这个截面上diatropic和paratropic的大小都差不多（毕竟苯分子没有sigma芳香性特征），所以键电流强度仅为0.007 a.u.。 5 总结 此文对分析磁感生电流的程序SYSMOIC做了介绍并给出了具体例子。由例子可见，通过SYSMOIC结合Gaussian绘制感生电流向量场图很方便，不需要借助第三方程序。SYSMOIC做键电流强度计算也很自动化，不需要像GIMIC那样谨慎地定义积分区域，而且键电流大小和方向还能直接通过箭头直观显示出来。SYSMOIC的计算速度也比较快。 本文只涉及了SYSMOIC最常用的功能，还另外有一些其它功能，比如对感生电流密度做拓扑分析，计算各种感生电流密度相关的场（ACID及其各向异性、磁屏蔽密度、vorticity张量的迹及其各向异性等）并绘制成等值面和填色平面图、计算任意格点上的感生电流张量和导数等等。 《详谈分子轨道对磁感生电流的贡献》（//www.umsyar.com/703）和本文有密切关联，非常建议阅读。 //www.umsyar.com/699 2024-02-09T17:09:00+08:00 谈谈有小虚频时热力学量的计算 On the calculation of thermodynamic quantities when there are small imaginary frequencies 文/Sobereva@北京科音  2024-Feb-9 0 前言 研究中在优化完极小点结构然后做振动分析时遇到虚频是家常便饭的烦人的问题。之前我专门写过博文《Gaussian中几何优化收敛后Freq时出现NO或虚频的原因和解决方法》（//www.umsyar.com/278）讲了怎么消虚频。除非程序存在bug，否则原理上不存在按上文处理完全解决不掉的虚频。一个现实问题是：虽然靠各种尝试、折腾，最终确实总能消掉虚频，但这个过程对于大体系来说过于浪费计算时间，也耗费研究者的精力。那么是否带有虚频的结构就一定不能拿来算能量相关的问题呢？是否虚频总是要不惜一切代价消除呢？此文就谈谈我个人的看法，给出理论分析，并且通过一个例子对我的观点予以验证。 本文牵扯到热力学数据计算的各种常识知识，没有这些常识的话先看《使用Shermo结合 程序方便地计算分子的各种热力学数据》（//www.umsyar.com/552），以及Shermo程序手册的附录，里面有计算原理的简洁易懂的概述。在北京科音基础（中级） 培训班（http://www.keinsci.com/workshop/KBQC_content.html）里对热力学量的计算有非常全面深入的介绍，欢迎关注和参加。 本文说的小虚频和过渡态必然存在的虚频没有丝毫联系，初学者绝对不要把本文内容胡乱套用到过渡态的情况。 1 导致小虚频的原因 为了避免误解，首先回顾一下优化极小点后振动分析时出现虚频的各种常见原因： (1)振动分析时用的影响势能面的设置与优化时的不同。这通常是没有基本常识的外行人才会犯的低级错误。注：JCTC, 17, 1701 (2021)中提出的single-point Hessian (SPH)方法通过对势能面添加偏移势，从而允许振动分析用的计算级别低于优化用的级别来节约时间，但这个做法非主流、经验性强，本文不涉及 (2)初始结构的点群高于实际极小点的点群，且如果优化后的结构维持住了初始的对称性，振动分析后就会发现破坏对称性的虚频，即按照虚频方式移动结构会导致结构的对称性下降。有时候这种有虚频的结构恰对应于过渡态 (3)程序bug。比如Gaussian 09 D.01的DFT-D3色散校正对解析Hessian的贡献的计算有bug，容易导致算势能面平缓体系出现虚频。这点在《DFT-D色散校正的使用》（//www.umsyar.com/210）我专门做了说明 (4)某些方法对势能面引入明显数值噪音。如Gaussian中的SMD溶剂模型的这个问题就比较明显 (5)几何优化收敛得不精确。几何优化收敛限过于宽松或Hessian的质量不太好会导致这种问题。 (6)积分格点质量不够好。这个问题对于M06-2X等对积分格点敏感的泛函尤为显著。关于 程序用的原子中心的积分格点参考《密度泛函计算中的格点积分方法》（//www.umsyar.com/69）。如果是CP2K等利用平面波的第一性原理程序，用的均匀分布的格点依赖于平面波截断能的设置。 以上前四种情况导致的虚频是一定要解决掉的，而且本来解决就非常简单。只要不犯低级错误自然就没有(1)的情况出现，解决(2)按照虚频模式调节结构后进一步优化即可，解决(3)就换成没bug的版本或其它计算程序，解决(4)可以换成其它的能在优化+振动分析中实现基本相同目的而没严重数值噪音的方法，如SMD换成IEFPCM，或者用SMD时结合SAS方式定义溶质孔洞表面。对于(5)和(6)情况带来的虚频，只要计算能力和精力允许，能消尽量要消，但如果由于计算资源或精力所限实在太难消掉、被迫要容忍虚频的存在，那么后续的能量相关问题的计算就是后文讨论的内容了。 对于柔性大分子、弱相互作用复合物这样势能面的某些维度非常平缓的体系，在做完几何优化和振动分析后很容易发现存在大小非常小的虚频，如波数为9.8i、27.3i cm^-1这种，普遍小于50，这一般是前述(5)或(6)情况所导致的。此类虚频的振动模式通常对应于甲基旋转运动、柔性骨架的变形运动、分子间相对运动。这类(5)或(6)情况所导致的虚频以下简称为“小虚频”。但要注意也不是所有波数小于50的虚频都属于这种情况，前述的其它四种情况也完全可能导致数值很小的虚频，所以绝对别光拿虚频大小来判断情况。 2 有小虚频时能量相关问题的计算 当遇到小虚频时，显然最完美的做法（具体选择视情况而定）是用诸如更严格的几何优化收敛限、更好的Hessian、更好的积分格点精度进一步做优化和振动分析，使得虚频完全消除后再计算能量相关问题。但是对于普通泛函算几百个原子有机体系，很可能算一次Hessian就要花一整天时间，几何优化一步也可能花费一个小时左右，最终把虚频完全消掉可能得花几天甚至一个礼拜，十分痛苦，甚至不可能实现。那么，在后续算能量相关问题时怎么考虑？下面根据被计算的量专门讨论。 2.1 计算电子能量 诸如Gaussian这样的主流 程序默认的几何优化收敛限是以受力和位移作为判据的。在默认收敛限下优化正常完成后，即便结构还没有与实际极小点特别精确一致，但通常相差很小。而且当前的最大和平均受力肯定已经很小了（因为收敛限考虑了它们），由于受力对应于能量梯度的负矢量，所以当前结构稍微再挪一点达到精确的极小点的过程中能量变化肯定非常小。即曰，在有小虚频的情况下，计算的电子能量一般是可以用的，也即结果和消掉虚频后再算的能量相差极小。 2.2 计算0 K下的内能U(0) U(0)是电子能量与零点振动能（ZPE）之和。在谐振近似模型下，每个振动模式对ZPE的贡献是1/2*h*ν，h是普朗克常数，ν是振动频率。频率很低的振动模式对ZPE的贡献是非常小的，每波数的振动频率对ZPE的贡献仅为5.98 J/mol，因此诸如30波数的振动模式只对ZPE有0.18 kJ/mol（0.04 kcal/mol）的贡献，数量级远远低于理论方法、基组、溶剂模型等因素带来的误差。顺带一提，Shermo计算热力学量的时候如果把settings.ini里的prtvib设为1，每个振动模式对热力学量的贡献都会明确输出。 主流 程序，以及Shermo程序在默认情况下，在计算振动对热力学量的贡献时都是不考虑虚频产生的贡献的，因此小虚频对热力学量的贡献会被直接忽略掉，这不会对ZPE的计算造成实际问题。因为如果把这个小虚频消了，从而产生一个实频，其频率肯定也是非常低的，对ZPE的贡献可忽略不计。也就是不管消不消这个虚频，算出的ZPE都没多大差别。不过，实际情况远没这么简单，因为在虚频消掉前后，由于二者所在势能面的位置的些许不同，或改用不同档次积分格点导致频率计算精度的差异，其它的实频模式也会多多少少受到影响，不过频率值受影响相对明显的主要是低频模式（可能影响几个甚至10个左右波数，而高频模式受影响一般小于2波数），它们对ZPE的贡献恰恰又很小，所以虚频消掉前后算的ZPE虽然未必差异小到诸如前面说的0.2 kJ/mol左右这种程度，但至少也不会大，比如对柔性大体系顶多也就影响个1-2 kJ/mol。 如上所述，鉴于是否小虚频对电子能量和ZPE的影响都不太明显，因此算电子能量和U(0)时一般可以容忍一个或多个小虚频的存在。 2.3 计算从0 K升温到当前温度对内能的影响 研究U(0)没太大实际意义，毕竟实际化学过程也不可能在0 K发生。计算有限温度，也就是高于0 K时候的热力学量，就需要考虑从0 K升温到当前温度对内能的改变量，其中振动产生的贡献下面用Uvib(0→T)表示。首先要知道一点，低频对Uvib(0→T)的贡献是明显高于高频模式的。这点从下面的J. Phys. Chem., 100, 16502 (1996)的图1可以明显看出来，图中ΔHvib(T)和这里说的Uvib(0→T)是一码事，图是对T=298.15K的情况绘制的。   虽然低频模式的贡献大，但是由图中放大的0-100波数的区间可见，在小几十波数范围以内，频率值对Uvib(0→T)的影响虽然不可忽视，但不算多大，都在2.0-2.4 kJ/mol范围内。这也就是说，由于存在一个小虚频导致少考虑一个实低频模式，会给Uvib(0→298.15K)带来-2.2 kJ/mol左右的误差，这虽然不算小，但还算可以凑合接受。如前所述消虚频后还可能会对其它低频模式产生不可忽视的影响，这对Uvib(0→T)的改变一般也不至于太大。 综上所述，若带着小虚频计算U(T)，会造成一定低估。对于有一个小虚频的情况，低估程度在0.5 kcal/mol左右。如果有多个小虚频，那还是别计算U(T)了，除非T远低于常温，或者结合后文说的把小虚频当成实频的做法。 2.4 计算熵、自由能 低频对熵的贡献远大于高频，而且在很低频区间内，频率越低贡献明显越大。在包括Gaussian在内，大多数 程序（ORCA除外）计算热力学量用的是传统的rigid-rotor harmonic oscillator (RRHO)模型，振动熵的计算完全基于谐振近似模型。此时在低频区间内，随着频率趋近于0，振动模式对熵的贡献猛增，如下面的Chem. Eur. J., 18, 9955 (2012)文中图2的虚线所示。Grimme在这篇文章中还提出了一种quasi-RRHO模型（QRRHO，后来也叫modified RRHO、mRRHO），以下称为QRRHO(Grimme)，它将自由转子和RRHO用的谐振近似模型算的振动熵做插值，频率越接近0自由转子的权重越大、频率越高谐振近似的权重越大，如下图实线所示。QRRHO(Grimme)对高频模式算的熵贡献和RRHO完全一样，而对于低频模式（特别是波数在100以下的）考虑得比RRHO明显更合理。由下图还可见QRRHO(Grimme)算的低频模式的熵比RRHO整体小得多，没有RRHO的严重高估的问题。此外，QRRHO(Grimme)下低频模式的振动熵对频率值的敏感性低于RRHO很多。   直接计算熵的用处不大，算熵的主要目的是用来算自由能，因为自由能里有-T*S项。显然低频对熵的贡献计算的准确度显著影响到自由能的计算准确度。假设体系有一个小虚频，并且在消掉之后会多一个20 cm^-1的实频，并且用的是QRRHO(Grimme)算298.15 K的熵，那么在有虚频的结构下由于缺失了这个实频，根据上图可知会造成熵的误差约为-18 J/mol/K，对自由能造成的误差约为-298.15*(-18)/1000=5.4 kJ/mol，略高于1 kcal/mol，这就不小了，都能赶上高质量ωB97M-V泛函结合大基组算普通有机反应的误差了。不过实际中小虚频带来的误差未必那么大，因为有和没有小虚频的结构下实低频模式的频率也往往有不可忽视的不同，例如下面第3节的例子在有虚频的结构下实低频模式的频率比精确极小点结构下的整体偏小（势能面形状原因所致），导致高估了不少实低频模式的熵，藉由-T*S项进而对自由能造成低估，于是和少考虑一个低实频对自由能的高估产生了很大程度的相互抵消，下面管这称为“熵的误差抵消”。 根据以上讨论，有小虚频的情况下计算熵和自由能都有一定风险。如果就有一个小虚频，算常温的情况往往还说得过去，而如果有多个小虚频，或计算高温的情况（注意T越大-T*S越大，熵的误差造成自由能的误差越大），我就不建议冒险带着虚频了。而且就算非要带着小虚频算，也至少应当用QRRHO模型，而切勿用RRHO模型，因为RRHO算的振动熵对低频频率值的敏感程度显著高于QRRHO，因而用有虚频的结构算的熵的误差上限比QRRHO的明显要大。 要注意QRRHO模型不止Grimme提出的这一种。如Shermo手册附录部分所介绍的，Truhlar的QRRHO是把所有小于一定阈值（通常为100 cm^-1，可以由Shermo的ravib参数控制）的低频统一提升为这个阈值再算它们对Uvib(0→T)和熵的贡献，这个做法的物理意义不甚明确，没有QRRHO(Grimme)的插值做法优雅。目前还有一种QRRHO是Minenkov等人在J. Comput. Chem., 44, 1807 (2023)中提出的，它对振动熵部分的考虑和QRRHO(Grimme)一样，但还同时将这种自由转子与谐振近似做插值的思想同时应用到了计算振动对内能的热校正量Ucorr上（此时没法单独计算ZPE），原理更理想，实际结果还略好一点点。这些QRRHO模型在最新版Shermo中都是支持的，通过ilowfreq参数控制。 2.5 将小虚频视为实频的做法 在Angew. Chem. Int. Ed. 2022, e202205735中的3.2节Grimme等人认为存在小虚频的情况下，只要使用QRRHO模型并且将小虚频视为实频处理，比如18i cm^-1直接当做18 cm^-1，则得到的熵就是可以用的、小虚频的存在就不再是个问题。这种做法表面看起来似乎有点莫名其妙、太任性了，但确实有一定道理，可以这么理解：在消除小虚频后，必定会得到相同数目的波数也很小的实频模式，它们都对熵有明显贡献。相对于完全忽略掉这些小虚频的贡献而导致明显低估熵，干脆把他们当做实频处理，从而“凑”出来相应数目的实低频是更好的做法，而且本身在QRRHO模型下，数值在十几、几十波数区间的频率对熵的贡献也都相差不悬殊，因此也用不着对这么人为凑的实频频率太较真。由于如前所述，Uvib(0→T)也对低频数值不敏感而低频的贡献又没小到可以随便忽略的程度，所以计算Uvib(0→298.15K)时也使用这种处理是有道理的。从Shermo 2.5版开始，可以通过settings.ini里的imagreal选项设置视为实频的虚频阈值，比如设为50，就代表大小在50 cm^-1以内的虚频都当做实频处理来计算振动对各种热力学量的贡献，此时从Shermo程序输出的频率信息中也会看到已经没有虚频了。 然而，这种将小虚频视为实频的做法目前很非主流，其思想在我来看过于主观和乐观、把问题想得太简单了，目前也缺乏系统性的检验，我个人不建议轻易使用。这种做法有一个明显问题：如前所述，有小虚频和消掉虚频的两种情况相差的不仅仅是一个小虚频转变成一个小实频，许多其它低频模式的频率也会有不小变化，而且有时是整体明显增加的。不做这个虚频到实频的转换的情况下能够有上一节说的熵的误差抵消效应，而这么做了之后，忽略虚频模式导致对熵的低估问题确实很大程度解决了，但在有小虚频结构下其它频率往往偏低而导致熵的高估带来的误差就没处抵消了。此外，如果当前体系的低频还有很多都是10 cm^-1 左右这种特别低的，这个范围附近即便是QRRHO(Grimme)或QRRHO(Minenkov)模型算的熵也对频率很敏感，问题会更大。 所以，将小虚频视为实频的做法虽然确实对一些情况更好地计算熵和自由能是有帮助的，但也很可能起到负面效果。所以这绝不是普适、稳健、理想的做法，更千万不能一看到存在小虚频就总是想靠这个做法来完全无视虚频。 3 有虚频时算热力学量的实例 笔者最近通过Gaussian 16在ωB97XD/6-311G*级别研究的一个以pi-pi堆积方式形成的共128个原子的三分子复合物，在优化和振动分析后发现存在6.44i cm^-1的小虚频，之后进一步优化时用opt=calcfc关键词提供精确的初始Hessian，优化完再做振动分析就发现没虚频了，此时最低的频率成了10.35 cm^-1。正好这个例子可以拿来考察消虚频前后计算的热力学量的差异，由此检验带着那个小虚频时以不同方式算的热力学量有多大误差。要强调的是，仅仅这么一个测试，显然不足以充分证明某种做法理想，因为这必须通过大量体系做统计分析，但至少误差较大的情况可以用来体现相应做法不够稳健和普适。 先看一下有虚频相对于无虚频时各个实频的频率变化情况。先把实频频率按照由低到高排列，用有虚频时候它们的频率减去无虚频时候的频率作为纵坐标值，而两种情况频率的平均值作为横坐标，作的图如下所示。可见，至少对于此例，有虚频时数值较低的实频是普遍明显偏低不少的，而>200 cm^-1的相对而言的中、高频则受到的影响很小。这种情形虽然不能说是普遍情况，但至少也是挺常见的一种。 此例有虚频相对于无虚频的情况计算的电子能量仅仅相差0.04 kcal/mol，小到完全可以忽略不计。用Shermo程序基于Gaussian振动分析输出文件计算的有虚频相对于无虚频时的ZPE和气态标况下的熵、内能、自由能热校正量的差异如以下表格所示。RRHO以及Truhlar、Grimme、Minenkov三种QRRHO分别对应Shermo的ilowfreq=0/1/2/3的情况，imagreal=0和50相当于忽略那个虚频以及把那个虚频当大小相同的实频考虑的情况。 首先看常规的忽略虚频贡献的情况。由表格可见： • ZPE：有虚频时和无虚频时ZPE相差仅-0.3 kcal/mol，属于通常可以接受的差异 • 熵：由于误差抵消的巧合，恰好RRHO下有小虚频时的熵和无虚频时的差异很小。QRRHO(Truhlar)的情况差异挺大，在于它把实低频都统一当成了100 cm^-1，这部分在有和没有小虚频情况之间几乎都抵消了，而有小虚频时相当于少考虑了一个100 cm^-1振动模式对熵的贡献，因此造成熵的严重低估。QRRHO(Grimme)和QRRHO(Minenkov)对熵的部分处理相同，由于前述的熵的误差抵消效应导致有和没有小虚频时熵的差异不算太大 • 内能的热校正量Ucorr：等于ZPE+Uvib(0→T)。由于Uvib(0→T)对低频的具体数值不太敏感，有无小虚频情况下实低频的贡献变化很小，而有小虚频时由于缺失了一个低实频对Uvib(0→T)的贡献，导致有小虚频时Uvib(0→T)偏低不少，Ucorr因此比ZPE偏低得明显更多。QRRHO(Minenkov)用的计算Ucorr的方式和其它模型不同，对于此例使得有小虚频时Ucorr偏低程度比其它模型更小，仅为-0.5 kcal/mol。因此在计算有小虚频情况下的U(T)时，用QRRHO(Minenkov)可能会比其它模型更有优势。 • 自由能的热校正量Gcorr：相对于其它模型，有小虚频相对于无虚频时Gcorr的差异在RRHO下是最大的，这在于有无小虚频时RRHO算的S的差异对此例恰好颇小，导致-T*S项的差异才0.03 kcal/mol，因此没能充分对Ucorr的差异产生抵消效应。而各种QRRHO模型下Gcorr的差异明显较小，这在于这种情况下S的差异不小，而且-T*S项的差异的符号和Ucorr的差异相反，因此相互抵消了很多。其中Gcorr差异最小的是QRRHO(Minenkov)，才-0.11 kcal/mol，体现出它可能适合作为有小虚频情况下算自由能的优先选择。 再来看将小虚频当实频处理的情况。根据以上表格可见，这种处理使得Ucorr在有和没有小虚频情况下计算的结果差异减小很多，因此对于算内能的目的，这种做法有一定意义。但是，除了QRRHO(Truhlar)外，这个做法却造成了有小虚频时算的S和Gcorr与无虚频时的差异变得显著更大，对RRHO最为严重，对QRRHO(Grimme)和QRRHO(Minenkov)问题也挺严重，这来自于有小虚频时实频模式偏低造成的对振动熵的高估无法与因为少考虑一个实低频导致的低估一定程度相抵消。把小虚频当实频处理时结合QRRHO(Truhlar)对此例倒是很不错，有小虚频相对于无虚频时的熵只相差0.27 cal/mol/K，这在于此种QRRHO将低频全都上拉到100 cm^-1，此时又把小虚频转化成了小低频，因此无虚频和有虚频时相当于有同等数目的100 cm^-1的实频，自然两种情况下算的熵很接近。 4 总结 根据以上讨论和实际体系的测试可知，对于因(5)和(6)情况造成有小虚频的情况，算电子能量和U(0)没大问题，不是非得消虚频。 对于计算Ucorr或与之相差RT的焓的热校正量Hcorr，用QRRHO(Minenkov)或QRRHO(Grimme)并同时把小虚频当实频处理是较好选择，不会因为存在小虚频的问题而带来太大误差。 有小虚频时最难搞的是熵，以及温度不很低时候Gcorr的计算，难点在于小虚频消掉前后实低频模式频率也可能受到不小影响且影响情况不易估计，然而即便用QRRHO(Grimme或Minenkov)时低频的具体数值对熵的影响也不算太小，所以明显不是把小虚频当成实频处理以避免少考虑一个实低频模式就能较好解决的，而且如第3节的测试可见，这种做法还可能造成熵和Gcorr在有小虚频时计算结果更差。所以当你对熵、Gcorr要求精度高时，小虚频要尽量消。实在搞不定的话，那就姑且用原理上最好的QRRHO(Minenkov)，并且建议不用小虚频当实频的处理，毕竟此做法非主流而且起到正面效果的可能性很有限。虽然把小虚频当实频并结合QRRHO(Truhlar)的情况下小虚频和无虚频时熵的差异很小，但在J, Comput. Chem., 44, 1807 (2023)的测试中QRRHO(Truhlar)算团簇反应自由能表现得远不如QRRHO(Minenkov)，甚至比RRHO没明显改进，所以我也不推荐。 //www.umsyar.com/698 2024-01-29T10:07:00+08:00 谈谈如何计算环张力能：以CPP和碳单环体系为例 On the calculation of ring strain energy: Taking CPP and cyclocarbon systems as examples 文/Sobereva@北京科音   2024-Jan-29 0 前言 环状化学体系具有不同程度的环张力，环张力和环状体系的诸多问题紧密相关，诸如生成焓、稳定性、合成难易程度等。环张力能（strain energy, SE）是衡量环张力大小的最直接的指标，它对应当前环状结构转变为假想的没有环张力的结构过程中能量变化的负值。 北京科音自然科学研究中心（www.keinsci.com）的卢天等人之前受邀发表了一篇专门精确考察碳环（cyclocarbon）体系环张力能的文章，文中对环张力的计算思想有详细介绍并给出了诸多讨论： Accurate theoretical evaluation of strain energy of all-carboatomic ring (cyclo[2n]carbon), boron nitride ring, and cyclic polyacetylene, Chin. Phys. B, 31, 126101 (2022) DOI: 10.1088/1674-1056/ac873a（不方便获取者也可以看预印版https://doi.org/10.26434/chemrxiv-2022-v8w9h，内容基本一样） 下面，将基于此论文中的内容简要说明环张力能的计算原理并给出具体例子。感兴趣的读者请务必仔细阅读论文原文以了解更多信息，也推荐以类似方式计算其它体系的环张力能时引用此论文。另外，读者如果对碳环体系感兴趣，很建议查看同作者之前发表的此类体系的各种研究文章和相关博文：//www.umsyar.com/carbon_ring.html。 1 同联反应法 设计同联（homodesmotic）反应是计算环张力能的很常用方法。这个方法设计一个化学反应，既使得环被打开使得环张力被完全释放，同时又尽可能不影响原本的成键情况，显然这个反应能的负值就对应了环张力能。 例如，要计算氧杂环丁烷的环张力能，就可以设计以下同联反应。红色部分对应具有显著环张力的氧杂环丁烷，将它插入到甲乙醚所示位置就变成了右边的产物，显然此时环张力就被完全释放掉了。与此同时，反应物和产物的成键特征并未发生变化，例如下图红字的第一个碳在反应前后始终是sp3杂化状态、始终连着一个氧原子和一个CH2，因此化学环境变化导致的成键特征的极轻微改变对这个反应能的影响可以忽略不计，也即这个反应能几乎只体现出环张力的释放。 同联反应法计算环张力能的缺点是反应的设计并不唯一，比如以上反应的甲乙醚也可以是甲醚、乙醚等等，结果或多或少会有点差异。显然反应应当设计得尽可能令各个键的特征在反应过程中不发生改变，但这并非总能很好实现，后文提到的碳环体系就是一例。 2 同联反应法计算实例：[6]CPP 这一节示例使用同联反应法计算[6]CPP的环张力能。此体系的几何结构，以及设计的同联反应如下所示。可见是把[6]CPP插入到了联苯中间变成了直线的八联苯，从而完全释放了环张力。 上图涉及的三个分子都在ωB97XD/6-31G*级别下用Gaussian 16程序进行了几何优化，并做振动分析确认了无虚频，输入输出文件都在//www.umsyar.com/attach/698/file.rar里。基于输出文件里最后一次输出的电子能量，用反应物能量之和减去产物能量，就得到了环张力能：627.51*(-463.1444911-1385.7232492+1849.0260518) = 99.3 kcal/mol。这个值比C-C单键的键能还要大一些，可见其环张力是相当强的。 如果想得到更精确结果，可以基于优化完的结构在更好的级别下算高精度单点再求差，比如基组可以改用明显更好的def2-TZVP。这对结果可能造成几kcal/mol程度的影响。 在《18个氮原子组成的环状分子长什么样？一篇文章全面揭示18氮环的特征！》（//www.umsyar.com/725）介绍的笔者的ChemPhysChem, 25, e202400377 (2024)研究文章中利用同联反应法计算了18氮环的环张力，推荐感兴趣的读者阅读了解细节。 3 外推法 如果被考察的环状体系是由重复单元构成的，比如前面的[6]CPP是苯环作为重复单元，就也可以用这一节介绍的外推法计算其环张力能。这个方法计算过程如下 (1)计算从小到大不同重复单元数（n）的环的每单元的能量E/n，这里E是整体能量。显然E/n是依赖于n的，n较小时环尺寸较小，自然环张力大，平均到每个单元的环张力能也大 (2)用Origin之类程序通过恰当的函数（一般是二次函数）拟合E/n与n之间的解析关系，并外推得到n无穷大时的E/n，这对应于没有环张力时候的每单元能量 (3)将当前考察的环状体系的E/n减去n无穷大时的E/n，就得到了当前体系的每单元环张力能（SE/n） (4)将当前体系的SE/n乘以当前的体系的n即是当前体系的环张力能 相较于同联反应法，外推法的优点： •原理非常严格 •没有同联反应法设计反应的任意性，并可以用于同联反应法不适合的体系 •可以顺带得到环张力能与n的解析关系，由此可以直接预测出任意重复单元数的环张力能 缺点： •需要计算从小到大不同尺寸的环状体系，总耗时明显高于同联反应法 •计算步骤略繁琐，得自己做拟合 如果你只是要考察当前一个环状体系的环张力能，而且可以设计出很合适的同联反应，那么就没必要用昂贵且略麻烦的外推法。 4 外推法计算实例：18碳环 2019年首次在凝聚相实验中观测到的18碳环分子由18个sp杂化的碳依次相连构成环状。由于sp杂化的碳明显倾向于形成180度键角，因此18碳环必然有显著的环张力，具体环张力能是多少非常值得研究，这正是前述Chin. Phys. B, 31, 126101 (2022)文章研究的主要内容。 此文中指出碳环类体系不适合同联反应法计算环张力能。碳环类体系是长、短C-C键交替构成，两个碳是一个重复单元，若用同联反应法计算化学式为C_2n的碳环的环张力能，可以设计如下反应式，相当于把碳环的每个重复单元插入到了一定长度的碳炔里。碳炔是sp杂化的碳形成的链状分子，两端由氢封端。PS：对碳炔感兴趣者建议阅读《氢封端碳链H-(C≡C)n-H (n = 3-9, 15)的电子光谱的尺寸依赖性：性质分析及对碳炔的预测》（//www.umsyar.com/679） 然而，如下图可见，将18碳环的重复单元插入到不同长度（通过m体现）的碳炔里，按照同联反应算出来的环张力能有明显区别，因此靠同联反应法难以得到严格的18碳环的环张力能。 碳炔和18碳环都具有C-C键长、短交替的特征。文中还发现，不管碳炔的m取多少，处于它中央的较长的C-C键和18碳环中的较长的C-C键之间的键长差异都不可忽视，对较短的C-C键也是如此。因此原理上没法设计一个同联反应使得18碳环的C-C键特征在转移到产物时几乎完全不发生改变，这也必然导致靠同联反应法不可能对此体系得到完全精确的环张力能，也即成键特征差异造成的能量改变会不可避免地掺入到计算出的环张力能中。 由于发现了以上问题，文中使用外推法计算碳环的环张力能。为了追求尽可能好的精度，文中用ORCA程序在精度很好的DLPNO-CCSD(T)/cc-pVTZ级别下对n=7到n=24的碳环都计算了单点能（基于ωB97XD/def2-TZVP优化的结构），ORCA输出文件在//www.umsyar.com/attach/698/file.rar里都提供了。每重复单元能量（E/n）与n的关系如下所示。由于n较小的碳环的电子结构特征有特殊性、和较大的碳环存在一定差异，因此文中只对n>=13的数据点做了二次函数的拟合，由下图的红色曲线可见拟合质量超极好，R^2近乎精确为1。 令拟合公式E/n = 0.8844/n^2 - 76.00997中的n为无穷大，就得到了无环张力对应的无穷大的碳环的每单元能量-76.00997 Hartree。每单元环张力能SE/n = E/n - 76.00997，于是有SE/n = 0.8844/n^2，也即碳环类体系的环张力能SE = 0.8844/n Hartree。18碳环的n=9，代进去就得到了其环张力能0.8844/9 Hartree = 61.7 kcal/mol。 值得一提的是，如果令前述同联反应法中碳炔的m设得很大，从而使得环张力能随m收敛，最终算出来的结果是59 kcal/mol，可见和严格的外推法的结果相比误差并不可忽略。 前述论文中也测试了用ωB97XD/def2-TZVP的能量结合外推法算碳环的环张力能，结果和上面用的昂贵的DLPNO-CCSD(T)/cc-pVTZ能量的结果没多大差别。这在一定程度上体现出ωB97XD/def2-TZVP级别算碳环类体系很不错。 文中还基于不同尺寸碳环的环张力能和键角之间的关系，推导出了适合碳环体系用的C-C-C键角力常数56.23 kcal/mol/rad^2。 此外，文中还对以下所示的18碳环的等电子体硼-氮环，以及具有单套pi共轭特征的环聚乙炔的环张力做了计算，发现硼-氮环类体系的环张力明显小于碳环类，而环聚乙炔类体系的环张力则略大于碳环类。文中利用Multiwfn程序（//www.umsyar.com/multiwfn）做波函数分析从电子结构角度清晰解释了它们的环张力与碳环之间存在差异的原因。文中还发现用常用的B3LYP泛函算全局pi共轭的环状聚合物的合理性较差，E/n随n的变化明显不合理。这些方面这里就不多提了，请读者务必阅读论文原文。 5 总结 环张力是环状分子的重要特性。本文简明扼要地介绍了两种最常用的计算环状分子的环张力能的方法，并以[6]CPP和18碳环作为具体例子进行了演示。同联反应法比较方便省事，但对于一些全局pi共轭的环状体系，诸如18碳环，则更适合外推法，明显严格得多，同时外推法还很有助于探究环尺寸与环张力之间的关系。还应当注意并不是所有环状体系都能用这两种方法之一计算环张力能，比如苯分子，它存在显著的和其环状结构紧密相关的六中心键，因而环张力本来就是难以被定义的。 //www.umsyar.com/696 2024-01-26T09:04:00+08:00 18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响 The unique aromaticity of isoelectronic molecule of cyclo[18]carbon, B6N6C6: Revealing the key role of the carbon atom bridging boron and nitrogen on electron delocalization 文/Sobereva@北京科音   2024-Jan-26 0 前言 由18个碳原子连接形成的环状体系18碳环（cyclo[18]carbon）具有十分特殊的几何和电子结构，自从它于2019年首次在凝聚相中被观测到后，18碳环及其衍生物就得到了化学界的广泛关注，并陆续有大量理论研究工作发表。特别是在Carbon, 165, 468-475 (2020)中，北京科音自然科学研究中心（http://www.keinsci.com）的卢天和江苏科技大学的刘泽玉等人证实18碳环具有明显的双pi芳香性特征，这是18碳环非常与众不同的一点。 18碳环的一种等电子体是B9N9，由B和N原子交替相连构成环状。此体系早已被研究过，被证实几乎不具备像18碳环那样的明显的芳香性。最近，卢天、刘泽玉和西班牙多诺斯蒂亚国际物理中心的Mesías Orozco-Ic发现，由6个(硼-碳-氮)重复单元依次相连构成的18碳环的另一种等电子体B6C6N6具有和18碳环很接近的芳香性。为什么在无芳香性的B9N9中掺杂进去碳原子后，或者说让碳原子桥联硼和氮后，B9N9的芳香性就能得到巨大提升？这是非常有意思的问题。于是以上研究者对B6C6N6的特征做了全面深入的理论研究，详细分析讨论了其芳香性特征和本质，成果近期在美国化学会的知名期刊Inorganic Chemistry上发表，非常推荐阅读： Exploring the Aromaticity Differences of Isoelectronic Species of Cyclo[18]carbon (C18), B6C6N6, and B9N9: The Role of Carbon Atoms as Connecting Bridges, Inorg. Chem., 62, 19986 (2023) https://doi.org/10.1021/acs.inorgchem.3c02675 此研究的研究意义不仅在于讨论了B6C6N6本身，还在于揭示了碳原子桥联硼和氮原子对电子离域特征的关键性影响，这对于通过碳掺杂的方式对纯硼、氮构成的体系（如硼-氮纳米管、二维层状h-BN体系等）进行改性提供了重要启示。此研究通过强大的Multiwfn程序（//www.umsyar.com/multiwfn）利用了各种波函数分析手段对B6C6N6进行了充分的考察并和18碳环、B9N9进行了对比，这也是充分运用波函数分析探究新颖体系的很好的范例。 下面本文就对Inorg. Chem., 62, 19986 (2023)这篇文章的主要研究思想、结论进行深入浅出的介绍，并对一些研究细节进行附加说明，以便于读者更好地理解此文，以及将此文的研究手段用于自己的研究。同作者之前还对18碳环及其衍生物的各方面特征做过充分的理论研究并得到了同行的广泛关注，成果汇总见//www.umsyar.com/carbon_ring.html。 1 B6N6C6的几何结构 做B6N6C6的各方面研究之前先得得到其可靠的几何结构。18碳环及衍生物用ωB97XD/def2-TZVP级别做几何优化是很稳妥的，在Carbon, 165, 468-475 (2020)以及//www.umsyar.com/carbon_ring.html里列举的此类体系的大量研究工作中都已经充分证明了这点。因此此文用Gaussian 16在此级别下优化了B6N6C6的结构。特别需要注意的一点是，B6N6C6基态是单重态，但是其闭壳层单重态波函数存在不稳定性，用DFT计算时需要作为对称破缺单重态来处理方可得到真正的基态，相关讨论见《谈谈片段组合波函数与自旋极化单重态》（//www.umsyar.com/82）。因此，此文对B6N6C6基态的各种研究都是基于对称破缺单重态的结构和波函数做的。如果误当成了闭壳层单重态来算，作者发现得到的几何结构将明显不合理，而且电子结构也明显错误，比如芳香性严重偏高。实际上之前有一篇论文也算过B6N6C6，但由于那些作者没注意这一点而误当成了闭壳层计算，导致结果是没有任何意义的。此例体现出DFT计算B6N6C6这种电子结构不同寻常的体系一定要注意做波函数稳定性测试。而ωB97XD/def2-TZVP级别下的B9N9、18碳环的稳定波函数都是单重态闭壳层的。 ωB97XD/def2-TZVP级别下优化得到的无虚频的B6N6C6、B9N9和18碳环的结构如下所示，红/灰/蓝分别对应硼/碳/氮，笛卡尔坐标在文章的SI里提供了。可见B6N6C6属于C6h点群，而B9N9属于D9h点群。 毕竟B6N6C6是一个新颖的结构，不排除静态相关很显著导致ωB97XD描述不理想的可能，故为了绝对稳妥，此文还用ORCA在CASSCF(12,12)/def2-TZVP级别下也做了几何优化。此活性空间设定下，最终会有12个占据数明显偏离整数的pi型自然轨道，而若活性空间进一步扩大到CAS(14,14)，则会再多出两个占据数接近整数的pi型自然轨道，这体现CAS(12,12)恰好足够充分描述B6N6C6的pi电子的静态相关效应。CASSCF(12,12)优化出的结构和ωB97XD的非常相近，体现出ωB97XD得到的此体系的几何结构靠谱。为了进一步体现ωB97XD计算出的波函数也合理，文中也顺带对比了一下这两个级别得到的单电子分布情况。文章按照《谈谈自旋密度、自旋布居以及在Multiwfn中的绘制和计算》（//www.umsyar.com/353）计算了ωB97XD波函数的自旋密度，而由于CASSCF没法得到自旋密度，故文章使用《使用Multiwfn计算odd electron density考察激发态单电子分布》（//www.umsyar.com/583）介绍的方法对CASSCF(12,12)波函数产生了单电子密度（ODD），等值面图对比如下所示。在不考虑自旋密度符号（体现在等值面颜色上，不同颜色对应单电子不同自旋方向）的情况下，可见二者展现出的单电子分布是基本吻合的，即单电子主要都分布在碳上。这在很大程度上体现出ωB97XD对称破缺计算产生的波函数也是合理的，适合用于之后的波函数分析。 2 B6N6C6的电荷分布 原子电荷是定量考察化学体系中电荷分布最常用的指标之一。18碳环由于所有的碳都是等价的，显然原子电荷都为0，而对B6N6C6和B9N9计算原子电荷，可以直接考察这俩体系中不同原子带的净电量、了解二者电荷分布的差异。原子电荷计算方法并不唯一，不同方法基于不同的思想、有不同的特点，在《原子电荷计算方法的对比》（DOI: 10.3866/PKU.WHXB2012281）中有大量介绍和对比。考虑到B6N6C6电子结构的特殊性，为了得到稳妥的结论，文中用Multiwfn计算了很常用的ADCH、Hirshfeld、Hirshfeld-I、CHELPG原子电荷（它们的计算例子看Multiwfn手册4.7节）、用NBO程序计算了常见的NPA原子电荷，结果如下所示。可见虽然不同原子电荷给出的具体数值不同，但可以得到共同的结论，就是氮明显带负电，硼明显带正电，这体现出了它们显著的电负性差异。而在B6N6C6中，硼和氮的净电荷差异明显小于它们在B9N9中，这体现出碳元素在硼、氮之间的插入明显促进了体系整体电荷分布的均衡性。而B6N6C6中碳自己的原子电荷则接近于0，即没怎么受周围化学环境的影响。 18碳环具有四类分子轨道：平面内和平面外两套π轨道，分别称为π-in和π-out，此外还有σ轨道，以及碳的内核原子轨道组成的内核分子轨道。B6C6N6和B9N9也有这四类轨道。为了深入考察各个原子上的电子是怎么分布在这些轨道上的，文中用Multiwfn基于Hirshfeld-I原子空间划分对三个体系计算了电子在这些轨道上的布居数，如下所示。具体计算方法我专门写了篇文章说明：《使用Multiwfn基于Hirshfeld-I划分计算特定类型电子在各个原子上的分布量》（//www.umsyar.com/697）。 从以上数据可以看到，在每个原子上，π-in和π-out电子分布量都是基本一致的，没有什么偏向性，体现出这些体系里π-in和π-out电子特征的近似等价性。B6C6N6和18碳环中的碳的电子结构特征非常相近，π-in和π-out电子数在两种体系里都基本为1.0，明显都处于sp杂化的状态。 3 B6N6C6的成键情况 键级是考察化学键特征非常常见的方式，参见《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）。文中用Multiwfn对B6C6N6、B9N9和18碳环计算了流行的Mayer键级，并且分解为了不同类型轨道的贡献，如下所示。Mayer键级体现了原子间等效的共享电子对数，由数据可见，三个体系里所有键都包含典型的单个σ键特征（即σ+core电子贡献的Mayer键级接近非常接近于1），同时π-in和π-out电子都构成了一定程度的π键特征，且二者贡献相近。B6C6N6中的π作用程度是N-B ≈ C-N > B-C。B9N9中所有B-N等价。18碳环由长、短两种C-C键构成，后者的π作用显著强于前者。文中还用Multiwfn计算了模糊键级以进一步确认了这些结论，数据见补充材料。 《使用IRI方法图形化考察化学体系中的化学键和弱相互作用》（//www.umsyar.com/598）中介绍了卢天提出的IRI-π函数，可以直观地展现化学体系中π电子相互作用情况。此文对B6C6N6和B9N9绘制了IRI-π=1.2的等值面图并用sign(λ2)ρ函数着色，如下所示。可见每个键上都出现了环状的等值面，非常直观地体现出这些键都存在明显的π作用，这和Mayer键级展现出的信息相吻合。 4 B6N6C6的芳香性 B6N6C6的芳香性是当前研究的重点。和18碳环一样，B6N6C6和B9N9的π-in和π-out电子皆为18个，都满足Hückel的4n+2芳香性规则，因此都有可能存在双π芳香性。衡量芳香性的方法非常多，在《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）里有充分的介绍。其中基于电子结构和基于磁性质的芳香性衡量方法的可靠性和普适性普遍较好，可以用于考察B6N6C6的芳香性并与18碳环和B9N9的进行对比。 AV1245*1000在《使用Multiwfn计算AV1245指数研究大环的芳香性》（//www.umsyar.com/519）里专门介绍过，是基于波函数定义的适合用于定量考察共轭大环体系芳香性的指数。文中首先使用Multiwfn计算了AV1245*1000，结果如下。除了总值外，还分别给出了π-out电子和其它电子的贡献（注：π-in、σ、core电子对AV1245的贡献无法精确分离故没有单独给出。但显然σ和core的离域性极低，不会对AV1245有什么贡献，故“其它电子贡献”实质上等于π-in电子的贡献）。由数据可见，B6C6N6的芳香性显著强于B9N9，并且π-out的芳香性比π-in还略强一点，这可能在于体系的π-in电子的共轭程度相对略低，源自于环的曲率使得平行于环的p原子轨道彼此交叠程度低于垂直于环平面的p原子轨道。单从AV1245*1000的数值来看，B6C6N6的芳香性甚至比18碳环还要更强一点。 NICS是非常流行和稳健的基于磁属性衡量芳香性的指标，它的不同形式在//www.umsyar.com/176里有专门介绍。其中NICS(0)zz对本文研究的这些环状共轭体系最为适合，结果如下所示，数值越负说明芳香性越强。为了更充分了解芳香性的来源，文中还按照《基于Gaussian的NMR=CSGT任务得到各个轨道对NICS贡献的方法》（//www.umsyar.com/670）介绍的方法分别计算了π-in、π-out和σ+core电子的贡献。NICS(0)zz体现出B6C6N6的芳香性接近18碳环，并且π-in和π-out的贡献相仿佛。而B9N9虽然形式上也满足4n+2规则，但由于其NICS(0)zz基本为0，因此可断定是非芳香性的，这必然是由于其π电子缺乏整体离域能力所致。虽然B9N9的各个B-N键如前述Mayer键级和IRI-π所示有明显π作用，但显然不意味着其π电子能够容易地离域。对这些体系，σ和core电子都由于其高度定域性而对芳香性基本没有影响。 《使用AICD程序研究电子离域性和磁感应电流密度》（//www.umsyar.com/147）和《使用AICD 2.0绘制磁感应电流图》（//www.umsyar.com/294）介绍的AICD图是使用得较普遍的图形化直观展现化学体系磁感生电流的方法。为了进一步从磁感生电流的角度展现B6C6N6的芳香性以及和B9N9的差异，文中对这些体系的π-in、π-out以及全部π电子分别绘制了AICD图，如下所示。外磁场方向垂直于环由下朝上施加。由于等值面上的箭头太小，为了看得清楚，图中还用大的橙色箭头做了清晰的标注。此图直观体现出，B6C6N6确实存在显著的π芳香性，因为感生电流方向符合左手定则并且连贯地环绕整个体系，这正是典型的π芳香性体系具备的特征。而缺乏芳香性的B9N9则明显不具备这样的特征，感生电流几乎只是绕着原子核或原子间区域转。 此文还按照《考察分子磁感生电流的程序GIMIC 2.0的使用（含24分钟演示视频）》（//www.umsyar.com/491）介绍的方法绘制了GIMIC磁感生电流的动画，如下所示，更进一步地体现出芳香性明显的B6C6N6能够形成整体感生环电流而非芳香性的B9N9不具备这样的特征。 B6C6N6：//www.umsyar.com/attach/696/B6C6N6_GIMIC.mp4 B9N9：//www.umsyar.com/attach/696/B9N9_GIMIC.mp4 文中还对穿越B6C6N6、B9N9和18碳环的化学键的感生电流进行了积分，数值分别为21.64、0.70、22.50 nA/T，这定量体现出B6C6N6和18碳环的芳香性相仿佛，而B9N9完全不具备芳香性。 在《通过Multiwfn绘制等化学屏蔽表面(ICSS)研究芳香性》（//www.umsyar.com/216）和《使用Multiwfn巨方便地绘制二维NICS平面图考察芳香性》（//www.umsyar.com/682）中介绍了ICSS函数，它直观地体现了化学体系的磁屏蔽和去屏蔽区域的位置和强度，是NICS的重要扩展，已被非常广泛应用于讨论芳香性问题。此研究使用Multiwfn对B6C6N6和B9N9都绘制了ICSS_ZZ = 10 ppm的等值面图以及垂直于环的截面的填色等值线图，如下所示。等值面图中橙色和蓝色分别对应屏蔽和去屏蔽区域。由图可见，B6C6N6的环内部完全是屏蔽区域，环外面完全是去屏蔽区域，这种特征和苯这样的典型芳香性体系完全相同，更进一步严格地证明了B6C6N6有显著的芳香性。而且B6C6N6的ICSS的这种特征和Carbon, 165, 468-475 (2020)报道的18碳环的高度一致，体现出二者芳香性的高度共性。而B9N9的ICSS等值面图的屏蔽和去屏蔽区域彼此交错，没有形成遍及整体的连贯的等值面，体现出此体系既没有芳香性也没有反芳香性，而是非芳香性。 5 B6N6C6的电子离域性 芳香性的本质来自于电子的离域，因此文中还特意对B6N6C6的电子离域性从多方面进行了考察。在《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）中介绍了Multiwfn可以计算和绘制的LOL-π函数，这对于考察π电子的离域情况极其有用，如今已经非常流行。下图对比了B6N6C6和B9N9的LOL-π=0.5等值面，明显可见B6N6C6的等值面比B9N9显得连贯得多，明显体现出B6N6C6的π电子有强得多的在18个原子上全局离域的能力，也一定程度解释了前者的芳香性比后者强得多的原因。 LOL-π函数比较适合直观定性考察，而《在Multiwfn中单独考察pi电子结构特征》里介绍的ELF-π函数的二分点数值则适合定量分析对比。下图是B6N6C6和B9N9的ELF-π函数等值面图，等值数值设为了令等值面刚好分裂时的数值，这被称为二分值，数值越大说明电子越容易跨越二分点区域发生离域。在图上也标注了二分点位置和二分值。由图可见，B6N6C6的π-in和π-out的二分值都明显大于B9N9的，这给B6N6C6的整体电子离域性更强又进一步提供了定量证据。 电子的费米穴是考察电子离域特征的一个重要的函数，它包含两个三维位置坐标，作图考察时通常将一个位置作为参考点，然后对另一个坐标绘图。Atoms-in-molecules理论的提出者Bader指出an electron can go where its hole goes and, if the Fermi hole is localized, then so is the electron，也就是假设一个电子出现在某个参考点位置，那么其费米穴函数明显分布在哪，这个电子就容易离域到哪去。因此，通过绘制费米穴图，可以更进一步了解体系处于不同位置的电子的离域能力。文中通过Multiwfn对B6N6C6、B9N9和18碳环分别绘制了费米穴的二维图，如下所示。参考点位置以绿色箭头标注，设在了不同原子的价层区域距离原子核特定距离的位置，此时图中的费米穴的分布就体现了这个位置的价电子容易离域到哪去。用Multiwfn绘制这种图的操作在《制作动画分析电子结构特征》（//www.umsyar.com/86）里有相关说明。 由上图等值线分布可见，B6N6C6的硼和碳的价电子的费米穴分布范围较广，体现出电子有往较远区域离域的能力，特别是其碳的费米穴分布特征和已知具有较显著全局离域性的18碳环的很接近。相比之下，B9N9的硼和氮的费米穴分布广度明显不及B6N6C6的相应元素的，等值线只分布在距离参考点较近的原子上，体现出B9N9的价电子的离域能力相对较差。 上面的等值线图适合定性讨论，为了能把以上信息转化成定量形式来更好地对比分析，文中定义了一个新的量叫做原子远程离域指数（atomic remote delocalization index, ARDI）。离域性指数（DI）在《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）里专门介绍过，是衡量特定两个原子间等效共享电子对数用的，而ARDI则是对单个原子定义的，定义为距离它两个化学键以上的原子与当前原子的所有DI之和，显然基于Multiwfn计算的所有原子间的DI可以简单地手算出来各个原子的ARDI。ARDI体现了特定原子上的电子的远程离域能力，数值越大说明远程离域能力越强。环上各个原子上的电子有显著的远程离域能力显然是这个环具备明显芳香性的前提。各原子ARDI如下所示，可见B6N6C6上的原子的ARDI都较大，平均值与18碳环的碳的相仿佛，而B9N9的硼和氮的ARDI都远小于B6N6C6相应元素的原子，这更充分体现了B9N9的电子的全局电子离域能力远不如B6N6C6，进一步展现了B9N9芳香性远低于B6N6C6的原因。 6 碳掺杂进B9N9使得芳香性巨幅提升的本质 根据前面的讨论，文章已经充分证明了芳香性的关系是B6C6N6 ≈ C18 ≫ B9N9。为什么将碳元素掺杂进入无芳香性的B9N9，就能使其芳香性巨幅提升、电子全局离域能力大幅增加？实际上原因并不难理解。尽管前面的Mayer键级体现出B9N9中的硼-氮之间有一定π作用，但由于硼和氮的电负性差异非常大，这无疑会导致电子倾向于在氮上定域而难以往相邻的硼上离域，至于往更远的原子上离域就更是难上加难了。而碳的电负性则正好介于硼和氮之间，显然硼-碳和氮-碳键的共价性显著强于硼-氮键，因此电子更容易离域过去，相当于碳给硼和氮之间架设了能够令电子容易离域过去的桥梁，充分调和了硼和氮之间电负性差异过大导致的电子跨键离域太难的矛盾。 为了更好地展现碳原子起到的价值，文中将垂直于环的各原子的2p原子轨道绘制了出来，π-out电子正是分布在这些p轨道上，π-out分子轨道也正是它们线性组合产生的。具体来说，由于NBO程序产生的预自然原子轨道（PNAO）可以较合理地描述分子环境中的原子轨道，因此文中按照《使用Multiwfn绘制NBO及相关轨道》（//www.umsyar.com/134）中的做法，用Multiwfn结合VMD将相应的PNAO轨道用等值面图形式做了展示，如下所示。等值面数值选为了适合对比的0.1。图中也把PNAO的能量（eV）标注在了上面。由图可清楚地看出，B9N9中的氮和硼的p轨道的空间延展程度差异很大，而且轨道能量差异也很大，势必电子难以在它们之间离域。而从B6C6N6的图中可见，碳的p轨道延展范围和轨道能量都恰好介于硼和氮之间，它的引入无疑能够显著帮助电子离域过去。 7 B6N6C6异构体的芳香性 以上研究的那种以B-C-N作为重复单元排列的B6N6C6结构是否真的是B6N6C6的能量最低结构？是否B6N6C6化学组成的体系还有芳香性更强的异构体？这是个值得考察的有趣的问题。为此，文中还另外考虑了两种B6N6C6构型，优化后得到的无虚频结构如下所示，可见B6N6C6'构型让所有碳尽可能连在一起，其它部分部分保持B9N9那样硼-氮交错的结构，而B6N6C6''构型则是由三个C-C-B-N-B-N重复单元构成的结构。这两个结构的基态是闭壳层的，能量分别比前面研究的B-C-N为重复单元构成的B6N6C6结构能量低280.0和250.2 kcal/mol。这体现出B6N6C6中B-C和C-N的平均键能比B-N和C-C是要更小的。虽然前面研究的B6N6C6结构能量不是这种化学组成的能量最低结构，但并不意味着它是不稳定的。根据《使用ORCA做从头算动力学(AIMD)的简单例子》（//www.umsyar.com/576）中介绍的方法，文中对此体系做了较高温度（500 K）下20 ps长度的基于wB97X-D3/6-311G*级别的从头算分子动力学模拟，用VMD每1 ps绘制一次的结构叠加图如下所示（颜色按照蓝-白-红变化），可见环状结构始终很好维持着，并没有发生异构化、解离等现象，体现出B6N6C6具有一定热力学稳定性。 文中还对B6N6C6'和B6N6C6''构型计算了NICS(0)zz，发现分别仅有1.8和1.4 ppm，可见这俩构型没有芳香性。这说明只有让碳插入到硼和氮中间，从而以B-C-N作为重复单元时，碳才能真正起到显著提升纯硼-氮环体系的电子离域性的作用。 8 B6N6C6最低三重态激发态的芳香性 前面考察的都是B6N6C6的单重态基态（S0）。最后，文章还研究了B6N6C6最低三重态激发态（T1态）的芳香性。计算发现S0-T1垂直和绝热激发能分别为1.62和1.17 eV。下图(a)是优化后的T1态B6N6C6的结构，可见对称性相对于S0态发生了下降。下图(b)是按照《使用Multiwfn作电子密度差图》（//www.umsyar.com/113）绘制的在B6N6C6的S0结构下T1态与S0态的密度差，蓝色和绿色分别为电子密度降低和增加部分。由图可见S0到T1态的电子激发是π-in → π-out的激发。由于垂直激发导致的电子密度变化破坏了B6N6C6基态的C6h对称性，这也是为什么在T1态势能面上从Franck-Condon点往极小点弛豫过程中几何结构对称性发生了下降。 文中对B6N6C6的T1态计算了NICS(0)zz，并分解为了不同类型电子的贡献。而且为了考察几何结构弛豫的影响，在S0极小点和T1极小点分别做了这种计算，结果如下所示。可见无论在哪种结构下，T1态由于NICS(0)zz远大于0，因此是显著反芳香性的，这和Baird规则所述情况一致。Baird规则指出单重态基态满足Hückel 4n+2规则的体系的最低三重态是反芳香性的。从NICS(0)zz还看到，T1态的反芳香性特征在T1极小点结构下比S0极小点结构下更弱，这体现出B6N6C6自发减小T1态的反芳香性可能是其S0→T1垂直激发后结构弛豫的重要驱动力。 9 总结 此文基于 计算和Multiwfn等程序做波函数分析，非常全面系统研究了18碳环的等电子体B6N6C6的几何结构和电子结构，从不同角度全面考察了其芳香性和电子离域性特征，并且和18碳环及其另外一个重要的等电子体B9N9进行了细致的对比。本文体现出以碳原子桥联硼和氮原子，可以显著提升纯硼-氮体系的电子离域性。这不仅对电子结构产生重要影响，必然也同时会影响其它性质，诸如光学性质、电子输运性质等。本文虽然研究的只是一维体系，可以预想电负性介于硼和氮之间的碳元素的掺杂必定也会对二维、三维纯硼-氮体系产生显著影响，而且掺杂方式的选取可能会对这些体系的性质起到一定调控作用。 如果读者对于18碳环及其相关体系感兴趣，非常推荐查看//www.umsyar.com/carbon_ring.html里汇总的本文作者的巨量其它研究工作，包含大量深入浅出的论文内容和研究思想介绍文章，以及对计算技术细节的附加说明。 //www.umsyar.com/693 2023-12-30T12:35:00+08:00 使用Dalton通过CC3方法极高精度计算激发态 Extremely high-precision calculation of excited states using the CC3 method in Dalton program 文/Sobereva@北京科音   2023-Dec-30 0 前言 耦合簇方法可以通过线性响应（linear response, LR）计算激发态。在常见范畴内，精度和耗时关系是LR-CC3>LR-CCSDR(3)≈EOM-CCSD(T)>LR-CCSD=EOM-CCSD>=CC2。LR-CC3如今被普遍视为是激发态计算的金标准（多参考特征很强以及双电子激发特征占主导的情况除外），激发能能精确到零点零几eV。但LR-CC3结合像样基组只能算得动个位数原子，对更大体系则可以用明显更便宜的LR-CCSDR(3)，精度也已经非常好了。 Dalton程序在耦合簇方面功能丰富，可以做LR-CC3、LR-CCSDR(3)、LR-CCSD、LR-CC2激发能计算，对于LR-CC2/CCSD/CC3还可以得到包括振子强度在内的诸多跃迁属性。Dalton的输入文件对于初学者来说比较复杂抽象，本文用尽可能简单易懂的语言完整介绍一下用免费的Dalton做这些激发态耦合簇计算的方法，使用有C2v对称性的非常典型的分子甲醛为例。读者如果对Dalton不熟悉，务必先阅读我之前写的《 程序Dalton的编译方法和运行方式简介》（//www.umsyar.com/463），我假定读者已经认真看过此文。 值得一提的是，虽然Dalton做LR-CC3计算如今在文献中很常见，但Dalton在这方面不是最快的，而且代码还没有并行化。如果追求更好的效率，可以用专门做耦合簇计算的e^T程序（https://www.etprogram.org），它的LR-CC3是所有程序中明显最快的。另外，LR-CCSD激发能和振子强度没任何必要非得用Dalton计算，常用的Gaussian和ORCA的EOM-CCSD也都做得很好，结果和Dalton的LR-CCSD是等同的。Dalton做耦合簇的激发能计算不支持考虑溶剂模型，而G16的EOM-CCSD则可以。还值得一提的是，Dalton的衍生程序LSDalton还额外支持RI-CC2，可以比Dalton的LR-CC2更高效率地计算较大体系，这不属于本文的范畴。 本文涉及的输入输出文件都可以在//www.umsyar.com/attach/693/file.zip里获得。计算用的是2023-Dec-22通过git下载的Dalton最新的开发版。 本文下面的例子首先要重复知名的TBE1激发能测试集（J. Chem. Phys., 128, 134110 (2008)）里的LR-CC3/def-TZVP算的甲醛的1-1A2、1-1B1、2-1A1三个单重态的垂直激发能。如文章表1所示，结果分别为3.95、9.18、10.45 eV。此文用的几何结构是MP2/6-31G*级别优化的，对应的结构文件是file文件包里的H2CO.xyz，本文的计算也将基于这个结构。 1 确定不可约表示顺序 为了指定各个不可约表示的激发态各算几个态，需要先知道C2v点群的各个不可约表示在Dalton程序中的顺序。最简单的方法是做个DFT单点计算，看一开始输出的不可约表示的顺序。为此，我们用Multiwfn（//www.umsyar.com/multiwfn）创建这个任务的输入文件。启动Multiwfn，然后依次输入 H2CO.xyz 100   //其它功能（Part 1） 2   //导出各种文件 19   //Dalton DFT.dal  //在当前目录下产生DFT.dal，对应B3LYP/6-31G*计算设置 H2CO.mol   //在当前目录下产生Dalton格式的体系定义文件H2CO.mol 由于当前计算要考虑对称性，因此把H2CO.mol里的Nosymmetry删掉。然后基于DFT.dal和H2CO.mol做计算，从输出文件DFT_H2CO.out里可以看到在SCF开始之前就输出了当前自动判断的点群，以及相应的各个不可约表示的顺序 @    Full point group is: C(2v)          @    Represented as:      C2v @  * The irrep name for each symmetry:    1: A1     2: B1     3: B2     4: A2 2 做LR-CC3激发能计算 写一个文本文件LRCC3.dal，内容如下（这是最简单写法，默认设置就已经适合的选项就没再做额外设置） **DALTON .RUN WAVE FUNCTIONS **WAVE FUNCTION .CC *CC INP .CC3 *CCEXCI .NCCEXCI 2 1 0 1 **END OF DALTON 对当前体系默认是做闭壳层HF计算得到参考态波函数。**DALTON控制任务类型，.RUN WAVE FUNCTIONS要求算单点并得到波函数。**WAVE FUNCTION设置波函数计算类型，.CC要求做耦合簇计算。*CC INP设置做什么形式耦合簇计算，.CC3要求做CC3计算。*CCEXCI要求做耦合簇的激发能计算，.NCCEXCI下面按照不可约表示的顺序指定各个不可约表示的能量最低激发态各计算几个。当前为了重复TBE1测试集的数据，要求A1算两个（由此得到1-1A1和2-1A1），B1算1个（得到1-1B1），B2不算，A2算1个（得到1-1A2）。 然后把前述的H2CO.mol复制为H2CO_TZVP.mol，手动把里面的6-31G*替换为Turbomole-TZVP，这是Dalton内置的def-TZVP基组的写法，对应Dalton目录下的basis子目录中的Turbomol-TZVP这个文件的名字。 然后基于LRCC3.dal和H2CO_TZVP.mol做计算，得到LRCC3_H2CO_TZVP.out。本节的文件都已经提供在了前述的file文件包里。打开输出文件，从里面可以看到如下激发能信息，2-1A1为10.44842 eV，1-1B1为9.18343 eV，1-1A2为3.94711 eV，可见和TBE1原文里的10.45、9.18、3.95 eV完全一致。表中||T1||是单激发贡献，和TBE1原文里给出的也是一致的。  +=============================================================================+  |  sym. | Exci.  |        CC3        Excitation energies            | ||T1||  |  |(spin, |        +------------------------------------------------------------+  | spat) |        |     Hartree    |       eV.      |     cm-1       |    %    |  +=============================================================================+  | ^1A1  |    1   |     0.3502215  |       9.53001  |     76864.734  |  91.00  |  | ^1A1  |    2   |     0.3839723  |      10.44842  |     84272.170  |  91.32  |  +-----------------------------------------------------------------------------+  | ^1B1  |    1   |     0.3374848  |       9.18343  |     74069.360  |  90.93  |  +-----------------------------------------------------------------------------+  | ^1A2  |    1   |     0.1450537  |       3.94711  |     31835.610  |  91.16  |  +=============================================================================+ 如果想要计算三重态激发态，就在.NCCEXCI下面的第二行定义。比如下面的写法，代表除了计算如上那些单重态激发态以外，还计算1-3A2和1-3A1三重态激发态 .NCCEXCI 2 1 0 1 1 0 0 1 如果不需要计算单重态激发态，就写 .NCCEXCI 0 0 0 0 1 0 0 1 笔者试了几个Dalton版本，包括2016、2018、2022开发版，算出来的三重态激发能都是错的，和TBE1表2里的1-3A2和1-3A1三重态激发能明显不符，而且||T1||严重偏低，我认为是bug。如果不利用对称性，即.mol文件里带着Nosymmetry，并且.NCCEXCI下面直接指定计算的激发态的总态数，则问题轻得多，可以看到||T1||都显著高于90%。 值得一提的是LR-CC3是对激发态一个一个独立计算的，故算的态数越多耗时明显越高。 如果只需要计算CC3基态能量，去掉*CCEXCI及下面的部分即可，相对于LR-CC3算激发态部分的耗时来说可以忽略不计。 3 做LR-CC3激发能+振子强度的计算 如果对上面算的激发态还要计算振子强度，就创建LRCC3mom.dal文件，写入以下内容，然后进行计算。其中**INTEGRAL指定要算基函数之间哪些类型的积分，.DIPLEN要求算长度形式的偶极积分，这是因为算振子强度要算跃迁偶极矩，会用到这类积分。*CCOPA模块用来计算耦合簇的基态到激发态的跃迁属性，对当前用的Dalton版本支持CCS、CC2、CCSD、CC3，其中CC3仅限单重态激发态。.DIPOLE要求计算长度形式的跃迁偶极矩及振子强度。 **DALTON .RUN WAVE FUNCTIONS **INTEGRAL .DIPLEN **WAVE FUNCTION .CC *CC INP .CC3 *CCEXCI .NCCEXCI 2 1 0 1 *CCOPA .DIPOLE **END OF DALTON 从输出文件（file文件包里的LRCC3mom_H2CO_TZVP.out）可以看到各个激发态的跃迁偶极矩和振子强度，比如2-1A1的如下所示，振子强度为0.34867845。这个值很接近TBE1原文表VII里这个态的LR-CC2和LR-CCSD振子强度（这篇文章里没给CC3的振子强度），分别为0.368和0.374。      Transition from ground state to:      number, multiplicity, symmetry :    2  ^1A1      frequency :   0.3839722612 a.u.    10.44842 e.V.     84272.2 cm^-1      +-----------+-----------------+-----------------+---------------------+      | operator  |   left moment   |  right moment   | transition strength |      +-----------+-----------------+-----------------+---------------------+      | XDIPLEN   |      0.00000000 |      0.00000000 |       0.00000000    |      | YDIPLEN   |      0.00000000 |      0.00000000 |       0.00000000    |      | ZDIPLEN   |     -0.83921304 |     -1.62309631 |       1.36212358    |      +-----------+-----------------+-----------------+---------------------+        oscillator strength (length gauge)   :      0.34867845 上面表格里的ZDIPLEN（偶极矩的Z分量算符）的transition strength值1.36212358对应的是跃迁偶极矩Z分量的平方。 如果你只需要跃迁偶极矩的特定分量，比如Z分量，可以把.DIPOLE改为 .OPERATOR ZDIPLEN 这样耗时会比计算所有分量时稍微低一点。 4 LR-CC2/CCSD/CCSDR(3)激发态计算 LR-CC2、LR-CCSD、LR-CCSDR(3)计算只需要把前面例子里的.CC3分别改成.CC2、.CCSD、.CCR(3)即可。注意LR-CCSDR(3)只能算激发能而无法算包括振子强度在内的跃迁属性。 这些级别的激发能计算的输入输出文件在file文件包里的other目录下都给了，耗时跟LR-CC3比都不值得一提。整体来说按LR-CC2、LR-CCSD、LR-CCSDR(3)的顺序，结果和金标准LR-CC3的越来越接近。比如它们的1-1B1激发能分别为9.348、9.257、9.186 eV，LR-CC3的为9.183 eV。所以LR-CC3激发能很难算得动的时候用LR-CCSDR(3)是很好的平替。 5 冻核设置 Dalton默认是不冻核的。显然恰当冻核可以节约耗时而几乎不影响精度。这里演示对H2CO做LR-CC3计算时冻结能量最低两个分子轨道（对应C和O的1s轨道）的方式。首先得得知它们的不可约表示，从之前B3LYP/6-31G*的输出文件中可以看到各个占据轨道能量的负值（Koopmans定理对应的电离能，相关知识参考《正确地认识分子的能隙(gap)、HOMO和LUMO》//www.umsyar.com/543）。 @ Summary of DFT KT ionization potentials [eV]: @ Symmetry 1 (A1 ) :   521.511     279.798      28.621      17.429      12.139 @ Symmetry 2 (B1 ) :    10.711 @ Symmetry 3 (B2 ) :    13.454       7.303 @ Symmetry 4 (A2 ) :  No IPs in this symmetry 显然能量最低两个轨道是A1。 写一个LR-CC3激发能计算输入文件LRCC3_FC.dal，内容如下。其中.FROIMP下面第一行指定对各个不可约表示冻结多少个能量最低占据轨道，下面第二行指定对各个不可约表示冻结多少个能量最高空轨道。当前冻结的是能量最低的两个A1轨道。 **DALTON .RUN WAVE FUNCTIONS **WAVE FUNCTION .CC *CC INP .CC3 .FROIMP 2 0 0 0 0 0 0 0 *CCEXCI .NCCEXCI 2 1 0 1 **END OF DALTON 从file文件包里的此任务的输出文件LRCC3_FC_H2CO_TZVP.out可看到，现在冻核后耗时是1分49秒，1-1B1激发能是9.18872 eV，之前没冻核时是3分32秒，1-1B1激发能是9.18343 eV。明显冻核对激发能的影响可完全忽略不计，而耗时大为降低，因而有重要实际意义。