Skip to content

Commit

Permalink
Refactor: remove GlobalV in charge extra
Browse files Browse the repository at this point in the history
  • Loading branch information
@YuLiu98
YuLiu98 committed Aug 13, 2024
1 parent 208cf7e commit 790545b
Show file tree
Hide file tree
Showing 7 changed files with 79 additions and 51 deletions.
63 changes: 35 additions & 28 deletions source/module_elecstate/module_charge/charge_extra.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,21 +20,25 @@ Charge_Extra::~Charge_Extra()
}
}

void Charge_Extra::Init_CE(const int& natom, const double& volume, const int& nrxx)
void Charge_Extra::Init_CE(const int& nspin,
const int& natom,
const double& volume,
const int& nrxx,
const std::string chg_extrap)
{
if(GlobalV::chg_extrap == "none")
if (chg_extrap == "none")
{
pot_order = 0;
}
else if(GlobalV::chg_extrap == "atomic")
else if (chg_extrap == "atomic")
{
pot_order = 1;
}
else if(GlobalV::chg_extrap == "first-order")
else if (chg_extrap == "first-order")
{
pot_order = 2;
}
else if(GlobalV::chg_extrap == "second-order")
else if (chg_extrap == "second-order")
{
pot_order = 3;
}
Expand All @@ -45,11 +49,12 @@ void Charge_Extra::Init_CE(const int& natom, const double& volume, const int& nr

if (pot_order > 1)
{
delta_rho1.resize(GlobalV::NSPIN, std::vector<double>(nrxx, 0.0));
delta_rho2.resize(GlobalV::NSPIN, std::vector<double>(nrxx, 0.0));
delta_rho1.resize(this->nspin, std::vector<double>(nrxx, 0.0));
delta_rho2.resize(this->nspin, std::vector<double>(nrxx, 0.0));
}

this->omega_old = volume;
this->nspin = nspin;

if(pot_order == 3)
{
Expand All @@ -68,7 +73,9 @@ void Charge_Extra::extrapolate_charge(
#endif
UnitCell& ucell,
Charge* chr,
Structure_Factor* sf)
Structure_Factor* sf,
std::ofstream& ofs_running,
std::ofstream& ofs_warning)
{
ModuleBase::TITLE("Charge_Extra","extrapolate_charge");
ModuleBase::timer::tick("Charge_Extra", "extrapolate_charge");
Expand Down Expand Up @@ -96,23 +103,23 @@ void Charge_Extra::extrapolate_charge(
if(rho_extr == 0)
{
sf->setup_structure_factor(&ucell, chr->rhopw);
GlobalV::ofs_running << " charge density from previous step !" << std::endl;
ofs_running << " charge density from previous step !" << std::endl;
return;
}


// if(lsda || noncolin) rho2zeta();

double** rho_atom = new double*[GlobalV::NSPIN];
for (int is = 0; is < GlobalV::NSPIN; is++)
double** rho_atom = new double*[this->nspin];
for (int is = 0; is < this->nspin; is++)
{
rho_atom[is] = new double[chr->rhopw->nrxx];
}
chr->atomic_rho(GlobalV::NSPIN, omega_old, rho_atom, sf->strucFac, ucell);
chr->atomic_rho(this->nspin, omega_old, rho_atom, sf->strucFac, ucell);
#ifdef _OPENMP
#pragma omp parallel for collapse(2) schedule(static, 512)
#endif
for (int is = 0; is < GlobalV::NSPIN; is++)
for (int is = 0; is < this->nspin; is++)
{
for (int ir = 0; ir < chr->rhopw->nrxx; ir++)
{
Expand All @@ -123,14 +130,14 @@ void Charge_Extra::extrapolate_charge(

if(rho_extr == 1)
{
GlobalV::ofs_running << " NEW-OLD atomic charge density approx. for the potential !" << std::endl;
ofs_running << " NEW-OLD atomic charge density approx. for the potential !" << std::endl;

if (pot_order > 1)
{
#ifdef _OPENMP
#pragma omp parallel for collapse(2) schedule(static, 512)
#endif
for (int is = 0; is < GlobalV::NSPIN; is++)
for (int is = 0; is < this->nspin; is++)
{
for (int ir = 0; ir < chr->rhopw->nrxx; ir++)
{
Expand All @@ -142,12 +149,12 @@ void Charge_Extra::extrapolate_charge(
// first order extrapolation
else if(rho_extr ==2)
{
GlobalV::ofs_running << " first order charge density extrapolation !" << std::endl;
ofs_running << " first order charge density extrapolation !" << std::endl;

#ifdef _OPENMP
#pragma omp parallel for collapse(2) schedule(static, 128)
#endif
for (int is = 0; is < GlobalV::NSPIN; is++)
for (int is = 0; is < this->nspin; is++)
{
for (int ir = 0; ir < chr->rhopw->nrxx; ir++)
{
Expand All @@ -160,18 +167,18 @@ void Charge_Extra::extrapolate_charge(
// second order extrapolation
else
{
GlobalV::ofs_running << " second order charge density extrapolation !" << std::endl;
ofs_running << " second order charge density extrapolation !" << std::endl;

find_alpha_and_beta(ucell.nat);
find_alpha_and_beta(ucell.nat, ofs_running, ofs_warning);

const double one_add_alpha = 1 + alpha;
const double beta_alpha = beta - alpha;

std::vector<std::vector<double>> delta_rho3(GlobalV::NSPIN, std::vector<double>(chr->rhopw->nrxx, 0.0));
std::vector<std::vector<double>> delta_rho3(this->nspin, std::vector<double>(chr->rhopw->nrxx, 0.0));
#ifdef _OPENMP
#pragma omp parallel for collapse(2) schedule(static, 64)
#endif
for (int is = 0; is < GlobalV::NSPIN; is++)
for (int is = 0; is < this->nspin; is++)
{
for (int ir = 0; ir < chr->rhopw->nrxx; ir++)
{
Expand All @@ -185,11 +192,11 @@ void Charge_Extra::extrapolate_charge(
}

sf->setup_structure_factor(&ucell, chr->rhopw);
chr->atomic_rho(GlobalV::NSPIN, ucell.omega, rho_atom, sf->strucFac, ucell);
chr->atomic_rho(this->nspin, ucell.omega, rho_atom, sf->strucFac, ucell);
#ifdef _OPENMP
#pragma omp parallel for collapse(2) schedule(static, 512)
#endif
for(int is=0; is<GlobalV::NSPIN; is++)
for (int is = 0; is < this->nspin; is++)
{
for(int ir=0; ir<chr->rhopw->nrxx; ir++)
{
Expand All @@ -200,7 +207,7 @@ void Charge_Extra::extrapolate_charge(

omega_old = ucell.omega;

for(int is=0; is<GlobalV::NSPIN; is++)
for (int is = 0; is < this->nspin; is++)
{
delete[] rho_atom[is];
}
Expand All @@ -209,7 +216,7 @@ void Charge_Extra::extrapolate_charge(
return;
}

void Charge_Extra::find_alpha_and_beta(const int& natom)
void Charge_Extra::find_alpha_and_beta(const int& natom, std::ofstream& ofs_running, std::ofstream& ofs_warning)
{
if(istep < 3) return;

Expand Down Expand Up @@ -244,7 +251,7 @@ void Charge_Extra::find_alpha_and_beta(const int& natom)
alpha = 0.0;
beta = 0.0;

ModuleBase::GlobalFunc::OUT(GlobalV::ofs_warning,"in find_alpha_and beta() det = ", det);
ModuleBase::GlobalFunc::OUT(ofs_warning, "in find_alpha_and beta() det = ", det);
}

if(det > 1e-20)
Expand All @@ -263,8 +270,8 @@ void Charge_Extra::find_alpha_and_beta(const int& natom)
}
}

GlobalV::ofs_running << " alpha = " << alpha << std::endl;
GlobalV::ofs_running << " beta = " << beta << std::endl;
ofs_running << " alpha = " << alpha << std::endl;
ofs_running << " beta = " << beta << std::endl;

return;
}
Expand Down
19 changes: 16 additions & 3 deletions source/module_elecstate/module_charge/charge_extra.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,11 +45,17 @@ class Charge_Extra
* But after ucell and Esolver are fully decoupled
* Init_CE will be removed and everything put back in the constructor
*
* @param nspin the number of spins
* @param natom the number of atoms
* @param volume the volume of the cell
* @param nrxx the number of grids
* @param chg_extrap the charge extrapolation method
*/
void Init_CE(const int& natom, const double& volume, const int& nrxx);
void Init_CE(const int& nspin,
const int& natom,
const double& volume,
const int& nrxx,
const std::string chg_extrap);

/**
* @brief charge extrapolation method
Expand All @@ -58,14 +64,18 @@ class Charge_Extra
* @param ucell the cell information
* @param chr the charge density
* @param sf the structure factor
* @param ofs_running the output stream
* @param ofs_warning the output stream
*/
void extrapolate_charge(
#ifdef __MPI
Parallel_Grid* Pgrid,
#endif
UnitCell& ucell,
Charge* chr,
Structure_Factor* sf);
Structure_Factor* sf,
std::ofstream& ofs_running,
std::ofstream& ofs_warning);

/**
* @brief update displacements
Expand All @@ -82,6 +92,7 @@ class Charge_Extra
int pot_order; ///< the specified charge extrapolation method
int rho_extr; ///< the actually used method
double omega_old; ///< the old volume of the last step
int nspin; ///< the number of spins

ModuleBase::Vector3<double>* dis_old1 = nullptr; ///< dis_old2 = pos_old1 - pos_old2
ModuleBase::Vector3<double>* dis_old2 = nullptr; ///< dis_old1 = pos_now - pos_old1
Expand All @@ -97,8 +108,10 @@ class Charge_Extra
* @brief determine alpha and beta
*
* @param natom the number of atoms
* @param ofs_running the output stream
* @param ofs_warning the output stream
*/
void find_alpha_and_beta(const int& natom);
void find_alpha_and_beta(const int& natom, std::ofstream& ofs_running, std::ofstream& ofs_warning);
};

#endif
34 changes: 18 additions & 16 deletions source/module_elecstate/test/charge_extra_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,29 +136,31 @@ TEST_F(ChargeExtraTest, InitCEWarningQuit)
{
GlobalV::chg_extrap = "wwww";
testing::internal::CaptureStdout();
EXPECT_EXIT(CE.Init_CE(ucell->nat, ucell->omega, charge.rhopw->nrxx), ::testing::ExitedWithCode(0), "");
EXPECT_EXIT(CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap),
::testing::ExitedWithCode(0),
"");
std::string output = testing::internal::GetCapturedStdout();
EXPECT_THAT(output, testing::HasSubstr("charge extrapolation method is not available"));
}

TEST_F(ChargeExtraTest, InitCECase1)
{
GlobalV::chg_extrap = "none";
CE.Init_CE(ucell->nat, ucell->omega, charge.rhopw->nrxx);
CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
EXPECT_EQ(CE.pot_order, 0);
}

TEST_F(ChargeExtraTest, InitCECase2)
{
GlobalV::chg_extrap = "atomic";
CE.Init_CE(ucell->nat, ucell->omega, charge.rhopw->nrxx);
CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
EXPECT_EQ(CE.pot_order, 1);
}

TEST_F(ChargeExtraTest, InitCECase3)
{
GlobalV::chg_extrap = "first-order";
CE.Init_CE(ucell->nat, ucell->omega, charge.rhopw->nrxx);
CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
EXPECT_EQ(CE.pot_order, 2);
EXPECT_NE(CE.delta_rho1.size(), 0);
EXPECT_NE(CE.delta_rho2.size(), 0);
Expand All @@ -167,7 +169,7 @@ TEST_F(ChargeExtraTest, InitCECase3)
TEST_F(ChargeExtraTest, InitCECase4)
{
GlobalV::chg_extrap = "second-order";
CE.Init_CE(ucell->nat, ucell->omega, charge.rhopw->nrxx);
CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
EXPECT_EQ(CE.pot_order, 3);
EXPECT_DOUBLE_EQ(CE.alpha, 1.0);
EXPECT_DOUBLE_EQ(CE.beta, 0.0);
Expand All @@ -181,12 +183,12 @@ TEST_F(ChargeExtraTest, InitCECase4)
TEST_F(ChargeExtraTest, ExtrapolateChargeCase1)
{
GlobalV::chg_extrap = "second-order";
CE.Init_CE(ucell->nat, ucell->omega, charge.rhopw->nrxx);
CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
CE.istep = 0;
CE.pot_order = 3;

GlobalV::ofs_running.open("log");
CE.extrapolate_charge(*ucell.get(), &charge, &sf);
CE.extrapolate_charge(*ucell.get(), &charge, &sf, GlobalV::ofs_running, GlobalV::ofs_warning);
GlobalV::ofs_running.close();

// Check the results
Expand All @@ -203,12 +205,12 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase1)
TEST_F(ChargeExtraTest, ExtrapolateChargeCase2)
{
GlobalV::chg_extrap = "second-order";
CE.Init_CE(ucell->nat, ucell->omega, charge.rhopw->nrxx);
CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
CE.istep = 1;
CE.pot_order = 3;

GlobalV::ofs_running.open("log");
CE.extrapolate_charge(*ucell.get(), &charge, &sf);
CE.extrapolate_charge(*ucell.get(), &charge, &sf, GlobalV::ofs_running, GlobalV::ofs_warning);
GlobalV::ofs_running.close();

// Check the results
Expand All @@ -225,12 +227,12 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase2)
TEST_F(ChargeExtraTest, ExtrapolateChargeCase3)
{
GlobalV::chg_extrap = "second-order";
CE.Init_CE(ucell->nat, ucell->omega, charge.rhopw->nrxx);
CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
CE.istep = 2;
CE.pot_order = 3;

GlobalV::ofs_running.open("log");
CE.extrapolate_charge(*ucell.get(), &charge, &sf);
CE.extrapolate_charge(*ucell.get(), &charge, &sf, GlobalV::ofs_running, GlobalV::ofs_warning);
GlobalV::ofs_running.close();

// Check the results
Expand All @@ -247,11 +249,11 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase3)
TEST_F(ChargeExtraTest, ExtrapolateChargeCase4)
{
GlobalV::chg_extrap = "second-order";
CE.Init_CE(ucell->nat, ucell->omega, charge.rhopw->nrxx);
CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
CE.istep = 3;

GlobalV::ofs_running.open("log");
CE.extrapolate_charge(*ucell.get(), &charge, &sf);
CE.extrapolate_charge(*ucell.get(), &charge, &sf, GlobalV::ofs_running, GlobalV::ofs_warning);
GlobalV::ofs_running.close();

// Check the results
Expand All @@ -269,7 +271,7 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase4)
TEST_F(ChargeExtraTest, UpdateAllDis)
{
GlobalV::chg_extrap = "second-order";
CE.Init_CE(ucell->nat, ucell->omega, charge.rhopw->nrxx);
CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
CE.istep = 3;
for (int i = 0; i < ucell->nat; ++i)
{
Expand All @@ -291,7 +293,7 @@ TEST_F(ChargeExtraTest, UpdateAllDis)
TEST_F(ChargeExtraTest, FindAlphaAndBeta)
{
GlobalV::chg_extrap = "second-order";
CE.Init_CE(ucell->nat, ucell->omega, charge.rhopw->nrxx);
CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
CE.istep = 3;
for (int i = 0; i < ucell->nat; ++i)
{
Expand All @@ -302,7 +304,7 @@ TEST_F(ChargeExtraTest, FindAlphaAndBeta)
}
}

CE.find_alpha_and_beta(ucell->nat);
CE.find_alpha_and_beta(ucell->nat, GlobalV::ofs_running, GlobalV::ofs_warning);

EXPECT_DOUBLE_EQ(CE.alpha, 1.0);
EXPECT_DOUBLE_EQ(CE.beta, 0.0);
Expand Down
2 changes: 1 addition & 1 deletion source/module_esolver/esolver_fp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,7 +105,7 @@ void ESolver_FP::before_all_runners(const Input_para& inp, UnitCell& cell)
this->print_rhofft(inp, GlobalV::ofs_running);

//! 5) initialize the charge extrapolation method if necessary
this->CE.Init_CE(GlobalC::ucell.nat, GlobalC::ucell.omega, this->pw_rhod->nrxx);
this->CE.Init_CE(GlobalV::NSPIN, GlobalC::ucell.nat, GlobalC::ucell.omega, this->pw_rhod->nrxx, inp.chg_extrap);

return;
}
Expand Down
Loading

0 comments on commit 790545b

Please sign in to comment.