diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 227d9a050a..33c30b3e36 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -24,14 +24,26 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, const int& mixing_ndim_in, const double& mixing_gg0_in, const bool& mixing_tau_in, - const double& mixing_beta_mag_in) + const double& mixing_beta_mag_in, + const double& mixing_gg0_mag_in, + const double& mixing_angle_in) { + // mixing parameters this->mixing_mode = mixing_mode_in; this->mixing_beta = mixing_beta_in; this->mixing_beta_mag = mixing_beta_mag_in; this->mixing_ndim = mixing_ndim_in; this->mixing_gg0 = mixing_gg0_in; this->mixing_tau = mixing_tau_in; + this->mixing_gg0_mag = mixing_gg0_mag_in; + this->mixing_angle = mixing_angle_in; + // nspin + this->_nspin = GlobalV::NSPIN; + this->_mixing_rho_unit_num = GlobalV::NSPIN; + this->_mixing_rho_type_num = 1; // only charge density + if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0 ) this->_mixing_rho_unit_num = 2; // {rho,|m|} + // will open for tau mixing in the future + // if (this->mixing_tau) this->_mixing_rho_type_num = 2; GlobalV::ofs_running<<"\n----------- Double Check Mixing Parameters Begin ------------"<mixing_mode < 0 ) - { - this->mixing->init_mixing_data(this->rho_mdata, - this->rhopw->npw * 2, - sizeof(std::complex)); - } - else - { - this->mixing->init_mixing_data(this->rho_mdata, - this->rhopw->npw * GlobalV::NSPIN, - sizeof(std::complex)); - } + this->mixing->init_mixing_data(this->rho_mdata, + this->rhopw->npw * this->_mixing_rho_type_num * this->_mixing_rho_unit_num, + sizeof(std::complex)); } else { - if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0 ) - { - this->mixing->init_mixing_data(this->rho_mdata, this->rhopw->nrxx * 2, sizeof(double)); - } - else - { - this->mixing->init_mixing_data(this->rho_mdata, this->rhopw->nrxx * GlobalV::NSPIN, sizeof(double)); - } + this->mixing->init_mixing_data(this->rho_mdata, this->rhopw->nrxx * this->_mixing_rho_type_num * this->_mixing_rho_unit_num, sizeof(double)); } #ifdef USE_PAW @@ -308,7 +304,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) auto inner_product_old = std::bind(&Charge_Mixing::inner_product_recip, this, std::placeholders::_1, std::placeholders::_2); - if (GlobalV::NSPIN == 1) + if (this->_nspin == 1) { rhog_in = chr->rhog_save[0]; rhog_out = chr->rhog[0]; @@ -317,17 +313,17 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) this->mixing->cal_coef(this->rho_mdata, inner_product_old); this->mixing->mix_data(this->rho_mdata, rhog_out); } - else if (GlobalV::NSPIN == 2) + else if (this->_nspin == 2) { // 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); + rhog_mag = new std::complex[npw * _nspin]; + rhog_mag_save = new std::complex[npw * _nspin]; + ModuleBase::GlobalFunc::ZEROS(rhog_mag, npw * _nspin); + ModuleBase::GlobalFunc::ZEROS(rhog_mag_save, npw * _nspin); // get rhog_mag[is*ngmc] and rhog_mag_save[is*ngmc] for (int ig = 0; ig < npw; ig++) { @@ -366,7 +362,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) 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++) + for (int is = 0; is < _nspin; is++) { ModuleBase::GlobalFunc::ZEROS(chr->rhog[is], npw); } @@ -379,7 +375,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) delete[] rhog_mag; delete[] rhog_mag_save; } - else if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE <= 0) + else if (this->_nspin == 4 && this->mixing_angle <= 0) { // normal broyden mixing for {rho, mx, my, mz} rhog_in = chr->rhog_save[0]; @@ -408,7 +404,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) this->mixing->cal_coef(this->rho_mdata, inner_product_old); this->mixing->mix_data(this->rho_mdata, rhog_out); } - else if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0) + else if (this->_nspin == 4 && this->mixing_angle > 0) { // special broyden mixing for {rho, |m|} proposed by J. Phys. Soc. Jpn. 82 (2013) 114706 // here only consider the case of mixing_angle = 1, which mean only change |m| and keep angle fixed @@ -489,13 +485,14 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) } // rhog to rho - if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0) + // |m| do not need FFT for angle-mixing method + if (this->_nspin == 4 && this->mixing_angle > 0) { chr->rhopw->recip2real(chr->rhog[0], chr->rho[0]); } else { - for (int is = 0; is < GlobalV::NSPIN; is++) + for (int is = 0; is < this->_nspin; is++) { chr->rhopw->recip2real(chr->rhog[is], chr->rho[is]); } @@ -555,7 +552,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr) { double* rhor_in; double* rhor_out; - if (GlobalV::NSPIN == 1) + if (this->_nspin == 1) { rhor_in = chr->rho_save[0]; rhor_out = chr->rho[0]; @@ -566,17 +563,17 @@ void Charge_Mixing::mix_rho_real(Charge* chr) this->mixing->cal_coef(this->rho_mdata, inner_product); this->mixing->mix_data(this->rho_mdata, rhor_out); } - else if (GlobalV::NSPIN == 2) + else if (this->_nspin == 2) { // 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); + 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); // get rho_mag[is*nnrx] and rho_mag_save[is*nnrx] for (int ir = 0; ir < nrxx; ir++) { @@ -616,7 +613,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr) 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++) + for (int is = 0; is < _nspin; is++) { ModuleBase::GlobalFunc::ZEROS(chr->rho[is], nrxx); //ModuleBase::GlobalFunc::ZEROS(rho_save[is], nrxx); @@ -630,7 +627,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr) delete[] rho_mag; delete[] rho_mag_save; } - else if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE <= 0) + else if (this->_nspin == 4 && this->mixing_angle <= 0) { // normal broyden mixing for {rho, mx, my, mz} rhor_in = chr->rho_save[0]; @@ -661,7 +658,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr) this->mixing->cal_coef(this->rho_mdata, inner_product); this->mixing->mix_data(this->rho_mdata, rhor_out); } - else if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0) + else if (this->_nspin == 4 && this->mixing_angle > 0) { // special broyden mixing for {rho, |m|} proposed by J. Phys. Soc. Jpn. 82 (2013) 114706 // here only consider the case of mixing_angle = 1, which mean only change |m| and keep angle fixed @@ -914,47 +911,50 @@ void Charge_Mixing::Kerker_screen_recip_new(std::complex* drhog) if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1) return; double fac, gg0, amin; + std::complex* drhog_tmp; - // consider a resize for mixing_angle - int resize_tmp = 1; - if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0) resize_tmp = 2; - - // implement Kerker for density and magnetization separately - for (int is = 0; is < GlobalV::NSPIN / resize_tmp; ++is) + // over all rho_type + for (int irho = 0; irho < this->_mixing_rho_type_num; ++irho) { - // new mixing method only support nspin=2 not nspin=4 - if (is >= 1) + // skip if other rho_type such as kinetic density is included + drhog_tmp = drhog + irho * this->_mixing_rho_unit_num * this->rhopw->npw; + // implement Kerker for density and magnetization separately + for (int is = 0; is < this->_mixing_rho_unit_num; ++is) { - if (GlobalV::MIXING_GG0_MAG <= 0.0001 || GlobalV::MIXING_BETA_MAG <= 0.1) + // kerker for magnetization + if (is >= 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; - //for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++) - //{ - // drhog[is * this->rhopw->npw + ig] *= 1; - //} - break; + //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 = this->mixing_gg0_mag; + amin = this->mixing_beta_mag; + } + else + { + fac = this->mixing_gg0; + amin = this->mixing_beta; } - fac = GlobalV::MIXING_GG0_MAG; - amin = GlobalV::MIXING_BETA_MAG; - } - else - { - fac = this->mixing_gg0; - amin = this->mixing_beta; - } - gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2); + gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2); #ifdef _OPENMP #pragma omp parallel for schedule(static, 512) #endif - 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); - drhog[is * this->rhopw->npw + ig] *= filter_g; + 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); + drhog_tmp[is * this->rhopw->npw + ig] *= filter_g; + } } } return; @@ -964,76 +964,96 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) { if (this->mixing_gg0 <= 0.0001 || this->mixing_beta <= 0.1) return; - // consider a resize for mixing_angle - int resize_tmp = 1; - if (GlobalV::NSPIN == 4 && GlobalV::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); - for (int is = 0; is < GlobalV::NSPIN / resize_tmp; ++is) + std::vector> drhog(this->rhopw->npw * this->_mixing_rho_type_num * this->_mixing_rho_unit_num); + std::vector drhor_filter(this->rhopw->nrxx * this->_mixing_rho_type_num * this->_mixing_rho_unit_num); + for (int irho = 0; irho < this->_mixing_rho_type_num; ++irho) { - // Note after this process some G which is higher than Gmax will be filtered. - // Thus we cannot use Kerker_screen_recip(drhog.data()) directly after it. - this->rhopw->real2recip(drhor + is * this->rhopw->nrxx, drhog.data() + is * this->rhopw->npw); + // skip if other rho_type such as kinetic density is included + int address_r = irho * this->_mixing_rho_unit_num * this->rhopw->nrxx; + int address_g = irho * this->_mixing_rho_unit_num * this->rhopw->npw; + for (int is = 0; is < this->_mixing_rho_unit_num; ++is) + { + // Note after this process some G which is higher than Gmax will be filtered. + // Thus we cannot use Kerker_screen_recip(drhog.data()) directly after it. + this->rhopw->real2recip(drhor + is * this->rhopw->nrxx + address_r, drhog.data() + is * this->rhopw->npw + address_g); + } } // implement Kerker for density and magnetization separately double fac, gg0, amin; - for (int is = 0; is < GlobalV::NSPIN / resize_tmp; is++) + std::complex* drhog_tmp; + double* drhor_tmp; + for (int irho = 0; irho < this->_mixing_rho_type_num; ++irho) { - - if (is >= 1) + // skip if other rho_type such as kinetic density is included + drhog_tmp = drhog + irho * this->_mixing_rho_unit_num * this->rhopw->npw; + for (int is = 0; is < this->_mixing_rho_unit_num; is++) { - if (GlobalV::MIXING_GG0_MAG <= 0.0001 || GlobalV::MIXING_BETA_MAG <= 0.1) + if (is >= 1) { + if (this->mixing_gg0_mag <= 0.0001 || this->mixing_beta_mag <= 0.1) + { #ifdef __DEBUG - assert(is == 1); // make sure break works + 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; - for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++) - { - drhog[is * this->rhopw->npw + ig] = 0; + double is_mag = this->_mixing_rho_unit_num - 1; + for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++) + { + drhog_tmp[is * this->rhopw->npw + ig] = 0; + } + break; } - break; + fac = this->mixing_gg0_mag; + amin = this->mixing_beta_mag; } - fac = GlobalV::MIXING_GG0_MAG; - amin = GlobalV::MIXING_BETA_MAG; - } - else - { - fac = this->mixing_gg0; - amin = this->mixing_beta; - } - - gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2); + else + { + fac = this->mixing_gg0; + amin = this->mixing_beta; + } + + gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2); #ifdef _OPENMP #pragma omp parallel for schedule(static, 512) #endif - for (int ig = 0; ig < this->rhopw->npw; ig++) - { - double gg = this->rhopw->gg[ig]; - // I have not decided how to handle gg=0 part, will be changed in future - //if (gg == 0) - //{ - // drhog[is * this->rhopw->npw + ig] *= 0; - // continue; - //} - double filter_g = std::max(gg / (gg + gg0), GlobalV::MIXING_GG0_MIN / amin); - drhog[is * this->rhopw->npw + ig] *= (1 - filter_g); + for (int ig = 0; ig < this->rhopw->npw; ig++) + { + double gg = this->rhopw->gg[ig]; + // I have not decided how to handle gg=0 part, will be changed in future + //if (gg == 0) + //{ + // drhog[is * this->rhopw->npw + ig] *= 0; + // continue; + //} + double filter_g = std::max(gg / (gg + gg0), GlobalV::MIXING_GG0_MIN / amin); + drhog_tmp[is * this->rhopw->npw + ig] *= (1 - filter_g); + } } } + // inverse FT - for (int is = 0; is < GlobalV::NSPIN / resize_tmp; ++is) + for (int irho = 0; irho < this->_mixing_rho_type_num; ++irho) { - this->rhopw->recip2real(drhog.data() + is * this->rhopw->npw, drhor_filter.data() + is * this->rhopw->nrxx); + // skip if other rho_type such as kinetic density is included + int address_r = irho * this->_mixing_rho_unit_num * this->rhopw->nrxx; + int address_g = irho * this->_mixing_rho_unit_num * this->rhopw->npw; + for (int is = 0; is < this->_mixing_rho_unit_num; ++is) + { + this->rhopw->recip2real(drhog.data() + is * this->rhopw->npw + address_g, drhor_filter.data() + is * this->rhopw->nrxx + address_r); + } } + for (int irho = 0; irho < this->_mixing_rho_type_num; ++irho) + { + // skip if other rho_type such as kinetic density is included + drhor_tmp = drhor + irho * this->_mixing_rho_unit_num * this->rhopw->nrxx; #ifdef _OPENMP #pragma omp parallel for schedule(static, 512) #endif - for (int ir = 0; ir < this->rhopw->nrxx * GlobalV::NSPIN / resize_tmp; ir++) - { - drhor[ir] -= drhor_filter[ir]; + for (int ir = 0; ir < this->rhopw->nrxx * this->_mixing_rho_unit_num; ir++) + { + drhor_tmp[ir] -= drhor_filter[ir]; + } } } @@ -1057,7 +1077,7 @@ double Charge_Mixing::inner_product_recip_new1(std::complex* rho1, std:: { double rnorm = 0.0; // consider a resize for mixing_angle - int resize_tmp = 1; +int resize_tmp = 1; if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0) resize_tmp = 2; #ifdef _OPENMP #pragma omp parallel for reduction(+ : rnorm) @@ -1187,14 +1207,11 @@ double Charge_Mixing::inner_product_recip_new2(std::complex* rhog1, std: 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; #ifdef _OPENMP #pragma omp parallel for reduction(+ : rnorm) #endif - for (int ir = 0; ir < this->rhopw->nrxx * GlobalV::NSPIN / resize_tmp; ++ir) + for (int ir = 0; ir < this->rhopw->nrxx * this->_mixing_rho_type_num * this->_mixing_rho_unit_num; ++ir) { rnorm += rho1[ir] * rho2[ir]; } diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 993c229057..c37cf5ecb6 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -56,7 +56,6 @@ class Charge_Mixing * */ void Kerker_screen_real(double* rho); - void Kerker_screen_real_test(double* rho); /** * @brief Inner product of two complex vectors @@ -65,7 +64,7 @@ class Charge_Mixing double inner_product_recip(std::complex* rho1, std::complex* rho2); double inner_product_recip_new1(std::complex* rho1, std::complex* rho2); double inner_product_recip_new2(std::complex* rho1, std::complex* rho2); - + /** * @brief Inner product of two double vectors * @@ -81,25 +80,17 @@ class Charge_Mixing * @param mixing_gg0_in mixing gg0 for Kerker screen * @param mixing_tau_in whether to use tau mixing * @param mixing_beta_mag_in mixing beta for magnetism + * @param mixing_gg0_mag_in mixing gg0 for magnetism + * @param mixing_angle_in mixing angle for non-colinear angle mixing */ void set_mixing(const std::string& mixing_mode_in, const double& mixing_beta_in, const int& mixing_ndim_in, const double& mixing_gg0_in, const bool& mixing_tau_in, - const double& mixing_beta_mag_in); // mohan add mixing_gg0_in 2014-09-27 - - // /** - // * @brief use auto set - // * - // */ - // void need_auto_set(); - - // /** - // * @brief auto set mixing gg0 and mixing_beta - // * - // */ - // void auto_set(const double& bandgap_in, const UnitCell& ucell_); + const double& mixing_beta_mag_in, + const double& mixing_gg0_mag_in, + const double& mixing_angle_in); /** * @brief Get the drho @@ -146,16 +137,19 @@ class Charge_Mixing double mixing_beta = 0.8; ///< mixing beta for density double mixing_beta_mag = 1.6; ///< mixing beta for magnetism int mixing_ndim = 8; ///< mixing ndim for broyden and pulay - double mixing_gg0 = 0.0; ///< mixing gg0 for Kerker screen + double mixing_gg0 = 0.0; ///< Kerker screen gg0 for density + double mixing_gg0_mag = 0.0; ///< Kerker screen gg0 for magnetism + double mixing_angle = -10.0; ///< mixing beta for non-colinear angle mixing bool mixing_tau = false; ///< whether to use tau mixing + int _nspin = 1; ///< number of spin + int _mixing_rho_type_num = 1; ///< how many types of rho stored in mixing_data. only charge density: 1; charge density and kinetic energy density: 2 + int _mixing_rho_unit_num = 1; ///< number of data_unit for each rho_block in mixing_data, depending on nspin and whether mixing_tau is used bool new_e_iteration = true; ModulePW::PW_Basis* rhopw = nullptr; ///< smooth grid ModulePW::PW_Basis* rhodpw = nullptr; ///< dense grid, same as rhopw for ncpp. - // bool autoset = false; - private: double rhog_dot_product(const std::complex* const* const rhog1, const std::complex* const* const rhog2) const; diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index 4233d0d99d..f4b39d123f 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -113,12 +113,12 @@ TEST_F(ChargeMixingTest, SetMixingTest) bool mixingtau = false; GlobalV::SCF_THR_TYPE = 1; std::string mode = "broyden"; - CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); + CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6, 0.0, 0.0); EXPECT_EQ(CMtest.rho_mdata.length, pw_basis.npw); GlobalV::SCF_THR_TYPE = 2; mode = "broyden"; - CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); + CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6, 0.0, 0.0); EXPECT_EQ(CMtest.rho_mdata.length, pw_basis.nrxx); EXPECT_EQ(CMtest.get_mixing_mode(), "broyden"); EXPECT_EQ(CMtest.get_mixing_beta(), 1.0); @@ -129,19 +129,19 @@ TEST_F(ChargeMixingTest, SetMixingTest) mixingtau = true; mode = "plain"; GlobalV::SCF_THR_TYPE = 1; - CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); + CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6, 0.0, 0.0); CMtest.mix_reset(); EXPECT_EQ(CMtest.tau_mdata.length, pw_basis.npw); GlobalV::SCF_THR_TYPE = 2; - CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); + CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6, 0.0, 0.0); CMtest.mix_reset(); EXPECT_EQ(CMtest.tau_mdata.length, pw_basis.nrxx); mode = "nothing"; std::string output; testing::internal::CaptureStdout(); - EXPECT_EXIT(CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6);, ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(CMtest.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6, 0.0, 0.0);, ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("This Mixing mode is not implemended yet,coming soon.")); } @@ -429,7 +429,7 @@ TEST_F(ChargeMixingTest, MixRhoTest) Charge_Mixing CMtest_recip; CMtest_recip.set_rhopw(&pw_basis, &pw_basis); GlobalV::SCF_THR_TYPE = 1; - CMtest_recip.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); + CMtest_recip.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6, 0.0, 0.0); CMtest_recip.mix_reset(); for(int i = 0 ; i < nspin * npw; ++i) { @@ -459,7 +459,7 @@ TEST_F(ChargeMixingTest, MixRhoTest) Charge_Mixing CMtest_real; GlobalV::SCF_THR_TYPE = 2; CMtest_real.set_rhopw(&pw_basis, &pw_basis); - CMtest_real.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); + CMtest_real.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6, 0.0, 0.0); CMtest_real.mix_reset(); for(int i = 0 ; i < nspin * nrxx; ++i) { @@ -545,7 +545,7 @@ TEST_F(ChargeMixingTest, MixDoubleGridRhoTest) Charge_Mixing CMtest_recip; CMtest_recip.set_rhopw(&pw_basis, &pw_dbasis); GlobalV::SCF_THR_TYPE = 1; - CMtest_recip.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6); + CMtest_recip.set_mixing(mode, beta, dim, gg0, mixingtau, 1.6, 0.0, 0.0); CMtest_recip.mix_reset(); for (int i = 0; i < nspin * npw; ++i) { diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index d9877a0969..0df1580bc5 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -79,7 +79,9 @@ namespace ModuleESolver GlobalV::MIXING_NDIM, GlobalV::MIXING_GG0, GlobalV::MIXING_TAU, - GlobalV::MIXING_BETA_MAG); + GlobalV::MIXING_BETA_MAG, + GlobalV::MIXING_GG0_MAG, + GlobalV::MIXING_ANGLE); // I use default value to replace autoset // using bandgap to auto set mixing_beta // if (std::abs(GlobalV::MIXING_BETA + 10.0) < 1e-6) diff --git a/source/module_io/test/for_testing_input_conv.h b/source/module_io/test/for_testing_input_conv.h index b444732606..597dbce282 100644 --- a/source/module_io/test/for_testing_input_conv.h +++ b/source/module_io/test/for_testing_input_conv.h @@ -271,7 +271,9 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, const int& mixing_ndim_in, const double& mixing_gg0_in, const bool& mixing_tau_in, - const double& mixing_beta_mag_in) + const double& mixing_beta_mag_in, + const double& mixing_gg0_mag_in, + const double& mixing_angle_in) { return; } diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp index 8b117c9f1e..77e0d1be42 100644 --- a/source/module_io/test/input_test.cpp +++ b/source/module_io/test/input_test.cpp @@ -158,6 +158,9 @@ TEST_F(InputTest, Default) EXPECT_DOUBLE_EQ(INPUT.mixing_beta,-10.0); EXPECT_EQ(INPUT.mixing_ndim,8); EXPECT_DOUBLE_EQ(INPUT.mixing_gg0,1.00); + EXPECT_DOUBLE_EQ(INPUT.mixing_gg0_mag,0.00); + EXPECT_DOUBLE_EQ(INPUT.mixing_beta_mag,-10.0); + EXPECT_DOUBLE_EQ(INPUT.mixing_angle,-10.0); EXPECT_EQ(INPUT.init_wfc,"atomic"); EXPECT_EQ(INPUT.mem_saver,0); EXPECT_EQ(INPUT.printe,100);