求助ORCA DLPNO-CCSD(t)-F12计算能量异常问题

pet · 发表于 Post on 6 day ago

本帖最后由 pet 于 2025-8-16 03:18 编辑

各位老师好，我在用 ORCA 计算 C₆H₁₂O₃ 阳离子体系的单点能时遇到了一些问题，请教一下各位老师。（实际上是为了计算同一分子不同构象的电离能，目前中性分子没有问题，所以这里简化为阳离子能量计算的问题）

流程如下：

先在 Gaussian 16 中用 M06-2X-D3/def2-TZVP 优化了不同构象。
再用 ORCA 计算单点能，关键词为：! DLPNO-CCSD(T)-F12 cc-pVDZ-F12 TightPNO cc-pVDZ-F12-CABS cc-pVTZ/C tightSCF SCFConvForced noautostart miniprint，其中一个构象（06i）的能量值明显异常。波函数稳定性测试显示The stability analysis shows that the wavefunction is stable。

图中给出了几个构象在 DFT 和 DLPNO-CCSD(T)-F12/cc-pVDZ-F12 下的能量对比，其中 06i 为异常结构。

另外，由于最近发现 M06-2X 优化结果有时不太理想，又用 B3LYP-D3BJ/6-311+G(d,p) 对 06i 重新优化，并计算单点（06i-b）。这两个结构确实有差异，DFT结果相差2.5kcal/mol，但是按说不至于导致DLPNO-CCSD(T)-F12 结果相差这么多吧。

附件中分别是m062x（06i-stable.out）和b3lyp（06i-b.out）优化完的结构在orca中计算单点的输出文件。请老师们帮忙分析一下可能的原因。

06i-stable.out (122.55 KB, 下载次数 Times of downloads: 4)

06i-b.out (117.44 KB, 下载次数 Times of downloads: 0)

sobereva · 发表于 Post on 6 day ago

多一事（DLPNO）不如少一事，是我的话就用CCSD(T)/cc-pVTZ，结合MP2估计的TZ-QZ基组校正量了，这样大小的体系算起来很轻松

我不推荐M06-2X优化。但凡B3LYP-D3(BJ)优化不理想的情况（下文），推荐优先考虑wB97XD（或者ORCA里的wB97X-D3，或组合方法wB97X-3c）
谈谈量子化学研究中什么时候用B3LYP泛函优化几何结构是适当的
http://sobereva.com/557（http://bbs.keinsci.com/thread-17899-1-1.html）

zjxitcc · 发表于 Post on 6 day ago

试试看加上

%basis
PCDTrimBas Overlap
PCDTrimAuxJ Coulomb
PCDTrimAuxJK Coulomb
PCDTrimAuxC Coulomb
PCDThresh -1
end

复制代码

结果有无改善

pet · 发表于 Post on 4 day ago

sobereva 发表于 2025-8-16 05:20
多一事（DLPNO）不如少一事，是我的话就用CCSD(T)/cc-pVTZ，结合MP2估计的TZ-QZ基组校正量了，这样大小的体 ...

谢谢sob老师的回答，关于CCSD(T)/cc-pVTZ结合MP2估计的TZ-QZ基组校正量这种外推方式我也有点疑问，我在ORCA中测试C3H6O3，关键词
! ExtrapolateEP2(3/4,cc,MP2) tightSCF noautostart miniprint
%maxcore 5000
%pal nprocs 16 end
竟然算了三个多小时，相同的体系我用CCSD(T)/jul-cc-pV(T+d)Z 也才算了40min，总觉得是不是哪里有问题，，，这个耗时我根本不敢去算C6H12O3阳离子的体系。

test.out (130.32 KB, 下载次数 Times of downloads: 0)

pet · 发表于 Post on 4 day ago

zjxitcc 发表于 2025-8-16 18:44
试试看加上
结果有无改善

谢谢老师的回复，结果还是一样.....

zjxitcc · 发表于 Post on 4 day ago

pet 发表于 2025-8-18 16:39
谢谢老师的回复，结果还是一样.....

是否可以上传一下06i-stable.out任务在加了3L建议的关键词后的计算输出文件？

pet · 发表于 Post on 4 day ago

zjxitcc 发表于 2025-8-18 16:42
是否可以上传一下06i-stable.out任务在加了3L建议的关键词后的计算输出文件？

好的，已传

06i-3.out (151.13 KB, 下载次数 Times of downloads: 2)

zjxitcc · 发表于 Post on 4 day ago

