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

计算化学公社

 找回密码
 现在注册!
查看: 1756|回复: 10

[Multiwfn资源与经验] Multiwfn的轨道定域化功能的使用以及与NBO、AdNDP分析的对比

[复制链接]

1万

帖子

25

威望

1万

eV
积分
31051

管理员

公社社长

发表于 2017-6-5 15:10:18 | 显示全部楼层 |阅读模式
Multiwfn的轨道定域化功能的使用以及与NBO、AdNDP分析的对比

文/Sobereva @北京科音
First release: 2017-Jun-5    Last update: 2018-Feb-19



摘要:本文首先介绍分子轨道定域化的含义,然后介绍在Multiwfn中做轨道定域化的基本步骤。然后通过5个例子,讲解如何基于定域化分子轨道讨论、分析体系的电子结构,并且与流行的AdNDP方法和NBO方法进行对比。本文的实例体现出轨道定域化的重要应用价值,也鲜明反映出NBO分析对于电子结构复杂、牵扯多中心离域的体系结果经常不靠谱,一定要慎用,最好不用。

1 谈谈分子轨道的定域化

我们做一般量子化学计算产生的分子轨道叫正则分子轨道(Canonical MO),它们是Fock或者KS算符的本征函数,本征值就是轨道的能量。下面说的MO都是指正则分子轨道。MO往往有很强的离域性,也就是说往往MO涉及很多原子,因此与化学键没法一一对应。所以,直接通过分子轨道讨论成键问题,多数情况毫无用处。有不少文章、毕业论文摆出来一堆MO图,但要么没怎么讨论(纯粹是凑数据),要么就是瞎讨论。例如,下面的图是丁烷的两个MO图,完全离域在整个分子。显然从这种MO图上根本没法直接对应众所周知的此体系存在C-C、C-H sigma键的事实。
1.png


要想将轨道与化学键直接联系起来,需要做轨道定域化,将MO变换成定域化分子轨道(localized MO, LMO)。轨道定域化相当于对占据的MO做酉变换,从而产生出相同数目的离域程度尽可能低的占据轨道。描述体系波函数的Slater行列式无论是由原先的MO所构成,还是由变换后的LMO所构成,所对应的体系各种可观测量,如体系能量、电子密度、静电势等等,都是不变的。还有很多量虽然是不可观测的,但也不受这种变换所影响,比如电子定域化函数(ELF)。

通常做轨道定域化只需要把占据轨道定域化,因为可观测性质、成键问题等等都是只与占据轨道相联系。如果也想把空轨道做定域化,那么也可以再把空轨道以相同方式做酉变换来得到离域程度尽量低的空轨道。将空轨道做定域化主要是用于一些基于定域轨道的电子相关方法,比如知名的但如今用处不大的Local MP2 (LMP2)。由于空轨道比占据轨道数目多得多,将它们也定域化耗时会增加甚巨,而且空轨道对我们分析体系的性质没影响,所以通常我们只把占据轨道定域化就够了。

轨道定域化有许多做法。比较知名的是
(1)Edmiston-Ruedenberg定域化:1963年提出。通过最大化轨道的自互斥能(亦即最小化轨道间的互斥能)来实现定域化。缺点是耗时很高,而且得算复杂的双电子积分,而且又没额外好处,故很不推荐。
(2)Boys定域化:1960年提出。应用比较广泛,通过最小化轨道涵盖的空间范围来定域化。需要利用偶极矩积分。
(3)Pipek-Mezey定域化:1989年提出,它通过通过最大化Mulliken原子布居数的平方和达到定域化目的,编程实现颇为简单,只需要重叠积分,而且耗时很低,被广泛使用。后来还有人提出基于其它布居分析方法来做Pipek-Mezey定域化。
(4)NLMO定域化:NLMO的N是Natural的意思。这是NBO开发者提出的定域化方法,它将高占据(Lewis型)的NBO和低占据(non-Lewis型)的NBO以一定方式混合,从而得到整数占据的定域化轨道。由于得到NLMO轨道得先产生NBO,故此方法只有NBO程序,以及也同样能做NBO分析的Molpro才支持。
以上方法计算过程需要迭代,迭代过程中不同轨道间不断混合,最后使得目标函数被最大化或最小化。

不同的轨道定域化方法得到的LMO图形有一定差异。对于双键的描述,不同方法差异很大,比如下图是乙烯的C=C键,Pipek-Mezey产生的LMO保留了sigma和pi分离的特征,然而Boys定域化则是以两个香蕉键来描述这个双键。


2.png

这两种双键的描述形式没法说谁对谁错,从物理角度上是等价的描述,但是显然Pipek-Mezey方法更好一些,结果更符合化学观念,更便于讨论。

