Common errors when using GROMACS

The vast majority of error messages generated by GROMACS are descriptive, informing the user where the exact error lies. Some errors that arise are noted below, along with more details on what the issue is and how to solve it.

Common errors during usage

Out of memory when allocating

The program has attempted to assign memory to be used in the calculation, but is unable to due to insufficient memory.

Possible solutions are:

  • reduce the scope of the number of atoms selected for analysis.
  • reduce the length of trajectory file being processed.
  • in some cases confusion between Ångström and nm may lead to users generating a pdb2gmx water box that is 103 times larger than what they think it is (e.g. gmx solvate).
  • use a computer with more memory.
  • install more memory in the computer.

The user should bear in mind that the cost in time and/or memory for various activities will scale with the number of atoms/groups/residues N or the simulation length T as order N, NlogN, or N2 (or maybe worse!) and the same for T, depending on the type of activity. If it takes a long time, have a think about what you are doing, and the underlying algorithm (see the Reference manual, man page, or use the -h flag for the utility), and see if there’s something sensible you can do that has better scaling properties.

Errors in pdb2gmx

Residue ‘XXX’ not found in residue topology database

This means that the force field you have selected while running pdb2gmx does not have an entry in the residue database for XXX. The residue database entry is necessary both for stand-alone molecules (e.g. formaldehyde) or a peptide (standard or non-standard). This entry defines the atom types, connectivity, bonded and non-bonded interaction types for the residue and is necessary to use pdb2gmx to build a top file. A residue database entry may be missing simply because the database does not contain the residue at all, or because the name is different.

For new users, this error appears because they are running pdb2gmx on a PDB file they have, without consideration of the contents of the file. A force field is not magical, it can only deal with molecules or residues (building blocks) that are provided in the residue database or included otherwise.

If you want to use pdb2gmx to automatically generate your topology, you have to ensure that the appropriate rtp entry is present within the desired force field and has the same name as the building block you are trying to use. If you call your molecule “HIS,” then pdb2gmx will try to build histidine, based on the [ HIS ] entry in the rtp file, so it will look for the exact atomic entries for histidine, no more no less.

If you want a topology for an arbitrary molecule, you cannot use pdb2gmx (unless you build the rtp entry yourself). You will have to build that entry by hand, or use another program (such as x2top or one of the scripts contributed by users) to build the top file.

If there is not an entry for this residue in the database, then the options for obtaining the force field parameters are:

  • see if there is a different name being used for the residue in the residue database and rename as appropriate,
  • parameterize the residue / molecule yourself (lots of work, even for an expert),
  • find a topology file for the molecule, convert it to an itp file and include it in your top file,
  • use another force field which has parameters available for this,
  • search the primary literature for publications for parameters for the residue that are consistent with the force field that is being used.

Long bonds and/or missing atoms

There are probably atoms missing earlier in the pdb file which makes pdb2gmx go crazy. Check the screen output of pdb2gmx, as it will tell you which one is missing. Then add the atoms in your pdb file, energy minimization will put them in the right place, or fix the side chain with e.g. the WHAT IF program.

Chain identifier ‘X’ was used in two non-sequential blocks

This means that within the coordinate file fed to pdb2gmx, the X chain has been split, possibly by the incorrect insertion of one molecule within another. The solution is simple: move the inserted molecule to a location within the file so that it is not splitting another molecule. This message may also mean that the same chain identifier has been used for two separate chains. In that case, rename the second chain to a unique identifier.

WARNING: atom X is missing in residue XXX Y in the pdb file

Related to the long bonds/missing atoms error above, this error is usually quite obvious in its meaning. That is, pdb2gmx expects certain atoms within the given residue, based on the entries in the force field rtp file. There are several cases to which this error applies:

  • Missing hydrogen atoms; the error message may be suggesting that an entry in the hdb file is missing. More likely, the nomenclature of your hydrogen atoms simply does not match what is expected by the rtp entry. In this case, use -ignh to allow pdb2gmx to add the correct hydrogens for you, or re-name the problematic atoms.
  • A terminal residue (usually the N-terminus) is missing H atoms; this usually suggests that the proper -ter option has not been supplied or chosen properly. In the case of the AMBER force fields, nomenclature is typically the problem. N-terminal and C-terminal residues must be prefixed by N and C, respectively. For example, an N-terminal alanine should not be listed in the pdb file as ALA, but rather NALA, as specified in the ffamber instructions.
  • Atoms are simply missing in the structure file provided to pdb2gmx; look for REMARK 465 and REMARK 470 entries in the pdb file. These atoms will have to be modeled in using external software. There is no GROMACS tool to re-construct incomplete models.

