new gint module (This is a big PR that replaces old integral codes) (#6296)

* version 0.0

* change gemm function

* enable gint_vl_gpu

* add new gemm function

* support rho calculation

* enable fvl calculation

* small change

* enable vlocal-metagga and tau calculation

* add set ddphi

* enable gint_vl_nspin4_gpu

* add gint_vl_metagga_nspin4_gpu

* enable fvl_meta_gpu

* optimize dgemm_vbatch

* set streams num

* small fix

* optimize gint_atom

* delete virtual function

* add cal_env to new gint module

* renew some notes

* simplify cuda_mem_wrapper.h

* update cuda_mem_wrapper.h

* fix get_wf_lcao.cpp

* fix assert error

* add gint_info init for lcao_others

* refactor output_dHR

* add gint_dvlocal to new gint module

* add namespace

* remove some shared_ptr

* rename phi_mul_phi_vldr3

* small modification

* remove some outdated comments

* small modification

* small modification

* renew gint interface in rdmft related code

* renew gint interface in get_pchg_lcao.cpp

* fully support new gint module

* reduce header dependency

* remove redundant ";"

* fix cmakelist

* add checkcudalasterror

* small modification

* add some timer

* Revert "change gemm function"

This reverts commit ab61eb4.

* fix a bug

* fix bugs

* delete unused declaration

* make new gint module the default compilation option

* add more specific debug info

* fix bug of gint_env_k.cpp

* fix a bug about mixing_dmr

* small modification

* fix bug of mixing_dmr

* modify comment

* replace pointer with std::vector in Numerical_Orbital

* change copy assignment function in intarray

* fix error in module_lr

* replace module with source

* include missing header

* fix a bug

* update some result.ref

* fix compilation error when MPI is disabled

* add makefile support

* fix compilation error

Author: dzzz2001

CMakeLists.txt

@@ -252,8 +252,8 @@ if(ENABLE_LCAO)
     add_compile_definitions(__PEXSI)
     set(CMAKE_CXX_STANDARD 14)
   endif()
-  if(NEW_GINT)
-    add_compile_definitions(__NEW_GINT)
+  if(OLD_GINT)
+    add_compile_definitions(__OLD_GINT)
   endif()
 else()
   set(ENABLE_MLALGO OFF)

source/Makefile.Objects

@@ -61,6 +61,7 @@ VPATH=./src_global:\
 ./module_hamilt_lcao/module_deltaspin:\
 ./module_hamilt_lcao/hamilt_lcaodft/operator_lcao:\
 ./module_hamilt_lcao/module_gint:\
+./module_hamilt_lcao/module_gint/temp_gint:\
 ./module_relax:\
 ./source_hamilt/module_vdw:\
 ./module_io:\
@@ -273,13 +274,13 @@ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       lcao_others.o\
       esolver_dm2rho.o\
 
-OBJS_GINT=gint.o\
+OBJS_GINT=gint_old.o\
       gint_gamma_env.o\
       gint_gamma_vl.o\
-      gint_fvl.o\
-      gint_rho.o\
-      gint_tau.o\
-      gint_vl.o\
+      gint_fvl_old.o\
+      gint_rho_old.o\
+      gint_tau_old.o\
+      gint_vl_old.o\
       gint_k_env.o\
       gint_k_sparse1.o\
       gint_k_pvpr.o\
@@ -298,6 +299,30 @@ OBJS_GINT=gint.o\
       cal_ddpsir_ylm.o\
       mult_psi_dmr.o\
       init_orb.o\
+      batch_biggrid.o\
+      big_grid.o\
+      biggrid_info.o\
+      divide_info.o\
+      gint_atom.o\
+      gint_common.o\
+      gint_dvlocal.o\
+      gint_env_gamma.o\
+      gint_env_k.o\
+      gint_fvl_meta.o\
+      gint_fvl.o\
+      gint_info.o\
+      gint_interface.o\
+      gint_rho.o\
+      gint_tau.o\
+      gint_vl_metagga_nspin4.o\
+      gint_vl_metagga.o\
+      gint_vl_nspin4.o\
+      gint_vl.o\
+      gint.o\
+      localcell_info.o\
+      phi_operator.o\
+      set_ddphi.o\
+      unitcell_info.o\
 
 OBJS_HAMILT=hamilt_pw.o\
     hamilt_sdft_pw.o\

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp

@@ -68,7 +68,7 @@ void Veff<OperatorLCAO<double, double>>::contributeHR()
     double* vr_eff1 = this->pot->get_effective_v(this->current_spin);
     double* vofk_eff1 = this->pot->get_effective_vofk(this->current_spin);
 
-#ifndef __NEW_GINT
+#ifdef __OLD_GINT
     if(XC_Functional::get_ked_flag())
     {
         Gint_inout inout(vr_eff1, vofk_eff1, Gint_Tools::job_type::vlocal_meta);
@@ -113,7 +113,7 @@ void Veff<OperatorLCAO<std::complex<double>, double>>::contributeHR()
     double* vr_eff1 = this->pot->get_effective_v(this->current_spin);
     double* vofk_eff1 = this->pot->get_effective_vofk(this->current_spin);
 
-#ifndef __NEW_GINT
+#ifdef __OLD_GINT
     // if you change the place of the following code,
     // rememeber to delete the #include
     if(XC_Functional::get_ked_flag())
@@ -155,7 +155,7 @@ void Veff<OperatorLCAO<std::complex<double>, std::complex<double>>>::contributeH
     ModuleBase::TITLE("Veff", "contributeHR");
     ModuleBase::timer::tick("Veff", "contributeHR");
 
-#ifndef __NEW_GINT
+#ifdef __OLD_GINT
     double* vr_eff1 = nullptr;
     double* vofk_eff1 = nullptr;
     for (int is = 0; is < 4; is++)
@@ -187,19 +187,15 @@ void Veff<OperatorLCAO<std::complex<double>, std::complex<double>>>::contributeH
         if(XC_Functional::get_ked_flag())
         {
             vofk_eff[is] = this->pot->get_effective_vofk(is);
-            if(is == 3)
-            {
-                ModuleGint::cal_gint_vl_metagga(vr_eff, vofk_eff, this->hR);
-            }
-        }
-        else
-        {
-            if(is == 3)
-            {
-                ModuleGint::cal_gint_vl(vr_eff, this->hR);
-            }
         }
     }
+    if(XC_Functional::get_ked_flag())
+    {
+        ModuleGint::cal_gint_vl_metagga(vr_eff, vofk_eff, this->hR);
+    } else
+    {
+        ModuleGint::cal_gint_vl(vr_eff, this->hR);
+    }
 #endif
 
     ModuleBase::timer::tick("Veff", "contributeHR");

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);
+#ifdef __OLD_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);
-
+#ifdef __OLD_GINT
         GG_in->initialize_pvpR(*ucell_in, GridD_in, nspin);
