24小时热门版块排行榜    

查看: 9319  |  回复: 0
当前主题已经存档。

wqd198686

铁杆木虫 (正式写手)

[交流] 总结:高斯量化计算总结(必看)

关于自旋多重度

定义: 多重度=2S+1, S=n*1/2,n为单电子数。所以,关键是单电子的数目是多少。

当有偶数个电子时,例如 O2,共有16个电子,那么单电子数目可能是0,即8个alpha和8个beta电子配对,对应单重态,但是也可能是有9个α电子和7个β电子,那么能成对的是7对,还剩2个α没有配对,于是n=2,对应的是多重度3。同理还可以有多重度5,7,9, ...一般而言,是多重度低的能量低,最稳定,所以,一般来说,偶数电子的体系多重度就是1。但是也有例外,例如O2就是一个大家都知道的例子,它的基态是三重态,其单重态反而是激发态。所以对于未知的体系,还是算几个保险一点,看哪个能量更低。

所以,总结一下,就是

电子数目是偶数,未成对电子数目n=0,2,4,6,...自旋多重度是1,3,5,7,...
电子数目是奇数,未成对电子数目n=1,3,5,7,...自旋多重度是2,4,6,8,...
多数情况是多重度低的能量低,有时(特别是有“磁”性的时候,例如顺磁的O2,以及Fe啊什么的),可能会高多重度的能量低,所以需要都算算,看哪个能量更低。
关于赝势:

简单来说,赝势就是不计算内层电子,而是把内层电子的贡献用一个势来描述,放在哈密顿里面。适用于重元素。

赝势基组,实际上包括 赝势和基组两个部分,内层电子采用赝势,即effective core potential (ECP),外层价电子采用一般的基组。 比如:
LanL2DZ: D95V on first row, Los Alamos ECP plus DZ on Na-Bi.
就是对第一行原子是 D95V (这个是非赝势基组),对Na-Bi是使用一个叫做Los Alamos的有效核势加上一个DZ基组。所以Lanl2dz就是对前面的原子全电子基组,对后面的原子是赝势基组。 (再次说明,量化里面,C,O那一行,算周期表的第一行)

使用赝势的3个原因:
1。没有相应的全电子基组。
2。减少计算量。
3。赝势可以包含重金属相对论效应的修正。

在高斯中,lanl2dz基组,在手册中可以查到其定义为:
LanL2DZ: D95V on first row , Los Alamos ECP plus DZ on Na-Bi,也就是说,对于C,O等元素来说(量化中,H大概算第0周期,C,O才是第一周期),Lanl2dz实际上还是全电子基组,而对于Na以后才是对内层电子用Los Alamos ECP赝势,外层电子用DZ 基组。

使用赝势的输入文件:

1.所有原子使用lanl2dz
--------------------------
#HF/lanl2dz opt

lanl2dz for all atoms

0  2
O  0.0 0.0 0.0
C  0.0 0.0 1.2
Cu 0.0 0.0 3.2

--------------------------

2.所有原子使用lanl2dz的另一种输入方法。
--------------------------
#HF/genecp opt

lanl2dz for all atoms

0  2
O  0.0 0.0 0.0
C  0.0 0.0 1.2
Cu 0.0 0.0 3.2

C O 0
lanl2dz
****
Cu 0
lanl2dz    //定义价电子的基组, C O 0 是碳,氧,零,其中零用作终止符号。
****

Cu 0
lanl2dz   //定义内层电子的赝势

--------------------------

3.混和基组,即有的使用全电子,有的使用Lanl2dz。格式同2。
--------------------------
#HF/genecp opt

lanl2dz for Cu, 6-31G(d) for C and O

0  2
O  0.0 0.0 0.0
C  0.0 0.0 1.2
Cu 0.0 0.0 3.2

C O 0
6-31G(d)   //另一种全电子基组
****
Cu 0
lanl2dz    //定义价电子的基组
****

Cu 0
lanl2dz   //定义内层电子的赝势

--------------------------
开壳层和闭壳层

闭壳层计算就是对于多重度是1的体系,此时α和β的电子数目相同,可以把α和β配对,成对的α和β使用同一个轨道,一个轨道上填充2个电子。
开壳层计算就是对α和β电子分别计算,一个轨道上只填充1个电子,一般来说,多重度是1时,开壳层计算和闭壳层计算会给出相同的结果。
限制开壳层计算是对多重度大于1的体系,此时α和β的电子数目不同,设有m个α和n个β电子,m>n,那么让前n个轨道上每个填充一个α和一个β,剩下的m-n个α电子再填充m-n个轨道。即前n个轨道是闭的(每个轨道2个电子),后m-n个轨道是开的(每个轨道1个电子)

