计算化学公社

 找回密码 Forget password
 注册 Register

Davidson 算法的preconditioner是如何推导出来的?

查看数: 6830 | 评论数: 23 | 收藏 Add to favorites 2
关灯 | 提示:支持键盘翻页<-左 右->
    组图打开中,请稍候......
发布时间: 2021-8-9 13:32

正文摘要:

本帖最后由 JohnCase 于 2021-9-1 04:09 编辑 求解大型稀疏对角矩阵特征值的Davidson算法,使用矩阵对角元素作为preconditioner, 以此计算出下一循环中的subspace basis,即 basis_new = (D - \lambda)^-1 *resid ...

回复 Reply

JohnCase 发表于 Post on 2025-9-3 18:18:06
SHENLIN 发表于 2025-9-3 16:54
是等于这个,就是要算的态的个数。为啥最好取这个,我看高斯前两圈取的都是nstate*3

感觉这是两段代码最不一样的地方
SHENLIN 发表于 Post on 2025-9-3 16:54:34
本帖最后由 SHENLIN 于 2025-9-3 17:03 编辑
JohnCase 发表于 2025-9-3 12:37
初猜数量是多少?算10个态的话,最好是设置前18个态作为初猜。

lowest等于size(residues, 2)吗

是等于这个,就是要算的态的个数。为啥最好取这个,我看高斯前两圈取的都是nstate*3
JohnCase 发表于 Post on 2025-9-3 12:37:01
本帖最后由 JohnCase 于 2025-9-3 14:23 编辑
SHENLIN 发表于 2025-9-2 10:50
再次感谢老师,不然这个bug不知道要找到什么时候。请问还有什么加速收敛的策略吗?我看每一圈计算norm的 ...

初猜数量是多少?算10个态的话,最好是设置前18个态作为初猜。

lowest等于size(residues, 2)吗
wzkchem5 发表于 Post on 2025-9-2 10:58:26
SHENLIN 发表于 2025-9-2 10:50
再次感谢老师,不然这个bug不知道要找到什么时候?请问还有什么加速收敛的策略吗?我看每一圈计算norm的 ...

大的特征值收敛慢是正常的,比如你算10个特征值,第10个特征值和第11个特征值的距离一般很接近,导致需要迭代很多轮才能把第10个特征值和第11个特征值之间的矩阵非对角元消干净。
迭代次数少一般是因为第一步求解比较多的特征值,收敛后再把高能量的解扔掉。比如第一步求解20个根,但从第二步开始只计算前10个根的residual vector,后10个根仍然放在Krylov子空间里,但不产生新的vector,也就是Krylov子空间按照20->30->40...的规律增加
SHENLIN 发表于 Post on 2025-9-2 10:50:57
本帖最后由 SHENLIN 于 2025-9-2 10:57 编辑
wzkchem5 发表于 2025-7-7 10:38
感觉大体还是有收敛趋势的,可能preconditioner实现得有问题?可以对角化一个对角矩阵,看是不是一步收敛 ...

再次感谢老师,不然这个bug不知道要找到什么时候。请问还有什么加速收敛的策略吗?我看每一圈计算norm的时候,前面小的特征值收敛比较快,后面大的收敛圈数要多得多,请问这里有什么改进的方法吗?另外我用高斯迭代的圈数好像要少很多,它是用了什么策略吗?
SHENLIN 发表于 Post on 2025-9-2 10:36:42
本帖最后由 SHENLIN 于 2025-9-2 10:56 编辑
JohnCase 发表于 2025-8-7 18:29
diag(i)是什么,i是nocc*vir这个维度吗

是这个维度,就是对角线元素。我还是没搞明白  之前preconditiner是这样写得
我看MRCC程序,把分母改成下面的这样

就没有问题了。我觉得是等价的啊

JohnCase 发表于 Post on 2025-8-7 18:29:19
SHENLIN 发表于 2025-7-31 13:48
do j = 1, lowest
            do i=1,nocc
                do a=1,nvir

