计算化学公社

标题: 透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍 [打印本页]

作者
Author: sobereva 时间: 2019-9-15 19:28
标题: 透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍

透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍

Thoroughly understanding the nature of hydrogen bonds and estimating hydrogen bond strength easily and reliably: An introduction to an important research article published in JCC in the 2019

文/Sobereva@北京科音 2019-Sep-15

0 前言

氢键是最重要的分子内和分子间相互作用类型之一，尽管已经有巨量文章对其从各个角度进行探讨，但即便到现在依然是一个研究的热点，尚有很多未探究充分、透彻的问题。2019年9月，北京科音自然科学研究中心（http://www.keinsci.com）连同伊朗Islamic Azad大学以及捷克科学院的研究者，发表了一篇名为Exploring Nature and Predicting Strength of Hydrogen Bonds: A Correlation Analysis Between Atoms-in-Molecules Descriptors, Binding Energy and SAPT theory的文章，见J. Comput. Chem., 40, 2868-2881 (2019) DOI: 10.1002/jcc.26068。此文章利用量子化学计算和波函数分析，充分讨论了氢键的本质，并提出了新的氢键分类和估计氢键强度的方法，文章内容非常具有普遍意义，内容丰富，值得关注。下面就依次对此文章的亮点和重点做简明直白的介绍，更详细内容请阅读原文。此文中的波函数分析部分全部使用Multiwfn（http://sobereva.com/multiwfn）程序完成。

文章下载： (, 下载次数 Times of downloads: 2849)

1 亮点1：构造了具有代表性的氢键作用体系

为了揭示氢键本质、讨论如何估计氢键强度，此文首先构造了一批氢键二聚体，如下所示

(, 下载次数 Times of downloads: 150)

这套氢键二聚体有几个特点：(1)样本数很大，使得后续分析讨论依据充分 (2)体系非常有代表性，常见的氢键类型X-H...Y (X,Y=O,N,F)都被纳入其中 (3)中性和带电荷氢键体系同时包含，结合能涵盖了从极弱到极强所有范围，因此可以提供关于氢键全面的视角 (4)样本都是由恰当选取的小分子构成，可以充分避免氢键以外的作用区域对讨论氢键时的干扰。以上体系都通过B3LYP-D3(BJ)/ma-TZVP级别做了几何优化，所得结构可以从文章的补充材料里获得。

如上表所示，此文通过高精度的CCSD(T)/jul-cc-pVTZ级别计算了结合能（Binding energy, BE），还通过PSI4程序通过高阶SAPT方法SAPT2+(3)δMP2结合aug-cc-pVTZ基组做了能量分解并同时得到结合能。文中数据体现这样的高阶SAPT方法不仅对中性体系结果和金标准CCSD(T)非常接近，即便对于带电荷体系也能得到合理结果，表明SAPT分析完全可以放心地用于带电荷氢键。在此文的补充材料里还证明，通过代入简单的线性回归公式，可以令CCSD(T)/jul-cc-pVTZ级别算的结合能很好地重现更准确但也昂贵得多的CCSD(T)/CBS（TZ->QZ外推）级别算的结合能。

考虑到本文构造的氢键数据集非常全面、有代表性，且提供了高精度结合能数据，因此可以在未来作为检验更廉价计算级别计算精度的测试集。

2 亮点2：通过能量分解充分展现了氢键的作用本质

SAPT是目前被使用最多、最被广为接受的能量分解方法，它能将结合能分解成不同的物理成分，从而更好地理解相互作用的本质。虽然SAPT早已被大量用于氢键体系的研究中，但是像此文这样使用SAPT非常系统、全面地考察大批量典型氢键的文章还非常少。文中对上述42个氢键二聚体分别给出了结合能当中的静电作用(Elst)、色散作用(Dis)、诱导(Ind)、交换互斥(Rep)部分各自的贡献，其中前三项对于氢键的形成起到正贡献，它们所占贡献的比例在下图给出

(, 下载次数 Times of downloads: 139)

