Fix part of 6608: Make IO output self-describing: running-log messages, cube data_desc, and ofs_running parameterization

docs/advanced/elec_properties/hs_matrix.md

@@ -120,7 +120,7 @@ Each file starts with a header:
  --- Ionic Step 1 ---
  # print H matrix in real space H(R)
  1 # number of spin directions
- 1 # spin index
+ 1 # spin index, E_Fermi = 10.3798 eV
  100 # number of localized basis
  50 # number of Bravais lattice vector R
 
@@ -138,6 +138,8 @@ Each file starts with a header:
  0 3 7 ...
 ```
 
+The spin index line may carry an optional `, E_Fermi = <value> eV` suffix when the Fermi energy is available (only for the H matrix). The value is printed in eV with 6 significant digits.
+
 The CSR format stores a sparse m × n matrix M in row form using three arrays (values, column indices, row pointers). According to Wikipedia:
 
 - The arrays **values** and **column indices** are of length NNZ (number of nonzero entries), and contain the non-zero values and the column indices of those values respectively.

source/source_esolver/esolver_dm2rho.cpp

@@ -68,25 +68,39 @@ void ESolver_DM2rho<TK, TR>::runner(BaseCell& basecell, const int istep)
     // this->pelec->psiToRho(*this->psi);
     LCAO_domain::dm2rho(this->dmat.dm->get_dmr_vec(), this->inp_->nspin, &this->chr, this->inp_->nelec, ucell.omega, false);
 
-    int nspin0 = this->inp_->nspin == 2 ? 2 : 1;
+    // Write every density channel (rho[0..3] for nspin=4) so the file set
+    // matches ctrl_output_fp: the init_chg=file reader reads chgs1..chgs<nspin>
+    // once chgs1.cube exists and aborts if any of them is missing.
+    const int nspin = this->inp_->nspin;
 
-    for (int is = 0; is < nspin0; is++)
+    for (int is = 0; is < nspin; is++)
     {
-        std::string fn = PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube";
+        // same naming as ctrl_output_fp and the init_chg=file reader:
+        // chg.cube for nspin=1, chgs<n>.cube otherwise
+        std::string fn = PARAM.globalv.global_out_dir + "chg";
+        if (nspin != 1)
+        {
+            fn += "s" + std::to_string(is + 1);
+        }
+        fn += ".cube";
 
         // write electron density
+        const std::string desc_chg = ModuleIO::make_data_desc("charge density", "magnetization density m", is, nspin);
+
         ModuleIO::write_vdata_palgrid(this->Pgrid,
                                       this->chr.rho[is],
                                       is,
-                                      this->inp_->nspin,
+                                      nspin,
                                       istep,
                                       fn,
                                       this->pelec->eferm.get_efval(is),
                                       &(ucell),
                                       3,
                                       1,
                                       PARAM.globalv.two_fermi,
-                                      false);
+                                      false,
+                                      GlobalV::ofs_running,
+                                      desc_chg);
     }
 
     ModuleBase::timer::end("ESolver_DM2rho", "runner");

source/source_esolver/esolver_fp.cpp

@@ -161,10 +161,10 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
     ModuleBase::TITLE("ESolver_FP", "after_scf");
 
     //! Output convergence information
-    ModuleIO::output_convergence_after_scf(conv_esolver, this->pelec->f_en.etot);
+    ModuleIO::output_convergence_after_scf(conv_esolver, this->pelec->f_en.etot, GlobalV::ofs_running);
 
     //! Write Fermi energy
-    ModuleIO::output_efermi(conv_esolver, this->pelec->eferm.ef);
+    ModuleIO::output_efermi(conv_esolver, this->pelec->eferm, GlobalV::ofs_running);
 
     //! Update delta_rho for charge extrapolation
     const module_charge::AtomicRhoCfg atomic_rho_cfg_after{
@@ -177,7 +177,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
 
     //! print out charge density, potential, elf, etc.
 	ModuleIO::ctrl_output_fp(ucell, *this->inp_, this->pelec, this->pw_big, this->pw_rhod, 
-			this->chr, this->solvent, this->Pgrid, istep); 
+			this->chr, this->solvent, this->Pgrid, istep, GlobalV::ofs_running); 
 
 }
 
@@ -262,7 +262,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
 
     //! output the initial charge density
     ModuleIO::write_chg_init(ucell, this->Pgrid, this->chr, this->pelec->eferm, istep,
-                             PARAM.globalv.global_out_dir, *this->inp_, PARAM.globalv.two_fermi);
+                             PARAM.globalv.global_out_dir, *this->inp_, PARAM.globalv.two_fermi, GlobalV::ofs_running);
 
     return;
 }

source/source_esolver/esolver_gets.cpp

@@ -208,7 +208,8 @@ void ESolver_GetS::runner(BaseCell& basecell, const int istep)
                         PARAM.globalv.global_matrix_dir,
                         PARAM.inp.calculation,
                         PARAM.inp.out_app_flag,
-                        PARAM.inp.nspin);
+                        PARAM.inp.nspin,
+                        GlobalV::ofs_running);
 
     if (this->inp_->out_mat_r[0])
     {
@@ -225,7 +226,8 @@ void ESolver_GetS::runner(BaseCell& basecell, const int istep)
                         PARAM.inp.calculation,
                         PARAM.inp.out_app_flag,
                         nlocal,
-                        PARAM.globalv.npol);
+                        PARAM.globalv.npol,
+                        GlobalV::ofs_running);
     }
 
     if (this->inp_->out_mat_ds[0])
@@ -251,7 +253,8 @@ void ESolver_GetS::runner(BaseCell& basecell, const int istep)
                              mat_R_options,
                              PARAM.globalv.gamma_only_local,
                              PARAM.globalv.npol,
-                             PARAM.globalv.nlocal);
+                             PARAM.globalv.nlocal,
+                             GlobalV::ofs_running);
     }
 
     ModuleBase::timer::end("ESolver_GetS", "runner");

source/source_esolver/esolver_ks.cpp

@@ -412,7 +412,7 @@ void ESolver_KS::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
                              istep);
 
     // 4) write band information to band.txt
-    ModuleIO::write_bands(*this->inp_, this->pelec->ekb, this->kv);
+    ModuleIO::write_bands(*this->inp_, this->pelec->ekb, this->kv, GlobalV::ofs_running);
 
 }
 

source/source_esolver/esolver_ks_lcao.cpp

@@ -338,7 +338,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(BaseCell& basecell)
 	    this->gd, this->psi, this->chr, hamilt_lcao,
 	    this->two_center_bundle_,
 	    this->orb_, this->pw_rho, this->pw_rhod,
-	    this->sf, this->locpp.vloc, this->exx_nao, this->exx_info_, this->solvent);
+	    this->sf, this->locpp.vloc, this->exx_nao, this->exx_info_, this->solvent, GlobalV::ofs_running);
 
     ModuleBase::timer::end("ESolver_KS_LCAO", "after_all_runners");
 }
@@ -649,7 +649,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
             this->orb_, this->pw_wfc, this->pw_rho, this->pw_big, this->sf,
             this->pw_rhod, this->locpp.vloc, this->solvent,
             this->rdmft_solver, this->deepks, this->exx_nao, this->exx_info_,
-            conv_esolver, this->scf_nmax_flag, istep);
+            conv_esolver, this->scf_nmax_flag, istep, GlobalV::ofs_running);
 
     //! 3) Clean up RA, which is used to serach for adjacent atoms
     if (!this->inp_->cal_force && !this->inp_->cal_stress)

source/source_esolver/esolver_ks_pw.cpp

@@ -333,7 +333,7 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
     pw::check_deltaspin_oscillation(iter, this->drho, this->p_chgmix, *this->inp_);
 
     // the output quantities
-    ModuleIO::ctrl_iter_pw(istep, iter, conv_esolver, this->stp, this->kv, this->pw_wfc, *this->inp_);
+    ModuleIO::ctrl_iter_pw(istep, iter, conv_esolver, this->stp, this->kv, this->pw_wfc, *this->inp_, GlobalV::ofs_running);
 }
 
 template <typename T, typename Device>
@@ -365,7 +365,8 @@ void ESolver_KS_PW<T, Device>::after_scf(UnitCell& ucell, const int istep, const
                                      this->pw_big,
                                      this->stp,
                                      this->Pgrid,
-                                     *this->inp_);
+                                     *this->inp_,
+                                     GlobalV::ofs_running);
 
     ModuleBase::timer::end("ESolver_KS_PW", "after_scf");
 }
@@ -463,7 +464,8 @@ void ESolver_KS_PW<T, Device>::after_all_runners(BaseCell& basecell)
                                         this->ppcell,
                                         this->solvent,
                                         this->Pgrid,
-                                        *this->inp_);
+                                        *this->inp_,
+                                        GlobalV::ofs_running);
 
     elecstate::teardown_estate_pw(this->pelec, this->vsep_cell);
 }

source/source_estate/init_scf.cpp

@@ -32,7 +32,7 @@ void init_scf(const UnitCell& ucell,
     pelec->pot->init_pot(pelec->charge);
 
     //! output the initial potential
-    ModuleIO::write_pot_init(ucell, pgrid, pelec, istep, out_dir, inp, PARAM.globalv.two_fermi);
+    ModuleIO::write_pot_init(ucell, pgrid, pelec, istep, out_dir, inp, PARAM.globalv.two_fermi, GlobalV::ofs_running);
 }
 
 } // namespace elecstate

source/source_estate/module_charge/chg_init.cpp

@@ -153,8 +153,8 @@ bool read_rho_file(Charge& chr,
  * @param rhopw [in] plane-wave basis for file decoding and Fourier transforms
  * @param cfg [in] file-reading configuration (suffix, dir, rank, logs)
  * @return true if a tau source was read; false if neither the restart binary
- *         nor the SPINX_TAU.cube file exists (caller falls back to the
- *         Thomas-Fermi tau init)
+ *         nor a tau.cube/taus1.cube (or legacy SPIN1_TAU.cube) file exists
+ *         (caller falls back to the Thomas-Fermi tau init)
  */
 bool read_kin_file(Charge& chr,
                    const ModulePW::PW_Basis& rhopw,
@@ -191,25 +191,59 @@ bool read_kin_file(Charge& chr,
         return true;
     }
 
-    // restart binary is absent; fall back to SPINX_TAU.cube files, probing on
-    // the parsing rank and broadcasting so all ranks agree to skip together.
+    // restart binary is absent; fall back to tau.cube/taus1.cube files, or to
+    // the legacy SPIN<n>_TAU.cube names written by 3.10-LTS and earlier.
+    // Probe on the parsing rank and broadcast so all ranks agree to skip
+    // together.
     std::stringstream ssc0;
-    ssc0 << readin_dir << "SPIN1_TAU.cube";
+    if (nspin == 1)
+    {
+        ssc0 << readin_dir << "tau.cube";
+    }
+    else
+    {
+        ssc0 << readin_dir << "taus1.cube";
+    }
+    const std::string legacy_name0 = readin_dir + "SPIN1_TAU.cube";
     bool cube_exists = false;
+    bool legacy_exists = false;
     if (rank == 0)
     {
         cube_exists = std::ifstream(ssc0.str()).good();
+        if (!cube_exists)
+        {
+            legacy_exists = std::ifstream(legacy_name0).good();
+        }
     }
     Parallel_Common::bcast_bool(cube_exists);
-    if (!cube_exists)
+    Parallel_Common::bcast_bool(legacy_exists);
+    if (!cube_exists && !legacy_exists)
     {
+        ofs_warning << " Kinetic energy density file " << ssc0.str()
+                    << " not found; tau will be initialized from the charge density." << std::endl;
         return false;
     }
+    if (legacy_exists)
+    {
+        ofs_warning << " Reading kinetic energy density from legacy file name " << legacy_name0
+                    << "; new outputs are named tau.cube (nspin=1) or taus<n>.cube." << std::endl;
+    }
 
     for (int is = 0; is < nspin; is++)
     {
         std::stringstream ssc;
-        ssc << readin_dir << "SPIN" << is + 1 << "_TAU.cube";
+        if (legacy_exists)
+        {
+            ssc << readin_dir << "SPIN" << is + 1 << "_TAU.cube";
+        }
+        else if (nspin == 1)
+        {
+            ssc << readin_dir << "tau.cube";
+        }
+        else
+        {
+            ssc << readin_dir << "taus" << is + 1 << ".cube";
+        }
         // mohan update 2012-02-10, sunliang update 2023-03-09
         ModuleIO::read_vdata_palgrid(pgrid,
                                      rank,

source/source_estate/module_charge/unittests/test_chg_extra.cpp

@@ -285,7 +285,6 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase4)
 
     EXPECT_EQ(output, expected_output);
     EXPECT_EQ(CE.get_rho_extr(), 3);
-    std::remove("./support/OLD2_SPIN1_CHG.cube");
 }
 
 TEST_F(ChargeExtraTest, UpdateAllDis)

source/source_estate/write_elecstat_pot.cpp

@@ -22,7 +22,8 @@ void write_elecstat_pot(
     const UnitCell* ucell,
     const double* v_eff,
     const surchem& solvent,
-    const int precision)
+    const int precision,
+    std::ofstream& ofs_running)
 {
     ModuleBase::TITLE("ModuleIO", "write_elecstat_pot");
     ModuleBase::timer::start("ModuleIO", "write_elecstat_pot");
@@ -84,7 +85,8 @@ void write_elecstat_pot(
                                   rho_basis->nxyz,
                                   rho_basis->nrxx,
                                   rho_basis->nplane,
-                                  rho_basis->startz_current);
+                                  rho_basis->startz_current,
+                                  ofs_running);
 
     //-------------------------------------------
     //! Write down the electrostatic potential
@@ -104,7 +106,9 @@ void write_elecstat_pot(
         precision,
         out_fermi,
         PARAM.globalv.two_fermi,
-        false);
+        false,
+        ofs_running,
+        "electrostatic potential");
 
     ModuleBase::timer::end("ModuleIO", "write_elecstat_pot");
     return;

source/source_estate/write_elecstat_pot.h

@@ -5,6 +5,7 @@
 #include "source_estate/module_charge/charge.h"
 #include "source_hamilt/module_surchem/surchem.h"
 
+#include <fstream>
 #include <string>
 
 namespace ModuleIO
@@ -33,7 +34,8 @@ void write_elecstat_pot(
     const UnitCell* ucell_,
     const double* v_eff_fixed,
     const surchem& solvent,
-    const int precision);
+    const int precision,
+    std::ofstream& ofs_running);
 
 } // namespace ModuleIO
 

source/source_estate/write_init.cpp

@@ -64,7 +64,8 @@ void ModuleIO::write_chg_init(
     const int istep,
     const std::string& out_dir,
     const Input_para& inp,
-    const bool two_fermi)
+    const bool two_fermi,
+    std::ofstream& ofs_running)
 {
     const int nspin = inp.nspin;
     assert(nspin == 1 || nspin == 2 || nspin == 4);
@@ -104,6 +105,8 @@ void ModuleIO::write_chg_init(
                     }
                 }
 
+                const std::string desc_chg = ModuleIO::make_data_desc("charge density", "magnetization density m", is, nspin);
+
                 ModuleIO::write_vdata_palgrid(para_grid,
                                               chr.rho[is],
                                               is,
@@ -115,7 +118,9 @@ void ModuleIO::write_chg_init(
                                               inp.out_chg[1],
                                               1,
                                               two_fermi,
-                                              false);
+                                              false,
+                                              ofs_running,
+                                              desc_chg);
             }
         }
     }
@@ -130,7 +135,8 @@ void ModuleIO::write_pot_init(
     const int istep,
     const std::string& out_dir,
     const Input_para& inp,
-    const bool two_fermi)
+    const bool two_fermi,
+    std::ofstream& ofs_running)
 {
     const int nspin = inp.nspin;
     assert(nspin == 1 || nspin == 2 || nspin == 4);
@@ -153,6 +159,8 @@ void ModuleIO::write_pot_init(
             {
                 std::string filename = gen_ini_filename("pot", out_dir, nspin, is, istep, include_geom_step);
 
+                const std::string desc_pot = ModuleIO::make_data_desc("effective potential", "effective magnetic field b", is, nspin);
+
                 ModuleIO::write_vdata_palgrid(para_grid,
                                               pelec->pot->get_eff_v(is),
                                               is,
@@ -164,7 +172,9 @@ void ModuleIO::write_pot_init(
                                               inp.out_pot[1],
                                               0,
                                               two_fermi,
-                                              false);
+                                              false,
+                                              ofs_running,
+                                              desc_pot);
             }
         }
     }