From 347be7e38b798bb9c5a88232db34049e3b73fc01 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Thu, 21 Dec 2023 16:55:55 +0800 Subject: [PATCH 01/12] delete Kerker_screen_real_test in charge_mixing.h --- source/module_elecstate/module_charge/charge_mixing.h | 1 - 1 file changed, 1 deletion(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 993c229057..32505177e9 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 From 31fa25243225e0d431ee08e19bf8ba5ec13ff8b3 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Thu, 21 Dec 2023 16:57:24 +0800 Subject: [PATCH 02/12] delete auto_set in charge_mixing.h --- .../module_elecstate/module_charge/charge_mixing.h | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 32505177e9..7162504022 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -88,18 +88,6 @@ class Charge_Mixing 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_); - /** * @brief Get the drho * From d595ed4da56b3834c80bbf25f06fc14c33b2f5c8 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Thu, 21 Dec 2023 17:50:36 +0800 Subject: [PATCH 03/12] refactore set_mixing() and UnitTest --- .../module_charge/charge_mixing.cpp | 35 ++++++++----------- .../module_charge/charge_mixing.h | 14 +++++--- .../test/charge_mixing_test.cpp | 16 ++++----- source/module_esolver/esolver_ks.cpp | 4 ++- .../module_io/test/for_testing_input_conv.h | 4 ++- 5 files changed, 38 insertions(+), 35 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index acece06a04..e8d8825cbe 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -24,14 +24,23 @@ 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_nspin = GlobalV::NSPIN; + if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0 ) this->mixing_nspin = 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_nspin, + 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_nspin, sizeof(double)); } #ifdef USE_PAW diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 7162504022..5ad8653844 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -80,13 +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 + const double& mixing_beta_mag_in, + const double& mixing_gg0_mag_in, + const double& mixing_angle_in); /** * @brief Get the drho @@ -133,16 +137,18 @@ 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_nspin = 1; ///< number of spin used in mixing 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 fb20b2519a..921e442b94 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -282,12 +282,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); @@ -298,19 +298,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.")); } @@ -632,7 +632,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) { @@ -662,7 +662,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) { @@ -748,7 +748,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 a97ebfd9dd..5638551bad 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; } From c64826f345f8114a47053e27a4bcf679cba2d7b2 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Thu, 28 Dec 2023 13:39:48 +0800 Subject: [PATCH 04/12] add local variables _nspin/_mixing_rho_type_num/_mixing_rho_unit_num --- .../module_charge/charge_mixing.cpp | 19 +++++++++++-------- .../module_charge/charge_mixing.h | 7 ++++--- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index e8d8825cbe..dca5c6dd9b 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -38,9 +38,12 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, this->mixing_gg0_mag = mixing_gg0_mag_in; this->mixing_angle = mixing_angle_in; // nspin - this->nspin = GlobalV::NSPIN; - this->mixing_nspin = GlobalV::NSPIN; - if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0 ) this->mixing_nspin = 2; + this->_nspin = GlobalV::NSPIN; + this->_mixing_rho_unit_num = GlobalV::NSPIN; + this->_mixing_rho_type_num = 1; + if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0 ) this->_mixing_rho_unit_num = 2; + // 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 <mixing->init_mixing_data(this->rho_mdata, - this->rhopw->npw * this->mixing_nspin, + this->rhopw->npw * this->_mixing_rho_type_num * this->_mixing_rho_unit_num, sizeof(std::complex)); } else { - this->mixing->init_mixing_data(this->rho_mdata, this->rhopw->nrxx * this->mixing_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 @@ -928,7 +931,7 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) 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 && 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; @@ -996,7 +999,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) @@ -1126,7 +1129,7 @@ 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 +// consider a resize for mixing_angle int resize_tmp = 1; if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0) resize_tmp = 2; diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 5ad8653844..c37cf5ecb6 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -64,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 * @@ -141,8 +141,9 @@ class Charge_Mixing 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_nspin = 1; ///< number of spin used in 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; From 3547d64c80ccfc9a66cbd3b17e44f11336e270af Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Thu, 28 Dec 2023 13:44:11 +0800 Subject: [PATCH 05/12] replace GlobalV::variables by local variables in mix_rho_real() --- source/module_elecstate/module_charge/charge_mixing.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index dca5c6dd9b..72b9646346 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -525,7 +525,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]; @@ -536,7 +536,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 == 2) + else if (this->_nspin == 2) { chr->allocate_rho_mag(); chr->get_rho_mag(); @@ -571,7 +571,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr) chr->get_rho_from_mag(); chr->destroy_rho_mag(); } - 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]; @@ -602,7 +602,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 From efdf84e7ba2bc608f6de3f844a084c7d06f739ce Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Thu, 28 Dec 2023 14:34:38 +0800 Subject: [PATCH 06/12] refactor Kerker_screen_real() using _mixing_rho_type_num/_mixing_rho_unit_num --- .../module_charge/charge_mixing.cpp | 110 ++++++++++-------- 1 file changed, 64 insertions(+), 46 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 72b9646346..30619dd57c 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -906,76 +906,94 @@ 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++) + 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 + int address_g = 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[is * this->rhopw->npw + ig + address_g] = 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[is * this->rhopw->npw + ig + address_g] *= (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 + int address_r = 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 = address_r; ir < address_r + this->rhopw->nrxx * this->_mixing_rho_unit_num; ir++) + { + drhor[ir] -= drhor_filter[ir]; + } } } From d08fed9e0d821ed570d043795fc0983619baf9f7 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Thu, 28 Dec 2023 14:44:55 +0800 Subject: [PATCH 07/12] refactor inner_product_real() using _mixing_rho_type_num/_mixing_rho_unit_num --- source/module_elecstate/module_charge/charge_mixing.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 30619dd57c..6da3672e03 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -1147,14 +1147,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]; } From 588daa428e8396e8324aa647a025e44496019d94 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Thu, 28 Dec 2023 14:52:46 +0800 Subject: [PATCH 08/12] add new EXPECT_EQ in input_test.cpp --- source/module_elecstate/module_charge/charge_mixing.cpp | 2 +- source/module_io/test/input_test.cpp | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 6da3672e03..927115fcf0 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -41,7 +41,7 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, this->_nspin = GlobalV::NSPIN; this->_mixing_rho_unit_num = GlobalV::NSPIN; this->_mixing_rho_type_num = 1; - if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0 ) this->_mixing_rho_unit_num = 2; + 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; diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp index 02e69b98b8..b9751bc5cd 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); From 633c67c3194436641301b00d97a61f8587d2ccfa Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Mon, 8 Jan 2024 11:01:54 +0800 Subject: [PATCH 09/12] replace GlobalV:NSPIN by _nspin --- .../module_elecstate/module_charge/charge_mixing.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 48adbb92f8..bf9ab66987 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -40,7 +40,7 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, // nspin this->_nspin = GlobalV::NSPIN; this->_mixing_rho_unit_num = GlobalV::NSPIN; - this->_mixing_rho_type_num = 1; + 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; @@ -569,10 +569,10 @@ void Charge_Mixing::mix_rho_real(Charge* chr) 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++) { @@ -612,7 +612,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); From e7df22a9d023d95b896a5173aab4b12f1883cba0 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Mon, 8 Jan 2024 11:15:03 +0800 Subject: [PATCH 10/12] replace GlobalV::variables by local variables in mix_rho_recip_new() --- .../module_charge/charge_mixing.cpp | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index bf9ab66987..eb7f1fe252 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -304,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]; @@ -313,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++) { @@ -362,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); } @@ -375,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]; @@ -404,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 @@ -485,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]); } From af8c70637fcffe543c357a08c8612ddb00286778 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Mon, 8 Jan 2024 11:31:04 +0800 Subject: [PATCH 11/12] refactor Kerker_screen_recip_new by using _mixing_rho_type_num & _mixing_rho_unit_num --- .../module_charge/charge_mixing.cpp | 60 ++++++++++--------- 1 file changed, 31 insertions(+), 29 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index eb7f1fe252..2d3e8fdb98 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -912,46 +912,48 @@ void Charge_Mixing::Kerker_screen_recip_new(std::complex* drhog) 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; - - // 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 + int address_g = 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[is * this->rhopw->npw + ig + address_g] *= filter_g; + } } } return; From cd4d7fae9942cb4baf933c9da16e451de85b0c70 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Mon, 8 Jan 2024 12:18:04 +0800 Subject: [PATCH 12/12] move pointer by drhog_tmp & drhor_tmp --- .../module_charge/charge_mixing.cpp | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 2d3e8fdb98..33c30b3e36 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -911,12 +911,13 @@ 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; // over all rho_type for (int irho = 0; irho < this->_mixing_rho_type_num; ++irho) { // skip if other rho_type such as kinetic density is included - int address_g = irho * this->_mixing_rho_unit_num * this->rhopw->npw; + 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) { @@ -952,7 +953,7 @@ void Charge_Mixing::Kerker_screen_recip_new(std::complex* drhog) { 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 + address_g] *= filter_g; + drhog_tmp[is * this->rhopw->npw + ig] *= filter_g; } } } @@ -980,10 +981,12 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) } // implement Kerker for density and magnetization separately double fac, gg0, amin; + std::complex* drhog_tmp; + double* drhor_tmp; for (int irho = 0; irho < this->_mixing_rho_type_num; ++irho) { // skip if other rho_type such as kinetic density is included - int address_g = irho * this->_mixing_rho_unit_num * this->rhopw->npw; + drhog_tmp = drhog + irho * this->_mixing_rho_unit_num * this->rhopw->npw; for (int is = 0; is < this->_mixing_rho_unit_num; is++) { if (is >= 1) @@ -996,7 +999,7 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) double is_mag = this->_mixing_rho_unit_num - 1; for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++) { - drhog[is * this->rhopw->npw + ig + address_g] = 0; + drhog_tmp[is * this->rhopw->npw + ig] = 0; } break; } @@ -1023,7 +1026,7 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) // continue; //} double filter_g = std::max(gg / (gg + gg0), GlobalV::MIXING_GG0_MIN / amin); - drhog[is * this->rhopw->npw + ig + address_g] *= (1 - filter_g); + drhog_tmp[is * this->rhopw->npw + ig] *= (1 - filter_g); } } } @@ -1043,13 +1046,13 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) for (int irho = 0; irho < this->_mixing_rho_type_num; ++irho) { // skip if other rho_type such as kinetic density is included - int address_r = irho * this->_mixing_rho_unit_num * this->rhopw->nrxx; + 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 = address_r; ir < address_r + this->rhopw->nrxx * this->_mixing_rho_unit_num; ir++) + for (int ir = 0; ir < this->rhopw->nrxx * this->_mixing_rho_unit_num; ir++) { - drhor[ir] -= drhor_filter[ir]; + drhor_tmp[ir] -= drhor_filter[ir]; } } }