O gás de rede

Neste capítulo vamos mostrar um uso de várias espécies de turtles em NetLogo, assim como a capacidade dos indivíduos de estas mudarem de espécie. Introduziremos também as listas e explicaremos como se usam. Fá-lo-emos usando o modelo do gás de rede para demonstrar algumas características da transição de fase líquido vapor, assim como da tensão superficial dos líquidos.

O gás de rede é um modelo onde os agentes são patches de dois tipos: –ocupados por uma partícula (pretos) e desocupados (brancos), os quais são referenciados por turtles de duas espécies: –ocupados e livres, as quais funcionam como etiquetas dos patches. Tal permite-nos aplicar código NetLogo somente aos patches com uma certa característica, sem estar constantemente a verificar se todos os patches satisfazem essa característica ou não (ask patches with [ variavel = valor ]), o que permite acelerar o programa.

Numa substância real as partículas interagem entre si, quer por pontes de hidrogénio, quer por forças de Van der Waals, no gás de rede modelamos as interacções por proximidade geométrica: os patches pretos ("partículas") interagem, se forem adjacentes na horizontal ou na vertical. Consequentemente cada patch preto pode interagir com quatro vizinhos, no máximo. Atribuímos a cada uma destas interacções a energia -1 (podemos interpretá-la como a energia média de interacção num sistema de unidades da energia conveniente).

O modelo

Com o gás de rede pretendemos modelar um fluido contido num recipiente fechado, com volume constante, imerso num banho que mantém a temperatura do fluido constante.

Nestas condições as partículas do fluido estão em equilíbrio térmico e estão constantemente a transitar entre os vários níveis de energia (-4, -3, -2, -1 e 0), cf. diagrama da direita. No entanto pela física estatística a transição entre os vários níveis de energia é feita com diferentes probabilidades:

  1. e-|ΔE|/KT, quando a partícula salta para um nível de energia mais elevado, onde -|ΔE| é a diferença de energia entre os dois níveis, T é a temperatura e K a constante de Boltzmann (a mesma que aparece na equação dos gases ideais: pV = NKT),
  2. 1 quando a partícula salta para um nível de energia mais baixo.

As regras anteriores, mediante a escolha conveniente de um sistema de unidades para a temperatura (aquele onde K = 1), traduzem-se no seguinte guião para os patches:

  1. A cada iteração um patch preto e um patch branco são escolhidos ao acaso.
  2. Se a transferência do patch preto para a posição do patch branco resultar numa diminuição da energia, a troca é feita e a partícula transita para o novo estado.
  3. Se pelo contrário resultar numa maior energia a troca é feita com probabilidade p = e-|ΔE|/T para o novo estado, onde |ΔE| é a diferença de energia entre os dois estados e T é a temperatura.

Interagindo com o modelo

Com este programa podemos variar a quantidade de partículas (ou densidade, número de partículas por unidade de área) e a temperatura do sistema de forma independente. Consequentemente este programa vai permitir investigar o diagrama de fases densidade-temperatura do gás de rede, o qual possui as mesmas características do diagrama de fases das substâncias simples correspondente às fases líquido e vapor.

Representação esquemática do diagrama de fases. Abaixo da temperatura crítica o estado de equilíbrio é o líquido a densidades elevadas e o gás a baixas densidades. Numa gama de densidades intermédia, delimitada pelas densidades de coexistência do gás e do líquido a essa temperatura (curva de coexistência), o estado de equilíbrio corresponde à coexistência entre os dois estados. Acima da temperatura crítica o estado estável é fluido com uma densidade que varia continuamente.

O leitor é agora convidado a interagir com este modelo através do seguinte applet:

Para correr este modelo no ambiente NetLogo, basta seguir o link gasderede, se tiver instalado o software.

