Skip to content

Commit

Permalink
Refactor: delete wg and ekb from class wavefunc, only elecstate::Elec…
Browse files Browse the repository at this point in the history 
…State have these two members now (#1478)

* Refactor: removed GlobalC::wf.wg and GlobalC::wf.ekb

* Fix: reset pelec_td from unique pointer to another name of pelec

* Fix: delete redundant wf.wg and wf.ekb

* Refactor: delete wg and ekb from class wavefunc

Co-authored-by: root <zhengdy@bjaisi.com>
  • Loading branch information
@dyzheng
dyzheng and root committed Nov 7, 2022
1 parent f145042 commit 100fd74
Show file tree
Hide file tree
Showing 100 changed files with 487 additions and 4,174 deletions.
85 changes: 43 additions & 42 deletions source/module_elecstate/elecstate.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -215,52 +215,53 @@ void ElecState::calEBand()

void ElecState::print_band(const int& ik, const int& printe, const int& iter)
{
// check the band energy.
//check the band energy.
bool wrong = false;
int nbands = this->ekb.nc;
for (int ib = 0; ib < nbands; ++ib)
{
if (abs(this->ekb(ik, ib)) > 1.0e10)
{
GlobalV::ofs_warning << " ik=" << ik + 1 << " ib=" << ib + 1 << " " << this->ekb(ik, ib) << " Ry"
<< std::endl;
wrong = true;
}
}
if (wrong)
for(int ib=0; ib<GlobalV::NBANDS; ++ib)
{
if( abs( this->ekb(ik, ib) ) > 1.0e10)
{
GlobalV::ofs_warning << " ik=" << ik+1 << " ib=" << ib+1 << " " << this->ekb(ik, ib) << " Ry" << std::endl;
wrong = true;
}
}
if(wrong)
{
ModuleBase::WARNING_QUIT("Threshold_Elec::print_eigenvalue", "Eigenvalues are too large!");
ModuleBase::WARNING_QUIT("Threshold_Elec::print_eigenvalue","Eigenvalues are too large!");
}

if (GlobalV::MY_RANK == 0)
{
// if( GlobalV::DIAGO_TYPE == "selinv" ) xiaohui modify 2013-09-02
if (GlobalV::KS_SOLVER == "selinv") // xiaohui add 2013-09-02
{
GlobalV::ofs_running << " No eigenvalues are available for selected inversion methods." << std::endl;
}
else
{
if (printe > 0 && ((iter + 1) % printe == 0))
{
// NEW_PART("ENERGY BANDS (Rydberg), (eV)");
GlobalV::ofs_running << std::setprecision(6);
GlobalV::ofs_running << " Energy (eV) & Occupations for spin=" << GlobalV::CURRENT_SPIN + 1
<< " K-point=" << ik + 1 << std::endl;
GlobalV::ofs_running << std::setiosflags(ios::showpoint);
for (int ib = 0; ib < nbands; ib++)
{
GlobalV::ofs_running << " " << std::setw(6) << ib + 1 << std::setw(15)
<< this->ekb(ik, ib) * ModuleBase::Ry_to_eV;
// for the first electron iteration, we don't have the energy
// spectrum, so we can't get the occupations.
GlobalV::ofs_running << std::setw(15) << this->wg(ik, ib);
GlobalV::ofs_running << std::endl;
}
}
}
}
return;


if(GlobalV::MY_RANK==0)
{
//if( GlobalV::DIAGO_TYPE == "selinv" ) xiaohui modify 2013-09-02
if(GlobalV::KS_SOLVER=="selinv") //xiaohui add 2013-09-02
{
GlobalV::ofs_running << " No eigenvalues are available for selected inversion methods." << std::endl;
}
else
{
if( printe>0 && ((iter+1) % printe == 0))
{
// NEW_PART("ENERGY BANDS (Rydberg), (eV)");


GlobalV::ofs_running << std::setprecision(6);
GlobalV::ofs_running << " Energy (eV) & Occupations for spin=" << GlobalV::CURRENT_SPIN+1 << " K-point=" << ik+1 << std::endl;
GlobalV::ofs_running << std::setiosflags(ios::showpoint);
for(int ib=0;ib<GlobalV::NBANDS;ib++)
{
GlobalV::ofs_running << " "<< std::setw(6) << ib+1
<< std::setw(15) << this->ekb(ik, ib) * ModuleBase::Ry_to_eV;
// for the first electron iteration, we don't have the energy
// spectrum, so we can't get the occupations.
GlobalV::ofs_running << std::setw(15) << this->wg(ik,ib);
GlobalV::ofs_running << std::endl;
}
}
}
}
return;
}

} // namespace elecstate
6 changes: 3 additions & 3 deletions source/module_elecstate/elecstate.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,12 +81,12 @@ class ElecState

