（更新5.0版）Gaussian优化过渡态的监控/诊断小脚本-tsmonitor

ggdh

和本脚本功能类似，但用于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：

4.显示了每一步的SCF圈数
5.对于普通的优化任务，也显示了本征值(如果都是0，那说明这步没有本征值）
6.改进排版


2.0版新特性（脚本介绍在后面，第一次接触的同学建议先看后面的脚本介绍）：
根据各位热心网友的建议增加了下面一些特性：
1.显示收敛情况，也就是Gaussian的4个收敛指标，如下图所示

这4列分别式Max Force，RMS Force，Max Displacement， RMS Displacement中的当前指标和收敛限的比值。
如果小于1，说明对应项达到了收敛指标，越大就离收敛指标越远。
用途：通过观察这一项，可以判断某些震荡体系离收敛还有多远的距离。

2. 支持读取柔性扫描的log文件，如下图所示

用途：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）的编号

这里优化的是一个还原消除反应的过渡态，这个反应中11-45和12-45键断裂，而11-12键生成，因此这步反应中的关键原子（AOI）被定为11 12 45

2. 运行脚本进行诊断
根据上面AOI的编号，输入命令，这里TS45a1.log是本例反应的log文件。

1. tsmonitor -a "11 12 45"  TS45a1.log

复制代码

可以得到如下结果：

下面就这个结果中每一部分对应的意思做下简介（能够直接看懂结果的可以略过此部分）：
第一行：这个脚本首先会检测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部分观察。

优化失败的例子：
下面再举一个优化失败的例子，和上面成功的例子是同一个过渡态，但是初猜不同，我们使用相同的命令进行分析

可以看到在三步的时候，虚频突然变得很小，而且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 改可执行权限（复制了上面代码的同学略过此步）


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

chemzhh

谢谢分享，正在找过渡态中，使用后再来汇报

granvia

支持！ 这个程序固然是好，但可能也失去了“意外惊喜发现”的机会 :D

杨开心

谢谢分享

sobereva

试了下，很好用。已经把这个程序结合我的实际例子加入到科音基础量化班ppt里进行推广

Shine剪水

谢谢钟老师！不仅是方便了监控，还进一步加深了对优化过程细节的认识！

409407227

谢谢分享

浪里白条

有个小问题想请教一下钟老师，程序里面eigenV=$(grep -A 1 ITU $1 | grep Eigenvalues | awk '{print $3" "$4" "$5" "$6"\n"}')，可我的log文件中没有“ ITU”关键词，因此运行脚本 没有得到Eigenvalues ， 请问只有我一个人是这样的吗？

ggdh

请问你计算用的关键词是什么啊?

浪里白条

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

mutron

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

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

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

zsu007

谢谢分享！

mutron

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

ggdh

浪里白条 发表于 2019-8-2 14:36
我的关键词是：#p b3lyp/genecp  fopt=(ts,noeigen,calcfc) optcyc=300 freq

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

ggdh

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

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