A novel grand canonical ensemble Monte Carlo (Pressure Controlled Grand Canonical ensemble Monte Carlo (PC-GCMC) method was developed for efficient estimation of solubility of molecules in polymer based on interatomic potential. This method contains novel algorithm, namely, the mu VT-CCB method, to simulate two different dynamics simultaneously, the insertion of a penetrant molecule into the polymer and the relaxation of the polymer chain. Therefore, this method is applicable to the estimation of solubility accompanied with polymer swelling or polymer plasticity. This method was validated by applying it to the estimation of solubility of various penetrants in polyethylene, and the results were compared with the experimental values. Calculation results indicated that this method could reproduce Henry and the Flory type isotherms, and obtained isotherms and swelling ratio were in good agreement with the reported experimental values. The required GCMC steps for calculation of solubility by PC-GCMC was ten times smaller than that by a conventional GCMC without the mu VT-CCB method. These results demonstrate the validity of PC-GCMC for efficient estimation of solubility in polymer.