Feature: LR-TDDFT support reading fxc from file or calculating fxc by a specified charge file

docs/advanced/input_files/input-main.md

@@ -423,6 +423,7 @@
     - [pexsi\_zero\_thr](#pexsi_zero_thr)
   - [Linear Response TDDFT](#linear-response-tddft)
     - [xc\_kernel](#xc_kernel)
+    - [lr\_init\_xc\_kernel](#lr_init_xc_kernel)
     - [lr\_solver](#lr_solver)
     - [lr\_thr](#lr_thr)
     - [nocc](#nocc)
@@ -3943,6 +3944,15 @@ These parameters are used to solve the excited states using. e.g. LR-TDDFT.
 Currently supported: `RPA`, `LDA`, `PBE`, `HSE`, `HF`.
 - **Default**: LDA
 
+### lr_init_xc_kernel
+
+- **Type**: String
+- **Description**: The method to initalize the xc kernel. 
+  - "default": Calculate xc kerenel ($f_\text{xc}$) from the ground-state charge density.
+  - "file": Read the xc kernel $f_\text{xc}$ on grid from the provided files. The following words should be the paths of ".cube" files, where the first 1 (*[nspin](#nspin)==1*) or 3 (*[nspin](#nspin)==2*, namely spin-aa, spin-ab and spin-bb) will be read in. The parameter [xc_kernel](#xc_kernel) will be invalid. Now only LDA-type kernel is supproted as the potential will be calculated by directly multiplying the transition density.
+  - "from_charge_file": Calculate fxc from the charge density read from the provided files. The following words should be the paths of ".cube" files, where the first [nspin]($nspin) files will be read in. 
+- **Default**: "default"
+
 ### lr_solver
 
 - **Type**: String

source/Makefile.Objects

@@ -731,7 +731,7 @@ OBJS_TENSOR=tensor.o\
     dmr_complex.o\
     operator_lr_hxc.o\
     operator_lr_exx.o\
-    kernel_xc.o\
+    xc_kernel.o\
     pot_hxc_lrtd.o\
     lr_spectrum.o\
     hamilt_casida.o\

source/module_io/read_input_item_tddft.cpp

@@ -309,6 +309,20 @@ void ReadInput::item_lr_tddft()
         read_sync_string(input.xc_kernel);
         this->add_item(item);
     }
+    {
+        Input_Item item("lr_init_xc_kernel");
+        item.annotation = "The method to initalize the xc kernel";
+        item.read_value = [](const Input_Item& item, Parameter& para) {
+            size_t count = item.get_size();
+            auto& ifxc = para.input.lr_init_xc_kernel;
+            for (int i = 0; i < count; i++) { ifxc.push_back(item.str_values[i]); }
+            };
+        item.reset_value = [](const Input_Item& item, Parameter& para) {
+            if (para.input.lr_init_xc_kernel.empty()) { para.input.lr_init_xc_kernel.push_back("default"); }
+            };
+        sync_stringvec(input.lr_init_xc_kernel, para.input.lr_init_xc_kernel.size(), "default");
+        this->add_item(item);
+    }
     {
         Input_Item item("lr_solver");
         item.annotation = "the eigensolver for LR-TDDFT";

source/module_io/test/read_input_ptest.cpp

@@ -420,6 +420,7 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_EQ(param.inp.nocc, param.inp.nbands);
     EXPECT_EQ(param.inp.nvirt, 1);
     EXPECT_EQ(param.inp.xc_kernel, "LDA");
+    EXPECT_EQ(param.inp.lr_init_xc_kernel[0], "default");
     EXPECT_EQ(param.inp.lr_solver, "dav");
     EXPECT_DOUBLE_EQ(param.inp.lr_thr, 1e-2);
     EXPECT_FALSE(param.inp.lr_unrestricted);

source/module_lr/CMakeLists.txt

@@ -17,13 +17,8 @@ if(ENABLE_LCAO)
     potentials/pot_hxc_lrtd.cpp
     lr_spectrum.cpp
     hamilt_casida.cpp
-    esolver_lrtd_lcao.cpp)
-
-    if(ENABLE_LIBXC)
-      list(APPEND objects
-        potentials/kernel_xc.cpp
-      )
-    endif()
+    esolver_lrtd_lcao.cpp
+    potentials/xc_kernel.cpp)
 
 add_library(
     lr

source/module_lr/esolver_lrtd_lcao.cpp

@@ -65,8 +65,8 @@ inline int cal_nupdown_form_occ(const ModuleBase::matrix& wg)
 template<typename T, typename TR>
 void LR::ESolver_LR<T, TR>::parameter_check()const
 {
-    std::set<std::string> lr_solvers = { "dav", "lapack" , "spectrum", "dav_subspace", "cg" };
-    std::set<std::string> xc_kernels = { "rpa", "lda", "pbe", "hf" , "hse" };
+    const std::set<std::string> lr_solvers = { "dav", "lapack" , "spectrum", "dav_subspace", "cg" };
+    const std::set<std::string> xc_kernels = { "rpa", "lda", "pwlda", "pbe", "hf" , "hse" };
     if (lr_solvers.find(this->input.lr_solver) == lr_solvers.end()) {
         throw std::invalid_argument("ESolver_LR: unknown type of lr_solver");
 }
@@ -104,7 +104,13 @@ void LR::ESolver_LR<T, TR>::set_dimension()
     GlobalV::ofs_running << "number of excited states to be solved: " << this->nstates << std::endl;
     if (input.ri_hartree_benchmark == "aims" && !input.aims_nbasis.empty())
     {
-        this->nbasis = [&]() -> int { int nbas = 0; for (int it = 0;it < ucell.ntype;++it) { nbas += ucell.atoms[it].na * input.aims_nbasis[it]; };return nbas;}();
+        // calculate total number of basis funcs, see https://en.cppreference.com/w/cpp/algorithm/inner_product
+        this->nbasis = std::inner_product(input.aims_nbasis.begin(), /* iterator1.begin */
+                                          input.aims_nbasis.end(),  /* iterator1.end */
+                                          ucell.atoms,  /* iterator2.begin */
+                                          0,  /* init value */
+                                          std::plus<int>(), /* iter op1 */
+                                          [](const int& a, const Atom& b) { return a * b.na; }); /* iter op2 */
         std::cout << "nbasis from aims: " << this->nbasis << std::endl;
     }
 }
@@ -592,11 +598,11 @@ void LR::ESolver_LR<T, TR>::init_pot(const Charge& chg_gs)
     {
         using ST = PotHxcLR::SpinType;
     case 1:
-        this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, ST::S1);
+        this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, ucell, chg_gs, GlobalC::Pgrid, ST::S1, input.lr_init_xc_kernel);
         break;
     case 2:
-        this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, openshell ? ST::S2_updown : ST::S2_singlet);
-        this->pot[1] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, openshell ? ST::S2_updown : ST::S2_triplet);
+        this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, ucell, chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_singlet, input.lr_init_xc_kernel);
+        this->pot[1] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, ucell, chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_triplet, input.lr_init_xc_kernel);
         break;
     default:
         throw std::invalid_argument("ESolver_LR: nspin must be 1 or 2");

source/module_lr/lr_spectrum.cpp

@@ -53,7 +53,7 @@ void LR::LR_Spectrum<double>::oscillator_strength(Grid_Driver& gd, const std::ve
 
             // 2. transition density
             double** rho_trans;
-            LR_Util::new_p2(rho_trans, 1, this->rho_basis.nrxx);
+            LR_Util::_allocate_2order_nested_ptr(rho_trans, 1, this->rho_basis.nrxx);
             this->cal_gint_rho(rho_trans, this->rho_basis.nrxx);
 
             // 3. transition dipole moment
@@ -67,7 +67,7 @@ void LR::LR_Spectrum<double>::oscillator_strength(Grid_Driver& gd, const std::ve
                 ModuleBase::Vector3<double> rc = rd * ucell.latvec * ucell.lat0; // real coordinate
                 transition_dipole_[istate] += rc * rho_trans[0][ir];
             }
-            LR_Util::delete_p2(rho_trans, 1);
+            LR_Util::_deallocate_2order_nested_ptr(rho_trans, 1);
         }
         transition_dipole_[istate] *= (ucell.omega / static_cast<double>(gint->get_ncxyz()));   // dv
         Parallel_Reduce::reduce_all(transition_dipole_[istate].x);
@@ -114,8 +114,8 @@ void LR::LR_Spectrum<std::complex<double>>::oscillator_strength(Grid_Driver& gd,
             // 2. transition density
             double** rho_trans_real;
             double** rho_trans_imag;
-            LR_Util::new_p2(rho_trans_real, 1, this->rho_basis.nrxx);
-            LR_Util::new_p2(rho_trans_imag, 1, this->rho_basis.nrxx);
+            LR_Util::_allocate_2order_nested_ptr(rho_trans_real, 1, this->rho_basis.nrxx);
+            LR_Util::_allocate_2order_nested_ptr(rho_trans_imag, 1, this->rho_basis.nrxx);
             // real part
             LR_Util::get_DMR_real_imag_part(DM_trans, DM_trans_real_imag, ucell.nat, 'R');
             this->gint->transfer_DM2DtoGrid(DM_trans_real_imag.get_DMR_vector());
@@ -140,8 +140,8 @@ void LR::LR_Spectrum<std::complex<double>>::oscillator_strength(Grid_Driver& gd,
                 ModuleBase::Vector3<std::complex<double>> rc_complex(rc.x, rc.y, rc.z);
                 transition_dipole_[istate] += rc_complex * std::complex<double>(rho_trans_real[0][ir], rho_trans_imag[0][ir]);
             }
-            LR_Util::delete_p2(rho_trans_real, 1);
-            LR_Util::delete_p2(rho_trans_imag, 1);
+            LR_Util::_deallocate_2order_nested_ptr(rho_trans_real, 1);
+            LR_Util::_deallocate_2order_nested_ptr(rho_trans_imag, 1);
         }
         transition_dipole_[istate] *= (ucell.omega / static_cast<double>(gint->get_ncxyz()));   // dv
         Parallel_Reduce::reduce_all(transition_dipole_[istate].x);

source/module_lr/operator_casida/operator_lr_hxc.cpp

@@ -61,15 +61,15 @@ namespace LR
         // \f[ \tilde{\rho}(r)=\sum_{\mu_j, \mu_b}\tilde{\rho}_{\mu_j,\mu_b}\phi_{\mu_b}(r)\phi_{\mu_j}(r) \f]
         double** rho_trans;
         const int& nrxx = this->pot.lock()->nrxx;
-        LR_Util::new_p2(rho_trans, 1, nrxx); // currently gint_kernel_rho uses PARAM.inp.nspin, it needs refactor
+        LR_Util::_allocate_2order_nested_ptr(rho_trans, 1, nrxx); // currently gint_kernel_rho uses PARAM.inp.nspin, it needs refactor
         ModuleBase::GlobalFunc::ZEROS(rho_trans[0], nrxx);
         Gint_inout inout_rho(rho_trans, Gint_Tools::job_type::rho, 1, false);
         this->gint->cal_gint(&inout_rho);
 
         // 3. v_hxc = f_hxc * rho_trans
         ModuleBase::matrix vr_hxc(1, nrxx);   //grid
-        this->pot.lock()->cal_v_eff(rho_trans, &GlobalC::ucell, vr_hxc, ispin_ks);
-        LR_Util::delete_p2(rho_trans, 1);
+        this->pot.lock()->cal_v_eff(rho_trans, GlobalC::ucell, vr_hxc, ispin_ks);
+        LR_Util::_deallocate_2order_nested_ptr(rho_trans, 1);
 
         // 4. V^{Hxc}_{\mu,\nu}=\int{dr} \phi_\mu(r) v_{Hxc}(r) \phi_\mu(r)
         Gint_inout inout_vlocal(vr_hxc.c, 0, Gint_Tools::job_type::vlocal);
@@ -103,18 +103,18 @@ namespace LR
                 double** rho_trans;
                 const int& nrxx = this->pot.lock()->nrxx;
 
-                LR_Util::new_p2(rho_trans, 1, nrxx); // nspin=1 for transition density
+                LR_Util::_allocate_2order_nested_ptr(rho_trans, 1, nrxx); // nspin=1 for transition density
                 ModuleBase::GlobalFunc::ZEROS(rho_trans[0], nrxx);
                 Gint_inout inout_rho(rho_trans, Gint_Tools::job_type::rho, 1, false);
                 this->gint->cal_gint(&inout_rho);
                 // print_grid_nonzero(rho_trans[0], nrxx, 10, "rho_trans");
 
                 // 3. v_hxc = f_hxc * rho_trans
                 ModuleBase::matrix vr_hxc(1, nrxx);   //grid
-                this->pot.lock()->cal_v_eff(rho_trans, &GlobalC::ucell, vr_hxc, ispin_ks);
+                this->pot.lock()->cal_v_eff(rho_trans, GlobalC::ucell, vr_hxc, ispin_ks);
                 // print_grid_nonzero(vr_hxc.c, this->poticab->nrxx, 10, "vr_hxc");
 
-                LR_Util::delete_p2(rho_trans, 1);
+                LR_Util::_deallocate_2order_nested_ptr(rho_trans, 1);
 
                 // 4. V^{Hxc}_{\mu,\nu}=\int{dr} \phi_\mu(r) v_{Hxc}(r) \phi_\mu(r)
                 Gint_inout inout_vlocal(vr_hxc.c, 0, Gint_Tools::job_type::vlocal);