Refactor charge density and density matrix modules

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 ######
 
 #==========================

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\
@@ -823,18 +824,21 @@ 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\
+    chg_atomic.o\
     charge_init.o\
-    charge_mpi.o\
+    chg_tools.o\
+    chg_parallel.o\
     charge_extra.o\
     charge_mixing.o\
     charge_mixing_dmr.o\
-    mix_resid.o\
-    mix_precond.o\
+    chg_drho.o\
+    chg_precond.o\
     charge_mixing_rho.o\
-    charge_mixing_uspp.o\
+    chg_symm.o\
+    chg_symm_detail.o\
+    chg_uspp.o\
     rhog_io.o\
     write_elecstat_pot.o\
     write_init.o\
@@ -892,8 +896,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\

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;

source/source_esolver/esolver_fp.cpp

@@ -2,7 +2,7 @@
 
 #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_symm.h"
 #include "source_cell/read_pp_ucell.h"
 #include "source_estate/param_update.h"
 #include "source_hamilt/module_ewald/h_ewald_pw.h"

source/source_esolver/esolver_ks.cpp

@@ -5,7 +5,8 @@
 #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_pw/module_pwdft/setup_pwwfc.h" // mohan add 20251018
 #include "source_hsolver/hsolver.h"
 #include "source_io/module_energy/write_eig_occ.h"
@@ -65,9 +66,24 @@ 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;
+    mix_cfg.mixing_tau = inp.mixing_tau;
+    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;
+    p_chgmix->set_mixing(mix_cfg, ucell.omega, ucell.tpiba);
     p_chgmix->init_mixing();
 
     //! 4) setup plane wave for electronic wave functions
@@ -97,7 +113,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 +131,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);
@@ -273,7 +291,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 

source/source_esolver/esolver_ks_lcao.cpp

@@ -7,7 +7,7 @@
 #include "source_lcao/setup_dftu_lcao.h"
 #include "source_pw/module_pwdft/dftu_base.h" // Plus_U_Base (PW and LCAO share it)
 #include "source_hamilt/hs_matrix_k.h"
-#include "source_estate/module_charge/symm_rho.h"
+#include "source_estate/module_charge/chg_symm.h"
 #include "source_lcao/lcao_domain.h" // need DeePKS_init
 #include "source_lcao/force_stress_lcao.h"
 #include "source_hamilt/module_gint/gint.h"
@@ -18,7 +18,7 @@
 #include "../source_lcao/module_ri/exx_opt_orb.h"
 #endif
 #include "source_lcao/module_rdmft/rdmft.h"
-#include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006
+#include "source_estate/module_charge/chg_routine.h" // use charge mixing, mohan add 20251006
 #include "source_estate/module_dm/init_dm.h" // init dm from electronic wave functions
 #include "source_io/module_restart/restart.h" // GlobalC::restart for load_exx_flag
 #include "source_io/module_ctrl/ctrl_runner_lcao.h" // use ctrl_runner_lcao() 
@@ -234,7 +234,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
 #endif
 
     // 16) the electron charge density should be symmetrized,
-    Symmetry_rho::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rho, ucell.symm);
+    module_charge::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rho, ucell.symm);
 
     // 17) update of RDMFT, added by jghan
     if (this->inp_->rdmft == true)
@@ -506,7 +506,7 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2rho_single(UnitCell& ucell, int istep, int
 #endif
 
     // 5) symmetrize the charge density
-    Symmetry_rho::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rho, ucell.symm);
+    module_charge::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rho, ucell.symm);
 
     // 6) calculate delta energy
     this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);

source/source_esolver/esolver_ks_lcao.h

@@ -6,7 +6,7 @@
 #include "source_lcao/record_adj.h" // adjacent atoms
 #include "source_basis/module_nao/two_center_bundle.h" // nao basis
 #include "source_hamilt/module_gint/gint_info.h"
-#include "source_estate/module_charge/gint_prec_ctrl.h"
+#include "source_hamilt/module_gint/gint_prec_ctrl.h"
 #include "source_lcao/setup_deepks.h" // for deepks, mohan add 20251008
 #include "source_lcao/setup_exx.h" // for exx, mohan add 20251008
 #include "source_lcao/module_rdmft/rdmft.h" // rdmft

source/source_esolver/esolver_ks_lcao_tddft.cpp

@@ -13,7 +13,7 @@
 #include "source_io/module_wf/read_wfc_nao.h"
 //------LCAO HSolver ElecState-------
 #include "source_estate/elecstate_tools.h"
