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.py

PySCF 默认的收敛判别标准非常严格，判别标准是能量收敛到 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.1

level shift 的单位是 hartree 。这个值一般不用特别大，0.2 左右对大多数体系应该就够用了。当体系的基组有较大的线性相关时，过大的 level shift 不利于数值稳定性。

3. 我们可以调整 DIIS 空间的大小来给予 DIIS 算法更大的自由度来搜索解，比如

mf.diis_space = 12

DIIS 的收敛还可以被一系列别的参数影响，更多的设定可以参考帮助文档 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 = .2
# Save 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-6
# Load 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)