# Variational minimization scheme for the one-particle reduced density matrix
functional theory in the ensemble N-representability domain

Matthieu Vladaj  ICGM, Université de Montpellier, CNRS, ENSCM, Montpellier,
France  Quentin Marécat  ICGM, Université de Montpellier, CNRS, ENSCM,
Montpellier, France  Bruno Senjean  ICGM, Université de Montpellier, CNRS,
ENSCM, Montpellier, France  Matthieu Saubanère<EMAIL_ADDRESS>ICGM, Université de Montpellier, CNRS, ENSCM, Montpellier, France  Univ.
Bordeaux, CNRS, LOMA, UMR 5798, F-33400 Talence, France

###### Abstract

The one-particle reduced density-matrix (1-RDM) functional theory is a
promising alternative to density-functional theory (DFT) that uses the 1-RDM
rather than the electronic density as a basic variable. However, long-standing
challenges such as the lack of Kohn–Sham scheme and the complexity of the pure
$N$-representability conditions are still impeding its wild utilization.
Fortunately, ensemble $N$-representability conditions derived in the natural
orbital basis are known and trivial, such that almost every functionals of the
1-RDM are actually natural orbital functionals which do not perform well for
all the correlation regimes. In this work, we propose a variational
minimization scheme in the ensemble $N$-representable domain that is not
restricted to the natural orbital representation of the 1-RDM. We show that
splitting the minimization into the diagonal and off-diagonal part of the
1-RDM can open the way toward the development of functionals of the orbital
occupations, which remains a challenge for the generalization of site-
occupation functional theory in chemistry. Our approach is tested on the
uniform Hubbard model using the Müller and the Töws–Pastor functionals, as
well as on the dihydrogen molecule using the Müller functional.

##  I Introduction

Over the last two decades, density functional theory (DFT) has played a major
role in the study of the physico-chemical properties of molecules and
materials. The key idea behind DFT, formulated by Hohenberg and Kohn Hohenberg
and Kohn (1964), is to replace the $N$-electron wavefunction by the electronic
density as the fundamental variable of the many-electron problem. Within the
Kohn–Sham (KS) decomposition scheme Kohn and Sham (1965), the electronic
repulsion effects are treated implicitly by a Hartree-exchange-correlation
(Hxc) functional of the density, thus leading to a good ratio between
computational cost and accuracy. However, despite the numerous Hxc functionals
developed in the past decades Marques  _et al._ (2012); Mardirossian and Head-
Gordon (2017), KS-DFT still faces issues when dealing with systems that
exhibit strong (or static) correlation effects Cramer and Truhlar (2009);
Cohen  _et al._ (2008).

Alternatively, replacing the density by the one-particle reduced density-
matrix (1-RDM) as the fundamental variable leads to the 1-RDM functional
theory (1-RDMFT) Löwdin (1955); Gilbert (1975); Coleman (1963); Haar (1961)
and has several advantages. First, kinetic and exchange energy contributions
are expressed exactly and analytically with respect to the 1-RDM. Second,
fractional occupations of the orbitals are naturally obtained within 1-RDMFT,
thus allowing for the proper description of bond dissociation where static
correlation effects are significant. Despite these advantages, 1-RDMFT does
not yet compete with DFT, partially due to the lack of an efficient KS scheme.
Indeed, the non-idempotent 1-RDM of an interacting system cannot be mapped to
its non-interacting equivalent that is idempotent by construction, as it
results from a single-electron Schrödinger equation Piris and Pernal (2017).
Besides, the proper minimization of the energy functional of the 1-RDM is
highly non-trivial as it requires maintaining the pure $N$-representability
conditions of the 1-RDM (also known as the generalized Pauli constraints
Altunbulak and Klyachko (2008); Schilling  _et al._ (2018); Castillo  _et al._
(2023)) which number increases significantly with the number of electrons and
orbitals Theophilou  _et al._ (2015a). Instead, the 1-RDM is usually written
in the natural orbital (NO) representation, and the minimization is performed
via a Langrangian-based system of differential equations where both the NOs
and their occupation numbers are optimized while satisfying the more trivial
ensemble $N$-representability conditions Coleman (1963, 1972); Valone (1980).
Because these conditions are known in the NO basis, Hxc functionals are almost
all NO-functionals (NOFs) Müller (1984); Csányi and Arias (2000); Csányi  _et
al._ (2002); Buijse and Baerends (2002); Sharma  _et al._ (2008); Lathiotakis
_et al._ (2009); Goedecker and Umrigar (1998); Holas (1999); Cioslowski and
Pernal (1999); Cioslowski  _et al._ (2023). The most recent and efficient ones
are the Piris natural orbitals functionals (PNOFs) Piris and Otto (2003);
Piris (2004); Leiva and Piris (2005); Piris (2006); Piris  _et al._ (2007);
Piris and Ugalde (2009); Piris  _et al._ (2010, 2011); Piris (2013); Piris and
Ugalde (2014); Lopez  _et al._ (2015) based on the cumulants Kutzelnigg and
Mukherjee (1999). Note that despite promising results, NOFs do not perform
well for all the correlation regimes Ramos-Cordoba  _et al._ (2015);
Cioslowski  _et al._ (2015); Mitxelena  _et al._ (2017); Rodriguez-Mayorga
_et al._ (2017); Mitxelena  _et al._ (2018), although recent progress are
worth mentioning Piris (2021).

Solving the aforementioned non-linear differential equations is
computationally demanding, so that alternative strategies have been developed
to define efficient and robust algorithms to perform the energy functional
minimization. Gilbert proposed to determine the NOs using a non-local
potential for functionals that are explicitly dependent on the 1-RDM Gilbert
(1975), whose existence was later shown by Pernal in the case of NOFs Pernal
(2005). In this approach, partially occupied orbitals are all degenerate. A
faster scheme for the NO optimization is based on a hermitian matrix built
from the Lagrange multiplier, as proposed by Piris and Ugalde Piris and Ugalde
(2009). One can also use the Hessian and couple the optimization of the NOs
and occupation numbers Castillo  _et al._ (2023). Another interesting scheme
called the local 1-RDMFT was developed by Lathiotakis and coworkers, where the
minimization of the NOs is restricted to a domain where the orbitals are
eigenfunctions of a single-electron Hamiltonian with a local potential, thus
bridging the 1-RDMFT and KS-DFT approaches Lathiotakis  _et al._ (2014);
Theophilou  _et al._ (2015b). Written in the NO representation, these methods
as well as methods based on an analogy with the temperature are particularly
suited for NOFs Goedecker and Colombo (1994); Baldsiefen and Gross (2013);
Saubanère  _et al._ (2016); Wang and Baerends (2022). Indeed, in that case the
energy does only depend explicitly on the occupation numbers while the
dependence on the NOs is implicitly accounted to in the Hartree and exchange
integrals. Alternatively, minimization schemes based on a different
representation than natural orbitals were also proposed, for instance using
projected gradient-based algorithms constrained to preserve ensemble
$N$-representability conditions Cancès and Pernal (2008).

In this work, we introduce a new variational minimization scheme so-called
density-matrix interpolation variational ansatz (DIVA) which relies on
remarkable topological properties of the ensemble $N$-representable set of
1-RDM. In contrast to the usual minimization schemes adapted to NOFs, DIVA
aims to minimize the energy functional of the 1-RDM expressed in any
orthonormal orbital basis, thus motivating the development of new 1-RDM
functionals beyond NOFs. By separating the minimization over the diagonal and
off-diagonal part of the 1-RDM, one can also make a connection between 1-RDMFT
and DFT in a lattice, also called site-occupation functional theory (SOFT)
Gunnarsson and Schönhammer (1986); Schönhammer  _et al._ (1995); Lima  _et
al._ (2002, 2003). We show how this separation can be used to estimate Hartree
Hxc potentials, thus paving the way toward the development of functionals of
the orbital occupations which is the current limitation for the implementation
of SOFT in chemistry Fromager (2015); Coe (2019).

The paper is organized as follows. In Sec. II, we give a brief introduction to
1-RDMFT. The separation of the minimization process into that bridges 1-RDMFT
and SOFT is then presented in Sec. III. Our variational minimization strategy
and its various implementations are described in Sec. IV, and tested on the
Hubbard model and the hydrogen molecule in Sec. V. Finally, conclusions and
perspectives are provided in Sec. VI. Note that in the following, we will use
‘1-RDM’ and ‘density matrix’ to refer to the same object, i.e. the one-
particle reduced density matrix.

##  II Reduced density-matrix functional theory

Let us consider the second-quantized electronic Hamiltonian written in an
orthonormal basis set, referred to as the computational basis set and composed
of $N$ spin-orbitals $\\{\varphi_{i}({\bf{x}})\\}$ that are different from the
set of natural spin-orbitals $\\{\phi_{i}({\bf{x}})\\}$,

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle\hat{T}+\hat{U}$ (1)
$\displaystyle=$
$\displaystyle\sum_{ij}^{N}t_{ij}\hat{c}_{i}^{\dagger}\hat{c}_{j}+\dfrac{1}{2}\sum_{ijkl}^{N}U_{ijkl}\hat{c}_{i}^{\dagger}\hat{c}_{j}^{\dagger}\hat{c}_{k}\hat{c}_{l},$

where $\hat{c}_{i}^{\dagger}$ ($\hat{c}_{i}$) corresponds to the creation
(annihilation) of an electron in $\varphi_{i}({\bf{x}})$. The one-body
integrals

$\displaystyle t_{ij}=\int{\rm
d}{\bf{x}}\varphi_{i}^{*}({\bf{x}})\left(-\dfrac{\nabla^{2}}{2}+\sum_{I}\dfrac{Z_{I}}{|\mathbf{r}-{\bf
R}_{I}|}\right)\varphi_{j}({\bf{x}})$ (2)

contain the kinetic and external potential contributions, while the two-body
integrals

$\displaystyle U_{ijkl}=\int{\rm d}{\bf{x}}_{1}{\rm
d}{\bf{x}}_{2}\dfrac{\varphi_{i}^{*}({\bf{x}}_{1})\varphi_{j}^{*}({\bf{x}}_{2})\varphi_{k}({\bf{x}}_{2})\varphi_{l}({\bf{x}}_{1})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}$
(3)

contain contributions from the electron-electron interactions. The summation
in Eq. (2) runs over the number of atoms at positions ${\mathbf{R}_{I}}$ and
${\bf{x}}_{i}=(\mathbf{r}_{i},\sigma_{i})$ corresponds to the position
$\mathbf{r}_{i}$ and spin $\sigma_{i}$ of the electron $i$.

In 1-RDMFT, the basic variable ${\bm{\gamma}}$ is the 1-RDM given by

$\displaystyle\gamma_{ij}=\langle\Psi|\hat{c}_{i}^{\dagger}\hat{c}_{j}|\Psi\rangle,$
(4)

with $|\Psi\rangle$ the ground-state wavefunction of $\hat{H}$. The energy can
in principle be written as a functional of the 1-RDM,

$\displaystyle E[{\bm{\gamma}}]=T[{\bm{\gamma}}]+W[{\bm{\gamma}}],$ (5)

where

$\displaystyle T[{\bm{\gamma}}]=\sum_{ij}t_{ij}\gamma_{ij}$ (6)

is an explicit functional of the 1-RDM and $W[{\bm{\gamma}}]$ is the energy
functional describing the electronic repulsion. The analytical form of this
functional remains unknown and is usually approximated, though it can be
written using the Levy–Lieb formalism as follows,

$\displaystyle
W[{\bm{\gamma}}]=\min_{\Psi\rightarrow{\bm{\gamma}}}\bra{\Psi}\hat{U}\ket{\Psi},$
(7)

where the notation $\Psi\rightarrow{\bm{\gamma}}$ refers to any wavefunction
giving the one-particle density matrix ${\bm{\gamma}}$. Within the variational
principle of 1-RDMFT, the ground-state energy $E_{0}$ can be obtained as a
saddle point,

$\displaystyle E_{\rm 0}=\min_{{\bm{\gamma}}}E[{\bm{\gamma}}].$ (8)

More precisely, as a corollary of the Hohenberg–Kohn theorem (later extended
by Gilbert Gilbert (1975) to non-local potential), there exists a one-to-one
mapping between the ground-state wavefunction and the 1-RDM, thus allowing to
determine the ground-state energy variationally. The conditions for
${\bm{\gamma}}$ to derive from a $N_{e}$-electron wavefunction – so-called
$N$-representability conditions – are highly non-trivial Altunbulak and
Klyachko (2008); Schilling  _et al._ (2018); Castillo  _et al._ (2023) and
cannot be used in practice to define the domain in which the minimization of
Eq. (8) takes place. However, the conditions for ${\bm{\gamma}}$ to be
ensemble $N$-representable are known and straightforward Coleman (1963) and
are given as follows,

$\begin{split}\sum_{k=1}^{N}\eta_{k}=N_{e},\;\;0\leq\eta_{k}\leq
1,\end{split}$ (9)

where $\eta_{k}$ is the occupation number of the $k$-th natural orbital.
Hence, it is common to use the natural orbital representation of the density
matrix, i.e. ${\bm{\gamma}}\equiv(\\{\eta_{k}\\},\\{\ket{\phi_{k}}\\})$. As
originally proposed by Gilbert Gilbert (1975), these conditions are fulfilled
during the minimization by introducing a Lagrangian function in the ensemble
$N$-representable domain,

$\displaystyle\mathcal{E}$ $\displaystyle=$ $\displaystyle
E[\\{\eta_{k}\\},\\{\ket{\phi_{k}}\\}]-\mu\bigg{(}\sum_{k}\cos^{2}(\theta_{k})-N_{e}\bigg{)}$
(10)
$\displaystyle-\sum_{jk}\lambda_{jk}\bigg{(}\langle\phi_{j}|\phi_{k}\rangle-\delta_{jk}\bigg{)},$

where $\mu$ and $\\{\lambda_{jk}\\}$ are Lagrange multipliers, the first term
($\eta_{k}=\cos^{2}(\theta_{k})$) fulfills Eq. (9) and the second term the
orthonormality of the natural orbitals $\\{|\phi_{k}\rangle\\}$. Minimizing
Eq. (10) with respect to the whole set of occupation numbers and NOs leads to
two non-linear differential equations,

$\displaystyle\dfrac{\partial\mathcal{E}}{\partial\theta_{k}}=\left(\mu-\dfrac{\partial
E}{\partial\eta_{k}}\right)\sin(2\theta_{k})=0,$ (11)

and

$\displaystyle\displaystyle\dfrac{\partial\mathcal{E}}{\partial\langle\phi_{k}|}=\dfrac{\partial
E}{\partial\langle\phi_{k}|}-\sum_{j}\lambda_{jk}|\phi_{j}\rangle=0,$ (12)

that are solved successively but with slow convergence Cartier and Giesbertz
(2024). This situation contrasts with the variational equations in KS-DFT and
SOFT that are solved efficiently through the definition of a single-particle
Schrödinger equation Kohn and Sham (1965); Gunnarsson and Schönhammer (1986).

##  III Link between 1-RDMFT and SOFT

Rather than decomposing the 1-RDM into natural orbitals and occupation
numbers, we decompose it in terms of its diagonal elements, i.e. orbital (or
site) occupations, ${\bf{n}}\equiv\\{\gamma_{ii}\\}_{i}$ and the $N(N-1)$ off-
diagonal elements denoted by ${\bm{\gamma}}^{\neq}$ such that
${\bm{\gamma}}={\bf{n}}\cdot\mathbb{1}+{\bm{\gamma}}^{\neq}$, where
$\mathbb{1}$ refers to the identity matrix. The energy functional in Eq. (5)
therefore reads

$E[{\bm{\gamma}}]=T[{\bf{n}},{\bm{\gamma}}^{\neq}]+W[{\bf{n}},{\bm{\gamma}}^{\neq}].$
(13)

Following the KS decomposition, Eq. (13) can be alternatively written as

$\displaystyle
E[{\bm{\gamma}}]=T_{s}[{\bf{n}}]+F[{\bf{n}},{\bm{\gamma}}^{\neq}],$ (14)

where $T_{s}[{\bf{n}}]$ represents the one-body (kinetic and potential) energy
contributions of the non-interacting system with orbital occupations
${\bf{n}}$ and $F[{\bf{n}},{\bm{\gamma}}^{\neq}]$ is a functional containing
all non trivial terms defined by

$\displaystyle
F[{\bf{n}},{\bm{\gamma}}^{\neq}]=T[{\bf{n}},{\bm{\gamma}}^{\neq}]-T_{s}[{\bf{n}}]+W[{\bf{n}},{\bm{\gamma}}^{\neq}].$
(15)

In analogy with the natural orbital representation, the minimization of the
energy leads to two saddle-point equations. The first saddle-point equation
reads

$\displaystyle\delta{\bf{n}}\left\\{\frac{\delta
T_{s}[{\bf{n}}]}{\delta{\bf{n}}}+\left.\frac{\delta
F[{\bf{n}},{\bm{\gamma}}^{\neq}]}{\delta{\bf{n}}}\right|_{{\bm{\gamma}}^{\neq}}\right\\}=0,$
(16)

and leads to a KS-like equation:

$\displaystyle\hat{h}_{\rm KS}[{\bf{n}}]|\Phi^{\rm
KS}\rangle=\left(\hat{T}+\hat{V}_{{\rm Hxc}}[{\bf{n}}]\right)|\Phi^{\rm
KS}\rangle=\mathcal{E}^{\rm KS}|\Phi^{\rm KS}\rangle$ (17)

with $V_{{\rm Hxc}}[{\bf{n}}]:=\left.V_{\rm
Hxc}[{\bf{n}}]\right|_{{\bm{\gamma}}^{\neq}}=\left.\delta
F[{\bf{n}},{\bm{\gamma}}^{\neq}]/\delta{\bf{n}}\right|_{{\bm{\gamma}}^{\neq}}$
being the Hxc potential defined at the saddle point of ${\bm{\gamma}}^{\neq}$
and $\hat{h}_{\rm KS}[{\bf{n}}]$ is the KS single-particle Hamiltonian. Eq.
(17) is determined concomitantly with the other saddle point equation:

$\displaystyle\delta{\bm{\gamma}}^{\neq}\left\\{\left.\frac{\delta
F[{\bf{n}},{\bm{\gamma}}^{\neq}]}{\delta{\bm{\gamma}}^{\neq}}\right|_{{\bf{n}}={\bf{n}}^{\rm
KS}}\right\\}=0,$ (18)

where ${\bf{n}}^{\rm KS}$ are the orbital occupations resulting from Eq. (17).
In practice, solving Eq. (17) appears numerically easy when $V_{{\rm
Hxc}}[{\bf{n}}]$ is known, however solving Eq. (18) remains difficult, in
particular due to the non-trivial representability conditions that
${\bm{\gamma}}^{\neq}$ has to fulfill.

Interestingly, one can also use a Levy–Lieb formalism to derive similar
equations,

$\displaystyle E_{\rm 0}$ $\displaystyle=$
$\displaystyle\min_{{\bf{n}}}\left\\{\min_{{\bm{\gamma}}\rightarrow{\bf{n}}}\left\\{E[{\bf{n}},{\bm{\gamma}}^{\neq}]\right\\}\right\\}$
(19) $\displaystyle=$
$\displaystyle\min_{{\bf{n}}}\left\\{\min_{{\bm{\gamma}}\rightarrow{\bf{n}}}\left\\{T_{s}[{\bf{n}}]+F[{\bf{n}},{\bm{\gamma}}^{\neq}]\right\\}\right\\}$
$\displaystyle=$
$\displaystyle\min_{{\bf{n}}}\left\\{T_{s}[{\bf{n}}]+\min_{{\bm{\gamma}}\rightarrow{\bf{n}}}\left\\{F[{\bf{n}},{\bm{\gamma}}^{\neq}]\right\\}\right\\},$

where ${\bm{\gamma}}\rightarrow{\bf{n}}$ corresponds to density matrices with
orbital occupations ${\bf{n}}$. In that framework, the Hxc potential is
defined as

$\displaystyle V_{{\rm
Hxc}}[{\bf{n}}]=\frac{\delta}{\delta{\bf{n}}}\min_{{\bm{\gamma}}\rightarrow{\bf{n}}}\left\\{F[{\bf{n}},{\bm{\gamma}}^{\neq}]\right\\},$
(20)

thus bridging 1-RDMFT and SOFT, where the Hxc potential functional of the
orbital occupations can be obtained after solving Eq. (18) for a fixed orbital
occupation ${\bf{n}}$.

##  IV Density-matrix Interpolation Variational ansatz

###  IV.1 Topological properties of the ensemble $N$-representable domain of
the 1-RDM

Let us consider the space $\mathcal{M}_{N}(\mathbb{R})$ of $N\times N$
symmetric and real matrices, $N$ being the number of spin-orbitals in the
system. In this section we denote by $\Omega_{\bm{\gamma}}$ the domain of
ensemble $N$-representable density matrices that is included in
$\mathcal{M}_{N}(\mathbb{R})$,
$\Omega_{\bm{\gamma}}\subset\mathcal{M}_{N}(\mathbb{R})$. A significant
advantage of working with natural orbitals is that the ensemble
$N$-representable conditions are straightforward, i.e. a matrix
${\bm{\gamma}}\in\mathcal{M}_{N}(\mathbb{R})$ is ensemble $N$-representable if
all the occupation numbers $\eta_{k}$ fulfill $0\leq\eta_{k}\leq 1$.
Similarly, it means that ${\bm{\gamma}}$ is ensemble $N$-representable if both
${\bm{\gamma}}$ and $\mathbb{1}-{\bm{\gamma}}$ are positive semi-definite,
i.e. both showing positive or null eigenvalues $\eta_{k}\geq 0$ and
$1-\eta_{k}\geq 0$ for all $1\leq k\leq N$. Hence, the domain of ensemble
$N$-representable density matrices reads

$\displaystyle\Omega_{{\bm{\gamma}}}=\left\\{{\bm{\gamma}}\in\mathcal{M}_{N}(\mathbb{R})\;|\;\eta_{k}\geq
0\;\&\;1-\eta_{k}\geq 0,\;\forall k\right\\}.$ (21)

Interesting topological properties for the ensemble $N$-representable 1-RDM
domain $\Omega_{{\bm{\gamma}}}$ follow from these definitions. The first
remarkable property of $\Omega_{{\bm{\gamma}}}$ is the convexity, see proofs
in Appendix A. It follows that given a set of ensemble $N$-representable
matrices $\left\\{{\bm{\gamma}}^{(i)}\in\Omega_{{\bm{\gamma}}}\right\\}$, all
matrices constructed as

${\bm{\gamma}}({\bf{z}})=\sum_{i}z_{i}{\bm{\gamma}}^{(i)}$ (22)

also belong to $\Omega_{{\bm{\gamma}}}$ if all $z_{i}\geq 0$ and
$\sum_{i}z_{i}=1$.

Let us now define $\mathring{\Omega}_{{\bm{\gamma}}}$ the ensemble of matrices
${\bm{\gamma}}$ such that both ${\bm{\gamma}}$ and $\mathbb{1}-{\bm{\gamma}}$
are positive definite (in contrast to Eq. (21) where they are positive semi-
definite), i.e.

$\displaystyle\mathring{\Omega}_{{\bm{\gamma}}}=\left\\{{\bm{\gamma}}\in\mathcal{M}_{N}(\mathbb{R})\;|\;\eta_{k}>0\;\&\;1-\eta_{k}>0,\;\forall
k\right\\}.$ (23)

The convexity of $\mathring{\Omega}_{{\bm{\gamma}}}$ follows the same proof as
in Appendix A for $\Omega_{{\bm{\gamma}}}$. Besides, we introduce

$\displaystyle\delta\Omega_{{\bm{\gamma}}}=\Omega_{{\bm{\gamma}}}\,\backslash\,\mathring{\Omega}_{{\bm{\gamma}}},$
(24)

such that $\delta\Omega_{{\bm{\gamma}}}$ contains the density matrices
${\bm{\gamma}}$ that exhibit at least one integer eigenvalue
$\eta_{k}=\\{0,1\\}$. Contrary to $\Omega_{{\bm{\gamma}}}$ and
$\mathring{\Omega}_{{\bm{\gamma}}}$, $\delta\Omega_{{\bm{\gamma}}}$ is not
convex. We show in Appendix B that $\mathring{\Omega}_{{\bm{\gamma}}}$
corresponds to the topological interior of $\Omega_{{\bm{\gamma}}}$,
$\delta\Omega_{{\bm{\gamma}}}$ is the topological boundary of
$\Omega_{{\bm{\gamma}}}$ and
$\Omega_{{\bm{\gamma}}}=\mathring{\Omega}_{{\bm{\gamma}}}\cup\delta\Omega_{{\bm{\gamma}}}$
is a topological closed set.

Finally, let us define
$\delta\Omega_{{\bm{\gamma}}}^{*}\subset\delta\Omega_{{\bm{\gamma}}}$ the set
of idempotent matrices, i.e. ${\bm{\gamma}}={{\bm{\gamma}}}^{2}$ if
${\bm{\gamma}}\in\delta\Omega_{{\bm{\gamma}}}^{*}$ and

$\displaystyle\delta\Omega_{{\bm{\gamma}}}^{*}=\left\\{{\bm{\gamma}}\in\Omega_{{\bm{\gamma}}}\,|\,\eta_{k}=\left\\{0,1\right\\},\,\forall
k\right\\},$ (25)

which is also non convex. $\delta\Omega_{{\bm{\gamma}}}^{*}$ corresponds to
the set of extreme points of $\Omega_{{\bm{\gamma}}}$, i.e. given
${\bm{\gamma}}_{*}\in\delta\Omega_{{\bm{\gamma}}}^{*}$ there is no
${\bm{\gamma}}^{A}$ and ${\bm{\gamma}}^{B}$ in $\Omega_{{\bm{\gamma}}}$ such
that ${\bm{\gamma}}_{*}=z{\bm{\gamma}}^{A}+(1-z){\bm{\gamma}}^{B}$ with
$0<z<1$. Following the Krein–Milman theorem, there exists a decomposition in
terms of idempotent matrices, i.e. for any
${\bm{\gamma}}\in\Omega_{\bm{\gamma}}$,

${\bm{\gamma}}({\bf{z}})=\sum_{i}z_{i}{\bm{\gamma}}^{(i)}_{*},$ (26)

with ${\bf{z}}=\\{z_{i}\\}$, $\sum_{i}z_{i}=1$ and
${\bm{\gamma}}^{(i)}_{*}\in\delta\Omega_{{\bm{\gamma}}}^{*}$. Given Eqs. (8)
and (26), it follows the variational ansatz:

$E_{\rm 0}=\min_{{\bf{z}}}E[{\bm{\gamma}}({\bf{z}})].$ (27)

Hence, the energy is minimized with respect to the weights ${\bf{z}}$ defining
the interpolation of ${\bm{\gamma}}({\bf{z}})$ in terms of idempotent 1-RDM.
We show in Appendix C that for each
${\bm{\gamma}}\in\mathring{\Omega}_{{\bm{\gamma}}}$, the sum in Eq. (26) runs
up to $2^{N}$ idempotent matrices. It follows that the number of variational
parameters rapidly becomes prohibitive when the size of the system increases.
To reduce the number of parameters, one could define the trial set as every
Slater determinants corresponding to the singly- and doubly-excited
determinants with respect to the Hartree–Fock state. The efficiency and
accuracy of such ansatz strongly depends on the size and relevance of the
trial set of $\\{{\bm{\gamma}}_{*}^{(i)}\\}$.

Rather than interpolating ${\bm{\gamma}}$ in terms of idempotent 1-RDM, we can
use the property that for each ${\bm{\gamma}}\in\Omega_{{\bm{\gamma}}}$ there
exists a decomposition in terms of boundary matrices belonging to
$\delta\Omega_{\bm{\gamma}}$ [see Eq. (24)]. Then, Eq. (26) turns to

${\bm{\gamma}}({\bf{z}})=\sum_{i}z_{i}{\bm{\gamma}}_{B}^{(i)}.$ (28)

