Skip to content

Commit e493841

Browse files
mohanchenmohanchen
committed
remove old psiToRho function in elecstate_lcao and replace with new rho_tau_lcao
1 parent b8e0d85 commit e493841

File tree

11 files changed

+53
-112
lines changed

11 files changed

+53
-112
lines changed

‎source/source_esolver/esolver_dm2rho.cpp‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@
1010
#include "source_io/cube_io.h"
1111
#include "source_io/io_npz.h"
1212
#include "source_io/print_info.h"
13+
#include "source_lcao/rho_tau_lcao.h" // mohan add 2025-10-24
1314

1415
namespace ModuleESolver
1516
{
@@ -45,9 +46,16 @@ void ESolver_DM2rho<TK, TR>::runner(UnitCell& ucell, const int istep)
4546

4647
ESolver_KS_LCAO<TK, TR>::before_scf(ucell, istep);
4748

49+
auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
50+
51+
if(!estate)
52+
{
53+
ModuleBase::WARNING_QUIT("ESolver_DM2rho::after_scf","pelec does not exist");
54+
}
55+
4856
// file name of DM
4957
std::string zipname = "output_DM0.npz";
50-
elecstate::DensityMatrix<TK, double>* dm = dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
58+
elecstate::DensityMatrix<TK, double>* dm = estate->get_DM();
5159

5260
// read DM from file
5361
ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(1)));
@@ -59,7 +67,9 @@ void ESolver_DM2rho<TK, TR>::runner(UnitCell& ucell, const int istep)
5967
ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(2)));
6068
}
6169

62-
this->pelec->psiToRho(*this->psi);
70+
// it's dangrous to design psiToRho function like this, mohan note 20251024
71+
// this->pelec->psiToRho(*this->psi);
72+
LCAO_domain::dm2rho(estate->DM->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
6373

6474
int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1;
6575

‎source/source_esolver/esolver_ks_lcao.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -469,7 +469,7 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2rho_single(UnitCell& ucell, int istep, int
469469
if (!skip_solve)
470470
{
471471
hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
472-
hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, skip_charge);
472+
hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, this->chr, PARAM.inp.nspin, skip_charge);
473473
}
474474

475475
// 4) EXX

‎source/source_esolver/esolver_ks_lcao_tddft.cpp‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -31,10 +31,9 @@
3131
#include "source_lcao/hamilt_lcao.h"
3232
#include "source_psi/psi.h"
3333

34-
//-----force& stress-------------------
3534
#include "source_lcao/FORCE_STRESS.h"
35+
#include "source_lcao/rho_tau_lcao.h" // mohan add 2025-10-24
3636

37-
//---------------------------------------------------
3837

3938
namespace ModuleESolver
4039
{
@@ -293,7 +292,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::hamilt2rho_single(UnitCell& ucell,
293292
{
294293
bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false;
295294
hsolver::HSolverLCAO<std::complex<double>> hsolver_lcao_obj(&this->pv, PARAM.inp.ks_solver);
296-
hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, skip_charge);
295+
hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, this->chr, PARAM.inp.nspin, skip_charge);
297296
}
298297
}
299298

@@ -574,8 +573,8 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::weight_dm_rho(const UnitCell& ucell)
574573
_pes->DM->cal_DMR();
575574
}
576575

577-
// get the real-space charge density
578-
this->pelec->psiToRho(this->psi[0]);
576+
// get the real-space charge density, mohan add 2025-10-24
577+
LCAO_domain::dm2rho(_pes->DM->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
579578
}
580579

581580
template class ESolver_KS_LCAO_TDDFT<double, base_device::DEVICE_CPU>;

‎source/source_estate/elecstate_lcao.cpp‎

