使用Multiwfn图形化研究弱相互作用

使用Multiwfn图形化研究弱相互作用

杨伟涛课题组在名为Revealing Noncovalent Interactions (JACS,132,6498-6506)一文中提出了一种新的可视化研究弱相互作用的方法，概念简单清晰，具有广泛意义，很有实用价值，本文目的之一在于介绍、推广这一方法，但与原文的角度有一些不同。使用这一方法需要计算空间上各点的RDG函数和sign(λ_2)ρ函数的值，虽然可以使用杨伟涛课题组开发的免费的NCIplot程序(http://www.chem.duke.edu/~yang/Software/softwareNCI.html)进行计算，但是使用起来有诸多不便。在Multiwfn 1.5程序中不仅可以实现NCIplot的所有功能，还做了很大改进，使这一分析方法更容易投入实际应用，本文目的之二就是结合实例介绍Multiwfn做这种分析时的操作方法。

这种分析方法的研究对象从距离上讲是中、长程相互作用，从类型上讲主要包括氢键、静电、范德华作用，也能展现位阻作用。虽然原文作者称这种研究方法的对象是非共价相互作用，笔者认为称为弱相互作用更为妥当，文中将使用这种称呼。这种分析方法也存在一些局限、弊端、随意性，这将在日后另文讨论，本文只以正面角度来讨论。本文所用的Multiwfn是一个免费、开源、方便、灵活的波函数分析程序，可以从http://multiwfn.codeplex.com下载最新版本，文中例子中的操作步骤是针对1.5版而言的，若以后版本中选项有所改动，根据程序中选项的含义就能明白怎么操作。

1. 用RDG函数等值面显示弱相互作用区域

此方法的主要目的在于凸显出体系中涉及弱相互作用的区域，以便直观地了解到分子中哪些区域与弱相互作用有关。也就是定义一个实空间函数，使其数值能够区分开体系中具有不同特征的区域。此方法使用的是约化密度梯度函数(Reduced density gradient, 下文称RDG)，其表达式为1/(2*(3*π^2)^(1/3))*|▽ρ(r)|/ρ(r)^(4/3)，其中▽是梯度算符，|▽ρ(r)|即电子密度梯度的模。实际上前面的常数项（其值为0.16162）可以忽略掉，只考察|▽ρ(r)|/ρ(r)^(4/3)部分。一个分子体系主要由以下四部分构成： [表1]

表格中的大、小的标准比较模糊，只是定性的。如果我们想要找弱相互作用区域，利用RDG函数的数值大小差异就可以将“原子核附近”和“分子边缘”区域去掉，但“弱相互作用区域”和“化学键附近”的RDG函数值、|▽ρ(r)|值都比较小，区分不开，但ρ(r)存在一定差异。所以，结合使用RDG函数和ρ(r)函数，就可以确定分子中哪些区域涉及弱相互作用。

如果设定立方网格，使网格中的点能够覆盖整个体系，做这些点上ρ(r) vs. RDG的散点图，就可以把上述概念图形化且定量地表述出来。下面来做甲烷二聚体的这种图。

首先要利用Gaussian生成体系的波函数文件(.wfn)。写一个甲烷二聚体的输入文件，route section写上#p opt mp2/aug-cc-pVDZ density out=wfn，坐标后面空一行写上.wfn文件的输出路径，附件里的methanedimer.gjf是已写好

的。用Gaussian03执行后得到methanedimer.wfn。

我建议在Gaussian输入文件中的Title Section部分写上与route section一致的内容（当然不要把#也写进去），这些内容将会自动传递到.wfn文件中的第一行，以免以后忘记.wfn是在什么条件下生成的。也可以不生成.wfn文件，因为Multiwfn可以读入.fch文件，所得结果将是一样的。要注意.wfn文件记录的波函数的基组角动量最高只支持到f函数（因为更高角动量函数在.wfn文件中没有标准的定义），如果用涉及到g角动量的基组，比如cc-pVQZ，就只能通过.fch文件把波函数信息传递给Multiwfn。

先把multiwfn.exe目录下的Settings.ini里的RDG_maxrho设为0.0（注意等号后面要留空格），其原因后文会解释。然后启动Multiwfn，依次输入： c:methanedimer.wfn //输入文件的路径

100 //功能100，其中包含Multiwfn中比较杂的功能 1 //绘制“函数1 vs. 函数2”散点图并生成相应格点文件

1,13 //输入函数1和函数2的序号，分别作为散点图的横轴和纵轴。在Multiwfn支持的函数中ρ(r)是第1号，RDG函数是第13号。

2 //用中等质量的网格，总共约512000个点，x,y,z方向的具体点数通过使x,y,z方向格点间距相等来自动确定。网格的区域自动往分子外延展6 bohr。

现在开始计算格点数据。格点数越多、体系所含Gauss函数越多，计算速度越慢，计算时间与二者都成正比。计算完毕后，输入

4 //设定散点图X轴 0,0.35 //X轴上下限的值 5 //设定散点图Y轴 0,2 //Y轴上下限的值 -1 //绘制散点图

很快ρ(r) vs. RDG的散点图就弹出来了，如图2左图所示。在图上点右键可以关闭，然后选功能1可以将图保存到当前目录下（即Multiwfn.exe所在目录下，后同）。如果对图的效果不满意，可以选功能2将数据导出到当前目录下output.txt，然后用origin、sigmaplot等程序做散点图。其中前三列是数据点的坐标，后两列是两个函数的值。

[图2]

按上述步骤绘制甲烷孤立状态的散点图得到图2右图（设定网格时选3，用高质量网格）。从图2可以看出，体系中存在与不存在弱相互作用时散点图最主要区别在于图中最左侧是否有一个竖条，在原文中这被称为spike。这个竖条上的点正是弱相互作用区域“RDG数值为0~中，ρ(r)^(4/3)数值为小”所对应的点；在右侧也有一个区域RDG 值接近0，这是C-H键区域“RDG数值为0~较小，ρ(r)^(4/3)数值为中”所对应的格点；图中坐标轴范围的更右侧就是原子核附近区域的点；图中左上角的尖峰再往上继续延伸就是分子边缘的区域，虽然离分子越远的地方|▽ρ(r)|和ρ(r)^(4/3)都越小，但后者比前者减小得更快，所以离分子越远RDG值越大，并直至无限大，可自行调整坐标轴观看。不同体系的散点图的成键、原子核附近、分子边缘区域都是类似的，一个体系中是否含有弱相互作用，就是看在ρ(r)较小区域是否有spike出现，这是此分析方法的要点。当然，网格不能太稀疏，如果在弱相互作用区域恰好没有点，spike也不会出现。

接下来，要用等值面确定这些对应于弱相互作用区域的点在实空间上的位置。在计算完格点后的那个菜单中，输入7，然后输入想看的等值面就可以观看第2个函数（即RDG函数）的等值面。其0.5的等值面如图3左图所示

[图3]

图上在两个甲烷中间出现了封闭的等值面，描绘的正是二者间的范德华作用，但是在分子上也出现了三角形封闭的等值面。这是因为成键区域和弱相互作用区域的RDG函数数值范围有很大程度的重叠，如果在图2左图上做一个y=0.5的直线，会发现这条直线不仅贯穿spike，还贯穿成键区域，所以相应的RDG=0.5的封闭等值面不止一个，而在成键区域附近也会出现。这也是前面所说，必须再通过ρ(r)函数区分开成键和弱相互作用区域。

屏蔽掉成键区域，也就是将ρ(r)值稍微大一些的区域，比如ρ(r)=0.05以上的RDG函数的数值设得很大，比如设成100，这样在散点图上y=0.5的直线就不会经过那个区域了，等值面也就只剩下弱相互作用区域。具体做法是在之前的菜单中选择-3，然后输入0,0.05，再输入100，这就表明将ρ(r)的范围在[0,0.05]以外区域的点的RDG函数数值设为100。重新绘制散点图，得到了图4的结果，等值面也变为了所期望的图3右图的情况。

[图4]

一般来讲，观看RDG函数一般观看的是0.5的等值面，这没有什么严格的物理意义，只是等值面大小比较适中。由于弱相互作用区域的ρ(r)一般不会越过0.05，散点图上y=0.5的直线在ρ(r)<=0.05的区域内也只与spike相交，所以每次作弱相互作用等值面图时没必要再考察散点图，只需直接将ρ(r)>0.05的点的RDG函数值设为较大数值就行了。由于这个步骤经常要做，为了方便，Multiwfn在settings.ini里面有一个RDG_maxrho参数，凡是涉及到计算RDG函数的功能，只要某个点的ρ(r)大于这个参数，这个点的RDG值就自动被设为100，这个参数默认被设定为0.05。所以用户就不需要再考虑屏蔽掉成键区域了，这已由Multiwfn自动完成。当然，在后文中作完整的散点图时、或者就是想通过RDG函数研究成键区域时，应当关闭这个功能，将RDG_maxrho设为0.0就代表关闭此功能。

2. 合理地设定网格

Multiwfn计算格点时默认将网格根据原子x,y,z坐标的最大值和最小值往外延展6 bohr，留出一定余地，避免等值面被截断。不过，对于通过RDG函数显示弱相互作用区域来说，这显得浪费了，因为弱相互作用区域是在整个体系内侧，这就导致很多格点白白用于描述没用的区域。如果格点质量不够高，作一些弱相互作用等值面还会有麻烦，比如直接用中等质量格点作苯二聚体的弱相互作用区域RDG=0.6的等值面会得到图5左侧结果，可见薄片状等值面千疮百孔，与原文中的图明显不同，这是因为这个区域格点太稀疏，对RDG函数描述得不够精细。

[图5]

如果将原本被浪费掉的格点利用起来，加强对分子内部区域的描述，将得到更好的等值面。下例将绘制苯二聚体弱相互作用等值面，并自定义延展距离。由于此例只对RDG函数感兴趣，用Multiwfn的功能5（计算格点文件并显示等值面）即可，不需要像前例中用功能100中的子功能1同时把ρ(r)也算出来。起动Multiwfn进行如下操作：

benzenedimer.wfn //已在附件中，几何结构来自原文补充材料，为B3LYP/6-31G*波函数 5 //功能5 13 //RDG函数 -10 //设定延展距离

0 //延展距离为0 bohr，即网格范围紧贴着体系。此时会看到功能-10条目上显示的current:变为了0。 2 //用中等质量格点 4 //设等值面数值 0.6 //等值面数值设为0.6 -1 //观看等值面

此时图像显示出来，如图5右侧所示，点击Show data range复选框可以用蓝线显示格点数据涵盖的区域。点Return关闭图像后，选功能2，高斯格点文件就会被输出到当前目录下的RDG.cub中。

由于总格点数没变，但涵盖的空间范围减小了，所以数据点更密，对弱相互作用区域描述得更精确，等值面的窟窿都不见了，好看了许多，很直观地表现出两个苯之间的π-π相互作用区域。如果点击界面右侧的Show data range，会用蓝框将网格包含的范围显示出来。虽然苯分子之间相互作用好看了，但是由于网格没有延展，导致苯环中间的体现位阻效应（见后文）的梭形的等值面被截断了一半。此例之所以观看的不是0.5的等值面，是因为0.5的等值面上也有窟窿，将等值面数值加大可以使等值面范围扩张，补上窟窿，使图像好看。

当然，绝不意味着有窟窿就说明是格点不够精细所致，因为等值面随数值由小到大的变化过程是：一堆点->一堆小等值面->带窟窿的大等值面->没窟窿的面，如果等值面数值取得较小，必然带窟窿。用Multiwfn绘制此体系的对称平面上的RDG函数填色图将易于理解这一点。在Multiwfn里输入以下命令即可绘制。为得到完整的图，先把RDG_maxrho设为0.0。

benzenedimer.wfn 4 //功能4，绘制平面图

-10 //设定延展距离。默认延展4.5 bohr，对于RDG函数偏大了 2 //改为只延展2 bohr 13 //RDG函数 1 //填色图

200,200 //X,Y方向格点数 1 //绘制XY平面 0 //XY平面的Z值为0

图像很快就弹出来了。如下图所示

[图6]

在分子边界以外没有弱相互作用的区域RDG值很大，远超过1，这样区域都显示为白色。图中央的区域正是描述苯二聚体弱相互作用区域扁片等值面的截面，如果等值面的数值不够大，等值面只能包围每个蓝色区域，显然彼此不相连，如果增大到对应绿色区域的值，孤立区域将会相连接，构成整个扁片等值面。如果继续增大到红色区域对应的值，则苯环中间代表位阻效应区域的等值面将与分子相连而无法区分。

由于原子间存在弱相互作用时（严格来说是指能被RDG函数等值面表现出来的作用）它们的距离一般不会太远，所以一般能猜到哪些区域可能有弱相互作用，而且有时人们只对诸多弱相互作用区域中的某个一感兴趣，此时网格只需要覆盖那个小区域即可，即便网格点数较少，由于密度大，所以也能描述得较精确，可以节省计算时间。然而确定网格空间位

置比较麻烦，Multiwfn在网格设置中提供了一个选项方便研究局部弱相互作用。例如图7中苯酚二聚体之间只有一小块区域相接触，若将网格中心设定在1号和14号原子之间，然后向四周延展一定距离，网格就能覆盖弱相互作用区域。

phenoldimer.wfn //已在附件中，几何结构来自原文补充材料，波函数为B3LYP/6-31G* 5 //功能5 13 //RDG函数 -10 //设定延展距离 3 //延展距离为3 bohr

7 //将两个原子的中点作为网格中心 1,14 //两个原子序号分别为1和14

40,40,40 //由于网格范围小，用较少的格点数就够了 4 //设等值面

0.5 //设等值面数值为0.5 -1 //观看等值面

此时得到图7的结果。网格中心也可以通过直接输入坐标来设定。

[图7]

3 判别弱相互作用的强度与类型

这种分析方法不仅可以指出哪里存在弱相互作用，还可以可视化地了解弱相互作用的强度与类型。

弱相互作用强度一般以相互作用能来衡量，但这是一个全局的量，应用到可视化分析中必须通过局域函数（实空间函数）。在AIM理论中，弱相互作用的临界点的ρ(r)是衡量相互作用强度的重要指标之一，其数值和键的强度存在正相关性，

因而也被用来定义键级。实际上，此文的分析方法在某种程度上可以视为AIM方法的扩展，RDG封闭的等值面一般包围着相应的临界点，如果某个弱相互作用在其临界点处ρ(r)较大，由于ρ(r)的连续性，一般在周围区域ρ(r)也会较大。所以，将ρ(r)的数值大小以不同的色彩映射到RDG等值面上，相互作用的强度就一目了然。

ρ(r)只能反映出强度，但类型需要由sign(λ_2)函数来反映，这个函数是电子密度Hessian矩阵的第二大的本征值λ_2的符号，在AIM理论中键临界点的sign(λ_2)=-1，环、笼临界点的sign(λ_2)=+1，在接近临界点的区域其值与临界点处一般相同。可以将sign(λ_2)函数用不同色彩投影到RDG等值面上，用来表现某一个区域的相互作用类型。

若将ρ(r)和sign(λ_2)函数相乘而得的sign(λ_2)ρ函数投影到RDG等值面上，则弱相互作用的位置、强度、类型都能一目了然地显现出来。

原文中色彩刻度被设为蓝->绿->红，色彩刻度一般设为[-0.04,0.02]，对于个别体系为了色彩效果更好，可以进行调整。不同颜色所代表的ρ(r)、λ_2数值以及所对应的相互作用类型可以用这个图来表示

[图8]

蓝色区域ρ(r)大、sign(λ_2)=-1，表现较强、起吸引作用的弱相互作用，符合这个特征的最常见的就是氢键，还包括较强的卤键等作用。当然如果把ρ(r)更大的，即成键区域也算进去，其等值面也是蓝色。绿色区域的ρ(r)很小，说明相互作用强度很弱，范德华作用区域符合这个特征。由于这样区域电子密度很小，λ_2的符号较为不稳定，所以可正可负。红色区域ρ(r)较大、sign(λ_2)=+1，对应于在环、笼中出现的较强的位阻效应区域（也被称为nonbonded overlap），产生张力，因而红色等值面周围原子间起互斥效应。

图9是甲酸二聚体的sign(λ_2)ρ vs. RDG的散点图和填色等值面图。如果将这个散点图的左边折叠到右边去，就还原为ρ(r) vs. RDG图。

[图13]

虽然能显示高斯类型格点文件的等值面的程序很多，但支持将一个函数数值用不同颜色填到另外一个函数的等值面上的

可视化程序比较有限，常用的GaussView虽然支持，但操作不便而且不够强大。VMD是观看、分析分子动力学结果最重要的软件之一，它在映射颜色和显示等值面方面也很好用，等值面又光滑又有光泽，填色的色彩变化细腻，调整等值面也比较容易，而且运行流畅。VMD可以免费在http://www.ks.uiuc.edu/Developme ... cgi?PackageName=VMD下载。

首先安装VMD，然后将func1.cub和func2.cub复制到VMD安装后的目录下，即vmd.exe所在路径。然后在此目录下编写一个文本文件，后缀为.vmd，比如ltwd.vmd。在里面填上如下内容： mol new func1.cub mol addfile func2.cub mol delrep 0 top

mol representation CPK 1.0 0.3 18.0 16.0 mol addrep top

mol representation Isosurface 0.50000 1 0 0 1 1 mol color Volume 0 mol addrep top

mol scaleminmax top 1 -0.04 0.02 color scale method BGR

保存文件后，启动VMD，选File-Load State，选择ltwd.vmd，图13下方的图就显示来了。若想把背景改成黑色，选Graphics-Colors-Display-Background-8 White。图中的弯曲的片状等值面边缘略有锯齿，可以通过增加此处格点密度来改善。从颜色上可看出，弯曲的片状等值面描述的是二聚体之间范德华作用，但部分区域也有微弱的位阻效应。圆片等值面只有中间呈蓝色，说明H与O之间形成了氢键，但并不像甲酸二聚体的氢键那么强。

VMD里面每个操作对应一条语句，载入.vmd本质上就是让.vmd文件中的语句全部执行，这就免得每次都手动执行载入文件、设定参数的一系列繁琐步骤。比如mol new func1.cub的含义就是读入当前路径下func1.cub文件（默认的当前路径一般就是vmd.exe所在路径），mol scaleminmax top 1 -0.04 0.02代表将1号representation（对应于显示填色等值面的那个图层）色彩刻度的下上限分别设为-0.04和0.02，color scale method BGR代表将色彩刻度由小到大设为蓝->绿->红。将背景设为白色的命令是color Display Background white，若将此命令添加到.vmd里面就能在载入.vmd文件时顺便执行，使背景自动设为白色。这些命令在VMD手册上都有解释。这些命令也可以在VMD的控制台下直接运行，控制台通过Extensions-Tk Console进入，比如想把色彩刻度改为从-0.05到0.05，就在控制台执行mol scaleminmax top 1 -0.05 0.05。有很多命令在VMD的GUI上有相应的选项，其中有些通过GUI操作会容易得多，比如调整等值面数值，可以在Graphics-Representations里面的列表中选定第二个显示模式（即Style为isosurface的那个），然后将下方Range的下上限分别设为比如0和1，点回车，之后拉动滑条就能在0~1范围内改变等值面。

另外，Chemcraft也可以实现相同功能，对初学者来说使用更简单，不过效果不如VMD，而且收费。使用方法： 先打开func2.cub，再点左下角Add cube，选择func1.cub。将Contour value设为0.5，敲回车，然后点Show isosurface。然后把Map other:选为2，再把Values range分别填-0.04和0.02，敲回车。

这种分析方法也可以用于分析晶体间内部弱相互作用。虽然Gaussian的PBC功能并不给力，但是简单的PBC计算还是可以胜任的。由于Gaussian的PBC计算是以高斯型函数作为基函数，所以波函数信息可以直接被Multiwfn读入并进行分析。这里将以金刚石晶体作为例子。

还是先获得波函数文件。Gaussian的输入文件就是压缩包里的pbc_diamond.gjf。运行之后，用formchk将.chk转化为.fch文件。当然，用.wfn作为Multiwfn的波函数输入也可以。由于计算的是素晶胞，波函数信息也只含有两个碳原子的，获得的等值面显然体现不出周期性，然而计算复晶胞又太费时。在Multiwfn里提供了一个晶胞波函数平移复制

功能，可以将这素晶胞波函数信息扩展为足够大的复晶胞波函数。复制次数越多，体系中高斯函数越多，计算格点时越慢，为避免计算格点时间太长，这里只向各方向平移复制两次。

启动Multiwfn后输入： c:pbc_diamond.fch 6 //修改波函数功能 14 //平移复制体系

4.7523,0,0 //第一个平移向量。平移向量在输出文件中的PBC vector段落，或者查看.fch中的Translation vectors字段。不要用输入文件中的平移向量，因为Gaussian可能自动改动坐标，平移向量也就变了。 1 //单位为bohr

2 //复制两次。接下来再对另外两个方向做相同的平移复制。 2.3762,4.1156,0 1 2

2.3762,1.3719,3.8802 1 2

0 //平移复制已完毕，保存当前波函数到当前目录下new.wfn，也就是压缩包中的pbc_diamond_dup.wfn。

由于得到的体系是斜着的，把所有原子都纳入立方网格中必将造成很多格点落在体系外而浪费掉。实际上只要取体系中间一个原子（28号）作为网格中心，然后向四周延展一些距离，所得等值面就足够体现周期性了。现在生成格点文件，依次输入

pbc_diamond_dup.wfn 100 1 15,13

7 //用两个原子的中点作为网格中心

28,28 //当输入的两个原子序号相同，则以此原子为中心向四周延展。Multiwfn默认延展距离是6 bohr，对此例子比较合适，不用修改。 80,80,80 //各方向格点数

算完后，用功能3保存格点文件，按照上例方法，用VMD显示结果如下。可见，由于碳原子间距离较近，原子空穴间充满红色等值面，体现很强的位阻效应。

[图14]

5 波函数质量产生的影响

在前面的例子中多数用的是B3LYP/6-31G*波函数，必定有人会认为用B3LYP/6-31G*计算以范德华作用为主的弱相互作用体系很不合理，但实际上相互作用能与理论方法关系更大，其值很差不代表相应的ρ(r)就会差，原文作者认为使用这种分析方法时用B3LYP/6-31G*就够了，这样计算量也小，而且改为MP2/6-311++G**后等值面并没有什么改变。甚至ρ(r)粗糙到只用自由电子密度叠加都能得到定性一致的结果，所以这种分析方法对电子密度的质量很不敏感，使之用于大体系成为了可能。

6 Promolecule近似及对结果产生的影响

使原子坐标保持在形成分子时的状态，将自由原子密度叠加得到的密度称为Promolecule密度，可以视为在形成分子前，电子密度尚未驰豫的电子密度。构建这种密度不需要量子化学计算，只要有分子结构文件和自由原子密度就能十分容易地构建。由于在弱相互作用区域Promolecule密度和实际密度差异并不像在成键区域那么显著，Promolecule密度在此分析方法中可以近似代替实际电子密度，等值面的基本形状不会有太大差异，但是定量上，也即等值面细部特征会有一定差异，因为电子密度驰豫后会在吸引作用区域聚集，尤其是会在位阻效应区域疏散以减小斥力，位阻效应越强，改变量越大。

原子在自由状态的真实电子密度是球对称的，但是对于比如基态氧原子，p轨道一个是双占据两个是单占据，量化算出

来的电子密度不是球对称的，这将导致分析结果出现不应有的取向性，所以需要将之球对称化。球对称化方法不是唯一的，比如可以简单地对三个p轨道占据数取平均，也有人利用GVB方法解决，原文的方法是用s型STO函数拟合B3LYP/6-31G*原子密度，第一、二、三周期的原子分别用1、2、3个STO拟合，以对应壳层数。这不仅解决了球对称问题，还有另一个好处，就是描述原子密度的函数大大减少了，6-31G*描述氧原子用28个GTO，而拟合后只用2个STO，这使得计算RDG函数、sign(λ_2)ρ的速度也大大加快。实际上STO用的少主要在于双电子积分时的困难，由于STO能正确地表现原子轨道波函数随r增大的收敛行为和原点处Cusp的特征，如果研究内容不涉及到双电子积分，比如这种只依赖ρ(r)的分析方法，用STO又便宜又好。

在Multiwfn中使用Promolecule密度计算RDG与sign(λ_2)ρ函数，只需在选择要计算的函数的界面中选择后面带着\的相应函数即可。注意此时体系中只能存在前三周期的原子。由于不涉及到波函数信息，只要将分子结构传递给Multiwfn就可以，所以既可以通过.wfn、.fch文件将分子结构传入Multiwfn，也可以通过Multiwfn支持的.pdb文件导入分子结构。在Promolecule近似下苯酚二聚体的散点图和等值面图如下所示

比RDG等值面用更大的函数数值。n=1时，由于远离分子时ρ(r)的收敛速度不像ρ(r)^(4/3)那样快，所以散点图中对应分子边缘区域的左上角的峰没那么明显。如果继续减小n，则弱相互作用与分子边缘区域就难以区分，数据点都一起缩到图中左下角。RDG函数原本是用在GGA密度泛函中表达电子密度梯度项，这也正是其名字的含义，之所以不是|▽ρ(r)|而要除以ρ(r)^(4/3)，是为了消掉|▽ρ(r)|的量纲，使梯度的模成为无量纲的量。而RDG函数，即n=4/3的情况用在分析弱相互作用上最合适，是属于巧合。

12. 总结

RDG填色等值面分析方法是十分方便、有用的分析方法，原理简单易懂、适用范围广，在Promolecule近似下只需要提供分子坐标就能很快地得到结果，值得大力推广。Multiwfn对这种分析方法投入实际应用提供了极大便利，操作傻瓜化，不涉及到任何复杂的理论、概念。此方法也存在一些不足，也因此可能有进一步改进或发展的余地，比如等值面数值的选取、色彩刻度的上下限设定有一定随意性，那么是否有可能将这种分析定量化？其等值面面积和包围的体积是否也有一定意义？此分析方法的提出纯粹是从可视化便利的角度出发的，是否能找出其背后与ELF/LOL等函数的理论联系，进而组合、衍生出一个能分析更广泛问题的函数？笔者曾在《电子定域性的图形分析》一文中浅谈了实空间函数图形化分析方法的优点，本文的RDG填色等值面分析方法进一步显现了这类分析方法的用处，这类方法在未来必将得到进一步发展，并逐渐流行起来，也希望读者多加实践

本文来源：https://www.bwwdw.com/article/53h6.html