Refactor the operators in source_lcao

source/Makefile.Objects

@@ -343,18 +343,28 @@ OBJS_HAMILT_OF=kedf_tf.o\
 OBJS_HAMILT_LCAO=hamilt_lcao.o\
     operator_lcao.o\
     ekinetic.o\
+    ekinetic_fs.o\
+    ekinetic_dh.o\
     nonlocal.o\
+    nonlocal_dh.o\
+    nonlocal_fs.o\
     overlap.o\
+    overlap_fs.o\
     td_ekinetic_lcao.o\
     td_nonlocal_lcao.o\
     td_pot_hybrid.o\
+    td_pot_hybrid_fs.o\
     veff_lcao.o\
+    veff_dh.o\
     meta_lcao.o\
     op_dftu_lcao.o\
     deepks_lcao.o\
     op_exx_lcao.o\
     dspin_lcao.o\
+    dspin_fs.o\
     dftu_lcao.o\
+    dftu_fs.o\
+    operator_fs_utils.o\
 
 OBJS_HCONTAINER=base_matrix.o\
     atom_pair.o\
@@ -721,7 +731,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     setup_pwrho.o\
     setup_pwwfc.o\
     update_cell_pw.o\
-    dftu_pw.o\
+    setup_dftu_pw.o\
     deltaspin_pw.o\
     forces.o\
     forces_us.o\
@@ -793,7 +803,7 @@ OBJS_DFTU=dftu.o\
       dftu_tools.o\
       dftu_occup.o\
       dftu_hamilt.o\
-      dftu_pw.o
+      setup_dftu_pw.o
 
 OBJS_DELTASPIN=basic_funcs.o\
       cal_mw_from_lambda.o\

source/source_esolver/esolver_ks_lcao.cpp

@@ -5,7 +5,7 @@
 #include "source_lcao/module_deltaspin/spin_constrain.h"
 #include "source_lcao/module_deltaspin/deltaspin_lcao.h"
 #include "source_lcao/dftu_lcao.h"
-#include "source_lcao/hs_matrix_k.hpp" // there may be multiple definitions if using hpp
+#include "source_hamilt/hs_matrix_k.h"
 #include "source_estate/module_charge/symmetry_rho.h"
 #include "source_lcao/LCAO_domain.h" // need DeePKS_init
 #include "source_lcao/FORCE_STRESS.h"

source/source_esolver/esolver_ks_pw.cpp

@@ -1,35 +1,27 @@
 #include "esolver_ks_pw.h"
 
-#include "source_cell/cal_ux.h"
 #include "source_estate/elecstate_pw.h"
 #include "source_estate/module_charge/symmetry_rho.h"
-#include "source_hamilt/module_xc/xc_functional.h" // use XC_Functional
 #include "source_hsolver/diago_iter_assist.h"
 #include "source_hsolver/diago_params.h"
 #include "source_hsolver/hsolver_pw.h"
-#include "source_hsolver/kernels/hegvd_op.h"
 #include "source_io/module_parameter/parameter.h"
-#include "source_lcao/module_deltaspin/spin_constrain.h"
-#include "source_lcao/module_dftu/dftu.h"
 #include "source_pw/module_pwdft/forces.h"
 #include "source_pw/module_pwdft/hamilt_pw.h"
-#include "source_pw/module_pwdft/onsite_proj.h"
 #include "source_pw/module_pwdft/stress_pw.h"
-#include "source_pw/module_pwdft/vsep_pw.h"
 
 #ifdef __DSP
 #include "source_base/kernels/dsp/dsp_connector.h"
 #endif
 
 #include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006
 #include "source_estate/setup_estate_pw.h"         // mohan add 20251005
-#include "source_estate/update_pot.h"              // mohan add 20251016
 #include "source_hamilt/module_xc/exx_info.h"      // use GlobalC::exx_info
 #include "source_io/module_ctrl/ctrl_output_pw.h"  // mohan add 20250927
 #include "source_pw/module_pwdft/deltaspin_pw.h"   // mohan add 20250309
-#include "source_pw/module_pwdft/dftu_pw.h"        // mohan add 20250309
 #include "source_pw/module_pwdft/setup_pot.h"      // mohan add 20250929
 #include "source_pw/module_pwdft/update_cell_pw.h" // mohan add 20250309
