专栏名称: 研之成理
夯实基础，让基础成就辉煌；传递思想，让思想改变世界。“研之成理科研平台”立足于科研基础知识与科研思想的传递与交流，旨在创建属于大家的科研乐园！主要内容包括文献赏析，资料分享，科研总结，论文写作，软件使用等。科研路漫漫，我们会一路陪伴你！
目录
相关文章推荐
国际新闻界  ·  2024年中国的新闻学研究 ·  昨天  
募格学术  ·  DeepSeek王炸登场！AI的神级应用被发 ... ·  昨天  
实验万事屋  ·  我一直觉得博士生发这10.7分的SCI期刊就 ... ·  昨天  
科研圈  ·  特朗普重拳打击美国高校：腰斩研究管理费，引发 ... ·  2 天前  
募格学术  ·  DeepSeek评中国最好的十大高校，和Ch ... ·  3 天前  
51好读  ›  专栏  ›  研之成理

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

研之成理  · 公众号  · 科研  · 2019-12-11 07:00

正文

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

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

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

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


  1. benzene = '''


  2. C 4.673795     6.280948    0.


  3. C 5.901190     5.572311    0.


  4. C  5.901190     4.155037    0.


  5. C  4.673795     3.446400    0.


  6. C  3.446400     4.155037    0.


  7. C  3.446400     5.572311    0.


  8. H 4.673795     7.376888    0.


  9. H 6.850301     6.120281    0.


  10. H  6.850301     3.607068    0.


  11. H  4.673795     2.350461    0.


  12. H 2.497289     3.607068    0.


  13. H  2.497289     6.120281    0.


  14. '''



  16. frompyscfimportgto, scf


  17. # Createan moleculeobject


  18. mol= gto.M(atom=benzene, basis='ccpvdz', verbose=4)


  19. # Createan mean-field object


  20. mf= scf.RHF(mol)




    



  21. # Solve the SCF problem


  22. mf.run()

然后在命令行里执行这个 Python 脚本就可以做一个最简单的平均场计算 


  1. python example.py

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


  1. mf.conv_tol = 1e-6

DIIS 方法

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

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

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


  1. mf.max_cycle = 200

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

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


  1. mf.level_shift = 0.1

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

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


  1. mf.diis_space = 12

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

初始猜测

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

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


  1. importnumpyas np


  2. frompyscfimportgto, scf


  3. mol =gto.M(atom='''


  4. H 0 0   0 ;H 0 0   0.5;


  5. H 0 0.5 0.5; H 0 0.5 0  ''', basis='sto3g')


  6. print('Default initial guess:')


  7. mf= scf.UHF(mol)


  8. mf.kernel()


  9. mf.analyze()



  11. print('\nAnti-ferromamagnetic initial guess:')


  12. dm =np.zeros((2,4,4))


  13. dm[0,0,0] = dm[0,2,2] = 1


  14. dm[1,1,1] = dm[1,3,3] = 1


  15. mf.kernel(dm)


  16. mf.analyze()

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


  1. Default initial guess:


  2. converged SCF energy = -0.521792378127812 = 4.6629367e-152S+1= 1


  3. **** MO energy ****


  4.                              alpha | beta                alpha | beta


  5. MO #1   energy= -1.19356543325823  | -1.19356543325823  occ= 1 | 1


  6. MO #2   energy= -0.0986706198268508 | -0.09867061982685 occ= 1 | 1


  7. MO #3   energy= 0.416828469148984  | 0.416828469148984 occ= 0 | 0


  8. MO #4   energy= 2.37630222058633   | 2.37630222058633   occ= 0 | 0


  9. ** Mulliken pop alpha/beta on meta-lowdinorthogonal AOs**


  10. ** Mulliken atomic charges( Nelec_alpha|Nelec_beta) **


  11. charge of 0H=      0.00000( 0.500000.50000)


  12. charge of  1H=     -0.00000( 0.500000.50000)


  13. charge of  2H=      0.00000( 0.500000.50000)


  14. charge of  3H=     -0.00000(     0.500000.50000)


  15. Dipole moment(X, Y, Z, Debye):  0.00000,  0.00000,  0.00000



  17. Anti-ferromamagnetic initial guess:


  18. converged SCF energy = -0.615622616250374 = 1.0086432S+1= 2.2437852


  19. **** MO energy ****


  20.                              alpha | beta                alpha | beta


  21. MO #1   energy= -1.19996932121733  | -1.19996932121733  occ= 1 | 1


  22. MO #2   energy= -0.192312470583276 | -0.192312470583277 occ= 1 | 1


  23. MO #3   energy= 0.510790949191351  | 0.510790949191353  occ= 0 | 0




    



  24. MO #4   energy= 2.38286361752655   | 2.38286361752655   occ= 0 | 0


  25. ** Mulliken pop alpha/beta on meta-lowdinorthogonal AOs**


  26. ** Mulliken atomic charges( Nelec_alpha|Nelec_beta) **


  27. charge of 0H=     -0.00000( 0.773240.22676)


  28. charge of 1H=     -0.00000( 0.226760.77324)


  29. charge of 2H=     -0.00000( 0.773240.22676)


  30. charge of 3H=      0.00000( 0.226760.77324)


  31. Dipole moment(X, Y, Z, Debye):  0.00000,  0.00000, -0.00000