图中前28个体系是中性氢键二聚体，后14个是带电荷氢键作用。图中从左到右总的排列顺序是氢键作用逐渐增强。由图可以一目了然发现，不同强度的氢键有着十分不同的本质。强度很低的氢键主要是通过静电吸引和色散吸引效果一起维持的。而中等强度氢键，比如水二聚体，静电吸引作用占绝对主导，而色散吸引和诱导作用一起起到辅助效果。非常强的氢键大多数是带电荷氢键，由图可见这种情况下静电和诱导作用共同起到关键性角色，而色散作用相对而言可以忽略不计。注意SAPT方法给出的“诱导项”本质上体现的是电荷转移、极化、轨道相互作用。

3 亮点3：一种新的氢键分类

根据SAPT分析数据和计算的氢键结合能，文中提出一种全新的氢键分类：

(, 下载次数 Times of downloads: 144)

这种分类方法相对于以往其他研究者提出的分类方法的关键好处在于，这种分类将氢键强度（由结合能BE体现）和氢键本质直接关联了起来。研究者们根据理论计算或实验观测到的氢键结合能，通过对照上表就能马上认识到这种氢键的主要本质大概是什么。由于此文考虑的样本较大，而且不是拘泥于特定类型体系，因此这个分类不仅可靠而且普适。再加上这个分类简单直白，所以很适合纳入教科书，对于学生正确理解氢键的本质很有好处。

4 亮点4：构建了氢键结合能与其能量成分之间的关系

按照常规逻辑，要想得到SAPT中的能量成分，肯定是要做SAPT能量分解的，但是研究者为此需要学习专门的程序的使用（PSI4或Molpro等），还需要付出计算代价。此文发现对于中性氢键作用体系，氢键结合能与各个能量成分间可以拟合出很好的线性关系，如下所示

(, 下载次数 Times of downloads: 159)

图中的式子中自变量BE是结合能。因此研究者只要通过一般方法计算出氢键结合能，直接代入上面的公式，就可以立马估计出各个SAPT能量成分了，省事极了！虽然不同的氢键体系都有一定个性，严格讨论时仍免不了需要直接做SAPT计算，但至少这套关系式足够给出定性正确的能量成分，用于粗略讨论足够了。注意用这套式子有个前提，就是你算的结合能必须确实只由某个氢键作用所主导，体系其它部分的作用相对而言可以忽略。

上面的图也体现出，氢键的所有能量项的绝对值都是随着氢键强度增加而增加的，对氢键形成产生正贡献的静电、色散、诱导作用增强时，起到负贡献的交换互斥作用也随之增加，很大程度上抵消掉结合作用。此外还可以看到，随着氢键增强，诱导作用增加的速度明显比色散作用增加的速度要快，即色散作用对氢键强度不是很敏感。对于强度只有<10 kJ/mol的弱氢键，色散作用不仅比诱导作用强，甚至和静电作用起到的地位不分伯仲，而当氢键强度进一步增加后，色散作用起到的地位就垫底了。

5 亮点5：证明了通过氢键的键临界点处的电子密度可以可靠地估计氢键强度

弱相互作用是波函数分析（wavefunction analysis）方法的常见研究对象。有很多研究者通过Atoms-in-molecules理论定义的键临界点(BCP)的属性（电子密度、电子密度的拉普拉斯函数、能量密度、势能密度等等）分析氢键的本质和估计氢键的强度。也有些人提出了一些指标，比如CVB指数、ΔΔVn指数来考察氢键强度。此前有一个知名的关系是Espinoa在Chem. Phys. Lett., 285, 170 (1998)提出来的BE=V(BCP)/2，即氢键的结合能约等于氢键的键临界点处势能密度的一半。然而，Espinoa的关系局限性非常大。首先，他当年考察氢键体系用的结合能的计算精度在如今来看非常差，而且他的文章只考虑了X-H...O (X=C,N,O)型氢键。

