E_pair discontinuity after restart and simulation crash

Dear lammps users,
I’m currently using the stable LAMMPS version from August 29th, 2024, to run simulations with a custom fix designed to create and delete bonds to form loops. My system consists of 50 bead-spring polymers, each containing 500 beads connected by FENE bonds. Additionally, “patchy” beads are added 0.5 units away from each bead center, with each bead and its corresponding patchy bead forming a rigid body. Regular beads interact via classic LJ repulsion, while patchy beads interact through an attractive Morse potential.
In several of these runs, using an NVE ensemble with a Langevin thermostat, I observe that after around 5e8 timesteps (with a timestep of 0.01 Brownian times), the E_pair stabilizes around -2. However, if I restart from a LAMMPS restart file saved at this point, the simulation crashes immediately, showing an E_pair around 5e+26, which is highly unexpected.

These are the pair potential currently used (and not changed after the restart):

pair_coeff * * lj/cut 0.0 0.0 0.0
pair_coeff 1 1 lj/cut 1.0 1.0 1.12246
pair_coeff 1 3 lj/cut 1.0 1.0 1.12246
pair_coeff 3 3 lj/cut 1.0 1.0 1.12246
pair_coeff 4 2 morse 25 14.0 0.0 0.20
pair_coeff 5 2 morse 25 14.0 0.0 0.20

Additionally, this happen either with or without my custom fix running from this restart point and even if the restart file is transformed into a data file. Do you have any insight on what could be causing this issue?

Thanks for your support,
Filippo

This is impossible to say with so little context information.

Best would be, if you can create a small test system that can be run quickly and reproduces the same behavior without required a long wait on a laptop or a local workstation. Then a developer can take this input deck and experiment with it and apply various modifications to determine what is at fault. There are many possibilities, and without even seeing the input, there is not chance of even speculating about it.