std::string classname = "none";

// print and check for band energy and occupations
void print_band(const int &ik, const int &printe, const int &iter);

protected:
// calculate ebands for all k points and all occupied bands
void calEBand();

// print and check for band energy and occupations
void print_band(const int &ik, const int &printe, const int &iter);
};

} // namespace elecstate
Expand Down
7 changes: 3 additions & 4 deletions source/module_elecstate/test/elecstate_lcao_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,8 +89,8 @@ namespace GlobalC

namespace WF_Local
{
int read_lowf(double** ctot, const int& is, const Parallel_Orbitals* ParaV, psi::Psi<double>*) {return 1;};
int read_lowf_complex(std::complex<double>** ctot, const int& ik, const Parallel_Orbitals* ParaV, psi::Psi<std::complex<double> >*) {return 1;}
int read_lowf(double** ctot, const int& is, const Parallel_Orbitals* ParaV, psi::Psi<double>*, elecstate::ElecState*) {return 1;};
int read_lowf_complex(std::complex<double>** ctot, const int& ik, const Parallel_Orbitals* ParaV, psi::Psi<std::complex<double> >*, elecstate::ElecState*) {return 1;}
void write_lowf(const std::string &name, double **ctot, const ModuleBase::matrix& ekb, const ModuleBase::matrix& wg) {}
void write_lowf_complex(const std::string &name, std::complex<double>** ctot, const int &ik, const ModuleBase::matrix& ekb, const ModuleBase::matrix& wg) {}
}
Expand Down Expand Up @@ -309,7 +309,6 @@ class ElecStateLCAOPrepare
void set_env()
{
GlobalV::NBANDS = nbands;
GlobalC::wf.wg = this->wg;

GlobalC::kv.nks = GlobalC::kv.nkstot = nk;
GlobalC::kv.isk.resize(nk,0);
Expand Down Expand Up @@ -505,7 +504,7 @@ class ElecStateLCAOPrepare
{
psik = (psi::Psi<std::complex<double>>*)(&(this->psi));
}
loc.allocate_dm_wfc(GlobalC::GridT.lgd, lowf, psigo, psik);
loc.allocate_dm_wfc(GlobalC::GridT.lgd, nullptr, lowf, psigo, psik);

elecstate::MockElecStateLCAO mesl(&GlobalC::CHR,&GlobalC::kv,nk,nbands,&loc,&uhm,&lowf,this->wg);

Expand Down
15 changes: 0 additions & 15 deletions source/module_elecstate/test/updaterhok_pw_test.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,29 +115,14 @@ Restart restart;
} // namespace GlobalC
Input INPUT;

/*
void Occupy::calculate_weights()
{
GlobalC::wf.wg(0,0)=2.0;
GlobalC::wf.wg(0,1)=0.0;
GlobalC::wf.wg(0,2)=0.0;
GlobalC::wf.wg(0,3)=0.0;
}
*/

void Restart::load_disk(const std::string mode, const int i) const {}


psi::Psi<complex<double>>* wavefunc::allocate(const int nks)
{
this->npwx = GlobalC::wfcpw->npwk_max;
this->wg.create(nks,GlobalV::NBANDS);
this->ekb = new double*[nks];
psi::Psi<std::complex<double>>* psi = new psi::Psi<std::complex<double>>(nks, GlobalV::NBANDS,npwx, nullptr);
for (int ik=0;ik<nks;ik++)
{
this->ekb[ik] = new double[GlobalV::NBANDS];
}
return psi;
}

Expand Down
67 changes: 19 additions & 48 deletions source/module_esolver/esolver_ks_lcao.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,12 +186,6 @@ void ESolver_KS_LCAO::Init(Input& inp, UnitCell_pseudo& ucell)
}
#endif

// Initialize the local wave functions.
// npwx, eigenvalues, and weights
// npwx may change according to cell change
// this function belongs to cell LOOP
GlobalC::wf.allocate_ekb_wg(GlobalC::kv.nks);

// Initialize the FFT.
// this function belongs to cell LOOP

Expand Down Expand Up @@ -222,6 +216,7 @@ void ESolver_KS_LCAO::cal_Force(ModuleBase::matrix& force)
GlobalV::TEST_FORCE,
GlobalV::TEST_STRESS,
this->LOC,
this->pelec,
this->psid,
this->psi,
this->UHM,
Expand Down Expand Up @@ -251,7 +246,7 @@ void ESolver_KS_LCAO::postprocess()
GlobalV::ofs_running << " !FINAL_ETOT_IS " << GlobalC::en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl;
GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl;

GlobalC::en.perform_dos(this->psid, this->psi, this->UHM);
GlobalC::en.perform_dos(this->psid, this->psi, this->UHM, this->pelec);
}

void ESolver_KS_LCAO::Init_Basis_lcao(ORB_control& orb_con, Input& inp, UnitCell_pseudo& ucell)
Expand Down Expand Up @@ -347,32 +342,16 @@ void ESolver_KS_LCAO::eachiterinit(const int istep, const int iter)
// mohan move it outside 2011-01-13
// first need to calculate the weight according to
// electrons number.
// mohan add iter > 1 on 2011-04-02
// because the GlobalC::en.ekb has not value now.
// so the smearing can not be done.
// if (iter > 1)Occupy::calculate_weights();

if (GlobalC::wf.init_wfc == "file")
{
if (iter == 1)
{
std::cout << " WAVEFUN -> CHARGE " << std::endl;

// The occupation should be read in together.
// Occupy::calculate_weights(); //mohan add 2012-02-15

// calculate the density matrix using read in wave functions
// and the ncalculate the charge density on grid.

// transform wg and ekb to elecstate first
for (int ik = 0; ik < this->pelec->ekb.nr; ++ik)
{
for (int ib = 0; ib < this->pelec->ekb.nc; ++ib)
{
this->pelec->ekb(ik, ib) = GlobalC::wf.ekb[ik][ib];
this->pelec->wg(ik, ib) = GlobalC::wf.wg(ik, ib);
}
}
if (this->psi != nullptr)
{
this->pelec->psiToRho(this->psi[0]);
Expand Down Expand Up @@ -479,7 +458,7 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr)
// print ekb for each k point and each band
for (int ik = 0; ik < GlobalC::kv.nks; ++ik)
{
GlobalC::en.print_band(ik);
this->pelec->print_band(ik, GlobalC::en.printe, iter);
}

#ifdef __MPI
Expand Down Expand Up @@ -633,14 +612,6 @@ void ESolver_KS_LCAO::eachiterfinish(int iter)

void ESolver_KS_LCAO::afterscf(const int istep)
{
for (int ik = 0; ik < this->pelec->ekb.nr; ++ik)
{
for (int ib = 0; ib < this->pelec->ekb.nc; ++ib)
{
GlobalC::wf.ekb[ik][ib] = this->pelec->ekb(ik, ib);
GlobalC::wf.wg(ik, ib) = this->pelec->wg(ik, ib);
}
}
// if (this->conv_elec || iter == GlobalV::SCF_NMAX)
// {
//--------------------------------------
Expand All @@ -655,7 +626,7 @@ void ESolver_KS_LCAO::afterscf(const int istep)
std::cout << "dim = " << GlobalC::chi0_hilbert.dim << std::endl;
// std::cout <<"oband = "<<GlobalC::chi0_hilbert.oband<<std::endl;
GlobalC::chi0_hilbert.wfc_k_grid = this->LOWF.wfc_k_grid;
GlobalC::chi0_hilbert.Chi();
GlobalC::chi0_hilbert.Chi(this->pelec->ekb);
}

for (int is = 0; is < GlobalV::NSPIN; is++)
Expand Down Expand Up @@ -705,7 +676,7 @@ void ESolver_KS_LCAO::afterscf(const int istep)

if (GlobalV::OUT_LEVEL != "m")
{
Threshold_Elec::print_eigenvalue(GlobalV::ofs_running);
Threshold_Elec::print_eigenvalue(GlobalV::ofs_running, this->pelec);
}

if (this->conv_elec)
Expand Down Expand Up @@ -739,22 +710,22 @@ void ESolver_KS_LCAO::afterscf(const int istep)
{
if (GlobalV::GAMMA_ONLY_LOCAL)
{
GlobalC::ld.save_npy_o(GlobalC::wf.ekb[0][nocc] - GlobalC::wf.ekb[0][nocc - 1], "o_tot.npy");
GlobalC::ld.save_npy_o(this->pelec->ekb(0, nocc) - this->pelec->ekb(0, nocc - 1), "o_tot.npy");
}
else
{
double homo = GlobalC::wf.ekb[0][nocc - 1];
double lumo = GlobalC::wf.ekb[0][nocc];
double homo = this->pelec->ekb(0, nocc - 1);
double lumo = this->pelec->ekb(0, nocc);
for (int ik = 1; ik < GlobalC::kv.nks; ik++)
{
if (homo < GlobalC::wf.ekb[ik][nocc - 1])
if (homo < this->pelec->ekb(ik, nocc - 1))
{
homo = GlobalC::wf.ekb[ik][nocc - 1];
homo = this->pelec->ekb(ik, nocc - 1);
GlobalC::ld.h_ind = ik;
}
if (lumo > GlobalC::wf.ekb[ik][nocc])
if (lumo > this->pelec->ekb(ik, nocc))
{
lumo = GlobalC::wf.ekb[ik][nocc];
lumo = this->pelec->ekb(ik, nocc);
GlobalC::ld.l_ind = ik;
}
}
Expand Down Expand Up @@ -806,7 +777,7 @@ void ESolver_KS_LCAO::afterscf(const int istep)
GlobalC::ld.save_npy_orbital_precalc(GlobalC::ucell.nat);

GlobalC::ld.cal_o_delta(dm_bandgap_gamma, *this->LOWF.ParaV);
GlobalC::ld.save_npy_o(GlobalC::wf.ekb[0][nocc] - GlobalC::wf.ekb[0][nocc - 1]
GlobalC::ld.save_npy_o(this->pelec->ekb(0, nocc) - this->pelec->ekb(0, nocc - 1)
- GlobalC::ld.o_delta,
"o_base.npy");
}
Expand Down Expand Up @@ -842,8 +813,8 @@ void ESolver_KS_LCAO::afterscf(const int istep)
GlobalC::ld.save_npy_orbital_precalc(GlobalC::ucell.nat);

GlobalC::ld.cal_o_delta_k(dm_bandgap_k, *this->LOWF.ParaV, GlobalC::kv.nks);
GlobalC::ld.save_npy_o(GlobalC::wf.ekb[GlobalC::ld.l_ind][nocc]
- GlobalC::wf.ekb[GlobalC::ld.h_ind][nocc - 1] - GlobalC::ld.o_delta,
GlobalC::ld.save_npy_o(this->pelec->ekb(GlobalC::ld.l_ind, nocc)
- this->pelec->ekb(GlobalC::ld.h_ind, nocc - 1) - GlobalC::ld.o_delta,
"o_base.npy");
}
}
Expand All @@ -855,13 +826,13 @@ void ESolver_KS_LCAO::afterscf(const int istep)
{
if (GlobalV::GAMMA_ONLY_LOCAL)
{
GlobalC::ld.save_npy_o(GlobalC::wf.ekb[0][nocc] - GlobalC::wf.ekb[0][nocc - 1],
GlobalC::ld.save_npy_o(this->pelec->ekb(0, nocc) - this->pelec->ekb(0, nocc - 1),
"o_base.npy"); // no scf, o_tot=o_base
}
else
{
GlobalC::ld.save_npy_o(GlobalC::wf.ekb[GlobalC::ld.l_ind][nocc]
- GlobalC::wf.ekb[GlobalC::ld.h_ind][nocc - 1],
GlobalC::ld.save_npy_o(this->pelec->ekb(GlobalC::ld.l_ind, nocc)
- this->pelec->ekb(GlobalC::ld.h_ind, nocc - 1),
"o_base.npy");
}
}
Expand Down Expand Up @@ -932,7 +903,7 @@ void ESolver_KS_LCAO::afterscf(const int istep)
{
ModuleRPA::DFT_RPA_interface rpa_interface(GlobalC::exx_global.info);
rpa_interface.rpa_exx_lcao().info.files_abfs = GlobalV::rpa_orbitals;
rpa_interface.out_for_RPA(*(this->LOWF.ParaV), *(this->psi), this->LOC);
rpa_interface.out_for_RPA(*(this->LOWF.ParaV), *(this->psi), this->LOC, this->pelec);
}
if (hsolver::HSolverLCAO::out_mat_hsR)
{
Expand Down

0 comments on commit 100fd74

Please sign in to comment.