# A force-based gradient descent method for ab initio atomic structure
relaxation

Yukuan Hu  State Key Laboratory of Scientific and Engineering Computing,
Academy of Mathematics and Systems Science, Chinese Academy of Sciences,
Beijing 100190, China  University of Chinese Academy of Sciences, Beijing
100049, China  Xingyu Gao<EMAIL_ADDRESS>Laboratory of Computational
Physics, Institute of Applied Physics and Computational Mathematics, Beijing
100088, China  Yafan Zhao  CAEP Software Center for High Performance Numerical
Simulation, Beijing 100088, China  Xin Liu<EMAIL_ADDRESS>State Key
Laboratory of Scientific and Engineering Computing, Academy of Mathematics and
Systems Science, Chinese Academy of Sciences, Beijing 100190, China
University of Chinese Academy of Sciences, Beijing 100049, China  Haifeng Song
<EMAIL_ADDRESS>Laboratory of Computational Physics, Institute of
Applied Physics and Computational Mathematics, Beijing 100088, China

###### Abstract

Force-based algorithms for ab initio atomic structure relaxation, such as
conjugate gradient methods, usually get stuck in the line minimization
processes along search directions, where expensive ab initio calculations are
triggered frequently to test trial positions before locating the next iterate.
We present a force-based gradient descent method, WANBB, that circumvents the
deficiency. At each iteration, WANBB enters the line minimization process with
a trial stepsize capturing the local curvature of the energy surface. The exit
is controlled by an unrestrictive criterion that tends to accept early trials.
These two ingredients streamline the line minimization process in WANBB. The
numerical simulations on nearly 80 systems with good universality demonstrate
the considerable compression of WANBB on the cost for the unaccepted trials
compared with conjugate gradient methods. We also observe across the board
significant and universal speedups as well as the superior robustness of WANBB
over several widely used methods. The latter point is theoretically
established. The implementation of WANBB is pretty simple, in that no a priori
physical knowledge is required and only two parameters are present without
tuning.

Atomic structure relaxation determines the ground-state atomic configuration
by searching for the local minimum in the energy landscape. It makes the
foundation and potential performance bottleneck for the search of global
structure minimum [1, 2, 3] and high-throughput calculations in material
design [4, 5, 6, 7, 8]. Conjugate gradient methods (CG) [9], the direct
inversion in the iterative subspace (DIIS) [10], and limited-memory Broyden-
Fletcher-Goldfarb-Shanno quasi-Newton methods (LBFGS) [11] are three widely
used methods to relax the atomic structure following the density functional
theory (DFT) calculation. While DIIS and LBFGS could be slightly faster in
some cases, CG is more stable in the situations where the initial atomic
configuration is far from equilibria. For large systems, hundreds of steps may
be required to converge to the typical atomic force tolerance of 0.01 eV/Å. So
it will be tremendously helpful even if the ab initio atomic relaxation could
be sped up to some degree.

One reason for the slow convergence of the local minimization goes to the
narrowly curved energy valley, which prevents the efficient execution of the
conventional CG. Thus, various preconditioners have been designed for
improving search directions via approximating the Hessian matrices. Zhao et
al. [12] notice that the Hessian matrix of large quantum dots can be
approximated by that of motifs comprising merely a few atoms. Some other works
construct general force fields based on a priori physical knowledge. Surrogate
potentials are used to either directly approximate the Hessian matrix [13, 14,
15], possibly with on-the-flight fittings [16, 17], or generate the
“preconditioned force” with the implicit inversion of the Hessian matrix [18].

Another reason is the inefficient execution of the line minimization (LM) to
which not much attention has been paid in the previous study. Actually,
gradient descent method and CG are two typical LM-based methods [9, 19, 20]. A
flowchart of the general LM-based method for the ab initio atomic relaxation
is depicted in Fig. 1.

StartptInput initial atomic positions $R_{0}$,force tolerance $\varepsilon>0$;
let $k:=0$ptForce tolerance satisfied?StopUpdate search direction
$D_{k}$ptCompute initial trial stepsize$\alpha_{k}^{\text{trial}}$;let scaling
factor $r_{k}:=1$ptUse an LM algorithm to update $r_{k}$;compute new trial
$R_{k}^{\text{trial}}:=R_{k}+r_{k}\alpha_{k}^{\text{trial}}D_{k}$ptLM
criterion accepts trial?Update atomic positions $R_{k+1}$; let
$k:=k+1$KSYesNoYesNoKSLM processExisting worksThis work Figure 1: A flowchart
of the general LM-based method, where $N\in\mathbb{N}$ refers to the number of
atoms, $R_{k}\in\mathbb{R}^{3\times N}$ stands for the Cartesian atomic
coordinates, $F_{k}\in\mathbb{R}^{3\times N}$ denotes the atomic forces
applied at $R_{k}$, and $D_{k}\in\mathbb{R}^{3\times N}$ is the search
direction. It should be noted that the white bold text “KS” refers to solving
the corresponding Kohn-Sham equations to obtain new energies and forces
whenever the atomic positions are updated. The LM process is marked out by the
shaded box.

The ab initio LM process, marked out by the shaded box in Fig. 1, is carried
out along search directions such as steepest descent or CG. During this
process, one begins with an initial trial stepsize and then repeatedly invokes
an LM algorithm to generate new trials until certain LM criterion is met.
Involving moving atoms, each trial step is followed by ab initio calculation,
e.g., solving the Kohn-Sham (KS) equation. But only the last accepted trial
step triggers the update of search direction. Taking the conventional CG for
instance, we test nearly 80 systems and find that each update of search
direction is followed by 2.4 calls of LM algorithms on average. In other
words, nearly 60% of the total computational overhead is paid for the
unaccepted trial steps in the LM process. Moreover, preconditioning with the
(fitted) classical force fields does not necessarily waive the cost in the LM
process [17]. As a result, the computational cost on the LM process is far
more than that on simply determining a search direction if the whole procedure
gets frequently trapped in the inner loop.

