The way of plotting electrostatic potential for local region of molecular surface using Multiwfn in combination with VMD

文/Sobereva@北京科音  2025-Aug-2

下面是笔者登记了的第19001-22000篇使用了极为流行的Multiwfn波函数分析程序（//www.umsyar.com/multiwfn）发表的文章的pdf合集。pdf都是我手动下载的，由于疏忽漏了十几篇文章对应的pdf文件。
注：根据Google学术统计，Multiwfn的引用次数目前已经达到约35000。由于Google学术向我推送的信息不完整，因此我迄今提供的pdf合集里的并不是目前全部引用了Multiwfn的文章。

下载链接：https://pan.baidu.com/s/17NWaECclsnTnY7rih22ebw
提取码：fktj
文件使用分卷压缩，共14.2 GB，都下载后对任意一个压缩包进行解压即可。

本次汇总的文章的列表见//www.umsyar.com/attach/747/19001-22000.txt，或者看//www.umsyar.com/multiwfn/all_citation4.html和//www.umsyar.com/multiwfn/all_citation5.html里面的第19001~22000号文章。更早的文章列表看Multiwfn主页（//www.umsyar.com/multiwfn）的Cited by页面里的其它页面。有的文章标题太长，为了避免文件路径超过Windows限制，网盘里分享的pdf的文件名只保留了前面一部分。

这些文章可以视为Multiwfn的例子库，便于用户了解Multiwfn在 、第一性原理、分子模拟实际研究中能发挥的重要价值。如果想找某一类应用实例，可以利用Acrobat对指定文件夹里全部文档进行搜索的功能搜索相关关键词。

由于本文件包里的文章都是它们刚在期刊网站上刊登的时候我就下载的，所以文档中大多没有页码、卷号，且多数都是还没经过proofing的原始状态，大家只要到Google学术搜索里搜索标题就能找到文章最终版链接。

之前的引用Multiwfn的第1~1010篇文章可以从这里下载：http://bbs.keinsci.com/thread-2850-1-1.html
之前的引用Multiwfn的第1011~2001篇文章可以从这里下载：http://bbs.keinsci.com/thread-6754-1-1.html
之前的引用Multiwfn的第2002~3000篇文章可以从这里下载：http://bbs.keinsci.com/thread-11023-1-1.html
之前的引用Multiwfn的第3001~4000篇文章可以从这里下载：http://bbs.keinsci.com/thread-14181-1-1.html
之前的引用Multiwfn的第4001~5000篇文章可以从这里下载：http://bbs.keinsci.com/thread-16840-1-1.html
之前的引用Multiwfn的第5001~6000篇文章可以从这里下载：http://bbs.keinsci.com/thread-19818-1-1.html
之前的引用Multiwfn的第6001~7000篇文章可以从这里下载：http://bbs.keinsci.com/thread-22597-1-1.html
之前的引用Multiwfn的第7001~8000篇文章可以从这里下载：http://bbs.keinsci.com/thread-25497-1-1.html
之前的引用Multiwfn的第8001~9000篇文章可以从这里下载：http://bbs.keinsci.com/thread-27978-1-1.html
之前的引用Multiwfn的第9001~10000篇文章可以从这里下载：http://bbs.keinsci.com/thread-31224-1-1.html
之前的引用Multiwfn的第10001~11000篇文章可以从这里下载：http://bbs.keinsci.com/thread-34124-1-1.html
之前的引用Multiwfn的第11001~12000篇文章可以从这里下载：http://bbs.keinsci.com/thread-36387-1-1.html
之前的引用Multiwfn的第12001~13000篇文章可以从这里下载：http://bbs.keinsci.com/thread-38881-1-1.html
之前的引用Multiwfn的第13001~14000篇文章可以从这里下载：http://bbs.keinsci.com/thread-41495-1-1.html
之前的引用Multiwfn的第14001~16000篇文章可以从这里下载：
http://bbs.keinsci.com/thread-46062-1-1.html
之前的引用Multiwfn的第16001~19000篇文章可以从这里下载：
http://bbs.keinsci.com/thread-50824-1-1.html

再次顺带强调关于恰当引用Multiwfn的问题。在Multiwfn程序主页的首页和Download页面、程序手册、程序启动时屏幕上显示的信息、程序包内的诸多文件（Multiwfn quick start.pdf文档、How to cite Multiwfn.pdf文档、LICENSE.txt文件）、Multiwfn入门贴（//www.umsyar.com/167）、Multiwfn FAQ（//www.umsyar.com/452）等所有用户可以轻易看到的所有地方都非常非常非常明确强调如何对Multiwfn进行正确的引用，但使用Multiwfn的文章中不引用或者胡乱引用的现象依然十分严重（特别是中国、印度、巴基斯坦、伊朗作者的文章尤甚）！第19001~22000篇文章里没按要求恰当引用Multiwfn的文章多达555篇！经常是只提及Multiwfn但不引用，甚至居然在引用Multiwfn的地方引用的是和Multiwfn根本没有任何直接联系的其他人的文章（极为恶劣！）还有不少文章里作者明明在研究中用了Multiwfn做分析，在论文里竟然连Multiwfn的名字都不提。不按明确声明的方式恰当引用，对免费而且没有任何经费支持的程序的发展十分不利，同时也是严重缺乏学术道德的行为。借此机会再次强烈呼吁用户重视对Multiwfn的恰当引用。最恰当的引用方式见Multiwfn可执行文件包内的How to cite Multiwfn.pdf文档的说明。没有恰当引用Multiwfn的文章在前述的引用列表中都已经明确注上了Multiwfn was not properly cited!的字样进行公示，并且在本文的文件包里放在了“未恰当引用Multiwfn的文章”子目录中。

文/Sobereva@北京科音  2025-Jun-26

原子、分子间的相互作用在化学体系中无处不在。化学体系中的相互作用可以按强度大体分为化学键相互作用和弱相互作用，也可以按作用特征大体分为共价和非共价相互作用。这些相互作用普遍在数值层面上进行研究，如研究相互作用能及其物理成份（能量分解）、几何参数、键级、AIM拓扑分析考察键临界点属性，等等。而如今，这些相互作用的可视化分析（visual analysis）日趋流行，在计算化学研究文章中被应用得越来越普遍、在文献中越来越常见。这些方法可以令化学家们十分直观生动地认识相互作用出现的位置、强度、本质，不仅应用十分便利，也对传统的基于数值的研究相互作用的方式提供了关键性的互补。这些方法在近年来还得到了快速发展，涌现了许多新的分析手段。显然，对这些可视化分析方法进行全面、系统的综述是极为重要且迫切的。

近期，北京科音自然科学研究中心（http://www.keinsci.com）的卢天应知名的Angew. Chem. Int. Ed.期刊的邀请，发表了名为Visualization Analysis of Covalent and Noncovalent Interactions in Real Space的综述，可在此访问：https://onlinelibrary.wiley.com/doi/10.1002/anie.202504895、https://pan.baidu.com/s/17x4Pm0LINmt9iHDBLXWhhw?pwd=86rv。若未订阅，也可以免费阅读此文在ChemRxiv上的预印版https://doi.org/10.26434/chemrxiv-2025-9t442（内容和正式版基本一样，写文章时请引用正式版）。

推荐使用Multiwfn程序做文中介绍的各种可视化分析时引用此综述！

此综述是到目前为止最完整和系统介绍各种化学体系中相互作用可视化分析方法的文章，面向广大化学家，特别是计算化学家。文章的写作注重令内容清晰易懂、深入浅出，且引发广泛读者们的兴趣。通过此文，化学研究者们可以快速且充分了解各种可视化分析方法的基本思想、适用范畴，并且结合文中提供的丰富的例子，读者能够充分体会到它们的重要实用价值。在看过本文有了基本的知识框架后，若想进一步深究物理、算法细节等方面的知识以及了解更多实际应用，推荐再阅读此综述中引用的各种方法的原文和相关文章。

此综述介绍的可视化分析方法十分全面，并且囊括了最新进展。主要介绍了NCI（亦称RDG）、IRI、IRI-pi、IGM、IGMH、mIGM、aNCI、amIGM、Hirshfeld表面分析、静电势分析、范德华势分析、ELF、LOL、LOL-pi、电子密度拉普拉斯函数、变形密度、价层电子密度。值得一提的是，其中IRI、IRI-pi、IGMH、mIGM、amIGM、范德华势分析、价层电子密度分析都是此文的作者所提出的。

综述中还简要提及了文章作者开发的十分流行的波函数分析程序Multiwfn（主页：//www.umsyar.com/multiwfn），这是唯一可以实现此文介绍的所有方法的程序，有很多方法还是此程序独家支持的。Multiwfn对于大量相互作用可视化分析方法在计算学领域的普及和推广起到了关键性的作用，此综述里几乎所有图都是基于Multiwfn计算的数据绘制的。Multiwfn的相关信息见《Multiwfn FAQ》（//www.umsyar.com/452）和介绍文章J. Chem. Phys., 161, 082503 (2024) DOI: 10.1063/5.0216272。

此外，//www.umsyar.com/667介绍的卢天在2024年发表的Visualization Analysis of Weak Interactions in Chemical Systems（Comprehensive Computational Chemistry Vol.2中的一章）和前述的综述有极为密切的联系，此文介绍的方法专注于弱相互作用，涉及的方法更少，而由于有明显更多的篇幅，因而对它们的介绍明显更深入细致，故和上面的综述有极大的互补性，强烈建议阅读！

还十分推荐阅读以下文章及其中引用的文章获取更多相关知识以及计算操作细节：
使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用
//www.umsyar.com/621（http://bbs.keinsci.com/thread-28147-1-1.html）
使用IRI方法图形化考察化学体系中的化学键和弱相互作用
//www.umsyar.com/598（http://bbs.keinsci.com/thread-23457-1-1.html）
使用Multiwfn做Hirshfeld surface分析直观展现分子晶体和复合物中的相互作用
//www.umsyar.com/701（http://bbs.keinsci.com/thread-43178-1-1.html）
谈谈范德华势以及在Multiwfn中的计算、分析和绘制
//www.umsyar.com/551（http://bbs.keinsci.com/thread-17271-1-1.html）
在Multiwfn中单独考察pi电子结构特征
//www.umsyar.com/432（http://bbs.keinsci.com/thread-10610-1-1.html）
ELF综述和重要文献小合集
http://bbs.keinsci.com/thread-2100-1-1.html
静电势与平均局部离子化能相关资料合集
http://bbs.keinsci.com/thread-219-1-1.html
使用Multiwfn研究分子动力学中的弱相互作用
//www.umsyar.com/186
使用Multiwfn做aNCI分析图形化考察动态过程中的蛋白-配体间的相互作用
//www.umsyar.com/591（http://bbs.keinsci.com/thread-21826-1-1.html）
通过电子定域化函数(ELF)、价层电子密度分析讨论亲核进攻的方向
//www.umsyar.com/606（http://bbs.keinsci.com/thread-23982-1-1.html）
通过价层电子密度分析展现分子电子结构
http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB201709252

此外，笔者的18碳环及相关体系的研究文章中充分利用了文中介绍的许多方法，是很好的应用实例，研究汇总见//www.umsyar.com/carbon_ring.html。

波函数分析与Multiwfn程序培训班（http://www.keinsci.com/WFN）对前述的综述中几乎所有方法都有十分全面、透彻、系统的介绍，并对各种分析方法在Multiwfn中的实际操作做了细致的演示并传授大量技巧，信息量远比综述里的大得多得多，因而很推荐对这些分析方法感兴趣者关注和参加。

下面是Visualization Analysis of Covalent and Noncovalent Interactions in Real Space综述中的图，供预览。

文/Sobereva@北京科音   2025-Jun-11

0 前言

NICS和磁感生电流是非常常用的磁响应性质一类的考察芳香性的方法，在《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）里有简要介绍并里面附了大量我写过的相关博文，而在“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/WFN）的“芳香性与离域性分析”一节里有极为全面深入的介绍。很多人对这两个方法在衡量芳香性方面的认识局限在：
(1)NICS(1)ZZ约为0说明环是非芳香性，明显为负说明环是芳香性且越负芳香性越强，明显为正说明环是反芳香性且越正反芳香性越强。注意NICS的形式很多，这里仅以最常用的NICS(1)ZZ举例
(2)感生电流的方向如果满足左手定则（diatropic电流），电流越大芳香性越强；如果与左手定则方向相反（paratropic电流），电流越大反芳香性越强。电流大小一般通过对键截面上的电流密度做积分来衡量，称为bond current strength (BCS)
但实际上，仅光凭以上知识讨论芳香性或反芳香性，在一些情况下会得到不正确甚至严重错误的结论。在本文笔者介绍一些大多数人都普遍没有注意到的用NICS和感生电流讨论芳香性方面的重要问题。把这些问题都彻底搞清楚了，基于磁属性研究芳香性时就不惧审稿人comment了。

本文中的计算用的几何结构，若无专门提及，皆是B3LYP/6-31G*下优化的无虚频结构。

1 NICS、感生电流与HOMO-LUMO gap的关系

感生电流和NICS计算公式都涉及到对占据轨道和空轨道的循环，并且空轨道与占据轨道能量之差出现在分母项中，这一点看《深入理解分子轨道对磁感生电流的贡献》（//www.umsyar.com/703）里的公式1、2、3便知。总的来说，HOMO-LUMO gap越大，空轨道与占据轨道整体能量差越大，导致NICS和感生电流越接近0。因此，如果两个体系的gap相差较大，不应直接凭NICS或BCS判断二者的芳香性强弱，还应该把gap差异对结果的影响考虑进去。在《深度揭示互为等电子体的苯、无机苯和carborazine的芳香性的显著差异》（//www.umsyar.com/731）里介绍的笔者的Chem. Eur. J., 30, e202403369 (2024)一文中就注意到了这一点，文中研究的carborazine的HOMO-LUMO gap（7.30 eV）明显小于苯（11.20 eV）和无机苯（12.18 eV），因此文中提到：It is noteworthy that the aromaticity of carborazine compared to benzene and borazine is more or less quantitatively overestimated by the value of BCS and NICS。

4n电子的Huckel反芳香性体系的HOMO-LUMO gap比4n+2的Huckel芳香性体系的整体更小，Chem. Soc. Rev., 44, 6597 (2015)认为这是前者的BCS比后者更大的原因。例如反芳香性的环丁二烯和芳香性的苯在B3LYP/6-31G*级别算的HOMO-LUMO gap分别为3.692和6.800 eV，使用《使用SYSMOIC程序绘制磁感生电流图和计算键电流强度》（//www.umsyar.com/702）介绍的SYSMOIC算出来的环丁二烯的BCS为-0.7941（长C-C键）和-0.7499 a.u.（短C-C键），而苯的C-C键的BCS大小则小得多，只有0.4216 a.u.。

如《正确地认识分子的能隙(gap)、HOMO和LUMO》（//www.umsyar.com/543）所说，HF成份越低的泛函算出来的HOMO-LUMO gap越低，因此若发现HF成份低的泛函算出来的NICS和BCS的绝对值较大，大概率也是上述原因（但并非HF成份越低的泛函算出来的总是越大，因为泛函的差异还影响到其它方面）。

2 NICS的环尺寸的依赖性

原理上，感生电流比NICS判断芳香性更合理，因为NICS对环尺寸有依赖性。虽然有的文章说NICS的环尺寸依赖性弱，但那往往只是对较小的环之间对比做的结论，再加上本身不同体系的芳香性就有所不同，令那种说法更不确切。

电流对特定位置产生的磁场可以用Biot–Savart方程得到。如果是环形电流，则在环中心产生的磁场为μ0*I/(2*R)，其中μ0是真空磁导率，I是电流大小，R是环的半径，更多信息见https://en.wikipedia.org/wiki/Biot%E2%80%93Savart_law。因此，对芳香性的情况，当感生电流相同时，对应芳香性的环形路径整体半径越大、原子距离环中心整体越远，在环中心产生的磁屏蔽效果就越弱、NICS的绝对值越小。如果认为感生电流I的大小是芳香性强弱的准确体现，那么显然用NICS对比环尺寸明显不同的体系的芳香性的大小就是不甚严格的，因为对越大的环会越低估芳香性。所以，如果要基于磁属性对比明显不同大小的环的芳香性，建议用感生电流而非NICS说事。对于反芳香性，也有和上面完全类似的情况。

3 感生电流对积分截面的依赖性

对应芳香性的环上，如果只考虑对芳香性有贡献的电子的感生电流，原理上来说在垂直于键（更确切来说是垂直于电流方向）的各个截面上积分电流密度得到的BCS都应该是相同的。但实际得到的BCS没有这么理想化，环上的各个键的BCS往往不完全相同，在于：
(1)积分的数值误差、有限的积分截面尺寸
(2)积分截面可能并不严格垂直于电流方向
(3)在一些键上存在方向彼此相反的感生电流的抵消作用，此效应对于多环芳烃等并环的情况尤为明显
(4)环周围的原子的感生电流产生的干扰，以及环上原子的定域化的电子（内核电子、sigma电子、孤对电子等）产生的局部感生电流添乱，这导致计算的BCS对键上的积分截面位置的选取可能存在不可忽视的依赖性
(5)基组不够完备，导致感生电流计算不够准确。对D2h点群的环丁二烯（反芳香性）极小点结构，用SYSMOIC在B3LYP/6-31G*下计算的短和长C-C键的BCS分别为0.7499和0.7941 a.u.，而基组用完备得多的def2-QZVP计算结果分别为0.7116和0.7111 a.u.，可见后者明显更等价

若发现要考察芳香性的环上的各个键的BCS不等价性不可忽视，可以取感生电流较明确、被上述的(4)干扰小的键的BCS作为环电流值衡量芳香性，或者对环上的键的BCS取平均作为环电流值（适合碳环、环丁二烯等情况）。此外，也可以尝试计算BCS时只考虑对芳香性有贡献的电子所占据的分子轨道，例如用感生电流研究pi芳香性时就只考虑pi轨道。

4 感生电流判断芳香性和反芳香性方面相对于NICS的优势

NICS很有可能给出关于芳香性错误的结论，不如观看和积分感生电流，以及观看《通过Multiwfn绘制等化学屏蔽表面(ICSS)研究芳香性》（//www.umsyar.com/216）里介绍的ICSS图来得靠谱，这一节给出几个例子。

JCTC, 6, 1131 (2010)发现NICS用于考察氟化氢三聚体(HF)3的芳香性有问题。在wB97XD/def-TZVP下优化出的此体系的结构如下所示。凭基本常识就知道这个体系是非芳香性的，但是B3LYP/def2-TZVP算出来的NICS(0)ZZ为13.52 ppm，仿佛有显著反芳香性似的。不过NICS(1)ZZ为-0.04 ppm，这倒是正确体现出此体系并没有反芳香性。

为什么(HF)3的NICS(0)ZZ那么正？这可以从此体系所在平面上的感生电流图上理解，如下图所示（磁场垂直于屏幕向外）。可见此体系并没有形成环绕整体的感生电流，但由于存在环绕各个HF的sigma键的dia电流，导致在环中央区域呈现出等效的para电流的特征（红色箭头所示），使得环中央是去屏蔽的，因而NICS(0)ZZ为正。而由于这种效应在垂直于体系平面方向上衰减很快，因此NICS(1)ZZ几乎为0。这也体现出，对于小环体系，由于sigma键的感生电流产生的磁场离环中心太近，NICS(0)ZZ用来衡量芳香性是极差的做法。但对于18碳环及其等电子体B6N6C6这样的大环倒不是什么问题，因为局域的电子对应的局部的感生电流造成的磁场在环中心处都基本衰减没了，这也是为什么《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（//www.umsyar.com/696）介绍的文章中用了NICS(0)ZZ讨论。

实际上对于苯分子也有类似情况。下图是苯分子的pi以外的电子产生的感生电流图，可见由于局域的sigma电子产生的局部dia电流的存在，在环中心也有等效的para电流（红色箭头所示）并使得环中心产生磁去屏蔽效应。不过由于苯的pi电子造成的全局的dia电流在环中心产生的磁屏蔽效应远强于para电流的磁去屏蔽效应，因此苯的NICS(0)ZZ还是负值。

我极其推荐阅读《使用Multiwfn绘制一维NICS曲线并通过积分衡量芳香性》（//www.umsyar.com/681）中的讨论，了解sigma和pi电子对NICS_ZZ的贡献在垂直于环方向上的变化特征。

《18个氮原子组成的环状分子长什么样？一篇文章全面揭示18氮环的特征！》（//www.umsyar.com/725）里介绍的我的ChemPhysChem, 25, e202400377 (2024)文章中，对18氮环的最稳定构型计算了NICS(0)ZZ，为3.42 ppm，也是表面上体现出了轻微的反芳香性。然而通过ICSS等值面图和AICD感生电流图，都充分确认了此体系没有芳香性或反芳香性特征，详见博文里的图。这也是个很好地体现NICS(0)ZZ展现（反）芳香性不靠谱的例子。

NICS亦有可能把非芳香性体系误判成芳香性的。Theor. Chem. Acc., 134, 8 (2015)指出NICS衡量很多过渡金属团簇的芳香性也是不合理的。其中一个例子是如下所示的Ti3(CO)3。基于PCCP, 13, 4576 (2011)的SI里给出的BP86/6-311+G(d)级别下优化的结构，我用PBE0/def2-TZVP算出的NICS(0)ZZ和NICS(1)ZZ分别是-141.1和-11.7 ppm，乍看起来明显芳香性。

下图左侧和右侧分别是Ti3(CO)3的分子平面上和分子平面上方1埃处的感生电流，数值过大的部分截掉了，粉色大圆球是Ti，磁场垂直屏幕向外。由图可见，无论是在分子平面上还是在其上方1埃处，都有显著的环绕Ti原子的para电流，很大程度等效带来了在环中心及其上方的dia电流的效果（红色箭头所示），即令环中心及其上方一定区域都是磁屏蔽的，这是NICS(0)ZZ和NICS(1)ZZ都为负值的原因。然而从图上并没有看到环绕整体的连贯的感生电流，因此Ti3(CO)3就是非芳香性的。可见此例不仅NICS(0)ZZ，连NICS(1)ZZ都误判了芳香性。

5 显著且连贯的dia感生电流未必能证明有明显的芳香性

从上一节看到，感生电流判断芳香性比个别位置的NICS更为严格且直观，但即便是感生电流，也未必总能正确衡量芳香性，这里以[N6H6]2+体系为例进行说明，更多讨论见ChemistrySelect, 2, 863 (2017)。下图展示了此体系的平面构型

此构型并不是极小点，而是有很多破坏平面的虚频，并且其能量甚至比其开环后的线型结构的还要高，因此明显没有芳香稳定化作用。用Multiwfn计算出它的多中心键级几乎为0（0.00035），也进一步体现出没有pi共轭。如果用《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）的做法绘制它的LOL-pi图，也会看到各个氮都有定域的孤对电子，根本没有像苯的pi电子一样在环上充分离域。然而，下图所示的[N6H6]2+平面构型的分子上方1埃处的感生电流（磁场垂直于屏幕向外）则是完全顺时针的dia电流，和具有显著芳香性的苯的情况看起来完全一样！

为什么[N6H6]2+的平面构型不是极小点？这是因为如下图的轨道占据情况所示，尽管此体系表面上有10个pi电子，看似满足Huckel的4n+2芳香性规则因而芳香性会使此构型被稳定化，但e2u对称性的轨道的反pi特征显著（对于苯来说是没有占据的），因此pi共轭作用很弱，再加上如此众多的孤对电子间的强烈互斥作用，显然令其极不稳定。

为什么[N6H6]2+的平面构型像普通芳香性体系一样有显著的贯穿整个环的dia电流？这需要利用《深入理解分子轨道对磁感生电流的贡献》（//www.umsyar.com/703）里介绍的知识，没看过者强烈建议认真阅读。此文说了，感生电流可视为由各种占据轨道向各种空轨道跃迁所贡献，上图蓝色箭头标注的跃迁是此结构下对感生电流贡献最大的跃迁（因为这种跃迁能最小），其中非占据的b2g轨道比占据的e2u轨道轨道多一个节面，这种跃迁是平动允许的，因此对dia电流有贡献。这和造成苯分子的dia电流的机制是完全相同的，只不过苯分子的dia电流主要是e1g到e2u跃迁所贡献的。

由此例可见，光靠是否存在全局的dia感生电流判断芳香性并不总是可靠！因为这种电流是否存在，和体系是否真的有芳香性特征并没有内在的必然联系，而与它关系真正最直接、最紧密的是轨道的对称性特征！

可能大多数人觉得，有芳香性的体系有满足左手定则的dia电流的原因在于：某个环有芳香性意味着环上有明显离域的电子，因此在外磁场的作用下会像经典电磁学描述的导体一样在环上产生dia电流。但实际上真正本质并非如此，不能忽视关键性的量子力学效应！一般的芳香性体系有dia电流实际上是一个表象，在于它们的轨道对称特征和跃迁正好能满足产生显著dia电流的条件，而能产生显著dia电流的并非只有芳香性体系。因此，要这样认为：有显著的贯穿整个环的dia电流是这个环具有芳香性的必要非充分条件。

由于磁屏蔽特征是感生电流的进一步衍生属性，[N6H6]2+的平面构型既然都和苯一样具有相同感生电流特征了，其自然NICS(0)ZZ和NICS(1)ZZ也都为负，分别为-178.1和-24.6 ppm。

此例彰显出判断芳香性最好不要光基于磁属性，即便这种判断芳香性的做法极为流行。同时结合多中心键级等基于电子结构的方法，以及结合芳香稳定化能或离域化能这种能量层面的方法，结论会更可靠、更有说服力。本节的[N6H6]2+的例子也体现出，形式上满足4n+2规则的体系未必真的有芳香性，若电子都没法在环上充分离域、不构成共轭体系，满足4n+2也仍是非芳香性的。

6 显著且连贯的para感生电流未必能证明有明显的反芳香性

J. Org. Chem., 88, 14831 (2023)通过一系列分子指出，即便环上存在显著且连贯的para电流，也未必能可靠地证明这个环就有明显的反芳香性。这一节就用其中最简单的体系，苯+1电荷的阳离子自由基（具有5个pi电子）和-1电荷的阴离子自由基（具有7个pi电子）来说明这一点，这种体系的（反）芳香性没法基于经典的Huckel规则判断。

