DFT+U Refactor

docs/advanced/elec_properties/band.md

@@ -8,7 +8,7 @@ out_chg 1
 ```
 
 With this input parameter, the converged charge density will be output in the files such as `chgs1.cube`, `chgs2.cube`, etc.
-Then, one can use the same `STRU` file, pseudopotential files and atomic orbital files (and the local density matrix file onsite.dm if DFT+U is used) to do a non-self-consistent (NSCF) calculation. In this example, the potential is constructed from the ground-state charge density from the proceeding calculation. Now the INPUT file is like:
+Then, one can use the same `STRU` file, pseudopotential files and atomic orbital files (and the local density matrix file dm_onsite.txt if DFT+U is used) to do a non-self-consistent (NSCF) calculation. In this example, the potential is constructed from the ground-state charge density from the proceeding calculation. Now the INPUT file is like:
 
 ```
 INPUT_PARAMETERS

docs/advanced/elec_properties/dos.md

@@ -10,7 +10,7 @@ out_chg              1
 ```
 
 this will produce the converged charge density, which is contained in the file SPIN1_CHG.cube.
-Then, use the same `STRU` file, pseudopotential file and atomic orbital file (and the local density matrix file onsite.dm if DFT+U is used) to do a non-self-consistent calculation. In this example, the potential is constructed from the ground-state charge density from the proceeding calculation. Now the INPUT file is like:
+Then, use the same `STRU` file, pseudopotential file and atomic orbital file (and the local density matrix file dm_onsite.txt if DFT+U is used) to do a non-self-consistent calculation. In this example, the potential is constructed from the ground-state charge density from the proceeding calculation. Now the INPUT file is like:
 
 ```
 INPUT_PARAMETERS

docs/advanced/input_files/input-main.md

@@ -3736,10 +3736,10 @@
 - **Type**: Integer
 - **Description**: The parameter controls the form of occupation matrix control used.
   - 0: No occupation matrix control is performed, and the onsite density matrix will be calculated from wavefunctions in each SCF step.
-  - 1: The first SCF step will use an initial density matrix read from a file named initial_onsite.dm, but for later steps, the onsite density matrix will be updated.
-  - 2: The same onsite density matrix from initial_onsite.dm will be used throughout the entire calculation.
+  - 1: The first SCF step will use an initial density matrix read from a file named dm_onsite_ini.txt, but for later steps, the onsite density matrix will be updated.
+  - 2: The same onsite density matrix from dm_onsite_ini.txt will be used throughout the entire calculation.
 
-  > Note: The easiest way to create initial_onsite.dm is to run a DFT+U calculation, look for a file named onsite.dm in the OUT.prefix directory, and make replacements there. The format of the file is rather straight-forward.
+  > Note: The easiest way to create dm_onsite_ini.txt is to run a DFT+U calculation with out_chg=1, look for a file named dm_onsite.txt in the OUT.prefix directory, copy and rename it to dm_onsite_ini.txt. The file dm_onsite_ini.txt should be placed in the directory specified by read_file_dir. The format of the file is rather straight-forward.
 - **Default**: 0
 
 ### onsite_radius

docs/parameters.yaml

@@ -4387,10 +4387,10 @@ parameters:
     description: |
       The parameter controls the form of occupation matrix control used.
       * 0: No occupation matrix control is performed, and the onsite density matrix will be calculated from wavefunctions in each SCF step.
-      * 1: The first SCF step will use an initial density matrix read from a file named initial_onsite.dm, but for later steps, the onsite density matrix will be updated.
-      * 2: The same onsite density matrix from initial_onsite.dm will be used throughout the entire calculation.
+      * 1: The first SCF step will use an initial density matrix read from a file named dm_onsite_ini.txt, but for later steps, the onsite density matrix will be updated.
+      * 2: The same onsite density matrix from dm_onsite_ini.txt will be used throughout the entire calculation.
 
-      [NOTE] The easiest way to create initial_onsite.dm is to run a DFT+U calculation, look for a file named onsite.dm in the OUT.prefix directory, and make replacements there. The format of the file is rather straight-forward.
+      [NOTE] The easiest way to create dm_onsite_ini.txt is to run a DFT+U calculation with out_chg=1, look for a file named dm_onsite.txt in the OUT.prefix directory, copy and rename it to dm_onsite_ini.txt. The file dm_onsite_ini.txt should be placed in the directory specified by read_file_dir. The format of the file is rather straight-forward.
     default_value: "0"
     unit: ""
     availability: ""

