| 查看: 10740 | 回复: 30 | |||||||||||||||||||||
| 【奖励】 本帖被评价25次，作者sobereva增加金币 19.8 个 | |||||||||||||||||||||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||||||||||||||||||||
[资源]
使用Multiwfn做电子密度、ELF、静电势、密度差等函数的盆分析
|
|||||||||||||||||||||
|
使用Multiwfn做电子密度、ELF、静电势、密度差等函数的盆分析 前言：在Multiwfn 3.0当中新加入了一个十分重要的功能，即盆分析功能（主功能17）。本文目的在于简要介绍盆分析的概念、算法，并结合几个应用实例介绍盆分析功能的基本操作方法并展现盆分析的基本用途。更多的信息参见Multiwfn手册3.20节和4.17节。 1 盆、吸引子和盆分析的基本概念 盆(basin)这个概念最早是由Bader在它的Atoms in molecules (AIM)理论中引入的，在AIM理论中盆分析是对于电子密度来进行的。随后这个概念被Savin、Silvi推广到用于分析电子定域化函数(ELF)。关于ELF的介绍可参见《电子定域性的图形分析》（http://hi.baidu.com/sobereva/item/4a8aeae2291d8d0f8c3ea8d8）和 物理化学学报,27,2786。实际上，盆这个概念对于任何实空间函数都适用，比如电子密度拉普拉斯函数、动能密度、静电势、电子密度差、福井函数等等。 盆是整个空间中的局部空间，而所有的盆所占空间总和就是整个空间。每个盆彼此间通过函数的零通量面(zero-flux surface)来划分，因而零通量面也叫做盆间面(interbasin surface)。具体来说，零通量面上的各个点上都满足函数的梯度在垂直于这个面方向上的分量为0的条件。每个盆当中唯一地包含一个“吸引子”(attractor)，它是这个盆空间内函数的极大点。吸引子和拓扑分析术语中的(3,-3)型临界点是一回事。 例如，下图是电子密度的盆的二维示意图，这种图经常出现在AIM分析中  棕色的圆点就是吸引子。由于电子密度的极大点总是很接近原子核（对于H可能偏得多一点。另外还有所谓的非核吸引子，但这里不提及），所以这些吸引子的位置基本就在原子核上。灰色细线是电子密度的梯度线，从各个方向向着吸引子汇集。粗的蓝线就是零通量面，将整个分子空间分割为一个个盆。这些盆和原子一一对应，所以AIM理论中将电子密度的盆称为原子盆。这样的盆被AIM理论认为描述了原子在分子中所占的空间，所以也被叫做AIM原子空间。上面的图只是盆在分子平面上的二维示意图，氢原子盆的三维图如下所示  注：虽然众多梯度线组成的面上所有点也都符合“垂直于这个面上函数的梯度分量为0”的条件，但是这样的面穿越吸引子，所以不特称为零通量面或盆间面。 将盆划分好之后，就可以在这些盆的空间里进行各种分析。通常有三类分析，它们都被Multiwfnd的盆分析模块所支持： (1)在盆里对指定的实空间函数进行积分。例如，积分电子密度，就得到了盆内电子的布居数。如果这个盆是AIM盆，那么得到的就是原子的电子布居数，令原子核电荷减去它就得到了AIM原子电荷。如果这个盆是ELF，且比如对应共价键，那么可以通过盆中的电子密度积分值讨论有多少电子与共价键的构成相关。还可以积分诸如自旋密度、能量密度等等。用于定义盆的实空间函数和在盆中被积分的实空间函数显然是可以不同的。Multiwfn可以让任意实空间函数在任意实空间函数定义的盆中被积分。 (2)计算盆内的电多极矩(electric multipole moments)，具体来讲是指单极矩、偶极矩和四极矩。这对于定量分析电子密度的分布十分有益。比如可以研究孤对电子对应的ELF盆对分子偶极矩的贡献（如J. Comp. Chem., 29, 1440）、考察形成复合物过程中电子密度被极化的过程、研究盆内电子分布在哪个方向较弥散、偏离球对称分布程度的多少等等。 (3)计算盆内的定域化指数(localization index, LI)和盆间的离域化指数(delocalization index, DI)。前者衡量平均来说有多少电子定域在盆里；后者衡量的是平均有多少电子离域在两个盆之间，或者说盆之间共享了多少电子。如果是AIM原子盆，那么它们之间的DI和共价键级是可以相对应的。 注意(2)和(3)方式的分析本质上也是需要在盆内进行积分才能实现的。所以，要想做盆分析，必须得通过某种算法将盆空间表示出来，并且能够在其中对特定实空间函数进行积分。这是下一节要讨论的。 2 生成盆的算法 由于AIM理论的重要性，好几十年前就开始有人不断提出各种各样的积分AIM原子盆的方法。这些方法都是解析方法，例如J.Comp.Chem.,21,1040和J.Phys.Chem.A,115,13169。遗憾的是，由于盆的边界复杂，以解析方式积分相当困难，非常耗时，且积分公式都很复杂，不容易实现。这些方法一般需要先搜索到电子密度的(3,-3)及(3,-1)临界点，然后利用(3,-1)临界点产生盆间面，然后以(3,-3)临界点为中心向外发散地做盆内的积分。这类解析方式的积分不仅慢、复杂、通常需要事先找临界点，而且有一个重要缺点就是基本上只能用于积分电子密度的盆，不适合用于其它实空间函数如ELF。这类方法的唯一好处就是积分精度很不错。经典的AIM程序如AIMPAC、AIM2000、AIMALL、Morphy之类用的都是这类解析积分方法。值得一提的是在J.Comp.Chem.,30,1082中作者提出一种利用类似Becke的DFT泛函积分的方法积分盆，虽然积分精度有所降低，但是积分速度快了很多，且算法复杂度降了不少。但是这种方法依然只适合于积分电子密度的盆。 Silvi他们组的专用于ELF盆分析的TopMod程序用的是一种基于立方格点的数值方式的盆积分方法。在此程序发表后，逐渐有一些其它的基于立方格点（实际上也完全可以用矩形格点）的盆积分方法和相应程序被提出。这类基于立方格点的数值积分方法速度比较快，最大的好处是可以用于任意类型的实空间函数，而缺点是积分精度逊于解析积分方法。或者说，用的格点数据的格点间距越小则积分精度越高，但是要想达到解析方法那么高的精度，那么计算格点数据的时间以及储存格点数据所需内存量将甚巨而难以实现。不过，只要对积分精度要求不苛刻，那么这类格点积分的精度一般也够了。值得一提的是这类方法产生盆之前不需要事先寻找临界点，吸引子会在产生盆的过程中自动获得，所以很方便。由于Multiwfn的盆分析模块强调的重点之一是普适性强，所以采用的是立方格点积分的方法，而非解析方式的积分盆。这类积分方法中目前最好的是near-grid方法（J.Phys.:Condens.Matter,21,084204），这也是Multiwfn默认的方法，而它的前身on-grid方法（Comput.Mat.Sci.,36,354）在Multiwfn中也支持。On-grid方法生成盆的速度更快，但明显看得出盆间面的划分准确度不如near-grid方法，因此积分值准确度也差一些。 积分不同的实空间函数对格点精度要求不同。比如电子密度拉普拉斯函数、源函数(source function)在原子核附近波动极大，所以就不适合用立方格点方法来积分。即便用了很高精度的格点，仍会发现在盆内积分拉普拉斯函数的结果偏离理想值0比较多（多级优化的积分方式可以缓解这一问题，见本文4.5节）。另外，不同区域的盆积分难度也不同。例如电子密度在原子核附近有尖峰，而在原子价层区域变化则较缓，所以在核电子对应的ELF盆里积分电子密度的精度就相对低一些，明显不如积分价层ELF盆那样容易得到较高精度。 near-grid方法得到吸引子位置、产生并积分盆的大体过程是： 首先提供要分析的实空间函数的格点数据。然后，依次从每个格点出发，按照梯度最大的方向沿着格点不断进行移动从而产生出一条轨迹。如果走到一个格点且发现它的数值比相邻的点都大，那么这个格点就会被当成吸引子，并且这条轨迹上所有点将被指认(assign)到这个吸引子。如果从某个点开始沿着梯度最大方向不断移动的过程中走到了一个已经被指认了吸引子的点，那么这条轨迹中的所有点也将被指认到这个吸引子。将所有格点循环一遍后，被指认到同一个吸引子的所有点就构成了这个吸引子对应的盆。盆内任意一实空间函数f的积分值可通过此式很容易地得到：∑f(r_i)*dv。其中加和符号循环盆内每个格点，编号以i表示，r_i即是格点的坐标。dv是空间微元，即每个小格子的体积。 3 现有的支持盆分析的程序 那些经典的利用解析方法做电子密度盆分析的AIM程序前面已经提到过了，即AIMPAC、AIM2000、AIMALL、Morphy等，它们没法用于对其它实空间函数进行盆分析。Bader程序（http://theory.cm.utexas.edu/bader/）是极少数基于立方格点方法进行电子密度盆分析的程序，操作还算方便，但是功能很弱，仅能计算AIM原子盆内的电子密度积分值。还有的程序如InteGriTy，虽然载入的是电子密度格点数据，但是实际上还是用传统解析方法进行盆积分，只不过任意处的电子密度不是直接计算而是通过格点数据插值得来。 普适型的基于立方格点方式做盆分析的程序最出名的是TopMod（http://www.lct.jussieu.fr/pagesperso/silvi/silvi_english.html），但是它实际上只支持电子密度和ELF的盆分析，也只能在这两类盆里计算LI/DI和积分电子密度。TopMod的操作超级繁琐，做个盆分析需要调用多个子程序（top_grid、top_bas、top_pop），得敲一大堆命令。而且程序诸多方面都显得莫名其妙，颇不好用，而且想可视化盆也非常困难。Multiwfn的盆分析功能的最初开发目标就是完胜TopMod并彻底取代它在ELF盆分析上的统治地位。TopGrid（http://www.lct.jussieu.fr/pagesperso/pilme/topchempage.html）是与TopMod关系甚密的另一个基于立方格点的盆分析程序，只能读入格点数据来做盆分析，生成盆的算法和TopMod有异，普适性号称比它要广，可以支持盆的电多极矩分析。但这程序甚至比TopMod还难用，而且计算结果笔者发现不准，比如观看它产生的盆会明显发现盆的形状和理想形状往往偏离不小。 相对于上述程序，Multiwfn的盆分析功能无论是普适性、计算速度、灵活性还是操作便利性上都有质的飞跃。Multiwfn可以对任意实空间函数划分盆，对任意实空间函数在盆内进行积分，还可以计算盆的电多极矩、LI和DI。实空间函数既可以直接由Multiwfn计算也可以从外部文件读取。分析结果能够直接可视化而不必再借助于第三方程序。同时，吸引子和盆都可以导出到外部文件，以便使用诸如VMD等可视化程序显示来得到更好的效果。由于Multiwfn本身代码的高效，以及作了充分的并行化，做盆分析速度很快，包括C60这样较大的体系的盆分析也不用花多少时间。Multiwfn盆分析模块操作极其方便，只需要按几个键，敲几下回车即可完成，自动化程度极高。而且输出信息简明易懂，各个功能、涉及到的原理在手册里都能方便地查阅到。可以说，Multiwfn 3.0的发布标志着实空间函数的盆分析进入了新纪元。 PS: 由于前面已经提到，立方格点方式的积分精度比解析方法弱一些，所以对于AIM盆分析，如果对精度要求苛刻（比如达到小数点后三位），那么我目前推荐用AIMALL去做积分。但如果要求不很严格，那么结合相对比较精细的格点，用Multiwfn得到的精度就足矣，这也能省大量的计算时间。 4 实例 以下五个实例涉及到盆分析模块的大部分操作，做过一遍后就能对盆分析有个基本的认识。本文所用Multiwfn为3.0版，一些选项在未来可能有变动，以屏幕上的实际提示为准。//后面的是注释。 4.1 HCN分子的电子密度盆（AIM盆）分析 启动Multiwfn并输入 examples\HCN.wfn 17 //盆分析 1 //生成盆并寻找吸引子 1 //让Multiwfn计算电子密度格点数据，用来得到它的盆和吸引子 2 //中等质量格点，也就是格点间距为0.1Bohr。格点数据的空间范围程序将自动确定，也就是从分子边界起往各个方向延展边界原子的1.7倍范德华半径的距离，如果想调这个倍数的话此时可以选-10（对于阴离子、存在带很大负电荷的原子的时候倍数应调大，以避免分子外围区域没有被完全地积分）。列表中还有其它很多定义格点设定的方式，这里暂时先不用 此时Multiwfn开始计算格点数据，然后生成盆并确定吸引子的位置。之后吸引子的坐标以及吸引子位置上的函数值会显示在屏幕上 Attractor X,Y,Z coordinate (Angstrom) Value 1 0.02047487 0.02047487 -1.54308028 0.34007350 2 0.02047487 0.02047487 -0.48472578 62.82169829 3 0.02047487 0.02047487 0.62654644 77.37068545 由于氢的密度低，很显然吸引子1就是与氢相对应的。如果之后还想重新看吸引子的信息，以及更详细的信息，可以选择选项-3。注意如果开启了Multiwfn的并行模式（默认是开启的），每次寻找出的吸引子的编号顺序都可能不同，以实际为准，后同。 现在选0来可视化刚找出来的吸引子和生成的盆。在弹出来的图形界面中，紫色的文字对应吸引子的标号。每个吸引子用绿色小球表示，但由于它们对于此例在原子核处，所以得取消选择show molecule选项才能看到它们。图形界面中各个控件的含义都很明确就不多提了，试试就明白了。界面右下角（如果是Linux版则在右上角）是盆列表，一个盆对应一个吸引子，盆的编号和吸引子编号是一致的，选择哪个盆就会显示出哪个盆。这里我们在盆列表中选1来让1号盆显示，并且选show basin interior让盆内部也显示，此时看到的图如下  位于盆边界的格点通过绿球显示，使得盆的轮廓被勾勒出来。如果不想显示盆，在盆列表里选None即可。 在盆列表里会看到"Unas"（意指Unassigned）和"Boun"（意指travelled to box boundary）。前者指的是产生盆过程中由于某些数值算法原因没有成功被指认所属吸引子的格点，后者是指产生盆的过程中按照梯度最大方向移动但是最终移动到格点数据范围边界的点。这两类格点一般没有什么物理意义，对分析结果影响甚微，通常出现在离体系较远的地方。值得一提的是，如果为了节省计算时间，没有让格点数据范围囊括整个分子而是只囊括了局部区域，那么可以得到只位于局部区域的盆，与此同时"Boun"类型的点会大量出现，依然无视它即可。在盆列表里选择"Unas"和"Bound"就可以分别看到这两类点分布在哪。对于当前体系，在产生盆之后从屏幕上可以看到 The number of unassigned grids: 0 The number of grids travelled to box boundary: 0 也就是说"Unas"和"boun"类型的格点数都为0，这是最普通的情况。也因此对于当前例子在盆里表里选了"Unas"或"Boun"之后什么都不会显示。 现在点击RETURN按钮关闭图形界面。 这里略微离题一下。如果你想得到更漂亮的图，或者想更充分地自定义图像效果，可以选-4将吸引子导出为当前目录下的attractor.pdb文件，并且选-5将盆导出为cube文件。选择-5后，会提示你输入要导出的盆的编号范围，比如想把此例的3个盆都输出就输入1,3。然后程序问你是否也把盆的内部区域输出到cube里，如果选n，接下来用cube文件绘制盆的时候将只能看见边界区域。之后，这三个盆会被输出到当前目录下的basin0001.cub、basin0002.cub和basin0003.cub。然后就可以用任意可以显示等值面的工具来绘制盆了。在这些cube文件中，盆所属区域的格点数值为1，其它的格点数值为0，所以只需显示isovalue为0.5（或其它的0~1之间的数）的等值面就能让盆显示出来了。比如我们用VMD将basin0002.cub载入，令分子结构以CPK方式显示，同时用透明方式显示出盆，即可看到如下效果  如果想让边界显得更圆滑，之前生成格点数据时可以选"High quality grid"。如果想让吸引子显示在图中，载入attractor.pdb文件然后用诸如point或者bead方式显示它即可。下面的图是通过VMD利用high quality grid的格点数据产生的，吸引子以紫色圆球描述，1和3号盆的边界以不同方式显示了出来  回归正题。在我们生成了盆以后，就可以通过盆分析界面里的选项2~5对盆进行各种分析了。先做哪种分析都无所谓。这里先选"2 Integrate real space functions in the basins"做盆积分。我们这里要在盆里积分的是电子密度，所以接下来可以选1，来让Multiwfn计算各个盆内的格点上的电子密度值以得到盆内的密度积分值。但是，由于前面已经生成过了一遍电子密度格点数据，这套数据仍然在内存里存着，因此就没必要让Multiwfn重新算每个点的密度了，直接调用内存里的格点数据的值就行了，这大大节省了时间。所以我们选"0 The values of the grid data stored in memory"，来让内存中存着的格点数据作为被积的函数。立刻结果就输出出来了 #Basin Integral Volume(Angstrom^3) 1 0.7567576964 298.34200000 2 5.2789080691 473.81000000 3 7.9481842894 732.50400000 Sum of above values: 13.98385005 除了每个盆的积分值还有盆的体积。不过盆的体积对于当前体系没什么意义，因为三个盆的总和就是整个格点数据空间范围，故格点数据空间范围取得越大盆的体积就越大。 第一个盆的积分值是0.757，也就是说氢的AIM原子空间的电子布居数是0.757，因此氢的AIM电荷是1.0-0.757=0.243。当前体系含有14个电子，但是积分值的总和与理想值有微小差距。如果想改善结果，建议在前面生成盆的时候选high quality grid，这样结果就和理想值相当符合了，当然，作为代价，计算时间会变长不少。 接下来，对盆进行电多极矩分析。电多极矩的计算公式在手册3.18.3节里有详细说明，只不过那些公式里的原子核坐标对于盆分析来说应当被替换为吸引子坐标，而积分范围应当由模糊空间替换为盆空间。 由于电子密度在原子核附近比较大的缘故，如果生成盆时的格点精度不够高，做盆的多极矩分析时由于在核附近积分误差较大，结果将不够精确。所以我们先把盆在较高格点精度下重新生成一遍。故输入 1 //重新生成盆并重新寻找吸引子 1 //计算新的格点数据，而非使用在内存中存着的格点数据来生成盆 1 //让Multiwfn计算电子密度格点数据 3 //高质量格点精度 重新生成盆并寻找吸引子之后，由于序号有可能会变化，所以建议再进入功能0里看看序号变没变。这里假定序号没变。然后选择功能3进入盆的多极矩分析界面。之后输入2，来得到2号盆（对应碳的AIM原子空间）的结果，如下所示 Basin electric monopole moment: -5.290550 Basin electric dipole moment: X= -0.006981 Y= -0.006960 Z= 1.076573 Magnitude= 1.076618 Basin electric quadrupole moment: QXX= -0.763834 QXY= 0.000114 QXZ= 0.002346 QYX= 0.000114 QYY= -0.763834 QYZ= 0.002294 QZX= 0.002346 QZY= 0.002294 QZZ= 1.527667 Its magnitude: 1.527667 盆的单极矩就是盆内电子密度积分值的负值。此盆的偶极矩Z分量为明显的正值，表明盆内的电子密度主要分布于相对于2号吸引子（碳原子核）位置的Z值更负的区域。四极矩的ZZ分量为正，而XX、YY分量为负，表明盆内的电子密度分布相对于球对称分布来说，在Z方向（分子轴方向）有所延展，而在其余方向有所收缩。四极矩的大小表现盆内电子密度偏离球对称的程度，此例的数值1.528不小，说明偏离球对称分布还是挺大的，这显然是因为形成了化学键所引起的。可见，通过盆的电多极矩分析，我们可以以定量的方式明确地讨论盆内电子的分布状态。实际上计算诸如盆内的八极矩、十六极矩也很容易做到，但是用处不是很大，所以Multiwfn最高给到四极矩。 最后，基于AIM盆，我们分析定域化指数(LI)和离域化指数(DI)。输入0退出多极矩分析界面后，选择功能4。等一会儿之后，Multiwfn就会输出盆之间的DI矩阵和盆内的LI值。如下所示 ************ Total delocalization index matrix ************ 1 2 3 1 0.97545421 0.90309943 0.07235478 2 0.90309943 3.52649949 2.62340006 3 0.07235478 2.62340006 2.69575485 Total localization index: 1: 0.262 2: 3.527 3: 6.612 比如DI(2,3)=2.623，表明2号盆（C的AIM原子空间）和3号盆（N的AIM原子空间）之间平均来说共享了2.623个电子，这和HCN中C与N的形式键级3.0挺接近。实际上，DI就是一种共价键级的衡量标准，随着键的极性的增强DI值一般会降低。DI矩阵的对角元是相应的列（或行）的矩阵元的加和，如果把DI就当成键级看待，那么对于闭壳层体系，DI矩阵的对角元便是原子价。因此，可以说此体系N原子的原子价是2.69。 1号盆（H的AIM原子空间）的定域化指数为0.262，说明平均有0.262个电子定域在这个盆内，这个数值远小于此盆内电子布居数，表明在HCN体系内，电子很容易从H的AIM空间中离域出去，外界的电子也容易离域进来，换句话说，盆内外的电子容易相互交换。 4.2 乙炔的ELF盆分析 启动Multiwfn并输入 examples\C2H2.wfn 17 1 //生成盆并寻找吸引子 9 //ELF 2 //中等质量格点 0 //观看结果。图像如下所示  理论上，在两个碳之间有一个环形的ELF吸引子，这是由于C-C之间存在的两个pi键所致。在Multiwfn中，吸引子都是以点来描述的。由于数值算法的原因，如图所见，Multiwfn会在这个环的位置上找出来一大堆数值差不多的吸引子，一起围成一个环形。显然不能让这一圈吸引子一个个地孤立着，否则物理意义就丧失了，所以Multiwfn会自动地将挨得近的且数值相仿佛的吸引子进行“归簇”(clustering)。这一圈吸引子也因此就被归为了一个简并的吸引子(degenerated attractor)，即图中所示的2号吸引子，原本那一圈吸引子就被称为2号吸引子的“成员吸引子”。而2号盆的空间也就是其成员吸引子对应的盆的总和。有少数时候，在默认设定下Multiwfn自动归簇的结果可能与预期的有差距，比如这一圈吸引子没有被归在一起而被归成了两个吸引子，这时应当在生成盆之前先用功能-6调节阈值，包括相对数值差异和距离间隔这两个条件阈值。如果任意两个吸引子之间的相对数值差异小于阈值，而且它们的距离也小于阈值的话，则它们会被归到一起。距离阈值是和格点间距有关的，假设定义的时候输入了一个值n，那么距离阈值就是n*sqrt(dx^2+dy^2+dz^2)，其中dx、dy和dz是三个方向的格点间距。如果将n设为0，那么Multiwfn就不会在寻找完吸引子之后对它们进行归簇。另外，在生成盆并搜索完吸引子之后，也可以用功能-6的选项3来自行指定将哪些吸引子归到一起，后面的例子将会用到这个功能。 根据ELF的符号法则，2号盆应该被指认为V(C1,C3)，这里V代表这个盆的电子是来自于C1和C3的价层电子。类似地，1号和5号盆应当被分别指认为V(C3,H4)和V(C1,H2)。3号和4号盆应当被分别指认为C(C3)和C(C1)，C代表Core，表明内核电子主要都在这个些盆里。在图形界面中取消Show molecule复选框不让分子结构显示出来就能看到这种盆了，如下所示  接下来计算电子密度在ELF盆中的积分值。关闭图形窗口后选择功能2，然后选择选项1（对应电子密度），就会看到如下结果 #Basin Integral Volume(Angstrom^3) 1 2.2160904060 460.52100000 2 5.3886831176 700.01500000 3 2.0849480277 0.81300000 4 2.0846030131 0.81700000 5 2.2103993998 459.13000000 Sum of above values: 13.98472396 C(C1)和C(C3)的积分值为2.08，很接近2，这也正对应了碳有两个内核电子。V(C3,H4)和V(C1,H2)的积分值为2.21，也接近2，对应于C-H键共享一对儿电子这个经典化学概念。值得一提的是，尽管此C(C1)和C(C3)是对称的、V(C3,H4)和V(C1,H2)是对称的，因此原则上它们的积分值应该相同，但实际上因为格点精度是有限，所以积分结果有一些数值差异，使用更高的格点精度能使结果与分子对称性满足得更好。按照经典化学键理论，乙炔的两个碳之间应该享有三对儿电子，共6个电子，但是V(C1,C3)的积分值5.39与6偏差略大。不过也没有理由要求这种基于实空间的分析的结果与经典化学概念完全一致。 然后选功能3再选-1，来得到每个盆的电多极矩分析结果。其中的4号盆，即C(C1)的结果为 Basin 4 Basin electric monopole moment: -2.084603 Basin electric dipole moment: X= 0.080293 Y= 0.080293 Z= -0.099329 Magnitude= 0.150865 Basin electric quadrupole moment: QXX= -0.000597 QXY= -0.004658 QXZ= 0.005628 QYX= -0.004658 QYY= -0.000597 QYZ= 0.005628 QZX= 0.005628 QZY= 0.005628 QZZ= 0.001194 Its magnitude: 0.001194 由于当前体系是轴对称的，而且轴在Z轴上，所以盆内的偶极矩的X和Y分量都应当为0。不过由于这次只用了medium quality grid，所以这两个分量不是很接近于0。但即便如此，当前的结果还是有定性意义的，起码看到偶极矩的大小非常小，而且四极矩的大小也很接近0，这说明原子的内核电子受分子环境影响甚小，基本接近于球对称分布。 输入0返回上一级菜单，然后选功能4计算LI和DI，结果如下 ************ Total delocalization index matrix ************ 1 2 3 4 5 1 1.31120729 1.03620018 0.16053897 0.02258264 0.09188551 2 1.03620018 2.70782782 0.31823886 0.32035614 1.03303264 3 0.16053897 0.31823886 0.50788903 0.00679996 0.02231124 4 0.02258264 0.32035614 0.00679996 0.50834742 0.15860868 5 0.09188551 1.03303264 0.02231124 0.15860868 1.30583807 Total localization index: 1: 1.560 2: 4.034 3: 1.825 4: 1.823 5: 1.557 DI(3,4)对应于C(C1)和C(C3)之间的离域化指数，其数值0.007十分接近于0，这表明在两个原子的内核区域之间电子很难相互离域。DI(1,3)=0.160和DI(2,3)=0.318虽然不大，但是还是有一定数值的，这说明C(C3)盆内的电子与和它相邻的仅有的两个盆V(C3,H4)和V(C1,C3)内的电子相互交换还是有一定几率的。DI(2,1)和DI(2,5)都接近于1，表明电子在C-C键区域和C-H键区域间离域很容易。从定域化指数上看，碳的内核的盆的LI为1.82，接近于盆内的电子布居数2.08，表明内核电子整体来说喜欢呆在内核区域，而不容易离域到外头去，而相应地外界电子也难以进来。而V(C3,H4)的LI为4.034，和它的电子布居数5.389差得挺大，表明V(C3,H4)内的电子容易与外界电子发生交换，这也是价层ELF盆的普遍情况。 4.3 水分子的静电势盆分析 虽然很早以前已经有人对静电势做过拓扑分析，但是对静电势进行盆分析的研究还相当少。 静电势同时存在原子核电荷主导的正值的区域，以及电子所主导的负值的区域。对于这种正负值同时存在的实空间函数，Multiwfn的盆分析模块在生成盆并搜索吸引子之前，会首先将负值区域的数值乘以-1来让符号颠倒，然后照常生成盆和寻找吸引子，弄完了之后再恢复原先的负值区域的数值。这样一来，负值区域也会找到吸引子并能产生对应的盆，但这些吸引子实际上是“排斥子”，即函数值的极小点。在Multiwfn程序中以及下文中，不管本质上是吸引子还是排斥子，一律都被叫做吸引子。 启动Multiwfn并输入 examples\H2O.fch //当然，也可以用.wfn或.wfx文件作为输入 17 1 //生成盆并搜索吸引子 12 //计算静电势格点数据 2 //中等质量格点 0 //查看结果 从显示出来的图中可以看到1、4、5号吸引子对应于原子核电荷产生的静电势极大值点。而2、3号吸引子对应于氧的孤对电子导致的静电势负值区域中的极小点。诸如这样的负值区域的极小点以及相应的盆都会用蓝色表示，以便与绿色表示的正值区域的极大点和盆进行区分。从盆列表里选择3号盆，可以看到如下图像，  注意如命令行窗口中的提示所显示的，在此例中生成盆的时候有800多个格点最后跑到了格点数据范围的边缘，在盆列表里选择"Boun"就能看到它们，如下图所示。显然这些点里分子都很远，所以这个问题可以无视。  盆分析模块提供了便利的计算吸引子、原子核间的几何参数的功能。进入功能-2后可以看到操作方法介绍。例如想测量2号吸引子与1号原子（氧原子）间的距离，就输入a1 c2，结果为2.260 Bohr。再比如，测量“吸引子2--氧原子核--吸引子3”的角度，就输入c2 a1 c3，结果为76.96度，这在某种程度上也可以看做是氧的孤对电子的夹角（孤对电子夹角还可以利用ELF、LOL来研究，结果与此有一定差异，见《使用Multiwfn做拓扑分析以及计算孤对电子角度》http://hi.baidu.com/sobereva/item/c84fa5089156806dd55a11c3）。然后输入q退到上一级菜单。 假设我们想手动把两个孤对电子对应的吸引子合并为一个吸引子，可以输入 -6 //设定归簇参数或自行对吸引子进行归簇 3 //将指定的吸引子归到一起 2,3 //2和3号吸引子将归为一个，同时盆也合并到一起 0 //返回 这时再选0打开图形界面，会看到吸引子的序号都变了，两个孤对电子对应的吸引子合并成了一个，序号为2。 选功能2，然后选1，来得到电子密度在静电势盆中的积分值，结果为 #Basin Integral Volume(Angstrom^3) 1 0.8682191150 239.18400000 2 0.7388522542 521.25800000 3 7.5266332598 26.99700000 4 0.8623600202 240.36100000 Sum of above values: 9.99606465 Integral of the grids travelled to box boundary: 0.00000011 虽然氧有两个孤对电子，但是在两个孤对电子对应的静电势盆，即目前的2号盆当中电子密度积分值仅有0.739。这是因为对任何分子，大部分离分子较近的空间总是由核电荷所主导，故而静电势为正。静电势为负值的区域都处在分子的较外侧，这些区域电子密度是比较有限的，所以电子密度积分值在负的静电势盆中不会很大。结果也显示，那些"Boun"类型的格点内的积分值几乎为0，这也是为什么前面说这些格点可以忽略。 在Multiwfn的盆分析模块中，允许对外部格点文件(.cub/.grd)记录的实空间函数产生盆并寻找吸引子。也可以将外部格点文件记录的函数作为被积函数得到它在各个盆内的积分值。这个设计使得盆分析模块的普适性、灵活性很强。比如要分析的实空间函数用Multiwfn没法直接计算，就可以考虑用第三方程序计算它的格点文件再用Multiwfn对它做盆分析。而对于静电势，用其它量子化学程序（如Gaussian自带的cubegen）计算会比用目前的Multiwfn要快（尤其是大体系），因此可以用这些程序代替Multiwfn计算静电势数据来节省静电势盆分析的总时间。在手册的4.17.4节的Part2中有个实例，介绍了利用cubegen生成的静电势和电子密度格点文件来做与本节相同的分析。 4.4 水分子的电子密度差的盆分析 在进入盆分析之前，首先要利用主功能5把电子密度差的格点数据算出来。电子密度差这里具体指分子密度与组成它的各个原子在孤立状态下的密度差。这需要涉及到生成原子波函数文件，为了省事，此例就直接用Multiwfn自带的原子波函数文件。做法也就是把examples目录下的atomwfn目录挪到上一级目录下（假设你是通过双击Multiwfn图标来运行Multiwfn的，否则就将atomwfn放到你调用Multiwfn时所在的目录下）。关于用Multiwfn做电子密度差的具体内容，见《使用Multiwfn作电子密度差图》（http://hi.baidu.com/sobereva/item/96eaeb5dd1b12309e7c4a5c1） 移动过atomwfn目录后，启动Multiwfn，然后输入 examples\H2O.fch 5 //计算格点数据 -2 //产生变形属性 1 //电子密度 3 //高质量格点。注意这里的“高质量格点”和前面用到的盆分析模块里的高质量格点不是一回事。这里定义的是总格点数而非格点间距，故而高质量只是对于较小分子来说的。由于密度差比电子密度、ELF变化都要复杂，为了盆分析可靠，建议用较好的格点质量。 0 //返回主菜单。刚才计算的格点数据会一直保存在内存里。 17 //盆分析 1 //生成盆并寻找吸引子 2 //基于保存在内存里的格点数据来产生盆和寻找吸引子 选择0进入图形界面观看电子密度差的吸引子和相应的盆。显示和不显示分子结构时的图分别如下面两图所示   电子密度差的正值和负值部分分别代表原子形成分子后，电子密度增加和减少的区域。上图中绿球对应正值部分的极大值，蓝球对应负值部分的极小值。如果搞不清楚为什么极值点这么分布，对照着电子密度差截面图就会弄明白。水分子垂直于分子平面和分子平面上的电子密度差图如下所示  对照平面图，就很容易理解3号和6号盆对应的是形成O-H键造成的电子密度增加区域，而4号和5号盆对应的是氧原子以sp3杂化后由于孤对电子的出现造成的密度增加的区域。 值得一提的是由于格点精度用得仍然不够高，导致7号吸引子的另一侧的等价吸引子没有被找出来。另外，8号吸引子偏离了理想位置。不过，这样的吸引子和盆都不是化学上感兴趣的，它们对应的是密度差变化的十分细微的结构，没有太大物理意义，所以可以忽视这些问题、 现在看看O-H键形成后造成了多少电子往键中间发生了挪动，也就是对3号或6号盆进行积分，被积函数也是电子密度差。于是，进入功能2，然后输入0，来将内存中保存着的电子密度差格点数据当做被积函数。结果立刻显示出来，3号和6号盆的积分值都是0.09。而对应于孤对电子的形成造成密度增加的4号或5号盆的积分值是0.21。 以前zhou2009在小木虫上发了很多关于分析电子密度差的帖子，据悉是通过excel手动对指定X,Y,Z空间范围内特定符号的电子密度差数据进行加和然后乘上空间微元得到积分值。这样的操作很耗费精力，而且只能处理很简单的体系。现在利用Multiwfn 3.0的盆分析功能，做这样的分析真是极其简单，而且更精确。 如果想要将等值面图和吸引子进行对比，可以选择-10退回Multiwfn主菜单，然后进入主功能13，再选-2，当前内存中的格点数据的等值面图就会显示出来，之前找出来的吸引子也会同时显示。下图是isovalue=0.05的图  PS:关于保存在内存里的格点数据这里多说几句。用.cub/.grd文件作为输入文件的话，其中的格点数据就会一直保存在内存里。另外用主功能5或者13计算或者对格点数据进行操作的话，得到的格点数据也会保存在内存里。当有新的格点数据生成时，旧的格点数据就会被覆盖。若在盆分析功能里的功能1当中让Multiwfn计算新的格点数据，那么原先内存里的格点数据就也被覆盖。 4.5 对乙烷进行基于AIM盆的源函数分析 源函数在Multiwfn手册2.6节有简要介绍。在The Quantum Theory of Atoms in Molecules-From Solid State to DNA and Drug Design一书的7.6章有个简明的综述建议读读。这个函数依赖于两个坐标，其中一个坐标是参考点。当源函数被用来讨论成键（包括弱相互作用）问题时，参考点一般取在BCP（键临界点）上。经常通过在AIM原子空间内对源函数进行积分，以研究原子对BCP处密度的贡献。我们这一节要分析乙烷的各部分对C-H键的BCP密度的贡献值。我们先要进行拓扑分析找出临界点位置。启动Multiwfn，依次输入 examples\ethane.wfn 2 // 拓扑分析 2 // 搜索核临界点 3 // 搜索BCP 0 // 观看结果 临界点分布和编号如下所示  这里我们打算将CP13这个C-H键的BCP取做参考点。从文本窗口能找到它的坐标，将坐标直接从屏幕上拷贝到系统的剪切板中（如果不知道怎么去做这步，见手册5.4节）。然后我们看看这个点的电子密度是多少，这一步不是必须的，但是由于对源函数进行全空间积分的结果原理上会精确等于参考点的电子密度值，所以我们现在获知了它的电子密度，就便于稍后检验积分得是否准确。选择7，然后输入临界点编号13，就能得到它的密度0.276277。 现在来设定参考点。一种方法是用settings.ini里的refx、refy和refz参数来设参考点的X,Y,Z坐标，但是这样设完了还得重启一下Multiwfn才能生效，这里我们用另一种方法来设。先选-10从拓扑分析界面返回到主菜单，然后输入1000，就进入了一个秘密界面，选择1，再把CP13的坐标粘贴进去并按回车，参考点坐标就设好了。 然后，依次输入以下命令生成AIM盆 17 // 盆分析 1 // 生成盆并寻找吸引子 1 // 电子密度 2 // 中等质量格点 进入功能0，会看到吸引子如下图分布  因为随后我们想看下方的甲基对BCP密度的贡献，为了免得稍后得手动去加和积分值，我们这里将1、2、3、4号盆都合并到一起成为甲基的AIM盆。故关闭图形界面后输入 -6 3 1,2,3,4 // 要归簇的吸引子编号 0 这时再进入图形界面，就会看到1、2、3、4号盆已合并成了1号盆。 接下来要对盆进行积分。这里特别要注意的是这次不能用功能2来获得积分值，因为源函数和拉普拉斯函数类似，在原子核附近变化特别剧烈，普通方式的立方格点积分精度太低，连定性正确的结果也得不到。所以Multiwfn特意为了在AIM空间下积分源函数、拉普拉斯函数，设计了功能5。它与功能2的区别是积分时候使用了多级优化，大意是指在电子密度越高的地方，每个格点就会分裂为越多个格点以使得格点间距变小，从而得到越高的积分精度。这使得原子核附近积分精度有了显著的提高，但是作为代价，计算量也变大了不少。所以一般的实空间函数积分用功能2就够了（控制优化方式的具体参数可以进入功能5后用选项-1先进行调节，这里暂不提及） 选择功能5，然后选择19（源函数），结果如下 #Basin Integral Volume(Angstrom^3) 1 0.2683862633 883.82800000 2 0.0027953503 235.13700000 3 0.0027789566 250.98000000 4 0.0061607390 124.72900000 5 0.0028977125 269.88600000 Sum of above values: 0.28301902 积分值总和0.28301902与参考点的密度0.276277十分接近，相对误差仅有2.44%，完全在可接受范围内。1号盆，也就是目前甲基的盆，对BCP密度的贡献为0.2684/0.2830*100%=94.8%，因此甲基是它自身的C-H键BCP密度的最主要贡献者，而其余原子只贡献5.2%。有人发现这个结论对于不同长度的烷烃都是十分一致的，这被用来说明化学基团的可移植性问题。 注意尽管利用功能5做源函数积分对于当前体系可以得到满意的精度，但是如果涉及到比较重的原子，比如Br，由于在内核部分函数波动过于剧烈，基于立方格点的积分方式无论怎么改进也得不到理想的结果。这不得不用解析的方法来积分，而解析方式的AIM盆积分或许会在未来Multiwfn版本中实现。 5 总结&其它 本文介绍了盆分析的意义、算法，结合很多例子对Multiwfn的盆分析操作进行了讲解。对于最一般的盆分析，在Multiwfn里其实就这么几步，简单得很： [文件名] // 一般用.wfn/.wfx/.fch文件 17 // 盆分析模块 1 // 生成盆并寻找吸引子 1 // 选择用什么实空间函数定义盆和吸引子，这里假设为电子密度 2 // 产生什么质量的格点数据，这里是中等质量的，适合用作预览 0 // 观看结果 然后根据要分析的内容，选择功能2~5即可。 新加入盆分析模块的Multiwfn 3.0的推出，对于实空间函数的盆分析领域，是具有革新意义的。Multiwfn在普适性、速度、灵活度、简便易用等各种方面相对以往程序都有着质的飞跃。特别是学术界既十分常用又超级难用的ELF盆分析程序TopMod，在Multiwfn 3.0推出后就没有任何存在价值了，Multiwfn将ELF盆分析的门槛降低到量化初学者也能轻松完成的程度。本文虽然用的例子都是小分子，但用于较大的分子也没有问题。除非是分子太大，这样的话可能会出现内存容量不足以储存格点数据的问题。对于这种情况，如果只对局部区域感兴趣，建议只将格点数据空间范围定义在相应的局部区域，来减少总格点数从而避免内存不足。 [ Last edited by zhou2009 on 2013-4-12 at 16:06 ] |
回复此楼» 收录本帖的淘贴专辑推荐
 纳米技术与能源及模拟 | 万卷阁之搜神记 | VASP | 新药研究 |
