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是甲酸二聚体的RDG vs. sign(λ_2)ρ的散点图和填色等值面图。如果将这个散点图的左边折叠到右边去，就还原为RDG vs. ρ(r)图。
[图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]
这里以苯酚二聚体为例介绍如何绘制sign(λ_2)ρ函数填色的RDG等值面图。
phenoldimer.wfn //文件名
100 //功能100
1 //绘制“函数1 vs. 函数2”散点图并生成相应格点文件
15,13 //sign(λ_2)ρ是第15号函数，RDG是第13号函数
-10 //设定延展距离
0 //延展距离为0 bohr
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]
虽然能显示高斯格点文件的等值面的程序很多，但支持将一个函数数值用不同颜色填到另外一个函数的等值面上的可视化程序比较有限，Multiwfn目前不支持，常用的Chemcraft也只支持投影到范德华表面上，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范围内改变等值面。
这种分析方法也可以用于分析晶体间内部弱相互作用。虽然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)ρ函数，只需在选择要计算的函数的界面中选择后面带着"with Promolecule approximation"的相应函数即可。注意此时体系中只能存在前三周期的原子。由于不涉及到波函数信息，只要将分子结构传递给Multiwfn就可以，所以既可以通过.wfn、.fch文件将分子结构传入Multiwfn，也可以通过Multiwfn支持的.pdb文件导入分子结构。在Promolecule近似下苯酚二聚体的散点图和等值面图如下所示
[图15]
与B3LYP/6-31G*波函数下的散点图相对比，在Promolecule近似下，主要区别是对应位阻效应的spike明显靠右，这说明芳香环中间的电子密度在形成分子后显著降低以减小互斥效应。并且Promolecule的散点图整体分布比较靠下。所以在观看Promolecule近似下的等值面就必须用与之前不同的等值面数值和电子密度屏蔽范围，如果仍只保留sign(λ_2)ρ在[-0.05,0.05]区域，则最右边的spike将会被截掉很大部分；如果观看的等值面仍为0.5，则y=0.5的直线在散点图中将与一大片范围相交。所以建议保留sign(λ_2)ρ在[-0.1,0.1]区域的点，并观看RDG=0.25的等值面，这样只有这四个spike对应的区域被保留了下来，并且y=0.25的直线只与它们相交。这样得到的填色等值面如图15下方所示，色彩刻度范围设为了[-0.04,0.03]，可见和B3LYP/6-31G*下的图没有太大差别，只是芳香环中间的梭形区域变成椭圆。
虽然靠实际电子密度分析时需要保留的sign(λ_2)ρ的范围有比较通用的值[-0.05,0.05]，但是在Promolecule近似下情况比较复杂，spike的位置变化范围较大，应保留的范围不好确定，所以最好每次都考察一下散点图确定保留范围。在Multiwfn中默认的保留范围是[-0.1,0.1]，可以通过修改settings.ini里的RDGprodens_maxrho参数来改变，当sign(λ_2)ρ小于-RDGprodens_maxrho或者大于RDGprodens_maxrho的点的RDG值将被自动设为100。若此参数设为0.0，则不自动作屏蔽。
[
Last edited by sobereva on 2010-10-2 at 21:58 ]
@通知他人 
引用回复
RSS订阅 ]
No comments:
Post a Comment