到底被使用广泛的Espinoa的估计氢键强度的公式准不准确？能不能找到更好的估计氢键强度的关系？这是一个极其重要的问题。因为，如果能确定一个简单且可靠的估计氢键强度的方法，那么就可以比较廉价地近似得到高精度量子化学方法计算的结果。更重要的是，利用这种关系可以估计一些不容易直接计算的氢键键能。比如说分子内氢键键能就不好算，因为这通常需要把体系进行截断、对截断处进行饱和，还得恰当对结构进行调节以避免严重位阻等问题，相当麻烦。再比如，有时候两个分子间形成多对氢键，一般的算分子间作用的方法只能给出总的作用能，想得到每个氢键各自的作用能很不容易实现。

此文将CCSD(T)算的42个氢键二聚体结合能与BCP处的电子密度ρ进行了线性拟合，如下所示

(, 下载次数 Times of downloads: 127)

由图可见，如果把带电荷和不带电荷的总共42个氢键体系直接拟合，虽然R^2看起来不差，但是MAPE（平均绝对百分误差）很大，因此这样拟合的关系在预测氢键作用能上没实际意义。然而，如上图(b)和(c)所示，若对中性和带电荷的氢键二聚体分别进行拟合，那么对二者都可以得到理想的线性关系，不仅R^2很高，同时MAPE很低，拟合出的公式很有实际意义。也就是说，我们只需要对氢键复合物先用量子化学程序进行几何结构优化，同时产生含有波函数的文件（如.wfn、.fch、.molden），然后用Multiwfn程序做电子密度的拓扑分析得到要考察的氢键的BCP处的电子密度（这个过程耗时几乎可忽略不计），再代入图中的公式，就可以直接得到较准确的氢键结合能。这样算出的结果准不准？假设有个中性体系的氢键的结合能是-5.0 kcal/mol（接近水二聚体的氢键强度），乘上MAPE=14.69%后可知，误差才0.7 kcal/mol而已，这样的误差实际上比起很多廉价方法（比如PM6-D3H4、PM7半经验方法）得到的氢键结合能准得多，很多时下算弱相互作用流行的DFT泛函也未必达得到这个精度。

为了考察是否能用BCP处的其它属性，或者一些与氢键有关的指标来准确预测氢键强度，文中也对它们与氢键结合能的关系做了回归分析，MAPE如下所示

(, 下载次数 Times of downloads: 135)

由表可见，用BCP处的能量密度（H）、bond degree（H/ρ）、势能密度的绝对值除以拉格朗日动能密度（|V|/G），都不能通过拟合关系很好预测氢键强度，至少比用电子密度来预测要差不少。Espinoa提出的用V/2来预测氢键结合能的做法也并不理想（对中性氢键MAPE达到23.35%），因此从此以后不应当再用这个关系估计氢键结合能。至于基于电子定域化函数定义的CVB指数，以及基于原子核位置静电势定义的ΔΔVn指数，虽然与氢键结合能确实有相关性，但从上面表格里的MAPE值来看也不算理想。因此此文有个发现值得强调：虽然BCP处的电子密度是最容易算的、形式最简单的量，但是基于它来预测氢键强度比用形式更复杂、看似更高级方式定义的量效果还好！

6 其它

在一些文章指出，SAPT的诱导项(E_ind)与共价作用程度有密切联系，甚至被一些文章直接拿来衡量共价作用程度。而在一些AIM研究文章里指出BCP处的|V|/G和H/ρ与共价作用也关系密切。考虑到这点，此文将|V|/G和H/ρ向E_ind进行了拟合。结果发现对于中性体系，|V|/G与E_ind之间可以通过二次函数很好地关联起来，暗示出BCP处|V|/G越大，氢键的共价成分越强，证明了将BCP处的|V|/G用于讨论不同中性体系的氢键的共价性是合理的。

想全面研究氢键本质，不仅要考虑平衡位置（优化后的复合物结构）下氢键的特征，还应当考察不同作用距离下氢键的特征。因此文章的补充材料里的S3一节里取了弱氢键、中等强度氢键以及带电荷的强氢键各一例，对SAPT总能量和各个分量随作用距离的变化情况进行了扫描和分析，如下所示

(, 下载次数 Times of downloads: 144)

