请教一个反应能垒使用CCSD(T)/aug-cc-pvtz和CBS-QB3计算结果相差巨大的原因

Melvin

各位老师好，最近计算了HCN+CH3→（TS）→CH3NCH这个反应，在CCSD(T)/aug-cc-pvtz//M062X/6-311G+(2df,2p)级别计算得到的能垒为19kcal/mol，而在CBS-QB3计算得到的能垒为7.6kcal/mol，两个结果相差甚远。
1.文献中有提到过“自旋污染会影响含 CN 自由基计算的准确性”，结果中CCSD(T)/aug-cc-pvtz下的T1诊断为0.0372，自旋污染也很大“S**2 before annihilation  0.9699, after0.7748”，是不是由于自旋污染导致的两个方法计算误差巨大？
2.如果需要通过需要消除自旋污染来获得准确的结果，具体应该怎么操作呢？
3.这种情况热力学组合方法的结果还有参考价值吗？

感谢各位老师的解答！！

wzkchem5

CBS-QB3和CCSD(T)都会受自旋污染影响，但我感觉两者影响的方向应该一致，而且不会有10kcal/mol这么大。其实我觉得更可能是其中一个计算没有收敛到正确的波函数，可以做个stability analysis看看

sobereva

这么小体系没必要CBS-QB3，建议用G4
用cc-pVQZ比aug-cc-pVTZ明显更适合当前问题
当前的差异和自旋污染的联系不大

Melvin

wzkchem5 发表于 2022-7-22 21:34
CBS-QB3和CCSD(T)都会受自旋污染影响，但我感觉两者影响的方向应该一致，而且不会有10kcal/mol这么大。其实 ...

感谢王老师，每次问题都得到了您的解答，真的受益匪浅。
按照您的建议，在计算CCSD(T)/aug-cc-pvtz时使用guess=read读取了M062x/6-311G+(2df,2p)的波函数计算，结果和CBS-QB3和G4的结果近似.

还想请教王老师一个小问题，波函数分析结果是The wavefunction is stable under the perturbations considered，但这个过渡态在CCSD(T)/aug-cc-pvtz水平下S2高达0.98，T1诊断为0.037，这种情况如果后续计算反应速率常数需要考虑消除自旋污染的影响吗？

Melvin

sobereva 发表于 2022-7-23 14:49
这么小体系没必要CBS-QB3，建议用G4
用cc-pVQZ比aug-cc-pVTZ明显更适合当前问题
当前的差异和自旋污染的 ...

感谢Sob老师的建议，后边我都换用G4和cc-pVQZ计算了。
还有一个小问题就是这个反应的过渡态S2为0.98，T1诊断为0.037，如果后续计算反应速率常数需要考虑消除自旋污染吗？

wzkchem5

Melvin 发表于 2022-7-27 09:20
感谢王老师，每次问题都得到了您的解答，真的受益匪浅。
按照您的建议，在计算CCSD(T)/aug-cc-pvtz时使 ...

自旋污染这个数问题不大，T1 diagnostic有点大但是也还可以（开壳层体系T1 diagnostic本来就偏大，有个别文献认为对于开壳层体系T1 diagnostic大于0.044才算大）。所以结果还是可以相信的。

Melvin

wzkchem5 发表于 2022-7-27 16:32
自旋污染这个数问题不大，T1 diagnostic有点大但是也还可以（开壳层体系T1 diagnostic本来就偏大，有个别 ...

太感谢了，谢谢王老师的解答！

xxtnenu

wzkchem5 发表于 2022-7-22 21:34
CBS-QB3和CCSD(T)都会受自旋污染影响，但我感觉两者影响的方向应该一致，而且不会有10kcal/mol这么大。其实 ...

请问老师，“S**2 before annihilation  0.9699, after0.7748” 这里S**2的值应该是0.9699还是0.7748.

Error on total polarization charges =  0.00966
SCF Done:  E(UB3LYP) =  -1489.65419124     A.U. after   24 cycles
            NFock= 24  Conv=0.36D-08     -V/T= 2.0018
<Sx>= 0.0000 <Sy>= 0.0000 <Sz>= 0.0000 <S**2>= 1.0058 S= 0.6206
<L.S>=  0.00000000000   
KE= 1.486982112146D+03 PE=-5.120539511139D+03 EE= 1.349334096389D+03
Annihilation of the first spin contaminant:
S**2 before annihilation     1.0058,   after     0.0480
Leave Link  502 at Tue Oct 14 16:20:28 2025, MaxMem= 10737418240 cpu:             531.3 elap:              49.6
(Enter /mnt/home/g16/l801.exe)
DoSCS=F DFT=T ScalE2(SS,OS)=  1.000000  1.000000

