计算化学公社

 找回密码 Forget password
 注册 Register
Views: 4914|回复 Reply: 7

[综合交流] ORCA结合Multiwfn计算RESP、RESP2和1.2*CM5原子电荷的懒人脚本

[复制链接 Copy URL]

4万

帖子

99

威望

4万

eV
积分
89946

管理员

公社社长+计算化学玩家

发表于 Post on 2022-3-8 06:09:54 | 显示全部楼层 Show all |阅读模式 Reading model
ORCA结合Multiwfn计算RESP、RESP2和1.2*CM5原子电荷的懒人脚本

文/Sobereva@北京科音
First release: 2022-Mar-15  Last update: 2022-Aug-6


之前笔者在以下文章中提供了三个Linux shell脚本,分别用来自动调用机子里的Gaussian和Multiwfn程序实现一键计算1.2*CM5、RESP和RESP2原子电荷,它们对于做经典力场的分子动力学非常重要。
计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒人脚本
http://sobereva.com/585http://bbs.keinsci.com/thread-21462-1-1.html
计算RESP原子电荷的超级懒人脚本(一行命令就算出结果)
http://sobereva.com/476http://bbs.keinsci.com/thread-12858-1-1.html
RESP2原子电荷的思想以及在Multiwfn中的计算
http://sobereva.com/531http://bbs.keinsci.com/thread-16190-1-1.html

