计算化学公社

标题: 为什么我使用ORCA5.0.3做对称破缺计算得到的OSS单点能，它是显著大于CSS单点能呢？ [打印本页]

作者
Author: 13799725819 时间: 2024-3-12 00:41
标题: 为什么我使用ORCA5.0.3做对称破缺计算得到的OSS单点能，它是显著大于CSS单点能呢？
本帖最后由 13799725819 于 2024-3-19 17:53 编辑

主要想研究分子中双自由基特征，或者分子中某单原子（氮）上相互垂直的两个轨道中的单电子自旋方向问题，下面进行了三个计算：（输出文件见附件）
1. 对称破缺计算（open shell singlet, OSS）输入文件：
%pal nprocs 12 end
%maxcore 2000
! RI-B2PLYP D3BJ RIJCOSX def2-TZVPD def2/J def2-TZVPD/C ExtremeSCF
%scf
HFTyp UHF
BrokenSym 1,1
end
* xyz 0 1
坐标。。。
*

------
输出文件：（附件：1.out）
<S**2>(High-Spin) = 2.1019
<S**2>(BrokenSym) = 0.3096
E(High-Spin)    =  -505.006873 Eh（也就是“TOTAL SCF ENERGY”，下同）
E(BrokenSym)    =  -505.015144 Eh
RI-MP2 CORRELATION ENERGY:    -0.479005150 Eh
Dispersion correction          -0.003429712
FINAL SINGLE POINT ENERGY    -505.497578519065

2. 对称破缺计算，结合波函数稳定性分析（open shell singlet, OSS）输入文件：
%pal nprocs 12 end
%maxcore 2000
! RI-B2PLYP D3BJ RIJCOSX def2-TZVPD def2/J def2-TZVPD/C ExtremeSCF
%scf
HFTyp UHF
BrokenSym 1,1
STABPerform true # default false
STABRestartUHFifUnstable true # restart the UHF-SCF if unstable
STABNRoots 3 # number of eigenpairs sought
STABMaxDim 3 # Davidson expansion space = MaxDim * NRoots
STABMaxIter 100 # maximum number of Davidson iterations
STABNGuess 4096 # size of initial guess matrix: 4096 x 4096
STABDTol 0.0001 # convergence criterion from iteration to iteration
STABRTol 0.0001 # convergence criterion max residual norm
STABlambda +0.5 # mixing parameter
STABORBWIN -1, -1, -1, -1, -1, -1, -1, -1 # defines the donor / acceptor spaces
# 4 parameters for RHF
# 8 paramters for UHF (4 alpha, 4 beta)
# orbital window, -1 refers to automatic determination
STABEWIN -5.0, 5.0 # lower and upper cutoff in Eh for automatic freezing
end
* xyz 0 1
相同坐标。。。
*

------
输出文件：（附件：2.out）
<S**2>(High-Spin) = 2.1019
<S**2>(BrokenSym) = 0.3096
E(High-Spin)    =  -505.006873 Eh
E(BrokenSym)    =  -505.015144 Eh
RI-MP2 CORRELATION ENERGY:    -0.479005150 Eh
Dispersion correction          -0.003429712
FINAL SINGLE POINT ENERGY    -505.497578519065

3. 闭壳层单重态（close shell singlet, CSS）计算输入文件：（程序默认计算RHF/RKS）
%pal nprocs 12 end
%maxcore 2000
! RI-B2PLYP D3BJ RIJCOSX def2-TZVPD def2/J def2-TZVPD/C ExtremeSCF
* xyz 0 1
相同坐标。。。
*

------
输出文件：（附件：3.out）
Total Energy    :       -505.01302386 Eh（也就是“TOTAL SCF ENERGY”）
RI-MP2 CORRELATION ENERGY:    -0.492645753 Eh
Dispersion correction          -0.003429712
FINAL SINGLE POINT ENERGY    -505.509099325939

ORCA5.0.3计算的OSS能量值比CSS高了 7.23 kcal/mol。

*******************************
下面使用G09.D01进行计算：(R/U)B2PLYPD3(BJ)/def2-TZVPD
close shell singlet：
SCF Done:  E(RB2PLYPD3) =  -505.016241133    A.U. after 17 cycles
E2(B2PLYPD3) = -0.4926741699D+00 E(B2PLYPD3) = -0.50550891530248D+03
输出文件总结部分@：
HF=-505.0162411\MP2=-505.5089153

比ORCA5.0.3计算的结果大了：0.12 kcal/mol

open shell singlet：
SCF Done:  E(UB2PLYPD3) =  -505.018047107    A.U. after 28 cycles
<S**2>= 0.3688
E2(B2PLYPD3) = -0.4844192704D+00 E(B2PLYPD3) = -0.50550246637785D+03
输出文件总结部分@：
HF=-505.0180471\MP2=-505.5024664

