<?xml version="1.0"?>
<feed xmlns="http://www.w3.org/2005/Atom" xml:lang="pt-BR">
	<id>http://fiscomp.if.ufrgs.br/index.php?action=history&amp;feed=atom&amp;title=C%C3%B3digo%3AModelo_de_Blume-Capel</id>
	<title>Código:Modelo de Blume-Capel - Histórico de revisão</title>
	<link rel="self" type="application/atom+xml" href="http://fiscomp.if.ufrgs.br/index.php?action=history&amp;feed=atom&amp;title=C%C3%B3digo%3AModelo_de_Blume-Capel"/>
	<link rel="alternate" type="text/html" href="http://fiscomp.if.ufrgs.br/index.php?title=C%C3%B3digo:Modelo_de_Blume-Capel&amp;action=history"/>
	<updated>2026-06-13T00:26:48Z</updated>
	<subtitle>Histórico de revisões para esta página neste wiki</subtitle>
	<generator>MediaWiki 1.43.0</generator>
	<entry>
		<id>http://fiscomp.if.ufrgs.br/index.php?title=C%C3%B3digo:Modelo_de_Blume-Capel&amp;diff=11544&amp;oldid=prev</id>
		<title>Misalocin: Criou página com &#039;&lt;syntaxhighlight lang=&quot;python&quot;&gt; import matplotlib.pyplot as plt import matplotlib.animation as animation import numpy as np  from IPython.display import HTML from numba import njit, prange &lt;/syntaxhighlight&gt;  &lt;syntaxhighlight lang=&quot;python&quot;&gt; # 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. @...&#039;</title>
		<link rel="alternate" type="text/html" href="http://fiscomp.if.ufrgs.br/index.php?title=C%C3%B3digo:Modelo_de_Blume-Capel&amp;diff=11544&amp;oldid=prev"/>
		<updated>2026-06-01T01:33:32Z</updated>

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