read fxc or charge for fxc from file

Author: maki49

docs/advanced/input_files/input-main.md

@@ -3921,6 +3921,15 @@ These parameters are used to solve the excited states using. e.g. LR-TDDFT.
 Currently supported: `RPA`, `LDA`, `PBE`, `HSE`, `HF`.
 - **Default**: LDA
 
+### init_fxc
+
+- **Type**: String
+- **Description**:  The method to initalize the xc kernel. 
+  - "gs": Calculate fxc from the ground-state charge density.
+  - "file_fxc": Read the xc kernel $f_\text{xc}$ on grid from the provided files. The following worlds should be the paths of ".cube" files, where the first 3 (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.
+  - "file_chg": 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**: "gs"
+
 ### lr_solver
 
 - **Type**: String

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("init_fxc");
+        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.init_fxc;
+            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.init_fxc.empty()) { para.input.init_fxc.push_back("gs"); }
+            };
+        sync_stringvec(input.init_fxc, para.input.init_fxc.size(), "gs");
+        this->add_item(item);
+    }
     {
         Input_Item item("lr_solver");
         item.annotation = "the eigensolver for LR-TDDFT";

source/module_io/test/read_input_ptest.cpp

@@ -417,6 +417,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.init_fxc[0], "gs");
     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/esolver_lrtd_lcao.cpp

@@ -603,11 +603,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.init_fxc);
         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.init_fxc);
+        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.init_fxc);
         break;
     default:
         throw std::invalid_argument("ESolver_LR: nspin must be 1 or 2");

source/module_lr/potentials/kernel.h

@@ -1,28 +1,27 @@
 #pragma once
 #include "module_basis/module_pw/pw_basis.h"
-#include "module_elecstate/module_charge/charge.h"
-#include "module_cell/unitcell.h"
-// #include <ATen/tensor.h>
 
 namespace LR
 {
     class KernelXC
     {
     public:
-        KernelXC() {};
+        KernelXC(const ModulePW::PW_Basis& rho_basis) : rho_basis_(rho_basis) {};
         ~KernelXC() {};
 
         // xc kernel for LR-TDDFT
 #ifdef USE_LIBXC
-        void f_xc_libxc(const int& nspin, const double& omega, const double& tpiba, const Charge* chg_gs);
+        void f_xc_libxc(const int& nspin, const double& omega, const double& tpiba, const double* const* const rho_gs, const double* const rho_core = nullptr);
 #endif
 
+        void set_kernel(const std::string& name, const std::vector<double>& vec) { this->kernel_set_[name] = vec; }
+        void set_kernel(const std::string& name, const std::vector<double>&& vec) { this->kernel_set_[name] = std::move(vec); }
         const std::vector<double>& get_kernel(const std::string& name) { return kernel_set_[name]; }
         const std::vector<ModuleBase::Vector3<double>>& get_grad_kernel(const std::string& name) { return grad_kernel_set_[name]; }
 
     protected:
 
-        const ModulePW::PW_Basis* rho_basis_ = nullptr;
+        const ModulePW::PW_Basis& rho_basis_;
         std::map<std::string, std::vector<double>> kernel_set_; // [kernel_type][nrxx][nspin]
         std::map<std::string, std::vector<ModuleBase::Vector3<double>>> grad_kernel_set_;// [kernel_type][nrxx][nspin],  intermediate terms for GGA
     };

source/module_lr/potentials/kernel_xc.cpp

@@ -7,44 +7,31 @@
 #include <xc.h>
 #include "module_hamilt_general/module_xc/xc_functional_libxc.h"
 
-void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const double& tpiba, const Charge* chg_gs)
+void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const double& tpiba, const double* const* const rho_gs, const double* const rho_core)
 {
     ModuleBase::TITLE("XC_Functional", "f_xc_libxc");
     ModuleBase::timer::tick("XC_Functional", "f_xc_libxc");
     // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/
 
+    assert(nspin == 1 || nspin == 2);
+
     std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(
         XC_Functional::get_func_id(), 
         (1 == nspin) ? XC_UNPOLARIZED : XC_POLARIZED);
-    int nrxx = chg_gs->nrxx;
+    const int& nrxx = rho_basis_.nrxx;
 
     // converting rho (extract it as a subfuntion in the future)
     // -----------------------------------------------------------------------------------
     std::vector<double> rho(nspin * nrxx);    // r major / spin contigous
-    std::vector<double> amag;
-    if (1 == nspin || 2 == nspin)
-    {
+
 #ifdef _OPENMP
 #pragma omp parallel for collapse(2) schedule(static, 1024)
 #endif
-        for (int is = 0; is < nspin; ++is)
-            for (int ir = 0; ir < nrxx; ++ir)
-                rho[ir * nspin + is] = chg_gs->rho[is][ir] + 1.0 / nspin * chg_gs->rho_core[ir];
-    }
-    else
+    for (int is = 0; is < nspin; ++is) { for (int ir = 0; ir < nrxx; ++ir) { rho[ir * nspin + is] = rho_gs[is][ir]; } }
+    if (rho_core)
     {
-        amag.resize(nrxx);
-#ifdef _OPENMP
-#pragma omp parallel for
-#endif
-        for (int ir = 0; ir < nrxx; ++ir)
-        {
-            const double arhox = std::abs(chg_gs->rho[0][ir] + chg_gs->rho_core[ir]);
-            amag[ir] = std::sqrt(std::pow(chg_gs->rho[1][ir], 2) + std::pow(chg_gs->rho[2][ir], 2) + std::pow(chg_gs->rho[3][ir], 2));
-            const double amag_clip = (amag[ir] < arhox) ? amag[ir] : arhox;
-            rho[ir * nspin + 0] = (arhox + amag_clip) / 2.0;
-            rho[ir * nspin + 1] = (arhox - amag_clip) / 2.0;
-        }
+        const double fac = 1.0 / nspin;
+        for (int is = 0; is < nspin; ++is) { for (int ir = 0; ir < nrxx; ++ir) { rho[ir * nspin + is] += fac * rho_core[ir]; } }
     }
 
     // -----------------------------------------------------------------------------------
@@ -81,7 +68,7 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl
 #endif
                 for (int ir = 0; ir < nrxx; ++ir) rhor[ir] = rho[ir * nspin + is];
                 gdr[is].resize(nrxx);
-                LR_Util::grad(rhor.data(), gdr[is].data(), *(chg_gs->rhopw), tpiba);
+                LR_Util::grad(rhor.data(), gdr[is].data(), rho_basis_, tpiba);
             }
             // 2. |\nabla\rho|^2
             sigma.resize(nrxx * ((1 == nspin) ? 1 : 3));

