【BAGEL使用心得】Gaussian、BAGEL联用做CASPT2几何优化的脚本

Freeman

【2022.07.31更新】


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

2022-7-31 16:56 上传 Uploaded

点击下载
Click to download
新脚本

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

案例：基态环丁二烯互变异构过渡态
这里用环丁二烯的例子讲述新的gau_bag.sh的用法。单重态环丁二烯的极小点是长方形，其长短边互变异构的过渡态是正方形。（图一）该过渡态是多参考体系，正好用CASPT2优化。

首先，创建并进入一个新文件夹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个满占据的非活性轨道
   
5. 
   
6. Title Card Required
   
7. 
   
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
    
17. 
    

复制代码

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





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


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的输出文件。

效果图


先运行，获得初始molden和archive，然后运行带有这个脚本的gjf。

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

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

Accelerator

感谢楼主的工作。

楼主是否有测试过用该脚本优化激发态？我在尝试使用这段脚本进行构型优化时，优化基态非常正常，但优化激发态则会在算完初始构型后退出，而生成的各输出文件内容都是正常的。以下为.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. }
    

复制代码

Freeman

Accelerator 发表于 2022-2-15 16:46
感谢楼主的工作。

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

95行改为awk 'NR==2{printf "%20.12e",$1}' $base/ENERGY.out > $outputfile
95行以下的所有FORCE_0改成FORCE_1
就好了

beefly

好象一些简单的优化算法遇到碳环就容易出问题。比如cfour，经常会在快要收敛的时候，突然把碳环上的一个配体掰到相反的方向

lijingbai2009

写BAGEL输入文件有个技巧，如果熟悉Python和json的话，就当Python脚本来写，用IDE也可以，这样括号，逗号不容易出错。

Freeman

重写更新了一下。冒个泡。

ABQTrap

老师您好！您提到的“一个办法是降低置信半径至0.1左右并固定，这样优化时就不会振荡，但是会增加无意义的耗时。” 解决优化时能量震荡问题，请问具体用哪个关键词设置？谢谢！

Freeman

ABQTrap 发表于 2023-3-8 14:44
老师您好！您提到的“一个办法是降低置信半径至0.1左右并固定，这样优化时就不会振荡，但是会增加无意义的 ...

https://nubakery.org/opt/optimize.html#keywords
检索“maxstep”

ABQTrap

Freeman 发表于 2023-3-8 17:00
https://nubakery.org/opt/optimize.html#keywords
检索“maxstep”

谢谢老师！还有一个问题：几何优化时，仅一轮优化，程序就自己停止了，都还没满足收敛条件。用自带的test文件也是这样，是编译的有问题吗？郁闷了

Freeman

ABQTrap 发表于 2023-3-8 17:15
谢谢老师！还有一个问题：几何优化时，仅一轮优化，程序就自己停止了，都还没满足收敛条件。用自带的test ...

我也遇到过。所以我在安装bagel的贴子说了，不要下载1.2.2的release，而要下载github上最新的版本。后者没有这个问题。

ABQTrap

Freeman 发表于 2023-3-8 17:18
我也遇到过。所以我在安装bagel的贴子说了，不要下载1.2.2的release，而要下载github上最新的版本。后者 ...

老师，试了github上的版本也是这样，几何优化一轮结束，要放弃了

ABQTrap

老师好，感谢分享gau_bag.sh的脚本。我在使用时遇到了问题，麻烦您看一下。reference.archive已经预先提供了，运行step 0时正常，step 1时报错了。错误信息和输入文件见附件。好像是step 1时找不到所需文件。另外一个问题是，环境变量里面需要设置BAGEL_NUM_THREADS和MKL_NUM_THREADS吗？比如我是一个cpu共8核，关闭了超线程，应该如何设置呢？万分感谢！

Freeman

step1的bagel任务正常结束了吗？有没有step0001.bag.archive？
两个都设为1。

ABQTrap

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基组不一样？谢谢

ABQTrap

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高很多，后两者能量结果基本接近。谢谢！