diag(i)是什么,i是nocc*vir这个维度吗
JohnCase 发表于 Post on 2025-8-7 18:19:19
SHENLIN 发表于 2025-7-31 16:11
好的,谢谢。一定要设定max space吗?直接用最大迭代圈数控制可以吗?

讲道理不需要设定max space。或者说设定max space不能大于A矩阵大小
hebrewsnabla 发表于 Post on 2025-7-31 16:25:45
SHENLIN 发表于 2025-7-31 16:11
好的,谢谢。一定要设定max space吗?直接用最大迭代圈数控制可以吗?

我不清楚你的程序有没有trial vector的space,也就是个数的概念,如果没有的话,大概是求n个根就只有n个trial vector?那效率应该是比较低的。
如果有space的概念,每一步会扩大space,但是没有max space(超过了就重来),也可能会有收敛问题,但也不一定。
SHENLIN 发表于 Post on 2025-7-31 16:11:52
本帖最后由 SHENLIN 于 2025-7-31 16:16 编辑
hebrewsnabla 发表于 2025-7-31 12:21
可以参考一下这个实现 https://github.com/pyscf/pyscf/b ... nalg_helper.py#L291

你的max space是多少 ...

好的,谢谢。一定要设定max space吗?直接用最大迭代圈数控制可以吗?
SHENLIN 发表于 Post on 2025-7-31 13:48:06
       do j = 1, lowest
            do i=1,nocc
                do a=1,nvir
                ia=(i-1)*nvir+a
                denom = -eorb(a+nocc)+eorb(i)+eigenvalues(j)           
                ! 处理分母过小的情况
                if (abs(denom) < t) then
                    denom = sign(t, denom)
                end if           
                ! 计算修正项
                correction(ia, j) = residues(ia, j) / denom
                end do
            end do
        end do   我这样写preconditioner就是正确的,能正常收敛
     do j = 1, size(residues, 2)
        do i = 1, size(correction, 1)
            D = eigenvalues(j) - diag(i)
            
            ! 处理分母过小的情况
            if (abs(D) < t) then
                D = sign(t, D)
            end if
            
            ! 计算修正项
            correction(i, j) = residues(i, j) / D
        end do
    end do这样写就是错误的,一直收敛不了。这两者不是等价的吗?好奇怪
hebrewsnabla 发表于 Post on 2025-7-31 12:21:25
本帖最后由 hebrewsnabla 于 2025-7-31 12:23 编辑

可以参考一下这个实现 https://github.com/pyscf/pyscf/b ... nalg_helper.py#L291

你的max space是多少?

SHENLIN 发表于 Post on 2025-7-31 09:21:19
本帖最后由 SHENLIN 于 2025-7-31 10:44 编辑

谢谢老师指导,确实是diag的原因,因为我如果把矩阵元都求出来,然后取对角线元素,上面的例子16步左右能收敛。但是我在写CIS的过程中,为了节约存储,不是不求出H(ia,jb),而是直接求sigma向量HV吗?所以我的diag取的都是HF轨道能级差项,忽略了后面电子积分项(我看很多程序好像都是这么做的)。初猜和preconditioner用的都是这个,需要怎么调整呢?
do i=1,nocc
     do a=1,nvir
           ia=(i-1)*nvir+a
            diag(ia)=eorb(a+nocc)-eorb(i)
       enddo
enddo
wzkchem5 发表于 Post on 2025-7-7 10:38:07
SHENLIN 发表于 2025-7-7 09:36
占据轨道5,虚轨道44,nstate=3,为啥我写得程序一直收敛不了?
Number of occupied orbitals   :     5
...

感觉大体还是有收敛趋势的,可能preconditioner实现得有问题?可以对角化一个对角矩阵,看是不是一步收敛,这样验证preconditioner有没有实现对

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2026-2-19 07:13 , Processed in 0.171426 second(s), 25 queries , Gzip On.

快速回复 返回顶部 返回列表 Return to list