Update out_chg 2 and out_pot 3, and elf

docs/advanced/output_files/output-specification.md

@@ -49,11 +49,18 @@ All output file naming conventions can be found in the online documentation (wit
 | `s12` | Non-collinear spin calculation |
 | `k#` | k-point index (e.g., `k1`, `k2`) |
 | `g#` | Ionic step index for relax/md (e.g., `g1`, `g2`) |
+| `ini` | Initial state (before electronic iteration), used before or after `g#` |
 
 **Important:**
 - All index numbers start from 1 (not 0)
 - For Gamma-only algorithm in LCAO, no `k` index is included
 - Overlap matrix `s` does not distinguish spin, so only one matrix is output
+- For initial charge density output (`out_chg = 2`):
+  - `out_freq_ion = 0`: `chg_ini.cube` (nspin=1) or `chgs{#}_ini.cube` (nspin=2/4) - output at every step (overwrite same file)
+  - `out_freq_ion > 0`: `chgg{#}_ini.cube` (nspin=1) or `chgs{#}g{#}_ini.cube` (nspin=2/4) - output every out_freq_ion steps
+- For initial potential output (`out_pot = 3`):
+  - `out_freq_ion = 0`: `pot_ini.cube` (nspin=1) or `pots{#}_ini.cube` (nspin=2/4) - output at every step (overwrite same file)
+  - `out_freq_ion > 0`: `potg{#}_ini.cube` (nspin=1) or `pots{#}g{#}_ini.cube` (nspin=2/4) - output every out_freq_ion steps
 
 ### 2.2 Examples
 
@@ -62,7 +69,22 @@ All output file naming conventions can be found in the online documentation (wit
 | `chgs1.cube` | Charge density, spin 1 |
 | `chgs2.cube` | Charge density, spin 2 |
 | `chgs3.cube` | Charge density, spin 3 (non-collinear with SOC) |
+| `chg_ini.cube` | Initial charge density (out_freq_ion=0, nspin=1) |
+| `chgs1_ini.cube` | Initial charge density (out_freq_ion=0, spin 1, nspin=2/4) |
+| `chgg1_ini.cube` | Initial charge density, geometry step 1 (nspin=1) |
+| `chgs1g1_ini.cube` | Initial charge density, spin 1, geometry step 1 (nspin=2/4) |
+| `pot_ini.cube` | Initial potential (out_freq_ion=0, nspin=1) |
+| `pots1_ini.cube` | Initial potential (out_freq_ion=0, spin 1, nspin=2/4) |
+| `potg1_ini.cube` | Initial potential, geometry step 1 (nspin=1) |
+| `pots1g1_ini.cube` | Initial potential, spin 1, geometry step 1 (nspin=2/4) |
 | `pots1.cube` | Local potential, spin 1 |
+| `elftot.cube` | ELF, total (nspin=1/4) |
+| `elfs1.cube` | ELF, spin 1 (nspin=2) |
+| `elfs2.cube` | ELF, spin 2 (nspin=2) |
+| `elftot.cube` | ELF, total (nspin=2) |
+| `elftotg1.cube` | ELF, total, geometry step 1 (nspin=1/4) |
+| `elfs1g1.cube` | ELF, spin 1, geometry step 1 (nspin=2) |
+| `elftotg1.cube` | ELF, total, geometry step 1 (nspin=2) |
 | `eig_occ.txt` | Eigenvalues and occupations |
 | `doss1g1_nao.txt` | DOS, spin 1, geometry step 1, NAO basis |
 | `wf_pw.dat` | Wavefunction, plane wave basis |
@@ -79,8 +101,25 @@ All output file naming conventions can be found in the online documentation (wit
 | `band.txt` | Band structure |
 | `chgs1.cube`, `chgs2.cube` | Charge density (spin 1, spin 2) |
 | `chg.cube` | Total charge density |
+| Initial charge density (out_chg=2) | |
+| `chg_ini.cube` | Initial charge density (out_freq_ion=0, nspin=1) |
+| `chgs{#}_ini.cube` | Initial charge density (out_freq_ion=0, spin {#}, nspin=2/4) |
+| `chgg{#}_ini.cube` | Initial charge density (out_freq_ion>0, geometry step {#}, nspin=1) |
+| `chgs{#}g{#}_ini.cube` | Initial charge density (out_freq_ion>0, spin {#}, geometry step {#}, nspin=2/4) |
+| Initial potential (out_pot=3) | |
+| `pot_ini.cube` | Initial potential (out_freq_ion=0, nspin=1) |
+| `pots{#}_ini.cube` | Initial potential (out_freq_ion=0, spin {#}, nspin=2/4) |
+| `potg{#}_ini.cube` | Initial potential (out_freq_ion>0, geometry step {#}, nspin=1) |
+| `pots{#}g{#}_ini.cube` | Initial potential (out_freq_ion>0, spin {#}, geometry step {#}, nspin=2/4) |
 | `taus1.cube`, `taus2.cube` | Kinetic energy density (tau) |
 | `pots1.cube`, `pots2.cube` | Local potential |
+| ELF (out_elf) | |
+| `elftot.cube` | ELF, total (nspin=1/4) |
+| `elfs1.cube`, `elfs2.cube` | ELF, spin 1/2 (nspin=2) |
+| `elftot.cube` | ELF, total (nspin=2) |
+| `elftotg{#}.cube` | ELF, total, geometry step {#} (nspin=1/4) |
+| `elfs{#}g{#}.cube` | ELF, spin {#}, geometry step {#} (nspin=2) |
+| `elftotg{#}.cube` | ELF, total, geometry step {#} (nspin=2) |
 
 ## 3. File Format Standards
 

examples/18_md/01_lcao_gamma_Si8/INPUT_1 → examples/18_md/02_LCAO_NVT_Si8/INPUT

@@ -3,6 +3,7 @@ INPUT_PARAMETERS
 suffix                 Si_nhc_nvt
 calculation            md
 nbands                 20
+nspin                  2
 symmetry               0
 pseudo_dir             ../../../tests/PP_ORB
 orbital_dir            ../../../tests/PP_ORB
@@ -15,7 +16,7 @@ scf_nmax               100
 #Parameters (3.Basis)
 basis_type             lcao
 ks_solver              genelpa
-gamma_only             1
+gamma_only             0 
 
 #Parameters (4.Smearing)
 smearing_method        gaussian
@@ -28,8 +29,15 @@ chg_extrap             second-order
 
 #Parameters (6.MD)
 md_type                nvt
-md_nstep               10
+md_nstep               5
 md_dt                  1
 md_tfirst              300
 md_tfreq               1
-md_tchain              1
+md_tchain              1
+
+# output initial charge density
+out_chg                2  4
+out_pot                3  4
+out_elf                1  4
+out_dipole             1
+out_freq_ion           2

examples/18_md/02_LCAO_NVT_Si8/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+2 2 2 0 0 0

examples/18_md/02_LCAO_NVT_Si8/STRU

@@ -0,0 +1,28 @@
+ATOMIC_SPECIES
+Si  28.085  Si_ONCV_PBE-1.0.upf
+
+NUMERICAL_ORBITAL
+Si_gga_8au_60Ry_2s2p1d.orb
+
+LATTICE_CONSTANT
+1.8897270    # 1 Angstrom = 1.8897270 bohr
+
+LATTICE_VECTORS
+5.43090      0.00000      0.00000
+0.00000      5.43090      0.00000
+0.00000      0.00000      5.43090
+
+ATOMIC_POSITIONS
+Direct 
+
+Si
+0.0
+8
+0.000	0.000	0.000	1	1	1
+0.000	0.500	0.500	1	1	1
+0.500	0.000	0.500	1	1	1
+0.500	0.500	0.000	1	1	1
+0.250	0.250	0.250	1	1	1
+0.250	0.750	0.750	1	1	1
+0.750	0.250	0.750	1	1	1
+0.750	0.750	0.250	1	1	1

examples/22_rt-tddft/02_H2_velocity_gauge/INPUT

@@ -20,7 +20,7 @@ smearing_method      gauss
 
 #Parameters (5.MD Parameters)
 md_type              nve
-md_nstep             1000
+md_nstep             5
 md_dt                0.005
 md_tfirst            0
 
@@ -46,6 +46,9 @@ out_chg              1
 out_efield           1
 out_dipole           1
 
+out_elf              1
+out_freq_ion         0
+
 
 ### [1] Energy cutoff determines the quality of numerical quadratures in your calculations.
 ###     So it is strongly recommended to test whether your result (such as converged SCF energies) is

source/source_esolver/esolver_dm2rho.cpp

@@ -79,7 +79,9 @@ void ESolver_DM2rho<TK, TR>::runner(UnitCell& ucell, const int istep)
                                       this->pelec->eferm.get_efval(is),
                                       &(ucell),
                                       3,
-                                      1);
+                                      1,
+                                      PARAM.globalv.two_fermi,
+                                      false);
     }
 
     ModuleBase::timer::end("ESolver_DM2rho", "runner");

source/source_esolver/esolver_double_xc.cpp

@@ -158,7 +158,8 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
 	}
 
     XC_Functional::set_xc_type(PARAM.inp.deepks_out_base);
-    this->pelec_base->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
+    elecstate::init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, istep, 
+		    PARAM.globalv.global_out_dir, PARAM.inp, this->pelec_base);
     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); 
 
     // DMR should be same size with Hamiltonian(R)

source/source_esolver/esolver_fp.cpp

@@ -175,9 +175,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
     //! set direction of magnetism, used in non-collinear case 
     elecstate::cal_ux(ucell);
 
-    //! output the initial charge density and potential
-    ModuleIO::write_chg_init(ucell, this->Pgrid, this->chr, this->pelec->eferm, istep, PARAM.inp);
-//    ModuleIO::write_pot_init(ucell, this->Pgrid, this->pelec, istep, PARAM.inp); 
+    //! output the initial charge density
+    ModuleIO::write_chg_init(ucell, this->Pgrid, this->chr, this->pelec->eferm, istep,
+                             PARAM.globalv.global_out_dir, PARAM.inp, PARAM.globalv.two_fermi);
 
     return;
 }

source/source_esolver/esolver_ks_lcao.cpp

@@ -195,7 +195,8 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
         this->dmat.dm->cal_DMR();
     }
     // 13.2) init_scf, should be before_scf? mohan add 2025-03-10
-    this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
+    elecstate::init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric,
+                          istep, PARAM.globalv.global_out_dir, PARAM.inp, this->pelec);
 
 #ifdef __MLALGO
     // 14) initialize DM2(R) of DeePKS, the DM2(R) is different from DM(R)

source/source_esolver/esolver_ks_lcao_tddft.cpp

@@ -177,7 +177,8 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)
                                         GlobalV::ofs_running,
                                         GlobalV::ofs_warning);
             this->exx_nao.before_scf(ucell, this->kv, this->orb_, this->p_chgmix, totstep, PARAM.inp);
-            this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
+            elecstate::init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, istep, 
+			    PARAM.globalv.global_out_dir, PARAM.inp, this->pelec);
 
             if (totstep <= PARAM.inp.td_tend + 1)
             {

source/source_esolver/esolver_ks_pw.cpp

@@ -165,7 +165,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
               this->chr, this->locpp, this->ppcell, this->dftu, this->vsep_cell,
               this->stp.template get_psi_t<T, Device>(), 
 	      this->p_hamilt, 
-	      this->pw_wfc, this->pw_rhod, PARAM.inp);
+	      this->pw_wfc, this->pw_rhod, PARAM.globalv.global_out_dir, PARAM.inp);
 
     // setup psi (electronic wave functions)
     this->stp.init(this->p_hamilt);

source/source_esolver/esolver_of.cpp

@@ -94,7 +94,9 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp)
     this->init_elecstate(ucell);
 
     // calculate the total local pseudopotential in real space
-    this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm); // atomic_rho, v_of_rho, set_vrs
+    const int istep=0;
+    elecstate::init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep, 
+		    PARAM.globalv.global_out_dir, PARAM.inp, this->pelec);
 
     // liuyu move here 2023-10-09
     // D in uspp need vloc, thus behind init_scf()
@@ -188,8 +190,6 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
     //! 1) call before_scf() of ESolver_FP
     ESolver_FP::before_scf(ucell, istep);
 
-
-
     if (ucell.cell_parameter_updated)
     {
         this->dV_ = ucell.omega / this->pw_rho->nxyz;
@@ -232,7 +232,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
         }
     }
 
-    this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm);
+    elecstate::init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep, PARAM.globalv.global_out_dir, PARAM.inp, this->pelec);
 
     Symmetry_rho::symmetrize_rho(PARAM.inp.nspin, this->chr, this->pw_rho, ucell.symm);
 

