| 查看: 5440 | 回复: 12 | ||||||||||
| 【奖励】 本帖被评价8次，作者sobereva增加金币 6.4 个 | ||||||||||
[资源]
使用Multiwfn分析Gaussian09的极化率、第一超极化率的输出
|
||||||||||
|
使用Multiwfn分析Gaussian09的极化率、第一超极化率的输出 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相差不算很大。 |
回复此楼» 收录本帖的淘帖专辑推荐
 量化软件学习 | SOB 集锦 | 理论化学计算知识 | Gaussian |
| 经验 | 计算-mulwiwfn | QuantChem | 自己感觉好的东西 |
» 猜你喜欢
全奖博士 英国利物浦大学 × 台湾清华大学 联合培养
已经有0人回复
-
寻博导
已经有13人回复
-
物理化学论文润色/翻译怎么收费?
已经有232人回复
-
基金正文30页指的是报告正文还是整个申请书
已经有9人回复
-
诚邀有志之士加入江苏大学环境学院合成生物学课题组！(长期接收学硕专硕、硕博连读)
已经有16人回复
-
一志愿南大理论与计算化学方向（070300）342分，求调剂
已经有1人回复
-
欢迎化学、化工、物理、材料、或计算机背景的考研咨询
已经有0人回复
-
哈尔滨工业大学（深圳）绿色化学团队招收2026年秋季入学博士生
已经有11人回复
-
澳大利亚昆士兰科技大学PhD博士奖学金
已经有0人回复
» 本主题相关商家推荐: (我也要在这里推广)
» 本主题相关价值贴推荐，对您同样有帮助:
- Multiwfn在AIM分析中的准确性如何？
已经有7人回复
- 利用Multiwfn计算Gaussian输出的激发态之间的跃迁偶极矩
已经有5人回复
- Multiwfn结合VMD绘制AIM拓扑分析图
已经有12人回复
- 求助高斯输出文件中坐标系怎么定义的？
已经有11人回复
- 使用Multiwfn结合VMD分析和绘制分子表面静电势分布
已经有14人回复
- 同一个结构2*1*1和1*1*1两个模型计算的静介电常数不同，介电函数曲线也不同
已经有4人回复
- 在gaussian计算中，先用小基组计算，在改成大基组计算出现的问题
已经有23人回复
- 关于高斯计算分子表面积等参数的问题
已经有10人回复
- 建议multiwfn开发GUI界面 文件打开什么的效率太低了
已经有10人回复
- 急急急！高斯生成不了wfn文件怎么办
已经有7人回复
- 使用Multiwfn绘制原子轨道图形、研究原子壳层结构及相对论效应的影响
已经有17人回复
- 使用Multiwfn绘制NBO及相关轨道
已经有30人回复
- 引用了multiwfn的论文。现在multiwfn软件已经有19次引用率了呢
已经有13人回复
- 使用Multiwfn做拓扑分析以及计算孤对电子角度
已经有18人回复
- 回答一些关于Multiwfn的疑问以及未来Multiwfn的发展打算
已经有37人回复
- 【求助】6-31G*中*的含义
已经有3人回复
- 【求助】如何读取输出文件中的键级
已经有4人回复
- 【求助】如何在以下高斯03输出文件中得到二阶超极化率
已经有15人回复
2楼2014-04-27 08:51:22
3楼2014-04-27 12:09:38
5楼2014-04-28 05:29:50
6楼2014-04-28 09:54:59
7楼2014-05-23 16:40:02
8楼2014-09-24 15:39:22
9楼2016-01-13 22:34:06
10楼2017-05-04 15:30:55
11楼2017-05-04 15:32:05
12楼2017-05-04 15:31:44
|
本帖内容被屏蔽 |
13楼2021-12-13 19:03:18
简单回复
安德4楼
2014-04-28 00:03
回复
五星好评 
50