使用Multiwfn计算CVB指数考察氢键强度
使用Multiwfn计算CVB指数考察氢键强度
Using Multiwfn to calculate CVB index and measure strength of hydrogen bonds
文/Sobereva@北京科音
First release: 2019-Feb-5 Last update: 2019-Feb-6
摘要：本文介绍基于电子定域化函数(ELF)定义的用于衡量氢键强度的CVB指数，并且演示怎么在Multiwfn程序（http://sobereva.com/multiwfn）中计算。掌握这个方法，在讨论氢键问题的时候相当于多了一个很有用的工具。
1 CVB指数的思想以及计算方法
core-valence bifurcation (CVB)指数是最初在Theor. Chem. Acc., 104,13 (2000)中由Silvi等人提出的一种考察氢键强度的方法，它利用了电子定域化函数(ELF)的拓扑分析的思想。如果对ELF不熟悉的话，请参阅《ELF综述和重要文献小合集》（http://bbs.keinsci.com/thread-2100-1-1.html）。后来有不少文章都使用了这个指数讨论氢键，比如Struct. Chem., 16, 203 (2005)、J. Phys. Chem. A, 115, 10078 (2011)，在一篇从理论角度分析氢键的综述Chem. Rev., 111, 2597 (2011)中也对CVB指数进行了介绍。
一般的氢键可以写为D-H...A的形式，其中D是氢键给体原子（donor），A是氢键受体原子（acceptor）。用ELF盆分析的语言来描述的话，这块区域由以下的ELF盆构成：
V(D,H)：D和与之成键的H共同构成的价层(valence)盆
C(D)和C(A)：D原子和A原子的内核(core)盆
V(A)：A原子的价层盆
相邻的盆之间都有所谓的二分点(bifurcation point)，对应于拓扑分析语言里的(3,-1)型ELF临界点。如果你不懂什么叫二分点、不懂怎么计算的话，可以看看《在Multiwfn中单独考察pi电子结构特征》（http://sobereva.com/432）一文中对ELF-pi函数考察其二分点的例子，本文第四节的例子也通过ELF拓扑分析考察了二分点。
CVB指数在Theor. Chem. Acc.那篇文章中的定义为：CVB index = ELF(C-V) - ELF(DH-A)
其中ELF(C-V)代表内核盆与价层盆之间的二分点数值，ELF(DH-A)代表V(D,H)与V(A)之间的二分点数值，也即H与A原子间(3,-1)型ELF临界点的ELF数值。CVB指数越负，通常氢键越强，这点在原文以后后续一些文章里都通过实例进行了验证。比较强的氢键的CVB指数一般明显为负；极强、带有明显共价特征的氢键可以达到很负的数值；中等强度氢键的CVB指数一般在0附近；较弱的氢键的CVB指数一般明显为正。其实为什么会这样也不难理解，因为氢键越强，H和A之间的距离通常越近，H...A作用的共价性也会越显著，势必ELF(DH-A)会越大、导致CVB index越负。其实说白了，CVB指数的名字虽然重在core-valence二分点数值，但真正关键的其实是ELF(DH-A)的数值，不同氢键的ELF(DH-A)的差异远比ELF(C-V)大得多得多。
在CVB指数的原文里，对CVB指数的计算细节没有说清楚，这导致之后不同文章在计算CVB指数的时候用的具体做法不同，而且CVB指数原文里的数据还存在错误，尤为明显的错误是把符号搞反了，而且在后来的某些相关文章中和CVB有关的叙述还往往是错的，总之被搞得一团糟，那些文章里的CVB数值也很难或根本无法重现出来。
笔者在Multiwfn里专门加入了CVB指数的自动计算功能，是主功能200的子功能1，只要输入D、H、A原子的序号即可计算。在Multiwfn中计算CVB指数的做法是笔者提出来的，笔者认为这种方式计算不仅耗时很低（基本上秒算）、操作极简单、实现极容易，而且意义明确。
在Multiwfn的这个功能中，CVB指数定义为：CVB index = ELF(C-V,D) - ELF(DH-A)
其中，ELF(DH-A)不是通过拓扑分析搜索出来的精确的V(D,H)与V(A)之间的二分点的数值，而是在H-A连线上计算ELF曲线，从中取H和A之间ELF极小点的数值，因为这样计算容易得多。而且由于精确的二分点位置几乎就在H-A的连线上，因此Multiwfn中这种近似方式得到的ELF(DH-A)与按照CVB指数原始定义计算的ELF(DH-A)几乎没有差异。
上式中ELF(C-V,D)的物理意义是D原子上的core-valence二分点数值。Multiwfn中虽然也可以通过其强大的拓扑分析功能（主功能2）得到准确的这种二分点位置和数值，但这种二分点在实际中往往不好考察，同一个原子可能有多个这样的二分点，在选取上有任意性。而且如果是第二周期之后的原子，二分点会更多，考察时更麻烦。而且笔者也认为，用精确的二分点也并没有太强的物理意义，因为它们往往都不在D-H连线方向上，因此和氢键特征也联系不紧密。考虑到这些，在Multiwfn自动计算CVB指数的功能中，ELF(C-V,D)是直接取D-H连线上相应的ELF极小点数值。
Multiwfn自动计算CVB指数的功能中会计算三个量，思想很简单，如下所示：
实际上你也完全可以通过主功能3来绘制ELF曲线（例子看Multiwfn手册4.3节），并让程序给出极小点位置和数值，由此手动得到CVB指数，只不过使用Multiwfn自动计算CVB指数的功能在操作上方便得多。上图中顺带搜索到的ELF(C-V,A)并没什么实际用处，但是可能有的用户对此感兴趣，所以程序也会顺带输出。
由于Theor. Chem. Acc.那篇文章里的数据有一定问题，而且对测试体系计算氢键结合能用的级别在现在来看非常烂（B3LYP/6-31G**），因此文中的数据不足为信。笔者用这篇文章里的一部分氢键二聚体，在可靠的B3LYP-D3(BJ)/def2-TZVP下优化并得到波函数文件，使用Multiwfn自动计算CVB指数的功能做了计算，并且用计算氢键结合能精度很高的MP2/jul-cc-pVQZ结合counterpoise校正计算了氢键结合能，结果如下
可见，CVB指数和氢键结合能之间的线性关系虽然不算完美，但线性关系还是很显著的。如果把OC...HF这个例外不计的话，那么线性关系就很理想了。
其实，CVB指数的两个组成部分自身也与氢键强度存在明显相关性（起码对这批体系而言），如下所示
由于“CVB指数 vs 结合能”的图和“ELF(DH-A) vs 结合能”很相似，所以CVB指数的重点是ELF(DH-A)，哪怕忽略ELF(C-V,D)其实也无大碍。ELF(C-V,D)与结合能的线性关系很好，但并不好从物理意义上解释原因，同时也可以看到对不同氢键其数值差异很小。总之，当你用CVB指数讨论氢键的时候，也可以顺带检验一下其两个组成部分和氢键强度的相关性。
值得一提的是，对这些氢键二聚体，笔者发现AIM理论里的键临界点(BCP)的电子密度ρ(BCP)和结合能的关系颇好，如下所示（笔者也考察了BCP处势能密度，发现与结合能的线性关系远没有这么好）。ρ(BCP)用Multiwfn可以非常容易地计算，见《使用Multiwfn做拓扑分析以及计算孤对电子角度》（http://sobereva.com/108）。AIM的相关知识看《AIM学习资料和重要文献合集》（http://bbs.keinsci.com/thread-362-1-1.html）。关于此项研究，笔者专门发表了论文，主要内容介绍见《透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍》（http://sobereva.com/513），此文对于氢键的研究意义很大，请读者务必一读。
2 CVB指数在Multiwfn中的计算一例：HF...HF
只有在2019-Feb-5及之后更新的Multiwfn中才有CVB指数自动计算功能。不了解Multiwfn的话看《Multiwfn入门tips》（http://sobereva.com/167），不知道怎么产生Multiwfn需要的输入文件的话看《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）。常见的.fch、.molden、.wfn等格式都可以作为CVB指数计算功能的输入文件。
下面以HF...HF作为例子演示在Multiwfn中计算CVB指数，由Gaussian 16 A.03在B3LYP-D3(BJ)/def2-TZVP下优化并产生的wfn文件是examples目录下的HF_HF.wfn，体系结构如下
启动Multiwfn，输入
examples\HF_HF.wfn
200 // 其它功能(part2)
1 // 计算CVB指数和相关的量
2,1,3 //氢键给体原子、氢原子和氢键受体原子的序号
输出信息如下
Core-valence bifurcation value at donor, ELF(C-V,D): 0.0936
Distance between corresponding minimum and the hydrogen: 0.743 Angstrom
Core-valence bifurcation value at acceptor, ELF(C-V,A): 0.1408
Distance between corresponding minimum and the hydrogen: 1.628 Angstrom
Bifurcation value at H-bond, ELF(DH-A): 0.0648
Distance between corresponding minimum and the hydrogen: 0.614 Angstrom
The CVB index, namely ELF(C-V,D) - ELF(DH-A): 0.028768
各项是什么含义，只要仔细看了上一节就肯定明白。输出中还顺便把相应的ELF曲线极小点与氢原子之间的距离显示了出来。
3 某些极强氢键的特殊考虑
对于某些极强的氢键，由于H和A之间的作用也非常强，存在显著共价特征，此时V(D,H)二分为V(D)和V(H)，例如下图所示的体系（由于O-H-F基本在一条直线上，所以就不把O-H和H-F的ELF曲线图单独绘制了，而是为了省事直接绘制O-F连线上的图）。
此体系中，由于F形式上是以F-阴离子形态存在，与带正电的氢的作用极强，这使H和F之间已经存在不可忽视的共价性（而一般氢键的主要本质是静电作用），因此在对应的二分点C位置，ELF数值远大于一般强度氢键的情况。H与F强烈吸引，导致水分子中的O3-H2键遭到很大破坏，拉长了O3-H2的距离，也同时导致V(H2)独立的ELF盆的出现，因此曲线上出现了极小点A。
像上面这种明显不具备典型D-H...A氢键特征的体系，前述的Multiwfn的CVB指数自动计算方法是明显不适用的，强行用的话结果是没意义的，必须通过主功能3作图并手动计算。这个体系可以认为存在两种氢键O-H...F和O...H-F，二者的结合能存在明显差异，其CVB指数在计算的时候也应当分别考虑：
考察O-H...F时，CVB index = ELF(B) - ELF(C)
考察O...H-F时，CVB index = ELF(D) - ELF(A)
由于ELF(B)和ELF(D)差不太多，而ELF(A)显著大于ELF(C)，因此O...H-F的CVB指数明显更负，体现O...H-F的结合强度比O-H...F肯定大得多。
在Multiwfn手册的3.200.1节当中还有关于CVB指数更多的信息，有兴趣的读者建议看看。
4 氢键受体不是单一原子情况时CVB指数的计算
有很多氢键的受体并不是特定的某个原子，比如pi-氢键，是以富集的pi电子云作为氢键受体。这种情况下，也没法用前述的Multiwfn中自动计算CVB指数的功能来计算，而需要根据实际情况手动进行计算。这里以HF...乙烯二聚体作为例子。B3LYP-D3(BJ)/def2-TZVP下优化并产生的波函数文件是examples目录下的C2H4_HF.wfn。其结构如下
先计算ELF(C-V,D)。按照之前用到的做法，在F7-H8连线方向上绘制ELF，查看对应的极小点即可。启动Multiwfn并输入
examples\C2H4_HF.wfn
3 //绘制曲线图
9 //绘制ELF
0 //设定延展距离
0 //延展距离设为0（使得图的两端正好是两个原子核位置）
1 //输入两个原子序号定义直线
7,8
把弹出来的图关闭，后处理菜单选择6，在屏幕上看到
Minimum X (Bohr): 0.356648 Value: 0.94377038E-01
这个极小点的数值0.0944即是ELF(C-V,D)。
对于当前体系的ELF(DH-A)，我们应当通过做ELF拓扑分析找到对应H8与乙烯pi区域之间的ELF二分点位置来得到其数值。接着在Multiwfn里输入
0 //返回主菜单
2 //拓扑分析
-11 //选择要做拓扑分析的函数
9 //ELF
6 //在特定圆球内随机撒初猜点的方式找临界点
4 //将圆球中心定义为三个原子的中点
1,4,8 //将C1,C4,H8的中心作为圆球的中点，这时圆球的区域内理应会包含当前要找的氢键处的ELF二分点
0 //开始在圆球内随机撒初猜点来找临界点（撒的点数可以事先通过选项11设定，圆球半径也可以自设）
-9 //返回拓扑分析主界面
0 //观看临界点分布
看到的图如下所示，为了清楚，这里只显示出(3,-1)和(3,-3)型ELF临界点，分别是桔球和紫球。
我们看到在H8与pi电子区域之间找到了ELF二分点，即5号临界点。点右上角的RETURN退出图形窗口，选择选项7考察临界点属性，再输入5，在屏幕上看到
Electron localization function (ELF): 0.1240943921E+00
即此体系氢键的ELF(DH-A)=0.1241。因此，此氢键的CVB指数为0.0944-0.1241=-0.0297。
由于这个体系对称性较高，其实不做拓扑分析也可以得到ELF(DH-A)，也就是用主功能3绘制H8与C1-C4的中点连线上的ELF曲线，然后从中取对应的极小点，数值和上面做拓扑分析得到的完全一致。但如果体系对称性不那么高，显然就没法这么通过ELF曲线图比较准确地得到实际二分点位置的ELF值了。