The problem may be related to time step.

The time step used in the source code can not ensure the conservation of energy , necessitating smaller time step .

Another possible reason for the absence of energy conservation is "bad" numerical integration scheme. For example, the 4th order Runge-Kutta method does not conserve energy. The simplest symplectic (energy conserving) method is leap-frog.

I would also suggest to start with small (or even zero) temperatures. It allows to avoid other possible issues like cut-off radius, etc.

Good luck,

Vitaly

The typical problem, when implementing MD simulations with any non-trivial interaction law, is to be sure that the expression for forces  is consistent with the formula for potential energy. With a non-pairwise interaction like Tersoff's potential, the first thing I would check is the above mentioned consistency. The simplest way  to do it, and that's the way I routinely use when I implement a new force field, is to evaluate numerically the gradient of the energy expression (for a minimal number of atoms)  and then comparing it to the explicit calculation  of forces.

Some help from algebraic manipulation software like Mathematica, Maple or Maxima may help as well.

Another possible source of trouble, may be   the interplay between the distance cut-off and periodic boundary conditions. You may try to see if things go better if your work without periodi boundary conditions  (like the case of a small cluster).