Skip to content

Refactor: remove set_matrix_grid #5558

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 1 commit into from
Nov 21, 2024
Merged

Refactor: remove set_matrix_grid #5558

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
1 change: 0 additions & 1 deletion source/Makefile.Objects
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -251,7 +251,6 @@ OBJS_ESOLVER=esolver.o\
OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
esolver_ks_lcao_tddft.o\
dpks_cal_e_delta_band.o\
set_matrix_grid.o\
lcao_before_scf.o\
esolver_gets.o\
lcao_others.o\
Expand Down
1 change: 0 additions & 1 deletion source/module_esolver/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,6 @@ if(ENABLE_LCAO)
esolver_ks_lcao.cpp
esolver_ks_lcao_tddft.cpp
dpks_cal_e_delta_band.cpp
set_matrix_grid.cpp
lcao_before_scf.cpp
esolver_gets.cpp
lcao_others.cpp
Expand Down
150 changes: 97 additions & 53 deletions source/module_esolver/lcao_before_scf.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,13 +39,94 @@ namespace ModuleESolver
{

template <typename TK, typename TR>
void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
void ESolver_KS_LCAO<TK, TR>::before_scf(const int istep)
{
ModuleBase::TITLE("ESolver_KS_LCAO", "beforesolver");
ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver");
ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf");

//! 1) call before_scf() of ESolver_FP
ESolver_FP::before_scf(istep);

if (GlobalC::ucell.ionic_position_updated)
{
this->CE.update_all_dis(GlobalC::ucell);
this->CE.extrapolate_charge(
#ifdef __MPI
&(GlobalC::Pgrid),
#endif
GlobalC::ucell,
this->pelec->charge,
&(this->sf),
GlobalV::ofs_running,
GlobalV::ofs_warning);
}

//----------------------------------------------------------
// about vdw, jiyy add vdwd3 and linpz add vdwd2
//----------------------------------------------------------
auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp, &(GlobalV::ofs_running));
if (vdw_solver != nullptr)
{
this->pelec->f_en.evdw = vdw_solver->get_energy();
}

// 1. prepare HS matrices, prepare grid integral
this->set_matrix_grid(this->RA);
// (1) Find adjacent atoms for each atom.
double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running,
PARAM.inp.out_level,
orb_.get_rcutmax_Phi(),
GlobalC::ucell.infoNL.get_rcutmax_Beta(),
PARAM.globalv.gamma_only_local);

atom_arrange::search(PARAM.inp.search_pbc,
GlobalV::ofs_running,
GlobalC::GridD,
GlobalC::ucell,
search_radius,
PARAM.inp.test_atom_input);

// (3) Periodic condition search for each grid.
double dr_uniform = 0.001;
std::vector<double> rcuts;
std::vector<std::vector<double>> psi_u;
std::vector<std::vector<double>> dpsi_u;
std::vector<std::vector<double>> d2psi_u;

Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, orb_, psi_u, dpsi_u, d2psi_u);

this->GridT.set_pbc_grid(this->pw_rho->nx,
this->pw_rho->ny,
this->pw_rho->nz,
this->pw_big->bx,
this->pw_big->by,
this->pw_big->bz,
this->pw_big->nbx,
this->pw_big->nby,
this->pw_big->nbz,
this->pw_big->nbxx,
this->pw_big->nbzp_start,
this->pw_big->nbzp,
this->pw_rho->ny,
this->pw_rho->nplane,
this->pw_rho->startz_current,
GlobalC::ucell,
GlobalC::GridD,
dr_uniform,
rcuts,
psi_u,
dpsi_u,
d2psi_u,
PARAM.inp.nstream);
psi_u.clear();
psi_u.shrink_to_fit();
dpsi_u.clear();
dpsi_u.shrink_to_fit();
d2psi_u.clear();
d2psi_u.shrink_to_fit();

// (2)For each atom, calculate the adjacent atoms in different cells
// and allocate the space for H(R) and S(R).
// If k point is used here, allocate HlocR after atom_arrange.
this->RA.for_2d(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());

