使用Multiwfn分析Gaussian的极化率、超极化率的输出

使用Multiwfn分析Gaussian的极化率、超极化率的输出
Using Multiwfn to analyze the polarizability and hyperpolarizability outputted by Gaussian

文/Sobereva @北京科音
First release: 2014-Apr-27     Last update: 2023-Aug-11


1 前言

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

本文内容对应的是Multiwfn官网上的最新版本,可以在其主页http://sobereva.com/multiwfn上免费下载,不要用老版本。Multiwfn是最为强大的波函数分析程序,此文介绍的只不过是它的一个附属功能而已。如果对Multiwfn不熟悉,建议阅读《Multiwfn波函数分析程序的意义、功能与用途》(http://sobereva.com/184)和《Multiwfn入门tips》(http://sobereva.com/167)。Multiwfn分析polar关键词输出的这个功能只支持Gaussian 09及之后版本,不支持更老版本。

Multiwfn的这个功能还支持计算与超瑞利散射(Hyper-Rayleigh Scattering, HRS)实验相关的各种量,在另一篇文章《使用Multiwfn计算与超瑞利散射(HRS)实验相关的量》(http://sobereva.com/499当中做了详细介绍。

如果你使用了本文介绍的Multiwfn的功能,请在发表文章时提及研究中使用了Multiwfn提取Gaussian输出文件中的相关数据并计算了相关的量,并应当根据Multiwfn程序启动时提示的信息恰当引用Multiwfn程序


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) 解析+数值混合导数。这就是基于低阶解析导数,再通过一次或多次有限差分求高阶导数以获得(超)极化率的方法。精度和速度介于全解析导数和全数值导数之间。

Gaussian在支持解析导数方面已经算是很好的,对于HF、DFT、半经验方法都能支持到三阶解析导数,因此可以全解析地得到静态和含频的α和β(而其它量化程序都很难对HF/DFT/半经验支持到这么高阶导数,通常也就能支持到二阶)。对这类支持三阶解析导数的方法,使用polar关键词就会直接给出α和β的结果。如果同时写上CPHF=RdFreq,并且在输入文件末尾空一行写上外场频率w,比如0.05 0.12 0.3或532nm 670nm(静态的情况总是会被计算,故不需要特意写),就会计算含频极化率α(-ω;ω)和含频第一超极化率β(-ω;ω,0)。如果要算的含频第一超极化率是SHG形式,即β(-2ω;ω,ω),就写polar=DCSHG即可,此时CPHF=RdFreq可以忽略不写,并且此时β(-ω;ω,0)也会照样计算。

对于MP2等Gaussian只能支持到二阶解析导数的方法,可以全解析地计算α,但是若想计算β,就得再做一次有限差分才行。Gaussian中,此时直接写polar只会计算α,若还想得到β就得写polar=Cubic,会自动基于解析二阶导数做有限差分来得到三阶导数。

对于CISD、QCISD、CCSD、MP3、MP4(SDQ)等Gaussian只支持一阶解析导数的方法,直接写polar关键词会基于解析一阶导数做一次有限差分得到二阶导数,即α。若还想得到β就得写polar=DoubleNumer(等价于polar=enonly),会基于解析一阶导数做两次有限差分得到三阶导数。

对于CCSD(T)、QCISD(T)、MP4(SDTQ)、MP5等Gaussian只支持计算单点能的方法,写polar关键词会对能量做两次有限差分并由此得到α。对于它们在Gaussian里没有办法直接给出β。

注意对于HF、DFT、半经验以外的方法,Gaussian都没法给出含频(超)极化率。

若想得到更高阶导数计算更高阶超极化率,比如DFT下计算第二超极化率γ,需要基于其解析三阶导数做一次有限差分。从G09 D.01开始,对于支持三阶解析导数的方法可以用polar=gamma来算静态γ和动态的γ(-2ω;ω,ω,0)、γ(-ω;ω,0,0),会基于解析三阶导数自动做有限差分来得到所需的四阶导数。末尾应空一行写上外场频率。


值得一提的是,利用Multiwfn,可以通过完全态求和方法基于Gaussian的CIS/TDHF/TDDFT计算的输出得到静态或动态的极化率和第一、二、三超极化率,见《使用Multiwfn基于完全态求和(SOS)方法计算极化率和超极化率》(http://sobereva.com/232)。Multiwfn还能图形化展现(超)极化率,见《使用Multiwfn通过单位球面表示法图形化考察(超)极化率张量》(http://sobereva.com/547)。


3 使用Multiwfn解析polar关键词的输出

对于上一节讨论的各种情况,polar关键词输出的信息都能被Multiwfn所解析。注意route section一定要写上#P,否则输出文件将无法被正确解析。

启动Multiwfn后,先输入Gaussian的polar任务的输出文件的路径,比如C:\CH3NH2_polar.out,然后进入主功能24,选1,之后会看到下面的菜单:
-3 Set unit in the output, current: a.u.
-2 Set output destination, current: Output to screen
-1 Toggle loading frequency-dependent result for options 1 and 7, current: No
0 Return
1 "Polar" + analytic 3-order deriv. (HF,pure/hybrid DFT,semi-empirical)
2 "Polar" + analytic 2-order deriv. (MP2,double-hybrid DFT,TDDFT,CIS...)
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...)
7 "Polar=gamma" + analytic 3-order deriv. (HF,DFT,semi-empirical)

