Skip to content

Commit 84ce690

Browse files
maki49maki49
committed
fix op-exx in Z-vector eq.
1 parent 5735afb commit 84ce690

File tree

5 files changed

+40
-25
lines changed

5 files changed

+40
-25
lines changed

‎source/module_lr/Grad/esolver_lr_grad.cpp‎

Lines changed: 24 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -266,39 +266,44 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
266266
}
267267

268268
ModuleBase::matrix force_hxc_dmtrans = lr_force.cal_force_hxc_dmtrans(dm_trans_real, *this->pot[ispin]);
269-
std::cout << "Force (Hxc-DMTrans term) of state " << istate << ": " << std::endl;
270-
LR_Util::print_value(force_hxc_dmtrans.c, ucell.nat, 3);
269+
if (PARAM.inp.test_force)
270+
ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "HXC DMTRANS FORCE (eV/Angstrom)", force_hxc_dmtrans, false);
271271

272272
const elecstate::DensityMatrix<T, double>& dm_gs = this->cal_dm_gs();
273273
ModuleBase::matrix force_hamiltgs_relaxed_diff = lr_force.cal_force_hamilt_gs_dm_relaxed_diff(relaxed_diff_dm_real, dm_gs, /*with_ewald=*/false);
274-
std::cout << "Force (GS-(T+Z) term) of state " << istate << ": " << std::endl;
275-
LR_Util::print_value(force_hamiltgs_relaxed_diff.c, ucell.nat, 3);
274+
if (PARAM.inp.test_force)
275+
ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "H_GS-(T+Z) FORCE (without EXX) (eV/Angstrom)", force_hamiltgs_relaxed_diff, false);
276276

277277
ModuleBase::matrix force_overlap_edm = lr_force.cal_force_overlap_edm(edm_real); // "-" sign has been included in the force factor
278-
std::cout << "Force (Overlap-EDM term) of state " << istate << ": " << std::endl;
279-
LR_Util::print_value(force_overlap_edm.c, ucell.nat, 3);
278+
if (PARAM.inp.test_force)
279+
ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "OVERLAP-EDM FORCE (eV/Angstrom)", force_overlap_edm, false);
280280

281281
#ifdef __EXX
282282
const double& alpha = this->exx_info.info_global.hybrid_alpha;
283283

284-
const auto& Ds_trans = get_exx_Ds_spin1(dm_trans, this->ucell, this->kv, this->paraMat_);
285-
ModuleBase::matrix force_exx_dmtrans = lr_force.cal_force_exx_dm_trans(Ds_trans, alpha);
286-
std::cout << "Force (EXX-DMTrans term) of state " << istate << ": " << std::endl;
287-
LR_Util::print_value(force_exx_dmtrans.c, ucell.nat, 3);
288-
289-
const auto& Ds_gs = get_exx_Ds_gs(dm_gs, this->ucell, this->kv, this->paraMat_, this->nspin);
290-
const auto& Ds_relaxed_diff = get_exx_Ds_spin1(relaxed_diff_dm, this->ucell, this->kv, this->paraMat_);
291-
ModuleBase::matrix force_exx_gs_diff = lr_force.cal_force_exx_gs_dm_relaxed_diff(Ds_gs, Ds_relaxed_diff, alpha);
292-
std::cout << "Force (EXX-GS-(T+Z) term) of state " << istate << ": " << std::endl;
293-
LR_Util::print_value(force_exx_gs_diff.c, ucell.nat, 3);
284+
if (LR::exx_kernel_list().count(xc_kernel))
285+
{
286+
const auto& Ds_trans = get_exx_Ds_spin1(dm_trans, this->ucell, this->kv, this->paraMat_);
287+
ModuleBase::matrix force_exx_dmtrans = lr_force.cal_force_exx_dm_trans(Ds_trans, alpha);
288+
if (PARAM.inp.test_force)
289+
ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "EXX DMTRANS FORCE (eV/Angstrom)", force_exx_dmtrans, false);
290+
force_hxc_dmtrans += force_exx_dmtrans;
294291

295-
// add exx force to the corresponding terms
296-
force_hxc_dmtrans += force_exx_dmtrans;
297-
force_hamiltgs_relaxed_diff += force_exx_gs_diff;
292+
}
293+
if (LR::exx_kernel_list().count(PARAM.inp.dft_functional))
294+
{
295+
const auto& Ds_gs = get_exx_Ds_gs(dm_gs, this->ucell, this->kv, this->paraMat_, this->nspin);
296+
const auto& Ds_relaxed_diff = get_exx_Ds_spin1(relaxed_diff_dm, this->ucell, this->kv, this->paraMat_);
297+
ModuleBase::matrix force_exx_gs_diff = lr_force.cal_force_exx_gs_dm_relaxed_diff(Ds_gs, Ds_relaxed_diff, alpha);
298+
if (PARAM.inp.test_force)
299+
ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "EXX GS-(T+Z) FORCE (eV/Angstrom)", force_exx_gs_diff, false);
300+
force_hamiltgs_relaxed_diff += force_exx_gs_diff;
301+
}
298302
#endif
299303
forces[istate] = force_hxc_dmtrans + force_hamiltgs_relaxed_diff + force_overlap_edm;
300304
}
301305
ModuleBase::timer::tick("ESolver_LR", "cal_force");
306+
// total force
302307
print_force(forces, std::cout);
303308
print_force(forces, GlobalV::ofs_running);
304309
return forces;