-#include "source_estate/module_charge/symm_rho.h"
+#include "source_estate/module_charge/chg_symm.h"
 #include "source_estate/module_dm/cal_dm_psi.h"
 #include "source_estate/module_dm/cal_edm_tddft.h"
 #include "source_estate/module_pot/h_tddft_pw.h"
@@ -378,7 +378,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::hamilt2rho_single(UnitCell& ucell, const
     // Symmetrize the charge density only for ground state
     if (istep <= 1)
     {
-        Symmetry_rho::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rho, ucell.symm);
+        module_charge::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rho, ucell.symm);
     }
 #ifdef __EXX
     if (this->exx_info_.info_ri.real_number)

source/source_esolver/esolver_ks_lcaopw.cpp

@@ -4,7 +4,7 @@
 #include <iostream>
 
 //--------------temporary----------------------------
-#include "source_estate/module_charge/symm_rho.h"
+#include "source_estate/module_charge/chg_symm.h"
 #include "source_estate/occupy.h"
 #include "source_hamilt/module_ewald/h_ewald_pw.h"
 //-----force-------------------
@@ -176,7 +176,7 @@ namespace ModuleESolver
         }
 #endif
 
-        Symmetry_rho::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rhod, ucell.symm);
+        module_charge::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rhod, ucell.symm);
 
         // deband is calculated from "output" charge density calculated
         // in sum_band

source/source_esolver/esolver_ks_pw.cpp

@@ -1,7 +1,7 @@
 #include "esolver_ks_pw.h"
 
 #include "source_estate/elecstate_pw.h"
-#include "source_estate/module_charge/symm_rho.h"
+#include "source_estate/module_charge/chg_symm.h"
 #include "source_hsolver/diago_iter_assist.h"
 #include "source_hsolver/diago_params.h"
 #include "source_hsolver/hsolver_pw.h"
@@ -14,7 +14,7 @@
 #include "source_base/kernels/dsp/dsp_connector.h"
 #endif
 
-#include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006
+#include "source_estate/module_charge/chg_routine.h" // use charge mixing, mohan add 20251006
 #include "source_estate/setup_estate_pw.h"         // mohan add 20251005
 #include "source_hamilt/module_xc/general_exx_info.h" // for General_Exx_Info type used via general_exx_info_
 #include "source_io/module_ctrl/ctrl_output_pw.h"  // mohan add 20250927
@@ -277,7 +277,7 @@ void ESolver_KS_PW<T, Device>::hamilt2rho_single(UnitCell& ucell, const int iste
     }
 
     // symmetrize the charge density
-    Symmetry_rho::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rhod, ucell.symm);
+    module_charge::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rhod, ucell.symm);
 
     ModuleBase::timer::end("ESolver_KS_PW", "hamilt2rho_single");
 }

source/source_esolver/esolver_of.cpp

@@ -3,7 +3,7 @@
 #include "source_io/module_parameter/parameter.h"
 //-----------temporary-------------------------
 #include "source_base/global_function.h"
-#include "source_estate/module_charge/symm_rho.h"
+#include "source_estate/module_charge/chg_symm.h"
 #include "source_hamilt/module_ewald/h_ewald_pw.h"
 #include "source_cell/cal_ux.h"
 #include "source_pw/module_pwdft/force_pw.h"
@@ -240,7 +240,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
     const int nspin = this->inp_->nspin;
     if (this->inp_->init_chg == "file")
     {
-        Symmetry_rho::symmetrize_rho(nspin, this->chr, this->pw_rho, ucell.symm);
+        module_charge::symmetrize_rho(nspin, this->chr, this->pw_rho, ucell.symm);
         for (int is = 0; is < nspin; ++is)
         {
             for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs)
@@ -385,19 +385,6 @@ void ESolver_OF::update_rho()
         }
     }
     // // ------------ turn on symmetry may cause instability in optimization ------------
