生成18碳环的itp拓扑时原子电荷的计算问题以及一些思考

不想飞的猫头鹰

      最近在回顾sob老师经典博文http://sobereva.com/441《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》的18碳环的RESP电荷计算部分时，我突然注意到了之前没有思考过的问题，18碳环中的所有碳原子按理说应该是化学环境等价的，因此需要设置等价性约束，这种情况下，常规的RESP电荷算出来必然为0，必须拟合出非原子中心的电荷，才能还原静电势分布。
1：然而，gromacs的itp分子拓扑文件中的原子电荷只能指定到实际的原子上，也就是只能在18碳环的itp文件中只能填入18个电荷数值。而按照Multiwfn手册的4.7.7.6节中的例子，产生了非原子中心位点的另外18个虚原子的电荷，这时候总共有36个位点，在itp中无法显示出来，这种情况下计算的RESP电荷又怎么用于分子动力学模拟？
2：我查看了社长对18碳环的研究中http://sobereva.com/674《8字形双环分子对18碳环的独特吸附行为的量子化学、波函数分析与分子动力学研究》中的原文https://doi.org/10.1039/D3CP01896B的SI补充信息时注意到有“Only van der Waals (vdW) interaction was taken into account because atomic charges of C18 are all close to zero and thus electrostatic interaction cannot be represented.”和“Because atomic charges of C18 are all zero due to its structure symmetry, electrostatic interaction was not taken into account in the simulation.”也就是说，在这种以范德华作用主导的超分子复合物体系中，直接把18碳环的碳原子设置为0没什么问题，因为文中已经用IGMH方法证实分子间相互作用的本质是范德华作用，静电势穿透图里穿透程度较低、范德华势更好地重合也说明了范德华作用的主导。然而，我又注意到在http://sobereva.com/572《探究18碳环独特的分子间相互作用与pi-pi堆积特征》研究18碳环二聚体堆积动力学模拟中，是跑的aimd，这应该是考虑到博文中提到的静电势互补的稳定化作用以及SAPT的结果里“静电作用相对次要，只贡献了19%”这百分之19仍不能直接忽略的影响，所以为了动态反映电子结构的变化才用的aimd。然而，博文中的静电势穿透程度仍然较小，碳环中心也有较大的范德华势相互重叠区域，虽然IGMH等值面仍然是显著的绿色来说明作用本质是范德华作用，但是EDA仍有19%的静电相互作用，我可不可以这样理解：EDA的数值只是代表能量的数值大小，但是并不能代表对应的能量项目的梯度也就是受力的本质，能量的数值跟势函数的形式有很大的关系，EDA的百分比并不能说明相互作用中受力种类的百分比，作用本质仍应以IGMH的结果为主。这种理解是否有问题？EDA的结果多大程度上说明了静电相互作用的不可忽略性？
3：在gromacs的itp中即使像研究OPP与碳环的吸附那样，指定碳环中的碳原子的电荷为0，这种情况下去跑两个18碳环的堆积动力学轨迹，由于范德华做作用，上下层碳原子也确实尽可能交替填充，但是失去了静电相互作用，这时候的分子拓扑是否严格？假如所研究的体系是含有极性很大的分子与18碳环，不能忽略静电相互作用，那么这时候就涉及到了指定电荷的问题，真想表现出来18碳环静电势正负交替的情况，碳原子的电荷就必然是有正有负，这时候计算RESP电荷时加上全部一样的等价性约束肯定不行了，应该设置什么样的等价性约束？对于wb97xd/tzvp级别优化的18碳环的波函数，我发现如果分别对奇数和偶数号的碳原子分别设置约束，就会是正负交替的，
采用MK方法时，奇偶数号碳原子电荷数值为-0.0000845294和0.0000845294，RMSE:    0.001870   RRMSE:    1.000000，
采用CHELPG方法时，奇偶数号碳原子电荷数值为0.0024948230和-0.0024948230，RMSE:    0.007195   RRMSE:    0.999994。
不采用任何约束时，基于MK方法的RESP电荷的计算结果是
     1(C )   0.0011125877
     2(C )   0.0006207184
     3(C )  -0.0022456861
     4(C )   0.0042459376
     5(C )  -0.0068684581
     6(C )   0.0068936429
     7(C )  -0.0047019671
     8(C )   0.0027870279
     9(C )  -0.0005592424
    10(C )  -0.0015922826
    11(C )   0.0032796759
    12(C )  -0.0036786178
    13(C )   0.0035855858
    14(C )  -0.0030841709
    15(C )   0.0023498013
    16(C )  -0.0023405057
    17(C )   0.0031048185
    18(C )  -0.0029088653
