使用Multiwfn做拓扑分析以及计算孤对电子角度
文/Sobereva 2011-Dec-1

Multiwfn程序(http://multiwfn.codeplex.com)从2.1版开始支持电子密度的拓扑分析，在刚刚发布的2.2版中进一步支持ELF、LOL、电子密度拉普拉斯函数、轨道波函数的拓扑分析。本文将简要介绍Multiwfn拓扑功能的基本使用方法。

电子密度函数大家都很熟悉，如果不熟悉ELF、LOL和拉普拉斯函数，可以参看《电子定域性的图形分析》一文，见http://hi.baidu.com/sobereva/blo ... e99b0c4a90a7d6.html。


1 拓扑分析简介

限于篇幅，这里只是对拓扑分析做极简单的说明，或者说是回顾。对于完全不了解AIM理论的人，不妨看看The Quantum Theory of Atoms in Molecules-From Solid State to DNA and Drug Design一书的第一章，对电子密度的拓扑分析有简介、清楚、易懂的介绍。

实空间函数（以三维空间坐标为变量的函数）的拓扑分析主要是指获取临界点，以及获得连接临界点的拓扑路径。电子密度的拓扑分析是由Bader首先提出的，是AIM理论中的重要组成部分，后来拓扑分析的思路又被Silvi、Savin等人用于分析ELF（电子定域化函数）。

临界点(Critical point, CP)是函数的梯度的模为0的点，分四类，以(3,X)表示。假设实空间函数的Hessian矩阵（3*3的二阶导数矩阵）的本征值有M个正值和N个负值，则X=M-N。
(3,-3)对应函数的局部极大点。对于电子密度函数，通常出现在离原子核很近的位置。对于ELF函数，通常出现在共价键区域、孤对电子区域、原子核处。
(3,-1)对应函数的二阶鞍点。函数在一个方向曲率为正，另两个方向为负。对于电子密度函数，通常出现在有相互作用的两个原子之间，也被称为键临界点(BCP)。对于ELF函数，位置不甚确定，经常偏离两个相关的(3,-3)的连线很多。
(3,+1)对应函数一阶鞍点，如同势能面上的过渡态。对于电子密度函数，通常出现在环体系平面中，如苯环的中心。对于ELF，(3,+1)和下面的(3,+3)通常较少讨论。
(3,+3)对应函数的局部极大点。对于电子密度函数，通常出现在笼状体系中，例如C60的中心。

寻找临界点通常使用牛顿法。给出一个初猜点，通过反复迭代，就能找到与之最邻近的临界点。但牛顿法不一定每次都成功，有可能一直不收敛，也可能迭代过程中碰到Hessian矩阵为奇矩阵的情况，或者收敛到已找到的临界点位置，此时就要尝试从下一个初猜点开始了。初猜点与期望的临界点位置越近越容易成功收敛到相应的临界点。

电子密度函数的结构很简单，一般都是在原子核处像一个山峰耸立，向着四周呈指数降低，典型体系的各种临界点位置都能大致猜出来，正如上面已经提到的。寻找(3,-3)点的时候初猜点都选在原子核位置。寻找(3,-1)点的时候，初猜点都选在两个原子或者两个已找到的(3,-3)点连线的中点，在这个过程中也可能顺便找到(3,+1)和(3,+3)点。另外也常常将三个原子或四个原子的几何中心作为初猜，这分别适合进一步寻找尚未找到的(3,+1)和(3,-3)点。假设利用这几种初猜方式还是没找到一些期望的临界点，但通过经验可以估计出期望的点出现的大致空间范围，那么就应当在这个空间范围内设定一大堆初猜点来尝试寻找那个临界点。

对于ELF、LOL和电子密度拉普拉斯函数，寻找对应原子内核处的临界点（对应高度定域的1s轨道电子）还是用原子核位置作为初猜。但是对于价层的ELF临界点，由于其分布较难预料，而且数目众多，最好的办法就是直接在每个原子附近都设一大堆初猜点来寻找各个临界点。不过，由于ELF临界点太多而且分布太复杂，想找全所有临界点绝非易事，所以，只要看到化学意义显著的临界点（所有的(3,-3)，以及部分(3,-1)）都已经找到了，临界点的搜索就可以停掉了。除了牛顿法，还有一种基于立方格点的临界点搜索方法，在TopMod中用的就是这种，但数值精度上令笔者比较鄙夷，我不打算在这里细谈（立方格点方法和牛顿法相结合倒是可能产生出高效且准确的方法）。

拓扑路径是连接临界点的路径。一般比较关注的是(3,-1)与(3,-3)之间的路径，沿着(3,-1)的Hessian矩阵本征值为正的本征向量方向分别向前和向后出发，顺着梯度方向不断前进就能找到两个(3,-3)，而经过的轨迹就是拓扑路径。对于电子密度函数，这样的拓扑路径近似于化学键，因此称为键径。对于ELF，通过生成拓扑路径则可以将(3,-1)与(3,-3)临界点之间的关系清晰地表现出来。


2 电子密度拓扑分析在Multiwfn中的操作

这里先介绍电子密度拓扑分析的操作，以Gaussian在HF/6-31G*下优化并生成波函数文件的多巴胺分子为例。先启动程序，然后依次输入

dopamine.wfn  //波函数文件
2  //拓扑分析功能
2  //以所有原子核作为初猜点，目的是寻找所有核临界点
3  //以每一对儿距离不很远的原子核连线中点作为初猜点，主要用来寻找键临界点
8  //生成键径
0  //观看结果
此时文本界面会输出所有临界点和键径的信息，最后提示你Fine, Poincare-Hopf relationship is satisfied，这是所有临界点都已找到的必要非充分条件，如果这个条件没满足，肯定临界点没找全。同时蹦出来一个图形界面，如下所示

紫色的球是(3,-3)，与原子核位置基本一致。桔黄色是键临界点，黄色是环临界点。棕色的线是键径。通过图形界面，可以选择显示哪些、不显示哪些，以及旋转缩放等操作。凭借化学直觉，以及Poincare-Hopf关系可知期望的临界点都已经找到了，就没必要再用选项4（三个原子中心作为初猜）、选项5（四个原子中心作为初猜）或选项6（在某个位置附近随机设一大堆点）来进一步找临界点了。

简而言之，对于一般的体系，找电子密度的临界点就是进入拓扑分析界面后，依次选2、3、4就行了，简单至极！

值得说明的是Multiwfn的代码经过充分优化，并且支持SMP并行，这使得其拓扑分析速度极快。在一般的主流四核机子上，上面的任务在Multiwfn里眨眼间就完成了。我曾比较过AIM2000和Multiwfn分析一个含900个高斯函数的C60体系，在四核并行的机子上，Multiwfn寻找临界点的速度是AIM2000的约100倍，生成键径速度是其将近30倍，而结果完全一致。AIMALL的速度比AIM2000要快，不过由于算法不很一致，也不好公平地比较，但从实际应用角度上来说，Multiwfn比AIMALL要快一个数量级。另外AIMALL的操作远不如Multiwfn这种交互式程序便利，其结果可视化程序aimstudio还要收费。可以说，想要操作简单、能直接可视化、免费而且速度飞快的AIM拓扑分析程序，Multiwfn是不二之选。

有些临界点比较特殊，在默认的设定下，用各种搜索方法可能都找不到临界点，这就需要人为调整一下参数了，在Multiwfn里提供了丰富的可调参数。比较典型的就是重原子的原子核的临界点找不到，这是因为原子核处电子密度极大值对应的峰太尖锐，难以定位。解决办法是适当降低收敛标准，也就是进入选项-1，将梯度和位移的收敛标准降低，然后重新用选项2以核为初猜来搜索临界点。另一种情况就是重原子和氢连线上的键临界点找不到，这就需要在参数设定界面中将Set scale factor for stepsize一项设低，建议设为0.5，然后再重新找临界点。其它选项在Multiwfn手册第三章的拓扑分析一节里都有说明。

Multiwfn还提供了一些选项用于分析、编辑结果。用选项7并且输入临界点的编号，就可以获知这个临界点的一切信息，包括Multiwfn支持的各种实空间函数的数值和电子密度的梯度向量和Hessian矩阵。利用-9可以方便地计算临界点间、原子间的距离、角度。进入-4则可以导入导出已找到的临界点、删除或人为地增加临界点。进入-5则可以显示拓扑路径的详细信息（路径轨迹的坐标、路径长度等）并进行编辑和导入导出。在Multiwfn里，拓扑路径实质上是以一连串的点来表示的。

除了寻找临界点和生成拓扑路径，Multiwfn的拓扑分析模块还能生成并绘制盆分界面，对于电子密度分析，就对应于原子间分界面(Interatom surface, IAS)。在IAS面上没有电子密度梯度线穿过，也称零通量面(zero-flux surface)，每个IAS都经过一个键临界点。IAS将整个分子区域分隔成一个个原子空间，也叫原子盆。从拓扑分析主界面进入选项10后，输入(3,-1)的编号，就能生成穿过这个(3,-1)的IAS。而在编号前加上负号，就说明将已经生成的穿过这个(3,-1)的IAS删掉。如果当前尚未生成任何IAS，输入0就可以生成所有IAS。如果已经有了IAS，输入0就会删掉所有IAS。下图中生成了部分IAS



3 ELF拓扑分析在Multiwfn中的操作

这一节谈谈ELF的拓扑分析操作，以苯酚为例。

phenol.wfn // HF/6-31G*下优化后得到的波函数
2 //进入拓扑分析模块
-11 //清空已有拓扑分析数据并选择新的实空间函数。默认的函数是电子密度
9 //选择ELF
6 //在指定半径和中心的球里随机分布指定数目的初猜点来搜索临界点
-1 //将每个原子核依次作为球中心。由于此分子有13个原子，如果提示中"Points:"此时显示的是比如1000，总共就会以13*1000=13000个初猜点来找临界点。在我的i7-2630QM机子上以四核并行模式运行，仅15秒钟就完成了搜索。
-9 //返回上一级菜单
0 //观看结果
你会看到结果远比电子密度分析时复杂，总共找出了126个临界点，很花哨。为了使结果更清楚，可以将意义不大的(3,+1)和(3,+3)临界点通过右侧复选框关掉。另外，原子球挡住了一些临界点，将Ratio of atomic size滑动条拖到最小，原子球就不显示了，分子此时相当于以棍棒模型显示。如下图所示。



依然是每一个紫球对应一个(3,-3)。可以看到每个原子核都有一个(3,-3)（对于氢原子会偏离些）。如果你发现原子核处没有紫球，说明它们搜索中漏过去了，可以用选项2以原子核位置为初猜再找一下，一般都能找到。C-C和C-O键之间都有紫球，这对应于共价键区域电子定域程度局部最高点。在氧原子旁边还有两个悬着的紫球，它们对应于氧的两个孤对电子的电子定域性局部最大点。从位置上看，按照传统杂化理论可以认为氧是处于sp3杂化状态。

在显示结果的时候，命令行窗口会提示Warning: Poincare-Hopf relationship is not satisfied, some CPs may missing，这表明这一次搜索并没有将所有临界点找全。不过没有关系，因为从图中可见所有有显著化学意义的临界点都找到了。如果你想彻底找全，那么就再次进入选项6，适当选择参数后重新进行临界点搜索，将会找出零星数个剩余的临界点。彻底找全所有类型ELF的临界点，对于不是很小的体系是很困难的事，一般没有必要去做。对于最重要的临界点类型，即(3,-3)，基本上做一次每个原子周围有1000个初猜点的搜索（即本例的操作）后就都能找到。

现在点RETURN关闭图形界面回到拓扑分析主菜单，选8来生成(3,-1)和(3,-3)间的拓扑路径，再次显示结果，如下图所示



原本看起来数目众多、很乱的桔黄色小球，即(3,-1)点，通过拓扑路径与(3,-3)相连后，此时关系变得比较清晰了。

对于个别体系，有时候会发现本该出现的(3,-1)与(3,-3)间的路径没有出现，这很可能是因为生成拓扑路径的默认参数不太适合。比较有效的解决方法是进入选项-2，将Stepsize设得小一些，比如将默认值减小一两倍。注意路径的最大长度是"步长*路径的最大点数"，如果生成路径过程中，迭代次数超过了最大点数还没遇见另外一个临界点，就说明这次路径生成失败。因此，步长设小了，那么路径最大长度就会减小，如果体系中有一些很长的路径，这些路径可能在此时就没法生成了。因此减小步长的同时也应当同时加大路径最大点数。步长越小越容易成功生成复杂的路径（比如曲率特别大的路径），但是由于需要加大路径最大点数，这就会造成计算耗时越多。默认的路径生成参数是折中考虑精度和效率的结果。

LOL和电子密度拉普拉斯函数与ELF有很大相似之处，此节的做法也同样适合这两种函数，只要在选项-11里选择相应的函数然后以本例同样的方法操作即可。另外还可以用这种方法来获得轨道波函数的临界点。


4 谈谈孤对电子角度怎么算

好几次看到有人问孤对电子夹角、二面角怎么算，由于Multiwfn的拓扑分析功能正好能解决这个问题，所以顺便谈谈。如果对此问题没兴趣，可以跳过这节。

在讨论怎么计算角度这个具体问题前，先要搞清楚孤对电子应该如何严格地描述。首先要明确的是，孤对电子本身不是可观测的，只是个概念，也没法唯一地去定义，只能说哪种方式的描述更有物理意义。一般意义上的孤对电子指的就是原子价层某局部空间内运动的两个自旋互相配对的电子。一种将这个概念具现化的方法是轨道定域化方法，它通过对占据的分子轨道（往往有很强离域特征）进行酉变换，获得能够与化学键、孤对电子相对应的轨道，很大程度上重现了lewis模型。这类方法最大缺点是物理意义不严格，酉变换过程是随意的，以什么标准去做酉变换因算法而异，不同的常用的定域化方法虽然能得到相似的图景，但是在多重键、孤对电子的描述上却可能很不相同，所以我并不建议用这类方式讨论孤对电子。NBO也是属于定域化轨道一类，虽然不是通过酉变换得到，但结果上看也通常差不多。从算法上NBO分析对孤对电子的描述在寡人看来也有点不合理性，在NBO分析中是先从密度矩阵的原子对角块搜索完CR、LP型NBO才开始寻找双中心NBO，等于说生成LP的优先级要高于双中心NBO，然而实际上孤对电子和化学键是同时形成的，应当同步考虑。这也导致了NBO的孤对电子轨道看起来和杂化理论显示的图景不符。

下面是HF/6-31G*下的苯酚的16号和22号NLMO轨道图形，以及两个轨道波函数相加和相减后的图形。顺带一提，NLMO是NBO分析方法框架内的一部分，它是由NBO进一步变换得到的，通常和NBO很相似，只是离域性比NBO稍大。NLMO也可以视为是直接通过对占据的分子轨道酉变换得到的。



16和22号NLMO对应氧的孤对电子，乍看上去好像氧是sp2杂化，一个孤对电子在平面上，一个垂直于平面。但实际上我们可以将这两个轨道人为地混合形成两个新的轨道，即16+22和16-22（为了省事未考虑归一化，这对定性分析没影响）。这两个新的轨道也可以说是对应孤对电子的，这样一来氧就成了sp3杂化了，与ELF分析的结果一致。或许会有人问，将16和22号轨道那么混合有道理么？这样混合是没问题的，因为NLMO间都是正交的，因此16+22以及16-22轨道也与其它轨道正交，这等于是对占据轨道再酉变换一下罢了，完全是在轨道定域化的框架内允许的操作。所以说定域化轨道模型有严重的任意性，连孤对电子轨道都不能唯一确定，更甭提再进一步计算孤对电子角度了。

物理意义比较严格的是通过实空间函数分析，具体来说是ELF或LOL函数，结果通常较相似。电子密度拉普拉斯函数也能分析，但它与电子定域化的关系只是间接的，函数行为上也不如ELF和LOL舒服。ELF较大的区域，表明相应区域内的电子定域性较强，可以推断电子的费米穴主要分布在这个区域内，即自旋相同电子间距离会较远，而自旋相反的电子相对来说就容易呆在一起，即说电子在这个区域内配对，这正是孤对电子概念严格的物理解释。如果对电子定域性问题想有进一步了解，可参阅笔者的《电子的定域性与相关穴》一文（http://hi.baidu.com/sobereva/blo ... 7954dca7ef3f8e.html）。因此，位置适当的ELF较大的区域就可以说是孤对电子的区域，那么通过拓扑分析给出的ELF较大区域当中的极大点，就可以说是孤对电子出现位置的最大程度的抽象，利用这个点就可以很容易地计算出孤对电子角度了。

比如，按照第三节的方法找完临界点后，我们想计算CP93-O12-CP83的夹角（标号见下图。在图形界面选择CP labels就能显示各个临界点的编号），就在拓扑分析主界面里面选选项-9以进入角度测量的界面，然后输入a c93 a12 c83，结果显示为88.228度。这里a代表angle，c代表critical point，a代表atom；想再测量CP83-O12-C4-C5的二面角就输入d c83 a12 a4 a5，结果是130.496度。按q可以返回主菜单。值得一提的是每次搜索临界点给出的临界点的编号可能不同，所以对一个体系重复分析时，每次都应在图形界面里看清楚各个临界点的当前编号。



我们也可以看看用LOL和拉普拉斯函数研究这个问题能得到什么结果。在拓扑分析的主菜单选择-11然后选择LOL函数，此时所有刚才ELF的拓扑分析会自动清空，然后重新搜索临界点，我们会看到各个(3,-3)的位置和ELF分析时都差不多。再次计算氧的两对儿孤对电子间的夹角，得到的结果是119.895度，明显大于ELF给出的数值。改用电子密度拉普拉斯函数再算一次，结果是131.244度，更大了。按照一般观念，孤对电子间相互互斥比较强，也应当大于sp3杂化的109度28分。结果在一定程度上显示ELF可能会低估孤对电子的夹角。实际上有时低估得更厉害，乃至于LOL的结论是有两对儿孤对电子出现的时候，ELF却认为两对儿孤对电子合在了一起。我的建议是使用LOL来研究孤对电子，LOL比ELF形式更简单，计算速度稍快，而且一些文章指出LOL也能给出比ELF更清楚的化学图景，而LOL在物理意义上也和ELF差不多。（LOL也并非各方面都比ELF好，等有机会再细谈）


5 临界点与等值面图

我们接下来将通过等值面图讨论ELF的域(Domain)、分歧点(Bifurcation point)与(3,-3)及(3,-1)临界点的关系，这有助于理解以及验证拓扑分析的结果。

先按照第三节的过程搜索完ELF临界点，然后直接退回到主菜单，并输入以下命令
5 //生成格点文件
9 //ELF函数
2 //中等精度格点，对于小体系定性分析这就够了。这时开始计算
-1 //显示等值面

由于已经搜索过了临界点，所以临界点也在等值面图上一起显示了出来。拉动等值面数值的滑动条，可以看到随着数值由小到大，等值面如下图变化（蓝色数字是等值面对应的数值）



由图可见，在数值很小时，等值面连成一个大域（或者叫做盆，是指等值面包含的封闭空间），也表明能超过这个定域性标准的区域非常广。随着数值增加，等值面开始逐渐减小。在0.52时，如红色箭头所指，恰好两个(3,+1)刚刚出现。它们的出现容易理解，因为(3,+1)代表这个点在两个方向上是函数值极小点，而另一个方向是极大点，由等值面随数值的变化过程图可以想象，那两个作为极小点的方向就是顺着等值面的两个正交方向，而作为极大点的方向是基本上垂直于等值面的方向。进一步增加等值面数值，等值面会进一步萎缩，并且很多地方等值面开始分离，联通的大域逐渐分解为小域。在0.66时，红色箭头所指的地方就是一个域开始分歧的地方，称为分歧点，凡是这种地方都对应于一个(3,-1)点，而拓扑分析时获得的相应(3,-1)点的ELF数值，也就正对应于ELF域在此处分歧时的数值。等值面数值增加到比较高的情况，如0.84，从图中可见ELF域已经充分分离了，对应每个成键区域以及氢的高定域性区域已经很清晰明朗，但是氧的两对儿孤对电子的高定域性区域还是没分开。继续增加到0.915时，发现两对儿孤对电子的域恰好分离开，红箭头指出了分歧点。如果进一步增加数值，则每个域都只会进一步缩小直至消失，但不会在这个过程中分歧成为更小的域，这种不能再分歧的域称为不可约域。每个不可约域里面都包含一个(3,-3)临界点，对应于域随着等值面数值增加越变越小最后消失的位置。

等值面图比较直观清楚，各个高定域性区域一目了然，通过不断调节等值面数值，还可以直接确定出(3,-1)和(3,-3)的位置和这些临界点的数值。但是等值面图不便于准确地定量分析，无法得到很精确的临界点的位置和函数数值。而且，这种找临界点的方法也比较累，尤其是通过观察等值面的变化来确定(3,+1)和(3,-3)临界点是很困难的。直接用拓扑分析虽然准确快速，但电子结构复杂时，临界点分布纷繁复杂，不能从图上直观地获知各个空间位置电子定域性的大小、高定域性区域的分布范围。所以，拓扑分析和等值面分析建议结合使用，前者用于精确定量分析，后者用于快速定性分析。另外，通过等值面图观察不可约域的位置，也有助于确定是否在拓扑分析中漏掉了一些(3,-3)临界点，如果某不可约域的范围内没找到(3,-3)，就应当在搜索设置中将球的中心设在那个区域，让大量在球体内随机分布的点作为初猜点来寻找那个(3,-3)。 返回小木虫查看更多