The prediction of the structural and thermodynamic properties of electrolyte solutions is critical for a huge range of practical situations where these solutions play a vital role. Theoretical models, such as the continuum solvent model, attempt to explain the behavior of solutions using a coarse-grained description of the interactions of species in the solution, whereas molecular simulations aim to directly compute the behavior of the solution, including the interactions between all ions and molecules in the system. Both methods have limitations: theoretical models are generally less accurate because they rely on assumptions, while molecular simulations require significant computational resources, particularly if higher accuracy is desired. To address these issues, we propose an affordable and effective method that combines the advantages of the modified Poisson–Boltzmann equation (MPBE) with classical molecular dynamics (MD) simulations to predict the radial distribution functions and thermodynamic properties of electrolyte solutions. We demonstrate a method of using the MPBE to compute the short-range potential of mean force (PMF) from the radial distribution functions (RDFs) and vice versa. Furthermore, we provide insights into the relationship between the RDFs and the short-range PMF based on the MPBE. Our analysis reveals that the effective short-range PMFs can be approximately calculated using low concentration simulations but the short-range PMFs are slightly concentration-dependent in simulations at higher concentrations. Additionally, we demonstrate that for concentrated solutions, osmotic coefficients can be calculated in agreement with experiment using a virial approach. This is based on the effective short-range PMFs and RDFs obtained from the MPBE method. Our proposed MPBE can therefore accelerate the calculation of the structural and thermodynamic properties of electrolyte solutions.