Cálculo numérico de la densidad de estados

Question

Estoy tratando de averiguar la interpretación numérica de la densidad de estados para un sistema fermiónico bajo un potencial periódico.

La ecuación de la densidad de estados es la siguiente

$$ DOS(E) = \sum_{k \in BZ, n} \delta(E-E_n(k)), $$

donde $E_n(k)$ son los valores propios de la matriz hamiltoniana particular que estoy resolviendo. Me gustaría utilizar la aproximación Cauchy/Lorentziana de la función Delta de manera que la primera ecuación se convierta en

$$ DOS(E) = \frac{1}{\pi} \lim_{\epsilon\to 0} \sum_{k \in BZ, n} \frac{\epsilon}{(E-E_n(k))^2+\epsilon^2}. $$

A partir de aquí, estoy confundido sobre cómo interpretar numéricamente la segunda ecuación. Tengo los respectivos valores propios del hamiltoniano, pero no sé cómo obtener el DOS utilizando $E$ . ¿Cómo puedo incluir $E$ en mi código? Discretizar E para mí significa que tomo una determinada ventana de energía alrededor de un determinado valor $E$ pero no sé cómo estructurarlo, o si debe ser un array, una rejilla... o algo más. Si E debe ser una cuadrícula, ¿debe ser una cuadrícula entre los valores mínimos y máximos de los valores propios de energía?

EDITAR: Hola a todos. Después de reflexionar sobre la respuesta de Murali se me ha ocurrido un pseudocódigo bastante malo pero me gustaría saber si voy en la dirección correcta.

Básicamente codifiqué una función para la función delta ampliada lorentziana así:

def delta_l(x):
  return (1/np.pi)*(epsilon/(epsilon**2 + x**2))

def dos(Egrid,Eigen):
    DOS = np.zeros((AllK,1))

    for j in range(Allk):
    DOS[j] = (1/AllK)*sum([delta(Egrid[j]-Eigen[i]) for i in range(np.shape(Eigen)[0])])

  return DOS

Aquí se le da a epsilon un valor de 0,1 sólo para probar. Los vectores propios del Hamiltoniano se obtuvieron introduciendo puntos de la FBZ:

AllK = len(np.arange(0, 1, 0.01)) * len(np.arange(0, 1, 0.01))

E  = np.zeros((AllK,4*n), float)

count = 0
for m in np.arange(0, 1, 0.01):
    for f in np.arange(0, 1, 0.01):
        kx = (m-f) * np.sqrt(3)/2
        ky = (m+f) * 3.0/2 - 1
        E[count] = Hamiltonian(kx*Kmag, ky*Kmag)
        count = count + 1

import pandas as pd
EinBZ = E.flatten()

Así que entonces obtengo todos los valores propios de la FBZ en esta matriz. ¿Voy en la dirección correcta?

Answer 1

Creo que no has entendido bien lo que significa la densidad de estados (DOS). La DOS es una función de densidad de probabilidad (PDF). Como ha señalado Andrew, toma la energía como entrada y devuelve el número de estados para una energía determinada.

No se puede discretizar $E$ ya que no son valores propios de ningún observable. Es el parámetro de entrada, y discretizarlo simplemente no tiene ningún sentido. Los valores de $E_n(k)$ se discretizan ya que son valores propios del hamiltoniano electrónico.

Si considera la primera ecuación de su pregunta,

$$ DOS(E) = \sum_{k \in BZ, n} \delta(E-E_n(k)), $$ Para las energías $E\neq E_n(k), DOS(E) = 0$ . Sin embargo, esto no es lo que se observa en los experimentos. Típicamente vemos algunas DOS finitas para energías no iguales a los estados propios debido a la incertidumbre de Heisenberg. Para tener en cuenta esto, añadimos un pequeño parámetro de ensanchamiento electrónico finito ( $\epsilon$ ) como se muestra en la segunda ecuación de su pregunta.

Para calcular el DOS, se toma E como parámetro que puede tomar cualquier valor y se fija $\epsilon$ , entonces calcula la doble suma sobre la zona de Brillouin (BZ) y las bandas (n) que son los valores propios de su hamiltoniano. La suma sobre BZ es simplemente la suma sobre k puntos en la zona de Brillouin y dividir la suma obtenida por el número total de k puntos. Elige una malla de k puntos razonable y asegúrate de que converge. Echa un vistazo al siguiente enlace ( http://www.iiserpune.ac.in/~smr2626/hands_on/week1/july1/bzsums_mastani.pdf ) si no tienes ni idea de la suma de BZ

Pseudocódigo:

def delta_l(x):
    return delta function(x)

def E(k):
    return Eigen values for Each k

def dos(E): (let us compute for some E value. This is very inefficient way. just writing for your understanding)
    sum = 0 
    for i in kpoints:
        for j in total_number_of_bands:
            sum = sum + delta_l(E-E(i)[j]) #where E(i)[j] is $j^{th}$ eigen value of $i^{th}$ kpoint 
    return sum/N # N is total number of kpoints