请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 909|回复: 16

[Gaussian/gview] (更新2.0版)Gaussian优化过渡态的监控/诊断小脚本-tsmonitor

[复制链接]

590

帖子

14

威望

2074

eV
积分
2944

Level 5 (御坂)

发表于 2019-7-30 11:02:16 | 显示全部楼层 |阅读模式
本帖最后由 ggdh 于 2019-8-15 16:25 编辑

2.0版新特性(脚本介绍在后面,第一次接触的同学建议先看后面的脚本介绍):
根据各位热心网友的建议增加了下面一些特性:
1.显示收敛情况,也就是Gaussian的4个收敛指标,如下图所示
TSMnew.png
这4列分别式Max Force,RMS Force,Max Displacement, RMS Displacement中的当前指标和收敛限的比值。
如果小于1,说明对应项达到了收敛指标,越大就离收敛指标越远。
用途:通过观察这一项,可以判断某些震荡体系离收敛还有多远的距离。

2. 支持读取柔性扫描的log文件,如下图所示
TMSnew2.png
用途:1,找过渡态之前往往会用柔性扫描找初始猜。心急的同学往往等不到扫描结束,只要出现最高点就可以中止任务,然后拿最高点,或者最高点的前一个点来做初猜。
2,通过最高点的能量,可以估计一下这个过渡态是否合理。


3.增加-c选项,比如
  1. tsmonitor -a "11 12 45"  -c a TS45a1.log
复制代码
可以输出所有的出现过的内坐标对应的原子序号。
用途:方便大家查看优化中间出现过哪些虚频

4.增加step列,方便查看当前优化了几步,或者是从中间取结构做续算
用途:1,某些震荡体系中间会出现能量最小值,取这个结构往往比取最后一帧有利。
2,TS优化过程中某一步很健康,其后健康状况恶化,可以取这一步加calcfc续算。

5.支持读取关键词以#p开头的激发态优化中基态和激发态的能量

6.增强排版,避免Eigenvector不对齐的情况出现,Eigenvalue 现在显示的值乘以了100,避免一行太长以及方便肉眼阅读。

欢迎大家找bug,提出建议!



下面开始脚本介绍:

前言:
在使用Gaussian优化激发态时,我们通常使用opt(ts,calcfc,noeigentest)关键词来寻找过渡态,参考:
简谈Gaussian里找过渡态的关键词opt=TS和QST2、3
如何判断是否得到了想要的过渡态,这里引用sob上面这个帖子里面的原话:
“怎么验证过渡态找没找对?如果优化出来的过渡态结构在期望的反应物和产物之间,同时有且只有一个虚频,一般能有6成把握认为过渡态找对了。如果通过gview查看,发现仅有的那个虚频的振动方向对应反应坐标方向,那8成可以认为过渡态找对了;如果走IRC还能连通自己期望的反应物和产物结构,那么说明100%找对了。”
在一些情况下,由于初猜不够合理,导致开始优化就没有走到正确的轨迹上,或者是优化到一半跑偏了(优化过渡态就像是玩某个即时跟踪类游戏,一旦跟丢目标,往往很难找回来),或者出现震荡,如果我们不闻不问,只等最后的结果,无疑会造成cpu机时和国家电力资源的巨大浪费。
其实这时候我们可以通过观察log文件中振动对应的Eigenvalues 和 Eigenvectors来进行早期诊断,如果发现情况不妙,可以及时修改结构初猜,中止没有前途的计算,避免浪费。
基于这个原理,这里编了一个bash脚本,用于诊断和监控Gaussian过渡态优化的情况。下面举例介绍一下用法。

脚本用法:

1. 确定关键原子AOI(Atom of Interest)的编号
tsm1.png
这里优化的是一个还原消除反应的过渡态,这个反应中11-45和12-45键断裂,而11-12键生成,因此这步反应中的关键原子(AOI)被定为11 12 45

2. 运行脚本进行诊断
根据上面AOI的编号,输入命令,这里TS45a1.log是本例反应的log文件。
  1. tsmonitor -a "11 12 45"  TS45a1.log
