Skip to content

Commit e8c24fb

Browse files
maki49maki49
committed
KEY fix: Hxc Hellmann-Feynman term with (T+DZ)
1 parent 61d8ddc commit e8c24fb

File tree

5 files changed

+98
-28
lines changed

5 files changed

+98
-28
lines changed

‎source/module_lr/Grad/esolver_lr_grad.cpp‎

Lines changed: 20 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -75,7 +75,7 @@ inline auto get_exx_Ds_spin1(const elecstate::DensityMatrix<T, T>& dm,
7575
return RI_2D_Comm::split_m2D_ktoR<T>(ucell, kv, DMk_trans_pointer, pmat, /*nspin=*/1)[0];
7676
}
7777
template <typename T>
78-
inline auto get_exx_Ds_gs(elecstate::DensityMatrix<T, double>& dm,
78+
inline auto get_exx_Ds_gs(const elecstate::DensityMatrix<T, double>& dm,
7979
const UnitCell& ucell, const K_Vectors& kv, const Parallel_Orbitals& pmat, const int nspin)
8080
-> std::map<int, std::map<std::pair<int, std::array<int, 3>>, RI::Tensor<T>>>
8181
{
@@ -167,7 +167,8 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
167167
const auto& c = LR_Util::get_psi_spin(*this->psi_ks, ispin, this->nk); // wavefunction coefficients of ground state
168168

169169
// calculate the force (the partial gradient of Lagrangian)
170-
LR_Force<T> lr_force(this->ucell, this->kv.kvec_d, this->paraMat_, *this->pw_rho, this->locpp, this->sf, this->gd, this->gint_, this->two_center_bundle_
170+
LR_Force<T> lr_force(this->ucell, this->kv.kvec_d, this->paraMat_,
171+
*this->pw_rhod, *this->pw_rho, this->locpp, this->sf, this->gd, this->gint_, this->two_center_bundle_
171172
#ifdef __EXX
172173
, std::weak_ptr<Exx_LRI<T>>(this->exx_lri), this->exx_info.info_global.hybrid_alpha
173174
#endif
@@ -245,7 +246,8 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
245246
std::cout << "Force (Hxc-DMTrans term) of state " << istate << ": " << std::endl;
246247
LR_Util::print_value(force_hxc_dmtrans.c, ucell.nat, 3);
247248

248-
ModuleBase::matrix force_hamiltgs_relaxed_diff = lr_force.cal_force_hamilt_gs_dm_relaxed_diff(relaxed_diff_dm_real, *pot_gs, /*with_ewald=*/false);
249+
const elecstate::DensityMatrix<T, double>& dm_gs = this->cal_dm_gs();
250+
ModuleBase::matrix force_hamiltgs_relaxed_diff = lr_force.cal_force_hamilt_gs_dm_relaxed_diff(relaxed_diff_dm_real, dm_gs, /*with_ewald=*/false);
249251
std::cout << "Force (GS-(T+Z) term) of state " << istate << ": " << std::endl;
250252
LR_Util::print_value(force_hamiltgs_relaxed_diff.c, ucell.nat, 3);
251253

@@ -261,8 +263,6 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
261263
std::cout << "Force (EXX-DMTrans term) of state " << istate << ": " << std::endl;
262264
LR_Util::print_value(force_exx_dmtrans.c, ucell.nat, 3);
263265

264-
elecstate::DensityMatrix<T, double> dm_gs(&this->paraMat_, this->nspin, this->kv.kvec_d, this->nk); //DX
265-
elecstate::cal_dm_psi(&this->paraMat_all_, this->wg_ks_all, *this->psi_ks_all, dm_gs); // nbands is important here
266266
const auto& Ds_gs = get_exx_Ds_gs(dm_gs, this->ucell, this->kv, this->paraMat_, this->nspin);
267267
const auto& Ds_relaxed_diff = get_exx_Ds_spin1(relaxed_diff_dm, this->ucell, this->kv, this->paraMat_);
268268
ModuleBase::matrix force_exx_gs_diff = lr_force.cal_force_exx_gs_dm_relaxed_diff(Ds_gs, Ds_relaxed_diff, alpha);
@@ -281,20 +281,28 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
281281
return forces;
282282
}
283283

284+
template<typename T, typename TR>
285+
elecstate::DensityMatrix<T, double> LR::ESolver_LR<T, TR>::cal_dm_gs()
286+
{
287+
elecstate::DensityMatrix<T, double> dm_gs(&this->paraMat_, this->nspin, this->kv.kvec_d, this->nk);
288+
elecstate::cal_dm_psi(&this->paraMat_all_, this->wg_ks_all, *this->psi_ks_all, dm_gs); // nbands is important here
289+
LR_Util::initialize_DMR(dm_gs, this->paraMat_, this->ucell, this->gd, this->orb_cutoff_); // nbands is not important here
290+
dm_gs.cal_DMR();
291+
return dm_gs;
292+
}
293+
284294
template<typename T, typename TR>
285295
void LR::ESolver_LR<T, TR>::test_force()
286296
{
287-
LR_Force<T> lr_force(this->ucell, this->kv.kvec_d, this->paraMat_, *this->pw_rho, this->locpp, this->sf, this->gd, this->gint_, this->two_center_bundle_
297+
LR_Force<T> lr_force(this->ucell, this->kv.kvec_d, this->paraMat_, *this->pw_rhod, *this->pw_rho,
298+
this->locpp, this->sf, this->gd, this->gint_, this->two_center_bundle_
288299
#ifdef __EXX
289300
, std::weak_ptr<Exx_LRI<T>>(this->exx_lri), this->exx_info.info_global.hybrid_alpha
290301
#endif
291302
);
292303

293-
elecstate::DensityMatrix<T, double> dm_gs(&this->paraMat_, this->nspin, this->kv.kvec_d, this->nk); //DX
294-
elecstate::cal_dm_psi(&this->paraMat_all_, this->wg_ks_all, *this->psi_ks_all, dm_gs); // nbands is important here
295-
LR_Util::initialize_DMR(dm_gs, this->paraMat_, this->ucell, this->gd, this->orb_cutoff_); // nbands is not important here
296-
dm_gs.cal_DMR();
297-
LR_Util::print_DMR(dm_gs, this->ucell.nat, "DM(R) of ground state");
304+
const elecstate::DensityMatrix<T, double>& dm_gs = this->cal_dm_gs();
305+
// LR_Util::print_DMR(dm_gs, this->ucell.nat, "DM(R) of ground state");
298306
///========================== test 1: reproduce the force of ground state =========================
299307
// energy density matrix of the ground state
300308
elecstate::DensityMatrix<T, double> edm_gs(&this->paraMat_, this->nspin, this->kv.kvec_d, this->nk); //DX
@@ -305,7 +313,7 @@ void LR::ESolver_LR<T, TR>::test_force()
305313
LR_Util::initialize_DMR(edm_gs, this->paraMat_, this->ucell, this->gd, this->orb_cutoff_);
306314
edm_gs.cal_DMR();
307315
// ground-state force
308-
ModuleBase::matrix force_gs = lr_force.reproduce_force_gs(dm_gs, edm_gs, *this->pot_gs);
316+
ModuleBase::matrix force_gs = lr_force.reproduce_force_gs(dm_gs, edm_gs);
309317
ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "Ground State FORCE (eV/Angstrom)", force_gs, false);
310318
/// ======================================= END test 1 =========================================
311319
///========================== test 2: reproduce the DX Hartree term =========================

‎source/module_lr/Grad/force/lr_force.cpp‎

Lines changed: 65 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -25,9 +25,35 @@ namespace LR
2525
return chr;
2626
}
2727

28+
template<typename TK>
29+
elecstate::Potential LR_Force<TK>::dm_to_hxc_potential(const elecstate::DensityMatrix<TK, double>& dm)
30+
{
31+
double etxc = 0.0, vtxc = 0.0;
32+
elecstate::Potential pot(&this->rhodpw_, &this->rhopw_, &this->ucell_,
33+
&this->locpp_.vloc, const_cast<Structure_Factor*>(&this->sf_),
34+
nullptr/*surchem*/, &etxc, &vtxc);
35+
PARAM.inp.vh_in_h ? pot.pot_register({ "hartree", "xc" }) : pot.pot_register({ "xc" });
36+
const Charge& charge = this->dm_to_charge(dm);
37+
pot.init_pot(0, &charge); // call update_from_charge inside
38+
return pot;
39+
}
40+
41+
template<typename TK>
42+
elecstate::Potential LR_Force<TK>::local_potential()
43+
{
44+
elecstate::Potential pot(&this->rhodpw_, &this->rhopw_, &this->ucell_,
45+
&this->locpp_.vloc, const_cast<Structure_Factor*>(&this->sf_),
46+
nullptr/*surchem*/, nullptr/*etxc*/, nullptr/*vtxc*/);
47+
pot.pot_register({ "local" });
48+
pot.init_pot(0, nullptr);
49+
return pot;
50+
}
51+
2852
template<typename TK>
2953
ModuleBase::matrix LR_Force<TK>::cal_force_hamilt_gs_dm_relaxed_diff(const elecstate::DensityMatrix<TK, double>& relax_diff_dm,
30-
const elecstate::Potential& pot_gs, const bool with_ewald)
54+
const elecstate::DensityMatrix<TK, double>& dm_gs,
55+
// const elecstate::Potential& pot_gs,
56+
const bool with_ewald)
3157
{
3258
const Charge chr_diff_relaxed = dm_to_charge(relax_diff_dm);
3359

@@ -39,11 +65,40 @@ namespace LR
3965
// 2. nonlocal pp (Hellmann-Feynman + Pulay)
4066
ModuleBase::matrix fvnl = cal_force_nonlocal(this->ucell_, this->kvec_d_, this->gd_, this->two_center_bundle_, relax_diff_dm);
4167

42-
// 3. local pp (Pulay) + Hartree + xc (grid integration)
68+
// // 3. local pp (Pulay) + Hartree + xc (grid integration)
69+
// ModuleBase::matrix fvl_dphi(this->ucell_.nat, 3);
70+
// ModuleBase::matrix stress_tmp; // no use now, only for passing into interfaces
71+
// PulayForceStress::cal_pulay_fs(relax_diff_dm.get_DMR_vector().size()/*nspin*/, fvl_dphi, stress_tmp,
72+
// relax_diff_dm, this->ucell_, &pot_gs, *this->gint_, true, false);
73+
74+
// 3.1. local pp (Pulay)
4375
ModuleBase::matrix fvl_dphi(this->ucell_.nat, 3);
4476
ModuleBase::matrix stress_tmp; // no use now, only for passing into interfaces
77+
elecstate::Potential pot_loc = this->local_potential();
4578
PulayForceStress::cal_pulay_fs(relax_diff_dm.get_DMR_vector().size()/*nspin*/, fvl_dphi, stress_tmp,
46-
relax_diff_dm, this->ucell_, &pot_gs, *this->gint_, true, false);
79+
relax_diff_dm, this->ucell_, &pot_loc, *this->gint_, true, false);
80+
81+
// 3.2. Hartree + xc (Pulay)
82+
// method 1
83+
// ModuleBase::matrix fgs_dphi(this->ucell_.nat, 3);
84+
// PulayForceStress::cal_pulay_fs(relax_diff_dm.get_DMR_vector().size()/*nspin*/, fgs_dphi, stress_tmp,
85+
// relax_diff_dm, this->ucell_, &pot_gs, *this->gint_, true, false);
86+
// ModuleBase::matrix fhxc_dphi = (fgs_dphi - fvl_dphi) * 0.5; // avoid double count of hxc Pulay term
87+
// method 2
88+
ModuleBase::matrix fhxc_dphi(this->ucell_.nat, 3);
89+
this->gint_->reset_DMRGint(dm_gs.get_DMR_vector().size());
90+
elecstate::Potential pot_hxc = this->dm_to_hxc_potential(dm_gs);
91+
this->gint_->reset_DMRGint(relax_diff_dm.get_DMR_vector().size());
92+
PulayForceStress::cal_pulay_fs(relax_diff_dm.get_DMR_vector().size()/*nspin*/, fhxc_dphi, stress_tmp,
93+
relax_diff_dm, this->ucell_, &pot_hxc, *this->gint_, true, false);
94+
fhxc_dphi *= 0.5; // avoid double count
95+
96+
// 3.3 Hartree + xc (Hellmann-Feynman)
97+
ModuleBase::matrix fhxc_dvhxc(this->ucell_.nat, 3);
98+
elecstate::Potential pot_hxc_relaxed_diff = this->dm_to_hxc_potential(relax_diff_dm);
99+
PulayForceStress::cal_pulay_fs(1/*nspin*/, fhxc_dvhxc, stress_tmp,
100+
relax_diff_dm, this->ucell_, &pot_hxc_relaxed_diff, *this->gint_, true, false);
101+
// fhxc_dvhxc *= 0.5; // avoid double count, but nspin=2 of ground-state dm cancels it here
47102

48103
// 4. kinetic (Pulay)
49104
std::vector<hamilt::HContainer<double>> dT = cal_hs_grad('T', this->ucell_, this->pv_, this->gd_, this->two_center_bundle_);
@@ -54,11 +109,13 @@ namespace LR
54109
ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "PW FORCE (eV/Angstrom)", f_pw, false);
55110
ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "NONLOCAL FORCE (eV/Angstrom)", fvnl, false);
56111
ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "KINETIC FORCE (eV/Angstrom)", ft_dphi, false);
57-
ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "LOCAL Pulay FORCE (pp+Hartree+XC) (eV/Angstrom)", fvl_dphi, false);
112+
ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "LOCAL-PP Pulay FORCE (eV/Angstrom)", fvl_dphi, false);
113+
ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "HARTREE+XC Pulay FORCE (eV/Angstrom)", fhxc_dphi, false);
114+
ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "HARTREE+XC Hellmann-Feynman FORCE (eV/Angstrom)", fhxc_dvhxc, false);
58115
}
59116

