diff --git a/docs/advanced/elec_properties/dos.md b/docs/advanced/elec_properties/dos.md index 3b42d6811f..47076b54cd 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/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/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..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 (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 (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/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_dos.cpp b/source/source_io/module_dos/cal_dos.cpp index bb6eba8b60..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,66 +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); - ene.push_back(curr_energy); - 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(" -4.995000 0.500000 0.500000")); + EXPECT_THAT(str, testing::HasSubstr(" -5.000000 0.500000 0.500000")); 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 4ad1091604..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 @@ -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,14 @@ 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(" -4.995000 0.500000 0.500000")); + EXPECT_THAT(str, testing::HasSubstr(" -5.000000 0.500000 0.500000")); EXPECT_THAT(str, testing::HasSubstr(" 18.000000 0.000000 8.000000")); + EXPECT_THAT(str, testing::HasSubstr("2.849588")); ifs.close(); - remove("dos.txt"); + remove("doss1g1.txt"); #ifdef __MPI } #endif @@ -154,7 +155,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 +179,14 @@ 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(" -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")); + EXPECT_THAT(str1, testing::HasSubstr("2.847771")); 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 6fec9b3757..80a5531946 100644 --- a/source/source_io/module_parameter/read_inp_out.cpp +++ b/source/source_io/module_parameter/read_inp_out.cpp @@ -341,7 +341,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 diff --git a/tools/02_postprocessing/plot-tools/README.md b/tools/02_postprocessing/plot-tools/README.md index 768db54648..3ed71aa070 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 | -| _pdosfile_ | `str` | Directory or path prefix containing the text PDOS files `pdoss1_*` (and `pdoss2_*` for spin-polarized runs) 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), 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`. 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 3387b60902..ea849bc88b 100644 --- a/tools/02_postprocessing/plot-tools/abacus_plot/dos.py +++ b/tools/02_postprocessing/plot-tools/abacus_plot/dos.py @@ -120,6 +120,17 @@ 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. 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"]) + else: + self.ax.set_xlabel("Energy(eV)", size=25) # notes if notes: @@ -162,20 +173,81 @@ 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)}") + # 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(): + 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 - :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 +257,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 +328,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"]) @@ -255,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. """ @@ -282,7 +389,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) @@ -384,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) @@ -575,11 +682,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 +716,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 @@ -682,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]} 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..1ee51f9459 100644 --- a/tools/02_postprocessing/plot-tools/abacus_plot/utils.py +++ b/tools/02_postprocessing/plot-tools/abacus_plot/utils.py @@ -237,6 +237,40 @@ def energy_minus_efermi(energy: Sequence, efermi: float) -> np.ndarray: return np.array(energy)-efermi +def gaussian_smearing(energy: Sequence, dos: Sequence, sigma: float, chunk: int = 1024) -> np.ndarray: + """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 + 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 + 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']