计算化学公社

标题: 【BAGEL使用心得】Gaussian、BAGEL联用做CASPT2几何优化的脚本 [打印本页]
作者
Author: Freeman    时间: 2021-12-23 15:25
标题: 【BAGEL使用心得】Gaussian、BAGEL联用做CASPT2几何优化的脚本
本帖最后由 Freeman 于 2022-7-31 17:00 编辑

【2022.07.31更新】


原来的脚本功能单一，代码不规范，应对不同的场景还要用户自己改代码，十分不便。现在gau_bag.sh升级重写了。细节如下。 (, 下载次数 Times of downloads: 50)

改进
1、原来的脚本只能应用于一阶导数的场景，现在可以应用于二阶导数的场景了。用GV观看CASPT2水平的振动模式，还能读各种热力学量，非常美滋滋。
2、原来BAGEL的输出，后一帧会把前一帧覆盖掉，没法复盘整个优化过程，现在每一帧的输出文件都会保存。
3、原来的脚本每次用都要改动以适应不同的核数、基组、活性空间、电子态的需求。现在这些信息都转到从gjf文件里的external的引号里读取了，用户就可以把脚本放到一个固定的地方而不用老是改了。
4、原来的Gaussian的输入文件里除了能量、梯度什么也没有，要看多态的能量和组态还要诉诸BAGEL的输出文件，很麻烦。现在BAGEL输出文件的重要信息也会打印到Gaussian输出文件里，用户就不用辗转于两个输出文件之间了。
5、原来的代码乱七八糟的，不规范，现在规范了很多，方便用户阅读修改。

案例：基态环丁二烯互变异构过渡态
这里用环丁二烯的例子讲述新的gau_bag.sh的用法。单重态环丁二烯的极小点是长方形，其长短边互变异构的过渡态是正方形。（图一）该过渡态是多参考体系，正好用CASPT2优化。
(, 下载次数 Times of downloads: 18)
首先，创建并进入一个新文件夹c4h4s，本节以下所有活动都在这个文件夹里进行。建模画出正方形的初猜结构，用BAGEL算个RHF单点，以提供后续CASSCF的活性空间。（文件c4h4s.bag.json）
  1. {"bagel":[
  2. {
  3.   "title":"molecule",
  4.   "basis":"cc-pvdz",
  5.   "df_basis":"cc-pvdz-jkfit",
  6.   "angstrom" : "true",
  7.   "geometry":[
  8.    { "atom" : "C", "xyz" : [        0.00000000,    1.01527700,    0.00000000]},
  9.    { "atom" : "C", "xyz" : [        1.01527700,    0.00000000,    0.00000000]},
  10.    { "atom" : "C", "xyz" : [       -1.01527700,    0.00000000,    0.00000000]},
  11.    { "atom" : "C", "xyz" : [        0.00000000,   -1.01527700,    0.00000000]},
  12.    { "atom" : "H", "xyz" : [        0.00000000,    2.09288200,    0.00000000]},
  13.    { "atom" : "H", "xyz" : [        2.09288200,    0.00000000,    0.00000000]},
  14.    { "atom" : "H", "xyz" : [        0.00000000,   -2.09288200,    0.00000000]},
  15.    { "atom" : "H", "xyz" : [       -2.09288200,    0.00000000,    0.00000000]}
  16.   ]
  17. },
  18. {
  19.   "title":"hf"
  20. },
  21. {
  22.   "title":"save_ref", // 保存波函数文件为c4h4s.bag.archive
  23.   "file":"c4h4s.bag"
  24. },
  25. {
  26.   "title":"print", // archive文件没法可视化，所以还要输出molden用来看轨道
  27.   "file":"c4h4s.bag.molden",
  28.   "orbitals":true
  29. }
  30. ]}
