Skip to content

Commit 97f2591

Browse files
maki49maki49
committed
redesign KernelXC: constructor sets all the members
fix compile fix a initialize bug
1 parent 69bedd8 commit 97f2591

File tree

7 files changed

+137
-123
lines changed

7 files changed

+137
-123
lines changed

‎source/module_lr/CMakeLists.txt‎

Lines changed: 2 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -17,13 +17,8 @@ if(ENABLE_LCAO)
1717
potentials/pot_hxc_lrtd.cpp
1818
lr_spectrum.cpp
1919
hamilt_casida.cpp
20-
esolver_lrtd_lcao.cpp)
21-
22-
if(ENABLE_LIBXC)
23-
list(APPEND objects
24-
potentials/xc_kernel.cpp
25-
)
26-
endif()
20+
esolver_lrtd_lcao.cpp
21+
potentials/xc_kernel.cpp)
2722

2823
add_library(
2924
lr

‎source/module_lr/esolver_lrtd_lcao.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -610,11 +610,11 @@ void LR::ESolver_LR<T, TR>::init_pot(const Charge& chg_gs)
610610
{
611611
using ST = PotHxcLR::SpinType;
612612
case 1:
613-
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, ST::S1, input.lr_init_xc_kernel);
613+
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, ucell, chg_gs, GlobalC::Pgrid, ST::S1, input.lr_init_xc_kernel);
614614
break;
615615
case 2:
616-
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);
617-
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);
616+
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);
617+
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);
618618
break;
619619
default:
620620
throw std::invalid_argument("ESolver_LR: nspin must be 1 or 2");

‎source/module_lr/operator_casida/operator_lr_hxc.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -68,7 +68,7 @@ namespace LR
6868

6969
// 3. v_hxc = f_hxc * rho_trans
7070
ModuleBase::matrix vr_hxc(1, nrxx); //grid
71-
this->pot.lock()->cal_v_eff(rho_trans, &GlobalC::ucell, vr_hxc, ispin_ks);
71+
this->pot.lock()->cal_v_eff(rho_trans, GlobalC::ucell, vr_hxc, ispin_ks);
7272
LR_Util::_deallocate_2order_nested_ptr(rho_trans, 1);
7373

7474
// 4. V^{Hxc}_{\mu,\nu}=\int{dr} \phi_\mu(r) v_{Hxc}(r) \phi_\mu(r)
@@ -111,7 +111,7 @@ namespace LR
111111

112112
// 3. v_hxc = f_hxc * rho_trans
113113
ModuleBase::matrix vr_hxc(1, nrxx); //grid
114-
this->pot.lock()->cal_v_eff(rho_trans, &GlobalC::ucell, vr_hxc, ispin_ks);
114+
this->pot.lock()->cal_v_eff(rho_trans, GlobalC::ucell, vr_hxc, ispin_ks);
115115
// print_grid_nonzero(vr_hxc.c, this->poticab->nrxx, 10, "vr_hxc");
116116

117117
LR_Util::_deallocate_2order_nested_ptr(rho_trans, 1);

‎source/module_lr/potentials/pot_hxc_lrtd.cpp‎

