special_bonds command in combination with the lj/cut/coul/long pair potential

I am not sure what could be wrong. But many things could be wrong. It may be a problem with your input.

  1. Check whether all the force parameters are right.
  2. When you say you are relaxing the system, do you have abnormally large bonds ?. Try visualizing your data file.
  3. What is your initial potential energy of the system ?. Do you have overlaps ?. Are you using periodic boundary conditions ?. Are you minimizing the system prior to equilibration ?
  4. When you say simulation runs fine for how long ?. what is the time step ?. How fast is the simulation box shrinking ? How are you assigning velocities to the atom ?

If the problem persists after checking the above, please post the output from lammps.