Contrary to what the error message says, the use of the option -missing is almost always inappropriate. The -missing option should only be used to generate specialized topologies for amino acid-like molecules to take advantage of rtp entries. If you find yourself using -missing in order to generate a topology for a protein or nucleic acid, don’t; the topology produced is likely physically unrealistic.

Atom X in residue YYY not found in rtp entry

If you are attempting to assemble a topology using pdb2gmx, the atom names are expected to match those found in the rtp file that define the building block(s) in your structure. In most cases, the problem arises from a naming mismatch, so simply re-name the atoms in your coordinate file appropriately. In other cases, you may be supplying a structure that has residues that do not conform to the expectations of the force field, in which case you should investigate why such a difference is occurring and make a decision based on what you find - use a different force field, manually edit the structure, etc.

No force fields found (files with name ‘forcefield.itp’ in subdirectories ending on ‘.ff’)

This means your environment is not configured to use GROMACS properly, because pdb2gmx cannot find its databases of forcefield information. This could happen because a GROMACS installation was moved from one location to another. Either follow the instructions about getting access to |Gromacs| after installation or re-install GROMACS before doing so.

Errors in grompp

Found a second defaults directive file

This is caused by the [defaults] directive appearing more than once in the topology or force field files for the system - it can only appear once. A typical cause of this is a second defaults being set in an included topology file, itp, that has been sourced from somewhere else. For specifications on how the topology files work, see the reference manual, Section 5.6.:

