文/Sobereva@北京科音  2024-Jul-31

0 前言

衡量化学体系的芳香性的方法非常多，见《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）。研究芳香性的文章数目甚巨，但大多数都是对分子、原子团簇这样的孤立体系研究的，主要在于很少有程序能够支持将衡量芳香性的方法用于周期性体系。如今Multiwfn已经可以将很多方法用于周期性体系了，包括多中心键级、AV1245、AVmin、PDI、FLU、FLU-pi、PLR、HOMA、Bird、ELF二分值、Shannon芳香性指数、芳香环的环临界点属性，等等，还可以通过绘制LOL-pi函数图像直观考察共轭情况。在本文中，将对其中大部分方法在Multiwfn中的操作进行演示。作为例子的体系是一个共价有机框架化合物（COF）的一层，此体系在《使用Multiwfn对周期性体系做键级分析和NAdO分析考察成键特征》（//www.umsyar.com/719）中已经作为例子分析过，结构如下所示，本文要通过不同方法对比绿色虚线里三个环的芳香性的大小。可见1号环是C3N3，原子序号为1,2,21,22,11,12。2号环和3号环相当于萘片段中的两个六元环，2号环序号为19,20,50,49,18,48，3号环序号为14,48,18,44,46,16。这里序号是按照成键关系排的。

笔者假定读者已经阅读过《衡量芳香性的方法以及在Multiwfn中的计算》充分了解了各种衡量芳香性方法的特征，也假定读者有了Multiwfn的基本使用常识，不了解者建议阅读《Multiwfn入门tips》（//www.umsyar.com/167）和《Multiwfn FAQ》（//www.umsyar.com/452）。Multiwfn可以从官网//www.umsyar.com/multiwfn免费下载，注意务必使用2024-Jul-1及以后更新的版本（注意看Multiwfn启动时的更新日期的提示），否则情况和本文所述不符。

本文涉及的多数芳香性分析方法都是基于波函数的，上面COF这个体系的由CP2K程序在PBE/DZVP-MOLOPT-SR-GTH级别计算得到的molden格式的波函数文件和//www.umsyar.com/719里用的是一致的，请从此文的链接里下载，是其中的COF-MOS-1_0.molden。此文件的产生方法在此文里也详细说了，不会用CP2K者必须仔细看。本文涉及到的诸如HOMA这样的芳香性指数是基于几何结构的，只要给Multiwfn提供的文件里包含结构信息就可以，如cif、pdb、xyz、gjf、mol2等常见格式都可以。如果被分析的区域是跨晶胞的，则输入文件必须能给Multiwfn提供晶胞信息，哪些格式能提供晶胞信息见《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）里的说明。

2  计算多中心键级

多中心键级是非常好、很严格的考察芳香性的指标。这一节计算一下当前体系中三个六元环的六中心键级，以此考察它们六中心共轭的强弱，这正比于它们的芳香性。多中心键级的概念参看《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）和《使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键》（//www.umsyar.com/138）。

启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
25  //离域性与芳香性分析
1  //多中心键级
1,2,21,22,11,12  //计算1号环。原子序号按原子连接关系输入
此时得到多中心键级结果为0.0399。类似地，再分别输入19,20,50,49,18,48 [回车]和14,48,18,44,46,16 [回车]计算2、3号环，多中心键级分别为0.0212和0.0281。

根据六中心键级大小可见此体系中的C3N3环（1号环）的芳香性显著强于六元碳环，而不与C3N3环连接的六元碳环（3号环）的芳香性比与C3N3环连接的六元碳环（2号环）更强。

3 计算AV1245和AVmin

AV1245的概念和计算方法在《使用Multiwfn计算AV1245指数研究大环的芳香性》（//www.umsyar.com/519）中介绍过，这里就不多说了。这里用AV1245算一下前述的三个环的芳香性。

启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
25  //离域性与芳香性分析
2  //AV1245
1,2,21,22,11,12  //计算1号环。原子序号按连接关系输入
看到以下结果
AV1245 times 1000 for the selected atoms is    7.35835761
AVmin times 1000 for the selected atoms is     7.358354 (   11   12    2   21)
即AV1245乘以1000和AVmin乘以1000都为7.358。

再输入19,20,50,49,18,48计算2号环，结果为
AV1245 times 1000 for the selected atoms is    4.46380673
AVmin times 1000 for the selected atoms is     3.441219 (   20   50   18   48)

再输入14,48,18,44,46,16计算3号环，结果为
AV1245 times 1000 for the selected atoms is    6.30980750
AVmin times 1000 for the selected atoms is     4.014836 (   46   16   48   18)

可见，无论是从AV1245还是AVmin上来看，芳香性都是1号环>3号环>2号环，这和多中心键级的结论完全一致。

4 计算PDI

这一节用PDI衡量芳香性。注意PDI只能用于六元环。Multiwfn中PDI是基于模糊空间定义的离域化指数算的，对周期性体系默认用的是Hirshfeld原子空间划分。

启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
15  //模糊空间分析
5  //PDI
[回车]  //用默认的格点间距计算原子重叠矩阵（AOM）。对较大体系可以用更大的比如0.35 Bohr格点间距明显节约这一步的计算时间，对结果影响甚微
1,2,21,22,11,12  //第一个环里的原子序号
结果如下

Delocalization index of     1(C )   --   22(N ):    0.104520
Delocalization index of     2(N )   --   11(C ):    0.104520
Delocalization index of    21(C )   --   12(N ):    0.104520
PDI value is    0.104520

再依次输入19,20,50,49,18,48 [回车]和14,48,18,44,46,16 [回车]分别计算第2、3个环的PDI，结果分别为0.083658和0.096507。PDI越大芳香性越强，可见PDI给出的芳香性大小的结论也是1号环>3号环>2号环。

5 计算FLU和FLU-pi

先计算FLU芳香性指数。启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
15  //模糊空间分析
6  //FLU
[回车]  //用默认的格点间距计算原子重叠矩阵（AOM）
现在从屏幕上可以看到做FLU计算用的参数。输入第一个环里的原子序号1,2,21,22,11,12，得到如下结果

        Atom pair         Contribution          DI
   1(C )  --    2(N ):       0.000800        1.474298
   2(N )  --   21(C ):       0.000312        1.508760
  21(C )  --   22(N ):       0.000800        1.474298
  22(N )  --   11(C ):       0.000312        1.508761
  11(C )  --   12(N ):       0.000800        1.474298
  12(N )  --    1(C ):       0.000312        1.508761
FLU value is    0.003335

再依次输入19,20,50,49,18,48 [回车]和14,48,18,44,46,16 [回车]分别计算第2、3个环的FLU，结果分别为0.017882和0.012199。由于FLU越小芳香性越强，因此FLU的芳香性大小的结论是1号环>3号环>2号环，和前述芳香性指标结论相同。

下面再来计算FLU-pi。算这个需要告诉Multiwfn当前体系的所有pi占据轨道序号。可以按照《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）肉眼一个个看来记录序号，但太麻烦。建议按照《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）介绍的做法自动指认。启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
100  //其它功能（Part 1）
22  //自动检测pi轨道
0  //当前轨道是离域的（分子/晶体轨道属于这类）
马上Multiwfn就显示了占据的pi轨道序号，为47,50,60,61,64,70-73,77,81,82,84-89,92-94。下面计算FLU-pi，接着输入
0  //对识别出的pi轨道什么都不做
0  //返回主菜单
15  //模糊原子空间分析
7  //FLU-pi
[回车]  //用默认的格点间距计算原子重叠矩阵（AOM）
47,50,60,61,64,70-73,77,81,82,84-89,92-94  //pi占据轨道的序号
之后输入第一个环里的原子序号1,2,21,22,11,12，得到如下结果

Average of DI-pi is    0.387267
        Atom pair         Contribution          DI
   1(C )  --    2(N ):       0.000191        0.375470
   2(N )  --   21(C ):       0.000191        0.399063
  21(C )  --   22(N ):       0.000191        0.375471
  22(N )  --   11(C ):       0.000191        0.399063
  11(C )  --   12(N ):       0.000191        0.375470
  12(N )  --    1(C ):       0.000191        0.399063
FLU-pi value is    0.001148

再依次输入19,20,50,49,18,48 [回车]和14,48,18,44,46,16 [回车]分别计算第2、3个环的FLU-pi，结果分别为0.028167和0.026396。由于FLU-pi越小芳香性越强，因此FLU-pi的芳香性大小的结论是1号环>3号环>2号环，和前述芳香性指标结论相同。

6 计算HOMA

HOMA是流行的基于环上的键长特征衡量芳香性的方法。给Multiwfn提供的输入文件里有结构信息就行了。由于COF-MOS-1_0.molden里也包含结构信息，所以此例还是用这个作为输入文件。启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
25  //芳香性分析
6  //HOMA和Bird芳香性指数
0  //开始计算HOMA（基于默认参数）
现在屏幕上显示了计算HOMA用的参数。输入第一个环里的原子序号1,2,21,22,11,12，得到如下结果

        Atom pair         Contribution  Bond length(Angstrom)
   1(C )  --    2(N ):      -0.003592        1.349180
   2(N )  --   21(C ):      -0.001191        1.342742
  21(C )  --   22(N ):      -0.003592        1.349181
  22(N )  --   11(C ):      -0.001191        1.342741
  11(C )  --   12(N ):      -0.003592        1.349181
  12(N )  --    1(C ):      -0.001191        1.342741
HOMA value is    0.985651

可见HOMA为0.985651，并且环上各个键的键长以及对HOMA的贡献量都给出了。再依次输入19,20,50,49,18,48 [回车]和14,48,18,44,46,16 [回车]分别计算第2、3个环的HOMA，结果分别为0.598160和0.781579。由于HOMA越接近1芳香性越强，因此HOMA的芳香性大小的结论是1号环>3号环>2号环，和前述芳香性指标结论相同。

Bird芳香性指数的计算方法和HOMA非常类似，只不过是在主功能25里的子功能6里选择2而非此例的0而已，这里就不演示了。

7 计算Shannon芳香性指数

Shannon芳香性指数是基于被考察的环上的各个键的键临界点的电子密度定义的，因此计算它之前必须先用Multiwfn搜索临界点。这方面的操作在《使用Multiwfn结合CP2K做周期性体系的atom-in-molecules (AIM)拓扑分析》（//www.umsyar.com/717）里有详细说明，这里就不对细节做具体介绍了。

启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
2  //拓扑分析
3  //用每一对原子连线的中点作为初猜点搜索临界点
0  //在图形窗口里查看结果

在图形界面里要求只显示(3,-1)临界点，也就是键临界点，并且要求把临界点序号显示出来，看到下图。可见第一个环上的临界点序号为72,76,67,59,55,61。

点击图形界面右上角的Return按钮，然后输入
20  //计算Shannon芳香性
72,76,67,59,55,61  //第1个环上的BCP序号
马上看到结果：

Electron density at CP 72:   0.3352748688  Local entropy:   0.2993061920
Electron density at CP 76:   0.3318277479  Local entropy:   0.2979424307
Electron density at CP 67:   0.3352750629  Local entropy:   0.2993062683
Electron density at CP 59:   0.3318276569  Local entropy:   0.2979423945
Electron density at CP 55:   0.3352752224  Local entropy:   0.2993063310
Electron density at CP 61:   0.3318279246  Local entropy:   0.2979425011
Total electron density:   2.0013084835
Total Shannon entropy:   1.7917461175
Expected maximum Shannon entropy:   1.7917594692

Shannon aromaticity index:    0.0000133518

即Shannon芳香性指数为0.0000133518，每个键临界点的电子密度和局部熵也都给出了，这些是计算Shannon芳香性指数的中间量。

类似地也输入第2、3个环上的键临界点的序号，分别为62,64,54,45,43,51和64,71,82,84,75,65，结果分别为0.0013759152和0.0010309007。由于Shannon芳香性指数越小芳香性越强，因此芳香性大小的结论是1号环>3号环>2号环，和前述芳香性指标结论相同。

8 基于环临界点属性判断芳香性

在某个环中央区域的环临界点位置上，若垂直于环的方向上电子密度曲率越负，说明这个环的芳香性越强。这种判断芳香性的方法需要利用AIM拓扑分析得到的环临界点的属性。为了用这种分析，和上一节一样先用主功能2做拓扑分析，然后在选项0里只显示(3,+1)临界点，即环临界点，看到的图如下。可见第1、2、3号环的环临界点序号分别为66、53、73。

关闭图形窗口后输入
21  //计算顺着某个方向的电子密度梯度和曲率
66  //第1个环的环临界点序号
1  //指定某个方向
0,0,1  //当前体系平行于XY平面，因此计算Z方向的电子密度的梯度和曲率
结果如下

Electron density is                            0.0287379670 a.u.
Electron density gradient is                   0.0000000000 a.u.
Electron density curvature is                 -0.0261250099 a.u.

可见，在第一个环的环临界点位置上垂直于这个环的电子密度曲率为-0.0261250099 a.u.。类似地再使用这个功能考察第2、3个环分别对应的第53、73号环临界点垂直于环平面的电子密度曲率，结果分别为-0.0154770974和-0.0159388407 a.u.。由于数值越负芳香性越强，因此芳香性大小的结论是1号环>3号环>2号环，和前述芳香性指标结论相同。

9 ELF-pi二分点考察芳香性

这一节通过各个环上ELF-pi函数的等值面首次发生二分的ELF-pi数值（ELF-pi二分值）来衡量芳香性的大小。这个值既可以通过反复微调等值面数值来获得，也可以通过对ELF-pi做拓扑分析来获得，前者更直观，后者更精确，相关信息见《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）。这一节演示前者的做法。

首先需要产生ELF-pi的格点数据。启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
100  //其它功能（Part 1）
22  //检测pi轨道
0  //轨道都是离域形式
2  //设置其它轨道占据数为0
0  //返回主菜单
5  //计算格点数据
9  //ELF
9  //利用晶胞平移矢量定义格点信息
[回车]  //坐标原点为(0,0,0)
[回车]  //盒子三个方向尺寸和相应晶胞边长一致
0.15  //用较精细的0.15 Bohr格点间距，使得通过观看等值面获得二分值能尽可能准确
2  //将格点数据导出为cube文件

现在当前目录下就有了ELF.cub，对应ELF-pi的格点数据。用VMD显示其等值面（不会操作的话用《在VMD里将cube文件瞬间绘制成效果极佳的等值面图的方法》//www.umsyar.com/483里提供的脚本）。选Display - Orthographic用正交视角，在Graphics - Representation里将等值面数值（isovalue）由小到大调节，会发现在等值面数值为0.761的时候正好1号环上等值面首次发生二分，位置如下图箭头所示。因此1号环的二分点数值为0.761。

由上图可见，在1号环ELF-pi等值面二分的时候，2、3号环早就发生了二分、等值面间间隔已经很大了，这说明1号环的芳香性远强于2、3号环。

再将等值面数值调节到0.623，会看到2、3号环上首次出现了ELF-pi的二分，如下图红色箭头所示。这个键是被两个六元环共享的，怎么区分哪个芳香性更强？这需要再看其它的键。在下图紫色箭头所示的2号环上的位置，只要等值面数值再增加一点点就会二分，而3号环就没有这个情况，说明3号环的芳香性比2号环更强、环上的电子共轭遇到的瓶颈更少。因此ELF-pi的分析结论和其它方法完全一致。

10 LOL-pi图形分析

作为前述分析的扩展和补充，这一节绘制LOL-pi图直观展示一下单层COF体系的pi共轭情况，这能十分清楚地让大家理解三个环上电子离域特征的差异。这种分析在前述的《在Multiwfn中单独考察pi电子结构特征》文中有大量介绍，在笔者的Theor. Chem. Acc, 139, 25 (2020)中也有很多例子，欢迎阅读和引用。

启动Multiwfn并载入COF-MOS-1_0.molden后，依次输入
100  //其它功能（Part 1）
22  //检测pi轨道
0  //轨道都是离域形式
2  //设置其它轨道占据数为0
0  //返回
4  //绘制平面图
10  //LOL
1  //填色图
[回车]  //用默认的格点数
1  //XY平面
4.25  //当前体系中的原子都在Z=3.25位置，绘制平面上方1 Bohr处的图像显然应该设Z=4.25 Bohr

利用后处理菜单对作图效果的一番调节，得到下图。如果不会调的话，仔细把Multiwfn手册4.4节的所有例子都仔细领会了，结合后处理菜单中提示得很清楚的选项名字去理解，自然就明白了。当前用的色彩刻度是0到0.65，超过上限的区域显示为白色。

由上图可见，在1号环上pi共轭很显著，LOL-pi在环上分布得比较均匀。而如图上我标注的箭头所示，2号环上有三个键的pi共轭显著低于其它地方，因此这个环上的六中心pi共轭作用明显受到了很大遏制。而3号环上的pi共轭整体相对好点，但也不是很均匀，而是在不同键上有大有小，所以芳香性只是介于1和2号环之间。

11 总结

Multiwfn支持极为丰富的考察化学体系芳香性的方法，本文对其中已支持周期性体系的大部分方法的操作过程进行了简明扼要的演示。从结果可见尽管这些方法思想差异很大，但对于当前研究的这个体系中的三个六元环，它们给出的芳香性顺序完全一致，互相印证了彼此的可靠性。当然对于一些特殊情况，由于一些方法原理上的局限性、稳健性和普适性的不足也有可能导致结果存在差异。本文示例的只是一个很标准、简单的体系，希望大家能充分举一反三，将Multiwfn应用于广泛体系的芳香性的研究中。

使用Multiwfn做任何分析在发表文章时都请务必记得按照程序启动时的提示恰当引用Multiwfn原文，给别人代算时也必须明确告知对方这一点。

文/Sobereva@北京科音  2024-Jul-12

0 前言

键级是考察化学键特征的一类重要方法，非常强大的波函数分析程序Multiwfn可以方便地计算很多类型的键级，在《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）的键级部分有具体的说明，Multiwfn手册4.9节有一些实际计算例子，在“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/workshop/WFN_content.html）里做了特别全面详细的讲解和演示。另外，Multiwfn还支持NAdO方法，可以将原子间的模糊键级以轨道形式图形化展现，对于了解成键本质极其有用，介绍和实例看《使用键级密度(BOD)和自然适应性轨道(NAdO)图形化研究化学键》（//www.umsyar.com/535）。

以上博文都是举的孤立体系的例子，而如今Multiwfn已经将上述方法扩展到了对周期性体系的分析上。将Multiwfn与CP2K做周期性DFT计算产生的波函数相结合，可以很容易地用上述方法考察周期性体系的成键情况，在本文将通过一个层状的共价有机框架化合物（COF）体系，以及一个Pd表面吸附苯的体系，对具体做法进行演示，希望读者举一反三将本文介绍的方法运用到实际研究中。使用本文的方法计算的结果若用于发表，记得需要按照Multiwfn启动时的提示恰当引用Multiwfn的原文。

读者务必使用2024-Jun-24及以后更新的Multiwfn版本，否则情况与本文所述不符。Multiwfn可以在主页//www.umsyar.com/multiwfn免费下载。不了解Multiwfn者建议阅读《Multiwfn入门tips》（//www.umsyar.com/167）和《Multiwfn FAQ》（//www.umsyar.com/452）。本文要用到非常流行、高效且免费的第一性原理程序CP2K，笔者假定读者已具备了相关常识，推荐不熟悉者通过“北京科音CP2K第一性原理计算培训班”（http://www.keinsci.com/workshop/KFP_content.html）系统性学习。本文要使用Multiwfn创建CP2K输入文件，这在《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）里有简要介绍。本文用的CP2K是2024.1版。

下面的例子涉及到的所有文件都可以在//www.umsyar.com/attach/719/file.rar里得到。

1 单层COF体系

这部分示例的COF体系的cif文件是本文文件包里的COF_16371N2.cif，用GaussView显示的结构如下，可见晶胞里有两层。

1.1 准备输入文件

此体系两层之间的相互作用主要是pi-pi堆积，对电子结构影响微乎其微，因此只计算一层就够了，可以节约时间，尤其是能显著节约NAdO分析过程中计算原子重叠矩阵（AOM）的时间。于是在GaussView里删除一层，然后保存为COF.gjf。由于当前体系不是导体且晶胞边长约15埃已经不小了，所以可以不用扩胞就直接做gamma点计算。注意CP2K导出molden文件只支持gamma点的情况，这在《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）里已明确说了。

启动Multiwfn，然后输入
COF.gjf
cp2k  //进入产生CP2K输入文件的界面
[回车]  //输入文件名用默认的COF.inp
-2  //要求产生molden文件
-7  //设置周期性
XY  //仅在平行于COF的方向考虑周期性，而垂直于COF的方向用非周期性
-9  //其它设置
16  //要求将Kohn-Sham矩阵导出到.csr文件里。此文件在Multiwfn计算NAdO轨道能量时会用到
0  //返回
0  //产生输入文件

现在当前目录下就有了COF.inp，对应PBE/DZVP-MOLOPT-SR-GTH级别的单点计算。用CP2K运行之，当前目录下会产生记录波函数信息的COF-MOS-1_0.molden和记录KS矩阵的COF-KS_SPIN_1-1_0.csr。按照《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》所述，在molden文件开头插入以下内容定义晶胞信息和当前赝势基组描述的各元素的价电子数
 [Cell]
 14.94670000     0.00000000     0.00000000
 -7.47335000    12.94422190     0.00000000
  0.00000000     0.00000000     8.00000000
 [Nval]
 C 4
 N 5
 H 1

