请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 1175|回复: 9

[综合交流] 计算RESP原子电荷的超级懒人脚本(一行命令就算出结果)

[复制链接]

1万

帖子

25

威望

2万

eV
积分
44694

管理员

公社社长

发表于 2019-4-17 20:51:35 | 显示全部楼层 |阅读模式
计算RESP原子电荷的超级懒人脚本(一行命令就算出结果)

文/Sobereva@北京科音
First release: 2019-Apr-17  Last update: 2019-May-8

RESP电荷主要优点是特别适合对柔性分子做分子动力学计算。之前笔者写过《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》(http://sobereva.com/441)一文,极其详细地介绍了RESP电荷,并结合实际例子充分展现了Multiwfn在计算RESP电荷上的强大性,可谓极度便利而且极度灵活,明显比Antechamber等其它能计算RESP的程序好用、方便得多得多得多。但在网上的答疑过程中,笔者也看到很多要计算RESP电荷的人完全是量子化学外行,甚至GaussView都不会用,Gaussian的关键词一点也不会写,需要一个极度傻瓜化的“一键”工具来计算RESP电荷,为此笔者在这里介绍一个基于Gaussian和Multiwfn的仅仅需要写一条命令就可以算出RESP电荷的Linux环境下的脚本,可谓将RESP电荷的计算极尽简化,不可能做到更简单了。

1 准备工作

在机子里安装Gaussian。不会装的话仔细看《Gaussian的安装方法及运行时的相关问题》(http://sobereva.com/439)。

http://sobereva.com/multiwfn下载Multiwfn,必须是2019-Apr-17及以后的版本。然后严格按照手册2.1.2节说的方式安装到Linux下面。

对Multiwfn目录下的examples\RESP\目录下的RESP.sh根据需要进行恰当的修改。此脚本默认调用的是g09,Gaussian16的用户应当将其中的g09替换为g16。此脚本默认为在B3LYP-D3(BJ)/def2-SVP级别下对结构进行优化,然后在B3LYP-D3(BJ)/def2-TZVP级别下计算单点能同时产生分子表面静电势数据,然后调用formchk把chk文件转化为fchk,最后调用Multiwfn产生RESP电荷。Gaussian计算过程中使用了隐式溶剂模型表现水环境。计算级别、溶剂设定信息在这个脚本里都体现得很清楚,大家可以根据需要对关键词进行修改。此脚本默认的计算级别算出的RESP电荷已经十分准确。几何优化部分如果算不动或者想显著降低耗时的话,可以把原有的关键词改为# PM6D3 opt,即使用PM6-D3半经验方法在真空中做优化,但这个级别只适合普通有机体系而且没有体系局部显离子特征的情况。单点任务部分算不动的话,可以把def2-TZVP降到6-311G**,如果是阴离子建议改为6-311+G**。有时候对柔性大体系的优化过程可能难收敛,可以尝试把opt改为诸如opt(loose,gdiis)之类使得达到收敛更容易的关键词,详见《量子化学计算中帮助几何优化收敛的常用方法》(http://sobereva.com/164)。


2 使用脚本

将Multiwfn目录下的examples\RESP\目录下的RESP.sh拷到当前目录,使用chmod +x ./RESP.sh给它增加可执行权限。把含有分子结构信息的文件也拷到当前目录,格式必须是Multiwfn认识的(可以是xyz/mol/mol2/pdb/fch/wfn/molden/gjf等等)。然后运行下面的命令:
./RESP.sh [文件名] [净电荷] [自旋多重度]
比如./RESP.sh H2O.xyz 0 1。如果净电荷和自旋多重度不写,分别默认为0和1。

计算完毕后,当前目录下就有了和输入文件同名的.chg文件。这是个文本文件,前四列是元素名和优化后的坐标,最后一列就是RESP电荷。以下是运算过程的输出例子:

[root@192 other2]# ./RESP.sh H2O.xyz
Net charge was not defined. Default to 0
Spin multiplicity was not defined. Default to 1
Running optimization task via Gaussian...
Running single point task via Gaussian...
Running formchk...
Running Multiwfn...
Finished! The optimized atomic coordinates with RESP charges (the last column) have been exported to H2O.chg in current folder


本文的做法只是以最一般的方式计算RESP电荷,如果需要更灵活的拟合方式,需要在Multiwfn的RESP模块的界面里进行设置,在《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》(http://sobereva.com/441)文中有充分、详细的介绍。

如果以本文的方式计算了RESP电荷,请别忘了引用Multiwfn原文。引用方式在《Multiwfn FAQ》(http://sobereva.com/452》里写明了。

评分

参与人数 7eV +33 收起 理由
ezez + 5 赞!
molu851 + 3 精品内容,十分有用,之前需要等待高斯优化.
azero + 5 赞!
yjmaxpayne + 5 好物!
fallleave + 5 好物!
k64_cc + 5 GJ!
zikuzi + 5 赞!

查看全部评分

北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

37

帖子

1

威望

1168

eV
积分
1225

Level 4 (黑子)

发表于 2019-4-18 15:17:55 | 显示全部楼层
报这个错误,未找到原因,求助
Net charge was not defined. Default to 0
Spin multiplicity was not defined. Default to 1
forrtl: severe (59): list-directed I/O syntax error, unit -5, file Internal List-Directed Read
Image              PC                Routine            Line        Source            
Multiwfn           000000000198D709  Unknown               Unknown  Unknown
Multiwfn           00000000019B7F3A  Unknown               Unknown  Unknown
Multiwfn           00000000019B64BB  Unknown               Unknown  Unknown
Multiwfn           00000000005E85B9  Unknown               Unknown  Unknown
Multiwfn           00000000005A8AF1  Unknown               Unknown  Unknown
Multiwfn           0000000000657E16  Unknown               Unknown  Unknown
Multiwfn           000000000040B41E  Unknown               Unknown  Unknown
libc-2.17.so       00002B7926F143D5  __libc_start_main     Unknown  Unknown
Multiwfn           000000000040B329  Unknown               Unknown  Unknown

1万

帖子

25

威望

2万

eV
积分
44694

管理员

公社社长

 楼主| 发表于 2019-4-18 17:36:24 | 显示全部楼层
wangzhehyd 发表于 2019-4-18 15:17
报这个错误,未找到原因,求助
Net charge was not defined. Default to 0
Spin multiplicity was n ...

可以修改下脚本,检查一下中间产生的Gaussian输入输出文件。注意一定得是最新的Multiwfn版本
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

37

帖子

1

威望

1168

eV
积分
1225

Level 4 (黑子)

发表于 2019-4-18 20:24:13 | 显示全部楼层
sobereva 发表于 2019-4-18 17:36
可以修改下脚本,检查一下中间产生的Gaussian输入输出文件。注意一定得是最新的Multiwfn版本

谢谢,已解决。请问Gaussian部分的计算是否可以通过TeraChem完成,然后再用multiwfn处理得到resp电荷?

1万

帖子

25

威望

2万

eV
积分
44694

管理员

公社社长

 楼主| 发表于 2019-4-19 03:35:36 | 显示全部楼层
wangzhehyd 发表于 2019-4-18 20:24
谢谢,已解决。请问Gaussian部分的计算是否可以通过TeraChem完成,然后再用multiwfn处理得到resp电荷?

我没有测试过Terachem产生的molden文件是否规矩,能否正确被Multiwfn读取(如果有问题通常载入文件的时候程序就会有警告)
如果能正确读取,把脚本稍微改写一下就可以实现Multiwfn+Terachem产生RESP电荷。如果不能正常读取,还得多加一步调用molden2aim把molden文件标准化的过程
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

42

帖子

0

威望

1678

eV
积分
1720

Level 5 (御坂)

发表于 2019-5-8 18:05:37 | 显示全部楼层
厉害,这个超方便,而且得到了结构优化的小分子,可以直接对接了。

另外请教一下sob老师,画出来的小分子用AM1半经验优化对分子对接结果应该无大影响吧?

1万

帖子

25

威望

2万

eV
积分
44694

管理员

公社社长

 楼主| 发表于 2019-5-8 18:28:30 | 显示全部楼层
azero 发表于 2019-5-8 18:05
厉害,这个超方便,而且得到了结构优化的小分子,可以直接对接了。

另外请教一下sob老师,画出来的小分 ...

虽然AM1也不是不能用,但已经过时了。如果为了优化时候图快,建议你改成Gaussian能支持的更好的PM6D3或PM7
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

42

帖子

0

威望

1678

eV
积分
1720

Level 5 (御坂)

发表于 2019-5-8 18:39:10 | 显示全部楼层
sobereva 发表于 2019-5-8 18:28
虽然AM1也不是不能用,但已经过时了。如果为了优化时候图快,建议你改成Gaussian能支持的更好的PM6D3或PM ...

这个样子,毕竟不是学量化的。
还是用sob老师你脚本优化出来的坐标省事。才几个小分子,花不了多少时间

62

帖子

1

威望

209

eV
积分
291

Level 3 能力者

发表于 2019-6-19 18:38:26 | 显示全部楼层
本帖最后由 wuzhiyi 于 2019-6-19 18:44 编辑

sob老师好,我是做生物出生的对量化不是很了解。传统上,我们一般用amber14SB来表示蛋白,然后用gaff来表示小分子的LJ,bonded/angle/dihedral。对于电荷的处理,人们一般要么用AM1-BCC,要么用HF/6-31G*做RESP,人们给出的理由是,因为LJ,bonded/angle/dihedral都是以partial charge是按HF/6-31G*来计算为基础来定义的所有不可以换。


想问一下sob老师对这种观点的看法如何?我稍微找了找,发现好像大部分时候,人们从新拟合电荷的时候都是在HF/6-31G*下拟合的。比如:
Development and benchmark to obtain AMBER parameters dataset for non-standard amino acids modified with 4-hydroxy-2-nonenal
Forcefield_PTM: Ab Initio Charge and AMBER Forcefield Parameters for Frequently Occurring Post-Translational Modifications
Reparametrization of Protein Force Field Nonbonded Interactions Guided by Osmotic Coefficient Measurements from Molecular Dynamics Simulations
Revised AMBER parameters for bioorganic phosphates


而且最新的FF19SB (https://chemrxiv.org/articles/ff ... in_Solution/8279681)仍然是沿用HF/6-31G*的电荷,然后在这个基础上重新拟合dihedral。
虽然ff15ipq (https://pubs.acs.org/doi/full/10.1021/acs.jctc.6b00567)是用IPolQ在MP2/cc-pVTZ下拟合的电荷,但好像没什么人用。想问一下sob老师能不能推荐一下用amberFF配合B3LYP-D3(BJ)/def2-TZVP来计算电荷的文献?

想问一下,用更好的级别比如B3LYP-D3(BJ)/def2-TZVP主要在哪一方面改善RESP的表现?我们一般用溶解自由能来作为一个很粗略的检测手段,想问一下,有没有这样的研究比如找一个小分子的溶解数据库(类似Mobley的FreeSolv:https://github.com/MobleyLab/FreeSolv,这其中有AM1-BCC和试验结果的对比)然后对比B3LYP-D3(BJ)/def2-TZVP和HF/6-31G*以及实验结果进行比较之类的。

另想问一下IPolQ这种改善隐式溶剂对水环境的表现的方法,sob老师的看法如何?

谢谢

1万

帖子

25

威望

2万

eV
积分
44694

管理员

公社社长

 楼主| 发表于 2019-6-19 21:28:30 | 显示全部楼层
wuzhiyi 发表于 2019-6-19 18:38
sob老师好,我是做生物出生的对量化不是很了解。传统上,我们一般用amber14SB来表示蛋白,然后用gaff来表示 ...

早在2003年提出AMBER03的时候原作者就已经把计算RESP电荷的水平提升到了B3LYP/cc-pVTZ结合IEFPCM隐式溶剂模型,还死咬早就过时的HF/6-31G*是完全匪夷所思的

稍有量化常识就知道,B3LYP-D3(BJ)/def2-TZVP结合隐式溶剂模型对波函数的描述精度比HF/6-31G*强百倍,后者那是当年计算条件差,不得已而为之,而且只能近似利用HF波函数质量特别恶心、高估极性这点来非常近似地体现水溶液中被极化的情况。显然如今这年头没有任何理由去用HF/6-31G*。

没有恰好现成的B3LYP-D3(BJ)/def2-TZVP结合隐式溶剂模型算RESP电荷的文献,因为这个级别是我的主张,而且显然是非常合理的。def2-TZVP的档次和cc-pVTZ接近,而对于DFT更划算。B3LYP加了D3去优化结构明显比起用B3LYP对于弱相互作用对体系结构影响显著的时候强得多。见
谈谈“计算时是否需要加DFT-D3色散校正?”
http://sobereva.com/413
B3LYP-D3(BJ)/def2-TZVP结合隐式溶剂模型算RESP电荷你大可放心去用,除非碰上缺乏量化常识的审稿人,否则绝对不会被怼,即便真碰上了,也有充分理由去回击。

计算级别问题我已经在此文1.5节里有非常详细的说明
RESP拟合静电势电荷的原理以及在Multiwfn中的计算
http://sobereva.com/441http://bbs.keinsci.com/thread-10880-1-1.html

IPolQ没怎么去特意关注,没什么意思。

FF19SB沿用HF/6-31G*是考虑和以往老版本蛋白质骨架参数的兼容性,然而自己算RESP电荷都是算的新的有机小分子,完全不需要考虑这种兼容性问题。
如果AMBER力场的开发者现在从头重新开发一套新的力场,我相信100%不会再去用HF/6-31G*算RESP电荷。

北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949-1号 )

GMT+8, 2019-8-21 09:44 , Processed in 0.169222 second(s), 25 queries .

快速回复 返回顶部 返回列表