From de24ce3d4d2ec51fd9ed5465a7a7f250615b391f Mon Sep 17 00:00:00 2001 From: sunliang98 <1700011430@pku.edu.cn> Date: Thu, 28 Nov 2024 14:48:57 +0800 Subject: [PATCH 1/2] Refactor: Move if (PARAM.inp.nspin == 4) into cal_ux(). --- source/module_cell/test/unitcell_test.cpp | 2 ++ source/module_elecstate/cal_ux.cpp | 6 ++++++ source/module_elecstate/cal_ux.h | 1 + source/module_esolver/esolver_ks_lcao.cpp | 10 ++-------- source/module_esolver/esolver_ks_lcao_tddft.cpp | 5 +---- source/module_esolver/esolver_ks_pw.cpp | 10 ++-------- source/module_esolver/esolver_of.cpp | 5 +---- source/module_esolver/esolver_of_tool.cpp | 10 ++-------- source/module_esolver/lcao_before_scf.cpp | 5 +---- source/module_esolver/lcao_others.cpp | 5 +---- source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp | 5 +---- .../module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp | 5 +---- 12 files changed, 21 insertions(+), 48 deletions(-) diff --git a/source/module_cell/test/unitcell_test.cpp b/source/module_cell/test/unitcell_test.cpp index a23d14a20f..ea46081986 100644 --- a/source/module_cell/test/unitcell_test.cpp +++ b/source/module_cell/test/unitcell_test.cpp @@ -1035,6 +1035,7 @@ TEST_F(UcellTest, CalUx1) ucell->atoms[0].m_loc_[0].set(0, -1, 0); ucell->atoms[1].m_loc_[0].set(1, 1, 1); ucell->atoms[1].m_loc_[1].set(0, 0, 0); + PARAM.input.nspin = 4; elecstate::cal_ux(*ucell); EXPECT_FALSE(ucell->magnet.lsign_); EXPECT_DOUBLE_EQ(ucell->magnet.ux_[0], 0); @@ -1051,6 +1052,7 @@ TEST_F(UcellTest, CalUx2) ucell->atoms[1].m_loc_[0].set(1, 1, 1); ucell->atoms[1].m_loc_[1].set(0, 0, 0); //(0,0,0) is also parallel to (1,1,1) + PARAM.input.nspin = 4; elecstate::cal_ux(*ucell); EXPECT_TRUE(ucell->magnet.lsign_); EXPECT_NEAR(ucell->magnet.ux_[0], 0.57735, 1e-5); diff --git a/source/module_elecstate/cal_ux.cpp b/source/module_elecstate/cal_ux.cpp index 40e9bdde93..441aacceae 100644 --- a/source/module_elecstate/cal_ux.cpp +++ b/source/module_elecstate/cal_ux.cpp @@ -1,8 +1,14 @@ #include "cal_ux.h" +#include "module_parameter/parameter.h" namespace elecstate { void cal_ux(UnitCell& ucell) { + if (PARAM.inp.nspin != 4) + { + return; + } + double amag, uxmod; int starting_it = 0; int starting_ia = 0; diff --git a/source/module_elecstate/cal_ux.h b/source/module_elecstate/cal_ux.h index 2d00770fb2..912c570bfa 100644 --- a/source/module_elecstate/cal_ux.h +++ b/source/module_elecstate/cal_ux.h @@ -5,6 +5,7 @@ namespace elecstate { + // Only for npsin = 4 void cal_ux(UnitCell& ucell); bool judge_parallel(double a[3], ModuleBase::Vector3 b); diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 437e634d21..4ce45f34c7 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -603,10 +603,7 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const // rho1 and rho2 are the same rho. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - if (PARAM.inp.nspin == 4) - { - elecstate::cal_ux(ucell); - } + elecstate::cal_ux(ucell); //! update the potentials by using new electron charge density this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); @@ -839,10 +836,7 @@ void ESolver_KS_LCAO::update_pot(UnitCell& ucell, const int istep, const if (!this->conv_esolver) { - if (PARAM.inp.nspin == 4) - { - elecstate::cal_ux(ucell); - } + elecstate::cal_ux(ucell); this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); } diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 3c60fb2526..ef231e84e9 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -238,10 +238,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, const int istep, const i // Calculate new potential according to new Charge Density if (!this->conv_esolver) { - if (PARAM.inp.nspin == 4) - { - elecstate::cal_ux(ucell); - } + elecstate::cal_ux(ucell); this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); } diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index f081f5596d..8bc6bf9fcc 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -297,10 +297,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) //! cal_ux should be called before init_scf because //! the direction of ux is used in noncoline_rho - if (PARAM.inp.nspin == 4) - { - elecstate::cal_ux(ucell); - } + elecstate::cal_ux(ucell); //! calculate the total local pseudopotential in real space this->pelec->init_scf(istep, this->sf.strucFac, ucell.symm, (void*)this->pw_wfc); @@ -470,10 +467,7 @@ void ESolver_KS_PW::update_pot(UnitCell& ucell, const int istep, cons { if (!this->conv_esolver) { - if (PARAM.inp.nspin == 4) - { - elecstate::cal_ux(ucell); - } + elecstate::cal_ux(ucell); this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); #ifdef __MPI diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index a5a68ba14a..ebf37a568c 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -314,10 +314,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) void ESolver_OF::update_potential(UnitCell& ucell) { // (1) get dL/dphi - if (PARAM.inp.nspin == 4) - { - elecstate::cal_ux(ucell); - } + elecstate::cal_ux(ucell); this->pelec->pot->update_from_charge(pelec->charge, &ucell); // Hartree + XC + external this->kinetic_potential(pelec->charge->rho, diff --git a/source/module_esolver/esolver_of_tool.cpp b/source/module_esolver/esolver_of_tool.cpp index 85c6b640fe..060a8f4f30 100644 --- a/source/module_esolver/esolver_of_tool.cpp +++ b/source/module_esolver/esolver_of_tool.cpp @@ -133,10 +133,7 @@ void ESolver_OF::cal_potential(double* ptemp_phi, double* rdLdphi, UnitCell& uce } } - if (PARAM.inp.nspin == 4) - { - elecstate::cal_ux(ucell); - } + elecstate::cal_ux(ucell); this->pelec->pot->update_from_charge(this->ptemp_rho_, &ucell); ModuleBase::matrix& vr_eff = this->pelec->pot->get_effective_v(); @@ -173,10 +170,7 @@ void ESolver_OF::cal_dEdtheta(double** ptemp_phi, Charge* temp_rho, UnitCell& uc { double* dphi_dtheta = new double[this->pw_rho->nrxx]; - if (PARAM.inp.nspin == 4) - { - elecstate::cal_ux(ucell); - } + elecstate::cal_ux(ucell); this->pelec->pot->update_from_charge(temp_rho, &ucell); ModuleBase::matrix& vr_eff = this->pelec->pot->get_effective_v(); diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index aa5d623dbc..4bfae34d8e 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -245,10 +245,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) // cal_ux should be called before init_scf because // the direction of ux is used in noncoline_rho //========================================================= - if (PARAM.inp.nspin == 4) - { - elecstate::cal_ux(ucell); - } + elecstate::cal_ux(ucell); // Peize Lin add 2016-12-03 #ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 75f4164e2f..9a41fe5657 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -246,10 +246,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) // cal_ux should be called before init_scf because // the direction of ux is used in noncoline_rho //========================================================= - if (PARAM.inp.nspin == 4) - { - elecstate::cal_ux(ucell); - } + elecstate::cal_ux(ucell); // pelec should be initialized before these calculations this->pelec->init_scf(istep, this->sf.strucFac, ucell.symm); diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp index 417c812059..588dbcc377 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp @@ -71,10 +71,7 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, } else { - if (PARAM.inp.nspin == 4) - { - elecstate::cal_ux(ucell_in); - } + elecstate::cal_ux(ucell_in); const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &ucell_in); // etxc = std::get<0>(etxc_vtxc_v); diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp index 898657febc..3e75afce36 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp @@ -64,10 +64,7 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, } else { - if(PARAM.inp.nspin==4) - { - elecstate::cal_ux(GlobalC::ucell); - } + elecstate::cal_ux(GlobalC::ucell); const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &GlobalC::ucell); // etxc = std::get<0>(etxc_vtxc_v); // may delete? // vtxc = std::get<1>(etxc_vtxc_v); // may delete? From 196dc803f63e2d3e59e2ab85cf532021a5b29a45 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Thu, 28 Nov 2024 08:17:53 +0000 Subject: [PATCH 2/2] [pre-commit.ci lite] apply automatic fixes --- source/module_esolver/esolver_ks_lcao.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 4ce45f34c7..388abb64bc 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -1052,7 +1052,8 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) if ( PARAM.inp.rdmft == true ) // rdmft, added by jghan, 2024-10-17 { ModuleBase::matrix occ_number_ks(this->pelec->wg); - for(int ik=0; ik < occ_number_ks.nr; ++ik) { for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik]; } + for(int ik=0; ik < occ_number_ks.nr; ++ik) { for(int inb=0; inb < occ_number_ks.nc; ++inb) { occ_number_ks(ik, inb) /= this->kv.wk[ik]; +}} this->rdmft_solver.update_elec(occ_number_ks, *(this->psi)); //initialize the gradients of Etotal on occupation numbers and wfc, and set all elements to 0.