A new approach to molecular diffusion is developed using density functionals for fluxes and the Metropolisalgorithm in the mass balance equation. This procedure results in a new equation for diffusion of interactingparticles which has multiple solutions and gives density distributions for coexisting and metastable phases. Itis shown that the diffusion of interacting molecules is driven by two variables: the density gradient of moleculesand the density of molecule−vacancy pairs (“pseudoparticles”).