Skip to content

Refactor: update logic of init_chg #5801

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
Jan 4, 2025
Merged

Refactor: update logic of init_chg #5801

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
119 changes: 74 additions & 45 deletions source/module_elecstate/module_charge/charge_init.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
this->pgrid = &pgrid;

bool read_error = false;
bool read_kin_error = false;
if (PARAM.inp.init_chg == "file" || PARAM.inp.init_chg == "auto")
{
GlobalV::ofs_running << " try to read charge from file" << std::endl;
Expand Down Expand Up @@ -101,72 +102,100 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
}
}

if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
if (read_error)
{
GlobalV::ofs_running << " try to read kinetic energy density from file" << std::endl;
// try to read charge from binary file first, which is the same as QE
std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->ngmc, {0.0, 0.0});
std::vector<std::complex<double>*> kin_g;
for (int is = 0; is < PARAM.inp.nspin; is++)
const std::string warn_msg
= " WARNING: \"init_chg\" is enabled but ABACUS failed to read charge density from file.\n"
" Please check if there is SPINX_CHG.cube (X=1,...) or {suffix}-CHARGE-DENSITY.restart in the "
"directory.\n";
std::cout << std::endl << warn_msg;
if (PARAM.inp.init_chg == "file")
{
kin_g.push_back(kin_g_space.data() + is * this->ngmc);
ModuleBase::WARNING_QUIT("Charge::init_rho",
"Failed to read in charge density from file.\nIf you want to use atomic "
"charge initialization, \nplease set init_chg to atomic in INPUT.");
}
}

std::stringstream binary;
binary << PARAM.globalv.global_readin_dir << PARAM.inp.suffix + "-TAU-DENSITY.restart";
if (ModuleIO::read_rhog(binary.str(), rhopw, kin_g.data()))
if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
{
// If the charge density is not read in, then the kinetic energy density is not read in either
if (!read_error)
{
GlobalV::ofs_running << " Read in the kinetic energy density: " << binary.str() << std::endl;
for (int is = 0; is < PARAM.inp.nspin; ++is)
GlobalV::ofs_running << " try to read kinetic energy density from file" << std::endl;
// try to read charge from binary file first, which is the same as QE
std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->ngmc, {0.0, 0.0});
std::vector<std::complex<double>*> kin_g;
for (int is = 0; is < PARAM.inp.nspin; is++)
{
rhopw->recip2real(kin_g[is], this->kin_r[is]);
kin_g.push_back(kin_g_space.data() + is * this->ngmc);
}
}
else
{
for (int is = 0; is < PARAM.inp.nspin; is++)

std::stringstream binary;
binary << PARAM.globalv.global_readin_dir << PARAM.inp.suffix + "-TAU-DENSITY.restart";
if (ModuleIO::read_rhog(binary.str(), rhopw, kin_g.data()))
{
std::stringstream ssc;
ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
// mohan update 2012-02-10, sunliang update 2023-03-09
if (ModuleIO::read_vdata_palgrid(
pgrid,
(PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_STOGROUP : GlobalV::MY_RANK),
GlobalV::ofs_running,
ssc.str(),
this->kin_r[is],
ucell.nat))
GlobalV::ofs_running << " Read in the kinetic energy density: " << binary.str() << std::endl;
for (int is = 0; is < PARAM.inp.nspin; ++is)
{
GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
rhopw->recip2real(kin_g[is], this->kin_r[is]);
}
}
else
{
for (int is = 0; is < PARAM.inp.nspin; is++)
{
std::stringstream ssc;
ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
// mohan update 2012-02-10, sunliang update 2023-03-09
if (ModuleIO::read_vdata_palgrid(
pgrid,
(PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_STOGROUP : GlobalV::MY_RANK),
GlobalV::ofs_running,
ssc.str(),
this->kin_r[is],
ucell.nat))
{
GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
}
else
{
read_kin_error = true;
std::cout << " WARNING: \"init_chg\" is enabled but ABACUS failed to read kinetic energy "
"density from file.\n"
" Please check if there is SPINX_TAU.cube (X=1,...) or "
"{suffix}-TAU-DENSITY.restart in the directory.\n"
<< std::endl;
break;
}
}
}
}
else
{
read_kin_error = true;
}
}
}
if (read_error)

if (PARAM.inp.init_chg == "atomic" || read_error) // mohan add 2007-10-17
{
const std::string warn_msg = " WARNING: \"init_chg\" is enabled but ABACUS failed to read charge density from file.\n"
" Please check if there is SPINX_CHG.cube (X=1,...) or {suffix}-CHARGE-DENSITY.restart in the directory.\n";
std::cout << std::endl << warn_msg;
if (PARAM.inp.init_chg == "auto")
if (read_error)
{
std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl << std::endl;
}
else if (PARAM.inp.init_chg == "file")
{
ModuleBase::WARNING_QUIT("Charge::init_rho", "Failed to read in charge density from file.\nIf you want to use atomic charge initialization, \nplease set init_chg to atomic in INPUT.");
std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl;
}
this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);
}

if (PARAM.inp.init_chg == "atomic" ||
(PARAM.inp.init_chg == "auto" && read_error)) // mohan add 2007-10-17
// wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
{
this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);

// liuyu 2023-06-29 : move here from atomic_rho(), which will be called several times in charge extrapolation
// wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
if (PARAM.inp.init_chg == "atomic" || read_kin_error)
{
if (read_kin_error)
{
std::cout << " Charge::init_rho: init kinetic energy density from rho." << std::endl;
}
const double pi = 3.141592653589790;
const double fact = (3.0 / 5.0) * pow(3.0 * pi * pi, 2.0 / 3.0);
for (int is = 0; is < PARAM.inp.nspin; ++is)
Expand Down
Loading