Weak interactions of graphene surface with reactive molecular impurities are the subject of many studies since noncovalent functionalization of surface via molecular doping is a powerful tool for tuning the electronic properties of graphene layers. In this work, the adsorption of diatomic halogen gas molecules, F, Cl, Br, I onto bilayer and multilayer pristine graphene surfaces were studied comparatively by Monte Carlo (MC) and molecular dynamics (MD) simulation techniques in canonical ensemble. The adsorption sites, adsorption capacity, coverage factors, adsorption isotherms, and adsorption kinetics were investigated and the adsorption energies were calculated for all adsorbates. Graphene was modeled as a two-dimensional layer of 200 carbon atoms in a honeycomb arrangement. The COMPASS force field was used in the simulations. The adsorption isotherms were obtained and fitted to Langmuir model. The kinetics of adsorption was studied and found to be first order. Both the monolayer and the multilayer adsorption of halogen molecules showed that van der Waals volumes of halogen molecules and also their polarizabilities display a competitive role in the saturation capacity and the strength of surface interactions.