Parallel Orbital-Updating Method

• The parallel orbital-updating method is an efficient strategy that replaces a large coupled eigenvalue problem with many independent orbital source problems and a small projection eigensolve.
• It leverages two-level parallelism by enabling both orbital-wise and intra-orbital computations, thereby overcoming scalability bottlenecks in electronic structure calculations.
• Its variants—spanning real-space finite-element, plane-wave, and direct energy minimization formulations—improve convergence and reduce global orthogonalization overhead.

The parallel orbital-updating method is an orbital/eigenfunction iteration based approach for electronic structure calculations and, more generally, for eigenvalue problems in which many eigenpairs are required. Its defining reformulation is to replace repeated solutions of a large coupled eigenvalue problem by many independent source or correction problems, one per orbital, followed by a small projected eigenvalue problem in the span of the updated orbitals. In the Kohn–Sham setting, this strategy is motivated by the single-particle structure of the equations and by the observation that large-scale orthogonalization and dense subspace operations are often the dominant scalability bottlenecks in conventional eigensolvers (Dai et al., 2014, Pan et al., 2017, Dai et al., 2024).

1. Historical emergence and problem class

The method emerged in the context of Kohn–Sham density functional theory (KS-DFT), where one seeks the lowest $N$ orthonormal orbitals solving a nonlinear eigenvalue problem of the form

$$\left(-\frac12\Delta + V_{\mathrm{eff}}(\rho)\right)\psi_i = \varepsilon_i \psi_i, \qquad \int \psi_i \psi_j = \delta_{ij}, \qquad \rho = \sum_{i=1}^N |\psi_i|^2.$$

After discretization and linearization inside self-consistent field (SCF) iteration, standard approaches repeatedly solve a generalized algebraic eigenproblem $A u=\lambda B u$ or, in plane-wave notation, $H\Psi=\varepsilon S\Psi$. The associated cost is severe because many such eigenproblems are solved across SCF iterations, and because $N_g\gg N$ in realistic calculations. The real-space formulation emphasized that the optimal complexity of traditional approaches is $\mathcal O(N^2N_g)$ when $A$ and $B$ are sparse and $\mathcal O(NN_g^2)$ when $A$ or $$\left(-\frac12\Delta + V_{\mathrm{eff}}(\rho)\right)\psi_i = \varepsilon_i \psi_i, \qquad \int \psi_i \psi_j = \delta_{ij}, \qquad \rho = \sum_{i=1}^N |\psi_i|^2.$$0 is dense (Dai et al., 2014).

The earliest formulation in the supplied literature is the real-space, finite-element, full-potential approach of Dai, Gong, Zhou, and Zhu, which described the method as an orbital iteration based parallel approach for electronic structure calculations. A later paper transferred the same philosophy to reciprocal-space plane-wave DFT, explicitly presenting a parallel orbital-updating based plane-wave basis method and two modified variants, with implementation in Quantum ESPRESSO. Subsequent work extended the framework to direct energy minimization on the Stiefel/Grassmann manifold, and later mathematical papers established convergence and error estimates for linear eigenvalue problems, including clustered eigenvalues and adaptive finite element discretizations (Dai et al., 2014, Pan et al., 2017, Dai et al., 2015, Dai et al., 27 Aug 2025).

2. Core formulation and algorithmic mechanism

The central algorithmic idea is to freeze the effective operator at the current iterate, update each orbital independently, and postpone mutual coupling and orthogonality restoration to a low-dimensional projection step. In the KS formulation of the real-space paper, with current orbitals $$\left(-\frac12\Delta + V_{\mathrm{eff}}(\rho)\right)\psi_i = \varepsilon_i \psi_i, \qquad \int \psi_i \psi_j = \delta_{ij}, \qquad \rho = \sum_{i=1}^N |\psi_i|^2.$$1, each orbital is updated by solving

$$\left(-\frac12\Delta + V_{\mathrm{eff}}(\rho)\right)\psi_i = \varepsilon_i \psi_i, \qquad \int \psi_i \psi_j = \delta_{ij}, \qquad \rho = \sum_{i=1}^N |\psi_i|^2.$$2

and the updated subspace is

$$\left(-\frac12\Delta + V_{\mathrm{eff}}(\rho)\right)\psi_i = \varepsilon_i \psi_i, \qquad \int \psi_i \psi_j = \delta_{ij}, \qquad \rho = \sum_{i=1}^N |\psi_i|^2.$$3

