计算化学公社

 找回密码 Forget password
 注册 Register
Views: 11741|回复 Reply: 11

[Gaussian/gview] 使用G16的external调用G16(无限套娃)

[复制链接 Copy URL]

112

帖子

1

威望

933

eV
积分
1065

Level 4 (黑子)

发表于 Post on 2020-8-13 01:10:21 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 ldatea 于 2020-8-24 18:50 编辑

目录
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减少耗时的意义就几乎没有了。

















gau_gau.7z

432.9 KB, 阅读权限: 10, 下载次数 Times of downloads: 31

评分 Rate

参与人数
Participants 5
威望 +1 eV +20 收起 理由
Reason
北大-陶豫 + 5 好物!
snljty + 5
ggdh + 5 不明觉厉
小范范1989 + 5 好物!
sobereva + 1

查看全部评分 View all ratings

877

帖子

36

威望

4805

eV
积分
6402

Level 6 (一方通行)

发表于 Post on 2020-8-14 11:49:25 | 显示全部楼层 Show all
本帖最后由 ggdh 于 2020-8-14 11:51 编辑

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

112

帖子

1

威望

933

eV
积分
1065

Level 4 (黑子)

 楼主 Author| 发表于 Post on 2020-8-14 12:18:02 | 显示全部楼层 Show all
ggdh 发表于 2020-8-14 11:49
能不能把这个脚本继续开发一下,实现使用gaussian 调用xtb 计算ONIOM啊。参考sob的这个帖子。。
http://bb ...

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

183

帖子

4

威望

1457

eV
积分
1720

Level 5 (御坂)

发表于 Post on 2020-8-14 12:42:04 | 显示全部楼层 Show all
本帖最后由 liuyuje714 于 2020-8-14 15:05 编辑

不错的东西。如果不考虑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.     NR != FNR {
  18.         if(FNR == 1) {mp3ene = $1}
  19.         if(FNR > 2) {
  20.             for(i = 1; i <= natom; i++) {
  21.                 gsub("D", "E", $1); fx3[i] = $1; getline
  22.                 gsub("D", "E", $1); fy3[i] = $1; getline
  23.                 gsub("D", "E", $1); fz3[i] = $1; getline
  24.             }
  25.         }
  26.     }

  27.     END {
  28.         enesum=(1-WeightedNum)*mp2ene+WeightedNum*mp3ene
  29.         printf("%20.12e%20.12e%20.12e%20.12e\n", enesum, 0, 0, 0)
  30.         for(j = 1; j <= natom; j++) {
  31.             fxsum=-(1-WeightedNum)*fx2[j]-WeightedNum*fx3[j]
  32.             fysum=-(1-WeightedNum)*fy2[j]-WeightedNum*fy3[j]
  33.             fzsum=-(1-WeightedNum)*fz2[j]-WeightedNum*fz3[j]
  34.             printf("%20.12e%20.12e%20.12e\n", fxsum, fysum, fzsum)
  35.         }
  36.     }
  37.    
  38. ' cfmp2.out cfmp3.out > $3
复制代码


评分 Rate

参与人数
Participants 1
eV +5 收起 理由
Reason
ldatea + 5 好物!

查看全部评分 View all ratings

112

帖子

1

威望

933

eV
积分
1065

Level 4 (黑子)

 楼主 Author| 发表于 Post on 2020-8-14 13:17:00 | 显示全部楼层 Show all
liuyuje714 发表于 2020-8-14 12:42
不错的东西。如果不考虑awk数据精度的限制和速度,直接就没必要用这个fortran接口了。把脚本中的最后那个./ ...

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

390

帖子

1

威望

2032

eV
积分
2442

Level 5 (御坂)

发表于 Post on 2020-8-18 01:52:34 | 显示全部楼层 Show all
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

1188

帖子

5

威望

2758

eV
积分
4046

Level 6 (一方通行)

发表于 Post on 2020-8-18 09:59:27 | 显示全部楼层 Show all
个人建议,遇到D还是用fortran写个小程序读一下,毕竟主流语言对D代表双精度这种(奇怪的)格式支持最好的也只有Fortran了。

112

帖子

1

威望

933

eV
积分
1065

Level 4 (黑子)

 楼主 Author| 发表于 Post on 2020-8-18 13:21:30 | 显示全部楼层 Show all
Daniel_Arndt 发表于 2020-8-18 01:52
我以前写awk脚本读取Gaussian 09的BOMD任务的代码中有一段如下:
/^ Total angular / {
    m = index( ...

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

112

帖子

1

威望

933

eV
积分
1065

Level 4 (黑子)

 楼主 Author| 发表于 Post on 2020-8-18 13:23:06 | 显示全部楼层 Show all
snljty 发表于 2020-8-18 09:59
个人建议,遇到D还是用fortran写个小程序读一下,毕竟主流语言对D代表双精度这种(奇怪的)格式支持最好的也 ...

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

52

帖子

0

威望

2402

eV
积分
2454

Level 5 (御坂)

发表于 Post on 2021-4-5 00:05:07 | 显示全部楼层 Show all
本帖最后由 lanthanum 于 2021-4-5 00:06 编辑

在运行原始的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参数后,输入和输出文件的格式有什么变化?
谢谢。

test0726_extmm_1.out

43.52 KB, 下载次数 Times of downloads: 0

197

帖子

5

威望

1689

eV
积分
1986

Level 5 (御坂)

发表于 Post on 2021-12-16 20:55:09 | 显示全部楼层 Show all
您好。我是fortran小白。请问3D20.12和4D20.12有区别吗?他们都可以用awk的%20.12e格式模仿吗?

4万

帖子

99

威望

4万

eV
积分
89975

管理员

公社社长+计算化学玩家

发表于 Post on 2021-12-24 02:14:03 | 显示全部楼层 Show all
Freeman 发表于 2021-12-16 20:55
您好。我是fortran小白。请问3D20.12和4D20.12有区别吗?他们都可以用awk的%20.12e格式模仿吗?

3D20.12是按D20.12输出三个值,4D20.12是输出4个
这么写的时候Fortran会用D代表科学计数法的E
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

本版积分规则 Credits rule

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

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

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