source/source_esolver/esolver_ks_lcao.cpp

@@ -370,7 +370,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     }
 #endif
 
-    init_dftu_lcao<TK>(istep, iter, PARAM.inp, &(this->dftu), this->dmat.dm, ucell, this->chr.rho, this->pw_rho->nrxx);
+    init_dftu_lcao<TK>(istep, iter, PARAM.inp.dft_plus_u, &(this->dftu), this->dmat.dm, ucell, this->chr.rho, this->pw_rho->nrxx);
 
 #ifdef __MLALGO
     // the density matrixes of DeePKS have been updated in each iter
@@ -481,7 +481,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
 	const std::vector<std::vector<TK>>& dm_vec = this->dmat.dm->get_DMK_vector();
 
     // 1) calculate the local occupation number matrix and energy correction in DFT+U
-    finish_dftu_lcao<TK>(iter, conv_esolver, PARAM.inp, &(this->dftu), ucell, dm_vec, this->kv, this->p_chgmix->get_mixing_beta(), hamilt_lcao);
+    finish_dftu_lcao<TK>(iter, conv_esolver, PARAM.inp.dft_plus_u, PARAM.inp.out_chg[0], &(this->dftu), ucell, dm_vec, this->kv, this->p_chgmix->get_mixing_beta(), hamilt_lcao, PARAM.globalv.global_out_dir, PARAM.inp.nspin, PARAM.globalv.npol);
 
     // 2) for deepks, calculate delta_e, output labels during electronic steps
     this->deepks.delta_e(ucell, this->kv, this->orb_, this->pv, this->gd, dm_vec, this->pelec->f_en, PARAM.inp);

source/source_io/module_hs/write_vxc.hpp

@@ -207,7 +207,7 @@ void write_Vxc(const int nspin,
         &vxcs_R_ao[0],ucell,/*for paraV*/ kv, Hexxd, Hexxc, hamilt::Add_Hexx_Type::k);
     std::vector<std::vector<double>> e_orb_exx; // orbital energy (EXX)
 #endif
-    hamilt::OperatorDFTU<hamilt::OperatorLCAO<TK, TR>> vdftu_op_ao(&vxc_k_ao, kv.kvec_d, nullptr, nullptr, kv.isk);
+    hamilt::OperatorDFTU<hamilt::OperatorLCAO<TK, TR>> vdftu_op_ao(&vxc_k_ao, kv.kvec_d, nullptr, nullptr, kv.isk, PARAM.globalv.npol);
 
     // 4. calculate and write the MO-matrix Exc
     Parallel_2D p2d;

source/source_io/module_parameter/read_input_item_exx_dftu.cpp

@@ -834,10 +834,10 @@ void ReadInput::item_dftu()
         item.type = "Integer";
         item.description = R"(The parameter controls the form of occupation matrix control used.
 * 0: No occupation matrix control is performed, and the onsite density matrix will be calculated from wavefunctions in each SCF step.
-* 1: The first SCF step will use an initial density matrix read from a file named initial_onsite.dm, but for later steps, the onsite density matrix will be updated.
-* 2: The same onsite density matrix from initial_onsite.dm will be used throughout the entire calculation.
+* 1: The first SCF step will use an initial density matrix read from a file named dm_onsite_ini.txt, but for later steps, the onsite density matrix will be updated.
+* 2: The same onsite density matrix from dm_onsite_ini.txt will be used throughout the entire calculation.
 
-[NOTE] The easiest way to create initial_onsite.dm is to run a DFT+U calculation, look for a file named onsite.dm in the OUT.prefix directory, and make replacements there. The format of the file is rather straight-forward.)";
+[NOTE] The easiest way to create dm_onsite_ini.txt is to run a DFT+U calculation with out_chg=1, look for a file named dm_onsite.txt in the OUT.prefix directory, copy and rename it to dm_onsite_ini.txt. The file dm_onsite_ini.txt should be placed in the directory specified by read_file_dir. The format of the file is rather straight-forward.)";
         item.default_value = "0";
         item.unit = "";
         item.availability = "";

source/source_lcao/FORCE_STRESS.cpp

@@ -420,7 +420,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
             std::vector<std::vector<double>>* dmk_d = nullptr;
             std::vector<std::vector<std::complex<double>>>* dmk_c = nullptr;
             assign_dmk_ptr<T>(dmat.dm, dmk_d, dmk_c, PARAM.globalv.gamma_only_local);
-            dftu.force_stress(ucell, gd, dmk_d, dmk_c, pv, fsr_dftu, force_u, stress_u, kv);
+            dftu.force_stress(ucell, gd, dmk_d, dmk_c, pv, fsr_dftu, force_u, stress_u, kv, PARAM.globalv.npol);
         }
         else
         {

source/source_lcao/LCAO_set.cpp

@@ -78,10 +78,22 @@ void LCAO_domain::set_pot(
           &(pelec->f_en.etxc), &(pelec->f_en.vtxc));
     }
 
-    //! 3) initialize DFT+U
     if (inp.dft_plus_u)
     {
-        dftu.init(ucell, &pv, kv.get_nks(), &orb);
+        dftu.init(ucell, &pv,
+                  PARAM.globalv.npol,
+                  inp.nspin, inp.orbital_corr, inp.yukawa_lambda,
+                  PARAM.globalv.global_readin_dir,
+                  PARAM.globalv.global_out_dir,
+                  inp.init_chg,
+                  pv.nloc,
+                  PARAM.globalv.gamma_only_local,
+                  inp.ks_solver,
+                  inp.cal_force,
+                  inp.cal_stress,
+                  inp.device,
+                  inp.kpar,
+                  &orb);
     }
 
     //! 4) init exact exchange calculations

source/source_lcao/dftu_lcao.cpp

