计算化学公社

 找回密码 Forget password
 注册 Register
Views: 52508|回复 Reply: 48

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

  [复制链接 Copy URL]

877

帖子

36

威望

4805

eV
积分
6402

Level 6 (一方通行)

发表于 Post on 2019-7-30 11:02:16 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 ggdh 于 2021-8-1 17:00 编辑

和本脚本功能类似,但用于orca结果监控的程序见ORCA几何优化的监控/诊断小脚本-orcamonitor

6.0版新特性(先立flag)
目前弄大文件的时候有点慢,6.0版应该会重构部分代码,加速(然后别人问我重构什么代码,我说重构shell代码。。。。
其他新特性暂时想不到,等待沙雕网友们补充

5.0版新特性
1. 测量每一步的距离,角度,二面角
用法如下
  1. tsmonitor5 -m "d1-2 a1-2-3 t1-2-3-4" tsopt.log
复制代码
其中d是distance,a是angle,t是torsion,可以一次测多个指标,中间用空格隔开,用引号引起来即可。对opt,柔性扫描,irc任务都适用
可能的使用场景:
a)找过渡态时,通过观察成键/断键的原子的距离,可以更加确定过渡态的健康程度
b)柔性扫描时,可以观察其他的一些结构参数伴随扫描坐标的改变,也可以直接测量扫描的坐标,然后和能量/玻尔兹曼分布那些一起拷贝出来方便作图

2. 从中途提取一个或者多个结构并保存为gjf文件
可能的使用场景:
a)找过渡态小连招
在进行经典的柔性扫描--找过渡态小连招过程中,我们需要取柔性扫描最高点附近的点作为过渡态的初始猜,这里假设我们先通过tsmonitor发现第9步是最高点,然后再运行:
  1. tsmonitor5 -e "s8-10"  -K  "opt(ts,noeigentest,calcall) m062x def2svp em(gd3)" -N 20 -M 20GB  relaxscan.log  
复制代码
这里-e "s8-10"是提取柔性扫描中8,9,10这三步的坐标,并在当前目录下产生优化过渡态的输入文件。
输出文件名是relaxscan_s8.gjf relaxscan_s9.gjf relaxscan_s10.gjf (因为输入文件名是relaxscan.log)
-K -N 和 -M分别设置guassian输入文件的关键词,核数和内存。
如果你计算资源充足,可能同时对这3个任务进行计算,这样成功找到过渡态的概率更大。
b)中途介入,挽救过渡态
有时候过渡态前4步还正常,第5步突然跟丢了,那这时候我们可以把第四步捞出来,减小其步长,增加calcall关键词,看能不能挽救:
  1. tsmonitor5 -e "o4"  -K  "opt(ts,noeigentest,calcall,maxstep=2,notrust) m062x def2svp em(gd3)"  tsopt.log  
复制代码
注意这里o4是指优化任务(opt)的第四步。
优化基态或者激发态也可以类似的操作,有时候优化任务时候会出现震荡,那我可以在最后的振荡结构中提出能量最低的那个结构,然后减少步长继续算,或者直接进行主观收敛。
c)结合Multiwfn制作炫酷动画
有时候我们需要把IRC或者是扫描的中间结构拿出来,算单点获得波函数,进而制作炫酷动画,参考:
通过键级曲线和ELF/LOL/RDG等值面动画研究化学反应过程
产生Gaussian的IRC和SCAN任务每个点的波函数文件的工具
用我这个工具也可以完成和sob那个工具类似的功能,操作如下:
  1. tsmonitor5 -e "r1-10 f1-10"  -K  "m062x def2svp"  irc.log  
复制代码
这里r1-10和f1-10对应irc的正反向各10步(注意这里面没用包含过渡态本身的结构),另外柔性扫描也能搞动画,比如:
  1. tsmonitor5 -e "s1-16"  -K  "td cam-b3lyp def2svp density scrf(pcm)"  tictscan.log  
复制代码
做一个TICT激发态势能面扫描,然后做激发态电子密度/空穴电子分布随角度的变化动画。
d)普通续算
有时候我们简单的需要最后一个结构进行后续计算,比如优化完了算个性质什么的:
  1. tsmonitor5 -e e opt.log  