60117
// from the formula, we do not need the non-ortho term (overlap*edm) here.
61-
return f_pw + fvnl + ft_dphi + fvl_dphi;
118+
return f_pw + fvnl + ft_dphi + fvl_dphi + fhxc_dphi + fhxc_dvhxc;
62119
}
63120
template<typename TK>
64121
ModuleBase::matrix LR_Force<TK>::cal_force_hxc_dmtrans(const elecstate::DensityMatrix<TK, double>& dm_trans, const PotHxcLR& pot_hxc)
@@ -85,12 +142,11 @@ namespace LR
85142
template<typename TK>
86143
ModuleBase::matrix LR_Force<TK>::reproduce_force_gs(
87144
const elecstate::DensityMatrix<TK, double>& dm_gs,
88-
const elecstate::DensityMatrix<TK, double>& edm_gs,
89-
const elecstate::Potential& pot_gs)
145+
const elecstate::DensityMatrix<TK, double>& edm_gs)
90146
{
91147
this->gint_->reset_DMRGint(PARAM.inp.nspin);
92-
// Hartree+xc term
93-
ModuleBase::matrix f_gs_hf_pulay = cal_force_hamilt_gs_dm_relaxed_diff(dm_gs, pot_gs); // pw+vnl+t_dphi+vl_dphi
148+
// local + Hartree + xc term, including Hellmann-Feynman and Pulay
149+
ModuleBase::matrix f_gs_hf_pulay = cal_force_hamilt_gs_dm_relaxed_diff(dm_gs, dm_gs); // pw+vnl+t_dphi+vl_dphi
94150
// edm term
95151
ModuleBase::matrix f_nonortho = cal_force_overlap_edm(edm_gs); // overlap
96152
this->gint_->reset_DMRGint(1);

‎source/module_lr/Grad/force/lr_force.h‎

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@ namespace LR
1717
LR_Force(const UnitCell& ucell,
1818
const std::vector<ModuleBase::Vector3<double>>& kvec_d,
1919
const Parallel_Orbitals& pv,
20+
const ModulePW::PW_Basis& rhodpw,
2021
const ModulePW::PW_Basis& rhopw,
2122
const pseudopot_cell_vl& locpp,
2223
const Structure_Factor& sf,
@@ -29,7 +30,7 @@ namespace LR
2930
#endif
3031
)///< for 2-center integrals
3132
: ucell_(ucell), kvec_d_(kvec_d), pv_(pv),
32-
rhopw_(rhopw), sf_(sf), locpp_(locpp), gd_(gd),
33+
rhodpw_(rhodpw), rhopw_(rhopw), sf_(sf), locpp_(locpp), gd_(gd),
3334
gint_(gint), two_center_bundle_(two_center_bundle)
3435
#ifdef __EXX
3536
, exx_lri_(exx_lri_in), alpha_(alpha)
@@ -39,7 +40,8 @@ namespace LR
3940

4041
/// 1. $Tr[H_{GS}^x * (T+D^Z)]$, where GS=groud state and $(T+D^Z)$ is the relaxed difference density matrix
4142
ModuleBase::matrix cal_force_hamilt_gs_dm_relaxed_diff(const elecstate::DensityMatrix<TK, double>& relaxed_diff_dm,
42-
const elecstate::Potential& pot_gs, const bool with_ewald = true);
43+
const elecstate::DensityMatrix<TK, double>& dm_gs, const bool with_ewald = true);
44+
// const elecstate::Potential& pot_gs, const bool with_ewald = true);
4345

4446
/// 2. $Tr[S^x * (EDM)]
4547
ModuleBase::matrix cal_force_overlap_edm(const elecstate::DensityMatrix<TK, double>& edm);
@@ -62,8 +64,8 @@ namespace LR
6264
// test functions
6365
/// reproduce the force of the ground state
6466
ModuleBase::matrix reproduce_force_gs(const elecstate::DensityMatrix<TK, double>& dm_gs,
65-
const elecstate::DensityMatrix<TK, double>& edm_gs,
66-
const elecstate::Potential& pot_gs);
67+
const elecstate::DensityMatrix<TK, double>& edm_gs);
68+
6769
/// repreduce the ground state local term
6870
ModuleBase::matrix reproduce_force_gs_loc(const elecstate::DensityMatrix<TK, double>& dm_gs,
6971
const elecstate::Potential& pot_gs);
@@ -72,6 +74,7 @@ namespace LR
7274
const UnitCell& ucell_;
7375
const std::vector<ModuleBase::Vector3<double>>& kvec_d_;
7476
const Parallel_Orbitals& pv_;
77+
const ModulePW::PW_Basis& rhodpw_;
7578
const ModulePW::PW_Basis& rhopw_;
7679
const pseudopot_cell_vl& locpp_;
7780
const Structure_Factor& sf_;
@@ -84,7 +87,8 @@ namespace LR
8487
#endif
8588

