11affiliationtext: CERMICS, Ecole des Ponts and Inria Paris, 6 & 8 avenue
Blaise Pascal, 77455 Marne-la-Vallée, France22affiliationtext: Dipartimento di
Chimica e Chimica Industriale, Università di Pisa, Via G. Moruzzi 13, 56124
Pisa, Italy

# Geometric Optimization of Restricted-Open and Complete Active Space Self-
Consistent Field Wavefunctions

Laurent Vidal  Tommaso Nottoli  Filippo Lipparini  Eric Cancès

###### Abstract

We explore Riemannian optimization methods for Restricted-Open-shell Hartree-
Fock (ROHF) and Complete Active Space Self-Consistent Field (CASSCF) methods.
After showing that ROHF and CASSCF can be reformulated as optimization
problems on so-called flag manifolds, we review Riemannian optimization basics
and their application to these specific problems. We compare these methods to
traditional ones and find robust convergence properties without fine-tuning of
numerical parameters. Our study suggests Riemannian optimization as a valuable
addition to orbital optimization for ROHF and CASSCF, warranting further
investigation.

##  1 Introduction

Orbital optimization is one of the most common task performed in quantum
chemistry calculations. It is the numerical problem associated with Hartree-
Fock (HF) [1] and Kohn-Sham Density Functional Theory (KS DFT) [2], as well as
a component of Complete Active Space Self-Consistent Field (CASSCF)
calculations [3, 4, 5], and is further encountered in orbital optimized post-
Hartree Fock methods [6, 7]. The various algorithms that have been proposed to
tackle this problem can be grouped into two families: fixed point methods,
such as Roothaan’s SCF algorithm [8, 9], and direct optimization methods, such
as quadratically convergent optimization strategies [10, 11]. For HF and DFT,
the former family is the most commonly employed, due to the existence of very
robust implementation that exploit convergence acceleration techniques such as
Pulay’s Direct Inversion in the Iterative Subspace [12, 13, 14] (DIIS),
constraint relaxation methods such as the Optimal Damping Algorithm [15, 16,
17] (ODA), or more sophisticated related techniques such as E-DIIS [18] or
A-DIIS [19]. Nevertheless, direct optimization techniques have received quite
some attention due to their robustness and due to the possibility of
implementing them avoiding dense linear algebra operations (e.g.,
diagonalization of the Fock matrix).

Direct minimization techniques for Restricted Open-Shell Hartree-Fock (ROHF)
calculations are relatively scarce compared to Self-Consistent Field (SCF)
methods, which predominantly feature in quantum chemistry software. Among
direct minimization approaches, noteworthy methods include Geometric Direct
Minimization (GDM) techniques [20] and the QC-SCF algorithm [21].
Additionally, the Second-Order SCF (SOSCF) algorithm [22, 23] and the DIIS-GDM
method [20, 24], which amalgamate aspects of both SCF and direct minimization
strategies, merit mention. Moreover, the CUHF method, as introduced by
Tsuchimochi and Scuseria [25], can be adapted for ROHF computations through
the utilization of a direct minimization procedure designed for Unrestricted
Hartree-Fock (UHF) calculations.

A difficulty associated with the formulation and implementation of direct
minimization techniques is due to the fact that the quantity that needs to be
optimized, such as the molecular orbitals (MO) coefficients, or the density
matrix, needs to satisfy nonlinear constraints. In other words, the
minimization set is not a vector space, but rather a differentiable manifold.

It is well-known that after discretization in a finite basis set, HF and KS
models can be formulated as optimization problems on Stiefel (molecular
orbital formalism) or Grassmann (density matrix formalism) manifolds [26].
These formulations lead to enlightening geometric interpretations of the
Hartree-Fock and Kohn-Sham equations, and to the design and convergence
analysis [27, 28, 29] of robust and efficient direct minimization algorithms.
The purpose of this article is to show that Restricted Open-Shell Hartree-Fock
(ROHF) and Complete Active Space Self-Consistent Field (CASSCF) methods can be
reformulated as optimization problems on so-called flag manifolds [30]. This
allows one to shed new light on the ROHF and CASSCF equations, and the direct
minimization algorithms used to solve these problems. While the work presented
in this article does not lead to game-changing improvements in orbital
optimization, we hope that it will provide the community with a set of
rigorous tools that can be used for further developments, as the ones recently
proposed by some of us for the extrapolation of the SCF density matrix in the
context of ab-initio molecular dynamics simulations [31, 32, 33].

This article is organized as follows. In Section 2, we briefly recall the
high-spin ROHF and CASSCF orbital optimization problem in terms of both
density-matrix (DM) and molecular orbitals (MO) formulations and provide a
simple geometric interpretation of the ROHF and CASSCF equations. In Section
3, we review the basic concepts of geometric optimization (Riemannian gradient
and Hessian, vector transport, affine connection, geodesic, retraction). In
Section 4, we discuss more specifically geometric optimization for ROHF and
CASSCF, providing also the tools to translate any algorithm formulated in the
MO formalism into the DM formalism and viceversa. We then provide geometric
interpretations of existing direct minimization algorithms for ROHF and
CASSCF, and propose new ones. We also introduce a direct optimization method
circumventing the use of virtual orbitals, which is useful for ROHF in large
basis sets (i.e. planewaves, finite elements, or wavelets). Numerical results
are reported in Section 5.

##  2 ROHF and CASSCF

In this section, we briefly recapitulate the orbital optimization problem for
ROHF and CASSCF and introduce the manifolds associated with the MO and DM
formalisms. ROHF and CASSCF methods indeed share common features. They both
involve

  * •

a set of $N_{I}$ doubly-occupied molecular orbitals, often called internal
orbitals;

  * •

a set of $N_{A}$ partially-occupied molecular orbitals, often called active
orbitals,

the latter being orthogonal to the former. Consider a molecular system with
$N$ electrons discretized in a basis set of size ${\mathcal{N}}_{\rm b}$. We
denote the electronic Hamiltonian by

$\widehat{H}_{N}=-\frac{1}{2}\sum_{i=1}^{N}\nabla_{{\mathbf{r}}_{i}}^{2}+\sum_{i=1}^{N}V_{\rm
nuc}({\mathbf{r}}_{i})+\sum_{1\leq i<j\leq
N}\frac{1}{|{\mathbf{r}}_{i}-{\mathbf{r}}_{j}|}.$ (1)

#### Notation

To avoid possible misunderstandings, we clarify here the notation that will be
adopted throughout the paper. Let us assume that we have discretized the
problem using an orthonormal set of atomic orbitals (AO)
$\\{\chi_{\mu}\\}_{\mu=1}^{{\mathcal{N}}_{\rm b}}$, obtained, for instance, by
Löwdin orthogonalization of a usual Gaussian-type or Slater-type basis. This
implies, for the overlap matrix:

$S_{\mu\nu}=\langle\chi_{\mu}|\chi_{\nu}\rangle=\delta_{\mu\nu}.$

We use greek letters $\mu,\nu,\ldots$ to label AOs. Molecular orbitals
$\\{\phi_{p}\\}_{p=1}^{{\mathcal{N}}_{\rm b}}$ are written as linear
combinations of atomic orbitals

$\phi_{p}=\sum_{\mu=1}^{{\mathcal{N}}_{\rm b}}C_{\mu p}\chi_{\mu}$

where the coefficient matrix $C\in{\mathbb{R}}^{{\mathcal{N}}_{\rm
b}\times{\mathcal{N}}_{\rm b}}$ is, due to our choice of orthogonal AOs, an
orthogonal matrix, i.e., $CC^{T}=C^{T}C=I_{{\mathcal{N}}_{\rm b}}$. We divide
the molecular orbitals into three sets: $N_{I}$ _internal_ orbitals, that are
always doubly occupied, $N_{A}$ _active_ orbitals, that are singly occupied in
ROHF and have varying occupation for CASSCF, and $N_{E}$ _external_ orbitals,
that are always empty. We use $i,j,\ldots$ to label internal orbitals,
$u,v,\ldots$ to label active orbitals, $a,b,\ldots$ to label external
orbitals, and $p,q,\ldots$ for generic ones. We call $\Pi^{I}$ and $\Pi^{A}$
the orthogonal projectors on the space spanned by internal and active
orbitals, respectively, in the AO basis, i.e.,

$\Pi^{I}_{\mu\nu}=\sum_{i=1}^{N_{I}}C_{\mu i}C_{\nu
i},\quad\Pi^{A}_{\mu\nu}=\sum_{u=N_{I}+1}^{N_{I}+N_{A}}C_{\mu u}C_{\nu u}$ (2)

$P\in{\mathbb{R}}^{{\mathcal{N}}_{\rm b}\times{\mathcal{N}}_{\rm b}}$ denotes
the one-body reduced density matrix (1-RDM) in the AO basis, while
$\gamma\in{\mathbb{R}}^{{\mathcal{N}}_{\rm b}\times{\mathcal{N}}_{\rm b}}$
denotes the one-body reduced density matrix in the MO basis. We call
$m_{\alpha}$ and $m_{\beta}$ the number of $\alpha$ and $\beta$ _active_
electrons, respectively, so that the total number of electrons is given by
$N=2N_{I}+m_{\alpha}+m_{\beta}$. In high-spin ROHF, $m_{\alpha}=N_{A}$ and
$m_{\beta}=0$. Using these conventions, the high-spin ROHF density matrix is
given by