之后就可以基于这个molden格式的波函数文件做各种周期性体系的分析了。

1.2 计算Mayer键级

首先计算非常常用的Mayer键级，它衡量原子间等效的共享电子的对数。启动Multiwfn，载入COF-MOS-1_0.molden，然后输入
9  //键级分析
1  //Mayer键级
Multiwfn首先计算重叠矩阵，然后立刻输出了Mayer键级，只有大于阈值的（由settings.ini里的bndordthres参数控制）才直接显示了出来。要想看所有原子间的Mayer键级可以再选择y导出键级矩阵然后查看相应的矩阵元。为了便于对照，下图显示了当前体系各个原子的序号，本文主要关注其中绿色虚线圈住的那些原子。

Multiwfn输出的一些有代表性的键的Mayer键级值如下
11(C )   19(C )    1.12454138
11(C )   22(N )    1.48268711
19(C )   48(C )    1.27625346
20(C )   50(C )    1.39369508
18(C )   48(C )    1.24722702
可见这些键的Mayer键级都不同程度地明显大于1.0，因此可以推测它们都不仅仅形成了sigma键，还有一定pi成份。在后文会计算pi键级和使用NAdO分析更进一步确认这一点。上面列出的C-C键键级可以体现出体系中哪种C-C键的强度更强。可见连接萘单元和C3N3单元的C11-C19的键是相对最弱的。

1.3 计算模糊键级（fuzzy bond order）

Mayer键级怕弥散函数，而模糊键级则没有这个问题，用于任意基组都可以。Mayer键级用于MOLOPT系列基组都是没问题的，不过这里作为例子，也计算一下模糊键级。由于计算模糊键级需要先构造AOM，对较大的体系很耗时，所以对于Mayer键级适用的情况建议优先用Mayer键级。注意对孤立体系Multiwfn计算模糊键级默认是基于Becke划分的，而对于周期性体系计算模糊键级默认是基于Hirshfeld划分的，因此默认设置下周期性体系的计算结果和孤立体系的没有严格的可比性。

在键级分析主功能的菜单里选择7，然后Multiwfn就会开始计算AOM，之后立刻给出模糊键级的计算结果，前述的那些键的模糊键级如下
11(C )   19(C )    1.05202861
11(C )   22(N )    1.50866214
19(C )   48(C )    1.22313834
20(C )   50(C )    1.38494183
18(C )   48(C )    1.17984661
虽然模糊键级与Mayer键级在定量上有一些差别，但不同的键的键级的大小顺序是完全一致的，没有结论上的差异。

1.4 计算pi电子贡献的Mayer键级

我在《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）中专门介绍了pi键级的计算方法，也即pi占据轨道对Mayer键级的贡献，没读过的话务必先阅读。这一节对单层COF这个周期性体系也做这个计算，看看各个键的pi作用程度的差异。由于当前体系是精确平行于XY平面的，所以让Multiwfn自动指认pi轨道前不需要先做轨道定域化。

启动Multiwfn并载入COF-MOS-1_0.molden后，依次输入
100  //其它功能（Part 1）
22  //检测pi轨道
0  //当前的轨道都是离域形式（如正则分子轨道）
2  //设置其它轨道占据数为0
0  //返回
9  //键级分析
1  //Mayer键级
前面提到的那些键的键级计算结果如下，这对应的就是pi键级。可见C3N3六元环中的C-N键的pi共享电子作用很强，比这里列出的C-C键的还要更显著。
11(C )   19(C )    0.14583409
11(C )   22(N )    0.41775145
19(C )   48(C )    0.28202057
20(C )   50(C )    0.38270736
18(C )   48(C )    0.27414145

使用此做法发表文章时，除了引用Multiwfn原文外还建议同时引用介绍Multiwfn做pi电子结构分析的文章Theor. Chem. Acc., 139, 25 (2020)。

1.5 NAdO分析

利用前述的《使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键》中介绍的NAdO分析，可以将sigma和pi作用分别以轨道形式呈现出来，非常有价值。例如在我的研究18氮环的ChemPhysChem (2024) https://doi.org/10.1002/cphc.202400377文中就用了NAdO分析非常清楚、严格地展示了较短的N-N键的sigma+pi作用特征。本节我们用NAdO分析展示一下C11-C19的sigma和pi轨道相互作用情况。

启动Multiwfn并载入COF-MOS-1_0.molden后，依次输入
15  //模糊空间分析
3  //计算AOM。对周期性体系默认是基于Hirshfeld原子空间计算的
[回车]  //用默认的0.2 Bohr格点间距。此设置对当前体系耗时不高，但对于明显更大体系为了节约耗时可以适当用更大的格点间距如0.35 Bohr，但会多多少少牺牲精度

计算完成后当前目录下出现了记录AOM的AOM.txt。在Multiwfn里接着输入
0  //返回
200  //其它功能（Part 2）
20  //BOD和NAdO分析
-1  //要求计算NAdO的能量
2  //从外部文件读取Fock矩阵
COF-KS_SPIN_1-1_0.csr  //输入此文件的实际路径
1  //基于AOM计算原子间的相互作用
[回车]  //读取当前目录下的AOM.txt
11,19  //分析11和19号原子间的相互作用
马上NAdO分析就做完了，屏幕上显示以下内容。这说明所有NAdO的本征值加和为1.052，和上一节得到的模糊键级一致。并且当前只有前两个NAdO的本征值明显大于0，它们是最值得考察的，其它的基本可以忽略

 Eigenvalues of NAdOs: (sum=   1.05204 )
   0.78955   0.19252   0.06784   0.02133   0.00807   0.00651   0.00127
   0.00096   0.00031   0.00011   0.00010   0.00005   0.00004   0.00000
...略

所有的NAdO轨道现在都已经导出到了当前目录下的NAdOs.mwfn中。现在输入y让Multiwfn载入此文件，之后退回到Multiwfn主菜单，进入主功能0查看NAdO轨道。不熟悉Multiwfn看轨道功能的话参考《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）。在图形界面里分别选择1号和2号轨道，可分别看到下面的图像，同时在文本窗口显示的NAdO本征值和轨道能量也标上了。这里等值面数值用的是0.08。

由上可见，C11-C19同时具有sigma和pi作用特征，前者远比后者显著得多。而且由于C19处在pi共轭区域，所以NAdO方法产生的主要描述C11-C19的pi作用的轨道并不完全定域在C11和C19上，也同时一定程度离域到了周围原子上。通过NAdO轨道能量还可以看到，C11和C19之间共享的sigma电子的能量明显比pi电子要低，这十分符合一般化学常识。

2 Pd(100)表面吸附苯

本节通过Pd(100)晶面吸附苯的体系，演示一下苯与Pd(100)基底这两个片段间Mayer键级的计算，以及对于像这样很大的体系怎么尽可能节约整个NAdO分析的时间。此体系的CP2K做结构优化的输入和输出文件，以及优化任务跑完后产生的restart文件都在本文文件包里提供了。优化完的结构如下，可见吸附作用很强，以至于Pd-C的作用都令苯产生弯曲了。

由于此例子的molden文件、FOM.txt和NAdOs.mwfn文件太大，本文的文件包里就没提供了。

2.1 准备输入文件

用Multiwfn载入本文文件包里的Pd+ben_opt-1.restart读取优化完的坐标和晶胞信息，然后输入以下命令创建一个做单点且产生molden文件的任务
cp2k  //产生CP2K输入文件
Pd+ben_SP.inp  //产生的输入文件名
-7  //设置周期性
XY
-2  //要求产生molden文件
6  //由于Pd基底是导体，开启smearing
0  //产生输入文件

为了帮助SCF收敛，把Pd+ben_SP.inp里的ALPHA设为0.15，NBROYDEN设为12。用CP2K计算Pd+ben_SP.inp，得到Pd+ben_SP-MOS-1_0.molden。然后在里面开头部分加入以下内容。

 [Cell]
 23.33880000     0.00000000     0.00000000
 0.00000000    23.33880000     0.00000000
 0.00000000     0.00000000    24.79750184
 [Nval]
 Pd 18
 C 4
 H 1

2.2 片段间Mayer键级

这一节将苯和Pd(100)基底分别定义为一个片段计算它们之间的Mayer键级。启动Multiwfn，载入Pd+ben_SP-MOS-1_0.molden，然后输入
9  //键级分析
-1  //定义片段
289-300  //苯的部分
1-288  //Pd(100)基底部分
1  //计算键级

Multiwfn开始产生重叠矩阵，之后给出了原子间键级，以下是其中一部分，可见苯中的C和有的Pd之间的Mayer键级较大，即共价作用很显著。

...略
 # 1412:       246(Pd)  270(Pd)    0.21887832
 # 1413:       246(Pd)  271(Pd)    0.21771061
 # 1414:       246(Pd)  289(C )    0.08497491
 # 1415:       246(Pd)  293(C )    0.09391700
 # 1416:       246(Pd)  294(C )    0.55236636
 # 1417:       247(Pd)  271(Pd)    0.23611928
 # 1418:       247(Pd)  272(Pd)    0.23841126
...略

老有人问怎么判断固体表面对分子是物理吸附还是化学吸附，判断方法很多，如结合能的数量级、原子间距离和原子半径的关系、电子密度差、IGMH（//www.umsyar.com/621）等等，而用Mayer或模糊键级是最能说明问题的。以上Mayer键级体现的C和Pd之间的明显的共价作用意味着二者之间形成了明显的化学键作用，显然Pd(100)对苯是化学吸附。如果键级只有比如零点零几的程度，那一般可以说是物理吸附，但也不排除形成典型离子键的可能，这通过原子电荷可以判断，计算方式可参考《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）。

最后Multiwfn还给出了片段间键级：
The bond order between fragment 1 and 2:    4.309784
即苯分子和Pd(100)之间等效共享电子对数多达4.3，显然太共价了，如果这都不叫化学吸附...

2.3 NAdO分析

这一节我们将Pd(100)与苯之间显著的共价作用通过NAdO方法以轨道形式直观展现。本节的做法和1.5节示例的常规的NAdO分析流程有两个明显的不同：
(1)对当前这样原子数又多、轨道数又多的大体系计算所有原子的AOM是相当耗时间的事情，而且导出、载入AOM.txt的耗时巨高而且巨占硬盘（正比于占据轨道数的平方乘以原子数）。为了极大程度节约时间，可以使用Multiwfn提供的一个特殊的策略，也就是在主功能15里用选项33计算感兴趣的两个片段间的片段重叠矩阵（FOM）并导出为FOM.txt，然后对两个片段之间做NAdO分析时直接从中读取要用的FOM，这样不仅省得导出和载入巨大的AOM.txt，而且产生FOM的过程中不需要计算片段里没有涉及的原子的AOM，当片段里的原子数占整个体系的原子数比例较少的时候这可以大幅节约时间。对于当前的例子，苯可以作为片段1，与苯相互作用最密切的那部分Pd原子适合作为片段2，这比起把所有Pd都定义为片段2能节约巨量耗时。但用户很难判断哪些Pd原子与苯作用密切，此时可以利用Multiwfn特意提供的一种定义片段2的方法，即片段1以外的原子中，与片段1任何一个原子之间Mayer键级大于指定阈值的原子都被定义为片段2。利用这个做法可以使得与苯的轨道相互作用最显著的Pd原子被定义为片段2，其数目远少于总的Pd数目。
(2)在计算FOM和做NAdO分析之前，需要先用主功能6的子功能38令电子根据轨道能量由低到高重新排列，使得轨道占据数都为整数。这是因为当前CP2K计算用了smearing，导致费米能级附近的轨道处于分数占据状态，这不是单行列式波函数应具有的特征，而Multiwfn的NAdO分析目前只支持单行列式波函数。重新排列后，当前的波函数就成了标准的单行列式波函数形式，从而可以做NAdO分析。另一方面，这么重排占据数后在计算AOM和FOM时可以只对占据轨道之间计算，也只有这部分是NAdO分析所实际需要的，只计算占据轨道之间的AOM/FOM的耗时、对内存的占用远远低于考虑所有轨道的情况，

启动Multiwfn，载入Pd+ben_SP-MOS-1_0.molden，然后输入
6  //修改波函数
38  //按照轨道能量由低到高重排占据数
-1  //返回
15  //模糊空间分析
33  //计算FOM
3  //定义两个片段并计算FOM，不属于片段1的原子若与片段1的任意原子的Mayer键级大于指定值则被定义为片段2
289-300  //苯的部分
0.001  //Mayer键级的阈值

之后看到以下信息，告诉了你哪些原子被作为了片段2，并且可见片段1和2之间的Mayer键级为4.114，和上一节看到的苯与整个Pd(100)表面的片段间Mayer键级4.310相差很小，这说明当前自动定义的片段2是合理的。如果你把片段2的原子在Multiwfn主功能0里用菜单栏的Other settings - Set atom highlighting功能高亮显示，会看到它们都是离苯不远的原子。

 Atoms in fragment 2:
 76,77,80-83,86,87,97,101-103,106-108,112,195-197,200-202,205-207,220,221,224-22
 6,230,231,241,245-247,250-252,256,265,266,269-272,274-277,280,281
 Mayer bond order between fragments 1 and 2:   4.114

之后输入0.35，Multiwfn就开始用0.35 Bohr间距的立方格点计算AOM，然后构造出FOM并导出到当前目录下的FOM.txt。在双路7R32机子上96核并行计算，整个过程花了26分钟（耗时与格点间距的三次方呈反比，耗时太高的话可以用比如0.5 Bohr格点间距，误差也还可以接受，此时只需要不到10分钟）。然后输入
0  //返回
200  //其它功能（Part 2）
20  //BOD/NAdO分析
4  //直接从当前目录下的FOM.txt中读取NAdO分析要用的片段1和2的FOM
之后开始了NAdO的计算，结果如下

 Eigenvalues of NAdOs: (sum=   6.95512 )
   0.79476   0.77393   0.71745   0.63686   0.62698   0.35965   0.22325
   0.16243   0.14986   0.14655   0.14545   0.13818   0.12618   0.12058
   0.11388   0.11078   0.10617   0.09335   0.09007   0.08979   0.08423
   0.08224   0.07942   0.07434   0.07391   0.07171   0.06475   0.05751
   0.05507   0.05478   0.05246   0.04964   0.04669   0.04030   0.03771
...略

当前显示的所有NAdO本征的加和6.955对应于苯和Pd(100)之间在Hirshfeld原子空间划分下算的模糊键级，和前面看到的片段间Mayer键级4.31有一定差异，这很正常，毕竟对原子空间的定义截然不同。由以上数据可见有不少NAdO轨道对模糊键级的贡献都显著，尤其是前六个。现在输入y载入新产生的NAdOs.mwfn，在主功能0里观看NAdO轨道，其中本征值最大的6个NAdO轨道如下所示，等值面数值用的0.03

从上图中可以看到这些轨道在苯与Pd(100)之间都是相位相同方式叠加的，必然都是起到明显成键作用的，所以本征值都不小。从轨道图形上可以看出这些轨道来自于苯上C原子的垂直于苯环的p轨道与Pd的原子轨道混合，而且有的图里直接就能清楚看出Pd用的是d原子轨道。例如NAdO(2)中和苯接触的Pd明显用的是dz2轨道，从Pd区域的轨道等值面形状就能看出这一点。感兴趣的话可以进一步按照《谈谈轨道成份的计算方法》（//www.umsyar.com/131）介绍的用SCPA方法做一下轨道成分分析。在Multiwfn主菜单里输入
8  //轨道成分分析
3  //SCPA方法
2  //2号NAdO轨道

Multiwfn马上就输出了轨道成份，下面列出的是Multiwfn返回的各个原子上各个角动量基函数产生的贡献。明显可以看到Pd主要都是用D基函数，而C用的都是P基函数。因此对NAdO(2)的轨道成份的分析体现了Pd的d原子轨道与C的p原子轨道的显著混合对Pd吸附苯有关键性贡献。可见以这种方式讨论成键特征能分析得巨清楚透彻！

Composition of each shell
Shell  1370 Type: D    in atom  196(Pd) :     0.50345 %
Shell  1405 Type: D    in atom  201(Pd) :     0.88530 %
Shell  1440 Type: D    in atom  206(Pd) :     0.50251 %
Shell  1538 Type: D    in atom  220(Pd) :     0.80192 %
Shell  1545 Type: D    in atom  221(Pd) :     0.76004 %
Shell  1570 Type: S    in atom  225(Pd) :     0.53572 %
Shell  1573 Type: D    in atom  225(Pd) :    10.40186 %
Shell  1577 Type: S    in atom  226(Pd) :     0.56745 %
Shell  1580 Type: D    in atom  226(Pd) :    10.21554 %
Shell  1608 Type: D    in atom  230(Pd) :     0.80255 %
Shell  1615 Type: D    in atom  231(Pd) :     0.76058 %
Shell  1713 Type: D    in atom  245(Pd) :     0.53415 %
Shell  1717 Type: S    in atom  246(Pd) :     0.96339 %
Shell  1720 Type: D    in atom  246(Pd) :    21.79157 %
Shell  1748 Type: D    in atom  250(Pd) :     0.53493 %
Shell  1752 Type: S    in atom  251(Pd) :     0.96056 %
Shell  1755 Type: D    in atom  251(Pd) :    21.80729 %
Shell  2019 Type: P    in atom  289(C ) :     1.48815 %
Shell  2024 Type: P    in atom  290(C ) :     1.47584 %
Shell  2029 Type: P    in atom  291(C ) :     6.18539 %
Shell  2034 Type: P    in atom  292(C ) :     1.44475 %
Shell  2039 Type: P    in atom  293(C ) :     1.41651 %
Shell  2044 Type: P    in atom  294(C ) :     6.20361 %

感兴趣的读者还可以用CDA方法从苯的分子轨道与Pd表面的晶体轨道的混合角度来考察Pd对苯吸附造成的电子转移和轨道相互作用的细节。详见《使用Multiwfn结合CP2K对周期性体系做电荷分解分析（CDA）》（//www.umsyar.com/716）。

3 总结

本文通过COF二维层状体系和Pd(100)表面吸附苯作为例子，演示了如何用Multiwfn程序计算原子间和片段间的Mayer键级和模糊键级，还演示了如何计算pi键级单独考察pi电子的贡献，以及讲解了如何做NAdO分析以轨道形式清楚直观地考察原子间或片段间的共价作用的内在特征。可见使用Multiwfn可以把成键特征情况展现得超级清楚透彻。本文的做法对其它情况的周期性体系，诸如三维原子晶体、过渡态结构等等，也都是完全适用的。更多细节和相关知识请参看本文提到的相关博文。

Multiwfn对周期性体系能算的键级不止本文示例的这些。对周期性体系Multiwfn还能算Mulliken键级和Wiberg键级，还能把Mulliken和Mayer键级分解成不同轨道的贡献，后者称为轨道占据数扰动的Mayer键级。具体介绍见《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471），做法见Multiwfn手册4.8节的相关例子，本文就不特意演示了。除了本文介绍的这些外Multiwfn还有很多其它考察成键的方法对周期性体系都可以用，比如AIM拓扑分析，见《使用Multiwfn结合CP2K做周期性体系的atom-in-molecules (AIM)拓扑分析》（//www.umsyar.com/717）。

使用Multiwfn计算的结果若用于发表，记得需要按照Multiwfn启动时的提示恰当引用Multiwfn的原文。给别人代算的话也必须明确告知对方这一点。

文/Sobereva@北京科音   2024-Jul-3

0 前言

Bader提出的非常知名的atom-in-molecules (AIM)理论在分析孤立体系和周期性体系的电子结构方面有非常广泛的应用，相关资料汇总见《AIM学习资料和重要文献合集》（http://bbs.keinsci.com/thread-362-1-1.html）。AIM分析框架中最常见的分析手段之一是电子密度的拓扑分析，也常被称为AIM拓扑分析，例如根据这种拓扑分析得到的键临界点（BCP）的位置可以用来讨论成键和弱相互作用特征，参考《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）和《Multiwfn支持的弱相互作用的分析方法概览》（//www.umsyar.com/252）里的相关介绍。非常强大的波函数分析程序Multiwfn（//www.umsyar.com/multiwfn）很早以前就已经支持了结合Gaussian、ORCA等 程序产生的波函数文件对孤立体系做AIM拓扑分析，见《使用Multiwfn做拓扑分析以及计算孤对电子角度》（//www.umsyar.com/108）和Multiwfn手册4.2节的大量例子。如今Multiwfn也已支持了对周期性体系做AIM拓扑分析，本文的目的是对此进行演示，以使得读者可以轻松地举一反三充分将这种手段应用于实际研究中。使用Multiwfn做AIM拓扑分析在发表文章时记得需要按照Multiwfn启动时的提示恰当引用Multiwfn的原文。

读者务必使用2024-Jun-16及以后更新的Multiwfn版本，否则情况与本文所述不符。Multiwfn可以在主页//www.umsyar.com/multiwfn免费下载。不了解Multiwfn者建议阅读《Multiwfn入门tips》（//www.umsyar.com/167）和《Multiwfn FAQ》（//www.umsyar.com/452）。本文要用到非常流行、高效且免费的第一性原理程序CP2K，不熟悉者推荐通过“北京科音CP2K第一性原理计算培训班”（http://www.keinsci.com/workshop/KFP_content.html）系统性学习。本文要使用Multiwfn创建CP2K输入文件，这在《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）里有简要介绍。本文用的CP2K是2024.1版。

