使用G16的external调用G16（无限套娃）

ldatea

目录
1.开发此功能想法的来源
2.附件说明
3.目前能实现的功能
4.一些细节的说明
5.目前遇到的困难

1.开发此功能想法的来源
         前段时间修改了一下sob老师的gau_orca   http://bbs.keinsci.com/thread-10141-1-1.html（《将Gaussian与ORCA联用搜索过渡态、产生IRC、做振动分析》），实现一些新的功能http://bbs.keinsci.com/thread-18782-1-1.html（《给sob老师的gau_orca添加额外的功能》），但是新的功能用处不大。这一次其实算是上一次的延续，由于不涉及ORCA，所以另开一贴。虽然没有使用ORCA，但是所有文件都是基于gau_orca修改而来的。这个模板可以说具有普遍意义，sob老师的拓荒非常关键。
     事实上，gaussian的IOP和external功能很强大，可以完成一些非常规操作。但是使用起来不是很方便，需要测试好多次，一些bug可能藏得很深。
     另外参考了g16自带的脚本extuffex（使用external调用MM程序做uff力场几何优化）题外话：这个脚本在win版的G09里也有，不管win还是Linux还有实现同样功能的.bat脚本，但是我在windows下用这个.bat没成功过，尝试修改也没试出来
这个脚本最奇怪的地方就是完全是sh脚本，与fortran没有关系，使用awk的printf(%20.12e格式）“模仿” 3D20.12或者4D20.12(具体细节，请在Gaussian安装目录下的bsd文件夹里面找。之前写错了，写成tests了，已经改成bsd)。最后的浮点数的指数是E+xx或者E-xx。Gaussain不管是D还是E的浮点数都是“认识”的。这实际上有不少好处，只需要一个shell脚本就可以完成所有的事，简单的数据提取shell脚本也不会比fortran难用，就是代码可读性差一些，一些提取字符串的功能不那么容易看懂。如果要做数学运算，shell脚本虽然也能实现，但是可读性就更差了。由于本人水平有限，目前仍然使用一个shell脚本和一个Fortran可执行文件。

2.附件说明
压缩包含有gau_gau.sh    ext_gau.f90  ext_gau    examples
gau_gau.sh由sob老师的orca.sh修改而来
ext_gau.f90由sob老师的ext_orca.sh修改而来
ext_gau.f90文件编译的可执行程序 ext_gau是用"gfortran extorca1111.f90 -o  extorca1111"在Unbuntu20.04LTS下生成的，由于使用了动态库，在其他系统下可能无法正常使用
examples含有一个例子，包含输入文件和部分输出文件，内容是用MP2.5/aug-cc-pvtz 优化H2—H2 (T构型)的弱相互作用体系。

3.目前能实现的功能
目前可以应用于有解析梯度的两种理论方法按一定比例组合的情况，通过分别计算解析梯度再按比例组合，来做几何优化。

注意：MP2.X以及SCS-MP3目前仅能做到可用(附件的功能只是MP2.5，请根据需要修改)，但是功能很烂，算是个半成品，很多时候这么做可能还不如选择合适的双杂化泛函。目前只是做到能分别调用g16计算MP2,MP3的能量和解析梯度，相当于MP2算了两遍。我接下来会尝试通过IOP减少没有意义的Link调用。

MP2.5是sob老师在http://sobereva.com/344（《Gaussian中非内置的理论方法和泛函的用法》）中提到的方法。这个方法显然没法直接做几何优化，但是gaussian可以分别计算MP2和MP3的解析梯度，所以需要external功能，就可以实现几何优化甚至振动分析了(MP3在Gaussian中没有解析Hessian,实现起来可能非常慢）不过我找了一个最简单的H2分子来看几何优化的效果。
此处MP2和MP2.5都很烂
注意：氢的弥散函数大多数情况下是没有必要的。此处使用的基组不一定合适，格点可能没必要这么多。本身对氢气几何优化没有必要加弥散函数，
这里用了弥散函数只是拿来对比一下。反正氢气特别小，随便折腾。
但是，如果要研究弱相互作用，弥散函数是有必要的，用多少层弥散函数也是需要考虑的。
（aug-cc-pvtz给氢加了S,P,D一共三层弥散函数）
H2 CCSD(T)/aug-cc-pv6z             0.74151  verytight
H2 CCSD(T)/cc-pv6z                 0.74148  verytight
H2 CCSD(T)/cc-pv5z                 0.74154  verytight
H2 CCSD(T)/cc-pvqz                 0.74187  verytight
H2 CCSD(T)/cc-pvtz                 0.74263  verytight
H2 CCSD(T)/aug-cc-pvtz             0.74299  verytight
H2 MP2/aug-cc-pvtz                 0.73744  tight
H2 MP2.5/aug-cc-pvtz               0.73864  tight
H2 DSD-PBEP86-D3(BJ)/aug-cc-pvtz   0.74215  (verytight)tight grid=400974
H2 B2PLYP/aug-cc-pvtz              0.74003  (verytight)tight grid=400974
H2 B2PLYP-D3(BJ)/aug-cc-pvtz       0.74005  (verytight)tight grid=400974

不过MP2.5之类还可以算弱相互作用，这个我还没测试过。(以后会简单测试一下，弱相互作用我没有经验，只能拿一些已经研究过的体系)
注意：MP3的解析梯度耗时远高于MP2的解析梯度,我算个H2—H2（T构型）的几何优化，MP3的解析梯度占的时间>90%。external默认的几何优化算法没有原版的MP2算法好，几何优化走的步数比较多，我将尝试改进。（这里纯粹测试优化算法，把ext_gau里面的WeightedNum调成0.，不过这个时候MP3就没必要算了，可以用任意一个PESMP3.out来代替，并且把gau_gau.sh中的MP3的部分注释掉）。

4.一些细节的说明

IOP(7/32=x)可以在L716运行结束后，以多种格式输出核坐标，能量，受力和hessian。gau_gau.sh使用的是IOP(7/32=3),由于输出是受力，而external要求一阶导数时，需要的是梯度，两者相差一个负号。ORCA的engrad关键词直接输出梯度，所以不需要变号。而Gaussian输出的总是受力，所以要加负号（ext_gau.f90梯度的计算都加了负号）
仔细看Gaussian自带的extuffex，也可以发现awk提取梯度时，带了一个负号。

G16可以在使用external功能时，同一个目录下的g16可执行文件可以同时处理两个输入文件，我猜只要两个任务不同时进同一个Link就没啥问题，external运行时，原始的任务总是会停在L402，L402不使用external或者MM程序时根本不会调用，不需要担心两个任务共存的问题，没有必要再复制一份新的G16。但是如果再运行一个新任务，则会报错。

5.目前遇到的困难
MP2和MP3在L801之前的计算路径是一模一样的，L801以后，MP3进L804再进L913,而MP2直接进L906。L804可以输出与L906一样的信息，包括MP2的E2，以及分别输出alpha-alpha,alpha-beta,beta-beta电子对E2的贡献。但是用把L906改成L804会导致MP2 force任务在L1002报错。MP3也会遇到类似的问题，替换以后也是L1002报错。如果路径中即有L804也有L906，则两次输出MP2的能量是不一样的（第一次是对的，第二次是错的）。
如果这两个都得分别调用，那么MP2.5通过IOP减少耗时的意义就几乎没有了。

ggdh

能不能把这个脚本继续开发一下，实现使用gaussian 调用xtb 计算ONIOM啊。参考sob的这个帖子。。
http://bbs.keinsci.com/thread-10106-1-1.html

ldatea

ggdh 发表于 2020-8-14 11:49
能不能把这个脚本继续开发一下，实现使用gaussian 调用xtb 计算ONIOM啊。参考sob的这个帖子。。
http://bb ...

不好意思，我不懂oniom，交界的地方需要怎么补足原子我不清楚。请给一个例子，并详细描述需求。

liuyuje714

不错的东西。如果不考虑awk数据精度的限制和速度，直接就没必要用这个fortran接口了。把脚本中的最后那个./ext_gau $3 $atoms $derivs直接换成下面的代码即可。当然对于这种较高精度运算，用fortran或者C等写接口才靠谱。

1. WeightedNum=0.500000000
   
2. awk '
   
3.     BEGIN {
   
4.         natom = "'$atoms'"+0; WeightedNum = "'$WeightedNum'"+0
   
5.     }
   
6.    
   
7.     NR == FNR {
   
8.         if(NR == 1) {mp2ene = $1}
   
9.         if(NR > 2) {
   
10.             for(i = 1; i <= natom; i++) {
    
11.                 gsub("D", "E", $1); fx2[i] = $1; getline
    
12.                 gsub("D", "E", $1); fy2[i] = $1; getline
    
13.                 gsub("D", "E", $1); fz2[i] = $1; getline
    
14.             }
    
15.         }
    
16.     }
    
17. 
    
18.     NR != FNR {
    
19.         if(FNR == 1) {mp3ene = $1}
    
20.         if(FNR > 2) {
    
21.             for(i = 1; i <= natom; i++) {
    
22.                 gsub("D", "E", $1); fx3[i] = $1; getline
    
23.                 gsub("D", "E", $1); fy3[i] = $1; getline
    
24.                 gsub("D", "E", $1); fz3[i] = $1; getline
    
25.             }
    
26.         }
    
27.     }
    
28. 
    
29.     END {
    
30.         enesum=(1-WeightedNum)*mp2ene+WeightedNum*mp3ene
    
31.         printf("%20.12e%20.12e%20.12e%20.12e\n", enesum, 0, 0, 0)
    
32.         for(j = 1; j <= natom; j++) {
    
33.             fxsum=-(1-WeightedNum)*fx2[j]-WeightedNum*fx3[j]
    
34.             fysum=-(1-WeightedNum)*fy2[j]-WeightedNum*fy3[j]
    
35.             fzsum=-(1-WeightedNum)*fz2[j]-WeightedNum*fz3[j]
    
36.             printf("%20.12e%20.12e%20.12e\n", fxsum, fysum, fzsum)
    
37.         }
    
38.     }
    
39.    
    
40. ' cfmp2.out cfmp3.out > $3

复制代码

ldatea

liuyuje714 发表于 2020-8-14 12:42
不错的东西。如果不考虑awk数据精度的限制和速度，直接就没必要用这个fortran接口了。把脚本中的最后那个./ ...

非常感谢大佬优化脚本，
我awk命令都是网上东拼西凑学的，现在只会最基础的简单功能，把for和if用在一起就不会了。
之前还一直没搞定awk读取含字母D的浮点数，本来按extuffex写也是可以的，为了稍微提高精度没有选择读.out里面的笛卡尔坐标下的力，而是iop(7/32=3)输出PES.out。awk直接读浮点数只认识含E的浮点数。我最后还是没搞成功，现在算是解决了。

Daniel_Arndt

ldatea 发表于 2020-8-14 13:17
非常感谢大佬优化脚本，
我awk命令都是网上东拼西凑学的，现在只会最基础的简单功能，把for和if用在一起 ...

我以前写awk脚本读取Gaussian 09的BOMD任务的代码中有一段如下：
/^ Total angular / {
    m = index($4,"D");
    j = substr($4,1,m - 1);
    s = substr($4,m + 1,1);
    j1 = substr($4,m + 2,1);
    j2 = substr($4,m + 3);
    if (j2 > 9)
    {
        printf("\nERROR1!\n");
        exit 1;
    }
    j2 = j1 * 10 + j2;
    if (s == "+")
        j = j * (10 ** j2);
    else if (s == "-")
        j = j / (10 ** j2);
    else
    {
        printf("\nERROR2!\n");
        exit 1;
    }
    printf("total angular momentum: %f h-bar\n", j);
}


后来社长说直接用字符串替换将“D”替换成“E”，就可以用awk读取了，不用像我这样手动写个parser。。。。。。

其实gawk从4.1.0版本开始支持MPFR，就是在调用awk脚本时加一个“-M”的参数，可以支持任意精度的浮点数，而不是默认的双精度的浮点数。目前只要你的Linux操作系统不是太老，默认的awk应该是gawk，而且肯定高于4.1.0版。但是有一个操作系统需要注意，openSUSE，我印象里它的gawk版本高于4.1.0，但是编译时没有支持MPFR。gawk里面的MPFR的具体使用方法可以参考 http://gawkextlib.sourceforge.net/mpfr/gawk-mpfr.html 。GNU官网上也可以查到一些内容。

针对你的情况，如果想用awk的话，可能还需要使用一个扩展函数：strtonum。你可以参考这个 https://byteofbio.com/archives/1.html 。

snljty

个人建议，遇到D还是用fortran写个小程序读一下，毕竟主流语言对D代表双精度这种(奇怪的)格式支持最好的也只有Fortran了。

ldatea

Daniel_Arndt 发表于 2020-8-18 01:52
我以前写awk脚本读取Gaussian 09的BOMD任务的代码中有一段如下：
/^ Total angular / {
    m = index( ...

感谢有价值的补充，写脚本免不了和一些棘手的问题纠缠，兼顾各方面的话，难度还是有的。有空我去看看python，不知道python在处理这种需求的时候有没有优秀之处。

ldatea

snljty 发表于 2020-8-18 09:59
个人建议，遇到D还是用fortran写个小程序读一下，毕竟主流语言对D代表双精度这种(奇怪的)格式支持最好的也 ...

求稳妥肯定是fortran，如果是经常要把脚本改来改去可能还是解释型语言的脚本比较方便。

lanthanum

在运行原始的test0726.gjf文件时，通过extuffex.bat脚本调用了高斯内置的uff分子力学方法，顺利执行成功。现在将关键词改成#p opt=nomicro nosymm external="C:\g16w\mm.exe -uff -force -external" ，
希望调用高斯工具包中的mm.exe utility实现同样功能，未能执行。输出文件见附件。
请问老师：
（1）如何写关键词，才能执行？
（2）mm.exe utility的输入和输出文件格式是什么？ 若用-external参数后，输入和输出文件的格式有什么变化？
谢谢。

Freeman

您好。我是fortran小白。请问3D20.12和4D20.12有区别吗？他们都可以用awk的%20.12e格式模仿吗？

sobereva

Freeman 发表于 2021-12-16 20:55
您好。我是fortran小白。请问3D20.12和4D20.12有区别吗？他们都可以用awk的%20.12e格式模仿吗？

3D20.12是按D20.12输出三个值，4D20.12是输出4个
这么写的时候Fortran会用D代表科学计数法的E