deepmodeling · mohanchen · Sep 21, 2026 · Sep 16, 2026 · Sep 16, 2026 · Sep 16, 2026
diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml
@@ -83,7 +83,7 @@ jobs:
 
         lcov --directory . --capture --output-file coverage.info
 
-        lcov --remove coverage.info '/usr/*' '*/test/*' '*/external/*' '*/build/*' --output-file coverage.filtered.info
+        lcov --remove coverage.info '/usr/*' '*/test/*' '*/unittests/*' '*/external/*' '*/build/*' --output-file coverage.filtered.info
 
         genhtml coverage.filtered.info --output-directory coverage-report
 

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
@@ -212,6 +212,13 @@ jobs:
         run: |
           ctest --test-dir build -V --timeout 1700 -R MODULE_ESTATE
 
+      - name: Module_Charge Unittests
+        env:
+          GTEST_COLOR: 'yes'
+          OMP_NUM_THREADS: '2'
+        run: |
+          ctest --test-dir build -V --timeout 1700 -R MODULE_CHARGE
+
       - name: Module_Hamilt Unittests
         env:
           GTEST_COLOR: 'yes'
@@ -343,4 +350,4 @@ jobs:
           GTEST_COLOR: 'yes'
           OMP_NUM_THREADS: '2'
         run: |