$P^{\rm ROHF}=2\Pi^{I}+\Pi^{A}$ (3)

where we note that the $\alpha$ and $\beta$ one-body spin density matrices
(1-SDM) are given by

$P^{\rm ROHF,\alpha}=\Pi^{I}+\Pi^{A},\quad P^{\rm ROHF,\beta}=\Pi^{I}$ (4)

For CASSCF, the MO 1-RDM $\gamma^{\Psi^{\rm CAS}}$ has a block structure, with

$\gamma^{\Psi^{\rm CAS}}_{ij}=2\delta_{ij},\quad\gamma^{\Psi^{\rm
CAS}}_{uv}=\langle\Psi^{\rm CAS}|\hat{E}_{uv}|\Psi^{\rm
CAS}\rangle,\quad\gamma^{\Psi^{\rm CAS}}_{ab}=0,$ (5)

where

$\hat{E}_{pq}=\hat{a}^{\dagger}_{p\alpha}\hat{a}_{q\alpha}+\hat{a}^{\dagger}_{p\beta}\hat{a}_{q\beta}$

is the spin-traced singlet excitation operator, and all other blocks
vanishing. In the AO basis, $P^{\Psi^{\rm CAS}}=C\gamma^{\Psi^{\rm CAS}}C^{T}$
satisfies thus the following relation:

$2\Pi^{I}\leq P^{\Psi^{\rm CAS}}\leq 2\Pi^{I}+\Pi^{A}.$ (6)

To define 1-SDM for CASSCF, we need to introduce the active spin densities
$\gamma_{uv}^{\Psi^{\rm CAS},\sigma}$, for $\sigma=\alpha,\beta$ which are
defined as follows:

$\gamma^{\Psi^{\rm CAS},\sigma}_{ij}=\delta_{ij},\quad\gamma^{\Psi^{\rm
CAS},\sigma}_{uv}=\langle\Psi^{\rm
CAS}|\hat{a}^{\dagger}_{u\sigma}\hat{a}_{v\sigma}|\Psi^{\rm CAS}\rangle,$ (7)

with all the other blocks (i.e., internal-active and all blocks with at least
one external index) vanishing. The AO spin densities are then given by

$P^{\Psi^{\rm CAS},\sigma}=C\gamma^{{\rm CAS},\sigma}C^{T}$ (8)

and it holds that

$\Pi^{I}\leq P^{\Psi^{\rm CAS},\sigma}\leq\Pi^{I}+\Pi^{A}.$ (9)

Direct optimization methods for ROHF and CASSCF can be divided into two
groups, depending on the degrees of freedom used for performing the
optimization. In the MO formalism, the main variable is the coefficient matrix
$C$. As mentioned previously, it is an orthonormal matrix, which can be seen
as a point of the orthogonal group $O({\mathcal{N}}_{b})$. In the DM
formulation, the main variable is the pair of orthogonal projectors $(\Pi^{\rm
I},\Pi^{\rm A})$, which can be identified with a point of the set

$\displaystyle{\mathcal{M}}_{\rm DM}:=\biggl{\\{}(\Pi^{\rm I},\Pi^{\rm A})$
$\displaystyle\in{\mathbb{R}}^{{\mathcal{N}}_{\rm b}\times{\mathcal{N}}_{\rm
b}}_{\rm sym}\times{\mathbb{R}}^{{\mathcal{N}}_{\rm b}\times{\mathcal{N}}_{\rm
b}}_{\rm sym}\mbox{ s.t. }(\Pi^{\rm I})^{2}=\Pi^{\rm I},\;(\Pi^{\rm
A})^{2}=\Pi^{\rm A}$ $\displaystyle\;{\rm Tr}(\Pi^{\rm I})=N_{I},\;{\rm
Tr}(\Pi^{\rm A})=N_{A},\mbox{ and }\Pi^{\rm I}\Pi^{\rm A}=0\biggr{\\}}.$ (10)

We will see later that the above set has a nice geometrical structure: it can
be canonically identified with the flag manifold ${\mathcal{M}}_{\rm
Flag}:={\rm Flag}(N_{\rm I},N_{\rm I}+N_{\rm
A};{\mathbb{R}}^{{\mathcal{N}}_{\rm b}})$.

The passage between MO and DM parameterization is done by the map

$\zeta:O({\mathcal{N}}_{b})\ni C\mapsto\left(\Pi^{\rm I},\Pi^{\rm
A}\right)\in{\mathcal{M}}_{\rm MO}\qquad\mbox{with }\Pi^{\rm I},\Pi^{\rm
A}\mbox{ given by Eq. \ref{eq:projectors}}.$ (11)

The dimension of the MO manifold $O({\mathcal{N}}_{b})$ is
$\frac{{\mathcal{N}}_{b}({\mathcal{N}}_{b}-1)}{2}$ (i.e. the number of degrees
of freedom in an orthogonal matrix), while the dimension of the DM manifold
${\mathcal{M}}_{\rm DM}$ can be shown to be
$N_{I}N_{A}+N_{I}N_{E}+N_{A}N_{E}$. The discrepancy comes from the fact that
rotations that mix orbitals of the same class do not affect the energy. In
mathematical terms, this can be formulated as follows: the DM manifold
${\mathcal{M}}_{\rm DM}$ can be identified with the quotient of the MO
manifold $O({\mathcal{N}}_{b})$ by the group $O(N_{I})\times O(N_{A})\times
O(N_{E})$. This identification has very practical consequences on the design
of direct optimization algorithms, as will be seen in Section 4.

The geometric structure described above corresponds to a well known fact in
quantum chemistry. In direct optimization implementations, changes in the
orbitals are parameterized via a rotation matrix

$U=e^{\kappa},$ (12)

where $\kappa\in{\mathbb{R}}^{{\mathcal{N}}_{\rm b}\times{\mathcal{N}}_{\rm
b}}$ is a skew-symmetric matrix with the following block structure:

$\kappa=\begin{pmatrix}0&\kappa_{IA}&\kappa_{IE}\\\
-\kappa^{T}_{IA}&0&\kappa_{AE}\\\
-\kappa^{T}_{IE}&-\kappa^{T}_{AE}&0\end{pmatrix}$ (13)

The vanishing diagonal blocks, that would mix orbitals belonging to the same
class, are the practical translation of the quotient process mentioned above.
We further note that the map $\kappa\mapsto Ce^{\kappa}$, with $\kappa$ as in
Eq. 13 provides a non-redundant local parametrization of the quotient manifold
$O({\mathcal{N}}_{b})/(O(N_{I})\times O(N_{A})\times O(N_{E}))$.

In standard quantum chemistry direct optimization implementations, a sequence
of MO coefficients $\\{C^{(k)}\\}_{k=0}^{N_{it}}$ is generated starting from
an initial guess $C^{(0)}$ such that the sequence of energies
$\\{E^{(k)}\\}_{k=0}^{N_{it}}$ is non increasing and hopefully converges to
the ground state energy. The passage from $C^{(k)}$ to $C^{(k+1)}$ is obtained
by

$C^{(k+1)}=C^{(k)}e^{\kappa^{(k)}}$ (14)

where $\kappa^{(k)}$ is the result of some optimization procedure (e.g.,
steepest descent, or Levenberg-Marquardt second order optimization). Eq. 14
amounts to changing the center of the local parameterization of the quotient
manifold. For the steepest descent and Newton optimization methods, the
calculation of $\kappa^{(k)}$ relies only on information relative to the point
$C^{(k)}$, i.e., it makes no use of the history. Employing an optimization
method that does, e.g., non-linear conjugate gradient (CG) or quasi-Newton
methods, comes with a complication. Let us consider non-linear CG as an
example. In the (flat) vector space ${\mathbb{R}}^{n}$, the CG descent
direction at a given iterate is computed by linearly combining the gradient at
the current iterate with the descent direct at the previous direction, see
Fig. 2 (left) and Eq. 18. On a Riemannian manifold, the gradient and descent
direction at a given iterate belong to the tangent space to the manifold at
this particular iterate, which changes from iteration to iteration. Therefore,
it is not possible to linearly combine tangent vectors at different points, as
required by CG, in an obvious way. The operation of correctly transferring a
vector quantity from the tangent space at a given point of the manifold to the
tangent space at another point of the manifold is called _transport_ , see Fig
2 (right) and Eq. 19.

&

$\scriptstyle x^{(k)}$$\mathcal{M}$$\scriptstyle x^{(k+1)}$$\scriptstyle
d^{(k-1)}$$\scriptstyle
d^{(k)}$$\scriptstyle-\nabla_{\mathfrak{g}}E\left(x^{(k)}\right)$$\scriptstyle\mathcal{T}_{t_{k-1}d^{(k-1)}}d^{(k-1)}$

Figure 1: Conjugate gradient algorithm in ${\mathbb{R}}^{n}$ (left).
Riemannian conjugate gradient (RCG) algorithm on a Riemannian manifold
$\mathcal{M}$ (right).