复制代码
-e e 就是用最后一个结构取产生gjf文件。

3. 修复bug若干
4. 输出排版进一步优化

4.0版新特性(脚本介绍在后面,第一次接触的同学建议先看后面的脚本介绍):
根据各位热心网友(尤其是冰释之川)的建议增加了下面一些特性:
1.IRC的输出信息更丰富了,增加scf cycle,相对上步能差(dE(kCal)),MAX force, RMS force,,
2.普通优化的任务,增加相对上一步能量差,方便快速判断能量收敛情况
3.柔性扫描任务,增加每个点的boltzmann分布(用-t选项设温度),方便看不同构象之间的比例
4.输出宽度超出窗口宽度时,不会自动换行,影响观看。

3.0版新特性(脚本介绍在后面,第一次接触的同学建议先看后面的脚本介绍):
根据各位热心网友的建议增加了下面一些特性:
1.对于有多个Link的Gaussian任务,会分别分析每个link
2.对于激发态优化,显示跃迁波长
3.对于IRC任务,优化了输出结果如下,其中F代表Forward,R代表Rerverse:
irc.png
4.显示了每一步的SCF圈数
5.对于普通的优化任务,也显示了本征值(如果都是0,那说明这步没有本征值)
6.改进排版


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, 下载次数 Times of downloads: 170

tsmonitor2

5.93 KB, 下载次数 Times of downloads: 140

tsmonitor3

7.33 KB, 下载次数 Times of downloads: 102

tsmonitor4

9.99 KB, 下载次数 Times of downloads: 126

tsmonitor5

21.37 KB, 下载次数 Times of downloads: 147

评分 Rate

参与人数
Participants 44
威望 +1 eV +191 收起 理由
Reason
lnf + 5 好物!
1983984613 + 5 谢谢分享
lurensan + 3 赞!
mfdsrax2 + 5 GJ!
北大-陶豫 + 5 好物!
丁越 + 5 钟叔大法好
hong1kun + 4 好物!
jqfeng + 5 精品内容
ene + 5
funok + 5 谢谢分享
panernie + 4
AceVolca + 5 太好用了
Satoru + 5 хорошо!
djjj148 + 5 好物!
朙天儿 + 5 赞!
RAL + 3 赞!
壹零壹室掃地僧 + 3 赞!
philartist + 5 赞!
typ_gauss + 4 好物!
cyx98 + 4 好物!

查看全部评分 View all ratings

22

帖子

0

威望

1887

eV
积分
1909

Level 5 (御坂)

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

1020

帖子

0

威望

3510

eV
积分
4530

Level 6 (一方通行)

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

评分 Rate

参与人数
Participants 1
eV +2 收起 理由
Reason
灰飞的旋律 + 2 我很赞同

查看全部评分 View all ratings

2

帖子

0

威望

17

eV
积分
19

Level 1 能力者

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

4万

帖子

99

威望

4万

eV
积分
89975

管理员

公社社长+计算化学玩家

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

评分 Rate

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

查看全部评分 View all ratings

北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

211

帖子

0

威望

2304

eV
积分
2515

Level 5 (御坂)

埋進雨裡聽海的呼吸

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

111

帖子

0

威望

1774

eV
积分
1885

Level 5 (御坂)

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

3

帖子

0

威望

218

eV
积分
221

Level 3 能力者

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

877

帖子

36

威望

4805

eV
积分
6402

Level 6 (一方通行)

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

3

帖子

0

威望

218

eV
积分
221

Level 3 能力者

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

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

88

帖子

0

威望

1941

eV
积分
2029

Level 5 (御坂)

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

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

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

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

评分 Rate

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

查看全部评分 View all ratings

307

帖子

1

威望

5739

eV
积分
6066

Level 6 (一方通行)

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

88

帖子

0

威望

1941

eV
积分
2029

Level 5 (御坂)

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

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

评分 Rate

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

查看全部评分 View all ratings

877

帖子

36

威望

4805

eV
积分
6402

Level 6 (一方通行)

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

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

877

帖子

36

威望

4805

eV
积分
6402

Level 6 (一方通行)

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

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

本版积分规则 Credits rule

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

GMT+8, 2023-2-7 03:49 , Processed in 0.251199 second(s), 31 queries .

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