计算化学公社

 找回密码 Forget password
 注册 Register
Views: 25820|回复 Reply: 24

[综合交流] 计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒人脚本

[复制链接 Copy URL]

4万

帖子

99

威望

4万

eV
积分
89863

管理员

公社社长+计算化学玩家

发表于 Post on 2021-1-28 20:41:00 | 显示全部楼层 Show all |阅读模式 Reading model
计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒人脚本

文/Sobereva@北京科音  2021-Jan-28


根据J. Phys. Chem. B., 121, 3864 (2017)中的测试,基于OPLS-AA力场的模拟很适合结合1.2*CM5原子电荷。虽然用1.2*CM5的时候水合自由能的计算精度稍逊于文中测试的另一种原子电荷1.14*CM1A-LBCC(因为这种原子电荷本来就是针对水合自由能训练的校正参数),但在计算其它属性上,如蒸发焓、密度,用1.2*CM5时的误差都更小,而且经验性更低、更有普适性。所以当大家用OPLS-AA力场研究新的小分子的时候,我比较推荐用1.2*CM5电荷。1.2*CM5电荷就是在Truhlar等人提出的原始的CM5电荷基础上乘上1.2,这相当于增大了原子电荷的数量级,等效地体现出了溶剂环境对溶质的极化作用。

