From 8959177daa83a34ac49300b3b1200428d0b45731 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Mon, 5 Feb 2024 16:23:07 +0800 Subject: [PATCH 1/5] rename inner_product_recip_new1 as inner_product_recip_simple --- source/module_elecstate/module_charge/charge_mixing.cpp | 2 +- source/module_elecstate/module_charge/charge_mixing.h | 2 +- source/module_elecstate/test/charge_mixing_test.cpp | 4 ++-- 3 files 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 505fac7772..d6a24fb4d0 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -1260,7 +1260,7 @@ double Charge_Mixing::inner_product_recip(std::complex* rho1, std::compl } // a simple inner product -double Charge_Mixing::inner_product_recip_new1(std::complex* rho1, std::complex* rho2) +double Charge_Mixing::inner_product_recip_simple(std::complex* rho1, std::complex* rho2) { double rnorm = 0.0; // consider a resize for mixing_angle diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 3c8cade931..fcc7e8e1f7 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -158,7 +158,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_simple(std::complex* rho1, std::complex* rho2); double inner_product_recip_new2(std::complex* rho1, std::complex* rho2); /** diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index 52eaaeb62f..7dc6a8804f 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -460,7 +460,7 @@ TEST_F(ChargeMixingTest, InnerDotNewTest) CMtest.set_rhopw(&pw_basis, &pw_basis); GlobalV::NSPIN = 1; - // inner_product_recip_new1 + // a simple sum for inner product std::vector> drhog1(pw_basis.npw); std::vector> drhog2(pw_basis.npw); for (int i = 0; i < pw_basis.npw; ++i) @@ -468,7 +468,7 @@ TEST_F(ChargeMixingTest, InnerDotNewTest) drhog1[i] = 1.0; drhog2[i] = double(i); } - double inner = CMtest.inner_product_recip_new1(drhog1.data(), drhog2.data()); + double inner = CMtest.inner_product_recip_simple(drhog1.data(), drhog2.data()); EXPECT_NEAR(inner, 0.5 * pw_basis.npw * (pw_basis.npw - 1), 1e-8); // inner_product_recip_new2 From de249f4de39dec2b135f35f3a658fab804d31190 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Mon, 5 Feb 2024 17:35:55 +0800 Subject: [PATCH 2/5] inner_product_recip_hartree support all cases --- .../module_charge/charge_mixing.cpp | 102 +++++++++++++----- .../module_charge/charge_mixing.h | 2 +- .../test/charge_mixing_test.cpp | 4 +- 3 files changed, 81 insertions(+), 27 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index d6a24fb4d0..b58625fd54 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -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); @@ -1259,9 +1253,12 @@ double Charge_Mixing::inner_product_recip(std::complex* rho1, std::compl return result; } -// a simple inner product +// 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; @@ -1276,19 +1273,48 @@ double Charge_Mixing::inner_product_recip_simple(std::complex* rho1, std #ifdef __MPI Parallel_Reduce::reduce_pool(rnorm); #endif + + ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_simple"); + return rnorm; } // a Hartree-like inner product -double Charge_Mixing::inner_product_recip_new2(std::complex* rhog1, std::complex* rhog2) +double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, std::complex* rhog2) { + 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; - if (GlobalV::NSPIN==2) + // a lambda function for summing the charge density + auto part_of_rho = [&]() { + double sum = 0.0; +#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; + return sum; + }; + + if (GlobalV::NSPIN==1) + { + sum += part_of_rho(); + } + else if (GlobalV::NSPIN==2) + { + // charge density part #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif @@ -1321,34 +1347,60 @@ double Charge_Mixing::inner_product_recip_new2(std::complex* rhog1, std: } 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; } - else if (GlobalV::NSPIN==4 && GlobalV::MIXING_ANGLE > 0) + else if (GlobalV::NSPIN==4) { if (!GlobalV::DOMAG && !GlobalV::DOMAG_Z) { + sum += part_of_rho(); + } + else if (this->mixing_angle <= 0) + { + // sum for tradtional mixing #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif - for (int ig = 0; ig < this->rhopw->npw; ++ig) + for (int ig = 0; ig < this->rhopw->npw; ig++) { - if (this->rhopw->gg[ig] < 1e-8) + 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 + 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) + { + 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 + 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 + else if (this->mixing_angle > 0) { - // another part with magnetization + // sum for angle mixing #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif @@ -1386,6 +1438,8 @@ double Charge_Mixing::inner_product_recip_new2(std::complex* rhog1, std: Parallel_Reduce::reduce_pool(sum); #endif + ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_hartree"); + sum *= GlobalC::ucell.omega * 0.5; return sum; diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index fcc7e8e1f7..059d28e6c7 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -159,7 +159,7 @@ class Charge_Mixing */ double inner_product_recip(std::complex* rho1, std::complex* rho2); double inner_product_recip_simple(std::complex* rho1, std::complex* rho2); - double inner_product_recip_new2(std::complex* rho1, std::complex* rho2); + double inner_product_recip_hartree(std::complex* rho1, std::complex* rho2); /** * @brief Inner product of two double vectors diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index 7dc6a8804f..63daae6d16 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -491,10 +491,10 @@ TEST_F(ChargeMixingTest, InnerDotNewTest) 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()); + 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_new2(drhog1_mag.data(), drhog2_mag.data()); + inner = CMtest.inner_product_recip_hartree(drhog1_mag.data(), drhog2_mag.data()); EXPECT_NEAR(inner, 236763.82650318215 * 2, 1e-8); } From 67423df487c0f7bc773464f39304958468654dab Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Tue, 6 Feb 2024 10:42:14 +0800 Subject: [PATCH 3/5] add a UnitTest to protect inner_product_real() --- .../test/charge_mixing_test.cpp | 72 +++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index 63daae6d16..783cb389dc 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -269,6 +269,78 @@ 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; + CMtest.set_rhopw(&pw_basis, &pw_basis); + GlobalV::NSPIN = 1; + + // a simple sum for inner product + 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_simple(drhog1.data(), drhog2.data()); + EXPECT_NEAR(inner, 0.5 * pw_basis.npw * (pw_basis.npw - 1), 1e-8); +} + TEST_F(ChargeMixingTest, KerkerScreenRecipTest) { Charge_Mixing CMtest; From 16ec0fb1fbc12a91f207925e6048deb51b7435c9 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Tue, 6 Feb 2024 11:31:27 +0800 Subject: [PATCH 4/5] add a UnitTest to protect inner_product_recip_hartree --- .../test/charge_mixing_test.cpp | 131 +++++++++++++++++- 1 file changed, 126 insertions(+), 5 deletions(-) diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index 783cb389dc..0d13cd05a4 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -326,19 +326,140 @@ TEST_F(ChargeMixingTest, InnerDotRealTest) 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 = 1; + GlobalV::NSPIN = 2; // a simple sum for inner product - std::vector> drhog1(pw_basis.npw); - std::vector> drhog2(pw_basis.npw); - for (int i = 0; i < pw_basis.npw; ++i) + 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 * (pw_basis.npw - 1), 1e-8); + 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, KerkerScreenRecipTest) From 382ede8dc02f292c0052e87ae5ffcaca9fb220e8 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Tue, 6 Feb 2024 13:01:33 +0800 Subject: [PATCH 5/5] rename inner_product_recip_rho and add a UnitTest to protect it --- .../module_charge/charge_mixing.cpp | 275 +++++++++--------- .../module_charge/charge_mixing.h | 10 +- .../test/charge_mixing_test.cpp | 189 +++++------- 3 files changed, 213 insertions(+), 261 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index b58625fd54..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 { @@ -1238,19 +1238,144 @@ 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]; + ModuleBase::TITLE("Charge_Mixing", "inner_product_recip_rho"); + ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_rho"); + + std::complex** rhog1 = new std::complex*[GlobalV::NSPIN]; + std::complex** rhog2 = 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; + rhog1[is] = rho1 + is * this->rhopw->npw; + rhog2[is] = rho2 + is * this->rhopw->npw; + } + + 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; + + auto part_of_noncolin = [&]() + { + double sum = 0.0; +#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]) * 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; + + if (GlobalV::GAMMA_ONLY_PW) + { + sum *= 2.0; + } + + // (2) Second part of density error. + // including |G|=0 term. + double sum2 = 0.0; + + sum2 += fac2 * (conj(rhog1[0][0] - rhog1[1][0]) * (rhog2[0][0] - rhog2[1][0])).real(); + + double mag = 0.0; +#ifdef _OPENMP +#pragma omp parallel for reduction(+ : mag) +#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 *= fac2; + + // if(GlobalV::GAMMA_ONLY_PW); + if (GlobalV::GAMMA_ONLY_PW) // Peize Lin delete ; 2020.01.31 + { + mag *= 2.0; + } + + // std::cout << " sum=" << sum << " mag=" << mag << std::endl; + sum2 += mag; + sum += sum2; + break; + } + case 4: + // non-collinear spin, added by zhengdy + if (!GlobalV::DOMAG && !GlobalV::DOMAG_Z) + sum += part_of_noncolin(); + else + { + // another part with magnetization +#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[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[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) + { + 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[1][ig]) * rhog2[1][ig]).real() + (conj(rhog1[2][ig]) * rhog2[2][ig]).real() + + (conj(rhog1[3][ig]) * rhog2[3][ig]).real()); + } + } + break; } - double result = this->rhog_dot_product(rho1_2d, rho2_2d); - delete[] rho1_2d; - delete[] rho2_2d; - return result; +#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; } // a simple inner product, now is not used anywhere. For test only. @@ -1465,136 +1590,6 @@ double Charge_Mixing::inner_product_real(double* rho1, double* rho2) return rnorm; } -double Charge_Mixing::rhog_dot_product(const std::complex* const* const rhog1, - const std::complex* const* const rhog2) const -{ - ModuleBase::TITLE("Charge_Mixing", "rhog_dot_product"); - ModuleBase::timer::tick("Charge_Mixing", "rhog_dot_product"); - 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; - - auto part_of_noncolin = [&]() // Peize Lin change goto to function at 2020.01.31 - { - double sum = 0.0; -#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]) * 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; - - if (GlobalV::GAMMA_ONLY_PW) - { - sum *= 2.0; - } - - // (2) Second part of density error. - // including |G|=0 term. - double sum2 = 0.0; - - sum2 += fac2 * (conj(rhog1[0][0] - rhog1[1][0]) * (rhog2[0][0] - rhog2[1][0])).real(); - - double mag = 0.0; -#ifdef _OPENMP -#pragma omp parallel for reduction(+ : mag) -#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 *= fac2; - - // if(GlobalV::GAMMA_ONLY_PW); - if (GlobalV::GAMMA_ONLY_PW) // Peize Lin delete ; 2020.01.31 - { - mag *= 2.0; - } - - // std::cout << " sum=" << sum << " mag=" << mag << std::endl; - sum2 += mag; - sum += sum2; - break; - } - case 4: - // non-collinear spin, added by zhengdy - if (!GlobalV::DOMAG && !GlobalV::DOMAG_Z) - sum += part_of_noncolin(); - else - { - // another part with magnetization -#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[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[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) - { - 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[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", "rhog_dot_product"); - - sum *= GlobalC::ucell.omega * 0.5; - - return sum; -} - 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 059d28e6c7..1bf3282359 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -155,9 +155,12 @@ 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_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); @@ -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 0d13cd05a4..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) @@ -462,6 +463,75 @@ TEST_F(ChargeMixingTest, InnerDotRecipHartreeTest) 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; @@ -578,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; - - // a simple sum for inner product - 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_simple(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_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); -} - TEST_F(ChargeMixingTest, MixRhoTest) { GlobalV::double_grid = false;