PySCF程序包平均场计算的一些收敛技巧
平均场计算是 PySCF 程序包里优化得比较并全面的模块之一。在平均场模块里，PySCF支持 RHF, UHF, ROHF, GHF, RKS, UKS, ROKS, GKS 等一系列方法来研究闭壳层体系、开壳层体系、复数哈密顿量体系、相对论效应、溶剂化效应。同时 PySCF 提供了大量的辅助功能来帮助平均场计算收敛。以下我们通过一些例子来演示在 PySCF 里收敛平均场计算的技巧。
以下的例子在 PySCF-1.5 以上发行版均可使用。
首先我们使用 PySCF 程序包实现一个最简单的平均场计算。在编辑器里创建以下 Python 代码
benzene = ''' C 4.673795 6.280948 0. C 5.901190 5.572311 0. C 5.901190 4.155037 0. C 4.673795 3.446400 0. C 3.446400 4.155037 0. C 3.446400 5.572311 0. H 4.673795 7.376888 0. H 6.850301 6.120281 0. H 6.850301 3.607068 0. H 4.673795 2.350461 0. H 2.497289 3.607068 0. H 2.497289 6.120281 0. '''from pyscf import gto, scf
Create an molecule object
mol = gto.M(atom=benzene, basis='ccpvdz', verbose=4)
Create an mean-field object
mf = scf.RHF(mol)
Solve the SCF problem
mf.run()
然后在命令行里执行这个 Python 脚本就可以做一个最简单的平均场计算
python example.pyPySCF 默认的收敛判别标准非常严格，判别标准是能量收敛到 1e-10 Eh，能量梯度的 norm 收敛到 1e-5 。收敛参数可以通过设定平均场对象的 .conv_tol, .conv_tol_grad 来调整。大多数时候，设定
mf.conv_tol = 1e-6的收敛标准已经足够化学问题的研究了。
DIIS 方法
DIIS 迭代方法是 PySCF 里默认的平均场收敛算法。在上面的例子里，DIIS 算法进行了 7 轮迭代得到了收敛结果。
当体系比较复杂时，DIIS 算法的默认参数不一定能得到收敛的平均场结果。此时，有一些常规的参数可以帮助 DIIS 算法收敛，比如
1. 我们可以增加迭代步数
mf.max_cycle = 200来等待 DIIS 的结果慢慢移动到收敛的结果上。
2. 我们可以调整 level shift 的大小来改变体系 HOMO-LUMO gap。较大的 HOMO-LUMO gap 一般能避免收敛过程的震荡。
mf.level_shift = 0.1level shift 的单位是 hartree 。这个值一般不用特别大，0.2 左右对大多数体系应该就够用了。当体系的基组有较大的线性相关时，过大的 level shift 不利于数值稳定性。
3. 我们可以调整 DIIS 空间的大小来给予 DIIS 算法更大的自由度来搜索解，比如
mf.diis_space = 12DIIS 的收敛还可以被一系列别的参数影响，更多的设定可以参考帮助文档 http://pyscf.org/pyscf/scf.html#hartree-fock
初始猜测
在平均场计算中，很多困难体系对初始猜测十分敏感，合适的初始猜测可以有效地帮助平均场收敛。PySCF 的平均场默认的初始猜测是基于 ANO-RCC 基组构造的原子密度矩阵的叠加。这个初始猜测对平衡结构，或较复杂的问题（如开壳层，非平衡结构，带有少量电荷的体系）都有不错的表现。除了 ANO 作为初始猜测以外，PySCF 还提供了一系列方法对平均场计算的初始猜测进行调整。
有一些体系需要特殊的初始猜测才能收敛到正确的态上，比如铁磁或反铁磁的初始猜测。此时，就需要对 PySCF 的初始猜测进行调整。PySCF 支持构造好一个密度矩阵，然后传给平均场对象作为初始猜测进行计算。例如
import numpy as np
from pyscf import gto, scf
mol = gto.M(atom='''
H 0 0 0 ; H 0 0 0.5;
H 0 0.5 0.5; H 0 0.5 0 ''', basis='sto3g')
print('Default initial guess:')
mf = scf.UHF(mol)
mf.kernel()
mf.analyze()
print('\nAnti-ferromamagnetic initial guess:')
dm = np.zeros((2,4,4))
dm[0,0,0] = dm[0,2,2] = 1
dm[1,1,1] = dm[1,3,3] = 1
mf.kernel(dm)
mf.analyze()
可以看到默认的初始猜测收敛到了能量较高的态上，而反铁磁的初始猜测收敛到能量较低的态上。
Default initial guess:
converged SCF energy = -0.521792378127812<S^2> = 4.6629367e-152S+1= 1
**** MO energy ****
alpha | beta alpha | beta
MO #1 energy= -1.19356543325823 | -1.19356543325823 occ= 1 | 1
MO #2 energy= -0.0986706198268508 | -0.09867061982685 occ= 1 | 1
MO #3 energy= 0.416828469148984 | 0.416828469148984 occ= 0 | 0
MO #4 energy= 2.37630222058633 | 2.37630222058633 occ= 0 | 0
** Mulliken pop alpha/beta on meta-lowdin orthogonal AOs**
** Mulliken atomic charges ( Nelec_alpha| Nelec_beta) **
charge of 0H= 0.00000( 0.500000.50000)
charge of 1H= -0.00000( 0.500000.50000)
charge of 2H= 0.00000( 0.500000.50000)
charge of 3H= -0.00000( 0.500000.50000)
Dipole moment(X, Y, Z, Debye): 0.00000, 0.00000, 0.00000
Anti-ferromamagnetic initial guess:
converged SCF energy = -0.615622616250374<S^2> = 1.0086432S+1= 2.2437852
**** MO energy ****
alpha | beta alpha | beta
MO #1 energy= -1.19996932121733 | -1.19996932121733 occ= 1 | 1
MO #2 energy= -0.192312470583276 | -0.192312470583277 occ= 1 | 1
MO #3 energy= 0.510790949191351 | 0.510790949191353 occ= 0 | 0
MO #4 energy= 2.38286361752655 | 2.38286361752655 occ= 0 | 0
** Mulliken pop alpha/beta on meta-lowdin orthogonal AOs**
** Mulliken atomic charges ( Nelec_alpha| Nelec_beta) **
charge of 0H= -0.00000( 0.773240.22676)
charge of 1H= -0.00000( 0.226760.77324)
charge of 2H= -0.00000( 0.773240.22676)
charge of 3H= 0.00000( 0.226760.77324)
Dipole moment(X, Y, Z, Debye): 0.00000, 0.00000, -0.00000
尽管构造密度矩阵初始猜测看起来比较麻烦，但这给平均场初始猜测的构造提供了最大的自由度。比如有一些计算中，我们希望对调 HOMO LUMO 轨道，这时可以很容易地通过修改轨道占据数 (mf.mo_occ) 或轨道系数 (mf.mo_coeff) 后计算密度矩阵实现。
from pyscf import gto, scf
mol = gto.M(atom='''
H 0 0 0 ; H 0 0 0.5;
H 0 0.5 0.5; H 0 0.5 0 ''', basis='sto3g')
mf = scf.UHF(mol)
mf.kernel()Exchange the occupancies of alpha HOMO and alpha LUMO
mf.mo_occ[0,1] = 0
mf.mo_occ[0,2] = 1
dm = mf.make_rdm1(mf.mo_coeff, mf.mo_occ)
mf.kernel(dm)
此外，PySCF 支持从 checkpoint 文件读取平均场计算的信息，然后根据当前计算的特征构造一个密度矩阵初始猜测。比如在一个先导的平均场计算里把平均场的信息保存在 checkpoint 文件里，
from pyscf import gto, scf
mol = gto.M(atom='H 0 0 0; F 0 0 2.5', basis='ccpvdz')
mf = scf.HF(mol)
mf.chkfile = 'hf.chk'
mf.mac_cycle = 2
mf.kernel()然后在另一个计算中读取 checkpoint 文件获得前面的平均场的信息并构造密度矩阵作为初始猜测
from pyscf import gto, scf
mol = gto.M(atom='H 0 0 0; F 0 0 2.5', basis='ccpvdz')
mf = scf.HF(mol)
dm = mf.from_chk('hf.chk')
mf.kernel(dm)先导计算里的分子结构和基组不必要和后续计算一样，PySCF 程序能自动识别结构和基组的差异，正确地把先导计算的结果投影到当前的计算里。通过这个特性，我们可以先对体系做一个小基组的计算来获得初始猜测，然后把初始猜测投影给大基组，从而减少大基组下计算的时间和不确定性。
更多初始猜测的例子参见 PySCF 源代码提供的例子 https://github.com/pyscf/pyscf/blob/master/examples/scf/14-restart.py 和 https://github.com/pyscf/pyscf/blob/master/examples/scf/15-initial_guess.py
二阶收敛方法
PySCF 提供了二阶收敛方法。在处理 DIIS 方法难以收敛的开壳层体系，非平衡结构体系时表现稳定，能够高效地收敛到与初始猜测非常接近的体系局域极小点。
相对 DIIS 方法，PySCF 二阶算法有如下优点：
1. DIIS 迭代方法容易卡在鞍点。二阶算法计算了能量的二阶导，有效地避免了鞍点问题。
2. 二阶算法直接优化能量梯度，不受 aufbau principle 的限制，对开壳层体系容易找到能量更低的态，即便这个态对应的轨道能违反了 aufbau principle。
3. 二阶算法连续地从初始态变换到终态，不会出现占据态和空态翻转引起体系波函数的剧烈变化，容易通过初试猜测来控制收敛的结果。
实际经验显示，碰到 DIIS 较难收敛的过渡金属化合物时，PySCF 二阶算法表现良好，相对 DIIS 算法经常能找到能量更低的态。当体系比较简单，比如平衡结构的闭壳层计算，DIIS 方法一般收敛迅速，二阶算法与 DIIS 表现相近，没有明显的优势。
PySCF 二阶算法可以通过 .newton() 方法调用。例如下面这个计算，使用 DIIS 算法很难收敛这个体系
from pyscf import gto, scf, dft
mol = gto.M(atom='''
Fe -3.301 -0.904 0.
N -3.078 0.552 0.
N -2.377 0.160 0.
''', basis='ccpvdz', verbose=4)
mf = dft.RKS(mol)
mf.xc = 'pbe'
mf.conv_tol = 1e-6
mf.level_shift = .2
mf.max_cycle = 200
mf.run()一些常规的设置，比如调整 level shift，DIIS space，或增加迭代的步数等都无法将体系的能量收敛到 1e-6 的精度。
...
HOMO = -0.128165038749955 LUMO = 0.0413831518440274
cycle= 198 E= -1371.88589574082 delta_E= 0.000406|g|= 0.123|ddm|= 0.0396
HOMO = -0.127982085807919 LUMO = 0.0414632150014509
cycle= 199 E= -1371.88540959507 delta_E= 0.000486|g|= 0.133|ddm|= 0.0422
HOMO = -0.127905703677111 LUMO = 0.0417282896132888
cycle= 200 E= -1371.88480160486 delta_E= 0.000608|g|= 0.144|ddm|= 0.0465
SCF not converged.当切换至二阶算法后
from pyscf import gto, scf, dft
mol = gto.M(atom='''
Fe -3.301 -0.904 0.
N -3.078 0.552 0.
N -2.377 0.160 0.
''', basis='ccpvdz', verbose=4)
mf = dft.RKS(mol)
mf.xc = 'pbe'
mf = mf.newton()
mf.run()程序很快就收敛到一个解，能量比上面未收敛的 DIIS 的解更低
...
macro= 2 E= -1371.88681267525 delta_E= -0.0754391|g|= 0.1776653 KF 13 JK
macro= 3 E= -1371.88884629902 delta_E= -0.00203362|g|= 0.001372923 KF 13 JK
macro= 4 E= -1371.88884629986 delta_E= -8.43329e-10|g|= 0.0002635921 KF 2 JK
Canonicalize SCF orbitals
macro X = 5 E=-1371.88884629986|g|= 0.000263592 total 14 KF 65 JK
converged SCF energy = -1371.88884629986如果一个体系 DIIS 迭代不收敛，很多时候我们可以把 DIIS 未收敛的结果作为初始猜测输入到二阶算法来获得一个解
from pyscf import gto, scf, dft
mol = gto.M(atom='''
Fe 3.301 -0.904 0.
N 3.078 0.552 0.
N 2.377 0.160 0.
''',
basis='ccpvdz', verbose=4, spin=0)
mf = dft.RKS(mol)
mf.xc = 'pbe'
mf.conv_tol = 1e-6
mf.level_shift = .2Save the DIIS results in the specified checkpoint file.
mf.chkfile = 'diis-iter.chk'
mf.run()Create a new DFT calculation
mf = dft.RKS(mol)
mf.xc = 'pbe'
mf.conv_tol = 1e-6Load the previous DIIS results into the new DFT calculation.
mf.dict.update(scf.chkfile.load('diis-iter.chk', 'scf'))
mf = mf.newton()
mf.run()
在这个例子里，我们需要把第一步的 DIIS 结果保存到一个 checkpoint 文件里以供另一个计算加载。在第二步使用二阶算法时，我们通过 scf.chkfile.load 这个函数读取 checkpoint 文件里的平均场结果，并通过 .dict.update() 这个方法把平均场结果载入到二阶方法的对象里。然后这个平均场结果会被用作后续二阶算法的初始猜测。
PySCF 二阶算法有很多参数可以调整。合适的参数设置可以有效地提高计算效率。如果一个计算对初始猜测没有特殊要求的话，PySCF 程序包提供了一种可能的配置，综合利用了各种优化技术，对很多体系都可以实现相对默认二阶方法 3 - 5 倍的速度提升。这个优化的配置可以通过函数 scf.fast_newton 来调用。
from pyscf import gto, scf, dft
mol = gto.M(atom='''
Fe -3.301 -0.904 0.
N -3.078 0.552 0.
N -2.377 0.160 0.
''', basis='ccpvdz', verbose=4)
mf = dft.RKS(mol)
mf.xc = 'pbe'
mf.conv_tol = 1e-6
scf.fast_newton(mf)