update grad and lapl

Author: maki49

source/Makefile.Objects

@@ -718,7 +718,6 @@ OBJS_TENSOR=tensor.o\
 
     OBJS_LR=lr_util.o\
     lr_util_hcontainer.o\
-    lr_util_xc.o\
     AX_parallel.o\
     AX_serial.o\
     dm_trans_parallel.o\

source/module_lr/CMakeLists.txt

@@ -7,7 +7,6 @@ if(ENABLE_LCAO)
   list(APPEND objects
     utils/lr_util.cpp
     utils/lr_util_hcontainer.cpp
-    utils/lr_util_xc.cpp
     AX/AX_parallel.cpp
     AX/AX_serial.cpp
     dm_trans/dm_trans_parallel.cpp

source/module_lr/potentials/pot_hxc_lrtd.cpp

@@ -6,6 +6,7 @@
 #include "module_hamilt_general/module_xc/xc_functional.h"
 #include <set>
 #include "module_lr/utils/lr_util.h"
+#include "module_lr/utils/lr_util_xc.hpp"
 #include "module_io/cube_io.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"    // tmp, for pgrid
 #define FXC_PARA_TYPE const double* const rho, ModuleBase::matrix& v_eff, const std::vector<int>& ispin_op = { 0,0 }

source/module_lr/potentials/xc_kernel.cpp

@@ -3,6 +3,7 @@
 #include "module_parameter/parameter.h"
 #include "module_base/timer.h"
 #include "module_lr/utils/lr_util.h"
+#include "module_lr/utils/lr_util_xc.hpp"
 #ifdef USE_LIBXC
 #include <xc.h>
 #include "module_hamilt_general/module_xc/xc_functional_libxc.h"
@@ -184,7 +185,7 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl
             //             + v2s2.at(ir * 6 + 3) * sigma[ir * 3];   //down-down
             //     }
             //     for (int isig = 0;isig < 3;++isig)
-            //         LR_Util::laplace(v2sigma2_sigma_r.data() + isig * nrxx, v2sigma2_sigma_r.data() + isig * nrxx, *(chg_gs->rhopw), tpiba2);
+            //         LR_Util::lapl(v2sigma2_sigma_r.data() + isig * nrxx, v2sigma2_sigma_r.data() + isig * nrxx, *(chg_gs->rhopw), tpiba2);
             //     // 3. $\nabla^2(v^\sigma)$
             //     std::vector<double> lap_vsigma(3 * nrxx);
             //     for (int ir = 0;ir < nrxx;++ir)
@@ -194,7 +195,7 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl
             //         lap_vsigma[2 * nrxx + ir] = vs.at(ir * 3 + 2) * 2.0;
             //     }
             //     for (int isig = 0;isig < 3;++isig)
-            //         LR_Util::laplace(lap_vsigma.data() + isig * nrxx, lap_vsigma.data() + isig * nrxx, *(chg_gs->rhopw), tpiba2);
+            //         LR_Util::lapl(lap_vsigma.data() + isig * nrxx, lap_vsigma.data() + isig * nrxx, *(chg_gs->rhopw), tpiba2);
             //     // add to v2rho2
             //     BlasConnector::axpy(3 * nrxx, 1.0, v2r2.data(), 1, to_mul_rho_.data(), 1);
             //     BlasConnector::axpy(3 * nrxx, -1.0, div_v2rhosigma_gdrho_r.data(), 1, to_mul_rho_.data(), 1);

source/module_lr/utils/lr_util.h

@@ -12,6 +12,15 @@
 using DAT = container::DataType;
 using DEV = container::DeviceType;
 
+#ifndef TO_COMPLEX_H
+#define TO_COMPLEX_H
+template <typename T> struct ToComplex;
+template <> struct ToComplex<double> { using type = std::complex<double>; };
+template <> struct ToComplex<std::complex<double>> { using type = std::complex<double>; };
+template <> struct ToComplex<float> { using type = std::complex<float>; };
+template <> struct ToComplex<std::complex<float>> { using type = std::complex<float>; };
+#endif
+
 namespace LR_Util
 {
     /// =====================PHYSICS====================
@@ -36,15 +45,7 @@ namespace LR_Util
     std::pair<ModuleBase::matrix, std::vector<std::pair<int, int>>>
         set_ix_map_diagonal(bool mode, int nc, int nv);
 
-    /// operators to calculate XC kernels
-    void grad(const double* rhor,
-        ModuleBase::Vector3<double>* gdr,
-        const ModulePW::PW_Basis& rho_basis,
-        const double& tpiba);
-    void laplace(const double* rhor,
-        double* lapn,
-        const ModulePW::PW_Basis& rho_basis,
-        const double& tpiba2);
+    // Operators to calculate xc kernel have been moved into lr_util_xc.hpp.
     /// =================ALGORITHM====================
 
     //====== newers and deleters========

source/module_lr/utils/lr_util_xc.cpp

@@ -0,0 +1,48 @@
+#pragma once
+#include "lr_util.h"
+namespace LR_Util
+{
+    template<typename T>
+    void grad(const T* rhor,
+        ModuleBase::Vector3<T>* gdr,
+        const ModulePW::PW_Basis& rho_basis,
+        const double& tpiba)
+    {
+        std::vector<typename ToComplex<T>::type> rhog(rho_basis.npw);
+        rho_basis.real2recip(rhor, rhog.data());
+        XC_Functional::grad_rho(rhog.data(), gdr, &rho_basis, tpiba);
+    }
+    template<typename T>
+    void grad(const std::vector<T>& rhor,
+        std::vector<ModuleBase::Vector3<T>>& gdr,
+        const ModulePW::PW_Basis& rho_basis,
+        const double& tpiba)
+    {
+        grad(rhor.data(), gdr.data(), rho_basis, tpiba);
+    }
+
+    template<typename T>
+    void lapl(const T* rhor, T* lapn,
+        const ModulePW::PW_Basis& rho_basis,
+        const double& tpiba2)
+    {
+        ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis.nrxx);
+        std::vector<typename ToComplex<T>::type> rhog(rho_basis.npw);
+        std::vector<T> tmp_rhor(rho_basis.nrxx);
+        rho_basis.real2recip(rhor, rhog.data());
+        for (int i = 0;i < 3;++i)
+        {
+            for (int ig = 0; ig < rho_basis.npw; ig++) { rhog[ig] *= pow(rho_basis.gcar[ig][i], 2); }
+            rho_basis.recip2real(rhog.data(), tmp_rhor.data());
+            for (int ir = 0; ir < rho_basis.nrxx; ir++) { lapn[ir] -= tmp_rhor[ir] * tpiba2; }
+        }
+    }
+    template<typename T>
+    void lapl(const std::vector<T>& rhor,
+        std::vector<T>& lapn,
+        const ModulePW::PW_Basis& rho_basis,
+        const double& tpiba2)
+    {
+        lapl(rhor.data(), lapn.data(), rho_basis, tpiba2);
+    }
+}