Skip to content

Commit

Permalink
Refactor: runner() of esolver_ks (#5445)
Browse files Browse the repository at this point in the history 
* Refactor: runner() of esolver_ks

* rename hamilt2density and diag as hamilt2density_single and hamilt2density
  • Loading branch information
@YuLiu98
YuLiu98 authored Nov 11, 2024
1 parent 32ba8d4 commit c03c879
Show file tree
Hide file tree
Showing 12 changed files with 332 additions and 403 deletions.
324 changes: 180 additions & 144 deletions source/module_esolver/esolver_ks.cpp
View file

Large diffs are not rendered by default.

8 changes: 6 additions & 2 deletions source/module_esolver/esolver_ks.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,12 +46,15 @@ class ESolver_KS : public ESolver_FP
//! Something to do after hamilt2density function in each iter loop.
virtual void iter_finish(const int istep, int& iter);

// calculate electron density from a specific Hamiltonian
virtual void hamilt2density(const int istep, const int iter, const double ethr);
// calculate electron density from a specific Hamiltonian with ethr
virtual void hamilt2density_single(const int istep, const int iter, const double ethr);

// calculate electron states from a specific Hamiltonian
virtual void hamilt2estates(const double ethr) {};

// calculate electron density from a specific Hamiltonian
void hamilt2density(const int istep, const int iter, const double ethr);

//! Something to do after SCF iterations when SCF is converged or comes to the max iter step.
virtual void after_scf(const int istep) override;

Expand Down Expand Up @@ -82,6 +85,7 @@ class ESolver_KS : public ESolver_FP
double scf_thr; // scf density threshold
double scf_ene_thr; // scf energy threshold
double drho; // the difference between rho_in (before HSolver) and rho_out (After HSolver)
double hsolver_error; // the error of HSolver
int maxniter; // maximum iter steps for scf
int niter; // iter steps actually used in scf
int out_freq_elec; // frequency for output
Expand Down
137 changes: 52 additions & 85 deletions source/module_esolver/esolver_ks_lcao.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -681,10 +681,17 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
SpinConstrain<TK, base_device::DEVICE_CPU>& sc = SpinConstrain<TK, base_device::DEVICE_CPU>::getScInstance();
sc.run_lambda_loop(iter - 1);
}

// save density matrix DMR for mixing
if (PARAM.inp.mixing_restart > 0 && PARAM.inp.mixing_dmr && this->p_chgmix->mixing_restart_count > 0)
{
elecstate::DensityMatrix<TK, double>* dm = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
dm->save_DMR();
}
}

//------------------------------------------------------------------------------
//! the 11th function of ESolver_KS_LCAO: hamilt2density
//! the 11th function of ESolver_KS_LCAO: hamilt2density_single
//! mohan add 2024-05-11
//! 1) save input rho
//! 2) save density matrix DMR for mixing
Expand All @@ -700,47 +707,16 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
//! 12) calculate delta energy
//------------------------------------------------------------------------------
template <typename TK, typename TR>
void ESolver_KS_LCAO<TK, TR>::hamilt2density(int istep, int iter, double ethr)
void ESolver_KS_LCAO<TK, TR>::hamilt2density_single(int istep, int iter, double ethr)
{
ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density");

// 1) save input rho
this->pelec->charge->save_rho_before_sum_band();

// 2) save density matrix DMR for mixing
if (PARAM.inp.mixing_restart > 0 && PARAM.inp.mixing_dmr && this->p_chgmix->mixing_restart_count > 0)
{
elecstate::DensityMatrix<TK, double>* dm = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
dm->save_DMR();
}
ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density_single");

// 3) solve the Hamiltonian and output band gap
{
// reset energy
this->pelec->f_en.eband = 0.0;
this->pelec->f_en.demet = 0.0;

hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, false);
// reset energy
this->pelec->f_en.eband = 0.0;
this->pelec->f_en.demet = 0.0;

if (PARAM.inp.out_bandgap)
{
if (!PARAM.globalv.two_fermi)
{
this->pelec->cal_bandgap();
}
else
{
this->pelec->cal_bandgap_updw();
}
}
}

// 4) print bands for each k-point and each band
for (int ik = 0; ik < this->kv.get_nks(); ++ik)
{
this->pelec->print_band(ik, PARAM.inp.printe, iter);
}
hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, false);

// 5) what's the exd used for?
#ifdef __EXX
Expand All @@ -754,59 +730,13 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2density(int istep, int iter, double ethr)
}
#endif

// 6) calculate the local occupation number matrix and energy correction in
// DFT+U
if (PARAM.inp.dft_plus_u)
{
// only old DFT+U method should calculated energy correction in esolver,
// new DFT+U method will calculate energy in calculating Hamiltonian
if (PARAM.inp.dft_plus_u == 2)
{
if (GlobalC::dftu.omc != 2)
{
const std::vector<std::vector<TK>>& tmp_dm
= dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMK_vector();
this->dftu_cal_occup_m(iter, tmp_dm);
}
GlobalC::dftu.cal_energy_correction(istep);
}
GlobalC::dftu.output();
}

// (7) for deepks, calculate delta_e
#ifdef __DEEPKS
if (PARAM.inp.deepks_scf)
{
const std::vector<std::vector<TK>>& dm
= dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMK_vector();

this->dpks_cal_e_delta_band(dm);
}
#endif

