PySCF程序包平均场计算的一些收敛技巧

PySCF程序包平均场计算的一些收敛技巧

平均场计算是 PySCF 程序包里优化得比较并全面的模块之一。在平均场模块里,PySCF支持 RHF, UHF, ROHF, GHF, RKS, UKS, ROKS, GKS 等一系列方法来研究闭壳层体系、开壳层体系、复数哈密顿量体系、相对论效应、溶剂化效应。同时 PySCF 提供了大量的辅助功能来帮助平均场计算收敛。以下我们通过一些例子来演示在 PySCF 里收敛平均场计算的技巧。

以下的例子在 PySCF-1.5 以上发行版均可使用。

首先我们使用 PySCF 程序包实现一个最简单的平均场计算。在编辑器里创建以下 Python 代码

代码语言:javascript
复制
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 脚本就可以做一个最简单的平均场计算

代码语言:javascript
复制
python example.py

PySCF 默认的收敛判别标准非常严格,判别标准是能量收敛到 1e-10 Eh,能量梯度的 norm 收敛到 1e-5 。收敛参数可以通过设定平均场对象的 .conv_tol, .conv_tol_grad 来调整。大多数时候,设定

代码语言:javascript
复制
mf.conv_tol = 1e-6

的收敛标准已经足够化学问题的研究了。

DIIS 方法

DIIS 迭代方法是 PySCF 里默认的平均场收敛算法。在上面的例子里,DIIS 算法进行了 7 轮迭代得到了收敛结果。

当体系比较复杂时,DIIS 算法的默认参数不一定能得到收敛的平均场结果。此时,有一些常规的参数可以帮助 DIIS 算法收敛,比如

1. 我们可以增加迭代步数

代码语言:javascript
复制
mf.max_cycle = 200

来等待 DIIS 的结果慢慢移动到收敛的结果上。

2. 我们可以调整 level shift 的大小来改变体系 HOMO-LUMO gap。较大的 HOMO-LUMO gap 一般能避免收敛过程的震荡。

代码语言:javascript
复制
mf.level_shift = 0.1

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

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

代码语言:javascript
复制
mf.diis_space = 12

DIIS 的收敛还可以被一系列别的参数影响,更多的设定可以参考帮助文档 http://pyscf.org/pyscf/scf.html#hartree-fock

初始猜测

在平均场计算中,很多困难体系对初始猜测十分敏感,合适的初始猜测可以有效地帮助平均场收敛。PySCF 的平均场默认的初始猜测是基于 ANO-RCC 基组构造的原子密度矩阵的叠加。这个初始猜测对平衡结构,或较复杂的问题(如开壳层,非平衡结构,带有少量电荷的体系)都有不错的表现。除了 ANO 作为初始猜测以外,PySCF 还提供了一系列方法对平均场计算的初始猜测进行调整。

有一些体系需要特殊的初始猜测才能收敛到正确的态上,比如铁磁或反铁磁的初始猜测。此时,就需要对 PySCF 的初始猜测进行调整。PySCF 支持构造好一个密度矩阵,然后传给平均场对象作为初始猜测进行计算。例如

代码语言:javascript
复制
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()

可以看到默认的初始猜测收敛到了能量较高的态上,而反铁磁的初始猜测收敛到能量较低的态上。

代码语言:javascript
复制
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) 后计算密度矩阵实现。

代码语言:javascript
复制
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 文件里,

代码语言:javascript
复制
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 文件获得前面的平均场的信息并构造密度矩阵作为初始猜测

代码语言:javascript
复制
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.pyhttps://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 算法很难收敛这个体系

代码语言:javascript
复制
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 的精度。

代码语言:javascript
复制
...
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.

当切换至二阶算法后

代码语言:javascript
复制
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 的解更低

代码语言:javascript
复制
...
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 未收敛的结果作为初始猜测输入到二阶算法来获得一个解

代码语言:javascript
复制
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 来调用。

代码语言:javascript
复制
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)