计算化学公社

标题: 关于为什么Multiwfn算的出RESP电荷与Antechamber的有所差异 [打印本页]

作者
Author: sobereva 时间: 2019-9-27 02:32
标题: 关于为什么Multiwfn算的出RESP电荷与Antechamber的有所差异

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

Regarding why the RESP charge calculated by Multiwfn is different from that of 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。此体系结构如下

(, 下载次数 Times of downloads: 90)

可见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节的例子。

作者
Author: 锁塘柳烟池 时间: 2019-9-27 05:07
最近也注意到了这个问题。antechamber强大之处在于能把这样一个杯芳烃 (, 下载次数 Times of downloads: 4) 的结构上的等价原子都自动约束上。对我来说这个分子写一个multiwfn的约束文件有点太麻烦了。

另，resp计算有两步，我的理解multiwfn加了约束就只能做一步计算。不知道antechamber自动做约束后是做一步还是两步呢

作者
Author: sobereva 时间: 2019-9-27 05:36

锁塘柳烟池发表于 2019-9-27 05:07
最近也注意到了这个问题。antechamber强大之处在于能把这样一个杯芳烃的结构上的等价原子都自动约束上。对 ...

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

如果想基于自设的等价性约束做两步拟合，只需要按照RESP那个博文末尾说的，基于自设的等价性约束做两次一步拟合就行了，这样过程清楚，不会有任何含糊、混乱、歧义。

作者
Author: Daniel_Arndt 时间: 2019-9-27 06:23
有一个纯属好奇的问题，如果用B3LYP-D3(BJ)/6-311G**并且带上隐式溶剂模型算RESP电荷，会被审稿人质疑吗？

作者
Author: sobereva 时间: 2019-9-27 07:43

Daniel_Arndt 发表于 2019-9-27 06:23
有一个纯属好奇的问题，如果用B3LYP-D3(BJ)/6-311G**并且带上隐式溶剂模型算RESP电荷，会被审稿人质疑吗？

没有可质疑的地方
顶多个别审稿人觉得基组还可以再提升提升，比如提升到def2-TZVP，对静电势的描述还是有可查觉的改进的

作者
Author: 乐嘻嘻嘻 时间: 2019-9-28 19:56
非常感谢老师的这篇帖子，想自己实践上尝试Multiwfn算的出RESP电荷与Antechamber的差异，但是在用高斯+amebr针对高斯的输出结果做resp电荷拟合时出现了如图所示的问题，麻烦问一下应该如何解决~

作者
Author: sobereva 时间: 2019-9-29 01:59

乐嘻嘻嘻发表于 2019-9-28 19:56
非常感谢老师的这篇帖子，想自己实践上尝试Multiwfn算的出RESP电荷与Antechamber的差异，但是在用高斯+ameb ...

从这上看不出来

作者
Author: sobereva 时间: 2019-10-19 13:10

锁塘柳烟池发表于 2019-9-27 05:07
最近也注意到了这个问题。antechamber强大之处在于能把这样一个杯芳烃的结构上的等价原子都自动约束上。对 ...

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

并且新版本Multiwfn已支持直接利用整体或子区域的点群对称性自动产生等价原子列表，比老版本方便多了，见此文新添加的3.5节的例子
RESP拟合静电势电荷的原理以及在Multiwfn中的计算
http://sobereva.com/441（http://bbs.keinsci.com/thread-10880-1-1.html）

欢迎光临计算化学公社 (http://bbs.keinsci.com/)