本帖最后由 zjxitcc 于 2025-8-18 20:08 编辑

这个体系从中性单重态失去1个电子变成阳离子二重态时，最可能失去1个电子的地方有两处：（1）C=O双键中那个O的孤电子对；（2）O-O过氧键中某一个O的孤电子对。因此当前二重态体系大概率会有2个UHF稳定解，代表两种物理含义：
（1）C=O双键中的O失去了1个电子，剩余1个单电子仍然在这个O附近。
（2）O-O过氧键中某个O失去了1个电子，剩余1个单电子仍然在这附近。
你上传的06i-stable.out和06i-3.out文件都是收敛到了情况（2），它的UHF能量远高于（1）的UHF能量，也就是说通常情况下我们应该取UHF解（1）进行后续计算。无论你用UHF DLPNO-CCSD(T)，还是UHF UCCSD(T)，由于它们都是从UHF出发，你都可能会碰到CCSD(T)能量异常，碰到问题的概率约等于50%。

我们先构造和使用4种SCF初猜，分别进行UHF/cc-pVDZ-F12计算，4个gjf文件如下

06i_uhf_gau.gjf (3.95 KB, 下载次数 Times of downloads: 4) ,

06i_frag1_uhf.gjf (4.36 KB, 下载次数 Times of downloads: 5) ,

06i_frag2_uhf.gjf (4.36 KB, 下载次数 Times of downloads: 6) ,

06i_frag3_uhf.gjf (4.36 KB, 下载次数 Times of downloads: 4)

为方便阅读，展示第一个gjf文件内容

%chk=06i_uhf_gau.chk
%mem=160GB
%nprocshared=64
#p UHF gen nosymm int=nobasistransform scf(xqc,maxcycle=500) stable=opt
UHF/cc-pVDZ-F12 sp calculation
1 2
C 2.20560700 3.00976100 -1.10959300
...

复制代码

我这里用Gaussian做UHF计算，比较顺手（世上所有使用高斯基组的量化程序都可以收敛出一样的解，用哪个程序都行）。算完后汇总，4个任务得到2种UHF结果

06i_uhf_gau.log: SCF Done: E(UHF) = -458.546387235 a.u. after 18 cycles
06i_frag1_uhf.log: SCF Done: E(UHF) = -458.482233098 A.U. after 31 cycles
06i_frag2_uhf.log: SCF Done: E(UHF) = -458.482233098 A.U. after 32 cycles
06i_frag3_uhf.log: SCF Done: E(UHF) = -458.546387234 A.U. after 23 cycles

复制代码

其中-458.482233 a.u.对应你06i-stable.out文件中的UHF能量。可以看到另一个UHF解-458.546387 a.u. 足足低了40.2 kcal/mol，应当取这个解进行后续计算。我们可以直接利用这个波函数，无需在ORCA中从零开始算起、无需在ORCA中多番尝试去找另一个UHF解。运行以下三行命令

formchk 06i_uhf_gau.chk 06i_uhf_gau.fch
fch2mkl 06i_uhf_gau.fch
orca_2mkl 06i_uhf_gau_o -gbw

复制代码

此举将会产生ORCA相关文件06i_uhf_orca.inp, 06i_uhf_orca.gbw, 06i_uhf_orca.mkl，其中已写好体系坐标，基组数据，UHF关键词和UHF轨道数据。formchk是Gaussian自带小程序，fch2mkl是MOKIT的小程序，orca_2mkl是ORCA自带小程序。我们只需打开inp文件将前几行修改为DLPNO-CCSD(T)关键词

%pal nprocs 32 end
%maxcore 5000
! UHF DLPNO-CCSD(T)-F12 TightPNO cc-pVDZ-F12-CABS cc-pVTZ/C TightSCF

复制代码

核数和内存可以根据你自己机器情况修改。这里不用写cc-pVDZ-F12，因为基组数据在fch2mkl传过来的时候已经自动写好了。提交ORCA任务，过几秒打开输出文件，可以看到UHF瞬间收敛到了更低的那个解