Neste applet podemos estudar como é que o estado final do fluido depende do estado inicial escolhendo um dos três estados (aleatorio, gota ou interface) no botão estado-inicial para diferentes temperaturas e densidades escolhidas nos cursores respectivos e premindo de seguida o botão preparar. Qualquer uma das três escolhas vai criar um número de patches ocupados com a densidade escolhida. Por exemplo para a temperatura de 0.30 e densidade de 50% obtemos respectivamente os seguintes estados iniciais (aleatório, gota e interface):

A evolução no tempo do sistema é mostrada através do gráfico EM.

O leitor tem à sua disposição mais dois botões densidade_automatica e temperatura_automatica os quais desenham respectivamente a energia média de equilíbrio por partícula em função da densidade para a temperatura seleccionada no cursor temperatura e a energia média por partícula em função da temperatura para a densidade seleccionada no cursor densidade. Para limpar os gráficos, basta premir o botão limpar.

Conclusões

O efeito do estado inicial

O estado inicial é irrelevante para a obtenção do estado de equilíbrio: –a baixas temperaturas o fluido condensa, a altas temperaturas o fluido evapora independentemente do estado inicial:

  • aleatorio i.e. estado gasoso,
  • gota ou interface i.e. estado líquido,

conforme atestam os seguintes exemplos:

Baixas temperaturas
Estado de equilíbrio para os estados iniciais aleatório, gota e interface para temperatura = .30 e densidade = 10%:

A baixas temperaturas o fluido está no estado líquido, como podemos ver pelo aparecimento de gotas.
No estado líquido é igualmente visível o aparecimento da tensão superficial, pois as gotas tendem a ter forma circular. Tal é consequência de a baixas temperaturas o sistema tender a minimizar a energia. A forma de o fazer é escolhendo a figura geométrica que tenha o menor perímetro para o número de partículas dado (o círculo), pois dessa forma obtém-se o maior número de interacções possível e consequentemente a menor energia. A forma de o fazer em dimensão três é escolhendo a forma geométrica com a menor superfície para o número de partículas dado (i.e. a esfera).
Note-se que a terceira imagem também só contem uma gota pois o mundo é toroidal (i.e. tem a forma de um doughnut), pelo que a gota da parte do topo é a continuação da gota da parte de baixo que aparece do outro lado do mundo.
Um outro efeito da forma toroidal do mundo é que para densidades mais altas a forma das gotas é rectangular (i.e. uma fatia de toro), como podemos ver para os seguintes parâmetros:
O valor da densidade que separa os dois comportamentos é dado por 2πr = 2×63, onde 63 é a largura / altura do toro neste modelo, o que corresponde à densidade 1 / π ≈ 32%.
Altas temperaturas
Estado de equilíbrio para os estados iniciais aleatório, gota e interface para temperatura = .70 e densidade = 10%:
Independentemente do estado inicial podemos observar que o estado de equilíbrio é o estado gasoso.

O efeito da densidade

Para temperaturas suficientemente altas, por forma a que o fluido esteja no estado gasoso, a mera variação da densidade permite fazer o sistema passar de um estado onde a energia por partícula é próxima de 0 (a baixas densidades) para um estado onde a energia média por partícula é próxima de -2 (densidades elevadas), cf. gráfico da direita para a temperatura = 0.70.

A forma do gráfico não é surpreendente, pois para temperaturas altas e densidades baixas as hipóteses de existirem partículas vizinhas são bastante baixas, por outro lado, para densidades altas, não há espaço para que as partículas estejam no estado "gasoso", pelo que as partículas são confinadas ao estado "condensado" , onde interagem com as vizinhas.

O efeito da temperatura

Os fluidos macroscópicos apresentam transições de fase. A água por exemplo passa do estado líquido ao estado gasoso a 100°C. A essa temperatura a energia interna do sistema fechado composto por água no estado líquido e no estado gasoso aumenta à medida que esta evapora, apesar da temperatura se manter constante. Este fenómeno encontra-se representado esquematicamente na curva a azul do gráfico da esquerda. Esta subida na vertical só é observada para sistemas macroscópicos, i.e. compostos por muitas partículas. Para sistemas compostos por relativamente poucas partículas, a subida da energia interna não é tão acentuada no ponto de ebulição, observando-se apenas uma inflecção no ponto de ebulição, cf. curva a vermelho no gráfico da esquerda. O gás de rede, para cada densidade, também apresenta um ponto de ebulição. Por exemplo para a densidade = 50%, o ponto de ebulição é por volta de temperatura = .55, como pode ser visto no gráfico da direita:

É igualmente visível a existência do ponto crítico quando densidade = 50% e temperatura = .55, como podemos ver no gráfico da direita, pois para esta densidade o sistema ferve à temperatura mais alta.

Neste gráfico estão representadas para temperatura = .55 as curvas cinzenta, vermelha e laranja que correspondem respectivamente à energia média por partícula em função da temperatura para as densidades 50%, 90% e 10%. A recta azul é a recta que corta temperatura = .55.

Para obter este gráfico seleccione densidade = 50% prima temperatura_automatica e deixe desenhar a curva cinzenta, depois seleccione densidade = 90% prima temperatura_automatica, mas não prima limpar e deixe desenhar a curva vermelha, por fim repita o mesmo processo para densidade = 10%.

Breve análise do código

O modelo está dividido nas seguintes secções:

Variáveis globais
tempocontém o número de iterações (tentativas de troca) já decorridas
energiacontém a energia média por partícula do sistema
particulascontém o número de partículas (patches pretos)
patches: ligacoespropriedade de cada patch que contem o número de ligações disponíveis, i.e. o número de patches pretos adjacentes na horizontal ou na vertical
Rotinas principais
prepararlimpa e inicia todas as variáveis, cria o número de partículas correspondente à densidade escolhida e coloca-as no estado inicial escolhido pelo utilizador (aleatório, gota, ou interface).
executaritera as regras de transformação dos patches, desenha o gráfico da evolução da energia média por partícula em função do tempo e mostra a energia média instantânea no monitor.
densidade_automaticaitera as regras de transformação dos patches e desenha o gráfico da energia média de equilíbrio por partícula em função da densidade para a temperatura fixada no cursor temperatura.
temperatura_automaticaitera as regras de transformação dos patches, desenha o gráfico da evolução da energia média de equilíbrio por partícula em função da temperatura para a densidade fixada no cursor densidade.
Rotinas auxiliares
prepara-aleatoriocria o estado inicial aleatorio.
prepara-gotacria o estado inicial gota.
prepara-interfacecria o estado inicial interface.
trocaritera as regras de transformação dos patches e é chamada pela rotina executar.
preparar-graficoinicia o gráfico da energia média.
graficomarca os pontos no gráfico da energia média.
calcula-energia-mediacalcula a energia média de equilíbrio por partícula e é chamada pelas rotinas densidade_automatica e temperatura_automatica.

O código

globals [ tempo energia particulas cor ]
breed [ ocupados ocupado ]
breed [ livres livre ]
patches-own [ ligacoes ]

to preparar
    locals [ patches-ocupados ]
    ct cp
    set-current-plot "EM"
    clear-plot
    set tempo 0
    ask patches [ 
        set pcolor white
        sprout-livres 1 [ ht ]
    ]
    if estado-inicial = "aleatorio" [ prepara-aleatorio ] 
    if estado-inicial = "gota" [ prepara-gota ]
    if estado-inicial = "interface" [ prepara-interface ]
    ask patches [ set ligacoes count neighbors4 with [ pcolor = black ] ]
    set patches-ocupados patches-from ocupados [ patch-here ]
    set energia sum values-from patches-ocupados [ligacoes] / (-2 * particulas)
    display
    preparar-grafico
    grafico
end

to prepara-aleatorio
    ask patches [
        if random-float 100 < densidade [ 
            set pcolor black
            ask livres-here [ set breed ocupados ]
        ]
    ]
    set particulas count ocupados
end