As ${\bm{\gamma}}_{B}^{(i)}$ can be seen as a linear combination of
${\bm{\gamma}}_{*}^{(i)}$, one can expect the amount of variational parameters
$\\{z_{i}\\}$ to be drastically reduced when too much
$\\{{\bm{\gamma}}_{*}^{(i)}\\}$ would be needed to reproduce ${\bm{\gamma}}$.
Indeed, any non-idempotent ${\bm{\gamma}}\in\mathring{\Omega}_{\bm{\gamma}}$
can in principle be obtained by interpolating only two boundary matrices.
Hence, while the set of idempotent matrices
$\left\\{{\bm{\gamma}}^{(i)}_{*}\right\\}$ in Eq. (26) is very large but
straightforward to construct, only a small relevant set of boundary matrices
$\left\\{{\bm{\gamma}}_{B}^{(i)}\right\\}$ would suffice but is non-trivial to
find. An efficient algorithm would consist in constructing this non-trivial
set of relevant $\left\\{{\bm{\gamma}}_{B}^{(i)}\right\\}$ on which the
minimization in Eq. (27) is performed. Eqs. (27) and (28) constitute the
density-matrix interpolation variational ansatz (DIVA). Note that by
construction, all density matrices of the form given by Eq. (28) are ensemble
$N$-representable.

Interestingly, the topological properties presented in
$\mathcal{M}_{N}(\mathbb{R})$ also hold if we restrain the considered metric
space to $\mathcal{M}^{N_{e}}_{N}(\mathbb{R})$ and
$\mathcal{M}^{{\bf{n}}}_{N}(\mathbb{R})$, corresponding to the spaces of
$N\times N$ symmetric and real matrices with a fixed trace (fixed number of
electron $N_{e}$) and fixed diagonal part ${\bf{n}}$, respectively.
Considering these restrictions should automatically reduce the number of trial
1-RDMs in Eq. (28).

###  IV.2 Iterative gradient-based implementations of DIVA

Let us now turn to the numerical implementations of DIVA. Starting from a
guess ${\bm{\gamma}}^{(0)}$ that already belongs to $\Omega_{{\bm{\gamma}}}$,
we want to construct iteratively a minimal trial set of boundary density
matrices to minimize the 1-RDMFT energy [Eq. (27)] according to Eq. (28),

${\bm{\gamma}}({\bf{z}})=z_{0}{\bm{\gamma}}^{(0)}+\sum_{i}z_{i}{\bm{\gamma}}_{B}^{(i)}.$
(29)

In the following, we highlight two ways of minimizing the energy functional
within DIVA. The first one consists in minimizing the energy functional with
respect to only one variational parameter $z_{i}$ at each iteration $i$ and is
called the mono-parameter DIVA. For the second one, the minimization is done
over the whole set of parameters $\\{z_{j}\\}_{j\leq i}$ at iteration $i$, and
is called the multi-parameter DIVA. These approaches are pictured in Figs. 1
and 2 and detailed below.

Figure 1: Illustration of the mono-parameter DIVA. The set
$\Omega_{{\bm{\gamma}}}$ is illustrated by the union of
$\delta\Omega_{{\bm{\gamma}}}$ (blue circle) and the topological interior. Red
arrows represent the direction of the energy gradient, used to find the
${\bm{\gamma}}_{B}^{(i)}$ belonging to the border of the representability
domain ($\delta\Omega_{\bm{\gamma}}$). Each ${\bm{\gamma}}^{(i)}$ is obtained
as an interpolation between ${\bm{\gamma}}^{(i-1)}$ and
${\bm{\gamma}}_{B}^{(i-1)}$. The star represents the ground-state 1-RDM
${\bm{\gamma}}_{0}$ (at the level of the functional approximation) that is
recovered after convergence.

Figure 2: Illustration of the multi-parameter DIVA. The set
$\Omega_{{\bm{\gamma}}}$ is illustrated by the union of
$\delta\Omega_{{\bm{\gamma}}}$ (blue circle) and the topological interior. Red
arrows represent the direction of the energy gradient, used to find the
${\bm{\gamma}}_{B}^{(i)}$ belonging to the border of the representability
domain ($\delta\Omega_{\bm{\gamma}}$). Each ${\bm{\gamma}}^{(i>2)}$ is a
linear combination of all ${\bm{\gamma}}_{B}^{(j)}$ with $j<i$. The star
represents the ground-state 1-RDM ${\bm{\gamma}}_{0}$ (at the level of the
functional approximation) that is recovered after convergence.

####  IV.2.1 Mono-parameter DIVA

The mono-parameter DIVA is implemented according to the following steps.

  1. 1.

Initialization (iteration $i=0$): set ${\bm{\gamma}}^{(0)}$ to be the
(idempotent) density matrix coming from an Hartree–Fock or KS-DFT calculation,
thus ensuring that ${\bm{\gamma}}^{(0)}$ belongs to $\Omega_{{\bm{\gamma}}}$,
and even to
$\delta\Omega_{\bm{\gamma}}^{*}\subset\delta\Omega_{\bm{\gamma}}\subset\Omega_{\bm{\gamma}}$.
To avoid divergence of the energy gradient with respect to
${\bm{\gamma}}^{(0)}$ having occupation numbers equal to 0, which happens for
many NOFs, the first gradient is set to
$\nabla_{ij}[{\bm{\gamma}}]=(1-\delta_{ij})\gamma_{ij}^{(0)}$.

  2. 2.

Search the boundary density matrix
${\bm{\gamma}}_{B}^{(i)}\in\delta\Omega_{{\bm{\gamma}}}$ following the
opposite direction of the energy gradient. This is done in two steps:

    1. (a)

search for ${\bm{\gamma}}_{\rm
out}={\bm{\gamma}}^{(i)}-\theta{\bm{\nabla}}[{\bm{\gamma}}^{(i)}]$ by
incrementing $\theta$ until ${\bm{\gamma}}_{\rm
out}\notin\Omega_{{\bm{\gamma}}}$

    2. (b)

Once ${\bm{\gamma}}_{\rm out}$ is found, vary $\theta$ using a bracketing
method until ${\bm{\gamma}}_{\rm out}\in\delta\Omega_{{\bm{\gamma}}}$ and set
${\bm{\gamma}}_{B}^{(i)}={\bm{\gamma}}_{\rm out}$

  3. 3.

start a new iteration, $i=i+1$

  4. 4.

Interpolate between ${\bm{\gamma}}^{(i-1)}$ and ${\bm{\gamma}}_{B}^{(i-1)}$ to
find
${\bm{\gamma}}^{(i)}=(1-z_{i}){\bm{\gamma}}^{(i-1)}+z_{i}{\bm{\gamma}}_{B}^{(i-1)}$
such that

$\dfrac{\partial E[{\bm{\gamma}}^{(i)}]}{\partial z_{i}}=0.$ (30)

  5. 5.

Compute the energy gradient ${\bm{\nabla}}[{\bm{\gamma}}^{(i)}]$. Note that
the number of electrons is conserved by setting the diagonal part
$\nabla_{ii}[{\bm{\gamma}}^{(i)}]=\nabla_{ii}[{\bm{\gamma}}^{(i)}]-\mu^{(i)}$
with $\mu^{(i)}=(1/N)\sum_{i}\nabla_{ii}$. If the orbital occupations should
be kept unchanged, set $\nabla_{ii}[{\bm{\gamma}}^{(i)}]=0$, $\forall i$.

  6. 6.

If $\delta E=|E^{(i-1)}-E^{(i)}|<\upsilon$ and
$\delta{\bm{\gamma}}=|{\bm{\gamma}}^{(i-1)}-{\bm{\gamma}}^{(i)}|<\varepsilon$
where $\upsilon$ and $\varepsilon$ are the convergence criteria, set
${\bm{\gamma}}^{\rm opt}={\bm{\gamma}}^{(i)}$, which should be equal to the
ground-state 1-RDM ${\bm{\gamma}}_{0}$ corresponding to the level of
approximation of the chosen 1-RDM functional (the star in Fig. 1). Else, go to
step 2.

At convergence, the density matrix fulfills Eq. (28),
${\bm{\gamma}}({\bf{z}})=\sum_{i}^{M}Z_{i}{\bm{\gamma}}_{B}^{(i)}$, with
$Z_{i}=z_{i+1}\prod_{i+2}^{M}(1-z_{i})$ and $M$ the number of iterations
needed to achieve convergence. Following Ref. Cancès and Pernal (2008), step 2
of the algorithm could be made more efficient by projecting
${\bm{\gamma}}^{(i)}$ to the boundary domain $\delta\Omega_{\bm{\gamma}}$.
However, it would lead to a ${\bm{\gamma}}_{B}^{(i)}$ that deviates from the
direction pointed by the gradient. Also, the opposite energy gradient
direction can be replaced by the conjugate gradient. This is what we
implemented in this work, as we observed a much faster convergence using it.

####  IV.2.2 Multi-parameter DIVA

Apart from step 4, the multi-parameter approach remains the same as the mono-
parameter one. Rather than doing the interpolation between
${\bm{\gamma}}^{(i-1)}$ and ${\bm{\gamma}}_{B}^{(i-1)}$ only, the
interpolation is done over the whole set of previously determined
${\bm{\gamma}}_{B}^{(j)}$ with $j<i$:

${\bm{\gamma}}^{(i)}=z^{(0)}{\bm{\gamma}}^{(0)}+\sum_{j=1}^{i-1}z_{j}{\bm{\gamma}}_{B}^{(j)}.$
(31)

The minimization is done over the whole set of parameters
${\bf{z}}=(z_{0},\dots,z_{(i-1)})$,

$\dfrac{\partial E[{\bm{\gamma}}^{(i)}]}{\partial{\bf{z}}}=0.$ (32)

In contrast to the previous implementation, we expect to achieve a much faster
convergence at the expense of more expensive optimizations of the ${\bf{z}}$
parameters.

####  IV.2.3 SOFT-DIVA

