★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ loovfnd(金币+3 ,VIP+0 ):SCF不收敛的标志什么的说清楚了最好了,呵呵,中秋快乐! 10-3 15:04 wuli8(金币+10 ,VIP+0 ):辛苦了。 10-20 19:37 yjcmwgk: 回帖置顶 2012-07-14 09:44:48
对GaussianFAQ的改进建议,请在http://muchong.com/bbs/viewthread.php?tid=1734563 回帖。
第二章:Gaussian使用:
(2.1)
版本一:高斯能不能做XXXXXXX方面的计算啊?
版本二:高斯的主要功能是什么?引用回帖: 分子能量和结构
过渡态的能量和结构
化学键以及反应的能量
分子轨道
偶极矩和多极矩
原子电荷和电势
振动频率
红外和拉曼光谱
核磁
极化率和超极化率
热力学性质
反应途径
(2.2)
请问2070是什么错误啊?
后续问题:什么是log/out文件的错误提示啊,在哪儿看啊?引用回帖: 2070错误
最基本的解释就是,你的系统在安装了g03/g09之后,运行时对内存调用出错.....
解决措施:
1 更换稳定版本的系统, 采用完整的正版系统.
2 使用linux平台,采用稳定版(非最新版), 如果linux下遇到Segementation Fault,在找别的方法处置, 一般只需一个指令,就可以一劳永逸....
对于一个新手,你不必管什么是2070错误。是什么错误应该看最后终止在哪一步,例如L1就是命令行里有错误,L301通常是自选多重度错误等,而且在结果文件中最后部分会给出错误的原因(就是error termination这一行【图7】),以便联查修改!
(2.3)
scf不收敛怎么办?几何优化不收敛怎么办?引用回帖: 虽然二者有联系,但实际上是两个相对独立的问题,要分开回答。
这里提供的只是解决方案,所有的方法都用上了还是不收敛得情况也是有的。
首先要分清scf不收敛和几何构型优化不收敛:scf不收敛指的是自洽场叠代不收敛(什么?没听说过什么叫自洽场?那还是回去学习些量化基础知识再开展计算吧!),可以认为是对指定结构的波函数不断优化的过程,是为了找到这个某个指定结构下能量最低的波函数,而几何构型优化是对结构的优化的过程,是为了找到某个指定的组分下能量极小结构(注意,不一定是能量最小结构)。在量子化学计算的几何构型优化中,每一步的几何构型优化都包含的很多次的scf计算。
1、scf不收敛的解决方案。
(1) 可以加大scf的循环次数,默认的循环次数是128次,通过scf=(maxcycle=n)来设置最大循环次数n。建议不要超过512,更多的循换没有必要。
(2) 如果加大循环次数不管用,在分子有对称性的情况下,使用scf=dsymm关键词来强制密度对称,有时可以收敛。另外,此关键词很多时候对"scf is confused”这种错误很管用。
(3) 使用scf=symm关键词,使用的前提同上,有时可以收敛。
(4) 如果(2)(3)两步都不行,可以将对称的分子中的某几个原子的位置微调,使分子丧失对称性。这等效于nosymm关键词,但个人经验,这种方式比nosymm好用的多。
(5) 如果还不行,只能拿出杀手锏了,就是使用qc,但不建议直接使用,而是使用xqc关键词,比如scf=(maxcycle=80,xqc),意思是如果scf正常计算(dc)在80个循环之内不收敛才进行昂贵的qc计算,因为scf不收敛多数在几个优化的过程中出现,无法判断哪一步优化的时候会出现scf不收敛,所以用xqc比纯粹使用qc要省时的多。
(6) 中级用户可以在输入文件的井号“#”开头那一行井号后面加上字母"p"来输出更多的信息,其中就有自洽场叠代的信息,分析原因可能会对采用什么方法提供指导。
(7) 前面有虫子提到一个有用的方案但没说清楚,我这里补充一下:如果用用小基组计算,scf可以收敛,那么保存好检查点文件,换成大基组的时候从检查点文件中读取初始猜测(使用guess=read关键词),有时可以算过去。
(8) 如果你计算中使用了含有弥散函数的基组(比如6-31+G*基组中的“+”,Aug-cc-pVTZ基组中的Aug等),可以先去掉弥散函数(比如例子中的两个基组变为6-31G*和cc-pVTZ),这样很容易收敛,需要的话,换成带有弥散函数的基组时,从检查点文件的读初始猜测。原理与(7)类似,有时候管用。
注意:(7)和(8)不适合很大的体系,此时,使用小基组或者不带弥散函数的基组与后来使用的大基组在基组数量上差得太多(如果你本来就遇到了scf收敛问题,而使用这两个方案时,机组数量差距超过512个,多数时候要出问题),反倒容易引起scf不收敛,因为读检查点文件时,对那些多出来的基组的猜测是空的,而直接计算时最起码还有个半经验的计算提供了相对合理的初始猜测。
(9) 上述方法有时可以组合使用,如果经过各种组合处理都不收敛,那么放弃吧,你的分子的电子结构太差了。
2、几何构型优化不收敛的解决方案。引起这个问题的原因比较多,可能也说不全,尽量做到全面。
(1) 如果是很小的分子(10原子以内),初始结构可能离平衡结构比较远,又在输出文件优化的最后一步判断是否收敛的位置(可以通过查找Threshold字段找到,它下面有四个判断项,都是YES才代表优化收敛)看到了“-- Number of steps exceeded, NStep=xxx”,可以通过加大优化的循环次数来解决问题。使用关键词opt=(maxcycle=n).
(2) 优化到如果四个收敛标准中前两项早就收敛了,而后两项尤其是第三项不收敛,这需要判断原因。如果每一步优化过程中的能量始终在下降,那么可以继续让他算,超过了最大步数停掉了的话把结构拿出来,重新提交就行了。如果能量忽大忽小出现跳跃,说明遇到了比较平坦的势能面,或者优化的算法不好。此时建议按如下顺序处理:
a. 使用opt=maxstep=n来缩小最大步长,即原子移动的距离,n的默认值是30,只能设置整数。
b. 使用关键词opt=gdiis,很多时候可以很快收敛。
c. 使用opt=calcfc关键词,在几何优化的第一步计算力常数,为优化定一个大方向,类似在第一个三岔路口选一条路。
d. 使用opt=calcall关键词,在几何优化计算的每一步都计算力常数,类似于在每一个三岔路口都来选一条路。这个关键词不到万不得已不要使用,非常耗时。
(3) 四个判断标准的后两项早就收敛了,但前两项似乎固定在某个数值上了。遇到这种情况,首先去检查你优化第一步的时候给出的每个轨道得本征值部分,也就是轨道的能量,如果HOMO和LUMO的能量相同或者非常接近,那么你的设置有问题,即自旋多重度不对。如果HOMO-LUMO能量没问题,那么将目前的计算停掉,提出结构,使用opt=gdiis关键词,很多时候管用。
(4) 如果是在几何构型优化的过程中出现因为scf不收敛导致无法优化,可以先用上面scf不收敛的解决方案来处理,如果都不行,可以通过在#后面添加p输出更多信息,如果看到每次scf能量都能收敛到10的-6次方以下然后开始出现能量的跳跃,就是收敛不到10的-8次方(这种情况还比较常见),可以使用scf=(conver=6)来让几何优化继续下去,等优化好了,再将此关键词去掉,我的经验,优化的结构一般不会有太大的差异。如果不到10的-6能量就开始跳越,不适合此方法。
(5) 如果几何优化过程中发现键变得异常短,而且此时几何优化总是不收敛。遇到这种情况,问题出在你的基组使用,肯定是使用了赝势基组却没有读入赝势。由于赝势基组将内层电子和原子核用一个有效核心势来表示,你没有读入赝势就说明没有读入这个有效核心势,那么相当于把原子的内层电子都给除掉了,巨大的核引力把原子之间的距离变得异常短也不足为奇。
比较通俗的解答请看2.10
(2.4)
Gaussian并行效率如何?引用回帖: 暂时的最佳答案
这个一般要看你的节点数 和 CPU 以及内存
集群测试环境:
操作系统:CentOS4
1个管理节点:2×CPU5410 16G内存 5×300GSAS硬盘做RAID5
8个计算节点:2×CPU5410 16G内存 2×143GSAS硬盘RAID0
以太网络链接:千兆网,TCP/IP协议
G03D01并行版(Linda7.1)
任务说明:主要测试并行SCF,共882个基函数,无对称性。
======================
CPU数 运行时间 效能
n xn e
--------------------
1 1920
2 1314 1.168
4 803 1.157
6 611 1.334
8 527 1.648
10 479 1.992
12 455 2.416
14 434 2.797
16 432 3.355
======================
效能e是我自己定义的:e=(n*xn-1*1920)/(1920-xn)
其中分子表示用n个CPU并行处理时比单个cpu多付出的cpu时间(计算资源而非实际时间),分母表示用n个CPU并行时节省的实际时间。
效能e体现的是以单个CPU运行的时间为基础,用n个cpu并行运行任务时,每节省1秒钟实际时间所需要付出的CPU时间代价即计算资源。其值越小越好。
e即考虑了实际时间的节省又包括了计算资源的消耗。
从上表可以看出,当使用两个计算节点4个CPU时,e最小,并行效率最好。因此推荐使用4个CPU并行计算单个高斯任务,兼顾速度和效益。
(2.5)
什么是自旋多重度?怎么设置?
后续问题:怎么判断单电子数啊?还有那个2|S|+1是怎么回事儿啊?S是什么东西?怎么算S啊?引用回帖: 量子化学中用到的自旋多重度是用来区分一组简并的波函数,这些波函数之间只存在自旋角动量的不同。【图8】
自旋多重度定义为2S+1,其中S是自旋角动量,它与体系内的单电子数(N)相关。S=N/2,所以说到底最简单的判断方法:自旋多重度等于单电子数+1。
怎们判断体系单电子数需要计算者对自己研究的体系有着深入的了解,这是开展科研的前提,不是偷懒的理由。
但有一点可以提醒Gaussian初学者,计算体系总电子数时需要加上体系的电荷,这样才不会算错单电子数。
如果实在无法判断到底有几个单电子,那么根据体系的总电子数,算各种可能情况,找那个能量最低的。
体系总电子数为偶数:自旋多重度可能是1,3,5,7....
体系总电子数为奇数:自旋多重度可能是2,4,6,8....
===============
现在介绍自旋多重度=2S+1的由来,没兴趣的略过,有兴趣可以参考。
大学的化学书中提到了决定原子中每个电子唯一性的n, l, m, 和m(s)四个量子数,其中这个m(s)就是自旋量子数,理论研究表明它只有两个值,即:+1/2和-1/2,由于现在绝大多数量化计算都是基于原子轨道线性组合分子轨道的方法,所以此套参数可以扩展到分子体系,只不过换了个说法而已。
如果一个原子轨道或者分子轨道的两个自旋轨道是充满的,因为其他三个量子数都一样,那么两个自旋量子数可以相互抵消,也就是这个轨道对体系的自旋角动量不产生贡献。
如果体系中所有的占据轨道的两个自旋轨道都是充满的,那么整个体系是没有自旋角动量,也就是只有一种自旋状态(没有自旋)。如果体系内有一个占据轨道之占据了一个自旋轨道,而另一个自旋轨道是空的,那么体系的就分成了两个部分,一部分是没有自旋的全部充满的轨道构成,另一部分就是这个单电子产生的+1/2或者-1/2的自旋,所以此时体系的自旋多重度是2,如果体系有两个以上的单电子,根据电子排布规则,这些成单的电子优先占据与其他单电子占据的自旋方向一致的自旋轨道,所以从这个角度来说,这些电子之间的自旋是不能抵消的,而只能累加。这样就导致体系分成了单电子数+1个部分,即全充满的所有轨道是一个部分,而单占据的每个自旋轨道又是一个独自的部分,体系的自旋多重度就是单电子数+1。
由于每个单电子产生|1/2|的自旋,那么根据自旋多重度=单电子数+1,自旋多重度也就等于总自旋2*|N*1/2|+1,而绝对值符号中的部分正是体系总自旋角动量的表达式,即S,这也就是自旋多重度=2S+1的由来。
(2.6)
什么是基组啊?什么是混合基组啊?怎么个写法啊?
后续问题:不好意思,楼上说的我还是不太明白。什么是你们所谓的Gauss基函啊?基组大小你们又是怎么知道的呢?
后续问题:请问6-311++g(3d2f,2pd)中各项的含义是什么啊?引用回帖: 基组是体系轨道的数学描述,对应着体系的波函。将其带入到薛定谔方程中,就可解出体系的本征值(也就是能量)。混合基组就是对不同的原子使用不同的基组!写法如下:
在命令行加入gen Pseudo=Read
========
C H 0
6-311++G**
****
Pt 0
LanL2DZ
****
Pt 0
LanL2DZ
========
Gaussian中常用的分子轨道是名为高斯轨道的高斯函数的线性组合。基组是体系内轨道的数学描述。大的基组由于对电子在空间上有小的限制而具有更大的精确性。
基组的大小:最小基组包含了描述轨道的最少的函数数量。例如STO-3G 是最小基组,每一个基本函数中含有三个高斯函数。增大基组的第一个方法就是增加每个原子基函数的数量。分裂基组,比如3-21G 和6-31G,对于价键轨道都用两个函数来进行描述。我个人认为大的基组就是加入限制,对体系的描述多。
6-311++g(3d2f,2pd)中表示:增加了更多的极化函数,包括三组分裂的价键基组,在重原子和氢原子上加弥散函数,在重原子上加的三组d 函数和两组f 函数,在氢原子上加的两组p 函数和一组d 函数。
如果你不知道你所选择的基组群的计算规模,你可以加入一个命令%kjob L??? 1,合适的Link下,添加一个kjob命令,来测试你的机组规模。
解释:%KJob LN [M]指示程序在执行Link N的第M个事件之后停止。例如%KJob L502 2表示执行模块502两次后停止。如果M值省略,默认值为1。
对于第二周期元素,几个常用基组的基组规模测试结果
basis functions primitive gaussians
sto-3g 5 15
3-21g 9 25
6-31g* 15 28
6-31+g** 19 32
6-311++g** 22 36
6-311++g(3df) 39 58
6-311++g(3d2f,2pd) 46 68
sddall 8 16
tzvp 19 35
aug-cc-pvdz 23 43
aug-cc-pvtz 46 73
aug-cc-pvqz 80 127
aug-cc-pv5z 127 209
aug-cc-pv6z 189 326
具体测试方法见
http://muchong.com/bbs/viewthread.php?tid=1498413&fpage=1 (2.7)
怎么选择泛函,怎么选择基组?引用回帖: 这不是一句话能说清楚的。
关于这个问题,有篇文章讲得特别好
JPCA 2007, 111, 10439-10452
-----------
以下均为个人经验,搞错了不负责
《什么样的泛函和基组的组合,适合做的事情》
这些都是我自己科研经验和以前读的文献的总结。也没写什么参考文献,就是随便扯扯,这就全凭自己的经验和记忆啦,想到哪儿写到哪儿。这里面肯定有不对的,甚至有荒谬之极的东西,还请大家指正。也希望自己抛砖引玉,其他高手也讲讲自己的经验。
很少见到高手们总结自己的科研经验和阅览的文献。小卒我这个屁股上还没褪干净蛋壳的小菜鸟,只好抛出一块又丑又粗糙的破砖,引引高手们的美玉啦。
以下探讨均不涉及cluster。主要是有机化合物、以及金属离子和有机配体组建成的配合物。其计算也都是基于单个分子的计算。周期体系不在考虑之列——这是一句无奈的词句,因为我的研究方向不在此处,所以我确实没有仔细思考过它们
1分子结构、键长和键角:
对于比较轻的元素,比如CHON之类,b3lyp/6-311++g**就很优秀,如果把基组加到了aug-cc-pv?z的程度,就没什么大的必要了。有人甚至说,Cl以前的元素用6-31g就好,如果有孤对电子加个d就行了,如果带负电荷就加个弥散。
对于一些中等重量的元素,LanL2DZ或者DZVP之类的基组就不错啦,没必要加更高的机组了,加了也是浪费。Fe的LanL2DZ有人换成SDD试过,换了以后得到不同基态结果 ,千万小心使用SDD。第一行过渡金属,如果只有一个,那么6-311+g*的计算量还是可以承受的;两个及两个以上,非赝势很难算。lanl2dz和sdd都是比较好的选择。除非做单点算,否则非要上mp2算就没太大必要了。
再重一些,比如稀土,LanL2DZ就不行了,M. Dolg先生搞出来的那个基组,斯图加特RECP的[5s4p3d]-GTO占了大部分市场份额。而泛函选择上,百分之四十用的是B3LYP,百分之四十用的是B3PW91,百分之五用的是MPn族,百分之五用的是CCSD(T)族,百分之五用的是PBE族,百分之五用的是经过改进的半经验。有时候还要自己选择,是用小核赝势呢,还是大核赝势呢。听说ADF用PBE/DZPZORA处理稀土,效果也不错,但是我没用过。
有金属的时候,据说tpss不错,但是我也没用过,不敢说。
这里我要批一批我偶然见到的一个文献:Int.J.Electrochem.Sci. 2009, 4, 295-307,伊朗的一群人做的工作。他们居然用6-31G*计算稀土元素,这简直是扯淡,大家不要学他。
半经验方法中,处理CHNO之类的元素,AM1和PM3都不错,PM6也很好,他们可以作为我们前期粗略优化的手段。不过要注意AM1处理共轭体系不太好。有人用Sparkle/AM1和Sparkle/PM3处理稀土元素。精度并不是多么好,但是极大地减轻了计算量,可以说还是很有成效的一个工作。不过这里不探讨半经验,仅仅探讨密度泛函,所以不多说了。
一个例子:[Fe(CN)6]3-,我会选择B3lyp/genecp来优化其分子结构,然后用PBE1PBE/genecp来计算其电子结构。genecp对应关系如下(其实,这个混合基组有自己的专门的名字,叫做LACVP+*)
Fe 0
LanL2DZ
****
C N 0
6-31+g*
****
Fe 0
LanL2DZ
分子结构计算结果的优劣可以用有2种方法加以检验,一是与可获得的实验数据进行对照(但不要忘了实验数据本身只是接近事实,而不是准确事实的本身,因为它们仍然是有误差的),二是在优化结构的基础上计算性质,同样的,还是要拿计算结果与更为有限的性质实验测定结果(它们同样有误差)进行对照。这充分体现了人类认识自然规律的有限性——我们总是只能获得相对真理。由于现有的计算理论本身在处理相关能和相对论效应时有很多困难(理论的和实际计算上的),所以,本人的研究都不涉及重金属,而只做轻元素体系,多为生物分子体系,大量的计算实验表明DZP加弥散基组与B3LYP的结合所得的结果可靠,所以一直使用。个人的精力是有限的,无法大量进行方法的尝试与鉴别,所以,只有大量阅读文献来弥补。实际上,就轻元素物质的结构而言,6-31++G**的结果已然非常不错了,但在计算性质时不妨对基组适当选择。如前所述,因为实验有误差,所以有时实验并非一定可靠,更重要的是,有的性质的测定,实验上是不可能进行的(至少目前是如此),因而,计算怎样走在实验的前面,恐怕值得思考。
RI-DFT、RI-MP2方法,能比常规方法快一个数量级,但精度接近。此外还有RI-CCSD、RI-CCSD(T)等,不过用得很少。Gaussian 0x的话可以试试纯GGA+密度拟合,不会比ADF的同类方法差太多。
2氢键
氢键的计算,还是以b3lyp、b97-?、pwb6k、b2plyp、b3p86、pbe1pbe、mpw1k这几族泛函比较好。基组的话,还是aug-cc-pv?z、mg3s好一些。其实劈裂价键基组也不错,但是必须加“大弥散”和“高角动量”,因为氢键还是比较弱一些,氢键对应的波函数必须加大离域程度。
其实氢键还算是一个比较好算的东西。流行的DFT中,除了OLYP和O3LYP之外,基本都有一定精度。氢键的计算,基组的选择影响大一些,大弥散基组是必须的。
虽然dft在研究氢键方面也一直有人在用,但现在大多数的审稿人总是对这种方法持怀疑的态度,因为对扩散相互作用,dft无法预测,而审稿人更为信任的是mp2,而且一般需要大基组来计算,如aug-cc-pvdz,aug-cc-pvtz等,但如果无法完成这些基组,一些低一些的基组来做计算也是可以被接受的,如6-31G++(d,p),6-311G++(d,p)等,但bsse校正是做氢键方面一个必须要做的一部分。如aug-cc-pvdz,aug-cc-pvtz等这些大基组,bsse似乎影响不大,可以不做。
3远程弱作用
远程作用就比较扯淡了。B3交换泛函彻底完蛋。虽然我主要也是B3LYP阵营的人,但是要我在这个领域夸夸B3LYP,我还真不能瞪着眼睛说瞎话。这个领域DFTB和DFTB+曾经试图在这个领域做些贡献,但是成就并不是很大。比较成功的是DFT-D,包括b97-d,b2plyp-d等等这些加上dispersion项的DFT了。
Truhlar的泛函 M05, M05-2X, M06, M06-2X等等可能只适用于某些体系,对于有些体系,照样给出糟糕的结果。LDA居然能出乎意料地给出vdW作用,而其他GGA泛函都不行。这应该是个偶然,可能是正负误差相消的结果。比如,用LDA优化出的石墨层间距和实验值符合得相当好。PBE-D方法总是倾向于高估vdW作用。
远程作用还有别的选择——尽管它会大大增加计算量
一种选择,使用MPn,还是不够准
另一种选择,是把MPn改造成SCF-MPn,计算量忒大
再一种选择,使用CCSD(t)之类的,……呃,还是算了吧,计算量我更受不了。
文献常见的高斯09的方法有:MPW1KCIS,MPW1B95,MPWB1K,jiewei先生推荐基组tzv(2d,2p)/QZV3P。MS的话一般用GGA/PEB。
基组问题:算远程作用必须用极大的基组——这是一句废话,因为你要让两个分子相互作用,那么它们的波函必须能“够到对方”才行。从TZV(d),到TZV(d,p),到TZV(2d,2p),到TZV(2df,2pd),到TZV(3d2f,3p2d),……继续加吧,能加多高,就加多高,算远程弱作用,您永远别吝啬你的资源。当然,我绝不反对你使用DZV族
另一个选择:从6-311++g(d,p),到6-311++g(2d,p),到6-311++g(2df,2pd),到6-311++g(3df,3pd),到6-311++g(3d2f,2p2d),……继续加呀加呀加呀加!
再来个选择:从aug-cc-pvtz,到aug-cc-pvdz,到aug-cc-pvqz,……子子孙孙无穷尽也
听说过一族巨型机组,叫ANO基组,当然还有它的一些变体。ANO规模很巨大,算的能量相当精确,但是对于分子性质没有那么好的效果。
4多加一句:关于核磁,给某些人提个醒儿
昨天,有人在qq上问我,他说核磁计算结果与实验值相差太大。我问他你用啥算的?他说是b3lyp/6-31g*。我直接就无语+仆街了。
很多人习惯用6-31g*。这确实是一个很优秀的基组。但是您居然扛着6-31g*跑去算核磁,误差不达到百分之一千才怪。算核磁最小的基组,至少要用6-31++g(d,p)。这已经是最低要求了,不能再低了,建议使用6-311++g(2df,2pd)或更高
(2.8)
NBO是什么?怎么用啊?引用回帖: (2.9)
消除虚频的方法和经验?引用回帖: 首先,搞懂什么是频率
中学的时候我们学过简谐振动,对应的回复力是f=-kx,对应的能量曲线,是一个开口向上的二次函数E=(k/2)x^2. 这样的振动,对应的x=0的点是能量极小值点(简单情况下也就是最小值点)。这时的振动频率我们也会求ω=(2π)sqrt(k/m)。显然它是一个正的频率,也就是通常意义下的振动频率。
然后,从频率的概念引出虚频的概念。
一维情况下,如果能量曲线是一个开口向下的二次曲线呢?首先,从能量上看,这是个不稳定的点,中学的物理书上称为“不稳平衡”。用现在的观点看,就是这个点的一阶导数是零(受力为0),且是能量极大值。如果套用上面的公式,“回复力”f=-k'x(实际上已经不是回复,而是让x越来越远了),这里k'是个负数,ω=2π sqrt(k'/m)显然就是一个虚数了,即所谓的虚频。Gaussian里面给出一个负的频率,就是对应这个虚频的。
实际情况下,分子的能量是一个高维的势能面,构型优化的时候,有时得到了极小值点,这样这个点的任意方向上,都可以近似为开口向上的二次函数,这样这里对应的振动频率就都是正的。对于极大值点,在每个方向都是开口向下的二次函数,那么频率就会都是负的——当然一般优化很少会遇到这样的情况。对于频率有正有负的情况,说明找到的点在某些方向上是极大值,有些方向上是极小值。如果要得到稳定的能量最低构型,显然需要通过微调分子的构型,消去所有的虚频。如何微调?要看虚频的振动方向。想象着虚频对应的就是开口向下的二次函数,显然,把分子坐标按照振动的方向移动一点点,分子应该就可以顺着势能面找到新的稳定点,但是也不能太小。而所谓的过渡态,则是连接反应物和产物之间的最低能量路径上的能量极大值。好比山谷中的A,B两点,它们之间的一个小土丘,就是过渡态,从A到B的反应,需要越过的是这个小土丘,而不是两边的高山。这样,过渡态就是在一个方向上是极大值,而在其它方向上都是极小值的点。因此,过渡态只有一个虚频。
在进行Gaussian计算时,经常会出现多个虚频个,每个都只有-30~-20左右,但是如果结构改动过大,则一不小心,虚频就会更大。消除虚频的主要方法是改变构型。其次在计算上还可以尝试:nosymm;加大循环次数;提高收敛度;iop(1/8=?)等。通常对于势能面很平,较小的虚频很难消除。为了解决以上的这种现象,按如下的步骤一步步调整:
1). 增加收敛标准
在优化时采用Scf(tight)的选项,增加收敛的标准。再去计算频率。如果还有虚频,参见下一步
2). 降低对称性
对称性的影响,很多情况下的虚频是由于分子本身的对称性造成的。这样,在优化时,如果必要,要将对称性降低,还有,输入文件有时是用内坐标。建议如果有虚频的话,将内坐标改成直角坐标优化。但是有时这一方法也并没有多大作用。
3). 坐标中和法
如果上述方法还有虚频,看一下虚频,找到强度较大的,将在频率中产生的原子的振动坐标加到相应的输入文件中。这样,重新计算,直到虚频没有。比较好的做法是分子结构输入用直角坐标,将这些产生虚频的振动坐标直接与原始坐标相加。加的时候不一定要将坐标直接相加,可以选择权重f=0.1-1.0相加等。
4). 其他方法
实际上,如果分子柔性较大,很难找到最低点,这是电子结构计算的问题,这种情况下,需要动力学的东西,用构象搜寻的办法解决。如:模拟退火,最陡下降法,淬火法等。将得到的能量最低的构象做一般的电子结构计算,这样,应当没有问题。
如果上面几步仍然无法解决问题,则有可能是因为计算还没有得到最稳定的结构,那么或者是分子结构有问题,或者是计算错了,再者就是游离出现在计算的范畴。
最后如果是寻找过渡态,如果只有一个虚频,而且该虚频对应某单键旋转震动,那么产生虚频的原因就是初始结构比较接近重叠构象。这时把优化过的结构绕振动方向旋转几十度再重新优化就有可能得到正确的结果。优化中等大小的有机分子(50-100原子)的时候,应该经常碰到这种情况。 当出现两个以上虚频时,可以先看一下振动模式。如果多余的虚频,和想要的振动模式发生在不同原子上,可以先固定过渡态模式所对应的原子,然后做opt到minimal,一般可以把多余的虚频去掉,再放开做全优化。或许会有所收获。
(2.10)
L9999错误是什么?如何解决?
L502错误是什么?如何解决?
这简直是初学者遇到的最频繁的两种错误了引用回帖: L9999错误是属于构型优化没有找到最后结果,即在规定的优化步骤内没有找到最终的合适的结构。(L9999错误)。
对于优化不收敛,即L9999错误,实际上是在规定的步数内没有完成优化,即还没有找到极小值点。(或者对于过渡态优化,还没有找到过渡态),可以增加优化步骤,或是拿最后的结构 接着进行优化。
此外这有几种可能性:
1. 看一下能量的收敛的情况,可能正在单调减小,眼看有收敛的趋势,这样的情况下,只要加大循环的步数opt(maxcycle=200),可能就可以解决问题了。
2. 加大循环步数还不能解决的(循环步数有人说超过200再不收敛,再加也不会有用了,这虽然不一定绝对正确,但200步应该也差不多了),有两种可能。一是查看能量,发现能量在振荡了,且变化已经很小了,这时可能重新算一下,或者构型稍微变一下,继续优化,就可以得到收敛的结果(当然也有麻烦的,看运气和经验了);二是构型变化太大,和你预计的差别过大,这很可能是你的初始构型太差了,优化不知道到哪里去了,这时最好检查一下初始构型,再从头优化。
3. 对于L9999快达到收敛时,考虑减小优化步长有时对于能量振荡的情况也是有用的,opt(maxstep=1).
一个建议是,对于大体系,难收敛体系,先用小基组,低精度算法优化一下,以得到较好的初始构型,再用高精度的计算接着算。如果前面的方法保留了chk文件,重新计算时需要使用 geom=allcheck 读入构型(就不必麻烦地写构型了), guess=read(读入初始波函数,可以加快第一步SCF收敛)。
L502错误对应的是SCF不收敛,通常有以下的解决方法:
a. 使用小基组,或低级算法计算,得到scf收敛的波函数,用guess=read读初始波函数。
b. 使用scf=qc,这个计算会慢,而且需要用stable关键字来测试结果是否波函数稳定。如果这个还不收敛,会提示L508错误。
c. 改变键长,一般是缩小一点,有时会有用。
d. 计算相同体系的其他电子态,比如相应的阴离子、阳离子体系或单重态体系,得到的收敛波函数作为初始猜测进行计算。
更专业的回答请参见(2.3)
(2.11)
writwa错误是什么?如何解决?引用回帖: 要么是内存问题。可以增加内存(%mem=?)。
要么或是硬盘问题,两种情况:
一种情况是,如果你使用的电脑是32位的操作系统,或者你使用的Gaussian软件是32位软件,两者满足其一,那么你的计算最多只能写16GB的硬盘空间,即使你的电脑剩余空间再多,需要的空间超过了16G也要出错,如果是这样的,在命令行尝试使用maxdisk=2000mw关键词,也许会解决问题,如果不行,你只能找64位的版本来做计算了。windows下使用老版本的gaussian作计算时还会出现单个文件2G的限制,这时需要在将这16G拆分成8个左右的磁盘读写文件。在命令行的上面添加
%rwf=name1.rwf,2gb,name2.rwf,2gb,name3.rwf,2gb,.......这里的文件名随便起
另一种情况是你的系统和软件都是64位的,这种情况下如果你没有定义最大硬盘,gaussian就会无限的写下去,直到满足要求为止,这时经常会出现硬盘被写满的情况,也就是你的硬盘没有空间了,此时同样需要定义maxdisk关键词,容量至少要比你的硬盘剩余空间小一些。
(2.12)优化震荡怎么解决?引用回帖: 1.尝试改变结构。首先略微减小键长,接下来略微增加键长,接下来再对结构作一点改变。
2.关闭DIIS外推(SCF=NoDIIS)或改进的外推(GDIIS)。同时进行更多的迭代( SCF=(MaxCycle=N) )。
3.使用强制的收敛方法。SCF=QC通常最佳,但在极少数情况下SCF=DM 更快。不要忘记给计算额外增加一千个左右的迭代。应当测试这个方法获得的波函,保证它最小,并且正好不是稳定点(使用Stable关键字)。
4.将两个迭代震荡构型相加,重新进行优化:
两个震荡构型
分别是
C x1 y1 z1
N x2 y2 z2
O x3 y3 z3
和
C x1' y1' z1'
N x2' y2' z2'
O x3' y3' z3'
那么你要用的新坐标就是
C (x1+x1')/2 (y1+y1')/2 (z1+z1')/2
N (x2+x2')/2 (y2+y2')/2 (z2+z2')/2
O (x3+x3')/2 (y3+y3')/2 (z3+z3')/2
5.减小积分步长,使用iop(1/8=?)单位是0.01埃。默认是30,你最小可以调整到1,但是要注意,这样很容易严重增大opt循环步数。
(2.13)chk文件转化成fchk文件时默认内存不足,或者从fchk文件生成cube文件时默认内存不足,怎么办?引用回帖: 解决方法:
(1)在winxp下,依次点击控制面板,系统,高级,环境变量,在上下两个输入框中分别加入一个新的环境变量——变量名GAUSS_MEMDEF,变量值200000000,然后重启电脑
(2)在Linux下,输入命令export GAUSS_MEMDEF=200000000,然后回车
这样就把formchk和cubegen可以使用的内存,从默认值6M增加到了200M
如果你想改到400M,GAUSS_MEMDEF这个环境变量你就改成400000000,以此类推。这个值只认比特,不认K\M\G这些符号
(2.14)溶剂化计算时如果溶剂不是高斯默认的话怎么处理?引用回帖: eps=? 溶剂的介电常数。
rsolv=? 溶剂半径,单位是埃。
density=? 溶剂的密度。
具体看高斯说明书中关于scrf关键字的解释
小卒版主说,这个问题的答案,以后有时间会继续扩充。先做个标记
(2.15)对一个有很高对称性的分子,计算中出现"Inaccurate quadrature in CalDSu."错误应该怎么处理?引用回帖: 这个错误并不是出现在优化的过程中,而是出现在第一次scf迭代时,不是scf不收敛,而是scf出了问题。往往没有什么好办法,只能稍微调整构型来降低了一些对称性(比如Ih转为D5d)才算过去。这个意义上讲,symmetry(pg=?)应该也可行。
guess=indo是一个可行的解决办法,但是只适用于很小的体系。大体系用这个关键词会出现不可克服的scf不收敛,因此他的问题可能出在体系太大,只能用降低对称性来处理。
小卒版主说,这个问题的答案,以后有时间会继续扩充。先做个标记
(2.16)什么是BSSE校正?为什么要做BSSE校正?怎么做BSSE校正?引用回帖: 计算A、B分子间的弱相互作用能,不能简单地通过E(interaction) = E(AB) - E(A) - E(B)来计算,因为E(AB)能量相对于E(A) + E(B)的降低来自两方面,一方面是真实的A、B分子间的相互作用能,这是我们要求的;另一方面来自于A、B分子的基函数在复合物体系中重叠,相当于增大了复合物的基组而使E(AB)能量降低,这个部分贡献如果也掺入E(interaction),则高估了相互作用能,所以要去掉,它称为Basis Set Superposition Error(BSSE)。所以双分子的相互作用能应该表述为E(interaction) = E(AB) - E(A) - E(B) + E(BSSE)。有时在分子内相互作用能计算时也要考虑BSSE。
计算E(BSSE)有多种方法,Gaussian03用的是目前广泛使用的Boys、Bernardi发展的counterpoise correction方法。设E(i)为第i个分子在自身基组下的能量,E(i)'为第i个分子在全部n个分子上的基函数都出现下的能量,则计算n个分子相互作用能中的E(BSSE) = ∑( E(i) - E(i)' ),这里E(BSSE)为正值。注意计算E(i)与E(i)'时的分子几何结构与处在复合物中时的一致。
要计算A、B两个分子的相互作用能,在Gaussian03中使用Counterpoise=2关键字,会计算5个体系,输出的能量按照如下顺序:
E(AB):A、B基组下AB复合物的能量
E(A,bAB):A、B基组下A的能量
E(B,bAB):A、B基组下B的能量
E(A):A基组下A的能量
E(B):B基组下B的能量
最后会输出corrected energy (E(corrected))和BSSE energy (E(BSSE)). 这里E(corrected)就是消除了因单体基组重叠造成的能量降低后的AB复合物能量,E(corrected) = E(AB) + E(BSSE)。E(BSSE) = (E(A) - E(A,bAB)) + (E(B) - E(B,bAB))。
BSSE校正后的真实的相互作用能可以这样计算:E(interaction) = E(corrected) - (E(A) + E(B))。也可以这样计算,是等价的:E(interaction) = E(AB) - E(A,bAB) - E(B,bAB)
计算过程中会输出类似这样的语句:"Counterpoise: doing DCBS calculation for fragment 1."这里就是说明接下来计算的是E(A,bAB)(假设A分子为fragment 1),其中DCBS代表dimer centered basis set,说明以A、B分子为中心的基函数都出现,但是计算中并不纳入B的电子和原子核,这称为计算A的能量时添加了B的ghost轨道;如果是"doing MCBS calculation for fragment 1",就是要计算E(A),MCBS代表monomer centered basis set,计算中只出现属于A分子的基函数。
若计算n个分子,则关键字为counterpoise=n,结果输出顺序与计算相互作用能的方法与双分子的情况是一样的。能量按如下顺序输出:E(AB),E(1)',E(2)'...E(n)',E(1),E(2)...E(n)。E(BSSE) = E(1) - E(1)' + E(2) - E(2)' + ... + E(n) - E(n)'。E(corrected) = E_complex + E(BSSE)。E(interaction) = E(corrected) - ( E(1) + E(2) + ... + E(n) )。计算过程中也用DCBS和MCBS来说明接下来将要计算的是哪项,但此时DCBS中的D的含义就不是具体指Dimer了,而是多分子复合物。
(2.17)一些特殊的DFT在Gaussian03中如何输入?引用回帖: 我所谓的特殊的指的是G03中不存在的DFT,但是可以通过用户自定义来得到的DFT。Gaussian03可以使用具有下列一般形式的任何模型
P1*[P3*ΔEx(non-local)+P4*Ex(Slater)]+P2*Ex(HF)+P5*ΔEc(non-local)+P6*Ec(local)
P1~P6的值由IOP(3)来指定
一套指定方式,支持四位小数
IOp(3/45=mmmmnnnn)
指定P1为mmmm/1000,P2为nnnn/1000.
IOp(3/46=mmmmnnnn)
指定P3为mmmm/1000,P4为nnnn/1000。
IOp(3/47=mmmmnnnn)
指定P5为mmmm/1000,P6为nnnn/1000。
另一套指定方式,支持五位小数
IOp(3/76=MMMMMNNNNN)
指定P1为MMMMM/10000,P2为NNNNN/10000.
IOp(3/77=MMMMMNNNNN)
指定P3为MMMMM/10000,P4为NNNNN/10000.
IOp(3/78=MMMMMNNNNN)
指定P5为MMMMM/10000,P6为NNNNN/10000.
举个例子
# BLYP IOp(3/45=10000200) IOp(3/46=07200800) IOp(3/47=08101000)
# BLYP IOp(3/76=1000002000) IOp(3/77=0720008000) IOp(3/78=0810010000)
这两个输入路径其实都相当于B3LYP
(1)mpw1k/6-31+G(d,p)
# mpwpw91/6-31+G(d,p) iop(3/76=0572004280)
(2)MPWB1K/6-31+G(d,p)
# mpwb95/6-31+G(d,p) IOp(3/76=0560004400)
(3)MPW1B95/6-31+G(d,p)
# mpwb95/6-31+G(d,p) IOp(3/76=0690003100)
(4)TPSS1KCIS/6-31+G(d,p)
# tpsskcis/6-31+G(d,p) IOp(3/76=0870001300)
(5)MPW3LYP/6-31+G(d,p)
# mpwlyp/6-31+G(d,p) IOp(3/76=1000002000) IOp(3/77=0720008000) IOp(3/78=0810010000)
(6)MPWKCIS1K/6-31+G(d,p)
# MPWKCIS/6-31+G(d,p) IOp(3/76= 0590004100)
(7)MPW1KCIS/6-31+G(d,p)
# MPWKCIS/6-31+G(d,p) IOp(3/76= 0850001500)
(8)PBE1KCIS/6-31+G(d,p)
# pbekcis/6-31+G(d,p) IOp(3/76=0780002200)
(9)PBE1W/6-31+G(d,p)
# pbepbe/6-31+G(d,p) IOp(3/78=0740010000)
(10)MPWLYP1W/6-31+G(d,p)
# mpwv5lyp/6-31+G(d,p) IOp(3/78=0880010000)
(11)PBELYP1W
# pbev5lyp/6-31+G(d,p) IOp(3/78=0540010000)
(12)TPSSLYP1W
# tpssv5lyp/6-31+G(d,p) IOp(3/78=0740010000)
(13)MPWLYP1M/6-31+G(d,p)
# mpwlyp/6-31+G(d,p) IOp(3/76=0950000500)
(14)MOHLYP/6-31+G(d,p)
# OV5LYP/6-31+G(d,p) IOp(3/77=1292010000) IOp(3/78=0500010000)
(15)MOHLYP2/6-31+G(d,p)
# OV5LYP/6-31+G(d,p) IOp(3/77=1849810515) IOp(3/78=0500005000) (for Gaussian03 B&C)
# OV5LYP/6-31+G(d,p) IOp(3/77=1292010000) IOp(3/78=0500005000) (for Gaussian03 D&E)
小卒版主说,这个问题的答案,以后有时间会继续扩充。先做个标记
(2.18)我们在用gauss寻找过渡态之前应该准备什么?(不断更新中)引用回帖:
过渡态的计算往往是Gaussian计算的难点之一,我做过多次过渡态的计算,甚至有两个课题就是专门研究反应过程的。所以给出一点浅见,请大家斧正。
过渡态的计算不再是黑匣子,需要人的参与很多。尤其是不能像基态优化那样“蛮干”了。我认为这是初学者算过渡态很难成功的原因之一。
个人建议,在做过渡态反应计算的之前,应该确认三件事情
(1)这个反应是基元反应吗?如果是,那这个反应会有几种可能的过渡态?如果不是,我考虑到了所有可能的中间态和过渡态了吗?不论它是不是基元反应,我考虑了它的全部可能的反应路径了吗?
(2)进攻位点对不对?
(3)进攻位点对了,那么我的进攻方向对不对?进攻分子和目标分子的相对位置(角度、距离)对不对?
比如,我在教我的师弟师妹的时候,举了一个过渡态计算的例子。【图8-1】过渡态的计算是很繁杂的。一个小小的DA反应,就要计算这么多东西——考虑各种中间态和过渡态,考虑各种自旋状态,考虑各种可能的进攻位点和进攻方向。见附件
初学者往往搞不懂过渡态和中间态的区别。我做一个很简单的解释:
(1)过渡态有虚频,中间态没有虚频
(2)如果一个反应路径中有中间态,那这个反应路径一定不是基元反应
(3)如果一个反应路径中只有一个过渡态,没有中间态,那这个反应路径一定是基元反应
(4)过渡态是反应物,生成物和中间态的联结纽带。反应物,生成物和中间态之间,必须由过渡态来连接
(请注意我在这个帖子中,严格区分了“反应”和“反应路径”两个概念。一个反应可能有多个反应路径。比如在附件中的DA反应,我就设计了3条反应路径,分别是:1-2-3-4-5,1-6-7-8-9-10-5,1-2-3-8-9-10-5)
每一个反应路径算完之后,都要给出irc验证。当然如果你能做二维甚或三维势能面扫描更好,一切看你能使用的计算资源。请注意,如果说irc计算量是N,那么M维势能面扫描所需的计算量是N^M
事实上JACS1999,121,4816-4826和JOC2004,69,3683-3692这两篇文章已经确定,附件中的1-2-3-4-5这个反应路径是最低能量反应路径,我只是拿DA反应的计算做个例子
===========
2009.9.23,majun04虫友补充道:
补充一点具体的问题:先把过渡态两边的结构(也可以说山谷的点)确定准确,再涉及过渡态的问题。过渡态两端链接的是反应物,产物,或者反应中间体。在一条路径上往往是多个过渡态的(也就是要翻过多个山峰到达终点)。在计算过渡态之前的反应物(或者中间体)的几何结构是非常重要的,对于这个结构是要注意的。你给出的结构是否是实验测量得到的,或者是文献报道的,如果是,证明你的“起点”具有一定的可靠性。如果没有文献报道,是你计算模拟出来的,这就要对这一点进行详细的阐述。需要从结构,配位,从电荷分布等方面阐述,具有分子间氢键的还要加上分子间氢键。如果是超分子结构,小分子要使用MP2或者HF辅助DFT方法进行验证。如果是大分子,需要考虑BSSE。把这些阐述明了,再进行过渡态的寻找。否则审稿人也会对你的结构提出质疑的。
谢谢你的具体问题补充。我写的东西确实是高度概括、流于笼统了
http://muchong.com/bbs/viewthread.php?tid=1491020&fpage=1 如果想进一步学习,请看这个帖子
(2.19)关于为什么使用QST3的时候必须优化反应物和产物?引用回帖: 今天QQ上有个虫友找我,说做QST3计算,但是没有优化反应物和产物的构型,就用这种未优化的初始构型和终止构型计算了TS构型。课题基本做完了,问我这样行不行。
我虽然很替他惋惜,但是我不得不说:“非常可惜,你的计算确实没有意义”
QST2和QST3是基于QST过渡态理论的计算方法。这要从QST计算方法的原理上来讲一讲了。
QST的原理如图【8-2】,A与B是势能面上两个“碗底”。我打算用QST方法找一个irc线。于是就先在A和B之间画一条直线,然后,逐步沿着这条直线上的垂线方向的力常数,开始偏移这条直线,直到找到一阶鞍点TS和合适的irc线。
如果A和B的构型本身就没有被优化,如A'和B'这两个点,那么QST也将忠实的按照这条直线上的垂线方向的力常数,开始偏移这条直线,直到找到irc'。但是,这个irc'与irc是势能面上不同的线。irc'和TS'并不是irc&TS。
如果你的产物构型和反应物构型搭建的很不错,很接近优化后的构型,那么,如果你不优化产物构型和反应物构型,你也能得到合理的TS。比如图中的A"和B"。但是,你得到的irc"仍旧是错的。
讲述完毕。心情沉痛,教训啊,很替他惋惜!大家一定要引以为戒
(2.20)什么是极化连续介质模型(PCM)?引用回帖: Tomasi和他的同事提出的极化连续介质模型(Polarizable Continuum Model, PCM)是一个经常用到的连续溶剂化方法,并且这些年来已经有很多的改进。PCM模型计算分子溶液中的分子自由能为三项的加和:
G(sol)=G(es)+G(dr)+G(cav)
这些成分代表了静电的(es)和散射-推斥对自由能(dr)的贡献,和空穴能(cav)。所有这三项都是由以原子位置为中心的连锁范德华球定义的空穴来计算的。反应场是通过位于分子空穴表面的点电荷(表观表面电荷模型)表示的。这里要讨论的PCM的特殊版本是一种用为Hartree-Fock的联合原子模型来构建空穴。在这种模型中,范德华表面是只由定位于重元素(非氢元素)的一些球构建的联合原子方法。每个原子的范德华半径是原子类型,连接性,分子总电荷,和连接氢原子数目的一个函数。在评估方程(1)中的三项时,这个空穴的用法稍微有些不同。【图8-3】
当计算空穴能Gcav时,采用范德华球定义的表面,溶剂可接近表面被用来计算散射-推斥对自由能(dr)的贡献。后一表面不同于前一表面,在后者中额外考虑了(理想化的)溶剂半径。在溶液中静电对自由能的贡献 Ges采用一个近似版本的溶剂排斥面,该排斥面通过用一个常数因子约化所有的半径,并在之后加入更多的不以原子为中心的球,以便得到稍微平滑些的表面。
定域化和表面电荷的计算是通过系统地将球表面分割成已知面积的镶嵌块计算每个表面元素的一个点电荷来达到的。
在Gaussian98中执行PCM/UAHF模型可以用SCRF关键词结合PCM专用修饰词。 溶剂可以用对SCRF关键词给出Solvent= modifier来指定,可接受的溶剂名称是Water(水), DMSO(二甲亚砜), NitroMethane(硝基甲烷), Methanol(甲醇), Ethanol(乙醇), Acetone(丙酮), DiChloroEthane(二氯乙烷), DiChloroMethane(二氯甲烷), TetraHydroFuran(四氢呋喃), Aniline(苯胺), ChloroBenzene(氯苯), Chloroform(氯仿), Ether(乙醚), Toluene(甲苯), Benzene(苯), CarbonTetrachloride(四氯化碳), Cyclohexane(环己烷), Hepaten, and Acetonitrile(乙腈)。 附加的选项可以在输入文件的末尾指定并用给SCRF关键词指定Read修饰词来读入。PCM溶剂化模型可以用于HF和DFT水平下的能量计算和梯度计算。 PCM计算产生的输出可以用DUMP选项显著地扩展。
下面的例子输入解释Cs对称性的乙醇的水溶剂化自由能的单电能计算(没有几何优化,【图8-4】):
#P B3LYP/6-31G(d) scf=tight int=finegrid SCRF=(PCM,Read,Solvent=Water)
pcm/b3lyp/6-31G(d) sp ethanol in water (Cs)
0 1
O1
C2 1 r2
C3 2 r3 1 a3
H4 3 r4 2 a4 1 180.0
H5 3 r5 2 a5 4 d5
H6 3 r5 2 a5 4 -d5
H7 2 r7 3 a7 1 d7
H8 2 r7 3 a7 1 -d7
H9 1 r9 2 a9 3 180.0
r2=1.42492915
r3=1.51965095
r4=1.09569807
r5=1.09496362
r7=1.10264669
r9=0.96904984
a3=107.81130783
a4=110.63999342
a5=110.37205263
a7=109.90077195
a9=107.87777748
d5=-120.23659087
d7=-121.12750852
DUMP
由PCM溶剂化模型引起的附加输出是由负责SCF计算的L502产生的:
-------------------------------------------------------------------
Solvent: WATER
Model : PCM/UAHF, Icomp = 4
Version: MATRIX INVERSION
Cavity : PENTAKISDODECAHEDRA with 60 initial tesserae
-------------------------------------------------------------------
Nord Group Hybr Charge Alpha Radius Bonded to
1 OH sp3 0.00 1.20 1.590 C2 [s]
2 CH2 sp3 0.00 1.20 1.860 O1 [s] C3 [s]
3 CH3 sp3 0.00 1.20 1.950 C2 [s]
-------------------------------------------------------------------
------------------------------------------------------
Dielectric Const = 78.39000
High.Fr.D.Const = 1.77600
d(Diel.Const.)/dT = -0.35620
Molar Volume = 18.07000
Therm.Exp.Coeff. = 0.00026
Radius = 1.38500
Absolute temper. = 298.00000
Number of spheres = 3
OMEGA = 40.00000
RET = 0.20000
FRO = 0.70000
Accuracy = 0.1D-05
------------------------------------------------------
头四行重复对水指定的设置或通常的PCM缺省设置。 溶质空穴由范德华球构成,范德华球是由规则的pentakisdodecahedra(五个十二面体——六十面体?)表示的,将每个球表面分成相同尺寸的60个单元。
后面的四行列出了UAHF分析的结果,确定只有三个中心(联合原子)。对于每个中心,假设的杂化和它的形式电荷,最后半径及溶剂专用约化参数Alpha一起列出。 后者通常用缺省值1.2, 但是也可以用选项
ALPHA=x.x
直接指定。
输出的最后部分列出了溶剂专用参数如介电常数和有效溶剂半径,和一些更多的缺省PCM设置,如初始球的数目和参数OMEGA, RET,和FRO的当前值。后面的这三个参数控制着渐入更多球(非以原子位置为中心)的过程,以便平滑表面。新的OMEGA值可以用
OMEGA=n,n
选项设置。 有意义的值在40.0-90.0之间(较高的值给出较少的加入球)。新的FRO值可以用
FRO=m,m
指定。有意义的值在0.7到0.2之间(较小的值给出较少的加入球)。RET指定加入新球的最小半径,新值可以用
RET=1.1
指定。增加的数值给出较少的附加球,非常大的值完全消除附加球。
PCM算法是首先进行一步气相能量计算,以便得到后面溶剂化自由能计算的参考点。在结束气相SCF循环后,列出了空穴产生的迭代过程的细节:
------------------------------------------------------
---------- CAVITY for ELECTROSTATIC term -----------
------------------------------------------------------
------- The SOLUTE is enclosed in ONE CAVITY -------
Total N.of Tesserae = 132
Surface Area (Ang**2) = 97.71529
Volume (Ang**3) = 84.60281
------------------------------------------------------
Original Sphere On Atom Re0 Alpha Surface
1 O1 1.590 1.200 24.08112
2 C2 1.860 1.200 27.26919
3 C3 1.950 1.200 46.36498
------------------------------------------------------
----------------------------------------------------------------------
AT CONVERGENCE
132 Tesserae over a maximum of 1500
Surface Area (Ang**2) = 97.71529
Volume (Ang**3) = 84.60281
Escaped Charge= 0.13334
Error on NUCLEAR pol.charges = 0.21898 Error on ELECTR. pol.charges =-0.33812
----------------------------------------------------------------------
------------------------------------------------------------------------------
dG(solv)/dEps (kcal/mol) = 0.00000
------------------------------------------------------
IN VACUO Dipole moment (Debye):
X= 0.0176 Y= 1.5625 Z= 0.0000 Tot= 1.5626
IN SOLUTION Dipole moment (Debye):
X= 0.1053 Y= 1.9429 Z= -0.0019 Tot= 1.9457
------------------------------------------------------
Tessera X Y Z QTot QSN QSE
1 2.84136 0.54870 3.42913 0.00354 -0.17977 0.18331
.
.
.
132 -2.70116 -2.26783 -4.13447 -0.00393 -0.25391 0.24997
在这个例子(一切正常的)中不需要附加球。 总的表面用132个镶嵌块("Tesserae" 表示。由于电子波函数的长尾(事实上,没有结尾),用当前表面定义的分子体积不包含体系所有电子密度,导致一些“逃逸电荷”。在每个表面元中心有一个表面电荷“QTot”,它含有一个从溶质核电荷而来的组分和一个从溶质电子电荷而来的组分。
在解电子薛定谔方程(包括附加的反应场效应)的求解过程得到的结果给出如下:
SCF Done: E(RB+HF-LYP) = -155.041616090 A.U. after 10 cycles
Convg = 0.1222D-08 -V/T = 2.0093
S**2 = 0.0000
KE= 1.536131750668D+02 PE=-5.252140797608D+02 EE= 1.349508863800D+02
------------------------------------------------------
-------------- VARIATIONAL PCM RESULTS -------------
------------------------------------------------------
(a.u.) = -155.033805
(a.u.) = -155.040760
(a.u.) = -155.041609
(a.u.) = -155.032883
(a.u.) = -155.041616
Total free energy in sol.
(with non electrost.terms) (a.u.) = -155.041162
------------------------------------------------------
(Unpol.Solute)-Solvent (kcal/mol) = -4.36
(Polar.Solute)-Solvent (kcal/mol) = -5.48
Solute Polarization (kcal/mol) = 0.58
Total Electrostatic (kcal/mol) = -4.90
------------------------------------------------------
Cavitation energy (kcal/mol) = 8.92
Dispersion energy (kcal/mol) = -11.39
Repulsion energy (kcal/mol) = 2.75
Total non electr. (kcal/mol) = 0.29
------------------------------------------------------
DeltaG (solv) (kcal/mol) = -4.62
------------------------------------------------------
这里列出的:
(a.u.) = -155.033805
是非微扰的气相SCF解,它被用作后面所有步骤的参考。后面描述为:
(a.u.) = -155.040760
的能量包括非极化的溶质和非极化的溶剂之间的相互作用。和气相参考能量对比得到相应的作用能:
(Unpol.Solute)-Solvent (kcal/mol) = -4.36
报告对应于非极化的溶质和极化的溶剂相互作用的总能量之后,下面的重要信息是关于极化的溶质的总能量。
(a.u.) = -155.032883
与非极化的气相总能量的能量差列出如下:
Solute Polarization (kcal/mol) = 0.58
而且永远应该是正的。极化的溶质和极化的溶剂的全部作用体系给出总能量:
(a.u.) = -155.041616
相应的静电作用能列出如下:
Total Electrostatic (kcal/mol) = -4.90
溶剂化能的非静电部分在同一个block里给出,以空穴化和散射-推斥能的加和结束:
Total non electr. (kcal/mol) = 0.29
非静电和静电贡献的加和给出了总的溶剂化自由能:
DeltaG (solv) (kcal/mol) = -4.62
然而,需要认识到,这里说的溶剂化自由能是指温度为0 K的气相静止体系。为了在特定温度下得到热力学有意义的溶剂化自由能,这些溶剂化自由能必须用气相热化学的标准处理进行修正。 对于计算溶剂化自由能自身,作为气相和溶液自由能的差值,这一步经常是省略的,PCM得出的数值被直接用来和实验值相比。对乙醇来说,实验溶剂化自由能已经被测量为-5.0kcal/mol。和这个值相比PCM预测的-4.6kcal/mol可以被认为是相当精确的。
当溶剂化能作为气相和溶液相自由能的差计算的时候,必须注意到它们各自标准态的定义。因为气相和液相浓度都是以mol值(mol/l)给出的,气相和液相的数据可以直接对比。然而,气相值经常是指分压为1atm。假设是理想气体行为,在298.15K这对应于1/24.46 mol/l。
改进的溶剂化自由能预测应该涵盖溶液中的结构松弛效应。用PCM模型的几何优化是可能的,但是比气相优化消耗多许多的时间。这不仅是由于每次的能量和提速计算需要更多的CPU时间,而且是由于优化过程的收敛缓慢和时常出现的振荡。有两个选项对于减轻一些收敛问题是有用的,它们是TSNUM和TSARE:
TSNUM指定每个球的表面元的数目。 PCM算法选择规则的多面体,该多面体的表面元的数目尽可能地接近TSNUM。除了缺省的数目60,一些更大的值64,80,或100或许对于减轻一些几何优化中的共振行为是有帮助的。
TSARE以单位(Angstrom)2指定表面元的面积。有意义的值的范围从0.4到0.2,越小的值导致越大的表面元数目。将表面元的尺寸设为特定的值导致等尺寸的表面元,而不考虑球的半径(对于TSNUM设置改变的情况并非如此)。
然而,在这两种情况中,总能量依赖于表面元的实际选择,对不同体系或不同异构体的对比只有在选项选择相同的时候才有意义。对于这里用到的乙醇的例子,几何优化在用缺省设置时在23圈仍没收敛,但是在用TSARE=0.3时在7圈内收敛,得到最后总的溶剂化能-155.043097702 au.和用气相几何作的PCM的单点能计算相比(总能量为-155.041371 au),这暗示一个-1.1kcal/mol的结构松弛能,因此对溶剂化能“改进的”的预测是-5.7kcal/mol。
在应用当前的PCM/UAHF模型于溶液中反应路径时的一个问题是为范德华半径的导出使用杂化和连结性的直接后果。因为杂化和连结性都不是平滑变化的,而是沿着反应路径突然地从一点到另一点,UAHF方法一定会在溶剂化自由能图上产生突然的断点。这些问题可以,从原理上,通过平滑地约化从一组半径到另外一组, 或者通过完全地避免UAHF方法,选择与连结性无关的半径。 对后者,一个通常的选择是Pauling半径,用选项RADII=Pauling可以做到。 (2.21)物化中学过反应的ΔG大于40千焦的话就认为不反应。但有时候可以发生的反应,过渡态的能量和反应物的能量算出来会差几十千卡(>40kj/mol),请问应该怎么理解?引用回帖: 首先要搞清楚热力学判据与动力学判据的区别。
所谓热力学判据,只考虑反应物和产物的平衡比例,不考虑具体实现过程,这个时候应用范特霍夫方程ΔG=-RTlnK,此时的ΔG是指产物与反应物的吉布斯自由能之差,这里大写的K是反应平衡常数。假设反应为简单的A->B,反应温度为298.15K(25℃),达到平衡时,假设B的浓度是A的一千万倍(这个时候我们可以肯定的说,A完全转换为B了吧。而且如果以B为反应物,我们说B->A是不反应的,也不会有人有意见了吧),那么,我们计算一下,ΔG=8.314*298.15*ln10000000=39954J/mol≈40kJ/mol≈10kcal/mol
如果放松点标准,生成B的浓度是A的浓度的100倍,那么此时ΔG=11.4kJ/mol=2.7kcal/mol
而动力学判据,考虑的是反应进行的具体路线,不同条件下反应机理是不一样的,反应速率是不一样的,动力学此时应用阿雷尼乌斯公式k=Aexp(-Ea/RT),此处的反应活化能Ea=ΔH≠(≠写成上角标)+nRT,由于一般化学反应的活化能或活化焓都在40-80kJ/mol,而室温下RT仅为2.5kJ/mol,所以多数时被忽略,这个时候活化能或者活化焓指的是过渡态与反应物之间的能量差或者焓变,而不是热力学上说的ΔG。由于指前因子A取决于反应的熵变(具体利用波尔兹曼分布推算,此略),所以也可以考虑把熵一并放到活化能里面,变成吉布斯自由能项,这个时候就可以指定和计算过渡态与反应物之间的ΔG,但这个时候一遍写做ΔG≠(≠写成上角标)。数值的变化不大(除非熵变特别明显),所以大概的范围还是20-30kcal/mol(低于5kcal则反应太快,没有办法可以检测到反应过程,超过40kcal则反应太慢,以至于我们看来相当稳定,几乎就是不反应了),你说计算出来的能量差在20+kcal,这个应该很正常——虽然你没有说到底是ΔH≠还是ΔG≠(注意能量单位,不要一会用千焦一会用千卡,自己习惯用哪个就一直用哪个就好)
再举一个简单的例子来说明热力学判据和动力学判据的差别:
氢气和氧气燃烧放热,典型的不可逆反应,这个用热力学判据可以解释。
然而你把氢气和氧气在常温常压下混合在一起,难道它们会自燃么?如果你不改变条件,放100年它也不反应生成水,这就是动力学上无法翻越那个势垒,只有改变了外界条件,它才会反应。
(2.22)什么是电子密度?什么是电子密度差?它们的物理含义是什么?如何得到?引用回帖: 电子密度其实就是ρ=Ψ*Ψ,ρ从Ψ的平方而来,它是从第一原理得到结论性数据,这一点与NBO、AIM、Mulliken等电荷分配有天渊之别。
电子密度的变化,就是电子密度差,Δρ=ρ2-ρ1,它表示的就是变化前后,空间中各点的电子密度分布的变化
如何得到电子密度和密度差?需要调用cubegen和cubman命令。Gaussian中cubman的操作方法,见下面的例子:
以两种不同状态的乙烯分子为例,介绍如何获得它们的电子密度分布的差别。大家都是用Gaussian的,我就不多解释了,直接上gjf文件吧。如果gjf文件看不懂,就别在这儿混了。
====gjf文件开始===
%chk=G:\sundries\C2H4\C2H4-L.chk
# hf/sto-3g
Low spin
0 1
C 0.00000000 0.00000000 0.67760000
H 0.92664718 0.00000000 1.21260000
H -0.92664718 0.00000000 1.21260000
C 0.00000000 0.00000000 -0.67760000
H 0.92664718 0.00000000 -1.21260000
H -0.92664718 0.00000000 -1.21260000
--Link1--
%chk=G:\sundries\C2H4\C2H4-H.chk
# hf/sto-3g
High spin
0 3
C 0.00000000 0.00000000 0.67760000
H 0.92664718 0.00000000 1.21260000
H -0.92664718 0.00000000 1.21260000
C 0.00000000 0.00000000 -0.67760000
H 0.92664718 0.00000000 -1.21260000
H -0.92664718 0.00000000 -1.21260000
====gjf文件结束====
运行这个gjf文件,获得C2H4-H.chk和C2H4-L.chk两个文件之后,打开DOS(这个DOS指的是Windows的命令行操作系统Disk Operation System)红色字体是我输入的,黑色字体是DOS的提示
====DOS操作开始====
Microsoft Windows XP [Version 5.1.2600]
(C) Copyright 1985-2001 Microsoft Corp.
C:\Documents and Settings\Administrator>d:
D:\>cd programfiles\gaussian03w
D:\ProgramFiles\gaussian03w>formchk G:\sundries\C2H4\C2H4-H.chk G:\sundries\C2H4\C2H4-H.fchk
Read checkpoint file G:\sundries\C2H4\C2H4-H.chk
Write formatted file G:\sundries\C2H4\C2H4-H.fchk
Rotating derivatives, DoTrsp=T IDiff= 1 LEDeriv= 353 LFDPrp= 0 LDFDPr= 0.
D:\ProgramFiles\gaussian03w>formchk G:\sundries\C2H4\C2H4-L.chk G:\sundries\C2H4\C2H4-L.fchk
Read checkpoint file G:\sundries\C2H4\C2H4-L.chk
Write formatted file G:\sundries\C2H4\C2H4-L.fchk
Rotating derivatives, DoTrsp=T IDiff= 1 LEDeriv= 353 LFDPrp= 0 LDFDPr= 0.
D:\ProgramFiles\gaussian03w>cubegen 0 density G:\sundries\C2H4\C2H4-H.fchk G:\sundries\C2H4\C2H4-H-D.cube 0 h
D:\ProgramFiles\gaussian03w>cubegen 0 density G:\sundries\C2H4\C2H4-L.fchk G:\sundries\C2H4\C2H4-L-D.cube 0 h
D:\ProgramFiles\gaussian03w>cubman
Action [Add, Copy, Difference, Properties, SUbtract, SCale, SQuare]? subtract
First input? G:\sundries\C2H4\C2H4-H-D.cube
Is it formatted [no,yes,old]? yes
Opened special file G:\sundries\C2H4\C2H4-H-D.cube.
Second input? G:\sundries\C2H4\C2H4-L-D.cube
Is it formatted [no,yes,old]? yes
Opened special file G:\sundries\C2H4\C2H4-L-D.cube.
Output file? G:\sundries\C2H4\C2H4-LH-Different.cube
Should it be formatted [no,yes,old]? yes
Opened special file G:\sundries\C2H4\C2H4-LH-Different.cube.
Input file titles:
High spin density
Electron density from Total SCF Density
Input file titles:
High spin densit
Electron density from Total SCF Density
SumAP= 12.0019420655 SumAN= 0.0000000000 SumA= 12.0019420655
CAMax= 2.1672300000 XYZ= -0.0360500000 -0.0360500000 1.3091610000
CAMin= 0.0000000000 XYZ=-9999.0000000000-9999.0000000000-9999.0000000000
SumBP= 12.0019319790 SumBN= 0.0000000000 SumB= 12.0019319790
CBMax= 2.1577300000 CBMin= 0.0000000000
SumOP= 0.1979999629 SumON= -0.1979898763 SumO= 0.0000100865
COMax= 0.0428610000 COMin= -0.0297300000
DipAE= 0.0002322172 0.0002783966 -0.0003923874
DipAN= 0.0000000000 0.0000000000 0.0000000000
DipA= 0.0002322172 0.0002783966 -0.0003923874
DipBE= 0.0002261708 0.0002769231 -0.0003922251
DipBN= 0.0000000000 0.0000000000 0.0000000000
DipB= 0.0002261708 0.0002769231 -0.0003922251
DipOE= 0.0000060464 0.0000014735 -0.0000001623
DipON= 0.0000000000 0.0000000000 0.0000000000
DipO= 0.0000060464 0.0000014735 -0.0000001623
D:\ProgramFiles\gaussian03w>
====DOS操作结束====
用GaussView打开C2H4-LH-Different.cube,生成附件中的三个图(图8-5),就OK了
(2.23)各方法的红外校正因子是多少引用回帖:
[ Last edited by yjcmwgk on 2010-8-21 at 14:26 ]