梯度的基组外推以及在结构优化中的实现

ionexchangeC

这篇文章主要是起一个抛砖引玉的作用，讨论比较浅显。本文中仅涉及HF/DFT的外推，后HF相关能的外推原理相似。
一、梯度外推公式推导
在社长的博文http://sobereva.com/172中，提到HF/DFT能量的外推公式如下：
E_SCF(L)=E_SCF(∞)+A*exp(-α*√L)
上式中一共有三个未知数：E_SCF(∞)、A和α，从原理上来讲它们都应当是分子坐标xyz的函数，但是其中α随分子坐标的变化幅度明显比另外两个小，所以不妨视它不随分子坐标变化而变化，即设对于分子中的每个原子来说都是dα/dx=dα/dy=dα/dz=0。这样不仅可以简化梯度的外推时的公式推导，而且还可以直接使用能量外推时拟合好的α值，这样就可以两点外推梯度了。
按博文中的做法消掉A，得：
E_SCF(∞) = [ E_SCF(M)*exp(-α*√N) - E_SCF(N)*exp(-α*√M) ] / [exp(-α*√N) - exp(-α*√M)]
等式两边同时对分子坐标求导（用x代表至少3N-6个的自变量）：
d(E_SCF(∞))/dx= [ d(E_SCF(M))/dx*exp(-α*√N) - d(E_SCF(N))/dx*exp(-α*√M) ] / [exp(-α*√N) - exp(-α*√M)]
这样就获得了梯度的外推公式。
二、Gaussian中的实现
要想在Gaussian中实现用外推的梯度优化，要用到external功能。这里我参考以下两篇文章：
http://bbs.keinsci.com/thread-10141-1-1.html
http://bbs.keinsci.com/thread-18901-1-1.html
写了.sh脚本和fortran程序来实现，实现的是def2-SVP/def2-TZVPP这两个基组外推CBS。本人计算机水平较差，修改代码时可能会留下一些不必要的语句，请大家见谅。三个文件在附件中给出，要是可执行程序用不了，就自行编译f90文件。泛函默认的是PBE0（写作PBE1PBE），想改的话直接改.sh文件就行。
三、补充
当然，用相同的思路还可以进行基组外推下的频率计算。也许以后会补充上。
对于一般的泛函来说，3-zeta基组下结构优化的结果已经很理想了，在基组外推下优化实际意义并不大。
想用三点外推来外推梯度的话，也是可以的，不过需要解三元方程，很麻烦，而且似乎不好消元。
使用该脚本时时输入文件写#p opt=nomicro external='./gau_gau.sh'
四、测试
测试2个体系：氟化氢和氟甲烷，泛函都是PBE0。测试结果如下：
氟化氢
基组        键长(埃)
def2-SVP        0.9202
def2-TZVPP        0.9178
def2-QZVPP        0.9176
外推CBS        0.9178
氟甲烷
优化前坐标：
C                 -1.13207554   -0.58780841    0.00000000
H                 -0.77542111   -1.59661841    0.00000000
H                 -0.77540270   -0.08341022    0.87365150
H                 -0.77540270   -0.08341022   -0.87365150
F                 -2.48207554   -0.58779178    0.00000000
基组        C-H键长        C-F键长        H-C-H键角        H-C-F键角        2-1-4-3二面角
def2-SVP        1.1028        1.3636        108.9295        110.0076        -118.6969
def2-TZVPP        1.0917        1.3759        109.8308        109.1092        -120.8918
def2-QZVPP        1.0911        1.3764        109.892        109.0471        -121.0457
外推CBS        1.0913        1.3764        109.8673        109.0722        -120.9835
观察数据发现，外推CBS的优化精度似乎介于3-zeta和4-zeta基组之间。

xizaokaiz

想请教一下，代码里定义的alpha=1.039是怎么来的呀，如果基组不同是不是要重新选择参数

ionexchangeC

2023.4.29更新：支持了计算解析频率 extrapolate(Freq).zip (7.71 KB, 下载次数 Times of downloads: 6)

2023-4-29 11:19 上传 Uploaded

点击下载
Click to download

ldatea

我记得很清楚，第一次看到sob老师那个Powell的帖子，第一个想到的就是优化基函数。但是水平非常有限，就打算拿H的基组试试，还在群里问sob老师那个f90文件怎么用。分别用用H2的能量和键长作为优化目标，结果发现很难搞。碰了一堆钉子就搁置了，没想到居然在这里能见到一样想法的人。

ionexchangeC

先用数值频率凑合了下，测试了一下氟甲烷的振动分析结果：
基组        频率1        频率2-3        频率4        频率5-6        频率7        频率8-9
def2-SVP        1140.37        1200.43        1493.14        1468.82        3038.03        3138.14
def2-TZVPP        1097.26        1198.78        1485.49        1489.73        3042.39        3129.65
def2-QZVPP        1091.44        1199.39        1483.27        1490.02        3041.11        3128.43
外推CBS        1095.50        1198.75        1485.24        1490.56        3042.56        3129.27

chrinide

对DFT意义确实不大，一般而言3z-4z基本就收敛了，用在波函数上，比如搞个中小标准测试集的基准的CCSD(T)-CBS几何结构数据库，看看疗效