move xc-dependent functions to a new file

Author: maki49

source/Makefile.Objects

@@ -718,6 +718,7 @@ 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,6 +7,7 @@ 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/utils/lr_util.cpp

@@ -62,8 +62,9 @@ namespace LR_Util
     template<>
     void matsym<double>(const double* in, const int n, const Parallel_2D& pmat, double* out)
     {
-        for (int i = 0;i < pmat.get_local_size();++i) {out[i] = in[i];
-}
+        for (int i = 0;i < pmat.get_local_size();++i) {            
+out[i] = in[i];
+        }
         const double alpha = 0.5, beta = 0.5;
         const int i1 = 1;
         pdtran_(&n, &n, &alpha, in, &i1, &i1, pmat.desc, &beta, out, &i1, &i1, pmat.desc);
@@ -72,17 +73,19 @@ namespace LR_Util
     void matsym<double>(double* inout, const int n, const Parallel_2D& pmat)
     {
         std::vector<double> tmp(n * n);
-        for (int i = 0;i < pmat.get_local_size();++i) {tmp[i] = inout[i];
-}
+        for (int i = 0;i < pmat.get_local_size();++i) {
+            tmp[i] = inout[i];
+        }
         const double alpha = 0.5, beta = 0.5;
         const int i1 = 1;
         pdtran_(&n, &n, &alpha, tmp.data(), &i1, &i1, pmat.desc, &beta, inout, &i1, &i1, pmat.desc);
     }
     template<>
     void matsym<std::complex<double>>(const std::complex<double>* in, const int n, const Parallel_2D& pmat, std::complex<double>* out)
     {
-        for (int i = 0;i < pmat.get_local_size();++i) {out[i] = in[i];
-}
+        for (int i = 0;i < pmat.get_local_size();++i) {
+            out[i] = in[i];
+        }
         const std::complex<double> alpha(0.5, 0.0), beta(0.5, 0.0);
         const int i1 = 1;
         pztranc_(&n, &n, &alpha, in, &i1, &i1, pmat.desc, &beta, out, &i1, &i1, pmat.desc);
@@ -91,8 +94,9 @@ namespace LR_Util
     void matsym<std::complex<double>>(std::complex<double>* inout, const int n, const Parallel_2D& pmat)
     {
         std::vector<std::complex<double>> tmp(n * n);
-        for (int i = 0;i < pmat.get_local_size();++i) {tmp[i] = inout[i];
-}
+        for (int i = 0;i < pmat.get_local_size();++i) {
+            tmp[i] = inout[i];
+        }
         const std::complex<double> alpha(0.5, 0.0), beta(0.5, 0.0);
         const int i1 = 1;
         pztranc_(&n, &n, &alpha, tmp.data(), &i1, &i1, pmat.desc, &beta, inout, &i1, &i1, pmat.desc);
@@ -101,88 +105,100 @@ namespace LR_Util
     container::Tensor mat2ten_double(ModuleBase::matrix& m)
     {
         container::Tensor t(DAT::DT_DOUBLE, DEV::CpuDevice, { m.nr, m.nc });
-        for (int i = 0;i < t.NumElements();++i) {t.data<double>()[i] = m.c[i];
-}
+        for (int i = 0;i < t.NumElements();++i) {            
+t.data<double>()[i] = m.c[i];
+        }
         return t;
     }
     std::vector<container::Tensor> mat2ten_double(std::vector<ModuleBase::matrix>& m)
     {
         std::vector<container::Tensor> t;
-        for (int i = 0;i < m.size();++i) { t.push_back(mat2ten_double(m[i]));
-}
+        for (int i = 0;i < m.size();++i) {
+            t.push_back(mat2ten_double(m[i]));
+        }
         return t;
     }
     ModuleBase::matrix ten2mat_double(container::Tensor& t)
     {
         ModuleBase::matrix m(t.shape().dims()[0], t.shape().dims()[1]);
-        for (int i = 0;i < t.NumElements();++i) {m.c[i] = t.data<double>()[i];
-}
+        for (int i = 0;i < t.NumElements();++i) {
+            m.c[i] = t.data<double>()[i];
+        }
         return m;
     }
     std::vector<ModuleBase::matrix> ten2mat_double(std::vector<container::Tensor>& t)
     {
         std::vector<ModuleBase::matrix> m;
-        for (int i = 0;i < t.size();++i) { m.push_back(ten2mat_double(t[i]));
-}
+        for (int i = 0;i < t.size();++i) {
+            m.push_back(ten2mat_double(t[i]));
+        }
         return m;
     }
     container::Tensor mat2ten_complex(ModuleBase::ComplexMatrix& m)
     {
         container::Tensor t(DAT::DT_COMPLEX_DOUBLE, DEV::CpuDevice, { m.nr, m.nc });
-        for (int i = 0;i < t.NumElements();++i) {t.data<std::complex<double>>()[i] = m.c[i];
-}
+        for (int i = 0;i < t.NumElements();++i) {
+            t.data<std::complex<double>>()[i] = m.c[i];
+        }
         return t;
     }
     std::vector<container::Tensor> mat2ten_complex(std::vector<ModuleBase::ComplexMatrix>& m)
     {
         std::vector<container::Tensor> t;
-        for (int i = 0;i < m.size();++i) { t.push_back(mat2ten_complex(m[i]));
-}
+        for (int i = 0;i < m.size();++i) {
+            t.push_back(mat2ten_complex(m[i]));
+        }
         return t;
     }
     ModuleBase::ComplexMatrix ten2mat_complex(container::Tensor& t)
     {
         ModuleBase::ComplexMatrix m(t.shape().dims()[0], t.shape().dims()[1]);
-        for (int i = 0;i < t.NumElements();++i) {m.c[i] = t.data<std::complex<double>>()[i];
-}
+        for (int i = 0;i < t.NumElements();++i) {
+            m.c[i] = t.data<std::complex<double>>()[i];
+        }
         return m;
     }
     std::vector<ModuleBase::ComplexMatrix> ten2mat_complex(std::vector<container::Tensor>& t)
     {
         std::vector<ModuleBase::ComplexMatrix> m;
-        for (int i = 0;i < t.size();++i) { m.push_back(ten2mat_complex(t[i]));
-}
+        for (int i = 0;i < t.size();++i) {
+            m.push_back(ten2mat_complex(t[i]));
+        }
         return m;
     }
 
     ModuleBase::matrix vec2mat(const std::vector<double>& v, const int nr, const int nc)
     {
         assert(v.size() == nr * nc);
         ModuleBase::matrix m(nr, nc, false);
-        for (int i = 0;i < v.size();++i) { m.c[i] = v[i];
-}
+        for (int i = 0;i < v.size();++i) {
+            m.c[i] = v[i];
+        }
         return m;
     }
     ModuleBase::ComplexMatrix vec2mat(const std::vector<std::complex<double>>& v, const int nr, const int nc)
     {
         assert(v.size() == nr * nc);
         ModuleBase::ComplexMatrix m(nr, nc, false);
-        for (int i = 0;i < v.size();++i) { m.c[i] = v[i];
-}
+        for (int i = 0;i < v.size();++i) {
+            m.c[i] = v[i];
+        }
         return m;
     }
     std::vector<ModuleBase::matrix> vec2mat(const std::vector<std::vector<double>>& v, const int nr, const int nc)
     {
         std::vector<ModuleBase::matrix> m(v.size());
-        for (int i = 0;i < v.size();++i) { m[i] = vec2mat(v[i], nr, nc);
-}
+        for (int i = 0;i < v.size();++i) {
+            m[i] = vec2mat(v[i], nr, nc);
+        }
         return m;
     }
     std::vector<ModuleBase::ComplexMatrix> vec2mat(const std::vector<std::vector<std::complex<double>>>& v, const int nr, const int nc)
     {
         std::vector<ModuleBase::ComplexMatrix> m(v.size());
-        for (int i = 0;i < v.size();++i) { m[i] = vec2mat(v[i], nr, nc);
-}
+        for (int i = 0;i < v.size();++i) {
+            m[i] = vec2mat(v[i], nr, nc);
+        }
         return m;
     }
 
@@ -270,32 +286,4 @@ namespace LR_Util
             vl.data(), &ldvl, vr.data(), &ldvr, work2.data(), &lwork, rwork.data(), &info);
         if (info) { std::cout << "ERROR: Lapack solver zgeev, info=" << info << std::endl; }
     }
-    void grad(const double* rhor,
-        ModuleBase::Vector3<double>* gdr,
-        const ModulePW::PW_Basis& rho_basis,
-        const double& tpiba)
-    {
-        std::vector<std::complex<double>> rhog(rho_basis.npw);
-        rho_basis.real2recip(rhor, rhog.data());
-        XC_Functional::grad_rho(rhog.data(), gdr, &rho_basis, tpiba);
-    }
-    void laplace(const double* rhor, double* lapn,
-        const ModulePW::PW_Basis& rho_basis,
-        const double& tpiba2)
-    {
-        ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis.nrxx);
-        std::vector<std::complex<double>> rhog(rho_basis.npw);
-        std::vector<double> 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;
-}
-        }
-    }
-}
+}

