Add WFN-history backend for Gamma-only NAO calculations

source/CMakeLists.txt

@@ -632,6 +632,8 @@ if(ENABLE_LCAO)
     ${ABACUS_BIN_NAME}
     PRIVATE
       hamilt_lcao
+      operator_ks_lcao
+      lcao_extrap
       tddft
       orb
       gint

source/source_esolver/esolver_fp.cpp

@@ -152,8 +152,21 @@ 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)
+        {
+            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);
+        }
     }
 
     //! 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"
@@ -27,6 +28,7 @@
 #include "source_lcao/LCAO_set.h" // mohan add 20251111
 #include "source_psi/setup_psi.h" // use Setup_Psi for deallocate_psi
 
+
 namespace ModuleESolver
 {
 
@@ -83,6 +85,8 @@ 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);
 
+    this->wf_history_lcao_.set_method(ModuleExtrap::wfc_extrap_method_from_string(inp.wfc_extrap));
+
     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);
@@ -194,12 +198,17 @@ 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 (!this->wf_history_lcao_.initialize_gamma_density(
+                *hamilt_lcao, this->pv, *(this->psi), this->pelec->wg, *(this->dmat.dm), this->chr,
+                PARAM.inp.nspin, PARAM.inp.ks_solver))
+        {
+            // 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,
@@ -550,6 +559,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->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

@@ -10,6 +10,7 @@
 #include "source_lcao/setup_exx.h" // for exx, mohan add 20251008
 #include "source_lcao/module_rdmft/rdmft.h" // rdmft
 #include "source_lcao/setup_dm.h" // mohan add 2025-10-30
+#include "source_lcao/module_extrap/wf_history_lcao.h"
 
 #include <memory>
 
@@ -81,6 +82,8 @@ class ESolver_KS_LCAO : public ESolver_KS
     //! Add density matrix class, mohan add 2025-10-30
     LCAO_domain::Setup_DM<TK> dmat;
 
+    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

@@ -50,6 +50,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

@@ -897,6 +897,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,15 @@
+add_library(
+    lcao_extrap
+    OBJECT
+    wf_history_lcao.cpp
+    wf_orthonormalize_lcao.cpp
+    wf_density_init_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
+)