Add WFN-history backend for Gamma-only NAO calculations

CMakeLists.txt

@@ -227,8 +227,11 @@ set(ABACUS_BIN_PATH ${CMAKE_CURRENT_BINARY_DIR}/${ABACUS_BIN_NAME})
 include_directories(${ABACUS_SOURCE_DIR})
 include_directories(${ABACUS_SOURCE_DIR}/source_base/module_container)
 
+set(ABACUS_MIN_CXX_STANDARD 14)
 if(NOT DEFINED CMAKE_CXX_STANDARD)
-  set(CMAKE_CXX_STANDARD 11)
+  set(CMAKE_CXX_STANDARD ${ABACUS_MIN_CXX_STANDARD})
+elseif(CMAKE_CXX_STANDARD LESS ABACUS_MIN_CXX_STANDARD)
+  message(FATAL_ERROR "ABACUS requires C++${ABACUS_MIN_CXX_STANDARD} or newer.")
 endif()
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
@@ -331,7 +334,6 @@ if(ENABLE_LCAO)
     target_link_libraries(${ABACUS_BIN_NAME} ${PEXSI_LIBRARY} ${SuperLU_DIST_LIBRARY} ${ParMETIS_LIBRARY} ${METIS_LIBRARY} pexsi)
     include_directories(${PEXSI_INCLUDE_DIR} ${ParMETIS_INCLUDE_DIR})
     add_compile_definitions(__PEXSI)
-    set(CMAKE_CXX_STANDARD 14)
   endif()
 else()
   set(ENABLE_MLALGO OFF)
@@ -452,7 +454,6 @@ if(USE_CUDA)
   # Always find CUDAToolkit to get CUDAToolkit_VERSION
   find_package(CUDAToolkit REQUIRED)
 
-  set_if_higher(CMAKE_CXX_STANDARD 14)
   if(CUDAToolkit_VERSION VERSION_GREATER_EQUAL "13.0")
     message(STATUS "CUDA ${CUDAToolkit_VERSION} detected. Setting CMAKE_CUDA_STANDARD to 17.")
     set_if_higher(CMAKE_CXX_STANDARD 17)
@@ -687,8 +688,6 @@ if(ENABLE_MLALGO OR DEFINED Torch_DIR)
   find_package(Torch REQUIRED)
   if(NOT Torch_VERSION VERSION_LESS "2.1.0")
     set_if_higher(CMAKE_CXX_STANDARD 17)
-  elseif(NOT Torch_VERSION VERSION_LESS "1.5.0")
-    set_if_higher(CMAKE_CXX_STANDARD 14)
   endif()
   include_directories(${TORCH_INCLUDE_DIRS})
   if(MKL_FOUND)
@@ -750,7 +749,6 @@ if(DEFINED LIBRI_DIR)
   set(ENABLE_LIBRI ON)
 endif()
 if(ENABLE_LIBRI)
-  set_if_higher(CMAKE_CXX_STANDARD 14)
   if(LIBRI_DIR)
   else()
     find_package(LibRI REQUIRED)
@@ -870,6 +868,8 @@ if(ENABLE_LCAO)
   target_link_libraries(
     ${ABACUS_BIN_NAME}
     hamilt_lcao
+    operator_ks_lcao
+    lcao_extrap
     tddft
     orb
     gint

source/source_esolver/esolver_fp.cpp

@@ -155,8 +155,20 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
     if (ucell.ionic_position_updated)
     {
         this->CE.update_all_dis(ucell);
-        this->CE.extrapolate_charge(&this->Pgrid, ucell, &this->chr, &this->sf,
-                                    GlobalV::ofs_running, GlobalV::ofs_warning);
+
+        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)
+        {
+            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);
+        }
     }
 
     //! calculate D2 or D3 vdW

source/source_esolver/esolver_ks_lcao.cpp

@@ -1,4 +1,5 @@
 #include "esolver_ks_lcao.h"
+#include "source_base/global_function.h"
 #include "source_base/module_external/blacs_connector.h"
 #include "source_cell/module_neighbor/sltk_atom_arrange.h"
 #include "source_estate/elecstate_tools.h"
@@ -17,8 +18,11 @@
 #include "../source_lcao/module_ri/exx_opt_orb.h"
 #endif
 #include "source_lcao/module_rdmft/rdmft.h"
+#include "source_lcao/module_extrap/wf_history_lcao.h"
+#include "source_lcao/module_operator_lcao/operator_lcao.h" // build overlap SR before WFN reorthonormalization
 #include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006
 #include "source_estate/module_dm/init_dm.h" // init dm from electronic wave functions
