diff --git a/source/Makefile.Objects b/source/Makefile.Objects
index 91e4f34bea..75ea579cdd 100644
--- a/source/Makefile.Objects
+++ b/source/Makefile.Objects
@@ -239,7 +239,6 @@ OBJS_ELECSTAT=elecstate.o\
     potential_types.o\
     pot_sep.o\
     pot_local.o\
-    pot_local_paw.o\
     H_Hartree_pw.o\
     H_TDDFT_pw.o\
     pot_xc.o\
diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp
index 1e80391dfc..cc5328c3c9 100644
--- a/source/source_esolver/esolver_double_xc.cpp
+++ b/source/source_esolver/esolver_double_xc.cpp
@@ -51,7 +51,6 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(UnitCell& ucell, const Input_p
         this->pelec_base = new elecstate::ElecStateLCAO<TK>(&(this->chr_base), // use which parameter?
                                                        &(this->kv),
                                                        this->kv.get_nks(),
-                                                       this->pw_rho,
                                                        this->pw_big);
     }    
 
@@ -87,7 +86,10 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(UnitCell& ucell, const Input_p
     this->dmat_base.allocate_dm(&this->kv, &this->pv, PARAM.inp.nspin);
 
     // 10) inititlize the charge density
-    this->chr_base.allocate(PARAM.inp.nspin);
+	this->chr_base.set_rhopw(this->pw_rhod); // mohan add 20251130
+	this->chr_base.allocate(PARAM.inp.nspin);
+	this->chr_base.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv);
+	this->chr_base.check_rho();
 
     // 11) initialize the potential
     if (this->pelec_base->pot == nullptr)
diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp
index afe1fd6e78..6a6f09ed98 100644
--- a/source/source_esolver/esolver_fp.cpp
+++ b/source/source_esolver/esolver_fp.cpp
@@ -41,24 +41,56 @@ void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp)
 {
     ModuleBase::TITLE("ESolver_FP", "before_all_runners");
 
-    //! read pseudopotentials
+    //! 1) read pseudopotentials
     elecstate::read_pseudo(GlobalV::ofs_running, ucell);
 
-    // setup pw_rho, pw_rhod, pw_big, sf, and read_pseudopotentials
+    //! 2) setup pw_rho, pw_rhod, pw_big, sf, and read_pseudopotentials
     pw::setup_pwrho(ucell, PARAM.globalv.double_grid, this->pw_rho_flag, 
-      this->pw_rho, this->pw_rhod, this->pw_big, 
-      this->classname, inp);
+      this->pw_rho, this->pw_rhod, this->pw_big, this->classname, inp);
 
-    // setup structure factors
+    //! 3) setup structure factors
     this->sf.set(this->pw_rhod, inp.nbspline);
 
-    // write geometry file
+    //! 4)  write geometry file
     ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif",
       ucell, "# Generated by ABACUS ModuleIO::CifParser", "data_?");
 
-    // init charge extrapolation
+    //! 5) init charge extrapolation
     this->CE.Init_CE(inp.nspin, ucell.nat, this->pw_rhod->nrxx, inp.chg_extrap);
 
+    //! 6) symmetry analysis should be performed every time the cell is changed
+    if (ModuleSymmetry::Symmetry::symm_flag == 1)
+    {
+        ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
+        ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
+    }
+
+    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
+
+    //! 7) setup k points in the Brillouin zone according to symmetry.
+    this->kv.set(ucell,ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
+    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
+
+    //! 8) print information
+    ModuleIO::print_parameters(ucell, this->kv, inp);
+
+    //! 9) parallel of FFT grid
+    this->Pgrid.init(this->pw_rhod->nx, this->pw_rhod->ny, this->pw_rhod->nz,
+            this->pw_rhod->nplane, this->pw_rhod->nrxx, pw_big->nbz, pw_big->bz);
+
+    //! 10) calculate the structure factor
+    this->sf.setup(&ucell, Pgrid, this->pw_rhod);
+
+    //! 11) setup the xc functional
+    XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
+    GlobalV::ofs_running<<XC_Functional::output_info()<<std::endl;
+
+    //! 11) initialize the charge density, we need to first set xc_type,
+    // then we can call chr.allocate()
+	this->chr.set_rhopw(this->pw_rhod); // mohan add 20251130
+	this->chr.allocate(inp.nspin); // mohan move this from setup_estate_pw, 20251128
+
+
     return;
 }
 
diff --git a/source/source_esolver/esolver_gets.cpp b/source/source_esolver/esolver_gets.cpp
index 4e79849367..108b148f7e 100644
--- a/source/source_esolver/esolver_gets.cpp
+++ b/source/source_esolver/esolver_gets.cpp
@@ -53,7 +53,6 @@ void ESolver_GetS::before_all_runners(UnitCell& ucell, const Input_para& inp)
         this->pelec = new elecstate::ElecStateLCAO<std::complex<double>>(&(this->chr), // use which parameter?
                                                                          &(this->kv),
                                                                          this->kv.get_nks(),
-                                                                         this->pw_rho,
                                                                          this->pw_big);
     }
 
diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp
index f19bbc8bdf..604c1569cc 100644
--- a/source/source_esolver/esolver_ks.cpp
+++ b/source/source_esolver/esolver_ks.cpp
@@ -62,46 +62,26 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
     // cell_factor
     this->ppcell.cell_factor = inp.cell_factor;
 
-    //! 3) setup charge mixing
+    //! 3) setup Exc for the first element '0' (all elements have same exc) 
+//    XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
+//    GlobalV::ofs_running<<XC_Functional::output_info()<<std::endl;
+
+    //! 4) setup charge mixing
     p_chgmix = new Charge_Mixing();
     p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod);
-
-    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
-
-    //! 4) setup Exc for the first element '0' (all elements have same exc) 
-    XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
-    GlobalV::ofs_running<<XC_Functional::output_info()<<std::endl;
-    
-    //! 5) setup the charge mixing parameters
     p_chgmix->set_mixing(inp.mixing_mode, inp.mixing_beta, inp.mixing_ndim,
       inp.mixing_gg0, inp.mixing_tau, inp.mixing_beta_mag, inp.mixing_gg0_mag,
       inp.mixing_gg0_min, inp.mixing_angle, inp.mixing_dmr, ucell.omega, ucell.tpiba);
-
     p_chgmix->init_mixing();
 
-    //! 6) symmetry analysis should be performed every time the cell is changed
-    if (ModuleSymmetry::Symmetry::symm_flag == 1)
-    {
-        ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
-        ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
-    }
-
-    //! 7) setup k points in the Brillouin zone according to symmetry.
-    this->kv.set(ucell,ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
-    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
-
-    //! 8) print information
-    ModuleIO::print_parameters(ucell, this->kv, inp);
-
-    //! 9) setup plane wave for electronic wave functions
+    //! 5) setup plane wave for electronic wave functions
     pw::setup_pwwfc(inp, ucell, *this->pw_rho, this->kv, this->pw_wfc);
 
-    //! 10) parallel of FFT grid 
-	Pgrid.init(this->pw_rhod->nx, this->pw_rhod->ny, this->pw_rhod->nz,
-			this->pw_rhod->nplane, this->pw_rhod->nrxx, pw_big->nbz, pw_big->bz);
-
-    //! 11) calculate the structure factor
-    this->sf.setup(&ucell, Pgrid, this->pw_rhod);
+    //! 6) read in charge density, mohan add 2025-11-28
+    //! Inititlize the charge density.
+    this->chr.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv, this->pw_wfc);
+    this->chr.check_rho(); // check the rho
+  
 }
 
 template <typename T, typename Device>
@@ -115,20 +95,14 @@ void ESolver_KS<T, Device>::hamilt2rho(UnitCell& ucell, const int istep, const i
     this->hamilt2rho_single(ucell, istep, iter, diag_ethr);
 
     // 2) for MPI: STOGROUP? need to rewrite
-    //<Temporary> It may be changed when more clever parallel algorithm is
-    // put forward.
-    // When parallel algorithm for bands are adopted. Density will only be
-    // treated in the first group.
-    //(Different ranks should have abtained the same, but small differences
-    // always exist in practice.)
+    //<Temporary> It may be changed when more clever parallel algorithm is put forward.
+    // When parallel algorithm for bands are adopted. Density will only be treated in the first group.
+    //(Different ranks should have abtained the same, but small differences always exist in practice.)
     // Maybe in the future, density and wavefunctions should use different
     // parallel algorithms, in which they do not occupy all processors, for
     // example wavefunctions uses 20 processors while density uses 10.
     if (PARAM.globalv.ks_run)
     {
-        // double drho = this->estate.caldr2();
-        // EState should be used after it is constructed.
-
         drho = p_chgmix->get_drho(&this->chr, PARAM.inp.nelec);
         hsolver_error = 0.0;
         if (iter == 1 && PARAM.inp.calculation != "nscf")
@@ -140,23 +114,16 @@ void ESolver_KS<T, Device>::hamilt2rho(UnitCell& ucell, const int istep, const i
             // so a more precise HSolver should be executed.
             if (hsolver_error > drho)
             {
-                diag_ethr = hsolver::reset_diag_ethr(GlobalV::ofs_running,
-                                                     PARAM.inp.basis_type,
-                                                     PARAM.inp.esolver_type,
-                                                     PARAM.inp.precision,
-                                                     hsolver_error,
-                                                     drho,
-                                                     diag_ethr,
-                                                     PARAM.inp.nelec);
+                diag_ethr = hsolver::reset_diag_ethr(GlobalV::ofs_running, PARAM.inp.basis_type,
+                            PARAM.inp.esolver_type, PARAM.inp.precision, hsolver_error,
+                            drho, diag_ethr, PARAM.inp.nelec);
 
                 this->hamilt2rho_single(ucell, istep, iter, diag_ethr);
 
                 drho = p_chgmix->get_drho(&this->chr, PARAM.inp.nelec);
 
                 hsolver_error = hsolver::cal_hsolve_error(PARAM.inp.basis_type,
-                                                          PARAM.inp.esolver_type,
-                                                          diag_ethr,
-                                                          PARAM.inp.nelec);
+                                PARAM.inp.esolver_type, diag_ethr, PARAM.inp.nelec);
             }
         }
     }
@@ -168,22 +135,17 @@ void ESolver_KS<T, Device>::runner(UnitCell& ucell, const int istep)
     ModuleBase::TITLE("ESolver_KS", "runner");
     ModuleBase::timer::tick(this->classname, "runner");
 
-    //----------------------------------------------------------------
     // 1) before_scf (electronic iteration loops)
-    //----------------------------------------------------------------
     this->before_scf(ucell, istep);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF");
 
-    //----------------------------------------------------------------
     // 2) SCF iterations
-    //----------------------------------------------------------------
     bool conv_esolver = false;
     this->niter = this->maxniter;
     this->diag_ethr = PARAM.inp.pw_diag_thr;
     this->scf_nmax_flag = false; // mohan add 2025-09-21
     for (int iter = 1; iter <= this->maxniter; ++iter)
 	{
-        // mohan add 2025-09-21
 		if(iter == this->maxniter)
 		{
 			this->scf_nmax_flag=true;
diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp
index 62f00adcb6..3a46d3b505 100644
--- a/source/source_esolver/esolver_ks_lcao.cpp
+++ b/source/source_esolver/esolver_ks_lcao.cpp
@@ -21,13 +21,6 @@
 #include "source_lcao/rho_tau_lcao.h" // mohan add 20251024
 #include "source_lcao/LCAO_set.h" // mohan add 20251111
 
-
-// tmp
-#include "source_psi/setup_psi.h" // use Setup_Psi
-#include "source_io/read_wfc_nao.h" // use read_wfc_nao
-#include "source_estate/elecstate_tools.h" // use fixed_weights
-
-
 namespace ModuleESolver
 {
 
@@ -61,7 +54,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     {
         // TK stands for double and std::complex<double>?
         this->pelec = new elecstate::ElecStateLCAO<TK>(&(this->chr), &(this->kv),
-          this->kv.get_nks(), this->pw_rho, this->pw_big);
+          this->kv.get_nks(), this->pw_big);
     }
 
     // 3) read LCAO orbitals/projectors and construct the interpolation tables.
@@ -458,8 +451,13 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
             }
             this->dftu.cal_energy_correction(ucell, istep);
         }
