Numerische Berechnung der Zustandsdichte

Ich versuche, die numerische Interpretation der Zustandsdichte für ein fermionisches System unter einem periodischen Potential herauszufinden.

Die Gleichung für die Zustandsdichte lautet

D Ö S ( E ) = ∑ k ∈ B Z , N δ ( E − E N ( k ) ) ,

Wo E N ( k ) sind die Eigenwerte der bestimmten Hamilton-Matrix, die ich löse. Ich möchte die Cauchy/Lorentzsche Näherung der Delta-Funktion verwenden, so dass die erste Gleichung jetzt wird

D Ö S ( E ) = 1 π lim ϵ → 0 ∑ k ∈ B Z , N ϵ ( E − E N ( k ) ) 2 + ϵ 2 .

Von hier an bin ich verwirrt, wie ich die zweite Gleichung numerisch interpretieren soll. Ich habe die jeweiligen Eigenwerte des Hamilton-Operators, aber ich weiß nicht, wie ich den DOS mit erhalten soll E . Wie füge ich ein E in meinem Code? E zu diskretisieren bedeutet für mich, dass ich um einen bestimmten Wert herum ein bestimmtes Energiefenster ergreife E , aber ich weiß nicht, wie ich es strukturieren soll, oder ob es ein Array, ein Raster ... oder etwas anderes sein muss. Wenn E ein Gitter sein soll, sollte es ein Gitter zwischen den minimalen und maximalen Werten der Energieeigenwerte sein?

EDIT: Hallo zusammen. Nachdem ich über Muralis Antwort nachgedacht habe, bin ich auf einen Pseudocode gekommen, der ziemlich schlecht ist, aber ich würde gerne wissen, ob ich in die richtige Richtung gehe.

Ich habe im Grunde eine Funktion für die erweiterte Lorentz-Delta-Funktion wie folgt codiert:


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

Hier wird Epsilon nur zum Testen ein Wert von 0,1 gegeben. Die Eigenvektoren des Hamiltonoperators wurden durch Eingabe von Punkten der FBZ erhalten:


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()

In diesem Array bekomme ich dann also alle Eigenwerte der FBZ. Gehe ich in die richtige Richtung?

Naiverweise hätte ich gedacht, Sie würden eine DOS-Funktion definieren wollen, die als Eingabe verwendet wird E und gibt die Zustandsdichte bei diesem Wert von zurück E . Dann können Sie beispielsweise die Zustandsdichte gegen die Energie aufzeichnen, indem Sie ein Array von übergeben E Werte zu dieser Funktion.
Sie müssen eine numerische Summierung/Integration durchführen k Und N für jeden Wert von E . Falls wo k eigentlich eine kontinuierliche Variable ist, ist es vielleicht besser, die Integration zuerst analytisch in Bezug auf die Wurzeln durchzuführen E N ( k ∗ ) = E und dann einfach die Beiträge summieren - das würde jede Unsicherheit aufgrund des Wertes vermeiden ϵ .
@RogerVadim k ist in diesem Fall nicht kontinuierlich. Was ich tue, ist, ein Gitter im Impulsraum so zu wählen, dass ich ungefähr 10.000 K-Punkte habe, die einen Teil der Brillouin-Zone abdecken.
@MadLad Ich bin mir nicht sicher, ob ich Sie verstehe: Sie diskretisieren es, um Computerberechnungen durchzuführen, aber anfangs ist es kontinuierlich (weil der Kristall unendlich ist)?
Richtig, tut mir leid, @RogerVadim. Wir beginnen mit einer irreduziblen Zelle im Impulsraum, und dann würden wir natürlich ein Integral in k machen, aber da ich die Numerik mache, nehme ich einfach einige Arrays in k, die den Bereich dieser Zelle abdecken.
Was ich vorgeschlagen habe, ist, ein Integral zu nehmen, bevor Sie Numerik machen - die Delta-Funktion verpflichtet Sie fast dazu, dies zu tun.

Antworten (1)

Ich glaube, Sie haben nicht richtig verstanden, was Zustandsdichte (DOS) bedeutet. DOS ist eine Wahrscheinlichkeitsdichtefunktion (PDF). Wie Andrew betonte, nimmt es Energie als Eingabe und gibt die Anzahl der Zustände für eine gegebene Energie zurück.

