计算化学公社

标题: 使用NEWUOA局部极小化算法对基组的指数进行优化 [打印本页]

作者
Author: ionexchangeC 时间: 2022-6-11 02:18
标题: 使用NEWUOA局部极小化算法对基组的指数进行优化
本帖最后由 ionexchangeC 于 2022-6-12 21:52 编辑

量子化学计算中基组的重要性不言而喻。虽然有很多现成的，已经搭配好的基组可供我们在计算中选用，不过总是会有现有基组不够用，不能满足对当前研究的需要的情况。遇到这种情况的解决办法之一就是自行为基组增加一些基函数。这些基函数可以由别的等级相似的基组中直接借用，但有一个问题就是这样移植来的基函数，指数并不一定是当前基组下最优的指数。如果对计算的精度要求高，可以尝试自行对基组的指数进行优化。
对基组的指数进行优化，可以从各种属性入手，例如有为算准极化率而优化的def2-nZVPD系列基组，为算准核磁共振磁屏蔽值的pcSseg系列基组，不过最常见策略的还是以能量计算为标准优化基组。一个非常简单的思路是：使用能量具有变分原理的方法（HF、CI等）进行能量计算，据此优化基函数的指数使得体系能量最低。这个体系可以采用最小的体系，也即一个原子。（当然如pc系列基组则是以一些小分子的能量计算结果来优化极化函数的指数）
为了优化能量随着基函数的指数变化的极小值点，需要一些优化算法。最容易想的就是牛顿法/赝牛顿法，但是这需要体系能量对基函数指数的导数，以我现有的水平显然是不知道它的解析表达式的，基于有限差分的话，要优化的基函数一多又显得特别麻烦。为此，我想到了sob老师的一篇博文《无需导数的局部极小化算法NEWUOA在Fortran中的使用简例》http://sobereva.com/536，决定用这种算法优化基函数指数。目前写的程序还非常死板，要用不同体系去优化不同的基组还得直接改源代码。这个程序我命名为basissetoptimizer，调用的程序是Gaussian16。核心的源代码如下：
PROGRAM test_newuoa
use newuoa_module
implicit real*8 (a-h,o-z)
real*8,allocatable :: X(:)
external CALFUN

nvar=3 !被优化的变量数量
allocate(X(nvar))
X(:)=0
X(1)=1.2
X(2)=0.4
X(3)=1.0 !基函数指数的初猜，要注意的是这种算法对初猜的要求还是比较高的
call newuoa_min(CALFUN, X, RHOBEG=0.1D0, RHOEND=1D-6, IPRINT=2, MAXFUN=50000) !可调节的参数
END PROGRAM

subroutine CALFUN(X,F)
real*8 :: X(:),F
real*8 ene
open(168,FILE="./1.gjf") !创建Gaussian输入文件
write(168,*) "%chk=1.chk" !以下均是输入文件的内容，文件名和计算配置可以自行修改
write(168,*) "%mem=6GB"
write(168,*) "%nproc=4"
write(168,*) "#p CISD/gen IOp(7/32=3) Force 5D 7F GFInput" !关键词，5D 7F是为了方便能量比较，GFInput输出基组信息。加Force关键词是为了能用IOp(7/32=3)以确定的格式输出能量
write(168,*) ""
write(168,*) "Title"
write(168,*) ""
write(168,*) "0 3"
write(168,*) "C 0. 0. 0."
write(168,*) ""
write(168,*) "C    0"
write(168,*) "6-311G"
write(168,*) "D 1 1.00"
write(168,*) "    ",X(1),"    1.0000000" !被优化的基函数指数，不希望被优化的话直接写数字，对于收缩型基函数，收缩系数也可以被优化
write(168,*) "D 1 1.00"
write(168,*) "    ",X(2),"    1.0000000"
write(168,*) "F 1 1.00"
write(168,*) "    ",X(3),"    1.0000000"
write(168,*) "****"
write(168,*) ""
write(168,*) ""
close(168)
call system('g16 <1.gjf> 1.out') !运行Gaussian
call system('./1.sh') !提取数据，1.sh在压缩包内
open(10,file="2.out",status="old")
read(10,*) ene
print "(D20.12)",ene
close(10)
call system('rm PES.out')
call system('rm 2.out')
F=ene
end subroutine
基于此程序，我进行了四个基组优化作为使用例。优化的过程详见压缩包里的result文本文件。
例一：为C的6-311G基组加入2个D极化函数和1个F极化函数，使得C原子在CISD理论方法下计算出的能量为最低。
优化结果是2个D极化函数的指数分别是1.061023D+00和3.155438D-01，F极化函数的指数是7.723239D-01。在Gaussian中的6-311G(2df)对于C来说，2个D极化函数的指数分别是0.1252000000D+01和0.3130000000D+00，F极化函数的指数是0.8000000000D+00，存在一定差异。二者分别在CISD方法下计算的能量为-37.7754986和-37.7751304，在自行优化的基组下计算能量确实更低些。

