Amesp中SCF不收敛的解决办法（修订版）

• Jul 2, 2023 修订：增加fch2amo, py2amesp和qchem2amesp介绍

在使用量子化学软件时基本上都需要进行自洽场（SCF）迭代计算，一些时候会遇到SCF不收敛的情况，在这里将详细介绍Amesp软件中解决SCF不收敛时的办法，其中大多数关键词都是在“>scf”模块中设置。

1 增大迭代圈数

在Amesp中默认的SCF迭代圈数为125圈，这在大多数情况下是足够的，而有一些体系未能在125圈内收敛，且有收敛的趋势时，可以采用增大迭代圈数的办法，例子如下所示：

! b3lyp def2-SVP d3bj
>scf
 maxcyc 300
end
>xyz 0 1
C      -1.09929085    0.22458629    0.00000000
H      -0.74263642   -0.78422372    0.00000000
H      -0.74261801    0.72898448    0.87365150
H      -0.74261801    0.72898448   -0.87365150
H      -2.16929085    0.22459947    0.00000000
end

注意，这个方法主要应用于有收敛趋势的体系，如果出现震荡以及没有收敛趋势的时候，切勿盲目使用，圈数也不宜过大，一般300以内即可，切勿几千几万的设置！

2 能级移动法

能级移动法的原理是拉大HOMO-LUMO gap，使得占据轨道和空轨道的混合减弱，对于大多数不收敛的情况笔者十分推荐首先尝试该方法，尤其是含有过渡金属的体系，具体的设置为:

>scf
 vshift 500
end

vhift后面的数值根据具体需要进行设定，设置的数值为整数n，一般来说，数值越大效果越好，但是也会使得迭代圈数增加，笔者推荐的值为n=300~1000之间。

3 使用其他初猜

SCF是否容易收敛，初猜是一个十分关键的因素。在Amesp中，默认的初猜是Harris，当使用默认的初猜SCF没法收敛时，可以使用其他初猜，在Amesp中提供的初猜有：harris，huckel，core，gwh，read。初猜的具体设置方式为：

>scf
 guess huckel
end

其中read为读取存储在mo文件中的波函数作为初猜，这个关键词可以用来实现小基组投影，因为小基组更容易收敛，因此可以先使用小基组收敛后再投影到大基组。比如先使用3-21G做收敛后，再将3-21G计算得到的波函数作为def2-SVP的初猜。除了使用小基组，也可以使用该体系的阳离子体系的波函数作为初猜，因为电子数目越少，SCF越容易收敛。另外也可以在原来几何结构的基础上稍微调整一点，比如小范围改变键长键角二面角等，收敛后再作为原来结构的初猜。

4 DIIS方法

DIIS的全称为Direct Inversion in the Iterative Subspace，是由Pulay发展的对于SCF收敛十分有效的方法，其原理是使用前n步的Fock矩阵线性组合成一个新的Fock矩阵，这些线性组合的系数通过极小化一个误差函数得到。除了组合Fock矩阵的CDIIS，在Amesp中也支持组合密度矩阵的EDIIS和ADIIS，默认的为前半段使用ADIIS，后半段使用CDIIS，通过修改相应的关键词，可以使用不同的组合，这里举个例子：

>scf
 diis ediis
end

diis后面的关键词可以书写adiis，ediis，cdiss，off。其中adiis表示SCF前半段使用ADIIS，后半段使用CDIIS。ediis表示SCF前半段使用EDIIS，后半段使用CDIIS。cdiis则是全程使用CDIIS。off则是关掉diis。当SCF不收敛的时候，可以尝试使用不同的组合，或者关掉diis。除了更换不同的组合，在Amesp中也可以设置子空间的大小，例子为：

>scf
 ndiis 10
end

其中默认的大小为18，当遇到不收敛的时候可以适当地增大或者减小这个值。在Amesp迭代初期并不会开启DIIS，这是因为要避免一开始距离结果较远的Fock矩阵（或者密度矩阵）的干扰，默认当误差函数的值小于diistol=0.1的时候开启DIIS，用户可以通过如下的设置来修改开启DIIS的时机：

>scf
 diistol 0.05
end

5 二次收敛方法