+#include "source_estate/module_dm/cal_dm_psi.h" // rebuild DMK from extrapolated WFN
 #include "source_io/module_ctrl/ctrl_runner_lcao.h" // use ctrl_runner_lcao() 
 #include "source_io/module_ctrl/ctrl_iter_lcao.h" // use ctrl_iter_lcao() 
 #include "source_io/module_ctrl/ctrl_scf_lcao.h" // use ctrl_scf_lcao()
@@ -27,6 +31,11 @@
 #include "source_lcao/LCAO_set.h" // mohan add 20251111
 #include "source_psi/setup_psi.h" // use Setup_Psi for deallocate_psi
 
+#include <iomanip>
+#include <sstream>
+#include <string>
+#include <type_traits>
+
 namespace ModuleESolver
 {
 
@@ -83,6 +92,14 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     LCAO_domain::set_psi_occ_dm_chg<TK>(this->kv, this->psi, this->pv, this->pelec,
       this->dmat, this->chr, inp);
 
+    const auto wfc_extrap_method = ModuleExtrap::wfc_extrap_method_from_string(inp.wfc_extrap);
+    if (wfc_extrap_method != ModuleExtrap::WfcExtrapMethod::None)
+    {
+        this->wf_history_lcao_ = std::make_unique<ModuleExtrap::WfHistoryLCAO<TK>>(wfc_extrap_method);
+        GlobalV::ofs_running << " WFN extrapolation method: "
+                             << ModuleExtrap::to_string(wfc_extrap_method) << std::endl;
+    }
+
     LCAO_domain::set_pot<TK>(ucell, this->kv, this->sf, *this->pw_rho, *this->pw_rhod,
       this->pelec, this->orb_, this->pv, this->locpp, this->dftu,
       this->solvent, this->exx_nao, this->deepks, inp);
@@ -176,6 +193,8 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
     this->dmat.dm->init_DMR(*hamilt_lcao->getHR());
 
+    bool initialized_by_wfc_extrap = false;
+
     // 13.1) decide the strategy for initializing DMR and HR
     if(istep == 0)//if the first scf step, readin DMR from file,
     {
@@ -194,12 +213,97 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
                 this->chr, PARAM.inp.ks_solver);
         }
     }
