Non-zero net charge at pdb2gmx leads to fatal error in grompp

User discussions
,
1

GROMACS version: 2022
GROMACS modification: No

Hello everyone. I’m trying to run a simple protein-in-water simulation using a heme-containing protein. I chose the CHARMM36 force field since it has the HEME group parametrization already and I’ve seen it used in MD simulations for this particular protein in the literature.

The problem I’m having is that pdb2gmx returns a non-zero net charge for the topology:

Total mass 44483.489 a.m.u.

Total charge -1.000 e

Writing topology

Now, when I proceed with the grompp command it ends with a fatal error and returns the following note and warning:

NOTE 1 [file topol.top, line 59559]:
  System has non-zero total charge: -1.000000
  Total charge should normally be an integer. See
  http://www.gromacs.org/Documentation/Floating_Point_Arithmetic
  for discussion on how close it should be to an integer.

WARNING 1 [file topol.top, line 59559]:
  You are using Ewald electrostatics in a system with net charge. This can
  lead to severe artifacts, such as ions moving into regions with low
  dielectric, due to the uniform background charge. We suggest to
  neutralize your system with counter ions, possibly in combination with a
  physiological salt concentration.

How can I solve this? Should I use the insert-molecules command to add a positive ion and reach neutral system? Is it that my initial structure is incorrect? And if so, how can I check what is causing the error?

In case it might be relevant, here is the heme information on the topology file:

residue 391 HEM rtp HEME q -2.0
  6197         FE    391    HEM     FE   1871       0.24     55.847
  6198        NPH    391    HEM     NA   1871      -0.18     14.007
  6199        NPH    391    HEM     NB   1871      -0.18     14.007
  6200        NPH    391    HEM     NC   1871      -0.18     14.007
  6201        NPH    391    HEM     ND   1871      -0.18     14.007
  6202        CPA    391    HEM    C1A   1871       0.12     12.011
  6203        CPB    391    HEM    C2A   1871      -0.06     12.011
  6204        CPB    391    HEM    C3A   1871      -0.06     12.011
  6205        CPA    391    HEM    C4A   1871       0.12     12.011
  6206        CPA    391    HEM    C1B   1871       0.12     12.011
  6207        CPB    391    HEM    C2B   1871      -0.06     12.011
  6208        CPB    391    HEM    C3B   1871      -0.06     12.011
  6209        CPA    391    HEM    C4B   1871       0.12     12.011
  6210        CPA    391    HEM    C1C   1871       0.12     12.011
  6211        CPB    391    HEM    C2C   1871      -0.06     12.011
  6212        CPB    391    HEM    C3C   1871      -0.06     12.011
  6213        CPA    391    HEM    C4C   1871       0.12     12.011
  6214        CPA    391    HEM    C1D   1871       0.12     12.011
  6215        CPB    391    HEM    C2D   1871      -0.06     12.011
  6216        CPB    391    HEM    C3D   1871      -0.06     12.011
  6217        CPA    391    HEM    C4D   1871       0.12     12.011
  6218        CPM    391    HEM    CHA   1872       -0.1     12.011
  6219         HA    391    HEM     HA   1872        0.1      1.008
  6220        CPM    391    HEM    CHB   1873       -0.1     12.011
  6221         HA    391    HEM     HB   1873        0.1      1.008
  6222        CPM    391    HEM    CHC   1874       -0.1     12.011
  6223         HA    391    HEM     HC   1874        0.1      1.008
  6224        CPM    391    HEM    CHD   1875       -0.1     12.011
  6225         HA    391    HEM     HD   1875        0.1      1.008
  6226        CT3    391    HEM    CMA   1876      -0.27     12.011
  6227        HA3    391    HEM   HMA1   1876       0.09      1.008
  6228        HA3    391    HEM   HMA2   1876       0.09      1.008
  6229        HA3    391    HEM   HMA3   1876       0.09      1.008
  6230        CT2    391    HEM    CAA   1877      -0.18     12.011
  6231        HA2    391    HEM   HAA1   1877       0.09      1.008
  6232        HA2    391    HEM   HAA2   1877       0.09      1.008
  6233        CT2    391    HEM    CBA   1878      -0.28     12.011
  6234        HA2    391    HEM   HBA1   1878       0.09      1.008
  6235        HA2    391    HEM   HBA2   1878       0.09      1.008
  6236         CC    391    HEM    CGA   1878       0.62     12.011
  6237         OC    391    HEM    O1A   1878      -0.76    15.9994
  6238         OC    391    HEM    O2A   1878      -0.76    15.9994
  6239        CT3    391    HEM    CMB   1879      -0.27     12.011
  6240        HA3    391    HEM   HMB1   1879       0.09      1.008
  6241        HA3    391    HEM   HMB2   1879       0.09      1.008
  6242        HA3    391    HEM   HMB3   1879       0.09      1.008
  6243        CE1    391    HEM    CAB   1880      -0.15     12.011
  6244        HE1    391    HEM    HAB   1880       0.15      1.008
  6245        CE2    391    HEM    CBB   1881      -0.42     12.011
  6246        HE2    391    HEM   HBB1   1881       0.21      1.008
  6247        HE2    391    HEM   HBB2   1881       0.21      1.008
  6248        CT3    391    HEM    CMC   1882      -0.27     12.011
  6249        HA3    391    HEM   HMC1   1882       0.09      1.008
  6250        HA3    391    HEM   HMC2   1882       0.09      1.008
  6251        HA3    391    HEM   HMC3   1882       0.09      1.008
  6252        CE1    391    HEM    CAC   1883      -0.15     12.011
  6253        HE1    391    HEM    HAC   1883       0.15      1.008
  6254        CE2    391    HEM    CBC   1884      -0.42     12.011
  6255        HE2    391    HEM   HBC1   1884       0.21      1.008
  6256        HE2    391    HEM   HBC2   1884       0.21      1.008
  6257        CT3    391    HEM    CMD   1885      -0.27     12.011
  6258        HA3    391    HEM   HMD1   1885       0.09      1.008
  6259        HA3    391    HEM   HMD2   1885       0.09      1.008
  6260        HA3    391    HEM   HMD3   1885       0.09      1.008
  6261        CT2    391    HEM    CAD   1886      -0.18     12.011
  6262        HA2    391    HEM   HAD1   1886       0.09      1.008
  6263        HA2    391    HEM   HAD2   1886       0.09      1.008
  6264        CT2    391    HEM    CBD   1887      -0.28     12.011
  6265        HA2    391    HEM   HBD1   1887       0.09      1.008
  6266        HA2    391    HEM   HBD2   1887       0.09      1.008
  6267         CC    391    HEM    CGD   1887       0.62     12.011
  6268         OC    391    HEM    O1D   1887      -0.76    15.9994
  6269         OC    391    HEM    O2D   1887      -0.76    15.9994   ; qtot -1
