计算化学公社

标题: 为什么GAUSSIAN中的restricted和unrestricted结果一样，而其他人却不是？ [打印本页]

作者
Author: 诸葛壹次心 时间: 2014-12-16 19:58
标题: 为什么GAUSSIAN中的restricted和unrestricted结果一样，而其他人却不是？
我算了一个小体系，算甲烷在一个C-H键变化过程中的能量变化（the potential energy curve was obtained by altering the length of a single C–H bond while keeping the remaining three C–H bonds at their equilibrium bond length ~1.086 Å and the HCH angles at the tetrahedral value，Full configuration interaction potential energy curves for breaking bonds to hydrogen:
An assessment of single-reference correlation methods，Antara Dutta and C. David Sherrill），但是发现了一个奇怪的现象。
如下图所示，无论是B3LYP/MP2/CCSD/CCSD(T)，我用gaussian算的restricted和unrestricted结果都一模一样，并且和文献中的restricted是一样的，但是他的unrestricted却有更好的结果，不知道是gaussian自己的处理有问题，还是我哪里设置得不对？贴上了我的一个log文件供检查。

另外，好像QCISD(T)的计算结果比CCSD(T)还好？这下为难了。。。

作者
Author: 诸葛壹次心 时间: 2014-12-16 20:03
gaussian处理unrestricted的奇怪结果讨论

作者
Author: sobereva 时间: 2014-12-16 20:07
你得获得对称破缺态才行。
直接算的话，无论你是用闭壳层来算，还是非限制性开壳层来算，结果和闭壳层都是一样的。必须用guess=mix或者片段初猜波函数或者stable=opt等方式获得对称破缺态（如果存在的话），这时非限制性开壳层的解和闭壳层才不同。只有键长超过一定长度后（所谓的不稳定点）用这些方法才有可能得到对称破缺态。

详见
谈谈片段组合波函数与自旋极化单重态
http://sobereva.com/82

作者
Author: jjjspring 时间: 2014-12-16 23:43
关注一下UHF的S^2值对于键长的变化关系，超过不稳定点之后应该是由0逐渐增加的。

作者
Author: 诸葛壹次心 时间: 2014-12-21 01:46

sobereva 发表于 2014-12-16 20:07
你得获得对称破缺态才行。
直接算的话，无论你是用闭壳层来算，还是非限制性开壳层来算，结果和闭壳层都是 ...

原来这篇帖子在这里用的！太神奇啦~谢谢sob！
我之前使用scan扫出来的曲线，发现加上guess=mix后没有用，猜想是不是scan的时候只对第一个键长用了而后面没有，只好用了一个笨办法，用modredundant或者内坐标限制住键长，然后加上guess=mix，每一个键长都写了一个单独的输入文件，加上一些尝试性的计算，大约算了1000个jobs，真是累死了。。。好在总算是搞定啦！
最后的图像是这样的：

UMP2/UQCISD(T)/UCCSD(T)总算是对上啦，可是UB3LYP还是和paper中的不一样，目前还没有找出原因（窃以为大概是作者放错数据啦，和他图片上画的能量差也不太一致，所以就不管啦）。但是进而我又遇到了这么一个问题：
键能应该是BDE(CH3-H)=E(CH3)+C(H)-E(CH4)，但是我看着这个图上的能量差的好像有点不太对劲。
我先查了一下BDE(CH3-H)的数据，有BDE(CH3-H, 298)=103.9 kcal/mol, BDE(CH3-H，460)=104.8 kcal/mol，我估计0K下应该也差不多是103左右，然后又用B3LYP/6-31G(D)和CCSD(T)/6-31G(D)算了一下（已加上ZPE），却发现B3LYP算出来103.19kcal/mol at 0K, 104.82kcal/mol at 298.15K, 而CCSD(T)算出来 94.80kcal/mol at 0K, 96.45 kcal/mol at 298.15K。CCSD(T)反而没有B3LYP准确了，而且误差很大，这是为什么呢？
另外一个小问题是CCSD(T)算频率好像比opt还慢的样子，这正常吗？CCSD(T)的频率计算和B3LYP相比哪个更准确?

作者
Author: 诸葛壹次心 时间: 2014-12-21 02:21
本帖最后由诸葛壹次心于 2014-12-21 02:22 编辑

jjjspring 发表于 2014-12-16 23:43
关注一下UHF的S^2值对于键长的变化关系，超过不稳定点之后应该是由0逐渐增加的。