由图明显可见不同的作用成份随氢键距离增加时衰减行为明显不同，详细分析讨论请看原文补充材料，在这里不再细谈。

7 总结

这篇JCC文章选取了非常有代表性的大批氢键体系，基于高精度的数据，从各个角度对氢键本质进行了深入的分析探讨，对于正确了解氢键本质非常有益，并且文中给出的氢键分类方法和估计氢键结合能的方法极具实用性，值得推广。这几点通过此文的Table of content图非常直观地体现了出来（其中电子密度ρ和结合能ΔE的单位分别用的是a.u.和kcal/mol）：

(, 下载次数 Times of downloads: 162)

8 相关阅读

本文提出的根据BCP处电子密度估计氢键键能的方法在Multiwfn手册4.2.1节以及《计算分子内氢键键能的几种方法》（http://sobereva.com/522）中给了实际例子，你会发现用起来非常简单，而且计算耗时极低。关于Multiwfn做AIM拓扑分析的更多信息，见《使用Multiwfn做拓扑分析以及计算孤对电子角度》（http://sobereva.com/108）和Multiwfn手册4.2节的各种例子。关于怎么用Gaussian等程序产生给Multiwfn用的输入文件，见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）。

上面提到的CVB指数的介绍和计算方法见：《使用Multiwfn计算CVB指数考察氢键强度》（http://sobereva.com/461）。

前述的ΔΔVn指数是怎么定义的，以及如何用Multiwfn计算，见Multiwfn手册4.1.2节。

用Multiwfn可以实现的各种研究弱相互作用的方法汇总见《Multiwfn支持的弱相互作用的分析方法概览》（http://sobereva.com/252），其中还对能量分解做了简介。

使用PSI4做SAPT计算，以及扫描各个SAPT能量项随作用距离的变化可参考《考察SAPT能量分解的能量项随分子二聚体间距变化的简单方法》（http://sobereva.com/469）。

分子间弱相互作用的本质简介，以及上文中提到的B3LYP-D3(BJ)，见《谈谈“计算时是否需要加DFT-D3色散校正？”》（http://sobereva.com/413）。

文中用到的ma-TZVP基组用法见《给def2以ma-方式加弥散函数的Gaussian格式的基组定义文件（含所有def2支持的元素）》（http://sobereva.com/509）。

文中提到的counterpoise校正的原理，以及为什么文中用一半Counterpoise校正，在此文有提及《谈谈BSSE校正与Gaussian对它的处理》（http://sobereva.com/46）。

作者
Author: pyscf 时间: 2019-9-15 22:39
这是个很有意义的工作，两点看法：
1. 42个体系数量很多但是谈不上全面. 2014年Weinhold提出了anti-electrostatic hydrogen bond在这个工作中一个都没有
比如一个简单的体系 [F...HOCO2]^2-
如果加上这种非典型的氢键，文章中的好几个拟合关系是否还成立？SAPT又会给出什么结论？

2. 这个工作只考虑了最简单的 A...B dimer类型的氢键体系并得到了很多拟合关系
但是对于稍微复杂的体系，比如水团簇结构，从简单体系得到的结论要用到这些复杂的体系，可不可行，文章中并没有给出可靠证据
实际的化学研究中，很少有人会只合成一个dimer来研究，更多的是复杂的结构 (比如蛋白)——dimer周围还有其他影响因素
因此一定程度上，这个工作有误导性

作者
Author: granvia 时间: 2019-9-15 23:50

pyscf 发表于 2019-9-15 22:39
这是个很有意义的工作，两点看法：
1. 42个体系数量很多但是谈不上全面. 2014年Weinhold提出了anti-elect ...

还有X-H...pi这类体系，不知道SAPT能否给出类似的结论？

作者
Author: beefly 时间: 2019-9-16 07:41
本帖最后由 beefly 于 2019-9-16 07:44 编辑

能量密度的正负有不同的物理意义，但是氢键分类里似乎未考虑。或者说，氢键临界点的能量密度是否应当都是负数（也许弱氢键除外）？