‎source/module_lr/Grad/multipliers/hamilt_zeq_left.h‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -49,15 +49,16 @@ namespace LR
4949
ModuleBase::TITLE("Z_vector_L", "Z_vector_L");
5050
this->DM_trans = LR_Util::make_unique<elecstate::DensityMatrix<T, T>>(&pmat, 1, kv.kvec_d, this->nk);
5151
LR_Util::initialize_DMR(*this->DM_trans, pmat, ucell, gd, orb_cutoff);
52-
// 1. $2\sum_bX_{ib}K_{ab}[D^X]-2\sum_jX_{ja}K_{ij}[D^X]$
52+
// $H_{ia}[D^Z]$, equals to $2K_{ab}[D^Z]$ when $D^Z$ is symmetrized
53+
// 1. diag term in A
5354
this->ops = new OperatorLRDiag<T>(eig_ks.c, pX[0], kv.get_nks() / nspin, nocc[0], nvirt[0]);
54-
// 2. $H_{ia}[D^Z]$, equals to $2K_{ab}[D^Z]$ when $D^Z$ is symmetrized
55+
// 2. Hessian (A+B) with GS XC kernel
5556
hamilt::Operator<T>* op_hz = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
5657
*this->DM_trans, gint, pot_hxc_gs, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
5758
{ 0 }, 2.0, ATYPE::CC_vo);
5859
this->ops->add(op_hz);
5960
#ifdef __EXX
60-
if (exx_kernel_list().count(xc_kernel))
61+
if (exx_kernel_list().count(PARAM.inp.dft_functional))
6162
{
6263
hamilt::Operator<T>* op_hz_exx = new OperatorLREXX<T>(nspin, naos, nocc[0], nvirt[0], ucell, psi_ks,
6364
*this->DM_trans, exx_lri, kv, pX[0], pc, pmat,

‎source/module_lr/Grad/multipliers/hamilt_zeq_right.h‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,7 @@ namespace LR
5656

5757
// note: calculation_type cannot repeated, or it will be ignored in ops->add()
5858
// 1. $2\sum_bX_{ib}K_{ab}[D^X]-2\sum_jX_{ja}K_{ij}[D^X]$
59+
// kernel: excited state
5960
this->ops = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
6061
*this->DM_trans, gint, pot, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
6162
{ 0 }, -2.0, ATYPE::CXC);
@@ -70,12 +71,13 @@ namespace LR
7071
}
7172
#endif
7273
// 2. $H_{ia}[T]$, equals to $2K_{ab}[T]$ when $T$ is symmetrized
74+
// kernel: ground state
7375
hamilt::Operator<T>* op_ht = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
7476
*this->DM_diff, gint, pot_hxc_gs, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
7577
{ 0 }, T(-2.0), ATYPE::CC_vo, hamilt::calculation_type::lr_dmdiff_hxc);
7678
this->ops->add(op_ht);
7779
#ifdef __EXX
78-
if (exx_kernel_list().count(xc_kernel))
80+
if (exx_kernel_list().count(PARAM.inp.dft_functional))
7981
{
8082
hamilt::Operator<T>* op_ht_exx = new OperatorLREXX<T>(nspin, naos, nocc[0], nvirt[0], ucell, psi_ks,
8183
*this->DM_diff, exx_lri, kv, pX[0], pc, pmat,

‎source/module_lr/esolver_lrtd_lcao.cpp‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -437,8 +437,12 @@ LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu
437437
this->gint_->initialize_pvpR(ucell, &this->gd, 1); // always use nspin=1 for transition density
438438

439439
// if EXX from scratch, init 2-center integral and calculate Cs, Vs
440+
// when:
441+
// 1. EXX xc_kernel
442+
// 2. cal_force with ground state with EXX functional
440443
#ifdef __EXX
441-
if ((xc_kernel == "hf" || xc_kernel == "hse") && this->input.lr_solver != "spectrum")
444+
if (((xc_kernel == "hf" || xc_kernel == "hse") && this->input.lr_solver != "spectrum")
445+
|| (PARAM.inp.cal_force && (PARAM.inp.dft_functional == "hf" || PARAM.inp.dft_functional == "hse")))
442446
{
443447
// set ccp_type according to the xc_kernel
444448
if (xc_kernel == "hf") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf; }

‎source/module_lr/operator_casida/operator_lr_exx.h‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -68,7 +68,10 @@ namespace LR
6868
this->BvK_cells = RI_Util::get_Born_von_Karmen_cells(period);
6969

7070
this->allocate_Ds_onebase();
71-
this->exx_lri.lock()->Hexxs.resize(1);
71+
if (!this->exx_lri.expired())
72+
{
73+
this->exx_lri.lock()->Hexxs.resize(1);
74+
}
7275
};
7376

7477
void init(const int ik_in) override {};

0 commit comments

Comments
 (0)