Skip to content

Refactor: rename Kerker_screen_recip() and add more UnitTests to protect Charge_Mixing #3609

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
WHUweiqingzhou merged 7 commits into from
Feb 7, 2024
Merged

Refactor: rename Kerker_screen_recip() and add more UnitTests to protect Charge_Mixing #3609

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
3055b34
delete old kerker_screen_recip()
WHUweiqingzhou Feb 6, 2024
8166de5
rename mix_rho_recip_new() as mix_rho_recip()
WHUweiqingzhou Feb 6, 2024
82f0cd3
rename Kerker_screen_recip_new() as Kerker_screen_recip()
WHUweiqingzhou Feb 6, 2024
71e63b4
use local vars in Kerker_screen_recip
WHUweiqingzhou Feb 6, 2024
3f1682e
use local vars in Kerker_screen_real
WHUweiqingzhou Feb 6, 2024
2005a46
use local vars in inner_product_real()
WHUweiqingzhou Feb 6, 2024
45ba106
polish UnitTests to protect Kerker_screen_recip() and Kerker_screen_r…
WHUweiqingzhou Feb 7, 2024
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
59 changes: 19 additions & 40 deletions source/module_elecstate/module_charge/charge_mixing.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -264,7 +264,7 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec)
return drho;
}

void Charge_Mixing::mix_rho_recip_new(Charge* chr)
void Charge_Mixing::mix_rho_recip(Charge* chr)
{
std::complex<double>* rhog_in = nullptr;
std::complex<double>* rhog_out = nullptr;
Expand All @@ -291,7 +291,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr)
{
rhog_in = rhogs_in;
rhog_out = rhogs_out;
auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, this, std::placeholders::_1);
auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip, this, std::placeholders::_1);
this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, true);
this->mixing->cal_coef(this->rho_mdata, inner_product);
this->mixing->mix_data(this->rho_mdata, rhog_out);
Expand Down Expand Up @@ -322,7 +322,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr)
rhog_in = rhog_mag_save;
rhog_out = rhog_mag;
//
auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, this, std::placeholders::_1);
auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip, this, std::placeholders::_1);
auto twobeta_mix
= [this, npw](std::complex<double>* out, const std::complex<double>* in, const std::complex<double>* sres) {
#ifdef _OPENMP
Expand Down Expand Up @@ -373,7 +373,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr)
rhog_in = rhogs_in;
rhog_out = rhogs_out;
const int npw = this->rhopw->npw;
auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, this, std::placeholders::_1); // use old one
auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip, this, std::placeholders::_1); // use old one
auto twobeta_mix
= [this, npw](std::complex<double>* out, const std::complex<double>* in, const std::complex<double>* sres) {
#ifdef _OPENMP
Expand Down Expand Up @@ -437,7 +437,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr)
//
rhog_in = rhog_magabs_save;
rhog_out = rhog_magabs;
auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, this, std::placeholders::_1); // use old one
auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip, this, std::placeholders::_1); // use old one
auto twobeta_mix
= [this, npw](std::complex<double>* out, const std::complex<double>* in, const std::complex<double>* sres) {
#ifdef _OPENMP
Expand Down Expand Up @@ -1027,7 +1027,7 @@ void Charge_Mixing::mix_rho(Charge* chr)
// --------------------Mixing Body--------------------
if (GlobalV::SCF_THR_TYPE == 1)
{
mix_rho_recip_new(chr);
mix_rho_recip(chr);
}
else if (GlobalV::SCF_THR_TYPE == 2)
{
Expand Down Expand Up @@ -1090,43 +1090,22 @@ void Charge_Mixing::mix_rho(Charge* chr)
}

void Charge_Mixing::Kerker_screen_recip(std::complex<double>* drhog)
{
if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1)
return;
const double fac = this->mixing_gg0;
const double gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2);
#ifdef _OPENMP
#pragma omp parallel for collapse(2) schedule(static, 512)
#endif
for (int is = 0; is < GlobalV::NSPIN; ++is)
{
for (int ig = 0; ig < this->rhopw->npw; ++ig)
{
double gg = this->rhopw->gg[ig];
double filter_g = std::max(gg / (gg + gg0), GlobalV::MIXING_GG0_MIN / this->mixing_beta);
drhog[is * this->rhopw->npw + ig] *= filter_g;
}
}
return;
}

void Charge_Mixing::Kerker_screen_recip_new(std::complex<double>* drhog)
{
if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1)
return;
double fac, gg0, amin;

// consider a resize for mixing_angle
int resize_tmp = 1;
if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0) resize_tmp = 2;
if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) resize_tmp = 2;