谢谢@jjjspring !
有道理，果然增加了！
在0.8-1.6的键长的时候都是0，在1.7开始就逐渐增加了，到3左右就增加到差不多1了，见下图。可是有几个小问题想请教一下：
1. log文件中的S**2就是S(S+1)吗？
2. 它有一个before annihilation和after annihilation，那么就是说单点能计算中已经自动消除了一部分的自旋污染（我用的UCCSD(T)/6-31G(d) guess=mix）？那么after annihilation到什么程度就认为可以接受了呢？我之前听Xiaosong Li(我不知道是李晓松还是李小松啊，反正是讲gaussian课的那个。。。)老师说差不多和预计的S(S+1)的10%以内就行了，可是这里S=0怎么破……我想S=1/2的时候是0.075, S=1的时候是0.2，那S=0不得至少0.05以下？可是这个算出来的after annihlation最后涨到0.06了，岂不是说算得不对？那该怎么办呢？
3. 我也试着用UCCSD(T)/6-31G(d)不加mix直接算了一下，s**2就一直是0了，就从数值上看起来好像很正常啊，那么我们可以从这个参数上发现计算存在问题么？难道是基于“超过不稳定点之后应该是由0逐渐增加的”的经验判断出计算不对？

作者
Author: 诸葛壹次心 时间: 2014-12-21 02:27
另外我其实还想问这样一个问题：为什么CCSD(T)算出来能量比UCCSD(T)更低但我们却偏要用unrestricted的计算结果？
我试着自己找了一下原因，对每一个键长做了下stable的计算，发现好像从s**2不等于0的时候就有instability了，所以是波函数不稳定所以结果不能采用么？

作者
Author: 诸葛壹次心 时间: 2014-12-21 02:32
还有一个一直困扰了我很久的这个体系的multiplicity究竟怎么取的问题，感觉总算是明白了一些，在键长到达一定长度之后原来计算结果是一样的！

作者
Author: 诸葛壹次心 时间: 2014-12-21 02:36
自旋极化单重态计算继续探讨

作者
Author: jjjspring 时间: 2014-12-21 03:21

诸葛壹次心发表于 2014-12-21 02:36
自旋极化单重态计算继续探讨

你这个体系很像双自由基，推荐你看下sob老师写的 "CASSCF计算双自由基以及双自由基特征的计算"里面UHF部分。有空你再回顾下RHF，UHF，和CI 对 H2 = H+H 的处理。
http://bbs.keinsci.com/forum.php?mod=viewthread&tid=322
在平衡位置或者靠近平衡位置，一般认为CCSD（T）优于B3LYP的，远离平衡的位置，这两个方法都不那么好。多组态方法，耦合簇方法包含到T甚至Q会更好。高斯没有CCSD（T）的解析梯度和解析力常数，用CCSD(T)算频率巨慢。
还有个建议，你算的这些小分子CH3， CH4的能量，频率之类的，结果自己拿不准对不对，可以上nist quantum chemistry database上比较一下，那里有常见分子不同方法不同基组下的结果，像B3LYP，CCSD（T）都包括，可以参考。

作者
Author: sobereva 时间: 2014-12-21 04:15
写guess(mix,always)再试，这才能确保每一步都重新以mix方式生成初猜。

CCSD(T)用那么小的基组没什么价值，此时不如B3LYP是很可能的。

CCSD(T)没解析梯度，优化本来就很慢了，频率要算二阶导数显然更慢。因此不建议用CCSD(T)优化和计算频率。

做谐振计算，绝不代表越高级别的方法算出来的频率就越准确，频率校正因子起到了关键性影响，见《谈谈谐振频率校正因子》（http://sobereva.com/221）。CCSD(T)算出来的频率不考虑频率校正因子的话还不如用HF在考虑的时候更准确。这也是为什么追求精度的热力学组合方法甚至用HF/6-31G*来计算频率。

平衡状态下是什么自旋多重度，scan计算的时候还是什么样。拉远之后，变成两个独立的碎片，它们自旋相同和自旋相反摆放能量一致，因此就简并了。

作者
Author: 诸葛壹次心 时间: 2014-12-21 14:59

jjjspring 发表于 2014-12-21 03:21
你这个体系很像双自由基，推荐你看下sob老师写的 "CASSCF计算双自由基以及双自由基特征的计算"里面UHF部 ...

