请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 720|回复: 0

[Multiwfn资源与经验] 使用Multiwfn超级方便地计算出概念密度泛函理论中定义的各种量

[复制链接]

1万

帖子

25

威望

2万

eV
积分
42584

管理员

公社社长

发表于 2019-5-21 12:15:20 | 显示全部楼层 |阅读模式
使用Multiwfn超级方便地计算出概念密度泛函理论中定义的各种量

文/Sobereva@北京科音  2019-May-21

1 前言

由Parr最初开始发展的概念密度泛函理论(Conceptual density functional theory,CDFT)也称密度泛函反应性理论(DFRT),是量子化学和波函数分析领域中的一个重要组成,专门用来预测和解释化学物质的反应活性、反应位点等问题,对于研究化学反应有重要意义。相关介绍以及笔者写过的相关博文见《概念密度泛函综述和重要文献合集》(http://bbs.keinsci.com/thread-384-1-1.html)。

CDFT框架里有大量概念,包含很多实空间函数、局部指数和全局指数,它们经常出现在各种文献里用来讨论反应问题。只要有N、N+1、N-1电子态的波函数和能量,就可以计算它们中的绝大部分。N指的是体系的电子数,比如某个分子最稳定状态是中性的,则N就是指中性时的电子数,N+1就是指带一个负电荷的阴离子状态的电子数,N-1就是指带一个正电荷的阳离子态的电子数。其实计算CDFT里的那些量并不复杂,手动算起来没什么难度,然而如果把所有量全都一一算出来,还是比较费事,而且容易弄错数据。另外,网上也经常有人问我怎么去算那些量,由于初学者的量子化学基础知识往往很薄弱,经常解释起来很费劲,甚至怎么说也说不明白。于是,笔者决定在Multiwfn里直接加入一个功能,通过极简单的步骤,就可以一次性把所有CDFT框架内常见的量全都算出来,过程完全傻瓜化。笔者相信这个功能有极高的实用性,经常要计算CDFT里涉及的量的读者务必要重视此功能。

这个功能从2019-May-21发布的Multiwfn 3.6正式版开始加入,是主功能100的子功能16。Multiwfn可以在官网http://sobereva.com/multiwfn免费下载。相关入门知识看《Multiwfn入门tips》(http://sobereva.com/167)、《Multiwfn FAQ》(http://sobereva.com/452)。


2 能计算的量

这个功能能直接计算的量主要包含下面这些,分为三大类,具体定义在这里就不写了,因为在Multiwfn手册3.100.16节都已经明确给出了。
• 全局指数(即对整个体系而言的数值)
垂直电离能(Vertical ionization potential,VIP)
垂直电子亲和能(Vertical electron affinity,VEA)
Mulliken电负性(Mulliken electronegativity)
化学势(Chemical potential)
电子的硬度(Hardness)
电子的软度(Softness)
亲电指数(Electrophilicity index)
亲核指数(Nucleophilicity index)
• 实空间函数(即变量是三维空间坐标的函数)
福井函数(Fukui function)
双描述符(Dual descriptor)
局部软度(Local softness)
局部亲电指数(Local electrophilicity index)
局部亲核指数(Local nucleophilicity index)
• 原子指数(即每个原子各有一个数值)
简缩福井函数(Condensed Fukui function)
简缩双描述符(Condensed dual descriptor)
简缩局部软度(Condensed local softness)
相对亲电指数(Relative electrophilicity index)
相对亲核指数(Relative nucleophilicity index)
简缩局部亲电指数(Condensed local electrophilicity index)
简缩局部亲核指数(Condensed local nucleophilicity index)

上述原子指数都可以写为基于原子电荷计算的形式。根据一些文章的测试,Hirshfeld电荷非常适合这种目的(不建议用NPA等其它电荷),所以Multiwfn自动计算它们都是用的Hirshfeld电荷。


3 实例:对苯酚计算CDFT定义的各种量

这个而功能的用法和细节在手册3.100.16节介绍了,本文就不再累述,下面直接看一个具体使用例子,对苯酚计算各种量。看完例子后建议再看一下手册这一节。

3.1 准备.wfn文件

用本文介绍的功能计算各种量之前需要先在当前目录下给出N.wfn、N+1.wfn和N-1.wfn,分别记录当前体系三个电子态的波函数信息。可以自行提供它们,但如果你是Gaussian用户,按照下面的流程操作省事得多。

计算前述所有量要用的几何结构都是相同的,都是N电子态的极小点结构。对于苯酚这个例子而言,我们需要对苯酚中性状态进行几何优化,这一步是用户自己手动做的,用什么程序、什么级别无所谓,靠谱就行。普通有机体系用常用的B3LYP-D3(BJ)/6-31G*级别一般就够了,如果想更便宜一些,用半经验、GFN-xTB之类也可以。在B3LYP/6-31G*下优化中性苯酚得到的结构已经直接给了,是Multiwfn目录下的examples目录中的phenol.xyz。

启动Multiwfn,依次输入
examples\phenol.xyz  //一开始载入的文件只要含有中性苯酚优化后的结构即可,其它格式如pdb/fch/gjf/mol/mol2等也都可以
100  //其它功能(Part 1)
16  //计算CDFT里定义的各种量的模块
[直接按回车]  //这一步让你输入gjf文件里对应的计算单点任务的关键词。按回车代表用默认的B3LYP/6-31G*,此级别虽然很便宜,但对于考察CDFT定义的量来说一般基本已可以满足需要
[直接按回车]  //这一步让你输入计算N、N+1和N-1态用的电荷和自旋多重度。按回车代表用对三个态分别用默认的0 1、-1 2和1 2

现在当前目录下已经有了N.gjf、N+1.gjf和N-1.gjf,是用于产生.wfn文件的单点任务的Gaussian输入文件。如果你不自己改gjf里的.wfn文件路径的话,把这些gjf都用Linux下的Gaussian执行后,N.wfn、N+1.wfn和N-1.wfn会产生在当前目录下,而如果用Windows下的Gaussian执行则会产生在临时文件目录下(比如D:\study\G16W\Scratch下)。如果你当前机子里已经装了Gaussian,而且之前Multiwfn的settings.ini里的gaupath已经被你设为了实际的Gaussian可执行文件路径,那么建议此时直接在Multiwfn窗口里输入y,这样Multiwfn就会自动调用Gaussian计算这个三个gjf文件,并在当前目录下产生相应的三个wfn文件,然后gjf和out文件会自动被删掉。
注:让Multiwfn顺利调用Windows版Gaussian需要做额外的设置,见Multiwfn手册Appendix 1。


3.2 计算各种指数

现在当前目录下已经有了N.wfn、N+1.wfn和N-1.wfn,可以开始算了。在当前模块里选择2 Calculate various quantitative indices,然后程序就会依次载入wfn文件,读取其中的能量信息和波函数信息,自动计算Hirshfeld电荷,最后给出各种指数。对于6-31G*基组下的苯酚,在Intel 4核机子下不出10秒钟就都算完了。算完的数据被输出到了当前目录下的CDFT.txt中,此文件内容如下:
Note: the E(HOMO) of TCE used for evaluating nucleophilicity index is the value evaluated at B3LYP/6-31G* level

Hirshfeld charges, condensed Fukui functions and condensed dual descriptors
Units used below are "e" (elementary charge)
     Atom     q(N)    q(N+1)   q(N-1)     f-       f+       f0      CDD
     1(C )  -0.0587  -0.1185   0.0852   0.1439   0.0598   0.1018  -0.0841
     2(C )  -0.0390  -0.1674   0.0268   0.0658   0.1284   0.0971   0.0626
     3(C )  -0.0597  -0.1873   0.0319   0.0916   0.1276   0.1096   0.0360
     4(C )   0.0737   0.0216   0.1739   0.1001   0.0522   0.0762  -0.0480
     5(C )  -0.0731  -0.1956   0.0090   0.0821   0.1225   0.1023   0.0404
     6(C )  -0.0415  -0.1730   0.0336   0.0751   0.1315   0.1033   0.0563
     7(H )   0.0417  -0.0048   0.0993   0.0576   0.0464   0.0520  -0.0112
     8(H )   0.0450  -0.0193   0.0904   0.0455   0.0642   0.0548   0.0187
     9(H )   0.0473  -0.0160   0.0953   0.0480   0.0633   0.0557   0.0153
    10(H )   0.0387  -0.0230   0.0854   0.0467   0.0617   0.0542   0.0150
    11(H )   0.0444  -0.0207   0.0914   0.0470   0.0651   0.0561   0.0181
    12(O )  -0.1977  -0.2443  -0.0517   0.1460   0.0465   0.0963  -0.0995
    13(H )   0.1789   0.1482   0.2294   0.0505   0.0307   0.0406  -0.0197

Condensed local electrophilicity/nucleophilicity index (e*eV)
     Atom              Electrophilicity          Nucleophilicity
     1(C )                  0.02576                  0.45535
     2(C )                  0.05533                  0.20827
     3(C )                  0.05496                  0.28982
     4(C )                  0.02249                  0.31688
     5(C )                  0.05280                  0.25977
     6(C )                  0.05665                  0.23773
     7(H )                  0.02001                  0.18234
     8(H )                  0.02767                  0.14394
     9(H )                  0.02729                  0.15195
    10(H )                  0.02659                  0.14779
    11(H )                  0.02807                  0.14875
    12(O )                  0.02005                  0.46203
    13(H )                  0.01325                  0.15965

Condensed local softnesses (Hartree*e) and relative electrophilicity/nucleophilicity (dimensionless)
     Atom         s-          s+          s0        s+/s-       s-/s+
     1(C )      0.3761      0.1562      0.2661      0.4154      2.4076
     2(C )      0.1720      0.3356      0.2538      1.9510      0.5126
     3(C )      0.2393      0.3333      0.2863      1.3926      0.7181
     4(C )      0.2617      0.1364      0.1990      0.5211      1.9190
     5(C )      0.2145      0.3202      0.2674      1.4925      0.6700
     6(C )      0.1963      0.3436      0.2700      1.7500      0.5714
     7(H )      0.1506      0.1213      0.1360      0.8057      1.2412
     8(H )      0.1189      0.1678      0.1433      1.4116      0.7084
     9(H )      0.1255      0.1655      0.1455      1.3187      0.7583
    10(H )      0.1221      0.1612      0.1416      1.3211      0.7570
    11(H )      0.1228      0.1702      0.1465      1.3858      0.7216
    12(O )      0.3816      0.1216      0.2516      0.3186      3.1384
    13(H )      0.1319      0.0804      0.1061      0.6094      1.6409

E(N):     -307.464860 Hartree
E(N+1):   -307.383614 Hartree
E(N-1):   -307.163438 Hartree
E_HOMO(N):     -0.218913 Hartree,   -5.9569 eV
E_HOMO(N+1):    0.161297 Hartree,    4.3891 eV
E_HOMO(N-1):   -0.464864 Hartree,  -12.6496 eV
Vertical IP:    0.301421 Hartree,    8.2021 eV
Vertical EA:   -0.081246 Hartree,   -2.2108 eV
Mulliken electronegativity:     0.110088 Hartree,    2.9956 eV
Chemical potential:            -0.110088 Hartree,   -2.9956 eV
Hardness (=fundamental gap):    0.382667 Hartree,   10.4129 eV
Softness:    2.613235 Hartree^-1,    0.0960 eV^-1
Electrophilicity index:    0.015835 Hartree,    0.4309 eV
Nucleophilicity index:     0.116285 Hartree,    3.1643 eV

上面的数据基本没有什么需要特别解释的,信息都非常易懂,如果有符号看不懂的,去对照手册3.100.16节的式子。唯一值得特别一提的是亲核指数(以及局部亲核指数),这个量是依赖于TCE(四氰基乙烯)的HOMO能量定义的,TCE的HOMO能量按理说是应当使用和当前计算级别完全一样的级别来计算的,这样得到的全局和局部亲核指数才是严格的,但Multiwfn在给出的时候是自动用笔者事先算好的B3LYP/6-31G*下的TCE的HOMO值-0.335198 Hartree算的。如果你当前用的不是B3LYP/6-31G*级别且追求更严格的结果,应当自行用当前级别去算TCE的HOMO值然后手动根据定义得到全局和局部亲核指数。

在Multiwfn手册里也有手动计算简缩福井函数和简缩双描述符的例子,见4.7.3节,用的也是Hirshfeld电荷和B3LYP/6-31G*级别,因此结果和上面自动算出来的完全一样。相比之下,用本文介绍的功能实在方便太多了,所有重要的量一口气全都输出了!


3.3 计算福井函数和双描述符

接下来计算福井函数和双描述符。在当前模块里选择3 Calculate grid data of Fukui function and dual descriptor,然后选择一个合适的格点设定,对于当前这样小体系选择Medium quality grid就够了(大体系应当用High quality grid或其它选项,详见http://sobereva.com/452中的Q39的讨论),然后程序会依次载入当前目录下的N.wfn、N+1.wfn和N-1.wfn并计算电子密度格点数据,之后会看到一个菜单。通过选择相应选项,可以把各种类型的福井函数以及双描述符直接显示成等值面图。比如此例我们依次选择选项1、2、3、4来把f+、f-、f0福井函数和双描述符都依次绘制出来,下图是把图像手动合并到一起的图,等值面数值都是用的0.007:
1.png

通过主功能5手动计算福井函数和双描述符的方法在Multiwfn手册4.5.4节已经详细介绍了,得到的图像和上图完全一样。使用本文介绍的步骤,明显远比手动操作方便得多得多得多!

值得一提的是,如果你想要更好的显示质量,可以借助VMD。比如我们选择6 Export grid data of f- as f-.cub in current folder,然后把当前目录下产生的f-.cub按照《在VMD里将cube文件瞬间绘制成效果极佳的等值面图的方法》(http://sobereva.com/483)里介绍的方法载入VMD并通过Tachyon渲染器渲染,只需要很简单几步就可以得到下图,可见效果极好!
2.png

另外,当前功能也可以对局部软度、局部亲电/亲核指数观看等值面和导出cube文件,因为它们的定义都是对特定类型的福井函数乘上一个数值,而Multiwfn当前功能里可以设乘上的系数。比如局部软度s-定义为全局软度乘以福井函数f-,我们根据前面算出的数据已知苯酚的局部软度是2.613235 Hartree^-1,因此若我们想导出局部软度cube文件,就选-1 Set the value multiplied to all grid data,输入2.613235,之后选6 Export grid data of f- as f-.cub in current folder导出的cube文件就对应于s-的格点数据了,若选2 Visualize isosurface of f-则看到的是s-的等值面。


4 总结

本文介绍了Multiwfn的超级便利的一口气算出来概念密度泛函理论里涉及到的各种量的功能,这使得初学者也可以非常顺利、快速且不犯错误地计算这些量。相信此功能对于将概念密度泛函理论广泛地应用于实际化学问题的研究有积极的促进作用,也鼓励大家在日常研究中充分使用此功能,以从繁复的手动操作中解脱出来。

评分

参与人数 1eV +5 收起 理由
cokoy + 5 好物!

查看全部评分

北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

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

GMT+8, 2019-6-19 01:56 , Processed in 0.166920 second(s), 28 queries .

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