Lines changed: 4 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -15,68 +15,6 @@
1515
namespace elecstate
1616
{
1717

18-
// multi-k case
19-
template <>
20-
void ElecStateLCAO<std::complex<double>>::psiToRho(const psi::Psi<std::complex<double>>& psi)
21-
{
22-
ModuleBase::TITLE("ElecStateLCAO", "psiToRho");
23-
ModuleBase::timer::tick("ElecStateLCAO", "psiToRho");
24-
25-
for (int is = 0; is < PARAM.inp.nspin; is++)
26-
{
27-
ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is],
28-
this->charge->nrxx); // mohan 2009-11-10
29-
}
30-
31-
//------------------------------------------------------------
32-
// calculate the charge density on real space grid.
33-
//------------------------------------------------------------
34-
35-
ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
36-
ModuleGint::cal_gint_rho(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->rho);
37-
38-
if (XC_Functional::get_ked_flag())
39-
{
40-
this->cal_tau(psi);
41-
}
42-
43-
this->charge->renormalize_rho();
44-
45-
ModuleBase::timer::tick("ElecStateLCAO", "psiToRho");
46-
return;
47-
}
48-
49-
// Gamma_only case
50-
template <>
51-
void ElecStateLCAO<double>::psiToRho(const psi::Psi<double>& psi)
52-
{
53-
ModuleBase::TITLE("ElecStateLCAO", "psiToRho");
54-
ModuleBase::timer::tick("ElecStateLCAO", "psiToRho");
55-
56-
for (int is = 0; is < PARAM.inp.nspin; is++)
57-
{
58-
ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is],
59-
this->charge->nrxx); // mohan 2009-11-10
60-
}
61-
62-
//------------------------------------------------------------
63-
// calculate the charge density on real space grid.
64-
//------------------------------------------------------------
65-
ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
66-
67-
ModuleGint::cal_gint_rho(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->rho);
68-
69-
if (XC_Functional::get_ked_flag())
70-
{
71-
this->cal_tau(psi);
72-
}
73-
74-
this->charge->renormalize_rho();
75-
76-
ModuleBase::timer::tick("ElecStateLCAO", "psiToRho");
77-
return;
78-
}
79-
8018
template <typename TK>
8119
void ElecStateLCAO<TK>::init_DM(const K_Vectors* kv, const Parallel_Orbitals* paraV, const int nspin)
8220
{
@@ -101,9 +39,9 @@ double ElecStateLCAO<std::complex<double>>::get_spin_constrain_energy()
10139

10240
#ifdef __PEXSI
10341
template <>
104-
void ElecStateLCAO<double>::dmToRho(std::vector<double*> pexsi_DM, std::vector<double*> pexsi_EDM)
42+
void ElecStateLCAO<double>::dm2Rho(std::vector<double*> pexsi_DM, std::vector<double*> pexsi_EDM)
10543
{
106-
ModuleBase::timer::tick("ElecStateLCAO", "dmToRho");
44+
ModuleBase::timer::tick("ElecStateLCAO", "dm2Rho");
10745

10846
int nspin = PARAM.inp.nspin;
10947
if (PARAM.inp.nspin == 4)
@@ -138,12 +76,12 @@ void ElecStateLCAO<double>::dmToRho(std::vector<double*> pexsi_DM, std::vector<d
13876

13977
this->charge->renormalize_rho();
14078

141-
ModuleBase::timer::tick("ElecStateLCAO", "dmToRho");
79+
ModuleBase::timer::tick("ElecStateLCAO", "dm2Rho");
14280
return;
14381
}
14482

14583
template <>
146-
void ElecStateLCAO<std::complex<double>>::dmToRho(std::vector<std::complex<double>*> pexsi_DM,
84+
void ElecStateLCAO<std::complex<double>>::dm2Rho(std::vector<std::complex<double>*> pexsi_DM,
14785
std::vector<std::complex<double>*> pexsi_EDM)
14886
{
14987
ModuleBase::WARNING_QUIT("ElecStateLCAO", "pexsi is not completed for multi-k case");

‎source/source_estate/elecstate_lcao.h‎

Lines changed: 1 addition & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -33,15 +33,6 @@ class ElecStateLCAO : public ElecState
3333
}
3434
}
3535

36-
// void init(Charge* chg_in):charge(chg_in){} override;
37-
38-
// interface for HSolver to calculate rho from Psi
39-
virtual void psiToRho(const psi::Psi<TK>& psi) override;
40-
// virtual void psiToRho(const psi::Psi<double>& psi) override;
41-
// return current electronic density rho, as a input for constructing Hamiltonian
42-
// const double* getRho(int spin) const override;
43-
virtual void cal_tau(const psi::Psi<TK>& psi) override;
44-
4536
// update charge density for next scf step
4637
// void getNewRho() override;
4738

@@ -65,20 +56,11 @@ class ElecStateLCAO : public ElecState
6556
* @param pexsi_EDM: pointers of energy-weighed density matrix (EDMK) calculated by pexsi, needed by MD, will be
6657
* stored in DensityMatrix::pexsi_EDM
6758
*/
68-
void dmToRho(std::vector<TK*> pexsi_DM, std::vector<TK*> pexsi_EDM);
59+
void dm2Rho(std::vector<TK*> pexsi_DM, std::vector<TK*> pexsi_EDM);
6960
#endif
7061

7162
DensityMatrix<TK, double>* DM = nullptr;
7263

73-
protected:
74-
// calculate electronic charge density on grid points or density matrix in real space
75-
// the consequence charge density rho saved into rho_out, preparing for charge mixing.
76-
// void updateRhoK(const psi::Psi<std::complex<double>>& psi) ;//override;
77-
// sum over all pools for rho and ebands
78-
// void parallelK();
79-
// calcualte rho for each k
80-
// void rhoBandK(const psi::Psi<std::complex<double>>& psi);
81-
8264
};
8365

8466
template <typename TK>

‎source/source_estate/elecstate_lcao_cal_tau.cpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
namespace elecstate
77
{
88

9+
/*
910
// calculate the kinetic energy density tau, multi-k case
1011
template <>
1112
void ElecStateLCAO<std::complex<double>>::cal_tau(const psi::Psi<std::complex<double>>& psi)
@@ -36,4 +37,6 @@ void ElecStateLCAO<double>::cal_tau(const psi::Psi<double>& psi)
3637
ModuleBase::timer::tick("ElecStateLCAO", "cal_tau");
3738
return;
3839
}
39-
}
40+
*/
41+
42+
}

‎source/source_hsolver/hsolver_lcao.cpp‎

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -34,13 +34,17 @@
3434
#include "source_hsolver/parallel_k2d.h"
3535
#include "source_io/module_parameter/parameter.h"
3636

37+
#include "source_lcao/rho_tau_lcao.h" // mohan add 20251024
38+
3739
namespace hsolver
3840
{
3941

4042
template <typename T, typename Device>
4143
void HSolverLCAO<T, Device>::solve(hamilt::Hamilt<T>* pHamilt,
4244
psi::Psi<T>& psi,
4345
elecstate::ElecState* pes,
46+
Charge &chr,
47+
const int nspin,
4448
const bool skip_charge)
4549
{
4650
ModuleBase::TITLE("HSolverLCAO", "solve");
@@ -97,9 +101,8 @@ void HSolverLCAO<T, Device>::solve(hamilt::Hamilt<T>* pHamilt,
97101

98102
if (!skip_charge)
99103
{
100-
// used in scf calculation
101-
// calculate charge by eigenpairs(eigenvalues and eigenvectors)
102-
pes->psiToRho(psi);
104+
// compute charge density from density matrix, mohan update 20251024
105+
LCAO_domain::dm2rho(_pes_lcao->DM->get_DMR_vector(), nspin, &chr);
103106
}
104107
else
105108
{
@@ -492,4 +495,4 @@ void HSolverLCAO<T, Device>::parakSolve_cusolver(hamilt::Hamilt<T>* pHamilt,
492495
template class HSolverLCAO<double>;
493496
template class HSolverLCAO<std::complex<double>>;
494497

495-
} // namespace hsolver
498+
} // namespace hsolver

‎source/source_hsolver/hsolver_lcao.h‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,8 @@
55
#include "source_hamilt/hamilt.h"
66
#include "source_basis/module_ao/parallel_orbitals.h"
77

8+
#include "source_estate/module_charge/charge.h" // mohan add 20251024
9+
810
namespace hsolver
911
{
1012

@@ -17,7 +19,9 @@ class HSolverLCAO
1719
void solve(hamilt::Hamilt<T>* pHamilt,
1820
psi::Psi<T>& psi,
1921
elecstate::ElecState* pes,
20-
const bool skip_charge);
22+
Charge &chr, // charge density
23+
const int nspin,
24+
const bool skip_charge);
2125

2226
private:
2327
void hamiltSolvePsiK(hamilt::Hamilt<T>* hm, psi::Psi<T>& psi, double* eigenvalue); // for kpar_lcao == 1
@@ -36,4 +40,4 @@ class HSolverLCAO
3640

3741
} // namespace hsolver
3842

39-
#endif
43+
#endif

‎source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -151,7 +151,8 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(int
151151
->update_lambda();
152152
}
153153
// diagonalization without update charge
154-
hsolver_t.solve(hamilt_t, psi_t[0], this->pelec, true);
154+
// mohan add two parameters charge and nspin, 2025-10-24
155+
hsolver_t.solve(hamilt_t, psi_t[0], this->pelec, *this->pelec->charge, PARAM.inp.nspin, true);
155156
elecstate::calculate_weights(this->pelec->ekb,
156157
this->pelec->wg,
157158
this->pelec->klist,

‎source/source_lcao/rho_tau_lcao.cpp‎

Lines changed: 10 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -4,19 +4,22 @@
44

55
void LCAO_domain::dm2rho(std::vector<hamilt::HContainer<double>*> &dmr,
66
const int nspin,
7-
Charge* chg)
7+
Charge* chr)
88
{
99
ModuleBase::TITLE("LCAO_domain", "dm2rho");
1010
ModuleBase::timer::tick("LCAO_domain", "dm2rho");
1111

1212
for (int is = 0; is < nspin; is++)
1313
{
14-
ModuleBase::GlobalFunc::ZEROS(chg->rho[is], chg->nrxx);
14+
ModuleBase::GlobalFunc::ZEROS(chr->rho[is], chr->nrxx);
1515
}
1616

17-
ModuleGint::cal_gint_rho(dmr, nspin, chg->rho);
17+
ModuleGint::cal_gint_rho(dmr, nspin, chr->rho);
1818

19-
chg->renormalize_rho();
19+
chr->renormalize_rho();
20+
21+
// should be moved somewhere else, mohan 20251024
22+
dm2tau(dmr, nspin, chr);
2023

2124
// symmetrize of charge density should be here, mohan 20251023
2225

@@ -27,7 +30,7 @@ void LCAO_domain::dm2rho(std::vector<hamilt::HContainer<double>*> &dmr,
2730

2831
void LCAO_domain::dm2tau(std::vector<hamilt::HContainer<double>*> &dmr,
2932
const int nspin,
30-
Charge* chg)
33+
Charge* chr)
3134
{
3235
ModuleBase::TITLE("LCAO_domain", "dm2tau");
3336
ModuleBase::timer::tick("LCAO_domain", "dm2tau");
@@ -36,12 +39,10 @@ void LCAO_domain::dm2tau(std::vector<hamilt::HContainer<double>*> &dmr,
3639
{
3740
for (int is = 0; is < nspin; is++)
3841
{
39-
ModuleBase::GlobalFunc::ZEROS(chg->kin_r[is], chg->nrxx);
42+
ModuleBase::GlobalFunc::ZEROS(chr->kin_r[is], chr->nrxx);
4043
}
41-
ModuleGint::cal_gint_tau(dmr, nspin, chg->kin_r);
44+
ModuleGint::cal_gint_tau(dmr, nspin, chr->kin_r);
4245
}
4346

4447
ModuleBase::timer::tick("LCAO_domain", "dm2tau");
4548
}
46-
47-

0 commit comments

Comments
 (0)