选项1~7对应上一节提到的不同情况,当前你在什么情况下计算的就选哪种,然后Multiwfn就会输出分析的结果。比如你用了CCSD结合polar=doublenumer就应当选5。之所以要对此进行选择是因为不同情况下Gaussian的输出格式相差极大,Multiwfn必须知道你在什么情况下计算的,才能相应地对输出内容进行解析。选项1~7都可以给出α,只有1、3、5才会给出β,只有7才会给出γ。

如果在HF/DFT/半经验计算时用了polar CPHF=rdfreq或polar=DCSHG做了含频(超)极化率计算,或者用的是polar=gamma,默认情况下只会解析静态的(超)极化率。如果想让Multiwfn分析含频的,应当先选-1将其状态切换为Yes。当Gaussian输入文件里定义了多个频率时,如0.1,0.22,0.35,用户可以选择分析哪个频率下的结果。对于β,用户也可以选择分析β(-ω;ω,0)还是β(-2ω;ω,ω),注意只有写了DCSHG时选择β(-2ω;ω,ω)才有意义。对于γ,可以选择分析γ(-2ω;ω,ω,0)还是γ(-ω;ω,0,0)。

默认情况下输出的单位是a.u.,如果选-3,也可以切换为SI或esu单位。默认情况下输出信息会显示到屏幕上,也可以选-2来更改为输出到当前目录下的polar.txt当中。

如果你使用的是DFT或MP2,且对当前体系又要做振动分析又要算极化率,那么没必要分别做freq和polar任务,只要做一次freq即可。用Multiwfn载入freq任务的输出文件,在Multiwfn此文介绍的界面里选择2 "Polar" + analytic 2-order deriv.即可提取极化率。


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.04620
XY=         -0.08830
YY=         28.89110
XZ=          0.00000
YZ=          0.00000
ZZ=         24.43580
Isotropic average polarizability:         26.12437
Isotropic average polarizability volume:       3.871231 Angstrom^3
Polarizability anisotropy (definition 1):          4.18643
Polarizability anisotropy (definition 2):          4.18363
Eigenvalues of polarizability tensor:      24.4358      25.0442      28.8931
Polarizability anisotropy (definition 3):          4.15314

如前所述,默认情况下所有输出的量都是以a.u.为单位。

极化率张量是个对称的3*3矩阵,故所以只有6个独立的矩阵元列了出来。每个体系的极化率的大小一般只用一个数值<α>来表示,也就是上面给出各向同性平均极化率,它是对XX、YY、ZZ取平均得到的。极化率也经常通过极化率体积来表示,上面也给出了。极化率并不是各向同性的,体系对不同方向射来的电场的响应是不同的。如上可见,Multiwfn也输出了极化率的各向异性值Δα,数值越大各向异性越强,球对称体系显然其值为0。各向异性可以以三种方式定义,分别如下所示,详见Multiwfn手册3.27.1节

ε代表极化率张量的本征值,它也被Multiwfn输出了出来。最常见的衡量极化率的各向异性程度是定义2,而定义3的形式一般是在衡量磁屏蔽张量各向异性时用的。

下面是第一超极化率部分的解析结果
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载入输出文件,并进入主功能24里的子功能1后,先选-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ω;ω,ω)情况的数据。


3.3 解析第二超极化率一例

这里我们对NH3计算第二超极化率γ。Gaussian输入文件是Multiwfn程序包里的examples\polar\NH3_gamma.gjf,内容如下

#p PBE1PBE/daug-cc-pVTZ polar=gamma
[空行]
pbe1pbe/aug-cc-pVTZ opted
[空行]
0 1
 N                  0.00000000    0.00000000    0.11367600
 H                  0.00000000    0.93902800   -0.26524400
 H                 -0.81322200   -0.46951400   -0.26524400
 H                  0.81322200   -0.46951400   -0.26524400
[空行]
532nm 680nm

想高精度计算小体系的γ,基组必须带非常充足的弥散函数,因此此例用了每个角动量加两层弥散的d-aug-cc-pVTZ,此基组非常昂贵。另外计算γ本身也很昂贵,所以哪怕是小体系,想高精度计算γ也需要给力的机子。polar=gamma关键词要求既计算静态的也计算动态的γ,如当前输入文件所示,对于动态的γ我们要计算的是532 nm和680 nm外场下的情况。

此例我们要解析γ(-2w;w,w,0)在680 nm下的结果。启动Multiwfn然后输入
examples\polar\NH3_gamma.out
24   //(超)极化率分析
1  //解析Gaussian的(超)极化率的计算输出
-1  //切换为载入动态(超)极化率(不选这个的话之后解析的将是静态的γ)
7  //对应polar=gamma的情况
2  //对应680 nm的情况
2  //解析γ(-2w;w,w,0)形式的动态第二超极化率

输出信息如下
 XXXX=     3091.740000
 YXXX=        0.000000
 ZXXX=        0.000040
 XYXX=XXYX=        0.000000
...略
 YZYZ=YYZZ=     3677.510000
 ZZYZ=ZYZZ=       -0.446631
 XZZZ=        0.000000
 YZZZ=       -0.400105
 ZZZZ=    18072.400000

 Magnitude of gamma:       5998.270291
 X component of gamma:     1574.323333
 Y component of gamma:     1574.284000
 Z component of gamma:     5569.774000
 Average of gamma (definition 1), gamma ||:     8718.3813
 Average of gamma (definition 2):     8737.028000
 gamma _|_:     2875.049333

上面每一项的定义笔者就不在这里介绍了,因为在Multiwfn手册3.27.1节里已经写得非常详细了,这些量都是考察γ时经常要涉及的。

如屏幕上的提示所示,上面给出的gamma信息对应的是输入朝向下的,而之前提到的Multiwfn给出的alpha和beta都是标准朝向下的,这个差异要注意。