下面的例子涉及到的所有文件都可以在//www.umsyar.com/attach/717/file.rar里得到。

1 实例1：NaCl晶体

这一节对NaCl做AIM拓扑分析，要使用CP2K产生的molden文件作为Multiwfn的输入文件。molden文件的产生方法在《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）里有详细说明。由于考虑k点时CP2K没法产生molden文件，因此必须使用足够大的超胞模型并只考虑gamma点。本文文件包里NaCl.cif是NaCl单胞的晶体结构文件，基于它我们创建CP2K算单点任务的超胞的输入文件，同时得到molden文件。

用Multiwfn载入NaCl.cif，然后输入
cp2k  //创建CP2K输入文件的界面
[回车]  //产生的输入文件名用默认的NaCl.inp
-2  //要求产生molden文件
-11  //几何操作界面
19  //构造超胞
2  //第1个方向平移复制为原先的2倍
2  //第2个方向平移复制为原先的2倍
2  //第3个方向平移复制为原先的2倍
-10  //返回
0  //产生输入文件

现在当前目录下就有了NaCl.inp，对应的是PBE/DZVP-MOLOPT-SR-GTH计算级别，对此体系是合适的。当前晶胞边长为11.28埃，对于绝缘体来说只考虑gamma点基本够用。用CP2K运行NaCl.inp，得到NaCl-MOS-1_0.molden。用文本编辑器打开NaCl.inp，将以下内容插入文件开头来定义晶胞信息从而使得Multiwfn在分析时能考虑周期性，并且定义Na和Cl的价电子数分别为9和7（不了解为什么这么设置的人看前述的《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》）。
 [Cell]
11.28112000     0.00000000     0.00000000
 0.00000000    11.28112000     0.00000000
 0.00000000     0.00000000    11.28112000
 [Nval]
 Na 9
 Cl 7

下面开始做AIM拓扑分析。启动Multiwfn，载入NaCl-MOS-1_0.molden，然后输入
2  //拓扑分析（默认是对电子密度做拓扑分析）
2  //以每个原子核位置为初猜搜索临界点

当前从屏幕上会看到找到了64个(3,-3)临界点，分别对应64个原子核位置的电子密度极大点。之所以虽然当前用的是赝势基组，在原子核位置还能找到电子密度极大点，是因为Multiwfn自动对这些原子指认了合适的EDF信息来描述内核电子，详见《在赝势下做波函数分析的一些说明》（//www.umsyar.com/156）。

然后再输入
3  //以每两个原子间连线中点作为搜索临界点的初猜位置
4  //以每三个原子的中心位置作为搜索临界点的初猜位置
0  //观看结果

此时文本窗口显示了各种临界点的数目，如下所示，并且提示Poincare-Hopf关系已经满足了，因此大概率所有临界点都已经找全了。
 The number of critical points of each type:
 (3,-3):    64,   (3,-1):   384,   (3,+1):   384,   (3,+3):    64
 Poincare-Hopf verification:    64  -   384  +   384  -    64  =     0
 Fine, Poincare-Hopf relationship is satisfied, all CPs may have been found

在Multiwfn的图形窗口可以看到下图，棕红色大球是Na，绿色大球是Cl，橙色、黄色、绿色小球分别是电子密度的(3,-1)、(3,+1)和(3,+3)临界点，而(3,-3)临界点都被原子球挡住了。在Multiwfn手册3.14节对Multiwfn的拓扑分析功能有完整、非常详细的说明，建议阅读。

下面再把AIM拓扑分析里常涉及的键径产生出来。在图形界面右上角点Return关闭之，然后选8产生(3,-3)和(3,-1)之间的拓扑路径，也即键径。在笔者的i9-13980HX CPU上24核并行很快就产生完了。之后再选0进入图形界面观看，看到下图，可见所有键径（棕色细线）都产生了出来，既连接相邻的Na和Cl，也连接每一对临近的两个Cl。

默认情况下Multiwfn把处于晶胞边缘的所有镜像原子和临界点也显示了出来，如果不想显示它们的话（也即只显示晶胞中唯一的那些），在图形界面顶端依次选择Other settings - Toggle showing all boundary atoms和Toggle showing all boundary CPs and paths。在文本窗口也会看到当前设置的状态。

下面来考察一下Na和Cl之间的BCP的属性。用选项0的图形界面右边的选项把其它临界点关闭显示而只显示(3,-1)临界点，选择CP labels复选框要求显示临界点序号，关闭键径显示并要求显示原子，然后就能看到下图，这是当前图像的一个边缘区域（点击save picture按钮可在当前目录下得到高分辨率图像文件，下图是剪切出的一部分）

可见189、512、238等序号的临界点都对应Na-Cl的BCP，考察哪个都可以，都是等价的。点击Return按钮关闭图形窗口，然后选择7进入临界点属性计算界面，然后再输入临界点序号189，此时就看到了这个临界点处各种函数的值，如下所示。这些函数在Multiwfn手册2.6节都有简要介绍，且在“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/workshop/WFN_content.html）里有十分全面深入详细的讲解。

CP type: (3,-1)
 Density of all electrons:  0.1144132659E-01
 Density of Alpha electrons:  0.5720663296E-02
 Density of Beta electrons:  0.5720663296E-02
 Spin density of electrons:  0.0000000000E+00
 Lagrangian kinetic energy G(r):  0.1168592222E-01
 G(r) in X,Y,Z:  0.7296255010E-03  0.1022667133E-01  0.7296253875E-03
 Hamiltonian kinetic energy K(r): -0.3025433987E-02
 Potential energy density V(r): -0.8660488230E-02
 Energy density E(r) or H(r):  0.3025433987E-02
 Laplacian of electron density:  0.5884598270E-01
 Electron localization function (ELF):  0.1993303578E-01
 Localized orbital locator (LOL):  0.1249064363E+00
...略

如我在全面提到的《Multiwfn支持的分析化学键的方法一览》文中所述，BCP位置上电子密度拉普拉斯函数和能量密度都为正是离子键的典型特征。当前BCP这两个值分别为0.0588 a.u.和0.00302 a.u.，都为正值，完全体现出Na-Cl是离子键的实事。此外，当前BCP的电子密度仅为0.01144 a.u.、ELF仅为0.0199，BCP处很小的电子密度和ELF值也是典型的离子键的普遍共性。

我在《使用Multiwfn+VMD快速地绘制高质量AIM拓扑分析图（含视频演示）》（//www.umsyar.com/445）介绍了一种极好的将Multiwfn和VMD联用快速、理想地绘制AIM拓扑分析图的做法，比在Multiwfn里直接显示的效果更好，而且视角可以自由旋转，没看过的话务必先看一下。对于当前体系，也可以像此文一样用Multiwfn文件包里的examples\scripts\AIM.bat和AIM.txt来作图，得到的图像如下。

上图效果明显不好，显得残缺不全，这是因为此时显示的原子、临界点、键径都只有晶胞里唯一的那些，而处于棱和面上的边界原子、临界点、键径的周期镜像都没有显示出来。为了避免这个问题，对于当前体系有原子、临界点、键径处于棱和面上的情况，应当改用examples\scripts\目录下的AIM+boundary.bat和AIM+boundary.txt来绘制，其用法与AIM.bat和AIM.txt的组合完全一致。只不过前者传递给Multiwfn的指令要求Multiwfn导出mol.pdb、CPs.pdb和paths.pdb文件时把边界的镜像也都纳入进去。利用AIM+boundary.bat和AIM+boundary.txt来绘制的结果如下，可见此时的效果就很好了，和前面在Multiwfn里看到的完全相同。注意此时这三个pdb文件里只有序号靠前的那些才是晶胞中唯一的那些。

2 层状COF

下面再举个例子，用Multiwfn对一个层状的共价有机框架化合物（COF）体系做AIM拓扑分析。此体系的cif文件是本文文件包里的COF-12000N2.cif，其结构图如下所示。

由于此体系的晶胞尺寸不小，计算时可以直接对原胞做只考虑gamma点的计算来得到molden文件。如《实验测定分子结构的方法以及将实验结构用于 计算需要注意的问题》（//www.umsyar.com/569）所说，由于X光衍射测的晶体结构中氢的原子位置通常不准确，最好做分析前先冻结重原子而优化一下氢原子的坐标，但此例为了省事就跳过这一步了。如果你要用AIM拓扑分析考察氢键，由于氢的位置是重中之重，那显然是绝对不能略过这一步的。

当前体系用PBE结合CP2K的GTH赝势基组计算完全可以，而此例作为演示，这回使用PBE结合pob-TZVP-rev2全电子基组计算。启动Multiwfn，载入COF-12000N2.cif，然后输入
cp2k  //创建CP2K输入文件的界面
[回车]  //产生的输入文件名用默认的COF-12000N2.inp
-2  //要求产生molden文件
2  //选择基组
15  //pob-TZVP-rev2
4  //开OT，对当前用的基组通常比对角化明显更容易收敛
0  //产生输入文件

用CP2K计算COF-12000N2.inp，然后将以下晶胞设置信息插入到计算得到的COF-12000N2-MOS-1_0.molden文件的开头
 [Cell]
 22.55600000     0.00000000     0.00000000
-11.27800000    19.53406901     0.00000000
  0.00000000     0.00000000     6.80000000

启动Multiwfn，载入COF-12000N2-MOS-1_0.molden，然后输入
2  //拓扑分析（默认是对电子密度做拓扑分析）
2  //以每个原子核位置为初猜位置搜索临界点
3  //以每两个原子间连线中点作为搜索临界点的初猜位置
4  //以每三个原子的中心位置作为搜索临界点的初猜位置
8  //产生(3,-3)到(3,-1)的拓扑路径
0  //观看结果
现在看到下图。可见临界点和键径都成功地得到了，不仅有化学键对应的BCP和键径，COF内氢键N-H...O作用以及层之间相互作用相对应的BCP和键径都有。而且不仅有晶胞内两层COF之间相互作用的BCP和键径，当前晶胞里的COF和相邻镜像晶胞里的COF相互作用对应的BCP和键径也都有，显然周期性完全充分考虑了。

虽然Multiwfn的文本窗口提示当前并未符合Poincare-Hopf关系，但由于感兴趣的临界点和键径都有了，所以就没必要进一步搜索了。

再利用前述名为AIM+boundary.bat的Windows批处理脚本结合AIM+boundary.txt里记录的指令，用Multiwfn+VMD绘制体系结构+临界点+键径图。然后再在VMD的文本窗口里输入pbc box命令显示盒子，VMD main界面里选择Display - Orthographic用正交视角。此时VMD显示的图像如下，可见非常理想！

顺带一提，对于当前体系要想完整展现层间相互作用的话，最理想的做法是用我提出的IGMH，见《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）和《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（//www.umsyar.com/667）。

3 对Ag晶胞基于电子密度格点数据做AIM拓扑分析

Multiwfn不仅可以基于CP2K产生的波函数文件做AIM拓扑分析，还可以基于从cube等文件里载入的电子密度的格点数据做拓扑分析，此时各个位置的电子密度、梯度和Hessian矩阵基于格点数据以B-spline方式插值得到。由于这个特性，Multiwfn做周期性体系的AIM拓扑分析不限于CP2K用户，诸如从VASP产生的CHGCAR文件里读取电子密度格点数据做AIM拓扑分析也可以。Multiwfn支持哪些格点数据文件格式见手册2.5节的说明。而且CP2K用户还可以在考虑k点的情况下对原胞计算得到电子密度cube文件结合Multiwfn做AIM拓扑分析，而不必在晶胞较小时非得弄成超胞模型再算。但基于电子密度格点数据做AIM拓扑分析的关键不足是对临界点只能得到电子密度及其导数信息以及一些与之直接相关的量，而依赖于波函数才能计算的如动能密度、ELF等函数就没法计算了。另外，如果电子密度格点数据是在赝势下计算产生的，基于它插值出的电子密度函数中，原子核处将没有电子密度极大点，也因此无法产生键径。

此例演示用Multiwfn对Ag的单胞基于电子密度cube文件做AIM拓扑分析的过程。使用的级别是PBE/TZVPP-MOLOPT-PBE-GTH-q19，这个基组是CP2K 2024.1的data目录下的BASIS_MOLOPT_UZH基组文件里的。之所以不用常用的BASIS_MOLOPT里的Ag的基组，在于那里面Ag的基组都是-q11的，即只描述11个价电子。使用赝势的情况下做AIM拓扑分析时，被描述的价电子越多，即使用越小核赝势，结果越准确、越可能接近全电子计算的结果。

启动Multiwfn，载入本文文件包里的Ag单胞的晶体结构文件Ag.cif，然后输入
cp2k  //创建CP2K输入文件的界面
[回车]  //产生的输入文件名用默认的Ag.inp
-3  //要求产生cube文件
1  //电子密度
6  //开smearing
8  //设置k点
10,10,10
0  //产生输入文件

把Ag.inp里BASIS_SET_FILE_NAME设为BASIS_MOLOPT_UZH，POTENTIAL_FILE_NAME设为POTENTIAL_UZH。Ag的BASIS_SET设为TZVPP-MOLOPT-PBE-GTH-q19，POTENTIAL设为GTH-PBE-q19。然后用CP2K计算，之后会得到记录晶胞内价电子密度的Ag-ELECTRON_DENSITY-1_0.cube。如果不了解cube格式的话，看《Gaussian型cube文件简介及读、写方法和简单应用》（//www.umsyar.com/125）。用文本编辑器打开这个cube文件，在第一行里写上box2cell。这样Multiwfn载入这个cube文件时，一旦发现第一行的内容包含box2cell，就会把格点数据的盒子信息转化为晶胞信息，从而在计算的时候能够考虑周期性（也可以不修改cube文件，而是载入cube文件后用隐藏的主功能1000里的子功能18来实现盒子信息向晶胞信息的转化）。

把Multiwfn的settings.ini里的iuserfunc设为-3，使得user-defined function对应基于载入的格点数据通过B-spline插值的函数。然后启动Multiwfn，载入Ag-ELECTRON_DENSITY-1_0.cube，之后输入
2  //拓扑分析
-11  //修改被分析的函数
100  //User-defined function，当前对应于插值方式产生的价层电子密度
3  //将每两个原子之间的中点作为搜索临界点的初猜
然后进入选项0，并把原子标签显示出来、把键设细，就会看到下图。可见临近原子之间的BCP都已经找到了

下面再把更多的临界点也找出来。当前情况不适合选择选项4、5分别用3、4个原子中心作为初猜点来搜索临界点。选项3用每两个原子之间中点当做初猜点时是考虑循环相邻晶胞中的镜像原子的，而选项4、5则只会循环当前晶胞里那些唯一的原子，因此此例用选项4、5搜索后还是会遗漏许多重要的临界点。当前最适合的方式是选择选项6里的子选项-1，此时会在每个原子附近特定半径内撒一批初猜点来搜索临界点。这么选过之后，会发现找到了一大堆新的临界点，进入选项0后看到的图如下所示。下图左侧把所有临界点都显示了出来，可以看到离原子核较近的区域临界点有不少，而且各种类型的都有。这是因为当前用了赝势，由于电子密度不是从原子核开始向四周单调下降的，因此会在价层区域出现各种各样的临界点，这些不是我们感兴趣的。把原子球设大掩盖这些没什么意义的临界点，并且旋转体系使得视角与晶轴平行后，就看到了下面的右图。可见在原子间相互作用区域找到的临界点的分布整齐、对称，因此可以认为有意义的临界点都找到了。

有一个概念叫平面性指数，在“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/workshop/WFN_content.html）里有一页幻灯片做了介绍，如下所示。我们这里来对Ag晶体计算一下这个平面性指数，看看是否和文章里说的情况吻合。注意这里说的平面性指数和《使用Multiwfn定量化和图形化考察分子的平面性（planarity）》（//www.umsyar.com/618）里介绍的我提出的衡量体系几何结构平面性的参数是两码事。

首先我们得确定要考察的笼临界点也即(3,+3)的临界点编号。下图左边只显示了笼临界点，可见当前体系有两类笼临界点，一种是诸如晶胞中心的那个（15号），它处于八面体中心。另一种是处于四面体中心的，比如63号临界点。下图右边只显示了所有BCP，此体系中所有的BCP都是等价的，我们考虑其中任意一个，比如1号。利用当前拓扑分析界面里的选项7依次对这三个临界点的属性进行考察，可以得到的这三个位置的基于格点数据插值产生的电子密度，也即User-defined real space function后面的值，结果如下（注意不要读Multiwfn直接输出的Density of electrons后面的值，那个是用孤立状态原子的电子密度叠加得到的准分子密度）
CP(15): 0.1572533472E-01 a.u.
CP(63): 0.2212960297E-01 a.u.
CP(1): 0.2706526462E-01 a.u.

可见两类笼临界点里电子密度最小值是0.01572 a.u.，BCP处的电子密度为0.02706 a.u.，因此平面性指数为0.01572/0.02706=0.58，在0.5左右，符合幻灯片里说的“其它金属及合金”的情况。

4 总结&其它

本文通过三个有代表性的例子详细介绍了在Multiwfn中对周期性体系做电子密度的拓扑分析，也即一般说的AIM拓扑分析的完整过程，所有要注意的细节都做了解释。通过例子可见，Multiwfn做周期性体系的拓扑分析的过程很简单，并不比分析孤立体系更复杂。读者如果实际操作一遍，还会感受到Multiwfn做这些分析的速度相当快。AIM拓扑分析对于考察周期性体系的电子结构有重要意义，推荐读者在实际研究中运用，使文章增光添彩。

很值得一提的是，Multiwfn的拓扑分析功能极其普适，绝不仅限于电子密度的拓扑分析。对于Multiwfn自身可以直接计算的函数，以及载入的格点数据里记录的函数，Multiwfn都可以做拓扑分析，对分子体系的例子参看比如《使用Multiwfn对静电势和范德华势做拓扑分析精确得到极小点位置和数值》（//www.umsyar.com/645）对静电势、范德华势的拓扑分析，以及《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）里面对ELF-pi函数的拓扑分析。对周期性体系，Multiwfn也完全可以以相同的方式对任意函数做拓扑分析。

使用本文的做法做AIM拓扑分析在发表文章时必须按照Multiwfn启动时的提示恰当引用Multiwfn，如果用于给别人代算也要明确告知对方这一点。

文/Sobereva@北京科音  2024-Jul-1

1 前言

笔者之前写过《使用Multiwfn做电荷分解分析(CDA)、绘制轨道相互作用图》（//www.umsyar.com/166）介绍了怎么用Multiwfn程序基于Gaussian等 程序产生的波函数文件对分子体系做电荷分解分析（CDA）。这种方法可以从片段轨道间相互作用的角度深入了解片段间是如何转移电子的。由于CDA的很高价值，Multiwfn做CDA分析已经得到非常广泛的使用。如果读者没看过此//www.umsyar.com/166的话一定要先仔细看一下，否则无法理解下文的内容。

从2024-Jun-5更新的Multiwfn开始，CDA模块已经支持了基于CP2K产生的molden文件做周期性体系的CDA分析，下面将通过一个简单的例子进行演示。Multiwfn可以在主页//www.umsyar.com/multiwfn免费下载，不了解此程序者参看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。我假定读者已经具备了CP2K的基本使用知识，不了解者推荐通过北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）完整系统地学习。

本文的例子是计算NaCl板吸附CO分子的体系，用CDA方法考察CO与NaCl板之间电子转移情况。这个体系优化后的结构如下，C和O分别是1和2号原子，3-110号原子是NaCl板。对这个体系，之前我在《使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线》（//www.umsyar.com/638）中通过密度差图的方式定性考察了电子转移，在《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）里通过片段电荷的方式定量考察了电子转移。而下面用CDA分析，可以以轨道角度在明显更深层次理解电子转移细节。

本文例子涉及的所有文件都在//www.umsyar.com/attach/716/file.7z里。CP2K用的版本是2024.1，Multiwfn是2024-Jun-5更新的版本，读者绝对不要用Multiwfn的更老版本。

2 产生输入文件

首先我们要用CP2K计算NaCl+CO、NaCl、CO各自的molden格式的波函数文件。如果你不知道怎么用CP2K产生molden文件的话，先仔细阅读《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651），我假定读者已经看过本文了。这里有几个要点：
(1)整体和片段的molden文件里的原子坐标必须对应
(2)整体和片段的计算必须使用严格相同的级别和设置，用的盒子也应当一致
(3)必须要求CP2K计算所有空轨道
(4)整体的molden文件里需要加入[Cell]字段定义晶胞信息，而片段的molden文件里加不加入不影响CDA结果。[Nval]字段可加入可不加入，也不影响CDA结果

本文文件包里的NaCl_CO.cif是以前我在PBE-D3(BJ)/DZVP-MOLOPT-SR-GTH下优化的NaCl板吸附CO的结构文件。用Multiwfn载入它，然后输入
cp2k  //创建CP2K输入文件
NaCl_CO.inp  //产生的输入文件名
-2  //要求产生molden文件
-9  //其它设置
12  //设置计算的空轨道数
-1  //计算所有空轨道
0  //返回
0  //产生输入文件

