我已经使用NumPy和Numba的JIT在Python中实现了2D Ising模型:
from timeit import default_timer as timer
import matplotlib.pyplot as plt
import numba as nb
import numpy as np
# TODO for Dict optimization.
# from numba import types
# from numba.typed import Dict
@nb.njit(nogil=True)
def initialstate(N):
'''
Generates a random spin configuration for initial condition
'''
state = np.empty((N,N),dtype=np.int8)
for i in range(N):
for j in range(N):
state[i,j] = 2*np.random.randint(2)-1
return state
@nb.njit(nogil=True)
def mcmove(lattice, beta, N):
'''
Monte Carlo move using Metropolis algorithm
'''
# # TODO* Dict optimization
# dict_param = Dict.empty(
# key_type=types.int64,
# value_type=types.float64,
# )
# dict_param = {cost : np.exp(-cost*beta) for cost in [-8, -4, 0, 4, 8] }
for _ in range(N):
for __ in range(N):
a = np.random.randint(0, N)
b = np.random.randint(0, N)
s = lattice[a, b]
dE = lattice[(a+1)%N,b] + lattice[a,(b+1)%N] + lattice[(a-1)%N,b] + lattice[a,(b-1)%N]
cost = 2*s*dE
if cost < 0:
s *= -1
#TODO* elif np.random.rand() < dict_param[cost]:
elif np.random.rand() < np.exp(-cost*beta):
s *= -1
lattice[a, b] = s
return lattice
@nb.njit(nogil=True)
def calcEnergy(lattice, N):
'''
Energy of a given configuration
'''
energy = 0
for i in range(len(lattice)):
for j in range(len(lattice)):
S = lattice[i,j]
nb = lattice[(i+1)%N, j] + lattice[i,(j+1)%N] + lattice[(i-1)%N, j] + lattice[i,(j-1)%N]
energy += -nb*S
return energy/2
@nb.njit(nogil=True)
def calcMag(lattice):
'''
Magnetization of a given configuration
'''
mag = np.sum(lattice, dtype=np.int32)
return mag
@nb.njit(nogil=True)
def ISING_model(nT, N, burnin, mcSteps):
"""
nT : Number of temperature points.
N : Size of the lattice, N x N.
burnin : Number of MC sweeps for equilibration (Burn-in).
mcSteps : Number of MC sweeps for calculation.
"""
T = np.linspace(1.2, 3.8, nT);
E,M,C,X = np.zeros(nT), np.zeros(nT), np.zeros(nT), np.zeros(nT)
n1, n2 = 1.0/(mcSteps*N*N), 1.0/(mcSteps*mcSteps*N*N)
for temperature in range(nT):
lattice = initialstate(N) # initialise
E1 = M1 = E2 = M2 = 0
iT = 1/T[temperature]
iT2= iT*iT
for _ in range(burnin): # equilibrate
mcmove(lattice, iT, N) # Monte Carlo moves
for _ in range(mcSteps):
mcmove(lattice, iT, N)
Ene = calcEnergy(lattice, N) # calculate the Energy
Mag = calcMag(lattice,) # calculate the Magnetisation
E1 += Ene
M1 += Mag
M2 += Mag*Mag
E2 += Ene*Ene
E[temperature] = n1*E1
M[temperature] = n1*M1
C[temperature] = (n1*E2 - n2*E1*E1)*iT2
X[temperature] = (n1*M2 - n2*M1*M1)*iT
return T,E,M,C,X
def main():
N = 32
start_time = timer()
T,E,M,C,X = ISING_model(nT = 64, N = N, burnin = 8 * 10**4, mcSteps = 8 * 10**4)
end_time = timer()
print("Elapsed time: %g seconds" % (end_time - start_time))
f = plt.figure(figsize=(18, 10)); #
# figure title
f.suptitle(f"Ising Model: 2D Lattice\nSize: {N}x{N}", fontsize=20)
_ = f.add_subplot(2, 2, 1 )
plt.plot(T, E, '-o', color='Blue')
plt.xlabel("Temperature (T)", fontsize=20)
plt.ylabel("Energy ", fontsize=20)
plt.axis('tight')
_ = f.add_subplot(2, 2, 2 )
plt.plot(T, abs(M), '-o', color='Red')
plt.xlabel("Temperature (T)", fontsize=20)
plt.ylabel("Magnetization ", fontsize=20)
plt.axis('tight')
_ = f.add_subplot(2, 2, 3 )
plt.plot(T, C, '-o', color='Green')
plt.xlabel("Temperature (T)", fontsize=20)
plt.ylabel("Specific Heat ", fontsize=20)
plt.axis('tight')
_ = f.add_subplot(2, 2, 4 )
plt.plot(T, X, '-o', color='Black')
plt.xlabel("Temperature (T)", fontsize=20)
plt.ylabel("Susceptibility", fontsize=20)
plt.axis('tight')
plt.show()
if __name__ == '__main__':
main()
当然，这是可行的:
我有两个主要问题:
- 还有什么需要优化的吗？我知道伊辛模型很难模拟，但看下表，我似乎遗漏了一些东西…
lattice size : 32x32
burnin = 8 * 10**4
mcSteps = 8 * 10**4
Simulation time = 365.98 seconds
lattice size : 64x64
burnin = 10**5
mcSteps = 10**5
Simulation time = 1869.58 seconds
- 我try 了另一个基于不使用字典一遍又一遍地计算指数的优化，但在我的测试中，它似乎更慢.我做错了什么？