使用QE计算ELF和DORI

欧拉角

本文原创非首发，未经同意不得转载内容。

引言

原子、离子或分子间结合成分子或晶体，其结合形式与物质的结构和物理化学性质有密切的联系[1]。共价键是通过原子间共享电子而形成，这些电子呈现出较强的定域性。那么什么是电子的定域性呢？电子的定域性并非指电子静止在某个空间点，而是指电子在某个区域内不断地运动，如共价键的电子通常就是在两个原子之间巡游[2]。为了定量展现电子的这种定域特性，Savin等提出了电子定域化函数 (electron localization function, ELF)[3]：

其中，D为精确动能密度与Weizsäcker动能密度之差，D0为无相互作用电子的动能密度。ELF为无量纲物理量，其数值范围由0到1，数值越大意味着空间中某处的电子的定域性越强。其中ELF等于1对应于电子的完全定域化；等于0则表示电子完全变为均匀电子气。我们往往关注的是ELF数值大于0.5的区域，其对应于成键区域、孤对电子、原子的壳层和内核。

另一种讨论原子间相互作用的方式是由Silva和Corminboeuf提出的密度重叠区域指示符(density overlap regions indicator, DORI)[4]：

DORI能够反映出原子间电子密度重叠的程度，它不仅能体现共价相互作用，同时也能将非共价相互作用区域呈现出来。DORI可以通过类似于上一篇文章中[5]所介绍的RDG那样进行sign(λ2)ρ填色，与ELF互补来进一步讨论所研究体系的相互作用类型。为了对上述两种函数在分析共价相互作用方面的应用进行说明，本文将展示若干Quantum ESPRESSO的计算实例。

准备工作

计算程序采用Quantum ESPRESSO 6.7[6]和VESTA 3.1.8[7]
测试体系为：
1. 来源于Materials Project [8]的结构GaAs(mp-2534)、NiO(mp-19009)
2. H2O2@PMCS来源于文献[9]
3. 在进行下一步计算前建议先阅读我的上一篇文章[5]，一些类似的操作本文将不再赘述。

计算参数

交换相关泛函	PBE
赝势	PSLibrary PAW
动能截断	60 Ry
密度截断	480 Ry
K网格	GaAs 8×8×8 NiO 10×10×10 H2O2@PMCS 2×2×1
展宽	0.02 Ry

计算过程和数据处理

A 计算过程
A1~A4计算过程可以参考[5]，只是在A3的过程中需要通过plot_num=8产生ELF格点，以及我们修改过的pp.x程序plot_num=123[注1]产生DORI格点。
B DORI填色图处理
B1~B4也是与[5]相同的，只是将B2过程中导入的RDG格点数据换成DORI格点数据，并且将DORI的等值面取为0.95。
C ELF切片填色图处理
C1 打开VESTA程序，将A2中out文件产生的xsf文件拖入VESTA窗体。

C2 选中我们关心的平面上的三个原子。

C3 在菜单栏中Edit->Lattice Planes->Add Lattice Planes->New建立新的切片，在New按钮的上方的Calculate the best plane for the selected atoms使得表面尽可能贴合选中的原子。

C4 在菜单栏中的Objects->Properties->Sections->Sections and slices->Saturation levels中的Min和Max分别改为0和1。

C5 调整其他显示部分，如原子的颜色、连接等。

案例

GaAs是一种原子晶体，通常作为太阳能电池材料，左下角为其晶体结构和ELF切片填色图。从图中我们可以发现，对于Ga-As之间的区域，电子是具有定域性的；右下角的图为NiO的晶体结构和ELF切片图，图中则体现出Ni-O之间的区域电子不具备定域性，因此我们可以认为GaAs中的Ga-As形成了共价键，而NiO中的Ni-O并不是共价键。