用CP2K运行NaCl_CO.inp，得到NaCl_CO-MOS-1_0.molden。然后将以下字段插入到molden文件的开头
[Cell]
16.92168000     0.00000000     0.00000000
 0.00000000    16.92168000     0.00000000
 0.00000000     0.00000000    25.00000000

将NaCl_CO.inp复制为NaCl.inp，再将NaCl.inp中&COORD里的CO部分删掉、PROJECT NaCl_CO改名为PROJECT NaCl，然后用CP2K运行之，得到NaCl-MOS-1_0.molden。将上面的[Cell]加入其中（这和CDA分析无关，加入这个的目的是为了之后在Multiwfn中能正确观看NaCl板的周期性的晶体轨道）。

将NaCl_CO.inp复制为CO.inp，再将CO.inp中&COORD里的NaCl部分删掉、PROJECT NaCl_CO改名为PROJECT CO，然后用CP2K运行之，得到CO-MOS-1_0.molden。这个可以不用刻意加上[Cell]，因为当前体系中CO离盒子边界很远，显示轨道时是否考虑周期性对看到的图像没影响。

接下来就可以进行CDA分析了。

3 进行CDA分析

启动Multiwfn，载入NaCl_CO-MOS-1_0.molden，然后依次输入
16  //CDA分析
2  //两个片段
CO-MOS-1_0.molden  //CO的波函数文件。注意要按整个体系里原子序号顺序载入片段
NaCl-MOS-1_0.molden  //NaCl板的波函数文件

现在屏幕上输出了一大堆复合物轨道的CDA结果，同时还输出了ECDA结果：
The net electrons obtained by frag. 2 = CT( 1-> 2) - CT( 2-> 1) =    0.1580
这个ECDA结果告诉你CO（片段1）往NaCl板（片段2）转移了0.158个电子，这和《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）里通过Mulliken方法计算的CO的片段电荷严格相同。

当前CDA分析输出的复合物轨道太多，一共437个（对应整个体系的占据轨道数），而绝大多数对d、b、r项的贡献微乎其微或完全为0，没有考察的意义。为了简化输出，在当前界面里输入
-3  //设置CDA输出时d、b、r项的阈值，它们之一的绝对值大于这个阈值的复合物轨道才会输出
0.005
0  //显示CDA和ECDA结果

当前结果如下
   Orb.      Occ.          d           b        d - b          r
    272    2.000000    0.006220   -0.000012    0.006232    0.000488
    275    2.000000    0.125862   -0.001874    0.127736    0.051155
    297    2.000000    0.001098    0.000305    0.000792   -0.005373
    318    2.000000    0.001053    0.000326    0.000727   -0.005720
    341    2.000000    0.000781    0.000200    0.000581   -0.006708
    352    2.000000    0.002291    0.000731    0.001560   -0.015705
    360    2.000000    0.001101    0.000298    0.000803   -0.009297
    414    2.000000    0.000635    0.000214    0.000421   -0.005229
    427    2.000000    0.000739    0.000219    0.000520   -0.005731
-------------------------------------------------------------------
Sum:     874.000000    0.153805    0.017732    0.136074   -0.037831
Note: The "Sum" includes all terms including those not printed above

CDA的总结果是CO向NaCl转移了d=0.154个电子，NaCl向CO反馈了b=0.018个电子，NaCl净获得d-b=0.136个电子。b项基本可以忽略不计，而d的主要贡献是由275号复合物轨道造成的。这个复合物轨道的贡献显然值得进一步探究，看看是由CO和NaCl板的哪些片段轨道混合所造成的。为此，输入
5  //将特定复合物轨道对CDA项的贡献进行分解
275  //要考察的复合物轨道序号
0.005  //输出阈值
现在看到如下结果

Occupation number of orbital   275 of the complex:  2.00000000
FragA Orb(Occ.)  FragB Orb(Occ.)      d           b        d - b          r
   5( 2.0000)     271( 2.0000)    0.000000    0.000000    0.000000    0.005183
   5( 2.0000)     289( 2.0000)    0.000000    0.000000    0.000000    0.005642
   5( 2.0000)     336( 2.0000)    0.000000    0.000000    0.000000    0.005127
   5( 2.0000)     341( 2.0000)    0.000000    0.000000    0.000000    0.008960
   5( 2.0000)     435( 0.0000)    0.009358    0.000000    0.009358    0.000000
   5( 2.0000)     439( 0.0000)    0.011934    0.000000    0.011934    0.000000
   5( 2.0000)     443( 0.0000)    0.010929    0.000000    0.010929    0.000000
   5( 2.0000)     447( 0.0000)    0.008032    0.000000    0.008032    0.000000
   5( 2.0000)     453( 0.0000)    0.009177    0.000000    0.009177    0.000000
   5( 2.0000)     465( 0.0000)    0.010755    0.000000    0.010755    0.000000
   5( 2.0000)     469( 0.0000)    0.013735    0.000000    0.013735    0.000000
   5( 2.0000)     473( 0.0000)    0.011282    0.000000    0.011282    0.000000
Sum of above terms:               0.085201    0.000000    0.085201    0.024911

可以看到CO向NaCl转移电子靠的都是CO的5号占据轨道和NaCl空轨道相互作用。NaCl板接收CO转移来的电子分散在其大量空轨道上，比如NaCl板的443号空轨道接受了0.011个电子。在当前的输出阈值下，上面输出的这些项的所有的d项加和为0.085，明显和275号复合物轨道的总d值0.153有很大差距，说明CO还向上面没输出的很多其它的NaCl板的空轨道转移了电子。

下图把275号复合物轨道、CO的5号占据轨道，以及NaCl接收CO贡献来的电子最多的两个轨道（439、469）展示在了一起，便于大家直观了解情况，等值面数值用的是0.02。复合物和片段的轨道可以用Multiwfn载入相应的波函数文件后用主功能0按照《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）说的观看。从此图可以清晰直观地认识到CO的孤对电子往吸附CO的那个Na位点的空轨道上转移了电子。还有大量也具有这种特征的NaCl板的空轨道也接收了CO的孤对电子，这里没全都画出来。

下面再看一下上图中275号复合物轨道是怎么由片段轨道构成的。按0返回到CDA菜单，然后进入选项2，输入275，看到以下信息

 Note: Only the fragment orbitals with contribution >  1.0 % will be shown below, the threshold can be changed by "compthresCDA" in settings.ini

 Occupation number of orbital   275 of the complex:  2.00000000
 Orbital     5 of fragment  1, Occ: 2.00000    Contribution:   87.42 %
 Sum of values shown above:     87.42 %

这告诉你275号复合物轨道由CO的5号占据轨道贡献了87.4%。默认情况下只有贡献大于1%的片段轨道才会输出，虽然NaCl板的一大堆非占据轨道都通过与CO的5号轨道混合而参与了这个复合物轨道的构成，但贡献都小于输出阈值，所以上面没看到。从当前信息可以认识到275号复合物轨道主要体现CO的5号轨道特征，这是为什么上面图中275号复合物轨道和CO的5号轨道看起来十分相似。

4 总结

此文通过一个简单的表面吸附的例子演示了怎么用Multiwfn的CDA功能分析周期性体系的电荷转移本质。讨论电子转移常见的套路就是绘制电子密度差、算算片段电荷，而从本文的例子看到Multiwfn可以从片段轨道相互作用对电子转移的本质进行深入的考察，掌握了本文介绍的方法明显可以给分析增光添彩。

使用Multiwfn做本文介绍的分析请记得在发表文章时引用Multiwfn程序启动时提示的Multiwfn原文，以及引用介绍CDA/GCDA方法的文章（Multiwfn用的实质上是我对CDA广义化后的GCDA形式），这在进入Multiwfn的CDA功能时屏幕上明确提示了。

文/Sobereva@北京科音   2024-Jun-8

1 前言

原子电荷（atomic charge）即原子所带的净电荷，对于讨论化学体系中原子的带电状态有重要意义，作为原子的一种描述符它也和原子的很多属性有密切联系。强烈建议阅读《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！》（//www.umsyar.com/714）里提到的笔者写的原子电荷的综述以全面了解原子电荷。

强大的波函数分析程序Multiwfn（//www.umsyar.com/multiwfn）支持非常丰富的原子电荷计算方法，其中多数已经支持周期性体系。Multiwfn支持的方法中有一类是基于模糊式原子空间划分，也就是每个原子对应一个平滑的权重函数，权重函数与体系的电子密度的乘积的全空间积分就是这个原子带的电子数（原子的布居数），原子的核电荷数减去它就是原子电荷。Multiwfn支持的这类方法包括Hirshfeld、Hirshfeld-I、MBIS、Becke，以及对Hirshfeld电荷做后校正得到的ADCH和CM5电荷。从2024-May-25更新的Multiwfn版本开始，Hirshfeld、CM5、Hirshfeld-I、MBIS都已经支持了周期性体系，本文将具体介绍怎么用Multiwfn与非常流行、高效且免费的第一性原理程序CP2K相结合非常方便快速地计算这些原子电荷。Multiwfn计算这些电荷的功能是普适的，并不限于结合CP2K，也可以基于其它第一性原理程序产生的体系的电子密度的cube文件来算，还支持VASP产生的记录电子密度的CHGCAR文件。
注：Multiwfn更老的一些版本也支持基于周期性波函数计算Hirshfeld和CM5电荷，但那时候用的算法对于周期性体系效率极低，本文介绍的做法与之有天壤之别。

下文第2节简要介绍一下Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷的基本特点，第3节介绍用Multiwfn+CP2K计算它们的具体方法，第4节给出具体例子，同时还会介绍怎么计算片段电荷。使用Multiwfn计算原子电荷在发表文章时记得需要按Multiwfn启动时的提示恰当引用Multiwfn。

Multiwfn可以在官网//www.umsyar.com/multiwfn免费下载，如果对Multiwfn不了解的话建议看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。如果你用CP2K创建给Multiwfn算原子电荷用的输入文件的话，笔者假定你对CP2K的使用已经有基本且正确的了解，不具备这些知识的话非常建议通过北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）完整、系统地学习一遍。本文的例子利用到Multiwfn创建CP2K的输入文件，相关介绍见《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）。

2 Hirshfeld、Hirshfeld-I、CM5和MBIS等原子电荷简介

Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷在Multiwfn手册3.9节的相应小节里都有很详细的介绍，因此这里不介绍细节，只是简单说一下它们的特点。

Theor. Chim. Acta (Berl.), 44, 129 (1977)中提出的Hirshfeld电荷是被广为使用的原子电荷，计算方式简单，物理意义清楚，在 研究孤立体系和第一性原理研究周期性体系的领域中都用得很多。Hirshfeld电荷为人诟病的地方是原子电荷数值整体严重偏小，对偶极矩、静电势重现性很差，这在我写的《原子电荷计算方法的对比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）里的对比测试中有充分体现。

JCP, 126, 144111 (2007)提出的Hirshfeld-I是在Hirshfeld基础上的改进，通过迭代方式不断更新原子权重函数直到收敛，这使得原子权重函数所描述的原子空间可以响应实际化学环境。其原子电荷数值显著大于Hirshfeld，对静电势的重现性也有明显改进。代价是Hirshfeld-I需要做迭代，往往几十轮，计算耗时明显高于Hirshfeld，而且需要事先提供体系中各种元素的各个氧化态的原子径向密度信息，实现起来较麻烦（为了实现Hirshfeld-I，我计算了360多个原子的径向密度文件并在Multiwfn中自带）。

一定要注意，虽然CP2K程序自己也支持Hirshfeld-I电荷的计算，但我测试发现其结果完全不靠谱，我检查了其源代码也确认了其实现方式明显不符合Hirshfeld-I的标准定义，因此切勿用CP2K自身的Hirshfeld-I电荷计算功能。CP2K算的Hirshfeld电荷在结合SHAPE_FUNCTION DENSITY选项时的结果不那么离谱，但结果和Multiwfn算的也存在不可忽视的差异，应当以Multiwfn的结果为准。

JCTC, 8, 527 (2012)中提出的CM5电荷和我在J. Theor. Comput. Chem., 11, 163 (2012)中提出的ADCH电荷一样都是对Hirshfeld电荷的后校正，对大部分情况都能使得它的原子电荷变得更大、更符合化学直觉，对静电势的重现性也明显变得更好。ADCH没有CM5那么强的经验性、不牵扯一堆经验参数，因而原理更为理想、更推荐使用。CM5最初是用于分子体系的，如今在应用于晶体方面也已经有了一定探索，例如JCTC, 16, 5884 (2020)将CM5与其它一些原子电荷对比考察了它们对一些固体的电荷分布的描述。CM5还被用于一些计算化学领域的方法，比如《比SMD算溶解自由能更好的溶剂模型uESE的使用》（//www.umsyar.com/593）里介绍的uESE溶剂模型依赖于CM5电荷，CM5乘以1.2后适合结合OPLS-AA力场做动力学模拟，见《计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒人脚本》（//www.umsyar.com/585）。

ADCH电荷在分子体系的研究方面已被广为使用，之前有人问我是否会把ADCH电荷扩展到周期性体系。我目前没这个打算，因为ADCH电荷的思想对周期性体系往往不适用。ADCH方法将原子偶极矩展开为周围原子的校正电荷，但对于诸如NaCl这样的体系，每个原子所处的环境是中心对称的，因此原子偶极矩为0、校正电荷都为0，故ADCH电荷会与Hirshfeld电荷完全相同。

MBIS (Minimal Basis Iterative Stockholder)电荷于JCTC, 12, 3894 (2016)中提出。类似于Hirshfeld-I，此方法也是通过迭代方式（因此耗时较高）优化原子空间直到收敛。MBIS比Hirshfeld-I的一大好处是不需要利用各个元素各个氧化态的原子径向密度信息，因此实现起来简单得多，也更优雅。原文里的测试体现出MBIS比Hirshfeld-I整体更具有优势。但MBIS经常收敛较慢，对于某些体系甚至需要200多轮才能充分收敛。目前Multiwfn里MBIS最高支持到Rn元素。

根据我的测试，对于大多数固体体系，原子电荷大小的关系是MBIS和Hirshfeld-I最大、谁大不一定，CM5大小介中，Hirshfeld电荷最小。如果你对数值的绝对大小不那么在乎，只是想体现相对大小关系的话，最“朴素”且很便宜的Hirshfeld也可以用。但如果对定量数值关心的话，建议用其它的。CM5、Hirshfeld-I和MBIS用于晶体哪个最理想没有定论，都可以算一算看看，如果倾向于得到较大数值的结果可以优先考虑Hirshfeld-I和MBIS。如果其中有的原子电荷明显不合理那就弃掉换别的方法，比如O的电荷如果算出来-2.34那就不应该用，因为显然O最多只能带-2电荷。

特别要强调的是，绝对不要把原子电荷与氧化态搞混，也切勿用原子电荷判断氧化态！氧化态怎么正确地判断看《使用Multiwfn结合CP2K计算晶体中原子的氧化态》（//www.umsyar.com/711）和《使用Multiwfn通过LOBA方法计算氧化态》（//www.umsyar.com/362）。

还有很多其它原子电荷计算方法，诸如Mulliken电荷、Lowdin电荷、SCPA电荷、AIM电荷（被个别文章称为Bader电荷）等，Multiwfn也都支持将它们用于周期性体系，这不属于本文的范畴。Mulliken电荷结合CP2K里常用的MOLOPT系列和pob系列基组用于周期性体系也说得过去。虽然它是基于希尔伯特空间划分的原子电荷中最原始、思想最朴素的一种，但实际结果一般还算定性正确，而且计算耗时极低。被用得很多的AIM电荷实际上很糟糕，往往出现违背化学常识的结果，这在前述的《原子电荷计算方法的对比》里已经充分体现了，不建议用。

3 Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷在Multiwfn中的计算方法

3.1 CP2K用户的情况

首先说对于CP2K用户怎么计算。Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷都是完全基于电子密度计算的，不直接牵扯到波函数，在Multiwfn中计算它们有两种可以用的输入文件：
(1)用CP2K产生的体系的电子密度cube文件作为输入文件。这种情况在CP2K计算时可以考虑k点，是我最推荐的，耗时非常低。如果不了解cube文件的格式的话，参看《Gaussian型cube文件简介及读、写方法和简单应用》（//www.umsyar.com/125）。注意CP2K有bug（起码直到我撰文时最新的CP2K 2024.1为止），产生的cube文件里的原子有效电荷都为0，这会导致Multiwfn无法正确计算原子电荷。这里说的原子的有效核电荷是指当前基组描述的原子的电子数，对于赝势基组它对应价电子数。为此，要么用文本编辑器手动改cube文件中每个原子的有效核电荷，要么把cube文件第一行写成各个元素的有效核电荷，例如N 5 B 3 Ti 12，这代表当前体系中的N、B、Ti的有效核电荷数分别为5、3、12
(2)用CP2K产生的记录了体系波函数的molden文件作为输入文件。详见《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）。这个方法的缺点在于没法考虑k点，对于小晶胞体系必须先扩胞到足够大，从而在计算时只考虑gamma点，这显然使得在CP2K中和Multiwfn中的计算耗时都会比用(1)的方式高很多

用Multiwfn载入上述输入文件其中一种后，进入主功能7，然后选择相应选项即可计算Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷。对于molden文件当输入文件时，Multiwfn还会让你输入计算用的格点间距，格点间距越小耗时越高而计算精度越高，通常用0.2 Bohr就可以，是精度和耗时的较好权衡。而使用cube文件做为输入文件时，计算原子电荷用的格点间距与此文件的格点间距直接对应，CP2K的CUTOFF设得越大则产生的cube文件里的格点间距会越小。

由于Multiwfn计算上述电荷用的是均匀分布的格点，不可能准确积分较重原子的内核区域电子密度，因此上述输入文件应当是CP2K做GPW计算得到的，也需要用赝势基组，此时cube文件记录的电子密度或者molden文件里记录的波函数只对应于价电子。

顺带一提，如果你要计算Mulliken电荷，用Multiwfn载入上述的molden文件，然后进主功能7，依次选择5、1即可得迅速到结果。Multiwfn中计算AIM电荷的方法笔者以后会另文说明。

3.2 其它程序用户的情况

如果你是VASP用户，可以用记录晶胞中价电子密度的CHGCAR文件给Multiwfn用于计算前述原子电荷，等同于3.1节说的第2类输入文件。只要文件名里包含CHGCAR字样就会被Multiwfn当做是CHGCAR文件来载入。注意为了让Multiwfn得知各类原子的有效核电荷数，CHGCAR文件的第一行需要以Nval为开头，后面依次写上各类元素的名字和有效核电荷，并以空格分隔。比如写为Nval Na 1 Cl 7，这代表Na和Cl的有效核电荷数分别为1和7。

如果你是其它第一性原理程序的用户，若程序可以产生记录晶胞中价电子密度的cube文件（即等同于3.1节说的第2类输入文件），也都可以用于给Multiwfn计算前述的原子电荷。

之后在Multiwfn中的计算过程和CP2K用户没任何区别。

4 实例

下面给出Multiwfn结合CP2K计算Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷的具体例子。相关的文件都可以在//www.umsyar.com/attach/712/file.zip中得到（4.2、4.3和4.5节的cube/molden文件除外，因为文件过大）。

4.1 TiO2晶体

此例对TiO2的金红石晶型计算前述的原子电荷。它的实验晶体结构的cif文件是本文文件包里的TiO2-Rutile.cif，结构图如下所示

先用CP2K计算出这个晶胞的电子密度的cube文件，用于之后Multiwfn做原子电荷计算。首先用Multiwfn创建CP2K输入文件，启动Multiwfn并载入TiO2-Rutile.cif，然后输入
cp2k  //创建CP2K的输入文件
[回车]  //产生的文件名用默认的TiO2-Rutile.inp
-3  //产生cube文件
1  //电子密度
8  //设置k点
4,4,6  //与各个晶胞边长相乘皆为18埃左右，够用了
0  //产生输入文件，计算级别为默认的PBE/DZVP-MOLOPT-SR-GTH

用CP2K运行刚得到的TiO2-Rutile.inp，马上就算完了，得到了TiO2-Rutile-ELECTRON_DENSITY-1_0.cube文件，在本文的文件包里已经提供了。用文本编辑器打开此文件，把第一行改为Ti 12 O 6并保存。这代表Ti和O的有效核电荷分别为12和6，直接对应于输入文件里它们分别用的DZVP-MOLOPT-SR-GTH-q12和DZVP-MOLOPT-SR-GTH-q6基组名里面的q后面的值。

启动Multiwfn，输入
TiO2-Rutile-ELECTRON_DENSITY-1_0.cube  //写此文件的实际路径
y  //从第一行载入各个元素有效核电荷
7  //布居分析与原子电荷计算
1  //Hirshfeld电荷

立刻看到如下Hirshfeld电荷计算结果。由于格点间距不是无穷小，因此等价的原子算出来的原子电荷也会有轻微的数值差异，可以自己手动取平均。如果用更大的CUTOFF，则cube文件的格点间距会更小，不等价性会更低。
Final atomic charges:
Atom    1(Ti):     0.59349459
Atom    2(O ):    -0.29461903
Atom    3(O ):    -0.29472601
Atom    4(Ti):     0.58363165
Atom    5(O ):    -0.29385236
Atom    6(O ):    -0.29385236
现在Multiwfn还问你是否把原子电荷导出为chg文件，这里选n不导出。

