No. of Coordinates in coordinate file does not match topology using Phosphate Ions

GROMACS version: 2021
GROMACS modification: No?

Dear all,
I am using the lysozyme tutorial with a PO4 3- Ion, instead of CL-. The error I get is when commanding:

gmx grompp -f minim.mdp -c 1AKI_solv_ions.gro -p topol.top -o em.tpr
number of coordinates in coordinate file (1AKI_solv_ions.gro, 33396)
             does not match topology (topol.top, 33830)

The walkthrough of commands was the same for the lysozyme tutorial except this line:

  • gmx genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname PO4 -nq 3 -conc 0.3

which added 186 NA ions and 62 PO4 ions. I believe this is where the coordinate issue arised from, most likely from the PO4 ions…

and adding lines within the topology for forcefields: #include “PO4ff.itp” and for the ions: #include “PO4.itp”.

I’m no advanced user at gromacs, and any help is appreciated.

The PO4ff.itp includes:

[ atomtypes ]
  opls_806  H806     1.0080     0.000    A    0.00000E+00   0.00000E+00
  opls_801  O801    15.9990     0.000    A    2.96000E-01   8.78640E-01
  opls_800  P800    30.9738     0.000    A    3.74000E-01   8.36800E-01
  opls_807  H807     1.0080     0.000    A    0.00000E+00   0.00000E+00
  opls_805  O805    15.9990     0.000    A    3.12000E-01   7.11280E-01
  opls_802  O802    15.9990     0.000    A    3.12000E-01   7.11280E-01
  opls_804  O804    15.9990     0.000    A    3.12000E-01   7.11280E-01
  opls_803  H803     1.0080     0.000    A    0.00000E+00   0.00000E+00

PO4 itp includes:

;
; GENERATED BY LigParGen Server
; Jorgensen Lab @ Yale University 
;
[ moleculetype ]
; Name               nrexcl
PO4                   3
[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  
     1   opls_800      1    UNL   P00      1     2.8980    30.9738 
     2   opls_801      1    UNL   O01      1    -1.2051    15.9990 
     3   opls_802      1    UNL   O02      1    -1.0619    15.9990 
     4   opls_803      1    UNL   H03      1     0.4728     1.0080 
     5   opls_804      1    UNL   O04      1    -1.0008    15.9990 
     6   opls_805      1    UNL   O05      1    -0.9984    15.9990 
     7   opls_806      1    UNL   H06      1     0.4508     1.0080 
     8   opls_807      1    UNL   H07      1     0.4445     1.0080 
[ bonds ]
    2     1     1      0.1480 439320.000
    3     1     1      0.1610 192464.000
    4     3     1      0.0945 462750.400
    5     1     1      0.1610 192464.000
    6     1     1      0.1610 192464.000
    7     5     1      0.0945 462750.400
    8     6     1      0.0945 462750.400

[ angles ]
;  ai    aj    ak funct            c0            c1            c2            c3 
    2     1     3     1    108.230    836.800
    1     3     4     1    108.500    460.240
    2     1     5     1    108.230    836.800
    2     1     6     1    108.230    836.800
    1     5     7     1    108.500    460.240
    1     6     8     1    108.500    460.240
    3     1     5     1    102.600    376.560
    5     1     6     1    102.600    376.560
    3     1     6     1    102.600    376.560

[ dihedrals ]
; IMPROPER DIHEDRAL ANGLES 
;  ai    aj    ak    al funct            c0            c1            c2            c3            c4            c5

[ dihedrals ]
; PROPER DIHEDRAL ANGLES
;  ai    aj    ak    al funct            c0            c1            c2            c3            c4            c5
    4    3    1    2        3       0.000   0.000   0.000  -0.000  -0.000   0.000
    7    5    1    2        3       0.000   0.000   0.000  -0.000  -0.000   0.000
    8    6    1    2        3       0.000   0.000   0.000  -0.000  -0.000   0.000
    7    5    1    3        3       0.000   0.000   0.000  -0.000  -0.000   0.000
    8    6    1    5        3       0.000   0.000   0.000  -0.000  -0.000   0.000
    8    6    1    3        3       0.000   0.000   0.000  -0.000  -0.000   0.000
    7    5    1    6        3       0.000   0.000   0.000  -0.000  -0.000   0.000
    5    1    3    4        3       0.000   0.000   0.000  -0.000  -0.000   0.000
    6    1    3    4        3       0.000   0.000   0.000  -0.000  -0.000   0.000

[ pairs ]
     2     4    1
     4     5    1
     2     7    1
     4     6    1
     3     7    1
     2     8    1
     3     8    1
     6     7    1
     5     8    1

genion can only be used to add monoatomic ions. It has no idea what the structure of PO4 is. You will have to add the PO4 molecules first, with gmx insert-molecules, then add solvent and the neutralizing Na+ ions.

Hi Jalemkul, thanks for replying!

I followed:

gmx pdb2gmx -f 1AKI_clean.pdb -o 1AKI_processed.gro -water spce
gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic
gmx insert-molecules -f 1AKI_newbox.gro -nmol 10 -ci po4.pdb -o 1AKI_boxpo4.gro

Then added #include po4ff.itp in forcefields, #include po4.itp before water topology and in the last like [molecules] PO4 10

Then

gmx solvate -cp 1AKI_boxpo4.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top
gmx grompp -f ions.mdp -c 1AKI_solv.gro -p topol.top -o ions.tpr

I still recieve a coordinate issue:

Fatal error:
number of coordinates in coordinate file (1AKI_solv.gro, 33846)
             does not match topology (topol.top, 33876)

Your topology is for fully protonated (neutral) phosphate. If you’re providing a fully deprotonated coordinate file, that’s the difference of 30 atoms (3 H on each, 10 molecules added).

There are two PO4 ions in the binding site of the protein and I was planning to keep them during the simulation. However, it seems pdb2gmx did not recognize them: Residue ‘PO4’ not found in residue topology database. Could I ask how to handle the PO4 already in PDB? I’m using charmm36-jul2022.ff. Thanks.

As far as I know, CHARMM does not support unsubstituted phosphate.

Okay, thanks. Could I ask if there is a way to parameterize it and add to the forcefield files?

The CHARMM parametrization protocols are well known and thoroughly described in the literature. I would suggest having a look at any of the protein force field and CGenFF papers.