[** Condition de calcul **] Dans la simulation, une grille de 20 $ \ fois 20 $ a été utilisée, et $ J = 1 $ et $ h = 0 $. Le nombre total de pas de Monte Carlo a été fixé à 40 000, et le nombre d'échantillons pour le calcul de la valeur proche thermique a été fixé à environ 8 000 (pas de Monte Carlo en moyenne d'environ 32 000).

En outre, l'interaction spin-spin des spins à la limite est [Condition aux limites périodique](https://ja.wikipedia.org/wiki/%E5%91%A8%E6%9C%9F%E7%9A%84% E5% A2% 83% E7% 95% 8C% E6% 9D% A1% E4% BB% B6) a été utilisé pour l'évaluation.

** Code de calcul **

Utilise Python 3.5.1. Opération implémentée et confirmée sur le notebook Jupyter 5.0.0.

"""
Simulation de Monte Carlo par la méthode Metropolis du modèle 2D Rising

"""

from random import random, randrange
import numpy as np
import matplotlib.pyplot as plt



#Définir une fonction pour calculer l'énergie pour un arrangement de spin donné alis
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 



#Génération d'état initial aléatoire:Le nombre de upspins et le nombre de downspins sont les mêmes. M=0
def Initial_rand(s, Nx,Ny):
    NN=int(Nx*Ny/2)
    for k in range(NN):
        i=randrange(Nx-1) #Un nombre aléatoire est choisi de 1 à N
        j=randrange(Ny-1) #Un nombre aléatoire est choisi de 1 à N
        s[i,j]=-1*s[i,j]
    
    return s


Nx= 20  #Divisé en 100 dans la direction x
Ny =20  #Divisé en 100 dans la direction y
Ntot=Nx*Ny


KBT_lis=np.linspace(0.001,6, 201) #Température 0.De 001 à 6(En unités kBT)Il fluctue par pas de 201.
ETplot=[]
MTplot=[]
CTplot=[]
for KBT in KBT_lis:
    
    J = 1 #Constante de couplage de rotation
    B = 0.0 #Champ magnétique externe
    steps =40000 #Étape MC
    avsteps=int(steps/5)
    
        
    #Génération de l'état initial:Placement de spin aléatoire
    s= np.ones([Nx,Ny], int)  #Nombre quantique pour N spins(Sz = +1 or -1)Définir comme 1
    s=Initial_rand(s, Nx,Ny)
 
    E = -J* Ecalc2(s, Nx,Ny) -B*np.sum(s) #(initiale)Calculez l'énergie.
    E2 = E**2  # E^Magasin 2
    
    #installer
    eplot = []  #Définir une liste pour stocker les valeurs énergétiques sous chaque température kBT
    Eavplot=[] #Définir une liste pour stocker la moyenne thermique de l'énergie
    
    e2plot = []#Définissez une liste pour stocker le carré de la valeur énergétique sous chaque température kBT. Utilisé pour calculer la chaleur spécifique.
    
    eplot.append(E/Ntot)
    Eavplot.append(E/Ntot)
    
    e2plot.append((E/Ntot)**2)
    
    Cvplot = [] ##Définir une liste pour stocker la chaleur spécifique sous chaque température kBT
    
    
    M = np.sum(s)  #Calcul de l'aimantation
    Mplot = []
    Mavplot=[]
    Mplot.append(M/Ntot) #Valeur moyenne thermique de l'aimantation
    Mavplot.append(M/Ntot) #Définir une liste pour stocker la moyenne thermique de l'aimantation

    #Principale
    for k in range(steps):
        i=randrange(Nx-1) #0 à N-Un nombre aléatoire est choisi jusqu'à 1
        j=randrange(Ny-1) #0 à N-Un nombre aléatoire est choisi jusqu'à 1
    
        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

        #Mise à jour du statut par la méthode Metropolis
        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)  #Calcul de la chaleur spécifique



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=Tracez une ligne à 0.
plt.xlabel("kBT")
plt.show()


plt.plot(KBT_lis,CTplot)

plt.ylabel("Heat capacity per spin")
plt.xlabel("kBT")
plt.show()

** Résultat (1) Evolution temporelle de l'énergie et de l'aimantation **

Ce qui suit est un exemple de la façon dont la quantité change avec l'étape de Monte Carlo pendant la simulation de Monte Carlo dans des conditions isothermes. Cela ne peut pas être obtenu directement à partir du code de calcul. Pour l'afficher, écrivez l'instruction de tracé suivante sans boucler la température dans le code de calcul (pour KBT dans 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=Tracez une ligne à 0.
plt.xlabel("kBT")
plt.show()

