Beside the analytic formula that you are using and Alrik commented on it, you can also use the numerical method that Wang-Landau presented to calculate the DOS (DOI: http://dx.doi.org/10.1103/PhysRevLett.86.2050), The algorithm is easy to implement. The only problem you might face is that you need to use the Hamiltonian matrix to provide the eigenvalues to the algorithm. If you are familiar with tight-binding method, that step is already taken.

Hi Ali, I am very interested by your explanations in this discussion,

because I am in bachelor degree of physics. Actually, the diagonalization,

the hamiltonian tight-binding and its eigenvalues are some new notions to me.

As I am also a new user of python programmation language, I would like you to

help me implementing the Wang-Landau algorithm (WLA), please. I have indeed read

this paper, but I don't really picture their implemention:

- how to define the transition probability and the moves acceptation or rejection?

- do I have to start from an initial energy value to obtain all the density of states values?

or do I have to write a loop for each energy value to get the corresponding density of states value?

If you can please help me, I will be very grateful to you. Below, I started to implement the WLA (in python)

supposing for example the eigenvalues of the Hamiltonian are in the ener array.

#code python

import numpy as np

ener = [

-6.494880428220E-01,

-3.623186537844E-01,

-3.622223934478E-01,

-3.622223934478E-01,

3.983706282684E-01,

3.985907117458E-01,

3.985907117458E-01,

5.990507050047E-01]

nsize = len(ener) # length of ener

emin = ener[0] # minimal eigenvalue

emax = ener[-1] # maximal eigenvalue

x = ener[5]

y = 0.0

gx = 1.0

gy = 1.0

fo = np.e # modification factor

f = fo

eps = np.exp(pow(10,-8)) # convergence factor

dos = np.ones(nsize)

hist = np.zeros(nsize)

k = 0 # iteration

while (f > eps):

k = k + 1 # iteration incrementation

y = x # initial state energy

gy = gx # dos of the initial state energy

x = emin + (emax-emin)*np.random.random() # random walk in energy ---> final state energy

f = np.sqrt(f) # reduction of the modification factor

gx = gx*f # dos of the final state energy

The effective density of states Nc in the conduction band or the valence band Nv is the density of electrons in the conduction band or holes in the valence band when the Fermi level coincides with the conduction band edge Ec or the valence band edge Ev.

For more details please see my book: Book Electronic Devices

If you find this book please recommend it for your colleagues and your students.

Best wishes

You can follow

http://link.springer.com/chapter/10.1007%2F978-3-540-74080-3_10)