+#include "source_pw/module_pwdft/setup_dftu_pw.h"  // mohan add 20250309
 
 namespace ModuleESolver
 {
@@ -150,7 +142,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
 
     if (ucell.cell_parameter_updated)
     {
-        this->stp.p_psi_init->prepare_init(PARAM.inp.pw_seed);
+        this->stp.p_psi_init->prepare_init(PARAM.inp.pw_seed, istep);
     }
 
     //! Init Hamiltonian (cell changed)

source/source_hamilt/CMakeLists.txt

@@ -9,6 +9,7 @@ endif()
 
 list(APPEND objects
     operator.cpp
+    hs_matrix_k.cpp
     module_ewald/H_Ewald_pw.cpp
     module_ewald/dnrm2.cpp
 )

source/source_hamilt/hs_matrix_k.cpp

@@ -0,0 +1,7 @@
+#include "source_hamilt/hs_matrix_k.h"
+
+namespace hamilt
+{
+    template class HS_Matrix_K<double>;
+    template class HS_Matrix_K<std::complex<double>>;
+}

source/source_lcao/hs_matrix_k.hpp → source/source_hamilt/hs_matrix_k.h

@@ -1,8 +1,9 @@
-#ifndef HS_MATRIX_K_HPP
-#define HS_MATRIX_K_HPP
+#ifndef HS_MATRIX_K_H
+#define HS_MATRIX_K_H
 
 #include "source_basis/module_ao/parallel_orbitals.h"
 
+#include <complex>
 #include <vector>
 namespace hamilt
 {
@@ -42,4 +43,4 @@ namespace hamilt
     };
 }
 
-#endif
+#endif

source/source_io/module_dhs/write_dH_terms.cpp

@@ -6,7 +6,7 @@
 #include "source_hamilt/module_hcontainer/output_hcontainer.h"
 #include "source_lcao/module_operator_lcao/ekinetic.h"
 #include "source_lcao/module_operator_lcao/nonlocal.h"
-#include "source_lcao/module_operator_lcao/operator_force_stress_utils.h"
+#include "source_lcao/module_operator_lcao/operator_fs_utils.h"
 #include "source_lcao/module_operator_lcao/veff_lcao.h"
 #include "source_hamilt/module_gint/gint_interface.h"
 #include "source_lcao/module_lr/utils/lr_util_xc.hpp"

source/source_io/module_energy/write_eband_terms.hpp

