计算化学公社

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

[其它] 关于为什么Multiwfn算的出RESP电荷与Antechamber的有所差异

[复制链接 Copy URL]

4万

帖子

99

威望

4万

eV
积分
89866

管理员

公社社长+计算化学玩家

发表于 Post on 2019-9-27 02:32:32 | 显示全部楼层 Show all |阅读模式 Reading model
关于为什么Multiwfn算的出RESP电荷与Antechamber有所差异

文/Sobereva@北京科音 2019-Sep-27

0 前言

RESP电荷的原理、细节以及计算在《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》(http://sobereva.com/441)有详细说明,Multiwfn是计算RESP电荷最方便易用和灵活的程序。有Multiwfn用户发现,Multiwfn算出的某些分子的RESP电荷与antechamber的不同,在此文我说一下原因。使用Multiwfn算RESP电荷的用户不是很有必要看此文,Multiwfn给的结果肯定是合理的。但如果想了解一些细节的话也可以看看。

简单来说导致差异的原因主要有两点:
(1)拟合点的分布不同
(2)等价约束设置不同

下文涉及的体系的相关文件可以在这里下:http://sobereva.com/attach/516/file.rar。本文用的是2019-Sep-26更新的Multiwfn 3.7(dev)版。


1 拟合点的分布不同带来的差异

Multiwfn在计算RESP电荷的时候,对于3.6、3.7版默认用的是MK方式设置拟合点(以后版本有可能会改变,届时看Multiwfn的RESP模块界面上的选项提示便知默认的是什么),拟合点位置由Multiwfn自己的代码设置,静电势由Multiwfn自己的代码计算或者调用cubegen计算。Antechamber计算RESP电荷时是从Gaussian的pop=MK IOp(6/33=2) IOp(6/42=6)输出文件里直接读取MK拟合点的坐标和静电势。虽然MK原文虽然定义了拟合点的分布规则,即每个原子上在不同原子半径倍数上有几层拟合点,但在具体确定拟合点坐标上,不同程序的细节不同,比如在球面上的拟合点以什么算法生成排布?拟合点的密度是多少?在原子接缝的地方是否自动去除过密的拟合点?等等。由于Multiwfn设置的拟合点的位置和Gaussian在这方面有一些不同,导致Multiwfn和Antechamber算的RESP电荷会有些许差别,但这个因素造成的差异通常不大,也没法说哪个程序更对,因为都是合理的。注:前述IOp(6/33=2)是让Gaussian把拟合点的位置和静电势数值打印到输出文件里的选项,而IOp(6/42=6)用来将拟合点的密度设为高于默认值。在Multiwfn的RESP电荷计算界面里也可以直接设拟合点的密度。

下面我们拿乙醇这个体系来对比一下。这里笔者直接用gview画了一个乙醇,保存成了.mol2文件。然后用antechamber产生了对应的gjf文件,里面的关键词是#HF/6-31G* SCF=tight Test Pop=MK iop(6/33=2) iop(6/42=6),原本带着的opt我给删了,对于本文的讨论没意义。然后用Gaussian 09 E.01进行计算,得到的chk文件我转成了fch文件。

用Multiwfn载入fch文件,进入RESP计算界面,直接选标准两步式RESP电荷计算,结果如下
Center      Charge
  1(C )   -0.250749
  2(H )    0.064449
  3(H )    0.064449
  4(H )    0.064449
  5(C )    0.547451
  6(H )   -0.091298
  7(H )   -0.091298
  8(O )   -0.714589
  9(H )    0.407135

在计算第二步的时候提示用了以下等价性约束
Constraint   1:    2(H )    3(H )    4(H )
Constraint   2:    6(H )    7(H )

再用antechamber计算RESP电荷,会出现一批文件,其中有两个RESP计算模块输出的文件ANTECHAMBER_RESP1.OUT和ANTECHAMBER_RESP2.OUT。第一个就是第一步拟合的输出,第二个就是第二步拟合的输出,最终结果看第二个里面的下面的部分

         Point Charges Before & After Optimization

    no.  At.no.    q(init)       q(opt)     ivary    d(rstr)/dq
    1   6      -0.229693      -0.229799      0       0.003990
    2   1       0.077688       0.059374      0       0.000000
    3   1       0.077682       0.059374      2       0.000000
    4   1       0.034340       0.059374      2       0.000000
    5   6       0.464072       0.533889      0       0.001841
    6   1      -0.057092      -0.086155      0       0.000000
    7   1      -0.057098      -0.086155      6       0.000000
    8   8      -0.720426      -0.720426    -99       0.001375
    9   1       0.410526       0.410526    -99       0.000000

q(init)是第一步拟合之后产生的电荷,q(opt)是经过第二步拟合最终得到的RESP电荷。ivary是等价性约束关系,比如no.为3和4的原子的ivary是2,就代表2、3、4在拟合过程中保持等价,这和前面看到的Multiwfn自动用的约束完全一样。8和9号原子的ivary为-99,代表它们在第二步拟合过程中电荷被约束为保持不变。

Antechamber这里产生的RESP电荷和前面给的Multiwfn算出来的有轻微差异,但顶多也就差0.01几。如果我们想让Multiwfn算出的电荷与Antechamber的精确一样,那就在选择开始计算RESP电荷之前先选一下选项8,然后再选1开始标准两步式拟合,按屏幕提示的要求把Gaussian输出文件路径输进去,这时Multiwfn就和Antechamber一样从Gaussian输出文件里读取拟合点的位置和静电势了,输出信息如下,和Antechamber给出的完全相同(个别原子仍有10的负六次方的差异是收敛限设置、舍入误差等数值细节导致的,完全可忽略不计)
Center      Charge
    1(C )   -0.229798
    2(H )    0.059373
    3(H )    0.059373
    4(H )    0.059373
    5(C )    0.533887
    6(H )   -0.086155
    7(H )   -0.086155
    8(O )   -0.720426
    9(H )    0.410526


2 等价约束设置不同带的差异

下面我们看苯甲醛这个例子。哪怕Multiwfn也直接从Gaussian输出文件里读取拟合点的信息,结果与Antechamber给出的仍有轻微差异,这是由于两个程序用的等价性约束设置不同所带来的。

我们先用Multiwfn基于Gaussian里的拟合点信息计算RESP电荷。启动Multiwfn,载入本文文件包里的C6H5CHO.fchk,进入RESP模块,先选择8,再选择1,再输入C6H5CHO.out的路径,得到的结果如下
   Center      Charge
     1(C )   -0.115868
     2(C )   -0.111841
     3(C )   -0.114944
     4(C )   -0.146385
     5(H )    0.137991
     6(H )    0.134593
     7(H )    0.153718
     8(H )    0.139930
     9(C )    0.391423
    10(H )    0.059972
    11(C )   -0.223448
    12(H )    0.155781
    13(C )    0.034578
    14(O )   -0.495500
Sum of charges:    0.000000
RMSE:    0.001373   RRMSE:    0.079595

从屏幕上的提示可见,这个体系的RESP电荷计算过程只做了第一步,因为按照RESP电荷定义的规则,这个体系本身不需要做第二步。在第一步计算的时候没有自动做任何等价性约束。

也用antechamber进行计算,从给出的ANTECHAMBER_RESP2.OUT中可以看到所有原子的ivary都为-99,即第二步相当于什么也没做,因此最终得到的RESP电荷来自于第一步。然后看ANTECHAMBER_RESP1.OUT,输出信息为
    no.  At.no.    q(init)       q(opt)     ivary    d(rstr)/dq
    1   6       0.000000      -0.101223      0       0.003514
    2   6       0.000000      -0.148518      0       0.002793
    3   6       0.000000      -0.139285      0       0.002916
    4   6       0.000000      -0.148518      2       0.002793
    5   1       0.000000       0.136078      0       0.000000
    6   1       0.000000       0.140390      0       0.000000
    7   1       0.000000       0.147385      0       0.000000
    8   1       0.000000       0.140390      6       0.000000
    9   6       0.000000       0.404823      0       0.001199
   10   1       0.000000       0.028701      0       0.000000
   11   6       0.000000      -0.139285      3       0.002916
   12   1       0.000000       0.147385      7       0.000000
   13   6       0.000000       0.001412      0       0.004999
   14   8       0.000000      -0.469734      0       0.001041


可见结果与Multiwfn产生的略有差异,特别是第11号原子。从ivary可见,第一步拟合的时候约束了4=2、8=6、11=3、12=7。此体系结构如下
1.png

可见antechamber自动把苯环上左右对称的原子电荷约束为了相同,而Multiwfn默认情况下没有做这个约束,因此俩程序结果不同。

如果想让Multiwfn也实现同样的约束,自己写个约束设置文件即可。写个文本文件比如叫eqs.txt,内容为
4,2
8,6
11,3
12,7


在Multiwfn计算RESP电荷前,选择设置等价性约束的选项,再选读取约束文件,输入eqs.txt的路径,然后选择2开始一步式RESP电荷拟合,此时自定义的约束会被利用上,结果为
  Center      Charge
    1(C )   -0.101222
    2(C )   -0.148519
    3(C )   -0.139284
    4(C )   -0.148519
    5(H )    0.136078
    6(H )    0.140390
    7(H )    0.147385
    8(H )    0.140390
    9(C )    0.404824
   10(H )    0.028701
   11(C )   -0.139284
   12(H )    0.147385
   13(C )    0.001411
   14(O )   -0.469734
Sum of charges:    0.000000
RMSE:    0.002020   RRMSE:    0.117100


可见此时结果和Antechamber给出的精确相同了。

对苯环上对称的原子设置等价性约束并不是RESP电荷原文里提出的标准做法,纯粹是Amber的开发者自己定义的规则。这个规则有一点道理,但带来的坏处是会使得静电势重现性误差增加。如上述信息所示,不做这个约束时RMSE是0.001373,做约束后就增加到了0.002020。对于Multiwfn用户,这种约束想施加就施加,不想施加就不施加。如果懒得手动一条一条写这种等价性约束设置,在Multiwfn中可以直接利用结构的对称性自动产生等价性约束设置文件,看《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》(http://sobereva.com/441)第3.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!

14

帖子

0

威望

1168

eV
积分
1182

Level 4 (黑子)

发表于 Post on 2019-9-27 05:07:59 | 显示全部楼层 Show all
最近也注意到了这个问题。antechamber强大之处在于能把这样一个杯芳烃 calix.mol2 (16.7 KB, 下载次数 Times of downloads: 4)

4万

帖子

99

威望

4万

eV
积分
89866

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2019-9-27 05:36:45 | 显示全部楼层 Show all
锁塘柳烟池 发表于 2019-9-27 05:07
最近也注意到了这个问题。antechamber强大之处在于能把这样一个杯芳烃的结构上的等价原子都自动约束上。对 ...

自动判断对称性并施加约束是未来版本预计加入的功能

如果想基于自设的等价性约束做两步拟合,只需要按照RESP那个博文末尾说的,基于自设的等价性约束做两次一步拟合就行了,这样过程清楚,不会有任何含糊、混乱、歧义。
北京科音自然科学研究中心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!

388

帖子

1

威望

2028

eV
积分
2436

Level 5 (御坂)

发表于 Post on 2019-9-27 06:23:50 | 显示全部楼层 Show all
有一个纯属好奇的问题,如果用B3LYP-D3(BJ)/6-311G**并且带上隐式溶剂模型算RESP电荷,会被审稿人质疑吗?

4万

帖子

99

威望

4万

eV
积分
89866

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2019-9-27 07:43:51 | 显示全部楼层 Show all
Daniel_Arndt 发表于 2019-9-27 06:23
有一个纯属好奇的问题,如果用B3LYP-D3(BJ)/6-311G**并且带上隐式溶剂模型算RESP电荷,会被审稿人质疑吗?

没有可质疑的地方
顶多个别审稿人觉得基组还可以再提升提升,比如提升到def2-TZVP,对静电势的描述还是有可查觉的改进的
北京科音自然科学研究中心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!

61

帖子

0

威望

436

eV
积分
497

Level 3 能力者

发表于 Post on 2019-9-28 19:56:24 | 显示全部楼层 Show all
非常感谢老师的这篇帖子,想自己实践上尝试Multiwfn算的出RESP电荷与Antechamber的差异,但是在用高斯+amebr针对高斯的输出结果做resp电荷拟合时出现了如图所示的问题,麻烦问一下应该如何解决~
LJ7]MF69CQ{}EP]9351V4SK.png

