Molybdenum (Mo) is a redox-sensitive metal element. Its isotope composition (delta 98Mo) has the potential to reconstruct oceanic oxygenation extent in deep time. To accurately reconstruct the bulk delta 98Mo value of ancient seawater from sediment records, it is necessary to clarify the Mo isotope fractionation mechanism during the transformation among thiomolybdates (MoOxS4-x2-, x = 0-4) in marine environments. For aqueous systems, predicting the equilibrium isotope fractionation factors theoretically needs careful treatment of solvent effects. The commonly used water-droplet or periodical boundary condition methods are not only troublesome but also time-consuming. We test the efficiency of the two-layer own N-layer integrated molecular orbital molecular mechanics (ONIOM) method in predicting the reduced partition function ratios (beta factors) of the MoOxS42- species in aqueous phases. Compared to other methods, the ONIOM method can produce beta factors with the same accuracy and considerably fewer computational resources. The predicted beta factors are in the order MoO4(aq)2- > MoO3S(aq)2- > MoO2S2(aq)2- > MoOS3(aq)2- > MoS4(aq)2-, which are in line with previous laboratory experiments and field observations. The new results provide the basic data for reconstructing delta 98Mo values in suboxic, ferruginous, and weakly euxinic sediments.