According to Sec. III, by splitting the variational minimization of the energy
into successive minimizations over the orbital occupations and the off-
diagonal elements of the 1-RDM, one can extract the Hxc potential, functional
of the orbital occupations used in SOFT [see Eq. (20)]. Everytime the second
saddle point equation in Eq. (18) is fulfilled, the Hxc potential of Eq. (17)
is computed. In practice Eq. (18) is solved via the DIVA algorithm while
setting the diagonal part of the energy gradient to zero, such that the
orbital occupations are fixed (see step 5 in Sec. IV.2.1). At convergence of
the DIVA algorithm, the Hxc potential can be extracted by aligning the
chemical potentials of the interacting and non-interacting (KS) systems as
follows

$\displaystyle V_{\rm Hxc}[{\bf{n}}]$ $\displaystyle=$
$\displaystyle\left.\dfrac{\delta
F[{\bf{n}},{\bm{\gamma}}^{\neq}]}{\delta{\bf{n}}}\right|_{{\bm{\gamma}}^{\neq}}$
(33) $\displaystyle=$ $\displaystyle\left.\dfrac{\delta
E[{\bf{n}},{\bm{\gamma}}^{\neq}]}{\delta{\bf{n}}}\right|_{{\bm{\gamma}}^{\neq}}-\dfrac{\delta
T_{s}[{\bf{n}}]}{\delta{\bf{n}}},$ $\displaystyle=$
$\displaystyle\left.\dfrac{\delta
E[{\bf{n}},{\bm{\gamma}}^{\neq}]}{\delta{\bf{n}}}\right|_{{\bm{\gamma}}^{\neq}}-\varepsilon_{\rm
H}[{\bf{n}}],$

where $\varepsilon_{\rm H}[{\bf{n}}]$ is the HOMO orbital energy of the non-
interacting system with fixed orbital occupations ${\bf{n}}$. Hence, the Hxc
potential is updated every time the minimizing off-diagonal elements of the
1-RDM are found. This Hxc potential is then plugged into the single-particle
Hamiltonian (or Fock matrix), thus generating a new idempotent 1-RDM with a
new set of orbital occupations. The off-diagonal elements will then be
optimized again using DIVA, and so on. Upon convergence, one obtains a
numerical Hxc potential, functional of the orbital occupations, with a level
of theory corresponding to the 1-RDM functional approximation used for the
calculation of the energy.

##  V Results and discussions

We propose to study the performance of DIVA on the one-dimensional (1D)
Hubbard model and on a simple molecular system. For the 1D Hubbard model, in
the spirit of Mitxelenal et al. Mitxelena  _et al._ (2017), we compare a NOF
(the Müller functional) with a functional derived in lattice representation
such as proposed by Töws and Pastor to approximate the correlation energy Töws
and Pastor (2011, 2012); Töws  _et al._ (2014). For molecular systems, only
NOFs are available and we also employ the Müller functional.

###  V.1 Hubbard model

The homogeneous Hubbard model can be defined by the following Hamiltonian,

$\hat{H}=-t\sum_{\langle
ij\rangle,\sigma}\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma}+\frac{U}{2}\sum_{i}\sum_{\sigma\sigma^{\prime}}\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma^{\prime}}^{\dagger}\hat{c}_{i\sigma^{\prime}}\hat{c}_{i\sigma}$
(34)

Figure 3:  Ground-state energy of the 1D Hubbard model as a function of the
electron filling $n$ at different values of the Coulomb integral $U/t$.
Results are given for the exact Bethe–Ansatz (gray dots), and using DIVA with
the Töws–Pastor functional (full orange lines) and the Müller functional
(dashed blue lines).

where the notation $\langle\cdot\rangle$ refers to first neighbor lattice
sites. In contrast to Eq. (1), the indices correspond to spatial orbitals (or
sites) and $\sigma=\\{\uparrow,\downarrow\\}$ denotes the spin. Focusing on
the homogeneous model allows us to bypass Eq. (17) and to solely focus on the
solution of Eq. (18) within the DIVA framework. In practice, chains with 202
sites have been considered with periodic boundary conditions, and the
bisection algorithm has been employed for bracketing and line minimization in
steps 2b and 4, respectively (see Sec. IV.2.1). We have considered the Müller
functional, defined as a function of the natural orbital occupation number
(Blöch function for periodic systems), and the Töws–Pastor functional that is
defined in the lattice site representation. Note that the Müller functional
has been investigated in the context of finite Hubbard chains Kamil  _et al._
(2016); Mitxelena  _et al._ (2017), and the Töws–Pastor functional has been
applied in the context of the Anderson model Töws and Pastor (2011, 2012);
Töws  _et al._ (2014). In terms of performance, we found that fast convergence
is achieved with DIVA (5 iteration max to reach an energy difference of
$\delta E/t<10^{-7}$) for both functionals. However, computational time is
smaller for the Töws–Pastor functional compared to the Müller one. Indeed, for
the Müller functional depending only on natural orbital occupation numbers, a
standard minimization algorithm in the natural orbital representation appears
well-suited Piris and Ugalde (2009), while a minimization scheme directly in
the lattice representation such as DIVA appears more adapted for the
Töws–Pastor functional. For the latter, the computation of the energy and
energy gradients (needed several time at each minimization step) do not
require to switch between lattice and natural orbital representations.

In Fig. 3, the ground-state energy is shown as a function of the electron
filling $n$. It shows that the Töws–Pastor functional is in fair agreement
with the exact Bethe–Ansatz (BA) solution Lieb and Wu (1968); Shiba (1972)
regardless of the electron filling and the strength of the Coulomb repulsion
$U/t$. On the contrary, the Müller functional leads to inaccurate energies
away from the half-band filling and especially for large $U/t$. Note that this
wrong behaviour goes together with an unphysical negative value of $\partial
E/\partial U$ corresponding to the number of double occupations (not shown).
These unphysical results for the Müller functional arises from equivalent spin
on-site exchange contributions to the energy, associated with the operator
$\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}\hat{c}_{i\sigma}$
that is supposed to be nil.

Figure 4:  Occupation numbers as a function of the Blöch wave number $k$ for
the half-band filled 1D Hubbard model at $U/t=1,4$ and $8$ (in ascending order
following the black arrow). The top panel corresponds to the half-band filling
($n=1$) case while the bottom panel corresponds to a three-quarter filling
case ($n=3/2$). Results are shown using the Töws–Pastor functional (full
orange lines) and the Müller functional (dashed blue lines).

In Fig. 4, the occupation numbers of the natural orbitals are represented as a
function of the plane-wave number $k$, at different electron fillings $n$,
different strengths of the Coulomb repulsion $U/t$, and for both the
Töws–Pastor and Müller functionals. At half-band filling ($n=1$), the
Töws–Pastor and the Müller functionals lead to very similar occupation numbers
with a $k$-dependency that is close to a step function, which is due the local
nature of the Coulomb interaction. The step-function-like behaviour also
explains the fast convergence of the algorithm as it corresponds to a quasi
homogeneous renormalization of all off-diagonal terms of the density matrix in
the lattice site representation, and thus to a renormalization of the kinetic
energy. Turning to the three-quarter band filling, this behavior remains,
however with large discrepancies between the Töws–Pastor and the Müller
functional. As already mentioned away from half-band filling, the Müller
functional leads to negative double occupations in particular when $U/t$
increases, which comes along with a drastic reduction of the kinetic energy
and thus drastic modification of the occupation numbers. In contrast,
occupation numbers obtained with the Töws–Pastor functional are less dependent
on $U/t$, which is consistent with the kinetic energy remaining finite even in
the strongly interacting regime away from half-band filling.

Figure 5: SOFT exchange and correlation potential $V_{xc}$ for the 1D Hubbard
model as a function the electron filling $n$ and at different values of the
Coulomb integral $U/t$. Results for the Töws–Pastor functional are compared
with BALDA and with BA. The embedded subplot shows the curve associated to
$U/t=1$.

Let us now investigate the SOFT exchange and correlation potential $V_{xc}$
obtained by Eq. (33). In Fig. 5, $V_{xc}$ obtained using the Töws–Pastor
functional is compared with the exact BA one and with the Bethe–Ansatz Local
Density Approximtion (BALDA) Lima  _et al._ (2002, 2003), with respect to the
electron filling $n$ and for different values of $U/t$. $V_{xc}$ is not shown
for the Müller functional as it leads to unphysical results away from half-
band filling. A major positive result for the Töws–Pastor functional is the
ability to fairly reproduce the $V_{xc}$ discontinuity at half-band filling,
enabling the description of the Mott–Hubbard transition from metal to
insulating regimes. This highlights the ability to design sound exchange and
correlation potentials for SOFT from 1-RDMFT approximations. Away from half-
filling ($0\leq n<1$ and $1<n\leq 2$), the Töws–Pastor functional leads to a
monotonous behavior of $V_{xc}$, in contrast to BALDA that shows unphysical
inflection points and even positive values for $V_{xc}$, particularly in the
weakly correlated regime away from the half-filled case (see the zoomed panel
in Fig. 5). This unphysical behaviour was already observed and shown in Eqs.
(31) and (32) in Ref. Senjean  _et al._ (2018). Finally, in contrast to BALDA
that has been specifically designed for the 1D Hubbard model, the Töws–Pastor
functional can a priori be applied to any lattice dimension.

###  V.2 Molecular sytems

In contrast to the Hubbard model, both Eqs. (17) and (18) have to be solved
for molecular systems. As a proof of concept, we apply our three different
DIVA implementations described in Sec. IV.2 to the hydrogen molecule H2 in the
6-31G basis, using the Müller functional. As DIVA is general and not tight to
any particular orbital representation, we decided to work in the orthogonal
atomic orbital (OAO) basis whose first orbital occupations are given by a DFT
calculation using a local spin density approximation functional (SVWN) from
the PSI4 program Smith  _et al._ (2020). Numerical energy gradients with
respect to ${\bm{\gamma}}^{\neq}$ and ${\bf{n}}$ have been computed with
finite difference.

