计算化学公社

标题: 关于较大分子脱CO进行NEVPT2计算时活性空间选择和单点能计算准确性 [打印本页]

作者
Author: mayang 时间: 2025-8-24 19:53
标题: 关于较大分子脱CO进行NEVPT2计算时活性空间选择和单点能计算准确性
本帖最后由 mayang 于 2025-8-24 19:55 编辑

各位老师们好，我目前计算的体系存在一些点的T1诊断值偏高,需要用多参考方法计算一下，这是其中一个点。我先用mokit对TS11和W10分别自动选择了活性空间，分别为（7e,7o）和（9e,9o），查看了自动选择的结果，发现并没有包括要断裂的C3-C4的sigma键，于是在 _asrot.fch中按照这篇帖子http://bbs.keinsci.com/forum.php?mod=viewthread&tid=53831&highlight=%B9%EC%B5%C0的方法进行了轨道的选择，包括了图片中的活性空间，对于TS11和W10都是（11e,11o）,最后算出的能垒为2.3kcal/mol(NEVPT2/cc-PVTZ，没有加零点能校正)，比CCSD(T)/cc-PVTZ低了6kcal/mol，请问各位老师我的计算结果是否合理？结构优化、频率计算和IRC计算方法是B3LYP/6-311

作者
Author: sobereva 时间: 2025-8-25 07:47
这个反应的静态相关凭经验看强不到哪去。当静态相关不是太强的情况，用多参考二阶微扰方法的结果还不如CCSD(T)（哪怕活性空间选择得当），很可能还不如revDSD-PBEP86-D3(BJ)、wB97M-V之类顶级泛函好，非必要时用多参考方法完全是走弯路

北京科音高级量子化学培训班（http://www.keinsci.com/KAQC）详细讲CASSCF和多参考计算部分的这页幻灯片专门强调了这一点：
(, 下载次数 Times of downloads: 0)

基组名必须严谨写完整，不得就写个6-311，否则谁知道你用的6-311G还是6-311G(3d2f,2pd)

作者
Author: wzkchem5 时间: 2025-8-25 10:22
T1诊断值的判断标准0.02只适用于闭壳层，因为提出T1诊断值的原文献只测了闭壳层体系。
开壳层体系有人提出过以0.044为判断标准，但用的人不多。不过总而言之，开壳层体系的T1诊断值>0.02未必表示CCSD(T)结果不可靠。

作者
Author: mayang 时间: 2025-8-25 10:23

wzkchem5 发表于 2025-8-25 10:22
T1诊断值的判断标准0.02只适用于闭壳层，因为提出T1诊断值的原文献只测了闭壳层体系。
开壳层体系有人提出 ...

我的是0.048

作者
Author: mayang 时间: 2025-8-25 10:27

sobereva 发表于 2025-8-25 07:47
这个反应的静态相关凭经验看强不到哪去。当静态相关不是太强的情况，用多参考二阶微扰方法的结果还不如CCSD ...

方法是b3lyp/6-311g(d,p)

作者
Author: mayang 时间: 2025-8-25 10:36

sobereva 发表于 2025-8-25 07:47
这个反应的静态相关凭经验看强不到哪去。当静态相关不是太强的情况，用多参考二阶微扰方法的结果还不如CCSD ...

老师好，这个过渡态的T1诊断值=0.048，审稿人认为能垒计算的值会有问题，如果不用多参考，是可以再跑几个单参考方法去回答他吗？

作者
Author: Uus/pMeC6H4-/キ 时间: 2025-8-25 11:00

wzkchem5 发表于 2025-8-25 10:22
T1诊断值的判断标准0.02只适用于闭壳层，因为提出T1诊断值的原文献只测了闭壳层体系。
开壳层体系有人提出 ...

这两个判断标准就是在http://bbs.keinsci.com/thread-25474-1-1.html提到的吧？确实http://bbs.keinsci.com/thread-5460-1-1.html测试的体系里T1诊断值没几个这么高的。

后面又看您在http://bbs.keinsci.com/thread-45602-1-1.html指出T1诊断只能表现HF轨道描述静态相关的能力，这样的话对此类“看似普通”的有机自由基反应，还有什么如http://bbs.keinsci.com/thread-19007-1-1.html提到的方法之类推荐用于多参考特征的判断呢？我有种刻板印象是论坛上很多人算的氧原子/单重态氧气/臭氧降解有机污染物的反应基本都要多参考方法，但是有很多做燃烧的用羟基自由基/过氧自由基/超氧自由基/超氧算提取氢的反应也用上了多参考方法……