复制代码
用Multiwfn看molden文件的轨道。选取四个pi轨道+四个电子作为活性空间，用BAGEL算CASSCF。（文件c4h4s+.bag.json）如果活性轨道不正好依次排在HOMO、LUMO附近，为了调转轨道顺序（相当于Gaussian的guess(read,alter)关键词），这一步就必须做。
  1. {"bagel":[
  2. {
  3.   "title":"load_ref", // 从c4h4s.bag.archive读取前一步的分子坐标、基组和轨道
  4.   "file":"c4h4s.bag"
  5. },
  6. {
  7.   "title":"casscf",
  8.   "nact":4, // 四个活性轨道
  9.   "nclosed":12, // 满占据的非活性轨道数
  10.   "active":[13,14,15,20] // 活性轨道的编号，从Multiwfn看molden文件而得来的
  11. },
  12. {
  13.   "title":"print",
  14.   "file":"c4h4s+.bag.molden",
  15.   "orbitals":true
  16. },
  17. {
  18.   "title":"save_ref",
  19.   "file":"c4h4s+.bag"
  20. }
  21. ]}
复制代码
看看molden文件，发现轨道现在正好排在前线了。复制c4h4s+.bag.archive为reference.archive。优化用的gjf文件如下。（文件c4h4s.bag.gjf）
  1. %nprocshared=1 // 只用一个核就可以了，因为密集的计算都是BAGEL做的
  2. %chk=c4h4s.bag.chk // 保存最后的频率信息到chk文件里
  3. # opt=(ts,noeigen,calcfc,nomicro) freq nosym // nomicro和nosym关键词一定要加
  4. external='bash ~/gau_bag.sh 9 cc-pvdz cc-pvdz-jkfit 2 0 4 12' // 调用在家目录下的gau_bag.sh脚本，给BAGEL用9个核，基组用cc-pvdz，辅助基组用cc-pvdz-jkfit，算两个（单重）态（其实算一个态就可以，这里是为了举例），几何优化针对第1个态（0是第一个态，1是第二个态，以此类推），有4个活性轨道、12个满占据的非活性轨道

  6. Title Card Required

  8. 0 1 // 以前的脚本是不从gjf读电荷数和电子多重度的，而是要直接改写脚本内容；现在的脚本是从这里读的，所以不要乱设了
  9. C                  0.00000000    1.01527700    0.00000000
  10. C                  1.01527700    0.00000000    0.00000000
  11. C                 -1.01527700    0.00000000    0.00000000
  12. C                  0.00000000   -1.01527700    0.00000000
  13. H                  0.00000000    2.09288200    0.00000000
  14. H                  2.09288200    0.00000000    0.00000000
  15. H                  0.00000000   -2.09288200    0.00000000
  16. H                 -2.09288200    0.00000000    0.00000000

复制代码
运行Gaussian，届时文件夹里会从step_0000开始出现每一帧对应的包含BAGEL输入输出文件的文件夹，上一步的step_xxxx.bag.archive会复制成reference.archive作为下一步的初猜。结果如图二，有优化路线，有频率，有热力学量。
(, 下载次数 Times of downloads: 18)




******************************************************************************


BAGEL自己的优化算法不太行，经常是好不容易优化到平衡位置旁边了，突然一抖，能量上升0.0几个Hartree，就前功尽弃了，对低血压有很强的治疗效果。一个办法是降低置信半径至0.1左右并固定，这样优化时就不会振荡，但是会增加无意义的耗时。Gaussian的优化算法倒是很好，因此笔者写了个将Gaussian和BAGEL联用的脚本，将BAGEL的CASPT2解析梯度和Gaussian的优化算法结合起来。

