From d14c9dc0002f30b864325f482aa911f285e60427 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 15:00:54 +0800 Subject: [PATCH 01/31] update 18_md examples --- examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT | 0 examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_0 | 0 examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_1 | 0 examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_2 | 0 examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_3 | 0 examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_4 | 0 examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_5 | 0 examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_6 | 0 examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_7 | 0 examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/KPT | 0 examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/README | 0 examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/STRU | 0 examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/run.sh | 0 13 files changed, 0 insertions(+), 0 deletions(-) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT (100%) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_0 (100%) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_1 (100%) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_2 (100%) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_3 (100%) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_4 (100%) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_5 (100%) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_6 (100%) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/INPUT_7 (100%) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/KPT (100%) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/README (100%) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/STRU (100%) rename examples/18_md/{01_lcao_gamma_Si8 => 01_LCAO_NVE_Si8}/run.sh (100%) 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_0 b/examples/18_md/01_LCAO_NVE_Si8/INPUT_0 similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/INPUT_0 rename to examples/18_md/01_LCAO_NVE_Si8/INPUT_0 diff --git a/examples/18_md/01_lcao_gamma_Si8/INPUT_1 b/examples/18_md/01_LCAO_NVE_Si8/INPUT_1 similarity index 100% rename from examples/18_md/01_lcao_gamma_Si8/INPUT_1 rename to examples/18_md/01_LCAO_NVE_Si8/INPUT_1 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 From a9b4ea93866d5d41a8bd4ab426025ae50e1dcc66 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 15:33:29 +0800 Subject: [PATCH 02/31] update out_chg 2 --- .../output_files/output-specification.md | 6 ++ examples/18_md/01_LCAO_NVE_Si8/INPUT_0 | 33 ------- .../INPUT_1 => 02_LCAO_NVT_Si8/INPUT} | 10 +- source/source_io/module_chgpot/write_init.cpp | 95 ++++++++++--------- source/source_io/module_chgpot/write_init.h | 19 +++- .../read_input_item_output.cpp | 7 +- 6 files changed, 83 insertions(+), 87 deletions(-) delete mode 100644 examples/18_md/01_LCAO_NVE_Si8/INPUT_0 rename examples/18_md/{01_LCAO_NVE_Si8/INPUT_1 => 02_LCAO_NVT_Si8/INPUT} (82%) diff --git a/docs/advanced/output_files/output-specification.md b/docs/advanced/output_files/output-specification.md index 27a9d8376ec..9bb5aa72b40 100644 --- a/docs/advanced/output_files/output-specification.md +++ b/docs/advanced/output_files/output-specification.md @@ -49,11 +49,13 @@ 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`), the naming convention is: `chgg{#}_ini.cube` (nspin=1) or `chgs{#}g{#}_ini.cube` (nspin=2/4) ### 2.2 Examples @@ -62,6 +64,8 @@ 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) | +| `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) | | `pots1.cube` | Local potential, spin 1 | | `eig_occ.txt` | Eigenvalues and occupations | | `doss1g1_nao.txt` | DOS, spin 1, geometry step 1, NAO basis | @@ -79,6 +83,8 @@ 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 | +| `chgg{#}_ini.cube` | Initial charge density, geometry step {#} (nspin=1) | +| `chgs{#}g{#}_ini.cube` | Initial charge density, spin {#}, geometry step {#} (nspin=2/4) | | `taus1.cube`, `taus2.cube` | Kinetic energy density (tau) | | `pots1.cube`, `pots2.cube` | Local potential | diff --git a/examples/18_md/01_LCAO_NVE_Si8/INPUT_0 b/examples/18_md/01_LCAO_NVE_Si8/INPUT_0 deleted file mode 100644 index 5827a8c0178..00000000000 --- a/examples/18_md/01_LCAO_NVE_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_NVE_Si8/INPUT_1 b/examples/18_md/02_LCAO_NVT_Si8/INPUT similarity index 82% rename from examples/18_md/01_LCAO_NVE_Si8/INPUT_1 rename to examples/18_md/02_LCAO_NVT_Si8/INPUT index c49ff60e925..76fb9c54d6a 100644 --- a/examples/18_md/01_LCAO_NVE_Si8/INPUT_1 +++ b/examples/18_md/02_LCAO_NVT_Si8/INPUT @@ -15,7 +15,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 +28,12 @@ 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 +out_freq_ion 2 diff --git a/source/source_io/module_chgpot/write_init.cpp b/source/source_io/module_chgpot/write_init.cpp index 8d41e617693..cfea40d0a9c 100644 --- a/source/source_io/module_chgpot/write_init.cpp +++ b/source/source_io/module_chgpot/write_init.cpp @@ -13,51 +13,58 @@ void ModuleIO::write_chg_init( const Input_para& inp) { const int nspin = inp.nspin; - assert(nspin == 1 || nspin ==2 || nspin == 4); + assert(nspin == 1 || nspin == 2 || nspin == 4); if (inp.out_chg[0] == 2) { - for (int is = 0; is < nspin; is++) + bool should_output = (istep == 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"; + int geom_step = istep + 1; - if(nspin==1) + for (int is = 0; is < nspin; is++) { - ss << "ini.cube"; - } - else if(nspin==2 || nspin==4) - { - ss << "s" << is + 1 << "ini.cube"; - } + std::stringstream ss; + ss << PARAM.globalv.global_out_dir << "chg"; - // 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; - } - } + if (nspin == 1) + { + ss << "g" << geom_step << "_ini.cube"; + } + else if (nspin == 2 || nspin == 4) + { + ss << "s" << is + 1 << "g" << geom_step << "_ini.cube"; + } - ModuleIO::write_vdata_palgrid(para_grid, - chr.rho[is], - is, - nspin, - istep, - ss.str(), - fermi_energy, - &(ucell), - inp.out_chg[1]); + 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]); + } } } return; @@ -71,9 +78,8 @@ void ModuleIO::write_pot_init( const int istep, const Input_para& inp) { - //! 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 (inp.out_pot[0] == 3) { @@ -82,11 +88,11 @@ void ModuleIO::write_pot_init( std::stringstream ss; ss << PARAM.globalv.global_out_dir << "pot"; - if(nspin==1) + if (nspin == 1) { ss << "ini.cube"; } - else if(nspin==2 || nspin==4) + else if (nspin == 2 || nspin == 4) { ss << "s" << is + 1 << "ini.cube"; } @@ -97,11 +103,10 @@ void ModuleIO::write_pot_init( nspin, istep, ss.str(), - 0.0, // efermi + 0.0, &(ucell), - 11, // precsion - 0); // out_fermi + 11, + 0); } } - } diff --git a/source/source_io/module_chgpot/write_init.h b/source/source_io/module_chgpot/write_init.h index 323cbe7af8f..705a3bccb16 100644 --- a/source/source_io/module_chgpot/write_init.h +++ b/source/source_io/module_chgpot/write_init.h @@ -1,14 +1,21 @@ #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 "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 { +// 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 (output at step 0 or every out_freq_ion steps). +// Output file naming convention: +// nspin=1: chgg{geom_step}_ini.cube (e.g., chgg1_ini.cube) +// nspin=2/4: chgs{spin}g{geom_step}_ini.cube (e.g., chgs1g1_ini.cube, chgs2g1_ini.cube) +// Note: geom_step starts from 1 (geom_step = istep + 1). void write_chg_init( const UnitCell& ucell, const Parallel_Grid ¶_grid, @@ -17,6 +24,9 @@ void write_chg_init( const int istep, const Input_para& inp); +// Write initial effective potential to cube file in real space. +// Triggered when inp.out_pot[0] == 3. +// Output: pot_ini.cube (nspin=1) or potsN_ini.cube (nspin=2/4) void write_pot_init( const UnitCell& ucell, const Parallel_Grid ¶_grid, @@ -25,6 +35,5 @@ void write_pot_init( const Input_para& inp); } - #endif 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..15a3d092e3c 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,12 @@ 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: + - 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`. + Here, {geom_step} denotes the geometry step index, starting from 1 (geom_step = istep + 1). + The output frequency is controlled by out_freq_ion (output at step 0 or every out_freq_ion steps). - -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. From c294b594aada23f668019af7120b6a7214382fff Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 15:46:22 +0800 Subject: [PATCH 03/31] update out_pot function --- .../output_files/output-specification.md | 5 ++ examples/18_md/02_LCAO_NVT_Si8/KPT | 4 + examples/18_md/02_LCAO_NVT_Si8/STRU | 28 +++++++ source/source_esolver/esolver_fp.cpp | 6 +- source/source_io/module_chgpot/write_init.cpp | 80 ++++++++++++++----- source/source_io/module_chgpot/write_init.h | 25 +++++- 6 files changed, 124 insertions(+), 24 deletions(-) create mode 100644 examples/18_md/02_LCAO_NVT_Si8/KPT create mode 100644 examples/18_md/02_LCAO_NVT_Si8/STRU diff --git a/docs/advanced/output_files/output-specification.md b/docs/advanced/output_files/output-specification.md index 9bb5aa72b40..4848c1c0863 100644 --- a/docs/advanced/output_files/output-specification.md +++ b/docs/advanced/output_files/output-specification.md @@ -56,6 +56,7 @@ All output file naming conventions can be found in the online documentation (wit - 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`), the naming convention is: `chgg{#}_ini.cube` (nspin=1) or `chgs{#}g{#}_ini.cube` (nspin=2/4) +- For initial potential output (`out_pot = 3`), the naming convention is: `potg{#}_ini.cube` (nspin=1) or `pots{#}g{#}_ini.cube` (nspin=2/4) ### 2.2 Examples @@ -66,6 +67,8 @@ All output file naming conventions can be found in the online documentation (wit | `chgs3.cube` | Charge density, spin 3 (non-collinear with SOC) | | `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) | +| `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 | | `eig_occ.txt` | Eigenvalues and occupations | | `doss1g1_nao.txt` | DOS, spin 1, geometry step 1, NAO basis | @@ -85,6 +88,8 @@ All output file naming conventions can be found in the online documentation (wit | `chg.cube` | Total charge density | | `chgg{#}_ini.cube` | Initial charge density, geometry step {#} (nspin=1) | | `chgs{#}g{#}_ini.cube` | Initial charge density, spin {#}, geometry step {#} (nspin=2/4) | +| `potg{#}_ini.cube` | Initial potential, geometry step {#} (nspin=1) | +| `pots{#}g{#}_ini.cube` | Initial potential, spin {#}, geometry step {#} (nspin=2/4) | | `taus1.cube`, `taus2.cube` | Kinetic energy density (tau) | | `pots1.cube`, `pots2.cube` | Local potential | 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/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 238c093580d..dbdd07d43f7 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -176,8 +176,10 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) 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); + ModuleIO::write_chg_init(ucell, this->Pgrid, this->chr, this->pelec->eferm, istep, + PARAM.globalv.global_out_dir, PARAM.inp); + ModuleIO::write_pot_init(ucell, this->Pgrid, this->pelec, istep, + PARAM.globalv.global_out_dir, PARAM.inp); return; } diff --git a/source/source_io/module_chgpot/write_init.cpp b/source/source_io/module_chgpot/write_init.cpp index cfea40d0a9c..f6dc697e101 100644 --- a/source/source_io/module_chgpot/write_init.cpp +++ b/source/source_io/module_chgpot/write_init.cpp @@ -1,5 +1,22 @@ +// ===================================================================== +// 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: +// - Initial charge density (out_chg = 2): chgg{#}_ini.cube or chgs{#}g{#}_ini.cube +// - Initial potential (out_pot = 3): potg{#}_ini.cube or pots{#}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 @@ -10,11 +27,17 @@ void ModuleIO::write_chg_init( const Charge &chr, const elecstate::Efermi &efermi, const int istep, + const std::string& out_dir, const Input_para& inp) { const int nspin = inp.nspin; 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) { bool should_output = (istep == 0) || @@ -27,7 +50,7 @@ void ModuleIO::write_chg_init( for (int is = 0; is < nspin; is++) { std::stringstream ss; - ss << PARAM.globalv.global_out_dir << "chg"; + ss << out_dir << "chg"; if (nspin == 1) { @@ -76,37 +99,52 @@ void ModuleIO::write_pot_init( const Parallel_Grid ¶_grid, elecstate::ElecState *pelec, const int istep, + const std::string& out_dir, const Input_para& inp) { const int nspin = inp.nspin; 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 = (istep == 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"; + int geom_step = istep + 1; - 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::stringstream ss; + ss << out_dir << "pot"; - ModuleIO::write_vdata_palgrid(para_grid, - pelec->pot->get_eff_v(is), - is, - nspin, - istep, - ss.str(), - 0.0, - &(ucell), - 11, - 0); + if (nspin == 1) + { + ss << "g" << geom_step << "_ini.cube"; + } + else if (nspin == 2 || nspin == 4) + { + ss << "s" << is + 1 << "g" << geom_step << "_ini.cube"; + } + + ModuleIO::write_vdata_palgrid(para_grid, + pelec->pot->get_eff_v(is), + is, + nspin, + istep, + ss.str(), + 0.0, + &(ucell), + 11, + 0); + } } } + return; } diff --git a/source/source_io/module_chgpot/write_init.h b/source/source_io/module_chgpot/write_init.h index 705a3bccb16..b0e7687214f 100644 --- a/source/source_io/module_chgpot/write_init.h +++ b/source/source_io/module_chgpot/write_init.h @@ -1,6 +1,23 @@ +// ===================================================================== +// 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: +// - Initial charge density (out_chg = 2): chgg{#}_ini.cube or chgs{#}g{#}_ini.cube +// - Initial potential (out_pot = 3): potg{#}_ini.cube or pots{#}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 #include "source_io/module_parameter/input_parameter.h" #include "source_estate/module_charge/charge.h" #include "source_estate/fp_energy.h" @@ -22,16 +39,22 @@ void write_chg_init( const Charge &chr, const elecstate::Efermi &efermi, const int istep, + const std::string& out_dir, const Input_para& inp); // Write initial effective potential to cube file in real space. // Triggered when inp.out_pot[0] == 3. -// Output: pot_ini.cube (nspin=1) or potsN_ini.cube (nspin=2/4) +// Output frequency is controlled by out_freq_ion (output at step 0 or every out_freq_ion steps). +// Output file naming convention: +// nspin=1: potg{geom_step}_ini.cube (e.g., potg1_ini.cube) +// nspin=2/4: pots{spin}g{geom_step}_ini.cube (e.g., pots1g1_ini.cube, pots2g1_ini.cube) +// 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 std::string& out_dir, const Input_para& inp); } From d676566940c9083a776b95a4ac417d11e76de746 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 15:59:34 +0800 Subject: [PATCH 04/31] =?UTF-8?q?feat(module=5Fio):=20=E4=BC=98=E5=8C=96?= =?UTF-8?q?=E5=88=9D=E5=A7=8B=E7=94=B5=E8=8D=B7=E5=AF=86=E5=BA=A6/?= =?UTF-8?q?=E5=8A=BF=E8=83=BD=E8=BE=93=E5=87=BA=EF=BC=8C=E6=94=AF=E6=8C=81?= =?UTF-8?q?=20out=5Ffreq=5Fion=20=E6=8E=A7=E5=88=B6=E5=92=8C=E5=8A=A8?= =?UTF-8?q?=E6=80=81=E6=96=87=E4=BB=B6=E5=90=8D?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 添加 gen_ini_filename() 辅助函数，统一生成初始电荷密度/势能文件名 - out_freq_ion=0 时输出单个固定名称文件（不带 g#） - out_freq_ion>0 时每个几何步输出独立文件（带 g#） - 更新文档，说明两种模式的区别 修改文件： - docs/advanced/output_files/output-specification.md - source/source_io/module_chgpot/write_init.cpp - source/source_io/module_chgpot/write_init.h - source/source_io/module_parameter/read_input_item_output.cpp --- .../output_files/output-specification.md | 26 +++++-- source/source_io/module_chgpot/write_init.cpp | 71 +++++++++++-------- source/source_io/module_chgpot/write_init.h | 22 ++++-- .../read_input_item_output.cpp | 24 +++++-- 4 files changed, 97 insertions(+), 46 deletions(-) diff --git a/docs/advanced/output_files/output-specification.md b/docs/advanced/output_files/output-specification.md index 4848c1c0863..30d2a1184ea 100644 --- a/docs/advanced/output_files/output-specification.md +++ b/docs/advanced/output_files/output-specification.md @@ -55,8 +55,12 @@ All output file naming conventions can be found in the online documentation (wit - 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`), the naming convention is: `chgg{#}_ini.cube` (nspin=1) or `chgs{#}g{#}_ini.cube` (nspin=2/4) -- For initial potential output (`out_pot = 3`), the naming convention is: `potg{#}_ini.cube` (nspin=1) or `pots{#}g{#}_ini.cube` (nspin=2/4) +- For initial charge density output (`out_chg = 2`): + - `out_freq_ion = 0`: `chg_ini.cube` (nspin=1) or `chgs{#}_ini.cube` (nspin=2/4) + - `out_freq_ion > 0`: `chgg{#}_ini.cube` (nspin=1) or `chgs{#}g{#}_ini.cube` (nspin=2/4) +- For initial potential output (`out_pot = 3`): + - `out_freq_ion = 0`: `pot_ini.cube` (nspin=1) or `pots{#}_ini.cube` (nspin=2/4) + - `out_freq_ion > 0`: `potg{#}_ini.cube` (nspin=1) or `pots{#}g{#}_ini.cube` (nspin=2/4) ### 2.2 Examples @@ -65,8 +69,12 @@ 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 | @@ -86,10 +94,16 @@ 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 | -| `chgg{#}_ini.cube` | Initial charge density, geometry step {#} (nspin=1) | -| `chgs{#}g{#}_ini.cube` | Initial charge density, spin {#}, geometry step {#} (nspin=2/4) | -| `potg{#}_ini.cube` | Initial potential, geometry step {#} (nspin=1) | -| `pots{#}g{#}_ini.cube` | Initial potential, spin {#}, geometry step {#} (nspin=2/4) | +| 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 | diff --git a/source/source_io/module_chgpot/write_init.cpp b/source/source_io/module_chgpot/write_init.cpp index f6dc697e101..f7eaa6ab88d 100644 --- a/source/source_io/module_chgpot/write_init.cpp +++ b/source/source_io/module_chgpot/write_init.cpp @@ -3,8 +3,8 @@ // to cube files in real space. It is part of the module_io package. // // Output files are named according to the following convention: -// - Initial charge density (out_chg = 2): chgg{#}_ini.cube or chgs{#}g{#}_ini.cube -// - Initial potential (out_pot = 3): potg{#}_ini.cube or pots{#}g{#}_ini.cube +// - 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: @@ -21,6 +21,41 @@ #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, @@ -45,21 +80,11 @@ void ModuleIO::write_chg_init( if (should_output) { - int geom_step = istep + 1; + bool include_geom_step = (inp.out_freq_ion > 0); for (int is = 0; is < nspin; is++) { - std::stringstream ss; - ss << out_dir << "chg"; - - if (nspin == 1) - { - ss << "g" << geom_step << "_ini.cube"; - } - else if (nspin == 2 || nspin == 4) - { - ss << "s" << is + 1 << "g" << geom_step << "_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) @@ -83,7 +108,7 @@ void ModuleIO::write_chg_init( is, nspin, istep, - ss.str(), + filename, fermi_energy, &(ucell), inp.out_chg[1]); @@ -117,28 +142,18 @@ void ModuleIO::write_pot_init( if (should_output) { - int geom_step = istep + 1; + bool include_geom_step = (inp.out_freq_ion > 0); for (int is = 0; is < nspin; is++) { - std::stringstream ss; - ss << out_dir << "pot"; - - if (nspin == 1) - { - ss << "g" << geom_step << "_ini.cube"; - } - else if (nspin == 2 || nspin == 4) - { - ss << "s" << is + 1 << "g" << geom_step << "_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(), + filename, 0.0, &(ucell), 11, diff --git a/source/source_io/module_chgpot/write_init.h b/source/source_io/module_chgpot/write_init.h index b0e7687214f..fd4c5506f8f 100644 --- a/source/source_io/module_chgpot/write_init.h +++ b/source/source_io/module_chgpot/write_init.h @@ -3,8 +3,8 @@ // to cube files in real space. It is part of the module_io package. // // Output files are named according to the following convention: -// - Initial charge density (out_chg = 2): chgg{#}_ini.cube or chgs{#}g{#}_ini.cube -// - Initial potential (out_pot = 3): potg{#}_ini.cube or pots{#}g{#}_ini.cube +// - 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: @@ -26,12 +26,22 @@ 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 (output at step 0 or every out_freq_ion steps). // Output file naming convention: -// nspin=1: chgg{geom_step}_ini.cube (e.g., chgg1_ini.cube) -// nspin=2/4: chgs{spin}g{geom_step}_ini.cube (e.g., chgs1g1_ini.cube, chgs2g1_ini.cube) +// 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, @@ -46,8 +56,8 @@ void write_chg_init( // Triggered when inp.out_pot[0] == 3. // Output frequency is controlled by out_freq_ion (output at step 0 or every out_freq_ion steps). // Output file naming convention: -// nspin=1: potg{geom_step}_ini.cube (e.g., potg1_ini.cube) -// nspin=2/4: pots{spin}g{geom_step}_ini.cube (e.g., pots1g1_ini.cube, pots2g1_ini.cube) +// 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, 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 15a3d092e3c..827a56e9af1 100644 --- a/source/source_io/module_parameter/read_input_item_output.cpp +++ b/source/source_io/module_parameter/read_input_item_output.cpp @@ -81,9 +81,14 @@ void ReadInput::item_output() - 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. The files are named as: - - 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`. + - 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`; + - 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`. Here, {geom_step} denotes the geometry step index, starting from 1 (geom_step = istep + 1). The output frequency is controlled by out_freq_ion (output at step 0 or every out_freq_ion steps). - -1: Disable the charge density auto-back-up file `{suffix}-CHARGE-DENSITY.restart`, useful for large systems. @@ -125,9 +130,16 @@ 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`; + * 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`. + Here, {geom_step} denotes the geometry step index, starting from 1 (geom_step = istep + 1). + The output frequency is controlled by out_freq_ion (output at step 0 or every out_freq_ion steps). The optional second integer controls the output precision. If not provided, the default precision is 8. From a48c058b769d2af36f81251d05067113f50b6883 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 16:06:47 +0800 Subject: [PATCH 05/31] =?UTF-8?q?fix(module=5Fio):=20=E4=BF=AE=E6=AD=A3=20?= =?UTF-8?q?out=5Ffreq=5Fion=3D0=20=E6=97=B6=E5=88=9D=E5=A7=8B=E7=94=B5?= =?UTF-8?q?=E8=8D=B7=E5=AF=86=E5=BA=A6/=E5=8A=BF=E8=83=BD=E8=BE=93?= =?UTF-8?q?=E5=87=BA=E9=80=BB=E8=BE=91?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - out_freq_ion=0 时，每个几何步都输出（覆盖同一个文件） - out_freq_ion>0 时，只在 istep 是 out_freq_ion 倍数时输出 - 更新所有相关文档和注释 修改文件： - docs/advanced/output_files/output-specification.md - source/source_io/module_chgpot/write_init.cpp - source/source_io/module_chgpot/write_init.h - source/source_io/module_parameter/read_input_item_output.cpp --- docs/advanced/output_files/output-specification.md | 8 ++++---- source/source_io/module_chgpot/write_init.cpp | 4 ++-- source/source_io/module_chgpot/write_init.h | 8 ++++++-- .../source_io/module_parameter/read_input_item_output.cpp | 6 ++++-- 4 files changed, 16 insertions(+), 10 deletions(-) diff --git a/docs/advanced/output_files/output-specification.md b/docs/advanced/output_files/output-specification.md index 30d2a1184ea..d677e3dba56 100644 --- a/docs/advanced/output_files/output-specification.md +++ b/docs/advanced/output_files/output-specification.md @@ -56,11 +56,11 @@ All output file naming conventions can be found in the online documentation (wit - 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) - - `out_freq_ion > 0`: `chgg{#}_ini.cube` (nspin=1) or `chgs{#}g{#}_ini.cube` (nspin=2/4) + - `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) - - `out_freq_ion > 0`: `potg{#}_ini.cube` (nspin=1) or `pots{#}g{#}_ini.cube` (nspin=2/4) + - `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 diff --git a/source/source_io/module_chgpot/write_init.cpp b/source/source_io/module_chgpot/write_init.cpp index f7eaa6ab88d..b8e6f3fbd0d 100644 --- a/source/source_io/module_chgpot/write_init.cpp +++ b/source/source_io/module_chgpot/write_init.cpp @@ -75,7 +75,7 @@ void ModuleIO::write_chg_init( if (inp.out_chg[0] == 2) { - bool should_output = (istep == 0) || + bool should_output = (inp.out_freq_ion == 0) || (inp.out_freq_ion > 0 && istep % inp.out_freq_ion == 0); if (should_output) @@ -137,7 +137,7 @@ void ModuleIO::write_pot_init( if (inp.out_pot[0] == 3) { - bool should_output = (istep == 0) || + bool should_output = (inp.out_freq_ion == 0) || (inp.out_freq_ion > 0 && istep % inp.out_freq_ion == 0); if (should_output) diff --git a/source/source_io/module_chgpot/write_init.h b/source/source_io/module_chgpot/write_init.h index fd4c5506f8f..7d5ecd5239b 100644 --- a/source/source_io/module_chgpot/write_init.h +++ b/source/source_io/module_chgpot/write_init.h @@ -38,7 +38,9 @@ std::string gen_ini_filename( // 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 (output at step 0 or every out_freq_ion steps). +// 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) @@ -54,7 +56,9 @@ void write_chg_init( // 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 (output at step 0 or every out_freq_ion steps). +// 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) 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 827a56e9af1..074ce3df679 100644 --- a/source/source_io/module_parameter/read_input_item_output.cpp +++ b/source/source_io/module_parameter/read_input_item_output.cpp @@ -85,12 +85,13 @@ void ReadInput::item_output() - 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). - The output frequency is controlled by out_freq_ion (output at step 0 or every out_freq_ion steps). - -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. @@ -134,12 +135,13 @@ In molecular dynamics simulations, the output frequency is controlled by out_fre * 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 output frequency is controlled by out_freq_ion (output at step 0 or every out_freq_ion steps). The optional second integer controls the output precision. If not provided, the default precision is 8. From 8378c0e5e1e39c71c740b1ae79264ec0614aae5e Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 16:17:29 +0800 Subject: [PATCH 06/31] fix a bug about out_pot --- source/source_io/module_ctrl/ctrl_output_fp.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index 710f4b87307..c23db5a990a 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -129,7 +129,7 @@ void ctrl_output_fp(UnitCell& ucell, pelec->pot->get_eff_v(is), is, nspin, - istep_in, + istep, fn, 0.0, // efermi &(ucell), From da4eefb4eedd56d2c6ade6f865fb7c309b7ec47c Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 16:40:37 +0800 Subject: [PATCH 07/31] fix bugs --- source/source_io/module_chgpot/write_init.cpp | 2 +- .../source_io/module_ctrl/ctrl_output_fp.cpp | 129 ++++++++++-------- 2 files changed, 71 insertions(+), 60 deletions(-) diff --git a/source/source_io/module_chgpot/write_init.cpp b/source/source_io/module_chgpot/write_init.cpp index b8e6f3fbd0d..b455c3833dd 100644 --- a/source/source_io/module_chgpot/write_init.cpp +++ b/source/source_io/module_chgpot/write_init.cpp @@ -156,7 +156,7 @@ void ModuleIO::write_pot_init( filename, 0.0, &(ucell), - 11, + inp.out_pot[1], 0); } } diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index c23db5a990a..67ff5e67de2 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 @@ -61,47 +60,53 @@ void ctrl_output_fp(UnitCell& ucell, // 4) write charge density if (PARAM.inp.out_chg[0] > 0) { - for (int is = 0; is < nspin; ++is) + // Only output when: + // - out_freq_ion == 0: output at every step (no g# in filename) but only one file + // - out_freq_ion > 0: output only when istep % out_freq_ion == 0 (with g# in filename), many files + if (PARAM.inp.out_freq_ion == 0 || istep % PARAM.inp.out_freq_ion == 0) { - std::string fn = PARAM.globalv.global_out_dir + "chg"; - - std::string spin_block; - if (nspin == 2 || nspin == 4) - { - spin_block = "s" + std::to_string(is + 1); - } - else if (nspin == 1) - { - // do nothing - } - - fn += spin_block + geom_block + ".cube"; - - ModuleIO::write_vdata_palgrid(para_grid, - chr.rho_save[is], - is, - nspin, - istep, // change istep_in to istep, mohan 20260222 - fn, - pelec->eferm.get_efval(is), - &(ucell), - PARAM.inp.out_chg[1], - 1); - - if (XC_Functional::get_ked_flag()) + for (int is = 0; is < nspin; ++is) { - fn = PARAM.globalv.global_out_dir + "tau"; + std::string fn = PARAM.globalv.global_out_dir + "chg"; + + std::string spin_block; + if (nspin == 2 || nspin == 4) + { + spin_block = "s" + std::to_string(is + 1); + } + else if (nspin == 1) + { + // do nothing + } fn += spin_block + geom_block + ".cube"; ModuleIO::write_vdata_palgrid(para_grid, - chr.kin_r_save[is], + chr.rho_save[is], is, nspin, - istep, + istep, // change istep_in to istep, mohan 20260222 fn, pelec->eferm.get_efval(is), - &(ucell)); + &(ucell), + PARAM.inp.out_chg[1], + 1); + + if (XC_Functional::get_ked_flag()) + { + fn = PARAM.globalv.global_out_dir + "tau"; + + fn += spin_block + geom_block + ".cube"; + + ModuleIO::write_vdata_palgrid(para_grid, + chr.kin_r_save[is], + is, + nspin, + istep, + fn, + pelec->eferm.get_efval(is), + &(ucell)); + } } } } @@ -109,32 +114,38 @@ void ctrl_output_fp(UnitCell& ucell, // 5) write potential if (PARAM.inp.out_pot[0] == 1 || PARAM.inp.out_pot[0] == 3) { - for (int is = 0; is < nspin; is++) + // Only output when: + // - out_freq_ion == 0: output at every step (no g# in filename) but only one file + // - out_freq_ion > 0: output only when istep % out_freq_ion == 0 (with g# in filename) but many files + if (PARAM.inp.out_freq_ion == 0 || istep % PARAM.inp.out_freq_ion == 0) { - std::string fn = PARAM.globalv.global_out_dir + "pot"; - - std::string spin_block; - if (nspin == 2 || nspin == 4) - { - spin_block = "s" + std::to_string(is + 1); - } - else if (nspin == 1) + for (int is = 0; is < nspin; is++) { - // do nothing - } + std::string fn = PARAM.globalv.global_out_dir + "pot"; + + std::string spin_block; + if (nspin == 2 || nspin == 4) + { + spin_block = "s" + std::to_string(is + 1); + } + else if (nspin == 1) + { + // do nothing + } + + fn += spin_block + geom_block + ".cube"; - fn += spin_block + geom_block + ".cube"; - - ModuleIO::write_vdata_palgrid(para_grid, - pelec->pot->get_eff_v(is), - is, - nspin, - istep, - fn, - 0.0, // efermi - &(ucell), - PARAM.inp.out_pot[1], // precision - 0); // out_fermi + ModuleIO::write_vdata_palgrid(para_grid, + pelec->pot->get_eff_v(is), + is, + nspin, + istep, + fn, + 0.0, // efermi + &(ucell), + PARAM.inp.out_pot[1], // precision + 0); // out_fermi + } } } else if (PARAM.inp.out_pot[0] == 2) @@ -154,7 +165,7 @@ 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 From 17177461f6244b0cfcb59abd27623c5844418d9e Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 17:51:48 +0800 Subject: [PATCH 08/31] update --- examples/18_md/02_LCAO_NVT_Si8/INPUT | 6 +- .../source_io/module_ctrl/ctrl_output_fp.cpp | 150 ++++++++---------- source/source_io/module_elf/write_elf.cpp | 27 ++-- source/source_io/module_elf/write_elf.h | 7 +- 4 files changed, 84 insertions(+), 106 deletions(-) diff --git a/examples/18_md/02_LCAO_NVT_Si8/INPUT b/examples/18_md/02_LCAO_NVT_Si8/INPUT index 76fb9c54d6a..b0c75b94e50 100644 --- a/examples/18_md/02_LCAO_NVT_Si8/INPUT +++ b/examples/18_md/02_LCAO_NVT_Si8/INPUT @@ -35,5 +35,9 @@ md_tfreq 1 md_tchain 1 # output initial charge density -out_chg 2 +out_chg 2 4 +out_pot 3 4 +out_elf 1 4 +out_dipole 1 out_freq_ion 2 + diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index 67ff5e67de2..a36dae61750 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -48,107 +48,93 @@ 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) { - // Only output when: - // - out_freq_ion == 0: output at every step (no g# in filename) but only one file - // - out_freq_ion > 0: output only when istep % out_freq_ion == 0 (with g# in filename), many files - if (PARAM.inp.out_freq_ion == 0 || istep % PARAM.inp.out_freq_ion == 0) + for (int is = 0; is < nspin; ++is) { - for (int is = 0; is < nspin; ++is) + std::string fn = PARAM.globalv.global_out_dir + "chg"; + + std::string spin_block; + if (nspin == 2 || nspin == 4) { - std::string fn = PARAM.globalv.global_out_dir + "chg"; - - std::string spin_block; - if (nspin == 2 || nspin == 4) - { - spin_block = "s" + std::to_string(is + 1); - } - else if (nspin == 1) - { - // do nothing - } + spin_block = "s" + std::to_string(is + 1); + } + else if (nspin == 1) + { + // do nothing + } + + fn += spin_block + geom_block + ".cube"; + + ModuleIO::write_vdata_palgrid(para_grid, + chr.rho_save[is], + is, + nspin, + istep, // change istep_in to istep, mohan 20260222 + fn, + pelec->eferm.get_efval(is), + &(ucell), + PARAM.inp.out_chg[1], + 1); + + if (XC_Functional::get_ked_flag()) + { + fn = PARAM.globalv.global_out_dir + "tau"; fn += spin_block + geom_block + ".cube"; ModuleIO::write_vdata_palgrid(para_grid, - chr.rho_save[is], + chr.kin_r_save[is], is, nspin, - istep, // change istep_in to istep, mohan 20260222 + istep, fn, pelec->eferm.get_efval(is), - &(ucell), - PARAM.inp.out_chg[1], - 1); - - if (XC_Functional::get_ked_flag()) - { - fn = PARAM.globalv.global_out_dir + "tau"; - - fn += spin_block + geom_block + ".cube"; - - ModuleIO::write_vdata_palgrid(para_grid, - chr.kin_r_save[is], - is, - nspin, - istep, - fn, - pelec->eferm.get_efval(is), - &(ucell)); - } + &(ucell)); } } } // 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) { - // Only output when: - // - out_freq_ion == 0: output at every step (no g# in filename) but only one file - // - out_freq_ion > 0: output only when istep % out_freq_ion == 0 (with g# in filename) but many files - if (PARAM.inp.out_freq_ion == 0 || istep % PARAM.inp.out_freq_ion == 0) + for (int is = 0; is < nspin; is++) { - for (int is = 0; is < nspin; is++) - { - std::string fn = PARAM.globalv.global_out_dir + "pot"; - - std::string spin_block; - if (nspin == 2 || nspin == 4) - { - spin_block = "s" + std::to_string(is + 1); - } - else if (nspin == 1) - { - // do nothing - } + std::string fn = PARAM.globalv.global_out_dir + "pot"; - fn += spin_block + geom_block + ".cube"; - - ModuleIO::write_vdata_palgrid(para_grid, - pelec->pot->get_eff_v(is), - is, - nspin, - istep, - fn, - 0.0, // efermi - &(ucell), - PARAM.inp.out_pot[1], // precision - 0); // out_fermi + std::string spin_block; + if (nspin == 2 || nspin == 4) + { + spin_block = "s" + std::to_string(is + 1); } + else if (nspin == 1) + { + // do nothing + } + + fn += spin_block + geom_block + ".cube"; + + ModuleIO::write_vdata_palgrid(para_grid, + pelec->pot->get_eff_v(is), + is, + nspin, + istep, + fn, + 0.0, // efermi + &(ucell), + PARAM.inp.out_pot[1], // precision + 0); // out_fermi } } - 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"; @@ -169,7 +155,7 @@ void ctrl_output_fp(UnitCell& ucell, } // 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; @@ -179,12 +165,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, @@ -192,12 +173,13 @@ void ctrl_output_fp(UnitCell& ucell, pw_rhod, para_grid, &(ucell), - PARAM.inp.out_elf[1]); + PARAM.inp.out_elf[1], + geom_block); } #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(), @@ -211,7 +193,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_elf/write_elf.cpp b/source/source_io/module_elf/write_elf.cpp index 347a119099b..6f95b056427 100644 --- a/source/source_io/module_elf/write_elf.cpp +++ b/source/source_io/module_elf/write_elf.cpp @@ -4,19 +4,16 @@ 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) { // 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. @@ -24,7 +21,7 @@ void write_elf( // 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. - int nspin_eff = (nspin == 4) ? 1 : nspin; + const int nspin_eff = (nspin == 4) ? 1 : nspin; std::vector> elf(nspin_eff, std::vector(rho_basis->nrxx, 0.)); // 1) calculate the kinetic energy density of vW KEDF @@ -100,7 +97,7 @@ 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) + const 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) for (int is = 0; is < nspin_eff; ++is) { for (int ir = 0; ir < rho_basis->nrxx; ++ir) @@ -118,12 +115,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 + "/elf" + geom_block + ".cube"; int is = -1; ModuleIO::write_vdata_palgrid(pgrid, @@ -141,11 +138,9 @@ void write_elf( { for (int is = 0; is < nspin; ++is) { - std::string fn_temp = out_dir + "/elf"; + std::string fn_temp = out_dir + "/elf" + std::to_string(is + 1) + geom_block + ".cube"; - fn_temp += std::to_string(is + 1) + ".cube"; - - int ispin = is + 1; + const int ispin = is + 1; ModuleIO::write_vdata_palgrid(pgrid, elf[is].data(), @@ -165,7 +160,7 @@ void write_elf( 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]); } - std::string fn = out_dir + "/elf.cube"; + std::string fn = out_dir + "/elf" + geom_block + ".cube"; int is = -1; ModuleIO::write_vdata_palgrid(pgrid, diff --git a/source/source_io/module_elf/write_elf.h b/source/source_io/module_elf/write_elf.h index 23a7e7c379e..7661a379b11 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,8 @@ 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); } #endif From 3137f0ddc1dd4b4aa2f118dd7d9b15872c985ab8 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 18:21:30 +0800 Subject: [PATCH 09/31] update ELF and add openmp parallel --- source/source_io/module_elf/write_elf.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/source/source_io/module_elf/write_elf.cpp b/source/source_io/module_elf/write_elf.cpp index 6f95b056427..050e99b11d7 100644 --- a/source/source_io/module_elf/write_elf.cpp +++ b/source/source_io/module_elf/write_elf.cpp @@ -1,5 +1,8 @@ #include "write_elf.h" #include "source_io/module_output/cube_io.h" +#ifdef _OPENMP +#include +#endif namespace ModuleIO { @@ -30,6 +33,7 @@ void write_elf( for (int is = 0; is < nspin_eff; ++is) { +#pragma omp parallel for schedule(static) default(none) shared(phi, rho, rho_basis, is) for (int ir = 0; ir < rho_basis->nrxx; ++ir) { phi[ir] = std::sqrt(std::abs(rho[is][ir])); @@ -44,6 +48,7 @@ void write_elf( std::complex img(0.0, 1.0); for (int j = 0; j < 3; ++j) { +#pragma omp parallel for schedule(static) default(none) shared(recip_gradient_phi, img, rho_basis, recip_phi, j) for (int ip = 0; ip < rho_basis->npw; ++ip) { recip_gradient_phi[ip] = img * rho_basis->gcar[ip][j] * recip_phi[ip] * rho_basis->tpiba; @@ -51,6 +56,7 @@ void write_elf( rho_basis->recip2real(recip_gradient_phi.data(), gradient_phi[j].data()); +#pragma omp parallel for schedule(static) default(none) shared(tau_vw, gradient_phi, is, j, rho_basis) for (int ir = 0; ir < rho_basis->nrxx; ++ir) { tau_vw[is][ir] += gradient_phi[j][ir] * gradient_phi[j][ir] / 2. * 2.; // convert Ha to Ry. @@ -65,6 +71,7 @@ void write_elf( * 2.0; // 10/3*(3*pi^2)^{2/3}, multiply by 2 to convert unit from Hartree to Ry, finally in Ry*Bohr^(-2) if (nspin == 1 || nspin == 4) { +#pragma omp parallel for schedule(static) default(none) shared(rho, tau_TF, c_tf, rho_basis) for (int ir = 0; ir < rho_basis->nrxx; ++ir) { if (rho[0][ir] > 0.0) @@ -82,6 +89,7 @@ void write_elf( // the spin-scaling law: 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) { +#pragma omp parallel for schedule(static) default(none) shared(rho, tau_TF, c_tf, is, rho_basis) for (int ir = 0; ir < rho_basis->nrxx; ++ir) { if (rho[is][ir] > 0.0) @@ -100,6 +108,7 @@ void write_elf( const 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) for (int is = 0; is < nspin_eff; ++is) { +#pragma omp parallel for schedule(static) default(none) shared(elf, tau, tau_vw, tau_TF, eps, is, rho_basis) for (int ir = 0; ir < rho_basis->nrxx; ++ir) { if (tau_TF[is][ir] > 1.0e-12) @@ -155,6 +164,7 @@ void write_elf( } std::vector elf_tot(rho_basis->nrxx, 0.0); +#pragma omp parallel for schedule(static) default(none) shared(elf_tot, tau, tau_vw, tau_TF, rho_basis) for (int ir = 0; ir < rho_basis->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]); From 94ea4ab3e2a7a919a6dd35615e6b18aeec3ab048 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 18:38:10 +0800 Subject: [PATCH 10/31] update elf --- examples/18_md/02_LCAO_NVT_Si8/INPUT | 1 + source/source_io/module_elf/write_elf.cpp | 124 ++++++++++++++-------- 2 files changed, 82 insertions(+), 43 deletions(-) diff --git a/examples/18_md/02_LCAO_NVT_Si8/INPUT b/examples/18_md/02_LCAO_NVT_Si8/INPUT index b0c75b94e50..1e97e86497c 100644 --- a/examples/18_md/02_LCAO_NVT_Si8/INPUT +++ 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 diff --git a/source/source_io/module_elf/write_elf.cpp b/source/source_io/module_elf/write_elf.cpp index 050e99b11d7..f7133671896 100644 --- a/source/source_io/module_elf/write_elf.cpp +++ b/source/source_io/module_elf/write_elf.cpp @@ -18,61 +18,78 @@ void write_elf( const int& precision, const std::string& geom_block) { - // 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. + // 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; - std::vector> elf(nspin_eff, std::vector(rho_basis->nrxx, 0.)); + const int nrxx = rho_basis->nrxx; + const int npw = rho_basis->npw; + + 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) { -#pragma omp parallel for schedule(static) default(none) shared(phi, rho, rho_basis, 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) { -#pragma omp parallel for schedule(static) default(none) shared(recip_gradient_phi, img, rho_basis, recip_phi, 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()); -#pragma omp parallel for schedule(static) default(none) shared(tau_vw, gradient_phi, is, j, rho_basis) - 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) { -#pragma omp parallel for schedule(static) default(none) shared(rho, tau_TF, c_tf, rho_basis) - 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) { @@ -86,15 +103,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) { -#pragma omp parallel for schedule(static) default(none) shared(rho, tau_TF, c_tf, is, rho_basis) - 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 { @@ -104,16 +125,20 @@ void write_elf( } } - // 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; // 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) { -#pragma omp parallel for schedule(static) default(none) shared(elf, tau, tau_vw, tau_TF, eps, is, rho_basis) - 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 @@ -141,13 +166,14 @@ void write_elf( ef_tmp, ucell_, precision, - out_fermi); + out_fermi); } else if (nspin == 2) { for (int is = 0; is < nspin; ++is) { - std::string fn_temp = out_dir + "/elf" + std::to_string(is + 1) + geom_block + ".cube"; + std::string fn_temp = out_dir + "/elf" + + std::to_string(is + 1) + geom_block + ".cube"; const int ispin = is + 1; @@ -160,15 +186,27 @@ void write_elf( ef_tmp, ucell_, precision, - out_fermi); + out_fermi); } - std::vector elf_tot(rho_basis->nrxx, 0.0); -#pragma omp parallel for schedule(static) default(none) shared(elf_tot, tau, tau_vw, tau_TF, rho_basis) - 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" + geom_block + ".cube"; From 097ac7ced43c5d84a1a49bdc1e3072efbec29e3c Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 18:50:44 +0800 Subject: [PATCH 11/31] update elf --- docs/advanced/output_files/output-specification.md | 14 ++++++++++++++ examples/18_md/02_LCAO_NVT_Si8/INPUT | 1 - source/source_io/module_elf/write_elf.cpp | 6 +++--- .../module_parameter/read_input_item_output.cpp | 10 ++++++---- 4 files changed, 23 insertions(+), 8 deletions(-) diff --git a/docs/advanced/output_files/output-specification.md b/docs/advanced/output_files/output-specification.md index d677e3dba56..834d18b6870 100644 --- a/docs/advanced/output_files/output-specification.md +++ b/docs/advanced/output_files/output-specification.md @@ -78,6 +78,13 @@ All output file naming conventions can be found in the online documentation (wit | `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 | @@ -106,6 +113,13 @@ All output file naming conventions can be found in the online documentation (wit | `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/02_LCAO_NVT_Si8/INPUT b/examples/18_md/02_LCAO_NVT_Si8/INPUT index 1e97e86497c..c07ca1f9f4a 100644 --- a/examples/18_md/02_LCAO_NVT_Si8/INPUT +++ b/examples/18_md/02_LCAO_NVT_Si8/INPUT @@ -41,4 +41,3 @@ out_pot 3 4 out_elf 1 4 out_dipole 1 out_freq_ion 2 - diff --git a/source/source_io/module_elf/write_elf.cpp b/source/source_io/module_elf/write_elf.cpp index f7133671896..9813a833876 100644 --- a/source/source_io/module_elf/write_elf.cpp +++ b/source/source_io/module_elf/write_elf.cpp @@ -154,7 +154,7 @@ void write_elf( if (nspin == 1 || nspin == 4) { - std::string fn = out_dir + "/elf" + geom_block + ".cube"; + std::string fn = out_dir + "/elftot" + geom_block + ".cube"; int is = -1; ModuleIO::write_vdata_palgrid(pgrid, @@ -172,7 +172,7 @@ void write_elf( { for (int is = 0; is < nspin; ++is) { - std::string fn_temp = out_dir + "/elf" + std::string fn_temp = out_dir + "/elf" + "s" + std::to_string(is + 1) + geom_block + ".cube"; const int ispin = is + 1; @@ -208,7 +208,7 @@ void write_elf( elf_tot[ir] = 0.0; } } - std::string fn = out_dir + "/elf" + geom_block + ".cube"; + std::string fn = out_dir + "/elftot" + geom_block + ".cube"; int is = -1; ModuleIO::write_vdata_palgrid(pgrid, 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 074ce3df679..a6b62d4947d 100644 --- a/source/source_io/module_parameter/read_input_item_output.cpp +++ b/source/source_io/module_parameter/read_input_item_output.cpp @@ -985,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. From bbe91cf0dbdb3299a76d2f1ec3b687c80585f72c Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 18:55:08 +0800 Subject: [PATCH 12/31] update example reference data --- tests/01_PW/nscf_out_pot/INPUT | 2 +- tests/01_PW/nscf_out_pot/pot.cube.ref | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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 From 1a20bb9f9cc05f87e303b6c485505f284024bad7 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 19:03:46 +0800 Subject: [PATCH 13/31] enable elf for rt-tddft, but results are wrong --- examples/22_rt-tddft/02_H2_velocity_gauge/INPUT | 5 ++++- source/source_io/module_parameter/read_input_item_output.cpp | 4 ++-- 2 files changed, 6 insertions(+), 3 deletions(-) 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_io/module_parameter/read_input_item_output.cpp b/source/source_io/module_parameter/read_input_item_output.cpp index a6b62d4947d..8bc0c895499 100644 --- a/source/source_io/module_parameter/read_input_item_output.cpp +++ b/source/source_io/module_parameter/read_input_item_output.cpp @@ -1012,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); From 2372ac56c86511d958cbff2ea0d0a33ab3b441b1 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Mon, 25 May 2026 08:22:50 +0800 Subject: [PATCH 14/31] fix elf test --- tests/01_PW/scf_out_elf/{refelf.cube => refelftot.cube} | 0 tests/integrate/tools/catch_properties.sh | 4 ++-- 2 files changed, 2 insertions(+), 2 deletions(-) rename tests/01_PW/scf_out_elf/{refelf.cube => refelftot.cube} (100%) 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/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 From 3717aa17fced2c1834bc92c37dcf79062520b298 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Mon, 25 May 2026 08:25:02 +0800 Subject: [PATCH 15/31] fix elf test in 03_NAO_multik --- tests/03_NAO_multik/scf_out_elf/{refelf.cube => refelftot.cube} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename tests/03_NAO_multik/scf_out_elf/{refelf.cube => refelftot.cube} (100%) 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 From 498180a2d439865e3bc7e4aecf6737c3e7d70838 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Mon, 25 May 2026 08:36:28 +0800 Subject: [PATCH 16/31] fix output of write_elf --- source/source_io/module_elf/write_elf.cpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/source/source_io/module_elf/write_elf.cpp b/source/source_io/module_elf/write_elf.cpp index 9813a833876..daf830db623 100644 --- a/source/source_io/module_elf/write_elf.cpp +++ b/source/source_io/module_elf/write_elf.cpp @@ -18,6 +18,7 @@ void write_elf( const int& precision, const std::string& geom_block) { + 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. @@ -32,6 +33,9 @@ void write_elf( const int nrxx = rho_basis->nrxx; const int npw = rho_basis->npw; + assert(nrxx>0); + assert(npw>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(nrxx, 0.)); @@ -154,7 +158,7 @@ void write_elf( if (nspin == 1 || nspin == 4) { - std::string fn = out_dir + "/elftot" + geom_block + ".cube"; + std::string fn = out_dir + "elftot" + geom_block + ".cube"; int is = -1; ModuleIO::write_vdata_palgrid(pgrid, @@ -172,7 +176,7 @@ void write_elf( { for (int is = 0; is < nspin; ++is) { - std::string fn_temp = out_dir + "/elf" + "s" + std::string fn_temp = out_dir + "elf" + "s" + std::to_string(is + 1) + geom_block + ".cube"; const int ispin = is + 1; @@ -208,7 +212,7 @@ void write_elf( elf_tot[ir] = 0.0; } } - std::string fn = out_dir + "/elftot" + geom_block + ".cube"; + std::string fn = out_dir + "elftot" + geom_block + ".cube"; int is = -1; ModuleIO::write_vdata_palgrid(pgrid, @@ -222,5 +226,6 @@ void write_elf( precision, out_fermi); } -} -} + ModuleBase::timer::end("ModuleIO", "write_elf"); +} // end write_elf +} // end namespace ModuleIO From 9df0aa24355d5c0c12e03ba108ba81e787327297 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Mon, 25 May 2026 21:07:38 +0800 Subject: [PATCH 17/31] fix bug --- tests/03_NAO_multik/scf_out_chg_pot1/pot.cube.ref | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From ea06c1453b93cc18c0662fe97357008ee74ffaaa Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Mon, 25 May 2026 21:08:22 +0800 Subject: [PATCH 18/31] update potential file, fix bug --- tests/03_NAO_multik/nscf_out_pot1/pot.cube.ref | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From dab9525109e032f357c601299a8b5801c49a62d2 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 26 May 2026 08:18:55 +0800 Subject: [PATCH 19/31] fix elf test in ofdft --- tests/07_OFDFT/24_OF_out_elf/{refelf.cube => refelftot.cube} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename tests/07_OFDFT/24_OF_out_elf/{refelf.cube => refelftot.cube} (100%) 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 From 8cc7d841562d41bdc81ddd353b99ff07f42a482f Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 29 May 2026 17:24:52 +0800 Subject: [PATCH 20/31] fix: Move write_pot_init to ElecState::init_scf for correct timing The write_pot_init was being called in ESolver_FP::before_scf before the effective potential was computed. This caused pot_ini.cube to contain: - All zeros for calculation=scf / first ionic step (istep=0) - Converged potential from previous ionic step for relax/md with istep>0 The fix moves write_pot_init to ElecState::init_scf, which is called after pot->init_pot(charge) computes the effective potential from the initial charge density. This ensures pot_ini.cube correctly contains the effective potential corresponding to the initial charge density. Changes: - Modified ElecState::init_scf signature to accept istep, out_dir, inp parameters - Added write_pot_init call after pot->init_pot() in init_scf - Updated pw::setup_pot to pass through the new parameters - Updated all callers (LCAO and PW) to provide the new parameters - Removed the premature write_pot_init call from ESolver_FP::before_scf --- source/source_esolver/esolver_fp.cpp | 4 +--- source/source_esolver/esolver_ks_lcao.cpp | 3 ++- source/source_esolver/esolver_ks_pw.cpp | 2 +- source/source_estate/elecstate.cpp | 9 ++++++++- source/source_estate/elecstate.h | 8 +++++++- source/source_pw/module_pwdft/setup_pot.cpp | 7 ++++++- source/source_pw/module_pwdft/setup_pot.h | 1 + 7 files changed, 26 insertions(+), 8 deletions(-) diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index dbdd07d43f7..6a13a4a5095 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -175,11 +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 + //! 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); - ModuleIO::write_pot_init(ucell, this->Pgrid, this->pelec, istep, - PARAM.globalv.global_out_dir, PARAM.inp); return; } diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 9d3906d6ebf..a4931243ac9 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); + this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm, + istep, PARAM.globalv.global_out_dir, PARAM.inp); #ifdef __MLALGO // 14) initialize DM2(R) of DeePKS, the DM2(R) is different from DM(R) 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_estate/elecstate.cpp b/source/source_estate/elecstate.cpp index dc2ebc7fd22..f6dbeef455f 100644 --- a/source/source_estate/elecstate.cpp +++ b/source/source_estate/elecstate.cpp @@ -3,6 +3,7 @@ #include "source_base/parallel_reduce.h" #include "source_base/tool_title.h" #include "occupy.h" +#include "source_io/module_chgpot/write_init.h" namespace elecstate { @@ -29,7 +30,10 @@ void ElecState::init_scf(const UnitCell& ucell, const ModuleBase::ComplexMatrix& strucfac, const bool* numeric, ModuleSymmetry::Symmetry& symm, - const void* wfcpw) + const void* wfcpw, + const int istep, + const std::string& out_dir, + const Input_para& inp) { //! core correction potential. this->charge->set_rho_core(ucell,strucfac, numeric); @@ -39,6 +43,9 @@ void ElecState::init_scf(const UnitCell& ucell, //! initialize the potential this->pot->init_pot(this->charge); + + //! output the initial potential + ModuleIO::write_pot_init(ucell, pgrid, this, istep, out_dir, inp); } diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h index 888b3b95dc0..49607f306fb 100644 --- a/source/source_estate/elecstate.h +++ b/source/source_estate/elecstate.h @@ -101,13 +101,19 @@ class ElecState * @param strucfac structure factor * @param symm symmetry * @param wfcpw PW basis for wave function if needed + * @param istep ionic step index + * @param out_dir output directory + * @param inp input parameters */ void init_scf(const UnitCell& ucell, const Parallel_Grid& pgrid, const ModuleBase::ComplexMatrix& strucfac, const bool* numeric, ModuleSymmetry::Symmetry& symm, - const void* wfcpw = nullptr); + const void* wfcpw = nullptr, + const int istep = 0, + const std::string& out_dir = "", + const Input_para& inp = PARAM.inp); std::string classname = "elecstate"; int iter = 0; ///< scf iteration diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp index 1073774b38c..1818d6fb438 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"); @@ -40,7 +41,7 @@ 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); + locpp.numeric, ucell.symm, (void*)pw_wfc, istep, out_dir, inp); //---------------------------------------------------------- //! 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 } From ab3c9215b8cc2a67c93e50e6fd3ba59581358296 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 29 May 2026 17:47:05 +0800 Subject: [PATCH 21/31] Remove unused parameters from ElecState::init_scf - Removed unused 'symm' and 'wfcpw' parameters from init_scf function - Updated all call sites to match the new signature - Simplified function interface by removing parameters not used in implementation --- source/source_esolver/esolver_double_xc.cpp | 2 +- source/source_esolver/esolver_ks_lcao.cpp | 2 +- source/source_esolver/esolver_ks_lcao_tddft.cpp | 2 +- source/source_esolver/esolver_of.cpp | 4 ++-- source/source_esolver/lcao_others.cpp | 2 +- source/source_estate/elecstate.cpp | 2 -- source/source_estate/elecstate.h | 2 -- source/source_estate/test/elecstate_base_test.cpp | 3 +-- source/source_hsolver/test/hsolver_supplementary_mock.h | 5 +++-- source/source_pw/module_pwdft/setup_pot.cpp | 2 +- 10 files changed, 11 insertions(+), 15 deletions(-) diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp index 14a86779edc..964dee667e5 100644 --- a/source/source_esolver/esolver_double_xc.cpp +++ b/source/source_esolver/esolver_double_xc.cpp @@ -158,7 +158,7 @@ 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); + this->pelec_base->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric); 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_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index a4931243ac9..caf8e7636cf 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -195,7 +195,7 @@ 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, + this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, istep, PARAM.globalv.global_out_dir, PARAM.inp); #ifdef __MLALGO diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index 033dbfa58d3..a1723e6723e 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -177,7 +177,7 @@ 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); + this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric); if (totstep <= PARAM.inp.td_tend + 1) { diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp index bf6f7ffde29..e697348d87c 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/source/source_esolver/esolver_of.cpp @@ -94,7 +94,7 @@ 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 + this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric); // atomic_rho, v_of_rho, set_vrs // liuyu move here 2023-10-09 // D in uspp need vloc, thus behind init_scf() @@ -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); + this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric); 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..629a680fc2c 100644 --- a/source/source_esolver/lcao_others.cpp +++ b/source/source_esolver/lcao_others.cpp @@ -172,7 +172,7 @@ 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); + this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric); // self consistent calculations for electronic ground state if (cal_type == "get_pchg") diff --git a/source/source_estate/elecstate.cpp b/source/source_estate/elecstate.cpp index f6dbeef455f..9b0d320ca45 100644 --- a/source/source_estate/elecstate.cpp +++ b/source/source_estate/elecstate.cpp @@ -29,8 +29,6 @@ void ElecState::init_scf(const UnitCell& ucell, const Parallel_Grid& pgrid, const ModuleBase::ComplexMatrix& strucfac, const bool* numeric, - ModuleSymmetry::Symmetry& symm, - const void* wfcpw, const int istep, const std::string& out_dir, const Input_para& inp) diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h index 49607f306fb..32904135b86 100644 --- a/source/source_estate/elecstate.h +++ b/source/source_estate/elecstate.h @@ -109,8 +109,6 @@ class ElecState const Parallel_Grid& pgrid, const ModuleBase::ComplexMatrix& strucfac, const bool* numeric, - ModuleSymmetry::Symmetry& symm, - const void* wfcpw = nullptr, const int istep = 0, const std::string& out_dir = "", const Input_para& inp = PARAM.inp); diff --git a/source/source_estate/test/elecstate_base_test.cpp b/source/source_estate/test/elecstate_base_test.cpp index e0d63df8220..4c3712d34dc 100644 --- a/source/source_estate/test/elecstate_base_test.cpp +++ b/source/source_estate/test/elecstate_base_test.cpp @@ -262,8 +262,7 @@ TEST_F(ElecStateTest, InitSCF) int istep = 0; ModuleBase::ComplexMatrix strucfac; elecstate->eferm = efermi; - ModuleSymmetry::Symmetry symm; - EXPECT_NO_THROW(elecstate->init_scf(ucell, pgrid, strucfac, nullptr, symm)); + EXPECT_NO_THROW(elecstate->init_scf(ucell, pgrid, strucfac, nullptr)); // delete elecstate->pot is done in the destructor of elecstate delete charge; } diff --git a/source/source_hsolver/test/hsolver_supplementary_mock.h b/source/source_hsolver/test/hsolver_supplementary_mock.h index 87155ecd5a7..d060ece0410 100644 --- a/source/source_hsolver/test/hsolver_supplementary_mock.h +++ b/source/source_hsolver/test/hsolver_supplementary_mock.h @@ -23,8 +23,9 @@ void ElecState::init_scf(const UnitCell& ucell, const Parallel_Grid& pgrid, const ModuleBase::ComplexMatrix& strucfac, const bool*, - ModuleSymmetry::Symmetry&, - const void*) + const int istep, + const std::string& out_dir, + const Input_para& inp) { return; } diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp index 1818d6fb438..899f21fcb65 100644 --- a/source/source_pw/module_pwdft/setup_pot.cpp +++ b/source/source_pw/module_pwdft/setup_pot.cpp @@ -41,7 +41,7 @@ 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, istep, out_dir, inp); + locpp.numeric, istep, out_dir, inp); //---------------------------------------------------------- //! 2) Symmetrize the charge density (rho) From 5a69d5cbd16baa6b70ff5eee051e553649aa50f7 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 29 May 2026 18:00:45 +0800 Subject: [PATCH 22/31] Fix missing io_basic library link in elecstate tests - Added io_basic library dependency to MODULE_ESTATE_elecstate_base test - Added io_basic library dependency to MODULE_ESTATE_elecstate_pw test - Fixes undefined reference to ModuleIO::write_pot_init --- source/source_estate/test/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/source_estate/test/CMakeLists.txt b/source/source_estate/test/CMakeLists.txt index aa69f4e29d8..27ab77690ba 100644 --- a/source/source_estate/test/CMakeLists.txt +++ b/source/source_estate/test/CMakeLists.txt @@ -43,14 +43,14 @@ AddTest( AddTest( TARGET MODULE_ESTATE_elecstate_base - LIBS parameter ${math_libs} base device + LIBS parameter ${math_libs} base device io_basic 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 ) AddTest( TARGET MODULE_ESTATE_elecstate_pw - LIBS parameter ${math_libs} planewave_serial base device + LIBS parameter ${math_libs} planewave_serial base device io_basic SOURCES elecstate_pw_test.cpp ../elecstate_pw.cpp ../elecstate_pw_cal_tau.cpp From 68350326278dea1bcf37824d0c45c043b2702d4d Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 29 May 2026 20:41:50 +0800 Subject: [PATCH 23/31] update init_scf --- source/source_esolver/esolver_double_xc.cpp | 2 +- source/source_esolver/esolver_ks_lcao_tddft.cpp | 3 ++- source/source_esolver/esolver_of.cpp | 7 +++---- source/source_esolver/lcao_others.cpp | 3 ++- source/source_estate/elecstate.h | 6 +++--- source/source_estate/test/elecstate_base_test.cpp | 3 ++- 6 files changed, 13 insertions(+), 11 deletions(-) diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp index 964dee667e5..5c7201904bf 100644 --- a/source/source_esolver/esolver_double_xc.cpp +++ b/source/source_esolver/esolver_double_xc.cpp @@ -158,7 +158,7 @@ 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); + this->pelec_base->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, 0, PARAM.globalv.global_out_dir, PARAM.inp); 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_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index a1723e6723e..0d151aee6b3 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); + this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, istep, + PARAM.globalv.global_out_dir, PARAM.inp); if (totstep <= PARAM.inp.td_tend + 1) { diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp index e697348d87c..f4ca985812c 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/source/source_esolver/esolver_of.cpp @@ -94,7 +94,8 @@ 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); // atomic_rho, v_of_rho, set_vrs + this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep, + PARAM.globalv.global_out_dir, PARAM.inp); // liuyu move here 2023-10-09 // D in uspp need vloc, thus behind init_scf() @@ -188,8 +189,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 +231,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) } } - this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric); + this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep, PARAM.globalv.global_out_dir, PARAM.inp); 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 629a680fc2c..cd6e253ecc5 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); + this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, + istep, PARAM.globalv.global_out_dir, PARAM.inp); // self consistent calculations for electronic ground state if (cal_type == "get_pchg") diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h index 32904135b86..39dd2fb9850 100644 --- a/source/source_estate/elecstate.h +++ b/source/source_estate/elecstate.h @@ -109,9 +109,9 @@ class ElecState const Parallel_Grid& pgrid, const ModuleBase::ComplexMatrix& strucfac, const bool* numeric, - const int istep = 0, - const std::string& out_dir = "", - const Input_para& inp = PARAM.inp); + const int istep, + const std::string& out_dir, + const Input_para& inp); std::string classname = "elecstate"; int iter = 0; ///< scf iteration diff --git a/source/source_estate/test/elecstate_base_test.cpp b/source/source_estate/test/elecstate_base_test.cpp index 4c3712d34dc..1024c77399e 100644 --- a/source/source_estate/test/elecstate_base_test.cpp +++ b/source/source_estate/test/elecstate_base_test.cpp @@ -262,7 +262,8 @@ TEST_F(ElecStateTest, InitSCF) int istep = 0; ModuleBase::ComplexMatrix strucfac; elecstate->eferm = efermi; - EXPECT_NO_THROW(elecstate->init_scf(ucell, pgrid, strucfac, nullptr)); + Input_para inp; + EXPECT_NO_THROW(elecstate->init_scf(ucell, pgrid, strucfac, nullptr, istep, PARAM.globalv.global_out_dir, inp)); // delete elecstate->pot is done in the destructor of elecstate delete charge; } From 04eb0f72c236e18d8aea2cf32b6709db02bc3fda Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 29 May 2026 20:43:43 +0800 Subject: [PATCH 24/31] fix --- source/source_esolver/esolver_double_xc.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp index 5c7201904bf..a62ed0cf804 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, 0, PARAM.globalv.global_out_dir, PARAM.inp); + this->pelec_base->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, istep, + PARAM.globalv.global_out_dir, PARAM.inp); XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); // DMR should be same size with Hamiltonian(R) From 6f3058257db20e645313f2a8c0d01aa0d4ba329c Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 29 May 2026 21:36:17 +0800 Subject: [PATCH 25/31] fix bug --- source/source_esolver/esolver_of.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp index f4ca985812c..c758d6991d5 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/source/source_esolver/esolver_of.cpp @@ -94,6 +94,7 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp) this->init_elecstate(ucell); // calculate the total local pseudopotential in real space + const int istep=0; this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep, PARAM.globalv.global_out_dir, PARAM.inp); From 84637598be9a160db84f1053b47d454706a920e1 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 29 May 2026 22:14:26 +0800 Subject: [PATCH 26/31] remove dependence of parameter for write_cube.cpp --- source/source_esolver/esolver_dm2rho.cpp | 3 ++- source/source_estate/test/CMakeLists.txt | 4 ++-- source/source_io/module_chgpot/get_pchg_lcao.cpp | 7 ++++--- source/source_io/module_chgpot/get_pchg_pw.h | 3 ++- .../module_chgpot/write_elecstat_pot.cpp | 3 ++- source/source_io/module_chgpot/write_init.cpp | 7 +++++-- source/source_io/module_ctrl/ctrl_output_fp.cpp | 14 ++++++++++---- source/source_io/module_dos/cal_ldos.cpp | 5 +++-- source/source_io/module_elf/write_elf.cpp | 12 ++++++++---- source/source_io/module_elf/write_elf.h | 3 ++- source/source_io/module_output/cube_io.h | 1 + source/source_io/module_output/write_cube.cpp | 4 ++-- source/source_io/module_wf/get_wf_lcao.cpp | 16 ++++++++++------ source/source_io/module_wf/get_wf_pw.h | 3 +++ source/source_io/test_serial/rho_io_test.cpp | 2 +- 15 files changed, 57 insertions(+), 30 deletions(-) diff --git a/source/source_esolver/esolver_dm2rho.cpp b/source/source_esolver/esolver_dm2rho.cpp index 0903afb29a9..8265f5493ea 100644 --- a/source/source_esolver/esolver_dm2rho.cpp +++ b/source/source_esolver/esolver_dm2rho.cpp @@ -79,7 +79,8 @@ void ESolver_DM2rho::runner(UnitCell& ucell, const int istep) this->pelec->eferm.get_efval(is), &(ucell), 3, - 1); + 1, + PARAM.sys.two_fermi); } ModuleBase::timer::end("ESolver_DM2rho", "runner"); diff --git a/source/source_estate/test/CMakeLists.txt b/source/source_estate/test/CMakeLists.txt index 27ab77690ba..a64989f88f7 100644 --- a/source/source_estate/test/CMakeLists.txt +++ b/source/source_estate/test/CMakeLists.txt @@ -43,14 +43,14 @@ AddTest( AddTest( TARGET MODULE_ESTATE_elecstate_base - LIBS parameter ${math_libs} base device io_basic + 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 ) AddTest( TARGET MODULE_ESTATE_elecstate_pw - LIBS parameter ${math_libs} planewave_serial base device io_basic + LIBS parameter ${math_libs} planewave_serial base device SOURCES elecstate_pw_test.cpp ../elecstate_pw.cpp ../elecstate_pw_cal_tau.cpp diff --git a/source/source_io/module_chgpot/get_pchg_lcao.cpp b/source/source_io/module_chgpot/get_pchg_lcao.cpp index 0c1b13d3b3d..3bf22cfe208 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.sys.two_fermi); } } } @@ -192,7 +193,7 @@ void Get_pchg_lcao::begin(double** rho, ssc.str(), ef_spin, ucell_in, - precision); + precision, 1, PARAM.sys.two_fermi); } } } @@ -244,7 +245,7 @@ void Get_pchg_lcao::begin(double** rho, ssc.str(), ef_spin, ucell_in, - precision); + precision, 1, PARAM.sys.two_fermi); } } } diff --git a/source/source_io/module_chgpot/get_pchg_pw.h b/source/source_io/module_chgpot/get_pchg_pw.h index 7609b30eb8a..1ca71fe4419 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.sys.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.sys.two_fermi); } } // 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..5bf026cd9bd 100644 --- a/source/source_io/module_chgpot/write_elecstat_pot.cpp +++ b/source/source_io/module_chgpot/write_elecstat_pot.cpp @@ -102,7 +102,8 @@ void write_elecstat_pot( ef_tmp, ucell, precision, - out_fermi); + out_fermi, + PARAM.sys.two_fermi); 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 b455c3833dd..1eaea729be8 100644 --- a/source/source_io/module_chgpot/write_init.cpp +++ b/source/source_io/module_chgpot/write_init.cpp @@ -111,7 +111,9 @@ void ModuleIO::write_chg_init( filename, fermi_energy, &(ucell), - inp.out_chg[1]); + inp.out_chg[1], + 1, + inp.two_fermi); } } } @@ -157,7 +159,8 @@ void ModuleIO::write_pot_init( 0.0, &(ucell), inp.out_pot[1], - 0); + 0, + inp.two_fermi); } } } diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index a36dae61750..2fbe1e5dcb9 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -83,7 +83,8 @@ void ctrl_output_fp(UnitCell& ucell, pelec->eferm.get_efval(is), &(ucell), PARAM.inp.out_chg[1], - 1); + 1, + PARAM.sys.two_fermi); if (XC_Functional::get_ked_flag()) { @@ -98,7 +99,10 @@ void ctrl_output_fp(UnitCell& ucell, istep, fn, pelec->eferm.get_efval(is), - &(ucell)); + &(ucell), + 11, // default precision + 1, // default out_fermi + PARAM.sys.two_fermi); } } } @@ -131,7 +135,8 @@ void ctrl_output_fp(UnitCell& ucell, 0.0, // efermi &(ucell), PARAM.inp.out_pot[1], // precision - 0); // out_fermi + 0, // out_fermi + PARAM.sys.two_fermi); } } else if (PARAM.inp.out_pot[0] == 2 && should_output) @@ -174,7 +179,8 @@ void ctrl_output_fp(UnitCell& ucell, para_grid, &(ucell), PARAM.inp.out_elf[1], - geom_block); + geom_block, + PARAM.sys.two_fermi); } #ifdef USE_LIBXC diff --git a/source/source_io/module_dos/cal_ldos.cpp b/source/source_io/module_dos/cal_ldos.cpp index 4463c4cf036..80aa6bcf7d0 100644 --- a/source/source_io/module_dos/cal_ldos.cpp +++ b/source/source_io/module_dos/cal_ldos.cpp @@ -90,7 +90,8 @@ void Cal_ldos::cal_ldos_lcao( 0, &ucell, precision, - 0); + 0, + PARAM.sys.two_fermi); // free memory delete[] ldos; @@ -162,7 +163,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.sys.two_fermi); } } diff --git a/source/source_io/module_elf/write_elf.cpp b/source/source_io/module_elf/write_elf.cpp index daf830db623..7b2c23bfb74 100644 --- a/source/source_io/module_elf/write_elf.cpp +++ b/source/source_io/module_elf/write_elf.cpp @@ -16,7 +16,8 @@ void write_elf( const Parallel_Grid& pgrid, const UnitCell* ucell_, const int& precision, - const std::string& geom_block) + const std::string& geom_block, + const bool two_fermi) { ModuleBase::timer::start("ModuleIO", "write_elf"); // For nspin = 4, we only calculate the total ELF using the @@ -170,7 +171,8 @@ void write_elf( ef_tmp, ucell_, precision, - out_fermi); + out_fermi, + two_fermi); } else if (nspin == 2) { @@ -190,7 +192,8 @@ void write_elf( ef_tmp, ucell_, precision, - out_fermi); + out_fermi, + two_fermi); } std::vector elf_tot(nrxx, 0.0); @@ -224,7 +227,8 @@ void write_elf( ef_tmp, ucell_, precision, - out_fermi); + out_fermi, + two_fermi); } 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 7661a379b11..ea5f1b788e7 100644 --- a/source/source_io/module_elf/write_elf.h +++ b/source/source_io/module_elf/write_elf.h @@ -17,7 +17,8 @@ void write_elf( const Parallel_Grid& pgrid, const UnitCell* ucell_, const int& precision, - const std::string& geom_block); + 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..23f369e210f 100644 --- a/source/source_io/module_output/cube_io.h +++ b/source/source_io/module_output/cube_io.h @@ -26,6 +26,7 @@ void write_vdata_palgrid(const Parallel_Grid& pgrid, const UnitCell* const ucell, const int precision = 11, const int out_fermi = 1, + const bool two_fermi, // whether to use two fermi energies const bool reduce_all_pool = false); // only reduce in the main pool as default /// read the full data from a cube 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_wf/get_wf_lcao.cpp b/source/source_io/module_wf/get_wf_lcao.cpp index cadf8d0edc0..8930f315b34 100644 --- a/source/source_io/module_wf/get_wf_lcao.cpp +++ b/source/source_io/module_wf/get_wf_lcao.cpp @@ -83,7 +83,10 @@ void Get_wf_lcao::begin(const UnitCell& ucell, 0, ss_out.str(), ef_tmp, - &(ucell)); + &(ucell), + 11, // default precision + 1, // default out_fermi + PARAM.sys.two_fermi); } } } @@ -126,12 +129,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.sys.two_fermi); // 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.sys.two_fermi); } } } @@ -245,7 +248,8 @@ void Get_wf_lcao::begin(const UnitCell& ucell, ef_tmp, &(ucell), 3, - 1); + 1, + PARAM.sys.two_fermi); // if (out_wfc_pw || out_wf_r) psi_g.fix_k(ik); @@ -316,12 +320,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.sys.two_fermi); // 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.sys.two_fermi); } } } diff --git a/source/source_io/module_wf/get_wf_pw.h b/source/source_io/module_wf/get_wf_pw.h index 9262ad181b6..f3f92de651b 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.sys.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.sys.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.sys.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..491b4ccdec7 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); 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()); From 92c9e0870747a17bd89f2d5c3f12f2b003627e4f Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 29 May 2026 22:20:34 +0800 Subject: [PATCH 27/31] fix bugs --- source/source_esolver/esolver_dm2rho.cpp | 2 +- source/source_estate/test/CMakeLists.txt | 6 ++++-- source/source_estate/test/elecstate_base_test.cpp | 2 +- source/source_io/module_chgpot/get_pchg_lcao.cpp | 6 +++--- source/source_io/module_chgpot/get_pchg_pw.h | 4 ++-- .../source_io/module_chgpot/write_elecstat_pot.cpp | 2 +- source/source_io/module_ctrl/ctrl_output_fp.cpp | 8 ++++---- source/source_io/module_dos/cal_ldos.cpp | 4 ++-- source/source_io/module_wf/get_wf_lcao.cpp | 12 ++++++------ source/source_io/module_wf/get_wf_pw.h | 6 +++--- 10 files changed, 27 insertions(+), 25 deletions(-) diff --git a/source/source_esolver/esolver_dm2rho.cpp b/source/source_esolver/esolver_dm2rho.cpp index 8265f5493ea..564be1201af 100644 --- a/source/source_esolver/esolver_dm2rho.cpp +++ b/source/source_esolver/esolver_dm2rho.cpp @@ -80,7 +80,7 @@ void ESolver_DM2rho::runner(UnitCell& ucell, const int istep) &(ucell), 3, 1, - PARAM.sys.two_fermi); + PARAM.globalv.two_fermi); } ModuleBase::timer::end("ESolver_DM2rho", "runner"); diff --git a/source/source_estate/test/CMakeLists.txt b/source/source_estate/test/CMakeLists.txt index a64989f88f7..16b14549505 100644 --- a/source/source_estate/test/CMakeLists.txt +++ b/source/source_estate/test/CMakeLists.txt @@ -43,14 +43,15 @@ AddTest( AddTest( TARGET MODULE_ESTATE_elecstate_base - LIBS parameter ${math_libs} base device + LIBS parameter ${math_libs} base device cell_info 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 + ../../source_io/module_chgpot/write_init.cpp ../../source_io/module_output/write_cube.cpp ) AddTest( TARGET MODULE_ESTATE_elecstate_pw - LIBS parameter ${math_libs} planewave_serial base device + LIBS parameter ${math_libs} planewave_serial base device cell_info SOURCES elecstate_pw_test.cpp ../elecstate_pw.cpp ../elecstate_pw_cal_tau.cpp @@ -60,6 +61,7 @@ AddTest( ../../source_psi/psi.cpp # ../../source_psi/kernels/psi_memory_op.cpp ../../source_base/module_device/memory_op.cpp + ../../source_io/module_chgpot/write_init.cpp ../../source_io/module_output/write_cube.cpp ) AddTest( diff --git a/source/source_estate/test/elecstate_base_test.cpp b/source/source_estate/test/elecstate_base_test.cpp index 1024c77399e..c40ce87f57f 100644 --- a/source/source_estate/test/elecstate_base_test.cpp +++ b/source/source_estate/test/elecstate_base_test.cpp @@ -187,7 +187,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; diff --git a/source/source_io/module_chgpot/get_pchg_lcao.cpp b/source/source_io/module_chgpot/get_pchg_lcao.cpp index 3bf22cfe208..e3faf934ce9 100644 --- a/source/source_io/module_chgpot/get_pchg_lcao.cpp +++ b/source/source_io/module_chgpot/get_pchg_lcao.cpp @@ -96,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, 1, PARAM.sys.two_fermi); + ssc.str(), ef_spin, ucell_in, precision, 1, PARAM.globalv.two_fermi); } } } @@ -193,7 +193,7 @@ void Get_pchg_lcao::begin(double** rho, ssc.str(), ef_spin, ucell_in, - precision, 1, PARAM.sys.two_fermi); + precision, 1, PARAM.globalv.two_fermi); } } } @@ -245,7 +245,7 @@ void Get_pchg_lcao::begin(double** rho, ssc.str(), ef_spin, ucell_in, - precision, 1, PARAM.sys.two_fermi); + precision, 1, PARAM.globalv.two_fermi); } } } diff --git a/source/source_io/module_chgpot/get_pchg_pw.h b/source/source_io/module_chgpot/get_pchg_pw.h index 1ca71fe4419..339cb179d98 100644 --- a/source/source_io/module_chgpot/get_pchg_pw.h +++ b/source/source_io/module_chgpot/get_pchg_pw.h @@ -153,7 +153,7 @@ void get_pchg_pw(const std::vector& out_pchg, ucell, 11, 1, - PARAM.sys.two_fermi, + PARAM.globalv.two_fermi, true); // reduce_all_pool is true } } @@ -236,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, 11, 1, PARAM.sys.two_fermi); + ModuleIO::write_vdata_palgrid(pgrid, rho_band[is].data(), is, nspin, 0, ssc.str(), 0.0, ucell, 11, 1, PARAM.globalv.two_fermi); } } // 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 5bf026cd9bd..4d55459271c 100644 --- a/source/source_io/module_chgpot/write_elecstat_pot.cpp +++ b/source/source_io/module_chgpot/write_elecstat_pot.cpp @@ -103,7 +103,7 @@ void write_elecstat_pot( ucell, precision, out_fermi, - PARAM.sys.two_fermi); + PARAM.globalv.two_fermi); ModuleBase::timer::end("ModuleIO", "write_elecstat_pot"); return; diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index 2fbe1e5dcb9..5e2eb8cf63d 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -84,7 +84,7 @@ void ctrl_output_fp(UnitCell& ucell, &(ucell), PARAM.inp.out_chg[1], 1, - PARAM.sys.two_fermi); + PARAM.globalv.two_fermi); if (XC_Functional::get_ked_flag()) { @@ -102,7 +102,7 @@ void ctrl_output_fp(UnitCell& ucell, &(ucell), 11, // default precision 1, // default out_fermi - PARAM.sys.two_fermi); + PARAM.globalv.two_fermi); } } } @@ -136,7 +136,7 @@ void ctrl_output_fp(UnitCell& ucell, &(ucell), PARAM.inp.out_pot[1], // precision 0, // out_fermi - PARAM.sys.two_fermi); + PARAM.globalv.two_fermi); } } else if (PARAM.inp.out_pot[0] == 2 && should_output) @@ -180,7 +180,7 @@ void ctrl_output_fp(UnitCell& ucell, &(ucell), PARAM.inp.out_elf[1], geom_block, - PARAM.sys.two_fermi); + PARAM.globalv.two_fermi); } #ifdef USE_LIBXC diff --git a/source/source_io/module_dos/cal_ldos.cpp b/source/source_io/module_dos/cal_ldos.cpp index 80aa6bcf7d0..400664a08c9 100644 --- a/source/source_io/module_dos/cal_ldos.cpp +++ b/source/source_io/module_dos/cal_ldos.cpp @@ -91,7 +91,7 @@ void Cal_ldos::cal_ldos_lcao( &ucell, precision, 0, - PARAM.sys.two_fermi); + PARAM.globalv.two_fermi); // free memory delete[] ldos; @@ -163,7 +163,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, PARAM.sys.two_fermi); + ModuleIO::write_vdata_palgrid(pgrid, ldos.data(), 0, PARAM.inp.nspin, 0, fn.str(), 0, &ucell, precision, 0, PARAM.globalv.two_fermi); } } diff --git a/source/source_io/module_wf/get_wf_lcao.cpp b/source/source_io/module_wf/get_wf_lcao.cpp index 8930f315b34..cc3bfd7f104 100644 --- a/source/source_io/module_wf/get_wf_lcao.cpp +++ b/source/source_io/module_wf/get_wf_lcao.cpp @@ -86,7 +86,7 @@ void Get_wf_lcao::begin(const UnitCell& ucell, &(ucell), 11, // default precision 1, // default out_fermi - PARAM.sys.two_fermi); + PARAM.globalv.two_fermi); } } } @@ -129,12 +129,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), 11, 1, PARAM.sys.two_fermi); + ModuleIO::write_vdata_palgrid(pgrid, wfc_real.data(), is, nspin, 0, ss_real.str(), ef_tmp, &(ucell), 11, 1, PARAM.globalv.two_fermi); // 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), 11, 1, PARAM.sys.two_fermi); + ModuleIO::write_vdata_palgrid(pgrid, wfc_imag.data(), is, nspin, 0, ss_imag.str(), ef_tmp, &(ucell), 11, 1, PARAM.globalv.two_fermi); } } } @@ -249,7 +249,7 @@ void Get_wf_lcao::begin(const UnitCell& ucell, &(ucell), 3, 1, - PARAM.sys.two_fermi); + PARAM.globalv.two_fermi); // if (out_wfc_pw || out_wf_r) psi_g.fix_k(ik); @@ -320,12 +320,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), 11, 1, PARAM.sys.two_fermi); + ModuleIO::write_vdata_palgrid(pgrid, wfc_real.data(), ispin, nspin, 0, ss_real.str(), ef_tmp, &(ucell), 11, 1, PARAM.globalv.two_fermi); // 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), 11, 1, PARAM.sys.two_fermi); + ModuleIO::write_vdata_palgrid(pgrid, wfc_imag.data(), ispin, nspin, 0, ss_imag.str(), ef_tmp, &(ucell), 11, 1, PARAM.globalv.two_fermi); } } } diff --git a/source/source_io/module_wf/get_wf_pw.h b/source/source_io/module_wf/get_wf_pw.h index f3f92de651b..0b9461ea68b 100644 --- a/source/source_io/module_wf/get_wf_pw.h +++ b/source/source_io/module_wf/get_wf_pw.h @@ -160,7 +160,7 @@ void get_wf_pw(const std::vector& out_wfc_norm, ucell, 11, 1, - PARAM.sys.two_fermi, + PARAM.globalv.two_fermi, true); // reduce_all_pool is true } } @@ -251,7 +251,7 @@ void get_wf_pw(const std::vector& out_wfc_norm, ucell, 11, 1, - PARAM.sys.two_fermi, + PARAM.globalv.two_fermi, true); // reduce_all_pool is true std::stringstream ss_imag; @@ -267,7 +267,7 @@ void get_wf_pw(const std::vector& out_wfc_norm, ucell, 11, 1, - PARAM.sys.two_fermi, + PARAM.globalv.two_fermi, true); // reduce_all_pool is true } } From a1b5a391dc9244caefd5359169aded59d78381a8 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 29 May 2026 22:33:48 +0800 Subject: [PATCH 28/31] fix bug --- source/source_esolver/esolver_dm2rho.cpp | 3 ++- source/source_esolver/esolver_fp.cpp | 2 +- source/source_estate/elecstate.cpp | 2 +- source/source_io/module_chgpot/get_pchg_lcao.cpp | 6 +++--- source/source_io/module_chgpot/get_pchg_pw.h | 2 +- .../source_io/module_chgpot/write_elecstat_pot.cpp | 3 ++- source/source_io/module_chgpot/write_init.cpp | 12 ++++++++---- source/source_io/module_chgpot/write_init.h | 6 ++++-- source/source_io/module_ctrl/ctrl_output_fp.cpp | 9 ++++++--- source/source_io/module_dos/cal_ldos.cpp | 5 +++-- source/source_io/module_elf/write_elf.cpp | 9 ++++++--- source/source_io/module_output/cube_io.h | 8 ++++---- source/source_io/module_wf/get_wf_lcao.cpp | 14 ++++++++------ source/source_io/test_serial/rho_io_test.cpp | 2 +- 14 files changed, 50 insertions(+), 33 deletions(-) diff --git a/source/source_esolver/esolver_dm2rho.cpp b/source/source_esolver/esolver_dm2rho.cpp index 564be1201af..1548d748321 100644 --- a/source/source_esolver/esolver_dm2rho.cpp +++ b/source/source_esolver/esolver_dm2rho.cpp @@ -80,7 +80,8 @@ void ESolver_DM2rho::runner(UnitCell& ucell, const int istep) &(ucell), 3, 1, - PARAM.globalv.two_fermi); + PARAM.globalv.two_fermi, + false); } ModuleBase::timer::end("ESolver_DM2rho", "runner"); diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 6a13a4a5095..1ddb636c8ab 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -177,7 +177,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, PARAM.inp); + PARAM.globalv.global_out_dir, PARAM.inp, PARAM.globalv.two_fermi); return; } diff --git a/source/source_estate/elecstate.cpp b/source/source_estate/elecstate.cpp index 9b0d320ca45..c06695cd3a5 100644 --- a/source/source_estate/elecstate.cpp +++ b/source/source_estate/elecstate.cpp @@ -43,7 +43,7 @@ void ElecState::init_scf(const UnitCell& ucell, this->pot->init_pot(this->charge); //! output the initial potential - ModuleIO::write_pot_init(ucell, pgrid, this, istep, out_dir, inp); + ModuleIO::write_pot_init(ucell, pgrid, this, istep, out_dir, inp, PARAM.globalv.two_fermi); } diff --git a/source/source_io/module_chgpot/get_pchg_lcao.cpp b/source/source_io/module_chgpot/get_pchg_lcao.cpp index e3faf934ce9..e0559c4e5f1 100644 --- a/source/source_io/module_chgpot/get_pchg_lcao.cpp +++ b/source/source_io/module_chgpot/get_pchg_lcao.cpp @@ -96,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, 1, PARAM.globalv.two_fermi); + ssc.str(), ef_spin, ucell_in, precision, 1, PARAM.globalv.two_fermi, false); } } } @@ -193,7 +193,7 @@ void Get_pchg_lcao::begin(double** rho, ssc.str(), ef_spin, ucell_in, - precision, 1, PARAM.globalv.two_fermi); + precision, 1, PARAM.globalv.two_fermi, false); } } } @@ -245,7 +245,7 @@ void Get_pchg_lcao::begin(double** rho, ssc.str(), ef_spin, ucell_in, - precision, 1, PARAM.globalv.two_fermi); + 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 339cb179d98..cee654aac83 100644 --- a/source/source_io/module_chgpot/get_pchg_pw.h +++ b/source/source_io/module_chgpot/get_pchg_pw.h @@ -236,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, 11, 1, PARAM.globalv.two_fermi); + 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 4d55459271c..da6c5b04df4 100644 --- a/source/source_io/module_chgpot/write_elecstat_pot.cpp +++ b/source/source_io/module_chgpot/write_elecstat_pot.cpp @@ -103,7 +103,8 @@ void write_elecstat_pot( ucell, precision, out_fermi, - PARAM.globalv.two_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 1eaea729be8..b7ffdc18893 100644 --- a/source/source_io/module_chgpot/write_init.cpp +++ b/source/source_io/module_chgpot/write_init.cpp @@ -63,7 +63,8 @@ void ModuleIO::write_chg_init( const elecstate::Efermi &efermi, const int istep, const std::string& out_dir, - const Input_para& inp) + const Input_para& inp, + const bool two_fermi) { const int nspin = inp.nspin; assert(nspin == 1 || nspin == 2 || nspin == 4); @@ -113,7 +114,8 @@ void ModuleIO::write_chg_init( &(ucell), inp.out_chg[1], 1, - inp.two_fermi); + two_fermi, + false); } } } @@ -127,7 +129,8 @@ void ModuleIO::write_pot_init( elecstate::ElecState *pelec, const int istep, const std::string& out_dir, - const Input_para& inp) + const Input_para& inp, + const bool two_fermi) { const int nspin = inp.nspin; assert(nspin == 1 || nspin == 2 || nspin == 4); @@ -160,7 +163,8 @@ void ModuleIO::write_pot_init( &(ucell), inp.out_pot[1], 0, - inp.two_fermi); + two_fermi, + false); } } } diff --git a/source/source_io/module_chgpot/write_init.h b/source/source_io/module_chgpot/write_init.h index 7d5ecd5239b..303f97e42ed 100644 --- a/source/source_io/module_chgpot/write_init.h +++ b/source/source_io/module_chgpot/write_init.h @@ -52,7 +52,8 @@ void write_chg_init( const elecstate::Efermi &efermi, const int istep, const std::string& out_dir, - const Input_para& inp); + 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. @@ -69,7 +70,8 @@ void write_pot_init( elecstate::ElecState *pelec, const int istep, const std::string& out_dir, - const Input_para& inp); + const Input_para& inp, + const bool two_fermi); } diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index 5e2eb8cf63d..87aa062bc0c 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -84,7 +84,8 @@ void ctrl_output_fp(UnitCell& ucell, &(ucell), PARAM.inp.out_chg[1], 1, - PARAM.globalv.two_fermi); + PARAM.globalv.two_fermi, + false); if (XC_Functional::get_ked_flag()) { @@ -102,7 +103,8 @@ void ctrl_output_fp(UnitCell& ucell, &(ucell), 11, // default precision 1, // default out_fermi - PARAM.globalv.two_fermi); + PARAM.globalv.two_fermi, + false); } } } @@ -136,7 +138,8 @@ void ctrl_output_fp(UnitCell& ucell, &(ucell), PARAM.inp.out_pot[1], // precision 0, // out_fermi - PARAM.globalv.two_fermi); + PARAM.globalv.two_fermi, + false); } } else if (PARAM.inp.out_pot[0] == 2 && should_output) diff --git a/source/source_io/module_dos/cal_ldos.cpp b/source/source_io/module_dos/cal_ldos.cpp index 400664a08c9..1e92e2e4e4b 100644 --- a/source/source_io/module_dos/cal_ldos.cpp +++ b/source/source_io/module_dos/cal_ldos.cpp @@ -91,7 +91,8 @@ void Cal_ldos::cal_ldos_lcao( &ucell, precision, 0, - PARAM.globalv.two_fermi); + PARAM.globalv.two_fermi, + false); // free memory delete[] ldos; @@ -163,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, PARAM.globalv.two_fermi); + 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 7b2c23bfb74..ea6e9b88050 100644 --- a/source/source_io/module_elf/write_elf.cpp +++ b/source/source_io/module_elf/write_elf.cpp @@ -172,7 +172,8 @@ void write_elf( ucell_, precision, out_fermi, - two_fermi); + two_fermi, + false); } else if (nspin == 2) { @@ -193,7 +194,8 @@ void write_elf( ucell_, precision, out_fermi, - two_fermi); + two_fermi, + false); } std::vector elf_tot(nrxx, 0.0); @@ -228,7 +230,8 @@ void write_elf( ucell_, precision, out_fermi, - two_fermi); + two_fermi, + false); } ModuleBase::timer::end("ModuleIO", "write_elf"); } // end write_elf diff --git a/source/source_io/module_output/cube_io.h b/source/source_io/module_output/cube_io.h index 23f369e210f..f97a639b1d0 100644 --- a/source/source_io/module_output/cube_io.h +++ b/source/source_io/module_output/cube_io.h @@ -24,10 +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 two_fermi, // whether to use two fermi energies - 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_wf/get_wf_lcao.cpp b/source/source_io/module_wf/get_wf_lcao.cpp index cc3bfd7f104..cd54b4cbdd5 100644 --- a/source/source_io/module_wf/get_wf_lcao.cpp +++ b/source/source_io/module_wf/get_wf_lcao.cpp @@ -86,7 +86,8 @@ void Get_wf_lcao::begin(const UnitCell& ucell, &(ucell), 11, // default precision 1, // default out_fermi - PARAM.globalv.two_fermi); + PARAM.globalv.two_fermi, + false); } } } @@ -129,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), 11, 1, PARAM.globalv.two_fermi); + 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), 11, 1, PARAM.globalv.two_fermi); + ModuleIO::write_vdata_palgrid(pgrid, wfc_imag.data(), is, nspin, 0, ss_imag.str(), ef_tmp, &(ucell), 11, 1, PARAM.globalv.two_fermi, false); } } } @@ -249,7 +250,8 @@ void Get_wf_lcao::begin(const UnitCell& ucell, &(ucell), 3, 1, - PARAM.globalv.two_fermi); + PARAM.globalv.two_fermi, + false); // if (out_wfc_pw || out_wf_r) psi_g.fix_k(ik); @@ -320,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), 11, 1, PARAM.globalv.two_fermi); + 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), 11, 1, PARAM.globalv.two_fermi); + 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/test_serial/rho_io_test.cpp b/source/source_io/test_serial/rho_io_test.cpp index 491b4ccdec7..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, false); + 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()); From 1e00db80bb40229717cc599f879c26b785514a1e Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 29 May 2026 23:15:39 +0800 Subject: [PATCH 29/31] add a new file init_scf --- source/source_esolver/esolver_double_xc.cpp | 4 +- source/source_esolver/esolver_ks_lcao.cpp | 4 +- .../source_esolver/esolver_ks_lcao_tddft.cpp | 4 +- source/source_esolver/esolver_of.cpp | 6 +-- source/source_esolver/lcao_others.cpp | 4 +- source/source_estate/CMakeLists.txt | 1 + source/source_estate/elecstate.cpp | 19 --------- source/source_estate/elecstate.h | 41 +++++++++++-------- source/source_estate/init_scf.cpp | 29 +++++++++++++ source/source_estate/test/CMakeLists.txt | 1 - source/source_pw/module_pwdft/setup_pot.cpp | 4 +- 11 files changed, 66 insertions(+), 51 deletions(-) create mode 100644 source/source_estate/init_scf.cpp diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp index a62ed0cf804..6a7cb86f679 100644 --- a/source/source_esolver/esolver_double_xc.cpp +++ b/source/source_esolver/esolver_double_xc.cpp @@ -158,8 +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, istep, - PARAM.globalv.global_out_dir, PARAM.inp); + 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_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index caf8e7636cf..08499df2bdf 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -195,8 +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, - istep, PARAM.globalv.global_out_dir, PARAM.inp); + 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 0d151aee6b3..66d7717f14f 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -177,8 +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, istep, - PARAM.globalv.global_out_dir, PARAM.inp); + 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_of.cpp b/source/source_esolver/esolver_of.cpp index c758d6991d5..b346df60fb8 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/source/source_esolver/esolver_of.cpp @@ -95,8 +95,8 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp) // calculate the total local pseudopotential in real space const int istep=0; - this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep, - PARAM.globalv.global_out_dir, PARAM.inp); + 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() @@ -232,7 +232,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) } } - this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep, PARAM.globalv.global_out_dir, PARAM.inp); + 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 cd6e253ecc5..57174c47ceb 100644 --- a/source/source_esolver/lcao_others.cpp +++ b/source/source_esolver/lcao_others.cpp @@ -172,8 +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, - istep, PARAM.globalv.global_out_dir, PARAM.inp); + 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 c06695cd3a5..1865396d83a 100644 --- a/source/source_estate/elecstate.cpp +++ b/source/source_estate/elecstate.cpp @@ -25,26 +25,7 @@ void ElecState::init_nelec_spin() } } -void ElecState::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) -{ - //! 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); - - //! output the initial potential - ModuleIO::write_pot_init(ucell, pgrid, this, istep, out_dir, inp, PARAM.globalv.two_fermi); -} void ElecState::init_ks(Charge* chr_in, // pointer for class Charge diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h index 39dd2fb9850..a4191df38d8 100644 --- a/source/source_estate/elecstate.h +++ b/source/source_estate/elecstate.h @@ -94,24 +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 - * @param istep ionic step index - * @param out_dir output directory - * @param inp input parameters - */ - 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); + + std::string classname = "elecstate"; int iter = 0; ///< scf iteration @@ -166,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 16b14549505..50ea26af79b 100644 --- a/source/source_estate/test/CMakeLists.txt +++ b/source/source_estate/test/CMakeLists.txt @@ -46,7 +46,6 @@ AddTest( LIBS parameter ${math_libs} base device cell_info 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 - ../../source_io/module_chgpot/write_init.cpp ../../source_io/module_output/write_cube.cpp ) AddTest( diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp index 899f21fcb65..6a4697ebaf9 100644 --- a/source/source_pw/module_pwdft/setup_pot.cpp +++ b/source/source_pw/module_pwdft/setup_pot.cpp @@ -40,8 +40,8 @@ void pw::setup_pot(const int istep, //---------------------------------------------------------- //! 1) Renew local pseudopotential //---------------------------------------------------------- - pelec->init_scf(ucell, para_grid, sf.strucFac, - locpp.numeric, istep, out_dir, inp); + elecstate::init_scf(ucell, para_grid, sf.strucFac, + locpp.numeric, istep, out_dir, inp, pelec); //---------------------------------------------------------- //! 2) Symmetrize the charge density (rho) From 09feb49cfe0903c8a114679c3090826c71d2b79e Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 30 May 2026 05:42:29 +0800 Subject: [PATCH 30/31] update estate tests --- source/source_estate/test/CMakeLists.txt | 10 ++++------ .../source_estate/test/elecstate_base_test.cpp | 18 ------------------ .../test/hsolver_supplementary_mock.h | 11 ----------- 3 files changed, 4 insertions(+), 35 deletions(-) diff --git a/source/source_estate/test/CMakeLists.txt b/source/source_estate/test/CMakeLists.txt index 50ea26af79b..9bd76ae2486 100644 --- a/source/source_estate/test/CMakeLists.txt +++ b/source/source_estate/test/CMakeLists.txt @@ -43,14 +43,14 @@ AddTest( AddTest( TARGET MODULE_ESTATE_elecstate_base - LIBS parameter ${math_libs} base device cell_info + 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 ) AddTest( TARGET MODULE_ESTATE_elecstate_pw - LIBS parameter ${math_libs} planewave_serial base device cell_info + LIBS parameter ${math_libs} planewave_serial base device SOURCES elecstate_pw_test.cpp ../elecstate_pw.cpp ../elecstate_pw_cal_tau.cpp @@ -58,9 +58,7 @@ 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 - ../../source_io/module_chgpot/write_init.cpp ../../source_io/module_output/write_cube.cpp ) AddTest( @@ -89,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 @@ -98,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 c40ce87f57f..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() @@ -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; - Input_para inp; - EXPECT_NO_THROW(elecstate->init_scf(ucell, pgrid, strucfac, nullptr, istep, PARAM.globalv.global_out_dir, inp)); - // 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 d060ece0410..7dc30d77233 100644 --- a/source/source_hsolver/test/hsolver_supplementary_mock.h +++ b/source/source_hsolver/test/hsolver_supplementary_mock.h @@ -19,17 +19,6 @@ void ElecState::init_nelec_spin() -void ElecState::init_scf(const UnitCell& ucell, - const Parallel_Grid& pgrid, - const ModuleBase::ComplexMatrix& strucfac, - const bool*, - const int istep, - const std::string& out_dir, - const Input_para& inp) -{ - return; -} - void ElecState::init_ks(Charge* chg_in, // pointer for class Charge const K_Vectors* klist_in, int nk_in, From cf38047a3a63afbe7546197906dffae458b6f8ab Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 30 May 2026 05:48:00 +0800 Subject: [PATCH 31/31] delete useless inclusion --- source/source_estate/elecstate.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/source/source_estate/elecstate.cpp b/source/source_estate/elecstate.cpp index 1865396d83a..5ed97e678c7 100644 --- a/source/source_estate/elecstate.cpp +++ b/source/source_estate/elecstate.cpp @@ -3,7 +3,6 @@ #include "source_base/parallel_reduce.h" #include "source_base/tool_title.h" #include "occupy.h" -#include "source_io/module_chgpot/write_init.h" namespace elecstate { @@ -26,8 +25,6 @@ void ElecState::init_nelec_spin() } - - void ElecState::init_ks(Charge* chr_in, // pointer for class Charge const K_Vectors* klist_in, int nk_in,