计算化学公社

 找回密码 Forget password
 注册 Register
Views: 7383|回复 Reply: 2

[ORCA] (新增一个ORCA保留scf每一步结果的独立脚本)给sob老师的gau_orca添加额外的功能

[复制链接 Copy URL]

112

帖子

1

威望

933

eV
积分
1065

Level 4 (黑子)

发表于 Post on 2020-8-1 13:57:54 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 ldatea 于 2020-8-14 13:00 编辑

文会根据需要,持续更注意:附件的程序使用的是opt=z-matrix,如果想要用其他关键词,直接在route section写上关键词是没用的,需要修改IOP产生新的计算路径我建议所有EF数值优化都使用opt=z-matrix,因其不仅可以使用纯粹z-matrix,还可以使用笛卡尔坐标和两者混合的坐标。
2020.8.4
更新ORCA的scf每一步的独立脚本gbwsave.sh,若收敛不正常,可以取出某一步的结果,以后将与其他的整合到一起。

今天早上传的版本有两个Bug,一个问题大点,一个问题小一点,目前暂未修复(这两个Bug解决需要大改代码,我现在的能力可能不够用,暂时没有更新计划,最后一次更新增加了输出部分,以便于观察运行情况及时发现bug)
波函数收敛后有一定概率触发,触发条件是:在orca_scf(或者orca_scf_mpi)还在运行的时候,如果.out文件写入比较慢,gbwsave.sh读取最后一行的命令读到了ORBITAL ENERGIES之后的输出,这个输出和scf部分的输出很像(前三个字符是空格+数字),
gbwsave.sh错误地将此部分的前三个字符当成scf收敛圈数,就错误地把某次的.gbw给覆盖了(或者造了一个本来不存在的)。在某个地方加一句sleep 或许可以减少发生概率。
我是初学者,尽可能详细的写了注释(怕我某天翻出来自己都看不懂),写的脚本可能出现很多复杂的效率低的部分或者bug,希望各位可以提出宝贵意见。


目录
1.简单的改进(直接贴代码)
2.用Gaussian的external功能调用ORCA做EF数值梯度几何优化。(修改IOP)
3.暂时未整合的独立脚本。
4.附件说明
5.写脚本时遇到的困难




1.简单的改进


(8.4更新)
一些很简单的功能,就直接贴代码了

如果需要获得每次调用ORCA得到的波函数(对IRC可能比较有用)
可以把sob老师的orca.sh最后三行前(清理目录)前添加
  1. gbwnumber=$(ls|grep mol.gbw |awk 'END{print NR}')                         #获得当前目录所有.gbw文件的数目
  2. #gbwnumber=$(echo 0000${gbwnumber})      #如果希望输出都按顺序排列,需要把行首注释去掉
  3. cp mol.gbw ${gbwname}mol.gbw
复制代码

这样每一次运行后的波函数都会被保留并且编号,但可能会产生一大堆.gbw,注意留足够的磁盘空间







2.用Gaussian的external功能调用ORCA做EF数值梯度几何优化。(修改IOP)

有关gau_orca的相关信息,见http://bbs.keinsci.com/thread-10141-1-1.html(《将Gaussian与ORCA联用搜索过渡态、产生IRC、做振动分析》)
ORCA的有不少功能是非常好用的,但是ORCA的几何优化不够理想,beefly老师专门吐槽过ORCA的数值梯度几何优化http://bbs.keinsci.com/thread-6860-1-1.html(不知道现在有没有改进)Gaussian的external功能挺好用的,就是遇到一些特殊情况非常恶心。(见下文)
我本来想用orca_的DLPNO-CCSD(T)给小分子几何优化,用gau_orca感觉挺不错。
ccsd(t)数值优化的计算路径
  1. 1/18=40,29=20000,38=1,172=1/1,14;
  2. 2/12=2,17=6,18=5,29=3,40=1/2;
  3. 3/11=9,25=1,30=1/1,2,3;
  4.                   4//1;
  5.                   5/5=2,38=5/2;
  6.                   8/6=4,9=120000,10=1/1,4;
  7.                   9/5=7,14=2/13;
  8.                   6/7=2,8=2,9=2,10=2/1;
  9. 1/18=40/14(2);
  10. 2/29=3/2;
  11. 99//99;
  12. 2/29=3/2;
  13. 3/11=9,25=1,30=1/1,2,3;
  14.                  4/5=5,16=3,69=1/1;
  15.                  5/5=2,38=5/2;
  16.                  8/6=4,9=120000,10=1/1,4;
  17.                  9/5=7,14=2/13;
  18. 1/18=40/14(-6);
  19. 2/29=3/2;
  20. 6/7=2,8=2,9=2,10=2/1;
  21. 99/9=10/99;
