diff --git a/source/module_elecstate/module_charge/charge.cpp b/source/module_elecstate/module_charge/charge.cpp index d94a9abb9f..709c49de47 100644 --- a/source/module_elecstate/module_charge/charge.cpp +++ b/source/module_elecstate/module_charge/charge.cpp @@ -258,113 +258,6 @@ void Charge::renormalize_rho(void) return; } -// for real space magnetic density -void Charge::get_rho_mag(void) -{ - ModuleBase::TITLE("Charge", "get_rho_tot_mag"); - - for (int ir = 0; ir < nrxx; ir++) - { - rho_mag[ir] = rho[0][ir] + rho[1][ir]; - rho_mag_save[ir] = rho_save[0][ir] + rho_save[1][ir]; - } - for (int ir = 0; ir < nrxx; ir++) - { - rho_mag[ir + nrxx] = rho[0][ir] - rho[1][ir]; - rho_mag_save[ir + nrxx] = rho_save[0][ir] - rho_save[1][ir]; - } - return; -} - -void Charge::get_rho_from_mag(void) -{ - ModuleBase::TITLE("Charge", "get_rho_from_mag"); - for (int is = 0; is < nspin; is++) - { - ModuleBase::GlobalFunc::ZEROS(rho[is], nrxx); - //ModuleBase::GlobalFunc::ZEROS(rho_save[is], nrxx); - } - for (int ir = 0; ir < nrxx; ir++) - { - rho[0][ir] = 0.5 * (rho_mag[ir] + rho_mag[ir+nrxx]); - rho[1][ir] = 0.5 * (rho_mag[ir] - rho_mag[ir+nrxx]); - } - - return; -} - -void Charge::allocate_rho_mag(void) -{ - rho_mag = new double[nrxx * nspin]; - rho_mag_save = new double[nrxx * nspin]; - - ModuleBase::GlobalFunc::ZEROS(rho_mag, nrxx * nspin); - ModuleBase::GlobalFunc::ZEROS(rho_mag_save, nrxx * nspin); - - return; -} - -void Charge::destroy_rho_mag(void) -{ - delete[] rho_mag; - delete[] rho_mag_save; - - return; -} - -// for reciprocal space magnetic density -void Charge::get_rhog_mag(void) -{ - ModuleBase::TITLE("Charge", "get_rhog_tot_mag"); - - for (int ig = 0; ig < ngmc; ig++) - { - rhog_mag[ig] = rhog[0][ig] + rhog[1][ig]; - rhog_mag_save[ig] = rhog_save[0][ig] + rhog_save[1][ig]; - } - for (int ig = 0; ig < ngmc; ig++) - { - rhog_mag[ig + ngmc] = rhog[0][ig] - rhog[1][ig]; - rhog_mag_save[ig + ngmc] = rhog_save[0][ig] - rhog_save[1][ig]; - } - return; -} - -void Charge::get_rhog_from_mag(void) -{ - ModuleBase::TITLE("Charge", "get_rhog_from_mag"); - for (int is = 0; is < nspin; is++) - { - ModuleBase::GlobalFunc::ZEROS(rhog[is], ngmc); - } - for (int ig = 0; ig < ngmc; ig++) - { - rhog[0][ig] = 0.5 * (rhog_mag[ig] + rhog_mag[ig+ngmc]); - rhog[1][ig] = 0.5 * (rhog_mag[ig] - rhog_mag[ig+ngmc]); - } - - return; -} - -void Charge::allocate_rhog_mag(void) -{ - rhog_mag = new std::complex[ngmc * nspin]; - rhog_mag_save = new std::complex[ngmc * nspin]; - - ModuleBase::GlobalFunc::ZEROS(rhog_mag, ngmc * nspin); - ModuleBase::GlobalFunc::ZEROS(rhog_mag_save, ngmc * nspin); - - return; -} - -void Charge::destroy_rhog_mag(void) -{ - delete[] rhog_mag; - delete[] rhog_mag_save; - - return; -} - //------------------------------------------------------- // superposition of atomic charges contained in the array // rho_at (read from pseudopotential files) diff --git a/source/module_elecstate/module_charge/charge.h b/source/module_elecstate/module_charge/charge.h index 55bf2ee369..58bc2ef615 100644 --- a/source/module_elecstate/module_charge/charge.h +++ b/source/module_elecstate/module_charge/charge.h @@ -36,12 +36,6 @@ class Charge double **rho = nullptr; double **rho_save = nullptr; - // for magnetic density, onyl support nspin=2 now - double *rho_mag = nullptr; - double *rho_mag_save = nullptr; - std::complex *rhog_mag = nullptr; - std::complex *rhog_mag_save = nullptr; - std::complex **rhog = nullptr; std::complex **rhog_save = nullptr; @@ -89,47 +83,6 @@ class Charge void renormalize_rho(void); - // for magnetic density - /** - * @brief allocate rho_mag[is*nnrx] and rho_mag_save[is*nnrx] - */ - void allocate_rho_mag(void); - - /** - * @brief destroy rho_mag[is*nnrx] and rho_mag_save[is*nnrx] - */ - void destroy_rho_mag(void); - - /** - * @brief get rho_mag[is*nnrx] - */ - void get_rho_mag(void); - - /** - * @brief get rho[is][nnrx] from rho_mag[is*nnrx] - */ - void get_rho_from_mag(void); - - /** - * @brief // allocate rhog_mag[is*ngmc] and rhog_mag_save[is*ngmc] - */ - void allocate_rhog_mag(void); - - /** - * @brief destroy rhog_mag[is*ngmc] and rhog_mag_save[is*ngmc] - */ - void destroy_rhog_mag(void); - - /** - * @brief get rhog_mag[is*ngmc] - */ - void get_rhog_mag(void); - - /** - * @brief get rhog[is][nnrx] from rhog_mag[is*ngmc] - */ - void get_rhog_from_mag(void); - double sum_rho(void) const; void save_rho_before_sum_band(void); diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index acece06a04..227d9a050a 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -319,11 +319,30 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) } else if (GlobalV::NSPIN == 2) { - chr->allocate_rhog_mag(); - chr->get_rhog_mag(); - rhog_in = chr->rhog_mag_save; - rhog_out = chr->rhog_mag; + // magnetic density + std::complex *rhog_mag = nullptr; + std::complex *rhog_mag_save = nullptr; const int npw = this->rhopw->npw; + // allocate rhog_mag[is*ngmc] and rhog_mag_save[is*ngmc] + rhog_mag = new std::complex[npw * GlobalV::NSPIN]; + rhog_mag_save = new std::complex[npw * GlobalV::NSPIN]; + ModuleBase::GlobalFunc::ZEROS(rhog_mag, npw * GlobalV::NSPIN); + ModuleBase::GlobalFunc::ZEROS(rhog_mag_save, npw * GlobalV::NSPIN); + // get rhog_mag[is*ngmc] and rhog_mag_save[is*ngmc] + for (int ig = 0; ig < npw; ig++) + { + rhog_mag[ig] = chr->rhog[0][ig] + chr->rhog[1][ig]; + rhog_mag_save[ig] = chr->rhog_save[0][ig] + chr->rhog_save[1][ig]; + } + for (int ig = 0; ig < npw; ig++) + { + rhog_mag[ig + npw] = chr->rhog[0][ig] - chr->rhog[1][ig]; + rhog_mag_save[ig + npw] = chr->rhog_save[0][ig] - chr->rhog_save[1][ig]; + } + // + rhog_in = rhog_mag_save; + rhog_out = rhog_mag; + // auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, this, std::placeholders::_1); auto twobeta_mix = [this, npw](std::complex* out, const std::complex* in, const std::complex* sres) { @@ -346,9 +365,19 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, twobeta_mix, true); this->mixing->cal_coef(this->rho_mdata, inner_product_new); this->mixing->mix_data(this->rho_mdata, rhog_out); + // get rhog[is][ngmc] from rhog_mag[is*ngmc] + for (int is = 0; is < GlobalV::NSPIN; is++) + { + ModuleBase::GlobalFunc::ZEROS(chr->rhog[is], npw); + } + for (int ig = 0; ig < npw; ig++) + { + chr->rhog[0][ig] = 0.5 * (rhog_mag[ig] + rhog_mag[ig+npw]); + chr->rhog[1][ig] = 0.5 * (rhog_mag[ig] - rhog_mag[ig+npw]); + } // delete - chr->get_rhog_from_mag(); - chr->destroy_rhog_mag(); + delete[] rhog_mag; + delete[] rhog_mag_save; } else if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE <= 0) { @@ -472,9 +501,6 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) } } - // renormalize rho in R-space would induce a error in K-space - //chr->renormalize_rho(); - // For kinetic energy density if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) { @@ -542,11 +568,29 @@ void Charge_Mixing::mix_rho_real(Charge* chr) } else if (GlobalV::NSPIN == 2) { - chr->allocate_rho_mag(); - chr->get_rho_mag(); - rhor_in = chr->rho_mag_save; - rhor_out = chr->rho_mag; + // magnetic density + double *rho_mag = nullptr; + double *rho_mag_save = nullptr; const int nrxx = this->rhopw->nrxx; + // allocate rho_mag[is*nnrx] and rho_mag_save[is*nnrx] + rho_mag = new double[nrxx * GlobalV::NSPIN]; + rho_mag_save = new double[nrxx * GlobalV::NSPIN]; + ModuleBase::GlobalFunc::ZEROS(rho_mag, nrxx * GlobalV::NSPIN); + ModuleBase::GlobalFunc::ZEROS(rho_mag_save, nrxx * GlobalV::NSPIN); + // get rho_mag[is*nnrx] and rho_mag_save[is*nnrx] + for (int ir = 0; ir < nrxx; ir++) + { + rho_mag[ir] = chr->rho[0][ir] + chr->rho[1][ir]; + rho_mag_save[ir] = chr->rho_save[0][ir] + chr->rho_save[1][ir]; + } + for (int ir = 0; ir < nrxx; ir++) + { + rho_mag[ir + nrxx] = chr->rho[0][ir] - chr->rho[1][ir]; + rho_mag_save[ir + nrxx] = chr->rho_save[0][ir] - chr->rho_save[1][ir]; + } + // + rhor_in = rho_mag_save; + rhor_out = rho_mag; 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) { @@ -571,9 +615,20 @@ void Charge_Mixing::mix_rho_real(Charge* chr) = std::bind(&Charge_Mixing::inner_product_real, this, std::placeholders::_1, std::placeholders::_2); this->mixing->cal_coef(this->rho_mdata, inner_product); this->mixing->mix_data(this->rho_mdata, rhor_out); + // get new rho[is][nrxx] from rho_mag[is*nrxx] + for (int is = 0; is < GlobalV::NSPIN; is++) + { + ModuleBase::GlobalFunc::ZEROS(chr->rho[is], nrxx); + //ModuleBase::GlobalFunc::ZEROS(rho_save[is], nrxx); + } + for (int ir = 0; ir < nrxx; ir++) + { + chr->rho[0][ir] = 0.5 * (rho_mag[ir] + rho_mag[ir+nrxx]); + chr->rho[1][ir] = 0.5 * (rho_mag[ir] - rho_mag[ir+nrxx]); + } // delete - chr->get_rho_from_mag(); - chr->destroy_rho_mag(); + delete[] rho_mag; + delete[] rho_mag_save; } else if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE <= 0) { @@ -668,7 +723,6 @@ void Charge_Mixing::mix_rho_real(Charge* chr) delete[] rho_magabs_save; } - chr->renormalize_rho(); double *taur_out, *taur_in; if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) { diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index fb20b2519a..4233d0d99d 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -23,175 +23,6 @@ Charge::Charge() { } -double Charge::sum_rho(void) const -{ - ModuleBase::TITLE("Charge", "sum_rho"); - - double sum_rho = 0.0; - int nspin0 = (nspin == 2) ? 2 : 1; - - for (int is = 0; is < nspin0; is++) - { - for (int ir = 0; ir < nrxx; ir++) - { - if(GlobalV::use_paw) - { - sum_rho += this->rho[is][ir] + this->nhat[is][ir]; - } - else - { - sum_rho += this->rho[is][ir]; - } - } - } - - // multiply the sum of charge density by a factor - sum_rho *= 500.0 / static_cast(this->rhopw->nxyz); - -#ifdef __MPI - Parallel_Reduce::reduce_pool(sum_rho); -#endif - - // mohan fixed bug 2010-01-18, - // sum_rho may be smaller than 1, like Na bcc. - //if (sum_rho <= 0.1) - //{ - //GlobalV::ofs_warning << " sum_rho=" << sum_rho << std::endl; - //ModuleBase::WARNING_QUIT("Charge::renormalize_rho", "Can't find even an electron!"); - //} - - return sum_rho; -} - -void Charge::renormalize_rho(void) -{ - ModuleBase::TITLE("Charge", "renormalize_rho"); - - const double sr = this->sum_rho(); - GlobalV::ofs_warning << std::setprecision(15); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_warning, "charge before normalized", sr); - const double normalize_factor = GlobalV::nelec / sr; - - for (int is = 0; is < nspin; is++) - { - for (int ir = 0; ir < nrxx; ir++) - { - rho[is][ir] *= normalize_factor; - } - } - - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_warning, "charge after normalized", this->sum_rho()); - - GlobalV::ofs_running << std::setprecision(6); - return; -} - -// for real space magnetic density -void Charge::get_rho_mag(void) -{ - ModuleBase::TITLE("Charge", "get_rho_tot_mag"); - - for (int ir = 0; ir < nrxx; ir++) - { - rho_mag[ir] = rho[0][ir] + rho[1][ir]; - rho_mag_save[ir] = rho_save[0][ir] + rho_save[1][ir]; - } - for (int ir = 0; ir < nrxx; ir++) - { - rho_mag[ir + nrxx] = rho[0][ir] - rho[1][ir]; - rho_mag_save[ir + nrxx] = rho_save[0][ir] - rho_save[1][ir]; - } - return; -} - -void Charge::get_rho_from_mag(void) -{ - ModuleBase::TITLE("Charge", "get_rho_from_mag"); - for (int is = 0; is < nspin; is++) - { - ModuleBase::GlobalFunc::ZEROS(rho[is], nrxx); - //ModuleBase::GlobalFunc::ZEROS(rho_save[is], nrxx); - } - for (int ir = 0; ir < nrxx; ir++) - { - rho[0][ir] = 0.5 * (rho_mag[ir] + rho_mag[ir+nrxx]); - rho[1][ir] = 0.5 * (rho_mag[ir] - rho_mag[ir+nrxx]); - } - - return; -} - -void Charge::allocate_rho_mag(void) -{ - rho_mag = new double[nrxx * nspin]; - rho_mag_save = new double[nrxx * nspin]; - - ModuleBase::GlobalFunc::ZEROS(rho_mag, nrxx * nspin); - ModuleBase::GlobalFunc::ZEROS(rho_mag_save, nrxx * nspin); - - return; -} - -void Charge::destroy_rho_mag(void) -{ - delete[] rho_mag; - delete[] rho_mag_save; - - return; -} - -// for reciprocal space magnetic density -void Charge::get_rhog_mag(void) -{ - ModuleBase::TITLE("Charge", "get_rhog_tot_mag"); - - for (int ig = 0; ig < ngmc; ig++) - { - rhog_mag[ig] = rhog[0][ig] + rhog[1][ig]; - rhog_mag_save[ig] = rhog_save[0][ig] + rhog_save[1][ig]; - } - for (int ig = 0; ig < ngmc; ig++) - { - rhog_mag[ig + ngmc] = rhog[0][ig] - rhog[1][ig]; - rhog_mag_save[ig + ngmc] = rhog_save[0][ig] - rhog_save[1][ig]; - } - return; -} - -void Charge::get_rhog_from_mag(void) -{ - ModuleBase::TITLE("Charge", "get_rhog_from_mag"); - for (int is = 0; is < nspin; is++) - { - ModuleBase::GlobalFunc::ZEROS(rhog[is], ngmc); - } - for (int ig = 0; ig < ngmc; ig++) - { - rhog[0][ig] = 0.5 * (rhog_mag[ig] + rhog_mag[ig+ngmc]); - rhog[1][ig] = 0.5 * (rhog_mag[ig] - rhog_mag[ig+ngmc]); - } - - return; -} - -void Charge::allocate_rhog_mag(void) -{ - rhog_mag = new std::complex[ngmc * nspin]; - rhog_mag_save = new std::complex[ngmc * nspin]; - - ModuleBase::GlobalFunc::ZEROS(rhog_mag, ngmc * nspin); - ModuleBase::GlobalFunc::ZEROS(rhog_mag_save, ngmc * nspin); - - return; -} - -void Charge::destroy_rhog_mag(void) -{ - delete[] rhog_mag; - delete[] rhog_mag_save; - - return; -} void Charge::set_rhopw(ModulePW::PW_Basis* rhopw_in) { this->rhopw = rhopw_in; @@ -315,40 +146,6 @@ TEST_F(ChargeMixingTest, SetMixingTest) EXPECT_THAT(output, testing::HasSubstr("This Mixing mode is not implemended yet,coming soon.")); } -/** -TEST_F(ChargeMixingTest, AutoSetTest) -{ - Charge_Mixing CMtest; - CMtest.set_rhopw(&pw_basis, &pw_basis); - GlobalV::SCF_THR_TYPE = 1; - GlobalV::NSPIN = 1; - - CMtest.set_mixing("broyden", 1.0, 1, 0.2, false); - CMtest.auto_set(0.0, GlobalC::ucell); - EXPECT_EQ(CMtest.mixing_beta, 1.0); - EXPECT_EQ(CMtest.mixing_gg0, 0.2); - - CMtest.need_auto_set(); - CMtest.auto_set(0.0, GlobalC::ucell); - EXPECT_EQ(CMtest.mixing_beta, 0.2); - EXPECT_EQ(CMtest.mixing->mixing_beta, 0.2); - EXPECT_EQ(CMtest.mixing_gg0, 1.0); - - CMtest.need_auto_set(); - CMtest.auto_set(1.0, GlobalC::ucell); - EXPECT_EQ(CMtest.mixing_beta, 0.7); - EXPECT_EQ(CMtest.mixing->mixing_beta, 0.7); - EXPECT_EQ(CMtest.mixing_gg0, 1.0); - - GlobalC::ucell.atoms = new Atom[1]; - GlobalC::ucell.ntype = 1; - GlobalC::ucell.atoms[0].ncpp.psd = "Sc"; - CMtest.need_auto_set(); - CMtest.auto_set(1.0, GlobalC::ucell); - EXPECT_EQ(CMtest.mixing_gg0, 1.0); -} -**/ - TEST_F(ChargeMixingTest, KerkerScreenRecipTest) { Charge_Mixing CMtest; diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 97da14c4e7..d9877a0969 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -438,6 +438,7 @@ namespace ModuleESolver // } p_chgmix->mix_rho(pelec->charge); + if (GlobalV::SCF_THR_TYPE == 2) pelec->charge->renormalize_rho(); // renormalize rho in R-space would induce a error in K-space //----------charge mixing done----------- } }