-    else if(PARAM.inp.esolver_type!="tddft")//if not, use the DMR calculated from last step
+    else if(PARAM.inp.esolver_type!="tddft")//initialize DMR from WFN history if required
     {
-        // 13.1.2) two cases are considered:
-        // 1. DMK in DensityMatrix is not empty (istep > 0), then DMR is initialized by DMK
-        // 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros
-        this->dmat.dm->cal_DMR();
+        if constexpr (std::is_same<TK, double>::value)
+        {
+            if (this->wf_history_lcao_ != nullptr)
+            {
+                // The newly-created HamiltLCAO has allocated S(R), but the actual
+                // overlap values are normally filled later by updateHk() inside the
+                // eigensolver.  WFN extrapolation needs S before entering SCF, so build
+                // only the overlap real-space matrix here and then fold it to S(Gamma).
+                auto* lcao_op = dynamic_cast<hamilt::OperatorLCAO<TK, TR>*>(hamilt_lcao->getOperator());
+                if (lcao_op == nullptr)
+                {
+                    ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_scf",
+                                             "Failed to access the LCAO operator chain for WFN extrapolation.");
+                }
+                ModuleBase::timer::start("WFN_Extrap", "prepare_overlap");
+                lcao_op->contributeHR();
+                // Refresh the current Gamma-only overlap matrix before reusing the previous WFN.
+                // The layout must follow the active LCAO solver because the WFN orthonormalizer
+                // assembles distributed local blocks using the same column/row-major convention.
+                const int sk_layout =
+                    ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver) ? 1 : 0;
+                hamilt_lcao->updateSk(0, sk_layout);
+                ModuleBase::timer::end("WFN_Extrap", "prepare_overlap");
+
+                ModuleBase::timer::start("WFN_Extrap", "apply");
+                const ModuleExtrap::WfExtrapApplyResult wfc_result
+                    = this->wf_history_lcao_->try_use_prev_wf_gamma(hamilt_lcao->getSk(),
+                                                                     this->pv,
+                                                                     *(this->psi),
+                                                                     this->pelec->wg);
+                ModuleBase::timer::end("WFN_Extrap", "apply");
+                if (wfc_result.ok())
+                {
+                    ModuleBase::timer::start("WFN_Extrap", "rebuild_density");
+                    elecstate::cal_dm_psi(this->dmat.dm->get_paraV_pointer(),
+                                          this->pelec->wg,
+                                          *(this->psi),
+                                          *(this->dmat.dm));
+                    this->dmat.dm->cal_DMR();
+                    LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
+                    ModuleBase::timer::end("WFN_Extrap", "rebuild_density");
+                    initialized_by_wfc_extrap = true;
+
+                    GlobalV::ofs_running << " WFN extrapolation: use_prev_wf from ionic step "
+                                         << wfc_result.snapshot_istep
+                                         << ", active bands = " << wfc_result.nactive_bands
+                                         << ", max |C^T S C - I| = "
+                                         << wfc_result.max_orthonormality_deviation << std::endl;
+                }
+                else
+                {
+                    std::ostringstream oss;
+                    oss << std::scientific << std::setprecision(16);
+                    oss << "WFN extrapolation failed with status: "
+                        << ModuleExtrap::to_string(wfc_result.status) << "\n\n"
+                        << " snapshot_istep=" << wfc_result.snapshot_istep << "\n"
+                        << " failed_state=" << wfc_result.failed_state << "\n"
+                        << " failed_pivot_index=" << wfc_result.failed_pivot_index << "\n"
+                        << " failed_pivot=" << wfc_result.failed_pivot << "\n"
+                        << " nstate=" << wfc_result.nstate << "\n"
+                        << " nbands=" << wfc_result.nbands << "\n"
+                        << " nbasis=" << wfc_result.nbasis << "\n"
+                        << " min_metric_diag=" << wfc_result.min_metric_diag << "\n"
+                        << " max_metric_diag=" << wfc_result.max_metric_diag << "\n"
+                        << " max_metric_abs=" << wfc_result.max_metric_abs << "\n"
+                        << " max_metric_asymmetry=" << wfc_result.max_metric_asymmetry << "\n"
+                        << " max_orthonormality_deviation="
+                        << wfc_result.max_orthonormality_deviation << "\n";
+                    const std::string message = oss.str();
+                    ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_scf", message);
+                }
+            }
+        }
+        else
+        {
+            if (this->wf_history_lcao_ != nullptr)
+            {
+                ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_scf",
+                                         "WFN extrapolation is currently supported only for the real Gamma-only NAO path.");
+            }
+        }
+
+        if (!initialized_by_wfc_extrap)
+        {
+            // 13.1.2) two cases are considered:
+            // 1. DMK in DensityMatrix is not empty (istep > 0), then DMR is initialized by DMK
+            // 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros
+            this->dmat.dm->cal_DMR();
+        }
     }
     // 13.2) init_scf, should be before_scf? mohan add 2025-03-10
     elecstate::init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric,
@@ -548,6 +652,11 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
             this->rdmft_solver, this->deepks, this->exx_nao,
             this->conv_esolver, this->scf_nmax_flag, istep);
 
+    if (conv_esolver && this->wf_history_lcao_ != nullptr && this->psi != nullptr)
+    {
+        this->wf_history_lcao_->update_after_scf(istep, *(this->psi), this->pelec->wg);
+    }
+
     //! 3) Clean up RA, which is used to serach for adjacent atoms
     if (!PARAM.inp.cal_force && !PARAM.inp.cal_stress)
     {

source/source_esolver/esolver_ks_lcao.h

@@ -21,6 +21,12 @@ template <typename T, typename TR>
 class ESolver_LR;
 }
 
+namespace ModuleExtrap
+{
+template <typename TK>
+class WfHistoryLCAO;
+}
+
 //-----------------------------------
 // ESolver for LCAO
 //-----------------------------------
@@ -81,6 +87,11 @@ class ESolver_KS_LCAO : public ESolver_KS
     //! Add density matrix class, mohan add 2025-10-30
     LCAO_domain::Setup_DM<TK> dmat;
 
+    //! Optional NAO wavefunction-history extrapolation state.
+    //! The object owns only copied wavefunction snapshots; transient objects such as
+    //! Psi, HamiltLCAO, DensityMatrix, and Charge are passed in only when needed.
+    std::unique_ptr<ModuleExtrap::WfHistoryLCAO<TK>> wf_history_lcao_;
+
 
     // For deepks method, mohan add 2025-10-08
     Setup_DeePKS<TK> deepks;

source/source_io/module_parameter/input_parameter.h

@@ -49,6 +49,7 @@ struct Input_para
     std::string init_chg = "atomic";    ///< "file","atomic"
     bool dm_to_rho = false;             ///< read density matrix from npz format and calculate charge density
     std::string chg_extrap = "default"; ///< xiaohui modify 2015-02-01