4万

帖子

99

威望

4万

eV
积分
89866

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2019-9-29 01:59:46 | 显示全部楼层 Show all
乐嘻嘻嘻 发表于 2019-9-28 19:56
非常感谢老师的这篇帖子,想自己实践上尝试Multiwfn算的出RESP电荷与Antechamber的差异,但是在用高斯+ameb ...

从这上看不出来
北京科音自然科学研究中心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
积分
89866

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2019-10-19 13:10:39 | 显示全部楼层 Show all
锁塘柳烟池 发表于 2019-9-27 05:07
最近也注意到了这个问题。antechamber强大之处在于能把这样一个杯芳烃的结构上的等价原子都自动约束上。对 ...

今日更新了Multiwfn,目前等价性约束、电荷约束对标准RESP两步式拟合的第一步也生效了。(对第二步没法利用,因为会和标准RESP拟合的第二步在原理上冲突)

并且新版本Multiwfn已支持直接利用整体或子区域的点群对称性自动产生等价原子列表,比老版本方便多了,见此文新添加的3.5节的例子
RESP拟合静电势电荷的原理以及在Multiwfn中的计算
http://sobereva.com/441http://bbs.keinsci.com/thread-10880-1-1.html
北京科音自然科学研究中心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-2 00:43 , Processed in 0.222607 second(s), 24 queries .

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