A small eigenproblem is then solved in $$\left(-\frac12\Delta + V_{\mathrm{eff}}(\rho)\right)\psi_i = \varepsilon_i \psi_i, \qquad \int \psi_i \psi_j = \delta_{ij}, \qquad \rho = \sum_{i=1}^N |\psi_i|^2.$$4, restoring orthonormality and extracting the next approximate eigenspace. In operator language, the large eigenproblem is replaced by $$\left(-\frac12\Delta + V_{\mathrm{eff}}(\rho)\right)\psi_i = \varepsilon_i \psi_i, \qquad \int \psi_i \psi_j = \delta_{ij}, \qquad \rho = \sum_{i=1}^N |\psi_i|^2.$$5 independent linear boundary value or source problems plus one $$\left(-\frac12\Delta + V_{\mathrm{eff}}(\rho)\right)\psi_i = \varepsilon_i \psi_i, \qquad \int \psi_i \psi_j = \delta_{ij}, \qquad \rho = \sum_{i=1}^N |\psi_i|^2.$$6-scale Rayleigh–Ritz type problem (Dai et al., 2014).

The plane-wave formulation uses the same structure inside SCF. At iteration $$\left(-\frac12\Delta + V_{\mathrm{eff}}(\rho)\right)\psi_i = \varepsilon_i \psi_i, \qquad \int \psi_i \psi_j = \delta_{ij}, \qquad \rho = \sum_{i=1}^N |\psi_i|^2.$$7, for each orbital $$\left(-\frac12\Delta + V_{\mathrm{eff}}(\rho)\right)\psi_i = \varepsilon_i \psi_i, \qquad \int \psi_i \psi_j = \delta_{ij}, \qquad \rho = \sum_{i=1}^N |\psi_i|^2.$$8, the baseline algorithm solves in parallel

$$\left(-\frac12\Delta + V_{\mathrm{eff}}(\rho)\right)\psi_i = \varepsilon_i \psi_i, \qquad \int \psi_i \psi_j = \delta_{ij}, \qquad \rho = \sum_{i=1}^N |\psi_i|^2.$$9

builds

$A u=\lambda B u$0

and then solves a reduced eigenproblem in that space. The modified formulations replace the source equation by a residual correction equation

$A u=\lambda B u$1

or by a projected correction equation

$A u=\lambda B u$2

In both modified methods the reduced space is enriched to dimension $A u=\lambda B u$3 by including both old orbitals and correction vectors (Pan et al., 2017).

A closely related but distinct variant arises in direct energy minimization. There, the full manifold gradient

$A u=\lambda B u$4

is approximated by retaining only diagonal entries of $A u=\lambda B u$5, which leads to independent orbital-wise directions

$A u=\lambda B u$6

followed by updates

$A u=\lambda B u$7

and periodic reconstruction by

$A u=\lambda B u$8

This formulation is presented as a parallel orbital-updating based optimization method rather than as an SCF eigensolver (Dai et al., 2015).

3. Discretization settings and major algorithmic variants

The literature supplied here contains three principal incarnations of the method in KS electronic structure, each preserving the same decomposition into independent orbital updates plus low-dimensional correction.

Setting	Paper	Defining feature
Full-potential adaptive FEM	(Dai et al., 2014)	$A u=\lambda B u$9 independent source problems plus small eigensolve in $H\Psi=\varepsilon S\Psi$0
Plane-wave KS-DFT	(Pan et al., 2017)	Source or correction equations in plane-wave space, then reduced eigenproblem of dimension $H\Psi=\varepsilon S\Psi$1 or $H\Psi=\varepsilon S\Psi$2
Direct energy minimization	(Dai et al., 2015)	Single orbital-updating approximation with periodic orthogonalization and occasional subspace diagonalization

In the real-space finite element formulation, the method was developed for full-potential molecular systems and combined naturally with adaptive refinement. The paper also introduced a stabilized variant that initializes $H\Psi=\varepsilon S\Psi$3 orbitals and performs the update/projection process on $H\Psi=\varepsilon S\Psi$4 functions, remarking that this variant is more stable in practice. Suggested initial guesses included Gaussian-type orbitals or Slater-type orbitals for full-potential calculations, and local plane-wave or local finite element/volume guesses for pseudopotentials (Dai et al., 2014).

