# 簡單高效的2D Ising model 蒙地卡羅程式實現

:::info
此文章主旨為展示如何使用程式進行簡短且高效的Ising model simulation
因此本文將不會對 Ising model 做過多的物理解釋
也不會使用較難的方式如使用 MPI
文中會展示的程式語言有 C++, Python, Julia
:::

## 1. Ising model 簡介
Ising model 是最簡單的相變模型，對於計算物理、統計物理及凝態物理都是較為重要且基本的模型
學習 Ising model simulation 對於瞭解這些領域及物理概念如：鐵磁性物質、相變等 有很大幫助

Ising model 有 Hamiltonian
\begin{equation}
H = E = -J\sum_{<i,j>}s_is_j+h \sum_i s_i
\end{equation}
其中$s$為自旋，有值$s=±1$；$J$為自旋-自旋耦合常數；$h$為外場，本文中$h=0, J=1$
$\sum_{<i,j>}$ 表示對所有點$i$及其鄰近點$j$做求和

在溫度$T$時，Ising model 滿足波茲曼分布(Boltzmann distribution)
\begin{equation}
P(σ)={1 \over Z}exp(-{E(σ) \over k_B T})
\end{equation}
其中$\sigma$為自旋組態，$Z$為配分函數，在模擬中波茲曼常數$k_B=1$

## 2. 模擬方法
在電腦上進行 Ising model 模擬的主要方法為馬可夫鍊蒙地卡羅 (MCMC)，其宗旨是在電腦上建立晶格，並透過 Metropolis 演算法使該晶格"演化"至平衡態，藉由統計在各平衡態上的觀測量，我們可以獲得系統在該溫度下的觀測量。
\
模擬的步驟分為以下幾個部分
* 設定參數：如溫度、晶格大小、演化步數等
* 初始化：初始化晶格、亂數等
* 演化
* 測量物理量：常見物理量有磁化強度、磁化率、熱容等
* 重複演化及測量直至有足夠多的統計數據

### Metroplis 演算法
Metropolis 演算法是一種高效的模擬演化方法，但要注意的是使用此方法的演化過程並沒有物理意義，物理意義僅在統計系統到達熱平衡時的物理量後才出現。
在 Metropolis 演算法中，狀態$\sigma_A$變為狀態$\sigma_B$的機率為
\begin{equation}
P(\sigma_A → \sigma_B)=min[1,exp(-{\Delta E_{AB} \over T})]
\end{equation}
*此方程式意為：若$\Delta E_{AB} < 0$則一定翻轉$(P=1)$，反之則有翻轉機率$P=exp(-{\Delta E_{AB} \over T})$

具體的演化做法為在晶格中選定一點自旋$s_i$(可以照順序選也可以隨機選)，計算其翻轉前後的能量差，由能量差可以得知翻轉機率$P$，取一在$[0,1)$上均勻分布的亂數$r$，若$r<P$則該點翻轉，之後進行下一次演化
\
注意，在一次僅選定一個自旋$s_i$的情況下，由於系統的總能量僅因該點翻轉而改變，因此我們不用考慮整個系統的能量，僅須關注該點造成的影響，此時我們有

$$
\begin{align}
\Delta E_{AB} &= E_B - E_A \\
&= \sum_{<j>}(-s_i)s_j-\sum_{<j>}s_is_j \\
&=-2\sum_{<j>}s_is_j \\
&=-2E_A
\end{align}
$$
則狀態$\sigma_A$變為狀態$\sigma_B$的機率可以寫為
\begin{equation}
P(\sigma_A → \sigma_B)=min[1,exp({2E_A \over T})]
\end{equation}

## 3. 程式碼實例
:::warning
此程式碼僅作簡單計算，結果有一定物理意義，但若要有更好更合理的結果則須考慮測量時系統是否已達到平衡、資料相關性等問題
:::

程式依照 C++, Python, Julia 順序提供範例，程式碼中有註解，在程式碼下會對內容進行詳細補充

*註：以下程式中變數名稱定義如下
$$
\begin{align}
L &= 晶格線性大小(邊長) \\
T &= 溫度 \\
lat &= 晶格本體 \\
i,j &= x,y\ 方向晶格編號 \\
m &= 磁化強度 \\
t &= 目前步數 \\
step &= 總步數 \\
dH &= 翻轉前後能量差 = 2E_A
\end{align}
$$

範例觀測之物理量為磁化強度(magnetization)，定義為
\begin{equation}
m={1 \over N} \sum_{i=0}^{N} s_i
\end{equation}
其中$N$為晶格總大小$L \times L$
### C++ 模擬實作

#include <iostream>
#include <cmath>
int main(){
    #定義參數
    constexpr int L=8;
    constexpr float T=2.0;
    constexpr int step=1000;
    #初始化
    int lat[L][L];
    float m[step];
    for(int i=0;i<L;i++){
        for(int j=0;j<L;j++){
            lat[i][j]=1;
        }
    }
    #開始演化
    for(int t=0;t<step;t++){
        float sum=0.0;
        #每一步演化皆將每個格點訪問一次
        for(int i=0;i<L;i++){
            for(int j=0;j<L;j++){
                #計算能量差(2E_A)
                float dH=-2.0*lat[i][j]*(lat[(i+1)%L][j]+lat[(i==0)?L-1:i-1][j]+lat[i][(j+1)%L]+lat[i][(j==0)?L-1:j-1]);
                float rf = 1.0 * rand() / (RAND_MAX + 1);
                #判斷是否翻轉
                if(rf < exp(dH/T)){ lat[i][j] = -lat[i][j]; }
                #進行測量
                sum+=1.0*lat[i][j];
            }
        }
        m[t]=sum/(L*L);
    }
    return 0;
}

* 使用 constexpr 能有較高的執行效率，減少程式運行時的計算
* `lat[i][j]=1`表示此範例在運行時採用冷啟動(固定自旋為1或-1)，如有需求也可以改用熱啟動(每個點隨機是1或-1)
*利用vector函式庫有助於程式優化，但這裡簡單演示，不過多引入非必要函式庫
* `[(i+1)%L]`及`[(i==0)?L-1:i-1]`其目的為建立週期性邊界條件，詳細運作原理看下面補充
* 由於`rf`範圍為$[0,1)<exp(-{\Delta E_{AB} \over T}), \forall \Delta E_{AB}\leq0$因此在實作中不用將$\Delta E_{AB} < 0$情況特別寫出
* 由於C++作圖不便，這裡沒有將我們想要的資料`m[t]`繪圖及存檔，在實際操作中需將該資料存下並用其他程式分析

#### 補充
為了實現週期性邊界條件，我們假設一晶格線性大小$L$，則其一邊可以表示為
\begin{equation}
[...,L-2,L-1],[0,1,...,L-2,L-1],[0,1,...]
\end{equation}
C++陣列編號由0開始，一個[ ]內為一個周期，可以看出$L-1$的右邊格$L$滿足$L\%L=0$，但因C\+\+中`(-1%L)=-1`無法透過`%`算符構成向左的週期性邊界條件，因此使用三元運算符`a ? b : c`代替，其邏輯上等價於

if(a){
    b;
}else{
    c;
}

即只有$i=0$時有左邊為$L-1$，其餘皆為$i-1$與原本相同

### Python 模擬實作

import numpy as np
import random as rd
import matplotlib.pyplot as plt
#設定參數
L=16
T=2.0
step=1000
#初始化
lat=np.ones((L,L))
m=np.zeros(step)
#開始演化
for t in range(step):
    for i in range(L):
        for j in range(L):
            dH=-2.0*lat[i,j]*(lat[(i+1)%L,j]+lat[(i-1)%L,j]+lat[i,(j+1)%L]+lat[i,(j-1)%L])
            #判斷是否翻轉
            if rd.random() < np.exp(dH/T):
                lat[i,j]=-lat[i,j]
    #進行測量
    m[t]=np.mean(lat)
#將結果繪圖
plt.plot(np.linspace(0,step-1,step),m)
plt.show()

* 利用`lat=np.ones((L,L))`直接初始化冷啟動晶格
* 利用`(i+1)%L`及`(i-1)%L`建立週期性邊界條件，需要注意的是與C++不同，Python中的取餘算符`%`回傳值與除數相同，如`(-1%3)=2`，因此可以直接使用`%`算符實現週期性邊界條件

由於 Python 運算速度較慢，若沒有需要使用已存在的函式庫，或使用平行運算、向量化等方式進行優化，不建議使用 Python 生產數據

### Julia 模擬實作

using Plots
function main(::Val{L}, ::Val{T}) where {L, T}
    #設定參數
    step = 1000
    #初始化
    lat = ones(Int8, L, L)
    m = zeros(step)
    #開始演化
    for t in 1:step
        for i in 1:L, j in 1:L
            dH = 2*lat[i, j]*(lat[i%L+1, j] + lat[(i==1) ? L : i-1, j] + lat[i, j%L+1] + lat[i, (j==1) ? L : j-1])
            #判斷是否翻轉
            rand(Float64) < exp(-dH/T) && (lat[i, j] = -lat[i, j])
        end
        #進行測量
        m[t] = sum(lat)/L^2
    end
    #將結果繪圖
    plot(1:step, m, framestyle=:box, xlabel = "step", ylabel = "m", label = "T = $T", size=(900, 600))
end
#呼叫函式以進行模擬
main(Val(16), Val(2.0))

* 在 Julia 中，將程式包裝為函式有助於編譯器優化，因此寫成函式會帶來極大的效能提升
* 將參數$L$及$T$用`::Val{}`包裝，有助於編譯器優化，能大幅提高效能，但會消耗較多的編譯資源
* 週期化邊界條件邏輯與C++相同，但因 Julia 陣列編號由1開始，因此詳細稍有不同
* 在 Julia 中，`a && b`邏輯等效於`if a; b; end`

## 4. 總結
由於 C++ 較為複雜且有些短處如：繪圖不易、語法困難，因此在初心者使用上若沒有特殊需求，作者強烈推薦使用 Julia 作為主要程式語言。Julia 有與 C++ 相近的運算速度及與 Python 相似的簡單語法，也包含 Plots、Dataframe 等工具，實作上不用為了進行分析而學習第二種程式語言，但由於該語言較新，社群較小、應用不如 C++ 般廣泛，在尋找資源時可能反而不如 C++ 方便