|
|
本帖最后由 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
|
评分 Rate
-
查看全部评分 View all ratings
|
发表于 Post on 2022-6-11 02:18:44
收藏 Add to favorites
楼主 Author