مدل آیزینگ#
ویدئوی جلسه
ویدئوی جلسه
مدل آیزینگ یکی از سادهترین و در عین حال قدرتمندترین مدلها در فیزیک آماری است. این مدل برای مطالعه سیستمهای ترمودینامیکی و پدیدههای مرتبط با فازهای مختلف ماده، به ویژه فرومغناطیس، استفاده میشود. در این بخش، به بررسی اهمیت مدل آیزینگ و ارتباط آن با فیزیک فرومغناطیس و دمای بحرانی میپردازیم.
۱. مدل آیزینگ چیست؟#
مدل آیزینگ یک مدل ریاضی است که برای توصیف رفتار سیستمهای متشکل از ذرات با اسپین (چرخش) استفاده میشود. در این مدل:
هر ذره (اسپین) میتواند دو حالت داشته باشد: ( +1 ) (اسپین بالا) یا ( -1 ) (اسپین پایین).
اسپینها بر روی یک شبکهی منظم (مانند شبکهی مربعی یا مکعبی) قرار میگیرند.
انرژی سیستم به صورت زیر تعریف میشود: \( E = -J \sum_{\langle i,j \rangle} \sigma_i.\sigma_j \)
\( \sigma_i \): اسپین ذرهی ( i )-ام.
( J ): ثابت برهمکنش بین اسپینها (اگر ( J > 0 )، برهمکنش فرومغناطیسی است).
مجموع \( \langle i,j \rangle \) نشاندهندهی جمع بر روی اسپینهای همسایه است.
۲. اهمیت مدل آیزینگ در مطالعه سیستمهای ترمودینامیکی#
مدل آیزینگ به عنوان یک مدل ساده، درک عمیقی از رفتار سیستمهای ترمودینامیکی ارائه میدهد. برخی از کاربردهای مهم آن عبارتند از:
الف) مطالعه انتقال فاز#
مدل آیزینگ برای مطالعه انتقال فاز بین فازهای مختلف ماده (مانند فازهای مغناطیسی و غیرمغناطیسی) استفاده میشود.
این مدل نشان میدهد که چگونه سیستم در دماهای مختلف از یک فاز به فاز دیگر منتقل میشود.
ب) رفتار جمعی سیستمها#
مدل آیزینگ رفتار جمعی ذرات در سیستمهای بزرگ را توصیف میکند.
این موضوع در مطالعه پدیدههایی مانند خودسازماندهی و تشکیل الگوها مفید است.
ج) شبیهسازیهای عددی#
مدل آیزینگ به عنوان یک مدل پایه برای شبیهسازیهای عددی در فیزیک آماری و محاسبات پیشرفته استفاده میشود.
۳. مدل آیزینگ و فیزیک فرومغناطیس#
مدل آیزینگ به طور خاص برای مطالعه فرومغناطیس (موادی که در دماهای پایین به طور خودبهخودی مغناطیده میشوند) استفاده میشود. در این زمینه:
الف) برهمکنش فرومغناطیسی#
در مدل آیزینگ، اگر ( J > 0 ) باشد، برهمکنش بین اسپینها فرومغناطیسی است.
این بدان معناست که اسپینها تمایل دارند همجهت شوند و سیستم به طور خودبهخودی مغناطیده شود.
ب) دمای بحرانی (Curie Temperature)#
دمای بحرانی (( T_c )) دمایی است که در آن سیستم از فاز فرومغناطیسی به فاز پارامغناطیسی منتقل میشود.
در دماهای پایینتر از ( T_c )، اسپینها همجهت میشوند و سیستم مغناطیده میشود.
در دماهای بالاتر از ( T_c )، اسپینها به طور تصادفی جهتگیری میکنند و سیستم غیرمغناطیسی میشود.
ج) رفتار بحرانی#
نزدیک به دمای بحرانی، سیستم رفتار بحرانی از خود نشان میدهد.
این رفتار با مقادیر بحرانی مانند نمای بحرانی (Critical Exponents) توصیف میشود.
۴. پیاده سازی مدل آیزینگ#
برای پیاده سازی مدل آیزینگ، کلاس Ising را نوشته ایم.
این کلاس باید توصیف کنندهی یک شبکه از اسپین ها باشد که بهترین روش برای نمایش آن، استفاده از یک ماتریس numpy است.
برای پیاده سازی آن، از ایدهی ارث بری یا inheritence که یکی از قابلیتهای بسیار مهم در برنامه نویسی شیء گرا میباشد بهره میبریم.
۱. ارثبری از np.ndarray و تابع __new__#
الف) ارثبری از np.ndarray#
ارثبری از np.ndarray به این معناست که کلاس شما تمامی ویژگیها و متدهای یک آرایهی NumPy را به ارث میبرد. این کار مزایای زیر را دارد:
استفاده از امکانات NumPy: میتوانید از توابع و عملیاتهای NumPy به طور مستقیم در کلاس خود استفاده کنید.
سادگی و کارایی: نیازی به تعریف مجدد بسیاری از توابع پایهای نیست، که باعث کاهش پیچیدگی کد میشود.
ب) تابع __new__#
تابع __new__ در پایتون برای ایجاد یک نمونه جدید از کلاس استفاده میشود. این تابع قبل از __init__ فراخوانی میشود و مسئول ایجاد شیء جدید است. در کلاس IsingModel، از __new__ برای ایجاد یک آرایهی دوبعدی با مقادیر تصادفی ( 1 ) یا ( -1 ) استفاده میشود. این آرایه سپس به عنوان یک نمونه از کلاس IsingModel بازگردانده میشود.
۲. محاسبهی انرژی#
الف) تابع convolve#
هدف محاسبهی جمع اسپین همسایههای هر کدام از خانههای شبکه است. برای این منظور اگر بتوانیم به وسیلهی یک ابزار در پایتون، ماتریس حالت اسپینها را تبدیل به یک ماتریس جدید کنیم که در هر خانه، به جای مقدار اسپین آن خانه، جمع همسایهها نوشته شود کار محاسباتی مان راحتتر خواهد بود. برای این منظور از متد convolve استفاده میکنیم.
تابع convolve از کتابخانهی scipy.ndimage برای اعمال یک فیلتر (کرنل) روی یک آرایه استفاده میشود. این تابع محاسباتی مانند جمع، ضرب یا میانگینگیری را روی همسایههای هر عنصر آرایه انجام میدهد.
ب) فیلتر همسایگی#
در مدل آیزینگ، فیلتر همسایگی به صورت یک ماتریس \( 3 \times 3 \) تعریف میشود که نشاندهندهی همسایههای هر اسپین است. این ماتریس معمولاً به شکل زیر است:
این فیلتر نشان میدهد که هر اسپین تنها با چهار همسایهی مستقیم خود (بالا، پایین، چپ، راست) در تعامل است.
ج) محاسبهی انرژی#
با اعمال فیلتر همسایگی روی آرایهی اسپینها، تابع convolve مجموع اسپینهای همسایه را برای هر اسپین محاسبه میکند. این کار به صورت برداری انجام میشود و بسیار بهینه است، زیرا نیاز به حلقههای تودرتو را کاهش میدهد.
۳. محاسبهی مغناطش#
به راحتی و با جمع کردن کل اسپینها میتوانیم مغناطش حالت را به دست آوریم.
import numpy as np
from scipy.ndimage import convolve
class IsingModel(np.ndarray):
def __new__(cls , N , *args , **kwargs ):
obj = np.random.choice([1 , -1] , size=[N,N] ).view(cls)
return obj
def __array_finalize__(self , obj):
if obj is None: return
def __init__(self , N):
self.N = N
self._filter = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]])
self._neighbors = convolve(self , self.filter , mode='wrap')
@property
def magnetization(self):
return float(np.sum(self))
@property
def filter(self):
return self._filter
@property
def neighbors(self):
return self._neighbors
@property
def energy(self):
return float(-0.5*np.sum( self*self.neighbors ))
def monte_carlo_step(self, temperature : float ) -> bool:
"""
Perform a single Monte Carlo sweep with sequential spin updates.
Returns the number of accepted spin flips.
"""
i = np.random.randint(0, self.N-1)
j = np.random.randint(0, self.N-1)
# Calculate the energy change if this spin is flipped
delta_E = 2 * self[i, j] * self.neighbors[i, j]
# Apply Metropolis criterion
if delta_E <= 0 or np.random.rand() < np.exp(-delta_E / temperature):
self[i, j] *= -1 # Flip the spin
# Update the n_sum matrix around the (i, j) element
self._neighbors[ i+1 % self.N , j] += 2 * self[i, j]
self._neighbors[ i-1 , j] += 2 * self[i, j]
self._neighbors[ i , j+1 % self.N] += 2 * self[i, j]
self._neighbors[ i , j-1] += 2 * self[i, j]
return True
return False
۴. فراهم کردن هنگرد#
در مجموع، برای یک شبکه اسپینی با ابعاد \(N \times N\) میتوان \(2^{N^2}\) حالت متفاوت در نظر گرفت. از فیزیک آماری میدانیم که در دمای T احتمال رخ دادن هر کدام از این حالات برابر با \(e^{-\frac{E}{kT}}\) میباشد. اگر توان محاسباتی ما اجازهی ساختن تمام این \(2^{N^2}\) حالت را به ما میداد، به سادگی میتوانستیم با محاسبهی رابطهی زیر مقدار متوسط مغناطش در هر دمایی را به دست آوریم
اما متاسفانه تعداد حالات آنقدر زیاد است که انجام این محاسبه ممکن نیست
یک محاسبهی ساده
با یک محاسبهی سرانگشتی، مقدار حافظهی لازم برای نگه داشتن تمام حالات ممکن یک شبکهی اسپینی با سایز \(10\times 10\) را حساب کنید. فرض کنید هر اسپین نیاز به یک بایت دارد.
آیا کامپیوتری میشناسید که این میزان حافظه داشته باشد؟؟
برای حل این مشکل، باید بتوانیم یک زیر مجموعه از تمام حالتها را پیدا کنیم که نمایندهی مناسبی از کل حالات باشند و با متوسط گیری روی آنها کمیت مورد نظر در دمای مورد نظر را به دست آوریم. این کار توسط روش متروپولیس ممکن شده است.
اثبات شرط متروپولیس#
شرط متروپولیس تضمین میکند که سیستم به تعادل ترمودینامیکی میرسد. این شرط به صورت زیر بیان میشود:
الف) تعادل جزئی (Detailed Balance)#
شرط تعادل جزئی بیان میکند که احتمال انتقال از حالت ( A ) به حالت ( B ) با احتمال انتقال از حالت ( B ) به حالت ( A ) متناسب است. به عبارت دیگر:
که در آن ( P(A) ) و ( P(B) ) احتمالات تعادلی حالتهای ( A ) و ( B ) هستند.
ب) احتمال انتقال در روش متروپولیس#
در روش متروپولیس، احتمال انتقال از حالت ( A ) به حالت ( B ) به صورت زیر تعریف میشود:
که در آن \( \Delta E = E_B - E_A \) تغییر انرژی و ( T ) دما است.
ج) اثبات شرط تعادل جزئی#
برای اثبات شرط تعادل جزئی، فرض کنید ( P(A) ) و ( P(B) ) احتمالات تعادلی حالتهای ( A ) و ( B ) باشند. این احتمالات به صورت زیر تعریف میشوند:
که در آن ( Z ) تابع پارتیشن است.
با جایگذاری این مقادیر در شرط تعادل جزئی، داریم:
با سادهسازی، به رابطهی زیر میرسیم:
این رابطه نشان میدهد که شرط متروپولیس شرط تعادل جزئی را برآورده میکند و سیستم به تعادل ترمودینامیکی میرسد.
۵ روش اجرا#
یک حالت رندوم میسازیم
با روش متروپولیس شروع به تغییر این حالت میکنیم تا به یک حالت متعادل برسیم
سپس روش متروپولیس را تکرار میکنیم و کمیت مورد نظر (مغناطش) را متوسط گیری میکنیم
Note
دقت کنید که حالتهایی که طبق شرط متروپولیس تغییر نمیکنند باید دو باره در متوسط گیری محاسبه شوند
این کار را برای هر دما چندین بار انجام میدهیم و متوسط مقادیر به دست آمده را میتوانیم به عنوان مغناطش در دمای مورد نظر گزارش کنیم.
۶ دمای بحرانی#
از فیزیک آماری میدانیم که مدل آیزینگ میتواند دمای بحرانی در تغییر مغناطش را نشان دهد. مقدار دمای بحرانی مدل آیزینگ \(T=2.26\) است. برای مشاهده و شبیه سازی این اثر، دما را در بازهی \((0.1, 4.1)\) تغییر میدهیم و متوسط مغناطش بر حسب دما را رسم میکنیم.
from matplotlib import pyplot as plt
import numpy as np
from concurrent.futures import ThreadPoolExecutor
import tqdm
import tqdm.contrib
import tqdm.contrib.concurrent
N_Size = 10
def simulate_ising_model(T , N=N_Size , steps_to_equilibrium = 2000 , steps_to_average = 10000):
E, M, A = np.zeros(steps_to_average), np.zeros(steps_to_average), np.zeros(steps_to_average)
E_ = np.zeros(steps_to_equilibrium)
model = IsingModel(N)
for i in range(steps_to_equilibrium):
model.monte_carlo_step(T)
E_[i] = model.energy
for i in range(steps_to_average):
res = model.monte_carlo_step(T)
E[i] = model.energy
M[i] = model.magnetization
A[i] = res
return T, E_ , E , M , A
temperatures = np.linspace(0.1, 4.1, 10)
N_Simulations = 20
with ThreadPoolExecutor() as executor:
#results = list(executor.map( simulate_ising_model , tqdm( np.repeat(temperatures, N_Simulations) ) ))
results = list(tqdm.contrib.concurrent.thread_map( simulate_ising_model , np.repeat(temperatures, N_Simulations) ) )
دقت کنید که در کد بالا از موازی سازی برای تسریع محاسبات استفاده شده است.
حال در بخش بعدی از محاسبات انجام شده پلاتهای مورد نظر را رسم میکنیم
M_avg = []
M_stdev = []
allEs = np.concatenate([ r[2] for r in results ])
e_bins = np.linspace(np.min(allEs) , np.max(allEs) , 20)
for T in temperatures:
res = [ r for r in results if r[0] == T ]
Es_ = [ r[1] for r in res ]
Es = [ r[2] for r in res ]
Ms = [ r[3] for r in res ]
As = [ r[4] for r in res ]
fig, axs = plt.subplots(2, 1, sharex=False, figsize=(10, 8))
fig.suptitle(f"T={T:.2f}")
boltzman = np.exp( -( e_bins[:-1] - e_bins[0] ) / T)
def ehist(E, T):
h , b = np.histogram(E , bins=e_bins)
hf = h.astype(float )
hf /= boltzman
return hf
for i, E in enumerate(Es):
axs[0].bar( e_bins[:-1] , ehist(E , T) , width=i/len(Es) , label=f'energy distribution {i}' , alpha=0.5)
axs[0].plot( e_bins[:-1] , ehist( np.concatenate(Es) , T )/N_Simulations , label='average')
#axs[0].plot( e_bins[:-1] , boltzman , label='Boltzmann')
axs[0].set_yscale('log')
axs[0].set_ylabel("# of states")
axs[0].set_xlabel("Energy")
axs[0].set_xlim( np.min(allEs) , np.max(allEs) )
#axs[0].legend()
for i, E_ in enumerate(Es_):
axs[1].plot(E_, label=f'energy over time {i}' , alpha=0.5)
axs[1].set_ylabel("Energy")
axs[1].set_xlabel("Step")
#axs[1].legend()
plt.show()
M_vals = [ np.mean( np.abs( M[::1] ) ) for M in Ms ]
M_avg.append( np.mean(M_vals) )
M_stdev.append( np.std(M_vals) )
#now plot magnetization with error bars
plt.errorbar(temperatures, M_avg , yerr=M_stdev , fmt='o' , label='Magnetization')
plt.legend()
plt.title("Magnetization vs Temperature")
plt.xlabel("Temperature")
plt.ylabel("Magnetization")
plt.show()
/tmp/ipykernel_41036/3420507827.py:19: RuntimeWarning: invalid value encountered in divide
hf /= boltzman
تکلیف#
تکلیف مدل آیزینگ
افزودن میدان مغناطیسی خارجی به مدل آیزینگ
برای افزودن یک میدان مغناطیسی خارجی ( h )، ترم جدیدی به هامیلتونی اضافه میشود:
\( h \): میدان مغناطیسی خارجی.
\( \sum_i S_i \): مجموع تمامی اسپینها در سیستم.
محاسبهی مغناطش متوسط (Magnetization)
مغناطش ( M ) به صورت میانگین مجموع اسپینها تعریف میشود:
که در آن ( N ) تعداد اسپینها در سیستم است.
محاسبهی پذیرفتاری مغناطیسی (Magnetic Susceptibility)
حساسیت مغناطیسی \(\chi\) نشاندهندهی پاسخ سیستم به تغییرات میدان مغناطیسی خارجی است و به صورت زیر تعریف میشود:
این کمیت را میتوان به صورت نوسانات مغناطش نیز بیان کرد:
که در آن:
\( k \): ثابت بولتزمان.
\( T \): دما.
\( \langle M \rangle \): میانگین مغناطش.
\( \langle M^2 \rangle \): میانگین مربع مغناطش.