$ k_BT = 1 $. On peut voir que la physique converge à mesure que le nombre de pas augmente (état d'atteindre l'état d'équilibre thermique).

Développement temporel de l'énergie par tour.

Evolution temporelle de l'aimantation par spin.

** Résultat (2) Dépendance à la température **

Le chiffre obtenu par la simulation est présenté ci-dessous. Si vous souhaitez obtenir des résultats plus précis, vous devez augmenter le nombre de grilles dans la simulation. La valeur théorique de $ T_c $, la température de transition de phase du ferromagnétique au paramagnétique, est $ T_c = 2,2691 $ [voir Addendum (4-1) (12)].

Énergie totale par tour. À mesure que la température augmente, des paires avec des spins opposés se produisent, augmentant l'énergie.

Dépendance à la température de l'aimantation. La magnétisation disparaît avec le chauffage (transition de phase du matériau ferromagnétique au matériau paramagnétique). Bien que cela soit difficile à voir sur la figure, à proprement parler, il change de manière discontinue à $ T = T_c $ (Transition de phase secondaire. % BB% A2% E7% A7% BB # .E7.AC.AC.E4.BA.8C.E7.A8.AE.E7.9B.B8.E8.BB.A2.E7.A7.BB)) y a-t-il.

Dépendance à la température de la chaleur spécifique. La transition magnétique est [transition de phase secondaire](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), et la chaleur spécifique se dissipe à $ T = T_c $. Cette tendance ressort également de cette figure.

Addenda

** Addendum (0): Evolution du temps de rotation **

Le changement de l'état de rotation pendant la simulation est illustré dans l'animation ci-dessous. Les états de rotation initiaux étaient tous à la hausse (upspin, $ s_ {i, j} = 1 $) et $ k_BT = 30 $. ** À cette température, l'état magnétique normal devient stable, et l'état de spin change de sorte que les spins ascendants et descendants sont à peu près les mêmes **.

Changement de l'état de rotation de la grille 100x100. ** Le rose représente upspin ($ s_ {i, j} = 1 ) et le bleu clair représente downspin ( s_ {i, j} = -1 $). La descente augmente avec le temps **.

** Addendum (1): Théorie **

Selon la mécanique statistique quantique, la valeur moyenne thermique de la quantité physique (opérateur observable) $ \ hat A $ dans l'état d'équilibre thermique de la température T est $ \ <A > $, où l'opérateur de Hamilton du système est $ \ hat H $. ，

\< A\> = Tr (\hat ρ \hat A) = e^{-\beta \hat H \hat A} /Z \tag{2}

Donné dans. Ici, Tr est la somme diagonale, et $ \ beta = 1 / (k_B T) $ est la quantité appelée température inverse.

$ \ Hat ρ $ est l 'opérateur de matrice de densité Et dans l'état d'équilibre thermique \hat ρ = e^{-\beta \hat H} /Z \tag{3} A été utilisé pour être.

Si vous utilisez l'affichage qui diagonale $ \ hat H $ et $ \ hat A $ en même temps, la somme diagonale est

\< A\> =\sum_{i=1}^{n} A_i \exp(\{ -\beta E_{i} \}), \ \ \ Z =\sum_{i=1}^{n} \exp(\{ -\beta E_{i} \})\tag{4}

Appel. $ E_ {i} $ est la $ i $ ème valeur propre d'énergie de $ \ hat H $, et $ Z $ est la Fonction de distribution E9% 85% 8D% E9% 96% A2% E6% 95% B0).

** Addendum (2): Simulation Monte Carlo par la méthode Metropolis **

À partir de l'Addendum (1), il a été trouvé que la valeur moyenne thermique des grandeurs physiques sous une température finie peut être théoriquement décrite. Ce système théorique sera l'un des monuments de la physique théorique. Cependant, du point de vue du ** calcul effectif **, il est généralement extrêmement difficile de calculer la valeur moyenne thermique. La raison en est que le nombre de solutions (nombre d'états) de l'équation de Schledinger est souvent énorme et ne peut pas être additionné pour les états de ** tous **.

** Par conséquent, plusieurs états appropriés sont échantillonnés à l'aide de nombres aléatoires, et $ \ <A > $ est approximativement calculé comme la moyenne d'entre eux (simulation de Monte Carlo). ** ** En supposant que le nombre d'échantillons est N, le $ \ <A_N > _ {MC} $ obtenu de cette manière devient toutes les paires lorsque N devient infini lorsque les conditions suivantes (1) et (2) sont satisfaites. On sait que la somme des angles s'approche progressivement de $ \ <A > $, qui se calcule exactement. Par conséquent, si un N suffisamment grand est pris, il est possible d'obtenir avec précision la valeur moyenne thermique de la grandeur physique.

\lim_{N\to\infty} \< A_N \>_{MC} = \< A\> \tag{5}

conditions (1): la simulation de Monte Carlo est de type ergod (2): ** Satisfaire aux "Conditions d'équilibrage détaillées" **

(1) est que l'état du système change avec l'évolution temporelle à partir d'un état initial arbitraire, mais il est possible de passer par tous les états possibles ([condition d'Ergod](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) définit la probabilité que l'état $ i $ soit réalisé dans l'état d'équilibre thermique $ \ pi (i) $ et la probabilité de transition de l'état $ i $ à l'état $ j $ (probabilité de transition) $ w (i, j). ) $ \pi(i)\ w(i,j)=\pi(j) \ w(j,i) Rencontrer. ** **

Le choix de $ w (i, j) $ est arbitraire. Dans la méthode metropolis, la probabilité de transition $ w (i, j) $ de l'état i à l'état j est définie comme suit.

w(i,j) =$ 1 (E_j \le E_i)\ $, $ \ e^{-\beta \ (E_j-E_i)} \ (E_j > E_i) \tag{6}$

** Cela crée plusieurs états les uns après les autres, et il a été prouvé qu'après un temps suffisant, les états possibles du système suivent la distribution de Boltzmann pour la population canonique en dynamique statistique [4]. La méthode de sélection de $ w (i, j) $ par la méthode de la métropole a donc une forte affinité avec la dynamique statistique et est indispensable en physique statistique computationnelle. ** **

** Addendum (3): Procédure de simulation Monte Carlo par la méthode Metropolis du modèle Zing **

Créez un tableau de spin arbitraire $ A_k $ 2. Ramassez au hasard le i-ème tour pour générer un nouveau placement de rotation $ A_ {k + 1} $ 3. Créez $ A_ {trial} avec la direction de rotation du spin i inversée (spin flip) pour créer un arrangement d'essai, et calculez $ energy $ E_ {tri} $. 4. Si $ E_i \ \ le E_ {essai} , acceptez le placement d'essai ( A_ {k + 1} $ = $ A_ {essai}) $

5. Si $ E_i \ \ gt E_ {essai} $, accepter avec probabilité $ P_ {tr} = exp [- (\ beta \ (E_ {essai} -E_i))] $ ($ A_ {k + 1} $ = $ A_ {essai} $)

Où 5 est 5-1. Générer un nombre aléatoire uniforme $ r $$ (0 \ le r \ le 1) $ 5-2. Comparez la relation de magnitude entre $ r $ et $ P_ {tr} $, et si $ P_ {tr} \ ge r $, alors $ A_ {k + 1} $ = $ A_ {essai} , sinon Si tel est le cas, rejetez le changement d'état (( A_ {k + 1} $ = $ A_k $))

Et. 6. Répétez les étapes 1 à 5 ci-dessus jusqu'à ce que la quantité physique ne soit pas trop élevée (jusqu'à ce que l'équilibre thermique soit atteint). Énergie E, aimantation M, chaleur spécifique Cv (volume égal), respectivement

E= -J \ \sum_{i=1}^{n} s_i s_{i+1}\tag{7} M= -J \ \sum_{i=1}^{n} s_i \tag{8} $C_v =\frac{E^2-<E>^2}{k_B T^2}\tag{9} $ Peut être donné. $ n $ est le nombre de particules. 7. Après avoir atteint l'équilibre thermique, calculez la moyenne des quantités physiques en plusieurs étapes. Le nombre d'étapes dépend de la température et du nombre de particules, mais plus la température est élevée, plus il peut être nécessaire d'étapes (car la fluctuation d'énergie devient plus importante).