// 2. density matrix extrapolation

Expand All @@ -63,12 +144,10 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
{
nsk = PARAM.inp.nspin;
ncol = this->pv.ncol_bands;
if (PARAM.inp.ks_solver == "genelpa"
|| PARAM.inp.ks_solver == "elpa"
|| PARAM.inp.ks_solver == "lapack"
|| PARAM.inp.ks_solver == "pexsi"
|| PARAM.inp.ks_solver == "cusolver"
|| PARAM.inp.ks_solver == "cusolvermp") {
if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "lapack"
|| PARAM.inp.ks_solver == "pexsi" || PARAM.inp.ks_solver == "cusolver"
|| PARAM.inp.ks_solver == "cusolvermp")
{
ncol = this->pv.ncol;
}
}
Expand All @@ -85,12 +164,11 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
}

// init wfc from file
if(istep == 0 && PARAM.inp.init_wfc == "file")
if (istep == 0 && PARAM.inp.init_wfc == "file")
{
if (! ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec))
if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec))
{
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO<TK, TR>::beforesolver",
"read wfc nao failed");
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO<TK, TR>::beforesolver", "read wfc nao failed");
}
}

Expand All @@ -116,10 +194,11 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
orb_,
DM
#ifdef __EXX
, istep
, GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step
, GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr
, GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs
,
istep,
GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step,
GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr,
GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs
#endif
);
}
Expand Down Expand Up @@ -169,41 +248,6 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
{
GlobalC::ucell.cal_ux();
}
ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver");
}

template <typename TK, typename TR>
void ESolver_KS_LCAO<TK, TR>::before_scf(const int istep)
{
ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf");

//! 1) call before_scf() of ESolver_FP
ESolver_FP::before_scf(istep);

if (GlobalC::ucell.ionic_position_updated)
{
this->CE.update_all_dis(GlobalC::ucell);
this->CE.extrapolate_charge(
#ifdef __MPI
&(GlobalC::Pgrid),
#endif
GlobalC::ucell,
this->pelec->charge,
&(this->sf),
GlobalV::ofs_running,
GlobalV::ofs_warning);
}

//----------------------------------------------------------
// about vdw, jiyy add vdwd3 and linpz add vdwd2
//----------------------------------------------------------
auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp, &(GlobalV::ofs_running));
if (vdw_solver != nullptr)
{
this->pelec->f_en.evdw = vdw_solver->get_energy();
}

this->beforesolver(istep);

// Peize Lin add 2016-12-03
#ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf
Expand Down
181 changes: 180 additions & 1 deletion source/module_esolver/lcao_others.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,186 @@ void ESolver_KS_LCAO<TK, TR>::others(const int istep)
return;
}

this->beforesolver(istep);
// 1. prepare HS matrices, prepare grid integral
// (1) Find adjacent atoms for each atom.
double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running,
PARAM.inp.out_level,
orb_.get_rcutmax_Phi(),
GlobalC::ucell.infoNL.get_rcutmax_Beta(),
PARAM.globalv.gamma_only_local);

atom_arrange::search(PARAM.inp.search_pbc,
GlobalV::ofs_running,
GlobalC::GridD,
GlobalC::ucell,
search_radius,
PARAM.inp.test_atom_input);

// (3) Periodic condition search for each grid.
double dr_uniform = 0.001;
std::vector<double> rcuts;
std::vector<std::vector<double>> psi_u;
std::vector<std::vector<double>> dpsi_u;
std::vector<std::vector<double>> d2psi_u;

Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, orb_, psi_u, dpsi_u, d2psi_u);

this->GridT.set_pbc_grid(this->pw_rho->nx,
this->pw_rho->ny,
this->pw_rho->nz,
this->pw_big->bx,
this->pw_big->by,
this->pw_big->bz,
this->pw_big->nbx,
this->pw_big->nby,
this->pw_big->nbz,
this->pw_big->nbxx,
this->pw_big->nbzp_start,
this->pw_big->nbzp,
this->pw_rho->ny,
this->pw_rho->nplane,
this->pw_rho->startz_current,
GlobalC::ucell,
GlobalC::GridD,
dr_uniform,
rcuts,
psi_u,
dpsi_u,
d2psi_u,
PARAM.inp.nstream);
psi_u.clear();
psi_u.shrink_to_fit();
dpsi_u.clear();
dpsi_u.shrink_to_fit();
d2psi_u.clear();
d2psi_u.shrink_to_fit();

