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

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

《使用Multiwfn图形化研究弱相互作用.docx》由会员分享，可在线阅读，更多相关《使用Multiwfn图形化研究弱相互作用.docx（36页珍藏版）》请在冰点文库上搜索。

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

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

杨伟涛课题组在名为RevealingNoncovalentInteractions（JACS,132,6498-6506）一文中提出了一种新的可视化研究弱相互作用的方法，概念简单清晰，具有广泛意义，很有实用价值，本文目的之一在于介绍、推广这一方法，但与原文的角度有一些不同。

使用这一方法需要计算空间上各点的RDG函数和sign（λ_2）ρ函数的值，虽然可以使用杨伟涛课题组开发的免费的NCIplot程序（http:

//www.chem.duke.edu/~yang/Software/softwareNCI.html）进行计算，但是使用起来有诸多不便。

在Multiwfn1.5程序中不仅可以实现NCIplot的所有功能，还做了很大改进，使这一分析方法更容易投入实际应用，本文目的之二就是结合实例介绍Multiwfn做这种分析时的操作方法。

这种分析方法的研究对象从距离上讲是中、长程相互作用，从类型上讲主要包括氢键、静电、范德华作用，也能展现位阻作用。

虽然原文作者称这种研究方法的对象是非共价相互作用，笔者认为称为弱相互作用更为妥当，文中将使用这种称呼。

这种分析方法也存在一些局限、弊端、随意性，这将在日后另文讨论，本文只以正面角度来讨论。

本文所用的Multiwfn是一个免费、开源、方便、灵活的波函数分析程序，可以从下载最新版本，文中例子中的操作步骤是针对1.5版而言的，若以后版本中选项有所改动，根据程序中选项的含义就能明白怎么操作。

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

此方法的主要目的在于凸显出体系中涉及弱相互作用的区域，以便直观地了解到分子中哪些区域与弱相互作用有关。

也就是定义一个实空间函数，使其数值能够区分开体系中具有不同特征的区域。

此方法使用的是约化密度梯度函数（Reduceddensitygradient,下文称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）。

写一个甲烷二聚体的输入文件，routesection写上#poptmp2/aug-cc-pVDZdensityout=wfn，坐标后面空一行写上.wfn文件的输出路径，附件里的methanedimer.gjf是已写好的。

用Gaussian03执行后得到methanedimer.wfn。

我建议在Gaussian输入文件中的TitleSection部分写上与routesection一致的内容（当然不要把#也写进去），这些内容将会自动传递到.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//绘制“函数1vs.函数2”散点图并生成相应格点文件

1,13//输入函数1和函数2的序号，分别作为散点图的横轴和纵轴。

在Multiwfn支持的函数中ρ（r）是第1号，RDG函数是第13号。

2//用中等质量的网格，总共约512000个点，x,y,z方向的具体点数通过使x,y,z方向格点间距相等来自动确定。

网格的区域自动往分子外延展6bohr。

现在开始计算格点数据。

格点数越多、体系所含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坐标的最大值和最小值往外延展6bohr，留出一定余地，避免等值面被截断。

不过，对于通过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//延展距离为0bohr，即网格范围紧贴着体系。

此时会看到功能-10条目上显示的current:

变为了0。

2//用中等质量格点

4//设等值面数值

0.6//等值面数值设为0.6

-1//观看等值面

此时图像显示出来，如图5右侧所示，点击Showdatarange复选框可以用蓝线显示格点数据涵盖的区域。

点Return关闭图像后，选功能2，高斯格点文件就会被输出到当前目录下的RDG.cub中。

由于总格点数没变，但涵盖的空间范围减小了，所以数据点更密，对弱相互作用区域描述得更精确，等值面的窟窿都不见了，好看了许多，很直观地表现出两个苯之间的π-π相互作用区域。

如果点击界面右侧的Showdatarange，会用蓝框将网格包含的范围显示出来。

虽然苯分子之间相互作用好看了，但是由于网格没有延展，导致苯环中间的体现位阻效应（见后文）的梭形的等值面被截断了一半。