Gaussian external接口脚本的基本思路
1、Gaussian把这一步的分子坐标、电荷、自旋多重度和要算的导数阶数（优化是一阶导数，频率是二阶导数）保存在一个文件里，文件名是Gau-[PID编号].EIn。
2、*读取EIn里的信息，生成其他程序的输入文件Gau-[PID编号].inp；
3、*设置环境变量，运行其他程序，获取其输出文件Gau-[PID编号].out；
4、*把out里的能量和梯度（或频率）信息复制到Gau-[PID编号].EOu里去；（BAGEL的能量和梯度信息可以另外分别保存在ENERGY.out和FORCE_0.out文件里，更方便读取）
5、Gaussian读取EOu里的信息，判断下一步的分子坐标该变成什么样；
6、如果是优化或IRC任务，重复1-5，直到收敛；如果是频率任务，就不用重复了。
以上打*的部分是接口脚本里要写的内容。每一步，Gaussian都会调用一次这个脚本，调用方式类似于“bash script.sh Gau-[PID编号].xxx（oniom里面用的，与本文无关） Gau-[PID编号].EIn Gau-[PID编号].EOu”，也就是说$2默认就是Gau-[PID编号].EIn，而$3默认就是Gau-[PID编号].EOu。

脚本使用流程
1、用现有的分子坐标编写一个用BAGEL的CASSCF算单点能的输入文件.json，如有必要就调整轨道顺序，注意要让它输出molden文件和archive文件，后者相当于Gaussian的chk文件；
2、用Multiwfn查看molden的轨道合不合心意，没问题的话就保存成.gjf格式；
3、改写.gjf文件，设置
  1. %nproc=1
  2. # opt(nomicro) external='bash gau_bagel.sh' nosym
复制代码
4、把gau_bagel.sh放到当前目录里，写好基组信息和CASSCF活性空间信息（脚本的65、19行；前一步已经调整好轨道顺序了，这里不用再调整了），添加BAGEL用到的环境变量（86-92行）；
5、运行Gaussian。
随后，脚本会创建一个名为“Gau-[PID编号]”的文件夹，然后把第一步的reference.archive文件复制进去，命名为Gau-[PID编号].archive。之后每走一步，都会读取上一步的archive作为轨道初猜，走完一步就更新archive里的信息，并且输出ENERGY.out和FORCE_0.out文件。最后在文件夹外面出现了Gaussian的输出文件。

效果图
(, 下载次数 Times of downloads: 27)

先运行[attach]43121[/attach]，获得初始molden和archive，然后运行带有这个脚本[attach]43120[/attach]的gjf。

P.S. BAGEL的输入文件格式非常变态，大家写输入文件的时候要注意一下方括号、大括号、逗号这些乱七八糟的东西有没有写全。

下期预告
测评一下包括BAGEL在内的常用来做多参考计算的程序在不同进程数下的速度。



作者
Author: Accelerator    时间: 2022-2-15 16:46
感谢楼主的工作。

楼主是否有测试过用该脚本优化激发态？我在尝试使用这段脚本进行构型优化时，优化基态非常正常，但优化激发态则会在算完初始构型后退出，而生成的各输出文件内容都是正常的。以下为.sh文件中的相应内容：
  1. cat>>$base/${base}.json<<EOF
  2. },
  3. {
  4. "title":"force",
  5. "target":1,
  6. "export":"true",
  7. "method":[{
  8. "title":"casscf",
  9. "nact":6,
  10. "nclosed":27,
  11. "nstate":3,
  12. "maxiter":100
  13. }]
  14. }
复制代码
作者
Author: Freeman    时间: 2022-2-15 17:14
Accelerator 发表于 2022-2-15 16:46
感谢楼主的工作。

楼主是否有测试过用该脚本优化激发态？我在尝试使用这段脚本进行构型优化时，优化基态 ...

