Fix TDOS/PDOS energy grid misalignment, unify DOS output file naming, and improve plot-tools TDOS

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.
 

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 |
 

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  :"* ) ]]

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  :"* ) ]]

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.

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)

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 <cmath>
+
 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<double> dos;
-    std::vector<double> ene;
-    std::vector<double> sum_elec;
-    std::vector<double> 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<double> 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<double> dos(npoints, 0.0);
+    std::vector<double> ene(npoints, 0.0);
+    std::vector<double> dos_smear(npoints, 0.0);
+    std::vector<double> 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<double>(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<double>(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]<dos_thr)
+            if (dos_smear[n] < dos_thr)
             {
-                 dos_smear[i]=0.0;
+                dos_smear[n] = 0.0;
             }
-            sum2 += dos_smear[i] * de_ev;
+            sum2 += dos_smear[n] * de_ev;
 
-            ofs_dos << std::setw(13) << ene[i]
-                 << " " << std::setw(12) << dos[i]
-                 << " " << std::setw(12) << sum_elec[i]
-                 << " " << std::setw(14) << dos_smear[i]
+            ofs_dos << std::setw(13) << ene[n]
+                 << " " << std::setw(12) << dos[n]
+                 << " " << std::setw(12) << sum_elec[n]
+                 << " " << std::setw(14) << dos_smear[n]
                  << " " << std::setw(14) << sum2 << std::endl;
         }
     }

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);
     }
 }

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);
     }
 }

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 <string>
+#include <sstream>
 #ifdef __MPI
 #include "mpi.h"
 #endif
@@ -63,8 +64,45 @@ TEST_F(DosTest,Dos)
 		ifs.open(dosp.fa.c_str());
 		std::string str((std::istreambuf_iterator<char>(ifs)),std::istreambuf_iterator<char>());
 		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"));
+    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