####  V.2.1 Convergence of the mono-parameter, multi-parameter, and SOFT-DIVA

Figure 6: 1-RDMFT energy with respect to the number of iterations, using SOFT-
DIVA, mono-parameter DIVA, and multi-parameter DIVA, for $R=1.0$ Å. The 6-31G
basis and the Müller functional were used.

In Fig. 6, we compare the convergence of the energy with respect to the number
of iterations, using either the mono- or multi-parameter DIVA implementations,
for an interatomic distance of $R=1.0$ Å. A convergence criteria of $10^{-8}$
Hartree for the energy and $10^{-5}$ for the Frobenius norm of the 1-RDM have
been used. As expected, the multi-parameter implementation converges in much
fewer iterations than the mono-parameter one, although it does not mean that
it is computationally faster than the mono-parameter DIVA. Indeed, each
iteration of mono-parameter DIVA has the same cost, i.e. the cost of the
optimization of a single variational parameter $z$. On the contrary, the
number of parameters to optimize in multi-parameter DIVA increases at each
additional iteration. Hence, there is a trade-off between the number of
iterations and the cost of optimizing many parameters at once. Note that in
this work, multi-parameter DIVA was always computationally much faster than
mono-parameter DIVA. Besides, the former always converges in around 20
iterations while thousands of iterations are sometimes required using mono-
parameter DIVA for some interatomic distances. This very slow convergence has
also been observed in other gradient-based algorithms and is attributed to the
fact that NOFs are not proper functionals of the 1RDM, since these functionals
are not invariant with respect to the rotation of two orbitals having
degenerate occupation numbers Cancès and Pernal (2008). We believe that the
convergence of mono-parameter DIVA could be significantly improved by using
analytical gradients. Alternatively, considering a simultaneous perturbation
stochastic approximation can reduce significantly the number of energy
evaluations needed to compute the numerical gradient, compared to standard
finite difference. This is left for future work.

Let us now consider the convergence of SOFT-DIVA, where the multi-parameter
implementation has been used for the optimization of the off-diagonal elements
of the 1-RDM [Eq. (18)]. As readily seen in Fig. 6, there is no numerical
advantage in splitting the variational principle into two parts. Indeed,
although constraining the orbital occupations fixed leads to faster
convergence when solving Eq. (18), one has to do it for each new orbital
occupations obtained when solving Eq. (17) with the previously obtained Hxc
potential, which is manifested by each peak of SOFT-DIVA in Fig. 6. However,
in practice it allows for the estimation of the Hxc potential, functional of
the orbital occupation, at the Müller functional level of approximation.

####  V.2.2 Energies and occupation numbers

Figure 7: 1-RDMFT energy with respect to the interatomic distance in the 6-31G
basis using the Müller functional. SOFT-DIVAn corresponds to the $n$-th
iteration of SOFT-DIVA, i.e. when the off-diagonal elements have been
optimized for a fixed orbital occupation.

Figure 8: Occupation numbers of each natural orbital $i$ denoted NOi with
respect to the interatomic distance. Dots correspond to the reference results
obtained with MOLGW.

In order to check the accuracy of our implementation, we compare the energy
and occupation numbers obtained by multi-parameter DIVA with the 1-RDMFT
implemented in MOLGW Bruneval  _et al._ (2016); Rodriguez-Mayorga (2022). As
readily seen in Figs. 7 and 8, the multi-parameter DIVA energy and occupation
numbers are on top of the reference ones obtained with MOLGW. One can also see
how the fractional occupation numbers of the natural orbitals in Fig. 8 are
slowly getting degenerate when stretching the bond, which is a very well known
result showing that 1-RDMFT can easily account for static correlation, in
contrast to KS-DFT for which the 1-RDM is idempotent by construction.

Turning to SOFT-DIVA, by stopping the algorithm at each new update of the Hxc
potential, one can see how much the change in orbital occupations impacts the
energy. In other words, one can look at the best 1-RDM functional energy
corresponding to a fixed orbital occupation (in the OAO basis), where only the
off-diagonal elements of the 1-RDM have been optimized. In Fig. 7, the energy
of the first iteration deviates from the reference one when stretching the
hydrogen bond, meaning that the orbital occupations are not optimal in this
case. This deviation is already corrected at the second iteration,
corresponding to the first update of the Hxc potential that gives sufficiently
accurate orbital occupations for the energy to be on top of the converged one.

Turning to the orbital occupations, the occupations of the OAO are given for
the first three iterations of SOFT-DIVA in Fig. 9. Note that the initial
orbital occupations of SOFT-DIVA are obtained by taking the diagonal (in the
OAO basis) of the Hxc potential from the converged DFT calculation, which
explains why they are not equal but close to the DFT occupations. According to
Fig. 9, the OAO occupations at the second iteration of SOFT-DIVA are indeed
much closer to the converged one, thus explaining why the energy in Fig. 7
matches the converged one. However, another iteration is required to get very
accurate OAO occupations.

Figure 9: OAO occupations with respect to the interatomic distance. For
symmetry reasons, there are only two different OAO occupation values
(corresponding to the 1s and 2s of the hydrogens), shown in blue and orange,
respectively. Converged DIVA occupations and DFT ones are depicted with
circles and squares, respectively. SOFT-DIVAn corresponds to the $n$-th
iteration of SOFT-DIVA.

##  VI Conclusions and perspectives

We have introduced a new variational minimization scheme for 1-RDMFT called
density-matrix interpolation variational ansatz (DIVA) that works in any
orbital basis representation while satisfying the ensemble
$N$-representability conditions. We have derived different implementations and
successfully tested them on the uniform Hubbard model with the Müller and the
Töws–Pastor functionals, as well as the hydrogen molecule with the Müller
functional. In summary, every implementations give the correct 1-RDM and
energy. The mono-parameter DIVA converges in a few iterations for the Hubbard
model, but exhibits very slow convergence for the hydrogen molecule. The
multi-parameter DIVA led to a much faster convergence, at the expense of more
involved optimizations of the variational parameters. Finally, by allowing the
estimation of the exchange and correlation potential for a given orbital
occupation vector, the SOFT-DIVA implementation contributes to the development
of orbital occupation functionals for chemistry. Besides, we stress that DIVA
is more adapted to the computational basis representation, as it would require
less swapping between the computational and the NO basis for the calculation
of energy and energy gradients, thus motivating the future development of
1-RDM functionals beyond NOFs. Work is in progress in these directions.

##  VII Acknowledgments

We thank Mauricio Rodriguez-Mayorga for fruitful discussions about 1-RDMFT in
general. QM and MS would like to thank the ANR (Grant No. ANR-19-CE29-0002
DESCARTES project) for funding.

##  Appendix A Convexity of the ensemble $N$-representable 1-RDM domain.

The ensemble $N$-representable domaine $\Omega_{{\bm{\gamma}}}$ is convex if,
for any ensemble $N$-representable 1-RDM ${\bm{\gamma}}^{A}$ and
${\bm{\gamma}}^{B}$, ${\bm{\gamma}}=z{\bm{\gamma}}^{A}+(1-z){\bm{\gamma}}^{B}$
with $z\in[0,1]$ is also ensemble $N$-representable.
Let us consider a $N$ dimensional vector ${\bf{x}}$. Thus we have

${\bf{x}}^{T}{\bm{\gamma}}{\bf{x}}=z{\bf{x}}^{T}{\bm{\gamma}}^{A}{\bf{x}}+(1-z){\bf{x}}^{T}{\bm{\gamma}}^{B}{\bf{x}}\geq
0$ (35)

since ${\bf{x}}^{T}{\bm{\gamma}}^{A}{\bf{x}}\geq 0$ and
${\bf{x}}^{T}{\bm{\gamma}}^{B}{\bf{x}}\geq 0$ (${\bm{\gamma}}^{A}$ and
${\bm{\gamma}}^{B}$ are positive semi-definite) and $0\leq z\leq 1$.
Consequently, ${\bm{\gamma}}$ is also positive semi-definite. Similarly,

${\bf{x}}^{T}\left(\mathbb{1}-{\bm{\gamma}}\right){\bf{x}}=z{\bf{x}}^{T}\left(\mathbb{1}-{\bm{\gamma}}^{A}\right){\bf{x}}+(1-z){\bf{x}}^{T}\left(\mathbb{1}-{\bm{\gamma}}^{B}\right){\bf{x}}\geq
0$ (36)

since $0\leq z\leq 1$ and
${\bf{x}}^{T}\left(\mathbb{1}-{\bm{\gamma}}^{A}\right){\bf{x}}\geq 0$ and
${\bf{x}}^{T}\left(\mathbb{1}-{\bm{\gamma}}^{B}\right){\bf{x}}\geq 0$) since
$\mathbb{1}-{\bm{\gamma}}^{A}$ and $\mathbb{1}-{\bm{\gamma}}^{B}$ are positive
semi-definite. Consequently, $\mathbb{1}-{\bm{\gamma}}$ is also semi positive
definite. Finally, ${\bm{\gamma}}$ and $\mathbb{1}-{\bm{\gamma}}$ being
positive semi-definite, ${\bm{\gamma}}$ is ensemble $N$-representable and
$\Omega_{{\bm{\gamma}}}$ is convex.

##  Appendix B Interior and boundary of $\Omega_{\gamma}$

We define the neighborhood $\mathcal{N}({\bm{\gamma}})$ of a matrix
${\bm{\gamma}}\in\mathcal{M}_{N}(\mathbb{R})$ as the open ball

