An efficient discontinuous molecular dynamics simulation algorithm is developed to simulate hydrogen-bonding systems. The efficiency of CPU time in the order of N 1092 N is stressed by comparisons to several other methods and algorithms. To evaluate this algorithm, special attention is given to the simulation of thermodynamic and transport properties as well as oligomer distribution of the hydrogen-bonding water system by use of DMD-1B and -4B site models. The results are compared to the results of Monte Carlo simulations and Wertheim’s first-order theory. The overall agreement is satisfactory within a standard deviation.