Refactor: new gint module

CMakeLists.txt

@@ -236,6 +236,9 @@ if(ENABLE_LCAO)
     add_compile_definitions(__PEXSI)
     set(CMAKE_CXX_STANDARD 14)
   endif()
+  if(NEW_GINT)
+    add_compile_definitions(__NEW_GINT)
+  endif()
 else()
   set(ENABLE_DEEPKS OFF)
   set(ENABLE_LIBRI OFF)

source/module_base/array_pool.h

@@ -40,7 +40,7 @@ namespace ModuleBase
         : nr(nr_in),
           nc(nc_in)
     {
-        this->ptr_1D = new T[nr * nc];
+        this->ptr_1D = new T[nr * nc]();
         this->ptr_2D = new T*[nr];
         for (int ir = 0; ir < nr; ++ir)
             this->ptr_2D[ir] = &this->ptr_1D[ir * nc];

source/module_base/vector3.h

@@ -365,6 +365,22 @@ template <class T> inline Vector3<T> cross(const Vector3<T> &u, const Vector3<T>
 //	         (u.z * (v.x * w.y - v.y * w.x)));
 // }
 
+// Overload the < operator for sorting
+template <class T> bool operator<(const Vector3<T> &u, const Vector3<T> &v)
+{
+    if (u.x < v.x)
+        return true;
+    if (u.x > v.x)
+        return false;
+    if (u.y < v.y)
+        return true;
+    if (u.y > v.y)
+        return false;
+    if (u.z < v.z)
+        return true;
+    return false;
+}
+
 // whether m1 != m2
 template <class T> inline bool operator!=(const Vector3<T> &u, const Vector3<T> &v)
 {

source/module_elecstate/elecstate_lcao.cpp

@@ -8,7 +8,8 @@
 #include "module_hamilt_lcao/module_gint/grid_technique.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_parameter/parameter.h"
-#include "elecstate_lcao_cal_tau.h"
+
+#include "module_hamilt_lcao/module_gint/new_gint/gint_interface.h"
 
 #include <vector>
 
@@ -59,13 +60,17 @@ void ElecStateLCAO<std::complex<double>>::psiToRho(const psi::Psi<std::complex<d
     //------------------------------------------------------------
 
     ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
+#ifndef __NEW_GINT
     this->gint_k->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
     Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho, PARAM.inp.nspin);
     this->gint_k->cal_gint(&inout);
+#else
+    ModuleGint::cal_gint_rho(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->rho);
+#endif
 
     if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
     {
-        elecstate::lcao_cal_tau_k(gint_k, this->charge);
+        this->cal_tau(psi);
     }
 
     this->charge->renormalize_rho();
@@ -92,15 +97,17 @@ void ElecStateLCAO<double>::psiToRho(const psi::Psi<double>& psi)
     //------------------------------------------------------------
     ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
 
+#ifndef __NEW_GINT 
     this->gint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
-
     Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho, PARAM.inp.nspin);
-
     this->gint_gamma->cal_gint(&inout);
+#else
+    ModuleGint::cal_gint_rho(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->rho);
+#endif
 
     if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
     {
-        elecstate::lcao_cal_tau_gamma(gint_gamma, this->charge);
+        this->cal_tau(psi);
     }
 
     this->charge->renormalize_rho();
@@ -158,17 +165,25 @@ void ElecStateLCAO<double>::dmToRho(std::vector<double*> pexsi_DM, std::vector<d
     }
 
     ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
+#ifndef __NEW_GINT
     this->gint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
     Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho, PARAM.inp.nspin);
     this->gint_gamma->cal_gint(&inout);
+#else
+    ModuleGint::cal_gint_rho(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->rho);
+#endif
     if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
     {
         for (int is = 0; is < PARAM.inp.nspin; is++)
         {
             ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[0], this->charge->nrxx);
         }
+#ifndef __NEW_GINT
         Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau);
         this->gint_gamma->cal_gint(&inout1);
+#else
+        ModuleGint::cal_gint_tau(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->kin_r);
+#endif
     }
 
     this->charge->renormalize_rho();

source/module_elecstate/elecstate_lcao.h

@@ -46,6 +46,7 @@ class ElecStateLCAO : public ElecState
     // virtual void psiToRho(const psi::Psi<double>& psi) override;
     //  return current electronic density rho, as a input for constructing Hamiltonian
     //  const double* getRho(int spin) const override;
+    virtual void cal_tau(const psi::Psi<TK>& psi) override;
 
     // update charge density for next scf step
     // void getNewRho() override;

source/module_elecstate/elecstate_lcao_cal_tau.cpp

@@ -1,56 +1,49 @@
 #include "elecstate_lcao.h"