source/source_esolver/lcao_others.cpp

@@ -172,7 +172,8 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
     elecstate::cal_ux(ucell);
 
     // pelec should be initialized before these calculations
-    this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
+    elecstate::init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, 
+		    istep, PARAM.globalv.global_out_dir, PARAM.inp, this->pelec);
 
     // self consistent calculations for electronic ground state
     if (cal_type == "get_pchg")

source/source_estate/CMakeLists.txt

@@ -1,5 +1,6 @@
 list(APPEND objects
     elecstate.cpp
+    init_scf.cpp
     elecstate_energy_terms.cpp
     elecstate_energy.cpp
     elecstate_exx.cpp

source/source_estate/elecstate.cpp

@@ -24,23 +24,6 @@ void ElecState::init_nelec_spin()
     }
 }
 
-void ElecState::init_scf(const UnitCell& ucell,
-                         const Parallel_Grid& pgrid,
-                         const ModuleBase::ComplexMatrix& strucfac, 
-                         const bool* numeric,
-                         ModuleSymmetry::Symmetry& symm, 
-                         const void* wfcpw)
-{
-    //! core correction potential.
-    this->charge->set_rho_core(ucell,strucfac, numeric);
-
-    //! renormalize the charge density
-    this->charge->renormalize_rho();
-
-    //! initialize the potential
-    this->pot->init_pot(this->charge);
-}
-
 
 void ElecState::init_ks(Charge* chr_in, // pointer for class Charge
                         const K_Vectors* klist_in,

source/source_estate/elecstate.h

@@ -94,20 +94,8 @@ class ElecState
         return;
     }
 