H2O2@PMCS是我们上一篇文章[5]所讨论过的例子，下图为其ELF切片填色图，切片分别选取了过氧化氢的两个氧原子和锌原子以及PMCS上的三个碳原子构成的两个平面。图中的PMCS上的C-N、C-C原子间和H2O2的O-O、O-H之间的区域都呈现较强的定域性，表明这些原子是通过共价作用结合的。H2O2的O原子不与氢连接的区域也显示出定域性，这是由于O原子带有孤对电子。Zn原子与其他原子间的区域为蓝色，意味着Zn原子与其他原子并非通过共价作用结合，但是无法区分这些作用是离子键还是弱相互作用。

我们通过DORI对H2O2@PMCS做进一步分析。下图中能够得到与RDG类似的结论，H2O2的O原子与PMCS上的Zn原子为静电作用，H2O2的O-H与Zn-N4的N原子为氢键作用。但不同于RDG，PMCS上的C-N、C-C原子间和H2O2的O-O、O-H之间的共价作用能够同时清晰地被展现出来，Zn-N之间的sign(λ2)ρ已经小于-0.1[注2]，我们可以认为Zn-N是离子键。图中可以看到，DORI并不能判断定域性，无法展现孤对电子也无法区分共价键和离子键，但能够得到与ELF互补的信息。

总结

本文通过QE和VESTA，以GaAs和NiO为例讨论了ELF判断共价作用的方法，又通过H2O2@PMCS为例，对ELF和DORI各自的局限性和互补性进行了讨论。我们建议可以通过ELF对共价作用或非共价作用进行初步判断，若为非共价作用再使用DORI进一步对离子键、氢键、vdW作用等进行进一步的区分。希望本文能够为大家的科研工作助力添彩，同时也欢迎大家共同参与讨论互相进步。

注释

[1] 这个版本是由我们自己修改增加了DORI的函数以及修正了一些bug，数值上的实现可以参考本文附录，修复bug并增加DORI的代码我们已经通过gitlab提交给QE官方，后续版本可能会加入。如果要使用DORI，我们暂时将plot_num指定为123。
[2] 按照杨伟涛教授的文章[10]所述，密度小于0.005的区域对应于弱作用，密度介于0.005~0.05的区域则为强非共价作用。

参考文献

[1] 黄昆. 固体物理学[M]. 北京: 高等教育出版社, 1988.
[2] https://www.reed.edu/chemistry/ROCO/Resonance/delocalized.html.
[3] Silvi B, Savin A. Classification of chemical bonds based on topological analysis of electron localization functions. Nature. 1994, 371, 683–686. DOI: 10.1038/371683a0.
[4] De Silva P, Corminboeuf C. Simultaneous visualization of covalent and noncovalent interactions using regions of density overlap. J. Chem. Theory Comput. 2014, 10, 3745–3756. DOI: 10.1021/ct500490b.
[5] 固体和表面的非共价相互作用分析
[6] https://github.com/QEF/q-e/releases.
[7] http://www.jp-minerals.org/vesta/en/download.html.
[8] https://www.materialsproject.org.
[9] Xu B, Wang H, Wang W, Gao L, Li S, Pan X, Wang H, Yang H, Meng X, Wu Q, Zheng L. A Single-atom nanozyme for wound disinfection applications. Angew. Chem. 2019, 131, 4965–4970. DOI: 10.1002/anie.201813994.
[10] Johnson E R, Keinan S, Mori-Sánchez P, et al. Revealing noncovalent interactions[J]. J. Am. Chem. Soc, 2010, 132(18): 6498-6506.
DOI: 10.1021/ja100936w
[11] http://sobereva.com/367
[12] http://sobereva.com/63

附录1 基于FFT算法的数值微分方法

假定f(r)是一个光滑的三维周期函数，存在Fourier变换对（忽略常数）

G是倒格矢[1]， ，将上述关系代入就能够将微分形式转换为一般的代数形式

如果用F和F^(-1)分别表示Fourier变换和Fourier逆变换，那么上面的公式则可以表达为一个简单的形式

