Gromacs+CP2K, cannot find atomic number

GROMACS version:2023
GROMACS modification: No
Force filed: Amber99SB-ildn, GAFF
CP2K version: 2023.1

Protein-ligand system, want to use QM for ligand and contacting residue, MM for the rest.

Command: gmx grompp -f md_qmmm.mdp -c md.gro -p topol.top -o md_qmmm.tpr -n index.ndx
Msg prompting:

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.

Here are some details
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

Some thoughts and description of the problem:

Where residues 444 and before are proteins, 501-505 are ions, and 505 and after are ligands etc. The one Connected to the ligand is Ca2+ of 503, so residues number 504 and 501 which are two calcium ions are not in the QM range, while atom 6125 of residue 503 is in the QM range. When gmx grompp can not find the atomic number is 504 (quite strange) because these three calcium ions are essentially only the difference in coordinates, and there is no problem between the various references between itp and topol.top. So only the last one of the three calcium ions has an atomic number problem, which is very puzzling to me.

Because the preorder, as a protein containing metal ions and small molecule ligands simulation, has been completed in accordance with the contents of the Gromacs training course, EM PR MD under the corresponding topology file, so the force field index when the packaging should not have problems.

Perhaps I have overlooked some details, and I will try to attach here as concisely as possible the relevant parts of the files:


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


Also, I test another system, whose .gro is:
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.

Can someone help me locate the problem?