为什么orca和gaussian的T1诊断不一致？

zealot_cdp86

大家好，请教一下T1诊断的问题

近期分别计算了一下萘(C10H8)分子和萘分子在氢提取反应中过渡态（C10H9）的T1。注后面用R2代表萘，R2-H代表萘分子在氢提取反应中过渡态，具体结构见附件。


通过orca（4.1.2）和gaussian（g16）分别计算R2-H的CCSD/cc-pvdz，结果发现orca的T1=0.0145，但是g16的T1=0.018。同时能量也有0.16kcal/mol的差别（g16-orca）。

正好手边上有R2的结果，不过是CCSD/cc-pvtz，结果发现orca的T1=0.010549，但是g16的T1=0.010547，两者非常接近，CCSD能量也很接近(<0.0005kcal/mol)。

请教一下，为什么orca和g16计算的R2-H的T1不一致？为什么单点能也有较小的误差？

谢谢！

Warm_Cloud

这明明是一致了。。。。两个软件收敛的标准不一样

sobereva

如置顶的新社员必读贴和论坛首页的公告栏所示，求助帖必须在帖子标题明确体现出“提问”、“求助”要素（仔细看http://bbs.keinsci.com/thread-9348-1-1.html），我已把你的不恰当标题“orca和gaussian的T1诊断不一致”改了，以后注意

zealot_cdp86

sobereva 发表于 2020-2-9 05:53
如置顶的新社员必读贴和论坛首页的公告栏所示，求助帖必须在帖子标题明确体现出“提问”、“求助”要素（仔 ...

谢谢！

zealot_cdp86

Warm_Cloud 发表于 2020-2-8 23:00
这明明是一致了。。。。两个软件收敛的标准不一样

R2-H单点能的差别可以接受，但是T1这个不是太确定。手上还有一个苯在氢提取反应的过渡态的结果，g16/orca/molpro的T1都不相同，差别还很大，具体如下：

g16, T1=0.045 (注：这个用了xqc,感觉得到的wavefunction应该是有问题的，能量比orca和molpro大~1.4 kcal/mol)
orca, T1=0.026 (注：单点能orca和molpro的差别为0.28 kcal/mol， orca结果大一点)
molpro, T1=0.014

如果换做这个例子，先只看orca和molpro的结果，orca的T1和molpro的T1定性不一致（T1建议标准0.02），那么以是以molpro的结果为准？还是orca?假设是收敛标准的问题，那么推广一点的说，什么样的收敛标准给出的T1的定性分析是正确的？

另外，谢谢提醒，目前在测试一下收敛标准的影响

sobereva

0.010549和0.010547是完全一样的，而0.045和0.026这才叫不同。
先确保相同基组下HF的单点能基本一样（比如Hartree为单位相符到小数点后五位），之后再说T1的事。如果都没收敛到相同波函数上，自然T1明显不同。
没有什么以谁为准之说，三个程序给出T1的算法都一样，正确操作情况下可以得到相同结果。

mutron

ORCA和gaussian的HF能量应该是精确一致（哪怕conver=6、normalscf对应1.0e-6 Eh，都小于0.001kcal/mol）

E(UHF)= -383.839248914704 Eh   <S**2>= 0.844184   T1 diagnostic= 0.017961059   ORCA tightscf（1.0e-8）
E(UHF)= -383.839600004203 Eh   <S**2>= 0.807829   T1 diagnostic= 0.014535132   ORCA normalscf（default、1.0e-6）
E(UHF)= -383.839254149      Eh   <S**2>= 0.8445       T1 diagnostic= 0.01799514     G16 scf=tihgt（default、1.0e-8）

收敛限都是1.0e-8的时候两者基本没有区别，问题出在ORCA默认SCF收敛限不太高的时候，有可能获取的是并没有收敛的波函数。把normalscf ORCA任务moread再做一次单点，还继续迭代了22次：E(UHF)= -383.839295799078 Eh   <S**2>=0.848578，  T1 diagnostic= 0.018347483 这时候结果就非常接近了。再重复一次只迭代了一遍。R2没有这个问题可能是闭壳层不容易出现这种情况

另外，G16、ORCA默认得出的这两个都不是稳定的波函数，需要stable=opt、STABPerform true STABRestartUHFifUnstable true获取稳定的波函数，再做几何优化。真正基态波函数对收敛限就不太敏感了：
E(UHF)= -383.859526509753 Eh   <S**2>= 1.956226   ORCA tightscf（1.0e-8）
E(UHF)= -383.859553794952 Eh   <S**2>= 1.956376   ORCA normalscf（default、1.0e-6）
E(UHF)= -383.859526503      Eh   <S**2>= 1.9562       G16 scf=tihgt（default、1.0e-8）

碰到一个问题想问问sober老师，G16的chk文件转为gbw作为ORCA初猜只迭代了一次就收敛了（必须int=NoBasisTransform）；但是反过来ORCA gbw转为molden-fchk-chk，G16却需要迭代二十多次，从Initial guess from the checkpoint file <S**2>= 0.7500看也不对，这是哪里出了问题？（R2闭壳层反应物没有这个问题）

zealot_cdp86

mutron 发表于 2020-2-10 22:37
ORCA和gaussian的HF能量应该是精确一致（哪怕conver=6、normalscf对应1.0e-6 Eh，都小于0.001kcal/mol）

...

特别感谢！！今天也做了点功课，得出的结论和你一致，是收敛标准的问题。测试结果如附图。另外所有结果可以从下面的link看到。
https://www.jianguoyun.com/p/DRxkwucQ-NToBhiHv9cC

算是上了一课，orca的默认收敛标准可能会导致得到的波函数不正确。

注：molpro那没软件，就没更新，不过，R1-H的molpro结果不对，猜想也是默认收敛标准的问题，导致得到的波函数不正确。

