diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 505fac7772..8d95fc8b0c 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -230,7 +230,7 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec) } ModuleBase::GlobalFunc::NOTE("Calculate the norm of the Residual std::vector: < R[rho] | R[rho_save] >"); - drho = this->inner_product_recip(drhog.data(), drhog.data()); + drho = this->inner_product_recip_rho(drhog.data(), drhog.data()); } else { @@ -282,13 +282,9 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) divide_data(chr->rhog[0], rhogs_out, rhoghf_out); } - // 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 - = std::bind(&Charge_Mixing::inner_product_recip, this, std::placeholders::_1, std::placeholders::_2); + // inner_product_recip_hartree is a hartree-like sum, unit is Ry + auto inner_product + = std::bind(&Charge_Mixing::inner_product_recip_hartree, this, std::placeholders::_1, std::placeholders::_2); // DIIS Mixing Only for smooth part, while high_frequency part is mixed by plain mixing method. if (GlobalV::NSPIN == 1) @@ -297,7 +293,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) rhog_out = rhogs_out; auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, 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_old); + this->mixing->cal_coef(this->rho_mdata, inner_product); this->mixing->mix_data(this->rho_mdata, rhog_out); } else if (GlobalV::NSPIN == 2) @@ -346,7 +342,7 @@ 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->cal_coef(this->rho_mdata, inner_product); 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++) @@ -397,7 +393,7 @@ 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_old); + this->mixing->cal_coef(this->rho_mdata, inner_product); this->mixing->mix_data(this->rho_mdata, rhog_out); } else if (GlobalV::NSPIN == 4 && GlobalV::MIXING_ANGLE > 0) @@ -461,9 +457,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr) } }; this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, twobeta_mix, true); - auto inner_product_tmp - = std::bind(&Charge_Mixing::inner_product_recip_new2, this, std::placeholders::_1, std::placeholders::_2); - this->mixing->cal_coef(this->rho_mdata, inner_product_tmp); + this->mixing->cal_coef(this->rho_mdata, inner_product); this->mixing->mix_data(this->rho_mdata, rhog_out); // get new |m| in real space using FT this->rhopw->recip2real(rhog_magabs + this->rhopw->npw, rho_magabs); @@ -1244,51 +1238,27 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) } } -double Charge_Mixing::inner_product_recip(std::complex* rho1, std::complex* rho2) +double Charge_Mixing::inner_product_recip_rho(std::complex* rho1, std::complex* rho2) { - std::complex** rho1_2d = new std::complex*[GlobalV::NSPIN]; - std::complex** rho2_2d = new std::complex*[GlobalV::NSPIN]; - for (int is = 0; is < GlobalV::NSPIN; is++) - { - rho1_2d[is] = rho1 + is * this->rhopw->npw; - rho2_2d[is] = rho2 + is * this->rhopw->npw; - } - double result = this->rhog_dot_product(rho1_2d, rho2_2d); - delete[] rho1_2d; - delete[] rho2_2d; - return result; -} + ModuleBase::TITLE("Charge_Mixing", "inner_product_recip_rho"); + ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_rho"); -// a simple inner product -double Charge_Mixing::inner_product_recip_new1(std::complex* rho1, std::complex* 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 ig = 0; ig < this->rhopw->npw * GlobalV::NSPIN / resize_tmp; ++ig) + std::complex** rhog1 = new std::complex*[GlobalV::NSPIN]; + std::complex** rhog2 = new std::complex*[GlobalV::NSPIN]; + for (int is = 0; is < GlobalV::NSPIN; is++) { - rnorm += (conj(rho1[ig]) * rho2[ig]).real(); + rhog1[is] = rho1 + is * this->rhopw->npw; + rhog2[is] = rho2 + is * this->rhopw->npw; } -#ifdef __MPI - Parallel_Reduce::reduce_pool(rnorm); -#endif - return rnorm; -} -// a Hartree-like inner product -double Charge_Mixing::inner_product_recip_new2(std::complex* rhog1, std::complex* rhog2) -{ static const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / GlobalC::ucell.tpiba2; static const double fac2 = ModuleBase::e2 * ModuleBase::FOUR_PI / (ModuleBase::TWO_PI * ModuleBase::TWO_PI); double sum = 0.0; - if (GlobalV::NSPIN==2) + auto part_of_noncolin = [&]() { + double sum = 0.0; #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif @@ -1296,7 +1266,28 @@ double Charge_Mixing::inner_product_recip_new2(std::complex* rhog1, std: { if (this->rhopw->gg[ig] < 1e-8) continue; - sum += (conj(rhog1[ig]) * (rhog2[ig])).real() / this->rhopw->gg[ig]; + sum += (conj(rhog1[0][ig]) * rhog2[0][ig]).real() / this->rhopw->gg[ig]; + } + sum *= fac; + return sum; + }; + + switch (GlobalV::NSPIN) + { + case 1: + sum += part_of_noncolin(); + break; + + case 2: { + // (1) First part of density error. +#ifdef _OPENMP +#pragma omp parallel for reduction(+ : sum) +#endif + for (int ig = 0; ig < this->rhopw->npw; ++ig) + { + if (this->rhopw->gg[ig] < 1e-8) + continue; + sum += (conj(rhog1[0][ig] + rhog1[1][ig]) * (rhog2[0][ig] + rhog2[1][ig])).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -1309,7 +1300,7 @@ double Charge_Mixing::inner_product_recip_new2(std::complex* rhog1, std: // including |G|=0 term. double sum2 = 0.0; - sum2 += fac2 * (conj(rhog1[0 + this->rhopw->npw]) * rhog2[0 + this->rhopw->npw]).real(); + sum2 += fac2 * (conj(rhog1[0][0] - rhog1[1][0]) * (rhog2[0][0] - rhog2[1][0])).real(); double mag = 0.0; #ifdef _OPENMP @@ -1317,7 +1308,7 @@ double Charge_Mixing::inner_product_recip_new2(std::complex* rhog1, std: #endif for (int ig = 0; ig < this->rhopw->npw; ig++) { - mag += (conj(rhog1[ig + this->rhopw->npw]) * rhog2[ig + this->rhopw->npw]).real(); + mag += (conj(rhog1[0][ig] - rhog1[1][ig]) * (rhog2[0][ig] - rhog2[1][ig])).real(); } mag *= fac2; @@ -1330,22 +1321,12 @@ double Charge_Mixing::inner_product_recip_new2(std::complex* rhog1, std: // std::cout << " sum=" << sum << " mag=" << mag << std::endl; sum2 += mag; sum += sum2; + break; } - else if (GlobalV::NSPIN==4 && GlobalV::MIXING_ANGLE > 0) - { + case 4: + // non-collinear spin, added by zhengdy if (!GlobalV::DOMAG && !GlobalV::DOMAG_Z) - { -#ifdef _OPENMP -#pragma omp parallel for reduction(+ : sum) -#endif - for (int ig = 0; ig < this->rhopw->npw; ++ig) - { - if (this->rhopw->gg[ig] < 1e-8) - continue; - sum += (conj(rhog1[ig]) * rhog2[ig]).real() / this->rhopw->gg[ig]; - } - sum *= fac; - } + sum += part_of_noncolin(); else { // another part with magnetization @@ -1356,14 +1337,15 @@ double Charge_Mixing::inner_product_recip_new2(std::complex* rhog1, std: { if (ig == this->rhopw->ig_gge0) continue; - sum += (conj(rhog1[ig]) * rhog2[ig]).real() / this->rhopw->gg[ig]; + sum += (conj(rhog1[0][ig]) * rhog2[0][ig]).real() / this->rhopw->gg[ig]; } sum *= fac; const int ig0 = this->rhopw->ig_gge0; if (ig0 > 0) { sum += fac2 - * ((conj(rhog1[ig0 + this->rhopw->npw]) * rhog2[ig0 + this->rhopw->npw]).real()); + * ((conj(rhog1[1][ig0]) * rhog2[1][ig0]).real() + (conj(rhog1[2][ig0]) * rhog2[2][ig0]).real() + + (conj(rhog1[3][ig0]) * rhog2[3][ig0]).real()); } double fac3 = fac2; if (GlobalV::GAMMA_ONLY_PW) @@ -1378,50 +1360,64 @@ double Charge_Mixing::inner_product_recip_new2(std::complex* rhog1, std: if (ig == ig0) continue; sum += fac3 - * ((conj(rhog1[ig + this->rhopw->npw]) * rhog2[ig + this->rhopw->npw]).real()); + * ((conj(rhog1[1][ig]) * rhog2[1][ig]).real() + (conj(rhog1[2][ig]) * rhog2[2][ig]).real() + + (conj(rhog1[3][ig]) * rhog2[3][ig]).real()); } } + break; } #ifdef __MPI Parallel_Reduce::reduce_pool(sum); #endif + ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_rho"); sum *= GlobalC::ucell.omega * 0.5; + delete[] rhog1; + delete[] rhog2; return sum; } -double Charge_Mixing::inner_product_real(double* rho1, double* rho2) +// a simple inner product, now is not used anywhere. For test only. +double Charge_Mixing::inner_product_recip_simple(std::complex* rho1, std::complex* rho2) { + ModuleBase::TITLE("Charge_Mixing", "inner_product_recip_simple"); + ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_simple"); + 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 ig = 0; ig < this->rhopw->npw * GlobalV::NSPIN / resize_tmp; ++ig) { - rnorm += rho1[ir] * rho2[ir]; + rnorm += (conj(rho1[ig]) * rho2[ig]).real(); } #ifdef __MPI Parallel_Reduce::reduce_pool(rnorm); #endif + + ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_simple"); + return rnorm; } -double Charge_Mixing::rhog_dot_product(const std::complex* const* const rhog1, - const std::complex* const* const rhog2) const +// a Hartree-like inner product +double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, std::complex* rhog2) { - ModuleBase::TITLE("Charge_Mixing", "rhog_dot_product"); - ModuleBase::timer::tick("Charge_Mixing", "rhog_dot_product"); + ModuleBase::TITLE("Charge_Mixing", "inner_product_recip_hartree"); + ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_hartree"); + static const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / GlobalC::ucell.tpiba2; static const double fac2 = ModuleBase::e2 * ModuleBase::FOUR_PI / (ModuleBase::TWO_PI * ModuleBase::TWO_PI); double sum = 0.0; + const int npw = this->rhopw->npw; - auto part_of_noncolin = [&]() // Peize Lin change goto to function at 2020.01.31 + // a lambda function for summing the charge density + auto part_of_rho = [&]() { double sum = 0.0; #ifdef _OPENMP @@ -1431,20 +1427,19 @@ double Charge_Mixing::rhog_dot_product(const std::complex* const* const { if (this->rhopw->gg[ig] < 1e-8) continue; - sum += (conj(rhog1[0][ig]) * rhog2[0][ig]).real() / this->rhopw->gg[ig]; + sum += (conj(rhog1[ig]) * rhog2[ig]).real() / this->rhopw->gg[ig]; } sum *= fac; return sum; }; - - switch (GlobalV::NSPIN) + + if (GlobalV::NSPIN==1) { - case 1: - sum += part_of_noncolin(); - break; - - case 2: { - // (1) First part of density error. + sum += part_of_rho(); + } + else if (GlobalV::NSPIN==2) + { + // charge density part #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif @@ -1452,7 +1447,7 @@ double Charge_Mixing::rhog_dot_product(const std::complex* const* const { if (this->rhopw->gg[ig] < 1e-8) continue; - sum += (conj(rhog1[0][ig] + rhog1[1][ig]) * (rhog2[0][ig] + rhog2[1][ig])).real() / this->rhopw->gg[ig]; + sum += (conj(rhog1[ig]) * (rhog2[ig])).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -1465,7 +1460,7 @@ double Charge_Mixing::rhog_dot_product(const std::complex* const* const // including |G|=0 term. double sum2 = 0.0; - sum2 += fac2 * (conj(rhog1[0][0] - rhog1[1][0]) * (rhog2[0][0] - rhog2[1][0])).real(); + sum2 += fac2 * (conj(rhog1[0 + this->rhopw->npw]) * rhog2[0 + this->rhopw->npw]).real(); double mag = 0.0; #ifdef _OPENMP @@ -1473,28 +1468,27 @@ double Charge_Mixing::rhog_dot_product(const std::complex* const* const #endif for (int ig = 0; ig < this->rhopw->npw; ig++) { - mag += (conj(rhog1[0][ig] - rhog1[1][ig]) * (rhog2[0][ig] - rhog2[1][ig])).real(); + mag += (conj(rhog1[ig + this->rhopw->npw]) * rhog2[ig + this->rhopw->npw]).real(); } mag *= fac2; - // if(GlobalV::GAMMA_ONLY_PW); - if (GlobalV::GAMMA_ONLY_PW) // Peize Lin delete ; 2020.01.31 + if (GlobalV::GAMMA_ONLY_PW) { mag *= 2.0; } - // std::cout << " sum=" << sum << " mag=" << mag << std::endl; sum2 += mag; sum += sum2; - break; } - case 4: - // non-collinear spin, added by zhengdy + else if (GlobalV::NSPIN==4) + { if (!GlobalV::DOMAG && !GlobalV::DOMAG_Z) - sum += part_of_noncolin(); - else { - // another part with magnetization + sum += part_of_rho(); + } + else if (this->mixing_angle <= 0) + { + // sum for tradtional mixing #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif @@ -1502,15 +1496,15 @@ double Charge_Mixing::rhog_dot_product(const std::complex* const* const { if (ig == this->rhopw->ig_gge0) continue; - sum += (conj(rhog1[0][ig]) * rhog2[0][ig]).real() / this->rhopw->gg[ig]; + sum += (conj(rhog1[ig]) * rhog2[ig]).real() / this->rhopw->gg[ig]; } sum *= fac; const int ig0 = this->rhopw->ig_gge0; if (ig0 > 0) { sum += fac2 - * ((conj(rhog1[1][ig0]) * rhog2[1][ig0]).real() + (conj(rhog1[2][ig0]) * rhog2[2][ig0]).real() - + (conj(rhog1[3][ig0]) * rhog2[3][ig0]).real()); + * ((conj(rhog1[ig0 + npw]) * rhog2[ig0 + npw]).real() + (conj(rhog1[ig0 + 2*npw]) * rhog2[ig0 + 2*npw]).real() + + (conj(rhog1[ig0 + 3*npw]) * rhog2[ig0 + 3*npw]).real()); } double fac3 = fac2; if (GlobalV::GAMMA_ONLY_PW) @@ -1525,22 +1519,77 @@ double Charge_Mixing::rhog_dot_product(const std::complex* const* const if (ig == ig0) continue; sum += fac3 - * ((conj(rhog1[1][ig]) * rhog2[1][ig]).real() + (conj(rhog1[2][ig]) * rhog2[2][ig]).real() - + (conj(rhog1[3][ig]) * rhog2[3][ig]).real()); + * ((conj(rhog1[ig + npw]) * rhog2[ig + npw]).real() + (conj(rhog1[ig + 2*npw]) * rhog2[ig + 2*npw]).real() + + (conj(rhog1[ig + 3*npw]) * rhog2[ig + 3*npw]).real()); + } + } + else if (this->mixing_angle > 0) + { + // sum for angle mixing +#ifdef _OPENMP +#pragma omp parallel for reduction(+ : sum) +#endif + for (int ig = 0; ig < this->rhopw->npw; ig++) + { + if (ig == this->rhopw->ig_gge0) + continue; + sum += (conj(rhog1[ig]) * rhog2[ig]).real() / this->rhopw->gg[ig]; + } + sum *= fac; + const int ig0 = this->rhopw->ig_gge0; + if (ig0 > 0) + { + sum += fac2 + * ((conj(rhog1[ig0 + this->rhopw->npw]) * rhog2[ig0 + this->rhopw->npw]).real()); + } + double fac3 = fac2; + if (GlobalV::GAMMA_ONLY_PW) + { + fac3 *= 2.0; + } +#ifdef _OPENMP +#pragma omp parallel for reduction(+ : sum) +#endif + for (int ig = 0; ig < this->rhopw->npw; ig++) + { + if (ig == ig0) + continue; + sum += fac3 + * ((conj(rhog1[ig + this->rhopw->npw]) * rhog2[ig + this->rhopw->npw]).real()); } } - break; } #ifdef __MPI Parallel_Reduce::reduce_pool(sum); #endif - ModuleBase::timer::tick("Charge_Mixing", "rhog_dot_product"); + + ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_hartree"); sum *= GlobalC::ucell.omega * 0.5; return sum; } +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) + { + rnorm += rho1[ir] * rho2[ir]; + } +#ifdef __MPI + Parallel_Reduce::reduce_pool(rnorm); +#endif + return rnorm; +} + void Charge_Mixing::divide_data(std::complex* data_d, std::complex*& data_s, std::complex*& data_hf) diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 3c8cade931..1bf3282359 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -155,11 +155,14 @@ 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 Actually, I am not sure if the definition of inner product for NSPIN=4 is correct, need to be checked. */ - 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); + double inner_product_recip_rho(std::complex* rho1, std::complex* rho2); + double inner_product_recip_simple(std::complex* rho1, std::complex* rho2); + double inner_product_recip_hartree(std::complex* rho1, std::complex* rho2); /** * @brief Inner product of two double vectors @@ -167,9 +170,6 @@ class Charge_Mixing */ double inner_product_real(double* rho1, double* rho2); - double rhog_dot_product(const std::complex* const* const rhog1, - const std::complex* const* const rhog2) const; - /** * @brief divide rho/tau to smooth and high frequency parts * @param data_d dense data diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index 52eaaeb62f..2a62d3c161 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -51,7 +51,7 @@ UnitCell ucell; /** * - Tested Functions: * - SetMixingTest: - * Charge_Mixing::set_mixing(mixing_mode_in,mixing_beta_in,mixing_ndim_in,mixing_gg0_in,mixing_tau_in) + * Charge_Mixing::set_mixing() * Charge_Mixing::init_mixing() * Charge_Mixing::set_rhopw(rhopw_in) * Charge_Mixing::get_mixing_mode() @@ -62,8 +62,9 @@ UnitCell ucell; * - KerkerScreenTest: Charge_Mixing::Kerker_screen_recip(drhog) * Charge_Mixing::Kerker_screen_real(drhog) * - screen drho with Kerker method - * - InnerDotTest: Charge_Mixing::inner_product_recip(rhog1, rhog2) - * Charge_Mixing::rhog_dot_product(rhog1, rhog2) + * - InnerDotTest: Charge_Mixing::inner_product_recip_hartree(rhog1, rhog2) + * Charge_Mixing::inner_product_recip_rho(rhog1, rhog2) + * Charge_Mixing::inner_product_recip_simple(rhog1, rhog2) * Charge_Mixing::inner_product_real(rho1, rho2) * - calculate the inner product of two vectors * - MixRhoTest: Charge_Mixing::mix_rho(chr) @@ -269,6 +270,268 @@ TEST_F(ChargeMixingTest, InitMixingTest) EXPECT_EQ(CMtest.rho_mdata.length, 2 * pw_basis.nrxx); } +TEST_F(ChargeMixingTest, InnerDotRealTest) +{ + Charge_Mixing CMtest; + // non mixing angle case + CMtest.set_mixing(GlobalV::MIXING_MODE, + GlobalV::MIXING_BETA, + GlobalV::MIXING_NDIM, + GlobalV::MIXING_GG0, + GlobalV::MIXING_TAU, + GlobalV::MIXING_BETA_MAG, + GlobalV::MIXING_GG0_MAG, + GlobalV::MIXING_GG0_MIN, + GlobalV::MIXING_ANGLE, + GlobalV::MIXING_DMR); + CMtest.set_rhopw(&pw_basis, &pw_basis); + GlobalV::NSPIN = 4; + + // a simple sum for inner product + std::vector drho1(pw_basis.nrxx * GlobalV::NSPIN); + std::vector drho2(pw_basis.nrxx * GlobalV::NSPIN); + for (int i = 0; i < pw_basis.nrxx * GlobalV::NSPIN; ++i) + { + drho1[i] = 1.0; + drho2[i] = double(i); + } + double inner = CMtest.inner_product_real(drho1.data(), drho2.data()); + EXPECT_NEAR(inner, 0.5 * pw_basis.nrxx * GlobalV::NSPIN * (pw_basis.nrxx * GlobalV::NSPIN - 1), 1e-8); + + // mixing angle case + GlobalV::MIXING_ANGLE = 1.0; + CMtest.set_mixing(GlobalV::MIXING_MODE, + GlobalV::MIXING_BETA, + GlobalV::MIXING_NDIM, + GlobalV::MIXING_GG0, + GlobalV::MIXING_TAU, + GlobalV::MIXING_BETA_MAG, + GlobalV::MIXING_GG0_MAG, + GlobalV::MIXING_GG0_MIN, + GlobalV::MIXING_ANGLE, + GlobalV::MIXING_DMR); + GlobalV::NSPIN = 4; + + // a simple sum for inner product + drho1.resize(pw_basis.nrxx * 2); + drho2.resize(pw_basis.nrxx * 2); + for (int i = 0; i < pw_basis.nrxx * 2; ++i) + { + drho1[i] = 1.0; + drho2[i] = double(i); + } + inner = CMtest.inner_product_real(drho1.data(), drho2.data()); + EXPECT_NEAR(inner, 0.5 * pw_basis.nrxx * 2 * (pw_basis.nrxx * 2 - 1), 1e-8); +} + +TEST_F(ChargeMixingTest, InnerDotRecipSimpleTest) +{ + Charge_Mixing CMtest; + // non mixing angle case + CMtest.set_mixing(GlobalV::MIXING_MODE, + GlobalV::MIXING_BETA, + GlobalV::MIXING_NDIM, + GlobalV::MIXING_GG0, + GlobalV::MIXING_TAU, + GlobalV::MIXING_BETA_MAG, + GlobalV::MIXING_GG0_MAG, + GlobalV::MIXING_GG0_MIN, + GlobalV::MIXING_ANGLE, + GlobalV::MIXING_DMR); + CMtest.set_rhopw(&pw_basis, &pw_basis); + GlobalV::NSPIN = 2; + + // a simple sum for inner product + std::vector> drhog1(pw_basis.npw * GlobalV::NSPIN); + std::vector> drhog2(pw_basis.npw * GlobalV::NSPIN); + for (int i = 0; i < pw_basis.npw * GlobalV::NSPIN; ++i) + { + drhog1[i] = 1.0; + drhog2[i] = double(i); + } + double inner = CMtest.inner_product_recip_simple(drhog1.data(), drhog2.data()); + EXPECT_NEAR(inner, 0.5 * pw_basis.npw * GlobalV::NSPIN * (pw_basis.npw * GlobalV::NSPIN - 1), 1e-8); +} + +TEST_F(ChargeMixingTest, InnerDotRecipHartreeTest) +{ + // REAL + Charge_Mixing CMtest; + CMtest.set_rhopw(&pw_basis, &pw_basis); + const int npw = pw_basis.npw; + const int nrxx = pw_basis.nrxx; + GlobalV::NSPIN = 1; + std::vector drhor1(pw_basis.nrxx); + std::vector drhor2(pw_basis.nrxx); + for (int i = 0; i < pw_basis.nrxx; ++i) + { + drhor1[i] = 1.0; + drhor2[i] = double(i); + } + double inner = CMtest.inner_product_real(drhor1.data(), drhor2.data()); + EXPECT_NEAR(inner, 0.5 * pw_basis.nrxx * (pw_basis.nrxx - 1), 1e-8); + + // RECIPROCAL NSPIN=1 + GlobalC::ucell.tpiba2 = 1.0; + GlobalC::ucell.omega = 2.0; + + GlobalV::NSPIN = 1; + std::vector> drhog1(pw_basis.npw); + std::vector> drhog2(pw_basis.npw); + for (int i = 0; i < pw_basis.nrxx; ++i) + { + drhor1[i] = 0.0; + } + drhor1[2] = 1.0; + pw_basis.real2recip(drhor1.data(), drhog1.data()); + pw_basis.real2recip(drhor2.data(), drhog2.data()); + + inner = CMtest.inner_product_recip_hartree(drhog1.data(), drhog2.data()); + EXPECT_NEAR(inner, -0.3 * ModuleBase::e2 * ModuleBase::FOUR_PI, 1e-8); + + // RECIPROCAL NSPIN=2 + GlobalV::NSPIN = 2; + drhog1.resize(pw_basis.npw * GlobalV::NSPIN); + drhog2.resize(pw_basis.npw * GlobalV::NSPIN); + std::vector> drhog1_mag(pw_basis.npw * GlobalV::NSPIN); + std::vector> drhog2_mag(pw_basis.npw * GlobalV::NSPIN); + for (int i = 0; i < pw_basis.npw * GlobalV::NSPIN; ++i) + { + drhog1[i] = std::complex(1.0, double(i)); + drhog2[i] = std::complex(1.0, 1.0); + } + // set mag + for (int i = 0; i < pw_basis.npw; ++i) + { + drhog1_mag[i] = drhog1[i] + drhog1[i+pw_basis.npw]; + drhog1_mag[i+pw_basis.npw] = drhog1[i] - drhog1[i+pw_basis.npw]; + drhog2_mag[i] = drhog2[i] + drhog2[i+pw_basis.npw]; + drhog2_mag[i+pw_basis.npw] = drhog2[i] - drhog2[i+pw_basis.npw]; + } + GlobalV::GAMMA_ONLY_PW = false; + inner = CMtest.inner_product_recip_hartree(drhog1_mag.data(), drhog2_mag.data()); + EXPECT_NEAR(inner, 236763.82650318215, 1e-8); + GlobalV::GAMMA_ONLY_PW = true; + inner = CMtest.inner_product_recip_hartree(drhog1_mag.data(), drhog2_mag.data()); + EXPECT_NEAR(inner, 236763.82650318215 * 2, 1e-8); + + // RECIPROCAL NSPIN=4 without mixing_angle + GlobalV::NSPIN = 4; + drhog1.resize(pw_basis.npw * GlobalV::NSPIN); + drhog2.resize(pw_basis.npw * GlobalV::NSPIN); + for (int i = 0; i < pw_basis.npw * GlobalV::NSPIN; ++i) + { + drhog1[i] = std::complex(1.0, double(i)); + drhog2[i] = std::complex(1.0, 1.0); + } + + GlobalV::DOMAG = false; + GlobalV::DOMAG_Z = false; + inner = CMtest.inner_product_recip_hartree(drhog1.data(), drhog2.data()); + EXPECT_NEAR(inner, 28260.091995611871, 1e-8); + GlobalV::GAMMA_ONLY_PW = true; + GlobalV::DOMAG = true; + GlobalV::DOMAG_Z = true; + inner = CMtest.inner_product_recip_hartree(drhog1.data(), drhog2.data()); + EXPECT_NEAR(inner, 110668.61166927818, 1e-8); + + // RECIPROCAL NSPIN=4 with mixing_angle + GlobalV::NSPIN = 4; + GlobalV::MIXING_ANGLE = 1.0; + CMtest.set_mixing(GlobalV::MIXING_MODE, + GlobalV::MIXING_BETA, + GlobalV::MIXING_NDIM, + GlobalV::MIXING_GG0, + GlobalV::MIXING_TAU, + GlobalV::MIXING_BETA_MAG, + GlobalV::MIXING_GG0_MAG, + GlobalV::MIXING_GG0_MIN, + GlobalV::MIXING_ANGLE, + GlobalV::MIXING_DMR); + drhog1.resize(pw_basis.npw * 2); + drhog2.resize(pw_basis.npw * 2); + for (int i = 0; i < pw_basis.npw * 2; ++i) + { + drhog1[i] = std::complex(1.0, double(i)); + drhog2[i] = std::complex(1.0, 1.0); + } + GlobalV::GAMMA_ONLY_PW = false; + inner = CMtest.inner_product_recip_hartree(drhog1.data(), drhog2.data()); + EXPECT_NEAR(inner, 36548.881431837777, 1e-8); + GlobalV::GAMMA_ONLY_PW = true; + inner = CMtest.inner_product_recip_hartree(drhog1.data(), drhog2.data()); + EXPECT_NEAR(inner, 44776.555369916401, 1e-8); +} + +TEST_F(ChargeMixingTest, InnerDotRecipRhoTest) +{ + // REAL + Charge_Mixing CMtest; + CMtest.set_rhopw(&pw_basis, &pw_basis); + GlobalV::NSPIN = 1; + std::vector drhor1(pw_basis.nrxx); + std::vector drhor2(pw_basis.nrxx); + for (int i = 0; i < pw_basis.nrxx; ++i) + { + drhor1[i] = 1.0; + drhor2[i] = double(i); + } + double inner = CMtest.inner_product_real(drhor1.data(), drhor2.data()); + EXPECT_NEAR(inner, 0.5 * pw_basis.nrxx * (pw_basis.nrxx - 1), 1e-8); + + // RECIPROCAL + GlobalC::ucell.tpiba2 = 1.0; + GlobalC::ucell.omega = 2.0; + + GlobalV::NSPIN = 1; + std::vector> drhog1(pw_basis.npw); + std::vector> drhog2(pw_basis.npw); + for (int i = 0; i < pw_basis.nrxx; ++i) + { + drhor1[i] = 0.0; + } + drhor1[2] = 1.0; + pw_basis.real2recip(drhor1.data(), drhog1.data()); + pw_basis.real2recip(drhor2.data(), drhog2.data()); + + inner = CMtest.inner_product_recip_rho(drhog1.data(), drhog2.data()); + EXPECT_NEAR(inner, -0.3 * ModuleBase::e2 * ModuleBase::FOUR_PI, 1e-8); + + GlobalV::NSPIN = 2; + drhog1.resize(pw_basis.npw * GlobalV::NSPIN); + drhog2.resize(pw_basis.npw * GlobalV::NSPIN); + for (int i = 0; i < pw_basis.npw * GlobalV::NSPIN; ++i) + { + drhog1[i] = std::complex(1.0, double(i)); + drhog2[i] = std::complex(1.0, 1.0); + } + GlobalV::GAMMA_ONLY_PW = false; + inner = CMtest.inner_product_recip_rho(drhog1.data(), drhog2.data()); + EXPECT_NEAR(inner, 236763.82650318215, 1e-8); + GlobalV::GAMMA_ONLY_PW = true; + inner = CMtest.inner_product_recip_rho(drhog1.data(), drhog2.data()); + EXPECT_NEAR(inner, 236763.82650318215 * 2, 1e-8); + + GlobalV::NSPIN = 4; + drhog1.resize(pw_basis.npw * GlobalV::NSPIN); + drhog2.resize(pw_basis.npw * GlobalV::NSPIN); + for (int i = 0; i < pw_basis.npw * GlobalV::NSPIN; ++i) + { + drhog1[i] = std::complex(1.0, double(i)); + drhog2[i] = std::complex(1.0, 1.0); + } + + GlobalV::DOMAG = false; + GlobalV::DOMAG_Z = false; + inner = CMtest.inner_product_recip_rho(drhog1.data(), drhog2.data()); + EXPECT_NEAR(inner, 28260.091995611871, 1e-8); + GlobalV::GAMMA_ONLY_PW = true; + GlobalV::DOMAG = true; + GlobalV::DOMAG_Z = true; + inner = CMtest.inner_product_recip_rho(drhog1.data(), drhog2.data()); + EXPECT_NEAR(inner, 110668.61166927818, 1e-8); +} + TEST_F(ChargeMixingTest, KerkerScreenRecipTest) { Charge_Mixing CMtest; @@ -385,119 +648,6 @@ TEST_F(ChargeMixingTest, KerkerScreenRecipNewTest) delete[] drhor_ref; } -TEST_F(ChargeMixingTest, InnerDotTest) -{ - // REAL - Charge_Mixing CMtest; - CMtest.set_rhopw(&pw_basis, &pw_basis); - GlobalV::NSPIN = 1; - std::vector drhor1(pw_basis.nrxx); - std::vector drhor2(pw_basis.nrxx); - for (int i = 0; i < pw_basis.nrxx; ++i) - { - drhor1[i] = 1.0; - drhor2[i] = double(i); - } - double inner = CMtest.inner_product_real(drhor1.data(), drhor2.data()); - EXPECT_NEAR(inner, 0.5 * pw_basis.nrxx * (pw_basis.nrxx - 1), 1e-8); - - // RECIPROCAL - GlobalC::ucell.tpiba2 = 1.0; - GlobalC::ucell.omega = 2.0; - - GlobalV::NSPIN = 1; - std::vector> drhog1(pw_basis.npw); - std::vector> drhog2(pw_basis.npw); - for (int i = 0; i < pw_basis.nrxx; ++i) - { - drhor1[i] = 0.0; - } - drhor1[2] = 1.0; - pw_basis.real2recip(drhor1.data(), drhog1.data()); - pw_basis.real2recip(drhor2.data(), drhog2.data()); - - inner = CMtest.inner_product_recip(drhog1.data(), drhog2.data()); - EXPECT_NEAR(inner, -0.3 * ModuleBase::e2 * ModuleBase::FOUR_PI, 1e-8); - - GlobalV::NSPIN = 2; - drhog1.resize(pw_basis.npw * GlobalV::NSPIN); - drhog2.resize(pw_basis.npw * GlobalV::NSPIN); - for (int i = 0; i < pw_basis.npw * GlobalV::NSPIN; ++i) - { - drhog1[i] = std::complex(1.0, double(i)); - drhog2[i] = std::complex(1.0, 1.0); - } - GlobalV::GAMMA_ONLY_PW = false; - inner = CMtest.inner_product_recip(drhog1.data(), drhog2.data()); - EXPECT_NEAR(inner, 236763.82650318215, 1e-8); - GlobalV::GAMMA_ONLY_PW = true; - inner = CMtest.inner_product_recip(drhog1.data(), drhog2.data()); - EXPECT_NEAR(inner, 236763.82650318215 * 2, 1e-8); - - GlobalV::NSPIN = 4; - drhog1.resize(pw_basis.npw * GlobalV::NSPIN); - drhog2.resize(pw_basis.npw * GlobalV::NSPIN); - for (int i = 0; i < pw_basis.npw * GlobalV::NSPIN; ++i) - { - drhog1[i] = std::complex(1.0, double(i)); - drhog2[i] = std::complex(1.0, 1.0); - } - - GlobalV::DOMAG = false; - GlobalV::DOMAG_Z = false; - inner = CMtest.inner_product_recip(drhog1.data(), drhog2.data()); - EXPECT_NEAR(inner, 28260.091995611871, 1e-8); - GlobalV::GAMMA_ONLY_PW = true; - GlobalV::DOMAG = true; - GlobalV::DOMAG_Z = true; - inner = CMtest.inner_product_recip(drhog1.data(), drhog2.data()); - EXPECT_NEAR(inner, 110668.61166927818, 1e-8); -} - -TEST_F(ChargeMixingTest, InnerDotNewTest) -{ - Charge_Mixing CMtest; - CMtest.set_rhopw(&pw_basis, &pw_basis); - GlobalV::NSPIN = 1; - - // inner_product_recip_new1 - std::vector> drhog1(pw_basis.npw); - std::vector> drhog2(pw_basis.npw); - for (int i = 0; i < pw_basis.npw; ++i) - { - drhog1[i] = 1.0; - drhog2[i] = double(i); - } - double inner = CMtest.inner_product_recip_new1(drhog1.data(), drhog2.data()); - EXPECT_NEAR(inner, 0.5 * pw_basis.npw * (pw_basis.npw - 1), 1e-8); - - // inner_product_recip_new2 - GlobalV::NSPIN = 2; - drhog1.resize(pw_basis.npw * GlobalV::NSPIN); - drhog2.resize(pw_basis.npw * GlobalV::NSPIN); - std::vector> drhog1_mag(pw_basis.npw * GlobalV::NSPIN); - std::vector> drhog2_mag(pw_basis.npw * GlobalV::NSPIN); - for (int i = 0; i < pw_basis.npw * GlobalV::NSPIN; ++i) - { - drhog1[i] = std::complex(1.0, double(i)); - drhog2[i] = std::complex(1.0, 1.0); - } - // set mag - for (int i = 0; i < pw_basis.npw; ++i) - { - drhog1_mag[i] = drhog1[i] + drhog1[i+pw_basis.npw]; - drhog1_mag[i+pw_basis.npw] = drhog1[i] - drhog1[i+pw_basis.npw]; - drhog2_mag[i] = drhog2[i] + drhog2[i+pw_basis.npw]; - drhog2_mag[i+pw_basis.npw] = drhog2[i] - drhog2[i+pw_basis.npw]; - } - GlobalV::GAMMA_ONLY_PW = false; - inner = CMtest.inner_product_recip_new2(drhog1_mag.data(), drhog2_mag.data()); - EXPECT_NEAR(inner, 236763.82650318215, 1e-8); - GlobalV::GAMMA_ONLY_PW = true; - inner = CMtest.inner_product_recip_new2(drhog1_mag.data(), drhog2_mag.data()); - EXPECT_NEAR(inner, 236763.82650318215 * 2, 1e-8); -} - TEST_F(ChargeMixingTest, MixRhoTest) { GlobalV::double_grid = false;