Experiments have suggested that ion correlation and fluctuation effects can be potentially important for multivalent ions in RNA folding. However, most existing computational methods for the ion electrostatics in RNA folding tend to ignore these effects. The previously reported Tightly Bound Ion (TBI) model can treat ion correlation and fluctuation but its applicability to biologically important RNAs is severely limited by the low computational efficiency. Here, based on Monte Carlo sampling for the many-body ion distribution, we develop a new computational model, Monte Carlo Tightly Bound Ion (MCTBI) model, for ion binding properties around an RNA. Due to an enhanced sampling algorithm for ion distribution, the model leads to significant improvement in computational efficiency. For example, for a 160-nt RNA, the model causes more than 10-fold increase in the computational efficiency, and the improvement in computational efficiency is more pronounced for larger systems. Furthermore, unlike the earlier model, which describes ion distribution using the number of bound ions around each nucleotide, the current MCTBI model is based on the three-dimensional coordinates of the ions. The higher efficiency of the model allows us to treat the ion effects for medium to large RNA molecules, RNA-ligand complexes, and RNA-protein complexes. This new model, together with proper RNA conformational sampling and energetics model, may serve as a starting point for further development for the ion effects in RNA folding and conformational changes and for large nucleic acids systems.