在高斯中,以HF为例,闭、开、限制开壳层计算分别是RHF,UHF,ROHF。如果只写HF,则按下面的方式取默认方法:
对多重度是1的体系,默认为RHF,对多重度大于1的体系,默认是UHF。
关于收敛问题 (L502, L508, L9999)

对于一个优化计算,它的过程是先做一个SCF计算,得到这个构型下的能量,然后优化构型,再做SCF,然后再优化构型。。。因此,会有两种不收敛的情况:一是在某一步的SCF不收敛(L502错误),或者构型优化没有找到最后结果(L9999错误)。

预备知识:计算时保存chk文件,可以在后续计算中使用guess=read读初始猜测.

对于SCF不收敛,通常有以下的解决方法:
1. 使用小基组,或低级算法计算,得到scf收敛的波函数,用guess=read读初始波函数。
2. 使用scf=qc,这个计算会慢,而且需要用stable关键字来测试结果是否波函数稳定。如果这个还不收敛,会提示L508错误。
3. 改变键长,一般是缩小一点,有时会有用。
4. 计算相同体系的其他电子态,比如相应的阴离子、阳离子体系或单重态体系,得到的收敛波函数作为初始猜测进行计算。
5. 待补充.

对于优化不收敛,即L9999错误,实际上是在规定的步数内没有完成优化,即还没有找到极小值点。(或者对于过渡态优化,还没有找到过渡态)

这有几种可能性:
1. 看一下能量的收敛的情况,可能正在单调减小,眼看有收敛的趋势,这样的情况下,只要加大循环的步数(opt(maxcycle=200)),可能就可以解决问题了。
2. 加大循环步数还不能解决的(循环步数有人说超过200再不收敛,再加也不会有用了,这虽然不一定绝对正确,但200步应该也差不多了),有两种可能。一是查看能量,发现能量在振荡了,且变化已经很小了,这时可能重新算一下,或者构型稍微变一下,继续优化,就可以得到收敛的结果(当然也有麻烦的,看运气和经验了);二是构型变化太大,和你预计的差别过大,这很可能是你的初始构型太差了,优化不知道到哪里去了,这时最好检查一下初始构型,再从头优化。
3. 对于L9999快达到收敛时,考虑减小优化步长有时对于能量振荡的情况也是有用的,opt(maxstep=1).(flyingheart )