接下来计算CM5电荷。在Multiwfn主功能7里选择16，马上就得到了结果
                      ---------- CM5 charges ----------
Atom:    1Ti  CM5 charge:    1.348393  Hirshfeld charge:    0.593495
Atom:    2O   CM5 charge:   -0.672093  Hirshfeld charge:   -0.294619
Atom:    3O   CM5 charge:   -0.672191  Hirshfeld charge:   -0.294726
Atom:    4Ti  CM5 charge:    1.338580  Hirshfeld charge:    0.583632
Atom:    5O   CM5 charge:   -0.671306  Hirshfeld charge:   -0.293852
Atom:    6O   CM5 charge:   -0.671306  Hirshfeld charge:   -0.293852
Summing up all CM5 charges:     0.00007649
可见CM5电荷的大小明显大于Hirshfeld电荷，看起来更合理，而Hirshfeld电荷值明显偏小。

下面再演示计算Hirshfeld-I电荷。把Multiwfn自带的examples目录下的atmrad目录挪到当前目录下，这样Multiwfn计算Hirshfeld-I电荷时就会自动用此目录下的各种元素各个价态的径向电子密度。之后在主功能7里选15，再选1用默认的设置进行Hirshfeld-I电荷计算，经过16轮迭代得到了结果：
Final atomic charges:
Atom    1(Ti):     3.01789464
Atom    2(O ):    -1.50888824
Atom    3(O ):    -1.50886131
Atom    4(Ti):     3.01787620
Atom    5(O ):    -1.50897240
Atom    6(O ):    -1.50897240
可见Hirshfeld-I电荷的大小不仅显著大于Hirshfeld电荷，比CM5电荷还要大很多。虽然数值很大，但完全在合理范围内，毕竟都没有大过原子的氧化态。

接下来再计算MBIS电荷。在主功能7里选择20，再选1用默认的设置进行MBIS电荷计算，经过19轮迭代得到了结果：
Final atomic charges:
Atom    1(Ti):     1.73473636
Atom    2(O ):    -0.86733169
Atom    3(O ):    -0.86731253
Atom    4(Ti):     1.73467041
Atom    5(O ):    -0.86734303
Atom    6(O ):    -0.86734303
对于当前例子，MBIS的数值介于CM5和Hirshfeld-I之间。但也有很多反例，MBIS也经常比Hirshfeld-I电荷还大不少，后文有例子。

以上原子电荷的计算顺序是随意的，不要误以为必须按上述顺序计算，你想算哪个就算哪个即可。

4.2 AlN晶体

这一节计算AlN晶体的原子电荷。完全可以按照和上一节相同的做法基于电子密度的cube文件来计算，但本节演示一下基于超胞的molden文件的计算过程。AlN原胞的cif文件是本文文件包里的AlN.cif，将之载入Multiwfn，然后输入以下命令创建CP2K输入文件
cp2k  //创建CP2K输入文件
AlN_443.inp  //产生的输入文件名
-11  //进入几何操作界面
19  //构造超胞
4  //a方向扩胞成原本的4倍
4  //b方向扩胞成原本的4倍
3  //c方向扩胞成原本的3倍
-10  //返回
-2  //要求产生molden文件
0  //产生输入文件

用CP2K运行刚产生的AlN_443.inp，得到AlN_443-MOS-1_0.molden。在此文件开头插入以下内容定义晶胞和各个元素的价电子数
[Cell]
12.44400000     0.00000000     0.00000000
-6.22200000    10.77682012     0.00000000
 0.00000000     0.00000000    14.93400000
[Nval]
Al 3
N 5

启动Multiwfn，载入AlN_443-MOS-1_0.molden，然后输入
7  //布居分析与原子电荷计算
16  //CM5电荷
[回车]  //使用均匀分布的格点计算
[回车]  //使用默认的0.2 Bohr格点间距。格点间距越小结果越精确但计算越耗时
在i9-13980HX CPU上24核并行，50秒就算完了，结果为

                      ---------- CM5 charges ----------
Atom:    1Al  CM5 charge:    0.467945  Hirshfeld charge:    0.302813
Atom:    2N   CM5 charge:   -0.467945  Hirshfeld charge:   -0.302812
Atom:    3Al  CM5 charge:    0.467945  Hirshfeld charge:    0.302812
Atom:    4N   CM5 charge:   -0.467946  Hirshfeld charge:   -0.302814
...略

Al的电负性很小，N的电负性很大，因此二者的原子电荷理应相差悬殊。当前的Hirshfeld电荷明显偏小。CM5电荷的数值虽然更大，但还是显得有些偏低。

下面再计算Hirshfeld-I电荷。在主功能7里输入
15  //Hirshfeld-I
1  //开始计算
[回车]  //使用默认的0.2 Bohr格点间距

经过37轮收敛，结果如下。可见Hirshfeld-I电荷的数量级比较合理，体现出AlN中成键的极强离子性
 Atom    1(Al):     1.80317384
 Atom    2(N ):    -1.80317397
 Atom    3(Al):     1.80317257
 Atom    4(N ):    -1.80317539
...略

最后再计算下MBIS电荷。在主功能7里输入
20  //MBIS
1  //开始计算
[回车]  //使用默认的0.2 Bohr格点间距

对此体系MBIS收敛很慢，经过190轮才达到收敛，在i9-13980HX上24核并行算了5分钟，在双路7R32 96核并行下算了不到一分钟。结果如下。可见对此体系，MBIS电荷的数值很大，和氧化态几乎正好是一样的。所以，如前所述，若你就是想算出来比较大的原子电荷，可以先考虑MBIS，但也有一些例外，后面会提到。
Atom    1(Al):     2.99599978
Atom    2(N ):    -2.99601488
Atom    3(Al):     2.99600129
Atom    4(N ):    -2.99599935

4.3 MOF-5晶体

MOF-5晶体含有Zn、O、C、H。其cif文件是本文文件包里的MOF-5.cif。由于此体系的晶胞边长约26埃，因此计算电子密度格点数据时不用考虑k点。按照3.1节的做法，基于PBE/DZVP-MOLOPT-SR-GTH级别的电子密度cube文件（第一行改为Zn 12 O 6 C 4 H 1），用Multiwfn计算各种原子电荷，结果如下。此体系中O有两类，分别给出，C则处于较多不同的化学环境，其电荷就不给出了。

Hirshfeld：
Zn=0.372 O=-0.332/-0.200 H=0.056

CM5：
Zn=0.644 O=-0.584/-0.305 H=0.116

Hirshfeld-I（36轮收敛，2*7R32机子上96核并行4分多钟算完）
Zn=1.568 O=-1.917/-0.746 H=0.139

MBIS（39轮收敛，2*7R32机子上96核并行3分多钟算完）
Zn=1.546 O=-1.662/-0.782 H=0.190

对这个体系，MBIS和Hirshfeld-I相符得很好，而且Zn的电荷大小很合理。Hirshfeld电荷再次数值最小、CM5比之大一些。

4.4 其它：MoS2、BaTiO3、MgF2、NaCl

注意MBIS电荷并非总是较大。例如《使用Multiwfn结合CP2K计算晶体中原子的氧化态》（//www.umsyar.com/711）里的单层MoS2的例子，在PBE/DZVP-MOLOPT-SR-GTH级别的波函数下，Mo的电荷的计算结果是Hirshfeld=0.330，CM5=0.827，Hirshfeld-I=1.987，MBIS=-0.059，Mulliken=-0.771。对这个体系来说CM5和Hirshfeld-I显得靠谱，MBIS电荷接近0有点说不通，Mulliken电荷明显有误导性。

BaTiO3在PBE/DZVP-MOLOPT-SR-GTH级别下的原子电荷计算结果如下。Hirshfeld还是严重偏小。Hirshfeld-I大得离谱，Ba的原子电荷都超过氧化态+2了。相较之下CM5和MBIS的结果比较靠谱，正好二者算的O的原子电荷差不多，差异主要在Ba和Ti的电荷的相对大小上。
Atom:    1Ba  CM5 charge:    1.103253  Hirshfeld charge:    0.404357
Atom:    2Ti  CM5 charge:    1.200467  Hirshfeld charge:    0.537212
Atom:    3O   CM5 charge:   -0.767963  Hirshfeld charge:   -0.310211
Atom:    4O   CM5 charge:   -0.767892  Hirshfeld charge:   -0.315693
Atom:    5O   CM5 charge:   -0.767892  Hirshfeld charge:   -0.315693

Hirshfeld-I
Atom    1(Ba):     2.37169239
Atom    2(Ti):     3.02878350
Atom    3(O ):    -1.78799416
Atom    4(O ):    -1.80625420
Atom    5(O ):    -1.80625434

MBIS：
Atom    1(Ba):     0.57391985
Atom    2(Ti):     1.69416825
Atom    3(O ):    -0.74678425
Atom    4(O ):    -0.76066530
Atom    5(O ):    -0.76066537

再看两个典型离子化合物的情况，都在PBE/DZVP-MOLOPT-SR-GTH级别下计算。对NaCl计算的Na的原子电荷为：Hirshfeld=0.217，CM5=0.443，Hirshfeld-I：1.053，MBIS：1.043，Mulliken：0.596
对MgF2计算的Mg的原子电荷为：
Hirshfeld=0.414，CM5=0.828，Hirshfeld-I：2.048，MBIS：1.888，Mulliken：1.439
可见对于强离子性化合物，Hirshfeld-I和MBIS电荷都在氧化态附近，并且由于方法定义的缺陷，往往电荷的大小还轻微大于氧化态的大小（但如果用大核赝势，对Na和Mg只描述最外层的3s价电子，则不会有这个问题）。其它原子电荷的大小关系为Mulliken>CM5>Hirshfeld。

以MgF2为例，这里顺带一提CP2K算的Hirshfeld-I电荷明显不对，千万不要用。在CP2K输入文件的&DFT/&PRINT里加上以下内容要求输出Hirshfeld-I电荷后，CP2K直接给出的Mg的Hirshfeld-I电荷为0.522，和Multiwfn算的相距甚远。
&HIRSHFELD
  SHAPE_FUNCTION DENSITY
  SELF_CONSISTENT T
&END HIRSHFELD

4.5 NaCl板吸附CO的片段电荷

NaCl板吸附CO是《使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线》（//www.umsyar.com/638）里的一个例子体系，这里基于这个模型计算一下原子电荷和CO片段的电荷。CP2K优化这个体系后产生的NaCl_CO-1.restart文件在本文的文件包里提供了，将之载入Multiwfn，然后输入
cp2k  //创建CP2K输入文件
NaCl_CO.inp  //产生的输入文件名
-7  //设置周期性
XY
-3  //要求产生cube文件
1  //电子密度
-2  //产生molden文件（用于之后算Mulliken电荷的目的）
0  //产生输入文件

计算完毕后，把NaCl_CO-ELECTRON_DENSITY-1_0.cube文件第一行改为Na 9 Cl 7 C 4 O 6。先用Hirshfeld方法各个原子的电荷和CO的片段电荷。启动Multiwfn，载入此cube文件，然后输入
7  //布居分析与原子电荷计算
-1  //定义片段
109,110  //CO的原子序号
1  //Hirshfeld电荷
结果为
...略
 Atom  109(C ):     0.13622804
 Atom  110(O ):    -0.03514734

Fragment charge:    0.10108070
Fragment population:    9.89891930

即CO的Hirshfeld片段电荷为0.101。

之后再选择CM5方法，结果为
...略
 Atom:  109C   CM5 charge:    0.147306  Hirshfeld charge:    0.136228
 Atom:  110O   CM5 charge:   -0.081972  Hirshfeld charge:   -0.035147
 Summing up all CM5 charges:     0.00510030

 Fragment charge:    0.06533361

之后再选择Hirshfeld-I方法，结果为
 Atom  109(C ):     0.09790054
 Atom  110(O ):    -0.08744136

 Fragment charge:    0.01045918

之后再选择MBIS方法，结果为
...略
 Atom  109(C ):     0.19642171
 Atom  110(O ):    -0.15937966

 Fragment charge:    0.03704205

再来算Mulliken电荷。用Multiwfn载入加入了恰当的[Cell]和[Nval]信息的CP2K计算后产生的NaCl_CO-MOS-1_0.molden文件，输入
7  //布居分析与原子电荷计算
-1  //定义片段
109,110  //CO的原子序号
5  //Mulliken电荷
1  //在屏幕上输出结果
结果为
...略
Atom   109(C )    Population:  3.60671601    Net charge:  0.39328399
Atom   110(O )    Population:  6.23484007    Net charge: -0.23484007
Total net charge:   -0.00000086

Fragment charge:    0.15844392

可见不同方法计算的结果在定量上有一定差异，但都有共性，即CO整体带少量正电，即电子从CO往NaCl板上转移了一些。并且C和O分别带正电和负电，且原子电荷的绝对值都是C大于O。

5 总结

本文简要介绍了四种知名的基于模糊式原子空间划分计算原子电荷的方法的特点，包括Hirshfeld、Hirshfeld-I、CM5和MBIS，并介绍了它们在Multiwfn中的计算方法，并且给了很多具体例子，以使得读者能了解如何使用Multiwfn结合CP2K或其它第一性原理程序很容易地计算它们。同时还说明了怎么直接计算片段电荷，这是定量讨论体系中两部分间电子转移的最严格的方式。推荐将这些原子电荷应用于实际问题的研究当中以考察固体与表面体系的电荷分布情况。记得届时应恰当引用Multiwfn启动时提示的Multiwfn原文以及相应的原子电荷的原文。

通过本文的一些体系的对比可以明显看到，Hirshfeld电荷原理上最简单，但明显缺点是整体严重偏小，CM5通常比它大不少，而且电荷的稳定性不错。虽然结合化学直觉来看，CM5电荷的绝对大小仍然往往能偏小，但至少也不会明显离谱，总比Hirshfeld电荷更建议使用。Hirshfeld-I和MBIS电荷普遍比较大，谁更大不一定，至少几乎总比CM5还大不少，从绝对大小来看往往比CM5更合理，但稳健性比CM5弱一些，对极个别体系可能会表现得不符合常识。对于离子性特征很强的化合物，Hirshfeld-I和MBIS电荷往往接近氧化态，由于方法并不完美，有时候原子电荷甚至比氧化态还要大一点点。当然，没有任何原子电荷计算方法是绝对完美的，本文介绍的这些原子电荷在实际当中都可以用，可以根据它们的特点和实际结果决定选用，在横向对比时方法必须统一。

顺带一提，对于计算固体表面或者多孔物质的原子电荷以用于基于力场的动力学模拟的目的，我目前最推荐CP2K算REPEAT电荷，北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）里讲“电子结构的分析”的部分做了专门的介绍。但REPEAT电荷需要在真空区定义拟合点，因而无法用于诸如本文所示的AlN、TiO2等致密的固体，而且对于远离真空区的原子也没法拟合出质量靠谱的原子电荷。所以REPEAT电荷有特定用处，但不是像本文介绍的那些原子电荷一样属于普适性的原子电荷计算方法。

文/Sobereva@北京科音  2024-May-23

1 前言

波函数分析程序Multiwfn可以基于电子波函数计算原子和片段的氧化态，这在笔者之前写的《使用Multiwfn通过LOBA方法计算氧化态》（//www.umsyar.com/362）一文中已做了详细介绍，没读过者一定要先阅读此文，否则无法理解下文的内容。此文介绍的做法已经被广泛用来准确地考察分子中的氧化态。如今Multiwfn已经支持将这种方法用于CP2K计算的周期性体系的波函数，从而考察晶体中原子的氧化态，这非常有实际价值，具体做法将在本文介绍。注意氧化态和原子电荷完全是两码事，周期性体系的原子电荷的计算看《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）。

Multiwfn可以在官网//www.umsyar.com/multiwfn免费下载，本文的读者必须使用2024年5月23日及以后更新的版本，否则没有本文介绍的特性。对Multiwfn不了解者看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。

2 用法

在Multiwfn中对周期性体系做LOBA或mLOBA分析氧化态的具体操作流程和对孤立体系做别无二致，都是先用主功能19产生定域化轨道（介绍见《Multiwfn的轨道定域化功能的使用以及与NBO、AdNDP分析的对比》//www.umsyar.com/380），然后进入主功能8（轨道成分分析），选择LOBA/mLOBA分析，之后输入一个LOBA方法判断定域化轨道电子归属的阈值，或者输入m让Multiwfn按照mLOBA的方式判断归属即可。定域化轨道的电子归属是什么意思，在前述的《使用Multiwfn通过LOBA方法计算氧化态》里已经专门说过了，这里不再累述。之后各个原子的氧化态就会立马输出出来。用户也可以自定义片段来得到片段的氧化态。

对周期性体系做LOBA/mLOBA分析的输入文件必须是CP2K产生的molden文件，并且需要自己在里面加入盒子信息。具体方法在《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）里有详细介绍。由于考虑k点时不能产生molden文件，因此原本晶胞较小的话需要扩胞成足够大的超胞，从而对其计算时只需要考虑gamma点。Multiwfn不支持Quantum ESPRESSO、VASP等其它第一性原理程序产生的波函数，也没有任何支持的必要，因为对于绝大多数情况，用Multiwfn创建CP2K单点输入文件并用CP2K计算是非常简单且快速的事，可参考《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）。

对于当前目的，CP2K计算时用赝势基组也行，用全电子基组也行，没有什么特别的讲究。CP2K计算的时候不要开smearing，否则轨道占据数不再精确为整数、不再是标准的单行列式波函数，此时也就没法做轨道定域化。如果有特殊原因非开smearing不可，可以在Multiwfn载入molden文件后进入主功能6，选择选项38，此时Multiwfn会按照轨道能量由低到高重排电子确定占据数，轨道占据数就都为整数了，之后退回到主菜单，再照常按前述步骤做LOBA/mLOBA分析即可。

3 例子

下面给出一系列通过Multiwfn使用LOBA/mLOBA计算固体的氧化态的例子。CP2K的输入输出文件在//www.umsyar.com/attach/711/file.zip里都给了。只有3.1节例子对应的molden文件我在里面直接给了，如果所有例子的molden文件都给的话压缩包就太大了。CP2K输入文件里的各种设置我这里就不细说了，在北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）里有非常全面系统的讲解。

3.1 单层氮化硼

氮化硼（BN）是知名的二维材料，这里完整介绍单层BN的氧化态计算的流程。本文文件包里的BN_bulk.cif是氮化硼晶体的结构文件，用GaussView打开它，可看到晶胞里有两层BN。删除其中的一层，然后保存为BN_ml.gjf。之后启动Multiwfn，载入BN_ml.gjf，输入
cp2k   //创建CP2K输入文件
[回车]  //产生的输入文件用默认的名字BN_ml.inp
-2  //要求导出molden文件
-11  //进入几何操作界面
19  //构造超胞
6  //第1个平行于层的方向扩胞成原先的6倍
6  //第2个平行于层的方向扩胞成原先的6倍
1   //垂直于层的方向不扩胞
-10  //返回
-7  //设置周期性
XY
0  //产生输入文件（计算级别为默认的PBE/DZVP-MOLOPT-SR-GTH）

用CP2K运行Multiwfn产生的BN_ml.inp，我这里用双路7R32 96核机子4秒钟就算完了。现在当前目录下有了BN_ml-MOS-1_0.molden，按照《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）所述，把以下信息插入进去，以定义晶胞信息和原子价电子数信息。改好的文件已经在本文的文件包里提供了，是BN_ml-MOS-1_0.molden。

 [Cell]
14.98944000     0.00000000     0.00000000
-7.49472000    12.98123580     0.00000000
 0.00000000     0.00000000     8.00000000
 [Nval]
B 3
N 5

启动Multiwfn，载入BN_ml-MOS-1_0.molden，然后输入
19  //轨道定域化
1  //只定域化占据轨道，用默认的Pipek-Mezey方法
8  //轨道成分分析
100  //LOBA/mLOBA分析

输入一个阈值，比如50，马上看到各个原子的以LOBA方法计算的氧化态，如下所示。可见B的氧化态为3、N为-3，整个体系所有原子氧化态加和为0，非常符合常识。这里用的阈值50%对于成键方式比较典型的情况是适合的，它的意思是如果有某原子对某个定域化轨道贡献超过50%（Multiwfn以Hirshfeld方法计算贡献值），则这个轨道的电子就完全归属于这个原子。

 Oxidation state of atom   1(B ) :  3
 Oxidation state of atom   2(N ) : -3
 Oxidation state of atom   3(B ) :  3
 Oxidation state of atom   4(N ) : -3
 Oxidation state of atom   5(B ) :  3
 Oxidation state of atom   6(N ) : -3
...略
 The sum of oxidation states:     0

我提出的modified LOBA (mLOBA)方法是对LOBA的重要改进，它将每个定域化轨道的电子完全归属到对它贡献最大的那个原子上，可以确保所有原子的氧化态加和总为0，而且也没有阈值选择的任意性（有的时候LOBA的结果对阈值敏感），比LOBA在原理上明显更好。在当前界面里输入m即得到mLOBA方法算的氧化态的结果，也是B和N的氧化态分别为3和-3。

3.2 BaTiO3

