RHF/RKS存在内部不稳定性/多个稳定解吗？

Accelerator

在检查波函数稳定性时，通常要么是RKS->UKS的不稳定性，要么是UKS->另一个UKS的内部不稳定性，那么是否存在RKS->另一个RKS的内部不稳定性呢？
比如，考虑两个距离很远、耦合很弱的片段构成的复合物A……B。如果是UKS的内部不稳定性，从自旋密度来考虑，有一类情况是对应单电子到底位于哪个片段上的问题；可能单电子在A上是一个稳定的波函数，单电子在B上是一个参数空间上的鞍点；也有可能两种情况都是稳定的UKS波函数。那么如果是RKS的情形，有些情况下，可能会有A2+……B2-、A……B、A2-……B2+这几种情况。在这种情况下，单参考态方法会如何表现呢？是否会存在波函数不稳定/多个稳定波函数的情况？
想到这一问题的背景是考虑双电子转移反应。对于伴随着单电子转移的反应，我们通过分析UKS波函数，或是RKS/UKS的交叉点，可以带来许多认识，而转移前后两个态耦合很弱的双电子转移反应如果要借助单参考态方法研究，就会涉及到单参考态方法会如何描述它们的RKS/RHF波函数的问题了。

wzkchem5

RKS->RKS内部不稳定性肯定是存在的，比如从HOMO激发两个电子到LUMO上的态，一般应该都有RKS->RKS内部不稳定性。
不过如果问有没有体系有RKS->RKS内部不稳定性却没有RKS->UKS不稳定性，就比较难说了。可以确认的是即使有也极其罕见，但是不知道能不能严格证明没有。
不管怎么说，因为RKS->UKS不稳定性的概率远大于RKS->RKS内部不稳定性的概率，所以估计大部分量化程序即使支持检查RKS->RKS内部不稳定性，也会先检查RKS->UKS不稳定性，然后如果发现不稳定，很可能就跳过后续检查，直接退出或者直接尝试优化稳定UKS波函数了。

Accelerator

wzkchem5 发表于 2022-11-12 21:33
RKS->RKS内部不稳定性肯定是存在的，比如从HOMO激发两个电子到LUMO上的态，一般应该都有RKS->RKS内部不稳定 ...

对于通过双电子激发构造出来的“不稳定波函数”有一个疑问，它是可以通过变分得到的吗？

另外感觉这个问题似乎跟自相互作用问题有关联，SIE或许会导致1L提到的几个态“融合”成一个电子比较离域的态？但是RHF/RKS参数空间的结构究竟是什么样，我目前没有一个清晰的感觉。

zjxitcc

我经常算键解离，共价键解离时容易出现RHF多个解，但这些解全是定性错误，所以平时一般不深究这其中哪个RHF解是内部稳定的（因为反正都用不着，都是选择用稳定的UHF波函数）。我想起来遇见过一个例子是wzkchem5说的双激发类型：RHF/STO-3G算氢气分子，输入文件如下

1. %chk=h2_rhf.chk
   
2. #p RHF/STO-3G nosymm int=nobasistransform
   
3. 
   
4. Title Card
   
5. 
   
6. 0 1
   
7. H   0.0   0.0   0.0
   
8. H   0.0   0.0   2.0
   
9. 
   
10. --Link1--
    
11. %chk=h2_rhf.chk
    
12. #p RHF chkbasis nosymm int=nobasistransform guess(read,alter) geom=allcheck
    
13. 
    
14. 1 2
    

复制代码

第一步是正常RHF计算，第二步中我们把HOMO, LUMO对调，这个两步任务会得到2个不同的RHF能量，均是1圈收敛

1. SCF Done:  E(RHF) = -0.783792654860     A.U. after    1 cycles
   
2. SCF Done:  E(RHF) = -0.541280622322     A.U. after    1 cycles

复制代码

说明都是极值点。接着传轨道到PySCF中检验各种波函数稳定性

1. bas_fch2py h2_rhf.fch

复制代码

去掉生成的.py文件末尾3行的注释符号，再加一句mf.stability(external=True)，提交，1圈收敛，结果显示