2

It is very common that proteins are not net neutral. Add a few Na and Cl ions , ideally at physiological concentration, such that the system becomes neutral.

3

Thank you for the reply. However, if I want to use the genion command to add the ions, I need a .tpr file as input, which is generated using grompp, correct? That’s the problem I have, the grompp command fails to run before I can get the .tpr file to use with genion.

4

These issues are addressed in basic tutorials. Please check out Lysozyme in Water for example.

5

Thanks for the answer. I was using the Lysozyme in water tutorial, however, going through it again in more detail I was able to figure out the problem.

Since I’m using the CHARMM36 force field I was trying to use PME electrostatics as recommended in the manual for CHARMM force fields. Looking at the tutorial ions.mdp, I switched the coulombtype option back to cutoff and was able to run the program with no issues.

Leaving the solution here in case someone has the same problem using the CHARMM force fields.

6

This has nothing to do with the CHARMM force field. The use of genion always requires an .mdp file that does not use PME. Since this input is not actually used for any real calculations, it does not need to be physically robust, merely syntactically correct.

7

Thank you very much for pointing that out and for the help.

8

I don’t know about a genion requirement of an mdp (tpr) file without PME and it seems to works fine for me with PME. Where does this non-PME requirement come in?

9

Perhaps I shouldn’t have phrased it as genion requiring that PME not be used, but grompp. There is no way generate a .tpr file for input to genion without either (1) using cutoff electrostatics or (2) using -maxwarn to override the warning about a non-zero charge in conjunction with PME. Requiring users to apply -maxwarn simply to generate an input to genion leads to bad habits, as I frequently see people haphazardly using this option in every invocation of grompp to make sure things “work.”

The only clean solutions for generating an input for genion are to use cutoff electrostatics or to create a special “integrator” or something that disables the PME net charge check internally. The system at this stage will always have a net charge; the check that is intended for preventing issues with actual dynamics aren’t relevant here, but since we co-opt the .tpr machinery to run genion, there’s always going to be some kind of conflict here.

10

That is quite annoying.
It would be good to find a better solution for this.