From b364da360615bc750174c8d852a2e578754665db Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Wed, 23 Oct 2024 01:18:39 +0800 Subject: [PATCH 01/23] read fxc or charge for fxc from file --- docs/advanced/input_files/input-main.md | 9 ++++ source/module_io/read_input_item_tddft.cpp | 14 ++++++ source/module_io/test/read_input_ptest.cpp | 1 + source/module_lr/esolver_lrtd_lcao.cpp | 6 +-- source/module_lr/potentials/kernel.h | 11 ++-- source/module_lr/potentials/kernel_xc.cpp | 33 ++++-------- source/module_lr/potentials/pot_hxc_lrtd.cpp | 53 +++++++++++++++++--- source/module_lr/potentials/pot_hxc_lrtd.h | 7 ++- source/module_lr/utils/lr_util.cpp | 2 - source/module_parameter/input_parameter.h | 1 + 10 files changed, 94 insertions(+), 43 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 40edbf5538..01c7ca34a8 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -3943,6 +3943,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 diff --git a/source/module_io/read_input_item_tddft.cpp b/source/module_io/read_input_item_tddft.cpp index ccd3190c54..2d6bc6c265 100644 --- a/source/module_io/read_input_item_tddft.cpp +++ b/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"; diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index 69d6d13a1c..07a34485c0 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/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.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); diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index aca88c38ae..e956e51bae 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -592,11 +592,11 @@ void LR::ESolver_LR::init_pot(const Charge& chg_gs) { using ST = PotHxcLR::SpinType; case 1: - this->pot[0] = std::make_shared(xc_kernel, this->pw_rho, &ucell, &chg_gs, ST::S1); + this->pot[0] = std::make_shared(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, ST::S1, input.init_fxc); break; case 2: - this->pot[0] = std::make_shared(xc_kernel, this->pw_rho, &ucell, &chg_gs, openshell ? ST::S2_updown : ST::S2_singlet); - this->pot[1] = std::make_shared(xc_kernel, this->pw_rho, &ucell, &chg_gs, openshell ? ST::S2_updown : ST::S2_triplet); + this->pot[0] = std::make_shared(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(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"); diff --git a/source/module_lr/potentials/kernel.h b/source/module_lr/potentials/kernel.h index 67ab8d0868..38da6087a4 100644 --- a/source/module_lr/potentials/kernel.h +++ b/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 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& vec) { this->kernel_set_[name] = vec; } + void set_kernel(const std::string& name, const std::vector&& vec) { this->kernel_set_[name] = std::move(vec); } const std::vector& get_kernel(const std::string& name) { return kernel_set_[name]; } const std::vector>& 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> kernel_set_; // [kernel_type][nrxx][nspin] std::map>> grad_kernel_set_;// [kernel_type][nrxx][nspin], intermediate terms for GGA }; diff --git a/source/module_lr/potentials/kernel_xc.cpp b/source/module_lr/potentials/kernel_xc.cpp index dee51ef3c5..96b0aa944f 100644 --- a/source/module_lr/potentials/kernel_xc.cpp +++ b/source/module_lr/potentials/kernel_xc.cpp @@ -7,44 +7,31 @@ #include #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 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 rho(nspin * nrxx); // r major / spin contigous - std::vector 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)); diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index ebdb8747c5..461ae97a76 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -6,28 +6,69 @@ #include "module_hamilt_general/module_xc/xc_functional.h" #include #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& 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& 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(this->rho_basis_); + std::set 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 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 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 diff --git a/source/module_lr/potentials/pot_hxc_lrtd.h b/source/module_lr/potentials/pot_hxc_lrtd.h index a04ab42bb9..f688c41aef 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.h +++ b/source/module_lr/potentials/pot_hxc_lrtd.h @@ -4,6 +4,8 @@ #include "kernel.h" #include #include + +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& 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& ispin_op = { 0,0 }); diff --git a/source/module_lr/utils/lr_util.cpp b/source/module_lr/utils/lr_util.cpp index eb5ec2bc43..fdb57a1e81 100644 --- a/source/module_lr/utils/lr_util.cpp +++ b/source/module_lr/utils/lr_util.cpp @@ -266,7 +266,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* gdr, const ModulePW::PW_Basis& rho_basis, @@ -295,5 +294,4 @@ namespace LR_Util } } } -#endif } diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index bfd9ae91cb..4e62414a03 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -299,6 +299,7 @@ struct Input_para // ============== #Parameters (10.lr-tddft) =========================== int lr_nstates = 1; ///< the number of 2-particle states to be solved + std::vector 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 From 1572d0b8ede10aa947e282a202b4c74aa58a4da3 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Sat, 2 Nov 2024 04:44:36 +0800 Subject: [PATCH 02/23] move xc-dependent functions to a new file --- source/Makefile.Objects | 1 + source/module_lr/CMakeLists.txt | 1 + source/module_lr/utils/lr_util.cpp | 101 +++++++-------------- source/module_lr/utils/lr_util.h | 2 - source/module_lr/utils/lr_util_xc.cpp | 28 ++++++ source/module_lr/utils/test/CMakeLists.txt | 4 +- 6 files changed, 67 insertions(+), 70 deletions(-) create mode 100644 source/module_lr/utils/lr_util_xc.cpp diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 351c8c340a..fd087180ea 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -724,6 +724,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\ diff --git a/source/module_lr/CMakeLists.txt b/source/module_lr/CMakeLists.txt index 4b816f09b3..f609c28b00 100644 --- a/source/module_lr/CMakeLists.txt +++ b/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 diff --git a/source/module_lr/utils/lr_util.cpp b/source/module_lr/utils/lr_util.cpp index fdb57a1e81..4ae65fcaea 100644 --- a/source/module_lr/utils/lr_util.cpp +++ b/source/module_lr/utils/lr_util.cpp @@ -80,74 +80,67 @@ namespace LR_Util void matsym>(const std::complex* in, const int n, const Parallel_2D& pmat, std::complex* out) { for (int i = 0;i < pmat.get_local_size();++i) { out[i] = in[i]; } - const std::complex 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); - } - template<> - void matsym>(std::complex* inout, const int n, const Parallel_2D& pmat) - { std::vector> tmp(pmat.get_local_size()); std::copy(inout, inout + pmat.get_local_size(), tmp.begin()); - const std::complex 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); - } -#endif - 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()[i] = m.c[i]; -} + for (int i = 0;i < t.NumElements();++i) { t.data()[i] = m.c[i]; } return t; } std::vector mat2ten_double(std::vector& m) { std::vector 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()[i]; -} + for (int i = 0;i < t.NumElements();++i) { + m.c[i] = t.data()[i]; + } return m; } std::vector ten2mat_double(std::vector& t) { std::vector 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>()[i] = m.c[i]; -} + for (int i = 0;i < t.NumElements();++i) { + t.data>()[i] = m.c[i]; + } return t; } std::vector mat2ten_complex(std::vector& m) { std::vector 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>()[i]; -} + for (int i = 0;i < t.NumElements();++i) { + m.c[i] = t.data>()[i]; + } return m; } std::vector ten2mat_complex(std::vector& t) { std::vector 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; } @@ -155,30 +148,34 @@ namespace LR_Util { 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>& 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 vec2mat(const std::vector>& v, const int nr, const int nc) { std::vector 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 vec2mat(const std::vector>>& v, const int nr, const int nc) { std::vector 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; } @@ -266,32 +263,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* gdr, - const ModulePW::PW_Basis& rho_basis, - const double& tpiba) - { - std::vector> 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> rhog(rho_basis.npw); - std::vector 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; -} - } - } -} +} \ No newline at end of file diff --git a/source/module_lr/utils/lr_util.h b/source/module_lr/utils/lr_util.h index 93f26be2a6..f6f70014db 100644 --- a/source/module_lr/utils/lr_util.h +++ b/source/module_lr/utils/lr_util.h @@ -36,7 +36,6 @@ namespace LR_Util std::pair>> 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* gdr, @@ -46,7 +45,6 @@ namespace LR_Util double* lapn, const ModulePW::PW_Basis& rho_basis, const double& tpiba2); -#endif /// =================ALGORITHM==================== //====== newers and deleters======== diff --git a/source/module_lr/utils/lr_util_xc.cpp b/source/module_lr/utils/lr_util_xc.cpp new file mode 100644 index 0000000000..cc8b3a1f6d --- /dev/null +++ b/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* gdr, + const ModulePW::PW_Basis& rho_basis, + const double& tpiba) + { + std::vector> 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> rhog(rho_basis.npw); + std::vector 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; } + } + } +} diff --git a/source/module_lr/utils/test/CMakeLists.txt b/source/module_lr/utils/test/CMakeLists.txt index 53d1db607e..fde01fa69b 100644 --- a/source/module_lr/utils/test/CMakeLists.txt +++ b/source/module_lr/utils/test/CMakeLists.txt @@ -1,7 +1,7 @@ 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 @@ -9,7 +9,7 @@ AddTest( 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 ) \ No newline at end of file From c027abb8dc9e9ff9ac19ef261658dcb37a1dc877 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Sat, 2 Nov 2024 13:54:51 +0800 Subject: [PATCH 03/23] minor fix --- docs/advanced/input_files/input-main.md | 6 +++--- source/module_lr/potentials/pot_hxc_lrtd.cpp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 01c7ca34a8..65e37d2e54 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -3946,9 +3946,9 @@ Currently supported: `RPA`, `LDA`, `PBE`, `HSE`, `HF`. ### 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. +- **Description**: The method to initalize the xc kernel. + - "gs": Calculate xc kerenel ($f_\text{xc}$) from the ground-state charge density. + - "file_fxc": 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 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" diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index 461ae97a76..baf68d185c 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -43,7 +43,7 @@ namespace LR double ef = 0.0; int prenspin = 1; std::vector v2rho2_tmp(nrxx); - ModuleIO::read_vdata_palgrid(GlobalC::Pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, init_fxc[is + 1], + ModuleIO::read_vdata_palgrid(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]; } } From d4f19862195e25bf6d8e34f85a66c1d52b7d5d91 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Sat, 2 Nov 2024 13:58:04 +0800 Subject: [PATCH 04/23] rename files --- source/Makefile.Objects | 2 +- source/module_lr/CMakeLists.txt | 2 +- source/module_lr/potentials/kernel_xc.cpp | 221 ------------------ source/module_lr/potentials/pot_hxc_lrtd.h | 2 +- source/module_lr/potentials/xc_kernel.cpp | 221 ++++++++++++++++++ .../potentials/{kernel.h => xc_kernel.h} | 0 6 files changed, 224 insertions(+), 224 deletions(-) delete mode 100644 source/module_lr/potentials/kernel_xc.cpp create mode 100644 source/module_lr/potentials/xc_kernel.cpp rename source/module_lr/potentials/{kernel.h => xc_kernel.h} (100%) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index fd087180ea..89fc9effd7 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -732,7 +732,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\ diff --git a/source/module_lr/CMakeLists.txt b/source/module_lr/CMakeLists.txt index f609c28b00..93b4f7ff9f 100644 --- a/source/module_lr/CMakeLists.txt +++ b/source/module_lr/CMakeLists.txt @@ -22,7 +22,7 @@ if(ENABLE_LCAO) if(ENABLE_LIBXC) list(APPEND objects - potentials/kernel_xc.cpp + potentials/xc_kernel.cpp ) endif() diff --git a/source/module_lr/potentials/kernel_xc.cpp b/source/module_lr/potentials/kernel_xc.cpp deleted file mode 100644 index 96b0aa944f..0000000000 --- a/source/module_lr/potentials/kernel_xc.cpp +++ /dev/null @@ -1,221 +0,0 @@ -#include "kernel.h" -#include "module_hamilt_general/module_xc/xc_functional.h" -#include "module_parameter/parameter.h" -#include "module_base/timer.h" -#include "module_lr/utils/lr_util.h" -#ifdef USE_LIBXC -#include -#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 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 funcs = XC_Functional_Libxc::init_func( - XC_Functional::get_func_id(), - (1 == nspin) ? XC_UNPOLARIZED : XC_POLARIZED); - const int& nrxx = rho_basis_.nrxx; - - // converting rho (extract it as a subfuntion in the future) - // ----------------------------------------------------------------------------------- - std::vector rho(nspin * nrxx); // r major / spin contigous - -#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] = rho_gs[is][ir]; } } - if (rho_core) - { - 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]; } } - } - - // ----------------------------------------------------------------------------------- - // for GGA - const bool is_gga = [&funcs]() - { - for (xc_func_type& func : funcs) - { - switch (func.info->family) - { - case XC_FAMILY_GGA: - case XC_FAMILY_HYB_GGA: - return true; - } - } - return false; - }(); - std::vector>> gdr; // \nabla \rho - std::vector sigma; // |\nabla\rho|^2 - std::vector sgn; // sgn for threshold mask - if (is_gga) - { - // 0. set up sgn for threshold mask - // in the case of GGA correlation for polarized case, - // a cutoff for grho is required to ensure that libxc gives reasonable results - - // 1. \nabla \rho - gdr.resize(nspin); - for (int is = 0; is < nspin; ++is) - { - std::vector rhor(nrxx); -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 1024) -#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(), rho_basis_, tpiba); - } - // 2. |\nabla\rho|^2 - sigma.resize(nrxx * ((1 == nspin) ? 1 : 3)); - if (1 == nspin) - { -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 1024) -#endif - for (int ir = 0; ir < nrxx; ++ir) - sigma[ir] = gdr[0][ir] * gdr[0][ir]; - } - else - { -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 256) -#endif - for (int ir = 0; ir < nrxx; ++ir) - { - sigma[ir * 3] = gdr[0][ir] * gdr[0][ir]; - sigma[ir * 3 + 1] = gdr[0][ir] * gdr[1][ir]; - sigma[ir * 3 + 2] = gdr[1][ir] * gdr[1][ir]; - } - } - } - // ----------------------------------------------------------------------------------- - //==================== XC Kernels (f_xc)============================= - this->kernel_set_.emplace("vrho", std::vector(nspin * nrxx)); - this->kernel_set_.emplace("v2rho2", std::vector(((1 == nspin) ? 1 : 3) * nrxx));//(nrxx* ((1 == nspin) ? 1 : 3)): 00, 01, 11 - if (is_gga) - { - this->kernel_set_.emplace("vsigma", std::vector(((1 == nspin) ? 1 : 3) * nrxx)); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 - this->kernel_set_.emplace("v2rhosigma", std::vector(((1 == nspin) ? 1 : 6) * nrxx)); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 - this->kernel_set_.emplace("v2sigma2", std::vector(((1 == nspin) ? 1 : 6) * nrxx)); //(nrxx* ((1 == nspin) ? 1 : 6)): 00, 01, 02, 11, 12, 22 - } - //MetaGGA ... - - for (xc_func_type& func : funcs) - { - const double rho_threshold = 1E-6; - const double grho_threshold = 1E-10; - - xc_func_set_dens_threshold(&func, rho_threshold); - - //cut off grho if not LDA (future subfunc) - - switch (func.info->family) - { - case XC_FAMILY_LDA: - xc_lda_vxc(&func, nrxx, rho.data(), this->kernel_set_["vrho"].data()); - xc_lda_fxc(&func, nrxx, rho.data(), this->kernel_set_["v2rho2"].data()); - break; - case XC_FAMILY_GGA: - case XC_FAMILY_HYB_GGA: - xc_gga_vxc(&func, nrxx, rho.data(), sigma.data(), - this->kernel_set_["vrho"].data(), this->kernel_set_["vsigma"].data()); - xc_gga_fxc(&func, nrxx, rho.data(), sigma.data(), - this->kernel_set_["v2rho2"].data(), this->kernel_set_["v2rhosigma"].data(), this->kernel_set_["v2sigma2"].data()); - break; - default: - throw std::domain_error("func.info->family =" + std::to_string(func.info->family) - + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); - break; - } - // some formulas for GGA - if (func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) - { - const std::vector& v2r2 = this->kernel_set_["v2rho2"]; - const std::vector& v2rs = this->kernel_set_["v2rhosigma"]; - const std::vector& v2s2 = this->kernel_set_["v2sigma2"]; - const std::vector& vs = this->kernel_set_["vsigma"]; - const double tpiba2 = tpiba * tpiba; - - if (1 == nspin) - { - // 0. drho - this->grad_kernel_set_.emplace("drho_gs", std::vector>(nrxx)); - for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["drho_gs"][ir] = gdr[0][ir]; - // 1. $2f^{\rho\sigma}*\nabla\rho$ - this->grad_kernel_set_.emplace("2_v2rhosigma_drho", std::vector>(nrxx)); - for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["2_v2rhosigma_drho"][ir] = gdr[0][ir] * v2rs.at(ir) * 2.; - // 2. $4f^{\sigma\sigma}*\nabla\rho$ - this->grad_kernel_set_.emplace("4_v2sigma2_drho", std::vector>(nrxx)); - for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["4_v2sigma2_drho"][ir] = sigma.at(ir) * v2s2.at(ir) * 4.; - } - // else if (2 == nspin) // wrong, to be fixed - // { - // // 1. $\nabla\cdot(f^{\rho\sigma}*\nabla\rho)$ - // std::vector div_v2rhosigma_gdrho_r(3 * nrxx); - // std::vector> v2rhosigma_gdrho_r(3 * nrxx); - // for (int ir = 0; ir < nrxx; ++ir) - // { - // v2rhosigma_gdrho_r[ir] = gdr[0][ir] * v2rs.at(ir * 6) * 4.0 - // + gdr[1][ir] * v2rs.at(ir * 6 + 1) * 2.0; //up-up - // v2rhosigma_gdrho_r[nrxx + ir] = gdr[0][ir] * (v2rs.at(ir * 6 + 3) * 2.0 + v2rs.at(ir * 6 + 1)) - // + gdr[1][ir] * (v2rs.at(ir * 6 + 2) * 2.0 + v2rs.at(ir * 6 + 4)); //up-down - // v2rhosigma_gdrho_r[2 * nrxx + ir] = gdr[1][ir] * v2rs.at(ir * 6 + 5) * 4.0 - // + gdr[0][ir] * v2rs.at(ir * 6 + 4) * 2.0; //down-down - // } - // for (int isig = 0;isig < 3;++isig) - // XC_Functional::grad_dot(v2rhosigma_gdrho_r.data() + isig * nrxx, div_v2rhosigma_gdrho_r.data() + isig * nrxx, chg_gs->rhopw, tpiba); - // // 2. $\nabla^2(f^{\sigma\sigma}*\sigma)$ - // std::vector v2sigma2_sigma_r(3 * nrxx); - // for (int ir = 0; ir < nrxx; ++ir) - // { - // v2sigma2_sigma_r[ir] = v2s2.at(ir * 6) * sigma[ir * 3] * 4.0 - // + v2s2.at(ir * 6 + 1) * sigma[ir * 3 + 1] * 4.0 - // + v2s2.at(ir * 6 + 3) * sigma[ir * 3 + 2]; //up-up - // v2sigma2_sigma_r[nrxx + ir] = v2s2.at(ir * 6 + 1) * sigma[ir * 3] * 2.0 - // + v2s2.at(ir * 6 + 4) * sigma[ir * 3 + 2] * 2.0 - // + (v2s2.at(ir * 6 + 2) * 4.0 + v2s2.at(ir * 6 + 3)) * sigma[ir * 3 + 1]; //up-down - // v2sigma2_sigma_r[2 * nrxx + ir] = v2s2.at(ir * 6 + 5) * sigma[ir * 3 + 2] * 4.0 - // + v2s2.at(ir * 6 + 4) * sigma[ir * 3 + 1] * 4.0 - // + 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); - // // 3. $\nabla^2(v^\sigma)$ - // std::vector lap_vsigma(3 * nrxx); - // for (int ir = 0;ir < nrxx;++ir) - // { - // lap_vsigma[ir] = vs.at(ir * 3) * 2.0; - // lap_vsigma[nrxx + ir] = vs.at(ir * 3 + 1); - // 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); - // // 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); - // BlasConnector::axpy(3 * nrxx, 1.0, v2sigma2_sigma_r.data(), 1, to_mul_rho_.data(), 1); - // BlasConnector::axpy(3 * nrxx, 1.0, lap_vsigma.data(), 1, to_mul_rho_.data(), 1); - // } - else - { - throw std::domain_error("nspin =" + std::to_string(nspin) - + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); - } - } - } // end for( xc_func_type &func : funcs ) - XC_Functional_Libxc::finish_func(funcs); - - if (1 == PARAM.inp.nspin || 2 == PARAM.inp.nspin) return; - // else if (4 == PARAM.inp.nspin) - else//NSPIN != 1,2,4 is not supported - { - throw std::domain_error("PARAM.inp.nspin =" + std::to_string(PARAM.inp.nspin) - + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); - } -} -#endif \ No newline at end of file diff --git a/source/module_lr/potentials/pot_hxc_lrtd.h b/source/module_lr/potentials/pot_hxc_lrtd.h index f688c41aef..00955e0a19 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.h +++ b/source/module_lr/potentials/pot_hxc_lrtd.h @@ -1,7 +1,7 @@ #pragma once #include "module_elecstate/potentials/pot_base.h" #include "module_elecstate/potentials/H_Hartree_pw.h" -#include "kernel.h" +#include "xc_kernel.h" #include #include diff --git a/source/module_lr/potentials/xc_kernel.cpp b/source/module_lr/potentials/xc_kernel.cpp new file mode 100644 index 0000000000..b901a752ce --- /dev/null +++ b/source/module_lr/potentials/xc_kernel.cpp @@ -0,0 +1,221 @@ +#include "xc_kernel.h" +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_parameter/parameter.h" +#include "module_base/timer.h" +#include "module_lr/utils/lr_util.h" +#ifdef USE_LIBXC +#include +#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 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 funcs = XC_Functional_Libxc::init_func( + XC_Functional::get_func_id(), + (1 == nspin) ? XC_UNPOLARIZED : XC_POLARIZED); + const int& nrxx = rho_basis_.nrxx; + + // converting rho (extract it as a subfuntion in the future) + // ----------------------------------------------------------------------------------- + std::vector rho(nspin * nrxx); // r major / spin contigous + +#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] = rho_gs[is][ir]; } } + if (rho_core) + { + 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]; } } + } + + // ----------------------------------------------------------------------------------- + // for GGA + const bool is_gga = [&funcs]() + { + for (xc_func_type& func : funcs) + { + switch (func.info->family) + { + case XC_FAMILY_GGA: + case XC_FAMILY_HYB_GGA: + return true; + } + } + return false; + }(); + std::vector>> gdr; // \nabla \rho + std::vector sigma; // |\nabla\rho|^2 + std::vector sgn; // sgn for threshold mask + if (is_gga) + { + // 0. set up sgn for threshold mask + // in the case of GGA correlation for polarized case, + // a cutoff for grho is required to ensure that libxc gives reasonable results + + // 1. \nabla \rho + gdr.resize(nspin); + for (int is = 0; is < nspin; ++is) + { + std::vector rhor(nrxx); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#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(), rho_basis_, tpiba); + } + // 2. |\nabla\rho|^2 + sigma.resize(nrxx * ((1 == nspin) ? 1 : 3)); + if (1 == nspin) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for (int ir = 0; ir < nrxx; ++ir) + sigma[ir] = gdr[0][ir] * gdr[0][ir]; + } + else + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int ir = 0; ir < nrxx; ++ir) + { + sigma[ir * 3] = gdr[0][ir] * gdr[0][ir]; + sigma[ir * 3 + 1] = gdr[0][ir] * gdr[1][ir]; + sigma[ir * 3 + 2] = gdr[1][ir] * gdr[1][ir]; + } + } + } + // ----------------------------------------------------------------------------------- + //==================== XC Kernels (f_xc)============================= + this->kernel_set_.emplace("vrho", std::vector(nspin * nrxx)); + this->kernel_set_.emplace("v2rho2", std::vector(((1 == nspin) ? 1 : 3) * nrxx));//(nrxx* ((1 == nspin) ? 1 : 3)): 00, 01, 11 + if (is_gga) + { + this->kernel_set_.emplace("vsigma", std::vector(((1 == nspin) ? 1 : 3) * nrxx)); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 + this->kernel_set_.emplace("v2rhosigma", std::vector(((1 == nspin) ? 1 : 6) * nrxx)); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 + this->kernel_set_.emplace("v2sigma2", std::vector(((1 == nspin) ? 1 : 6) * nrxx)); //(nrxx* ((1 == nspin) ? 1 : 6)): 00, 01, 02, 11, 12, 22 + } + //MetaGGA ... + + for (xc_func_type& func : funcs) + { + const double rho_threshold = 1E-6; + const double grho_threshold = 1E-10; + + xc_func_set_dens_threshold(&func, rho_threshold); + + //cut off grho if not LDA (future subfunc) + + switch (func.info->family) + { + case XC_FAMILY_LDA: + xc_lda_vxc(&func, nrxx, rho.data(), this->kernel_set_["vrho"].data()); + xc_lda_fxc(&func, nrxx, rho.data(), this->kernel_set_["v2rho2"].data()); + break; + case XC_FAMILY_GGA: + case XC_FAMILY_HYB_GGA: + xc_gga_vxc(&func, nrxx, rho.data(), sigma.data(), + this->kernel_set_["vrho"].data(), this->kernel_set_["vsigma"].data()); + xc_gga_fxc(&func, nrxx, rho.data(), sigma.data(), + this->kernel_set_["v2rho2"].data(), this->kernel_set_["v2rhosigma"].data(), this->kernel_set_["v2sigma2"].data()); + break; + default: + throw std::domain_error("func.info->family =" + std::to_string(func.info->family) + + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); + break; + } + // some formulas for GGA + if (func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) + { + const std::vector& v2r2 = this->kernel_set_["v2rho2"]; + const std::vector& v2rs = this->kernel_set_["v2rhosigma"]; + const std::vector& v2s2 = this->kernel_set_["v2sigma2"]; + const std::vector& vs = this->kernel_set_["vsigma"]; + const double tpiba2 = tpiba * tpiba; + + if (1 == nspin) + { + // 0. drho + this->grad_kernel_set_.emplace("drho_gs", std::vector>(nrxx)); + for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["drho_gs"][ir] = gdr[0][ir]; + // 1. $2f^{\rho\sigma}*\nabla\rho$ + this->grad_kernel_set_.emplace("2_v2rhosigma_drho", std::vector>(nrxx)); + for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["2_v2rhosigma_drho"][ir] = gdr[0][ir] * v2rs.at(ir) * 2.; + // 2. $4f^{\sigma\sigma}*\nabla\rho$ + this->grad_kernel_set_.emplace("4_v2sigma2_drho", std::vector>(nrxx)); + for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["4_v2sigma2_drho"][ir] = sigma.at(ir) * v2s2.at(ir) * 4.; + } + // else if (2 == nspin) // wrong, to be fixed + // { + // // 1. $\nabla\cdot(f^{\rho\sigma}*\nabla\rho)$ + // std::vector div_v2rhosigma_gdrho_r(3 * nrxx); + // std::vector> v2rhosigma_gdrho_r(3 * nrxx); + // for (int ir = 0; ir < nrxx; ++ir) + // { + // v2rhosigma_gdrho_r[ir] = gdr[0][ir] * v2rs.at(ir * 6) * 4.0 + // + gdr[1][ir] * v2rs.at(ir * 6 + 1) * 2.0; //up-up + // v2rhosigma_gdrho_r[nrxx + ir] = gdr[0][ir] * (v2rs.at(ir * 6 + 3) * 2.0 + v2rs.at(ir * 6 + 1)) + // + gdr[1][ir] * (v2rs.at(ir * 6 + 2) * 2.0 + v2rs.at(ir * 6 + 4)); //up-down + // v2rhosigma_gdrho_r[2 * nrxx + ir] = gdr[1][ir] * v2rs.at(ir * 6 + 5) * 4.0 + // + gdr[0][ir] * v2rs.at(ir * 6 + 4) * 2.0; //down-down + // } + // for (int isig = 0;isig < 3;++isig) + // XC_Functional::grad_dot(v2rhosigma_gdrho_r.data() + isig * nrxx, div_v2rhosigma_gdrho_r.data() + isig * nrxx, chg_gs->rhopw, tpiba); + // // 2. $\nabla^2(f^{\sigma\sigma}*\sigma)$ + // std::vector v2sigma2_sigma_r(3 * nrxx); + // for (int ir = 0; ir < nrxx; ++ir) + // { + // v2sigma2_sigma_r[ir] = v2s2.at(ir * 6) * sigma[ir * 3] * 4.0 + // + v2s2.at(ir * 6 + 1) * sigma[ir * 3 + 1] * 4.0 + // + v2s2.at(ir * 6 + 3) * sigma[ir * 3 + 2]; //up-up + // v2sigma2_sigma_r[nrxx + ir] = v2s2.at(ir * 6 + 1) * sigma[ir * 3] * 2.0 + // + v2s2.at(ir * 6 + 4) * sigma[ir * 3 + 2] * 2.0 + // + (v2s2.at(ir * 6 + 2) * 4.0 + v2s2.at(ir * 6 + 3)) * sigma[ir * 3 + 1]; //up-down + // v2sigma2_sigma_r[2 * nrxx + ir] = v2s2.at(ir * 6 + 5) * sigma[ir * 3 + 2] * 4.0 + // + v2s2.at(ir * 6 + 4) * sigma[ir * 3 + 1] * 4.0 + // + 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); + // // 3. $\nabla^2(v^\sigma)$ + // std::vector lap_vsigma(3 * nrxx); + // for (int ir = 0;ir < nrxx;++ir) + // { + // lap_vsigma[ir] = vs.at(ir * 3) * 2.0; + // lap_vsigma[nrxx + ir] = vs.at(ir * 3 + 1); + // 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); + // // 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); + // BlasConnector::axpy(3 * nrxx, 1.0, v2sigma2_sigma_r.data(), 1, to_mul_rho_.data(), 1); + // BlasConnector::axpy(3 * nrxx, 1.0, lap_vsigma.data(), 1, to_mul_rho_.data(), 1); + // } + else + { + throw std::domain_error("nspin =" + std::to_string(nspin) + + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); + } + } + } // end for( xc_func_type &func : funcs ) + XC_Functional_Libxc::finish_func(funcs); + + if (1 == PARAM.inp.nspin || 2 == PARAM.inp.nspin) return; + // else if (4 == PARAM.inp.nspin) + else//NSPIN != 1,2,4 is not supported + { + throw std::domain_error("PARAM.inp.nspin =" + std::to_string(PARAM.inp.nspin) + + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); + } +} +#endif \ No newline at end of file diff --git a/source/module_lr/potentials/kernel.h b/source/module_lr/potentials/xc_kernel.h similarity index 100% rename from source/module_lr/potentials/kernel.h rename to source/module_lr/potentials/xc_kernel.h From 85b1db054accc9c27eaff9a456c3235447b3cda1 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 14:28:03 +0800 Subject: [PATCH 05/23] rename init_fxc as init_xc_kernel --- docs/advanced/input_files/input-main.md | 10 +++++----- source/module_io/read_input_item_tddft.cpp | 8 ++++---- source/module_io/test/read_input_ptest.cpp | 2 +- source/module_lr/esolver_lrtd_lcao.cpp | 8 ++++---- source/module_lr/potentials/pot_hxc_lrtd.cpp | 19 ++++++++++--------- source/module_lr/potentials/pot_hxc_lrtd.h | 2 +- source/module_parameter/input_parameter.h | 2 +- 7 files changed, 26 insertions(+), 25 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 65e37d2e54..5c933fab75 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -3943,14 +3943,14 @@ 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 +### init_xc_kernel - **Type**: String - **Description**: The method to initalize the xc kernel. - - "gs": Calculate xc kerenel ($f_\text{xc}$) from the ground-state charge density. - - "file_fxc": 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 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" + - "from_chg_groundstate": 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 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. + - "from_chg_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**: "from_chg_groundstate" ### lr_solver diff --git a/source/module_io/read_input_item_tddft.cpp b/source/module_io/read_input_item_tddft.cpp index 2d6bc6c265..b407827222 100644 --- a/source/module_io/read_input_item_tddft.cpp +++ b/source/module_io/read_input_item_tddft.cpp @@ -310,17 +310,17 @@ void ReadInput::item_lr_tddft() this->add_item(item); } { - Input_Item item("init_fxc"); + Input_Item item("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.init_fxc; + auto& ifxc = para.input.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.init_fxc.empty()) { para.input.init_fxc.push_back("gs"); } + if (para.input.init_xc_kernel.empty()) { para.input.init_xc_kernel.push_back("from_chg_groundstate"); } }; - sync_stringvec(input.init_fxc, para.input.init_fxc.size(), "gs"); + sync_stringvec(input.init_xc_kernel, para.input.init_xc_kernel.size(), "from_chg_groundstate"); this->add_item(item); } { diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index 07a34485c0..c674ec51c3 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -420,7 +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.init_fxc[0], "gs"); + EXPECT_EQ(param.inp.init_xc_kernel[0], "from_chg_groundstate"); EXPECT_EQ(param.inp.lr_solver, "dav"); EXPECT_DOUBLE_EQ(param.inp.lr_thr, 1e-2); EXPECT_FALSE(param.inp.lr_unrestricted); diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index e956e51bae..3276f03bae 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -66,7 +66,7 @@ template void LR::ESolver_LR::parameter_check()const { std::set lr_solvers = { "dav", "lapack" , "spectrum", "dav_subspace", "cg" }; - std::set xc_kernels = { "rpa", "lda", "pbe", "hf" , "hse" }; + std::set 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"); } @@ -592,11 +592,11 @@ void LR::ESolver_LR::init_pot(const Charge& chg_gs) { using ST = PotHxcLR::SpinType; case 1: - this->pot[0] = std::make_shared(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, ST::S1, input.init_fxc); + this->pot[0] = std::make_shared(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, ST::S1, input.init_xc_kernel); break; case 2: - this->pot[0] = std::make_shared(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(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_triplet, input.init_fxc); + this->pot[0] = std::make_shared(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_singlet, input.init_xc_kernel); + this->pot[1] = std::make_shared(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_triplet, input.init_xc_kernel); break; default: throw std::invalid_argument("ESolver_LR: nspin must be 1 or 2"); diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index baf68d185c..a3094ae127 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -14,7 +14,7 @@ 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 Parallel_Grid& pgrid, - const SpinType& st_in, const std::vector& init_fxc) + const SpinType& st_in, const std::vector& init_xc_kernel) :xc_kernel(xc_kernel_in), tpiba_(ucell_in->tpiba), spin_type_(st_in), xc_kernel_components_(*rho_basis_in) { @@ -24,26 +24,27 @@ namespace LR this->pot_hartree = LR_Util::make_unique(this->rho_basis_); - std::set local_xc = { "lda", "pbe", "hse" }; + std::set local_xc = { "lda", "pwlda", "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()); this->set_integral_func(this->spin_type_, this->xc_type_); - if (init_fxc[0] == "file_fxc") + if (init_xc_kernel[0] == "file") { - assert(this->xc_kernel == "lda"); + std::set lda_xc = { "lda", "pwlda" }; + assert(lda_xc.count(this->xc_kernel)); const int spinsize = (1 == nspin) ? 1 : 3; std::vector v2rho2(spinsize * nrxx); // read fxc adn add to xc_kernel_components - assert(init_fxc.size() >= spinsize + 1); + assert(init_xc_kernel.size() >= spinsize + 1); for (int is = 0;is < spinsize;++is) { double ef = 0.0; int prenspin = 1; std::vector v2rho2_tmp(nrxx); - ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, init_fxc[is + 1], + ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, init_xc_kernel[is + 1], v2rho2_tmp.data(), ucell_in->nat); for (int ir = 0;ir < nrxx;++ir) { v2rho2[ir * spinsize + is] = v2rho2_tmp[ir]; } } @@ -52,16 +53,16 @@ namespace LR } #ifdef USE_LIBXC - if (init_fxc[0] == "file_chg") + if (init_xc_kernel[0] == "from_chg_file") { - assert(init_fxc.size() >= 2); + assert(init_xc_kernel.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]; + const std::string file = init_xc_kernel[init_xc_kernel.size() > nspin ? 1 + is : 1]; ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, file, rho_for_fxc[is], ucell_in->nat); } diff --git a/source/module_lr/potentials/pot_hxc_lrtd.h b/source/module_lr/potentials/pot_hxc_lrtd.h index 00955e0a19..6dee1d05b0 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.h +++ b/source/module_lr/potentials/pot_hxc_lrtd.h @@ -19,7 +19,7 @@ namespace LR /// 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 Parallel_Grid& pgrid, - const SpinType& st_in = SpinType::S1, const std::vector& init_fxc = { "gs" }); + const SpinType& st_in = SpinType::S1, const std::vector& init_xc_kernel = { "from_chg_groundstate" }); ~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& ispin_op = { 0,0 }); diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index 4e62414a03..d564289345 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -299,7 +299,7 @@ struct Input_para // ============== #Parameters (10.lr-tddft) =========================== int lr_nstates = 1; ///< the number of 2-particle states to be solved - std::vector init_fxc = {}; ///< The method to initalize the xc kernel + std::vector init_xc_kernel = {}; ///< 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 From e2f7b636b50c749e206b4543e4a0abb009004dbe Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 14:44:28 +0800 Subject: [PATCH 06/23] update grad and lapl --- source/Makefile.Objects | 1 - source/module_lr/CMakeLists.txt | 1 - source/module_lr/potentials/pot_hxc_lrtd.cpp | 1 + source/module_lr/potentials/xc_kernel.cpp | 5 +- source/module_lr/utils/lr_util.h | 19 ++++---- source/module_lr/utils/lr_util_xc.cpp | 28 ------------ source/module_lr/utils/lr_util_xc.hpp | 48 ++++++++++++++++++++ 7 files changed, 62 insertions(+), 41 deletions(-) delete mode 100644 source/module_lr/utils/lr_util_xc.cpp create mode 100644 source/module_lr/utils/lr_util_xc.hpp diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 89fc9effd7..d34d08242e 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -724,7 +724,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\ diff --git a/source/module_lr/CMakeLists.txt b/source/module_lr/CMakeLists.txt index 93b4f7ff9f..47a8874aef 100644 --- a/source/module_lr/CMakeLists.txt +++ b/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 diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index a3094ae127..054aaa6a09 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -6,6 +6,7 @@ #include "module_hamilt_general/module_xc/xc_functional.h" #include #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& ispin_op = { 0,0 } diff --git a/source/module_lr/potentials/xc_kernel.cpp b/source/module_lr/potentials/xc_kernel.cpp index b901a752ce..caa0e6dfdd 100644 --- a/source/module_lr/potentials/xc_kernel.cpp +++ b/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 #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 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); diff --git a/source/module_lr/utils/lr_util.h b/source/module_lr/utils/lr_util.h index f6f70014db..1c47f6fff0 100644 --- a/source/module_lr/utils/lr_util.h +++ b/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 struct ToComplex; +template <> struct ToComplex { using type = std::complex; }; +template <> struct ToComplex> { using type = std::complex; }; +template <> struct ToComplex { using type = std::complex; }; +template <> struct ToComplex> { using type = std::complex; }; +#endif + namespace LR_Util { /// =====================PHYSICS==================== @@ -36,15 +45,7 @@ namespace LR_Util std::pair>> set_ix_map_diagonal(bool mode, int nc, int nv); - /// operators to calculate XC kernels - void grad(const double* rhor, - ModuleBase::Vector3* 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======== diff --git a/source/module_lr/utils/lr_util_xc.cpp b/source/module_lr/utils/lr_util_xc.cpp deleted file mode 100644 index cc8b3a1f6d..0000000000 --- a/source/module_lr/utils/lr_util_xc.cpp +++ /dev/null @@ -1,28 +0,0 @@ -#include "lr_util.h" -namespace LR_Util -{ - void grad(const double* rhor, - ModuleBase::Vector3* gdr, - const ModulePW::PW_Basis& rho_basis, - const double& tpiba) - { - std::vector> 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> rhog(rho_basis.npw); - std::vector 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; } - } - } -} diff --git a/source/module_lr/utils/lr_util_xc.hpp b/source/module_lr/utils/lr_util_xc.hpp new file mode 100644 index 0000000000..229cbf7335 --- /dev/null +++ b/source/module_lr/utils/lr_util_xc.hpp @@ -0,0 +1,48 @@ +#pragma once +#include "lr_util.h" +namespace LR_Util +{ + template + void grad(const T* rhor, + ModuleBase::Vector3* gdr, + const ModulePW::PW_Basis& rho_basis, + const double& tpiba) + { + std::vector::type> rhog(rho_basis.npw); + rho_basis.real2recip(rhor, rhog.data()); + XC_Functional::grad_rho(rhog.data(), gdr, &rho_basis, tpiba); + } + template + void grad(const std::vector& rhor, + std::vector>& gdr, + const ModulePW::PW_Basis& rho_basis, + const double& tpiba) + { + grad(rhor.data(), gdr.data(), rho_basis, tpiba); + } + + template + 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::type> rhog(rho_basis.npw); + std::vector 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 + void lapl(const std::vector& rhor, + std::vector& lapn, + const ModulePW::PW_Basis& rho_basis, + const double& tpiba2) + { + lapl(rhor.data(), lapn.data(), rho_basis, tpiba2); + } +} From e87340be5453abdf99ccbae22240935b5436cb15 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 15:04:07 +0800 Subject: [PATCH 07/23] delete useless matrix transformer --- source/module_lr/utils/lr_util.cpp | 108 ++++------------------------- source/module_lr/utils/lr_util.h | 14 ---- 2 files changed, 12 insertions(+), 110 deletions(-) diff --git a/source/module_lr/utils/lr_util.cpp b/source/module_lr/utils/lr_util.cpp index 4ae65fcaea..1418764126 100644 --- a/source/module_lr/utils/lr_util.cpp +++ b/source/module_lr/utils/lr_util.cpp @@ -62,7 +62,7 @@ namespace LR_Util template<> void matsym(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]; } + std::copy(in, in + pmat.get_local_size(), out); 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); @@ -79,105 +79,21 @@ namespace LR_Util template<> void matsym>(const std::complex* in, const int n, const Parallel_2D& pmat, std::complex* out) { - for (int i = 0;i < pmat.get_local_size();++i) { out[i] = in[i]; } + std::copy(in, in + pmat.get_local_size(), out); + const std::complex alpha(0.5, 0.0), beta(0.5, 0.0); const int i1 = 1; - std::vector> tmp(pmat.get_local_size()); - std::copy(inout, inout + pmat.get_local_size(), tmp.begin()); - container::Tensor t(DAT::DT_DOUBLE, DEV::CpuDevice, { m.nr, m.nc }); - for (int i = 0;i < t.NumElements();++i) { t.data()[i] = m.c[i]; } - return t; - } - std::vector mat2ten_double(std::vector& m) - { - std::vector t; - 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()[i]; - } - return m; - } - std::vector ten2mat_double(std::vector& t) - { - std::vector m; - 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>()[i] = m.c[i]; - } - return t; - } - std::vector mat2ten_complex(std::vector& m) - { - std::vector t; - 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>()[i]; - } - return m; + pztranc_(&n, &n, &alpha, in, &i1, &i1, pmat.desc, &beta, out, &i1, &i1, pmat.desc); } - std::vector ten2mat_complex(std::vector& t) - { - std::vector m; - for (int i = 0;i < t.size();++i) { - m.push_back(ten2mat_complex(t[i])); - } - return m; - } - - ModuleBase::matrix vec2mat(const std::vector& 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]; - } - return m; - } - ModuleBase::ComplexMatrix vec2mat(const std::vector>& 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]; - } - return m; - } - std::vector vec2mat(const std::vector>& v, const int nr, const int nc) - { - std::vector m(v.size()); - for (int i = 0;i < v.size();++i) { - m[i] = vec2mat(v[i], nr, nc); - } - return m; - } - std::vector vec2mat(const std::vector>>& v, const int nr, const int nc) + template<> + void matsym>(std::complex* inout, const int n, const Parallel_2D& pmat) { - std::vector m(v.size()); - for (int i = 0;i < v.size();++i) { - m[i] = vec2mat(v[i], nr, nc); - } - return m; + std::vector> tmp(pmat.get_local_size()); + std::copy(inout, inout + pmat.get_local_size(), tmp.begin()); + const std::complex 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); } +#endif // for the first matrix in the commutator void setup_2d_division(Parallel_2D& pv, int nb, int gr, int gc) diff --git a/source/module_lr/utils/lr_util.h b/source/module_lr/utils/lr_util.h index 1c47f6fff0..892a42ed7e 100644 --- a/source/module_lr/utils/lr_util.h +++ b/source/module_lr/utils/lr_util.h @@ -74,20 +74,6 @@ namespace LR_Util template void matsym(T* inout, const int n, const Parallel_2D& pmat); #endif - ///======== Tensor - Matrix transformer========== - container::Tensor mat2ten_double(ModuleBase::matrix& m); - std::vector mat2ten_double(std::vector& m); - ModuleBase::matrix ten2mat_double(container::Tensor& t); - std::vector ten2mat_double(std::vector& t); - container::Tensor mat2ten_complex(ModuleBase::ComplexMatrix& m); - std::vector mat2ten_complex(std::vector& m); - ModuleBase::ComplexMatrix ten2mat_complex(container::Tensor& t); - std::vector ten2mat_complex(std::vector& t); - - ModuleBase::matrix vec2mat(const std::vector& v, const int nr, const int nc); - ModuleBase::ComplexMatrix vec2mat(const std::vector>& v, const int nr, const int nc); - std::vector vec2mat(const std::vector>& v, const int nr, const int nc); - std::vector vec2mat(const std::vector>>& v, const int nr, const int nc); ///===================Psi wrapper================= /// psi(nk=1, nbands=nb, nk * nbasis) -> psi(nb, nk, nbasis) without memory copy From 68429a5dde0cdf27c2c7e7c12c5007a048db1a18 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 15:31:59 +0800 Subject: [PATCH 08/23] update comments --- source/module_lr/potentials/pot_hxc_lrtd.h | 1 + source/module_lr/potentials/xc_kernel.h | 4 +++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/source/module_lr/potentials/pot_hxc_lrtd.h b/source/module_lr/potentials/pot_hxc_lrtd.h index 6dee1d05b0..c37cdd216b 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.h +++ b/source/module_lr/potentials/pot_hxc_lrtd.h @@ -14,6 +14,7 @@ namespace LR /// S1: K^Hartree + K^xc /// S2_singlet: 2*K^Hartree + K^xc_{upup} + K^xc_{updown} /// S2_triplet: K^xc_{upup} - K^xc_{updown} + /// S2_updown: K^Hartree + (K^xc_{upup}, K^xc_{updown}, K^xc_{downup} or K^xc_{downdown}), according to `ispin_op` (for spin-polarized systems) 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 diff --git a/source/module_lr/potentials/xc_kernel.h b/source/module_lr/potentials/xc_kernel.h index 38da6087a4..097c91be9d 100644 --- a/source/module_lr/potentials/xc_kernel.h +++ b/source/module_lr/potentials/xc_kernel.h @@ -3,17 +3,19 @@ namespace LR { + /// @brief Calculate the exchange-correlation (XC) kernel ($f_{xc}=\delta^2E_xc/\delta\rho^2$) and store its components. class KernelXC { public: KernelXC(const ModulePW::PW_Basis& rho_basis) : rho_basis_(rho_basis) {}; ~KernelXC() {}; - // xc kernel for LR-TDDFT #ifdef USE_LIBXC + /// @brief Calculate the XC kernel using libxc. 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 + // Setters and getters for certain components, such as v2rho2, v2rhosigma, v2sigma2, or other custom combinations. void set_kernel(const std::string& name, const std::vector& vec) { this->kernel_set_[name] = vec; } void set_kernel(const std::string& name, const std::vector&& vec) { this->kernel_set_[name] = std::move(vec); } const std::vector& get_kernel(const std::string& name) { return kernel_set_[name]; } From 8cdd9c5fbfacfa29762a6d843ac10cdfb73fa12c Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 15:33:14 +0800 Subject: [PATCH 09/23] rename: lr_init_xc_kernel --- docs/advanced/input_files/input-main.md | 2 +- source/module_io/read_input_item_tddft.cpp | 8 ++++---- source/module_io/test/read_input_ptest.cpp | 2 +- source/module_lr/esolver_lrtd_lcao.cpp | 6 +++--- source/module_lr/potentials/pot_hxc_lrtd.cpp | 14 +++++++------- source/module_lr/potentials/pot_hxc_lrtd.h | 2 +- source/module_parameter/input_parameter.h | 2 +- 7 files changed, 18 insertions(+), 18 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 5c933fab75..d2f1c30211 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -3943,7 +3943,7 @@ These parameters are used to solve the excited states using. e.g. LR-TDDFT. Currently supported: `RPA`, `LDA`, `PBE`, `HSE`, `HF`. - **Default**: LDA -### init_xc_kernel +### lr_init_xc_kernel - **Type**: String - **Description**: The method to initalize the xc kernel. diff --git a/source/module_io/read_input_item_tddft.cpp b/source/module_io/read_input_item_tddft.cpp index b407827222..ce71273dc7 100644 --- a/source/module_io/read_input_item_tddft.cpp +++ b/source/module_io/read_input_item_tddft.cpp @@ -310,17 +310,17 @@ void ReadInput::item_lr_tddft() this->add_item(item); } { - Input_Item item("init_xc_kernel"); + 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.init_xc_kernel; + 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.init_xc_kernel.empty()) { para.input.init_xc_kernel.push_back("from_chg_groundstate"); } + if (para.input.lr_init_xc_kernel.empty()) { para.input.lr_init_xc_kernel.push_back("from_chg_groundstate"); } }; - sync_stringvec(input.init_xc_kernel, para.input.init_xc_kernel.size(), "from_chg_groundstate"); + sync_stringvec(input.lr_init_xc_kernel, para.input.lr_init_xc_kernel.size(), "from_chg_groundstate"); this->add_item(item); } { diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index c674ec51c3..9cf702c2ac 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -420,7 +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.init_xc_kernel[0], "from_chg_groundstate"); + EXPECT_EQ(param.inp.lr_init_xc_kernel[0], "from_chg_groundstate"); EXPECT_EQ(param.inp.lr_solver, "dav"); EXPECT_DOUBLE_EQ(param.inp.lr_thr, 1e-2); EXPECT_FALSE(param.inp.lr_unrestricted); diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 3276f03bae..d741752bcb 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -592,11 +592,11 @@ void LR::ESolver_LR::init_pot(const Charge& chg_gs) { using ST = PotHxcLR::SpinType; case 1: - this->pot[0] = std::make_shared(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, ST::S1, input.init_xc_kernel); + this->pot[0] = std::make_shared(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(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_singlet, input.init_xc_kernel); - this->pot[1] = std::make_shared(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_triplet, input.init_xc_kernel); + this->pot[0] = std::make_shared(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(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"); diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index 054aaa6a09..39b90775e6 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -15,7 +15,7 @@ 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 Parallel_Grid& pgrid, - const SpinType& st_in, const std::vector& init_xc_kernel) + const SpinType& st_in, const std::vector& lr_init_xc_kernel) :xc_kernel(xc_kernel_in), tpiba_(ucell_in->tpiba), spin_type_(st_in), xc_kernel_components_(*rho_basis_in) { @@ -32,20 +32,20 @@ namespace LR this->xc_type_ = XCType(XC_Functional::get_func_type()); this->set_integral_func(this->spin_type_, this->xc_type_); - if (init_xc_kernel[0] == "file") + if (lr_init_xc_kernel[0] == "file") { std::set lda_xc = { "lda", "pwlda" }; assert(lda_xc.count(this->xc_kernel)); const int spinsize = (1 == nspin) ? 1 : 3; std::vector v2rho2(spinsize * nrxx); // read fxc adn add to xc_kernel_components - assert(init_xc_kernel.size() >= spinsize + 1); + assert(lr_init_xc_kernel.size() >= spinsize + 1); for (int is = 0;is < spinsize;++is) { double ef = 0.0; int prenspin = 1; std::vector v2rho2_tmp(nrxx); - ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, init_xc_kernel[is + 1], + ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, lr_init_xc_kernel[is + 1], v2rho2_tmp.data(), ucell_in->nat); for (int ir = 0;ir < nrxx;++ir) { v2rho2[ir * spinsize + is] = v2rho2_tmp[ir]; } } @@ -54,16 +54,16 @@ namespace LR } #ifdef USE_LIBXC - if (init_xc_kernel[0] == "from_chg_file") + if (lr_init_xc_kernel[0] == "from_chg_file") { - assert(init_xc_kernel.size() >= 2); + assert(lr_init_xc_kernel.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_xc_kernel[init_xc_kernel.size() > nspin ? 1 + is : 1]; + const std::string file = lr_init_xc_kernel[lr_init_xc_kernel.size() > nspin ? 1 + is : 1]; ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, file, rho_for_fxc[is], ucell_in->nat); } diff --git a/source/module_lr/potentials/pot_hxc_lrtd.h b/source/module_lr/potentials/pot_hxc_lrtd.h index c37cdd216b..11d15474f2 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.h +++ b/source/module_lr/potentials/pot_hxc_lrtd.h @@ -20,7 +20,7 @@ namespace LR /// 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 Parallel_Grid& pgrid, - const SpinType& st_in = SpinType::S1, const std::vector& init_xc_kernel = { "from_chg_groundstate" }); + const SpinType& st_in = SpinType::S1, const std::vector& lr_init_xc_kernel = { "from_chg_groundstate" }); ~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& ispin_op = { 0,0 }); diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index d564289345..ac4d74fca7 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -299,7 +299,7 @@ struct Input_para // ============== #Parameters (10.lr-tddft) =========================== int lr_nstates = 1; ///< the number of 2-particle states to be solved - std::vector init_xc_kernel = {}; ///< The method to initalize the xc kernel + std::vector lr_init_xc_kernel = {}; ///< 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 From b4b1cb6bb50e2ba1c48ef5d53ff7b65be4078874 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 15:53:53 +0800 Subject: [PATCH 10/23] add const before std::set --- source/module_lr/esolver_lrtd_lcao.cpp | 4 ++-- source/module_lr/potentials/pot_hxc_lrtd.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index d741752bcb..4a37f3b99d 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -65,8 +65,8 @@ inline int cal_nupdown_form_occ(const ModuleBase::matrix& wg) template void LR::ESolver_LR::parameter_check()const { - std::set lr_solvers = { "dav", "lapack" , "spectrum", "dav_subspace", "cg" }; - std::set xc_kernels = { "rpa", "lda", "pwlda", "pbe", "hf" , "hse" }; + const std::set lr_solvers = { "dav", "lapack" , "spectrum", "dav_subspace", "cg" }; + const std::set 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"); } diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index 39b90775e6..98fc2c7218 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -25,7 +25,7 @@ namespace LR this->pot_hartree = LR_Util::make_unique(this->rho_basis_); - std::set local_xc = { "lda", "pwlda", "pbe", "hse" }; + const std::set local_xc = { "lda", "pwlda", "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 @@ -34,7 +34,7 @@ namespace LR if (lr_init_xc_kernel[0] == "file") { - std::set lda_xc = { "lda", "pwlda" }; + const std::set lda_xc = { "lda", "pwlda" }; assert(lda_xc.count(this->xc_kernel)); const int spinsize = (1 == nspin) ? 1 : 3; std::vector v2rho2(spinsize * nrxx); From 31364df57e4ef1a0539c60bfef74c2e7857ff721 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 16:00:04 +0800 Subject: [PATCH 11/23] rename variables in PotHxcLR --- source/module_lr/potentials/pot_hxc_lrtd.cpp | 26 ++++++++++---------- source/module_lr/potentials/pot_hxc_lrtd.h | 8 +++--- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index 98fc2c7218..a7189c4ba8 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -13,29 +13,29 @@ 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, + PotHxcLR::PotHxcLR(const std::string& xc_kernel, const ModulePW::PW_Basis* rho_basis, const UnitCell* ucell, const Charge* chg_gs/*ground state*/, const Parallel_Grid& pgrid, - const SpinType& st_in, const std::vector& lr_init_xc_kernel) - :xc_kernel(xc_kernel_in), tpiba_(ucell_in->tpiba), spin_type_(st_in), - xc_kernel_components_(*rho_basis_in) + const SpinType& st, const std::vector& lr_init_xc_kernel) + :xc_kernel_(xc_kernel), tpiba_(ucell->tpiba), spin_type_(st), + xc_kernel_components_(*rho_basis) { - this->rho_basis_ = rho_basis_in; + this->rho_basis_ = rho_basis; 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(this->rho_basis_); const std::set local_xc = { "lda", "pwlda", "pbe", "hse" }; - if (local_xc.find(this->xc_kernel) != local_xc.end()) + 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 + 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()); this->set_integral_func(this->spin_type_, this->xc_type_); if (lr_init_xc_kernel[0] == "file") { const std::set lda_xc = { "lda", "pwlda" }; - assert(lda_xc.count(this->xc_kernel)); + assert(lda_xc.count(this->xc_kernel_)); const int spinsize = (1 == nspin) ? 1 : 3; std::vector v2rho2(spinsize * nrxx); // read fxc adn add to xc_kernel_components @@ -46,7 +46,7 @@ namespace LR int prenspin = 1; std::vector v2rho2_tmp(nrxx); ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, lr_init_xc_kernel[is + 1], - v2rho2_tmp.data(), ucell_in->nat); + v2rho2_tmp.data(), ucell->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)); @@ -65,12 +65,12 @@ namespace LR { const std::string file = lr_init_xc_kernel[lr_init_xc_kernel.size() > nspin ? 1 + is : 1]; ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, file, - rho_for_fxc[is], ucell_in->nat); + rho_for_fxc[is], ucell->nat); } - this->xc_kernel_components_.f_xc_libxc(nspin, ucell_in->omega, ucell_in->tpiba, rho_for_fxc, chg_gs->rho_core); + this->xc_kernel_components_.f_xc_libxc(nspin, ucell->omega, ucell->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 { this->xc_kernel_components_.f_xc_libxc(nspin, ucell->omega, ucell->tpiba, chg_gs->rho, chg_gs->rho_core); } #else ModuleBase::WARNING_QUIT("KernelXC", "to calculate xc-kernel in LR-TDDFT, compile with LIBXC"); #endif @@ -96,7 +96,7 @@ namespace LR break; } // XC - if (xc_kernel == "rpa" || xc_kernel == "hf") { return; } // no xc + if (this->xc_kernel_ == "rpa" || this->xc_kernel_ == "hf") { return; } // no xc #ifdef USE_LIBXC this->kernel_to_potential_[spin_type_](rho[0], v_eff, ispin_op); #else diff --git a/source/module_lr/potentials/pot_hxc_lrtd.h b/source/module_lr/potentials/pot_hxc_lrtd.h index 11d15474f2..59d8d71dbe 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.h +++ b/source/module_lr/potentials/pot_hxc_lrtd.h @@ -18,9 +18,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 Parallel_Grid& pgrid, - const SpinType& st_in = SpinType::S1, const std::vector& lr_init_xc_kernel = { "from_chg_groundstate" }); + PotHxcLR(const std::string& xc_kernel, const ModulePW::PW_Basis* rho_basis, + const UnitCell* ucell, const Charge* chg_gs/*ground state*/, const Parallel_Grid& pgrid, + const SpinType& st = SpinType::S1, const std::vector& lr_init_xc_kernel = { "from_chg_groundstate" }); ~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& ispin_op = { 0,0 }); @@ -33,7 +33,7 @@ namespace LR /// GGA: v2rho2, v2rhosigma, v2sigma2 /// meta-GGA: v2rho2, v2rhosigma, v2sigma2, v2rholap, v2rhotau, v2sigmalap, v2sigmatau, v2laptau, v2lap2, v2tau2 KernelXC xc_kernel_components_; - const std::string xc_kernel; + const std::string xc_kernel_; const double& tpiba_; const SpinType spin_type_ = SpinType::S1; XCType xc_type_ = XCType::None; From 504e1594dc67555abed918d62e4f6979044158ac Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 16:19:07 +0800 Subject: [PATCH 12/23] rename gdr as gradrho --- source/module_lr/potentials/xc_kernel.cpp | 32 +++++++++++------------ source/module_lr/utils/lr_util_xc.hpp | 8 +++--- 2 files changed, 20 insertions(+), 20 deletions(-) diff --git a/source/module_lr/potentials/xc_kernel.cpp b/source/module_lr/potentials/xc_kernel.cpp index caa0e6dfdd..fb3bd6aa7a 100644 --- a/source/module_lr/potentials/xc_kernel.cpp +++ b/source/module_lr/potentials/xc_kernel.cpp @@ -50,7 +50,7 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl } return false; }(); - std::vector>> gdr; // \nabla \rho + std::vector>> gradrho; // \nabla \rho std::vector sigma; // |\nabla\rho|^2 std::vector sgn; // sgn for threshold mask if (is_gga) @@ -60,7 +60,7 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl // a cutoff for grho is required to ensure that libxc gives reasonable results // 1. \nabla \rho - gdr.resize(nspin); + gradrho.resize(nspin); for (int is = 0; is < nspin; ++is) { std::vector rhor(nrxx); @@ -68,8 +68,8 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl #pragma omp parallel for schedule(static, 1024) #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(), rho_basis_, tpiba); + gradrho[is].resize(nrxx); + LR_Util::grad(rhor.data(), gradrho[is].data(), rho_basis_, tpiba); } // 2. |\nabla\rho|^2 sigma.resize(nrxx * ((1 == nspin) ? 1 : 3)); @@ -79,7 +79,7 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl #pragma omp parallel for schedule(static, 1024) #endif for (int ir = 0; ir < nrxx; ++ir) - sigma[ir] = gdr[0][ir] * gdr[0][ir]; + sigma[ir] = gradrho[0][ir] * gradrho[0][ir]; } else { @@ -88,9 +88,9 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl #endif for (int ir = 0; ir < nrxx; ++ir) { - sigma[ir * 3] = gdr[0][ir] * gdr[0][ir]; - sigma[ir * 3 + 1] = gdr[0][ir] * gdr[1][ir]; - sigma[ir * 3 + 2] = gdr[1][ir] * gdr[1][ir]; + sigma[ir * 3] = gradrho[0][ir] * gradrho[0][ir]; + sigma[ir * 3 + 1] = gradrho[0][ir] * gradrho[1][ir]; + sigma[ir * 3 + 2] = gradrho[1][ir] * gradrho[1][ir]; } } } @@ -146,10 +146,10 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl { // 0. drho this->grad_kernel_set_.emplace("drho_gs", std::vector>(nrxx)); - for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["drho_gs"][ir] = gdr[0][ir]; + for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["drho_gs"][ir] = gradrho[0][ir]; // 1. $2f^{\rho\sigma}*\nabla\rho$ this->grad_kernel_set_.emplace("2_v2rhosigma_drho", std::vector>(nrxx)); - for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["2_v2rhosigma_drho"][ir] = gdr[0][ir] * v2rs.at(ir) * 2.; + for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["2_v2rhosigma_drho"][ir] = gradrho[0][ir] * v2rs.at(ir) * 2.; // 2. $4f^{\sigma\sigma}*\nabla\rho$ this->grad_kernel_set_.emplace("4_v2sigma2_drho", std::vector>(nrxx)); for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["4_v2sigma2_drho"][ir] = sigma.at(ir) * v2s2.at(ir) * 4.; @@ -161,12 +161,12 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl // std::vector> v2rhosigma_gdrho_r(3 * nrxx); // for (int ir = 0; ir < nrxx; ++ir) // { - // v2rhosigma_gdrho_r[ir] = gdr[0][ir] * v2rs.at(ir * 6) * 4.0 - // + gdr[1][ir] * v2rs.at(ir * 6 + 1) * 2.0; //up-up - // v2rhosigma_gdrho_r[nrxx + ir] = gdr[0][ir] * (v2rs.at(ir * 6 + 3) * 2.0 + v2rs.at(ir * 6 + 1)) - // + gdr[1][ir] * (v2rs.at(ir * 6 + 2) * 2.0 + v2rs.at(ir * 6 + 4)); //up-down - // v2rhosigma_gdrho_r[2 * nrxx + ir] = gdr[1][ir] * v2rs.at(ir * 6 + 5) * 4.0 - // + gdr[0][ir] * v2rs.at(ir * 6 + 4) * 2.0; //down-down + // v2rhosigma_gdrho_r[ir] = gradrho[0][ir] * v2rs.at(ir * 6) * 4.0 + // + gradrho[1][ir] * v2rs.at(ir * 6 + 1) * 2.0; //up-up + // v2rhosigma_gdrho_r[nrxx + ir] = gradrho[0][ir] * (v2rs.at(ir * 6 + 3) * 2.0 + v2rs.at(ir * 6 + 1)) + // + gradrho[1][ir] * (v2rs.at(ir * 6 + 2) * 2.0 + v2rs.at(ir * 6 + 4)); //up-down + // v2rhosigma_gdrho_r[2 * nrxx + ir] = gradrho[1][ir] * v2rs.at(ir * 6 + 5) * 4.0 + // + gradrho[0][ir] * v2rs.at(ir * 6 + 4) * 2.0; //down-down // } // for (int isig = 0;isig < 3;++isig) // XC_Functional::grad_dot(v2rhosigma_gdrho_r.data() + isig * nrxx, div_v2rhosigma_gdrho_r.data() + isig * nrxx, chg_gs->rhopw, tpiba); diff --git a/source/module_lr/utils/lr_util_xc.hpp b/source/module_lr/utils/lr_util_xc.hpp index 229cbf7335..a270482ba7 100644 --- a/source/module_lr/utils/lr_util_xc.hpp +++ b/source/module_lr/utils/lr_util_xc.hpp @@ -4,21 +4,21 @@ namespace LR_Util { template void grad(const T* rhor, - ModuleBase::Vector3* gdr, + ModuleBase::Vector3* gradrho, const ModulePW::PW_Basis& rho_basis, const double& tpiba) { std::vector::type> rhog(rho_basis.npw); rho_basis.real2recip(rhor, rhog.data()); - XC_Functional::grad_rho(rhog.data(), gdr, &rho_basis, tpiba); + XC_Functional::grad_rho(rhog.data(), gradrho, &rho_basis, tpiba); } template void grad(const std::vector& rhor, - std::vector>& gdr, + std::vector>& gradrho, const ModulePW::PW_Basis& rho_basis, const double& tpiba) { - grad(rhor.data(), gdr.data(), rho_basis, tpiba); + grad(rhor.data(), gradrho.data(), rho_basis, tpiba); } template From fb7b93e21c462e65efc727923330536e6e3655cb Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 17:09:06 +0800 Subject: [PATCH 13/23] rename spinsize as n_component --- source/module_lr/potentials/pot_hxc_lrtd.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index a7189c4ba8..9a2885e6df 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -36,18 +36,18 @@ namespace LR { const std::set lda_xc = { "lda", "pwlda" }; assert(lda_xc.count(this->xc_kernel_)); - const int spinsize = (1 == nspin) ? 1 : 3; - std::vector v2rho2(spinsize * nrxx); + const int n_component = (1 == nspin) ? 1 : 3; + std::vector v2rho2(n_component * nrxx); // read fxc adn add to xc_kernel_components - assert(lr_init_xc_kernel.size() >= spinsize + 1); - for (int is = 0;is < spinsize;++is) + assert(lr_init_xc_kernel.size() >= n_component + 1); + for (int is = 0;is < n_component;++is) { double ef = 0.0; int prenspin = 1; std::vector v2rho2_tmp(nrxx); ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, lr_init_xc_kernel[is + 1], v2rho2_tmp.data(), ucell->nat); - for (int ir = 0;ir < nrxx;++ir) { v2rho2[ir * spinsize + is] = v2rho2_tmp[ir]; } + for (int ir = 0;ir < nrxx;++ir) { v2rho2[ir * n_component + is] = v2rho2_tmp[ir]; } } this->xc_kernel_components_.set_kernel("v2rho2", std::move(v2rho2)); return; From 955d5c6674091d499a15852f36426db320bb8753 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 17:10:37 +0800 Subject: [PATCH 14/23] use something else to replace the for-loop --- source/module_lr/potentials/xc_kernel.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/source/module_lr/potentials/xc_kernel.cpp b/source/module_lr/potentials/xc_kernel.cpp index fb3bd6aa7a..be087f5b11 100644 --- a/source/module_lr/potentials/xc_kernel.cpp +++ b/source/module_lr/potentials/xc_kernel.cpp @@ -144,15 +144,17 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl if (1 == nspin) { + using V3 = ModuleBase::Vector3; // 0. drho - this->grad_kernel_set_.emplace("drho_gs", std::vector>(nrxx)); - for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["drho_gs"][ir] = gradrho[0][ir]; + this->grad_kernel_set_.emplace("drho_gs", gradrho[0]); // 1. $2f^{\rho\sigma}*\nabla\rho$ - this->grad_kernel_set_.emplace("2_v2rhosigma_drho", std::vector>(nrxx)); - for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["2_v2rhosigma_drho"][ir] = gradrho[0][ir] * v2rs.at(ir) * 2.; + this->grad_kernel_set_.emplace("2_v2rhosigma_drho", std::vector(nrxx)); + std::transform(gradrho[0].begin(), gradrho[0].end(), v2rs.begin(), this->grad_kernel_set_["2_v2rhosigma_drho"].begin(), + [](const V3& a, const V3& b) {return a * b * 2.; }); // 2. $4f^{\sigma\sigma}*\nabla\rho$ - this->grad_kernel_set_.emplace("4_v2sigma2_drho", std::vector>(nrxx)); - for (int ir = 0; ir < nrxx; ++ir)this->grad_kernel_set_["4_v2sigma2_drho"][ir] = sigma.at(ir) * v2s2.at(ir) * 4.; + this->grad_kernel_set_.emplace("4_v2sigma2_drho", std::vector(nrxx)); + std::transform(sigma.begin(), sigma.end(), v2s2.begin(), this->grad_kernel_set_["4_v2sigma2_drho"].begin(), + [](const V3& a, const V3& b) {return a * b * 4.; }); } // else if (2 == nspin) // wrong, to be fixed // { From 441e2fa7cb44379a82848e49a5c322bab82eb50b Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Tue, 5 Nov 2024 10:03:44 +0000 Subject: [PATCH 15/23] [pre-commit.ci lite] apply automatic fixes --- source/module_lr/potentials/pot_hxc_lrtd.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index 9a2885e6df..d7b18af527 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -110,7 +110,7 @@ namespace LR { auto& funcs = this->kernel_to_potential_; auto& fxc = this->xc_kernel_components_; - if (xc == XCType::LDA) switch (s) + if (xc == XCType::LDA) { switch (s) { case SpinType::S1: funcs[s] = [this, &fxc](FXC_PARA_TYPE)->void @@ -153,7 +153,7 @@ namespace LR + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); break; } - else if (xc == XCType::GGA || xc == XCType::HYB_GGA) switch (s) + } else if (xc == XCType::GGA || xc == XCType::HYB_GGA) { switch (s) { case SpinType::S1: funcs[s] = [this, &fxc](FXC_PARA_TYPE)->void @@ -204,7 +204,7 @@ namespace LR + std::string(__FILE__) + " line " + std::to_string(__LINE__)); break; } - else + } else { throw std::domain_error("GlobalV::XC_Functional::get_func_type() =" + std::to_string(XC_Functional::get_func_type()) + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); From 1e2ab248c33a59f0fa5dba224902d7db21c16676 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 18:35:58 +0800 Subject: [PATCH 16/23] flatten the map --- source/module_lr/potentials/pot_hxc_lrtd.cpp | 23 ++++++------ source/module_lr/potentials/xc_kernel.cpp | 38 ++++++++++---------- source/module_lr/potentials/xc_kernel.h | 19 +++++----- 3 files changed, 40 insertions(+), 40 deletions(-) diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index d7b18af527..394de59a6b 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -37,7 +37,7 @@ namespace LR const std::set lda_xc = { "lda", "pwlda" }; assert(lda_xc.count(this->xc_kernel_)); const int n_component = (1 == nspin) ? 1 : 3; - std::vector v2rho2(n_component * nrxx); + this->xc_kernel_components_.v2rho2.resize(n_component * nrxx); // read fxc adn add to xc_kernel_components assert(lr_init_xc_kernel.size() >= n_component + 1); for (int is = 0;is < n_component;++is) @@ -47,9 +47,8 @@ namespace LR std::vector v2rho2_tmp(nrxx); ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, lr_init_xc_kernel[is + 1], v2rho2_tmp.data(), ucell->nat); - for (int ir = 0;ir < nrxx;++ir) { v2rho2[ir * n_component + is] = v2rho2_tmp[ir]; } + for (int ir = 0;ir < nrxx;++ir) { xc_kernel_components_.v2rho2[ir * n_component + is] = v2rho2_tmp[ir]; } } - this->xc_kernel_components_.set_kernel("v2rho2", std::move(v2rho2)); return; } @@ -115,7 +114,7 @@ namespace LR case SpinType::S1: funcs[s] = [this, &fxc](FXC_PARA_TYPE)->void { - for (int ir = 0;ir < nrxx;++ir) { v_eff(0, ir) += ModuleBase::e2 * fxc.get_kernel("v2rho2").at(ir) * rho[ir]; } + for (int ir = 0;ir < nrxx;++ir) { v_eff(0, ir) += ModuleBase::e2 * fxc.v2rho2.at(ir) * rho[ir]; } }; break; case SpinType::S2_singlet: @@ -125,7 +124,7 @@ namespace LR { const int irs0 = 3 * ir; const int irs1 = irs0 + 1; - v_eff(0, ir) += ModuleBase::e2 * (fxc.get_kernel("v2rho2").at(irs0) + fxc.get_kernel("v2rho2").at(irs1)) * rho[ir]; + v_eff(0, ir) += ModuleBase::e2 * (fxc.v2rho2.at(irs0) + fxc.v2rho2.at(irs1)) * rho[ir]; } }; break; @@ -136,7 +135,7 @@ namespace LR { const int irs0 = 3 * ir; const int irs1 = irs0 + 1; - v_eff(0, ir) += ModuleBase::e2 * (fxc.get_kernel("v2rho2").at(irs0) - fxc.get_kernel("v2rho2").at(irs1)) * rho[ir]; + v_eff(0, ir) += ModuleBase::e2 * (fxc.v2rho2.at(irs0) - fxc.v2rho2.at(irs1)) * rho[ir]; } }; break; @@ -145,7 +144,7 @@ namespace LR { assert(ispin_op.size() >= 2); const int is = ispin_op[0] + ispin_op[1]; - for (int ir = 0;ir < nrxx;++ir) { v_eff(0, ir) += ModuleBase::e2 * fxc.get_kernel("v2rho2").at(3 * ir + is) * rho[ir]; } + for (int ir = 0;ir < nrxx;++ir) { v_eff(0, ir) += ModuleBase::e2 * fxc.v2rho2.at(3 * ir + is) * rho[ir]; } }; break; default: @@ -180,17 +179,17 @@ namespace LR std::vector> e_drho(nrxx); for (int ir = 0;ir < nrxx;++ir) { - e_drho[ir] = -(fxc.get_grad_kernel("2_v2rhosigma_drho").at(ir) * rho[ir] - + fxc.get_grad_kernel("4_v2sigma2_drho").at(ir) * (fxc.get_grad_kernel("drho_gs").at(ir) * drho.at(ir)) - + drho.at(ir) * fxc.get_kernel("vsigma").at(ir) * 2.); + e_drho[ir] = -(fxc.v2rhosigma_2drho.at(ir) * rho[ir] + + fxc.v2sigma2_4drho.at(ir) * (fxc.drho_gs.at(ir) * drho.at(ir)) + + drho.at(ir) * fxc.vsigma.at(ir) * 2.); } XC_Functional::grad_dot(e_drho.data(), vxc_tmp.data(), this->rho_basis_, this->tpiba_); // 2. $f^{\rho\rho}\rho_1+2f^{\rho\sigma}\nabla\rho\cdot\nabla\rho_1$ for (int ir = 0;ir < nrxx;++ir) { - vxc_tmp[ir] += (fxc.get_kernel("v2rho2").at(ir) * rho[ir] - + fxc.get_grad_kernel("2_v2rhosigma_drho").at(ir) * drho.at(ir)); + vxc_tmp[ir] += (fxc.v2rho2.at(ir) * rho[ir] + + fxc.v2rhosigma_2drho.at(ir) * drho.at(ir)); } BlasConnector::axpy(nrxx, ModuleBase::e2, vxc_tmp.data(), 1, v_eff.c, 1); }; diff --git a/source/module_lr/potentials/xc_kernel.cpp b/source/module_lr/potentials/xc_kernel.cpp index be087f5b11..f6ae538ee1 100644 --- a/source/module_lr/potentials/xc_kernel.cpp +++ b/source/module_lr/potentials/xc_kernel.cpp @@ -96,13 +96,13 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl } // ----------------------------------------------------------------------------------- //==================== XC Kernels (f_xc)============================= - this->kernel_set_.emplace("vrho", std::vector(nspin * nrxx)); - this->kernel_set_.emplace("v2rho2", std::vector(((1 == nspin) ? 1 : 3) * nrxx));//(nrxx* ((1 == nspin) ? 1 : 3)): 00, 01, 11 + this->vrho.resize(nspin * nrxx); + this->v2rho2.resize(((1 == nspin) ? 1 : 3) * nrxx);//(nrxx* ((1 == nspin) ? 1 : 3)): 00, 01, 11 if (is_gga) { - this->kernel_set_.emplace("vsigma", std::vector(((1 == nspin) ? 1 : 3) * nrxx)); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 - this->kernel_set_.emplace("v2rhosigma", std::vector(((1 == nspin) ? 1 : 6) * nrxx)); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 - this->kernel_set_.emplace("v2sigma2", std::vector(((1 == nspin) ? 1 : 6) * nrxx)); //(nrxx* ((1 == nspin) ? 1 : 6)): 00, 01, 02, 11, 12, 22 + this->vsigma.resize(((1 == nspin) ? 1 : 3) * nrxx);//(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 + this->v2rhosigma.resize(((1 == nspin) ? 1 : 6) * nrxx); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 + this->v2sigma2.resize(((1 == nspin) ? 1 : 6) * nrxx); //(nrxx* ((1 == nspin) ? 1 : 6)): 00, 01, 02, 11, 12, 22 } //MetaGGA ... @@ -118,15 +118,13 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl switch (func.info->family) { case XC_FAMILY_LDA: - xc_lda_vxc(&func, nrxx, rho.data(), this->kernel_set_["vrho"].data()); - xc_lda_fxc(&func, nrxx, rho.data(), this->kernel_set_["v2rho2"].data()); + xc_lda_vxc(&func, nrxx, rho.data(), vrho.data()); + xc_lda_fxc(&func, nrxx, rho.data(), v2rho2.data()); break; case XC_FAMILY_GGA: case XC_FAMILY_HYB_GGA: - xc_gga_vxc(&func, nrxx, rho.data(), sigma.data(), - this->kernel_set_["vrho"].data(), this->kernel_set_["vsigma"].data()); - xc_gga_fxc(&func, nrxx, rho.data(), sigma.data(), - this->kernel_set_["v2rho2"].data(), this->kernel_set_["v2rhosigma"].data(), this->kernel_set_["v2sigma2"].data()); + xc_gga_vxc(&func, nrxx, rho.data(), sigma.data(), vrho.data(), vsigma.data()); + xc_gga_fxc(&func, nrxx, rho.data(), sigma.data(), v2rho2.data(), v2rhosigma.data(), v2sigma2.data()); break; default: throw std::domain_error("func.info->family =" + std::to_string(func.info->family) @@ -136,24 +134,24 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl // some formulas for GGA if (func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) { - const std::vector& v2r2 = this->kernel_set_["v2rho2"]; - const std::vector& v2rs = this->kernel_set_["v2rhosigma"]; - const std::vector& v2s2 = this->kernel_set_["v2sigma2"]; - const std::vector& vs = this->kernel_set_["vsigma"]; + const std::vector& v2r2 = this->v2rho2; + const std::vector& v2rs = this->v2rhosigma; + const std::vector& v2s2 = this->v2sigma2; + const std::vector& vs = this->vsigma; const double tpiba2 = tpiba * tpiba; if (1 == nspin) { using V3 = ModuleBase::Vector3; // 0. drho - this->grad_kernel_set_.emplace("drho_gs", gradrho[0]); + this->drho_gs = gradrho[0]; // 1. $2f^{\rho\sigma}*\nabla\rho$ - this->grad_kernel_set_.emplace("2_v2rhosigma_drho", std::vector(nrxx)); - std::transform(gradrho[0].begin(), gradrho[0].end(), v2rs.begin(), this->grad_kernel_set_["2_v2rhosigma_drho"].begin(), + this->v2rhosigma_2drho.resize(nrxx); + std::transform(gradrho[0].begin(), gradrho[0].end(), v2rs.begin(), this->v2rhosigma_2drho.begin(), [](const V3& a, const V3& b) {return a * b * 2.; }); // 2. $4f^{\sigma\sigma}*\nabla\rho$ - this->grad_kernel_set_.emplace("4_v2sigma2_drho", std::vector(nrxx)); - std::transform(sigma.begin(), sigma.end(), v2s2.begin(), this->grad_kernel_set_["4_v2sigma2_drho"].begin(), + this->v2sigma2_4drho.resize(nrxx); + std::transform(sigma.begin(), sigma.end(), v2s2.begin(), v2sigma2_4drho.begin(), [](const V3& a, const V3& b) {return a * b * 4.; }); } // else if (2 == nspin) // wrong, to be fixed diff --git a/source/module_lr/potentials/xc_kernel.h b/source/module_lr/potentials/xc_kernel.h index 097c91be9d..f0022e3a70 100644 --- a/source/module_lr/potentials/xc_kernel.h +++ b/source/module_lr/potentials/xc_kernel.h @@ -15,17 +15,20 @@ namespace LR 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 - // Setters and getters for certain components, such as v2rho2, v2rhosigma, v2sigma2, or other custom combinations. - void set_kernel(const std::string& name, const std::vector& vec) { this->kernel_set_[name] = vec; } - void set_kernel(const std::string& name, const std::vector&& vec) { this->kernel_set_[name] = std::move(vec); } - const std::vector& get_kernel(const std::string& name) { return kernel_set_[name]; } - const std::vector>& get_grad_kernel(const std::string& name) { return grad_kernel_set_[name]; } + // See https://libxc.gitlab.io/manual/libxc-5.1.x/ for the naming convention of the following members. + // std::map> kernel_set_; // [kernel_type][nrxx][nspin] + std::vector vrho; + std::vector vsigma; + std::vector v2rho2; + std::vector v2rhosigma; + std::vector v2sigma2; + // std::map>> grad_kernel_set_;// [kernel_type][nrxx][nspin], intermediate terms for GGA + std::vector> drho_gs; + std::vector> v2rhosigma_2drho; ///< $f^{\rho\sigma}*\nabla\rho *2$ + std::vector> v2sigma2_4drho; ///< $f^{\sigma\sigma}*\nabla\rho *4$ protected: - const ModulePW::PW_Basis& rho_basis_; - std::map> kernel_set_; // [kernel_type][nrxx][nspin] - std::map>> grad_kernel_set_;// [kernel_type][nrxx][nspin], intermediate terms for GGA }; } From 800a02439187fd1705755f94404e49a574be6bf4 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 19:21:30 +0800 Subject: [PATCH 17/23] use std::any_of --- source/module_lr/potentials/xc_kernel.cpp | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/source/module_lr/potentials/xc_kernel.cpp b/source/module_lr/potentials/xc_kernel.cpp index f6ae538ee1..ea04740c99 100644 --- a/source/module_lr/potentials/xc_kernel.cpp +++ b/source/module_lr/potentials/xc_kernel.cpp @@ -37,19 +37,8 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl // ----------------------------------------------------------------------------------- // for GGA - const bool is_gga = [&funcs]() - { - for (xc_func_type& func : funcs) - { - switch (func.info->family) - { - case XC_FAMILY_GGA: - case XC_FAMILY_HYB_GGA: - return true; - } - } - return false; - }(); + const bool is_gga = std::any_of(funcs.begin(), funcs.end(), [](const xc_func_type& f) { f.info->family == XC_FAMILY_GGA || f.info->family == XC_FAMILY_HYB_GGA; }); + std::vector>> gradrho; // \nabla \rho std::vector sigma; // |\nabla\rho|^2 std::vector sgn; // sgn for threshold mask From 3d2009f986ff6f52a50b1c58488b80c1478e6fa9 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 19:24:55 +0800 Subject: [PATCH 18/23] rename new/delete_p2 as new/delete_pointer_dim2 --- source/module_lr/lr_spectrum.cpp | 12 ++++++------ source/module_lr/operator_casida/operator_lr_hxc.cpp | 8 ++++---- source/module_lr/potentials/pot_hxc_lrtd.cpp | 4 ++-- source/module_lr/utils/gint_move.hpp | 4 ++-- source/module_lr/utils/lr_util.h | 4 ++-- source/module_lr/utils/lr_util.hpp | 4 ++-- 6 files changed, 18 insertions(+), 18 deletions(-) diff --git a/source/module_lr/lr_spectrum.cpp b/source/module_lr/lr_spectrum.cpp index 5cd9acfdb3..85616704e5 100644 --- a/source/module_lr/lr_spectrum.cpp +++ b/source/module_lr/lr_spectrum.cpp @@ -53,7 +53,7 @@ void LR::LR_Spectrum::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::new_pointer_dim2(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::oscillator_strength(Grid_Driver& gd, const std::ve ModuleBase::Vector3 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::delete_pointer_dim2(rho_trans, 1); } transition_dipole_[istate] *= (ucell.omega / static_cast(gint->get_ncxyz())); // dv Parallel_Reduce::reduce_all(transition_dipole_[istate].x); @@ -114,8 +114,8 @@ void LR::LR_Spectrum>::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::new_pointer_dim2(rho_trans_real, 1, this->rho_basis.nrxx); + LR_Util::new_pointer_dim2(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>::oscillator_strength(Grid_Driver& gd, ModuleBase::Vector3> rc_complex(rc.x, rc.y, rc.z); transition_dipole_[istate] += rc_complex * std::complex(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::delete_pointer_dim2(rho_trans_real, 1); + LR_Util::delete_pointer_dim2(rho_trans_imag, 1); } transition_dipole_[istate] *= (ucell.omega / static_cast(gint->get_ncxyz())); // dv Parallel_Reduce::reduce_all(transition_dipole_[istate].x); diff --git a/source/module_lr/operator_casida/operator_lr_hxc.cpp b/source/module_lr/operator_casida/operator_lr_hxc.cpp index 31e64a9acf..41bf6d130e 100644 --- a/source/module_lr/operator_casida/operator_lr_hxc.cpp +++ b/source/module_lr/operator_casida/operator_lr_hxc.cpp @@ -61,7 +61,7 @@ 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::new_pointer_dim2(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); @@ -69,7 +69,7 @@ namespace LR // 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); + LR_Util::delete_pointer_dim2(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,7 +103,7 @@ 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::new_pointer_dim2(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); @@ -114,7 +114,7 @@ namespace LR 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::delete_pointer_dim2(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); diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index 394de59a6b..0c74943ad5 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -57,7 +57,7 @@ namespace LR { assert(lr_init_xc_kernel.size() >= 2); double** rho_for_fxc; - LR_Util::new_p2(rho_for_fxc, nspin, nrxx); + LR_Util::new_pointer_dim2(rho_for_fxc, nspin, nrxx); double ef = 0.0; int prenspin = 1; for (int is = 0;is < nspin;++is) @@ -67,7 +67,7 @@ namespace LR rho_for_fxc[is], ucell->nat); } this->xc_kernel_components_.f_xc_libxc(nspin, ucell->omega, ucell->tpiba, rho_for_fxc, chg_gs->rho_core); - LR_Util::delete_p2(rho_for_fxc, nspin); + LR_Util::delete_pointer_dim2(rho_for_fxc, nspin); } else { this->xc_kernel_components_.f_xc_libxc(nspin, ucell->omega, ucell->tpiba, chg_gs->rho, chg_gs->rho_core); } #else diff --git a/source/module_lr/utils/gint_move.hpp b/source/module_lr/utils/gint_move.hpp index 7d21ca6037..ee1c15b129 100644 --- a/source/module_lr/utils/gint_move.hpp +++ b/source/module_lr/utils/gint_move.hpp @@ -12,11 +12,11 @@ using D2 = void(*) (T**, size_t); // template // using D3 = void(*) (T***, size_t, size_t); // template -// D2 d2 = LR_Util::delete_p2; +// D2 d2 = LR_Util::delete_pointer_dim2; // template // D3 d3 = LR_Util::delete_p3; // Change to C++ 11 -D2 d2 = LR_Util::delete_p2; +D2 d2 = LR_Util::delete_pointer_dim2; // D3 d3 = LR_Util::delete_p3; diff --git a/source/module_lr/utils/lr_util.h b/source/module_lr/utils/lr_util.h index 892a42ed7e..79cb5ff95a 100644 --- a/source/module_lr/utils/lr_util.h +++ b/source/module_lr/utils/lr_util.h @@ -51,10 +51,10 @@ namespace LR_Util //====== newers and deleters======== /// @brief delete 2d pointer template - void delete_p2(T** p2, size_t size); + void delete_pointer_dim2(T** p2, size_t size); /// @brief new 2d pointer template - void new_p2(T**& p2, size_t size1, size_t size2); + void new_pointer_dim2(T**& p2, size_t size1, size_t size2); template ct::Tensor newTensor(const ct::TensorShape& shape) { diff --git a/source/module_lr/utils/lr_util.hpp b/source/module_lr/utils/lr_util.hpp index 97a5113eb7..8c5ec1b871 100644 --- a/source/module_lr/utils/lr_util.hpp +++ b/source/module_lr/utils/lr_util.hpp @@ -32,7 +32,7 @@ namespace LR_Util /// @param size1 /// @param size2 template - void new_p2(T**& p2, size_t size1, size_t size2) + void new_pointer_dim2(T**& p2, size_t size1, size_t size2) { p2 = new T * [size1]; for (size_t i = 0; i < size1; ++i) @@ -46,7 +46,7 @@ namespace LR_Util /// @param p2 /// @param size template - void delete_p2(T** p2, size_t size) + void delete_pointer_dim2(T** p2, size_t size) { if (p2 != nullptr) { From 2743d59e24d06096ece91447a4eef4cddf8cb285 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Tue, 5 Nov 2024 20:26:39 +0800 Subject: [PATCH 19/23] a usage of std::inner_product --- source/module_lr/esolver_lrtd_lcao.cpp | 1 + source/module_lr/potentials/xc_kernel.cpp | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 4a37f3b99d..c45f911aa9 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -105,6 +105,7 @@ void LR::ESolver_LR::set_dimension() 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;}(); + this->nbasis = std::inner_product(input.aims_nbasis.begin(), input.aims_nbasis.end(), ucell.atoms, 0, std::plus(), [](const int& a, const Atom& b) { return a * b.na; }); std::cout << "nbasis from aims: " << this->nbasis << std::endl; } } diff --git a/source/module_lr/potentials/xc_kernel.cpp b/source/module_lr/potentials/xc_kernel.cpp index ea04740c99..d883fb7a22 100644 --- a/source/module_lr/potentials/xc_kernel.cpp +++ b/source/module_lr/potentials/xc_kernel.cpp @@ -37,7 +37,7 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl // ----------------------------------------------------------------------------------- // for GGA - const bool is_gga = std::any_of(funcs.begin(), funcs.end(), [](const xc_func_type& f) { f.info->family == XC_FAMILY_GGA || f.info->family == XC_FAMILY_HYB_GGA; }); + const bool is_gga = std::any_of(funcs.begin(), funcs.end(), [](const xc_func_type& f)->bool { return f.info->family == XC_FAMILY_GGA || f.info->family == XC_FAMILY_HYB_GGA; }); std::vector>> gradrho; // \nabla \rho std::vector sigma; // |\nabla\rho|^2 From ba5b170b0066075c42c217d61d1d60ee3b904abd Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Wed, 6 Nov 2024 11:24:42 +0800 Subject: [PATCH 20/23] rename 2order_nested_ptr --- source/module_lr/lr_spectrum.cpp | 12 ++++++------ source/module_lr/operator_casida/operator_lr_hxc.cpp | 8 ++++---- source/module_lr/potentials/pot_hxc_lrtd.cpp | 4 ++-- source/module_lr/utils/gint_move.hpp | 4 ++-- source/module_lr/utils/lr_util.h | 4 ++-- source/module_lr/utils/lr_util.hpp | 4 ++-- 6 files changed, 18 insertions(+), 18 deletions(-) diff --git a/source/module_lr/lr_spectrum.cpp b/source/module_lr/lr_spectrum.cpp index 85616704e5..e38d81609a 100644 --- a/source/module_lr/lr_spectrum.cpp +++ b/source/module_lr/lr_spectrum.cpp @@ -53,7 +53,7 @@ void LR::LR_Spectrum::oscillator_strength(Grid_Driver& gd, const std::ve // 2. transition density double** rho_trans; - LR_Util::new_pointer_dim2(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::oscillator_strength(Grid_Driver& gd, const std::ve ModuleBase::Vector3 rc = rd * ucell.latvec * ucell.lat0; // real coordinate transition_dipole_[istate] += rc * rho_trans[0][ir]; } - LR_Util::delete_pointer_dim2(rho_trans, 1); + LR_Util::_deallocate_2order_nested_ptr(rho_trans, 1); } transition_dipole_[istate] *= (ucell.omega / static_cast(gint->get_ncxyz())); // dv Parallel_Reduce::reduce_all(transition_dipole_[istate].x); @@ -114,8 +114,8 @@ void LR::LR_Spectrum>::oscillator_strength(Grid_Driver& gd, // 2. transition density double** rho_trans_real; double** rho_trans_imag; - LR_Util::new_pointer_dim2(rho_trans_real, 1, this->rho_basis.nrxx); - LR_Util::new_pointer_dim2(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>::oscillator_strength(Grid_Driver& gd, ModuleBase::Vector3> rc_complex(rc.x, rc.y, rc.z); transition_dipole_[istate] += rc_complex * std::complex(rho_trans_real[0][ir], rho_trans_imag[0][ir]); } - LR_Util::delete_pointer_dim2(rho_trans_real, 1); - LR_Util::delete_pointer_dim2(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(gint->get_ncxyz())); // dv Parallel_Reduce::reduce_all(transition_dipole_[istate].x); diff --git a/source/module_lr/operator_casida/operator_lr_hxc.cpp b/source/module_lr/operator_casida/operator_lr_hxc.cpp index 41bf6d130e..9a297bf0be 100644 --- a/source/module_lr/operator_casida/operator_lr_hxc.cpp +++ b/source/module_lr/operator_casida/operator_lr_hxc.cpp @@ -61,7 +61,7 @@ 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_pointer_dim2(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); @@ -69,7 +69,7 @@ namespace LR // 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_pointer_dim2(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); @@ -103,7 +103,7 @@ namespace LR double** rho_trans; const int& nrxx = this->pot.lock()->nrxx; - LR_Util::new_pointer_dim2(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); @@ -114,7 +114,7 @@ namespace LR 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_pointer_dim2(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); diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index 0c74943ad5..0e13ddd1bc 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -57,7 +57,7 @@ namespace LR { assert(lr_init_xc_kernel.size() >= 2); double** rho_for_fxc; - LR_Util::new_pointer_dim2(rho_for_fxc, nspin, nrxx); + LR_Util::_allocate_2order_nested_ptr(rho_for_fxc, nspin, nrxx); double ef = 0.0; int prenspin = 1; for (int is = 0;is < nspin;++is) @@ -67,7 +67,7 @@ namespace LR rho_for_fxc[is], ucell->nat); } this->xc_kernel_components_.f_xc_libxc(nspin, ucell->omega, ucell->tpiba, rho_for_fxc, chg_gs->rho_core); - LR_Util::delete_pointer_dim2(rho_for_fxc, nspin); + LR_Util::_deallocate_2order_nested_ptr(rho_for_fxc, nspin); } else { this->xc_kernel_components_.f_xc_libxc(nspin, ucell->omega, ucell->tpiba, chg_gs->rho, chg_gs->rho_core); } #else diff --git a/source/module_lr/utils/gint_move.hpp b/source/module_lr/utils/gint_move.hpp index ee1c15b129..8b31a09450 100644 --- a/source/module_lr/utils/gint_move.hpp +++ b/source/module_lr/utils/gint_move.hpp @@ -12,11 +12,11 @@ using D2 = void(*) (T**, size_t); // template // using D3 = void(*) (T***, size_t, size_t); // template -// D2 d2 = LR_Util::delete_pointer_dim2; +// D2 d2 = LR_Util::_deallocate_2order_nested_ptr; // template // D3 d3 = LR_Util::delete_p3; // Change to C++ 11 -D2 d2 = LR_Util::delete_pointer_dim2; +D2 d2 = LR_Util::_deallocate_2order_nested_ptr; // D3 d3 = LR_Util::delete_p3; diff --git a/source/module_lr/utils/lr_util.h b/source/module_lr/utils/lr_util.h index 79cb5ff95a..7eed4f3062 100644 --- a/source/module_lr/utils/lr_util.h +++ b/source/module_lr/utils/lr_util.h @@ -51,10 +51,10 @@ namespace LR_Util //====== newers and deleters======== /// @brief delete 2d pointer template - void delete_pointer_dim2(T** p2, size_t size); + void _deallocate_2order_nested_ptr(T** p2, size_t size); /// @brief new 2d pointer template - void new_pointer_dim2(T**& p2, size_t size1, size_t size2); + void _allocate_2order_nested_ptr(T**& p2, size_t size1, size_t size2); template ct::Tensor newTensor(const ct::TensorShape& shape) { diff --git a/source/module_lr/utils/lr_util.hpp b/source/module_lr/utils/lr_util.hpp index 8c5ec1b871..9fda5ccb8d 100644 --- a/source/module_lr/utils/lr_util.hpp +++ b/source/module_lr/utils/lr_util.hpp @@ -32,7 +32,7 @@ namespace LR_Util /// @param size1 /// @param size2 template - void new_pointer_dim2(T**& p2, size_t size1, size_t size2) + void _allocate_2order_nested_ptr(T**& p2, size_t size1, size_t size2) { p2 = new T * [size1]; for (size_t i = 0; i < size1; ++i) @@ -46,7 +46,7 @@ namespace LR_Util /// @param p2 /// @param size template - void delete_pointer_dim2(T** p2, size_t size) + void _deallocate_2order_nested_ptr(T** p2, size_t size) { if (p2 != nullptr) { From 4633457f3508f950da6acfd9ed528132961ac3e2 Mon Sep 17 00:00:00 2001 From: LUNASEA <33978601+maki49@users.noreply.github.com> Date: Wed, 6 Nov 2024 11:26:05 +0800 Subject: [PATCH 21/23] add annotation and minor changes Co-authored-by: kirk0830 <67682086+kirk0830@users.noreply.github.com> --- source/module_lr/esolver_lrtd_lcao.cpp | 8 +++++++- source/module_lr/potentials/pot_hxc_lrtd.cpp | 2 +- source/module_lr/potentials/xc_kernel.cpp | 2 +- 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index c45f911aa9..75560a69dc 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -105,7 +105,13 @@ void LR::ESolver_LR::set_dimension() 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;}(); - this->nbasis = std::inner_product(input.aims_nbasis.begin(), input.aims_nbasis.end(), ucell.atoms, 0, std::plus(), [](const int& a, const Atom& b) { return a * b.na; }); + // 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(), /* iter op1 */ + [](const int& a, const Atom& b) { return a * b.na; }); /* iter op2 */ std::cout << "nbasis from aims: " << this->nbasis << std::endl; } } diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index 0e13ddd1bc..93a22a482f 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -36,7 +36,7 @@ namespace LR { const std::set lda_xc = { "lda", "pwlda" }; assert(lda_xc.count(this->xc_kernel_)); - const int n_component = (1 == nspin) ? 1 : 3; + const int n_component = (1 == nspin) ? 1 : 3; // spin components of fxc: (uu, ud=du, dd) when nspin=2 this->xc_kernel_components_.v2rho2.resize(n_component * nrxx); // read fxc adn add to xc_kernel_components assert(lr_init_xc_kernel.size() >= n_component + 1); diff --git a/source/module_lr/potentials/xc_kernel.cpp b/source/module_lr/potentials/xc_kernel.cpp index d883fb7a22..dfa44d227d 100644 --- a/source/module_lr/potentials/xc_kernel.cpp +++ b/source/module_lr/potentials/xc_kernel.cpp @@ -37,7 +37,7 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl // ----------------------------------------------------------------------------------- // for GGA - const bool is_gga = std::any_of(funcs.begin(), funcs.end(), [](const xc_func_type& f)->bool { return f.info->family == XC_FAMILY_GGA || f.info->family == XC_FAMILY_HYB_GGA; }); + const bool is_gga = std::any_of(funcs.begin(), funcs.end(), [](const xc_func_type& f) { return f.info->family == XC_FAMILY_GGA || f.info->family == XC_FAMILY_HYB_GGA; }); std::vector>> gradrho; // \nabla \rho std::vector sigma; // |\nabla\rho|^2 From 83266daf128d30d7ee273e9a6ddfed204a50f5e7 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Wed, 6 Nov 2024 17:26:18 +0800 Subject: [PATCH 22/23] redesign KernelXC: constructor sets all the members fix compile fix a initialize bug --- source/module_lr/CMakeLists.txt | 9 +- source/module_lr/esolver_lrtd_lcao.cpp | 7 +- .../operator_casida/operator_lr_hxc.cpp | 4 +- source/module_lr/potentials/pot_hxc_lrtd.cpp | 79 +++----------- source/module_lr/potentials/pot_hxc_lrtd.h | 22 ++-- source/module_lr/potentials/xc_kernel.cpp | 103 ++++++++++++++---- source/module_lr/potentials/xc_kernel.h | 37 +++++-- 7 files changed, 137 insertions(+), 124 deletions(-) diff --git a/source/module_lr/CMakeLists.txt b/source/module_lr/CMakeLists.txt index 47a8874aef..a672e05c92 100644 --- a/source/module_lr/CMakeLists.txt +++ b/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/xc_kernel.cpp - ) - endif() + esolver_lrtd_lcao.cpp + potentials/xc_kernel.cpp) add_library( lr diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 75560a69dc..6eecea6b4d 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -104,7 +104,6 @@ void LR::ESolver_LR::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 */ @@ -599,11 +598,11 @@ void LR::ESolver_LR::init_pot(const Charge& chg_gs) { using ST = PotHxcLR::SpinType; case 1: - this->pot[0] = std::make_shared(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, ST::S1, input.lr_init_xc_kernel); + this->pot[0] = std::make_shared(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(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(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_triplet, input.lr_init_xc_kernel); + this->pot[0] = std::make_shared(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(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"); diff --git a/source/module_lr/operator_casida/operator_lr_hxc.cpp b/source/module_lr/operator_casida/operator_lr_hxc.cpp index 9a297bf0be..890b03398a 100644 --- a/source/module_lr/operator_casida/operator_lr_hxc.cpp +++ b/source/module_lr/operator_casida/operator_lr_hxc.cpp @@ -68,7 +68,7 @@ namespace LR // 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); 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) @@ -111,7 +111,7 @@ namespace LR // 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::_deallocate_2order_nested_ptr(rho_trans, 1); diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index 93a22a482f..1508b87301 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -1,5 +1,4 @@ #include "pot_hxc_lrtd.h" -#include "module_elecstate/potentials/pot_base.h" #include "module_parameter/parameter.h" #include "module_elecstate/potentials/H_Hartree_pw.h" #include "module_base/timer.h" @@ -7,76 +6,24 @@ #include #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& ispin_op = { 0,0 } namespace LR { // constructor for exchange-correlation kernel - PotHxcLR::PotHxcLR(const std::string& xc_kernel, const ModulePW::PW_Basis* rho_basis, const UnitCell* ucell, - const Charge* chg_gs/*ground state*/, const Parallel_Grid& pgrid, + PotHxcLR::PotHxcLR(const std::string& xc_kernel, const ModulePW::PW_Basis& rho_basis, const UnitCell& ucell, + const Charge& chg_gs/*ground state*/, const Parallel_Grid& pgrid, const SpinType& st, const std::vector& lr_init_xc_kernel) - :xc_kernel_(xc_kernel), tpiba_(ucell->tpiba), spin_type_(st), - xc_kernel_components_(*rho_basis) + :xc_kernel_(xc_kernel), tpiba_(ucell.tpiba), spin_type_(st), rho_basis_(rho_basis), nrxx_(chg_gs.nrxx), + nspin_(PARAM.inp.nspin == 1 || (PARAM.inp.nspin == 4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z) ? 1 : 2), + pot_hartree_(LR_Util::make_unique(&rho_basis)), + xc_kernel_components_(rho_basis, ucell, chg_gs, pgrid, nspin_, xc_kernel, lr_init_xc_kernel), //call XC_Functional::set_func_type and libxc + xc_type_(XCType(XC_Functional::get_func_type())) { - this->rho_basis_ = rho_basis; - 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(this->rho_basis_); - - const std::set local_xc = { "lda", "pwlda", "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()); - this->set_integral_func(this->spin_type_, this->xc_type_); - - if (lr_init_xc_kernel[0] == "file") - { - const std::set lda_xc = { "lda", "pwlda" }; - assert(lda_xc.count(this->xc_kernel_)); - const int n_component = (1 == nspin) ? 1 : 3; // spin components of fxc: (uu, ud=du, dd) when nspin=2 - this->xc_kernel_components_.v2rho2.resize(n_component * nrxx); - // read fxc adn add to xc_kernel_components - assert(lr_init_xc_kernel.size() >= n_component + 1); - for (int is = 0;is < n_component;++is) - { - double ef = 0.0; - int prenspin = 1; - std::vector v2rho2_tmp(nrxx); - ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, lr_init_xc_kernel[is + 1], - v2rho2_tmp.data(), ucell->nat); - for (int ir = 0;ir < nrxx;++ir) { xc_kernel_components_.v2rho2[ir * n_component + is] = v2rho2_tmp[ir]; } - } - return; - } - -#ifdef USE_LIBXC - if (lr_init_xc_kernel[0] == "from_chg_file") - { - assert(lr_init_xc_kernel.size() >= 2); - double** rho_for_fxc; - LR_Util::_allocate_2order_nested_ptr(rho_for_fxc, nspin, nrxx); - double ef = 0.0; - int prenspin = 1; - for (int is = 0;is < nspin;++is) - { - const std::string file = lr_init_xc_kernel[lr_init_xc_kernel.size() > nspin ? 1 + is : 1]; - ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, file, - rho_for_fxc[is], ucell->nat); - } - this->xc_kernel_components_.f_xc_libxc(nspin, ucell->omega, ucell->tpiba, rho_for_fxc, chg_gs->rho_core); - LR_Util::_deallocate_2order_nested_ptr(rho_for_fxc, nspin); - } - else { this->xc_kernel_components_.f_xc_libxc(nspin, ucell->omega, ucell->tpiba, chg_gs->rho, chg_gs->rho_core); } -#else - ModuleBase::WARNING_QUIT("KernelXC", "to calculate xc-kernel in LR-TDDFT, compile with LIBXC"); -#endif - } + if (std::set({ "lda", "pwlda", "pbe", "hse" }).count(xc_kernel)) { this->set_integral_func(this->spin_type_, this->xc_type_); } } - void PotHxcLR::cal_v_eff(double** rho, const UnitCell* ucell, ModuleBase::matrix& v_eff, const std::vector& ispin_op) + void PotHxcLR::cal_v_eff(double** rho, const UnitCell& ucell, ModuleBase::matrix& v_eff, const std::vector& ispin_op) { ModuleBase::TITLE("PotHxcLR", "cal_v_eff"); ModuleBase::timer::tick("PotHxcLR", "cal_v_eff"); @@ -86,10 +33,10 @@ namespace LR switch (this->spin_type_) { case SpinType::S1: case SpinType::S2_updown: - v_eff += elecstate::H_Hartree_pw::v_hartree(*ucell, const_cast(this->rho_basis_), 1, rho); + v_eff += elecstate::H_Hartree_pw::v_hartree(ucell, const_cast(&this->rho_basis_), 1, rho); break; case SpinType::S2_singlet: - v_eff += 2 * elecstate::H_Hartree_pw::v_hartree(*ucell, const_cast(this->rho_basis_), 1, rho); + v_eff += 2 * elecstate::H_Hartree_pw::v_hartree(ucell, const_cast(&this->rho_basis_), 1, rho); break; default: break; @@ -171,7 +118,7 @@ namespace LR // std::cout << std::endl;}; std::vector> drho(nrxx); // transition density gradient - LR_Util::grad(rho, drho.data(), *(this->rho_basis_), this->tpiba_); + LR_Util::grad(rho, drho.data(), this->rho_basis_, this->tpiba_); std::vector vxc_tmp(nrxx, 0.0); @@ -183,7 +130,7 @@ namespace LR + fxc.v2sigma2_4drho.at(ir) * (fxc.drho_gs.at(ir) * drho.at(ir)) + drho.at(ir) * fxc.vsigma.at(ir) * 2.); } - XC_Functional::grad_dot(e_drho.data(), vxc_tmp.data(), this->rho_basis_, this->tpiba_); + XC_Functional::grad_dot(e_drho.data(), vxc_tmp.data(), &this->rho_basis_, this->tpiba_); // 2. $f^{\rho\rho}\rho_1+2f^{\rho\sigma}\nabla\rho\cdot\nabla\rho_1$ for (int ir = 0;ir < nrxx;++ir) diff --git a/source/module_lr/potentials/pot_hxc_lrtd.h b/source/module_lr/potentials/pot_hxc_lrtd.h index 59d8d71dbe..f7b78355d7 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.h +++ b/source/module_lr/potentials/pot_hxc_lrtd.h @@ -1,14 +1,12 @@ #pragma once -#include "module_elecstate/potentials/pot_base.h" #include "module_elecstate/potentials/H_Hartree_pw.h" #include "xc_kernel.h" #include #include -class Parallel_Grid; namespace LR { - class PotHxcLR : public elecstate::PotBase + class PotHxcLR { public: /// S1: K^Hartree + K^xc @@ -16,23 +14,25 @@ namespace LR /// S2_triplet: K^xc_{upup} - K^xc_{updown} /// S2_updown: K^Hartree + (K^xc_{upup}, K^xc_{updown}, K^xc_{downup} or K^xc_{downdown}), according to `ispin_op` (for spin-polarized systems) enum SpinType { S1 = 0, S2_singlet = 1, S2_triplet = 2, S2_updown = 3 }; + /// XCType here is to determin the method of integration from kernel to potential, not the way calculating the kernel enum XCType { None = 0, LDA = 1, GGA = 2, HYB_GGA = 4 }; /// constructor for exchange-correlation kernel - PotHxcLR(const std::string& xc_kernel, const ModulePW::PW_Basis* rho_basis, - const UnitCell* ucell, const Charge* chg_gs/*ground state*/, const Parallel_Grid& pgrid, + PotHxcLR(const std::string& xc_kernel, const ModulePW::PW_Basis& rho_basis, + const UnitCell& ucell, const Charge& chg_gs/*ground state*/, const Parallel_Grid& pgrid, const SpinType& st = SpinType::S1, const std::vector& lr_init_xc_kernel = { "from_chg_groundstate" }); ~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& ispin_op = { 0,0 }); - int nrxx; + void cal_v_eff(double** rho, const UnitCell& ucell, ModuleBase::matrix& v_eff, const std::vector& ispin_op = { 0,0 }); + const int& nrxx = nrxx_; private: - int nspin; - std::unique_ptr pot_hartree; + const ModulePW::PW_Basis& rho_basis_; + const int nspin_ = 1; + const int nrxx_ = 1; + std::unique_ptr pot_hartree_; /// different components of local and semi-local xc kernels: /// LDA: v2rho2 /// GGA: v2rho2, v2rhosigma, v2sigma2 /// meta-GGA: v2rho2, v2rhosigma, v2sigma2, v2rholap, v2rhotau, v2sigmalap, v2sigmatau, v2laptau, v2lap2, v2tau2 - KernelXC xc_kernel_components_; + const KernelXC xc_kernel_components_; const std::string xc_kernel_; const double& tpiba_; const SpinType spin_type_ = SpinType::S1; diff --git a/source/module_lr/potentials/xc_kernel.cpp b/source/module_lr/potentials/xc_kernel.cpp index dfa44d227d..021c77ba3b 100644 --- a/source/module_lr/potentials/xc_kernel.cpp +++ b/source/module_lr/potentials/xc_kernel.cpp @@ -4,10 +4,69 @@ #include "module_base/timer.h" #include "module_lr/utils/lr_util.h" #include "module_lr/utils/lr_util_xc.hpp" +#include +#include "module_io/cube_io.h" #ifdef USE_LIBXC #include #include "module_hamilt_general/module_xc/xc_functional_libxc.h" +#endif + +LR::KernelXC::KernelXC(const ModulePW::PW_Basis& rho_basis, + const UnitCell& ucell, + const Charge& chg_gs, + const Parallel_Grid& pgrid, + const int& nspin, + const std::string& kernel_name, + const std::vector& lr_init_xc_kernel) :rho_basis_(rho_basis) +{ + if (!std::set({ "lda", "pwlda", "pbe", "hse" }).count(kernel_name)) { return; } + XC_Functional::set_xc_type(kernel_name); // for hse, (1-alpha) and omega are set here + const int& nrxx = rho_basis.nrxx; + if (lr_init_xc_kernel[0] == "file") + { + const std::set lda_xc = { "lda", "pwlda" }; + assert(lda_xc.count(kernel_name)); + const int n_component = (1 == nspin) ? 1 : 3; // spin components of fxc: (uu, ud=du, dd) when nspin=2 + this->v2rho2_.resize(n_component * nrxx); + // read fxc adn add to xc_kernel_components + assert(lr_init_xc_kernel.size() >= n_component + 1); + for (int is = 0;is < n_component;++is) + { + double ef = 0.0; + int prenspin = 1; + std::vector v2rho2_tmp(nrxx); + ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, lr_init_xc_kernel[is + 1], + v2rho2_tmp.data(), ucell.nat); + for (int ir = 0;ir < nrxx;++ir) { this->v2rho2_[ir * n_component + is] = v2rho2_tmp[ir]; } + } + return; + } + +#ifdef USE_LIBXC + if (lr_init_xc_kernel[0] == "from_chg_file") + { + assert(lr_init_xc_kernel.size() >= 2); + double** rho_for_fxc; + LR_Util::_allocate_2order_nested_ptr(rho_for_fxc, nspin, nrxx); + double ef = 0.0; + int prenspin = 1; + for (int is = 0;is < nspin;++is) + { + const std::string file = lr_init_xc_kernel[lr_init_xc_kernel.size() > nspin ? 1 + is : 1]; + ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, file, + rho_for_fxc[is], ucell.nat); + } + this->f_xc_libxc(nspin, ucell.omega, ucell.tpiba, rho_for_fxc, chg_gs.rho_core); + LR_Util::_deallocate_2order_nested_ptr(rho_for_fxc, nspin); + } + else { this->f_xc_libxc(nspin, ucell.omega, ucell.tpiba, chg_gs.rho, chg_gs.rho_core); } +#else + ModuleBase::WARNING_QUIT("KernelXC", "to calculate xc-kernel in LR-TDDFT, compile with LIBXC"); +#endif +} + +#ifdef USE_LIBXC 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"); @@ -85,13 +144,13 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl } // ----------------------------------------------------------------------------------- //==================== XC Kernels (f_xc)============================= - this->vrho.resize(nspin * nrxx); - this->v2rho2.resize(((1 == nspin) ? 1 : 3) * nrxx);//(nrxx* ((1 == nspin) ? 1 : 3)): 00, 01, 11 + this->vrho_.resize(nspin * nrxx); + this->v2rho2_.resize(((1 == nspin) ? 1 : 3) * nrxx);//(nrxx* ((1 == nspin) ? 1 : 3)): 00, 01, 11 if (is_gga) { - this->vsigma.resize(((1 == nspin) ? 1 : 3) * nrxx);//(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 - this->v2rhosigma.resize(((1 == nspin) ? 1 : 6) * nrxx); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 - this->v2sigma2.resize(((1 == nspin) ? 1 : 6) * nrxx); //(nrxx* ((1 == nspin) ? 1 : 6)): 00, 01, 02, 11, 12, 22 + this->vsigma_.resize(((1 == nspin) ? 1 : 3) * nrxx);//(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 + this->v2rhosigma_.resize(((1 == nspin) ? 1 : 6) * nrxx); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 + this->v2sigma2_.resize(((1 == nspin) ? 1 : 6) * nrxx); //(nrxx* ((1 == nspin) ? 1 : 6)): 00, 01, 02, 11, 12, 22 } //MetaGGA ... @@ -107,13 +166,13 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl switch (func.info->family) { case XC_FAMILY_LDA: - xc_lda_vxc(&func, nrxx, rho.data(), vrho.data()); - xc_lda_fxc(&func, nrxx, rho.data(), v2rho2.data()); + xc_lda_vxc(&func, nrxx, rho.data(), vrho_.data()); + xc_lda_fxc(&func, nrxx, rho.data(), v2rho2_.data()); break; case XC_FAMILY_GGA: case XC_FAMILY_HYB_GGA: - xc_gga_vxc(&func, nrxx, rho.data(), sigma.data(), vrho.data(), vsigma.data()); - xc_gga_fxc(&func, nrxx, rho.data(), sigma.data(), v2rho2.data(), v2rhosigma.data(), v2sigma2.data()); + xc_gga_vxc(&func, nrxx, rho.data(), sigma.data(), vrho_.data(), vsigma_.data()); + xc_gga_fxc(&func, nrxx, rho.data(), sigma.data(), v2rho2_.data(), v2rhosigma_.data(), v2sigma2_.data()); break; default: throw std::domain_error("func.info->family =" + std::to_string(func.info->family) @@ -123,24 +182,24 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl // some formulas for GGA if (func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) { - const std::vector& v2r2 = this->v2rho2; - const std::vector& v2rs = this->v2rhosigma; - const std::vector& v2s2 = this->v2sigma2; - const std::vector& vs = this->vsigma; + const std::vector& v2r2 = this->v2rho2_; + const std::vector& v2rs = this->v2rhosigma_; + const std::vector& v2s2 = this->v2sigma2_; + const std::vector& vs = this->vsigma_; const double tpiba2 = tpiba * tpiba; if (1 == nspin) { using V3 = ModuleBase::Vector3; // 0. drho - this->drho_gs = gradrho[0]; + this->drho_gs_ = gradrho[0]; // 1. $2f^{\rho\sigma}*\nabla\rho$ - this->v2rhosigma_2drho.resize(nrxx); - std::transform(gradrho[0].begin(), gradrho[0].end(), v2rs.begin(), this->v2rhosigma_2drho.begin(), + this->v2rhosigma_2drho_.resize(nrxx); + std::transform(gradrho[0].begin(), gradrho[0].end(), v2rs.begin(), this->v2rhosigma_2drho_.begin(), [](const V3& a, const V3& b) {return a * b * 2.; }); // 2. $4f^{\sigma\sigma}*\nabla\rho$ - this->v2sigma2_4drho.resize(nrxx); - std::transform(sigma.begin(), sigma.end(), v2s2.begin(), v2sigma2_4drho.begin(), + this->v2sigma2_4drho_.resize(nrxx); + std::transform(sigma.begin(), sigma.end(), v2s2.begin(), v2sigma2_4drho_.begin(), [](const V3& a, const V3& b) {return a * b * 4.; }); } // else if (2 == nspin) // wrong, to be fixed @@ -158,7 +217,7 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl // + gradrho[0][ir] * v2rs.at(ir * 6 + 4) * 2.0; //down-down // } // for (int isig = 0;isig < 3;++isig) - // XC_Functional::grad_dot(v2rhosigma_gdrho_r.data() + isig * nrxx, div_v2rhosigma_gdrho_r.data() + isig * nrxx, chg_gs->rhopw, tpiba); + // XC_Functional::grad_dot(v2rhosigma_gdrho_r.data() + isig * nrxx, div_v2rhosigma_gdrho_r.data() + isig * nrxx, chg_gs.rhopw, tpiba); // // 2. $\nabla^2(f^{\sigma\sigma}*\sigma)$ // std::vector v2sigma2_sigma_r(3 * nrxx); // for (int ir = 0; ir < nrxx; ++ir) @@ -174,7 +233,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::lapl(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 lap_vsigma(3 * nrxx); // for (int ir = 0;ir < nrxx;++ir) @@ -184,8 +243,8 @@ 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::lapl(lap_vsigma.data() + isig * nrxx, lap_vsigma.data() + isig * nrxx, *(chg_gs->rhopw), tpiba2); - // // add to v2rho2 + // 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); // BlasConnector::axpy(3 * nrxx, 1.0, v2sigma2_sigma_r.data(), 1, to_mul_rho_.data(), 1); diff --git a/source/module_lr/potentials/xc_kernel.h b/source/module_lr/potentials/xc_kernel.h index f0022e3a70..0db7f3d616 100644 --- a/source/module_lr/potentials/xc_kernel.h +++ b/source/module_lr/potentials/xc_kernel.h @@ -1,33 +1,46 @@ #pragma once #include "module_basis/module_pw/pw_basis.h" - +#include "module_cell/unitcell.h" +#include "module_hamilt_pw/hamilt_pwdft/parallel_grid.h" +#include "module_elecstate/module_charge/charge.h" +#define CREF(x) const std::vector& x = x##_; +#define CREF3(x) const std::vector>& x = x##_; namespace LR { /// @brief Calculate the exchange-correlation (XC) kernel ($f_{xc}=\delta^2E_xc/\delta\rho^2$) and store its components. class KernelXC { public: - KernelXC(const ModulePW::PW_Basis& rho_basis) : rho_basis_(rho_basis) {}; + KernelXC(const ModulePW::PW_Basis& rho_basis, + const UnitCell& ucell, + const Charge& chg_gs, + const Parallel_Grid& pgrid, + const int& nspin, + const std::string& kernel_name, + const std::vector& lr_init_xc_kernel); ~KernelXC() {}; + // const references + CREF(vrho);CREF(vsigma); CREF(v2rho2); CREF(v2rhosigma); CREF(v2sigma2); + CREF3(drho_gs); CREF3(v2rhosigma_2drho); CREF3(v2sigma2_4drho); + + private: #ifdef USE_LIBXC /// @brief Calculate the XC kernel using libxc. 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 - // See https://libxc.gitlab.io/manual/libxc-5.1.x/ for the naming convention of the following members. // std::map> kernel_set_; // [kernel_type][nrxx][nspin] - std::vector vrho; - std::vector vsigma; - std::vector v2rho2; - std::vector v2rhosigma; - std::vector v2sigma2; + std::vector vrho_; + std::vector vsigma_; + std::vector v2rho2_; + std::vector v2rhosigma_; + std::vector v2sigma2_; // std::map>> grad_kernel_set_;// [kernel_type][nrxx][nspin], intermediate terms for GGA - std::vector> drho_gs; - std::vector> v2rhosigma_2drho; ///< $f^{\rho\sigma}*\nabla\rho *2$ - std::vector> v2sigma2_4drho; ///< $f^{\sigma\sigma}*\nabla\rho *4$ + std::vector> drho_gs_; + std::vector> v2rhosigma_2drho_; ///< $f^{\rho\sigma}*\nabla\rho *2$ + std::vector> v2sigma2_4drho_; ///< $f^{\sigma\sigma}*\nabla\rho *4$ - protected: const ModulePW::PW_Basis& rho_basis_; }; } From 11ab3d63f60b84a32b340771ecc0d3a57e8d546e Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Fri, 8 Nov 2024 03:20:57 +0800 Subject: [PATCH 23/23] rename and update doc --- docs/advanced/input_files/input-main.md | 9 +++++---- source/module_io/read_input_item_tddft.cpp | 4 ++-- source/module_io/test/read_input_ptest.cpp | 2 +- source/module_lr/potentials/pot_hxc_lrtd.h | 2 +- source/module_lr/potentials/xc_kernel.cpp | 2 +- 5 files changed, 10 insertions(+), 9 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index d2f1c30211..8e15de5b6b 100644 --- a/docs/advanced/input_files/input-main.md +++ b/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) @@ -3947,10 +3948,10 @@ Currently supported: `RPA`, `LDA`, `PBE`, `HSE`, `HF`. - **Type**: String - **Description**: The method to initalize the xc kernel. - - "from_chg_groundstate": 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 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. - - "from_chg_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**: "from_chg_groundstate" + - "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 diff --git a/source/module_io/read_input_item_tddft.cpp b/source/module_io/read_input_item_tddft.cpp index ce71273dc7..bfd0542729 100644 --- a/source/module_io/read_input_item_tddft.cpp +++ b/source/module_io/read_input_item_tddft.cpp @@ -318,9 +318,9 @@ void ReadInput::item_lr_tddft() 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("from_chg_groundstate"); } + 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(), "from_chg_groundstate"); + sync_stringvec(input.lr_init_xc_kernel, para.input.lr_init_xc_kernel.size(), "default"); this->add_item(item); } { diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index 9cf702c2ac..7e17e2ae37 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -420,7 +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], "from_chg_groundstate"); + 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); diff --git a/source/module_lr/potentials/pot_hxc_lrtd.h b/source/module_lr/potentials/pot_hxc_lrtd.h index f7b78355d7..c165999ffb 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.h +++ b/source/module_lr/potentials/pot_hxc_lrtd.h @@ -19,7 +19,7 @@ namespace LR /// constructor for exchange-correlation kernel PotHxcLR(const std::string& xc_kernel, const ModulePW::PW_Basis& rho_basis, const UnitCell& ucell, const Charge& chg_gs/*ground state*/, const Parallel_Grid& pgrid, - const SpinType& st = SpinType::S1, const std::vector& lr_init_xc_kernel = { "from_chg_groundstate" }); + const SpinType& st = SpinType::S1, const std::vector& lr_init_xc_kernel = { "default" }); ~PotHxcLR() {} void cal_v_eff(double** rho, const UnitCell& ucell, ModuleBase::matrix& v_eff, const std::vector& ispin_op = { 0,0 }); const int& nrxx = nrxx_; diff --git a/source/module_lr/potentials/xc_kernel.cpp b/source/module_lr/potentials/xc_kernel.cpp index 021c77ba3b..f5fa4f8fc8 100644 --- a/source/module_lr/potentials/xc_kernel.cpp +++ b/source/module_lr/potentials/xc_kernel.cpp @@ -44,7 +44,7 @@ LR::KernelXC::KernelXC(const ModulePW::PW_Basis& rho_basis, } #ifdef USE_LIBXC - if (lr_init_xc_kernel[0] == "from_chg_file") + if (lr_init_xc_kernel[0] == "from_charge_file") { assert(lr_init_xc_kernel.size() >= 2); double** rho_for_fxc;