|
|
这篇文章主要是起一个抛砖引玉的作用，讨论比较浅显。本文中仅涉及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基组之间。 |
评分 Rate
-
查看全部评分 View all ratings
|
发表于 Post on 2022-4-30 16:16:23
收藏 Add to favorites
楼主 Author