加速SCF收敛可以尝试的一些手段

Eudaimonia

先祝大家新年快乐！
起因：本人在使用cp2k对Co2NiS4(001)晶面，3*3*3超晶胞，15Å真空层，共112个原子模型进行扩胞后的几何优化时(输入文件由社长的Multiwfn生成)，发现SCF迭代的收敛速度奇慢无比。先反复确定了磁性等设置后，最后输入文件在在Broy./Diag.默认迭代上限128步后，距收敛线还有约2个数量级的差异，于是打开.inp文件进行了一番尝试；在本文测试中，由于是单台工作站且人在工作站前，懒得生成Eenergy计算文件，直接使用Geo_OPT的第一次SCF迭代结果作为判定标准，完成第一次迭代后手动暂停进行下一个测试；
本文使用了两个计算模型：
a.Co2NiS4(001)晶面，3*3*3超晶胞，15Å真空层，共112个原子,该模型在默认参数下很难收敛；
b.单原子催化剂的Fe-N4-C模型，该模型在默认参数下很容易收敛；
Note: 由于本人基础较差，部分解释可能是胡诌，希望各位老师不吝赐教！

在inp文件中，控制SCF迭代主要由&DIAGONALIZATION及&MIXING下的关键词进行控制，社长贴心的在关键词后添加了注释方便理解；
&DIAGONALIZATION部分控制对角化算法的选择，此处不涉及尝试，提示中有大体系可使用DAVIDSON
&MIXING字段下有Method，Alpha及NBROYDEN三个参数；

尝试1：更改method：
Method是选择SCF迭代过程中，新旧矩阵的替代方法，由Multiwfn生成的inp文件默认使用BROYDEN_MIXING，本人尝试过PULAY_MIXING，在该体系下，该方法收敛更加困难。

尝试2：更改NBROYDEN：
NBROYDEN后注释为Number of previous steps stored for the actual mixing scheme，我理解为是综合参考之前X步的矩阵来混合生成新矩阵，cp2k程序默认值为4，社长认为设置较大值有利于SCF迭代，故Multiwfn生成的输入文件该处默认值为8，我继续增大了NBROYDEN的值，进行测试；
测试结果如下所示，在测试结果中，以Coverge steps为评价标准，越小的Coverge steps表明在SCF迭代过程中使用更少的步数到达了收敛限(3E-6)，最后列出计算出的总能量，以对比计算结果是否一致；

可以看到，无论是Co2NiS4模型还是Fe-N4-C模型，进一步增大NBROYDEN能够加速收敛，且计算出的能量基本一致，在设置到24时，能加速约15%（1/6~1/7）；

尝试3：更改Alpha值；
在增大NBROYDEN的基础上，本人又对该参数进行了测试；Alpha后的注释为Mixing xx of new density matrix with the old one，我理解为是新生成的SCF矩阵中由多少是新生成的，越大说明该SCF矩阵新生成的部分越多；查阅手册发现此处实际还有个Beta参数，但实际功能与Alpha相似，新旧矩阵混合由rho_mix(g) =rho_in(g) + alpha*g^2/(g^2 + beta^2)*(rho_out(g)-rho_in(g))决定，Alpha在分子上，Beta在分母上；
测试结果如下所示：



可以看到，计算出的能量基本是一致的，对于之前收敛较为困难的Co2NiS4体系，适当降低Alpha的值有利于SCF的收敛，但Alpha值降的过低反而有害；而对于本身就容易收敛的Fe-N4-C体系，无论增大或是降低Alpha的值都不利于SCF收敛；

总结：
对于目前的两个算例中
1.进一步增加NBROYDEN的值总是有利于SCF的收敛，cp2k默认的NBROYDEN值为4，Multiwfn生成的inp文件默认为8，可能8还是较为保守，大家可以在进行计算前再进一步增大NBROYDEN的值，观察是否能进一步加速收敛，或者激进一点直接倍增NBROYDEN的值（需要更多的算例支持，建议还是以测试结果为准）；
2.更改Alpha值需要慎重进行测试，在收敛困难的体系中可能可以倍增收敛速度，但是一定要先测试再改，不当的设置结果会引起收敛速度的倒退！

其他：
请确认你的SCF收敛是收敛速度很慢，存在收敛趋势时进行测试，在没有收敛趋势时（在某个值反复震荡时），一定要先排查kpoint，自旋多重度，磁矩，smearing等是否设置恰当！
两个测试的原始输入文件在附件中；

卡开发发

