[PYTHON] Génération Spinmap du modèle ing par la méthode Metropolis Hastings (algorithme M-H)

introduction

La carte de spin du modèle de ging en treillis carré bidimensionnel a été créée pour chaque température par la méthode de précipitation de la métropole, elle est donc résumée. Le modèle Zing n'est qu'une explication simplifiée, donc si vous voulez en savoir plus, veuillez vous référer au manuel de dynamique statistique. Vous pouvez également créer des cartes tournantes telles que des grilles triangulaires autres que des grilles carrées par cette méthode, alors essayez-la. (Il suffit de changer le changement hamiltonien)

Algorithme de Metropolis-Hastings

La méthode Metropolis Hastings est la méthode Markov Chain Monte Carlo (MCMC) pour obtenir des séquences d'échantillons aléatoires à partir de distributions de probabilités difficiles à échantillonner directement, et est souvent utilisée en statistique et en physique statistique. Cette séquence peut être utilisée pour approximer la distribution, intégrer et calculer les valeurs attendues.

algorithme

Définissez un état initial arbitraire, répétez la mise à jour par la méthode suivante et adoptez l'état final comme état plausible. Supposons qu'un état $ S $ soit adopté. Génère l'état $ \ aigue {S} $ selon la fonction de densité de probabilité $ g (\ aigue {S} \ mid S) $ comme candidat pour la mise à jour. $ g (\ aigue {S} \ mid S) $ est une fonction de distribution de la probabilité que l'état $ S $ passe à l'état $ \ aigu {S} $. Ici, la vraisemblance (probabilité) de l'état est $ P (S) , et les nouvelles fonctions $ a_1 = \ frac {P (\ sharp {S})} {P (S)}, a_2 = \ frac Introduisez {g (S \ mi \ aigu {S})} {g (\ aigu {S} \ mid S)} $ $ a = a_1a_2 $$. $ a_1 $ est le rapport de vraisemblance (rapport de probabilité) et $ a_2 $ est le rapport de la fonction de densité de probabilité. L'état $ S $ utilise $ a $ pour mettre à jour l'état selon les règles suivantes. ・ When $ a \ ge 1 $ Mettez à jour l'état $ S $ à l'état $ \ aigu {S} $ et passez à l'étape suivante. ・ Quand $ a \ lt 1 $ Met à jour l'état $ S $ à l'état $ \ aigu {S} $ avec la probabilité $ a $ et passe à l'étape suivante, et ensuite sans mettre à jour l'état avec la probabilité $ 1-a $ Passez à l'étape de. La mise à jour est répétée jusqu'à ce que l'état ne soit pas complètement mis à jour.

Modèle Ising

Le modèle Rising est un modèle principalement utilisé comme modèle d'un matériau ferromagnétique. (Ce n'est pas un modèle propre aux matériaux ferromagnétiques, mais fonctionne comme un modèle pour divers phénomènes.) Ce modèle est une théorie qui ne traite que de l'interaction entre les points d'intérêt de la grille et leurs points de grille se rejoignant. Bien que ce soit une théorie si simple, c'est un excellent modèle qui peut décrire le phénomène de transition de phase. L'hamiltonien de ce modèle est $ sans force externe H = -J\sum_{\}\sigma_i\sigma_j Obtenu par $. Où J est l'énergie d'interaction d'échange, $ \ <i, j > $ est la combinaison de points de grille adjacents, et $ \ sigma_k $ est la valeur de spin au point de grille k. La vraisemblance du modèle de Zing est exprimée en utilisant l'hamiltonien, où $ H $ pour l'hamiltonien du système et $ \ frac {\ exp (\ frac {H} {k_BT})} {Z_ \ beta lorsque la température est T. Représenté comme} $. En d'autres termes, si le déplacement hamiltonien dû au changement d'état est $ \ Delta H $, alors $ a_1 = \ exp (\ frac {\ Delta H} {k_BT}) $. En outre, si le candidat de mise à jour est un état dans lequel le spin d'un point de réseau sélectionné au hasard est inversé, la probabilité est constante quelle que soit la transition de tout état à tout état. En d'autres termes, puisque la densité de probabilité ne change pas en fonction de la direction de la transition, $ a_2 = 1 $, donc $ a = a_1 = \ exp (\ frac {\ Delta H} {k_BT}) $.