Sie können nicht diskretisieren E da sie keine Eigenwerte irgendeiner Observablen sind. Es ist der Eingabeparameter, und es macht einfach keinen Sinn, ihn zu diskretisieren. Die Werte von E N ( k ) sind diskretisiert, da sie Eigenwerte des elektronischen Hamiltonoperators sind.

Wenn Sie die 1. Gleichung in Ihrer Frage berücksichtigen,

D Ö S ( E ) = ∑ k ∈ B Z , N δ ( E − E N ( k ) ) ,
Für Energien E ≠ E N ( k ) , D Ö S ( E ) = 0 . Dies ist jedoch nicht das, was Sie in Experimenten beobachten. Typischerweise sehen wir aufgrund der Heisenberg-Unsicherheit eine endliche DOS für Energien, die nicht den Eigenzuständen entsprechen. Um dies zu berücksichtigen, fügen wir einen kleinen endlichen elektronischen Verbreiterungsparameter hinzu ( ϵ ) wie in der zweiten Gleichung Ihrer Frage gezeigt.

Um den DOS zu berechnen, nehmen Sie E als Parameter, der jeden Wert annehmen und fixieren kann ϵ , berechnen Sie dann die doppelte Summierung über die Brillouin-Zone (BZ) und die Bänder (n), die die Eigenwerte Ihres Hamilton-Operators sind. Das Summieren über BZ ist einfach das Summieren über k Punkte in der Brillouin-Zone und Dividieren der erhaltenen Summe durch die Gesamtzahl von k Punkten. Wählen Sie ein vernünftiges k-Punkt-Gitter und stellen Sie sicher, dass es konvergiert. Werfen Sie einen Blick auf den folgenden Link ( http://www.iiserpune.ac.in/~smr2626/hands_on/week1/july1/bzsums_mastani.pdf ), wenn Sie keine Ahnung von BZ-Summierung haben

Pseudocode:

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
Mir ist bewusst, dass ich über K-Punkte in der Brillouin-Zone und auch über die Bänder summieren muss. "Sie nehmen E als Parameter, der jeden Wert annehmen kann". Damit bin ich meistens verwechselt. Sollte E also ein Array in Inkrementen sein, die von der minimalen Energie zur maximalen Energie im System gehen? Der Lorentzian würde mir in diesem Fall Werte nahe 1 für jede Energie geben, die nahe an einer Bandenergie liegt, und daher wird DOS höher sein, wenn Bänder an bestimmten Kx, Ky-Punkten eine Energie nahe diesem E-Wert haben.
Der Lorentzian ist nicht 1 bei E = E(k). Es explodiert (1/epsilon). Ja, E kann jeden Wert annehmen, nicht nur zwischen Min und Max. Ich werde ein Analogon geben. Betrachten Sie eine 1D-Gaußsche Verteilung (~e^-{kx^2}). Unabhängig von Ihrem x erhalten Sie einen Wert für die Verteilung (was Sie finden, ist die Wahrscheinlichkeitsdichte bei x). Es ist die Verteilung, die eine Funktion von x ist. Forts.
Wenn Sie weiter zum Mittelwert gehen, fallen Ihre Werte sehr schnell ab. Ebenso möchten Sie die Anzahl der Zustände für eine bestimmte Energie. Wenn es keine Zustände mit dieser Energie gibt, ist Ihr DOS 0 (zB wenn Sie E in der Bandlücke wählen, ist Ihr DOS = 0, da es keine Zustände in der Bandlücke gibt). Natürlich werden Ihre Dosen höher sein, wenn sie mit den Energien der Band übereinstimmen.
Ich denke, Sie haben mir eine bessere Vorstellung davon gegeben, wie ich darüber nachdenken und es berechnen kann. Ich werde einiges ausprobieren, danke! Ich werde alles kommentieren, wenn ich Zweifel habe. @Murali
@MadLad: Ich habe einen Pseudocode hinzugefügt, schau dir an. Dies nur zum Verständnis. Sie können einen viel besseren Code schreiben.
Numerische Berechnung der Zustandsdichte
AntwortenHierDatenschutz-Bestimmungen1y, 0v