为什么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)还好？这下为难了。。。

诸葛壹次心

gaussian处理unrestricted的奇怪结果讨论

sobereva

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

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

jjjspring

关注一下UHF的S^2值对于键长的变化关系，超过不稳定点之后应该是由0逐渐增加的。

诸葛壹次心

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相比哪个更准确?

诸葛壹次心

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逐渐增加的”的经验判断出计算不对？

诸葛壹次心

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

诸葛壹次心

还有一个一直困扰了我很久的这个体系的multiplicity究竟怎么取的问题，感觉总算是明白了一些，在键长到达一定长度之后原来计算结果是一样的！

诸葛壹次心

自旋极化单重态计算继续探讨

jjjspring

诸葛壹次心 发表于 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）都包括，可以参考。

sobereva

写guess(mix,always)再试，这才能确保每一步都重新以mix方式生成初猜。

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

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

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

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

诸葛壹次心

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

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

诸葛壹次心

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

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

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

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

谢谢sob！

诸葛壹次心

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输入文件如下：

1. ***,CH4
   
2. memory,3000,M
   
3. print,basis,orbitals
   
4. angstrom
   
5. geometry={                !define the nuclear coordinates
   
6. C;
   
7. H, 1, R;
   
8. H, 1, 1.08600001, 2, 109.47120259;
   
9. H, 1, 1.08600000, 2, 109.47120275, 3, 120.00001386;
   
10. H, 1, 1.08600000, 2, 109.47120275, 4, 119.99997229;
    
11. }
    
12. basis=vdz                  !triple-zeta basis set
    
13. R=0.8
    
14. do i=1,43,1
    
15. R_total(i)=R
    
16. uhf;wf,spin=0;rotate,5.1,6.1,45
    
17. {casscf
    
18. closed,4
    
19. occ,6
    
20. }
    
21. escf(i)=energy
    
22. R=R+0.1
    
23. enddo
    
24. {table,R_total,escf
    
25. head,R_total,escf
    
26. save,CH4_29_mol_cas22_opt_080.tab,new
    
27. title,Results for CH4_29_mol_cas22_opt_080, basis $basis
    
28. sort,1,2}
    
29. ---
    

复制代码

我试图使用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的结果作为初猜就好啦~如下

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

复制代码

诸葛壹次心

双自由基计算进一步讨论