Remove GlobalC::dftu (#6694)

* add update_pot in source_estate

* split update_pot in rt-TDDFT into two functions, one is the original update_pot, which is replaced by update_pot function in estate, the other is save2 function used in after_iter

* fix bugs

* update

* update esolver_ks

* small updates

* small bug fixed

* update E_bandgap to E_gap(k)

* delete rho_restart file

* change setup_parameters to print_parameters

* move setup_pw.cpp to setup_pwwfc.cpp in module_pwdft

* update esolver, small things

* delete some old_gint codes in esolver

* fix bugs, update outputs

* fix bug in tests

* change efermi to Efermi in fp_energy.h

* fix bugs

* fix exd and exc in ctrl_runner_lcao

* fix cmake

* fix bug

* delete lcao_after_scf and lcao_before_scf

* modify codes including read_wfc_nao to delete pelec, add setup_psi_lcao codes

* fix bugs

* delete grid integral in FORCE_STRESS

* move read psi to esolver, not in psi

* fix system bug

* fix bug in read_wfc_nao.h

* change PARAM.inp to inp

* fix bugs

* delete omega in pelec, remove some deepks codes to setup_deepks.cpp

* remove exx redundant parameters in Hamilt_LCAO

* Fix circular dependencies of header files.

* fix bugs in MLALGO

* refactor exx lcao

* add rho_tau_lcao files

* fix bug in setup_deepks

* move rho_tau_lcao to source_lcao

* remove old psiToRho function in elecstate_lcao and replace with new rho_tau_lcao

* update estate

* fix bug in ELF

* fix bug

* move DM outside of pelec

* deal with init_dm functions

* remove pelec in dftu, replace with dmk_d(ouble) and dmk_c(omplex)

* fix cal_ldos with DM

* update delta_spin

* update delta spin run_lambda_loop

* the big issue to remove dm may be caused by cal_mw_from_lambda

* remove some DM

* keep updating DM

* refactor some DM in double_xc

* now all functions can be compiled

* fix it

* update again

* add and compile setup_dm

* update fix bugs

* simplify dm

* one more fix

* now I can compile the code

* fix some bugs in EXX version

* fix FORCE_STRESS

* fix 2

* fix 3

* fix 4

* fix 5

* improve esolver_ks_lcao

* improve again

* comment out psi2rho in init_dm

* clean esolver_ks_lcao

* update structure_factor

* update sf

* update

* fix sf

* fix

* fix structure factor

* update DFTU, remove GlobalC of DFTU

* remove some GlobalC

* remove GlobalC, rename ModuleDFTU::DFTU class to Plus_U

* keep updating DFTU

* fix +U forces and stress

* fix dm bug in PEXSI

* fix PEXSI

* Update edm.cpp

* fix edm

* Update hsolver_lcao.cpp

* Update hsolver_lcao.cpp

* keep removing dftu

* update dftu, remove GlobalC

* keep removing GlobalC::DFTU

* remove all GlobalC::dftu

* fix gpu +U

* fix tests related to removal of GlobalC::dftu

* fix again

* fix bug that ignoring pass dftu pointer

* update +U

* add output of DFT+U energy term

---------

Co-authored-by: Levi Zhou <[email protected]>

Author: mohanchen

source/source_esolver/esolver_double_xc.cpp

@@ -146,7 +146,8 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
             this->kv,
             this->two_center_bundle_,
             this->orb_,
-            this->dmat_base.dm,
+			this->dmat_base.dm,
+			&this->dftu,
 			this->deepks,
 			istep,
 			this->exx_nao);
@@ -396,6 +397,7 @@ void ESolver_DoubleXC<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
                        this->kv,
                        this->pw_rho,
                        this->solvent,
+                       this->dftu,
                        this->deepks,
 					   this->exx_nao,
 					   &ucell.symm);

source/source_esolver/esolver_gets.h

@@ -37,7 +37,7 @@ class ESolver_GetS : public ESolver_KS<std::complex<double>>
 
     TwoCenterBundle two_center_bundle_;
 
-    // temporary introduced during removing GlobalC::ORB
+    // temporary introduced
     LCAO_Orbitals orb_;
 };
 } // namespace ModuleESolver

source/source_esolver/esolver_ks.cpp

@@ -15,6 +15,7 @@
 #include "source_estate/elecstate_print.h" // print_etot
 #include "source_io/print_info.h" // print_parameters
 #include "source_psi/setup_psi.h" // mohan add 20251009
+#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-07
 
 namespace ModuleESolver
 {
@@ -281,9 +282,19 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
                                        this->pelec->nelec_spin.data());
 
     // 2.2) charge mixing 
