Skip to content

Commit 35fe565

Browse files
YuLiu98YuLiu98
authored
Refactor: update logic of init_chg (#5801)
1 parent d061719 commit 35fe565

File tree

1 file changed

+74
-45
lines changed

1 file changed

+74
-45
lines changed

‎source/module_elecstate/module_charge/charge_init.cpp‎

Lines changed: 74 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
3737
this->pgrid = &pgrid;
3838

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

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

115-
std::stringstream binary;
116-
binary << PARAM.globalv.global_readin_dir << PARAM.inp.suffix + "-TAU-DENSITY.restart";
117-
if (ModuleIO::read_rhog(binary.str(), rhopw, kin_g.data()))
120+
if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
121+
{
122+
// If the charge density is not read in, then the kinetic energy density is not read in either
123+
if (!read_error)
118124
{
119-
GlobalV::ofs_running << " Read in the kinetic energy density: " << binary.str() << std::endl;
120-
for (int is = 0; is < PARAM.inp.nspin; ++is)
125+
GlobalV::ofs_running << " try to read kinetic energy density from file" << std::endl;
126+
// try to read charge from binary file first, which is the same as QE
127+
std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->ngmc, {0.0, 0.0});
128+
std::vector<std::complex<double>*> kin_g;
129+
for (int is = 0; is < PARAM.inp.nspin; is++)
121130
{
122-
rhopw->recip2real(kin_g[is], this->kin_r[is]);
131+
kin_g.push_back(kin_g_space.data() + is * this->ngmc);
123132
}
124-
}
125-
else
126-
{
127-
for (int is = 0; is < PARAM.inp.nspin; is++)
133+
134+
std::stringstream binary;
135+
binary << PARAM.globalv.global_readin_dir << PARAM.inp.suffix + "-TAU-DENSITY.restart";
136+
if (ModuleIO::read_rhog(binary.str(), rhopw, kin_g.data()))
128137
{
129-
std::stringstream ssc;
130-
ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
131-
// mohan update 2012-02-10, sunliang update 2023-03-09
132-
if (ModuleIO::read_vdata_palgrid(
133-
pgrid,
134-
(PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_STOGROUP : GlobalV::MY_RANK),
135-
GlobalV::ofs_running,
136-
ssc.str(),
137-
this->kin_r[is],
138-
ucell.nat))
138+
GlobalV::ofs_running << " Read in the kinetic energy density: " << binary.str() << std::endl;
139+
for (int is = 0; is < PARAM.inp.nspin; ++is)
139140
{
140-
GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
141+
rhopw->recip2real(kin_g[is], this->kin_r[is]);
141142
}
142143
}
144+
else
145+
{
146+
for (int is = 0; is < PARAM.inp.nspin; is++)
147+
{
148+
std::stringstream ssc;
149+
ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
150+
// mohan update 2012-02-10, sunliang update 2023-03-09
151+
if (ModuleIO::read_vdata_palgrid(
152+
pgrid,
153+
(PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_STOGROUP : GlobalV::MY_RANK),
154+
GlobalV::ofs_running,
155+
ssc.str(),
156+
this->kin_r[is],
157+
ucell.nat))
158+
{
159+
GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
160+
}
161+
else
162+
{
163+
read_kin_error = true;
164+
std::cout << " WARNING: \"init_chg\" is enabled but ABACUS failed to read kinetic energy "
165+
"density from file.\n"
166+
" Please check if there is SPINX_TAU.cube (X=1,...) or "
167+
"{suffix}-TAU-DENSITY.restart in the directory.\n"
168+
<< std::endl;
169+
break;
170+
}
171+
}
172+
}
173+
}
174+
else
175+
{
176+
read_kin_error = true;
143177
}
144178
}
145179
}
146-
if (read_error)
180+
181+
if (PARAM.inp.init_chg == "atomic" || read_error) // mohan add 2007-10-17
147182
{
148-
const std::string warn_msg = " WARNING: \"init_chg\" is enabled but ABACUS failed to read charge density from file.\n"
149-
" Please check if there is SPINX_CHG.cube (X=1,...) or {suffix}-CHARGE-DENSITY.restart in the directory.\n";
150-
std::cout << std::endl << warn_msg;
151-
if (PARAM.inp.init_chg == "auto")
183+
if (read_error)
152184
{
153-
std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl << std::endl;
154-
}
155-
else if (PARAM.inp.init_chg == "file")
156-
{
157-
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.");
185+
std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl;
158186
}
187+
this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);
159188
}
160189

161-
if (PARAM.inp.init_chg == "atomic" ||
162-
(PARAM.inp.init_chg == "auto" && read_error)) // mohan add 2007-10-17
190+
// wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
191+
if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
163192
{
164-
this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);
165-
166-
// liuyu 2023-06-29 : move here from atomic_rho(), which will be called several times in charge extrapolation
167-
// wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
168-
if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
193+
if (PARAM.inp.init_chg == "atomic" || read_kin_error)
169194
{
195+
if (read_kin_error)
196+
{
197+
std::cout << " Charge::init_rho: init kinetic energy density from rho." << std::endl;
198+
}
170199
const double pi = 3.141592653589790;
171200
const double fact = (3.0 / 5.0) * pow(3.0 * pi * pi, 2.0 / 3.0);
172201
for (int is = 0; is < PARAM.inp.nspin; ++is)

0 commit comments

Comments
 (0)