diff --git a/Manuscript/eDFT.tex b/Manuscript/eDFT.tex
index 6981e5c..ae2e3f7 100644
--- a/Manuscript/eDFT.tex
+++ b/Manuscript/eDFT.tex
@@ -50,16 +50,18 @@
 \newcommand{\EPT}{E_\text{PT2}}
 \newcommand{\EFCI}{E_\text{FCI}}
 
-% matrices
+% matrices/operator
 \newcommand{\br}{\bm{r}}
 \newcommand{\bw}{\bm{w}}
 \newcommand{\bG}{\bm{G}}
 \newcommand{\bS}{\bm{S}}
 \newcommand{\bGamma}[1]{\bm{\Gamma}^{#1}}
+\newcommand{\opGamma}[1]{\hat{\Gamma}^{#1}}
 \newcommand{\bHc}{\bm{h}}
 \newcommand{\bF}[1]{\bm{F}^{#1}}
 \newcommand{\Ex}[1]{\Omega^{#1}}
 
+
 % elements
 \newcommand{\ew}[1]{w_{#1}}
 \newcommand{\eG}[1]{G_{#1}}
@@ -172,7 +174,154 @@ Atomic units are used throughout.
 
 \subsection{GOK-DFT}
 
+The GOK ensemble energy is defined as follows:
+\beq
+E^{\bw}=\sum_{K\geq0}w_KE^{(K)},
+\eeq
+where the $K$th energy level $E^{(K)}$ [$K=0$ refers to the ground state] is the eigenvalue of the electronic Hamiltonian $\hat{H}=\hat{h}
+%\sum^N_{i=1}\hat{h}(i) 
++\hat{W}_{\rm ee}$, where $\hat{h}\equiv \sum^N_{i=1}\left(-\frac{1}{2}\nabla_{\br_i}^2+v_{\rm ne}(\br_i)\right)$ is the one-electron operator describing kinetic and nuclear attraction energies, and $\hat{W}_{\rm ee}$ is the electron repulsion operator. 
+The (positive) ensemble weights $w_K$ decrease with increasing index $K$. They are normalized, i.e.
+\beq
+w_0=1-\sum_{K>0}w_K,
+\eeq
+so that only the weights $\bw\equiv\left(w_1,w_2,\ldots w_K,\ldots\right)$ assigned to the excited states can vary independently.  
+For simplicity we will assume in the following that the energies are not degenerate. Note that the theory can be extended to multiplets simply by assigning the same ensemble weight to all degenerate states~\cite{}. In GOK-DFT, the ensemble energy is determined variationally as follows:
+\beq\label{eq:var_ener_gokdft}
+E^{{\bw}}=
+\underset{\opGamma{{\bw}}}{\rm min}\left\{
+{\rm
+Tr}\left[\opGamma{{\bw}}\hat{h}\right]
++
+{E}^{{\bw}}_{\rm
+Hx}\left[n_{\opGamma{\bw}}\right]
++
+{E}^{{\bw}}_{\rm
+c}\left[n_{\opGamma{\bw}}\right]
+\right\},
+\eeq
+where ${\rm
+Tr}$ denotes the trace and the trial ensemble density matrix operator reads
+\beq
+\opGamma{{\bw}}=\sum_{K\geq 0}w_K\ket{\Phi^{(K)}}\bra{\Phi^{(K)}}.
+\eeq
+The determinants (or configuration state functions) $\Phi^{(K)}$ are all constructed from the same set of (ensemble Kohn--Sham) orbitals that is optimized variationally and the trial ensemble density is simply the weighted sum of the individual densities: 
+\beq
+n_{\opGamma{\bw}}(\br)=\sum_{K\geq 0}w_Kn_{\Phi^{(K)}}(\br).
+\eeq
+As readily seen from Eq.~(\ref{eq:var_ener_gokdft}), both Hartree-exchange and 
+correlation energies are described with density functionals that are {\it weight-dependent}. We focus here on the (exact) Hx part which is defined as follows:
+\beq
+{E}^{{\bw}}_{\rm
+Hx}[n]=\sum_{K\geq 0}w_K\bra{\Phi^{(K)}[n]}\hat{W}_{\rm
+ee}\ket{\Phi^{(K)}[n]}
+\eeq
+where the KS wavefunctions fulfill the ensemble density constraint
+\beq
+\sum_{K\geq 0}w_Kn_{\Phi^{(K)}[n]}(\br)=n(\br).
+\eeq 
+The (approximate) description of the correlation part is discussed in Sec.~\ref{sec:eDFA}.\\
+
+In practice, one is not much interested in ensemble energies but rather in excitation energies or individual energy levels (for geometry optimizations, for example). The latter can be extracted exactly as follows~\cite{}: 
+\beq\label{eq:indiv_ener_from_ens}
+E^{(I)}&=&E^{{\bw}}+\sum_{K>0}\left(\delta_{IK}-w_K\right)\dfrac{\partial
+E^{{\bw}}}{\partial w_K},
+\eeq
+where, according to the variational ensemble energy expression of Eq.~(\ref{eq:var_ener_gokdft}),  
+%\subsection{Hybrid GOK-DFT}
+
+
+\subsection{Approximations}
+
+In order to compute (approximate) energy levels within generalized
+GOK-DFT we use a two-step procedure. The first step consists in
+optimizing variationally the ensemble density matrix according to
+Eq.~(\ref{eq:var_princ_Gamma_ens}) with an approximate Hxc ensemble
+functional where (i) the ghost-interaction correction functional $\overline{E}^{{\bw}}_{\rm
+Hx}[n]$ in
+Eq.~(\ref{eq:exact_GIC}) is
+neglected, for simplicity, and (ii) the weight-dependent correlation
+energy is described at the local density level of approximation. 
+At this
+level of approximation, the two correlation functionals $\overline{E}^{{\bw}}_{\rm
+c}[n]$ and ${E}^{{\bw}}_{\rm
+c}[n]$ are actually identical and can be expressed as
+\beq\label{eq:eLDA_corr_fun}
+{E}^{{\bw}}_{\rm
+c}[n]=\int d\br\;n(\br)\epsilon_{c}^{\bw}(n(\br)). 
+\eeq
+More
+details about the construction of such a functional will be given in the
+following. In order to remove ghost interactions from the variational energy
+expression used in the first step, we then employ the (in-principle-exact)
+expression in Eq.~(\ref{eq:exact_ind_ener_OEP-like}). In this second
+step, the response of the individual density matrices to weight
+variations (last term on the right-hand side of
+Eq.~(\ref{eq:exact_ind_ener_OEP-like})) is neglected. The complete GIC
+procedure can be summarized as follows,    
+\beq
+{\bmg}^{\bw}\approx\argmin_{{\bm\gamma}^{\bw}}
+\Big\{
+{\rm
+Tr}\left[{\bm \gamma}^{{\bw}}{\bm h}\right]+W_{\rm
+HF}\left[{\bm\gamma}^{\bw}\right]
++
+{E}^{{\bw}}_{\rm
+c}\left[n_{\bm\gamma^{\bw}}\right]
+%+E^{{\bw}}_{\rm c}\left[n_{\hat{\Gamma}^{{\bw}}}\right]
+\Big\},
+\nonumber\\
+\eeq
+and
+\beq
+E^{(I)}&&\approx{\rm
+Tr}\left[{\bmg}^{(I)}{\bm h}\right]
++\frac{1}{2} \Tr(\bmg^{(I)} \, \bG \,
+\bmg^{(I)})
+\nonumber\\
+&&+{E}^{{\bw}}_{\rm
+c}\left[n_{\bmg^{\bw}}\right]
++\int d\br\,\dfrac{\delta {E}^{{\bw}}_{\rm
+c}\left[n_{\bmg^{\bw}}\right]}{\delta
+n({\br})}\left(n_{\bmg^{(I)}}(\br)-n_{\bmg^{\bw}}(\br)\right)
+\nonumber\\
+&&+\sum_{K>0}\left(\delta_{IK}-w_K\right)\left.       \dfrac{\partial {E}^{{\bw}}_{\rm
+c}\left[n\right]}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}}
+,
+\eeq
+thus leading to the final implementable expression [see Eq.~(\ref{eq:eLDA_corr_fun})]
+\beq
+E^{(I)}&&\approx{\rm
+Tr}\left[{\bmg}^{(I)}{\bm h}\right]
++\frac{1}{2} \Tr(\bmg^{(I)} \, \bG \,
+\bmg^{(I)})
+\nonumber\\
+&&+\int d\br\,
+{\epsilon}^{{\bw}}_{\rm
+c}(n_{\bmg^{\bw}}(\br))\,n_{\bmg^{(I)}}(\br)
+\nonumber\\
+&&
++\int d\br\,\left.\dfrac{\partial {\epsilon}^{{\bw}}_{\rm
+c}(n)}{\partial n}\right|_{n=n_{\bmg^{\bw}}(\br)}n_{\bmg^{\bw}}(\br)\left(n_{\bmg^{(I)}}(\br)-n_{\bmg^{\bw}}(\br)\right)
+\nonumber\\
+&&
++\int d\br\,\sum_{K>0}\left(\delta_{IK}-w_K\right)n_{\bmg^{\bw}}(\br)\left.
+\dfrac{\partial {\epsilon}^{{\bw}}_{\rm
+c}(n)}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}(\br)}.
+\eeq
+
+\blue{$================================$}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Theory (old)}
+\label{sec:eDFT_old}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Kohn--Sham formulation of GOK-DFT}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \subsection{Hybrid GOK-DFT}
+\label{sec:geKS}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 Since Hartree and exchange energy contributions cannot be separated in
 the one-dimensional case, we introduce in the following an alternative
@@ -348,7 +497,7 @@ ${\bmg}^{\bw}$ and ${\bm h}={\bm t}+{\bm v}_{\rm ext}$ is the total one-electron
 Hamiltonian matrix representation. When the minimum is reached, the
 ensemble energy and its derivatives can be used to extract individual
 ground- and excited-state energies as follows:\cite{Deur_2018b} 
-\beq\label{eq:indiv_ener_from_ens}
+\beq
 E^{(I)}&=&E^{{\bw}}+\sum_{K>0}\left(\delta_{IK}-w_K\right)\dfrac{\partial
 E^{{\bw}}}{\partial w_K}.
 \eeq
@@ -489,120 +638,6 @@ Hxc}(n)}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}(\br)}.
 \eeq
 }
 
-\subsection{Approximations}
-
-In order to compute (approximate) energy levels within generalized
-GOK-DFT we use a two-step procedure. The first step consists in
-optimizing variationally the ensemble density matrix according to
-Eq.~(\ref{eq:var_princ_Gamma_ens}) with an approximate Hxc ensemble
-functional where (i) the ghost-interaction correction functional $\overline{E}^{{\bw}}_{\rm
-Hx}[n]$ in
-Eq.~(\ref{eq:exact_GIC}) is
-neglected, for simplicity, and (ii) the weight-dependent correlation
-energy is described at the local density level of approximation. 
-At this
-level of approximation, the two correlation functionals $\overline{E}^{{\bw}}_{\rm
-c}[n]$ and ${E}^{{\bw}}_{\rm
-c}[n]$ are actually identical and can be expressed as
-\beq\label{eq:eLDA_corr_fun}
-{E}^{{\bw}}_{\rm
-c}[n]=\int d\br\;n(\br)\epsilon_{c}^{\bw}(n(\br)). 
-\eeq
-More
-details about the construction of such a functional will be given in the
-following. In order to remove ghost interactions from the variational energy
-expression used in the first step, we then employ the (in-principle-exact)
-expression in Eq.~(\ref{eq:exact_ind_ener_OEP-like}). In this second
-step, the response of the individual density matrices to weight
-variations (last term on the right-hand side of
-Eq.~(\ref{eq:exact_ind_ener_OEP-like})) is neglected. The complete GIC
-procedure can be summarized as follows,    
-\beq
-{\bmg}^{\bw}\approx\argmin_{{\bm\gamma}^{\bw}}
-\Big\{
-{\rm
-Tr}\left[{\bm \gamma}^{{\bw}}{\bm h}\right]+W_{\rm
-HF}\left[{\bm\gamma}^{\bw}\right]
-+
-{E}^{{\bw}}_{\rm
-c}\left[n_{\bm\gamma^{\bw}}\right]
-%+E^{{\bw}}_{\rm c}\left[n_{\hat{\Gamma}^{{\bw}}}\right]
-\Big\},
-\nonumber\\
-\eeq
-and
-\beq
-E^{(I)}&&\approx{\rm
-Tr}\left[{\bmg}^{(I)}{\bm h}\right]
-+\frac{1}{2} \Tr(\bmg^{(I)} \, \bG \,
-\bmg^{(I)})
-\nonumber\\
-&&+{E}^{{\bw}}_{\rm
-c}\left[n_{\bmg^{\bw}}\right]
-+\int d\br\,\dfrac{\delta {E}^{{\bw}}_{\rm
-c}\left[n_{\bmg^{\bw}}\right]}{\delta
-n({\br})}\left(n_{\bmg^{(I)}}(\br)-n_{\bmg^{\bw}}(\br)\right)
-\nonumber\\
-&&+\sum_{K>0}\left(\delta_{IK}-w_K\right)\left.       \dfrac{\partial {E}^{{\bw}}_{\rm
-c}\left[n\right]}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}}
-,
-\eeq
-thus leading to the final implementable expression [see Eq.~(\ref{eq:eLDA_corr_fun})]
-\beq
-E^{(I)}&&\approx{\rm
-Tr}\left[{\bmg}^{(I)}{\bm h}\right]
-+\frac{1}{2} \Tr(\bmg^{(I)} \, \bG \,
-\bmg^{(I)})
-\nonumber\\
-&&+\int d\br\,
-{\epsilon}^{{\bw}}_{\rm
-c}(n_{\bmg^{\bw}}(\br))\,n_{\bmg^{(I)}}(\br)
-\nonumber\\
-&&
-+\int d\br\,\left.\dfrac{\partial {\epsilon}^{{\bw}}_{\rm
-c}(n)}{\partial n}\right|_{n=n_{\bmg^{\bw}}(\br)}n_{\bmg^{\bw}}(\br)\left(n_{\bmg^{(I)}}(\br)-n_{\bmg^{\bw}}(\br)\right)
-\nonumber\\
-&&
-+\int d\br\,\sum_{K>0}\left(\delta_{IK}-w_K\right)n_{\bmg^{\bw}}(\br)\left.
-\dfrac{\partial {\epsilon}^{{\bw}}_{\rm
-c}(n)}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}(\br)}.
-\eeq
-
-\blue{$================================$}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{Theory (old)}
-\label{sec:eDFT_old}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Kohn--Sham formulation of GOK-DFT}
-
-Let us start from the analog for ensembles of Levy's universal
-functional,   
-\beq\label{eq:ens_LL_func}
-F^{\bw}[n]&=&
-\underset{\hat{\gamma}^{{\bw}}\rightarrow n}{\rm min}\left\{{\rm
-Tr}\left[\hat{\gamma}^{{\bw}}\left(\hat{T}+\hat{W}_{\rm
-ee}\right)\right]\right\}
-\eeq
-where ${\rm
-Tr}$ denotes the trace. The minimization over trial ensemble density matrix operators
-$\hat{\gamma}^{{\bw}}=\sum_{K\geq0}w_K\vert\Psi^{(K)}\rangle\langle\Psi^{(K)}\vert$
-is performed under the following density constraint:
-\beq
-{\rm
-Tr}\left[\hat{\gamma}^{{\bw}}\hat{n}(\br)\right]=\sum_{K\geq0}w_Kn_{\Psi^{(K)}}(\br)=n(\br),
-\eeq
-where $\hat{n}(\br)$ is the density operator, $n_{\Psi}$ denotes the
-density of wavefunction $\Psi$, and
-$\bw\equiv\left(w^{(1)},w^{(2)},\ldots\right)$ is the collection of
-(decreasing) ensemble weights assigned to the excited states. Note that
-$w_0=1-\sum_{K>0}w_K\geq 0$.
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Hybrid GOK-DFT}
-\label{sec:geKS}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-
 \subsection{Exact ensemble exchange in hybrid GOK-DFT}