I am trying to calculate the density of states for a 2D system given certain band structure data organized in a specific format. Below is an example of my procedure:

Beginning with a simple 2 x 2 Hamiltonian in which kx, ky represent the momenta in the x and y directions, respectively:

ε1 = -2*t*Cos[kx]; 
ε2 = -2*t*Cos[ky]; 
ε12 = 0.3*t*Sin[kx]*Sin[ky]; 
t = 1; 
λ = 2.5;   
h = {{ε1, ε12 + I*λ}, {ε12 - I*λ, ε2}}; 

I wrote a loop which diagonalizes the above Hamiltonian at different momenta and calculates the energy (eigenvalues) at each kx/ky pair, appending a 3-vector consisting of a particular pair of momenta and a corresponding energy (i.e {px, py, Max[e1, e1]}) to a list (dos1 or dos2) at each iteration:

n = 100; 
dos1 = {}; 
dos2 = {}; 
For[i = 0, i < n, i++,
  For[j = 0, j < n, j++, 
    px = -Pi + 2*Pi*(i/n) + 2*(Pi/(2*n)); 
    py = -Pi + 2*Pi*(j/n) + 2*(Pi/(2*n)); 
    e1 = Eigenvalues[h /. {kx -> px, ky -> py}][[1]]; 
    e2 = Eigenvalues[h /. {kx -> px, ky -> py}][[2]]; 
    dos1 = Append[dos1, {px, py, Max[e1, e2]}];
    dos2 = Append[dos2, {px, py, Min[e1, e2]}];];]; 

Everything works fine up to here. I can plot the band structure (a ListPlot3D plot of dos1 and dos2 resulting from the previous loop),

I want to calculate and plot the density of states in the following way:

I wish to define a small energy increment (i.e., an interval within dostot[[m, 3]] where dostot = Union[dos1, dos2] and m is just a counter) and for each energy increment over the whole spectrum, count the number of momenta corresponding to the energies in that interval. I'm thinking something like this (counting only x-direction momenta).

dostot = Union[dos1, dos2];     
dosx = {}; 
For[m = 1, m <= 20000, m++,
  de = Interval[dostot[[m,3]] - 0.2, dostot[[m,3]] + 0.2]; 
  For[i = 1, i <= 20000, i++, 
    dosx = Append[dosx, Count[de, dostot[[i,1]]]];];]; 

Obviously this can't work because dostot[[i, 1]] is a kx momentum which is not an energy like dostot[[m, 3]], so how can I count the momentum corresponding to each energy increment de, and how can I avoid appending all the zeroes which result when a momentum is not inside an energy interval associated with a particular m?

The final plot would have energy (in the specified increments) along the horizontal axis and the density of states (number of momenta in each energy range) on the vertical axis.

  • 2
    $\begingroup$ ClearAll is not a command that clears anything. Standing alone as you used it, It is just a symbol without any value. Perhaps you want Clear["Global`*"] $\endgroup$
    – m_goldberg
    Commented Jun 16, 2016 at 23:54

1 Answer 1


I think this is what you're looking for:

dostot = Union[dos1, dos2][[All, 3]];

{emin, emax} = MinMax[dostot];

de = .25;

Histogram[dostot, {emin, emax, de}, "ProbabilityDensity"]


The definition of dostot differs from yours in that I discard all the coordinate values that would only be needed to make 3D plots. The relevant energy is stored as the third component of each entry, so I only keep that part by doing Union[dos1, dos2][[All, 3]]. The rest is the implementation of your counting idea by using Histogram with a manual specification of the bin width de corresponding to the energy increment.

Here is what the plot looks like when the bin width is made smaller:

de = .05;

Histogram[dostot, {emin, emax, de}, "ProbabilityDensity"]


  • $\begingroup$ This is what I was looking for, thank you very much for your excellent answer! $\endgroup$
    – e.g
    Commented Jun 17, 2016 at 0:39
  • $\begingroup$ And how about unidimensional DOS? $\endgroup$ Commented Dec 26, 2018 at 15:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.