+    // SCF will continue if U is not converged for uramping calculation
+	bool converged_u = true;
+	// to avoid unnecessary dependence on dft+u, refactor is needed
+#ifdef __LCAO
+	if (PARAM.inp.dft_plus_u)
+	{
+		converged_u = this->dftu.u_converged();
+	}
+#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, PARAM.inp);
+      this->scf_thr, this->scf_ene_thr, converged_u, PARAM.inp);
 
     // 2.3) Update potentials (should be done every SF iter)
     elecstate::update_pot(ucell, this->pelec, this->chr, conv_esolver);

source/source_esolver/esolver_ks.h

@@ -7,6 +7,7 @@
 #include "source_estate/module_charge/charge_mixing.h" // use charge mixing
 #include "source_psi/psi.h" // use electronic wave functions
 #include "source_hamilt/hamilt.h" // use Hamiltonian
+#include "source_lcao/module_dftu/dftu.h" // mohan add 20251107
 
 namespace ModuleESolver
 {
@@ -61,6 +62,9 @@ class ESolver_KS : public ESolver_FP
     //! Electronic wavefunctions
     psi::Psi<T>* psi = nullptr;
 
+    //! DFT+U method, mohan add 2025-11-07
+    Plus_U dftu;
+
     std::string basisname;      //! esolver_ks_lcao.cpp
     double esolver_KS_ne = 0.0; //! number of electrons
     double diag_ethr;           //! the threshold for diagonalization

source/source_esolver/esolver_ks_lcao.cpp

@@ -4,7 +4,6 @@
 #include "source_lcao/hs_matrix_k.hpp" // there may be multiple definitions if using hpp
 #include "source_estate/module_charge/symmetry_rho.h"
 #include "source_lcao/LCAO_domain.h" // need DeePKS_init
-#include "source_lcao/module_dftu/dftu.h"
 #include "source_lcao/FORCE_STRESS.h"
 #include "source_estate/elecstate_lcao.h"
 #include "source_lcao/hamilt_lcao.h"
@@ -99,8 +98,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     // 9) initialize DFT+U
     if (inp.dft_plus_u)
     {
-        auto* dftu = ModuleDFTU::DFTU::get_instance();
-        dftu->init(ucell, &this->pv, this->kv.get_nks(), &orb_);
+        this->dftu.init(ucell, &this->pv, this->kv.get_nks(), &orb_);
     }
 
     // 10) init local pseudopotentials
@@ -189,7 +187,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     {
         this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
             ucell, this->gd, &this->pv, this->pelec->pot, this->kv,
-            two_center_bundle_, orb_, this->dmat.dm, this->deepks, istep, exx_nao);
+            two_center_bundle_, orb_, this->dmat.dm, &this->dftu, this->deepks, istep, exx_nao);
     }
 
     // 9) for each ionic step, the overlap <phi|alpha> must be rebuilt
@@ -274,7 +272,7 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
                        this->gd, this->pv, this->pelec, this->dmat, this->psi,
                        two_center_bundle_, orb_, force, this->scs,
                        this->locpp, this->sf, this->kv,
-                       this->pw_rho, this->solvent, this->deepks,
+                       this->pw_rho, this->solvent, this->dftu, this->deepks,
                        this->exx_nao, &ucell.symm);
 
     // delete RA after cal_force
@@ -338,7 +336,8 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     // call iter_init() of ESolver_KS
     ESolver_KS<TK>::iter_init(ucell, istep, iter);
 
-    module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, this->dmat.dm->get_DMR_pointer(1)->get_nnr(), PARAM.inp); 
+    module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, this->dftu, 
+      this->dmat.dm->get_DMR_pointer(1)->get_nnr(), PARAM.inp); 
 
     // mohan update 2012-06-05
     this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(ucell);
@@ -377,10 +376,10 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     {
         if (istep != 0 || iter != 1)
         {
-            GlobalC::dftu.set_dmr(this->dmat.dm);
+            this->dftu.set_dmr(this->dmat.dm);
         }
         // Calculate U and J if Yukawa potential is used
-        GlobalC::dftu.cal_slater_UJ(ucell, this->chr.rho, this->pw_rho->nrxx);
+        this->dftu.cal_slater_UJ(ucell, this->chr.rho, this->pw_rho->nrxx);
     }
 
 #ifdef __MLALGO
