计算化学公社
标题: 使用IRI方法图形化考察化学体系中的化学键和弱相互作用 [打印本页]
作者Author: sobereva 时间: 2021-6-1 11:09
标题: 使用IRI方法图形化考察化学体系中的化学键和弱相互作用
使用IRI方法图形化考察化学体系中的化学键和弱相互作用
Using IRI method to graphically study chemical bonds and weak interactions in chemical systems
文/Sobereva@北京科音 First release: 2021-May-31 Last update: 2024-Mar-29
1 前言
原子间的相互作用在化学体系中无处不在。把相互作用图形化展现无疑对考察化学问题非常有益，可以令化学家快速直观了解体系中什么位置有什么样的相互作用。目前图形化展现相互作用的方法不少，有很多都已经非常流行，在波函数分析程序Multiwfn（http://sobereva.com/multiwfn）里全都支持，这里先简单回顾一下。Multiwfn中支持的RDG方法（也称NCI方法）可以通过着色等值面图直观地展现弱相互作用，见这些博文以及其中引用的其它相关博文：《使用Multiwfn图形化研究弱相互作用》（http://sobereva.com/68）、《使用Multiwfn结合CP2K通过NCI和IGM方法图形化考察固体和表面的弱相互作用》（http://sobereva.com/588）、《使用Multiwfn做aNCI分析图形化考察动态过程中的蛋白-配体间的相互作用》（http://sobereva.com/59）。Multiwfn支持的Hirshfeld surface分析可以通过着色的表面展现分子间相互作用位置和强度，见Multiwfn手册4.12.5和4.12.6节的例子以及演示视频《使用Multiwfn结合VMD绘制Hirshfeld surface图分析分子晶体中的弱相互作用》（https://www.bilibili.com/video/av35738671/）。也有一些方法可以图形化展现共价作用，例如Multiwfn支持的ELF、LOL、SCI、价层电子密度、电子密度拉普拉斯函数、变形密度等，相关博文和文章汇总见《ELF综述和重要文献小合集（http://bbs.keinsci.com/thread-2100-1-1.html）、《使用Multiwfn作电子密度差图》（http://sobereva.com/113）、《通过价层电子密度分析展现分子电子结构》（http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB20170925）。为了将弱相互作用和共价作用共同展现，有人提出将RDG和ELF联用，见比如此文的例子：《通过键级曲线和ELF/LOL/RDG等值面动画研究化学反应过程》（http://sobereva.com/200）。此外，Multiwfn支持的DORI函数可以通过其等值面直接将弱相互作用和化学键作用一起展现，见《使用DORI函数同时考察共价和非共价相互作用》（http://sobereva.com/367）。Multiwfn支持的IGMH方法也可以将化学键和弱相互作用都予以展现，但对不同强度的相互作用需要用不同等值面数值才能展现得比较鲜明，见《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（http://sobereva.com/621）。另外Multiwfn还支持其它一些图形化的分析手段，如流行的atoms-in-molecules (AIM)、定域化轨道、键级密度（BOD）等，见《Multiwfn支持的分析化学键的方法一览》（http://sobereva.com/471）和《Multiwfn支持的弱相互作用的分析方法概览》（http://sobereva.com/252）中的相关信息。
虽然如上所示，已有的图形化展现相互作用的方法已经较为丰富，但之前始终没有一个能够很好地同时展现各种类型相互作用的方法。虽然DORI看似能做到，但其函数不仅定义很复杂（因此计算耗时也高、编程实现也不便），图像效果还特别差，等值面看着很杂乱。IGMH虽然也能展现各类相互作用，但没法在一张图里（即同一个等值面数值下）同时把不同强度的相互作用都展现清楚。RDG与ELF联用虽然图像效果还不错，但由于需要同时考虑两个函数，不仅绘图、分析起来很麻烦，计算耗时也高。
笔者通过对RDG函数的特征进行分析，发现只需要稍加修改，就可以通过其等值面非常清楚地把所有类型相互作用在一幅图里都展现出来，这个新的函数笔者命名为interaction region indicator (IRI)，可以译为“相互作用区域指示函数”。包含此函数的原理介绍、与其它方法的分析对比以及诸多应用实例已于近期发表在专注于化学领域方法学发展的新刊Chemistry-Methods上，欢迎阅读（此期刊是开放访问的，可免费下载），也请记得引用：
Tian Lu,* Qinxue Chen, Interaction Region Indicator (IRI): A Simple Real Space Function Clearly Revealing Both Chemical Bonds and Weak Interactions, Chemistry—Methods, 1, 231-239 (2021) https://chemistry-europe.onlinelibrary.wiley.com/doi/10.1002/cmtd.202100007
也非常推荐阅读和引用此文介绍的笔者的综述文章：《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（http://sobereva.com/667），其中对IRI有充分的介绍并给了诸多应用例子。
在下面第2节，笔者将对IRI函数进行简要的介绍，更多细节和讨论请看上面的论文。然后在第3、4节将分别通过分子体系和周期性体系的实例演示如何通过Multiwfn波函数分析程序与VMD可视化程序相结合绘制IRI图。之后在第5节，将介绍笔者在上面论文里连带提出的IRI的变体IRI-π，此函数可以直观地展现化学体系中的π相互作用位置、类型和强度。由于IRI和IRI-pi计算方便快捷，有免费易用的计算程序，图像效果好，普适性强，预期在未来会得到广泛应用。
后文计算实例中用到的文件都可以在本文的文件包http://sobereva.com/attach/598/file.zip里找到。
2 IRI方法简介
2.1 IRI的定义
IRI函数的定义特别简单，如下所示：
(, 下载次数 Times of downloads: 133)
可见，只要有电子密度ρ及其梯度即可计算IRI。做完常规量子化学计算后，就有了波函数信息，就可以很容易地计算电子密度及其梯度。由于IRI不依赖于波函数，因此如果有很高解析度的实验测定的电子密度分布数据的话也可以用来计算IRI。
IRI中的参数a非常关键。IRI标准定义中的a取1.1，这是笔者对许多体系进行测试经验选取的最佳数值，此时IRI可以最理想地同时展现化学键和弱相互作用。如果a=4/3，那么IRI和RDG就仅相差前面一个常数系数而已。因此对于已经支持了RDG分析的程序，程序开发者仅仅需要改一行代码即可实现IRI。
和RDG分析一样，也可以把sign(λ2)ρ函数通过不同颜色投影到IRI等值面上来区分不同区域的作用强度和特征。这里λ2代表电子密度Hessian矩阵第二大本征值，sign()代表取符号，sign(λ2)ρ的含义请看IRI方法的原文，或者笔者对RDG方法的介绍文章《使用Multiwfn图形化研究弱相互作用》（http://sobereva.com/68）。在IRI原文中笔者建议使用下面的色彩刻度对IRI等值面通过sign(λ2)ρ来着色，这个图是Multiwfn程序的examples目录下的IRI_colorbar.png，用IRI图发文章的时候可以直接用，其中单位是a.u.
(, 下载次数 Times of downloads: 196)
上面这种着色方式意味着如果某IRI等值面基本上是绿色，就说明此处是范德华作用区域（也可以是比如极弱的氢键，这种程度的氢键是色散作用占主导，详见《透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍》http://sobereva.com/513）。如果等值面颜色明显偏红，说明这里存在一定位阻作用，若是鲜红则说明位阻很强。如果等值面颜色明显偏蓝，说明存在显著的吸引作用，比如一般强度的氢键、卤键等。如果等值面完全是蓝色的，说明此处要么是相对来说很强的弱相互作用，因此作用区域的电子密度能达到>=0.04 a.u.的程度，要么是化学键作用，成键区域电子密度通常显著大于0.04 a.u.。
2.2 IRI与其它方法的对比
对于展现弱相互作用区域，IRI和RDG图像其实没有显著区别，但如果也要同时展现化学键作用，那么IRI相对于RDG的图像效果方面的优势很显著。下图是一个非常典型的含有分子间、分子内弱相互作用以及化学键作用的苯酚二聚体。
(, 下载次数 Times of downloads: 127)
如上可见，作RDG图并且不设密度截断的话（平时做RDG分析弱相互作用时通常都设密度截断，从而令密度较大的地方不显示RDG等值面），如下所示，当等值面数值选取得适合展现分子间弱相互作用时（RDG=0.5），化学键作用就展现不清楚，等值面都连成一片；而当等值面数值设为RDG=0.25，此时化学键作用展现得不错，但弱相互作用区域的等值面就非常小了，很难观察。而IRI则选取了一个十分恰当的a参数，在IRI=1.1等值面下同时把共价键、环内位阻区域、分子间氢键和范德华作用区域都清晰地展现了出来。根据笔者经验，对多数体系IRI展现各种相互作用比较合适的等值面数值是0.9~1.1，Multiwfn提供的绘制IRI图的脚本里用的是1.0，对绝大多数情况效果都令人满意。大家也可以自行根据实际情况微调等值面数值来令等值面大小更合适、图像更美观、更易于分析。
IRI和DORI函数的目的一致，都是同时展现化学键和弱相互作用区域。但相比于IRI，DORI的图像效果实在太差了，计算还更耗时、定义还更复杂，所以以后DORI就没有任何存在意义了。下面是一个对比，考察的是实验上已合成出来的两个金刚烷之间由C-C键连接的体系。可见IRI图把范德华作用、位阻作用和化学键都展现得清清楚楚，图像很干净；而DORI图则显得噪音很明显，有许多细碎的等值面，而且描述化学键的等值面还都没有分离开，描述两个金刚烷之间作用的等值面边缘还特粗糙。DORI图像效果对等值面数值的选取特别敏感，但无论怎么调，图像效果始终都很烂。
(, 下载次数 Times of downloads: 121)
再看下面的G-C碱基对体系。(a)、(b)分别是sign(λ2)ρ着色的IRI和DORI等值面图，可见前者的等值面干干净净，把所有化学键、碱基间的氢键和范德华作用都展现了出来，而DORI图还是显得很粗糙，等值面形状很难看，许多化学键也区分不开。DORI的等值面效果为啥这么烂，可以从体系平面上的填色图来更充分考察。下图(c)、(d)分别是体系平面上的IRI和DORI填色图，由前者可以看到IRI函数变化比较光滑，从这幅图也可以直接清楚地看出在这个平面上此体系的各种相互作用出现的区域。而从DORI的填色平面图上看DORI的变化非常复杂、凌乱，自然其等值面图效果也很差，这本质上是因为DORI函数本身定义得就很不好。
(, 下载次数 Times of downloads: 126)
AIM理论通常用键临界点（BCP）、键径展现原子间的相互作用。但是这种方式对许多类型的相互作用区域展现得不够全面、直观，特别是有些明显存在氢键相互作用的地方还会被漏掉，即找不到相应的BCP。在IRI的原文中笔者就这个问题做了详细的理论分析，充分证明了为什么有些氢键只能靠IRI或RDG等值面才能展现，而相应的BCP则不存在。文中的讨论使用了下图的Ni(NH3)2(OH)2做为例子，从图上可见N-H...O的氢键没有对应的BCP，而在IRI图上则显示出蓝色的等值面，体现出这是绝对不可忽视的、强度不很弱的内氢键。
(, 下载次数 Times of downloads: 120)
2.3 IRI应用于一些体系的效果
下面是IRI原文里绘制的各种类型体系的着色IRI图，笔者就不再一一分析了，原文里讨论得很详细，请自行阅读。总之从这些图可以看到，不管是普通的共价键、离子键，还是配位键，不管是周期表靠前还是靠后的元素的成键，IRI等值面一律都能将它们展现得非常清楚。与此同时，分子内各种弱相互作用也都能通过IRI图看得非常清楚，和化学直觉完全一致。值得一提的是，对于下面图中的二茂铁，不仅Fe-C作用通过等值面上的蓝色部分展现了出来，笼状区域以及C-Fe-C三元环内区域的强烈位阻作用也都以红色鲜明地体现了出来。
(, 下载次数 Times of downloads: 130)
(, 下载次数 Times of downloads: 123)
与RDG+ELF联用明显不同的是，IRI是单一且随核坐标改变能够平滑变化的函数，因此不仅可以用于考察势能面的极小点结构，还可以用于过渡态，乃至IRC路径上的各个点，从而生动鲜活地展现出整个化学过程中原子间相互作用是如何变化的，这对于深入理解反应的特征和本质非常有好处。下面是IRI原文的例子，是Diels-Alder加成过程，可见最初乙烯和1,3-丁二烯之间是范德华相互作用，随着反应的进行，两个片段间的等值面形状一点点变化，C-C之间的区域一点点变蓝，充分体现出共价键在一点点形成、逐渐变强，到最后变成了普通的C-C键。
(, 下载次数 Times of downloads: 133)
与上图对应的动图如下，可见此动画对反应中相互作用的变化展现得极为生动，也很能激发外行人的兴趣。通过批处理脚本，就可以自动调用Multiwfn对IRC的各个点依次计算、用VMD产生出IRC每个点的图像，最后再用ffmpeg或ImageMagick等程序就可以合并成动画。参考《通过键级曲线和ELF/LOL/RDG等值面动画研究化学反应过程》（http://sobereva.com/200）里的过程并结合Multiwfn里做IRI分析的操作步骤稍微举一反三就能实现。
下面是OH- + CH3CH2Br → CH3CH2OH + Br-的SN2取代反应的IRI的动画，可见也把新键的形成、旧键的断裂、相互作用特征的变化展现得特别直观、灵动。
3 分子体系的IRI分析实例
3.1 苯酚二聚体的等值面图
下面以苯酚二聚体为例展示怎么通过Multiwfn+VMD绘制IRI图。读者务必使用2021年2月及以后更新的Multiwfn版本，Multiwfn最新版本可以在主页http://sobereva.com/multiwfn免费下载。如果对Multiwfn缺乏基本了解的话强烈建议看看《Multiwfn入门tips》（http://sobereva.com/167）和《Multiwfn FAQ》（http://sobereva.com/452）以补充常识知识。本文的VMD使用1.9.3版，可以在http://www.ks.uiuc.edu/Research/vmd/免费下载。Multiwfn+VMD绘制IRI图和绘制RDG图的过程极为相似，如果follow下文的例子遇到困难，建议看看绘制RDG图的演示视频（https://www.bilibili.com/video/av71561024）。能把RDG图画出来自然也就能把IRI图顺利画出来。
做IRI计算需要提供一个含有Multiwfn可以读入波函数信息的文件，用比如wfn/wfx/fch/molden/mwfn等格式都可以，格式介绍以及生成方法见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）。本例就使用Gaussian 16在B3LYP-D3(BJ)/6-311G**级别下对苯酚二聚体做优化任务产生的fch文件当输入文件，这是本文文件包里的PhenolDimer.fchk。
启动Multiwfn，然后依次输入
PhenolDimer.fchk //此处输入PhenolDimer.fchk的实际路径
20 //弱相互作用可视化分析
4 //IRI分析
3 //高质量格点。格点设置显著影响计算耗时和最终得到的图像中等值面的光滑程度，务必看一下《用Multiwfn+VMD做RDG分析时的一些要点和常见问题》（http://sobereva.com/291）中关于格点设置的说明。当前选的“高质量格点”的总格点数是固定的，对中、小体系够用，但对于很大体系时格点间距会较大，等值面还是会有锯齿、不光滑。对很大体系建议手动输入格点间距
用不了多久就算完了，在后处理界面选3 Output cube files to func1.cub and func2.cub in current folder，在当前目录下就出现了func1.cub和func2.cub，分别是记录sign(λ2)ρ和IRI格点数据的cube文件。把这两个文件移动到VMD目录下。然后把Multiwfn目录下的examples子目录下的IRIfill.vmd文件拷到VMD目录下，这是VMD作图脚本。
启动VMD，在文本窗口输入source IRIfill.vmd来执行这个脚本，VMD就会自动载入两个cub文件并且修改作图设置，然后就会看到下图，这就是IRI原文里的图2
(, 下载次数 Times of downloads: 146)
如果想修改IRI等值面数值的话，进入VMD main界面的Graphics - Representation，在Isovalue文本框里输入想设的值然后按回车即可。
顺带一提，如果你之前没看过《用Multiwfn+VMD做RDG分析时的一些要点和常见问题》（http://sobereva.com/291）的话强烈建议看一下，里面有很多重要的作图要点对于绘制IRI图也同样适用。
做过RDG分析的人都知道有个重要的RDG vs sign(λ2)ρ散点图。使用Multiwfn也可以绘制IRI vs sign(λ2)ρ散点图，也就是在Multiwfn计算IRI完毕后的后处理菜单中选-1 Draw scatter graph，然后就会看到下图
(, 下载次数 Times of downloads: 128)
此图里sign(λ2)ρ在-0.04~0.03 a.u.范围内有好几个spike，对应的是苯酚间的弱相互作用，在小于-0.25 a.u.的部分还有好多spike，对应的是化学键作用。关于散点图的意义和分析请看前文里提到的那些RDG/NCI分析文章，这里就不多说了。
如果把散点数据导出成output.txt文件，并且使用Multiwfn文件包里的examples\scripts\IRIscatter.gnu脚本，还可以得到彩色的散点图，不仅明显更好看，通过颜色明显更便于对照指认各个spike和等值面的对应关系。例子见Multiwfn手册4.20.4节。
之前有人问绘制IRI图的时候能否只显示特定部分，从而避免其它区域的等值面在分析时产生干扰。这是可以的，有几种做法：
(1)Multiwfn载入波函数文件后，在做IRI分析前，先进入主功能6，用选项-3或-4可以分别设置只保留哪些原子的贡献、扣除哪些原子的贡献。之后照常做IRI分析即可。
(2)当Multiwfn让你设置格点的时候，恰当定义要计算的空间范围，使计算格点数据的区域只涉及你感兴趣的区域。Multiwfn提供了丰富、灵活的格点设置方式，见Multiwfn手册3.6节的介绍，特别是其中10 Set box of grid data visually using a GUI window的设置方式可以在图形界面里直接设置盒子的位置和大小，非常直观。
(3)用Multiwfn的主功能13中的格点数据屏蔽功能可以只保留特定区域，参见Multiwfn手册4.13.4节的实例，在《使用Multiwfn结合CP2K通过NCI和IGM方法图形化考察固体和表面的弱相互作用》（http://sobereva.com/588）的3.6节也有这种做法的例子。
如果你只希望IRI展现弱相互作用区域，而不想把成键区域也表现出来，可以在Multiwfn做完IRI分析后选择屏幕上的选项9 Screen out covalent bond regions (set IRI to 100 for regions with sign(lambda2)rho < -0.1 a.u.)（此选项在2024-Mar-29及以后更新的Multiwfn里才有），之后sign(λ2)ρ小于-0.1 a.u.的区域内的IRI值就会被设为一个很大的值100。这样做的理由是形成化学键的区域通常电子密度都较大，因而sign(λ2)ρ比-0.1 a.u.更负，此区域的IRI设成很大的值后再绘制等值面图时就看不到这块区域了。之后你再导出func1.cub和func2.cub、用IRIfill.vmd脚本作图，就会发现成键区域确实已经看不见了，对苯酚的例子如下所示
(, 下载次数 Times of downloads: 79)
3.2 G-C二聚体的IRI平面填色图
这一节演示一下怎么获得前文中的G-C碱基对平面上的IRI填色图。注意读者务必用2022-Jul-4及以后更新的版本。
启动Multiwfn，载入本文文件包里的GC.fchk（这是Gaussian 16在B3LYP-D3(BJ)/6-311G**级别下优化G-C二聚体得到的）。在主功能0里可以看到此体系完全处于Z=0的XY平面上，因此我们将对这个平面作图。依次输入
4 //绘制平面图
24 //IRI
1 //填色图
直接按回车 //使用默认的格点数
0 //设置体系向四周的延展距离（用于定义作图范围）
2 //2 Bohr
1 //XY平面
0 //Z=0
此时图像蹦出来了。在上面点右键关闭图像，然后输入下面的内容改进作图效果
8 //显示键
14 //棕色
19 //设置色彩变化方式
2 //Reversed rainbow
4 //设置原子标签颜色
1 //红色
1 //设置色彩刻度下限和上限
0,2
-1 //重新绘制
此时看到下图，效果很理想。其中橙色+绿色区域对应存在明显相互作用的区域。
(, 下载次数 Times of downloads: 133)
4 周期性体系的IRI分析实例
Multiwfn也可以基于CP2K第一性原理程序产生的周期性体系的波函数绘制IRI图。读者请先仔细阅读《使用Multiwfn结合CP2K通过NCI和IGM方法图形化考察固体和表面的弱相互作用》（http://sobereva.com/588）了解怎么用CP2K+Multiwfn+VMD绘制周期性体系的RDG图，然后只需要把选择RDG分析的地方改成选择IRI分析、把source RDGfill.vmd改成source IRIfill.vmd即可获得周期性体系的IRI图。
这里举一个简单的绘制例子，对重构之后的硅(001)表面绘制IRI图考察成键情况。如果不知道表面重构是怎么回事，看下图
(, 下载次数 Times of downloads: 123)
本文文件包里Si_surface.inp是CP2K的输入文件，是对表面重构后的结构在PBE/pob-TZVP级别下算单点并得到molden文件的任务。文件包里的Si_surface.molden是这个任务产生的molden文件，笔者已按照http://sobereva.com/588所述的方式将晶胞信息作为[Cell]字段手动加入了其中。
启动Multiwfn，载入Si_surface.molden。在开始IRI计算之前，我们可以先看一下此体系的原子位置和晶胞，便于我们确定最佳的格点设定。输入0进入主功能0，然后选择菜单上的Other settings - Toggle showing cell frame，此时可以看到下图
(, 下载次数 Times of downloads: 138)
上图紫色的边框展现了CP2K做计算的时候的晶胞范围。结合上图的坐标轴，以及Multiwfn文本窗口显示的信息，我们可知此体系Z坐标最大的原子大约是Z=13 Bohr。
点击Return按钮关闭图形窗口回到Multiwfn主菜单，然后接着输入
20 //弱相互作用可视化分析
4 //IRI分析
9 //基于晶胞的平移矢量定义格点
按回车 //代表用(0,0,0)位置作为格点数据的原点（从前面的图可见这也正是当前晶胞的原点）
0,0,16 //格点所在的盒子的三个方向的边长。前两个0代表第1、2方向的边长恰等于晶胞的第1、2方向的边长，16代表在晶胞的第3个方向（此例对应Z方向）取16 Bohr作为盒子的第三个方向边长。此时的盒子足矣涵盖所有原子所在的区域，不可能造成等值面被截断，同时又避免了计算真空区浪费耗时
0.15 //格点间距。数值越小要算的点数越多，耗时越高，而图像越平滑
算完了之后还是像前面的例子一样选择3导出func1.cub和func2.cub，然后把它们挪到VMD目录下，把Multiwfn目录下的examples\IRIfill.vmd文件挪到VMD目录下，然后启动VMD并在文本窗口输入source IRIfill.vmd，马上就看到了图像。为了让效果更好，做如下调整
(1)在VMD的文本窗口输入pbc box把盒子边框显示出来
(2)对于当前这个体系，默认的IRI=1.0等值面不太合适，等值面偏大。因此进入Graphics - Representation，在Isovalue文本框里输入一个小一些的值0.7（可反复尝试来确定），此时图像效果较好
(3)为了让图像更有层次感，选Display - Depth cueing开启景深效果。然后选Display - Display Settings，把Cue Mode设为Linear，把Cue Start和Cue End分别设为1.75和2.50
最终的图像如下，两个视角的图一起给出了
(, 下载次数 Times of downloads: 115)
从上面的图可以看到每个Si-Si键都有相应的蓝色的IRI等值面出现，体现出这是化学键作用。在比较大的笼状区域里有很多绿色的IRI等值面，体现出这些区域中硅原子间形成了范德华作用。还有一些相对较小的笼和五元环状区域中IRI等值面颜色为明显的红色，直观地体现出这些区域存在显著的位阻作用。注意这种相对较小的笼、环是在硅表面重构形成新的Si-Si键后才出现的，体现出形成新键令能量降低是重构的驱动力，而在这个过程中也要同时克服造成新的位阻导致的能量升高因素。
5 IRI-π分析
在Multiwfn中有一套专门的分析π电子结构的流程，见《在Multiwfn中单独考察pi电子结构特征》（http://sobereva.com/432）以及笔者的专门的介绍论文Theor. Chem. Acc. 139, 25 (2020) https://doi.org/10.1007/s00214-019-2541-z。笔者将这种分析思路与IRI相结合，就诞生了IRI-π。简单来说IRI-π这个函数就相当于在计算IRI的时候只考虑π电子，或者说在Multiwfn里先把π轨道以外的轨道占据数清零然后照常计算IRI。IRI原文2.6节对IRI-π的特性做了专门的讨论，这里笔者只是简单说一下。
下图(a)和(b)分别是乙烯和乙炔的IRI-π等值面图。由于这个函数是专门展现π作用的，和位阻区域完全无关，所以按照下图的色彩刻度根据ρ着色即可，而没必要用sign(λ2)ρ来着色。由IRI-π等值面形状可见，乙烯这样一套π作用和乙炔这样两套π作用对应的等值面形状完全不同，前者只在π平面上、下方有等值面，而后者则是绕着键轴出现了环形等值面。因此根据等值面形状可以判断π作用类型。
(, 下载次数 Times of downloads: 122)
上图的(c)是电子结构很特殊的18碳环（cyclo[18]carbon）的ρ着色的IRI-π等值面图。关于这个体系笔者曾做过大量研究并发表了一系列研究文章，汇总见http://sobereva.com/carbon_ring.html。此体系有个特点是其中C-C键一长一短交替出现，由上图可见每个C-C键的IRI-π等值面都是环状的，说明每个C-C键都是两套π作用。但两类C-C键的π作用强度明显不同，较短的C-C键的环状IRI-π等值面颜色比较长的C-C键的更蓝，很直观地说明较短的C-C键的π作用区域电子密度更大，也因此π作用更强。
再看噻吩的例子。下图(a)和(b)是噻吩在不同等值面数值下的IRI-π等值面图，使用和上图一样的着色方式。由图(a)的等值面颜色可见Cα-Cβ的π作用比Cβ-Cβ的明显更强，这和一般常识以及Multiwfn可以计算的π电子的Mayer键级结论相符。对于图(b)展现的IRI-π=1.7等值面，哪怕不看颜色，光看等值面形状也能看出两种C-C键的显著区别。Cα-Cβ的IRI-π等值面都延伸到接近Cα和Cβ原子核的位置了，而Cβ-Cβ的IRI-π等值面只处在π作用区域处。因此，此例说明IRI-π函数在特定等值面数值下的形状也往往能判断π作用的强弱。
(, 下载次数 Times of downloads: 140)
Multiwfn具有强大的盆分析功能，对Multiwfn支持的任意实空间函数都可以获得其极大、极小点位置以及这些位置的具体函数值，见《使用Multiwfn做电子密度、ELF、静电势、密度差等函数的盆分析》（http://sobereva.com/179）。之后还可以用Multiwfn的拓扑分析功能对盆分析给出的极值点位置进一步refine得到更准确的结果。这种操作对于IRI-π也照样适用。上图(c)里每个小绿球是Multiwfn搜索出的IRI-π的极小点位置（可视为π作用区域最有代表性的位置），相应位置的π电子密度（ρ-π）和IRI-π数值都给出了。可见Cα-Cβ的IRI-π的极小点位置的密度显著大于Cβ-Cβ的，在定量上进一步体现出见Cα-Cβ的π作用更强。
IRI-π在Multiwfn中的计算也很容易，限于本文的篇幅在这里就不详细介绍了，大家可以直接看IRI分析的英文教程http://sobereva.com/multiwfn/res/IRI_tutorial.zip里的相应部分，写得非常详细易懂。
后来笔者发表了研究18碳环衍生物C18-(CO)n的成键特征和电子离域性的文章，见Chem. Eur. J., 28, e202103815 (2022)，在其中使用了IRI-π分析，所得图像如下所示，充分体现了IRI-π在研究pi作用特征、横向对比不同体系pi作用强度方面的价值。此文的深入浅出的介绍见《深入揭示18碳环的重要衍生物C18-(CO)n的电子结构和光学特性》（http://sobereva.com/640），非常欢迎大家阅读。此外，在《不寻常的环[18]碳前驱体C18Br6的电子结构和芳香性》（http://sobereva.com/664）里介绍的Chem. Eur. J., e202300348 (2023)文中还通过IRI-π考察了C18Br6的pi电子的特征、在《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（http://sobereva.com/696）里介绍的Inorg. Chem., 62, 19986 (2023)文中还通过IRI-π考察了18碳环的等电子体B6N6C6和B9N9的pi作用特征。非常推荐在引用Multiwfn原文、Chemistry Methods上的IRI原文的同时引用这几篇文章作为IRI-π的应用范例。
(, 下载次数 Times of downloads: 107)
6 总结
本文介绍了近来新提出的IRI函数。此函数是图形化展现化学体系中各种相互作用很理想的函数，其定义简单，计算快速，编程实现容易，图像效果好，普适性强，有公开且免费的计算程序，作图简单方便，对孤立体系和周期性体系都能用，因此IRI有着独特的意义和显著的实用性，是计算化学家的波函数分析工具箱中的有用的新工具，欢迎大家在未来的研究中广泛使用。另外，IRI的变体IRI-π给分析π相互作用又增添了与以往截然不用的新方法，很值得在未来研究π相互作用的场合使用。
作者Author: captain 时间: 2021-6-3 21:02
本帖最后由 captain 于 2021-6-3 21:03 编辑
看来大神的CP2K第一性原理程序培训班准备就绪了
作者Author: 蕊蕊花儿杰杰 时间: 2021-6-7 10:26
您好，我用GROMACS跑混合物的动力学，，可以用这个展现其中两个分子间的相互作用吗？可以展现过程的种相互作用随时间的变化吗？谢谢社长。新手小白，刚接触，不是很懂。
作者Author: snljty 时间: 2021-6-7 10:39
参考卢老师的博文《使用Multiwfn做aNCI分析图形化考察动态过程中的蛋白-配体间的相互作用》 http://bbs.keinsci.com/thread-21826-1-1.html。看了一下手册目前版本似乎不支持promolecular的IRI分析。
作者Author: sobereva 时间: 2021-6-15 19:09
考察这个用promolecular的aRDG，见
使用Multiwfn做aNCI分析图形化考察动态过程中的蛋白-配体间的相互作用
http://sobereva.com/591（http://bbs.keinsci.com/thread-21826-1-1.html）
使用Multiwfn研究分子动力学中的弱相互作用
http://sobereva.com/186
如果只需要考察一帧当中的分子间作用就够了，用promolecular的RDG，或者IGM：
通过独立梯度模型(IGM)考察分子间弱相互作用
http://sobereva.com/407（http://bbs.keinsci.com/thread-9472-1-1.html）
IRI目前不支持promolecular的形式。
作者Author: yinlina 时间: 2021-9-18 15:13
老师您好，我将分子拖进去输入指令后，不显示有4IRI分析这个选项，而且multiwfn里面没有相应的IRI脚本，请问是我的版本问题吗？
作者Author: sobereva 时间: 2021-9-19 05:23
文中都特意说了
读者务必使用2021年2月及以后更新的Multiwfn版本
作者Author: zhaoxy008 时间: 2021-11-27 18:32
请问老师，作用区域的电子密度从那里能得到数据？
作者Author: sobereva 时间: 2021-11-28 14:17
直接看IRI等值面颜色对照色彩刻度条便知。当然Multiwfn里也有一大堆其它方式可以衡量，比如做AIM拓扑分析看BCP处的电子密度、绘制原子连线间电子密度曲线图等等
作者Author: zhaoxy008 时间: 2022-2-16 17:09
请教Sob老师：画红圈的对应关系从哪里看呢？
作者Author: sobereva 时间: 2022-2-17 08:32
用Multiwfn获得各个等值面里面的BCP的电子密度数值，和散点图的spike位置对照便知
使用Multiwfn做拓扑分析以及计算孤对电子角度
http://sobereva.com/108
作者Author: sun666 时间: 2022-3-16 19:16
小白请教老师，静电作用也属于弱相互作用。对于阳离子与阴离子之间的IRI分析，好像只有范德华作用而没有静电。那么如何判断静电主导还是范德华力主导呢。
作者Author: sobereva 时间: 2022-3-17 03:43
凭常识就知道是静电作用主导
IRI只是体现作用区域的特征，但整体带电荷的阴、阳离子间，静电作用是长程的，并不仅体现在接触区域上
可以做能量分解来说明静电和范德华作用各自具体贡献多少，而且下文的EDA-FF能量分解还可以告诉你各个原子具体的贡献
使用Multiwfn做基于分子力场的能量分解分析
http://sobereva.com/442（http://bbs.keinsci.com/thread-10907-1-1.html）
作者Author: 赵云跳槽 时间: 2022-3-22 19:16
本帖最后由 赵云跳槽 于 2022-3-22 19:40 编辑
(, 下载次数 Times of downloads: 109)
Sob老师好，我看了教程，有几个点没看懂？
1. 2号位置的吸引子，为何重新处理后，数值变化那么大？
原则上应该有变化，但从原来的0.1348到0.233E-5，这个跨度有点大
2. 通过分析显示，3号位置的值为0.3911，为何老师你那给的是0.050？
3. 如果我获得了吸引子位置，我直接新启动Multiwfn，使用主功能1，然后输入坐标，这样也给出了IRI-pi和ρ-pi的值，这样是不是更简单点？
作者Author: sobereva 时间: 2022-3-23 10:35
1 原理上，如果2号吸引子位置是精确的，其IRI-pi应当精确为0，对应此处pi电子密度的梯度为0。盆分析时候定位精度取决于格点间距，只可能落在与精确位置最近的格点上，故注定不可能定位特别精确，肯定偏离实际位置一些，这会导致IRI-pi不那么为0。但pi电子密度变化没有IRI-pi那么剧烈，所以用盆分析得到的相对粗糙的IRI-pi吸引子位置计算pi电子密度的准确性还可以接受，前提是格点间距不能用得太大。
2 IRI原文里用的格点是high quality grid，教程是medium quality grid。无所谓，3号点没什么实际意义，而且其IRI-pi值也不是关注的对象。文中仅仅是标了一下而已，都没有去提。
3 不是。用拓扑分析对IRI-pi吸引子的位置做refine，对于更精确地得到pi电子密度有意义。如果你的要求很低，用主功能1直接算盆分析给出的粗糙的吸引子位置的pi电子密度也能凑合。但如果好几个pi作用区域，pi作用强度相仿佛，你想根据各自的pi电子密度进行对比，那最好还是用拓扑文件做refine得到的IRI-pi临界点位置来算pi电子密度，否则可能造成误判。
作者Author: 忘不掉的回忆 时间: 2023-3-13 16:32
老师好，我想请教一个问题。
您在本文中提到“使用Multiwfn文件包里的examples\scripts\IRIscatter.gnu脚本，还可以得到彩色的散点图，”
我做出了这个彩色的散点图，为了更明确的说明spike的范围，我想着和黑白散点图一样，
在彩色散点图的横坐标轴上的主刻度之间添加一个分刻度。
我自己试着在IRIscatter.gnu脚本中添加，但是没成功。
所以我想请问老师，应该如何更改？谢谢老师，谢谢老师。
(, 下载次数 Times of downloads: 59)
(, 下载次数 Times of downloads: 65)
作者Author: sobereva 时间: 2023-3-13 21:49
还不如把主刻度间隔改小
作者Author: 忘不掉的回忆 时间: 2023-3-14 08:34
好的，谢谢老师。
作者Author: flying~ 时间: 2023-4-13 04:56
卢老师好，关于IRI方法，我想请教一下，如果是用其他DFT软件计算的电子密度，比如FHI-aims，能不能通过您引入的电子密度及其微分的IRI公式进行计算呢，可不可以理解为DFT获得电子密度后，进一步计算可以通过数学软件实现？这方面您有没有经验或者有可以实现的code没？谢谢！
作者Author: sobereva 时间: 2023-4-13 21:57
如果程序只能产生电子密度格点数据，你可以尝试用三维B-spline插值通过有限差分计算各个位置的电子密度的导数，然后代进IRI公式。如果程序直接就能产生电子密度梯度矢量的格点数据，那做很简单的运算就能得到IRI。
作者Author: flying~ 时间: 2023-4-13 23:12
好的，谢谢卢老师指点，我再和计算合作者沟通一下，看看他们会不会操作，发表文章时候肯定要引用您的这篇IRI的文章。
作者Author: xptracy 时间: 2023-6-3 13:25
(, 下载次数 Times of downloads: 59)
请问下vmd画IRI，导出pdb的时候，出现libpng warning: iCCP: known incorrect sRGB profile，导不出pdb是什么问题？
作者Author: sobereva 时间: 2023-6-4 02:40
没见过。重启VMD重新操作一遍再试，或者把保存路径设到诸如D盘再试
作者Author: Melvin 时间: 2024-2-22 20:55
想请教sob老师一个问题，比如研究十几个原子的有机分子的几个构象。用IRI看有些构象只有共价作用，有些构象分子内还存在范德华作用，但是存在范德华作用的构象是不是也不一定是最低能的构象，是否能用IRI说明或者解释不同构象之间能量差异的原因呢？
作者Author: sobereva 时间: 2024-2-22 21:47
得把所有因素全都考虑进去，诸如静电吸引/排斥作用、位阻作用等等。如果牵扯到溶剂模型，诸如带显著电荷部分是露在外面还是包埋在里面还会对溶解自由能有很大影响。
例如下面这篇文章在讨论片段间作用能的时候就把不同因素做了对比分析
直观解释分子间相互作用如何影响不对称催化：Nature Chemistry上一个很好的IGMH分析范例
http://sobereva.com/700（http://bbs.keinsci.com/thread-43122-1-1.html）
作者Author: Melvin 时间: 2024-2-22 23:07
本帖最后由 Melvin 于 2024-2-22 23:09 编辑
感谢sob老师的解答。是的，我就是看到了您最近的这篇讲解Nature Chemistry的博文学到了很多但产生的一点疑惑，您博文的案例中是IGMH分别对两个过渡态统计了相互作用说明能量差异。
我想详细再请教一下不同的地方是，如果是一个分子的两个构象，结构差别较大，比如气相条件下乙醚扭转120°C-O键优化后形成的不同构象，扭转后的结构虽然所有的共价作用一样但是静电/位阻/范德华作用都发生了变化，是不是没法像您博文里这个情况一样仅靠IGMH或IRI图很好的去说明呢？IGMH和IRI只能说明共价/位阻/范德华作用,因此还需要补充静电作用的分析？
作者Author: sobereva 时间: 2024-2-22 23:16
不同情况适合的分析方式不同，如果分子内相互作用看起来没显著差异，也可以考虑用刘述斌的能量分解，参考Multiwfn手册4.21.2节的例子和里面引用的相关文章
作者Author: Melvin 时间: 2024-2-22 23:19
谢谢sob老师！
作者Author: sobereva 时间: 2024-3-29 20:45
今日更新的Multiwfn给IRI分析后处理菜单中加入了选项9 Screen out covalent bond regions (set IRI to 100 for regions with sign(lambda2)rho < -0.1 a.u.)。选了这个再导出格点数据、照常绘制IRI图后，共价键的区域就看不到了。相应地在本文3.1节末尾举了例子。
作者Author: pcw94121 时间: 2024-7-9 10:49
苯酚二聚体的等值面图
4 //绘制平面图
24 //IRI
1 //填色图
直接按回车 //使用默认的格点数
0 //设置体系向四周的延展距离（用于定义作图范围）
2 //2 Bohr
1 //XY平面
0 //Z=0
此时图像蹦出来了。在上面点右键关闭图像，然后输入下面的内容改进作图效果
8 //显示键
14 //棕色
19 //设置色彩变化方式
2 //Reversed rainbow
4 //设置原子标签颜色
1 //红色
1 //设置色彩刻度下限和上限
0,2
-1 //重新绘制
为什么我做到设置色彩刻度下限和上限输入0，2之后就闪退
作者Author: sobereva 时间: 2024-7-9 13:55
分清楚半角和全角逗号
作者Author: 赵云跳槽 时间: 2024-10-8 21:19
IRI-pi绘图，有两种着色方式
1.根据ρ着色
2.用sign(λ2)ρ着色
两者都是主功能100, 1，前者是键入1,24，而后者键入15,24 ？不知道理解对不对？
作者Author: 虎王 时间: 2024-10-8 23:15
IRI是否使用于分析第一激发三重态、第二激发三重态的分子间作用？
作者Author: sobereva 时间: 2024-10-9 04:21
能得到相应激发态的电子密度即可分析
作者Author: sobereva 时间: 2024-10-9 04:23
对
作者Author: 虎王 时间: 2024-10-9 15:47
通过在TDDFT计算中，于命令行加入density指令，从而把目标激发态的电子密度存进chk吗？
作者Author: sobereva 时间: 2024-10-10 06:29
仔细看下文。必须波函数文件里有对应激发态的自然轨道才行
详谈Multiwfn支持的输入文件类型、产生方法以及相互转换
http://sobereva.com/379（http://bbs.keinsci.com/thread-6020-1-1.html）
最简单的做法，TDDFT的时候用out=wfn关键词产生wfn文件，里面会记录root指定的激发态的自然轨道，此wfn文件给Multiwfn用，做的分析就是对激发态波函数分析的
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |
