diff --git a/docs/advanced/output_files/output-specification.md b/docs/advanced/output_files/output-specification.md index 27a9d8376ec..834d18b6870 100644 --- a/docs/advanced/output_files/output-specification.md +++ b/docs/advanced/output_files/output-specification.md @@ -49,11 +49,18 @@ All output file naming conventions can be found in the online documentation (wit | `s12` | Non-collinear spin calculation | | `k#` | k-point index (e.g., `k1`, `k2`) | | `g#` | Ionic step index for relax/md (e.g., `g1`, `g2`) | +| `ini` | Initial state (before electronic iteration), used before or after `g#` | **Important:** - All index numbers start from 1 (not 0) - For Gamma-only algorithm in LCAO, no `k` index is included - Overlap matrix `s` does not distinguish spin, so only one matrix is output +- For initial charge density output (`out_chg = 2`): + - `out_freq_ion = 0`: `chg_ini.cube` (nspin=1) or `chgs{#}_ini.cube` (nspin=2/4) - output at every step (overwrite same file) + - `out_freq_ion > 0`: `chgg{#}_ini.cube` (nspin=1) or `chgs{#}g{#}_ini.cube` (nspin=2/4) - output every out_freq_ion steps +- For initial potential output (`out_pot = 3`): + - `out_freq_ion = 0`: `pot_ini.cube` (nspin=1) or `pots{#}_ini.cube` (nspin=2/4) - output at every step (overwrite same file) + - `out_freq_ion > 0`: `potg{#}_ini.cube` (nspin=1) or `pots{#}g{#}_ini.cube` (nspin=2/4) - output every out_freq_ion steps ### 2.2 Examples @@ -62,7 +69,22 @@ All output file naming conventions can be found in the online documentation (wit | `chgs1.cube` | Charge density, spin 1 | | `chgs2.cube` | Charge density, spin 2 | | `chgs3.cube` | Charge density, spin 3 (non-collinear with SOC) | +| `chg_ini.cube` | Initial charge density (out_freq_ion=0, nspin=1) | +| `chgs1_ini.cube` | Initial charge density (out_freq_ion=0, spin 1, nspin=2/4) | +| `chgg1_ini.cube` | Initial charge density, geometry step 1 (nspin=1) | +| `chgs1g1_ini.cube` | Initial charge density, spin 1, geometry step 1 (nspin=2/4) | +| `pot_ini.cube` | Initial potential (out_freq_ion=0, nspin=1) | +| `pots1_ini.cube` | Initial potential (out_freq_ion=0, spin 1, nspin=2/4) | +| `potg1_ini.cube` | Initial potential, geometry step 1 (nspin=1) | +| `pots1g1_ini.cube` | Initial potential, spin 1, geometry step 1 (nspin=2/4) | | `pots1.cube` | Local potential, spin 1 | +| `elftot.cube` | ELF, total (nspin=1/4) | +| `elfs1.cube` | ELF, spin 1 (nspin=2) | +| `elfs2.cube` | ELF, spin 2 (nspin=2) | +| `elftot.cube` | ELF, total (nspin=2) | +| `elftotg1.cube` | ELF, total, geometry step 1 (nspin=1/4) | +| `elfs1g1.cube` | ELF, spin 1, geometry step 1 (nspin=2) | +| `elftotg1.cube` | ELF, total, geometry step 1 (nspin=2) | | `eig_occ.txt` | Eigenvalues and occupations | | `doss1g1_nao.txt` | DOS, spin 1, geometry step 1, NAO basis | | `wf_pw.dat` | Wavefunction, plane wave basis | @@ -79,8 +101,25 @@ All output file naming conventions can be found in the online documentation (wit | `band.txt` | Band structure | | `chgs1.cube`, `chgs2.cube` | Charge density (spin 1, spin 2) | | `chg.cube` | Total charge density | +| Initial charge density (out_chg=2) | | +| `chg_ini.cube` | Initial charge density (out_freq_ion=0, nspin=1) | +| `chgs{#}_ini.cube` | Initial charge density (out_freq_ion=0, spin {#}, nspin=2/4) | +| `chgg{#}_ini.cube` | Initial charge density (out_freq_ion>0, geometry step {#}, nspin=1) | +| `chgs{#}g{#}_ini.cube` | Initial charge density (out_freq_ion>0, spin {#}, geometry step {#}, nspin=2/4) | +| Initial potential (out_pot=3) | | +| `pot_ini.cube` | Initial potential (out_freq_ion=0, nspin=1) | +| `pots{#}_ini.cube` | Initial potential (out_freq_ion=0, spin {#}, nspin=2/4) | +| `potg{#}_ini.cube` | Initial potential (out_freq_ion>0, geometry step {#}, nspin=1) | +| `pots{#}g{#}_ini.cube` | Initial potential (out_freq_ion>0, spin {#}, geometry step {#}, nspin=2/4) | | `taus1.cube`, `taus2.cube` | Kinetic energy density (tau) | | `pots1.cube`, `pots2.cube` | Local potential | +| ELF (out_elf) | | +| `elftot.cube` | ELF, total (nspin=1/4) | +| `elfs1.cube`, `elfs2.cube` | ELF, spin 1/2 (nspin=2) | +| `elftot.cube` | ELF, total (nspin=2) | +| `elftotg{#}.cube` | ELF, total, geometry step {#} (nspin=1/4) | +| `elfs{#}g{#}.cube` | ELF, spin {#}, geometry step {#} (nspin=2) | +| `elftotg{#}.cube` | ELF, total, geometry step {#} (nspin=2) | ## 3. File Format Standards diff --git a/examples/18_md/01_lcao_gamma_Si8/INPUT b/examples/18_md/01_LCAO_NVE_Si8/INPUT similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/INPUT rename to examples/18_md/01_LCAO_NVE_Si8/INPUT diff --git a/examples/18_md/01_lcao_gamma_Si8/INPUT_2 b/examples/18_md/01_LCAO_NVE_Si8/INPUT_2 similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/INPUT_2 rename to examples/18_md/01_LCAO_NVE_Si8/INPUT_2 diff --git a/examples/18_md/01_lcao_gamma_Si8/INPUT_3 b/examples/18_md/01_LCAO_NVE_Si8/INPUT_3 similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/INPUT_3 rename to examples/18_md/01_LCAO_NVE_Si8/INPUT_3 diff --git a/examples/18_md/01_lcao_gamma_Si8/INPUT_4 b/examples/18_md/01_LCAO_NVE_Si8/INPUT_4 similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/INPUT_4 rename to examples/18_md/01_LCAO_NVE_Si8/INPUT_4 diff --git a/examples/18_md/01_lcao_gamma_Si8/INPUT_5 b/examples/18_md/01_LCAO_NVE_Si8/INPUT_5 similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/INPUT_5 rename to examples/18_md/01_LCAO_NVE_Si8/INPUT_5 diff --git a/examples/18_md/01_lcao_gamma_Si8/INPUT_6 b/examples/18_md/01_LCAO_NVE_Si8/INPUT_6 similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/INPUT_6 rename to examples/18_md/01_LCAO_NVE_Si8/INPUT_6 diff --git a/examples/18_md/01_lcao_gamma_Si8/INPUT_7 b/examples/18_md/01_LCAO_NVE_Si8/INPUT_7 similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/INPUT_7 rename to examples/18_md/01_LCAO_NVE_Si8/INPUT_7 diff --git a/examples/18_md/01_lcao_gamma_Si8/KPT b/examples/18_md/01_LCAO_NVE_Si8/KPT similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/KPT rename to examples/18_md/01_LCAO_NVE_Si8/KPT diff --git a/examples/18_md/01_lcao_gamma_Si8/README b/examples/18_md/01_LCAO_NVE_Si8/README similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/README rename to examples/18_md/01_LCAO_NVE_Si8/README diff --git a/examples/18_md/01_lcao_gamma_Si8/STRU b/examples/18_md/01_LCAO_NVE_Si8/STRU similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/STRU rename to examples/18_md/01_LCAO_NVE_Si8/STRU diff --git a/examples/18_md/01_lcao_gamma_Si8/run.sh b/examples/18_md/01_LCAO_NVE_Si8/run.sh similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/run.sh rename to examples/18_md/01_LCAO_NVE_Si8/run.sh diff --git a/examples/18_md/01_lcao_gamma_Si8/INPUT_0 b/examples/18_md/01_lcao_gamma_Si8/INPUT_0 deleted file mode 100644 index 5827a8c0178..00000000000 --- a/examples/18_md/01_lcao_gamma_Si8/INPUT_0 +++ /dev/null @@ -1,33 +0,0 @@ -INPUT_PARAMETERS -#Parameters (1.General) -suffix Si_nve -calculation md -nbands 20 -symmetry 0 -pseudo_dir ../../../tests/PP_ORB -orbital_dir ../../../tests/PP_ORB - -#Parameters (2.Iteration) -ecutwfc 30 -scf_thr 1e-5 -scf_nmax 100 - -#Parameters (3.Basis) -basis_type lcao -ks_solver genelpa -gamma_only 1 - -#Parameters (4.Smearing) -smearing_method gaussian -smearing_sigma 0.001 - -#Parameters (5.Mixing) -mixing_type broyden -mixing_beta 0.3 -chg_extrap second-order - -#Parameters (6.MD) -md_type nve -md_nstep 10 -md_dt 1 -md_tfirst 300 diff --git a/examples/18_md/01_lcao_gamma_Si8/INPUT_1 b/examples/18_md/02_LCAO_NVT_Si8/INPUT similarity index 73% rename from examples/18_md/01_lcao_gamma_Si8/INPUT_1 rename to examples/18_md/02_LCAO_NVT_Si8/INPUT index c49ff60e925..c07ca1f9f4a 100644 --- a/examples/18_md/01_lcao_gamma_Si8/INPUT_1 +++ b/examples/18_md/02_LCAO_NVT_Si8/INPUT @@ -3,6 +3,7 @@ INPUT_PARAMETERS suffix Si_nhc_nvt calculation md nbands 20 +nspin 2 symmetry 0 pseudo_dir ../../../tests/PP_ORB orbital_dir ../../../tests/PP_ORB @@ -15,7 +16,7 @@ scf_nmax 100 #Parameters (3.Basis) basis_type lcao ks_solver genelpa -gamma_only 1 +gamma_only 0 #Parameters (4.Smearing) smearing_method gaussian @@ -28,8 +29,15 @@ chg_extrap second-order #Parameters (6.MD) md_type nvt -md_nstep 10 +md_nstep 5 md_dt 1 md_tfirst 300 md_tfreq 1 -md_tchain 1 \ No newline at end of file +md_tchain 1 + +# output initial charge density +out_chg 2 4 +out_pot 3 4 +out_elf 1 4 +out_dipole 1 +out_freq_ion 2 diff --git a/examples/18_md/02_LCAO_NVT_Si8/KPT b/examples/18_md/02_LCAO_NVT_Si8/KPT new file mode 100644 index 00000000000..f5f7f4ec34c --- /dev/null +++ b/examples/18_md/02_LCAO_NVT_Si8/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 2 2 0 0 0 diff --git a/examples/18_md/02_LCAO_NVT_Si8/STRU b/examples/18_md/02_LCAO_NVT_Si8/STRU new file mode 100644 index 00000000000..795530e8cf2 --- /dev/null +++ b/examples/18_md/02_LCAO_NVT_Si8/STRU @@ -0,0 +1,28 @@ +ATOMIC_SPECIES +Si 28.085 Si_ONCV_PBE-1.0.upf + +NUMERICAL_ORBITAL +Si_gga_8au_60Ry_2s2p1d.orb + +LATTICE_CONSTANT +1.8897270 # 1 Angstrom = 1.8897270 bohr + +LATTICE_VECTORS +5.43090 0.00000 0.00000 +0.00000 5.43090 0.00000 +0.00000 0.00000 5.43090 + +ATOMIC_POSITIONS +Direct + +Si +0.0 +8 +0.000 0.000 0.000 1 1 1 +0.000 0.500 0.500 1 1 1 +0.500 0.000 0.500 1 1 1 +0.500 0.500 0.000 1 1 1 +0.250 0.250 0.250 1 1 1 +0.250 0.750 0.750 1 1 1 +0.750 0.250 0.750 1 1 1 +0.750 0.750 0.250 1 1 1 \ No newline at end of file diff --git a/examples/22_rt-tddft/02_H2_velocity_gauge/INPUT b/examples/22_rt-tddft/02_H2_velocity_gauge/INPUT index b388f7dfa21..878e4f25e7e 100644 --- a/examples/22_rt-tddft/02_H2_velocity_gauge/INPUT +++ b/examples/22_rt-tddft/02_H2_velocity_gauge/INPUT @@ -20,7 +20,7 @@ smearing_method gauss #Parameters (5.MD Parameters) md_type nve -md_nstep 1000 +md_nstep 5 md_dt 0.005 md_tfirst 0 @@ -46,6 +46,9 @@ out_chg 1 out_efield 1 out_dipole 1 +out_elf 1 +out_freq_ion 0 + ### [1] Energy cutoff determines the quality of numerical quadratures in your calculations. ### So it is strongly recommended to test whether your result (such as converged SCF energies) is diff --git a/source/source_esolver/esolver_dm2rho.cpp b/source/source_esolver/esolver_dm2rho.cpp index 0903afb29a9..1548d748321 100644 --- a/source/source_esolver/esolver_dm2rho.cpp +++ b/source/source_esolver/esolver_dm2rho.cpp @@ -79,7 +79,9 @@ void ESolver_DM2rho::runner(UnitCell& ucell, const int istep) this->pelec->eferm.get_efval(is), &(ucell), 3, - 1); + 1, + PARAM.globalv.two_fermi, + false); } ModuleBase::timer::end("ESolver_DM2rho", "runner"); diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp index 14a86779edc..6a7cb86f679 100644 --- a/source/source_esolver/esolver_double_xc.cpp +++ b/source/source_esolver/esolver_double_xc.cpp @@ -158,7 +158,8 @@ void ESolver_DoubleXC::before_scf(UnitCell& ucell, const int istep) } XC_Functional::set_xc_type(PARAM.inp.deepks_out_base); - this->pelec_base->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm); + elecstate::init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, istep, + PARAM.globalv.global_out_dir, PARAM.inp, this->pelec_base); XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); // DMR should be same size with Hamiltonian(R) diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 238c093580d..1ddb636c8ab 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -175,9 +175,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) //! set direction of magnetism, used in non-collinear case elecstate::cal_ux(ucell); - //! output the initial charge density and potential - ModuleIO::write_chg_init(ucell, this->Pgrid, this->chr, this->pelec->eferm, istep, PARAM.inp); -// ModuleIO::write_pot_init(ucell, this->Pgrid, this->pelec, istep, PARAM.inp); + //! output the initial charge density + ModuleIO::write_chg_init(ucell, this->Pgrid, this->chr, this->pelec->eferm, istep, + PARAM.globalv.global_out_dir, PARAM.inp, PARAM.globalv.two_fermi); return; } diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 9d3906d6ebf..08499df2bdf 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -195,7 +195,8 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) this->dmat.dm->cal_DMR(); } // 13.2) init_scf, should be before_scf? mohan add 2025-03-10 - this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm); + elecstate::init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, + istep, PARAM.globalv.global_out_dir, PARAM.inp, this->pelec); #ifdef __MLALGO // 14) initialize DM2(R) of DeePKS, the DM2(R) is different from DM(R) diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index 033dbfa58d3..66d7717f14f 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -177,7 +177,8 @@ void ESolver_KS_LCAO_TDDFT::runner(UnitCell& ucell, const int istep) GlobalV::ofs_running, GlobalV::ofs_warning); this->exx_nao.before_scf(ucell, this->kv, this->orb_, this->p_chgmix, totstep, PARAM.inp); - this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm); + elecstate::init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, istep, + PARAM.globalv.global_out_dir, PARAM.inp, this->pelec); if (totstep <= PARAM.inp.td_tend + 1) { diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index 6714821d02f..57d0a1100be 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -165,7 +165,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) this->chr, this->locpp, this->ppcell, this->dftu, this->vsep_cell, this->stp.template get_psi_t(), this->p_hamilt, - this->pw_wfc, this->pw_rhod, PARAM.inp); + this->pw_wfc, this->pw_rhod, PARAM.globalv.global_out_dir, PARAM.inp); // setup psi (electronic wave functions) this->stp.init(this->p_hamilt); diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp index bf6f7ffde29..b346df60fb8 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/source/source_esolver/esolver_of.cpp @@ -94,7 +94,9 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp) this->init_elecstate(ucell); // calculate the total local pseudopotential in real space - this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm); // atomic_rho, v_of_rho, set_vrs + const int istep=0; + elecstate::init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep, + PARAM.globalv.global_out_dir, PARAM.inp, this->pelec); // liuyu move here 2023-10-09 // D in uspp need vloc, thus behind init_scf() @@ -188,8 +190,6 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) //! 1) call before_scf() of ESolver_FP ESolver_FP::before_scf(ucell, istep); - - if (ucell.cell_parameter_updated) { this->dV_ = ucell.omega / this->pw_rho->nxyz; @@ -232,7 +232,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) } } - this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm); + elecstate::init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep, PARAM.globalv.global_out_dir, PARAM.inp, this->pelec); Symmetry_rho::symmetrize_rho(PARAM.inp.nspin, this->chr, this->pw_rho, ucell.symm); diff --git a/source/source_esolver/lcao_others.cpp b/source/source_esolver/lcao_others.cpp index b3ad0c71499..57174c47ceb 100644 --- a/source/source_esolver/lcao_others.cpp +++ b/source/source_esolver/lcao_others.cpp @@ -172,7 +172,8 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) elecstate::cal_ux(ucell); // pelec should be initialized before these calculations - this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm); + elecstate::init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, + istep, PARAM.globalv.global_out_dir, PARAM.inp, this->pelec); // self consistent calculations for electronic ground state if (cal_type == "get_pchg") diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt index b1d4d9014f0..1bd1aff06d2 100644 --- a/source/source_estate/CMakeLists.txt +++ b/source/source_estate/CMakeLists.txt @@ -1,5 +1,6 @@ list(APPEND objects elecstate.cpp + init_scf.cpp elecstate_energy_terms.cpp elecstate_energy.cpp elecstate_exx.cpp diff --git a/source/source_estate/elecstate.cpp b/source/source_estate/elecstate.cpp index dc2ebc7fd22..5ed97e678c7 100644 --- a/source/source_estate/elecstate.cpp +++ b/source/source_estate/elecstate.cpp @@ -24,23 +24,6 @@ void ElecState::init_nelec_spin() } } -void ElecState::init_scf(const UnitCell& ucell, - const Parallel_Grid& pgrid, - const ModuleBase::ComplexMatrix& strucfac, - const bool* numeric, - ModuleSymmetry::Symmetry& symm, - const void* wfcpw) -{ - //! core correction potential. - this->charge->set_rho_core(ucell,strucfac, numeric); - - //! renormalize the charge density - this->charge->renormalize_rho(); - - //! initialize the potential - this->pot->init_pot(this->charge); -} - void ElecState::init_ks(Charge* chr_in, // pointer for class Charge const K_Vectors* klist_in, diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h index 888b3b95dc0..a4191df38d8 100644 --- a/source/source_estate/elecstate.h +++ b/source/source_estate/elecstate.h @@ -94,20 +94,8 @@ class ElecState return; } - /** - * @brief Init rho_core, init rho, renormalize rho, init pot - * - * @param ucell unit cell - * @param strucfac structure factor - * @param symm symmetry - * @param wfcpw PW basis for wave function if needed - */ - void init_scf(const UnitCell& ucell, - const Parallel_Grid& pgrid, - const ModuleBase::ComplexMatrix& strucfac, - const bool* numeric, - ModuleSymmetry::Symmetry& symm, - const void* wfcpw = nullptr); + + std::string classname = "elecstate"; int iter = 0; ///< scf iteration @@ -162,5 +150,26 @@ class ElecState bool skip_weights = false; }; +/** + * @brief Init rho_core, init rho, renormalize rho, init pot + * + * @param ucell unit cell + * @param pgrid parallel grid + * @param strucfac structure factor + * @param numeric numeric flag + * @param istep ionic step index + * @param out_dir output directory + * @param inp input parameters + * @param pelec pointer to ElecState + */ +void init_scf(const UnitCell& ucell, + const Parallel_Grid& pgrid, + const ModuleBase::ComplexMatrix& strucfac, + const bool* numeric, + const int istep, + const std::string& out_dir, + const Input_para& inp, + ElecState* pelec); + } // namespace elecstate #endif diff --git a/source/source_estate/init_scf.cpp b/source/source_estate/init_scf.cpp new file mode 100644 index 00000000000..a86c82b7b0c --- /dev/null +++ b/source/source_estate/init_scf.cpp @@ -0,0 +1,29 @@ +#include "elecstate.h" +#include "source_io/module_chgpot/write_init.h" + +namespace elecstate +{ + +void init_scf(const UnitCell& ucell, + const Parallel_Grid& pgrid, + const ModuleBase::ComplexMatrix& strucfac, + const bool* numeric, + const int istep, + const std::string& out_dir, + const Input_para& inp, + ElecState* pelec) +{ + //! core correction potential. + pelec->charge->set_rho_core(ucell, strucfac, numeric); + + //! renormalize the charge density + pelec->charge->renormalize_rho(); + + //! initialize the potential + 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); +} + +} // namespace elecstate diff --git a/source/source_estate/test/CMakeLists.txt b/source/source_estate/test/CMakeLists.txt index aa69f4e29d8..9bd76ae2486 100644 --- a/source/source_estate/test/CMakeLists.txt +++ b/source/source_estate/test/CMakeLists.txt @@ -43,7 +43,7 @@ AddTest( AddTest( TARGET MODULE_ESTATE_elecstate_base - LIBS parameter ${math_libs} base device + LIBS parameter ${math_libs} base device SOURCES elecstate_base_test.cpp ../elecstate.cpp ../elecstate_tools.cpp ../occupy.cpp ../../source_psi/psi.cpp ../../source_base/module_fft/fft_bundle.cpp ../../source_base/module_fft/fft_cpu.cpp ) @@ -58,7 +58,6 @@ AddTest( ../occupy.cpp ../module_charge/charge_mpi.cpp ../../source_psi/psi.cpp - # ../../source_psi/kernels/psi_memory_op.cpp ../../source_base/module_device/memory_op.cpp ) @@ -88,7 +87,7 @@ AddTest( AddTest( TARGET MODULE_ESTATE_charge_mixing - LIBS parameter base ${math_libs} psi device planewave_serial cell_info + LIBS parameter base ${math_libs} psi device planewave_serial cell_info SOURCES charge_mixing_test.cpp ../module_charge/charge_mixing.cpp ../module_charge/charge_mixing_dmr.cpp ../module_charge/charge_mixing_residual.cpp ../module_charge/charge_mixing_preconditioner.cpp ../module_charge/charge_mixing_rho.cpp @@ -97,7 +96,7 @@ AddTest( AddTest( TARGET MODULE_ESTATE_charge_extra - LIBS parameter ${math_libs} base device cell_info + LIBS parameter ${math_libs} base device cell_info SOURCES charge_extra_test.cpp ../module_charge/charge_extra.cpp ../../source_io/module_output/read_cube.cpp ../../source_io/module_output/write_cube.cpp ../../source_io/module_output/output.cpp ../../source_base/module_fft/fft_bundle.cpp ../../source_base/module_fft/fft_cpu.cpp ) diff --git a/source/source_estate/test/elecstate_base_test.cpp b/source/source_estate/test/elecstate_base_test.cpp index e0d63df8220..beb5a0a2e16 100644 --- a/source/source_estate/test/elecstate_base_test.cpp +++ b/source/source_estate/test/elecstate_base_test.cpp @@ -117,9 +117,6 @@ void Charge::check_rho() * - elecstate::ElecState::psiToRho() * - elecstate::ElecState::print_psi() * - elecstate::ElecState::getNewRho() - * - InitSCF: elecstate::ElecState::init_scf() - * - trivial calling to elecstate::ElecState::init_scf() - * - the production function is init charge and pot for scf calculation * - FixedWeights: elecstate::ElecState::fixed_weights() * - fix wg using external weights: ocp_kb * - CalEBand: elecstate::ElecState::cal_eband() @@ -187,7 +184,7 @@ TEST_F(ElecStateTest, Constructor) elecstate::ElecState* elecstate_new = new elecstate::ElecState(charge, rhopw, bigpw); EXPECT_EQ(elecstate_new->charge, charge); EXPECT_EQ(elecstate_new->bigpw, bigpw); - EXPECT_EQ(elecstate_new->eferm.two_efermi, PARAM.sys.two_fermi); + EXPECT_EQ(elecstate_new->eferm.two_efermi, PARAM.globalv.two_fermi); delete elecstate_new; delete bigpw; delete rhopw; @@ -253,21 +250,6 @@ TEST_F(ElecStateTest, VirtualBaseFuncs) EXPECT_NO_THROW(elecstate->getNewRho()); } -TEST_F(ElecStateTest, InitSCF) -{ - Charge* charge = new Charge; - elecstate->charge = charge; - elecstate->pot = new elecstate::Potential; - elecstate::Efermi efermi; - int istep = 0; - ModuleBase::ComplexMatrix strucfac; - elecstate->eferm = efermi; - ModuleSymmetry::Symmetry symm; - EXPECT_NO_THROW(elecstate->init_scf(ucell, pgrid, strucfac, nullptr, symm)); - // delete elecstate->pot is done in the destructor of elecstate - delete charge; -} - TEST_F(ElecStateTest, FixedWeights) { EXPECT_EQ(PARAM.input.nbands, 6); diff --git a/source/source_hsolver/test/hsolver_supplementary_mock.h b/source/source_hsolver/test/hsolver_supplementary_mock.h index 87155ecd5a7..7dc30d77233 100644 --- a/source/source_hsolver/test/hsolver_supplementary_mock.h +++ b/source/source_hsolver/test/hsolver_supplementary_mock.h @@ -19,16 +19,6 @@ void ElecState::init_nelec_spin() -void ElecState::init_scf(const UnitCell& ucell, - const Parallel_Grid& pgrid, - const ModuleBase::ComplexMatrix& strucfac, - const bool*, - ModuleSymmetry::Symmetry&, - const void*) -{ - return; -} - void ElecState::init_ks(Charge* chg_in, // pointer for class Charge const K_Vectors* klist_in, int nk_in, diff --git a/source/source_io/module_chgpot/get_pchg_lcao.cpp b/source/source_io/module_chgpot/get_pchg_lcao.cpp index 0c1b13d3b3d..e0559c4e5f1 100644 --- a/source/source_io/module_chgpot/get_pchg_lcao.cpp +++ b/source/source_io/module_chgpot/get_pchg_lcao.cpp @@ -1,6 +1,7 @@ #include "get_pchg_lcao.h" #include "source_io/module_output/cube_io.h" +#include "source_io/module_parameter/parameter.h" #include "source_estate/module_charge/symmetry_rho.h" #include "source_estate/module_dm/cal_dm_psi.h" #include "source_lcao/module_gint/gint_interface.h" @@ -95,7 +96,7 @@ void Get_pchg_lcao::begin(double** rho, double ef_spin = ef_all_spin[is]; ModuleIO::write_vdata_palgrid(pgrid, rho_save[is].data(), is, nspin, 0, - ssc.str(), ef_spin, ucell_in, precision); + ssc.str(), ef_spin, ucell_in, precision, 1, PARAM.globalv.two_fermi, false); } } } @@ -192,7 +193,7 @@ void Get_pchg_lcao::begin(double** rho, ssc.str(), ef_spin, ucell_in, - precision); + precision, 1, PARAM.globalv.two_fermi, false); } } } @@ -244,7 +245,7 @@ void Get_pchg_lcao::begin(double** rho, ssc.str(), ef_spin, ucell_in, - precision); + precision, 1, PARAM.globalv.two_fermi, false); } } } diff --git a/source/source_io/module_chgpot/get_pchg_pw.h b/source/source_io/module_chgpot/get_pchg_pw.h index 7609b30eb8a..cee654aac83 100644 --- a/source/source_io/module_chgpot/get_pchg_pw.h +++ b/source/source_io/module_chgpot/get_pchg_pw.h @@ -153,6 +153,7 @@ void get_pchg_pw(const std::vector& out_pchg, ucell, 11, 1, + PARAM.globalv.two_fermi, true); // reduce_all_pool is true } } @@ -235,7 +236,7 @@ void get_pchg_pw(const std::vector& out_pchg, std::stringstream ssc; ssc << global_out_dir << "pchgi" << ib + 1 << "s" << is + 1 << ".cube"; - ModuleIO::write_vdata_palgrid(pgrid, rho_band[is].data(), is, nspin, 0, ssc.str(), 0.0, ucell); + ModuleIO::write_vdata_palgrid(pgrid, rho_band[is].data(), is, nspin, 0, ssc.str(), 0.0, ucell, 11, 1, PARAM.globalv.two_fermi, false); } } // else if_separate_k is false } // end of ib loop over nbands diff --git a/source/source_io/module_chgpot/write_elecstat_pot.cpp b/source/source_io/module_chgpot/write_elecstat_pot.cpp index be3cd77b89e..da6c5b04df4 100644 --- a/source/source_io/module_chgpot/write_elecstat_pot.cpp +++ b/source/source_io/module_chgpot/write_elecstat_pot.cpp @@ -102,7 +102,9 @@ void write_elecstat_pot( ef_tmp, ucell, precision, - out_fermi); + out_fermi, + PARAM.globalv.two_fermi, + false); ModuleBase::timer::end("ModuleIO", "write_elecstat_pot"); return; diff --git a/source/source_io/module_chgpot/write_init.cpp b/source/source_io/module_chgpot/write_init.cpp index 8d41e617693..b7ffdc18893 100644 --- a/source/source_io/module_chgpot/write_init.cpp +++ b/source/source_io/module_chgpot/write_init.cpp @@ -1,63 +1,122 @@ +// ===================================================================== +// This module handles the output of initial charge density and potential +// to cube files in real space. It is part of the module_io package. +// +// Output files are named according to the following convention: +// - out_freq_ion = 0: chg_ini.cube/pot_ini.cube, chgs1_ini.cube/chgs2_ini.cube +// - out_freq_ion > 0: chgg{#}_ini.cube/potg{#}_ini.cube, chgs{#}g{#}_ini.cube +// - Geometry step index starts from 1 (geom_step = istep + 1) +// +// Usage: +// ModuleIO::write_chg_init(ucell, para_grid, chr, efermi, istep, out_dir, inp); +// ModuleIO::write_pot_init(ucell, para_grid, pelec, istep, out_dir, inp); +// +// Module: module_io/module_chgpot +// ===================================================================== + #include "source_io/module_chgpot/write_init.h" #include "source_io/module_output/cube_io.h" +#include "source_base/tool_quit.h" #include #include +std::string ModuleIO::gen_ini_filename( + const std::string& prefix, + const std::string& out_dir, + const int nspin, + const int is, + const int istep, + const bool include_geom_step) +{ + std::stringstream ss; + ss << out_dir << prefix; + + if (nspin == 1) + { + if (include_geom_step) + { + ss << "g" << (istep + 1) << "_ini.cube"; + } + else + { + ss << "_ini.cube"; + } + } + else if (nspin == 2 || nspin == 4) + { + ss << "s" << (is + 1); + if (include_geom_step) + { + ss << "g" << (istep + 1); + } + ss << "_ini.cube"; + } + + return ss.str(); +} + void ModuleIO::write_chg_init( const UnitCell& ucell, const Parallel_Grid ¶_grid, const Charge &chr, const elecstate::Efermi &efermi, const int istep, - const Input_para& inp) + const std::string& out_dir, + const Input_para& inp, + const bool two_fermi) { const int nspin = inp.nspin; - assert(nspin == 1 || nspin ==2 || nspin == 4); + assert(nspin == 1 || nspin == 2 || nspin == 4); + + if (istep < 0) + { + ModuleBase::WARNING_QUIT("write_chg_init", "istep must be >= 0"); + } if (inp.out_chg[0] == 2) { - for (int is = 0; is < nspin; is++) + bool should_output = (inp.out_freq_ion == 0) || + (inp.out_freq_ion > 0 && istep % inp.out_freq_ion == 0); + + if (should_output) { - std::stringstream ss; - ss << PARAM.globalv.global_out_dir << "chg"; + bool include_geom_step = (inp.out_freq_ion > 0); - if(nspin==1) - { - ss << "ini.cube"; - } - else if(nspin==2 || nspin==4) + for (int is = 0; is < nspin; is++) { - ss << "s" << is + 1 << "ini.cube"; - } + std::string filename = gen_ini_filename("chg", out_dir, nspin, is, istep, include_geom_step); + + double fermi_energy = 0.0; + if (nspin == 1 || nspin == 4) + { + fermi_energy = efermi.ef; + } + else if (nspin == 2) + { + if (is == 0) + { + fermi_energy = efermi.ef_up; + } + else if (is == 1) + { + fermi_energy = efermi.ef_dw; + } + } - // mohan add 2025-10-18 - double fermi_energy = 0.0; - if(nspin == 1 || nspin ==4) - { - fermi_energy = efermi.ef; - } - else if(nspin == 2) - { - if(is==0) - { - fermi_energy = efermi.ef_up; - } - else if(is==1) - { - fermi_energy = efermi.ef_dw; - } - } - - ModuleIO::write_vdata_palgrid(para_grid, - chr.rho[is], - is, - nspin, - istep, - ss.str(), - fermi_energy, - &(ucell), - inp.out_chg[1]); + ModuleIO::write_vdata_palgrid(para_grid, + chr.rho[is], + is, + nspin, + istep, + filename, + fermi_energy, + &(ucell), + inp.out_chg[1], + 1, + two_fermi, + false); + } } } return; @@ -69,39 +128,45 @@ void ModuleIO::write_pot_init( const Parallel_Grid ¶_grid, elecstate::ElecState *pelec, const int istep, - const Input_para& inp) + const std::string& out_dir, + const Input_para& inp, + const bool two_fermi) { - //! output total local potential of the initial charge density const int nspin = inp.nspin; - assert(nspin == 1 || nspin ==2 || nspin == 4); + assert(nspin == 1 || nspin == 2 || nspin == 4); + + if (istep < 0) + { + ModuleBase::WARNING_QUIT("write_pot_init", "istep must be >= 0"); + } if (inp.out_pot[0] == 3) { - for (int is = 0; is < nspin; is++) + bool should_output = (inp.out_freq_ion == 0) || + (inp.out_freq_ion > 0 && istep % inp.out_freq_ion == 0); + + if (should_output) { - std::stringstream ss; - ss << PARAM.globalv.global_out_dir << "pot"; + bool include_geom_step = (inp.out_freq_ion > 0); - if(nspin==1) - { - ss << "ini.cube"; - } - else if(nspin==2 || nspin==4) + for (int is = 0; is < nspin; is++) { - ss << "s" << is + 1 << "ini.cube"; - } + std::string filename = gen_ini_filename("pot", out_dir, nspin, is, istep, include_geom_step); - ModuleIO::write_vdata_palgrid(para_grid, - pelec->pot->get_eff_v(is), - is, - nspin, - istep, - ss.str(), - 0.0, // efermi - &(ucell), - 11, // precsion - 0); // out_fermi + ModuleIO::write_vdata_palgrid(para_grid, + pelec->pot->get_eff_v(is), + is, + nspin, + istep, + filename, + 0.0, + &(ucell), + inp.out_pot[1], + 0, + two_fermi, + false); + } } } - + return; } diff --git a/source/source_io/module_chgpot/write_init.h b/source/source_io/module_chgpot/write_init.h index 323cbe7af8f..303f97e42ed 100644 --- a/source/source_io/module_chgpot/write_init.h +++ b/source/source_io/module_chgpot/write_init.h @@ -1,30 +1,78 @@ +// ===================================================================== +// This module handles the output of initial charge density and potential +// to cube files in real space. It is part of the module_io package. +// +// Output files are named according to the following convention: +// - out_freq_ion = 0: chg_ini.cube/pot_ini.cube, chgs1_ini.cube/chgs2_ini.cube +// - out_freq_ion > 0: chgg{#}_ini.cube/potg{#}_ini.cube, chgs{#}g{#}_ini.cube +// - Geometry step index starts from 1 (geom_step = istep + 1) +// +// Usage: +// ModuleIO::write_chg_init(ucell, para_grid, chr, efermi, istep, out_dir, inp); +// ModuleIO::write_pot_init(ucell, para_grid, pelec, istep, out_dir, inp); +// +// Module: module_io/module_chgpot +// ===================================================================== + #ifndef WRITE_INIT_H #define WRITE_INIT_H -#include "source_io/module_parameter/input_parameter.h" // use inp -#include "source_estate/module_charge/charge.h" // use chg -#include "source_estate/fp_energy.h" // use efermi -#include "source_estate/elecstate.h" // use pelec +#include +#include "source_io/module_parameter/input_parameter.h" +#include "source_estate/module_charge/charge.h" +#include "source_estate/fp_energy.h" +#include "source_estate/elecstate.h" namespace ModuleIO { +// Generate initial data file name. +// prefix: "chg" or "pot" +std::string gen_ini_filename( + const std::string& prefix, + const std::string& out_dir, + const int nspin, + const int is, + const int istep, + const bool include_geom_step); + +// Write initial charge density to cube file in real space. +// Triggered when inp.out_chg[0] == 2. +// Output frequency is controlled by out_freq_ion: +// out_freq_ion = 0: every step output (overwrite same file) +// out_freq_ion > 0: output every out_freq_ion steps +// Output file naming convention: +// out_freq_ion = 0: chg_ini.cube (nspin=1), chgs1_ini.cube/chgs2_ini.cube (nspin=2/4) +// out_freq_ion > 0: chgg{geom_step}_ini.cube (nspin=1), chgs{spin}g{geom_step}_ini.cube (nspin=2/4) +// Note: geom_step starts from 1 (geom_step = istep + 1). void write_chg_init( const UnitCell& ucell, const Parallel_Grid ¶_grid, const Charge &chr, const elecstate::Efermi &efermi, const int istep, - const Input_para& inp); + const std::string& out_dir, + const Input_para& inp, + const bool two_fermi); +// Write initial effective potential to cube file in real space. +// Triggered when inp.out_pot[0] == 3. +// Output frequency is controlled by out_freq_ion: +// out_freq_ion = 0: every step output (overwrite same file) +// out_freq_ion > 0: output every out_freq_ion steps +// Output file naming convention: +// out_freq_ion = 0: pot_ini.cube (nspin=1), pots1_ini.cube/pots2_ini.cube (nspin=2/4) +// out_freq_ion > 0: potg{geom_step}_ini.cube (nspin=1), pots{spin}g{geom_step}_ini.cube (nspin=2/4) +// Note: geom_step starts from 1 (geom_step = istep + 1). void write_pot_init( const UnitCell& ucell, const Parallel_Grid ¶_grid, elecstate::ElecState *pelec, const int istep, - const Input_para& inp); + const std::string& out_dir, + const Input_para& inp, + const bool two_fermi); } - #endif diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index 710f4b87307..87aa062bc0c 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -1,10 +1,9 @@ #include "ctrl_output_fp.h" // use ctrl_output_fp() - -#include "../module_output/cube_io.h" // use write_vdata_palgrid -#include "../module_dipole/dipole_io.h" // use write_dipole +#include "../module_output/cube_io.h" // use write_vdata_palgrid +#include "../module_dipole/dipole_io.h" // use write_dipole #include "source_estate/module_charge/symmetry_rho.h" // use Symmetry_rho #include "source_hamilt/module_xc/xc_functional.h" // use XC_Functional -#include "source_io/module_chgpot/write_elecstat_pot.h" // use write_elecstat_pot +#include "source_io/module_chgpot/write_elecstat_pot.h" // use write_elecstat_pot #include "source_io/module_elf/write_elf.h" #ifdef USE_LIBXC @@ -49,17 +48,15 @@ void ctrl_output_fp(UnitCell& ucell, } std::string geom_block; - if (istep_in == -1) - { - // do nothing - } - else if (istep_in >= 0) + bool should_output = (PARAM.inp.out_freq_ion == 0); + if (istep_in >= 0) { geom_block = "g" + std::to_string(istep + 1); + should_output = true; } // 4) write charge density - if (PARAM.inp.out_chg[0] > 0) + if (PARAM.inp.out_chg[0] > 0 && should_output) { for (int is = 0; is < nspin; ++is) { @@ -86,7 +83,9 @@ void ctrl_output_fp(UnitCell& ucell, pelec->eferm.get_efval(is), &(ucell), PARAM.inp.out_chg[1], - 1); + 1, + PARAM.globalv.two_fermi, + false); if (XC_Functional::get_ked_flag()) { @@ -101,13 +100,17 @@ void ctrl_output_fp(UnitCell& ucell, istep, fn, pelec->eferm.get_efval(is), - &(ucell)); + &(ucell), + 11, // default precision + 1, // default out_fermi + PARAM.globalv.two_fermi, + false); } } } // 5) write potential - if (PARAM.inp.out_pot[0] == 1 || PARAM.inp.out_pot[0] == 3) + if ((PARAM.inp.out_pot[0] == 1 || PARAM.inp.out_pot[0] == 3) && should_output) { for (int is = 0; is < nspin; is++) { @@ -129,15 +132,17 @@ void ctrl_output_fp(UnitCell& ucell, pelec->pot->get_eff_v(is), is, nspin, - istep_in, + istep, fn, 0.0, // efermi &(ucell), PARAM.inp.out_pot[1], // precision - 0); // out_fermi + 0, // out_fermi + PARAM.globalv.two_fermi, + false); } } - else if (PARAM.inp.out_pot[0] == 2) + else if (PARAM.inp.out_pot[0] == 2 && should_output) { std::string fn = PARAM.globalv.global_out_dir + "potes"; fn += geom_block + ".cube"; @@ -154,11 +159,11 @@ void ctrl_output_fp(UnitCell& ucell, &(ucell), pelec->pot->get_fixed_v(), solvent, - PARAM.inp.out_pot[1]); + PARAM.inp.out_pot[1]); } // 6) write ELF - if (PARAM.inp.out_elf[0] > 0) + if (PARAM.inp.out_elf[0] > 0 && should_output) { chr.cal_elf = true; Symmetry_rho srho; @@ -168,12 +173,7 @@ void ctrl_output_fp(UnitCell& ucell, } std::string out_dir = PARAM.globalv.global_out_dir; - ModuleIO::write_elf( -#ifdef __MPI - pw_big->bz, - pw_big->nbz, -#endif - out_dir, + ModuleIO::write_elf(out_dir, istep, nspin, chr.rho, @@ -181,12 +181,14 @@ void ctrl_output_fp(UnitCell& ucell, pw_rhod, para_grid, &(ucell), - PARAM.inp.out_elf[1]); + PARAM.inp.out_elf[1], + geom_block, + PARAM.globalv.two_fermi); } #ifdef USE_LIBXC // 7) write xc(r) - if (PARAM.inp.out_xc_r[0] >= 0) + if (PARAM.inp.out_xc_r[0] >= 0 && should_output) { ModuleIO::write_libxc_r(PARAM.inp.out_xc_r[0], XC_Functional::get_func_id(), @@ -200,7 +202,7 @@ void ctrl_output_fp(UnitCell& ucell, #endif // 8) write dipole moment - if (PARAM.inp.out_dipole == 1) + if (PARAM.inp.out_dipole == 1 && should_output) { for (int is = 0; is < nspin; ++is) { diff --git a/source/source_io/module_dos/cal_ldos.cpp b/source/source_io/module_dos/cal_ldos.cpp index 4463c4cf036..1e92e2e4e4b 100644 --- a/source/source_io/module_dos/cal_ldos.cpp +++ b/source/source_io/module_dos/cal_ldos.cpp @@ -90,7 +90,9 @@ void Cal_ldos::cal_ldos_lcao( 0, &ucell, precision, - 0); + 0, + PARAM.globalv.two_fermi, + false); // free memory delete[] ldos; @@ -162,7 +164,7 @@ void stm_mode_pw(const elecstate::ElecStatePW>* pelec, << ".cube"; const int precision = PARAM.inp.out_ldos[1]; - ModuleIO::write_vdata_palgrid(pgrid, ldos.data(), 0, PARAM.inp.nspin, 0, fn.str(), 0, &ucell, precision, 0); + ModuleIO::write_vdata_palgrid(pgrid, ldos.data(), 0, PARAM.inp.nspin, 0, fn.str(), 0, &ucell, precision, 0, PARAM.globalv.two_fermi, false); } } diff --git a/source/source_io/module_elf/write_elf.cpp b/source/source_io/module_elf/write_elf.cpp index 347a119099b..ea6e9b88050 100644 --- a/source/source_io/module_elf/write_elf.cpp +++ b/source/source_io/module_elf/write_elf.cpp @@ -1,74 +1,100 @@ #include "write_elf.h" #include "source_io/module_output/cube_io.h" +#ifdef _OPENMP +#include +#endif namespace ModuleIO { void write_elf( -#ifdef __MPI - const int& bz, - const int& nbz, -#endif const std::string& out_dir, const int& istep_in, const int& nspin, const double* const* rho, const double* const* tau, - ModulePW::PW_Basis* rho_basis, + ModulePW::PW_Basis* const rho_basis, const Parallel_Grid& pgrid, const UnitCell* ucell_, - const int& precision) + const int& precision, + const std::string& geom_block, + const bool two_fermi) { - // For nspin = 4, we only calculate the total ELF using the rho_total and tau_total, - // containing in the first channel of rho and tau. - // What's more, we have not introduced the U(1) and SU(2) gauge invariance corrections - // proposed by Desmarais J K, Vignale G, Bencheikh K, et al. Physical Review Letters, 2024, 133(13): 136401, - // where the current density is also included in the ELF calculation. + ModuleBase::timer::start("ModuleIO", "write_elf"); + // For nspin = 4, we only calculate the total ELF using the + // rho_total and tau_total, containing in the first channel of + // rho and tau. + // What's more, we have not introduced the U(1) and SU(2) gauge + // invariance corrections proposed by Desmarais J K, Vignale G, + // Bencheikh K, et al. Physical Review Letters, 2024, + // 133(13): 136401, where the current density is also included + // in the ELF calculation. + + const int nspin_eff = (nspin == 4) ? 1 : nspin; + + const int nrxx = rho_basis->nrxx; + const int npw = rho_basis->npw; - int nspin_eff = (nspin == 4) ? 1 : nspin; + assert(nrxx>0); + assert(npw>0); - std::vector> elf(nspin_eff, std::vector(rho_basis->nrxx, 0.)); + std::vector> elf(nspin_eff, std::vector(nrxx, 0.)); // 1) calculate the kinetic energy density of vW KEDF - std::vector> tau_vw(nspin_eff, std::vector(rho_basis->nrxx, 0.)); - std::vector phi(rho_basis->nrxx, 0.); // phi = sqrt(rho) + std::vector> tau_vw(nspin_eff, std::vector(nrxx, 0.)); + std::vector phi(nrxx, 0.); for (int is = 0; is < nspin_eff; ++is) { - for (int ir = 0; ir < rho_basis->nrxx; ++ir) +#pragma omp parallel for schedule(static) \ + default(none) firstprivate(nrxx) shared(phi, rho, rho_basis, is) + for (int ir = 0; ir < nrxx; ++ir) { phi[ir] = std::sqrt(std::abs(rho[is][ir])); } - - std::vector> gradient_phi(3, std::vector(rho_basis->nrxx, 0.)); - std::vector> recip_phi(rho_basis->npw, 0.0); - std::vector> recip_gradient_phi(rho_basis->npw, 0.0); - + + std::vector> gradient_phi(3, std::vector(nrxx, 0.)); + std::vector> recip_phi(npw, 0.0); + std::vector> recip_gradient_phi(npw, 0.0); + rho_basis->real2recip(phi.data(), recip_phi.data()); - + std::complex img(0.0, 1.0); for (int j = 0; j < 3; ++j) { - for (int ip = 0; ip < rho_basis->npw; ++ip) +#pragma omp parallel for schedule(static) \ + default(none) firstprivate(img, npw) \ + shared(recip_gradient_phi, rho_basis, recip_phi, j) + for (int ip = 0; ip < npw; ++ip) { - recip_gradient_phi[ip] = img * rho_basis->gcar[ip][j] * recip_phi[ip] * rho_basis->tpiba; + recip_gradient_phi[ip] + = img * rho_basis->gcar[ip][j] + * recip_phi[ip] * rho_basis->tpiba; } rho_basis->recip2real(recip_gradient_phi.data(), gradient_phi[j].data()); - for (int ir = 0; ir < rho_basis->nrxx; ++ir) +#pragma omp parallel for schedule(static) \ + default(none) firstprivate(nrxx) \ + shared(tau_vw, gradient_phi, is, j, rho_basis) + for (int ir = 0; ir < nrxx; ++ir) { - tau_vw[is][ir] += gradient_phi[j][ir] * gradient_phi[j][ir] / 2. * 2.; // convert Ha to Ry. + tau_vw[is][ir] += gradient_phi[j][ir] + * gradient_phi[j][ir] / 2. * 2.; } } } // 2) calculate the kinetic energy density of TF KEDF - std::vector> tau_TF(nspin_eff, std::vector(rho_basis->nrxx, 0.)); + std::vector> tau_TF(nspin_eff, std::vector(nrxx, 0.)); + const double c_tf = 3.0 / 10.0 * std::pow(3 * std::pow(M_PI, 2.0), 2.0 / 3.0) - * 2.0; // 10/3*(3*pi^2)^{2/3}, multiply by 2 to convert unit from Hartree to Ry, finally in Ry*Bohr^(-2) + * 2.0; // convert unit from Hartree to Ry if (nspin == 1 || nspin == 4) { - for (int ir = 0; ir < rho_basis->nrxx; ++ir) +#pragma omp parallel for schedule(static) \ + default(none) firstprivate(c_tf, nrxx) \ + shared(rho, tau_TF, rho_basis) + for (int ir = 0; ir < nrxx; ++ir) { if (rho[0][ir] > 0.0) { @@ -82,14 +108,19 @@ void write_elf( } else if (nspin == 2) { - // the spin-scaling law: tau_TF[rho_up, rho_dn] = 1/2 * (tau_TF[2*rho_up] + tau_TF[2*rho_dn]) + // spin-scaling: tau_TF[rho_up,rho_dn] + // = 1/2 * (tau_TF[2*rho_up] + tau_TF[2*rho_dn]) for (int is = 0; is < nspin; ++is) { - for (int ir = 0; ir < rho_basis->nrxx; ++ir) +#pragma omp parallel for schedule(static) \ + default(none) firstprivate(c_tf, nrxx) \ + shared(rho, tau_TF, is, rho_basis) + for (int ir = 0; ir < nrxx; ++ir) { if (rho[is][ir] > 0.0) { - tau_TF[is][ir] = 0.5 * c_tf * std::pow(2.0 * rho[is][ir], 5.0 / 3.0); + tau_TF[is][ir] = 0.5 * c_tf + * std::pow(2.0 * rho[is][ir], 5.0 / 3.0); } else { @@ -99,15 +130,20 @@ void write_elf( } } - // 3) calculate the enhancement factor F = (tau_KS - tau_vw) / tau_TF, and then ELF = 1 / (1 + F^2) - double eps = 1.0e-5; // suppress the numerical instability in LCAO (Ref: Acta Phys. -Chim. Sin. 2011, 27(12), 2786-2792. doi: 10.3866/PKU.WHXB20112786) + // 3) calculate the enhancement factor F = (tau_KS - tau_vw) / + // tau_TF, and then ELF = 1 / (1 + F^2) + const double eps = 1.0e-5; for (int is = 0; is < nspin_eff; ++is) { - for (int ir = 0; ir < rho_basis->nrxx; ++ir) +#pragma omp parallel for schedule(static) \ + default(none) firstprivate(eps, nrxx) \ + shared(elf, tau, tau_vw, tau_TF, is, rho_basis) + for (int ir = 0; ir < nrxx; ++ir) { if (tau_TF[is][ir] > 1.0e-12) { - elf[is][ir] = (tau[is][ir] - tau_vw[is][ir] + eps) / tau_TF[is][ir]; + elf[is][ir] = (tau[is][ir] - tau_vw[is][ir] + eps) + / tau_TF[is][ir]; elf[is][ir] = 1. / (1. + elf[is][ir] * elf[is][ir]); } else @@ -118,12 +154,12 @@ void write_elf( } // 4) output the ELF = 1 / (1 + F^2) to cube file - double ef_tmp = 0.0; - int out_fermi = 0; + const double ef_tmp = 0.0; + const int out_fermi = 0; if (nspin == 1 || nspin == 4) { - std::string fn = out_dir + "/elf.cube"; + std::string fn = out_dir + "elftot" + geom_block + ".cube"; int is = -1; ModuleIO::write_vdata_palgrid(pgrid, @@ -135,17 +171,18 @@ void write_elf( ef_tmp, ucell_, precision, - out_fermi); + out_fermi, + two_fermi, + false); } else if (nspin == 2) { for (int is = 0; is < nspin; ++is) { - std::string fn_temp = out_dir + "/elf"; - - fn_temp += std::to_string(is + 1) + ".cube"; + std::string fn_temp = out_dir + "elf" + "s" + + std::to_string(is + 1) + geom_block + ".cube"; - int ispin = is + 1; + const int ispin = is + 1; ModuleIO::write_vdata_palgrid(pgrid, elf[is].data(), @@ -156,16 +193,31 @@ void write_elf( ef_tmp, ucell_, precision, - out_fermi); + out_fermi, + two_fermi, + false); } - std::vector elf_tot(rho_basis->nrxx, 0.0); - for (int ir = 0; ir < rho_basis->nrxx; ++ir) + std::vector elf_tot(nrxx, 0.0); +#pragma omp parallel for schedule(static) \ + default(none) firstprivate(eps, nrxx) \ + shared(elf_tot, tau, tau_vw, tau_TF, rho_basis) + for (int ir = 0; ir < nrxx; ++ir) { - elf_tot[ir] = (tau[0][ir] + tau[1][ir] - tau_vw[0][ir] - tau_vw[1][ir]) / (tau_TF[0][ir] + tau_TF[1][ir]); - elf_tot[ir] = 1. / (1. + elf_tot[ir] * elf_tot[ir]); + if (tau_TF[0][ir] + tau_TF[1][ir] > 1.0e-12) + { + elf_tot[ir] = (tau[0][ir] + tau[1][ir] + - tau_vw[0][ir] - tau_vw[1][ir] + eps) + / (tau_TF[0][ir] + tau_TF[1][ir]); + elf_tot[ir] = 1. + / (1. + elf_tot[ir] * elf_tot[ir]); + } + else + { + elf_tot[ir] = 0.0; + } } - std::string fn = out_dir + "/elf.cube"; + std::string fn = out_dir + "elftot" + geom_block + ".cube"; int is = -1; ModuleIO::write_vdata_palgrid(pgrid, @@ -177,7 +229,10 @@ void write_elf( ef_tmp, ucell_, precision, - out_fermi); + out_fermi, + two_fermi, + false); } -} -} + ModuleBase::timer::end("ModuleIO", "write_elf"); +} // end write_elf +} // end namespace ModuleIO diff --git a/source/source_io/module_elf/write_elf.h b/source/source_io/module_elf/write_elf.h index 23a7e7c379e..ea5f1b788e7 100644 --- a/source/source_io/module_elf/write_elf.h +++ b/source/source_io/module_elf/write_elf.h @@ -8,10 +8,6 @@ namespace ModuleIO { void write_elf( -#ifdef __MPI - const int& bz, - const int& nbz, -#endif const std::string& out_dir, const int& istep_in, const int& nspin, @@ -20,7 +16,9 @@ void write_elf( ModulePW::PW_Basis* rho_basis, const Parallel_Grid& pgrid, const UnitCell* ucell_, - const int& precision); + const int& precision, + const std::string& geom_block, + const bool two_fermi); } #endif diff --git a/source/source_io/module_output/cube_io.h b/source/source_io/module_output/cube_io.h index ce4fc69f82f..f97a639b1d0 100644 --- a/source/source_io/module_output/cube_io.h +++ b/source/source_io/module_output/cube_io.h @@ -24,9 +24,10 @@ void write_vdata_palgrid(const Parallel_Grid& pgrid, const std::string& fn, const double ef, const UnitCell* const ucell, - const int precision = 11, - const int out_fermi = 1, - const bool reduce_all_pool = false); // only reduce in the main pool as default + const int precision, + const int out_fermi, + const bool two_fermi, + const bool reduce_all_pool); /// read the full data from a cube file bool read_cube(const std::string& file, diff --git a/source/source_io/module_output/write_cube.cpp b/source/source_io/module_output/write_cube.cpp index 39c5454d4b4..e7ae67a7631 100644 --- a/source/source_io/module_output/write_cube.cpp +++ b/source/source_io/module_output/write_cube.cpp @@ -2,7 +2,6 @@ #include "source_base/parallel_comm.h" #include "source_pw/module_pwdft/parallel_grid.h" #include "source_io/module_output/cube_io.h" -#include "source_io/module_parameter/parameter.h" #include @@ -22,6 +21,7 @@ void ModuleIO::write_vdata_palgrid(const Parallel_Grid& pgrid, const UnitCell* const ucell, const int precision, const int out_fermi, + const bool two_fermi, const bool reduce_all_pool) { ModuleBase::TITLE("ModuleIO", "write_vdata_palgrid"); @@ -73,7 +73,7 @@ void ModuleIO::write_vdata_palgrid(const Parallel_Grid& pgrid, ss << nspin << " # number of spin directions "; if (out_fermi == 1) { - if (PARAM.globalv.two_fermi) + if (two_fermi) { if (is == 0) { diff --git a/source/source_io/module_parameter/read_input_item_output.cpp b/source/source_io/module_parameter/read_input_item_output.cpp index 4881cb7d48a..8bc0c895499 100644 --- a/source/source_io/module_parameter/read_input_item_output.cpp +++ b/source/source_io/module_parameter/read_input_item_output.cpp @@ -80,7 +80,18 @@ void ReadInput::item_output() - nspin = 1: `tau.cube`; - nspin = 2: `taus1.cube`, and `taus2.cube`; - nspin = 4: `taus1.cube`, `taus2.cube`, `taus3.cube`, and `taus4.cube`; - - 2: On top of 1, also output the initial charge density files with a suffix name as '_ini', such as `taus1_ini.cube`, etc. + - 2: On top of 1, also output the initial charge density files. The files are named as: + - out_freq_ion = 0: + - nspin = 1: `chg_ini.cube`; + - nspin = 2: `chgs1_ini.cube` and `chgs2_ini.cube`; + - nspin = 4: `chgs1_ini.cube`, `chgs2_ini.cube`, `chgs3_ini.cube`, and `chgs4_ini.cube`; + - output at every step (overwrite same file) + - out_freq_ion > 0: + - nspin = 1: `chgg{geom_step}_ini.cube` (e.g., `chgg1_ini.cube`); + - nspin = 2: `chgs1g{geom_step}_ini.cube` and `chgs2g{geom_step}_ini.cube`; + - nspin = 4: `chgs1g{geom_step}_ini.cube`, `chgs2g{geom_step}_ini.cube`, `chgs3g{geom_step}_ini.cube`, and `chgs4g{geom_step}_ini.cube`. + - output every out_freq_ion steps + Here, {geom_step} denotes the geometry step index, starting from 1 (geom_step = istep + 1). - -1: Disable the charge density auto-back-up file `{suffix}-CHARGE-DENSITY.restart`, useful for large systems. The second integer controls the precision of the charge density output. If not given, `3` is used as default. For restarting from this file and other high-precision calculations, `10` is recommended. @@ -120,9 +131,17 @@ In molecular dynamics simulations, the output frequency is controlled by out_fre * nspin = 4: pots1.cube, pots2.cube, pots3.cube, and pots4.cube * 2: Output the electrostatic potential on real space grids into OUT.{suffix}/pot_es.cube. The Python script named tools/average_pot/aveElecStatPot.py can be used to calculate the average electrostatic potential along the z-axis and outputs it into ElecStaticPot_AVE. Please note that the total local potential refers to the local component of the self-consistent potential, excluding the non-local pseudopotential. The distinction between the local potential and the electrostatic potential is as follows: local potential = electrostatic potential + XC potential. * 3: Apart from 1, also output the total local potential of the initial charge density. The files are named as: - * nspin = 1: pots1_ini.cube; - * nspin = 2: pots1_ini.cube and pots2_ini.cube; - * nspin = 4: pots1_ini.cube, pots2_ini.cube, pots3_ini.cube, and pots4_ini.cube + * out_freq_ion = 0: + * nspin = 1: `pot_ini.cube`; + * nspin = 2: `pots1_ini.cube` and `pots2_ini.cube`; + * nspin = 4: `pots1_ini.cube`, `pots2_ini.cube`, `pots3_ini.cube`, and `pots4_ini.cube`; + * output at every step (overwrite same file) + * out_freq_ion > 0: + * nspin = 1: `potg{geom_step}_ini.cube` (e.g., `potg1_ini.cube`); + * nspin = 2: `pots1g{geom_step}_ini.cube` and `pots2g{geom_step}_ini.cube`; + * nspin = 4: `pots1g{geom_step}_ini.cube`, `pots2g{geom_step}_ini.cube`, `pots3g{geom_step}_ini.cube`, and `pots4g{geom_step}_ini.cube`. + * output every out_freq_ion steps + Here, {geom_step} denotes the geometry step index, starting from 1 (geom_step = istep + 1). The optional second integer controls the output precision. If not provided, the default precision is 8. @@ -966,12 +985,14 @@ If EXX(exact exchange) is calculated (i.e. dft_fuctional==hse/hf/pbe0/scan0 or r item.type = R"(Integer \[Integer\](optional))"; item.description = R"(Whether to output the electron localization function (ELF) in the folder `OUT.${suffix}`. The files are named as * nspin = 1: - * elf.cube: ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i}|^2} - \frac{|\nabla\rho|^2}{8\rho}}{\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}}$; + * elftot.cube: ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i}|^2} - \frac{|\nabla\rho|^2}{8\rho}}{\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}}$; * nspin = 2: - * elf1.cube, elf2.cube: ${\rm{ELF}}_\sigma = \frac{1}{1+\chi_\sigma^2}$, $\chi_\sigma = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i,\sigma}|^2} - \frac{|\nabla\rho_\sigma|^2}{8\rho_\sigma}}{\frac{3}{10}(6\pi^2)^{2/3}\rho_\sigma^{5/3}}$; - * elf.cube: ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i,\sigma}{f_i |\nabla\psi_{i,\sigma}|^2} - \sum_{\sigma}{\frac{|\nabla\rho_\sigma|^2}{8\rho_\sigma}}}{\sum_{\sigma}{\frac{3}{10}(6\pi^2)^{2/3}\rho_\sigma^{5/3}}}$; + * elfs1.cube, elfs2.cube: ${\rm{ELF}}_\sigma = \frac{1}{1+\chi_\sigma^2}$, $\chi_\sigma = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i,\sigma}|^2} - \frac{|\nabla\rho_\sigma|^2}{8\rho_\sigma}}{\frac{3}{10}(6\pi^2)^{2/3}\rho_\sigma^{5/3}}$; + * elftot.cube: ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i,\sigma}{f_i |\nabla\psi_{i,\sigma}|^2} - \sum_{\sigma}{\frac{|\nabla\rho_\sigma|^2}{8\rho_\sigma}}}{\sum_{\sigma}{\frac{3}{10}(6\pi^2)^{2/3}\rho_\sigma^{5/3}}}$; * nspin = 4 (noncollinear): - * elf.cube: ELF for total charge density, ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i}|^2} - \frac{|\nabla\rho|^2}{8\rho}}{\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}}$ + * elftot.cube: ELF for total charge density, ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i}|^2} - \frac{|\nabla\rho|^2}{8\rho}}{\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}}$ + +When `out_freq_ion > 0`, a geometry step suffix `g{#}` is appended to the file names (e.g., `elftotg1.cube`, `elfs1g1.cube`). The second integer controls the precision of the kinetic energy density output, if not given, will use 3 as default. For purpose restarting from this file and other high-precision involved calculation, recommend to use 10. @@ -991,9 +1012,9 @@ In molecular dynamics calculations, the output frequency is controlled by out_fr } }; item.check_value = [](const Input_Item& item, const Parameter& para) { - if (para.input.out_elf[0] > 0 && para.input.esolver_type != "ksdft" && para.input.esolver_type != "ofdft") + if (para.input.out_elf[0] > 0 && para.input.esolver_type != "ksdft" && para.input.esolver_type != "ofdft" && para.input.esolver_type != "tddft") { - ModuleBase::WARNING_QUIT("ReadInput", "ELF is only aviailable for ksdft and ofdft"); + ModuleBase::WARNING_QUIT("ReadInput", "ELF is only available for ksdft, ofdft and tddft"); } }; sync_intvec(input.out_elf, 2, 0); diff --git a/source/source_io/module_wf/get_wf_lcao.cpp b/source/source_io/module_wf/get_wf_lcao.cpp index cadf8d0edc0..cd54b4cbdd5 100644 --- a/source/source_io/module_wf/get_wf_lcao.cpp +++ b/source/source_io/module_wf/get_wf_lcao.cpp @@ -83,7 +83,11 @@ void Get_wf_lcao::begin(const UnitCell& ucell, 0, ss_out.str(), ef_tmp, - &(ucell)); + &(ucell), + 11, // default precision + 1, // default out_fermi + PARAM.globalv.two_fermi, + false); } } } @@ -126,12 +130,12 @@ void Get_wf_lcao::begin(const UnitCell& ucell, // Output real part std::stringstream ss_real; ss_real << global_out_dir << "wfi" << ib + 1 << "s" << is + 1 << "re.cube"; - ModuleIO::write_vdata_palgrid(pgrid, wfc_real.data(), is, nspin, 0, ss_real.str(), ef_tmp, &(ucell)); + ModuleIO::write_vdata_palgrid(pgrid, wfc_real.data(), is, nspin, 0, ss_real.str(), ef_tmp, &(ucell), 11, 1, PARAM.globalv.two_fermi, false); // Output imaginary part std::stringstream ss_imag; ss_imag << global_out_dir << "wfi" << ib + 1 << "s" << is + 1 << "im.cube"; - ModuleIO::write_vdata_palgrid(pgrid, wfc_imag.data(), is, nspin, 0, ss_imag.str(), ef_tmp, &(ucell)); + ModuleIO::write_vdata_palgrid(pgrid, wfc_imag.data(), is, nspin, 0, ss_imag.str(), ef_tmp, &(ucell), 11, 1, PARAM.globalv.two_fermi, false); } } } @@ -245,7 +249,9 @@ void Get_wf_lcao::begin(const UnitCell& ucell, ef_tmp, &(ucell), 3, - 1); + 1, + PARAM.globalv.two_fermi, + false); // if (out_wfc_pw || out_wf_r) psi_g.fix_k(ik); @@ -316,12 +322,12 @@ void Get_wf_lcao::begin(const UnitCell& ucell, ss_real << global_out_dir << "wfi" << ib + 1 << "s" << ispin + 1 << "k" << ik0 + 1 << "re.cube"; const double ef_tmp = this->pes_->eferm.get_efval(ispin); - ModuleIO::write_vdata_palgrid(pgrid, wfc_real.data(), ispin, nspin, 0, ss_real.str(), ef_tmp, &(ucell)); + ModuleIO::write_vdata_palgrid(pgrid, wfc_real.data(), ispin, nspin, 0, ss_real.str(), ef_tmp, &(ucell), 11, 1, PARAM.globalv.two_fermi, false); // Output imaginary part std::stringstream ss_imag; ss_imag << global_out_dir << "wfi" << ib + 1 << "s" << ispin + 1 << "k" << ik0 + 1 << "im.cube"; - ModuleIO::write_vdata_palgrid(pgrid, wfc_imag.data(), ispin, nspin, 0, ss_imag.str(), ef_tmp, &(ucell)); + ModuleIO::write_vdata_palgrid(pgrid, wfc_imag.data(), ispin, nspin, 0, ss_imag.str(), ef_tmp, &(ucell), 11, 1, PARAM.globalv.two_fermi, false); } } } diff --git a/source/source_io/module_wf/get_wf_pw.h b/source/source_io/module_wf/get_wf_pw.h index 9262ad181b6..0b9461ea68b 100644 --- a/source/source_io/module_wf/get_wf_pw.h +++ b/source/source_io/module_wf/get_wf_pw.h @@ -160,6 +160,7 @@ void get_wf_pw(const std::vector& out_wfc_norm, ucell, 11, 1, + PARAM.globalv.two_fermi, true); // reduce_all_pool is true } } @@ -250,6 +251,7 @@ void get_wf_pw(const std::vector& out_wfc_norm, ucell, 11, 1, + PARAM.globalv.two_fermi, true); // reduce_all_pool is true std::stringstream ss_imag; @@ -265,6 +267,7 @@ void get_wf_pw(const std::vector& out_wfc_norm, ucell, 11, 1, + PARAM.globalv.two_fermi, true); // reduce_all_pool is true } } diff --git a/source/source_io/test_serial/rho_io_test.cpp b/source/source_io/test_serial/rho_io_test.cpp index 1a96ecc5627..20f5a843e68 100644 --- a/source/source_io/test_serial/rho_io_test.cpp +++ b/source/source_io/test_serial/rho_io_test.cpp @@ -118,7 +118,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); + ModuleIO::write_vdata_palgrid(pgrid, rho[0], 0, nspin, 0, "test_write_vdata_palgrid.cube", 0.461002, ucell, 11, 1, false, false); 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_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp index 1073774b38c..6a4697ebaf9 100644 --- a/source/source_pw/module_pwdft/setup_pot.cpp +++ b/source/source_pw/module_pwdft/setup_pot.cpp @@ -23,6 +23,7 @@ void pw::setup_pot(const int istep, hamilt::HamiltBase* p_hamilt, // hamiltonian ModulePW::PW_Basis_K *pw_wfc, // pw for wfc const ModulePW::PW_Basis *pw_rhod, // pw for rhod + const std::string& out_dir, const Input_para& inp) // input parameters { ModuleBase::TITLE("pw", "setup_pot"); @@ -39,8 +40,8 @@ void pw::setup_pot(const int istep, //---------------------------------------------------------- //! 1) Renew local pseudopotential //---------------------------------------------------------- - pelec->init_scf(ucell, para_grid, sf.strucFac, - locpp.numeric, ucell.symm, (void*)pw_wfc); + elecstate::init_scf(ucell, para_grid, sf.strucFac, + locpp.numeric, istep, out_dir, inp, pelec); //---------------------------------------------------------- //! 2) Symmetrize the charge density (rho) @@ -141,6 +142,7 @@ template void pw::setup_pot, base_device::DEVICE_CPU>( hamilt::HamiltBase* p_hamilt, // hamiltonian ModulePW::PW_Basis_K *pw_wfc, // pw for wfc const ModulePW::PW_Basis *pw_rhod, // pw for rhod + const std::string& out_dir, const Input_para& inp); // input parameters @@ -160,6 +162,7 @@ template void pw::setup_pot, base_device::DEVICE_CPU>( hamilt::HamiltBase* p_hamilt, // hamiltonian ModulePW::PW_Basis_K *pw_wfc, // pw for wfc const ModulePW::PW_Basis *pw_rhod, // pw for rhod + const std::string& out_dir, const Input_para& inp); // input parameters #if ((defined __CUDA) || (defined __ROCM)) @@ -180,6 +183,7 @@ template void pw::setup_pot, base_device::DEVICE_GPU>( hamilt::HamiltBase* p_hamilt, // hamiltonian ModulePW::PW_Basis_K *pw_wfc, // pw for wfc const ModulePW::PW_Basis *pw_rhod, // pw for rhod + const std::string& out_dir, const Input_para& inp); // input parameters template void pw::setup_pot, base_device::DEVICE_GPU>( @@ -198,6 +202,7 @@ template void pw::setup_pot, base_device::DEVICE_GPU>( hamilt::HamiltBase* p_hamilt, // hamiltonian ModulePW::PW_Basis_K *pw_wfc, // pw for wfc const ModulePW::PW_Basis *pw_rhod, // pw for rhod + const std::string& out_dir, const Input_para& inp); // input parameters #endif diff --git a/source/source_pw/module_pwdft/setup_pot.h b/source/source_pw/module_pwdft/setup_pot.h index a292579a6c2..ecb89bc213a 100644 --- a/source/source_pw/module_pwdft/setup_pot.h +++ b/source/source_pw/module_pwdft/setup_pot.h @@ -30,6 +30,7 @@ void setup_pot(const int istep, hamilt::HamiltBase* p_hamilt, // hamiltonian ModulePW::PW_Basis_K *pw_wfc, // pw for wfc const ModulePW::PW_Basis *pw_rhod, // pw for rhod + const std::string& out_dir, const Input_para& inp); // input parameters } diff --git a/tests/01_PW/nscf_out_pot/INPUT b/tests/01_PW/nscf_out_pot/INPUT index 63fd05c8cc4..3afecb93938 100644 --- a/tests/01_PW/nscf_out_pot/INPUT +++ b/tests/01_PW/nscf_out_pot/INPUT @@ -5,7 +5,7 @@ calculation nscf nbands 8 symmetry 1 -out_pot 1 +out_pot 1 3 #Parameters (2.Iteration) ecutwfc 5 diff --git a/tests/01_PW/nscf_out_pot/pot.cube.ref b/tests/01_PW/nscf_out_pot/pot.cube.ref index b7de1a905ef..9490ba892b8 100644 --- a/tests/01_PW/nscf_out_pot/pot.cube.ref +++ b/tests/01_PW/nscf_out_pot/pot.cube.ref @@ -1,4 +1,4 @@ -Ionic_Step 0 Cubefile created from ABACUS. Inner loop is z, followed by y and x +Ionic_Step 1 Cubefile created from ABACUS. Inner loop is z, followed by y and x 1 # number of spin directions 1 0.0 0.0 0.0 9 0.000000 0.495556 0.495556 diff --git a/tests/01_PW/scf_out_elf/refelf.cube b/tests/01_PW/scf_out_elf/refelftot.cube similarity index 100% rename from tests/01_PW/scf_out_elf/refelf.cube rename to tests/01_PW/scf_out_elf/refelftot.cube diff --git a/tests/03_NAO_multik/nscf_out_pot1/pot.cube.ref b/tests/03_NAO_multik/nscf_out_pot1/pot.cube.ref index 9b45c89cf68..6bff38268fd 100644 --- a/tests/03_NAO_multik/nscf_out_pot1/pot.cube.ref +++ b/tests/03_NAO_multik/nscf_out_pot1/pot.cube.ref @@ -1,4 +1,4 @@ -Ionic_Step 0 Cubefile created from ABACUS. Inner loop is z, followed by y and x +Ionic_Step 1 Cubefile created from ABACUS. Inner loop is z, followed by y and x 1 # number of spin directions 2 0.0 0.0 0.0 12 0.000000 0.425000 0.425000 diff --git a/tests/03_NAO_multik/scf_out_chg_pot1/pot.cube.ref b/tests/03_NAO_multik/scf_out_chg_pot1/pot.cube.ref index 2b7ef788960..6dddbcfc4cd 100644 --- a/tests/03_NAO_multik/scf_out_chg_pot1/pot.cube.ref +++ b/tests/03_NAO_multik/scf_out_chg_pot1/pot.cube.ref @@ -1,4 +1,4 @@ -Ionic_Step 0 Cubefile created from ABACUS. Inner loop is z, followed by y and x +Ionic_Step 1 Cubefile created from ABACUS. Inner loop is z, followed by y and x 1 # number of spin directions 2 0.0 0.0 0.0 8 0.662500 0.000000 0.000000 diff --git a/tests/03_NAO_multik/scf_out_elf/refelf.cube b/tests/03_NAO_multik/scf_out_elf/refelftot.cube similarity index 100% rename from tests/03_NAO_multik/scf_out_elf/refelf.cube rename to tests/03_NAO_multik/scf_out_elf/refelftot.cube diff --git a/tests/07_OFDFT/24_OF_out_elf/refelf.cube b/tests/07_OFDFT/24_OF_out_elf/refelftot.cube similarity index 100% rename from tests/07_OFDFT/24_OF_out_elf/refelf.cube rename to tests/07_OFDFT/24_OF_out_elf/refelftot.cube diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index 859d35309fc..c6070a8fc44 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -219,8 +219,8 @@ fi # echo $out_elf #------------------------------- if ! test -z "$out_elf" && [ $out_elf == 1 ]; then - elf1ref=refelf.cube - elf1cal=OUT.autotest/elf.cube + elf1ref=refelftot.cube + elf1cal=OUT.autotest/elftot.cube python3 $COMPARE_SCRIPT $elf1ref $elf1cal 3 echo "ComparePot1_pass $?" >>$1 fi