计算化学公社

 找回密码 Forget password
 注册 Register
Views: 1413|回复 Reply: 3

[量化理论] 梯度的基组外推以及在结构优化中的实现

[复制链接 Copy URL]

133

帖子

2

威望

801

eV
积分
974

Level 4 (黑子)

发表于 Post on 2022-4-30 16:16:23 | 显示全部楼层 Show all |阅读模式 Reading model
这篇文章主要是起一个抛砖引玉的作用,讨论比较浅显。本文中仅涉及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基组之间。

extrapolate.zip

6.26 KB, 下载次数 Times of downloads: 9

评分 Rate

参与人数
Participants 4
eV +19 收起 理由
Reason
丁越 + 5 赞!
sobereva + 8
wzkchem5 + 3
chrinide + 3 赞!

查看全部评分 View all ratings

272

帖子

0

威望

3941

eV
积分
4213

Level 6 (一方通行)

发表于 Post on 2022-4-30 17:31:19 | 显示全部楼层 Show all
对DFT意义确实不大,一般而言3z-4z基本就收敛了,用在波函数上,比如搞个中小标准测试集的基准的CCSD(T)-CBS几何结构数据库,看看疗效

133

帖子

2

威望

801

eV
积分
974

Level 4 (黑子)

 楼主 Author| 发表于 Post on 2022-5-1 00:48:04 | 显示全部楼层 Show all
先用数值频率凑合了下,测试了一下氟甲烷的振动分析结果:
基组        频率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

112

帖子

1

威望

933

eV
积分
1065

Level 4 (黑子)

发表于 Post on 2022-7-5 20:37:02 | 显示全部楼层 Show all
我记得很清楚,第一次看到sob老师那个Powell的帖子,第一个想到的就是优化基函数。但是水平非常有限,就打算拿H的基组试试,还在群里问sob老师那个f90文件怎么用。分别用用H2的能量和键长作为优化目标,结果发现很难搞。碰了一堆钉子就搁置了,没想到居然在这里能见到一样想法的人。

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-5 08:51 , Processed in 0.425477 second(s), 26 queries .

快速回复 返回顶部 返回列表 Return to list