例二：为Br的lanl08基组加入2个D极化函数和1个F极化函数，使得Br原子在CISD理论方法下计算出的能量为最低。
本例注意对于Br来说lanl08是赝势基组，不要忘记输入文件末尾写赝势。本例中初猜比较重要，这里我给了D极化函数指数0.4、0.15，F极化函数指数0.35。优化后的结果为：
D极化函数指数：6.262004D-01 2.466382D-01
F极化函数指数：5.693046D-01

例三：为S的def2-TZVP基组加入1个描述核电子的大指数S型基函数和1个描述核电子的大指数P型基函数，使得S原子在HF理论方法下计算出的能量为局部极低。
本例中体系能量在优化过程中变化不大，最终优化的结果为S基函数指数为9.208457D+00，P基函数指数为6.918190D+00。

例四：对F的月份基组jun-cc-pVTZ的三个弥散函数的指数进行优化，使得氟离子在CISD理论方法下计算出的能量最低。
本例的初猜直接选用aug-cc-pVTZ的SPD角动量弥散函数的指数即可。最终的结果是：
优化后的jun-cc-pVTZ基组算F-：
S基函数指数0.1177508704D+00
P基函数指数0.7310888554D-01
D基函数指数0.2913912155D+00
能量-99.7234903
未优化的jun-cc-pVTZ基组算F-：
S基函数指数0.9158000000D-01
P基函数指数0.7361000000D-01
D基函数指数0.2920000000D+00
能量-99.7234509

作者
Author: wzkchem5 时间: 2022-6-11 03:24
我最近一段时间也在做一些derivative free优化问题，用的是simplex方法。不知道有没有什么文章比较过simplex和NEWUOA哪个好。
另外最近我发现，当被优化函数可以表示为很多函数的平方和f(x)=f1(x)^2+...+fn(x)^2的时候（比如优化一个泛函的参数，让泛函在训练集上的计算结果尽量接近高精度理论计算值，f(x)可以取为训练集所有分子的误差的平方和），可以用DIIS加速simplex方法，具体来说就是以(f1(x),...,fn(x))作为DIIS的error vector，如果某一步DIIS失败就改成做一步simplex。效果相比simplex有极其显著的提升，50个变量的函数跑个200步就差不多收敛了。但是这个方法是利用了f(x)的结构，所以不适用于f(x)是黑箱子的情形（比如SCF能量）。
不过其实也可以思考一下，量化软件其实不止打印了SCF能量，还有单电子项、双电子项、XC项等等，利用这些信息不知道能不能加速收敛（乃至于可以做一个能量分解分析，把能量分解成各个AO、AO对的贡献，这样就有更多项了）。换句话说已知f(x)=f1(x)+...+fn(x)，但是各项有正有负，能不能设计出一种算法，相比只知道f(x)的值的情况下收敛得更快

作者
Author: chrinide 时间: 2022-6-11 10:31
本帖最后由 chrinide 于 2022-6-11 12:43 编辑

wzkchem5 发表于 2022-6-11 03:24
我最近一段时间也在做一些derivative free优化问题，用的是simplex方法。不知道有没有什么文章比较过simple ...

可以试试这个包 PDFO (Powell's Derivative-Free Optimization solvers) package. https://www.pdfo.net 主流的derivative free优化算法都包含了

作者
Author: ionexchangeC 时间: 2022-6-13 21:06
做了一些修改，省掉了1.sh和无意义地解一次CPSCF方程
(, 下载次数 Times of downloads: 5)

作者
Author: ionexchangeC 时间: 2022-7-14 20:42
7.14补充：
部分情况下，使用该种算法时会导致基组的指数被优化成负数，之后导致任务失败。此时有2个解决方法：
1.RHOBEG设小一点，避免优化时一步的步长太大走到负数；
2.将输入文件段的X(1)改成诸如Y(1)，令Y(1)=X(1)^2，保证了非负性。

作者
Author: Warm_Cloud 时间: 2022-7-15 07:56

ionexchangeC 发表于 2022-7-14 20:42
7.14补充：
部分情况下，使用该种算法时会导致基组的指数被优化成负数，之后导致任务失败。此时有2个解决 ...

指数为负的问题，我直接用的abs(X(1))

欢迎光临计算化学公社 (http://bbs.keinsci.com/)