Multiwfn计算的一个晶体分子的RESP电荷与Antechamber不太一致求助

Sunrise

老师好，大家好，
我最近学习使用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等文件见附件。

sobereva

这是因为你设的约束不合理。从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 )