source/module_lr/utils/lr_util.h

@@ -36,7 +36,6 @@ namespace LR_Util
     std::pair<ModuleBase::matrix, std::vector<std::pair<int, int>>>
         set_ix_map_diagonal(bool mode, int nc, int nv);
 
-#ifdef  USE_LIBXC
     /// operators to calculate XC kernels
     void grad(const double* rhor,
         ModuleBase::Vector3<double>* gdr,
@@ -46,7 +45,6 @@ namespace LR_Util
         double* lapn,
         const ModulePW::PW_Basis& rho_basis,
         const double& tpiba2);
-#endif
     /// =================ALGORITHM====================
 
     //====== newers and deleters========

source/module_lr/utils/lr_util_xc.cpp

@@ -0,0 +1,28 @@
+#include "lr_util.h"
+namespace LR_Util
+{
+    void grad(const double* rhor,
+        ModuleBase::Vector3<double>* gdr,
+        const ModulePW::PW_Basis& rho_basis,
+        const double& tpiba)
+    {
+        std::vector<std::complex<double>> rhog(rho_basis.npw);
+        rho_basis.real2recip(rhor, rhog.data());
+        XC_Functional::grad_rho(rhog.data(), gdr, &rho_basis, tpiba);
+    }
+    void laplace(const double* rhor, double* lapn,
+        const ModulePW::PW_Basis& rho_basis,
+        const double& tpiba2)
+    {
+        ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis.nrxx);
+        std::vector<std::complex<double>> rhog(rho_basis.npw);
+        std::vector<double> 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; }
+        }
+    }
+}

source/module_lr/utils/test/CMakeLists.txt

@@ -1,15 +1,15 @@
 remove_definitions(-DUSE_LIBXC)
 AddTest(
   TARGET lr_util_phys_test
-  LIBS parameter base ${math_libs} device container
+  LIBS parameter base ${math_libs} device container planewave #for FFT
   SOURCES lr_util_physics_test.cpp  ../lr_util.cpp
   ../../../module_basis/module_ao/parallel_2d.cpp
   ../../../module_io/orb_io.cpp
 )
 
 AddTest(
   TARGET lr_util_algo_test
-  LIBS parameter base ${math_libs} device psi container
+  LIBS parameter base ${math_libs} device psi container planewave #for FFT
   SOURCES lr_util_algorithms_test.cpp ../lr_util.cpp 
   ../../../module_basis/module_ao/parallel_2d.cpp
 )