如今用免费的ORCA量子化学程序的人也很多。为了便于那些主要做分子动力学模拟、没买Gaussian又不太懂ORCA程序使用的人也能便利地计算上述原子电荷,笔者写了能自动调用ORCA和Multiwfn的计算那些电荷的Linux shell脚本。这些脚本提供在了Multiwfn程序文件包里(可以在Multiwfn主页http://sobereva.com/multiwfn免费下载),必须是2022年3月8日及以后更新的Multiwfn版本里才有,包括:
examples\RESP\RESP_ORCA.sh:计算RESP电荷的脚本
examples\RESP\RESP2_ORCA.sh:计算RESP2电荷的脚本
examples\scripts\1.2CM5_ORCA.sh:计算1.2*CM5电荷的脚本

这些脚本的用法和前述帖子里介绍的基于Gaussian的脚本精确一致,需要留意的地方也都相同,只不过脚本中调用Gaussian的地方变成了调用ORCA而已,故不再累述用法。

这些脚本运行之前记得用文本编辑器打开,把ORCA=和orca_2mkl=后面的内容分别改为当前机子里实际的ORCA和orca_2mkl工具的路径。并把nprocs=后面的值改为计算时要调用的CPU核心数。如果不会装ORCA的话看《量子化学程序ORCA的安装方法》(http://sobereva.com/451)。

RESP_ORCA.sh和RESP2_ORCA.sh里默认用的优化级别和基于Gaussian的脚本(RESP.sh和RESP2.sh)有所不同,这里用的是ORCA才支持的B97-3c,因为这个级别做优化很快,结果准确度也不错。由于这个差异,以及ORCA和Gaussian在溶剂模型的实现上有所差异,所以基于Gaussian和基于ORCA的脚本得到的RESP或RESP2电荷可能有零点零几的差别,这点没必要在意,都是合理的。

关于使用脚本时哪些溶剂可以直接用、溶剂名怎么写,请在ORCA手册里搜“solvents in the SMD library”查看内置的溶剂名列表。注意ORCA里有些溶剂名是带空格的,对这种情况要把溶剂名用双引号扩住,例如./RESP2_ORCA.sh HF.pdb 0 1 "ETHYL ETHANOATE"。

如果你的输入的结构文件里的结构就已经足够好,不想让脚本自动再做优化浪费时间,可以用examples\RESP\目录下的RESP_ORCA_noopt.sh和RESP2_ORCA_noopt.sh分别代替前述的RESP_ORCA.sh和RESP2_ORCA.sh,它们的用法完全一样,只不过带_noopt后缀的不做优化步骤。

使用这些脚本计算原子电荷发表文章的话请在文中恰当引用Multiwfn和ORCA。

评分 Rate

参与人数
Participants 7
eV +30 收起 理由
Reason
ezez + 5 赞!
Frozen-Penguin + 4 GJ!
ggdh + 5 谢谢
kay + 5 赞!
tjuptz + 3 懒人福音
zsu007 + 5 赞!
牧生 + 3 期待已久

查看全部评分 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!

52

帖子

2

威望

2374

eV
积分
2466

Level 5 (御坂)

发表于 Post on 2022-3-13 20:30:19 | 显示全部楼层 Show all
貌似比用gaussian快挺多

52

帖子

2

威望

2374

eV
积分
2466

Level 5 (御坂)

发表于 Post on 2022-3-14 10:39:13 | 显示全部楼层 Show all
以下情况在脚本里面不兼容一下么

Atomic radii used:
Element:C      vdW radius (Angstrom): 1.500
Element:N      vdW radius (Angstrom): 1.500
Element:S      vdW radius (Angstrom): 1.750
vdW radius used in fitting for Br is missing, input the radius (Angstrom)
Hint: If pressing ENTER button directly, corresponding UFF radii scaled by 1/1.2 will be used, this is a usually reasonable choice
forrtl: severe (59): list-directed I/O syntax error, unit -5, file Internal List-Directed Read
Image              PC                Routine            Line        Source            
Multiwfn           00000000021432BB  Unknown               Unknown  Unknown
Multiwfn           00000000021705A2  Unknown               Unknown  Unknown
Multiwfn           000000000216F328  Unknown               Unknown  Unknown
Multiwfn           000000000089E25E  Unknown               Unknown  Unknown
Multiwfn           00000000008A53EF  Unknown               Unknown  Unknown
Multiwfn           00000000008799D6  Unknown               Unknown  Unknown
Multiwfn           000000000082724A  Unknown               Unknown  Unknown
Multiwfn           000000000040AA22  Unknown               Unknown  Unknown
libc-2.28.so       000014AE61C3B493  __libc_start_main     Unknown  Unknown
Multiwfn           000000000040A92E  Unknown               Unknown  Unknown

4万

帖子

99

威望

4万

eV
积分
89946

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2022-3-15 06:02:22 | 显示全部楼层 Show all
wangzhehyd 发表于 2022-3-14 10:39
以下情况在脚本里面不兼容一下么

Atomic radii used:

我今日更新了Multiwfn程序和里面的脚本,最新版在遇到没有半径的时候会自动用UFF radii scaled by 1/1.2。

脚本里在Multiwfn运行时加上了-ispecial 1选项,这就是要求遇到没半径的时候不问用户而自动用推荐值。
北京科音自然科学研究中心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!

52

帖子

2

威望

2374

eV
积分
2466

Level 5 (御坂)

发表于 Post on 2022-3-15 09:47:37 | 显示全部楼层 Show all
sobereva 发表于 2022-3-15 06:02
我今日更新了Multiwfn程序和里面的脚本,最新版在遇到没有半径的时候会自动用UFF radii scaled by 1/1.2 ...

昨天我还改源代码去解决这个问题,今天就有官方支持版本了,高效

此外,对于”对于第四周期之后是赝势基组。由于molden格式有不能记录实际核电荷数的缺陷,因此若你的体系包含第四周期之后的原子,最后给出的原子电荷将是错误的。对这样的体系请手动用ORCA计算得到molden文件并进行修改,然后自行用Multiwfn进行计算”这个问题,比如碘原子,我尝试在脚本中增加了一句(最后一句)
  1. #Convert to .molden file
  2. echo Running orca_2mkl...
  3. $orca_2mkl SP -molden > /dev/null
  4. mv SP.molden.input SP.molden
  5. sed -i '/^I/s/ 53 / 25 /' SP.molden
复制代码


测试了几个分子都可以通过,但是不清楚实际核电荷数会不会随分子不同发生变化。
如果会发生变化,有没有办法在Multiwfn处理molden文件前先获得这个实际的核电荷数,加一步自动处理molden文件

期待官方解决方案。

4万

帖子

99

威望

4万

eV
积分
89946

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2022-3-15 11:18:17 | 显示全部楼层 Show all
wangzhehyd 发表于 2022-3-15 09:47
昨天我还改源代码去解决这个问题,今天就有官方支持版本了,高效

此外,对于”对于第四周期 ...

用def2的时候,核电荷数只取决于元素,与分子无关
在脚本里自动对molden做替换确实是个不错的主意
北京科音自然科学研究中心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!

4万

帖子

99

威望

4万

eV
积分
89946

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2022-3-15 12:41:21 | 显示全部楼层 Show all
我刚刚更新了官网上Multiwfn最新版里的本帖提及的脚本,对于def2支持的所有元素(前六周期)都能用了,不再限于前四周期
北京科音自然科学研究中心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!

52

帖子

2

威望

2374

eV
积分
2466

Level 5 (御坂)

发表于 Post on 2022-4-7 13:44:45 | 显示全部楼层 Show all
本帖最后由 wangzhehyd 于 2022-4-7 13:46 编辑

基于1949个分子(随机挑选),比较用gaussian和orca算出来的resp原子电荷的差异。
gas情况下
by_element.png

solv(water)情况下
by_element.png


本版积分规则 Credits rule

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

GMT+8, 2023-2-6 05:40 , Processed in 0.429891 second(s), 26 queries .

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