| 高斯/MOPAC | 量化软件学习 | SOB 集锦 | 喜欢的 |
| 量化计算 | 计算化学-辅助工具 | 第一性原理计算 | 量子化学-第一性原理 |
| 量化经验 | 量化软件-Multiwfn | nature and science | 专业科研软件 |
| 高斯专辑 | 科研感悟 | 路哥的淘贴 |
» 猜你喜欢
电催化或能源材料
已经有2人回复
-
关于电解水稳定性测试的电解液问题
已经有0人回复
-
物理化学论文润色/翻译怎么收费?
已经有284人回复
-
电催化或能源材料合作
已经有0人回复
-
[博士招聘]伍伦贡大学 University of Wollongong 王彩云老师课题组博士招聘
已经有0人回复
-
山东科技大学（青岛）化工学院招收还有调剂名额（仅限化学专业），速联系
已经有15人回复
-
硅碳负极材料问题
已经有0人回复
» 本主题相关商家推荐: (我也要在这里推广)
» 本主题相关价值贴推荐，对您同样有帮助:
- 衡量芳香性的方法以及在Multiwfn中的计算
已经有12人回复
- 使用Multiwfn绘制原子轨道图形、研究原子壳层结构及相对论效应的影响
已经有17人回复
- 使用Multiwfn作电子密度差图
已经有37人回复
- 使用Multiwfn做拓扑分析以及计算孤对电子角度
已经有18人回复
- 回答一些关于Multiwfn的疑问以及未来Multiwfn的发展打算
已经有37人回复
- 【sobereva个人文集】电子定域性的图形分析
已经有56人回复
24楼2013-04-10 19:29:05
3楼2013-03-27 07:59:54
5楼2013-03-27 08:37:34
7楼2013-03-27 08:45:20
90