上述混合算法其实也是大多数固体的计算程序所采用的方法，对VASP、DMol3等诸多程序其实也适用，这里只给出上述两个程序的对照，其他程序大家可以自行对照手册（或者欢迎讨论）。由于学艺不精，所以只好对这些参数大致上做一些解释：

1、alpha其实也就是density mixing的参数（对照VASP的AMIX和DMol3的density mixing->charge）。alpha通常取小有利于减小SCF震荡，但会导致收敛缓慢。一般小的alpha有利于gap比较小的体系，能够避免能级顺序频繁交换，过小的alpha可能会导致一些数值问题。
2、beta出现在kerker或broyden混合中（对照VASP的BMIX和DMol3的preconditioner->q0），能明显看到公式的形式就是不同的G分量的贡献引入的不同的混合，一般来说低频分量收敛会比高频分量来的慢，对于金属性体系（比如某些表面）的电荷震荡(charge sloshing)问题会有较好的帮助。
3、nbroyden，其实也就是DIIS的size（对照VASP的MAXMIX和DMol3的DIIS size），就是混合之前的步数，有些平面波程序甚至能用到20～40，比较大的情况有利于收敛，但是对一些SCF比较病态的体系可能有时候反而小一些（甚至不开）才行，大了会引起震荡。

最后要说的就是，除了（个人经验上）磁矩设置不合理导致迭代不到基态引起的能量与mixing参数呈现依赖，一般而言这些迭代的参数不会影响计算结果，只影响SCF收敛快慢。当SCF出现难收敛的情况时，这通常是最后的方案，建议优先检查模型结构、磁性、涉及到精度的参数。通常也不必做特别复杂的测试，保障模型合理精度合理的情况下遇到不收敛再根据上述经验适当调整为好。上述算法的原理部分应该可以在Richard M Martin的书中找到，了解一些原理总还是要比胡乱调整更为合适。

sobereva

根据我的一些测试，一般NBROYDEN设8差不多够了，再大一般也不会令收敛明显更快，所以Multiwfn默认用了相对保守的8。但如果碰上难收敛的情况可以尝试设比如10、12试试，如果还不好收敛，也没太大必要尝试设更大了，而应该考虑别的法子。ALPHA一般0.3或0.4就行，之前试过更大或更小通常都不好。

Eudaimonia

卡开发发 发表于 2022-2-1 20:42
上述混合算法其实也是大多数固体的计算程序所采用的方法，对VASP、DMol3等诸多程序其实也适用，这里只给出 ...

谢谢卡开发发老师的解释，请问一下Richard M Martin的书指的是：
Interacting Electrons: Theory and Computational Approaches 还是 Electronic Structure: Basic Theory and Practical Methods？

Eudaimonia

sobereva 发表于 2022-2-1 21:18
根据我的一些测试，一般NBROYDEN设8差不多够了，再大一般也不会令收敛明显更快，所以Multiwfn默认用了相对 ...

嗯嗯，对于我做过的绝大多数体系来说，Multiwfn默认的参数设置都是合理的

卡开发发

Eudaimonia 发表于 2022-2-2 09:38
谢谢卡开发发老师的解释，请问一下Richard M Martin的书指的是：
Interacting Electrons: Theory and Co ...

这块看Electronic Structrure，这本书侧重于基本原理，比如不同的基组算法等；而Interacting Electrons这本偏向于理论方面，主要讲凝聚态的GW DMFT和QMC等计算方法。

Eudaimonia

卡开发发 发表于 2022-2-2 16:05
这块看Electronic Structrure，这本书侧重于基本原理，比如不同的基组算法等；而Interacting Electrons这 ...

谢谢老师，我去找一下相关资源，仔细研读一下

asd9998876

刚刚在cp2k的google group上面看到了回复，还可以增加smearing的温度，不知道这一点的影响怎么样呢

You should further decrease ALPHA and increase ELECTRONIC_TEMPERATURE in &SMEAR. Only activate SURFACE_DIPOLE_CORRECTION once the SCF is converged without it and you can use a restart wavefunction.

With ELECTRONIC_TEMPERATURE 3000 and ALPHA 0.05 the first SCF converges.

asd9998876

刚刚在cp2k的google group上面看到了回复，还可以增加smearing的温度，不知道这一点的影响怎么样呢

You should further decrease ALPHA and increase ELECTRONIC_TEMPERATURE in &SMEAR. Only activate SURFACE_DIPOLE_CORRECTION once the SCF is converged without it and you can use a restart wavefunction.

With ELECTRONIC_TEMPERATURE 3000 and ALPHA 0.05 the first SCF converges.