计算化学公社
标题: 高性价比热力学组合方法G4(MP2)-6X的计算方法 [打印本页]
作者Author: sobereva 时间: 2018-5-27 00:04
标题: 高性价比热力学组合方法G4(MP2)-6X的计算方法
高性价比热力学组合方法G4(MP2)-6X的计算方法
Cost-effective thermodynamics composite method G4(MP2)-6X
文/Sobereva @北京科音 First release: 2018-May-26 Last update: 2019-Jul-30
在2011年，J. Chem. Theory Comput., 7, 112-120中，几个非Gn系列方法的官方人员提出了G4(MP2)-6X热力学组合方法，号称耗时与G4(MP2)相仿佛，而精度很接近昂贵得多的G4。其相对于G4(MP2)最主要变化是把MP2换成了SCS-MP2，把CCSD(T)的CCSD和(T)部分的相关能乘了系数，把优化和振动分析的泛函从B3LYP改为BMK。此方法竟然一直到当下最新的G16 B.01里都仍然没有被加入，好在文章的补充材料里给出了基于Gaussian做此方法计算的Perl脚本，使用很简单，这里简单说一下用法。原作者给的脚本只能输出H(0)和H(T)，没法输出指定温度下的自由能和内能，因此我对脚本进行了一些修改使之能够输出。
首先需要编辑模板.gjf文件，内容如下，也可以直接从这里下：
(, 下载次数 Times of downloads: 93)
。需要将此文件中的坐标、电荷和自旋多重度改成自己分子的情况
%chk=mol.chk
# BMK/6-31+G(2df,p) Opt
A molecule G4(MP2)-6X calculation
0 1
C 0.00000000 0.00000000 -0.56221066
H 0.00000000 -0.92444767 -1.10110537
H -0.00000000 0.92444767 -1.10110537
O 0.00000000 0.00000000 0.69618930
--Link1--
%chk=mol.chk
# Geom=AllCheck Guess=Read BMK/6-31+G(2df,p) Freq
--Link1--
%chk=mol.chk
# Geom=AllCheck Guess=Read CCSD(T,FrzG4)/GTBas1
--Link1--
%chk=mol.chk
# Geom=AllCheck Guess=Read MP2(FrzG4)/GTMP2LargeXP
--Link1--
%chk=mol.chk
# Geom=AllCheck Guess=Read HF/GFHFB3
--Link1--
%chk=mol.chk
# Geom=AllCheck Guess=Read HF/GFHFB4
用Gaussian运行此脚本，产生比如G4MP2_6x.out。G09和G16经测试都可以用。
下载此Perl脚本：
(, 下载次数 Times of downloads: 81)
。把Gaussian输出文件文件和这个.pl文件都拷到Linux下，运行G4MP2_6x.pl G4MP2_6x.out，这个Perl脚本就会自动把相关数据从Gaussian输出文件中提取出来并进行处理，默认是在标况下算的。输出信息例子如下
Temperature (K) 298.15
Pressure (atm) 1
NImag 0
E_ele -114.40450436
H(0K) -114.37806386
H(T) -114.37424742
U(T) -114.37519161
G(T) -114.39906286
Nimag就是虚频数目，E_ele就是电子能量，其它的都不言自明，单位是Hartree。
如果需要计算别的温度和大气压的情况，用文本编辑器打开.pl脚本，修改开头的$temp和$pres即可。
注：后来此方法的作者又提出了G4(MP2)-XK，把G4(MP2)-6X用的Pople基组改为了def2系列，使得此方法可以用于H~Rn的主族体系，对前四周期精度和G4(MP2)-6X相仿佛。在其原文DOI: 10.1021/acs.jctc.9b00449的补充材料里提供了相应的结合Gaussian使用的Perl脚本。
作者Author: wuzhiyi 时间: 2019-7-25 15:35
谢谢sob老师的分享 有一问题请教 我的小分子液相和气相优化结构有较大的差异 想问能否用热力学的方式来算？
我的想法是这样
%chk=mol.chk
# BMK/6-31+G(2df,p) Opt SCRF=SMD
A molecule G4(MP2)-6X calculation
0 1
C 0.00000000 0.00000000 -0.56221066
H 0.00000000 -0.92444767 -1.10110537
H -0.00000000 0.92444767 -1.10110537
O 0.00000000 0.00000000 0.69618930
--Link1--
%chk=mol.chk
# Geom=AllCheck Guess=Read BMK/6-31+G(2df,p) Freq SCRF=SMD
--Link1--
%chk=mol.chk
# Geom=AllCheck Guess=Read CCSD(T,FrzG4)/GTBas1
--Link1--
%chk=mol.chk
# Geom=AllCheck Guess=Read MP2(FrzG4)/GTMP2LargeXP
--Link1--
%chk=mol.chk
# Geom=AllCheck Guess=Read HF/GFHFB3
--Link1--
%chk=mol.chk
# Geom=AllCheck Guess=Read HF/GFHFB4
--Link1--
%chk=mol.chk
# Geom=AllCheck Guess=Read M052X/6-31G*
--Link1--
%chk=mol.chk
# Geom=AllCheck Guess=Read M052X/6-31G* SCRF=SMD
最后加上溶解自由能
作者Author: sobereva 时间: 2019-7-25 21:18
不是“热力学的方式”而是“热力学组合方法的方式”
这么做从实际角度上可行，但原理上不十分严格，毕竟热力学组合方法里的参数是对于气相优化的情况拟合的
作者Author: wuzhiyi 时间: 2019-7-30 16:02
本帖最后由 wuzhiyi 于 2019-7-30 17:05 编辑
https://pubs.acs.org/doi/10.1021/acs.jctc.9b00449
这里好像有人提出了一个将G4(MP2)-6X拓展到第五行以内的所有的主族元素
作者说在有些测试集表现比原本的G4(MP2)-6X好
作者Author: wuzhiyi 时间: 2019-9-18 19:12
想请教社长一个比较细节的问题，对于HLC的计算，作者分成了三种情况考虑，闭壳，开壳，原子和单电子对。
在作者的perl脚本中，为了区分这几种情况，作者考虑了四个变量$na，$mult，$nh和$ex
$nh和$ex是根据$xyz[2]这个变量生产的
在gaussian输出中对应Atomic Type
在我的输出中永远是0.
我想问一下这Atomic Type到底是啥，我翻manual好像是onion的东西，但肯定这里做的是纯QM的，和onion无关。
作者Author: sobereva 时间: 2019-9-19 16:25
量化计算这个始终都是0，不用管（对于Bq原子是1000）
是ONIOM不是ONION
作者Author: wuzhiyi 时间: 2019-9-19 16:45
谢谢社长 我发现在perl中$xyz[2]对应的是第二个值不是第三个值
所以这应该对应原子序号
作者Author: zmm0418 时间: 2019-10-8 17:08
请问社长，想用这个方法计算单点能，不知道精度怎么样呀？
作者Author: sobereva 时间: 2019-10-9 08:37
精度很理想
作者Author: 1355447622 时间: 2020-5-13 09:26
本帖最后由 1355447622 于 2020-5-13 09:29 编辑
sobereva 老师您好，我想用这个方法计算一个单独硫原子的热力学数据，但是算完后我在Liunx里运行脚本时他没有数据，我算的其他分子倒是挺正常的，请问怎么办呢？
C:\Users\admin\Desktop\233[scu@localhost ~]$ cd wzk/[scu@localhost wzk]$ perl G4MP2_6x.pl S.out
Can't take log of 0 at G4MP2_6x.pl line 298, <LOG> line 1749.
S.out [scu@localhost wzk]$
作者Author: linq5203 时间: 2020-5-13 22:27
本帖最后由 linq5203 于 2020-5-13 23:26 编辑
1355447622 发表于 2020-5-13 09:26
sobereva 老师您好，我想用这个方法计算一个单独硫原子的热力学数据，但是算完后我在Liunx里运行脚本时他没 ...
单原子连频率都没有，当然报错。详情可看Shermo手册或自行复习物化（不是阉割版的话）
作者Author: sobereva 时间: 2020-5-14 15:58
单个原子虽然能做振动分析（此时只有平动贡献），但没法优化，显然脚本调用Gaussian时会报错。你可以修改脚本，把优化的步骤去掉。
作者Author: 一条君 时间: 2020-7-4 11:36
请问老师此方法如何在windows下输出能量呢？谢谢
作者Author: sobereva 时间: 2020-7-5 02:03
何必非要用windows，本来这种耗时高的任务就严重不适合windows版Gaussian
不想装实体linux就用虚拟机就完了
作者Author: Melvin 时间: 2022-5-25 11:47
sobereva老师请问这个“可以修改脚本，把优化的步骤去掉”这个怎么实现呢（比如用的CBS-QB3），是在gjf文件里加什么关键字就可以实现吗？
作者Author: sobereva 时间: 2022-5-26 05:44
自己把perl脚本内容看一遍就知道改哪
当前脚本是做G4(MP2)-6X的，跟CBS-QB3没关系。你若要用CBS-QB3而且不优化，在G16里写CBS-QB3=noopt就完了。
作者Author: Melvin 时间: 2022-5-26 14:18
本帖最后由 Melvin 于 2022-5-26 14:30 编辑
受教了，感谢Sob老师，还有就是noopt只能G16实现吗，G09这个指令会报错，有别的方法实现吗
作者Author: sobereva 时间: 2022-5-27 03:14
是
没有简单的办法
作者Author: Melvin 时间: 2022-5-27 14:19
谢谢Sob老师
作者Author: 黑色桃花 时间: 2023-9-16 15:22
卢老师，感谢您之前的回帖，学生现在想对过渡态及初始及最终产物想进行单点能（算吉布斯自由能），但是涉及到离子（价态不为0，自旋多重度为1），那么这种情况怎么采用G4MP2-6X计算，或许说能用这种方法算吗，如果能，我该怎么改，请老师赐教，谢谢卢老师~~
作者Author: sobereva 时间: 2023-9-18 09:06
恰当设置净电荷，该怎么算还怎么算
作者Author: 黑色桃花 时间: 2023-9-18 21:09
好的，感谢卢老师
作者Author: 黑色桃花 时间: 2023-9-18 21:19
卢老师，学生再问一个“不恰当”的问题，这个方法在2023年算高级吗，或许说这个方法会不会被审稿人质疑方法落后，请您不吝赐教~谢谢卢老师
作者Author: sobereva 时间: 2023-9-18 23:50
不落后。除非将来理论方法有什么显著突破，否则在未来10年内都妥妥属于高精度范畴。
作者Author: 黑色桃花 时间: 2023-9-24 15:37
老师，我想请问一下，采用G4MP2-6X该怎么像shermo一样读取它的热校正值呢，有办法吗
作者Author: sobereva 时间: 2023-9-25 01:56
读懂了脚本自然清楚哪些变量是干什么的，输出出来就完了
作者Author: 黑色桃花 时间: 2023-9-25 13:48
哈哈，好的，谢谢卢老师
作者Author: 黑色桃花 时间: 2023-9-28 13:07
老师，如果算的这个体系是离子体系（0，1），上述的基组都包含弥散函数吗，可以用吗
作者Author: sobereva 时间: 2023-9-28 16:55
说清楚什么叫上述
作者Author: 黑色桃花 时间: 2023-9-29 16:31
好的老师，我就是看了G4MP2-6X热力学组合方法输入文件中含有很多的方法基组，那么这些基组对于需要考虑弥散函数的离子体系（电荷0，自旋多重度1），还适用嘛，换句话说就是输入文件的那些基组考虑弥散函数了嘛，谢谢老师的回复
作者Author: sobereva 时间: 2023-9-30 05:47
电荷为0明明不是离子体系
涉及到什么基组自行看方法原文便知
作者Author: sobereva 时间: 2023-9-30 14:47
先把方法原文看一遍，搞清楚原理，再把这perl脚本从头到尾读一遍，不懂perl语法就随手google查一下，撑死了2个小时也能搞明白脚本是怎么运作的了，自然也就知道怎么改。
作者Author: 林小浪 时间: 2023-9-30 14:48
sob老师您好，我想请教您一个问题。我最近在用G4MP2-XK方法计算CH4分子的H的时候，计算出的H数值总是不正确。
下图是我在G4MP2-XK方法原文J. Chem. Theory Comput. 2019, 15, 8, 4478–4484中下载的计算输入文件，根据自己的计算情况进行了修改，将原文自带的基组传到g09/basis目录下，计算了一个CH4分子。
(, 下载次数 Times of downloads: 19)
使用将该方法自带的perl文件，对计算生成的log文件进行计算，得到的H（298k）如下。得到H（298K）为0.04730458，这个数值显然是有问题的。
(, 下载次数 Times of downloads: 21)
之前使用G4MP2-6X计算的甲烷分子的H（298K）为 -40.42744909，这个数值是比较合理的。由于我使用的g09.d01版本在处理大一点的体系（比如TNT）时，会有l906报错，看了您的帖子，只有更换g09.e01以上版本能解决这个问题，所以想用G4MP2-XK代替G4MP2-6X。现在遇到这个问题，还希望sob老师能指点一二。
(, 下载次数 Times of downloads: 19)
作者Author: sobereva 时间: 2023-9-30 14:50
建议发邮件联系脚本的作者
作者Author: 林小浪 时间: 2023-9-30 14:57
好的，谢谢sob老师，节日快乐
作者Author: 黑色桃花 时间: 2023-9-30 16:27
好的，感谢老师的回复，老师您人真的太好了（我问了太多问题，您都不厌其烦地回复），真心谢谢您~
作者Author: tangdu 时间: 2024-11-12 21:55
您好，请问输出能量的指令是什么？我怎么也输出不出来
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |