Refactor DFT+U codes

source/Makefile.Objects

@@ -741,6 +741,9 @@ OBJS_SRCPW=h_ewald_pw.o\
     setup_pwrho.o\
     setup_pwwfc.o\
     update_cell_pw.o\
+    dftu_base.o\
+    dftu_output.o\
+    dftu_pw.o\
     setup_dftu_pw.o\
     deltaspin_pw.o\
     force_pw.o\
@@ -809,7 +812,6 @@ OBJS_DFTU=dftu.o\
       dftu_force.o\
       dftu_yukawa.o\
       dftu_folding.o\
-      dftu_io.o\
       dftu_tools.o\
       dftu_occup.o\
       dftu_hamilt.o\

source/source_lcao/module_dftu/CMakeLists.txt

@@ -3,11 +3,9 @@ list(APPEND objects
       dftu_force.cpp
       dftu_yukawa.cpp
       dftu_folding.cpp
-      dftu_io.cpp
       dftu_tools.cpp
       dftu_occup.cpp
       dftu_hamilt.cpp
-      dftu_pw.cpp
 )
 
 add_library(

source/source_lcao/module_dftu/dftu.cpp

@@ -21,22 +21,8 @@
 #include <vector>
 
  // mohan add 2025-11-06
-double Plus_U::energy_u = 0.0;
-
-std::vector<double> Plus_U::U = {}; // U (Hubbard parameter U)
-
-std::vector<double> Plus_U::U0 = {}; // U0 (target Hubbard parameter U0)
-
-std::vector<int> Plus_U::orbital_corr = {}; //
-
-double Plus_U::uramping = 0.0; // increase U by uramping, default is -1.0
-
-int Plus_U::omc=0; // occupation matrix control
-
-int Plus_U::mixing_dftu=0; //whether to mix locale
-int Plus_U::nspin=0;
-
-bool Plus_U::Yukawa=false; // whether to use Yukawa potential
+// Static member definitions moved to dftu_base.cpp (Plus_U_Base::)
+// Plus_U inherits these from Plus_U_Base.
 
 Plus_U::Plus_U()
 {}
@@ -68,14 +54,9 @@ void Plus_U::init(UnitCell& cell,
 {
     ModuleBase::TITLE("Plus_U", "init");
 
-#ifndef __MPI
-    std::cout << "DFT+U module is only accessible in mpi versioin" << std::endl;
-    exit(0);
-#endif
-
     this->paraV = pv;
 
-#ifdef __LCAO    
+#ifdef __LCAO
     ptr_orb_ = orb;
     if(ptr_orb_ != nullptr)
     {
@@ -84,16 +65,6 @@ void Plus_U::init(UnitCell& cell,
     ucell = &cell;
 #endif
 
-    Plus_U::nspin = nspin;
-    Plus_U::orbital_corr = orbital_corr;
-    Plus_U::Yukawa = yukawa_potential;
-    this->yukawa_lambda = yukawa_lambda;
-
-    this->global_readin_dir = global_readin_dir;
-    this->global_out_dir = global_out_dir;
-    this->init_chg = init_chg;
-    this->npol = npol;
-
     if (pv != nullptr)
     {
         const int global_rows = pv->get_global_row_size();
@@ -107,204 +78,23 @@ void Plus_U::init(UnitCell& cell,
             ModuleBase::WARNING_QUIT("Plus_U::init", "nlocal does not match global matrix dimension");
         }
     }
-    this->nlocal = nlocal;
-
-    this->gamma_only_local = gamma_only_local;
-    this->ks_solver = ks_solver;
-    this->cal_force = cal_force;
-    this->cal_stress = cal_stress;
-    this->device = device;
-    this->kpar = kpar;
-
-    // mohan update 2025-11-06
-    Plus_U::energy_u = 0.0;
-
-    this->locale.resize(cell.nat);
-    this->locale_save.resize(cell.nat);
-    // only for PW base
-    this->eff_pot_pw_index.resize(cell.nat);
-    int pot_index = 0;
-
-    this->iatlnmipol2iwt.resize(cell.nat);
-
-    int num_locale = 0;
-    // it:index of type of atom
-    for (int it = 0; it < cell.ntype; ++it)
-    {
-        for (int ia = 0; ia < cell.atoms[it].na; ia++)
-        {
-            // ia:index of atoms of this type
-            // determine the size of locale
-            const int iat = cell.itia2iat(it, ia);
-
-            locale[iat].resize(cell.atoms[it].nwl + 1);
-            locale_save[iat].resize(cell.atoms[it].nwl + 1);
-
-            // initialize the arrry iatlnm2iwt[iat][l][n][m]
-            this->iatlnmipol2iwt[iat].resize(cell.atoms[it].nwl + 1);
-
-            if(!has_correlated_orbital(it))
-            {
-                continue;
-            }
-
-            const int tlp1_npol = (get_orbital_corr(it)*2+1)*npol;
-            const int tlp1 = 2 * get_orbital_corr(it) + 1;
-            const int elem_size = tlp1 * tlp1;
-    // eff_pot_pw_index: per-atom offset into eff_pot_pw (and uom_array)
-    //
-    // nspin=1: offset = sum(tlp1^2 for preceding atoms), total = sum(all tlp1^2)
-    // nspin=2: same per-spin-channel offset; after the loop, pot_index *= 2
-    //          to create split layout: [all_spin_up | all_spin_down]
-    //          spin-up  at eff_pot_pw[eff_pot_pw_index[iat] + mm]
-    //          spin-down at eff_pot_pw[size/2 + eff_pot_pw_index[iat] + mm]
-    // nspin=4: offset = sum(tlp1_npol^2) where tlp1_npol = (2l+1)*npol = 2*(2l+1)
-    //          each atom occupies (2*tlp1)^2 = 4*tlp1^2 entries for 4 Pauli blocks
-            if(nspin == 4)
-            {
-                this->eff_pot_pw_index[iat] = pot_index;
-                pot_index += tlp1_npol * tlp1_npol;
-            }
-            else // nspin=1 or nspin=2: one tlp1^2 block per atom per spin channel
-            {
-                this->eff_pot_pw_index[iat] = pot_index;
-                pot_index += elem_size;
-            }
-
-            for (int l = 0; l <= cell.atoms[it].nwl; l++)
-            {
-                const int N = cell.atoms[it].l_nchi[l];
-
-                locale[iat][l].resize(N);
-                locale_save[iat][l].resize(N);
-
-                for (int n = 0; n < N; n++)
-                {
-                    if (nspin == 1 || nspin == 2)
-                    {
-                        locale[iat][l][n].resize(2);
-                        locale_save[iat][l][n].resize(2);
-
-                        locale[iat][l][n][0].create(2 * l + 1, 2 * l + 1);
-                        locale[iat][l][n][1].create(2 * l + 1, 2 * l + 1);
-
-                        locale_save[iat][l][n][0].create(2 * l + 1, 2 * l + 1);
-                        locale_save[iat][l][n][1].create(2 * l + 1, 2 * l + 1);
-                        num_locale += (2 * l + 1) * (2 * l + 1) * 2;
-                    }
-                    else if (nspin == 4) // SOC
-                    {
-                        locale[iat][l][n].resize(1);
-                        locale_save[iat][l][n].resize(1);
 
-                        locale[iat][l][n][0].create((2 * l + 1) * npol, (2 * l + 1) * npol);
-                        locale_save[iat][l][n][0].create((2 * l + 1) * npol, (2 * l + 1) * npol);
-                        num_locale += (2 * l + 1) * (2 * l + 1) * npol * npol;
-                    }
-                }
-            }
-
-            // initialize the arrry iatlnm2iwt[iat][l][n][m]
-            this->iatlnmipol2iwt[iat].resize(cell.atoms[it].nwl + 1);
-            for (int L = 0; L <= cell.atoms[it].nwl; L++)
-            {
-                this->iatlnmipol2iwt[iat][L].resize(cell.atoms[it].l_nchi[L]);
-
-                for (int n = 0; n < cell.atoms[it].l_nchi[L]; n++)
-                {
-                    this->iatlnmipol2iwt[iat][L][n].resize(2 * L + 1);
-
-                    for (int m = 0; m < 2 * L + 1; m++)
-                    {
-                        this->iatlnmipol2iwt[iat][L][n][m].resize(npol);
-                    }
-                }
-            }
-
-            for (int iw = 0; iw < cell.atoms[it].nw * npol; iw++)
-            {
-                int iw0 = iw / npol;
-                int ipol = iw % npol;
-                int iwt = cell.itiaiw2iwt(it, ia, iw);
-                int l = cell.atoms[it].iw2l[iw0];
-                int n = cell.atoms[it].iw2n[iw0];
-                int m = cell.atoms[it].iw2m[iw0];
-
-                this->iatlnmipol2iwt[iat][l][n][m][ipol] = iwt;
-            }
-        }
-    }
-    // allocate memory for eff_pot_pw
-    // nspin=2: split layout [all_spin_up | all_spin_down], double the size
-    // nspin=4: each atom already has 4*tlp1^2 (tlp1_npol^2) entries for Pauli blocks
-    if (nspin == 2) pot_index *= 2;
-
-    this->eff_pot_pw.resize(pot_index, 0.0);
-    this->uom_array.resize(pot_index, 0.0);
-    this->uom_save.resize(pot_index, 0.0);
-
-    if (Yukawa)
-    {
-        this->Fk.resize(cell.ntype);
-
-        this->U_Yukawa.resize(cell.ntype);
-        this->J_Yukawa.resize(cell.ntype);
-
-        for (int it = 0; it < cell.ntype; it++)
-        {
-            const int NL = cell.atoms[it].nwl + 1;
-
-            this->Fk[it].resize(NL);
-            this->U_Yukawa[it].resize(NL);
-            this->J_Yukawa[it].resize(NL);
-
-            for (int l = 0; l < NL; l++)
-            {
-                int N = cell.atoms[it].l_nchi[l];
-
-                this->Fk[it][l].resize(N);
-                for (int n = 0; n < N; n++)
-                {
-                    this->Fk[it][l][n].resize(l + 1, 0.0);
-                }
-
-                this->U_Yukawa[it][l].resize(N, 0.0);
-                this->J_Yukawa[it][l].resize(N, 0.0);
-            }
-        }
-    }
-
-    if (omc != 0)
-    {
-        std::stringstream sst;
-        sst << this->global_readin_dir << "dm_onsite_ini.txt";
-        this->read_occup_m(cell, sst.str(), this->init_chg, nspin, npol);
-#ifdef __MPI
-        this->local_occup_bcast(cell, nspin, npol);
-#endif
-
-        mark_locale_initialized();
-        this->copy_locale(cell);
-    }
-    else
-    {
-        if (this->init_chg == "file")
-        {
-            std::stringstream sst;
-            sst << this->global_readin_dir << "dm_onsite.txt";
-            this->read_occup_m(cell, sst.str(), this->init_chg, nspin, npol);
-#ifdef __MPI
-            this->local_occup_bcast(cell, nspin, npol);
-#endif
-            mark_locale_initialized();
-        }
-        else
-        {
-            this->zero_locale(cell);
-        }
-    }
-
-    ModuleBase::Memory::record("Plus_U::locale", sizeof(double) * num_locale);
+    this->init_base(cell,
+                    npol,
+                    nspin,
+                    orbital_corr,
+                    yukawa_potential,
+                    yukawa_lambda,
+                    global_readin_dir,
+                    global_out_dir,
+                    init_chg,
+                    nlocal,
+                    gamma_only_local,
+                    ks_solver,
+                    cal_force,
+                    cal_stress,
+                    device,
+                    kpar);
     return;
 }
 
@@ -465,45 +255,8 @@ void Plus_U::cal_energy_correction(const UnitCell& ucell,
 
 #endif
 
-void Plus_U::uramping_update()
-{
-    // Yukawa calculates U directly every iteration, no need for ramping
-    if (Yukawa) {
-        return;
-    }
-    // if uramping < 0.1, use the original U
-    if (this->uramping < 0.01) {
-        return;
-}
-    // loop to change U
-    for (int i = 0; i < this->U0.size(); i++)
-    {
-        if (this->U[i] + this->uramping < this->U0[i])
-        {
-            this->U[i] += this->uramping;
-        }
-        else
-        {
-            this->U[i] = this->U0[i];
-        }
-    }
-}
-
-bool Plus_U::u_converged()
-{
-    // Yukawa calculates U directly every iteration, always considered converged
-    if (Yukawa) {
-        return true;
-    }
-    for (int i = 0; i < this->U0.size(); i++)
-    {
-        if (this->U[i] != this->U0[i])
-        {
-            return false;
-        }
-    }
-    return true;
-}
+// uramping_update() and u_converged() are now implemented in
+// dftu_base.cpp as Plus_U_Base methods (inherited by Plus_U).
 
 #ifdef __LCAO