使用Multiwfn结合CP2K计算晶体中原子的氧化态

使用Multiwfn结合CP2K计算晶体中原子的氧化态
Using Multiwfn in combination with CP2K to calculate oxidation state for atoms in crystals

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

1 前言

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

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

2 用法

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

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

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

3 例子

下面给出一系列通过Multiwfn使用LOBA/mLOBA计算固体的氧化态的例子。CP2K的输入输出文件在http://sobereva.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格式的波函数文件》（http://sobereva.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的结果对阈值敏感），比mLOBA在原理上明显更好。在当前界面里输入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观看分子轨道》（http://sobereva.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进行正确的引用。