Gromacs

Blowing Up

    Table of contents
    1. 1. Diagnosing an Unstable System

    Version as of 14:17, 30 Nov 2021

    to this version.

    Return to Version archive.

    View current version

    "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. 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,
    • you are using inappropriate temperature coupling, perhaps on inappropriate groups, or
    • you have a single water molecule somewhere within the system that is isolated from the other water molecules.

    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:
      • 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.
    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.

    Page last modified 12:34, 19 Jul 2012 by mabraham