作者
Author: sobereva 时间: 2025-8-25 11:25

mayang 发表于 2025-8-25 10:36
老师好，这个过渡态的T1诊断值=0.048，审稿人认为能垒计算的值会有问题，如果不用多参考，是可以再跑几个 ...

给出CCSD自然轨道占据数，有助于判断是否有必要多参考方法。不应只用T1一个指标说事。很多T1不小的体系用CCSD(T)照样能算得很好

用前面我说的泛函可以算算，单参考方法里也可以用OO-CCD(T)算算（包含了轨道优化过程，对静态相关的描述能力强于CCSD(T)。相当于Gaussian里的BD(T)）

另外，NEVPT4(SD)也可以试试，不过这个比NEVPT2贵非常非常多，可能得考虑减小活性空间、去掉非必要的活性轨道才能算得动（本身它对活性空间的要求倒也比NEVPT2更低）。

作者
Author: mayang 时间: 2025-8-25 12:33

sobereva 发表于 2025-8-25 11:25
给出CCSD自然轨道占据数，有助于判断是否有必要多参考方法。不应只用T1一个指标说事。很多T1不小的体系用 ...

好的，谢谢老师，我去算一算，如果用NEVPT4(SD)，活性空间应该怎么选呢，只涉及C-C的sigma和C-O键的pi键可以吗

作者
Author: wzkchem5 时间: 2025-8-25 13:28

Uus/pMeC6H4-/キ发表于 2025-8-25 11:00
这两个判断标准就是在http://bbs.keinsci.com/thread-25474-1-1.html提到的吧？确实http://bbs.keinsci.c ...

对。多参考特征多强算强，也和需要的精度以及计算量有关，如果臭氧要和一个100原子的分子反应，且有好几个可能的反应位点，那即便明知多参考特征强，也不得不用DFT算。反过来，对于特别小的体系，MRCISD+Q能跑得动的，即使CCSD(T)结果可靠，为求万无一失还是往往会上MRCISD+Q

作者
Author: mayang 时间: 2025-8-25 14:06

sobereva 发表于 2025-8-25 11:25
给出CCSD自然轨道占据数，有助于判断是否有必要多参考方法。不应只用T1一个指标说事。很多T1不小的体系用 ...

老师，这是自然轨道占据数，您可以帮我看看需不需要用多参考方法吗？
2.000015 2.000011 2.000010 2.000010 2.000009 2.000008
2.000006 2.000005 2.000003 1.984913 1.980282 1.977935
1.974843 1.973732 1.969959 1.969321 1.968002 1.967169
1.966534 1.963607 1.962308 1.959100 1.958046 1.954980
1.953239 1.952291 1.946929 1.943658 1.932464 1.927319
1.905124 0.997299 0.083411 0.061614 0.056387 0.048978
0.044904 0.029141 0.028020 0.026946 0.024497 0.023615
0.023007 0.022308 0.021996 0.021862 0.020899 0.020234
0.019751 0.018040 0.017020 0.016595 0.011624 0.011276
0.009914 0.009830 0.009368 0.008682 0.008467 0.008297
0.008051 0.007755 0.007418 0.007199 0.006851 0.006606
0.006484 0.006313 0.006257 0.005995 0.005613 0.005524
0.005324 0.005228 0.005013 0.004743 0.004545 0.004380
0.004318 0.004136 0.003987 0.003768 0.003646 0.003456
0.003321 0.003063 0.002863 0.002707 0.002632 0.002439
0.002352 0.002282 0.002142 0.002101 0.001968 0.001905
0.001861 0.001833 0.001810 0.001722 0.001680 0.001617
0.001588 0.001551 0.001495 0.001431 0.001330 0.001297
0.001270 0.001223 0.001185 0.001127 0.001036 0.001008
0.000988 0.000984 0.000947 0.000894 0.000880 0.000827
0.000802 0.000755 0.000735 0.000642 0.000616 0.000587
0.000550 0.000512 0.000490 0.000474 0.000452 0.000391
0.000377 0.000347 0.000332 0.000327 0.000313 0.000309
0.000222 0.000175 0.000134 0.000112 0.000107 0.000099
0.000088 0.000084 0.000069 0.000059 0.000054

作者
Author: wjc404 时间: 2025-8-25 17:41

mayang 发表于 2025-8-25 14:06
老师，这是自然轨道占据数，您可以帮我看看需不需要用多参考方法吗？
2.000015 2.000011 2.000 ...

谢谢分享。看上去是用不着使用MR方法的。单点能的话，如果是高斯可以用BD(T)，如果是ORCA可以用OOCCD(T)，最大的双激发振幅<0.1就OK。

