From d14c9dc0002f30b864325f482aa911f285e60427 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 24 May 2026 15:00:54 +0800 Subject: [PATCH 01/19] 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/19] 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/19] 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/19] =?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/19] =?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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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