使用Multiwfn计算(超)极化率密度

使用Multiwfn计算(超)极化率密度

文/Sobereva(3)

First release:2015-Sep-13   Last update:2015-Nov-12



前一阵子有国际友人问怎么用Multiwfn算第二超极化率密度,并且提供了一篇算这个的文献Chem. Phys. Lett., 254, 158-164 (1996)。另外J. Comput. Chem., 26, 1543这篇文章也专门讨论此问题。实际上用Multiwfn的函数自定义运算功能直接就能计算,都不需要任何功能的扩展,这里就说说什么是(超)极化率密度,以及具体怎么去计算。将超极化率密度作出图来放到论文里,可以给(超)极化率计算文章增色不少,讨论也能够更深入。本文用的是2015年9月13日更新的Multiwfn 3.3.8(dev)版,可在http://sobereva.com/multiwfn上免费下载,入门必读见《Multiwfn入门tips》(http://sobereva.com/167)。计算程序使用Gaussian09W D.01。

本文先介绍基本概念,然后用个很简单的小分子甲醛来演示怎么计算第二超极化率密度。


1 基本概念


不了解(超)极化率及其在Gaussian中的计算的话可以先看这里的简单介绍《使用Multiwfn分析Gaussian的极化率、第一超极化率的输出》(http://sobereva.com/231)。

对体系偶极矩向外电场进行Taylor展开得下式



对电子密度函数也可以向外电场进行Taylor展开,其中ρ(n)体现出电子密度对外场的n阶响应。其中ρ(0)就相当于无外场时的电子密度分布。



基于偶极矩(电子贡献的部分)与电子密度的关系,进行对比我们就知道各阶(超)极化率和ρ(n)是怎么对应的。通过ρ(n)可以考察体系不同区域对(超)极化率数值的贡献,这被称为(超)极化率密度。



计算不同的(超)极化率密度过程大同小异,本文我们只讨论计算和绘制上面列出的阶数最高的,即第二超极化率密度。搞明白这个后绘制极化率密度、第一超极化率密度也就是小菜了,举一反三即可,计算公式在文末给出了。

第二超极化率密度ρ(3)是三阶张量函数,第二超极化率γ是个四阶张量,我们不可能把每个分量都讨论一遍,否则太费劲。我们这里只关注最重要的一个分量,zzzz。具体写出来是这样:



我们一旦把ρ(3)zzz这个函数图形绘制出来,就能考察体系不同区域对γzzzz的贡献了。计算这个密度对外场z分量的三阶导数我们需要用有限差分的方法。具体计算公式如下所示

图中开头的是有限差分求二阶导数的公式,我们将它代入有限差分求一阶导数的公式就得到了有限差分求三阶导数的公式。最后的红框就是我们实际通过做密度差以及密度和求ρ(3)zzz的公式了。括号里的Fz代表在z方向外加微量电场(数值大小对应差分步长)时的结果。差分步长过大过小都会造成结果精度的下降,根据测试用0.003 a.u.左右比较不错。

利用Multiwfn的实空间函数自定义操作功能,结合Gaussian等程序产生的波函数文件,就可以很容易地根据上式得到ρ(3)zzz图像了。下面我们用甲醛作为实例来演示计算。


2 计算γzzzz


讨论ρ(3)zzz前,我们先用Gaussian把γzzzz给算出来,当然这一步不是必需的。既然讨论的是分量,就必须把分子朝向先明确了。本文用的示例分子甲醛的朝向如下,结构已在B3LYP/TZVP下优化。



算γ的输入文件如下。这会基于HF的解析三阶导数做一次有限差分得到γ对应的四阶导数。

# HF/aug-cc-pVTZ polar=gamma

B3LYP/TZVP opted

0 1
 C                  0.00000000    0.00000000   -0.52710800
 H                  0.00000000    0.93885600   -1.11413900
 H                  0.00000000   -0.93885600   -1.11413900
 O                  0.00000000    0.00000000    0.67386600

0.0

末尾的0.0代表外场为0,即我们计算的是静态的γ。显然用HF算超极化率肯定不准,但本文只是示意而已所以无所谓。从结果末尾可见γ的zzzz分量为0.145571D+04 a.u.,即1455.71 a.u.。


3 产生波函数文件


下面我们获得不同大小外电场下的波函数文件。

虽然一般加电场计算的时候都要用nosymm避免Gaussian把分子自动调整到标准朝向下(否则实际外加电场方向将和期望的不符),但当前例子里的坐标已经是在标准朝向下了,所以不用写nosymm。

我们总共要得到四个波函数文件,分别是外场为-2Fz、-Fz、Fz、2Fz的情况。这里我们用0.003a.u.外场为步长。我们先计算-2Fz,于是外场就是-0.006a.u.,对应field=z+60关键词(提醒一下,Gaussian中field关键词定义的外场方向和一般定义是反过来的)。输入文件如下:

# HF/aug-cc-pVTZ field=Z+60 out=wfx

B3LYP/TZVP opted

0 1
 C                  0.00000000    0.00000000   -0.52710800
 H                  0.00000000    0.93885600   -1.11413900
 H                  0.00000000   -0.93885600   -1.11413900
 O                  0.00000000    0.00000000    0.67386600

C:\H2CO-2FZ.wfx

算完上面这个任务后,改成field=Z+30再次计算得到H2CO-FZ.wfx,改成field=Z-30得到H2CO+FZ.wfx,改成field=Z-60得到H2CO+2FZ.wfx。这样所需的四个外场下的波函数文件就都齐了。

实际上Gaussian能产生的.wfn、.wfx、.fch都可以作为Multiwfn的输入文件,结果都一样。但是当外场很小时(特别是小到0.001 a.u.程度时),可能是因为.wfn格式的数据记录有效位数相对有限的原因,结果和.wfx偏差较大,所以这里用数值记录精度更高的.wfx格式。


4 通过自定义运算获得ρ(3)zzz


Multiwfn的主功能3、4、5,即绘制曲线图、平面图、等值面图的功能都可以设定对实空间函数的自定义运算,作密度差是其最常见的应用,有兴趣可以看《使用Multiwfn作电子密度差图》(http://sobereva.com/113)。本文算ρ(3)zzz是自定义运算的一个很特殊也很能证明其价值的应用。我们按照ρ(3)zzz计算公式的分子定义自定义操作的运算顺序,然后除以分母就行了。

启动Multiwfn,然后依次输入
c:\H2CO+2FZ.wfx   //首先载入ρ(3)zzz计算公式的分子中的第一项ρ(2Fz)
5      //计算格点数据并显示等值面
0      //自定义操作
5      //接下来对5个文件进行操作
-,c:\H2CO+FZ.wfx     //公式中的-2ρ(Fz)项,通过减两次来表现
-,c:\H2CO+FZ.wfx
+,c:\H2CO-FZ.wfx     //公式中的+2ρ(Fz)项,通过加两次来表现
+,c:\H2CO-FZ.wfx
-,c:\H2CO-2FZ.wfx    //公式中的ρ(-2Fz)项
1    //电子密度
2    //中等质量格点
然后程序开始计算密度数据并按照上述定义的规则进行运算,算完之后,我们要让格点数据除以分母,即2(Fz)^3项。对于当前情况,这一项是2*0.003^3=0.000000054。接着输入
6    //除以某个值
0.000000054
-1   //观看等值面
在蹦出的图形窗口中将isovalue设大一些,这里设为0.5,如下所示



这就是ρ(3)zzz啦!第一节已经给出了ρ(3)zzz和γzzzz的关系,即各个位置的ρ(3)zzz乘上相应位置的-Z再全空间积分就是γzzzz,因此我们容易想象体系不同位置对γzzzz的贡献。比如在原点附近(大约C=O中央),Z接近于0,所以这部分区域即便等值面挺大,即ρ(3)zzz数值不小,在乘了-Z后对γzzzz的贡献也不会太大。而O的上方有一大块蓝色(负值)等值面,Z的数值也不小,和-Z相乘后肯定那部分电子对γzzzz有较大正贡献。在两个H下方有一大块绿色(正值)区域,那部分Z为不小的负值,故乘上-Z后也对γzzzz有较大正贡献。

这么考察还稍微有点费劲,我们可以把已经得到的这个ρ(3)zzz乘上-Z之后再考察图形,这样就能更直观看到不同位置的贡献了,因为直接一积分就是γzzzz了。我们可以用Multiwfn的格点数据处理功能实现。我们关闭图形窗口,然后输入
0    //返回主菜单
13   //处理格点数据
11   //对格点数据进行操作
20   //乘上坐标变量(此功能在2015年9月13日及之后更新的Multiwfn版本中才有)
Z    //乘上Z坐标
11
5    //乘上个常数
-1   //乘以-1。此时格点数据就相当于-Z*ρ(3)zzz的格点数据了
-2   //观看等值面
数值为1.0的等值面如下所示



这下考察起来更容易了,也就是图中绿色等值面越大的区域对γzzzz的贡献越正,蓝色等值面越大的区域对γzzzz的贡献越负。可见,在分子价层区域,基本上都是负贡献,而在O、H外围区域有着很大的正贡献,并且显然明显超越了负值区域,这是为什么当前体系的γzzzz是一个很大正值(前面计算结果为1455.71 a.u.)。对于一个较大体系,通过这么张图,立刻就能解释清楚γzzzz为什么大或者为什么小,为什么正或者为什么负,各个区域贡献一目了然,太爽啦!放到论文里自然会加分。

这张图也引出一个很重要的问题,也就是计算γ的基组一定要有充足的弥散函数!!!如果没有弥散函数,那么可以想象,上图中分子边缘弥散程度较高的绿色区域就没法合理描述,此时γzzzz肯定会被低估,甚至为负值!为了证明这点,我们在cc-pVDZ下进行计算,结果仅为26.4 a.u.,而在6-31G*下计算,结果为-108.85 a.u.,连符号都错了。而加入弥散函数成为6-31++G*后,结果为726.27 a.u.,虽然定量上很烂,但趋势基本对了。

不同的基组下、不同的理论方法下,γzzzz不同,显然-Z*ρ(3)zzz也相应地不同,因为后者全空间积分就等于前者。因此,我们也可以通过比较-Z*ρ(3)zzz图来考察为什么不同基组、不同的理论方法下结果有异。我们还可以用一个高精度级别下计算的-Z*ρ(3)zzz和一个较烂级别下计算的-Z*ρ(3)zzz作差值图(用Multiwfn主功能13的子功能11做两个格点数据间的相减运算即可),马上就能知道较烂的级别是因为在什么区域描述不好而导致结果不好。此时误差就不再仅仅是一个数值了,我们可以洞察内在真相。

-Z*ρ(3)zzz的积分是γzzzz,我们也可以很容易地通过数值验证这一点。关闭刚才的图形界面,然后进入选项17,再输入1,就得到了当前格点数据在全空间中的统计信息,其中
Integral of all data:                    1429.5133014493
就是基于立方格点积分的结果。这个值和我们第二节通过导数方法得到的γzzzz值1455.71 a.u.很接近。如果我们之前在计算格点数据时在设定格点的界面中将延展距离设大到10 Bohr,并且用更好的格点(high quality grid),则结果为1459.43 a.u.,非常接近导数方法计算的结果了。


5 绘制ρ(3)zzz和-Z*ρ(3)zzz的平面图


我们也可以利用Multiwfn的主功能4十分容易地绘制ρ(3)zzz和-Z*ρ(3)zzz的平面图,自定义运算的设置方式完全一样。我们先绘制甲醛分子平面上的ρ(3)zzz等值线图。

c:\H2CO+2FZ.wfx
4      //绘制平面图
0      //自定义操作
5
-,c:\H2CO+FZ.wfx
-,c:\H2CO+FZ.wfx
+,c:\H2CO-FZ.wfx
+,c:\H2CO-FZ.wfx
-,c:\H2CO-2FZ.wfx
1   //电子密度
2   //等值线图
[回车]   //用默认的格点设定
0   //调整延展距离
10   //设为10 Bohr
3   //YZ平面
0   //X=0
关闭图像,然后选
-7  //将数据乘上个值
18518518.5185185   //即分母0.000000054的倒数
-1  //重新绘图

结果如下,实线和虚线分别代表正值和负值部分。后处理菜单还有大量选项可以调节绘图效果。



再来绘制-Z*ρ(3)zzz的等值线图。关掉Multiwfn,把settings.ini的iuserfunc设为23,因为从手册2.6节可以看到设为23时用户自定义函数等价于-Z*ρ(r),因此我们直接对用户自定义函数进行自定义操作就直接得到了-Z*ρ(3)zzz。绘图步骤和上面一样,唯一不同的是选择实空间函数的时候选100(用户自定义函数)。结果如下



也可以作成填色图+等值线图,绘图类型选绘制填色图然后后处理菜单再把等值线打开就行了:


附录


评论已关闭