NEO density functional theory (NEO-DFT)^{Pak:2007, Chakraborty:2008, Chakraborty:2009} is an extension of DFT to multicomponent systems within the NEO framework. The Hohenberg-Kohn theorems have been extended to multicomponent systems, where the reference is expressed as the product of electronic and nuclear Slater determinants composed of Kohn-Sham orbitals. The NEO-DFT total energy is

$E[{\rho}^{\text{e}},{\rho}^{\text{p}}]={E}_{\text{ext}}[{\rho}^{\text{e}},{\rho}^{\text{p}}]+{E}_{\text{ref}}[{\rho}^{\text{e}},{\rho}^{\text{p}}]+{E}_{\text{exc}}[{\rho}^{\text{e}}]+{E}_{\text{pxc}}[{\rho}^{\text{p}}]+{E}_{\text{epc}}[{\rho}^{\text{e}},{\rho}^{\text{p}}]$ | (13.43) |

In this equation, ${E}_{\text{ext}}[{\rho}^{\text{e}},{\rho}^{\text{p}}]$ is the interaction of the electronic and protonic densities with the external potential created by the classical nuclei, and ${E}_{\text{ref}}[{\rho}^{\text{e}},{\rho}^{\text{p}}]$ contains the electron-electron, proton-proton, and electron-proton classical Coulomb energies, as well as the noninteracting kinetic energies of the quantum particles. The terms ${E}_{\text{exc}}[{\rho}^{\text{e}}]$, ${E}_{\text{pxc}}[{\rho}^{\text{p}}]$, and ${E}_{\text{epc}}[{\rho}^{\text{e}},{\rho}^{\text{p}}]$ are the electron-electron exchange-correlation functional, the proton-proton exchange-correlation functional, and the electron-proton correlation functional, respectively. The quantities ${\rho}^{\text{e}}({\mathbf{r}}_{1}^{\text{e}})=2{\sum}_{i=1}^{{N}_{\text{e}}/2}{|{\psi}_{i}^{\text{e}}({\mathbf{r}}_{1}^{\text{e}})|}^{2}$ and ${\rho}^{\text{p}}({\mathbf{r}}_{1}^{\text{p}})={\sum}_{I=1}^{{N}_{\text{p}}}{|{\psi}_{I}^{\text{p}}({\mathbf{r}}_{1}^{\text{p}})|}^{2}$ are the electron and proton densities, respectively, and ${\psi}_{i}^{\text{e}}({\mathbf{r}}_{1}^{\text{e}})$ and ${\psi}_{I}^{\text{p}}({\mathbf{r}}_{1}^{\text{p}})$ are the electronic and protonic Kohn-Sham spatial orbitals, respectively. These orbitals are obtained by solving two sets of coupled Kohn-Sham equations for the electrons and quantum protons:

$\left(-{\displaystyle \frac{1}{2}}{\nabla}^{2}+{v}_{\text{eff}}^{\text{e}}({\mathbf{r}}_{1}^{\text{e}})\right){\psi}_{i}^{\text{e}}={\u03f5}_{i}^{\text{e}}{\psi}_{i}^{\text{e}}$ | (13.44) | ||

$\left(-{\displaystyle \frac{1}{2{m}_{\text{p}}}}{\nabla}^{2}+{v}_{\text{eff}}^{\text{p}}({\mathbf{r}}_{1}^{\text{p}})\right){\psi}_{I}^{\text{p}}={\u03f5}_{I}^{\text{p}}{\psi}_{I}^{\text{p}}.$ | (13.45) |

The effective potentials ${v}_{\text{eff}}$ and ${v}_{\text{eff}}$ are obtained by taking the derivative of the total energy expression in Eq. (13.44) with respect to electron density and proton density, respectively. Analogous to NEO-HF, these electronic and protonic Kohn-Sham orbitals are expanded as linear combinations of electronic or protonic Gaussian basis functions (${\varphi}_{\mu}^{\text{e}}({\mathbf{r}}_{\text{e}})$ and ${\varphi}_{{\mu}^{\prime}}^{\text{p}}({\mathbf{r}}_{\text{p}})$). The extension to open-shell electron systems is analogous to the NEO-UHF method.

