DFT+U Refactor (#7538)

* format dftu_io.cpp

* update output formats of DFT+U

* update dftu, remove PARAM

* fix read_occup_m function

* DFT+U I/O: Fix multiple issues with onsite.dm format, precision, and MPI

This commit addresses several issues in the DFT+U I/O code related to the
onsite.dm file handling:

1. Format compatibility fix: Updated read_occup_m() to parse the new output format
   with labels Atom=, L=, ORBITAL=, and spin= instead of the old tokens atoms, L,
   zeta, and spin. The writer was already using the new format but the reader
   was not updated, causing new onsite.dm files to be unreadable.

2. Off-by-one fix: The writer outputs 1-based indices (iat+1, is+1) for human
   readability, but read_occup_m() was using these values directly as array
   indices. Added iat -= 1 and spin -= 1 to convert back to 0-based indices.

3. Precision fix: Added std::setprecision(8) << std::fixed for the collinear
   (nspin=1,2) diag=false output path in write_occup_m(). The eigenvalues path
   and SOC path already had this, but the matrix values path was missing it.
   This ensures stable numeric precision for restart data.

4. out_chg logic fix: The out_chg parameter is supposed to control whether
   onsite.dm is written, but the implementation was still checking if(!ofdftu)
   even when out_chg == false. Moved the file-open check and write operations
   inside the out_chg && MY_RANK == 0 block.

5. MPI fix: Only rank 0 opens the onsite.dm file, but all ranks were executing
   the if(!ofdftu) check and write_occup_m() call. Non-root ranks would see an
   unopened stream and potentially fail. Now all file operations are rank-0-only.

6. Header fix: Added missing #include <string> and #include <vector> to
   dftu_lcao.h. The header was using std::string and std::vector in function
   signatures but relying on transitive includes.

7. Documentation fix: Updated stale comments in dftu_lcao.h:
   - Changed @PARAM inp to @PARAM dft_plus_u in init_dftu_lcao()
   - Added documentation for global_out_dir, nspin, and npol parameters in
     finish_dftu_lcao()

8. Error handling: Replaced all exit(0) calls with ModuleBase::WARNING_QUIT()
   for consistent error handling across the codebase.

Files modified:
- source/source_lcao/module_dftu/dftu_io.cpp
- source/source_lcao/dftu_lcao.h
- source/source_lcao/dftu_lcao.cpp

* DFT+U I/O: Rename onsite.dm to dm_onsite.txt and adjust file paths

This commit renames the DFT+U occupation matrix files and adjusts their read/write paths:

1. File name changes:
   - onsite.dm → dm_onsite.txt (output from DFT+U calculations)
   - initial_onsite.dm → dm_onsite_ini.txt (user-provided initial occupation matrix)

2. Path adjustments:
   - dm_onsite.txt: Both written and read from global_out_dir (OUT.prefix)
     - Previously read from global_readin_dir
     - This ensures the file is always in the output directory
   - dm_onsite_ini.txt: Read from global_readin_dir (set via read_file_dir parameter)
     - Previously read from global_out_dir
     - This allows users to place initial files in any directory

3. Code changes:
   - dftu_io.cpp: Updated output filename and error messages
   - dftu.cpp: Updated read paths for both files
   - dftu_lcao.h and module_operator_lcao/dftu_lcao.cpp: Updated comments

4. Documentation updates:
   - read_input_item_exx_dftu.cpp: Updated omc parameter description
   - parameters.yaml: Updated omc parameter description
   - input-main.md: Updated omc parameter description
   - band.md and dos.md: Updated file references
   - tests/17_DS_DFTU/README.md: Updated test documentation
   - tests/17_DS_DFTU/run_scf_nscf.sh: Updated copy logic for dm_onsite.txt

5. Backward compatibility:
   - The new format with labels Atom=, L=, ORBITAL=, spin= was already in the writer
   - This commit updates the reader to parse the new format
   - 1-based indices in output are converted to 0-based when reading

Files modified:
- source/source_lcao/module_dftu/dftu_io.cpp
- source/source_lcao/module_dftu/dftu.cpp
- source/source_lcao/dftu_lcao.h
- source/source_lcao/module_operator_lcao/dftu_lcao.cpp
- source/source_io/module_parameter/read_input_item_exx_dftu.cpp
- docs/parameters.yaml
- docs/advanced/input_files/input-main.md
- docs/advanced/elec_properties/band.md
- docs/advanced/elec_properties/dos.md
- tests/17_DS_DFTU/README.md
- tests/17_DS_DFTU/run_scf_nscf.sh

* remove some usage of PARAM

* remove C++17 codes, reduce PARAM usage in dftu

* update

* remove all PARAM in dftu modules

* fix a bug

* fix bug about Yukawa potential

* update

* try fixing the bug?

* update the directory

---------

Co-authored-by: abacus_fixer <mohanchen@pku.eud.cn>

Author: mohanchen

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

@@ -3737,10 +3737,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

@@ -4389,10 +4389,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: ""

interfaces/ASE_interface/abacuslite/core.py

@@ -529,7 +529,8 @@ def fixed_density(self,
                                              'symmetry': 0,
                                              'out_band': 1,
                                              'kspacing': 0.0,       # overwrite
-                                             'gamma_only': False})  # overwrite
+                                             'gamma_only': False,
+                                             'read_file_dir': 'OUT.ABACUS'})  # overwrite
 
         profile = self.profile if profile is None else profile
 

source/source_esolver/esolver_ks_lcao.cpp

@@ -372,7 +372,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
@@ -483,7 +483,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/input_conv.cpp

@@ -212,7 +212,6 @@ void Input_Conv::Convert()
 
     if (PARAM.inp.dft_plus_u)
     {
-        Plus_U::Yukawa = PARAM.inp.yukawa_potential;
         Plus_U::omc = PARAM.inp.omc;
         Plus_U::orbital_corr = PARAM.inp.orbital_corr;
         Plus_U::uramping = PARAM.globalv.uramping;

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
         {