-          ctest --test-dir build -V --timeout 1700 -E 'integrate_test|01_PW|02_NAO_Gamma|03_NAO_multik|04_FF|05_rtTDDFT|06_SDFT|07_OFDFT|08_EXX|09_DeePKS|10_others|11_PW_GPU|12_NAO_Gamma_GPU|13_NAO_multik_GPU|15_rtTDDFT_GPU|16_SDFT_GPU|17_DS_DFTU|MODULE_BASE|MODULE_IO|MODULE_HSOLVER|MODULE_CELL|MODULE_MD|MODULE_PSI|MODULE_ESTATE|MODULE_RI|MODULE_HAMILT|MODULE_PW|MODULE_LCAO|MODULE_AO|MODULE_NAO|MODULE_RELAX|MODULE_LR'
+          ctest --test-dir build -V --timeout 1700 -E 'integrate_test|01_PW|02_NAO_Gamma|03_NAO_multik|04_FF|05_rtTDDFT|06_SDFT|07_OFDFT|08_EXX|09_DeePKS|10_others|11_PW_GPU|12_NAO_Gamma_GPU|13_NAO_multik_GPU|15_rtTDDFT_GPU|16_SDFT_GPU|17_DS_DFTU|MODULE_BASE|MODULE_IO|MODULE_HSOLVER|MODULE_CELL|MODULE_MD|MODULE_PSI|MODULE_ESTATE|MODULE_CHARGE|MODULE_RI|MODULE_HAMILT|MODULE_PW|MODULE_LCAO|MODULE_AO|MODULE_NAO|MODULE_RELAX|MODULE_LR'
diff --git a/source/Makefile b/source/Makefile
@@ -259,14 +259,6 @@ ${OBJ_DIR}/parse_args.o: $(ABACUS_ROOT)source_io/parse_args.cpp $(BUILD_INFO_DIR
 	@mkdir -p $(dir $@)
 	${CXX} ${OPTS} ${OPTS_MPI} -c ${HONG} $< -o $@
 
-# Explicit rule: source_cell/module_symmetry and source_estate/module_charge
-# both contain a symm_rho.cpp. VPATH resolves symm_rho.o to the source_cell
-# one (listed first), so the Symmetry_rho implementation needs an explicit
-# object name to be compiled from the correct source.
-${OBJ_DIR}/symm_rho_charge.o: $(ABACUS_ROOT)source_estate/module_charge/symm_rho.cpp
-	@mkdir -p $(dir $@)
-	${CXX} ${OPTS} ${OPTS_MPI} -c ${HONG} $< -o $@
-
 ###### END of ABACUS INFO PART ######
 
 #==========================

diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -330,6 +330,7 @@ OBJS_GINT=batch_biggrid.o\
       gint_fvl.o\
       gint_info.o\
       gint_interface.o\
+      gint_prec_ctrl.o\
       gint_rho.o\
       gint_drho.o\
       gint_tau.o\
@@ -825,19 +826,25 @@ OBJS_SRCPW=h_ewald_pw.o\
     vnl_pw_grad.o\
     vnl_pw_init_vnl.o\
     vnl_pw_qrad.o\
-    chgmixing.o\
-    gint_prec_ctrl.o\
+    chg_routine.o\
     charge.o\
-    charge_init.o\
-    charge_mpi.o\
-    charge_extra.o\
-    charge_mixing.o\
-    charge_mixing_dmr.o\
-    mix_resid.o\
-    mix_precond.o\
-    charge_mixing_rho.o\
-    charge_mixing_uspp.o\
-    rhog_io.o\
+    chg_atomic.o\
+    chg_atomic_inner.o\
+    chg_init.o\
+    chg_tools.o\
+    chg_parallel.o\
+    chg_extra.o\
+    chg_mix.o\
+    chg_mix_rho.o\
+    chg_dmr.o\
+    chg_drho.o\
+    chg_drho_inner.o\
+    chg_precond.o\
+    chg_tau.o\
+    chg_symm.o\
+    chg_symm_detail.o\
+    chg_uspp.o\
+    chg_rhog_io.o\
     write_elecstat_pot.o\
     write_init.o\
     fp_energy.o\
@@ -894,8 +901,6 @@ OBJS_SRCPW=h_ewald_pw.o\
     stress_pw.o\
     of_stress_pw.o\
     of_print_info.o\
-    symm_rho_charge.o\
-    symm_rhog.o\
     setup_psi_pw.o\
     setup_psi.o\
     psi_prepare.o\

diff --git a/source/source_base/parallel_grid.cpp b/source/source_base/parallel_grid.cpp
@@ -352,7 +352,12 @@ void Parallel_Grid::reduce(double* rhotot, const double* const rhoin, const bool
         return;
     }
 
-    assert(rhoin != nullptr);
+    // A rank may own zero real-space grid points (nrxx == 0) when the grid is
+    // decomposed across more processes than it has slabs. In that case the
+    // source buffer is legitimately null: MPI_Gatherv is called with
+    // sendcount 0 below and ignores the send buffer. Only a null buffer with a
+    // non-zero nrxx is a genuine bug.
+    assert(rhoin != nullptr || this->nrxx == 0);
     assert(this->nrxx == this->ncxy * this->nczp);
 
     int pool_size = 0;

diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp
@@ -45,7 +45,7 @@ class XC_First_Order_FDM : public ModuleDFPT::XC_First_Order
         fdm_ = new elecstate::PotXC_FDM(rho_basis, chg0, ucell);
         chg1_ = new Charge();
         chg1_->set_rhopw(rho_basis);
-        chg1_->allocate(chg0->nspin, false);
+        chg1_->allocate(chg0->nspin, false, false, PARAM.inp.test_charge);
         veff_1_.create(chg0->nspin, chg0->nrxx);
     }
 

diff --git a/source/source_esolver/esolver_dm2rho.cpp b/source/source_esolver/esolver_dm2rho.cpp
@@ -66,7 +66,7 @@ void ESolver_DM2rho<TK, TR>::runner(BaseCell& basecell, const int istep)
 
     // it's dangerous to design psiToRho function like this, mohan note 20251024
     // this->pelec->psiToRho(*this->psi);
-    LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), this->inp_->nspin, &this->chr);
+    LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), this->inp_->nspin, &this->chr, this->inp_->nelec, ucell.omega, false);
 
     int nspin0 = this->inp_->nspin == 2 ? 2 : 1;
 

diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp
@@ -12,6 +12,8 @@
 //-----HSolver ElecState Hamilt--------
 #include "source_estate/elecstate_lcao.h"
 #include "source_estate/elecstate_tools.h"
+#include "source_estate/module_charge/chg_init.h"
+#include "source_estate/module_charge/chg_tools.h"
 #include "source_hsolver/hsolver_lcao.h"
 #include "source_io/module_parameter/parameter.h"
 #include "source_io/module_restart/restart.h" // GlobalC::restart for load_exx_flag
@@ -90,11 +92,26 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(BaseCell& basecell, const Inpu
     this->dmat_base.allocate_dm(&this->kv, &this->pv, this->inp_->nspin);
 
     // 10) inititlize the charge density
+    module_charge::InitRhoCfg init_rho_cfg;
+    init_rho_cfg.init_chg = this->inp_->init_chg;
+    init_rho_cfg.suffix = this->inp_->suffix;
+    init_rho_cfg.esolver_type = this->inp_->esolver_type;
+    init_rho_cfg.global_readin_dir = PARAM.globalv.global_readin_dir;
+    init_rho_cfg.nelec = this->inp_->nelec;
+    init_rho_cfg.nbands = this->inp_->nbands;
+    init_rho_cfg.test_charge = this->inp_->test_charge;
+    init_rho_cfg.domag = PARAM.globalv.domag;
+    init_rho_cfg.domag_z = PARAM.globalv.domag_z;
+    init_rho_cfg.npol = PARAM.globalv.npol;
+    init_rho_cfg.meta_gga = XC_Functional::get_ked_flag();
     this->chr_base.set_rhopw(this->pw_rhod);           // mohan add 20251130
-    const bool kin_den = this->chr_base.kin_density(); // mohan add 20251202
-    this->chr_base.allocate(this->inp_->nspin, kin_den);
-    this->chr_base.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv);
-    this->chr_base.check_rho();
+    const bool kin_den = XC_Functional::get_ked_flag() || (this->inp_->out_elf[0] > 0); // mohan add 20251202
+    this->chr_base.allocate(this->inp_->nspin, kin_den, XC_Functional::get_ked_flag(),
+                            this->inp_->test_charge);
+    this->chr_base.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv, nullptr, init_rho_cfg);
+    module_charge::check_rho(this->chr_base.rho, this->chr_base.nspin,
+                             this->chr_base.rhopw->nrxx, ucell.omega,
+                             this->chr_base.rhopw->nxyz, this->inp_->nelec);
 
     // 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
@@ -2,14 +2,15 @@
 
 #include "source_base/tool_quit.h"
 #include "source_cell/cal_ux.h"
-#include "source_estate/module_charge/symm_rho.h"
+#include "source_estate/module_charge/chg_atomic.h"
+#include "source_estate/module_charge/chg_symm.h"
 #include "source_cell/read_pp_ucell.h"
 #include "source_estate/param_update.h"
 #include "source_hamilt/module_ewald/h_ewald_pw.h"
 #include "source_hamilt/module_vdw/vdw.h"
 #include "source_io/module_output/output_log.h"
 #include "source_io/module_output/print_info.h"
-#include "source_estate/rhog_io.h"
+#include "source_estate/module_charge/chg_rhog_io.h"
 #include "source_io/module_parameter/parameter.h"
 
 #include "source_pw/module_pwdft/setup_pwrho.h" // mohan 20251005
@@ -147,8 +148,8 @@ void ESolver_FP::before_all_runners(BaseCell& basecell, const Input_para& inp)
     //! 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
-    const bool kin_den = this->chr.kin_density(); // mohan add 20251202
-	this->chr.allocate(inp.nspin, kin_den); // mohan move this from setup_estate_pw, 20251128
+    const bool kin_den = XC_Functional::get_ked_flag() || (inp.out_elf[0] > 0); // mohan add 20251202
+	this->chr.allocate(inp.nspin, kin_den, XC_Functional::get_ked_flag(), inp.test_charge); // mohan move this from setup_estate_pw, 20251128
 
 
     return;