作者
Author: yjmaxpayne 时间: 2019-9-16 10:47
社长威武，朝着定量化分析氢键及其他弱相关作用又迈进了一大步。

作者
Author: yflchx 时间: 2019-9-16 11:04

granvia 发表于 2019-9-15 23:50
还有X-H...pi这类体系，不知道SAPT能否给出类似的结论？

X-H...π大概属于Weak to medium类（−2.5 to −14.0 kcal/mol），本质应该是Dispersion + Electrostatics类别。

作者
Author: ShuangfeiZhu 时间: 2019-9-16 11:07
原文的优化基组是ma-TZVPP，笔误？

作者
Author: sobereva 时间: 2019-9-17 02:59

pyscf 发表于 2019-9-15 22:39
这是个很有意义的工作，两点看法：
1. 42个体系数量很多但是谈不上全面. 2014年Weinhold提出了anti-elect ...

1 你若感兴趣可以测试

2 这不是文本的考察范畴。氢键类型极多，pi氢键、二氢键、共振辅助氢键、协同氢键等等，不可能一篇文章把所有氢键类型全都纳入其中。网罗各种氢键，属于大型的氢键综述文章考虑的范畴
如果各种类型氢键体系全都混杂到一起计算讨论，内容会相当繁复、文章冗长，反倒没法对于最常见、典型的氢键的本质得到清晰、明确的结论。相对于本文考察的体系而言，本文中的论述没有误导性

作者
Author: sobereva 时间: 2019-9-17 03:01

ShuangfeiZhu 发表于 2019-9-16 11:07
原文的优化基组是ma-TZVPP，笔误？

没笔误。有何笔误？

作者
Author: sobereva 时间: 2019-9-17 03:08

granvia 发表于 2019-9-15 23:50
还有X-H...pi这类体系，不知道SAPT能否给出类似的结论？

虽然没有纳入文章的分析（要不然还得分门别类讨论内容太多了），但预计文中的结论也适用。

作者
Author: sobereva 时间: 2019-9-17 03:12

beefly 发表于 2019-9-16 07:41
能量密度的正负有不同的物理意义，但是氢键分类里似乎未考虑。或者说，氢键临界点的能量密度是否应当都是负 ...

可以看原文的表5，里面的H(r)是考察的体系的BCP处的能量密度。可见较弱、中等偏弱的都是正值，后面更强的都是负值。
顺带一提，这也再次反映出一个常见现象，也就是用H(BCP)符号是否为负值判断是否是共价作用过于宽松了。

作者
Author: chemhou 时间: 2019-9-17 04:35
工作很有意义

作者
Author: granvia 时间: 2019-9-17 10:18

sobereva 发表于 2019-9-17 03:08
虽然没有纳入文章的分析（要不然还得分门别类讨论内容太多了），但预计文中的结论也适用。

这个还有待探讨吧？

作者
Author: granvia 时间: 2019-9-17 10:19

yflchx 发表于 2019-9-16 11:04
X-H...π大概属于Weak to medium类（2.5 to 14.0 kcal/mol），本质应该是Dispersion + Ele ...

有文献支持你后面一条论断吗？

作者
Author: yflchx 时间: 2019-9-17 12:16

granvia 发表于 2019-9-17 10:19
有文献支持你后面一条论断吗？

Jan Rezac, Kevin E. Riley, and Pavel Hobza, J. Chem. Theory Comput. 2011, 7, 2427–2438.

请查看Table 2 others部分。

作者
Author: Yaqi 时间: 2019-10-13 17:14
谢谢分享，对于理解弱相互作用很有帮助！

作者
Author: Yaqi 时间: 2020-1-10 19:41
请教卢老师：帖子中，第二幅能量分解的柱状图，是用什么作图软件实现这种层叠的柱状图效果？

作者
Author: sobereva 时间: 2020-1-10 22:12

Yaqi 发表于 2020-1-10 19:41
请教卢老师：帖子中，第二幅能量分解的柱状图，是用什么作图软件实现这种层叠的柱状图效果？

那个图是合作者Saeed画的，应该用的是Origin

