Herein, a new theoretical method is presented for predicting osmotic equilibria and activities, where a bulk liquid and its corresponding vapor phase are simulated by means of molecular dynamics using explicit polarization. Calculated time averaged number density profiles provide the amount of evaporated molecules present in the vapor phase and consequently the vapor-phase density. The activity of the solvent and the corresponding osmotic coefficient are determined by the vapor density at different solute concentrations with respect to the reference vapor, density of the pure solvent. With the extended Debye Hiickel equation for the activity coefficient along with the corresponding Gibbs Duhem relation, the activity coefficients of the solutes are calculated by fitting the osmotic coefficients. A simple model based on the combination of Poisson processes and Maxwell Boltzmann velocity distributions is introduced to interpret statistical phenomena observed during the simulations, which are related to evaporation and recondensation. This method is applied to aqueous dysprosium nitrate [Dy(NO3)(3)] solutions at different concentrations. The obtained densities of the liquid.bulk and the osmotic and activity coefficients are in good agreement with the experimental results for concentrated and saturated solutions. Density profiles of the liquid vapor interface at different concentrations provide detailed insight into the spatial distributions of all compounds.