Lines changed: 13 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -1,82 +1,29 @@
11
#include "pot_hxc_lrtd.h"
2-
#include "module_elecstate/potentials/pot_base.h"
32
#include "module_parameter/parameter.h"
43
#include "module_elecstate/potentials/H_Hartree_pw.h"
54
#include "module_base/timer.h"
65
#include "module_hamilt_general/module_xc/xc_functional.h"
76
#include <set>
87
#include "module_lr/utils/lr_util.h"
98
#include "module_lr/utils/lr_util_xc.hpp"
10-
#include "module_io/cube_io.h"
119
#include "module_hamilt_pw/hamilt_pwdft/global.h" // tmp, for pgrid
1210
#define FXC_PARA_TYPE const double* const rho, ModuleBase::matrix& v_eff, const std::vector<int>& ispin_op = { 0,0 }
1311
namespace LR
1412
{
1513
// constructor for exchange-correlation kernel
16-
PotHxcLR::PotHxcLR(const std::string& xc_kernel, const ModulePW::PW_Basis* rho_basis, const UnitCell* ucell,
17-
const Charge* chg_gs/*ground state*/, const Parallel_Grid& pgrid,
14+
PotHxcLR::PotHxcLR(const std::string& xc_kernel, const ModulePW::PW_Basis& rho_basis, const UnitCell& ucell,
15+
const Charge& chg_gs/*ground state*/, const Parallel_Grid& pgrid,
1816
const SpinType& st, const std::vector<std::string>& lr_init_xc_kernel)
19-
:xc_kernel_(xc_kernel), tpiba_(ucell->tpiba), spin_type_(st),
20-
xc_kernel_components_(*rho_basis)
17+
:xc_kernel_(xc_kernel), tpiba_(ucell.tpiba), spin_type_(st), rho_basis_(rho_basis), nrxx_(chg_gs.nrxx),
18+
nspin_(PARAM.inp.nspin == 1 || (PARAM.inp.nspin == 4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z) ? 1 : 2),
19+
pot_hartree_(LR_Util::make_unique<elecstate::PotHartree>(&rho_basis)),
20+
xc_kernel_components_(rho_basis, ucell, chg_gs, pgrid, nspin_, xc_kernel, lr_init_xc_kernel), //call XC_Functional::set_func_type and libxc
21+
xc_type_(XCType(XC_Functional::get_func_type()))
2122
{
22-
this->rho_basis_ = rho_basis;
23-
this->nrxx = chg_gs->nrxx;
24-
this->nspin = (PARAM.inp.nspin == 1 || (PARAM.inp.nspin == 4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z)) ? 1 : 2;
25-
26-
this->pot_hartree = LR_Util::make_unique<elecstate::PotHartree>(this->rho_basis_);
27-
28-
const std::set<std::string> local_xc = { "lda", "pwlda", "pbe", "hse" };
29-
if (local_xc.find(this->xc_kernel_) != local_xc.end())
30-
{
31-
XC_Functional::set_xc_type(this->xc_kernel_); // for hse, (1-alpha) and omega are set here
32-
this->xc_type_ = XCType(XC_Functional::get_func_type());
33-
this->set_integral_func(this->spin_type_, this->xc_type_);
34-
35-
if (lr_init_xc_kernel[0] == "file")
36-
{
37-
const std::set<std::string> lda_xc = { "lda", "pwlda" };
38-
assert(lda_xc.count(this->xc_kernel_));
39-
const int n_component = (1 == nspin) ? 1 : 3; // spin components of fxc: (uu, ud=du, dd) when nspin=2
40-
this->xc_kernel_components_.v2rho2.resize(n_component * nrxx);
41-
// read fxc adn add to xc_kernel_components
42-
assert(lr_init_xc_kernel.size() >= n_component + 1);
43-
for (int is = 0;is < n_component;++is)
44-
{
45-
double ef = 0.0;
46-
int prenspin = 1;
47-
std::vector<double> v2rho2_tmp(nrxx);
48-
ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, lr_init_xc_kernel[is + 1],
49-
v2rho2_tmp.data(), ucell->nat);
50-
for (int ir = 0;ir < nrxx;++ir) { xc_kernel_components_.v2rho2[ir * n_component + is] = v2rho2_tmp[ir]; }
51-
}
52-
return;
53-
}
54-
55-
#ifdef USE_LIBXC
56-
if (lr_init_xc_kernel[0] == "from_chg_file")
57-
{
58-
assert(lr_init_xc_kernel.size() >= 2);
59-
double** rho_for_fxc;
60-
LR_Util::_allocate_2order_nested_ptr(rho_for_fxc, nspin, nrxx);
61-
double ef = 0.0;
62-
int prenspin = 1;
63-
for (int is = 0;is < nspin;++is)
64-
{
65-
const std::string file = lr_init_xc_kernel[lr_init_xc_kernel.size() > nspin ? 1 + is : 1];
66-
ModuleIO::read_vdata_palgrid(pgrid, GlobalV::MY_RANK, GlobalV::ofs_running, file,
67-
rho_for_fxc[is], ucell->nat);
68-
}
69-
this->xc_kernel_components_.f_xc_libxc(nspin, ucell->omega, ucell->tpiba, rho_for_fxc, chg_gs->rho_core);
70-
LR_Util::_deallocate_2order_nested_ptr(rho_for_fxc, nspin);
71-
}
72-
else { this->xc_kernel_components_.f_xc_libxc(nspin, ucell->omega, ucell->tpiba, chg_gs->rho, chg_gs->rho_core); }
73-
#else
74-
ModuleBase::WARNING_QUIT("KernelXC", "to calculate xc-kernel in LR-TDDFT, compile with LIBXC");
75-
#endif
76-
}
23+
if (std::set<std::string>({ "lda", "pwlda", "pbe", "hse" }).count(xc_kernel)) { this->set_integral_func(this->spin_type_, this->xc_type_); }
7724
}
7825

79-
void PotHxcLR::cal_v_eff(double** rho, const UnitCell* ucell, ModuleBase::matrix& v_eff, const std::vector<int>& ispin_op)
26+
void PotHxcLR::cal_v_eff(double** rho, const UnitCell& ucell, ModuleBase::matrix& v_eff, const std::vector<int>& ispin_op)
8027
{
8128
ModuleBase::TITLE("PotHxcLR", "cal_v_eff");
8229
ModuleBase::timer::tick("PotHxcLR", "cal_v_eff");
@@ -86,10 +33,10 @@ namespace LR
8633
switch (this->spin_type_)
8734
{
8835
case SpinType::S1: case SpinType::S2_updown:
89-
v_eff += elecstate::H_Hartree_pw::v_hartree(*ucell, const_cast<ModulePW::PW_Basis*>(this->rho_basis_), 1, rho);
36+
v_eff += elecstate::H_Hartree_pw::v_hartree(ucell, const_cast<ModulePW::PW_Basis*>(&this->rho_basis_), 1, rho);
9037
break;
9138
case SpinType::S2_singlet:
92-
v_eff += 2 * elecstate::H_Hartree_pw::v_hartree(*ucell, const_cast<ModulePW::PW_Basis*>(this->rho_basis_), 1, rho);
39+
v_eff += 2 * elecstate::H_Hartree_pw::v_hartree(ucell, const_cast<ModulePW::PW_Basis*>(&this->rho_basis_), 1, rho);
9340
break;
9441
default:
9542
break;
@@ -171,7 +118,7 @@ namespace LR
171118
// std::cout << std::endl;};
172119

173120
std::vector<ModuleBase::Vector3<double>> drho(nrxx); // transition density gradient
174-
LR_Util::grad(rho, drho.data(), *(this->rho_basis_), this->tpiba_);
121+
LR_Util::grad(rho, drho.data(), this->rho_basis_, this->tpiba_);
175122

176123
std::vector<double> vxc_tmp(nrxx, 0.0);
177124

@@ -183,7 +130,7 @@ namespace LR
183130
+ fxc.v2sigma2_4drho.at(ir) * (fxc.drho_gs.at(ir) * drho.at(ir))
184131
+ drho.at(ir) * fxc.vsigma.at(ir) * 2.);
185132
}
186-
XC_Functional::grad_dot(e_drho.data(), vxc_tmp.data(), this->rho_basis_, this->tpiba_);
133+
XC_Functional::grad_dot(e_drho.data(), vxc_tmp.data(), &this->rho_basis_, this->tpiba_);
187134

188135
// 2. $f^{\rho\rho}\rho_1+2f^{\rho\sigma}\nabla\rho\cdot\nabla\rho_1$
189136
for (int ir = 0;ir < nrxx;++ir)

‎source/module_lr/potentials/pot_hxc_lrtd.h‎

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,38 +1,38 @@
11
#pragma once
2-
#include "module_elecstate/potentials/pot_base.h"
32
#include "module_elecstate/potentials/H_Hartree_pw.h"
43
#include "xc_kernel.h"
54
#include <unordered_map>
65
#include <memory>
76

8-
class Parallel_Grid;
97
namespace LR
108
{
11-
class PotHxcLR : public elecstate::PotBase
9+
class PotHxcLR
1210
{
1311
public:
1412
/// S1: K^Hartree + K^xc
1513
/// S2_singlet: 2*K^Hartree + K^xc_{upup} + K^xc_{updown}
1614
/// S2_triplet: K^xc_{upup} - K^xc_{updown}
1715
/// 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)
1816
enum SpinType { S1 = 0, S2_singlet = 1, S2_triplet = 2, S2_updown = 3 };
17+
/// XCType here is to determin the method of integration from kernel to potential, not the way calculating the kernel
1918
enum XCType { None = 0, LDA = 1, GGA = 2, HYB_GGA = 4 };
2019
/// constructor for exchange-correlation kernel
21-
PotHxcLR(const std::string& xc_kernel, const ModulePW::PW_Basis* rho_basis,
22-
const UnitCell* ucell, const Charge* chg_gs/*ground state*/, const Parallel_Grid& pgrid,
20+
PotHxcLR(const std::string& xc_kernel, const ModulePW::PW_Basis& rho_basis,
21+
const UnitCell& ucell, const Charge& chg_gs/*ground state*/, const Parallel_Grid& pgrid,
2322
const SpinType& st = SpinType::S1, const std::vector<std::string>& lr_init_xc_kernel = { "from_chg_groundstate" });
2423
~PotHxcLR() {}
25-
void cal_v_eff(const Charge* chg/*excited state*/, const UnitCell* ucell, ModuleBase::matrix& v_eff) override {};
26-
void cal_v_eff(double** rho, const UnitCell* ucell, ModuleBase::matrix& v_eff, const std::vector<int>& ispin_op = { 0,0 });
27-
int nrxx;
24+
void cal_v_eff(double** rho, const UnitCell& ucell, ModuleBase::matrix& v_eff, const std::vector<int>& ispin_op = { 0,0 });
25+
const int& nrxx = nrxx_;
2826
private:
29-
int nspin;
30-
std::unique_ptr<elecstate::PotHartree> pot_hartree;
27+
const ModulePW::PW_Basis& rho_basis_;
28+
const int nspin_ = 1;
29+
const int nrxx_ = 1;
30+
std::unique_ptr<elecstate::PotHartree> pot_hartree_;
3131
/// different components of local and semi-local xc kernels:
3232
/// LDA: v2rho2
3333
/// GGA: v2rho2, v2rhosigma, v2sigma2
3434
/// meta-GGA: v2rho2, v2rhosigma, v2sigma2, v2rholap, v2rhotau, v2sigmalap, v2sigmatau, v2laptau, v2lap2, v2tau2
35-
KernelXC xc_kernel_components_;
35+
const KernelXC xc_kernel_components_;
3636
const std::string xc_kernel_;
3737
const double& tpiba_;
3838
const SpinType spin_type_ = SpinType::S1;

0 commit comments

Comments
 (0)