In the plane-wave formulation, the novelty was to adapt the real-space orbital-updating philosophy to reciprocal-space electronic structure, where conventional solvers are widely used but can become bottlenecked by FFT communication, large-scale orthogonalization, and synchronization-heavy dense linear algebra. The paper explicitly stated that Algorithms 2 and 3 inherit the main advantages of Algorithm 1 but increase the reduced eigenproblem dimension from $H\Psi=\varepsilon S\Psi$5 to $H\Psi=\varepsilon S\Psi$6. It also reported that Algorithm 3 gave results similar to Algorithm 2, so only Algorithm 2, denoted MParO, was reported numerically in detail (Pan et al., 2017).

In the optimization-based formulation, the method is expressed on the Stiefel manifold $H\Psi=\varepsilon S\Psi$7 and the Grassmann manifold $H\Psi=\varepsilon S\Psi$8, reflecting the rotational invariance of the occupied subspace. The practical algorithm used Barzilai–Borwein step sizes, a nonmonotone Armijo-type condition

$H\Psi=\varepsilon S\Psi$9

periodic orthogonalization, and periodic diagonalization of the projected Hamiltonian matrix $N_g\gg N$0 in order to keep the diagonal approximation effective (Dai et al., 2015).

4. Parallel structure, implementation, and computational behavior

A major theme across the literature is two-level parallelization. The first level is orbital-level parallelism: once the operator or Hamiltonian is frozen, the source or correction problems for different orbitals are independent and can be assigned to different processors or processor groups. The second level is intra-orbital parallelism: each source problem is itself a large linear system or PDE solve, so standard parallel mechanisms—multigrid, domain decomposition, parallel FFTs, or parallel matrix-vector operations—can be used within each orbital solve. The method therefore does not eliminate the underlying parallelization of the discretization; it adds a new layer above it (Pan et al., 2017, Dai et al., 2024).

In the plane-wave implementation in Quantum ESPRESSO, using norm-conserving pseudopotentials and modified Broyden density mixing, the comparison baseline was the code’s conjugate-gradient-like band-by-band diagonalization. The paper emphasized that the orbital-updating method changes the communication pattern favorably by confining orthogonality and eigenvalue coupling to a very small reduced space. Numerical experiments on MgO, Al, and Si supercells showed that the basic method ParO and especially MParO become advantageous as the system grows. For MgO supercells at $N_g\gg N$1-point and $N_g\gg N$2 Ry, the 512-atom case gave $N_g\gg N$3 s for CG, $N_g\gg N$4 s for ParO, and $N_g\gg N$5 s for MParO. For a $N_g\gg N$6 Si supercell with 1000 atoms and 2000 orbitals, MParO required $N_g\gg N$7 s on 80 processors and $N_g\gg N$8 s on 640 processors, whereas CG required $N_g\gg N$9 s and $\mathcal O(N^2N_g)$0 s, respectively. The Al tests also showed a reliability effect: for 256 atoms at $\mathcal O(N^2N_g)$1, CG failed to converge while ParO and MParO converged (Pan et al., 2017).

In the real-space direct minimization implementation in OCTOPUS 4.0.1, with LDA, LCAO initial guesses, and Troullier–Martins norm-conserving pseudopotential, the method was tested on systems ranging from benzene to protein-scale and carbon-cluster examples. Problem sizes reached $\mathcal O(N^2N_g)$2 orbitals, real-space grids up to $\mathcal O(N^2N_g)$3, and core counts up to 512. The proposed method, denoted Opt-Par, obtained essentially the same converged energies as the Opt-Z3W baseline, while exposing a large orbital-parallelizable fraction of the runtime. With $\mathcal O(N^2N_g)$4 and $\mathcal O(N^2N_g)$5, the orbital-parallelizable fraction was reported as 75% for $\mathcal O(N^2N_g)$6, 70% for 2JMO, 66% for FAS2, and 67% for $\mathcal O(N^2N_g)$7. With $\mathcal O(N^2N_g)$8 and $\mathcal O(N^2N_g)$9, the fraction increased to 82% for $A$0, 75% for 2JMO, and 73% for FAS2, while the total time for $A$1 decreased from $A$2 s for the baseline to $A$3 s for the orbital-updating method (Dai et al., 2015).