-#include "elecstate_lcao_cal_tau.h"
+#include "module_hamilt_lcao/module_gint/new_gint/gint_interface.h"
+
 #include "module_base/timer.h"
 
 namespace elecstate
 {
 
 // calculate the kinetic energy density tau, multi-k case
-void lcao_cal_tau_k(Gint_k* gint_k, 
-                    Charge* charge)
+template <>
+void ElecStateLCAO<std::complex<double>>::cal_tau(const psi::Psi<std::complex<double>>& psi)
 {
     ModuleBase::timer::tick("ElecStateLCAO", "cal_tau");
 
     for (int is = 0; is < PARAM.inp.nspin; is++)
     {
-        ModuleBase::GlobalFunc::ZEROS(charge->kin_r[is], charge->nrxx);
+        ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx);
     }
-    Gint_inout inout1(charge->kin_r, Gint_Tools::job_type::tau, PARAM.inp.nspin);
-    gint_k->cal_gint(&inout1);
-
+#ifndef __NEW_GINT
+    Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau, PARAM.inp.nspin);
+    this->gint_k->cal_gint(&inout1);
+#else
+    ModuleGint::cal_gint_tau(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->kin_r);
+#endif
     ModuleBase::timer::tick("ElecStateLCAO", "cal_tau");
     return;
 }
 
 // calculate the kinetic energy density tau, gamma-only case
-void lcao_cal_tau_gamma(Gint_Gamma* gint_gamma,
-                        Charge* charge)
+template <>
+void ElecStateLCAO<double>::cal_tau(const psi::Psi<double>& psi)
 {
     ModuleBase::timer::tick("ElecStateLCAO", "cal_tau");
 
     for (int is = 0; is < PARAM.inp.nspin; is++)
     {
-        ModuleBase::GlobalFunc::ZEROS(charge->kin_r[is], charge->nrxx);
+        ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx);
     }
-    Gint_inout inout1(charge->kin_r, Gint_Tools::job_type::tau, PARAM.inp.nspin);
-    gint_gamma->cal_gint(&inout1);
+#ifndef __NEW_GINT
+    Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau, PARAM.inp.nspin);
+    this->gint_gamma->cal_gint(&inout1);
+#else
+    ModuleGint::cal_gint_tau(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->kin_r);
+#endif
 
     ModuleBase::timer::tick("ElecStateLCAO", "cal_tau");
     return;
 }
-template <> 
-void lcao_cal_tau<double>(Gint_Gamma* gint_gamma, 
-                  Gint_k* gint_k, 
-                  Charge* charge)
-{
-    lcao_cal_tau_gamma(gint_gamma, charge);
-}
-template <> 
-void lcao_cal_tau<complex<double>>(Gint_Gamma* gint_gamma, 
-                    Gint_k* gint_k, 
-                    Charge* charge)
-{
-    lcao_cal_tau_k(gint_k, charge);
-}
-
-} // namespace elecstate
+}

source/module_esolver/esolver_ks_lcao.cpp

@@ -34,7 +34,6 @@
 #include "module_base/global_function.h"
 #include "module_cell/module_neighbor/sltk_grid_driver.h"
 #include "module_elecstate/cal_ux.h"
-#include "module_elecstate/elecstate_lcao_cal_tau.h"
 #include "module_elecstate/module_charge/symmetry_rho.h"
 #include "module_elecstate/occupy.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need DeePKS_init
@@ -933,7 +932,8 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
     // 1) calculate the kinetic energy density tau, sunliang 2024-09-18
     if (PARAM.inp.out_elf[0] > 0)
     {
-        elecstate::lcao_cal_tau<TK>(&(this->GG), &(this->GK), this->pelec->charge);
+        assert(this->psi != nullptr);
+        this->pelec->cal_tau(*(this->psi));
     }
 
     //! 2) call after_scf() of ESolver_KS

source/module_esolver/esolver_ks_lcao.h

@@ -5,6 +5,8 @@
 // for grid integration
 #include "module_hamilt_lcao/module_gint/gint_gamma.h"
 #include "module_hamilt_lcao/module_gint/gint_k.h"
+#include "module_hamilt_lcao/module_gint/new_gint/gint_info.h"
+#include "module_hamilt_lcao/module_gint/new_gint/gint.h"
 #ifdef __DEEPKS
 #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h"
 #endif

source/module_esolver/lcao_before_scf.cpp

@@ -115,6 +115,26 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
                              dpsi_u,
                              d2psi_u,
                              PARAM.inp.nstream);
+
+#ifdef __NEW_GINT
+    auto gint_info = std::make_shared<ModuleGint::GintInfo>(
+        this->pw_big->nbx,
+        this->pw_big->nby,
+        this->pw_big->nbz,
+        this->pw_rho->nx,
+        this->pw_rho->ny,
+        this->pw_rho->nz,
+        0,
+        0,
+        this->pw_big->nbzp_start,
+        this->pw_big->nbx,
+        this->pw_big->nby,
+        this->pw_big->nbzp,
+        orb_.Phi,
+        ucell,
+        this->gd);
+    ModuleGint::Gint::init_gint_info(gint_info);
+#endif
     psi_u.clear();
     psi_u.shrink_to_fit();
     dpsi_u.clear();