作者
Author: mayang 时间: 2025-8-25 18:44

wjc404 发表于 2025-8-25 17:41
谢谢分享。看上去是用不着使用MR方法的。单点能的话，如果是高斯可以用BD(T)，如果是ORCA可以用OOCCD(T) ...

谢谢老师

作者
Author: mayang 时间: 2025-9-1 21:13

wjc404 发表于 2025-8-25 17:41
谢谢分享。看上去是用不着使用MR方法的。单点能的话，如果是高斯可以用BD(T)，如果是ORCA可以用OOCCD(T) ...

老师，您是怎么从自然占据轨道看出来的不用多参考方法呢

作者
Author: zjxitcc 时间: 2025-9-11 17:02
本帖最后由 zjxitcc 于 2025-9-13 15:25 编辑

如果我没理解错，你这里能垒是指TS11与W10的相对能量？能展示一下TS11和W10的UCCSD(T)/cc-pVTZ具体能量数值吗？

作者
Author: mayang 时间: 2025-9-17 19:50
本帖最后由 mayang 于 2025-9-17 19:53 编辑

zjxitcc 发表于 2025-9-11 17:02
如果我没理解错，你这里能垒是指TS11与W10的相对能量？能展示一下TS11和W10的UCCSD(T)/cc-pVTZ具体能量数值 ...

ccsd(t)
TS11:-383.4820859 Hartree
W10:-383.4953641 Hartree

nevpt2(11e，11o):
TS11:-383.50238 Hartree
W10:-383.50602 Hartree

作者
Author: mayang 时间: 2025-9-17 19:55

zjxitcc 发表于 2025-9-11 17:02
如果我没理解错，你这里能垒是指TS11与W10的相对能量？能展示一下TS11和W10的UCCSD(T)/cc-pVTZ具体能量数值 ...

ccsd(t):
TS11 -383.4820859;w10 -383.4953641
nevpt2(11e,11o):
TS11:-383.5023788
w10:-383.50601683

作者
Author: zjxitcc 时间: 2025-10-5 11:29
本帖最后由 zjxitcc 于 2025-10-5 11:35 编辑

这个问题比较复杂。可以分为两部分来理解

（1）UCCSD(T)计算
W10这个体系有两个UHF解，在cc-pVTZ基组下电子能量分别为 (单位: a.u.)

	solution 1	solution 2
UHF	-381.911153	-381.922008
UCCSD(T)	-383.495364	-383.490805

第一个UHF解是不稳定的，第二个是稳定的。你用了前者进行UCCSD(T)计算，这通常是不合理的。有一个经验规律是，不稳定的UHF解通常会导致更大的相关能，所以尽管你的UHF能量偏高，但你的UCCSD(T)能量却偏低（这不是什么好事）。可能你是用Gaussian/ORCA“写一行关键词，提交1个任务获得结果”这种黑箱方式做的计算，以为这样就可以了，没有检验UHF波函数稳定性。对于电子结构较为复杂的体系，即使是单参考计算也要加倍小心，无论使用哪个传统量化软件都要仔细检查计算过程和计算结果。建议先做UHF计算，检验UHF波函数稳定性，确保波函数稳定，然后读取UHF波函数进行UCCSD(T)计算，分成两三步。

而你给出的TS11分子只有一个UHF稳定解，UHF能量是-381.910078 a.u., UCCSD(T)能量是-383.482087 a.u., 所以这一基元步骤的UCCSD(T)电子能垒是5.47 kcal/mol。这一步如果弄错了，后面对比NEVPT2能垒就容易困惑

最后算出的能垒为2.3kcal/mol(NEVPT2/cc-PVTZ，没有加零点能校正)，比CCSD(T)/cc-PVTZ低了6kcal/mol

（2）NEVPT2计算
这一基元步骤所需的最小活性空间为CAS(3,3)，活性轨道包含目标C-C键成键轨道、反键轨道，以及体系自身二重态对应的1个单占据轨道；对应的活性电子也是3个。

(, 下载次数 Times of downloads: 0)

上图是W10的3个重要GVB轨道，可以用GaussView/Multiwfn打开文件W10_cc-pVTZ_uhf_uno_asrot2gvb22_s.fch进行可视化。在W10分子中目标C-C键还未断裂，成键轨道占据数为1.99，非常接近2.0，因此本段所说的“重要”不是根据活性轨道占据数来判断的，而是根据基元反应机理判断的。目标C-C键轨道不在32号轨道附近，需要用permute_orb()函数将其调换至31号、33号轨道位置。

理论上讲CAS(3,3)就足够了，可能你会疑惑NEVPT2(3,3)是否足够准确，这的确需要测试。我们可以使用更高阶的多参考方法例如FIC-MRCISD+Q或FIC-NEVPT4(SD)，这两个在ORCA里都有，在活性空间合理的情况下往往可以获得靠谱的结果。若活性空间增大，由于基组cc-pVTZ下空轨道数目较多，FIC-NEVPT4(SD)大概率算不动。这是我算的结果

(, 下载次数 Times of downloads: 0)

可以看到，基于CASSCF(3,3)的FIC-NEVPT2和FIC-NEVPT4(SD)能垒与UCCSD(T)能垒很接近。就算FIC-NEVPT2结果是凑巧，还有FIC-NEVPT4(SD)可以说明问题。如果一定要追求更大活性空间，需要小心地增加的活性轨道和活性电子。使用MOKIT的automr对W10和TS11分别进行自动CASSCF/cc-pVTZ计算，自动确定的活性空间分别是(7,7)和(9,9)，显然我们需要增大/减小活性空间使二者一致。

注意，无论看xxx_gvbN_s.fch文件（在这里N=21或22）还是xxx_CASSCF_NO.fch文件，都可以看出活性轨道主要成份，前者是后者的初始轨道。但如果要调换轨道实现活性空间的一致，则需要调换第一个文件，这是因为GVB pair轨道是局域在每一根化学键上的，化学意义十分清晰，而CASSCF自然轨道是离域在整个体系中的，往往作为结果查看，而不作为可调换的过程。GVB是automr做自动CASSCF计算默认的一个中间步骤，它也很重要，不能忽略它直接去看CASSCF自然轨道。待我们看完这两个文件中的轨道
W10_cc-pVTZ_uhf_uno_asrot2gvb22_s.fch
TS11_cc-pVTZ_uhf_uno_asrot2gvb21_s.fch
可以发现乙烯基取代五元环上3根C-C pi键是比较重要的轨道，占据数偏离2.0/0.0最多，同时它们也跟自由基轨道可能共轭

(, 下载次数 Times of downloads: 0)

所以优先将这3对轨道加入活性空间。注意每根化学键对应一对轨道（成键和反键轨道）、2个电子，因此在CAS(3,3)基础上扩增为CAS(9,9)，电子能垒变为4.38 kcal/mol（见上面表格数据），还算合理。若还想进一步扩大活性空间，应当优先考虑C-O pi键（因为它在GVB轨道里是偏离2.0/0.0次多的），从而可以构造CAS(11,11)。但这里有个额外的问题，如果进行CASSCF(11,11)轨道优化，W10体系计算会很顺利（能量丝滑下降，下降幅度小，很快收敛），但对TS11会发现CASSCF轨道变化较大，能量下降较多，这是因为在W10的GVB和CASSCF计算中C-O上始终只有一根C-O pi键和一根C-O sigma键，O的其他价电子是俩孤电子对，在TS11的GVB计算中也是这样；但在TS11的CASSCF计算中有两根C-O pi键、一根C-O sigma键、一对O孤电子，意思就是有一对O孤电子转变成了O->C 配位pi键，成键方式发生较大变化，需要进一步扩大活性空间，将O的孤电子对、C-O sigma等加入活性空间，问题趋于复杂。有一个折中的奇技淫巧是，既然轨道比较合理，我们可以不做轨道优化，直接进行CASCI(11,11)-NEVPT2计算，把O的俩孤电子对留在双占据空间里，不管它们怎么变化，这就是上述表格中FIC-NEVPT2(11,11) CIonly的来历，能垒为5.56 kcal/mol，也比较合理。

总结：该有机反应在单参考和多参考计算均合理的前提下，二者能垒非常接近（有1.0 kcal/mol左右的区别是正常的，但不会有6 kcal/mol的区别）。
PS1: NEVPT2方法有两个变种SC-NEVPT2和FIC-NEVPT2，后者稍微准确一点点，本回答用的也是FIC-NEVPT2方法。ORCA支持两种，PySCF目前仅支持SC-NEVPT2。
PS2: 如果使用MOKIT以及本文技巧发表文章，请记得引用MOKIT程序，点击此处可以查看引用MOKIT的53篇已发表文献。

计算涉及的输入文件和脚本见附件（没有提供GVB和CASSCF的fch文件以及输出文件，那些文件也比较大，不适合上传论坛。请阅读后自行练习计算，才能较好地理解）
(, 下载次数 Times of downloads: 2)

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