笔者之前在《计算RESP原子电荷的超级懒人脚本(一行命令就算出结果)》(http://sobereva.com/476)和《RESP2原子电荷的思想以及在Multiwfn中的计算》(http://sobereva.com/531)中分别给出了超级懒人计算RESP和RESP2电荷的Linux脚本,完全不会用Gaussian的人都能轻松计算,已经有不少人在用。最近碰到完全不会用Gaussian的人试图计算1.2*CM5电荷向我求助,我遂又写了个类似的计算1.2*CM5电荷的懒人Linux脚本,在这里说一下。

这个脚本是Multiwfn文件包里的examples\scripts\1.2CM5.sh,在2021-Jan-28及之后更新的Multiwfn中才有,Multiwfn可以在其主页http://sobereva.com/multiwfn免费下载。

此脚本使用很简单:先确保Gaussian在当前机子里已经装了,见《Gaussian的安装方法及运行时的相关问题》(http://sobereva.com/439),也确保将Multiwfn按照手册2.1.2节的说明在机子里装了。之后假设1.2CM5.sh和一个结构文件phenol.xyz都在当前目录下,只需要运行./1.2CM5.sh phenol.xyz,之后在屏幕上就会看到脚本的运行过程:
Net charge was not defined. Default to 0
Spin multiplicity was not defined. Default to 1
Running optimization task via Gaussian...
Done!
Running formchk...
Running Multiwfn...
Finished! The optimized atomic coordinates with 1.2*CM5 charges (the last column) have been exported to phenol.chg in current folder


最终得到的phenol.chg用文本编辑器打开后可见,其中2、3、4列是优化后的XYZ坐标(埃),最后一列就是1.2*CM5电荷了。

脚本原理:这个脚本会将输入文件里的结构作为初猜结构用Gaussian在B3LYP-D3(BJ)/def2-SVP下做几何优化,然后调用formchk将得到的chk转化为fchk文件,然后调用Multiwfn计算CM5电荷,最后给出1.2*CM5电荷。

几点相关事项:
• 用的输入文件可以是任意含有结构信息而且Multiwfn支持的文件格式,比如pdb/mol/mol2/xyz/fch/gjf/wfn等等等等,见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》(http://sobereva.com/379)。
• 默认情况下电荷当成0,自旋多重度当成1。如果是比如不带电的二重态体系,应当写比如./1.2CM5.sh phenol.xyz 0 2。
• 如果运行脚本时提示没可执行权限,先运行chmod +x ./1.2CM5.sh
• 此脚本调用的是Gaussian 16,如果当前机子里装的是诸如Gaussian 09,需要将文件的第22行改为Gaussian=g09。
• 从脚本里的关键词可见,计算是在真空下进行的,这是刻意而为之,不要自行改成在溶剂模型下进行。
• 如果脚本运行失败,有这么几种可能:(1)Gaussian根本没恰当安装,应当确保能手动调用Gaussian运行一个最简单的计算任务 (2)Multiwfn没按照手册2.1.2节恰当安装。其中有些图形库不方便装的话,可以用Multiwfn的noGUI版,此时不需要装图形库 (3)Gaussian计算失败,可能是版本太老而不支持em=GD3BJ关键词,可将脚本中此关键词去掉。也可能几何优化或SCF不收敛,注意看计算中途产生的gau.out文件内容判断原因。
• 如果你手头已经有Gaussian或其它程序产生的Multiwfn支持的波函数文件了(fch/wfn/wfx/molden/mwfn等),就没必要再用这个脚本了。直接载入Multiwfn,依次输入7、16、1就能得到CM5电荷,再手动乘上1.2即可。
• 前述J. Phys. Chem. B.文章中只是考虑了中性体系,没有考虑带电体系。带电体系用什么原子电荷结合OPLS-AA没有统一说法,一般情况下不能用1.2*CM5电荷,因为此时总电荷都不是整数了。笔者建议此时用Multiwfn算RESP2电荷,虽然与OPLS-AA兼容性没有充分测试,但原理上问题不大。
• 如果使用此脚本计算1.2*CM5电荷,请按照Multiwfn程序包里How to cite Multiwfn.pdf文档的说明恰当引用Multiwfn。

如果你没买Gaussian,也可以用免费的ORCA量子化学程序结合Multiwfn算1.2*CM5原子电荷,笔者也提供了相应的傻瓜式脚本,见《ORCA结合Multiwfn计算RESP、RESP2和1.2*CM5原子电荷的懒人脚本》(http://sobereva.com/637)。

评分 Rate

参与人数
Participants 2
eV +10 收起 理由
Reason
zsu007 + 5
朙天儿 + 5

查看全部评分 View all ratings

北京科音自然科学研究中心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!

861

帖子

0

威望

2462

eV
积分
3323

Level 5 (御坂)

发表于 Post on 2021-1-28 21:39:20 | 显示全部楼层 Show all
本帖最后由 牧生 于 2021-1-29 07:38 编辑

日期有笔误,应该是
2021-Jan-28及之后更新的Multiwfn中才有。

另有一个疑问,“带电体系用什么原子电荷结合OPLS-AA没有统一说法,一般情况下不能用1.2*CM5电荷,因为此时总电荷都不是整数了。笔者建议此时用Multiwfn算RESP2电荷,虽然与OPLS-AA兼容性没有充分测试,但原理上问题不大。”

我的理解,如果做表面活性剂自组装,比如十二烷基磺酸钠或别的表面活性剂,想用opls力场,此时应该是1.2*CM5电荷也可以(十二烷基磺酸根就用1.2*CM5脚本计算,得到总电荷就是-1.2了,钠离子就用ions自带的,但手动改为电荷+1.2,行吗??),使用RESP2电荷也可以的??

如果做不带电的体系油水分离,如甲苯和水,想要用opls力场,1.2*CM5和RESP2哪个会更好呢?

如果我要做聚合物的模拟,如聚丙烯酸钠,要用opls力场,那么,应该用1.2*CM5电荷计算三聚体,然后手动重复建立聚合物链,此时还能不能用RESP2电荷?





4万

帖子

99

威望

4万

eV
积分
89863

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2021-1-29 19:34:37 | 显示全部楼层 Show all
牧生 发表于 2021-1-28 21:39
日期有笔误,应该是
2021-Jan-28及之后更新的Multiwfn中才有。

笔误改了

Na用+1.2明显说不通。用RESP2明显有意义得多

用1.2*CM5,毕竟有JPCB的文献支撑。水的原子电荷不要自己算,直接用水模型里定义的。

用1.2*CM5
北京科音自然科学研究中心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!

108

帖子

0

威望

2713

eV
积分
2821

Level 5 (御坂)

发表于 Post on 2021-2-20 17:19:45 | 显示全部楼层 Show all
牧生 发表于 2021-1-28 21:39
日期有笔误,应该是
2021-Jan-28及之后更新的Multiwfn中才有。

我觉得做聚合物用原力场电荷参数就挺好的。

312

帖子

1

威望

3106

eV
积分
3438

Level 5 (御坂)

发表于 Post on 2021-3-2 08:23:14 | 显示全部楼层 Show all
请问楼主和社长研究小分子多聚现象的时候,GAFF小分子+amber99sb-ildn水好?还是OPLS-AA好?

116

帖子

0

威望

1387

eV
积分
1503

Level 5 (御坂)

发表于 Post on 2021-3-2 19:49:47 | 显示全部楼层 Show all
exity 发表于 2021-3-2 08:23
请问楼主和社长研究小分子多聚现象的时候,GAFF小分子+amber99sb-ildn水好?还是OPLS-AA好?

GAFF2+amber14 SB ,RESP2电荷

312

帖子

1

威望

3106

eV
积分
3438

Level 5 (御坂)

发表于 Post on 2021-3-3 08:30:48 | 显示全部楼层 Show all
nianbin 发表于 2021-3-2 19:49
GAFF2+amber14 SB ,RESP2电荷

谢谢

861

帖子

0

威望

2462

eV
积分
3323

Level 5 (御坂)

发表于 Post on 2021-3-7 10:23:41 | 显示全部楼层 Show all
本帖最后由 牧生 于 2021-3-7 10:49 编辑

https://doi.org/10.1021/acs.jctc.9b00947

针对带电的SDS,RESP电荷比CM5好一些(文中没有使用1.2*CM5)

4万

帖子

99

威望

4万

eV
积分
89863

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2021-3-8 00:39:04 | 显示全部楼层 Show all
牧生 发表于 2021-3-7 10:23
https://doi.org/10.1021/acs.jctc.9b00947

针对带电的SDS,RESP电荷比CM5好一些(文中没有使用1.2*CM5 ...

那篇文章对拟合静电势电荷的计算方式不好。MP2直接在气相下算没有考虑溶剂的等效极化效应,而在SMD模型下算时,又没有像RESP2(0.5)那样避免高估极化程度。而且体系是阴离子,也没带弥散函数。鉴于原子电荷对结果影响极大,应当测试得更严谨才对。
北京科音自然科学研究中心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!

861

帖子

0

威望

2462

eV
积分
3323

Level 5 (御坂)

发表于 Post on 2021-3-8 07:49:27 | 显示全部楼层 Show all
本帖最后由 牧生 于 2021-3-8 07:51 编辑
sobereva 发表于 2021-3-8 00:39
那篇文章对拟合静电势电荷的计算方式不好。MP2直接在气相下算没有考虑溶剂的等效极化效应,而在SMD模型下 ...

所以,考虑使用RESP2才是更好的选择吧

4万

帖子

99

威望

4万

eV
积分
89863

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2021-3-9 04:30:10 | 显示全部楼层 Show all
牧生 发表于 2021-3-8 07:49
所以,考虑使用RESP2才是更好的选择吧

原理上是
北京科音自然科学研究中心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!

66

帖子

0

威望

381

eV
积分
447

Level 3 能力者

发表于 Post on 2021-7-22 22:22:15 | 显示全部楼层 Show all
本帖最后由 Sasuke 于 2021-7-22 22:23 编辑

社长 我用了这个脚本 高斯版本是gaussian09 Multiwfn是目前最新版3.8 然后高斯软件是可以运行的因为成功试算了里面自带的算例.Multiwfn也是按照手册2.1.2装的。然后我用MS画了苯乙烯三聚体导出mol2文件作为这个脚本的输入文件,结果运行报错如图所示。gaussian09的out文件里面显示输入文件有误 A syntax error was detected in the input line. 是高斯版本的问题吗?

高斯输出文件

高斯输出文件

脚本运行后报错

脚本运行后报错

4万

帖子

99

威望

4万

eV
积分
89863

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2021-7-23 03:41:33 | 显示全部楼层 Show all
Sasuke 发表于 2021-7-22 22:22
社长 我用了这个脚本 高斯版本是gaussian09 Multiwfn是目前最新版3.8 然后高斯软件是可以运行的因为成功试 ...

Gaussian版本太老,不支持DFT-D3。去掉em=GD3BJ
北京科音自然科学研究中心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!

71

帖子

0

威望

315

eV
积分
386

Level 3 能力者

发表于 Post on 2021-8-13 11:11:59 | 显示全部楼层 Show all
请问下社长,我的分子是个聚合物,原子数目有点多1050个,这个用56个核的双路服务器Gaussian还能算得动吗?
我提交了两次似乎都在运行一段时间后报错,错误信息是:
  1. Error in internal coordinate system.
  2. Error termination via Lnk1e in ***/g16/l103.exe at Thu Aug 12 10:50:23 2021.
  3. Job cpu time:       5 days 16 hours 46 minutes 17.4 seconds.
  4. Elapsed time:       0 days  2 hours 27 minutes 48.4 seconds.
  5. File lengths (MBytes):  RWF=  56574 Int=      0 D2E=      0 Chk=   2389 Scr=      2
复制代码

请问下这个错误我应该如何解决呢?非常感谢!

4万

帖子

99

威望

4万

eV
积分
89863

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2021-8-13 22:48:45 | 显示全部楼层 Show all
mgqqlwq 发表于 2021-8-13 11:11
请问下社长,我的分子是个聚合物,原子数目有点多1050个,这个用56个核的双路服务器Gaussian还能算得动吗? ...

算不动
另外,这种报错可以在脚本里的opt关键词里加入cartesian关键词解决。

聚合物通常都是对重复单元来计算原子电荷,然后拼接成整体,而非直接对整个聚合物直接算。
北京科音自然科学研究中心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-1 23:00 , Processed in 0.469506 second(s), 25 queries .

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