如下所示，JOC文中给出了苯阳离子、中性苯、苯阴离子的分子平面上方1 Bohr处的感生电流，磁场垂直于屏幕向外，计算用的是对各自状态优化的几何结构，电场强度越大箭头越蓝、越小越红。可见，苯阳离子和苯阴离子在碳环上方的电流都是逆时针的para电流，与中性苯的情况相反，乍看起来苯的阴、阳离子都是显著反芳香性的。

然而，JOC这篇文章计算了苯的阴、阳离子的芳香稳定化能（用到的反应见文章的SI），发现这两个体系是有芳香稳定化作用的，因此从能量角度来说是有芳香性的。我还通过Multiwfn在文中给的苯离子的极小点结构下在B3LYP/6-31G*级别计算了六中心键级，苯的阳离子和阴离子的结果分别为0.044和0.045，只有中性苯的0.086的约一半。具体来说，苯阳离子和阴离子的六中心键级分别几乎来自于alpha和beta电子的贡献（且和中性苯的相应自旋电子的贡献差不多），而另一个自旋几乎没贡献。也就是说，相对于苯，离子状态下电子数改变1的那个自旋的电子对芳香性没明显贡献，而电子数未变的自旋的电子对芳香性的贡献基本保持不变。

由于从能量和电子结构角度都确认了苯的阴、阳离子是芳香性的，那么为什么它们的感生电流是para的，显得像普通反芳香性体系一样？JOC文中给了下面的图对感生电流进行了解释。由于Jahn-Teller效应，阳离子（5π）和阴离子（7π）体系的对应于苯的二重简并的pi轨道发生了轻微分裂。图中蓝色大箭头标注了对感生电流产生主导的跃迁。对苯离子体系的情况，由于能级分裂造成的能量差较小，明显小于1→2的跃迁，因此对阳离子和阴离子分别是1→1和2→2跃迁对感生电流产生了主导性的贡献。由于这种跃迁涉及的两个轨道的节面数相同，是转动允许的，因此如《深入理解分子轨道对磁感生电流的贡献》所述，贡献的是para电流。

由本节和上一节的例子可见，虽然感生电流对大多数体系判断芳香性是合理的，而且比起只看个别位置的NICS_ZZ可靠不少，但在极个别体系上还是能暴露出其局限性。因此若有可能，除了基于磁属性分析外，还是建议再结合一两种靠谱的手段，特别是可靠又容易计算的多中心键级，来衡量芳香性，减少误判芳香性、被审稿人怼的可能（本来芳香性这方面的“说法”和“观点”就相当多）。

但也没必要盲目用过多的方法衡量芳香性，否则一方面太冗余，另一方面本身有的芳香性指标就不那么靠谱、普适，纳入讨论反倒无益于揭示芳香性的真相。有文章认为键的均衡化、出现dia和para的感生电流，属于芳香性的二阶（次级）特征，而能量的稳定化和电子的离域，才算是芳香性的一阶（即最关乎本质的）特征，这个观点我是认可的。

文/Sobereva@北京科音  2025-Apr-5

1 前言

Dalton这个 程序虽然用起来复杂，但有很多独特的亮点，例如可以在SOPPA(CCSD)和CASSCF级别下计算各种磁属性（笔者最近在Chem. Eur. J., e202404138 (2025)文中就用Dalton在CASSCF下算了NICS），可以用响应方式较准确地计算磷光寿命，可以算双光子吸收（北京科音高级 培训班http://www.keinsci.com/KAQC里专门有一节讲了)，可以用高级别耦合簇方法算激发态（见《使用Dalton通过CC3方法极高精度计算激发态》//www.umsyar.com/463），等等。可惜Dalton的SCF收敛方面做得一般，比不上Gaussian，而且单靠Dalton无法产生一些重要的轨道，比如UHF计算产生的UNO轨道经常被用于CASSCF的初猜，但Dalton甚至连UHF都做不了，因为Dalton里开壳层只支持RO形式。

为解决以上问题，从2025-Apr-5更新的Multiwfn开始，Multiwfn产生Dalton输入文件的功能支持了将Multiwfn中当前的轨道波函数写入到Dalton的输入文件.dal里作为初猜波函数。因此借助Multiwfn，Dalton用户可以用Gaussian、ORCA、GAMESS-US等诸多主流的 程序产生的波函数当初猜，从而解决Dalton遇到SCF不收敛的问题。而且，Dalton用户还可以直接用Multiwfn产生的UNO、定域化分子轨道（LMO，见//www.umsyar.com/380）等当CASSCF初猜。

下面就通过例子简单演示一下具体实现方法，本文中涉及的文件都可以在//www.umsyar.com/attach/740/file.zip下载。Multiwfn可以在//www.umsyar.com/multiwfn免费下载，相关常识看《Multiwfn FAQ》（//www.umsyar.com/452）。如果读者还不知道Dalton怎么安装、运行、产生输入文件的话，看《 程序Dalton的编译方法和运行方式简介》（//www.umsyar.com/463）。本文用的Multiwfn是2025-Apr-5更新的版本，Gaussian是G16 B.01，Dalton是2022版。

如果本文介绍的Multiwfn的功能给你的研究带来了帮助，请在论文中引用Multiwfn刚启动时提示的Multiwfn程序的原文，这是对Multiwfn开发和维护最好的支持！

2 把Gaussian收敛的HF波函数当Dalton做CCSD(T)的初猜

这一节以CH3NH2为例演示把Gaussian收敛的HF/cc-pVTZ波函数当Dalton做CCSD(T)/cc-pVTZ的初猜。首先用Gaussian做一个普通的单点计算，输入文件是本文文件包里的CH3NH2.gjf，内容如下所示。算完后把得到的CH3NH2.chk用formchk转成CH3NH2.fch。

%chk=D:\CH3NH2.chk
# HF/cc-pVTZ int=NoBasisTransform
[空行]
test
[空行]
0 1
 C                  0.05159500    0.70381800    0.00000000
 H                  0.59439800    1.06209400    0.88166000
 H                  0.59439800    1.06209400   -0.88166000
 H                 -0.94293100    1.18451300    0.00000000
 N                  0.05159500   -0.76078400    0.00000000
 H                 -0.45830300   -1.10306000    0.81258800
 H                 -0.45830300   -1.10306000   -0.81258800
[空行]
[空行]

之后启动Multiwfn，载入CH3NH2.fch，依次输入
100  //其它功能（Part 1）
2  //产生输入文件
19  //产生Dalton输入文件
CCSDpT.dal  //要产生的输入文件名
[回车]  //要产生的.mol文件名用默认的CH3NH2.mol
y  //由于当前内存里有基函数信息，Multiwfn问你是否把轨道展开系数写入到.dal文件里
1  //使用Dalton标准格式（4F18.14）方式写入。即四个值换一行，每个值占18列、小数位数14个

此时当前目录下就有了CCSDpT.dal和CH3NH2.mol。把CH3NH2.mol里的基组名都替换为cc-pVTZ_emsl。目前CCSDpT.dal里的计算设置对应B3LYP算单点，自行把关键词改为CCSD(T)计算用的，此时CCSDpT.dal的内容如下
**DALTON INPUT
.RUN WAVE FUNCTIONS
**WAVE FUNCTIONS
.CC
*CC INPUT
.CC(T)
*ORBITAL
.MOSTART     //代表从**MOLORB部分读取轨道展开系数作为初猜
 FORM18
**MOLORB (Punched by Multiwfn)
 -0.00001750284370  0.00060444251600  0.00013849165900 -0.00001024308040
  0.00002442752580 -0.00001726827680 -0.00014237439200  0.00004058047020
...略
**END OF INPUT

用Dalton进行计算，可以看到SCF一轮就收敛了，达到了我们的目的：
@ *** DIIS converged in   1 iterations !
@     Converged SCF energy, gradient:    -95.252488149534    1.44D-07
    - total time used in SIRFCK :              0.00 seconds
之后CCSD(T)计算也顺利跑完了。

一些注意事项和相关信息：
• 对于广义或部分广义收缩基组，都应当用int=NoBasisTransform关键词，拿不准就始终带上。详见《将Gaussian等程序收敛的波函数作为ORCA的初猜波函数的方法》（//www.umsyar.com/517）里对此关键词的解释。
• 把Gaussian的cc-pVDZ、cc-pVTZ的波函数当Dalton的初猜时，Dalton要用cc-pVDZ_emsl、cc-pVTZ_emsl基组，而不要用不带_emsl后缀的，否则还是需要迭代很多次才能收敛（甚至可能不收敛），但带不带_emsl的版本收敛后的结果是相同的。而如果是def、def2系列基组，没什么特别要注意的，比如Gaussian里def-TZVP和def2-TZVP做的计算，在Dalton里基组名分别写成Turbomole-TZVP、def2_tzvp即可，注意对大小写有要求，必须对应Dalton目录下的basis子目录里的基组文件名。
• Dalton默认是基于球谐型基函数做的计算，因此如果在Gaussian里使用的是比如6-31G*等D壳层默认是笛卡尔型基函数的基组，应当写上5d关键词强行要求用球谐型，详见《谈谈5d、6d型d壳层基函数与它们在Gaussian中的标识》（//www.umsyar.com/51）。虽然Dalton也可以用笛卡尔型基函数，但Multiwfn往Dalton输入文件里写入初猜轨道系数的功能不支持此情况。
• 只要一个 程序能输出Multiwfn可以识别的带有基函数信息的波函数文件，并且用的是球谐型基函数，则波函数都可以通过以上方式作为Dalton的初猜。支持的文件详见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）。例如ORCA做计算产生的molden文件也可以载入Multiwfn后像上面那样产生带初猜波函数的Dalton输入文件。

3 把Gaussian收敛的DFT波函数当Dalton做DFT的初猜

此例以顺铂体系为例，把Gaussian收敛的PBE0波函数当Dalton做PBE0计算的初猜，对Pt用Stuttgart小核赝势，对配体用6-311G*。Gaussian输入文件如下（本文文件包里的Ptcoord.gjf）：

%chk=D:\Ptcoord.chk
#P PBE1PBE/genecp int=NoBasisTransform
[空行]
b3lyp/def2TZVP opted    //众所周知这是纯给自己看的标题行，别问为什么写成这个
[空行]
0 1
[坐标略...]
[空行]
N H Cl
6-311G*
****
Pt
SDD
****
[空行]
Pt
SDD
[空行]
[空行]

按照和上一节例子相同的过程，把chk转成fch，再用Multiwfn产生.dal和.mol文件。.dal里的泛函名改成PBE0，然后设置.mol里的基组和赝势成为下面这样，当前对Pt用的赝势和赝势基组和Gaussian里的情况一样（若不确定的话，可以Gaussian里用gfinput关键词输出定义，并和Dalton的基组、赝势目录下的相应文件对照）。

ATOMBASIS
test molecule
Generated by Multiwfn
Atomtypes=4 Angstrom Nosymmetry charge=0
Charge=1.0 Atoms=6 Basis=6-311G*
H1   -0.82596600    1.64390200    2.14978600
H2    0.00000000    2.40772900    0.93575300
H3    0.82596600    1.64390200    2.14978600
H4   -0.82596600   -1.64390200    2.14978600
H5   -0.00000000   -2.40772900    0.93575300
H6    0.82596600   -1.64390200    2.14978600
Charge=7.0 Atoms=2 Basis=6-311G*
N1    0.00000000    1.59755500    1.56108400
N2   -0.00000000   -1.59755500    1.56108400
Charge=17.0 Atoms=2 Basis=6-311G*
Cl1    0.00000000    1.70827400   -1.36819100
Cl2   -0.00000000   -1.70827400   -1.36819100
Charge=78.0 Atoms=1 Basis=stuttgart_rsc_1997_ecp ECP=stuttgart_rsc_1997_ecp
Pt1   -0.00000000   -0.00000000    0.18195700

用Dalton计算后，经过4轮SCF就轻易收敛了，结果为-1152.597056，和Gaussian给出的-1152.596983非常接近。而如果此例Gaussian和Dalton都用HF的话，会和上例一样，结果精确相同且SCF一轮可收敛。

4 把UNO轨道当Dalton做CASSCF的初猜

之前在《CASSCF计算双自由基以及双自由基特征的计算》（//www.umsyar.com/264）里示例过用Gaussian基于UNO初猜轨道做CASSCF(2,2)计算的方法，不了解这种计算和UNO概念的话建议先看看。这一节演示使用Gaussian产生UHF波函数，读入Multiwfn并产生UNO轨道，然后写入.dal文件用于Dalton做CASSCF(2,2)计算的初猜轨道的流程。示例体系还是264号博文里的C4H8双自由基，基组还是此文里的6-31G*。

首先用以下Gaussian输入文件（本文文件包里的C4H8.gjf）对C4H8做UHF计算得到对称破缺波函数。不了解为什么用guess=mix的话看《谈谈片段组合波函数与自旋极化单重态》（//www.umsyar.com/82）。

%chk=D:\C4H8.chk
# UHF/6-31G* guess=mix 5d
[空行]
ub3lyp/6-31g(d) opted
[空行]
0 1
C     -0.74742092    1.77656753    0.00000000
H     -0.62438907    2.32965189    0.92333358
H     -0.62438907    2.32965189   -0.92333358
C     -0.74742092    0.30933780    0.00000000
H     -1.24978876   -0.09122321    0.88294562
H     -1.24978876   -0.09122321   -0.88294562
C      0.74742092   -0.30933780    0.00000000
H      1.24978876    0.09122321   -0.88294562
H      1.24978876    0.09122321    0.88294562
C      0.74742092   -1.77656753    0.00000000
H      0.62438907   -2.32965189   -0.92333358
H      0.62438907   -2.32965189    0.92333358
[空行]
[空行]