对啊，其实我主要想研究的是barrierless的这一类反应。那么什么样的反应具有barrierless的特征呢？我一直以为barrierless的都是这种双自由基的反应，或者是R+O2这样的multiplicity-2的反应，但上次听一个同学说H+烯烃也是，感觉好奇怪啊，这方面有什么规律么？
sob的那篇帖子其实之前就看过，只不过当时好多地方没有看懂。。。感觉好惭愧。。。这次又看了以便果然收获多多~
CCCBDB的网站也很有用，果断收藏之~
谢谢jjj@jjjspring ！

作者
Author: 诸葛壹次心 时间: 2014-12-21 15:38

sobereva 发表于 2014-12-21 04:15
写guess(mix,always)再试，这才能确保每一步都重新以mix方式生成初猜。

CCSD(T)用那么小的基组没什么价 ...

啊，原来还有一个这么好用的关键词！果然科学技术才是第一生产力！好在我之前也是写了小脚本批量生成，不然就亏大发了。。。

果然CCSD(T)换上了cc-pvtz一下子键能就靠谱啦！

谢谢sob！

作者
Author: 诸葛壹次心 时间: 2015-2-1 15:51
本帖最后由诸葛壹次心于 2015-2-3 02:55 编辑

sobereva 发表于 2014-12-21 04:15
写guess(mix,always)再试，这才能确保每一步都重新以mix方式生成初猜。

CCSD(T)用那么小的基组没什么价 ...

我继续使用molpro计算这个体系（最终目标是学会用caspt2计算这种双自由基的体系啦~），但是有个奇怪的现象。我试图按照http://bbs.keinsci.com/forum.php?mod=viewthread&tid=322中那样使用
#P UHF/6-31G* guess=mix
#P CAS(2,2,UNO)/6-31G* guess=read
但是molpro中好像不能直接mix，我就按网上说的改成了rotate,5.1,6.1,45，但是好像不能得到和gaussian一样的结果，不知道为什么？
贴上我的molpro输入文件如下：

***,CH4
memory,3000,M
print,basis,orbitals
angstrom
geometry={ !define the nuclear coordinates
C;
H, 1, R;
H, 1, 1.08600001, 2, 109.47120259;
H, 1, 1.08600000, 2, 109.47120275, 3, 120.00001386;
H, 1, 1.08600000, 2, 109.47120275, 4, 119.99997229;
}
basis=vdz !triple-zeta basis set
R=0.8
do i=1,43,1
R_total(i)=R
uhf;wf,spin=0;rotate,5.1,6.1,45
{casscf
closed,4
occ,6
}
escf(i)=energy
R=R+0.1
enddo
{table,R_total,escf
head,R_total,escf
save,CH4_29_mol_cas22_opt_080.tab,new
title,Results for CH4_29_mol_cas22_opt_080, basis $basis
sort,1,2}
---

复制代码

我试图使用uhf;wf,spin=0;rotate,5.1,6.1,45达到#P UHF/6-31G* guess=mix的效果，然后用
{casscf
closed,4
occ,6}来读取和计算，可是好像结果没有变化。（还没有学会怎么在molpro中用cas(2,2,uno)，不过好像加不加uno结果应该差不多？）
以下是结果

----------------------------
搞了一整天终于找到原因啦，没人回答我就来自问自答一下啦~
主要是没有加上想要的初始猜测，molpro中casscf默认寻找初猜的优先级是mcscf，然后scf
所以按照上面的写法就相当于scan的时候只在第一步使用了scf的结果，后面都是用的casscf自己。
所以加上了start,2200.2调用uhf的结果作为初猜就好啦~如下

***,CH4_4
memory,3000,M
print,basis,orbitals
angstrom
geometry={ !define the nuclear coordinates
C;
H, 1, R;
H, 1, 1.08600001, 2, 109.47120259;
H, 1, 1.08600000, 2, 109.47120275, 3, 120.00001386;
H, 1, 1.08600000, 2, 109.47120275, 4, 119.99997229;
}
basis=vdz !triple-zeta basis set
R=0.8
do i=1,43,1
R_total(i)=R
uhf;wf,spin=0;rotate,5.1,6.1,45
{casscf
closed,4
occ,6
wf,10,1,0
start,2200.2 !感觉在uhf结束后强制保存到2200.2会更保险，但是这个case结果都一样就没加啦~
}
escf(i)=energy
R=R+0.1
enddo
{table,R_total,escf
head,R_total,escf
save,CH4_4.tab,new
title,Results for CH4_4, basis $basis
sort,1,2}
---