Sum of charges:  -0.0000000000
RMSE:    0.001870   RRMSE:    0.999877
基于CHELPG方法的RESP电荷的计算结果是
     1(C )  -0.0050911386
     2(C )  -0.0001456799
     3(C )   0.0055277056
     4(C )  -0.0062412199
     5(C )   0.0065780819
     6(C )  -0.0086746366
     7(C )   0.0094692001
     8(C )  -0.0113333929
     9(C )   0.0181280227
    10(C )  -0.0191679231
    11(C )   0.0079715884
    12(C )   0.0026762764
    13(C )  -0.0080403489
    14(C )   0.0069590338
    15(C )  -0.0034468686
    16(C )   0.0042173591
    17(C )  -0.0067151617
    18(C )   0.0073291022
Sum of charges:  -0.0000000000
RMSE:    0.007193   RRMSE:    0.999717
这里又产生一个问题，MK和CHELPG方法产生的正负顺序不一致，这时候应该相信谁的结果？
另外，同样可以用于分子动力学模拟的ADCH电荷计算时不用设置等价性约束，这时候的计算结果是
Atom    1(C ):     0.00000239
Atom    2(C ):     0.00000239
Atom    3(C ):     0.00008788
Atom    4(C ):    -0.00005450
Atom    5(C ):    -0.00005140
Atom    6(C ):     0.00002691
Atom    7(C ):    -0.00005513
Atom    8(C ):     0.00000018
Atom    9(C ):     0.00005848
Atom   10(C ):    -0.00001482
Atom   11(C ):    -0.00001482
Atom   12(C ):     0.00005848
Atom   13(C ):     0.00000018
Atom   14(C ):    -0.00005513
Atom   15(C ):     0.00002691
Atom   16(C ):    -0.00005140
Atom   17(C ):    -0.00005450
Atom   18(C ):     0.00008788
可见，ADCH电荷布居分析时并没有把碳原子判断为等价性的。
4.这个问题的根源在于，gromacs的itp没法反映非原子中心的电荷的能力，或者说，这是分子力学没有办法通过原子电荷来反应电四极矩之类的更细致的电子结构的问题。那么，在当前的条件下，如果硬要用基于力场的分子动力学去跑像18碳环这样的体系的分子动力学模拟，应该用什么电荷？还是索性把像是18碳环这样的范德华作用主导的分子中的原子电荷本身就很小的原子电荷直接设置为0？

感谢指点和不吝赐教！

sobereva

C18的那点静电作用，对多数情况的常温下的动力学行为不会有显著影响，对于下面两篇文章研究的情况就是如此，起码造成的误差比其它方面（如范德华作用描述误差、C18结构柔性描述误差、OPP部分描述的误差等）还小。如果说C18都必须靠虚拟点来描述其静电分布不均衡的特征的话，那文中的OPP部分也该这么考虑，那就太麻烦了
理论设计新颖的基于18碳环构成的双马达超分子体系
http://sobereva.com/684（http://bbs.keinsci.com/thread-39132-1-1.html）
8字形双环分子对18碳环的独特吸附行为的量子化学、波函数分析与分子动力学研究
http://sobereva.com/674（http://bbs.keinsci.com/thread-38312-1-1.html）

只有对动力学精度要求较高、需要描述得很精细，比如Carbon, 171, 514-523 (2021)这篇文章里通过动力学考察18碳环二聚体，那才需要考虑静电作用的影响。如果不描述出静电作用，两个18碳环都不会呈现长键对着短键的交错结构。


所有那些基于数值格点的方法计算的原子电荷不可能算出来18个碳的电荷正好相同，因为实际用的格点不是无穷精细。增加格点密度可以增加等价性。

不想飞的猫头鹰

sobereva 发表于 2024-4-21 04:19
C18的那点静电作用，对多数情况的常温下的动力学行为不会有显著影响，对于下面两篇文章研究的情况就是如此 ...

感谢sob老师，之前没手写过虚拟位点，用的都是直接下载好的文件中的虚拟位点比如martini2粗粒化力场中的胆固醇中的虚拟位点来做铰链结构，我今天白天不断尝试所以费了不少时间，以下是我的尝试记录和附带的文件：
首先从gromacs官方手册里找到虚拟位点相关的说明：
https://manual.gromacs.org/current/reference-manual/topologies/topology-file-formats.html#id36
https://manual.gromacs.org/current/reference-manual/functions/interaction-methods.html#virtualsites