$\displaystyle\mathcal{N}({\bm{\gamma}})=\left\\{{\bm{\gamma}}^{\prime}\in\mathcal{M}_{N}(\mathbb{R})\;|\;||{\bm{\gamma}}-{\bm{\gamma}}^{\prime}||<\epsilon,\epsilon\rightarrow
0\right\\}$ (37)

where $||\cdot||$ refers to the Frobenius norm. We consider
${\bm{\gamma}}^{\prime}\in\mathcal{N}({\bm{\gamma}})$ and
$\delta{\bm{\gamma}}=\left({\bm{\gamma}}^{\prime}-{\bm{\gamma}}\right)/\lambda$
with $\lambda=||{\bm{\gamma}}^{\prime}-{\bm{\gamma}}||<\epsilon$, such that

$\displaystyle{\bm{\gamma}}^{\prime}={\bm{\gamma}}+\lambda\delta{\bm{\gamma}}.$
(38)

Following a perturbation approach with
${\bm{\gamma}}^{\prime}={\bm{\gamma}}+\lambda\delta{\bm{\gamma}}$ and
$\lambda\rightarrow 0$, we propose to develop the eigenvalues
${\bm{\eta}}^{\prime}={\bm{\eta}}+\lambda\delta{\bm{\eta}}$ and eigenvectors
${\bf{u}}^{\prime}={\bf{u}}+\lambda\delta{\bf{u}}$ as a power of infinitesimal
$\lambda$ where the ${\bm{\eta}}$ and ${\bf{u}}$ correspond to the eigenvalues
diagonal matrix and eigenvector matrix of ${\bm{\gamma}}$, respectively,

${\bm{\gamma}}{\bf{u}}^{\dagger}={\bf{u}}^{\dagger}{\bm{\eta}}.$ (39)

The unitary transformation ${\bf{u}}^{\prime}$ that diagonalizes
${\bm{\gamma}}^{\prime}$ corresponds to

${\bm{\gamma}}^{\prime}{\bf{u}}^{\prime\dagger}={\bf{u}}^{\prime\dagger}{\bm{\eta}}^{\prime},$
(40)

and Eq. (39) is recovered at the zero order development. At first order it
reads

$\delta{\bm{\gamma}}{\bf{u}}^{\dagger}+{\bm{\gamma}}\delta{\bf{u}}^{\dagger}-\delta{\bf{u}}^{\dagger}{\bm{\eta}}-{\bf{u}}^{\dagger}\delta{\bm{\eta}}=0.$
(41)

Multiplying by ${\bf{u}}$ on the left side leads to

$\displaystyle{\bf{u}}\delta{\bm{\gamma}}{\bf{u}}^{\dagger}+{\bf{u}}{\bm{\gamma}}{\delta{\bf{u}}}^{\dagger}-{\bf{u}}{\delta{\bf{u}}}^{\dagger}{\bm{\eta}}-{\bf{u}}{\bf{u}}^{\dagger}\delta{\bm{\eta}}=0$
$\displaystyle{\bf{u}}\delta{\bm{\gamma}}{{\bf{u}}}^{\dagger}-\delta{\bm{\eta}}+\left[{\bm{\eta}},{\bf{u}}{\delta{\bf{u}}}^{\dagger}\right]=0.$
(42)

On the diagonal part of Eq. (42), the last term on the right hand side
cancels, such that

$\delta{\bm{\eta}}={\bf{u}}\delta{\bm{\gamma}}{{\bf{u}}}^{\dagger}.$ (43)

Let us define

$\displaystyle d({\bm{\gamma}})={\rm min}\left(\eta_{1},1-\eta_{N}\right),$
(44)

where the eigenvalues are sorted in ascending order. $d({\bm{\gamma}})$
appears as a pseudo-distance from $\delta\Omega_{\bm{\gamma}}$ as
$d({\bm{\gamma}})\geq 0$ for ${\bm{\gamma}}\in\Omega_{\bm{\gamma}}$,
$d({\bm{\gamma}})>0$ for ${\bm{\gamma}}\in\mathring{\Omega}_{{\bm{\gamma}}}$,
$d({\bm{\gamma}})=0$ for ${\bm{\gamma}}\in\delta\Omega_{\bm{\gamma}}$ and
finally, $d({\bm{\gamma}})<0$ for non-representable matrices
${\bm{\gamma}}\notin\Omega_{\bm{\gamma}}$. It follows that

$\displaystyle d({\bm{\gamma}}^{\prime})={\rm
min}\left(\eta_{1}+\delta\eta_{1},1-(\eta_{N}+\delta\eta_{N})\right),$ (45)

