IRI描述符:在IRI框架内定量描述化学体系中的相互作用特征
IRI描述符:在IRI框架内定量描述化学体系中的相互作用特征
IRI descriptors: Quantitatively characterizing interactions in chemical systems within the IRI framework
文/Sobereva@北京科音 2026-Nov-6
1 前言
我在2021年提出的相互作用区域指示函数(interaction region indicator, IRI)是直观展现化学体系中相互作用非常好用也很流行的方法(截止到撰文时已经被2000多篇文章使用),在笔者的以下文章里有充分的介绍,如果不了解IRI的话一定要先仔细看看:
IRI原文:Chemistry—Methods, 1, 231-239 (2021) https://doi.org/10.1002/cmtd.202100007
使用IRI方法图形化考察化学体系中的化学键和弱相互作用
http://sobereva.com/598(http://bbs.keinsci.com/thread-23457-1-1.html)
Angew. Chem.上发表了全面介绍各种共价和非共价相互作用可视化分析方法的综述
http://sobereva.com/746(http://bbs.keinsci.com/thread-54361-1-1.html)
一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表!
http://sobereva.com/667(http://bbs.keinsci.com/thread-37629-1-1.html)
IRI以及其它许多相关方法(IGMH、mIGM、amIGM、IGM、RDG/NCI等),都使用某种函数通过其等值面展现相互作用区域,并将sign(λ2)ρ按照数值以不同颜色映射到等值面上用于直观区分相互作用的特征。对于IRI,标准的sign(λ2)ρ色彩刻度条以及含义如下所示
注:也不要太教条地理解不同sign(λ2)ρ数值范围对应的相互作用特征,始终要结合相互作用的物理本质理解。比如一个净电荷分布较为分散、较大的阳离子和阴离子形成的二聚体,二者之间显然有显著的长程静电吸引作用,但由于二者间接触区域的电子密度(ρ)较低,sign(λ2)ρ也会在0附近,因而着色主要为绿色,就不能简单认为是范德华作用主导它们之间的相互作用了。
一般说的IRI图,即sign(λ2)ρ着色的IRI函数的等值面图,虽然能直观展现出在哪里有什么类型的相互作用,但不便于确切地定量讨论和定量横向对比相互作用情况。对于弱相互作用,如《Multiwfn支持的弱相互作用的分析方法概览》(http://sobereva.com/252)所示,有一些可以定量讨论相互作用的方法,但都各有明显局限性,例如:
• AIM理论定义的键临界点位置上的属性虽然可以较好地用于相互作用分类、对比强度(例如如http://sobereva.com/513所述,氢键键能正比于BCP位置的电子密度),但是有不少不可忽视的弱相互作用并没有对应的BCP,就完全没法用这种方法表征了。而且BCP仅仅是三维空间中的一个点,光靠这个位置的属性在很多情况下并不能比较全面地反映相互作用情况,比如作用区域往往很大的范德华吸引和pi-pi堆积。
• 对IRI做拓扑分析然后考察其临界点的属性,例子见Multiwfn手册4.2.11节。这种做法虽然可以对很多没有BCP而能被IRI展现出的相互作用进行定量考察,但能反映的也仅仅是作用区域的个别点的属性,局限性还是明显。
• 对IRI、IGMH、NCI等方法有了解的人都知道格点数据可以以散点图形式描述。对于IRI来说横坐标是sign(λ2)ρ值,纵坐标是IRI值,图中的spike体现明显的相互作用的存在,根据spike的横坐标的具体位置可以定量对比相互作用。但当体系中相互作用出现位点较多时,这种图会很乱,难以指认各个spike对应空间中什么区域的相互作用。在我来看散点图是没什么价值的,尤其是有了本文介绍的IRI描述符后。
• 包括《使用sobEDA和sobEDAw方法做非常准确、快速、方便、普适的能量分解分析》(http://sobereva.com/685)介绍的sobEDA在内的绝大多数能量分解方法,可以给出片段间相互作用能和物理成份来定量对比,但难以单独考察体系中不同位点的相互作用(有时候靠恰当构造模型体系能实现,但比较麻烦)。
之前有不少Multiwfn的用户曾问我怎么定量对比IRI图展现的不同位置的相互作用。容易理解,若在不同位点对应的IRI等值面内进行恰当的积分,显然有可能获得一些有意义的量用于做这种定量对比。实际上之前已经有一些文章对NCI方法使用的RDG函数的等值面内的某些函数进行积分来讨论问题,只不过那些方法各有一些局限性或不合理性。受到那些文章的启发,并结合IRI函数的特征,我提出“IRI描述符”(IRI descriptors),里面定义了一系列具体的量,将在下面专门介绍一下,然后会给出用Multiwfn的计算实例。我实测发现这些IRI描述符对于定量分析讨论不同相互作用颇有益处,推荐大家在实际研究中充分运用。
注:IRI描述符尚未发表在期刊上,以后会专门撰文发表,届时会更新本文并提及。目前使用此方法研究时,请引用IRI原文、Multiwfn程序原文及Multiwfn手册3.200.14节(包含了IRI描述符和做这个分析的域分析模块的说明)。也很推荐同时引用http://sobereva.com/746和http://sobereva.com/667中提到的介绍了IRI方法的笔者的综述。
2 IRI描述符的定义
IRI等值面通常取1.0(也是Multiwfn自带的VMD的IRI作图脚本IRIfill.vmd的默认值),每个这样的IRI等值面内包围的区域叫做IRI域(IRI domain)。按照sign(λ2)ρ的数值(a.u.),我经验性地划分以下区间,这使得一个IRI域可能包含多个不同的子域:
• strong attraction(satt,强吸引):sign(λ2)ρ < -0.05
• moderate attraction(matt,中等吸引):-0.05 <= sign(λ2)ρ < -0.013
• weak attraction(watt,弱吸引):-0.013 <= sign(λ2)ρ < +0.009
• repulsion(rep,互斥):+0.009 <= sign(λ2)ρ
• attraction(att):是strong attraction、moderate attraction、weak attraction的并集
• total(tot):是attraction和repulsion的并集
satt子域部分主要对应共价键,以及弱相互作用中相对来说非常强的那些,它们伴随着显著的共价性,可参见http://sobereva.com/513里介绍的我的专门研究各种氢键特征的文章J. Comput. Chem., 40, 2868 (2019)中的讨论。matt主要对应常规强度的氢键、卤键、二氢键等显著弱相互作用,watt主要对应色散这种很弱的相互作用。注意体积较大的阴-阳离子间尽管相互作用能挺大的,但由于作用区域电子密度较低,也可能落在matt甚至watt里。
在不同区域内都可以对特定函数进行积分。如果被积分的是处处为1的常数,则积分结果就是相应区域的体积(V)。被积分对象如果是电子密度,得到的就是参与相应作用的电子数(q),和作用强度有密切关联。积分其它函数,如动能密度、势能密度、能量密度、变形密度等,原理上也可以,但本文暂不考虑。这样就有了许多IRI描述符。例如某个IRI域的V_watt就是其中对应弱相互作用部分的体积,q_matt是描述中等吸引作用区域内的电子数,q_tot显然对应整个IRI域内包含的电子数。以此类推。
V_att/V_tot用来衡量一个IRI域内体现吸引作用的体积占总体积的比例,q_att/q_tot用来衡量一个IRI域内起吸引作用的电子数占总电子数的比例。
上述的q和V都是广延量。为了能描述单位体积对应的作用强度,我给IRI描述符中还定义了以下三项,量纲都是电子密度。相对于AIM分析中经常只考察BCP这一个点的电子密度,这些描述符展现的信息明显更全面、充分
• q_att/V_att:IRI域里体现吸引作用的区域中的平均电子密度
• q_rep/V_rep:IRI域里体现互斥作用的区域中的平均电子密度
• q_tot/V_tot:整个IRI域里平均电子密度
这些量很有意义。比如有一个中等强度氢键对应的较小的IRI等值面和一个pi-pi堆积对应的特别广阔的IRI等值面,由于前者的体积远小于后者的,导致前者的q_att可能比后者更小,且实际吸引作用强度可能没后者大;但前者的q_att/V_att必定比后者大得多,体现出单位大小的局部区域内前者强得多。
IRI域里sign(λ2)ρ最负值和最正值也算是描述符,分别记为sl2r_low和sl2r_high。
《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》(http://sobereva.com/621)介绍的笔者提出的IGMH方法是和IRI同等重要、同等流行的相互作用可视化方法,但各有所长。肯定有人会想能不能对IGMH方法定义的δg_inter等值面也计算上述描述符,实际上是不合适的。其中很关键的一个原因是IRI=1.0 a.u.等值面通常可以把所有相互作用区域全都较好展现出来,而δg_inter则根据等值面数值的大小区分展现不同强度的相互作用。
3 IRI描述符在Multiwfn中的计算
这一节介绍下IRI描述符的计算实现。笔者开发的强大的波函数分析程序Multiwfn很久以前的版本就包含了域分析(domain analysis)功能,作为主功能200的子功能14。域分析功能可以对任意函数的特定等值面内包围的空间(域)进行函数的积分和计算体积,任意Multiwfn支持的函数,甚至靠格点数据文件记录的外部函数,都可以作为被积函数,并且还支持对等值面内sign(λ2)ρ处于特定数值范围的子域进行积分,利用这些特征显然能直接实现IRI描述符的计算。为了让大家计算前述的IRI描述符尽可能方便,Multiwfn的2026.11.3版对域分析功能做了很大改进,直接提供了现成选项能够一次性把所有用户选定的各个域的所有IRI描述符全都算出来,顺便还给出各个域结果的总和。域分析功能的详细介绍见Multiwfn手册3.200.14节。
Multiwfn的域分析功能是基于均匀分布的格点的,用它计算IRI描述符的流程是:首先对体系划定一块三维空间区域(盒子),对里面均匀分布的各个格点位置计算IRI函数值,然后用笔者提出的一种算法产生出各个域,每个域对应一批格点。之后对每个域中的各个格点计算出电子密度和sign(λ2)ρ值,就可以进一步积分得到各种IRI描述符。显然,格点间距越小,积分越精确,但格点数越多,计算越耗时。为了精度足够,建议格点间距用0.15 Bohr及以下。
下面就举一些典型例子,示例利用域分析功能计算IRI描述符,并展现IRI描述符的实际价值。Multiwfn程序可以在官网http://sobereva.com/multiwfn免费下载,读者请使用2026.11.3及以后版本。如果对Multiwfn不了解,请阅读《Multiwfn FAQ》(http://sobereva.com/452)、《Multiwfn入门tips》(http://sobereva.com/167)、《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》(http://sobereva.com/379)。
4 例:苯酚二聚体
下图是对Multiwfn自带的苯酚二聚体波函数文件examples\phenoldimer.wfn按照http://sobereva.com/598所示的标准流程绘制的sign(λ2)ρ着色的IRI等值面图,使用Multiwfn自带的IRIfill.vmd脚本绘制。后文若未注明,用的都是基态极小点结构、IRI等值面数值为1.0 a.u.、sign(λ2)ρ着色范围为-0.04到0.02 a.u.按蓝-绿-红变化,且0对应正绿色。
由上图可清楚看到各种作用区域。两个羟基间的青蓝色等值面对应明显的氢键,绿色大等值面对应色散作用,深蓝色对应共价键作用,苯环中央的红色等值面对应位阻作用。下面我们对对应不同相互作用的IRI等值面内进行积分获得IRI描述符。
启动Multiwfn,然后输入
examples\phenoldimer.wfn
200 //其它功能(Part 2)
14 //域分析
从屏幕上的选项文字可见当前默认设置:
2 Select real space function to be calculated for option 1, current: 24
3 Set criterion for defining domain, current: <1.0
即默认就会把Multiwfn中的第24号实空间函数(对应于IRI)数值小于1.0的部分作为被分析的域。利用这两个选项,你可以把其它任意函数大于或小于某个值的部分设成被分析的域,十分灵活和普适。
接着输入
1 //开始计算格点数据并产生域
-10 //设置延展距离,即各个方向上盒子边界距离体系最边缘的原子的距离
1 //延展距离为1 Bohr。故意设置一个较小的数值是为了节约盒子尺寸,从而减少要算的格点数。由于从当前体系等值面图可以看到等值面不会出现在体系外围,所以留一个较小的延展距离就够了
2 //Medium-quality grid(中等质量格点),格点间距对应0.1 Bohr
关于格点设置这里专门提一下。如果你对结果精度要求不高,或者只是试算的目的,或者对大体系发现用中等质量格点耗时高得不可接受或者由于内存不够而崩溃(内存消耗量正比于格点数),也可以选对应格点间距为0.2 Bohr的low-quality grid,耗时和内存消耗量只有medium-quality grid的1/8。
马上在屏幕上看到了域的列表,如下所示,有好几十个,各个域包含的格点数、体积、几何中心都给出了
Domain N_Grids Volume(Angs^3) X,Y,Z of domain center (Angs)
1 1 0.0001 3.52412 -2.80817 -4.16865
2 326 0.0483 3.79007 -2.32688 -3.44626
3 1 0.0001 1.32412 1.19183 -2.96865
4 1 0.0001 -6.97588 2.59183 -2.66865
...略
现在看到一个功能丰富的后处理菜单,从文字提示可见,可以对域进行积分、可视化、按照域的体积排序、计算IRI描述符、把域导出成cube或txt文件,等等。我们现在要做的是确认我们感兴趣的域的序号。选择3 Visualize domains,会看到如下观看域的界面。此图上每个浅绿色小圆球是一个域的几何中心,上面的文字是这个域的序号。在右下角的列表里可以选择显示哪个域,下图中我选了32号,可见它正对应于苯酚间色散作用的IRI域,这个域里面的格点都直观显示了出来。从图中显示的域的中心的位置结合之前的IRI等值面图还可以立刻判断出20号域对应氢键,18号域对应苯环的位阻,16号对应O-C化学键,等等。
点击右上角RETURN按钮关闭图形窗口,然后选择6 Calculate IRI interaction descriptors for domains,再输入域的序号就可以对指定的域计算IRI描述符了。这里我们取前面提到的几个典型的域来计算,输入32,20,18,16之后马上看到以下信息,每一列的含义提示得非常明白
Meaning of output:
satt: Strong attraction (sign(lambda2)rho< -0.05 a.u.)
matt: Moderate attraction (-0.05 <= sign(lambda2)rho < -0.013 a.u.)
watt: Weak attraction (-0.013 <= sign(lambda2)rho < +0.009 a.u.)
rep: Repulsion (+0.009 a.u. <= sign(lambda2)rho)
att (attraction) = satt + matt + watt
tot (total) = att + rep
Number of electrons (q):
Domain q_satt q_matt q_watt q_att q_rep q_tot q_att/q_tot
32 0.00000 0.00000 0.01524 0.01524 0.00000 0.01524 1.00000
20 0.00000 0.00325 0.00000 0.00325 0.00000 0.00325 1.00000
18 0.00000 0.00000 0.00000 0.00000 0.00203 0.00203 0.00000
16 0.08597 0.00000 0.00000 0.08597 0.00000 0.08597 1.00000
Total 0.08597 0.00325 0.01524 0.10447 0.00203 0.10650 0.98098
Volumes (Angstrom^3):
Domain V_satt V_matt V_watt V_att V_rep V_tot V_att/V_tot
32 0.00000 0.00000 0.50353 0.50353 0.00000 0.50353 1.00000
20 0.00000 0.01956 0.00000 0.01956 0.00000 0.01956 1.00000
18 0.00000 0.00000 0.00000 0.00000 0.01556 0.01556 0.00000
16 0.04668 0.00000 0.00000 0.04668 0.00000 0.04668 1.00000
Total 0.04668 0.01956 0.50353 0.56977 0.01556 0.58533 0.97342
Ratios between q and V (Bohr^3) correspond to average electron density (a.u.)
sl2r_low and sl2r_high: Lowest and highest sign(lambda2)rho (a.u.)
Domain q_att/V_att q_rep/V_rep q_tot/V_tot sl2r_low sl2r_high
32 0.00449 NaN 0.00449 -0.00622 0.00628
20 0.02465 NaN 0.02465 -0.02703 -0.02189
18 NaN 0.01929 0.01929 0.01535 0.02129
16 0.27293 NaN 0.27293 -0.35250 -0.23603
从q和V的信息可以看到,表现典型色散作用的32号域中的电子完全归属于watt,表现普通氢键的20号域中的电子完全归属于matt,表现典型位阻的18号域中的电子完全归属于rep,对应典型共价键的16号域中的电子完全归属于satt。这充分体现了我根据sign(λ2)ρ的数值范围所做的watt、matt、satt、rep的子域的划分非常合理。
再更具体对比一下不同域的IRI描述符的差异。描述32号域中起吸引作用的电子数q_att=0.01524比20号域的q_att=0.00325大得多,为什么色散作用区域里的电子数比氢键的大还那么多?关键在于前者的域体积很大,V_att=0.50353 Å^3,而后者只有0.01956 Å^3,悬殊的体积差异从前面给出的IRI等值面图上就可以很直观地看出来。而32号域里表现吸引作用区域的平均电子密度为q_att/V_att=0.00449 a.u.,比20号域的0.02465 a.u.小得多得多,非常明确地体现出单位作用区域中氢键是远远强于色散作用的。注:按照《使用Multiwfn做拓扑分析以及计算孤对电子角度》(http://sobereva.com/108)的方法用Multiwfn算出来的这个氢键的BCP处的电子密度是0.02628 a.u.,与氢键域的q_att/V_att很相近。
描述C-O化学键的16号域的q_att/V_att=0.27293 a.u.,比氢键的域大了一个数量级,体现出化学键和弱相互作用的强度差异很悬殊。这个键的BCP处的电子密度是0.28985 a.u.,也和q_att/V_att很相近。
完全描述吸引作用的32、20、16号域里面att和tot是相同的,因为它们都没有表现互斥作用的rep部分存在,也因此它们的q_rep/V_rep项显示的是NaN(not a number),因为作为分母的描述互斥部分的体积V_rep为0。
完全描述六元环位阻的18号域由于体积很小(仅V_rep=V_tot=0.01556 Å^3),因此里面的电子数也只有q_rep=q_tot=0.00203。其平均电子密度q_rep/V_rep=q_tot/V_tot=0.01929 a.u.和对应氢键作用的20号域的0.02465 a.u.相仿佛。
从输出信息中sl2r_low和sl2r_high还可以清楚地看到各个域里sign(λ2)ρ的下限和上限是多少。比如对应位阻的18号域的sign(λ2)ρ处在0.01535到0.02129 a.u.之间。顺带一提,JCC, 46, e70268 (2025)文中把RDG等值面(在弱相互作用区域中它和IRI等值面很相似)当中sign(λ2)ρ在-0.02到+0.02 a.u.的部分都算作范德华作用区域了,这明显非常不妥、范围太广了,这会导致典型、鲜明的位阻作用区域里有一大块都被误算成了范德华作用。而我根据经验把watt区域定义成sign(λ2)ρ处于-0.013到+0.009 a.u.则很合适,能把范德华作用充分纳入又不至于把位阻作用算进去。
此例反映出靠IRI域积分得到IRI描述符的方式很好用。只要有IRI等值面出现,就能计算里面包围的域的特征来描述相互作用。
5 例:氢键二聚体CH4...OH2和HOH...OH2
这一节通过IRI描述符对比同样是氢键作用的两个体系,CH4...OH2和HOH...OH2。这两个体系分别是http://sobereva.com/513介绍的我的J. Comput. Chem., 40, 2868 (2019)文章中的体系,在原文优化的结构下我在B3LYP/def-TZVP下用Gaussian 16算单点后得到的fch文件是http://sobereva.com/attach/786/file.zip里的CH4_OH2.fch和HOH_OH2.fch。按照常规流程绘制的它们的IRI图如下所示,从颜色上已经可以看出HOH...OH2的氢键更强,这里我们通过IRI描述符把这个差异定量化表达一下。
启动Multiwfn,载入CH4_OH2.fch,然后依次输入
200 //其它功能(Part 2)
14 //域分析
1 //开始计算格点数据和产生域,默认对IRI<1.0来做
在现在看到的设定格点的界面里,你可以像上个例子一样去操作。但由于氢键等值面挺小,如果你对域积分精度要求比较高,建议这里选High-quality grid,格点间距对应很小的0.06 Bohr,足够积分得很准确。
后处理菜单中进入选项3检查一下域的图像是否正常以及域的序号,再进入选项6计算IRI描述符,输入对应氢键的域的序号,下面看到以下信息
Number of electrons (q):
Domain q_satt q_matt q_watt q_att q_rep q_tot q_att/q_tot
9 0.00000 0.00000 0.00150 0.00150 0.00000 0.00150 1.00000
Total 0.00000 0.00000 0.00150 0.00150 0.00000 0.00150 1.00000
Volumes (Angstrom^3):
Domain V_satt V_matt V_watt V_att V_rep V_tot V_att/V_tot
9 0.00000 0.00000 0.03537 0.03537 0.00000 0.03537 1.00000
Total 0.00000 0.00000 0.03537 0.03537 0.00000 0.03537 1.00000
Ratios between q and V (Bohr^3) correspond to average electron density (a.u.)
sl2r_low and sl2r_high: Lowest and highest sign(lambda2)rho (a.u.)
Domain q_att/V_att q_rep/V_rep q_tot/V_tot sl2r_low sl2r_high
9 0.00629 NaN 0.00629 -0.00697 -0.00551
类似地,对HOH...OH2的氢键也这么计算IRI描述符。结果为:
q_att=q_tot=0.00284
V_matt=V_tot=0.01764 Å^3
q_tot/V_tot=0.02382 a.u.
可以看到,HOH...OH2氢键比CH4...OH2氢键的IRI域体积小了约一倍,但电子数却多了90%,并由此导致IRI域里的平均电子密度前者是后者的3.8倍,显然前者的氢键明显更强。
这个例子和其他人的JCC, 46, e70268 (2025)文中的发现相呼应。那篇文章计算了一系列只有一个明显作用位点的弱相互作用体系,氢键、卤键、磷键、pi-氢键体系都有,发现规律:对于越强的弱相互作用,静电作用贡献比例越大而色散作用贡献比例越小(这点和我的JCC, 40, 2868 (2019)文章中由弱到强不同氢键的研究完全吻合),并且RDG等值面里电子积分值(q_NCI)越大、等值面内体积(V_NCI)越小、q_NCI/V_NCI越大,因而作用区域的电子越呈现空间上定域的分布。故q_NCI/V_NCI可以与作用类型相挂钩。本节的例子体现出IRI等值面里的积分值也完全有相同的规律。故读者可以根据IRI描述符q_att/V_att区分作用类型,数值越大暗示静电作用相对于色散作用的比例越大,且作用越强(除非等值面体积相差过于悬殊)。顺带一提,我的JCC, 40, 2868 (2019)报道的CH4...OH2氢键相互作用能为-0.61 kcal/mol,静电和色散分别贡献-0.79和-0.84 kcal/mol;HOH...OH2氢键相互作用能为-4.96 kcal/mol,静电和色散分别贡献-8.26和-2.46 kcal/mol,可见二者的静电/色散贡献比例差异极大。
6 例:二醇内氢键
一些明显的相互作用有可能是没有BCP的,因而靠AIM拓扑分析完全无法研究其特征,但是通过IRI则总是可以将之展现出来,并且可以通过IRI描述符定量讨论。这一节就以两个二醇内氢键体系为例对比一下。这两个体系分别是乙二醇和丙二醇,在B3LYP/6-311G**下用Gaussian 16做opt freq后得到的fch文件是http://sobereva.com/attach/786/file.zip里的ethanediol.fch和propanediol.fch。用标准方式绘制IRI图,如下所示。
上图能明显看到两个体系都有内氢键。氢键是有取向性的,由于乙二醇的内氢键键角偏离直线太大,导致氢键较弱,氢键的吸引作用区域颜色只是稍微发青,而丙二醇内氢键键角明显更接近直线,因而氢键更强,氢键的吸引作用区域比较蓝。此外,还能看到丙二醇的氢键等值面旁边有一块橙红色体现位阻的等值面,而乙二醇的这两部分结合为了一个等值面的不同区域。对丙二醇可以用Multiwfn的拓扑分析功能搜索到内氢键的BCP以及对应位阻区域的RCP,但是对乙二醇却无法搜索到它们,这不是程序的问题,而是电子密度分布特征使然。
下面我们对内氢键起吸引作用的区域,以及伴生的分子内互斥作用的区域,通过IRI描述符来定量考察一下。按照上一节的例子操作,结果如下:
*乙二醇(11号是对应氢键+位阻的域)
Domain q_satt q_matt q_watt q_att q_rep q_tot q_att/q_tot
11 0.00000 0.00174 0.00015 0.00189 0.00207 0.00397 0.47773
Domain V_satt V_matt V_watt V_att V_rep V_tot V_att/V_tot
11 0.00000 0.01783 0.00176 0.01959 0.01946 0.03905 0.50164
Domain q_att/V_att q_rep/V_rep q_tot/V_tot sl2r_low sl2r_high
11 0.01433 0.01577 0.01505 -0.01626 0.01829
可以看到这个域里面处于吸引部分和处于互斥部分的电子数分别是0.00189和0.00207,吸引和互斥部分的体积分别是0.01959和0.01946 Å^3,都相仿佛。这和前面IRI等值面图上看到的情况差不多。
*丙二醇(9号是对应氢键的域,15号是对应位阻的域)
Domain q_satt q_matt q_watt q_att q_rep q_tot q_att/q_tot
9 0.00000 0.00361 0.00000 0.00361 0.00000 0.00361 1.00000
15 0.00000 0.00000 0.00000 0.00000 0.00265 0.00265 0.00000
Domain V_satt V_matt V_watt V_att V_rep V_tot V_att/V_tot
9 0.00000 0.02346 0.00000 0.02346 0.00000 0.02346 1.00000
15 0.00000 0.00000 0.00000 0.00000 0.03054 0.03054 0.00000
Domain q_att/V_att q_rep/V_rep q_tot/V_tot sl2r_low sl2r_high
9 0.02281 NaN 0.02281 -0.02543 -0.01939
15 NaN 0.01288 0.01288 0.00946 0.01536
现在对比一下丙二醇和乙二醇的IRI描述符。可以看到丙二醇的氢键域V_att=0.02346 Å^3比乙二醇对应氢键作用区域的V_att=0.01959 Å^3体积略大一点,而前者的电子数0.00361则远多于后者的0.00189,因而前者的平均密度q_att/V_att=0.02281 a.u.比后者的0.01433 a.u.大得多,明显丙二醇的内氢键更强。
再对比一下互斥部分。丙二醇15号域体现互斥部分的平均电子密度q_rep/V_rep=0.01288 a.u.比乙二醇11号域的q_rep/V_rep=0.01577 a.u.更小,清晰地体现出丙二醇的分子内位阻作用更弱。这一点和IRI等值面图相对应,明显可以看到丙二醇中体现互斥的等值面的橙色没有乙二醇的那么深。
可见,IRI描述符充分体现出由于结构特征差异,丙二醇的内氢键比乙二醇的更强的同时,分子内位阻效应还更弱。
7 例:6-螺烯
http://sobereva.com/attach/786/file.zip里的6-helicene.fchk是B3LYP/6-31G*级别下优化的6-螺烯。其标准IRI图如下,可见螺烯中央有一块弱相互作用等值面,既有描述分子内色散作用的绿色区域,也有描述弱位阻作用的棕色区域。
如果你的体系比较大因而盒子不太小的话,格点数会很多,导致高耗时。像本例这种情况,感兴趣的相互作用只出现在体系很小的局部,因此为了节约格点,在设定格点的界面里选10 Set box of grid data visually using a GUI window,通过界面右侧滑杆让蓝色的盒子范围大致能框住感兴趣的等值面出现区域,并且在右下角把格点间距设为0.06 Bohr,如下所示,这样是最划算的盒子。建议阅读《Multiwfn FAQ》(http://sobereva.com/452)中的Q39以及Multiwfn手册3.6节以对格点设置有更充分了解。之后点右上角的RETURN按钮退出此界面,Multiwfn就会开始计算格点数据并产生域了。
在后处理菜单对螺烯中央的那个IRI域计算IRI描述符,结果如下。可见att和rep部分电子数比例为0.02842/0.01645=1.73,比较符合化学直觉,且这两个子域里平均电子密度差不多。
q_satt q_matt q_watt q_att q_rep q_tot q_att/q_tot
0.00000 0.00000 0.02842 0.02842 0.01645 0.04488 0.63338
V_satt V_matt V_watt V_att V_rep V_tot V_att/V_tot
0.00000 0.00000 0.45470 0.45470 0.24857 0.70328 0.64655
q_att/V_att q_rep/V_rep q_tot/V_tot sl2r_low sl2r_high
0.00926 0.00981 0.00946 -0.01244 0.01220
8 例:SN2反应过程
在IRI原文里专门给出了用IRI考察典型的SN2反应[OH]- + CH3CH2Br -> CH3CH2OH + Br-的图,如下所示,横坐标是IRC上的反应坐标,路径上一些有代表性的结构给出了IRI=0.95 a.u.的等值面图。
为了进一步展现IRI描述符的价值、从定量角度理解SN2反应过程中O...C作用的变化以及C...Br作用的变化,这一节演示一下对整个IRC路径上各个点计算IRI描述符,然后分析它们是如何随着反应过程变化的。http://sobereva.com/attach/786/file.zip当中SN2目录下有SN2_0001.wfn到SN2_0080.wfn共80个wfn文件,对应80个IRC点,在Gaussian跑完IRC任务后按照《产生Gaussian的IRC和SCAN任务每个点的波函数文件的工具:IRCsplit和SCANsplit》(http://sobereva.com/199)介绍的做法即可得到它们。为了自动对每个wfn文件计算IRI描述符,参考《详谈Multiwfn的命令行方式运行和批量运行的方法》(http://sobereva.com/612)里介绍的方法,自己写一个Bash shell脚本,内容如下
#!/bin/bash
for inf in *.wfn
do
echo doing $inf..
Multiwfn $inf << EOF > /dev/null
200
14
1
-10
1
3
6b
a 5,8
0
0
q
EOF
done
把这个脚本放到那些wfn文件所在的目录下并执行,就会依次对各个wfn文件用Multiwfn按照指定的指令进行计算,并在当前目录下产生和输入文件同名但后缀是txt的文件。脚本运行原理只要大家认真读过http://sobereva.com/612都能明白,这里仅对传入Multiwfn的关键命令解释一下。6b和之前例子用的选项6的用处是相同的,只不过6b不是在屏幕上显示IRI描述符而是往文本文件里输出它们,让你输入文件名时直接按回车(对应a 5,8下面那个空行)就代表输出的文件名是把输入文件名的后缀替换为txt。a 5,8的意思是把几何中心距离O5和C8原子的中点最近的域作为这次被积分的域。通过这个设置可以确保被积分的就是我们感兴趣的那个域。由于不同IRC点对应O5-C8的域的序号往往不一致、序号不是固定的,所以必须得靠这个设置实现目的。
脚本跑完后当前目录下就有了SN2_0001.txt到SN2_0080.txt共80个txt文件了,用古法手写数据提取脚本比较费时间,如今最省事的做法莫过于直接让AI agent干这种杂事,比如我使用Google的Antigravity挂Gemini 3.8 Flash就可以精准地提取数据。你就告诉agent,这个目录里有一批txt文件对应不同体系,要从中提取q_att、q_rep、q_tot、V_att、V_rep、V_tot、q_att/V_att、q_tot/V_tot属性,要求提取结果中每一行对应一个体系,每一列对应一个属性。agent很快就提取完了,然后放到Origin等程序里作图即可。
先来看看SN2反应过程中从弱相互作用开始逐渐变成化学键的O...C作用对应的IRI域的描述符变化,应结合前面的IRI等值面图来理解。可见黑线所示的q_att从一开始是整体不断上升的,说明OH-和CH3CH2Br的弱吸引作用不断增强,然后再迅速上升,说明从弱相互作用很快变成了共价键这种强相互作用。在键形成后就不再发生改变了,曲线微小的抖动来自于数值积分噪音。红线所示的q_rep从反应一开始不断下降,主要在于体系结构逐渐改变令位阻效应有所削弱。到了反应坐标大约1.4 √amu*Bohr时,q_rep突然降到0,这可以从这个点的IRI等值面图理解,在下图中给出了。可见到这个点的时候,对应O...C作用的等值面和上方的等值面刚好分离开了,使得O...C作用的等值面里不再体现任何互斥作用,只表现O-C吸引作用了。有趣的是,经过这个点后q_att的斜率显著变大、增速明显加快,可以把这个点视为是O...C从弱相互作用变化成化学键的临界点!
下面左图是V_att、V_rep、V_tot的变化,可见O...C作用的IRI域的体积随反应进行迅速减小,体现出O...C相互作用涉及的电子的空间分布从弥散很快变得局域。下图右边的q_att/V_att曲线随反应进行呈现出指数型的增长,体现用于形成O-C共价键的电子密度迅速累积起来,直到反应坐标大约5 √amu*Bohr时达到稳定,此时化学键形成完毕。值得一提的是,大家如果看前面给出的IRC能量变化曲线,会发现也是在大约5 √amu*Bohr处开始能量下降的速度显著减缓,说明令能量迅速下降的最主要驱动力,即O-C化学键的形成过程,就此基本完毕。
在反应过程中C-Br从化学键逐渐削弱成CH3CH2OH和Br-间的弱相互作用。下图是C-Br键的q_att、q_rep、q_tot变化,可见随着这旧键的断裂,描述吸引作用的电子数q_att先迅速下降,下降到低位后进一步的下降显得十分缓和,说明已经到了弱相互作用的程度。起初q_rep为0,体现没有互斥作用,但随着Br-的脱离,中途一度出现了位阻作用(从之前给的IRI等值面图上也能看到红色区域),最后位阻作用又消失了,只剩下弱吸引作用。
C-Br键的V_tot是整体上升的、q_att/V_att整体是指数型下降的,和O-C的变化趋势基本上反过来,这里就不展示和细说了。
9 总结&其它
本文介绍了笔者提出的IRI描述符,这把原本用于图形化定性分析的IRI方法扩展到了定量分析范畴,便于准确分析和横向对比不同的相互作用。本文给出的一系列例子充分展现出IRI描述符能够把相互作用的相对强度、类型、互斥和吸引的成分等方面进行很好的说明。本文还展现出IRI描述符能够用于化学反应等几何结构变化的过程中来对电子结构的改变进行细致的探究。这种分析手段比起流行的AIM拓扑分析仅考察临界点的属性能更全面、充分阐释相互作用情况。
在IRI原文和《使用IRI方法图形化考察化学体系中的化学键和弱相互作用》(http://sobereva.com/598)中还介绍了IRI的变体,专门用于考察pi作用的IRI-pi。对IRI-pi等值面包围的空间同样可以用本文的方法进行积分,用来定量对比不同pi作用的特征,这可以称为IRI-pi描述符,限于篇幅本文就不再示例和讨论了。用域分析模块积分之前把pi以外的轨道占据数清零即可积分IRI-pi域,《在Multiwfn中单独考察pi电子结构特征》(http://sobereva.com/432)里面说了怎么清零。
顺带一提,为了灵活性,Multiwfn的域分析功能的后处理菜单中还提供了选项2b,这里面你可以要求对特定的域内sign(λ2)ρ处于特定范围的区域积分特定的实空间函数,比如你可以要求积分-0.05 < sign(λ2)ρ < -0.015 a.u.区域内的能量密度。这使得用户可以根据实际需要定义额外的描述符。