deletions | additions       diff --git a/RDMFT.tex b/RDMFT.tex  index e398ad1..c546ec6 100644  --- a/RDMFT.tex  +++ b/RDMFT.tex 
...
\section{Real-space reduced density matrix functional theory}  An alternative approach to DFT, that can model electrons using a single-particle framework is reduced density matrix functional theory (RDMFT) \cite{Gilbert_1975}. Here we present the current results of an ongoing effort to develop a real-space version of RDMFT and to implement it in the Octopus code.  Within RDMFT the total energy of a system is given as a functional of the one body reduced density matrix (1-RDM) %  \begin{equation}  \gamma(\mathbf{r},\mathbf{r'})=N\int\cdots\int d\mathbf{r_2}...d\mathbf{r_N} \Psi^*(\mathbf{r'},\mathbf{r_2}...\mathbf{r_N})\Psi(\mathbf{r},\mathbf{r_2}...\mathbf{r_N}) 
...
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} \begin{equation}  \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. \end{equation}  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} \begin{equation}  \langle \phi_{i} | \phi_{j}\rangle = \delta_{ij}. \label{eqorth}  \end{eqnarray} \end{equation}  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_i\} ,\{\phi_i(\mathbf{r})\})= E - \mu \left(\sum_{i=1}^\infty 2\sin^2( 2\pi\vartheta_i)-N\right)-\sum_{i,j=1}^\infty \lambda_{ji}\left(\langle\phi_i|\phi_j\rangle-\delta_{ij}\right) 
...
% \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 and, accordingly,  for a fixed set of occupations the energy functional is minimized with respect to variations of the orbitals until overall convergence 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} \begin{equation}  \lambda_{ji} = n_i\langle\phi_j|-\frac{\nabla^2}{2}+v_{ext}|\phi_i\rangle  +\int d\mathbf{r} \frac{\delta E_{Hxc}}{\delta \phi_i^{*}(\mathbf{r})}\phi_k^{*}(\mathbf{r}).  \end{eqnarray} \end{equation}  At the extremum, the matrix of the Lagrange multiplyers must be Hermitian, i.e.   \begin{eqnarray}  \lambda_{ji}-\lambda_{ij}^{*}=0\label{eqlambda} 
...
\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, this option is not good mainly because 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 the correspondning orbitals look like those for a particle in a box. Using the exact-exchange approximation (EXX) in an optimized effective potential framework results in a larger number of bound states than HF or the local density approximation (LDA) due to the EXX functional being self-interaction free for both occupied and unoccupied orbitals and the resulting effective potential having the correct asymptotic behavior. Using HF or LDA orbitals to start a RDMFT calculation, the natural orbitals did not converge to any reasonable shape but even when starting from EXX one needs to further localize the unoccupied states. Thus, we multiply each unoccupied orbital with gaussian exponential functions which are centered at the positions of the atoms. The unbound states need to be more localized than the bound ones thus we choose a larger exponent for them. All orbitals are orthogonalized after this multiplication process. The multiplications with gaussian exponentials has to be done such that one does not break energetic degeneracies, i.e.\ when two orbitals are energetically degenerate or 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 RDMFT in Octopus and compare it with results obtained by a gaussian the Gaussian  basis set RDMFT  code named HIPPO created by N.N. Lathiotakis \cite{ML_Benchmark}. HIPPO~\cite{ML_Benchmark}.  For the octopus Octopus  calculation we kept thirtheen 13  natural orbitals with the smallest occupation number being of the order of $10^{-5}$ after the RDMFT calculation had converged. The HIPPO calculation was performed using $30$ natural orbitals. The RDMFT curve obtained with octopus looks similar to the one from HIPPO and other gaussian Gaussian  implementations of RDMFT \cite{GPB2005} keeping the nice feature of not diverging strongly in the dissociation limit. However,especially  for internuclear distances $R$ bigger than $1 a.u$ our curve the real-space energy  lies above the HIPPO one. We believe that the remaining difference can be removed by further improving the initial guess for the orbitals that we use in octopus Octopus,  because a trial calculation using HF orbitals from a Gaussian implementation showed a curve almost identical to the one from the HIPPO code (not shown in the figure). In the future, we plan to include an support for  open-shell version of RDMFT systems  and additional approximate xc  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.