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

ldatea

本文会根据需要，持续更新注意:附件的程序使用的是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. 
   
3. 1/18=40,29=20000,38=1,172=1/1,14;
   
4. 2/12=2,17=6,18=5,29=3,40=1/2;
   
5. 3/11=9,25=1,30=1,41=9900000,43=2/1;
   
6. 4/20=17,24=3,113=1,114=1/2;
   
7. 6/7=2,8=2,9=2,10=2/1;
   
8. 1/18=40/14(2);
   
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/1;
    
13. 4/20=17,24=3,113=1,114=1/2;
    
14. 6/7=2,8=2,9=2,10=2/1;            
    
15. 1/18=40/14(-4);
    
16. 2/29=3/2;
    
17. 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脚本里这***被搞成了什么鬼样子。
但不管怎么说，这套代码会漏掉几次收敛，不过一般也就漏一次，也不会连续漏。如果是因为每一圈太快而漏，那也没太大必要保留中间波函数，直接算也费不了多少时间。

ggdh

弱弱的问下 你的体系试过直接用orca的opt优化没有？
一般能用上CCSD(T)的数值优化的体系应该非常小把
这张情况下不同的opt算法之间的误差还会有这么大么

ldatea

ggdh 发表于 2020-8-9 12:25
弱弱的问下 你的体系试过直接用orca的opt优化没有？
一般能用上CCSD(T)的数值优化的体系应该非常小把
这 ...

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


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