sobereva

mutron 发表于 2020-2-10 22:37
ORCA和gaussian的HF能量应该是精确一致（哪怕conver=6、normalscf对应1.0e-6 Eh，都小于0.001kcal/mol）

...

三重态那个问题，原因在于虽然Multiwfn导出的fch对应非限制性波函数，但unfchk的时候beta轨道被丢掉了，导致产生的chk里面只有alpha部分，因此相当于用了ORCA的alpha轨道同时作为Gaussian初始的alpha和beta轨道（这点可以通过guess(read,only)输出的初猜轨道系数来验证）。
至于怎么让unfchk将Multiwfn导出的fch正确识别成非限制性的从而把beta轨道在转换成chk时也保留下来目前我还不清楚，应该是fch文件里一些杂七杂八的信息会影响unfchk的判断。

2021-Oct-1注：以上问题在2021-Oct-1及以后更新的Multiwfn版本中已经解决，具体情况见本帖22楼

目前可以用Multiwfn的主功能100的子功能2里的保存gjf文件的功能，程序会问是否把当前的波函数直接写入gjf里面，这样可以令非限制性波函数的alpha和beta轨道都给Gaussian作为初猜。缺点是这样的gjf文件对于大体系、用高质量基组时会较大。

mutron

sobereva 发表于 2020-2-13 18:48
三重态那个问题，原因在于虽然Multiwfn导出的fch对应非限制性波函数，但unfchk的时候beta轨道被丢掉了， ...

谢谢老师的解答

guess=(read,only)读取gbw转为chk文件的轨道系数的确alpha和beta居然是一样的（不过和ORCA的alpha和beta都不一样），molden转为fchk里面还有beta轨道，unfchk得到的chk再转为fchk里面的确就剩alpha轨道了。对比了高斯任务chk转fchk后发现，gbw转fchk里面缺了有关ILSW的那一段导致unfchk的时候没有了beta轨道（想问一下ILSW是什么信息）

用gjf好像会出错，读取不了轨道。由multiwfn转出来的gjf里面第一行指定格式的是5(E16.9)，我看到高斯手册里面的格式是(*E**.*)，改成(5E16.9)可以顺利读取迭代4次收敛，<S**2>也非常一致，不过我看初猜（guess=(read,only)）的轨道能量和ORCA不太一样，是不是这样改会有点问题？

sobereva

mutron 发表于 2020-2-15 02:11
谢谢老师的解答

guess=(read,only)读取gbw转为chk文件的轨道系数的确alpha和beta居然是一样的（不过和 ...

ILSW是Gaussian内部的信息，我不晓得unfchk是怎么从中判断的，暂时也没时间研究。
5(E16.9)用windows版能正常读取，linux下失败。我已更新官网上的Multiwfn改成了(5E16.9)。
初始能量和ORCA不同应当是积分格点不同的问题，也可能是Gaussian在初期自动先降低积分精度之类导致的

mutron

sobereva 发表于 2020-2-16 09:43
ILSW是Gaussian内部的信息，我不晓得unfchk是怎么从中判断的，暂时也没时间研究。
5(E16.9)用windows版 ...

谢谢老师，新版没问题了

zealot_cdp86

mutron 发表于 2020-2-10 22:37
ORCA和gaussian的HF能量应该是精确一致（哪怕conver=6、normalscf对应1.0e-6 Eh，都小于0.001kcal/mol）

...

您好，再请教一下关于自旋的问题。如你计算的stable波函数结果显示，最后得到的的<S**>=1.956，这比理想自旋结果要高很多，这能说明得到的波函数不是真实的么？还是说如果能通过STABPerform true STABRestartUHFifUnstable true测试，就说明波函数是没问题的？无需管自旋问题？谢谢！

mutron

zealot_cdp86 发表于 2020-2-18 23:39
您好，再请教一下关于自旋的问题。如你计算的stable波函数结果显示，最后得到的的=1.956，这比理想自旋结 ...

<S**2>与理想值偏差较大也就是波函数和真实波函数偏差大，只要alpha、beta轨道不唯一匹配（自旋极化）就会出现偏差（自旋污染），但这并不代表基于UHF算出来的CCSD能量不准（基组得加大）；stable是判断波函数是否稳定，和自旋污染两个问题

zealot_cdp86

mutron 发表于 2020-2-19 02:17
与理想值偏差较大也就是波函数和真实波函数偏差大，只要alpha、beta轨道不唯一匹配（自旋极化）就会出现 ...

UHF波函数问题有些细节还是比较模糊想请教一下。

1. 一直以为UHF最终会得到“真实”的稳定波函数。但是对于这个讨论例子而言，<S**2>明显较大，说明与真实波函数差别较大，那为什么后面CCSD的计算还有意义么？为什么说不代表CCSD能量不准？另外，以前见过很多大个10-20%的值，这是第一次见大这么多的。
2. 对于UHF后面的CCSD（T）计算，我可不可以这样理解，UHF收敛的波函数需要是当前结构下能量最低的并且稳定的波函数，自旋污染不是主要问题。

我看sober老师的基础量化培训班讲义（page128, slide77）最下面有一句话：
“对于双杂化和后HF，如果其参考态波函数是稳定的，则最终得到的结果也相当于对应基态波函数的”

那么如果有一个体系通过不同的数值方法得到了两个不同能量的stable的UHF波函数，同时假设体系自旋布居看起来都无法说是有定性错误，那么应该取哪一个来做后面的CCSD(T）？能量最低的？还是说这两个最后给出的CCSD（T）的结果“一定”会差不多？正好有个例子正在测试中，想看看最后CCSD(T）结果差别大不大。。。

另外，基组得加大是为什么？这里pvdz只是为了测试方便。