# Cria uma tabela de referência com o mesmo tamanho
# da matriz de estado. Essa matriz terá, para uma partícula de
# índice 1, uma lista com os índices referente á um passo para cada
# lado, usado para acelerar a simulação.
@njit
def genLookup(side):
index = np.arange(side**2)
up = (index - side) % (side**2)
down = (index + side) % (side**2)
left = np.where(index % side == 0, index - 1 + side, index - 1)
right = np.where((index + 1) % side == 0, index + 1 - side, index + 1)
return np.stack((up, down, left, right))
# Função que simula a mudança de um número predefinido de spins aleatoriamente
# Aceita a mudança se ela diminuir a energia, recusa se ela
# não ultrapassar o parâmetro de corte.
@njit
def timeSteps(state, lookup, steps = 1, J = 1, D = 1, T = 1):
N = state.size
varE = 0.0
# Repete a função para N passos
for i in range(steps):
# A cada passo, escolhe uma partícula e um spin aleatórios
index = np.random.randint(0, N)
oldSpin = state[index]
newSpin = np.random.randint(-1, 2)
# Pula se os spins forem iguais
if newSpin == oldSpin:
continue
# Soma o estado dos vizinhos através da tabela de referência
neiSum = state[lookup[0][index]] +\
state[lookup[1][index]] +\
state[lookup[2][index]] +\
state[lookup[3][index]]
# Calcula a mudança de energia do passo
dE = -J * (newSpin - oldSpin) * neiSum + D * (newSpin**2 - oldSpin**2)
# Verifica se o passo foi ou não aceito
if np.random.random() < np.exp(-dE / T):
state[index] = newSpin
varE += dE
# Retorna o estado depois de n passos e a variação de energia total
return state, varE
# Função que cria a animação
def aniSpaceEvol(states, D, temp, ene, step=10):
steps = np.shape(states)[0]
time = np.linspace(0,steps,steps)
side = int(np.sqrt(np.shape(states[0]))[0])
layout = """
SS0
SS1
SS2
"""
fig, axes = plt.subplot_mosaic(layout,
figsize=(6, 4),
dpi=72,
layout="constrained")
mag = np.abs( np.mean(states , axis=1))
qua = np.mean(np.abs( states), axis=1)
sysAx = axes['S']
eneAx = axes['0']
magAx = axes['1']
quaAx = axes['2']
sysLine = sysAx.imshow(
states[0].reshape((side,side)),
cmap=plt.get_cmap('plasma',3),
vmin=-1,
vmax=1,
)
cbar = fig.colorbar(sysLine, ax=sysAx)
cbar.set_ticks([-2./3, 0, 2./3])
cbar.set_ticklabels(['Spin -1', 'Spin 0', 'Spin +1'])
eneLine, = eneAx.plot([],[])
magLine, = magAx.plot([],[])
quaLine, = quaAx.plot([],[])
sysAx.set_xlim(-0.5,side-0.5)
eneAx.set_xlim(time[0],time[-1])
magAx.set_xlim(time[0],time[-1])
quaAx.set_xlim(time[0],time[-1])
magAx.sharex(eneAx)
quaAx.sharex(eneAx)
sysAx.set_ylim(-0.5,side-0.5)
eneAx.set_ylim(np.min(ene),np.max(ene))
magAx.set_ylim(0, 1)
quaAx.set_ylim(0, 1)
sysAx.set_title("Distribuição de spin")
eneAx.set_ylabel("Energia")
magAx.set_ylabel("Magnetização")
quaAx.set_ylabel("Quadripolo")
eneAx.tick_params(labelbottom=False)
magAx.tick_params(labelbottom=False)
quaAx.set_xlabel("Passos (MCS)")
sysAx.set_axis_off()
eneAx.grid()
magAx.grid()
quaAx.grid()
def update(frame):
sysLine.set_array(states[frame].reshape((side,side)))
eneLine.set_data(time[:frame], ene[:frame])
magLine.set_data(time[:frame], mag[:frame])
quaLine.set_data(time[:frame], qua[:frame])
sysAx.set_title(f"Temperatura = {temp:1.3f}; D = {D:.3f}")
return [sysLine, eneLine, magLine, quaLine]
anim = animation.FuncAnimation(
fig,
update,
frames=range(0, states.shape[0], step),
interval=50,
blit=True
)
plt.close(fig)
return anim
# Função utilizada para gerar os dados do espaço de fases.
# Utiliza o njit pela complexidade matemática.
# Recebe listas com os valores de D e T que se que simular
# além do tamanho da grade da simulação interna e quantos MCS se quer realizar
@njit
def genPhaseSpace(TStates, DStates, simSize, mcs, J=1.0):
N = simSize ** 2
nT = TStates.size
nD = DStates.size
M = np.zeros((steps, nT, nD))
Q = np.zeros((steps, nT, nD))
lookup = genLookup(simSize)
# utriliza-se o prange para aproveitar dos mútiplos núcleos da CPU
# Itera sobre os valores de T e D que se quer realizar a simulação
for i in prange(nT):
for j in range(nD):
# Cria um estrado inicial para o sitema, igual para todos
state = np.ones(N, dtype=np.int64)
# Calcula-se a Magnetização e o quadrupolo médio do sistema inicial
M[0, i, j] = np.sum(state) / N
Q[0, i, j] = np.sum(np.abs(state)) / N
# Realiza um número recebido de passos de monte carlo no sistema
for k in range(1,mcs):
# A energia após cada passo é descartada, pois não nos interessa aqui
state, _ = timeSteps(
state, lookup,
steps=N,
J=J,
D=DStates[j],
T=TStates[i]
)
# Salva os valors de M e Q no tempo k
M[k, i, j] = np.sum(state) / N
Q[k, i, j] = np.sum(np.abs(state)) / N
return M, Q
def aniPhaseSpace(M, Q, T, D, points = [], step=5):
N = M.shape[0]
extent = [D[0], D[-1], T[0], T[-1]]
frames = range(1, N + 1, step)
MFrames = [np.mean(M[:k], axis=0) for k in frames]
QFrames = [np.mean(Q[:k], axis=0) for k in frames]
MQFrames = [np.abs(m - q) for m, q in zip(MFrames, QFrames)]
fig, axes = plt.subplots(1, 3, figsize=(14, 5), dpi=75, sharey=True)
imshow_kwargs = dict(origin='lower', aspect='auto', extent=extent,
interpolation='nearest', cmap="viridis", vmin = 0, vmax=1)
im0 = axes[0].imshow(MFrames[0], **imshow_kwargs)
im1 = axes[1].imshow(QFrames[0], **imshow_kwargs)
im2 = axes[2].imshow(MQFrames[0], **imshow_kwargs)
axes[0].set(ylabel="Temperatura T")
plt.colorbar(im2, ax=axes[-1], fraction=0.04, pad=0.04)
imTitle = fig.suptitle('', fontsize=12)
for ax, im, title in zip(axes,
[im0, im1, im2],
[
'<|M|> (Magnetização)',
'<Q> (Momento quadripolar)',
'|<M> - <Q>| (Diferença entre pontos)'
]):
for i, point in enumerate(points):
ax.plot(point[0], point[1], 'bo', label='Ponto simulados' if i == 0 else None)
ax.plot(1.965, 0.61, 'ro', label='PTC (teórico)')
ax.set(xlabel='D', title=title)
ax.legend()
def update(idx):
im0.set_array(MFrames[idx])
im1.set_array(QFrames[idx])
im2.set_array(MQFrames[idx])
imTitle.set_text(f'Tempo {list(frames)[idx]}/{N}')
return [im0, im1, im2, imTitle]
anim = animation.FuncAnimation(
fig, update,
frames=len(list(frames)),
interval=80,
blit=True
)
plt.tight_layout()
plt.close(fig)
return anim
gridSize = 20
mcs = 30_000
states = np.zeros((mcs,gridSize**2),dtype=int)
D, T, J = (1.965, 0.61, 1)
lookup = genLookup(gridSize)
states[0], E_i = timeSteps(
states[0].copy(),
lookup,
D = D,
T = T,
steps=50*gridSize**2
)
E = np.zeros(mcs)
side = int(np.sqrt(states[0].size))
E[0] = D * np.sum(states[0]**2)
E[0] += - J * np.sum(states[0] * (states[0][lookup[1]] + states[0][lookup[2]]))
for i in range(1,mcs):
states[i], E_i = timeSteps(
states[i-1].copy(),
lookup,
D = D,
T = T,
steps=gridSize**2
)
E[i] = E[i-1] + E_i
anim = aniSpaceEvol(states, step = 100, D=D, temp = T, ene = E/(gridSize**2))
anim.save('TCP.gif', writer='ffmpeg')
simSize = 20
gridSize = 50
mcs = 500
J = 1.0
TStates = np.linspace(0.1, 5.1, gridSize)
DStates = np.linspace(0.0, 5.0, gridSize)
M,Q = genPhaseSpace(TStates, DStates, simSize, mcs)
anim = aniPhaseSpace(
M, Q,
TStates, DStates,
points = [(1.965,0.91), (1.665,0.61), (2.265,0.61)],
step = 5
)
HTML(anim.to_jshtml())
anim.save('phaseSpace.gif', writer='ffmpeg')