diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 251b81280a..102ee6e246 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -38,9 +38,9 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, GlobalV::ofs_running<<"mixing_beta: "<< this->mixing_beta <mixing_gg0 <mixing_beta_mag <mixing_ndim <* rhog_in = nullptr; std::complex* rhog_out = nullptr; - if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4) + if (GlobalV::NSPIN == 1) { rhog_in = chr->rhog_save[0]; rhog_out = chr->rhog[0]; @@ -316,8 +315,38 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) }; this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, twobeta_mix, true); } + else if (GlobalV::NSPIN == 4) + { + // I do not merge this part with nspin=2 because the method of nspin=2 is almost done, + // while nspin=4 is not finished yet. I will try more methods for nspin=4 in the future. + rhog_in = chr->rhog_save[0]; + rhog_out = chr->rhog[0]; + const int npw = this->rhopw->npw; + auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, 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 +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = 0; i < npw; ++i) + { + out[i] = in[i] + this->mixing_beta * sres[i]; + } + // magnetism, mx, my, mz +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = npw; i < 4 * npw; ++i) + { + out[i] = in[i] + this->mixing_beta_mag * sres[i]; + } + }; + this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, twobeta_mix, true); + } // can choose inner_product_recip_new1 or inner_product_recip_new2 + // inner_product_recip_new1 is a simple sum + // inner_product_recip_new2 is a hartree-like sum, unit is Ry auto inner_product_new = std::bind(&Charge_Mixing::inner_product_recip_new2, this, std::placeholders::_1, std::placeholders::_2); auto inner_product_old @@ -326,7 +355,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) { this->mixing->cal_coef(this->rho_mdata, inner_product_new); } - else if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4) + else if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4) // nspin=4 can use old inner_product { this->mixing->cal_coef(this->rho_mdata, inner_product_old); } @@ -344,6 +373,8 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) { chr->rhopw->recip2real(chr->rhog[is], chr->rho[is]); } + + // renormalize rho in R-space would induce a error in K-space //chr->renormalize_rho(); // For kinetic energy density @@ -400,7 +431,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr) { double* rhor_in; double* rhor_out; - if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4) + if (GlobalV::NSPIN == 1) { rhor_in = chr->rho_save[0]; rhor_out = chr->rho[0]; @@ -420,19 +451,47 @@ void Charge_Mixing::mix_rho_real(Charge* chr) #ifdef _OPENMP #pragma omp parallel for schedule(static, 256) #endif - for (int i = 0; i < nrxx; ++i) - { - out[i] = in[i] + this->mixing_beta * sres[i]; - } + for (int i = 0; i < nrxx; ++i) + { + out[i] = in[i] + this->mixing_beta * sres[i]; + } // magnetism #ifdef _OPENMP #pragma omp parallel for schedule(static, 256) #endif - for (int i = nrxx; i < 2 * nrxx; ++i) - { - out[i] = in[i] + this->mixing_beta_mag * sres[i]; - } - }; + for (int i = nrxx; i < 2 * nrxx; ++i) + { + out[i] = in[i] + this->mixing_beta_mag * sres[i]; + } + }; + this->mixing->push_data(this->rho_mdata, rhor_in, rhor_out, screen, twobeta_mix, true); + } + else if (GlobalV::NSPIN == 4) + { + // I do not merge this part with nspin=2 because the method of nspin=2 is almost done, + // while nspin=4 is not finished yet. I will try more methods for nspin=4 in the future. + rhor_in = chr->rho_save[0]; + rhor_out = chr->rho[0]; + const int nrxx = this->rhopw->nrxx; + auto screen = std::bind(&Charge_Mixing::Kerker_screen_real, this, std::placeholders::_1); + auto twobeta_mix + = [this, nrxx](double* out, const double* in, const double* sres) { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = 0; i < nrxx; ++i) + { + out[i] = in[i] + this->mixing_beta * sres[i]; + } + // magnetism, mx, my, mz +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = nrxx; i < 4 * nrxx; ++i) + { + out[i] = in[i] + this->mixing_beta_mag * sres[i]; + } + }; this->mixing->push_data(this->rho_mdata, rhor_in, rhor_out, screen, twobeta_mix, true); } @@ -642,14 +701,18 @@ void Charge_Mixing::Kerker_screen_recip_new(std::complex* drhog) for (int is = 0; is < GlobalV::NSPIN; ++is) { // new mixing method only support nspin=2 not nspin=4 - if (is == 1 && GlobalV::NSPIN == 2) + if (is >= 1) { if (GlobalV::MIXING_GG0_MAG <= 0.0001 || GlobalV::MIXING_BETA_MAG <= 0.1) { - for (int ig = 0; ig < this->rhopw->npw; ig++) - { - drhog[is * this->rhopw->npw + ig] *= 1; - } +#ifdef __DEBUG + assert(is == 1); // make sure break works +#endif + double is_mag = GlobalV::NSPIN - 1; + //for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++) + //{ + // drhog[is * this->rhopw->npw + ig] *= 1; + //} break; } fac = GlobalV::MIXING_GG0_MAG; @@ -692,13 +755,17 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) for (int is = 0; is < GlobalV::NSPIN; is++) { - if (is == 1 && GlobalV::NSPIN == 2) + if (is >= 1) { if (GlobalV::MIXING_GG0_MAG <= 0.0001 || GlobalV::MIXING_BETA_MAG <= 0.1) { - for (int ig = 0; ig < this->rhopw->npw; ig++) +#ifdef __DEBUG + assert(is == 1); // make sure break works +#endif + double is_mag = GlobalV::NSPIN - 1; + for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++) { - drhog[is * this->rhopw->npw + ig] *= 0; + drhog[is * this->rhopw->npw + ig] = 0; } break; } diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 327714a5af..95178c3b71 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -301,7 +301,7 @@ void Input::Default(void) mixing_beta = -10; mixing_ndim = 8; mixing_gg0 = 1.00; // use Kerker defaultly - mixing_beta_mag = -10.0; // only set when nspin == 2 + mixing_beta_mag = -10.0; // only set when nspin == 2 || nspin == 4 mixing_gg0_mag = 0.0; // defaultly exclude Kerker from mixing magnetic density mixing_gg0_min = 0.1; // defaultly minimum kerker coefficient mixing_tau = false; @@ -2975,6 +2975,13 @@ void Input::Default_2(void) // jiyy add 2019-08-04 scf_thr_type = 1; } } + + // set nspin with noncolin + if (noncolin || lspinorb) + { + nspin = 4; + } + // mixing parameters if (mixing_beta < 0.0) { @@ -2990,12 +2997,14 @@ void Input::Default_2(void) // jiyy add 2019-08-04 } else if (nspin == 4) // I will add this { - mixing_beta = 0.2; + mixing_beta = 0.4; + mixing_beta_mag = 1.6; + mixing_gg0_mag = 0.0; } } else { - if (nspin == 2 && mixing_beta_mag < 0.0) + if ((nspin == 2 || nspin == 4) && mixing_beta_mag < 0.0) { if (mixing_beta <= 0.4) { @@ -3003,15 +3012,10 @@ void Input::Default_2(void) // jiyy add 2019-08-04 } else { - mixing_beta_mag = 1.6; + mixing_beta_mag = 1.6; // 1.6 can be discussed } } } - // set nspin with noncolin - if (noncolin || lspinorb) - { - nspin = 4; - } } #ifdef __MPI void Input::Bcast() diff --git a/tests/integrate/104_PW_NC_magnetic/result.ref b/tests/integrate/104_PW_NC_magnetic/result.ref index 298931edca..384804bb8d 100644 --- a/tests/integrate/104_PW_NC_magnetic/result.ref +++ b/tests/integrate/104_PW_NC_magnetic/result.ref @@ -1,6 +1,6 @@ -etotref -5866.174730432665 -etotperatomref -2933.0873652163 +etotref -5866.197307043378 +etotperatomref -2933.0986535217 pointgroupref O_h spacegroupref O_h nksibzref 1 -totaltimeref 0.77367 +totaltimeref 2.23 diff --git a/tests/integrate/204_NO_KP_NC_deltaspin/INPUT b/tests/integrate/204_NO_KP_NC_deltaspin/INPUT index c8b5a3e677..63fcf21927 100644 --- a/tests/integrate/204_NO_KP_NC_deltaspin/INPUT +++ b/tests/integrate/204_NO_KP_NC_deltaspin/INPUT @@ -25,6 +25,7 @@ smearing_sigma 0.07 mixing_type pulay mixing_beta 0.3 +mixing_beta_mag 0.3 mixing_gg0 0.0 ks_solver genelpa