In this work, we present an LM-based gradient descent method, called WANBB,
which avoids trapping in the LM process during atomic structure relaxation. We
achieve this by (i) calculating an initial trial stepsize that captures the
local curvature of potential energy surface (PES) and (ii) devising an
unrestrictive LM criterion that tends to accept early trials. Nearly 80
systems have been tested, including organic molecules, metallic systems,
semiconductors, surface-molecule adsorption systems, ABX3 perovskites, etc.
The set of benchmark tests favors the universality of everyday ab initio
atomic relaxation. The unaccepted trial steps in the LM process account for
only about 1% of the total solved KS equations on average in WANBB, compared
with nearly 60% in CG. The considerable compression of the LM process leads to
a prominent saving on running time. The average speedup factors of WANBB over
CG, DIIS, and LBFGS are about 1.5, 1.2, and 1.2, respectively. The robustness
of the atomic relaxation method, i.e., convergence to equilibria regardless of
initial configurations, is also investigated. While CG, DIIS, and LBFGS fail
on some systems, WANBB manages to converge across the benchmark. This
robustness is theoretically established and helpful for the cases with local
lattice distortion or reconstruction. Last but not least, WANBB is pretty
simple. Free of a priori physical knowledge, it can work with only two
parameters present but no tuning required.

In what follows, we are ready to deliver the algorithmic development. To
facilitate narration, we collect some notations beforehand. We denote by
$E(R)\in\mathbb{R}$ and $F(R)\in\mathbb{R}^{3\times N}$, respectively, the
potential energy and atomic forces evaluated at $R$. When describing
algorithms, we use subscripts for abbreviation; e.g., both $E_{k}$ and $F_{k}$
are evaluated at $R_{k}$. We denote by $\mathbf{I}$ the identity matrix in
$\mathbb{R}^{3\times N}$. The operator $\left\langle\cdot,\cdot\right\rangle$
represents the Frobenius inner product of two $3$-by-$N$ matrices

$\left\langle
M,N\right\rangle:=\mathrm{Tr}(M^{\top}N),\quad\forall~{}M,N\in\mathbb{R}^{3\times
N},$

whereas $\|\cdot\|_{\mathrm{F}}$ yields the Frobenius norm as
$\|M\|_{\mathrm{F}}:=\sqrt{\left\langle M,M\right\rangle}$ for any
$M\in\mathbb{R}^{3\times N}$. We further denote by $\|M\|_{2,\infty}$ the
$\ell_{2,\infty}$-norm of matrix $M$, namely, the maximum among the
$\ell_{2}$-norms of columns in $M$. In particular, $\|F\|_{2,\infty}$ gives
the maximum $\ell_{2}$-norm among the atomic forces. Throughout iterations, we
take the atomic forces $F_{k}$ as the search direction $D_{k}$.

The whole procedure shown in Fig. 1 can fall into the trap of the LM process
owing to two factors. One is the poorly chosen initial trial stepsize, which
degrades the starting states of the LM process. The other is the restrictive
LM criterion, which acts as a stringent stopping rule for the LM algorithm.
The ideal situation is that the procedure always passes the test of the LM
criterion with just the original initial stepsize. Intuitively, we require
both good initial stepsizes and a loose but safe criterion to achieve this
target. Here, “safe” means that the convergence property is maintained with
the LM criterion.

At first sight, a good initial stepsize should grasp, as much as possible, the
curvature of PES around the current configuration. As early as 1988, Barzilai
and Borwein [21] propose the following renowned Barzilai-Borwein (BB)
stepsizes: for any $k\geq 1$,

$\alpha^{\text{BB1}}_{k}:=\frac{\left\langle
S_{k-1},S_{k-1}\right\rangle}{\left\langle
S_{k-1},Y_{k-1}\right\rangle},\quad\alpha^{\text{BB2}}_{k}:=\frac{\left\langle
S_{k-1},Y_{k-1}\right\rangle}{\left\langle Y_{k-1},Y_{k-1}\right\rangle},$

where $S_{k-1}:=R_{k}-R_{k-1}$ and $Y_{k-1}:=F_{k-1}-F_{k}$ represent the
difference of successive atomic positions and forces, respectively.
Remarkably, these two stepsizes satisfy

$\displaystyle\alpha^{\text{BB1}}_{k}=$
$\displaystyle\operatorname*{arg\,min}_{\alpha}\|\alpha^{-1}S_{k-1}-Y_{k-1}\|_{\mathrm{F}},$
$\displaystyle\alpha^{\text{BB2}}_{k}=$
$\displaystyle\operatorname*{arg\,min}_{\alpha}\|\alpha
Y_{k-1}-S_{k-1}\|_{\mathrm{F}}.$

In other words, $\left(\alpha^{\text{BB1}}_{k}\right)^{-1}\mathbf{I}$ and
$\left(\alpha^{\text{BB2}}_{k}\right)^{-1}\mathbf{I}$ approximate the Hessian
matrix around $R_{k}$ along the last search direction. They capture the local
curvature of PES with merely two $3$-by-$N$ matrices. The alternating BB (ABB)
stepsize [22] is presented as follows

