  $H_2$ is the most abundant interstellar species. Its deuterated forms ($HD$
and $D_2$) are also significantly abundant. Huge abundances of these molecules
could be explained by considering the chemistry occurring on the interstellar
dust. Because of its simplicity, Rate equation method is widely used to study
the formation of grain-surface species. However, since recombination efficiency
of formation of any surface species are heavily dependent on various physical
and chemical parameters, Monte Carlo method would be best method suited to take
care of randomness of the processes. We perform Monte Carlo simulation to study
the formation of $H_2$, $HD$ and $D_2$ on interstellar ices. Adsorption
energies of surface species are the key inputs for the formation of any species
on interstellar dusts but binding energies of deuterated species are yet to
known with certainty. A zero point energy correction exists between
hydrogenated and deuterated species which should be considered while modeling
the chemistry on the interstellar dusts. Following some earlier studies, we
consider various sets of adsorption energies to study the formation of these
species in diverse physical circumstances. As expected, noticeable difference
in these two approaches (Rate equation method and Monte Carlo method) is
observed for production of these simple molecules on interstellar ices. We
introduce two factors, namely, $S_f$ and $\beta$ to explain these
discrepancies: $S_f$ is a scaling factor, which could be used to correlate
discrepancies between Rate equation and Monte Carlo methods. $\beta$ factor
indicates the formation efficiency under various circumstances. Higher values
of $\beta$ indicates a lower production efficiency. We found that $\beta$
increases with a decrease in rate of accretion from gas phase to grain phase.
