deletions | additions       diff --git a/RDMFT1.tex b/RDMFT1.tex  index 0250b01..a7075c8 100644  --- a/RDMFT1.tex  +++ b/RDMFT1.tex 
...
% \sum_i \int d\mathbf{r} \delta \phi_i(\mathbf{r})\Big[\frac{\delta E}{\delta \phi_{i}(\mathbf{r})}-\sum_{k}\lambda_{ki}\phi_{k}^{*}(\mathbf{r})\Big] = 0  %\end{eqnarray}  The variation of $\Omega$ is done in two steps: for a fixed set of orbitals the energy functional is minimized with respect to occupation numbers and accordingly for a fixed set of occupations the energy functional is minimized with respect to variations of the orbitals until overall convergenceof the absolute value of the energy  is achieved. As a starting point we use results from a Hartree-Fock calculation and first optimize the occupation numbers. Since the correct $\mu$ is not known it is determined via bisection: for every $\mu$ the objective functional is minimized with respect to $\vartheta_i$ until the condition \ref{eqsumocc} is satisfied.\par Due to the dependence on the occupation numbers, the natural orbital minimization does not lead to an eigenvalue equation like in Hartree-Fock. The implementation of the natural orbital minimization follows the method by Piris and Ugalde (\cite{Piris}). As one can show for fixed occupation numbers one obtains  \begin{eqnarray}  \lambda_{ji} = n_i\langle\phi_j|-\frac{\nabla^2}{2}+v_{ext}|\phi_i\rangle 
...
\begin{eqnarray}  F_{ji}=\theta(i-j)(\lambda_{ji}-\lambda^{*}_{ij})+\theta(j-i)(\lambda^{*}_{ij}-\lambda_{ji}),\label{eqF}  \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 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 overall convergence criteria is the convergence of the relative total energy and the matrix elements of $\mathbf{F}$ \label{eqF} to be close enough 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. In contrast to this the weakly occupied natural orbitals are localized around the atoms. As a result, using HF or DFT orbitals to start a RDMFT calculation, the natural orbitals did not converge to any reasonable shape (some of the weakly occupied natural orbitals were still unbound). In order to check that the only problem of the RDMFT implementation in octopus is the appropriate choice of an initial guess, we performed an RDMFT calculation using HF orbitals from a Gaussian basis set code as our starting point. This provides more localized unoccupied orbitals since they have to be linear combination of Gaussians. With this initial guess we found natural orbitals that looked reasonable. The process of finding an appropriate initial guess for the weakly occupied orbitals without having to rely on a different code is still ongoing. However, we are sure this is the only remaining problem in the implementation. In Fig. \ref{fig:rdmft_h2} we show the results obtained with the initial guess from the Gaussian implementation. 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 weakest occupied having occupancy of the order of $10^-6$ after RDMFT had converged. The curve looks the same as in other implementations of RDMFT (cite..) 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.