尽管构造密度矩阵初始猜测看起来比较麻烦，但这给平均场初始猜测的构造提供了最大的自由度。 比如有一些计算中，我们希望对调 HOMO LUMO 轨道，这时可以很容易地通过修改轨道占据数 (mf.mo_occ) 或轨道系数 (mf.mo_coeff) 后计算密度矩阵实现。 


  1. from pyscf import gto, scf


  2. mol =gto.M(atom='''


  3. H 0 0   0  ; H 0 0   0.5;


  4. H 0 0.5 0.5; H 0 0.5 0  ''', basis='sto3g')


  5. mf = scf.UHF(mol)


  6. mf.kernel()



  8. # Exchange the occupancies of alpha HOMO and alpha LUMO


  9. mf.mo_occ[0,1] = 0


  10. mf.mo_occ[0,2] = 1


  11. dm = mf.make_rdm1(mf.mo_coeff, mf.mo_occ)


  12. mf.kernel(dm)

此外，PySCF 支持从 checkpoint 文件读取平均场计算的信息，然后根据当前计算的特征构造一个密度矩阵初始猜测。 比如在一个先导的平均场计算里把平均场的信息保存在 checkpoint 文件里 


  1. from pyscf import gto, scf


  2. mol = gto.M(atom='H 0 0 0; F 0 0 2.5', basis='ccpvdz')


  3. mf = scf.HF(mol)


  4. mf.chkfile = 'hf.chk'


  5. mf.mac_cycle = 2


  6. mf.kernel()

然后在另一个计算中读取 checkpoint 文件获得前面的平均场的信息并构造密度矩阵作为初始猜测 


  1. from pyscf import gto, scf


  2. mol = gto.M(atom='H 0 0 0; F 0 0 2.5', basis='ccpvdz')


  3. mf = scf.HF(mol)


  4. dm = mf.from_chk('hf.chk')


  5. 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 算法很难收敛这个体系 


  1. frompyscfimportgto, scf,dft




    



  2. mol =gto.M(atom='''


  3. Fe  -3.301 -0.904 0.


  4. N   -3.078 0.552 0.


  5. N   -2.377  0.160  0.


  6. ''', basis='ccpvdz', verbose=4)


  7. mf = dft.RKS(mol)


  8. mf.xc = 'pbe'


  9. mf.conv_tol = 1e-6


  10. mf.level_shift = .2


  11. mf.max_cycle = 200


  12. mf.run()

一些常规的设置，比如调整 level shift，DIIS space，或增加迭代的步数等都无法将体系的能量收敛到 1e -6 的精度。 


  1. ...


  2.   HOMO = -0.128165038749955  LUMO = 0.0413831518440274


  3. cycle= 198 E= -1371.88589574082  delta_E= 0.000406|g|= 0.123|ddm|= 0.0396


  4.   HOMO = -0.127982085807919  LUMO = 0.0414632150014509


  5. cycle= 199 E= -1371.88540959507  delta_E= 0.000486|g|= 0.133|ddm|= 0.0422


  6.   HOMO = -0.127905703677111  LUMO = 0.0417282896132888


  7. cycle= 200 E= -1371.88480160486  delta_E= 0.000608|g|= 0.144|ddm|= 0.0465


  8. SCF not converged.

当切换至二阶算法后 


  1. from pyscf import gto, scf, dft


  2. mol = gto.M(atom='''


  3. Fe  -3.301 -0.904  0.


  4. N   -3.078  0.552  0.


  5. N   -2.377  0.160  0.


  6. ''', basis='ccpvdz', verbose=4)


  7. mf = dft.RKS(mol)


  8. mf.xc = 'pbe'


  9. mf = mf.newton()


  10. mf.run()






请到「今天看啥」查看全文


推荐文章
国际新闻界  ·  2024年中国的新闻学研究
昨天
募格学术  ·  DeepSeek王炸登场！AI的神级应用被发现，一天时间帮你提升科研技能与实战应用能力
昨天
实验万事屋  ·  我一直觉得博士生发这10.7分的SCI期刊就是水刊！但看了这四川大学中性粒细胞氨基酸代谢的SCI文章，竟然意外地不错……
昨天
科研圈  ·  特朗普重拳打击美国高校：腰斩研究管理费，引发大学公愤；谷歌最新AI超越奥数金牌选手平均水平 | 科研圈日报
2 天前
募格学术  ·  DeepSeek评中国最好的十大高校，和ChatGPT给出的答案大不同
3 天前
国家人文历史  ·  被诅咒的乱世黄金：那些拿到沙皇黄金的，没一个有好下场
7 年前
读者  ·  夜听 | 愿你阳光下像个孩子，风雨里像个大人
7 年前
班主任家园  ·  高考成绩出来了，请不要打扰别人的幸福！
7 年前
电影铺子  ·  这才是最厉害的小三（不要再骂吴越了）
7 年前
彬彬有理  ·  遇到一个差不多先生，你嫁吗？
7 年前
Sov5搜索   ·   小百科   ·   今天看啥   ·   移动版
51好读 - 好文章就要读起来!