Decomposition of Interaction Energies in LAMMPS

The Anatomy of a Simulation "Blow-Up"

When a simulation becomes unstable, it's typically because two or more atoms have moved too close to one another, resulting in an unphysically large repulsive force and a correspondingly massive potential energy. The total energy of the system spikes, integration algorithms fail, and the simulation crashes.

Your LAMMPS output might show the total energy (TotEng) rapidly increasing, but this single number doesn't tell you which interaction is responsible. The total potential energy, $E_{\text{pot}}$, is a sum of several terms:

Our first task is to see which of these components is the source of the runaway energy. You can customize your thermodynamic output in LAMMPS to print these individual components at each step. By adding a line like this to your input script, you can monitor the breakdown:

thermo_style custom step temp epair ecoul elong ebond eangle edihed press vol

In your system, which involves rigid SPC/E water molecules, the ebond, eangle, and edihed terms within water will be handled by constraints (e.g., SHAKE) and shouldn't be the primary source of instability between the nanoparticle and the water. The instability is most likely coming from the pairwise interactions, captured by epair (short-range) and elong (long-range, or kspace). If you see epair or elong skyrocket just before the crash, you've confirmed the problem lies in the non-bonded interactions.

Pinpointing the Problem Atoms

Knowing that the pair energy is the problem is a good start, but it doesn't tell you which atoms are involved. To find the specific atoms with extremely high energy, we use the compute pe/atom command. This command calculates the potential energy contribution for each atom individually.

The LAMMPS documentation provides the details for this command. Pay close attention to the optional keywords which allow you to calculate the per-atom energy arising from only the pair interactions.

By using this compute, you can tell LAMMPS to calculate the pair potential energy for every atom and then write this information to a dump file. Atoms with huge potential energies are your culprits.

Here is a practical example of how you would implement this in your LAMMPS script:

# Calculate the per-atom potential energy from pair interactions only
compute per_atom_pe all pe/atom pair




# Dump atom ID, type, position, and the calculated per-atom PE
# This dump should be frequent enough to capture the state just before a crash
dump 1 all custom 100 dump.instability_check id type x y z c_per_atom_pe

After running the simulation, you can inspect the dump.instability_check file. By sorting this file by the last column (c_per_atom_pe), you will immediately find the atoms with the highest potential energy. These are the atoms involved in the unphysically close contact.

The image below shows what these potential energy curves look like. When your simulation blows up, it means at least one pair of atoms has moved to a very small separation distance $r$, climbing to an extremely high point on the repulsive part of the curve. The compute pe/atom command helps you find which atoms are doing the climbing.

Isolating Problematic Group Interactions

Now you know which atoms have high energy. Given your system, it's highly likely they are at the interface between the maghemite nanoparticle and the surrounding water. The final diagnostic step is to confirm this by calculating the interaction energy specifically between these two groups of atoms. For this, LAMMPS provides the compute group/group command.

This command is ideal for quantifying the interaction between subsets of your system, such as a nanoparticle and the solvent.

To use this, you first need to define two groups: one for the atoms in your maghemite nanoparticle and one for the atoms in the water molecules. You can do this based on atom type. For example, if your maghemite Fe is type 1, maghemite O is type 2, water O is type 3, and water H is type 4, your script would include:

# Define groups based on atom types
group nanoparticle type 1 2
group water      type 3 4




# Compute the interaction energy between the two groups
# We include both short-range (pair) and long-range (kspace) interactions
compute NP_H2O_E nanoparticle group/group water pair yes kspace yes




# You can then output this value using thermo_style or a variable
thermo_style custom step temp c_NP_H2O_E ...

A large, positive, and rapidly increasing value for c_NP_H2O_E would be definitive proof that the instability originates from the cross-interactions between your nanoparticle and the water.

Summary and Next Steps

In this lesson, we have established a three-step protocol for diagnosing numerical instability in your LAMMPS simulation:

1. System-Level Diagnosis: Use thermo_style custom to monitor the different components of the total energy (epair, elong, etc.) and identify which type of interaction is becoming unstable.
2. Atom-Level Diagnosis: Use compute pe/atom combined with dump custom to identify the specific atoms that are acquiring massive potential energies just before the simulation crashes.
3. Group-Level Diagnosis: Use group and compute group/group to quantify the interaction energy between your nanoparticle and the water, confirming that the instability lies at the interface.

By following this workflow, you can move beyond simply knowing that your simulation failed and instead gather concrete data that points to the exact atoms and interaction types at the heart of the problem.

In our next lesson, we will take the atoms you've identified with this method and use visualization and further analysis to see why their interaction is so problematic. We will look for atomic overlaps at the water-oxide interface, which is the most likely physical reason for the high energies you are diagnosing.