计算完毕后把C4H8.chk转成C4H8.fch，载入Multiwfn，然后依次输入
200  //其它功能（Part 2）
16  //基于fch文件里的密度矩阵产生自然轨道。这个功能在《在Multiwfn中基于fch产生自然轨道的方法与激发态波函数、自旋自然轨道分析实例》（//www.umsyar.com/403）里有专门的介绍
SCF  //读取fch里的SCF类型的密度矩阵
1  //产生空间自然轨道
此时从屏幕上看到如下信息，这便是UNO轨道的占据数。其中最高占据自然轨道（HONO）和最低非占据自然轨道（LUNO）的占据数分别是1.123225和0.876775，都偏离整数显著，很满足期望。

 Occupation numbers:
    2.000000    2.000000    1.999999    1.999999    1.999904    1.999898
    1.999881    1.999835    1.999309    1.999120    1.999004    1.998977
    1.996829    1.996536    1.996053    1.123225    0.876775    0.003947
    0.003464    0.003171    0.001023    0.000996    0.000880    0.000691
...略

接着输入
n  //不把自然轨道导出成new.mwfn文件并重新载入，因为当前不需要如此（如果你想看一下当前的UNO轨道图像的话，应选y，之后用Multiwfn主功能0观看轨道）
0  //返回主菜单
100  //其它功能（Part 1）
2  //产生输入文件
19  //产生Dalton输入文件
CASSCF.dal  //要产生的输入文件名
[回车]  //要产生的.mol文件名用默认的C4H8.mol
y  //把当前内存里的轨道（即UNO）的展开系数写入到.dal文件里
1  //使用Dalton标准格式（4F18.14）方式写入

当前C4H8.mol里的基组默认就是6-31G*，不用改。修改CASSCF.dal，把以下两行
.DFT
 B3LYPg
替换为CASSCF(2,2)计算的设置：
.MCSCF
*CONFIGURATION INPUT
.CAS SPACE
2
.ELECTRONS
2
.INACTIVE
15

Dalton计算完毕后，活性空间内自然轨道占据数为1.2679和0.7320，能量为-156.024938 Hartree，和Gaussian在同级别做CASSCF精确一致。你若想看Dalton做CASSCF产生的自然轨道，把计算产生的.tar.gz文件里的molden.inp解压出来载入Multiwfn，并效仿《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）的方式观看即可。

5 关于线性依赖基函数的问题

较大基组（尤其是带弥散函数的），计算原子分布较致密的体系很容易出现线性依赖问题，Gaussian和Dalton此时都会自动去除线性依赖基函数，导致最终求解出来的轨道数会少于基函数数目。这种情况下，Gaussian输出信息中的NBsUse（实际用的基函数数目，等同于实际求解出来的轨道数）小于NBasis，Dalton输出文件中的Total number of orbitals小于Number of basis functions。这给上文的传轨道带来了麻烦。有两个解决方法：
(1)不去除线性依赖基函数：Gaussian计算时带着IOp(3/32=2)关键词、Dalton计算时在**WAVE FUNCTIONS里面加上以下内容，从而让两个程序都不去除线性依赖基函数
*ORBITAL
.AO DELETE
1E-100
这种做法的缺点是可能导致数值不稳定性。
(2)自动去除线性依赖基函数，但让Dalton少读取轨道：比如Gaussian计算显示NBsUse=50、NBasis=54，因此自动去除了4个线性依赖基函数。因此当Multiwfn载入fch文件后，最后4个轨道的展开系数都会为0，导出到.dal文件里也是如此。此时Dalton若照常读取全部轨道，计算一开始会报错*** ORTHO-FATAL ERROR ***。解决方法是在*ORBITAL里写上
.DELETE
4
这样Dalton就不会读取最后4个轨道了，计算可以正常进行。

6 .dal里的轨道展开系数存在星号而无法读取的问题

有时候Multiwfn写入到.dal里的轨道展开系数会存在一堆星号，导致Dalton读取轨道时提示input conversion error错误。这不是Multiwfn的问题，而是Dalton自身设计不周。Dalton定义的**MOLORB字段中记录展开系数的格式在前面说了，每个数值都是F18.14格式，因此当有的轨道展开系数特别大，超过了其整数位记录上限，就会被记录成一堆星号。这种问题在基组较弥散而又不让自动去除线性依赖基函数的情况下出现概率大。这是为什么Multiwfn导出Dalton输入文件的时候特意让你选择记录的格式，前面的例子都是选1用4F18.14格式，还可以选2用4E20.12格式，这是使用科学计数法输出，因此就不会出现数值过大记录成星号的问题了。相应地，必须自己修改Dalton的读取轨道展开系数的源代码。对于Dalton 2022版，需要修改Dalton源代码目录下的DALTON/sirius/sirgp.F文件，把PFMT = '(4F18.14)'替换为PFMT = '(4E20.12)'，之后重新编译即可。

7 其它

Multiwfn就相当一个大熔炉，能读取诸多 程序产生的波函数，又能导出波函数作为诸多 程序的初猜，这在《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）里专门说了。对于传轨道到ORCA的情况我之前还专门写过《将Gaussian等程序收敛的波函数作为ORCA的初猜波函数的方法》（//www.umsyar.com/517）。

MOKIT也是很好的可以实现不同程序间传轨道的工具，fch到.dal和.mol的转换可以用里面的fch2dal子程序。相对于MOKIT，本文介绍的用Multiwfn传轨道有几个额外的好处：
(1)Multiwfn不需要安装MOKIT所需的Python环境。尤其是对于Windows用户，Multiwfn解压即用明显更为方便
(2)Multiwfn特意考虑了.dal里记录轨道展开系数的格式的问题并给了解决办法，上一节说了
(3)Multiwfn自身可以产生诸多类型的轨道（定域化分子轨道、空间/自旋自然轨道、NAdO、AdNDP、NTO、双正交化轨道等等）并输出作为Dalton的初猜
而fch2dal也有额外的优点，它会把fch里的基组、赝势的具体定义直接写到输出的.mol里，因此如果Gaussian计算时用的基组、赝势在Dalton里没有自带也能直接计算而不需要自定义。并且fch2dal支持用笛卡尔型基函数的情况。

文/Sobereva@北京科音   2025-Mar-31

0 前言

使用计算化学程序过程中普遍离不开可视化问题，对如今已经非常流行、十分强大、巨快、开源免费的第一性原理程序CP2K自然也不例外。最近有人在思想家公社QQ群里问“cp2k可视化一般搭配哪个程序呢？”，针对这个问题，笔者简要罗列一下我推荐的在CP2K使用过程中的可视化程序，主要面向初学者。本文只是粗略概述，这些程序在笔者讲授的北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/KFP）中都会反复利用和详细演示具体操作。除本文说的外，还有大量其它可视化程序，比如Avogadro、gabedit、Chemcraft等，我认为对于CP2K用户来说，用好本文提到的这些就已经完全足够了。文中提到的Multiwfn可以在官网//www.umsyar.com/multiwfn下载，相关信息见《Multiwfn FAQ》（//www.umsyar.com/452）。

1 建模阶段的可视化

能够构建用于做第一性原理计算的周期性体系的程序很多，其中我十分推荐GaussView，这是绝大多数内行做 计算的人都会用的可视化程序，其图形界面设计得很好，所见即所得，大多数情况用这一个程序就可以很好完成CP2K计算前的建模。GaussView可以直接读取记录周期性体系原子坐标和晶胞信息的最常见的cif文件。在GaussView里构建好结构后，可以保存成Gaussian程序的输入文件（gjf），里面有原子坐标信息和晶胞的平移矢量信息（Tv开头的行），然后用Multiwfn载入这样的gjf文件后，就可以按照《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）文中说的创建CP2K输入文件了，整个过程相当方便。

许多人误以为GaussView只是结合Gaussian使用的，或者至多是用于 程序算孤立体系建模的，这是大错特错！由于Gaussian本来就有计算周期性的功能（但由于功能性和效率的原因，极少有人用Gaussian算周期性体系），所以Gaussian的御用图形界面GaussView也支持周期性体系，可以直接对周期性体系添加/删除/替换原子、调整键长/键角/二面角、平移/旋转片段、把其它体系粘贴到当前体系里、测量几何参数，等等，对周期性体系和对分子体系操作一样非常灵活和顺手。并且GaussView专门有个PBC editor界面，在里面可以定义晶胞参数、重新定义晶胞、修改周期性、晶轴变换、扩胞、切晶面、判断空间群、居中体系、删除晶胞外的原子、把晶胞外的原子卷入晶胞，等等，常用的操作几乎都提供了。GaussView对于周期性体系建模相当好用，然而其这方面的价值却被很多人低估，甚至有人明明不了解还胡乱贬损。

点击下图红色箭头处的按钮就可以进入前述的PBC editor界面。红、绿、蓝边框分别是晶胞的a、b、c轴。

还值得一提的是，GaussView的菜单栏的Tools中的Atom selection界面很有用。用滑框选择工具选择一批原子成为黄色后，Atom selection界面里就可以看到选中的原子序号，格式和Multiwfn程序里要求的完全一致。例如在Multiwfn里创建几何优化任务的输入文件时，若你选择要冻结某些原子，就可以把序号从这里复制出来并直接粘贴到Multiwfn的窗口中，这样实现诸如slab边缘原子的冻结设置巨方便。

Multiwfn还可以载入CP2K的输入文件、Quantum ESPRESSO的输入文件、VASP的POSCAR等格式，主功能100的子功能2里选择保存为gjf或者cif文件后，也可以用GaussView载入并观看和编辑。顺带一提，在建模方面Multiwfn也提供了很多实用功能，很多还是GaussView没有的，见《Multiwfn中非常实用的几何操作和坐标变换功能介绍》（//www.umsyar.com/610）。

2 几何优化结果的可视化

如果要看CP2K做几何优化得到的最终结构，最简单的方法是用Multiwfn载入此任务产生的restart文件，之后进入主功能0直接就看到了，例如下面的聚噻吩体系。界面上方的菜单以及图形界面右侧有一大堆选项可以调整效果，诸如视图的旋转和缩放、键的粗细、成键判断阈值、原子序号是否显示/字号/颜色、原子球大小、高亮特定原子。菜单栏的Tools列表里还有一些辅助工具，比如测量几何参数、导出所有内坐标、输出笛卡尔或分数坐标，等等。

很多时候我们还关心优化的过程，尤其是当体系结构复杂时，不观看优化轨迹的话，往往都不好判断优化过程中哪里发生了何种变化，心里没数。观看优化轨迹最好的程序是VMD，在http://www.ks.uiuc.edu/Research/vmd/可免费下载，强烈建议用VMD 1.9.3（撰此文时1.9.4还没发布正式版，其测试版的bug巨多，强烈不建议用）。CP2K的几何优化任务默认会输出.xyz文件，这是多帧的轨迹文件，里面记录了优化过程的每一帧的坐标，不熟悉xyz格式的话看《谈谈记录化学体系结构的xyz文件》（//www.umsyar.com/477）。把xyz文件载入VMD就可以观看优化轨迹了。即便任务没跑完，也可以载入此文件看看最新一步的结构。

xyz文件里不记录晶胞信息，因此在VMD里没法显示出晶胞边框，给观看周期性体系造成不便。对于不变胞的优化，可以用Multiwfn载入restart文件，此时屏幕上不仅显示晶胞信息，还会显示在VMD里定义晶胞的命令，如下图所示。

 把这命令复制到VMD的命令行窗口运行，之后再输入pbc box命令，晶胞在VMD里就显示出来了，如下所示。

如果做的是变胞优化任务，晶胞在优化过程中不是固定不变的，此时建议让CP2K输出pdb格式的优化轨迹，里面记录了每一帧的晶胞信息。此文件载入VMD并显示盒子后，会看到随着轨迹的播放，盒子也实时变化。

上述说明不仅适用于优化极小点，对于dimer方法优化过渡态也是一样的。

VESTA（http://jp-minerals.org/vesta/en/）也是一个很好且免费的主要面向周期性体系的可视化工具，默认情况下的显示效果往往比VMD甚至还更好些。若想用VESTA看CP2K优化的结果，可以把restart文件载入Multiwfn，主功能100的子功能2里选择导出cif或VASP的输入文件（POSCAR），然后就能载入VESTA了。

3 反应路径的可视化

CP2K做NEB类任务需要对给定的初始结构之间插点。CP2K虽然能够自动插点（点=image），但没法在计算前进行预览。利用《sobNEB：产生CP2K的NEB的插点的方便的工具》（//www.umsyar.com/660）里介绍的笔者的sobNEB程序插点的话，可以直接产生traj.xyz轨迹文件，放到VMD里通过播放动画或者多帧叠加显示，就可以直观判断初始的NEB的image是否合理了。

NEB类任务计算过程中是好多副本一起算的，每个副本都会在计算过程中输出它所负责的那个image优化过程的xyz轨迹文件。要想观看最终收敛的NEB轨迹，需要自己把这些副本输出的轨迹的最后一帧合并在一起作为新的xyz轨迹，放到VMD里就可以观看最终的反应路径动画了。手动做比较麻烦，建议用脚本实现。北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/KFP）里专门给了个shell脚本自动做这件事情，运行后就会在当前目录下产生[项目名]_traj.xyz文件，里面记录了NEB最新一步的轨迹（即便任务还没跑完，也可以看最新的NEB轨迹是什么样），同时还会产生[项目名]_ene.txt，里面记录了最新一步的各个点的能量。下面是培训里的一个实例，对载入的[项目名]_traj.xyz文件做多帧叠加显示，直观显示了Na的迁移轨迹

4 振动分析的可视化

CP2K做完振动分析后往往需要观看振动矢量了解振动的特征。推荐做法是用《使CP2K计算的振动模式可以被GaussView观看的程序：MfakeG》（//www.umsyar.com/656）里介绍的笔者开发的MfakeG工具，把记录振动模式信息的.mol后缀的molden文件转化成仿Gaussian输出文件，之后载入GaussView就可以用Results - Vibrations选项观看振动模式了，例如下图的草酸晶体的例子

还有一种做法是使用免费的可视化程序Jmol（https://jmol.sourceforge.net），也可以载入CP2K产生的molden文件观看振动模式。不过我还是觉得GaussView用起来舒服。

如果要基于CP2K的振动分析观看振动光谱，推荐用Multiwfn。Multiwfn具有非常强大的绘制各种光谱的功能，见《使用Multiwfn绘制红外、拉曼、UV-Vis、ECD、VCD和ROA光谱图》（//www.umsyar.com/224）和《使用Multiwfn绘制NMR谱》（//www.umsyar.com/565）。Multiwfn可以载入CP2K做振动分析产生的输出文件按博文所述绘制红外光谱，如果要求计算了拉曼还可以绘制拉曼光谱。此外，Multiwfn还可以载入常规TDDFT和XAS-TDDFT任务的输出文件分别绘制UV-Vis光谱和X光吸收光谱，还可以考虑旋轨耦合效应，例如下图是CP2K+Multiwfn绘制的NaAlO2晶体的Al的K-edge XAS。Multiwfn还可以载入CP2K的NMR任务的输出文件绘制NMR谱。这些在北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/KFP）里有非常系统的讲解和丰富的例子。

5 动力学轨迹可视化

CP2K的第一性原理动力学（叫从头算动力学也行）是CP2K的杰出强项之一。和GROMACS、AMBER、NAMD等经典力场动力学程序的情况一样，VMD也是观看其动力学轨迹的不二的选择。例如培训里有个模拟质子轰击石墨烯层的例子，VMD载入轨迹后并多帧叠加显示、根据帧号着色，直观展现了模拟过程中质子的运动轨迹：

还可以自己写VMD脚本从CP2K产生的记录原子速度的xyz文件中把速度信息读入VMD，并根据速度着色，展现出质子的动能是如何传递到石墨烯板上并扩散开来的。下图黄色是打入的质子，越蓝的石墨烯原子的速度越大。

VMD不仅可以观看动力学过程中结构的变化，还可以绘制等值面图观看电子结构的变化，例如下图这个例子，直观展现了水合电子是怎么在模拟过程中出现的。

通过写VMD tcl脚本，还可以实现很复杂的分析，如下面幻灯片里讲的高温下正癸烷的裂解产物分析。PS：做动力学的，不会写分析脚本的话，稍微做深一点就会寸步难行。

北京科音分子动力学与GROMACS培训班（http://www.keinsci.com/KGMX）专门用约100页幻灯片特别完整、详细讲VMD的使用，并且还额外用近100页幻灯片深入讲VMD脚本的编写，如果想系统学习VMD的话这是极佳的途径。

笔者偶尔看到有人用OVITO程序可视化CP2K产生的轨迹，我没用过那个程序，也完全不理解为什么有人不用VMD而用OVITO，明明VMD已经可以完美地满足一切需求。笔者在答疑时看到有一些CP2K用户还被OVITO坑了：OVITO根据边缘原子位置自动确定了盒子边框，有人以为那就是实际的晶胞边界，导致对结果产生了严重误解。

6 电子结构的可视化和分析

Multiwfn可以基于CP2K的输出文件绘制效果非常理想的能带结构图，如CrO2体系：

《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）中介绍了怎么用CP2K产生molden文件。Multiwfn将之载入后，就可以在《使用Multiwfn绘制态密度(DOS)图考察电子结构》（//www.umsyar.com/482）讲的基础上举一反三绘制效果很好的DOS、PDOS、OPDOS、LDOS图，下图是WO3体系：

Multiwfn还能基于CP2K的molden文件观看轨道，做法可参考《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）。还支持对特定k点看轨道：

Multiwfn还能对周期性体系做超级丰富的波函数分析，其中很多都是以图形方式展现的，比如可以基于《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）和《使用Multiwfn考察周期性体系的芳香性》（//www.umsyar.com/722）讲的做法计算LOL-pi格点数据，之后可以在Multiwfn里直接观看等值面；也可以导出成cube文件后，载入VMD或VESTA程序绘制，如下所示，极为生动直观地展现了一个COF体系的pi电子共轭路径

再比如下图是Multiwfn载入CP2K对硅表面计算产生的molden文件，做轨道定域化后直接显示出轨道等值面图，直观展现了体系不同位置的电子结构特征。相关信息见《Multiwfn的轨道定域化功能的使用以及与NBO、AdNDP分析的对比》（//www.umsyar.com/380）。

为避免此章太长，Multiwfn可以针对周期性体系做的巨量的可视化分析这里就不一一提及，很多分析我都写过博文，感兴趣的读者请阅读：
使用Multiwfn结合CP2K的波函数模拟周期性体系的隧道扫描显微镜（STM）图像
//www.umsyar.com/671（http://bbs.keinsci.com/thread-37740-1-1.html）
使用Multiwfn对周期性体系做键级分析和NAdO分析考察成键特征
//www.umsyar.com/719（http://bbs.keinsci.com/thread-47176-1-1.html）
使用Multiwfn结合CP2K做周期性体系的atom-in-molecules (AIM)拓扑分析
//www.umsyar.com/717（http://bbs.keinsci.com/thread-46927-1-1.html）
使用Multiwfn结合CP2K对周期性体系做电荷分解分析（CDA）
//www.umsyar.com/716（http://bbs.keinsci.com/thread-46878-1-1.html）
使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用
//www.umsyar.com/621（http://bbs.keinsci.com/thread-28147-1-1.html）
使用IRI方法图形化考察化学体系中的化学键和弱相互作用
//www.umsyar.com/598（http://bbs.keinsci.com/thread-23457-1-1.html）
使用Multiwfn结合CP2K通过NCI和IGM方法图形化考察固体和表面的弱相互作用
//www.umsyar.com/588（http://bbs.keinsci.com/thread-21742-1-1.html）
使用Multiwfn绘制分子和固体表面的距离投影图
//www.umsyar.com/589（http://bbs.keinsci.com/thread-21754-1-1.html）
使用Multiwfn图形化展现原子对色散能的贡献以及色散密度
//www.umsyar.com/705（http://bbs.keinsci.com/thread-44723-1-1.html）
使用Multiwfn做Hirshfeld surface分析直观展现分子晶体和复合物中的相互作用
//www.umsyar.com/701（http://bbs.keinsci.com/thread-43178-1-1.html）
使用CP2K结合Multiwfn对周期性体系模拟UV-Vis光谱和考察电子激发态
//www.umsyar.com/634（http://bbs.keinsci.com/thread-28006-1-1.html）
使用Multiwfn计算分子和晶体中孔洞的直径
//www.umsyar.com/643（http://bbs.keinsci.com/thread-30696-1-1.html）
使用Multiwfn计算晶体结构中自由区域的体积、图形化展现自由区域
//www.umsyar.com/617（http://bbs.keinsci.com/thread-25241-1-1.html）
使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线
//www.umsyar.com/638（http://bbs.keinsci.com/thread-28225-1-1.html）

此外，用VMD载入CP2K产生的Hartree势（静电势的负值），还可以实现对晶体表面静电势的可视化，一个COF体系的例子如下所示。静电势在研究静电主导的相互作用方面有重要意义，参见《静电势与平均局部离子化能相关资料合集》（http://bbs.keinsci.com/thread-219-1-1.html）里面的资料。

和静电势同等重要的是笔者提出的范德华势，适用于考察范德华作用主导的相互作用，可以由Multiwfn计算，见《谈谈范德华势以及在Multiwfn中的计算、分析和绘制》（//www.umsyar.com/551）。

文/Sobereva@北京科音   2025-Feb-21

1 前言

《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）里介绍的笔者提出的图形化展现相互作用的方法IGMH已被广为使用，之前我写的《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》（//www.umsyar.com/727）和《谈谈pi-pi相互作用》（//www.umsyar.com/737）里涉及到了IGMH的等值面面积，至少对pi-pi作用来说它和相互作用强度有密切的正相关性。有不少读者都问我怎么得到面积，本文就专门说一下。计算面积的同时还会顺带得到等值面内包围的体积。此文说的IGMH等值面具体是指IGMH方法里定义的δg_inter函数的等值面。如果读者不熟悉IGMH，应先把//www.umsyar.com/621看了。读者应使用Multiwfn官网//www.umsyar.com/multiwfn上的最新版本以免和本文的情况不符。不了解Multiwfn者参看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。

下面介绍两种做法，第一种方法是使用Multiwfn的定量分子表面分析（主功能12）计算IGMH等值面面积和体积，这种情况只适合存在一个等值面，且这个等值面就是你要研究的等值面的情况。另一种方法更为普适，需要借助免费的ChimeraX程序载入Multiwfn产生的cub文件显示IGMH等值面，可以测量图中任意一个等值面的面积和体积，对于《使用IRI方法图形化考察化学体系中的化学键和弱相互作用》（//www.umsyar.com/598）介绍的IRI函数、《谈谈范德华势以及在Multiwfn中的计算、分析和绘制》（//www.umsyar.com/551）介绍的范德华势等各种函数也都可以这么测量。

2 使用Multiwfn的定量分子表面分析功能计算IGMH等值面的面积和体积

这一节以18碳环与C60富勒烯的复合物为例演示怎么直接用Multiwfn得到IGMH等值面的面积和体积。前述的《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》介绍的笔者的Chem. Eur. J., 30, e202402227 (2024)研究中得到的此体系的波函数文件C60-C18.wfn在//www.umsyar.com/attach/738/file.rar里提供了。此体系的IGMH图如下所示（δg_inter函数等值面为0.002 a.u.）

首先将Multiwfn目录下的settings.ini里的iuserfunc设为91，这代表把用户自定义函数（user-defined function）设为IGMH方法的δg_inter函数。在Multiwfn手册2.7节里有可用的用户自定义函数的完整列表。之后启动Multiwfn，载入C60-C18.wfn，之后依次输入
1000   //隐藏功能
16  //定义片段
2  //定义两个片段
1-18  //18碳环里的原子序号范围
19-78   //富勒烯里的原子序号范围
12  //定量分子表面分析功能
1  //设置用于定义表面的函数
2  //某个实空间函数
100  //用户自定义函数
0.002  //定义表面用的等值面数值
6  //开始分析，不考虑被映射的函数
接下来程序开始计算δg_inter的格点数据，过一会儿，在屏幕上看到以下信息
 Volume:    69.13431 Bohr^3  (  10.24465 Angstrom^3)
 Estimated density according to mass and volume (M/V):  151.8505 g/cm^3
 Overall surface area:         323.36643 Bohr^2  (  90.55182 Angstrom^2)

在后处理菜单选-3可以观看当前考察的δg_inter函数的0.002 a.u.等值面，如下所示，确实就是前面给出的笔者的论文Chem. Eur. J., 30, e202402227 (2024)里的那个等值面。其体积是上面显示的10.24 Å^3，面积是90.55 Å^2，和文中报道的一致。

顺带一提，在后处理菜单中还可以选-2将当前算出来的δg_inter的格点数据导出为surf.cub，之后可以被第三方程序可视化和分析。

还值得一提的是进入主功能12的时候可以看到选项3 Spacing of grid points for generating molecular surface用来设置格点间距，默认是0.25 Bohr，数值设得越小计算耗时越高而统计精度越高。

Multiwfn对各种实空间函数（包括从外部的.cub等格点数据文件读入的）都可以利用定量分子表面分析功能计算其等值面的面积和体积，另一个使用例子见《使用Multiwfn计算轨道的体积》（//www.umsyar.com/734）。

3 使用Multiwfn结合ChimeraX获得IGMH等值面的面积和体积

这一节以2-pyridoxine和2-aminopyridine的二聚体为例演示利用ChimeraX程序得到特定的IGMH等值面的面积和体积。examples\2-pyridoxine_2-aminopyridine.wfn是Multiwfn程序包里自带的这个体系的波函数文件，在0.01 a.u.的δg_inter等值面数值下分子间的IGMH等值面图如下所示，可见有两个等值面，此例分别获得它们的面积和体积

首先照常对这个二聚体做IGMH分析。启动Multiwfn，然后依次输入
examples\2-pyridoxine_2-aminopyridine.wfn
20  //弱相互作用可视化分析
11  //IGMH分析
2  //两个片段
1-12  //第1个片段
c  //其它原子作为第2个片段
4  //设置格点间距（格点分布覆盖整个体系）
0.2  //格点间距为0.2 Bohr
3  //导出格点数据
当前目录下有了dg_inter.cub，是δg_inter的cub文件。