这意味着，f(r)的梯度计算可以依次对不同分量先进行Fourier变换，然后与该分量对应的倒格矢乘积运算后进行Fourier逆变换即可。实际编写代码时Fourier正逆变换可以直接使用如FFTW或FFTPACK等数学库，上述算法在QE中对应于QE目录下\Modules\gradutils.f90中的fft_gradient_r2r子程序。▽▽f(r)也可以按上述过程进行类似推导，此处不再赘述，其对应的是同一个源码文件中的fft_laplacian子程序。上述算法不仅能用于QE，对于其他程序也是可行的，包括cube或CHGCAR格式等，例如SIESTA程序的Util\Grid文件夹下有个cdf_laplacian.F90，这段代码也采用了相同的方法。

附录2 DORI的数值实现

按照定义，

我们先对分号的分子进行展开，这个只需要利用链式求导法则即可

其中▽ρ(r)为电子密度在r处的梯度，▽▽ρ(r)则是电子密度在r处的Hessian矩阵，可以使用附录1的方法进行计算。接着将分子分母同时乘以ρ(r)^6，就能得到以下公式：

实际计算时，我们在分母的|▽ρ(r)|^2加1e-5，即在不影响定性结果的基础上，来防止在体系边缘上因分母更快趋于0而使得DORI为1。

卡开发发

DORI的这个处理办法好，我之前没有展开直接做噪声还是比较大的。

sobereva

DORI这个函数其实很不好。最近我提出了IRI，形式比DORI简单得多得多，实现超级容易（只要已经实现了RDG，只需要改一行代码就得到IRI），而且图像效果好得多，使得DORI已经没有使用价值了。见
Tian Lu,* Qinxue Chen, Interaction Region Indicator (IRI): A Very Simple Real Space Function Clearly Revealing Both Chemical Bonds and Weak Interactions
预印版：ChemRxiv (2021) DOI: 10.26434/chemrxiv.13591142
（文章已经投出去了）

对比可见DORI效果很差，噪音非常大，不同的键之间都没有干净区分，两个片段间的范德华作用区域形状很丑。而IRI非常优雅地展现出所有相互作用区域，非共价作用区域保持了原有RDG的特征，化学键区域被平滑的等值面漂亮地展现了出来。



顺带一提，目前最新的Multiwfn已可以对CP2K产生的.molden文件计算ELF、LOL、IRI、DORI、RDG等各种依赖于波函数的函数，支持非正交盒子

例如氧化石墨烯+水的IRI（近期我会写个帖子介绍计算细节）


ELF





后记：IRI方法的论文已正式发表于Chemistry-Methods期刊，并且专门写了个博文做了详细介绍，见：
使用IRI方法图形化考察化学体系中的化学键和弱相互作用
http://sobereva.com/598（http://bbs.keinsci.com/thread-23457-1-1.html）

欧拉角

sobereva 发表于 2021-2-5 17:53
DORI这个函数其实很不好。最近我提出了IRI，形式比DORI简单得多得多，实现超级容易（只要已经实现了RDG，只 ...

欢迎卢老师一块参与讨论。IRI可能确实会是个好方法，如果我们测试下来合适也会考虑加入程序。本文的目的除了是提供这些分析方法，其实也是对数值实现方面简要做一些讨论，如您所述其实您的方法只需要对代码做很简单的修改就能够实现。其实无论是DORI、RDG，还是您提到的IRI都可能各自有其局限性，比如孤对电子方面无法体现等，最终还是要结合其他方法综合讨论，这不是很严重的问题。

sobereva

欧拉角 发表于 2021-2-5 18:19
欢迎卢老师一块参与讨论。IRI可能确实会是个好方法，如果我们测试下来合适也会考虑加入程序。本文的目的 ...

IRI的目的在于展现相互作用区域，孤对电子是ELF、LOL等描述电子定域性的函数的适合范畴

欧拉角