-    // if (ModuleSymmetry::Symmetry::symm_flag == 1)
-    // {
-    //     Symmetry_rho srho;
-    //     for (int is = 0; is < this->inp_->nspin; is++)
-    //     {
-    //         srho.begin(is, *(this->chr), this->pw_rho, Pgrid, ucell.symm);
-    //         for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs)
-    //         {
-    //             this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]);
-    //         }
-    //     }
-    // }
-    // // --------------------------------------------------------------------------------
 }
 
 /**

source/source_esolver/esolver_of_tddft.cpp

@@ -3,7 +3,7 @@
 #include "source_io/module_parameter/parameter.h"
 //-----------temporary-------------------------
 #include "source_base/global_function.h"
-#include "source_estate/module_charge/symm_rho.h"
+#include "source_estate/module_charge/chg_symm.h"
 #include "source_hamilt/module_ewald/h_ewald_pw.h"
 #include "source_cell/cal_ux.h"
 //-----force-------------------

source/source_esolver/esolver_sdft_pw.cpp

@@ -3,7 +3,7 @@
 #include "source_base/global_variable.h"
 #include "source_base/memory_recorder.h"
 #include "source_base/parallel_comm.h"
-#include "source_estate/module_charge/symm_rho.h"
+#include "source_estate/module_charge/chg_symm.h"
 #include "source_hsolver/diago_iter_assist.h"
 #include "source_hsolver/diago_params.h"
 #include "source_io/module_parameter/parameter.h"
@@ -193,7 +193,7 @@ void ESolver_SDFT_PW<T, Device>::hamilt2rho_single(UnitCell& ucell, int istep, i
 
     if (PARAM.globalv.ks_run)
     {
-        Symmetry_rho::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rho, ucell.symm);
+        module_charge::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rho, ucell.symm);
         this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);
     }
     else

source/source_esolver/lcao_others.cpp

@@ -5,7 +5,7 @@
 #include "source_cell/module_neighbor/sltk_grid_driver.h"
 #include "source_esolver/esolver_ks_lcao.h"
 #include "source_estate/elecstate_lcao.h"
-#include "source_estate/module_charge/symm_rho.h"
+#include "source_estate/module_charge/chg_symm.h"
 #include "source_estate/module_dm/cal_dm_psi.h"
 #include "source_hamilt/module_gint/gint.h"
 #include "source_io/module_chgpot/get_pchg_lcao.h"

source/source_estate/CMakeLists.txt

@@ -27,20 +27,21 @@ list(APPEND objects
     module_pot/td_field.cpp
     module_pot/td_field_profiles.cpp
     module_pot/td_field_manager.cpp
-    module_charge/chgmixing.cpp
-    module_charge/gint_prec_ctrl.cpp
+    module_charge/chg_routine.cpp
     module_charge/charge.cpp
+    module_charge/chg_atomic.cpp
     module_charge/charge_init.cpp
-    module_charge/charge_mpi.cpp
+    module_charge/chg_tools.cpp
+    module_charge/chg_parallel.cpp
     module_charge/charge_extra.cpp
     module_charge/charge_mixing.cpp
     module_charge/charge_mixing_dmr.cpp
-    module_charge/mix_resid.cpp
-    module_charge/mix_precond.cpp
+    module_charge/chg_drho.cpp
+    module_charge/chg_precond.cpp
     module_charge/charge_mixing_rho.cpp
-    module_charge/charge_mixing_uspp.cpp
-    module_charge/symm_rho.cpp
-    module_charge/symm_rhog.cpp
+    module_charge/chg_symm.cpp
+    module_charge/chg_uspp.cpp
+    module_charge/chg_symm_detail.cpp
     rhog_io.cpp
     write_elecstat_pot.cpp
     write_init.cpp

source/source_estate/elecstate_pw.cpp

@@ -5,6 +5,7 @@
 #include "source_base/module_device/device.h"
 #include "source_base/parallel_reduce.h"
 #include "source_base/timer.h"
+#include "source_estate/module_charge/chg_parallel.h"
 #include "source_estate/uspp_density.h"
 #include "source_hamilt/module_xc/xc_functional.h"
 #include "source_io/module_parameter/parameter.h"
@@ -171,7 +172,7 @@ template <typename T, typename Device>
 void ElecStatePW<T, Device>::parallelK()
 {
 #ifdef __MPI
-    this->charge->rho_mpi();
+    module_charge::rho_mpi(*this->charge);
 #endif
 }
 

source/source_estate/elecstate_pw_cal_tau.cpp

@@ -1,5 +1,7 @@
 #include "elecstate_pw.h"
 
+#include "source_estate/module_charge/chg_parallel.h"
+
 namespace elecstate {
 
 template<typename T, typename Device>
@@ -55,7 +57,7 @@ void ElecStatePW<T, Device>::cal_tau(const psi::Psi<T, Device>& psi)
         }
     }
 #ifdef __MPI
-    this->charge->kin_r_mpi();
+    module_charge::kin_r_mpi(*this->charge);
 #endif
     ModuleBase::TITLE("ElecStatePW", "cal_tau");
 }