+    std::string wfc_extrap = "none";   ///< wavefunction extrapolation method for LCAO calculations
     bool init_vel = false;              ///< read velocity from STRU or not  liuyu 2021-07-14
 
     std::string input_file = "INPUT";   ///< input file name

source/source_io/module_parameter/read_input_item_system.cpp

@@ -859,6 +859,47 @@ Available options are:
         };
         this->add_item(item);
     }
+    {
+        Input_Item item("wfc_extrap");
+        item.annotation = "none; use_prev_wf";
+        item.category = "System variables";
+        item.type = "String";
+        item.description = R"(Wavefunction extrapolation method for LCAO calculations.
+
+* none: Disable wavefunction-based extrapolation.
+* use_prev_wf: Use the previous ionic step wavefunctions as the initial guess.
+
+This option is currently limited to Gamma-only LCAO calculations.
+The k-point, ASPC, and GExt_PROJ paths will be enabled by later updates.)";
+        item.default_value = "none";
+        read_sync_string(input.wfc_extrap);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.wfc_extrap != "none" && para.input.wfc_extrap != "use_prev_wf")
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "wfc_extrap should be none or use_prev_wf");
+            }
+            if (para.input.wfc_extrap == "use_prev_wf")
+            {
+                if (para.input.basis_type != "lcao")
+                {
+                    ModuleBase::WARNING_QUIT("ReadInput", "WFN-based extrapolation is available only for LCAO");
+                }
+                if (!para.input.gamma_only)
+                {
+                    ModuleBase::WARNING_QUIT("ReadInput", "WFN-based extrapolation is not implemented for k-points");
+                }
+                if (para.input.esolver_type == "tddft")
+                {
+                    ModuleBase::WARNING_QUIT("ReadInput", "WFN-based extrapolation is not implemented for tddft");
+                }
+                if (para.input.nspin == 4)
+                {
+                    ModuleBase::WARNING_QUIT("ReadInput", "WFN-based extrapolation is not implemented for nspin = 4");
+                }
+            }
+        };
+        this->add_item(item);
+    }
     {
         Input_Item item("ecutwfc");
         item.annotation = "energy cutoff for wave functions";

source/source_lcao/CMakeLists.txt

@@ -6,23 +6,9 @@ add_subdirectory(module_deltaspin)
 
 if(ENABLE_LCAO)
     add_subdirectory(module_operator_lcao)
-    list(APPEND objects
+    add_subdirectory(module_extrap)
+    set(hamilt_lcao_sources
         hamilt_lcao.cpp
-        module_operator_lcao/operator_lcao.cpp
-        module_operator_lcao/veff_lcao.cpp
-        module_operator_lcao/meta_lcao.cpp
-        module_operator_lcao/op_dftu_lcao.cpp
-        module_operator_lcao/deepks_lcao.cpp
-        module_operator_lcao/op_exx_lcao.cpp
-        module_operator_lcao/overlap.cpp
-        module_operator_lcao/ekinetic.cpp
-        module_operator_lcao/nonlocal.cpp
-        module_operator_lcao/td_ekinetic_lcao.cpp
-        module_operator_lcao/td_nonlocal_lcao.cpp
-        module_operator_lcao/td_pot_hybrid.cpp
-        module_operator_lcao/dspin_lcao.cpp
-        module_operator_lcao/dftu_lcao.cpp
-        module_operator_lcao/operator_force_stress_utils.cpp
         dftu_lcao.cpp
         pulay_fs_center2.cpp
         FORCE_STRESS.cpp
@@ -58,7 +44,7 @@ if(ENABLE_LCAO)
     add_library(
         hamilt_lcao
         OBJECT
-        ${objects}
+        ${hamilt_lcao_sources}
     )
 
     if(ENABLE_COVERAGE)

source/source_lcao/module_extrap/CMakeLists.txt

@@ -0,0 +1,14 @@
+add_library(
+    lcao_extrap
+    OBJECT
+    wf_history_lcao.cpp
+    wf_orthonormalize_lcao.cpp
+)
+
+if(BUILD_TESTING)
+  add_subdirectory(test)
+endif()
+
+if(ENABLE_COVERAGE)
+  add_coverage(lcao_extrap)
+endif()

source/source_lcao/module_extrap/test/CMakeLists.txt

@@ -0,0 +1,9 @@
+AddTest(
+  TARGET MODULE_LCAO_extrap_test
+  LIBS parameter ${math_libs} psi base device
+  SOURCES
+    wf_extrap_test.cpp
+    ../wf_history_lcao.cpp
+    ../wf_orthonormalize_lcao.cpp
+    ${ABACUS_SOURCE_DIR}/source_basis/module_ao/parallel_orbitals.cpp
+)