之后退回到Multiwfn主菜单，输入xyz后按回车，再输入2-pyridoxine_2-aminopyridine.xyz，当前目录下就得到了记录当前结构的2-pyridoxine_2-aminopyridine.xyz文件。

去https://www.rbvi.ucsf.edu/chimerax/download.html下载ChimeraX并安装。本文用的是ChimeraX 1.9。

启动ChimeraX，将2-pyridoxine_2-aminopyridine.xyz和dg_inter.cub依次拖入程序界面载入，然后左右随便拖动一下下图红箭头所示的竖杠，激活这个等值面的设置，然后再在下图蓝箭头所示的文本框里输入要用的等值面数值0.01，之后看到的等值面就是下图这样

之后选择窗口上方的Tools - Volume Data - Measure and Color Blobs，之后按住Alt键并点击你要考察面积和体积的那个等值面，那个等值面就被自动着色了，并且在界面右下角显示了其体积和面积，如下所示，分别为0.31 Å^3和3.28 Å^2。

值得一提的是，由于ChimeraX和Multiwfn的定量分子表面分析功能产生等值面的算法不同，因此得到的等值面的面积和体积会有轻微差异。例如用ChimeraX对上一节的18碳环和富勒烯之间的δg_inter=0.002 a.u.的等值面进行测量，得到的面积是89.90 Å^2，体积是10.65 Å^3，面积和Multiwfn给出的90.55 Å^2相差0.7%。

笔者之前还录过一个视频《使用Multiwfn和ChimeraX绘制自定义着色的电子定域化函数(ELF)等值面图》（https://youtu.be/vC48iEB8PwI、https://www.bilibili.com/video/av85684420）演示ChimeraX里的和等值面显示、测量相关的操作，感兴趣的读者可以看看。

文/Sobereva@北京科音

First release: 2025-Feb-18   Last update: 2025-Feb-19

0 前言

pi-pi相互作用（pi-pi interaction）是化学体系中很常见而且很重要的一种弱相互作用，本文全面谈谈这种相互作用的各方面特征，以令读者对其有充分的认识、能在自己的研究中正确分析讨论。我在互联网上广泛答疑时也时不时看到有关于pi-pi作用的非常错误的理解，比如有人被一些文章误导，竟以为pi-pi作用来自于轨道相互作用或静电相互作用，而且这种说法还广泛以讹传讹，笔者希望通过此文以正视听。此外，氢键、范德华作用等概念在结构化学书里普遍都有，但pi-pi作用鲜有提及，我也推荐相关书籍作者和结构化学教师参考此文的内容将pi-pi作用纳入书籍和教学。本文中大量涉及非常流行的波函数分析程序Multiwfn，相关信息见《Multiwfn FAQ》（//www.umsyar.com/452）。

我有不少研究文章都和pi-pi作用有密切联系，是pi-pi作用研究的典型范例，非常推荐读者阅读并欢迎引用：
• Theoretical Insight into Complexation Between Cyclocarbons and C60 Fullerene, Chem. Eur. J., 30, e202402227 (2024)。在《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》（//www.umsyar.com/727）中有详细介绍，专门研究了不同尺寸的碳环和富勒烯之间的pi-pi作用
• Intermolecular interaction characteristics of the all-carboatomic ring, cyclo[18]carbon: Focusing on molecular adsorption and stacking, Carbon, 171, 514-523 (2021)。在《全面探究18碳环独特的分子间相互作用与pi-pi堆积特征》（//www.umsyar.com/572）中有详细介绍，其中涉及两个18碳环之间的pi-pi作用
• Molecular assembly with a figure-of-eight nanohoop as a strategy for the collection and stabilization of cyclo[18]carbon, Phys. Chem. Chem. Phys., 25, 16707 (2023)。在《8字形双环分子对18碳环的独特吸附行为的 、波函数分析与分子动力学研究》（//www.umsyar.com/674）中有详细介绍，其中涉及18碳环与OPP双环分子的pi-pi作用
• Comment on “18 and 12 – Member carbon rings (cyclo[n]carbons) – A density functional study”, Mat. Sci. Eng.: B, 273, 115425 (2021)。在《18碳环（cyclo[18]carbon）与石墨烯的相互作用：基于簇模型的研究一例》（//www.umsyar.com/615）中有详细介绍，其中涉及18碳环与石墨烯的pi-pi作用

1 pi-pi作用的基本特征

pi-pi作用的定义和说法比较乱，这里给出我认为最严格准确的定义：“pi-pi作用是相距较近的两个片段上彼此朝向相对的pi电子之间的独特的色散吸引作用”。这里做几个注释，结合后文的实例读者会了解得更充分：

(1)pi-pi作用可以是分子间的也可以是分子内的，满足以上定义即可
(2)诸如苯分子里面的pi电子之间的作用不叫pi-pi作用，那属于共享电子作用
(3)两套pi电子的分布必须近乎彼此相对才可能算pi-pi作用。而诸如两套pi电子近乎肩并肩挨着就不能算pi-pi作用，像是不能说环丁二烯里面两套近乎定域的pi电子之间是分子内pi-pi作用
(4)“相距较近”一般是在4.0-4.5埃以内。也不是说更远距离就完全没有pi-pi作用了，只不过由于色散作用随作用距离r呈-1/r^6快速衰减行为，因此距离稍微一远pi-pi作用就非常弱了，就不太值得一提了。但相互作用的片段间距离也不能太近，若显著小于相接触的原子的范德华半径之和，则显著的位阻互斥作用会远大于pi-pi吸引作用，使得pi-pi作用也相对不值得一提
(5)pi-pi作用最常出现在碳原子间，因为碳最容易带显著的pi电子。碳的Bondi和CSD范德华半径分别为1.70和1.77埃。如果没有特殊的因素影响相互作用距离的话，在平衡结构（势能面极小点结构）下，C-C间的pi-pi作用出现在3.4-3.6埃左右是最常见的
(6)pi-pi作用属于弱相互作用和非共价相互作用范畴。虽然名义上算作弱相互作用，但实际强度可大可小，直接取决于作用面积，详见后文
(7)有一个词叫pi-pi堆积（pi-pi stacking），这个词往往和pi-pi作用混用。但在我来看，这个词更适合用来形容由于pi-pi作用的吸引效果使得相互作用的两个部分发生紧密结合从而产生的彼此堆积的结构特征

下图就是一个再典型不过的存在显著pi-pi作用的体系，是晕苯的二聚体。高精度理论计算的极小点结构下平面间距大约在3.45埃

苯二聚体有不同构型，彼此垂直的T型二聚体是能量最低构型，还有一种能量与之非常接近的局部极小点构型是平行位错构型，属于pi-pi堆积结构。按照《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）的做法绘制出的几何结构结合pi电子等值面图如下所示，两个分子的pi电子的分布非常清楚直观，可以看到彼此相对。注意并非必须两个分子的所有pi电子都彼此对着才算pi-pi作用，像下图这样哪怕只要一部分对着也同样算pi-pi作用。而如果完全没有对着、完全错开了，那就不算pi-pi作用了。更多的展现pi电子分布的图见笔者的Theor. Chem. Acc., 139, 25 (2020)一文。

2 pi-pi作用的物理本质

我在《谈谈“计算时是否需要加DFT-D3色散校正？” 》（//www.umsyar.com/413）提到过，原子间相互作用从物理本质上可以基本分为静电、色散、交换、Pauli互斥、轨道相互作用；我在《谈谈“计算时是否需要加DFT-D3色散校正？”》（//www.umsyar.com/413）里也说过物理本质和相互作用表象之间的关系。pi-pi作用类似于氢键、卤键等，是属于表象层面的，用来描述一种具有特定特征的相互作用，其内在本质就是前面说的来自于相距较近的pi电子之间的色散作用，这一点在Grimme的一篇经典的讨论pi-pi作用本质的文章Angew. Chem. Int. Ed., 47, 3430 (2008)中通过严格的理论计算和分析对比已经论证得十分严格充分，而且也早已是 领域内行人的共识。

有的文章试图从轨道相互作用解释pi-pi堆积出现的本质，这是极其有误导性的，却还广泛以讹传讹。2007年的时候wiki上的词条说pi-pi作用来自于pi共轭体系间p轨道的重叠，前述Angew文章的脚注中专门对此进行了批判，指出没有任何 计算能支持这种说法。这也不难理解，本来两个分子的pi轨道之间重叠程度极低（可以用《分子间轨道重叠的图形显示和计算》//www.umsyar.com/163介绍的方法考察），远低于形成共价键时原子轨道间的重叠程度，因此有效重叠不足。而且两个分子的pi占据轨道之间混合会产生全占据的成键和反键轨道，二者效果充分抵消了，对成键没贡献；而若一个分子的pi占据轨道和另一个分子的pi非占据轨道混合，意味着要出现分子间pi->pi*电子转移，这通常又没有明显能使体系能量变得更低的内在驱动力（不像常见的较近距离的n->pi*超共轭往往可以明显降低体系能量）。此外，按照《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）说的，用Multiwfn对pi-pi堆积的两个分子间去计算衡量两个部分之间有效共享电子对数的Mayer键级、模糊键级或离域化指数（DI），会发现数值很接近0，比如《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》（//www.umsyar.com/727）介绍的文章里对18碳环与富勒烯之间算的Mayer键级仅为可忽略不计的0.04，也充分体现出pi-pi作用中的轨道相互作用成份基本可忽略不计。所以无论从哪个角度，pi-pi堆积的出现都不是轨道相互作用驱动的，顶多是在色散作用驱动下，pi-pi堆积结构出现后顺带出现了点可忽略的轨道相互作用而已。所以pi-pi作用在主流学术界被归为非共价相互作用范畴，这是完全正确的。然而有一些论文试图从轨道相互作用分析pi-pi作用的特征，还摆出一堆轨道图，比如Phys. Chem. Chem. Phys., 15, 9397 (2013)就是典型文章，切勿把这类文章太当回事，很多都是瞎讨论、尬讨论。

pi-pi堆积的体系有一个普遍特点是在极小点结构下，pi-pi作用区域的原子间通常不是正好对着，而是相互错位。下图是晕苯二聚体的俯视图，其中一个晕苯用红色显示，可以很清楚看到错位特征，即一层的碳对着另一层的六元环中心。一层层堆叠的石墨中每一层之间也同样是这样错位的。出现自发错位有不同的解释并且存在一定争议，有人认为是因为错位的结构下比严格面对面的结构下的原子间的Pauli互斥更低，例如Chem. Sci., 11, 6758 (2020)持这种观点。也有人认为是因为错位的结构下的静电相互作用能更负（静电吸引作用更强），在Phys. Chem. Chem. Phys., 24, 8979 (2022)中作者通过能量分解论证了这一点并反驳了Chem. Sci.那篇文章。静电作用引起位错的原因从逻辑上容易理解：形成pi-pi堆积的两部分都有丰富的pi电子，面对面构型下pi电子间静电互斥会较大，错开后互斥自然就没那么强了。虽然笔者更同意静电效应是位错出现的主导，但我也不否认降低Pauli互斥也可能是产生位错的另一个驱动力，尽管相对次要。

顺带一提，不止是pi-pi堆积二聚体，还有很多其它体系都是色散吸引作用驱动分子间结合，而静电作用进一步影响几何结构。《静电效应主导了氢气、氮气二聚体的构型》（//www.umsyar.com/209）介绍的笔者的研究就是很好的例子，H2、N2的各种二聚体的相互作用能通过能量分解分析可以发现都是色散作用为主，但不同构型下的静电作用的差异则直接影响不同构型的稳定性乃至存在与否，并且笔者发现通过静电势互补原理可以给出很直观的解释。这种静电势互补的分析方法也在《全面探究18碳环独特的分子间相互作用与pi-pi堆积特征》（//www.umsyar.com/572）介绍的笔者的论文中也用来解释为什么18碳环的pi-pi堆积二聚体是错位的。平面环状结构的18碳环具有独特的平面内和平面外两套pi电子，此文中确认了18碳环能够形成分子间彼此平行的pi-pi堆积二聚体，下图是按照《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图》（//www.umsyar.com/443）绘制的二聚体极小点结构下的单体的静电势着色的分子范德华表面的叠加图，可以看到两个碳环表面静电势为正和为负的区域在错位的结构下是以互补的方式重叠的，很清晰地解释了错位的产生原因，即这种构型下的静电吸引作用比一个个碳原子彼此精确对着的时候更强。

