使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用
PerformingIGMHanalysisbyMultiwfntoshowinteractionsinchemicalsystemsveryclearlyandintuitively
Firstreleas:2022-Mar-7Lastupdate:2024-Oct-29
下文将对IGMH的各个方面进行介绍，包括以下内容第1节：对IGMH的基本原理做简要说明，使读者快速了解此方法的思想第2节：给出IGMH用于各类体系的效果展示，使读者明白IGMH分析都能以什么方式研究什么问题第3节：将IGMH和IRI、NCI、IGM方法进行对比，便于读者了解什么时候适合用什么图形化分析方法第4节：对Multiwfn中的IGMH功能和所需的输入文件等方面进行说明第5节：给出使用Multiwfn做IGMH分析的具体例子第6节：对本文内容做总结
IGMH的详细介绍请读者仔细阅读前述的IGMH方法的原文，本节只是把关键信息做简要说明。
IGMH的基本思想和其前身IGM是共通的，都是使用三维函数δg来把原子间相互作用区域直观展现出来。为了便于说明，我们首先看一个最简单的情况，氢气分子H2。在H-H键轴上，两个氢的电子密度曲线图如下图(a)所示
上图(b)中的黑色曲线对应的g代表gradient，是两个H原子电子密度梯度之和。在两个H正中央g精确为0，这是因为在此处两个H的电子密度梯度大小相同但符号相反，正好抵消了。上图(b)中的绿色曲线对应的g_IGM是IGM（independentgradientmodel，独立梯度模型）型电子密度梯度，相当于两个H原子的电子密度梯度的绝对值加和，因此不会因为符号差异而出现相互抵消，仿佛两个原子是彼此“独立”的。显然g_IGM是g的上限。上图(b)中的蓝色曲线是δg，对应于g_IGM减去g。由图可见在差不多X=1.2埃处，δg达到了极大点。δg在原子间相互作用区域能有这么个峰出现，来自于两个原子的密度存在重叠，也暗示原子间相互作用的存在。在分子的平衡结构下，原子间密度重叠越大，这个峰往往越高，也暗示出原子间相互作用越强因而使得原子间距离越近。如果俩原子间都没有可查觉的相互作用，g和g_IGM就不会有差异，δg也就不会出现峰。
上图(b)中画了一个Y=0.2的虚线，和δg的交点用黑色箭头标注了。这体现出，对这个体系，若把δg绘制成等值面图，在0.2等值面数值（isovalue）的情况下就会出现一个等值面，包裹的区域就是上图中两个箭头之间的区域。可见δg函数的等值面可以用来把原子间存在明显相互作用的区域图形化展示出来，而且等值面数值取得越小，越弱的相互作用区域就越可能也出现相应等值面。而等值面数值取得较大时，则只会有较强的相互作用区域被展现出来。
上面给出的情况是一维情形，实际研究的当然都是三维情况。对于三维情况以如下方式定义g、g_IGM和δg。粗体r是坐标矢量，i循环所有原子，||符号代表取里面矢量的模，▽是矢量微分算符，ρ_i代表i原子的电子密度。
δg可以展现当前体系中所有原子间的相互作用。为了能专注于展现特定片段间（interfragment）和片段内（intrafragment）的相互作用，可分别定义δg_inter和δg_intra，如下所示。其中A循环用户定义的各个片段
IGM和IGMH方法不仅能展现相互作用区域，还能将sign(λ2)ρ函数通过不同颜色投影到δg、δg_inter、δg_intra的等值面上展现相互作用类型和强度。λ2是电子密度Hessian矩阵第二大本征值，sign()是取符号的意思。相互作用区域里某个位置的λ2若小于0体现出在此处原子间有吸引作用，大于0则体现互斥作用（但这只是个对多数相互作用区域适用的粗糙的、经验性的判据，切勿教条化理解）。ρ是当前体系实际电子密度，在相互作用区域此值越大暗示相互作用越强。sign(λ2)和ρ相乘得到的sign(λ2)ρ显然有能力对相互作用类型和强度进行区分。在IGMH原文里笔者给出了下图，可以作为IGMH图中sign(λ2)ρ的标准着色方式，不同颜色对应的常规解释也都给出了
为了解决IGM的这些问题，IGMH使用Hirshfeld原子空间划分的方式从体系实际的电子密度中划分出各个原子的电子密度，由于考虑了实际电子结构和化学环境，从物理意义上就比IGM强。这个改进还同时使得IGM等值面太肥厚的问题被理想地解决，IGMH方法的δg/inter/intra等值面总是比较薄且美观，能够纯粹地展现出相互作用区域的特征。IGMH的实现细节在IGMH原文2.2节有具体说明，本文就不再过多叙述了。IGMH和IGM还有个差异在于，IGMH分析中用的sign(λ2)ρ是基于实际电子密度计算的，而IGM分析中的sign(λ2)ρ是基于粗糙的准分子（promolecular）密度算的，准分子密度只是各个原子孤立状态电子密度的简单叠加，没有体现原子间相互作用对电子密度的改变。
由于IGMH要基于实际电子密度做分析，因此输入文件必须包含波函数信息，这需要做量子化学或第一性原理计算得到。另外，IGMH计算形式比IGM更为复杂，计算量更大。显然，做IGMH分析的总耗时是明显高于IGM的。不过以如今主流的几十核双路服务器的计算能力来说，IGMH分析用于几百个原子体系完全没问题。所以，除非真是碰到巨大体系用IGMH确实算不动，否则*强烈*不建议用IGM，而总应当用IGMH。实际上，整个IGM/IGMH分析流程中真正最耗时的是体系的几何优化（一般都是用DFT方法），IGMH和IGM分析本身的耗时通常只占几何优化的零头。从这个角度来说，IGM和IGMH分析能处理的体系尺度并没有实质性区别（除非你用很便宜的分子力场或半经验层面的方法做几何优化，此时的IGM、IGMH分析的准确性无疑会打一定折扣）。
为了能让IGMH在分析片段间相互作用时也能做一些定量分析讨论，笔者定义了一些指数，它们对于IGM也是同样适用的。
原子对δg指数（atomicpairδgindex）记为δG_pair，定义如下。其中A和B是要被考察相互作用的两个片段，i和j分别是两个片段中的原子
原子对的δG_pair(%)定义为它的δG_pair占两个片段间所有δG_pair总和的百分比，可以粗略体现原子对对片段间相互作用的贡献程度。
原子δg指数记为δG_atom，它可以衡量各个原子对片段间作用的贡献程度，其定义如下，其中i是一个片段里的原子，j循环另一个片段里的所有原子。相应地还可以定义δG_atom(%)用于衡量某原子对片段间作用贡献的百分比。
这一节展示一些IGMH图，以令读者对IGMH有直观的认识、充分了解其重要实用价值。图片都来自IGMH原文，原文里面有详细的说明和讨论，这一节里笔者对每个例子只是简单地提几句。这些图用的sign(λ2)ρ的色彩刻度基本都和前文1.1节给出的一致，有些图为了颜色更便于观察，将色彩刻度下限和上限做了适当的调节。为了图像效果尽可能好，这些图的等值面数值都经过反复调节，选取了最能充分、清楚展现感兴趣的相互作用的值。所有结构都用恰当的级别做过几何优化。
下图是一个三重联锁共价有机笼体系，由两个四面体笼嵌套构成，要想分离两个单体需要破坏共价键才行。IGMH分析中将两个单体各做为一个片段。下面(a)图是IGMH着色等值面图，(b)图是将原子按照δG_atom进行着色的图。可见IGMH将两个单体之间色散主导的作用区域，特别是两个分子中央的pi-pi堆积区域充分展现了出来，这也验证了这个体系的合成者在其原文里基于几何结构对pi-pi堆积存在性所作的推断。显然，光从几何结构上凭经验和感觉去分析分子间相互作用，远不如做个IGMH分析直观、严格、有说服力。从δG_atom着色图上可以进一步看出，对两个单体相互作用贡献最大的就是pi-pi堆积涉及的那些原子，其次是周围的原子，而整个体系顶端和底端蓝色区域原子的δG_atom几乎为0，因此对结合的重要性相对很低，这也是因为它们距离另一个单体原子相对较远。
IGMH分析不仅可以用于分子，也可以用于周期性体系。例如下图是Multiwfn基于CP2K第一性原理程序对冰计算的周期性波函数得到的IGMH图，任意选取的一个水和其它的水分别被定义为两个片段。可见中间的水与周围的水之间的四个氢键被蓝色的等值面展现得非常清楚。
下图通过IGMH展现了沸石与其中吸附的甲苯之间的相互作用，沸石和甲苯各定义为一个片段，结构使用CP2K做优化得到。由图可见甲苯与周围的沸石原子有大面积的相互作用。由于等值面都是绿色，说明吸附是色散作用主导的物理吸附。
下图通过IGMH展现了一种多孔共价有机框架化合物（COF）的层间相互作用。可见两层之间有无限延展的大面积的绿色等值面，无疑说明这种COF层间靠广阔的pi-pi堆积作用相结合。由图也可以看到等值面在晶胞边缘处被合理地截断，也体现出Multiwfn中的IGMH可以完美地考虑周期边界条件。
前面的例子主要是通过IGMH展现弱相互作用，而实际上IGMH对于很强的，甚至化学键程度的作用也同样可以很好地展现。例如下图是冠醚衍生物结合锂离子形成的体系，IGMH分析时将锂离子和冠醚定义为两个片段，由图可见确实IGMH图将锂离子与周围带负电的氧之间的高度离子性相互作用区域很好地展现了出来。另外，此图里还把键临界点和键径一起画了出来，每个临界点的电子密度都在图上给出了，通过这些数据便于定量讨论不同的氧与锂离子作用强度的差异。这有明显益处，毕竟光是从等值面颜色上判断容易有视觉误差、颜色差异较小时难以从视觉上区分，而且看到的颜色还多多少少受可视化程序里光源设置的影响。
下图是P4和Li6两种原子团簇结构。这回我们不划分片段，直接用sign(λ2)ρ着色的δg函数来展现体系中的所有相互作用。从下图左上角的IGMH方法的δg图上可以看到P-P化学键区域（等值面为蓝色的部分），以及三个P之间的环形位阻区域和四个P之间的笼形位阻区域（等值面为红色的部分），都被正确地展现了出来。IGM方法的δg函数的图在下图也给出了，可见效果较烂，位阻作用区域描述得乱七八糟，而且本来四面体笼状结构有很强的笼张力，成键区域的电子理应往外凸才对，结果图中描述化学键的蓝色等值面主体出现区域却反倒相对于P-P键轴往里收缩，可见IGM方法在此例表现得很差。下图右上角是P4的IRI图，可见它虽然和IGMH的δg图有异，但也把各种特征作用区域理想地展现了出来，而且相对来说IRI的图还显得更简单干净一些。
IGMH：如果你的目就是分析特定片段间或片段内的相互作用，很明确地知道片段该怎么定义，那么IGMH绝对是最理想的选择。从上一节的大量例子可看出，IGMH对各类体系表现得都很理想，图像很美观。
如果你之前已经会用Multiwfn做IGM分析，那么IGMH都不用学，因为IGMH和IGM分析操作上的差异仅仅在于进入主功能20之后选主功能11还是10，以及输入文件是否含有波函数信息，其它地方都是完全相同的。
如果体系很大，用DFT优化结构太耗时，用xtb等程序通过很廉价的半经验层面的GFN-xTB方法做优化也可以。xtb产生的molden文件也可以用于IGMH分析，结果大多数情况可以接受，但是由于GFN-xTB级别的波函数质量还是明显比DFT要差，基于它做的IGMH图和基于DFT波函数得到的图还是有可查觉的差距。所以最好还是基于像样的半经验方法或GFN-xTB方法优化的结构，用DFT算个单点任务得到用于IGMH分析的波函数文件。
常有人问怎么IGMH图作完了发现上面没有等值面，这大概率是因为片段间相互作用很弱，而默认的等值面数值不够小所致的。把等值面的数值手动改小到比如0.005甚至0.003，肯定能看见。如果此时还看不见等值面，那就说明片段间的相互作用可以完全忽略了，片段间的距离应该很远。
这里说的簇模型是指把晶体中的一个分子和与之有直接接触的一圈分子提取出来作为一个簇，中心分子和周围分子分别定义为IGMH分析的两个片段，这样分析结果就能完整展现晶体中的分子与周围环境分子的相互作用了。虽然基于CP2K计算产生的周期性波函数也能做分子晶体的IGMH分析，但使用簇模型来做有几个优点：(1)计算只考虑中心分子和邻近的一圈分子，可避免计算更外围、和中心分子没直接作用的分子导致浪费计算量。而且由于此时的IGMH分析过程中不需要考虑周期性，计算耗时比考虑周期性时低得多(2)避免更外围分子的结构显示在图中令图像混乱（虽然作为周期性体系分析时也可以在VMD中恰当设置来避免显示无关分子，但要多花费点步骤）(3)不需要用第一性原理程序做周期性计算，而只需要用量子化学程序就可以做，因此对量子化学研究者们更方便
Multiwfn瞬间就构造出了团簇，从屏幕上的提示可看到这个簇有88个原子，并且屏幕上还巨贴心地把中心分子中的原子序号给了出来，当前为1-3,25,32-34,69。把这个序号记下来，之后IGMH分析时要用到。
现在可以选菜单中的选项0看一眼新构造出的团簇是什么样，如下所示，可见非常理想，确实是中心分子被周围一圈分子所围绕
点图形窗口右上角的RETURN按钮关闭窗口，然后选选项“-3Outputsystemto.gjffile”，再输入cluster.gjf，当前团簇结构就保存为了当前目录下的cluster.gjf，这是Gaussian的输入文件格式。Multiwfn现在可以关了。
%chk=cluster.chk#PB3LYP/6-31G**em=GD3BJopt=readopt[空行]GeneratedbyMultiwfn[空行]01[坐标部分略...][空行]noatomsatoms=H[空行][空行]
下面就是做IGMH计算了。启动Multiwfn，然后输入cluster.fch20//弱相互作用可视化分析11//IGMH2//定义两个片段1-3,25,32-34,69//第一个片段里的序号，此例对应中心分子序号c//所有其余原子作为第二个片段，等同于输入4-24,26-31,35-68,70-8811//格点设置方式为：选中一批原子，在其四周扩展一定距离定义盒子（即格点分布区域），并指定格点间距。这种模式设置格点最适合当前情况1-3,25,32-34,69//盒子所围绕的中心分子中的原子序号2A//在中心分子四周扩展2埃定义盒子，通常这足够大了（如果最终从IGMH图中发现有些等值面被截断，可以再加大）0.15//格点间距（Bohr）。越小等值面越光滑，而计算耗时越高。鉴于当前体系较小，耗时必定不高，所以用了这样很精细的格点以得到尽可能理想的图像
算完后在后处理菜单选择3将格点数据导出成当前目录下的cub文件，然后将其中的dg_inter.cub和sl2r.cub挪到VMD目录下，也把Multiwfn文件包里的examples目录下的IGM_inter.vmd脚本挪到VMD目录下。启动VMD，然后输入sourceIGM_inter.vmd，即看到下图中δg_inter为0.01的图，此图中等值面上的蓝色部分将所有的氢键作用区域都很好地展现了出来。若进入Graphics-Representation，点击其中第二项，把Isovalue从0.1改小为0.004，看到的就是下图右侧的图，可见此图中绿色扁片部分把相对较弱的范德华作用区域也给展现了出来。
由此例可见，依靠Multiwfn以簇模型方式考察分子晶体中的弱相互作用特别便利。之后大家还可以进一步调整图像效果，比如把中心分子和周围分子用不同显示方式进行区分以凸显中间分子，还可以把键临界点和键径也画上去，如下所示
本文对笔者提出的IGMH方法的原理和不少细节做了介绍，给出了大量应用范例，充分体现出了IGMH的图像效果十分理想，在分析化学体系中的相互作用方面具有极高的实用性和很好的普适性。IGMH和IRI、NCI、IGM等图形化相互作用分析方法之间的差异和选用也在文中做了讨论。本文还介绍了强大的Multiwfn波函数分析程序中的IGMH分析的功能，并给出了具体的做IGMH分析的操作示例。可以看到IGMH分析在Multiwfn中使用特别方便，步骤简单，很容易上手。由于IGMH拥有诸多优点，又有不可替代的重要价值，并且有高效、易用的Multiwfn可以对各类体系实现IGMH分析的程序，IGMH显然会非常流行。