We present a molecular dynamics method for the simulation of reactions on metal surfaces. The method employs reactive potentials obtained from our modified version of the bond order conservation method of Shustorovich. Although we have used these potentials in solving classical equations of motion, quantum mechanical calculations could also employ them. Example calculations in four and six dimensions are given for the dissociation of H-2 over Ni(111) and H-SCH3 over Au(111). It is demonstrated that the calculated quantities representing reaction rate constants show the correct Arrhenius behavior. We show that the coupling of the reaction coordinate to "rotational" degrees of freedom contains important effects in agreement with existing work in the literature.