For RHF and RKS, this problem is not apparent because the optimization takes
place in a Grassmann manifold, for which the parallel transport map is trivial
in the right parameterization (see Eq. 23). This is not the case for the flag
manifolds on which ROHF and CASSCF optimization problems are set (see Eq. 27).

Given $(\Pi^{\rm I},\Pi^{\rm A})\in{\mathcal{M}}_{\rm DM}$, there exists a
unique (up to an irrelevant global phase) normalized ROHF wavefunction
$\Phi_{(\Pi^{\rm I},\Pi^{\rm A})}^{\rm ROHF}$ with maximal spin polarization
$S=S_{z}=\frac{N_{A}}{2}$ associated with $(\Pi^{\rm I},\Pi^{\rm A})$: it is
the Slater determinant whose spin-1-RDM in the AO basis $(\chi_{\mu})$ is
given by Eq. 4. The ROHF energy functional

$E^{\rm ROHF}(\Pi^{\rm I},\Pi^{\rm A}):=\langle\Phi_{(\Pi^{\rm I},\Pi^{\rm
A})}^{\rm ROHF}|\widehat{H}_{N}|\Phi_{(\Pi^{\rm I},\Pi^{\rm A})}^{\rm
ROHF}\rangle$ (15)

is therefore a well-defined function of $(\Pi^{\rm I},\Pi^{\rm A})$, and in
fact a quadratic function in $\Pi^{\rm I}$ and $\Pi^{\rm A}$. From a
geometrical point of view, the ROHF problem is therefore a smooth optimization
problem on a flag manifold, for which the energy is quadratic in the density-
matrix formalism.

CASSCF can be also seen as an optimization problem on a flag manifold. In the
spin-collinear approximation, the corresponding CASSCF energy functional can
be written as E^CAS(Π^I,Π^A) = min_Ψ∈W_Π^I,Π^A^CAS ⟨Ψ—^H_N—Ψ⟩, with
W_Π^I,Π^A^CAS := { Ψ s.t. ∥ Ψ∥=1, Π^I ≤P^Ψ,σ ≤Π^I+Π^A, Tr(P^Ψ,σ)=N_I+m_σ,
σ=α,β}. Recall that for $A,B\in{\mathbb{R}}^{{\mathcal{N}}_{\rm
b}\times{\mathcal{N}}_{\rm b}}_{\rm sym}$, $A\leq B$ means that $X^{T}AX\leq
X^{T}BX$ for all $X\in{\mathbb{R}}^{{\mathcal{N}}_{\rm b}}$.

A very appealing feature of quotient manifolds is that if closed form
expressions for parallel transport and geodesics on ${\mathcal{M}}$ are known,
then closed form expressions for parallel transports on ${\mathcal{M}}/G$ can
be derived from the ones on ${\mathcal{M}}$. The manifold
$O({\mathcal{N}}_{\rm b})$ is in fact a Lie group. For this reason, closed
form expressions for parallel transport and geodesics on $O({\mathcal{N}}_{\rm
b})$ can be constructed from the exponential map.

##  3 Optimization on Riemannian manifolds

Riemannian optimization (i.e. optimization on manifold endowed with a
Riemannian metric) is a major field of computational mathematics with many
applications in various areas of science and technology. Several Riemannian
optimization libraries are available, in which the most common Riemannian
optimization methods are implemented. One of the advantages of using an
optimization library is that this allows one to test and compare many
different optimization methods with limited development effort. For
optimization in the flat space ${\mathbb{R}}^{n}$, the user of an optimization
library is just asked to provide the code returning the value of the function
and its gradient at an input point $x\in{\mathbb{R}}^{n}$ (and possibly also,
for some methods, a preconditioner and/or the Hessian at $x$). For
optimization on a Riemannian manifold ${\mathcal{M}}$, the user is asked to
provide four pieces of codes returning respectively:

  1. 1.

the value of the function and its Riemannian gradient at an input point
$x\in{\mathcal{M}}$, (and possibly also, for some methods, a preconditioner
and/or the Riemannian Hessian at $x$);

  2. 2.

the value of $\mathfrak{g}_{x}(p_{x},q_{x})\in{\mathbb{R}}$, where
$\mathfrak{g}_{x}$ is the Riemannian metric, for an input point
$x\in{\mathcal{M}}$ and two tangent vectors $p_{x},q_{x}\in
T_{x}{\mathcal{M}}$ at point $x$;

  3. 3.

the value of ${\mathcal{R}}_{x}(p_{x})\in{\mathcal{M}}$, where ${\mathcal{R}}$
is the chosen retraction, for an input point $x\in{\mathcal{M}}$ and a tangent
vector $p_{x}\in T_{x}{\mathcal{M}}$ at point $x$;

  4. 4.

the value of ${\mathcal{T}}_{p_{x}}q_{x}\in
T_{{\mathcal{R}}_{x}(p_{x})}{\mathcal{M}}$, where ${\mathcal{T}}$ is the
chosen transport, for an input point $x\in{\mathcal{M}}$ and two tangent
vectors $p_{x},q_{x}\in T_{x}{\mathcal{M}}$ at point $x$.

Let us first recall the role of $\mathfrak{g}$, ${\mathcal{R}}$ and
${\mathcal{T}}$ in Riemannian optimization algorithms, and illustrate these
concepts on the simple example of optimization on the orthogonal group
$O({\mathcal{N}}_{b})$. It is well-known that the tangent space to
$O({\mathcal{N}}_{b})$ at some $C\in O({\mathcal{N}}_{b})$ is given by
T_CO(N_b) = { CA, A ∈R^N_b ×N_b_antisym }, where
$A\in{\mathbb{R}}^{{\mathcal{N}}_{b}\times{\mathcal{N}}_{b}}_{\rm antisym}$ is
the vector space of ${\mathcal{N}}_{b}\times{\mathcal{N}}_{b}$ real
antisymmetric matrix. The Frobenius inner product ⟨M,N ⟩_F :=Tr(M^TN) =
∑_μ,ν=1^N_b M_μν N_μν on
${\mathbb{R}}^{{\mathcal{N}}_{b}\times{\mathcal{N}}_{b}}$ induces a Riemannian
metric on $O({\mathcal{N}}_{b})$ defined by g_C(CA,CA’)=Tr(A^T A’)=-Tr(AA’)
for all  C ∈O(N_b)  and  CA,CA’ ∈T_CO(N_b). Consider a smooth function
$E:{\mathbb{R}}^{{\mathcal{N}}_{b}\times{\mathcal{N}}_{b}}\to{\mathbb{R}}$.
The gradient of $E$ at some point
$C\in{\mathbb{R}}^{{\mathcal{N}}_{b}\times{\mathcal{N}}_{b}}$ is the unique
matrix $\nabla E(C)\in{\mathbb{R}}^{{\mathcal{N}}_{b}\times{\mathcal{N}}_{b}}$
such that E(C+δC) = E(C) + ⟨∇E(C),δC ⟩_F + o(δC). If $C\in
O({\mathcal{N}}_{b})$, the Riemannian gradient of $E$ at $C$ is the matrix
$\nabla_{\mathfrak{g}}E(C)\in T_{C}O({\mathcal{N}}_{b})$ obtained by
orthogonally projecting $\nabla E(C)$ on $T_{C}O({\mathcal{N}}_{b})$ for the
Frobenius inner product. Its expression is given by ∇_gE(C) = 12 C ( C^T ∇E(C)
- ∇E(C)^T C ).

A retraction ${\mathcal{R}}$ on a manifold ${\mathcal{M}}$ is a map
${\mathcal{R}}:T{\mathcal{M}}\to{\mathcal{M}}$ such that for all
$x\in{\mathcal{M}}$ the restriction
${\mathcal{R}}_{x}:T_{x}{\mathcal{M}}\to{\mathcal{M}}$ of ${\mathcal{R}}$ to
$T_{x}{\mathcal{M}}$ satisfies for $p_{x}\in T_{x}{\mathcal{M}}$

${\mathcal{R}}_{x}(p_{x})=x+p_{x}+o(p_{x}).$ (16)

Among other things (see below), retractions are used to map straight lines, or
more generally paths, drawn on the vector space $T_{x}{\mathcal{M}}$ onto
paths drawn on the curved manifold ${\mathcal{M}}$. As a matter of example,
the iterates of the fixed-step gradient (also called steepest descent)
algorithm are defined by

$d^{(k)}=-\nabla_{\mathfrak{g}}E(x^{(k)}),\quad
x^{(k+1)}={\mathcal{R}}_{x^{(k)}}(td^{(k)}),\qquad\mbox{(fixed-step steepest
descent)}$ (17)

for a chosen fixed step $t>0$. In words, starting from a point
$x^{(k)}\in{\mathcal{M}}$, the descent direction $d^{(k)}$ is chosen equal to
the opposite of the gradient, which is the steepest descent direction for
infinitesimal length steps, a step $td^{(k)}$ is made in this direction, and
finally, the vector $td^{(k)}\in T_{x^{(k)}}{\mathcal{M}}$ is mapped back to a
point of the manifold thanks to the retraction (see Fig. 2).

