deletions | additions       diff --git a/RDMFT.tex b/RDMFT.tex  index 08c3690..8913f2b 100644  --- a/RDMFT.tex  +++ b/RDMFT.tex 
...
\section{Real-space reduced density matrix density-matrix  functional theory} An alternative approach to DFT, 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 RDMFT,  the total energy of a system is given as a functional of the one body one-body  reduced density matrix density-matrix  (1-RDM) %  \begin{equation}  \gamma(\vec{r},\vec{r'})=N\int\cdots\int \mathrm{d}\vec{r_2}\ldots\vec{d}\vec{r}_N\, \Psi^*(\vec{r}',\vec{r}_2...\vec{r}_N)\Psi(\vec{r},\vec{r}_2...\vec{r}_N) 
...
\begin{multline}  \label{eqenergy}  E=-\sum_{i=1}^\infty n_i \int \mathrm{d}\vec{r} \phi^{*}_{i}(\vec{r})\frac{\nabla^2}{2}\phi_{i}(\vec{r})+\sum_{i=1}^\infty n_i \int \mathrm{d}\vec{r}\, v_{\mathrm{ext}}(\vec{r})|\phi_{i}(\vec{r})|^{2}\\  +\frac{1}{2}\sum_{i,j=1}^\infty n_{i} n_{j}\int \mathrm{d}\vec{r} \mathrm{d}\vec{r}' \frac{|\phi_{i}(\vec{r})|^{2} |\phi_{j}(\vec{r})|^{2}}{|\vec{r}-\vec{r}'|} + E_{xc}\left[\{n_{j}\},\{\phi_{j}\}\right]\ E_{\rm xc}\left[\{n_{j}\},\{\phi_{j}\}\right]\  . \end{multline}  The third term is the Hartree energy, $E_H$, $E_{\rm H}$,  and the fourth the xc energy, $E_{xc}$. $E_{\rm xc}$.  As in DFT, the exact functional of RDMFT is unknown. However, the part that needs to be approximated, $E_{xc}[\gamma]$, $E_{\rm xc}[\gamma]$,  comes, contrary to DFT, only from the electron-electron interaction, as the interacting kinetic energy can be explicitly 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{GPB2005,ML2008,P2013}. A common approximation for $E_{xc}$ $E_{\rm xc}$  is the M\"uller functional~\cite{Mueller_1984}, which has the form \begin{multline}  E_{xc}\left[\{n_j\},\{\phi_j\}\right]=\\ E_{\rm xc}\left[\{n_j\},\{\phi_j\}\right]=\\  -\frac{1}{2}\sum_{i,j=1}^\infty \sqrt{n_{i} n_{j}}\iint \mathrm{d}\vec{r} \mathrm{d}\vec{r}' \frac{\phi_{i}^{*}(\vec{r})\phi_i(\vec{r}') \phi_{j}^{*}(\vec{r}')\phi_j(\vec{r})}{|\vec{r}-\vec{r'}|}  \end{multline}  and is the only $E_{xc}$ $E_{\rm xc}$  implemented in Octopus for the moment. For closed-shell systems systems,  the necessary and sufficient conditions for the 1-RDM to be $N$-representable \cite{Coleman_1963}, i.e.\ \textit{i.e.}\  to correspond to a $N$-electron wavefunction, is that $ 0 \leq n_{i} \leq 2$ and \begin{equation}  \sum_{i=1}^{\infty}n_{i}=N.\label{eqsumocc}  \end{equation} 
...
\Omega\left[N,\{\vartheta_i\} ,\{\phi_i(\vec{r})\}\right]= 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)  \end{multline}  which has to be stationary with respect to variations in $\{\vartheta_i\}$, $\{\phi_i(\vec{r})\}$ and $\{\phi_i^{*}(\vec{r})\}$. In any practical calculation the infinite sums have to be truncated including only a finite number of occupation numbers and natural orbitals. However, since the occupation numbers $n_j$ decay very quickly for $j>N$ $j>N$,  this is not problematic. The variation of $\Omega$ is done in two steps: for a fixed set of orbitals 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 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. Due to the dependence on the occupation numbers, the natural orbital natural-orbital  minimization does not lead to an eigenvalue equation like in DFT or Hartree-Fock. The implementation of the natural orbital minimization follows the method by Piris and Ugalde~\cite{Piris}. Varying $\Omega$ with respect to the orbitals for fixed occupation numbers one obtains \begin{equation}  \lambda_{ji} = n_i\left\langle\phi_j\left|-\frac{\nabla^2}{2}+v_{ext}\right|\phi_i\right\rangle n_i\left\langle\phi_j\left|-\frac{\nabla^2}{2}+v_{\rm ext}\right|\phi_i\right\rangle  +\int \mathrm{d}\vec{r} \frac{\delta E_{Hxc}}{\delta E_{\rm Hxc}}{\delta  \phi_i^{*}(\vec{r})}\phi_j^{*}(\vec{r}). \end{equation}  At the extremum, the matrix of the Lagrange multipliers must be Hermitian, \textit{i.e.}   \begin{equation} 
...
\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 Heaviside function. We initialize the whole matrix as $F_{ji}=(\lambda_{ji}+\lambda_{ij}^{*})/2$. In every iteration we diagonalize $\mathbf{F}$ $\mathbf{F}$,  keeping the diagonal elements for the next iteration 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~\cite{Piris}. The criterion to exit the natural orbital natural-orbital  optimization is that the difference in the total energies calculated in two successive $\mathbf{F}$ diagonalizations is smaller than a threshold. 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}), $\mathbf{F}$  are close to zero. Asalready  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. Unfortunately, there are unbound states among the HF/DFT unoccupied states which are a bad starting point for the weakly occupied natural orbitals. When calculated in a finite grid these orbitals are essentially the eigenstates of a particle in a box. Using the exact-exchange approximation (EXX) in an optimized effective potential 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 self-interaction-free  for both occupied and unoccupied orbitals. Using HF or LDA orbitals to start a RDMFT calculation, the natural orbitals do not converge to any reasonable shape shape,  but even when starting from EXX one needs to further localize the unoccupied states. Thus, we have found that in order to improve the starting point for our calculation we can multiply each unoccupied orbital by a set of Gaussian functions centered at the positions of the atom. atoms.  As the unbound states are initially more delocalized than the bound ones ones,  we choose a larger exponent for them. 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 the Gaussian basis set Gaussian-basis-set  RDMFT code HIPPO~\cite{ML_Benchmark}. For the Octopus calculation calculation,  we kept 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 implementations of RDMFT \cite{GPB2005} \cite{GPB2005},  keeping the nice feature of not diverging strongly in the dissociation limit. However, for internuclear distances $R$ bigger than $1~a.u.$ 1 a.u.,  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, 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 support for open-shell systems and additional xc functionals.