Different dynamical Monte Carlo algorithms to investigate molecule formation on surfaces are developed, evaluated and compared with the deterministic approach based on reaction-rate equations. These include a null event algorithm, the n-fold way / BKL algorithm and an "hybrid" variant of the latter. NO 2 formation by NO oxidation on Pyrex and O recombination on silica with the formation of O 2 are taken as case studies. It is shown that dynamical Monte Carlo schemes are flexible, simple to implement, describe easily elementary processes that are not straightforward to include in deterministic simulations, can run very efficiently if appropriately chosen and give highly reliable results. Moreover, the present approach provides a relatively simple procedure to describe fully coupled surface and gas phase chemistries. The influence of the grid size on the CPU calculation time and the accuracy of the results, the effect of back diffusion and its inclusion in a deterministic formulation, and the influence of Langmuir-Hinsehlwood recombination involving two physisorbed atoms, are investigated and discussed.