$x^{(k)}$$t_{k}d^{(k)}$$T_{x}\mathcal{M}$$\mathcal{M}$$x^{(k+1)}=\mathcal{R}_{x^{(k)}}\left(t_{k}d^{(k)}\right)$$-\nabla_{\mathfrak{g}}E$$\mathcal{R}_{x^{(k)}}$$x^{(k+1)}$
Figure 2: Riemannian steepest descent (RSD) algorithm.

It follows from (16) that in the limit of small step lengths, we have
x^(k+1)=R_x^(k)(t d^(k)) = x^(k) + t d^(k) + o(t) where the remainder term
$o(t)$ can be interpreted as a correction due to the curvature of the manifold
${\mathcal{M}}$.

Among all possible retractions on a Riemannian manifold ${\mathcal{M}}$, one
is canonical: it is the one defined from the geodesics, called the exponential
map, and denoted by ${\rm Exp}$. In the case of $O({\mathcal{N}}_{b})$, the
exponential map has a simple closed expression and is related to the usual
exponential of matrices:

${\rm Exp}_{C}(CA)=Ce^{A},\quad
A\in{\mathbb{R}}^{{\mathcal{N}}_{b}\times{\mathcal{N}}_{b}}_{\rm
antisym}\quad\mbox{(exponential map on $O({\mathcal{N}}_{b})$).}$

In Riemannian optimization, vector transports are used in particular to
combine the descent directions of previous iterates. Let us further elaborate
on this point, that was qualitatively discussed in Section 2. In the standard
conjugate gradient algorithm in ${\mathbb{R}}^{n}$, the descent direction
$d^{(k)}$ at iterate $x^{(k)}$ is a linear combination of $-\nabla E(x^{(k)})$
and $d^{(k-1)}$, the previous descent direction:

$\left|\begin{array}[]{l}d^{(k)}=-g^{(k)}+\beta_{k}d^{(k-1)},\\\
x^{(k+1)}=x^{(k)}+t_{k}d^{(k)},\end{array}\right.\quad\mbox{(nonlinear CG
algorithm in ${\mathbb{R}}^{n}$)}$ (18)

with $g^{(k)}=\nabla E(x^{(k)})$ and either

$\displaystyle\beta_{k}^{\rm FR}:=\frac{\|g^{(k)}\|^{2}}{\|g^{(k-1)}\|^{2}}$
$\displaystyle\mbox{(Fletcher-Reeves)},$ $\displaystyle\beta_{k}^{\rm
PR}:=\frac{{g^{(k)}}^{T}(g^{(k)}-g^{(k-1)})}{\|g^{(k-1)}\|^{2}}$
$\displaystyle\mbox{(Polak-Ribi\\`{e}re)},$ $\displaystyle\beta_{k}^{\rm
HS}:=\frac{{g^{(k)}}^{T}(g^{(k)}-g^{(k-1)})}{(g^{(k)}-g^{(k-1)})^{T}d^{(k)}}$
$\displaystyle\mbox{(Hestenes-Stiefel)}.$

The step length $t_{k}\in{\mathbb{R}}$ is obtained by a line search technique
such as Armijo, Wolfe, or Hager-Zhang [34] linesearch.

This idea cannot be directly used in optimization on manifolds, because
$-\nabla E(x^{(k)})$ and $d^{(k-1)}$ belong to different vector spaces, namely
$T_{x^{(k)}}{\mathcal{M}}$ and $T_{x^{(k-1)}}{\mathcal{M}}$ respectively.
Before being combined with $-\nabla E(x^{(k)})$ to form the new descent
direction $d^{(k)}$, the vector $d^{(k-1)}$ must be transported from
$T_{x^{(k-1)}}{\mathcal{M}}$ to $T_{x^{(k)}}{\mathcal{M}}$, using a transport
map ${\mathcal{T}}$. A transport map takes as input two vectors $p_{x},q_{x}$
of the tangent space $T_{x}{\mathcal{M}}$ at point $x$, and returns a vector
${\mathcal{T}}_{p_{x}}q_{x}$ of the tangent space at point
${\mathcal{R}}_{x}(p_{x})$ (compatibility condition with the retraction
${\mathcal{R}}$). The map $(p_{x},q_{x})\mapsto{\mathcal{T}}_{p_{x}}q_{x}$ is
linear in the variable $q_{x}$, and satisfies the consistency relation
${\mathcal{T}}_{0}q_{x}=q_{x}$. We thus have

$\left|\begin{array}[]{l}d^{(k)}=-g^{(k)}+\beta_{k}{\mathcal{T}}_{t_{k-1}d^{(k-1)}}d^{(k-1)},\\\
x^{(k+1)}={\mathcal{R}}_{x^{(k)}}(t_{k}d^{(k)}),\end{array}\right.\quad\mbox{(Riemannian
CG algorithm)}$ (19)

with either

$\displaystyle\beta_{k}^{\rm
RFR}:=\frac{{\mathfrak{g}}_{x^{(k)}}(g^{(k)},g^{(k)})}{\mathfrak{g}_{x^{(k-1)}}(g^{(k-1)},g^{(k-1)})}$
$\displaystyle\mbox{(Riemannian Fletcher-Reeves)},$
$\displaystyle\beta_{k}^{\rm
RPR}:=\frac{\mathfrak{g}_{x^{(k)}}(g^{(k)},g^{(k)}-{\mathcal{T}}_{t_{k-1}d^{(k-1)}}(g^{(k-1)}))}{\mathfrak{g}_{x^{(k)}}(g^{(k-1)},g^{(k-1)})}$
$\displaystyle\mbox{(Riemannian Polak-Ribi\\`{e}re)},$
$\displaystyle\beta_{k}^{\rm
RHS}:=\frac{\mathfrak{g}_{x^{(k)}}(g^{(k)},g^{(k)}-{\mathcal{T}}_{t_{k-1}d^{(k-1)}}(g^{(k-1)}))}{\mathfrak{g}_{x^{(k)}}(g^{(k)},{\mathcal{T}}_{t_{k-1}d^{(k-1)}}(d^{(k-1)}))-g_{x^{(k-1)}}(g^{(k-1)},d^{(k-1)})}$
$\displaystyle\mbox{(Riemannian Hestenes-Stiefel)}.$

Together with $x^{(k)}={\mathcal{R}}_{x^{(k-1)}}(t_{k-1}d^{(k-1)})$, the fact
that ${\mathcal{T}}_{p_{x}}q_{x}\in T_{{\mathcal{R}}_{x}(p_{x})}{\mathcal{M}}$
(comptabilitly with the retraction) ensures that
${\mathcal{T}}_{t_{k-1}d_{(k-1)}}d^{(k-1)}$ belongs to
$T_{x^{(k)}}{\mathcal{M}}$.

Among all transports compatible with the exponential map ${\rm Exp}$
associated with the metric $\mathfrak{g}$, one is canonical: it is the
parallel transport associated with the Levi-Civita connection of the metric
$\mathfrak{g}$. For the example of $O({\mathcal{N}}_{b})$, this parallel
transport has an extremely simple form T_CA(CB) = Ce^A B (parallel transport
on $O({\mathcal{N}}_{b})$).

##  4 Optimization on Grassmann and flag manifolds

In RHF and RKS models, the state of the system is described by a point of the
Grassmann manifold Gr(N,N_b) ≅⏟{P ∈R^N_b ×N_b_sym  s.t.  P^2=P, Tr(P)=N }_DM
formalism ≅⏟O(N_b)/(O(N) ×O(N_b-N))_MO formalism. In the DM formalism, the
Grassmann manifold is parameterized by the matrix $P$ of the orthogonal
projector on the vector space spanned by the doubly-occupied MO. In the MO
formalism, it is represented by an orthogonal matrix $C\in
O({\mathcal{N}}_{b})$, the first $N$ columns of $C$ corresponding to the $N$
doubly-occupied orbitals, and the last ${\mathcal{N}}_{b}-N$ ones to the
virtual orbitals. The gauge invariance in the MO formulation is taken into
account by quotienting $O({\mathcal{N}}_{b})$ by the group $O(N)\times
O({\mathcal{N}}_{b}-N)$ (occupied-occupied and virtual-virtual rotations,
repsectiely).

Likewise, in the ROHF model and the outer CASSCF minimization problem, the
state is represented by a point of the flag manifold Flag(N_I,N_I+N_A,N_b)
≅M_DM ≅⏟O(N_b)/(O(N_I) ×O(N_A) ×O(N_E))_MO formalism.

In both cases, the MO formalism involves the quotient of the orthogonal group
$O({\mathcal{N}}_{b})$ by a closed subgroup ($O(N)\times
O({\mathcal{N}}_{b}-N)$ for RHF/RKS, $O(N)\times O(N_{I})\times O(N_{A})\times
O(N_{E})$ for ROHF/CASSCF). As a consequence, the closed form expressions for
the canonical retraction and parallel transport on $O({\mathcal{N}}_{b})$ can
be translated into closed form expressions for canonical retraction and
parallel transport on the quotient manifold [26, 30, 35, 36], leading to the
following formulae:

  * •

RHF/RKS setting:

