Ising (伊辛)模型為：

這里要用到Metropolis采樣，可看這篇文章：

代碼主要參照參考資料[1]， 是采用XY Ising模型。自己有做了些改動和注釋，看起來會更容易些。代碼如下：

import random

import matplotlib.pyplot as plt

import numpy as np

import copy

import math

import time

def main():

size = 30 # 體系大小

T = 2 # 溫度

ising = get_one_sample(sizeOfSample=size, temperature=T) # 得到符合玻爾茲曼分布的ising模型

plot(ising) # 畫圖

energy_s = []

for i in range(size):

for j in range(size):

energy_s0 = getEnergy(i=i, j=j, S=ising) # 獲取格點的能量

energy_s = np.append(energy_s, [energy_s0], axis=0)

plt.hist(energy_s, bins=50, density=1, facecolor='red', alpha=0.7) # 畫出格點能量分布 # bins是分布柱子的個數，density是歸一化，后面兩個參數是管顏色的

plt.show()

def get_one_sample(sizeOfSample, temperature):

S = 2 * np.pi * np.random.random(size=(sizeOfSample, sizeOfSample)) # 隨機初始狀態，角度是0到2pi

print('體系大小：', S.shape)

initialEnergy = calculateAllEnergy(S) # 計算隨機初始狀態的能量

print('系統的初始能量是:', initialEnergy)

newS = np.array(copy.deepcopy(S))

for i00 in range(1000): # 循環一定次數，得到平衡的抽樣分布

newS = Metropolis(newS, temperature) # Metropolis方法抽樣，得到玻爾茲曼分布的樣品體系

newEnergy = calculateAllEnergy(newS)

if np.mod(i00, 100) == 0:

print('循環次數%i, 當前系統能量是:' % (i00+1), newEnergy)

print('循環次數%i, 當前系統能量是:' % (i00 + 1), newEnergy)

return newS

def Metropolis(S, T): # S是輸入的初始狀態， T是溫度

delta_max = 0.5 * np.pi # 角度最大的變化度數，默認是90度，也可以調整為其他

k = 1 # 玻爾茲曼常數

for i in range(S.shape[0]):

for j in range(S.shape[0]):

delta = (2 * np.random.random() - 1) * delta_max # 角度變化在-90度到90度之間

newAngle = S[i, j] + delta # 新角度

energyBefore = getEnergy(i=i, j=j, S=S, angle=None) # 獲取該格點的能量

energyLater = getEnergy(i=i, j=j, S=S, angle=newAngle) # 獲取格點變成新角度時的能量

alpha = min(1.0, math.exp(-(energyLater - energyBefore)/(k * T))) # 這個接受率對應的是玻爾茲曼分布

if random.uniform(0, 1) <= alpha:

S[i, j] = newAngle # 接受新狀態

else:

pass # 保持為上一個狀態

return S

def getEnergy(i, j, S, angle=None): # 計算(i,j)位置的能量，為周圍四個的相互能之和

width = S.shape[0]

height = S.shape[1]

top_i = i - 1 if i > 0 else width - 1 # 用到周期性邊界條件

bottom_i = i + 1 if i < (width - 1) else 0

left_j = j - 1 if j > 0 else height - 1

right_j = j + 1 if j < (height - 1) else 0

environment = [[top_i, j], [bottom_i, j], [i, left_j], [i, right_j]]

energy = 0

if angle == None:

for num_i in range(4):

energy += -np.cos(S[i, j] - S[environment[num_i][0], environment[num_i][1]])

else:

for num_i in range(4):

energy += -np.cos(angle - S[environment[num_i][0], environment[num_i][1]])

return energy

def calculateAllEnergy(S): # 計算整個體系的能量

energy = 0

for i in range(S.shape[0]):

for j in range(S.shape[1]):

energy += getEnergy(i, j, S)

return energy/2 # 作用兩次要減半

def plot(S): # 畫圖

X, Y = np.meshgrid(np.arange(0, S.shape[0]), np.arange(0, S.shape[0]))

U = np.cos(S)

V = np.sin(S)

plt.figure()

plt.quiver(X, Y, U, V, units='inches')

plt.show()

if __name__ == '__main__':

main()

計算結果，得到某個樣品(溫度

)：

能量分布為：

收斂情況為：

溫度

的情況，結果為：

溫度

的情況，結果為：

長和寬改為60， 溫度為

情況，結果為：

把程序修改為只有spin-up, spin-down兩種狀態，計算不同溫度下的總磁矩，可得到轉變溫度。程序為：

import random

import matplotlib.pyplot as plt

import numpy as np

import copy

import math

import time

def main():

size = 30 # 體系大小

for T in np.linspace(0.02, 5, 100):

ising, magnetism = get_one_sample(sizeOfSample=size, temperature=T)

print('溫度=', T, ' 磁矩=', magnetism, ' 總能量=', calculateAllEnergy(ising))

plt.plot(T, magnetism, 'o')

plt.show()

def get_one_sample(sizeOfSample, temperature):

newS = np.zeros((sizeOfSample, sizeOfSample)) # 初始狀態

magnetism = 0

for i00 in range(100):

newS = Metropolis(newS, temperature)

magnetism = magnetism + abs(sum(sum(np.cos(newS))))/newS.shape[0]**2

magnetism = magnetism/100

return newS, magnetism

def Metropolis(S, T): # S是輸入的初始狀態， T是溫度

k = 1 # 玻爾茲曼常數

for i in range(S.shape[0]):

for j in range(S.shape[0]):

newAngle = np.random.randint(-1, 1)*np.pi

energyBefore = getEnergy(i=i, j=j, S=S, angle=None) # 獲取該格點的能量

energyLater = getEnergy(i=i, j=j, S=S, angle=newAngle) # 獲取格點變成新角度時的能量

alpha = min(1.0, math.exp(-(energyLater - energyBefore)/(k * T))) # 這個接受率對應的是玻爾茲曼分布

if random.uniform(0, 1) <= alpha:

S[i, j] = newAngle # 接受新狀態

else:

pass # 保持為上一個狀態

return S

def getEnergy(i, j, S, angle=None): # 計算(i,j)位置的能量，為周圍四個的相互能之和

width = S.shape[0]

height = S.shape[1]

top_i = i - 1 if i > 0 else width - 1 # 用到周期性邊界條件

bottom_i = i + 1 if i < (width - 1) else 0

left_j = j - 1 if j > 0 else height - 1

right_j = j + 1 if j < (height - 1) else 0

environment = [[top_i, j], [bottom_i, j], [i, left_j], [i, right_j]]

energy = 0

if angle == None:

for num_i in range(4):

energy += -np.cos(S[i, j] - S[environment[num_i][0], environment[num_i][1]])

else:

for num_i in range(4):

energy += -np.cos(angle - S[environment[num_i][0], environment[num_i][1]])

return energy

def calculateAllEnergy(S): # 計算整個體系的能量

energy = 0

for i in range(S.shape[0]):

for j in range(S.shape[1]):

energy += getEnergy(i, j, S)

return energy/2 # 作用兩次要減半

if __name__ == '__main__':

main()

結果為(磁矩隨溫度的變化關系)：

參考資料：

+2

總結

以上是生活随笔為你收集整理的python蒙特卡洛仿真_蒙特卡洛模拟Ising模型（附Python代码）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。