// (2)For each atom, calculate the adjacent atoms in different cells
// and allocate the space for H(R) and S(R).
// If k point is used here, allocate HlocR after atom_arrange.
this->RA.for_2d(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());

// 2. density matrix extrapolation

// set the augmented orbitals index.
// after ParaO and GridT,
// this information is used to calculate
// the force.

// init psi
if (this->psi == nullptr)
{
int nsk = 0;
int ncol = 0;
if (PARAM.globalv.gamma_only_local)
{
nsk = PARAM.inp.nspin;
ncol = this->pv.ncol_bands;
if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "lapack"
|| PARAM.inp.ks_solver == "pexsi" || PARAM.inp.ks_solver == "cusolver"
|| PARAM.inp.ks_solver == "cusolvermp")
{
ncol = this->pv.ncol;
}
}
else
{
nsk = this->kv.get_nks();
#ifdef __MPI
ncol = this->pv.ncol_bands;
#else
ncol = PARAM.inp.nbands;
#endif
}
this->psi = new psi::Psi<TK>(nsk, ncol, this->pv.nrow, nullptr);
}

// init wfc from file
if (istep == 0 && PARAM.inp.init_wfc == "file")
{
if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec))
{
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO<TK, TR>::beforesolver", "read wfc nao failed");
}
}

// prepare grid in Gint
LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, orb_, *this->pw_rho, *this->pw_big);

// init Hamiltonian
if (this->p_hamilt != nullptr)
{
delete this->p_hamilt;
this->p_hamilt = nullptr;
}
if (this->p_hamilt == nullptr)
{
elecstate::DensityMatrix<TK, double>* DM = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr,
PARAM.globalv.gamma_only_local ? nullptr : &(this->GK),
&this->pv,
this->pelec->pot,
this->kv,
two_center_bundle_,
orb_,
DM
#ifdef __EXX
,
istep,
GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step,
GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr,
GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs
#endif
);
}

#ifdef __DEEPKS
// for each ionic step, the overlap <psi|alpha> must be rebuilt
// since it depends on ionic positions
if (PARAM.globalv.deepks_setorb)
{
const Parallel_Orbitals* pv = &this->pv;
// build and save <psi(0)|alpha(R)> at beginning
GlobalC::ld.build_psialpha(PARAM.inp.cal_force,
GlobalC::ucell,
orb_,
GlobalC::GridD,
*(two_center_bundle_.overlap_orb_alpha));

if (PARAM.inp.deepks_out_unittest)
{
GlobalC::ld.check_psialpha(PARAM.inp.cal_force, GlobalC::ucell, orb_, GlobalC::GridD);
}
}
#endif
if (PARAM.inp.sc_mag_switch)
{
spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();
sc.init_sc(PARAM.inp.sc_thr,
PARAM.inp.nsc,
PARAM.inp.nsc_min,
PARAM.inp.alpha_trial,
PARAM.inp.sccut,
PARAM.inp.sc_drop_thr,
GlobalC::ucell,
&(this->pv),
PARAM.inp.nspin,
this->kv,
PARAM.inp.ks_solver,
this->p_hamilt,
this->psi,
this->pelec);
}
//=========================================================
// cal_ux should be called before init_scf because
// the direction of ux is used in noncoline_rho
//=========================================================
if (PARAM.inp.nspin == 4)
{
GlobalC::ucell.cal_ux();
}

// pelec should be initialized before these calculations
this->pelec->init_scf(istep, this->sf.strucFac, GlobalC::ucell.symm);
// self consistent calculations for electronic ground state
Expand Down
Loading
Loading