复制代码
可以得到如下结果:
tsm2.png
下面就这个结果中每一部分对应的意思做下简介(能够直接看懂结果的可以略过此部分):
第一行:这个脚本首先会检测Gaussian运算是否正常结束,第一行如果是Normal termination...就是正常结束,如果是Error termination...就是错误结束,如果没有第一行,那么可能是异常中止或者正在运行。
A部分:A部分列出了相对第一步的能量(单位kCal/mol),主要用于判断收敛情况以及是否有震荡。
B部分:B部分列出了前4个最小的振动本征值,一个健康的过渡态优化过程中应该始终保持一个负的绝对值较大(小于-0.01)的本征值,最后的频率大概是这个值乘以10000,最后收敛时应只剩一个负值。
C部分:C部分列出了过渡态虚频对应的内坐标,以及这些坐标中是否含有AOI,比如D22-12,45的意思是虚频对应的第22号二面角中包含了AOI中的12 45号原子。一个健康的的过渡态优化过程中,AOI应该尽可能多的参与到虚频对应的内坐标中,在本例中,我们定义的11 12 45号原子,始终参与到了虚频对应的内坐标中。
D部分:列出最后一步虚频内坐标对应的原子编号。有时候AOI不出现在内坐标当中,我们又想知道这些内坐标到底对应哪些原子,可以分别复制这些原子编号到gview中的atom selection部分观察。

优化失败的例子:
下面再举一个优化失败的例子,和上面成功的例子是同一个过渡态,但是初猜不同,我们使用相同的命令进行分析
TSM3.png
可以看到在三步的时候,虚频突然变得很小,而且AOI从内坐标中消失。这说明这个过渡态优化已经“病入膏肓”,或者可以理解为gaussian“跟丢了”过渡态,可以“放弃治疗”了。果然,后续的优化都无法回到正轨上了。

后记:
  • 根据我非常有限的过渡态优化经验,一个健康的过渡态的Eigenvalue 应该小于-0.01(绝对值大于0.01),并且应该有大于1个的AOI参与到虚频对应的内坐标中。这里仅供参考。实际使用时,大家可以根据实际过程的不健康程度来决定是否放弃治疗。另外,值得一提的是,如果不使用calcall,那么Eigenvalue的值可能并不准。
  • 本脚本默认显示前5个虚频对应的内坐标,使用脚本时可以用-v选项来控制显示的内坐标的数量,范围是1到10个。
  • 运行脚本时可以不用-a选项指定AOI,这时候可以根据Eigenvalue大小来进行初略判断,同时,如果前5个内坐标全部是D打头的,那这个状态可能也不太健康。
  • 可以使用下面命令把屏幕上的输出同时保存到文件中
    1. tsmonitor -a “11 13  45”  TSopt.log | tee TSopt-mon.log
    复制代码

  • 对于普通优化的任务,也可也使用本脚本查看能量相对第一步的变化以及相对上一步的变化,单位是kCal/mol,比grep Done稍微友好一点。

脚本的安装方法:
1,不支持bash3 版本,使用前请运行bash --version 确定你的bash是4.x版本
2,拷贝到PATH中的某个路径下(不知道怎么操作的同学,首先确定你的当前目录下有tsmonitor这个文件,然后的复制下面代码到terminal即可)
  1. mkdir ~/scripts
  2. cp tsmonitor ~/scripts/tsmonitor
  3. chmod +x ~/scripts/tsmonitor
  4. echo 'export PATH=~/scripts:$PATH' >> ~/.bashrc
  5. source ~/.bashrc
复制代码

3,chmod +x tsmonitor 改可执行权限(复制了上面代码的同学略过此步)


最后所有苦哈哈找过渡态的同学们共勉,我们太难了。。。大家有啥意见和建议,欢迎留言


tsmonitor

3.35 KB, 下载次数: 56

tsmonitor2

5.93 KB, 下载次数: 7

评分

参与人数 14威望 +1 eV +63 收起 理由
lidocaine39 + 5 谢谢
wfmf1994 + 5 超有用,谢谢
shenchaoren + 4 好物!
浪里白条 + 5 牛!
王二葛 + 5 好物!
终_焉 + 5 赞!
复前行79 + 5 GJ!
Shine剪水 + 5 谢谢~
Novice + 5 好物!
zsu007 + 5 牛!
sobereva + 1
ABetaCarw + 5 好物!
req + 4 とてもいい!
小范范1989 + 5 好物!

查看全部评分

18

帖子

0

威望

775

eV
积分
793

Level 4 (黑子)

发表于 2019-7-30 11:51:24 | 显示全部楼层
谢谢分享,正在找过渡态中,使用后再来汇报

497

帖子

0

威望

911

eV
积分
1408

Level 4 (黑子)

发表于 2019-7-30 12:01:03 来自手机 | 显示全部楼层
支持! 这个程序固然是好,但可能也失去了“意外惊喜发现”的机会 :D

2

帖子

0

威望

17

eV
积分
19

Level 1 能力者

发表于 2019-7-30 15:19:15 | 显示全部楼层
谢谢分享

1万

帖子

25

威望

2万

eV
积分
44706

管理员

公社社长

发表于 2019-7-30 18:43:50 | 显示全部楼层
试了下,很好用。已经把这个程序结合我的实际例子加入到科音基础量化班ppt里进行推广