// implement Kerker for density and magnetization separately
for (int is = 0; is < GlobalV::NSPIN / resize_tmp; ++is)
{
// new mixing method only support nspin=2 not nspin=4
if (is >= 1)
{
if (GlobalV::MIXING_GG0_MAG <= 0.0001 || GlobalV::MIXING_BETA_MAG <= 0.1)
if (this->mixing_gg0_mag <= 0.0001 || this->mixing_beta_mag <= 0.1)
{
#ifdef __DEBUG
assert(is == 1); // make sure break works
Expand All @@ -1138,8 +1117,8 @@ void Charge_Mixing::Kerker_screen_recip_new(std::complex<double>* drhog)
//}
break;
}
fac = GlobalV::MIXING_GG0_MAG;
amin = GlobalV::MIXING_BETA_MAG;
fac = this->mixing_gg0_mag;
amin = this->mixing_beta_mag;
}
else
{
Expand All @@ -1154,7 +1133,7 @@ void Charge_Mixing::Kerker_screen_recip_new(std::complex<double>* drhog)
for (int ig = 0; ig < this->rhopw->npw; ++ig)
{
double gg = this->rhopw->gg[ig];
double filter_g = std::max(gg / (gg + gg0), GlobalV::MIXING_GG0_MIN / amin);
double filter_g = std::max(gg / (gg + gg0), this->mixing_gg0_min / amin);
drhog[is * this->rhopw->npw + ig] *= filter_g;
}
}
Expand All @@ -1167,7 +1146,7 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
return;
// consider a resize for mixing_angle
int resize_tmp = 1;
if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0) resize_tmp = 2;
if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) resize_tmp = 2;
//
std::vector<std::complex<double>> drhog(this->rhopw->npw * GlobalV::NSPIN / resize_tmp);
std::vector<double> drhor_filter(this->rhopw->nrxx * GlobalV::NSPIN / resize_tmp);
Expand All @@ -1184,21 +1163,21 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)

if (is >= 1)
{
if (GlobalV::MIXING_GG0_MAG <= 0.0001 || GlobalV::MIXING_BETA_MAG <= 0.1)
if (this->mixing_gg0_mag <= 0.0001 || this->mixing_beta_mag <= 0.1)
{
#ifdef __DEBUG
assert(is == 1); // make sure break works
#endif
double is_mag = GlobalV::NSPIN - 1;
if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0) is_mag = 1;
if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) is_mag = 1;
for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++)
{
drhog[is * this->rhopw->npw + ig] = 0;
}
break;
}
fac = GlobalV::MIXING_GG0_MAG;
amin = GlobalV::MIXING_BETA_MAG;
fac = this->mixing_gg0_mag;
amin = this->mixing_beta_mag;
}
else
{
Expand All @@ -1219,7 +1198,7 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
// drhog[is * this->rhopw->npw + ig] *= 0;
// continue;
//}
double filter_g = std::max(gg / (gg + gg0), GlobalV::MIXING_GG0_MIN / amin);
double filter_g = std::max(gg / (gg + gg0), this->mixing_gg0_min / amin);
drhog[is * this->rhopw->npw + ig] *= (1 - filter_g);
}
}
Expand Down Expand Up @@ -1387,7 +1366,7 @@ double Charge_Mixing::inner_product_recip_simple(std::complex<double>* rho1, std
double rnorm = 0.0;
// consider a resize for mixing_angle
int resize_tmp = 1;
if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0) resize_tmp = 2;
if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) resize_tmp = 2;
#ifdef _OPENMP
#pragma omp parallel for reduction(+ : rnorm)
#endif
Expand Down Expand Up @@ -1575,7 +1554,7 @@ double Charge_Mixing::inner_product_real(double* rho1, double* rho2)
double rnorm = 0.0;
// consider a resize for mixing_angle
int resize_tmp = 1;
if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0) resize_tmp = 2;
if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) resize_tmp = 2;

#ifdef _OPENMP
#pragma omp parallel for reduction(+ : rnorm)
Expand Down
5 changes: 2 additions & 3 deletions source/module_elecstate/module_charge/charge_mixing.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,7 @@ class Charge_Mixing
* @brief charge mixing for reciprocal space
* @param chr pointer of Charge object
*/
void mix_rho_recip_new(Charge* chr);
void mix_rho_recip(Charge* chr);

/**
* @brief charge mixing for real space
Expand All @@ -145,7 +145,6 @@ class Charge_Mixing
* @param rhog charge density in reciprocal space
*/
void Kerker_screen_recip(std::complex<double>* rhog);
void Kerker_screen_recip_new(std::complex<double>* rhog);

/**
* @brief Kerker screen method for real space
Expand All @@ -157,7 +156,7 @@ class Charge_Mixing
* @brief Inner product of two complex vectors
* @brief inner_product_recip_rho is used for charge, like get_drho()
* @brief inner_product_recip_hartree and inner_product_recip_simple are used for charge mixing
* @brief inner_product_recip_simple is used for test
* @brief inner_product_recip_simple is only used for test
* @brief Actually, I am not sure if the definition of inner product for NSPIN=4 is correct, need to be checked.
*/
double inner_product_recip_rho(std::complex<double>* rho1, std::complex<double>* rho2);
Expand Down
Loading