Ising (伊辛)模型為：

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

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

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) # 獲取格點(diǎn)的能量

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

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

plt.show()

def get_one_sample(sizeOfSample, temperature):

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

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

initialEnergy = calculateAllEnergy(S) # 計(jì)算隨機(jī)初始狀態(tài)的能量

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

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

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

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

newEnergy = calculateAllEnergy(newS)

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

print('循環(huán)次數(shù)%i, 當(dāng)前系統(tǒng)能量是:' % (i00+1), newEnergy)

print('循環(huán)次數(shù)%i, 當(dāng)前系統(tǒng)能量是:' % (i00 + 1), newEnergy)

return newS

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

delta_max = 0.5 * np.pi # 角度最大的變化度數(shù)，默認(rèn)是90度，也可以調(diào)整為其他

k = 1 # 玻爾茲曼常數(shù)

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) # 獲取該格點(diǎn)的能量

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

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

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

S[i, j] = newAngle # 接受新狀態(tài)

else:

pass # 保持為上一個(gè)狀態(tài)

return S

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

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): # 計(jì)算整個(gè)體系的能量

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()

計(jì)算結(jié)果，得到某個(gè)樣品(溫度

)：

能量分布為：

收斂情況為：

溫度

的情況，結(jié)果為：

溫度

的情況，結(jié)果為：

長(zhǎng)和寬改為60， 溫度為

情況，結(jié)果為：

把程序修改為只有spin-up, spin-down兩種狀態(tài)，計(jì)算不同溫度下的總磁矩，可得到轉(zhuǎn)變溫度。程序?yàn)?#xff1a;

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)) # 初始狀態(tài)

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ài)， T是溫度

k = 1 # 玻爾茲曼常數(shù)

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) # 獲取該格點(diǎn)的能量

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

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

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

S[i, j] = newAngle # 接受新狀態(tài)

else:

pass # 保持為上一個(gè)狀態(tài)

return S

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

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): # 計(jì)算整個(gè)體系的能量

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()

結(jié)果為(磁矩隨溫度的變化關(guān)系)：

參考資料：

+2

總結(jié)

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

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。