文/Sobereva@北京科音 2025-Dec-9
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。
这里假定读者已经了解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中实现。
下面列举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)图的等值面范围更大,可以进一步看到配体在很多区域和蛋白质间也有显著的范德华主导的相互作用,对应绿色等值面区域。
这一节专门说一下用于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关闭此策略。
本节将会具体演示如何绘制出来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)。因为这些例子中计算格点数据的区域都没有离模拟用的盒子边界太近,因此不需要考虑周期性。
这个例子的轨迹文件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,都可以自己改。色彩刻度范围设得越小,等值面上不同特征的相互作用区域的颜色差异越大。
这个例子的通过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节所示的这个体系的图一样。
这个例子的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节给出的此体系的图像。
这个例子的动力学模拟是用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节给出的此体系的图像。
本文充分介绍了笔者在Struct. Bond., 190, 297 (2026)中提出的amIGM分析方法,并通过四个例子演示了amIGM在Multiwfn中的分析操作。由本文可见,amIGM分析可以非常生动直观地展现动力学过程中的分子间的平均的弱相互作用,而且分析步骤操作简单,因此非常推荐大家将之应用于实际问题的研究中。amIGM的图像效果远好于aNCI,而耗时相仿佛,因此aNCI就不需要再考虑了。
]]>文/Sobereva@北京科音 2025-Dec-7
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里。
在演示使用严格的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 //同时显示空穴和电子的等值面
至少截止到笔者撰文时(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%
文/Sobereva@北京科音 2025-Nov-24
《使用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,既足够安全,节约耗时的效果又足够显著。
在《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的设置不影响其耗时和结果。
]]>文/Sobereva@北京科音 2025-Nov-24
弱相互作用的图形化分析方法日趋流行,其中笔者提出的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,内容和正式版没明显区别。引用时请引用正式版
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/免费下载。
为了做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号原子了。
碲晶体是由一条条碲原子链有序堆积构成的。作为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,因此暗示了很轻微的共价作用。
本文介绍并示例了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中计算的例子,非常欢迎阅读!]]>
文/Sobereva@北京科音 2025-Nov-19
很多人在对周期性体系绘制三维函数(如电子密度差)的等值面时,容易遇到一个情况是发现感兴趣的等值面在晶胞的边缘,导致等值面被截断,无法观看完整。虽然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)。
北京科音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观看等值面,现在看到的图如下。可见等值面已完全在盒子中央了,非常容易考察水合电子的形态。
此功能对任何函数的格点数据都可以用。格点数据可以是从cube、CHGCAR等Multiwfn支持的记录格点数据的文件中读取的,也可以是Multiwfn的主功能5等功能基于波函数文件直接计算出来并存在内存中的。
本文介绍的功能对于盒子是非正交的情况也可以照常用,平移的方向对应于实际三个轴的方向,只不过Multiwfn目前无法正确显示这种情况的等值面,应当用VMD、VESTA等程序显示。
前面的例子中,格点数据对应的盒子(即均匀分布的格点所处的范围)和做周期性第一性原理计算对应的盒子(对晶体体系来说也相当于晶胞)是完全一致的。如果格点数据的盒子范围和周期性计算用的盒子范围不对应,比如计算格点数据的区域只是整个体系中的一小块,那么使用本文的功能就没任何意义。
如果你的体系本来就没周期性,本文的功能虽然也能使用,但并没有任何实际意义。
如果本文的功能给你的研究带来了便利,发表文章时请引用Multiwfn启动时提示的原文。
]]>文/Sobereva@北京科音 2025-Oct-31
碳原子形成环状的体系如今受到越来越多的关注,自2019年合成的18碳环(cyclo[18]carbon)开始,目前已经有越来越多的碳单环体系被合成,小到C6,大到C48、C50碳环的合成都已经被报道,碳环的复合物、衍生物在未来极可能成为有重要实际应用价值的一类体系。北京科音自然科学研究中心(http://www.keinsci.com)的卢天和江苏科技大学的刘泽玉等人近年来对碳环类体系做了大量的理论研究工作,汇总见//www.umsyar.com/carbon_ring.html,近期发表的综述见Acc. Mater. Res., 6, 1220 (2025) https://doi.org/10.1021/accountsmr.5c00131。
超碱原子是具有M3O化学组成的三角形的簇,M是碱金属原子。M3O的电离能比碱金属原子还要更小,早已被用于设计具有显著非线性光学性质的体系,典型的如卢天、丁迅雷等人在J. Comput. Chem., 38, 1574 (2017)中构造的M3O@Si12C12型体系。将M3O和碳环相结合能够获得具有何种特征的体系,无疑十分值得探索。最近卢天和刘泽玉等人共同理论研究了C2n(n=5-10)碳环与超碱原子M3O(M=Li,Na,K)形成的复合物,十分全面系统地考察了这些复合物的几何结构、电子结构、分子间相互作用、光学吸收、非线型光学性质。研究成果已经发表在英国皇家化学会RSC旗下的J. Mater. Chem. C期刊上,欢迎阅读和引用:
Wenwen Zhao, Jiaojiao Wang, Xiufen Yan, Tian Lu,* Zeyu Liu,* Obtaining excellent optical molecules by screening superalkali-doped cyclo[2n]carbons, M3O@C2n (M = Li, Na, and K, n = 5–10), J. Mater. Chem. C , 13, 17862 (2025) http://doi.org/10.1039/d5tc01675d
此研究表明M3O与碳环可以组成静电吸引作用主导的盐型复合物[M3O]+@[C2n]-,体系的极化率随着原子序数的增加、碳环尺寸的增加而增大,并且具有显著的各向异性。而第一超极化率的变化不具有很强的规律性,筛选发现其中Li3O@C20具有显著的第一超极化率,且只在>200 nm的近紫外和可见光范围具有吸收。本研究对基于碳环构造具有特殊非线型光学特征的体系给予了新的启示,充分体现了将超碱原子置入碳环或其衍生物是一种可行的设计非线性光学材料的思路!
具体内容请阅读原文,以下仅对文中的部分图片进行展示。
不同超碱原子M3O与不同碳环形成的复合物的势能面极小点结构:
//www.umsyar.com/621介绍的IGMH方法展现的M3O与C20碳环之间的相互作用特征、不同M3O@C2n复合物的相互作用能,以及//www.umsyar.com/685介绍的sobEDAw方法做的能量分解
不同M3O@C2n复合物的各向同性极化率,以及按照//www.umsyar.com/547介绍的方法图形化展现的Li3O@C20的第一超极化率张量
M3O@C20的极化率和第一超极化率随碱金属原子的变化,以及极化率和第一超极化率随外场频率的变化
Li3O@C20的电子光谱、可见光范围的关键激发态S0-S8的轨道跃迁特征,以及按照//www.umsyar.com/434介绍的方法展现的空穴和电子分布。
文/Sobereva@北京科音 2025-Sep-20
Dalton程序有不可替代的价值,例如北京科音高级 培训班(http://www.keinsci.com/KAQC)里笔者会讲到使用Dalton算磷光速率、旋轨耦合矩阵元(包括激发态间的)、双光子吸收,这都是Dalton的特色功能。网上流传着一些错误显著错误的关于Dalton的溶剂模型的说法,甚至有的文中说什么“必须自定义空腔半径,否则计算无法继续下去”,明显误人子弟。本文专门说说用户应该了解的Dalton的PCM隐式溶剂模型的几个关键知识,其中不少内容在手册里都找不到。本文的情况对应Dalton 2018-2022版,对其它版本可能适用也可能不适用。
Dalton的隐式溶剂模型只支持PCM,不支持SMD等其它的。以下是个Dalton做B3LYP级别的单点+偶极矩计算的最简单的.dal文件例子,使用PCM模型描述水环境。
**DALTON INPUT
.RUN PROPERTIES
*PCM
.SOLVNT
WATER
*PCMCAV
**WAVE FUNCTIONS
.DFT
B3LYPg
**END OF DALTON INPUT
.SOLVNT下面是溶剂名,自带18种溶剂,包括(写化学式还是溶剂名都可以):
H2O WATER
CH3OH METHANOL
C2H5OH ETHANOL
CHCL3 CLFORM
CH2CL2 METHYLCL
C2H4CL2 12DICLET
CCL4 TETRACLC
C6H6 BENZENE
C6H5CH3 TOLUENE
C6H5CL CLBENZ
CH3NO2 NITROMET
C7H16 N-EPTANE
C6H12 CYCLOHEX
C6H5NH2 ANILINE
CH3COCH3 ACETONE
THF TETHYDFU
DMSO DIMETSOX
CH3CN ACETONIT
输出文件中可以看到当前的溶剂信息
** LOOKING UP INTERNALLY STORED DATA FOR SOLVENT = WATER **
Optical and physical constants:
EPS= 78.390; EPSINF= 1.776; RSOLV= 1.385 A; VMOL= 18.070 ML/MOL;
*PCMCAV字段用于定义溶质孔洞,必须出现,如果下面没做额外设置,代表以默认的方式产生溶质的孔洞。Dalton用的是 程序中普遍使用的靠原子范德华球叠加的方式构造孔洞。
*PCM里加入.NEQRSP可以使用非平衡溶剂,计算对应吸收的电子激发问题用得着。
如果要自定义新溶剂,比如静态介电常数为23.0、动态介电常数为1.84的溶剂,*PCM部分改为:
*PCM
.SOLVNT
WATER
.EPS
23
.EPSINF
1.84
Dalton带着PCM计算时,输出信息里可以看到诸如以下内容,给出的是各个原子的以Bohr为单位的X、Y、Z坐标和以埃为单位的构造溶质孔洞用的原子半径。
********SPHERES IN PCMSPHGEN************
INDEX X Y Z R
1 2.1499166192D+00 -1.9894550417D+00 1.6741801837D+00 1.2000000000D+00
2 4.0096560598D+00 2.5073910458D-01 0.0000000000D+00 1.2000000000D+00
3 2.1499166192D+00 -1.9894550417D+00 -1.6741801837D+00 1.2000000000D+00
...略
这些原子半径是怎么来的?实际上在源代码目录下的sirius/sircav.F中可以发现以下内容,可见用的是Bondi原子半径,但个C、N、O的半径是U. Pisa修改后的。对没定义半径的元素,半径直接当成0。
! A.Bondi, J.Phys.Chem. 68: 441-451(1964) gives alternate
! values, and a few transition metals.
allocate(rvdw(99))
rvdw = (/ 1.20d0, 1.22d0, 0.00d0, 0.00d0, 2.08d0, 1.85d0,
& 1.54d0, 1.40d0, 1.35d0, 1.60d0, 2.31d0, 0.00d0,
& 2.05d0, 2.00d0, 1.90d0, 1.85d0, 1.81d0, 1.91d0,
& 2.31d0, 13*0.0d0, 2.00d0, 2.00d0, 1.95d0, 1.98d0,
& 2.44d0, 13*0.0d0, 2.20d0, 2.20d0, 2.15d0, 0.00d0,
& 2.62d0, 27*0.0d0, 2.40d0, 16*0.0d0 /)
! override the above table with U. Pisa's experience
! as to what works best for singly bonded C,N,O
rvdw(6) = 1.70d0
rvdw(7) = 1.60d0
rvdw(8) = 1.50d0
所以,Dalton用的原子半径不是胡来的,是有依据的,不存在无条件的“必须自定义空腔半径”。仅当你的体系牵扯到部分元素没有内置半径时(而且那些元素的原子与溶剂有所接触时)才绝对需要补充半径。
本节给一个例子,B3LYP/def2-SVP计算水中的乙醇,通过在Dalton中自定义原子半径使得其计算结果与Gaussian 16默认的IEFPCM的相吻合,从而加深读者对Dalton的PCM模型的理解和信心。
以下是Gaussian的输入文件,用scrf关键词默认的IEFPCM溶剂模型表现默认的水环境
#P B3LYP/def2SVP nosymm scrf
[空行]
Generated by Multiwfn
[空行]
0 1
C 1.17229118 -0.41192328 0.00000000
H 1.13768688 -1.05277427 0.88593800
H 2.12181861 0.13268542 0.00000000
H 1.13768688 -1.05277427 -0.88593800
C -0.00000000 0.55479430 0.00000000
H 0.05413938 1.20787442 -0.88657584
H 0.05413938 1.20787442 0.88657584
O -1.19922622 -0.21255820 0.00000000
H -1.94540849 0.40035379 0.00000000
输出文件在此://www.umsyar.com/attach/752/G16_scrf.out。从其中以下内容里可以看到溶剂信息和原子半径信息。其中Re0是原子半径(埃),Gaussian 16对IEFPCM模型默认用的是UFF力场定义的范德华半径。Alpha是给各个原子半径乘的系数。
Solvent : Water, Eps= 78.355300 Eps(inf)= 1.777849
------------------------------------------------------------------------------
Spheres list:
ISph on Nord Re0 Alpha Xe Ye Ze
1 C 1 1.9255 1.100 1.172291 -0.411923 0.000000
2 H 2 1.4430 1.100 1.137687 -1.052774 0.885938
3 H 3 1.4430 1.100 2.121819 0.132685 0.000000
4 H 4 1.4430 1.100 1.137687 -1.052774 -0.885938
5 C 5 1.9255 1.100 0.000000 0.554794 0.000000
6 H 6 1.4430 1.100 0.054139 1.207874 -0.886576
7 H 7 1.4430 1.100 0.054139 1.207874 0.886576
8 O 8 1.7500 1.100 -1.199226 -0.212558 0.000000
9 H 9 1.4430 1.100 -1.945408 0.400354 0.000000
------------------------------------------------------------------------------
当前Gaussian 16计算的偶极矩为1.8592 D。
下面是Dalton的.dal文件。.ICESPH设2代表要自定义原子半径,.NESFP设9代表要读入9个原子的半径(乙醇有9个原子,此例全都自定义)。.INA下面写上各个要修改半径的原子的序号。.RIN下面是各个原子的半径(埃),这里改成和Gaussian用的相同的UFF半径。注意顺序要按照.mol文件里的写,本例的.mol见//www.umsyar.com/attach/752/ethanol.mol。.ALPHA定义各个原子的alpha值,都设成和Gaussian相同的1.1。Dalton里的B3LYP关键词对应的并不是Gaussian里的B3LYP,因此这里写B3LYPg使用和Gaussian相同的定义。
**DALTON INPUT
.RUN PROPERTIES
*PCM
.SOLVNT
WATER
.ICESPH
2
.NESFP
9
*PCMCAV
.INA
1
2
3
4
5
6
7
8
9
.RIN
1.4430
1.4430
1.4430
1.4430
1.4430
1.4430
1.9255
1.9255
1.7500
.ALPHA
1.1
1.1
1.1
1.1
1.1
1.1
1.1
1.1
1.1
**WAVE FUNCTIONS
.DFT
B3LYPg
**END OF DALTON INPUT
Dalton的输出文件是//www.umsyar.com/attach/752/PCMmodcav_ethanol.out,可见其中的偶极矩的计算结果为1.8599 D,和Gaussian 16给出的1.8592完美吻合!
再顺带一提怎么让Gaussian得到Dalton默认孔洞设置下的结果。Dalton对H、C、O设的半径分别是1.2、1.7、1.5埃,alpha默认用的是1.2,因此Gaussian输入文件要写成下面这样与之对应:
#P B3LYP/def2SVP nosymm scrf=read
[空行]
Generated by Multiwfn
[空行]
0 1
C 1.17229118 -0.41192328 0.00000000
H 1.13768688 -1.05277427 0.88593800
H 2.12181861 0.13268542 0.00000000
H 1.13768688 -1.05277427 -0.88593800
C -0.00000000 0.55479430 0.00000000
H 0.05413938 1.20787442 -0.88657584
H 0.05413938 1.20787442 0.88657584
O -1.19922622 -0.21255820 0.00000000
H -1.94540849 0.40035379 0.00000000
[空行]
alpha=1.2
modifysph
[空行]
H 1.2
C 1.7
O 1.5
[空行]
[空行]
Gaussian 16的偶极矩计算结果为1.9136 D,Dalton在默认孔洞下计算的输出文件为//www.umsyar.com/attach/752/PCM_ethanol.out,结果为1.9139 D,可见再次吻合得极好!
直接按照前面说的方式在Dalton中用PCM模型做计算,最终会出现警告
576: Warning, element 34 263 of SI too big: set to zero
577: Warning, element 50 249 of SI too big: set to zero
578: Warning, element 52 249 of SI too big: set to zero
579: Warning, element 56 247 of SI too big: set to zero
...
很多人看到这警告就大惊失色了。实际上这些Warning可以直接无视,如果想不显示这些Warning,在*PCM里加上以下内容即可(在Dalton自带的涉及PCM的测试文件里总是带着这个)
.NPCMMT
0
这里说一下原因。.NPCMMT选项在手册里找不到,但在源代码sirius/sirpcm.F里能看到其说明
NPCMMT = 0 No correction of the DI, SI and C matrices
NPCMMT = 1 Correction of DI and SI (default)
NPCMMT = 2 Correction of DI, SI and C
进而看sirief.F,会发现NPCMMT为1、2的时候都会自动做矩阵的检查,并可能导致那些warning的出现。是否会出现warning和是否自定义原子半径并没必然关系。NPCMMT为1、2时做的校正其实没什么意义,和为0时的结果差异微乎其微。对于前例,为1时偶极矩为1.913867、单点能为-154.9289657907 Ha,为0时分别为1.913182和-154.9289638876 Ha。
计算时还有可能出现** WARNING ** A VERY POOR TESSELATION HAS BEEN CHOSEN,这通过设NPCMMT为0也不会避免。根据我的测试,这个warning一般也是无害的,结果还是正常的。实在不放心的话可以按照上文的方式,把Gaussian中的原子半径设成与Dalton相同后对Dalton的结果做验证,也可以尝试把Dalton的原子半径设成和Gaussian默认的一致看看是否warning能消除,且与Gaussian结果相符。
]]>
文/Sobereva@北京科音 2025-Aug-31
D-pi-A型分子普遍具有较大的第一超极化率,其中D(donor)是电子给体,A(acceptor)是电子受体,pi是指具有全局pi共轭特征的片段,可以称为linker或者bridge。2019年首次通过STM观测到的由18个碳原子构成的环状体系,18碳环(cyclo[18]carbon),同时具有平面内和平面外两套全局pi共轭特征,这点在《18碳环及衍生物的十分全面系统的研究综述已在Acc. Mater. Res.期刊发表!》(//www.umsyar.com/749)介绍的综述中有充分说明,如果你对18碳环缺乏了解十分推荐阅读此文。显然,18碳环的这种特征使得它也很有希望作为pi-linker构建具有显著第一超极化率的独特的D-pi-A型分子。基于这个idea,北京科音自然科学研究中心的卢天和江苏科技大学的刘泽玉等人对18碳环与不同donor、acceptor相连接构成的D-pi-A型分子的特征,尤其是光学性质,做了充分全面的理论研究。此工作近期已经发表在Phys. Chem. Chem. Phys.期刊上,欢迎阅读和引用:
Jingbo Xu, Jiaojiao Wang, Xiaohui Chen, Wenwen Zhao, Xiufen Yan, Zeyu Liu,* Tian Lu,* Aihua Yuan,* Design of donor–π–acceptor type cyclo[18]carbon derivatives for infrared nonlinear optical materials: a theoretical perspective, Phys. Chem. Chem. Phys., 27, 11993 (2025) https://doi.org/10.1039/d5cp00736d
下文将会对这篇研究文章的主要内容进行介绍并添加一些附加信息,使读者更好地理解文章的研究思想和主要结果,图片、表格来自于文章正文或补充材料。同作者之前发表过大量其它的和18碳环及其衍生物、复合物有关的研究工作,十分推荐阅览//www.umsyar.com/carbon_ring.html里的汇总,其中有很多篇研究工作都与18碳环及相关体系的光学性质有关。
此文研究的将18碳环作为pi-linker构建的示意图如下所示,共考虑四种情况:(1)两侧都连接氢 (2)一侧连氢一侧连硝基 (3)一侧连氨基一侧连氢 (4)一侧连氨基一侧连硝基。其中(4)对应最典型的D-pi-A形式的体系。而连氢的话,在电子激发时氢不会起到明显的给/吸电子效应,它的作用主要是让相连的碳从原本在碳环中的sp杂化状态变成sp2杂化。
此文用Carbon, 165, 468 (2020)论证过描述碳环十分可靠的ωB97XD/def2-TZVP级别对上述体系做了几何优化,得到的无虚频的笛卡尔坐标在文中的补充材料里都提供了。各体系的碳环部分的长、宽(上图的d1和d2)在下表给出了,可见无论碳环连的是什么基团,环的形状都没什么区别,都是椭圆形。
各个体系的结构图如下(补充材料图S2),键长以埃为单位标注了,并根据键长对键进行了着色。由图可见,这些18碳环衍生物的C-C键长度是交替变化的,这个特征和孤立状态的碳环一致。
前面的表格里QM和Qn分别是体系的各个部分的Mulliken和NPA电荷,原子电荷的相关知识见《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表!》(//www.umsyar.com/714)介绍的笔者写的综述。可见无论哪种原子电荷计算方法,H-C18-H和NH2-C18-H中,18碳环都能获得电子,故带有明显的负电。而在H-C18-NO2和NH2-C18-NO2中,由于吸电子能力更强的-NO2从18碳环上吸了一些电子,导致18碳环部分的净电荷恰好接近0。
前面表格里的μ0是体系的基态的永久偶极矩。有趣的是H-C18-H的偶极矩并不为0,因为此体系并不是精确中心对称的(哪怕用很严收敛限优化也是如此),仔细看前面图中标注的键长可以看出这点来。H-C18-NO2和NH2-C18-H的偶极矩都较大,而D-pi-A特征最鲜明的NH2-C18-NO2具有特别大的偶极矩,高达15.5 Debye。
按照《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图》(//www.umsyar.com/443)绘制的4个18碳环衍生物的范德华表面静电势填色图如下所示,由此可以清楚直观地看出它们的电荷分布特征的差异。颜色越红、越蓝分别体现出其附近的原子带正电、负电越显著。并且表面静电势最大点和最小点的位置和数值也都标注了出来。可见虽然NH2-C18-NO2的碳环中部并不怎么带净电荷,但左右两端由于连接的基团的给/吸电子特征存在巨大差异,因而形成了显著的正负电荷分离现象,这也是为什么它有最大的偶极矩。
下图是基于《使用Multiwfn绘制态密度(DOS)图考察电子结构》(//www.umsyar.com/482)介绍的方法绘制的18碳环衍生物的TDOS和PDOS图,HOMO、LUMO位置以及相应的gap都标注了。在相同计算级别下,单独的18碳环的HOMO和LUMO能级分别为-8.45和-1.70 eV,对应HOMO-LUMO gap为6.75 eV,而当前研究的18碳环衍生物的gap在4.62-4.9 eV范围。可见接上基团,令两端的碳从sp杂化变成sp2杂化从而一定程度破坏平面内的芳香性会使得gap有所缩窄。PDOS曲线体现出HOMO、LUMO轨道几乎都是C18部分贡献的,但氨基,特别是硝基,对于其它前线轨道也有明显贡献,这暗示出一些低阶电子激发态会牵扯到它们。
本研究使用TDDFT在ωB97XD/def2-TZVP级别下做了电子激发计算,算了最低50个态,并模拟了电子光谱,如下所示。可见当前研究的18碳环衍生物的吸收分为三个部分:500多nm的可见光区的弱吸收、300多nm的近紫外区的弱吸收、200nm左右的更远紫外区的强吸收。如Carbon, 165, 461 (2020)的计算研究所示,独立的18碳环是没有可见光吸收的,而当前研究表面给它恰当接上H、氨基、硝基后能够令体系显色。基团的种类一定程度影响最大峰位置和峰高度,但影响程度很有限。
为了更好地理解这些体系的吸收光谱的内在本质,文中按照《使用Multiwfn绘制电荷转移光谱(CTS)直观分析电子光谱内在特征》(//www.umsyar.com/628)介绍的方法绘制了电荷转移光谱,如下所示。可见不管接上什么基团,>300 nm区域的吸收几乎完全都是18碳环片段上的局域激发,而氨基和硝基则对于200nm左右的吸收有少量贡献,使得对应的电子激发有一些电荷转移特征,但仍以18碳环区域的激发为主。
能量最低同时振子强度又显著的电子激发态被称为crucial state(关键态),对于当前研究的体系来说对应于可见光区的吸收。各个体系的关键态的信息如下所示,包括激发能、振子强度、激发态偶极矩相对于基态偶极矩的变化的大小
为了更好地了解这些激发的本质特征,文中按照《使用Multiwfn做空穴-电子分析全面考察电子激发特征》(//www.umsyar.com/434)介绍的方法对各个体系的关键态绘制了空穴(蓝色等值面)、电子(绿色等值面)分布图和热图,如下所示,图中还标注了两种衡量电子激发特征的定量指标,在434博文里也都详细介绍了。可见尽管这些激发主要都来自于碳环部分的pi电子,而结合定量指标D和Sr,还是能看出差别的。四个体系中,NH2-C18-NO2的关键态具有最显著的电荷转移特征,而且从蓝色和绿色等值面上可以看到是NH2向NO2方向整体一定程度转移电子。而H-C18-H的电荷转移特征几乎没有,空穴和电子几乎都是相对于体系中心对称分布的。H-C18-NO2的情况较接近H-C18-H,而NH2-C18-H的情况较接近NH2-C18-NO2。
电子激发态特征和第一超极化率之间有密切联系,通过SOS公式关联起来,在《使用Multiwfn基于完全态求和(SOS)方法计算极化率和超极化率》(//www.umsyar.com/232)里有详细介绍。SOS简化后可以得到双、三能级公式,见《使用Multiwfn对第一超极化率做双能级和三能级模型分析》(//www.umsyar.com/512)和《谈谈计算第一超极化率的双能级公式》(//www.umsyar.com/361)。双能级公式用来讨论影响第一超极化率的本质因素十分有用且流行。双能级公式指出静态第一超极化率与关键态的振子强度呈正比、与激发态相对于基态的偶极矩变化的大小呈正比、与激发能的三次方呈反比。根据前面表格里的关键态信息可见,相对于其它三个18碳环衍生物,NH2-C18-NO2的关键态具有最大的振子强度、最大的偶极矩变化,而激发能和它们则差不多,根据双能级公式可以预期NH2-C18-NO2有最大的第一超极化率。但双能级公式只适合用来解释和定性预测第一超极化率,更严格的结论、准确的数值还是需要依赖于用严格的方法定量计算,见下一节。
此文使用Gaussian 16在ωB97XD/aug-cc-pVTZ(-f,-d)级别下以CPKS方法解析地计算了前述四种18碳环衍生物的二次谐波生成(SHG)形式的静态和动态的第一超极化率,β(-2ω;ω,ω)。aug-cc-pVTZ(-f,-d)基组是常用的aug-cc-pVTZ基组去掉非氢原子的f和氢原子的d极化函数,根据笔者以前的充分的测试,这种简化可以节约很多计算时间而并不会令计算精度下降太多。文中具体考察的量包括第一超极化率的总大小(βtot)、平行于偶极矩的分量(βvec)、各笛卡尔分量(βX、βY、βZ),这些量都可以通过Multiwfn基于Gaussian输出文件直接产生,见《使用Multiwfn分析Gaussian的极化率、超极化率的输出》(//www.umsyar.com/231)。
文中计算出的静态(∞ nm)和1460、1907 nm下的动态第一超极化率的总大小如下图所示。可见H-C18-H的β相当小,光是把一个H替换成硝基也没什么变化,但替换成氨基则能令β显著增大,而同时连上氨基和硝基构成典型D-pi-A体系后β变得颇大,这体现出氨基和硝基在提升电场响应性质方面的显著的协同作用。从下图还可以看到随着外场波长的减小,即外场频率的增大,对应的β有大幅提升,即这些体系具有很显著的频率-色散效应。
NH2-C18-NO2的静态的βtot为4462 a.u.,算是颇大了,同样是D-pi-A型分子的p-nitroaniline (NH2-C6H4-NO2)的实验气相值只有1400 a.u.。可见18碳环作为pi-linker有独特价值,是个颇好的pi-linker。
文中还使用《使用Multiwfn通过单位球面表示法图形化考察(超)极化率张量》(//www.umsyar.com/547)介绍的方法直观考察了各个体系的第一超极化率的各向异性,并基于《使用Multiwfn计算(超)极化率密度》(//www.umsyar.com/305)介绍的方法考察了不同区域、片段对第一超极化率的贡献,请感兴趣的读者阅读原文中的讨论。
为了考察溶剂效应对光学性质的影响,文中还在IEFPCM隐式溶剂模型描述的环己烷环境下也做了计算,发现环己烷会导致光谱略微红移、振子强度略微下降,但整体特征不变。环己烷还会造成绝大多数情况β的显著增加,但不同体系的β大小的趋势并未改变。
本文介绍的Phys. Chem. Chem. Phys., 27, 11993 (2025)文章首次对18碳环与H、氨基、硝基构成的分子的几何结构、电子结构和光学性质做了充分的理论研究,研究证实18碳环可以作为很理想的pi-linker构造典型的D-pi-A体系以获得具有显著非线性光学性质的分子,文中构建的NH2-C18-NO2的第一超极化率已达到了很大的值。此文的研究工作显著拓宽了化学家们对碳环的衍生物的认识,对于未来利用碳环作为骨架构造非线性光学材料提供了重要的参考。值得一提的是,目前从小到大很多碳环皆已实验合成,以不同尺寸的碳环作为pi-linker,有望获得具有不同非线性光学性能的分子。
]]>The way of plotting electrostatic potential for local region of molecular surface using Multiwfn in combination with VMD