Skip to content

Commit 1a56745

Browse files
maki49maki49
committed
fix dm-to-charge, ground state Hxc potential and remove the redundant Ewald term
1 parent 0d7373f commit 1a56745

File tree

7 files changed

+38
-23
lines changed

7 files changed

+38
-23
lines changed

‎source/module_lr/Grad/esolver_lr_grad.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -76,8 +76,8 @@ ct::Tensor LR::ESolver_LR<T, TR>::solve_zvector_eqation(const int ispin)
7676
#ifdef __EXX
7777
std::weak_ptr<Exx_LRI<T>>(this->exx_lri), this->exx_info.info_global.hybrid_alpha,
7878
#endif
79-
this->gint_, std::weak_ptr<PotHxcLR>(this->pot[ispin]), this->kv,
80-
this->paraX_, this->paraC_, this->paraMat_, this->spin_types[ispin]);
79+
this->gint_, std::weak_ptr<PotHxcLR>(this->pot[ispin]), std::weak_ptr<PotHxcLR>(this->pot_hxc_gs),
80+
this->kv, this->paraX_, this->paraC_, this->paraMat_, this->spin_types[ispin]);
8181
ModuleBase::timer::tick("ESolver_LR", "solve_zvector_eqation");
8282
return Z;
8383
}

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

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@ class ForcePWTerms
1616
const ModulePW::PW_Basis& rhopw,
1717
const pseudopot_cell_vl& locpp,
1818
const Structure_Factor& sf,
19+
const bool with_ewald = true,
1920
const elecstate::ElecState* pelec = nullptr)
2021
{
2122
if (PARAM.inp.nspin == 4) { throw std::runtime_error("ForcePWTerms: nspin=4 is not supported."); }
@@ -30,8 +31,10 @@ class ForcePWTerms
3031
//--------------------------------------------------------
3132
// ewald force: use plane wave only.
3233
//--------------------------------------------------------
33-
f_pw.cal_force_ew(ucell, fewalds, &rhopw, &sf);
34-
34+
if (with_ewald)
35+
{
36+
f_pw.cal_force_ew(ucell, fewalds, &rhopw, &sf);
37+
}
3538
//--------------------------------------------------------
3639
// force due to core correlation.
3740
//--------------------------------------------------------

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

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include "pulay_force_hcontainer.h"
44
#include "module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h" // only for gint terms
55
#include "module_lr/utils/lr_util.h"
6+
// #include "module_lr/utils/lr_util_hcontainer.h"
67
namespace LR
78
{
89
template<typename TK>
@@ -20,11 +21,7 @@ namespace LR
2021
this->gint_->transfer_DM2DtoGrid(dm.get_DMR_vector()); // 2d block to grid
2122
Gint_inout inout_rho(chr.rho, Gint_Tools::job_type::rho, nspin_dm, false);
2223
this->gint_->cal_gint(&inout_rho);
23-
// point spin2 charge to spin1
24-
if (nspin_dm == 1 && nspin_global == 2)
25-
{
26-
chr.rho[1] = chr.rho[0];
27-
}
24+
// if (nspin_dm == 1 && nspin_global == 2), chr.rho[1][irxx]=0 has been set in Charge::allocate()
2825
return chr;
2926
}
3027

@@ -34,7 +31,7 @@ namespace LR
3431
const Charge chr_diff_relaxed = dm_to_charge(relax_diff_dm);
3532

3633
// 1. local pp (Hellmann-Feynman)(fvl_dvl) + ewald + core correction (+ self-consistent charge)
37-
ModuleBase::matrix f_pw = ForcePWTerms<double>()(this->ucell_, chr_diff_relaxed, this->rhopw_, this->locpp_, this->sf_);
34+
ModuleBase::matrix f_pw = ForcePWTerms<double>()(this->ucell_, chr_diff_relaxed, this->rhopw_, this->locpp_, this->sf_, /*with_ewald=*/false);
3835

3936
// 2. nonlocal pp (Hellmann-Feynman + Pulay)
4037
ModuleBase::matrix fvnl = cal_force_nonlocal(this->ucell_, this->kvec_d_, this->gd_, this->two_center_bundle_, relax_diff_dm);
@@ -71,7 +68,9 @@ namespace LR
7168
{
7269
// const double* dS[3] = { dSloc_x, dSloc_y, dSloc_z };
7370
std::vector<hamilt::HContainer<double>> dS = cal_hs_grad('S', this->ucell_, this->pv_, this->gd_, this->two_center_bundle_);
74-
// cal_pulay_fs should support HContainer
71+
// test: output dS
72+
// std::cout << "dS in 3 directions:\n";
73+
// for (int i = 0;i < 3;++i) { LR_Util::print_HR(dS.at(i), this->ucell_.nat, "dS" + std::to_string(i)); }
7574
ModuleBase::matrix foverlap = PulayForceStress::cal_pulay_fs(edm, this->ucell_, dS, -1.);
7675
if (PARAM.inp.test_force)
7776
{

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

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -72,7 +72,7 @@ namespace LR
7272
const double& exx_alpha,
7373
#endif
7474
TGint* gint,
75-
std::weak_ptr<PotHxcLR> pot,
75+
std::weak_ptr<PotHxcLR> pot_hxc_gs,
7676
const K_Vectors& kv,
7777
const std::vector<Parallel_2D>& px,
7878
const Parallel_2D& pc,
@@ -90,9 +90,9 @@ namespace LR
9090
elecstate::DensityMatrix<T, T> DM_diff_relaxed(&pmat, 1, kv.kvec_d, nk); //T+DZ
9191
LR_Util::initialize_DMR(DM_diff_relaxed, pmat, ucell, gd, orb_cutoff);
9292
/// operators
93-
// 1. 0.5$H_{ia}[T]$, equals to $K_{ab}[T]$ when $T$ is symmetrized
93+
// 1. 0.5$H_{ia}[T+Z]$, equals to $K_{ab}[T+Z]$ when $(T+Z)$ is symmetrized
9494
OperatorLRHxc<T> op_ht(nspin, naos, nocc, nvirt, psi_ks,
95-
DM_diff_relaxed, gint, pot, ucell, orb_cutoff, gd, kv, p_occ_occ, pc, pmat,
95+
DM_diff_relaxed, gint, pot_hxc_gs, ucell, orb_cutoff, gd, kv, p_occ_occ, pc, pmat,
9696
{ 0 }, T(1.0), ATYPE::CC_oo);
9797
// 2. $2\sum_{jb,kc} g^{xc}_{ia, jb, kc}X_{jb}X_{kc}$
9898
// use pointer here for polymorphism
@@ -101,7 +101,7 @@ namespace LR
101101
// `weak_ptr=shared_ptr` is automatically called in the constructor of OperatorLRHxc, so we don't need to do it manually
102102
// if `pot_grad` is passed into a function rather than a class, we need to write `weak_ptr=shared_ptr` explicitly
103103
std::shared_ptr<PotGradXCLR> pot_grad =
104-
std::make_shared<PotGradXCLR>(pot.lock()->xc_kernel_components, pot.lock()->get_rho_basis(), ucell, pot.lock()->nrxx);
104+
std::make_shared<PotGradXCLR>(pot_hxc_gs.lock()->xc_kernel_components, pot_hxc_gs.lock()->get_rho_basis(), ucell, pot_hxc_gs.lock()->nrxx);
105105
OperatorLRHxc<T> op_gxc(nspin, naos, nocc, nvirt, psi_ks,
106106
DM_trans, gint, pot_grad, ucell, orb_cutoff, gd, kv, p_occ_occ, pc, pmat,
107107
{ 0 }, T(-2.0), ATYPE::CC_oo);
@@ -150,7 +150,7 @@ namespace LR
150150
cal_dm_trans(0, X); // transition density matrix DX
151151
cal_dm_diff_relaxed(0, X, Z); // relaxed difference density matrix T+DZ
152152
// the 3 terms
153-
op_ht.act(/*nband=*/1, ld_oo, /*npol=*/1, X, W);
153+
op_ht.act(/*nband=*/1, ld_oo, /*npol=*/1, X, W); //comment out this line to test H[T+Z]=0
154154
if (has_local_xc) { op_gxc.act(/*nband=*/1, ld_oo, /*npol=*/1, X, W); }
155155
std::cout << "W (H[T+Z]) + W(gxc) terms: " << std::endl;
156156
LR_Util::print_value(W, nk, p_occ_occ[0].get_col_size(), p_occ_occ[0].get_row_size());

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

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@ namespace LR
2828
const double& exx_alpha,
2929
#endif
3030
TGint* gint,
31-
std::weak_ptr<PotLRBase> pot,
31+
std::weak_ptr<PotLRBase> pot_hxc_gs,
3232
const K_Vectors& kv,
3333
const std::vector<Parallel_2D>& pX,
3434
const Parallel_2D& pc,
@@ -38,16 +38,16 @@ namespace LR
3838
#ifdef __EXX
3939
exx_lri, exx_alpha,
4040
#endif
41-
gint, pot, kv, pX, pc, pmat, spin_type)
41+
gint, pot_hxc_gs, kv, pX, pc, pmat, spin_type)
4242
{
4343
ModuleBase::TITLE("Z_vector_L", "Z_vector_L");
4444
this->DM_trans = LR_Util::make_unique<elecstate::DensityMatrix<T, T>>(&pmat, 1, kv.kvec_d, this->nk);
4545
LR_Util::initialize_DMR(*this->DM_trans, pmat, ucell, gd, orb_cutoff);
4646
// 1. $2\sum_bX_{ib}K_{ab}[D^X]-2\sum_jX_{ja}K_{ij}[D^X]$
4747
this->ops = new OperatorLRDiag<T>(eig_ks.c, pX[0], kv.get_nks() / nspin, nocc[0], nvirt[0]);
48-
// 2. $H_{ia}[T]$, equals to $2K_{ab}[T]$ when $T$ is symmetrized
48+
// 2. $H_{ia}[D^Z]$, equals to $2K_{ab}[D^Z]$ when $D^Z$ is symmetrized
4949
OperatorLRHxc<T>* op_hz = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
50-
*this->DM_trans, gint, pot, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
50+
*this->DM_trans, gint, pot_hxc_gs, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
5151
{ 0 }, 2.0, ATYPE::CC_vo);
5252
this->ops->add(op_hz);
5353
#ifdef __EXX

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

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ namespace LR
2929
#endif
3030
TGint* gint,
3131
std::weak_ptr<PotHxcLR> pot,
32+
std::weak_ptr<PotHxcLR> pot_hxc_gs,
3233
const K_Vectors& kv,
3334
const std::vector<Parallel_2D>& pX,
3435
const Parallel_2D& pc,
@@ -51,7 +52,7 @@ namespace LR
5152
{ 0 }, -2.0, ATYPE::CXC);
5253
// 2. $H_{ia}[T]$, equals to $2K_{ab}[T]$ when $T$ is symmetrized
5354
OperatorLRHxc<T>* op_ht = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
54-
*this->DM_diff, gint, pot, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
55+
*this->DM_diff, gint, pot_hxc_gs, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
5556
{ 0 }, T(-2.0), ATYPE::CC_vo);
5657
this->ops->add(op_ht);
5758
// 3. $2\sum_{jb,kc} g^{xc}_{ia, jb, kc}X_{jb}X_{kc}$
@@ -60,6 +61,11 @@ namespace LR
6061
*this->DM_trans, gint, this->pot_grad, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
6162
{ 0 }, T(-2.0), ATYPE::CC_vo);
6263
this->ops->add(op_gxc);
64+
// // test: op_ht only
65+
// delete this->ops;
66+
// this->ops = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
67+
// *this->DM_diff, gint, pot_hxc_gs, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
68+
// { 0 }, T(-2.0), ATYPE::CC_vo);
6369
#ifdef __EXX
6470
// add EXX operators here
6571
#endif
@@ -85,6 +91,10 @@ namespace LR
8591
std::vector<ct::Tensor> dm_diff_2d = cal_dm_diff_blas(X, psi_ks_is, naos, nocc[is], nvirt[is]);
8692
#endif
8793
for (int ik = 0;ik < this->nk;++ik) { this->DM_diff->set_DMK_pointer(ik, dm_diff_2d[ik].data<T>()); }
94+
// std::cout << "difference density matrix" << std::endl;
95+
// for (int ik = 0;ik < this->nk;++ik) { LR_Util::print_value(dm_diff_2d[ik].data<T>(), naos, naos); }
96+
// std::cout << "test: set dm_diff to zero" << std::endl;
97+
// for (int ik = 0;ik < this->nk;++ik) { dm_diff_2d[ik].zero(); }
8898
};
8999
}
90100
virtual void hPsi(const T* const psi, T* const hpsi, const int ld_psi, const int nband) const override

