计算化学公社

标题: 为什么orca和gaussian的T1诊断不一致？ [打印本页]

作者
Author: zealot_cdp86 时间: 2020-2-8 21:10
标题: 为什么orca和gaussian的T1诊断不一致？

大家好，请教一下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不一致？为什么单点能也有较小的误差？

谢谢！

作者
Author: Warm_Cloud 时间: 2020-2-8 23:00
这明明是一致了。。。。两个软件收敛的标准不一样

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

作者
Author: zealot_cdp86 时间: 2020-2-9 10:43

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

谢谢！

作者
Author: zealot_cdp86 时间: 2020-2-9 10:55

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的定性分析是正确的？

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

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

作者
Author: mutron 时间: 2020-2-10 22:37
本帖最后由 mutron 于 2020-2-11 01:43 编辑

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闭壳层反应物没有这个问题）

作者
Author: zealot_cdp86 时间: 2020-2-10 23:34

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结果不对，猜想也是默认收敛标准的问题，导致得到的波函数不正确。

作者
Author: sobereva 时间: 2020-2-13 18:48

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文件对于大体系、用高质量基组时会较大。

作者
Author: mutron 时间: 2020-2-15 02:11
本帖最后由 mutron 于 2020-2-15 06:57 编辑

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不太一样，是不是这样改会有点问题？

作者
Author: sobereva 时间: 2020-2-16 09:43

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在初期自动先降低积分精度之类导致的

作者
Author: mutron 时间: 2020-2-16 21:33

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

谢谢老师，新版没问题了

作者
Author: zealot_cdp86 时间: 2020-2-18 23:39

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测试，就说明波函数是没问题的？无需管自旋问题？谢谢！

作者
Author: mutron 时间: 2020-2-19 02:17

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

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

作者
Author: zealot_cdp86 时间: 2020-2-19 10:49

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只是为了测试方便。

作者
Author: zjxitcc 时间: 2020-2-19 10:59
本帖最后由 zjxitcc 于 2020-2-19 11:01 编辑

zealot_cdp86 发表于 2020-2-19 10:49
UHF波函数问题有些细节还是比较模糊想请教一下。

1. 一直以为UHF最终会得到“真实”的稳定波函数。但 ...

你问了个不错的问题。UHF有两个（或更多）稳定解时，选哪个解做CC其实是有争议的问题，没有定论。如果是新手，建议就取最低的UHF解做CC计算。如果是懂电子结构、能分辨出两个解有啥不同、哪个更有实际意义的话，可以选那个解，但是这仅限于你懂两个解有啥物理意义/化学图像。
不过，这时候实际上更明智的做法是脱离单参考框架，做CASSCF和CASPT2/NEVPT2。

如果UHF能得到你所谓的“真实”的波函数，那你后面还做CC干啥？你想说的是自旋纯态吧，都U了，当然几乎不太可能是自旋纯态了。

作者
Author: sobereva 时间: 2020-2-19 13:26

zealot_cdp86 发表于 2020-2-19 10:49
UHF波函数问题有些细节还是比较模糊想请教一下。

1. 一直以为UHF最终会得到“真实”的稳定波函数。但 ...

UHF只是提供参考态轨道，之后做CCSD又产生了新的波函数，UHF波函数自旋污染大不代表UCCSD波函数自旋污染就大（而且CCSD不像UMP那样，UHF自旋污染大的话UMP波函数也逃不过）。显然最后能量的意义、好坏不能拿UHF的说事。
而且自旋污染大小和能量准确性并没有绝对的关系，诸如对称破缺方式算双自由基不可能没有显著的自旋污染。

作者
Author: mutron 时间: 2020-2-19 21:36

zealot_cdp86 发表于 2020-2-19 10:49
UHF波函数问题有些细节还是比较模糊想请教一下。

1. 一直以为UHF最终会得到“真实”的稳定波函数。但 ...

小小补充一下，要是碰到有不同稳定波函数解的，我个人觉得应该在当前电子结构下做几何优化（再在大基组下获得比较准确电子能量），根据玻尔兹曼分布看哪一个成分比较大，再根据具体反应讨论

作者
Author: zealot_cdp86 时间: 2020-3-5 22:54
抱歉，最近处理文章耽搁了这部分。
附件是R2-H体系的两个解，HF都是stable的，HF能量差了5.4kcal/mol，最后算的CCSD(T)-F12/cc-pvdz结果只差了0.89kcal/mol。不知道这个结果有没有可扩展性，但是似乎CCSD的能量差相比HF会缩小。

继续请教个问题，为什么说这个时候应该用多参考框架？你是指如果HF可以converge到两个stable波函数，这个现象暗示这个体系用多参考更好？还是T1诊断的结果？

对多参考框架理解非常有限，想请教一下，如果要验证你的假设，使用多参考方法，对于上面的体系，你推荐哪一个？如果是CASSCF，如何选活性空间？（2,2）? 谢谢

我一直的理解是HF提供唯一的参考态波函数（"真实波函数"），然后CCSD计算的结果取决于HF波函数的质量。

作者
Author: zealot_cdp86 时间: 2020-3-5 22:55

sobereva 发表于 2020-2-19 13:26
UHF只是提供参考态轨道，之后做CCSD又产生了新的波函数，UHF波函数自旋污染大不代表UCCSD波函数自旋污染 ...

谢谢！

作者
Author: zealot_cdp86 时间: 2020-3-5 23:03

mutron 发表于 2020-2-19 21:36
小小补充一下，要是碰到有不同稳定波函数解的，我个人觉得应该在当前电子结构下做几何优化（再在大基组下 ...

是这样，一般在计算过程中，我想大家都是得到一个stable就可以了，很少有特意去寻找有没有不同稳定的波函数，所以很难碰到多个稳定波函数的情况。如果有一些诊断方式，如T1诊断（多参考），那会给我们一个提醒，这样碰到的就会多一些。

另外，这个体系是在M06-2X/6-311G(d,p)下优化出来的，这里主要考虑体系一致性（DFT+basis set），同时有计算效率的考量（例如多环分子的优化，太大算不动）。

BTW，谢谢！

作者
Author: sobereva 时间: 2021-10-1 14:13

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

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

今日更新了Multiwfn，使用最新版本Multiwfn产生记录UHF/UKS轨道的fch文件，在unfchk将之转换成chk后，Gaussian可以通过guess=read正常读取里面的alpha和beta轨道了

经反复研究，发现关键在于fch里必须有这部分信息

Num ILSW I 1
ILSW I N= 1
1

复制代码

否则Gaussian在guess=read时不会从中读beta轨道。之前版本Multiwfn产生的fch未考虑这一点，如今考虑了。ILSW第一个值代表波函数类型，1代表real UHF，更多含义从Gaussian程序员手册里可以查到。

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