[PYTHON] Spinmap-Generierung eines Modells nach der Metropolis-Hastings-Methode (M-H-Algorithmus)

Einführung

Die Spin-Karte des zweidimensionalen quadratischen Gitter-Ging-Modells wurde für jede Temperatur nach der Metropolis-Hasting-Methode erstellt und daher zusammengefasst. Das Zing-Modell ist nur eine vereinfachte Erklärung. Wenn Sie mehr wissen möchten, lesen Sie bitte das Lehrbuch der statistischen Dynamik. Mit dieser Methode können Sie auch Spin-Maps wie andere Dreiecksgitter als quadratische Gitter erstellen. Probieren Sie es also aus. (Ändern Sie einfach die Hamilton-Änderung)

Metropolis-Hastings-Algorithmus

Die Metropolis Hastings-Methode ist die Markov-Ketten-Monte-Carlo-Methode (MCMC) zum Erhalten von Sequenzen von Zufallsstichproben aus Wahrscheinlichkeitsverteilungen, die schwer direkt abzutasten sind, und wird häufig in der Statistik und der statistischen Physik verwendet. Diese Sequenz kann verwendet werden, um Verteilungen zu approximieren, erwartete Werte zu integrieren und zu berechnen.

Algorithmus

Legen Sie einen beliebigen Anfangszustand fest, wiederholen Sie die Aktualisierung mit der folgenden Methode und übernehmen Sie den Endzustand als plausiblen Zustand. Angenommen, ein Zustand $ S $ wird angenommen. Generieren Sie den Zustand $ \ akut {S} $ gemäß der Wahrscheinlichkeitsdichtefunktion $ g (\ akut {S} \ mid S) $ als Kandidat für die Aktualisierung. $ g (\ akut {S} \ mid S) $ ist eine Verteilungsfunktion der Wahrscheinlichkeit, dass sich der Zustand $ S $ in den Zustand $ \ akut {S} $ bewegt. Hier ist die Wahrscheinlichkeit des Zustands $ P (S) $ und die neuen Funktionen $ a_1 = \ frac {P (\ akut {S})} {P (S)}, a_2 = \ frac Führen Sie {g (S \ mid \ akut {S})} {g (\ akut {S} \ mid S)} $ $ a = a_1a_2 $ ein. $ a_1 $ ist das Wahrscheinlichkeitsverhältnis (Wahrscheinlichkeitsverhältnis) und $ a_2 $ ist das Verhältnis der Wahrscheinlichkeitsdichtefunktion. Der Status $ S $ verwendet $ a $, um den Status gemäß den folgenden Regeln zu aktualisieren. ・ Wenn $ a \ ge 1 $ Aktualisieren Sie den Status $ S $ auf den Status $ \ akut {S} $ und fahren Sie mit dem nächsten Schritt fort. ・ Wenn $ a \ lt 1 $ den Zustand $ S $ mit der Wahrscheinlichkeit $ a $ auf den Zustand $ \ akut {S} $ aktualisieren und mit dem nächsten Schritt fortfahren, und als nächstes, ohne den Zustand mit der Wahrscheinlichkeit $ 1-a $ zu aktualisieren Fahren Sie mit dem Schritt von fort. Die Aktualisierung wird wiederholt, bis der Status nicht vollständig aktualisiert ist.

Ising Modell

Das Rising-Modell ist ein Modell, das hauptsächlich als Modell eines ferromagnetischen Materials verwendet wird. (Es ist kein Modell, das ferromagnetischen Materialien eigen ist, sondern dient als Modell für verschiedene Phänomene.) Dieses Modell ist eine Theorie, die sich nur mit der Wechselwirkung zwischen den interessierenden Gitterpunkten und ihren wiederverbindenden Gitterpunkten befasst. Obwohl es sich um eine so einfache Theorie handelt, ist es ein hervorragendes Modell, das das Phänomen des Phasenübergangs beschreiben kann. Der Hamiltonianer dieses Modells ist $ ohne äußere Kraft H = -J\sum_{\}\sigma_i\sigma_j Erhalten von $. Wobei J die Wechselwirkungsenergie ist, $ \ <i, j > $ die Kombination benachbarter Gitterpunkte ist und $ \ sigma_k $ der Spinwert am Gitterpunkt k ist. Die Wahrscheinlichkeit des Zing-Modells wird unter Verwendung des Hamilton-Operators ausgedrückt, wobei $ H $ für den Hamilton-Operator des Systems und $ \ frac {\ exp (\ frac {H} {k_BT})} {Z_ \ beta, wenn die Temperatur T ist. Dargestellt als} $. Mit anderen Worten, wenn die Hamilton-Verschiebung aufgrund der Zustandsänderung $ \ Delta H $ ist, dann ist $ a_1 = \ exp (\ frac {\ Delta H} {k_BT}) $. Wenn der Aktualisierungskandidat ein Zustand ist, in dem der Spin eines zufällig ausgewählten Gitterpunkts invertiert ist, ist die Wahrscheinlichkeit unabhängig vom Übergang von einem Zustand zu einem Zustand konstant. Mit anderen Worten, da sich die Wahrscheinlichkeitsdichte in Abhängigkeit von der Richtung des Übergangs nicht ändert, ist $ a_2 = 1 $, also $ a = a_1 = \ exp (\ frac {\ Delta H} {k_BT}) $.

