Als Anwendung der Monte-Carlo-Methode auf die statistische Mechanik wird häufig eine Simulation des Spinmagnetismus unter Verwendung eines ansteigenden Modells durchgeführt. ** In diesem Artikel wird die Monte-Carlo-Simulation eines zweidimensionalen Anstiegsmodells mit Python durchgeführt. ** **.
Der Überblick über Theorie und Methode ist in Addendum zusammengefasst. Wenn Sie interessiert sind, lesen Sie ihn bitte. ..
$ N_x \ times N_y $ Der Hamilton-Operator des zweidimensionalen Ising-Modells auf dem Gitter ist unten angegeben [1,2].
** Der Zweck dieser Arbeit ist es, die Temperaturabhängigkeit der Wärmeenergie, Magnetisierung und spezifischen Wärme dieses zweidimensionalen ansteigenden Spinsystems durch Monte-Carlo-Simulation nach der Metropolis-Methode zu bestimmen. </ font> **
[** Berechnungsbedingung **] In der Simulation wurde ein Raster von $ 20 \ mal 20 $ verwendet und $ J = 1 $ und $ h = 0 $. Die Gesamtzahl der Monte-Carlo-Schritte wurde auf 40.000 eingestellt, und die Anzahl der Proben für die Berechnung des thermischen Nahwerts wurde auf ungefähr 8000 eingestellt (Monte-Carlo-Schritte gemittelt von ungefähr 32000).
Darüber hinaus ist die Spin-Spin-Wechselwirkung von Spins an der Grenze Periodische Randbedingung E5% A2% 83% E7% 95% 8C% E6% 9D% A1% E4% BB% B6) wurde zur Bewertung verwendet.
Verwendet Python 3.5.1. Implementierung und Überprüfung des Vorgangs auf Jupyter Notebook 5.0.0.
"""
Monte-Carlo-Simulation nach der Metropolis-Methode des 2D-Rising-Modells
"""
from random import random, randrange
import numpy as np
import matplotlib.pyplot as plt
#Definieren Sie eine Funktion zur Berechnung der Energie für eine bestimmte Spinanordnung
def Ecalc2(alis2, Nx,Ny):
dum=0
for i in range(-1,Nx):
for j in range(0,Ny):
l=i
m=j
if l==-1:
l=Nx
if l == Nx:
l=0
if m ==-1:
m=Ny
if m== Ny:
m=0
ll=i+1
if ll == Nx:
ll=0
dum+=alis2[l, m]*alis2[ll, m]
for i in range(0,Nx):
for j in range(-1,Ny):
l=i
m=j
if l==-1:
l=Nx
if l == Nx:
l=0
if m ==-1:
m=Ny
if m== Ny:
m=0
mm=j+1
if mm == Ny:
mm=0
dum+=alis2[l, m]*alis2[l, mm]
return dum
#Zufällige Erzeugung des Anfangszustands:Die Anzahl der Upspins und die Anzahl der Downspins sind gleich. M.=0
def Initial_rand(s, Nx,Ny):
NN=int(Nx*Ny/2)
for k in range(NN):
i=randrange(Nx-1) #Eine Zufallszahl wird von 1 bis N gewählt
j=randrange(Ny-1) #Eine Zufallszahl wird von 1 bis N gewählt
s[i,j]=-1*s[i,j]
return s
Nx= 20 #In x-Richtung in 100 geteilt
Ny =20 #In y-Richtung in 100 geteilt
Ntot=Nx*Ny
KBT_lis=np.linspace(0.001,6, 201) #Temperatur 0.Von 001 bis 6(In kBT-Einheiten)Es schwankt in Schritten von 201.
ETplot=[]
MTplot=[]
CTplot=[]
for KBT in KBT_lis:
J = 1 #Spin-Kopplungskonstante
B = 0.0 #Externes Magnetfeld
steps =40000 #MC Schritt
avsteps=int(steps/5)
#Generierung des Anfangszustands:Zufällige Spinplatzierung
s= np.ones([Nx,Ny], int) #Quantenzahl für N Spins(Sz = +1 or -1)Als 1 einstellen
s=Initial_rand(s, Nx,Ny)
E = -J* Ecalc2(s, Nx,Ny) -B*np.sum(s) #(Initiale)Energie berechnen.
E2 = E**2 # E^Speichern 2
#installieren
eplot = [] #Definieren Sie eine Liste zum Speichern von Energiewerten unter jeder Temperatur kBT
Eavplot=[] #Definieren Sie eine Liste, um den thermischen Durchschnitt der Energie zu speichern
e2plot = []#Definieren Sie eine Liste, um das Quadrat des Energiewerts unter jeder Temperatur kBT zu speichern. Wird zur Berechnung der spezifischen Wärme verwendet.
eplot.append(E/Ntot)
Eavplot.append(E/Ntot)
e2plot.append((E/Ntot)**2)
Cvplot = [] ##Definieren Sie eine Liste, um spezifische Wärme unter jeder Temperatur kBT zu speichern
M = np.sum(s) #Berechnung der Magnetisierung
Mplot = []
Mavplot=[]
Mplot.append(M/Ntot) #Thermischer Mittelwert der Magnetisierung
Mavplot.append(M/Ntot) #Definieren Sie eine Liste, um den thermischen Durchschnitt der Magnetisierung zu speichern
#Maine
for k in range(steps):
i=randrange(Nx-1) #0 bis N.-Eine Zufallszahl wird bis zu 1 gewählt
j=randrange(Ny-1) #0 bis N.-Eine Zufallszahl wird bis zu 1 gewählt
s_trial=s.copy()
s_trial[i,j]= -1*s[i,j]
delta_E=2*s_trial[i,j]*-1*J*(s[i+1,j]+s[i-1,j]+s[i,j+1]+s[i,j-1])-B*(s_trial[i,j]-s[i,j])
E_trial =E+ delta_E
#Statusaktualisierung nach der Metropolis-Methode
if E_trial < E :
s = s_trial
E = E_trial
else :
if random() < np.exp(-(delta_E)/KBT):
s = s_trial
E = E_trial
eplot.append(E/Ntot)
e2plot.append((E/Ntot)**2)
M = np.sum(s)
Mplot.append(M/Ntot)
ETplot.append(np.sum(eplot[-avsteps:])/avsteps)
MTplot.append(np.sum(Mplot[-avsteps:])/avsteps)
CTplot.append((np.sum(e2plot[-avsteps:])/avsteps-((np.sum(eplot[-avsteps:]))/avsteps)**2)/KBT**2) #Berechnung der spezifischen Wärme
plt.plot(KBT_lis,ETplot)
plt.ylabel("Totla energy per spin")
plt.xlabel("kBT")
plt.show()
plt.plot(KBT_lis,MTplot)
plt.ylabel("Total magnetic moment per spin")
plt.hlines([0], 0, KBT_lis[-1], linestyles="-") # y=Zeichnen Sie eine Linie bei 0.
plt.xlabel("kBT")
plt.show()
plt.plot(KBT_lis,CTplot)
plt.ylabel("Heat capacity per spin")
plt.xlabel("kBT")
plt.show()
Das Folgende ist ein Beispiel dafür, wie sich die Menge mit dem Monte-Carlo-Schritt während der Monte-Carlo-Simulation unter isothermen Bedingungen ändert. Dies kann nicht direkt aus dem Berechnungscode abgerufen werden. Um es anzuzeigen, schreiben Sie die folgende Plotanweisung, ohne die Temperatur im Berechnungscode zu schleifen (für KBT in KBT_lis :).
plt.plot(eplot)
plt.ylabel("Totla energy per spin")
plt.xlabel("kBT")
plt.show()
plt.plot(Mplot)
plt.ylabel("Total magnetic moment per spin")
#plt.hlines([0], 0, steps, linestyles="-") # y=Zeichnen Sie eine Linie bei 0.
plt.xlabel("kBT")
plt.show()
$ k_BT = 1 $. Es ist ersichtlich, dass die Physik mit zunehmender Anzahl von Schritten konvergiert (der Zustand des Erreichens des thermischen Gleichgewichtszustands).
Zeitliche Entwicklung der Energie pro Spin.
Zeitliche Entwicklung der Magnetisierung pro Spin.
Die durch die Simulation erhaltene Zahl ist unten gezeigt. Wenn Sie genauere Ergebnisse erhalten möchten, können Sie die Anzahl der Gitter in der Simulation erhöhen. Der theoretische Wert von $ T_c $, der Phasenübergangstemperatur von ferromagnetisch zu paramagnetisch, beträgt $ T_c = 2,2691 $ [siehe Anhang (4-1) (12)].
Gesamtenergie pro Drehung. Mit steigender Temperatur treten Paare mit entgegengesetzten Spins auf, wodurch die Energie zunimmt.
Temperaturabhängigkeit der Magnetisierung. Die Magnetisierung verschwindet beim Erhitzen (Phasenübergang von ferromagnetischem Material zu paramagnetischem Material). Obwohl es in der Abbildung streng genommen schwer zu erkennen ist, ändert es sich diskontinuierlich bei $ T = T_c $ (Übergang der Sekundärphase % BB% A2% E7% A7% BB # .E7.AC.AC.E4.BA.8C.E7.A8.AE.E7.9B.B8.E8.BB.A2.E7.A7.BB) Gibt es.
Temperaturabhängigkeit der spezifischen Wärme. Der magnetische Übergang ist [Sekundärphasenübergang](https://ja.wikipedia.org/wiki/%E7%9B%B8%E8%BB%A2%E7%A7%BB#.E7.AC.AC.E4 .BA.8C.E7.A8.AE.E7.9B.B8.E8.BB.A2.E7.A7.BB), und die spezifische Wärme wird bei $ T = T_c $ abgeführt. Diese Tendenz ist auch aus dieser Abbildung ersichtlich.
Die Änderung des Spinzustands während der Simulation wird in der folgenden Animation gezeigt. Die anfänglichen Spinzustände waren alle nach oben (Upspin, $ s_ {i, j} = 1 $) und $ k_BT = 30 $. ** Bei dieser Temperatur wird der normale magnetische Zustand stabil und der Spinzustand ändert sich, so dass die Aufwärts- und Abwärtsdrehungen ungefähr gleich sind **.
Änderung des Spinzustands des 100x100-Gitters. ** Pink steht für Upspin ($ s_ {i, j} = 1
Gemäß der quantenstatistischen Mechanik beträgt der thermische Durchschnittswert der physikalischen Größe (beobachtbarer Operator) $ \ hat A $ im thermischen Gleichgewichtszustand bei der Temperatur T $ \ <A > $, wobei der Hamilton-Operator des Systems $ \ hat H $ ist. ,
Gegeben in. Hier ist Tr die diagonale Summe und $ \ beta = 1 / (k_B T) $ ist die Größe, die als Umkehrtemperatur bezeichnet wird.
$ \ Hat ρ $ ist der Density Matrix Operator Und im thermischen Gleichgewichtszustand
Wenn Sie die Anzeige verwenden, die gleichzeitig $ \ hat H $ und $ \ hat A $ diagonalisiert, ist die diagonale Summe
Anruf. $ E_ {i} $ ist der $ i $ -te Energieeigenwert von $ \ hat H $ und $ Z $ ist die Verteilungsfunktion E9% 85% 8D% E9% 96% A2% E6% 95% B0).
Aus Anhang (1) wurde herausgefunden, dass der thermische Mittelwert physikalischer Größen unter einer endlichen Temperatur theoretisch beschrieben werden kann. Dieses theoretische System wird eines der Denkmäler der theoretischen Physik sein. Unter dem Gesichtspunkt der ** tatsächlichen Berechnung ** ist es jedoch im Allgemeinen äußerst schwierig, den thermischen Durchschnittswert zu berechnen. Der Grund ist, dass die Anzahl der Lösungen (Anzahl der Zustände) der Schledinger-Gleichung oft enorm ist und nicht für die Zustände von ** all ** summiert werden kann.
** Daher werden mehrere geeignete Zustände unter Verwendung von Zufallszahlen abgetastet, und $ \ <A > $ wird ungefähr als Durchschnitt von ihnen berechnet (Monte-Carlo-Simulation). ** **. Unter der Annahme, dass die Anzahl der Abtastwerte N ist, werden die auf diese Weise erhaltenen $ \ <A_N > _ {MC} $ zu allen Paaren, da N unendlich wird, wenn die folgenden Bedingungen (1) und (2) erfüllt sind. Es ist bekannt, dass sich die Summe der Winkel allmählich $ \ <A > $ nähert, was genau berechnet wird. Wenn daher ein ausreichend großes N genommen wird, ist es möglich, den thermischen Durchschnittswert der physikalischen Größe genau zu erhalten.
Bedingungen (1): Die Monte-Carlo-Simulation ist ergodisch (2): ** "Detaillierte Ausgleichsbedingungen" erfüllen **
(1) ist, dass sich der Zustand des Systems mit der Zeitentwicklung ab einem beliebigen Anfangszustand ändert, es jedoch möglich ist, alle möglichen Zustände zu durchlaufen ([Ergod-Bedingung](https: //ja.wikipedia). .org / wiki /% E3% 82% A8% E3% 83% AB% E3% 82% B4% E3% 83% BC% E3% 83% 89% E7% 90% 86% E8% AB% 96)).
** (2) ist die Wahrscheinlichkeit, dass der Zustand $ i $ im thermischen Gleichgewichtszustand $ \ pi (i) $ realisiert wird, und die Wahrscheinlichkeit des Übergangs vom Zustand $ i $ in den Zustand $ j $ (Übergangswahrscheinlichkeit) beträgt $ w (i, j). ) $
Die Wahl von $ w (i, j) $ ist beliebig. Bei der Metropolenmethode wird die Übergangswahrscheinlichkeit $ w (i, j) $ vom Zustand i zum Zustand j wie folgt eingestellt.
** Dadurch entstehen nacheinander verschiedene Zustände, und es wurde nachgewiesen, dass die möglichen Zustände des Systems nach einer ausreichenden Zeit der Boltzmann-Verteilung für die kanonische Bevölkerung in der statistischen Dynamik folgen [4]. Die Methode zur Auswahl von $ w (i, j) $ nach der Metropolenmethode hat daher eine hohe Affinität zur statistischen Dynamik und ist in der rechnergestützten statistischen Physik unverzichtbar. ** **.
Erstellen Sie ein beliebiges Spin-Array $ A_k $
Nehmen Sie zufällig den i-ten Spin auf, um eine neue Spinplatzierung $ A_ {k + 1} $ zu generieren
Erstellen Sie $ A_ {trial} mit invertierter Spinrichtung von Spin i (Spin Flip), um eine Testanordnung zu erstellen, und berechnen Sie $ energy $ E_ {tri} $.
Wenn $ E_i \ \ le E_ {trial}
Wenn $ E_i \ \ gt E_ {trial} $, akzeptieren Sie mit Wahrscheinlichkeit $ P_ {tr} = exp [- (\ beta \ (E_ {trial} -E_i))] $ ($ A_ {k + 1} $ = $ A_ {trial} $)
Wo 5 ist
5-1 Generieren Sie eine einheitliche Zufallszahl $ r $$ (0 \ le r \ le 1) $
5-2. Vergleichen Sie die Größenbeziehung zwischen $ r $ und $ P_ {tr} $, und wenn $ P_ {tr} \ ge r $, dann $ A_ {k + 1} $ = $ A_ {trial}
Und. 6. Wiederholen Sie die obigen Schritte 1 bis 5, bis die physikalische Größe nicht zu groß ist (bis das thermische Gleichgewicht erreicht ist). Energie E, Magnetisierung M (gleiches Volumen) spezifische Wärme Cv
Es gibt eine genaue Lösung im zweidimensionalen Ising-Modell ohne Magnetfeld (h = 0).
Für den zweidimensionalen Hamilton-Operator nach Gleichung (1) sind die allmählichen Lösungen von Energie E, spezifischer Wärme Cv und Magnetisierung M bei Erhöhung der Anzahl der Gitterpunkte wie folgt [3].
Der Magnetismus M ist $ T \ lt T_c (= \ frac {J} {0,4407 \ k_B}) $
Wobei j, $ \ kappa $, $ \ kappa '$ sind
(Berechnungen in Python finden Sie in [5]).
Es liegt im Verhalten von "Spin", einem der quantenmechanischen Freiheitsgrade von Elektronen, und ist das Ergebnis einer Spin-Spin-Konkurrenz (Spin-Spin-Wechselwirkung). Daher ist die Quantenmechanik erforderlich, um die Essenz von Magneten zu verstehen [1]. Die maßgebliche Gleichung ist die Schrödinger-Gleichung (in der nicht-relativistischen Grenze). Daher ist (A) ** Um den Magnetismus einer Substanz zu simulieren, ist es notwendig, die Schrödinger-Gleichung unter Berücksichtigung der Spin-Spin-Wechselwirkung zu lösen. ** **.
Die statistische Mechanik ist die Untersuchung der thermodynamischen Eigenschaften von Materialien bei endlichen Temperaturen. Demnach sind quantitative Informationen über den Magnetismus bei endlicher Temperatur (B) **, die alle Lösungen der (polymorphen) Schrödinger-Gleichung mit vielen Drehungen kennen und mit etwas Gewicht gemittelt werden ** Du kannst es haben.
Leider ist es extrem schwierig, (A) und (B) zu machen. Die folgenden drei Punkte sind die Hauptgründe.
** (C) Die Substanz der Spin-Spin-Wechselwirkung wird Austauschwechselwirkung genannt, die aus den statistischen Eigenschaften von Elektronen abgeleitet wird, und die scheinbare funktionelle Form dieser Wechselwirkung ist unbekannt **,
** (D) Es ist schwierig, die polymorphe Schrödinger-Gleichung zu lösen, selbst wenn die funktionelle Form der Wechselwirkung bekannt ist **
** (E) Selbst wenn die Schrödinger-Gleichung gelöst werden kann, ist es sehr schwierig, den thermischen statistischen Mittelwert (Berechnung der Verteilungsfunktion usw.) zu berechnen **
Weil.
Daher wird bei der Untersuchung des Magnetismus die ** Modellinteraktion (Modell Hamiltonian) für (C) und (D) festgelegt, und die Schledinger-Gleichung wird durch numerische Berechnung gelöst. Für (E) wird in diesem Artikel Es wurde oft lange Zeit getan, um den Rechenaufwand zu reduzieren, indem eine ** effiziente Stichprobe ** durchgeführt wurde, wie in behandelt.
Phase durch Mehrkörperinteraktion (elektronische Korrelation) Das Schreiben von Übergängen ist immer noch eine sehr herausfordernde Aufgabe. Das Zing-Modell und das Heisenberg-Modell (zweidimensionales XY-Modell), das es quantenmechanisch handhabt, haben eine wichtige Rolle bei der Diskussion der Thermodynamik (Phasenübergang usw.) magnetischer Materialien zusammen mit der (Quanten-) Monte-Carlo-Methode gespielt [2]. ..
Kei Yoshida, ["Magnetic"](https://www.amazon.co.jp/%E7%A3%81%E6%80%A7-%E8%8A%B3%E7%94%B0-% E5% A5% 8E / dp / 4000054422), Iwanami Publishing, 1991.
Seiji Miyashita, ["Thermische und statistische Dynamik"](https://www.amazon.co.jp/%E7%86%B1-%E7%B5%B1%E8%A8%88%E5%8A % 9B% E5% AD% A6-% E7% 89% A9% E7% 90% 86% E5% AD% A6% E5% 9F% BA% E7% A4% 8E% E3% 82% B7% E3% 83% AA% E3% 83% BC% E3% 82% BA-% E5% AE% AE% E4% B8% 8B-% E7% B2% BE% E4% BA% 8C / dp / 4563023841 / ref = sr_1_6? S = Bücher & dh = UTF8 & qid = 1502402477 & sr = 1-6 & Schlüsselwörter =% E5% AE% AE% E4% B8% 8B% E7% B2% BE% E4% BA% 8C), Bakufukan, 1993.
Akira Morita und Yutaka Okabe, ["Simulation Statistical Physics"](https://www.amazon.co.jp/Windows%E3%83%A6%E3%83%BC%E3%82%B6%E3% 83% BC% E3% 81% AE% E3% 81% 9F% E3% 82% 81% E3% 81% AE-% E3% 82% B7% E3% 83% 9F% E3% 83% A5% E3% 83 % AC% E3% 83% BC% E3% 82% B7% E3% 83% A7% E3% 83% B3% E7% B5% B1% E8% A8% 88% E7% 89% A9% E7% 90% 86 -% E6% A3% AE% E7% 94% B0-% E7% AB% A0 / dp / 4621042033 / ref = sr_1_1? S = Bücher & dh = UTF8 & qid = 1502402499 & sr = 1-1 & Schlüsselwörter =% E3% 82% B7% E3% 83% 9F% E3% 83% A5% E3% 83% AC% E3% 83% BC% E3% 82% B7% E3% 83% A7% E3% 83% B3% E7% B5% B1% E8% A8% 88% E7% 89% A9% E7% 90% 86), Maruzen Co., Ltd., 1996.
Qiitas Artikel von kaityo256 "Detaillierte Ausgleichsbedingungen und Bestimmung der Übergangswahrscheinlichkeit bei der Monte-Carlo-Methode" ist sehr einfach zu verstehen und hilfreich. Es war. Empfehlung.
Qiita-Artikel "[Wissenschaftliche / technische Berechnung durch Python] Liste der Verwendung von (speziellen) Funktionen, die in der Physik unter Verwendung von scipy verwendet werden" von Python Ich schreibe über die Berechnung des perfekten elliptischen Integrals.
*** "Wie Sie vielleicht wissen, verwenden spinbasierte Monte-Carlo-Simulationen im Allgemeinen Cluster-Aktualisierungsmethoden wie die Swendsen-Wang- und Wolff-Algorithmen. Diese Methode ist im Vergleich zu Single-Spin-Flips sehr leistungsfähig. Wie es geht. Ich habe einen einfachen Artikel über Swendsen-Wang geschrieben, also hoffe ich, dass Sie ihn hilfreich finden. ***
http://qiita.com/kaityo256/items/6539261993e282edc5aa "
Recommended Posts