deletions | additions
diff --git a/RDMFT1.tex b/RDMFT1.tex
index d347487..8a51787 100644
--- a/RDMFT1.tex
+++ b/RDMFT1.tex
...
In RDMFT the total energy is given by
\begin{eqnarray}
E=\sum_{i=1}^\infty\int d\mathbf{r} n_{i}\phi^{*}_{i}(\mathbf{r})\left(-\frac{\nabla^2}{2}\right) \phi_{i}(\mathbf{r})+\sum_{i=1}^\infty \int d\mathbf{r} V_{\mathrm{ext}}(\mathbf{r})n_{i}|\phi_{i}(\mathbf{r})|^{2}\\
+\frac{1}{2}\sum_{i,j=1}^\infty n_{i} n_{j}\int d\mathbf{r} d\mathbf{r'} \frac{|\phi_{i}(\mathbf{r})|^{2} |\phi_{j}(\mathbf{r})|^{2}}{|\mathbf{r}-\mathbf{r'}|} +
E_{xc}\left[\{n_{j}\},\{\phi_{j}\}\right].\label{eqenergy} E_{xc}\left[\{n_{j}\},\{\phi_{j}\}\right].
\label{eqenergy}
\end{eqnarray}
The third term is the Hartree energy, E_H and the forth the exchange-correlation energy $E_{xc}$. As in DFT, the exact functional of RDMFT is unknown. However, the only part that needs to be approximated $E_{xc}[\gamma]$ comes, contrary to DFT, only from the electron-electron interaction, as the interacting kinetic energy can be explicitely expressed in terms of $\gamma$.
Different approximate functionals are employed and minimized with respect to the 1-RDM in order to find the ground state energy (cite BBC, Piris, and Nek here). A common approximation for $E_{xc}$ is the M\"uller functional \cite{Mueller_1984}, which has the following form
...
and is the only $E_{xc}$ implemented in octopus for the moment.
\par
For closed-shell systems the necessary and sufficient conditions for the 1-RDM to be $N$-representable \cite{Coleman_1963}, i.e.\ to correspond to a $N$-electron wavefunction is that $ 0 \leq n_{i} \leq 2$ and
\begin{eqnarray}
\sum_{i=1}^{\infty}n_{i}=N.\label{eqsumocc}
\end{eqnarray} %\begin{eqnarray}
%\sum_{i=1}^{\infty}n_{i}=N.\label{eqsumocc}
%\end{eqnarray}
Note that within the RDMFT implementation in octopus only closed-shell systems are treated at the moment. Minimization of the energy functional of Eq. (\ref{eqenergy}) is performed under the $N$-representability constraints and the orthonormality requierements of the natural orbitals,
\begin{eqnarray}
\langle %\begin{eqnarray}
%\langle \phi_{i} | \phi_{j}\rangle = \delta_{ij}. \label{eqorth}
\end{eqnarray} %\end{eqnarray}
The bounds on the occupation numbers are automatically satisfied by setting $n_{i}=2\sin^2(2\pi\vartheta_i)$ and varying $\vartheta_{i}$ without constraints. The conditions (\ref{eqsumocc}) and (\ref{eqorth}) are taken into account via Lagrange multipliers $\mu$ and $\lambda_{ij}$, respectively. Then, one can define the following functional
\begin{eqnarray}
\Omega(N,\{\vartheta\} ,\{\phi_i(\mathbf{r})\})= E - \mu \left(\sum_{i=1}^\infty 2sin^2( 2\pi\vartheta_i)-N\right)-\sum_{i,j=1}^\infty \lambda_{ji}(\langle\phi_i|\phi_j\rangle-\delta_{ij})
...
\end{eqnarray}
where $\theta$ is the unit-step Heavside function. We initialize the whole matrix as $F_{ji}=(\lambda_{ji}+\lambda_{ij}^{*})/2$ diagonalize it and keep 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).
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.\ 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. 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.
