Solvation effect might have a tremendous influence on chemical reactions. However, precise quantum chemistry calculations are most often done either in vacuum neglecting the role of the solvent or using continuum solvent model ignoring its molecular nature. We propose a new method coupling a quantum description of the solute using electronic density functional theory with a classical grand-canonical treatment of the solvent using molecular density functional theory. Unlike previous work, both densities are minimised self consistently, accounting for mutual polarisation of the molecular solvent and the solute. The electrostatic interaction is accounted using the full electron density of the solute rather than fitted point charges. The introduced methodology represents a good compromise between the two main strategies to tackle solvation effect in quantum calculation. It is computationally more effective than a direct quantum-mechanics/molecular mechanics coupling, requiring the exploration of many solvent configurations. Compared to continuum methods it retains the full molecular-level description of the solvent. We validate this new framework onto two usual benchmark systems: a water solvated in water and the symmetrical nucleophilic substitution between chloromethane and chloride in water. The prediction for the free energy profiles are not yet fully quantitative compared to experimental data but the most important features are qualitatively recovered. The method provides a detailed molecular picture of the evolution of the solvent structure along the reaction pathway.