-        this->dftu.output(ucell);
-    }
+		this->dftu.output(ucell);
+		// use the converged occupation matrix for next MD/Relax SCF calculation
+		if (conv_esolver)
+		{
+			this->dftu.initialed_locale = true;
+		}
+	}
 
     // 2) for deepks, calculate delta_e, output labels during electronic steps
     this->deepks.delta_e(ucell, this->kv, this->orb_, this->pv, this->gd, dm_vec, this->pelec->f_en, PARAM.inp);
@@ -487,12 +485,6 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
         }
     }
 
-    // use the converged occupation matrix for next MD/Relax SCF calculation
-    if (PARAM.inp.dft_plus_u && conv_esolver)
-    {
-        this->dftu.initialed_locale = true;
-    }
-
     // control the output related to the finished iteration
     ModuleIO::ctrl_iter_lcao<TK, TR>(ucell, PARAM.inp, this->kv, this->pelec, *this->dmat.dm,
       this->pv, this->gd, this->psi, this->chr, this->p_chgmix, 
diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp
index 6c19a6c5a1..61ee2e10bc 100644
--- a/source/source_esolver/esolver_ks_pw.cpp
+++ b/source/source_esolver/esolver_ks_pw.cpp
@@ -80,7 +80,7 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
     //! Call before_all_runners() of ESolver_KS
     ESolver_KS<T, Device>::before_all_runners(ucell, inp);
 
-    //! setup and allocation for pelec, charge density, potentials, etc. 
+    //! setup and allocation for pelec, potentials, etc. 
     elecstate::setup_estate_pw<T, Device>(ucell, this->kv, this->sf, this->pelec, this->chr,
       this->locpp, this->ppcell, this->vsep_cell, this->pw_wfc, this->pw_rho,
       this->pw_rhod, this->pw_big, this->solvent, inp);
@@ -286,7 +286,7 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
     // pp projectors, liuyu 2023-10-24
     if (PARAM.globalv.use_uspp)
     {
-        ModuleBase::matrix veff = this->pelec->pot->get_effective_v();
+        ModuleBase::matrix veff = this->pelec->pot->get_eff_v();
         this->ppcell.cal_effective_D(veff, this->pw_rhod, ucell);
     }
 
diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp
index f609b3e2c1..fef5e5fb14 100644
--- a/source/source_esolver/esolver_of.cpp
+++ b/source/source_esolver/esolver_of.cpp
@@ -73,40 +73,15 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp)
 
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
 
-    XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
+//  XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
     int func_type = XC_Functional::get_func_type();
     if (func_type > 2)
     {
         ModuleBase::WARNING_QUIT("esolver_of", "meta-GGA and Hybrid functionals are not supported by OFDFT.");
     }
 
-    // symmetry analysis should be performed every time the cell is changed
-    if (ModuleSymmetry::Symmetry::symm_flag == 1)
-    {
-        ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
-        ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
-    }
-
-    // Setup the k points according to symmetry.
-    kv.set(ucell,ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
-    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
-
-    // print information
-    // mohan add 2021-01-30
-    ModuleIO::print_parameters(ucell, kv, inp);
-
-    // initialize the real-space uniform grid for FFT and parallel
-    // distribution of plane waves
-    Pgrid.init(pw_rho->nx,
-                        pw_rho->ny,
-                        pw_rho->nz,
-                        pw_rho->nplane,
-                        pw_rho->nrxx,
-                        pw_big->nbz,
-                        pw_big->bz); // mohan add 2010-07-22, update 2011-05-04
-    // Calculate Structure factor
-    sf.setup(&ucell, Pgrid, pw_rho);
-    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS");
+    this->chr.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv);
+    this->chr.check_rho(); // check the rho
 
     // initialize local pseudopotential
     this->locpp.init_vloc(ucell,pw_rho);
@@ -122,7 +97,7 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp)
     // liuyu move here 2023-10-09
     // D in uspp need vloc, thus behind init_scf()
     // calculate the effective coefficient matrix for non-local pseudopotential projectors
-    ModuleBase::matrix veff = this->pelec->pot->get_effective_v();
+    ModuleBase::matrix veff = this->pelec->pot->get_eff_v();
 
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT POTENTIAL");
 
@@ -214,6 +189,8 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
     //! 1) call before_scf() of ESolver_FP
     ESolver_FP::before_scf(ucell, istep);
 