95行改为awk 'NR==2{printf "%20.12e",$1}' $base/ENERGY.out > $outputfile
95行以下的所有FORCE_0改成FORCE_1
就好了
作者
Author: beefly    时间: 2022-5-5 21:31
好象一些简单的优化算法遇到碳环就容易出问题。比如cfour，经常会在快要收敛的时候，突然把碳环上的一个配体掰到相反的方向
作者
Author: lijingbai2009    时间: 2022-5-11 01:52
写BAGEL输入文件有个技巧，如果熟悉Python和json的话，就当Python脚本来写，用IDE也可以，这样括号，逗号不容易出错。
作者
Author: Freeman    时间: 2022-7-31 17:09
重写更新了一下。冒个泡。
作者
Author: ABQTrap    时间: 2023-3-8 14:44
老师您好！您提到的“一个办法是降低置信半径至0.1左右并固定，这样优化时就不会振荡，但是会增加无意义的耗时。” 解决优化时能量震荡问题，请问具体用哪个关键词设置？谢谢！
作者
Author: Freeman    时间: 2023-3-8 17:00
ABQTrap 发表于 2023-3-8 14:44
老师您好！您提到的“一个办法是降低置信半径至0.1左右并固定，这样优化时就不会振荡，但是会增加无意义的 ...

https://nubakery.org/opt/optimize.html#keywords
检索“maxstep”
作者
Author: ABQTrap    时间: 2023-3-8 17:15
Freeman 发表于 2023-3-8 17:00
https://nubakery.org/opt/optimize.html#keywords
检索“maxstep”

谢谢老师！还有一个问题：几何优化时，仅一轮优化，程序就自己停止了，都还没满足收敛条件。用自带的test文件也是这样，是编译的有问题吗？郁闷了
作者
Author: Freeman    时间: 2023-3-8 17:18
ABQTrap 发表于 2023-3-8 17:15
谢谢老师！还有一个问题：几何优化时，仅一轮优化，程序就自己停止了，都还没满足收敛条件。用自带的test ...

我也遇到过。所以我在安装bagel的贴子说了，不要下载1.2.2的release，而要下载github上最新的版本。后者没有这个问题。
作者
Author: ABQTrap    时间: 2023-3-9 16:49
Freeman 发表于 2023-3-8 17:18
我也遇到过。所以我在安装bagel的贴子说了，不要下载1.2.2的release，而要下载github上最新的版本。后者 ...

老师，试了github上的版本也是这样，几何优化一轮结束，要放弃了
作者
Author: ABQTrap    时间: 2023-3-23 00:05
老师好，感谢分享gau_bag.sh的脚本。我在使用时遇到了问题，麻烦您看一下。reference.archive已经预先提供了，运行step 0时正常，step 1时报错了。错误信息和输入文件见附件。好像是step 1时找不到所需文件。另外一个问题是，环境变量里面需要设置BAGEL_NUM_THREADS和MKL_NUM_THREADS吗？比如我是一个cpu共8核，关闭了超线程，应该如何设置呢？万分感谢！

作者
Author: Freeman    时间: 2023-3-23 01:45
step1的bagel任务正常结束了吗？有没有step0001.bag.archive？
两个都设为1。
作者
Author: ABQTrap    时间: 2023-3-23 08:25
本帖最后由 ABQTrap 于 2023-3-23 09:15 编辑
Freeman 发表于 2023-3-23 01:45
step1的bagel任务正常结束了吗？有没有step0001.bag.archive？
两个都设为1。

老师好，附件是step 1的输出文件，信息是ERROR: EXCEPTION RAISED:  changing geometry and basis set at the same time is not allowed. 但是我没有改变基组设定，都是cc-pvdz。step 1没有看到step0001.bag.archive。我的初始reference.archive是这么来的：保持坐标相同的情况下用ORCA在cc-pvdz水平得到gbw文件，转换为molden作为初始波函数，用于BAGEL计算。因为BAGEL总是找不到正确的波函数（C2H4的过渡态例子来自论坛http://bbs.keinsci.com/forum.php ... 9165&highlight=mrcc）。是不是这方面有影响？ORCA与BAGEL对cc-pvdz基组不一样？谢谢
作者
Author: ABQTrap    时间: 2023-3-23 10:55
本帖最后由 ABQTrap 于 2023-3-23 10:56 编辑
Freeman 发表于 2023-3-23 01:45
step1的bagel任务正常结束了吗？有没有step0001.bag.archive？
两个都设为1。

