p450酶体系判断波函数稳定性问题

sylar

p450酶体系，包含一个铁卟啉环，是四重态简并到二重态，用pbe0/def2-svp速度略慢，所以试着使用r2scan-3c算，
初步解决收敛问题以后进行了波函数稳定性检测，首先试了二重态，输入关键部分如下：

!QMMM UKS r2scan-3c notrah verytightscf soscf
%pal nprocs 12 end
%maxcore 2000
%scf
AutoStart true
#BrokenSym 2,1
DirectResetFreq 1
DIIS Start 0.1 MaxIt 5 MaxEq 20 BFac 1.2 MaxC 15.0 end
SOSCFStart 0.0005
MaxIter 500
STABPerform true
end

得到如下结果，

      *** CONVERGENCE OF RESIDUAL NORM REACHED ***

-----------------------------
SCF STABILITY ANALYSIS RESULT
-----------------------------

UHF/UKS->UHF/UKS - singlet - internal

   Root    Eigenvalue (au)
     0      -0.000218
     1       0.006383
     2       0.014730

Stability Analysis indicates an UNSTABLE HF/KS wave function with one negative eigenvalue.

之后程序自动试图找到正确的波函数，出现以下结果，能量大幅震荡升高无法收敛

SCF ITERATIONS (restarted)
--------------------------
ITER       Energy         Delta-E        Max-DP      RMS-DP        Damp
  0  -4052.7736790308   0.000000000000 0.00450332  0.00002393 0.7000
  1  -4052.7734694511   0.000209579739 0.00622190  0.00003025 0.7000
  2  -4052.7726029985   0.000866452558 0.01454909  0.00005487 0.7000
  3  -4052.7675201063   0.005082892203 0.03534604  0.00014970 0.7000
  4  -4052.7314742913   0.036045815008 0.09607929  0.00033530 0.7000
  5  -4052.3062532910   0.425221000341 0.41520735  0.00162750 0.7000
  6  -4051.6304294008   0.675823890144 1.00571460  0.00149026 0.7000
  7  -4050.5839545901   1.04647481068228.91864806  0.06345000 0.7000
  8  -4030.5015474341  20.08240715602724.77646097  0.08446787 0.7000
  9  -3942.9818102237  87.519737210358936.62185775  1.94423335 0.7000
10  -3676.9187348417 266.063075382021287.12902571  0.59422251 0.7000
11  -3826.5017140618-149.582979220096347.89514449  0.84992277 0.7000
12  -3830.2612043399  -3.759490278107409.78491112  0.75552369 0.7000
13  -3789.3521073956  40.909096944322550.95593315  1.15240079 0.7000
14  -3607.7590637882 181.593043607428256.63755962  0.74372994 0.7000
15  -3737.7630277705-130.0039639822981115.59257333  1.95462190 0.7000
16  -3741.2749898141  -3.511962043629386.38975004  0.77957209 0.7000
17  -3827.7562264399 -86.481236625792330.39254056  0.55386042 0.7000
18  -3817.9756269938   9.780599446116749.75870352  1.54843944 0.7000
19  -3659.8650443401 158.110582653670290.08599715  0.64999569 0.7000
20  -3788.0274174555-128.162373115380256.64214134  0.62699408 0.7000
21  -3802.8090369334 -14.781619477938530.58767934  0.86925583 0.7000
22  -3772.8729205477  29.936116385730556.75873793  0.99573357 0.7000

之后尝试了四重态，输入关键词一样，得到以下结果，波函数在64轮后收敛，且稳定

      *** CONVERGENCE OF RESIDUAL NORM REACHED ***

-----------------------------
SCF STABILITY ANALYSIS RESULT
-----------------------------

UHF/UKS->UHF/UKS - singlet - internal

   Root    Eigenvalue (au)
     0       0.011282
     1       0.016302
     2       0.017150

Stability Analysis indicates a stable HF/KS wave function.

以此为初猜，无法收敛到二重态，同样能量升高并震荡


想请教一下各位老师，

1.结果中-0.000218这个值的意思是存在比当前波函数能量低 -0.000218 au的波函数吗？（已解决）
2. -0.000218 au换算下来才0.125kcal/mol,是不是说这个波函数可以用？（已解决）
3.有文献提到铁卟啉环CPDI体系用杂化泛函如b3lyp与CASSCF计算的自旋密度较好吻合，r2scan-3c是不是不适合这种体系？
4.主要目的是加快计算速度，大家有什么建议给我

谢谢大家

wzkchem5

不是，这是orbital Hessian的本征值，它和”稳定波函数的能量比当前波函数低多少“没有任何联系。只要orbital Hessian的本征值小于0，不管和0多接近，稳定波函数的能量都有可能显著低于当前波函数的能量，低多少都有可能。

北大-陶豫

类比一下结构优化的虚频，虚频数值的大小跟该结构还能再优化下来多少能量没有关系

sylar

wzkchem5 发表于 2022-2-21 17:07
不是，这是orbital Hessian的本征值，它和”稳定波函数的能量比当前波函数低多少“没有任何联系。只要orbit ...

明白了，谢谢！我补充了一些尝试的结果，老师你有空看一下有什么建议

sylar

北大-陶豫 发表于 2022-2-21 18:43
类比一下结构优化的虚频，虚频数值的大小跟该结构还能再优化下来多少能量没有关系

谢谢，您这么说理解更直观