to prepara-gota
    locals [ lado ]
    set particulas (count patches) * densidade / 100
    set lado int ((sqrt particulas - 1) / 2)
    ask patches with [ abs pxcor <= lado and abs pycor <= lado ][
        set pcolor black
        ask livres-here [ set breed ocupados ]
    ]
    ask n-of (particulas - (2 * lado + 1) ^ 2) patches with [ 
        (abs pxcor = lado + 1 and abs pycor <= lado + 1) or 
        (abs pxcor <= lado + 1 and abs pycor = lado + 1)
    ][ 
        set pcolor black
        ask livres-here [ set breed ocupados ]
    ] 
end

to prepara-interface
    locals [ altura ]
    set particulas (count patches) * densidade / 100
    set altura int (particulas / world-width)
    ask patches with [ pycor < altura + min-pycor ][ 
        set pcolor black
        ask livres-here [ set breed ocupados ]
    ]
    ask n-of (particulas - altura * world-width) patches with [
        pycor = altura + min-pycor
    ][ 
        set pcolor black 
        ask livres-here [ set breed ocupados ]
    ]
end

to executar
    locals [ ocupado1 livre1 patches-ocupados ]
    no-display
    set ocupado1 one-of ocupados
    set livre1 one-of livres
    trocar ocupado1 livre1
    every 0.1 [ 
        display
        set patches-ocupados patches-from ocupados [ patch-here ]
        set energia sum values-from patches-ocupados [ligacoes] / (-2 * particulas)
        grafico
    ]
end

to trocar [ ocupado1 livre1 ]
    locals [ diferenca ]
    set diferenca ligacoes-of livre1 - ligacoes-of ocupado1
    if random-float 1.0 < exp (diferenca / temperatura) [ 
        ask livre1 [ 
            ask neighbors4 [ set ligacoes ligacoes + 1 ] 
            ask patch-here [ set pcolor black ]
            set breed ocupados
        ]
        ask ocupado1 [ 
            ask neighbors4 [ set ligacoes ligacoes - 1 ] 
            ask patch-here [ set pcolor white ]
            set breed livres
        ] 
    ]
    set tempo tempo + 1
end

to preparar-grafico
    set-current-plot "EM"
    set-current-plot-pen "eixo"
    auto-plot-off
    plotxy 0 -2
    plotxy 100000000 -2
    auto-plot-on
end


to grafico
    set-current-plot "EM"
    set-current-plot-pen "epp"
    plotxy tempo energia
end

to temperatura_automatica
    set estado-inicial "gota"
    set-current-plot "energia / temperatura"
    set-plot-pen-color blue
    ppu
    plotxy .55 -.5
    ppd
    plotxy .55 -2
    if cor = 0 [ set cor 5 ]
    set-plot-pen-color cor
    set temperatura 0.45
    while [ temperatura < 0.65 ][
        preparar
        calcula-energia-media
        set-current-plot "energia / temperatura"
        if temperatura = .45 [
            ppu
            plotxy temperatura energia
            ppd
        ]
        plotxy temperatura energia
        set temperatura temperatura + .01
    ]
    set cor cor + 10
end

to densidade_automatica
    set estado-inicial "gota"
    set-current-plot "energia / densidade"
    clear-plot
    set densidade 1
    while [ densidade < 100 ][
        preparar
        calcula-energia-media
        set-current-plot "energia / densidade"
        plotxy densidade energia
        set densidade densidade + 1
    ]
end

to calcula-energia-media
    locals [ lista ]
    set lista []
    while [ length lista < 50 ][
        executar
        every 0.1 [ 
            display 
            ifelse length lista < 30 [
                set lista lput energia lista
            ][
                ifelse (energia > sum lista / 30) and (length lista = 30) [
                    set lista but-first lista
                    set lista lput energia lista    
                ][
                    set lista lput energia lista
                ]
            ]
        ]      
    ]
    set lista sublist lista 30 49
    set energia sum lista / 20
end

Vejamos em detalhe as componentes lógicas novas do código:

Rotinas principais
preparar 
to preparar
    locals [ patches-ocupados ]
    ct cp
    set-current-plot "EM"
    clear-plot
    set tempo 0
    ask patches [ 
        set pcolor white
        sprout-livres 1 [ ht ]
    ]
    if estado-inicial = "aleatorio" [ prepara-aleatorio ] 
    if estado-inicial = "gota" [ prepara-gota ]
    if estado-inicial = "interface" [ prepara-interface ]
    ask patches [ set ligacoes count neighbors4 with [ pcolor = black ] ]
    set patches-ocupados patches-from ocupados [ patch-here ]
    set energia sum values-from patches-ocupados [ligacoes] / (-2)
    set energia energia / particulas
    display
    preparar-grafico
    grafico
end

Esta rotina começa por limpar turtles e patches que possam existir de execuções anteriores do programa com os comandos ct e cp e por limpar o gráfico da energia média (EM) com os comandos: 

set-current-plot "EM"
clear-plot

De seguida inicia todos os patches com cor branca (desocupados) e associa a cada um uma turtle da espécie livres que funciona como etiqueta dos patches livres e por essa razão escondemo-la: 

ask patches [ 
    set pcolor white
    sprout-livres 1 [ ht ]
]

usamos esta artimanha, pois o NetLogo não nos permite criar espécies de patches.

Depois cria o estado inicial escolhido pelo utilizador: 

set particulas (count patches) * densidade / 100
if estado-inicial = "aleatorio" [ prepara-aleatorio ] 
if estado-inicial = "gota" [ prepara-gota ]
if estado-inicial = "interface" [ prepara-interface ]

De seguida a cada patch atribui o número de ligacoes ("interacções") possíveis nesse patch, número esse que corresponde ao número de patches ocupados e adjacentes na horizontal ou na vertical (neighbors4 with [ pcolor = black ]): 

ask patches [ set ligacoes count neighbors4 with [ pcolor = black ] ]

Para tornarmos o código mais legível criamos o conjunto dos patches-ocupados com o comando set patches-ocupados patches-from ocupados [ patch-here ], que corresponde aos patches referenciados pelo repórter [ patch-here ] aplicado às turtles da espécie ocupados.

De seguida calcula a energia do sistema, a qual é metade do simétrico da soma do número de pares de patches pretos adjacentes na horizontal ou na vertical através do seguinte comando:

set energia sum values-from patches-ocupados [ligacoes] / (-2).

Vemos aqui o primeiro exemplo do uso de listas, pois o comando values-from patches-ocupados [ligacoes] cria uma lista formada pelos valores da variável ligacoes relativos ao conjunto de agentes patches-ocupados, valores esses que são depois somados com a instrução sum.

Dividimos por 2, pois contámos a "ligações" duas vezes: –quando chegámos ao patch A ligado ao patch B, contámos uma "ligação", "ligação" esta que tornou a ser contada quando chegámos ao patch B (que obviamente está ligado ao patch A). Tomamos para a energia o simétrico do número de "ligações", pois as energias potenciais são convencionalmente representadas por números negativos.

De seguida a energia média por partícula é calculada com o comando set energia energia / particulas e por fim activa-se a actualização da janela gráfica (com a instrução display, ver à frente), que pode ter sido desactivada por uma execução anterior do modelo, e inicia-se o gráfico da energia média (EM) marcando o ponto correspondente à configuração inicial: 

display
preparar-grafico
grafico
to executar 
to executar
    locals [ ocupado1 livre1 patches-ocupados ]
    no-display
    set ocupado1 one-of ocupados
    set livre1 one-of livres
    Trocar ocupado1 livre1
    every 0.1 [ 
        display
        set patches-ocupados patches-from ocupados [ patch-here ]
        set energia sum values-from patches-ocupados [ligacoes] / (-2)
        set energia energia / particulas
        grafico
    ]
end

Esta é a rotina responsável pelo guião dos patches. Os comandos 

set ocupado1 one-of ocupados
set livre1 one-of livres