作者
Author: Yaqi 时间: 2020-1-11 20:04

sobereva 发表于 2020-1-10 22:12
那个图是合作者Saeed画的，应该用的是Origin

好的，我试一下，谢谢老师

作者
Author: snljty 时间: 2020-2-27 23:47
请问比如氨基酸水溶液中两性离子，分子内羧基负离子和氨基阳离子之间的氢键，适合用这个BCP处电子密度的方法估计么，两边有明显电荷，用中性氢键体系预测公式是否合理呢？谢谢

作者
Author: sobereva 时间: 2020-2-28 03:05

snljty 发表于 2020-2-27 23:47
请问比如氨基酸水溶液中两性离子，分子内羧基负离子和氨基阳离子之间的氢键，适合用这个BCP处电子密度的方 ...

文中预测中性氢键键能的公式，对于氢键作用位置明确的（不显著牵扯二级作用等情况）的中性氢键很普适
但是你的情况，COO-上负电荷同时明显分布在两个氧上，比较复杂，原理上肯定不可能光靠一个位置的BCP靠谱地预测氢键作用。应当用下文里其它做法计算
计算分子内氢键键能的几种方法
http://sobereva.com/522

作者
Author: snljty 时间: 2020-2-28 08:25

sobereva 发表于 2020-2-28 03:05
文中预测中性氢键键能的公式，对于氢键作用位置明确的（不显著牵扯二级作用等情况）的中性氢键很普适
但 ...

了解，谢谢老师！

作者
Author: granvia 时间: 2020-7-30 00:06
有点觉得B3LYP-D3BJ优化很弱的氢键复合物（如结合能只有-0.5 kcal/mol左右）似乎未必可靠吧（毕竟 B3LYP不胜任弱相互作用，且DFT-D3属于经验方法）？由于体系不是很大，为什么不用CCSD(T)甚或MP2优化个别体系测试一下呢？

作者
Author: sobereva 时间: 2020-7-30 03:09

granvia 发表于 2020-7-30 00:06
有点觉得B3LYP-D3BJ优化很弱的氢键复合物（如结合能只有-0.5 kcal/mol左右）似乎未必可靠吧（毕竟 B3LYP不 ...

B3LYP-D3(BJ)算弱相互作用能精度已经不错，而优化对于计算级别要求明显更低，所以这级别就够了。我和合作者大量的计算经验、前人的文章都可以证明这个级别对于优化完全够用了。
另外更关键的是，当前文章给出的拟合关系是希望其他人能拿来解决实际问题、用于预测他们自己体系的作用能。若比如当前文章用CCSD/def2-TZVP做优化并且得到波函数信息，虽然名义上结构和电子密度质量会更准确一点点，但这导致拟合公式缺乏实用性，毕竟稍微大一点的体系拿这个级别都根本算不动。

作者
Author: granvia 时间: 2020-7-30 12:29

sobereva 发表于 2020-7-30 03:09
B3LYP-D3(BJ)算弱相互作用能精度已经不错，而优化对于计算级别要求明显更低，所以这级别就够了。我和合作 ...

我只是觉得对于特别弱的氢键体系，如结合能只有-0.6kcal/mol的H2O-CH4体系，DFT再精确也无法达到这样的精度吧。要知道，CCSD(T)的单点计算都是基于DFT势能面的基础上的，所以应该对于这类个别体系用CCSD(T)验证一下平衡构型是否可靠，毕竟小于1kcal的能量差都可视为柔性体系了，与分子轻微的相对转动能量相当了吧

作者
Author: sobereva 时间: 2020-7-31 01:23

granvia 发表于 2020-7-30 12:29
我只是觉得对于特别弱的氢键体系，如结合能只有-0.6kcal/mol的H2O-CH4体系，DFT再精确也无法达到这样的精 ...

B3LYP-D3(BJ)优化这个不会有什么问题
就算因为势能面平缓，可能对这个体系B3LYP-D3(BJ)结构误差比其它的稍微大一丁点，但由于本身势能面平缓使得高级别算的作用能对结构不敏感，所以用当前计算级别也没什么问题

