deletions | additions       diff --git a/RDMFT.tex b/RDMFT.tex  index 9952b47..afa25bb 100644  --- a/RDMFT.tex  +++ b/RDMFT.tex 
...
\end{eqnarray}  where $\theta$ is the unit-step Heavside function. We initialize the whole matrix as $F_{ji}=(\lambda_{ji}+\lambda_{ij}^{*})/2$. In every iteration we diagonalize $\mathbf{F}$ keeping the diagonal elements for the next iteration while changing the off-diagonal ones to (\ref{eqF}). At the solution all off-diagonal elements of this matrix vanish, hence, the matrices $\mathbf{F}$ and $\gamma$ can be brought simultaneously to a diagonal form. Thus, the $\{\phi_i\}$ which are the solutions of Eq. (\ref{eqlambda}) can be found by diagonalization of $\mathbf{F}$ in an iterative manner (see \cite{Piris} for details). The criterium to exit the natural orbital optimization is that the difference in the total energies calculated in two successive $\mathbf{F}$ diagonalizations is smaller than a threshhold. Overall convergence is achieved when the difference in the total energies in two successive occupation-number optimizations and the non-diagonal matrix elements of $\mathbf{F}$, see Eq.\ (\ref{eqF}), are close to zero.  As already mentioned above, one needs an initial guess for the natural orbitals both for the first step of occupation-number optimization but also for the optimization of the natural orbitals. A rather obvious choice would be the occupied and a few unoccupied orbitals resulting from a DFT or HF calculation. Unfortunatelly, as we found out, there are unbound states among the HF (or DFT) unoccupied states which provide a very bad starting point for the weakly occupied natural orbitals. This is not surprising as especially the higher unoccupied HF or DFT orbitals see in practice no potential. In an implementation on a finite grid they act like particles in a box. OEP gives more bound states than HF or DFT. This is because all unoccupied HF states suffer from self interaction which makes them quiet unlocalized and DFT effective potential decays exponentially giving more unbound states than OEP. As a result, using HF or DFT orbitals to start a RDMFT calculation, the natural orbitals did not converge to any reasonable shape but even starting with OEP one needs to further localize the unoccupied states. Thus what we use at the moment as initial guess for the RDMFT natural orbitals is OEP orbitals which are further localized with a gaussian exponential function. The unbound states need to be more localized thus we use a bigger exponent than for the bound ones. One has to be carefull not to break degeneracies, i.e. when two orbitals are nearly degenerate they should be localized using the same exponent.In Fig. \ref{fig:rdmft_h2} we show the dissociation curve of H$_{2}$ obtained with Octopus and we compare it with RDMFT results using a gaussian basis set code named HIPPO created by N.N. Lathiotakis \cite{ML2008}. \cite{ML_Benchmark}.  It shows the energy curve for the hydrogen molecule as a function of the distance $R$ between the two hydrogen atoms. For this calculation we kept thirtheen natural orbitals with the smalles occupation number being of the order of $10^{-6}$ after the RDMFT calculation had converged. The curve looks the same as in other implementations of RDMFT \cite{GPB2005} keeping the nice feature of not diverging in the dissociation limit. Therefore, apart from the problem of the initial guess for the weakly occupied orbitals the RDMFT implementation is working. Hopefully, this last remaining problem will be solved in the very near future. In addition, we plan to include an open-shell version of RDMFT and additional approximate functionals. Since many available functionals are of very similar form as the M\"uller functional, only differing in their dependence on the occupation numbers, their implementation will be rather straightforward.