+#endif
     }
 
     ~Veff<OperatorLCAO<TK, TR>>(){};

source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_gint.hpp

@@ -20,7 +20,7 @@ namespace PulayForceStress
     {
         const int nspin = PARAM.inp.nspin;
 
-#ifndef __NEW_GINT
+#ifdef __OLD_GINT
         if (set_dmr_gint) { gint.transfer_DM2DtoGrid(dm.get_DMR_vector()); }    // 2d block to grid
         for (int is = 0; is < nspin; ++is)
         {

source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.cpp

@@ -2,6 +2,7 @@
 
 #include "module_parameter/parameter.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
+#include "module_hamilt_lcao/module_gint/temp_gint/gint_interface.h"
 #include <vector>
 
 void sparse_format::cal_dS(const UnitCell& ucell,
@@ -49,7 +50,6 @@ delete[] fsr_dh.DHloc_fixedR_y;
 delete[] fsr_dh.DHloc_fixedR_z;
 return;
 }
-
 void sparse_format::cal_dH(const UnitCell& ucell,
                            const Parallel_Orbitals& pv,
                            LCAO_HS_Arrays& HS_Arrays,
@@ -58,6 +58,7 @@ void sparse_format::cal_dH(const UnitCell& ucell,
                            const LCAO_Orbitals& orb,
                            const int& current_spin,
                            const double& sparse_thr,
+                           const ModuleBase::matrix& v_eff,
                            Gint_k& gint_k)
 {
     ModuleBase::TITLE("sparse_format", "cal_dH");
@@ -106,8 +107,38 @@ void sparse_format::cal_dH(const UnitCell& ucell,
     delete[] fsr_dh.DHloc_fixedR_y;
     delete[] fsr_dh.DHloc_fixedR_z;
 
-    gint_k.cal_dvlocal_R_sparseMatrix(current_spin, sparse_thr, HS_Arrays, &pv, ucell, grid);
-
+    if(PARAM.inp.nspin==2)
+    {
+#ifdef __OLD_GINT
+        gint_k.allocate_pvdpR();
+        // note: some MPI process will not have grids when MPI cores are too
+        // many, v_eff in these processes are empty
+        const double* vr_eff1
+            = v_eff.nc * v_eff.nr > 0 ? &(v_eff(current_spin, 0)) : nullptr;
+
+        if (!PARAM.globalv.gamma_only_local) 
+        {
+            if (PARAM.inp.vl_in_h) 
+            {
+                Gint_inout inout(vr_eff1,
+                                 current_spin,
+                                 Gint_Tools::job_type::dvlocal);
+                gint_k.cal_gint(&inout);
+            }
+        }
+        gint_k.cal_dvlocal_R_sparseMatrix(current_spin, sparse_thr, HS_Arrays, &pv, ucell, grid);
+        gint_k.destroy_pvdpR();
+#else
+        const double* vr_eff1
+            = v_eff.nc * v_eff.nr > 0 ? &(v_eff(current_spin, 0)) : nullptr;
+        if (!PARAM.globalv.gamma_only_local) 
+        {
+            ModuleGint::cal_dvlocal_R_sparseMatrix(
+                PARAM.inp.nspin, PARAM.globalv.npol, current_spin, PARAM.globalv.nlocal,
+                sparse_thr, vr_eff1, pv, ucell, grid, HS_Arrays);
+        }
+#endif
+    }
     return;
 }
 

source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.h

@@ -19,6 +19,7 @@ void cal_dH(const UnitCell& ucell,
             const LCAO_Orbitals& orb,
             const int& current_spin,
             const double& sparse_thr,
+            const ModuleBase::matrix& v_eff,
             Gint_k& gint_k);
 
 // calculated the derivative of the overlap matrix: <phi|dphi>

source/module_hamilt_lcao/module_gint/CMakeLists.txt

@@ -2,13 +2,13 @@
 if(ENABLE_LCAO)
 
 list(APPEND objects
-    gint.cpp
+    gint_old.cpp
     gint_gamma_env.cpp
     gint_gamma_vl.cpp
-    gint_fvl.cpp
-    gint_rho.cpp
-    gint_tau.cpp
-    gint_vl.cpp
+    gint_fvl_old.cpp
+    gint_rho_old.cpp
+    gint_tau_old.cpp
+    gint_vl_old.cpp
     gint_k_env.cpp
     gint_k_sparse1.cpp
     gint_k_pvpr.cpp
@@ -29,7 +29,7 @@ list(APPEND objects
     init_orb.cpp
 )
 
-if(NEW_GINT)
+if(NOT DEFINED OLD_GINT)
   list(APPEND objects
       temp_gint/biggrid_info.cpp
       temp_gint/big_grid.cpp
@@ -45,13 +45,34 @@ if(NEW_GINT)
       temp_gint/gint_tau.cpp
       temp_gint/gint_fvl.cpp
       temp_gint/gint_fvl_meta.cpp
+      temp_gint/gint_env_gamma.cpp
+      temp_gint/gint_env_k.cpp
+      temp_gint/gint_dvlocal.cpp
       temp_gint/localcell_info.cpp
       temp_gint/phi_operator.cpp
       temp_gint/set_ddphi.cpp
       temp_gint/unitcell_info.cpp
       temp_gint/gint_common.cpp
       temp_gint/gint_interface.cpp
       )
+  if(USE_CUDA)
+    list(APPEND objects
+        temp_gint/kernel/gint_gpu_vars.cpp
+        temp_gint/kernel/phi_operator_gpu.cu
+        temp_gint/kernel/phi_operator_kernel.cu
+        temp_gint/kernel/set_const_mem.cu
+        temp_gint/batch_biggrid.cpp
+        temp_gint/gint_vl_gpu.cpp
+        temp_gint/gint_rho_gpu.cpp
+        temp_gint/gint_fvl_gpu.cpp
+        temp_gint/gint_vl_metagga_gpu.cpp
+        temp_gint/gint_vl_nspin4_gpu.cpp
+        temp_gint/gint_vl_metagga_nspin4_gpu.cpp
+        temp_gint/gint_tau_gpu.cpp
+        temp_gint/gint_fvl_meta_gpu.cpp
+        temp_gint/kernel/dgemm_vbatch.cu
+    )
+  endif()
 endif()
 
 if(USE_CUDA)

source/module_hamilt_lcao/module_gint/gint_fvl.cpp renamed to source/module_hamilt_lcao/module_gint/gint_fvl_old.cpp

@@ -337,7 +337,6 @@ void Gint_k::cal_dvlocal_R_sparseMatrix(const int& current_spin,
     std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>>
         pvdpRz_soc_sparseMatrix;
 
-    int lgd = 0;
     double temp_value_double;
     std::complex<double> temp_value_complex;