$\displaystyle\mbox{tangent space at
$C$}\quad\cong\left\\{C\kappa,\quad\kappa=\begin{pmatrix}0&\kappa_{OV}\\\
-\kappa^{T}_{OV}&0\end{pmatrix}\right\\},$ (20)
$\displaystyle\mbox{metric:}\qquad\mathfrak{g}_{C}(C\kappa,C\kappa^{\prime})={\rm
Tr}(\kappa^{T}\kappa^{\prime})=2{\rm
Tr}(\kappa_{OV}^{T}\kappa_{OV}^{\prime}),$ (21)
$\displaystyle\mbox{exponential map (canonical
retraction):}\qquad{\mathcal{R}}_{C}(C\kappa)=Ce^{\kappa},$ (22)
$\displaystyle\mbox{parallel
transport:}\qquad{\mathcal{T}}_{C\kappa}(C\kappa^{\prime})=Ce^{\kappa}\kappa^{\prime},$
(23)

  * •

ROHF/CASSCF setting:

$\displaystyle\mbox{tangent space at
$C$}\quad\cong\left\\{C\kappa,\quad\kappa=\begin{pmatrix}0&\kappa_{IA}&\kappa_{IE}\\\
-\kappa_{IA}^{T}&0&\kappa_{AE}\\\
-\kappa_{IE}^{T}&-\kappa_{AE}^{T}&0\end{pmatrix}\right\\},$ (24)
$\displaystyle\mbox{metric:}\qquad\mathfrak{g}_{C}(C\kappa,C\kappa^{\prime})={\rm
Tr}(\kappa^{T}\kappa^{\prime})$ (25) $\displaystyle\mbox{exponential map
(canonical retraction):}\qquad{\mathcal{R}}_{C}(C\kappa)=Ce^{\kappa},$ (26)
$\displaystyle\mbox{parallel
transport:}\qquad{\mathcal{T}}_{C\kappa}(C\kappa^{\prime})=Ce^{\kappa}e^{-\phi_{\kappa}}(\kappa^{\prime}),$
(27)

where $\phi_{\kappa}:\mathfrak{K}\to\mathfrak{K}$ is the linear operator on
the vector space K= { κ= (0 κIAκIE\- κIAT0 κAE\- κIET\- κAET0 ) } defined by

$\displaystyle\phi_{\kappa}(\kappa^{\prime})$ $\displaystyle=\frac{1}{2}{\rm
Proj}_{\mathfrak{K}}\left([\kappa,\kappa^{\prime}]\right)$
$\displaystyle=\frac{1}{2}\begin{pmatrix}0&-\kappa_{IE}[\kappa_{AE}^{\prime}]^{T}+\kappa_{IE}^{\prime}\kappa_{AE}^{T}&\kappa_{IA}\kappa_{AE}^{\prime}-\kappa_{IA}^{\prime}\kappa_{AE}\\\
\kappa_{AE}^{\prime}\kappa_{IE}^{T}-\kappa_{AE}[\kappa_{IE}^{\prime}]^{T}&0&-\kappa_{IA}^{T}\kappa_{IE}^{\prime}+[\kappa_{IA}^{\prime}]^{T}\kappa_{IE}\\\
-[\kappa_{AE}^{\prime}]^{T}\kappa_{IA}^{T}+\kappa_{AE}^{T}[\kappa_{IA}^{\prime}]^{T}&[\kappa_{IE}^{\prime}]^{T}\kappa_{IA}-\kappa_{IE}^{T}\kappa_{IA}^{\prime}&0\end{pmatrix},$

and

$e^{-\phi_{\kappa}}=\sum_{n=0}^{+\infty}\frac{(-1)^{n}}{n!}\underbrace{(\phi_{\kappa}\circ\cdots\circ\phi_{\kappa})}_{n\;\rm
times}.$ (28)

In computational codes, it is convenient to represent a tangent vector
$C\kappa$ by the block
$\kappa_{OV}\in{\mathbb{R}}^{N\times({\mathcal{N}}_{b}-N)}$ for RHS/RKS, and
the blocks $(\kappa_{IA},\kappa_{IE},\kappa_{AE})\in{\mathbb{R}}^{N_{I}\times
N_{A}}\times{\mathbb{R}}^{N_{I}\times N_{E}}{\mathbb{R}}^{N_{A}\times N_{E}}$
for ROHF/CASSCF. It follows from Eq. 23, that in this representation the
parallel transport for RHF/RKS is the identity operator. This is not the case
for ROHF/CASSCF where the transport of the block matrix $\kappa^{\prime}$ is
done by the map $e^{-\phi_{\kappa}}(\kappa^{\prime})$, which transforms and
mixes the IA/IE/AE blocks of $\kappa^{\prime}$. Let us note however that in
the special case when the transported vector $C\kappa^{\prime}$ is collinear
to the vector $C\kappa$ along which it is transported, then the transport
formula Eq. 27 dramatically simplifies. Indeed, we then have
$[\kappa,\kappa^{\prime}]=0$, and therefore
$e^{-\phi_{\kappa}}(\kappa^{\prime})=\kappa^{\prime}$. This occurs for the
Riemannian conjugate gradient method (see Eq. 19), but not for quasi-Newton
methods such as BFGS.

##  5 Numerical Results

In this section, we analyze the performance of Riemannian optimization
algorithms for solving the ROHF and CASSCF minimization problems for a few
selected test cases. Let us first provide some implementation details.

#### General implementation.

We focus specifically on Riemannian steepest descent (RSD), nonlinear
conjugate gradient (RCG) and low-memory Broyden-Fletcher-Goldfarb-Shanno
(R-LBFGS) methods, all endowed with preconditioning. We refer to Refs [35, 36]
for general introductions to Riemannian optimization methods. Our code is
structured as follows: first, we implemented the RSD, RCG, and R-LBFGS
optimization routines in the MO formalism within a Julia [37] package which is
then interfaced with PySCF [38] for ROHF and CFOUR [39] for CASSCF
calculations. These software handle the operations specific to the ROHF and
CASSCF models, including the generation of AO basis sets and initial guess
MOs, the computation of electronic integrals, and the evaluation of energies
and Frobenius gradients.

In our Julia package, we use for all methods the exponential retraction (26)
and parallel transport (27) as outlined in the previous section. For parallel
transport, the exponential operator is computed by truncating the series (28)
so that the Frobenius norm of the last term falls below numerical precision.

Our implementation of RCG is based on Algorithm 1 in Boumal et al. [40] with
Polak-Ribière coefficient $\beta^{\rm RPR}$ as above. For R-LBFGS, we
implemented Algorithm 2 in Huang et al. [41]. For CASSCF, we use the inverse
diagonal of the Hessian as preconditioner. In the case of ROHF, we tested two
different preconditioners. The first one is the modified inverse diagonal
Hessian as discussed in Ref [21]. The other is detailed in Appendix. Our
results are presented for the second choice of preconditioner that showcased
the best performance for our test case.

All methods use Hager-Zhang [34] linesearch as implemented in the
LineSearches.jl [42] Julia package. Computations are considered to have
reached convergence when the Frobenius norm of the Riemannian gradient reaches
$10^{-5}$. Comprehensive implementation details are available in our publicly
accessible GitHub repository.111https://github.com/LaurentVidal95/ROHFToolkit

#### Details specific to the R-LBFGS implementation.

At each iteration, the R-LBFGS method constructs an approximation $B$ of the
inverse Hessian using a certain number $m$ of vectors stored in memory from
previous iterations, through an iterative procedure [41]. In addition to
preconditioning, it is important to note that the performance of R-LBFGS is
influenced by the selection of the maximum depth $m_{\rm max}$, the initial
guess $B_{0}$ for the approximate inverse Hessian in the iterative process and
the choice of restart strategies, which determines the iterations at which the
history is reset. For both ROHF and CASSCF, we define $B_{0}=\gamma\,{\rm Id}$
with $\gamma$ as in Ref [41]. If $P$ is the preconditioner, the preconditioned
version of R-LBFGS is obtained by replacing $B_{0}$ with $\hat{\gamma}P$, with
$\hat{\gamma}$ as in Ref [43].

We experimented two restart strategies which depend on the preconditioning.
For the first one, called “dynamic” R-LBFGS, the diagonal Hessian used for
preconditioning is updated at each iteration. The history is reset whenever
the direction obtained from the R-LBFGS quasi-newton system is not a descent
direction. For the second method called “fixed” R-LBFGS, we use the same
preconditioner $P={\rm diag}({\rm Hess}^{(0)})^{-1}$, corresponding to the
inverse diagonal Hessian for the guess orbitals, at each iteration until the
inverse diagonal Hessian at current point, ${\rm diag}({\rm
Hess}^{(k)})^{-1}$, deviates too much from $P$. When this happens, the history
is reset, and the procedure is reinitialized with $P={\rm diag}({\rm
Hess}^{(k)})^{-1}$. When using the preconditioner described in appendix for
ROHF, we applied the dynamic strategy.

###  5.1 ROHF

For ROHF we tested the three aforementioned Riemannian optimization methods on
Ti2O4 in its $D_{\rm 2h}$ geometry, using Dunning’s cc-pVTZ basis set [44,
45]. This system is employed as a template for addressing SCF convergence
issues222See https://www.scm.com/doc/ADF/Examples/SCF_Ti2O4.html in the
Amsterdam Density Functional (ADF) quantum-chemistry package [46]. In order to
compare the performance of Riemannian algorithms in different convergence
regimes, calculations were started from both a core initial guess (Fig. 3) and
a guess closer to a minimum (Fig. 4). The second guess is obtained by a
standard SCF+DIIS method for ROHF, with Guest and Saunders coefficients [47],
stopped when the Frobenius norm of the gradient reaches $10^{-1}$.