一个建议是,对于大体系,难收敛体系,先用小基组,低精度算法优化一下,以得到较好的初始构型,再用高精度的计算接着算。如果前面的方法保留了chk文件,重新计算时需要使用 geom=allcheck 读入构型(就不必麻烦地写构型了), guess=read(读入初始波函数,可以加快第一步SCF收敛)。
关于对称性:
(原贴 http://210.34.15.126/cgi-bin/topic.cgi?forum=3&topic=1736

刚刚在zixia上看到zhuoliu的大作。结合自己以前知道的东西,总结一下。

1。gaussian中输入什么对称性,一般优化的结果仍然还是那个对称性,比如CO2,如果初始两个CO键长输入不是完全相等(比如一个1.214,一个1.215),那么程序就会判断为C∞v 对称,那么优化结果虽然键长几乎相等,但仍然认为是C∞v ,这个从振动频率或者分子轨道对称性上可以看出来。--我们知道,CO2实际上是直线的两边对称的构型,其对称性应该是D∞v 。因此,为了得到高的对称性,必须输入的时候,精确地输入数值,比如sqrt(2),就要保留很多的小数点,180.0角,就不能写成179.9。

2。有时计算过程中对称性会变化,比如做过渡态的时候,这时需要用 IOP(2/16=3),否则计算会出错退出。

3。比如用直角坐标输入一个正三角形构型,其对称性应该是D3h,但是如果输入的小数点后面的数字不够多,那么常常得到的是C2v或其它。为了消除输入文件中坐标的有效位数的影响,得到较高的对称性,可以降低对称性判断的严格性。一般可以用symm=loose,这等价于 IOP (2/17=4, 2/18=3)。还可以减小这4和3这两个数值,使得更加loose,但不能过小,否则会出错。symm=loose只是在第一步判断输入构型的对称性时用到。
此外,也可以用gaussview来调整设置初始构型的对称性。

4。如果要降低对称性,那么可以用symm(PG=C3v)等等来做。使判断出来的对称性为C3v的一个子群。即由PG来限制最高对称性。

附用到的IOP的详细解释。

IOp(2/16)
action taken if the point group changes during an optimization. IOP(2/16=0)Abort the job. IOP(2/16=1)Keep going. IOP(2/16=2)Keep going and leave symmetry on, using the old symmetry. IOP(2/16=3)Keep going and leave symmetry on, using the new symmetry

IOp(17)
Tolerance for distance comparisons in symmetry determination.
0  Default (determined in the symmetry package, currently 1.d-8).
N>0  10-N.
N<0 10N, use the same tolerance for orientation.

IOp(18)
Tolerance for non-distance comparisons in symmetry determination.
0  Default (determined in the symmetry package, currently 1.d-7).
N>0  10 -N.
N<0  10 N, use the same tolerance for orientation.

(supi)
附加一个我的:有时候分子对称性太高经常无法优化收敛,这时可以加一个nosymm,有时挺管用的。

关于频率计算 关键字 Freq

频率计算,需要求能量对坐标的二阶导数,得到力常数,然后除以原子质量,求得振动频率。得到的振动频率,对应于实验上的红外和拉曼光谱(根据对称性,判断是否有红外或拉曼活性)。默认计算红外振动的强度,不计算拉曼强度。

频率计算必须在优化好的构型下进行。因此,必须使用和优化一样的方法和基组,进行频率计算。

有些方法,使用解析方法求解二次导数,例如DFT,HF等;而有些方法只能用数值方法,例如CCSD等。数值方法计算量很大,但是如果中断了,可以用freq(restart)接着计算,而解析方法则不能restart。

频率修正因子,是一个经验的东西,别人也是计算了一系列分子以后,总结(或者说拟合)出来的。所以有多准,是不是适合你在算的体系,就不敢保证了。
一些修正因子:
http://srdata.nist.gov/cccbdb/vsf.asp
http://www.bsc.ustc.edu.cn/~dxl/gaussian/VFS.mht (从前一个地方下载回来,便于国内访问)
关于虚频
(原贴:http://210.34.15.126/cgi-bin/top ... amp;show=0&man=

首先,什么是频率。
中学的时候我们学过简谐振动,对应的回复力是f=-kx,对应的能量曲线,是一个开口向上的二次函数E=kx^2/2. 这样的振动,对应的x=0的点是能量极小值点(简单情况下也就是最小值点)。这时的振动频率我们也会求:ω=2π sqrt(k/m)。显然它是一个正的频率,也就是通常意义下的振动频率。

那么,一维情况下,如果能量曲线是一个开口向下的二次曲线呢?首先,从能量上看,这是个不稳定的点,中学的物理书上称为“不稳平衡”。用现在的观点看,就是这一点导数是零(受力为0),且是能量极大值。如果套用上面的公式,“回复力”f=-k'x(实际上已经不是回复,而是让x越来越远了),这里k'是个负数,ω=2π sqrt(k'/m)显然就是一个虚数了,即所谓的虚频。Gaussian里面给出一个负的频率,就是对应这个虚频的。

实际情况下,分子的能量是一个高维的势能面,构型优化的时候,有时得到了极小值点,这样这个点的任意方向上,都可以近似为开口向上的二次函数,这样这里对应的振动频率就都是正的。对于极大值点,在每个方向都是开口向下的二次函数,那么频率就会都是负的——当然一般优化很少会遇到这样的情况。对于频率有正有负的情况,说明找到的点在某些方向上是极大值,有些方向上是极小值。如果要得到稳定的能量最低构型,显然需要通过微调分子的构型,消去所有的虚频。如何微调?要看虚频的振动方向。想象着虚频对应的就是开口向下的二次函数,显然,把分子坐标按照振动的方向移动一点点,分子应该就可以顺着势能面找到新的稳定点,但是也不能太小。而所谓的过渡态,则是连接反应物和产物之间的最低能量路径上的能量极大值。好比山谷中的A,B两点,它们之间的一个小土丘,就是过渡态,从A到B的反应,需要越过的是这个小土丘,而不是两边的高山。这样,过渡态就是在一个方向上是极大值,而在其它方向上都是极小值的点。因此,过渡态只有一个虚频。

优化得到虚频,消虚频的方法,就是根据虚频对应的振动位移(输出文件每一个频率下面都有一个类似坐标的3列数值,给出的就是每个原子在这个振动频率时候的移动方式),把这些位移,乘以一个适当大小的因子,直接和原来平衡构型相加,就得到新的构型了。

举例说明:

#HF/3-21G opt freq

test

0 2
H 0.0 0.0  0.0
H 0.0 0.0  1.0
H 0.0 0.0 -1.0

上面的计算,是优化一个 H-H-H的直线构型然后算频率。我们知道,3个H的稳定构型,应该是一个H2分子,加上一个H原子。但是由于我们输入的时候,中间H原子两边的H的键长相等,因此,在这个对称性的限制下,结果给出的“稳定构型”是:
   H           0.000000    0.000000    0.000000
   H           0.000000    0.000000    0.934091
   H           0.000000    0.000000   -0.934091
频率分析,有一个虚频:

                    1                      2                      3
                   SGU                    PIU                    PIU
Frequencies -- -2291.1909              1121.1486              1121.1486
Red. masses --     1.0078                 1.0078                 1.0078
Frc consts  --     3.1172                 0.7464                 0.7464
IR Inten    --    42.9706                 8.8167                 8.8167
Raman Activ --     0.0000                 0.0000                 0.0000
Depolar (P) --     0.0000                 0.0000                 0.0000
Depolar (U) --     0.0000                 0.0000                 0.0000
Atom AN      X      Y      Z        X      Y      Z        X      Y      Z
  1   1     0.00   0.00   0.82     0.82   0.03   0.00    -0.03   0.82   0.00
  2   1     0.00   0.00  -0.41    -0.41  -0.02   0.00     0.02  -0.41   0.00
  3   1     0.00   0.00  -0.41    -0.41  -0.02   0.00     0.02  -0.41   0.00

第一个是虚频。可以看出,其振动模式是:中间的H的z坐标变大,两边的两个H的z坐标变小,(然后是中间的H的z坐标变小,两边的H的z坐标变大,完成一次振动),这个如果还不清楚,画一下图就知道了。

现在,我们就把原来的坐标,加上振动模式对应的坐标:

   H           0.000000    0.000000    0.000000
   H           0.000000    0.000000    0.934091  +
   H           0.000000    0.000000   -0.934091

   0.00   0.00   0.82   
   0.00   0.00  -0.41      0.00   0.00  -0.41   直接加的结果,是
   H           0.000000    0.000000    0.820000
   H           0.000000    0.000000    0.524091
   H           0.000000    0.000000   -1.344091

画图知道,显然对原来的构型,变化太多了。对这个例子,如果取一个系数0.1乘以振动坐标,再求和,那么结果就是:

   H           0.000000    0.000000    0.082000
   H           0.000000    0.000000    0.893091
   H           0.000000    0.000000   -0.975091

显然就合理多了。用这个坐标去重新优化,就可以得到真正的没有虚频的稳定结构了。
这个系数0.1,只是个经验的数值,对于不同体系,可以自己设定。设置过小,会得到同样的虚频,设置过大,可能会得到别的构型(当然也可能碰巧得到更好的构型)。
  



关于BSSE

设计算AB体系的BSSE矫正,过程如下:
1.首先正常计算AB体系的能量E(AB)。
2.计算在AB基组下,A的能量E(A/AB)。
3.计算在AB基组下,B的能量E(B/AB)。
4.正常计算A的能量E(A)。
5.正常计算B的能量E(B)。
这样,可以找到5个scf计算,和相应的5个能量。

然后程序会自动给出BSSE修正能Counterpoise: BSSE energy =E(A)+E(B)-E(A/AB)-E(B/AB),以及经过BSSE修正以后体系的总能量Counterpoise: corrected energy=E(AB)+BSSE energy  (记为E6).

一般来说,基组变大的结果是能量更低,因此,一般BSSE energy>0,修正后的体系能量比原来稍高一点。

不考虑BSSE,我们计算的结合能是用的 Eb=E(AB)-E(A0)-E(B0),注意,这里A0和B0分别是A,B的各自的优化得到的稳定构型的能量。那么考虑了BSSE以后,就是(Eb)'= Eb+BSSE energy.
当然也就是 E6- E(A0)-E(B0).
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 wqd198686 的主题更新
普通表情 高级回复 (可上传附件)
信息提示
请填处理意见