Merge branch 'develop' into orb-energy-term

source/Makefile.Objects

@@ -185,16 +185,16 @@ OBJS_CELL=atom_pseudo.o\
     print_cif.o\
 
 OBJS_DEEPKS=LCAO_deepks.o\
-        LCAO_deepks_fgamma.o\
-        LCAO_deepks_fk.o\
+        deepks_fgamma.o\
+        deepks_fk.o\
         LCAO_deepks_odelta.o\
         LCAO_deepks_io.o\
         LCAO_deepks_mpi.o\
         LCAO_deepks_pdm.o\
         LCAO_deepks_psialpha.o\
         LCAO_deepks_torch.o\
         LCAO_deepks_vdelta.o\
-        LCAO_deepks_hmat.o\
+        deepks_hmat.o\
         LCAO_deepks_interface.o\
         orbital_precalc.o\
         orbital_precalc_k.o\
@@ -240,12 +240,15 @@ OBJS_ESOLVER=esolver.o\
     esolver_of.o\
     esolver_of_tool.o\
     esolver_of_interface.o\
+    print_funcs.o\
 
 OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       esolver_ks_lcao_elec.o\
       esolver_ks_lcao_tddft.o\
-      esolver_ks_lcao_tmpfunc.o\
+      dpks_cal_e_delta_band.o\
+      dftu_cal_occup_m.o\
       io_npz.o\
+      set_matrix_grid.o\
 
 OBJS_GINT=gint.o\
       gint_gamma_env.o\

source/module_esolver/CMakeLists.txt

@@ -10,14 +10,17 @@ list(APPEND objects
     esolver_of.cpp
     esolver_of_interface.cpp
     esolver_of_tool.cpp
+    print_funcs.cpp
 )
 if(ENABLE_LCAO)
   list(APPEND objects
       esolver_ks_lcao.cpp
       esolver_ks_lcao_elec.cpp
       esolver_ks_lcao_tddft.cpp
-      esolver_ks_lcao_tmpfunc.cpp
+      dpks_cal_e_delta_band.cpp
       io_npz.cpp
+      set_matrix_grid.cpp
+      dftu_cal_occup_m.cpp
   )
 endif()
 

source/module_esolver/dftu_cal_occup_m.cpp

@@ -0,0 +1,50 @@
+#include "esolver_ks_lcao.h"
+#include "module_hamilt_lcao/module_dftu/dftu.h"
+
+namespace ModuleESolver
+{
+
+    using namespace std;
+
+    //! dftu occupation matrix for gamma only using dm(double)
+    template <>
+    void ESolver_KS_LCAO<double, double>::dftu_cal_occup_m(
+         const int& iter, 
+         const vector<vector<double>>& dm)const
+    {
+		GlobalC::dftu.cal_occup_m_gamma(
+				iter, 
+				dm, 
+				this->p_chgmix->get_mixing_beta(),
+                this->p_hamilt);
+    }
+
+    //! dftu occupation matrix for multiple k-points using dm(complex)
+    template <>
+    void ESolver_KS_LCAO<complex<double>, double>::dftu_cal_occup_m(
+       const int& iter, 
+       const vector<vector<complex<double>>>& dm)const
+    {
+        GlobalC::dftu.cal_occup_m_k(
+                 iter, 
+                 dm, 
+                 this->kv, 
+                 this->p_chgmix->get_mixing_beta(), 
+                 this->p_hamilt);
+    }
+
+    //! dftu occupation matrix
+    template <>
+    void ESolver_KS_LCAO<complex<double>, complex<double>>::dftu_cal_occup_m(
+       const int& iter, 
+       const vector<vector<complex<double>>>& dm)const
+    {
+        GlobalC::dftu.cal_occup_m_k(
+                 iter, 
+                 dm, 
+                 this->kv, 
+                 this->p_chgmix->get_mixing_beta(), 
+                 this->p_hamilt);
+    }
+
+}

source/module_esolver/esolver_ks_lcao_tmpfunc.cpp → source/module_esolver/dpks_cal_e_delta_band.cpp

@@ -1,5 +1,4 @@
 #include "esolver_ks_lcao.h"
