计算化学公社

标题: Multiwfn计算的一个晶体分子的RESP电荷与Antechamber不太一致求助 [打印本页]

作者
Author: Sunrise 时间: 2019-12-23 05:19
标题: Multiwfn计算的一个晶体分子的RESP电荷与Antechamber不太一致求助
本帖最后由 Sunrise 于 2019-12-23 05:39 编辑

老师好，大家好，
我最近学习使用Multiwfn计算RESP电荷，测试了两个药物小分子发现与Antechamber严格一致。但是在测试WIMBAV13（见下图）分子的时候，发现不能严格的重复出来，也看了sobereva老师关于计算RESP电荷的原理，知道Multiwfn的计算是合理的。因为我后续要用GAFF力场算一些晶体的性质，而力场参数是基于Antechamber计算的RESP拟合的，目前发现用不同的RESP电荷计算的一些晶面的性质有所差异，所以想学习一下是否能调试的和Antechamber一致。

Multiwfn和Antechamber计算的WIMBAV13的RESP电荷见下图，个别原子差别相对较大比如O48差别为0.08多点，H49差别为0.074，C10和C16为0.1多。我按照Antechamber的计算结果对一些原子进行了等价约束（见附件），计算采用的选项是7 ，18 ，8， 5， 1， 2，先输入的是fchk文件，期间输入了eqvcons路径和.log高斯输出文件路径。高斯输入文件采用的关键词是：#B3LYP/6-31G* SCF=Tight Pop=MK IOp(6/33=2,6/41=10,6/42=17)。相关的高斯输入、输出以及Antechamber计算的含RESP电荷的mol2等文件见附件。

作者
Author: sobereva 时间: 2019-12-23 18:01
这是因为你设的约束不合理。从antechamber产生的ANTECHAMBER_RESP1.OUT可以看到
      Point Charges Before & After Optimization

no.  At.no. q(init)    q(opt)    ivary d(rstr)/dq
1  16    0.000000    -0.277823    0    0.001693
2 6    0.000000    0.101853    0    0.003503
3 1    0.000000    0.073943    0    0.000000
4 7    0.000000    -0.551781    0    0.000892
5 1    0.000000    0.400528    0    0.000000
6 6    0.000000    0.417748    0    0.001164
7 6    0.000000    0.090544    0    0.003706
8 1    0.000000    0.104999    0    0.000000
9 1    0.000000    0.089557    0    0.000000
10 6    0.000000    -0.004008    0    0.004996
11 1    0.000000    0.060069    0    0.000000
12 1    0.000000    0.036329    0    0.000000
13 6    0.000000    -0.091853    0    0.003682
14 6    0.000000    -0.084755    0    0.003814
15 1    0.000000    0.110882    0    0.000000
16 6    0.000000    0.347442    0    0.001383
17 6    0.000000    -0.084755    14    0.003814
18 1    0.000000    0.110882    15    0.000000
19 6    0.000000    -0.300552    0    0.001578
20 6    0.000000    -0.300552    19    0.001578
21 6    0.000000    0.586801    0    0.000840
22 6    0.000000    -0.197408    0    0.002259
23 1    0.000000    0.026830    0    0.000000
24 1    0.000000    0.028316    0    0.000000
25 1    0.000000    0.031959    0    0.000000
26 6    0.000000    -0.333745    0    0.001435
27 1    0.000000    0.081791    0    0.000000
28 1    0.000000    0.073719    0    0.000000
29 1    0.000000    0.057435    0    0.000000
30 6    0.000000    -0.285126    0    0.001655
31 1    0.000000    0.051811    0    0.000000
32 1    0.000000    0.052418    0    0.000000
33 1    0.000000    0.056105    0    0.000000
34 6    0.000000    0.586801    21    0.000840
35 6    0.000000    -0.378127    0    0.001278
36 1    0.000000    0.087847    0    0.000000
37 1    0.000000    0.059704    0    0.000000
38 1    0.000000    0.098862    0    0.000000
39 6    0.000000    -0.379703    0    0.001273
40 1    0.000000    0.090020    0    0.000000
41 1    0.000000    0.093439    0    0.000000
42 1    0.000000    0.057755    0    0.000000
43 6    0.000000    -0.438870    0    0.001111
44 1    0.000000    0.093003    0    0.000000
45 1    0.000000    0.091099    0    0.000000
46 1    0.000000    0.113017    0    0.000000
47 8    0.000000    -0.498543    0    0.000983
48 8    0.000000    -0.637149    0    0.000775
49 1    0.000000    0.481244    0    0.000000

