Gromacs

Blowing Up

    Table of contents
    1. 1. Diagnosing an Unstable System

    "Blowing up" is a highly technical term used to describe a common sort of simulation failure. In brief, it describes a failure typically due to an unacceptably large force that ends up resulting in a failure of the integrator.

    To give a bit more background, it's important to remember that molecular dynamics numerically integrates Newton's equations of motion by taking small, discrete timesteps, and using these timesteps to determine new velocities and positions from velocities, positions, and forces at the previous timestep. If forces become too large at one timestep, this can result in extremely large changes in velocity/position when going to the next timestep. Typically, this will result in a cascade of errors: one atom experiences a very large force one timestep, and thus goes shooting across the system in an uncontrolled way in the next timestep, overshooting its preferred location or landing on top of another atom or something similar. This then results in even larger forces the next timestep, more uncontrolled motions, and so on. Ultimately, this will cause the simulation package to crash in some way, since it can't cope with such situations. In simulations with constraints, the first symptom of this will usually be some LINCS or SHAKE warning or error -- not because the constraints are the source of the problem, but just because they're the first thing to crash. Similarly, in simulations with domain decomposition, you may see messages about particles being more than a cell length out of the domain decomposition cell of their charge group, which are symptomatic of your underlying problem, and not the DD algorithm itself. Likewise for warnings about tabulated or 1-4 interactions being outside the distance supported by the table. This can happen on one computer system while another resulted in a stable simulation because of the impossibility of numerical reproducibility of these calculations on different computer systems.

    Possible causes include:

    • you didn't minimize well enough,
    • you have a bad starting structure, perhaps with steric clashes,
    • you are using too large a timestep (particularly given your choice of constraints),
    • you are doing particle insertion in free energy calculations without using soft core,
    • you are using inappropriate pressure coupling (e.g. when you are not in equilibrium, Berendsen can be best while relaxing the volume, but you will need to switch to a more accurate pressure-coupling algorithm later),
    • you are using inappropriate temperature coupling, perhaps on inappropriate groups, or
    • your position restraints are to coordinates too different from those present in the system, or
    • you have a single water molecule somewhere within the system that is isolated from the other water molecules, or
    • you are experiencing a bug in mdrun.

    Because blowing up is due, typically, to forces that are too large for a particular timestep size, there are a couple of basic solutions:

    1. Make sure the forces don't get that large, or
    2. use a smaller timestep.

    Better system preparation can help with 1, if the problems are occurring near the beginning of a simulation.

     

    Diagnosing an Unstable System

    Troubleshooting a system that is blowing up can be challenging, especially for an inexperienced user.  Here are a few general tips that one may find useful when addressing such a scenario:

    1. If the crash is happening relatively early (within a few steps), set nstxout (or nstxtcout) to 1, capturing all possible frames.  Watch the resulting trajectory to see which atoms/residues/molecules become unstable first.
    2. Simplify the problem to try to establish a cause:
      • If you have a new box of solvent, try minimizing and simulating a single molecule to see if the instability is due to some inherent problem with the molecule's topology or if instead there are clashes in your starting configuration.
      • If you have a protein-ligand system, try simulating the protein alone in the desired solvent.  If it is stable, simulate the ligand in vacuo to see if its topology gives stable configurations, energies, etc.
      • Remove the use of fancy algorithms, particularly if you haven't equilibrated thoroughly first
    3. Monitor various components of the system's energy using g_energy.  If an intramolecular term is spiking, that may indicate improper bonded parameters, for example.
    4. Make sure you haven't been ignoring error messages (missing atoms when running pdb2gmx, mismatching names when running grompp, etc) or using work-arounds (like using grompp -maxwarn when you shouldn't be) to make sure your topology is intact and being interpreted correctly.
    5. Make sure you are using appropriate settings in your .mdp file for the force field you have chosen and the type of system you have.  Particularly important settings are treatment of cutoffs, proper neighbor searching interval (nstlist), and temperature coupling.  Improper settings can lead to a breakdown in the model physics, even if the starting configuration of the system is reasonable.

    If using implicit solvation, starting your equilibration with a smaller time step than your production run can help energy equipartition more stably.

    There are several common situations in which instability frequently arises, usually in the introduction of new species (ligands or other molecules) into the system.  To determine the source of the problem, simplify the system (e.g. the case of a protein-ligand complex):

    1. Does the protein (in water) minimize adequately by itself? This is a test of the integrity of the coordinates and system preparation.  If this fails, something probably went wrong when running pdb2gmx (see below), or maybe genion placed an ion very close to the protein (it is random, after all).
    2. Does the ligand minimize in vacuo? This is a test of the topology.  If it does not, check your parameterization of the ligand and any implementation of new parameters in force field files.
    3. (If 2 is successful) Does the ligand minimize in water, and/or does a short simulation of the ligand in water succeed?

    Other sources of possible problems are in the biomolecule topology itself.

    1. Did you use -missing when running pdb2gmx?  If so, don't.  Reconstruct missing coordinates rather than ignoring them.
    2. Did you override long/short bond warnings by changing the lengths?  If so, don't.  You probably have missing atoms or some terrible input geometry.
    Page last modified 20:04, 13 Dec 2014 by JLemkul?