Ionic solvation in a polar solvent is modeled by a central charge surrounded by dipolar molecules posted on two-dimensional distorted lattice sites with simple rotational dynamics. Density of states is calculated by applying the Wang-Landau algorithm to both the energy and polarization states. The free energy landscapes of solvent molecules as a function of polarization are depicted to explore the competition between the thermal fluctuation and solvation energy. Without a central charge, for temperatures higher than the energy scale of the dipole-dipole interactions, the energy landscape for the small polarization region exhibits a parabolic shape as predicted by Marcus [Rev. Mod. Phys. 65, 599 (1993)] for electron transfer reaction, while there is an additional quartic contribution to the landscape for the large polarization region. When the temperature drops, the simulated free energy landscapes are no longer smooth due to the presence of multiple local minima arising from the frustrated interaction among the dipoles. The parabolic contribution becomes negligible and the energy landscape becomes quartic in shape. For a strong central charge, the energy landscape exhibits an asymmetric profile due to the contributions of linear and cubic terms that arise from the charge-dipole interactions.