Phase A: extract output/write_occup_m as free functions in dftu_io namespace

- New files: source_pw/module_pwdft/dftu_output.h, dftu_output.cpp
  - declare and implement dftu_io::output and dftu_io::write_occup_m
  - first parameter is const Plus_U_Base&, accesses state via public getters
  - JacobiRotate and CalculateEigenvalues inline helpers migrated here
- dftu_base.h: added public getters get_U_Yukawa / get_J_Yukawa;
  removed output/write_occup_m member declarations
- dftu_base.cpp: removed output/write_occup_m implementations and the
  two Jacobi inline helpers (no longer referenced here)
- Callers updated:
  - setup_dftu_pw.cpp: dftu.output(...) -> dftu_io::output(dftu, ...)
  - setup_dftu_lcao.cpp: dftu_ptr->output(...) -> dftu_io::output(*dftu_ptr, ...)
- CMakeLists.txt and Makefile.Objects: added dftu_output.cpp/dftu_output.o

All 4 dftu unit tests pass.

Author: abacus_fixer

source/Makefile.Objects

@@ -742,6 +742,7 @@ OBJS_SRCPW=h_ewald_pw.o\
     setup_pwwfc.o\
     update_cell_pw.o\
     dftu_base.o\
+    dftu_output.o\
     dftu_pw.o\
     setup_dftu_pw.o\
     deltaspin_pw.o\

source/source_lcao/setup_dftu_lcao.cpp

@@ -1,5 +1,6 @@
 #include "setup_dftu_lcao.h"
 #include "source_lcao/module_dftu/dftu.h"
+#include "source_pw/module_pwdft/dftu_output.h" // mohan add 2025-11-08
 #include "source_estate/module_dm/density_matrix.h"
 #include "source_lcao/hamilt_lcao.h"
 
@@ -67,7 +68,7 @@ void finish_dftu_lcao(const int iter,
         }
         dftu_ptr->cal_energy_correction(ucell, iter);
     }
-    dftu_ptr->output(ucell, out_chg, global_out_dir, nspin, npol);
+    dftu_io::output(*dftu_ptr, ucell, out_chg, global_out_dir, nspin, npol);
 
     /// use the converged occupation matrix for next MD/Relax SCF calculation
     if (conv_esolver)

source/source_pw/module_pwdft/CMakeLists.txt

@@ -13,6 +13,7 @@ list(APPEND objects
     op_pw_exx_ace.cpp
     op_pw_exx_pot.cpp
     dftu_base.cpp
+    dftu_output.cpp
     dftu_pw.cpp
     setup_pot.cpp
     setup_pwrho.cpp

source/source_pw/module_pwdft/dftu_base.cpp

@@ -14,72 +14,9 @@
 #include <sstream>
 #include <vector>
 
-// local inline helpers for eigenvalue calculation (mirrors dftu_io.cpp)
-inline void JacobiRotate(std::vector<std::vector<double>>& A, int p, int q, int n)
-{
-    if (std::abs(A[p][q]) > 1e-10)
-    {
-        double r = (A[q][q] - A[p][p]) / (2.0 * A[p][q]);
-        double t = 0.0;
-        if (r >= 0)
-        {
-            t = 1.0 / (r + sqrt(1.0 + r * r));
-        }
-        else
-        {
-            t = -1.0 / (-r + sqrt(1.0 + r * r));
-        }
-        double c = 1.0 / sqrt(1.0 + t * t);
-        double s = t * c;
-
-        A[p][p] -= t * A[p][q];
-        A[q][q] += t * A[p][q];
-        A[p][q] = A[q][p] = 0.0;
-
-        for (int k = 0; k < n; k++)
-        {
-            if (k != p && k != q)
-            {
-                double Akp = c * A[k][p] - s * A[k][q];
-                double Akq = s * A[k][p] + c * A[k][q];
-                A[k][p] = A[p][k] = Akp;
-                A[k][q] = A[q][k] = Akq;
-            }
-        }
-    }
-}
-
-inline std::vector<double> CalculateEigenvalues(std::vector<std::vector<double>>& A, int n)
-{
-    std::vector<double> eigenvalues(n);
-    while (true)
-    {
-        int p = 0, q = 1;
-        for (int i = 0; i < n; i++)
-        {
-            for (int j = i + 1; j < n; j++)
-            {
-                if (std::abs(A[i][j]) > std::abs(A[p][q]))
-                {
-                    p = i;
-                    q = j;
-                }
-            }
-        }
-
-        if (std::abs(A[p][q]) < 1e-10)
-        {
-            for (int i = 0; i < n; i++)
-            {
-                eigenvalues[i] = A[i][i];
-            }
-            break;
-        }
-
-        JacobiRotate(A, p, q, n);
-    }
-    return eigenvalues;
-}
+// local inline helpers for eigenvalue calculation (JacobiRotate, CalculateEigenvalues)
+// have been migrated to dftu_output.cpp, where they are used by dftu_io::write_occup_m.
+// mohan refactored 2025-11-08
 
 // static member definitions (mohan add 2025-11-06)
 double Plus_U_Base::energy_u = 0.0;
@@ -584,237 +521,6 @@ void Plus_U_Base::set_locale_flat(const int iat, const int l, const int spin,
 }
 
 
-void Plus_U_Base::output(const UnitCell& ucell,
-                         bool out_chg,
-                         const std::string& global_out_dir,
-                         int nspin,
-                         int npol)
-{
-    ModuleBase::TITLE("Plus_U_Base", "output");
-
-    GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
-    GlobalV::ofs_running << " | #DFT+U INFORMATION# |" << std::endl;
-    GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
-
-    for (int T = 0; T < ucell.ntype; T++)
-    {
-        const int NL = ucell.atoms[T].nwl + 1;
-
-        for (int L = 0; L < NL; L++)
-        {
-            const int N = ucell.atoms[T].l_nchi[L];
-
-            if (L >= get_orbital_corr(T) && has_correlated_orbital(T))
-            {
-                if (L != get_orbital_corr(T))
-                {
-                    continue;
-                }
-
-                if (!Yukawa)
-                {
-                    GlobalV::ofs_running << " Type=" << T+1 << " L=" << L << " ORBITAL=" << 0
-                                         << " U=" << this->U[T] * ModuleBase::Ry_to_eV << " eV" << std::endl;
-                }
-                else
-                {
-                    for (int n = 0; n < N; n++)
-                    {
-                        if (n != 0)
-                        {
-                            continue;
-                        }
-                        double Ueff = (this->U_Yukawa[T][L][n] - this->J_Yukawa[T][L][n]) * ModuleBase::Ry_to_eV;
-                        GlobalV::ofs_running << " Type=" << T+1 << " L=" << L << "  ORBITAL=" << n
-                                             << " U=" << this->U_Yukawa[T][L][n] * ModuleBase::Ry_to_eV << " eV"
-                                             << " J=" << this->J_Yukawa[T][L][n] * ModuleBase::Ry_to_eV << " eV"
-                                             << std::endl;
-                    }
-                }
-            }
-        }
-    }
-
-    GlobalV::ofs_running << " Local Occupation Matrices for each atom" << std::endl;
-    this->write_occup_m(ucell, GlobalV::ofs_running, true, nspin, npol);
-
-    // Write dm_onsite.txt
-    if (out_chg && GlobalV::MY_RANK == 0)
-    {
-        std::ofstream ofdftu;
-        ofdftu.open(global_out_dir + "dm_onsite.txt");
-        if (!ofdftu)
-        {
-            ModuleBase::WARNING_QUIT("Plus_U_Base::output", "Can't create file dm_onsite.txt");
-        }
-        this->write_occup_m(ucell, ofdftu, false, nspin, npol);
-        ofdftu.close();
-    }
-
-    GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
-    GlobalV::ofs_running << " | # END DFT+U INFO    |" << std::endl;
-    GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl << std::endl;
-
-    return;
-}
-
-
-void Plus_U_Base::write_occup_m(const UnitCell& ucell,
-                                std::ofstream& ofs,
-                                bool diag,
-                                int nspin,
-                                int npol)
-{
-    ModuleBase::TITLE("Plus_U_Base", "write_occup_m");
-
-    if (GlobalV::MY_RANK != 0)
-    {
-        return;
-    }
-
-    for (int T = 0; T < ucell.ntype; T++)
-    {
-        if (!has_correlated_orbital(T))
-        {
-            continue;
-        }
-        const int NL = ucell.atoms[T].nwl + 1;
-        const int LC = get_orbital_corr(T);
-
-        for (int I = 0; I < ucell.atoms[T].na; I++)
-        {
-            const int iat = ucell.itia2iat(T, I);
-
-            for (int l = 0; l < NL; l++)
-            {
-                if (l != get_orbital_corr(T))
-                {
-                    continue;
-                }
-
-                const int N = ucell.atoms[T].l_nchi[l];
-
-                for (int n = 0; n < N; n++)
-                {
-                    if (n != 0)
-                    {
-                        continue;
-                    }
-
-                    ofs << "\n Atom= " << iat+1;
-                    ofs << " L= " << l;
-                    ofs << " ORBITAL= " << n << std::endl;
-
-                    if (nspin == 1 || nspin == 2)
-                    {
-                        double sum0[2];
-                        for (int is = 0; is < 2; is++)
-                        {
-                            if (diag)
-                            {
-                                std::vector<std::vector<double>> A(2 * l + 1, std::vector<double>(2 * l + 1));
-                                for (int m0 = 0; m0 < 2 * l + 1; m0++)
-                                {
-                                    for (int m1 = 0; m1 < 2 * l + 1; m1++)
-                                    {
-                                        A[m0][m1] = locale[iat][l][n][is](m0, m1);
-                                    }
-                                }
-                                std::vector<double> eigenvalues = CalculateEigenvalues(A, 2 * l + 1);
-                                sum0[is] = 0.0;
-                                ofs << " Eigenvalues for spin=" << is+1 << std::endl;
-                                ofs << std::setprecision(8) << std::fixed;
-                                for (int i = 0; i < 2 * l + 1; i++)
-                                {
-                                    ofs << std::setw(12) << eigenvalues[i];
-                                    sum0[is] += eigenvalues[i];
-                                }
-                                ofs << std::endl;
-                                ofs << " sum is " << std::setw(12) << sum0[is] << std::endl;
-                            }
-                            ofs << " spin= " << is+1 << std::endl;
-                            ofs << std::setprecision(8) << std::fixed;
-                            for (int m0 = 0; m0 < 2 * l + 1; m0++)
-                            {
-                                for (int m1 = 0; m1 < 2 * l + 1; m1++)
-                                {
-                                    ofs << std::setw(12)
-                                        << locale[iat][l][n][is](m0, m1);
-                                }
-                                ofs << std::endl;
-                            }
-                        }
-                        if (diag)
-                        {
-                            ofs << std::setw(12) << std::setprecision(8)
-                                << std::fixed << " Magnetism for atom " << iat+1 << ": " << sum0[0] - sum0[1]
-                                << std::endl;
-                        }
-                    }
-                    else if (nspin == 4) // SOC
-                    {
-                        if (diag)
-                        {
-                            double sum0[4];
-                            std::vector<std::vector<double>> A(2 * l + 1, std::vector<double>(2 * l + 1));
-                            int index = 0;
-                            for (int is = 0; is < 4; is++)
-                            {
-                                for (int m0 = 0; m0 < 2 * l + 1; m0++)
-                                {
-                                    for (int m1 = 0; m1 < 2 * l + 1; m1++)
-                                    {
-                                        A[m0][m1] = locale[iat][l][n][0].c[index];
-                                        index++;
-                                    }
-                                }
-                                std::vector<double> eigenvalues = CalculateEigenvalues(A, 2 * l + 1);
-                                sum0[is] = 0.0;
-                                ofs << " Eigenvalues for is=" << is << std::endl;
-                                ofs << std::setprecision(8) << std::fixed;
-                                for (int i = 0; i < 2 * l + 1; i++)
-                                {
-                                    ofs << std::setw(12) << eigenvalues[i];
-                                    sum0[is] += eigenvalues[i];
-                                }
-                                ofs << std::endl;
-                                ofs << " sum is " << std::setw(12) << sum0[is] << std::endl;
-                            }
-                            ofs << std::setw(12) << std::setprecision(8)
-                                << std::fixed << " Magnetism for atom " << iat + 1 << ": "
-                                << sum0[1] << " " << sum0[2] << " " << sum0[3] << std::endl;
-                        }
-                        else
-                        {
-                            for (int m0 = 0; m0 < 2 * l + 1; m0++)
-                            {
-                                for (int ipol0 = 0; ipol0 < npol; ipol0++)
-                                {
-                                    const int m0_all = m0 + (2 * l + 1) * ipol0;
-
-                                    for (int m1 = 0; m1 < 2 * l + 1; m1++)
-                                    {
-                                        for (int ipol1 = 0; ipol1 < npol; ipol1++)
-                                        {
-                                            int m1_all = m1 + (2 * l + 1) * ipol1;
-                                            ofs << std::setw(12) << std::setprecision(8) << std::fixed
-                                                << locale[iat][l][n][0](m0_all, m1_all);
-                                        }
-                                    }
-                                    ofs << std::endl;
-                                }
-                            }
-                        }
-                    }
-                } // n
-            }     // l
-        }         // I
-    }             // T
-
-    return;
-}
-
-
 void Plus_U_Base::read_occup_m(const UnitCell& ucell,
                                const std::string& fn,
                                const std::string& init_chg,

source/source_pw/module_pwdft/dftu_base.h

@@ -55,6 +55,10 @@ class Plus_U_Base
     static bool has_correlated_orbital(int it) { return orbital_corr[it] != -1; }
     static const int* get_orbital_corr_data() { return orbital_corr.data(); }
 
+    // mohan add 2025-11-08 for dftu_io::output free function
+    double get_U_Yukawa(int it, int l, int n) const { return U_Yukawa[it][l][n]; }
+    double get_J_Yukawa(int it, int l, int n) const { return J_Yukawa[it][l][n]; }
+
     static double get_energy() { return energy_u; }
     static void set_energy(const double &e) { energy_u = e; }
     static void set_double_energy() { energy_u *= 2.0; }
@@ -187,22 +191,12 @@ class Plus_U_Base
     std::vector<std::vector<std::vector<std::vector<ModuleBase::matrix>>>> locale_save;
 
     //=============================================================
-    // In dftu_io.cpp
-    // For reading/writing/broadcasting/copying relevant data structures
+    // output() and write_occup_m() have been extracted to free functions
+    // in source_pw/module_pwdft/dftu_output.cpp as dftu_io::output and
+    // dftu_io::write_occup_m. They access Plus_U_Base via public getters.
+    // mohan refactored 2025-11-08
     //=============================================================
-  public:
-    void output(const UnitCell& ucell,
-                bool out_chg,
-                const std::string& global_out_dir,
-                int nspin,
-                int npol);
-
   protected:
-    void write_occup_m(const UnitCell& ucell,
-                       std::ofstream& ofs,
-                       bool diag,
-                       int nspin,
-                       int npol);
     void read_occup_m(const UnitCell& ucell,
                       const std::string& fn,
                       const std::string& init_chg,