复制代码
但是写上#p opt(nomicro,ef,z-matrix) external='orca.sh',结果调用ORCA的时候关键词含有engard。
这是因为external使用ef方法用的是L113(功能:使用解析导数进行EF优化),所以给orca.sh传的第2个参数是1(要求一阶解析导数)
这些通过查看计算路径就可以看出来
  1. 1/18=40,26=1,29=20000,38=1,172=1/1,13;          多了26=1,这是个控制函数优化精度的选项,1相当于默认值
  2. 2/12=2,17=6,18=5,29=3,40=1/2;                一样
  3. 3/11=9,25=1,30=1,41=9900000,43=2,71=1/1;  多了后三条  41用9900000表示用external  43用2表示不考虑外部电荷(scrf)   71=1是导数识别标识(用于scrf)1阶导数
  4. 4/20=17,22=1,24=3,113=1,114=1/2;                    20=17表示用external而不是那些分子力学程序,22=1要求外部程序提供一阶导,24=3更新轨道特征值和其他一些东西
  5.    113=1 使用file 747 中的第 1 条命令 ,  114=1表示onionm的真实系统
  6. 6/7=2,8=2,9=2,10=2/1;   一样
  7. 7/44=-1/16;       多了L716  
  8. 1/18=40,26=1/13(2);     这是有解析梯度的EF
  9. 2/29=3/2;    一样
  10. 99//99;        一样
  11. 2/29=3/2;       一样
  12. 3/11=9,25=1,30=1,41=9900000,43=2,71=1/1;  和上一次调用L301一样
  13. 4/16=2,20=17,22=1,24=3,113=1,114=1/2;       比上一次调用L402多了16=2,这个选项控制是否转换从初猜中读取的波函数,2表示转换为当前原子坐标
  14.                  7/44=-1/16;                                               和前面一样
  15.                  1/18=40,26=1/13(-4);                                这是有解析梯度的EF
  16. 2/29=3/2;                       最后三个一样
  17. 6/7=2,8=2,9=2,10=2/1;              
  18. 99/9=10/99;
