Add WFN-history backend for Gamma-only NAO calculations

[Growl1234]

This PR introduces an optional WFN-based initialization path for NAO calculations.

A new wfc_extrap use_prev_wf path is added. It stores the converged wavefunction from the previous ionic step, rebuilds the current overlap matrix, reorthonormalizes the previous wavefunction with the current overlap matrix, and then reconstructs the density matrix and charge density before entering the normal SCF cycle. When this path is enabled, the existing charge-density extrapolation is skipped, so the initial density is generated from the WFN-history backend instead of being extrapolated directly in real space.

The implementation currently focuses on the Gamma-only backend, including MPI-distributed matrices. This is intended as a common backend for future improvements and more WFN-based predictors such as ASPC and GExt_PROJ.

Here's a simple test that shows the accuracy and effectiveness of the WFN-based path: test.zip. The regtesting is left for a next commit.

Still missing:

• Try take a look if we can further reduce the complication of code.
• Documentation mention.

What's planned in the future:

• Support for k-points: the case would be more complicated because the backend has to handle complex per-k wavefunctions, k-dependent overlap matrices $S(k)$, local/global k-point ownership, density-matrix reconstruction from all k-points, and PBC phase or convention corrections between ionic steps.

Details

Linked Issue

None

Unit Tests and/or Case Tests for my changes

• Unittest is added in the new minor module
• Regtest: ‎02_NAO_Gamma/relax_wfc_extrap

Exact Verification Performed

• Commands run:
• Result summary:
• Checks not run, with reason:

What's changed?

• See above.

Governance Checklist

• Global dependencies: no net increase in GlobalV, GlobalC, or PARAM code references, or exception requested below with reason, scope, risk, and cleanup plan.
• Default parameters: no new default arguments added to existing interfaces, or exception requested below.
• Headers: no unnecessary header dependencies or .hpp propagation, or rationale provided below.
• Line endings: text files use LF; only .bat and .cmd use CRLF.
• Build linkage: new source files are listed in the relevant source/CMakeLists.txt, or rationale provided below.
• Documentation: behavior/interface changes include documentation updates [TODO]
• CodeRabbit: if automatic review has not started and the repository has CodeRabbit installed, request @coderabbitai review.

INPUT Parameter Changes

• Parameters added/removed/changed: wfc_extrap none/use_prev_wf
• docs/parameters.yaml updated: not yet
• docs/advanced/input_files/input-main.md updated: not yet
• If not updated, explain why no INPUT documentation update is required:

Core Module Impact

• Affected core modules: module_lcao
• Risk summary: No risks since this is a new minor module
• Compatibility or performance impact:

Governance Exception

• Rule:
• Reason:
• Scope:
• User or maintenance risk:
• Why the normal rule cannot be followed now:
• Follow-up cleanup plan:
• Requested approver:

[Copilot]

Pull request overview

Note

Copilot was unable to run its full agentic suite in this review.

Adds Gamma-only LCAO wavefunction-history extrapolation (use_prev_wf) by snapshotting the previous step’s coefficients and reorthonormalizing them against the current overlap matrix before SCF.

Changes:

• Introduces LCAO wavefunction snapshot/history helpers and a Gamma-only reorthonormalization routine.
• Wires extrapolation into the LCAO ESolver flow and adds wfc_extrap input validation/handling (including skipping charge extrapolation when enabled).
• Updates CMake to enforce a minimum C++ standard and adds a new module_extrap subdirectory target.

Reviewed changes

Copilot reviewed 14 out of 14 changed files in this pull request and generated 6 comments.

Show a summary per file

File	Description
source/source_lcao/module_extrap/wf_snapshot_lcao.h	Stores/restores full Psi coefficient snapshots plus occupations.
source/source_lcao/module_extrap/wf_orthonormalize_lcao.h	Declares Gamma-only reorthonormalization API + diagnostic result struct.
source/source_lcao/module_extrap/wf_orthonormalize_lcao.cpp	Implements MPI-aware assembly + Cholesky-based reorthonormalization.
source/source_lcao/module_extrap/wf_history_lcao.h	Defines wavefunction history class and apply-result diagnostics.
source/source_lcao/module_extrap/wf_history_lcao.cpp	Implements history management and use_prev_wf apply path.
source/source_lcao/module_extrap/wf_extrap_method.h	Adds extrapolation method/status enums and string conversions.
source/source_lcao/module_extrap/CMakeLists.txt	Adds lcao_extrap object library target for new module sources.
source/source_lcao/CMakeLists.txt	Enables module_extrap subdir and adds new sources to LCAO objects list.
source/source_io/module_parameter/read_input_item_system.cpp	Adds wfc_extrap input item, docs, and validation checks.
source/source_io/module_parameter/input_parameter.h	Adds wfc_extrap parameter with default "none".
source/source_esolver/esolver_ks_lcao.h	Adds wf_history_lcao_ member to ESolver.
source/source_esolver/esolver_ks_lcao.cpp	Initializes history, prepares overlap early, applies extrapolation, rebuilds density, updates history after SCF.
source/source_esolver/esolver_fp.cpp	Skips charge extrapolation when wfc_extrap is enabled (LCAO, istep>0).
CMakeLists.txt	Enforces a minimum C++ standard globally (now 14).

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[mohanchen]

