Finding mathematically the ground state density in DFT

Question

To find the ground state density in DFT, you set the following Lagrangian:

$$L = E[\rho(\vec r)] - \mu\left(\int \rho(\vec r) \mathrm{d}\vec r - N\right)$$

While minimising with respect to the density, we get:

$$\frac{\partial L}{\partial \rho} = \frac{\partial E[\rho(\vec r)]}{\partial \rho} - \mu = 0$$

My problem is with the explicit derivation of the energy functional. Let's consider the Thomas-Fermi-Weizsacker kinetic energy functional:

$$T_s(\rho) = C_F \int \rho^{5/3}(\vec r)\mathrm{d}\vec r + C\int \frac{| \nabla \rho (\vec r)|^2}{\rho (\vec r)}\mathrm{d}\vec r$$

For example, how do you derive this with respect to the density?

$$\frac{\partial T_s(\rho)}{\partial \rho} = C_F \frac{\partial}{\partial \rho}\int \rho^{5/3}(\vec r)\mathrm{d}\vec r + C\frac{\partial}{\partial \rho}\int \frac{| \nabla \rho (\vec r)|^2}{\rho (\vec r)}\mathrm{d}\vec r$$

For the first term, you could use Leibniz integral rule, since the density very far away from the nucleus is zero. However, I really have no idea how to tackle the second term. Moreover, what is the algorithm afterwards to find the ground state density? How does one figure out this density from the equation?

Answer 1

What we have to do here is a functional derivative. Lets consider the following method of finding functional derivartives:

$$F[\rho+\delta \rho]-F[\rho] = \int_\Omega \delta\rho \frac{\delta F[\rho]}{\delta\rho}\mathrm{d}\mathbf{r}+\mathscr{O}(\delta\rho^2)\tag{1}\label{FuncDer}$$

Lets now consider the Thomas-Fermi kinetic energy functional using the above equation. The Thomas-Fermi functional is given as:

$$T_\mathrm{TF}[\rho] = C_\mathrm{TF}\int_\Omega\rho^{5/3}\mathrm{d}\mathbf{r}\tag{2}$$

Lets now try to evaluate $T_\mathrm{TF}[\rho+\delta\rho]$:

$$T_\mathrm{TF}[\rho+\delta\rho] = C_\mathrm{TF}\int_\Omega\left(\rho+\delta\rho\right)^{5/3}\mathrm{d}\mathbf{r}\tag{3}$$

Now lets do Taylor expansion of the integrand:

$$T_\mathrm{TF}[\rho+\delta\rho] = C_\mathrm{TF}\int_\Omega\rho^{5/3}+\frac{5}{3}\rho^{2/3}\delta\rho+\mathscr{O}(\delta \rho^2)\mathrm{d}\mathbf{r}\tag{4.1}$$ $$T_\mathrm{TF}[\rho+\delta\rho]= C_\mathrm{TF}\int_\Omega\rho^{5/3}\mathrm{d}\mathbf{r}+\int_\Omega C_\mathrm{TF}\frac{5}{3}\rho^{2/3}\delta\rho\mathrm{d}\mathbf{r}+\mathscr{O}(\delta \rho^2)\tag{4.2}$$

Now subtracting $T_\mathrm{TF}[\rho]$:

$$T_\mathrm{TF}[\rho+\delta\rho] - T_\mathrm{TF}[\rho]=\int_\Omega C_\mathrm{TF}\frac{5}{3}\rho^{2/3}\delta\rho\mathrm{d}\mathbf{r}+\mathscr{O}(\delta \rho^2)\tag{5}$$

Thus by comparing with Eq. (\ref{FuncDer}) it can be seen that:

$$\frac{\delta T_\mathrm{TF}[\rho]}{\delta\rho}=C_\mathrm{TF}\frac{5}{3}\rho^{2/3}\tag{6}$$

Now for the second part of the kinetic energy functional in your question, i.e. the Weizsacker kinetic energy functional:

$$T_\mathrm{W}[\rho]=C\int_\Omega\frac{\nabla\rho\cdot\nabla\rho}{\rho}\mathrm{d}\mathbf{r}\tag{7}$$

Again lets evaluate $T_\mathrm{W}[\rho+\delta\rho]$:

$$T_\mathrm{W}[\rho+\delta\rho]=C\int_\Omega\frac{\nabla(\rho+\delta\rho)\cdot\nabla(\rho+\delta\rho)}{\rho+\delta\rho}\mathrm{d}\mathbf{r}\tag{8}$$

Now doing a Taylor expansion again:

$$T_\mathrm{W}[\rho+\delta\rho] = C\int_\Omega \frac{\nabla\rho\cdot\nabla\rho}{\rho} + \frac{-2\nabla^2\rho+\frac{\nabla\rho\cdot\nabla\rho}{\rho}}{\rho}\delta\rho+\mathscr{O}(\delta\rho^2)\mathrm{d}\mathbf{r}\tag{9.1}$$

$$T_\mathrm{W}[\rho+\delta\rho] = C\int_\Omega \frac{\nabla\rho\cdot\nabla\rho}{\rho}\mathrm{d}\mathbf{r} + \int_\Omega C\left(-\frac{2\nabla^2\rho}{\rho}+\frac{\nabla\rho\cdot\nabla\rho}{\rho^2}\right)\delta\rho\mathrm{d}\mathbf{r}+\mathscr{O}(\delta\rho^2)\tag{9.2}$$

Now subtracting $T_\mathrm{W}[\rho]$:

$$T_\mathrm{W}[\rho+\delta\rho] - T_\mathrm{W}[\rho] = \int_\Omega C\left(-\frac{2\nabla^2\rho}{\rho}+\frac{\nabla\rho\cdot\nabla\rho}{\rho^2}\right)\delta\rho\mathrm{d}\mathbf{r}+\mathscr{O}(\delta\rho^2)\tag{10}$$

Thus by comparing with Eq. (\ref{FuncDer}) it can be seen that:

$$\frac{\delta T_\mathrm{W}[\rho]}{\delta\rho} = C\left(-\frac{2\nabla^2\rho}{\rho}+\frac{\nabla\rho\cdot\nabla\rho}{\rho^2}\right)\tag{11}$$

Now for the second part of your question. For molecular calculations involving kinetic energy functionals you have to look into orbital-free density functional theory (OF-DFT). The kinetic energy functionals are not used in Kohn-Sham density functional theory (KS-DFT). In KS-DFT the kinetic energy are calculated from the Kohn-Sham orbitals instead of the density.

Often when working with OF-DFT the equations can be solved numerically using finte element method (FEM) for a nice introduction to this take a look at Joel Davidsson master thesis [1]. Keep in mind that doing OF-DFT calculations is not a trivial task!

[1] : https://www.diva-portal.org/smash/get/diva2:864857/FULLTEXT01.pdf