作者
Author: granvia 时间: 2020-7-31 12:11

sobereva 发表于 2020-7-31 01:23
B3LYP-D3(BJ)优化这个不会有什么问题
就算因为势能面平缓，可能对这个体系B3LYP-D3(BJ)结构误差比其它的 ...

明白了，多谢解答

作者
Author: 牧生 时间: 2022-4-7 21:18
本帖最后由牧生于 2022-4-7 21:50 编辑

请教这样一个让我困惑的疑问，以具体实例进行问题描述：

1、以冰晶胞和水，少量NaCl和CO2，低温高压得到水合物以后，取某个水分子及附近的几个水分子为A体系，该体系净电荷为0

2、取CL原子附近的几个分子为B体系，该体系净电荷为-1

3、删掉虚原子，整理好原子顺序，使用orca得到波函数信息（r2scan-3c，极为简单的体系，应该足够）

4、A体系和B体系的IGMH图分别如下
(, 下载次数 Times of downloads: 70) (, 下载次数 Times of downloads: 81)

5、从上图中看出，A体系中，四个等值面看起来都一样的。而在B体系中，一个氯原子与周围5个水分子形成了5个氢键，其中三个等值面都是蓝色，为典型的氢键作用；另外两个氢键以绿色为主。那么有个疑问，这两个绿色的部分，由于分子间距离和角度符合氢键的特征，被归为氢键，但等值面又是绿色的，也可以说是范德华相互作用。
问题1：如果非要进行一个总结，这两个绿色等值面的地方，到底应该归结为范德华相互作用，还是氢键呢？

问题2：我该如何定量去描述这几个氢键（a~i）的强度，如何进行对比。

自我发挥的方法如下：
如果都按照体系为电中性计算，参考此文http://sobereva.com/522 套用公式，得到的结果如下

A (neutral)				B (charged)
a	b	c	d	e	f	g	h	i
-9.36	-6.61	-6.99	-10.1	-6.52	-6.77	-3.98	-2.26	-6.05

a~d的键能有所区别，是因为水分子有部分参与形成五边形，部分参与形成六边形，所以键能不一样，这个我能理解。但是对于e~i，虽然数据值看起来合理，但是由于体系带电，上表中直接套用不带电的公式得到的数值应该是不可靠的？？
如果套用带电的那个公式，那么得到的值则明显低的多，表明键的强度更大，如下表

B (charged)
e	f	g	h	i
-11.88970222	-12.25455046	-8.102392287	-5.531298989	-11.18637682

这就出现了我不能理解的矛盾。等值面颜色来看，蓝色的强度>绿色，但套用带电的那个公式，数值明显不靠谱。
如果用CVB，得到的结果是

neutral				charged
a	b	c	d	e	f	g	h	i
-0.075187	-0.033891	-0.032233	-0.059152	-0.068066	-0.069687	-0.016576	0.028026	-0.059147

归根结底，我该如何用数值，定量的去区分这几个氢键的强度呢？什么时候用CVB，什么时候用结合能呢？

作者
Author: sobereva 时间: 2022-4-8 16:02

牧生发表于 2022-4-7 21:18
请教这样一个让我困惑的疑问，以具体实例进行问题描述：

1、以冰晶胞和水，少量NaCl和CO2，低温高压得到 ...

g和i应当算氢键，通过结构和常识就能判断
几个氢键相对强弱可以用BCP的密度讨论
当前涉及到带电和不带电的情况对比，具体定量强度应当用高精度量子化学方法算相互作用能来得到和对比

作者
Author: 牧生 时间: 2022-4-10 15:37
本帖最后由牧生于 2022-4-10 15:48 编辑

sobereva 发表于 2022-4-8 16:02
g和i应当算氢键，通过结构和常识就能判断
几个氢键相对强弱可以用BCP的密度讨论
当前涉及到带电和不带 ...

重新进行了计算，使用orca得到输入文件。为了表明都是同一个级别，且为了照顾阴离子体系B，所以，都选用了-2加了弥散函数，直接选择了7b得到orca的输入文件，得到的结果好一些，具体如下：