Pipek-Mezey方法的一个缺点是不适合有弥散函数的情况,因为Mulliken布居在有弥散函数的时候结果缺乏意义。这点在《原子电荷计算方法的对比》(http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml)、《分子轨道成分的计算》(http://sioc-journal.cn/Jwk_hxxb/CN/abstract/abstract340458.shtml)里面都谈过。不过近年来也有人提出基于IAO (Intrinsic Atomic Orbitals)代替原始基函数来做Pipek-Mezey定域化的做法(得到的轨道称为IBO),这使得弥散函数不会对结果有什么影响,因为IAO的意义和NBO框架中的极小集NAO很相似,从原始基函数变换成IAO的过程中已经把弥散函数的不良影响充分去除了。实际上,将其它不怕弥散函数的布居分析方法(如Hirshfeld、Becke)结合Pipek-Mezey方法使用也可以使定域化结果不会因为有弥散函数而受到不良影响。

众所周知,NBO轨道的离域程度非常低,这点很像LMO。实际上从等值面图形上会看到NBO比LMO的定域程度其实更高,但NBO的占据数却不像LMO那样是整数,因此NBO轨道不能算做LMO中的一种。实际上,将高占据的NBO和低占据的NBO混合成前述的NLMO之后,就会发现NLMO比对应的NBO的离域程度要高一些,而占据数也变成了LMO要求的整数。产生NBO的过程和Edmiston-Ruedenberg、Boys、Pipek-Mezey定域化方法有本质的不同,NBO轨道是搜索出来的。对于寻找双中心、三中心轨道,NBO程序会循环不同的两个或三个原子组合,检验对应的块密度矩阵的本征值是否超过阈值来判断算不算找到了它们之间的成键轨道,NBO程序里的具体实现其实比这里说的复杂得多,有兴趣者可以看看NBO分析资料汇总里的文章:http://bbs.keinsci.com/forum.php?mod=viewthread&tid=102。需要特别强调的是,NBO的这种搜索方式非常不natural,对于成键方式典型的有机体系还行,但是对于电子结构复杂的体系严重不科学!试图用NBO讨论这些体系的电子结构特征完全是坑爹,经常会得到离谱的结论!有很多人拿着NBO输出的荒谬的结果还各种胡乱讨论一番,纯粹是以讹传讹。后文我们对比LMO和NBO轨道就会明显看到这一点。

Multiwfn支持的AdNDP分析方法在《使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键》(http://sobereva.com/138)一文中做过充分介绍,强烈建议作者阅读,本文不再详述。AdNDP所做的事也是搜索出离域程度比MO更低的轨道,从而便于展现和讨论体系中不同区域共享电子的情况。AdNDP既可以搜索出NBO那样1~3中心轨道,也可以搜索出更多中心轨道,对这类轨道我管它叫半定域化轨道。AdNDP轨道占据数和NBO轨道一样也不是严格整数,因此不属于LMO的一种。下文的例子中,我们也会对比Pipek-Mezey方法产生的LMO与AdNDP分析结果的异同。

总的来说,按照轨道离域程度的上限由低到高,有这样一个关系:NBO-LMO-AdNDP-MO。之所以说“上限”,是因为即便MO也并非都是离域的。比如算个水,就会看到能量最低的MO完全定域在氧上。

虽然LMO、NBO、AdNDP轨道都不是Fock或KS算符的本征函数,但是我们在单电子近似框架下仍可以讨论这些轨道的能量,也就是把Fock或KS矩阵从原始基函数表象变换到相应轨道的表象下,此时对角元就是各个轨道能量,例如第i轨道的能量就是第i对角元,表达式是<φi|f|φi>,这里f代表Fock或KS算符。我们做NBO分析时看到的程序输出的NAO、NBO、NLMO等轨道的能量,就是NBO对量化程序产生的Fock/KS矩阵做这么一个简单的表象变换得到的。


2 在Multiwfn中做分子轨道定域化

Multiwfn从3.4版开始支持轨道定域化方法方法,到了3.5版已经支持基于Mulliken和Lowdin布居的Pipek-Mezey方法,以及Foster-Boys方法。其中PM-Mulliken方法速度最快(对于只定域化占据轨道,处理200个原子一般都没问题),所以是默认的,但是它和PM-Lowdin都不适合有弥散函数的情况,此时应该改为Foster-Boys方法。不过Foster-Boys方法对于中、大体系会比Pipek-Mezey方法慢两三个数量级,而且收敛更慢,因此没有弥散函数时不建议使用。更多介绍和讨论参见手册3.21节。

Multiwfn可以在其主页http://sobereva.com/multiwfn免费下载。不熟悉此程序的请参看《Multiwfn入门tips》(http://sobereva.com/167)、《Multiwfn波函数分析程序的意义、功能与用途》(http://sobereva.com/184)。

把含有基函数信息的输入文件(如fch、molden等,不清楚的话看《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》http://sobereva.com/379)载入到Multiwfn里,然后进入主功能19就可以做轨道定域化。可以选择只定域化占据轨道,还是把占据轨道和空轨道分别定域化。对于非限制性开壳层体系,程序会把alpha和beta部分分别做定域化。此功能不支持限制性开壳层波函数。

定域化过程中会显示迭代细节,当Delta P降到设定的阈值以下时即宣告收敛。默认的阈值,以及默认的循环次数上限在这个功能的界面里会直接看到,可以自己设,一般不需要改。

默认设定下,当迭代收敛后,程序会通过非常便宜的SCPA轨道成分计算方法(不熟悉者《谈谈轨道成份的计算方法》http://sobereva.com/131)计算各个定域化轨道的成份,然后对轨道特征进行指认,判断是单中心、双中心还是更多中心的轨道。然后程序自动把波函数导出到当前目录下的new.fch中。默认情况下,程序会再把这个文件载入,此时内存里的轨道波函数就变成了LMO了,可以直接做进一步分析,如考察轨道图形、计算轨道成份等。注意此时LMO的占据数是对的,但能量还是原先MO的能量。

如果要得到LMO的能量,应当在轨道定域化界面里选择-4 If calculate and print orbital energies将状态切换到Yes,此时Multiwfn就会让你输入一个含有Fock/KS矩阵元数据的文本文件,Multiwfn会从中读取矩阵,并产按照如前所述的方式计算LMO的能量,并且把轨道按照能量由低到高排序储存。输入的文件的格式要求在手册3.21节已经说明了。对于Gaussian用户,可以用IOp(5/33=3)把Fock/KS矩阵提出来,但是每轮SCF迭代都会输出一堆矩阵,输出文件会很大。最简单的方式是把Fock/KS矩阵从Gaussian产生的.47文件(GENNBO输入文件)中提取出来。具体做法是写pop=nboread,然后输入文件末尾空一行写比如$NBO archive file=C:\YOSORO $END,则计算完成后就会产生C:\YOSORO.47文件,把其中$FOCK到$END之间的所有数据拷到一个文本文件里,这部分信息就是Fock/KS矩阵按照下三角矩阵顺序记录的矩阵元,把这个文件作为给Multiwfn的输入就行了。

下面我们通过几个实际体系,演示一下Multiwfn轨道定域化的使用,并对结果进行讨论,与此同时还会与NBO或AdNDP轨道进行一些对比。在《使用Multiwfn通过LOBA方法计算氧化态》(http://sobereva.com/362)这篇文章中使用的LOBA方法是基于定域化轨道的计算氧化态的方法,文中已经利用到了Multiwfn的这个轨道定域化功能,建议大家也看看。

下文涉及的各种文件都可以在此处下载: orbloc.rar (6.22 MB, 下载次数: 25)

评分

参与人数 3eV +15 收起 理由
winterzen + 5 好物!
librakitty + 5 谢谢!
小范范1989 + 5 sob老师又放一个大招。赞

查看全部评分

北京科音自然科学研究中心:http://www.keinsci.com  不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训
思想家公社的门口Blog:http://sobereva.com  发布大量原创计算化学相关博文
Multiwfn量子化学波函数分析程序主页:http://sobereva.com/multiwfn
计算化学公社论坛:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加。

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

2294

帖子

9

威望

3549

eV
积分
6023

Level 6 (一方通行)

首席卖萌官

发表于 2017-6-5 15:51:20 | 显示全部楼层
大博士连续出大招,需仰视得见
尔曹身与名俱灭,不废江河万古流

121

帖子

0

威望

1727

eV
积分
1848

Level 5 (御坂)

发表于 2017-6-6 18:57:12 | 显示全部楼层
谢谢老师分享

107

帖子

1

威望

1100

eV
积分
1227

Level 4 (黑子)

发表于 2017-6-14 03:49:23 | 显示全部楼层
"Pipek-Mezey方法虽然对小体系速度很快,几十原子体系耗时也不高,但是超过50原子体系就挺慢了"

我曾经也做过大概200个原子的体系的PM和Boys局域化的对比,确实PM要慢很多,可是从理论上说两者都是N^3的标度,请问sob老师PM比Boys慢的原因是什么?

150

帖子

3

威望

1119

eV
积分
1329

Level 4 (黑子)

发表于 2018-2-19 08:49:42 | 显示全部楼层
社长那部分PM的代码写的是4次方的。。。。需要把一些量提前计算存储才是3次。。。。

1万

帖子

25

威望

1万

eV
积分
31051

管理员

公社社长

 楼主| 发表于 2018-2-19 10:10:50 | 显示全部楼层
Warm_Cloud 发表于 2018-2-19 08:49
社长那部分PM的代码写的是4次方的。。。。需要把一些量提前计算存储才是3次。。。。

哪些量?
北京科音自然科学研究中心:http://www.keinsci.com  不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训
思想家公社的门口Blog:http://sobereva.com  发布大量原创计算化学相关博文
Multiwfn量子化学波函数分析程序主页:http://sobereva.com/multiwfn
计算化学公社论坛:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加。

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

150

帖子

3

威望

1119

eV
积分
1329

Level 4 (黑子)

发表于 2018-2-19 10:29:57 | 显示全部楼层
!------ Calculate Q value used in Pipek-Mezey localization with Mulliken population
!Qval is also used for determining convergence for other localization methods
real*8 function Qval(Cmat,imo,jmo,iatm)
use defvar
real*8 Cmat(nbasis,nbasis)
integer imo,jmo,iatm,ibas
Qval=0
do ibas=basstart(iatm),basend(iatm)
        Qval=Qval+sum( (Cmat(:,imo)*Cmat(ibas,jmo)+Cmat(ibas,imo)*Cmat(:,jmo)) *Sbas(ibas,:) )
end do
Qval=Qval/2D0
end function
这个地方社长你是现用现算的,虽然是一重循环,但是实际上是2次方,里面的数组操作计算量很大。
计算这里的时候:

用矩阵做分部乘法就可以到3次。

1万

帖子

25

威望

1万

eV
积分
31051

管理员

公社社长

 楼主| 发表于 2018-2-19 11:08:47 | 显示全部楼层
Warm_Cloud 发表于 2018-2-19 10:29
!------ Calculate Q value used in Pipek-Mezey localization with Mulliken population
!Qval is also u ...


你的图我看不到,图是本地路径
我晓得你的意思
PM-lowdin是完美的三次方标度。有了这个PM-Mullliken就再无用武之地了,直接写
      Qij=sum(Cmat(istart:iend,imo)*Cmat(istart:iend,jmo))
      Qii=sum(Cmat(istart:iend,imo)*Cmat(istart:iend,imo))
      Qjj=sum(Cmat(istart:iend,jmo)*Cmat(istart:iend,jmo))
就够了


北京科音自然科学研究中心:http://www.keinsci.com  不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训
思想家公社的门口Blog:http://sobereva.com  发布大量原创计算化学相关博文
Multiwfn量子化学波函数分析程序主页:http://sobereva.com/multiwfn
计算化学公社论坛:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加。

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

150

帖子

3

威望

1119

eV
积分
1329

Level 4 (黑子)

发表于 2018-2-19 11:23:15 | 显示全部楼层
sobereva 发表于 2018-2-19 11:08
你的图我看不到,图是本地路径
我晓得你的意思
PM-lowdin是完美的三次方标度。有了这个PM-Mullliken ...

哈哈,这倒是,我这边的情况是boys的比pm难收敛。

1万

帖子

25

威望

1万

eV
积分
31051

管理员

公社社长

 楼主| 发表于 2018-2-19 12:08:24 | 显示全部楼层
Warm_Cloud 发表于 2018-2-19 11:23
哈哈,这倒是,我这边的情况是boys的比pm难收敛。

是这样的,特别是空轨道
北京科音自然科学研究中心:http://www.keinsci.com  不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训
思想家公社的门口Blog:http://sobereva.com  发布大量原创计算化学相关博文
Multiwfn量子化学波函数分析程序主页:http://sobereva.com/multiwfn
计算化学公社论坛:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加。

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

1万

帖子

25

威望

1万

eV
积分
31051

管理员

公社社长

 楼主| 发表于 2018-2-20 12:16:24 | 显示全部楼层
对这个功能又做了改进,轨道定域化完直接会输出每个轨道的主要特征,明显便于找感兴趣的轨道
无标题.png

另外,鉴于比较常用,新版本里轨道定域化功能作为主功能19出现了

北京科音自然科学研究中心:http://www.keinsci.com  不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训
思想家公社的门口Blog:http://sobereva.com  发布大量原创计算化学相关博文
Multiwfn量子化学波函数分析程序主页:http://sobereva.com/multiwfn
计算化学公社论坛:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加。

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

本版积分规则

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

GMT+8, 2018-5-27 15:55 , Processed in 0.276918 second(s), 28 queries .

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