Three dimensional implementations of liquid state theories offer an efficient alternative to computer simulations for the atomic-level description of aqueous solutions in complex environments. In this context, we present a (classical) molecular density functional theory (MDFT) of water that is derived from first principles and is based on two classical density fields, a scalar one, the particle density, and a vectorial one, the multipolar polarization density. Its implementation requires as input the partial charge distribution of a water molecule and three measurable bulk properties, namely the structure factor and the k-dependent longitudinal and transverse dielectric constants. It has to be complemented by a solute-solvent three-body term that reinforces tetrahedral order at short range.The approach is shown to provide the correct three-dimensional microscopic solvation profile around various molecular solutes, possibly possessing H-bonding sites, at a computer cost two-three orders of magnitude lower than with explicit simulations. [22][23][24][25]. They can yield reliable predictions for both the microscopic structure and the thermodynamic properties of molecular fluids in bulk, interfacial, or confined conditions at a much more modest computational cost than molecular-dynamics or Monte-Carlo simulations. A current challenge concerns their implementation in three dimensions in order to describe molecular liquids, solutions, and mixtures in complex environments such as atomistically-resolved solid interfaces or biomolecular media. There have been a number of recent efforts in that direction using 3D-RISM [26][27][28][29][30][31], lattice field [32,33] or Gaussian field [20,34] theories. Recently, a molecular density functional theory (MDFT) approach to solvation has been introduced. [35][36][37][38][39][40] It relies on the definition of a free-energy functional depending on the full six-dimensional position and orientation solvent density. In the so-called homogeneous reference fluid (HRF) approximation, the (unknown) excess free energy can be inferred from the angular-dependent direct correlation function of the bulk solvent, that can be predetermined from molecular simulations of the pure solvent. Compared to reference molecular dynamics calculations, such approximation was shown to be accurate for polar, non-hydrogen bonded fluids [35,[38][39][40][41], but to require some corrections for water [39,[42][43][44].In this paper we introduce a simplified version of MDFT that can be derived rigorously for SPC-or TIPnP-like representations of water, involving a single LennardJones interaction site and distributed partial charges. In that case we will seek a simpler functional form expressed in terms of the particle density n(r) and site-distributed polarisation density P(r), and requiring as input simpler physical quantities than the full position and orientation-dependent direct correlation function. This functional may be considered as a multipolar generalization of the generic dipolar fluid free-ene...