8689
Charge dm_to_charge(const elecstate::DensityMatrix<TK, double>& dm);
87-
90+
elecstate::Potential dm_to_hxc_potential(const elecstate::DensityMatrix<TK, double>& dm);
91+
elecstate::Potential local_potential();
8892
// probably move frome the ground state?
8993
// void build_dHS()
9094
};

‎source/module_lr/esolver_lrtd_lcao.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -148,6 +148,7 @@ namespace LR
148148
ct::Tensor solve_zvector_eqation(const int ispin);
149149
std::vector<ModuleBase::matrix> cal_force(const int ispin);
150150
void test_force(); // test: reproduce the force of ground state
151+
elecstate::DensityMatrix<T, double> cal_dm_gs(); ///< ground-state density matrix
151152

152153
#ifdef __EXX
153154
/// Tdata of Exx_LRI is same as T, for the reason, see operator_lr_exx.h

‎source/module_lr/operator_casida/operator_lr_exx.cpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,7 @@ namespace LR
4141
{
4242
// term1: Co -> CvX, i.e. [CvX] Cv
4343
DMBand<double> dm_band1(ucell, pmat, this->kv.kvec_c, this->BvK_cells, this->cvx_full, this->psi_ks_full);
44-
dm_band1.eval(iv, io, ik);
44+
dm_band1.eval(io, iv, ik);
4545
// term2: Cv -> CoX^T, i.e. Co [CoX^T]
4646
DMBand<double> dm_band2(ucell, pmat, this->kv.kvec_c, this->BvK_cells, this->psi_ks_full, this->coxt_full);
4747
dm_band2.eval(io, iv, ik);
@@ -149,7 +149,8 @@ namespace LR
149149

150150
// 3. set [AX]_iak = DM_onbase * Hexxs for each occ-virt pair and each k-point
151151
// caution: parrallel
152-
if (PARAM.inp.cal_force) this->cal_coxt_cvx(psi_in);
152+
if (this->dm_pq_ == MO_TO_AO_TYPE::CXC || this->dm_pq_ == MO_TO_AO_TYPE::CXC_o)
153+
this->cal_coxt_cvx(psi_in);
153154
for (int io = 0;io < this->nocc;++io)
154155
{
155156
for (int iv = 0;iv < this->nvirt;++iv)

0 commit comments

Comments
 (0)