这是学生某次计算的输出文件最后一次scf done的内容，关键词为stable ub3lyp/6-311+g(d,p) scrf=(solvent=water) nosymm em=gd3bj integral=ultrafine，如果想知道当前方法基组下的S**2的值是看 <S**2>= 1.0058 还是看 after     0.0480呢？在输出文件的最后还看到***********************************************************************
Stability analysis using <AA,BB:AA,BB> singles matrix:
***********************************************************************
1PDM for each excited state written to RWF  633
Ground to excited state transition densities written to RWF  633

Eigenvectors of the stability matrix:

Excited state symmetry could not be determined.
Eigenvector   1:  1.088-?Sym  Eigenvalue= 0.0877741  <S**2>=0.046
     52A -> 53A        0.99976

Excited state symmetry could not be determined.
Eigenvector   2:  2.283-?Sym  Eigenvalue= 0.0956305  <S**2>=1.053
     48A -> 53A       -0.20511
     49A -> 53A       -0.14914
     51A -> 53A        0.94810

Excited state symmetry could not be determined.
Eigenvector   3:  2.341-?Sym  Eigenvalue= 0.1075253  <S**2>=1.120
     49A -> 53A        0.80285
     50A -> 53A       -0.44417
     51A -> 53A        0.14027
     49B -> 56B        0.10473
     52B -> 53B       -0.13654
The wavefunction is stable under the perturbations considered.
Leave Link  914 at Tue Oct 14 16:21:16 2025, MaxMem= 10737418240 cpu:            1729.7 elap:              46.5
(Enter /mnt/home/g16/l601.exe)
Copying SCF densities to generalized density rwf, IOpCl= 1 IROHF=0.
这里的s2则是<S**2>= 1.120，请问老师此处只是计算过程而非最终的s2？ 谢谢老师

zjxitcc

xxtnenu 发表于 2025-10-15 07:09
请问老师，“S**2 before annihilation  0.9699, after0.7748” 这里S**2的值应该是0.9699还是0.7748.

...

（1）不用看S**2 before annihilation xxx, after xxx，直接看输出文件最后一次SCF Done下方1~2行处的<S**2>，当然，它也等于before annihilation紧跟的数字，在你的例子里就是1.0058。
（2）波函数稳定性检验步骤中
Eigenvector ...  Eigenvalue ... <S**2>=...
这里<S**2>没有意义，之所以会输出是因为波函数稳定性功能通常会借用CIS和TDDFT计算代码和打印结果的代码，<S**2>作为一个不想要的附带输出。还是那句话：直接看输出文件最后一次SCF Done下方1~2行处的<S**2>。

xxtnenu

zjxitcc 发表于 2025-10-15 09:00
（1）不用看S**2 before annihilation xxx, after xxx，直接看输出文件最后一次SCF Done下方1~2行处的， ...

谢谢老师解答！

xxtnenu

zjxitcc 发表于 2025-10-15 09:00
（1）不用看S**2 before annihilation xxx, after xxx，直接看输出文件最后一次SCF Done下方1~2行处的， ...

还请问下老师，问什么有时候结果输出scf done后没有输出s**2值呢？
关键词为：
%oldchk=opt.chk
%chk=stable.chk

#p stable=opt guess=read rb3lyp/6-311+g(d,p) scrf=(solvent=water) nosymm em=gd3bj integral=ultrafine

zjxitcc

xxtnenu 发表于 2025-10-16 07:23
还请问下老师，问什么有时候结果输出scf done后没有输出s**2值呢？
关键词为：
%oldchk=opt.chk

UHF/UKS计算会输出<S**2>，而RKS计算不会输出，因为RKS的<S**2>一定严格为0，理论特点决定，无需在输出文件中打印该值。有一种特殊情况，是在stable=opt过程中发现RKS波函数不稳定，继而程序自动去优化UKS波函数，才会在UKS能量底下1~2行出现<S**2>；你在11L展示的例子是RKS稳定的。

如果想系统性学习 单行列式波函数的<S^2>是什么含义、如何计算，让基础更扎实一点点，可以阅读Szabo的《Modern Quantum Chemistry》一书前3章。

xxtnenu

zjxitcc 发表于 2025-10-16 08:54
UHF/UKS计算会输出，而RKS计算不会输出，因为RKS的一定严格为0，理论特点决定，无需在输出文件中打印该值 ...

收到，谢谢老师指点，学生这就看一下教材内容