Thank you for your pull request. This PR presents a creative idea, and we are open to merging it, yet several issues need to be discussed first:
First, has this algorithm been formally published? Are there other software packages that implement similar algorithms? Could you also document the relevant algorithm details within this issue thread?
Second, we need to confirm whether this algorithm delivers tangible performance benefits. Integrating it will require substantial modifications to the code framework and complicate the runtime logic of the LCAO module. We only intend to merge acceleration algorithms with proven practical value. If the feature ends up unused after merging, there is no justification for introducing extra complexity to the codebase.
Third, we prefer minimal edits to esolver files for all new feature implementations. For instance, the current implementation adds over one hundred lines of code to esolver_ks_lcao.cpp, which is not acceptable. Modifications limited to fewer than ten lines would be reasonable. The core design principle of the esolver module is to maintain a streamlined workflow; it is not meant to host detailed implementations of specific algorithms.

[Growl1234]

Thank you for reviewing!

the current implementation adds over one hundred lines of code to esolver_ks_lcao.cpp

Those code is most likely for debugging and finding errors. Once further validation from maintainer's side is done, they could be considered dropped.

[Growl1234]

First, has this algorithm been formally published? Are there other software packages that implement similar algorithms? Could you also document the relevant algorithm details within this issue thread?

This is not intended to introduce a "creative algorithm" or any "new method". It just implements the basic, lowest-order wavefunction-history initialization referred to as use_prev_wf here: the converged wavefunction from the previous ionic step is retained, reorthonormalized with respect to the current overlap matrix, and then used to reconstruct the density matrix and charge density before the normal SCF cycle.

This basic idea is already used, in various forms, in many quantum-chemistry and electronic-structure packages, including Gaussian, ORCA, and CP2K. In particular, the CP2K implementation was the direct motivation for introducing this capability in ABACUS. Its USE_PREV_WF path, together with methods built on top of the same infrastructure, especially ASPC, often performs well in practical geometry optimizations and molecular-dynamics calculations and is one of the useful strengths of CP2K in this area.

Moreover, this low-level capability is also a necessary foundation for introducing formally published higher-order wavefunction-extrapolation methods, such as ASPC. These methods do not merely "extrapolate" MO coefficient matrix; they require a consistent wavefunction-history backend, including the reorthonormalization of predicted orbitals in the current basis and the subsequent reconstruction of the density matrix and charge density.

[Growl1234]

Second, we need to confirm whether this algorithm delivers tangible performance benefits. Integrating it will require substantial modifications to the code framework and complicate the runtime logic of the LCAO module. We only intend to merge acceleration algorithms with proven practical value. If the feature ends up unused after merging, there is no justification for introducing extra complexity to the codebase.

For the current Gamma-only NAO benchmark, the WFN-based initialization reduces the total number of SCF iterations from 44 to 35 and reduces the total runtime from about 32.7 s to 28.2 s, while preserving the final result. This is the result from small-cell test, but if tested with large cells, the improvements should be more significant. That said, since it is not yet implemented with k-points, this PR should not be regarded as broadly applicable or as demonstrating general performance improvements across all ABACUS workflows.

[Growl1234]

Third, we prefer minimal edits to esolver files for all new feature implementations. For instance, the current implementation adds over one hundred lines of code to esolver_ks_lcao.cpp, which is not acceptable. Modifications limited to fewer than ten lines would be reasonable. The core design principle of the esolver module is to maintain a streamlined workflow; it is not meant to host detailed implementations of specific algorithms.

I agree with this point. The current implementation leaves too much WFN-specific logic in ESolver_KS_LCAO, while the esolver should only coordinate the general SCF workflow. I will attempt to move the parts regarding WFN-based extrapolation path into module_extrap later. This should reduce the feature-specific changes in the esolver files to a few calls and keep the detailed algorithmic implementation self-contained in module_extrap.

Btw, the diagnostics part in the implementation most likely can be dropped.

[Growl1234]

