diff --git a/docs/advanced/elec_properties/hs_matrix.md b/docs/advanced/elec_properties/hs_matrix.md index 06cdc78786f..5adf14e3971 100644 --- a/docs/advanced/elec_properties/hs_matrix.md +++ b/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 = 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. diff --git a/source/source_esolver/esolver_dm2rho.cpp b/source/source_esolver/esolver_dm2rho.cpp index dda7b13d064..77b7825cc34 100644 --- a/source/source_esolver/esolver_dm2rho.cpp +++ b/source/source_esolver/esolver_dm2rho.cpp @@ -68,17 +68,29 @@ void ESolver_DM2rho::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 + // 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.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), @@ -86,7 +98,9 @@ void ESolver_DM2rho::runner(BaseCell& basecell, const int istep) 3, 1, PARAM.globalv.two_fermi, - false); + false, + GlobalV::ofs_running, + desc_chg); } ModuleBase::timer::end("ESolver_DM2rho", "runner"); diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 5a661ce4a35..773fd35b5f5 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/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; } diff --git a/source/source_esolver/esolver_gets.cpp b/source/source_esolver/esolver_gets.cpp index 5c27a14a18f..82000ddd105 100644 --- a/source/source_esolver/esolver_gets.cpp +++ b/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"); diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index f4826fca64d..eac97871445 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/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); } diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index f8c89b038ef..3236ea7d414 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -338,7 +338,7 @@ void ESolver_KS_LCAO::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::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) diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index 0a49659661a..0ed5011a4ee 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -333,7 +333,7 @@ void ESolver_KS_PW::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 @@ -365,7 +365,8 @@ void ESolver_KS_PW::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::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); } diff --git a/source/source_estate/init_scf.cpp b/source/source_estate/init_scf.cpp index eb9cf262262..45887b8ea69 100644 --- a/source/source_estate/init_scf.cpp +++ b/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 diff --git a/source/source_estate/module_charge/chg_init.cpp b/source/source_estate/module_charge/chg_init.cpp index 3960f467950..cc3f7dffc6e 100644 --- a/source/source_estate/module_charge/chg_init.cpp +++ b/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_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.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, diff --git a/source/source_estate/module_charge/unittests/test_chg_extra.cpp b/source/source_estate/module_charge/unittests/test_chg_extra.cpp index 962960fa726..7cac11d626f 100644 --- a/source/source_estate/module_charge/unittests/test_chg_extra.cpp +++ b/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) diff --git a/source/source_estate/write_elecstat_pot.cpp b/source/source_estate/write_elecstat_pot.cpp index 15a5363ed62..8ba9768a941 100644 --- a/source/source_estate/write_elecstat_pot.cpp +++ b/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; diff --git a/source/source_estate/write_elecstat_pot.h b/source/source_estate/write_elecstat_pot.h index bee575b95c3..8a1418e50ad 100644 --- a/source/source_estate/write_elecstat_pot.h +++ b/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 #include 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 diff --git a/source/source_estate/write_init.cpp b/source/source_estate/write_init.cpp index 58cdb534287..b58fe8ce986 100644 --- a/source/source_estate/write_init.cpp +++ b/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); } } } diff --git a/source/source_estate/write_init.h b/source/source_estate/write_init.h index 303f97e42ed..e475aa67cb8 100644 --- a/source/source_estate/write_init.h +++ b/source/source_estate/write_init.h @@ -18,6 +18,7 @@ #define WRITE_INIT_H #include +#include #include "source_io/module_parameter/input_parameter.h" #include "source_estate/module_charge/charge.h" #include "source_estate/fp_energy.h" @@ -53,7 +54,8 @@ void 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); // Write initial effective potential to cube file in real space. // Triggered when inp.out_pot[0] == 3. @@ -71,7 +73,8 @@ void 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); } diff --git a/source/source_io/module_chgpot/get_pchg_lcao.cpp b/source/source_io/module_chgpot/get_pchg_lcao.cpp index a7c1d307d37..3ececff4b1f 100644 --- a/source/source_io/module_chgpot/get_pchg_lcao.cpp +++ b/source/source_io/module_chgpot/get_pchg_lcao.cpp @@ -78,9 +78,9 @@ void Get_pchg_lcao::begin_gamma(const UnitCell& ucell, std::stringstream ssc; ssc << global_out_dir << "pchgi" << ib + 1 << "s" << is + 1 << ".cube"; - ofs_running << " Writing cube file " << ssc.str() << std::endl; + const std::string desc_pchg = ModuleIO::make_data_desc("partial charge", "partial magnetization m", is, nspin_); - ModuleIO::write_vdata_palgrid(pgrid, rho[is].data(), is, nspin_, 0, ssc.str(), 0.0, &ucell, precision, 0, false, false); + ModuleIO::write_vdata_palgrid(pgrid, rho[is].data(), is, nspin_, 0, ssc.str(), 0.0, &ucell, precision, 0, false, false, ofs_running, desc_pchg); } } } @@ -170,9 +170,9 @@ void Get_pchg_lcao::begin_k(const ModulePW::PW_Basis& rho_pw, std::stringstream ssc; ssc << global_out_dir << "pchgi" << ib + 1 << "s" << is + 1 << "k" << ik + 1 << ".cube"; - ofs_running << " Writing cube file " << ssc.str() << std::endl; + const std::string desc_pchg = ModuleIO::make_data_desc("partial charge", "partial magnetization m", is, nspin_); - ModuleIO::write_vdata_palgrid(pgrid, rho[is].data(), is, nspin_, 0, ssc.str(), 0.0, &ucell, precision, 0, false, false); + ModuleIO::write_vdata_palgrid(pgrid, rho[is].data(), is, nspin_, 0, ssc.str(), 0.0, &ucell, precision, 0, false, false, ofs_running, desc_pchg); } } } @@ -210,9 +210,9 @@ void Get_pchg_lcao::begin_k(const ModulePW::PW_Basis& rho_pw, std::stringstream ssc; ssc << global_out_dir << "pchgi" << ib + 1 << "s" << is + 1 << ".cube"; - ofs_running << " Writing cube file " << ssc.str() << std::endl; + const std::string desc_pchg = ModuleIO::make_data_desc("partial charge", "partial magnetization m", is, nspin_); - ModuleIO::write_vdata_palgrid(pgrid, rho[is].data(), is, nspin_, 0, ssc.str(), 0.0, &ucell, precision, 0, false, false); + ModuleIO::write_vdata_palgrid(pgrid, rho[is].data(), is, nspin_, 0, ssc.str(), 0.0, &ucell, precision, 0, false, false, ofs_running, desc_pchg); } } } diff --git a/source/source_io/module_chgpot/get_pchg_pw.cpp b/source/source_io/module_chgpot/get_pchg_pw.cpp index bb738965377..10561069885 100644 --- a/source/source_io/module_chgpot/get_pchg_pw.cpp +++ b/source/source_io/module_chgpot/get_pchg_pw.cpp @@ -116,7 +116,8 @@ void Get_pchg_pw::begin(UnitCell* ucell, const std::vector& out_pchg, const std::string& global_out_dir, const bool if_separate_k, - const bool noncolin) const + const bool noncolin, + std::ofstream& ofs_running) const { // Resolve global band ownership collectively before validating the selection. const Parallel::ParaBandOutput band_output(psi_.get_nbands(), global_nbands_, Parallel::make_band_world()); @@ -141,11 +142,11 @@ void Get_pchg_pw::begin(UnitCell* ucell, std::fill(work.becsum.begin(), work.becsum.end(), 0.0); if (if_separate_k) { - write_separate(band, *ucell, pgrid, kv, global_out_dir, noncolin, band_output, &work); + write_separate(band, *ucell, pgrid, kv, global_out_dir, noncolin, band_output, ofs_running, &work); } else { - write_summed(band, ucell, pgrid, kv, global_out_dir, noncolin, band_output, &work); + write_summed(band, ucell, pgrid, kv, global_out_dir, noncolin, band_output, ofs_running, &work); } } } @@ -236,6 +237,7 @@ void Get_pchg_pw::write_separate(const int band, const std::string& out_dir, const bool noncolin, const Parallel::ParaBandOutput& band_output, + std::ofstream& ofs_running, Workspace* work) const { // Collinear spin channels share the same physical k-point numbering in file names. @@ -258,7 +260,7 @@ void Get_pchg_pw::write_separate(const int band, const int component_end = work->is_spinor ? 4 : spin_index + 1; for (int component = component_begin; component < component_end; ++component) { - write_cube(band, component, k_number, ucell, pgrid, out_dir, true, work->density[component]); + write_cube(band, component, k_number, ucell, pgrid, out_dir, true, work->density[component], ofs_running); } } } @@ -271,6 +273,7 @@ void Get_pchg_pw::write_summed(const int band, const std::string& out_dir, const bool noncolin, const Parallel::ParaBandOutput& band_output, + std::ofstream& ofs_running, Workspace* work) const { for (int ik = 0; ik < kv.get_nks(); ++ik) @@ -288,7 +291,7 @@ void Get_pchg_pw::write_summed(const int band, symmetrize(ucell, work); for (int is = 0; is < nspin_; ++is) { - write_cube(band, is, 0, *ucell, pgrid, out_dir, false, work->density[is]); + write_cube(band, is, 0, *ucell, pgrid, out_dir, false, work->density[is], ofs_running); } } @@ -460,7 +463,8 @@ void Get_pchg_pw::write_cube(const int band, const Parallel_Grid& pgrid, const std::string& out_dir, const bool separate_k, - const std::vector& values) const + const std::vector& values, + std::ofstream& ofs_running) const { std::stringstream filename; filename << out_dir << "pchgi" << band + 1 << "s" << component + 1; @@ -469,7 +473,8 @@ void Get_pchg_pw::write_cube(const int band, filename << "k" << k_number; } filename << ".cube"; - ModuleIO::write_vdata_palgrid(pgrid, values.data(), component, nspin_, 0, filename.str(), 0.0, &ucell, 11, 0, false, separate_k); + const std::string desc_pchg = ModuleIO::make_data_desc("partial charge", "partial magnetization m", component, nspin_); + ModuleIO::write_vdata_palgrid(pgrid, values.data(), component, nspin_, 0, filename.str(), 0.0, &ucell, 11, 0, false, separate_k, ofs_running, desc_pchg); } // Explicit instantiation emits both precisions for each supported device from this .cpp file. diff --git a/source/source_io/module_chgpot/get_pchg_pw.h b/source/source_io/module_chgpot/get_pchg_pw.h index dd9860b16f5..9713c8bdbfd 100644 --- a/source/source_io/module_chgpot/get_pchg_pw.h +++ b/source/source_io/module_chgpot/get_pchg_pw.h @@ -46,7 +46,8 @@ class Get_pchg_pw const std::vector& out_pchg, const std::string& global_out_dir, const bool if_separate_k, - const bool noncolin) const; + const bool noncolin, + std::ofstream& ofs_running) const; private: const psi::Psi& psi_; @@ -76,6 +77,7 @@ class Get_pchg_pw const std::string& out_dir, const bool noncolin, const Parallel::ParaBandOutput& band_output, + std::ofstream& ofs_running, Workspace* work) const; void write_summed(const int band, UnitCell* ucell, @@ -84,6 +86,7 @@ class Get_pchg_pw const std::string& out_dir, const bool noncolin, const Parallel::ParaBandOutput& band_output, + std::ofstream& ofs_running, Workspace* work) const; void calc_density(const int spin_index, const double weight, const bool noncolin, const bool accumulate, Workspace* work) const; @@ -104,7 +107,8 @@ class Get_pchg_pw const Parallel_Grid& pgrid, const std::string& out_dir, const bool separate_k, - const std::vector& values) const; + const std::vector& values, + std::ofstream& ofs_running) const; }; } // namespace ModuleIO diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index ca8f5173358..a426bb8f3b2 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -26,7 +26,8 @@ void ctrl_output_fp(UnitCell& ucell, Charge& chr, surchem& solvent, Parallel_Grid& para_grid, - const int istep) + const int istep, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "ctrl_output_fp"); ModuleBase::timer::start("ModuleIO", "ctrl_output_fp"); @@ -80,6 +81,8 @@ void ctrl_output_fp(UnitCell& ucell, fn += spin_block + geom_block + ".cube"; + 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, @@ -91,7 +94,9 @@ void ctrl_output_fp(UnitCell& ucell, inp.out_chg[1], 1, PARAM.globalv.two_fermi, - false); + false, + ofs_running, + desc_chg); if (XC_Functional::get_ked_flag()) { @@ -99,6 +104,8 @@ void ctrl_output_fp(UnitCell& ucell, fn += spin_block + geom_block + ".cube"; + const std::string desc_tau = ModuleIO::make_data_desc("kinetic energy density", "kinetic energy density magnetization m", is, nspin); + ModuleIO::write_vdata_palgrid(para_grid, chr.kin_r[is], is, @@ -110,7 +117,9 @@ void ctrl_output_fp(UnitCell& ucell, 11, // default precision 1, // default out_fermi PARAM.globalv.two_fermi, - false); + false, + ofs_running, + desc_tau); } } } @@ -134,6 +143,8 @@ void ctrl_output_fp(UnitCell& ucell, fn += spin_block + geom_block + ".cube"; + 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, @@ -145,7 +156,9 @@ void ctrl_output_fp(UnitCell& ucell, inp.out_pot[1], // precision 0, // out_fermi PARAM.globalv.two_fermi, - false); + false, + ofs_running, + desc_pot); } } else if (inp.out_pot[0] == 2 && should_output) @@ -165,7 +178,8 @@ void ctrl_output_fp(UnitCell& ucell, &(ucell), pelec->pot->get_fixed_v(), solvent, - inp.out_pot[1]); + inp.out_pot[1], + ofs_running); } // 6) write ELF @@ -181,16 +195,17 @@ void ctrl_output_fp(UnitCell& ucell, std::string out_dir = PARAM.globalv.global_out_dir; ModuleIO::write_elf(out_dir, - istep, - nspin, - chr.rho, - chr.kin_r, - pw_rhod, - para_grid, - &(ucell), - inp.out_elf[1], - geom_block, - PARAM.globalv.two_fermi); + istep, + nspin, + chr.rho, + chr.kin_r, + pw_rhod, + para_grid, + &(ucell), + inp.out_elf[1], + geom_block, + PARAM.globalv.two_fermi, + ofs_running); } #ifdef __LIBXC @@ -215,7 +230,7 @@ void ctrl_output_fp(UnitCell& ucell, { std::stringstream ss_dipole; ss_dipole << global_out_dir << "dipole_s" << is + 1 << ".txt"; - ModuleIO::write_dipole(ucell, chr.rho[is], pw_rhod, istep, ss_dipole.str(), GlobalV::ofs_running); + ModuleIO::write_dipole(ucell, chr.rho[is], pw_rhod, istep, ss_dipole.str(), ofs_running); } } diff --git a/source/source_io/module_ctrl/ctrl_output_fp.h b/source/source_io/module_ctrl/ctrl_output_fp.h index 10c4c755620..59352509a3f 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.h +++ b/source/source_io/module_ctrl/ctrl_output_fp.h @@ -1,6 +1,8 @@ #ifndef CTRL_OUTPUT_FP_H #define CTRL_OUTPUT_FP_H +#include + struct Input_para; class UnitCell; class Charge; @@ -29,7 +31,8 @@ void ctrl_output_fp(UnitCell& ucell, Charge& chr, surchem& solvent, Parallel_Grid& para_grid, - const int istep); + const int istep, + std::ofstream& ofs_running); } #endif diff --git a/source/source_io/module_ctrl/ctrl_output_pw.cpp b/source/source_io/module_ctrl/ctrl_output_pw.cpp index 42f14ab6393..6fa96d9e1ef 100644 --- a/source/source_io/module_ctrl/ctrl_output_pw.cpp +++ b/source/source_io/module_ctrl/ctrl_output_pw.cpp @@ -26,7 +26,8 @@ void ModuleIO::ctrl_iter_pw(const int istep, Setup_Psi_pw& stp, const K_Vectors& kv, const ModulePW::PW_Basis_K* pw_wfc, - const Input_para& inp) + const Input_para& inp, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "ctrl_iter_pw"); ModuleBase::timer::start("ModuleIO", "ctrl_iter_pw"); @@ -76,7 +77,7 @@ void ModuleIO::ctrl_iter_pw(const int istep, stp.psi_cpu[0], kv, pw_wfc, - GlobalV::ofs_running); + ofs_running); } ModuleBase::timer::end("ModuleIO", "ctrl_iter_pw"); @@ -96,7 +97,8 @@ void ModuleIO::ctrl_scf_pw(const int istep, const ModulePW::PW_Basis_Big* pw_big, Setup_Psi_pw& stp, const Parallel_Grid& para_grid, - const Input_para& inp) + const Input_para& inp, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "ctrl_scf_pw"); ModuleBase::timer::start("ModuleIO", "ctrl_scf_pw"); @@ -167,7 +169,7 @@ void ModuleIO::ctrl_scf_pw(const int istep, inp.out_app_flag, inp.bndpar, PARAM.globalv.global_out_dir, - GlobalV::ofs_running); + ofs_running); } } @@ -180,14 +182,14 @@ void ModuleIO::ctrl_scf_pw(const int istep, // Use the solver's native wavefunction precision for FFT and projector overlaps. ModuleIO::Get_pchg_pw output(*stp.template get_psi_t(), *pw_wfc, *pw_rho, *pw_rhod, ppcell, inp.nspin, inp.nbands); - output.begin(&ucell, para_grid, kv, inp.out_pchg, PARAM.globalv.global_out_dir, inp.if_separate_k, inp.noncolin); + output.begin(&ucell, para_grid, kv, inp.out_pchg, PARAM.globalv.global_out_dir, inp.if_separate_k, inp.noncolin, ofs_running); } if (inp.out_wfc_norm.size() > 0 || inp.out_wfc_re_im.size() > 0) { // Use the solver's native wavefunction precision for FFT. ModuleIO::Get_wf_pw output(*stp.template get_psi_t(), *pw_wfc, *pw_rho, *pw_rhod, inp.nspin, inp.nbands); - output.begin(ucell, para_grid, kv, inp.out_wfc_norm, inp.out_wfc_re_im, PARAM.globalv.global_out_dir); + output.begin(ucell, para_grid, kv, inp.out_wfc_norm, inp.out_wfc_re_im, PARAM.globalv.global_out_dir, ofs_running); } //------------------------------------------------------------------ @@ -232,7 +234,7 @@ void ModuleIO::ctrl_scf_pw(const int istep, { spinconstrain::SpinConstrain>& sc = spinconstrain::SpinConstrain>::getScInstance(); spinconstrain::pw::cal_mi_pw(sc.state_, sc.psi, sc.pelec); - spinconstrain::print_Mag_Force(sc, GlobalV::ofs_running); + spinconstrain::print_Mag_Force(sc, ofs_running); } //------------------------------------------------------------------ @@ -263,7 +265,8 @@ void ModuleIO::ctrl_runner_pw(UnitCell& ucell, pseudopot_cell_vnl& ppcell, surchem& solvent, Parallel_Grid& para_grid, - const Input_para& inp) + const Input_para& inp, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "ctrl_runner_pw"); ModuleBase::timer::start("ModuleIO", "ctrl_runner_pw"); @@ -295,7 +298,8 @@ void ModuleIO::ctrl_runner_pw(UnitCell& ucell, PARAM.globalv.dos_setemin, inp.dos_emin_ev, inp.dos_sigma, - inp.ldos_line); + inp.ldos_line, + ofs_running); } //---------------------------------------------------------- @@ -317,7 +321,7 @@ void ModuleIO::ctrl_runner_pw(UnitCell& ucell, Numerical_Basis numerical_basis; numerical_basis.output_overlap(stp.psi_cpu[0], sf, kv, pw_wfc, ucell, i); } - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "BASIS OVERLAP (Q and S) GENERATION."); + ModuleBase::GlobalFunc::DONE(ofs_running, "BASIS OVERLAP (Q and S) GENERATION."); } } @@ -357,7 +361,7 @@ void ModuleIO::ctrl_runner_pw(UnitCell& ucell, inp.of_ml_kernel_file, ucell.omega, pw_rho, - GlobalV::ofs_running); + ofs_running); write_mlkedf_desc.generateTrainData_KS(PARAM.globalv.global_mlkedf_descriptor_dir, stp.template get_psi_t(), @@ -386,7 +390,8 @@ template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_CPU const ModulePW::PW_Basis_Big* pw_big, Setup_Psi_pw& stp, const Parallel_Grid& para_grid, - const Input_para& inp); + const Input_para& inp, + std::ofstream& ofs_running); // complex + CPU template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_CPU>(const int nstep, @@ -401,7 +406,8 @@ template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_CP const ModulePW::PW_Basis_Big* pw_big, Setup_Psi_pw& stp, const Parallel_Grid& para_grid, - const Input_para& inp); + const Input_para& inp, + std::ofstream& ofs_running); #if ((defined __CUDA) || (defined __ROCM)) // complex + GPU @@ -417,7 +423,8 @@ template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_GPU const ModulePW::PW_Basis_Big* pw_big, Setup_Psi_pw& stp, const Parallel_Grid& para_grid, - const Input_para& inp); + const Input_para& inp, + std::ofstream& ofs_running); // complex + GPU template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_GPU>(const int nstep, @@ -432,7 +439,8 @@ template void ModuleIO::ctrl_scf_pw, base_device::DEVICE_GP const ModulePW::PW_Basis_Big* pw_big, Setup_Psi_pw& stp, const Parallel_Grid& para_grid, - const Input_para& inp); + const Input_para& inp, + std::ofstream& ofs_running); #endif // complex + CPU @@ -448,7 +456,8 @@ template void ModuleIO::ctrl_runner_pw, base_device::DEVICE_ pseudopot_cell_vnl& ppcell, surchem& solvent, Parallel_Grid& para_grid, - const Input_para& inp); + const Input_para& inp, + std::ofstream& ofs_running); // complex + CPU template void ModuleIO::ctrl_runner_pw, base_device::DEVICE_CPU>(UnitCell& ucell, @@ -463,7 +472,8 @@ template void ModuleIO::ctrl_runner_pw, base_device::DEVICE pseudopot_cell_vnl& ppcell, surchem& solvent, Parallel_Grid& para_grid, - const Input_para& inp); + const Input_para& inp, + std::ofstream& ofs_running); #if ((defined __CUDA) || (defined __ROCM)) // complex + GPU @@ -479,7 +489,8 @@ template void ModuleIO::ctrl_runner_pw, base_device::DEVICE_ pseudopot_cell_vnl& ppcell, surchem& solvent, Parallel_Grid& para_grid, - const Input_para& inp); + const Input_para& inp, + std::ofstream& ofs_running); // complex + GPU template void ModuleIO::ctrl_runner_pw, base_device::DEVICE_GPU>(UnitCell& ucell, @@ -494,5 +505,6 @@ template void ModuleIO::ctrl_runner_pw, base_device::DEVICE pseudopot_cell_vnl& ppcell, surchem& solvent, Parallel_Grid& para_grid, - const Input_para& inp); + const Input_para& inp, + std::ofstream& ofs_running); #endif diff --git a/source/source_io/module_ctrl/ctrl_output_pw.h b/source/source_io/module_ctrl/ctrl_output_pw.h index 62518166a11..9491afbf787 100644 --- a/source/source_io/module_ctrl/ctrl_output_pw.h +++ b/source/source_io/module_ctrl/ctrl_output_pw.h @@ -6,6 +6,8 @@ #include "source_psi/psi.h" // define psi #include "source_psi/setup_psi_pw.h" // use Setup_Psi class +#include + class pseudopot_cell_vnl; namespace ModuleIO @@ -18,7 +20,8 @@ void ctrl_iter_pw(const int istep, Setup_Psi_pw& stp, const K_Vectors& kv, const ModulePW::PW_Basis_K* pw_wfc, - const Input_para& inp); + const Input_para& inp, + std::ofstream& ofs_running); // print out information in 'after_scf' in ESolver_KS_PW template @@ -34,7 +37,8 @@ void ctrl_scf_pw(const int istep, const ModulePW::PW_Basis_Big* pw_big, Setup_Psi_pw& stp, const Parallel_Grid& para_grid, - const Input_para& inp); + const Input_para& inp, + std::ofstream& ofs_running); // print out information in 'after_all_runners' in ESolver_KS_PW template @@ -50,7 +54,8 @@ void ctrl_runner_pw(UnitCell& ucell, pseudopot_cell_vnl& ppcell, surchem& solvent, Parallel_Grid& para_grid, - const Input_para& inp); + const Input_para& inp, + std::ofstream& ofs_running); // print out information in 'after_all_runners' in ESolver_KS_PW (runtime version) template @@ -67,7 +72,8 @@ void ctrl_runner_pw(UnitCell& ucell, surchem& solvent, const base_device::DeviceContext* ctx, Parallel_Grid& para_grid, - const Input_para& inp); + const Input_para& inp, + std::ofstream& ofs_running); } // namespace ModuleIO #endif diff --git a/source/source_io/module_ctrl/ctrl_runner_lcao.cpp b/source/source_io/module_ctrl/ctrl_runner_lcao.cpp index 4966af3924b..5ceec5f9165 100644 --- a/source/source_io/module_ctrl/ctrl_runner_lcao.cpp +++ b/source/source_io/module_ctrl/ctrl_runner_lcao.cpp @@ -37,7 +37,8 @@ void ctrl_runner_lcao(UnitCell& ucell, // unitcell ModuleBase::matrix &vloc, // local pseudopotential Exx_NAO &exx_nao, const Exx_Info& exx_info, - surchem &solvent) // solvent model + surchem &solvent, // solvent model + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "ctrl_runner_lcao"); ModuleBase::timer::start("ModuleIO", "ctrl_runner_lcao"); @@ -54,7 +55,7 @@ void ctrl_runner_lcao(UnitCell& ucell, // unitcell ModuleIO::Cal_ldos::cal_ldos_lcao(pelec->eferm, chr, dmat, kv, pelec->ekb, pelec->wg, psi[0], pgrid, pv, gd, ucell, inp.stm_bias, inp.nspin, PARAM.globalv.global_out_dir, - PARAM.globalv.two_fermi, inp.out_ldos[1]); + PARAM.globalv.two_fermi, inp.out_ldos[1], ofs_running); } // 3) print out exchange-correlation potential @@ -128,7 +129,8 @@ void ctrl_runner_lcao(UnitCell& ucell, // unitcell hexxc_ptr #endif , - sparse_thr + sparse_thr, + ofs_running ); } @@ -188,7 +190,8 @@ template void ctrl_runner_lcao(UnitCell& ucell, // unitcell ModuleBase::matrix &vloc, // local pseudopotential Exx_NAO &exx_nao, const Exx_Info& exx_info, - surchem &solvent); // solvent model + surchem &solvent, // solvent model + std::ofstream& ofs_running); // TK: complex TR: double template void ctrl_runner_lcao, double>(UnitCell& ucell, // unitcell @@ -210,7 +213,8 @@ template void ctrl_runner_lcao, double>(UnitCell& ucell, ModuleBase::matrix &vloc, // local pseudopotential Exx_NAO> &exx_nao, const Exx_Info& exx_info, - surchem &solvent); // solvent model + surchem &solvent, // solvent model + std::ofstream& ofs_running); // TK: complex TR: complex template void ctrl_runner_lcao, std::complex>(UnitCell& ucell, // unitcell @@ -232,6 +236,7 @@ template void ctrl_runner_lcao, std::complex>(UnitC ModuleBase::matrix &vloc, // local pseudopotential Exx_NAO> &exx_nao, const Exx_Info& exx_info, - surchem &solvent); // solvent model + surchem &solvent, // solvent model + std::ofstream& ofs_running); } // end namespace diff --git a/source/source_io/module_ctrl/ctrl_runner_lcao.h b/source/source_io/module_ctrl/ctrl_runner_lcao.h index 6f51f9ae273..96434852343 100644 --- a/source/source_io/module_ctrl/ctrl_runner_lcao.h +++ b/source/source_io/module_ctrl/ctrl_runner_lcao.h @@ -10,6 +10,8 @@ #include "source_lcao/setup_exx.h" // for exx, mohan add 20251018 #include "source_lcao/allocate_dm.h" // for density matrix, mohan add 20251103 +#include + namespace ModuleIO { @@ -33,7 +35,8 @@ void ctrl_runner_lcao(UnitCell& ucell, // unitcell ModuleBase::matrix &vloc, // local pseudopotential Exx_NAO &exx_nao, const Exx_Info& exx_info, - surchem &solvent); // solvent model + surchem &solvent, // solvent model + std::ofstream& ofs_running); } diff --git a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp index e3f6bd6dd1f..dddf2576c24 100644 --- a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp +++ b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp @@ -9,6 +9,7 @@ #include "source_lcao/hamilt_lcao.h" // use hamilt::HamiltLCAO #include +#include // functions #include "../module_unk/berryphase.h" // use berryphase @@ -94,7 +95,8 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, const Exx_Info& exx_info, const bool conv_esolver, const bool scf_nmax_flag, - const int istep) + const int istep, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "ctrl_scf_lcao"); ModuleBase::timer::start("ModuleIO", "ctrl_scf_lcao"); @@ -167,7 +169,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, PARAM.globalv.two_fermi, inp.bndpar, PARAM.globalv.global_out_dir, - GlobalV::ofs_running); + ofs_running); } //------------------------------------------------------------------ @@ -178,7 +180,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, const int precision = inp.out_dmr[1]; ModuleIO::write_dmr(dm->get_dmr_vec(), &ucell, precision, pv, out_app_flag, - ucell.get_iat2iwt(), ucell.nat, istep); + ucell.get_iat2iwt(), ucell.nat, istep, ofs_running); } //------------------------------------------------------------------ @@ -193,7 +195,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, } const int precision = inp.out_dmk[1]; - ModuleIO::write_dmk(dm->get_dmk_vec(), kv, precision, efermis, &(ucell), pv, global_out_dir, istep); + ModuleIO::write_dmk(dm->get_dmk_vec(), kv, precision, efermis, &(ucell), pv, global_out_dir, istep, ofs_running); } //------------------------------------------------------------------ @@ -219,7 +221,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, PARAM.globalv.nlocal, PARAM.inp.ks_solver, GlobalV::DRANK, - GlobalV::ofs_running); + ofs_running); } //------------------------------------------------------------------ @@ -264,7 +266,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, -1, // -1 when called in after scf true, // no used when after scf GlobalV::MY_RANK, - GlobalV::ofs_running); + ofs_running); #endif //------------------------------------------------------------------ @@ -278,7 +280,8 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, ModuleIO::write_hsr(hr_vec, sr, &ucell, inp.out_hsr[0], precision, pv, out_app_flag, gamma_only, ucell.get_iat2iwt(), ucell.nat, istep, - PARAM.globalv.global_out_dir); + PARAM.globalv.global_out_dir, pelec->eferm, + ofs_running); } //------------------------------------------------------------------ @@ -339,7 +342,8 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, gd, kv, p_ham_tk, - &dftu); + &dftu, + ofs_running); //------------------------------------------------------------------ //! 7c) Output atomic dH components (dT/dτ, dV^NL/dτ, dV^L/dτ, dV^H/dτ, dV^XC/dτ), only for nspin =1, 2 now @@ -396,6 +400,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, dh_params.istep = istep; dh_params.gamma_only = gamma_only; dh_params.append = out_app_flag; + dh_params.ofs_running = &ofs_running; if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) { // per-spin DM (1-indexed): nspin=1 -> {spin0}, nspin=2 -> {spin-up, spin-down}. @@ -452,6 +457,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, h_params.global_matrix_dir = PARAM.globalv.global_matrix_dir; h_params.ks_solver = PARAM.inp.ks_solver; h_params.drank = GlobalV::DRANK; + h_params.ofs_running = &ofs_running; if (inp.out_mat_h_t[0]) { ModuleIO::write_h_t(h_params); @@ -526,6 +532,11 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, pv, GlobalV::DRANK, PARAM.inp.ks_solver); + + if (GlobalV::DRANK == 0) + { + ofs_running << " Write T(k) matrix in NAO basis to file: " << t_fn << std::endl; + } } delete ekinetic; @@ -545,7 +556,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, PARAM.globalv.search_pbc, PARAM.inp.out_level, PARAM.globalv.gamma_only_local, - &GlobalV::ofs_running, + &ofs_running, GlobalV::MY_RANK); mylcalculator.calculate(inp.suffix, global_out_dir, ucell, inp.out_mat_l[1], GlobalV::MY_RANK, istep_in); } @@ -574,8 +585,8 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, { spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); sc.cal_mi_lcao(istep); - spinconstrain::print_Mi(sc, GlobalV::ofs_running); - spinconstrain::print_Mag_Force(sc, GlobalV::ofs_running); + spinconstrain::print_Mi(sc, ofs_running); + spinconstrain::print_Mag_Force(sc, ofs_running); } //------------------------------------------------------------------ @@ -717,7 +728,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, inp.orbital_dir, &ucell, kv.kvec_d, - GlobalV::ofs_running, + ofs_running, GlobalV::MY_RANK, GlobalV::NPROC); tqo.calculate(); @@ -787,7 +798,8 @@ template void ModuleIO::ctrl_scf_lcao( const Exx_Info& exx_info, const bool conv_esolver, const bool scf_nmax_flag, - const int istep); + const int istep, + std::ofstream& ofs_running); // For multiple k-points template void ModuleIO::ctrl_scf_lcao, double>( @@ -816,7 +828,8 @@ template void ModuleIO::ctrl_scf_lcao, double>( const Exx_Info& exx_info, const bool conv_esolver, const bool scf_nmax_flag, - const int istep); + const int istep, + std::ofstream& ofs_running); template void ModuleIO::ctrl_scf_lcao, std::complex>( UnitCell& ucell, @@ -844,4 +857,5 @@ template void ModuleIO::ctrl_scf_lcao, std::complex const Exx_Info& exx_info, const bool conv_esolver, const bool scf_nmax_flag, - const int istep); + const int istep, + std::ofstream& ofs_running); diff --git a/source/source_io/module_ctrl/ctrl_scf_lcao.h b/source/source_io/module_ctrl/ctrl_scf_lcao.h index 87f1aa9dd53..013a25901da 100644 --- a/source/source_io/module_ctrl/ctrl_scf_lcao.h +++ b/source/source_io/module_ctrl/ctrl_scf_lcao.h @@ -17,6 +17,7 @@ #include "source_pw/module_pwdft/stru_fac.h" // use Structure_Factor #include +#include namespace ModuleIO { @@ -47,6 +48,7 @@ void ctrl_scf_lcao(UnitCell& ucell, const Exx_Info& exx_info, const bool conv_esolver, const bool scf_nmax_flag, - const int istep); + const int istep, + std::ofstream& ofs_running); } // namespace ModuleIO #endif diff --git a/source/source_io/module_dhs/write_dh.cpp b/source/source_io/module_dhs/write_dh.cpp index 3fca7a9fa1e..9a61cf3dee9 100644 --- a/source/source_io/module_dhs/write_dh.cpp +++ b/source/source_io/module_dhs/write_dh.cpp @@ -4,6 +4,7 @@ #include "source_base/timer.h" #include "source_io/module_hs/hs_dense_io.h" #include "source_io/module_hs/hsr_writer.h" +#include "source_io/module_output/spin_tag.h" #include "source_cell/ucell_io.h" #include "source_io/module_parameter/parameter.h" #include "source_hamilt/module_hcontainer/hcontainer_funcs.h" @@ -26,6 +27,7 @@ void write_dh_perI(WriteDHParams& params, const std::string& rprefix, const std::string& kprefix, const std::string& label, + const std::string& term_name, std::array*>, 3>& g, const std::vector& atom_filter) { @@ -78,13 +80,21 @@ void write_dh_perI(WriteDHParams& params, #endif { std::string fr = r_dir + ModuleIO::dhr_gen_fname(rprefix + tag, ispin, params.append, params.istep); + // dH/dR is not a Hamiltonian: no Fermi energy annotation in the header + const double no_efermi = 0.0; #ifdef __MPI ModuleIO::write_hcontainer_csr( - fr, &ucell, 8, &hR_s, params.istep, ispin, nspin, label, ""); + fr, &ucell, 8, &hR_s, params.istep, ispin, nspin, label, "", no_efermi, false); #else ModuleIO::write_hcontainer_csr( - fr, &ucell, 8, hR, params.istep, ispin, nspin, label, ""); + fr, &ucell, 8, hR, params.istep, ispin, nspin, label, "", no_efermi, false); #endif + if (params.ofs_running != nullptr) + { + const std::string spin_tag = ModuleIO::make_spin_tag(ispin, nspin); + *params.ofs_running << " Write dH/dR (" << term_name << " term)" << spin_tag + << " matrix in NAO basis to file: " << fr << std::endl; + } } } @@ -157,12 +167,6 @@ void write_dH_components(WriteDHParams& params, const Exx_Info& exx_info) } #endif - GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - GlobalV::ofs_running << " | |" << std::endl; - GlobalV::ofs_running << " | #Print out dH/dR components# |" << std::endl; - GlobalV::ofs_running << " | |" << std::endl; - GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - if (PARAM.inp.out_mat_dh[0]) { write_dH_sum(params); diff --git a/source/source_io/module_dhs/write_dh.h b/source/source_io/module_dhs/write_dh.h index 1258daad4fb..e4076efb50f 100644 --- a/source/source_io/module_dhs/write_dh.h +++ b/source/source_io/module_dhs/write_dh.h @@ -11,6 +11,7 @@ #include #include +#include #include ///for lack of make_unique in c++11 @@ -51,6 +52,7 @@ struct WriteDHParams // terms: V^H needs the total density (sum over spins), V^XC the spin-resolved densities. std::vector*> dmR; const Charge* chg = nullptr; // ground-state charge for XC Hellmann-Feynman (FDM) + std::ofstream* ofs_running = nullptr; // running log for the "Write ... to file" messages }; // Returns 0-based atom indices to output (converted from the 1-based user-facing values stored at param[2+]); @@ -74,6 +76,7 @@ void write_dh_perI(WriteDHParams& params, const std::string& rprefix, const std::string& kprefix, const std::string& label, + const std::string& term_name, std::array*>, 3>& g, const std::vector& atom_filter = {}); diff --git a/source/source_io/module_dhs/write_dh_terms.cpp b/source/source_io/module_dhs/write_dh_terms.cpp index 77dbf2a97f5..f3d78f6ed69 100644 --- a/source/source_io/module_dhs/write_dh_terms.cpp +++ b/source/source_io/module_dhs/write_dh_terms.cpp @@ -175,6 +175,7 @@ bool write_dH_veff_term(WriteDHParams& params, const std::string& rprefix, const std::string& kprefix, const std::string& label, + const std::string& term_name, const std::vector& atom_filter = {}) { const UnitCell& ucell = *params.ucell; @@ -194,7 +195,7 @@ bool write_dH_veff_term(WriteDHParams& params, fill_dH_veff(params, pot, hf_type, ispin, c); - ModuleIO::write_dh_perI(params, ispin, rprefix, kprefix, label, c.g, atom_filter); + ModuleIO::write_dh_perI(params, ispin, rprefix, kprefix, label, term_name, c.g, atom_filter); } return true; } @@ -218,7 +219,7 @@ bool write_dH_t(WriteDHParams& params) fill_dH_t(params, c); - ModuleIO::write_dh_perI(params, ispin, "dtr", "dtk", "dT", c.g, af_t); + ModuleIO::write_dh_perI(params, ispin, "dtr", "dtk", "dT", "kinetic", c.g, af_t); } ModuleBase::timer::end("ModuleIO", "write_dH_t"); @@ -241,7 +242,7 @@ bool write_dH_vnl(WriteDHParams& params) fill_dH_vnl(params, c); - ModuleIO::write_dh_perI(params, ispin, "dvnlr", "dvnlk", "dV^NL", c.g, af_vnl); + ModuleIO::write_dh_perI(params, ispin, "dvnlr", "dvnlk", "dV^NL", "nonlocal", c.g, af_vnl); } ModuleBase::timer::end("ModuleIO", "write_dH_vnl"); @@ -254,7 +255,7 @@ bool write_dH_vl(WriteDHParams& params) ModuleBase::timer::start("ModuleIO", "write_dH_vl"); const std::vector af_vl = dh_atom_filter(PARAM.inp.out_mat_dh_vl); - const bool ok = write_dH_veff_term(params, params.pot_vl, "vl", "dvlr", "dvlk", "dV^L", af_vl); + const bool ok = write_dH_veff_term(params, params.pot_vl, "vl", "dvlr", "dvlk", "dV^L", "local", af_vl); ModuleBase::timer::end("ModuleIO", "write_dH_vl"); return ok; @@ -266,7 +267,7 @@ bool write_dH_vh(WriteDHParams& params) ModuleBase::timer::start("ModuleIO", "write_dH_vh"); const std::vector af_vh = dh_atom_filter(PARAM.inp.out_mat_dh_vh); - const bool ok = write_dH_veff_term(params, params.pot_vh, "hartree", "dvhr", "dvhk", "dV^H", af_vh); + const bool ok = write_dH_veff_term(params, params.pot_vh, "hartree", "dvhr", "dvhk", "dV^H", "hartree", af_vh); ModuleBase::timer::end("ModuleIO", "write_dH_vh"); return ok; @@ -278,7 +279,7 @@ bool write_dH_vh_pulay(WriteDHParams& params) ModuleBase::timer::start("ModuleIO", "write_dH_vh_pulay"); const std::vector af_vh_pulay = dh_atom_filter(PARAM.inp.out_mat_dh_vh); - const bool ok = write_dH_veff_term(params, params.pot_vh, "none", "dvhr_pulay_", "dvhk_pulay_", "dV^H (Pulay)", af_vh_pulay); + const bool ok = write_dH_veff_term(params, params.pot_vh, "none", "dvhr_pulay_", "dvhk_pulay_", "dV^H (Pulay)", "hartree Pulay", af_vh_pulay); ModuleBase::timer::end("ModuleIO", "write_dH_vh_pulay"); return ok; @@ -292,7 +293,7 @@ bool write_dH_vxc(WriteDHParams& params) const std::vector af_vxc = dh_atom_filter(PARAM.inp.out_mat_dh_vxc); const bool ok = write_dH_veff_term(params, params.pot_vxc, params.chg ? "xc" : "none", - "dvxcr", "dvxck", "dV^XC", af_vxc); + "dvxcr", "dvxck", "dV^XC", "xc", af_vxc); ModuleBase::timer::end("ModuleIO", "write_dH_vxc"); return ok; @@ -304,7 +305,7 @@ bool write_dH_vxc_pulay(WriteDHParams& params) ModuleBase::timer::start("ModuleIO", "write_dH_vxc_pulay"); const std::vector af_vxc_pulay = dh_atom_filter(PARAM.inp.out_mat_dh_vxc); - const bool ok = write_dH_veff_term(params, params.pot_vxc, "none", "dvxcr_pulay_", "dvxck_pulay_", "dV^XC (Pulay)", af_vxc_pulay); + const bool ok = write_dH_veff_term(params, params.pot_vxc, "none", "dvxcr_pulay_", "dvxck_pulay_", "dV^XC (Pulay)", "xc Pulay", af_vxc_pulay); ModuleBase::timer::end("ModuleIO", "write_dH_vxc_pulay"); return ok; @@ -372,7 +373,7 @@ bool write_dH_sum(WriteDHParams& params) accumulate(c); } - ModuleIO::write_dh_perI(params, ispin, "dhr", "dhk", "dH", sum.g, dh_atom_filter(PARAM.inp.out_mat_dh)); + ModuleIO::write_dh_perI(params, ispin, "dhr", "dhk", "dH", "total", sum.g, dh_atom_filter(PARAM.inp.out_mat_dh)); } ModuleBase::timer::end("ModuleIO", "write_dH_sum"); diff --git a/source/source_io/module_dm/test/write_dmk_test.cpp b/source/source_io/module_dm/test/write_dmk_test.cpp index 42a4c181f9c..4d6a2a8f1bc 100644 --- a/source/source_io/module_dm/test/write_dmk_test.cpp +++ b/source/source_io/module_dm/test/write_dmk_test.cpp @@ -164,8 +164,9 @@ TEST(DMKTest,WriteDMK) { kv.kc_done = true; kv.kd_done = true; - ModuleIO::write_dmk(dmk, kv, 3, efs, ucell, pv, out_dir, istep); - ModuleIO::write_dmk(dmk_multik, kv, 3, efs, ucell, pv, out_dir, istep); + std::ofstream ofs_log; + ModuleIO::write_dmk(dmk, kv, 3, efs, ucell, pv, out_dir, istep, ofs_log); + ModuleIO::write_dmk(dmk_multik, kv, 3, efs, ucell, pv, out_dir, istep, ofs_log); std::ifstream ifs; diff --git a/source/source_io/module_dm/write_dmk.cpp b/source/source_io/module_dm/write_dmk.cpp index eac2ca021d5..beec822652a 100644 --- a/source/source_io/module_dm/write_dmk.cpp +++ b/source/source_io/module_dm/write_dmk.cpp @@ -1,4 +1,5 @@ #include "source_io/module_dm/write_dmk.h" +#include "source_io/module_output/spin_tag.h" #include "source_base/parallel_common.h" #include "source_base/module_external/scalapack_connector.h" @@ -186,7 +187,8 @@ void ModuleIO::write_dmk(const std::vector>& dmk, const UnitCell* ucell, const Parallel_2D& pv, const std::string& dmk_dir, - const int istep) + const int istep, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "write_dmk"); ModuleBase::timer::start("ModuleIO", "write_dmk"); @@ -250,7 +252,8 @@ void ModuleIO::write_dmk(const std::vector>& dmk, } else { - //std::cout << " Write the density matrix to file " << fn << std::endl; + const std::string spin_tag = ModuleIO::make_spin_tag(ispin, nspin); + ofs_running << " Write DM(k)" << spin_tag << " matrix in NAO basis to file: " << fn << std::endl; } @@ -343,7 +346,8 @@ template void ModuleIO::write_dmk(const std::vector> const UnitCell* ucell, const Parallel_2D& pv, const std::string& dmk_dir, - const int istep); + const int istep, + std::ofstream& ofs_running); template void ModuleIO::write_dmk>(const std::vector>>& dmk, const K_Vectors &kv, @@ -352,5 +356,6 @@ template void ModuleIO::write_dmk>(const std::vector #include #include @@ -76,7 +77,8 @@ void write_dmk(const std::vector>& dmk, const UnitCell* ucell, const Parallel_2D& pv, const std::string& dmk_dir, - const int istep); + const int istep, + std::ofstream& ofs_running); } // namespace ModuleIO diff --git a/source/source_io/module_dm/write_dmr.cpp b/source/source_io/module_dm/write_dmr.cpp index 33306a43e53..6b95e726a31 100644 --- a/source/source_io/module_dm/write_dmr.cpp +++ b/source/source_io/module_dm/write_dmr.cpp @@ -1,4 +1,5 @@ #include "write_dmr.h" +#include "source_io/module_output/spin_tag.h" #include "source_io/module_parameter/parameter.h" #include "source_hamilt/module_hcontainer/hcontainer_funcs.h" @@ -84,7 +85,8 @@ void write_dmr(const std::vector*> dmr, const bool append, const int* iat2iwt, const int nat, - const int istep) + const int istep, + std::ofstream& ofs_running) { const int nspin = dmr.size(); assert(nspin > 0); @@ -110,6 +112,8 @@ void write_dmr(const std::vector*> dmr, const int out_type = 1; std::string fname = PARAM.globalv.global_out_dir + dmr_gen_fname(out_type, ispin, append, istep); write_dmr_csr(fname, ucell, precision, &dm_serial, istep, ispin, nspin); + const std::string spin_tag = ModuleIO::make_spin_tag(ispin, nspin); + ofs_running << " Write DM(R)" << spin_tag << " matrix in NAO basis to file: " << fname << std::endl; } } } diff --git a/source/source_io/module_dm/write_dmr.h b/source/source_io/module_dm/write_dmr.h index 037d6efec61..644e2828b52 100644 --- a/source/source_io/module_dm/write_dmr.h +++ b/source/source_io/module_dm/write_dmr.h @@ -5,6 +5,7 @@ #include "source_hamilt/module_hcontainer/hcontainer.h" #include "source_cell/unitcell.h" +#include #include namespace ModuleIO @@ -62,7 +63,8 @@ void write_dmr(const std::vector*> dmr, const bool append, const int* iat2iwt, const int nat, - const int istep); + const int istep, + std::ofstream& ofs_running); } // namespace ModuleIO #endif diff --git a/source/source_io/module_dos/cal_ldos.cpp b/source/source_io/module_dos/cal_ldos.cpp index d43d843ae61..63f0e2c3f5b 100644 --- a/source/source_io/module_dos/cal_ldos.cpp +++ b/source/source_io/module_dos/cal_ldos.cpp @@ -29,7 +29,8 @@ void Cal_ldos::cal_ldos_lcao( const int nspin, const std::string& global_out_dir, const bool two_fermi, - const int out_ldos_precision) + const int out_ldos_precision, + std::ofstream& ofs_running) { for (int ie = 0; ie < stm_bias[2]; ie++) { @@ -85,7 +86,7 @@ void Cal_ldos::cal_ldos_lcao( // write ldos to cube file std::stringstream fn; - fn << global_out_dir << "LDOS_" << en << "eV" + fn << global_out_dir << "ldos_" << en << "ev" << ".cube"; const int precision = out_ldos_precision; @@ -100,7 +101,9 @@ void Cal_ldos::cal_ldos_lcao( precision, 0, two_fermi, - false); + false, + ofs_running, + "local DOS"); // free memory delete[] ldos; @@ -132,21 +135,22 @@ void cal_ldos_pw(const elecstate::ElecStatePW>* pelec, const bool dos_setemin, const double dos_emin_ev, const double dos_sigma, - const std::vector& ldos_line) + const std::vector& ldos_line, + std::ofstream& ofs_running) { if (out_ldos[0] == 1 || out_ldos[0] == 3) { const int out_ldos_precision = out_ldos[1]; ModuleIO::stm_mode_pw(pelec, psi, ctx, pgrid, ucell, stm_bias, nspin, global_out_dir, two_fermi, - out_ldos_precision); + out_ldos_precision, ofs_running); } if (out_ldos[0] == 2 || out_ldos[0] == 3) { ModuleIO::ldos_mode_pw(pelec, psi, ctx, pgrid, ucell, global_out_dir, nbands, dos_edelta_ev, dos_scale, dos_setemax, dos_emax_ev, dos_setemin, dos_emin_ev, - two_fermi, dos_sigma, ldos_line); + two_fermi, dos_sigma, ldos_line, ofs_running); } } @@ -160,7 +164,8 @@ void stm_mode_pw(const elecstate::ElecStatePW>* pelec, const int nspin, const std::string& global_out_dir, const bool two_fermi, - const int out_ldos_precision) + const int out_ldos_precision, + std::ofstream& ofs_running) { for (int ie = 0; ie < stm_bias[2]; ie++) { @@ -221,11 +226,11 @@ void stm_mode_pw(const elecstate::ElecStatePW>* pelec, } std::stringstream fn; - fn << global_out_dir << "LDOS_" << en << "eV" + fn << global_out_dir << "ldos_" << en << "ev" << ".cube"; const int precision = out_ldos_precision; - ModuleIO::write_vdata_palgrid(pgrid, ldos.data(), 0, nspin, 0, fn.str(), 0, &ucell, precision, 0, two_fermi, false); + ModuleIO::write_vdata_palgrid(pgrid, ldos.data(), 0, nspin, 0, fn.str(), 0, &ucell, precision, 0, two_fermi, false, ofs_running, "local DOS"); } } @@ -245,12 +250,13 @@ void ldos_mode_pw(const elecstate::ElecStatePW>* pelec, const double dos_emin_ev, const bool two_fermi, const double dos_sigma, - const std::vector& ldos_line) + const std::vector& ldos_line, + std::ofstream& ofs_running) { double emax = 0.0; double emin = 0.0; - prepare_dos(GlobalV::ofs_running, + prepare_dos(ofs_running, pelec->eferm, pelec->ekb, pelec->klist->get_nks(), @@ -490,7 +496,8 @@ template void cal_ldos_pw(const elecstate::ElecStatePW< const bool dos_setemin, const double dos_emin_ev, const double dos_sigma, - const std::vector& ldos_line); + const std::vector& ldos_line, + std::ofstream& ofs_running); #if defined(__CUDA) || defined(__ROCM) template void cal_ldos_pw(const elecstate::ElecStatePW>* pelec, const psi::Psi, base_device::DEVICE_GPU>& psi, @@ -510,7 +517,8 @@ template void cal_ldos_pw(const elecstate::ElecStatePW< const bool dos_setemin, const double dos_emin_ev, const double dos_sigma, - const std::vector& ldos_line); + const std::vector& ldos_line, + std::ofstream& ofs_running); #endif } // namespace ModuleIO diff --git a/source/source_io/module_dos/cal_ldos.h b/source/source_io/module_dos/cal_ldos.h index 8888f2d7008..6b6380ee711 100644 --- a/source/source_io/module_dos/cal_ldos.h +++ b/source/source_io/module_dos/cal_ldos.h @@ -11,6 +11,8 @@ #include "source_cell/module_neighbor/sltk_grid_driver.h" // Grid_Driver #include "source_base/matrix.h" // matrix +#include + namespace ModuleIO { template @@ -36,7 +38,8 @@ class Cal_ldos const int nspin, const std::string& global_out_dir, const bool two_fermi, - const int out_ldos_precision); + const int out_ldos_precision, + std::ofstream& ofs_running); }; // namespace Cal_ldos @@ -59,7 +62,8 @@ void cal_ldos_pw(const elecstate::ElecStatePW>* pelec, const bool dos_setemin, const double dos_emin_ev, const double dos_sigma, - const std::vector& ldos_line); + const std::vector& ldos_line, + std::ofstream& ofs_running); template void stm_mode_pw(const elecstate::ElecStatePW>* pelec, @@ -71,7 +75,8 @@ void stm_mode_pw(const elecstate::ElecStatePW>* pelec, const int nspin, const std::string& global_out_dir, const bool two_fermi, - const int out_ldos_precision); + const int out_ldos_precision, + std::ofstream& ofs_running); template void ldos_mode_pw(const elecstate::ElecStatePW>* pelec, @@ -89,7 +94,8 @@ void ldos_mode_pw(const elecstate::ElecStatePW>* pelec, const double dos_emin_ev, const bool two_fermi, const double dos_sigma, - const std::vector& ldos_line); + const std::vector& ldos_line, + std::ofstream& ofs_running); /* * @brief Get grid points and shifts for interpolation. diff --git a/source/source_io/module_elf/write_elf.cpp b/source/source_io/module_elf/write_elf.cpp index ea6e9b88050..fe6164c0107 100644 --- a/source/source_io/module_elf/write_elf.cpp +++ b/source/source_io/module_elf/write_elf.cpp @@ -17,7 +17,8 @@ void write_elf( const UnitCell* ucell_, const int& precision, const std::string& geom_block, - const bool two_fermi) + const bool two_fermi, + std::ofstream& ofs_running) { ModuleBase::timer::start("ModuleIO", "write_elf"); // For nspin = 4, we only calculate the total ELF using the @@ -173,7 +174,9 @@ void write_elf( precision, out_fermi, two_fermi, - false); + false, + ofs_running, + "electron localization function"); } else if (nspin == 2) { @@ -182,11 +185,12 @@ void write_elf( std::string fn_temp = out_dir + "elf" + "s" + std::to_string(is + 1) + geom_block + ".cube"; - const int ispin = is + 1; + // nspin==4 never reaches this branch (only elftot is written), the magnetic label is a placeholder. + const std::string desc_elf = ModuleIO::make_data_desc("electron localization function", "electron localization function magnetization m", is, nspin); ModuleIO::write_vdata_palgrid(pgrid, elf[is].data(), - ispin, + is, nspin, istep_in, fn_temp, @@ -195,7 +199,9 @@ void write_elf( precision, out_fermi, two_fermi, - false); + false, + ofs_running, + desc_elf); } std::vector elf_tot(nrxx, 0.0); @@ -231,7 +237,9 @@ void write_elf( precision, out_fermi, two_fermi, - false); + false, + ofs_running, + "electron localization function"); } ModuleBase::timer::end("ModuleIO", "write_elf"); } // end write_elf diff --git a/source/source_io/module_elf/write_elf.h b/source/source_io/module_elf/write_elf.h index ea5f1b788e7..970843e6526 100644 --- a/source/source_io/module_elf/write_elf.h +++ b/source/source_io/module_elf/write_elf.h @@ -5,6 +5,8 @@ #include "source_basis/module_pw/pw_basis.h" #include "source_estate/module_charge/charge.h" +#include + namespace ModuleIO { void write_elf( @@ -18,7 +20,8 @@ void write_elf( const UnitCell* ucell_, const int& precision, const std::string& geom_block, - const bool two_fermi); + const bool two_fermi, + std::ofstream& ofs_running); } #endif diff --git a/source/source_io/module_energy/write_bands.cpp b/source/source_io/module_energy/write_bands.cpp index 6727b15a7f7..f796b19f07a 100644 --- a/source/source_io/module_energy/write_bands.cpp +++ b/source/source_io/module_energy/write_bands.cpp @@ -7,13 +7,15 @@ #include "source_base/module_parallel/para_bridge.h" #include "source_base/timer.h" #include "source_base/tool_title.h" +#include "source_io/module_output/spin_tag.h" #include "source_io/module_parameter/parameter.h" // mohan add 20250911 #ifdef __MPI #include #endif -void ModuleIO::write_bands(const Input_para& inp, const ModuleBase::matrix& ekb, const K_Vectors& kv) +void ModuleIO::write_bands(const Input_para& inp, const ModuleBase::matrix& ekb, const K_Vectors& kv, + std::ofstream& ofs_running) { // write band information to band.txt if (inp.out_band[0]) @@ -38,9 +40,13 @@ void ModuleIO::write_bands(const Input_para& inp, const ModuleBase::matrix& ekb, ss << ".txt"; + const std::string eig_file = ss.str(); + const std::string spin_tag = ModuleIO::make_spin_tag(is, nspin0); + ofs_running << " Write eigenvalues" << spin_tag << " to file: " << eig_file << std::endl; + const double eshift = 0.0; nscf_bands(is, - ss.str(), + eig_file, inp.nbands, eshift, inp.out_band[1], // precision @@ -73,8 +79,6 @@ void ModuleIO::nscf_bands(const int& is, GlobalV::ofs_running << "\n"; */ - GlobalV::ofs_running << " Write eigenvalues to file: " << eig_file << std::endl; - // number of k points without spin; // nspin = 1,2, nkstot = nkstot_np * nspin; // nspin = 4, nkstot = nkstot_np diff --git a/source/source_io/module_energy/write_bands.h b/source/source_io/module_energy/write_bands.h index c87e050653f..09c2135aff3 100644 --- a/source/source_io/module_energy/write_bands.h +++ b/source/source_io/module_energy/write_bands.h @@ -5,13 +5,16 @@ #include "source_base/matrix.h" #include "source_cell/klist.h" +#include + namespace ModuleIO { void write_bands(const Input_para& inp, const ModuleBase::matrix& ekb, - const K_Vectors& kv); + const K_Vectors& kv, + std::ofstream& ofs_running); /** * @brief calculate the band structure @@ -25,7 +28,7 @@ void write_bands(const Input_para& inp, * @param kv klist */ void nscf_bands(const int& is, - const std::string &eig_file, + const std::string &eig_file, const int& nband, const double& fermie, const int& precision, diff --git a/source/source_io/module_hs/angmom_op_mat.cpp b/source/source_io/module_hs/angmom_op_mat.cpp index d8e3fd2441a..16f8dc7e44d 100644 --- a/source/source_io/module_hs/angmom_op_mat.cpp +++ b/source/source_io/module_hs/angmom_op_mat.cpp @@ -225,17 +225,7 @@ ModuleIO::Angmom_op::Angmom_op( this->ofs_ = &this->fallback_ofs_; } - *ofs_ << "\n\n\n\n"; - *ofs_ << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - *ofs_ << " | |" << std::endl; - *ofs_ << " | Angular momentum expectation value calculation: |" << std::endl; - *ofs_ << " | This is a post-processing step. The expectation value of operator |" << std::endl; - *ofs_ << " | Lx, Ly, Lz (, in which a and b are ABACUS numerical atomic |" << std::endl; - *ofs_ << " | orbitals) will be calculated. |" << std::endl; - *ofs_ << " | The result will be printed to file with name l{x,y,z}[g{step}]_nao.txt |" << std::endl; - *ofs_ << " | |" << std::endl; - *ofs_ << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; - *ofs_ << "\n\n\n\n"; + *ofs_ << " Calculate angular momentum expectation values Lx, Ly, Lz () in NAO basis." << std::endl; int ntype_ = ucell.ntype; Parallel_Common::bcast_int(ntype_); @@ -427,6 +417,7 @@ void ModuleIO::Angmom_op::calculate( ofout << title; this->kernel(&ofout, ucell, d, precision); ofout.close(); + *ofs_ << " Write L(R) (" << d << " component) matrix in NAO basis to file: " << fn << std::endl; } ModuleBase::timer::end("Angmom_op", "calculate"); } diff --git a/source/source_io/module_hs/dhs_sparse_writer.cpp b/source/source_io/module_hs/dhs_sparse_writer.cpp index 9eca747934d..5c4a049546f 100644 --- a/source/source_io/module_hs/dhs_sparse_writer.cpp +++ b/source/source_io/module_hs/dhs_sparse_writer.cpp @@ -6,6 +6,7 @@ #include "source_base/global_function.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" +#include "source_io/module_output/spin_tag.h" #include #include @@ -25,7 +26,8 @@ void ModuleIO::save_dH_sparse(const int& istep, const std::string& calculation, const bool out_app_flag, const int nspin, - const int nlocal) { + const int nlocal, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "save_dH_sparse"); ModuleBase::timer::start("ModuleIO", "save_dH_sparse"); SparseWriteOptions lat_r_options; @@ -149,8 +151,10 @@ void ModuleIO::save_dH_sparse(const int& istep, } else if (!binary) { - GlobalV::ofs_running << " " << label << " data are in file: " - << comp.fname[ispin].str() << std::endl; + const char* matrix_name = (fileflag == "s") ? "dS/dR" : "dH/dR"; + const std::string spin_tag = ModuleIO::make_spin_tag(ispin, nspin); + ofs_running << " Write " << matrix_name << " (" << comp.axis << " component)" << spin_tag + << " matrix in NAO basis to file: " << comp.fname[ispin].str() << std::endl; } comp.ofs[ispin].open(comp.fname[ispin].str().c_str(), mode); detail::check_output_file_open(comp.ofs[ispin], comp.fname[ispin].str(), "ModuleIO::save_dH_sparse"); diff --git a/source/source_io/module_hs/dhs_sparse_writer.h b/source/source_io/module_hs/dhs_sparse_writer.h index 6b84a05c884..8aab9c8007e 100644 --- a/source/source_io/module_hs/dhs_sparse_writer.h +++ b/source/source_io/module_hs/dhs_sparse_writer.h @@ -3,6 +3,7 @@ #include "source_lcao/lcao_hs_arrays.h" +#include #include class Parallel_Orbitals; @@ -21,7 +22,8 @@ void save_dH_sparse(const int& istep, const std::string& calculation, const bool out_app_flag, const int nspin, - const int nlocal); + const int nlocal, + std::ofstream& ofs_running); } // namespace ModuleIO #endif diff --git a/source/source_io/module_hs/hsr_legacy.cpp b/source/source_io/module_hs/hsr_legacy.cpp index 40e37cbdfe1..0f1e37daa3e 100644 --- a/source/source_io/module_hs/hsr_legacy.cpp +++ b/source/source_io/module_hs/hsr_legacy.cpp @@ -28,7 +28,8 @@ void ModuleIO::output_dSR(const int& istep, const MatROutputOptions& options, const bool gamma_only_local, const int npol, - const int nlocal) + const int nlocal, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "output_dSR"); ModuleBase::timer::start("ModuleIO", "output_dSR"); @@ -41,7 +42,7 @@ void ModuleIO::output_dSR(const int& istep, ModuleIO::save_dH_sparse(istep, pv, HS_Arrays, options.sparse_threshold, options.binary, fileflag_s, options.precision, options.global_out_dir, options.global_matrix_dir, options.calculation, options.out_app_flag, - options.nspin, nlocal); + options.nspin, nlocal, ofs_running); sparse_format::destroy_dH_R_sparse(HS_Arrays, options.nspin); @@ -60,7 +61,8 @@ void ModuleIO::output_dHR(const int& istep, const MatROutputOptions& options, const bool gamma_only_local, const int npol, - const int nlocal) + const int nlocal, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "output_dHR"); ModuleBase::timer::start("ModuleIO", "output_dHR"); @@ -86,7 +88,7 @@ void ModuleIO::output_dHR(const int& istep, ModuleIO::save_dH_sparse(istep, pv, HS_Arrays, options.sparse_threshold, options.binary, fileflag_h, options.precision, options.global_out_dir, options.global_matrix_dir, options.calculation, options.out_app_flag, - options.nspin, nlocal); + options.nspin, nlocal, ofs_running); sparse_format::destroy_dH_R_sparse(HS_Arrays, options.nspin); @@ -106,12 +108,13 @@ void ModuleIO::output_SR(Parallel_Orbitals& pv, const std::string& global_matrix_dir, const std::string& calculation, const bool out_app_flag, - const int nspin) + const int nspin, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "output_SR"); ModuleBase::timer::start("ModuleIO", "output_SR"); - GlobalV::ofs_running << " Overlap matrix file is in " << SR_filename << std::endl; + ofs_running << " Write S(R) matrix in NAO basis to file: " << SR_filename << std::endl; LCAO_HS_Arrays HS_Arrays; @@ -164,7 +167,8 @@ void ModuleIO::output_TR(const int istep, const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, const std::string& TR_filename, - const MatROutputOptions& options) + const MatROutputOptions& options, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "output_TR"); ModuleBase::timer::start("ModuleIO", "output_TR"); @@ -182,8 +186,7 @@ void ModuleIO::output_TR(const int istep, { sst << options.global_out_dir << TR_filename; } - GlobalV::ofs_running << " T(R) data are in file: " << sst.str() << std::endl; - + ofs_running << " Write T(R) matrix in NAO basis to file: " << sst.str() << std::endl; sparse_format::cal_TR(ucell, pv, HS_Arrays, grid, two_center_bundle, orb, options.sparse_threshold); ModuleIO::SparseWriteOptions sparse_options; sparse_options.filename = sst.str(); @@ -218,7 +221,8 @@ template void ModuleIO::output_SR(Parallel_Orbitals& pv, const std::string& global_matrix_dir, const std::string& calculation, const bool out_app_flag, - const int nspin); + const int nspin, + std::ofstream& ofs_running); template void ModuleIO::output_SR>(Parallel_Orbitals& pv, const Grid_Driver& grid, hamilt::Hamilt>* p_ham, @@ -230,4 +234,5 @@ template void ModuleIO::output_SR>(Parallel_Orbitals& pv, const std::string& global_matrix_dir, const std::string& calculation, const bool out_app_flag, - const int nspin); + const int nspin, + std::ofstream& ofs_running); diff --git a/source/source_io/module_hs/hsr_legacy.h b/source/source_io/module_hs/hsr_legacy.h index 9c4c61f8b38..334374125e1 100644 --- a/source/source_io/module_hs/hsr_legacy.h +++ b/source/source_io/module_hs/hsr_legacy.h @@ -9,6 +9,7 @@ #include "source_hamilt/hamilt.h" #include "source_lcao/lcao_hs_arrays.h" +#include #include // Legacy LCAO_HS_Arrays-based sparse matrix output path (dH/dR, dS/dR, T(R), S(R)). @@ -43,7 +44,8 @@ void output_dHR(const int& istep, const MatROutputOptions& options, const bool gamma_only_local, const int npol, - const int nlocal); + const int nlocal, + std::ofstream& ofs_running); void output_dSR(const int& istep, const UnitCell& ucell, @@ -55,7 +57,8 @@ void output_dSR(const int& istep, const MatROutputOptions& options, const bool gamma_only_local, const int npol, - const int nlocal); + const int nlocal, + std::ofstream& ofs_running); void output_TR(const int istep, const UnitCell& ucell, @@ -65,7 +68,8 @@ void output_TR(const int istep, const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, const std::string& TR_filename, - const MatROutputOptions& options); + const MatROutputOptions& options, + std::ofstream& ofs_running); template void output_SR(Parallel_Orbitals& pv, @@ -79,7 +83,8 @@ void output_SR(Parallel_Orbitals& pv, const std::string& global_matrix_dir, const std::string& calculation, const bool out_app_flag, - const int nspin); + const int nspin, + std::ofstream& ofs_running); } // namespace ModuleIO diff --git a/source/source_io/module_hs/hsr_writer.cpp b/source/source_io/module_hs/hsr_writer.cpp index 7a83b2c57cb..10b5218fa69 100644 --- a/source/source_io/module_hs/hsr_writer.cpp +++ b/source/source_io/module_hs/hsr_writer.cpp @@ -1,4 +1,5 @@ #include "hsr_writer.h" +#include "source_io/module_output/spin_tag.h" #include "source_base/module_out/sparse_matrix.h" #include "source_base/tool_quit.h" @@ -6,6 +7,7 @@ #include "source_cell/unitcell.h" #include "source_base/global_variable.h" #include "source_basis/module_ao/parallel_orbitals.h" +#include "source_estate/fp_energy.h" #include "source_hamilt/module_hcontainer/hcontainer.h" #include "source_hamilt/module_hcontainer/hcontainer_funcs.h" #include "source_hamilt/module_hcontainer/output_hcontainer.h" @@ -81,7 +83,9 @@ void ModuleIO::write_hcontainer_csr(const std::string& fname, const int ispin, const int nspin, const std::string& label, - const std::string& representation_note) + const std::string& representation_note, + const double efermi_eV, + const bool has_efermi) { std::ofstream ofs; if (istep <= 0) @@ -101,7 +105,12 @@ void ModuleIO::write_hcontainer_csr(const std::string& fname, ofs << " --- Ionic Step " << istep + 1 << " ---" << std::endl; ofs << " # print " << label << " matrix in real space " << label << "(R)" << std::endl; ofs << " " << nspin << " # number of spin directions" << std::endl; - ofs << " " << ispin + 1 << " # spin index" << std::endl; + ofs << " " << ispin + 1 << " # spin index"; + if (has_efermi) + { + ofs << ", E_Fermi = " << std::setprecision(6) << efermi_eV << " eV"; + } + ofs << std::endl; ofs << " " << mat_serial->get_nbasis() << " # number of localized basis" << std::endl; ofs << " " << mat_serial->size_R_loop() << " # number of Bravais lattice vector R" << std::endl; ofs << std::endl; @@ -279,7 +288,9 @@ void ModuleIO::write_hsr(const std::vector*>& hr_vec, const int* iat2iwt, const int nat, const int istep, - const std::string& global_out_dir) + const std::string& global_out_dir, + const elecstate::Efermi& eferm, + std::ofstream& ofs_running) { if (out_type != 1 && out_type != 2) { @@ -317,9 +328,12 @@ void ModuleIO::write_hsr(const std::vector*>& hr_vec, } else { + const double efermi_eV = eferm.get_efval(ispin) * ModuleBase::Ry_to_eV; write_hcontainer_csr( - fname, ucell, precision, &hr_serial, istep, ispin, nspin, "H", representation_note); + fname, ucell, precision, &hr_serial, istep, ispin, nspin, "H", representation_note, efermi_eV, true); } + const std::string spin_tag = ModuleIO::make_spin_tag(ispin, nspin); + ofs_running << " Write H(R)" << spin_tag << " matrix in NAO basis to file: " << fname << std::endl; } } @@ -348,8 +362,9 @@ void ModuleIO::write_hsr(const std::vector*>& hr_vec, else { write_hcontainer_csr( - fname, ucell, precision, &sr_serial, istep, 0, 1, "S", representation_note); + fname, ucell, precision, &sr_serial, istep, 0, 1, "S", representation_note, 0.0, false); } + ofs_running << " Write S(R) matrix in NAO basis to file: " << fname << std::endl; } } } @@ -357,10 +372,10 @@ void ModuleIO::write_hsr(const std::vector*>& hr_vec, // Explicit instantiations template void ModuleIO::write_hcontainer_csr( const std::string&, const UnitCell*, const int, - hamilt::HContainer*, const int, const int, const int, const std::string&, const std::string&); + hamilt::HContainer*, const int, const int, const int, const std::string&, const std::string&, const double, const bool); template void ModuleIO::write_hcontainer_csr>( const std::string&, const UnitCell*, const int, - hamilt::HContainer>*, const int, const int, const int, const std::string&, const std::string&); + hamilt::HContainer>*, const int, const int, const int, const std::string&, const std::string&, const double, const bool); template void ModuleIO::write_hcontainer_csr_binary( const std::string&, hamilt::HContainer*, const int, const bool); @@ -372,10 +387,10 @@ template void ModuleIO::write_hsr( const hamilt::HContainer*, const UnitCell*, const int, const int, const Parallel_2D&, const bool, const bool, const int*, const int, const int, - const std::string&); + const std::string&, const elecstate::Efermi&, std::ofstream&); template void ModuleIO::write_hsr>( const std::vector>*>&, const hamilt::HContainer>*, const UnitCell*, const int, const int, const Parallel_2D&, const bool, const bool, const int*, const int, const int, - const std::string&); + const std::string&, const elecstate::Efermi&, std::ofstream&); diff --git a/source/source_io/module_hs/hsr_writer.h b/source/source_io/module_hs/hsr_writer.h index b4f836035da..3fb8e6bc5a3 100644 --- a/source/source_io/module_hs/hsr_writer.h +++ b/source/source_io/module_hs/hsr_writer.h @@ -9,6 +9,11 @@ class UnitCell; +namespace elecstate +{ +struct Efermi; +} + namespace hamilt { template @@ -44,6 +49,8 @@ std::string dhr_gen_fname(const std::string& prefix, const int istep); /// Write a single HContainer to CSR file with header. +/// @param efermi_eV the Fermi energy in eV for this spin channel +/// @param has_efermi whether to append the Fermi energy to the spin-index line template void write_hcontainer_csr(const std::string& fname, const UnitCell* ucell, @@ -53,7 +60,9 @@ void write_hcontainer_csr(const std::string& fname, const int ispin, const int nspin, const std::string& label, - const std::string& representation_note); + const std::string& representation_note, + const double efermi_eV, + const bool has_efermi); /// Write one HContainer record in the native binary CSR format. template @@ -75,7 +84,9 @@ void write_hsr(const std::vector*>& hr_vec, const int* iat2iwt, const int nat, const int istep, - const std::string& global_out_dir); + const std::string& global_out_dir, + const elecstate::Efermi& eferm, + std::ofstream& ofs_running); } // namespace ModuleIO diff --git a/source/source_io/module_hs/hterm_writer.cpp b/source/source_io/module_hs/hterm_writer.cpp index 30c1921b923..a4be8399b56 100644 --- a/source/source_io/module_hs/hterm_writer.cpp +++ b/source/source_io/module_hs/hterm_writer.cpp @@ -1,4 +1,5 @@ #include "hterm_writer.h" +#include "source_io/module_output/spin_tag.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" @@ -67,6 +68,7 @@ static void setup_veff_hcontainer(hamilt::HContainer& hR, static void gather_and_write(const std::string& prefix, const std::string& label, + const std::string& term_name, hamilt::HContainer& hR, const UnitCell& ucell, const Parallel_Orbitals& pv, @@ -79,7 +81,8 @@ static void gather_and_write(const std::string& prefix, const std::string& calculation, const bool out_app_flag, const std::string& global_out_dir, - const std::string& global_matrix_dir) + const std::string& global_matrix_dir, + std::ofstream* ofs_running) { const int nbasis = hR.get_nbasis(); #ifdef __MPI @@ -94,11 +97,20 @@ static void gather_and_write(const std::string& prefix, const bool md_no_append = (calculation == "md") && !out_app_flag; const std::string& out_dir = md_no_append ? global_matrix_dir : global_out_dir; const std::string fname = out_dir + hsr_gen_fname(prefix, ispin, append, istep); + // H(R) component terms (T, V^NL, V^L, V^H, V^XC, V^EXX) are parts of + // the Hamiltonian, not the full H(R): they do not carry the Fermi annotation + const double no_efermi = 0.0; #ifdef __MPI - write_hcontainer_csr(fname, &ucell, 8, &hr_serial, istep, ispin, nspin, label, ""); + write_hcontainer_csr(fname, &ucell, 8, &hr_serial, istep, ispin, nspin, label, "", no_efermi, false); #else - write_hcontainer_csr(fname, &ucell, 8, &hR, istep, ispin, nspin, label, ""); + write_hcontainer_csr(fname, &ucell, 8, &hR, istep, ispin, nspin, label, "", no_efermi, false); #endif + if (ofs_running != nullptr) + { + const std::string spin_tag = ModuleIO::make_spin_tag(ispin, nspin); + *ofs_running << " Write H(R) (" << term_name << " term)" << spin_tag << " matrix in NAO basis to file: " << fname + << std::endl; + } } } @@ -163,6 +175,7 @@ struct HTermSpec std::string k_prefix; std::string r_prefix; std::string r_label; + std::string term_name; // Fills hR_tmp for spin channel ispin. setup_veff_hcontainer has already // been run when use_veff_layout is true; the builder adds the term values. std::function&, int)> build; @@ -201,9 +214,9 @@ static void write_h_term(WriteHParams& params, const HTermSpec& spec) if (also_hR) { - gather_and_write(spec.r_prefix, spec.r_label, hR_tmp, ucell, pv, nspin, ispin, istep, append, - iat2iwt, nat, params.calculation, params.out_app_flag, - params.global_out_dir, params.global_matrix_dir); + gather_and_write(spec.r_prefix, spec.r_label, spec.term_name, hR_tmp, ucell, pv, nspin, ispin, + istep, append, iat2iwt, nat, params.calculation, params.out_app_flag, + params.global_out_dir, params.global_matrix_dir, params.ofs_running); } } @@ -225,6 +238,7 @@ void write_h_t(WriteHParams& params) spec.k_prefix = "tk"; spec.r_prefix = "t"; spec.r_label = "T"; + spec.term_name = "kinetic"; spec.build = [&](hamilt::HContainer& hR_tmp, int /*ispin*/) { hamilt::EKinetic> tmp_ekinetic(nullptr, kv.kvec_d, &hR_tmp, &ucell, orb_cutoff, &gd, @@ -248,6 +262,7 @@ void write_h_vnl(WriteHParams& params) spec.k_prefix = "vnlk"; spec.r_prefix = "vnl"; spec.r_label = "V^NL"; + spec.term_name = "nonlocal"; spec.build = [&](hamilt::HContainer& hR_tmp, int /*ispin*/) { hamilt::Nonlocal> tmp_nonlocal(nullptr, kv.kvec_d, &hR_tmp, &ucell, orb_cutoff, &gd, @@ -273,6 +288,7 @@ void write_h_vl(WriteHParams& params) spec.k_prefix = "vlk"; spec.r_prefix = "vl"; spec.r_label = "V^L"; + spec.term_name = "local"; spec.build = [&](hamilt::HContainer& hR_tmp, int /*ispin*/) { setup_veff_hcontainer(hR_tmp, ucell, gd, pv, orb_cutoff); ModuleGint::cal_gint_vl(v_local, &hR_tmp); @@ -298,6 +314,7 @@ void write_h_vh(WriteHParams& params) spec.k_prefix = "vhk"; spec.r_prefix = "vh"; spec.r_label = "V^H"; + spec.term_name = "hartree"; spec.build = [&](hamilt::HContainer& hR_tmp, int ispin) { setup_veff_hcontainer(hR_tmp, ucell, gd, pv, orb_cutoff); ModuleGint::cal_gint_vl(&v_h(ispin, 0), &hR_tmp); @@ -334,6 +351,7 @@ void write_h_vxc(WriteHParams& params) spec.k_prefix = "vxck"; spec.r_prefix = "vxc"; spec.r_label = "V^XC"; + spec.term_name = "xc"; spec.build = [&](hamilt::HContainer& hR_tmp, int ispin) { setup_veff_hcontainer(hR_tmp, ucell, gd, pv, orb_cutoff); ModuleGint::cal_gint_vl(&v_xc(ispin, 0), &hR_tmp); @@ -363,7 +381,8 @@ static void write_h_exx_impl(const UnitCell& ucell, const bool out_app_flag, const std::string& ks_solver, const int drank, - const Exx_Info& exx_info) + const Exx_Info& exx_info, + std::ofstream* ofs_running) { const auto& Hexxs = ex->get_Hexxs(); // vector over spin of map> const int nspin_out = (nspin == 2 ? 2 : 1); @@ -382,8 +401,9 @@ static void write_h_exx_impl(const UnitCell& ucell, if (also_hR) { - gather_and_write("vexx", "V^EXX", hR_tmp, ucell, pv, nspin, ispin, istep, append, iat2iwt, nat, - calculation, out_app_flag, global_out_dir, global_matrix_dir); + gather_and_write("vexx", "V^EXX", "exx", hR_tmp, ucell, pv, nspin, ispin, istep, append, + iat2iwt, nat, calculation, out_app_flag, global_out_dir, global_matrix_dir, + ofs_running); } } } @@ -413,7 +433,8 @@ void write_h_exx(WriteHParams& params, const Exx_Info& exx_info) write_h_exx_impl(ucell, pv, params.exd, kv, nspin, istep, append, iat2iwt, nat, also_hR, params.npol, params.nlocal, params.gamma_only_local, params.global_out_dir, params.global_matrix_dir, - params.calculation, params.out_app_flag, params.ks_solver, params.drank, exx_info); + params.calculation, params.out_app_flag, params.ks_solver, params.drank, exx_info, + params.ofs_running); } } else @@ -423,7 +444,8 @@ void write_h_exx(WriteHParams& params, const Exx_Info& exx_info) write_h_exx_impl(ucell, pv, params.exc, kv, nspin, istep, append, iat2iwt, nat, also_hR, params.npol, params.nlocal, params.gamma_only_local, params.global_out_dir, params.global_matrix_dir, - params.calculation, params.out_app_flag, params.ks_solver, params.drank, exx_info); + params.calculation, params.out_app_flag, params.ks_solver, params.drank, exx_info, + params.ofs_running); } } diff --git a/source/source_io/module_hs/hterm_writer.h b/source/source_io/module_hs/hterm_writer.h index 64e5572d55f..7aa2f721fc6 100644 --- a/source/source_io/module_hs/hterm_writer.h +++ b/source/source_io/module_hs/hterm_writer.h @@ -12,6 +12,7 @@ #include "source_hamilt/module_xc/exx_info.h" #include +#include #include #include @@ -51,6 +52,7 @@ struct WriteHParams std::string global_matrix_dir; std::string ks_solver; int drank = 0; + std::ofstream* ofs_running = nullptr; // running log for the "Write ... to file" messages #ifdef __EXX // The gamma-only (TK==double) exx interfaces used by the EXX term. // Deliberately NOT templated on TK, because it would force WriteHParams and diff --git a/source/source_io/module_hs/output_mat_sparse.cpp b/source/source_io/module_hs/output_mat_sparse.cpp index d8d737d8956..5179bd5477e 100644 --- a/source/source_io/module_hs/output_mat_sparse.cpp +++ b/source/source_io/module_hs/output_mat_sparse.cpp @@ -17,7 +17,8 @@ void output_mat_sparse(const MatSparseOutputOptions& options, const Grid_Driver& grid, const K_Vectors& kv, hamilt::Hamilt* p_ham, - Plus_U_Base* p_dftu) + Plus_U_Base* p_dftu, + std::ofstream& ofs_running) { LCAO_HS_Arrays HS_Arrays; // store sparse arrays @@ -52,7 +53,8 @@ void output_mat_sparse(const MatSparseOutputOptions& options, two_center_bundle, orb, tr_filename, - mat_R_options); + mat_R_options, + ofs_running); } //! generate a file containing the derivatives of the Hamiltonian matrix (in Ry/Bohr) @@ -70,7 +72,8 @@ void output_mat_sparse(const MatSparseOutputOptions& options, mat_R_options, gamma_only_local, npol, - nlocal); + nlocal, + ofs_running); } //! generate a file containing the derivatives of the overlap matrix (in Ry/Bohr) if (options.out_mat_ds) @@ -86,7 +89,8 @@ void output_mat_sparse(const MatSparseOutputOptions& options, mat_R_options, gamma_only_local, npol, - nlocal); + nlocal, + ofs_running); } // add by jingan for out r_R matrix 2019.8.14 @@ -106,7 +110,8 @@ void output_mat_sparse(const MatSparseOutputOptions& options, calculation, out_app_flag, nlocal, - npol); + npol, + ofs_running); } return; @@ -122,7 +127,8 @@ template void output_mat_sparse(const MatSparseOutputOptions& options, const Grid_Driver& grid, const K_Vectors& kv, hamilt::Hamilt* p_ham, - Plus_U_Base* p_dftu); + Plus_U_Base* p_dftu, + std::ofstream& ofs_running); template void output_mat_sparse>(const MatSparseOutputOptions& options, const int& istep, @@ -134,6 +140,7 @@ template void output_mat_sparse>(const MatSparseOutputOptio const Grid_Driver& grid, const K_Vectors& kv, hamilt::Hamilt>* p_ham, - Plus_U_Base* p_dftu); + Plus_U_Base* p_dftu, + std::ofstream& ofs_running); } // namespace ModuleIO diff --git a/source/source_io/module_hs/output_mat_sparse.h b/source/source_io/module_hs/output_mat_sparse.h index 2d0b698c8bb..86e780758cc 100644 --- a/source/source_io/module_hs/output_mat_sparse.h +++ b/source/source_io/module_hs/output_mat_sparse.h @@ -8,6 +8,8 @@ #include "source_cell/module_neighbor/sltk_grid_driver.h" #include "source_pw/module_pwdft/dftu_base.h" // mohan add 20251107 +#include + namespace ModuleIO { struct MatSparseOutputOptions @@ -36,7 +38,8 @@ void output_mat_sparse(const MatSparseOutputOptions& options, const Grid_Driver& grid, const K_Vectors& kv, hamilt::Hamilt* p_ham, - Plus_U_Base* p_dftu); + Plus_U_Base* p_dftu, + std::ofstream& ofs_running); } // namespace ModuleIO #endif // OUTPUT_MAT_SPARSE_H diff --git a/source/source_io/module_hs/pos_op_mat.cpp b/source/source_io/module_hs/pos_op_mat.cpp index b02e30eb79c..04ec8d45e85 100644 --- a/source/source_io/module_hs/pos_op_mat.cpp +++ b/source/source_io/module_hs/pos_op_mat.cpp @@ -91,7 +91,8 @@ void Position_op::out_rR(const UnitCell& ucell, const std::string& calculation, const bool out_app_flag, const int nlocal, - const int npol) + const int npol, + std::ofstream& ofs_running) { ModuleBase::TITLE("Position_op", "out_rR"); ModuleBase::timer::start("Position_op", "out_rR"); @@ -107,7 +108,8 @@ void Position_op::out_rR(const UnitCell& ucell, nlocal, npol, sparse_threshold, - binary); + binary, + ofs_running); ModuleBase::timer::end("Position_op", "out_rR"); } diff --git a/source/source_io/module_hs/pos_op_mat.h b/source/source_io/module_hs/pos_op_mat.h index 63def4d9bec..96a6513c23e 100644 --- a/source/source_io/module_hs/pos_op_mat.h +++ b/source/source_io/module_hs/pos_op_mat.h @@ -78,7 +78,8 @@ class Position_op const std::string& calculation, const bool out_app_flag, const int nlocal, - const int npol); + const int npol, + std::ofstream& ofs_running); private: std::unique_ptr basis_; diff --git a/source/source_io/module_hs/pos_op_writer.cpp b/source/source_io/module_hs/pos_op_writer.cpp index be1d1dc2a42..15e1fa201ea 100644 --- a/source/source_io/module_hs/pos_op_writer.cpp +++ b/source/source_io/module_hs/pos_op_writer.cpp @@ -26,7 +26,8 @@ void PosOpWriter::out_lat_r(const UnitCell& ucell, int nlocal, int npol, double sparse_threshold, - bool binary) + bool binary, + std::ofstream& ofs_running) { ModuleBase::TITLE("PosOpWriter", "out_lat_r"); ModuleBase::timer::start("PosOpWriter", "out_lat_r"); @@ -236,6 +237,7 @@ void PosOpWriter::out_lat_r(const UnitCell& ucell, binary, open_in_append, "PosOpWriter::out_lat_r"); + ofs_running << " Write r(R) matrix in NAO basis to file: " << ssr.str() << std::endl; std::remove(tem1.str().c_str()); } diff --git a/source/source_io/module_hs/pos_op_writer.h b/source/source_io/module_hs/pos_op_writer.h index 52c74e6de7a..1fd22bcfc81 100644 --- a/source/source_io/module_hs/pos_op_writer.h +++ b/source/source_io/module_hs/pos_op_writer.h @@ -39,7 +39,8 @@ class PosOpWriter int nlocal, int npol, double sparse_threshold, - bool binary); + bool binary, + std::ofstream& ofs_running); private: const PosOpBasis& basis_; diff --git a/source/source_io/module_hs/unittests/CMakeLists.txt b/source/source_io/module_hs/unittests/CMakeLists.txt index af2c1ff79f9..05bc0102d85 100644 --- a/source/source_io/module_hs/unittests/CMakeLists.txt +++ b/source/source_io/module_hs/unittests/CMakeLists.txt @@ -70,6 +70,7 @@ AddTest( SOURCES test_hsr_writer.cpp ../hsr_writer.cpp + ../../../source_estate/fp_energy.cpp ../../../source_cell/ucell_io.cpp ../../../source_basis/module_ao/parallel_orbitals.cpp ../../../source_hamilt/module_hcontainer/test/tmp_mocks.cpp diff --git a/source/source_io/module_hs/unittests/test_dhs_sparse_writer.cpp b/source/source_io/module_hs/unittests/test_dhs_sparse_writer.cpp index e2335454cc7..3885c7c9472 100644 --- a/source/source_io/module_hs/unittests/test_dhs_sparse_writer.cpp +++ b/source/source_io/module_hs/unittests/test_dhs_sparse_writer.cpp @@ -16,6 +16,11 @@ #include "source_io/module_hs/dhs_sparse_writer.h" +#include + +// Dummy output stream for running-log writes in tests; discards output. +static std::ofstream ofs_running; + #ifdef __MPI #include #endif @@ -35,7 +40,7 @@ TEST(DhsSparseWriter, TextCountsOnlyValuesAboveThreshold) arrays.dHRx_sparse[0][r_vector][1][0] = 0.0; arrays.dHRx_sparse[0][r_vector][1][1] = -2.0; - ModuleIO::save_dH_sparse(5, pv, arrays, 1e-10, false, "h", 8, "./", "./", "scf", false, 1, 2); + ModuleIO::save_dH_sparse(5, pv, arrays, 1e-10, false, "h", 8, "./", "./", "scf", false, 1, 2, ofs_running); const std::vector lines = read_lines("dhrxs1g6_nao.csr"); ASSERT_GE(lines.size(), 7); @@ -86,7 +91,7 @@ TEST(DhsSparseWriter, BinaryCountsOnlyValuesAboveThreshold) arrays.dHRx_sparse[0][r_vector][1][0] = 0.0; arrays.dHRx_sparse[0][r_vector][1][1] = -2.0; - ModuleIO::save_dH_sparse(6, pv, arrays, 1e-10, true, "h", 8, "./", "./", "scf", false, 1, 2); + ModuleIO::save_dH_sparse(6, pv, arrays, 1e-10, true, "h", 8, "./", "./", "scf", false, 1, 2, ofs_running); std::ifstream ifs("dhrxs1g7_nao.csr", std::ios::binary); ASSERT_TRUE(ifs.is_open()); @@ -122,7 +127,7 @@ TEST(DhsSparseWriter, SocWritesAllDirections) arrays.dHRy_soc_sparse[r_vector][0][1] = std::complex(2.0, -1.0); arrays.dHRz_soc_sparse[r_vector][1][1] = std::complex(-3.0, 0.5); - ModuleIO::save_dH_sparse(7, pv, arrays, 1e-10, false, "s", 8, "./", "./", "scf", false, 4, 2); + ModuleIO::save_dH_sparse(7, pv, arrays, 1e-10, false, "s", 8, "./", "./", "scf", false, 4, 2, ofs_running); const std::string x_output = read_file("dsrxs1g8_nao.csr"); const std::string y_output = read_file("dsrys1g8_nao.csr"); diff --git a/source/source_io/module_hs/unittests/test_hsr_writer.cpp b/source/source_io/module_hs/unittests/test_hsr_writer.cpp index f4a9e6c90d8..b2f0079681f 100644 --- a/source/source_io/module_hs/unittests/test_hsr_writer.cpp +++ b/source/source_io/module_hs/unittests/test_hsr_writer.cpp @@ -11,6 +11,8 @@ #include "csr_test_helpers.h" #include "source_base/module_out/csr_reader.h" +#include "source_base/constants.h" +#include "source_estate/fp_energy.h" #include "source_hamilt/module_hcontainer/hcontainer_funcs.h" #include "source_io/module_hs/hsr_writer.h" @@ -18,6 +20,7 @@ #include #include #include +#include #include #include @@ -158,13 +161,15 @@ TEST(HsrWriterIo, HContainerCsrHeaderKeepsCurrentFormat) double values[4] = {1.0, 0.0, 0.5, 2.0}; fill_matrix(matrix, pv, values); - ModuleIO::write_hcontainer_csr(filename, &ucell, 5, &matrix, 0, 0, 1, "H", ""); + // label "H": header carries the Fermi energy of this spin channel + const double efermi_eV = 5.4321; + ModuleIO::write_hcontainer_csr(filename, &ucell, 5, &matrix, 0, 0, 1, "H", "", efermi_eV, true); const std::string output = read_file(filename); EXPECT_THAT(output, testing::HasSubstr(" --- Ionic Step 1 ---\n")); EXPECT_THAT(output, testing::HasSubstr(" # print H matrix in real space H(R)\n")); EXPECT_THAT(output, testing::HasSubstr(" 1 # number of spin directions\n")); - EXPECT_THAT(output, testing::HasSubstr(" 1 # spin index\n")); + EXPECT_THAT(output, testing::HasSubstr(" 1 # spin index, E_Fermi = 5.4321 eV\n")); EXPECT_THAT(output, testing::HasSubstr(" 2 # number of localized basis\n")); EXPECT_THAT(output, testing::HasSubstr(" 1 # number of Bravais lattice vector R\n")); EXPECT_THAT(output, testing::HasSubstr(" user_defined_lattice\n")); @@ -195,13 +200,17 @@ TEST(HsrWriterIo, GammaFoldedHeaderKeepsCsrReadable) double values[4] = {1.0, 0.0, 0.5, 2.0}; fill_matrix(matrix, pv, values); + // label "H" without Fermi energy available: pass has_efermi = false + const double no_efermi = 0.0; ModuleIO::write_hcontainer_csr( - filename, &ucell, 5, &matrix, 0, 0, 1, "H", representation_note); + filename, &ucell, 5, &matrix, 0, 0, 1, "H", representation_note, no_efermi, false); const std::string output = read_file(filename); EXPECT_THAT(output, testing::HasSubstr("# representation: " + representation_note + "\n")); EXPECT_THAT(output, testing::HasSubstr(" 1 # number of Bravais lattice vector R\n")); EXPECT_THAT(output, testing::HasSubstr(" 0 0 0 3\n")); + // has_efermi = false: no Fermi annotation even though label is "H" + EXPECT_THAT(output, testing::Not(testing::HasSubstr("E_Fermi"))); ModuleIO::csrFileReader reader(filename); ASSERT_EQ(reader.getNumberOfR(), 1); @@ -224,18 +233,51 @@ TEST(HsrWriterIo, HContainerCsrAppendKeepsCurrentStepSections) double values[4] = {1.0, 0.0, 0.0, 1.0}; fill_matrix(matrix, pv, values); - ModuleIO::write_hcontainer_csr(filename, &ucell, 4, &matrix, 0, 0, 1, "S", ""); - ModuleIO::write_hcontainer_csr(filename, &ucell, 4, &matrix, 1, 0, 1, "S", ""); + // label "S": header carries no Fermi energy (has_efermi = false) + const double ignored_efermi = 0.0; + ModuleIO::write_hcontainer_csr(filename, &ucell, 4, &matrix, 0, 0, 1, "S", "", ignored_efermi, false); + ModuleIO::write_hcontainer_csr(filename, &ucell, 4, &matrix, 1, 0, 1, "S", "", ignored_efermi, false); const std::string output = read_file(filename); EXPECT_EQ(count_substr(output, " --- Ionic Step "), 2); EXPECT_THAT(output, testing::HasSubstr(" --- Ionic Step 1 ---\n")); EXPECT_THAT(output, testing::HasSubstr(" --- Ionic Step 2 ---\n")); EXPECT_EQ(count_substr(output, " # print S matrix in real space S(R)\n"), 2); + // S(R) files must not contain a Fermi energy annotation + EXPECT_THAT(output, testing::Not(testing::HasSubstr("E_Fermi"))); std::remove(filename.c_str()); } +TEST(HsrWriterIo, HContainerCsrHeaderCarriesPerChannelFermi) +{ + UnitCell ucell; + init_unitcell(ucell); + Parallel_Orbitals pv; + init_serial_orbitals(pv); + hamilt::HContainer matrix(&pv); + double values[4] = {1.0, 0.0, 0.5, 2.0}; + fill_matrix(matrix, pv, values); + + // Each spin channel is written to its own file (write mode, istep=0) + const double efermi_up_eV = 5.4321; + const double efermi_dw_eV = 3.2109; + const std::string filename_up = "write_hs_r_header_fermi_up.csr"; + const std::string filename_dw = "write_hs_r_header_fermi_dw.csr"; + std::remove(filename_up.c_str()); + std::remove(filename_dw.c_str()); + ModuleIO::write_hcontainer_csr(filename_up, &ucell, 5, &matrix, 0, 0, 2, "H", "", efermi_up_eV, true); + ModuleIO::write_hcontainer_csr(filename_dw, &ucell, 5, &matrix, 0, 1, 2, "H", "", efermi_dw_eV, true); + + const std::string output_up = read_file(filename_up); + const std::string output_dw = read_file(filename_dw); + EXPECT_THAT(output_up, testing::HasSubstr(" 1 # spin index, E_Fermi = 5.4321 eV\n")); + EXPECT_THAT(output_dw, testing::HasSubstr(" 2 # spin index, E_Fermi = 3.2109 eV\n")); + + std::remove(filename_up.c_str()); + std::remove(filename_dw.c_str()); +} + // --------------------------------------------------------------------------- // write_hcontainer_csr_binary: native binary layout // --------------------------------------------------------------------------- @@ -407,8 +449,13 @@ TEST(HsrWriterIo, HContainerBinaryMpiGatherWritesCompleteFiles) init_sparse_output_globals(); std::vector*> hr_vec(1, &hr_parallel); + elecstate::Efermi eferm; + eferm.two_efermi = false; + eferm.ef = 0.4; // Ry; used only for the H(R) header + std::ofstream ofs_running_null; // not opened; test does not inspect the running log ModuleIO::write_hsr( - hr_vec, &sr_parallel, &ucell, 2, 8, parallel_pv, true, true, iat2iwt, 1, 0, "./"); + hr_vec, &sr_parallel, &ucell, 2, 8, parallel_pv, true, true, iat2iwt, 1, 0, "./", eferm, + ofs_running_null); MPI_Barrier(MPI_COMM_WORLD); if (mpi_rank == 0) @@ -438,6 +485,82 @@ TEST(HsrWriterIo, HContainerBinaryMpiGatherWritesCompleteFiles) #endif } +// --------------------------------------------------------------------------- +// write_hsr: text CSR (out_type=1) threads per-spin Fermi energy from Efermi +// --------------------------------------------------------------------------- + +TEST(HsrWriterIo, WriteHsrTextCsrCarriesPerSpinFermiFromEfermi) +{ + // Exercise the fix path inside write_hsr: eferm.get_efval(ispin) * Ry_to_eV + // for a two-Fermi (nspin=2) case. The only existing write_hsr test uses + // out_type=2 (binary), which skips the Fermi branch entirely. + // append=true suppresses the ionic-step suffix in the generated names. + const std::string hr_up_filename = "hrs1_nao.csr"; + const std::string hr_dw_filename = "hrs2_nao.csr"; + const std::string sr_filename = "sr_nao.csr"; + std::remove(hr_up_filename.c_str()); + std::remove(hr_dw_filename.c_str()); + std::remove(sr_filename.c_str()); + + UnitCell ucell; + init_unitcell(ucell); + // The __MPI path in write_hsr calls set_atomic_trace and gatherParallels, + // which require a real atom-to-orbital map; nullptr with nat=0 throws on + // that path. iat2it is released by UnitCell's Statistics member. + ucell.iat2it = new int[1]; + ucell.iat2it[0] = 0; + ucell.set_iat2iwt(1); + const int* iat2iwt_ptr = ucell.get_iat2iwt(); + Parallel_Orbitals pv; + init_serial_orbitals(pv); + // gatherParallels reaches get_indexes_row(iat) on this source pv, which + // dereferences iat2iwt_; init_serial_orbitals does not set it, so the + // atomic trace must be installed explicitly. + pv.set_atomic_trace(iat2iwt_ptr, 1, 2); + + hamilt::HContainer hr_up(&pv); + hamilt::HContainer hr_dw(&pv); + hamilt::HContainer sr(&pv); + double values[4] = {1.0, 0.0, 0.5, 2.0}; + double sr_values[4] = {1.0, 0.0, 0.0, 1.0}; + fill_matrix(hr_up, pv, values); + fill_matrix(hr_dw, pv, values); + fill_matrix(sr, pv, sr_values); + + init_sparse_output_globals(); + + elecstate::Efermi eferm; + eferm.two_efermi = true; + eferm.ef_up = 0.5; // Ry + eferm.ef_dw = 0.3; // Ry + const double ef_up_eV = eferm.ef_up * ModuleBase::Ry_to_eV; // 6.802849 + const double ef_dw_eV = eferm.ef_dw * ModuleBase::Ry_to_eV; // 4.0817094 + + std::vector*> hr_vec; + hr_vec.push_back(&hr_up); + hr_vec.push_back(&hr_dw); + + std::ofstream ofs_running_null; // not opened; test does not inspect the running log + ModuleIO::write_hsr( + hr_vec, &sr, &ucell, 1, 8, pv, true, true, iat2iwt_ptr, 1, 0, "./", eferm, + ofs_running_null); + + const std::string output_up = read_file(hr_up_filename); + const std::string output_dw = read_file(hr_dw_filename); + + // write_hsr multiplies each spin's Fermi by Ry_to_eV and passes has_efermi=true + std::ostringstream expected_up; + expected_up << " 1 # spin index, E_Fermi = " << std::setprecision(6) << ef_up_eV << " eV\n"; + std::ostringstream expected_dw; + expected_dw << " 2 # spin index, E_Fermi = " << std::setprecision(6) << ef_dw_eV << " eV\n"; + EXPECT_THAT(output_up, testing::HasSubstr(expected_up.str())); + EXPECT_THAT(output_dw, testing::HasSubstr(expected_dw.str())); + + std::remove(hr_up_filename.c_str()); + std::remove(hr_dw_filename.c_str()); + std::remove(sr_filename.c_str()); +} + int main(int argc, char** argv) { #ifdef __MPI diff --git a/source/source_io/module_hs/vxc_op_r.cpp b/source/source_io/module_hs/vxc_op_r.cpp index 38345a7be3a..c1ca15afc4a 100644 --- a/source/source_io/module_hs/vxc_op_r.cpp +++ b/source/source_io/module_hs/vxc_op_r.cpp @@ -1,4 +1,5 @@ #include "vxc_op_r.h" +#include "source_io/module_output/spin_tag.h" #include "source_io/module_hs/hs_sparse_io.h" #include "source_base/module_out/filename.h" @@ -40,7 +41,8 @@ void write_Vxc_R(const int nspin, const std::vector>>>>* Hexxc #endif , - const double sparse_thr) + const double sparse_thr, + std::ofstream& ofs_running) { ModuleBase::TITLE("ModuleIO", "write_Vxc_R"); @@ -127,6 +129,8 @@ void write_Vxc_R(const int nspin, all_R_coor, *pv, options); + const std::string spin_tag = ModuleIO::make_spin_tag(is, nspin); + ofs_running << " Write Vxc(R)" << spin_tag << " matrix in NAO basis to file: " << options.filename << std::endl; } } @@ -140,7 +144,7 @@ template void write_Vxc_R( , const std::vector>>>*, const std::vector>>>>* #endif - , const double); + , const double, std::ofstream&); template void write_Vxc_R, double>( const int, const Parallel_Orbitals*, const UnitCell&, Structure_Factor&, surchem&, @@ -151,7 +155,7 @@ template void write_Vxc_R, double>( , const std::vector>>>*, const std::vector>>>>* #endif - , const double); + , const double, std::ofstream&); template void write_Vxc_R, std::complex>( const int, const Parallel_Orbitals*, const UnitCell&, Structure_Factor&, surchem&, @@ -162,6 +166,6 @@ template void write_Vxc_R, std::complex>( , const std::vector>>>*, const std::vector>>>>* #endif - , const double); + , const double, std::ofstream&); } // namespace ModuleIO diff --git a/source/source_io/module_hs/vxc_op_r.h b/source/source_io/module_hs/vxc_op_r.h index a51559f32b1..85f42be1ae5 100644 --- a/source/source_io/module_hs/vxc_op_r.h +++ b/source/source_io/module_hs/vxc_op_r.h @@ -11,6 +11,7 @@ #include "source_hamilt/module_xc/exx_info.h" #include +#include #include #include #include @@ -50,7 +51,8 @@ void write_Vxc_R(const int nspin, const std::vector>>>>* Hexxc #endif , - const double sparse_thr); + const double sparse_thr, + std::ofstream& ofs_running); } // namespace ModuleIO diff --git a/source/source_io/module_output/cube_io.h b/source/source_io/module_output/cube_io.h index d6e9b247731..6f554fef949 100644 --- a/source/source_io/module_output/cube_io.h +++ b/source/source_io/module_output/cube_io.h @@ -1,12 +1,48 @@ #ifndef CUBE_IO_H #define CUBE_IO_H #include "source_cell/unitcell.h" +#include "spin_tag.h" +#include #include +#include class Parallel_Grid; namespace ModuleIO { +/// Build a descriptive label for the data channel at spin index is. +/// For nspin==2 and is>=0: appends " (spin up )" or " (spin down)". +/// For nspin==4 and is>0: returns mag_base with the axis letter appended, +/// e.g. mag_base "magnetization density m" yields "magnetization density mx". +/// For nspin==1, nspin==4 with is==0, or spin-summed calls (is<0): no suffix. +/// +/// @param base label of the scalar channel (is == 0), e.g. "charge density". +/// @param mag_base label of the nspin==4 magnetic channels (is = 1..3), +/// ending with the component letter: "m" for magnetization densities +/// (rho[1..3]), "b" for magnetic fields (the B_xc components of the +/// effective potential). +inline std::string make_data_desc(const std::string& base, const std::string& mag_base, const int is, const int nspin) +{ + std::string desc = base + make_spin_tag(is, nspin); + if (nspin == 4 && is > 0) + { + desc = mag_base; + if (is == 1) + { + desc += "x"; + } + else if (is == 2) + { + desc += "y"; + } + else if (is == 3) + { + desc += "z"; + } + } + return desc; +} + /// read volumetric data from .cube file into the parallel distributed grid. bool read_vdata_palgrid(const Parallel_Grid& pgrid, const int my_rank, @@ -27,7 +63,9 @@ void write_vdata_palgrid(const Parallel_Grid& pgrid, const int precision, const int out_fermi, const bool two_fermi, - const bool reduce_all_pool); + const bool reduce_all_pool, + std::ofstream& ofs_running, + const std::string& data_desc); /// read the full data from a cube file bool read_cube(const std::string& file, diff --git a/source/source_io/module_output/output_log.cpp b/source/source_io/module_output/output_log.cpp index a642f8b7322..8220e16f710 100644 --- a/source/source_io/module_output/output_log.cpp +++ b/source/source_io/module_output/output_log.cpp @@ -1,5 +1,6 @@ #include "output_log.h" +#include "source_estate/fp_energy.h" #include "source_io/module_parameter/parameter.h" #include "source_base/constants.h" #include "source_base/formatter.h" @@ -54,13 +55,23 @@ void output_after_relax(bool conv_ion, bool conv_esolver, std::ofstream& ofs_run } } -void output_efermi(const bool &convergence, double& efermi, std::ofstream& ofs_running) +void output_efermi(const bool &convergence, const elecstate::Efermi& eferm, std::ofstream& ofs_running) { if (convergence && PARAM.inp.out_level != "m") { -// mohan comment out 2025-06-22 -// ofs_running << std::setprecision(16); -// ofs_running << " EFERMI = " << std::setprecision(11) << efermi * ModuleBase::Ry_to_eV << " eV" << std::endl; + if (eferm.two_efermi) + { + ofs_running << " #FERMI ENERGY# " + << std::setprecision(11) << eferm.ef_up * ModuleBase::Ry_to_eV + << " " << eferm.ef_dw * ModuleBase::Ry_to_eV + << " eV # two fermi energies (up, down)" << std::endl; + } + else + { + ofs_running << " #FERMI ENERGY# " + << std::setprecision(11) << eferm.ef * ModuleBase::Ry_to_eV + << " eV" << std::endl; + } } } diff --git a/source/source_io/module_output/output_log.h b/source/source_io/module_output/output_log.h index 8b6c84cbe96..100e6cbf222 100644 --- a/source/source_io/module_output/output_log.h +++ b/source/source_io/module_output/output_log.h @@ -9,6 +9,11 @@ class MDCell; +namespace elecstate +{ +struct Efermi; +} + namespace ModuleIO { @@ -16,19 +21,19 @@ namespace ModuleIO /// @param convergence if is convergence /// @param energy the total energy in Ry /// @param ofs_running the output stream -void output_convergence_after_scf(const bool&convergence, double& energy, std::ofstream& ofs_running = GlobalV::ofs_running); +void output_convergence_after_scf(const bool&convergence, double& energy, std::ofstream& ofs_running); /// @brief output after relaxation /// @param conv_ion if is convergence for ions /// @param conv_esolver if is convergence for electrons /// @param ofs_running the output stream -void output_after_relax(bool conv_ion, bool conv_esolver, std::ofstream& ofs_running = GlobalV::ofs_running); +void output_after_relax(bool conv_ion, bool conv_esolver, std::ofstream& ofs_running); /// @brief output the fermi energy /// @param convergence if is convergence -/// @param efermi +/// @param eferm the fermi energy structure (supports single and two fermi energies) /// @param ofs_running the output stream -void output_efermi(const bool &convergence, double& efermi, std::ofstream& ofs_running = GlobalV::ofs_running); +void output_efermi(const bool &convergence, const elecstate::Efermi& eferm, std::ofstream& ofs_running); /// @brief calculate and output the vacuum level /// We first determine the vacuum direction, then get the vacuum position based on the minimum of charge density, @@ -48,7 +53,7 @@ void output_vacuum_level(const UnitCell* ucell, const int& nrxx, const int& nplane, const int& startz_current, - std::ofstream& ofs_running = GlobalV::ofs_running); + std::ofstream& ofs_running); /// @brief output atomic forces /// @param ofs the output stream diff --git a/source/source_io/module_output/spin_tag.h b/source/source_io/module_output/spin_tag.h new file mode 100644 index 00000000000..19562a55fe7 --- /dev/null +++ b/source/source_io/module_output/spin_tag.h @@ -0,0 +1,27 @@ +#ifndef MODULEIO_SPIN_TAG_H +#define MODULEIO_SPIN_TAG_H + +#include + +namespace ModuleIO +{ + +/// Running-log suffix naming the collinear spin channel of an output. +/// Returns " (spin up )" for nspin==2 and is==0, " (spin down)" for nspin==2 +/// and is==1, and an empty string otherwise (nspin==1, nspin==4, or is<0). +inline std::string make_spin_tag(const int is, const int nspin) +{ + if (nspin != 2 || is < 0) + { + return ""; + } + if (is == 0) + { + return " (spin up )"; + } + return " (spin down)"; +} + +} // namespace ModuleIO + +#endif diff --git a/source/source_io/module_output/write_cube.cpp b/source/source_io/module_output/write_cube.cpp index 1771a28026b..5c3d8e9af54 100644 --- a/source/source_io/module_output/write_cube.cpp +++ b/source/source_io/module_output/write_cube.cpp @@ -22,7 +22,9 @@ void ModuleIO::write_vdata_palgrid(const Parallel_Grid& pgrid, const int precision, const int out_fermi, const bool two_fermi, - const bool reduce_all_pool) + const bool reduce_all_pool, + std::ofstream& ofs_running, + const std::string& data_desc) { ModuleBase::TITLE("ModuleIO", "write_vdata_palgrid"); @@ -173,6 +175,8 @@ void ModuleIO::write_vdata_palgrid(const Parallel_Grid& pgrid, data_xyz_full, precision); + ofs_running << " Write " << data_desc << " to file: " << fn << std::endl; + end = time(nullptr); ModuleBase::GlobalFunc::OUT_TIME("write_vdata_palgrid", start, end); } @@ -220,9 +224,6 @@ void ModuleIO::write_cube(const std::string& file, std::ofstream ofs(file); - // mohan add 2025-09-10 - GlobalV::ofs_running << " Write data to file: " << file << std::endl; - for (int i = 0; i < 2; ++i) { ofs << comment[i] << "\n"; diff --git a/source/source_io/module_wf/get_wf_lcao.cpp b/source/source_io/module_wf/get_wf_lcao.cpp index ea1189acaee..70b1e46784f 100644 --- a/source/source_io/module_wf/get_wf_lcao.cpp +++ b/source/source_io/module_wf/get_wf_lcao.cpp @@ -3,6 +3,7 @@ #include "source_hamilt/module_gint/gint_env_gamma.h" #include "source_hamilt/module_gint/gint_env_k.h" #include "source_io/module_output/cube_io.h" +#include "source_io/module_output/spin_tag.h" #include #include @@ -46,6 +47,11 @@ void Get_wf_lcao::begin_gamma(const UnitCell& ucell, for (int is = 0; is < nspin_; ++is) { psi_gamma_->fix_k(is); + // Collinear (nspin==2) wavefunction cubes carry a spin tag. + const std::string spin_tag = ModuleIO::make_spin_tag(is, nspin_); + const std::string desc_wfc_norm = "wave function (norm)" + spin_tag; + const std::string desc_wfc_real = "wave function (real)" + spin_tag; + const std::string desc_wfc_imag = "wave function (imag)" + spin_tag; ModuleGint::Gint_env_gamma gint_env(psi_gamma_->get_pointer(), ¶_orb_, nbands_, nlocal, wfc_gamma.data()); for (int ib = 0; ib < nbands_; ++ib) { @@ -70,12 +76,7 @@ void Get_wf_lcao::begin_gamma(const UnitCell& ucell, std::stringstream ss_out; ss_out << global_out_dir << ss_file.str(); - std::stringstream ss_info; - ss_info << "Wave func. " << ib + 1 << " spin " << is + 1 << " saved in"; - - ModuleBase::GlobalFunc::OUT(ofs_running, ss_info.str(), ss_file.str()); - - ModuleIO::write_vdata_palgrid(pgrid, wfc_norm.data(), is, nspin_, 0, ss_out.str(), 0.0, &ucell, precision, 0, false, false); + ModuleIO::write_vdata_palgrid(pgrid, wfc_norm.data(), is, nspin_, 0, ss_out.str(), 0.0, &ucell, precision, 0, false, false, ofs_running, desc_wfc_norm); } if (re_im_bands_picked[ib]) @@ -93,7 +94,9 @@ void Get_wf_lcao::begin_gamma(const UnitCell& ucell, precision, 0, false, - false); + false, + ofs_running, + desc_wfc_real); std::stringstream ss_imag; ss_imag << global_out_dir << "wfi" << ib + 1 << "s" << is + 1 << "k1im.cube"; @@ -108,7 +111,9 @@ void Get_wf_lcao::begin_gamma(const UnitCell& ucell, precision, 0, false, - false); + false, + ofs_running, + desc_wfc_imag); } } } @@ -179,10 +184,8 @@ void Get_wf_lcao::begin_k(const UnitCell& ucell, std::stringstream ss_out; ss_out << global_out_dir << ss_file.str(); - std::stringstream ss_info; - ss_info << "Wave func. " << ib + 1 << " spin " << spin_index + 1 << " k-point " << k_number << " saved in"; - - ModuleBase::GlobalFunc::OUT(ofs_running, ss_info.str(), ss_file.str()); + // Only collinear (nspin==2) runs get a spin tag; nspin==4 writes the spinor modulus. + const std::string desc_wfc_norm = "wave function (norm)" + ModuleIO::make_spin_tag(spin_index, nspin_); ModuleIO::write_vdata_palgrid(pgrid, wfc_norm.data(), @@ -195,11 +198,17 @@ void Get_wf_lcao::begin_k(const UnitCell& ucell, precision, 0, false, - false); + false, + ofs_running, + desc_wfc_norm); } if (re_im_bands_picked[ib]) { + // Only collinear (nspin==2) channels get a spin tag; nspin==4 spinor components have none. + const std::string spin_tag = ModuleIO::make_spin_tag(spin_index, nspin_); + const std::string desc_wfc_real = "wave function (real)" + spin_tag; + const std::string desc_wfc_imag = "wave function (imag)" + spin_tag; for (int ipol = 0; ipol < npol; ++ipol) { // Spinors write their two components separately; collinear states use the spin channel. @@ -225,7 +234,9 @@ void Get_wf_lcao::begin_k(const UnitCell& ucell, precision, 0, false, - false); + false, + ofs_running, + desc_wfc_real); std::stringstream ss_imag; ss_imag << global_out_dir << "wfi" << ib + 1 << "s" << component_index + 1 << "k" << k_number << "im.cube"; @@ -240,7 +251,9 @@ void Get_wf_lcao::begin_k(const UnitCell& ucell, precision, 0, false, - false); + false, + ofs_running, + desc_wfc_imag); } } } diff --git a/source/source_io/module_wf/get_wf_pw.cpp b/source/source_io/module_wf/get_wf_pw.cpp index f6c5e54d875..c3d3401c2d9 100644 --- a/source/source_io/module_wf/get_wf_pw.cpp +++ b/source/source_io/module_wf/get_wf_pw.cpp @@ -6,6 +6,7 @@ #include "source_base/module_parallel/para_bridge.h" #include "source_base/tool_quit.h" #include "source_io/module_output/cube_io.h" +#include "source_io/module_output/spin_tag.h" #include #include @@ -83,7 +84,8 @@ void Get_wf_pw::begin(const UnitCell& ucell, const K_Vectors& kv, const std::vector& out_wfc_norm, const std::vector& out_wfc_re_im, - const std::string& global_out_dir) const + const std::string& global_out_dir, + std::ofstream& ofs_running) const { // Resolve global band ownership collectively before validating the selection. const Parallel::ParaBandOutput band_output(psi_.get_nbands(), global_nbands_, Parallel::make_band_world()); @@ -106,7 +108,7 @@ void Get_wf_pw::begin(const UnitCell& ucell, { std::fill(work.values[is].begin(), work.values[is].end(), 0.0); } - write_norm(band, ucell, pgrid, kv, global_out_dir, band_output, &work); + write_norm(band, ucell, pgrid, kv, global_out_dir, band_output, ofs_running, &work); } for (int band = 0; band < global_nbands_; ++band) { @@ -119,7 +121,7 @@ void Get_wf_pw::begin(const UnitCell& ucell, std::fill(work.values[is].begin(), work.values[is].end(), 0.0); std::fill(work.imag[is].begin(), work.imag[is].end(), 0.0); } - write_complex(band, ucell, pgrid, kv, global_out_dir, band_output, &work); + write_complex(band, ucell, pgrid, kv, global_out_dir, band_output, ofs_running, &work); } } @@ -208,6 +210,7 @@ void Get_wf_pw::write_norm(const int band, const K_Vectors& kv, const std::string& out_dir, const Parallel::ParaBandOutput& band_output, + std::ofstream& ofs_running, Workspace* work) const { // Collinear spin channels share the same physical k-point numbering in file names. @@ -221,7 +224,7 @@ void Get_wf_pw::write_norm(const int band, // Wavefunction amplitudes carry the inverse square root of the cell volume. const double scale = std::sqrt(1.0 / ucell.omega); calc_norm(spin_index, scale, work); - write_cube(band, spin_index, k_number, "", ucell, pgrid, out_dir, work->values[spin_index]); + write_cube(band, spin_index, k_number, "", ucell, pgrid, out_dir, work->values[spin_index], ofs_running); } } @@ -232,6 +235,7 @@ void Get_wf_pw::write_complex(const int band, const K_Vectors& kv, const std::string& out_dir, const Parallel::ParaBandOutput& band_output, + std::ofstream& ofs_running, Workspace* work) const { // Collinear spin channels share the same physical k-point numbering in file names. @@ -252,8 +256,8 @@ void Get_wf_pw::write_complex(const int band, { const int wfc_component = work->is_spinor && component == 1 ? 1 : 0; calc_component(work->wfcr[wfc_component], work->phase, scale, &work->values[component], &work->imag[component]); - write_cube(band, component, k_number, "re", ucell, pgrid, out_dir, work->values[component]); - write_cube(band, component, k_number, "im", ucell, pgrid, out_dir, work->imag[component]); + write_cube(band, component, k_number, "re", ucell, pgrid, out_dir, work->values[component], ofs_running); + write_cube(band, component, k_number, "im", ucell, pgrid, out_dir, work->imag[component], ofs_running); } } } @@ -313,11 +317,28 @@ void Get_wf_pw::write_cube(const int band, const UnitCell& ucell, const Parallel_Grid& pgrid, const std::string& out_dir, - const std::vector& values) const + const std::vector& values, + std::ofstream& ofs_running) const { std::stringstream filename; filename << out_dir << "wfi" << band + 1 << "s" << component + 1 << "k" << k_number << part << ".cube"; - ModuleIO::write_vdata_palgrid(pgrid, values.data(), component, nspin_, 0, filename.str(), 0.0, &ucell, 11, 0, false, true); + + std::string data_desc; + if (part == "re") + { + data_desc = "wave function (real)"; + } + else if (part == "im") + { + data_desc = "wave function (imag)"; + } + else + { + data_desc = "wave function (norm)"; + } + // collinear spin channel (nspin=2); spinor components (nspin=4) are not spin channels + data_desc += ModuleIO::make_spin_tag(component, nspin_); + ModuleIO::write_vdata_palgrid(pgrid, values.data(), component, nspin_, 0, filename.str(), 0.0, &ucell, 11, 0, false, true, ofs_running, data_desc); } // Explicit instantiation emits both precisions for each supported device from this .cpp file. diff --git a/source/source_io/module_wf/get_wf_pw.h b/source/source_io/module_wf/get_wf_pw.h index bb370d1732b..5d738f2e7e1 100644 --- a/source/source_io/module_wf/get_wf_pw.h +++ b/source/source_io/module_wf/get_wf_pw.h @@ -44,7 +44,8 @@ class Get_wf_pw const K_Vectors& kv, const std::vector& out_wfc_norm, const std::vector& out_wfc_re_im, - const std::string& global_out_dir) const; + const std::string& global_out_dir, + std::ofstream& ofs_running) const; private: const psi::Psi& psi_; @@ -72,6 +73,7 @@ class Get_wf_pw const K_Vectors& kv, const std::string& out_dir, const Parallel::ParaBandOutput& band_output, + std::ofstream& ofs_running, Workspace* work) const; void write_complex(const int band, const UnitCell& ucell, @@ -79,6 +81,7 @@ class Get_wf_pw const K_Vectors& kv, const std::string& out_dir, const Parallel::ParaBandOutput& band_output, + std::ofstream& ofs_running, Workspace* work) const; void calc_norm(const int spin_index, const double scale, Workspace* work) const; @@ -96,7 +99,8 @@ class Get_wf_pw const UnitCell& ucell, const Parallel_Grid& pgrid, const std::string& out_dir, - const std::vector& values) const; + const std::vector& values, + std::ofstream& ofs_running) const; }; } // namespace ModuleIO diff --git a/source/source_io/test/CMakeLists.txt b/source/source_io/test/CMakeLists.txt index 107ff7de930..41ee440a3d7 100644 --- a/source/source_io/test/CMakeLists.txt +++ b/source/source_io/test/CMakeLists.txt @@ -11,6 +11,11 @@ AddTest( ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp ) +AddTest( + TARGET MODULE_IO_spin_tag + SOURCES test_spin_tag.cpp +) + if(NOT ENABLE_MPI) return() endif() diff --git a/source/source_io/test/outputlog_test.cpp b/source/source_io/test/outputlog_test.cpp index e7fc8886d47..71aa620ab8d 100644 --- a/source/source_io/test/outputlog_test.cpp +++ b/source/source_io/test/outputlog_test.cpp @@ -8,6 +8,7 @@ #include "source_base/constants.h" #include "source_base/global_variable.h" +#include "source_estate/fp_energy.h" #include "source_io/module_output/output_log.h" #ifdef __MPI @@ -105,9 +106,10 @@ TEST(OutputAfterRelaxTest, TestConvergence) TEST(OutputEfermiTest, TestNotConvergence) { bool convergence = false; - double efermi = 1.0; + elecstate::Efermi eferm; + eferm.ef = 1.0; std::ofstream ofs_running("test_output_efermi_noconvergence.txt"); - ModuleIO::output_efermi(convergence, efermi, ofs_running); + ModuleIO::output_efermi(convergence, eferm, ofs_running); ofs_running.close(); std::ifstream ifs_running("test_output_efermi_noconvergence.txt"); @@ -124,10 +126,11 @@ TEST(OutputEfermiTest, TestNotConvergence) { TEST(OutputEfermiTest, TestMOutputLevel) { bool convergence = true; - double efermi = 1.0; + elecstate::Efermi eferm; + eferm.ef = 1.0; TestParameters::input().out_level = "m"; // Setting output level to "m" std::ofstream ofs_running("test_output_efermi_m_outputlevel.txt"); - ModuleIO::output_efermi(convergence, efermi, ofs_running); + ModuleIO::output_efermi(convergence, eferm, ofs_running); ofs_running.close(); std::ifstream ifs_running("test_output_efermi_m_outputlevel.txt"); @@ -142,6 +145,55 @@ TEST(OutputEfermiTest, TestMOutputLevel) { std::remove("test_output_efermi_m_outputlevel.txt"); } +TEST(OutputEfermiTest, TestSingleFermi) { + bool convergence = true; + elecstate::Efermi eferm; + eferm.two_efermi = false; + eferm.ef = 0.4; // in Ry + TestParameters::input().out_level = "i"; + std::ofstream ofs_running("test_output_efermi_single.txt"); + ModuleIO::output_efermi(convergence, eferm, ofs_running); + ofs_running.close(); + + std::ifstream ifs_running("test_output_efermi_single.txt"); + std::stringstream ss; + ss << ifs_running.rdbuf(); + std::string file_content = ss.str(); + ifs_running.close(); + + EXPECT_THAT(file_content, testing::HasSubstr("#FERMI ENERGY# ")); + EXPECT_THAT(file_content, testing::HasSubstr("eV")); + // 0.4 Ry * 13.605698 eV/Ry = 5.4422792 eV + EXPECT_THAT(file_content, testing::HasSubstr("5.4422792")); + std::remove("test_output_efermi_single.txt"); +} + +TEST(OutputEfermiTest, TestTwoFermi) { + bool convergence = true; + elecstate::Efermi eferm; + eferm.two_efermi = true; + eferm.ef_up = 0.3; // in Ry + eferm.ef_dw = 0.5; // in Ry + TestParameters::input().out_level = "i"; + std::ofstream ofs_running("test_output_efermi_two.txt"); + ModuleIO::output_efermi(convergence, eferm, ofs_running); + ofs_running.close(); + + std::ifstream ifs_running("test_output_efermi_two.txt"); + std::stringstream ss; + ss << ifs_running.rdbuf(); + std::string file_content = ss.str(); + ifs_running.close(); + + EXPECT_THAT(file_content, testing::HasSubstr("#FERMI ENERGY# ")); + EXPECT_THAT(file_content, testing::HasSubstr("two fermi energies (up, down)")); + // 0.3 Ry * 13.605698 eV/Ry = 4.0817094 eV + EXPECT_THAT(file_content, testing::HasSubstr("4.0817094")); + // 0.5 Ry * 13.605698 eV/Ry = 6.802849 eV + EXPECT_THAT(file_content, testing::HasSubstr("6.802849")); + std::remove("test_output_efermi_two.txt"); +} + UnitCell::UnitCell() { ntype = 1; diff --git a/source/source_io/test/test_spin_tag.cpp b/source/source_io/test/test_spin_tag.cpp new file mode 100644 index 00000000000..78146dd64b9 --- /dev/null +++ b/source/source_io/test/test_spin_tag.cpp @@ -0,0 +1,28 @@ +#include "source_io/module_output/spin_tag.h" + +#include "gtest/gtest.h" + +#include + +/** + * Tests for ModuleIO::make_spin_tag, the running-log suffix shared by the + * H(R), DM, Vxc(R), band and cube writers. + */ + +TEST(SpinTagTest, CollinearSpinChannels) +{ + EXPECT_EQ(ModuleIO::make_spin_tag(0, 2), " (spin up )"); + EXPECT_EQ(ModuleIO::make_spin_tag(1, 2), " (spin down)"); +} + +TEST(SpinTagTest, NoTagForSingleSpinOrNoncollinear) +{ + EXPECT_EQ(ModuleIO::make_spin_tag(0, 1), ""); + EXPECT_EQ(ModuleIO::make_spin_tag(0, 4), ""); + EXPECT_EQ(ModuleIO::make_spin_tag(3, 4), ""); +} + +TEST(SpinTagTest, NoTagForSpinSummedCall) +{ + EXPECT_EQ(ModuleIO::make_spin_tag(-1, 2), ""); +} diff --git a/source/source_io/test/write_bands_test.cpp b/source/source_io/test/write_bands_test.cpp index a892f015f7c..4fdcbdbb7e5 100644 --- a/source/source_io/test/write_bands_test.cpp +++ b/source/source_io/test/write_bands_test.cpp @@ -59,6 +59,9 @@ namespace using BandTopology = std::pair; +// Dummy output stream for running-log writes in tests; discards output. +std::ofstream ofs_running; + // GoogleTest evaluates this generator during InitGoogleTest, after MPI_Init. std::vector band_topologies() { @@ -301,7 +304,7 @@ class BandOutputTest : public testing::TestWithParam { SCOPED_TRACE(distributed ? "distributed" : "replicated"); prepare(nspin, 3, distributed); - ModuleIO::write_bands(input_, ekb_, kv_); + ModuleIO::write_bands(input_, ekb_, kv_, ofs_running); if (nspin == 2) { expect_output("bands1.txt", 0, 3, 0.0, 8, 1); @@ -363,11 +366,11 @@ TEST_P(BandOutputTest, OverwriteAndAppend) prepare(1, 3, true); seed_file("band.txt", "stale output\n"); Parallel::barrier(*world_); - ModuleIO::write_bands(input_, ekb_, kv_); + ModuleIO::write_bands(input_, ekb_, kv_, ofs_running); expect_output("band.txt", 0, 3, 0.0, 8, 1); Parallel::barrier(*world_); parameters_->set_output(directory_, true); - ModuleIO::write_bands(input_, ekb_, kv_); + ModuleIO::write_bands(input_, ekb_, kv_, ofs_running); expect_output("band.txt", 0, 3, 0.0, 8, 2); } @@ -375,11 +378,11 @@ TEST_P(BandOutputTest, OutputDisabled) { prepare(1, 3, true); input_.out_band[0] = 0; - ModuleIO::write_bands(input_, ekb_, kv_); + ModuleIO::write_bands(input_, ekb_, kv_, ofs_running); expect_missing("band.txt"); seed_file("band.txt", "keep this content\n"); Parallel::barrier(*world_); - ModuleIO::write_bands(input_, ekb_, kv_); + ModuleIO::write_bands(input_, ekb_, kv_, ofs_running); if (world_->rank() == 0) { EXPECT_EQ(read_text("band.txt"), "keep this content\n"); @@ -410,11 +413,32 @@ TEST_P(BandOutputTest, EmptyBandShard) { EXPECT_EQ(ekb_.nc, 0); } - ModuleIO::write_bands(input_, ekb_, kv_); + ModuleIO::write_bands(input_, ekb_, kv_, ofs_running); expect_output("bands1.txt", 0, 1, 0.0, 8, 1); expect_output("bands2.txt", 1, 1, 0.0, 8, 1); } +TEST_P(BandOutputTest, RunningLogNamesSpinChannels) +{ + prepare(2, 3, false); + const std::string log_name = "write_bands_running_p" + std::to_string(world_->size()) + "_r" + + std::to_string(world_->rank()) + ".log"; + std::ofstream log(log_name); + EXPECT_TRUE(log.is_open()); + ModuleIO::write_bands(input_, ekb_, kv_, log); + log.close(); + + std::ifstream in(log_name); + const std::string text((std::istreambuf_iterator(in)), std::istreambuf_iterator()); + in.close(); + std::remove(log_name.c_str()); + + const std::string up_line = " Write eigenvalues (spin up ) to file: " + directory_ + "bands1.txt"; + const std::string down_line = " Write eigenvalues (spin down) to file: " + directory_ + "bands2.txt"; + EXPECT_NE(text.find(up_line), std::string::npos) << text; + EXPECT_NE(text.find(down_line), std::string::npos) << text; +} + INSTANTIATE_TEST_SUITE_P(Topology, BandOutputTest, testing::ValuesIn(band_topologies()), topology_name); } // namespace diff --git a/source/source_io/test_serial/rho_io_test.cpp b/source/source_io/test_serial/rho_io_test.cpp index 7bc3471c7f8..38b226fc994 100644 --- a/source/source_io/test_serial/rho_io_test.cpp +++ b/source/source_io/test_serial/rho_io_test.cpp @@ -97,7 +97,7 @@ TEST_F(RhoIOTest, Write) ucell->atoms[0].ncpp.zv = 4; Parallel_Grid pgrid(nx, ny, nz, nz, nrxx, nz, 1); ModuleIO::read_vdata_palgrid(pgrid, my_rank, ofs_running, "support/chg.cube", rho[0], ucell->nat); - ModuleIO::write_vdata_palgrid(pgrid, rho[0], 0, nspin, 0, "test_write_vdata_palgrid.cube", 0.461002, ucell, 11, 1, false, false); + ModuleIO::write_vdata_palgrid(pgrid, rho[0], 0, nspin, 0, "test_write_vdata_palgrid.cube", 0.461002, ucell, 11, 1, false, false, ofs_running, "charge density"); std::ifstream ifs1("test_write_vdata_palgrid.cube", std::ifstream::binary | std::ifstream::ate); std::ifstream ifs2("support/chg.cube", std::ifstream::binary | std::ifstream::ate); EXPECT_EQ(ifs1.tellg(), ifs2.tellg()); diff --git a/source/source_lcao/module_lr/utils/exciton_plotter.cpp b/source/source_lcao/module_lr/utils/exciton_plotter.cpp index b7c4483198f..037ba0586d7 100644 --- a/source/source_lcao/module_lr/utils/exciton_plotter.cpp +++ b/source/source_lcao/module_lr/utils/exciton_plotter.cpp @@ -496,8 +496,8 @@ void ExcitonPlotter::plot_average_density(const int istate, const std::string for (int is = 0; is < this->nspin_x; ++is) { - const std::string filename = this->output_dir_ + "Exciton_avg_" + type + "_state" + std::to_string(istate) - + "_spin" + std::to_string(is) + ".cube"; + const std::string filename = this->output_dir_ + "exc_" + type + "_st" + std::to_string(istate + 1) + + "_s" + std::to_string(is + 1) + ".cube"; ModuleIO::write_vdata_palgrid(this->Pgrid, rho_result[is], is, @@ -509,7 +509,9 @@ void ExcitonPlotter::plot_average_density(const int istate, const std::string 8, 0, false, /*two_fermi*/ - false); + false, + GlobalV::ofs_running, + "exciton density"); } LR_Util::_deallocate_2order_nested_ptr(rho_result, this->nspin_x); } @@ -586,7 +588,7 @@ void ExcitonPlotter::plot_average_slice(const int istate, } const std::string filename - = this->output_dir_ + "Exciton_avg_" + type + "_slice_state" + std::to_string(istate) + ".dat"; + = this->output_dir_ + "exc_" + type + "_slice_st" + std::to_string(istate + 1) + ".dat"; write_slice_data(this->ucell, geom, density, @@ -791,7 +793,7 @@ void ExcitonPlotter::plot_cond_slice(const int istate, std::cout << "Bloch-sum evaluation complete (cached " << cell_res << "x" << cell_res << " home-cell positions)." << std::endl; const std::string filename - = this->output_dir_ + "Exciton_cond_" + type + "_slice_state" + std::to_string(istate) + ".dat"; + = this->output_dir_ + "exc_cond_" + type + "_slice_st" + std::to_string(istate + 1) + ".dat"; write_slice_data(this->ucell, geom, density, diff --git a/tests/01_PW/scf_out_ldos/README b/tests/01_PW/scf_out_ldos/README index 239ee25dae8..6236cd82e7e 100644 --- a/tests/01_PW/scf_out_ldos/README +++ b/tests/01_PW/scf_out_ldos/README @@ -1,3 +1,3 @@ 1. SG15 ONCV pseudopotential 2. symmetry=on -3. test LDOS with 2 eV, output LDOS_2eV.cube +3. test LDOS with 2 eV, output ldos_2ev.cube diff --git a/tests/01_PW/scf_out_ldos/LDOS.cube.ref b/tests/01_PW/scf_out_ldos/ldos.cube.ref similarity index 100% rename from tests/01_PW/scf_out_ldos/LDOS.cube.ref rename to tests/01_PW/scf_out_ldos/ldos.cube.ref diff --git a/tests/01_PW/scf_out_ldos/result.ref b/tests/01_PW/scf_out_ldos/result.ref index 60735dc1920..588715b328c 100644 --- a/tests/01_PW/scf_out_ldos/result.ref +++ b/tests/01_PW/scf_out_ldos/result.ref @@ -1,6 +1,6 @@ etotref 13.59535887413501 etotperatomref 6.7976794371 -LDOS.cube_pass 0 +ldos.cube_pass 0 pointgroupref T_d spacegroupref O_h nksibzref 1 diff --git a/tests/02_NAO_Gamma/scf_out_hk_spin2/hrs1_nao.csr.ref b/tests/02_NAO_Gamma/scf_out_hk_spin2/hrs1_nao.csr.ref index 50e89733b66..698669eeffd 100644 --- a/tests/02_NAO_Gamma/scf_out_hk_spin2/hrs1_nao.csr.ref +++ b/tests/02_NAO_Gamma/scf_out_hk_spin2/hrs1_nao.csr.ref @@ -1,7 +1,7 @@ --- Ionic Step 1 --- # print H matrix in real space H(R) 2 # number of spin directions - 1 # spin index + 1 # spin index, E_Fermi = 14.2079 eV 8 # number of localized basis 1 # number of Bravais lattice vector R @@ -35,3 +35,4 @@ 0 4 1 5 2 6 3 7 0 4 1 5 2 6 3 7 # CSR row pointers 0 2 4 6 8 10 12 14 16 + diff --git a/tests/02_NAO_Gamma/scf_out_hk_spin2/hrs2_nao.csr.ref b/tests/02_NAO_Gamma/scf_out_hk_spin2/hrs2_nao.csr.ref index 41706c30819..86c0a285611 100644 --- a/tests/02_NAO_Gamma/scf_out_hk_spin2/hrs2_nao.csr.ref +++ b/tests/02_NAO_Gamma/scf_out_hk_spin2/hrs2_nao.csr.ref @@ -1,7 +1,7 @@ --- Ionic Step 1 --- # print H matrix in real space H(R) 2 # number of spin directions - 2 # spin index + 2 # spin index, E_Fermi = 14.2079 eV 8 # number of localized basis 1 # number of Bravais lattice vector R @@ -35,3 +35,4 @@ 0 4 1 5 2 6 3 7 0 4 1 5 2 6 3 7 # CSR row pointers 0 2 4 6 8 10 12 14 16 + diff --git a/tests/03_NAO_multik/CASES_CPU.txt b/tests/03_NAO_multik/CASES_CPU.txt index bd57f5d9d02..a85a0938ec3 100644 --- a/tests/03_NAO_multik/CASES_CPU.txt +++ b/tests/03_NAO_multik/CASES_CPU.txt @@ -43,6 +43,7 @@ scf_out_mul_spin4 scf_out_mul_nupdw scf_out_elf scf_out_qo +scf_out_hsr_spin2 nscf_out_dos nscf_out_band_pband nscf_out_pot1 diff --git a/tests/03_NAO_multik/CASES_GPU.txt b/tests/03_NAO_multik/CASES_GPU.txt index 1453c8cf282..ef50c4fff3e 100644 --- a/tests/03_NAO_multik/CASES_GPU.txt +++ b/tests/03_NAO_multik/CASES_GPU.txt @@ -42,6 +42,7 @@ scf_out_mul_spin4 scf_out_mul_nupdw scf_out_elf scf_out_qo +scf_out_hsr_spin2 nscf_out_dos nscf_out_band_pband nscf_out_pot1 diff --git a/tests/03_NAO_multik/nscf_out_hsr_tr_rr/hrs1_nao.csr.ref b/tests/03_NAO_multik/nscf_out_hsr_tr_rr/hrs1_nao.csr.ref index 8e9584bfd1a..f4779ab1e67 100644 --- a/tests/03_NAO_multik/nscf_out_hsr_tr_rr/hrs1_nao.csr.ref +++ b/tests/03_NAO_multik/nscf_out_hsr_tr_rr/hrs1_nao.csr.ref @@ -1,7 +1,7 @@ --- Ionic Step 1 --- # print H matrix in real space H(R) 1 # number of spin directions - 1 # spin index + 1 # spin index, E_Fermi = -4.67919 eV 6 # number of localized basis 11 # number of Bravais lattice vector R diff --git a/tests/03_NAO_multik/scf_out_hsr/hrs1_nao.csr.ref b/tests/03_NAO_multik/scf_out_hsr/hrs1_nao.csr.ref index ff316f89108..84cdba91759 100644 --- a/tests/03_NAO_multik/scf_out_hsr/hrs1_nao.csr.ref +++ b/tests/03_NAO_multik/scf_out_hsr/hrs1_nao.csr.ref @@ -1,7 +1,7 @@ --- 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.4252 eV 8 # number of localized basis 13 # number of Bravais lattice vector R diff --git a/tests/03_NAO_multik/scf_out_hsr_spin2/hrs1_nao.csr.ref b/tests/03_NAO_multik/scf_out_hsr_spin2/hrs1_nao.csr.ref index de17e412592..f6c84569230 100644 --- a/tests/03_NAO_multik/scf_out_hsr_spin2/hrs1_nao.csr.ref +++ b/tests/03_NAO_multik/scf_out_hsr_spin2/hrs1_nao.csr.ref @@ -1,7 +1,7 @@ --- Ionic Step 1 --- # print H matrix in real space H(R) 2 # number of spin directions - 1 # spin index + 1 # spin index, E_Fermi = 10.3798 eV 8 # number of localized basis 13 # number of Bravais lattice vector R diff --git a/tests/03_NAO_multik/scf_out_hsr_spin2/hrs2_nao.csr.ref b/tests/03_NAO_multik/scf_out_hsr_spin2/hrs2_nao.csr.ref index d4aa39c5b68..c0d22344976 100644 --- a/tests/03_NAO_multik/scf_out_hsr_spin2/hrs2_nao.csr.ref +++ b/tests/03_NAO_multik/scf_out_hsr_spin2/hrs2_nao.csr.ref @@ -1,7 +1,7 @@ --- Ionic Step 1 --- # print H matrix in real space H(R) 2 # number of spin directions - 2 # spin index + 2 # spin index, E_Fermi = 10.3798 eV 8 # number of localized basis 13 # number of Bravais lattice vector R diff --git a/tests/03_NAO_multik/scf_out_hsr_spin4/hrs1_nao.csr.ref b/tests/03_NAO_multik/scf_out_hsr_spin4/hrs1_nao.csr.ref index 06b276cc3ae..53e9c725198 100644 --- a/tests/03_NAO_multik/scf_out_hsr_spin4/hrs1_nao.csr.ref +++ b/tests/03_NAO_multik/scf_out_hsr_spin4/hrs1_nao.csr.ref @@ -1,7 +1,7 @@ --- Ionic Step 1 --- # print H matrix in real space H(R) 1 # number of spin directions - 1 # spin index + 1 # spin index, E_Fermi = -4.70347 eV 26 # number of localized basis 7 # number of Bravais lattice vector R diff --git a/tests/integrate/validation_tools/CompareFile.py b/tests/integrate/validation_tools/CompareFile.py index b8d55d2957e..4aa6cb13547 100755 --- a/tests/integrate/validation_tools/CompareFile.py +++ b/tests/integrate/validation_tools/CompareFile.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -import os,sys +import sys usage = ''' python CompareFile.py file1 file2 [accuracy] [-abs 0] [-com_type 1] @@ -54,10 +54,13 @@ def is_two_data_diff(data1,data2,epsilon,compareAbsoluteDiff): return False def ReadFile(file1,lines): - if os.path.isfile(file1): - with open(file1) as f1: + # Open the path directly instead of pre-checking with os.path.isfile, + # so that non-regular files such as /dev/fd/ pipes produced by bash + # process substitution <(...) are also accepted. + try: + with open(file1) as f1: for line in f1.readlines(): lines.append(line) - else: + except OSError: print("Error: can not find file %s" % file1) sys.exit(1) diff --git a/tests/integrate/validation_tools/props_cube.sh b/tests/integrate/validation_tools/props_cube.sh index 593dc2d469a..fb048df4f93 100644 --- a/tests/integrate/validation_tools/props_cube.sh +++ b/tests/integrate/validation_tools/props_cube.sh @@ -76,8 +76,8 @@ fi #--------------------------------------- if ! test -z "$has_ldos" && [ $has_ldos == 1 ]; then stm_bias=$(get_input_key_value "stm_bias" "OUT.autotest/INPUT.info") - python3 $COMPARE_SCRIPT LDOS.cube.ref OUT.autotest/LDOS_"$stm_bias"eV.cube 8 - echo "LDOS.cube_pass $?" >> $props_result_file + python3 $COMPARE_SCRIPT ldos.cube.ref OUT.autotest/ldos_"$stm_bias"ev.cube 8 + echo "ldos.cube_pass $?" >> $props_result_file fi #--------------------------------------- diff --git a/tests/integrate/validation_tools/props_mat.sh b/tests/integrate/validation_tools/props_mat.sh index f114b5bbbc4..1ec823fc41a 100644 --- a/tests/integrate/validation_tools/props_mat.sh +++ b/tests/integrate/validation_tools/props_mat.sh @@ -14,6 +14,29 @@ # run_mat_props() - everything from S(k) through dH(k) term matrices # run_mat_dm_props() - out_dm block (between wfc_lcao and mulliken) +# Compare an H(R) CSR file with its reference. +# The "# spin index, E_Fermi = eV" header annotation is printed with +# setprecision(6), so the matrix body is compared at 1e-8 with the annotation +# stripped, and the E_Fermi values are compared separately at 1e-3 eV. +# A missing file is reported as a failure (process substitution would +# otherwise hand CompareFile.py an empty stream and pass silently). +compare_hr_csr(){ + local ref=$1 + local cal=$2 + local strip='s/, E_Fermi = [0-9.eE+-]* eV//' + local pick='s/.*E_Fermi = \([0-9.eE+-]*\) eV.*/\1/p' + if ! test -f "$ref" || ! test -f "$cal"; then + echo "Error: can not find file $ref or $cal" + return 1 + fi + python3 $COMPARE_SCRIPT <(sed "$strip" "$ref") <(sed "$strip" "$cal") 8 || return 1 + if [ "$(sed -n "$pick" "$ref" | wc -l)" -ne "$(sed -n "$pick" "$cal" | wc -l)" ]; then + echo "Error: E_Fermi annotation count differs between $ref and $cal" + return 1 + fi + python3 $COMPARE_SCRIPT <(sed -n "$pick" "$ref") <(sed -n "$pick" "$cal") 3 +} + run_mat_dm1_props(){ #------------------------------- @@ -234,10 +257,10 @@ fi #----------------------------------- #echo $has_hs2 if ! test -z "$has_hs2" && [ $has_hs2 == 1 ]; then - python3 $COMPARE_SCRIPT hrs1_nao.csr.ref OUT.autotest/hrs1_nao.csr 8 + compare_hr_csr hrs1_nao.csr.ref OUT.autotest/hrs1_nao.csr echo "CompareHR_pass $?" >>$props_result_file if ! test -z "$nspin" && [ "$nspin" -eq 2 ]; then - python3 $COMPARE_SCRIPT hrs2_nao.csr.ref OUT.autotest/hrs2_nao.csr 8 + compare_hr_csr hrs2_nao.csr.ref OUT.autotest/hrs2_nao.csr echo "CompareHR2_pass $?" >>$props_result_file fi python3 $COMPARE_SCRIPT sr_nao.csr.ref OUT.autotest/sr_nao.csr 8 diff --git a/tools/02_postprocessing/stm/stm.py b/tools/02_postprocessing/stm/stm.py index 8bc3df71378..2e978ea71f2 100644 --- a/tools/02_postprocessing/stm/stm.py +++ b/tools/02_postprocessing/stm/stm.py @@ -24,7 +24,7 @@ def read_ldos(self, bias): if(abs(bias) < 1e-5): bias = 0.0 - filename = f'{self.dirname}/LDOS_{bias:.5g}eV.cube' + filename = f'{self.dirname}/ldos_{bias:.5g}ev.cube' print('read in ' + filename) self.ldos, self.atoms = read_cube_data(filename) self.cell = self.atoms.cell