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

计算化学公社

 找回密码
 现在注册!
查看: 2285|回复: 4

[综合交流] 使用NWChem做分数占据数的DFT计算

[复制链接]

2万

帖子

25

威望

2万

eV
积分
46973

管理员

公社社长

发表于 2017-2-25 03:42:35 | 显示全部楼层 |阅读模式
使用NWChem做分数占据数的DFT计算

文/Sobereva @北京科音   2017-Feb-25


1 前言

分数占据数(Fractional occupation number, FON)计算是指被计算的体系有的轨道占据数不为整数,也因此允许整个体系的电子数不为整数。这看似有点超现实,实际中体系电子数不可能为非整数,毕竟电子不可分割,但对于DFT(或HF)计算是可以实现的,而且在程序实现上超简单。大家知道,常规DFT计算的时候都是已知电子数,根据每一轮迭代时解出的KS轨道能量从低往高填充,然后用占据轨道构建用于下一轮迭代用的密度矩阵以及计算体系能量。而FON计算时,有的轨道依然是整数占据,有的则是自设的分数占据数(这有点像自然轨道),构建密度矩阵时考虑所有这些占据数不为0的轨道即可,除此以外不需要对程序做任何修改。FON计算既可以是限制性计算,轨道占据数允许在[0,2]范围,也可以是非限制性计算,alpha/beta轨道占据数允许在[0,1]范围。

做FON计算有一些特别的用处,很重要之一就是考察各种近似的交换-相关泛函存在的自相互作用误差(self-interaction error, SIE)导致的离域性误差(delocalization error)问题。这里引用doi: 10.3866/PKU.WHXB201605301中的一张图,图中ΔN指相对于中性状态时电子数的改变,纵坐标是相对于中性状态时体系能量的变化。
0.png

理论上,如果交换-相关泛函是精确,那么体系能量随电子数的变化应当如图中绿线所示,是精确的折线。但是HF成份较低的泛函存在较大的SIE问题,导致电子离域性过强,所以从右图上可看到PBE、B3LYP等泛函的曲线都往外凸,而且HF成份越低凸得越厉害。HF交换势则完全没有SIE问题,但是在图中看到其曲线往里凹,高估了电子定域性,也不很合理。近年来比较火的w调控方法则可以对范围分离泛函的HF交换势掺入的方式优化到一个较理想的情形,它使得HOMO轨道能量和体系电离能尽可能接近(对于精确泛函二者则严格相同),从而极大程度解决SIE问题,使得对超极化率、中/大体系激发能等方面精度提升很多,图中可见经过w调控的泛函LC-PBE0*的曲线确实相当平直。更多和w调控的讨论见《优化长程校正泛函w参数的简便工具optDFTw》(http://sobereva.com/346)。

如今大量使用w调控的文章都纷纷绘制如上的图像讨论不同泛函的SIE,本文就介绍怎么在NWChem中做FON计算并绘制如上的图像。


2 在NWChem中做分数占据数计算

NWChem是为数不多的能支持FON计算的量化程序,开源免费。编译方法参看《NWChem 6.6编译方法》(http://sobereva.com/270)。

NWChem中支持分数占据数的DFT计算。不是光指定一个总电子数就完了,而必须在DFT段落中用fon关键词来指定占据方式。下面的语句是个用于闭壳层FON计算的设定例子:
fon partial 2 electrons 0.8 filled 5
代表有两个轨道是分数占据的,0.8个电子均分给这两个轨道。另外有5个双占据轨道。因此此体系一共有0.8+5*2=10.8个电子。

对于开壳层体系,需要分别对alpha和beta轨道定义占据数,比如
fon alpha partial 2 electrons 0.8 filled 6
fon beta partial 1 electrons 1 filled 5
ODFT
这代表有6个占据数为1.0的alpha轨道,另有0.8个电子均分给两个alpha轨道。占据的beta轨道有6个,占据数都为1.0(含义上等价于fon beta partial 0 electrons 0 filled 6,但不能这么写否则程序报错)。ODFT必须写,指明做开壳层DFT计算(O代表Open)。

在NWChem中普通的计算都是通过charge关键词指定体系净电荷以及在DFT字段中用mult关键词指定自旋多重度来让程序知道电子是怎么排布的。而做FON的时候我们是自行指定电子排布方式的,因此不需要再写charge和mult了。


3 对吡咯做分数占据数计算

这里我们对吡咯体系做FON计算。吡咯总共36个电子,因此中性状态下alpha和beta轨道各有18个占据轨道。假设我们当前研究体系有36.7个电子的情形,则多出来的0.7个电子应当被视为处在一个alpha轨道上(可以姑且视作是处在中性状态下的alpha LUMO轨道上,但注意36.7个电子计算和36个电子计算得到的轨道形状和能量都是不同的)。B3LYP/6-31G*时输入文件如下,结构是事先在B3LYP/6-31G*下对中性状态优化的。
GEOMETRY
C                  0.00000000    1.12552000    0.33154500
C                  0.00000000    0.71275300   -0.98347600
C                  0.00000000   -0.71275300   -0.98347600
C                  0.00000000   -1.12552000    0.33154500
N                  0.00000000    0.00000000    1.12233100
H                  0.00000000    0.00000000    2.13032900
H                  0.00000000    2.11411400    0.76778100
H                  0.00000000    1.36093100   -1.84951800
H                  0.00000000   -1.36093100   -1.84951800
H                  0.00000000   -2.11411400    0.76778100
END
BASIS
* library 6-31G*
END
DFT
fon alpha partial 1 electrons 0.7 filled 18
fon beta partial 1 electrons 1 filled 17
odft
XC b3lyp
END
TASK DFT ENERGY

输出信息和