escolhem ao acaso uma etiqueta (turtle) ocupada e uma livre e os patches executam de seguida as regras definidas no guião (Trocar ocupado1 livre1).

O comando no-display é usado para que a janela gráfica não mostre a actualização do mundo, pois essa actualização consome muitos recursos de processamento. Este comando é usado em conjunto com o comando display o qual actualiza o mundo e está activo a cada 0.1 segundos, em consequência do comando every 0.1 [ comandos ]. Tal contribui significativamente para que o modelo corra a uma velocidade aceitável.

densidade_automatica 
to densidade_automatica
    set estado-inicial "gota"
    set-current-plot "energia / densidade"
    clear-plot
    set densidade 1
    while [ densidade < 100 ][
        preparar
        calcula-energia-media
        set-current-plot "energia / densidade"
        plotxy densidade energia
        set densidade densidade + 1
    ]
end

Esta é a rotina que desenha o gráfico da energia média em função da densidade a temperatura constante (dada pelo cursor com o mesmo nome). No gráfico a densidade varia entre 0 e 100%, em intervalos de 1%.

temperatura_automatica 
to temperatura_automatica
    set estado-inicial "gota"
    set-current-plot "energia / temperatura"
    set-plot-pen-color blue
    ppu
    plotxy .55 -.5
    ppd
    plotxy .55 -2
    if cor = 0 [ set cor 5 ]
    set-plot-pen-color cor
    set temperatura 0.45
    while [ temperatura < 0.65 ][
        preparar
        calcula-energia-media
        set-current-plot "energia / temperatura"
        if temperatura = .45 [
            ppu
            plotxy temperatura energia
            ppd
        ]
        plotxy temperatura energia
        set temperatura temperatura + .01
    ]
    set cor cor + 10
end

Esta é a rotina que desenha o gráfico da energia média em função da temperatura a densidade constante (dada pelo cursor com o mesmo nome). No gráfico a temperatura varia entre 0.45 e 0.65, em intervalos de 0.01. Esta rotina desenha também uma recta vertical sobre o ponto correspondente à temperatura = 0.55.

Rotinas auxiliares
prepara-aleatorio 
to prepara-aleatorio
    ask patches [
        if random-float 100 < densidade [ 
            set pcolor black
            ask livres-here [ set breed ocupados ]
        ]
    ]
    set particulas count ocupados
end

Esta é a rotina que cria o estado inicial aleatorio (estado gasoso). A qual ocupa cada patch com a probabilidade densidade / 100 pintando o patch de preto (set pcolor black) e mudando-lhe a etiqueta para ocupados (ask livres-here [ set breed ocupados ]).

prepara-gota 
to prepara-gota
    locals [ lado ]
    set particulas (count patches) * densidade / 100
    set lado int ((sqrt particulas - 1) / 2)
    ask patches with [ abs pxcor <= lado and abs pycor <= lado ][
        set pcolor black
        ask livres-here [ set breed ocupados ]
    ]
    ask n-of (particulas - (2 * lado + 1) ^ 2) patches with [ 
        (abs pxcor = lado + 1 and abs pycor <= lado + 1) or 
        (abs pxcor <= lado + 1 and abs pycor = lado + 1)
    ][ 
        set pcolor black
        ask livres-here [ set breed ocupados ]
    ] 
end

Esta rotina é semelhante à anterior, o que varia é a disposição das partículas.

prepara-interface 
to prepara-interface
    locals [ altura ]
    set particulas (count patches) * densidade / 100
    set altura int (particulas / world-width)
    ask patches with [ pycor < altura + min-pycor ][
        set pcolor black
        ask livres-here [ set breed ocupados ]
    ]
    ask n-of (particulas - altura * world-width) patches with [
        pycor = altura + min-pycor
    ][ 
        set pcolor black 
        ask livres-here [ set breed ocupados ]
    ]
end

Esta rotina é semelhante à anterior, o que varia é a disposição das partículas.