Modèle de ging en treillis carré bidimensionnel

Un modèle dans lequel des électrons sont présents aux points de réseau de l'espace étiré dans un motif de grille est appelé un modèle de ging en réseau carré bidimensionnel. L'hamiltonien de ce modèle est $$, en supposant $ J = 1 $ H = -\sum_{<i,j>}\sigma_i\sigma_j $$ Peut être exprimé comme. (Cela ne change dans aucun réseau.) En utilisant cela, le changement d'hamiltonien lorsque la valeur de spin d'un point $ (k, l) $ est inversée est $. \Delta H = -\sum_{i=0}^{1}(\sigma_{k,l}\sigma_{k+2i-1,l}+\sigma_{k,l}\sigma_{k,l+2i-1}) $ Sera. (Cela dépend de la grille.) $ \ Sigma_ {N, j} = \ sigma_ {0, j}, \ sigma_ {-1, j} = \ sigma_ {N-1, j}, \ sigma_ {i, N} = \ sigma_ {i, 0} et \ sigma_ {i, -1} = \ sigma_ {i, N-1} $ sont fournis. 例として左図のような系を考え、格子点(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})$の簡単な形式で表せる。

la mise en oeuvre

L'implémentation a été réalisée en Python. Le code entier ressemble à ceci:

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)

Tout d'abord, je vais vous expliquer GenerateSpinState.py. Soit la taille et la température du treillis comme des variables d'instance. Il génère également l'état initial de la carte de rotation.

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

Trouvez le changement hamiltonien $ \ Delta H $ lorsque la valeur de spin de tout point de grille $ (x, y) $ est inversée. En tant que caractéristique du modèle Zing, seuls les points de réseau adjacents sont reflétés dans le calcul, il est donc seulement nécessaire d'extraire et de calculer uniquement les points de réseau autour des points de réseau inversés.

# 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

Un point de grille est sélectionné au hasard, et la carte de rotation avec le point inversé est enregistrée sous new_state jusqu'à ce qu'il soit décidé de la mettre à jour ou non. N'oubliez pas «.copy» lorsque vous déplacez «self.state» vers «new_state» comme une carte de rotation avant l'inversion. Sans .copy, la mise à jour de new_state mettra également à jour self.state et la carte rotative continuera à se mettre à jour.

# 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

Après avoir appelé les deux fonctions ci-dessus, l'état est mis à jour en adoptant new_state avec une probabilité de $ \ exp {\ frac {\ Delta {H}} {k_B \ beta}} $. S'il n'est pas adopté, le statut ne sera pas mis à jour. Le nombre de répétitions est fixé à 10000, mais il n'y a pas de base physique ou mathématique, et c'était une bonne valeur, alors je l'ai définie. En conséquence, j'ai obtenu le résultat le plus probable par hasard, alors je l'ai adopté. (Il n'a pas été complètement mis à jour après 1000 fois.)

# 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 qui stocke la température que vous voulez vérifier en premier self.T J'ai créé des cartes tournantes à différentes températures.

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)

Ceci conclut l'explication de GenerateSpinState.py. Ensuite, pour main.py, la taille de la grille et la disposition de la température sont définies, et GenerateSpinState.py est appelé pour obtenir une carte de rotation après la température. Dans ce cas, la taille du réseau est de 36 et l'agencement de température est de 0,1 à 5,5 divisé en 100 parties égales. La séquence de température a été créée en considérant que la température de transition de phase théorique est de 2,26. J'ai créé l'image gif suivante avec create_spinmap_gif. (Carte de rotation adoptée à chaque température)

Les références

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