5. Mathematical analysis, clustered spectra, and adaptive extensions

The later numerical analysis papers place the method in an abstract Hilbert-space eigenproblem

$A$4

with symmetric bilinear forms and a discrete spectrum containing eigenvalue clusters with multiplicities $A$5. Their central analytical contribution is the identification of a quasi-orthogonality mechanism: even before explicit orthogonalization, independently updated orbitals are sufficiently close to mutually orthogonal exact eigendirections that one can recover an orthogonal basis with controlled error. This is the key step that allows rigorous analysis for clustered and multiple eigenvalues, where individual eigenvectors are not the correct primary object and distances between subspaces must be used instead (Dai et al., 2024).

For the projected ParO step, the analysis proves finite-step error estimates. If the updated component spaces $A$6 are close enough to the corresponding discrete eigenspaces $A$7, then after solving the small projected eigenproblem there exists a $A$8-orthonormal basis of each $A$9 such that

$B$0

and

$B$1

The same paper analyzes simplified and practical shifted-inverse based ParO algorithms. For the simplified method, the grouped eigenspace error contracts with an asymptotic factor

$B$2

analogous to classical shifted-inverse iteration. For the practical adaptive-shift method, the asymptotic factor becomes

$B$3

and in the exact discrete case, where $B$4, the paper derives cubic convergence (Dai et al., 2024).

The adaptive finite element analysis extends the method to an outer Solve $B$5 Estimate $B$6 Mark $B$7 Refine loop. On each mesh, ParO iterations are performed only until the algebraic error is of the same order as the current FE discretization error; the paper states explicitly that one should not oversolve the inner ParO iterations on coarse meshes. Residual-type estimators are constructed for both FE eigenspaces and ParO eigenspaces, and the analysis proves that the two are tightly related: $B$8 Under the stated assumptions, the adaptive ParO approximations satisfy

$B$9

with $\mathcal O(NN_g^2)$0 independent of mesh size. This provides a mathematical justification for adaptive ParO on clustered eigenvalue problems of linear elliptic operators (Dai et al., 27 Aug 2025).

6. Scope of the term and related but distinct usages

A recurring source of confusion is that parallel orbital-updating has been used for more than one computational task. In the KS-DFT and eigenvalue literature, the term refers to updating approximate orbitals or eigenfunctions separately by source or correction solves, then performing a reduced eigenspace correction. In contrast, the QMC paper “GPGPU for orbital function evaluation with a new updating scheme” used the phrase in a different sense: the central algorithmic innovation was quasi-simultaneous updating (Q.S.), defined by

$\mathcal O(NN_g^2)$1

so that all trial electron moves are tested independently against the same old configuration. This exposed $\mathcal O(NN_g^2)$2 independent orbital-evaluation tasks for CUDA execution and accelerated the orbital-evaluation bottleneck by about 30.67× for solid TiO$\mathcal O(NN_g^2)$3 with $\mathcal O(NN_g^2)$4 electrons on a GeForce GTX 480. That paper concerns ab-initio QMC orbital spline evaluation, not KS eigensolver reformulation (Uejima et al., 2012).

Another boundary case is the 2026 work on constrained optimization algorithms for orbital optimization in quantum chemistry. That paper is directly relevant to orbital updating in the broad sense because it formulates orbital optimization as a Stiefel-manifold constrained minimization driven by one- and two-particle reduced density matrices from MP2, CASCI, or DMRG. However, it explicitly does not propose a truly parallel orbital-updating algorithm in the sense of independently solvable orbital subproblems executed concurrently. Its relevance is architectural rather than direct: it separates the correlated solver from the orbital optimizer through an RDM interface, but it does not introduce orbital-wise parallel updates, block-parallel decomposition, or asynchronous orbital scheduling (Zhang et al., 16 Jun 2026).

Within the principal KS/eigenvalue lineage, the method is therefore best understood as a family of algorithms that replace the dominant large-scale eigensolve by many independent orbital updates plus a small reduced eigenproblem, thereby reducing global orthogonalization cost and exposing two-level parallelism. The real-space, plane-wave, optimization, and adaptive finite element variants differ in discretization, update equations, and convergence control, but they share the same structural principle and the same computational objective: better scalability and better robustness when many orbitals are required (Dai et al., 2014, Pan et al., 2017, Dai et al., 27 Aug 2025).