@@ -489,18 +488,18 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
     // 1) calculate the local occupation number matrix and energy correction in DFT+U
     if (PARAM.inp.dft_plus_u)
     {
-        // only old DFT+U method should calculated energy correction in esolver,
-        // new DFT+U method will calculate energy in calculating Hamiltonian
+        // old DFT+U method calculates energy correction in esolver,
+        // new DFT+U method calculates energy in Hamiltonian
         if (PARAM.inp.dft_plus_u == 2)
         {
-            if (GlobalC::dftu.omc != 2)
+            if (this->dftu.omc != 2)
             {
-                ModuleDFTU::dftu_cal_occup_m(iter, ucell, dm_vec, this->kv,
-                  this->p_chgmix->get_mixing_beta(), hamilt_lcao);
+                dftu_cal_occup_m(iter, ucell, dm_vec, this->kv,
+                  this->p_chgmix->get_mixing_beta(), hamilt_lcao, this->dftu);
             }
-            GlobalC::dftu.cal_energy_correction(ucell, istep);
+            this->dftu.cal_energy_correction(ucell, istep);
         }
-        GlobalC::dftu.output(ucell);
+        this->dftu.output(ucell);
     }
 
     // 2) for deepks, calculate delta_e, output labels during electronic steps
@@ -532,7 +531,7 @@ 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)
     {
-        GlobalC::dftu.initialed_locale = true;
+        this->dftu.initialed_locale = true;
     }
 
     // control the output related to the finished iteration
@@ -566,7 +565,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
     //! 2) output of lcao every few ionic steps
     ModuleIO::ctrl_scf_lcao<TK, TR>(ucell,
             PARAM.inp, this->kv, this->pelec, this->dmat.dm, this->pv,
-            this->gd, this->psi, hamilt_lcao, this->two_center_bundle_,
+            this->gd, this->psi, hamilt_lcao, this->dftu, this->two_center_bundle_,
             this->orb_, this->pw_wfc, this->pw_rho, this->pw_big, this->sf,
             this->rdmft_solver, this->deepks, this->exx_nao, 
             this->conv_esolver, this->scf_nmax_flag, istep);

source/source_esolver/esolver_ks_lcao.h

@@ -69,32 +69,34 @@ class ESolver_KS_LCAO : public ESolver_KS<TK>
     //! GintInfo: used to store some basic infomation about module_gint
     std::unique_ptr<ModuleGint::GintInfo> gint_info_;
 
+    //! NAO: store related information 
+    LCAO_Orbitals orb_;
+
     //! NAO orbitals: two-center integrations
     TwoCenterBundle two_center_bundle_;
 
     //! Add density matrix class, mohan add 2025-10-30
     LCAO_domain::Setup_DM<TK> dmat;
 
+
+    // For deepks method, mohan add 2025-10-08
+    Setup_DeePKS<TK> deepks;
+
+    // For exact-exchange energy, mohan add 2025-10-08
+    Exx_NAO<TK> exx_nao;
+
     //! For RDMFT calculations, added by jghan, 2024-03-16 
     rdmft::RDMFT<TK, TR> rdmft_solver;
 
-    //! NAO: store related information 
-    LCAO_Orbitals orb_;
+    //! For linear-response TDDFT
+    friend class LR::ESolver_LR<double, double>;
+    friend class LR::ESolver_LR<std::complex<double>, double>;
 
     // Temporarily store the stress to unify the interface with PW,
     // because it's hard to seperate force and stress calculation in LCAO.
     ModuleBase::matrix scs;
     bool have_force = false;
 
-    // deepks method, mohan add 2025-10-08
-    Setup_DeePKS<TK> deepks;
-
-    // exact-exchange energy, mohan add 2025-10-08
-    Exx_NAO<TK> exx_nao;
-
-    friend class LR::ESolver_LR<double, double>;
-    friend class LR::ESolver_LR<std::complex<double>, double>;
-
 
   public:
     const Record_adj & get_RA() const { return RA; }

source/source_esolver/esolver_ks_pw.cpp

@@ -55,7 +55,13 @@ ESolver_KS_PW<T, Device>::~ESolver_KS_PW()
 template <typename T, typename Device>
 void ESolver_KS_PW<T, Device>::allocate_hamilt(const UnitCell& ucell)
 {
-    this->p_hamilt = new hamilt::HamiltPW<T, Device>(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell, &ucell);
+	this->p_hamilt = new hamilt::HamiltPW<T, Device>(
+			this->pelec->pot, 
+			this->pw_wfc, 
+			&this->kv, 
+			&this->ppcell, 
+			&this->dftu,
+			&ucell);
 }
 
 template <typename T, typename Device>
@@ -133,8 +139,10 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
     this->allocate_hamilt(ucell);
 
     //! Setup potentials (local, non-local, sc, +U, DFT-1/2)
+    // note: init DFT+U is done here for pw basis for every scf iteration, however, 
+    // init DFT+U is done in "before_all_runners" in LCAO basis. This should be refactored, mohan note 2025-11-06
     pw::setup_pot(istep, ucell, this->kv, this->sf, this->pelec, this->Pgrid,
-              this->chr, this->locpp, this->ppcell, this->vsep_cell,
+              this->chr, this->locpp, this->ppcell, this->dftu, this->vsep_cell,
               this->stp.psi_t, this->p_hamilt, this->pw_wfc, this->pw_rhod, PARAM.inp);
 
     // setup psi (electronic wave functions)
