使用Multiwfn分析Gaussian09的极化率、第一超极化率的输出

使用Multiwfn分析Gaussian09的极化率、第一超极化率的输出

文/Sobereva(3)
First release: 2014-Apr-27     Last update: 2015-Feb-8


1 前言

Gaussian中polar关键词专门用来计算极化率α和第一超极化率β的(不能直接计算第二超极化率γ)。Gaussian输出的很多信息总是令初学者很难看懂,(超)极化率的输出更是如此,内容杂乱,而且不同理论方法下输出的信息格式还差异很大,分析起来很不方便。Multiwfn的主功能200里的子功能7就是专门用来解析Gaussian09的(超)极化率的输出的,使之简洁易懂。与此同时,Multiwfn还会计算一些对于分析(超)极化率十分重要的量,比如<α>、α的各向异性值、β在分子偶极矩方向的分量等等,这使得通过Gaussian研究(超)极化率变得方便得多。此文的目的就是介绍如何使用Multiwfn分析Gaussian09的(超)极化率的输出。

本文使用的Multiwfn为3.3.7版,可以在其主页http://Multiwfn.codeplex.com上免费下载,不要用老版本,老版本无此功能。Multiwfn是最为强大的波函数分析程序,此文介绍的只不过是它的一个附属功能而已。如果对Multiwfn不熟悉,建议阅读《Multiwfn波函数分析程序的意义、功能与用途》(http://sobereva.com/184)和《Multiwfn入门tips》(http://sobereva.com/167)。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可以忽略不写,并且此时β(-ω;ω,0)也会照样计算。

对于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 D.01开始,对于支持三阶解析导数的方法可以用polar=gamma来算静态γ和动态的γ(-2ω;ω,ω,0)、γ(-ω;ω,0,0),会基于解析三阶导数自动做有限差分来得到所需的四阶导数。末尾应空一行写上外场频率。

值得一提的是,利用Multiwfn,可以通过完全态求和方法基于Gaussian的CIS/TDHF/TDDFT计算的输出得到静态或动态的极化率和第一、二、三超极化率,见《使用Multiwfn基于完全态求和(SOS)方法计算极化率和超极化率》(http://sobereva.com/232)。


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,用户可以选择分析哪个频率下的结果。用户也可以选择分析β(-ω;ω,0)还是β(-2ω;ω,ω),注意只有写了DCSHG时选择β(-2ω;ω,ω)才有意义。

选项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,至少直到G09 D.01版依然存在这个问题。有文章专门指出了这点,见《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,一些文中也写为β_vec)。另一个文献中常见的量是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和0.1 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 0.1


令Multiwfn载入输出文件,并进入主功能200里的子功能7后,先选-1切换成解析动态值,然后再选1。然后程序问载入哪个频率下的结果,此例我们选2,即w=0.07的情况。然后我们得到了w=0.07下的极化率输出。之后程序问分析β(-ω;ω,0)还是β(-2ω;ω,ω)的结果。由于没写DCSHG,我们只能选择β(-ω;ω,0)。输出如下

 Frequency-dependent first hyperpolarizability Beta(-w;w,0)
 XXX=            16.334900
 XYX= YXX=        7.369050
 YYX=            16.593300
 XZX= ZXX=        0.000000
 YZX= ZYX=        0.000000
 ZZX=            -2.918320
 XXY=             8.947060
 XYY= YXY=       15.473800
 YYY=           -40.341800
 XZY= ZXY=        0.000000
 YZY= ZYY=        0.000000
 ZZY=             1.240500
 XXZ=             0.000000
 XYZ= YXZ=        0.000000
 YYZ=             0.000000
 XZZ= ZXZ=       -2.796870
 YZZ= ZYZ=        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

和上一例相比,可见在含频外场下β值增大了。各项的含义和上例一致,不再累述。值得一提的是对于beta(-ω;ω,0)这种情况,只有i,j,k三个标号中的i,j可以互换,比如xyy=yxy≠yyx,张量的对称性没有静态β的那么高,故有18个独立的元素。不过,当外场频率不是很大时,Kleinman对称性是可以近似满足的,从上面的结果可见,诸如xyx和xxy相差不算很大。

假设我们要看w=0.1时的β(-ω;ω,0),那么再次进入此功能,然后选1-3-1即可。

基于当前输出文件没法查看β(-2ω;ω,ω),因为当前例子没写DCSHG。要看它的话应把polar后面加上=DCSHG重算一次,然后再进入此功能,并选择分析β(-2ω;ω,ω)情况的数据。

仅有一条评论

  1. Guga

    请问Sob老师,G09输出的gamma有这个符号的问题吗?

添加新评论