The practical implementation of the NEO-DFT method requires an electron-electron exchange-correlation functional, a proton-proton exchange-correlation functional, and an electron-proton correlation functional. Any conventional electron-electron exchange-correlation functional can be used within the NEO-DFT framework.^{Brorsen:2018} Because the proton-proton exchange and correlation are negligible in molecular systems, only the exchange at the NEO-Hartree-Fock level is included to eliminate self-interaction error in the NEO-DFT method. A suitable electron-proton correlation functional is essential for obtaining an accurate proton densities and energies, and the epc17-2^{Yang:2017, Brorsen:2017} and epc19^{Tao:2019} functionals are designed to achieve this goal. These two functionals are based on the multicomponent extension of the Colle-Salvetti formalism. The epc17-2 functional is of the local density approximation (LDA) type with the functional form:

${E}_{\text{epc}}[{\rho}^{\text{e}},{\rho}^{\text{p}}]=-{\displaystyle \int \mathit{d}\mathbf{r}\frac{{\rho}^{\text{e}}(\mathbf{r}){\rho}^{\text{p}}(\mathbf{r})}{a-b{[{\rho}^{\text{e}}(\mathbf{r}){\rho}^{\text{p}}(\mathbf{r})]}^{1/2}+c{\rho}^{\text{e}}(\mathbf{r}){\rho}^{\text{p}}(\mathbf{r})}}.$ | (13.46) |

The epc19 functional is its multicomponent generalized gradient approximation (GGA) extension that depends on the electron and proton density gradients and is of the form:

${E}_{\text{epc}}[{\rho}^{\text{e}},{\rho}^{\text{p}},\nabla {\rho}^{\text{e}},\nabla {\rho}^{\text{p}}]=-{\displaystyle \int}d\mathbf{r}{\displaystyle \frac{{\rho}^{\text{e}}(\mathbf{r}){\rho}^{\text{p}}(\mathbf{r})}{a-b{[{\rho}^{\text{e}}(\mathbf{r}){\rho}^{\text{p}}(\mathbf{r})]}^{1/2}+c{\rho}^{\text{e}}(\mathbf{r}){\rho}^{\text{p}}(\mathbf{r})}}\times $ | (13.47) | |||

$\left\{1-d\left({\displaystyle \frac{{[{\rho}^{\text{e}}(\mathbf{r}){\rho}^{\text{p}}(\mathbf{r})]}^{-1/3}}{{(1+{m}_{\text{p}})}^{2}}}\left[{m}_{\text{p}}^{2}{\displaystyle \frac{{\nabla}^{2}{\rho}^{\text{e}}(\mathbf{r})}{{\rho}^{\text{e}}(\mathbf{r})}}-2{m}_{\text{p}}{\displaystyle \frac{\nabla {\rho}^{\text{e}}(\mathbf{r})\cdot \nabla {\rho}^{\text{p}}(\mathbf{r})}{{\rho}^{\text{e}}(\mathbf{r}){\rho}^{\text{p}}(\mathbf{r})}}+{\displaystyle \frac{{\nabla}^{2}{\rho}^{\text{p}}(\mathbf{r})}{{\rho}^{\text{p}}(\mathbf{r})}}\right]\right)\text{exp}\left[{\displaystyle \frac{-k}{{[{\rho}^{\text{e}}(\mathbf{r}){\rho}^{\text{p}}(\mathbf{r})]}^{1/6}}}\right]\right\}$ |

In addition to the parameters $a$, $b$, and $c$ in the epc17-2 functional,^{Brorsen:2017} the epc19 functional^{Tao:2019} has the additional $d$ and $k$ parameters and also depends on the proton mass ${m}_{\text{p}}$.
Analogous to the NEO-HF analytic energy gradients, the NEO-DFT analytic gradients are also available for these two functionals.