sobereva 发表于 2021-2-5 18:21
IRI的目的在于展现相互作用区域，孤对电子是ELF、LOL等描述电子定域性的函数的适合范畴

RDG和DORI或者IRI只要使用其中某种配合ELF或LOL其实就足以展现体系内相互作用情况了，如果还是不充分其实还可以补充别的分析方法，所以我认为也并不是非要用某一种方式去分析才行的。

卡开发发

sobereva 发表于 2021-2-5 17:53
DORI这个函数其实很不好。最近我提出了IRI，形式比DORI简单得多得多，实现超级容易（只要已经实现了RDG，只 ...

话说回来，cp2k功能是挺多，貌似不是很友好，稳定性也不是很好，我记得我以前看到过（http://bbs.keinsci.com/thread-4391-1-1.html），迄今也不知道这些bug到底解决了没，相比之下QE毕竟测试的数据还比较多（https://molmod.ugent.be/deltacodesdft）。另外对于一些小的晶胞或者一些金属乃至于金属氧化物表面得考虑多k点的问题，局限性会比较大。所以有个QE的能用的方法我认为还是很有价值的。

风飞

老师，你好，请问怎样修改pp. x文件呢，或者怎样获得修改后的pp. x文件呢？

欧拉角

风飞 发表于 2021-2-6 11:42
老师，你好，请问怎样修改pp. x文件呢，或者怎样获得修改后的pp. x文件呢？

我们已经在gitlab上向QE提交了，现在也是邮件沟通中，我们后续还可能会加入一些其他的功能。QE的官方也应该也会花时间对代码和注释进行调整，您要是不是太着急可以等QE下个版本，我们还不确定是否能够更新到官方版本中。如果您需要的是ELF，原始的ELF未考虑自旋极化形式，虽然可能定性不会有太大差异。如果您需要的是DORI来展现非共价作用，如同我在前面提到的，RDG也是可以替代的。在这里我主要还是希望提供一种解决问题的思路，不过有问题我欢迎大家一块讨论。

gog

sobereva 发表于 2021-2-5 17:53
DORI这个函数其实很不好。最近我提出了IRI，形式比DORI简单得多得多，实现超级容易（只要已经实现了RDG，只 ...

不同软件添加不同方法的代码实现的。比较的标准不一致。
从效果看，cp2k IRI更剩一筹。

卡开发发

gog 发表于 2021-8-4 17:04
不同软件添加不同方法的代码实现的。比较的标准不一致。
从效果看，cp2k IRI更剩一筹。

IRI的量纲太魔幻了，我想不出其物理意义，量纲这件事不说清楚的话，导致的结果就是在做单位换算的时候变得异常麻烦。

sobereva

卡开发发 发表于 2021-8-4 17:50
IRI的量纲太魔幻了，我想不出其物理意义，量纲这件事不说清楚的话，导致的结果就是在做单位换算的时候变 ...

密度和密度梯度都用a.u.代入，IRI也写a.u.即可
这个函数的初衷只是为了恰当的图形化展现，物理意义不是重点

卡开发发

sobereva 发表于 2021-8-6 08:01
密度和密度梯度都用a.u.代入，IRI也写a.u.即可
这个函数的初衷只是为了恰当的图形化展现，物理意义不是 ...

物理意义其次，量纲我认为是不应忽视的问题。比如对长度为Angstrom的单位如何转换呢？在原子单位下很多记作a.u.既可以是长度又可以是能量或其他，若不同程序使用的长度单位，显示后将出现不统一的情况（注：并不是所有图形显示程序都采用a.u.单位）。无量纲的RDG和DORI因为是无量纲的，在不同长度单位下不受影响，而Laplacian应当理解为能量密度相同的量纲，长度单位转换过程应该与密度处理相同才对；但是如果简单把Laplacian看成是密度的导数，则量纲看上去是长度^(-5)，如果这样转换再用某些以Angstrom为单位的显示程序显示可能就会得到怪异的结果。