BaTiO3晶体的cif文件是本文文件包里的BaTiO3.cif。用Multiwfn创建其4*4*4超胞在PBE/DZVP-MOLOPT-SR-GTH下做单点的CP2K计算的输入文件。当前的超胞在每个方向有16埃，足够大了。CP2K计算后得到molden文件，在里面加入恰当的[Cell]和[Nval]字段后，启动Multiwfn并将之载入，然后输入
19  //轨道定域化
-9  //设置产生定域化轨道后自动计算轨道成份的方法
0  //不输出轨道成份。这些信息不是当前感兴趣的，当前体系又大，不输出可以节约不少时间
1  //只定域化占据轨道，用默认的Pipek-Mezey方法
由于当前体系大、轨道数很多，因此轨道定域化过程耗时较高。在EPYC 7R32机子上花了约8分钟。之后输入
8  //轨道成分分析
100  //LOBA/mLOBA分析
m  //显示mLOBA方法计算的原子氧化态

马上看到以下结果。可见Ba、Ti、O的氧化态都和常识完全一致

 Oxidation state of atom   1(Ba) :  2
 Oxidation state of atom   2(Ti) :  4
 Oxidation state of atom   3(O ) : -2
 Oxidation state of atom   4(O ) : -2
 Oxidation state of atom   5(O ) : -2
...略

3.3 MOF-5

MOF-5是知名的金属有机框架化合物，介绍见https://en.wikipedia.org/wiki/MOF-5。其中含有Zn，这里算算Zn的氧化态是多少。MOF-5晶体的cif文件在本文的文件包里提供了，是立方晶胞，边长25.86埃。由于边长已经大到可以只考虑gamma点做计算，因此CP2K计算时不需要扩胞，直接对原胞算单点得到molden文件即可。相应的CP2K输入文件已提供在了本文的文件包里。氧化态计算过程同前，7R32机子上几分钟可以算完整个过程。氧化态结果为Zn=2、O=-2、H=1、C=2/0/2。可见Zn、O、H的氧化态都完全符合预期。当前体系中C的氧化态有三种，和所处的化学环境差异有关。

3.4 单层MoS2

MoS2是知名的二维材料体系，它的体相结构文件是本文文件包里的MoS2_bulk.cif。完全效仿3.1节氮化硼的例子，取单层并基于5*5*1的超胞模型做它的氧化态的计算。如果以LOBA方法计算氧化态并输入50作为阈值，结果为

 Oxidation state of atom   1(Mo) :  6
 Oxidation state of atom   2(S ) : -2
 Oxidation state of atom   3(Mo) :  6
 Oxidation state of atom   4(S ) : -2
 Oxidation state of atom   5(S ) : -2
 Oxidation state of atom   6(S ) : -2
 Oxidation state of atom   7(Mo) :  6
...略
 The sum of oxidation states:   50

虽然Mo的价电子有6个，看似Mo的氧化态为6也不是不可能，但实际上明显是错的，因为当前所有原子的氧化态加和为100。考虑到S的氧化态不可能比-2更负，明显当前的结果高估了Mo的氧化态。即便尝试各种阈值，也始终得不到能令LOBA方法算的氧化态加和为0的情况。因此LOBA对此体系完全失败。

输入m使用mLOBA计算氧化态，结果为

 Oxidation state of atom   1(Mo) :  4
 Oxidation state of atom   2(S ) : -2
 Oxidation state of atom   3(S ) : -2
 Oxidation state of atom   4(Mo) :  4
 Oxidation state of atom   5(S ) : -2
 Oxidation state of atom   6(S ) : -2
 Oxidation state of atom   7(Mo) :  2
 Oxidation state of atom   8(S ) : -2
 Oxidation state of atom   9(S ) : -2
...略
 The sum of oxidation states:     0

可见mLOBA方法算的所有原子氧化态加和为期望的0，但是Mo有的氧化态是2有的是4，和当前体系中所有Mo等价的事实不符，原因后面会说。解决方法是输入-1定义片段，然后输入所有Mo对应的序号。序号的快速查询方法是在Multiwfn主功能0的菜单里选择Tools - Get atom indices of a given element，然后输入Mo，点OK，就会返回下面的内容，将之粘贴到Multiwfn定义片段的窗口里即可。
1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58,61,64,67,70,73
片段定义好后，再次输入m，可以看到Oxidation state of the fragment: 100，即这个片段总的氧化态为100。由于当前体系里有25个完全等价的Mo，因此Mo的氧化态为100/25=4，是完全正确的值。

为什么mLOBA算出来的明明等价的Mo会有不同的氧化态？这可以通过考察定域化轨道的特征结合mLOBA的原理予以了解。默认情况下，Multiwfn的主功能19做完轨道定域化后会自动输出各个定域化轨道的成份，可以看到除了单中心、双中心的定域化轨道外，还有大量三中心的：

 More delocalized LMOs: (Three largest contributions are printed)
  301:   64(Mo) 24.8%   46(Mo) 24.8%   61(Mo) 24.8%
  302:   46(Mo) 24.8%   49(Mo) 24.8%   31(Mo) 24.8%
  303:   16(Mo) 24.8%   34(Mo) 24.8%   31(Mo) 24.8%
...略

从轨道成份上明显可知这些三中心轨道都是等价的。轨道定域化做完后按照《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）说的，在Multiwfn的主功能0里随便看一个三中心轨道，如下所示

可见这个轨道体现的是三个Mo之间的三中心作用。按照mLOBA的思想，这个轨道上的两个电子会完全归属到对它贡献最大的原子上，但由于三个Mo的贡献本质上相同，实际的贡献量由于数值计算误差会有极其轻微的差异，所以这个轨道上的电子归属到哪个Mo是有随机性的，这是为什么Mo的氧化态有的为2有的为4。可见，只要把轨道定域化和mLOBA的思想都理解了，遇到可疑的结果自然就能通过具体分析搞明白是怎么回事、知道怎么恰当解决。

4 总结

本文详细介绍了怎么用Multiwfn基于CP2K对固体做DFT计算得到的波函数通过LOBA和mLOBA方法得到原子的氧化态，给了很多例子。可见本文的方法操作简单，结果非常合理，这对于讨论原子在固体中的状态、对体系进行分类很有价值。作为练习，读者可以对NaCl、MgO、AlN、SiO2晶体也都算算氧化态。值得一提的，如本文的例子所体现的，如有了我后来提出的mLOBA方法，原本的LOBA就没什么使用价值了，我建议总是默认用mLOBA。

使用本文的方法计算在发表结果时务必记得按照Multiwfn程序启动时的提示对Multiwfn进行正确的引用。

Key points and simple examples of hybrid functional calculations with CP2K

文/Sobereva@北京科音

 First release: 2023-Nov-7   Last update: 2024-Mar-9

0 前言

对周期性体系，杂化泛函的计算耗时众所周知远高于纯泛函。著名的第一性原理程序CP2K不光纯泛函的计算超级快，对大体系做杂化泛函计算也非常快（尽管还是远比纯泛函贵），远远快于Quantum ESPRESSO (QE)、VASP等纯粹基于平面波的程序。CP2K结合像样的基组做杂化泛函的单点计算，在一般双路服务器上算到上千原子都没问题。然而，使用CP2K做高效的杂化泛函计算有很多非常关键性的要点必须知道，远远不是像Gaussian程序那样光写个杂化泛函名字就了事的。然而，笔者经常在网上回答CP2K问题时看到很多人完全不具备这些最基本常识，胡用瞎用CP2K做杂化泛函计算，导致耗时巨高、始终出不来结果、任务崩溃、SCF完全不收敛等各种问题。因此笔者觉得很有必要写个小文将CP2K做杂化泛函计算的最关键常识和要点挨个提一下。Multiwfn程序（//www.umsyar.com/multiwfn）具有极其便利的创建CP2K输入文件的功能，在本文的最后会还用Multiwfn产生C60晶体的CP2K杂化泛函计算的输入文件作为实际例子。

本文的内容至少对CP2K 2024.1，以及2023-Nov-7及以后更新的Multiwfn是适用的。

笔者在北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）的“能量的计算及相关问题”部分里讲杂化泛函计算远比本文全面、深入得多，并列举诸多经验技巧和例子，还详细介绍HSE06、CAM-B3LYP、wB97X等各种范围分离形式的杂化泛函以及revDSD-PBEP86-D3(BJ)等各种双杂化泛函的用法和要点，使用杂化泛函结合k点计算的要点，欢迎想学习更多知识者关注和参加。

1 CP2K杂化泛函计算的特性

---------
2024-Feb-18更新：以下之前写的关于k点的说法已经过时。从CP2K 2024.1版开始支持了杂化泛函计算时考虑k点，称为RI-HFXk算法，并且可以无缝结合ADMM降低耗时。此时需要以RI近似计算HF交换部分。不仅能算能量还能做优化（包括变胞）和能带计算。从第3届的“北京科音CP2K第一性原理计算培训班”开始已经加入了RI-HFXk的全面深入的介绍和具体例子，参考http://bbs.keinsci.com/thread-43683-1-1.html。

首先要知道CP2K做杂化泛函计算的一个关键不足是不支持考虑k点，因此如果被计算的是周期性体系，晶胞必须大到只考虑gamma点就足够。如果你要计算的体系的晶胞原本比较小，显然必须先扩胞到足够大才行。如果不知道扩到多大够用，严格来说需要做你感兴趣的属性的计算结果随扩胞倍数的收敛性测试。虽然CP2K能够对含有非常多原子的超胞做得动杂化泛函计算，但耗时会明显高于QE等程序对原胞在考虑k点的情况下做杂化泛函计算。

最能充分体现CP2K杂化泛函计算效率优势的情况是算一些较大分子晶体、MOF、COF等一些大晶胞体系，以及做第一性原理动力学等情况。对小晶胞体系如上所述经过扩胞后也能算，只不过不体现CP2K的长处。

CP2K显然没法用杂化泛函算能带，因为不支持k点，这种目的只能用QE等程序代替。必须用CP2K算的话，应当用CP2K支持的纯泛函里算带隙相对比较好的比如HLE17。
---------

由于杂化泛函远比普通泛函昂贵，尤其是用于几何优化、NEB、动力学等涉及到结构变化的任务。因此如果你算的问题用纯泛函就足以描述得不错、用杂化泛函不会带来明显好处的话，显然应该用纯泛函。记住，非必要甭考虑杂化泛函，纯泛函的耗时仅有杂化泛函的一个微不足道的零头。对几何优化、振动分析、NEB等高耗时任务用纯泛函，算势垒、反应能的时候用更好的杂化泛函是可以接受且常见的策略。

周期性体系的杂化泛函计算不仅耗时远高于纯泛函计算，对内存要求更是远远高于纯泛函计算。经常用CP2K做杂化泛函任务的服务器的内存容量绝对是多多益善，选择服务器的时候要注意这一点。也别天真地试图拿个人计算机将就着跑CP2K的杂化泛函计算。

CP2K做杂化泛函计算的前提是编译时必须带着libint电子积分库，否则没法算杂化泛函中HF交换项中涉及的双电子积分。CP2K的编译方法见//www.umsyar.com/586。

杂化泛函用的赝势基组和赝势就用常见的给PBE等纯泛函优化的就可以。从CP2K 9.1开始，data目录下的BASIS_MOLOPT_UZH和POTENTIAL_UZH文件里还分别给出了专门给PBE0优化的赝势基组和赝势，如C的DZVP-MOLOPT-PBE0-GTH-q4和GTH-PBE0-q4，原理上来说它们用于杂化泛函计算会更理想一些（但实际未必有什么优势）。

2 库仑截断

杂化泛函计算之所以比纯泛函昂贵得多，在于计算杂化泛函的HF交换项时要计算大量的双电子积分（two-electron integral, ERI），尤其是对于周期性体系来说，不光要计算晶胞内基函数之间的ERI，还要计算中心晶胞与周围镜像晶胞基函数之间的ERI。再加上由于只能算gamma点而晶胞必然较大，原子数、基函数通常非常多，无疑要算的ERI是巨量的。

由于周期边界条件，原理上来说要计算的ERI是无穷的，显然这会导致计算无法进行。因此实际计算时必须使用库仑截断（Coulomb truncation），也就是将ERI中描述电子库仑相互作用的1/r算符截断在某个距离（截断半径），超过截断半径时1/r视为0，相应地ERI就不需要计算了，这个策略使得ERI的计算数目是有限的。如果你计算的是非周期性体系，由于本身ERI数目就是有限的，因此虽然也可以用库仑截断，但这就不是必须的了。

杂化泛函的库仑截断的半径设置显然既影响耗时也影响结果。截断半径越小，要计算的ERI越少，耗时就越低，但结果越偏离泛函原本的精确结果。显然不能为了一味地节约时间而把截断半径设得太小，否则结果会很垃圾。6埃通常是耗时和精度的较好权衡，设得明显更小会令精度有明显损失，设得更大并不会令精度有明显提升却会增加耗时不少。如果对精度要求很高也可以用8埃，设得更大就完全没必要了。

要注意像6、8埃这种程度的截断半径其实还没有大到令绝对能量充分收敛，但能量在相互求差时，由截断半径带来的系统性误差能很大程度抵消。可以认为用不同的截断半径相当于用不同的理论方法，直接对比能量的情况必须保持截断半径的设置相统一！

值得一提的是，不一定截断半径小结果就一定差，有可能截断半径造成的对ERI忽略的误差和泛函本身的误差相抵消导致恰好对某些体系某些问题的计算比用更大截断半径时结果还更好。

截断半径应当小于晶胞最短尺寸的一半以免导致得到缺乏物理意义的结果。这里所谓的最短尺寸，是(001)晶面间距、(010)晶面间距、(100)晶面间距三者中的那个最小值。也因此为了能够使用6埃截断半径，晶胞最短尺寸应当大于12埃，如果晶胞不够大则应当扩胞（即便抛开库仑截断问题不谈，单从只能计算gamma点这点来说，晶胞也不能过小）。

Multiwfn创建的杂化泛函输入文件中，若周期性设成了NONE，则不使用库仑截断，若是周期性，Multiwfn会自动加上库仑截断的设置。如果晶胞最短尺寸够大，Multiwfn会自动把截断半径设为6埃，如果不够大，则会根据当前晶胞最短尺寸自动设成能设的最大半径。

以下是库仑截断的设置方式，加入到&DFT/&XC里面。POTENTIAL_TYPE TRUNCATED代表用库仑截断的1/r算符，CUTOFF_RADIUS设置截断半径（埃）。T_C_G_DATA t_c_g.dat是默认的可以不写，这指定的是库仑截断用的gamma函数的数据文件t_c_g.dat，在CP2K的data目录下。
&HF
  &INTERACTION_POTENTIAL
    POTENTIAL_TYPE TRUNCATED
    CUTOFF_RADIUS 6.0
    T_C_G_DATA t_c_g.dat
  &END INTERACTION_POTENTIAL
&END HF

3 积分屏蔽

积分屏蔽（integral screening）是指尽可能在不明显牺牲精度的前提下减少要计算的ERI的量，这对于CP2K杂化泛函计算的耗时和内存使用量有至关重要的影响！上一节的库仑截断本身就是一种积分屏蔽策略，CP2K里还有另外两种非常重要的积分屏蔽方法，也是普遍应用在Gaussian等 程序中的：
(1)利用Schwarz不等式做积分屏蔽：利用Schwarz不等式关系，可以根据数目较少的二指标ERI估计出数目甚巨的三、四指标ERI中哪些数值肯定小于阈值从而可以忽略掉，这可以大幅减少ERI要算的数目（原本总共要算的ERI中占比最大的就是四个基函数序号不同的那些ERI）。阈值通过&DFT/&XC/&HF/&SCREENING中EPS_SCHWARZ设置，此值越大，被忽略的ERI数目越多、计算精度越差。默认的1E-10虽然足够保证精度，但此时太贵，没必要。通常设为1E-6就可以了，是精度和耗时的较好权衡，想要精度更好点用1E-8。
(2)结合密度矩阵做积分屏蔽：ERI对能量的贡献取决于它与特定密度矩阵元的乘积。因此哪怕ERI不是很小，但与之相乘的密度矩阵元很小，那么这样的ERI也完全不需要计算。CP2K定义了特定方法快速估计一个ERI与相应密度矩阵元乘积值的上限，如果小于上述的EPS_SCHWARZ，则此ERI就会忽略掉。这种依赖于密度矩阵的积分屏蔽默认是关闭的，如果开启则把&DFT/&XC/&HF/&SCREENING里的SCREEN_ON_INITIAL_P设为T。由于这种积分屏蔽能令耗时巨幅降低，因此我强烈建议总是开启，也因此Multiwfn产生的杂化泛函的输入文件里都自动把SCREEN_ON_INITIAL_P设成T。

基于密度矩阵做积分屏蔽有一个极其关键的要点是初猜的密度矩阵质量不能太差！然而很多人对此一无所知！默认情况下初始的密度矩阵是CP2K自动做初猜产生的，由于它和SCF收敛时的密度矩阵相差很大、质量很糙，因此此时积分屏蔽通常做得极烂，不仅可以忽略的ERI可能没忽略，还会把一些重要的ERI给忽略掉，由此可能导致KS矩阵质量很烂、SCF极难收敛。而且忽略哪些ERI是在SCF过程一开始决定好的，即便之后随着SCF迭代密度矩阵质量不断变好，也不会由此使得积分屏蔽逐步变得合理，这导致即便最后SCF收敛了，得到的能量的精度也可能很烂。因此当SCREEN_ON_INITIAL_P设T时，必须读取一个基本合理的波函数当初猜。通常是在杂化泛函计算前先用相同的基组做一个便宜的纯泛函计算，比如PBE0计算前先用PBE做一个单点计算，得到记录了PBE收敛的波函数信息的wfn文件，然后PBE0计算时用SCF_GUESS RESTART，并把WFN_RESTART_FILE_NAME指定成那个wfn文件，以读取其作为初猜。这不仅使得基于密度矩阵的积分屏蔽可以合理进行，还有一个额外的好处是由于初猜波函数质量较好（至少比默认自动产生的好得多），使得杂化泛函计算时达到SCF收敛所需的迭代次数会比较少，这进一步节约了时间（注意即便抛开ERI计算耗时不谈，杂化泛函在构造KS矩阵时花的时间也明显高于纯泛函，因为涉及到密度矩阵元与ERI的大量的乘加操作）。

一般的 程序虽然也基于密度矩阵做积分屏蔽，但不是非得像CP2K这样得读取一个像样的初猜波函数，这是因为它们会根据当前最新的密度矩阵判断这一轮SCF用的ERI哪些是能被忽略的。CP2K的做法和它们不同，主要是因为CP2K倾向于让用户做incore形式的SCF而非一般量化程序一般用的direct形式的SCF，见后文。

当SCREEN_ON_INITIAL_P设T时，对于杂化泛函的几何优化、分子动力学等任务的续算，要注意光有.restart文件还不行，还必须提供上一步的.wfn文件当初猜波函数，否则也是会由于基于自动初猜的密度矩阵做积分屏蔽非常恶心而令任务无法正常跑下去。

如果你嫌麻烦而不想先做纯泛函计算产生收敛的波函数再做杂化泛函计算，那就不要写SCREEN_ON_INITIAL_P T。对于小体系、小基组特别是非周期性的计算，不基于密度矩阵做积分屏蔽时耗时也不高。而对于计算量较大的杂化泛函计算，则强烈建议不要偷这个懒而多花巨量不必要的时间。

4 ADMM

CP2K还支持其开发者独创的ADMM (Auxiliary Density Matrix Methods)方法减少要算的ERI数目从而进一步显著加快杂化泛函的计算，其原理见以下北京科音CP2K第一性原理计算培训班的ppt。

ppt里提到的主基组就是指平时实际用的基组，诸如DZVP-MOLOPT-SR-GTH。对于不同的主基组，有一些专门构造的与之匹配的ADMM辅助基组，比如MOLOPT系列基组适合搭配的是CP2K的data目录下的BASIS_ADMM_UZH文件里的那些辅助基组，从小到大依次是admm-dz、admm-dzp、admm-tzp、admm-tz2p。admm-dzp是2-zeta结合极化函数的辅助基组，通常算是保底质量。DZVP-MOLOPT-SR-GTH主基组适合搭配admm-dzp，更好的主基组如TZVP-MOLOPT-GTH结合admm-dzp也可以接受，若结合admm-tzp精度会更好些，这都比起不用ADMM能巨幅节约计算时间和内存。辅助基组的大小与主基组越接近，ADMM给杂化泛函HF交换项带来的误差越小，但需要计算越多ERI因而越耗时。如果辅助基组和主基组设为相同的，那么ADMM不会造成任何误差，相应地也不会节约任何耗时（反倒还会增加一些耗时）。

data目录下的BASIS_ADMM文件中的FIT3、pFIT3等辅助基组都比较过时了，对元素涵盖也远不如BASIS_ADMM_UZH里的完整，因此如今不再建议使用（它们出现在不少官方例子文件和文献中，如今不要效仿）。目前版本的Multiwfn产生利用ADMM的杂化泛函的输入文件时默认设的是admm-dzp。

注意前述这些辅助基组都是对GTH赝势基组用的，千万别做全电子计算时也用它们。要想全电子计算时用ADMM，可以用(aug-)pcseg系列基组作为主基组，在《在线基组和赝势数据库一览》（//www.umsyar.com/309）中介绍的BSE基组数据库中有和它匹配的(aug-)admm辅助基组，诸如pcseg-2适合搭配admm-2。pcseg是一类构造得很好的对DFT计算非常划算的基组，《谈谈 中基组的选择》（//www.umsyar.com/336）里简要提及了，在北京科音中级 培训班（http://www.keinsci.com/workshop/KBQC_content.html）里讲基组的地方做了具体介绍。

