Grompp -r --does not match the number of atoms in the topology

Ok! This should be straight forward but it’s driving me nuts. I haven’t done a lot of system building since earlier versions of gromacs. Now, I’m fighting with position restraints and the -r parameter in grompp.

My system topology file defines each individual component e.g.:

…lots of lipid, water and ion stuff…
#include “/home/ubuntu/DE_NOVO_PROTEINS/EPOC_0/0/Protein_A.itp”
#ifdef STRONG
#include “/home/ubuntu/DE_NOVO_PROTEINS/EPOC_0/0/STRONG_posre.itp”
#endif

[ system ]
Some name

[ molecules ]
…lipids, water and ions stuff…
Protein_A 4

The file Protein_A.itp defines the topology of a CG helix, 71 beads in length (1…71).
The file STRONG_posre.itp defines the x, y, and z position restraints of Protein_A.itp - also, 57 beads in length (1…71).

I am struggling with the wording in the documentation concerning the -r parameter from Grompp.

When using position restraints, a file with restraint coordinates must be supplied with -r (can be the same file as supplied for -c).

I will guess that when it says “can be the same file as supplied for -c” I’m assuming that’s on the condition that you are restraining all atoms/beads in the structure -c file? Surely, you can’t just supply the whole structure file of a protein, lipid and water when the generated restraint file only contains the atoms/beads for the protein?

So I’m running with my assumption that supplying grompp -r with a .gro structure file that only contains the restrained helix beads and numbered 1…71, as per its topology and its position restraint file. But, sadly, when I do that I get the following Warning:

WARNING 2 [file system.top, line 28]:
  The number of atoms in
  /home/ubuntu/DE_NOVO_PROTEINS/EPOC_0/0/helix_minimisation_SD.gro (71)
  does not match the number of atoms in the topology (9793). Will assume
  that the first 71 atoms in the topology and
  /home/ubuntu/DE_NOVO_PROTEINS/EPOC_0/0/helix_minimisation_SD.gro match.

This is the command:

gmx grompp -f /home/ubuntu/DE_NOVO_PROTEINS/NPT_short_eq_STRONG.mdp -c /home/ubuntu/DE_NOVO_PROTEINS/EPOC_0/0/minimization.gro -p /home/ubuntu/DE_NOVO_PROTEINS/EPOC_0/0/system.top -o /home/ubuntu/DE_NOVO_PROTEINS/EPOC_0/0/npt_eq.tpr -maxwarn 1 -r /home/ubuntu/DE_NOVO_PROTEINS/EPOC_0/0/helix_minimisation_SD.gro

Any help is much appreciated. I’ve been on Gromacs all day and I think my eyes are getting tired.
Thanks
Anthony

Hi,

position restraint files are assumed to have the same order of atoms as the original configuration. If you supply one with only 71 atoms, as the warning says it will assume that those 71 atoms are the same as the first 71 atoms in the original configuration.

If those atoms are the ones you are restraining, this is correct. If your restrained atoms are not the first 71 atoms, however, this will not work.

The way to fix this warning is to create a position restraints (-r) file with the same atoms as the original configuration, and setting the position restraint coordinates of your restrained atoms in that file. Only the restrained atom positions are important in this file, the rest are ignored.

While slightly inconvenient this does work in some way as a consistency test, that you are supplying a position restraints file that does correspond to the original.

Regards,
Petter

Hi Petter,

Thanks for the reply.

Just to confirm some details. The position restraint file (.itp) is 71 beads long (has the 71 x, y, z -fc entries). The corresponding topology it restrains (protein_A.itp) is also 71 beads long, and finally, the accompanying -r structure.gro file, also 71 beads long. All numbering 1… 71, with the corresponding entries in my system.top file:

#include protein_A.itp
#include position_restraint_file.itp

I fail to see what I’m doing wrong. Perhaps the complete system.top might reveal any errors I’ve made:

#include "../../martini_3.0.b.3.2/martini_v3.0.b.3.2.itp"
#include "../../martini_3.0.b.3.2/DPPC_M3.itp"
#include "../../martini_3.0.b.3.2/POPC_M3.itp"
#include "../../martini_3.0.b.3.2/martini_v3.0_solvents.itp"
#include "/home/ubuntu/DE_NOVO_PROTEINS/EPOC_0/0/Protein_A.itp"
#ifdef STRONG
#include "/home/ubuntu/DE_NOVO_PROTEINS/EPOC_0/0/STRONG_posre.itp"
#endif
#ifdef MEDIUM
#include "/home/ubuntu/DE_NOVO_PROTEINS/EPOC_0/0/MEDIUM_posre.itp"
#endif
#ifdef WEAK
#include "/home/ubuntu/DE_NOVO_PROTEINS/EPOC_0/0/WEAK_posre.itp"
#endif
#include "../../martini_3.0.b.3.2/martini_v3.0_ions.itp"

[ system ]
; name
KRSNMPYWFWGLLLGLLIYYLNSRK