复制代码

作者
Author: 诸葛壹次心 时间: 2015-2-1 15:54
双自由基计算进一步讨论

作者
Author: helpme 时间: 2015-2-1 15:56
我想问一句，换了基组以后，多重度3的和 Restrict多重度1的，到底哪个能量低（在键长比较大的地方）？

作者
Author: 诸葛壹次心 时间: 2015-2-1 16:18

helpme 发表于 2015-2-1 15:56
我想问一句，换了基组以后，多重度3的和 Restrict多重度1的，到底哪个能量低（在键长比较大的地方）？

你是指换成什么基组啊？
8#的图是在6-31g(d)下算的，这时Restricted多重度1的低于多重度3的（键长较大时）。你是指换成cc-pvtz么？

作者
Author: helpme 时间: 2015-2-2 09:38

诸葛壹次心发表于 2015-2-1 16:18
你是指换成什么基组啊？
8#的图是在6-31g(d)下算的，这时Restricted多重度1的低于多重度3的（键长较大时 ...

是的，就是想知道换成cc-pvtz时，到底RCCSD(T)和UCCSD(T)哪个低？
因为看你说 “果然CCSD(T)换上了cc-pvtz一下子键能就靠谱啦！”

作者
Author: 诸葛壹次心 时间: 2015-2-3 02:48

helpme 发表于 2015-2-2 09:38
是的，就是想知道换成cc-pvtz时，到底RCCSD(T)和UCCSD(T)哪个低？
因为看你说 “果然CCSD(T)换上了cc-pv ...

换成cc-pvtz之后见下图，三重态和单重态（单重态中uccsd(t)，guess=mix）在键长较长时是一样的能量。原因还是sob之前说的。但是RCCSD(t)稍微变了一点，之前是偏低，现在是偏高。

作者
Author: helpme 时间: 2015-2-3 19:34
“现在是偏高”-------我想要的就是这句话。

作者
Author: 诸葛壹次心 时间: 2015-2-4 14:22

helpme 发表于 2015-2-3 19:34
“现在是偏高”-------我想要的就是这句话。

为什么想要这个呢？能详细说明一下原因么？谢谢~

作者
Author: sobereva 时间: 2015-2-4 15:03

诸葛壹次心发表于 2015-2-4 14:22
为什么想要这个呢？能详细说明一下原因么？谢谢~

RHF不能正确表现解离曲线，拉远后能量会比两个片段之和要高，见此贴讨论
小谈几种量子化学方法对氢分子解离过程的描述
http://hi.baidu.com/sobereva/item/fa2772ad3197ac17a9cfb7d8
相应地，由于RHF波函数作为参考态不正确，RCCSD(T)也没法正确描述解离曲线。

作者
Author: 诸葛壹次心 时间: 2015-2-4 17:13

sobereva 发表于 2015-2-4 15:03
RHF不能正确表现解离曲线，拉远后能量会比两个片段之和要高，见此贴讨论
小谈几种量子化学方法对氢分子 ...

经典啊！原来是这样！谢谢sob~
那么这篇文中提到的“势能面交叉”如果存在的话，是否就需要用态平均（state average）和多态（multi state）的方法来解决呢？（我前几天刚看到这个，但是还没搞太明白）另外一般什么体系下可能会出现势能面交叉呢？

作者
Author: sobereva 时间: 2015-2-4 17:47

诸葛壹次心发表于 2015-2-4 17:13
经典啊！原来是这样！谢谢sob~
那么这篇文中提到的“势能面交叉”如果存在的话，是否就需要用态平均（st ...

如果是自旋多重度不同的态的交叉，DFT实际上就可以。
对于自旋多重度相同的态的交叉（一般是避免交叉），离交叉点比较远的地方，用TDDFT之类普通研究激发态的方法都可以研究，但离交叉点近的区域，态平均CASSCF是标准方法。如果在CASSCF基础上进一步用二阶微扰，如CASPT2，就应该用其多态版本，单态版本只适合研究某个态。

势能面交叉是很普遍的现象，特别是对有机分子来说。随便找找光化学的资料都会有很多例子。

作者
Author: 诸葛壹次心 时间: 2015-2-4 18:43

sobereva 发表于 2015-2-4 17:47
如果是自旋多重度不同的态的交叉，DFT实际上就可以。
对于自旋多重度相同的态的交叉（一般是避免交叉 ...

