Achieving ultra-deep hydrodesulfurization means enabling removal of the last fractions of sulfur, contained in refractory molecules, from oil. Improving the state-of-the-art Co-promoted MoS₂ (CoMoS) catalyst or the development of novel catalysts is crucial for this. Improving CoMoS requires more insight in the way sulfur containing molecules interact with it. Herein, we model the adsorption of sulfur containing molecules on the S-edge, M-edge, corner and basal plane of CoMoS using density functional theory. The obtained adsorption configurations and energies point to a preference towards physisorption at the S-edge and chemisorption in vacancies at the M-edge and corner. Smaller molecules, such as thiophene and methylthiol, were found to prefer vacancies when adsorbing while larger, sterically hindered molecules as 4,6-dimethyldibenzothiophene prefer physisorption on the brim of the edges or the basal plane through van der Waals interactions. Hydrogenation generally leads to a preference towards adsorption at vacancies for thiophene and dibenzothiophene while for 4,6-dimethyldibenzothiophene hydrogenation leads to preferential adsorption on the S-edge brim, possibly explaining why 4,6-dimethyldibenzothiophene does not get desulfurized directly but follows a hydrogenation route. Thiolate formation energies were also calculated for the different molecules and used to predict which sites are most likely to be involved in breaking carbon-sulfur bonds. The thiolate formation energies show the inert nature of the basal plane towards breaking carbon-sulfur and sulfur-hydrogen bonds. Additionally, activation energies for thiophene and dibenzothiophene carbon-sulfur bond scission indicate that both molecules follow the direct desulfurization route on under-coordinated sites or vacancies.