// 8) for delta spin
if (PARAM.inp.sc_mag_switch)
{
SpinConstrain<TK, base_device::DEVICE_CPU>& sc = SpinConstrain<TK, base_device::DEVICE_CPU>::getScInstance();
sc.cal_MW(iter, this->p_hamilt);
}

// 9) use new charge density to calculate energy
this->pelec->cal_energies(1);

// 10) symmetrize the charge density
Symmetry_rho srho;
for (int is = 0; is < PARAM.inp.nspin; is++)
{
srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::ucell.symm);
}

// 11) compute magnetization, only for spin==2
GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx,
this->pelec->charge->nxyz,
this->pelec->charge->rho,
this->pelec->nelec_spin.data());

// 12) calculate delta energy
this->pelec->f_en.deband = this->pelec->cal_delta_eband();
}
Expand Down Expand Up @@ -923,6 +853,43 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(const int istep, int& iter)
{
ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish");

// 6) calculate the local occupation number matrix and energy correction in
// DFT+U
if (PARAM.inp.dft_plus_u)
{
// only old DFT+U method should calculated energy correction in esolver,
// new DFT+U method will calculate energy in calculating Hamiltonian
if (PARAM.inp.dft_plus_u == 2)
{
if (GlobalC::dftu.omc != 2)
{
const std::vector<std::vector<TK>>& tmp_dm
= dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMK_vector();
this->dftu_cal_occup_m(iter, tmp_dm);
}
GlobalC::dftu.cal_energy_correction(istep);
}
GlobalC::dftu.output();
}

// (7) for deepks, calculate delta_e
#ifdef __DEEPKS
if (PARAM.inp.deepks_scf)
{
const std::vector<std::vector<TK>>& dm
= dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMK_vector();

this->dpks_cal_e_delta_band(dm);
}
#endif

// 8) for delta spin
if (PARAM.inp.sc_mag_switch)
{
SpinConstrain<TK, base_device::DEVICE_CPU>& sc = SpinConstrain<TK, base_device::DEVICE_CPU>::getScInstance();
sc.cal_MW(iter, this->p_hamilt);
}

// call iter_finish() of ESolver_KS
ESolver_KS<TK>::iter_finish(istep, iter);

Expand Down
4 changes: 1 addition & 3 deletions source/module_esolver/esolver_ks_lcao.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,9 +50,7 @@ class ESolver_KS_LCAO : public ESolver_KS<TK> {

virtual void iter_init(const int istep, const int iter) override;

virtual void hamilt2density(const int istep,
const int iter,
const double ethr) override;
virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override;

virtual void update_pot(const int istep, const int iter) override;

Expand Down
51 changes: 18 additions & 33 deletions source/module_esolver/esolver_ks_lcao_tddft.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,10 +118,8 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(const Input_para& inp, UnitCell&
this->pelec_td = dynamic_cast<elecstate::ElecStateLCAO_TDDFT*>(this->pelec);
}

void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, const double ethr)
void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int iter, const double ethr)
{
pelec->charge->save_rho_before_sum_band();

if (wf.init_wfc == "file")
{
if (istep >= 1)
Expand Down Expand Up @@ -171,11 +169,23 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons
hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec_td, false);
}
}
// else
// {
// ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "HSolver has not been initialed!");
// }

// symmetrize the charge density only for ground state
if (istep <= 1)
{
Symmetry_rho srho;
for (int is = 0; is < PARAM.inp.nspin; is++)
{
srho.begin(is, *(pelec->charge), pw_rho, GlobalC::ucell.symm);
}
}

// (7) calculate delta energy
this->pelec->f_en.deband = this->pelec->cal_delta_eband();
}

void ESolver_KS_LCAO_TDDFT::iter_finish(const int istep, int& iter)
{
// print occupation of each band
if (iter == 1 && istep <= 2)
{
Expand All @@ -201,32 +211,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons
<< std::endl;
}

for (int ik = 0; ik < kv.get_nks(); ++ik)
{
this->pelec_td->print_band(ik, PARAM.inp.printe, iter);
}

// using new charge density.
this->pelec->cal_energies(1);

// symmetrize the charge density only for ground state
if (istep <= 1)
{
Symmetry_rho srho;
for (int is = 0; is < PARAM.inp.nspin; is++)
{
srho.begin(is, *(pelec->charge), pw_rho, GlobalC::ucell.symm);
}
}

// (6) compute magnetization, only for spin==2
GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx,
this->pelec->charge->nxyz,
this->pelec->charge->rho,
pelec->nelec_spin.data());

// (7) calculate delta energy
this->pelec->f_en.deband = this->pelec->cal_delta_eband();
ESolver_KS_LCAO<std::complex<double>, double>::iter_finish(istep, iter);
}

void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)
Expand Down
4 changes: 3 additions & 1 deletion source/module_esolver/esolver_ks_lcao_tddft.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,10 +30,12 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO<std::complex<double>, doubl
int td_htype = 1;

protected:
virtual void hamilt2density(const int istep, const int iter, const double ethr) override;
virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override;

virtual void update_pot(const int istep, const int iter) override;

virtual void iter_finish(const int istep, int& iter) override;

virtual void after_scf(const int istep) override;

void cal_edm_tddft();
Expand Down
Loading

0 comments on commit c03c879

Please sign in to comment.