我按照Multiwfn手册中的计算18碳环非原子中心的RESP电荷的方法，得知了短键中心处的虚原子带负电，长键中心处的虚原子带正电，碳原子本身带微量正电。我的wb97xd/tzvp的fchk文件中，1-2号是短键，2-3号是长键，以此类推，19 21 23...35号虚原子放到短键中心，itp中该类虚原子类型设置为MA，20 22 24...36号虚原子放到长键中心，itp中该类虚原子类型设置为MB，虚原子的原子名均加上V前缀，序号和电荷组都依照数字顺延，均属于同一残基。最后用Multiwfn转换gif为gro并替换原子名。这样，[ atoms ]里就有了36个可分配电荷的位点，我就可以把sob老师例子里的电荷指定上了，下面是我的[ atoms ]字段：
[ atoms ]
;  Index   type   residue  resname   atom        cgnr     charge       mass
     1     c1         1      MOL     C1            1    0.0642410906   12.010736
     2     c1         1      MOL     C2            2    0.0642410906   12.010736
     3     c1         1      MOL     C3            3    0.0642410906   12.010736
     4     c1         1      MOL     C4            4    0.0642410906   12.010736
     5     c1         1      MOL     C5            5    0.0642410906   12.010736
     6     c1         1      MOL     C6            6    0.0642410906   12.010736
     7     c1         1      MOL     C7            7    0.0642410906   12.010736
     8     c1         1      MOL     C8            8    0.0642410906   12.010736
     9     c1         1      MOL     C9            9    0.0642410906   12.010736
    10     c1         1      MOL     C10          10    0.0642410906   12.010736
    11     c1         1      MOL     C11          11    0.0642410906   12.010736
    12     c1         1      MOL     C12          12    0.0642410906   12.010736
    13     c1         1      MOL     C13          13    0.0642410906   12.010736
    14     c1         1      MOL     C14          14    0.0642410906   12.010736
    15     c1         1      MOL     C15          15    0.0642410906   12.010736
    16     c1         1      MOL     C16          16    0.0642410906   12.010736
    17     c1         1      MOL     C17          17    0.0642410906   12.010736
    18     c1         1      MOL     C18          18    0.0642410906   12.010736
    19     MA         1      MOL     V19          19    -0.503822021   0.00000
    20     MB         1      MOL     V20          20    0.3753398397   0.00000
    21     MA         1      MOL     V21          21    -0.503822021   0.00000
    22     MB         1      MOL     V22          22    0.3753398397   0.00000
    23     MA         1      MOL     V23          23    -0.503822021   0.00000
    24     MB         1      MOL     V24          24    0.3753398397   0.00000
    25     MA         1      MOL     V25          25    -0.503822021   0.00000
    26     MB         1      MOL     V26          26    0.3753398397   0.00000
    27     MA         1      MOL     V27          27    -0.503822021   0.00000
    28     MB         1      MOL     V28          28    0.3753398397   0.00000
    29     MA         1      MOL     V29          29    -0.503822021   0.00000
    30     MB         1      MOL     V30          30    0.3753398397   0.00000
    31     MA         1      MOL     V31          31    -0.503822021   0.00000
    32     MB         1      MOL     V32          32    0.3753398397   0.00000
    33     MA         1      MOL     V33          33    -0.503822021   0.00000
    34     MB         1      MOL     V34          34    0.3753398397   0.00000
    35     MA         1      MOL     V35          35    -0.503822021   0.00000
    36     MB         1      MOL     V36          36    0.3753398397   0.00000

然后，根据上面的虚拟位点的添加方式，在itp中添加了[ virtual_sites2 ]字段，把19-36这18个虚拟位点添加上使用funct1函数，指定和相邻原子的距离，参数a是从gview里量取短键和长键的键长后，除以10再除以2变成gromacs纳米单位的键长的一半，下面是我的写法：
[ virtual_sites2 ]
; Vsite from           funct   a      
19       1       2     1       0.06106
20       2       3     1       0.067187
21       3       4     1       0.06106
22       4       5     1       0.067187
23       5       6     1       0.06106
24       6       7     1       0.067187
25       7       8     1       0.06106
26       8       9     1       0.067187
27       9       10     1       0.06106
28       10       11     1       0.067187
29       11       12     1       0.06106
30       12       13     1       0.067187
31       13       14     1       0.06106
32       14       15     1       0.067187
33       15       16     1       0.06106
34       16       17     1       0.067187
35       17       18     1       0.06106
36       18       1     1       0.067187