|
---|---
Figure 3: Calculations from core initial guess. On the left, energy difference with respect to the converged energy along the iterations. On the right, Frobenius norm of the Riemannian gradient along the iterations. Only the first 500 iterations of RSD are shown on the graph for readability. |
---|---
Figure 4: Calculations from a better initial guess, 1 $E_{h}$ from the
expected energy. On the left, energy difference with respect to the converged
energy along the iterations. On the right, Frobenius norm of the Riemannian
gradient along the iterations.

In both cases, all three methods provide stable convergence toward a local
minimum of the energy. In the optimal scenario, a finely tuned SCF+DIIS method
outperforms the Riemannian optimization methods we have tested. However, the
performance and stability of SCF routines for ROHF are notoriously sensitive
to the choice of method and acceleration parameters, as illustrated in Fig. 5.
On the other hand, the direct minimization methods described in this paper
have the advantage of offering robust convergence, which is a valuable feature
in terms of user’s time and effort.

We were unfortunately not able to make a direct comparison to the GDM
algorithm [20] due to our lack of access to the code or to the fine details of
the implementation. Nevertheless, a simple-minded test performed with the free
trial version of Q-Chem [48] using default parameters showed that GDM and
GDM+DIIS exhibit similar performances to our RCG implementation on Ti2O4.

Figure 5: Comparison of standard SCF+DIIS methods for ROHF as listed in Ref
[47], for Ti2O4 in cc-pVTZ basis set, starting from our good initial guess
about 1 Hartree away from the expected energy. The vertical axis shows the
Frobenius norm of the Riemannian gradient along the iterations. Only two
methods yield convergence.

An important point that needs to be discussed here is what minimum the various
algorithms converged to. Our proof-of-concept implementation does not enforce
point-group symmetry, so our calculations were performed in the $C_{\rm 1}$
group.

Our calculations converged to two different minima, one at
$-1996.191285~{}E_{h}$, which was systematically obtained when starting from
the core guess, and one at $-1996.179398~{}E_{h}$, obtained when starting from
the better guess. We also looked at the lowest triplets for each Irrep
enforcing symmetry using a quadratically convergent ROHF implementation [21],
and we determined that the lowest triplet is the $B_{\rm 2u}$ state at
$-1996.142005$ $E_{h}$. The stability analysis of such solution reveals
however that a lower energy, symmetry-broken solution exists. We therefore
conclude that the two solutions found with our Riemannian optimization
algorithms are two lower-energy symmetry-broken solutions. Whether a lower-
energy, symmetry broken-solution is desirable or not depends on the aims of
the study, and ultimately on the user; however, we note that our algorithms
can be generalized to enforce point group symmetry, which we plan to do in the
future.

###  5.2 CASSCF

The CASSCF method has been tested by running calculations on a subgroup of the
benchmark set used in Ref. [49] and Ref. [50] using Pople’s 6-31G* basis set
[51]. Convergence properties of direct minimization algorithms were compared
against two well-established CASSCF optimization algorithms namely Super CI
(SCI) [52, 53, 54, 55, 56] and the norm-extended optimization (NEO) [57, 58],
the latter one being a genuine second-order algorithm. All computations were
carried out using two different choices, to simulate, respectively, a
troublesome scenario where the calculation starts relatively far away from the
converged result, and an ideal starting point that should be close to the
final minimum. For the former scenario, we use canonical restricted Hartree-
Fock (RHF) orbitals, while as a good starting point we exploit unrestricted
natural orbitals (UNO) [59, 60]

We report in Tab. 1 for each algorithm the average number of iterations
required to converge CASSCF. The number of iterations for each tested system
is reported in the supporting information. As expected, the values related to
the RHF guess are systematically higher than the ones related to the UNO
guess. Moreover, we notice that the numbers related to the RHF guess show a
high variability as indicated by the large standard deviation, thus being
strongly system dependent. The average number of iterations for all direct
minimization methods with the exception of RSD is comparable with and in some
cases outperforms the ones of SCI.

Algorithm | $\langle\text{It.}\rangle^{\text{RHF}}$ | $\sigma^{\text{RHF}}$ | $\langle\text{It.}\rangle^{\text{UNO}}$ | $\sigma^{\text{UNO}}$
---|---|---|---|---
RSD | 115.9 | 112.9 | 32.3 | 12.5
RCG | 31.9 | 11.6 | 15.4 | 3.0
R-LBFGS(dyn) | 37.9 | 11.3 | 19.2 | 3.5
R-LBFGS(fix) | 35.6 | 15.2 | 19.1 | 3.7
SCI(DIIS) | 38.4 | 26.0 | 14.9 | 6.2
SCI(no DIIS) | 61.5 | 25.1 | 21.6 | 9.9
NEO | 12.2 | 1.8 | 5.1 | 0.5
Table 1: Average number of iterations ($\langle\text{It.}\rangle$) and
standard deviation ($\sigma$) for each tested algorithm starting with two
different guess orbitals, namely restricted Hartree-Fock (RHF) and
unrestricted natural orbitals (UNO).

We conclude this section by looking more in detail at one example, namely,
pyridine using a standard CAS(6,6) wavefunction. In Fig. 6, we compare the
convergence behavior of the R-LBFGS as implemented in the present study with a
naive implementation that simply translates the gradient from previous points,
without parallel transport. In both cases, we start from canonical orbitals.
We note that while both implementations get stuck for a while on a plateau,
the R-LBFGS overcomes it in a few iterations and then converges smoothly. On
the contrary, the naive LBFGS implementation exhibits worse convergence, with
very small and even positive slope steps at the beginning of the optimization.
This demonstrates the importance of properly accounting for parallel
transport.

Figure 6: On the left, energy difference with respect to the converged energy
along the iterations for a proper R-LBFGS implementation that makes use of the
parallel transport (green curve) and a more basic one (orange curve). On the
right, Frobenius norm of the gradient along the iterations.

Another important point concerns the actual solution obtained. When starting
from a poor guess such as the one given by canonical orbitals, without any
kind of manual selection for the active space, the optimizer may get stuck in
local minima that are ultimately related by orbitals swapping between the
inactive and active domains. Using once again pyridine as a test case, we
observe that all the algorithms converged to the same local minimum
(-246.766857 $E_{h}$) with the exception of SCI that converged to a higher
minimum (-246.756489 $E_{h}$).

These two minima are characterized by different converged active orbitals.
This difference can be assessed simply by visual inspection or by checking the
singular values of the difference between the active part of the one-body
density matrix in the AO basis for the two calculations. If the converged
active orbitals were (almost) the same, we would observe (almost) vanishing
singular values. On the contrary, as depicted in Fig. 7, two orbitals are
completely different between the two results.

Figure 7: First ten singular values of the difference matrix between the
active AO-based one-body density matrices stemming from two calculations that
reached different local minima.

We also note that by manually tuning the DIIS parameters, we managed to
achieve convergence to the lower-energy solution using SCI. Again, this
underlines the robustness of the Riemannian optimization algorithms, which is
valuable in terms of user time and effort.

##  6 Conclusions and perspectives

In this contribution, we have explored the use of Riemannian optimization
methods on the flag manifold to optimize restricted-open and complete active
space self-consistent field wavefunctions. After discussing the geometry of
the problem, we have reviewed the general aspects of Riemannian optimization
and its application to the aforementioned chemical problem. We have then
compared various algorithms to traditional ones. The Riemannian optimization
methods illustrated in this work all show robust convergence properties, and
do so without requiring the user to finely tune the parameters that control
the optimization. Even in the naive implementation presented here, they
demonstrate that they can be competitive with other traditional
implementations in terms of number of iterations, and thus overall
computational cost. Nevertheless, this is just a proof-of-concept study, for
which several further developments are required. First, the overall
underwhelming performance of the Riemannian Quasi-Newton L-BFGS method for
CASSCF, which is expected to outperform conjugate gradient, as observed for
ROHF when starting from a good guess (see Fig. 4), can be explained by the
basic preconditioner and initialization of the inverse Hessian we use. Finding
better preconditioners and inverse Hessian approximations is not a
straightforward task, and requires further attention. We also have not
investigated optimal parameters, including exploring different line-search and
restarting strategies, which would require to run extensive numerical tests,
but that can greatly improve the overall performance of the methods.

In conclusion, we believe that Riemannian optimization is a valuable addition
to the SCF optimization toolbox for ROHF and CASSCF, and that further
exploration of the use of such techniques is worthy of attention.

## Acknowledgements

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation program (grant
agreement EMC2 No. 810367). F.L. and T.N. acknowledge financial support from
ICSC-Centro Nazionale di Ricerca in High Performance Computing, Big Data, and
Quantum Computing, funded by the European Union – Next Generation EU – PNRR,
Missione 4 Componente 2 Investimento 1.4.

