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

Freeman

【2022.07.31更新】


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

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

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