通过独立梯度模性(IGM)考察分子内和分子间的弱相互作用

本文尚未写完,仅供心急的人预览


通过独立梯度模性(IGM)考察分子内和分子间的弱相互作用

文/Sobereva@北京科音  2018-Jan-24


0 前言

早在2010年,笔者写了一篇《使用Multiwfn图形化研究弱相互作用》(http://sobereva.com/68),文中详细介绍了杨伟涛等人通过约化密度梯度(RDG)图形化考察弱相互作用的方法,以及在Multiwfn中的实现。这种分析方法在文献中也常被叫做NCI(Noncovalent interaction)方法。此文贴出来后,不断有大量使用Multiwfn做RDG分析的文章涌现,而RDG分析方法也被很多研究者进一步发展,例如:
(1)aRDG方法,用分析动力学过程中的弱相互作用,见《使用Multiwfn研究分子动力学中的弱相互作用》(http://sobereva.com/186
(2)将RDG与ELF联用,同时考察弱相互作用和共价相互作用,例子见《通过键级曲线和ELF/LOL/RDG等值面动画研究化学反应过程》(http://sobereva.com/200
(3)DORI方法,可以通过一个函数同时展现体系中各种类型的相互作用,见《使用DORI函数同时考察共价和非共价相互作用》(http://sobereva.com/367
(4)对RDG等值面内部空间的属性进行积分,从而能够定量讨论相互作用强度。见Multiwfn手册4.200.14节的例子。

2017年,在Phys. Chem. Chem. Phys., 19, 17928中新提出一种叫做独立梯度模型(Independent Gradient Model, IGM)的方法。这个方法受到RDG方法的很大启发,但是底层思想不同。和RDG一样,IGM也是重在展现弱相互作用区域及其特征,但是有很多优势,简要列举如下,后文会详细说:
(a)给出的信息会明确划分成片段间和片段内两套数据,因此想考察分子间相互作用时,不会被分子内相互作用所干扰。(实际上,用RDG方法时,如果结合Multiwfn手册4.13.4节示例的格点数据屏蔽功能,其实也能达到相同目的)
(b)计算时只依赖于原子坐标,而且计算很快。因为此方法计算时只需要自由原子密度,并不需要体系的波函数信息,而自由原子密度是内置于IGM分析程序当中的。(实际上,RDG也可以基于初分子(promolecular)密度来近似地快速计算)
(c)IGM方法给出的等值面图比RDG更为饱满一些,不像RDG图在格点稀疏的时候容易出现难看的锯齿
(d)可以考察每一对原子间的相互作用程度
(e)可以给出每个原子对片段间相互作用程度,而且可以利用VMD通过填色直观地展现
(f)在分子间相互作用区域,IGM方法里定义的dg函数的大小与相互作用强度有正相关性

2018-Jan-23及之后更新的Multiwfn 3.5版已经全面支持IGM分析,本文的目的一方面是介绍IGM方法,一方面通过实际体系示例怎么通过Multiwfn做IGM分析。值得一提的是,在DOI: 10.1002/cphc.201701325中,作者又提出了基于波函数计算的IGM分析方法,不过基于波函数的版本计算复杂,又没什么额外好处,所以Multiwfn不支持,也不打算做介绍。阅读本文前,强烈建议仔细阅读前述博文,至少弄清楚RDG分析的基本原理以及在Multiwfn中的操作,否则在理解本文内容时将遇到困难。

Multiwfn可以免费在http://sobereva.com/multiwfn下载。如果对Multiwfn不了解,看《Multiwfn入门tips》(http://sobereva.com/167)、《Multiwfn波函数分析程序的意义、功能与用途》(http://sobereva.com/184)。除本文提到的方法外,Multiwfn还支持一大票弱相互作用分析方法,看《Multiwfn支持的弱相互作用的分析方法概览》(http://sobereva.com/252)。本文使用的VMD是1.9.3版,可在http://www.ks.uiuc.edu/Research/vmd/免费下载。


1 IGM方法的原理

IGM的详细介绍看其原文Phys. Chem. Chem. Phys., 19, 17928 (2017),这里主要通过语言解释一下IGM最核心的思想,看懂本文后有兴趣可以再看看原文的数学表达。我们先看一个很简单的情况,氢分子。在氢分子的轴线上,H1和H2的电子密度分布如下所示,这里用的原子的电子密度是原子在孤立状态下计算的

从上图可看出,在两个原子间,两个原子的梯度方向正好是相反的。H1的电子密度的梯度为负(即往左边密度增大),而H2的电子密度的梯度为正(即往右边密度增大)。因此,在原子间,两个原子的密度梯度是相互抵消的,而且在两原子正中央正好精确抵消,导致初分子密度(两个原子自由状态密度叠加构建的密度)的梯度为0,此处也正是与初分子密度相对应的键临界点(BCP)的位置。

一般方式计算的初分子密度的梯度等于各个原子密度梯度的加和。而IGM型密度梯度,则是把每个原子的梯度先取绝对值再加和。二者求差,就得到了所谓的δg函数。对上面的氢分子体系,图示如下

图中黑线是一般方式计算的初分子密度梯度,在原子核和原子连线中点都为0。绿线是IGM型密度梯度,由于计算时对两个原子密度梯度都取了绝对值,所以不会相互抵消,在原子之间其数值也很大。IGM密度梯度是普通方式计算的密度梯度的上限,在相互原子间相互作用区域二者相差最大,即δg最大。δg对应图中深蓝色曲线,可见利用δg这个函数可以明确地把原子间相互作用区域展现出来。而且,原子间相互作用越强,相互作用区域的δg就会最大。这是因为如果原子间相互作用较弱的话,那么优化后的结构两个原子距离就会较远,它们相交的区域都是各自密度梯度已经比较小的部分,此时IGM型密度梯度和普通密度梯度相差也就没那么大了。

之所以当前这个方法叫做IGM(独立梯度模型),是因为常规方式初分子密度梯度时,相当于原子密度梯度之间发生了“干涉”,影响了“相位”。而在计算IGM型密度梯度时,纯粹是原子的密度梯度大小的简单叠加,不存在干涉现象,因此这些原子的密度梯度就彼此“独立”了。

δg是个三维实空间函数,有不同的考察方式,如果图形化展现的话,一般是像RDG那样绘制等值面图。比如对上面的氢分子体系,如果绘制δg=0.2的等值面图,就肯定会看到两个氢之间出现了等值面,因为上图中对应Y=0.2的虚线和δg曲线在黑色箭头标注的位置出现了相交。

δg对应于Multiwfn中第37号自定义函数。如果将settings.ini里的iuserfunc设为37,再用主功能4对第100号函数绘制填色平面图的话,就可以看到下面的效果。下图是GC碱基二聚体平面上的情况

由图可见,所有成键的原子间δg都较大,最大值处都超过了图中色彩刻度上限0.3因此显示为白色。而弱相互作用区域,比如N-H...N的氢键区域,只有一块深蓝色,说明这样的地方δg相对较小,作用较弱,但终究比没有相互作用的地方δg值要大。IGM原文里通过测试还发现,二聚体的结合能,和两个单体间对应弱相互作用的BCP位置上的δg值有很好的正相关性。因此δg是个很有用的函数。

δg可以划分为用于讨论片段内相互作用的δg_intra和讨论片段间相互作用的δg_inter。在说这个之前,这里先把δg的比较严谨的数学定义给出来。其中i是原子序号,▽ρ是梯度矢量,abs(▽ρ)代表对里面▽ρ矢量的每个分量都取绝对值(还是保持矢量形式),| |代表对矢量取模。

IGM原文中计算δg_intra和δg_inter只考虑到两个片段的情况,但笔者将之推广到了多个片段的情况,可表达为:

式中A是片段编号。计算时不要求所有片段的并集必须对应整个体系,但不对应时,在计算δg的时候必须只考虑这些片段涉及的原子。只要想明白了δg是怎么定义的,就很容易理解δg_inter为什么是这样定义的。

由上式可见,δg_intra+δg_inter=δg。δg衡量的是体系中所有原子间相互作用,稍加变通就成了只体现片段间相互作用的δg_inter,再把这部分从δg中扣除就成了描述片段内相互作用的δg_intra。IGM的这个思想相当不错。

对δg_intra和δg_inter分别绘制等值面图,就可以清楚看到分子内和分子间存在的相互作用区域。还可以像RDG分析时那样,把sign(lambda2)rho函数以不同颜色投影到这两个函数的等值面上,从而能够清楚地判断相互作用区域是吸引还是互斥作用,以及强度如何。另外,观看δg_intra等值面的时候还可以像RDG分析那样,把电子密度较大的区域的δg_intra值设为0,这样等值面就只体现出分子内弱相互作用区域了,而化学键区域就都被屏蔽掉了。

还可以定义原子对δg指数,用于要考察相互作用的两个片段间的每对原子的相互作用大小。这个指标在IGM原文里没有提出,是我把IGM实现进Multiwfn的时候顺带想出来的。计算某两个原子间的这个指数,就是把这两个原子各自作为一个片段来计算它们之间的δg_inter函数,然后再对这个函数值在全空间进行积分,数值越大显然说明两个原子间作用越强。然后,可以再计算原子δg指数,就是把这个原子与另一个片段中所有原子的原子对δg指数进行加和。之后,可以在VMD里将每个原子按照原子δg指数用不同颜色显示出来,哪些原子对片段间相互作用起主要贡献就一目了然了。

评论已关闭