Ctrl+K

Debugging simulation crashes#

One of the most common crashes you will encounter using ACEMD is

openmm.OpenMMException: Particle coordinate is NaN.

This crash can be caused by various issues in the setup of your system or the ACEMD input file.

Possible causes#

  1. Clashes of atoms in the starting coordinates

  2. Incorrect parameterization of small molecules

  3. Incorrect box dimensions for a solvated system

  4. Overly strong external forces (like restraints)

  5. NNP instabilities

General advice for debbugging issues#

If your simulation crashes directly after minimization you can try increasing the minimization steps. If atom clashes are the cause of the crash, a larger number of minimization steps might give the system the chance to converge to a lower energy solution where the clashes are resolved.

If the simulation still crashes despite longer minimization, or during a production run, you should try setting the following options in your input.yaml file:

This will make ACEMD write out the starting coordinates in the XTC file (stepzero: true) and will make ACEMD write into the trajectory file at every single simulation step (trajectoryperiod: 1).

Warning

These options will significantly decrease simulation speed, as the writing of the trajectory files at every step is orders of magnitude slower than the MD step time. After fixing your issue, remove these settings from the input file.

With these settings you should be able to see the steps right before the crash of the simulation. If you visualize the simulation using pymol or vmd or any other MD visualization tool, you should be able to see some unnatural bond stretchings or an atom jumping very rapidly across the simulation area. You should take a closer look at the atoms or bonds which behave abnormally and see if there are any clashes with other atoms around them.

Note

Don’t be confused by atoms “leaving” the initial simulation box. Simulations produced by ACEMD are written as unwrapped trajectories, despite the periodic boundary conditions being used. See our tutorial on wrapping if you are confused about this https://software.acellera.com/htmd/tutorials/WrappingTutorial.html

Finally if the cause of the crash is still not apparent to the eye you can try also setting the following option in the input.yaml file:

trajforcefile: "forces.xtc"

This will make ACEMD write out the forces applied to each atom at every trajectoryperiod steps. The XTC file format is normally used to store atom coordinates but it can be abused for storing forces as well. One way of visualizing the forces is using the view_forces function from the acemd module. This function requires VMD to be installed on your system and available in the system PATH.

from acemd.utils import view_forces
from moleculekit.molecule import Molecule

mol = Molecule("structure.prmtop")
mol.read("output.xtc")
view_forces(mol, "forces.xtc", step=0)

This will visualize the first frame of the simulation and the forces which were applied on each atom, to see other simulation steps set the step= value accordingly.

If you wait a little bit for VMD to stop drawing all the arrows, you will see the force vectors as green arrows originating from each atom. The larger the force, the larger the arrow. If you find the atoms with the larger forces and look in their environment you should be able to figure out the cause of the explosion.

Starting coordinate clashes#

System building is a delicate procedure and many builders like CHARMM psfgen or AMBER tleap can produce structures with clashes on atoms if someone doesn’t pay too close attention. The above tips you help you identify the cause of the crash.

Incorrect parameterization of small molecules#

This mistake should manifest itself as extreme forces applied mostly to the parameterized small molecules (meaning residues which are not part of the chosen forcefield you are using but were provided additionally during the building procedure). If you follow the above tips and the strong forces seems to step from the small molecule, one possible cause could be that your molecule has been wrongly parameterized and you should take a look at the parameters which were produced.

Incorrect box dimensions#

If you are simulating a solvated system, the box dimensions should be set correctly in the input file. If you are running an equilibration simulation (or any simulation with barostat: true), the box dimensions will adapt during the simulation to the target pressure. However if the initial box dimensions are very far from the initial box size in the starting coordinates, the barostat might not be able to save the simulation from crashing.

These mistakes can usually be detected by looking at the simulation logs and monitoring value of the box volume. If it has very large fluctuations it could be that the cause is wrong starting dimensions. Some fluctuation is to be expected though if the barostat is on.

If you are running production simulations (or any simulation with barostat: false), the box dimensions will not adapt anymore and therefore should be set to the exact correct box dimensions. Therefore we also generally only perform production runs after the equilibration runs and use the final box size of the equilibrarion as the box size of the production runs.

Strong external forces#

Another cause for simulation crashes are very strong external forces applied to the atoms. By external forces we mean positional restraints as well as NNP forces. You should take a look at your input.yaml file and the extforces section to see what external forces are being applied. If you are applying positional restraints, try removing them from the file and running a short simulation. If the simulation works, this means that the restraints are the cause of the crash. This could mean that either they are too strong (check the k values) or the restraints are applied on the wrong atoms (check the sel value). Fixing the atom selection or decreasing the force constants should resolve the issue.

NNP instabilities#

Neural network potentials have many advantages over classical MD simulations but sometimes they can produce poor extrapolations on systems which are far from their training data. In those cases the forces which they predict for your system might be wrong and lead to crashes. One way to increase their stability is by reducing the timestep of the simulation in the input.yaml file. This will possibly lead to fewer extrapolations but will accordingly decrease your simulation speed. Another way to solve this issue is to use a different NNP model which is better suited for your system. Our AceFF NNP models are continuously being developed and improved so take a look if a newer version of the AceFF forcefield is available.

On this page