  Energy density and energy flux was introduced along Takesue's method.
Particle energies were localized at particle positions using Dirac delta
function. The energy density was connected with the energy flux by continuity
equation. New method was proposed to deal with multi--body interaction.
Consequently, the energy current between particles could be calculated even if
the multi--body interaction appeared. Moreover, An application to the molecular
dynamics simulation of hydrogen adsorption on a graphene with modified Brenner
reactive empirical bond order potential was demonstrated.
