引言

约化密度梯度(reduced density gradient, RDG)是密度泛函理论中用来描述电子非均匀性的无量纲参量[注1],其表达式为

2010年,杨伟涛教授提出基于RDG的非共价相互作用(noncovalent interactions, NCI)[1]分析方法,可以用来归属原子间或分子间的相互作用。虽然电子定域化函数(electron localization function,ELF)和分子中的原子(atoms-in-molecules, AIM)理论等能够用来分析化学键,但是对氢键、空间位阻、π-π堆叠等非共价相互作用的讨论有一定的局限性,而NCI分析则可以作为揭示这些非共价相互作用强有力的工具。因为分子的电子密度呈现指数衰减,所以在远离分子处RDG呈现很大值;相反,在相互作用的区域RDG的值则很小,其中非共价相互作用区域的电子密度和RDG都比较低。

由于氢键和空间位阻区域的电子密度和RDG很接近而难以进行区分,所以需要再借助通过电子密度Hessian矩阵的第二大本征值的正负号(sign(λ2))来进一步判断[注2];电子密度大小则可以反映出相互作用的强度,一般认为vdW作用的区域ρ<0.005,因此可以通过对RDG等值面进行sign(λ2)ρ填色或者做RDG vs sign(λ2)ρ的散点图来讨论所研究体系的相互作用类型。

在Sobereva的博文中已经介绍了RDG相关的概念和使用Multiwfn对非周期边界系进行NCI分析的详细过程[2]。但固体和表面以及牵扯过渡金属的体系的NCI讨论较少。因此,本文将通过若干实例的计算和相应的分析,对这些体系使用Quantum ESPRESSO(QE)计算RDG的过程和NCI分析方法进行介绍。


准备工作

本文所采用的计算程序如下:

1. Quantum ESPRESSO 6.7 [3] 

实际上QE在5.1.0版本就已经支持NCI分析,所以即使正在使用旧版本的用户也不必担心无法完成本文的工作。

2. VESTA 3.1.8 [4] 

个别新版本有bug,所以本文采用旧版本。

测试体系为:

1. 来源于Materials Project [5]的结构

MoS2(mp-2815)和Ni(OH)2(mp-27912)

2. [email protected]来源于文献[6]


计算参数

交换相关泛函

vdW-DF3-opt2

赝势

超软

GBRV

动能截断

40 Ry

密度截断

480 Ry

k网格

MoS2

9×9×3

Ni(OH)2(Ni U=6.2 eV)

9×9×3

[email protected]

2×2×1

展宽

Gaussian

0.02 Ry

计算参数选择的依据:从文献[7]测试结果来看,vdW-DF3-opt2在S22×5和S66×8的表现都还不错,原文中所使用的是GBRV超软赝势[8],该赝势本身也有很好的可靠性和可移植性[9,10]。建议密度截断(ecutrho)选取充分以保证计算结果可靠以及避免RDG计算过程产生噪声。


计算过程和数据处理

A1 使用PW模块

进行结构优化,并进行自洽计算产生波函数以及电子密度

A2 通过PW/tools/pwo2xsf将out文件转换为xsf格式[注3],如

A3 使用PP模块进行后处理,通过plot_num=19产生RDG格点,以及 output_format=5来产生xsf文件

A4 类似过程A3,使用PP模块通过plot_num=20产生sign(λ2)ρ的xsf文件,但建议提前将文件夹拷贝一份或者通过fileout更改xsf输出格点文件的名称避免上一步的结果被覆盖。

B1 打开VESTA程序,将A2中out文件产生的xsf文件拖入VESTA窗体

B2菜单栏中Edit->Edit Data->Volumetric Data,分别点击Isosurface和Surface coloring旁的Import分别导入RDG网格和sign(λ2)ρ的xsf文件。

B3 菜单栏中Objects->Properties->Isosurfaces,将Isosurfaces中的Isosurface level改为0.5,以及将Surface coloring中的Saturation level的Max改为0.02,Min改为-0.04,这样能够和文献[1]较好符合。

B4 调整其他显示部分,比如菜单中去除勾选Objects-> Volumetric Data->Show Sections,以及在Objects->Properties->Atoms对原子颜色进行调整,在Edit->Bonds对显示成键的判断进行调整,使得作图更加美观,但具体调整的细节这里不进行赘述。如果对操作过程有问题也欢迎大家讨论。


案例

MoS2是层状结构材料,MoS2层间出现了RDG等值面围成的区域,且层间S……S连线的区域呈现淡绿色(ρ<0.005),证明了层间存在相互作用,且通过vdW作用相结合并主要是以层间S……S的相互作用为主导;而层内的三个Mo原子之间显示为红色,表现出位阻效应。

Ni(OH)2与MoS2情况类似,但Ni(OH)2层间作用以O-H……O-H的氢键作用为主导;层内的Ni原子之间也呈现出位阻效应。在Ni-O键中形成的蓝色环形区域是因为Ni-O之间相互作用比较强,而pp.x程序在处理RDG的时候将ρ>0.05的网格直接设置为100[注4],导致共价键或离子键在图中不被显示。

[email protected]这一类结构常见于单原子催化过程中,也是我们计算较多的一类结构。文献[6]给出的吸附能为-0.45 eV,仅从吸附能来看可能是物理吸附。尽管本文采用的方法与文献不同,但给出了-0.47 eV的结果与文献[6]十分接近。

为了进一步确定H2O2与基底相互作用的类型,我们进行了NCI分析,结果表明过氧化氢与Zn-N4的相互作用却是是非共价相互作用,其中H2O2中一个O-H与Zn-N4中的N原子形成了氢键;而另一个O-H的O与Zn原子产生了静电作用,这种静电作用从颜色上看要更强于氢键作用,因此可以认为H2O2在PMCS的吸附作用是由这种静电作用所主导的物理吸附。