用MOLOPT系列基组做周期性体系的杂化泛函计算几乎是一定要结合ADMM的，否则极难算得动！因为MOLOPT是完全广义收缩基组，而CP2K利用的计算源高斯函数（PGTF）之间ERI的libint库又没有专门考虑广义收缩的情况，因此造成等效的PGTF数目甚巨，ERI的数目又与PGTF数目的四次方近似成正比，数目显然多得恐怖。如果你用片段收缩基组，比如def2-SVP、6-31G**、pcseg-1、pob-DZVP-rev2之类去做杂化泛函计算，由于需要考虑的PGTF数目远比档次与之相仿佛的DZVP-MOLOPT-SR-GTH要少，因此ADMM方法不是非用不可，但不用的话耗时会显著高于DZVP-MOLOPT-SR-GTH搭配admm-dzp的情况。

不能有的原子用ADMM加速计算而有的原子不用。因此如果你用混合基组，其中有的原子用的主基组有ADMM辅助基组，有的没有，你又想用ADMM，那么没有ADMM辅助基组的那些原子的辅助基组就只能设定成其主基组。虽然这些原子没法被ADMM技术加速，但至少其它那部分原子能被加速。

前面提到的各种节约耗时的策略是彼此完全兼容的，同时利用可以最大程度节约时间。即周期性体系杂化泛函计算时一般要用库仑截断、利用Schwarz不等式做积分屏蔽、结合密度矩阵做积分屏蔽、开启ADMM。所有这些设定都会对能量产生影响，因此写文章的时候最好注明，横向对比时应统一。

启用ADMM需要在&KIND里加上BASIS_SET AUX_FIT [辅助基组名]，同时增加一行BASIS_SET_FILE_NAME [辅助基组文件名]。并且在&DFT里加入以下内容
&AUXILIARY_DENSITY_MATRIX_METHOD
&END AUXILIARY_DENSITY_MATRIX_METHOD
如果用的是对角化而非OT，还需要在以上字段里加入ADMM_PURIFICATION_METHOD NONE要求不做纯化。

5 ERI的储存和内存使用量的控制

SCF有两种常见形式，incore SCF是指SCF计算一开始将所有要用到的ERI全都算出来并储存在内存里，之后SCF每一轮在构造KS矩阵时会直接从内存里读取ERI而不再重新计算。direct SCF是指SCF每一轮都重新计算要用到的ERI，这也叫on-the-fly方式计算ERI。显然，incore比direct总耗时低得多，但代价是机子的必须内存非常大，因为要存的ERI通常非常多。在CP2K的杂化泛函计算中，如果你给的内存足够存得下所有ERI，那么只有SCF第一轮的时候由于要计算所有ERI所以耗时很高，而之后SCF每一轮的耗时就很低了。如果给的内存不够储存所有ERI，那么能在内存里存多少ERI就存多少，存的那部分在之后每一轮SCF中会直接读取，存不下的ERI在每一轮SCF中会重新算，相当于介于incore和direct之间的情况。显然，做大体系杂化泛函计算时应当尽量在大内存机子上做、分配给CP2K尽可能多的内存，以让尽可能多的ERI能存到内存中、尽可能减少每轮SCF过程中要on-the-fly计算的ERI的量来节约耗时。

常有人问为什么他做杂化泛函计算时老是卡在SCF刚开始的位置不动，这是因为SCF第一轮肯定是要计算所有需要考虑的ERI，耗时显然往往很高，对大体系、大基组的情况花一、两个小时都是常事，看上去好像一直停在这里似的。如果你给的内存较少而导致大部分ERI都得on-the-fly方式去算，那么之后SCF每一轮都得花将近这么多时间，整个任务跑完可能得半天甚至一天时间。如果PRINT_LEVEL设的是MEDIUM，在SCF第一轮算完后就会输出Number of sph. ERI's calculated on the fly和Number of sph. ERI's stored in-core，前者是之后每一轮SCF都需要重新算的ERI数，后者是已存到内存中而之后不需要SCF每一轮重新算的ERI数。显然前者相对于后者越小，之后每一轮SCF比第一轮时间少得越多。

CP2K用户通常默认用popt版本，每个MPI进程在计算HF交换项时最大内存使用量由&DFT/&XC/&HF/&MEMORY里的MAX_MEMORY控制，单位为MB，扣除掉储存零碎数据的占用量外，其余都用来储存ERI。MAX_MEMORY与MPI并行进程数的乘积必须小于当前机子空余物理内存量，显然应当尽量给大以尽可能减少on-the-fly方式算的ERI量，但也不要太顶着头分配，因为进入HF交换项计算模块之前CP2K还会占一定量的内存，故要有适当余量。在SCF第一轮不断往内存里写入ERI、内存占用量不断增大的过程中，若空余内存已用完时CP2K还在继续往里存ERI，CP2K就会马上崩溃，甚至还可能导致计算机暂时失去响应。所以MAX_MEMORY应当设得恰到好处。

如果空余物理内存很有限，为了能让MAX_MEMORY能设得比较大，可以减小MPI进程数。为了不因此浪费CPU运算能力，此时可以用psmp版，此时每个MPI进程下属的OpenMP线程都可以共享这个MPI进程储存的ERI。

6 杂化泛函的定义方式

这一节简要说一下杂化泛函怎么正确在输入文件里定义。对于非周期性体系，比如用PBE0，直接在&DFT里写上以下内容即可。此时不用库仑截断而用完整的1/r算符，EPS_SCHWARZ为默认的1E-10、SCREEN_ON_INITIAL_P为默认的F。
&XC
  &XC_FUNCTIONAL PBE0
  &END XC_FUNCTIONAL
&END XC

显然对周期性体系绝对不能简单写成上面那样，至少也得自己去指定为使用库仑截断。结合前面的讲解，周期性体系PBE0计算的泛函定义部分一般写为下面这样。PBE0泛函由75%的PBE交换项+25%的HF交换项+100%的PBE相关项构成，所以需要恰当地定义&PBE里面的参数，并在&HF里用FRACTION指定HF成份。
&XC
  &XC_FUNCTIONAL
    &PBE
      SCALE_X 0.75
      SCALE_C 1.0
    &END PBE
  &END XC_FUNCTIONAL
  &HF
    FRACTION 0.25
    &SCREENING
      EPS_SCHWARZ 1E-6
      SCREEN_ON_INITIAL_P T
    &END SCREENING
    &INTERACTION_POTENTIAL
      POTENTIAL_TYPE TRUNCATED
      CUTOFF_RADIUS 6.0
    &END INTERACTION_POTENTIAL
    &MEMORY
      MAX_MEMORY 3000
    &END MEMORY
  &END HF
&END XC

有的泛函在&XC_FUNCTIONAL里有特定的字段，可以省得自己去定义细节参数，但HF成份依然必须自己定义，否则相当于只做了残缺的纯泛函计算，结果毫无意义。比如HF成分为10%的TPSSh泛函可以写为
&XC
  &XC_FUNCTIONAL   
    &HYB_MGGA_XC_TPSSH
    &END HYB_MGGA_XC_TPSSH
  &END XC_FUNCTIONAL
  &HF
    FRACTION 0.1
    ...同上
  &END HF
&END XC

更多的泛函定义方式可以直接看Multiwfn产生的相应杂化泛函的输入文件里的写法，都是严格正确的。

7 一个常见警告：The Kohn Sham matrix is not 100% occupied

做杂化泛函计算，特别是用了ADMM时，很容易在SCF开始后看到The Kohn Sham matrix is not 100% occupied警告。具体原因在https://www.cp2k.org/faq:hfx_eps_warning有介绍，这里就不多说了。对于实际计算来说，看到这个警告时，先把&QS里的EPS_PGF_ORB设为比默认明显更小的1E-12，如果计算正常，即便还有这个警告也可直接无视。EPS_PGF_ORB设得越小，越无法利用重叠矩阵和KS矩阵的稀疏性节约时间，但这对耗时的影响相对于杂化泛函计算的总耗时来说微不足道。目前版本的Multiwfn产生的杂化泛函的输入文件里自动就把EPS_PGF_ORB设成了1E-12。

对于晶胞较小的情况，本身矩阵稀疏性就差，开发者建议上来就在&QS里设MIN_PAIR_LIST_RADIUS -1从而完全不利用矩阵的稀疏性节约时间，此时绝对不会再出现那个警告，但比起用EPS_PGF_ORB 1E-12时耗时会高一些。

8 简单例子：C60分子晶体

最后，这里以C60分子晶体的单点任务作为例子演示一个标准的周期性杂化泛函的计算。如果你还不了解Multiwfn创建CP2K输入文件的功能，建议先阅读《CP2K第一性原理程序在CentOS中的简易安装方法》（//www.umsyar.com/586）。

此体系的cif文件和下面涉及的输入文件都可以在//www.umsyar.com/attach/690/file.zip中找到。此体系结构如下所示，一共240个C原子。晶胞边长是14.07埃，对于分子晶体来说不扩胞的情况下只考虑gamma点可以接受。

启动Multiwfn，输入C60.cif的路径载入之，然后输入
cp2k
PBE.inp   //要产生的PBE计算的输入文件名
0   //基于默认设置产生输入文件
现在当前目录下就有了PBE.inp，对应PBE/DZVP-MOLOPT-SR-GTH单点任务。然后接着在Multiwfn里输入
cp2k
PBE0.inp   //要产生的PBE0计算的输入文件名
1   //选择理论方法
-6   //PBE0结合ADMM
0   //产生输入文件
由PBE0.inp内容可见，Multiwfn自动指定的计算设置完全合适，对应的是PBE0/DZVP-MOLOPT-SR-GTH结合admm-dzp辅助基组、6埃库仑截断、EPS_SCHWARZ 1E-6、SCREEN_ON_INITIAL_P T、开启OT。

手动把PBE0.inp中的WFN_RESTART_FILE_NAME前头的#去掉并设为PBE-RESTART.wfn，把SCF_GUESS RESTART开头的#去掉，从而使得PBE0计算时读取PBE收敛的波函数当初猜。然后根据你的机子的空余物理内存量和要用的MPI进程数设置合适的MAX_MEMORY。

把PBE.inp和PBE0.inp放到同一个目录下，先运行PBE.inp，算完后当前目录下就出现了PBE-RESTART.wfn，然后再运行PBE0.inp。在《淘宝上购买的双路EPYC 7R32 96核服务器的使用感受和杂谈》（//www.umsyar.com/653）介绍的笔者的2*7R32 512GB双路服务器上用96核并行、MAX_MEMORY设5000，用CP2K 2023.2 popt版，PBE计算花了8秒，PBE0计算花了178秒。PBE0任务的SCF迭代过程如下所示，可见第一轮耗时较高，由于当前给的内存足够储存所有ERI，因而ERI没有on-the-fly计算的，所以之后每一轮SCF的耗时远远低于第一轮。

9 总结

本文把强大的CP2K高效地做杂化泛函计算的各方面要点进行了介绍，以避免初学者对各种关键信息毫不知情而凭感觉胡用瞎用。本文的内容也充分体现出用Multiwfn产生CP2K输入文件非常容易，令研究者从复杂的输入文件编写中充分解放，但这绝不意味着CP2K的使用者就可以对关键性技术细节毫不知情，否则会各种踩坑、白浪费时间、算出无意义的结果。CP2K用户必须掌握的知识、要领远多于大部分其它计算程序。笔者开设的北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）可以使得学员完整系统透彻掌握全部这些知识。

Simulating scanning tunneling microscope (STM) images of periodic systems using Multiwfn in combination with wavefunction produced by CP2K

文/Sobereva@北京科音    2023-Jun-4

1 前言

强大的波函数分析程序Multiwfn（//www.umsyar.com/multiwfn）有非常方便好用的模拟隧道扫描显微镜（STM）图像的功能，模拟的原理以及Multiwfn中此功能的用法我在《使用Multiwfn模拟扫描隧道显微镜(STM)图像》（//www.umsyar.com/549）里有专门说明，但那篇博文是以孤立体系为例。Multiwfn支持对CP2K产生的周期性体系波函数进行大量的分析，在北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）里我做了很多介绍，其中包括模拟周期性体系的STM图。本文的目的是举一个完整的例子，能令甚至之前没用过CP2K的人都能轻易使用Multiwfn模拟周期性体系的STM图。没读过//www.umsyar.com/549的读者必须要先仔细阅读此文。不了解Multiwfn的话建议看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。值得一提的是，CP2K自己也有模拟STM的功能，但远没有Multiwfn来得好用和方便，所以这里就不多说了。

2D Mater., 6, 015005 (2019)文中给出了几层厚度的黑磷薄片表面的0.7 V偏压下的常电流模式的STM图，如下所示，本文将对黑磷计算同样条件下的STM图并与之对照。

本文为了节约时间只用单层黑磷，而且图省事直接从三维晶体中截出来一层结构直接做单点计算得到波函数。严格来说应该用多层模型，并且固定最下层黑磷并令上层的自发弛豫。

本文的例子涉及的文件都可以在//www.umsyar.com/attach/671/file.rar中得到。读者务必使用2023-Jun-4及以后更新的Multiwfn版本，否则情况与本文所述不同。

2 得到单层黑磷的molden文件

如《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）所述，用Multiwfn分析CP2K的波函数必须让CP2K产生记录波函数信息的molden文件。这一节我们就对单层黑磷做一个单点计算来得到此文件。

本文文件包里的Phosphorus-black.cif是黑磷的晶体结构，用GaussView打开，只保留中间的单层部分而删除其它原子，如下所示，然后保存为Phosphorus-black.gjf。

Multiwfn有极其便利的创建CP2K输入文件的功能，见《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587），这里我们用Multiwfn创建产生molden文件的单点任务文件。由于CP2K没法产生考虑k点时的molden文件，因此必须构造超胞用gamma点计算。当前的晶胞参数是a=3.31埃、b=4.38埃，在a方向复制为7倍、b方向复制为5倍后两个方向都有>20埃，此时只考虑gamma点是可以接受的。当前模拟STM用的偏压为正，此时能量最低的一个或多个空轨道对STM会产生贡献，因此必须要求CP2K求解出一些空轨道才行，几十个就够了。模拟STM用的偏压越大应当计算的空轨道越多，因为越高的空轨道会被涉及。

启动Multiwfn，载入Phosphorus-black.gjf，然后输入
cp2k   //产生CP2K输入文件
SP.inp  //将要产生的CP2K输入文件名
-7  //设置周期性
XY  //XY二维周期性
-2  //产生molden文件
-11  //进入结构编辑界面
19  //产生超胞
7  //a方向复制的倍数
5  //b方向复制的倍数
1  //c方向（垂直于表面的方向）不变
-10  //返回
-9  //其它设置
12   //计算空轨道
20  //算最低20个空轨道
0  //返回
0  //产生输入文件

此时SP.inp就产生在当前目录下了，在本文的文件包里提供了，计算级别是Multiwfn默认的PBE/DZVP-MOLOPT-SR-GTH。用CP2K运行SP.inp，在我的2*EPYC 7R32双路服务器上不到10秒钟就算完了。之后按照《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）所述，把晶胞信息和价电子信息手动加入.molden文件，即在开头插入以下字段。本文文件包里的SP-MOS-1_0.molden是已经改好的。

 [Cell]
 23.17000000     0.00000000     0.00000000
  0.00000000    21.90000000     0.00000000
  0.00000000     0.00000000    12.11600000
 [Nval]
 P 5

3 用Multiwfn模拟黑磷常距离模式的STM

首先模拟一下常距离模式的STM。启动Multiwfn，载入SP-MOS-1_0.molden，然后输入
300  //其它功能 Part 3
4  //模拟STM
2  //设置偏压
0.7  //0.7V，和文献里的相同
0   //开始计算

默认计算的是Z最大值原子上方0.7埃的XY平面，X和Y方向计算范围和当前晶胞的X和Y方向的跨度一致。瞬间就算完了，在新出现的作图设置菜单里输入0绘图，图像如下所示

可以再做一些调整，只让最上层的P用粗体字标注，下层的用细体字显示以作区分，并且用和文献里相似的配色，并且调整坐标轴让刻度整齐。为此，把图像关闭后，接着输入
4  //修改显示原子标签的距离阈值
0.71 A //在图上显示距离作图平面0.71埃以内的原子的标签
y  //将距离作图平面更远的原子的标签以细体字显示
9  //修改色彩变化方式
13   //黑-橙-黄
7  //修改标签间隔
3,3,0.0005
8  //修改色彩刻度
0,0.002  //下限和上限
1  //保存图像

此时的图像如下所示

此图的效果已经很好了。和实验STM图像展现出的信息近似一致。

4 用Multiwfn模拟黑磷常电流模式的STM

这次模拟常电流模式的STM，这与文献里的STM图对应。启动Multiwfn，载入SP-MOS-1_0.molden，然后输入
300  //其它功能 Part 3
4  //模拟STM
2  //设置偏压
0.7  //0.7V，和文献里的相同
1  //切换模式为常电流模式
0   //开始计算
3   //绘制STM图像
0.0009   //常电流值（反复尝试，取一个图像效果较好、和实验图像较接近的即可）
4  //修改显示原子标签的距离阈值
2.6 A  //在图上显示距离被计算区域顶端2.6埃以内的原子的标签（计算的区域顶端的z坐标是计算开始之前界面里7 Set range in Z direction选项后面显示的第二个值。当前设2.6埃可以只把顶层磷原子纳入进来）
y  //将距离更远的原子的标签以细体字显示
9  //修改色彩变化方式
13   //黑-橙-黄
7  //修改标签间隔
3,3,0.1
1  //保存图像
现在看到下图，和常高模式的STM图展现出的信息基本一致

5 其它

本文的例子充分体现出用Multiwfn+CP2K绘制周期性体系表面的STM真是又快又方便，鼓励大家应用在实际当中。请别忘了恰当引用Multiwfn。作为练习，大家可以绘制一下单层MoS2和石墨炔的STM图。

对于无gap的体系，如金属表面（包括吸附小分子后），绘制STM有一点需要注意。这种情况一般都是开smearing的，此时费米能级附近的轨道占据数不为整数，导出的molden文件里的轨道占据数也是如此，此时Multiwfn没法直接绘制STM（Multiwfn会以为是记录自然轨道的多组态波函数）。解决方法是先进入主功能300的子功能9（计算费米能级功能），此功能会将电子按照轨道能量由低到高以整数方式填充来修改轨道占据数，相当于0 K的占据方式。此外，在此功能里输入当前的温度，程序会给你费米能级值，将之记下来。退出此功能后，就可以照常进入STM绘制界面，届时再把费米能级值输入进选项3 Set Fermi level即可。另外值得一提的是，用Multiwfn创建CP2K输入文件时开启smearing的话，产生的输入文件里自动就会有计算空轨道的设置，就不用自己再单独设置了。

Solution when SCF is difficult to converge in CP2K

文/Sobereva@北京科音

First release: 2023-May-15   Last update: 2024-Feb-22

0 前言

经常有人问CP2K中遇到SCF难收敛怎么解决，这是很家常便饭的问题。这里就基于笔者讲授的“北京科音CP2K第一性原理培训班”（http://www.keinsci.com/workshop/KFP_content.html）里面相关部分的内容专门写个博文，供CP2K用户参考。这里假定读者已经了解CP2K计算的基本的背景知识（如k点、CUTOFF、对角化、OT、&BS、smearing等等），如果不了解相关名词、算法的话，十分推荐参加这个培训全面、系统学习一遍。

笔者之前专门写过Gaussian用户遇到SCF不收敛的解决方法的文章《解决SCF不收敛问题的方法》（//www.umsyar.com/61），本文和那篇文章一样，已经把所有可能有用的方法和要注意的问题说得尽可能全面、充分了，读者遇到难收敛/不收敛时需要结合实际情况一个个尝试。切勿问笔者“这文章我看了，但还是解决不了”之类的，我没有任何可额外补充的（至多是给出当前计算的尽可能详细的信息，比如给出输入输出文件，然后我能回答哪些问题需要优先考虑、哪些方法值得优先尝试）。

下面第1节先说说遇到SCF不收敛时一般情况需要考虑的问题和解决策略，然后第2、3节再分别说使用对角化和OT时分别可以进一步考虑的。最后再说一下几何优化、分子动力学时与SCF收敛有关的问题。在此之前先看培训班ppt了解一下CP2K是怎么判断SCF收敛的：

1 一般情况时考虑的

• 检查结构合理性。SCF是否容易收敛和结构的合理性关系极大，化学意义越强的结构通常越容易收敛。因此遇到SCF很难收敛时，需要结合结构化学、固体物理常识检查结构是否靠谱，如模型是否有意义、是否漏了原子、是否某些原子所处的化学环境和实际明显不符（如过渡金属本该与配体成键的地方却光秃秃地悬着）等等。且要结合图像肉眼检查，确保不存在不该出现的结构严重扭曲、盒子边缘原子和周期镜像原子发生过近接触等问题。若有这些问题，即便SCF收敛了，结果显然也毫无意义。