复制代码
sob老师的gau_orca是专门用来传递一阶解析导数和二阶解析导数,没有数值差分几何优化的功能。需要自己写代码添加(我写了一个样板,见附件,代码可能有冗余,.f90文件编译的extorca1111可执行程序是用"gfortran extorca1111.f90 -o  extorca1111"在Unbuntu20.04LTS下生成的,由于使用了动态库,所以文件比sob老师的extorca小了很多,在其他系统下可能无法正常使用

gaussian的external功能不太灵活,想用L114得费好大劲。
我先尝试用
  1. #p external='orca1111.sh'    nonstd
  2. .......
复制代码

结果是L1错误,因为用了nonstd关键词没法再自己加关键词,把external='orca.sh' 写在nonstd后面也没有用。
external='orca.sh'去掉则会导致L402无法调用orca.sh而失败,我真的有点受不了这破烂功能。
我想了好久
只好用一种更加奇怪的办法
  1. #p  external='orca1111.sh' iop(99/5=0,99/9=10) skip=ov9 extraoverlays

  2. 1/18=40,29=20000,38=1,172=1/1,14;
  3. 2/12=2,17=6,18=5,29=3,40=1/2;
  4. 3/11=9,25=1,30=1,41=9900000,43=2/1;
  5. 4/20=17,24=3,113=1,114=1/2;
  6. 6/7=2,8=2,9=2,10=2/1;
  7. 1/18=40/14(2);
  8. 2/29=3/2;
  9. 99//99;
  10. 2/29=3/2;
  11. 3/11=9,25=1,30=1,41=9900000,43=2/1;
  12. 4/20=17,24=3,113=1,114=1/2;
  13. 6/7=2,8=2,9=2,10=2/1;            
  14. 1/18=40/14(-4);
  15. 2/29=3/2;
  16. 6/7=2,8=2,9=2,10=2/1;
复制代码
4/113=1 使用file 747 中的第 1 条命令,我没找到这个文件,或许改掉这个文件的内容可以不用这么奇怪的办法。

我测试的时候用的是HF以缩减测试时间。
目前还有一点小问题,用这种方法得到的gaussian输出文件archive段落是空的,而ef数值差分优化,优化完成的几何结构不是在末尾,要翻到前面去找。考虑写个啥来显示一下最终结果。

3.一些未整合的脚本

gbwsave.sh
这是一个备份ORCA每一步SCF迭代结果的shell脚本
目前的版本写的不太好,有小bug,下次更新要修改的部分比较多,我由于不太熟练,只能不定期更新。
这个部分以后会增加,而且有朝一日要整合进一个脚本里


4.附件说明:
orca1111.sh和sob老师的orca.sh几乎一模一样,改动了写入mol.inp的内容,并且在关键词里写上了关闭了incremental Fock matrix的选项,在运行结束后不清理目录以便发现bug,把运行extorca改成了运行extorca1111
extorca1111.f90保留了extorca.f90原有的代码,添加了gaussian给"orca1111.sh"第二个参数传递0时,只从ORCA输出文件读取能量的代码。
由于不太会用fortran,提取ORCA输出文件的能量使用的是aaa.sh的shell脚本,shell脚本只需要一句就可以了,将ORCA的单点能直接写到Energy.out文件里。由于只有一个浮点数,其他冗余信息都没了,fortran直接read就完了。关于fortran定位字符串的实现,http://bbs.keinsci.com/thread-14848-1-1.html(《Fortran有快速识别文件内容的工具吗?》)sob老师的回答有详细描述。(但我太菜了,搞不定,直接复制代码就会报错,还要处理一下)
aaa.sh是extorca1111通过system函数调用命令行直接运行aaa.sh
虽然shell脚本有点多,但是这样也不算坏,因为编译.f90文件要额外花时间,debug还麻烦,一旦写错了得重新编译。shell脚本只要修改了就能用,你想读取其他程序产生的能量也很方便,灵活度高。只要把orca1111.sh把运行ORCA的命令换掉,aaa.sh把grep查找的内容换掉就可以了。至少读能量这种简单的数据,fortran的好处没有显示出来,甚至gaussian官方某个自带的external程序也是用shell脚本处理外部程序的输出的。
2020.8.4更新
新增ORCA保留scf每一步结果的独立脚本gbwsave此版本有小bug,不是特别严重,触发条件已经明确,具体情况见本贴开头。如果关键词没有“nofinalgrid”,则收敛后用更高精度的格点计算的结果不会备份。多步任务没有测试过,可能会出bug。


还能想到的,可能有用的方法(调用ORCA,甚至调用Gaussian都行)
1.用CCSD(T)/小基组-MP2/小基组+MP2/大基组,获得CCSD(T)/大基组近似能量,然后几何优化。MP2的解析导数在这里不知道用处大不大,最省事的方法就是全数值导数,反正CCSD(T)只能老老实实用数值导数。
2.CCSD在变量少的时候几何优化,Gaussian对于有解析梯度的理论方法强制使用解析梯度,由于L1110求CCSD梯度耗时远高于单点,变量少的时候完全可以考虑数值优化。
3.用能量基组外推,获得更高精度的能量,只需要简单数值运算即可,sob老师还有一个小程序enecalc,见 http://sobereva.com/308 (《Gaussian单点能自动相互运算工具enecalc》)。改一下拿来算这个也不错。

5.写脚本时遇到的困难(有几个我不敢拿定主意的地方)
1.对于几何优化,ORCA的输出文件可能变的很长,如果通过grep查找两个特定字符串("SCF ITERATIONS"和“”“SCF CONVERGED”)的行数来判断orca_scf是否在运行,会不会消耗较多时间?我目前是通过ps aux|||(一堆管道)来得知scf_orca是否在运行。
2.把迭代结果输出到.out和把迭代结果写入.gbw的先后顺序是怎么样的?会不会出现计算中途备份的.gbw是写了一半的.gbw?(如果有这种情况,读取写了一半的.gbw可能会导致波函数不收敛)。这个我没搞明白。
3.又查出一个莫名其妙的bug
  1. Run ORCA
  2. Please wait.   orca_scf (or orca_scf_mpi) is not running!
  3. orca_scf (or orca_scf_mpi) is running! Number of running is 1.(not ITERATIONS numbers)
  4. move reactant.gbw ****  ITERATIONS   0  ****
  5. move reactant.gbw ****  ITERATIONS   1  ****
  6. move reactant.gbw ****  ITERATIONS   2  ****
  7. move reactant.gbw ****  ITERATIONS   3  ****
  8. move reactant.gbw ****  ITERATIONS   4  ****
  9. move reactant.gbw ****  ITERATIONS   6  ****
  10. move reactant.gbw ****  ITERATIONS   7  ****
  11. move reactant.gbw ****  ITERATIONS   8  ****
  12. move reactant.gbw ****  ITERATIONS   9  ****
  13. move reactant.gbw ****  ITERATIONS  10  ****
  14. move reactant.gbw ****  ITERATIONS  11  ****
  15. move reactant.gbw ****  ITERATIONS  12  ****
  16. move reactant.gbw ****  ITERATIONS  13  ****
  17. move reactant.gbw ****  ITERATIONS  14  ****
  18. Please wait.   orca_scf (or orca_scf_mpi) is not running!
  19. ORCA running finished! ORCA terminated normally
复制代码

这是一个看起来还算正常的输出,但是如果加入
  1. echo $(tail -1 $output )   #这句是新加的
  2. result=$(tail -1 $output |cut -c 1-3)    #获取最后一行输出  这也改了一下,因为我一开始怀疑是"$(awk 'END {print}' $output)" 有问题
  3. result2=$(expr match "$result" "[0-9, ][0-9, ][0-9]"  )  #这句没变
复制代码


立马出现妖魔鬼怪 (把我这个目录里所有的文件名都给输出了)(红字部分是echo的正常输出)
Run ORCA
Please wait.   orca_scf (or orca_scf_mpi) is not running!
orca_scf (or orca_scf_mpi) is running! Number of running is 1.(not ITERATIONS numbers)
bug.docx bug.txt gbwsave.sh gbwsave.sh混乱了 gbwsaveDone2020-08-05_16:07:44 gbwsaveDone2020-08-05_16:34:27 gbwsaveNotfinished2020-08-05_16:36:11 gbwsave修改.sh gbwsave框架和一些额外注释.sh gbwsave第一版.sh part_test.sh reactant.E.tmp reactant.H.tmp reactant.K.tmp reactant.S.tmp reactant.T.tmp reactant.dummy.tmp reactant.engrad reactant.gbw reactant.ges reactant.inp reactant.int.tmp reactant.opt reactant.out reactant.prop reactant.xyz reactant_property.txt reactant_trj.xyz runall sudosu ~$bug.docx 新建文本文档.txt
# of grid points (after initial pruning) ... 51950 ( 0.0 sec)
Grid point division into batches done ... 0.5 sec
Now organizing SCF variables ... done

bug.docx bug.txt gbwsave.sh gbwsave.sh混乱了 gbwsaveDone2020-08-05_16:07:44 gbwsaveDone2020-08-05_16:34:27 gbwsaveNotfinished2020-08-05_16:36:11 gbwsave修改.sh gbwsave框架和一些额外注释.sh gbwsave第一版.sh part_test.sh reactant.E.tmp reactant.H.tmp reactant.K.tmp reactant.S.tmp reactant.SM12.tmp reactant.SP12.tmp reactant.T.tmp reactant.diis.tmp reactant.dummy.tmp reactant.engrad reactant.gbw reactant.ges reactant.grid.tmp reactant.inp reactant.int.tmp reactant.opt reactant.out reactant.prop reactant.xyz reactant_property.txt reactant_trj.xyz runall sudosu ~$bug.docx 新建文本文档.txt Starting incremental Fock matrix formation bug.docx bug.txt gbwsave.sh gbwsave.sh混乱了 gbwsaveDone2020-08-05_16:07:44 gbwsaveDone2020-08-05_16:34:27 gbwsaveNotfinished2020-08-05_16:36:11 gbwsave修改.sh gbwsave框架和一些额外注释.sh gbwsave第一版.sh part_test.sh reactant.E.tmp reactant.H.tmp reactant.K.tmp reactant.S.tmp reactant.SM12.tmp reactant.SP12.tmp reactant.T.tmp reactant.diis.tmp reactant.dummy.tmp reactant.engrad reactant.gbw reactant.ges reactant.grid.tmp reactant.inp reactant.int.tmp reactant.opt reactant.out reactant.prop reactant.xyz reactant_property.txt reactant_trj.xyz runall sudosu ~$bug.docx 新建文本文档.txt
0 -690.0547450139 0.000000000000 0.13747907 0.00989124 0.4483969 0.7000
move reactant.gbw ****  ITERATIONS   0  ****
1 -690.1937901532 -0.139045139222 0.10276143 0.00678357 0.1734227 0.7000
move reactant.gbw ****  ITERATIONS   1  ****
2 -690.2390744860 -0.045284332811 0.11312586 0.00872183 0.0285016 0.0000
move reactant.gbw ****  ITERATIONS   2  ****
3 -690.3219913604 -0.082916874455 0.07690480 0.00480596 0.0996510 0.0000
move reactant.gbw ****  ITERATIONS   3  ****
4 -690.3369223803 -0.014931019892 0.01777994 0.00130864 0.0121746 0.0000
move reactant.gbw ****  ITERATIONS   4  ****
(重复前面的那一长串)Restarting incremental Fock matrix formation(重复前面的那一长串)
move reactant.gbw ****  ITERATIONS   6  ****
7 -690.33789890 -0.0000217078 0.000205 0.001757 0.001472 0.000150
move reactant.gbw ****  ITERATIONS   7  ****
8 -690.33789915 -0.0000002467 0.000243 0.000552 0.000807 0.000060
move reactant.gbw ****  ITERATIONS   8  ****
10 -690.33790082 0.0000000386 0.000066 0.000089 0.000230 0.000014
move reactant.gbw ****  ITERATIONS  10  ****
11 -690.33790093 -0.0000001064 0.000006 0.000012 0.000033 0.000002
move reactant.gbw ****  ITERATIONS  11  ****
13 -690.33790093 -0.0000000014 0.000001 0.000001 0.000004 0.000000
move reactant.gbw ****  ITERATIONS  13  ****
14 -690.33790093 0.0000000000 0.000001 0.000001 0.000002 0.000000
move reactant.gbw ****  ITERATIONS  14  ****
Grid point re-assignment to atoms done ... 0.1 sec
Final grid set up in 4.1 sec
Total SCF time: 0 days 0 hours 0 min 52 sec

Please wait.   orca_scf (or orca_scf_mpi) is not running!
ORCA running finished! ORCA terminated normally
第五次迭代也出事了。
不管是tail还是awk都是一模一样的错误
我现在非常非常懵逼。输入文件reactant.inp已经传上来了。(为了缩减测试时间,基组用了个没法得到有意义的6-31g)

大概猜出来是什么东西的问题了,去掉一部分输出,直接输命令行就是。
tail -1 reactant.out
               *** Restarting incremental Fock matrix formation ***
万恶之源肯定是***,不知道在shell脚本里这***被搞成了什么鬼样子。
但不管怎么说,这套代码会漏掉几次收敛,不过一般也就漏一次,也不会连续漏。如果是因为每一圈太快而漏,那也没太大必要保留中间波函数,直接算也费不了多少时间。

aaa.sh

108 Bytes, 阅读权限: 10, 下载次数 Times of downloads: 6

extorca1111.f90

1.62 KB, 阅读权限: 10, 下载次数 Times of downloads: 11

extorca1111

21.33 KB, 阅读权限: 10, 下载次数 Times of downloads: 3

orca1111.sh

1.03 KB, 阅读权限: 10, 下载次数 Times of downloads: 7

H2-cc-tz.gjf

583 Bytes, 阅读权限: 10, 下载次数 Times of downloads: 5

gbwsave.sh

6.28 KB, 阅读权限: 10, 下载次数 Times of downloads: 8

reactant.inp

595 Bytes, 阅读权限: 10, 下载次数 Times of downloads: 0

评分 Rate

参与人数
Participants 7
eV +32 收起 理由
Reason
li447fan + 3 好物!
哇哇吐 + 5 赞!
喵星大佬 + 3 牛!
thanhtam + 3
ggdh + 5 太强了!
sobereva + 10 精品内容
shalene + 3 你太可爱

查看全部评分 View all ratings

877

帖子

36

威望

4805

eV
积分
6402

Level 6 (一方通行)

发表于 Post on 2020-8-9 12:25:20 | 显示全部楼层 Show all
弱弱的问下 你的体系试过直接用orca的opt优化没有?
一般能用上CCSD(T)的数值优化的体系应该非常小把
这张情况下不同的opt算法之间的误差还会有这么大么

112

帖子

1

威望

933

eV
积分
1065

Level 4 (黑子)

 楼主 Author| 发表于 Post on 2020-8-11 00:51:26 | 显示全部楼层 Show all
本帖最后由 ldatea 于 2020-8-11 09:45 编辑
ggdh 发表于 2020-8-9 12:25
弱弱的问下 你的体系试过直接用orca的opt优化没有?
一般能用上CCSD(T)的数值优化的体系应该非常小把
这 ...

讲真,现在回过头来看这个确实不好用。DLPNO的加速效果在分子小的时候没明显优势,使用ORCA必须放弃对称性加速,这挺尴尬的。
Gaussian在保持对称性上做的不错,有对称性的体系可以通过进一步减少优化的变量数来降低时间消耗。但是在ORCA上不能利用对称性加速,致命一击。


beefly老师提到的虚原子的坐标还会改变的问题,已经可以解决了。
我不清楚ORCA是不是只能做全优化,看beefly老师的意思我感觉是没法冻结变量Gaussian的优化更容易上手,在收敛限的设置上比ORCA稳健。
这个工具最大的用处大概就是帮助对Gaussian有足够了解,但对ORCA很不熟悉的新人了。
写这些的主要目的还是学习写shell脚本,以及我比较懒暂时不想看ORCA的手册。


本版积分规则 Credits rule

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

GMT+8, 2023-2-7 03:47 , Processed in 0.264373 second(s), 25 queries .

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