diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 8d95fc8b0c..9e2df119be 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -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* rhog_in = nullptr; std::complex* rhog_out = nullptr; @@ -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); @@ -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* out, const std::complex* in, const std::complex* sres) { #ifdef _OPENMP @@ -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* out, const std::complex* in, const std::complex* sres) { #ifdef _OPENMP @@ -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* out, const std::complex* in, const std::complex* sres) { #ifdef _OPENMP @@ -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) { @@ -1090,27 +1090,6 @@ void Charge_Mixing::mix_rho(Charge* chr) } void Charge_Mixing::Kerker_screen_recip(std::complex* 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* drhog) { if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1) return; @@ -1118,7 +1097,7 @@ void Charge_Mixing::Kerker_screen_recip_new(std::complex* drhog) // 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) @@ -1126,7 +1105,7 @@ void Charge_Mixing::Kerker_screen_recip_new(std::complex* drhog) // 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 @@ -1138,8 +1117,8 @@ void Charge_Mixing::Kerker_screen_recip_new(std::complex* drhog) //} break; } - fac = GlobalV::MIXING_GG0_MAG; - amin = GlobalV::MIXING_BETA_MAG; + fac = this->mixing_gg0_mag; + amin = this->mixing_beta_mag; } else { @@ -1154,7 +1133,7 @@ void Charge_Mixing::Kerker_screen_recip_new(std::complex* 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; } } @@ -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> drhog(this->rhopw->npw * GlobalV::NSPIN / resize_tmp); std::vector drhor_filter(this->rhopw->nrxx * GlobalV::NSPIN / resize_tmp); @@ -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 { @@ -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); } } @@ -1387,7 +1366,7 @@ double Charge_Mixing::inner_product_recip_simple(std::complex* 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 @@ -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) diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 1bf3282359..bc6f74e4d1 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -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 @@ -145,7 +145,6 @@ class Charge_Mixing * @param rhog charge density in reciprocal space */ void Kerker_screen_recip(std::complex* rhog); - void Kerker_screen_recip_new(std::complex* rhog); /** * @brief Kerker screen method for real space @@ -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* rho1, std::complex* rho2); diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index 2a62d3c161..113a77c9b5 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -536,77 +536,168 @@ TEST_F(ChargeMixingTest, KerkerScreenRecipTest) { Charge_Mixing CMtest; CMtest.set_rhopw(&pw_basis, &pw_basis); - GlobalV::NSPIN = 1; GlobalC::ucell.tpiba = 1.0; - + // nspin = 1 + GlobalV::NSPIN = 1; + std::complex* drhog = new std::complex[GlobalV::NSPIN*pw_basis.npw]; + std::complex* drhog_old = new std::complex[GlobalV::NSPIN*pw_basis.npw]; + for (int i = 0; i < GlobalV::NSPIN*pw_basis.npw; ++i) + { + drhog_old[i] = drhog[i] = std::complex(1.0, 1.0); + } + // no kerker CMtest.mixing_gg0 = 0.0; - std::complex* drhog = new std::complex[pw_basis.npw]; - std::complex* drhog_old = new std::complex[pw_basis.npw]; - double* drhor = new double[pw_basis.nrxx]; - double* drhor_ref = new double[pw_basis.nrxx]; + CMtest.Kerker_screen_recip(drhog); + for (int i = 0; i < GlobalV::NSPIN*pw_basis.npw; ++i) + { + EXPECT_EQ(drhog[i], drhog_old[i]); + } + // kerker + CMtest.mixing_gg0 = 1.0; + CMtest.Kerker_screen_recip(drhog); + double gg0 = std::pow(0.529177, 2); for (int i = 0; i < pw_basis.npw; ++i) + { + double gg = this->pw_basis.gg[i]; + double ref = std::max(gg / (gg + gg0), 0.1 / CMtest.mixing_beta); + EXPECT_NEAR(drhog[i].real(), ref, 1e-10); + EXPECT_NEAR(drhog[i].imag(), ref, 1e-10); + } + delete[] drhog; + delete[] drhog_old; + + // nspin = 2 + GlobalV::NSPIN = 2; + CMtest.mixing_beta = 0.4; + CMtest.mixing_beta_mag = 1.6; + drhog = new std::complex[GlobalV::NSPIN*pw_basis.npw]; + drhog_old = new std::complex[GlobalV::NSPIN*pw_basis.npw]; + for (int i = 0; i < GlobalV::NSPIN*pw_basis.npw; ++i) { drhog_old[i] = drhog[i] = std::complex(1.0, 1.0); } + // mixing_gg0 = 0.0 + CMtest.mixing_gg0 = 0.0; CMtest.Kerker_screen_recip(drhog); - for (int i = 0; i < pw_basis.npw; ++i) + for (int i = 0; i < GlobalV::NSPIN*pw_basis.npw; ++i) { EXPECT_EQ(drhog[i], drhog_old[i]); } - - // RECIPROCAL + // mixing_gg0 = 1.0, mixing_gg0_mag = 0.0 CMtest.mixing_gg0 = 1.0; CMtest.Kerker_screen_recip(drhog); - const double gg0 = std::pow(0.529177, 2); + gg0 = std::pow(0.529177, 2); for (int i = 0; i < pw_basis.npw; ++i) { - std::complex ration = drhog[i] / drhog_old[i]; double gg = this->pw_basis.gg[i]; - double ration_ref = std::max(gg / (gg + gg0), 0.1 / CMtest.mixing_beta); - EXPECT_NEAR(ration.real(), ration_ref, 1e-10); - EXPECT_NEAR(ration.imag(), 0, 1e-10); + double ref = std::max(gg / (gg + gg0), 0.1 / CMtest.mixing_beta); + // rho + EXPECT_NEAR(drhog[i].real(), ref, 1e-10); + EXPECT_NEAR(drhog[i].imag(), ref, 1e-10); + // mag + EXPECT_NEAR(drhog[i+pw_basis.npw].real(), 1.0, 1e-10); + EXPECT_NEAR(drhog[i+pw_basis.npw].imag(), 1.0, 1e-10); } + delete[] drhog; + delete[] drhog_old; - // REAL - pw_basis.recip2real(drhog, drhor_ref); - pw_basis.recip2real(drhog_old, drhor); - + // nspin = 4 + GlobalV::NSPIN = 4; + drhog = new std::complex[GlobalV::NSPIN*pw_basis.npw]; + drhog_old = new std::complex[GlobalV::NSPIN*pw_basis.npw]; + for (int i = 0; i < GlobalV::NSPIN*pw_basis.npw; ++i) + { + drhog_old[i] = drhog[i] = std::complex(1.0, 1.0); + } + // mixing_gg0 = 0.0 CMtest.mixing_gg0 = 0.0; - // nothing happens - CMtest.Kerker_screen_real(drhor); - + CMtest.Kerker_screen_recip(drhog); + for (int i = 0; i < GlobalV::NSPIN*pw_basis.npw; ++i) + { + EXPECT_EQ(drhog[i], drhog_old[i]); + } + // mixing_gg0 = 1.0, mixing_gg0_mag = 0.0 CMtest.mixing_gg0 = 1.0; - CMtest.Kerker_screen_real(drhor); - for (int i = 0; i < pw_basis.nrxx; ++i) + CMtest.Kerker_screen_recip(drhog); + gg0 = std::pow(0.529177, 2); + for (int i = 0; i < pw_basis.npw; ++i) { - EXPECT_NEAR(drhor[i], drhor_ref[i], 1e-8); + double gg = this->pw_basis.gg[i]; + double ref = std::max(gg / (gg + gg0), 0.1 / CMtest.mixing_beta); + // rho + EXPECT_NEAR(drhog[i].real(), ref, 1e-10); + EXPECT_NEAR(drhog[i].imag(), ref, 1e-10); + } + for (int i = 0; i < 3*pw_basis.npw; ++i) + { + EXPECT_NEAR(drhog[i + pw_basis.npw].real(), 1.0, 1e-10); + EXPECT_NEAR(drhog[i + pw_basis.npw].imag(), 1.0, 1e-10); + } + // mixing_gg0 = 1.0, mixing_gg0_mag = 2.0 + CMtest.mixing_gg0 = 1.0; + CMtest.mixing_gg0_mag = 2.0; + CMtest.Kerker_screen_recip(drhog); + double gg1 = std::pow(1.0 * 0.529177, 2); + double gg2 = std::pow(2.0 * 0.529177, 2); + for (int i = 0; i < pw_basis.npw; ++i) + { + double gg = this->pw_basis.gg[i]; + double ref = std::max(gg / (gg + gg1), 0.1 / CMtest.mixing_beta); + // rho + EXPECT_NEAR(drhog[i].real(), ref * ref, 1e-10); + EXPECT_NEAR(drhog[i].imag(), ref * ref, 1e-10); + } + for (int i = 0; i < pw_basis.npw; ++i) + { + double gg = this->pw_basis.gg[i]; + double ref = std::max(gg / (gg + gg2), 0.1 / CMtest.mixing_beta_mag); + // rho + for (int j = 1; j < GlobalV::NSPIN; ++j) + { + EXPECT_NEAR(drhog[i + pw_basis.npw * j].real(), ref, 1e-10); + EXPECT_NEAR(drhog[i + pw_basis.npw * j].imag(), ref, 1e-10); + } } - delete[] drhog; delete[] drhog_old; - delete[] drhor; - delete[] drhor_ref; } -TEST_F(ChargeMixingTest, KerkerScreenRecipNewTest) +TEST_F(ChargeMixingTest, KerkerScreenRealTest) { Charge_Mixing CMtest; CMtest.set_rhopw(&pw_basis, &pw_basis); GlobalC::ucell.tpiba = 1.0; - // for new kerker + // nspin = 1 + GlobalV::NSPIN = 1; + double* drhor = new double[GlobalV::NSPIN*pw_basis.nrxx]; + double* drhor_ref = new double[GlobalV::NSPIN*pw_basis.nrxx]; + for (int i = 0; i < GlobalV::NSPIN*pw_basis.nrxx; ++i) + { + drhor_ref[i] = drhor[i] = 1.0; + } + // no kerker + CMtest.mixing_gg0 = 0.0; + CMtest.Kerker_screen_real(drhor); + for (int i = 0; i < GlobalV::NSPIN*pw_basis.nrxx; ++i) + { + EXPECT_EQ(drhor[i], drhor_ref[i]); + } + delete[] drhor; + delete[] drhor_ref; + + // nspin = 2 GlobalV::NSPIN = 2; CMtest.mixing_gg0 = 0.0; - GlobalV::MIXING_GG0_MAG = 0.0; std::complex* drhog = new std::complex[GlobalV::NSPIN*pw_basis.npw]; std::complex* drhog_old = new std::complex[GlobalV::NSPIN*pw_basis.npw]; - double* drhor = new double[GlobalV::NSPIN*pw_basis.nrxx]; - double* drhor_ref = new double[GlobalV::NSPIN*pw_basis.nrxx]; + drhor = new double[GlobalV::NSPIN*pw_basis.nrxx]; + drhor_ref = new double[GlobalV::NSPIN*pw_basis.nrxx]; for (int i = 0; i < GlobalV::NSPIN*pw_basis.npw; ++i) { drhog_old[i] = drhog[i] = std::complex(1.0, 1.0); } - CMtest.Kerker_screen_recip_new(drhog); + CMtest.Kerker_screen_recip(drhog); // no kerker for (int i = 0; i < GlobalV::NSPIN*pw_basis.npw; ++i) { EXPECT_EQ(drhog[i], drhog_old[i]); @@ -615,7 +706,7 @@ TEST_F(ChargeMixingTest, KerkerScreenRecipNewTest) // RECIPROCAL CMtest.mixing_gg0 = 1.0; GlobalV::MIXING_GG0_MAG = 0.0; - CMtest.Kerker_screen_recip_new(drhog); + CMtest.Kerker_screen_recip(drhog); const double gg0 = std::pow(0.529177, 2); for (int i = 0; i < pw_basis.npw; ++i) { @@ -646,6 +737,7 @@ TEST_F(ChargeMixingTest, KerkerScreenRecipNewTest) delete[] drhog_old; delete[] drhor; delete[] drhor_ref; + } TEST_F(ChargeMixingTest, MixRhoTest)