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

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_k_sparse1.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;
 

source/module_hamilt_lcao/module_gint/kernels/cuda/cuda_tools.cu

@@ -4,27 +4,7 @@
 
 #include "cuda_tools.cuh"
 
-cudaError_t check(cudaError_t result, const char *const func, const char *const file, const int line)
-{
-    if (result != cudaSuccess)
-    {
-        fprintf(stderr, "CUDA Runtime Error at %s:%d code=%s \"%s\" \n", file, line, cudaGetErrorString(result), func);
-        exit(EXIT_FAILURE);
-    }
-    return result;
-}
-cudaError_t __checkCudaLastError(const char *file, const int line)
-{
-    cudaError_t result = cudaGetLastError();
-    if (result != cudaSuccess)
-    {
-        fprintf(stderr, "%s(%i) : getLastCudaError():%s\n", file, line, cudaGetErrorString(result));
-        assert(result == cudaSuccess);
-    }
-    return result;
-}
-
-void dump_cuda_array_to_file(double* cuda_array,
+void dump_cuda_array_to_file(const double* cuda_array,
                              int width,
                              int hight,
                              const std::string& filename)

source/module_hamilt_lcao/module_gint/kernels/cuda/cuda_tools.cuh

@@ -9,21 +9,47 @@
 #include <iostream>
 #include <sstream>
 
-#define checkCuda(val) check(val, #val, __FILE__, __LINE__)
-#define checkCudaLastError() __checkCudaLastError(__FILE__, __LINE__)
+#define checkCuda(val) check((val), #val, __FILE__, __LINE__)
+#define checkCudaLastError() __getLastCudaError(__FILE__, __LINE__)
 
-cudaError_t check(cudaError_t result, const char *const func, const char *const file, const int line);
-cudaError_t __checkCudaLastError(const char *file, const int line);
+inline void check(cudaError_t result, char const *const func, const char *const file,
+           int const line) {
+  if (result) {
+    fprintf(stderr, "CUDA error at %s:%d code=%d(%s) \"%s\" \n", file, line,
+            static_cast<unsigned int>(result), cudaGetErrorString(result), func);
+    exit(EXIT_FAILURE);
+  }
+}
+
+inline void __getLastCudaError(const char *file,
+                               const int line) 
+{
+  cudaError_t err = cudaGetLastError();
+
+  if (cudaSuccess != err) {
+    fprintf(stderr,
+            "%s(%i) : getLastCudaError() CUDA error :"
+            " (%d) %s.\n",
+            file, line, static_cast<int>(err),
+            cudaGetErrorString(err));
+    exit(EXIT_FAILURE);
+  }
+}
 
-void dump_cuda_array_to_file(double* cuda_array,
+static inline int ceildiv(int x, int y)
+{
+    return (x + y - 1) / y;
+}
+
+void dump_cuda_array_to_file(const double* cuda_array,
                              int width,
                              int hight,
                              const std::string& filename);
 
-inline int ceil_div(int a, int b)
-{
-    return (a + b - 1) / b;
-}
+// inline int ceil_div(int a, int b)
+// {
+//     return (a + b - 1) / b;
+// }
 
 /*
  * @brief: A simple wrapper for cudaMalloc and cudaFree, sync and async CUDA

source/module_hamilt_lcao/module_gint/kernels/cuda/vbatch_matrix_mul.cuh

@@ -320,11 +320,6 @@ static __global__ void vbatched_gemm_kernel(int* M,
                                           alpha_tmp);
 }
 
-static inline int ceildiv(int x, int y)
-{
-    return (x + y - 1) / y;
-}
-
 /**
  * Performs a batched matrix multiplication using the vbatched_gemm_impl
  * function.