//www.umsyar.com/category/QC //www.umsyar.com/734 2024-12-30T02:00:00+08:00 使用Multiwfn计算轨道的体积 Using Multiwfn to calculate volume of orbitals 文/Sobereva@北京科音  2025-Jan-1 1 前言 有人在思想家公社QQ群（//www.umsyar.com/QQrule.html）里问怎么计算分子轨道的体积，并且贴了一张轨道等值面图。只要计算出这个等值面里包围的体积就可以达到这个目的（虽然也可以有其它方式的定义，但没这么简单直观）。这里我介绍一下怎么利用Multiwfn程序的定量分子表面分析功能来实现这个目的。定量分子表面分析功能之前在《使用Multiwfn的定量分子表面分析功能预测反应位点、分析分子间相互作用》（//www.umsyar.com/159）有简要介绍，结合Multiwfn的极度灵活的设计，这个功能还能做很多其它的分析，正如本文所展示的。 如果读者不熟悉Multiwfn，应当阅读《Multiwfn FAQ》（//www.umsyar.com/452）、《Multiwfn入门tips》（//www.umsyar.com/167）、《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）。Multiwfn可以在//www.umsyar.com/multiwfn免费下载，不要用上古版本。 2 实例 这里以Multiwfn计算下面这个D-pi-A体系的HOMO轨道的体积为例，图中的轨道等值面数值为0.05，此体系的波函数文件是Multiwfn目录下的examples\excit\D-pi-A.fchk。注意等值面包围的体积直接依赖于等值面数值的选取，横向对比时必须统一。 启动Multiwfn，然后输入 examples\excit\D-pi-A.fchk 5  //计算格点数据 4  //轨道波函数 h  //HOMO 4  //设置格点间距 0.25  //0.25 Bohr是精度和耗时的较好权衡。体系大的话也可以用大一些的格点间距来节约耗时但会损失一些精度 0  //返回主菜单 13  //处理格点数据 11  //格点数据运算 13  //取绝对值。轨道波函数有正有负，这一步使之都变成正值 -1  //返回主菜单 12  //定量分子表面分析 1  //选择定义表面的方式 11  //用内存中的格点数据定义等值面 0.05  //等值面数值 6  //开始表面分析且不考虑被映射的函数 从屏幕上可以看到以下结果，即体积为14.2 Angstrom^3，顺带着轨道等值面的面积69.0 Angstrom^2也显示出来了。  Volume:    95.86111 Bohr^3  (  14.20515 Angstrom^3)  Estimated density according to mass and volume (M/V):   25.0417 g/cm^3  Overall surface area:         246.50998 Bohr^2  (  69.02983 Angstrom^2) 感兴趣的话，还可以选择选项-3 Visualize the surface看一下当前等值面对应的图像，体积就是这个等值面里面包围的部分 对比一下，看看下面这个第18号分子轨道的体积如何 退回到主菜单，把下面这一串命令复制到Multiwfn窗口里即可，结果是7.3 Angstrom^3，可见只有HOMO的一半，和肉眼看上去的轨道体积差异相一致。 5 4 18 4 0.25 0 13 11 13 -1 12 1 11 0.05 6 3 批量计算 在《详谈Multiwfn的命令行方式运行和批量运行的方法》（//www.umsyar.com/612）里专门讲了怎么用Multiwfn批处理计算。这里提供一个脚本，可以便捷地用前文的方法批量计算特定序号范围的轨道的体积。//www.umsyar.com/attach/734/orbvol.sh是个Bash shell脚本，用来计算examples\excit\D-pi-A.fchk这个体系所有占据轨道（1到56号）的体积，内容如下所示。其中istart和iend是起始和终止的轨道序号，这俩值和被处理的波函数文件路径都可以根据实际需要修改。运行之前，Multiwfn可执行文件所在目录需要添加到PATH环境变量下使得能够通过Multiwfn命令启动之。 #!/bin/bash istart=1 iend=56 for ((i=$istart;i<=$iend;i=i+1)) do cat << EOF >> orbvol.txt 5 4 $i 4 0.25 0 13 11 13 -1 12 1 11 0.05 6 -1 -1 EOF done echo "q" >> orbvol.txt echo "Running Multiwfn..." Multiwfn examples/excit/D-pi-A.fchk < orbvol.txt >> result.txt grep "Volume:" result.txt | nl -v$istart rm ./orbvol.txt result.txt 运行后，屏幕上看到各个轨道的体积：      1   Volume:     2.88332 Bohr^3  (   0.42726 Angstrom^3)      2   Volume:     2.89404 Bohr^3  (   0.42885 Angstrom^3)      3   Volume:     2.42980 Bohr^3  (   0.36006 Angstrom^3)      4   Volume:     2.40622 Bohr^3  (   0.35656 Angstrom^3) ...略     54   Volume:    95.10824 Bohr^3  (  14.09359 Angstrom^3)     55   Volume:    95.86240 Bohr^3  (  14.20534 Angstrom^3)     56   Volume:    95.86111 Bohr^3  (  14.20515 Angstrom^3) 将轨道体积相对于轨道序号作柱形图，如下所示。可见前16号轨道的体积非常小，这是因为它们都是内核轨道。 笔者之前写过《通过轨道离域指数(ODI)衡量轨道的空间离域程度》（//www.umsyar.com/525），里面用我提出的ODI指数考察了与本文相同的D-pi-A体系的各个轨道的离域程度，并且绘制了ODI与轨道序号之间的柱形图。将上图与ODI的柱形图相对比，可以看到两个图有明显的互补性，即ODI越低（体现轨道越离域、越能同时覆盖很多原子），轨道体积倾向于越大。如果你的目的就是想衡量轨道的离域程度，ODI相对更严格，毕竟不依赖于轨道等值面数值的选取（轨道等值面数值设得很大的话，轨道等值面就完全看不到了，体积为0）。而用轨道体积来说事的好处是可以直接与等值面图相对应，结合图像讨论非常清楚直观。   4 其它 本文的做法不限于计算考察分子轨道的体积，对任意轨道都可以考察，如定域化轨道、NTO轨道、双正交化轨道、NAdO轨道、AdNDP轨道等等，载入含有相应轨道的波函数文件即可，产生方法在下文说了。 Multiwfn的轨道定域化功能的使用以及与NBO、AdNDP分析的对比 //www.umsyar.com/380 用于非限制性开壳层波函数的双正交化方法的原理与应用 //www.umsyar.com/448 使用键级密度(BOD)和自然适应性轨道(NAdO)图形化研究化学键 //www.umsyar.com/535 使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键 //www.umsyar.com/138   效仿本文的做法还可以计算其它任意实空间函数的等值面内包围的体积。比如计算自旋密度的等值面内的体积，只需要在主功能5里选择被计算的函数的那一步选择自旋密度即可。如果是像ELF那样处处为正的函数，计算其等值面内的体积前无需像前文那样需要先进入主功能13将其取绝对值；反之如果函数处处为负，则也需要先取绝对值。 //www.umsyar.com/733 2024-12-30T01:58:00+08:00 谈谈分子间结合能的构成以及分解分析思想 On the components of intermolecular binding energy and idea of decomposition analysis 文/Sobereva@北京科音  2024-Dec-31 在我网上回答大量 的问题的过程中，经常看到有些人搞不清楚对分子间相互作用有贡献的量之间的关系、不知道怎么将结合能分解成不同部分考察影响结果的内在因素。这使我决定作一张图，把分子间结合自由能的构成一次性梳理清楚，同时展现我的对分子间结合能的分解分析的思想，如下所示： 下面我会把图中各项进行简要的解释，便于读者理解每项是怎么回事、知道怎么计算，务必结合上图看。上图的源头是“溶剂中的结合自由能”，因为这个是包含的成份最多，内涵最复杂的量。本文仅考虑两个单体分子A、B间结合成为复合物AB的情况。 • 溶剂中的结合自由能ΔG_bind(sol)：这是指在溶剂（sol）环境下，两个分子从各自孤立状态相互接近并形成二聚体过程中的自由能的总变化。即ΔG_bind(sol) = G_AB(sol) - G_A(sol) - G_B(sol)。如《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》（//www.umsyar.com/327）所述，溶剂下的体系的自由能G(sol)可以写为气相自由能G(gas)和溶解自由能ΔG_solv的加和，因此可以分解成两个部分： • 溶解自由能的贡献Δ(ΔG_solv)：这体现溶剂效应对分子间结合产生的影响，即复合物的溶解自由能减去每个分子的溶解自由能。这部分可正可负，为正代表溶剂效应不利于结合，通常出现在两个单体间以静电作用为主的情况，即溶剂效应对单体的总稳定化程度比对复合物的稳定化程度更显著；为负代表溶剂效应有利于分子间结合，例如《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》（//www.umsyar.com/727）所讲的情况，溶剂的疏水效应能显著促进富勒烯和碳环之间的结合。 溶解自由能可分为极性和非极性两部分，因此讨论溶剂效应对分子间结合产生的影响时可以具体讨论这两部分各自的影响是多大。隐式溶剂模型下溶解自由能的标准计算方式在前述的《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》里说了，是溶剂模型下的电子能量减去气相下的电子能量。Gaussian做考虑SMD溶剂模型的单点计算时可以直接给出非极性部分对溶解自由能的贡献，把溶解自由能减去这部分就是极性部分对溶解自由能的贡献。复合物的溶解自由能的极性（非极性）贡献减去每个单体的这部分，就是溶剂效应的极性（非极性）部分对结合自由能的贡献。 • 气相结合自由能ΔG_bind(gas)：通过G_AB(gas) - G_A(gas) - G_B(gas)得到。气相下G的计算方式见《使用Shermo结合 程序方便地计算分子的各种热力学数据》（//www.umsyar.com/552）。注意，如果你本身就是研究溶剂环境下的结合，而且溶剂效应对几何结构的影响不可忽视，那么复合物和单体的opt和freq任务都应当在溶剂模型下做，即计算气相自由能用的几何结构是在溶剂模型下优化的，自由能热校正量也是基于溶剂模型下振动分析信息得到的。 由于G=H-TS，故ΔG_bind = ΔH_bind - TΔS_bind，遂有了下面两部分。 • ΔH_bind：ΔG_bind(gas)中的焓效应部分，即复合物的焓减去各个单体的焓 • -TΔS_bind：ΔG_bind(gas)中的熵效应部分，即复合物的熵减去各个单体的熵然后乘以-T。由于分子间结合会令熵大幅降低，ΔS_bind为负，因此这部分总是不利于分子间结合，称为熵罚，而且温度（T）越大明显越不利结合。如果分子间无法形成足够强烈的相互作用令ΔH_bind足够负，那么它很可能被-TΔS_bind抵消掉导致分子间无法结合。我在北京科音初级 培训班（http://www.keinsci.com/KEQC）里给了相关计算例子，体现出在气相标况下，水二聚体在热力学上是无法自发形成的，但是能在分子间形成两个显著的氢键的甲酰胺二聚体是能形成的。 温度对ΔH_bind的影响较小，而对-TΔS_bind影响巨大，因此分子间能否结合极大程度取决于温度。在《8字形双环分子对18碳环的独特吸附行为的 、波函数分析与分子动力学研究》（//www.umsyar.com/674）介绍的笔者的文章Phys. Chem. Chem. Phys., 25, 16707 (2023)里通过Shermo程序做了碳环、8字形分子OPP以及二者的复合物的热力学量相对于温度的扫描，并进而给出了ΔG_bind、ΔH_bind、-TΔS_bind随温度的变化，准确指出了分子间能结合的临界温度。前述的《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》介绍的笔者的Chem. Eur. J., 30, e202402227 (2024)也做了这样的分析。 如Shermo程序手册附录部分所讲的热力学数据计算常识所示，熵是由体系的平动、转动、振动以及电子（跃迁或简并）共同贡献的，Shermo可以分别输出，而且Shermo还能给出各个振动模式对熵的贡献。通过这些信息可以深入了解熵的构成，并对各项在复合物和单体之间求差以进一步了解ΔS_bind的内在本质。 由于焓等于电子能量与焓的热校正量之和，即H = E + H_corr，故ΔH_bind又可以分为以下两部分： • 气相结合能ΔE_bind：即E(AB)-E(A)-E(B)，其中E(AB)、E(A)、E(B)分别是AB、A、B的气相的（即不带溶剂模型的）电子能量，并且计算用的几何结构是对它们各自分别优化后的。ΔE_bind体现出在气相下，A、B从孤立状态下结合为AB过程中的电子能量的变化。如果你不清楚什么叫电子能量，看《谈谈该从Gaussian输出文件中的什么地方读电子能量》（//www.umsyar.com/488）。注意，平时说结合能的时候可以指ΔG_bind也可以指ΔE_bind，只有前者才能体现热力学上结合的可能性，后者只能用来讨论分子间内在的作用强度、结合的驱动力。没有前提的情况下说结合能的时候习俗上指的是更容易算的ΔE_bind，本文也是这个习俗。 • ΔH_corr：即复合物的H_corr减去各个单体的。这体现出分子间结合过程中热力学效应对焓变的贡献。理想气体近似下H=U+RT，故双分子结合的情况有ΔH_corr = ΔU_corr - RT，其中ΔU_corr是结合过程中内能的热校正量的变化。类似于熵，U_corr也是能够通过Shermo输出平动、转动、振动（包括ZPE在内）、电子跃迁的贡献从而了解ΔU_corr的内在本质的。 ΔE_bind = ΔE_int + ΔE_deform，即气相结合能可以分解为气相相互作用能和变形能，下面专门说一下这两项。 • 气相相互作用能ΔE_int：即复合物的电子能量减去各个单体的电子能量，计算用的几何结构都是优化后的复合物结构（单体的结构直接从优化后的复合物结构中抠出来，不再做优化）。 ΔE_int可以通过笔者提出的高效、普适、很流行的sobEDAw方法分解为分子间的静电相互作用、交换互斥作用、轨道相互作用、色散作用，细节这里就不说了，在《使用sobEDA和sobEDAw方法做非常准确、快速、方便、普适的能量分解分析》（//www.umsyar.com/685）里有极其详细的说明。 其中，轨道相互作用部分还可以通过《使用Multiwfn通过ETS-NOCV方法深入分析片段间的轨道相互作用》（//www.umsyar.com/609）介绍的ETS-NOCV方法进一步将其本质讨论得很清楚。 《使用Multiwfn做基于分子力场的能量分解分析》（//www.umsyar.com/442）介绍的笔者提出的EDA-FF能量分解方法可以将静电、交换互斥、色散作用都分解为原子对之间的贡献和原子的贡献，从而能对它们的本质和来源有确切的了解。但注意EDA-FF得到的这三项和sobEDAw给出的有所不同，毕竟计算原理相差很大，前者基于力场计算，而后者基于色散校正的DFT理论计算。另外，色散部分还可以按《使用Multiwfn图形化展现原子对色散能的贡献以及色散密度》（//www.umsyar.com/705）介绍的方式深入地分析。 • 变形能ΔE_deform：在分子形成复合物的整个过程中，单体分子会从孤立状态的结构（单独做优化得到的势能面极小点结构）变成它在复合物中的结构，这会造成分子的能量增加，称为分子的变形能，可以很容易地手动自行计算。结合能中的变形能ΔE_deform是各个分子的变形能的总和，显然可以讨论各个分子各自的贡献。分子的柔性越强、分子间相互作用越强，其变形能倾向于越大。分子变形过程对应各个内坐标的变化，因此变形能可以分解为不同内坐标的贡献来具体讨论。还可以用类似《使用Dushin分解重组能和计算Huang-Rhys因子》（//www.umsyar.com/330）中提到的分解重组能的方式，把变形能分解为极小点结构下的不同振动模式的贡献来讨论。 以上算是把分子间结合自由能里的所有构成做了详细的拆分介绍。把这些原理搞清楚了，在讨论分子间结合强弱的时候就容易分析得比较深入透彻。如果你由于缺乏基本常识，看了上文之后还是不会计算里面的项，推荐通过北京科音初级 培训班（http://www.keinsci.com/KEQC）系统性学习一遍，自然就明白了。 //www.umsyar.com/731 2024-12-12T08:51:00+08:00 深度揭示互为等电子体的苯、无机苯和carborazine的芳香性的显著差异 Deeply revealing the significant differences in the aromaticity of benzene, inorganic benzene and carborazine 文/Sobereva@北京科音  2024-Dec-12 1 前言 众所周知，苯具有极其显著的芳香性特征。它的重要的等电子体borazine（硼吖嗪）也被称为无机苯，具有B3N3H6的化学组成，是把苯的六元碳环上的碳替换为交替相连的氮和硼。无机苯的芳香性虽然早就被很多文章研究过，都指出其芳香性明显不及苯，但不同文章里关于它的芳香性强度说法不一，甚至出入很大，还有的文章使用不严谨的方法对其芳香性进行表征得到了误导性的结论。New J. Chem., 39, 2483 (2015)提出了carborazine（我将它译为碳硼吖嗪）分子，是苯的另一种等电子体，它的结构介于苯与无机苯之间，如下所示，虽然此文指出了它具有芳香性，但对芳香性研究得并不深入和严格。 为了真正全面、透彻阐释苯及其上述两种等电子体的芳香性特征、揭示芳香性差异的内在本质，以及了解carborazine的特点，北京科音自然科学研究中心的卢天和江苏科技大学的刘泽玉等人近期在知名的Chem. Eur. J.（欧洲化学）期刊上发表了专门的研究文章充分研究了上述体系的电子结构和芳香性，非常欢迎大家阅读和引用： Yang Wu, Xiufen Yan, Zeyu Liu,* Tian Lu,* Mengdi Zhao, Jingbo Xu, Jiaojiao Wang, Aromaticity in Isoelectronic Analogues of Benzene, Carborazine and Borazine, from Electronic Structure and Magnetic Property, Chem. Eur. J., 30, e202403369 (2024) DOI: 10.1002/chem.202403369 此文可通过以下链接免费在线阅读： https://onlinelibrary.wiley.com/share/author/PYSFYG86FNTHBBWXITBI?target=10.1002/chem.202403369 这篇文章还被选为了Chem. Eur. J.的封面文章（DOI: 10.1002/chem.202486603）： 下文将对这篇文章的主要内容做深入浅出的介绍，便于读者快速顺利地了解文章的研究结论，并且同时附加一些额外信息，使得读者可以更充分了解研究思想和计算细节，从而能够参考本文的内容更好地做自己的研究。 本文充分利用了强大的Multiwfn程序（//www.umsyar.com/multiwfn）做波函数分析讨论芳香性，是波函数分析方法和Multiwfn程序的很好且十分典型的应用范例。Multiwfn能做的芳香性分析在《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）里有简略介绍，在北京科音 波函数分析与Multiwfn程序培训班（http://www.keinsci.com/WFN）里有多达160页以上的幻灯片极为完整、全面、系统地讲授芳香性的各种概念和全部分析手段，非常推荐参加！ 与本文的研究有密切相关的是《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（//www.umsyar.com/696）介绍的Inorg. Chem., 62, 19986 (2023)一文，十分建议看完下文后认真阅读。 2 苯、无机苯和carborazine的结构和稳定性 文中首先用ωB97XD/def2-TZVP对上述三个体系做了几何优化，这是相当稳的级别，也被用于//www.umsyar.com/carbon_ring.html列举的18碳环及其各种衍生物的研究中，后文的各种分析也都是在这个级别下做的。得到的苯的C-C键长1.387埃和电子衍射实验测定的1.399埃相符极好，得到的无机苯的B-N键长1.426埃也完全在X光衍射测定的1.422-1.429埃范围内。所有三个结构都是纯平面的，苯、无机苯和carborazine分别属于D6h、D3h和C2h点群。 尚未合成的carborazine的稳定性并不得而知，跑从头算动力学（AIMD）是一种有效的检验手段，例如在《18个氮原子组成的环状分子长什么样？一篇文章全面揭示18氮环的特征！》（//www.umsyar.com/725）介绍的文章里通过AIMD直接证明了18氮环非常不稳定。本文对carborazine和borazine用ORCA程序在ωB97XD/6-311G*级别下做了50 ps AIMD，由于AIMD的昂贵耗时使得模拟时间很有限，文中故意用了1000 K的较高温度令势能面的采样尽量充分。文中还确认了当前温度下的AIMD过程的平均总能量是高于极小点结构下的ZPE的，因此当前的模拟不会低估实际的振动基态能量。这种模拟在北京科音高级 培训班（http://www.keinsci.com/KAQC）里有超详细的讲解，在《使用ORCA做从头算动力学(AIMD)的简单例子》（//www.umsyar.com/576）里讲了一些粗略的信息。模拟轨迹多帧叠加图如下所示，每1 ps绘制一次，按照模拟时间颜色按蓝-白-红变化。可见模拟过程中carborazine和无机苯一样始终保持结构稳定，并未解离或异构化，而且结构波动范围相仿佛，这很大程度体现了carborazine像无机苯一样是稳定的物质，在未来很有可能被成功合成出来。 文中还在ωB97XD/def2-TZVP下计算了原子化能，carborazine的原子化能-1200.9 kcal/mol和无机苯的-1225.6 kcal/mol非常接近，进一步体现出carborazine具有显著的稳定性。 值得一提的是carborazine的结构并非是其B2C2N2H6的化学组成的所有的可能异构体中能量最低的，这从文中在SI里给出的各种可能的异构体相对于carborazine的能量可以看出，如下所示，可见让两个碳挨在一起的能量更低。但由于原子的迁移要破坏原先的C-N和C-B键，决速步的能垒必定巨高，因此这样的自发异构化并不会在现实中发生。 3 苯、无机苯和carborazine的电荷分布 原子电荷的概念在《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！》（//www.umsyar.com/714）介绍的笔者的文章里有充分介绍，这是定量衡量体系电荷分布最常用的指标。本文通过流行的ADCH方法以及NPA、Hirshfeld和Hirshfeld-I、CHELPG方法都对前述三个体系计算了原子电荷，结果如下所示。可见虽然不同方法计算的原子电荷存在定量差异，但结论都是共通的，即硼和氮由于电负性一个很小一个极大，导致分别带显著的正电荷和负电荷，而在carborazine中它们的差异远小于在无机苯中，这体现出在硼和氮之间插入的碳可以极大地均衡体系中的电荷分布。另外，carborazine中的碳的原子电荷的大小很小，体现出旁边的氮和硼并没有对碳的净电荷带来明显影响。 《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）介绍了Multiwfn对pi电子的分析功能。本文使用Multiwfn通过Mulliken方法计算了前述体系的非氢原子的pi布居数以考察它们pi电子的分布情况。结果是苯(C=0.995)、无机苯(B=0.429, N=1.558)、carborazine (B=0.586, C=0.985, N=1.411)。此结果清楚地体现了相对于无机苯，carborazine中在硼和氮之间引入的碳原子可以平衡它们的pi电子数的差异。而carborazine中的碳的pi布居数和它在苯中几乎没区别。 为什么要讨论环上的pi电子分布？因为它的分布量和分布均衡性与芳香性有关键的联系。较强的pi芳香性必须同时满足环上的pi电子丰富、pi电子分布较均衡，以及pi电子在环上的离域性较强这三个条件。虽然无机苯的pi电子数和苯一样，都是6个电子，但环上的原子的pi布居数的最小值，即硼的pi布居数，仅有0.429，因此姑且不论其电子的离域性与苯的差距有多大，光凭这点，其芳香性至多也只有苯的约40%。反之，carborazine六元环上硼的pi布居数达到0.586，因此可以认为其芳香性的上限约有苯的60%。 4 苯、无机苯和carborazine的pi轨道 pi分子轨道的特征与pi电子的芳香性有密切联系，所以此文也考察了pi轨道的特征。前述三个体系的HOMO、LUMO轨道都是pi轨道，它们的能级和等值面图如下所示，其它两个占据的pi轨道也一起显示了。可见carborazine的HOMO-LUMO gap比苯和无机苯明显更小，一定程度暗示可能其稳定性相对更弱一些。 Multiwfn可以非常便利地计算轨道成份，见《谈谈轨道成份的计算方法》（//www.umsyar.com/131）。使用Mulliken方法计算的上面图中展示的占据的pi轨道的成份如下所示。可见无机苯的所有pi轨道全都是由氮原子贡献所主导的，这也必定导致无机苯的pi电子几乎都集中且定域在氮原子上，因此在这点上它就不可能有显著的芳香性。而虽然carborazine的HOMO-3和HOMO-5也是以氮原子的贡献为主（因为氮的pz原子轨道的能量比碳和硼的都要低，故偏向于对能量较低的分子轨道贡献较多），但还不至于起到绝对的主导性，而其HOMO的主要贡献则来自硼和碳。由于carborazine的硼、碳、氮都同时明显参与了pi占据轨道，它势必有比无机苯强得多的六中心pi电子离域。 5 苯、无机苯和carborazine的成键特征 键级的概念在《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）中做了简要介绍。文中用Multiwfn对上述三个体系计算了Mayer键级，它体现两个原子间等效共享的电子对儿数，结果如下表左侧所示，并且还按照《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）所述的方法将sigma和pi电子的贡献分别给出了。可见虽然B-N、C-N、C-B的pi键级都比C-C键的更小，但数值也都不低，所以它们的pi作用都是不容忽视的，也因此切勿以为B-N键只能算作纯sigma键。carborazine六元环上的各个键的pi键级都大于无机苯的B-N键，这点已在一定程度上能体现出carborazine的整体pi共轭是强于无机苯的。 文中还对前述体系用Multiwfn做了流行的atom-in-molecules (AIM)拓扑分析，以从键临界点（BCP）处的实空间函数角度考察体系的成键特征，结果如上表右侧所示。AIM的相关知识和Multiwfn里的操作见AIM学习资料和重要文献合集（http://bbs.keinsci.com/thread-362-1-1.html）里提到的内容。共价键在BCP处的能量密度（H）应当为负，并且eta指数（η）大于1是最典型共价键的特征。当前的结果体现出苯的C-C键和carborazine的C-N键都具有最典型的共价键特征，而carborazine和无机苯的N-B的共价性相对最低、离子性相对最强，因此H没那么负且η显著小于1，这无疑是硼和氮的巨大电负性差异所带来的。 要注意切勿把B-N当做离子键。有些文章盲目使用BCP位置的电子密度拉普拉斯函数（▽2ρ）的正负作为共价键的判据，如Inorg. Chim. Acta., 360, 619 (2007)以这个为理由说B-N是离子键，这是很误导的。之前我在《AIM键临界点处电子密度拉普拉斯值符号判断相互作用类型失败原因的图形分析》（//www.umsyar.com/161）中专门说过类似情况。对B-N键，如前所示Mayer键级显著大于1，而且BCP处的能量密度明显为负，就已经充分体现出这属于共价键了。文中更进一步在补充材料里还用Multiwfn绘制了电子定域化函数（ELF）和变形密度图，分别如下图左侧和右侧所示，可见B-N之间有显著的电子高定域性区域，并且从孤立状态变成分子过程中在B-N键上有显著的电子密度增加（图中红色实线为正值部分），这都进一步证明B-N键本质上是共价键，只不过极性很强而已。关于变形密度的知识见《使用Multiwfn作电子密度差图》（//www.umsyar.com/113）和《通过价层电子密度分析展现分子电子结构》（http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB201709252），关于ELF的知识见《ELF综述和重要文献小合集》（http://bbs.keinsci.com/thread-2100-1-1.html）和《电子定域化函数的含义与函数形式》（http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB20112786），绘制方法参考Multiwfn手册4.4节，在 波函数分析与Multiwfn程序培训班（http://www.keinsci.com/WFN）里有更丰富的例子、更细致的讲解和各种技巧的传授。 6 苯、无机苯和carborazine的电子离域特征 如《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）和笔者的Theor. Chem. Acc., 139, 25 (2020)所提到的，pi电子的局域轨道指示函数（LOL），即LOL-pi，对于展现pi电子的离域性极为直观。对前述三个体系用Multiwfn计算并用VMD绘制的LOL-pi等值面图如下所示，可见苯的六中心共轭作用最强，等值面均匀、理想地贯穿六元环。虽然carborazine和无机苯的LOL-pi等值面都没有分布得很均匀，特别是在B-N键处受到了阻碍，算是电子共轭的“瓶颈”，但是carborazine的每个B-C-N单元上的等值面较为均匀、连贯，这是无机苯不具有的，因此从LOL-pi展现的电子离域性图上可以期望carborazine有明显更强的芳香性。 文中还绘制了EDDB_P-π等值面图，详见原文。LOL-pi展现的是不同区域电子的离域程度，而EDDB_P-π展现的是不同区域离域的pi电子量，二者体现的信息存在互补性。EDDB_P-π等值面图及其布居分析展现出无论是carborazine还是无机苯，硼上的离域的pi电子都比环上的其它原子更少，原因之一是硼上的pi电子数本来就非常少，这在前面的pi布居分析中已经体现了。相比之下，carborazine的硼的离域的pi电子为0.360，明显多于无机苯的0.184，这进一步体现出carborazine具有更强的pi芳香性。并且0.360这个值是苯的碳原子的离域pi电子数0.885的41%，这说明carborazine至少有苯41%的芳香性。 进一步，文章还使用Multiwfn计算了多中心键级（MCBO），以及sigma+core和pi电子各自的贡献量，相关知识参看《使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键》（//www.umsyar.com/138）和《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）。结果如下所示，可见六中心电子离域性是苯 > carborazine > 无机苯，差异非常鲜明，这和上述其它分析给出的结论非常一致。从MCBO的数值大小来说，carborazine具有苯的大约一半的芳香性，而无机苯则连carborazine芳香性的一半都不到。 7 苯、无机苯和carborazine的磁屏蔽特征 有大量的方法从体系对外磁场的屏蔽角度衡量芳香性，参考笔者写过的《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）、《使用Multiwfn绘制一维NICS曲线并通过积分衡量芳香性》（//www.umsyar.com/681）、《使用Multiwfn计算FiPC-NICS芳香性指数》（//www.umsyar.com/724）、《通过Multiwfn绘制等化学屏蔽表面(ICSS)研究芳香性》（//www.umsyar.com/216）、《使用Multiwfn巨方便地绘制二维NICS平面图考察芳香性》（//www.umsyar.com/682），在《 波函数分析与Multiwfn程序培训班》（http://www.keinsci.com/WFN）里有更全面系统的讲授。本文从磁屏蔽角度对苯、无机苯和carborazine的芳香性特征进行了定量的考察以及直观的展现。 NICS及其变体在定量衡量芳香性强弱方面极为常用，文中对苯、无机苯和carborazine计算了NICS最理想的形式之一NICS(1)ZZ，并通过《基于Gaussian的NMR=CSGT任务得到各个轨道对NICS贡献的方法》（//www.umsyar.com/670）介绍的做法分解成了sigma+core和pi电子各自的贡献，数值越负说明芳香性越强。由数据可见，carborazine的pi芳香性虽然不及苯，但也相当显著，而无机苯的芳香性极为微弱，接近非芳香性。sigma电子对芳香性的贡献可忽略不计。 如《使用Multiwfn绘制一维NICS曲线并通过积分衡量芳香性》所述，对NICS_ZZ在垂直于环的方向上积分，即∫NICS_ZZ或称integrated NICS (INICS)，是比NICS(1)ZZ还更为严格的基于磁屏蔽特征衡量芳香性的方法，因为它的结果不依赖于计算磁屏蔽张量的位置的选取。从上面表格的数据可见∫NICS_ZZ的结论与NICS(1)ZZ相一致，更进一步体现了之前结论的可靠性。∫NICS_ZZ计算过程中利用到的穿过环中心且垂直于环方向上的NICS_ZZ曲线可以绘制出来并进行直观对比，如下所示，横坐标是相对于环中心（零点）的位置。根据pi电子的NICS_ZZ曲线，可以清楚地看出carborazine的曲线与苯比较接近，这将carborazine显著的pi芳香性特征体现得非常确切，而无机苯的pi芳香性相比之下显得微乎其微，其曲线在全范围都很接近0。三个体系的sigma电子的NICS_ZZ曲线几乎完全重合，很严格地体现出它们的sigma芳香性基本没有任何区别。 下图左列是Multiwfn结合VMD绘制的ICSS_ZZ等值面图，等值面数值为3 ppm，橙色和深蓝色分别是对外磁场屏蔽和去屏蔽区域。从这种图上能非常完整地了解体系的磁屏蔽特征。可见carborazine和苯一样都是在分子外围形成了环状的去屏蔽区域，又一次体现出了carborazine类高度似于苯的显著的芳香性，而无机苯的去屏蔽区域等值面则明显显得不够连续，体现出其芳香性远比carborazine弱得多。 上图还给出了Multiwfn直接绘制的填色等值线图，中列和右列分别是分子平面以及分子的侧截面，颜色越红和越蓝分别对应相应位置的磁屏蔽和去屏蔽效应越强。可以看到carborazine在分子中心是磁屏蔽的，颜色和苯一样是偏红的，而borazine在分子中心却是去屏蔽的，颜色偏蓝，这充分反映了carborazine和苯在芳香性上是类似物，而无机苯可近似归为非芳香性分子一类。 8 苯、无机苯和carborazine的感生电流 在外磁场作用下在环状区域内离域的电子可以形成感生环电流，这是衡量电子离域性和芳香性的非常直观的方法，其中AICD展现方式尤为常见，详见《使用AICD 2.0绘制磁感应电流图》（//www.umsyar.com/294）。本文对苯、无机苯和carborazine分别绘制了总的、sigma电子的和pi电子的AICD图，如下所示，从AICD(pi)图中箭头标注的感生电流上可以清楚地看出carborazine形成了和苯一样的贯穿整个六元环的显著的感生电流，而无机苯则完全不具备这样的特征。AICD图进一步有力地体现出，carborazine和苯一样都属于芳香性物质，而无机苯是近乎非芳香性的，这和ICSS分析的结论一致，彼此相互印证。 此文还使用《使用SYSMOIC程序绘制磁感生电流图和计算键电流强度》（//www.umsyar.com/702）中介绍的做法用CTOCD-DZ2方法绘制了分子平面上方1 Bohr处的感生电流图，如下所示，灰、红、蓝分别对应碳、硼、氮原子。由图可见CTOCD-DZ2方法展现的感生电流图和pi电子的AICD图的现象很类似，而这种平面图看得更清楚一些。此图体现carborazine分子平面上方pi区域的感生电流的强度虽不像苯那么均匀，在硼上方相对弱一些，但终究还是形成了明显的全局感生电流。而无机苯的感生电流则主要在氮原子周围打转，体现出pi电子主要都是定域在各个氮的局部原子空间里。 为了便于定量对比，此文还计算了六元环上各种键的电流强度（BCS）值，即在键中心位置垂直于键的平面上的感生电流的积分，结果如上图标注的数值所示，单位为nA/T。可见carborazine的感生电流强度与苯相差并不悬殊，而无机苯则相当微弱。这从定量角度上进一步验证了前面基于肉眼观看感生电流图所做出的结论。 值得一提的是，在《深入理解分子轨道对磁感生电流的贡献》（//www.umsyar.com/703）中我给了计算感生电流的具体公式，假定只有HOMO-LUMO跃迁对感生电流有显著贡献，那么感生电流大小与HOMO-LUMO gap之间会有反比关系。在前面说过，carborazine的HOMO-LUMO gap明显小于苯和无机苯，因此基于磁属性和感生电流得到的关于carborazine的芳香性的强度会有一定程度的高估。但即便如此也不影响它具有明显芳香性的结论，而且前面从多中心键级、LOL-pi等分析角度也全都体现了carborazine有远强于无机苯的芳香性。 9 carborazine中的碳原子对硼和氮原子桥联的重要意义 前面从各个角度极为充分、全面、严格地论证了carborazine具有无机苯不具备的显著的芳香性，而二者的差异仅在于前者给后者的氮和硼之间引入了碳原子起到桥联作用。为什么碳原子的桥联对于提升芳香性有如此关键性的影响？此文将Multiwfn、NBO、VMD相结合绘制了前述三个体系的环上的各原子的对应2pz原子轨道的预正交自然原子轨道（PNAO）等值面的叠加图，过程可参照《使用Multiwfn绘制NBO及相关轨道》（//www.umsyar.com/134）举一反三，如下所示。可见2pz的空间延展程度是硼 > 碳 > 氮，这和电负性增大的顺序相一致。氮和硼的轨道重叠颇差，电子明显难以离域过去。而给二者之间插入的碳原子，增加离域路径上的原子轨道重叠程度、降低相连原子的电负性不平衡性，明显可以起到了帮助电子离域过去的“中介”，从而使得carborazine能够具有无机苯不具备的显著芳香性。之前笔者在《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（//www.umsyar.com/696）介绍的文章中也使用了类似的角度进行分析。 文章里还研究了其它4种苯的等电子体，以考察碳的数目和位置对芳香性的影响。它们的结构如下所示，NICS(1)ZZ也标在了图上，注意carborazine的值如前所述是-23.1 ppm。通过对比可见，环上的碳原子数多不意味着芳香性就强，例如下图第1、3个分子虽然有4个碳原子，但单从NICS(1)ZZ来看其芳香性还没只有两个碳的carborazine强。下图第4个分子和carborazine的碳原子数一样，但碳没有插入在硼和氮之间，芳香性也远不如carborazine强。而只有下图第2个分子，碳桥联了氮和硼，同时碳原子数比carborazine还更多，这才带来了稍微更强的芳香性。这个对比体现出对只有硼、氮的环上引入碳，而且能起到对它们的桥联作用时，才能显著提升芳香性。 10 总结 本文介绍的Chem. Eur. J. 2024, e202403369一文十分全面透彻探讨了互为等电子体的苯、无机苯（borazine，B3N3H6）和尚未实验合成的carborazine（B2N2C2H6）的电子结构和芳香性的差异。此文的重点是综合运用各种分析手段，充分论证了含有碳原子桥联硼、氮原子特征的carborazine具有明显的芳香性（尽管强度一定程度弱于苯），同时还严格地证明了无机苯的芳香性极弱，甚至很大程度上体现出非芳香性分子的特征。文章对此现象的本质从电子结构层面进行了详细的分析和解释。此文的研究工作对于化学家们更好地了解由碳、氮、硼构成的物质的电子离域性和芳香性有显著的意义。 值得一提的是，此文的研究也是使用Multiwfn等程序做波函数分析深入剖析化学物质特征的非常典型的文章，可以作为波函数分析的很好的范例作为参考。此研究涉及的研究方法在前面提及的各种博文里有一些粗略、初步的介绍，通过 波函数分析与Multiwfn程序培训班（http://www.keinsci.com/WFN）可以一次性全面深入透彻掌握文中用到的一切波函数分析方法的背景知识和使用具体程序计算的方法，从而能顺利地复现此文的研究并举一反三探索更多的体系。 //www.umsyar.com/730 2024-12-08T05:19:00+08:00 Molclus FAQ 文/Sobereva@北京科音  2024-Dec-8 笔者开发的免费的分子构象和团簇构型搜索程序Molclus（http://www.keinsci.com/research/molclus.html）如今已经非常流行，在网上经常有人提问此程序的使用问题。这里我把一些比较常见的问题进行统一回答。日后还会对本文内容有进一步追加。 Q1：Molclus是干什么？ A：如果你没有任何 背景知识，先把《辨析计算化学中的任务类型和理论方法》（//www.umsyar.com/680）看了。多原子分子、原子/分子团簇体系的势能面上一般有不止一个极小点，直接基于自己猜的初始结构去用Gaussian等程序做几何优化，极可能会优化到不是能量最低的极小点，这会导致之后做的分析讨论没实际意义、结果有误导性。Molclus可以帮助你获得能量最低的极小点，即全局最小点的结构，这是体系最有意义的结构。 此外，柔性分子、团簇的势能面上通常有一批能量都很低、能量相差不大的极小点结构，在实际研究的温度下都有一定出现比例（参见《根据Boltzmann分布计算分子不同构象所占比例》//www.umsyar.com/165），计算构象权重的光谱（参见《使用Multiwfn绘制构象权重平均的光谱）//www.umsyar.com/383和《使用Multiwfn绘制NMR谱》//www.umsyar.com/565）以及计算构象权重的热力学量（参见《使用Shermo结合 程序方便地计算分子的各种热力学数据》//www.umsyar.com/552）需要得到这些低能量极小点的结构和出现百分比，Molclus也可以很好地实现这个目的。 Q2：怎么引用Molclus程序？ A：计算化学程序的引用在《写计算化学文章时引用理论方法、基组、程序时应注意的问题》（//www.umsyar.com/370）里明确强调了。使用Molclus发文章（包括给别人代算）必须按照主页http://www.keinsci.com/research/molclus.html上明确写明的方式引用，引用必须在正文里体现。Molclus完全免费，按要求正确引用Molclus是对此程序开发和维护的最好的支持。 Q3：使用Molclus程序遇到问题如何求助？ A：在计算化学公社论坛（http://bbs.keinsci.com）的“ ”版块发帖，发帖时选择Molclus分类。也可以在笔者创建的思想家公社QQ群里提问，群介绍和群号见//www.umsyar.com/QQrule.html。我看到后都会尽量回答。提问时必须把使用细节交代得尽可能清楚，避免含糊其辞或引起歧义，注意看《在网上求助计算化学问题的时候必须把问题描述得详细、具体、准确、清楚》（//www.umsyar.com/620）。 Q4：Molclus程序怎么上手？ A：Molclus主页http://www.keinsci.com/research/molclus.html上明确给了具体教程，写得极为详细，认真一个字一个字阅读、照着操作即可。 需要注意的是，Molclus虽然非常易于使用、灵活方便，但不代表这是完全黑箱程序，使用者不能一点计算化学常识都没有。Molclus需要调用第三方程序来做单点、几何优化、振动分析等任务。Molclus调用哪个程序，用户就应当有相应程序的最基本常识。例如使用Molclus调用Gaussian做优化和单点能计算是最常涉及的情况，因此用户就得懂Gaussian最基本的使用，要不然连Gaussian怎么安装和运行、Gaussian模板文件里怎么写关键词、最常见的Gaussian报错都不知道怎么解决可不行。而且仅当你具备基本的Gaussian计算常识的情况下，才能知道怎么针对当前的体系根据精度要求和计算资源选择合适的级别用Molclus做构型或构象搜索。不具备相关知识的话非常建议通过北京科音初级 培训班（http://www.keinsci.com/KEQC）快速且系统地学习。 Q5：Molclus做构型/构象搜索的基本流程是什么？ Molclus的最简单的使用过程是： (1)用Molclus的gentor或genmer或第三方程序产生一批初始构型/构象记录在traj.xyz文件里，这是多帧的xyz格式的文件，见《谈谈记录化学体系结构的xyz文件》（//www.umsyar.com/477） (2)用Molclus对当前目录下的traj.xyz里的构型/构象依次做优化，得到的结构和能量都会储存在isomers.xyz里 (3)用Molclus里的isostat工具对isomers.xyz中的结构去重、按能量排序，得到cluster.xyz，里面第一个就是全局能量最小点结构，后面的是能量次低的结构 更准确而且也是实际中特别常用的流程如下，这里以普通有机体系为例 (1)如上所述先产生一批初始构型/构象 (2)用Molclus结合xtb程序在很便宜但也较糙的GFN-xTB级别下对这些构型/构象依次做很快速的优化（对柔性大体系这一步可能要对几百甚至上千个结构做计算) (3)用isostat对优化后的结构去重、按能量排序得到cluster.xyz (4)将cluster.xyz改名为traj.xyz，用Molclus调用Gaussian在精度够用而且也不贵的B3LYP-D3(BJ)/6-31G*级别下对之前能量最低的几kcal/mol的结构（通常在几十个以内）进行更准确的优化以及做振动分析（后者是为了检验虚频情况以及获得自由能热校正量），并且之后自动在M06-2X/def2-TZVP级别下做较高精度的单点能计算。此时isomers.xyz里存的是各个结构及其较准确的自由能 (5)用isostat对isomers.xyz里的结构去重、按自由能排序得到cluster.xyz，并且可以输入温度让isostat直接算出来各个结构的Boltzmann出现比率 记住如果是做溶液情况下的构型/构象搜索，单点计算必须带溶剂模型（建议用SMD），优化和振动分析最好也带上以图稳妥（建议用IEFPCM）。如果嫌振动分析贵也可以不做，最后isostat会将按照电子能量而非自由能排序，给出的Boltzmann分布的准确性可能会差不少。 没有哪种构型/构象搜索流程对所有情况都是最理想的，精度和效率不可能同时特别理想，上面说的只是对一般有机体系的精度和效率比较均衡的搜索流程。使用者可根据实际情况随机应变。具体操作细节在Molclus官网上的教程里都有体现。 上述提到的Gaussian是收费的，因版权问题实在不方便用Gaussian的话也可以让Molclus调用免费的ORCA 程序来代替，再加上xtb和Molclus都是免费的，因此利用Molclus做构象搜索可以全程免费，只不过ORCA在优化+振动分析方面的整体稳健性方面不如Gaussian。用ORCA的话，上述优化+振动分析部分建议用B97-3c，单点计算部分建议用wB97M-V/def2-TZVP结合gCP形式的BSSE校正。ORCA和xtb的使用在《北京科音高级 培训班》（http://www.keinsci.com/KAQC）里有极其全面的介绍。 Q6：哪种方式产生初始构型/构象最适合我？ A：不同方式产生初始构型/构象各有长处，在Molclus官网上写明了，如下所示。有不清楚的可以在公社论坛或QQ群里向我咨询 • Molclus自带的genmer程序结合Molclus：用于分子团簇或原子团簇构型搜索，使用简便。对于分子团簇搜索这通常是首选。唯一局限性是对分子的构象考虑不够充分，因为genmer将分子当成刚性来堆积，因此不适合单分子柔性很强的情况，这种情况更适合后面说的做动力学程序结合Molclus。相关介绍和示例看《genmer：生成团簇初始构型结合molclus做团簇结构搜索的超便捷工具》（http://bbs.keinsci.com/thread-2369-1-1.html） • Molclus自带的gentor程序结合Molclus：用于分子构象搜索。使用简便，但不支持环状区域构象搜索（环状区域是指诸如环己烷这样有多种构象的环状区域）。对于可旋转的键不很多的分子的构象搜索这是首选，可控性最强。相关介绍和示例看《gentor：扫描方式做分子构象搜索的便捷工具》（http://bbs.keinsci.com/thread-2388-1-1.html） • 第三方的Confab或Frog2构象生成程序结合Molclus：用于分子构象搜索。使用最为傻瓜化，但不支持环状区域构象搜索、不支持普通有机分子以外的情况。相关介绍和示例看《将Confab或Frog2与Molclus联用对有机体系做构象搜索》http://bbs.keinsci.com/thread-20063-1-1.html） • xtb程序跑动力学轨迹结合Molclus：普适性极强，用于构象搜索、原子/分子团簇构型搜索皆可，使用容易。但由于模拟时间长度有限，有动力学采样不充分导致遗漏构型/构象的风险（风险随动力学模拟时间增长而减小）。相关介绍和示例看《使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索》（http://bbs.keinsci.com/thread-16255-1-1.html） • GROMACS等经典力场分子动力学程序结合Molclus：用于分子团簇搜索和分子构象搜索。使用稍繁琐，因为得创建拓扑文件，且被动力学模拟的体系必须有恰当的力场，好处是动力学模拟耗时极低，因此动力学采样不充分导致遗漏构型/构象的风险很小。相关介绍和示例看《使用molclus程序做团簇构型搜索和分子构象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）。不了解GROMACS使用的话推荐通过《北京科音分子动力学与GROMACS培训班》（http://www.keinsci.com/KGMX）学习相关使用知识。 Q7：为什么运行Molclus的时候闪退/崩溃？ 只有两个可能： (1)被Molclus所读取的traj.xyz的格式有问题导致程序崩溃。要按照《谈谈记录化学体系结构的xyz文件》（//www.umsyar.com/477）说的用文本编辑器检查其内容看是否符合标准格式。也可以试图往免费的VMD可视化程序里载入，如果载入不了也说明格式有问题 (2)被Molclus所读取的settings.ini文件有问题。典型情况是用户修改里面设置的时候把格式改乱了，导致无法被读取。也有的用户干了多余的事，自作主张地把他觉得不需要的内容给删掉了，导致Molclus没法读取要读取的内容。还有的用户直接把老版本Molclus的settings.ini放到新版本Molclus目录下用，由于新、老版本此文件里面包含的选项不同，导致新版本Molclus无法正常读取之（所以不要干这种事，要改什么设置一定要在当前版本Molclus自带的settings.ini基础上按照原本的格式改） 死活搞不明白的话，到计算化学公社论坛发帖，把整个Molclus目录压缩后作为附件上传。 Q8：为什么运行Molclus自带的genmer或gentor后程序闪退/崩溃？ A：要么是涉及的xyz文件格式有问题，自行检查；要么是genmer的genmer.ini、gentor的gentor.ini里的格式没写对或者写的内容有bug。死活搞不明白的话，到计算化学公社论坛发帖，把整个genmer或gentor目录压缩后作为附件上传。 Q9：为什么产生的isomers.xyz内容为空？ A：这说明连一个结构都没跑成功。一定要仔细看Molclus运行窗口中的提示，Molclus在执行什么、运行状态如何，都提示得清清楚楚。导致此问题常见的原因如下，是哪种情况根据Molclus输出的提示就能很容易地判断（实在看不懂的话在提问时把完整截图贴出来） (1)Molclus调用Gaussian、xtb等第三方程序失败。通常是第三方程序在机子里没装好，或者settings.ini里调用它们的命令没写对 (2)连一个任务都还没跑完。如果你确信你的计算设置没有离谱的地方，就继续等着 (3)所有任务都运行失败了。如果把settings.ini里的ibkout设为1，每个任务的被调用的第三方程序的输出文件都会被保留，应当根据里面的信息结合相应程序的使用常识判断当前是什么原因所致。常见的问题包括：模板文件里关键词写错了、净电荷或自旋多重度设错了、SCF/几何优化不收敛。 Q10：怎么Molclus老也算不完？ A：如果根据Molclus的提示能确信任务正在算着，这种情况无非就是以下原因 (1)CPU忒烂导致耗时太高。一些初学者对耗时没基本常识，在个位数核心的CPU上就想让Molclus调用Gaussian或ORCA用DFT优化上百原子体系、对其做振动分析，耗时不高才怪。牵扯到大批量上百原子体系的DFT计算都应该有个像样的性能较好的服务器，要不然做构象搜索可能得跑N天甚至一周以上。也可能使用者的CPU并不差，但没恰当并行，或者非要用Windows版的Gaussian，导致CPU性能无法充分发挥 (2)用的计算级别过于昂贵。诸如让Molclus调用Gaussian对挺大的体系用def2-TZVP这么贵的基组做优化、对其用双杂化泛函做单点计算甚至优化，耗时不高才怪 (3)遇到了难收敛。应根据被调用的程序的使用常识，根据其临时输出文件内容判断当前是什么情况。Molclus调用第三方程序用的命令在屏幕上提示得清清楚楚，被调用的程序产生的临时输出文件是哪个也体现得明明白白，比如Gaussian就是gau.out。对于调用Gaussian做几何优化，若一个任务老跑不完，显然应该想到用GaussView打开gau.out，看看是遇到了震荡难收敛还是什么其它情况，以及弄清楚当前跑多少步了、是不是每一步的耗时太高，等等。如果真遇上几何优化难收敛，Gaussian用户可以参考《 计算中帮助几何优化收敛的常用方法》（//www.umsyar.com/164）在模板文件里加上适当的帮助收敛的关键词。对于Molclus调用ORCA也是类似的处理思路 Q11：Molclus对各个结构优化时有个别没成功，结果可以用么？ A：我只能说没优化成功的任务占总结构数目比例越小，遗漏重要结构的概率越低。如果你用ibkout=3把失败的任务的输出文件都备份了出来，可以看看其优化最后一帧的结构和isostat排序后能量最低一批结构的差异，如果看上去优化失败的结构继续跑下去也不太可能收敛到能量较低的结构去，那这样的结构就可以舍弃。而如果你怀疑它继续优化有可能能得到一个不可忽视的能量很低的结构，那就把他回炉结合适当的关键词重新优化。值得注意的是settings.ini里可以通过ngeom指定对traj.xyz里哪些帧进行处理，iappend设1可以把新得到的结果在之前的isomers.xyz上续写，因此可以方便地将回炉的结果累加到之前的上面。 Q12：Molclus能在Windows下用么？ A：能是能，但在Linux下用更理想。xtb程序虽然有Windows版，但有些bug，Windows下使用时容易莫名其妙中断。因此笔者在Windows主机下用Molclus调用xtb时都是在VMware虚拟机里装的Linux客户机里用。而且Gaussian的Windows版的计算效率不及Linux版，尤其是CPU核数较多时相差很悬殊，这点在《Gaussian的安装方法及运行时的相关问题》（//www.umsyar.com/439）里专门说了。 Q13：为什么Windows下用Molclus调用Gaussian计算时蹦出Gaussian图形窗口后卡住不动？ A：那是因为你把settings.ini里gaussian_path的路径设成了比如g16w.exe（对于Gaussian 16来说），带w后缀的只是Gaussian 16图形界面这层壳，g16.exe才是真正的Gaussian 16程序本体，gaussian_path必须指定成它的路径。 Q14：为什么Molclus调用xtb不能运行，屏幕上提示诸如：'xtb'不是内部或外部命令，也不是可运行的程序 A：显然是xtb程序没在机子里装好，倘若装好了就直接可以用xtb命令调用xtb程序，自然不会有这种提示。仔细看《将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析》（//www.umsyar.com/421）了解xtb怎么安装。 Q15：为什么Windows下用Molclus调用Gaussian计算时提示No executable for file l1.exe？ A：没有把Gaussian的文件夹加入GAUSS_EXEDIR环境变量所致。这点在《使用molclus程序做团簇构型搜索和分子构象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）里专门说了，在文中搜GAUSS_EXEDIR。或者你加了但路径没写对，也可能环境变量尚未生效（尝试重启机子再试）。 Q16：genmer产生构型时怎么考虑单分子的构象？是否需要先对单分子做构象搜索？ A：genmer把每个分子当做刚性看待，所以没法考虑。如果单分子的柔性较强，因此形成复合物的时候构象可能发生极大的变化，那么建议使用分子动力学方式产生团簇初始结构的traj.xyz文件。对于柔性分子先对单分子做构象搜索意义不大，毕竟分子间相互作用可能导致复合物中的分子的构象相对于在孤立状态下发生很大变化。 Q17：genmer产生多少个初始构型合适？ A：产生的初始构型越多，显然经过批量优化后遗漏掉能量最低构型的可能性就越低。对于分子较小、分子数较少的团簇的构型搜索，比如水的四聚体，初始构型产生大几十个就足够了。分子越大、分子数越多的团簇，由于势能面维度更高、极小点更多，原理上应该产生越多的初始构型。鉴于调用xtb程序做GFN-xTB级别的初筛耗时非常低（调用GFN-FF分子力场的话耗时更低），优化大几百、上千个都完全不是难事，因此在计算量可以接受的情况下对较大的体系产生的初始构型原则上多多益善。也可以根据之前得到的结果凭直觉判断是不是还有遗漏能量更低构型的可能，如果怀疑有这种可能性，可以让genmer再产生一批构型，让Molclus再次批量优化并续写到之前的isomers.xyz上。 Q18：gentor该旋转哪些键？ A：旋转的键越多，产生的初始构象数显著越多，导致构象搜索耗时明显越高，因此不能太多。只有柔性的、现实环境下易于旋转的键才应该在gentor里要求旋转来产生各种初始构象，这凭借基本常识就能判断。现实中不易旋转的键有：(1)有单套pi共轭作用的键，比如肽键就有pi共轭使其C-N键的旋转势垒巨高、在常温下不可能自发旋转 (2)旋转后会严重破坏之前存在的显著的吸引性的分子内弱相互作用，比如导致破坏好多氢键 (3)旋转后会出现严重的位阻作用而且无法靠旋转其它部分回避。 也不是所有易于旋转的键都应当纳入考虑，这凭常识就能判断，比如甲基虽然在常温下极易旋转，但旋转前后的结构是相同的，自然旋转它毫无意义。 如果你实在缺乏常识，索性用confab产生初始构象，但不像gentor那样有可控性。 Q19：isostat的能量和结构去重阈值怎么设置合适？ A：若无特殊情况就用默认设置即可。阈值多大合适和计算级别的精度无关，而取决于几何优化的收敛精度。默认的去重阈值适合几何优化收敛限偏严的情况（如Gaussian默认的收敛限）。当你用了较宽松的几何优化收敛限，去重阈值也需要相应地放宽。 如果发现得到的cluster.xyz中有肉眼看上去特别接近的结构，说明去重阈值偏小，导致本应该只保留其中一个的情况下多个结构都被保留了，此时可以稍微加大阈值（如变成默认的两倍）重新尝试。但阈值也不能太大，否则能量和结构比较接近的对应实际不同极小点的构象会只保留一个。 Q20：isostat判断重复结构的原理是什么？ A：在《使用molclus程序做团簇构型搜索和分子构象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）里专门说了，仔细看相关部分。 Q21：怎么对过渡态做构象搜索？ A：不小的体系的化学反应只发生在体系的一部分区域，可以叫反应区域，整个反应路径中只有这部分区域的原子涉及的内坐标有特别显著变化、牵扯化学键的改变，基于这点可以做过渡态的构象搜索。首先照常搜索过渡态，然后用gentor对非反应区域做各种二面角的旋转产生traj.xyz，在Molclus做批量优化时用settings.ini里的freeze设置将反应区域的原子冻住，之后照常用Molclus做剩下的构象搜索步骤即可。当然，反应区域和非反应区域会存在一定耦合，非反应区域处于不同构象时反应区域的结构也会多多少少有所不同。因此为了严谨起见，建议以上述构象搜索得到的能量最低一批结构作为初猜再次优化过渡态，并取其中能量最低的过渡态。 Q22：为什么用OpenBabel通过Confab方法产生一批构象时失败？ A：和OpenBabel版本可能有关，读者请下载《将Confab或Frog2与Molclus联用对有机体系做构象搜索》（http://bbs.keinsci.com/thread-20063-1-1.html）里直接提供的OpenBabel版本，并且先确保对帖子里的体系能正常产生构象然后再做自己的体系。并且要参考帖子确保输入的OpenBabel命令合理，包括选项的大小写。OpenBabel不要装到默认的C盘的目录下，免得由于操作系统可能的权限限制导致OpenBabel无法成功创建新文件。也要确保提供给OpenBabel的mol2文件格式合理，参考《谈谈记录化学体系结构的mol2文件》（//www.umsyar.com/655）。 Q23：Molclus怎么在超算上用？ A：笔者不用超算，也并未特意考虑Molclus在超算上用的情况。但Molclus是完全可以在超算上跑的，参考以下帖子 《consearch：一键提交slurm的molclus构象搜索脚本》（http://bbs.keinsci.com/thread-43932-1-1.html） 《分享molclus配套PBS作业提交脚本》（http://bbs.keinsci.com/thread-15598-1-1.html） 《分解合并Molclus任务的小脚本splitMC》（http://bbs.keinsci.com/thread-14361-1-1.html） //www.umsyar.com/727 2024-10-28T20:02:00+08:00 全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！ Comprehensively revealing the unique pi-pi interactions between different cyclocarbon and fullerene! 文/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版，对其它版本不一定适用，请结合相应版本程序提示和手册随机应变。 《深度揭示互为等电子体的苯、无机苯和carborazine的芳香性的显著差异》（//www.umsyar.com/731》介绍的笔者的Chem. Eur. J., 30, e202403369 (2024)文章中使用本文介绍的方法清晰直观地展现和讨论了苯、无机苯和carborazine芳香性的差异，非常推荐阅读并将此文作为例子引用。 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）和本文有密切关联，非常建议阅读。