Skip to content

Refactor: runner() of esolver_ks #5445

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
Nov 11, 2024
Merged

Refactor: runner() of esolver_ks #5445

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
398915e
Refactor: runner() of esolver_ks
YuLiu98 Nov 8, 2024
cb7fc58
Merge branch 'develop' of https://github.com/deepmodeling/abacus-deve…
YuLiu98 Nov 9, 2024
65815c4
rename hamilt2density and diag as hamilt2density_single and hamilt2de…
YuLiu98 Nov 10, 2024
7d5807d
Merge branch 'develop' of https://github.com/deepmodeling/abacus-deve…
YuLiu98 Nov 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
324 changes: 180 additions & 144 deletions source/module_esolver/esolver_ks.cpp
View file
Open in desktop

Large diffs are not rendered by default.

8 changes: 6 additions & 2 deletions source/module_esolver/esolver_ks.h
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading