This paper presents the nuclear-electronic orbital density functional theory [NEO-DFT(ee)] method for including electron-electron correlation and nuclear quantum effects self-consistently in quantum chemical calculations. The NEO approach is designed to treat a relatively small number of nuclei quantum mechanically, while the remaining nuclei are treated classically. In the NEO-DFT(ee) approach, the correlated electron density is used to obtain the nuclear molecular orbitals, and the resulting nuclear density is used to obtain the correlated electron density during an iterative procedure that continues until convergence of both the nuclear and electronic densities. This approach includes feedback between the correlated electron density and the nuclear wavefunction. The application of this approach to bihalides and acetylene indicates that the nuclear quantum effects do not significantly impact the electron correlation energy, but the quantum nuclear energy is enhanced in the NEO-DFT(ee) B3LYP method. The excellent agreement of the NEO-DFT(ee)-optimized bihalide structures with the vibrationally averaged geometries from grid-based quantum dynamical methods provides validation for the NEO-DFT(ee) approach. Electron-proton correlation could be included by the development of an electron-nucleus correlation functional. Alternatively, explicit electron-proton correlation could be included directly into the NEO self-consistent-field framework with Gaussian-type geminal functions.