In this work, a method to estimate solid elemental sulfur solubility in pure and gas mixtures using Monte Carlo (MC) molecular simulation is proposed. This method is based on Isobaric-Isothermal (NPT) ensemble and the Widom insertion technique for the gas phase and a continuum model for the solid phase. This method avoids the difficulty of having to deal with high rejection rates that are usually encountered when simulating using Gibbs ensemble. The application of this method is tested with a system made of pure hydrogen sulfide gas (H2S) and solid elemental sulfur. However, this technique may be used for other solid-vapor systems provided the fugacity of the solid phase is known (e.g., through experimental work). Given solid fugacity at the desired pressure and temperature, the mole fraction of the solid dissolved in gas that would be in chemical equilibrium with the solid phase might be obtained. In other words a set of MC molecular simulation experiments is conducted on a single box given the pressure and temperature and for different mole fractions of the solute. The fugacity of the gas mixture is determined using the Widom insertion method and is compared with that predetermined for the solid phase until one finds the mole fraction which achieves the required fugacity. In this work, several examples of MC have been conducted and compared with experimental data. The Lennard-Jones parameters related to the sulfur molecule model (ɛ, σ) have been optimized to achieve better match with the experimental work.