Promolecule近似

Promolecule近似就是将球对称的自由原子密度直接进行叠加,即

Promolecule近似下的密度通常可以作为自洽场的初始猜测,也能够用于大体系的定性计算。在QE中也能实现Promolecule近似的计算,这是因为QE中UPF格式的赝势文件中已经存有自由原子的电子密度。在自洽场计算过程中,我们只进行初始猜测而不做自洽过程,因此将electron_maxstep设置为0,且startingwfc设置为’atomic’,就可以得到Promolecule近似下的密度。

由于Promolecule近似做NCI分析也能给出定性可靠的结果[1],因此在处理更大的体系时可以使用这种方法。这里仍然以[email protected]体系为例,展示Promolecule近似下NCI分析的结果,下图表明Promolecule近似(等值面取0.035)与自洽计算后的差距并不大。


散点图

由于非共价相互作用区域电子密度和RDG都比较低,且λ2的符号为负,通常在散点图中呈现为一个或多个尖峰(spike)。如[email protected]的散点图中,左下角横坐标位于-0.04~0.02的spike就对应于H2O2中OH与N的氢键。

而promolecule近似下的散点图与上图的结果较为接近

那么如何绘制这样的散点图呢?我在这里提供了一个自己写的脚本nci_scatter.py能够直接绘制散点图。需要预先安装python 3.x,调用到的python组件有numpy和matplotlib,只要通过pip install就能够进行安装。只要执行nci_scatter.py脚本按照每一步的提示用鼠标拖入对应的xsf文件就能够直接画出上面的散点图。

“nci_scatter.py”

链接:

https://pan.baidu.com/s/1rtUFsfvlFPtLiqROc2qNHw

提取码:csgo

请大家请自行下载~


总结

Quantum ESPRESSO是一款功能十分强大且开源免费的第一性原理计算程序,NCI分析也是一种强大的非共价相互作用分析方法。QE已经内置了RDG的计算,但却很少有人写QE做NCI分析的教程,本文以固体和表面中的一些常见体系为案例,对实现过程和分析进行了一些叙述讨论,更多的理论基础和分析细节还是建议大家多研读文献。希望本文能够为大家的科研工作助力添彩,同时也欢迎大家共同参与讨论互相进步。


注释

[1] RDG是无量纲参量,这是因为

[2] 这些概念来源于AIM理论,电子密度的Hessian矩阵是一个3×3的矩阵,对角化后得到三个本征值λ1>λ2>λ3。若λ2<0意味着λ3也小于0,对应于键临界点(bond critical point,BCP);若λ2>0意味着λ1也大于0,对应于环临界点(ring critical point, RCP)。

[3] pwo2xsf的作用是将pw计算得到的out文件转化为xsf文件,可选选项有

-ic 初始结构

-lc 最终结构

-oc 最优结构

-a 将轨迹转化为AXSF

[4] 严谨起见,我们检查了QE的pp.x程序的源代码,发现在QE 6.x版本分析磁性体系的sign(λ2)ρ和ELF时,电子密度存在bug(这是因为自旋极化的计算中,QE 6.x版本的电子密度不是按照spin up和spin down进行存储的,而是按照total和spin进行存储的,但sign(λ2)ρ和ELF这两段代码忘记进行修改)。这点我们已经向QE官方提交了这个bug并得到了确认,所以先建议大家不要分析磁性体系,后续版本会进行改进。本文中的案例Ni(OH)2的结果展示是我们使用修改后的程序,因此结果可靠。


参考资料

[1] Johnson E R, Keinan S, Mori-Sánchez P, et al. Revealing noncovalent interactions[J]. Journal of the American Chemical Society, 2010, 132(18): 6498-6506.

DOI: 10.1021/ja100936w

[2] http://sobereva.com/68

[3] https://github.com/QEF/q-e/releases

[4] http://www.jp-minerals.org/vesta/en/download.html

[5] https://www.materialsproject.org

[6] Xu B, Wang H, Wang W, et al. A Single‐Atom Nanozyme for Wound Disinfection Applications[J]. Angewandte Chemie, 2019, 131(15): 4965-4970.

DOI: 10.1002/anie.201813994

[7] Chakraborty D, Berland K, Thonhauser T. Next-Generation Nonlocal van der Waals Density Functional[J]. Journal of Chemical Theory and Computation, 2020, 16(9): 5893-5911.

DOI: 10.1021/acs.jctc.0c00471

[8] Garrity K F, Bennett J W, Rabe K M, et al. Pseudopotentials for high-throughput DFT calculations[J]. Computational Materials Science, 2014, 81: 446-452.

DOI: 10.1016/j.commatsci.2013.08.053

[9] Marsman M, Marzari N, Nitzsche U, et al. Reproducibility in density functional theory calculations of solids[J]. Science, 2016, 351: 6280.

DOI: 10.1126/science.aad3000

[10] https://molmod.ugent.be/deltacodesdft


学术交流微信群


测试狗同步辐射咨询预约群 / 同步辐射GIWAXS/GISAXS交流群 / OPV钙钛矿制备与表征交流群 / 测试GO-XRD精修交流群 / 催化理论与实验交流群 / 分子模拟应用交流群 / 电催化交流群 / 测试狗推文创作群/CSG学术干货分享群 / CSG能源材料交流
加微信群方式:添加GO小妹微信15680750806,
请备注“昵称—学校/机构”


学术交流QQ群


测试狗EBSD学习交流QQ群:597369268
机器学习+催化计算QQ群:966186610
Lammps分子模拟实践交流群:1137524036
材料测试&模拟计算QQ群:776761012
科研字幕组(科学10分钟)QQ群:703443384
测试GO软件资源共享QQ群:567425397