使用Multiwfn分析Gaussian09的极化率、第一超极化率的输出
文/Sobereva(3) 2014-Apr-27
1 前言
Gaussian中polar关键词专门用来计算极化率α和第一超极化率β的(不能直接计算第二超极化率γ)。Gaussian输出的很多信息总是令初学者很难看懂,(超)极化率的输出更是如此,内容杂乱,而且不同理论方法下输出的信息格式还差异很大,分析起来很不方便。Multiwfn的主功能200里的子功能7就是专门用来解析Gaussian09的(超)极化率的输出的,使之简洁易懂。与此同时,Multiwfn还会计算一些对于分析(超)极化率十分重要的量,比如<α>、α的各向异性值、β在分子偶极矩方向的分量等等,这使得通过Gaussian研究(超)极化率变得方便得多。此文的目的就是介绍如何使用Multiwfn分析Gaussian09的(超)极化率的输出。
本文使用的Multiwfn为3.3.1版,可以在其主页http://Multiwfn.codeplex.com上免费下载,不要用老版本,老版本无此功能。Multiwfn是最为强大的波函数分析程序,此文介绍的只不过是它的一个附属功能而已。如果对Multiwfn不熟悉,建议阅读《Multiwfn波函数分析程序的意义、功能与用途》(http://hi.baidu.com/sobereva/item/b35581418db96cea1f19bc09)和《Multiwfn入门tips》(http://hi.baidu.com/sobereva/item/896ee3a19f6d7d3c020a4d76)。Multiwfn分析polar关键词输出的这个功能只支持G09,不支持G03。
2 Gaussian中(超)极化率的计算
这一节介绍在Gaussian中计算(超)极化率的基本原理和所需要的关键词。
将能量E对均匀外电场F进行Taylor展开得到如下式子
μ0就是分子的永久偶极矩矢量(无外场下的偶极矩)。α是分子的极化率,也称线性光学系数,是个3*3二阶张量。β是第一超极化率,也称分子的二阶非线性光学(NLO)系数,是个3*3*3的三阶张量。γ是第二超极化率,亦称分子的三阶NLO系数,研究得相对少一些。而更高的δ极少被讨论。(超)极化率数值和电场频率是相关的,外场频率为0的情况被称为静态(超)极化率;如果外场频率不是0,比如某一频率的光,对应的就是动态(超)极化率,或者叫含频(超)极化率。
Gaussian的polar关键词求这些量是通过导数方法,其原理很明确。如上式可见,求极化率要令能量对外场做两次导数,求第一超极化率要做三次导数。
求能量的导数分为三种情况:
(1)解析导数。这种方法求导数又快又精确,但是编程实现很困难,特别是对于高级的后HF而言。这种方式求能量导数时需要用到分子轨道系数对外场的导数,这通过求解CPHF方程来实现。并且利用含频CPHF方程,可以得到动态(超)极化率。这种求(超)极化率的方式也被称为CPHF方法。
(2)数值导数。这种方式通过有限差分方式得到能量的导数,也称为有限场(FF)方法。有限差分在所有数值算法书里都会介绍,比如求一个函数f在0处的导数可表示为f(0)=[f(0.001)-f(-0.001)/(2*0.001)。通过这种方法求(超)极化率非常费时间,需要在不同数值、不同方向的外场下计算很多很多次单点,而且求高阶导数精度很差,因为每求一次数值导数都会由于数值噪音而累积误差,随着求导阶数增加误差迅速放大,做二阶以上导数时精度就很烂了(除非通过额外修正、外推等特殊方式处理)。另外,通过数值导数也没法求得动态(超)极化率。
(3) 解析+数值混合导数。这就是基于低阶解析导数,再通过一次或多次有限差分求高阶导数以获得(超)极化率的方法。精度和速度介于全解析导数和全数值导数之间。
G09在支持解析导数方面已经算是很好的,对于HF、DFT、半经验方法都能支持到三阶解析导数,因此可以全解析地得到静态和含频的α和β(而其它量化程序都很难对HF/DFT/半经验支持到这么高阶导数,通常也就能支持到二阶)。对这类支持三阶解析导数的方法,使用polar关键词就会直接给出α和β的结果。如果同时写上CPHF=RdFreq,并且在输入文件末尾空一行写上外场频率w,比如0.05 0.12 0.3(静态的情况总是会被计算,故不需要写0.0),就会计算含频极化率α(-ω;ω)和含频第一极化率β(-ω;ω,0)。如果要算的含频第一极化率是SHG形式,即β(-2ω;ω,ω),就写polar=DCSHG即可,此时CPHF=RdFreq可以忽略不写。
对于MP2等G09只能支持到二阶解析导数的方法,可以全解析地计算α,但是若想计算β,就得再做一次有限差分才行。G09中,此时直接写polar只会计算α,若还想得到β就得写polar=Cubic,会自动基于解析二阶导数做有限差分来得到三阶导数。
对于CISD、QCISD、CCSD、MP3、MP4(SDQ)等G09只支持一阶解析导数的方法,直接写polar关键词会基于解析一阶导数做一次有限差分得到二阶导数,即α。若还想得到β就得写polar=DoubleNumer(等价于polar=enonly),会基于解析一阶导数做两次有限差分得到三阶导数。
对于CCSD(T)、QCISD(T)、MP4(SDTQ)、MP5等G09只支持计算单点能的方法,写polar关键词会对能量做两次有限差分并由此得到α。对于它们在G09里没有办法直接给出β。
注意对于HF、DFT、半经验以外的方法,G09都没法给出含频(超)极化率。
若想得到更高阶导数计算更高阶超极化率,比如DFT下计算第二超极化率γ,需要基于其解析三阶导数做一次有限差分,但G09没有关键词能自动去做,因此就得自己手写或者通过自编脚本生成一批在不同外电场下的Gaussian输入文件来实现了。
值得一提的是,利用Multiwfn,可以通过完全态求和方法基于Gaussian的CIS/TDHF/TDDFT计算的输出得到静态或动态的极化率和第一、二、三超极化率,见《使用Multiwfn基于完全态求和(SOS)方法计算极化率和超极化率》(http://hi.baidu.com/sobereva/item/f4960b3314d400be134b1480)。
3 使用Multiwfn解析polar关键词的输出
对于上一节讨论的各种情况,polar关键词输出的信息都能被Multiwfn所解析。注意route section一定要写上#P,否则输出文件将无法被正确解析。
启动Multiwfn后,先输入G09的polar任务的输出文件的路径,比如c:\CH3NH2_polar.out,然后进入主功能200,选7,之后会看到下面的菜单:
-3 Set the unit in the output, current: a.u.
-2 Set the output destination, current: Output to screen
-1 Toggle if load frequency-dependent result for option 1, current: No
0 Return
1 "Polar" + analytic 3-order deriv. (HF/DFT/Semi-empirical)
2 "Polar" + analytic 2-order deriv. (MP2...)
3 "Polar=Cubic" + analytic 2-order deriv.
4 "Polar" + analytic 1-order deriv. (CISD,QCISD,CCSD,MP3,MP4(SDQ)...)
5 "Polar=DoubleNumer" or "Polar=EnOnly" + analytic 1-order deriv.
6 "Polar" + energy only (CCSD(T),QCISD(T),MP4(SDTQ),MP5...)
默认情况下输出的单位是a.u.,如果选-3,也可以切换为SI或esu单位。
默认情况下输出信息会显示到屏幕上,也可以选-2来更改为输出到当前目录下的polar.txt当中。
如果在HF/DFT/半经验计算时用了polar CPHF=rdfreq或polar=DCSHG做了含频(超)极化率计算,默认情况下只会解析静态的α和β。如果想让Multiwfn分析含频的α和β,应当先选-1。注意如果在Gaussian输入文件里定义了多个频率值,如0.1,0.2,0.3,则Multiwfn只会去解析ω=0.1下的结果。
选项1~6对应上一节涉及的6种情况,当前你在什么情况下计算的就选哪种,然后Multiwfn就会输出分析的结果。比如你用了CCSD结合polar=doublenumer就应当选5。之所以要对此进行选择是因为不同情况下G09的输出格式相差极大,Multiwfn必须知道你在什么情况下计算的,才能相应地对输出内容进行解析。如前面所讨论的,这6种情况中只有1、3、5才会给出β,其它的只给出α。
3.1 静态(超)极化率一例
下面就以分析甲胺在#P PBE1PBE/aug-cc-pVTZ polar关键词下得到的输出文件为例进行说明。
当前例子用的是PBE1PBE,是DFT中常用的一种泛函,故对应的是前述第一种情况,因此当我们看到上面的菜单时应选1,然后马上会看到α和β的分析结果,下面将对各项的含义进行说明。我们先看α部分
Dipole moment:
X,Y,Z= -0.506266 0.150290 0.000000 Norm= 0.528103
Static polarizability:
XX= 25.046200
XY= -0.088300
YY= 28.891100
XZ= 0.000000
YZ= 0.000000
ZZ= 24.435800
Isotropic average polarizability: 26.124367
Polarizability anisotropy (definition 1): 4.186426
Eigenvalues of polarizability tensor: 24.43580 25.04417 28.89313
Polarizability anisotropy (definition 2): 4.153140
如前所述,默认情况下所有输出的量都是以a.u.为单位。
极化率张量是个对称的3*3矩阵,故所以只有6个独立的矩阵元列了出来。每个体系的极化率的大小一般只用一个数值<α>来表示,也就是上面给出各项同性平均极化率,它是对XX、YY、ZZ取平均得到的。极化率并不是各项同性的,体系对不同方向射来的电场的响应是不同的。如上可见,Multiwfn也输出了极化率的各向异性值Δα,数值越大各向异性越强,球对称体系显然其值为0。各向异性可以以两种方式定义,定义1和定义2分别如下所示
ε代表极化率张量的本征值,它也被Multiwfn输出了出来。通常衡量极化率的各向异性程度用的是定义1,而定义2的形式一般是在衡量磁屏蔽张量各向异性时用的。
下面是第一超极化率部分的解析结果
Static first hyperpolarizability:
XXX= 13.173100
XXY= 4.937890
XYY= 13.647700
YYY= -35.484000
XXZ= 0.000000
XYZ= 0.000000
YYZ= 0.000000
XZZ= -2.763570
YZZ= 1.342190
ZZZ= 0.000000
Beta_X= 24.05723 Beta_Y= -29.20392 Beta_Z= 0.00000
Magnitude of first hyperpolarizability: 37.836745
Projection of beta on dipole moment: -31.373491
Beta || : -18.824094
Beta ||(z) : 0.000000
Beta _|_(z) : 0.000000
Gaussian的β的输出有个众所周知的问题,也就是符号是反的,即每个β元素的数值都应该乘以-1。有文章专门指出了这点,见《Gaussian程序计算的一阶超极化率的符号问题》(四川师范大学学报(自然科学版),33,228)。这个问题在Multiwfn中已经被考虑进去了,Multiwfn输出的β的符号都是正确的。
β有3*3*3=27个元素。静态β满足Kleinman对称性,三个标号可以随意置换而完全不影响结果,比如xyy=yxy=yyx,只有10个元素是唯一的,所以Multiwfn也只输出了10个元素。
Beta_X、Beta_Y、Beta_Z衡量的是β在X、Y、Z方向上的分量。Magnitude of first hyperpolarizability体现了β的整体的大小。其定义分别如下所示
只有β在偶极矩方向上的投影值才是可以通过电场诱导二次谐波产生实验(EFISH)测量出来的,和实验值对比通常比的都是这个值。Multiwfn给出了β在偶极矩方向上的投影值(β_prj)。另一个文献中常见的量是Beta ||,||代表平行于偶极矩,它正比于β_prj
文献中常见的Beta ||(z)和Beta _|_(z)的定义如下所示,分别衡量的是平行和垂直于Z方向上的β的值
当偶极矩方向和Z方向一致时,Beta_Z=β_prj,Beta ||(z)=Beta ||。
值得一提的是,默认情况下,Gaussian会将输入文件里体系的坐标调整到标准朝向下,这通常会令体系发生旋转。Gaussian计算出的,也即Multiwfn所解析的,都是标准朝向下的结果,因此X、Y、Z可能和你的输入文件里的笛卡尔轴方向不对应。如果想避免自动调整到标准朝向,用nosymm关键词。
3.2 动态(超)极化率一例
此例使用Multiwfn解析polar关键词产生的动态(超)极化率。输入文件如下,电磁场的能量为0.07 a.u.
#p PBE1PBE/aug-cc-pVTZ polar CPHF=RdFreq
test
0 1
C -0.55391731 0.43227932 0.14513431
H -1.29461381 -0.13844703 0.71478821
H -0.40071988 1.37810560 0.67504812
H -0.99858078 0.66864399 -0.83680924
N 0.70741247 -0.31205725 0.11165308
H 0.57479710 -1.18908235 -0.38659422
H 1.39864521 0.20863272 -0.42322026
0.07
令Multiwfn载入输出文件,并进入主功能200里的子功能7后,先选-1切换成解析动态值,然后再选1,结果如下(极化率部分没什么好说的,略过)
Frequency-dependent first hyperpolarizability beta(-w;w,0)
XXX= 16.334900
XYX= 7.369050
YYX= 16.593300
XZX= 0.000000
YZX= 0.000000
ZZX= -2.918320
XXY= 8.947060
XYY= 15.473800
YYY= -40.341800
XZY= 0.000000
YZY= 0.000000
ZZY= 1.240500
XXZ= 0.000000
XYZ= 0.000000
YYZ= 0.000000
XZZ= -2.796870
YZZ= 1.609830
ZZZ= 0.000000
Beta_X= 29.34451 Beta_Y= -30.96003 Beta_Z= 0.00000
Magnitude of first hyperpolarizability: 42.657048
Projection of beta on dipole moment: -36.941911
Beta || : -22.165146
Beta ||(z) : 0.000000
Beta _|_(z) : 0.000000
和上一例相比,可见在含频外场下β值增大了。各项的含义和上例一致,不再累述。值得一提的是这次Multiwfn输出了beta(-ω;ω,0)的18个元素,而不是上例那样10个。原因在于,对于beta(-ω;ω,0)这种情况,只有i,j,k三个标号中的i,j可以互换,比如xyy=yxy≠yyx,张量的对称性没有静态β的那么高,故有18个独立的元素。不过,当外场频率不是很大时,Kleinman对称性是可以近似满足的,从上面的结果可见,诸如xyx和xxy相差不算很大。 |