然后我想，由于这时候指定上了很多较为明显的净电荷，而且碳环中实原子和虚原子的位置很近，如果计算了分子内非键相互作用，就容易在模拟中由于键振动时导致静电相互作用过大使模拟崩溃（我之后尝试了50ps真空中单碳环300k的nvt模拟，1fs步长，跑到后面果然是由于势能NaN崩溃了），所以应该排除分子内静电相互作用的计算。这时候我就想仿照粗粒化力场拓扑martini_v2.0_CHOL_02.itp中的写法，添加[ exclusions ]项，想来想去，必须把所有原子彼此间的项都加上，因此我写成了以下这样：
[ exclusions ]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
23 24 25 26 27 28 29 30 31 32 33 34 35 36
24 25 26 27 28 29 30 31 32 33 34 35 36
25 26 27 28 29 30 31 32 33 34 35 36
26 27 28 29 30 31 32 33 34 35 36
27 28 29 30 31 32 33 34 35 36
28 29 30 31 32 33 34 35 36
29 30 31 32 33 34 35 36
30 31 32 33 34 35 36
31 32 33 34 35 36
32 33 34 35 36
33 34 35 36
34 35 36
35 36

然后，我生成了tpr，并进行模拟，发现一开始好像还不错，但后期振动模式越来越诡异，甚至原子坍塌成一团了。不清楚这是什么原因导致的？我自己分析有两点，要么是我的mdp文件中哪里写的有问题，要么是上面exclusion字段虽然避免了分子内静电相互作用的计算，但是也导致了nrexcl以上的分子内范德华相互作用的计算被忽略，而这些范德华力对于这样一个不大的分子环张力的表现时不可忽略的。思来想去，好像没有办法在gromacs中只计算分子内范德华相互作用而不计算分子内非0电荷原子间的静电相互作用。
没有别的方法了，然后我想，按照sob老师的提醒，以常规方式计算resp电荷或者adch电荷或者1.2*cm5电荷的数值都那么小了，而且在md过程中热波动和各种随机效应的影响下，就算没法用18个实原子的点电荷反映出36个实原子加虚拟位点混合的静电势，但是动力学行为应该基本差不多，即使有些原子的吸附位点不那么准也会被动态过程模糊掉，所以要是实在没有更好的加上虚拟位点又能够让模拟稳定的方法，直接索性拿旋转不变性更好的CHELPG，配合拟合格点间距更小的无等价性约束的RESP电荷分配给这18个碳原子，是否合理？我把格点间距直接索性设置到了0.05Bohr，2696v3 18核算了将近两个小时，结果如下：
   Center       Charge
     1(C )  -0.0001012264
     2(C )   0.0000797468
     3(C )  -0.0000459640
     4(C )   0.0000209014
     5(C )   0.0000157849
     6(C )  -0.0000430119
     7(C )   0.0000596121
     8(C )  -0.0000626730
     9(C )   0.0000589037
    10(C )  -0.0000545314
    11(C )   0.0000370691
    12(C )  -0.0000309145
    13(C )   0.0000244121
    14(C )   0.0000001841
    15(C )  -0.0000334498
    16(C )   0.0000646071
    17(C )  -0.0000933236
    18(C )   0.0001038733
Sum of charges:  -0.0000000000
RMSE:    0.007211   RRMSE:    1.000000

可见常规方式计算的RESP电荷数值确实不大，而非原子中心拟合的方法得到的数值大了3个数量级以上，把虚原子的电荷加到itp里确实容易导致静电相互作用太大。

yxdd98

不想飞的猫头鹰 发表于 2024-4-21 20:14
感谢sob老师，之前没手写过虚拟位点，用的都是直接下载好的文件中的虚拟位点比如martini2粗粒化力场中的 ...

主要问题是虚拟位点的定义，virtual_sites2，使用funct 1时，a指的是两个连线原子的权重系数，比如在你的体系里，虚拟位点位于碳碳键中部的话，a就取0.5就行，你按实际距离来设置虚拟点，就会导致重叠了。
跑了一下，发现18碳环会发生各种奇怪的扭曲构象，不清楚itp里的键参数是否有改进的地方

不想飞的猫头鹰

yxdd98 发表于 2024-7-18 10:17
主要问题是虚拟位点的定义，virtual_sites2，使用funct 1时，a指的是两个连线原子的权重系数，比如在你的 ...

感谢纠正谬误，我刚才把a都调到0.5了，又发现之前上传的nvtmd.mdp写法有点问题，后面是冻结分子的设置，手动删去并把控温组改成system，消除平动，其他设置不变。我之前试过了sobtop中不同的基于hessian生成力常数的方法，结合着虚拟位点这些电荷都会产生很诡异的运动模式。我刚才想到，既然想要维持碳环的刚性和平面性，而不希望它扭曲，同时体现出局部电荷特征，我就索性把所有键角和二面角力常数调大了10倍，然后再在真空中跑nvtmd，步长仍然是2fs,a使用原来的数值或者改正后的0.5都尝试了，发现这次带着虚拟位点电荷也不会使结构有不稳定的问题了。虽然有些随意性和人为性，但是至少能带着电荷模拟了。