评分

参与人数 3eV +13 收起 理由
阿锋001 + 5 期待!
ggdh + 5 哈哈哈哈,开心啊!
Novice + 3 好物!

查看全部评分

北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

146

帖子

0

威望

1076

eV
积分
1222

Level 4 (黑子)

埋進雨裡聽海的呼吸

发表于 2019-7-30 20:28:51 | 显示全部楼层
谢谢钟老师!不仅是方便了监控,还进一步加深了对优化过程细节的认识!
奔跑吧 驕傲的少年

64

帖子

0

威望

451

eV
积分
515

Level 4 (黑子)

发表于 2019-7-31 07:47:59 | 显示全部楼层
谢谢分享

3

帖子

0

威望

168

eV
积分
171

Level 3 能力者

发表于 2019-8-1 11:07:31 | 显示全部楼层
有个小问题想请教一下钟老师,程序里面eigenV=$(grep -A 1 ITU $1 | grep Eigenvalues | awk '{print $3" "$4" "$5" "$6"\n"}'),可我的log文件中没有“ ITU”关键词,因此运行脚本 没有得到Eigenvalues , 请问只有我一个人是这样的吗?

590

帖子

14

威望

2074

eV
积分
2944

Level 5 (御坂)

 楼主| 发表于 2019-8-1 11:27:07 | 显示全部楼层
请问你计算用的关键词是什么啊?

3

帖子

0

威望

168

eV
积分
171

Level 3 能力者

发表于 2019-8-2 14:36:50 | 显示全部楼层
本帖最后由 浪里白条 于 2019-8-2 14:43 编辑

我的关键词是:#p b3lyp/genecp  fopt=(ts,noeigen,calcfc) optcyc=300 freq

58

帖子

0

威望

1618

eV
积分
1676

Level 5 (御坂)

发表于 2019-8-4 17:25:38 | 显示全部楼层
本帖最后由 mutron 于 2019-8-4 20:24 编辑

不知道可否把D部分修改为:所有值小于-0.01对应的内坐标原子编号,以及,如果一行里面只剩下一个负值,输出其内坐标对应原子编号。

因为原D部分输出除了虚频内坐标,其它好像没太多作用。如果碰到跟丢过渡态的情况,我们比较感兴趣的是到底跟错了哪个虚频,以及最初最大那几个负值Eigenvalue有没有包含感兴趣的方向

当然内坐标对应原子编号可以在log文件里面找,但是把感兴趣的虚频内坐标一起输出来看会方便很多,不用来回切换找

评分

参与人数 1eV +5 收起 理由
ggdh + 5 好建议,我更新了以后回复你

查看全部评分

159

帖子

1

威望

3199

eV
积分
3378

Level 5 (御坂)

发表于 2019-8-4 17:53:29 | 显示全部楼层
谢谢分享!

58

帖子

0

威望

1618

eV
积分
1676

Level 5 (御坂)

发表于 2019-8-4 19:59:31 | 显示全部楼层
本帖最后由 mutron 于 2019-8-4 20:14 编辑

写错编辑掉
可以考虑左边加上优化的步数序号,有时候看到哪一步开始走偏了定位会比较快,优化步数多了不用一行行数下去
0.png

评分

参与人数 1eV +5 收起 理由
ggdh + 5 该功能已添加

查看全部评分

590

帖子

14

威望

2074

eV
积分
2944

Level 5 (御坂)

 楼主| 发表于 2019-8-5 22:11:48 | 显示全部楼层
浪里白条 发表于 2019-8-2 14:36
我的关键词是:#p b3lyp/genecp  fopt=(ts,noeigen,calcfc) optcyc=300 freq

我试了你的关键词 其中optcyc=300不认,说我语法错误。是不是你的Gaussian版本太老?
我建议你使用g09的版本。
把这个关键词删掉以后可以正常运行脚本。

590

帖子

14

威望

2074

eV
积分
2944

Level 5 (御坂)

 楼主| 发表于 2019-8-5 22:17:18 | 显示全部楼层
mutron 发表于 2019-8-4 17:25
不知道可否把D部分修改为:所有值小于-0.01对应的内坐标原子编号,以及,如果一行里面只剩下一个负值,输出 ...

需要知道的是,所有列出的Eigenvector 都对应第一个,也是绝对值最大的那个虚频。
而且gaussian 的log里面也不会列出其他虚频对应的本征矢量。
在2.0版本中,加了-c 选项。 运行命令的时候输入 -c a 可以打印出所有出现过的内坐标及其对应的原子列表。
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949-1号 )

GMT+8, 2019-8-21 20:46 , Processed in 0.182893 second(s), 28 queries .

快速回复 返回顶部 返回列表