** Addendum (4): Autre **

** Addendum (4-1) Solution exacte du modèle 2D Rising **

Il existe une solution exacte dans le modèle Ising bidimensionnel sans champ magnétique (h = 0). Pour l'hamiltonien bidimensionnel donné par l'équation (1), les solutions graduelles de l'énergie E, de la chaleur spécifique Cv et de l'aimantation M lorsque le nombre de points du réseau est augmenté sont les suivantes [3]. \frac{E}{N}=-J\ coth(j)\ [1+\frac{2}{\pi}\kappa'\ K(\kappa)] \tag{10}

\frac{Cv}{Nk_B}=\frac{1}{\pi} \frac{j}{sinh(j)}[2K(\kappa)-\frac{2}{\pi}-cosh^2(j) \ E(\kappa)+\kappa' sinh^2(j)\ K(\kappa)]\tag{11}

Le magnétisme M est $ T \ lt T_c (= \ frac {J} {0.4407 \ k_B}) $

M=[\frac{cosh^2(j)}{sinh^4(j)}\ (sinh^2j-1)]^{\frac{1}{8}} \tag{12} Et dans $ T \ gt T_c $, $ M = 0 $. ** Ce $ T_c $ est la température de transition de phase du ferromagnétique au paramagnétique, ce qui montre simplement que le modèle bidimensionnel peut décrire la transition de phase magnétique. ** Dans une dimension, $ M $ est une fonction continue et il n'y a pas de transition magnétique (même chose, mais cela peut être décrit comme $ T_c $ = 0K $).

Où j, $ \ kappa $, $ \ kappa '$ sont respectivement

j=\frac{2J}{k_BT}, \kappa=2\frac{sinh(j)}{cosh^2(j)}, \kappa'=2\frac{sinh^2(j)-1}{cosh^2(j)} \tag{13} Donné dans. De plus, $ K (\ kappa) $ et $ E (\ kappa ') $ sont respectivement des intégrales elliptiques parfaites de première et de deuxième classe.

{K(\kappa)=\int_0^{\pi/2}{\frac{1}{\sqrt{1-k^2\sin^2\theta}}}d\theta}\tag{13}

{E(k)=\int_0^{\pi/2}{\sqrt{1-k^2\sin^2\theta}}d\theta}\tag{14}

(Voir [5] pour les calculs en Python).

** Addendum (4-2) Matériau magnétique **

C'est dans le comportement du "spin", qui est l'un des degrés de liberté de la mécanique quantique des électrons, et est le résultat de la compétition spin-spin (interaction spin-spin). Par conséquent, la mécanique quantique est nécessaire pour comprendre l'essence des aimants [1]. L'équation dominante est l'équation de Schrödinger (dans la limite non relativiste). Par conséquent, (A) ** Pour simuler le magnétisme d'une substance, il est nécessaire de résoudre l'équation de Schrödinger en considérant l'interaction spin-spin. ** **

La mécanique statistique est l'étude des propriétés thermodynamiques des matériaux à des températures finies. Selon lui, les informations quantitatives sur le magnétisme à température finie sont (B) ** connaissant toutes les solutions de l'équation (polymorphe) de Schrödinger avec de nombreux spins, et faisant la moyenne avec un certain poids ** Tu peux l'avoir.

Malheureusement, il est extrêmement difficile de faire (A) et (B). Les trois points suivants en sont les principales raisons.

** (C) La substance de l'interaction spin-spin est appelée interaction d'échange dérivée des propriétés statistiques des électrons, et la forme fonctionnelle apparente de cette interaction est inconnue **,

** (D) Difficile de résoudre l'équation de Schledinger multi-corps même si la forme fonctionnelle de l'interaction est connue **

** (E) Même si l'équation de Schrödinger peut être résolue, il est très difficile de calculer la valeur moyenne statistique thermique (calcul de la fonction de distribution, etc.) **

Parce que.

Par conséquent, dans l'étude du magnétisme, ** l'interaction du modèle (hamiltonien modèle) est définie pour (C) et (D), et l'équation de Schrödinger en la considérant est résolue par un calcul numérique. Cela a souvent été fait pendant longtemps pour réduire la quantité de calcul en effectuant un ** échantillonnage efficace ** comme indiqué dans.

Phase par interaction multi-corps (corrélation électronique) Décrire la transition est toujours une tâche très difficile. Le modèle de Zing et le modèle de Heisenberg (modèle XY bidimensionnel) qui le gère de manière quantique ont joué un rôle important dans la discussion de la thermodynamique (transition de phase, etc.) des matériaux magnétiques avec la méthode (quantique) de Monte Carlo [2]. ..

Les références

1. 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.

2. Seiji Miyashita, ["Thermal and Statistical Dynamics"](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 = livres & ie = UTF8 & qid = 1502402477 & sr = 1-6 & keywords =% E5% AE% AE% E4% B8% 8B% E7% B2% BE% E4% BA% 8C), Bakufukan, 1993.

3. Akira Morita et 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 = livres & ie = UTF8 & qid = 1502402499 & sr = 1-1 & mots-clés =% 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.

4. L'article de Qiita par kaityo256 "Conditions d'équilibrage détaillées et comment déterminer la probabilité de transition dans la méthode de Monte Carlo" est très facile à comprendre et utile. C'était. recommandation.

5. Article Qiita "[Calcul scientifique / technique par Python] Liste des utilisations des fonctions (spéciales) utilisées en physique en utilisant scipy" par Python J'écris sur le calcul de l'intégrale elliptique parfaite.

Postscript

12 août 2017 Nous avons reçu les précieux commentaires suivants de M. kaityo256. Ce sera utile. Merci beaucoup.

*** "Comme vous le savez peut-être, les simulations de Monte Carlo basées sur le spin utilisent généralement des méthodes de mise à jour de cluster telles que les algorithmes Swendsen-Wang et Wolff. Elles sont super puissantes comparées aux retournements de spin unique. Comment faire. J'ai écrit un article simple sur Swendsen-Wang, donc j'espère que vous le trouverez utile. ***

http://qiita.com/kaityo256/items/6539261993e282edc5aa "