Multicomponent Density Functional Theory (MDFT) is a promising methodology to incorporate nuclear quantum effects, such as zero-point energy or tunneling, or to simulate other types of particles such as muons or positrons using particle densities as basic quantities. As for standard electronic DFT, a still ongoing challenge is to achieve the most efficient implementations. We introduce a multicomponent DFT implementation within the framework of auxiliary DFT, focusing on molecular systems comprising electrons and quantum protons. We introduce a dual variational procedure to determine auxiliary electron and proton densities which leads to a succession of approximate energy expressions. Electronic and protonic fitted densities are employed in (i) electron-electron, proton-proton, and electron-proton classical Coulomb interactions and (ii) electron exchange-correlation, proton-proton exchange, and electron-proton correlation (EPC) potentials. If needed, exact exchange among electrons or among protons is computed by the variational fitting of the corresponding Fock potential. The implementation is carried out in deMon2k. We test various electron proton correlation functionals on proton affinities. We find that auxiliary densities can be safely used in electron-electron, proton-proton, and electron-proton classical Coulomb interactions, as well as in EPC, albeit with some precautions related to the choice of the electronic auxiliary basis set that must be flexible enough. Computational tests reported indicate that introduction of density fitting in MDFT is clearly advantageous in terms of computational effort with good scaling properties with respect to the number of electrons and protons treated at the DFT level.