trocar 
to trocar [ ocupado1 livre1 ]
    locals [ diferenca ]
    set diferenca ligacoes-of livre1 - ligacoes-of ocupado1
    if random-float 1.0 < exp (diferenca / temperatura) [ 
        ask livre1 [ 
            ask neighbors4 [ set ligacoes ligacoes + 1 ] 
            ask patch-here [ set pcolor black ]
            set breed ocupados
        ]
        ask ocupado1 [ 
            ask neighbors4 [ set ligacoes ligacoes - 1 ] 
            ask patch-here [ set pcolor white ]
            set breed livres
        ] 
    ]
    set tempo tempo + 1
end

Esta rotina recebe dois parâmetros uma turtle ocupada e uma livre. Se a energia potencial da partícula no patch livre for maior que a energia potencial da partícula no patch ocupado, a troca é feita com probabilidade exp (diferenca / temperatura), caso contrário a troca é sempre feita, pois nesse caso exp (diferenca / temperatura) ≥ 1.

A troca consiste em actualizar o número de ligacoes resultante da transferência da partícula, em mudar a cor do patch onde a turtle (etiqueta) está, através da instrução ask patch-here, e em mudar a espécie da turtle.

preparar-grafico 
to preparar-grafico
    set-current-plot "EM"
    set-current-plot-pen "eixo"
    auto-plot-off
    plotxy 0 -2
    plotxy 100000000 -2
    auto-plot-on
end

Esta rotina prepara o gráfico de energia média, desenhando uma linha de base no valor -2.0, que é o valor mínimo possível para a energia média. Este valor só é atingido se todos os patches estiverem ocupados (densidade = 100%).

grafico 
to grafico
    set-current-plot "EM"
    set-current-plot-pen "epp"
    plotxy tempo energia
end

Esta rotina simplesmente acrescenta ao gráfico da energia média (EM) o novo ponto (tempo, energia).

calcula-energia-media 
to calcula-energia-media
    locals [ lista ]
    set lista []
    while [ length lista < 50 ][
        executar
        every 0.1 [ 
            display 
            ifelse length lista < 30 [
                set lista lput energia lista
            ][
                ifelse (energia > sum lista / 30) and (length lista = 30) [
                    set lista but-first lista
                    set lista lput energia lista    
                ][
                    set lista lput energia lista
                ]
            ]
        ]      
    ]
    set lista sublist lista 30 49
    set energia sum lista / 20
end

Esta é a rotina que calcula a energia média por partícula. Para o fazer deixa a energia por partícula estabilizar em torno do valor de equilíbrio e de seguida pega em vinte valores consecutivos e calcula a média.

A parte aparentemente complicada é o "estabilizar em torno do ponto de equilíbrio", mas tal até que acaba por ser fácil, pois quando se começa com o estado-inicial gota a energia vai crescendo com o tempo até estabilizar em torno do ponto de equilíbrio. Precisamos apenas de um critério para saber que a energia já não está a crescer mais.

As médias móveis, bastante usadas na análise técnica de acções, fornecem esse critério. Para calcularmos uma média móvel pegamos nos últimos n valores da energia e dividimos por n. Neste modelo escolhemos n = 30 para suavizar as flutuações na energia e dizemos que esta já não está a crescer mais se o último valor for inferior à média móvel.

Esta rotina faz um extenso uso de listas. Começa por iniciar a variável lista com uma lista vazia (set lista []), de seguida aumenta, em cada iteração, o tamanho da lista enquanto este for menor que trinta (length lista < 30) adicionando a última energia ao fim da lista com o seguinte comando: 

set lista lput energia lista

Quando a lista tiver trinta elementos, remove a energia mais antiga da lista e adiciona a mais recente usando os seguintes comandos: 

set lista but-first lista
set lista lput energia lista

até que a energia por partícula estabilize. Quando tal acontecer junta à lista as vinte energias seguintes e constrói uma lista com essas últimas vinte energias (set lista sublist lista 30 49) e depois calcula a média.