求助QM/MM，在Gromacs调用CP2K时金属离子找不到原子序数atomic number

blah75862

执行命令：gmx grompp -f md_qmmm.mdp -c md.gro -p topol.top -o md_qmmm.tpr -n index.ndx

报错出现：
-------------------------------------------------------
Program:     gmx grompp, version 2023-EasyBuild_4.7.2.dev0_r5ec768ef07c86f5c99d1670192971a80bd72412f
Source file: src/gromacs/applied_forces/qmmm/qmmmtopologypreprocessor.cpp (line 282)

Fatal error:
Atoms 6126 does not have atomic number needed for QMMM. Check atomtypes
section in your topology or forcefield.

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------


目的是解决这个报错，以下是详情。

详细信息：
md.gro --部分：
  444GLY    OC1 6120   5.943   7.349   7.315  0.3482  0.1903 -0.2700
  444GLY    OC2 6121   6.111   7.447   7.197 -0.5262 -0.4131 -0.2818
  500ZN      ZN 6122   4.594   6.198   5.794 -0.3117  0.0016  0.4411
  501ZN      ZN 6123   3.831   4.887   6.583 -0.3727 -0.1098 -0.1493
  502CA      CA 6124   3.950   4.982   5.326  0.0528  0.1777  0.1857
  503CA      CA 6125   4.622   4.244   5.181  0.4271  0.0206 -0.5234
  504CA      CA 6126   4.298   5.365   7.396 -0.3069 -0.0691  0.0228
  505UNL      N 6127   3.807   4.195   4.600 -0.2368 -0.1051 -0.8978
  505UNL     N1 6128   4.257   4.045   4.451  0.4084 -0.1262  0.8145


一些问题的思考与描述：
其中444号残基及之前的是蛋白，501-505是离子，505及之后的是配体。其中与配体相连的是503的Ca2+，所以残基序数为504和501这两个钙离子并不在QM范围中，而503号残基的6125号原子在QM范围中。当gmx grompp找不到原子序数的是504（挺奇怪的）因为这三个钙离子本质上只有坐标的差别，并不存在在itp和topol.top之间的各种引用之间出问题。所以三个钙离子中只有最后一个出现原子序数问题就让我非常费解。
因为前序的时候，作为含有金属离子的蛋白质与小分子配体模拟的话，已经按照Gromacs培训班的内容，在对应拓扑文件下完成了EM PR MD，所以力场索引时的打包不应该出现问题。

也许是我忽略了某些细节，在此尽量简洁地附上有关的文件之局部：

-------------------------------------------------------

topol.top:
; Include forcefield parameters
#include "amber99sb-ildn.ff/forcefield.itp"

; Include chain topologies
#include "ligand_GMX.itp"
#include "protein_Protein_chain_A.itp"
#include "protein_Ion_chain_A2.itp"

; Include topology for ions
#include "amber99sb-ildn.ff/ions.itp"


[ molecules ]
; Compound        #mols
Protein_chain_A     1
Ion_chain_A2        1
ligand             1         

-------------------------------------------------------
-------------------------------------------------------
protein_Ion_chain_A2.itp：


[ moleculetype ]
; Name            nrexcl
Ion_chain_A2        3

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB      massB
; residue 500 ZN  rtp ZN   q +2.0
     1         Zn    500     ZN     ZN      1          2       65.4   ; qtot 2
; residue 501 ZN  rtp ZN   q +2.0
     2         Zn    501     ZN     ZN      2          2       65.4   ; qtot 4
; residue 502 CA  rtp CA   q +2.0
     3         C0    502     CA     CA      3          2      40.08   ; qtot 6
; residue 503 CA  rtp CA   q +2.0
     4         C0    503     CA     CA      4          2      40.08   ; qtot 8
; residue 504 CA  rtp CA   q +2.0
     5         C0    504     CA     CA      5          2      40.08   ; qtot 10

-------------------------------------------------------
-------------------------------------------------------
ions.itp：

[ moleculetype ]
; molname       nrexcl
CA              1

[ atoms ]
; id    at type         res nr  residu name     at name  cg nr  charge
1       C0              1       CA              CA       1      2.00000


-------------------------------------------------------
-------------------------------------------------------
ffnonbonded.itp：

[ atomtypes ]
; name      at.num  mass     charge ptype  sigma      epsilon
C0          20      40.08    0.0000  A   3.05240e-01  1.92376e+00



-------------------------------------------------------

求助关于这一找原子序数的解决办法。谢谢。

PS 看到同问题的另一位发帖博主遇到的问题似乎还不一样，那位朋友转去纯用CP2K了，所以也就没有再PM)



4.13 12:08 追加：

我又测试了另外一个体系，.gro文件的相关段落如下：
  242GLY      C 4760   5.517   6.424   8.200 -0.0493 -0.1042 -0.9216
  242GLY    OC1 4761   5.584   6.510   8.143  0.8850  0.2375  0.4917
  242GLY    OC2 4762   5.468   6.441   8.315  0.5356 -0.3495  0.1248
  996CA      CA 4763   4.875   3.319   2.432 -0.2133  0.2644 -0.0913
  997CA      CA 4764   2.858   3.920   2.379  0.3885 -0.0312 -0.2469
  998ZN      ZN 4765   3.949   3.718   1.950 -0.2158 -0.1123  0.0862
  999ZN      ZN 4766   3.382   3.195   2.857  0.1799 -0.0062 -0.0012
  996CA      CA 4767   4.676   7.417   5.963 -0.2843  0.0728  0.2189
  997CA      CA 4768   5.774   6.037   4.583 -0.1622 -0.1721  0.0064
  998ZN      ZN 4769   5.475   7.231   5.263  0.1066  0.0747 -0.3438
  999ZN      ZN 4770   5.772   6.296   5.737 -0.2014  0.0324  0.1067
1000UNK    H10 4771   3.372   4.284   4.956  1.7591 -1.2583  0.6062
1000UNK    H22 4772   4.072   3.899   4.915 -1.9854 -0.4021 -0.8280
1000UNK    H02 4773   3.672   4.634   4.565 -1.0338 -0.2016  2.0067

-------------------------------------------------------
Program:     gmx grompp, version 2023-EasyBuild_4.7.2.dev0_r5ec768ef07c86f5c99d1670192971a80bd72412f
Source file: src/gromacs/applied_forces/qmmm/qmmmtopologypreprocessor.cpp (line 282)

Fatal error:
Atoms 4770 does not have atomic number needed for QMMM. Check atomtypes
section in your topology or forcefield.

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------

因为金属离子出问题的一直是主体蛋白质部分的最后一个，倾向于认为是程序判定时处理最后一个离子文件出了问题。

kexin

遇到了类似的问题，请问楼主最后怎么解决的？