neutral 直接套用不带电公式				charged 直接套用带电公式
a	b	c	d	e	f	g	h	i
-11.80013081	-7.454337147	-8.103973693	-9.889485485	-11.83344969	-12.17984292	-11.14618356	-8.002261914	-5.418346706

数值上看起来比之前合理多了，但IGMH图形上出现了我不好解释的地方，以体系A为例，如图所示
(, 下载次数 Times of downloads: 77)

同一个等值面，从外向里看都是蓝色，从里往外看有很明显的红圈，这个该如何解释呢？

以下为计算5个水分子的inp文件，请帮忙看一下，精度是否足够了。

! wB97X-2 D3 ma-def2-QZVPP autoaux RIJCOSX tightSCF noautostart miniprint nopop
%maxcore  1000
%pal nprocs 4 end
# DFT-D3(BJ) parameter for wB97X-2, see SI of PCCP, 20, 23175 (2018)
%method
D3S6 0.547
D3A1 3.520
D3S8 0.0
D3A2 7.795
end
* xyz 0 1
O    13.45000000 9.69000000 36.79000000
H    13.69000000 8.86000000 37.20000000
H    13.87000000 10.36000000 37.34000000
O    16.35000000 7.55000000 34.63000000
H    16.29000000 6.96000000 33.87000000
H    17.26000000 7.85000000 34.62000000
O    14.76000000 9.53000000 34.55000000
H    14.25000000 9.39000000 35.35000000
H    15.47000000 8.89000000 34.61000000
O    16.24000000 11.82000000 34.48000000
H    15.83000000 12.34000000 33.78000000
H    15.69000000 11.05000000 34.56000000
O    13.45000000 9.48000000 32.21000000
H    13.81000000 9.48000000 33.09000000
H    12.52000000 9.29000000 32.32000000
*

作者
Author: wzkchem5 时间: 2022-4-10 16:28

牧生发表于 2022-4-10 08:37
重新进行了计算，使用orca得到输入文件。为了表明都是同一个级别，且为了照顾阴离子体系B，所以，都选用 ...

虽然这个跟你观察到的IGMH图的问题估计关系不大，但是如果用双杂化泛函，建议考虑MP2项对密度的贡献（%mp2 density relaxed end；如果算不动，把relaxed改成unrelaxed，但精度会差一点。作图时注意不要用SCF密度，要用MP2密度），否则算出来的密度的精度可能还不如普通杂化泛函，因为漏掉了一些correlation的贡献

作者
Author: sobereva 时间: 2022-4-11 00:25

牧生发表于 2022-4-10 15:37
重新进行了计算，使用orca得到输入文件。为了表明都是同一个级别，且为了照顾阴离子体系B，所以，都选用 ...

因为等值面不是无限薄，而函数也不是处处相同，因此两个方向看会有所不同。如果把等值面数值设大一点，等值面会有所收缩，两个方向看的差异会更小。反正怎么看中间都是蓝色，就没比要刻意区分了。

wB97X-2是双杂化泛函，默认情况下产生的波函数是它的杂化泛函部分的，然而这还不如用B3LYP等普通泛函算得好。要么产生双杂化泛函的波函数（需要额外时间），要么用B3LYP的波函数，建议后者。因为我的JCC 2019年的文章里的参数是在B3LYP/ma-TZVP下拟合的，套用预测公式时波函数也应当在那个级别下计算。

作者
Author: 牧生 时间: 2022-4-11 07:22

sobereva 发表于 2022-4-11 00:25
因为等值面不是无限薄，而函数也不是处处相同，因此两个方向看会有所不同。如果把等值面数值设大一点，等 ...

OK，问题解决，感谢解答
进行了B3LYP/ma-TZVP和wB97X-2 D3 ma-def2-QZVPP得到的结果对比，带电体系，得到的数值几乎完全精确相同，差异不超过0.1%；对于不带电的体系，大约有1~2%的差异。

欢迎光临计算化学公社 (http://bbs.keinsci.com/)