From f33a883cbb5b82287f0caffe112a1651133cc3bf Mon Sep 17 00:00:00 2001 From: DanielZhao0432 Date: Tue, 13 Oct 2026 13:48:44 +0800 Subject: [PATCH 1/8] fix(source_io): label TDOS histogram bins with left edge to align with PDOS grid cal_dos() recorded each energy bin [E, E+de) at its right edge, so the TDOS energy grid started at emin+de while the PDOS text grid starts at emin. Both files now share the same emin + n*de grid, allowing point-by-point comparison of TDOS and PDOS output. Only the energy column shifts by one step; dos/dos_int/dos_smear values are unchanged. Verified: ctest --test-dir build -V -R MODULE_DOS (3/3 passed). --- source/source_io/module_dos/cal_dos.cpp | 4 +++- source/source_io/module_dos/unittests/test_cal_dos.cpp | 4 ++-- .../source_io/module_dos/unittests/test_write_dos_pw.cpp | 8 ++++---- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/source/source_io/module_dos/cal_dos.cpp b/source/source_io/module_dos/cal_dos.cpp index bb6eba8b60..6d0a1553b1 100644 --- a/source/source_io/module_dos/cal_dos.cpp +++ b/source/source_io/module_dos/cal_dos.cpp @@ -210,7 +210,9 @@ bool ModuleIO::cal_dos(const int& is, // index for spin if (GlobalV::MY_RANK == 0) { dos.push_back(nstates); - ene.push_back(curr_energy); + // label the bin [e_old, curr_energy) with its left edge so that the + // TDOS grid (emin + n*de_ev) matches the PDOS grid point-by-point + ene.push_back(e_old); sum_elec.push_back(sum); } } diff --git a/source/source_io/module_dos/unittests/test_cal_dos.cpp b/source/source_io/module_dos/unittests/test_cal_dos.cpp index 62d8dd73bc..47581b79e3 100644 --- a/source/source_io/module_dos/unittests/test_cal_dos.cpp +++ b/source/source_io/module_dos/unittests/test_cal_dos.cpp @@ -63,8 +63,8 @@ TEST_F(DosTest,Dos) ifs.open(dosp.fa.c_str()); std::string str((std::istreambuf_iterator(ifs)),std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("4801 # number of points")); - EXPECT_THAT(str, testing::HasSubstr(" -4.995000 0.500000 0.500000")); - EXPECT_THAT(str, testing::HasSubstr(" 18.000000 0.000000 8.000000")); + EXPECT_THAT(str, testing::HasSubstr(" -5.000000 0.500000 0.500000")); + EXPECT_THAT(str, testing::HasSubstr(" 17.995000 0.000000 8.000000")); ifs.close(); remove("doss1_pw.txt"); #ifdef __MPI diff --git a/source/source_io/module_dos/unittests/test_write_dos_pw.cpp b/source/source_io/module_dos/unittests/test_write_dos_pw.cpp index 4ad1091604..6f99720567 100644 --- a/source/source_io/module_dos/unittests/test_write_dos_pw.cpp +++ b/source/source_io/module_dos/unittests/test_write_dos_pw.cpp @@ -103,8 +103,8 @@ TEST_F(DosPWTest,Dos1) ifs.open("dos.txt"); std::string str((std::istreambuf_iterator(ifs)),std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("4801 # number of points")); - EXPECT_THAT(str, testing::HasSubstr(" -4.995000 0.500000 0.500000")); - EXPECT_THAT(str, testing::HasSubstr(" 18.000000 0.000000 8.000000")); + EXPECT_THAT(str, testing::HasSubstr(" -5.000000 0.500000 0.500000")); + EXPECT_THAT(str, testing::HasSubstr(" 17.995000 0.000000 8.000000")); ifs.close(); remove("dos.txt"); #ifdef __MPI @@ -181,8 +181,8 @@ TEST_F(DosPWTest,Dos2) ifs.open("dos.txt"); std::string str1((std::istreambuf_iterator(ifs)),std::istreambuf_iterator()); EXPECT_THAT(str1, testing::HasSubstr("2526 # number of points")); - EXPECT_THAT(str1, testing::HasSubstr(" -4.997500 0.500000 0.500000")); - EXPECT_THAT(str1, testing::HasSubstr(" 7.567500 0.000000 8.000000")); + EXPECT_THAT(str1, testing::HasSubstr(" -5.002500 0.500000 0.500000")); + EXPECT_THAT(str1, testing::HasSubstr(" 7.562500 0.000000 8.000000")); ifs.close(); remove("dos.txt"); #ifdef __MPI From e48e7ba41762bd756aeb05499ba9246ec8274134 Mon Sep 17 00:00:00 2001 From: DanielZhao0432 Date: Tue, 13 Oct 2026 13:57:51 +0800 Subject: [PATCH 2/8] refactor(source_io): unify PW/LCAO DOS output naming to doss{spin}g{geom}.txt - write_dos_pw: dos.txt/dos1/dos2.txt -> doss{spin}g{geom}.txt, the same scheme as LCAO; the geometry index is now always tagged by passing the real istep in ctrl_output_pw - write_dos_lcao / write_pdos_text: drop the _nao/_pw basis suffix (pdoss{spin}g{geom}.txt); the now-unused basis parameter of write_pdos_text is removed - docs (parameters.yaml, input-main.md, dos.md, output-specification.md), examples/13_dos and tests/performance/catch_properties.sh updated accordingly Verified: ctest --test-dir build -V -R MODULE_DOS (3/3 passed). --- docs/advanced/elec_properties/dos.md | 6 +++-- docs/advanced/input_files/input-main.md | 2 +- .../output_files/output-specification.md | 6 ++--- docs/parameters.yaml | 2 +- examples/13_dos/01_pw_Al/run.sh | 2 +- examples/13_dos/02_lcao_Si2/run.sh | 2 +- examples/13_dos/README | 2 +- .../source_io/module_ctrl/ctrl_output_pw.cpp | 4 ++- .../source_io/module_dos/cal_pdos_gamma.cpp | 2 +- .../source_io/module_dos/cal_pdos_multik.cpp | 2 +- .../unittests/test_write_dos_pw.cpp | 12 ++++----- .../unittests/test_write_pdos_text.cpp | 25 +++++++++---------- .../source_io/module_dos/write_dos_lcao.cpp | 2 +- source/source_io/module_dos/write_dos_pw.cpp | 11 +++----- .../source_io/module_dos/write_pdos_text.cpp | 3 +-- source/source_io/module_dos/write_pdos_text.h | 3 +-- .../module_parameter/read_inp_out.cpp | 2 +- tests/performance/catch_properties.sh | 4 +-- 18 files changed, 45 insertions(+), 47 deletions(-) diff --git a/docs/advanced/elec_properties/dos.md b/docs/advanced/elec_properties/dos.md index d9f78bc4cc..acbbf60348 100644 --- a/docs/advanced/elec_properties/dos.md +++ b/docs/advanced/elec_properties/dos.md @@ -68,7 +68,9 @@ Gamma 8 8 8 0 0 0 ``` -Run the program, and you will see a file named doss1g1_nao.txt in the output directory. The columns are: energy(eV), dos (number of states in each energy bin), dos_int (integrated DOS), dos_smear(1/eV) (Gaussian-smeared DOS), dos_smear_int (integrated smeared DOS). Plot the file with graphing software, and you'll get the DOS. +Run the program, and you will see a file named doss1g1.txt in the output directory. The columns are: energy(eV), dos (number of states in each energy bin), dos_int (integrated DOS), dos_smear(1/eV) (Gaussian-smeared DOS), dos_smear_int (integrated smeared DOS). Plot the file with graphing software, and you'll get the DOS. + +The `dos` column is a raw per-bin histogram, while `dos_smear` is the physically meaningful density of states in 1/eV. The `abacus-plot` tool plots `dos_smear` by default; set `"smearing": false` in its config.json to plot the raw histogram instead, or apply your own smearing to the raw column with `abacus_plot.utils.gaussian_smearing`. For nspin=2, two files (doss1*/doss2*) are written, one per spin channel; the tool pairs them automatically whether you pass either file, a directory, or a list of both files. For nspin=4, a single spinor-summed file is written. ``` # energy(eV) dos dos_int dos_smear(1/eV) dos_smear_int @@ -80,7 +82,7 @@ Run the program, and you will see a file named doss1g1_nao.txt in the output dir ## PDOS -Along with the DOS files, we also produce the projected density of states (PDOS) in files named pdoss{spin}g{geom}_{basis}.txt (e.g., pdoss1g1_nao.txt). +Along with the DOS files, we also produce the projected density of states (PDOS) in files named pdoss{spin}g{geom}.txt (e.g., pdoss1g1.txt). The PDOS file uses a plain-text format. Each row corresponds to one energy point and one atom. Columns: energy(eV), atom (1-based), species, then pdos values ordered as s(1), p(3), d(5), f(7), etc. Zeta components are summed, and values below 1e-6 are zeroed out. diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 366c98bbd4..b7b800bcd7 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -2063,7 +2063,7 @@ - **Description**: Whether to output the density of states (DOS) and projected density of states (PDOS). For more information, refer to the dos.md. - 0: no output - 1: output the density of states (DOS) and projected density of states (PDOS, LCAO only) - - nspin=1 or 4: doss1g{geom}_{basis}.txt and pdoss1g{geom}_{basis}.txt, where geom is the geometry index when cell changes or ions move while basis is either pw or nao. + - nspin=1 or 4: doss1g{geom}.txt and pdoss1g{geom}.txt (PDOS is LCAO only), where geom is the geometry index when cell changes or ions move. PW and LCAO calculations use the same naming. - nspin=2: doss1/doss2 and pdoss1/pdoss2 files for two spin channels. Note: values 2 and 3 are no longer supported. Setting out_dos to 2 or 3 will raise an error. - **Default**: 0 diff --git a/docs/advanced/output_files/output-specification.md b/docs/advanced/output_files/output-specification.md index 834d18b687..71e5502686 100644 --- a/docs/advanced/output_files/output-specification.md +++ b/docs/advanced/output_files/output-specification.md @@ -86,7 +86,7 @@ All output file naming conventions can be found in the online documentation (wit | `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 | +| `doss1g1.txt` | DOS, spin 1, geometry step 1 | | `wf_pw.dat` | Wavefunction, plane wave basis | | `sr.csr` | Overlap matrix in real space (no spin index) | @@ -95,7 +95,7 @@ All output file naming conventions can be found in the online documentation (wit | File Name | Description | |-----------|-------------| | `running_scf.log` | SCF iteration log | -| `dos.txt` | Density of states | +| `doss1g1.txt` | Density of states | | `eig_occ.txt` | Eigenvalues and occupations | | `mulliken.txt` | Mulliken population analysis | | `band.txt` | Band structure | @@ -308,7 +308,7 @@ void OutputMyData::write(const std::string& filename, int precision) | File | Location | Strengths | |------|----------|-----------| | `chgs1.cube` | `tests/03_NAO_multik/scf_out_chg_tau/OUT.autotest/` | Short name, spin index convention | -| `dos.txt` | `tests/03_NAO_multik/scf_out_dos_spin2/OUT.autotest/` | Clear header with metadata | +| `doss1g1.txt` | `tests/03_NAO_multik/scf_out_dos_spin2/OUT.autotest/` | Clear header with metadata | | `eig_occ.txt` | Same as above | Clear block structure for spin/k-point | | `mulliken.txt` | `tests/03_NAO_multik/scf_out_mul/OUT.autotest/` | Clear atom separators | diff --git a/docs/parameters.yaml b/docs/parameters.yaml index a6d202867b..4ab89fd3cc 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -3449,7 +3449,7 @@ parameters: Whether to output the density of states (DOS) and projected density of states (PDOS). For more information, refer to the dos.md. * 0: no output * 1: output the density of states (DOS) and projected density of states (PDOS, LCAO only) - * nspin=1 or 4: doss1g{geom}_{basis}.txt and pdoss1g{geom}_{basis}.txt, where geom is the geometry index when cell changes or ions move while basis is either pw or nao. + * nspin=1 or 4: doss1g{geom}.txt and pdoss1g{geom}.txt (PDOS is LCAO only), where geom is the geometry index when cell changes or ions move. PW and LCAO calculations use the same naming. * nspin=2: doss1/doss2 and pdoss1/pdoss2 files for two spin channels. Note: values 2 and 3 are no longer supported. Setting out_dos to 2 or 3 will raise an error. default_value: "0" diff --git a/examples/13_dos/01_pw_Al/run.sh b/examples/13_dos/01_pw_Al/run.sh index 618ea42960..626f7008c3 100755 --- a/examples/13_dos/01_pw_Al/run.sh +++ b/examples/13_dos/01_pw_Al/run.sh @@ -18,7 +18,7 @@ if [[ ! -f scf.output ]] || [[ ! -f nscf.output ]] || [[ ! -f OUT.ABACUS/running_scf.log ]] || [[ ! -f OUT.ABACUS/running_nscf.log ]] || - [[ ! -f OUT.ABACUS/dos.txt ]] || + [[ ! -f OUT.ABACUS/doss1g1.txt ]] || [[ ! -f OUT.ABACUS/chg.cube ]] || [[ ! ( "$(tail -1 OUT.ABACUS/running_scf.log)" == " Total Time :"* ) ]] || [[ ! ( "$(tail -1 OUT.ABACUS/running_nscf.log)" == " Total Time :"* ) ]] diff --git a/examples/13_dos/02_lcao_Si2/run.sh b/examples/13_dos/02_lcao_Si2/run.sh index fa91c0046c..626f7008c3 100755 --- a/examples/13_dos/02_lcao_Si2/run.sh +++ b/examples/13_dos/02_lcao_Si2/run.sh @@ -18,7 +18,7 @@ if [[ ! -f scf.output ]] || [[ ! -f nscf.output ]] || [[ ! -f OUT.ABACUS/running_scf.log ]] || [[ ! -f OUT.ABACUS/running_nscf.log ]] || - [[ ! -f OUT.ABACUS/doss1g1_nao.txt ]] || + [[ ! -f OUT.ABACUS/doss1g1.txt ]] || [[ ! -f OUT.ABACUS/chg.cube ]] || [[ ! ( "$(tail -1 OUT.ABACUS/running_scf.log)" == " Total Time :"* ) ]] || [[ ! ( "$(tail -1 OUT.ABACUS/running_nscf.log)" == " Total Time :"* ) ]] diff --git a/examples/13_dos/README b/examples/13_dos/README index fb4cbe244f..fa353354d7 100644 --- a/examples/13_dos/README +++ b/examples/13_dos/README @@ -16,7 +16,7 @@ It requires two steps of calculation: - set 'out_dos' to '1' - set 'pw_diag_thr' to a suitable value when basis_type is 'pw' -Finish! Then you can see the DOS file in OUT.ABACUS (dos.txt for PW basis, doss1g{geom}_nao.txt for LCAO basis). The columns are: energy(eV), dos(1/eV), dos_int (integrated DOS), dos_smear(1/eV), dos_smear_int (integrated smeared DOS). Plot it! +Finish! Then you can see the DOS file in OUT.ABACUS (doss1g{geom}.txt for both PW and LCAO basis; for nspin=2 there are doss1*/doss2* files, one per spin channel). The columns are: energy(eV), dos(1/eV), dos_int (integrated DOS), dos_smear(1/eV), dos_smear_int (integrated smeared DOS). Plot it! Note: 1)If the "out_dos" is used, the DOS file will appear in the directory. diff --git a/source/source_io/module_ctrl/ctrl_output_pw.cpp b/source/source_io/module_ctrl/ctrl_output_pw.cpp index 6fa96d9e1e..9a107bd090 100644 --- a/source/source_io/module_ctrl/ctrl_output_pw.cpp +++ b/source/source_io/module_ctrl/ctrl_output_pw.cpp @@ -116,7 +116,9 @@ void ModuleIO::ctrl_scf_pw(const int istep, { bool out_dos_tmp = false; - int istep_in = -1; + // always tag the DOS file with the geometry index, so that PW and + // LCAO share the same doss{spin}g{geom}.txt naming + int istep_in = istep; // default value of out_freq_ion is 0 if (inp.out_freq_ion == 0) diff --git a/source/source_io/module_dos/cal_pdos_gamma.cpp b/source/source_io/module_dos/cal_pdos_gamma.cpp index 1bb5273b0e..3b9ee1b1dd 100644 --- a/source/source_io/module_dos/cal_pdos_gamma.cpp +++ b/source/source_io/module_dos/cal_pdos_gamma.cpp @@ -148,7 +148,7 @@ void ModuleIO::cal_pdos( if (GlobalV::MY_RANK == 0) { write_pdos_text(ucell, pdos.data(), nspin, nlocal, npoints, - emin, dos_edelta_ev, global_out_dir, "nao", istep); + emin, dos_edelta_ev, global_out_dir, istep); ModuleIO::write_orb_info(&ucell); } } diff --git a/source/source_io/module_dos/cal_pdos_multik.cpp b/source/source_io/module_dos/cal_pdos_multik.cpp index c20668f492..95ee8c7ba1 100644 --- a/source/source_io/module_dos/cal_pdos_multik.cpp +++ b/source/source_io/module_dos/cal_pdos_multik.cpp @@ -188,7 +188,7 @@ void ModuleIO::cal_pdos( if (GlobalV::MY_RANK == 0) { write_pdos_text(ucell, pdos.data(), nspin, nlocal, npoints, - emin, dos_edelta_ev, global_out_dir, "nao", istep); + emin, dos_edelta_ev, global_out_dir, istep); ModuleIO::write_orb_info(&ucell); } } diff --git a/source/source_io/module_dos/unittests/test_write_dos_pw.cpp b/source/source_io/module_dos/unittests/test_write_dos_pw.cpp index 6f99720567..12b533e080 100644 --- a/source/source_io/module_dos/unittests/test_write_dos_pw.cpp +++ b/source/source_io/module_dos/unittests/test_write_dos_pw.cpp @@ -76,7 +76,7 @@ TEST_F(DosPWTest,Dos1) dosp.wg, *kv, dosp.nbands, - -1, // istep_in + 0, // istep_in fermi_energy, dosp.de_ev, dos_scale, @@ -100,13 +100,13 @@ TEST_F(DosPWTest,Dos1) { #endif std::ifstream ifs; - ifs.open("dos.txt"); + ifs.open("doss1g1.txt"); std::string str((std::istreambuf_iterator(ifs)),std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("4801 # number of points")); EXPECT_THAT(str, testing::HasSubstr(" -5.000000 0.500000 0.500000")); EXPECT_THAT(str, testing::HasSubstr(" 17.995000 0.000000 8.000000")); ifs.close(); - remove("dos.txt"); + remove("doss1g1.txt"); #ifdef __MPI } #endif @@ -154,7 +154,7 @@ TEST_F(DosPWTest,Dos2) dosp.wg, *kv, dosp.nbands, - -1, // istep_in + 0, // istep_in fermi_energy, dosp.de_ev, dos_scale, @@ -178,13 +178,13 @@ TEST_F(DosPWTest,Dos2) { #endif std::ifstream ifs; - ifs.open("dos.txt"); + ifs.open("doss1g1.txt"); std::string str1((std::istreambuf_iterator(ifs)),std::istreambuf_iterator()); EXPECT_THAT(str1, testing::HasSubstr("2526 # number of points")); EXPECT_THAT(str1, testing::HasSubstr(" -5.002500 0.500000 0.500000")); EXPECT_THAT(str1, testing::HasSubstr(" 7.562500 0.000000 8.000000")); ifs.close(); - remove("dos.txt"); + remove("doss1g1.txt"); #ifdef __MPI } #endif diff --git a/source/source_io/module_dos/unittests/test_write_pdos_text.cpp b/source/source_io/module_dos/unittests/test_write_pdos_text.cpp index 2ab6f9009b..00681a3cb5 100644 --- a/source/source_io/module_dos/unittests/test_write_pdos_text.cpp +++ b/source/source_io/module_dos/unittests/test_write_pdos_text.cpp @@ -69,7 +69,7 @@ class WritePdosTextTest : public ::testing::Test } }; -// Spin-unpolarized: writes pdoss1_.txt with s and 3 p columns. +// Spin-unpolarized: writes pdoss1.txt with s and 3 p columns. TEST_F(WritePdosTextTest, BasicOutput) { GlobalV::ofs_running.open("test_write_pdos_text.log"); @@ -80,7 +80,6 @@ TEST_F(WritePdosTextTest, BasicOutput) const double emin = -1.0; const double dos_edelta_ev = 0.5; const std::string out_dir = "./"; - const std::string basis = "lcao"; const int istep = -1; // pdos matrix: nlocal x npoints; value = orbital * 100 + point @@ -94,9 +93,9 @@ TEST_F(WritePdosTextTest, BasicOutput) } ModuleIO::write_pdos_text(*ucell, &pdos, nspin, nlocal, npoints, - emin, dos_edelta_ev, out_dir, basis, istep); + emin, dos_edelta_ev, out_dir, istep); - const std::string str = read_file("./pdoss1_lcao.txt"); + const std::string str = read_file("./pdoss1.txt"); // header lines EXPECT_THAT(str, testing::HasSubstr("# npoints: 5")); @@ -109,7 +108,7 @@ TEST_F(WritePdosTextTest, BasicOutput) // last data row: energy=1.0, n=4 -> s=4, p=104,204,304 EXPECT_THAT(str, testing::HasSubstr(" 1.000000 1 Si 4.000000 104.000000 204.000000 304.000000")); - remove("./pdoss1_lcao.txt"); + remove("./pdoss1.txt"); GlobalV::ofs_running.close(); remove("test_write_pdos_text.log"); } @@ -128,12 +127,12 @@ TEST_F(WritePdosTextTest, IstepFileName) pdos.zero_out(); ModuleIO::write_pdos_text(*ucell, &pdos, nspin, nlocal, npoints, - emin, dos_edelta_ev, "./", "lcao", 2); + emin, dos_edelta_ev, "./", 2); - const std::string str = read_file("./pdoss1g3_lcao.txt"); + const std::string str = read_file("./pdoss1g3.txt"); EXPECT_THAT(str, testing::HasSubstr("# istep: 3")); - remove("./pdoss1g3_lcao.txt"); + remove("./pdoss1g3.txt"); GlobalV::ofs_running.close(); remove("test_write_pdos_text.log"); } @@ -164,17 +163,17 @@ TEST_F(WritePdosTextTest, SpinPolarized) } ModuleIO::write_pdos_text(*ucell, pdos_arr, nspin, nlocal, npoints, - emin, dos_edelta_ev, "./", "lcao", -1); + emin, dos_edelta_ev, "./", -1); - const std::string s1 = read_file("./pdoss1_lcao.txt"); - const std::string s2 = read_file("./pdoss2_lcao.txt"); + const std::string s1 = read_file("./pdoss1.txt"); + const std::string s2 = read_file("./pdoss2.txt"); // spin 1: s orbital w=0, n=0 -> 1000 EXPECT_THAT(s1, testing::HasSubstr("1000.000000")); // spin 2: s orbital w=0, n=0 -> 2000 EXPECT_THAT(s2, testing::HasSubstr("2000.000000")); - remove("./pdoss1_lcao.txt"); - remove("./pdoss2_lcao.txt"); + remove("./pdoss1.txt"); + remove("./pdoss2.txt"); GlobalV::ofs_running.close(); remove("test_write_pdos_text.log"); } diff --git a/source/source_io/module_dos/write_dos_lcao.cpp b/source/source_io/module_dos/write_dos_lcao.cpp index 77ecb3c903..926d94bd4e 100644 --- a/source/source_io/module_dos/write_dos_lcao.cpp +++ b/source/source_io/module_dos/write_dos_lcao.cpp @@ -67,7 +67,7 @@ void write_dos_lcao( ss << "g" << istep+1; } - ss << "_nao.txt"; + ss << ".txt"; ModuleIO::cal_dos(is, ss.str(), diff --git a/source/source_io/module_dos/write_dos_pw.cpp b/source/source_io/module_dos/write_dos_pw.cpp index 502b398474..1ce707f59d 100644 --- a/source/source_io/module_dos/write_dos_pw.cpp +++ b/source/source_io/module_dos/write_dos_pw.cpp @@ -50,16 +50,13 @@ void ModuleIO::write_dos_pw( for (int is = 0; is < nspin0; ++is) { // DOS_ispin contains not smoothed dos + // unified naming with the LCAO DOS/PDOS output: doss{spin}g{geom}.txt std::stringstream ss; - ss << global_out_dir << "dos"; + ss << global_out_dir << "doss" << is + 1; - if(nspin0==2) + if(istep_in>=0) { - ss << "s" << is + 1; - } - else - { - // do nothing; + ss << "g" << istep_in+1; } ss << ".txt"; diff --git a/source/source_io/module_dos/write_pdos_text.cpp b/source/source_io/module_dos/write_pdos_text.cpp index 0afc5ae57c..f0b3331dcf 100644 --- a/source/source_io/module_dos/write_pdos_text.cpp +++ b/source/source_io/module_dos/write_pdos_text.cpp @@ -17,7 +17,6 @@ void ModuleIO::write_pdos_text( const double& emin, const double& dos_edelta_ev, const std::string& out_dir, - const std::string& basis, const int istep) { ModuleBase::TITLE("ModuleIO", "write_pdos_text"); @@ -33,7 +32,7 @@ void ModuleIO::write_pdos_text( { ss << "g" << istep + 1; } - ss << "_" << basis << ".txt"; + ss << ".txt"; std::ofstream ofs(ss.str().c_str()); diff --git a/source/source_io/module_dos/write_pdos_text.h b/source/source_io/module_dos/write_pdos_text.h index e7b7f3a5ad..dcb9a18462 100644 --- a/source/source_io/module_dos/write_pdos_text.h +++ b/source/source_io/module_dos/write_pdos_text.h @@ -8,7 +8,7 @@ namespace ModuleIO { /// Write PDOS in plain text long-table format. -/// File name: pdoss{spin}g{geom}_{basis}.txt +/// File name: pdoss{spin}g{geom}.txt /// Columns: energy(eV) atom species l z m pdos(1/eV) /// For nspin=4, the two spinor components are summed. void write_pdos_text( @@ -20,7 +20,6 @@ void write_pdos_text( const double& emin, const double& dos_edelta_ev, const std::string& out_dir, - const std::string& basis, const int istep); } diff --git a/source/source_io/module_parameter/read_inp_out.cpp b/source/source_io/module_parameter/read_inp_out.cpp index 8ac0f3a23f..7a9645af08 100644 --- a/source/source_io/module_parameter/read_inp_out.cpp +++ b/source/source_io/module_parameter/read_inp_out.cpp @@ -327,7 +327,7 @@ Also controled by out_freq_ion and out_app_flag. item.description = R"(Whether to output the density of states (DOS) and projected density of states (PDOS). For more information, refer to the dos.md. * 0: no output * 1: output the density of states (DOS) and projected density of states (PDOS, LCAO only) - * nspin=1 or 4: doss1g{geom}_{basis}.txt and pdoss1g{geom}_{basis}.txt, where geom is the geometry index when cell changes or ions move while basis is either pw or nao. + * nspin=1 or 4: doss1g{geom}.txt and pdoss1g{geom}.txt (PDOS is LCAO only), where geom is the geometry index when cell changes or ions move. PW and LCAO calculations use the same naming. * nspin=2: doss1/doss2 and pdoss1/pdoss2 files for two spin channels. Note: values 2 and 3 are no longer supported. Setting out_dos to 2 or 3 will raise an error.)"; item.default_value = "0"; diff --git a/tests/performance/catch_properties.sh b/tests/performance/catch_properties.sh index d8b3794c41..9ee3b44a8e 100644 --- a/tests/performance/catch_properties.sh +++ b/tests/performance/catch_properties.sh @@ -89,10 +89,10 @@ fi #echo $total_charge if ! test -z "$has_dos" && [ $has_dos -eq 1 ]; then - total_dos=`cat OUT.autotest/doss1g1_nao.txt | awk 'END {print}' | awk '{print $3}'` + total_dos=`cat OUT.autotest/doss1g1.txt | awk 'END {print}' | awk '{print $3}'` echo "totaldosref $total_dos" >> $1 fi -# smearing_dos=`sum_file OUT.autotest/doss1g1_nao.txt` +# smearing_dos=`sum_file OUT.autotest/doss1g1.txt` # echo "totaldossmearing $smearing_dos" >> $1 From c0b7bbed7f1b5d884511b25bd6f2483cc21cfa6c Mon Sep 17 00:00:00 2001 From: DanielZhao0432 Date: Tue, 13 Oct 2026 13:58:07 +0800 Subject: [PATCH 3/8] fix(tools): plot-tools TDOS: smeared column by default, nspin=2 spin pairing, x-axis label - TDOS plots the Gaussian-smeared column (dos_smear, unit 1/eV) by default instead of the raw per-bin histogram; a new smearing switch (plot(..., smearing=False) or "smearing": false in config.json) restores the raw column, and utils.gaussian_smearing() lets users apply their own smearing width to it - TDOS accepts per-spin file pairs of the new text format: tdosfile may be a list of both spin files, a single spin file (paired in both directions, so doss2* works too), or a directory; energy grids are checked before merging. nspin=4 single spinor-summed files work as a single channel - the x-axis label is set centrally in DOSPlot._set_figure (bottom-most subplot only), so TDOS figures now show 'Energy(eV)' and the duplicated set_xlabel calls in PDOS branches are removed --- tools/02_postprocessing/plot-tools/README.md | 3 +- .../plot-tools/abacus_plot/dos.py | 153 +++++++++++++----- .../plot-tools/abacus_plot/main.py | 3 +- .../plot-tools/abacus_plot/utils.py | 20 +++ 4 files changed, 141 insertions(+), 38 deletions(-) diff --git a/tools/02_postprocessing/plot-tools/README.md b/tools/02_postprocessing/plot-tools/README.md index 768db54648..bd17549f22 100644 --- a/tools/02_postprocessing/plot-tools/README.md +++ b/tools/02_postprocessing/plot-tools/README.md @@ -118,11 +118,12 @@ First, prepare a file named 'config.json' in json format: If you only want to plot total DOS, you can modify `pdosfile` to `tdosfile` and do not set `species` and `pdosfig`. | Property | Type | Note | | :------------: | :-------------------------------------------------------------------------: | :-------------------------------------------------------------------------------------------------------------------------------------------------------------------: | -| _tdosfile_ | `str` | Total DOS data file output from ABACUS | +| _tdosfile_ | `str` or `List[str]` | Total DOS data file(s) output from ABACUS. For `nspin=2`, pass one spin-channel file (`doss1*` or `doss2*`; the counterpart is paired automatically), a directory containing them, or a list of both files. For `nspin=4`, ABACUS writes a single spinor-summed file | | _pdosfile_ | `str` | Directory or path prefix containing the text PDOS files `pdoss1_*` (and `pdoss2_*` for spin-polarized runs) output from ABACUS | | _efermi_ | `float` | Fermi level in eV | | _energy_range_ | `list` | Range of energy in eV | | _shift_ | `bool` | If set `'true'`, it will evaluate band gap. Default: `'false'` | +| _smearing_ | `bool` | If `'true'` (default), plot the Gaussian-smeared DOS column (`dos_smear`, unit 1/eV); if `'false'`, plot the raw per-bin histogram column. The raw column is also accessible as `TDOS(...).dos` for custom smearing, e.g. with `abacus_plot.utils.gaussian_smearing` | | _index_ | `List[int]` or `Dict[str, List[int]]` or `Dict[str, Dict[str, List[int]]]]` | Extract PDOS of each atom e.g. [index_0, index_1...] or {index*0:[l_0, l_1, ...], ...}, [index_0:{l_0:[m_0, m_1, ...], ...}, ...] | | \_atom_index* | `List[int]` or `Dict[str, List[int]]` or `Dict[str, Dict[str, List[int]]]]` | Extract PDOS of each atom with same atom*index e.g. [atom_index_0, atom_index_1...] or {atom_index_0:[l_0, l_1, ...], ...}, [atom_index_0:{l_0:[m_0, m_1, ...], ...}, ...] | | \_species* | `List[str]` or `Dict[str, List[int]]` or `Dict[str, Dict[str, List[int]]]]` | Extract PDOS of each atom with same species e.g. [elem_0, elem_1...] or {elem*0:[l_0, l_1, ...], ...}, [elem_0:{l_0:[m_0, m_1, ...], ...}, ...] | diff --git a/tools/02_postprocessing/plot-tools/abacus_plot/dos.py b/tools/02_postprocessing/plot-tools/abacus_plot/dos.py index 3387b60902..527681b0a8 100644 --- a/tools/02_postprocessing/plot-tools/abacus_plot/dos.py +++ b/tools/02_postprocessing/plot-tools/abacus_plot/dos.py @@ -120,6 +120,16 @@ def _set_figure(self, energy_range: Sequence = [], dos_range: Sequence = [], not # x-axis if energy_range: self.ax.set_xlim(energy_range[0], energy_range[1]) + # label the x-axis only on the bottom-most subplot when several + # axes share the x-axis; plain axes have no subplotspec and are + # always labeled + spec = self.ax.get_subplotspec() + if spec is None or spec.is_last_row(): + if "xlabel_params" in self.plot_params.keys(): + self.ax.set_xlabel("Energy(eV)", ** + self.plot_params["xlabel_params"]) + else: + self.ax.set_xlabel("Energy(eV)", size=25) # notes if notes: @@ -162,20 +172,66 @@ def _set_figure(self, energy_range: Sequence = [], dos_range: Sequence = [], not class TDOS(DOS): - """Parse total DOS data""" + """Parse total DOS data + + ``tdosfile`` may be a single TDOS file, a list of one or two per-spin + files, or a directory containing them. ABACUS writes one 5-column text + file per spin channel (``doss1*``/``doss2*``; legacy PW names + ``dos1.txt``/``dos2.txt``). A single spin-channel file is paired with + its counterpart when present, searching in both directions, so passing + either spin file works. For nspin=4 ABACUS writes a single + spinor-summed file, which is handled like nspin=1. + """ - def __init__(self, tdosfile: PathLike = None) -> None: + def __init__(self, tdosfile: Union[PathLike, Sequence[PathLike]] = None) -> None: self.tdosfile = tdosfile self._read() super().__init__(self.nspin) - def _read(self) -> tuple: - """Read total DOS data file + def _find_spin_files(self) -> List[str]: + """Locate the per-spin TDOS file(s) from a file, a list, or a directory.""" + if isinstance(self.tdosfile, (list, tuple)): + files = [str(f) for f in self.tdosfile] + if not 1 <= len(files) <= 2: + raise ValueError( + f"tdosfile accepts 1 or 2 files, got {len(files)}") + return files + + path = Path(self.tdosfile) + if path.is_dir(): + candidates = [] + for spin in (1, 2): + matches = sorted(glob.glob(str(path / f"doss{spin}*.txt"))) + if not matches: + matches = sorted(glob.glob(str(path / f"dos{spin}.txt"))) + if matches: + candidates.append(matches[0]) + if not candidates: + raise FileNotFoundError( + f"No TDOS file found in directory '{path}'") + return candidates - :params tdosfile: string of TDOS data file + if path.is_file(): + name = path.name + for one, two in (("doss1", "doss2"), ("dos1", "dos2")): + if one in name or two in name: + spin1 = path.parent / name.replace(two, one) + spin2 = path.parent / name.replace(one, two) + files = [str(p) for p in (spin1, spin2) if p.is_file()] + return files + # any other file (e.g. legacy TDOS.dat): parse as is + return [str(path)] + + def _read_one(self, fname: str) -> Dict[str, Any]: + """Read one TDOS file; the returned keys map 1:1 to the file columns. + + Returns a dict with ``energy`` of shape (npoints, 1) and dos/dos_int/ + dos_smear/dos_smear_int arrays whose second dimension counts the spin + channels stored in this file (None where the column is absent, e.g. + the legacy 2/3-column formats carry no smeared DOS). """ rows = [] - with open(self.tdosfile) as f: + with open(fname) as f: for line in f: if line.startswith("#"): continue @@ -185,36 +241,70 @@ def _read(self) -> tuple: continue rows.append([float(x) for x in line.split()]) if not rows: - raise ValueError(f"Empty TDOS file: {self.tdosfile}") + raise ValueError(f"Empty TDOS file: {fname}") # New format has two header lines ("ionic step" / "number of points") # with a single column; data lines have 5 columns per spin channel. # Old format (DOS1_smearing.dat) has only energy + dos columns. ncols = len(rows[-1]) data = np.asarray([r for r in rows if len(r) == ncols]) + out = {"energy": None, "dos": None, "dos_int": None, + "dos_smear": None, "dos_smear_int": None} if ncols == 2: - self.nspin = 1 - self.energy, self.dos = np.hsplit(data, 2) + out["energy"], out["dos"] = np.hsplit(data, 2) elif ncols == 3: - self.nspin = 2 - self.energy, dos_up, dos_dw = np.hsplit(data, 3) - self.dos = np.hstack((dos_up, dos_dw)) + energy, dos_up, dos_dw = np.hsplit(data, 3) + out["energy"] = energy + out["dos"] = np.hstack((dos_up, dos_dw)) elif ncols == 5: - self.nspin = 1 - self.energy, self.dos, self.dos_int, self.dos_smear, self.dos_smear_int = np.hsplit(data, 5) + out["energy"], out["dos"], out["dos_int"], out["dos_smear"], out["dos_smear_int"] = np.hsplit(data, 5) elif ncols == 9: - self.nspin = 2 - self.energy, dos_up, dos_dw, dos_int_up, dos_int_dw, smear_up, smear_dw, smear_int_up, smear_int_dw = np.hsplit(data, 9) - self.dos = np.hstack((dos_up, dos_dw)) - self.dos_int = np.hstack((dos_int_up, dos_int_dw)) - self.dos_smear = np.hstack((smear_up, smear_dw)) - self.dos_smear_int = np.hstack((smear_int_up, smear_int_dw)) + energy, dos_up, dos_dw, int_up, int_dw, sm_up, sm_dw, sint_up, sint_dw = np.hsplit(data, 9) + out["energy"] = energy + out["dos"] = np.hstack((dos_up, dos_dw)) + out["dos_int"] = np.hstack((int_up, int_dw)) + out["dos_smear"] = np.hstack((sm_up, sm_dw)) + out["dos_smear_int"] = np.hstack((sint_up, sint_dw)) else: raise ValueError(f"Unrecognized TDOS format: {ncols} columns") + return out - def _shift_energy(self, efermi: float = 0, shift: bool = False, prec: float = 0.01): + def _read(self) -> None: + spin_files = self._find_spin_files() + parsed = [self._read_one(f) for f in spin_files] + + keys = ("dos", "dos_int", "dos_smear", "dos_smear_int") + self.energy = parsed[0]["energy"] + if len(parsed) == 1: + for key in keys: + setattr(self, key, parsed[0][key]) + self.nspin = parsed[0]["dos"].shape[1] + else: + if parsed[0]["dos"].shape[1] != 1 or parsed[1]["dos"].shape[1] != 1: + raise ValueError( + "Each TDOS file already contains two spin channels; " + "do not pass two files") + if not np.allclose(parsed[0]["energy"], parsed[1]["energy"]): + raise ValueError( + "Spin-up and spin-down TDOS files have different energy " + f"grids: '{spin_files[0]}' vs '{spin_files[1]}'") + self.nspin = 2 + for key in keys: + up, dw = parsed[0][key], parsed[1][key] + setattr(self, key, None if up is None or dw is None + else np.hstack((up, dw))) + + def _selected_dos(self, smearing: bool) -> np.ndarray: + """DOS array to plot or analyze: the smeared column (unit 1/eV) when + available and requested, otherwise the raw histogram column.""" + if smearing and self.dos_smear is not None: + return self.dos_smear + return self.dos + + def _shift_energy(self, efermi: float = 0, shift: bool = False, prec: float = 0.01, smearing: bool = True): + dos = self._selected_dos(smearing) if shift: vb, cb = self.set_vcband( - energy_minus_efermi(self.energy, efermi), self.dos, prec) + energy_minus_efermi(self.energy, efermi), dos, prec) self.info(vb, cb) energy_f = np.concatenate((vb.band, cb.band)) else: @@ -222,12 +312,13 @@ def _shift_energy(self, efermi: float = 0, shift: bool = False, prec: float = 0. return energy_f - def plot(self, fig: Figure, ax: Union[axes.Axes, Sequence[axes.Axes]], efermi: float = 0, energy_range: Sequence[float] = [], dos_range: Sequence[float] = [], shift: bool = False, prec: float = 0.01, **kwargs): + def plot(self, fig: Figure, ax: Union[axes.Axes, Sequence[axes.Axes]], efermi: float = 0, energy_range: Sequence[float] = [], dos_range: Sequence[float] = [], shift: bool = False, prec: float = 0.01, smearing: bool = True, **kwargs): - energy_f = self._shift_energy(efermi, shift, prec) + energy_f = self._shift_energy(efermi, shift, prec, smearing) + dos = self._selected_dos(smearing) dosplot = DOSPlot(fig, ax, self.nspin, **kwargs) - dosplot.ax = self._plot(dosplot, energy_f, self.dos, "TDOS") + dosplot.ax = self._plot(dosplot, energy_f, dos, "TDOS") if "notes" in dosplot.plot_params.keys(): dosplot._set_figure(energy_range, dos_range, notes=dosplot.plot_params["notes"]) @@ -282,7 +373,7 @@ def _find_spin_files(self) -> List[str]: if path.is_file(): name = path.name if name.startswith("pdoss"): - # A concrete spin-channel file e.g. pdoss1g1_nao.txt: use it + # A concrete spin-channel file e.g. pdoss1g1.txt: use it # directly as the spin-1 channel and look for the matching # spin-2 file by replacing "s1" with "s2" in the filename. spin1 = str(path) @@ -575,11 +666,6 @@ def _parial_plot(self, if isinstance(species, (list, tuple)): dosplot = DOSPlot(fig, ax[0], self.nspin, **kwargs) - if "xlabel_params" in dosplot.plot_params.keys(): - dosplot.ax.set_xlabel("Energy(eV)", ** - dosplot.plot_params["xlabel_params"]) - else: - dosplot.ax.set_xlabel("Energy(eV)", size=25) self._plot(dosplot, energy_f, tdos, "TDOS") for elem in dos.keys(): self._plot(dosplot, energy_f, dos[elem], elem) @@ -614,11 +700,6 @@ def _parial_plot(self, else: dosplot._set_figure(energy_range, dos_range) dosplots.append(dosplot) - if "xlabel_params" in dosplot.plot_params.keys(): - dosplot.ax.set_xlabel("Energy(eV)", ** - dosplot.plot_params["xlabel_params"]) - else: - dosplot.ax.set_xlabel("Energy(eV)", size=25) return dosplots diff --git a/tools/02_postprocessing/plot-tools/abacus_plot/main.py b/tools/02_postprocessing/plot-tools/abacus_plot/main.py index 123f62b114..e565e9fa76 100644 --- a/tools/02_postprocessing/plot-tools/abacus_plot/main.py +++ b/tools/02_postprocessing/plot-tools/abacus_plot/main.py @@ -70,9 +70,10 @@ def show_cmdline(cls, args): figsize = text.pop("figsize", (12, 10)) fig, ax = plt.subplots(figsize=figsize) prec = text.pop("prec", 0.01) + smearing = text.pop("smearing", True) tdos = TDOS(tdosfile) tdos.plot(fig, ax, efermi=efermi, shift=shift, - energy_range=energy_range, dos_range=dos_range, prec=prec, **text) + energy_range=energy_range, dos_range=dos_range, prec=prec, smearing=smearing, **text) dpi = text.pop("dpi", 400) fig.savefig(text.pop("tdosfig", "tdos.png"), dpi=dpi) diff --git a/tools/02_postprocessing/plot-tools/abacus_plot/utils.py b/tools/02_postprocessing/plot-tools/abacus_plot/utils.py index b42610c1bb..8f2f8e1c16 100644 --- a/tools/02_postprocessing/plot-tools/abacus_plot/utils.py +++ b/tools/02_postprocessing/plot-tools/abacus_plot/utils.py @@ -237,6 +237,26 @@ def energy_minus_efermi(energy: Sequence, efermi: float) -> np.ndarray: return np.array(energy)-efermi +def gaussian_smearing(energy: Sequence, dos: Sequence, sigma: float) -> np.ndarray: + """Apply Gaussian smearing to a raw (per-bin histogram) DOS. + + Follows the same convention as the smearing ABACUS applies internally: + ``dos_smear[j] = sum_i dos[i] * G(E_j - E_i)`` with a normalized + Gaussian G of width parameter ``sigma`` (the ``dos_sigma`` input). + + :params energy: energy grid in eV + :params dos: raw DOS (number of states per energy bin), shape + (npoints,) or (npoints, nspin) + :params sigma: Gaussian width parameter in eV + """ + energy = np.asarray(energy, dtype=float).flatten() + dos = np.asarray(dos, dtype=float) + width = np.sqrt(2.0) * sigma + delta = energy[:, None] - energy[None, :] + gauss = np.exp(-(delta / width) ** 2) / (np.sqrt(np.pi) * width) + return gauss @ dos + + angular_momentum_label = ['s', 'p', 'd', 'f', 'g'] From 73fad1f5349dd6c7b9be629153d1bd12669829da Mon Sep 17 00:00:00 2001 From: DanielZhao0432 Date: Wed, 14 Oct 2026 22:08:51 +0800 Subject: [PATCH 4/8] fix(tools): address PR review: matplotlib compat guard and TDOS spin-pairing robustness - guard get_subplotspec()/is_last_row() behind getattr/hasattr so plotting still works on matplotlib < 3.4 (setup.py declares no minimum matplotlib version) - canonicalize a two-file tdosfile list so the spin-1 channel always comes first (a reversed list no longer swaps the up/down labels) - directory discovery now pairs spin files by identical name suffix (doss1g1.txt only pairs doss2g1.txt, never doss2g2.txt) and falls back to the legacy PW dos.txt - docs: the raw dos column is described as states-per-bin (not 1/eV); the smearing switch note now covers its effect on the shift=true gap estimate; PDOS file patterns corrected to pdoss1*/pdoss2* Verified with abacus_dos_test cases: reversed list equals forward list, mixed-geometry directory stays single-spin, legacy dos.txt directory loads, TDOS/PDOS plots unchanged. --- examples/13_dos/README | 2 +- tools/02_postprocessing/plot-tools/README.md | 4 +- .../plot-tools/abacus_plot/dos.py | 48 ++++++++++++------- 3 files changed, 35 insertions(+), 19 deletions(-) diff --git a/examples/13_dos/README b/examples/13_dos/README index fa353354d7..3ac7578126 100644 --- a/examples/13_dos/README +++ b/examples/13_dos/README @@ -16,7 +16,7 @@ It requires two steps of calculation: - set 'out_dos' to '1' - set 'pw_diag_thr' to a suitable value when basis_type is 'pw' -Finish! Then you can see the DOS file in OUT.ABACUS (doss1g{geom}.txt for both PW and LCAO basis; for nspin=2 there are doss1*/doss2* files, one per spin channel). The columns are: energy(eV), dos(1/eV), dos_int (integrated DOS), dos_smear(1/eV), dos_smear_int (integrated smeared DOS). Plot it! +Finish! Then you can see the DOS file in OUT.ABACUS (doss1g{geom}.txt for both PW and LCAO basis; for nspin=2 there are doss1*/doss2* files, one per spin channel). The columns are: energy(eV), dos (number of states in each energy bin), dos_int (integrated DOS), dos_smear(1/eV), dos_smear_int (integrated smeared DOS). Plot it! Note: 1)If the "out_dos" is used, the DOS file will appear in the directory. diff --git a/tools/02_postprocessing/plot-tools/README.md b/tools/02_postprocessing/plot-tools/README.md index bd17549f22..3ed71aa070 100644 --- a/tools/02_postprocessing/plot-tools/README.md +++ b/tools/02_postprocessing/plot-tools/README.md @@ -119,11 +119,11 @@ If you only want to plot total DOS, you can modify `pdosfile` to `tdosfile` and | Property | Type | Note | | :------------: | :-------------------------------------------------------------------------: | :-------------------------------------------------------------------------------------------------------------------------------------------------------------------: | | _tdosfile_ | `str` or `List[str]` | Total DOS data file(s) output from ABACUS. For `nspin=2`, pass one spin-channel file (`doss1*` or `doss2*`; the counterpart is paired automatically), a directory containing them, or a list of both files. For `nspin=4`, ABACUS writes a single spinor-summed file | -| _pdosfile_ | `str` | Directory or path prefix containing the text PDOS files `pdoss1_*` (and `pdoss2_*` for spin-polarized runs) output from ABACUS | +| _pdosfile_ | `str` | Directory or path prefix containing the text PDOS files `pdoss1*` (and `pdoss2*` for spin-polarized runs), e.g. `pdoss1g1.txt`, output from ABACUS | | _efermi_ | `float` | Fermi level in eV | | _energy_range_ | `list` | Range of energy in eV | | _shift_ | `bool` | If set `'true'`, it will evaluate band gap. Default: `'false'` | -| _smearing_ | `bool` | If `'true'` (default), plot the Gaussian-smeared DOS column (`dos_smear`, unit 1/eV); if `'false'`, plot the raw per-bin histogram column. The raw column is also accessible as `TDOS(...).dos` for custom smearing, e.g. with `abacus_plot.utils.gaussian_smearing` | +| _smearing_ | `bool` | If `'true'` (default), plot the Gaussian-smeared DOS column (`dos_smear`, unit 1/eV); if `'false'`, plot the raw per-bin histogram column. The raw column is also accessible as `TDOS(...).dos` for custom smearing, e.g. with `abacus_plot.utils.gaussian_smearing`. This switch also selects the data used by the `shift=true` band-gap estimate; with the smeared column, Gaussian tails can narrow small apparent gaps, so the reported gap is a threshold-based estimate that depends on `dos_sigma` and `prec` | | _index_ | `List[int]` or `Dict[str, List[int]]` or `Dict[str, Dict[str, List[int]]]]` | Extract PDOS of each atom e.g. [index_0, index_1...] or {index*0:[l_0, l_1, ...], ...}, [index_0:{l_0:[m_0, m_1, ...], ...}, ...] | | \_atom_index* | `List[int]` or `Dict[str, List[int]]` or `Dict[str, Dict[str, List[int]]]]` | Extract PDOS of each atom with same atom*index e.g. [atom_index_0, atom_index_1...] or {atom_index_0:[l_0, l_1, ...], ...}, [atom_index_0:{l_0:[m_0, m_1, ...], ...}, ...] | | \_species* | `List[str]` or `Dict[str, List[int]]` or `Dict[str, Dict[str, List[int]]]]` | Extract PDOS of each atom with same species e.g. [elem_0, elem_1...] or {elem*0:[l_0, l_1, ...], ...}, [elem_0:{l_0:[m_0, m_1, ...], ...}, ...] | diff --git a/tools/02_postprocessing/plot-tools/abacus_plot/dos.py b/tools/02_postprocessing/plot-tools/abacus_plot/dos.py index 527681b0a8..ea849bc88b 100644 --- a/tools/02_postprocessing/plot-tools/abacus_plot/dos.py +++ b/tools/02_postprocessing/plot-tools/abacus_plot/dos.py @@ -121,10 +121,11 @@ def _set_figure(self, energy_range: Sequence = [], dos_range: Sequence = [], not if energy_range: self.ax.set_xlim(energy_range[0], energy_range[1]) # label the x-axis only on the bottom-most subplot when several - # axes share the x-axis; plain axes have no subplotspec and are - # always labeled - spec = self.ax.get_subplotspec() - if spec is None or spec.is_last_row(): + # axes share the x-axis. get_subplotspec()/is_last_row() require + # matplotlib >= 3.4; fall back to always labeling on older versions + get_spec = getattr(self.ax, "get_subplotspec", None) + spec = get_spec() if get_spec is not None else None + if spec is None or not hasattr(spec, "is_last_row") or spec.is_last_row(): if "xlabel_params" in self.plot_params.keys(): self.ax.set_xlabel("Energy(eV)", ** self.plot_params["xlabel_params"]) @@ -195,20 +196,35 @@ def _find_spin_files(self) -> List[str]: if not 1 <= len(files) <= 2: raise ValueError( f"tdosfile accepts 1 or 2 files, got {len(files)}") + # canonicalize spin order: the spin-1 channel always comes first + if len(files) == 2: + first, second = Path(files[0]).name, Path(files[1]).name + if ("doss2" in first or "dos2" in first) and \ + ("doss1" in second or "dos1" in second): + files = files[::-1] return files path = Path(self.tdosfile) if path.is_dir(): - candidates = [] - for spin in (1, 2): - matches = sorted(glob.glob(str(path / f"doss{spin}*.txt"))) - if not matches: - matches = sorted(glob.glob(str(path / f"dos{spin}.txt"))) - if matches: - candidates.append(matches[0]) - if not candidates: + matches1 = sorted(glob.glob(str(path / "doss1*.txt"))) + if not matches1: + matches1 = sorted(glob.glob(str(path / "dos1.txt"))) + if not matches1 and (path / "dos.txt").is_file(): + matches1 = [str(path / "dos.txt")] # legacy PW nspin=1 name + if not matches1: raise FileNotFoundError( f"No TDOS file found in directory '{path}'") + spin1 = Path(matches1[0]) + candidates = [str(spin1)] + # pair with the spin-2 file carrying the same name suffix, e.g. + # doss1g1.txt -> doss2g1.txt, so different geometry steps are + # never mixed + for one, two in (("doss1", "doss2"), ("dos1", "dos2")): + if one in spin1.name: + spin2 = spin1.parent / spin1.name.replace(one, two) + if spin2.is_file(): + candidates.append(str(spin2)) + break return candidates if path.is_file(): @@ -346,8 +362,8 @@ class PDOS(DOS): continuous ``l = 0..nwl``. Zeta components are already summed by ABACUS. ``pdosfile`` may point to a single spin-channel file - (``pdoss1_*.txt``), to the spin-up file of a spin-polarized pair, or to a - path/prefix from which the ``pdoss1_*`` and (optionally) ``pdoss2_*`` + (``pdoss1*.txt``), to the spin-up file of a spin-polarized pair, or to a + path/prefix from which the ``pdoss1*`` and (optionally) ``pdoss2*`` files are discovered automatically. """ @@ -475,7 +491,7 @@ def _read(self): :params pdosfile: path of a text PDOS file, a spin-up file of a spin-polarized pair, or a path/prefix used to discover the - ``pdoss1_*``/``pdoss2_*`` files. + ``pdoss1*``/``pdoss2*`` files. """ spin_files = self._find_spin_files() self.nspin = len(spin_files) @@ -763,7 +779,7 @@ def plot(self, # energy_range=energy_range, dos_range=dos_range, notes={'s': '(a)'}) # fig.savefig("tdos.png") - # directory or prefix containing the text PDOS files pdoss1_*[/pdoss2_*] + # directory or prefix containing the text PDOS files pdoss1*[/pdoss2*] pdosfile = "./pdos" pdos = PDOS(pdosfile) #species = {"Ag": [2], "Cl": [1], "In": [0]} From 2f78590092415599f1479cc52ea66ce90f51778c Mon Sep 17 00:00:00 2001 From: DanielZhao0432 Date: Wed, 14 Oct 2026 22:49:27 +0800 Subject: [PATCH 5/8] perf(tools): chunk gaussian_smearing to bound memory usage Accumulate the smeared DOS in row chunks (default 1024) instead of allocating the dense npoints x npoints Gaussian matrix, which would take ~3.2 GB at 20000 energy points. Numerical results are unchanged (identical to machine precision), verified against the on-disk smeared column and with 20000-point stress input (~0.7 GB peak RSS). --- .../plot-tools/abacus_plot/utils.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/tools/02_postprocessing/plot-tools/abacus_plot/utils.py b/tools/02_postprocessing/plot-tools/abacus_plot/utils.py index 8f2f8e1c16..d3c526155a 100644 --- a/tools/02_postprocessing/plot-tools/abacus_plot/utils.py +++ b/tools/02_postprocessing/plot-tools/abacus_plot/utils.py @@ -237,24 +237,35 @@ def energy_minus_efermi(energy: Sequence, efermi: float) -> np.ndarray: return np.array(energy)-efermi -def gaussian_smearing(energy: Sequence, dos: Sequence, sigma: float) -> np.ndarray: +def gaussian_smearing(energy: Sequence, dos: Sequence, sigma: float, chunk: int = 1024) -> np.ndarray: """Apply Gaussian smearing to a raw (per-bin histogram) DOS. Follows the same convention as the smearing ABACUS applies internally: ``dos_smear[j] = sum_i dos[i] * G(E_j - E_i)`` with a normalized Gaussian G of width parameter ``sigma`` (the ``dos_sigma`` input). + The result is accumulated in row chunks, so memory usage stays + O(chunk * npoints) instead of O(npoints**2) for the dense + Gaussian matrix. + :params energy: energy grid in eV :params dos: raw DOS (number of states per energy bin), shape (npoints,) or (npoints, nspin) :params sigma: Gaussian width parameter in eV + :params chunk: number of output rows computed per block """ energy = np.asarray(energy, dtype=float).flatten() dos = np.asarray(dos, dtype=float) width = np.sqrt(2.0) * sigma - delta = energy[:, None] - energy[None, :] - gauss = np.exp(-(delta / width) ** 2) / (np.sqrt(np.pi) * width) - return gauss @ dos + norm = np.sqrt(np.pi) * width + n = energy.size + out = np.empty((n,) + dos.shape[1:], dtype=float) + for start in range(0, n, chunk): + stop = min(start + chunk, n) + delta = energy[start:stop, None] - energy[None, :] + gauss = np.exp(-(delta / width) ** 2) / norm + out[start:stop] = gauss @ dos + return out angular_momentum_label = ['s', 'p', 'd', 'f', 'g'] From 7ffda6ab4d2608e2cfa0e98509b45f9ae5d9f7f8 Mon Sep 17 00:00:00 2001 From: DanielZhao0432 Date: Thu, 15 Oct 2026 16:51:06 +0800 Subject: [PATCH 6/8] refactor(source_io): node-sample the TDOS smeared column, aligned with PDOS cal_dos() now uses a fixed node grid of floor((emax-emin)/de)+1 points, the same formula as cal_pdos, so TDOS and PDOS files always have the same number of rows and point-by-point identical energy grids (the previous while-loop could drop or add the last row depending on floating-point accumulation). Bin boundaries accumulate from a single variable so adjacent bins share the exact same floating-point edge and no eigenvalue is double-counted. The dos_smear column is now evaluated directly from the exact eigenvalue positions with the same normalized-Gaussian convention as cal_pdos, instead of convolving the histogram: the sum of all PDOS orbitals now equals the TDOS smeared column point by point (verified to ~1e-5, i.e. the 6-decimal output precision, on Si2/Fe/SOC cases; electron counts from dos_int are unchanged: 16/54/38). This also replaces the O(npoints^2) histogram convolution with O(npoints x nstates) accumulation and cuts the per-bin MPI reductions to two vector reductions. Verified: ctest --test-dir build -V -R MODULE_DOS (3/3 passed). --- source/source_io/module_dos/cal_dos.cpp | 112 ++++++++---------- .../module_dos/unittests/test_cal_dos.cpp | 5 +- .../unittests/test_write_dos_pw.cpp | 4 +- .../plot-tools/abacus_plot/utils.py | 13 +- 4 files changed, 63 insertions(+), 71 deletions(-) diff --git a/source/source_io/module_dos/cal_dos.cpp b/source/source_io/module_dos/cal_dos.cpp index 6d0a1553b1..7975c5e6de 100644 --- a/source/source_io/module_dos/cal_dos.cpp +++ b/source/source_io/module_dos/cal_dos.cpp @@ -5,6 +5,8 @@ #include "source_base/global_variable.h" #include "source_base/parallel_reduce.h" +#include + void ModuleIO::prepare_dos(std::ofstream& ofs_running, const elecstate::Efermi &energy_fermi, const ModuleBase::matrix& ekb, @@ -126,15 +128,6 @@ bool ModuleIO::cal_dos(const int& is, // index for spin ofs_dos << istep+1 << " # ionic step" << std::endl; } - std::vector dos; - std::vector ene; - std::vector sum_elec; - std::vector dos_smear; // dos_smearing - dos.clear(); - ene.clear(); - sum_elec.clear(); - dos_smear.clear(); - #ifdef __MPI MPI_Barrier(MPI_COMM_WORLD); #endif @@ -169,17 +162,29 @@ bool ModuleIO::cal_dos(const int& is, // index for spin << std::setw(15) << "dos_smear_int" << std::endl; } - std::vector e_mod(npoints, 0.0); - - double sum = 0.0; - double curr_energy = emin_ev; - double e_old = 0.0; - - while (curr_energy < emax_ev) + // Node-sampled grid shared with the PDOS text output: en = emin + n*de. + // Each node carries the histogram count of bin [en, en+de) and the + // Gaussian-smeared DOS evaluated directly at en from the exact + // eigenvalue positions, using the same normalized-Gaussian convention + // as cal_pdos (width parameter bcoeff), so that the sum of all PDOS + // orbitals equals this smeared column point by point. + std::vector dos(npoints, 0.0); + std::vector ene(npoints, 0.0); + std::vector dos_smear(npoints, 0.0); + std::vector sum_elec(npoints, 0.0); + + const double b = std::sqrt(2.0) * bcoeff; + const double gauss_norm = std::sqrt(ModuleBase::PI) * b; + + // en accumulates by de_ev so that adjacent bins [en, en+de) share the + // exact same floating-point boundary and no eigenvalue is counted twice + // or skipped; the printed labels still coincide with the PDOS nodes + // emin + n*de at the 1e-6 eV output precision. + double en = emin_ev; + for (int n = 0; n < npoints; ++n) { - double nstates = 0.0; - e_old = curr_energy; - curr_energy += de_ev; + ene[n] = en; + const double bin_hi = en + de_ev; // nks is the number of k-points in the 'pool' for (int ik = 0; ik < nks; ik++) @@ -190,68 +195,47 @@ bool ModuleIO::cal_dos(const int& is, // index for spin // band index for (int ib = 0; ib < nbands; ib++) { - // compare et and e_old(curr_energy) in ev unit. - if (ekb(ik, ib) * ModuleBase::Ry_to_eV >= e_old - && ekb(ik, ib) * ModuleBase::Ry_to_eV < curr_energy) + const double e = ekb(ik, ib) * ModuleBase::Ry_to_eV; + if (e >= en && e < bin_hi) { - nstates += wk[ik] * nkstot; + dos[n] += wk[ik] * nkstot; } + const double de = en - e; + dos_smear[n] += wk[ik] * std::exp(-de * de / b / b) / gauss_norm; } } } + en = bin_hi; + } #ifdef __MPI - const int npool = GlobalV::KPAR * bndpar; - Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, nstates); + const int npool = GlobalV::KPAR * bndpar; + Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, dos.data(), npoints); + Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, dos_smear.data(), npoints); #endif - nstates = nstates / static_cast(nkstot); - sum += nstates; - if (GlobalV::MY_RANK == 0) - { - dos.push_back(nstates); - // label the bin [e_old, curr_energy) with its left edge so that the - // TDOS grid (emin + n*de_ev) matches the PDOS grid point-by-point - ene.push_back(e_old); - sum_elec.push_back(sum); - } - } + double sum = 0.0; + double sum2 = 0.0; + const double dos_thr = 1.0e-12; - // Use Gaussian smearing to smooth the DOS if (GlobalV::MY_RANK == 0) { - dos_smear.resize(dos.size()); - - double b = sqrt(2.0) * bcoeff; - for (int i = 0; i < dos.size() ; i++) + for (int n = 0; n < npoints; ++n) { - double Gauss = 0.0; - - for (int j = 0; j < dos.size(); j++) - { - double denergy = ene[j] - ene[i]; - double de2 = denergy * denergy; - Gauss = exp(-de2 / b / b) / sqrt(ModuleBase::PI) / b; - dos_smear[j] += dos[i] * Gauss; - } - } + dos[n] /= static_cast(nkstot); + sum += dos[n]; + sum_elec[n] = sum; - // mohan add 2025-06-08 - const double dos_thr = 1.0e-12; - double sum2 = 0.0; - - for (int i = 0; i < dos.size(); i++) - { - if(dos_smear[i](ifs)),std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("4801 # number of points")); EXPECT_THAT(str, testing::HasSubstr(" -5.000000 0.500000 0.500000")); - EXPECT_THAT(str, testing::HasSubstr(" 17.995000 0.000000 8.000000")); + EXPECT_THAT(str, testing::HasSubstr(" 18.000000 0.000000 8.000000")); + // node-sampled smeared DOS: Gaussian evaluated at the exact eigenvalue + // position, e.g. node -5.000 eV gets 0.5*G(0) = 0.5/(sqrt(2*pi)*0.07) + EXPECT_THAT(str, testing::HasSubstr("2.849588")); ifs.close(); remove("doss1_pw.txt"); #ifdef __MPI diff --git a/source/source_io/module_dos/unittests/test_write_dos_pw.cpp b/source/source_io/module_dos/unittests/test_write_dos_pw.cpp index 12b533e080..37819b5c1e 100644 --- a/source/source_io/module_dos/unittests/test_write_dos_pw.cpp +++ b/source/source_io/module_dos/unittests/test_write_dos_pw.cpp @@ -104,7 +104,8 @@ TEST_F(DosPWTest,Dos1) std::string str((std::istreambuf_iterator(ifs)),std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("4801 # number of points")); EXPECT_THAT(str, testing::HasSubstr(" -5.000000 0.500000 0.500000")); - EXPECT_THAT(str, testing::HasSubstr(" 17.995000 0.000000 8.000000")); + EXPECT_THAT(str, testing::HasSubstr(" 18.000000 0.000000 8.000000")); + EXPECT_THAT(str, testing::HasSubstr("2.849588")); ifs.close(); remove("doss1g1.txt"); #ifdef __MPI @@ -183,6 +184,7 @@ TEST_F(DosPWTest,Dos2) EXPECT_THAT(str1, testing::HasSubstr("2526 # number of points")); EXPECT_THAT(str1, testing::HasSubstr(" -5.002500 0.500000 0.500000")); EXPECT_THAT(str1, testing::HasSubstr(" 7.562500 0.000000 8.000000")); + EXPECT_THAT(str1, testing::HasSubstr("2.847771")); ifs.close(); remove("doss1g1.txt"); #ifdef __MPI diff --git a/tools/02_postprocessing/plot-tools/abacus_plot/utils.py b/tools/02_postprocessing/plot-tools/abacus_plot/utils.py index d3c526155a..1ee51f9459 100644 --- a/tools/02_postprocessing/plot-tools/abacus_plot/utils.py +++ b/tools/02_postprocessing/plot-tools/abacus_plot/utils.py @@ -238,11 +238,14 @@ def energy_minus_efermi(energy: Sequence, efermi: float) -> np.ndarray: def gaussian_smearing(energy: Sequence, dos: Sequence, sigma: float, chunk: int = 1024) -> np.ndarray: - """Apply Gaussian smearing to a raw (per-bin histogram) DOS. - - Follows the same convention as the smearing ABACUS applies internally: - ``dos_smear[j] = sum_i dos[i] * G(E_j - E_i)`` with a normalized - Gaussian G of width parameter ``sigma`` (the ``dos_sigma`` input). + """Apply Gaussian smearing to a raw (per-bin histogram) DOS column. + + This mirrors the histogram-based convention: ``dos_smear[j] = sum_i + dos[i] * G(E_j - E_i)`` with a normalized Gaussian G of width parameter + ``sigma`` (the ``dos_sigma`` input). Note that the ``dos_smear`` column + written by ABACUS is instead computed by smearing the eigenvalues at + their exact positions; the two differ by the O(bin-width) binning + error of the histogram, which is expected. The result is accumulated in row chunks, so memory usage stays O(chunk * npoints) instead of O(npoints**2) for the dense From a13ddb4756a5dede67c3b02fda2a220d2a84d808 Mon Sep 17 00:00:00 2001 From: ABACUS-PR-BOT <342302420+ABACUS-PR-BOT@users.noreply.github.com> Date: Fri, 16 Oct 2026 07:37:38 +0800 Subject: [PATCH 7/8] Fix DOS plotting review follow-ups and add focused regressions --- .../module_dos/unittests/test_cal_dos.cpp | 35 +++++ .../unittests/test_write_dos_pw.cpp | 36 +++++ tools/02_postprocessing/plot-tools/README.md | 11 ++ .../plot-tools/abacus_plot/dos.py | 7 +- .../plot-tools/abacus_plot/main.py | 5 +- .../plot-tools/abacus_plot/utils.py | 4 + .../plot-tools/examples/Si/README | 2 +- .../plot-tools/tests/test_dos.py | 127 ++++++++++++++++++ .../plot-tools/tests/test_main.py | 57 ++++++++ .../plot-tools/tests/test_utils.py | 40 ++++++ 10 files changed, 321 insertions(+), 3 deletions(-) create mode 100644 tools/02_postprocessing/plot-tools/tests/test_dos.py create mode 100644 tools/02_postprocessing/plot-tools/tests/test_main.py create mode 100644 tools/02_postprocessing/plot-tools/tests/test_utils.py diff --git a/source/source_io/module_dos/unittests/test_cal_dos.cpp b/source/source_io/module_dos/unittests/test_cal_dos.cpp index eec4f7dcd7..abcf942bad 100644 --- a/source/source_io/module_dos/unittests/test_cal_dos.cpp +++ b/source/source_io/module_dos/unittests/test_cal_dos.cpp @@ -3,6 +3,7 @@ #include "../cal_dos.h" #include "source_base/global_variable.h" #include +#include #ifdef __MPI #include "mpi.h" #endif @@ -68,6 +69,40 @@ TEST_F(DosTest,Dos) // node-sampled smeared DOS: Gaussian evaluated at the exact eigenvalue // position, e.g. node -5.000 eV gets 0.5*G(0) = 0.5/(sqrt(2*pi)*0.07) EXPECT_THAT(str, testing::HasSubstr("2.849588")); + std::istringstream rows(str); + std::string line; + std::getline(rows, line); // ionic step + std::getline(rows, line); // declared grid size + int declared_points = 0; + std::istringstream(line) >> declared_points; + EXPECT_EQ(declared_points, 4801); + std::getline(rows, line); // column names + int row_count = 0; + while (std::getline(rows, line)) + { + std::istringstream fields(line); + double energy = 0.0; + double raw = 0.0; + double integral = 0.0; + double smeared = 0.0; + double smeared_integral = 0.0; + ASSERT_TRUE(fields >> energy >> raw >> integral >> smeared >> smeared_integral); + EXPECT_NEAR(energy, -6.0 + row_count * 0.005, 5.0e-7); + if (row_count == 200) + { + EXPECT_DOUBLE_EQ(energy, -5.0); + EXPECT_DOUBLE_EQ(raw, 0.5); + EXPECT_DOUBLE_EQ(integral, 0.5); + EXPECT_NEAR(smeared, 2.849588, 5.0e-7); + } + if (row_count == 4800) + { + EXPECT_DOUBLE_EQ(energy, 18.0); + EXPECT_DOUBLE_EQ(integral, 8.0); + } + ++row_count; + } + EXPECT_EQ(row_count, declared_points); ifs.close(); remove("doss1_pw.txt"); #ifdef __MPI diff --git a/source/source_io/module_dos/unittests/test_write_dos_pw.cpp b/source/source_io/module_dos/unittests/test_write_dos_pw.cpp index 37819b5c1e..384af45ac6 100644 --- a/source/source_io/module_dos/unittests/test_write_dos_pw.cpp +++ b/source/source_io/module_dos/unittests/test_write_dos_pw.cpp @@ -1,6 +1,7 @@ #include "gtest/gtest.h" #include "gmock/gmock.h" #include "../write_dos_pw.h" +#include #ifdef __MPI #include "mpi.h" #endif @@ -185,6 +186,41 @@ TEST_F(DosPWTest,Dos2) EXPECT_THAT(str1, testing::HasSubstr(" -5.002500 0.500000 0.500000")); EXPECT_THAT(str1, testing::HasSubstr(" 7.562500 0.000000 8.000000")); EXPECT_THAT(str1, testing::HasSubstr("2.847771")); + std::istringstream rows(str1); + std::string line; + std::getline(rows, line); // ionic step + std::getline(rows, line); // declared grid size + int declared_points = 0; + std::istringstream(line) >> declared_points; + EXPECT_EQ(declared_points, 2526); + std::getline(rows, line); // column names + int row_count = 0; + while (std::getline(rows, line)) + { + std::istringstream fields(line); + double energy = 0.0; + double raw = 0.0; + double integral = 0.0; + double smeared = 0.0; + double smeared_integral = 0.0; + ASSERT_TRUE(fields >> energy >> raw >> integral >> smeared >> smeared_integral); + EXPECT_NEAR(energy, -5.0625 + row_count * 0.005, 5.0e-7); + if (row_count == 12) + { + // The -5 eV state is 0.0025 eV above this automatic-grid node. + EXPECT_DOUBLE_EQ(energy, -5.0025); + EXPECT_DOUBLE_EQ(raw, 0.5); + EXPECT_DOUBLE_EQ(integral, 0.5); + EXPECT_NEAR(smeared, 2.847771, 5.0e-7); + } + if (row_count == 2525) + { + EXPECT_DOUBLE_EQ(energy, 7.5625); + EXPECT_DOUBLE_EQ(integral, 8.0); + } + ++row_count; + } + EXPECT_EQ(row_count, declared_points); ifs.close(); remove("doss1g1.txt"); #ifdef __MPI diff --git a/tools/02_postprocessing/plot-tools/README.md b/tools/02_postprocessing/plot-tools/README.md index 3ed71aa070..f1bb9962cd 100644 --- a/tools/02_postprocessing/plot-tools/README.md +++ b/tools/02_postprocessing/plot-tools/README.md @@ -94,6 +94,17 @@ Then, the following command will output parsed partial DOS to directory `PBAND*_ abacus-plot -d -o ``` +## Testing + +From this directory, run the focused DOS reader, plotting and smearing regressions: + +```shell +MPLBACKEND=Agg PYTHONPATH=. python3 -m unittest discover -s tests -v +``` + +These Python tests use the package dependencies from `setup.py` and are discovered +by `unittest`; they do not add a dependency to the C++ CMake build. + ### DOS First, prepare a file named 'config.json' in json format: diff --git a/tools/02_postprocessing/plot-tools/abacus_plot/dos.py b/tools/02_postprocessing/plot-tools/abacus_plot/dos.py index ea849bc88b..6c310eeb1c 100644 --- a/tools/02_postprocessing/plot-tools/abacus_plot/dos.py +++ b/tools/02_postprocessing/plot-tools/abacus_plot/dos.py @@ -299,7 +299,12 @@ def _read(self) -> None: raise ValueError( "Each TDOS file already contains two spin channels; " "do not pass two files") - if not np.allclose(parsed[0]["energy"], parsed[1]["energy"]): + energy_up = parsed[0]["energy"] + energy_dw = parsed[1]["energy"] + # ABACUS prints energies to 1e-6 eV. Keep the comparison below + # that precision, without a tolerance that grows with |energy|. + if energy_up.shape != energy_dw.shape or not np.allclose( + energy_up, energy_dw, rtol=0, atol=1e-8): raise ValueError( "Spin-up and spin-down TDOS files have different energy " f"grids: '{spin_files[0]}' vs '{spin_files[1]}'") diff --git a/tools/02_postprocessing/plot-tools/abacus_plot/main.py b/tools/02_postprocessing/plot-tools/abacus_plot/main.py index e565e9fa76..859a49094e 100644 --- a/tools/02_postprocessing/plot-tools/abacus_plot/main.py +++ b/tools/02_postprocessing/plot-tools/abacus_plot/main.py @@ -95,8 +95,11 @@ def show_cmdline(cls, args): fig, ax = plt.subplots( len(atom_index), 1, sharex=True, figsize=figsize) if species: + # Lists overlay components on one axes; dictionaries use + # one row per component for angular-momentum-resolved plots. fig, ax = plt.subplots( - len(species), 1, sharex=True, figsize=figsize) + len(species) if isinstance(species, dict) else 1, 1, + sharex=True, figsize=figsize) prec = text.pop("prec", 0.01) pdos = PDOS(pdosfile) pdos.plot(fig=fig, ax=ax, index=index, atom_index=atom_index, species=species, efermi=efermi, shift=shift, diff --git a/tools/02_postprocessing/plot-tools/abacus_plot/utils.py b/tools/02_postprocessing/plot-tools/abacus_plot/utils.py index 1ee51f9459..0559e9d6ab 100644 --- a/tools/02_postprocessing/plot-tools/abacus_plot/utils.py +++ b/tools/02_postprocessing/plot-tools/abacus_plot/utils.py @@ -257,6 +257,10 @@ def gaussian_smearing(energy: Sequence, dos: Sequence, sigma: float, chunk: int :params sigma: Gaussian width parameter in eV :params chunk: number of output rows computed per block """ + if np.ndim(sigma) != 0 or not np.isfinite(sigma) or sigma <= 0: + raise ValueError("sigma must be finite and positive") + if isinstance(chunk, (bool, np.bool_)) or not isinstance(chunk, (int, np.integer)) or chunk <= 0: + raise ValueError("chunk must be a positive integer") energy = np.asarray(energy, dtype=float).flatten() dos = np.asarray(dos, dtype=float) width = np.sqrt(2.0) * sigma diff --git a/tools/02_postprocessing/plot-tools/examples/Si/README b/tools/02_postprocessing/plot-tools/examples/Si/README index 24b1038ee0..494c3665bb 100644 --- a/tools/02_postprocessing/plot-tools/examples/Si/README +++ b/tools/02_postprocessing/plot-tools/examples/Si/README @@ -11,7 +11,7 @@ NiteL (1) The input files of this example is chosen from ABACUS/examples/02(c-d). Extra parameter "out_proj_band" for PBANDS calculation is needed in 02c. -(2) The PDOS is now read from the text files written by ABACUS (`pdoss1_*.txt`, and `pdoss2_*.txt` for +(2) The PDOS is now read from the text files written by ABACUS (`pdoss1*.txt`, and `pdoss2*.txt` for spin-polarized calculations) instead of the old XML file. Set `pdosfile` in config.json to the directory or prefix containing those files (e.g. the `OUT.*` directory, or the path prefix ending in `pdos`). diff --git a/tools/02_postprocessing/plot-tools/tests/test_dos.py b/tools/02_postprocessing/plot-tools/tests/test_dos.py new file mode 100644 index 0000000000..3655ad5b0c --- /dev/null +++ b/tools/02_postprocessing/plot-tools/tests/test_dos.py @@ -0,0 +1,127 @@ +"""Regressions for TDOS formats, spin pairing and plotted data.""" + +import tempfile +import unittest +from pathlib import Path +from types import SimpleNamespace +from unittest import mock + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +from abacus_plot.dos import DOSPlot, TDOS + + +class TDOSTest(unittest.TestCase): + def setUp(self): + self.temp = tempfile.TemporaryDirectory() + self.addCleanup(self.temp.cleanup) + self.addCleanup(plt.close, "all") + self.root = Path(self.temp.name) + self.energy = np.array([-1., 0., 1.]) + + def write_spin(self, name, energy, scale): + path = self.root / name + raw = scale * np.arange(1., len(energy) + 1.) + data = np.column_stack((energy, raw, np.cumsum(raw), + raw + 10., np.cumsum(raw + 10.))) + np.savetxt(path, data, header="1 # ionic step\n{} # number of points".format(len(energy)), + comments="", fmt="%.12f") + return path + + def test_formats_and_default_raw_selection(self): + up = np.array([1., 2., 3.]) + down = np.array([4., 5., 6.]) + smear_up = up + 10. + smear_down = down + 20. + formats = [ + (np.column_stack((self.energy, up)), up[:, None], None), + (np.column_stack((self.energy, up, down)), np.column_stack((up, down)), None), + (np.column_stack((self.energy, up, up.cumsum(), smear_up, smear_up.cumsum())), + up[:, None], smear_up[:, None]), + (np.column_stack((self.energy, up, down, up.cumsum(), down.cumsum(), + smear_up, smear_down, smear_up.cumsum(), smear_down.cumsum())), + np.column_stack((up, down)), np.column_stack((smear_up, smear_down))), + ] + for data, raw, smeared in formats: + with self.subTest(columns=data.shape[1]): + path = self.root / "format.txt" + np.savetxt(path, data) + tdos = TDOS(path) + np.testing.assert_array_equal(tdos.energy[:, 0], self.energy) + np.testing.assert_array_equal(tdos.dos, raw) + self.assertEqual(tdos.nspin, raw.shape[1]) + if smeared is None: + self.assertIsNone(tdos.dos_smear) + else: + np.testing.assert_array_equal(tdos.dos_smear, smeared) + for smearing in (True, False): + fig, ax = plt.subplots() + tdos.plot(fig, ax, smearing=smearing) + expected = smeared if smearing and smeared is not None else raw + np.testing.assert_array_equal(ax.lines[0].get_ydata(), expected[:, 0]) + if tdos.nspin == 2: + np.testing.assert_array_equal(ax.lines[1].get_ydata(), -expected[:, 1]) + self.assertEqual(ax.get_xlabel(), "Energy(eV)") + plt.close(fig) + + def test_shift_uses_the_selected_column(self): + path = self.write_spin("single.txt", self.energy, 1.) + tdos = TDOS(path) + bands = (SimpleNamespace(band=tdos.energy[:1], value=0.), + SimpleNamespace(band=tdos.energy[1:], value=1.)) + for smearing in (True, False): + with mock.patch.object(tdos, "set_vcband", return_value=bands) as split: + with mock.patch.object(tdos, "info"): + fig, ax = plt.subplots() + tdos.plot(fig, ax, shift=True, smearing=smearing) + expected = tdos.dos_smear if smearing else tdos.dos + np.testing.assert_array_equal(split.call_args[0][1], expected) + np.testing.assert_array_equal(ax.lines[0].get_ydata(), expected[:, 0]) + + def test_recognized_pairs_preserve_spin_order(self): + up = self.write_spin("doss1g1.txt", self.energy, 1.) + down = self.write_spin("doss2g1.txt", self.energy, 7.) + for source in (up, down, self.root, [up, down], [down, up]): + with self.subTest(source=str(source)): + tdos = TDOS(source) + self.assertEqual(tdos.nspin, 2) + np.testing.assert_array_equal(tdos.dos[:, 0], [1., 2., 3.]) + np.testing.assert_array_equal(tdos.dos[:, 1], [7., 14., 21.]) + + def test_directory_does_not_mix_geometries(self): + self.write_spin("doss1g1.txt", self.energy, 1.) + self.write_spin("doss2g2.txt", self.energy, 7.) + tdos = TDOS(self.root) + self.assertEqual(tdos.nspin, 1) + np.testing.assert_array_equal(tdos.dos[:, 0], [1., 2., 3.]) + + def test_spin_grids_require_equal_shape_and_absolute_agreement(self): + energy = np.array([100., 100.0001, 100.0002]) + up = self.write_spin("doss1g1.txt", energy, 1.) + for other in (energy + 0.0001, energy[:2], energy[:1]): + with self.subTest(other=other): + down = self.write_spin("doss2g1.txt", other, 7.) + with self.assertRaisesRegex(ValueError, "different energy grids"): + TDOS([up, down]) + for other in (energy, energy + 5e-9): + down = self.write_spin("doss2g1.txt", other, 7.) + self.assertEqual(TDOS([up, down]).nspin, 2) + + def test_plain_axes_and_missing_subplot_apis(self): + fig = plt.figure() + plain = fig.add_axes([0.1, 0.1, 0.8, 0.8]) + DOSPlot(fig, plain, xlabel_params={"size": 13})._set_figure() + self.assertEqual(plain.get_xlabel(), "Energy(eV)") + self.assertEqual(plain.xaxis.label.get_fontsize(), 13) + for get_spec in (None, lambda: object()): + fig, ax = plt.subplots() + with mock.patch.object(ax, "get_subplotspec", get_spec, create=True): + DOSPlot(fig, ax)._set_figure() + self.assertEqual(ax.get_xlabel(), "Energy(eV)") + + +if __name__ == "__main__": + unittest.main() diff --git a/tools/02_postprocessing/plot-tools/tests/test_main.py b/tools/02_postprocessing/plot-tools/tests/test_main.py new file mode 100644 index 0000000000..164d46a2f8 --- /dev/null +++ b/tools/02_postprocessing/plot-tools/tests/test_main.py @@ -0,0 +1,57 @@ +"""Exercise the actual CLI selector-to-PDOS plotting path.""" + +import json +import tempfile +import unittest +from pathlib import Path +from unittest import mock + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +from abacus_plot.main import main + + +class PDOSCommandTest(unittest.TestCase): + def test_list_and_dictionary_layouts(self): + with tempfile.TemporaryDirectory() as temp: + root = Path(temp) + pdos = root / "pdoss1g1.txt" + pdos.write_text("# npoints: 3\n" + "-1 1 Si 1\n-1 2 O 2\n" + "0 1 Si 2\n0 2 O 3\n" + "1 1 Si 3\n1 2 O 4\n") + selectors = {"species": ["Si", "O"], "atom_index": [1, 2], "index": [1, 2]} + for key, values in selectors.items(): + for as_dict in ((False, True) if key == "species" else (True,)): + with self.subTest(key=key, dictionary=as_dict): + plt.close("all") + image = root / "{}-{}.png".format(key, as_dict) + config = {"pdosfile": str(pdos), "pdosfig": str(image), + "figsize": [4, 3], "dpi": 80, + "xlabel_params": {"size": 13}, + key: {str(value): [0] for value in values} if as_dict else values} + config_path = root / "config.json" + config_path.write_text(json.dumps(config)) + with mock.patch("sys.argv", ["abacus-plot", "-d", "-p", "-f", str(config_path)]): + main() + fig = plt.gcf() + self.assertEqual(len(fig.axes), 2 if as_dict else 1) + self.assertEqual(fig.axes[-1].get_xlabel(), "Energy(eV)") + self.assertEqual(fig.axes[-1].xaxis.label.get_fontsize(), 13) + if as_dict: + self.assertEqual(fig.axes[0].get_xlabel(), "") + self.assertEqual(fig.axes[0].lines[0].get_ydata().tolist(), [1., 2., 3.]) + self.assertEqual(fig.axes[1].lines[0].get_ydata().tolist(), [2., 3., 4.]) + else: + self.assertEqual(len(fig.axes[0].get_legend_handles_labels()[1]), 3) + self.assertEqual(fig.axes[0].lines[0].get_ydata().tolist(), [3., 5., 7.]) + self.assertEqual(fig.axes[0].lines[1].get_ydata().tolist(), [1., 2., 3.]) + self.assertEqual(fig.axes[0].lines[2].get_ydata().tolist(), [2., 3., 4.]) + self.assertGreater(image.stat().st_size, 0) + plt.close("all") + + +if __name__ == "__main__": + unittest.main() diff --git a/tools/02_postprocessing/plot-tools/tests/test_utils.py b/tools/02_postprocessing/plot-tools/tests/test_utils.py new file mode 100644 index 0000000000..a2bf819ae0 --- /dev/null +++ b/tools/02_postprocessing/plot-tools/tests/test_utils.py @@ -0,0 +1,40 @@ +"""Gaussian helper validation and analytic histogram-smearing checks.""" + +import math +import unittest + +import numpy as np + +from abacus_plot.utils import gaussian_smearing + + +class GaussianSmearingTest(unittest.TestCase): + def test_invalid_width_and_chunk_are_rejected(self): + for sigma in (0., -1., float("nan"), float("inf"), -float("inf"), [1., 2.]): + with self.subTest(sigma=sigma): + with self.assertRaisesRegex(ValueError, "sigma"): + gaussian_smearing([-1., 0., 1.], [0., 1., 0.], sigma) + for chunk in (0, -1, 1.5, True, np.bool_(True)): + with self.subTest(chunk=chunk): + with self.assertRaisesRegex(ValueError, "chunk"): + gaussian_smearing([-1., 0., 1.], [0., 1., 0.], 1., chunk=chunk) + + def test_analytic_two_spin_result_is_chunk_independent(self): + energy = [-1., 0., 1.] + counts = np.array([[1., 0.], [0., 2.], [1., 0.]]) + norm = math.sqrt(2. * math.pi) + expected = np.array([[(1. + math.exp(-2.)) / norm, 2. * math.exp(-0.5) / norm], + [2. * math.exp(-0.5) / norm, 2. / norm], + [(1. + math.exp(-2.)) / norm, 2. * math.exp(-0.5) / norm]]) + for chunk in (1, 2, 1024, np.int64(2)): + with self.subTest(chunk=chunk): + actual = gaussian_smearing(energy, counts, 1., chunk=chunk) + np.testing.assert_allclose(actual, expected, rtol=1e-14, atol=0.) + self.assertTrue(np.isfinite(actual).all()) + self.assertTrue((actual >= 0.).all()) + np.testing.assert_allclose(gaussian_smearing(energy, counts[:, 0], 1., chunk=chunk), + expected[:, 0], rtol=1e-14, atol=0.) + + +if __name__ == "__main__": + unittest.main() From 24bec28d1c7d12ce21812c948316791bb3e0a57d Mon Sep 17 00:00:00 2001 From: ABACUS-PR-BOT <342302420+ABACUS-PR-BOT@users.noreply.github.com> Date: Fri, 16 Oct 2026 07:41:32 +0800 Subject: [PATCH 8/8] Exercise default DOS selection and document helper validation --- tools/02_postprocessing/plot-tools/abacus_plot/utils.py | 4 ++-- tools/02_postprocessing/plot-tools/tests/test_dos.py | 5 ++++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/tools/02_postprocessing/plot-tools/abacus_plot/utils.py b/tools/02_postprocessing/plot-tools/abacus_plot/utils.py index 0559e9d6ab..5cdcd4b426 100644 --- a/tools/02_postprocessing/plot-tools/abacus_plot/utils.py +++ b/tools/02_postprocessing/plot-tools/abacus_plot/utils.py @@ -254,8 +254,8 @@ def gaussian_smearing(energy: Sequence, dos: Sequence, sigma: float, chunk: int :params energy: energy grid in eV :params dos: raw DOS (number of states per energy bin), shape (npoints,) or (npoints, nspin) - :params sigma: Gaussian width parameter in eV - :params chunk: number of output rows computed per block + :params sigma: finite positive Gaussian width parameter in eV + :params chunk: positive integer number of output rows computed per block """ if np.ndim(sigma) != 0 or not np.isfinite(sigma) or sigma <= 0: raise ValueError("sigma must be finite and positive") diff --git a/tools/02_postprocessing/plot-tools/tests/test_dos.py b/tools/02_postprocessing/plot-tools/tests/test_dos.py index 3655ad5b0c..db1107c7cc 100644 --- a/tools/02_postprocessing/plot-tools/tests/test_dos.py +++ b/tools/02_postprocessing/plot-tools/tests/test_dos.py @@ -59,7 +59,10 @@ def test_formats_and_default_raw_selection(self): np.testing.assert_array_equal(tdos.dos_smear, smeared) for smearing in (True, False): fig, ax = plt.subplots() - tdos.plot(fig, ax, smearing=smearing) + if smearing: + tdos.plot(fig, ax) + else: + tdos.plot(fig, ax, smearing=False) expected = smeared if smearing and smeared is not None else raw np.testing.assert_array_equal(ax.lines[0].get_ydata(), expected[:, 0]) if tdos.nspin == 2: