From 93badfa871a42ccb09906f486d3b7fd616d63b8c Mon Sep 17 00:00:00 2001 From: Yu Liu <77716030+YuLiu98@users.noreply.github.com> Date: Wed, 28 Aug 2024 10:18:22 +0800 Subject: [PATCH] Fix: bug in charge extrapolation (#5007) * Fix: bug in charge extrapolation * Test: update result.ref --- .../module_charge/charge_extra.cpp | 91 ++++++++++--------- .../module_charge/charge_extra.h | 18 ++-- .../test/charge_extra_test.cpp | 22 ++--- source/module_esolver/esolver_fp.cpp | 5 +- tests/integrate/108_PW_RE_NEW/result.ref | 11 ++- tests/integrate/108_PW_RE_PINT_RKS/result.ref | 10 +- 6 files changed, 87 insertions(+), 70 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_extra.cpp b/source/module_elecstate/module_charge/charge_extra.cpp index 22b1d6cdee..fc40a2d494 100644 --- a/source/module_elecstate/module_charge/charge_extra.cpp +++ b/source/module_elecstate/module_charge/charge_extra.cpp @@ -20,11 +20,7 @@ Charge_Extra::~Charge_Extra() } } -void Charge_Extra::Init_CE(const int& nspin, - const int& natom, - const double& volume, - const int& nrxx, - const std::string chg_extrap) +void Charge_Extra::Init_CE(const int& nspin, const int& natom, const int& nrxx, const std::string chg_extrap) { if (chg_extrap == "none") { @@ -47,13 +43,13 @@ void Charge_Extra::Init_CE(const int& nspin, ModuleBase::WARNING_QUIT("Charge_Extra","charge extrapolation method is not available !"); } - this->omega_old = volume; this->nspin = nspin; - if (pot_order > 1) + if (pot_order > 0) { delta_rho1.resize(this->nspin, std::vector(nrxx, 0.0)); delta_rho2.resize(this->nspin, std::vector(nrxx, 0.0)); + delta_rho3.resize(this->nspin, std::vector(nrxx, 0.0)); } if(pot_order == 3) @@ -110,39 +106,18 @@ void Charge_Extra::extrapolate_charge( // if(lsda || noncolin) rho2zeta(); - double** rho_atom = new double*[this->nspin]; - for (int is = 0; is < this->nspin; is++) - { - rho_atom[is] = new double[chr->rhopw->nrxx]; - } - chr->atomic_rho(this->nspin, omega_old, rho_atom, sf->strucFac, ucell); -#ifdef _OPENMP -#pragma omp parallel for collapse(2) schedule(static, 512) -#endif - for (int is = 0; is < this->nspin; is++) - { - for (int ir = 0; ir < chr->rhopw->nrxx; ir++) - { - chr->rho[is][ir] -= rho_atom[is][ir]; - chr->rho[is][ir] *= omega_old; - } - } - if(rho_extr == 1) { ofs_running << " NEW-OLD atomic charge density approx. for the potential !" << std::endl; - if (pot_order > 1) - { #ifdef _OPENMP -#pragma omp parallel for collapse(2) schedule(static, 512) +#pragma omp parallel for collapse(2) schedule(static, 128) #endif - for (int is = 0; is < this->nspin; is++) + for (int is = 0; is < this->nspin; is++) + { + for (int ir = 0; ir < chr->rhopw->nrxx; ir++) { - for (int ir = 0; ir < chr->rhopw->nrxx; ir++) - { - delta_rho1[is][ir] = chr->rho[is][ir]; - } + chr->rho[is][ir] = delta_rho1[is][ir]; } } } @@ -158,8 +133,6 @@ void Charge_Extra::extrapolate_charge( { for (int ir = 0; ir < chr->rhopw->nrxx; ir++) { - delta_rho2[is][ir] = delta_rho1[is][ir]; - delta_rho1[is][ir] = chr->rho[is][ir]; chr->rho[is][ir] = 2 * delta_rho1[is][ir] - delta_rho2[is][ir]; } } @@ -174,7 +147,6 @@ void Charge_Extra::extrapolate_charge( const double one_add_alpha = 1 + alpha; const double beta_alpha = beta - alpha; - std::vector> delta_rho3(this->nspin, std::vector(chr->rhopw->nrxx, 0.0)); #ifdef _OPENMP #pragma omp parallel for collapse(2) schedule(static, 64) #endif @@ -182,9 +154,6 @@ void Charge_Extra::extrapolate_charge( { for (int ir = 0; ir < chr->rhopw->nrxx; ir++) { - delta_rho3[is][ir] = delta_rho2[is][ir]; - delta_rho2[is][ir] = delta_rho1[is][ir]; - delta_rho1[is][ir] = chr->rho[is][ir]; chr->rho[is][ir] = one_add_alpha * delta_rho1[is][ir] + beta_alpha * delta_rho2[is][ir] - beta * delta_rho3[is][ir]; } @@ -192,6 +161,11 @@ void Charge_Extra::extrapolate_charge( } sf->setup_structure_factor(&ucell, chr->rhopw); + double** rho_atom = new double*[this->nspin]; + for (int is = 0; is < this->nspin; is++) + { + rho_atom[is] = new double[chr->rhopw->nrxx]; + } chr->atomic_rho(this->nspin, ucell.omega, rho_atom, sf->strucFac, ucell); #ifdef _OPENMP #pragma omp parallel for collapse(2) schedule(static, 512) @@ -205,8 +179,6 @@ void Charge_Extra::extrapolate_charge( } } - omega_old = ucell.omega; - for (int is = 0; is < this->nspin; is++) { delete[] rho_atom[is]; @@ -296,4 +268,41 @@ void Charge_Extra::update_all_dis(const UnitCell& ucell) assert(iat == ucell.nat); } return; +} + +void Charge_Extra::update_delta_rho(const UnitCell& ucell, const Charge* chr, const Structure_Factor* sf) +{ + if (pot_order == 0) + { + return; + } + + // obtain the difference between chr->rho and atomic_rho + double** rho_atom = new double*[this->nspin]; + for (int is = 0; is < this->nspin; is++) + { + rho_atom[is] = new double[chr->rhopw->nrxx]; + } + chr->atomic_rho(this->nspin, ucell.omega, rho_atom, sf->strucFac, ucell); + +#ifdef _OPENMP +#pragma omp parallel for collapse(2) schedule(static, 512) +#endif + for (int is = 0; is < this->nspin; is++) + { + for (int ir = 0; ir < chr->rhopw->nrxx; ir++) + { + delta_rho3[is][ir] = delta_rho2[is][ir]; + delta_rho2[is][ir] = delta_rho1[is][ir]; + delta_rho1[is][ir] = chr->rho[is][ir] - rho_atom[is][ir]; + delta_rho1[is][ir] *= ucell.omega; + } + } + + for (int is = 0; is < this->nspin; is++) + { + delete[] rho_atom[is]; + } + delete[] rho_atom; + return; } \ No newline at end of file diff --git a/source/module_elecstate/module_charge/charge_extra.h b/source/module_elecstate/module_charge/charge_extra.h index 4e84e17658..50e87c35c1 100644 --- a/source/module_elecstate/module_charge/charge_extra.h +++ b/source/module_elecstate/module_charge/charge_extra.h @@ -47,15 +47,10 @@ class Charge_Extra * * @param nspin the number of spins * @param natom the number of atoms - * @param volume the volume of the cell * @param nrxx the number of grids * @param chg_extrap the charge extrapolation method */ - void Init_CE(const int& nspin, - const int& natom, - const double& volume, - const int& nrxx, - const std::string chg_extrap); + void Init_CE(const int& nspin, const int& natom, const int& nrxx, const std::string chg_extrap); /** * @brief charge extrapolation method @@ -87,11 +82,19 @@ class Charge_Extra */ void update_all_dis(const UnitCell& ucell); + /** + * @brief update the difference of charge density + * + * @param ucell the cell information + * @param chr the charge density + * @param sf the structure factor + */ + void update_delta_rho(const UnitCell& ucell, const Charge* chr, const Structure_Factor* sf); + private: int istep = 0; ///< the current step int pot_order; ///< the specified charge extrapolation method int rho_extr; ///< the actually used method - double omega_old; ///< the old volume of the last step int nspin; ///< the number of spins ModuleBase::Vector3* dis_old1 = nullptr; ///< dis_old2 = pos_old1 - pos_old2 @@ -100,6 +103,7 @@ class Charge_Extra std::vector> delta_rho1; ///< the last step difference of rho and atomic_rho std::vector> delta_rho2; ///< the second last step difference of rho and atomic_rho + std::vector> delta_rho3; ///< the third last step difference of rho and atomic_rho double alpha; ///< parameter used in the second order extrapolation double beta; ///< parameter used in the second order extrapolation diff --git a/source/module_elecstate/test/charge_extra_test.cpp b/source/module_elecstate/test/charge_extra_test.cpp index d263e5491c..14b98716f4 100644 --- a/source/module_elecstate/test/charge_extra_test.cpp +++ b/source/module_elecstate/test/charge_extra_test.cpp @@ -136,7 +136,7 @@ TEST_F(ChargeExtraTest, InitCEWarningQuit) { GlobalV::chg_extrap = "wwww"; testing::internal::CaptureStdout(); - EXPECT_EXIT(CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap), + EXPECT_EXIT(CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap), ::testing::ExitedWithCode(0), ""); std::string output = testing::internal::GetCapturedStdout(); @@ -146,21 +146,21 @@ TEST_F(ChargeExtraTest, InitCEWarningQuit) TEST_F(ChargeExtraTest, InitCECase1) { GlobalV::chg_extrap = "none"; - CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap); + CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap); EXPECT_EQ(CE.pot_order, 0); } TEST_F(ChargeExtraTest, InitCECase2) { GlobalV::chg_extrap = "atomic"; - CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap); + CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap); EXPECT_EQ(CE.pot_order, 1); } TEST_F(ChargeExtraTest, InitCECase3) { GlobalV::chg_extrap = "first-order"; - CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap); + CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap); EXPECT_EQ(CE.pot_order, 2); EXPECT_NE(CE.delta_rho1.size(), 0); EXPECT_NE(CE.delta_rho2.size(), 0); @@ -169,7 +169,7 @@ TEST_F(ChargeExtraTest, InitCECase3) TEST_F(ChargeExtraTest, InitCECase4) { GlobalV::chg_extrap = "second-order"; - CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap); + CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap); EXPECT_EQ(CE.pot_order, 3); EXPECT_DOUBLE_EQ(CE.alpha, 1.0); EXPECT_DOUBLE_EQ(CE.beta, 0.0); @@ -183,7 +183,7 @@ TEST_F(ChargeExtraTest, InitCECase4) TEST_F(ChargeExtraTest, ExtrapolateChargeCase1) { GlobalV::chg_extrap = "second-order"; - CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap); + CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap); CE.istep = 0; CE.pot_order = 3; @@ -205,7 +205,7 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase1) TEST_F(ChargeExtraTest, ExtrapolateChargeCase2) { GlobalV::chg_extrap = "second-order"; - CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap); + CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap); CE.istep = 1; CE.pot_order = 3; @@ -227,7 +227,7 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase2) TEST_F(ChargeExtraTest, ExtrapolateChargeCase3) { GlobalV::chg_extrap = "second-order"; - CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap); + CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap); CE.istep = 2; CE.pot_order = 3; @@ -249,7 +249,7 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase3) TEST_F(ChargeExtraTest, ExtrapolateChargeCase4) { GlobalV::chg_extrap = "second-order"; - CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap); + CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap); CE.istep = 3; GlobalV::ofs_running.open("log"); @@ -271,7 +271,7 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase4) TEST_F(ChargeExtraTest, UpdateAllDis) { GlobalV::chg_extrap = "second-order"; - CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap); + CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap); CE.istep = 3; for (int i = 0; i < ucell->nat; ++i) { @@ -293,7 +293,7 @@ TEST_F(ChargeExtraTest, UpdateAllDis) TEST_F(ChargeExtraTest, FindAlphaAndBeta) { GlobalV::chg_extrap = "second-order"; - CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap); + CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap); CE.istep = 3; for (int i = 0; i < ucell->nat; ++i) { diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index 0391df06eb..e8b5d92f41 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -106,7 +106,7 @@ void ESolver_FP::before_all_runners(const Input_para& inp, UnitCell& cell) this->print_rhofft(inp, GlobalV::ofs_running); //! 5) initialize the charge extrapolation method if necessary - this->CE.Init_CE(GlobalV::NSPIN, GlobalC::ucell.nat, GlobalC::ucell.omega, this->pw_rhod->nrxx, inp.chg_extrap); + this->CE.Init_CE(GlobalV::NSPIN, GlobalC::ucell.nat, this->pw_rhod->nrxx, inp.chg_extrap); return; } @@ -120,6 +120,9 @@ void ESolver_FP::after_scf(const int istep) // 1) write fermi energy ModuleIO::output_efermi(this->conv_elec, this->pelec->eferm.ef); + // 2) update delta rho for charge extrapolation + CE.update_delta_rho(GlobalC::ucell, &(this->chr), &(this->sf)); + if (istep % PARAM.inp.out_interval == 0) { // 3) write charge density diff --git a/tests/integrate/108_PW_RE_NEW/result.ref b/tests/integrate/108_PW_RE_NEW/result.ref index 4926b7bbe2..333605dd79 100644 --- a/tests/integrate/108_PW_RE_NEW/result.ref +++ b/tests/integrate/108_PW_RE_NEW/result.ref @@ -1,5 +1,6 @@ -etotref -253.6810444502974633 -etotperatomref -63.4202611126 -totalforceref 0.094780 -totalstressref 2.603647 -totaltimeref 3.99 +etotref -253.6810827494790033 +etotperatomref -63.4202706874 +totalforceref 0.078384 +totalstressref 2.163967 +STRU_ION1_D_pass 0 +totaltimeref 4.77 \ No newline at end of file diff --git a/tests/integrate/108_PW_RE_PINT_RKS/result.ref b/tests/integrate/108_PW_RE_PINT_RKS/result.ref index eb3c924859..2bc7d3f797 100644 --- a/tests/integrate/108_PW_RE_PINT_RKS/result.ref +++ b/tests/integrate/108_PW_RE_PINT_RKS/result.ref @@ -1,5 +1,5 @@ -etotref -253.6810443355996085 -etotperatomref -63.4202610839 -totalforceref 0.094440 -totalstressref 2.593312 -totaltimeref 5.96 +etotref -253.6810827448143186 +etotperatomref -63.4202706862 +totalforceref 0.078204 +totalstressref 2.142930 +totaltimeref 6.65 \ No newline at end of file