-    /**
-     * @brief Init rho_core, init rho, renormalize rho, init pot
-     *
-     * @param ucell unit cell
-     * @param strucfac structure factor
-     * @param symm symmetry
-     * @param wfcpw PW basis for wave function if needed
-     */
-    void init_scf(const UnitCell& ucell,
-                  const Parallel_Grid& pgrid,
-                  const ModuleBase::ComplexMatrix& strucfac,
-                  const bool* numeric,
-                  ModuleSymmetry::Symmetry& symm,
-                  const void* wfcpw = nullptr);
+
+
     std::string classname = "elecstate";
 
     int iter = 0;                                  ///< scf iteration
@@ -162,5 +150,26 @@ class ElecState
     bool skip_weights = false;
 };
 
+/**
+ * @brief Init rho_core, init rho, renormalize rho, init pot
+ *
+ * @param ucell unit cell
+ * @param pgrid parallel grid
+ * @param strucfac structure factor
+ * @param numeric numeric flag
+ * @param istep ionic step index
+ * @param out_dir output directory
+ * @param inp input parameters
+ * @param pelec pointer to ElecState
+ */
+void init_scf(const UnitCell& ucell,
+              const Parallel_Grid& pgrid,
+              const ModuleBase::ComplexMatrix& strucfac,
+              const bool* numeric,
+              const int istep,
+              const std::string& out_dir,
+              const Input_para& inp,
+              ElecState* pelec);
+
 } // namespace elecstate
 #endif