+
+
     if (ucell.cell_parameter_updated)
     {
         this->dV_ = ucell.omega / this->pw_rho->nxyz;
@@ -317,10 +294,10 @@ void ESolver_OF::update_potential(UnitCell& ucell)
     this->kedf_manager_->get_potential(this->chr.rho,
                                        this->pphi_,
                                        this->pw_rho,
-                                       this->pelec->pot->get_effective_v()); // KEDF potential
+                                       this->pelec->pot->get_eff_v()); // KEDF potential
     for (int is = 0; is < PARAM.inp.nspin; ++is)
     {
-        const double* vr_eff = this->pelec->pot->get_effective_v(is);
+        const double* vr_eff = this->pelec->pot->get_eff_v(is);
         for (int ir = 0; ir < this->pw_rho->nrxx; ++ir)
         {
             this->pdEdphi_[is][ir] = vr_eff[ir];
@@ -516,9 +493,9 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell, const bool conv_eso
     this->kedf_manager_->get_potential(this->chr.rho,
                                        this->pphi_,
                                        this->pw_rho,
-                                       this->pelec->pot->get_effective_v()); // KEDF potential
+                                       this->pelec->pot->get_eff_v()); // KEDF potential
         
-        const double* vr_eff = this->pelec->pot->get_effective_v(0);
+        const double* vr_eff = this->pelec->pot->get_eff_v(0);
         for (int ir = 0; ir < this->pw_rho->nrxx; ++ir)
         {
             this->pdEdphi_[0][ir] = vr_eff[ir];
diff --git a/source/source_esolver/esolver_of_tool.cpp b/source/source_esolver/esolver_of_tool.cpp
index 6e7ebbd7fb..62aa9aa0b4 100644
--- a/source/source_esolver/esolver_of_tool.cpp
+++ b/source/source_esolver/esolver_of_tool.cpp
@@ -20,7 +20,6 @@ void ESolver_OF::init_elecstate(UnitCell& ucell)
     if (this->pelec == nullptr)
     {
         this->pelec = new elecstate::ElecState((Charge*)(&chr), this->pw_rho, pw_big);
-        this->chr.allocate(PARAM.inp.nspin);
     }
 
     delete this->pelec->pot;
@@ -139,7 +138,7 @@ void ESolver_OF::cal_potential(double* ptemp_phi, double* rdLdphi, UnitCell& uce
 
     elecstate::cal_ux(ucell);
     this->pelec->pot->update_from_charge(this->ptemp_rho_, &ucell);
-    ModuleBase::matrix& vr_eff = this->pelec->pot->get_effective_v();
+    ModuleBase::matrix& vr_eff = this->pelec->pot->get_eff_v();
 
     this->kedf_manager_->get_potential(this->ptemp_rho_->rho,
                                        temp_phi,
@@ -179,7 +178,7 @@ void ESolver_OF::cal_dEdtheta(double** ptemp_phi, Charge* temp_rho, UnitCell& uc
 
     elecstate::cal_ux(ucell);
     this->pelec->pot->update_from_charge(temp_rho, &ucell);
-    ModuleBase::matrix& vr_eff = this->pelec->pot->get_effective_v();
+    ModuleBase::matrix& vr_eff = this->pelec->pot->get_eff_v();
 
     this->kedf_manager_->get_potential(temp_rho->rho,
                                        ptemp_phi,
diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt
index bd7be0e782..9789df523c 100644
--- a/source/source_estate/CMakeLists.txt
+++ b/source/source_estate/CMakeLists.txt
@@ -13,7 +13,6 @@ list(APPEND objects
     module_pot/H_Hartree_pw.cpp
     module_pot/pot_xc.cpp
     module_pot/pot_local.cpp
-    module_pot/pot_local_paw.cpp
     module_pot/potential_new.cpp
     module_pot/potential_types.cpp
     module_pot/pot_sep.cpp
diff --git a/source/source_estate/elecstate.cpp b/source/source_estate/elecstate.cpp
index 58fdf7344a..1831ed8fce 100644
--- a/source/source_estate/elecstate.cpp
+++ b/source/source_estate/elecstate.cpp
@@ -40,11 +40,13 @@ void ElecState::init_scf(const int istep,
 
     //! other effective potentials need charge density,
     // choose charge density from ionic step 0.
+/*
     if (istep == 0)
     {
         this->charge->init_rho(this->eferm,ucell, pgrid, strucfac, symm, (const void*)this->klist, wfcpw);
         this->charge->check_rho(); // check the rho
     }
+*/
 
     //! renormalize the charge density
     this->charge->renormalize_rho();
@@ -57,12 +59,11 @@ void ElecState::init_scf(const int istep,
 void ElecState::init_ks(Charge* chr_in, // pointer for class Charge
                         const K_Vectors* klist_in,
                         int nk_in,
-                        ModulePW::PW_Basis* rhopw_in,
                         const ModulePW::PW_Basis_Big* bigpw_in)
 {
     this->charge = chr_in;
     this->klist = klist_in;
-    this->charge->set_rhopw(rhopw_in);
+//    this->charge->set_rhopw(rhopw_in); // mohan comment out 20251130
     this->bigpw = bigpw_in;
     // init nelec_spin with nelec and nupdown
     this->init_nelec_spin();
diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h
index 4ac3cb8f8a..1ed769a875 100644
--- a/source/source_estate/elecstate.h
+++ b/source/source_estate/elecstate.h
@@ -35,7 +35,6 @@ class ElecState
     void init_ks(Charge* chr_in, // pointer for class Charge
                  const K_Vectors* klist_in,
                  int nk_in, // number of k points
-                 ModulePW::PW_Basis* rhopw_in,
                  const ModulePW::PW_Basis_Big* bigpw_in);
 
     // return current electronic density rho, as a input for constructing Hamiltonian
diff --git a/source/source_estate/elecstate_energy.cpp b/source/source_estate/elecstate_energy.cpp
index b6079df0aa..bbdc269fec 100644
--- a/source/source_estate/elecstate_energy.cpp
+++ b/source/source_estate/elecstate_energy.cpp
@@ -105,15 +105,15 @@ void ElecState::cal_bandgap_updw()
 /// @brief calculate deband
 double ElecState::cal_delta_eband(const UnitCell& ucell) const
 {
-    // out potentials from potential mixing
-    // total energy and band energy corrections
+	ModuleBase::timer::tick("ElecState", "cal_delta_eband");
+	// out potentials from potential mixing
+	// total energy and band energy corrections
     double deband0 = 0.0;
-
     double deband_aux = 0.0;
 
     // only potential related with charge is used here for energy correction
-    // on the fly calculate it here by v_effective - v_fixed
-    const double* v_eff = this->pot->get_effective_v(0);
+    // on the fly calculate it here by v_eff - v_fixed
+    const double* v_eff = this->pot->get_eff_v(0);
     const double* v_fixed = this->pot->get_fixed_v();
     const double* v_ofk = nullptr;
     const bool v_ofk_flag = (XC_Functional::get_ked_flag());
@@ -122,10 +122,11 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
     {
         deband_aux -= this->charge->rho[0][ir] * (v_eff[ir] - v_fixed[ir]);
     }
+
     if (v_ofk_flag)
     {
-        v_ofk = this->pot->get_effective_vofk(0);
-        // cause in the get_effective_vofk, the func will return nullptr
+        v_ofk = this->pot->get_eff_vofk(0);
+        // cause in the get_eff_vofk, the func will return nullptr
         if (v_ofk == nullptr && this->charge->rhopw->nrxx > 0)
         {
             ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband", "v_ofk is nullptr");
@@ -138,14 +139,14 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
 
     if (PARAM.inp.nspin == 2)
     {
-        v_eff = this->pot->get_effective_v(1);
+        v_eff = this->pot->get_eff_v(1);
         for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
         {
             deband_aux -= this->charge->rho[1][ir] * (v_eff[ir] - v_fixed[ir]);
         }
         if (v_ofk_flag)
         {
-            v_ofk = this->pot->get_effective_vofk(1);
+            v_ofk = this->pot->get_eff_vofk(1);
             if (v_ofk == nullptr && this->charge->rhopw->nrxx > 0)
             {
                 ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband", "v_ofk is nullptr");
@@ -160,7 +161,7 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
     {
         for (int is = 1; is < 4; is++)
         {
-            v_eff = this->pot->get_effective_v(is);
+            v_eff = this->pot->get_eff_v(is);
             for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
             {
                 deband_aux -= this->charge->rho[is][ir] * v_eff[ir];
@@ -178,13 +179,16 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
 
     // \int rho(r) v_{exx}(r) dr = 2 E_{exx}[rho]
     deband0 -= 2 * this->f_en.exx; // Peize Lin add 2017-10-16
+
+	ModuleBase::timer::tick("ElecState", "cal_delta_eband");
     return deband0;
 }
 
 /// @brief calculate descf
 double ElecState::cal_delta_escf() const
 {
-    ModuleBase::TITLE("energy", "delta_escf");
+    ModuleBase::TITLE("ElecState", "cal_delta_escf");
+	ModuleBase::timer::tick("ElecState", "cal_delta_escf");
     double descf = 0.0;
 
     // now rho1 is "mixed" charge density
@@ -192,21 +196,21 @@ double ElecState::cal_delta_escf() const
     // because in "deband" the energy is calculated from "output" charge density,
     // so here is the correction.
     // only potential related with charge is used here for energy correction
-    // on the fly calculate it here by v_effective - v_fixed
-    const double* v_eff = this->pot->get_effective_v(0);
+    // on the fly calculate it here by v_eff - v_fixed
+    const double* v_eff = this->pot->get_eff_v(0);
     const double* v_fixed = this->pot->get_fixed_v();
     const double* v_ofk = nullptr;
 
     if (XC_Functional::get_ked_flag())
     {
-        v_ofk = this->pot->get_effective_vofk(0);
+        v_ofk = this->pot->get_eff_vofk(0);
     }
     for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
     {
         descf -= (this->charge->rho[0][ir] - this->charge->rho_save[0][ir]) * (v_eff[ir] - v_fixed[ir]);
         if (XC_Functional::get_ked_flag())
         {
-            // cause in the get_effective_vofk, the func will return nullptr
+            // cause in the get_eff_vofk, the func will return nullptr
             assert(v_ofk != nullptr);
             descf -= (this->charge->kin_r[0][ir] - this->charge->kin_r_save[0][ir]) * v_ofk[ir];
         }
@@ -214,10 +218,10 @@ double ElecState::cal_delta_escf() const
 
     if (PARAM.inp.nspin == 2)
     {
-        v_eff = this->pot->get_effective_v(1);
+        v_eff = this->pot->get_eff_v(1);
         if (XC_Functional::get_ked_flag())
         {
-            v_ofk = this->pot->get_effective_vofk(1);
+            v_ofk = this->pot->get_eff_vofk(1);
         }
         for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
         {
@@ -232,7 +236,7 @@ double ElecState::cal_delta_escf() const
     {
         for (int is = 1; is < 4; is++)
         {
-            v_eff = this->pot->get_effective_v(is);
+            v_eff = this->pot->get_eff_v(is);
             for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
             {
                 descf -= (this->charge->rho[is][ir] - this->charge->rho_save[is][ir]) * v_eff[ir];
@@ -247,6 +251,14 @@ double ElecState::cal_delta_escf() const
     assert(this->charge->rhopw->nxyz > 0);
 
     descf *= this->charge->rhopw->omega / this->charge->rhopw->nxyz;
+
+// mohan move the code here, 2025-11-28
+#ifdef __MPI
+        MPI_Bcast(&descf, 1, MPI_DOUBLE, 0, BP_WORLD);
+#endif
+
+
+	ModuleBase::timer::tick("ElecState", "cal_delta_escf");
     return descf;
 }
 
@@ -311,7 +323,7 @@ void ElecState::cal_energies(const int type)
     }
     else
     {
-        ModuleBase::WARNING_QUIT("cal_energies", "The form of total energy functional is unknown!");
+        ModuleBase::WARNING_QUIT("ElecState::cal_energies", "The form of total energy functional is unknown!");
     }
 }
 
diff --git a/source/source_estate/elecstate_lcao.h b/source/source_estate/elecstate_lcao.h
index 3cf06641e3..bf1f11e1f7 100644
--- a/source/source_estate/elecstate_lcao.h
+++ b/source/source_estate/elecstate_lcao.h
@@ -18,10 +18,9 @@ class ElecStateLCAO : public ElecState
     ElecStateLCAO(Charge* chr_in,
                   const K_Vectors* klist_in,
                   int nks_in,
-                  ModulePW::PW_Basis* rhopw_in,
                   ModulePW::PW_Basis_Big* bigpw_in)
     {
-        init_ks(chr_in, klist_in, nks_in, rhopw_in, bigpw_in);
+        init_ks(chr_in, klist_in, nks_in, bigpw_in);
         this->classname = "ElecStateLCAO";
     }
 
diff --git a/source/source_estate/elecstate_pw.cpp b/source/source_estate/elecstate_pw.cpp
index 247227149e..2d24321f2a 100644
--- a/source/source_estate/elecstate_pw.cpp
+++ b/source/source_estate/elecstate_pw.cpp
@@ -17,7 +17,6 @@ ElecStatePW<T, Device>::ElecStatePW(ModulePW::PW_Basis_K* wfc_basis_in,
                                     K_Vectors* pkv_in,
                                     UnitCell* ucell_in,
                                     pseudopot_cell_vnl* ppcell_in,
-                                    ModulePW::PW_Basis* rhodpw_in,
                                     ModulePW::PW_Basis* rhopw_in,
                                     ModulePW::PW_Basis_Big* bigpw_in)
     : basis(wfc_basis_in)
@@ -26,7 +25,7 @@ ElecStatePW<T, Device>::ElecStatePW(ModulePW::PW_Basis_K* wfc_basis_in,
     this->rhopw_smooth = rhopw_in;
     this->ppcell = ppcell_in;
     this->ucell = ucell_in;
-    this->init_ks(chr_in, pkv_in, pkv_in->get_nks(), rhodpw_in, bigpw_in);
+    this->init_ks(chr_in, pkv_in, pkv_in->get_nks(), bigpw_in);
 }
 
 template<typename T, typename Device>
diff --git a/source/source_estate/elecstate_pw.h b/source/source_estate/elecstate_pw.h
index 10b54676e6..f674beed89 100644
--- a/source/source_estate/elecstate_pw.h
+++ b/source/source_estate/elecstate_pw.h
@@ -24,7 +24,6 @@ class ElecStatePW : public ElecState
                 K_Vectors* pkv_in,
                 UnitCell* ucell_in,
                 pseudopot_cell_vnl* ppcell_in,
-                ModulePW::PW_Basis* rhodpw_in,
                 ModulePW::PW_Basis* rhopw_in,
                 ModulePW::PW_Basis_Big* bigpw_in);
 
diff --git a/source/source_estate/elecstate_pw_sdft.cpp b/source/source_estate/elecstate_pw_sdft.cpp
index 4c04fbd886..1df5545cb5 100644
--- a/source/source_estate/elecstate_pw_sdft.cpp
+++ b/source/source_estate/elecstate_pw_sdft.cpp
@@ -42,4 +42,4 @@ template class ElecStatePW_SDFT<std::complex<double>, base_device::DEVICE_CPU>;
 #if ((defined __CUDA) || (defined __ROCM))
 template class ElecStatePW_SDFT<std::complex<double>, base_device::DEVICE_GPU>;
 #endif
-} // namespace elecstate
\ No newline at end of file
+} // namespace elecstate
diff --git a/source/source_estate/elecstate_pw_sdft.h b/source/source_estate/elecstate_pw_sdft.h
index 506cabdd89..0ffa00efc4 100644
--- a/source/source_estate/elecstate_pw_sdft.h
+++ b/source/source_estate/elecstate_pw_sdft.h
@@ -12,11 +12,10 @@ class ElecStatePW_SDFT : public ElecStatePW<T, Device>
                      K_Vectors* pkv_in,
                      UnitCell* ucell_in,
                      pseudopot_cell_vnl* ppcell_in,
-                     ModulePW::PW_Basis* rhodpw_in,
                      ModulePW::PW_Basis* rhopw_in,
                      ModulePW::PW_Basis_Big* bigpw_in)
         : ElecStatePW<T,
-                      Device>(wfc_basis_in, chr_in, pkv_in, ucell_in, ppcell_in, rhodpw_in, rhopw_in, bigpw_in)
+                      Device>(wfc_basis_in, chr_in, pkv_in, ucell_in, ppcell_in, rhopw_in, bigpw_in)
     {
         this->classname = "ElecStatePW_SDFT";
     }
diff --git a/source/source_estate/module_charge/charge.cpp b/source/source_estate/module_charge/charge.cpp
index be86a04800..9dfba3782a 100644
--- a/source/source_estate/module_charge/charge.cpp
+++ b/source/source_estate/module_charge/charge.cpp
@@ -80,10 +80,17 @@ void Charge::destroy()
 void Charge::allocate(const int& nspin_in)
 {
     ModuleBase::TITLE("Charge", "allocate");
+
+	if (this->rhopw == nullptr)
+	{
+		ModuleBase::WARNING_QUIT("Charge::allocate","rhopw is nullptr.");
+	}
+
     this->nrxx = this->rhopw->nrxx;
     this->nxyz = this->rhopw->nxyz;
     this->ngmc = this->rhopw->npw;
 
+
     if (allocate_rho == true)
     {
         this->destroy();
@@ -634,6 +641,8 @@ void Charge::save_rho_before_sum_band()
 
 double Charge::cal_rho2ne(const double* rho_in) const
 {
+    assert(this->rhopw->nxyz > 0); // mohan add 2025-12-02
+
     double ne = 0.0;
     for (int ir = 0; ir < this->rhopw->nrxx; ir++)
     {
diff --git a/source/source_estate/module_charge/charge.h b/source/source_estate/module_charge/charge.h
index 479eb8dc66..a9190676a0 100644
--- a/source/source_estate/module_charge/charge.h
+++ b/source/source_estate/module_charge/charge.h
@@ -7,7 +7,7 @@
 #include "source_base/parallel_global.h"
 #include "source_basis/module_pw/pw_basis.h"
 #include "source_cell/module_symmetry/symmetry.h"
-#include "source_estate/fp_energy.h"
+// #include "source_estate/fp_energy.h"
 #include "source_pw/module_pwdft/parallel_grid.h"
 
 //a forward declaration of UnitCell
@@ -77,8 +77,7 @@ class Charge
      * @param klist [in] k points list if needed
      * @param wfcpw [in] PW basis for wave function if needed
      */
-    void init_rho(elecstate::Efermi& eferm_iout,
-                  const UnitCell& ucell,
+    void init_rho(const UnitCell& ucell,
                   const Parallel_Grid& pgrid,
                   const ModuleBase::ComplexMatrix& strucFac,
                   ModuleSymmetry::Symmetry& symm,
@@ -97,8 +96,6 @@ class Charge
                       const ModuleBase::ComplexMatrix& structure_factor, 
                       const bool* numeric);
 
-    void set_rho_core_paw();
-
     void renormalize_rho();
 
     double sum_rho() const;
diff --git a/source/source_estate/module_charge/charge_init.cpp b/source/source_estate/module_charge/charge_init.cpp
index b3a034c3d2..95f675eb90 100644
--- a/source/source_estate/module_charge/charge_init.cpp
+++ b/source/source_estate/module_charge/charge_init.cpp
@@ -18,8 +18,7 @@
 #include "source_io/rhog_io.h"
 #include "source_io/read_wf2rho_pw.h"
 
-void Charge::init_rho(elecstate::Efermi& eferm_iout,
-                      const UnitCell& ucell,
+void Charge::init_rho(const UnitCell& ucell,
                       const Parallel_Grid& pgrid,
                       const ModuleBase::ComplexMatrix& strucFac,
                       ModuleSymmetry::Symmetry& symm,
@@ -372,11 +371,6 @@ void Charge::set_rho_core(const UnitCell& ucell,
     return;
 } // end subroutine set_rhoc
 
-void Charge::set_rho_core_paw()
-{
-    ModuleBase::TITLE("Charge","set_rho_core_paw");
-}
-
 
 void Charge::non_linear_core_correction
 (
@@ -397,7 +391,7 @@ void Charge::non_linear_core_correction
 
 	double gx = 0.0;
     double rhocg1 = 0.0;
-    double *aux;
+    double *aux = nullptr;
 
     // here we compute the fourier transform is the charge in numeric form
     if (numeric)
diff --git a/source/source_estate/module_dm/cal_dm_psi.cpp b/source/source_estate/module_dm/cal_dm_psi.cpp
index 7f68838c94..bd14f21774 100644
--- a/source/source_estate/module_dm/cal_dm_psi.cpp
+++ b/source/source_estate/module_dm/cal_dm_psi.cpp
@@ -15,8 +15,8 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,
                        const psi::Psi<double>& wfc,
                        elecstate::DensityMatrix<double, double>& DM)
 {
-    ModuleBase::TITLE("elecstate", "cal_dm");
-    ModuleBase::timer::tick("elecstate", "cal_dm");
+    ModuleBase::TITLE("elecstate", "cal_dm_psi");
+    ModuleBase::timer::tick("elecstate", "cal_dm_psi");
 
     // dm.resize(wfc.get_nk(), ParaV->ncol, ParaV->nrow);
     const int nbands_local = wfc.get_nbands();
@@ -33,13 +33,10 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,
         // dm[ik].create(ParaV->ncol, ParaV->nrow);
         //  wg_wfc(ib,iw) = wg[ib] * wfc(ib,iw);
 
-        psi::Psi<double> wg_wfc(1, 
-                                wfc.get_nbands(), 
-                                wfc.get_nbasis(),
-                                wfc.get_nbasis(),
-                                true);
-        wg_wfc.set_all_psi(wfc.get_pointer(), wg_wfc.size());
+        psi::Psi<double> wg_wfc(1, wfc.get_nbands(), 
+          wfc.get_nbasis(), wfc.get_nbasis(), true);
 
+        wg_wfc.set_all_psi(wfc.get_pointer(), wg_wfc.size());
 
         int ib_global = 0;
         for (int ib_local = 0; ib_local < nbands_local; ++ib_local)
@@ -50,7 +47,6 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,
                 if (ib_global >= wg.nc)
                 {
                     break;
-                    ModuleBase::WARNING_QUIT("ElecStateLCAO::cal_dm", "please check global2local_col!");
                 }
             }
 			if (ib_global >= wg.nc)
@@ -70,7 +66,7 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,
         psiMulPsi(wg_wfc, wfc, dmk_pointer);
 #endif
     }
-    ModuleBase::timer::tick("elecstate", "cal_dm");
+    ModuleBase::timer::tick("elecstate", "cal_dm_psi");
 
     return;
 }
@@ -80,8 +76,8 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,
                        const psi::Psi<std::complex<double>>& wfc,
                        elecstate::DensityMatrix<std::complex<double>, double>& DM)
 {
-    ModuleBase::TITLE("elecstate", "cal_dm");
-    ModuleBase::timer::tick("elecstate", "cal_dm");
+    ModuleBase::TITLE("elecstate", "cal_dm_psi");
+    ModuleBase::timer::tick("elecstate", "cal_dm_psi");
 
     // dm.resize(wfc.get_nk(), ParaV->ncol, ParaV->nrow);
     const int nbands_local = wfc.get_nbands();
@@ -104,6 +100,7 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,
         
         const std::complex<double>* pwfc = wfc.get_pointer();
         std::complex<double>* pwg_wfc = wg_wfc.get_pointer();
+
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 1024)
 #endif
@@ -139,8 +136,9 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,
         if (PARAM.inp.ks_solver == "cg_in_lcao")
         {
             psiMulPsi(wg_wfc, wfc, dmk_pointer);
-        } else 
-        {
+		} 
+		else 
+		{
             psiMulPsiMpi(wg_wfc, wfc, dmk_pointer, ParaV->desc_wfc, ParaV->desc);
         }
 #else
@@ -148,7 +146,7 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,
 #endif
     }
 
-    ModuleBase::timer::tick("elecstate", "cal_dm");
+    ModuleBase::timer::tick("elecstate", "cal_dm_psi");
     return;
 }
 
diff --git a/source/source_estate/module_dm/density_matrix.cpp b/source/source_estate/module_dm/density_matrix.cpp
index 8dd87e73c4..7282226a84 100644
--- a/source/source_estate/module_dm/density_matrix.cpp
+++ b/source/source_estate/module_dm/density_matrix.cpp
@@ -29,7 +29,7 @@ template <typename TK, typename TR>
 DensityMatrix<TK, TR>::DensityMatrix(const Parallel_Orbitals* paraV_in, const int nspin, const std::vector<ModuleBase::Vector3<double>>& kvec_d, const int nk)
     : _paraV(paraV_in), _nspin(nspin), _kvec_d(kvec_d), _nk((nk > 0 && nk <= _kvec_d.size()) ? nk : _kvec_d.size())
 {
-    ModuleBase::TITLE("DensityMatrix", "DensityMatrix-MK");
+    ModuleBase::TITLE("DensityMatrix", "resize_DMK");
     const int nks = _nk * _nspin;
     this->_DMK.resize(nks);
     for (int ik = 0; ik < nks; ik++)
@@ -42,7 +42,7 @@ DensityMatrix<TK, TR>::DensityMatrix(const Parallel_Orbitals* paraV_in, const in
 template <typename TK, typename TR>
 DensityMatrix<TK, TR>::DensityMatrix(const Parallel_Orbitals* paraV_in, const int nspin) :_paraV(paraV_in), _nspin(nspin), _kvec_d({ ModuleBase::Vector3<double>(0,0,0) }), _nk(1)
 {
-    ModuleBase::TITLE("DensityMatrix", "DensityMatrix-GO");
+    ModuleBase::TITLE("DensityMatrix", "resize_gamma");
     this->_DMK.resize(_nspin);
     for (int ik = 0; ik < this->_nspin; ik++)
     {
diff --git a/source/source_estate/module_dm/init_dm.cpp b/source/source_estate/module_dm/init_dm.cpp
index 669360ee1e..c45190e01e 100644
--- a/source/source_estate/module_dm/init_dm.cpp
+++ b/source/source_estate/module_dm/init_dm.cpp
@@ -2,6 +2,7 @@
 #include "source_estate/module_dm/cal_dm_psi.h"
 #include "source_estate/elecstate_tools.h"
 #include "source_estate/cal_ux.h"
+#include "source_lcao/rho_tau_lcao.h" // mohan add 2025-11-12
 
 template <typename TK>
 void elecstate::init_dm(UnitCell& ucell,
@@ -16,26 +17,15 @@ void elecstate::init_dm(UnitCell& ucell,
 
 	if (iter == 1 && exx_two_level_step == 0)
 	{
-		std::cout << " WAVEFUN -> CHARGE " << std::endl;
-
-		// calculate the density matrix using read in wave functions
-		// and then calculate the charge density on grid.
-
-		pelec->skip_weights = true;
-		elecstate::calculate_weights(pelec->ekb,
-				pelec->wg,
-				pelec->klist,
-				pelec->eferm,
-				pelec->f_en,
-				pelec->nelec_spin,
-				pelec->skip_weights);
+		std::cout << " LCAO WAVEFUN -> CHARGE " << std::endl;
 
 		elecstate::calEBand(pelec->ekb, pelec->wg, pelec->f_en);
+
 		elecstate::cal_dm_psi(dmat.dm->get_paraV_pointer(), pelec->wg, *psi, *dmat.dm);
 		dmat.dm->cal_DMR();
 
-//		pelec->psiToRho(*psi); // I found this sentence is useless, mohan add 2025-11-04
-		pelec->skip_weights = false;
+        // mohan add 2025-11-12, use density matrix to calculate the charge density
+        LCAO_domain::dm2rho(dmat.dm->get_DMR_vector(), PARAM.inp.nspin, &chr);
 
 		elecstate::cal_ux(ucell);
 
diff --git a/source/source_estate/module_pot/pot_local_paw.cpp b/source/source_estate/module_pot/pot_local_paw.cpp
deleted file mode 100644
index 03dab7d59e..0000000000
--- a/source/source_estate/module_pot/pot_local_paw.cpp
+++ /dev/null
@@ -1,25 +0,0 @@
-#include "pot_local_paw.h"
-
-#include "source_base/timer.h"
-#include "source_base/tool_title.h"
-
-#include <complex>
-
-namespace elecstate
-{
-
-//==========================================================
-// This routine computes the local potential in real space
-//==========================================================
-void PotLocal_PAW::cal_fixed_v(double *vl_pseudo // store the local pseudopotential
-)
-{
-    ModuleBase::TITLE("PotLocal_PAW", "cal_fixed_v");
-    ModuleBase::timer::tick("PotLocal_PAW", "cal_fixed_v");
-
-    // GlobalV::ofs_running <<" set local pseudopotential done." << std::endl;
-    ModuleBase::timer::tick("PotLocal_PAW", "cal_fixed_v");
-    return;
-}
-
-} // namespace elecstate
\ No newline at end of file
diff --git a/source/source_estate/module_pot/pot_local_paw.h b/source/source_estate/module_pot/pot_local_paw.h
deleted file mode 100644
index 7399827c33..0000000000
--- a/source/source_estate/module_pot/pot_local_paw.h
+++ /dev/null
@@ -1,24 +0,0 @@
-#ifndef POTLOCALPAW_H
-#define POTLOCALPAW_H
-
-#include "source_base/matrix.h"
-#include "pot_base.h"
-
-namespace elecstate
-{
-
-class PotLocal_PAW : public PotBase
-{
-  public:
-    PotLocal_PAW()
-    {
-      this->fixed_mode = true;
-      this->dynamic_mode = false;
-    }
-
-    void cal_fixed_v(double* vl_pseudo) override;
-};
-
-} // namespace elecstate
-
-#endif
\ No newline at end of file
diff --git a/source/source_estate/module_pot/potential_new.cpp b/source/source_estate/module_pot/potential_new.cpp
index b1968339aa..ba59b5e8ed 100644
--- a/source/source_estate/module_pot/potential_new.cpp
+++ b/source/source_estate/module_pot/potential_new.cpp
@@ -109,10 +109,10 @@ void Potential::allocate()
 		return;
 	}
 
-    this->v_effective_fixed.resize(nrxx);
+    this->v_eff_fixed.resize(nrxx);
     ModuleBase::Memory::record("Pot::veff_fix", sizeof(double) * nrxx);
 
-    this->v_effective.create(nspin, nrxx);
+    this->v_eff.create(nspin, nrxx);
     ModuleBase::Memory::record("Pot::veff", sizeof(double) * nspin * nrxx);
 
     this->veff_smooth.create(nspin, nrxx_smooth);
@@ -120,7 +120,7 @@ void Potential::allocate()
 
     if (XC_Functional::get_ked_flag())
     {
-        this->vofk_effective.create(nspin, nrxx);
+        this->vofk_eff.create(nspin, nrxx);
         ModuleBase::Memory::record("Pot::vofk", sizeof(double) * nspin * nrxx);
 
         this->vofk_smooth.create(nspin, nrxx_smooth);
@@ -162,11 +162,11 @@ void Potential::update_from_charge(const Charge*const chg, const UnitCell*const
 
     if (!this->fixed_done)
     {
-        this->cal_fixed_v(this->v_effective_fixed.data());
+        this->cal_fixed_v(this->v_eff_fixed.data());
         this->fixed_done = true;
     }
 
-    this->cal_v_eff(chg, ucell, this->v_effective);
+    this->cal_v_eff(chg, ucell, this->v_eff);
 
     // interpolate potential on the smooth mesh if necessary
     this->interpolate_vrs();
@@ -202,7 +202,7 @@ void Potential::cal_fixed_v(double* vl_pseudo)
     ModuleBase::TITLE("Potential", "cal_fixed_v");
     ModuleBase::timer::tick("Potential", "cal_fixed_v");
 
-    this->v_effective_fixed.assign(this->v_effective_fixed.size(), 0.0);
+    this->v_eff_fixed.assign(this->v_eff_fixed.size(), 0.0);
     for (size_t i = 0; i < this->components.size(); i++)
     {
         if (this->components[i]->fixed_mode)
@@ -219,10 +219,10 @@ void Potential::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Mo
     ModuleBase::TITLE("Potential", "cal_veff");
     ModuleBase::timer::tick("Potential", "cal_veff");
 
-    const int nspin_current = this->v_effective.nr;
-    const int nrxx = this->v_effective.nc;
-    // first of all, set v_effective to zero.
-    this->v_effective.zero_out();
+    const int nspin_current = this->v_eff.nr;
+    const int nrxx = this->v_eff.nc;
+    // first of all, set v_eff to zero.
+    this->v_eff.zero_out();
 
     // add fixed potential components
     // nspin = 2, add fixed components for all
@@ -231,11 +231,11 @@ void Potential::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Mo
     {
         if (i == 0 || nspin_current == 2)
         {
-            ModuleBase::GlobalFunc::COPYARRAY(this->v_effective_fixed.data(), this->get_effective_v(i), nrxx);
+            ModuleBase::GlobalFunc::COPYARRAY(this->v_eff_fixed.data(), this->get_eff_v(i), nrxx);
         }
     }
 
-    // cal effective by every components
+    // cal eff by every components
     for (size_t i = 0; i < this->components.size(); i++)
     {
         if (this->components[i]->dynamic_mode)
@@ -265,14 +265,14 @@ void Potential::init_pot(int istep, const Charge*const chg)
 void Potential::get_vnew(const Charge* chg, ModuleBase::matrix& vnew)
 {
     ModuleBase::TITLE("Potential", "get_vnew");
-    vnew.create(this->v_effective.nr, this->v_effective.nc);
-    vnew = this->v_effective;
+    vnew.create(this->v_eff.nr, this->v_eff.nc);
+    vnew = this->v_eff;
 
     this->update_from_charge(chg, this->ucell_);
     //(used later for scf correction to the forces )
     for (int iter = 0; iter < vnew.nr * vnew.nc; ++iter)
     {
-        vnew.c[iter] = this->v_effective.c[iter] - vnew.c[iter];
+        vnew.c[iter] = this->v_eff.c[iter] - vnew.c[iter];
     }
 
     return;
@@ -296,7 +296,7 @@ void Potential::interpolate_vrs(void)
         ModuleBase::ComplexMatrix vrs(nspin, rho_basis_->npw);
         for (int is = 0; is < nspin; is++)
         {
-            rho_basis_->real2recip(&v_effective(is, 0), &vrs(is, 0));
+            rho_basis_->real2recip(&v_eff(is, 0), &vrs(is, 0));
             rho_basis_smooth_->recip2real(&vrs(is, 0), &veff_smooth(is, 0));
         }
 
@@ -305,15 +305,15 @@ void Potential::interpolate_vrs(void)
             ModuleBase::ComplexMatrix vrs_ofk(nspin, rho_basis_->npw);
             for (int is = 0; is < nspin; is++)
             {
-                rho_basis_->real2recip(&vofk_effective(is, 0), &vrs_ofk(is, 0));
+                rho_basis_->real2recip(&vofk_eff(is, 0), &vrs_ofk(is, 0));
                 rho_basis_smooth_->recip2real(&vrs_ofk(is, 0), &vofk_smooth(is, 0));
             }
         }
     }
     else
     {
-        this->veff_smooth = this->v_effective;
-        this->vofk_smooth = this->vofk_effective;
+        this->veff_smooth = this->v_eff;
+        this->vofk_smooth = this->vofk_eff;
     }
 
     ModuleBase::timer::tick("Potential", "interpolate_vrs");
diff --git a/source/source_estate/module_pot/potential_new.h b/source/source_estate/module_pot/potential_new.h
index ec2688ed9b..6ed7bdc93d 100644
--- a/source/source_estate/module_pot/potential_new.h
+++ b/source/source_estate/module_pot/potential_new.h
@@ -18,7 +18,7 @@ namespace elecstate
  * 2. Func init_pot()
  *     a. need istep for update_for_tddft();
  *     b. need Charge for update_from_charge();
- *     c. it will reset fixed_done to false, v_effective_fixed will be calculated;
+ *     c. it will reset fixed_done to false, v_eff_fixed will be calculated;
  *     d. it should be called after Charge is initialized;
  *     e. it can only be called once in one SCF loop
  * 3. Func pot_register() and components
@@ -30,8 +30,8 @@ namespace elecstate
  *     f. "efield", PotEfield introduces electronic field including dipole correction part of potentials;
  *     g. "gatefield", PotGate introduces gate field part of potentials;
  * 4. Func update_from_charge()
- *     a. regenerate v_effective
- *     b. if Meta-GGA is choosed, it will regenerate vofk_effective
+ *     a. regenerate v_eff
+ *     b. if Meta-GGA is choosed, it will regenerate vofk_eff
  * 5. Func update_for_tddft()
  *     a. in principle, it should be added to components, but it related to real time(istep)
  *     b. it should be called after update_from_charge() as a compensation;
@@ -41,9 +41,9 @@ namespace elecstate
  *         2. use the final delta_V_eff for calculating force correction
  * 7. Func write_potential()
  * 8. Func write_elecstat_pot()
- * 9. interfaces for v_effective_fixed/v_effective/vofk_effective
+ * 9. interfaces for v_eff_fixed/v_eff/vofk_eff
  * 10. Func interpolate_vrs()
- *    a. interpolate v_effective on the smooth mesh
+ *    a. interpolate v_eff on the smooth mesh
  */
 class Potential : public PotBase
 {
@@ -53,7 +53,7 @@ class Potential : public PotBase
     // In constructor, size of every potential components should be allocated
     // rho_basis_in is the dense grids, rho_basis_smooth_in is the smooth grids in USPP
     // charge density and potential are defined on dense grids,
-    // but effective potential needs to be interpolated on smooth grids in order to compute Veff|psi>
+    // but eff potential needs to be interpolated on smooth grids in order to compute Veff|psi>
     // Note: rho_basis_in and rho_basis_smooth_in are the same in NCPP
     Potential(const ModulePW::PW_Basis* rho_basis_in,
               const ModulePW::PW_Basis* rho_basis_smooth_in,
@@ -78,61 +78,61 @@ class Potential : public PotBase
     PotBase* get_pot_type(const std::string& pot_type);
 
     // interfaces to get values
-    ModuleBase::matrix& get_effective_v()
+    ModuleBase::matrix& get_eff_v()
     {
-        return this->v_effective;
+        return this->v_eff;
     }
-    const ModuleBase::matrix& get_effective_v() const
+    const ModuleBase::matrix& get_eff_v() const
     {
-        return this->v_effective;
+        return this->v_eff;
     }
 
-    double* get_effective_v(int is)
+    double* get_eff_v(int is)
     {
-        if (this->v_effective.nc > 0)
+        if (this->v_eff.nc > 0)
         {
-            return &(this->v_effective(is, 0));
+            return &(this->v_eff(is, 0));
         }
         else
         {
             return nullptr;
         }
     }
-    const double* get_effective_v(int is) const
+    const double* get_eff_v(int is) const
     {
-        if (this->v_effective.nc > 0)
+        if (this->v_eff.nc > 0)
         {
-            return &(this->v_effective(is, 0));
+            return &(this->v_eff(is, 0));
         }
         else
         {
             return nullptr;
         }
     }
-    ModuleBase::matrix& get_effective_vofk()
+    ModuleBase::matrix& get_eff_vofk()
     {
-        return this->vofk_effective;
+        return this->vofk_eff;
     }
-    const ModuleBase::matrix& get_effective_vofk() const
+    const ModuleBase::matrix& get_eff_vofk() const
     {
-        return this->vofk_effective;
+        return this->vofk_eff;
     }
-    double* get_effective_vofk(int is)
+    double* get_eff_vofk(int is)
     {
-        if (this->vofk_effective.nc > 0)
+        if (this->vofk_eff.nc > 0)
         {
-            return &(this->vofk_effective(is, 0));
+            return &(this->vofk_eff(is, 0));
         }
         else
         {
             return nullptr;
         }
     }
-    const double* get_effective_vofk(int is) const
+    const double* get_eff_vofk(int is) const
     {
-        if (this->vofk_effective.nc > 0)
+        if (this->vofk_eff.nc > 0)
         {
-            return &(this->vofk_effective(is, 0));
+            return &(this->vofk_eff(is, 0));
         }
         else
         {
@@ -166,11 +166,11 @@ class Potential : public PotBase
 
     double* get_fixed_v()
     {
-        return this->v_effective_fixed.data();
+        return this->v_eff_fixed.data();
     }
     const double* get_fixed_v() const
     {
-        return this->v_effective_fixed.data();
+        return this->v_eff_fixed.data();
     }
     const ModulePW::PW_Basis *get_rho_basis() const
     {
@@ -195,18 +195,20 @@ class Potential : public PotBase
 
     void allocate();
 
-    std::vector<double> v_effective_fixed;
-    ModuleBase::matrix v_effective;
+    std::vector<double> v_eff_fixed;
+    ModuleBase::matrix v_eff;
 
     ModuleBase::matrix veff_smooth; // used in uspp liuyu 2023-10-12
     ModuleBase::matrix vofk_smooth; // used in uspp liuyu 2023-10-12
 
     ModuleBase::matrix v_xc; // if PAW is used, vxc must be stored separately
 
-    float *s_veff_smooth = nullptr, *s_vofk_smooth = nullptr;
-    double *d_veff_smooth = nullptr, *d_vofk_smooth = nullptr;
+    float *s_veff_smooth = nullptr;
+    float *s_vofk_smooth = nullptr;
+    double *d_veff_smooth = nullptr;
+    double *d_vofk_smooth = nullptr;
 
-    ModuleBase::matrix vofk_effective;
+    ModuleBase::matrix vofk_eff;
 
     bool fixed_done = false;
 
diff --git a/source/source_estate/module_pot/potential_types.cpp b/source/source_estate/module_pot/potential_types.cpp
index f62a34aa75..ba74f5143e 100644
--- a/source/source_estate/module_pot/potential_types.cpp
+++ b/source/source_estate/module_pot/potential_types.cpp
@@ -13,7 +13,6 @@
 #include "pot_xc.h"
 #include "potential_new.h"
 #include "pot_sep.h"
-#include "pot_local_paw.h"
 #ifdef __LCAO
 #include "H_TDDFT_pw.h"
 #endif
@@ -34,13 +33,13 @@ PotBase* Potential::get_pot_type(const std::string& pot_type)
     }
     else if (pot_type == "xc")
     {
-        return new PotXC(this->rho_basis_, this->etxc_, this->vtxc_, &(this->vofk_effective));
+        return new PotXC(this->rho_basis_, this->etxc_, this->vtxc_, &(this->vofk_eff));
     }
     else if (pot_type == "surchem")
     {
         return new PotSurChem(this->rho_basis_,
                               this->structure_factors_,
-                              this->v_effective_fixed.data(),
+                              this->v_eff_fixed.data(),
                               this->solvent_);
     }
     else if (pot_type == "efield")
diff --git a/source/source_estate/setup_estate_pw.cpp b/source/source_estate/setup_estate_pw.cpp
index d384b3a171..5daa7be1b5 100644
--- a/source/source_estate/setup_estate_pw.cpp
+++ b/source/source_estate/setup_estate_pw.cpp
@@ -28,20 +28,15 @@ void elecstate::setup_estate_pw(UnitCell& ucell, // unitcell
         {
             //! SDFT only supports double precision currently
             pelec = new elecstate::ElecStatePW_SDFT<std::complex<double>, Device>(pw_wfc,
-                &chr, &kv, &ucell, &ppcell,
-                pw_rhod, pw_rho, pw_big);
+                &chr, &kv, &ucell, &ppcell, pw_rho, pw_big);
         }
         else
         {
             pelec = new elecstate::ElecStatePW<T, Device>(pw_wfc,
-                &chr, &kv, &ucell, &ppcell,
-                pw_rhod, pw_rho, pw_big);
+                &chr, &kv, &ucell, &ppcell, pw_rho, pw_big);
         }
     }
 
-    //! Inititlize the charge density.
-    chr.allocate(inp.nspin);
-
     //! Initialize DFT-1/2
     if (PARAM.inp.dfthalf_type > 0)
     {
diff --git a/source/source_estate/test/elecstate_base_test.cpp b/source/source_estate/test/elecstate_base_test.cpp
index 24ec4a5fff..8cc19b63b4 100644
--- a/source/source_estate/test/elecstate_base_test.cpp
+++ b/source/source_estate/test/elecstate_base_test.cpp
@@ -80,11 +80,7 @@ void ModulePW::PW_Basis::distribute_r()
 void Charge::set_rho_core(const UnitCell& ucell, ModuleBase::ComplexMatrix const&, const bool*)
 {
 }
-void Charge::set_rho_core_paw()
-{
-}
-void Charge::init_rho(elecstate::Efermi&,
-                      const UnitCell&,
+void Charge::init_rho(const UnitCell&,
                       const Parallel_Grid&,
                       ModuleBase::ComplexMatrix const&,
                       ModuleSymmetry::Symmetry& symm,
@@ -201,11 +197,10 @@ TEST_F(ElecStateTest, Constructor)
 TEST_F(ElecStateTest, InitKS)
 {
     Charge* charge = new Charge;
-    ModulePW::PW_Basis* rhopw = new ModulePW::PW_Basis;
     ModulePW::PW_Basis_Big* bigpw = new ModulePW::PW_Basis_Big;
     K_Vectors* klist = new K_Vectors;
     int nk = 1;
-    EXPECT_NO_THROW(elecstate->init_ks(charge, klist, nk, rhopw, bigpw));
+    EXPECT_NO_THROW(elecstate->init_ks(charge, klist, nk, bigpw));
     EXPECT_EQ(elecstate->charge, charge);
     EXPECT_EQ(elecstate->bigpw, bigpw);
     EXPECT_EQ(elecstate->klist, klist);
@@ -215,14 +210,12 @@ TEST_F(ElecStateTest, InitKS)
     EXPECT_EQ(elecstate->wg.nc, PARAM.input.nbands);
     delete klist;
     delete bigpw;
-    delete rhopw;
     delete charge;
 }
 
 TEST_F(ElecStateTest, GetRho)
 {
     Charge* charge = new Charge;
-    ModulePW::PW_Basis* rhopw = new ModulePW::PW_Basis;
     ModulePW::PW_Basis_Big* bigpw = new ModulePW::PW_Basis_Big;
     K_Vectors* klist = new K_Vectors;
     int nk = 1;
@@ -236,7 +229,7 @@ TEST_F(ElecStateTest, GetRho)
             charge->rho[i][j] = 1.0;
         }
     }
-    elecstate->init_ks(charge, klist, nk, rhopw, bigpw);
+    elecstate->init_ks(charge, klist, nk, bigpw);
     EXPECT_EQ(elecstate->getRho(0), &(charge->rho[0][0]));
     EXPECT_EQ(elecstate->getRho(0)[nrxx - 1], 1.0);
     for (int i = 0; i < PARAM.input.nspin; ++i)
@@ -246,7 +239,6 @@ TEST_F(ElecStateTest, GetRho)
     delete[] charge->rho;
     delete klist;
     delete bigpw;
-    delete rhopw;
     delete charge;
 }
 
diff --git a/source/source_estate/test/elecstate_pw_test.cpp b/source/source_estate/test/elecstate_pw_test.cpp
index 4128d12812..72fa56375b 100644
--- a/source/source_estate/test/elecstate_pw_test.cpp
+++ b/source/source_estate/test/elecstate_pw_test.cpp
@@ -130,11 +130,7 @@ Fcoef::~Fcoef()
 void Charge::set_rho_core(const UnitCell& ucell, ModuleBase::ComplexMatrix const&, const bool*)
 {
 }
-void Charge::set_rho_core_paw()
-{
-}
-void Charge::init_rho(elecstate::Efermi&,
-                      const UnitCell&,
+void Charge::init_rho(const UnitCell&,
                       const Parallel_Grid&,
                       ModuleBase::ComplexMatrix const&,
                       ModuleSymmetry::Symmetry& symm,
@@ -244,7 +240,6 @@ TEST_F(ElecStatePWTest, ConstructorDouble)
                                                                                                klist,
                                                                                                ucell,
                                                                                                ppcell,
-                                                                                               rhodpw,
                                                                                                rhopw,
                                                                                                bigpw);
     EXPECT_EQ(elecstate_pw_d->classname, "ElecStatePW");
@@ -260,7 +255,6 @@ TEST_F(ElecStatePWTest, ConstructorSingle)
                                                                                               klist,
                                                                                               ucell,
                                                                                               ppcell,
-                                                                                              rhodpw,
                                                                                               rhopw,
                                                                                               bigpw);
     EXPECT_EQ(elecstate_pw_s->classname, "ElecStatePW");
@@ -278,7 +272,6 @@ TEST_F(ElecStatePWTest, InitRhoDataDouble)
                                                                                                klist,
                                                                                                ucell,
                                                                                                ppcell,
-                                                                                               rhodpw,
                                                                                                rhopw,
                                                                                                bigpw);
     elecstate_pw_d->init_rho_data();
@@ -298,7 +291,6 @@ TEST_F(ElecStatePWTest, InitRhoDataSingle)
                                                                                               klist,
                                                                                               ucell,
                                                                                               ppcell,
-                                                                                              rhodpw,
                                                                                               rhopw,
                                                                                               bigpw);
     elecstate_pw_s->init_rho_data();
@@ -315,7 +307,6 @@ TEST_F(ElecStatePWTest, ParallelKDouble)
                                                                                                klist,
                                                                                                ucell,
                                                                                                ppcell,
-                                                                                               rhodpw,
                                                                                                rhopw,
                                                                                                bigpw);
     EXPECT_NO_THROW(elecstate_pw_d->parallelK());
@@ -329,7 +320,6 @@ TEST_F(ElecStatePWTest, ParallelKSingle)
                                                                                               klist,
                                                                                               ucell,
                                                                                               ppcell,
-                                                                                              rhodpw,
                                                                                               rhopw,
                                                                                               bigpw);
     EXPECT_NO_THROW(elecstate_pw_s->parallelK());
diff --git a/source/source_estate/test/potential_new_test.cpp b/source/source_estate/test/potential_new_test.cpp
index b874be7b39..ea260d248b 100644
--- a/source/source_estate/test/potential_new_test.cpp
+++ b/source/source_estate/test/potential_new_test.cpp
@@ -170,9 +170,9 @@ TEST_F(PotentialNewTest, ConstructorCPUDouble)
     pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
-    EXPECT_EQ(pot->v_effective_fixed.size(), 100);
-    EXPECT_EQ(pot->v_effective.nr, PARAM.input.nspin);
-    EXPECT_EQ(pot->v_effective.nc, 100);
+    EXPECT_EQ(pot->v_eff_fixed.size(), 100);
+    EXPECT_EQ(pot->v_eff.nr, PARAM.input.nspin);
+    EXPECT_EQ(pot->v_eff.nc, 100);
 }
 
 TEST_F(PotentialNewTest, ConstructorCPUSingle)
@@ -182,9 +182,9 @@ TEST_F(PotentialNewTest, ConstructorCPUSingle)
     pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
-    EXPECT_EQ(pot->v_effective_fixed.size(), 100);
-    EXPECT_EQ(pot->v_effective.nr, PARAM.input.nspin);
-    EXPECT_EQ(pot->v_effective.nc, 100);
+    EXPECT_EQ(pot->v_eff_fixed.size(), 100);
+    EXPECT_EQ(pot->v_eff.nr, PARAM.input.nspin);
+    EXPECT_EQ(pot->v_eff.nc, 100);
 }
 
 TEST_F(PotentialNewTest, ConstructorNRXX0)
@@ -203,11 +203,11 @@ TEST_F(PotentialNewTest, ConstructorXC3)
     pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
-    EXPECT_EQ(pot->v_effective_fixed.size(), 100);
-    EXPECT_EQ(pot->v_effective.nr, PARAM.input.nspin);
-    EXPECT_EQ(pot->v_effective.nc, 100);
-    EXPECT_EQ(pot->vofk_effective.nr, PARAM.input.nspin);
-    EXPECT_EQ(pot->vofk_effective.nc, 100);
+    EXPECT_EQ(pot->v_eff_fixed.size(), 100);
+    EXPECT_EQ(pot->v_eff.nr, PARAM.input.nspin);
+    EXPECT_EQ(pot->v_eff.nc, 100);
+    EXPECT_EQ(pot->vofk_eff.nr, PARAM.input.nspin);
+    EXPECT_EQ(pot->vofk_eff.nc, 100);
 }
 
 TEST_F(PotentialNewTest, ConstructorGPUDouble)
@@ -218,9 +218,9 @@ TEST_F(PotentialNewTest, ConstructorGPUDouble)
     pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
-    EXPECT_EQ(pot->v_effective_fixed.size(), 100);
-    EXPECT_EQ(pot->v_effective.nr, PARAM.input.nspin);
-    EXPECT_EQ(pot->v_effective.nc, 100);
+    EXPECT_EQ(pot->v_eff_fixed.size(), 100);
+    EXPECT_EQ(pot->v_eff.nr, PARAM.input.nspin);
+    EXPECT_EQ(pot->v_eff.nc, 100);
 }
 
 TEST_F(PotentialNewTest, ConstructorGPUSingle)
@@ -232,9 +232,9 @@ TEST_F(PotentialNewTest, ConstructorGPUSingle)
     pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
-    EXPECT_EQ(pot->v_effective_fixed.size(), 100);
-    EXPECT_EQ(pot->v_effective.nr, PARAM.input.nspin);
-    EXPECT_EQ(pot->v_effective.nc, 100);
+    EXPECT_EQ(pot->v_eff_fixed.size(), 100);
+    EXPECT_EQ(pot->v_eff.nr, PARAM.input.nspin);
+    EXPECT_EQ(pot->v_eff.nc, 100);
 }
 
 TEST_F(PotentialNewTest, Getters)
@@ -289,9 +289,9 @@ TEST_F(PotentialNewTest, CalFixedV)
     }
     double* vl_pseudo = new double[1000];
     pot->cal_fixed_v(vl_pseudo);
-    for (int i = 0; i < pot->v_effective_fixed.size(); i++)
+    for (int i = 0; i < pot->v_eff_fixed.size(); i++)
     {
-        EXPECT_DOUBLE_EQ(pot->v_effective_fixed[i], 0.0);
+        EXPECT_DOUBLE_EQ(pot->v_eff_fixed[i], 0.0);
     }
     delete[] vl_pseudo;
 }
@@ -319,9 +319,9 @@ TEST_F(PotentialNewTest, CalVeff)
     ModuleBase::matrix v_eff;
     v_eff.create(2, 100);
     pot->cal_v_eff(chg,this->ucell,v_eff);
-    for (int i = 0; i < pot->v_effective_fixed.size(); i++)
+    for (int i = 0; i < pot->v_eff_fixed.size(); i++)
     {
-        EXPECT_DOUBLE_EQ(pot->v_effective_fixed[i], 0.0);
+        EXPECT_DOUBLE_EQ(pot->v_eff_fixed[i], 0.0);
     }
     delete chg;
 }
@@ -417,8 +417,8 @@ TEST_F(PotentialNewTest, GetEffectiveVmatrix)
     rhopw->nrxx = 100;
     pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
-    ModuleBase::matrix v_eff_tmp = pot->get_effective_v();
-    const ModuleBase::matrix v_eff_tmp_const = pot->get_effective_v();
+    ModuleBase::matrix v_eff_tmp = pot->get_eff_v();
+    const ModuleBase::matrix v_eff_tmp_const = pot->get_eff_v();
     EXPECT_EQ(v_eff_tmp.nr, PARAM.input.nspin);
     EXPECT_EQ(v_eff_tmp.nc, 100);
     EXPECT_EQ(v_eff_tmp_const.nr, PARAM.input.nspin);
@@ -427,8 +427,8 @@ TEST_F(PotentialNewTest, GetEffectiveVmatrix)
     {
         for (int ic = 0; ic < v_eff_tmp.nc; ic++)
         {
-            EXPECT_DOUBLE_EQ(v_eff_tmp(ir, ic), pot->v_effective(ir, ic));
-            EXPECT_DOUBLE_EQ(v_eff_tmp_const(ir, ic), pot->v_effective(ir, ic));
+            EXPECT_DOUBLE_EQ(v_eff_tmp(ir, ic), pot->v_eff(ir, ic));
+            EXPECT_DOUBLE_EQ(v_eff_tmp_const(ir, ic), pot->v_eff(ir, ic));
         }
     }
 }
@@ -439,24 +439,24 @@ TEST_F(PotentialNewTest, GetEffectiveVarray)
     rhopw->nrxx = 100;
     pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
-    double* v_eff_tmp = pot->get_effective_v(0);
-    const double* v_eff_tmp_const = pot->get_effective_v(0);
+    double* v_eff_tmp = pot->get_eff_v(0);
+    const double* v_eff_tmp_const = pot->get_eff_v(0);
     for (int ic = 0; ic < rhopw->nrxx; ic++)
     {
-        EXPECT_DOUBLE_EQ(v_eff_tmp[ic], pot->v_effective(0, ic));
-        EXPECT_DOUBLE_EQ(v_eff_tmp_const[ic], pot->v_effective(0, ic));
+        EXPECT_DOUBLE_EQ(v_eff_tmp[ic], pot->v_eff(0, ic));
+        EXPECT_DOUBLE_EQ(v_eff_tmp_const[ic], pot->v_eff(0, ic));
     }
     v_eff_tmp[0] = 1.0;
-    EXPECT_DOUBLE_EQ(pot->v_effective(0, 0), 1.0);
+    EXPECT_DOUBLE_EQ(pot->v_eff(0, 0), 1.0);
     EXPECT_DOUBLE_EQ(v_eff_tmp_const[0], 1.0);
 }
 
 TEST_F(PotentialNewTest, GetEffectiveVarrayNullptr)
 {
     pot = new elecstate::Potential;
-    EXPECT_EQ(pot->v_effective.nc, 0);
-    double* v_eff_tmp = pot->get_effective_v(0);
-    const double* v_eff_tmp_const = pot->get_effective_v(0);
+    EXPECT_EQ(pot->v_eff.nc, 0);
+    double* v_eff_tmp = pot->get_eff_v(0);
+    const double* v_eff_tmp_const = pot->get_eff_v(0);
     EXPECT_EQ(v_eff_tmp, nullptr);
     EXPECT_EQ(v_eff_tmp_const, nullptr);
 }
@@ -469,8 +469,8 @@ TEST_F(PotentialNewTest, GetEffectiveVofkmatrix)
     rhopw->nrxx = 100;
     pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
-    ModuleBase::matrix vofk_eff_tmp = pot->get_effective_vofk();
-    const ModuleBase::matrix vofk_eff_tmp_const = pot->get_effective_vofk();
+    ModuleBase::matrix vofk_eff_tmp = pot->get_eff_vofk();
+    const ModuleBase::matrix vofk_eff_tmp_const = pot->get_eff_vofk();
     EXPECT_EQ(vofk_eff_tmp.nr, PARAM.input.nspin);
     EXPECT_EQ(vofk_eff_tmp.nc, 100);
     EXPECT_EQ(vofk_eff_tmp_const.nr, PARAM.input.nspin);
@@ -479,8 +479,8 @@ TEST_F(PotentialNewTest, GetEffectiveVofkmatrix)
     {
         for (int ic = 0; ic < vofk_eff_tmp.nc; ic++)
         {
-            EXPECT_DOUBLE_EQ(vofk_eff_tmp(ir, ic), pot->vofk_effective(ir, ic));
-            EXPECT_DOUBLE_EQ(vofk_eff_tmp_const(ir, ic), pot->vofk_effective(ir, ic));
+            EXPECT_DOUBLE_EQ(vofk_eff_tmp(ir, ic), pot->vofk_eff(ir, ic));
+            EXPECT_DOUBLE_EQ(vofk_eff_tmp_const(ir, ic), pot->vofk_eff(ir, ic));
         }
     }
 }
@@ -491,24 +491,24 @@ TEST_F(PotentialNewTest, GetEffectiveVofkarray)
     rhopw->nrxx = 100;
     pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
-    double* vofk_eff_tmp = pot->get_effective_vofk(0);
-    const double* vofk_eff_tmp_const = pot->get_effective_vofk(0);
+    double* vofk_eff_tmp = pot->get_eff_vofk(0);
+    const double* vofk_eff_tmp_const = pot->get_eff_vofk(0);
     for (int ic = 0; ic < rhopw->nrxx; ic++)
     {
-        EXPECT_DOUBLE_EQ(vofk_eff_tmp[ic], pot->vofk_effective(0, ic));
-        EXPECT_DOUBLE_EQ(vofk_eff_tmp_const[ic], pot->vofk_effective(0, ic));
+        EXPECT_DOUBLE_EQ(vofk_eff_tmp[ic], pot->vofk_eff(0, ic));
+        EXPECT_DOUBLE_EQ(vofk_eff_tmp_const[ic], pot->vofk_eff(0, ic));
     }
     vofk_eff_tmp[0] = 1.0;
-    EXPECT_DOUBLE_EQ(pot->vofk_effective(0, 0), 1.0);
+    EXPECT_DOUBLE_EQ(pot->vofk_eff(0, 0), 1.0);
     EXPECT_DOUBLE_EQ(vofk_eff_tmp_const[0], 1.0);
 }
 
 TEST_F(PotentialNewTest, GetEffectiveVofkarrayNullptr)
 {
     pot = new elecstate::Potential;
-    EXPECT_EQ(pot->v_effective.nc, 0);
-    double* vofk_eff_tmp = pot->get_effective_vofk(0);
-    const double* vofk_eff_tmp_const = pot->get_effective_vofk(0);
+    EXPECT_EQ(pot->v_eff.nc, 0);
+    double* vofk_eff_tmp = pot->get_eff_vofk(0);
+    const double* vofk_eff_tmp_const = pot->get_eff_vofk(0);
     EXPECT_EQ(vofk_eff_tmp, nullptr);
     EXPECT_EQ(vofk_eff_tmp_const, nullptr);
 }
@@ -519,14 +519,14 @@ TEST_F(PotentialNewTest, GetFixedV)
     pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
-    EXPECT_EQ(pot->v_effective_fixed.size(), 100);
+    EXPECT_EQ(pot->v_eff_fixed.size(), 100);
     double* v_eff_fixed_tmp = pot->get_fixed_v();
     const double* v_eff_fixed_tmp_const = pot->get_fixed_v();
     for (int ic = 0; ic < rhopw->nrxx; ic++)
     {
         v_eff_fixed_tmp[ic] = ic;
-        EXPECT_DOUBLE_EQ(v_eff_fixed_tmp[ic], pot->v_effective_fixed[ic]);
-        EXPECT_DOUBLE_EQ(v_eff_fixed_tmp_const[ic], pot->v_effective_fixed[ic]);
+        EXPECT_DOUBLE_EQ(v_eff_fixed_tmp[ic], pot->v_eff_fixed[ic]);
+        EXPECT_DOUBLE_EQ(v_eff_fixed_tmp_const[ic], pot->v_eff_fixed[ic]);
     }
 }
 
@@ -594,12 +594,12 @@ TEST_F(PotentialNewTest, InterpolateVrsDoubleGrids)
 
      pot = new elecstate::Potential(rhodpw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
 
-     for (int ir = 0; ir < pot->v_effective.nr; ir++)
+     for (int ir = 0; ir < pot->v_eff.nr; ir++)
      {
-         for (int ic = 0; ic < pot->v_effective.nc; ic++)
+         for (int ic = 0; ic < pot->v_eff.nc; ic++)
          {
-             pot->v_effective(ir, ic) = ir + ic;
-             pot->vofk_effective(ir, ic) = ir + 2 * ic;
+             pot->v_eff(ir, ic) = ir + ic;
+             pot->vofk_eff(ir, ic) = ir + 2 * ic;
          }
      }
 
@@ -655,12 +655,12 @@ TEST_F(PotentialNewTest, InterpolateVrsSingleGrids)
 
      pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
 
-     for (int ir = 0; ir < pot->v_effective.nr; ir++)
+     for (int ir = 0; ir < pot->v_eff.nr; ir++)
      {
-         for (int ic = 0; ic < pot->v_effective.nc; ic++)
+         for (int ic = 0; ic < pot->v_eff.nc; ic++)
          {
-             pot->v_effective(ir, ic) = ir + ic;
-             pot->vofk_effective(ir, ic) = ir + 2 * ic;
+             pot->v_eff(ir, ic) = ir + ic;
+             pot->vofk_eff(ir, ic) = ir + 2 * ic;
          }
      }
 
diff --git a/source/source_estate/update_pot.cpp b/source/source_estate/update_pot.cpp
index 8a49f66e0c..f0f0ef861a 100644
--- a/source/source_estate/update_pot.cpp
+++ b/source/source_estate/update_pot.cpp
@@ -12,9 +12,6 @@ void elecstate::update_pot(UnitCell& ucell, // unitcell
         elecstate::cal_ux(ucell);
         pelec->pot->update_from_charge(&chr, &ucell);
         pelec->f_en.descf = pelec->cal_delta_escf();
-#ifdef __MPI
-        MPI_Bcast(&(pelec->f_en.descf), 1, MPI_DOUBLE, 0, BP_WORLD);
-#endif
     }
     else
     {
diff --git a/source/source_hsolver/test/hsolver_supplementary_mock.h b/source/source_hsolver/test/hsolver_supplementary_mock.h
index b6beb9cb2d..3c3c500194 100644
--- a/source/source_hsolver/test/hsolver_supplementary_mock.h
+++ b/source/source_hsolver/test/hsolver_supplementary_mock.h
@@ -33,7 +33,6 @@ void ElecState::init_scf(const int istep,
 void ElecState::init_ks(Charge* chg_in, // pointer for class Charge
                         const K_Vectors* klist_in,
                         int nk_in,
-                        ModulePW::PW_Basis* rhopw_in,
                         const ModulePW::PW_Basis_Big* bigpw_in)
 {
     return;
@@ -45,7 +44,6 @@ ElecStatePW<T, Device>::ElecStatePW(ModulePW::PW_Basis_K* wfc_basis_in,
                                     K_Vectors* pkv_in,
                                     UnitCell* ucell_in,
                                     pseudopot_cell_vnl* ppcell_in,
-                                    ModulePW::PW_Basis* rhodpw_in,
                                     ModulePW::PW_Basis* rhopw_in,
                                     ModulePW::PW_Basis_Big* bigpw_in)
     : basis(wfc_basis_in)
diff --git a/source/source_hsolver/test/test_hsolver_sdft.cpp b/source/source_hsolver/test/test_hsolver_sdft.cpp
index 4118ba7bb7..7ab1464195 100644
--- a/source/source_hsolver/test/test_hsolver_sdft.cpp
+++ b/source/source_hsolver/test/test_hsolver_sdft.cpp
@@ -262,7 +262,7 @@ namespace ModulePW {
 class TestHSolverPW_SDFT : public ::testing::Test
 {
   public:
-    TestHSolverPW_SDFT() : stoche(8, 1, 0, 0), elecstate_test(nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr)
+    TestHSolverPW_SDFT() : stoche(8, 1, 0, 0), elecstate_test(nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr)
     {
     }
     ModulePW::PW_Basis_K pwbk;
@@ -408,4 +408,4 @@ int main(int argc, char** argv)
     return result;
 }
 
-#endif
\ No newline at end of file
+#endif
diff --git a/source/source_io/ctrl_output_fp.cpp b/source/source_io/ctrl_output_fp.cpp
index 0cde6e6011..bbe0fea792 100644
--- a/source/source_io/ctrl_output_fp.cpp
+++ b/source/source_io/ctrl_output_fp.cpp
@@ -121,7 +121,7 @@ void ctrl_output_fp(UnitCell& ucell,
             fn += spin_block + geom_block + ".cube";
 
             ModuleIO::write_vdata_palgrid(para_grid,
-                    pelec->pot->get_effective_v(is),
+                    pelec->pot->get_eff_v(is),
                     is,
                     nspin,
                     istep_in,
diff --git a/source/source_io/ctrl_output_pw.cpp b/source/source_io/ctrl_output_pw.cpp
index 34d57caa01..decd1f680a 100644
--- a/source/source_io/ctrl_output_pw.cpp
+++ b/source/source_io/ctrl_output_pw.cpp
@@ -366,7 +366,7 @@ void ModuleIO::ctrl_runner_pw(UnitCell& ucell,
                                                pw_wfc,
                                                pw_rho,
                                                ucell,
-                                               pelec->pot->get_effective_v(0));
+                                               pelec->pot->get_eff_v(0));
     }
 #endif
 
diff --git a/source/source_io/ctrl_scf_lcao.cpp b/source/source_io/ctrl_scf_lcao.cpp
index 889e50f926..9439f15576 100644
--- a/source/source_io/ctrl_scf_lcao.cpp
+++ b/source/source_io/ctrl_scf_lcao.cpp
@@ -216,7 +216,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
 			inp.out_mat_t,
 			inp.out_mat_r,
 			istep,
-			pelec->pot->get_effective_v(),
+			pelec->pot->get_eff_v(),
 			pv,
 			two_center_bundle,
 			orb,
diff --git a/source/source_io/write_elecstat_pot.h b/source/source_io/write_elecstat_pot.h
index 76c90bf70d..bd3f2a8734 100644
--- a/source/source_io/write_elecstat_pot.h
+++ b/source/source_io/write_elecstat_pot.h
@@ -18,7 +18,7 @@ namespace ModuleIO
 /// @param rho_basis
 /// @param chr
 /// @param ucell_
-/// @param v_effective_fixed
+/// @param v_eff_fixed
 void write_elecstat_pot(
 #ifdef __MPI
     const int& bz,
@@ -29,7 +29,7 @@ void write_elecstat_pot(
     ModulePW::PW_Basis* rho_basis,
     const Charge* const chr,
     const UnitCell* ucell_,
-    const double* v_effective_fixed,
+    const double* v_eff_fixed,
     const surchem& solvent);
 
 } // namespace ModuleIO
diff --git a/source/source_io/write_init.cpp b/source/source_io/write_init.cpp
index 8881647023..ba435790f6 100644
--- a/source/source_io/write_init.cpp
+++ b/source/source_io/write_init.cpp
@@ -91,7 +91,7 @@ void ModuleIO::write_pot_init(
             }
 
             ModuleIO::write_vdata_palgrid(para_grid,
-                                          pelec->pot->get_effective_v(is),
+                                          pelec->pot->get_eff_v(is),
                                           is,
                                           nspin,
                                           istep,
diff --git a/source/source_lcao/LCAO_set.cpp b/source/source_lcao/LCAO_set.cpp
index 2c87852fde..ebb296ffd2 100644
--- a/source/source_lcao/LCAO_set.cpp
+++ b/source/source_lcao/LCAO_set.cpp
@@ -39,9 +39,6 @@ void LCAO_domain::set_psi_occ_dm_chg(
     //! 4) init DMK, but DMR is constructed in before_scf()
     dmat.allocate_dm(&kv, &pv, inp.nspin);
 
-    //! 5) init charge density
-    chr.allocate(inp.nspin);
-
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "CHARGE");
 
     return;
diff --git a/source/source_lcao/module_operator_lcao/veff_lcao.cpp b/source/source_lcao/module_operator_lcao/veff_lcao.cpp
index 0df6ed33a5..f664844544 100644
--- a/source/source_lcao/module_operator_lcao/veff_lcao.cpp
+++ b/source/source_lcao/module_operator_lcao/veff_lcao.cpp
@@ -65,8 +65,8 @@ void Veff<OperatorLCAO<double, double>>::contributeHR()
     //(1) prepare data for this k point.
     // copy the local potential from array.
     //-----------------------------------------
-    double* vr_eff1 = this->pot->get_effective_v(this->current_spin);
-    double* vofk_eff1 = this->pot->get_effective_vofk(this->current_spin);
+    double* vr_eff1 = this->pot->get_eff_v(this->current_spin);
+    double* vofk_eff1 = this->pot->get_eff_vofk(this->current_spin);
 
     if(XC_Functional::get_ked_flag())
     {
@@ -95,8 +95,8 @@ void Veff<OperatorLCAO<std::complex<double>, double>>::contributeHR()
     //(1) prepare data for this k point.
     // copy the local potential from array.
     //-----------------------------------------
-    double* vr_eff1 = this->pot->get_effective_v(this->current_spin);
-    double* vofk_eff1 = this->pot->get_effective_vofk(this->current_spin);
+    double* vr_eff1 = this->pot->get_eff_v(this->current_spin);
+    double* vofk_eff1 = this->pot->get_eff_vofk(this->current_spin);
 
     if(XC_Functional::get_ked_flag())
     {
@@ -126,16 +126,17 @@ void Veff<OperatorLCAO<std::complex<double>, std::complex<double>>>::contributeH
     std::vector<const double*> vofk_eff(4, nullptr);
     for (int is = 0; is < 4; is++)
     {
-        vr_eff[is] = this->pot->get_effective_v(is);
+        vr_eff[is] = this->pot->get_eff_v(is);
         if(XC_Functional::get_ked_flag())
         {
-            vofk_eff[is] = this->pot->get_effective_vofk(is);
+            vofk_eff[is] = this->pot->get_eff_vofk(is);
         }
     }
     if(XC_Functional::get_ked_flag())
     {
         ModuleGint::cal_gint_vl_metagga(vr_eff, vofk_eff, this->hR);
-    } else
+    } 
+    else
     {
         ModuleGint::cal_gint_vl(vr_eff, this->hR);
     }
@@ -150,4 +151,4 @@ template class Veff<OperatorLCAO<double, double>>;
 template class Veff<OperatorLCAO<std::complex<double>, double>>;
 
 template class Veff<OperatorLCAO<std::complex<double>, std::complex<double>>>;
-}
\ No newline at end of file
+}
diff --git a/source/source_lcao/pulay_fs_gint.hpp b/source/source_lcao/pulay_fs_gint.hpp
index 9603a28c34..46040ce340 100644
--- a/source/source_lcao/pulay_fs_gint.hpp
+++ b/source/source_lcao/pulay_fs_gint.hpp
@@ -24,8 +24,8 @@ namespace PulayForceStress
         {
             for (int is = 0; is < nspin; ++is)
             {
-                vr_eff[is] = pot->get_effective_v(is);
-                vofk_eff[is] = pot->get_effective_vofk(is);
+                vr_eff[is] = pot->get_eff_v(is);
+                vofk_eff[is] = pot->get_eff_vofk(is);
             }
             ModuleGint::cal_gint_fvl_meta(nspin, vr_eff, vofk_eff, dm.get_DMR_vector(), isforce, isstress, &f, &s);
         }
@@ -33,7 +33,7 @@ namespace PulayForceStress
         {
             for(int is = 0; is < nspin; ++is)
             {
-                vr_eff[is] = pot->get_effective_v(is);
+                vr_eff[is] = pot->get_eff_v(is);
             }
             ModuleGint::cal_gint_fvl(nspin, vr_eff, dm.get_DMR_vector(), isforce, isstress, &f, &s);
         }
diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp
index f4aaea2cec..05a058671d 100644
--- a/source/source_pw/module_ofdft/evolve_ofdft.cpp
+++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp
@@ -25,10 +25,10 @@ void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec,
     }
 
     pelec->pot->update_from_charge(&chr, &ucell); // Hartree + XC + external
-    this->cal_tf_potential(chr.rho, pw_rho, pelec->pot->get_effective_v()); // TF potential
+    this->cal_tf_potential(chr.rho, pw_rho, pelec->pot->get_eff_v()); // TF potential
     if (PARAM.inp.of_cd)
     {
-        this->cal_CD_potential(psi_, pw_rho, pelec->pot->get_effective_v(), PARAM.inp.of_mCD_alpha); // CD potential
+        this->cal_CD_potential(psi_, pw_rho, pelec->pot->get_eff_v(), PARAM.inp.of_mCD_alpha); // CD potential
     }
 
 #ifdef _OPENMP
@@ -36,7 +36,7 @@ void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec,
 #endif
     for (int is = 0; is < PARAM.inp.nspin; ++is)
     {
-        const double* vr_eff = pelec->pot->get_effective_v(is);
+        const double* vr_eff = pelec->pot->get_eff_v(is);
         for (int ir = 0; ir < pw_rho->nrxx; ++ir)
         {
             Hpsi[is * pw_rho->nrxx + ir] = vr_eff[ir]*psi_[is * pw_rho->nrxx + ir];
@@ -188,14 +188,17 @@ void Evolve_OFDFT::cal_CD_potential(std::vector<std::complex<double>> psi_,
         pw_rho->real2recip(rCurrent_z,recipCurrent_z);
         for (int ik = 0; ik < pw_rho->npw; ++ik)
         {
-            recipCDPotential[ik]=recipCurrent_x[ik]*pw_rho->gcar[ik].x+recipCurrent_y[ik]*pw_rho->gcar[ik].y+recipCurrent_z[ik]*pw_rho->gcar[ik].z;
+            recipCDPotential[ik]=recipCurrent_x[ik]*pw_rho->gcar[ik].x
+              +recipCurrent_y[ik]*pw_rho->gcar[ik].y
+              +recipCurrent_z[ik]*pw_rho->gcar[ik].z;
             recipCDPotential[ik]*=imag/pw_rho->gg[ik];
         }
         pw_rho->recip2real(recipCDPotential,rCDPotential);
 
         for (int ir = 0; ir < pw_rho->nrxx; ++ir)
         {
-            rpot(0, ir) -= mCD_para*2.0*std::real(rCDPotential[ir])*std::pow(ModuleBase::PI,3) / (2.0*std::pow(std::real(kF_r[ir]),2));
+            rpot(0, ir) -= mCD_para*2.0*std::real(rCDPotential[ir])*std::pow(ModuleBase::PI,3) 
+              / (2.0*std::pow(std::real(kF_r[ir]),2));
         }
         delete[] recipCurrent_x;
         delete[] recipCurrent_y;
@@ -278,9 +281,11 @@ void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec,
         for (int ir = 0; ir < pw_rho->nrxx; ++ir)
         {
             K4[is * nrxx + ir]=-1.0*K4[is * nrxx + ir]*dt*imag;
-            pphi_[is * nrxx + ir]+=1.0/6.0*(K1[is * nrxx + ir]+2.0*K2[is * nrxx + ir]+2.0*K3[is * nrxx + ir]+K4[is * nrxx + ir]);
+            pphi_[is * nrxx + ir]+=1.0/6.0*(K1[is * nrxx + ir]
+              +2.0*K2[is * nrxx + ir]+2.0*K3[is * nrxx + ir]
+              +K4[is * nrxx + ir]);
         }
     }
 
     ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi");
-}
\ No newline at end of file
+}
diff --git a/source/source_pw/module_pwdft/forces_us.cpp b/source/source_pw/module_pwdft/forces_us.cpp
index dd2fbc5f8d..f5a6a4444e 100644
--- a/source/source_pw/module_pwdft/forces_us.cpp
+++ b/source/source_pw/module_pwdft/forces_us.cpp
@@ -31,7 +31,7 @@ void Forces<FPTYPE, Device>::cal_force_us(ModuleBase::matrix& forcenl,
 
     ModuleBase::matrix forceq(ucell.nat, 3);
 
-    ModuleBase::matrix veff = elec.pot->get_effective_v();
+    ModuleBase::matrix veff = elec.pot->get_eff_v();
     ModuleBase::ComplexMatrix vg(PARAM.inp.nspin, npw);
     // fourier transform of the total effective potential
     for (int is = 0; is < PARAM.inp.nspin; is++)
@@ -143,4 +143,4 @@ void Forces<FPTYPE, Device>::cal_force_us(ModuleBase::matrix& forcenl,
 template class Forces<double, base_device::DEVICE_CPU>;
 #if ((defined __CUDA) || (defined __ROCM))
 template class Forces<double, base_device::DEVICE_GPU>;
-#endif
\ No newline at end of file
+#endif
diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp
index e4bc159c80..7a3b0c711a 100644
--- a/source/source_pw/module_pwdft/setup_pot.cpp
+++ b/source/source_pw/module_pwdft/setup_pot.cpp
@@ -66,7 +66,7 @@ void pw::setup_pot(const int istep,
 	//! D in uspp need vloc, thus behind init_scf()
 	//! calculate the effective coefficient matrix
 	//! for non-local pseudopotential projectors
-	ModuleBase::matrix veff = pelec->pot->get_effective_v();
+	ModuleBase::matrix veff = pelec->pot->get_eff_v();
 
 	ppcell.cal_effective_D(veff, pw_rhod, ucell);
 
diff --git a/source/source_pw/module_pwdft/stress_func_us.cpp b/source/source_pw/module_pwdft/stress_func_us.cpp
index adb8745901..1881dacafe 100644
--- a/source/source_pw/module_pwdft/stress_func_us.cpp
+++ b/source/source_pw/module_pwdft/stress_func_us.cpp
@@ -26,7 +26,7 @@ void Stress_PW<FPTYPE, Device>::stress_us(ModuleBase::matrix& sigma,
 
     ModuleBase::matrix stressus(3, 3);
 
-    ModuleBase::matrix veff = this->pelec->pot->get_effective_v();
+    ModuleBase::matrix veff = this->pelec->pot->get_eff_v();
     ModuleBase::ComplexMatrix vg(PARAM.inp.nspin, npw);
     // fourier transform of the total effective potential
     for (int is = 0; is < PARAM.inp.nspin; is++)
diff --git a/source/source_pw/module_pwdft/stress_pw.cpp b/source/source_pw/module_pwdft/stress_pw.cpp
index ee069ad898..f79bd6349c 100644
--- a/source/source_pw/module_pwdft/stress_pw.cpp
+++ b/source/source_pw/module_pwdft/stress_pw.cpp
@@ -92,7 +92,7 @@ void Stress_PW<FPTYPE, Device>::cal_stress(ModuleBase::matrix& sigmatot,
         this->stress_mgga(ucell,
                           sigmaxc,
                           this->pelec->wg,
-                          this->pelec->pot->get_effective_vofk(),
+                          this->pelec->pot->get_eff_vofk(),
                           pelec->charge,
                           p_kv,
                           wfc_basis,