@@ -29,8 +29,8 @@ void write_eband_terms(const int nspin,
                        const TwoCenterBundle& two_center_bundle
 #ifdef __EXX
                        ,
-                       std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd = nullptr,
-                       std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc = nullptr
+                       std::vector<std::map<int, std::map<hamilt::TAC, RI::Tensor<double>>>>* Hexxd = nullptr,
+                       std::vector<std::map<int, std::map<hamilt::TAC, RI::Tensor<std::complex<double>>>>>* Hexxc = nullptr
 #endif
 )
     {

source/source_io/module_hs/write_H_terms.cpp

@@ -14,7 +14,7 @@
 #include "source_hamilt/module_hcontainer/output_hcontainer.h"
 #include "source_lcao/module_operator_lcao/ekinetic.h"
 #include "source_lcao/module_operator_lcao/nonlocal.h"
-#include "source_lcao/module_operator_lcao/operator_force_stress_utils.h"
+#include "source_lcao/module_operator_lcao/operator_fs_utils.h"
 #ifdef __EXX
 #include "source_lcao/module_operator_lcao/op_exx_lcao.h"
 #include "source_lcao/module_ri/Exx_LRI_interface.h"

source/source_io/module_hs/write_vxc.hpp

@@ -156,8 +156,8 @@ void write_Vxc(const int nspin,
                bool cal_exx
 #ifdef __EXX
                ,
-               std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd = nullptr,
-               std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc = nullptr
+               std::vector<std::map<int, std::map<hamilt::TAC, RI::Tensor<double>>>>* Hexxd = nullptr,
+               std::vector<std::map<int, std::map<hamilt::TAC, RI::Tensor<std::complex<double>>>>>* Hexxc = nullptr
 #endif
 )
 {

source/source_io/module_hs/write_vxc_r.hpp

@@ -41,8 +41,8 @@ void write_Vxc_R(const int nspin,
                  bool real_number
 #ifdef __EXX
                  ,
-                 const std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* const Hexxd,
-                 const std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* const Hexxc
+                 const std::vector<std::map<int, std::map<hamilt::TAC, RI::Tensor<double>>>>* const Hexxd,
+                 const std::vector<std::map<int, std::map<hamilt::TAC, RI::Tensor<std::complex<double>>>>>* const Hexxc
 #endif
                  ,
                  const double sparse_thr = 1e-10)

source/source_lcao/CMakeLists.txt

@@ -9,19 +9,28 @@ if(ENABLE_LCAO)
         hamilt_lcao.cpp
         module_operator_lcao/operator_lcao.cpp
         module_operator_lcao/veff_lcao.cpp
+        module_operator_lcao/veff_dh.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/overlap_fs.cpp
         module_operator_lcao/ekinetic.cpp
+        module_operator_lcao/ekinetic_fs.cpp
+        module_operator_lcao/ekinetic_dh.cpp
         module_operator_lcao/nonlocal.cpp
+        module_operator_lcao/nonlocal_dh.cpp
+        module_operator_lcao/nonlocal_fs.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/td_pot_hybrid_fs.cpp
         module_operator_lcao/dspin_lcao.cpp
+        module_operator_lcao/dspin_fs.cpp
         module_operator_lcao/dftu_lcao.cpp
-        module_operator_lcao/operator_force_stress_utils.cpp
+        module_operator_lcao/dftu_fs.cpp
+        module_operator_lcao/operator_fs_utils.cpp
         dftu_lcao.cpp
         pulay_fs_center2.cpp
         FORCE_STRESS.cpp

source/source_lcao/hamilt_lcao.h

@@ -5,7 +5,7 @@
 #include "source_cell/klist.h"
 #include "source_cell/module_neighbor/sltk_grid_driver.h"
 #include "source_hamilt/hamilt.h"
-#include "source_lcao/hs_matrix_k.hpp"
+#include "source_hamilt/hs_matrix_k.h"
 #include "source_hamilt/module_hcontainer/hcontainer.h"
 
 #include <memory>

source/source_lcao/module_deepks/test/deepks_test_e_deltabands.cpp

@@ -1,6 +1,6 @@
 #include "deepks_test_runner.h"
 
-#include "source_lcao/hs_matrix_k.hpp"
+#include "source_hamilt/hs_matrix_k.h"
 #include "source_lcao/module_operator_lcao/deepks_lcao.h"
 
 #include <complex>

source/source_lcao/module_operator_lcao/CMakeLists.txt

@@ -5,16 +5,25 @@ add_library(
     op_dftu_lcao.cpp
     meta_lcao.cpp
     veff_lcao.cpp
+    veff_dh.cpp
     deepks_lcao.cpp
     overlap.cpp
+    overlap_fs.cpp
     ekinetic.cpp
+    ekinetic_fs.cpp
+    ekinetic_dh.cpp
     nonlocal.cpp
+    nonlocal_dh.cpp
+    nonlocal_fs.cpp
     td_ekinetic_lcao.cpp
     td_nonlocal_lcao.cpp
     td_pot_hybrid.cpp
+    td_pot_hybrid_fs.cpp
     dspin_lcao.cpp
+    dspin_fs.cpp
     dftu_lcao.cpp
-    operator_force_stress_utils.cpp
+    dftu_fs.cpp
+    operator_fs_utils.cpp
 )
 
 if(ENABLE_COVERAGE)

source/source_lcao/module_operator_lcao/dftu_force_stress.hpp → source/source_lcao/module_operator_lcao/dftu_fs.cpp

@@ -1,4 +1,3 @@
-#pragma once
 #include "dftu_lcao.h"
 #include "source_base/parallel_reduce.h"
 #include "source_base/timer.h"
@@ -426,4 +425,77 @@ void DFTU<OperatorLCAO<TK, TR>>::cal_stress_IJR(const int& iat1,
     }
 }
 
+// explicit member function instantiations
+template void DFTU<OperatorLCAO<double, double>>::cal_force_stress(
+    const bool cal_force, const bool cal_stress,
+    ModuleBase::matrix& force, ModuleBase::matrix& stress);
+template void DFTU<OperatorLCAO<std::complex<double>, double>>::cal_force_stress(
+    const bool cal_force, const bool cal_stress,
+    ModuleBase::matrix& force, ModuleBase::matrix& stress);
+template void DFTU<OperatorLCAO<std::complex<double>, std::complex<double>>>::cal_force_stress(
+    const bool cal_force, const bool cal_stress,
+    ModuleBase::matrix& force, ModuleBase::matrix& stress);
+
+template void DFTU<OperatorLCAO<double, double>>::cal_force_IJR(
+    const int& iat1, const int& iat2,
+    const Parallel_Orbitals* paraV,
+    const std::unordered_map<int, std::vector<double>>& nlm1_all,
+    const std::unordered_map<int, std::vector<double>>& nlm2_all,
+    const std::vector<double>& vu_in,
+    const hamilt::BaseMatrix<double>** dmR_pointer,
+    const int nspin,
+    double* force1, double* force2);
+template void DFTU<OperatorLCAO<std::complex<double>, double>>::cal_force_IJR(
+    const int& iat1, const int& iat2,
+    const Parallel_Orbitals* paraV,
+    const std::unordered_map<int, std::vector<double>>& nlm1_all,
+    const std::unordered_map<int, std::vector<double>>& nlm2_all,
+    const std::vector<double>& vu_in,
+    const hamilt::BaseMatrix<double>** dmR_pointer,
+    const int nspin,
+    double* force1, double* force2);
+template void DFTU<OperatorLCAO<std::complex<double>, std::complex<double>>>::cal_force_IJR(
+    const int& iat1, const int& iat2,
+    const Parallel_Orbitals* paraV,
+    const std::unordered_map<int, std::vector<double>>& nlm1_all,
+    const std::unordered_map<int, std::vector<double>>& nlm2_all,
+    const std::vector<double>& vu_in,
+    const hamilt::BaseMatrix<double>** dmR_pointer,
+    const int nspin,
+    double* force1, double* force2);
+
+template void DFTU<OperatorLCAO<double, double>>::cal_stress_IJR(
+    const int& iat1, const int& iat2,
+    const Parallel_Orbitals* paraV,
+    const std::unordered_map<int, std::vector<double>>& nlm1_all,
+    const std::unordered_map<int, std::vector<double>>& nlm2_all,
+    const std::vector<double>& vu_in,
+    const hamilt::BaseMatrix<double>** dmR_pointer,
+    const int nspin,
+    const ModuleBase::Vector3<double>& dis1,
+    const ModuleBase::Vector3<double>& dis2,
+    double* stress);
+template void DFTU<OperatorLCAO<std::complex<double>, double>>::cal_stress_IJR(
+    const int& iat1, const int& iat2,
+    const Parallel_Orbitals* paraV,
+    const std::unordered_map<int, std::vector<double>>& nlm1_all,
+    const std::unordered_map<int, std::vector<double>>& nlm2_all,
+    const std::vector<double>& vu_in,
+    const hamilt::BaseMatrix<double>** dmR_pointer,
+    const int nspin,
+    const ModuleBase::Vector3<double>& dis1,
+    const ModuleBase::Vector3<double>& dis2,
+    double* stress);
+template void DFTU<OperatorLCAO<std::complex<double>, std::complex<double>>>::cal_stress_IJR(
+    const int& iat1, const int& iat2,
+    const Parallel_Orbitals* paraV,
+    const std::unordered_map<int, std::vector<double>>& nlm1_all,
+    const std::unordered_map<int, std::vector<double>>& nlm2_all,
+    const std::vector<double>& vu_in,
+    const hamilt::BaseMatrix<double>** dmR_pointer,
+    const int nspin,
+    const ModuleBase::Vector3<double>& dis1,
+    const ModuleBase::Vector3<double>& dis2,
+    double* stress);
+
 } // namespace hamilt

source/source_lcao/module_operator_lcao/dftu_lcao.cpp

@@ -687,8 +687,6 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::cal_v_of_u(const std::vector<do
     }
 }
 
-#include "dftu_force_stress.hpp"
-
 template class hamilt::DFTU<hamilt::OperatorLCAO<double, double>>;
 template class hamilt::DFTU<hamilt::OperatorLCAO<std::complex<double>, double>>;
 template class hamilt::DFTU<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>>;