fully support new gint module

Author: dzzz2001

source/module_esolver/lcao_before_scf.cpp

@@ -69,6 +69,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     Gint_Tools::init_orb(dr_uniform, rcuts, ucell, orb_, psi_u, dpsi_u, d2psi_u);
 
     //! 5) set periodic boundary conditions
+#ifndef __NEW_GINT
     this->GridT.set_pbc_grid(this->pw_rho->nx,
                              this->pw_rho->ny,
                              this->pw_rho->nz,
@@ -92,9 +93,17 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
                              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();
+    LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell, orb_, *this->pw_rho, *this->pw_big);
 
     //! 6) prepare grid integral
-#ifdef __NEW_GINT
+#else
     auto gint_info = std::make_shared<ModuleGint::GintInfo>(
         this->pw_big->nbx,
         this->pw_big->nby,
@@ -114,22 +123,12 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     ModuleGint::Gint::set_gint_info(gint_info);
 #endif
 
-    psi_u.clear();
-    psi_u.shrink_to_fit();
-    dpsi_u.clear();
-    dpsi_u.shrink_to_fit();
-    d2psi_u.clear();
-    d2psi_u.shrink_to_fit();
-
     // 7) 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.
     this->RA.for_2d(ucell, this->gd, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
 
-    // 8) after ions move, prepare grid in Gint
-    LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell, orb_, *this->pw_rho, *this->pw_big);
-
-    // 9) initialize the Hamiltonian operators
+    // 8) initialize the Hamiltonian operators
     // if atom moves, then delete old pointer and add a new one
     if (this->p_hamilt != nullptr)
     {
@@ -169,7 +168,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
 
 
 #ifdef __MLALGO
-    // 10) for each ionic step, the overlap <phi|alpha> must be rebuilt
+    // 9) for each ionic step, the overlap <phi|alpha> must be rebuilt
     // since it depends on ionic positions
     if (PARAM.globalv.deepks_setorb)
     {
@@ -198,7 +197,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
 #endif
 
-    // 11) prepare sc calculation
+    // 10) prepare sc calculation
     if (PARAM.inp.sc_mag_switch)
     {
         spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();
@@ -217,7 +216,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
                    this->pelec);
     }
 
-    // 12) set xc type before the first cal of xc in pelec->init_scf
+    // 11) set xc type before the first cal of xc in pelec->init_scf
     // Peize Lin add 2016-12-03
 #ifdef __EXX
     if (PARAM.inp.calculation != "nscf")
@@ -233,10 +232,10 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
 #endif
 
-    // 13) init_scf, should be before_scf? mohan add 2025-03-10
+    // 12) init_scf, should be before_scf? mohan add 2025-03-10
     this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
 
-    // 14) initalize DMR
+    // 13) initalize DMR
     // DMR should be same size with Hamiltonian(R)
     dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)
         ->get_DM()
@@ -247,26 +246,26 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     this->ld.init_DMR(ucell, orb_, this->pv, this->gd);
 #endif
 
-    // 15) two cases are considered:
+    // 14) 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
     if (istep > 0)
     {
         dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->cal_DMR();
     }
 
-    // 16) the electron charge density should be symmetrized,
+    // 15) the electron charge density should be symmetrized,
     // here is the initialization
     Symmetry_rho srho;
     for (int is = 0; is < PARAM.inp.nspin; is++)
     {
         srho.begin(is, this->chr, this->pw_rho, ucell.symm);
     }
 
-    // 17) why we need to set this sentence? mohan add 2025-03-10
+    // 16) why we need to set this sentence? mohan add 2025-03-10
     this->p_hamilt->non_first_scf = istep;
 
-    // 18) update of RDMFT, added by jghan
+    // 17) update of RDMFT, added by jghan
     if (PARAM.inp.rdmft == true)
     {
         // necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp

source/module_esolver/lcao_others.cpp

@@ -93,14 +93,14 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
                          PARAM.inp.test_atom_input);
 
     // (3) Periodic condition search for each grid.
+#ifndef __NEW_GINT
     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, ucell, orb_, psi_u, dpsi_u, d2psi_u);
