Exibir código-fonte para Código:Modelo de Blume-Capel

<syntaxhighlight lang="python">
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np

from IPython.display import HTML
from numba import njit, prange
</syntaxhighlight>

<syntaxhighlight lang="python">
# 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
</syntaxhighlight>

<syntaxhighlight lang="python">
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')
</syntaxhighlight>

<syntaxhighlight lang="python">
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')
</syntaxhighlight>