• 确认关键词、设置没明显不合理的。CP2K的关键词很复杂，对使用者的知识水平要求很高，若用了明显不当的关键词和设置，SCF很难收敛或收敛到离谱结果是极其正常的事，然而很多要点是初学者难以注意到的。比如计算自旋极化体系的时候必须得在&DFT里写UKS（或等价的LSD），否则即便自旋多重度设的是>1，程序也会愣当成是闭壳层算。再比如用wavelet的Poisson solver算孤立体系的时候，如果体系有明显电子密度分布的区域跨越了盒子，结果就无意义。

CP2K属于那种是完全手动挡还很个性的车。如果经验不够丰富的话，笔者建议总是用Multiwfn（//www.umsyar.com/multiwfn）产生其输入文件，见《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）。Multiwfn的CP2K输入文件创建功能将CP2K尽最大可能包装成黑箱，不仅用起来极其方便省事，对初学者还能大大减小关键词使用不当的可能性。比如将自旋多重度设成>1时自动就会添加UKS、用wavelet做非周期性计算时会自动设置足够大又不太浪费的盒子尺寸并且加入自动让体系在盒子里居中的设置。

• 确保体系的净电荷、自旋多重度的设置合理。这俩如果设置不当的话，相当于计算了无意义或激发态的电子结构，SCF极难收敛也是很可能出现的事。特别是对于磁性体系，自旋多重度的设置对SCF收敛往往影响极大。

• 如果SCF有收敛趋势，可增大SCF迭代次数上限再试，即&SCF里的MAX_SCF，其默认的50对于很多情况明显太小。Multiwfn产生的CP2K输入文件里默认用128（128是为了优雅，2的7次方。和Gaussian程序一致），绝大多数情况都够了，但也可能碰上极个别体系仍需要加大。

笔者超级反感一碰到SCF不收敛，也不观察收敛趋势，上来就增加迭代次数上限的做法，因为这是最典型的菜鸟思维，绝大多数情况毫无用处而白浪费时间。因为如果SCF过程都发生显著震荡了、没收敛的苗头，再怎么迭代下去都是白搭，纯粹是毫无意义地浪费时间和计算资源。

• 如果发现用其它泛函或基组下计算能收敛的话，可以读取其收敛的波函数作为初猜，这至少比自动产生的初猜波函数要好、更容易收敛。也即&SCF里用SCF_GUESS RESTART，且&DFT里WFN_RESTART_FILE_NAME指定成之前任务的.wfn文件（若考虑了k点是.kp文件）。注意两次的基组特征必须定性一致，比如之前计算用的大核赝势，而现在用小核赝势或全电子基组，这时候读之前的波函数就没意义了。

通常较小的基组下计算比大基组下计算更容易收敛，比如直接用TZV2P档次的MOLOPT基组不收敛的话，可以尝试用DZVP档次的算，若收敛了，读它收敛的波函数为初猜。

往往HF成份越大的泛函越容易收敛，很大程度是因为其HOMO-LUMO gap更大、SCF过程中占据和空轨道混合程度相对更缓和所致。因此遇上SCF难收敛时可以尝试用HF成份比实际要用的泛函更大的泛函，若能收敛，读取其收敛的波函数当初猜。HF成份可以参考《不同DFT泛函的HF成份一览》（//www.umsyar.com/282）。

如果原本用的是纯泛函，想尝试杂化泛函但用不动，也可以用DFT+U试试，耗时和纯泛函基本一致，但此时不支持k点，此外怎么获得+U用的Ueff参数也是个问题。如果实际计算不打算+U，也可以尝试3或4 eV的+U看看能否令SCF收敛，如果能，则不+U时读取其收敛的波函数当初猜。

• 做杂化泛函计算时（特别是对于周期性体系）通常会在&DFT / &XC / &HF / &SCREENING里开启SCREEN_ON_INITIAL_P从而基于初始的密度矩阵做双电子积分的屏蔽来巨幅节约耗时，Multiwfn产生的杂化泛函计算的输入文件里默认也开了这个。若此时没有读取之前收敛的波函数当初猜（通常是便宜得多的纯泛函计算的），则CP2K会基于自动初猜的波函数对应的密度矩阵做积分屏蔽，此时结果通常无意义，SCF能收敛的概率也非常低！然而这个杂化泛函计算的关键要点却经常被初学者所忽略！笔者强烈建议在做杂化泛函计算之前仔细把此文认真读一遍：《CP2K做杂化泛函计算的关键要点和简单例子》（//www.umsyar.com/690）。

• 确保恰当考虑了k点、k点数够多。如果晶胞边长不够大，显然必须考虑k点而不能只用默认的Gamma点。而且k点设得还得够多，否则SCF难收敛极为正常。k点如果设置不当，就算SCF能收敛了也没意义。众所周知，算能量问题，原则上k点应当大到足以令能量随k点数增加而基本收敛的程度。

• 尝试增大&DFT / &MGRID里的CUTOFF、REL_CUTOFF。这俩截断能设置控制格点精度，不够大的话不仅结果不准，甚至可能造成SCF难收敛。如果当前明显已经够大了，就没必要再盲目尝试再增加了。注：还有种说法是尝试&DFT / &XC / &XC_GRID里写USE_FINER_GRID T将计算交换-相关泛函部分的格点精度提升为CUTOFF的四倍，但根据我的经验这并没有什么实际用处，而耗时暴增。

• CP2K默认的产生初猜波函数的方法是孤立原子密度的叠加，在整个任务一开始时会对每个KIND对应的原子在孤立状态下做计算，用PRINT_LEVEL MEDIUM时可以看到具体细节。计算孤立原子所用的电子组态可以通过&KIND里的MAGNETIZATION（定义alpha电子数减beta电子数）或者&BS（可以具体控制各个主量子数的各个角量子数上的alpha和beta电子数）来设置。初猜波函数中的孤立原子的组态越接近整个体系收敛的波函数中原子的实际组态则SCF越容易收敛、SCF不收敛的可能性越低。

如果体系有磁性，遇到SCF不收敛时一定别忘了尝试用MAGNETIZATION或&BS尽可能令初猜波函数中原子组态与实际情况接近，这对SCF是否能够顺利收敛有关键性影响。如果当前体系在文献（可用google搜索来找）或高通量数据库（如Material Project）中已经被计算过，可以参考那里面报道的原子自旋磁矩或原子组态来设。有时同一种元素在体系里处于不止一种化学环境，差异很大时（比如明显处于不同价态）可分成多个&KIND分别设置。

即便体系没有磁性，利用&BS使初猜波函数中的原子带电状态使之尽可能符合实际也往往是有益的。比如对于NaCl这样高度离子性体系计算时SCF不收敛的话，可以尝试令初猜中Na和Cl的原子组态分别对应Na+和Cl-的状态（默认情况下分别对应中性的Na和Cl）。曾经我用PBEsol/pob-TZVP-rev2计算10层的NaCl板时通过这种做法成功解决了SCF不收敛问题。

• 尝试将&QS / EPS_DEFAULT减小几个数量级，默认为1E-10。EPS_DEFAULT这个参数控制一堆涉及数值计算精度的参数，如果设得太大，有可能造成SCF难收敛。若这个数值当前已经很小了，比如都1E-12了，尝试继续减小则没必要。

• 改用其它泛函、基组。如果原先的计算级别就是极难收敛，而有其它级别性价比、精度相仿佛，且发现能收敛，那就别死磕当前的了。

• 在GPW与GAPW之间切换。如果之前用的GPW结合GTH赝势基组不收敛，尝试用GPAW结合一个全电子基组计算试试，反之亦然。

• 修改默认产生波函数的方式。参见手册里&SCF / SCF_GUESS设置。但其实这没什么用，除了默认的ATOMIC（即前述的基于孤立原子密度叠加产生的）以及RESTART（读取之前的波函数当初猜），其它的一个能打的都没有。

用ATOMIC产生初猜时，CP2K计算每个KIND对应的孤立原子也是需要SCF迭代的，上限是固定的50轮（没法靠关键词改），在PRINT_LEVEL MEDIUM时其迭代过程会在一开始输出。有时候我发现到了50轮时收敛精度还很差，故有可能在某些情况下由于这个问题导致初猜质量烂，不利于实际体系的SCF收敛。我发现可以通过atom_kind_orbitals.F里atom%optimization%max_iter修改迭代次数上限，比如改成300然后重新编译，一般孤立原子的SCF就都能收敛了。

• 如果加了外电场时SCF很难收敛，可能是外电场过大所致，可尝试更小的数值。如果就是要用较大的外电场，可以试试读取较小电场或不带电场时收敛的波函数当初猜。有的体系在很强外电场下本来就不可能稳定存在，这时候SCF不收敛自然很正常，是物理上决定的，也别盲目去计算。

• 实在没辙时干脆放宽SCF收敛限，即加大&SCF / EPS_SCF，默认是1E-5（这其实已经挺宽松了），这是走投无路时的做法。或者换程序，如免费且和CP2K同样流行的Quantum ESPRESSO，或者尝试CP2K里的纯粹基于平面波的SIRIUS模块。某些体系用CP2K的quickstep模块就是很难收敛，如某些d族纯金属的slab，这和算法、基函数特征有关。

2 用对角化的情况可专门考虑的

• 对无gap（如金属、石墨烯）或gap很小的体系（非常窄带半导体）应当使用smearing。即加入&SCF / &SMEAR / METHOD FERMI_DIRAC，且在&SCF里用ADDED_MOS要求求解一些空轨道从而能被热激发的电子所占据。ADDED_MOS数目够用就行，设太大了会增加耗时而且也没用处。Multiwfn产生CP2K输入文件时如果开了smearing会自动做ADDED_MOS设置，并根据体系原子数估计一个通常够用的值。用较小电子温度（如200K、300K）时若不收敛可尝试再提升温度，即&SMEAR / ELECTRONIC_TEMPERATURE设的。温度设高时ADDED_MOS往往需要相应地增加，否则热激发的电子没法按Fermi-Dirac分布占该占的空轨道（此时计算时会有警告）。

我老在网上看到有人不理解smearing原理是什么就胡用、乱试smearing以试图解决SCF不收敛，所以这里强调一下。一定要注意，对于gap不是很小的体系（比如氧化锂等绝缘体），根本没必要尝试smearing，本来这在原理上就不是smearing能起到帮助SCF收敛的情况，而且本文还有其它一堆方法可以尝试用于帮助SCF收敛，完全轮不到考虑用smearing。

用smearing还会对结果产生一定影响，结果会依赖于电子温度T的设置，用Fermi-Dirac函数做smearing时最终给出的能量相当于是T温度下电子的自由能。如果你就是要得到基态（0K）的电子能量，原则上最好做T→0的外推。如果你不做外推，在能量横向对比时至少也应当用相同的smearing设置。关于smearing的物理意义和实际意义，以及做外推的方法，我在北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）里有专门详细的讲解。

• 对gap较小的体系可尝试使用能级移动，即&SCF / LEVEL_SHIFT，设比如0.3、0.4左右，默认单位是Hartree。这使得SCF过程中人为加大空轨道与占据轨道之间的能级差，从而减轻彼此间混合的剧烈程度而有时令SCF更容易收敛。

• 修改&SCF / &MIXING / METHOD设置SCF过程中将旧密度矩阵与新密度矩阵混合得到实际每一轮密度矩阵的方法。默认的DIRECT_P_MIXING往往表现极差。通常用BROYDEN_MIXING，这也是Multiwfn产生的输入文件里默认的。亦可尝试PULAY_MIXING，但一般不比BROYDEN_MIXING好。

• 对于使用BROYDEN_MIXING的情况，尝试修改&MIXING中的ALPHA（新密度矩阵混入旧的的比例，一般用默认的0.4，即混入40%。难收敛时可尝试诸如0.1、0.2、0.3）和NBROYDEN（默认的4明显太小，改大到8通常很有助于收敛，这是Multiwfn产生的输入文件默认的。也可以尝试进一步改大到比如12）

• 对于使用PULAY_MIXING的情况，尝试修改&MIXING中的NPULAY （默认的4偏小，改大往往很有助于收敛）。注：NPULAY、NBROYDEN、NBUFFER三个关键词实际上等价。

• 对于使用DIRECT_P_MIXING的情况，收敛到&SCF / EPS_DIIS值后会切换到加速收敛的DIIS算法。但默认的EPS_DIIS为0.1，有点偏小，可以尝试改大以早点启用DIIS（但借助DIIS通常仍不如BROYDEN_MIXING）。还可以把控制最大DIIS矢量的&SCF / MAX_DIIS从默认的4改大到诸如7。

• 在对角化框架下各种方法都不行的话，尝试改用OT。如果OT能收敛的话，而且当前由于特殊目的又必须用对角化，也可以用OT收敛的波函数当对角化计算的初猜。由于OT不支持k点，因此如果之前对角化是结合k点算的较小晶胞情况，用OT时就只能靠超胞模型了。OT也不支持用ADDED_MOS设置在SCF计算过程中求解空轨道，因此必须用smearing的导体体系就别指望OT了。如果是用GFN1-xTB，总应当用OT，这也是Multiwfn默认的。

3 用OT的情况可专门考虑的

• 如果之前没用outer SCF则开启之，这对解决SCF不收敛非常有效，也即加入&SCF / &OUTER_SCF字段。由于其重要性，Multiwfn产生的输入文件里开OT后默认就是启用outer SCF的。此时每次到达&SCF / MAX_SCF所设的inner SCF的步数上限时，就会进入outer SCF过程更新OT的preconditioner，然后进入下一次的inner SCF迭代过程。每次更新preconditioner后经常会看到比之前以更快的速度趋近收敛限。原理上来说，OT每一步都更新preconditioner是最理想、SCF最容易收敛的，但由于这样太昂贵，因此默认情况下（即不用outer SCF时）只会在最开始计算一次preconditioner。

开启outer SCF时inner SCF的迭代次数上限应当设得比平时更小（否则经历很多次SCF迭代才做一次outer SCF，无法有效达到用outer SCF的目的），15至35是比较合适的范围。outer SCF迭代次数上限通过&SCF / &OUTER_SCF / MAX_SCF控制，总SCF迭代次数上限相当于inner SCF次数上限和(outer SCF次数上限+1)的乘积。

&SCF / EPS_SCF设的是inner SCF收敛限，&SCF / &OUTER_SCF / EPS_SCF设的是outer SCF收敛限，后者应当小于或等于前者（一般设为相同）。

• 尝试其它的preconditioner。FULL_ALL通常最好，但最贵。大体系可以尝试明显更便宜的FULL_SINGLE_INVERSE和FULL_KINETIC。个别时候FULL_KINETIC反倒比FULL_ALL更容易收敛。GAPW计算总应当用FULL_ALL。

• 修改OT算法，把&OT里的ALGORITHM从默认的STRICT改为IRAC，SCF经常会收敛得更稳健，非常值得尝试！特别是我发现用pob系列全电子基组的时候，用IRAC往往容易收敛得多，因此Multiwfn产生基于pob系列基组的输入文件时若用户选了OT，则默认用IRAC。而且用约束性DFT（CDFT）时，我发现用IRAC也会令SCF容易收敛得多得多。

• &OT / MINIMIZER控制OT用的极小化方法，若之前用的是较快的DIIS（也是Multiwfn产生的用OT的输入文件里默认的），可改为更稳健但耗时更高的CG（共轭梯度法），也可以尝试BROYDEN。

• MINIMIZER用CG时，尝试用&OT / LINESEARCH选项把线搜索方式设为比默认2PNT更贵但也更稳的3PNT。

4 几何优化和分子动力学过程中的SCF收敛问题

在几何优化过程中，有可能中途偶尔出现几步SCF不收敛。此时程序还会继续算下去。只要最终几何优化收敛时SCF也是收敛的状态就行了。（注：从CP2K 2024.1开始，优化中途SCF不收敛会导致任务自动停止，&SCF里写IGNORE_CONVERGENCE_FAILURE才会在SCF不收敛的情况下也继续优化）

有时候在几何优化初期SCF一直不收敛，不要急于把任务停掉，可以先凑合跑跑看看情况。初始几何结构相对于实际的极小点结构而言肯定是存在多多少少不合理性的，此时SCF也肯定比极小点结构处更难收敛。在SCF不收敛、受力也较糙的情况下凑合跑一定步数后，结构一般会变得比一开始更为合理，此时SCF也往往变得更容易收敛。

在几何优化、分子动力学过程中，程序会自动利用之前一些结构的波函数外推出当前步的初猜波函数，这比起直接产生的初猜波函数质量通常好得多。外推方式可以通过&DFT / &QS / EXTRAPOLATION来具体设置，默认的ASPC通常是最好的选择。由于这一点，当几何优化、分子动力学跑一定步数之后，SCF达到收敛限所需的迭代次数明显减小、出现SCF不收敛的可能性也大为降低。

注意当使用k点时CP2K不支持外推产生初猜波函数，因此会在几何优化、分子动力学等涉及结构变化的任务中每一步重新产生初猜波函数，即EXTRAPOLATION USE_GUESS的情况。我一般建议改为EXTRAPOLATION USE_PREV_P，代表直接使用上一步SCF收敛的波函数当初猜。

如果想减小MD中途出现SCF不收敛的可能，可以尝试减小步长（&MD / TIMESTEP）。如果想减小几何优化中途出现SCF不收敛的可能，可减小优化的步长上限，即置信半径（如通过&MOTION / &GEO_OPT / &BFGS / TRUST_RADIUS）。这样做之后由于相邻两个结构间差异变小了，ASPC或USE_PREV_P给出的初始波函数对于当前结构而言就越理想，自然SCF收敛也会更容易，且有利于让SCF能收敛的状态保持下去。

sobNEB: A convenient tool of generating interpolation points for NEB task of CP2K

文/Sobereva@北京科音  2023-Feb-28

1 介绍

CP2K有很好用的NEB功能可以得到反应路径，靠CI-NEB还可以同时得到过渡态结构，用法在北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）里讲得很详细。虽然CP2K的NEB功能自己就能自动线性插点，但是偶有bug（有的初始间距是其它的两倍），而且没法人为控制自己提供的中间结构对应哪个点，也不便于预览初始点。

我开发的sobNEB在插点时可控性更强。虽然也是线性插点，但有以下好处
(1)可以直接控制自己提供的结构对应哪个初始点，其它点都在相邻的提供的结构之间线性插入
(2)可直接产生对应CP2K的&BAND部分中的信息便于写NEB输入文件
(3)会将用户提供的点和所有额外插的点的结构一起写入traj.xyz文件以在VMD中观看所有初始点构成的轨迹，以便判断插点是否合适。

sobNEB的下载地址：//www.umsyar.com/soft/sobNEB_1.0.zip。有.exe后缀的是Windows版，无后缀的是Linux版。

注意sobNEB不会对输入的结构自动消除平动转动。

2 用法

用户只需编辑sobNEB目录下的sobNEB.ini，指定自行提供的各个xyz文件的路径以及对应的点的序号（之间以冒号分隔），在启动程序后就会开始处理。只有始、末端的结构是必须自己提供的，中间的结构可以不提供也可以提供一个或者多个。

下面是sobNEB.ini文件的例子，涉及到的三个xyz文件在sobNEB目录下都提供了。conf1.xyz和conf2.xyz是分子的两个构象，分别作为始、末端。TS_guess.xyz是自己摆的过渡态初猜结构，其点号位于始末、端点号的正中央。

conf1.xyz : 1
TS_guess.xyz  : 11
conf2.xyz : 21

启动sobNEB后会看到以下信息

 sobNEB: Generate interpolated .xyz files and CP2K input for NEB calculation
 Programmed by Tian Lu (sobereva[at]sina.com)
 Version 1.0, release date: 2023-Feb-19

 Loading sobNEB.ini...
  3 structures are given

 Generating interpolated structures between given systems  1 and  2
 Generating point   2 : 2.xyz
 Generating point   3 : 3.xyz
 Generating point   4 : 4.xyz
 Generating point   5 : 5.xyz
 Generating point   6 : 6.xyz
 Generating point   7 : 7.xyz
 Generating point   8 : 8.xyz
 Generating point   9 : 9.xyz
 Generating point  10 : 10.xyz

 Generating interpolated structures between given systems  2 and  3
 Generating point  12 : 12.xyz
 Generating point  13 : 13.xyz
 Generating point  14 : 14.xyz
 Generating point  15 : 15.xyz
 Generating point  16 : 16.xyz
 Generating point  17 : 17.xyz
 Generating point  18 : 18.xyz
 Generating point  19 : 19.xyz
 Generating point  20 : 20.xyz

 Done! .xyz file of each point has been generated in current folder
 Corresponding &BAND field of NEB task of CP2K has been exported to BAND.txt in
current folder
 traj.xyz is the trajectory file containing all points for previewing purpose

可见产生了2-10号、12-20号，共18个插点结构，当前目录下出现了相应序号的xyz文件。

当前目录下出现的BAND.txt内容如下，可以直接粘贴到&BAND字段里
    &REPLICA
      COORD_FILE_NAME conf1.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 2.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 3.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 4.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 5.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 6.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 7.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 8.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 9.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 10.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME TS_guess.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 12.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 13.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 14.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 15.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 16.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 17.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 18.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 19.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 20.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME conf2.xyz
    &END REPLICA

可以用VMD检查当前目录下同时产生的traj.xyz。里面21帧同时叠加显示时如下所示，可见初始点的分布是合理的，在两个构象间变化，可以用于NEB计算。