-
     this->GridT.set_pbc_grid(this->pw_rho->nx,
                              this->pw_rho->ny,
                              this->pw_rho->nz,
@@ -124,7 +124,16 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
                              dpsi_u,
                              d2psi_u,
                              PARAM.inp.nstream);
-    #ifdef __NEW_GINT
+        
+    psi_u.clear();
+    psi_u.shrink_to_fit();
+    dpsi_u.clear();
+    dpsi_u.shrink_to_fit();
+    d2psi_u.clear();
+    d2psi_u.shrink_to_fit();
+    // prepare grid in Gint
+    LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell, orb_, *this->pw_rho, *this->pw_big);
+#else
     auto gint_info = std::make_shared<ModuleGint::GintInfo>(
         this->pw_big->nbx,
         this->pw_big->nby,
@@ -142,14 +151,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
         ucell,
         this->gd);
     ModuleGint::Gint::set_gint_info(gint_info);
-    #endif
-    
-    psi_u.clear();
-    psi_u.shrink_to_fit();
-    dpsi_u.clear();
-    dpsi_u.shrink_to_fit();
-    d2psi_u.clear();
-    d2psi_u.shrink_to_fit();
+#endif
 
     // (2)For each atom, calculate the adjacent atoms in different cells
     // and allocate the space for H(R) and S(R).
@@ -206,9 +208,6 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
         }
     }
 