source/module_lr/potentials/pot_hxc_lrtd.cpp

@@ -6,28 +6,69 @@
 #include "module_hamilt_general/module_xc/xc_functional.h"
 #include <set>
 #include "module_lr/utils/lr_util.h"
-
+#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 }
 namespace LR
 {
     // constructor for exchange-correlation kernel
     PotHxcLR::PotHxcLR(const std::string& xc_kernel_in, const ModulePW::PW_Basis* rho_basis_in, const UnitCell* ucell_in,
-        const Charge* chg_gs/*ground state*/, const SpinType& st_in)
-        :xc_kernel(xc_kernel_in), tpiba_(ucell_in->tpiba), spin_type_(st_in)
+        const Charge* chg_gs/*ground state*/, const Parallel_Grid& pgrid,
+        const SpinType& st_in, const std::vector<std::string>& init_fxc)
+        :xc_kernel(xc_kernel_in), tpiba_(ucell_in->tpiba), spin_type_(st_in),
+        xc_kernel_components_(*rho_basis_in)
     {
         this->rho_basis_ = rho_basis_in;
         this->nrxx = chg_gs->nrxx;
         this->nspin = (PARAM.inp.nspin == 1 || (PARAM.inp.nspin == 4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z)) ? 1 : 2;
 
         this->pot_hartree = LR_Util::make_unique<elecstate::PotHartree>(this->rho_basis_);
+
         std::set<std::string> local_xc = { "lda", "pbe", "hse" };
         if (local_xc.find(this->xc_kernel) != local_xc.end())
         {
             XC_Functional::set_xc_type(this->xc_kernel);    // for hse, (1-alpha) and omega are set here
             this->xc_type_ = XCType(XC_Functional::get_func_type());
-#ifdef USE_LIBXC
             this->set_integral_func(this->spin_type_, this->xc_type_);
-            this->xc_kernel_components_.f_xc_libxc(nspin, ucell_in->omega, ucell_in->tpiba, chg_gs);
+
+            if (init_fxc[0] == "file_fxc")
+            {
+                assert(this->xc_kernel == "lda");
+                const int spinsize = (1 == nspin) ? 1 : 3;
+                std::vector<double> v2rho2(spinsize * nrxx);
+                // read fxc adn add to xc_kernel_components
+                assert(init_fxc.size() >= spinsize + 1);
+                for (int is = 0;is < spinsize;++is)
+                {
+                    double ef = 0.0;
+                    int prenspin = 1;
+                    std::vector<double> v2rho2_tmp(nrxx);
+                    ModuleIO::read_vdata_palgrid(GlobalC::Pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, init_fxc[is + 1],
+                        v2rho2_tmp.data(), ucell_in->nat);
+                    for (int ir = 0;ir < nrxx;++ir) { v2rho2[ir * spinsize + is] = v2rho2_tmp[ir]; }
+                }
+                this->xc_kernel_components_.set_kernel("v2rho2", std::move(v2rho2));
+                return;
+            }
+
+#ifdef USE_LIBXC
+            if (init_fxc[0] == "file_chg")
+            {
+                assert(init_fxc.size() >= 2);
+                double** rho_for_fxc;
+                LR_Util::new_p2(rho_for_fxc, nspin, nrxx);
+                double ef = 0.0;
+                int prenspin = 1;
+                for (int is = 0;is < nspin;++is)
+                {
+                    const std::string file = init_fxc[init_fxc.size() > nspin ? 1 + is : 1];
+                    ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, file,
+                        rho_for_fxc[is], ucell_in->nat);
+                }
+                this->xc_kernel_components_.f_xc_libxc(nspin, ucell_in->omega, ucell_in->tpiba, rho_for_fxc, chg_gs->rho_core);
+                LR_Util::delete_p2(rho_for_fxc, nspin);
+            }
+            else { this->xc_kernel_components_.f_xc_libxc(nspin, ucell_in->omega, ucell_in->tpiba, chg_gs->rho, chg_gs->rho_core); }
 #else 
             ModuleBase::WARNING_QUIT("KernelXC", "to calculate xc-kernel in LR-TDDFT, compile with LIBXC");
 #endif
@@ -63,7 +104,6 @@ namespace LR
         ModuleBase::timer::tick("PotHxcLR", "cal_v_eff");
     }
 
-#ifdef USE_LIBXC
     void PotHxcLR::set_integral_func(const SpinType& s, const XCType& xc)
     {
         auto& funcs = this->kernel_to_potential_;
@@ -168,5 +208,4 @@ namespace LR
                 + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__));
         }
     }
-#endif
 } // namespace LR

source/module_lr/potentials/pot_hxc_lrtd.h

@@ -4,6 +4,8 @@
 #include "kernel.h"
 #include <unordered_map>
 #include <memory>
+
+class Parallel_Grid;
 namespace LR
 {
     class PotHxcLR : public elecstate::PotBase
@@ -15,8 +17,9 @@ namespace LR
         enum SpinType { S1 = 0, S2_singlet = 1, S2_triplet = 2, S2_updown = 3 };
         enum XCType { None = 0, LDA = 1, GGA = 2, HYB_GGA = 4 };
         /// constructor for exchange-correlation kernel
-        PotHxcLR(const std::string& xc_kernel_in, const ModulePW::PW_Basis* rho_basis_in, const UnitCell* ucell_in, const Charge* chg_gs/*ground state*/,
-            const SpinType& st_in = SpinType::S1);
+        PotHxcLR(const std::string& xc_kernel_in, const ModulePW::PW_Basis* rho_basis_in,
+            const UnitCell* ucell_in, const Charge* chg_gs/*ground state*/, const Parallel_Grid& pgrid,
+            const SpinType& st_in = SpinType::S1, const std::vector<std::string>& init_fxc = { "gs" });
         ~PotHxcLR() {}
         void cal_v_eff(const Charge* chg/*excited state*/, const UnitCell* ucell, ModuleBase::matrix& v_eff) override {};
         void cal_v_eff(double** rho, const UnitCell* ucell, ModuleBase::matrix& v_eff, const std::vector<int>& ispin_op = { 0,0 });

source/module_lr/utils/lr_util.cpp

@@ -270,7 +270,6 @@ 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; }
     }
-#ifdef USE_LIBXC
     void grad(const double* rhor,
         ModuleBase::Vector3<double>* gdr,
         const ModulePW::PW_Basis& rho_basis,
@@ -299,5 +298,4 @@ namespace LR_Util
 }
         }
     }
-#endif
 }

source/module_parameter/input_parameter.h

@@ -296,6 +296,7 @@ struct Input_para
 
     // ==============   #Parameters (10.lr-tddft) ===========================
     int lr_nstates = 1; ///< the number of 2-particle states to be solved
+    std::vector<std::string> init_fxc = {};    ///< The method to initalize the xc kernel
     int nocc = -1;      ///< the number of occupied orbitals to form the 2-particle basis
     int nvirt = 1;      ///< the number of virtual orbitals to form the 2-particle basis (nocc + nvirt <= nbands)
     std::string xc_kernel = "LDA"; ///< exchange correlation (XC) kernel for LR-TDDFT