[ molecules ]
; name  number
DPPC             132
POPC             102
WN               6465
NA               110
CL               126
Protein_A        4

Many thanks
Anthony

The position restraints file (-r) should correspond to the full system topology, not only the restrained residues.

Petter

Hi Petter,

Thank you for clarifying that point.

I’ve just grompp’d using the full topology for -r with no problems. However, I’ll do some digging as to why this was necessary in the first place, the reason escapes me.

Many thanks
Anthony

Hi @acnash,
I’m having the same sort of issue where it says, I want to restraint some part of my protein-RNA complex and have generated position restraint files for each chain individually as I have almost 11 chains and also included them in my topology description:

typ;
;	File 'topol.top' was generated
;	By user: psg (1000)
;	On host: pslab
;	At date: Wed Aug 18 17:09:36 2021
;
;	This is a standalone topology file
;
;	Created by:
;	                     :-) GROMACS - gmx pdb2gmx, 2019 (-:
;	
;	Executable:   /usr/local/gromacs/bin/gmx
;	Data prefix:  /usr/local/gromacs
;	Working dir:  /home/psg/Desktop/new_trial
;	Command line:
;	  gmx pdb2gmx -f rib_U_mutated.pdb -o rib.gro -ter
;	Force field was read from the standard GROMACS share directory.
;


; Include forcefield parameters
#include "charmm36-mar2019_psg.ff/forcefield.itp"

; include paromomycin parameters
#include "par.prm"

; Include chain topologies
#include "topol_Ion_chain_7.itp"
#include "topol_RNA_chain_A.itp"
#include "topol_RNA_chain_A2.itp"
#include "topol_RNA_chain_a.itp"
#include "topol_RNA_chain_a2.itp"
#include "topol_RNA_chain_a3.itp"
#include "topol_RNA_chain_a4.itp"
#include "topol_RNA_chain_a5.itp"
#include "topol_RNA_chain_a6.itp"
#include "topol_Protein_chain_c.itp"
#include "topol_Protein_chain_e.itp"
#include "topol_Protein_chain_i.itp"
#include "topol_Protein_chain_k.itp"
#include "topol_Protein_chain_l.itp"
#include "topol_RNA_chain_v.itp"
#include "topol_RNA_chain_x.itp"
#include "topol_RNA_chain_y.itp"

#ifdef POSRES_RES
#include "posre_RNA_chain_A_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_a_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_a2_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_a3_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_a4_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_a5_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_a6_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_a7_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_a8_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_a9_.itp"
#endif

#ifdef POSRES_RES
#include "posre_Protein_chain_c.itp"
#endif

#ifdef POSRES_RES
#include "posre_Protein_chain_e_.itp"
#endif

#ifdef POSRES_RES
#include "posre_Protein_chain_l_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_v_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_v2_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_x_.itp"
#endif

#ifdef POSRES_RES
#include "posre_RNA_chain_y_.itp"
#endif

; Include water topology
#include "charmm36-mar2019_psg.ff/tip3p.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

; Include topology for paramomycin
#include "par.itp"

; Include topology for ions
#include "charmm36-mar2019_psg.ff/ions.itp"

[ system ]
; Name
Protein in water

[ molecules ]
; Compound        #mols
Ion_chain_7         1
RNA_chain_A         1
RNA_chain_A2        1
RNA_chain_a         1
RNA_chain_a2        1
RNA_chain_a3        1
RNA_chain_a4        1
RNA_chain_a5        1
RNA_chain_a6        1
Protein_chain_c     1
Protein_chain_e     1
Protein_chain_i     1
Protein_chain_k     1
Protein_chain_l     1
RNA_chain_v         1
RNA_chain_x         1
RNA_chain_y         1
PAR		1
SOL                 9
SOL         67458
NA               224

When I issue the grompp command I get such error:
The bey32.pdb has the description of the atoms I want to restraint.

$ gmx grompp -f nvt.mdp -c em.gro -r bey32.pdb -p topol.top -o nvt.tpr
                       :-) GROMACS - gmx grompp, 2019 (-:

                            GROMACS is written by:
     Emile Apol      Rossen Apostolov      Paul Bauer     Herman J.C. Berendsen
    Par Bjelkmar      Christian Blau   Viacheslav Bolnykh     Kevin Boyd    
 Aldert van Buuren   Rudi van Drunen     Anton Feenstra       Alan Gray     
  Gerrit Groenhof     Anca Hamuraru    Vincent Hindriksen  M. Eric Irrgang  
  Aleksei Iupinov   Christoph Junghans     Joe Jordan     Dimitrios Karkoulis
    Peter Kasson        Jiri Kraus      Carsten Kutzner      Per Larsson    
  Justin A. Lemkul    Viveca Lindahl    Magnus Lundborg     Erik Marklund   
    Pascal Merz     Pieter Meulenhoff    Teemu Murtola       Szilard Pall   
    Sander Pronk      Roland Schulz      Michael Shirts    Alexey Shvetsov  
   Alfons Sijbers     Peter Tieleman      Jon Vincent      Teemu Virolainen 
 Christian Wennberg    Maarten Wolf   
                           and the project leaders:
        Mark Abraham, Berk Hess, Erik Lindahl, and David van der Spoel

Copyright (c) 1991-2000, University of Groningen, The Netherlands.
Copyright (c) 2001-2018, The GROMACS development team at
Uppsala University, Stockholm University and
the Royal Institute of Technology, Sweden.
check out http://www.gromacs.org for more information.

GROMACS is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License
as published by the Free Software Foundation; either version 2.1
of the License, or (at your option) any later version.

GROMACS:      gmx grompp, version 2019
Executable:   /usr/local/gromacs/bin/gmx
Data prefix:  /usr/local/gromacs
Working dir:  /home/psg/Desktop/new_trial
Command line:
  gmx grompp -f nvt.mdp -c em.gro -r bey32.pdb -p topol.top -o nvt.tpr

Ignoring obsolete mdp entry 'title'
Setting the LD random seed to 543556490
Generated 100032 of the 100128 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 1
Generated 65937 of the 100128 1-4 parameter combinations
Excluding 3 bonded neighbours molecule type 'Ion_chain_7'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_A'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_A2'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_a'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_a2'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_a3'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_a4'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_a5'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_a6'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'Protein_chain_c'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'Protein_chain_e'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'Protein_chain_i'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'Protein_chain_k'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'Protein_chain_l'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_v'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_x'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_y'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'PAR'
turning H bonds into constraints...
Excluding 2 bonded neighbours molecule type 'SOL'
turning H bonds into constraints...
Excluding 2 bonded neighbours molecule type 'SOL'
Excluding 1 bonded neighbours molecule type 'NA'
turning H bonds into constraints...

NOTE 1 [file topol.top, line 156]:
  System has non-zero total charge: -0.279996
  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.
  


Setting gen_seed to -1819690577
Velocities were taken from a Maxwell distribution at 300 K
Removing all charge groups because cutoff-scheme=Verlet
WARNING: all CONECT records are ignored

WARNING 1 [file topol.top, line 156]:
  The number of atoms in bey32.pdb (1944) does not match the number of
  atoms in the topology (214158). Will assume that the first 1944 atoms in
  the topology and bey32.pdb match.

WARNING: all CONECT records are ignored

WARNING 2 [file topol.top, line 156]:
  The number of atoms in bey32.pdb (1944) does not match the number of
  atoms in the topology (214158). Will assume that the first 1944 atoms in
  the topology and bey32.pdb match.


NOTE 2 [file topol.top, line 156]:
  In moleculetype 'Ion_chain_7' 32 atoms are not bound by a potential or
  constraint to any other atom in the same moleculetype. Although
  technically this might not cause issues in a simulation, this often means
  that the user forgot to add a bond/potential/constraint or put multiple
  molecules in the same moleculetype definition by mistake. Run with -v to
  get information for each atom.


NOTE 3 [file topol.top, line 156]:
  The bond in molecule-type RNA_chain_A2 between atoms 308 C2 and 309 O2
  has an estimated oscillational period of 1.9e-02 ps, which is less than
  10 times the time step of 2.0e-03 ps.
  Maybe you forgot to change the constraints mdp option.

Analysing residue names:
There are:   256        Ion residues
There are:   297      Other residues
There are:    17        RNA residues
There are:    76    Protein residues
There are: 67467      Water residues
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
Analysing Protein...
Number of degrees of freedom in T-Coupling group Protein is 3152.98
Number of degrees of freedom in T-Coupling group non-Protein is 432771.03
Determining Verlet buffer for a tolerance of 0.005 kJ/mol/ps at 300 K
Calculated rlist for 1x1 atom pair-list as 1.128 nm, buffer size 0.128 nm
Set rlist, assuming 4x4 atom pair-list, to 1.104 nm, buffer size 0.104 nm
Note that mdrun will redetermine rlist based on the actual pair-list setup

NOTE 4 [file nvt.mdp]:
  You are using a plain Coulomb cut-off, which might produce artifacts.
  You might want to consider using PME electrostatics.


This run will generate roughly 512 Mb of data

There were 4 notes

There were 2 warnings

-------------------------------------------------------
Program:     gmx grompp, version 2019
Source file: src/gromacs/gmxpreprocess/grompp.cpp (line 2309)

Fatal error:
Too many warnings (2).
If you are sure all warnings are harmless, use the -maxwarn option.

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

How did you solve the issue? Kindly help

How to go about doing that?
Say, I want to restrain the residues of the outer part of a protein-RNA complex where I have a ligand in the center and I select the residues beyond 32 Angstrom of that ligand, hence in any molecular visualisation program I select those residues and save it as a pdb file and use it as a restraint coordinate file, I’ll get the same sort of error as @acnash reported.
And in my case I have almost 7500 atoms in my original file and almost 1944 atoms which I want to restraint. How to go about creating the position restraints (-r) file with same atoms as the original configuration of my atoms? if possible please do share the trick on how to go about doing it?