根据《关于为什么Multiwfn算的出RESP电荷与Antechamber的有所差异》http://sobereva.com/516一文的说明，可知在RESP第一步的时候只用了如下等价性约束
14,17
15,18
19,20
21,34

然而在你的eqvcons.txt里还添加了对氢的等价性约束，这是明显不当的。
如《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》（http://sobereva.com/441）一文所述，Multiwfn在两步式RESP电荷计算时，自定义等价性约束的话只会对第一步生效，而第二步还是按照标准定义去设置等价性约束。你额外定义的那些对氢的等价性约束是多余的，因为在第二步的时候自动就会去设置。

因此对当前体系若想得到和antechamber相同的结果，只需把eqvcons.txt里对氢的约束删掉，只保留
14,17
15,18
19,20
21,34

之后Multiwfn就会得到如下RESP电荷：
Center    Charge
   1(S )  -0.2778231890
   2(C ) 0.1037280767
   3(H ) 0.0718497242
   4(N )  -0.5517814749
   5(H ) 0.4005279266
   6(C ) 0.4177479925
   7(C ) 0.0872780055
   8(H ) 0.1012870788
   9(H ) 0.1012870788
10(C ) 0.0056180143
11(H ) 0.0439777580
12(H ) 0.0439777580
13(C )  -0.0918524698
14(C )  -0.0847561650
15(H ) 0.1108826540
16(C ) 0.3474409913
17(C )  -0.0847561650
18(H ) 0.1108826540
19(C )  -0.3005514376
20(C )  -0.3005514376
21(C ) 0.6273590263
22(C )  -0.2299804582
23(H ) 0.0341987928
24(H ) 0.0341987928
25(H ) 0.0341987928
26(C )  -0.3334592116
27(H ) 0.0672903316
28(H ) 0.0672903316
29(H ) 0.0672903316
30(C )  -0.3059007791
31(H ) 0.0550201219
32(H ) 0.0550201219
33(H ) 0.0550201219
34(C ) 0.4803339332
35(C )  -0.3596544838
36(H ) 0.0822477494
37(H ) 0.0822477494
38(H ) 0.0822477494
39(C )  -0.4774512365
40(H ) 0.1089353227
41(H ) 0.1089353227
42(H ) 0.1089353227
43(C )  -0.0881758764
44(H ) 0.0146290059
45(H ) 0.0146290059
46(H ) 0.0146290059
47(O )  -0.4985435057
48(O )  -0.6371490901
49(H ) 0.4812443355

可见和Antechamber的完全一样。

而且在第二步RESP拟合时你会看到Multiwfn直接自动就合理地对H做等价性约束了，用户完全没必要干预：
**** Stage 2: RESP fitting under strong hyperbolic penalty
Atom equivalence constraint imposed in this fitting stage:
Constraint 1: 8(H ) 9(H )
Constraint 2: 11(H ) 12(H )
Constraint 3: 23(H ) 24(H ) 25(H )
Constraint 4: 27(H ) 28(H ) 29(H )
Constraint 5: 31(H ) 32(H ) 33(H )
Constraint 6: 36(H ) 37(H ) 38(H )
Constraint 7: 40(H ) 41(H ) 42(H )
Constraint 8: 44(H ) 45(H ) 46(H )

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