好的，谢谢sob~我赶紧去看看相关的例子~

PS: 看到“很普遍”这几个字顿时有种欲哭无泪的感觉，这简直就是一个无底洞啊。。。我最开始只是想研究下barrierless反应的。。。

作者
Author: sobereva 时间: 2015-2-4 19:13

诸葛壹次心发表于 2015-2-4 18:43
好的，谢谢sob~我赶紧去看看相关的例子~

PS: 看到“很普遍”这几个字顿时有种欲哭无泪的感觉，这简直 ...

之所以有著名的kasha rule，正因为激发态能量相近，交叉频繁，所以...太普遍了

PS：最近恰有一篇用了Multiwfn的研究含铂化合物的单-三重态交叉的文章有兴趣可以看看：http://www.sciencedirect.com/science/article/pii/S1566119915000312

作者
Author: 诸葛壹次心 时间: 2015-2-5 01:49

sobereva 发表于 2015-2-4 19:13
之所以有著名的kasha rule，正因为激发态能量相近，交叉频繁，所以...太普遍了

PS：最近恰有一篇用 ...

好的，谢谢sob~

作者
Author: helpme 时间: 2015-2-5 09:30

诸葛壹次心发表于 2015-2-4 14:22
为什么想要这个呢？能详细说明一下原因么？谢谢~

能量是计算化学里面的“黄金规则”，绝大多数情况下，以能量低的为准，能量低的为基态。

作者
Author: 诸葛壹次心 时间: 2015-2-5 13:29

helpme 发表于 2015-2-5 09:30
能量是计算化学里面的“黄金规则”，绝大多数情况下，以能量低的为准，能量低的为基态。

这样子啊，我之前算的那个RCCSD(T)/6-31g(d)的case检验了stability，发现是不稳定的，我想这时虽然能量更低应该也不能用作基态是吗？

作者
Author: sobereva 时间: 2015-2-5 13:46

诸葛壹次心发表于 2015-2-5 13:29
这样子啊，我之前算的那个RCCSD(T)/6-31g(d)的case检验了stability，发现是不稳定的，我想这时虽然能量更 ...

不稳定的话，用stable=opt，通常会产生出真正基态波函数，给出的能量应当比先前的更低。

作者
Author: 诸葛壹次心 时间: 2015-2-5 16:42

sobereva 发表于 2015-2-5 13:46
不稳定的话，用stable=opt，通常会产生出真正基态波函数，给出的能量应当比先前的更低。

这。。。更低的话岂不是偏得更远了？看来ccsd(t)/6-31g(d)这样的搭配果然很有问题呀。但我想起来好像CBS-QB3的方法中有一步就是CCSD(T)/6-31+G(d')，貌似两个基组也差不了多少，是不是说明键长较长时使用CBS-QB3的方法也要比较小心啊？

作者
Author: sobereva 时间: 2015-2-5 16:53

诸葛壹次心发表于 2015-2-5 16:42
这。。。更低的话岂不是偏得更远了？看来ccsd(t)/6-31g(d)这样的搭配果然很有问题呀。但我想起来好像CBS- ...

我不清楚你说的是什么偏离更远。

CBS-QB3里CCSD(T)/6-31+G(d')不是直接作为能量的，而只是用来获得高阶的相关能的校正。这和经常用的一种近似有一定类似：
CCSD(T)/大基组 ≈ CCSD(T)/小基组 - MP2/小基组 + MP2/大基组

作者
Author: 诸葛壹次心 时间: 2015-2-5 18:39

sobereva 发表于 2015-2-5 16:53
我不清楚你说的是什么偏离更远。

CBS-QB3里CCSD(T)/6-31+G(d')不是直接作为能量的，而只是用来获得高 ...

这样啊，好的，谢谢sob！

PS: 我指的偏离更远是因为之前在CCSD(T)/6-31G(d)下RCCSD(T)比UCCSD(T)的能量偏低，而按照“一般能量更低的作为基态或者稳定态”的说法，事实上应该偏高才对，但如果使用stable=opt之后RCCSD(T)的结果更低，就会与偏高的事实相去更远。

作者
Author: sobereva 时间: 2015-2-5 19:25
那你看看RHF/6-31G*和对称破缺的UHF/6-31G*谁低。理应肯定是后者低，这样的话就等于RCCSD(T)的相关能比UCCSD(T)相关能更大才导致了这种情况。

欢迎光临计算化学公社 (http://bbs.keinsci.com/)