Density Functional Theory with Generalized Gradient Approxmations (GGA) - What Happens if the Electron Density Shows a "Kink" at the Ion Position?

My question may be stupid, so please correct me if you find anything which is obviously erroneous. 


In the following I will place a question mark (?) besides points/steps I consider doubtful.


My observations:


1. In density functional theory (DFT) (at zero temperature), the energy of the system $E$ is broken up as $E = \int\_{V} d v \left[ e\_{K}(n(\vec{r})) + e\_{H}(n(\vec{r})) + e\_{XC}(n(\vec{r})) \right]$, where $V$ is the system size, $e$ the energy density per unit volume, $n(\vec{r})$ the electronic density, subscript "$K$" for "kinetic", "$H$" for "Hartree" and "$XC$" for "exchange-correlation".
2. In the generalized gradient approximation (gga), $e$ also depends on $\nabla n(\vec{r})$. Thus write $e(n(\vec{r}), \nabla n(\vec{r}))$.


Assumptions:
1. spherical symmetry (closed-shell)


By spherical symmetry $e\_{XC}(n(\vec{r}), \nabla n(\vec{r})) = e(n(r), \nabla n(r))$. Of course in writing this we are restricting the variation in density $\delta n$ in the **radial** direction only. I'm doubtful about the validity of this step. 


Anyway, proceed as usual. The exchange-correlation potential is defined as the functional derivative of $e\_{XC}$ w.r.t. $\delta n$. Now please point out if there are any mistakes in the following steps:


$$\begin{eqnarray} \delta \int dv e\_{XC}(n(r), n^{\prime}(r)) &=& \int dv \left[ \frac{\partial e\_{XC}}{\partial n(r)} \delta n(r) + \frac{\partial e\_{XC}}{\partial n^{\prime}(r)} \delta n^{\prime}(r) \right] \end{eqnarray}$$,


where $n^{\prime}(r) = \partial\_r n(r)$. The first term is done; consider the second term. $\delta n^{\prime} (r) = \partial\_r \delta n(r)$. Thus we integrate the second term on the right hand side by parts, and obtain


$4 \pi \int\_0^{r\_0} dr \partial\_r \left[ r^2 \frac{\partial e\_{XC}}{\partial n^{\prime}} \delta n(r) \right] - 4 \pi \int\_0^{r\_0} dr \left[ \partial\_r \left( r^2 \frac{e\_{XC}}{n^{\prime}}\right) \delta n(r) \right]$. The first term is the boundary term; ignore it for the time being. The second term seems to give the correction to the exchange correlation potential due to the introduction of gradient-dependence. Note that the second term is equivalent to
$\int\_0^{r\_0} dr 4 \pi r^2 \left[ - \frac{1}{r^2} \partial\_r \left( r^2 \frac{\partial e\_{XC}}{\partial n^{\prime}} \right)\right]$, the square-bracketed term seems to give the exchange-correlation potential (?).


Now $- \frac{1}{r^2} \partial\_r \left( r^2 \frac{\partial e\_{XC}}{\partial n^{\prime}} \right) = - \frac{2}{r} \frac{\partial e\_{XC}}{\partial n^{\prime}} - \partial\_r \frac{\partial e\_{XC}}{\partial n^{\prime}}$. What I worry about is the first term, which seems to diverge at $r=0$ (but it shouldn't; or is it cancelled by a piece in the second term?)


A expression for $e\_{XC}$ may be the becke88 functional (exchange only) (PRA, 38, 3098). $$ \begin{eqnarray}e^{Becke}\_{X} &=& e^{LDA}\_{X}(n(r)) - \beta \sum\_{\sigma} \int n\_{\sigma}^{4/3} \frac{x\_{\sigma}^2}{1+6 \beta x\_{\sigma} \sinh^{-1}x\_{\sigma}} d\vec{r} \end{eqnarray}$$, where $\sigma$ is the spin-index (For closed-shell system $n\_{+} = n\_{-} = n/2$ (?)), $\beta$ is a parameter (=0.0042 a.u.), and $x\_{\sigma} = | \nabla n\_{\sigma}|/n\_{\sigma}^{4/3}$. In spherically symmetric system, can we change $| \nabla n\_{\sigma}|$ to $|n^{\prime}(r)|$?. Now I am not sure about the following claims


1. $\frac{\partial e\_{X}^{Becke}}{\partial n^{\prime}(r)}$ is non-zero at the origin. (?) The reason is that, while it seems that $n^{\prime}(0)$ should be zero, it may not be. Consider helium. When solving the schroedinger eq. for the 1s state under arbitrary potential which tends to $2/r$ at $r=0$, the slope of the radial wave function (excluding the $r$-factor from the volume element) is non-zero (Consider for example the radial wave function of the hydrogen ground state, which behaves as $e^{-r}$). This means $n^{\prime}(0) \neq 0$. But then $\frac{\partial e\_{X}^{Becke}}{\partial n^{\prime}(r)}$ will be non-zero at the origin and thus $- \frac{2}{r} \frac{\partial e\_{XC}}{\partial n^{\prime}} $ will tend to infinity at $r=0$.


This should not be the case because the hartree (i.e. "classical" contribution by all electrons) contribution to the potential is small as $r \rightarrow 0$, and the exchange contribution should cancel the "self-action" in the hartree part.


Now, sorry for the long post. Note that I might have done something very wrong, such that the whole argument does not stand. Thanks for your time.