-    // prepare grid in Gint
-    LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell, orb_, *this->pw_rho, *this->pw_big);
-
     // init Hamiltonian
     if (this->p_hamilt != nullptr)
     {

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h

@@ -50,7 +50,9 @@ class Veff<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
         this->cal_type = calculation_type::lcao_gint;
 
         this->initialize_HR(ucell_in, GridD_in);
+#ifndef __NEW_GINT
         GK_in->initialize_pvpR(*ucell_in, GridD_in, nspin);
+#endif
     }
     /**
      * @brief Construct a new Veff object for Gamma-only calculation
@@ -69,8 +71,9 @@ class Veff<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
     {
         this->cal_type = calculation_type::lcao_gint;
         this->initialize_HR(ucell_in, GridD_in);
-
+#ifndef __NEW_GINT
         GG_in->initialize_pvpR(*ucell_in, GridD_in, nspin);
+#endif
     }
 
     ~Veff<OperatorLCAO<TK, TR>>(){};

source/module_lr/esolver_lrtd_lcao.cpp

@@ -397,6 +397,7 @@ LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu
     this->gint_->gridt = &this->gt_;
 
     // (3) Periodic condition search for each grid.
+#ifndef __NEW_GINT
     double dr_uniform = 0.001;
     std::vector<double> rcuts;
     std::vector<std::vector<double>> psi_u;
@@ -451,7 +452,25 @@ LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu
         &ucell,
         &orb);
     this->gint_->initialize_pvpR(ucell, &this->gd, 1); // always use nspin=1 for transition density
-
+#else
+    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::set_gint_info(gint_info);
+#endif
     // if EXX from scratch, init 2-center integral and calculate Cs, Vs 
 #ifdef __EXX
     if ((xc_kernel == "hf" || xc_kernel == "hse") && this->input.lr_solver != "spectrum")

source/module_lr/esolver_lrtd_lcao.h

@@ -17,6 +17,7 @@
 #include "module_elecstate/module_dm/density_matrix.h"
 #include "module_lr/potentials/pot_hxc_lrtd.h"
 #include "module_lr/hamilt_casida.h"
+#include "module_hamilt_lcao/module_gint/temp_gint/gint.h"
 #ifdef __EXX
 // #include <RI/physics/Exx.h>
 #include "module_ri/Exx_LRI.h"
@@ -34,6 +35,9 @@ namespace LR
         ESolver_LR(const Input_para& inp, UnitCell& ucell);
         ~ESolver_LR() {
             delete this->psi_ks;
+#ifdef __NEW_GINT
+            ModuleGint::Gint::set_gint_info(nullptr);
+#endif
         }
 
         ///input: input, call, basis(LCAO), psi(ground state), elecstate

source/module_lr/lr_spectrum.cpp

@@ -6,6 +6,7 @@
 #include "module_lr/utils/lr_util.h"
 #include "module_lr/utils/lr_util_hcontainer.h"
 #include "module_lr/utils/lr_util_print.h"
+#include "module_hamilt_lcao/module_gint/temp_gint/gint_interface.h"
 
 template <typename T>
 elecstate::DensityMatrix<T, T> LR::LR_Spectrum<T>::cal_transition_density_matrix(const int istate, const T* X_in, const bool need_R)
@@ -34,14 +35,6 @@ elecstate::DensityMatrix<T, T> LR::LR_Spectrum<T>::cal_transition_density_matrix
     return DM_trans;
 }
 
-template<typename T>
-void LR::LR_Spectrum<T>::cal_gint_rho(double** rho, const int& nrxx)
-{
-    ModuleBase::GlobalFunc::ZEROS(rho[0], nrxx);
-    Gint_inout inout_rho(rho, Gint_Tools::job_type::rho, 1, false);
-    this->gint->cal_gint(&inout_rho);
-}
-
 inline void check_sum_rule(const double& osc_tot)
 {
     if (std::abs(osc_tot - 1.0) > 1e-3) {
@@ -59,12 +52,16 @@ ModuleBase::Vector3<double> LR::LR_Spectrum<double>::cal_transition_dipole_istat
     const elecstate::DensityMatrix<double, double>& DM_trans = this->cal_transition_density_matrix(istate);
     for (int is = 0;is < this->nspin_x;++is)
     {
-        this->gint->transfer_DM2DtoGrid({ DM_trans.get_DMR_vector().at(is) });
-
         // 2. transition density
         double** rho_trans;
         LR_Util::_allocate_2order_nested_ptr(rho_trans, 1, this->rho_basis.nrxx);
+#ifndef __NEW_GINT
+        this->gint->transfer_DM2DtoGrid({ DM_trans.get_DMR_vector().at(is) });
         this->cal_gint_rho(rho_trans, this->rho_basis.nrxx);
+#else
+        ModuleBase::GlobalFunc::ZEROS(rho_trans[0], this->rho_basis.nrxx);
+        ModuleGint::cal_gint_rho({ DM_trans.get_DMR_vector().at(is) }, 1, rho_trans);
+#endif
 
         // 3. transition dipole moment
         for (int ir = 0; ir < rho_basis.nrxx; ++ir)
@@ -108,14 +105,24 @@ ModuleBase::Vector3<std::complex<double>> LR::LR_Spectrum<std::complex<double>>:
 
         // real part
         LR_Util::get_DMR_real_imag_part(DM_trans, DM_trans_real_imag, ucell.nat, 'R');
+#ifndef __NEW_GINT
         this->gint->transfer_DM2DtoGrid(DM_trans_real_imag.get_DMR_vector());
         this->cal_gint_rho(rho_trans_real, this->rho_basis.nrxx);
+#else
+        ModuleBase::GlobalFunc::ZEROS(rho_trans_real[0], this->rho_basis.nrxx);
+        ModuleGint::cal_gint_rho(DM_trans_real_imag.get_DMR_vector(), 1, rho_trans_real);
+#endif
         // LR_Util::print_grid_nonzero(rho_trans_real[0], this->rho_basis.nrxx, 10, "rho_trans");
 
         // imag part
         LR_Util::get_DMR_real_imag_part(DM_trans, DM_trans_real_imag, ucell.nat, 'I');
+#ifndef __NEW_GINT
         this->gint->transfer_DM2DtoGrid(DM_trans_real_imag.get_DMR_vector());
         this->cal_gint_rho(rho_trans_imag, this->rho_basis.nrxx);
+#else
+        ModuleBase::GlobalFunc::ZEROS(rho_trans_imag[0], this->rho_basis.nrxx);
+        ModuleGint::cal_gint_rho(DM_trans_real_imag.get_DMR_vector(), 1, rho_trans_imag);
+#endif
         // LR_Util::print_grid_nonzero(rho_trans_imag[0], this->rho_basis.nrxx, 10, "rho_trans");
 
         // 3. transition dipole moment