and $d({\bm{\gamma}}^{\prime})$ is contained in the interval
$d({\bm{\gamma}}^{\prime})$ $\in$
$]d({\bm{\gamma}})-\lambda,d({\bm{\gamma}})+\lambda[$ since
$|\delta{\bm{\eta}}|\leq\lambda|{\bf{u}}\delta{\bm{\gamma}}{\bf{u}}^{\dagger}|\leq\lambda$.

Consequently, for each ${\bm{\gamma}}\in\mathring{\Omega}_{{\bm{\gamma}}}$
there is a neighborhood with $\epsilon<d({\bm{\gamma}})$ such that
$\mathcal{N}({\bm{\gamma}})\subseteq\Omega_{{\bm{\gamma}}}$, thus making
$\mathring{\Omega}_{{\bm{\gamma}}}$ the topological interior of
$\Omega_{{\bm{\gamma}}}$. On the contrary, for
${\bm{\gamma}}\in\delta\Omega_{{\bm{\gamma}}}$, since $d({\bm{\gamma}})=0$ it
becomes trivial to find matrices in $\mathcal{N}({\bm{\gamma}})$ that do not
belong to $\Omega_{{\bm{\gamma}}}$. Hence, $\delta\Omega_{{\bm{\gamma}}}$ is
the topological boundary of $\Omega_{{\bm{\gamma}}}$.

##  Appendix C Density-matrix decomposition

We show that for any non-idempotent ${\bm{\gamma}}$, there exists at least one
decomposition in terms of (up to $2^{N}$) idempotent matrices such that
${\bm{\gamma}}=\sum_{i}z_{i}{\bm{\gamma}}_{*}^{(i)}$ with $\sum_{i}z_{i}=1$.
To that aim, we propose the following iterative algorithm. Consider
${\bm{\gamma}}\in\Omega_{{\bm{\gamma}}}\backslash\delta\Omega_{{\bm{\gamma}}}^{*}$
that has $F({\bm{\gamma}})$ fractional eigenvalues.

  1. 1.

Choose one among the $F({\bm{\gamma}})$ fractional eigenvalues $\eta_{k}$

  2. 2.

Construct, in the diagonal representation of ${\bm{\gamma}}$,

$\displaystyle{\bm{\gamma}}^{\prime}(\lambda)={\bm{\gamma}}+\lambda{\bm{\gamma}}^{(k)}$
(46) $\displaystyle\gamma^{(k)}_{ij}=0\,\forall\,(i,j)\neq(k,k),\,{\rm
and}\,\gamma^{(k)}_{kk}=1$ (47)

with $\lambda$ a scalar and
${\bm{\gamma}}^{\prime}(\lambda)\in\Omega_{{\bm{\gamma}}}$ for
$\lambda\in[-\lambda^{A},\lambda^{B}]$ with $\lambda^{A}=\eta_{k}$ and
$\lambda^{B}=1-\eta_{k}$. By construction,
$F\left({\bm{\gamma}}^{\prime}(-\lambda^{A})\right)=F\left({\bm{\gamma}}^{\prime}(\lambda^{B})\right)=F({\bm{\gamma}})-1$.

  3. 3.

Decompose ${\bm{\gamma}}$ in terms of matrices with higher number of integer
eigenvalues,

$\displaystyle{\bm{\gamma}}=(1-\eta_{k}){\bm{\gamma}}^{\prime}(-\lambda^{A})+\eta_{k}{\bm{\gamma}}^{\prime}(\lambda^{B})$
(48)

Repeat steps 1 and 2 for ${\bm{\gamma}}^{\prime}(-\lambda^{A})$ and
${\bm{\gamma}}^{\prime}(\lambda^{B})$ and insert their decomposition into Eq.
(48). Continue the algorithm until

$\displaystyle{\bm{\gamma}}=\sum_{i}z_{i}{\bm{\gamma}}_{*}^{(i)},$ (49)

with ${\bm{\gamma}}_{*}^{(i)}\in\delta\Omega_{\bm{\gamma}}^{*}$ and
$\sum_{i}z_{i}=1$. Each iteration of the algorithm multiplies by two the
number of matrices in the decomposition, making the sum in Eq. (49) running
ultimately up to $2^{F({\bm{\gamma}})}$ elements.

## References

  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. Sham, Phys. Rev. 140, A1133 (1965).
  * Marques  _et al._ (2012) M. A. L. Marques, M. J. T. Oliveira, and T. Burnus, Comput. Phys. Commun. 183, 2272 (2012).
  * Mardirossian and Head-Gordon (2017) N. Mardirossian and M. Head-Gordon, Mol. Phys. 115, 2315 (2017).
  * Cramer and Truhlar (2009) C. J. Cramer and D. G. Truhlar, Phys. Chem. Chem. Phys. 11, 10757 (2009).
  * Cohen  _et al._ (2008) A. J. Cohen, P. Mori-Sánchez, and W. Yang, Science 321, 792 (2008).
  * Löwdin (1955) P.-O. Löwdin, Phys. Rev. 97, 1474 (1955).
  * Gilbert (1975) T. L. Gilbert, Phys. Rev. B 12, 2111 (1975).
  * Coleman (1963) A. J. Coleman, Rev. Mod. Phys. 35, 668 (1963).
  * Haar (1961) D. T. Haar, Rep. Prog. Phys. 24, 304 (1961).
  * Piris and Pernal (2017) M. Piris and K. Pernal, Phys. Rev. A 96, 046501 (2017).
  * Altunbulak and Klyachko (2008) M. Altunbulak and A. Klyachko, Commun. Math. Phys. 282, 287 (2008).
  * Schilling  _et al._ (2018) C. Schilling, M. Altunbulak, S. Knecht, A. Lopes, J. D. Whitfield, M. Christandl, D. Gross, and M. Reiher, Phys. Rev. A 97, 052503 (2018).
  * Castillo  _et al._ (2023) F. Castillo, J.-P. Labbé, J. Liebert, A. Padrol, E. Philippe, and C. Schilling, in  _Annales Henri Poincaré_, Vol. 24 (Springer, 2023) pp. 2241–2321.
  * Theophilou  _et al._ (2015a) I. Theophilou, N. N. Lathiotakis, M. A. L. Marques, and N. Helbig, J. Chem. Phys. 142 (2015a).
  * Coleman (1972) A. J. Coleman, J. Math. Phys. 13, 214 (1972).
  * Valone (1980) S. M. Valone, J. Chem. Phys. 73, 1344 (1980).
  * Müller (1984) A. M. K. Müller, Phys. Lett. A 105, 446 (1984).
  * Csányi and Arias (2000) G. Csányi and T. A. Arias, Phys. Rev. B 61, 7348 (2000).
  * Csányi  _et al._ (2002) G. Csányi, S. Goedecker, and T. A. Arias, Phys. Rev. A 65, 032510 (2002).
  * Buijse and Baerends (2002) M. A. Buijse and E. J. Baerends, Mol. Phys. 100, 401 (2002).
  * Sharma  _et al._ (2008) S. Sharma, J. K. Dewhurst, N. N. Lathiotakis, and E. K. U. Gross, Phys. Rev. B 78, 201103 (2008).
  * Lathiotakis  _et al._ (2009) N. N. Lathiotakis, N. Helbig, A. Zacarias, and E. K. U. Gross, J. Chem. Phys. 130, 064109 (2009).
  * Goedecker and Umrigar (1998) S. Goedecker and C. J. Umrigar, Phys. Rev. Lett. 81, 866 (1998).
  * Holas (1999) A. Holas, Phys. Rev. A 59, 3454 (1999).
  * Cioslowski and Pernal (1999) J. Cioslowski and K. Pernal, J. Chem. Phys. 111, 3396 (1999).
  * Cioslowski  _et al._ (2023) J. Cioslowski, C. Schilling, and R. Schilling, J. Chem. Phys. 158 (2023).
  * Piris and Otto (2003) M. Piris and P. Otto, Int. J. Quantum Chem. 94, 317 (2003).
  * Piris (2004) M. Piris, Recent Res. Devel. Quantum Chem. 4, 1 (2004).
  * Leiva and Piris (2005) P. Leiva and M. Piris, J. Mol. Struct. THEOCHEM 719, 63 (2005).
  * Piris (2006) M. Piris, Int. J. Quantum Chem. 106, 1093 (2006).
  * Piris  _et al._ (2007) M. Piris, X. Lopez, and J. M. Ugalde, J. Chem. Phys. 126, 214103 (2007).
  * Piris and Ugalde (2009) M. Piris and J. M. Ugalde, J. Comput. Chem. 30, 2078 (2009).
  * Piris  _et al._ (2010) M. Piris, J. M. Matxain, X. Lopez, and J. M. Ugalde, J. Chem. Phys. 132 (2010).
  * Piris  _et al._ (2011) M. Piris, X. Lopez, F. Ruipérez, J. M. Matxain, and J. M. Ugalde, J. Chem. Phys. 134, 164102 (2011).
  * Piris (2013) M. Piris, Int. J. Quantum Chem. 113, 620 (2013).
  * Piris and Ugalde (2014) M. Piris and J. M. Ugalde, Int. J. Quantum Chem. 114, 1169 (2014).
  * Lopez  _et al._ (2015) X. Lopez, M. Piris, F. Ruipeérez, and J. M. Ugalde, J. Phys. Chem. A 119, 6981 (2015).
  * Kutzelnigg and Mukherjee (1999) W. Kutzelnigg and D. Mukherjee, J. Chem. Phys. 110, 2800 (1999).
  * Ramos-Cordoba  _et al._ (2015) E. Ramos-Cordoba, X. Lopez, M. Piris, and E. Matito, J. Chem. Phys. 143, 164112 (2015).
  * Cioslowski  _et al._ (2015) J. Cioslowski, M. Piris, and E. Matito, J. Chem. Phys. 143 (2015).
  * Mitxelena  _et al._ (2017) I. Mitxelena, M. Piris, and M. Rodriguez-Mayorga, J. Phys. Condens. Matter 29, 425602 (2017).
  * Rodriguez-Mayorga  _et al._ (2017) M. Rodriguez-Mayorga, E. Ramos-Cordoba, M. Via-Nadal, M. Piris, and E. Matito, Phys. Chem. Chem. Phys. 19, 24029 (2017).
  * Mitxelena  _et al._ (2018) I. Mitxelena, M. Rodriguez-Mayorga, and M. Piris, Eur. Phys. J. B 91, 109 (2018).
  * Piris (2021) M. Piris, Phys. Rev. Lett. 127, 233001 (2021).
  * Pernal (2005) K. Pernal, Phys. Rev. Lett. 94, 233002 (2005).
  * Lathiotakis  _et al._ (2014) N. N. Lathiotakis, N. Helbig, A. Rubio, and N. I. Gidopoulos, Phys. Rev. A 90, 032511 (2014).
  * Theophilou  _et al._ (2015b) I. Theophilou, N. N. Lathiotakis, N. I. Gidopoulos, A. Rubio, and N. Helbig, J. Chem. Phys. 143, 054106 (2015b).
  * Goedecker and Colombo (1994) S. Goedecker and L. Colombo, Phys. Rev. Lett. 73, 122 (1994).
  * Baldsiefen and Gross (2013) T. Baldsiefen and E. K. U. Gross, Comput. Theor. Chem. 1003, 114 (2013).
  * Saubanère  _et al._ (2016) M. Saubanère, M. B. Lepetit, and G. M. Pastor, Phys. Rev. B 94, 045102 (2016).
  * Wang and Baerends (2022) J. Wang and E. J. Baerends, Phys. Rev. Lett. 128, 013001 (2022).
  * Cancès and Pernal (2008) E. Cancès and K. Pernal, J. Chem. Phys. 128, 134108 (2008).
  * Gunnarsson and Schönhammer (1986) O. Gunnarsson and K. Schönhammer, Phys. Rev. Lett. 56, 1968 (1986).
  * Schönhammer  _et al._ (1995) K. Schönhammer, O. Gunnarsson, and R. Noack, Phys. Rev. B 52, 2504 (1995).
  * Lima  _et al._ (2002) N. Lima, L. Oliveira, and K. Capelle, Europhys. Lett. 60, 601 (2002).
  * Lima  _et al._ (2003) N. A. Lima, M. F. Silva, L. N. Oliveira, and K. Capelle, Phys. Rev. Lett. 90, 146402 (2003).
  * Fromager (2015) E. Fromager, Mol. Phys. 113, 419 (2015).
  * Coe (2019) J. P. Coe, Phys. Rev. B 99, 165118 (2019).
  * Cartier and Giesbertz (2024) N. G. Cartier and K. J. H. Giesbertz, J. Chem. Theory Comput. (2024).
  * Töws and Pastor (2011) W. Töws and G. M. Pastor, Phys. Rev. B 83, 235101 (2011).
  * Töws and Pastor (2012) W. Töws and G. M. Pastor, Phys. Rev. B 86, 245123 (2012).
  * Töws  _et al._ (2014) W. Töws, M. Saubanère, and G. M. Pastor, Theor. Chem. Acc. 133, 1 (2014).
  * Kamil  _et al._ (2016) E. Kamil, R. Schade, T. Pruschke, and P. E. Blöchl, Phys. Rev. B 93, 085141 (2016).
  * Lieb and Wu (1968) E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968).
  * Shiba (1972) H. Shiba, Phys. Rev. B 6, 930 (1972).
  * Senjean  _et al._ (2018) B. Senjean, N. Nakatani, M. Tsuchiizu, and E. Fromager, Phys. Rev. B 97, 235105 (2018).
  * Smith  _et al._ (2020) D. G. A. Smith, L. A. Burns, A. C. Simmonett, R. M. Parrish, M. C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. Di Remigio, A. Alenaizan, A. M. James, S. Lehtola, J. P. Misiewicz, M. Scheurer, R. A. Shaw, J. B. Schriber, Y. Xie, Z. L. Glick, D. A. Sirianni, J. S. O’Brien, J. M. Waldrop, A. Kumar, E. G. Hohenstein, B. P. Pritchard, B. R. Brooks, H. F. Schaefer, A. Y. Sokolov, K. Patkowski, A. E. DePrince, U. Bozkaya, R. A. King, F. A. Evangelista, J. M. Turney, T. D. Crawford, and C. D. Sherrill, J. Chem. Phys. 152, 184108 (2020).
  * Bruneval  _et al._ (2016) F. Bruneval, T. Rangel, S. M. Hamed, M. Shao, C. Yang, and J. B. Neaton, Comput. Phys. Commun. 208, 149 (2016).
  * Rodriguez-Mayorga (2022) M. Rodriguez-Mayorga, Standalone NOFT module (1.0). Zenodo (2022), 10.5281/zenodo.7352733.