老师，目前已经确定是使用了ORCA的molden文件用于生成archive文件导致的错误。把生成的archive文件再用BAGEL算一遍CASSCF能量输出archive就好了。运行顺利了。想问下，对于自旋极化单重态的问题，有什么好的方法帮助收敛吗？目前采用了每一步都计算hessian（calcall），还是很难收敛，不像broken-symmetry DFT那么容易收敛。而且BAGEL总是找不到正确波函数，算出来的能量比ORCA、Gaussian高很多，后两者能量结果基本接近。谢谢！
作者
Author: Freeman    时间: 2023-3-23 17:58
本帖最后由 Freeman 于 2023-3-23 18:00 编辑
ABQTrap 发表于 2023-3-23 10:55
老师，目前已经确定是使用了ORCA的molden文件用于生成archive文件导致的错误。把生成的archive文件再用BA ...

我记得bagel不能读取molden波函数作为初猜。
既然用到bagel了，就不要考虑自旋极化单重态波函数了。你可以直接用bagel给出的rhf轨道当初猜算casscf，毕竟这个体系casscf波函数本来就是两个闭壳层多项式的线性叠加。
优化不收敛有没有可能是因为活性空间没选好呢？

作者
Author: ABQTrap    时间: 2023-3-23 21:18
Freeman 发表于 2023-3-23 17:58
我记得bagel不能读取molden波函数作为初猜。
既然用到bagel了，就不要考虑自旋极化单重态波函数了。你可 ...

老师好，我用的乙烯过渡态例子来自于论坛http://bbs.keinsci.com/forum.php?mod=viewthread&tid=29165&highlight=fic+mrcc，并且用ORCA和高斯重现过结果，应该没有问题，但是bagel几何优化总是产生奇怪的结果。正常过渡态虚频应该在1700-1800cm-1左右，bagel的CASSCF算出来的总是不对，所有很纳闷。可能输入文件有地方不对。
作者
Author: Freeman    时间: 2023-3-23 21:52
ABQTrap 发表于 2023-3-23 21:18
老师好，我用的乙烯过渡态例子来自于论坛http://bbs.keinsci.com/forum.php?mod=viewthread&tid=29165&hi ...

你检查过活性空间没有？
作者
Author: ABQTrap    时间: 2023-3-27 21:19
Freeman 发表于 2023-3-23 21:52
你检查过活性空间没有？

老师好，您提醒我了。因为之前用ORCA和高斯的时候用的都是CASSCF（2electrons，2orbitals），所以在BAGEL里面计算时想当然照用了。今天尝试了BAGEL的CASSCF（4electrons，4orbitals），BAGEL结果都对了，与ORCA和高斯很接近。至于为什么BAGEL需要更大的活性空间，我还不了解，请老师解答。感谢老师这几天的详细回复！
作者
Author: Freeman    时间: 2023-3-28 03:05
ABQTrap 发表于 2023-3-27 21:19
老师好，您提醒我了。因为之前用ORCA和高斯的时候用的都是CASSCF（2electrons，2orbitals），所以在BAGEL ...

感谢啥啊。我明明没帮上什么忙。我看到你到bagel的github主页上也问问题了，给出的输入文件，一看就知道完全没理睬我第16、18楼的建议。如果听不懂，哪怕问一问呢？谁曾经还不是个菜鸟呢？但你偏不，就自己乱搞，结果白花了这5天时间，至今不知道问题出在哪里。我想帮忙也帮不上。现在歪打正着，别感谢我，感谢运气吧。阿弥陀佛，善哉善哉。
作者
Author: ABQTrap    时间: 2023-3-28 08:11
Freeman 发表于 2023-3-28 03:05
感谢啥啊。我明明没帮上什么忙。我看到你到bagel的github主页上也问问题了，给出的输入文件，一看就知道 ...

还是要感谢老师提醒检查活性空间
作为菜鸟觉得CASSCF太难了
目前只是测试和熟悉软件，所以用了简单的例子，没想到这么曲折，后面算自己的东西估计更头疼，到时再请教老师。



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3