此例之所以观看的不是0.5的等值面，是因为0.5的等值面上也有窟窿，将等值面数值加大可以使等值面范围扩张，补上窟窿，使图像好看。

当然，绝不意味着有窟窿就说明是格点不够精细所致，因为等值面随数值由小到大的变化过程是：

一堆点->一堆小等值面->带窟窿的大等值面->没窟窿的面，如果等值面数值取得较小，必然带窟窿。

用Multiwfn绘制此体系的对称平面上的RDG函数填色图将易于理解这一点。

在Multiwfn里输入以下命令即可绘制。

为得到完整的图，先把RDG_maxrho设为0.0。

benzenedimer.wfn

4//功能4，绘制平面图

-10//设定延展距离。

默认延展4.5bohr，对于RDG函数偏大了

2//改为只延展2bohr

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//延展距离为3bohr

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，对应于在环、笼中出现的较强的位阻效应区域（也被称为nonbondedoverlap），产生张力，因而红色等值面周围原子间起互斥效应。

图9是甲酸二聚体的sign（λ_2）ρvs.RDG的散点图和填色等值面图。

如果将这个散点图的左边折叠到右边去，就还原为ρ（r）vs.RDG图。

[图9]

散点图左边的spike的sign（λ_2）ρ很负，对应很强的氢键，因而相应的等值面为蓝色。

这个spike在散点图上看起来像是一条一条地有规律地组成的，并不致密，这是因为落在这个空间区域的格点偏少，由于格点是均匀、规则地以立方形式排列的，所以函数值变化起来比较有规律，如果增加这个区域格点密度，这个spike会更为致密。

右侧的spike的sign（λ_2）ρ为较小正值，对应于图中棕色圆片等值面，体现了微弱的位阻效应。

比较下面的例子会看到，这种靠弱相互作用结合的复合物，即便之间有位阻效应出现也不会太强，否则将不足以被弱相互作用抵消掉。

至于只靠范德华这种很弱作用力结合的复合物，在平衡状态下不会有位阻区域产生，除非是很强的范德华作用，如π-π堆叠，才可能有微弱的对应位阻效应的区域出现。

下面的例子是邻氯苯酚

[图10]

在苯环中间有红色梭形区域，体现较强位阻效应，对应散点图最右边的spike。

羟基与氯原子之间RDG等值面一小半橘红色，一大半绿色，在散点图上分别对应着spike尖端x值约为+0.06和-0.017的spike。

如果横坐标不是sign（λ_2）ρ而只是ρ（r），则这两个spike是合并在一起的，无法区分究竟代表什么类型的作用，而引入sign（λ_2）使其本质一目了然。

这个等值面说明羟基与氯原子间既存在着位阻效应，也存在着弱氢键作用，互斥和吸引效应并存。

如果做AIM分析，会发现橘红色区域里面是一个（3,+1）临界点，绿色区域里面是一个（3,-1）临界点，两个临界点扩展后连成一个等值面，但各自区域的特征仍然能够靠颜色分辨。

估计会有人存在疑问，羟基与氯原子之间有一大半区域是绿色，从色彩刻度条上看应该对应范德华作用，为何说是弱氢键？

一方面，从散点图上看，范德华作用的spike尖端的ρ（r）不会达到这么大，在原文中作者建议将是否ρ（r）小于0.005作为相互作用是否属于范德华作用的评判标准。

另外O-H----Cl这样的构成也符合形成氢键的条件。

这种情况实际上理应显示成淡蓝色，但由于色彩刻度上下限的设定有一定随意性，导致同一个色彩刻度范围未必对每个体系都很合适。

如果不同体系不用同一个上下限，在横向比较作用强弱的时候又会缺乏基准，而且需要每次手动调整，略微麻烦，一般来说还是按照原文，色彩刻度统一使用[-0.04,0.02]范围为好。

但遇到颜色与期望的差异较大时，最好还是看看散点图上相应spike的位置来确认。

色彩刻度的随意性是这种分析方法的一个弊端，不同色彩刻度下得到的此体系的等值面如下图所示，差异是很明显的。