@@ -162,7 +170,7 @@ void ESolver_KS_PW<T, Device>::iter_init(UnitCell& ucell, const int istep, const
     ESolver_KS<T, Device>::iter_init(ucell, istep, iter);
 
     // 2) perform charge mixing for KSDFT using pw basis
-    module_charge::chgmixing_ks_pw(iter, this->p_chgmix, PARAM.inp);
+    module_charge::chgmixing_ks_pw(iter, this->p_chgmix, this->dftu, PARAM.inp);
 
     // 3) mohan move harris functional here, 2012-06-05
     // use 'rho(in)' and 'v_h and v_xc'(in)
@@ -172,14 +180,13 @@ void ESolver_KS_PW<T, Device>::iter_init(UnitCell& ucell, const int istep, const
     // should before lambda loop in DeltaSpin
     if (PARAM.inp.dft_plus_u && (iter != 1 || istep != 0))
     {
-        auto* dftu = ModuleDFTU::DFTU::get_instance();
         // only old DFT+U method should calculate energy correction in esolver,
         // new DFT+U method will calculate energy when evaluating the Hamiltonian
-        if (dftu->omc != 2)
+        if (this->dftu.omc != 2)
         {
-            dftu->cal_occ_pw(iter, this->stp.psi_t, this->pelec->wg, ucell, PARAM.inp.mixing_beta);
+            this->dftu.cal_occ_pw(iter, this->stp.psi_t, this->pelec->wg, ucell, PARAM.inp.mixing_beta);
         }
-        dftu->output(ucell);
+        this->dftu.output(ucell);
     }
 }
 
@@ -389,7 +396,7 @@ void ESolver_KS_PW<T, Device>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
 
     // Calculate forces
     ff.cal_force(ucell, force, *this->pelec, this->pw_rhod, &ucell.symm,
-                 &this->sf, this->solvent, &this->locpp, &this->ppcell, 
+                 &this->sf, this->solvent, &this->dftu, &this->locpp, &this->ppcell, 
                  &this->kv, this->pw_wfc, this->stp.psi_d);
 }
 
@@ -401,7 +408,7 @@ void ESolver_KS_PW<T, Device>::cal_stress(UnitCell& ucell, ModuleBase::matrix& s
     // mohan add 2025-10-12
     this->stp.update_psi_d();
 
-    ss.cal_stress(stress, ucell, this->locpp, this->ppcell, this->pw_rhod,
+    ss.cal_stress(stress, ucell, this->dftu, this->locpp, this->ppcell, this->pw_rhod,
                   &ucell.symm, &this->sf, &this->kv, this->pw_wfc, this->stp.psi_d);
 
     // external stress

source/source_esolver/esolver_of.cpp

@@ -575,7 +575,10 @@ double ESolver_OF::cal_energy()
 void ESolver_OF::cal_force(UnitCell& ucell, ModuleBase::matrix& force)
 {
     Forces<double> ff(ucell.nat);
-    ff.cal_force(ucell, force, *pelec, this->pw_rho, &ucell.symm, &sf, this->solvent, &this->locpp);
+ 
+    // here nullptr is for DFT+U, which may cause bugs, mohan note 2025-11-07
+    // solvent can be used? mohan ask 2025-11-07
+    ff.cal_force(ucell, force, *pelec, this->pw_rho, &ucell.symm, &sf, this->solvent, nullptr, &this->locpp);
 }
 
 /**

source/source_esolver/lcao_others.cpp

@@ -182,6 +182,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
             two_center_bundle_,
             orb_,
             this->dmat.dm,
+            &this->dftu,
 			this->deepks,
 			istep,
 			this->exx_nao);

source/source_estate/elecstate_energy_terms.cpp

@@ -4,7 +4,7 @@
 #include "source_estate/module_pot/gatefield.h"
 #include "source_lcao/module_deepks/LCAO_deepks.h"
 #include "source_lcao/module_deltaspin/spin_constrain.h"
-#include "source_lcao/module_dftu/dftu.h"
+#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-06
 
 namespace elecstate
 {
@@ -36,7 +36,7 @@ double ElecState::get_solvent_model_Acav()
 
 double ElecState::get_dftu_energy()
 {
-    return GlobalC::dftu.get_energy();
+    return Plus_U::get_energy();
 }
 
 double ElecState::get_local_pp_energy()
@@ -52,4 +52,4 @@ double ElecState::get_local_pp_energy()
     return local_pseudopot_energy;
 }
 
-} // namespace elecstate
+} // namespace elecstate