1. <class 'pyscf.scf.hf.RHF'> wavefunction is stable in the internal stability analysis
   
2. <class 'pyscf.scf.hf.RHF'> wavefunction has a real -> complex instability
   
3. <class 'pyscf.scf.hf.RHF'> wavefunction has a RHF/RKS -> UHF/UKS instability.

复制代码

我挺惊讶-0.5412806 a.u.这个高能RHF解居然是内部稳定的。

另外一个问题是这些波函数（尤其是能量更高的稳定解）能否通过变分得到。这个我也不知道咋说比较好，我个人有多种理解：推导出FC=SCE公式的过程我觉得是用了变分的。但是具体解FC=SCE的过程却未必是变分的，常见的SCF迭代算法（无论是DIIS、轨道旋转一阶梯度、轨道旋转二阶梯度）我认为都不能称为变分的，只有2*2 Jacob轨道旋转才能称为变分的轨道优化算法，其每圈能量必然下降（或不变），且必然收敛到极小值点，不会收敛到鞍点上，在这种变分的轨道优化算法下，不可能得到 不稳定解，顶多只能得到高能的稳定解。

如果换个相近的问题：能否通过常见的SCF收敛手段得到这些解？那应该是可以的，就是尝试片段组合波函数构建初猜、scf=qc等。如果实在难以构造，那得借助MOM/delta SCF方法使其收敛到想要的态，然后再检验轨道梯度是否为零、波函数是否稳定。

wzkchem5

zjxitcc 发表于 2022-11-12 17:26
我经常算键解离，共价键解离时容易出现RHF多个解，但这些解全是定性错误，所以平时一般不深究这其中哪个RHF ...

2x2 Jacobi旋转为什么必然会收敛到极小值点？
假设一个SCF已经收敛的体系的orbital Hessian的对角元全是正值，但是orbital Hessian的最小本征值是负值。那么按我理解，这个体系无论怎么做2x2 Jacobi旋转，都会使能量升高，但是沿着orbital Hessian最小本征值的本征矢方向扰动波函数仍然可以使得能量降低。

zjxitcc

wzkchem5 发表于 2022-11-13 01:05
2x2 Jacobi旋转为什么必然会收敛到极小值点？
假设一个SCF已经收敛的体系的orbital Hessian的对角元全是 ...

我的理解是2x2旋转一定会找到下降的那个方向，比如对于RHF/GVB波函数，我推导过电子能量是轨道旋转角的一元四次（三角）函数，对一元四次方程我可以求出所有解析解，然后取能量下降最多的那个旋转角，所以我感觉“orbital Hessian最小本征值的本征矢方向”一定逃不出解析解，一定不会2x2旋转漏掉，不知道我的理解是否正确。
（当然实际体系不止一个轨道旋转角，每次只能求1个最优角度，多元四次方程是耦合在一起的）

wzkchem5

zjxitcc 发表于 2022-11-12 18:14
我的理解是2x2旋转一定会找到下降的那个方向，比如对于RHF/GVB波函数，我推导过电子能量是轨道旋转角的一 ...

嗯，如果允许同时旋转不止一对轨道，那么我同意一定能找到一个能量极小值点。
但是Jacobi rotation一般指的是每次只转一对轨道吧？此时就有可能出现这样的情况，按某种特定的方式同时旋转若干对轨道可以让能量降低，但是单独旋转其中的任意一对都会导致能量升高。也就是不同轨道之间的耦合项导致了能量可以降低

zjxitcc

wzkchem5 发表于 2022-11-13 01:59
嗯，如果允许同时旋转不止一对轨道，那么我同意一定能找到一个能量极小值点。
但是Jacobi rotation一般 ...

这确实有可能

hebrewsnabla

补充一下，高斯也是可以计算RHF internal instability的，stable关键词有int, RRHF这些选项。但是不如pyscf清晰，pyscf可以对每一种HF/KS指定计算internal, external, real2complex。

zjxitcc

hebrewsnabla 发表于 2022-11-13 18:34
补充一下，高斯也是可以计算RHF internal instability的，stable关键词有int, RRHF这些选项。但是不如pyscf ...

我没注意到高斯也可以检测RHF内部稳定性一直以为高斯不行