Código:Modelo de Blume-Capel

De Física Computacional
Ir para navegação Ir para pesquisar 
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np

from IPython.display import HTML
from numba import njit, prange

# 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')
Disponível em “http://fiscomp.if.ufrgs.br/index.php?title=Código:Modelo_de_Blume-Capel&oldid=11544”