$\alpha^{\text{ABB}}_{k}:=\left\\{\begin{array}[]{ll}\alpha^{\text{BB1}}_{k},&\text{if}~{}\mathrm{mod}(k,2)=1,\\\
\alpha^{\text{BB2}}_{k},&\text{if}~{}\mathrm{mod}(k,2)=0,\end{array}\right.$
(1)

which is shown to overwhelm BB1 and BB2 in many contexts (e.g., [22, 23]). It
is not difficult to prove that the implementation of (A)BB stepsizes requires
the same computational and storage complexity as computing the conjugate
parameter in CG.

Despite the simplicity, the (A)BB stepsizes turn out to be significant
benefits in solving various optimization problems. The most related
application goes to solving KS equations as orthogonality constrained
optimization problems, where the (A)BB stepsizes play as important local
acceleration components; see, e.g., [23, 24, 25, 26, 27]. However, the
application of the (A)BB stepsizes on the atomic structure relaxation problem
remains to be investigated. In this work, we take the (A)BB stepsizes as the
initial trial stepsizes to exploit their merits.

Without LM, the energy values produced by gradient descent methods with the
(A)BB stepsizes do not necessarily decrease monotonically. This phenomenon
reminds us of the possible divergence when starting the algorithms far from
equilibria. Nevertheless, the monotone LM criterion, adopted by the
conventional CG, can ruin the merits of the (A)BB stepsizes, since frequent
calls of LM algorithms are entailed to assure the monotonic decrease of the
energy and the local curvature of PES may be lost. Therefore, it is more
sensible to design a nonmonotone LM (NLM) criterion, accepting the original
(A)BB stepsizes in most cases.

The existing NLM criteria primarily fall into two categories: the “max” type
(MNLM) [28] and the “average” type (ANLM) [29]. In the $(k+1)$-th iteration,
the MNLM criterion asks the LM algorithm to return a factor $r_{k}>0$ such
that

$E\left(R_{k}+r_{k}\alpha_{k}^{\text{trial}}F_{k}\right)\leq\bar{E}_{k}-c\cdot
r_{k}\alpha_{k}^{\text{trial}}\|F_{k}\|_{\mathrm{F}}^{2},$

where $\bar{E}_{k}:=\max_{i=k-m(k)+1}^{k}\left\\{E_{i}\right\\}$,
$m(k)\in[0,M]$, both $c$ and $M$ are given positive constants, i.e., there is
a sufficient reduction with respect to (w.r.t.) the largest one of the past
$m(k)$ energies. The ANLM criterion requires $r_{k}$ to fulfill

$\displaystyle E\left(R_{k}+r_{k}\alpha_{k}^{\text{trial}}F_{k}\right)$
$\displaystyle\leq C_{k}-c\cdot
r_{k}\alpha_{k}^{\text{trial}}\|F_{k}\|_{\mathrm{F}}^{2},$ (2)
$\displaystyle\left\langle
F\left(R_{k}+r_{k}\alpha_{k}^{\text{trial}}F_{k}\right),F_{k}\right\rangle$
$\displaystyle\leq\sigma\|F_{k}\|_{\mathrm{F}}^{2},$ (3)

where $\sigma\in(0,1)$, the monitoring sequence $\left\\{C_{k}\right\\}$ is
updated through

$C_{k+1}:=\frac{\eta_{k}Q_{k}C_{k}+E_{k+1}}{\eta_{k}Q_{k}+1},~{}Q_{k+1}:=\eta_{k}Q_{k}+1$
(4)

with $C_{0}:=E_{0}$, $Q_{0}:=1$, and $\eta_{k}\in[0,1]$. Roughly in the ANLM
criterion, a sufficient reduction w.r.t. $C_{k}$ is obligatory. It is
straightforward to deduce that $C_{k}$ is a weighted average of the past
energies $\\{E_{j}\\}_{j=0}^{k-1}$ and the current $E_{k}$. The construction
of the monitoring sequence is completely free of a priori physical knowledge.
Clearly, entailing reduction w.r.t. $\bar{E}_{k}$ and $C_{k}$ respectively
instead of $E_{k}$, both MNLM and ANLM criteria are less restrictive than the
monotone criterion adopted by the conventional CG because $\bar{E}_{k}$,
$C_{k}\geq E_{k}$. Moreover, as noted in [29], the MNLM criterion is sensitive
to the choice of $M$ in some contexts, and is often outperformed by the ANLM
criterion.

Figure 2: The energy curves when relaxing the silicene surface with one
adsorbed H2O molecule using gradient descent method equipped with the ABB
stepsizes and ANLM ($\eta_{k}\equiv 0.85$222This value is recommended by
[29].) or WANLM ($\mu_{k}\equiv 0.05$) criterion. “#KS” refers to the number
of solved KS equations. We plot out the history with $\text{\\#KS}\leq 30$.
The blue solid line with square markers and orange line with triangle markers
stand for energies produced by “ABB + ANLM” and “ABB + WANLM”, respectively.
The blue dashed line with square markers and orange line with triangle markers
represent $\\{C_{k}\\}$ in ANLM and $\\{B_{k}\\}$ in WANLM, respectively. The
black arrows indicate the ab initio calculations dedicated to the trials
rejected by the ANLM criterion.

However, preliminary simulations show that the ANLM criterion is not loose
enough for the (A)BB stepsizes during the atomic relaxation. We take the
gradient descent direction armed with the ABB stepsizes and ANLM criterion to
relax the silicene surface with one adsorbed H2O molecule. From the short
relaxation history in Fig. 2, one can observe that the monitoring sequence
$\\{C_{k}\\}$ (the blue dashed line with square markers) decreases fast, stays
close to the energies $\\{E_{k}\\}$, and does not tolerate somewhat large
increments. As a result, the combination “ABB + ANLM” solves 3 KS equations
for the unaccepted trials during the short history, which are marked out by
the black arrows in Fig. 2. For example, the initial trial step at the 8-th
iteration gives an energy value of -99.92 eV, which is higher than $C_{7}=$
-100.03 eV and thus is rejected by the ANLM criterion.

To this end, we devise a reweighted ANLM (WANLM) criterion, which differs from
the ANLM criterion in the weights assigned to the past and current energies.
Specifically, we denote the monitoring sequence in WANLM by
$\left\\{B_{k}\right\\}$ and update it via

$B_{k+1}:=\frac{B_{k}+\mu_{k}P_{k}E_{k+1}}{1+\mu_{k}P_{k}},~{}P_{k+1}:=1+\mu_{k}P_{k}$
(5)

with $B_{0}:=E_{0}$, $P_{0}:=1$, and $\mu_{k}\in[0,1]$. Then WANLM needs
$r_{k}$ to satisfy both Eqs. (2) and (3) with $C_{k}$, $Q_{k}$, and $\eta_{k}$
replaced by $B_{k}$, $P_{k}$, and $\mu_{k}$, respectively.

The looseness of the WANLM criterion can be understood via the comparison with
the ANLM criterion. Due to Eq. (4) and $\eta_{0}\in[0,1]$, the coefficient of
$E_{1}$ is $1/(1+\eta_{0})\geq 0.5$, indicating a large portion of $E_{1}$ in
$C_{1}$ and a small gap $C_{1}-E_{1}$. Since the monitoring sequence
$\\{C_{k}\\}$ decreases monotonically (see the appendix for a rigorous proof),
it would then stay close to the energy values $\\{E_{k}\\}$, rendering the
ANLM criterion stringent for the (A)BB stepsizes as reflected in Fig. 2. In
contrast, from Eq. (5) and $\mu_{0}\in[0,1]$, the coefficient of $E_{1}$ is
$\mu_{0}/(1+\mu_{0})\leq 0.5$. Therefore, the monitoring sequence
$\\{B_{k}\\}$ in the WANLM criterion is relatively far away from the energies
$\\{E_{k}\\}$ at the beginning and gradually approaches the latter one,
leaving sufficient room to accept the (A)BB stepsizes. As shown in Fig. 2, the
first 7 energy values obtained by “ABB + ANLM” and “ABB + WANLM” are the same.
But one can find a much larger distance from $\\{E_{k}\\}$ to $\\{B_{k}\\}$
than $\\{C_{k}\\}$, which enables the WANLM criterion to accept the 8-th
initial trial step without invoking extra ab initio calculations and retain
the merits of the ABB stepsizes. All the subsequent iterations of “ABB +
WANLM” benefit from this characteristic and progress to lower energy values
faster. Regarding safety, it can be verified that the sufficient reduction in
terms of the atomic force norm holds for the monitoring sequence
$\\{B_{k}\\}$. One can then establish the convergence to equilibria,
regardless of initial configurations, after noticing the lower boundedness of
$\\{B_{k}\\}$. For more details, please refer to the appendix.

The gradient descent method combining the (A)BB stepsizes with the WANLM
criterion is referred to as WANBB. We set the initial trial stepsize to be
$4.8\times 10^{-2}$ Å2/eV at the first iteration, and take the alternating
strategy described in Eq. (1) afterwards because it leads to better
performance than BB1 and BB2 in our simulations. Before invoking an LM
algorithm, we take absolute value in case that $\left\langle
S_{k},Y_{k}\right\rangle<0$ and make truncation for safeguard:

$\alpha_{k}^{\text{ABB}}:=\min\left\\{\left|\alpha_{k}^{\text{ABB}}\right|,\max\left\\{-\log_{10}(\|F_{k}\|_{2,\infty}),1.0\right\\}\right\\}.$
(6)

As is shown in the appendix, the convergence of WANBB to equilibria can be
attributed to both Eqs. (2) and (3). For the sake of implementation, it is
more efficient to fullfill only Eq. (2). In WANLM, we set $\mu_{k}\equiv
0.05$, $c=10^{-4}$. The LM algorithm used by WANBB forms quadratic or cubic
interpolations based on recent information of PES; for details see [19,
Chapter 3].

We have implemented WANBB in the in-house plane-wave code CESSP [30, 31, 32,
33]. The exchange-correlation energy is described by the generalized gradient
approximation [34]. Electron-ion interactions are treated with the projector
augmented-wave (PAW) potentials based on the open-source ABINIT Jollet-
Torrent-Holzwarth dataset library in the PAW-XML format [35]. The total
energies are calculated using the Monkhorst-Pack mesh [36] with the $k$-mesh
spacing of 0.2 Å-1. The typical plane-wave energy cutoffs are around 500~600
eV. The KS equations corresponding to the updated atomic positions are solved
by the preconditioned self-consistent field (SCF) iteration [32].

In addition to WANBB, we test the performances of three relaxation methods,
including the conventional CG [19, 9], LBFGS [11], and DIIS [10]. By the way,
DIIS is the so-called quasi-Newton method in some software such as VASP [37,
38], which can be derived from Broyden’s method [39]. The implementations of
CG and DIIS follow those in [19, 11], and LBFGS is realized by the transition
state tool VTST [40]. CG employs Brent’s method [41] as the LM algorithm.

The benchmark for performance test contains nearly 80 systems from various
categories, including organic molecules, metallic systems, semiconductors,
surface-molecule adsorption systems, and ABX3 perovskites. The number of atoms
ranges from 2 to 429. Some of them are available online, e.g., in Materials
Project [42] and Organic Materials Database [43]. In the beginning, these
systems are at their ideal crystal positions for defects, heterostructures,
and substitutional alloys, or simply place a molecule on top of a surface.
These are the typical initial configurations used by researchers in their
simulations. More information on the benchmark can be found in Supplementary
Material 333See Supplemental Material at [URL will be inserted by publisher]
for information on the benchmark set and detailed numerical results.. The
atomic structure relaxation is terminated if either the maximum atomic force
falls below 0.01 eV/Å or the number of the solved KS equations (#KS) arrives
at 1000. The convergence criterion for the SCF iteration is $10^{-5}$ eV.

We evaluate the performances of the relaxation methods in two ways. One is the
number of the solved KS equations (#KS) needed to reach the 0.01 eV/Å force
tolerance. The other is the running time (CPU) of the atomic structure
relaxation. We adopt Dolan and Moré’s performance profile [45] for an overall
comparison among CG, WANBB, DIIS, and LBFGS. To this end, for any method
indicated by $m\in\left\\{\text{CG,~{}WANBB,~{}DIIS,~{}LBFGS}\right\\}$, we
denote its #KS and CPU for relaxing the $n$-th system by
$\text{\\#KS}_{n,m},~{}\text{CPU}_{n,m}$, respectively 444For visualization,
we reset $\text{\\#KS}_{n,m}=2000$ and $\text{CPU}_{n,m}=864000$ when the
method $m$ fails on the $n$-th system., define ratios

$r_{n,m}^{\text{KS}}:=\frac{\text{\\#KS}_{n,m}}{\min_{m}\text{\\#KS}_{n,m}},\quad
r_{n,m}^{\text{CPU}}:=\frac{\text{CPU}_{n,m}}{\min_{m}\text{CPU}_{n,m}},$

and for any $\omega\geq 0$, let

$\displaystyle\pi_{m}^{\text{KS}}(\omega)$ $\displaystyle:=\frac{\text{the
number of systems satisfying }r_{n,m}^{\text{KS}}\leq\omega}{\text{the number
of systems}},$ $\displaystyle\pi_{m}^{\text{CPU}}(\omega)$
$\displaystyle:=\frac{\text{the number of systems satisfying
}r_{n,m}^{\text{CPU}}\leq\omega}{\text{the number of systems}}.$

The quantities $\pi_{m}^{\text{KS}}(\omega)$ and
$\pi_{m}^{\text{CPU}}(\omega)$ estimate respectively the probabilities of
$r_{n,m}^{\text{KS}}\leq\omega$ and $r_{n,m}^{\text{CPU}}\leq\omega$ over the
benchmark 555To give a fair and expressive numerical comparison, we (i) filter
out the the systems where the differences between the maximum and the minimum
of the converged energies per atom by the four solvers are larger than 1 meV;
(ii) retain the systems at which certain solvers fail. The final number of
systems used for Fig. 3 is 68.. For example,
$\pi_{\text{WANBB}}^{\text{CPU}}(1)$ yields the portion of the systems where
WANBB is the fastest. And $\pi_{\text{WANBB}}^{\text{CPU}}(2)$ means the
estimated probability that, for one system, the running time of WANBB is not
more than twice the shortest running time needed by CG, WANBB, DIIS, and
LBFGS. In terms of the #KS, $\pi_{m}^{\text{KS}}(\omega)$ can be interpreted
in a similar way. We plot $\pi_{m}^{\text{KS}}(\omega)$ and
$\pi_{m}^{\text{CPU}}(\omega)$ as a function of $\omega$ for each method in
Fig. 3. Basically, the larger the area under the curve, the better the overall
performance of the method. It is easy to see that the overall performance of
WANBB ranks the best. More specifically, one can tell from the supplementary
Tables I and II [44] that the average CPU speedup factors of WANBB over CG,
DIIS, and LBFGS are 1.51, 1.21 and 1.16, respectively.

Figure 3: The #KS (left) and CPU (right) performance profiles of CG, WANBB,
DIIS, and LBFGS.

As shown in Fig. 3, though sometimes requiring solving more KS equations than
LBFGS, WANBB turns out to be faster than LBFGS in most tested systems. LBFGS
may make larger atomic displacements due to constant stepsize, which results
in underestimating the quality of initial wave functions and slowing down the
SCF iteration. Later, we will explain it in more detail using the relaxation
of the PuO2 surface with one adsorbed H2O molecule (see Fig. 9).

In addition to the efficiency, the robustness is also worth the whistle. As
shown in the supplementary Tables I and II [44], WANBB converges in all the
benchmark systems, while CG fails in 1 case due to the breakdown of Brent’s
method, DIIS and LBFGS fail in 8 cases and 1 case, respectively, due to
divergence. This kind of robustness is theoretically guaranteed by the
convergence property of WANBB proved in the appendix.

Taking the performance of CG as the baseline, we illustrate the acceleration
obtained by WANBB for each system, as depicted in Fig. 4. We have an average
CPU speedup factor around 1.51, and a factor above 2.0 on about 20% of
systems. According to the supplementary Tables I and II [44], CG takes on
average 117.68 calls of solving the KS equations to reach the force tolerance.
Thus a speedup factor of 1.5~2 can mean a significant saving of computational
overhead. This saving is mainly attributed to the enhanced ab initio LM
process of WANBB since the distribution of #KS speedup is consistent with that
of CPU speedup in Fig. 4. Moreover, we also count the percentage of the LM
steps invoked by WANBB and CG for unaccepted trials and draw the frequency
distribution over the benchmark in Fig. 5. The average LM step percentages for
unaccepted trials are 59.09% and 1.47% for CG and WANBB, respectively. It is
obvious from these statistics that WANBB escapes from the LM process quickly,
usually without invoking LM algorithms at all.

Figure 4: The speedup factors of #KS and the #KS of CG (the upper two), the
speedup factors of CPU and the CPU of CG (the lower two). The 1.0 baselines
are marked out by bold black dashed lines. The system where CG fails is not
included. Figure 5: The frequency distributions of the LM step percentages of
WANBB and CG for unaccepted trials.

The following discussion is concentrated on the performance comparison between
CG and WANBB for each category. In Fig. 6, we plot the mean speedup factors by
category and the number of atoms, respectively. All description below is
supported by Tables I and II in Supplementary Material [44]. Whenever speaking
of “speedup” below, we are talking about the CPU speedup if not specified.

Figure 6: Mean factors of speedup by category (left) and the number of atoms
(right).

Organic molecules. This class of systems includes various alkanes and amino
acids. The amino acids contain common organic atomic types, e.g., oxygen,
nitrogen, sulfur, etc. The number of atoms in this class ranges from 26 to
106. Most of the initial configurations are available in Organic Materials
Database [43]. Compared with CG, WANBB achieves a mean speedup factor of 1.81
for this class and speedup factors greater than 2 on over 20% of systems. In
the case of Glutamate, the acceleration factor even reaches 3.18, which means
a saving of nearly 70% computational overhead.

Metallic systems. We consider Ag16 cluster, Cu(111) surface with adsorbed
organic molecules, and some alloys including Ni3Al based intermetallic
compound as well as high-entropy alloys (HEAs) such as FeCrNiCoAl. The speedup
factors on these systems w.r.t. CG are mostly around 1.5~3. In particular, for
Cu(111) surface with adsorbed organic molecules, CG needs to solve over 400 KS
equations for atomic relaxations, while WANBB cuts the required #KS roughly in
half. It is also worth mentioning that WANBB can efficiently treat HEAs with
potential strong local lattice distortion [48]. Actually, local lattice
distortion is an essential factor in studying HEAs [49], and can be described
by the supercell method with ab initio atomic structure relaxation. In Fig. 7,
the local lattice distortion obtained by relaxing the body-centered cubic
supercell of FeCrNiCoAl is illustrated by the smearing radial distribution
function (RDF), in sharp contrast to the peaks at coordination shells in the
case of the ideal lattice. On average in our benchmark, WANBB achieves a
speedup of 1.54 over CG for HEAs.

Figure 7: RDF of FeCrNiCoAl with 6 Å cutoff.

Semiconductors & heterojunction systems. We test semiconductors (CdSe, GaAs,
Si, etc.) and heterojunction systems (GaAs-InAs). Moreover, we consider
semiconductors with defects by constructing the bulk GaAs and Si supercells
with some vacancies. In most of these systems, the acceleration factors are
above 1.5, and sometimes exceed 3. We also investigate the performances of
WANBB compared with CG on Si supercells with increasing system sizes; see Fig.
8. As the system size increases, WANBB is more stable in terms of #KS, and its
CPU advantage becomes more prominent.

Figure 8: Left: the average #KS and CPU speedup factors of WANBB over CG on
$N\times 1\times 1$ Si supercells containing from 8 to 256 atoms, where each
size has 10 random samples and the shift on each atom away from equilibrium is
no larger than 0.1 Å. Right: the #KS of WANBB and CG on $N\times 1\times 1$ Si
supercells, where the solid lines represent the average values, the upper and
lower boundaries of the shaded area stand for the maximum and minimum #KS,
respectively.

ABX3 perovskite systems. This type of systems is fairly popular in recent
solar cell studies. We consider the cases where “A” = CH3NH3, CH(NH2)2, or Cs,
“B” = Pb or Sn, “X” = Br, Cl, or I. Their initial configurations come from
Materials Project [42]. Due to the existence of the organic molecules and
their relative weak interaction with the inorganic framework, the atomic
relaxations of such systems are quite challenging. Taking CH3NH3SnI3 for
instance, CG solves 581 KS equations, while WANBB entails only about one fifth
of the expenditure.

Surface-molecule adsorption systems. Several surface-molecule adsorption
systems have been tested, including a 2D silicene surface with one adsorbed
H2O molecule (total 21 atoms), a 2D MoS2 surface with one adsorbed NH3
molecule (total 16 atoms), 2D PbS(001) and (111) surfaces with one adsorbed
oleic acid molecule ((001) total 254 atoms, (111) total 429 atoms), and a 2D
PuO2 surface with one adsorbed H2O molecule (total 195 atoms). The results
show that CG often undergoes quite a long way until convergence in this
category, and even fails in one case due to the breakdown of Brent’s method.
On the contrary, WANBB converges in all cases and usually costs less time; it
is even twice as fast as CG on PbS(001) surface with one adsorbed oleic acid
molecule. One interesting observation is that WANBB occasionally takes equal
or a few more #KS than CG, even though it is less time-costly. This is also
observed when comparing WANBB and LBFGS. To explain, we plot the cumulative
SCF iteration numbers of CG, WANBB, and LBFGS when relaxing the PuO2 surface
with one adsorbed H2O molecule; see Fig. 9.

Figure 9: Cumulative SCF iteration numbers of CG, WANBB, and LBFGS on PuO2
surface with one adsorbed H2O molecule.

One could observe that CG and LBFGS often consume more SCF iterations per KS
equation than WANBB. In fact, CG and LBFGS make larger atomic displacements,
resulting in degrading initial wave functions and slowing down the SCF
iteration. On the other hand, Fig. 9 also reveals that CG and LBFGS produce
more efficient search directions, yielding fewer iteration numbers. This
weakness, however, does not prevent WANBB from being a more favorable
alternative in the context of ab initio simulations, because the cut down on
the LM process obviously brings much more benefits.

2D and other systems. The 2D materials in our benchmark include C96, Si32, and
Germanene on hexagonal boron-nitride. The average speedup factors of WANBB
over CG on these systems are 1.26. We also include systems of other type,
e.g., cerium-iron oxide and potassium chloride. The former one is of practical
usage for NO${}_{\text{x}}$ reduction. WANBB reduces the numbers of solved KS
equations by half on both of them.

We recognize that the inefficiency of the LM process is an important obstacle
to the high-performance execution of ab initio atomic structure relaxation.
With the most widely used CG, nearly 60% of the computational overhead on
average is attributed to the unaccepted trials in the LM process. To this end,
we propose a force-based gradient descent method called WANBB, where the
initial trial (A)BB stepsizes grasp well the local curvature of PES and the
unrestrictive WANLM criterion often accepts the initial trials without calling
extra ab initio calculations. The robustness of WANBB, i.e., the convergence
to equilibria regardless of initial configurations, is theoretically
guaranteed. The numerical simulations on the benchmark containing nearly 80
various systems reveal the average 1.47% computational cost for unaccepted
trials in WANBB and demonstrate its universal speedups as well as superior
robustness over CG, DIIS, and LBFGS.

Our work concentrates on the LM process. Nevertheless, efficient search
directions together with subtly designed LM can be another great boon. This
topic merits futural study.

Acknowledgements

The authors acknowledge Jun Fang for the help with CESSP, and thank Zhen Yang
and Lifang Wang for the help on the construction of the benchmark. The
computations were (partly) done on the high performance computers of State Key
Laboratory of Scientific and Engineering Computing, Chinese Academy of
Sciences. This work was partially supported by the National Natural Science
Foundation of China under Grant Nos. 12125108, 11971466, 11991021, 11991020,
12021001, and 12288201, Key Research Program of Frontier Sciences, Chinese
Academy of Sciences (No. ZDBS-LY-7022), the CAS AMSS-PolyU Joint Laboratory in
Applied Mathematics, and the Science Challenge Project under Grant No.
TZ2018002.

*

##  Appendix A Convergence Analysis of WANBB

In this part, we analyze the convergence property of WANBB. Before moving on,
we define the level set

$\mathcal{L}:=\left\\{R:E(R)\leq E_{0}\right\\}.$

For the ease of reference, we restate the nonmonotone conditions in WANLM as
follows: in the $(k+1)$-th iteration, $R_{k+1}$ satisfies

$\displaystyle E_{k+1}$ $\displaystyle\leq
B_{k}-c\cdot\alpha_{k}\|F_{k}\|_{\mathrm{F}}^{2},$ (7)
$\displaystyle\left\langle F_{k+1},F_{k}\right\rangle$
$\displaystyle\leq\sigma\|F_{k}\|_{\mathrm{F}}^{2}.$ (8)

By the arguments in [29, Lemma 1.1], we know that the stepsize $\alpha_{k}$
satisfying both Eqs. (7) and (8) always exists and that $B_{k}$ is a convex
combination of $\left\\{E_{i}\right\\}_{i=0}^{k}$.

We first show that the atomic position sequence generated by WANBB remains in
$\mathcal{L}$.

######  Lemma 1.

Let $\left\\{R_{k}\right\\}$ be the atomic position sequence generated by
WANBB with $\mu_{k}\in[\mu_{\min},\mu_{\max}]\subseteq[0,1]$. Then we have
$\left\\{R_{k}\right\\}\subseteq\mathcal{L}$.

######  Proof.

We prove by contradiction. Suppose there exists $l_{1}>0$ such that
$E_{l_{1}}>E_{0}$. By the condition (7),

$E_{l_{1}}\leq
B_{l_{1}-1}-c\cdot\alpha_{l_{1}-1}\|F_{l_{1}-1}\|_{\mathrm{F}}^{2}\leq
B_{l_{1}-1}.$

Since $B_{l_{1}-1}$ is a convex combination of
$\left\\{E_{0},\ldots,E_{l_{1}-1}\right\\}$, there must exist some $l_{2}\leq
l_{1}-1$ for which $E_{l_{1}}\leq E_{l_{2}}$ holds. We can repeat the above
arguments to obtain a nonnegative sequence $\left\\{l_{i}\right\\}$ which is
finite for its strict monotonicity and ends up with 0. Therefore, we get
$E_{0}<E_{l_{1}}\leq E_{l_{2}}\leq\cdots\leq E_{0}$, a contradiction. ∎

From the update formula (5), one could obtain that the surrogate sequence
$\\{B_{k}\\}$ is always lower bounded by the potential energy sequence
$\\{E_{k}\\}$.

######  Lemma 2.

Let $\\{R_{k}\\}$ be the atomic position sequence generated by WANBB with
$\mu_{k}\in[\mu_{\min},\mu_{\max}]\subseteq[0,1]$. Then we have $E_{k}\leq
B_{k}$ for any $k\geq 0$.

######  Proof.

The claim is true for $k=0$ by definition. For $k\geq 1$, by the update
formula (5) of $B_{k+1}$, one has

$B_{k}=\frac{B_{k-1}+\mu_{k-1}Q_{k-1}E_{k}}{1+\mu_{k-1}Q_{k-1}}\geq\frac{E_{k}+\mu_{k-1}Q_{k-1}E_{k}}{1+\mu_{k-1}Q_{k-1}}=E_{k},$

where the inequality follows from the condition (7). The proof is complete. ∎

With Lemmas 1 and 2 in place, we are ready to show the convergence of WANBB.
The proof relies on the conditions (7) and (8) as well as the update formula
(5). Specifically, the condition (8) ensures a uniform positive lower bound of
stepsizes, while Eqs. (7) and (5) guide us to a sufficient reduction over
$\\{B_{k}\\}$.

######  Theorem 1.

Suppose the potential energy is bounded from below and $F$ is Lipschitz
continuous, with modulus $L>0$, on $\mathcal{L}$, namely,

$\|F(R)-F(\bar{R})\|_{\mathrm{F}}\leq
L\|R-\bar{R}\|_{\mathrm{F}},~{}\forall~{}R,\bar{R}\in\mathcal{L}.$

Let $\left\\{R_{k}\right\\}$ be the atomic position sequence generated by
WANBB with $\mu_{k}\in[\mu_{\min},\mu_{\max}]\subseteq(0,1]$. Then we have

$\lim_{k\to\infty}\|F_{k}\|_{\mathrm{F}}=0.$

That is to say, the atoms will approach equilibria, regardless of initial
configurations.

######  Proof.

Substracting $\|F_{k}\|_{\mathrm{F}}^{2}$ from both side of Eq. (8), we have

$\left\langle
F_{k}-F_{k+1},F_{k}\right\rangle\geq\left(1-\sigma\right)\|F_{k}\|_{\mathrm{F}}^{2}.$

Combining Lemma 1, the Lipschitz continuity of $F$ in $\mathcal{L}$, and the
fact $\sigma\in(0,1)$, we can further derive

$\displaystyle(1-\sigma)\|F_{k}\|_{\mathrm{F}}^{2}$
$\displaystyle\leq\|F_{k}-F_{k+1}\|_{\mathrm{F}}\cdot\|F_{k}\|_{\mathrm{F}}$
$\displaystyle\leq L\|R_{k}-R_{k+1}\|\cdot\|F_{k}\|$
$\displaystyle=\alpha_{k}L\|F_{k}\|_{\mathrm{F}}^{2},$

and hence a lower bound for stepsize whenever $F_{k}$ does not vanish:

$\alpha_{k}\geq\frac{1-\sigma}{L}.$ (9)

Plugging Eq. (9) into condition (7), we obtain

$E_{k+1}\leq B_{k}-c\cdot\alpha_{k}\|F_{k}\|_{\mathrm{F}}^{2}\leq
B_{k}-\frac{c(1-\sigma)}{L}\|F_{k}\|_{\mathrm{F}}^{2}.$ (10)

Let $\beta:=c(1-\sigma)/L$. By the update formula (5) of $B_{k+1}$, we obtain
a recursion for the surrogate sequence using Eq. (10):

$\displaystyle B_{k+1}$
$\displaystyle=\frac{B_{k}+\mu_{k}Q_{k}E_{k+1}}{1+\mu_{k}Q_{k}}$ (11)
$\displaystyle\leq\frac{B_{k}+\mu_{k}Q_{k}\left(B_{k}-\beta\|F_{k}\|_{\mathrm{F}}^{2}\right)}{Q_{k+1}}$
$\displaystyle=B_{k}-\beta\frac{\mu_{k}Q_{k}}{Q_{k+1}}\|F_{k}\|_{\mathrm{F}}^{2}.$

A byproduct of Eq. (11) is the monotonicity of $\left\\{B_{k}\right\\}$. Since
$E$ is bounded from below on $\mathcal{L}$ and $B_{k}\geq E_{k}$ for any
$k\geq 0$ (by Lemma 2), $\left\\{B_{k}\right\\}$ is lower bounded as well.
Consequently, by Eq. (11),

$\sum_{k=0}^{\infty}\frac{\mu_{k}Q_{k}\|F_{k}\|_{\mathrm{F}}^{2}}{Q_{k+1}}\leq\frac{1}{\beta}\sum_{k=0}^{\infty}\left(B_{k}-B_{k+1}\right)<+\infty.$
(12)

Following from the update formula of $\left\\{Q_{k}\right\\}$ and the fact
$\mu_{\max}\geq\mu_{k}\geq\mu_{\min}>0$, $Q_{k}\geq 1$, we can derive

$\frac{\mu_{k}Q_{k}}{Q_{k+1}}\geq\frac{\mu_{\min}Q_{k}}{\mu_{k}Q_{k}+1}\geq\frac{\mu_{\min}}{\mu_{k}+1/Q_{k}}\geq\frac{\mu_{\min}}{\mu_{\max}+1},$

which, together with Eq. (12), yields the summability of $\|F_{k}\|$. Hence,
$\|F_{k}\|\to 0$ as $k\to\infty$ as desired. ∎

## References

  * Oganov and Glass [2006] A. R. Oganov and C. W. Glass, J. Chem. Phys. 124, 244704 (2006).
  * Wang  _et al._ [2010] Y. Wang, J. Lv, L. Zhu, and Y. Ma, Phys. Rev. B 82, 094116 (2010).
  * Chen  _et al._ [2017a] Z. Chen, W. Jia, X. Jiang, S. Li, and L.-W. Wang, Comput. Phys. Commun. 219, 35 (2017a).
  * Goedecker  _et al._ [2005] S. Goedecker, W. Hellmann, and T. Lenosky, Phys. Rev. Lett. 95, 055501 (2005).
  * Vilhelmsen and Hammer [2012] L. B. Vilhelmsen and B. Hammer, Phys. Rev. Lett. 108, 126101 (2012).
  * Hu  _et al._ [2013] C. H. Hu, A. R. Oganov, Q. Zhu, G. R. Qian, G. Frapper, A. O. Lyakhov, and H. Y. Zhou, Phys. Rev. Lett. 110, 165504 (2013).
  * Li  _et al._ [2013] Q. Li, D. Zhou, W. Zheng, Y. Ma, and C. Chen, Phys. Rev. Lett. 110, 136403 (2013).
  * Zhang  _et al._ [2013] L. Zhang, J. Luo, A. Saraiva, B. Koiller, and A. Zunger, Nat. Commun. 4, 1 (2013).
  * Shewchuk [1994] J. Shewchuk,  _An Introduction to the Conjugate Gradient Method without the Agonizing Pain_, Tech. Rep. (Carnegie Mellon University, Pittsburgh, 1994).
  * Pulay [1980] P. Pulay, Chem. Phys. Lett. 73, 393 (1980).
  * Liu and Nocedal [1989] D. Liu and J. Nocedal, Math. Program. 45, 503 (1989).
  * Zhao  _et al._ [2006] Z. Zhao, L.-W. Wang, and J. Meza, Phys. Rev. B 73, 193309 (2006).
  * Goedecker  _et al._ [2001] S. Goedecker, F. Lançon, and T. Deutsch, Phys. Rev. B 64, 161102(R) (2001).
  * Fernández-Serra  _et al._ [2003] M. V. Fernández-Serra, E. Artacho, and J. M. Soler, Phys. Rev. B 67, 100101(R) (2003).
  * Packwood  _et al._ [2016] D. Packwood, J. Kermode, L. Mones, N. Bernstein, J. Woolley, N. Gould, C. Ortner, and G. Csányi, J. Chem. Phys. 144, 164109 (2016).
  * Chen  _et al._ [2014] Z. Chen, J. Li, S. Li, and L.-W. Wang, Phys. Rev. B 89, 144110 (2014).
  * Chen  _et al._ [2017b] Z. Chen, J. Li, S. Li, and L.-W. Wang, Phys. Rev. B 96, 115141 (2017b).
  * Liu  _et al._ [2018] L. Liu, Z. Chen, Y. Yao, and L.-W. Wang, Phys. Rev. B 98, 144109 (2018).
  * Nocedal and Wright [2006] J. Nocedal and S. Wright,  _Numerical Optimization 2nd Edition_ (Springer, New York, NY, 2006).
  * Press  _et al._ [2007] W. Press, H. William, S. Teukolsky, W. Vetterling, A. Saul, and B. Flannery,  _Numerical Recipes 3rd Edition: the Art of Scientific Computing_ (Cambridge University Press, 2007).
  * Barzilai and Borwein [1988] J. Barzilai and J. Borwein, IMA J. Numer. Anal. 8, 141 (1988).
  * Dai and Fletcher [2005] Y. Dai and R. Fletcher, Numer. Math. 100, 21 (2005).
  * Gao  _et al._ [2019] B. Gao, X. Liu, and Y. Yuan, SIAM J. Sci. Comput. 41, A1949 (2019).
  * Gao  _et al._ [2020] B. Gao, G. Hu, Y. Kuang, and X. Liu, An orthogonalization-free parallelizable framework for all-electron calculations in density functional theory, eprint arXiv:physics.comp-ph/2007.14228 (2020).
  * [25] B. Gao, X. Liu, L. Wang, N. Xiao, and Y. Yuan, STOP – a toolbox for Stiefel manifold optimization.
  * Xiao  _et al._ [2020] N. Xiao, X. Liu, and Y. Yuan, Optim. Methods Softw. 10.1080/10556788.2020.1852236 (2020).
  * Xiao  _et al._ [2021] N. Xiao, X. Liu, and Y. Yuan, A penalty-free infeasible approach for a class of nonsmooth optimization problems over the Stiefel manifold, eprint arXiv:math.OC/2103.03514 (2021).
  * Grippo  _et al._ [1986] L. Grippo, F. Lampariello, and S. Lucidi, SIAM J. Numer. Anal. 23, 707 (1986).
  * Zhang and Hager [2004] H. Zhang and W. Hager, SIAM J. Optim. 14, 1043 (2004).
  * Fang  _et al._ [2016] J. Fang, X. Gao, H. Song, and H. Wang, J. Chem. Phys. 144, 244103 (2016).
  * Gao  _et al._ [2017] X. Gao, Z. Mo, J. Fang, H. Song, and H. Wang, Comput. Phys. Commun. 211, 54 (2017).
  * Zhou  _et al._ [2018] Y. Zhou, H. Wang, Y. Liu, X. Gao, and H. Song, Phys. Rev. E 97, 033305 (2018).
  * Fang  _et al._ [2019] J. Fang, X. Gao, and H. Song, Commun. Comput. Phys. 26, 1196 (2019).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Jollet  _et al._ [2014] F. Jollet, M. Torrent, and N. Holzwarth, Comput. Phys. Commun. 185, 1246 (2014).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
  * Johnson [1988] D. D. Johnson, Phys. Rev. B 38, 12807 (1988).
  * [40] Jónsson group and Henkelman group, Transition state tools for VASP.
  * Brent [1973] R. P. Brent,  _Algorithms for Minimization without Derivatives_ (Prentice-Hall, Englewood Cliffs, New Jersey, 1973).
  * Jain  _et al._ [2013] A. Jain, S. Ong, G. Hautier, W. Chen, W. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. Persson, APL Mater. 1, 011002 (2013).
  * Borysov  _et al._ [2017] S. Borysov, R. Geilhufe, and A. Balatsky, PLoS ONE 12, e0171501 (2017).
  * Note [1] See Supplemental Material at [URL will be inserted by publisher] for information on the benchmark set and detailed numerical results.
  * Dolan and Moré [2002] E. Dolan and J. Moré, Math. Program. 91, 201 (2002).
  * Note [2] For visualization, we reset $\text{\\#KS}_{n,m}=2000$ and $\text{CPU}_{n,m}=864000$ when the method $m$ fails on the $n$-th system.
  * Note [3] To give a fair and expressive numerical comparison, we (i) filter out the the systems where the differences between the maximum and the minimum of the converged energies per atom by the four solvers are larger than 1 meV; (ii) retain the systems at which certain solvers fail. The final number of systems used for Fig. 3 is 68.
  * Wu  _et al._ [2020] J. Wu, Z. Yang, J. Xian, X. Gao, D. Lin, and H. Song, Front. Mater. , 396 (2020).
  * Song  _et al._ [2017] H. Song, F. Tian, Q. M. Hu, L. Vitos, Y. Wang, J. Shen, and N. Chen, Phys. Rev. Mater. 1, 023404 (2017).