[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1       1       no       1.0       1.0

One solution is to simply comment out (or delete) the lines of code out in the file where it is included for the second time i.e.,:

;[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
;1       1       no       1.0       1.0

A better approach to finding a solution is to re-think what you are doing. The [defaults] directive should only be appearing at the top of your top file where you choose the force field. If you are trying to mix two force fields, then you are asking for trouble. If a molecule itp file tries to choose a force field, then whoever produced it is asking for trouble.

Invalid order for directive xxx

The directives in the .top and .itp files have rules about the order in which they can appear, and this error is seen when the order is violated. Consider the examples and discussion in chapter 5 of the reference manual, and/or from tutorial material. The include file mechanism cannot be used to #include a file in just any old location, because they contain directives and these have to be properly placed.

In particular, Invalid order for directive defaults is a result of defaults being set in the topology or force field files in the inappropriate location; the [defaults] section can only appear once and must be the first directive in the topology. The [defaults] directive is typically present in the force field file (forcefield.itp), and is added to the topology when you #include this file in the system topology.

If the directive in question is [atomtypes] (which is the most common source of this error) or any other bonded or nonbonded [*types] directive, typically the user is adding some non-standard species (ligand, solvent, etc) that introduces new atom types or parameters into the system. As indicated above, these new types and parameters must appear before any [moleculetype] directive. The force field has to be fully constructed before any molecules can be defined.

Atom index n in position_restraints out of bounds

A common problem is placing position restraint files for multiple molecules out of order. Recall that a position restraint itp file containing a [ position_restraints ] block can only belong to the [ moleculetype ] block that contains it. For example:

WRONG:

#include "topol_A.itp"
#include "topol_B.itp"
#include "ligand.itp"

#ifdef POSRES
#include "posre_A.itp"
#include "posre_B.itp"
#include "ligand_posre.itp"
#endif

RIGHT:

#include "topol_A.itp"
#ifdef POSRES
#include "posre_A.itp"
#endif

#include "topol_B.itp"
#ifdef POSRES
#include "posre_B.itp"
#endif

#include "ligand.itp"
#ifdef POSRES
#include "ligand_posre.itp"
#endif

Further, the atom index of each [position_restraint] must be relative to the [moleculetype], not relative to the system (because the parsing has not reached [molecules] yet, there is no such concept as “system”). So you cannot use the output of a tool like genrestr blindly (as genrestr -h warns).

System has non-zero total charge

Notifies you that counter-ions may be required for the system to neutralize the charge or there may be problems with the topology.

If the charge is not very close to an integer, then this indicates that there is a problem with the topology. If pdb2gmx has been used, then look at the right-hand comment column of the atom listing, which lists the cumulative charge. This should be an integer after every residue (and/or charge group where applicable). This will assist in finding the residue where things start departing from integer values. Also check the terminal capping groups that have been used.

If the charge is already close to an integer, then the difference is caused by rounding errors and not a major problem.

Note for PME users: It is possible to use a uniform neutralizing background charge in PME to compensate for a system with a net background charge. This may however, especially for non-homogeneous systems, lead to unwanted artifacts, as shown in 181 (http://pubs.acs.org/doi/abs/10.1021/ct400626b). Nevertheless, it is a standard practice to actually add counter-ions to make the system net neutral.

Incorrect number of parameters

Look at the topology file for the system. You’ve not given enough parameters for one of the bonded definitions. Sometimes this also occurs if you’ve mangled the Include File Mechanism or the topology file format (see: reference manual Chapter 5) when you edited the file.

Number of coordinates in coordinate file does not match topology

This is pointing out that, based on the information provided in the topology file, top, the total number of atoms or particles within the system does not match exactly with what is provided within the coordinate file, often a gro or a pdb.

The most common reason for this is simply that the user has failed to update the topology file after solvating or adding additional molecules to the system, or made a typographical error in the number of one of the molecules within the system. Ensure that the end of the topology file being used contains something like the following, that matches exactly with what is within the coordinate file being used, in terms of both numbers and order of the molecules:

[ molecules ]
; Compound   #mol
Protein      1
SOL          10189
NA+          10

Fatal error: No such moleculetype XXX

Each type of molecule in your [ molecules ] section of your top file must have a corresponding [ moleculetype ] section defined previously, either in the top file or an included itp file. See the reference manual section 5.6.1 for the syntax description. Your top file doesn’t have such a definition for the indicated molecule. Check the contents of the relevant files, how you have named your molecules, and how you have tried to refer to them later. Pay attention to the status of #ifdef and / or #include statements.

T-Coupling group XXX has fewer than 10% of the atoms

It is possible to specify separate thermostats (temperature coupling groups) for every molecule type within a simulation. This is a particularly bad practice employed by many new users to molecular dynamics simulations. Doing so is a bad idea, as you can introduce errors and artifacts that are hard to predict. In some cases it is best to have all molecules within a single group, using the default System group. If separate coupling groups are required to avoid the hot-solvent, cold-solute problem, then ensure that they are of sufficient size and combine molecule types that appear together within the simulation. For example, for a protein in water with counter-ions, one would likely want to use Protein and Non-Protein.

The cut-off length is longer than half the shortest box vector or longer than the smallest box diagonal element. Increase the box size or decrease rlist

This error is generated in the cases as noted within the message. The dimensions of the box are such that an atom will interact with itself (when using periodic boundary conditions), thus violating the minimum image convention. Such an event is totally unrealistic and will introduce some serious artefacts. The solution is again what is noted within the message, either increase the size of the simulation box so that it is at an absolute minimum twice the cut-off length in all three dimensions (take care here if are using pressure coupling, as the box dimensions will change over time and if they decrease even slightly, you will still be violating the minimum image convention) or decrease the cut-off length (depending on the force field utilised, this may not be an option).

Atom index (1) in bonds out of bounds

This kind of error looks like:

Fatal error:
[ file spc.itp, line 32 ]
Atom index (1) in bonds out of bounds (1-0).
This probably means that you have inserted topology
section "settles" in a part belonging to a different
molecule than you intended to. in that case move the
"settles" section to the right molecule.

This error is fairly self-explanatory. You should look at your top file and check that all of the [molecules] sections contain all of the data pertaining to that molecule, and no other data. That is, you cannot #include another molecule type (itp file) before the previous [moleculetype] has ended. Consult the examples in chapter 5 of the reference manual for information on the required ordering of the different [sections]. Pay attention to the contents of any files you have included with #include directives.

This error can also arise if you are using a water model that is not enabled for use with your chosen force field by default. For example, if you are attempting to use the SPC water model with an AMBER force field, you will see this error. The reason is that, in spc.itp, there is no #ifdef statement defining atom types for any of the AMBER force fields. You can either add this section yourself, or use a different water model.

XXX non-matching atom names

This error usually indicates that the order of the topology file does not match that of the coordinate file. When running grompp, the program reads through the topology, mapping the supplied parameters to the atoms in the coordinate file. If there is a mismatch, this error is generated. To remedy the problem, make sure that the contents of your [ molecules ] directive matches the exact order of the atoms in the coordinate file.

In a few cases, the error is harmless. Perhaps you are using a coordinate file that has the old (pre-4.5) ion nomenclature. In this case, allowing grompp to re-assign names is harmless. For just about any other situation, when this error comes up, it should not be ignored. Just because the -maxwarn option is available does not mean you should use it in the blind hope of your simulation working. It will undoubtedly blow up.

The sum of the two largest charge group radii (X) is larger than rlist - rvdw/rcoulomb

This error warns that some combination of settings will result in poor energy conservation at the longest cutoff, which occurs when charge groups move in or out of pair list range. The error can have two sources:

  • Your charge groups encompass too many atoms. Most charge groups should be less than 4 atoms or less.
  • Your mdp settings are incompatible with the chosen algorithms. For switch or shift functions, rlist must be larger than the longest cutoff (rvdw or rcoulomb) to provide buffer space for charge groups that move beyond the neighbor searching radius. If set incorrectly, you may miss interactions, contributing to poor energy conservation.

A similar error (“The sum of the two largest charge group radii (X) is larger than rlist”) can arise under two following circumstances:

  • The charge groups are inappropriately large or rlist is set too low.
  • Molecules are broken across periodic boundaries, which is not a problem in a periodic system. In this case, the sum of the two largest charge groups will correspond to a value of twice the box vector along which the molecule is broken.

Invalid line in coordinate file for atom X

This error arises if the format of the gro file is broken in some way. The most common explanation is that the second line in the gro file specifies an incorrect number of atoms, causing grompp to continue searching for atoms but finding box vectors.

Errors in mdrun

Stepsize too small, or no change in energy. Converged to machine precision, but not to the requested Fmax

This may not be an error as such. It is simply informing you that during the energy minimization process mdrun reached the limit possible to minimize the structure with your current parameters. It does not mean that the system has not been minimized fully, but in some situations that may be the case. If the system has a significant amount of water present, then an Epot of the orderof -105 to -106 (in conjunction with an Fmax between 10 and 1000 kJ mol-1 nm-1) is typically a reasonable value for starting most MD simulations from the resulting structure. The most important result is likely the value of Fmax, as it describes the slope of the potential energy surface, i.e. how far from an energy minimum your structure lies. Only for special purposes, such as normal mode analysis type of calculations, it may be necessary to minimize further. Further minimization may be achieved by using a different energy minimization method or by making use of double precision-enabled GROMACS.

Energy minimization has stopped because the force on at least one atom is not finite

This likely indicates that (at least) two atoms are too close in the input coordinates, and the forces exerted on each other are greater in magnitude than can be expressed to the extent of the precision of GROMACS, and therefore minimization cannot proceed. It is sometimes possible to minimize systems that have infinite forces with the use of soft-core potentials, which scale down the magnitude of Lennard-Jones interactions with the use of the GROMACS free energy code. This approach is an accepted workflow for equilibration of some coarse-grained systems such as Martini.

LINCS/SETTLE/SHAKE warnings

Sometimes, when running dynamics, mdrun may suddenly stop (perhaps after writing several pdb files) after a series of warnings about the constraint algorithms (e.g. LINCS, SETTLE or SHAKE) are written to the log file. These algorithms often used to constrain bond lengths and/or angles. When a system is blowing up (i.e. exploding due to diverging forces), the constraints are usually the first thing to fail. This doesn’t necessarily mean you need to troubleshoot the constraint algorithm. Usually it is a sign of something more fundamentally wrong (physically unrealistic) with your system. See also the advice here about diagnosing unstable systems.

1-4 interaction not within cut-off

Some of your atoms have moved so two atoms separated by three bonds are separated by more than the cut-off distance. This is BAD. Most importantly, do not increase your cut-off! This error actually indicates that the atoms have very large velocities, which usually means that (part of) your molecule(s) is (are) blowing up. If you are using LINCS for constraints, you probably also already got a number of LINCS warnings. When using SHAKE this will give rise to a SHAKE error, which halts your simulation before the 1-4 not within cutoff error can appear.

There can be a number of reasons for the large velocities in your system. If it happens at the beginning of the simulation, your system might be not equilibrated well enough (e.g. it contains some bad contacts). Try a(nother) round of energy minimization to fix this. Otherwise you might have a very high temperature, and/or a timestep that is too large. Experiment with these parameters until the error stops occurring. If this doesn’t help, check the validity of the parameters in your topology!

Simulation running but no output

Not an error as such, but mdrun appears to be chewing up CPU time but nothing is being written to the output files. There are a number of reasons why this may occur:

  • Your simulation might simply be (very) slow, and since output is buffered, it can take quite some time for output to appear in the respective files. If you are trying to fix some problems and you want to get output as fast as possible, you can set the environment variable GMX_LOG_BUFFER to 0.
  • Something might be going wrong in your simulation, causing e.g. not-a-numbers (NAN) to be generated (these are the result of e.g. division by zero). Subsequent calculations with NAN’s will generate floating point exceptions which slow everything down by orders of magnitude.
  • You might have all nst* parameters (see your mdp file) set to 0, this will suppress most output.
  • Your disk might be full. Eventually this will lead to mdrun crashing, but since output is buffered, it might take a while for mdrun to realize it can’t write.

Can not do Conjugate Gradients with constraints

This means you can’t do energy minimization with the conjugate gradient algorithm if your topology has constraints defined. Please check the reference manual.

Pressure scaling more than 1%

This error tends to be generated when the simulation box begins to oscillate (due to large pressures and / or small coupling constants), the system starts to resonate and then crashes. This can mean that the system isn’t equilibrated sufficiently before using pressure coupling. Therefore, better / more equilibration may fix the issue.

It is recommended to observe the system trajectory prior and during the crash. This may indicate if a particular part of the system / structure is the problem.

In some cases, if the system has been equilibrated sufficiently, this error can mean that the pressure coupling constant, tau-p, is too small (particularly when using the Berendsen weak coupling method). Increasing that value will slow down the response to pressure changes and may stop the resonance from occurring. You are also more likely to see this error if you use Parrinello-Rahman pressure coupling on a system that is not yet equilibrated - start with the much more forgiving Berendsen method first, then switch to other algorithms.

This error can also appear when using a timestep that is too large, e.g. 5 fs, in the absence of constraints and / or virtual sites.

Range Checking error

This usually means your simulation is blowing up. Probably you need to do better energy minimization and/or equilibration and/or topology design.

X particles communicated to PME node Y are more than a cell length out of the domain decomposition cell of their charge group

This is another way that mdrun tells you your system is blowing up. If you have particles that are flying across the system, you will get this fatal error. The message indicates that some piece of your system is tearing apart (hence out of the “cell of their charge group”). Refer to the Blowing Up page for advice on how to fix this issue.

A charge group moved too far between two domain decomposition steps.

See information above.

Software inconsistency error: Some interactions seem to be assigned multiple times

See information above

There is no domain decomposition for n ranks that is compatible with the given box and a minimum cell size of x nm

This means you tried to run a parallel calculation, and when mdrun tried to partition your simulation cell into chunks, it couldn’t. The minimum cell size is controlled by the size of the largest charge group or bonded interaction and the largest of rvdw, rlist and rcoulomb, some other effects of bond constraints, and a safety margin. Thus it is not possible to run a small simulation with large numbers of processors. So, if grompp warned you about a large charge group, pay attention and reconsider its size. mdrun prints a breakdown of how it computed this minimum size in the log file, so you can perhaps find a cause there.

If you didn’t think you were running a parallel calculation, be aware that from 4.5, GROMACS uses thread-based parallelism by default. To prevent this, give mdrun the -ntmpi 1 command line option. Otherwise, you might be using an MPI-enabled GROMACS and not be aware of the fact.