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/test/vector3_test.cpp

@@ -64,6 +64,9 @@
  *   - VneV
  *     - overload operator "!=" to assert
  *     - the inequality between two 3d vectors
+ *	 - VltV
+ *     - overload operator "<" to sort
+ * 	   - the "less than" relationship between two 3d vectors
  *   - StdOutV
  *     - overload operator "<<" to print out
  *     - a 3d vectors on standard output
@@ -703,6 +706,22 @@ TEST_F(Vector3Test,VneV)
 	EXPECT_TRUE(wp != w);
 }
 
+TEST_F(Vector3Test, VltV)
+{
+	ModuleBase::Vector3<double> u, up;
+	u.set(da, db, dc);
+	up.set(dc, db, da);
+	EXPECT_TRUE(u < up);
+	ModuleBase::Vector3<float> v, vp;
+	v.set(fa, fb, fc);
+	vp.set(fa, fb, fc);
+	EXPECT_FALSE(v < vp);
+	ModuleBase::Vector3<int> w, wp;
+	w.set(ia, ib, ic);
+	wp.set(ib, ib, ic);
+	EXPECT_TRUE(w < wp);
+}
+
 TEST_F(Vector3Test,StdOutV)
 {
 	// double Vector3

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_basis/module_ao/ORB_atomic.h

@@ -66,6 +66,7 @@ class Numerical_Orbital
 
 	const inline Numerical_Orbital_Lm& PhiLN( const int &L, const int &N)const
 	{ 	
+		assert(this->phiLN != nullptr);
 		return this->phiLN[ this->find_chi(L, N) ];
 	}
 

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/temp_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/temp_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/temp_gint/gint_info.h"
+#include "module_hamilt_lcao/module_gint/temp_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();