@@ -165,7 +166,13 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
     ModuleIO::output_efermi(conv_esolver, this->pelec->eferm.ef);
 
     //! Update delta_rho for charge extrapolation
-    CE.update_delta_rho(ucell, &(this->chr), &(this->sf));
+    const module_charge::AtomicRhoCfg atomic_rho_cfg_after{
+        PARAM.inp.nelec,
+        PARAM.inp.test_charge,
+        PARAM.globalv.domag,
+        PARAM.globalv.domag_z,
+        GlobalV::ofs_warning};
+    CE.update_delta_rho(ucell, &(this->chr), *this->pw_rhod, &(this->sf), atomic_rho_cfg_after);
 
     //! print out charge density, potential, elf, etc.
 	ModuleIO::ctrl_output_fp(ucell, *this->inp_, this->pelec, this->pw_big, this->pw_rhod, 
@@ -218,8 +225,15 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
     if (ucell.ionic_position_updated)
     {
         this->CE.update_all_dis(ucell);
-        this->CE.extrapolate_charge(&this->Pgrid, ucell, &this->chr, &this->sf,
-                                    GlobalV::ofs_running, GlobalV::ofs_warning);
+        const module_charge::AtomicRhoCfg atomic_rho_cfg_before{
+            PARAM.inp.nelec,
+            PARAM.inp.test_charge,
+            PARAM.globalv.domag,
+            PARAM.globalv.domag_z,
+            GlobalV::ofs_warning};
+        this->CE.extrapolate_charge(&this->Pgrid, ucell, &this->chr, *this->pw_rhod,
+                                    &this->sf, GlobalV::ofs_running, GlobalV::ofs_warning,
+                                    atomic_rho_cfg_before);
     }
 
     //! Evaluate the vdW correction once for this ionic configuration.
@@ -268,7 +282,7 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool&
             // Only pool 0 writes the rhog file (rhog is identical across pools).
             if (GlobalV::MY_POOL == 0)
             {
-                elecstate::write_rhog(PARAM.globalv.global_out_dir + this->inp_->suffix + "-CHARGE-DENSITY.restart",
+                module_charge::write_rhog(PARAM.globalv.global_out_dir + this->inp_->suffix + "-CHARGE-DENSITY.restart",
                                      PARAM.globalv.gamma_only_pw,
                                      this->pw_rhod,
                                      this->inp_->nspin,
@@ -289,7 +303,7 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool&
                 }
                 if (GlobalV::MY_POOL == 0)
                 {
-                    elecstate::write_rhog(PARAM.globalv.global_out_dir + this->inp_->suffix + "-TAU-DENSITY.restart",
+                    module_charge::write_rhog(PARAM.globalv.global_out_dir + this->inp_->suffix + "-TAU-DENSITY.restart",
                                          PARAM.globalv.gamma_only_pw,
                                          this->pw_rhod,
                                          this->inp_->nspin,

diff --git a/source/source_esolver/esolver_fp.h b/source/source_esolver/esolver_fp.h
@@ -5,7 +5,7 @@
 #include "source_base/timer_wrapper.h"
 #include "source_basis/module_pw/pw_basis.h"          // plane wave basis
 #include "source_estate/elecstate.h"                  // electronic states
-#include "source_estate/module_charge/charge_extra.h" // charge extrapolation
+#include "source_estate/module_charge/chg_extra.h" // charge extrapolation
 #include "source_hamilt/module_surchem/surchem.h"     // solvation model
 #include "source_base/parallel_grid.h"     // Parallel_Grid (value member below)
 #include "source_pw/module_pwdft/stru_fac.h"  // structure factor
@@ -81,7 +81,7 @@ class ESolver_FP : public ESolver
     pseudopot_cell_vl locpp;
 
     //! charge extrapolation method
-    Charge_Extra CE;
+    module_charge::Charge_Extra CE;
 
     //! solvent model
     surchem solvent;

diff --git a/source/source_esolver/esolver_gets.h b/source/source_esolver/esolver_gets.h
@@ -1,6 +1,7 @@
 #ifndef ESOLVER_GETS_H
 #define ESOLVER_GETS_H
 
+#include "source_basis/module_ao/parallel_orbitals.h"
 #include "source_basis/module_nao/two_center_bundle.h"
 #include "source_cell/unitcell.h"
 #include "source_esolver/esolver_ks.h"

diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp
@@ -5,7 +5,10 @@
 #include "source_io/module_json/output_info.h"
 
 #include "source_estate/update_pot.h" // mohan add 20251016
-#include "source_estate/module_charge/chgmixing.h" // mohan add 20251018
+#include "source_estate/module_charge/chg_routine.h" // mohan add 20251018
+#include "source_estate/module_charge/chg_drho.h" // module_charge::cal_drho/cal_dkin
+#include "source_estate/module_charge/chg_init.h" // module_charge::InitRhoCfg
+#include "source_estate/module_charge/chg_tools.h" // module_charge::check_rho
 #include "source_pw/module_pwdft/setup_pwwfc.h" // mohan add 20251018
 #include "source_hsolver/hsolver.h"
 #include "source_io/module_energy/write_eig_occ.h"
@@ -65,18 +68,48 @@ void ESolver_KS::before_all_runners(BaseCell& basecell, const Input_para& inp)
     //! 3) setup charge mixing
     p_chgmix = new Charge_Mixing();
     p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod);
-    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);
+    MixingConfig mix_cfg;
+    mix_cfg.mixing_mode = inp.mixing_mode;
+    mix_cfg.mixing_beta = inp.mixing_beta;
+    mix_cfg.mixing_ndim = inp.mixing_ndim;
+    mix_cfg.mixing_gg0 = inp.mixing_gg0;
+    // tau mixing is only meaningful for kinetic-energy-density functionals
+    mix_cfg.mixing_tau = inp.mixing_tau && XC_Functional::get_ked_flag();
+    mix_cfg.mixing_beta_mag = inp.mixing_beta_mag;
+    mix_cfg.mixing_gg0_mag = inp.mixing_gg0_mag;
+    mix_cfg.mixing_gg0_min = inp.mixing_gg0_min;
+    mix_cfg.mixing_angle = inp.mixing_angle;
+    mix_cfg.mixing_dmr = inp.mixing_dmr;
+    mix_cfg.nspin = inp.nspin;
+    mix_cfg.scf_thr_type = inp.scf_thr_type;
+    mix_cfg.double_grid = PARAM.globalv.double_grid;
+    mix_cfg.gamma_only_pw = PARAM.globalv.gamma_only_pw;
+    mix_cfg.domag = PARAM.globalv.domag;
+    mix_cfg.domag_z = PARAM.globalv.domag_z;
+    mix_cfg.scf_nmax = inp.scf_nmax;
+    p_chgmix->set_mixing(mix_cfg, ucell.omega, ucell.tpiba);
     p_chgmix->init_mixing();
 
     //! 4) setup plane wave for electronic wave functions
     pw::setup_pwwfc(inp, ucell, *this->pw_rho, this->kv, this->pw_wfc);
 
     //! 5) 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
+    module_charge::InitRhoCfg init_rho_cfg;
+    init_rho_cfg.init_chg = inp.init_chg;
+    init_rho_cfg.suffix = inp.suffix;
+    init_rho_cfg.esolver_type = inp.esolver_type;
+    init_rho_cfg.global_readin_dir = PARAM.globalv.global_readin_dir;
+    init_rho_cfg.nelec = inp.nelec;
+    init_rho_cfg.nbands = inp.nbands;
+    init_rho_cfg.test_charge = inp.test_charge;
+    init_rho_cfg.domag = PARAM.globalv.domag;
+    init_rho_cfg.domag_z = PARAM.globalv.domag_z;
+    init_rho_cfg.npol = PARAM.globalv.npol;
+    init_rho_cfg.meta_gga = XC_Functional::get_ked_flag();
+    this->chr.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv, this->pw_wfc, init_rho_cfg);
+    module_charge::check_rho(this->chr.rho, this->chr.nspin, this->chr.rhopw->nrxx, ucell.omega,
+                             this->chr.rhopw->nxyz, inp.nelec); // check the rho
 
 }
 
@@ -97,7 +130,8 @@ void ESolver_KS::hamilt2rho(UnitCell& ucell, const int istep, const int iter, co
     // example wavefunctions uses 20 processors while density uses 10.
     if (PARAM.globalv.ks_run)
     {
-        drho = p_chgmix->get_drho(&this->chr, this->inp_->nelec);
+        drho = module_charge::cal_drho(&this->chr, this->inp_->nelec, *this->pw_rho,
+                                       p_chgmix->get_mixing_config(), ucell.omega, ucell.tpiba);
         hsolver_error = 0.0;
         if (iter == 1 && this->inp_->calculation != "nscf")
         {
@@ -114,7 +148,8 @@ void ESolver_KS::hamilt2rho(UnitCell& ucell, const int istep, const int iter, co
 
                 this->hamilt2rho_single(ucell, istep, iter, diag_ethr);
 
-                drho = p_chgmix->get_drho(&this->chr, this->inp_->nelec);
+                drho = module_charge::cal_drho(&this->chr, this->inp_->nelec, *this->pw_rho,
+                                               p_chgmix->get_mixing_config(), ucell.omega, ucell.tpiba);
 
                 hsolver_error = hsolver::cal_hsolve_error(this->inp_->basis_type,
                                 this->inp_->esolver_type, diag_ethr, this->inp_->nelec);
@@ -255,9 +290,20 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &
     }
 #endif
 
-    module_charge::chgmixing_ks(iter, ucell, this->pelec, this->chr, this->p_chgmix, 
-      this->pw_rhod->nrxx, this->drho, this->oscillate_esolver, conv_esolver, hsolver_error, 
-      this->scf_thr, this->scf_ene_thr, converged_u, *this->inp_);
+    module_charge::ScfMixingCtx ctx;
+    ctx.hsolver_error = hsolver_error;
+    ctx.scf_thr = this->scf_thr;
+    ctx.scf_ene_thr = this->scf_ene_thr;
+    ctx.converged_u = converged_u;
+    ctx.ks_run = PARAM.globalv.ks_run;
+    ctx.drho = this->drho;
+    ctx.oscillate_esolver = this->oscillate_esolver;
+    ctx.conv_esolver = conv_esolver;
+    module_charge::chgmixing_ks(iter, ucell, this->pelec, this->chr,
+        *this->chr.rhopw, this->p_chgmix, ctx, *this->inp_);
+    this->drho = ctx.drho;
+    this->oscillate_esolver = ctx.oscillate_esolver;
+    conv_esolver = ctx.conv_esolver;
 
     // 2.3) Update potentials (should be done every SF iter)
     elecstate::update_pot(ucell, this->pelec, this->chr, conv_esolver);
@@ -277,7 +323,8 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &
     double dkin = 0.0; // for meta-GGA
     if (XC_Functional::get_ked_flag())
     {
-        dkin = p_chgmix->get_dkin(&this->chr, this->inp_->nelec);
+        dkin = module_charge::cal_dkin(&this->chr, this->inp_->nelec, *this->pw_rho,
+                                       p_chgmix->get_mixing_config(), ucell.omega);
     }
 
     // Iter finish