## References

  * [1] D. R. Hartree.  The calculation of atomic structures.  Wiley, London, 1957.
  * [2] W. Kohn and L. J. Sham.  Self-consistent equations including exchange and correlation effects.  Phys. Rev., 140:A1133–A1138, Nov 1965.
  * [3] Hans-Joachim Werner.  Matrix-Formulated Direct Multiconfiguration Self-Consistent Field and Multiconfiguration Reference Configuration-Interaction Methods.  Adv. Chem. Phys., 69:1–62, 1987.
  * [4] Ron Shepard.  The Multiconfiguration Self-Consistent Filed Method.  Adv. Chem. Phys., 69:63–200, 1987.
  * [5] Björn O. Roos.  The Complete Active Space Self-Consistent Field Method and its Applications in Electronic Structure Calculations.  Adv. Chem. Phys., 69:399–445, 1987.
  * [6] Gustavo E Scuseria and Henry F Schaefer III.  The optimization of molecular orbitals for coupled cluster wavefunctions.  Chemical physics letters, 142(5):354–358, 1987.
  * [7] C David Sherrill, Anna I Krylov, Edward FC Byrd, and Martin Head-Gordon.  Energies and analytic gradients for a coupled-cluster doubles model using variational brueckner orbitals: Application to symmetry breaking in o 4+.  The Journal of chemical physics, 109(11):4171–4181, 1998.
  * [8] C. C. J. Roothaan.  New developments in molecular orbital theory.  Rev. Mod. Phys., 23:69–89, Apr 1951.
  * [9] C. C.J. Roothaan.  Self-consistent field theory for open shells of electronic systems.  Rev. of Mod. Phys., 32:179–185, 1960.
  * [10] G.B. Bacskay.  A quadratically convergent Hartree—Fock (QC-SCF) method. application to closed shell systems.  Chemical Physics, 61(3):385–404, oct 1981.
  * [11] George B. Bacskay.  A quadratically convergent Hartree-Fock (QC-SCF) method. Application to open shell orbital optimization and coupled perturbed Hartree-Fock calculations.  Chem. Phys., 65:383–396, 1982.
  * [12] Péter Pulay.  Convergence acceleration of iterative sequences - the case of SCF iteration.  Chem. Phys. Lett., 73:393–398, 1980.
  * [13] P Pulay.  Improved SCF convergence acceleration.  J. Comput. Chem., 3:556–560, 1982.
  * [14] Tracy P. Hamilton and Peter Pulay.  Direct inversion in the iterative subspace (DIIS) optimization of open-shell, excited-state, and small multiconfiguration SCF wave functions.  J. Chem. Phys., 84:5728–5734, 1985.
  * [15] E. Cancès and C. Le Bris.  Can we outperform the DIIS approach for electronic structure calculations?  International Journal of Quantum Chemistry, 79(2):82–90, 2000.
  * [16] Eric Cancès.  SCF algorithms for HF electronic calculations, pages 17–43.  Springer, Berlin, Heidelberg, 2000.
  * [17] E. Cancès.  Self-consistent field algorithms for Kohn–Sham models with fractional occupation numbers.  The Journal of Chemical Physics, 114(24):10616–10622, 2001.
  * [18] Konstantin N. Kudin, Gustavo E. Scuseria, and Eric Cancès.  A black-box self-consistent field convergence algorithm: One step closer.  J. Chem. Phys., 116:8255–8261, 2002.
  * [19] Xiangqian Hu and Weitao Yang.  Accelerating self-consistent field convergence with the augmented Roothaan-Hall energy function.  The Journal of Chemical Physics, 132(5):054109, 02 2010.
  * [20] B.D. Dunietz, T. Van Voorhis, and M. Head-Gordon.  Geometric direct minimization of Hartree-Fock calculations involving open shell wavefunctions with spin restricted orbitals.  Journal of Theoretical and Computational Chemistry, 01(02):255–261, oct 2002.
  * [21] Tommaso Nottoli, Jürgen Gauss, and Filippo Lipparini.  A black-box, general purpose quadratic self-consistent field code with and without Cholesky Decomposition of the two-electron integrals.  Molecular Physics, 119(21-22):e1974590, November 2021.
  * [22] G. Chaban, M.W. Schmidt, and M.S. Gordon.  Approximate second order method for orbital optimization of SCF and MCSCF wavefunctions.  Theoretical Chemistry Accounts, 97(1):88–95, 1997.
  * [23] F. Neese.  Approximate second-order SCF convergence for spin unrestricted wavefunctions.  Chemical Physics Letters, 325(1-3):93–98, 2000.
  * [24] T. Van Voorhis and M. Head-Gordon.  A geometric approach to direct minimization.  Molecular Physics, 100(11):1713–1721, jun 2002.
  * [25] Takashi Tsuchimochi and Gustavo E Scuseria.  Communication: Rohf theory made simple.  The Journal of chemical physics, 133(14), 2010.
  * [26] A. Edelman, T.A. Arias, and S.T. Smith.  The geometry of algorithms with orthonormality constraints.  J. Matrix Anal. Appl., 20:303–353, 1998.
  * [27] E. Cancès and C. Le Bris.  On the convergence of SCF algorithms for the Hartree-Fock equations.  ESAIM: Mathematical Modelling and Numerical Analysis, 34(4):749–774, 2000.
  * [28] A. Levitt.  Convergence of gradient-based algorithms for the Hartree-Fock equations.  ESAIM: Mathematical Modelling and Numerical Analysis, 46(6):1321–1336, 2012.
  * [29] Éric Cancès, Gaspard Kemlin, and Antoine Levitt.  Convergence analysis of direct minimization and self-consistent iterations.  SIAM Journal on Matrix Analysis and Applications, 42(1):243–274, 2021.
  * [30] Ke Ye, Ken Sze-Wai Wong, and Lek-Heng Lim.  Optimization on flag manifolds.  Mathematical Programming, 194:621–660, 2022.
  * [31] É Polack, A. Mikhalev, G. Dusson, B. Stamm, and F. Lipparini.  An approximation strategy to compute accurate initial density matrices for repeated self-consistent field calculations at different geometries.  Mol. Phys., 118(19-20):e1779834, 2020.
  * [32] É Polack, G. Dusson, B. Stamm, and F. Lipparini.  Grassmann extrapolation of density matrices for born–oppenheimer molecular dynamics.  J. Chem. Theory Comput., 17(11):6965–6973, 2021.
  * [33] Federica Pes, Etienne Polack, Patrizia Mazzeo, Genevieve Dusson, Benjamin Stamm, and Filippo Lipparini.  A quasi time-reversible scheme based on density matrix extrapolation on the grassmann manifold for born–oppenheimer molecular dynamics.  The Journal of Physical Chemistry Letters, 14(43):9720–9726, 2023\.  PMID: 37879072.
  * [34] William W Hager and Hongchao Zhang.  Algorithm 851: Cg_descent, a conjugate gradient method with guaranteed descent.  ACM Transactions on Mathematical Software (TOMS), 32(1):113–137, 2006.
  * [35] P-A Absil, Robert Mahony, and Rodolphe Sepulchre.  Optimization algorithms on matrix manifolds.  Princeton University Press, 2008.
  * [36] Nicolas Boumal.  An introduction to optimization on smooth manifolds.  Cambridge University Press, 2023.
  * [37] J. Bezanson, A. Edelman, S. Karpinski, and V.B. Shah.  Julia: A fresh approach to numerical computing.  SIAM Review, 59(1):65–98, 2017.
  * [38] Qiming Sun, Xing Zhang, Samragni Banerjee, Peng Bao, Marc Barbry, Nick S Blunt, Nikolay A Bogdanov, George H Booth, Jia Chen, Zhi-Hao Cui, et al.  Recent developments in the pyscf program package.  The Journal of chemical physics, 153(2), 2020.
  * [39] Devin A Matthews, Lan Cheng, Michael E Harding, Filippo Lipparini, Stella Stopkowicz, Thomas-C Jagau, Péter G Szalay, Jürgen Gauss, and John F Stanton.  Coupled-cluster techniques for computational chemistry: The cfour program package.  The Journal of Chemical Physics, 152(21), 2020.
  * [40] Nicolas Boumal and P-A Absil.  Low-rank matrix completion via preconditioned optimization on the grassmann manifold.  Linear Algebra and its Applications, 475:200–239, 2015.
  * [41] Wen Huang, Kyle A Gallivan, and P-A Absil.  A broyden class of quasi-newton methods for riemannian optimization.  SIAM Journal on Optimization, 25(3):1660–1685, 2015.
  * [42] P Mogensen and A Riseth.  Optim: A mathematical optimization package for julia.  Journal of Open Source Software, 3(24), 2018.
  * [43] Hans De Sterck and Alexander JM Howse.  Nonlinearly preconditioned l-bfgs as an acceleration mechanism for alternating least squares with application to tensor decomposition.  Numerical Linear Algebra with Applications, 25(6):e2202, 2018.
  * [44] Thom H. Dunning, Jr.  Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen.  J. Chem. Phys., 1007:4572–4585, 1989.
  * [45] Rick A. Kendall, Thom H. Dunning, and Robert J. Harrison.  Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions.  J. Chem. Phys., 96:6796–6806, 1992.
  * [46] G t Te Velde, Friedrich Matthias Bickelhaupt, Evert Jan Baerends, C Fonseca Guerra, Stan JA van Gisbergen, Jaap G Snijders, and Tom Ziegler.  Chemistry with adf.  Journal of Computational Chemistry, 22(9):931–967, 2001.
  * [47] Boris N Plakhutin and Ernest R Davidson.  Canonical form of the hartree-fock orbitals in open-shell systems.  The Journal of Chemical Physics, 140(1), 2014.
  * [48] Yihan Shao, Zhengting Gan, Evgeny Epifanovsky, Andrew TB Gilbert, Michael Wormit, Joerg Kussmann, Adrian W Lange, Andrew Behn, Jia Deng, and Xintian Feng.  Advances in molecular quantum chemistry contained in the Q-Chem 4 program package.  Molecular Physics, 113(2):184–215, 2015.
  * [49] Filipe Menezes, Daniel Kats, and Hans-Joachim Werner.  Local complete active space second-order perturbation theory using pair natural orbitals (PNO-CASPT2).  J. Chem. Phys., 145:124115, 2016.
  * [50] Tommaso Nottoli, Jürgen Gauss, and Filippo Lipparini.  Second-order CASSCF algorithm with the Cholesky decomposition of the two-electron integrals.  J. Chem. Theory Comput., 17(11):6819–6831, 2021.
  * [51] Warren J Hehre, Robert Ditchfield, and John A Pople.  Self—consistent molecular orbital methods. xii. further extensions of gaussian—type basis sets for use in molecular orbital studies of organic molecules.  The Journal of Chemical Physics, 56(5):2257–2261, 1972.
  * [52] Björn O Roos.  The complete active space scf method in a Fock-matrix-based super-ci formulation.  Int. J. Quant. Chem., 18(S14):175–189, 1980.
  * [53] Per EM Siegbahn, Jan Almlöf, Anders Heiberg, and Björn O Roos.  The complete active space scf (casscf) method in a newton–raphson formulation with application to the hno molecule.  J. Chem. Phys., 74(4):2384–2396, 1981.
  * [54] Per Åke Malmqvist, Alistair Rendell, and Björn O Roos.  The restricted active space self-consistent-field method, implemented with a split graph unitary group approach.  J. Phys. Chem., 94(14):5477–5482, 1990.
  * [55] Christian Kollmar, Kantharuban Sivalingam, Benjamin Helmich-Paris, Celestino Angeli, and Frank Neese.  A perturbation-based super-ci approach for the orbital optimization of a casscf wave function.  J. Comp. Chem., 40(14):1463–1470, 2019.
  * [56] Celestino Angeli, Stefano Evangelisti, Renzo Cimiraglia, and Daniel Maynau.  A novel perturbation-based complete active space–self-consistent-field algorithm: Application to the direct calculation of localized orbitals.  J. Chem. Phys., 117(23):10525–10533, 2002.
  * [57] Hans Jørgen Aa Jensen and Poul Jørgensen.  A direct approach to second‐order MCSCF calculations using a norm extended optimization scheme.  J. Chem. Phys., 80(3):1204–1214, 1984.
  * [58] Hans Jørgen Aa Jensen and Hans Agren.  A Direct, restricted-step, second-order MC SCF program for large scale ab initio calculations.  Chem. Phys., 104:229–250, 1986.
  * [59] Peter Pulay and Tracy P. Hamilton.  UHF natural orbitals for defining and starting MC-SCF calculations.  J. Chem. Phys., 88:4926–4933, 1988.
  * [60] Zsuzsanna Tóth and Peter Pulay.  Comparison of Methods for Active Orbital Selection in Multiconfigurational Calculations.  J. Chem. Theory Comput., 16(12):7328–7341, 2020.

