Kohn-Sham equations from variational principle I'm trying to understand how the Kohn-Sham equations arise from the variational principle, failing.
I think my problem is the inability to apply the variational principle. Or, I lack some crucial understanding about many body physics. (Background: I've never studied many body physics, but basic quantum mechanics and classical electrodynamics, yes. I run into DFT while searching for ways to model biochemical reactions.)
So the Kohn Sham equations are the time-independent Schrödinger equations (one for each electron)
\begin{equation}
\big[ -\frac{\nabla^2}{2} + v_{\mathrm{eff}}(r) \big] \phi_i = \epsilon_i \phi_i.
\end{equation}
In density functional theory the energy functional can be written (for non-interacting electrons that have the same electron density as the exact system)
\begin{align}
E_v[n(\vec{r})] =& T_s[n_s(\vec{r})] + \int \mathrm{d}^3r \, v_{\mathrm{eff}}(\vec{r}) n_s(\vec{r}),
\end{align}
where $n(\vec{r})$ is the (trial) electron density, $T_s$ is the kinetic energy and the integral is the potential energy of the electrons moving in a potential $v(r)$. 
With a fixed number of $N$ electrons, the method of Lagrange multipliers ought to have the form (as far as I've understood)
\begin{equation}
\frac{ \partial E[n] }{ \partial n(\vec{r}) } = \lambda \frac{\partial}{\partial n(\vec{r})} \int \mathrm{d}^3r \, n(\vec{r})
\end{equation}
I would start to expand the right hand side as 
\begin{align}
\frac{\partial}{\partial n(\vec{r})} \int \mathrm{d}^3r \, n(\vec{r}) =& \int \mathrm{d}^3r \frac{ \partial n(\vec{r}) }{ \partial n(\vec{r}) } \\
=& \int \mathrm{d}^3 r 
\end{align}
And the left hand side as
\begin{align}
& \int \mathrm{d}^3r \, \bigg\{ \big(\frac{\partial}{\partial n}v(\vec{r})\big) n(\vec{r})  + v(\vec{r}) \frac{\partial}{\partial n}n(\vec{r}) \bigg\} + \frac{\partial}{\partial n}T_s[n] \\
=& \int \mathrm{d}^3 r \, v(\vec{r}) + \frac{\partial}{\partial n}T_s[n] 
\end{align}
What am I doing wrong? What should I be doing?
 A: After studying variational calculus for a while, here is what I gathered. 
The variation of a functional $J[L(x,f,f')] = \int \, L(x,f,f')\mathrm{d}x$ can be written as
\begin{equation}
\delta J[L] = \int \frac{ \delta J[f] }{ \delta f(x) } \delta f(x) \mathrm{d}x
\end{equation}
with
\begin{equation}
\frac{\delta J}{\delta f(x)} = \frac{\partial L}{\partial f} - \frac{\mathrm{d}}{\mathrm{d}x} \frac{\partial L}{\partial f'}
\end{equation}
Here, $L(x,f,f') = E_v[n(r)] - \lambda\big( \int n(r) \mathrm{d}^3 r - N \big)$, a function of r and n(r) (not n'). 
In my question I tried to apply the required but not sufficient condition: $\nabla_{r,n(r)} E_v[n(\vec{r})] = \lambda \nabla_{r,n(r)} \big(\int n(\vec{r}) \mathrm{d}^3r - N \big)$.
Instead, if I take the variation of the Lagrangian $E_v[n] - \lambda \big(\int n(\vec{r}) \mathrm{d}^3r - N \big)$ of the fictious system:
\begin{align}
&\delta \bigg(E_v[n(r)] - \lambda \big(\int n(r) \mathrm{d}^3r - N \big) \bigg) \\
=& \int \frac{\delta E_v[n]}{\delta n(r)} \delta n(r) \mathrm{d}^3r - \lambda \int \bigg( \frac{\delta (\int n(r) \mathrm{d}^3 r)}{\delta n(r)} \bigg) \delta n(r) \mathrm{d}^3 r \\
=& \int \frac{\delta E_v[n]}{\delta n(r)} \delta n(r) \mathrm{d}^3r - \lambda \int \bigg( \frac{\partial n(r)}{\partial n(r)} \bigg) \delta n(r) \mathrm{d}^3 r \\
\end{align}
Furthermore, after expanding $E_v[n(r)] = \int v(r)n(r) \mathrm{d}^3 r + T_s[n(\vec{r})]$ I get
\begin{align}
& \int \frac{\partial \big(v(r) n(r)\big)}{\partial n(r)} \delta n(r)\mathrm{d}^3 r + \int \frac{\delta T_s[n(r)]}{\delta n(r)} \delta n(r) \mathrm{d}^3 r - \lambda \bigg( \int \delta n(r) \mathrm{d}^3r \bigg) \\
=& \int \big[n(r)\frac{\partial v(r)}{\partial n(r)} + v(r)\frac{\partial n(r)}{\partial n(r)} \big]\delta n(r)\mathrm{d}^3 r + \int \frac{\delta T_s[n(r)]}{\delta n(r)} \delta n(r) \mathrm{d}^3 r - \lambda \bigg( \int \delta n(r) \mathrm{d}^3r \bigg) \\
=& \int \delta n(r) \bigg\{ v(r) + \frac{\delta T_s[n(r)]}{\delta n(r)} - \lambda \bigg\} \mathrm{d}^3 r = 0
\end{align}
Which is what I wanted to get. 
The next question is why do we arrive at the 1D time-independent Schrödinger equation from this?
First, the integral is zero only when the kernel is zero (except for a set of points with zero measure), since the energy functional (which is the expectation value of the system's energy) is positive (by convention) and $\lambda$ is a constant, so that the sign of the integrand is always the same. This means that the integrand must be zero.
After that, we can multiply the part in $\{  \}$-brackets with some function $\phi$, and notice that this is the Schrödinger equation we were looking for. After this, we can define the electron density to be $n(r) = \sum_i |\phi_i|^2$. After iterating enough times, you should get the same $n(r)$ which was your initail guess.
Although I must admit that I'm not a 100% satisfied with that last part, yet. I'm not sure if that's the correct mindset, or if it was known from elsewhere that $n(r) = \sum_i \phi_i$.
