//www.umsyar.com/category/Multiwfn Multiwfn //www.umsyar.com/759 2025-12-08T23:42:00+08:00 使用amIGM方法图形化直观展现动态过程中的平均弱相互作用 Using amIGM method to graphically and intuitively represent averaged weak interactions in a dynamic process 文/Sobereva@北京科音   2025-Dec-9 0 前言 mIGM是笔者在Struct. Bond., 190, 297 (2026) DOI: 10.1007/430_2025_95中提出的能够非常快速、直观地展现自定义片段间相互作用的方法，实用性极强，在《使用mIGM方法基于几何结构快速图形化展现弱相互作用》（//www.umsyar.com/755）中做了专门介绍并演示了怎么通过Multiwfn程序实现，如果没看过此文的话一定要先看看。在同一篇论文里，基于mIGM方法笔者还进一步提出了amIGM方法，全称为averaged mIGM，此方法可以展现动态环境中的片段间的平均弱相互作用。amIGM方法远远比我以前在《使用Multiwfn研究分子动力学中的弱相互作用》（//www.umsyar.com/186）中介绍的与它用处类似的aNCI方法好用，因此aNCI方法完全可以弃了！ 下面第1节将介绍amIGM并与其它方法对比，第2节介绍一些计算要点，第3节将会示例如何通过Multiwfn程序做amIGM分析重现amIGM原文中的图。读者请务必使用2025-Dec-8及以后更新的Multiwfn版本。如果你对Multiwfn不了解，强烈建议看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。 使用Multiwfn做amIGM分析发表文章时应同时引用上面提到的amIGM原文以及Multiwfn启动后提示的程序原文，也推荐一起引用我发表的包含amIGM在内的相互作用可视化分析方法综述Angew. Chem. Int. Ed., 137, e202504895 (2025)，此文的介绍见//www.umsyar.com/746。 1 amIGM方法介绍 1.1 amIGM方法的特点 这里假定读者已经了解mIGM方法，也看过了《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（//www.umsyar.com/667）和《Angew. Chem.上发表了全面介绍各种共价和非共价相互作用可视化分析方法的综述》（//www.umsyar.com/746）这两篇综述以及《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）了解了sign(λ2)ρ着色的δg_inter等值面是怎么一回事。mIGM、IGMH、IRI、NCI等方法都是对于单一几何结构做的分析，然而现实中分子是在不断运动的，因此分子间相互作用也是随时间变化的，有很多情况光是用一帧或几帧结构进行分析是明显无法全面描述动态环境中出现的相互作用的情况的。amIGM方法将mIGM扩展到了动态环境中的弱相互作用的分析上，对分子动力学模拟得到的几百甚至几千帧结构进行计算，其等值面能够展现模拟过程中存在的分子间的平均的相互作用，具体实现细节可以看前述Struct. Bond.原文中的第4节。 amIGM与mIGM间的关系高度类似于aNCI（也叫aRDG）和NCI间的关系，amIGM相较于aNCI有以下关键性优点： (1)amIGM可以自定义片段、只展现特定片段间的相互作用。而aNCI会同时把所有弱相互作用都展示出来，经常导致图像非常混乱。虽然可以靠Multiwfn的格点数据屏蔽功能做后处理来一定程度上屏蔽掉不想要的aNCI等值面，但不仅麻烦而且原理明显不如amIGM优雅 (2)可以通过调节amIGM的等值面数值决定是只展现较显著的弱相互作用，还是也把不太显著的弱相互作用一同展现。而aNCI没法做这个区分 (3)amIGM的等值面明显不像aNCI那么容易在边缘出现很严重的难看的锯齿 (4)aNCI经常会出现一些没任何意义、不对应实际弱相互作用的垃圾等值面，而amIGM没这个问题 相比之下，aNCI没有继续用的价值。Multiwfn里还支持我之前提出的aIGM，是把IGM扩展到研究动态环境，但其等值面和IGM一样非常臃肿，远不如amIGM好，因此也没有使用价值。原理上也可以把依赖于波函数的IGMH方法扩展成aIGMH，但用起来必定十分昂贵（无论是AIMD产生巨量数目的波函数文件还是aIGMH分析过程），所以笔者并没有对其效果进行检验以及在Multiwfn中实现。 1.2 amIGM的实际效果 下面列举amIGM的一些应用例子，令读者能充分了解amIGM的价值，图都来自于amIGM原文。 下面的amIGM和aIGM图是对标况的水盒子进行动力学模拟后得到的平均的sign(λ2)ρ着色的平均的δg_inter等值面图，用于展示一个水分子和周围的水分子间的相互作用，周围的不断运动的水的结构没有画出来。这里sign(λ2)ρ使用标准的色彩变化方式和范围，和《使用mIGM方法基于几何结构快速图形化展现弱相互作用》（//www.umsyar.com/755）里第一张图相同，后同。由下图可见amIGM图很好地把这个水分子作为氢键给体和作为氢键受体的特征展现了出来，等值面精确出现在氢键会出现的范围，非常对称，而且形状优雅、易于观看，并且淡蓝色还明确体现出其作用强度比普通范德华作用更大。aIGM的等值面虽然和amIGM相似，但过于鼓囊、肥大，丑多了，而且O-H冲着的等值面还蓝得过头了（显得化学键作用似的）。   上面的mIGM图是对动力学过程中随便取的一帧算的，可见中心水分子周围四个等值面虽然各对应一个氢键，但与水分子的C2v对称性明显不符，显然无法像amIGM那样如实地描述液态水中的平均的相互作用。 作为对比，下面给出aNCI方法的图，做法在//www.umsyar.com/186里有。下图左边部分是aNCI直接出的图，可见有一大堆红蓝相间的零碎的等值面，严重扰乱视觉。用Multiwfn做格点数据屏蔽后得到右图，依然明显不理想。虚线标注的那一大坨红蓝相间的等值面完全意义不明，箭头所示的绿色分叉的等值面的存在也难以理解，说不清楚是什么玩意儿。可见aNCI远不如amIGM。   下面这张图是amIGM方法展示的水中的苯酚与周围的水之间的相互作用，三张图分别对应不同数值的δg_inter等值面。可见0.005的图把苯酚与水之间的氢键作用区域很清楚地展现了出来，数值更小的0.003的图还同时把水与苯酚之间形成pi-氢键的位点展现了出来，数值最小的0.0018的图还同时把苯酚与水之间范德华作用为主的区域明确地展现了出来。   范德华势是笔者之前提出的弱相互作用的重要分析方法，在《谈谈范德华势以及在Multiwfn中的计算、分析和绘制》（//www.umsyar.com/551）中对螺烯吸附He原子体系用范德华势等值面图的分析方法很好地解释了动力学模拟得到的He的密度分布。下图(b)是模拟过程中He的位置的叠加图（每隔一定帧数绘制一次，按照帧号由小到大颜色按照蓝-白-红变化）。下图(a)是amIGM等值面图，由于螺烯与He的范德华作用极弱，所以平均的δg_inter的等值面数值用的是非常小的值。可见amIGM图很清晰直观地展现了在模拟过程中螺烯在哪些区域与He原子有相对明显的相互作用，和(b)图能很好地对应上，充分说明了amIGM的合理性。下图(c)图是aNCI图，可见效果十分不堪入目，乱七八糟！   下面是富勒烯在碳纳米管体系中的动力学模拟，(a)是体系示意图，(b)是模拟过程中富勒烯质心位置，可见由于动能-势能反复交换，富勒烯在碳管中反复穿梭做活塞运动。(c)是对这个动力学过程绘制的amIGM图，两种视角都给出了。可见等值面均匀、完整覆盖了碳管内壁，充分体现出整个模拟过程中富勒烯与碳管在这些区域的作用相当充分和均匀。   我在《使用Multiwfn做aNCI分析图形化考察动态过程中的蛋白-配体间的相互作用》（//www.umsyar.com/591）中曾示例如何用aNCI方法考察苄脒阳离子配体与胰蛋白酶间的动态平均的弱相互作用，下图是这个体系的amIGM分析的结果。(a)是0.006 a.u.等值面，(b)是更小的0.003 a.u.等值面，注意蛋白质部分的结构是动力学轨迹中随便取的一帧绘制的。由(a)图可清晰地看到配体的两个氨基与蛋白质有四处鲜明的氢键作用（作用中心区域的等值面明显发蓝）。(b)图的等值面范围更大，可以进一步看到配体在很多区域和蛋白质间也有显著的范德华主导的相互作用，对应绿色等值面区域。   2 amIGM分析的输入文件和一些要点 这一节专门说一下用于Multiwfn做amIGM分析用的输入文件如何合理地准备，以及计算必须知道的要点。 做amIGM需要用户提供跑分子动力学得到的xyz轨迹文件，对此格式不了解的话看《谈谈记录化学体系结构的xyz文件》（//www.umsyar.com/477）。VMD可以载入GROMACS、AMBER、CP2K、LAMMPS等诸多程序跑动力学得到的轨迹文件，选择File - Save coordinate并且选xyz格式，就可以得到xyz轨迹文件。 标准的xyz文件里每个原子的名字都是元素名，给Multiwfn做amIGM分析的文件应当满足这一点，因为只有判断对了元素，amIGM分析的结果才是合理的。如果原子的名字不是元素名，Multiwfn会根据名字猜元素（自动去掉其中的数字，然后试图匹配各个元素名以确定元素，比如N2会被判断为氮，CA会被判断为钙），若猜错了则准分子密度无法正确构造，会导致amIGM结果虚假甚至离谱。自己用文本编辑器打开xyz文件看一下便知是否记录的是元素名，也可以看Multiwfn载入xyz文件后看屏幕上提示的化学组成是什么（即各种元素都是多少个原子）确认Multiwfn是否判断对了元素。 注：pdb格式专门定义了一列用于记录元素名，当前目录下如果有和载入的xyz文件名相同但后缀是pdb的文件，Multiwfn就会优先从pdb文件记录元素的那一列里读取元素名。如果pdb里对某些原子没有提供元素名，对这些原子Multiwfn仍会根据xyz文件里的原子名去猜。 做amIGM分析和mIGM一样需要定义片段。一般只定义两个片段，这种情况下，第1个片段必须对应轨迹中坐标被固定的部分，amIGM图会展示它与第2个组（运动的部分）的相互作用。例如之前展示的例子中，纯水盒子的模拟中一个水要被固定并作为第1个组；苯酚+水的轨迹中苯酚要被固定并作为第1个组；配体+蛋白质的轨迹中配体是被固定的并作为第1个组；富勒烯+碳纳米管的轨迹中碳纳米管要被固定并作为第1个组...再比如，如果研究固体表面吸附分子，表面要被固定并作为第1个组。被固定的部分在动力学模拟过程中可以用冻结（freeze）将坐标严格固定住，也可以用位置限制势结合足够大的力常数将坐标基本固定住。也可以动力学时先不做固定，等动力学跑完之后用VMD的Extensions - RMSD trajectory tool插件或者用GROMACS的trjconv命令结合-fit rot+trans做叠合（align）将第1个组对应的部分消除平动转动（但如果在align处理后第1个组的结构在整个轨迹中还是有较大变化，比如某个要考察相互作用的基团反复翻转，就不适合做amIGM分析了，至少对这个区域而言）。 模拟的体系里面的原子数不要太多，够用就行了，要不然amIGM分析起来慢，轨迹文件也大。比如考察水中的苯酚与水的相互作用，只要模拟的盒子比苯酚稍微大一圈，往盒子里充满水，进行动力学模拟就行了，不要把盒子尺寸弄得过大而导致填充的水太多。做amIGM分析用的动力学轨迹也可以是从完整轨迹中抠出来的一部分。比如考察水中的蛋白质与配体间的相互作用，做动力学的时候肯定是蛋白质+配体+溶液的模型，而做amIGM分析用的轨迹文件只需要包含蛋白质+配体部分即可，水的部分应该去除以避免严重浪费时间。当蛋白质较大的时候，为了避免耗时过高，还应当借助VMD的选择语句（见《VMD里原子选择语句的语法和例子》//www.umsyar.com/504）只把配体以及与配体有明显相互作用的氨基酸残基抠出来成为一个簇并保存成轨迹文件。 用于amIGM分析的凝聚相体系（如水中的苯酚）的动力学模拟最好在NVT下进行（以避免控压导致的原子位置改变对结果可能产生的不良影响），在此之前应先做NPT平衡相模拟使得盒子尺寸达到充分的平衡。 amIGM分析用的轨迹文件记录的帧显然应该涵盖你希望考察弱相互作用的模拟阶段。比如你希望考察稳定状态下蛋白质与配体的相互作用，显然应该取结构已经达到稳定的阶段。而如果模拟过程中配体反复在两种明显不同的构象A-B之间变化，我建议对轨迹做簇分析，把处于A构象的那些帧和B构象的那些帧分别提取成两个轨迹文件，分别做amIGM分析，从而能考察这两个构象与蛋白质相互作用的区别。 amIGM分析耗时和考虑的轨迹帧数呈线性正比。显然考虑的帧数越多，amIGM的结果原理上越好、其图像越能如实表现出模拟过程中的片段间平均的弱相互作用。一般建议至少考虑200帧。如果发现效果不够好，比如等值面不够连贯、不能如预期般充分完整展现平均相互作用，则应当考虑更多帧。这一点类似于计算空间分布函数（sdf），采样越充分越好。 再说一下amIGM分析用的格点。读者请先阅读《Multiwfn FAQ》（//www.umsyar.com/452）的Q39了解格点设置的基本知识。格点数据所处的盒子范围固定的情况下，格点间距越小，amIGM等值面越平滑、格点数越多、耗时越高、Multiwfn导出的cube文件越大，通常0.2 Bohr格点间距就已经够用了，要更好一点可以用0.15 Bohr。盒子的范围应该囊括并尽量刚好囊括感兴趣的弱相互作用可能出现的区域，避免该有的等值面不出现或被截断，以及避免浪费格点用于描述无意义的区域。 Multiwfn的amIGM分析过程中会自动使用节约耗时的策略，如果一个格点与第1个片段在第1帧中的任意一个原子的距离在其scaled范德华半径内，而且这个格点与其它片段的任意原子的距离在整个轨迹中至少有一次在其scaled范德华半径内，这个格点才会被考虑，否则在计算过程中会被忽视。这里scaled范德华半径是指范德华半径与Multiwfn的settings.ini中的amIGMvdwscl参数的乘积。amIGMvdwscl默认为2.0，足够稳妥又能显著节约耗时。但如果极个别情况下发现等值面看起来像是被截断了，应将amIGMvdwscl设得更大或者设为0关闭此策略。 3 amIGM计算实例 本节将会具体演示如何绘制出来1.2节展示的各种amIGM图，用到的xyz轨迹文件都给出了，除了螺烯+He原子的体系外都是我用GROMACS程序跑出来的。GROMACS做动力学模拟的方法在北京科音分子动力学与GROMACS培训班（http://www.keinsci.com/KGMX）里讲授得极为详细，学过一遍后这些轨迹文件肯定都能自己跑出来。螺烯+He原子的体系是我用xtb程序跑的，在北京科音高级 培训班（http://www.keinsci.com/KAQC）的“从头算分子动力学”部分专门全面系统讲了怎么用xtb跑动力学并给了很多例子，其中就包括这个例子。 值得一提的是，虽然amIGM分析可以支持周期性，但这些xyz文件里并没有加入Multiwfn可以识别的盒子信息（如果加入的话，需要在注释行写入盒子的三个矢量，诸如Tv_1: 7.426 0.0 0.0 Tv_2: 3.6 6 6.40 0.0 Tv_3: 0.0 0.0 10.0）。因为这些例子中计算格点数据的区域都没有离模拟用的盒子边界太近，因此不需要考虑周期性。 3.1 水中的苯酚与水的相互作用 这个例子的轨迹文件phenol_wat.xyz在此压缩包里：//www.umsyar.com/attach/759/phenol_wat.7z。体系共1342个原子，在模拟过程中苯酚在盒子中央被固定。NVT下模拟了1 ns（之前已经过了平衡相模拟），每1 ps保存一帧，故xyz轨迹文件共1001帧。模拟控温在298.15 K。第1帧的结构如下所示，苯酚原子（1-13号原子）用CPK方式显示，水用透明棍状方式显示。   启动Multiwfn，载入phenol_wat.xyz，然后输入 20  //弱相互作用可视化分析 -12  //amIGM分析 2  //定义两个片段 1-13  //作为第1个片段的苯酚的原子序号 c  //其它原子作为第2个片段 [回车]  //考虑所有帧 11  //围绕接下来指定的片段往各方向延展一定距离来定义计算格点数据的空间范围 1-13  //苯酚片段的原子序号 3 A  //延展距离为3埃，对此例够大了 [回车]  //用默认的0.2 Bohr格点间距。这个格点间距对应的图像质量足够好了，想要等值面的边缘锯齿更少可以降低到0.15 Bohr，但耗时多一倍多 在双路7R32服务器上96核并行计算的总共耗时约100秒。在后处理菜单选3将平均的sign(λ2)ρ和平均的δg_inter格点数据分别导出为当前目录下的avgsl2r.cub和avgdg_inter.cub，然后放到VMD目录下。再把Multiwfn的examples目录下的aIGM.vmd脚本放到VMD目录下。之后启动VMD，在命令行窗口输入source aIGM.vmd执行脚本，就会马上显示出amIGM图。在VMD里进入Graphics - Representation界面可以看到两个显示方式（简称为Rep）。选择Isosurface对应的Rep后可以在界面里改Isovalue值，在Trajectory标签页里可以改色彩刻度范围。将Isovalue设为比如0.003，并且将CPK显示方式的Rep的Selected Atoms设为serial 1 to 13后，看到的图像就和1.2节展示的这个体系的图像一样了，一些箭头、文字自行ps上去即可。 如aIGM.vmd脚本的内容可以看到，默认的等值面数值是0.008，默认的色彩刻度范围是-0.05到0.05，都可以自己改。色彩刻度范围设得越小，等值面上不同特征的相互作用区域的颜色差异越大。 3.2 螺烯与被吸附的He原子的相互作用 这个例子的通过xtb程序用GFN0-xTB方法跑出来的轨迹文件helicene-He.xyz在此压缩包里：//www.umsyar.com/attach/759/helicene-He.7z。体系共43个原子，控温在10 K做了2500 ps动力学，从中均匀取了1000帧。轨迹跑完后，用VMD自带的RMSD Trajectory tool工具对螺烯部分做了align来消除其平动和转动，因此虽然在模拟过程中没有对螺烯做固定，经过这么处理后，整个轨迹中螺烯的位置、朝向都没变。并且本身螺烯是刚性的，结构波动程度甚微。 启动Multiwfn，载入helicene-He.xyz，然后输入 20  //弱相互作用可视化分析 -12  //amIGM分析 2  //定义两个片段 1-42  //作为第1个片段的螺烯的原子序号 c  //其它原子作为第2个片段 [回车]  //考虑所有帧 11  //围绕接下来指定的片段往各方向延展一定距离来定义计算格点数据的空间范围 1-42  //螺烯片段的原子序号 2 A  //延展距离设2埃 [回车]  //用默认的0.2 Bohr格点间距 在双路7R32服务器上96核并行计算的总共耗时约40秒。按前例在VMD中显示出等值面后，再把等值面数值设为很小的0.00004，图像就和第1.2节所示的这个体系的图一样。 3.3 富勒烯与碳纳米管的相互作用 这个例子的GROMACS程序结合GROMOS力场跑出来的轨迹文件CNT_C60.xyz在此压缩包里：//www.umsyar.com/attach/759/CNT_C60.7z。体系共456个原子，一开始把富勒烯放在距离碳纳米管管口一定距离处。模拟过程控温在200 K。模拟共100 ps，每0.5 ps保存一帧，故轨迹共201帧。轨迹跑完后，用VMD的RMSD Trajectory tool工具对碳纳米管部分做了align来消除其整体运动。 启动Multiwfn，载入CNT_C60.xyz，然后输入 20  //弱相互作用可视化分析 -12  //amIGM分析 2  //定义两个片段 1-396  //作为第1个片段的碳纳米管的原子序号 c  //其它原子作为第2个片段 [回车]  //考虑所有帧 11  //围绕接下来指定的片段往各方向延展一定距离来定义计算格点数据的空间范围 1-396  //碳纳米管片段的原子序号 0  //延展距离设为0，因为等值面只在纳米管内壁区域出现 0.15  //用0.15 Bohr格点间距 按照之前的做法在VMD里作图，等值面数值设0.0001，就可以看到1.2节给出的此体系的图像。 3.4 配体与蛋白质的相互作用 这个例子的动力学模拟是用GROMACS程序结合GROMOS力场对完整蛋白质+配体+溶剂做的，在NPT预平衡后在NVT系综下模拟了1 ns，每1 ps保存一帧，共1001帧，模拟过程中将配体部分做了冻结。随后用VMD将轨迹中配体与相邻的蛋白质残基截出来构成团簇，共161原子，轨迹文件cluster.xyz在此压缩包里：//www.umsyar.com/attach/759/cluster.7z。由于这个xyz文件里的原子名不是元素名，因此此压缩包里还提供了cluster.pdb，里面记录了各个原子的元素信息，从而使得Multiwfn能正确识别元素。这些文件的产生方式在《使用Multiwfn做aNCI分析图形化考察动态过程中的蛋白-配体间的相互作用》（//www.umsyar.com/591）的第2、3节有详细说明（那篇文章里文件名叫cluster），因此这里不再累述。 把cluster.xyz和cluster.pdb放在相同目录下。启动Multiwfn，载入cluster.xyz，如屏幕上的提示所示元素信息从cluster.pdb中读取了，显示的化学组成H72 C53 N16 O18 S2完全正确。然后依次输入 20  //弱相互作用可视化分析 -12  //amIGM分析 2  //定义两个片段 144-161  //作为第1个片段的配体的原子序号 c  //其它原子作为第2个片段 [回车]  //考虑所有帧 11  //围绕接下来指定的片段往各方向延展一定距离来定义计算格点数据的空间范围 144-161  //配体的原子序号 3 A  //延展距离设为3埃 0.15  //用0.15 Bohr格点间距 在双路7R32服务器上96核并行计算的总共耗时110秒。按照之前的做法在VMD里作图，等值面数值设0.006或0.003，并恰当在VMD里设置配体和周围残基的显示方式，就可以看到1.2节给出的此体系的图像。 4 总结 本文充分介绍了笔者在Struct. Bond., 190, 297 (2026)中提出的amIGM分析方法，并通过四个例子演示了amIGM在Multiwfn中的分析操作。由本文可见，amIGM分析可以非常生动直观地展现动力学过程中的分子间的平均的弱相互作用，而且分析步骤操作简单，因此非常推荐大家将之应用于实际问题的研究中。amIGM的图像效果远好于aNCI，而耗时相仿佛，因此aNCI就不需要再考虑了。 //www.umsyar.com/758 2025-12-07T14:21:00+08:00 Multiwfn结合ORCA的TDDFT计算做空穴-电子等分析的方法 Method of performing hole-electron and relevant analyses via Multiwfn in combination with TDDFT calculation of ORCA 文/Sobereva@北京科音  2025-Dec-7 0 前言 Multiwfn是电子激发分析的十分强大的武器库，在《Multiwfn支持的电子激发分析方法一览》（//www.umsyar.com/437）里有全面盘点。其主功能18中有很多分析方法依赖于参考态轨道波函数和激发态组态系数，典型的如《使用Multiwfn做空穴-电子分析全面考察电子激发特征》（//www.umsyar.com/434）介绍的空穴-电子分析和δr指数、《在Multiwfn中通过IFCT方法计算电子激发过程中任意片段间的电子转移量》（//www.umsyar.com/433）介绍的IFCT分析、《使用Multiwfn做自然跃迁轨道(NTO)分析》（//www.umsyar.com/377）介绍的NTO分析，等等一大堆。它们通常结合目前研究电子激发最常用的TDDFT方法使用。在这些文章里我都明确示例了怎么结合Gaussian的TDDFT计算实现这些分析，而如今免费又强大的 程序ORCA的用户也很多，因此在本文专门示例一下怎么结合ORCA来实现。本文只使用如今特别流行的空穴-电子分析作为例子，其它的需要轨道+组态系数的分析方法所用的文件与空穴-电子分析完全一样，就不一一举例了。 这里我假定读者已经看过上述博文了解了空穴-电子分析的原理和用法。如果读者不了解Multiwfn的话，看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。读者请务必使用2025-Dec-7及以后更新的Multiwfn版本，否则情况与本文说的不符。Multiwfn可以在官网//www.umsyar.com/multiwfn免费下载。读者务必使用ORCA 6.1.1及以后的版本，本文用的是ORCA 6.1.1。我假定读者已了解ORCA的各种常识和做TDDFT的方法，在《北京科音高级 培训班（http://www.keinsci.com/KAQC）中对ORCA的各方面使用有极为全面系统的讲解，ORCA的安装方法见《 程序ORCA的安装方法》（//www.umsyar.com/451）。 本文的例子使用与《使用Multiwfn做空穴-电子分析全面考察电子激发特征》中相同的D-pi-A体系，只不过基组从6-31G*改为了ORCA用户更常用的def2-SVP。所有例子的输入输出文件都在//www.umsyar.com/attach/758/file.rar里。 1 TDA-DFT情况下的空穴-电子分析 在演示使用严格的TDDFT做电子激发分析之前，这一节先说ORCA默认用的TDA近似下的TDDFT的情况，即TDA-DFT的情况，这种情况最为简单。 使用ORCA运行以下输入文件，对应本文文件包里的ORCA_TDA_D-pi-A\TDA.inp，算完后得到输出文件TDA.out。还同时在当前目录下产生了一堆其它文件，除了TDA.gbw外都删掉。 ! CAM-B3LYP def2-SVP pal16 tightSCF miniprint %tddft   nroots 5   tprint 1E-8 end %maxcore 3000 * xyz 0 1  C                  3.55863800   -1.13874700    0.39983600  C                  2.17031600   -1.13391800    0.39359500 ...略 * 注意%tddft中的tprint 1E-8很重要，代表对电子激发贡献大于1E-8的组态都输出系数，相当于Gaussian的IOp(9/40=4)。想让结果更精确也可以设得更小，但改进甚微，而输出文件尺寸会增大、Multiwfn的分析耗时会增加。 运行orca_2mkl TDA -molden命令，会基于TDA.gbw转换出TDA.molden.input，其中记录了参考态轨道波函数，类似于Gaussian的fch文件。 启动Multiwfn，载入TDA.molden.input，然后输入 18  //电子激发分析 1  //空穴-电子分析 [回车]  //这一步让用户输入ORCA的输出文件路径，直接按回车代表载入与输入文件同路径的同名但后缀为out的文件（TDA.out） 2  //以分析第2激发态为例 现在Multiwfn就从TDA.out中载入了组态系数。之后可照常做后面的分析，比如输入 1  //可视化空穴、电子、跃迁密度等函数 2  //中等质量格点 3  //同时显示空穴和电子的等值面 2 TDDFT情况下的空穴-电子分析 至少截止到笔者撰文时（2025-Dec-7）的最新版ORCA 6.1.1来说，ORCA尚无法在TDDFT计算时将激发组态系数和去激发组态系数分别输出。为了使得Multiwfn基于ORCA的TDDFT计算做空穴-电子分析可行，需要导出记录所有组态系数的json文件并令Multiwfn从其中读取。 ！！注意！！由于笔者撰文时最新版ORCA对于参考态为开壳层的情况无法在json文件中以正确格式输出组态系数，因此本节的做法目前只适用于参考态为闭壳层的情况！等以后ORCA修正了这个bug，Multiwfn也会使本节的做法支持开壳层参考态。 运行以下输入文件，对应本文文件包里的ORCA_TDDFT_D-pi-A\TDDFT.inp，算完后得到输出文件TDDFT.out。还同时在当前目录下产生了一堆其它文件，除了TDDFT.gbw外都删掉。 ! CAM-B3LYP def2-SVP pal16 tightSCF miniprint %tddft   nroots 5   tda false end %maxcore 3000 * xyz 0 1  C                  3.55863800   -1.13874700    0.39983600  C                  2.17031600   -1.13391800    0.39359500 ...略 * 运行orca_2mkl TDDFT -molden命令，会基于TDDFT.gbw转换出TDDFT.molden.input。 创建一个文本文件叫TDDFT.json.conf，在里面写入如下内容 { "CIS": true, "CISNRoots": true } 然后运行orca_2json TDDFT.gbw，就会在当前目录下产生TDDFT.json。 说明：TDDFT.json.conf是orca_2json的控制文件，控制orca_2json在处理TDDFT.gbw时用的设置，当前内容要求把所有激发的组态系数都写入json文件，而默认不会写入。若控制文件的名字为orca.json.conf，则对于当前目录下任意名字的gbw在用orca_2json处理时都会用其中的设置。 启动Multiwfn，载入TDDFT.molden.input，然后输入 18  //电子激发分析 1  //空穴-电子分析 [回车]  //载入与输入文件在同目录下同名的out文件（TDDFT.out） 2  //以分析第2激发态为例 此时如屏幕上的提示所示，由于Multiwfn发现在TDDFT.out的同目录下有名字相同但后缀为json的文件（TDDFT.json），Multiwfn就自动从json文件中载入了组态系数。json文件里记录的是全部组态系数，Multiwfn默认只载入绝对值大于1E-4的组态系数以节约之后的分析时间，这个阈值可以通过settings.ini中的ORCAloadcoeff修改。 现在就进入了空穴-电子分析界面，照常做分析即可，结果和同级别下Multiwfn+Gaussian的分析结果差异可忽略不计。 技巧：如果你把Multiwfn的settings.ini里的orca_2mklpath设为了orca_2mkl的实际路径，则以上例子中不需要手动用orca_2mkl转换gbw文件成molden.input文件，而是可以直接用Multiwfn载入gbw文件，此时Multiwfn会自动将之转换为.molden.input文件并载入，然后再删掉molden.input。 顺带再举个例子，对ORCA的TDDFT计算使用《使用Multiwfn便利地查看所有激发态中的主要轨道跃迁贡献》（//www.umsyar.com/529）介绍的功能。启动Multiwfn后载入TDDFT.out，然后进入主功能18的子功能15，就得到了如下结果，所有5个态的组态系数也都是从json文件中读取的。  #   1   3.8520 eV    321.87 nm   f=  0.01007   Spin multiplicity= 1:    H-4 -> L 82.3%, H-4 -> L+2 13.0%  #   2   4.0470 eV    306.36 nm   f=  0.64110   Spin multiplicity= 1:    H -> L 87.4%, H-3 -> L 5.2%  #   3   4.3730 eV    283.52 nm   f=  0.00009   Spin multiplicity= 1:    H-6 -> L 84.4%, H-6 -> L+2 12.7%  #   4   4.7440 eV    261.35 nm   f=  0.01900   Spin multiplicity= 1:    H -> L+1 38.7%, H-2 -> L 25.3%, H -> L+3 18.0%, H-1 -> L 5.8%  #   5   4.8350 eV    256.43 nm   f=  0.00409   Spin multiplicity= 1:    H -> L+3 43.3%, H-2 -> L 41.1%, H-3 -> L+1 5.1% //www.umsyar.com/756 2025-11-24T11:30:00+08:00 通过格点屏蔽巨幅降低IGMH可视化分析片段间相互作用的耗时 Significantly reducing computational cost of IGMH visual analysis for interfragment interactions via grid screening 文/Sobereva@北京科音  2025-Nov-24 1 说明 《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）中介绍的笔者提出的IGMH方法已经被广泛用于直观分析考察化学体系中的相互作用。这种分析最主要的组成部分是计算δg、δg_inter、δg_intra、sign(λ2)ρ的格点数据，这也是这种分析最耗时的过程。对于大体系，特别是在CPU比较弱的机子上，IGMH计算格点数据的分析耗时往往挺高。 大多数人用IGMH方法主要是作图分析考察片段间相互作用，这种情况实际上只需要对片段间重叠的区域计算δg_inter和sign(λ2)ρ格点数据就行了，其它区域的格点完全都不需要计算。考虑到这一点，从2025-Nov-24更新的Multiwfn版本开始，支持了一种巨幅节约IGMH图形化考察片段间相互作用耗时的做法，具体来说是这样的：在Multiwfn的settings.ini文件中现在有一个参数名为IGMvdwscl，一个片段的表面是由其中各个原子的范德华半径乘以这个系数对应的原子球叠加构成的表面。在做IGMH或mIGM或IGM分析过程中计算格点数据的那一步，如果某个格点出现在两个或多个片段的表面之间的重叠区域内，则这个格点就会被计算，否则会被直接忽略。利用这个策略，往往可以令IGMH格点数据计算过程耗时降低好几倍甚至更多。能节约百分之多少耗时具体取决于你设的格点数据分布的盒子范围包含了多少可以被忽略的格点。显然，用了这个策略后就只能对δg_inter函数的等值面作图了，而不能对δg和δg_intra等值面作图，因为只有δg_inter的等值面才肯定会落在片段间重叠区域。 IGMvdwscl默认为0，代表不启用这种节约耗时的策略。如果要启用此策略，就恰当设其数值，设得越大降低耗时效果越低，设得越小越有可能截断δg_inter等值面。一般建议用2.0，既足够安全，节约耗时的效果又足够显著。 2 实例 在《8字形双环分子对18碳环的独特吸附行为的 、波函数分析与分子动力学研究》（//www.umsyar.com/674）文中介绍了我研究的OPP双环分子与两个18碳环形成的超分子复合物，文中给出了IGMH图，OPP、第一个碳环、第二个碳环各定义为一个片段。以这个体系为例，我们对比一下用和不用前述加速策略的耗时，在拥有24个物理核心的i9-13980HX CPU上Multiwfn用24核并行计算。这个体系的波函数文件2C18_OPP.wfn在这个压缩包里：//www.umsyar.com/attach/756/2C18_OPP.7z。体系总共260个原子，5504个GTF（碳环用的6-311G*基组，OPP双环分子用的6-31G*）。 启动Multiwfn，载入2C18_OPP.wfn，然后输入 20  //弱相互作用可视化分析 11  //IGMH分析 3  //定义三个片段 243-260   //第1个片段原子序号 225-242  //第2个片段原子序号 c  //其它原子作为第3个片段 4  //自定义格点间距或格点数，格点覆盖整个体系 0.2  //0.2 Bohr 格点数据计算总耗时达到600秒，虽然不算太长，但也得在屏幕前等好一阵，对于明显更大的体系或者更弱的CPU没准要花一个小时。对导出的格点数据用VMD作图，0.002等值面的图如下所示。 注意在Multiwfn计算完格点数据后在屏幕上还显示了三种函数的全空间积分值，对应相应格点数据数值的总和乘上格子体积：  Integral of delta-g over whole space:        231.384179 a.u.  Integral of delta-g_inter over whole space:    1.184764 a.u.  Integral of delta-g_intra over whole space:  230.199414 a.u. 下面看看使用节约耗时策略后的情况。把settings.ini里的IGMvdwscl设为2，重新启动Multiwfn并重复前面的操作，计算格点数据之前会有以下提示，说明90.96%的格点都被忽略掉了。  Prescreening grids with IGMvdwscl parameter:  2.00  Percent of screened grids:     90.96% 计算总共耗时仅为90秒，只有之前耗时的15%！当前看到的格点数据积分值如下，可见δg_inter函数的积分值几乎没变，充分说明了当前没有忽略掉δg_inter主要分布区域的格点。  Integral of delta-g over whole space:         41.955223 a.u.  Integral of delta-g_inter over whole space:    1.179150 a.u.  Integral of delta-g_intra over whole space:   40.776072 a.u. 用导出的格点数据再次绘制sign(λ2)ρ着色的δg_inter图，会看到得到的图和上面的没有任何差别！另外，使用这种降低耗时的策略也完全不会影响δg_inter与sign(λ2)ρ之间的散点图。 由此例可见，将IGMvdwscl设为2是非常安全的巨幅加速IGMH可视化展现片段间相互作用的策略，推荐使用！这个策略对《使用mIGM方法基于几何结构快速图形化展现弱相互作用》（//www.umsyar.com/755）介绍的mIGM方法同样奏效，只不过mIGM方法本身耗时就远远低于IGMH，对于本文的例子一瞬间就算完，因此这个策略对mIGM不会带来可查觉的收益。 至于IGMH、mIGM、IGM分析的后处理菜单中的计算δG_atom和δG_pair指数的功能，和本文介绍的加速策略无关，IGMvdwscl的设置不影响其耗时和结果。 //www.umsyar.com/755 2025-11-24T11:26:00+08:00 使用mIGM方法基于几何结构快速图形化展现弱相互作用 Using mIGM method to rapidly graphically represent weak interactions based on geometry 文/Sobereva@北京科音   2025-Nov-24 0 前言 弱相互作用的图形化分析方法日趋流行，其中笔者提出的IRI和IGMH方法如今已经被使用得相当普遍，另外还有NCI、aNCI、IGM等方法。笔者专门有两篇综述文章对其进行了非常全面、完整的介绍，如果还没读过的话强烈建议阅读： Angew. Chem.上发表了全面介绍各种共价和非共价相互作用可视化分析方法的综述 //www.umsyar.com/746 一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！ //www.umsyar.com/667 近期笔者提出了一种新的弱相互作用可视化方法称为mIGM（modified independent gradient model）并已实现在了Multiwfn程序中。如下一节所示，mIGM有重要、独特的价值。mIGM的原文如下，非常欢迎阅读和引用。使用Multiwfn做mIGM分析时请将此文与Multiwfn程序启动时提示的程序原文一起引用。 Tian Lu, Graphically revealing weak interactions in dynamic environments using amIGM method, Struct. Bond., 190, 297 (2026) DOI: 10.1007/430_2025_95 注1：此文同时作为Computational Methods for the Analysis of Non-Covalent Interactions书的一章出版 注2：mIGM方法是我提出amIGM方法时的重要副产物，所以文章标题上写的是amIGM 注3：如果没有权限访问，也可以看ChemRxiv上的预印版https://doi.org/10.26434/chemrxiv-2025-zts13，内容和正式版没明显区别。引用时请引用正式版 1 mIGM方法介绍 IGMH在前述综述里以及《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）中都有充分的介绍，这里就不再多说了。mIGM方法是IGMH的变体，把IGMH中基于波函数计算的实际电子密度替换为准分子密度（promolecular density），准分子密度只需要元素和原子坐标信息就可以由Multiwfn构造出来。 mIGM分析的优点是仅依赖于几何结构而不需要波函数信息，且速度远快于IGMH。IGM方法虽然也有这方面的优点，但mIGM的图像效果明显好于IGM。因此IGM方法可以完全被淘汰了。 mIGM显然不如基于波函数计算的IGMH严格，但对于图形化考察弱相互作用的目的，多数情况下其结果不比IGMH差多少（但对于片段带显著净电荷的情况有可能差得多一点）。当然，对于不大的体系、不难获得波函数文件的情况，还是建议优先用IGMH。 由于以上特点，mIGM很适合这三种情况： (1)需要快速得到结果的场合，尤其适用于很大体系和大批量分析 (2)作为IGMH分析前的快速预览目的 (3)不便于得到Multiwfn支持的波函数文件的情况。例如Multiwfn不支持QE、VASP、M$等完全基于平面波的第一性原理程序产生的波函数文件，这些程序的用户可以把优化完的几何结构保存成cif格式作为Multiwfn做mIGM分析的输入文件）。甚至一点都不会理论计算的人也可以直接拿较准确实验得到的pdb、cif文件做mIGM分析。 这里以苯基丙氨酸三聚体为例对比一下不同弱相互作用可视化方法的结果，先来看IGMH、IGM、NCI方法的情况。下面这张图是mIGM原文里的图，给出了sign(λ2)ρ着色的不同方法定义的等值面，IGMH和IGM分析中把每个分子被定义为了一个片段。可见IGMH图像很理想，等值面光滑、形状优雅，很好地展现出了普通范德华作用为主和氢键作用为主的相互作用区域。IGM的图很糟，等值面过于肥大，往往离原子核太近，并进而导致在箭头所示的局部区域出现不合理的橙色着色。NCI方法无法区分分子内和分子间相互作用，而且等值面显得零散、稀碎，在一些地方还有很难看的锯齿（把格点间距降到很小可以避免，但巨幅增加耗时和格点数据尺寸）。 下面再来看mIGM的图。这里给出的是sign(λ2)ρ着色的mIGM的δg_inter函数分别为0.006、0.01、0.03 a.u.的等值面，sign(λ2)ρ的色彩变化方式同前。可见效果很好，与IGMH图的效果十分接近。并且mIGM和IGMH一样，等值面的数值设得越大可以着重展现越强的相互作用。0.03的等值面中只出现了相对较强的分子间氢键作用，数值中等的0.01图中还能看到较显著的范德华作用区域，数值最小的0.006图中还把特别弱的范德华作用区域也展现了出来。 接下来就给出Multiwfn做mIGM分析的两个简单例子。第一个例子重复上面苯基丙氨酸三聚体的图，第二个例子以碲晶体为例演示将mIGM用于周期性体系。实际上mIGM和IGMH、IGM分析过程几乎完全相同，唯一差异仅仅是在Multiwfn主功能20里选择mIGM而非IGMH或IGM。因此如果你之前就看过//www.umsyar.com/621会了IGMH分析，或看过《通过独立梯度模型(IGM)考察分子间弱相互作用》（//www.umsyar.com/407）会了IGM分析，用mIGM其实都没有需要额外学的。 读者务必使用2025年11月23日及以后更新的Multiwfn版本，否则情况和下文所述不符。Multiwfn可以在官网//www.umsyar.com/multiwfn免费下载。不了解Multiwfn者建议看《Multiwfn入门tips》（//www.umsyar.com/167）和《Multiwfn FAQ》（//www.umsyar.com/452）。本文例子用的VMD是1.9.3版，可以在http://www.ks.uiuc.edu/Research/vmd/免费下载。 2 mIGM分析实例1：苯基丙氨酸三聚体 为了做mIGM分析，需要给Multiwfn提供记录了几何优化后的苯基丙氨酸三聚体结构信息的文件。Multiwfn支持什么格式详见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）。常见的xyz、pdb、cif、mol2等大量格式都可以作为输入文件。Multiwfn目录下的examples\phenylalanineresiduestrimer.xyz是优化后的三聚体的结构文件。顺带一提，如果想方便地得到里面各个分子的序号的话，可以将之载入Multiwfn，进入主功能0，选菜单栏的Tools - Select fragment，然后输入某分子中任意原子序号，则这个分子中所有原子序号都会给出，可以复制出来之后粘贴到Multiwfn窗口里用来定义片段。 启动Multiwfn，依次输入 examples\phenylalanineresiduestrimer.xyz 20 //弱相互作用可视化分析 -10 //mIGM分析 3 //定义3个片段 1-29 //第1个片段的原子序号（对应第1个分子） 30-40,52,53,56,57,63-65,77-87 //第2个片段的原子序号（对应第2个分子） c //其它原子作为第3个片段（对应第3个分子） 4 //自定义格点间距。如果你对于定义格点的意义和方式不了解的话，建议参看《Multiwfn FAQ》（//www.umsyar.com/452）中的Q39 0.2 //格点间距设0.2 Bohr，这足以得到足够光滑的mIGM图像，想要效果更好点可以用0.15 一瞬间就算完了，之后选3把格点数据导出成cube文件，当前目录下就产生了dg_inter.cub、dg_intra.cub、dg.cub、sl2r.cub。对于展现片段间相互作用，只需要保留分别对应sign(λ2)ρ和δg_inter格点数据的sl2r.cub和dg_inter.cub文件即可，其它的可以删掉。 将sl2r.cub和dg_inter.cub挪到VMD目录下，并且把Multiwfn自带的examples目录下的IGM_inter.vmd脚本文件挪到VMD目录下。启动VMD，之后在文本窗口输入source IGM_inter.vmd，脚本就会被执行，然后看到下图 当前默认用的δg_inter等值面数值是0.01，如果想改的话，在VMD Main窗口进入Graphics - Representation，在下面标注的文本框里输入数值后按回车即可 如果想把图像效果设得像上一节展示的mIGM原文的图那样，在Graphics - Representation窗口中选中显示体系结构的那个Representation（后文简写为Rep），把Drawing Method从CPK改为Licorice，并把Bond Radius恰当设小。然后选中显示等值面的那个Rep，把Material改成AOShiny。 在Multiwfn做mIGM分析的后处理菜单还可以选择选项6在mIGM框架内计算分别衡量原子和原子对儿对特定两个片段间相互作用贡献程度的δG_atom和δG_pair指数。这两个指数在//www.umsyar.com/621和IGMH原文中都有介绍，在IGMH分析的官方教程//www.umsyar.com/multiwfn/res/IGMH_tutorial.zip中有计算例子。mIGM和IGMH框架下计算的这些指数显然在定量数值上会有所不同，但能说明的问题一致。 这里对前面的图中中央和上方两个分子之间的相互作用计算δG_atom和δG_pair指数。在之前的Multiwfn的mIGM分析的后处理菜单中依次输入 6  //计算δG_atom和δG_pair指数 2,3  //对片段2、3之间计算 3  //用Ultrafine格点计算这些指数 马上就算完了，从屏幕上的提示可见δG_atom和δG_pair指数已经被导出到了当前目录下的atmdg.txt文件中。然后再输入y导出atmdg.pdb，此文件中原子的beta因子对应的是δG_atom指数，可以在VMD中通过对原子着色来直观显示。 把atmdg.pdb载入VMD，背景设成白色，Graphics - Colors - Color Scale里把Method设为BWR（蓝-白-红方式变化），Graphics - Representation里把Coloring Method设为Beta，Trajectory标签页里把Color Scale Data Range下限和上限分别设为-5和5，此时看到的图如下。完全白色对应于δG_atom为0，暗示是对于相互作用没贡献的原子。越红对应δG_atom越大，暗示对分子间相互作用的重要性越大。提醒：δG_atom只是以粗糙方式定义的，切勿视为对片段间相互作用能的确切贡献量，它只适合直观区分不同原子对相互作用可能的重要性、令你能快速判断哪些原子很值得关注。真正要确切得到各个原子对相互作用能的贡献，应当用Multiwfn中的EDA-FF方法，介绍见《使用Multiwfn做基于分子力场的能量分解分析》（//www.umsyar.com/442）。 在atmdg.txt中可以看到片段2、3上原子的δG_atom指数和归一化后的δG_atom%指数，由高到低排序：  Atomic delta-g indices of fragment  2 and percentage contributions  Atom   57 :    0.458233  (  22.76 % )  Atom   87 :    0.215283  (  10.69 % )  Atom   79 :    0.187531  (   9.31 % )  Atom   39 :    0.157183  (   7.81 % ) ...略  Atomic delta-g indices of fragment  3 and percentage contributions  Atom   58 :    0.314871  (  15.64 % )  Atom   69 :    0.297674  (  14.78 % )  Atom   67 :    0.227089  (  11.28 % )  Atom   44 :    0.178482  (   8.86 % ) ...略 还可以看到两个片段的原子之间的由大到小排序的δG_pair指数和归一化后的δG_pair%指数：  Atomic pair delta-g indices and percentage contributions (zero terms are not shown)    57   67 :    0.153236  (   7.61 % )    87   58 :    0.142887  (   7.10 % )    57   55 :    0.077699  (   3.86 % )    53   58 :    0.071819  (   3.57 % ) ...略 可以把上面的一些δG_pair%指数自行通过ps标注在图上，如下所示。 提示：当前atmdg.pdb里只记录了2、3片段的原子，为了确认比如δG_pair指数最大的57、67号原子对应图上哪两个，可以同时把完整结构文件phenylalanineresiduestrimer.xyz载入VMD，选择语句输入serial 57 67并设为VDW方式显示，并把Sphere Scale改小到0.3，这样就能直接在图上看到57和67号原子了。 3 mIGM分析实例2：碲晶体 碲晶体是由一条条碲原子链有序堆积构成的。作为mIGM方法分析周期性体系的例子，这一节使用mIGM直观展现碲晶体中的一条链与周围的链之间的相互作用。 //www.umsyar.com/attach/755/Te.cif是碲晶体的原胞。由于如下所示，其原胞太小，无法在图上展现清楚链之间的相互作用，因此做mIGM分析之前需要先把它扩成足够大的超胞。这会用到《Multiwfn中非常实用的几何操作和坐标变换功能介绍》（//www.umsyar.com/610）里介绍的功能。  启动Multiwfn，载入Te.cif，然后输入 300  //其它功能（Part 3） 7  //几何相关操作 19  //扩胞 4  //第1个轴方向扩成原先的4倍 4  //第2个轴方向扩成原先的4倍 3  //第3个轴方向扩成原先的3倍 此时可以选择0观看一下当前的结构，会看到晶胞大小已经足够大了。我们要把其中最接近晶胞中央的一条链作为第1个片段。为了获得其原子序号，在Multiwfn图形窗口顶端选择Tools - Select fragment，并且输入中央一条链上任意一个原子的序号比如98，就会返回整条链上的原子序号，并且整条链都被高亮显示了，如下所示 点右上角的Return按钮关闭图形窗口，接着在Multiwfn中输入 -10  //从几何操作界面中返回 0  //返回到主菜单 20  //弱相互作用可视化分析 -10  //mIGM分析 2  //定义两个片段 48,51,54,91,92,94,95,97,98  //第1个片段的原子序号 c  //其它原子作为第2个片段 2  //中等质量格点。对于周期性来说，这个选项代表格点间距约0.25 Bohr，格点均匀分布在整个晶胞中。如果想让等值面更平滑，可以选择高质量格点，对应0.15 Bohr格点间距，此时产生的cub文件会是中等质量格点的1/(0.15/0.25)^3=4.6倍 很快就算完了，后处理菜单选2导出cube文件，再把sl2r.cub和dg_inter.cub用IGM_inter.vmd脚本在VMD中作图，默认的0.01等值面如下所示，两种视角都给出了。为了让中间的链和其它的链在图上便于区分，对中间的链设置了两个Rep，选择语句都是serial 48 51 54 91 92 94 95 97 98，着色方式都设为ColorID并选iceblue颜色，其中一个Rep用CPK风格显示，另一个Rep用DynamicBonds风格显示并且把Distance Cutoff改大到3以令碲原子间的键能被显示出来。之后再创建一个DynamicBonds风格的Rep，选择语句写not serial 48 51 54 91 92 94 95 97 98，用Name着色，材质设EdgyGlass。此外，在文本窗口中输入pbc box把格点数据的盒子边框（当前对应晶胞边框）显示出来。 由上图可清晰地看出中间的碲链的每个原子都在好几个方向与周围的链有相对明显的相互作用。还可以把等值面数值改小到0.003以让相互作用展现得更全面，此时图像如下所示，可以看到在有的区域等值面颜色相对发蓝，暗示这些区域的相互作用比绿色区域的更为显著，甚至有可能有一定共价作用成份。这点可以用《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）介绍的Mayer键级在一定程度上论证。如mIGM原文中所提到，用CP2K在PBE/DZVP-MOLOPT-SR-GTH级别下对此体系计算并产生波函数文件，然后用Multiwfn计算原子间Mayer键级，上图中相对近距离接触的两个碲（二者连线穿越上图中偏蓝的等值面）之间的键级为0.137，不是非常接近0，因此暗示了很轻微的共价作用。 4 总结 本文介绍并示例了Struct. Bond., 190, 297 (2026)中提出并在相互作用可视化分析综述Angew. Chem. Int. Ed., 2025, 64, e202504895 (2025)中讲到的mIGM方法。使用Multiwfn做mIGM分析速度极快，只需要原子坐标信息就能算，对于分析弱相互作用来说它通常是非常流行的IGMH方法的很好的近似，因此很有实用性，推荐大家在IGMH算不动或者波函数文件不便于得到的时候使用mIGM。此外，mIGM的提出使得IGM完全失去使用价值了，不用再考虑了。 mIGM分析是基于准分子密度的，这样的近似密度描述分子间相互作用大多能接受，因为分子间相互作用区域电子密度颇低，而且弱相互作用的共价性大多微弱因此对作用区域密度影响较小。而mIGM如果用于描述化学键作用的话就太糙了，因为化学键的形成会导致实际电子密度在成键区域显著偏离准分子密度，此时应当用IGMH（若不需要区分片段内和片段间相互作用则应当用IRI）。 虽然mIGM分析已经很快了，但如果体系巨大或者格点间距设得很小导致格点数很多，计算还是会花一定时间。此时推荐使用《通过格点屏蔽巨幅降低IGMH可视化分析片段间相互作用的耗时》（//www.umsyar.com/756）介绍的方法，对mIGM也同样适用，往往可以节约好几倍的格点数据计算耗时。 由于mIGM非常快，因此对不太大的体系算分子动力学过程产生的几百、几千帧结构都比较容易。基于这一点，我在Struct. Bond., 190, 297 (2026)中还将mIGM扩展到了用于可视化研究动态环境中的平均弱相互作用的情况，称为amIGM方法。这个方法远远比《使用Multiwfn研究分子动力学中的弱相互作用》（//www.umsyar.com/186）中介绍的目的相似的aNCI方法好用，在《使用amIGM方法图形化直观展现动态过程中的平均弱相互作用》（//www.umsyar.com/759）中对amIGM做了非常详细的介绍并且给出了丰富的在Multiwfn中计算的例子，非常欢迎阅读！   //www.umsyar.com/754 2025-11-20T01:04:06+08:00 Multiwfn的格点数据平移功能介绍 Introduction to grid data translation function in Multiwfn 文/Sobereva@北京科音  2025-Nov-19 1 前言 很多人在对周期性体系绘制三维函数（如电子密度差）的等值面时，容易遇到一个情况是发现感兴趣的等值面在晶胞的边缘，导致等值面被截断，无法观看完整。虽然VMD程序在显示格点数据等值面的时候可以要求把周期镜像显示出来，从而试图让位于晶胞边缘的等值面看起来完整，但等值面在晶胞边界的位置会不连贯，有个突变，因而还是不好看。还有一个解决方法是先用《Multiwfn中非常实用的几何操作和坐标变换功能介绍》（//www.umsyar.com/610）里提供的功能令晶胞进行平移使得感兴趣的区域在盒子中央，然后再重新计算格点数据，但这个过程不仅麻烦而且还花费重算一次的时间。 为了完美地解决以上问题，在2025-Nov-19更新的Multiwfn中增加了一个新功能，可以令格点数据连带着原子坐标在盒子的第1、2、3个轴上分别平移特定百分比，从而将感兴趣的等值面移到便于观察的地方，整个过程瞬间完成。下面就通过两个实际例子演示一下。Multiwfn可以在官网//www.umsyar.com/multiwfn免费下载，不了解Multiwfn者建议看《Multiwfn FAQ》（//www.umsyar.com/452）。记录格点数据最常用的格式是cube，不了解的话建议看《Gaussian型cube文件简介及读、写方法和简单应用》（//www.umsyar.com/125）。 2 实例 北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/KFP）里我讲CP2K做水合电子的AIMD模拟的幻灯片中，给了三个时刻的自旋密度等值面的图片，如下所示。由图可见从2000 fs时开始形成水合电子，到了2100 fs时水合电子已经完全形成。然而2100 fs时的等值面在盒子最边上，看起来很不舒服。这个问题靠前述的Multiwfn的格点数据平移功能即可解决。 2100 fs时刻CP2K产生的这个体系的自旋密度的cube文件MD-SPIN_DENSITY-1_2100.cube在//www.umsyar.com/attach/754/file.rar中。我们先看一下等值面图，用VMD、VESTA、Multiwfn等观看都可以。此例用Multiwfn载入cube文件后，进入主功能0，把等值面数值设为0.001，并且点击show data range复选框要求显示格点数据盒子边框后就可以看到下图 显然，为了让等值面位于盒子中央便于观看，应该对体系在第1个轴（当前对应X轴）的负方向平移30%左右；在第2个轴（当前对应Y轴）的正方向平移约50%；在第3个轴（当前对应Z轴）上用不着平移。 关闭图形窗口回到Multiwfn主菜单，依次输入 13  //处理内存中的格点数据的功能 19  //平移格点数据 -0.3  //在第1个轴负方向平移30% 0.5  //在第2个轴正方向平移50% [回车]  //不在第3个轴方向上平移 此时屏幕上看到如下提示，显示了在各个方向上平移了多少个格点，以及平移矢量。注意这里都是按正值显示（由于周期性，前面输入-0.3等同于输入了0.7）。  Translate along the 1st axis by   62 grids  Translate along the 2nd axis by   45 grids  Translate along the 3rd axis by    0 grids  Translation vector:    8.630404    6.264003    0.000000 Angstrom 当前程序检测到在平移后有些原子露在了盒子外面，问你是否把它们卷到盒子里，这里输入y要求卷入。 现在格点数据就处理好了。可以选择选项0把格点数据导出为新的.cub文件，也可以直接选选项-2观看等值面，现在看到的图如下。可见等值面已完全在盒子中央了，非常容易考察水合电子的形态。 3 注意事项 此功能对任何函数的格点数据都可以用。格点数据可以是从cube、CHGCAR等Multiwfn支持的记录格点数据的文件中读取的，也可以是Multiwfn的主功能5等功能基于波函数文件直接计算出来并存在内存中的。 本文介绍的功能对于盒子是非正交的情况也可以照常用，平移的方向对应于实际三个轴的方向，只不过Multiwfn目前无法正确显示这种情况的等值面，应当用VMD、VESTA等程序显示。 前面的例子中，格点数据对应的盒子（即均匀分布的格点所处的范围）和做周期性第一性原理计算对应的盒子（对晶体体系来说也相当于晶胞）是完全一致的。如果格点数据的盒子范围和周期性计算用的盒子范围不对应，比如计算格点数据的区域只是整个体系中的一小块，那么使用本文的功能就没任何意义。 如果你的体系本来就没周期性，本文的功能虽然也能使用，但并没有任何实际意义。 如果本文的功能给你的研究带来了便利，发表文章时请引用Multiwfn启动时提示的原文。 //www.umsyar.com/750 2025-08-31T07:27:00+08:00 使用Multiwfn结合VMD绘制分子局部区域表面静电势的方法 The way of plotting electrostatic potential for local region of molecular surface using Multiwfn in combination with VMD 文/Sobereva@北京科音  2025-Aug-31 《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图》（//www.umsyar.com/443）介绍的方法已经是如今非常主流的绘制分子表面静电势的方法。最近有人问“请问用Multiwfn画分子表面的静电势分布时，怎么画某个区域或某个原子的静电势分布呢？”，针对这个局部区域表面静电势的绘制问题，我在本文介绍一下做法。阅读本文之前建议先阅读《使用Multiwfn结合VMD分析和绘制分子表面静电势分布》（//www.umsyar.com/196）了解利用Multiwfn的输出文件手动在VMD中绘制整个体系表面静电势的基本操作过程，里面的步骤在下文里也会用到。 读者务必使用2025-Aug-14及以后更新的Multiwfn版本，否则不具有本文提到的一些特性。如果对Multiwfn不了解，看《Multiwfn FAQ》（//www.umsyar.com/452）了解相关知识。本文用的VMD是1.9.3版，用其它版本后果自负。 本文以非常常见的分子苯酚为例进行演示，将要分别绘制对应苯环的所有碳原子区域和氧原子区域的分子表面静电势图。不了解给Multiwfn用的波函数文件怎么产生的话看《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）。 启动Multiwfn，载入苯酚的波函数文件，即Multiwfn自带的examples目录下的phenol.wfn，然后输入 12  //定量分子表面分析 0  //开始分析。当前的分子表面默认对应电子密度为0.001 a.u.的等值面，即Bader定义的真空下的范德华表面 分析结束后可以看到整个分子表面的静电势的最小值和最大值分别为-0.042242和0.083100。之后作图我们用-0.05到0.05 a.u.的色彩刻度范围。  Global surface minimum: -0.042242 a.u. at   1.464068   3.342877  -0.000761 Ang  Global surface maximum:  0.083100 a.u. at  -1.948624   3.093517   0.021792 Ang 在Multiwfn的后处理菜单中输入 12  //对特定片段做定量分子表面分析 1-6  //苯环上碳原子的序号 现在屏幕上显示了对这个局部分子表面区域做表面静电势分析得到的各种指标。关于Multiwfn是怎么把分子表面划分出对应特定原子片段的，看Multiwfn手册3.15.2.2节 y  //导出locsurf.pqr文件 5  //把分子结构导出为pdb文件 mol.pdb  //导出的文件名 pqr格式是知名的pdb格式的变体，locsurf.pqr里的每个原子对应于分子表面上的一个顶点。如Multiwfn在屏幕上的提示所示，pqr文件的倒数第三列按照格式定义本来是用来记录原子电荷的，而当前这一列用来记录被映射的函数值，即以a.u.为单位的静电势。locsurf.pqr中的残基号那一列为1和为0分别对应这个表面顶点属于和不属于自定义的片段。 将mol.pdb载入VMD，VMD Main窗口中Graphics - Representation里把Drawing Method设成CPK。再把locsurf.pqr载入VMD，Graphics - Representation里把显示它用的选择语句设成resid 1，然后Drawing Method设Points（如果之前开了GLSL记得关闭），Size设10左右，Coloring Method用Charge，在Trajectory标签页里把色彩刻度范围设成-0.05到0.05。Graphics - Colors - Color Scale里把Method设BWR，现在看到下图。如果想弄出来色彩刻度轴，参考//www.umsyar.com/443的说明。 由上图可见，苯环的分子范德华表面中碳原子对应的区域是淡蓝色，说明此处静电势略负，这来自于丰富的pi电子对静电势的负贡献。 下面绘制对应于氧原子的局部范德华表面的静电势。在Multiwfn主功能12的后处理菜单选择11 Output surface properties of each atom，然后后输入y导出locsurf.pqr。这个文件内容和之前的locsurf.pqr唯一差异是现在残基号直接对应于相应表面顶点对应的原子序号。 重新按上文方式绘制，只不过把之前用的locsurf.pqr改成本次的locsurf.pqr，并且把显示表面顶点的那个representation的选择语句改成serial 12，因为氧原子是12号原子。现在看到下图。可见这部分区域中间部分的蓝色很深，这来自于氧原子带显著的负电荷，同时其孤对电子离这部分区域较近。 下文是按照//www.umsyar.com/443的做法绘制的整个分子表面的静电势填色图。相较之下，明显上文的做法可以把自己真正感兴趣的区域的静电势专门展现清楚，而不受到其它区域的视觉干扰。 最后提醒一下，按照本文的方式绘制局部分子表面的静电势图，发文章时记得要按照Multiwfn启动时的提示对Multiwfn的原文进行正确引用。 //www.umsyar.com/747 2025-08-02T15:59:43+08:00 使用了Multiwfn发表的第19001~22000篇文章打包下载 文/Sobereva@北京科音  2025-Aug-2 下面是笔者登记了的第19001-22000篇使用了极为流行的Multiwfn波函数分析程序（//www.umsyar.com/multiwfn）发表的文章的pdf合集。pdf都是我手动下载的，由于疏忽漏了十几篇文章对应的pdf文件。 注：根据Google学术统计，Multiwfn的引用次数目前已经达到约35000。由于Google学术向我推送的信息不完整，因此我迄今提供的pdf合集里的并不是目前全部引用了Multiwfn的文章。 下载链接：https://pan.baidu.com/s/17NWaECclsnTnY7rih22ebw 提取码：fktj 文件使用分卷压缩，共14.2 GB，都下载后对任意一个压缩包进行解压即可。 本次汇总的文章的列表见//www.umsyar.com/attach/747/19001-22000.txt，或者看//www.umsyar.com/multiwfn/all_citation4.html和//www.umsyar.com/multiwfn/all_citation5.html里面的第19001~22000号文章。更早的文章列表看Multiwfn主页（//www.umsyar.com/multiwfn）的Cited by页面里的其它页面。有的文章标题太长，为了避免文件路径超过Windows限制，网盘里分享的pdf的文件名只保留了前面一部分。 这些文章可以视为Multiwfn的例子库，便于用户了解Multiwfn在 、第一性原理、分子模拟实际研究中能发挥的重要价值。如果想找某一类应用实例，可以利用Acrobat对指定文件夹里全部文档进行搜索的功能搜索相关关键词。 由于本文件包里的文章都是它们刚在期刊网站上刊登的时候我就下载的，所以文档中大多没有页码、卷号，且多数都是还没经过proofing的原始状态，大家只要到Google学术搜索里搜索标题就能找到文章最终版链接。 之前的引用Multiwfn的第1~1010篇文章可以从这里下载：http://bbs.keinsci.com/thread-2850-1-1.html 之前的引用Multiwfn的第1011~2001篇文章可以从这里下载：http://bbs.keinsci.com/thread-6754-1-1.html 之前的引用Multiwfn的第2002~3000篇文章可以从这里下载：http://bbs.keinsci.com/thread-11023-1-1.html 之前的引用Multiwfn的第3001~4000篇文章可以从这里下载：http://bbs.keinsci.com/thread-14181-1-1.html 之前的引用Multiwfn的第4001~5000篇文章可以从这里下载：http://bbs.keinsci.com/thread-16840-1-1.html 之前的引用Multiwfn的第5001~6000篇文章可以从这里下载：http://bbs.keinsci.com/thread-19818-1-1.html 之前的引用Multiwfn的第6001~7000篇文章可以从这里下载：http://bbs.keinsci.com/thread-22597-1-1.html 之前的引用Multiwfn的第7001~8000篇文章可以从这里下载：http://bbs.keinsci.com/thread-25497-1-1.html 之前的引用Multiwfn的第8001~9000篇文章可以从这里下载：http://bbs.keinsci.com/thread-27978-1-1.html 之前的引用Multiwfn的第9001~10000篇文章可以从这里下载：http://bbs.keinsci.com/thread-31224-1-1.html 之前的引用Multiwfn的第10001~11000篇文章可以从这里下载：http://bbs.keinsci.com/thread-34124-1-1.html 之前的引用Multiwfn的第11001~12000篇文章可以从这里下载：http://bbs.keinsci.com/thread-36387-1-1.html 之前的引用Multiwfn的第12001~13000篇文章可以从这里下载：http://bbs.keinsci.com/thread-38881-1-1.html 之前的引用Multiwfn的第13001~14000篇文章可以从这里下载：http://bbs.keinsci.com/thread-41495-1-1.html 之前的引用Multiwfn的第14001~16000篇文章可以从这里下载： http://bbs.keinsci.com/thread-46062-1-1.html 之前的引用Multiwfn的第16001~19000篇文章可以从这里下载： http://bbs.keinsci.com/thread-50824-1-1.html 再次顺带强调关于恰当引用Multiwfn的问题。在Multiwfn程序主页的首页和Download页面、程序手册、程序启动时屏幕上显示的信息、程序包内的诸多文件（Multiwfn quick start.pdf文档、How to cite Multiwfn.pdf文档、LICENSE.txt文件）、Multiwfn入门贴（//www.umsyar.com/167）、Multiwfn FAQ（//www.umsyar.com/452）等所有用户可以轻易看到的所有地方都非常非常非常明确强调如何对Multiwfn进行正确的引用，但使用Multiwfn的文章中不引用或者胡乱引用的现象依然十分严重（特别是中国、印度、巴基斯坦、伊朗作者的文章尤甚）！第19001~22000篇文章里没按要求恰当引用Multiwfn的文章多达555篇！经常是只提及Multiwfn但不引用，甚至居然在引用Multiwfn的地方引用的是和Multiwfn根本没有任何直接联系的其他人的文章（极为恶劣！）还有不少文章里作者明明在研究中用了Multiwfn做分析，在论文里竟然连Multiwfn的名字都不提。不按明确声明的方式恰当引用，对免费而且没有任何经费支持的程序的发展十分不利，同时也是严重缺乏学术道德的行为。借此机会再次强烈呼吁用户重视对Multiwfn的恰当引用。最恰当的引用方式见Multiwfn可执行文件包内的How to cite Multiwfn.pdf文档的说明。没有恰当引用Multiwfn的文章在前述的引用列表中都已经明确注上了Multiwfn was not properly cited!的字样进行公示，并且在本文的文件包里放在了“未恰当引用Multiwfn的文章”子目录中。 //www.umsyar.com/746 2025-06-26T16:54:00+08:00 Angew. Chem.上发表了全面介绍各种共价和非共价相互作用可视化分析方法的综述 文/Sobereva@北京科音  2025-Jun-26 原子、分子间的相互作用在化学体系中无处不在。化学体系中的相互作用可以按强度大体分为化学键相互作用和弱相互作用，也可以按作用特征大体分为共价和非共价相互作用。这些相互作用普遍在数值层面上进行研究，如研究相互作用能及其物理成份（能量分解）、几何参数、键级、AIM拓扑分析考察键临界点属性，等等。而如今，这些相互作用的可视化分析（visual analysis）日趋流行，在计算化学研究文章中被应用得越来越普遍、在文献中越来越常见。这些方法可以令化学家们十分直观生动地认识相互作用出现的位置、强度、本质，不仅应用十分便利，也对传统的基于数值的研究相互作用的方式提供了关键性的互补。这些方法在近年来还得到了快速发展，涌现了许多新的分析手段。显然，对这些可视化分析方法进行全面、系统的综述是极为重要且迫切的。 近期，北京科音自然科学研究中心（http://www.keinsci.com）的卢天应知名的Angew. Chem. Int. Ed.期刊的邀请，发表了名为Visualization Analysis of Covalent and Noncovalent Interactions in Real Space的综述，可在此访问：https://onlinelibrary.wiley.com/doi/10.1002/anie.202504895、https://pan.baidu.com/s/17x4Pm0LINmt9iHDBLXWhhw?pwd=86rv。若未订阅，也可以免费阅读此文在ChemRxiv上的预印版https://doi.org/10.26434/chemrxiv-2025-9t442（内容和正式版基本一样，写文章时请引用正式版）。 推荐使用Multiwfn程序做文中介绍的各种可视化分析时引用此综述！ 此综述是到目前为止最完整和系统介绍各种化学体系中相互作用可视化分析方法的文章，面向广大化学家，特别是计算化学家。文章的写作注重令内容清晰易懂、深入浅出，且引发广泛读者们的兴趣。通过此文，化学研究者们可以快速且充分了解各种可视化分析方法的基本思想、适用范畴，并且结合文中提供的丰富的例子，读者能够充分体会到它们的重要实用价值。在看过本文有了基本的知识框架后，若想进一步深究物理、算法细节等方面的知识以及了解更多实际应用，推荐再阅读此综述中引用的各种方法的原文和相关文章。 此综述介绍的可视化分析方法十分全面，并且囊括了最新进展。主要介绍了NCI（亦称RDG）、IRI、IRI-pi、IGM、IGMH、mIGM、aNCI、amIGM、Hirshfeld表面分析、静电势分析、范德华势分析、ELF、LOL、LOL-pi、电子密度拉普拉斯函数、变形密度、价层电子密度。值得一提的是，其中IRI、IRI-pi、IGMH、mIGM、amIGM、范德华势分析、价层电子密度分析都是此文的作者所提出的。 综述中还简要提及了文章作者开发的十分流行的波函数分析程序Multiwfn（主页：//www.umsyar.com/multiwfn），这是唯一可以实现此文介绍的所有方法的程序，有很多方法还是此程序独家支持的。Multiwfn对于大量相互作用可视化分析方法在计算学领域的普及和推广起到了关键性的作用，此综述里几乎所有图都是基于Multiwfn计算的数据绘制的。Multiwfn的相关信息见《Multiwfn FAQ》（//www.umsyar.com/452）和介绍文章J. Chem. Phys., 161, 082503 (2024) DOI: 10.1063/5.0216272。 此外，//www.umsyar.com/667介绍的卢天在2024年发表的Visualization Analysis of Weak Interactions in Chemical Systems（Comprehensive Computational Chemistry Vol.2中的一章）和前述的综述有极为密切的联系，此文介绍的方法专注于弱相互作用，涉及的方法更少，而由于有明显更多的篇幅，因而对它们的介绍明显更深入细致，故和上面的综述有极大的互补性，强烈建议阅读！ 还十分推荐阅读以下文章及其中引用的文章获取更多相关知识以及计算操作细节： 使用mIGM方法基于几何结构快速图形化展现弱相互作用 //www.umsyar.com/755（http://bbs.keinsci.com/thread-57066-1-1.html） 使用amIGM方法图形化直观展现动态过程中的平均弱相互作用 //www.umsyar.com/759（http://bbs.keinsci.com/thread-57317-1-1.html） 使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用 //www.umsyar.com/621（http://bbs.keinsci.com/thread-28147-1-1.html） 使用IRI方法图形化考察化学体系中的化学键和弱相互作用 //www.umsyar.com/598（http://bbs.keinsci.com/thread-23457-1-1.html） 使用Multiwfn做Hirshfeld surface分析直观展现分子晶体和复合物中的相互作用 //www.umsyar.com/701（http://bbs.keinsci.com/thread-43178-1-1.html） 谈谈范德华势以及在Multiwfn中的计算、分析和绘制 //www.umsyar.com/551（http://bbs.keinsci.com/thread-17271-1-1.html） 在Multiwfn中单独考察pi电子结构特征 //www.umsyar.com/432（http://bbs.keinsci.com/thread-10610-1-1.html） ELF综述和重要文献小合集 http://bbs.keinsci.com/thread-2100-1-1.html 静电势与平均局部离子化能相关资料合集 http://bbs.keinsci.com/thread-219-1-1.html 使用Multiwfn研究分子动力学中的弱相互作用 //www.umsyar.com/186 使用Multiwfn做aNCI分析图形化考察动态过程中的蛋白-配体间的相互作用 //www.umsyar.com/591（http://bbs.keinsci.com/thread-21826-1-1.html） 通过电子定域化函数(ELF)、价层电子密度分析讨论亲核进攻的方向 //www.umsyar.com/606（http://bbs.keinsci.com/thread-23982-1-1.html） 通过价层电子密度分析展现分子电子结构 http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB201709252 此外，笔者的18碳环及相关体系的研究文章中充分利用了文中介绍的许多方法，是很好的应用实例，研究汇总见//www.umsyar.com/carbon_ring.html。 波函数分析与Multiwfn程序培训班（http://www.keinsci.com/WFN）对前述的综述中几乎所有方法都有十分全面、透彻、系统的介绍，并对各种分析方法在Multiwfn中的实际操作做了细致的演示并传授大量技巧，信息量远比综述里的大得多得多，因而很推荐对这些分析方法感兴趣者关注和参加。 下面是Visualization Analysis of Covalent and Noncovalent Interactions in Real Space综述中的图，供预览。                                             //www.umsyar.com/743 2025-06-01T10:25:46+08:00 使用NICS和磁感生电流考察芳香性时的一些易被忽视的重要问题 文/Sobereva@北京科音   2025-Jun-11 0 前言 NICS和磁感生电流是非常常用的磁响应性质一类的考察芳香性的方法，在《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）里有简要介绍并里面附了大量我写过的相关博文，而在“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/WFN）的“芳香性与离域性分析”一节里有极为全面深入的介绍。很多人对这两个方法在衡量芳香性方面的认识局限在： (1)NICS(1)ZZ约为0说明环是非芳香性，明显为负说明环是芳香性且越负芳香性越强，明显为正说明环是反芳香性且越正反芳香性越强。注意NICS的形式很多，这里仅以最常用的NICS(1)ZZ举例 (2)感生电流的方向如果满足左手定则（diatropic电流），电流越大芳香性越强；如果与左手定则方向相反（paratropic电流），电流越大反芳香性越强。电流大小一般通过对键截面上的电流密度做积分来衡量，称为bond current strength (BCS) 但实际上，仅光凭以上知识讨论芳香性或反芳香性，在一些情况下会得到不正确甚至严重错误的结论。在本文笔者介绍一些大多数人都普遍没有注意到的用NICS和感生电流讨论芳香性方面的重要问题。把这些问题都彻底搞清楚了，基于磁属性研究芳香性时就不惧审稿人comment了。 本文中的计算用的几何结构，若无专门提及，皆是B3LYP/6-31G*下优化的无虚频结构。 1 NICS、感生电流与HOMO-LUMO gap的关系 感生电流和NICS计算公式都涉及到对占据轨道和空轨道的循环，并且空轨道与占据轨道能量之差出现在分母项中，这一点看《深入理解分子轨道对磁感生电流的贡献》（//www.umsyar.com/703）里的公式1、2、3便知。总的来说，HOMO-LUMO gap越大，空轨道与占据轨道整体能量差越大，导致NICS和感生电流越接近0。因此，如果两个体系的gap相差较大，不应直接凭NICS或BCS判断二者的芳香性强弱，还应该把gap差异对结果的影响考虑进去。在《深度揭示互为等电子体的苯、无机苯和carborazine的芳香性的显著差异》（//www.umsyar.com/731）里介绍的笔者的Chem. Eur. J., 30, e202403369 (2024)一文中就注意到了这一点，文中研究的carborazine的HOMO-LUMO gap（7.30 eV）明显小于苯（11.20 eV）和无机苯（12.18 eV），因此文中提到：It is noteworthy that the aromaticity of carborazine compared to benzene and borazine is more or less quantitatively overestimated by the value of BCS and NICS。 4n电子的Huckel反芳香性体系的HOMO-LUMO gap比4n+2的Huckel芳香性体系的整体更小，Chem. Soc. Rev., 44, 6597 (2015)认为这是前者的BCS比后者更大的原因。例如反芳香性的环丁二烯和芳香性的苯在B3LYP/6-31G*级别算的HOMO-LUMO gap分别为3.692和6.800 eV，使用《使用SYSMOIC程序绘制磁感生电流图和计算键电流强度》（//www.umsyar.com/702）介绍的SYSMOIC算出来的环丁二烯的BCS为-0.7941（长C-C键）和-0.7499 a.u.（短C-C键），而苯的C-C键的BCS大小则小得多，只有0.4216 a.u.。 如《正确地认识分子的能隙(gap)、HOMO和LUMO》（//www.umsyar.com/543）所说，HF成份越低的泛函算出来的HOMO-LUMO gap越低，因此若发现HF成份低的泛函算出来的NICS和BCS的绝对值较大，大概率也是上述原因（但并非HF成份越低的泛函算出来的总是越大，因为泛函的差异还影响到其它方面）。 2 NICS的环尺寸的依赖性 原理上，感生电流比NICS判断芳香性更合理，因为NICS对环尺寸有依赖性。虽然有的文章说NICS的环尺寸依赖性弱，但那往往只是对较小的环之间对比做的结论，再加上本身不同体系的芳香性就有所不同，令那种说法更不确切。 电流对特定位置产生的磁场可以用Biot–Savart方程得到。如果是环形电流，则在环中心产生的磁场为μ0*I/(2*R)，其中μ0是真空磁导率，I是电流大小，R是环的半径，更多信息见https://en.wikipedia.org/wiki/Biot%E2%80%93Savart_law。因此，对芳香性的情况，当感生电流相同时，对应芳香性的环形路径整体半径越大、原子距离环中心整体越远，在环中心产生的磁屏蔽效果就越弱、NICS的绝对值越小。如果认为感生电流I的大小是芳香性强弱的准确体现，那么显然用NICS对比环尺寸明显不同的体系的芳香性的大小就是不甚严格的，因为对越大的环会越低估芳香性。所以，如果要基于磁属性对比明显不同大小的环的芳香性，建议用感生电流而非NICS说事。对于反芳香性，也有和上面完全类似的情况。 3 感生电流对积分截面的依赖性 对应芳香性的环上，如果只考虑对芳香性有贡献的电子的感生电流，原理上来说在垂直于键（更确切来说是垂直于电流方向）的各个截面上积分电流密度得到的BCS都应该是相同的。但实际得到的BCS没有这么理想化，环上的各个键的BCS往往不完全相同，在于： (1)积分的数值误差、有限的积分截面尺寸 (2)积分截面可能并不严格垂直于电流方向 (3)在一些键上存在方向彼此相反的感生电流的抵消作用，此效应对于多环芳烃等并环的情况尤为明显 (4)环周围的原子的感生电流产生的干扰，以及环上原子的定域化的电子（内核电子、sigma电子、孤对电子等）产生的局部感生电流添乱，这导致计算的BCS对键上的积分截面位置的选取可能存在不可忽视的依赖性 (5)基组不够完备，导致感生电流计算不够准确。对D2h点群的环丁二烯（反芳香性）极小点结构，用SYSMOIC在B3LYP/6-31G*下计算的短和长C-C键的BCS分别为0.7499和0.7941 a.u.，而基组用完备得多的def2-QZVP计算结果分别为0.7116和0.7111 a.u.，可见后者明显更等价 若发现要考察芳香性的环上的各个键的BCS不等价性不可忽视，可以取感生电流较明确、被上述的(4)干扰小的键的BCS作为环电流值衡量芳香性，或者对环上的键的BCS取平均作为环电流值（适合碳环、环丁二烯等情况）。此外，也可以尝试计算BCS时只考虑对芳香性有贡献的电子所占据的分子轨道，例如用感生电流研究pi芳香性时就只考虑pi轨道。 4 感生电流判断芳香性和反芳香性方面相对于NICS的优势 NICS很有可能给出关于芳香性错误的结论，不如观看和积分感生电流，以及观看《通过Multiwfn绘制等化学屏蔽表面(ICSS)研究芳香性》（//www.umsyar.com/216）里介绍的ICSS图来得靠谱，这一节给出几个例子。 JCTC, 6, 1131 (2010)发现NICS用于考察氟化氢三聚体(HF)3的芳香性有问题。在wB97XD/def-TZVP下优化出的此体系的结构如下所示。凭基本常识就知道这个体系是非芳香性的，但是B3LYP/def2-TZVP算出来的NICS(0)ZZ为13.52 ppm，仿佛有显著反芳香性似的。不过NICS(1)ZZ为-0.04 ppm，这倒是正确体现出此体系并没有反芳香性。   为什么(HF)3的NICS(0)ZZ那么正？这可以从此体系所在平面上的感生电流图上理解，如下图所示（磁场垂直于屏幕向外）。可见此体系并没有形成环绕整体的感生电流，但由于存在环绕各个HF的sigma键的dia电流，导致在环中央区域呈现出等效的para电流的特征（红色箭头所示），使得环中央是去屏蔽的，因而NICS(0)ZZ为正。而由于这种效应在垂直于体系平面方向上衰减很快，因此NICS(1)ZZ几乎为0。这也体现出，对于小环体系，由于sigma键的感生电流产生的磁场离环中心太近，NICS(0)ZZ用来衡量芳香性是极差的做法。但对于18碳环及其等电子体B6N6C6这样的大环倒不是什么问题，因为局域的电子对应的局部的感生电流造成的磁场在环中心处都基本衰减没了，这也是为什么《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（//www.umsyar.com/696）介绍的文章中用了NICS(0)ZZ讨论。   实际上对于苯分子也有类似情况。下图是苯分子的pi以外的电子产生的感生电流图，可见由于局域的sigma电子产生的局部dia电流的存在，在环中心也有等效的para电流（红色箭头所示）并使得环中心产生磁去屏蔽效应。不过由于苯的pi电子造成的全局的dia电流在环中心产生的磁屏蔽效应远强于para电流的磁去屏蔽效应，因此苯的NICS(0)ZZ还是负值。   我极其推荐阅读《使用Multiwfn绘制一维NICS曲线并通过积分衡量芳香性》（//www.umsyar.com/681）中的讨论，了解sigma和pi电子对NICS_ZZ的贡献在垂直于环方向上的变化特征。 《18个氮原子组成的环状分子长什么样？一篇文章全面揭示18氮环的特征！》（//www.umsyar.com/725）里介绍的我的ChemPhysChem, 25, e202400377 (2024)文章中，对18氮环的最稳定构型计算了NICS(0)ZZ，为3.42 ppm，也是表面上体现出了轻微的反芳香性。然而通过ICSS等值面图和AICD感生电流图，都充分确认了此体系没有芳香性或反芳香性特征，详见博文里的图。这也是个很好地体现NICS(0)ZZ展现（反）芳香性不靠谱的例子。 NICS亦有可能把非芳香性体系误判成芳香性的。Theor. Chem. Acc., 134, 8 (2015)指出NICS衡量很多过渡金属团簇的芳香性也是不合理的。其中一个例子是如下所示的Ti3(CO)3。基于PCCP, 13, 4576 (2011)的SI里给出的BP86/6-311+G(d)级别下优化的结构，我用PBE0/def2-TZVP算出的NICS(0)ZZ和NICS(1)ZZ分别是-141.1和-11.7 ppm，乍看起来明显芳香性。   下图左侧和右侧分别是Ti3(CO)3的分子平面上和分子平面上方1埃处的感生电流，数值过大的部分截掉了，粉色大圆球是Ti，磁场垂直屏幕向外。由图可见，无论是在分子平面上还是在其上方1埃处，都有显著的环绕Ti原子的para电流，很大程度等效带来了在环中心及其上方的dia电流的效果（红色箭头所示），即令环中心及其上方一定区域都是磁屏蔽的，这是NICS(0)ZZ和NICS(1)ZZ都为负值的原因。然而从图上并没有看到环绕整体的连贯的感生电流，因此Ti3(CO)3就是非芳香性的。可见此例不仅NICS(0)ZZ，连NICS(1)ZZ都误判了芳香性。   5 显著且连贯的dia感生电流未必能证明有明显的芳香性 从上一节看到，感生电流判断芳香性比个别位置的NICS更为严格且直观，但即便是感生电流，也未必总能正确衡量芳香性，这里以[N6H6]2+体系为例进行说明，更多讨论见ChemistrySelect, 2, 863 (2017)。下图展示了此体系的平面构型   此构型并不是极小点，而是有很多破坏平面的虚频，并且其能量甚至比其开环后的线型结构的还要高，因此明显没有芳香稳定化作用。用Multiwfn计算出它的多中心键级几乎为0（0.00035），也进一步体现出没有pi共轭。如果用《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）的做法绘制它的LOL-pi图，也会看到各个氮都有定域的孤对电子，根本没有像苯的pi电子一样在环上充分离域。然而，下图所示的[N6H6]2+平面构型的分子上方1埃处的感生电流（磁场垂直于屏幕向外）则是完全顺时针的dia电流，和具有显著芳香性的苯的情况看起来完全一样！   为什么[N6H6]2+的平面构型不是极小点？这是因为如下图的轨道占据情况所示，尽管此体系表面上有10个pi电子，看似满足Huckel的4n+2芳香性规则因而芳香性会使此构型被稳定化，但e2u对称性的轨道的反pi特征显著（对于苯来说是没有占据的），因此pi共轭作用很弱，再加上如此众多的孤对电子间的强烈互斥作用，显然令其极不稳定。   为什么[N6H6]2+的平面构型像普通芳香性体系一样有显著的贯穿整个环的dia电流？这需要利用《深入理解分子轨道对磁感生电流的贡献》（//www.umsyar.com/703）里介绍的知识，没看过者强烈建议认真阅读。此文说了，感生电流可视为由各种占据轨道向各种空轨道跃迁所贡献，上图蓝色箭头标注的跃迁是此结构下对感生电流贡献最大的跃迁（因为这种跃迁能最小），其中非占据的b2g轨道比占据的e2u轨道轨道多一个节面，这种跃迁是平动允许的，因此对dia电流有贡献。这和造成苯分子的dia电流的机制是完全相同的，只不过苯分子的dia电流主要是e1g到e2u跃迁所贡献的。 由此例可见，光靠是否存在全局的dia感生电流判断芳香性并不总是可靠！因为这种电流是否存在，和体系是否真的有芳香性特征并没有内在的必然联系，而与它关系真正最直接、最紧密的是轨道的对称性特征！ 可能大多数人觉得，有芳香性的体系有满足左手定则的dia电流的原因在于：某个环有芳香性意味着环上有明显离域的电子，因此在外磁场的作用下会像经典电磁学描述的导体一样在环上产生dia电流。但实际上真正本质并非如此，不能忽视关键性的量子力学效应！一般的芳香性体系有dia电流实际上是一个表象，在于它们的轨道对称特征和跃迁正好能满足产生显著dia电流的条件，而能产生显著dia电流的并非只有芳香性体系。因此，要这样认为：有显著的贯穿整个环的dia电流是这个环具有芳香性的必要非充分条件。 由于磁屏蔽特征是感生电流的进一步衍生属性，[N6H6]2+的平面构型既然都和苯一样具有相同感生电流特征了，其自然NICS(0)ZZ和NICS(1)ZZ也都为负，分别为-178.1和-24.6 ppm。 此例彰显出判断芳香性最好不要光基于磁属性，即便这种判断芳香性的做法极为流行。同时结合多中心键级等基于电子结构的方法，以及结合芳香稳定化能或离域化能这种能量层面的方法，结论会更可靠、更有说服力。本节的[N6H6]2+的例子也体现出，形式上满足4n+2规则的体系未必真的有芳香性，若电子都没法在环上充分离域、不构成共轭体系，满足4n+2也仍是非芳香性的。 6 显著且连贯的para感生电流未必能证明有明显的反芳香性 J. Org. Chem., 88, 14831 (2023)通过一系列分子指出，即便环上存在显著且连贯的para电流，也未必能可靠地证明这个环就有明显的反芳香性。这一节就用其中最简单的体系，苯+1电荷的阳离子自由基（具有5个pi电子）和-1电荷的阴离子自由基（具有7个pi电子）来说明这一点，这种体系的（反）芳香性没法基于经典的Huckel规则判断。 如下所示，JOC文中给出了苯阳离子、中性苯、苯阴离子的分子平面上方1 Bohr处的感生电流，磁场垂直于屏幕向外，计算用的是对各自状态优化的几何结构，电场强度越大箭头越蓝、越小越红。可见，苯阳离子和苯阴离子在碳环上方的电流都是逆时针的para电流，与中性苯的情况相反，乍看起来苯的阴、阳离子都是显著反芳香性的。   然而，JOC这篇文章计算了苯的阴、阳离子的芳香稳定化能（用到的反应见文章的SI），发现这两个体系是有芳香稳定化作用的，因此从能量角度来说是有芳香性的。我还通过Multiwfn在文中给的苯离子的极小点结构下在B3LYP/6-31G*级别计算了六中心键级，苯的阳离子和阴离子的结果分别为0.044和0.045，只有中性苯的0.086的约一半。具体来说，苯阳离子和阴离子的六中心键级分别几乎来自于alpha和beta电子的贡献（且和中性苯的相应自旋电子的贡献差不多），而另一个自旋几乎没贡献。也就是说，相对于苯，离子状态下电子数改变1的那个自旋的电子对芳香性没明显贡献，而电子数未变的自旋的电子对芳香性的贡献基本保持不变。 由于从能量和电子结构角度都确认了苯的阴、阳离子是芳香性的，那么为什么它们的感生电流是para的，显得像普通反芳香性体系一样？JOC文中给了下面的图对感生电流进行了解释。由于Jahn-Teller效应，阳离子（5π）和阴离子（7π）体系的对应于苯的二重简并的pi轨道发生了轻微分裂。图中蓝色大箭头标注了对感生电流产生主导的跃迁。对苯离子体系的情况，由于能级分裂造成的能量差较小，明显小于1→2的跃迁，因此对阳离子和阴离子分别是1→1和2→2跃迁对感生电流产生了主导性的贡献。由于这种跃迁涉及的两个轨道的节面数相同，是转动允许的，因此如《深入理解分子轨道对磁感生电流的贡献》所述，贡献的是para电流。   由本节和上一节的例子可见，虽然感生电流对大多数体系判断芳香性是合理的，而且比起只看个别位置的NICS_ZZ可靠不少，但在极个别体系上还是能暴露出其局限性。因此若有可能，除了基于磁属性分析外，还是建议再结合一两种靠谱的手段，特别是可靠又容易计算的多中心键级，来衡量芳香性，减少误判芳香性、被审稿人怼的可能（本来芳香性这方面的“说法”和“观点”就相当多）。 但也没必要盲目用过多的方法衡量芳香性，否则一方面太冗余，另一方面本身有的芳香性指标就不那么靠谱、普适，纳入讨论反倒无益于揭示芳香性的真相。有文章认为键的均衡化、出现dia和para的感生电流，属于芳香性的二阶（次级）特征，而能量的稳定化和电子的离域，才算是芳香性的一阶（即最关乎本质的）特征，这个观点我是认可的。 //www.umsyar.com/740 2025-04-05T06:31:00+08:00 利用Multiwfn令Dalton计算时使用其它程序产生的轨道作为初猜 Using Multiwfn to take orbitals generated by other programs as initial guess in Dalton calculations 文/Sobereva@北京科音  2025-Apr-5 1 前言 Dalton这个 程序虽然用起来复杂，但有很多独特的亮点，例如可以在SOPPA(CCSD)和CASSCF级别下计算各种磁属性（笔者最近在Chem. Eur. J., e202404138 (2025)文中就用Dalton在CASSCF下算了NICS），可以用响应方式较准确地计算磷光寿命，可以算双光子吸收（北京科音高级 培训班http://www.keinsci.com/KAQC里专门有一节讲了)，可以用高级别耦合簇方法算激发态（见《使用Dalton通过CC3方法极高精度计算激发态》//www.umsyar.com/463），等等。可惜Dalton的SCF收敛方面做得一般，比不上Gaussian，而且单靠Dalton无法产生一些重要的轨道，比如UHF计算产生的UNO轨道经常被用于CASSCF的初猜，但Dalton甚至连UHF都做不了，因为Dalton里开壳层只支持RO形式。 为解决以上问题，从2025-Apr-5更新的Multiwfn开始，Multiwfn产生Dalton输入文件的功能支持了将Multiwfn中当前的轨道波函数写入到Dalton的输入文件.dal里作为初猜波函数。因此借助Multiwfn，Dalton用户可以用Gaussian、ORCA、GAMESS-US等诸多主流的 程序产生的波函数当初猜，从而解决Dalton遇到SCF不收敛的问题。而且，Dalton用户还可以直接用Multiwfn产生的UNO、定域化分子轨道（LMO，见//www.umsyar.com/380）等当CASSCF初猜。 下面就通过例子简单演示一下具体实现方法，本文中涉及的文件都可以在//www.umsyar.com/attach/740/file.zip下载。Multiwfn可以在//www.umsyar.com/multiwfn免费下载，相关常识看《Multiwfn FAQ》（//www.umsyar.com/452）。如果读者还不知道Dalton怎么安装、运行、产生输入文件的话，看《 程序Dalton的编译方法和运行方式简介》（//www.umsyar.com/463）。本文用的Multiwfn是2025-Apr-5更新的版本，Gaussian是G16 B.01，Dalton是2022版。 如果本文介绍的Multiwfn的功能给你的研究带来了帮助，请在论文中引用Multiwfn刚启动时提示的Multiwfn程序的原文，这是对Multiwfn开发和维护最好的支持！ 2 把Gaussian收敛的HF波函数当Dalton做CCSD(T)的初猜 这一节以CH3NH2为例演示把Gaussian收敛的HF/cc-pVTZ波函数当Dalton做CCSD(T)/cc-pVTZ的初猜。首先用Gaussian做一个普通的单点计算，输入文件是本文文件包里的CH3NH2.gjf，内容如下所示。算完后把得到的CH3NH2.chk用formchk转成CH3NH2.fch。 %chk=D:\CH3NH2.chk # HF/cc-pVTZ int=NoBasisTransform [空行] test [空行] 0 1  C                  0.05159500    0.70381800    0.00000000  H                  0.59439800    1.06209400    0.88166000  H                  0.59439800    1.06209400   -0.88166000  H                 -0.94293100    1.18451300    0.00000000  N                  0.05159500   -0.76078400    0.00000000  H                 -0.45830300   -1.10306000    0.81258800  H                 -0.45830300   -1.10306000   -0.81258800 [空行] [空行] 之后启动Multiwfn，载入CH3NH2.fch，依次输入 100  //其它功能（Part 1） 2  //产生输入文件 19  //产生Dalton输入文件 CCSDpT.dal  //要产生的输入文件名 [回车]  //要产生的.mol文件名用默认的CH3NH2.mol y  //由于当前内存里有基函数信息，Multiwfn问你是否把轨道展开系数写入到.dal文件里 1  //使用Dalton标准格式（4F18.14）方式写入。即四个值换一行，每个值占18列、小数位数14个 此时当前目录下就有了CCSDpT.dal和CH3NH2.mol。把CH3NH2.mol里的基组名都替换为cc-pVTZ_emsl。目前CCSDpT.dal里的计算设置对应B3LYP算单点，自行把关键词改为CCSD(T)计算用的，此时CCSDpT.dal的内容如下 **DALTON INPUT .RUN WAVE FUNCTIONS **WAVE FUNCTIONS .CC *CC INPUT .CC(T) *ORBITAL .MOSTART     //代表从**MOLORB部分读取轨道展开系数作为初猜  FORM18 **MOLORB (Punched by Multiwfn)  -0.00001750284370  0.00060444251600  0.00013849165900 -0.00001024308040   0.00002442752580 -0.00001726827680 -0.00014237439200  0.00004058047020 ...略 **END OF INPUT 用Dalton进行计算，可以看到SCF一轮就收敛了，达到了我们的目的： @ *** DIIS converged in   1 iterations ! @     Converged SCF energy, gradient:    -95.252488149534    1.44D-07     - total time used in SIRFCK :              0.00 seconds 之后CCSD(T)计算也顺利跑完了。 一些注意事项和相关信息： • 对于广义或部分广义收缩基组，都应当用int=NoBasisTransform关键词，拿不准就始终带上。详见《将Gaussian等程序收敛的波函数作为ORCA的初猜波函数的方法》（//www.umsyar.com/517）里对此关键词的解释。 • 把Gaussian的cc-pVDZ、cc-pVTZ的波函数当Dalton的初猜时，Dalton要用cc-pVDZ_emsl、cc-pVTZ_emsl基组，而不要用不带_emsl后缀的，否则还是需要迭代很多次才能收敛（甚至可能不收敛），但带不带_emsl的版本收敛后的结果是相同的。而如果是def、def2系列基组，没什么特别要注意的，比如Gaussian里def-TZVP和def2-TZVP做的计算，在Dalton里基组名分别写成Turbomole-TZVP、def2_tzvp即可，注意对大小写有要求，必须对应Dalton目录下的basis子目录里的基组文件名。 • Dalton默认是基于球谐型基函数做的计算，因此如果在Gaussian里使用的是比如6-31G*等D壳层默认是笛卡尔型基函数的基组，应当写上5d关键词强行要求用球谐型，详见《谈谈5d、6d型d壳层基函数与它们在Gaussian中的标识》（//www.umsyar.com/51）。虽然Dalton也可以用笛卡尔型基函数，但Multiwfn往Dalton输入文件里写入初猜轨道系数的功能不支持此情况。 • 只要一个 程序能输出Multiwfn可以识别的带有基函数信息的波函数文件，并且用的是球谐型基函数，则波函数都可以通过以上方式作为Dalton的初猜。支持的文件详见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）。例如ORCA做计算产生的molden文件也可以载入Multiwfn后像上面那样产生带初猜波函数的Dalton输入文件。 3 把Gaussian收敛的DFT波函数当Dalton做DFT的初猜 此例以顺铂体系为例，把Gaussian收敛的PBE0波函数当Dalton做PBE0计算的初猜，对Pt用Stuttgart小核赝势，对配体用6-311G*。Gaussian输入文件如下（本文文件包里的Ptcoord.gjf）： %chk=D:\Ptcoord.chk #P PBE1PBE/genecp int=NoBasisTransform [空行] b3lyp/def2TZVP opted    //众所周知这是纯给自己看的标题行，别问为什么写成这个 [空行] 0 1 [坐标略...] [空行] N H Cl 6-311G* **** Pt SDD **** [空行] Pt SDD [空行] [空行] 按照和上一节例子相同的过程，把chk转成fch，再用Multiwfn产生.dal和.mol文件。.dal里的泛函名改成PBE0，然后设置.mol里的基组和赝势成为下面这样，当前对Pt用的赝势和赝势基组和Gaussian里的情况一样（若不确定的话，可以Gaussian里用gfinput关键词输出定义，并和Dalton的基组、赝势目录下的相应文件对照）。 ATOMBASIS test molecule Generated by Multiwfn Atomtypes=4 Angstrom Nosymmetry charge=0 Charge=1.0 Atoms=6 Basis=6-311G* H1   -0.82596600    1.64390200    2.14978600 H2    0.00000000    2.40772900    0.93575300 H3    0.82596600    1.64390200    2.14978600 H4   -0.82596600   -1.64390200    2.14978600 H5   -0.00000000   -2.40772900    0.93575300 H6    0.82596600   -1.64390200    2.14978600 Charge=7.0 Atoms=2 Basis=6-311G* N1    0.00000000    1.59755500    1.56108400 N2   -0.00000000   -1.59755500    1.56108400 Charge=17.0 Atoms=2 Basis=6-311G* Cl1    0.00000000    1.70827400   -1.36819100 Cl2   -0.00000000   -1.70827400   -1.36819100 Charge=78.0 Atoms=1 Basis=stuttgart_rsc_1997_ecp ECP=stuttgart_rsc_1997_ecp Pt1   -0.00000000   -0.00000000    0.18195700 用Dalton计算后，经过4轮SCF就轻易收敛了，结果为-1152.597056，和Gaussian给出的-1152.596983非常接近。而如果此例Gaussian和Dalton都用HF的话，会和上例一样，结果精确相同且SCF一轮可收敛。 4 把UNO轨道当Dalton做CASSCF的初猜 之前在《CASSCF计算双自由基以及双自由基特征的计算》（//www.umsyar.com/264）里示例过用Gaussian基于UNO初猜轨道做CASSCF(2,2)计算的方法，不了解这种计算和UNO概念的话建议先看看。这一节演示使用Gaussian产生UHF波函数，读入Multiwfn并产生UNO轨道，然后写入.dal文件用于Dalton做CASSCF(2,2)计算的初猜轨道的流程。示例体系还是264号博文里的C4H8双自由基，基组还是此文里的6-31G*。 首先用以下Gaussian输入文件（本文文件包里的C4H8.gjf）对C4H8做UHF计算得到对称破缺波函数。不了解为什么用guess=mix的话看《谈谈片段组合波函数与自旋极化单重态》（//www.umsyar.com/82）。 %chk=D:\C4H8.chk # UHF/6-31G* guess=mix 5d [空行] ub3lyp/6-31g(d) opted [空行] 0 1 C     -0.74742092    1.77656753    0.00000000 H     -0.62438907    2.32965189    0.92333358 H     -0.62438907    2.32965189   -0.92333358 C     -0.74742092    0.30933780    0.00000000 H     -1.24978876   -0.09122321    0.88294562 H     -1.24978876   -0.09122321   -0.88294562 C      0.74742092   -0.30933780    0.00000000 H      1.24978876    0.09122321   -0.88294562 H      1.24978876    0.09122321    0.88294562 C      0.74742092   -1.77656753    0.00000000 H      0.62438907   -2.32965189   -0.92333358 H      0.62438907   -2.32965189    0.92333358 [空行] [空行] 计算完毕后把C4H8.chk转成C4H8.fch，载入Multiwfn，然后依次输入 200  //其它功能（Part 2） 16  //基于fch文件里的密度矩阵产生自然轨道。这个功能在《在Multiwfn中基于fch产生自然轨道的方法与激发态波函数、自旋自然轨道分析实例》（//www.umsyar.com/403）里有专门的介绍 SCF  //读取fch里的SCF类型的密度矩阵 1  //产生空间自然轨道 此时从屏幕上看到如下信息，这便是UNO轨道的占据数。其中最高占据自然轨道（HONO）和最低非占据自然轨道（LUNO）的占据数分别是1.123225和0.876775，都偏离整数显著，很满足期望。  Occupation numbers:     2.000000    2.000000    1.999999    1.999999    1.999904    1.999898     1.999881    1.999835    1.999309    1.999120    1.999004    1.998977     1.996829    1.996536    1.996053    1.123225    0.876775    0.003947     0.003464    0.003171    0.001023    0.000996    0.000880    0.000691 ...略 接着输入 n  //不把自然轨道导出成new.mwfn文件并重新载入，因为当前不需要如此（如果你想看一下当前的UNO轨道图像的话，应选y，之后用Multiwfn主功能0观看轨道） 0  //返回主菜单 100  //其它功能（Part 1） 2  //产生输入文件 19  //产生Dalton输入文件 CASSCF.dal  //要产生的输入文件名 [回车]  //要产生的.mol文件名用默认的C4H8.mol y  //把当前内存里的轨道（即UNO）的展开系数写入到.dal文件里 1  //使用Dalton标准格式（4F18.14）方式写入 当前C4H8.mol里的基组默认就是6-31G*，不用改。修改CASSCF.dal，把以下两行 .DFT  B3LYPg 替换为CASSCF(2,2)计算的设置： .MCSCF *CONFIGURATION INPUT .CAS SPACE 2 .ELECTRONS 2 .INACTIVE 15 Dalton计算完毕后，活性空间内自然轨道占据数为1.2679和0.7320，能量为-156.024938 Hartree，和Gaussian在同级别做CASSCF精确一致。你若想看Dalton做CASSCF产生的自然轨道，把计算产生的.tar.gz文件里的molden.inp解压出来载入Multiwfn，并效仿《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）的方式观看即可。 5 关于线性依赖基函数的问题 较大基组（尤其是带弥散函数的），计算原子分布较致密的体系很容易出现线性依赖问题，Gaussian和Dalton此时都会自动去除线性依赖基函数，导致最终求解出来的轨道数会少于基函数数目。这种情况下，Gaussian输出信息中的NBsUse（实际用的基函数数目，等同于实际求解出来的轨道数）小于NBasis，Dalton输出文件中的Total number of orbitals小于Number of basis functions。这给上文的传轨道带来了麻烦。有两个解决方法： (1)不去除线性依赖基函数：Gaussian计算时带着IOp(3/32=2)关键词、Dalton计算时在**WAVE FUNCTIONS里面加上以下内容，从而让两个程序都不去除线性依赖基函数 *ORBITAL .AO DELETE 1E-100 这种做法的缺点是可能导致数值不稳定性。 (2)自动去除线性依赖基函数，但让Dalton少读取轨道：比如Gaussian计算显示NBsUse=50、NBasis=54，因此自动去除了4个线性依赖基函数。因此当Multiwfn载入fch文件后，最后4个轨道的展开系数都会为0，导出到.dal文件里也是如此。此时Dalton若照常读取全部轨道，计算一开始会报错*** ORTHO-FATAL ERROR ***。解决方法是在*ORBITAL里写上 .DELETE 4 这样Dalton就不会读取最后4个轨道了，计算可以正常进行。 6 .dal里的轨道展开系数存在星号而无法读取的问题 有时候Multiwfn写入到.dal里的轨道展开系数会存在一堆星号，导致Dalton读取轨道时提示input conversion error错误。这不是Multiwfn的问题，而是Dalton自身设计不周。Dalton定义的**MOLORB字段中记录展开系数的格式在前面说了，每个数值都是F18.14格式，因此当有的轨道展开系数特别大，超过了其整数位记录上限，就会被记录成一堆星号。这种问题在基组较弥散而又不让自动去除线性依赖基函数的情况下出现概率大。这是为什么Multiwfn导出Dalton输入文件的时候特意让你选择记录的格式，前面的例子都是选1用4F18.14格式，还可以选2用4E20.12格式，这是使用科学计数法输出，因此就不会出现数值过大记录成星号的问题了。相应地，必须自己修改Dalton的读取轨道展开系数的源代码。对于Dalton 2022版，需要修改Dalton源代码目录下的DALTON/sirius/sirgp.F文件，把PFMT = '(4F18.14)'替换为PFMT = '(4E20.12)'，之后重新编译即可。 7 其它 Multiwfn就相当一个大熔炉，能读取诸多 程序产生的波函数，又能导出波函数作为诸多 程序的初猜，这在《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）里专门说了。对于传轨道到ORCA的情况我之前还专门写过《将Gaussian等程序收敛的波函数作为ORCA的初猜波函数的方法》（//www.umsyar.com/517）。 MOKIT也是很好的可以实现不同程序间传轨道的工具，fch到.dal和.mol的转换可以用里面的fch2dal子程序。相对于MOKIT，本文介绍的用Multiwfn传轨道有几个额外的好处： (1)Multiwfn不需要安装MOKIT所需的Python环境。尤其是对于Windows用户，Multiwfn解压即用明显更为方便 (2)Multiwfn特意考虑了.dal里记录轨道展开系数的格式的问题并给了解决办法，上一节说了 (3)Multiwfn自身可以产生诸多类型的轨道（定域化分子轨道、空间/自旋自然轨道、NAdO、AdNDP、NTO、双正交化轨道等等）并输出作为Dalton的初猜 而fch2dal也有额外的优点，它会把fch里的基组、赝势的具体定义直接写到输出的.mol里，因此如果Gaussian计算时用的基组、赝势在Dalton里没有自带也能直接计算而不需要自定义。并且fch2dal支持用笛卡尔型基函数的情况。