Enge Bindungsdichte von Zuständen genauer erhalten

Entschuldigen Sie die späte Antwort, aber hoffentlich kann dies für jemand anderen nützlich sein!

Sie können das Rauschen mit einem elliptischen Integral reduzieren.

$$D(\varepsilon)=\frac{1}{4 \pi ^3 t}\int_{-\pi }^{\pi }d\phi K\left( \sqrt{1-\left(\frac{\varepsilon +2 t \cos\phi }{4 t}\right)^2}\right)$$

Wobei K das vollständige elliptische Integral erster Art ist: http://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html .

Es ist nicht trivial, hierher zu kommen. Und selbst aus diesem Ausdruck muss das Integral numerisch mit Sorgfalt durchgeführt werden (es hat Singularitäten für viele Werte von$\varepsilon$), sollte aber bessere Ergebnisse liefern. Wenn ich sechs Sekunden lang in Mathematica laufe, bekomme ich (mit$t=\frac{1}{2\sqrt{3}}$):

TL;DR

Mit Matlaboder Mathematikakann ich nicht helfen, aber in Pythongibt es eine Implementierung: sc_dos Here$D$ist die halbe Bandbreite$D = 6t$.

import numpy as np
import gftool as gt

eps = np.linspace(-1.2, 1.2, num=6001)
dos = gt.sc_dos(eps, half_bandwidth=1)

Die Auswertung des DOS dauert bei mir ~100 ms.

---

Sie haben bereits die richtigen Ausdrücke angegeben. Wir haben die Green-Funktion

$$G(z) = \frac{1}{N} \sum_{\boldsymbol{k}} \frac{1}{z-\epsilon_{\boldsymbol{k}}}$$ und die (normalisierte) Zustandsdichte (DOS) $$D(\epsilon) = \frac{1}{N} \sum_{\boldsymbol{k}} \delta(\epsilon - \epsilon_{\boldsymbol{k}}) = -\frac{1}{\pi} \Im G(\epsilon+i0^+),$$ Wo$\epsilon$ist eine reelle Energievariable. Die zweite Gleichheit ist Sokhotski–Plemelj .

Eine naive Summe über die Punkte ist extrem anspruchsvoll, da eine ungeheure Anzahl von$\boldsymbol{k}$Punkte sind in 3 Dimensionen erforderlich.

Um die Funktion zu glätten, können wir die Green'sche Funktion auf einer Kontur parallel zur reellen Achse auswerten, die um ein Endliches verschoben ist$\eta>0$in die obere komplexe Halbebene:

$$D(\epsilon) \approx -\frac{1}{\pi} \Im G(\epsilon+i\eta).$$

Je größer wir wählen$\eta$desto reibungsloser wird die Funktion, aber andererseits verlieren wir Funktionen.

Da uns nur die thermodynamische Grenze interessiert$N\rightarrow \infty$, ein klügerer Ansatz als nur Sampling$\boldsymbol{k}$, ist die Summe durch das Integral zu ersetzen. Für Integrale haben wir mehr oder weniger effiziente Algorithmen.

Rechnen wir also

$$G(z) = \int \frac{d^3 k}{{2\pi}^3} \frac{1}{z-\epsilon_{\boldsymbol{k}}}$$ stattdessen. Der$\epsilon_\boldsymbol{k}$, ist für alle symmetrisch$k_{x_i}$:$\epsilon(k_{x_i}) = \epsilon(-k_{x_i})$, es genügt also, über ein Achtel der Brillouin-Zone zu integrieren.

Und schließlich können wir analytische Ergebnisse für die Integrale verwenden. Wir stellen fest, dass wir die 3D-Green-Funktion in Form bekannter Ergebnisse der 1D- und 2D-Green-Funktion ausdrücken können, wie wir sie haben

$$\epsilon_{\boldsymbol{k}} = \epsilon_{k_x, k_y, k_z} = \epsilon^{2D}_{k_x, k_y} - 2t \cos(k_z) = \epsilon^{1D}_{k_x} - 2t[\cos(k_y) + \cos(k_z)]$$ und deshalb $$G(z) = \int \frac{dk_z}{2\pi} G^{2D} (z - 2t\cos(k_z)) = \int \frac{dk_y}{2\pi} \int \frac{dk_z}{2\pi} G^{1D}(z - 2t[\cos(k_y) + \cos(k_z)]).$$

Die eindimensionale Greensche Funktion$G^{1D}(z)$leicht ausgewertet werden kann, ist die zweidimensionale Green'sche Funktion$G^{2D}(z)$kann durch das vollständige elliptische Integral erster Art ausgedrückt werden (das in Standardlehrbüchern zu finden ist). Verwenden$G^{2D}(z)$ist im Grunde das Ergebnis von bRost03 .

Ein sehr kluger Kerl namens Joyce hat sogar einen Ausdruck für gefunden$G^{3D}(z)$im Jahr 1973. Die Gleichungen sind etwas langwierig und kompliziert, daher vermeide ich es, sie hier zu kopieren. Aber wir haben sie in einem PythonModul implementiert gftool>=0.8.0, siehe sc_dos . Dort finden Sie auch die entsprechenden Referenzen.

Ich hatte das gleiche Problem. Ich habe die Formel verwendet$\rho(\epsilon)=\sum_{\vec{k}}\delta(\epsilon - \epsilon(\vec{k}))$die Zustandsdichte numerisch zu berechnen. Ich habe eine Zusammenfassung gemacht$100$ $k-$Werte für jede Komponente und verwendet eine Gaußsche Verteilung mit$\sigma =0.1$für die Delta-Funktion, um das folgende Diagramm zu erhalten. Größer verwenden$\sigma$glättet schließlich die Singularität in der Ableitung bei ungefähr$\epsilon=\pm 2$.

Der Code wurde in C++ geschrieben und lief ungefähr$60$Sekunden.$y$Achse ist$\rho(\epsilon)$,$x$Achse ist$\epsilon$Und$t=1$.

PS: Ich habe die Summierung über die halbe Brillouin-Zone gemacht, was ich für meine Anwendung brauchte.