Skip to content

Commit

Permalink
redesign KernelXC: constructor sets all the members
Browse files Browse the repository at this point in the history 
fix compile

fix a initialize bug
  • Loading branch information
@maki49
maki49 committed Nov 7, 2024
1 parent 69bedd8 commit ac9b788
Show file tree
Hide file tree
Showing 7 changed files with 137 additions and 124 deletions.
9 changes: 2 additions & 7 deletions source/module_lr/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
7 changes: 3 additions & 4 deletions source/module_lr/esolver_lrtd_lcao.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,6 @@ void LR::ESolver_LR<T, TR>::set_dimension()
GlobalV::ofs_running << "number of excited states to be solved: " << this->nstates << std::endl;
if (input.ri_hartree_benchmark == "aims" && !input.aims_nbasis.empty())
{
this->nbasis = [&]() -> int { int nbas = 0; for (int it = 0;it < ucell.ntype;++it) { nbas += ucell.atoms[it].na * input.aims_nbasis[it]; };return nbas;}();
// calculate total number of basis funcs, see https://en.cppreference.com/w/cpp/algorithm/inner_product
this->nbasis = std::inner_product(input.aims_nbasis.begin(), /* iterator1.begin */
input.aims_nbasis.end(), /* iterator1.end */
Expand Down Expand Up @@ -610,11 +609,11 @@ void LR::ESolver_LR<T, TR>::init_pot(const Charge& chg_gs)
{
using ST = PotHxcLR::SpinType;
case 1:
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, ST::S1, input.lr_init_xc_kernel);
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, ucell, chg_gs, GlobalC::Pgrid, ST::S1, input.lr_init_xc_kernel);
break;
case 2:
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_singlet, input.lr_init_xc_kernel);
this->pot[1] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_triplet, input.lr_init_xc_kernel);
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, ucell, chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_singlet, input.lr_init_xc_kernel);
this->pot[1] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, ucell, chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_triplet, input.lr_init_xc_kernel);
break;
default:
throw std::invalid_argument("ESolver_LR: nspin must be 1 or 2");
Expand Down
4 changes: 2 additions & 2 deletions source/module_lr/operator_casida/operator_lr_hxc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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);
Expand Down
79 changes: 13 additions & 66 deletions source/module_lr/potentials/pot_hxc_lrtd.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,82 +1,29 @@
#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"
#include "module_hamilt_general/module_xc/xc_functional.h"
#include <set>
#include "module_lr/utils/lr_util.h"
#include "module_lr/utils/lr_util_xc.hpp"
#include "module_io/cube_io.h"
#include "module_hamilt_pw/hamilt_pwdft/global.h" // tmp, for pgrid
#define FXC_PARA_TYPE const double* const rho, ModuleBase::matrix& v_eff, const std::vector<int>& ispin_op = { 0,0 }
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<std::string>& 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<elecstate::PotHartree>(&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<elecstate::PotHartree>(this->rho_basis_);

const std::set<std::string> 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<std::string> 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<double> 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<std::string>({ "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<int>& ispin_op)
void PotHxcLR::cal_v_eff(double** rho, const UnitCell& ucell, ModuleBase::matrix& v_eff, const std::vector<int>& ispin_op)
{
ModuleBase::TITLE("PotHxcLR", "cal_v_eff");
ModuleBase::timer::tick("PotHxcLR", "cal_v_eff");
Expand All @@ -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<ModulePW::PW_Basis*>(this->rho_basis_), 1, rho);
v_eff += elecstate::H_Hartree_pw::v_hartree(ucell, const_cast<ModulePW::PW_Basis*>(&this->rho_basis_), 1, rho);
break;
case SpinType::S2_singlet:
v_eff += 2 * elecstate::H_Hartree_pw::v_hartree(*ucell, const_cast<ModulePW::PW_Basis*>(this->rho_basis_), 1, rho);
v_eff += 2 * elecstate::H_Hartree_pw::v_hartree(ucell, const_cast<ModulePW::PW_Basis*>(&this->rho_basis_), 1, rho);
break;
default:
break;
Expand Down Expand Up @@ -171,7 +118,7 @@ namespace LR
// std::cout << std::endl;};

std::vector<ModuleBase::Vector3<double>> 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<double> vxc_tmp(nrxx, 0.0);

Expand All @@ -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)
Expand Down
22 changes: 11 additions & 11 deletions source/module_lr/potentials/pot_hxc_lrtd.h
View file
Original file line number Diff line number Diff line change
@@ -1,38 +1,38 @@
#pragma once
#include "module_elecstate/potentials/pot_base.h"
#include "module_elecstate/potentials/H_Hartree_pw.h"
#include "xc_kernel.h"
#include <unordered_map>
#include <memory>

class Parallel_Grid;
namespace LR
{
class PotHxcLR : public elecstate::PotBase
class PotHxcLR
{
public:
/// 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 };
/// 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<std::string>& 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<int>& ispin_op = { 0,0 });
int nrxx;
void cal_v_eff(double** rho, const UnitCell& ucell, ModuleBase::matrix& v_eff, const std::vector<int>& ispin_op = { 0,0 });
const int& nrxx = nrxx_;
private:
int nspin;
std::unique_ptr<elecstate::PotHartree> pot_hartree;
const ModulePW::PW_Basis& rho_basis_;
const int nspin_ = 1;
const int nrxx_ = 1;
std::unique_ptr<elecstate::PotHartree> 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;
Expand Down
Loading

0 comments on commit ac9b788

Please sign in to comment.