-#include "module_hamilt_lcao/module_dftu/dftu.h"
 #ifdef __DEEPKS
 #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
@@ -10,47 +9,6 @@ namespace ModuleESolver
 
     using namespace std;
 
-    //! dftu occupation matrix for gamma only using dm(double)
-    template <>
-    void ESolver_KS_LCAO<double, double>::dftu_cal_occup_m(
-         const int& iter, 
-         const vector<vector<double>>& dm)const
-    {
-		GlobalC::dftu.cal_occup_m_gamma(
-				iter, 
-				dm, 
-				this->p_chgmix->get_mixing_beta(),
-                this->p_hamilt);
-    }
-
-    //! dftu occupation matrix for multiple k-points using dm(complex)
-    template <>
-    void ESolver_KS_LCAO<complex<double>, double>::dftu_cal_occup_m(
-       const int& iter, 
-       const vector<vector<complex<double>>>& dm)const
-    {
-        GlobalC::dftu.cal_occup_m_k(
-                 iter, 
-                 dm, 
-                 this->kv, 
-                 this->p_chgmix->get_mixing_beta(), 
-                 this->p_hamilt);
-    }
-
-    //! dftu occupation matrix
-    template <>
-    void ESolver_KS_LCAO<complex<double>, complex<double>>::dftu_cal_occup_m(
-       const int& iter, 
-       const vector<vector<complex<double>>>& dm)const
-    {
-        GlobalC::dftu.cal_occup_m_k(
-                 iter, 
-                 dm, 
-                 this->kv, 
-                 this->p_chgmix->get_mixing_beta(), 
-                 this->p_hamilt);
-    }
-
 #ifdef __DEEPKS
     template<>
     void ESolver_KS_LCAO<double, double>::dpks_cal_e_delta_band(

source/module_esolver/esolver_ks.cpp

@@ -24,6 +24,8 @@
 #endif
 #include "module_io/json_output/output_info.h"
 
+#include "print_funcs.h" // mohan add 2024-07-27
+
 namespace ModuleESolver
 {
 
@@ -258,7 +260,7 @@ void ESolver_KS<T, Device>::before_all_runners(const Input_para& inp, UnitCell&
 
     this->pw_wfc->collect_local_pw(inp.erf_ecut, inp.erf_height, inp.erf_sigma);
 
-    this->print_wfcfft(inp, GlobalV::ofs_running);
+    Print_functions::print_wfcfft(inp, *this->pw_wfc, GlobalV::ofs_running);
 
     //! 10) initialize the real-space uniform grid for FFT and parallel
     //! distribution of plane waves
@@ -383,79 +385,6 @@ void ESolver_KS<T, Device>::hamilt2density(const int istep, const int iter, cons
     ModuleBase::timer::tick(this->classname, "hamilt2density");
 }
 
-//------------------------------------------------------------------------------
-//! the 6th function of ESolver_KS: print_wfcfft
-//! mohan add 2024-05-11
-//------------------------------------------------------------------------------
-template <typename T, typename Device>
-void ESolver_KS<T, Device>::print_wfcfft(const Input_para& inp, std::ofstream& ofs)
-{
-    ofs << "\n\n\n\n";
-    ofs << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
-           ">>>>"
-        << std::endl;
-    ofs << " |                                                                 "
-           "   |"
-        << std::endl;
-    ofs << " | Setup plane waves of wave functions:                            "
-           "   |"
-        << std::endl;
-    ofs << " | Use the energy cutoff and the lattice vectors to generate the   "
-           "   |"
-        << std::endl;
-    ofs << " | dimensions of FFT grid. The number of FFT grid on each "
-           "processor   |"
-        << std::endl;
-    ofs << " | is 'nrxx'. The number of plane wave basis in reciprocal space "
-           "is   |"
-        << std::endl;
-    ofs << " | different for charege/potential and wave functions. We also set "
-           "   |"
-        << std::endl;
-    ofs << " | the 'sticks' for the parallel of FFT. The number of plane wave "
-           "of  |"
-        << std::endl;
-    ofs << " | each k-point is 'npwk[ik]' in each processor                    "
-           "   |"
-        << std::endl;
-    ofs << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"
-           "<<<<"
-        << std::endl;
-    ofs << "\n\n\n\n";
-    ofs << "\n SETUP PLANE WAVES FOR WAVE FUNCTIONS" << std::endl;
-
-    double ecut = inp.ecutwfc;
-    if (std::abs(ecut - this->pw_wfc->gk_ecut * this->pw_wfc->tpiba2) > 1e-6)
-    {
-        ecut = this->pw_wfc->gk_ecut * this->pw_wfc->tpiba2;
-        ofs << "Energy cutoff for wavefunc is incompatible with nx, ny, nz and "
-               "it will be reduced!"
-            << std::endl;
-    }
-    ModuleBase::GlobalFunc::OUT(ofs, "energy cutoff for wavefunc (unit:Ry)", ecut);
-    ModuleBase::GlobalFunc::OUT(ofs,
-                                "fft grid for wave functions",
-                                this->pw_wfc->nx,
-                                this->pw_wfc->ny,
-                                this->pw_wfc->nz);
-    ModuleBase::GlobalFunc::OUT(ofs, "number of plane waves", this->pw_wfc->npwtot);
-    ModuleBase::GlobalFunc::OUT(ofs, "number of sticks", this->pw_wfc->nstot);
-
-    ofs << "\n PARALLEL PW FOR WAVE FUNCTIONS" << std::endl;
-    ofs << " " << std::setw(8) << "PROC" << std::setw(15) << "COLUMNS(POT)" << std::setw(15) << "PW" << std::endl;
-
-    for (int i = 0; i < GlobalV::NPROC_IN_POOL; ++i)
-    {
-        ofs << " " << std::setw(8) << i + 1 << std::setw(15) << this->pw_wfc->nst_per[i] << std::setw(15)
-            << this->pw_wfc->npw_per[i] << std::endl;
-    }
-
-    ofs << " --------------- sum -------------------" << std::endl;
-    ofs << " " << std::setw(8) << GlobalV::NPROC_IN_POOL << std::setw(15) << this->pw_wfc->nstot << std::setw(15)
-        << this->pw_wfc->npwtot << std::endl;
-    ModuleBase::GlobalFunc::DONE(ofs, "INIT PLANEWAVE");
-}
-
 //------------------------------------------------------------------------------
 //! the 7th function of ESolver_KS: run
 //! mohan add 2024-05-11

source/module_esolver/esolver_ks_lcao.cpp

@@ -878,8 +878,11 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
             // calculation, noted by zhengdy-soc
             if (this->psi != nullptr && (istep % PARAM.inp.out_interval == 0))
             {
-                hamilt::MatrixBlock<TK> h_mat, s_mat;
+                hamilt::MatrixBlock<TK> h_mat;
+                hamilt::MatrixBlock<TK> s_mat;
+
                 this->p_hamilt->matrix(h_mat, s_mat);
+
                 if (hsolver::HSolverLCAO<TK>::out_mat_hs[0])
                 {
                     ModuleIO::save_mat(istep,
@@ -908,7 +911,7 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
 #ifdef __DEEPKS
                 if(GlobalV::deepks_out_labels && GlobalV::deepks_v_delta)
                 {
-                    GlobalC::ld.save_h_mat(h_mat.p,this->ParaV.nloc);
+                    DeePKS_domain::save_h_mat(h_mat.p, this->ParaV.nloc);
                 }
 #endif
             }

source/module_esolver/esolver_ks_lcao_elec.cpp

@@ -37,82 +37,6 @@
 namespace ModuleESolver
 {
 
-template <typename TK, typename TR>
-void ESolver_KS_LCAO<TK, TR>::set_matrix_grid(Record_adj& ra)
-{
-    ModuleBase::TITLE("ESolver_KS_LCAO", "set_matrix_grid");
-    ModuleBase::timer::tick("ESolver_KS_LCAO", "set_matrix_grid");
-
-    // (1) Find adjacent atoms for each atom.
-    GlobalV::SEARCH_RADIUS = atom_arrange::set_sr_NL(GlobalV::ofs_running,
-                                                     GlobalV::OUT_LEVEL,
-                                                     GlobalC::ORB.get_rcutmax_Phi(),
-                                                     GlobalC::ucell.infoNL.get_rcutmax_Beta(),
-                                                     GlobalV::GAMMA_ONLY_LOCAL);
-
-    atom_arrange::search(PARAM.inp.search_pbc,
-                         GlobalV::ofs_running,
-                         GlobalC::GridD,
-                         GlobalC::ucell,
-                         GlobalV::SEARCH_RADIUS,
-                         GlobalV::test_atom_input);
-
-    // ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"SEARCH ADJACENT
-    // ATOMS");
-
-    // (3) Periodic condition search for each grid.
-    double dr_uniform = 0.001;
-    std::vector<double> rcuts;
-    std::vector<std::vector<double>> psi_u;
-    std::vector<std::vector<double>> dpsi_u;
-    std::vector<std::vector<double>> d2psi_u;
-
-    Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, psi_u, dpsi_u, d2psi_u);
-
-    this->GridT.set_pbc_grid(this->pw_rho->nx,
-                             this->pw_rho->ny,
-                             this->pw_rho->nz,
-                             this->pw_big->bx,
-                             this->pw_big->by,
-                             this->pw_big->bz,
-                             this->pw_big->nbx,
-                             this->pw_big->nby,
-                             this->pw_big->nbz,
-                             this->pw_big->nbxx,
-                             this->pw_big->nbzp_start,
-                             this->pw_big->nbzp,
-                             this->pw_rho->ny,
-                             this->pw_rho->nplane,
-                             this->pw_rho->startz_current,
-                             GlobalC::ucell,
-                             dr_uniform,
-                             rcuts,
-                             psi_u,
-                             dpsi_u,
-                             d2psi_u,
-                             PARAM.inp.nstream);
-    psi_u.clear();
-    psi_u.shrink_to_fit();
-    dpsi_u.clear();
-    dpsi_u.shrink_to_fit();
-    d2psi_u.clear();
-    d2psi_u.shrink_to_fit();
-    // (2)For each atom, calculate the adjacent atoms in different cells
-    // and allocate the space for H(R) and S(R).
-    // If k point is used here, allocate HlocR after atom_arrange.
-    Parallel_Orbitals* pv = &this->ParaV;
-    ra.for_2d(*pv, GlobalV::GAMMA_ONLY_LOCAL);
-    if (!GlobalV::GAMMA_ONLY_LOCAL)
-    {
-        // need to first calculae lgd.
-        // using GridT.init.
-        this->GridT.cal_nnrg(pv);
-    }
-
-    ModuleBase::timer::tick("ESolver_KS_LCAO", "set_matrix_grid");
-    return;
-}
-
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
 {

source/module_esolver/esolver_ks_pw.cpp

@@ -595,6 +595,12 @@ void ESolver_KS_PW<T, Device>::hamilt2density(const int istep,
                            GlobalV::use_uspp,
                            GlobalV::RANK_IN_POOL,
                            GlobalV::NPROC_IN_POOL,
+
+                           hsolver::DiagoIterAssist<T, Device>::SCF_ITER,
+                           hsolver::DiagoIterAssist<T, Device>::need_subspace,
+                           hsolver::DiagoIterAssist<T, Device>::PW_DIAG_NMAX,
+                           hsolver::DiagoIterAssist<T, Device>::PW_DIAG_THR,
+
                            false);
 
 
@@ -1080,6 +1086,12 @@ void ESolver_KS_PW<T, Device>::hamilt2estates(const double ethr) {
                            GlobalV::use_uspp,
                            GlobalV::RANK_IN_POOL,
                            GlobalV::NPROC_IN_POOL,
+
+                           hsolver::DiagoIterAssist<T, Device>::SCF_ITER,
+                           hsolver::DiagoIterAssist<T, Device>::need_subspace,
+                           hsolver::DiagoIterAssist<T, Device>::PW_DIAG_NMAX,
+                           hsolver::DiagoIterAssist<T, Device>::PW_DIAG_THR,
+
                            true);
     } else {
         ModuleBase::WARNING_QUIT("ESolver_KS_PW",