Zweidimensionales quadratisches Gittermodell

Ein Modell, in dem Elektronen an den Gitterpunkten in einem gitterartigen Raum vorhanden sind, wird als zweidimensionales quadratisches Gittermodell bezeichnet. Hamiltonian dieses Modells ist $$, vorausgesetzt $ J = 1 $ H = -\sum_{<i,j>}\sigma_i\sigma_j $$ Kann ausgedrückt werden als. (Dies ändert sich in keinem Gitter.) Wenn Sie dies verwenden, beträgt die Änderung in Hamilton, wenn der Spinwert eines Punktes $ (k, l) $ umgekehrt wird, $. \Delta H = -\sum_{i=0}^{1}(\sigma_{k,l}\sigma_{k+2i-1,l}+\sigma_{k,l}\sigma_{k,l+2i-1}) $ Wird sein. (Dies hängt vom Raster ab.) $ \ Sigma_ {N, j} = \ sigma_ {0, j}, \ sigma_ {-1, j} = \ sigma_ {N-1, j}, \ sigma_ {i, N} = \ sigma_ {i, 0} und \ sigma_ {i, -1} = \ sigma_ {i, N-1} $ werden bereitgestellt. 例として左図のような系を考え、格子点(i,j)を反転させたときハミルトニアンはどれほど変化するかを求める。先述の通りイジングモデルでは最隣接格子点のみ考えるので、格子点(i,j)が反転することでスピンの商が変わる格子点の組み合わせは(i,j)と(i+1,j)、(i,j)と((i-1,j)、(i,j)と((i,j+1)、(i,j)と((i,j-1)それぞれ4つのみである。幸いスピンの値は1、-1のみを取るので、ハミルトニアン変化は$\Delta{H}=2(\sigma_{i,j}\sigma_{i+1,j}+\sigma_{i,j}\sigma_{i-1,j}+\sigma_{i,j}\sigma_{i,j+1}+\sigma_{i,j}\sigma_{i,j-1})$の簡単な形式で表せる。

Implementierung

Die Implementierung erfolgte in Python. Der gesamte Code sieht folgendermaßen aus:

GenerateSpinState.py

import numpy as np

class GSS:
    def __init__(self, N, T):
        # lattice size
        self.N = N
        # temperature list
        self.T = T
        # generate N*N spin state by random
        self.state = 2 * np.random.randint(0, 2, (self.N, self.N)) - 1

    # calc hamiltonian displacement　which is square lattice, after flipping one site(x,y)
    def calc_Hamiltonian_square(self, x, y):
        delta_H = 2 * self.state[x, y] * (self.state[(x+1)%self.N, y] + self.state[x, (y+1)%self.N] + self.state[(x-1)%self.N, y] + self.state[x, (y-1)%self.N])
        return delta_H

    # flip random spin site
    def flip_spin(self):
        new_state = self.state.copy()
        [x, y] = np.random.randint(0, self.N, (2))
        new_state[x, y] *= -1
        return new_state, x, y

    # calc specious spin state by metropolis method
    def calc_spin_state(self, t):
        n = 10000
        for i in range(n):
            # get flip site
            new_state, x, y = self.flip_spin()
            # calc hamiltonian displacement
            deltaH = self.calc_Hamiltonian_square(x, y)
            # decide whether to adopt
            if np.random.rand() <= np.exp(-deltaH/t):
                self.state = new_state

    def calc_each_temperature(self):
        # save spin state
        X = []
        for t in self.T:
            # init spin state
            self.state = 2 * np.random.randint(0, 2, (self.N, self.N)) - 1
            # generate spin state which temperature t
            self.calc_spin_state(t)
            # append generate spin state which temperature t
            X.append(self.state)
        return np.array(X)

main.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation

import GenerateSpinState

# save path
# train test data
npy_path = "./result//"
# spin map
spin_map_path = "./result/"


# lattice size
N = 32

# number of data (square lattice)
count = 100

# temperaturelist/Significant digits is 3
T_list = np.linspace(0.1, 5.5, count).round(3)

# make gif which write spin state transition
def create_spinmap_gif(data, T, save_path):
    save_path = save_path + "squarelattice.gif"
    fig = plt.figure()

    def Plot(num):
        plt.cla()
        plt.imshow(data[num])
        t_len = len(str(T[num]))
        title = "temperature:" + str(T[num]) + "0"*(5-t_len)
        plt.title(title)

    anim = FuncAnimation(fig, Plot, frames=100)
    anim.save(save_path, writer='imagemagick')


# call module GSS
GSS = GenerateSpinState.GSS(N, T_list)

# create spin map
X_train = GSS.calc_each_temperature()

# create gif
create_spinmap_gif(X_train, T_list, spin_map_path)

Zuerst werde ich GenerateSpinState.py erklären. Die Größe und Temperatur des Gitters seien Instanzvariablen. Es erzeugt auch den Anfangszustand der Spin-Map.

def __init__(self, N, T):
    # lattice size
    self.N = N
    # temperature list
    self.T = T
    # generate N*N spin state by random
    self.state = 2 * np.random.randint(0, 2, (self.N, self.N)) - 1

Finden Sie die Hamilton-Änderung $ \ Delta H $, wenn der Spinwert eines beliebigen Gitterpunkts $ (x, y) $ invertiert ist. Als Merkmal des Zing-Modells werden nur die wieder benachbarten Gitterpunkte in die Berechnung einbezogen, sodass nur die Gitterpunkte um die invertierten Gitterpunkte extrahiert und berechnet werden müssen.

# calc hamiltonian displacement　which is square lattice, after flipping one site(x,y)
def calc_Hamiltonian_square(self, x, y):
    delta_H = 2 * self.state[x, y] * (self.state[(x+1)%self.N, y] + self.state[x, (y+1)%self.N] + self.state[(x-1)%self.N, y] + self.state[x, (y-1)%self.N])
    return delta_H

Ein Gitterpunkt wird zufällig ausgewählt und die Spin-Map mit dem invertierten Punkt wird als "new_state" gespeichert, bis entschieden wird, ob sie aktualisiert werden soll oder nicht. Vergessen Sie nicht ".copy", wenn Sie "self.state" vor der Inversion als Spinmap nach "new_state" verschieben. Ohne .copy aktualisiert das Aktualisieren von new_state auch self.state und die Spinmap wird weiterhin aktualisiert.

# flip random spin site
def flip_spin(self):
    new_state = self.state.copy()
    [x, y] = np.random.randint(0, self.N, (2))
    new_state[x, y] *= -1
    return new_state, x, y

Nach dem Aufrufen der beiden oben genannten Funktionen wird der Status aktualisiert, indem "new_state" mit einer Wahrscheinlichkeit von $ \ exp {\ frac {\ Delta {H}} {k_B \ beta}} $ übernommen wird. Wenn es nicht übernommen wird, wird der Status nicht aktualisiert. Die Anzahl der Wiederholungen ist auf 10000 eingestellt, aber es gibt keine physikalische oder mathematische Grundlage, und es war ein guter Wert, also habe ich ihn eingestellt. Infolgedessen habe ich das wahrscheinlichste Ergebnis zufällig erhalten, also habe ich es übernommen. (Es wurde nach 1000 Mal nicht vollständig aktualisiert.)

# calc specious spin state by metropolis method
def calc_spin_state(self, t):
    n = 10000
    for i in range(n):
        # get flip site
        new_state, x, y = self.flip_spin()
        # calc hamiltonian displacement
        deltaH = self.calc_Hamiltonian_square(x, y)
        # decide whether to adopt
        if np.random.rand() <= np.exp(-deltaH/t):
            self.state = new_state

Liste, in der die Temperatur gespeichert ist, die Sie zuerst überprüfen möchten. Ich habe Spin-Maps bei verschiedenen Temperaturen erstellt.

def calc_each_temperature(self):
    # save spin state
    X = []
    for t in self.T:
        # init spin state
        self.state = 2 * np.random.randint(0, 2, (self.N, self.N)) - 1
        # generate spin state which temperature t
        self.calc_spin_state(t)
        # append generate spin state which temperature t
        X.append(self.state)
    return np.array(X)

Damit ist die Erläuterung von GenerateSpinState.py abgeschlossen. Als nächstes werden für main.py die Gittergröße und das Temperaturarray festgelegt, und GenerateSpinState.py wird aufgerufen, um nach der Temperatur eine Spin-Map zu erhalten. In diesem Fall beträgt die Größe des Gitters 36 und die Temperaturanordnung 0,1 bis 5,5, aufgeteilt in 100 gleiche Teile. Die Temperatursequenz wurde unter Berücksichtigung der theoretischen Phasenübergangstemperatur von 2,26 erstellt. Ich habe das folgende GIF-Bild mit create_spinmap_gif erstellt. (Spin Map bei jeder Temperatur übernommen)

Verweise

Metropolis–Hastings algorithm Hastings, W.K. (1970). "Monte Carlo Sampling Methods Using Markov Chains and Their Applications".