----------------------------------------D-I-I-S--------------------------------------------
Iteration Energy (Eh) Delta-E RMSDP MaxDP DIISErr Damp Time(sec)
-------------------------------------------------------------------------------------------
*** Starting incremental Fock matrix formation ***
*** Initializing SOSCF ***
---------------------------------------S-O-S-C-F--------------------------------------
Iteration Energy (Eh) Delta-E RMSDP MaxDP MaxGrad Time(sec)
--------------------------------------------------------------------------------------
1 -458.5463872330064987 0.00e+00 8.74e-09 7.14e-08 6.49e-09 3.1
*** Restarting incremental Fock matrix formation ***
2 -458.5463872330067261 -2.27e-13 5.89e-09 3.35e-08 3.19e-08 2.9
*** Gradient check signals convergence ***
*****************************************************
* SUCCESS *
* SCF CONVERGED AFTER 2 CYCLES *
*****************************************************

复制代码

迅速进入DLPNO-CCSD(T)计算步骤，后面比较耗机时，任务我kill掉了。你自己阅读、理解我的解答后可自行计算。

PS1. 若想查看两个UHF解的单电子分别在哪里，可运行以下Python脚本

from mokit.lib.gaussian import uno
uno('06i_uhf_gau.fch')

复制代码

产生文件06i_uhf_gau_UNO.fch，用GaussView/Multiwfn打开此fch文件看单占据轨道。

PS2. 虽然该二重态体系在UHF下有两个解，但它在UM062X下可能有两个解、也可能只有一个解，结论不可简单照搬。原则上你之前算的01i~05i也应当做这种检查。由于此体系电子结构略为复杂，我不建议你采用基组外推的方式进行计算，那会包含多个SCF步骤然后把自己绕晕进去。

PS3. 上述技巧多次用到了MOKIT，若您使用这些技巧并发表文章，除引用ORCA外也建议引用MOKIT。引用MOKIT的50篇已发表文章可以看https://jeanwsr.gitlab.io/mokit-doc-mdbook/citing.html

pet · 发表于 Post on 4 day ago

zjxitcc 发表于 2025-8-18 20:00
这个体系从中性单重态失去1个电子变成阳离子二重态时，最可能失去1个电子的地方有两处：（1）C=O双键中那 ...

非常感谢老师您细致的回复！一定认真学习并正确引用~

我能够理解整个流程，但是关于您的gjf文件还有两个疑问：
1. 06i_frag1_uhf.gjf 和 06i_frag2_uhf.gjf 应该是对应您提到的第二种失电子情况：从O-O过氧键中某一个O的孤电子对上失去电子，但是为什么gjf的片段2里只包含OH呢？
2. 我尝试把OOH包含进frag2，设置自旋多重度为1 2 1 1 0 2以及1 2 0 2 1 1 , 我理解的第二种情况对应OOH上失电子，第一种情况对应frag1失电子，但是这两种算出来都是收敛到了高能量的情况（结果上看与您的 06i_frag1_uhf.gjf 和 06i_frag2_uhf.gjf完全一致），，，这是为什么呢？为什么第一种已经排除了OOH失电子，也还是不能收敛到羰基O失电子的情况去呢？

zjxitcc · 发表于 Post on 3 day ago

本帖最后由 zjxitcc 于 2025-8-19 17:24 编辑

pet 发表于 2025-8-18 21:47
非常感谢老师您细致的回复！一定认真学习并正确引用~

我能够理解整个流程，但是关于您的gjf文件还有 ...

存在06i_frag1_uhf.gjf和06i_frag2_uhf.gjf这两个文件，是因为当时我对“O-O过氧键中某个O失去了1个电子”没有100%的把握，我怀疑它可以进一步细分为12号O失去1个电子、13号O失去1个电子两种不同情况，所以片段2只含OH两个原子。等我算完才发现这里没有两种细分，实际上是12号O失去1个电子，说成“O-O过氧键中某个O失去了1个电子”也没问题。所以这两个gjf文件对应的都是我在8L回答中说的情况（2）。

06i_frag3_uhf.gjf对应的是我在8L回答中说的情况（1）。而06i_uhf_gau.gjf这个文件，就是想看看Gaussian默认初猜能收敛出啥来，一个常规计算，没想到它直接收敛到了较低能量的UHF解。你说把OOH划分成一个片段当然也可以，也对应情况（2）。只有像我给出的06i_frag3_uhf那样划分才对应情况（1）。

pet · 发表于 Post on 前天 20:17

zjxitcc 发表于 2025-8-19 15:50
存在06i_frag1_uhf.gjf和06i_frag2_uhf.gjf这两个文件，是因为当时我对“O-O过氧键中某个O失去了1个电子 ...

好的，非常感谢您

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[ORCA] 求助ORCA DLPNO-CCSD(t)-F12计算能量异常问题

评分 Rate