要注意的是，有些文章作者极度夸大了静电作用对于pi-pi堆积形成的作用，不仅认为平行位错构型的出现是因为此时静电吸引作用最有利，还鼓吹pi-pi堆积结构的形成的本质就是静电吸引作用，这是大错特错！最典型的这种文章就是Chem. Sci.,3 , 2191 (2012)，这篇文章的发表严重误导了很多读者。此文虽然是在前述的Grimme的pi-pi作用的研究之后写的，而且文中还引了那篇文章，但给我的感觉是此文的作者完全没好好看那篇文章，似乎也完全不懂什么叫电子相关作用，而依然基于古老且过于简单（忽略了pi电子结构和穿透效应）的Hunter和Sanders的四极矩观点盲目、主观、武断、信誓旦旦地解释pi-pi堆积的位错结构形成的本质，而完全没有可信、严格的理论和计算数据作为依据，各种自说自话，并且还错误地解释一些文献里的理论计算数据，甚至还不承认pi-pi作用的概念。作者还执拗地认为非得是存在精确面对面堆积才能说明pi-pi作用的存在，并因为实际pi-pi堆积体系的极小点结构不是面对面堆积就被他用来否认存在pi-pi作用。此文文中给出了下图，令一些初学者不明觉厉，误以为不仅清晰解释了为什么平行位错结构比精确面对面结构更有利，还误以为这种二聚体的形成主要靠的就是静电吸引作用，完全没分清楚pi-pi堆积二聚体的形成中色散作用和静电作用的主次关系。还有不少其它文章如J. Chem. Soc., Dalton Trans., 2000, 3885也是类似地对pi-pi作用存在严重误解、给出误导性的示意图，读者看到这种文章时切勿当回事！

能量分解分析在pi-pi作用的研究中应用得非常普遍，对于认识其本质非常有帮助。这类方法将总的相互作用能分解成不同物理成份。例如将《使用sobEDA和sobEDAw方法做非常准确、快速、方便、普适的能量分解分析》（//www.umsyar.com/685）介绍的笔者提出的流行的sobEDAw能量分解方法用于苯的两种二聚体构型可得到如下结果（计算级别：B3LYP-D3(BJ)/6-311+G(2d,p)并考虑counterpoise，单位为kcal/mol），可见色散项ΔE_disp远比静电作用项ΔE_els更负，轨道相互作用项ΔE_orb更是微乎其微，故色散作用对苯的分子间结合起到主导作用。特别是在平行位错的结构下的ΔE_disp比T构型下的明显负得多（同时ΔE_els的大小还变小），充分体现出平行位错结构下才有的pi-pi作用的本质是色散作用，这和上文对pi-pi作用本质的讨论相一致。平行位错结构下，色散作用对总吸引作用项的贡献百分比是：-7.6/(-7.6-0.8-1.6)*100=76%。

顺带一提，pi-pi堆积的结构在现实环境中一般很容易发生分子间相对滑移（除非有额外的位阻效应等阻碍），这是由于滑移导致能量变化很小。Carbon, 171, 514-523 (2021)文中做18碳环二聚体的势能面扫描、Phys. Chem. Chem. Phys., 24, 8979 (2022)中做的苯分子平行二聚体的势能面扫描都体现了这一点。并且由于滑移方向的势能面很缓，哪怕室温下也非常容易出现滑移。在Carbon, 171, 514-523 (2021)文中我给出了18碳环二聚体在100K下的4 ps的从头算动力学的轨迹叠加图（消除了下方的18碳环的整体运动以着重表现相对运动），如下所示，分子位置根据模拟时间按照蓝-白-红变化，可见在模拟过程中的确出现了显著的滑移

3 pi-pi作用的强度

pi-pi堆积的极小点结构下的每一对近距离接触的原子的pi-pi作用都是非常弱的，就是普通范德华作用的程度，但是当涉及pi-pi作用的原子很多时，就可以达到很高的强度，甚至能达到近乎化学键的程度。最强的pi-pi作用就是从微观来看无限延展的石墨中的层间pi-pi作用了。下面看一些实例，以令读者更好地理解pi-pi作用的强度特征。

Angew. Chem. Int. Ed., 47, 3430 (2008)文中计算了T型苯二聚体（下图a）、平行位错苯二聚体（下图b）、环己烷二聚体（下图c和d是两个视角）的相互作用能，并且还将它们的环数n从1逐渐加到4（分别对应下图e、f、g的结构），以考察了三种作用类型的相互作用能随环数的变化，分别对应下图的aromatic T-shaped、aromatic stacked和saturated stacked三条曲线。由图可见，当pi-pi作用涉及的原子很少时，即苯二聚体，由于pi-pi作用还不够强，因此平行位错二聚体的能量比静电吸引作用更强的苯T型二聚体略微更高，相互作用强度比起无pi-pi作用特征的环己烷二聚体没优势。但随着环数增加，pi-pi作用强度不断增大，使得具有pi-pi堆积结构的二聚体的相互作用变得显著强于不具备这种特征的二聚体。四并苯的二聚体的相互作用能已达到-16 kcal/mol（67 kJ/mol），这已经是水二聚体这种普通强度氢键作用能的3倍左右了。

可见，前面提到的Chem. Sci.,3 , 2191 (2012)那篇文章拿两个苯及其衍生物在液体和晶体中普遍缺少平行位错构型的出现，以及平行位错的苯二聚体不是能量最低构型，以此鼓吹pi-pi作用不是一种客观存在的作用，这明显是极其狭隘的。

对本文第1节的两个晕苯的二聚体，在DLPNO-CCSD(T)级别下算出来的相互作用能达到-20 kcal/mol，也是相当强了。

前面给出的18碳环二聚体的相互作用能，在Carbon, 171, 514-523 (2021)中我用ωB97X-V/def2-QZVPP结合counterpoise校正的计算结果为-9.2 kcal/mol（-38.5 kJ/mol），达到平行位错的苯二聚体的相互作用能的大约三倍，无疑也是非常显著的pi-pi作用了。

不是只有平面区域之间才可能有pi-pi作用，例如下面三个体系的pi-pi作用区域都是有显著曲率的。下图左边的结构是Org. Lett., 17, 5292 (2015)中合成出的主-客体复合物晶体的一部分，由于客体分子与富勒烯之间的pi-pi作用面积颇大，B97D/QZVP级别计算的相互作用能达到了-50 kcal/mol，已经是弱化学键的强度了。下图右侧是碗烯的pi-pi堆积二聚体，相互作用的研究见J. Phys. Chem. A, 116, 11920 (2012)。

4 图形化展现pi-pi作用区域

《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）以及《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（//www.umsyar.com/667）中的综述详细介绍的笔者提出的IGMH方法用于图形化展现片段间相互作用效果极佳，已被广为使用。IGMH用于展现自定义片段间的pi-pi作用、判断pi-pi作用是否存在尤为好用，因此在此给出一些简单例子予以体现。其方法的原理、细节、计算操作见上面的文章。如果你不知道pi-pi作用可能出现在哪，不方便定义片段，则应当用《使用IRI方法图形化考察化学体系中的化学键和弱相互作用》（//www.umsyar.com/598）介绍的笔者的IRI方法，可以把体系中所有相互作用区域全都展现出来，也包括化学键作用区域。

不是只有纯碳的pi区域之间才可能有pi-pi作用，pi-pi作用也可以涉及到其它原子。众所周知，DNA结构中相邻层的碱基与碱基之间是平行堆叠的，再加上碱基区域有大量pi电子，所以这也属于典型的pi-pi作用。下图是从DNA中截取的堆叠的GCGC碱基对，两层分子各定义为一个做IGMH分析的片段，两层之间的绿色的IGMH方法定义的delta_g_inter函数的0.004 a.u.等值面非常清晰、优雅地展现出了pi-pi作用的主体区域。图中还有个面积较小的等值面出现在O...O之间，虽然C=O键的O有pi电子，但这俩O的pi电子不是彼此对着的，所以这就只能算是普通的范德华作用区域了。

下图是《8字形双环分子对18碳环的独特吸附行为的 、波函数分析与分子动力学研究》（//www.umsyar.com/674）介绍的笔者的论文中的一张图，研究的是OPP双环分子结合两个18碳环后的复合物。图中绿色等值面直观地展现出在什么区域18碳环与OPP主体分子间有显著的pi-pi作用。可见pi-pi作用主要出现在OPP和18碳环近距离接触的体系的两端，而OPP靠中间区域离18碳环较远因此缺乏pi-pi作用。顺带一提，18碳环是极少有的同时具有平面内（in-plane）和平面外（out-of-plane）两套全局离域的pi电子的体系，当前体系中18碳环主要凭借的是平行于环方向分布的pi电子与OPP大环上的垂直于苯环的pi电子产生pi-pi作用，这是极为新颖、独特的pi-pi作用形式。如果读者对碳环类体系感兴趣，非常建议进入//www.umsyar.com/carbon_ring.html观看笔者发表的所有碳环相关的理论研究和对应的深入浅出的评述文章。

上图中还对利用Multiwfn的基于力场的能量分解（EDA-FF）功能计算了各个原子对分子间相互作用的贡献并对原子进行了着色，展现出了更丰富的信息。这种分析见《使用Multiwfn做基于分子力场的能量分解分析》（//www.umsyar.com/442）。

在《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》（//www.umsyar.com/727）介绍的文章中，笔者全面研究了碳环与富勒烯的相互作用，其中给出了下图，展示了不同尺寸的两个碳环与富勒烯形成的三聚体。并且为了直观区分富勒烯-碳环和碳环-碳环两种pi-pi作用，分别将对应的等值面用绿色和青色着色。可见此图极好地将体系中两种pi-pi作用出现的主要位置展现了出来，比起在文中只是给出个几何结构图信息丰富得多。

标准的IGMH图是通过sign(λ2)ρ函数对等值面着色的，ρ是电子密度。按照IGMH方法的标准的色彩刻度着色的话，ρ接近0的区域的等值面都是绿色。由于pi-pi作用区域、普通范德华作用、极弱的氢键（如C-H作为氢键给体时）等情况的作用区域的电子密度都很低（0.01 a.u.以内），因此相应的标准方式着色的IGMH等值面都会是绿色或很接近绿色。因此解释IGMH图的等值面时必须结合不同特征的弱相互作用的基本定义，不能单纯根据颜色去乱解释。比如网上有人问我“下面这个图里红圈部分是pi-pi作用吗？”，这实在是匪夷所思的问题。跨越那个等值面和右下方苯环的pi电子区域直接作用的是一个氢原子，氢原子又没pi电子，怎么可能是pi-pi作用！？而且右下方的苯环也明显不可能和左上方的苯环有pi-pi作用，一方面二者离得老远，另一方面二者的pi电子区域也根本不相互对着，显然无论如何也不可能算是pi-pi作用。

那么上图中红圈里的等值面算什么作用？明显这是pi-氢键，C-H作为氢键给体，其中氢带一定正电荷（这一点用H的原子电荷就能体现，原子电荷的相关知识见《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！》//www.umsyar.com/714这篇综述），右下角的苯环上的pi电子令苯环上方显负电性而作为氢键受体（用静电势可直接体现，参见http://bbs.keinsci.com/thread-219-1-1.html里汇总的静电势相关资料和我的相关博文）。

下面再举一个IRI图展现pi-pi作用区域的例子。Multiwfn的原文之一J. Chem. Phys., 161, 082503 (2024)里给出了144-轮烯的LOL-pi函数的等值面图，如下图左侧所示，此图清晰直观地展现出了这个全局大共轭体系的pi电子的主要离域路径。根据此图展现的pi电子的分布情况，凭直觉就知道这个体系里肯定存在pi-pi作用。下图右侧是对这个体系的其中一个局部区域绘制的IRI等值面图，蓝色的等值面展现出了化学键作用区域，绿色的等值面清楚直观地展现出了pi-pi作用区域。可见这个体系非常有趣，pi-pi作用区域绵延不断贯通整个体系！也正是有分子内pi-pi作用的存在，此体系才能形成螺旋状结构，要不然就散了（如同用不能描述色散作用的理论方法优化DNA结构时结构会散掉）。

5 衡量pi-pi作用强度的方法

这一节说一下如何衡量pi-pi作用强度。最简单的方法就是计算pi-pi堆积二聚体的相互作用能，作用能越负说明作用越强。但如果要讨论二聚体的热力学稳定性，则需要计算结合自由能，溶剂环境中还得考虑溶剂模型。相关知识参见《谈谈分子间结合能的构成以及分解分析思想》（//www.umsyar.com/733）。复合物AB在特定结构下，A和B的相互作用能的最常规的计算方法就是用E(AB)-E(A)-E(B)方式手动计算，做sobEDAw能量分解（//www.umsyar.com/685）时也会顺带给你相互作用能。

以上述方法算相互作用能的一个问题是，如果两个分子之间不仅仅有pi-pi作用，还有其它作用（如氢键），那么得到的只是总相互作用能。如果想只得到pi-pi作用能，有几个办法可以用：
(1)用《使用Multiwfn做基于分子力场的能量分解分析》（//www.umsyar.com/442）介绍的EDA-FF能量分解方法，将两个分子的pi作用区域定义为两个片段，让Multiwfn给出基于分子力场的这两个部分的相互作用能，取其中的范德华作用能（即色散作用和交换-互斥作用之和）。如果你只需要色散部分，还有另一种做法，见《使用Multiwfn图形化展现原子对色散能的贡献以及色散密度》（//www.umsyar.com/705）。这两种做法还都可以给出具体原子产生的贡献，并可以对原子进行着色以便通过图像直观展现和考察
(2)如果分子间同时有pi-pi作用和氢键，且其它作用可忽略不计，可以按《透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍》（//www.umsyar.com/513）介绍的方法使用Multiwfn做拓扑分析估计出氢键作用能，从总相互作用能中扣掉之作为pi-pi作用能
(3)对体系进行恰当改造，基本保留每个分子参与pi-pi作用的部分，然后将总相互作用能近似当成pi-pi作用能。

如果是分子内的pi-pi作用能的估计，还可以参考《计算分子内氢键键能的几种方法》（//www.umsyar.com/522）里的说明举一反三处理。