It seems moving specific logic from ESolver_KS_LCAO to module_extrap didn't reduce the complication of code but raised the risk of compilation failure? Unfortunately I won't have time to take a deeper look at this in the next weeks, so I'll shelve this PR for now.

[Copilot]

Pull request overview

Copilot reviewed 26 out of 26 changed files in this pull request and generated 1 comment.

Suppressed comments (4)

source/source_lcao/module_extrap/wf_density_init_lcao.cpp:115

• Rule: Do not increase cross-layer control through GlobalV, GlobalC, or PARAM
  Severity: warning
  Location: source/source_lcao/module_extrap/wf_density_init_lcao.cpp:113-115
  Reason: New GlobalV::ofs_running usage is introduced in the new extrapolation backend, increasing reliance on global state in a reusable module.
  Suggested action: Propagate success diagnostics (snapshot_istep, nactive_bands, max deviation) back to the ESolver layer and perform logging there, or accept and document a governance exception with reason/scope/risk/cleanup plan.
  Exception: human approval required

    GlobalV::ofs_running << " WFN extrapolation: use_prev_wf from ionic step " << result.snapshot_istep
                         << ", active bands = " << result.nactive_bands
                         << ", max |C^T S C - I| = " << result.max_orthonormality_deviation << std::endl;

source/source_lcao/module_extrap/wf_history_lcao.cpp:168

• reorthonormalize_gamma_lcao() uses the provided occupation matrix to decide how many bands to reorthonormalize (and leaves the rest untouched). Passing snapshot.wg (previous ionic step occupations) can be wrong when occupations change between steps (e.g., metallic systems / smearing), potentially leaving currently-occupied bands non-orthonormal and corrupting the DM rebuilt with wg_now.

    const WfOrthonormalizeResult orth_result = reorthonormalize_gamma_lcao(current_overlap,
                                                                           pv,
                                                                           psi_trial,
                                                                           snapshot.wg,
                                                                           1.0e-12,
                                                                           pivot_threshold,
                                                                           check_tolerance);

source/source_esolver/esolver_ks_lcao.cpp:3

• source_base/global_function.h appears unused in this translation unit (no ModuleBase::GlobalFunc usage). Keeping it increases compile dependencies and can introduce avoidable coupling to global headers.

#include "esolver_ks_lcao.h"
#include "source_base/global_function.h"
#include "source_base/module_external/blacs_connector.h"

source/source_esolver/esolver_fp.cpp:198

• Skipping charge-density extrapolation solely based on wfc_extrap/istep can leave chr/sf in a stale state when the WFN-history path is unavailable (e.g., previous SCF not converged so no snapshot was stored, or history was cleared). In that case initialize_gamma_density() falls back and returns false, but charge extrapolation has already been skipped, so the SCF may start from an inconsistent density for the updated ionic positions.

        const bool skip_charge_extrap_for_wfc = PARAM.inp.basis_type == "lcao"
                                                && PARAM.inp.wfc_extrap != "none"
                                                && istep > 0;
        if (skip_charge_extrap_for_wfc)
        {
            this->sf.setup(&ucell, this->Pgrid, this->pw_rhod);
            GlobalV::ofs_running << " charge density extrapolation is skipped because wfc_extrap = "
                                 << PARAM.inp.wfc_extrap << "." << std::endl;
        }
        else
        {
            this->CE.extrapolate_charge(&this->Pgrid, ucell, &this->chr, &this->sf,
                                        GlobalV::ofs_running, GlobalV::ofs_warning);
        }

[Copilot]

Pull request overview

Copilot reviewed 29 out of 29 changed files in this pull request and generated 1 comment.

Suppressed comments (1)

source/source_lcao/module_extrap/wf_orthonormalize_lcao.cpp:332

• Similar to overlap_global, coeff_global is materialized as a dense nactive_bands×nbasis_global array on every rank and combined via MPI_Allreduce. This can become a large O(nbands·nbasis) memory+bandwidth cost (in addition to the overlap gather) and may negate the runtime benefit of WFN-based initialization on bigger systems. Consider using a distributed representation/workflow (or adding a size guard and clear diagnostics) so wfc_extrap fails fast rather than risking memory exhaustion.

    TimerGuard timer("WFN_Extrap", "assemble_coeff");
    const int nbasis_global = pv.get_global_row_size();
    if (coeff_local == nullptr || nbasis_global <= 0 || ncoeff_local <= 0 || nbasis_local <= 0
        || nactive_bands < 0 || pv.get_row_size() != nbasis_local)
    {
        return WfcExtrapStatus::DimensionMismatch;
    }

    coeff_global.assign(static_cast<std::size_t>(nactive_bands) * static_cast<std::size_t>(nbasis_global), 0.0);

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.