‎source/module_lr/Grad/multipliers/zeq_solver.hpp‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -81,6 +81,7 @@ namespace LR
8181
#endif
8282
TGint* gint,
8383
std::weak_ptr<PotHxcLR> pot,
84+
std::weak_ptr<PotHxcLR> pot_hxc_gs,
8485
const K_Vectors& kv,
8586
const std::vector<Parallel_2D>& px,
8687
const Parallel_2D& pc,
@@ -97,7 +98,7 @@ namespace LR
9798
#ifdef __EXX
9899
exx_lri, exx_alpha,
99100
#endif
100-
gint, pot, kv, px, pc, pmat, spin_type);
101+
gint, pot, pot_hxc_gs, kv, px, pc, pmat, spin_type);
101102
ModuleBase::timer::tick("Z_vector", "Z_vector_R");
102103
ops_R.hPsi(X, R.data<T>(), nk * px[0].get_local_size(), nstates); // act each operators on X
103104
ModuleBase::timer::tick("Z_vector", "Z_vector_R");
@@ -109,11 +110,13 @@ namespace LR
109110
#ifdef __EXX
110111
exx_lri, exx_alpha,
111112
#endif
112-
gint, pot, kv, px, pc, pmat, spin_type);
113+
gint, pot_hxc_gs, kv, px, pc, pmat, spin_type);
113114

114115
// 3. solve Z-vector equation
115116
std::cout << "The right side of the Z-vector equation:" << std::endl;
116117
LR_Util::print_value(R.data<T>(), nstates, nloc_per_band);
117118
solve_Z_CG(Z, R.data<T>(), nloc_per_band, nstates, std::bind(&HamiltLR<T>::hPsi, &ops_L, std::placeholders::_1, std::placeholders::_2, nloc_per_band, nstates));
119+
// test: set Z to 0
120+
// for (int i = 0; i < nstates * nloc_per_band; ++i) { Z[i] = T(0.0); }
118121
}
119122
}

0 commit comments

Comments
 (0)