另外屏幕对比度、可视化程序中分子与光源的相对朝向等诸多问题都可能影响色彩，给定量比较带来些麻烦。

[图11]

下图的四个体系分别是：

1.环氧乙烷与氟化氯通过卤键形成的复合物2.二环[2,2,1]庚烷3.gauche构象的苯乙胺4.直立构象的甲基环己烷。

这些留给读者自行分析，在原文的补充材料中也有很多例子，值得一看。

[图12]

4实例

这里以苯酚二聚体为例介绍如何绘制sign（λ_2）ρ函数填色的RDG等值面图。

phenoldimer.wfn//文件名

100//功能100

1//绘制“函数1vs.函数2”散点图并生成相应格点文件

15,13//sign（λ_2）ρ是第15号函数，RDG是第13号函数

-10//设定延展距离

0//延展距离为0bohr

2//用中等质量格点

现在开始计算，由于体系稍大，所以计算稍微耗时。

由于计算sign（λ_2）需要计算电子密度的Hessian矩阵，会比计算RDG函数要费时不少。

经过约5~10分钟（视CPU速度而定）计算完毕。

当选取的函数1和函数2分别为sign（λ_2）ρ和RDG函数时，散点图的x,y坐标轴范围会分别自动调到[-0.05,0.05]和[0.0,2.0]，所以直接用功能-1就能看到合适的散点图了，如下图上方所示。

之后选择功能3将函数1和函数2的高斯类型格点文件输出到当前目录下func1.cub和func2.cub。

[图13]

虽然能显示高斯类型格点文件的等值面的程序很多，但支持将一个函数数值用不同颜色填到另外一个函数的等值面上的可视化程序比较有限，常用的GaussView虽然支持，但操作不便而且不够强大。

VMD是观看、分析分子动力学结果最重要的软件之一，它在映射颜色和显示等值面方面也很好用，等值面又光滑又有光泽，填色的色彩变化细腻，调整等值面也比较容易，而且运行流畅。

VMD可以免费在http:

//www.ks.uiuc.edu/Developme...cgi?

PackageName=VMD下载。

首先安装VMD，然后将func1.cub和func2.cub复制到VMD安装后的目录下，即vmd.exe所在路径。

然后在此目录下编写一个文本文件，后缀为.vmd，比如ltwd.vmd。

在里面填上如下内容：

molnewfunc1.cub

moladdfilefunc2.cub

moldelrep0top

molrepresentationCPK1.00.318.016.0

moladdreptop

molrepresentationIsosurface0.5000010011

molcolorVolume0

moladdreptop

molscaleminmaxtop1-0.040.02

colorscalemethodBGR

保存文件后，启动VMD，选File-LoadState，选择ltwd.vmd，图13下方的图就显示来了。

若想把背景改成黑色，选Graphics-Colors-Display-Background-8White。

图中的弯曲的片状等值面边缘略有锯齿，可以通过增加此处格点密度来改善。

从颜色上可看出，弯曲的片状等值面描述的是二聚体之间范德华作用，但部分区域也有微弱的位阻效应。

圆片等值面只有中间呈蓝色，说明H与O之间形成了氢键，但并不像甲酸二聚体的氢键那么强。

VMD里面每个操作对应一条语句，载入.vmd本质上就是让.vmd文件中的语句全部执行，这就免得每次都手动执行载入文件、设定参数的一系列繁琐步骤。

比如molnewfunc1.cub的含义就是读入当前路径下func1.cub文件（默认的当前路径一般就是vmd.exe所在路径），molscaleminmaxtop1-0.040.02代表将1号representation（对应于显示填色等值面的那个图层）色彩刻度的下上限分别设为-0.04和0.02，colorscalemethodBGR代表将色彩刻度由小到大设为蓝->绿->红。

将背景设为白色的命令是colorDisplayBackgroundwhite，若将此命令添加到.vmd里面就能在载入.vmd文件时顺便执行，使背景自动设为白色。

这些命令在VMD手册上都有解释。

这些命令也可以在VMD的控制台下直接运行，控制台通过Extensions-TkConsole进入，比如想把色彩刻度改为从-0.05到0.05，就在控制台执行molscal