比ORCA5.0.3计算的结果大了：-3.07 kcal/mol

G09.D01计算的OSS能量值比CSS高了 4.05 kcal/mol。

作者
Author: wjc404 时间: 2024-3-12 08:12
本帖最后由 wjc404 于 2024-3-12 08:37 编辑

对称破缺计算不是黑箱，需要有一定的经验。你的体系期望对称破缺的原因是啥（过渡金属反铁磁耦合/双自由基体系/反芳香性共轭/...）?

在对称破缺的计算中，可以尝试让程序生成spin natural orbitals或unrestricted natural orbitals（翻翻orca手册看怎么生成），然后用molden或multiwfn等可视化软件，查看占据数绝对值接近1的轨道是否合适。

作者
Author: zjxitcc 时间: 2024-3-12 09:47
本帖最后由 zjxitcc 于 2024-3-12 09:50 编辑

提供的信息不足，无法下结论。还没到看轨道、轨道占据数的时候，先自己想一下回答这些问题
（1）你用的是双杂化泛函B2PLYP，因此1个输出文件里至少有2个能量，第一个是SCF收敛处的自洽泛函的能量，第二个是真正的双杂化泛函能量。你没说你对两个任务比较的是哪个能量。
（2）没说自洽部分的波函数到底稳不稳定。
（3）没说到底差多少（不仅要有数值，还要有单位）。最好列个简单的表格，把自洽泛函的能量、自洽泛函的<S**2>和双杂化泛函能量都列出来。
这些都理清楚了，就搞清楚自己在算什么了。

即使你准确地回答了上述三点，即使你的计算全部合理，也有可能因为体系存在较强多参考特征，导致RDFT-B2PLYP的微扰部分失效，相关能异常（理论方法局限所致），此时可以低于其他任何方法的能量，不足为奇。

作者
Author: 13799725819 时间: 2024-3-12 14:14

wjc404 发表于 2024-3-12 08:12
对称破缺计算不是黑箱，需要有一定的经验。你的体系期望对称破缺的原因是啥（过渡金属反铁磁耦合/双自由基 ...

打算研究双自由基体系，或者分子中某个原子的垂直轨道中单电子的自旋问题；可以看“分子中某个原子的垂直轨道中单电子的自旋”相反的情况吗？问题又做了补充，麻烦再给看一下，谢谢！

作者
Author: 13799725819 时间: 2024-3-12 14:21

zjxitcc 发表于 2024-3-12 09:47
提供的信息不足，无法下结论。还没到看轨道、轨道占据数的时候，先自己想一下回答这些问题
（1）你用的是 ...

问题又做了补充，麻烦再给看一下，谢谢！
（1）我取的“FINAL SINGLE POINT ENERGY”应该是B2PLYP最终的能量值；要关注的是这个值。
（2）2.中做了波函数稳定性分析，我使用的是ORCA手册中的案例方法，这样做可行吗？结果显示是稳定的。
（3）这里计算的OSS和CSS的能量值相差了 > 5 kcal/mol，更新的问题中列出了这些值，麻烦再给瞧瞧，谢谢！
体系可能会有显著的多参考特征。

作者
Author: wjc404 时间: 2024-3-12 16:46
本帖最后由 wjc404 于 2024-3-15 08:26 编辑

你这个体系是三氟乙酰基氮宾，如果是单重态，显然有较强的multireference特征。你用的又是ORCA软件，强烈建议用NEVPT2来算，比KSDFT更好，在开启RI的情况下用def2-TZVP也不昂贵。至于活性空间的确定，可以借助https://gitlab.com/jxzou/mokit这个软件。

我用ORCA-5.0.3软件算了一下，在FIC-NEVPT2(4,4)/def2-QZVP水平下做几何优化（得用数值梯度），优化后的坐标如下：
  C    0.079167 -0.094701 0.002454
  N    1.180981 0.451291 0.310224
  O    -0.435414 0.996482 -0.336702
  C    -0.491641 -1.497699 0.005018
  F    0.412289 -2.336459 0.492546
  F    -1.581753 -1.541329 0.757217
  F    -0.804661 -1.860958 -1.230756

4e4o的活性空间用multiwfn绘制轨道，如下：
(, 下载次数 Times of downloads: 12)

作者
Author: 13799725819 时间: 2024-3-12 23:19

wjc404 发表于 2024-3-12 16:46
你这个体系是三氟乙酰基氮宾，如果是单重态，显然有较强的multireference特征。你用的又是ORCA软件，强烈建 ...

谢谢你的建议！

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