使用最簡單的逐點翻轉法逐點模擬
import time
import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt
def init_state(n):
''' generates a random spin configuration for initial condition'''
# np.random.randint(2, size=(n,n)) 表示產生n*n隨機數,返回為列表
state = 2*np.random.randint(2, size=(n,n))-1
return state
def flipping(grid, beta):
'''monte carlo move using metropolis algorithm '''
n = len(grid)
for i in range(n):
for j in range(n):
# 產生0到n-1的隨機數
a = np.random.randint(0, n)
b = np.random.randint(0, n)
s = grid[a][b]
e = grid[(a+1)%n][b] + grid[a][(b+1)%n] + grid[(a-1)%n][b] + grid[a][(b-1)%n]
cost = 2*s*e
# 如果能量降低接受翻轉
if cost < 0:
s *= -1
# 在0到1產生隨機數，如果概率小於exp(-e/(kt))翻轉
elif rand() < np.exp(-cost*beta):
s *= -1
grid[a][b] = s
return grid
def calculate_energy(grid):
'''energy of a given configuration'''
energy = 0
n = len(grid)
for i in range(n):
for j in range(n):
s = grid[i][j]
e = grid[(i+1)%n][j] + grid[i][(j+1)%n] + grid[(i-1)%n][j] + grid[i][(j-1)%n]
# 負號來自哈密頓量
energy += -e*s
# 最近鄰4個格點
return energy/4
def calculate_magnetic(grid):
'''magnetization of a given configuration'''
mag = np.sum(grid)
return mag
nt = 2**8 # 溫度取點數量
n = 2**4 # 點陣尺寸, n x n
eqsteps = 2**10 # mc方法平衡步數
mcsteps = 2**10 # mc方法計算步數
energy = # 內能
magnetization = # 磁矩
specificheat = # 比熱容/溫度的漲落
susceptibility = # 磁化率/磁矩的漲落
t= np.linspace(1.2,3.8,nt)
t =list(t)
n1 = 1/mcsteps*n*n
n2 = 1/mcsteps**2*n*n
time_start=time.time()
j = 0
for t in t:
# 初始構型
e = 0 # 每一溫度下的能量
m = 0 # 每一溫度下的磁矩
cv = 0 # 每一溫度下的熱容
k = 0 # 每一溫度下的磁化率
config = init_state(n)
# 熱浴，達到平衡態
for i in range(eqsteps):
flipping(config, 1/t)
# 抽樣計算
for i in range(mcsteps):
flipping(config, 1/t)
e = calculate_energy(config)
m= calculate_magnetic(config)
e += e
m += m
cv += e*e
k += m*m
cv = (n1*cv - n2*e*e)/(t*t)
k = (n1*k - n2*m*m)/t
j +=1
if j%10==0:
print("已完成第%d步模擬" %j)
time_end=time.time()
print('totally cost',time_end-time_start)
magnetization = np.array(magnetization)
f = plt.figure(figsize=(18, 10)); # plot the calculated values
sp = f.add_subplot(2, 2, 1 );
plt.plot(t, energy, 'o', color="red");
plt.xlabel("temperature (t)", fontsize=20);
plt.ylabel("energy ", fontsize=20);
sp = f.add_subplot(2, 2, 2 );
plt.plot(t, abs(magnetization), 'o', color="blue");
plt.xlabel("temperature (t)", fontsize=20);
plt.ylabel("magnetization ", fontsize=20);
sp = f.add_subplot(2, 2, 3 );
plt.plot(t, specificheat, 'o', color="red");
plt.xlabel("temperature (t)", fontsize=20);
plt.ylabel("specific heat ", fontsize=20);
sp = f.add_subplot(2, 2, 4 );
plt.plot(t, susceptibility, 'o', color="blue");
plt.xlabel("temperature (t)", fontsize=20);
plt.ylabel("susceptibility", fontsize=20);
Excel必備工具箱 方方格仔
今天無意間發現這個方方格仔，功能 小技巧挺多的，挺好玩的，特整理分享下 希望經常處理資料的朋友可以減少工作時間，提高效率。方方格仔excel工具箱技巧分享 方方格仔官網 方方格仔安裝指引 2 然後在官網裡點選 方方格仔工具箱 6 然後在解壓後的檔案中點選方方格仔的安裝程式 setup結尾 7 在安裝...
方方格仔access 有什麼好用的辦公軟體推薦？
辦公軟體有個人使用和團隊使用之分，個人使用的辦公軟體 it無憂 答主已經說得很全面了 往上翻 這裡我就不再贅述，我這裡給大家來聊一款團隊辦公協作的軟體。在進入到正題之前，我們先想另乙個問題，除了日程安排 請假考勤等一般性需求之外，企業日常辦公協需要解決的最大問題是什麼？有調查結果顯示，企業中普通員工...
數格仔算面積的方法 方格法計算面積 格仔與面積
例題 圖1是某廣告公司為某種商品設計的商標圖案，若圖中每個小正方形的面積都是l，則陰影部分的面積是 a 6 b 6.5 c 7 d 7.5 此題為某地中考試題，我們先不談題目的解答，而來聊一聊這個有點古怪的圖形。有這麼一類圖形，它們是由多邊形和乙個封閉的網格組成，並且。多邊形的所有頂點都在網格的橫 ...