##  Appendix A A variant of the diagonal Hessian preconditioner for ROHF

Let $C\in O({\mathcal{N}}_{b})$ and $\kappa\in\mathfrak{K}$ such that
$C\kappa\in T_{C}O({\mathcal{N}}_{b})$. The ROHF Hessian applied to $C\kappa$
is given by

$\mathcal{L}_{C}(C\kappa)=C\begin{pmatrix}0&X&Y\\\ -X^{T}&0&Z\\\
-Y^{T}&-Z^{T}&0\end{pmatrix}$ (29)

where the matrices $X\in\mathbb{R}^{N_{I}\times N_{A}}$,
$Y\in\mathbb{R}^{N_{I}\times N_{E}}$ and $Z\in\mathbb{R}^{N_{A}\times N_{E}}$
are defined by

$\begin{array}[]{l}X=2(\kappa_{IA}(F_{I}-F_{A})_{AA}-(F_{I}-F_{A})_{II}\kappa_{IA})+\kappa_{IE}(2F_{I}-F_{A})_{EA}+(F_{I}-2F_{A})_{IE}\kappa_{AE}^{T}\\\
\hskip
170.71652pt+(2J(\lambda_{1})-K(\lambda_{1}))_{IA}+J(\lambda_{2})_{IA}\\\
Y=\kappa_{IA}(4F_{I}-2F_{A})_{AE}+4(\kappa_{IE}(F_{I})_{EE}-(F_{I})_{II}\kappa_{IE})-2((F_{I})_{IA}+(F_{A})_{IA})\kappa_{AE}\\\
\hskip
170.71652pt+4(2J(\lambda_{1})-K(\lambda_{1}))_{IE}+2(2J(\lambda_{2})-K(\lambda_{2}))_{IE}\\\
Z=\kappa_{IA}^{T}(2F_{I}-F_{A})_{IE}-2(F_{I}+F_{A})_{AI}\kappa_{IE}+4(\kappa_{AE}(F_{A})_{EE}-(F_{A})_{AA}\kappa_{AE})\\\
\hskip
170.71652pt+2(2J(\lambda_{1})-K(\lambda_{1}))_{AE}+2(2J(\lambda_{2})-K(\lambda_{2}))_{AE}.\end{array}$
(30)

In the above expression, we adopted the following conventions: the operators
$J$ and $K$ are the standard exchange and Coulomb operators, $F_{I}$ and
$F_{A}$ are the internal and active Fock matrices, defined for all
$\Pi_{I}=C_{I}C_{I}^{T}$ and $\Pi_{A}=C_{A}C_{A}^{T}$ by

$\begin{array}[]{l}F_{I}=\displaystyle
h+2J(\Pi_{I})+J(\Pi_{A})-K(\Pi_{I})-\frac{1}{2}K(\Pi_{A})\\\\[4.26773pt]
F_{A}=\displaystyle\frac{1}{2}(h+2J(\Pi_{I})+J(\Pi_{A})-K(\Pi_{I})-K(\Pi_{A}))\end{array}$

and the matrices
$\lambda_{1},\lambda_{2}\in{\mathbb{R}}^{{\mathcal{N}}_{b}\times{\mathcal{N}}_{b}}_{\rm
sym}$ are given by

$\lambda_{1}=\begin{pmatrix}0&\kappa_{IA}&\kappa_{IE}\\\
\kappa_{IA}^{T}&0&0\\\
\kappa_{IE}^{T}&0&0\end{pmatrix}\quad\mbox{and}\quad\lambda_{2}=\begin{pmatrix}0&-\kappa_{IA}&0\\\
-\kappa_{IA}^{T}&0&\kappa_{AE}\\\ 0&\kappa_{AE}^{T}&0\end{pmatrix}.$ (31)

Each block $X$, $Y$ and $Z$ can be decomposed as the sum of two terms, the
first one, denoted by $(\widetilde{X},\widetilde{Y},\widetilde{Z})$, being
simple to compute in terms of internal and active Fock operators, and the
second one, denoted by $(\Omega_{X},\Omega_{Y},\Omega_{Z})$, being more costly
to compute:

$\begin{array}[]{l}X=2(\kappa_{IA}(F_{I}-F_{A})_{AA}-(F_{I}-F_{A})_{II}\kappa_{IA})+\Omega_{X}:=\widetilde{X}+\Omega_{X}\\\
Y=4(\kappa_{IE}(F_{I})_{EE}-(F_{I})_{II}\kappa_{IE})+\Omega_{Y}:=\widetilde{Y}+\Omega_{Y}\\\
Z=4(\kappa_{AE}(F_{A})_{EE}-(F_{A})_{AA}\kappa_{AE})+\Omega_{Z}:=\widetilde{Z}+\Omega_{Z}\end{array}.$
(32)

In our implementation, we define a preconditioned direction $C\kappa_{\rm
prec}$ as a solution to the linear system

$\widetilde{L}_{C}(C\kappa_{\rm prec})=C\kappa.$ (33)

involving the approximate Hessian

$\widetilde{\mathcal{L}}_{C}(C\kappa)=C\begin{pmatrix}0&\widetilde{X}&\widetilde{Y}\\\
-\widetilde{X}^{T}&0&\widetilde{Z}\\\
-\widetilde{Y}^{T}&-\widetilde{Z}^{T}&0\end{pmatrix}.$ (34)

The advantage of formulation (32) is that the lowest eigenvalue of
$\widetilde{L}_{C}$ can be estimated with respect to $F_{I}$ and $F_{A}$,
which allows to apply a shift when $\widetilde{\mathcal{L}}$ is not positive
definite (which is expected when starting far from a minimum). In addition,
the system (33) reads as three Sylvester matrix equations, that can be solved
using standard LAPACK optimized routines.