请选择 进入手机版 | 继续访问电脑版
第9届北京科音分子动力学与GROMACS培训班将于4月17~20日于北京举办,请点击此链接查看培训详情,欢迎参加和相互转告!

计算化学公社

 找回密码
 现在注册!
查看: 1897|回复: 24

[Quantum ESPRESSO] 固体和表面的非共价相互作用分析

[复制链接]

21

帖子

0

威望

204

eV
积分
225

Level 3 能力者

发表于 2021-1-22 13:30:27 | 显示全部楼层 |阅读模式
管理员Sobereva 2021-Feb-27注:本论坛现在还有一篇介绍Multiwfn+CP2K联用便利地实现NCI分析的教程,见《使用Multiwfn结合CP2K通过NCI和IGM方法图形化考察固体和表面的弱相互作用》(http://bbs.keinsci.com/thread-21742-1-1.html


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

引言

约化密度梯度(reduced density gradient, RDG)是密度泛函理论中用来描述电子非均匀性的无量纲参量[注1],其表达式为
公式01.jpg
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. H2O2@PMCS来源于文献[6]
计算参数选择的依据:从文献[7]测试结果来看,vdW-DF3-opt2在S22×5和S66×8的表现都还不错,原文中所使用的是GBRV超软赝势[8],该赝势本身也有很好的可靠性和可移植性[9,10]。建议密度截断(ecutrho)选取充分以保证计算结果可靠以及避免RDG计算过程产生噪声。

计算参数

交换相关泛函vdW-DF3-opt2
赝势
GBRV
动能截断40 Ry
密度截断480 Ry
k网格MoS2    9×9×3
Ni(OH)2(Ni U=6.2 eV)    9×9×3
H2O2@PMCS    2×2×1
展宽 Gaussian    0.02 Ry


计算过程和数据处理

A1 使用PW模块
进行结构优化,并进行自洽计算产生波函数以及电子密度
  1. mpirun -np QE安装目录/bin/pw.x < MoS2.pw.in > MoS2.pw.out &
复制代码
A2 通过PW/tools/pwo2xsf将out文件转换为xsf格式[注3],如
  1. QE安装目录/PW/tools/pwo2xsf -oc MoS2.pw.out > MoS2.xsf
复制代码
A3 使用PP模块进行后处理,通过plot_num=19产生RDG格点,以及 output_format=5来产生xsf文件
  1. mpirun -np QE安装目录/bin/pp.x < MoS2.pp.in > MoS2.pp.out &
复制代码
A4 类似过程A3,使用PP模块通过plot_num=20产生sign(λ2)ρ的xsf文件,但建议提前将文件夹拷贝一份或者通过fileout更改xsf输出格点文件的名称避免上一步的结果被覆盖。
B1 打开VESTA程序,将A2中out文件产生的xsf文件拖入VESTA窗体
1.png
B2菜单栏中Edit->Edit Data->Volumetric Data,分别点击Isosurface和Surface coloring旁的Import分别导入RDG网格和sign(λ2)ρ的xsf文件。
2.png
B3 菜单栏中Objects->Properties->Isosurfaces,将Isosurfaces中的Isosurface level改为0.5,以及将Surface coloring中的Saturation level的Max改为0.02,Min改为-0.04,这样能够和文献[1]较好符合。
3.png
B4 调整其他显示部分,比如菜单中去除勾选Objects-> Volumetric Data->Show Sections,以及在Objects->Properties->Atoms对原子颜色进行调整,在Edit->Bonds对显示成键的判断进行调整,使得作图更加美观,但具体调整的细节这里不进行赘述。如果对操作过程有问题也欢迎大家讨论。

案例
4.png
MoS2是层状结构材料,MoS2层间出现了RDG等值面围成的区域,且层间S……S连线的区域呈现淡绿色(ρ<0.005),证明了层间存在相互作用,且通过vdW作用相结合并主要是以层间S……S的相互作用为主导;而层内的三个Mo原子之间显示为红色,表现出位阻效应。
5.png
Ni(OH)2与MoS2情况类似,但Ni(OH)2层间作用以O-H……O-H的氢键作用为主导;层内的Ni原子之间也呈现出位阻效应。在Ni-O键中形成的蓝色环形区域是因为Ni-O之间相互作用比较强,而pp.x程序在处理RDG的时候将ρ>0.05的网格直接设置为100[注4],导致共价键或离子键在图中不被显示。
6.png
H2O2@PMCS这一类结构常见于单原子催化过程中。文献[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近似就是将球对称的自由原子密度直接进行叠加,即
公式05.jpg
Promolecule近似下的密度通常可以作为自洽场的初始猜测,也能够用于大体系的定性计算。在QE中也能实现Promolecule近似的计算,这是因为QE中UPF格式的赝势文件中已经存有自由原子的电子密度。在自洽场计算过程中,我们只进行初始猜测而不做自洽过程,因此将electron_maxstep设置为0,且startingwfc设置为’atomic’,就可以得到Promolecule近似下的密度。
7.png
由于Promolecule近似做NCI分析也能给出定性可靠的结果[1],因此在处理更大的体系时可以使用这种方法。这里仍然以H2O2@PMCS体系为例,展示Promolecule近似下NCI分析的结果,下图表明Promolecule近似(等值面取0.035)与自洽计算后的差距并不大。

散点图

由于非共价相互作用区域电子密度和RDG都比较低,且λ2的符号为负,通常在散点图中呈现为一个或多个尖峰(spike)。如H2O2@PMCS的散点图中,左下角横坐标位于-0.04~0.02的spike就对应于H2O2中OH与N的氢键。
8.png
而promolecule近似下的散点图与上图的结果较为接近
9.png
那么如何绘制这样的散点图呢?我在这里提供了一个自己写的脚本nci_scatter.py nci_scatter.py (1.61 KB, 下载次数: 10)

评分

参与人数 16eV +71 收起 理由
anrushan + 2 谢谢分享
yjmaxpayne + 5 好物!
ChemG + 5 谢谢分享
biogon + 5
Henryugg + 4 好物!
ghifi37 + 3 精品内容
alonewolfyang + 5 好物!
冰释之川 + 5 精品内容
ggdh + 5 精品内容
朙天儿 + 5 好物!
ABetaCarw + 5 精品内容
newple + 5 好物!
asdf + 5 精品内容
978142355 + 5 666
ene + 5 好!
卡开发发 + 2 хорошо!

查看全部评分

学无止境。

2265

帖子

3

威望

8626

eV
积分
10951

Level 6 (一方通行)

Ab Initio Amateur

发表于 2021-1-22 13:47:09 | 显示全部楼层
这活整的有意思,也算是造福大家了,毕竟之前对固体方面讨论RDG的帖子不多见。说不上感谢,有些问题我也是请教@978142355。

实际上那天讨论完我发现VMD也是可以直接绘制RDG的等值面填色图,VMD自身就支持xsf文件的读取,另外需要指定一个pbc box -color black也能画出来差不多的图,不过实际上看上去效果差距不是很大。

另外,QE的sign(lambda2)rho和ELF问题我已经和QE官方做了反映,看看下个版本也许可能会改进。

评分

参与人数 2eV +7 收起 理由
Penson + 2 赞!
ene + 5 好!

查看全部评分

满招损,谦受益。热衷于理论和方法研究水平不高但欢迎讨论。

791

帖子

0

威望

2145

eV
积分
2936

Level 5 (御坂)

发表于 2021-1-22 14:10:22 来自手机 | 显示全部楼层
支持VASP吗?毕竟VASP用户也是个大族群

21

帖子

0

威望

204

eV
积分
225

Level 3 能力者

 楼主| 发表于 2021-1-22 14:16:38 | 显示全部楼层
granvia 发表于 2021-1-22 14:10
支持VASP吗?毕竟VASP用户也是个大族群

这个工作是用QE的pp.x实现的。VASP要实现需要自己编写程序或者修改代码,虽然能做但是在这边公开分享不是太合适。
学无止境。

3万

帖子

25

威望

3万

eV
积分
65560

管理员

公社社长+计算化学玩家

发表于 2021-1-22 17:41:37 | 显示全部楼层
Multiwfn可以载入CP2K的molden文件,回头我把代码改改支持考虑PBC,使得Multiwfn可以对CP2K的波函数做RDG分析而不需要考虑buffer区域

顺带一提,对于任意第一性原理程序优化出来的周期性结构,Multiwfn可以基于结构文件(xyz、pdb、mol2...)直接做基于promolecular近似的RDG分析

2021-Feb-27补充:相应帖子已经写好,见
使用Multiwfn结合CP2K通过NCI和IGM方法图形化考察固体和表面的弱相互作用
http://sobereva.com/588http://bbs.keinsci.com/thread-21742-1-1.html

评分

参与人数 1eV +5 收起 理由
nianbin + 5 谢谢

查看全部评分

北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395(已满),2号:466017436(已满)。3号:764390338(可加),合计8000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

701

帖子

23

威望

3226

eV
积分
4387

Level 6 (一方通行)

发表于 2021-1-22 18:40:41 | 显示全部楼层
sobereva 发表于 2021-1-22 17:41
Multiwfn可以载入CP2K的molden文件,回头我把代码改改支持考虑PBC,使得Multiwfn可以对CP2K的波函数做RDG分 ...

如果某个作用跨胞了,那Multiwfn默认是不支持这种作用的RDG分析吧,
这种情况下应该建立超胞然后用Multiwfn?

2265

帖子

3

威望

8626

eV
积分
10951

Level 6 (一方通行)

Ab Initio Amateur

发表于 2021-1-22 18:40:53 | 显示全部楼层
sobereva 发表于 2021-1-22 17:41
Multiwfn可以载入CP2K的molden文件,回头我把代码改改支持考虑PBC,使得Multiwfn可以对CP2K的波函数做RDG分 ...

非正交的格子可能还是这样的方式方便一些

评分

参与人数 1eV +2 收起 理由
ene + 2

查看全部评分

满招损,谦受益。热衷于理论和方法研究水平不高但欢迎讨论。

791

帖子

0

威望

2145

eV
积分
2936

Level 5 (御坂)

发表于 2021-1-22 20:55:47 来自手机 | 显示全部楼层
对于DFT-D级别的计算,色散作用也属于非共价作用,这个NCI分析应该没法处理了吧?

21

帖子

0

威望

204

eV
积分
225

Level 3 能力者

 楼主| 发表于 2021-1-22 22:01:00 | 显示全部楼层
granvia 发表于 2021-1-22 20:55
对于DFT-D级别的计算,色散作用也属于非共价作用,这个NCI分析应该没法处理了吧?

我们测试过PBE-D3,对结果的影响并不是很显著。
学无止境。

3万

帖子

25

威望

3万

eV
积分
65560

管理员

公社社长+计算化学玩家

发表于 2021-1-23 06:41:45 | 显示全部楼层
granvia 发表于 2021-1-22 20:55
对于DFT-D级别的计算,色散作用也属于非共价作用,这个NCI分析应该没法处理了吧?

一般意义的NCI是基于电子密度做的,DFT-D不影响密度分布
至于诸如VV10那种可以在迭代过程中影响密度的色散校正,对密度影响微乎其微
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395(已满),2号:466017436(已满)。3号:764390338(可加),合计8000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

3万

帖子

25

威望

3万

eV
积分
65560

管理员

公社社长+计算化学玩家

发表于 2021-1-23 06:43:21 | 显示全部楼层
ggdh 发表于 2021-1-22 18:40
如果某个作用跨胞了,那Multiwfn默认是不支持这种作用的RDG分析吧,
这种情况下应该建立超胞然后用Multi ...

目前如此。我近期打算在计算密度时显式地考虑PBC,免得先扩胞了
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395(已满),2号:466017436(已满)。3号:764390338(可加),合计8000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

49

帖子

0

威望

271

eV
积分
320

Level 3 能力者

发表于 2021-1-23 08:17:24 | 显示全部楼层
太好了,非常详细的教程,欧拉老师能否提供一下 MoS2.pw.in  文件,好让我们这些想入门的人能够按照教程操练一下?

21

帖子

0

威望

204

eV
积分
225

Level 3 能力者

 楼主| 发表于 2021-1-23 09:17:13 | 显示全部楼层
renzhogn424 发表于 2021-1-23 08:17
太好了,非常详细的教程,欧拉老师能否提供一下 MoS2.pw.in  文件,好让我们这些想入门的人能够按照教程操 ...


需要注意的是,vdW-DF3-opt2应该只有QE新版本能使用,如上面所述实际用PBE-D3BJ差异也不太大。 MoS2.pw.in (1.42 KB, 下载次数: 10)
学无止境。

791

帖子

0

威望

2145

eV
积分
2936

Level 5 (御坂)

发表于 2021-1-23 11:50:00 来自手机 | 显示全部楼层
sobereva 发表于 2021-1-23 06:41
一般意义的NCI是基于电子密度做的,DFT-D不影响密度分布
至于诸如VV10那种可以在迭代过程中影响密度的色 ...

那是不是说基于电子密度的分析并不能用于解释色散作用?

3万

帖子

25

威望

3万

eV
积分
65560

管理员

公社社长+计算化学玩家

发表于 2021-1-23 15:18:07 | 显示全部楼层
granvia 发表于 2021-1-23 11:50
那是不是说基于电子密度的分析并不能用于解释色散作用?

基于电子密度做RDG分析,是可以图形化展现弱相互作用区域和强度的(只是“图形化展现”,不是“解释”本质)

从DFT原理上来讲,通过基态电子密度也是可以严格计算色散作用的,不过这和RDG不是一个范畴的事
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395(已满),2号:466017436(已满)。3号:764390338(可加),合计8000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )

GMT+8, 2021-4-15 06:51 , Processed in 0.279566 second(s), 32 queries .

快速回复 返回顶部 返回列表