@@ -9,14 +9,14 @@ namespace ModuleESolver
 template <typename TK>
 void init_dftu_lcao(const int istep,
                      const int iter,
-                     const Input_para& inp,
+                     int dft_plus_u,
                      void* dftu,
                      void* dm,
                      const UnitCell& ucell,
                      double** rho,
                      const int nrxx)
 {
-    if (!inp.dft_plus_u)
+    if (!dft_plus_u)
     {
         return;
     }
@@ -36,15 +36,19 @@ void init_dftu_lcao(const int istep,
 template <typename TK>
 void finish_dftu_lcao(const int iter,
                        const bool conv_esolver,
-                       const Input_para& inp,
+                       int dft_plus_u,
+                       bool out_chg,
                        void* dftu,
                        const UnitCell& ucell,
                        const std::vector<std::vector<TK>>& dm_vec,
                        const K_Vectors& kv,
                        const double mixing_beta,
-                       void* hamilt_lcao)
+                       void* hamilt_lcao,
+                       const std::string& global_out_dir,
+                       int nspin,
+                       int npol)
 {
-    if (!inp.dft_plus_u)
+    if (!dft_plus_u)
     {
         return;
     }
@@ -54,7 +58,7 @@ void finish_dftu_lcao(const int iter,
 
     /// old DFT+U method calculates energy correction in esolver,
     /// new DFT+U method calculates energy in Hamiltonian
-    if (inp.dft_plus_u == 2)
+    if (dft_plus_u == 2)
     {
         if (dftu_ptr->omc != 2)
         {
@@ -63,7 +67,7 @@ void finish_dftu_lcao(const int iter,
         }
         dftu_ptr->cal_energy_correction(ucell, iter);
     }
-    dftu_ptr->output(ucell);
+    dftu_ptr->output(ucell, out_chg, global_out_dir, nspin, npol);
 
     /// use the converged occupation matrix for next MD/Relax SCF calculation
     if (conv_esolver)
@@ -75,15 +79,16 @@ void finish_dftu_lcao(const int iter,
 /// Template instantiation
 template void init_dftu_lcao<double>(const int istep,
                                       const int iter,
-                                      const Input_para& inp,
+                                      int dft_plus_u,
                                       void* dftu,
                                       void* dm,
                                       const UnitCell& ucell,
                                       double** rho,
                                       const int nrxx);
+
 template void init_dftu_lcao<std::complex<double>>(const int istep,
                                                     const int iter,
-                                                    const Input_para& inp,
+                                                    int dft_plus_u,
                                                     void* dftu,
                                                     void* dm,
                                                     const UnitCell& ucell,
@@ -92,21 +97,30 @@ template void init_dftu_lcao<std::complex<double>>(const int istep,
 
 template void finish_dftu_lcao<double>(const int iter,
                                         const bool conv_esolver,
-                                        const Input_para& inp,
+                                        int dft_plus_u,
+                                        bool out_chg,
                                         void* dftu,
                                         const UnitCell& ucell,
                                         const std::vector<std::vector<double>>& dm_vec,
                                         const K_Vectors& kv,
                                         const double mixing_beta,
-                                        void* hamilt_lcao);
+                                        void* hamilt_lcao,
+                                        const std::string& global_out_dir,
+                                        int nspin,
+                                        int npol);
+
 template void finish_dftu_lcao<std::complex<double>>(const int iter,
                                                       const bool conv_esolver,
-                                                      const Input_para& inp,
+                                                      int dft_plus_u,
+                                                      bool out_chg,
                                                       void* dftu,
                                                       const UnitCell& ucell,
                                                       const std::vector<std::vector<std::complex<double>>>& dm_vec,
                                                       const K_Vectors& kv,
                                                       const double mixing_beta,
-                                                      void* hamilt_lcao);
+                                                      void* hamilt_lcao,
+                                                      const std::string& global_out_dir,
+                                                      int nspin,
+                                                      int npol);
 
 } // namespace ModuleESolver

source/source_lcao/dftu_lcao.h

@@ -1,9 +1,10 @@
 #ifndef DFTU_LCAO_H
 #define DFTU_LCAO_H
 
+#include <string>
+#include <vector>
 #include "source_cell/unitcell.h"
 #include "source_cell/klist.h"
-#include "source_io/module_parameter/input_parameter.h"
 
 namespace ModuleESolver
 {
@@ -16,7 +17,7 @@ namespace ModuleESolver
  *
  * @param istep Current ionic step
  * @param iter Current SCF iteration
- * @param inp Input parameters
+ * @param dft_plus_u DFT+U mode (0=disabled, 1=old, 2=new)
  * @param dftu DFT+U object
  * @param dm Density matrix
  * @param ucell Unit cell
@@ -26,7 +27,7 @@ namespace ModuleESolver
 template <typename TK>
 void init_dftu_lcao(const int istep,
                      const int iter,
-                     const Input_para& inp,
+                     int dft_plus_u,
                      void* dftu,
                      void* dm,
                      const UnitCell& ucell,
@@ -41,24 +42,32 @@ void init_dftu_lcao(const int istep,
  *
  * @param iter Current SCF iteration
  * @param conv_esolver Whether ESolver has converged
- * @param inp Input parameters
+ * @param dft_plus_u DFT+U mode (0=disabled, 1=old, 2=new)
+ * @param out_chg Whether to output dm_onsite.txt
  * @param dftu DFT+U object
  * @param ucell Unit cell
  * @param dm_vec Density matrix vector
  * @param kv K-vectors
  * @param mixing_beta Mixing beta parameter
  * @param hamilt_lcao Hamiltonian LCAO object
+ * @param global_out_dir Output directory for dm_onsite.txt
+ * @param nspin Number of spin channels (1, 2, or 4)
+ * @param npol Number of polarizations
  */
 template <typename TK>
 void finish_dftu_lcao(const int iter,
                        const bool conv_esolver,
-                       const Input_para& inp,
+                       int dft_plus_u,
+                       bool out_chg,
                        void* dftu,
                        const UnitCell& ucell,
                        const std::vector<std::vector<TK>>& dm_vec,
                        const K_Vectors& kv,
                        const double mixing_beta,
-                       void* hamilt_lcao);
+                       void* hamilt_lcao,
+                       const std::string& global_out_dir,
+                       int nspin,
+                       int npol);
 
 } // namespace ModuleESolver
 

source/source_lcao/hamilt_lcao.cpp

@@ -225,9 +225,10 @@ HamiltLCAO<TK, TR>::HamiltLCAO(const UnitCell& ucell,
             {
                 plus_u = new OperatorDFTU<OperatorLCAO<TK, TR>>(this->hsk,
                                                               this->kv->kvec_d,
-															  this->hR, // no explicit call yet
-															  p_dftu, // mohan add 2025-11-07
-															  this->kv->isk);
+															  this->hR,
+															  p_dftu,
+															  this->kv->isk,
+															  PARAM.globalv.npol);
             }
             else
             {
@@ -381,9 +382,10 @@ HamiltLCAO<TK, TR>::HamiltLCAO(const UnitCell& ucell,
             {
                 plus_u = new OperatorDFTU<OperatorLCAO<TK, TR>>(this->hsk,
                                                               this->kv->kvec_d,
-															  this->hR, // no explicit call yet
-															  p_dftu, // mohan add 2025-11-07
-                                                              this->kv->isk);
+															  this->hR,
+															  p_dftu,
+                                                              this->kv->isk,
+                                                              PARAM.globalv.npol);
             }
             else
             {