常规的SCF求解是使用迭代法求解，迭代法的形式为x_{i+1}=f(x_{i})，除了迭代法，还有直接极小化的方法：min f(x)，量化软件中的几何结构优化就是这么做的，是通过改变自变量x（对于几何结构优化来说就是原子坐标，对于SCF来说就是波函数）来使得f(x)处在极值点。在Amesp中，可以使用能量对波函数的一阶导数以及二阶导数来搜索极小值。具体设置方式为：

>scf
 soscf 1st
end

soscf后面的关键词可以写1st，2nd，off。其中1st是使用一阶导数，这与ORCA中的soscf类似，而2nd是使用二阶导数，与Gaussian中的QC以及ORCA中的TRAH方法类似。笔者推荐使用1st，效果不错（尤其是对开壳层体系，U和RO）且每次循环的耗时与迭代法基本一致，而2nd则耗时更高，它每次循环都需要额外迭代求解一个线性方程，在Amesp中是调用了CP-SCF的代码来实现的。

6 增加积分精度

在SCF过程中，数值精度也会影响到收敛的情况。在Amesp中可以通过增加电子积分的精度以及DFT格点的精度来提高数值稳定性从而改善收敛情况。具体的例子为：

! b3lyp def2-SVP d3bj grid5
>scf
 accint 12
 incfock off
end
>xyz 0 1
C       -1.09929085    0.22458629    0.00000000
H       -0.74263642   -0.78422372    0.00000000
H       -0.74261801    0.72898448    0.87365150
H       -0.74261801    0.72898448   -0.87365150
H       -2.16929085    0.22459947    0.00000000
end

其中accint是设置电子积分精度的，默认为10^(-10)，即筛选掉小于10^(-10)的积分，上述的设置为10^(-12)。grid5是设置DFT格点的，默认为grid3，对于不太好收敛的情况，可以尝试设置grid4或者grid5。而incfock off是关掉Incremental Fock。以上的三种设置都能增加数值精度，能够增大SCF收敛的概率，但是也会相应地增加耗时。

7 阻尼和温度展宽

在Amesp迭代初期，会使用阻尼方法，其原理为混合第n步和n-1步的密度矩阵作为第n步的密度矩阵来计算Fock矩阵：

D_{n} = damp*D_{n-1} + (1-damp)*D_{n}

Amesp中damp默认为0.7，当启动DIIS时会关闭阻尼，修改damp的操作为：

>scf
 damp 0.65
end

温度展宽的关键词为:

>scf
 fermit 500
end

其中fermit后面的数字（整数）为温度（K），默认为0，一般设置300-1000即可。

8 从其他量化程序传轨道给AMESP

借助擅长传轨道的MOKIT程序，可从常见的量子化学程序传轨道给Amesp，达到SCF 1圈极速收敛。以Gaussian为例，假设我们完成了一个复杂体系的UHF计算，得到了一个hard.fch文件，执行

fch2amo hard.fch
a2m hard.amo

第一行是生成Amesp输入文件hard.aip和波函数文件hard.amo；第二行是用Amesp自带的a2m小程序将amo转化为mo文件（类似高斯的chk文件）。然后提交Amesp计算任务即可。

也支持从PySCF程序直接导出相关文件。假设读者已在PySCF中做完了一个SCF计算，结果保存在mf对象中，在.py脚本中添加以下两行

from mokit.lib import py2amesp
py2amesp(mf, aipname='hard.aip')

便可导出hard.aip和hard.amo文件。

另外，Q-Chem程序也可以产生fch文件，不过格式与高斯的fch文件有些不同，MOKIT采用了一个专门的API来处理，即在Python中运行

from mokit.lib import qchem2amesp
qchem2amesp('hard.fch', 'hard.aip')

便可产生hard.aip和hard.amo文件。

9 波函数稳定性分析

收敛后的波函数不一定是稳定的，这种情况会影响能量以及HF/DFT后续的计算，因此很多情况下需要进行波函数稳定性分析，Amesp支持R, U和RO三种情况的波函数稳定性分析，其关键词为：

>scf
 stable on
end

10 总结

当出现SCF不收敛的情况时，可以尝试本文中所提到的方法，1~7多种方式可结合使用，但不能保证SCF能100%收敛，当实在无法收敛时，可考虑更换方法基组；或采用方式8从其他量化程序中读入相同计算级别下的轨道。