在特定条件下，pi-pi作用强度和IGMH的等值面面积有密切的正相关性。《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》（//www.umsyar.com/727）介绍的文章中给出了下图，是不同原子数的碳环与富勒烯之间的相互作用能和描述pi-pi作用的IGMH等值面面积的关系，可见有挺理想的线性关系。通过这种关系，可以直接根据pi-pi作用区域的等值面面积估计对应的pi-pi相互作用能。这种面积的计算方法见《计算IGMH等值面的面积和体积的方法》（//www.umsyar.com/738）。

Chem. Commun., 48, 9239 (2012)提出的LOLIPOP方法值得一提。基于单个分子的波函数文件，就可以用Multiwfn非常容易地计算体系中的各个六元环的LOLIPOP指数，此值越小说明这个环发生pi-pi堆积的能力越强，原文对不同体系以苯分子作为探针分子进行了测试证实了这一点。Multiwfn手册3.100.14节有LOLIPOP的详细介绍，4.100.14节有计算实例。

6 不同理论方法描述pi-pi作用的能力

色散作用的更深层本质是电子的库仑相关作用，因此做 或第一性原理计算时，若想描述好pi-pi作用，用的方法必须能描述好库仑相关，显然高精度后HF方法都没问题，比如CCSD(T)在计算pi-pi作用上可以算是金标准。至于低级别的后HF方法MP2则倾向于明显高估pi-pi作用，绝对不要用。

对于特别常用的DFT，描述pi-pi作用的能力基本等同于描述普通色散作用的能力，如果泛函原本描述色散作用烂（如PBE、PBE0）或者完全不能描述（如B3LYP），就必须带色散校正，参考《谈谈“计算时是否需要加DFT-D3色散校正？” 》（//www.umsyar.com/413）、《DFT-D色散校正的使用》（//www.umsyar.com/210）、《谈谈 研究中什么时候用B3LYP泛函优化几何结构是适当的》（//www.umsyar.com/557）。本来就带色散校正的wB97M-V、wB97X-D3等直接就能很好描述pi-pi作用。M06-2X描述pi-pi作用虽然定性正确但不算太好，加了零阻尼DFT-D3色散校正后对pi-pi作用的描述有明显改进，但还是不如B3LYP-D3(BJ)，对比测试见考虑了很多pi-pi作用体系的L7测试集的原文（J. Chem. Theory Comput., 9, 3364 (2013)）。双杂化泛函由于带有MP2项，所以都有描述pi-pi作用的能力，但一般在考虑了色散校正后才会变得足够好，如revDSD-PBEP86-D3(BJ)。

实际上DFT-D那种形式的色散校正在原理上对于描述pi-pi作用并非很理想，因为pi电子不是绕着原子核球对称分布的，而DFT-D校正能的公式依赖的只是原子间距离（这里不考虑三体校正项），没体现出pi电子在原子核周围的具体分布特征。不过这倒也不是明显问题，常用的DFT-D3、DFT-D4色散校正对于描述pi-pi作用从实际效果上来看并没有什么问题。

在半经验方法层面，专门考虑了对色散作用描述的GFN2-xTB、PM6-D3H4X'、PM6-ML在表现pi-pi作用方面优秀，见J. Chem. Theory Comput., 21, 678 (2025)里3图基于L7测试集的测试，PM6-D3和PM7只能算是定性正确。至于没专门考虑对色散作用描述的诸如PM6、AM1等方法则完全失败。

主流的分子力场，如GAFF、AMBER、CHARMM、OPLS-AA、MMFF94等，对pi-pi作用的描述虽然跟像样的 方法比算不上出色，但至少也算定性正确，因为它们都有描述色散作用的能力。J. Chem. Inf. Model., 49, 944 (2009)的测试专门体现了这点（不过这篇文章也有不少漏洞和槽点）。

7 疏水作用与pi-pi作用的关系

众所周知，疏水效应使得水环境下非极性物质倾向于发生聚集，本质是溶剂的熵效应。两个石墨烯片段在水中会自发堆积在一起，这算疏水作用还是pi-pi作用？实际上二者都有，对于堆积结构的形成起到协同作用。疏水效应更为普遍，无论两个溶质的接触区域是否有pi电子，只要溶质是基本无极性的，在水中都有疏水作用促使它们发生结合。而对于pi电子区域暴露的两个溶质，疏水效应则在pi-pi作用的基础上进一步促进了它们的pi-pi堆积结构的出现。如《全面揭示各种碳环与富勒烯之间独特的pi-pi相互作用！》（//www.umsyar.com/727）介绍的文章的理论计算所示，在水环境下碳环与富勒烯之间的结合自由能的大小显著大于在真空下，充分体现了这一点。

顺带一提，前述的Chem. Sci.,3 , 2191 (2012)一文居然误以为溶剂环境下pi共轭体系间出现堆积结构仅仅是因为溶剂效应，并似乎试图靠这个否认真空环境下也存在pi-pi相互作用，甚至说the terms "pi-stacking" or "pi–pi interactions" do not describe any physically meaningful interaction，实在是难以理喻！！！PS：这样的文章能通过Chem. Sci.的审稿真是离谱。

8 总结&判断pi-pi作用的标准

本文系统地对pi-pi作用的各个方面进行了介绍，包括其基本特征、物理本质、强度范围、图形化展现方法、考察强度的方法、理论计算方法的精度、疏水作用与它的关系。通过本文，读者应该已经对pi-pi作用有了较全面的了解，能够进行正确的分析讨论，并且能认识到哪些文章或书籍里的说法是有误导性的。我也推荐读者接下来阅读本文一开始提到的笔者的一系列和pi-pi作用有关的研究文章和对应的介绍博文。

最后再总结一下pi-pi作用的常规判断标准，便于读者能确切判断哪些作用算是pi-pi作用。通过下面第1、2条就可以进行粗略判断，3、4、5可以作为进一步检验。一般意义的pi-pi作用应当能同时满足所有条件
(1)相互作用的两部分都有pi电子且其分布彼此相互对着。如果拿不准有没有pi电子分布、分布朝向如何，可以按照《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）介绍的方法直接把pi电子密度等值面图画出来，一目了然
(2)相互作用的pi电子之间离得不太远。比如明显超过5埃的直接就可以忽略了
(3)用Multiwfn做IGMH分析（如果得不到波函数文件或体系太大难以算得动的话，可以改用极便宜且只依赖于原子坐标的mIGM），在等值面数值调到诸如0.003 a.u.这样较小的值时，在预期出现pi-pi作用的区域能明显看到等值面
(4)两部分之间的Mayer键级或模糊键级或离域化指数非常小（远小于0.1），故而轨道相互作用可基本忽略。注：如《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）所述，Multiwfn在主功能9里计算键级之前可以用选项-1定义两个片段，之后做键级计算时会给出两个部分之间的总键级，即片段间每一对原子的键级的总和
(5)使用sobEDAw等能量分解方法分析，色散项应当占所有吸引项总和的大部分，轨道相互作用项应当占比微小

文/Sobereva@北京科音  2025-Jan-1

1 前言

有人在思想家公社QQ群（//www.umsyar.com/QQrule.html）里问怎么计算分子轨道的体积，并且贴了一张轨道等值面图。只要计算出这个等值面里包围的体积就可以达到这个目的（虽然也可以有其它方式的定义，但没这么简单直观）。这里我介绍一下怎么利用Multiwfn程序的定量分子表面分析功能来实现这个目的。定量分子表面分析功能之前在《使用Multiwfn的定量分子表面分析功能预测反应位点、分析分子间相互作用》（//www.umsyar.com/159）有简要介绍，结合Multiwfn的极度灵活的设计，这个功能还能做很多其它的分析，正如本文所展示的。

如果读者不熟悉Multiwfn，应当阅读《Multiwfn FAQ》（//www.umsyar.com/452）、《Multiwfn入门tips》（//www.umsyar.com/167）、《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）。Multiwfn可以在//www.umsyar.com/multiwfn免费下载，不要用上古版本。

2 实例

这里以Multiwfn计算下面这个D-pi-A体系的HOMO轨道的体积为例，图中的轨道等值面数值为0.05，此体系的波函数文件是Multiwfn目录下的examples\excit\D-pi-A.fchk。注意等值面包围的体积直接依赖于等值面数值的选取，横向对比时必须统一。

启动Multiwfn，然后输入
examples\excit\D-pi-A.fchk
5  //计算格点数据
4  //轨道波函数
h  //HOMO
4  //设置格点间距
0.25  //0.25 Bohr是精度和耗时的较好权衡。体系大的话也可以用大一些的格点间距来节约耗时但会损失一些精度
0  //返回主菜单
13  //处理格点数据
11  //格点数据运算
13  //取绝对值。轨道波函数有正有负，这一步使之都变成正值
-1  //返回主菜单
12  //定量分子表面分析
1  //选择定义表面的方式
11  //用内存中的格点数据定义等值面
0.05  //等值面数值
6  //开始表面分析且不考虑被映射的函数

从屏幕上可以看到以下结果，即体积为14.2 Angstrom^3，顺带着轨道等值面的面积69.0 Angstrom^2也显示出来了。

 Volume:    95.86111 Bohr^3  (  14.20515 Angstrom^3)
 Estimated density according to mass and volume (M/V):   25.0417 g/cm^3
 Overall surface area:         246.50998 Bohr^2  (  69.02983 Angstrom^2)

感兴趣的话，还可以选择选项-3 Visualize the surface看一下当前等值面对应的图像，体积就是这个等值面里面包围的部分

对比一下，看看下面这个第18号分子轨道的体积如何

退回到主菜单，把下面这一串命令复制到Multiwfn窗口里即可，结果是7.3 Angstrom^3，可见只有HOMO的一半，和肉眼看上去的轨道体积差异相一致。
5
4
18
4
0.25
0
13
11
13
-1
12
1
11
0.05
6

3 批量计算

在《详谈Multiwfn的命令行方式运行和批量运行的方法》（//www.umsyar.com/612）里专门讲了怎么用Multiwfn批处理计算。这里提供一个脚本，可以便捷地用前文的方法批量计算特定序号范围的轨道的体积。//www.umsyar.com/attach/734/orbvol.sh是个Bash shell脚本，用来计算examples\excit\D-pi-A.fchk这个体系所有占据轨道（1到56号）的体积，内容如下所示。其中istart和iend是起始和终止的轨道序号，这俩值和被处理的波函数文件路径都可以根据实际需要修改。运行之前，Multiwfn可执行文件所在目录需要添加到PATH环境变量下使得能够通过Multiwfn命令启动之。

#!/bin/bash
istart=1
iend=56
for ((i=$istart;i<=$iend;i=i+1))
do
cat << EOF >> orbvol.txt
5
4
$i
4
0.25
0
13
11
13
-1
12
1
11
0.05
6
-1
-1
EOF
done
echo "q" >> orbvol.txt
echo "Running Multiwfn..."
Multiwfn examples/excit/D-pi-A.fchk < orbvol.txt >> result.txt
grep "Volume:" result.txt | nl -v$istart
rm ./orbvol.txt result.txt

运行后，屏幕上看到各个轨道的体积：
     1   Volume:     2.88332 Bohr^3  (   0.42726 Angstrom^3)
     2   Volume:     2.89404 Bohr^3  (   0.42885 Angstrom^3)
     3   Volume:     2.42980 Bohr^3  (   0.36006 Angstrom^3)
     4   Volume:     2.40622 Bohr^3  (   0.35656 Angstrom^3)
...略
    54   Volume:    95.10824 Bohr^3  (  14.09359 Angstrom^3)
    55   Volume:    95.86240 Bohr^3  (  14.20534 Angstrom^3)
    56   Volume:    95.86111 Bohr^3  (  14.20515 Angstrom^3)

将轨道体积相对于轨道序号作柱形图，如下所示。可见前16号轨道的体积非常小，这是因为它们都是内核轨道。

笔者之前写过《通过轨道离域指数(ODI)衡量轨道的空间离域程度》（//www.umsyar.com/525），里面用我提出的ODI指数考察了与本文相同的D-pi-A体系的各个轨道的离域程度，并且绘制了ODI与轨道序号之间的柱形图。将上图与ODI的柱形图相对比，可以看到两个图有明显的互补性，即ODI越低（体现轨道越离域、越能同时覆盖很多原子），轨道体积倾向于越大。如果你的目的就是想衡量轨道的离域程度，ODI相对更严格，毕竟不依赖于轨道等值面数值的选取（轨道等值面数值设得很大的话，轨道等值面就完全看不到了，体积为0）。而用轨道体积来说事的好处是可以直接与等值面图相对应，结合图像讨论非常清楚直观。

4 其它

本文的做法不限于计算考察分子轨道的体积，对任意轨道都可以考察，如定域化轨道、NTO轨道、双正交化轨道、NAdO轨道、AdNDP轨道等等，载入含有相应轨道的波函数文件即可，产生方法在下文说了。
Multiwfn的轨道定域化功能的使用以及与NBO、AdNDP分析的对比
//www.umsyar.com/380
用于非限制性开壳层波函数的双正交化方法的原理与应用
//www.umsyar.com/448
使用键级密度(BOD)和自然适应性轨道(NAdO)图形化研究化学键
//www.umsyar.com/535
使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键
//www.umsyar.com/138