O pacote LAMMPS como uma ferramenta para o ensino de transição de fases

The LAMMPS package as a tool for teaching phase transition

Resumos

A interpretação e descrição dos fenômenos naturais, intrínseca à Ciência, demanda abstração conceitual, utilizando observações e construções teóricas. Conceitos como átomos e moléculas, junto a equações matemáticas, foram desenvolvidos para descrever interações e evolução temporal das propriedades naturais. A análise criteriosa desses modelos, considerando parâmetros termodinâmicos, viabiliza a investigação de fenômenos complexos. Ferramentas computacionais, como a Dinâmica Molecular, permitem modelar sistemas atômico-moleculares, auxiliando na compreensão das transformações da matéria. O potencial de Lennard Jones, um modelo isotrópico proposto nos anos 1920, é eficaz na descrição termodinâmica de fluidos simples. Neste trabalho, revisaremos este modelo para construir o diagrama de fases da temperatura versus densidade ($T\times \rho$) usando simulações de Dinâmica Molecular com o LAMMPS, disponível gratuitamente online. Obteremos a temperatura e densidade críticas, assim como a região de coexistência das fases líquida e vapor. Isso proporcionará aos professores de física, em todos os níveis educacionais, uma abordagem computacional simples para ensinar sobre transições de fase.

Palavras-chave:
Termodinâmica; Diagrama de Fases; Dinâmica Molecular; Lennard-Jones

The interpretation and description of natural phenomena, intrinsic to Science, require conceptual abstraction, using observations and theoretical constructs. Concepts like atoms and molecules, along with mathematical equations, have been developed to describe interactions and the temporal evolution of natural properties. Careful analysis of these models, considering thermodynamic parameters, enables the investigation of complex phenomena. Computational tools, such as Molecular Dynamics, allow for the modeling of atomic-molecular systems, aiding in the understanding of matter transformations. The Lennard-Jones potential, an isotropic model proposed in the 1920s, is effective in the thermodynamic description of simple fluids. In this work, we will review this model to construct the temperature and density phase diagram ($T\times \rho$) using Molecular Dynamics simulations with LAMMPS, available freely online. We will obtain critical temperature and density, as well as the coexistence region of the liquid and vapor phases. This will provide physics teachers at all educational levels with a simple computational approach to teach about phase transitions.

Keywords:
Thermodynamics; Phase Diagrams; Molecular Dynamics; Lennard-Jones

1. Introdução

A Termodinâmica, como um ramo fundamental da Física, investiga as relações entre calor, trabalho e energia nos mais variados sistemas compostos de um número muito grande de constituintes. Além disso, busca entender como a energia é transferida e transformada entre diversas formas durante os processos termodinâmicos [¹, ²]. Tais sistemas são compostos por muitos corpos (da ordem do Número de Avogrado, $\sim {10}^{23}$ partículas) e tem sido objeto de estudo da Termodinâmica desde o seu surgimento em meados do século XVIII, onde o interesse pelo papel do calor nos processos físicos se tornou primordial, inclusive, ao desenvolvimento tecnológico e econômico [³, ⁴]. Desde então, sistemas completamente diferentes (átomos, moléculas em fluidos, sólidos ou plasmas; gases quânticos em semicondutores ou metais; anãs brancas, sistemas de partículas de leptons, quarks e bósons do universo primordial) tem se mostrado fieis a leis físicas bastante gerais e de natureza fenomenológica envolvendo algumas poucas variáveis macroscópicas [⁵]. Enquanto teoria relacionando as propriedades macroscópicas de um sistema e que analisa como a energia é transformada e como a entropia evolui, a Termodinâmica reconheceu os aspectos aleatório e probabilístico dos processos naturais, afirmando-se como um fundamento para a energia e a entropia que governam o mundo [⁶].

As condições termodinâmicas, descritas através de suas variáveis de estado, nas quais uma determinada substância existirá em suas diferentes fases (como gás, líquida e sólida) é sumarizada em um diagrama de fases [⁷]. Tais diagramas permitem representar graficamente as possíveis transições de fase, como, por exemplo, a fusão, a evaporação e a condensação, além de fornecer uma compreensão valiosa sobre o comportamento de materiais em diferentes regimes termodinâmicos, contribuindo para a compreensão de fenômenos complexos que cercam o nosso dia a dia, e também aplicações práticas em diversas áreas da ciência e engenharia [⁸]. Dada a importância da água para a nossa sobrevivência, seu diagrama de fases é um exemplo amplamente utilizado para fins pedagógicos em todos os níveis de ensino [⁹]. Neste, ilustram-se os distintos estados de agregação (sólido, líquido e gasoso), através dos quais a água pode existir de acordo com a temperatura e pressão. Através do diagrama, pode-se estudar as curvas de aquecimento e resfriamento, bem como a relação entre pressão e ponto de ebulição/ condensação [¹⁰]. Tais conceitos estabelecem uma base importante para compreender os comportamentos térmicos da matéria [¹¹].

É possível aprimorar significativamente a obtenção de diagramas de fases por meio da aplicação de simulações computacionais. Nelas, a evolução dos sistemas físicos pode ser analisada tanto de forma determinística (através da técnica de Dinâmica Molecular [¹²]) quanto probabilística (via técnica de Monte Carlo [¹³], por exemplo), gerando médias termodinâmicas que representam cada ponto no gráfico dos diagramas [¹⁴]. Desta forma, o computador se torna um poderoso laboratório, tanto de pesquisa quanto de ensino, permitindo a elucidação de conceitos físicos a todos os níveis de instrução, particularmente aqueles relacionados às transições de fase, que desempenham um papel significativo em muitos aspectos da vida cotidiana, desde a culinária até a climatização de ambientes. Seja no ensino médio ou no ensino superior, a exploração do diagrama de fase une a teoria termodinâmica à realidade prática, enriquecendo a compreensão dos alunos sobre as complexidades do mundo físico que nos rodeia [¹⁵, ¹⁶].

Simulações computacionais vêm sendo utilizadas como ferramentas didáticas no ensino de ciências há bastante tempo [¹⁷, ¹⁸, ¹⁹, ²⁰, ²¹]. Com a crescente integração tecnológica na educação, especialmente após a pandemia de Covid-19, a adoção de simulações computacionais como uma ferramenta para aprimorar a compreensão de conceitos científicos é agora amplamente aceita por educadores em todo o mundo. Como resultado, uma diversidade de aplicações computacionais está sendo desenvolvida e implementada no ensino de Física [²², ²³, ²⁴]. Os efeitos da integração de simulações computacionais no ensino de Física também foram analisados por uma extensa revisão realizada por Banda e Nzabahimana (2021), que ratificaram a recomendação da integração das simulações interativas para aprimorar a compreensão conceitual da física dos alunos [²⁵]. No contexto das transições de fase, softwares computacionais tendem a otimizar o aprendizado (e a discussão) de conceitos mais abstratos a partir, por exemplo, do comportamento das funções termodinâmicas na região de criticalidade, evidenciando os expoentes críticos a partir da generalização das transições de fase já conhecidas.

Neste artigo, buscamos estimular professores de Física (e de Ciências em geral) a utilizarem outro tipo de simulação computacional – aquela geralmente voltada à pesquisa científica, e disponível em softwares gratuitos que nos permitem caracterizar termodinamicamente as mais variadas substâncias. Assim sendo, apresentaremos uma técnica de construção de diagrama de fases a partir de simulações de Dinâmica Molecular para um modelo de fluido simples, composto por apenas um tipo de partícula e que interage via potencial de Lennard-Jones [²⁶]. Utilizaremos o pacote de simulações gratuito LAMMPS, acrônimo do termo em inglês Large-scale Atomic/Molecular Massively Parallel Simulator – um código clássico de dinâmica molecular amplamente utilizado no contexto de pesquisas. Esse trabalho está organizado da seguinte forma: na Seção ², iremos apresentar um resumo sobre Transições de Fase e Curva de Coexistência; na Seção ³, discutiremos a técnica de Dinâmica Molecular, o modelo de Lennard-Jones, o método computacional utilizado para construção do diagrama de fases, além de técnicas usadas para análise dos dados gerados; na Seção ⁴, iremos apresentar o pacote de simulação LAMMPS e o código utilizado na execução desse trabalho; na Seção ⁵, apresentaremos os resultados obtidos e sua discussão; e, finalmente, na Seção ⁶, iremos apresentar nossas considerações finais.

2. Transições de Fase e Curva de Coexistência

As transições de fase são fenômenos físicos que ocorrem quando uma substância passa por uma alteração em seu estado de agregação, como da fase líquida para a fase gasosa (por exemplo, ao fervermos a água), ou da fase sólida para a fase líquida (quando o gelo derrete). Tais mudanças de fase são governadas por propriedades termodinâmicas e ocorrem em temperaturas específicas para diferentes valores de pressão, pelos quais o sistema está sujeito. Para entender esse conceito, é útil explorar a transição e coexistência líquido-gás da água, um fenômeno que encontramos com frequência em nosso cotidiano.

Evidenciando a transição de fase líquido-gás, a qual ocorre quando a água muda de seu estado líquido para o estado gasoso, vê-se que tal fenômeno é influenciado pela temperatura e pressão, de acordo com os seguintes passos:

1. 1.

   Aquecimento da água: inicialmente, a água está no estado líquido. Quando aquecida, suas moléculas ganham energia cinética, o que as faz vibrar mais rapidamente, separando-se umas das outras.

2. 2.

   Atingindo o Ponto de Ebulição: para o caso de sistemas fechados, como por exemplo a água dentro de uma panela com uma tampa bem vedada, à medida que a temperatura aumenta, a energia cinética das moléculas da água também aumenta até atingir uma temperatura específica, o equilíbrio é atingido quando a pressão da fase líquida se iguala a pressão de vapor, conhecida como ponto de ebulição. Como estamos injetando calor ao sistema continuamente, essa energia adicionada é utilizada para superar as forças intermoleculares de atração que mantêm as moléculas de água juntas como um líquido, aumentando a distância entre as mesmas. Portanto, a energia adicionada durante a fervura contribui tanto para a energia cinética das moléculas (à medida que se movem mais rapidamente) quanto para a energia potencial do sistema (à medida que as forças intermoleculares são superadas), desta forma ocorre a transição da substância de um estado líquido para um estado gasoso. Especificamente, para a água ao nível do mar e a uma pressão atmosférica normal, o ponto de ebulição é de 100 graus Celsius (ou 373K). Neste ponto, é importante relembrar que para sistemas abertos, como é o caso de uma poça de água ou roupa molhada secando sob o sol, ocorre o processo de evaporação, abaixo desse ponto de ebulição. Na evaporação apenas as moléculas na superfície do líquido escapam para o ar e, mesmo a temperaturas muito baixas como em dias frios de inverno, algumas moléculas na superfície da água ainda têm energia suficiente para evaporar, tal que taxa de evaporação também é influenciada por fatores como a umidade relativa do ar, a superfície de contato, o movimento do ar sobre a poça (ventilação), e a pressão atmosférica. Em condições de baixa umidade, alta ventilação e pressão atmosférica normal, a evaporação ocorre mais rapidamente.

Durante a transição de fase, há um equilíbrio entre a fase líquida e a fase gasosa. Isso significa que, mesmo que parte da água tenha se transformado em vapor, ainda haverá uma parte remanescente no estado líquido. Tal coexistência é fundamental, sendo representada pela curva de coexistência – uma representação gráfica das condições de temperatura e pressão (ou pressão e densidade, ou temperatura e densidade, etc) em que duas fases diferentes de uma substância podem coexistir em equilíbrio. No caso da água, a curva de coexistência entre a fase líquida e a fase gasosa é bem conhecida, como mostrado na Figura 1. Os pontos que formam essa curva, e que separam a região líquida da região gasosa, representam as condições em que a água está passando por uma transição de fase. Se você estiver abaixo dessa linha no diagrama de fases a água estará no estado gasoso; se estiver acima, estará no estado líquido. Isso permite entender como a água pode coexistir nas duas fases, dependendo da temperatura e pressão, sobre os pontos da linha de transição (ou de coexistência). Essa linha termina no Ponto Crítico (PC): acima destes valores de temperatura e pressão, as duas fases se tornam indistinguíveis.

Figura 1
Diagrama pressão vs temperatura da água. As linhas tracejadas representam as curvas de coexistência. A esfera inferior indica o Ponto Triplo (PT), onde as fases gasosa, líquida e sólida coexistem, e esfera superior indica a ocorrência do Ponto Crítico (PC).

Voltando à linha de coexistência em um diagrama Pressão P e Temperatura T (conforme ilustrado na Figura 1), temos que, para cada ponto contido nessa linha, haverá dois valores de volume, ou densidade, para um mesmo valor de temperatura e pressão. Outra forma útil de representar o sistema é, então, através de isotermas em um diagrama Pressão P e volume V, como ilustrado na Figura 2(a). Neste diagrama, ao invés de haver uma curva de coexistência, teremos uma região de coexistência delimitada pela curva em cor vermelha, como se pode ver na Figura 2(a). As curvas, que se iniciam a altas pressões e baixos volumes correspondendo à fase líquida (ver seta preta à esquerda), atravessam a região de coexistência e, terminam a baixas pressões e altos volumes correspondendo à fase de gás (ver seta preta à direita). Essas curvas representam isotermas, ou seja, os pontos contidos nelas são de mesma temperatura. Além disso, a região de coexistência é delimitada acima por um ponto de inflexão (ver seta vermelha), ou o ponto crítico. Como se pode ver através da Figura 2(a), na região de coexistência (ver área em vermelho), o sistema apresenta patamares no espaço p-V¹. Tais patamares, de pressão e temperatura constantes, indicam os dois valores de volume que o sistema ocupa: um menor, correspondente à fase líquida, e outro maior, correspondente à fase gasosa. Assim, nessa região, teremos ambas as fases líquido e gás coexistindo para uma mesma dada temperatura e pressão.

Figura 2
(a) Diagrama pressão P – volume V com todas as isotermas até o ponto crítico – representado pela seta em vermelho. (b) Diagrama Temperatura T – densidade $\rho$ para o Neônio, Monóxido de Carbono, Nitrogênio e Oxigênio. Figura adaptada da Ref. [²⁷].

O ponto crítico, como vimos, caracteriza-se por uma inflexão no diagrama P–V, assinalando o término da região de coexistência. É interessante examinar a maneira como as densidades do líquido e do vapor se aproximam uma da outra, tornando-se idênticas no ponto crítico. Desta forma poderemos introduzir a noção de expoentes críticos, e para isso iremos assumir que a diferença entre as densidades líquido e vapor se anulem no ponto crítico de acordo com a lei de potência:

(1) $${\rho }_l-{\rho }_v\approx {(T_c-T)}^{\beta },$$

onde ${\rho }_l$, ${\rho }_v$, T_c e T são a densidade da fase líquida, densidade da fase vapor, temperatura crítica e temperatura, respectivamente [¹, ²]. Ao supor essa aproximação, podemos reescrever a equação acima multiplicando ${(T_c-T)}^{\beta }$ por uma constante de proporcionalidade:

(2) $${\rho }_l-{\rho }_v=A{(T_c-T)}^{\beta }.$$

A análise do comportamento ao redor do ponto crítico, onde as densidades do líquido e do vapor se aproximam, tornando-se idênticas, foi apresentada por Guggenheim, em 1945, quando publicou originalmente o diagrama de onde retiramos os dados apresentados na Figura 2(b) [²⁷]. Tal resultado é um diagrama Temperatura T e densidade $\rho$ e, assim como o diagrama P–V, ele apresenta uma região de coexistência das fases líquidas e vapor quando analisamos a variação da temperatura com a densidade ao longo de curvas de pressão constante (isóbaras). Ajustando a equação 2 a estes resultados, e utilizando o expoente crítico experimental $\beta =0.35$ [²], pode-se estimar a temperatura crítica T_c. Por outro lado, é possível também obter a densidade crítica; para isso, usa-se a regra de média Caillet-Mathias [²], ou a lei do diâmetro retilíneo, segundo a qual, ao redor de um ponto crítico, a média aritmética entre ambas as densidades, líquido e vapor, é linear com a temperatura crítica. Desta forma, podemos escrever,

(3) $$\frac{{\rho }_l+{\rho }_v}{2}={\rho }_c+B(T_c-T),$$

onde ${\rho }_c$ e B são a densidade crítica e um parâmetro ajustável, respectivamente. Ajustando os resultados ilustrados na Figura 2 (b) à equação acima, pode-se, então, estimar o valor da densidade crítica. Vale ressaltar que essa regra deixa de ser válida muito próxima do ponto crítico e outros métodos devem ser utilizados [²].

3. Modelos e Métodos Computacionais

3.1 A Dinâmica Molecular Clássica

A proposta deste trabalho é utilizar simulações de Dinâmica Molecular (DM) Clássica para o ensino de Termodinâmica. A DM permite estudar problemas que envolvem sistemas com muitas partículas a partir do conhecimento das energias potenciais de interação entre elas, sendo possível analisar tanto as propriedades de transporte quanto as de equilíbrio do sistema de interesse [²⁸].

De forma simplificada, em uma simulação de DM deve-se resolver as equações clássicas de movimento de Newton para cada partícula do sistema,

(4) $$m_i{\ddot{\mathbf{r}}}_i={\mathbf{f}}_i,$$

onde m_i é a massa da i-ésima partícula e

(5) $${\mathbf{f}}_i=-{\nabla }_{{\mathbf{r}}_i}U,$$

é a força resultante sobre a partícula i dada pela sua interação com as outras partículas que compõem o fluido. $U({\mathbf{r}}_i)$ é a energia potencial que representa as interações do sistema de N partículas e que define o Campo de Forças necessário à obtenção da força resultante. Na ampla maioria das situações, visando-se reduzir o custo computacional das simulações, o Campo de Forças será um potencial efetivo entre pares de partículas que levará em conta todas as demais interações (entre tripletos, entre quadrupletos, etc), desde que consiga reproduzir os dados experimentais [²⁹]. Neste sentido, o potencial de Lennard Jones é capaz de descrever o comportamento de diversos sistemas [²⁸]. Contudo, sistemas iônicos possuem termos coulombianos no seu Campo de Forças (por conta da interação eletrostática) e a modelagem de fluidos poliatômicos exige a definição de potenciais de ligação entre as partículas que compõem cada molécula. Portanto, simulações de DM implicam em resolver a equação diferencial (4) para cada partícula que compõe o sistema, resultando em 3N equações para um sistema tridimensional.

O algoritmo padrão de um programa de Dinâmica Molecular (DM) é representado no diagrama da Figura 3. Após a leitura da configuração inicial, o próximo passo é calcular a força resultante atuante sobre cada partícula. A natureza dessa força depende do Campo de Forças adotado e, no nosso caso, utilizamos o potencial de Lennard-Jones. Para sistemas mais complexos, modelos mais elaborados são necessários. Uma vez conhecida a força resultante, podemos determinar a aceleração de acordo com a equação de Movimento (4). A partir deste ponto, iniciamos a integração das equações de movimento, avançando o tempo até alcançar o tempo máximo de simulação, t_max. A cada incremento infinitesimal de tempo $\delta t$, conhecido como o passo de simulação, as equações de movimento são resolvidas para calcular as novas posições e velocidades de cada partícula nos instantes t e $t+\delta t$. Existem várias técnicas numéricas para realizar essa integração, sendo o algoritmo de velocity – Verlet um dos mais amplamente empregados [²⁹]. Para fins de simplicidade, não entraremos em detalhes sobre os aspectos específicos da integração aqui. Após essa integração das equações de movimento, um programa de DM normalmente calcula várias propriedades de interesse, como a energia do sistema, pressão, temperatura, entalpia e uma gama de outras propriedades. Além disso, é possível armazenar as trajetórias das partículas para análises posteriores e pós-processamento: através de softwares específicos (como ovito [³⁰]), várias análises podem ser realizadas sobre as propriedades do sistema, a saber: possíveis transições de fase, ordem translacional e rotacional, comportamento dinâmico.

Figura 3
Fluxograma esquemático de algoritmo padrão de um programa de Dinâmica Molecular (MD) [²⁹].

3.2 O Potencial de Lennard-Jones

O principal potencial² de pares utilizado para representar a interação entre as partículas de um sistema em DM é o potencial de Lennard-Jones [²⁶],

(6) $$U(r)=4ϵ\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right],$$

onde $ϵ$ e $\sigma$ são a profundidade do poço do potencial e o diâmetro da partícula, respectivamente. Esse potencial simples é capaz de reproduzir a força atrativa que surge devido às forças de London/forças de van der Waals que surgem pelos dipolos induzidos momentaneamente quando dois átomos começam a se aproximar e, a força repulsiva que surge a pequenas distâncias, menores que $2^{1/6}\sigma$ – como indicado na Figura 4. Essa repulsão ocorre devido à sobreposição das nuvens eletrônicas dos átomos e, é relacionada à força de troca do princípio de exclusão de Pauli [³¹].

Figura 4
Potencial e Força de Lennard-Jones em função da distância r. Setas indicam a distância $\sigma$, a distância de equilíbrio r_eq e a distância r_Fmax, onde a força atrativa, em módulo, é máxima.

Para facilitar a execução dos passos de DM, é conveniente estabelecer um raio de corte $r_c=2.5\sigma$ para o potencial, assim, se as moléculas estiverem a uma distância de separação igual ou maior que r_c, a interação entre ambas as moléculas não será levada em consideração durante a evolução temporal do sistema, ou seja, $U(r\ge r_c)=0$. Este potencial de interação está ilustrado pela curva de cor preta na Figura 4 juntamente com a força de interação (curva pontilhada vermelha) dada pelo gradiente negativo do potencial $F(r)=-{\nabla }_{r}U$. A primeira observação que pode ser feita acerca desse potencial é que $\sigma$ é a distância na qual $U(r)=0$ e $ϵ$ é a profundidade do potencial de interação, ambos estão indicados no gráfico. Por outro lado, se fizermos $F(r)=0$, podemos encontrar a distância de equilíbrio r_eq (indicada no gráfico), de forma que se duas partículas encontram-se a uma distância de separação menor que r_eq, elas sentirão uma repulsão, enquanto se estiverem separadas por uma distância maior, estarão sujeitas a uma atração. Outra indicação no gráfico é a distância onde a força é máxima r_Fmax e, nota-se que para distâncias maiores que esta, a força decresce continuamente até alcançar a distância de corte r_c.

Não somente, é interessante descrever o sistemas em unidades reduzidas (ou adimensionais). Como podemos ver pela Figura 2(b), a curva de coexistência líquido-gás de diversas substâncias caem sobre uma mesma curva. Em simulações de DM, é conveniente a utilização destas unidades, já que estas permitem a descrição qualitativa de uma ampla faixa de pressões e temperaturas. Através delas, pode-se escrever propriedades, tais como: a densidade, a temperatura, a pressão, assim como outras [²⁸]. Esse procedimento facilita, por exemplo, a execução de simulações de DM, já que lidar com constantes muito pequenas, como a constante de Boltzmann k_B, ao mesmo tempo com unidades muito grandes, pode levar a erros matemáticos. Além desse fato, temos uma motivação maior para tratar sistemas em unidades reduzidas, que é o fato de permitir múltiplas combinações de densidade $\rho$, temperatura T e parâmetros de Lennard-Jones, como a profundidade do potencial $ϵ$ e o diâmetro da partícula $\sigma$, que correspondem ao mesmo estado termodinâmico [²⁸]. Desta forma, definimos as unidades reduzidas como

(7) $${\rho }^{*}\equiv \rho {\sigma }^3,T^{*}\equiv T\frac{k_B}{ϵ}e{p}^{*}\equiv \frac{p{\sigma }^3}{ϵ}.$$

Por exemplo, para uma dada simulação com densidade ${\rho }^{*}=1.0$ e temperatura $T^{*}=1.0$ e $p^{*}=1.0$, se convertermos essas unidades reduzidas em unidades reais (usando as Equações 7), como mostrado na Tabela 1 para o caso do elemento Argônio, temos um sistema a uma densidade $\rho =1680$ kg/m³, uma temperatura $T=119.8$ K e pressão $p=41.9$ MPa. Dessa forma, confere-se aos resultados obtidos nas simulações de DM uma dada universalidade. Por simplicidade, iremos omitir em nossas notações o símbolo $*$ (indicativo de unidades reduzidas), sem qualquer perda de generalidade.

Tabela 1
Conversão das unidades reduzidas para unidades reais do elemento Argônio fazendo uso da equação 7 e os parâmetros $ϵ/k_B=119.8$ K, $\sigma =3.405\times {10}^{-10}$ m, respectivamente [²⁸].

Propriedade	Unidades Reduzidas		Unidades Reais
Temperatura	$T^{*}=1.0$	$\leftrightarrow$	$T=119.8$ K
Densidade	${\rho }^{*}=1.0$	$\leftrightarrow$	$\rho =1680$ kg/m3
Pressão	$p^{*}=1.0$	$\leftrightarrow$	$p=41.9$ MPa

3.3 Obtendo as propriedades do sistema

Com o objetivo de construir o diagrama de fase $T\times \rho$ (Temperatura T em função da densidade $\rho$) do sistema, iremos utilizar a metodologia descrita em Watanabe et al. [³²], Blas et al. [³³] e Silmore et al. [³⁴]. O primeiro estudou o diagrama de fases de um sistema de partículas monoatômicas, enquanto os outros dois, polímeros de diferentes comprimentos. Portanto, a metodologia consiste em introduzir em uma caixa de simulação de tamanho fixo L, um número N de partículas, de forma a ter uma densidade $\rho$ líquida alta a uma dada temperatura T fixa não tão grande de maneira a ultrapassar a temperatura crítica T_c do sistema (aquela temperatura a partir da qual as fases líquida e gasosa serão indistinguíveis). Em seguida, aumenta-se um dos lados da caixa cúbica, normalmente na direção z, centralizando a fase líquida entre $-L_z/2e{L}_z/2$, de forma a se ter uma caixa paralelepípeda, como ilustrado na Figura 5. Além disso, iremos utilizar ao longo de nossas simulações condições de contorno periódicas, ou seja, quando uma partícula atravessa um dos lados da caixa, ela é substituída instantaneamente no lado oposto por outra com a mesma velocidade a fins de se manter a continuidade do sistema.

Figura 5
Representação esquemática do sistema simulado do plano yz, com o eixo x perpendicular à folha.

Com o sistema pronto, realiza-se equilibrações a T fixa, permitindo que algumas partículas se libertem da fase líquida e comecem a povoar os espaços vazios da caixa. Como a região densa está centralizada na origem da caixa, ela ocupa os espaços entre $-10\sigma$ e $10\sigma$ na direção x, entre $-10\sigma$ e $10\sigma$ na direção y, e entre $-10\sigma$ e $10\sigma$ na direção z. A partir do momento que partículas passam a ocupar posições maiores que $10\sigma$ e menores que $-10\sigma$ na direção z, é possível observar a coexistência das fases vapor e líquido. Para estimar as densidades líquido ${\rho }_l$ e vapor ${\rho }_v$, calcula-se o perfil de densidade ao longo do eixo z da caixa; assim, se definem as regiões correspondentes às fases de baixa e alta densidade. Em seguida, faz-se um ajuste com a seguinte função sigmoide [³³],

(8) $$\rho (z)=\frac{1}{2}({\rho }_l+{\rho }_v)-\frac{1}{2}({\rho }_l-{\rho }_v)\tanh \left(\frac{z-z_0}{d}\right),$$

onde ${\rho }_l$, ${\rho }_v$, z₀ e d são parâmetros ajustáveis. Os dois últimos parâmetros correspondem à posição da interface vapor-líquido e à espessura da mesma, respectivamente. Vale ressaltar que a espessura interfacial d é uma terminologia cuja referência se faz à distância entre dois pontos específicos, normalmente associados a uma interface ou à fronteira entre duas fases. Neste caso, é uma medida com dimensão de distância $\sigma$. Para melhor ilustrar, pode-se utilizar o exemplo da interface entre água e óleo. Se temos um tanque com tal mistura, ao se medir a espessura da interface onde ambos se encontram, estaremos medindo a distância entre os pontos onde a concentração de óleo é de $10\mathrm{\%}$ e $90\mathrm{\%}$ do total. Tal medida nos ajuda a entender o quão misturadas estão as substâncias. Se a espessura da fronteira é grande, isso significa que o óleo e a água não estão bem misturados, entretanto se a espessura é pequena, isso indica uma boa mistura. Então, em resumo, quando falamos sobre a espessura interfacial, estamos simplesmente medindo o quão espessa é a camada onde duas substâncias se encontram e, para entender o quão bem misturadas elas estão. Voltando à equação 8, após definidas as densidades da fase líquida e da fase de vapor, podemos utilizar a diferença entre elas, usando a equação 2, e ajustar nossos resultados para obter a temperatura crítica T_c do sistema. Além disso, através da lei dos diâmetros retilíneos dado pela equação 3, podemos estimar o valor da densidade crítica ${\rho }_c$ através de mais um ajuste de curva. Seguindo todos esses passos, obtendo os valores de ${\rho }_l$, ${\rho }_v$, T_c e ${\rho }_c$, pode-se construir o diagrama de fases $T\times \rho$.

Neste trabalho, simulações foram realizadas a uma T fixa com número fixo de partículas N = 6080, e o controle da temperatura foi realizada através do termostato de Langevin [²⁸], o qual adiciona um termo aleatório ao movimento das partículas. Este termostato simula a influência de um banho térmico por meio da adição de duas forças na equação 5: uma força de viscosidade proporcional à velocidade das partículas e uma força aleatória de natureza gaussiana, representando as colisões com o banho térmico. A inclusão das forças viscosa e aleatória, como no caso da dinâmica de Langevin, é necessária para modelar um sistema à temperatura constante devido à natureza aleatória do processo de equilíbrio térmico. Estas forças simulam o efeito do ambiente térmico sobre as partículas do sistema, garantindo que as mesmas mantenham uma distribuição de velocidades típica de um ensemble canônico – a distribuição de Maxwell-Boltzmann [³⁵]. No entanto, existem outros termostatos que podem ser empregados para manter a temperatura constante em simulações de DM, são eles:

• •

  Termostato de Berendsen [³⁶]: aplica uma escala de velocidade a todas as partículas do sistema com base na diferença entre a temperatura atual do sistema e a temperatura desejada. No entanto, ele não reproduz fielmente a distribuição de Maxwell-Boltzmann.

• •

  Termostato de Nose-Hoover [³⁷, ³⁸]: introduz variáveis dinâmicas adicionais que evoluem ao longo do tempo e são acopladas às partículas do sistema. Ele ajusta as velocidades das partículas para manter a temperatura constante, permitindo uma distribuição de velocidades mais realista.

Esses são apenas alguns exemplos de termostatos que podem ser utilizados para manter a temperatura constante em simulações de DM. A escolha do termostato adequado depende das propriedades específicas do sistema que está sendo simulado e dos objetivos da simulação.

Para a produção do diagrama de fases, diversas simulações a diferentes temperaturas (no intervalo entre 0.5 e 1.2) foram realizadas, onde se equilibrou o sistema com $1\times {10}^6$ passos temporais, permitindo que algumas partículas se libertassem da fase líquida e começassem a povoar os espaços vazios da caixa paralelepípeda. Visando assegurar que o sistema estava em equilíbrio, o comportamento da Energia Cinética e Energia Potencial do sistema foram analisadas (Figura 6). Nota-se que ao fim da equilibração, ambas as energias não variavam demasiadamente, indicando que o sistema se encontra em equilíbrio termodinâmico. Após esses passos, realizamos mais $1\times {10}^6$ passos de simulação para a produção dos resultados. Os códigos para realizações das simulações e análise dos dados em DM podem ser encontrados no repositório de códigos https://github.com/thiagopuccinelli/slab-rbef. Na próxima seção, apresentaremos o pacote de simulações usado neste trabalho.

Figura 6
Energia potencial e cinética do sistema em função dos passos de simulação para o sistema, cuja temperatura é $T=0.5$.

4. Apresentação do Pacote de Simulação LAMMPS e script de entrada

O LAMMPS é um código clássico de DM [³⁹] que modela sistemas atômicos, poliméricos, biológicos, de estado sólido (metais, cerâmicas, óxidos), granulares, coarse-grained ou macroscópicos, utilizando uma variedade de potenciais interatômicos (campos de força) e condições de contorno. Ele pode modelar sistemas bidimensionais ou tridimensionais com tamanhos que vão desde apenas algumas partículas até bilhões.

No sentido mais geral, o LAMMPS integra as equações de movimento de Newton para um conjunto de partículas que interagem entre si. Uma partícula pode ser um átomo, uma molécula ou um elétron, um aglomerado de átomos coarse-grained, ou um aglomerado de material mesoscópico ou macroscópico. Os modelos de interação que o LAMMPS inclui são, na sua maioria, de curto alcance; Entretanto, também estão incluídos alguns modelos de longo alcance. Esse pacote de simulações é distribuído gratuitamente pelos laboratórios Sandia, através do link (https://www.lammps.org/). Nesta seção, descreveremos algumas porções do código em LAMMPS usados para a realização deste trabalho.

Quando se inicia uma simulação em DM, faz-se necessário definir algumas propriedades do sistema, por exemplo: as unidades utilizadas, a dimensão do sistema, as condições de contorno (afinal estamos resolvendo equações diferenciais) e o tipo de partículas e/ou moléculas encontradas no sistema. Tais características estão ilustradas na Figura 7. Temos, nas linhas 2 a 6, os comandos que estão transmitindo as informações ao LAMMPS sobre as unidades, dimensões, condições de contorno (neste caso, periódicas) e o tipo atômico molecular, respectivamente.

Figura 7
Introdução do código de simulações.

Em seguida, deve-se informar ao simulador as informações a respeito da caixa de simulação, como os tamanhos laterias L_x, L_y e L_z, o número de partículas (e como as inserir na caixa). Este procedimento está ilustrado na Figura 8. Como se pode observar, as informações a respeito da caixa de simulação estão no comando na linha 1, seguida do comando para criação dessa caixa com as devidas dimensões na linha 2. Estamos centrando a origem da caixa na origem $(\mathrm{0,\; 0,\; 0})$, ou seja, as dimensões para o comando são $-L_x/2$, $L_x/2$, $-L_y/2$, $L_y/2$, $-L_z/2$ e $L_z/2$. Desta forma, teremos uma caixa cúbica com lado L = 20. Outra nota importante: optamos por inserir as 6080 partículas de forma randômica dentro da caixa (na linha 3). As partículas poderiam também ser inseridas, na caixa, em uma rede cristalina pré-definida.

Figura 8
Bloco de características da caixa de simulação e inserção de partículas no sistema.

Agora, devemos informar ao LAMMPS o tipo de interação entre as partículas do sistema. Existe um número grande de formas de interação já incluídas no LAMMPS, como é o caso da interação de Lennard-Jones, cuja a forma é dada pela Eq. 6. Os comandos para esta definição estão ilustrados na Figura 9. Definimos, portanto, uma interação de Lennard-Jones com raio de corte $r_c=2.5\sigma$ através do comando na linha 1, bem como definimos os parâmetros da interação $ϵ$ e $\sigma$ através do comando apresentado na linha 2. Nestas simulações, estamos usando como $ϵ$ e $\sigma$, em unidades reduzidas, como sendo iguais a $1.0$.

Figura 9
Bloco de definição do potencial de interação.

Na Figura 10, encontramos alguns comandos que definem algumas características técnicas da simulação, como por exemplo, a otimização da simulação e parâmetro de incremento de tempo (também chamado de passo temporal). Como otimização, estamos utilizando, nas linhas 1 e 2, os comandos que definem a lista de vizinhos e como esta é modificada ao longo da simulação, respectivamente. Esta forma de otimização cria uma lista de partículas vizinhas para cada par de partículas. Todos os pares de átomos dentro de uma distância de corte de vizinhança igual ao raio de corte (da força) mais uma dada distância, são armazenados na lista (para o nosso caso esta distância é dada por $2.5\sigma +0.5\sigma$). Tipicamente, quanto maior essa dada distância, menos frequentemente as listas de vizinhos precisam ser construídas, contudo mais pares de partículas devem ser verificados para possíveis interações de força a cada passo temporal. O comando de modificação está definindo que essa lista será checada a cada 10 passos de simulação e sem atrasos. Além disso, definimos como incremento de tempo, um passo de $0.001$ em unidades reduzidas de tempo através do comando inserido na linha 3.

Figura 10
Bloco destinado a características técnicas e de otimização da simulação.

Em seguida, definimos algumas características adicionais ao sistema ilustradas na Figura 11. Estamos agrupando as partículas do sistema ao grupo “tracer” através do comando na linha 1. Além de definir a massa das partículas $m=1.0$, em unidades reduzidas, através do comando na linha 2. Como também, minimizamos a energia do sistema através do comando dado pela linha 3, com o objetivo de evitar sobreposições de partículas, visto que elas foram introduzidas aleatoriamente dentro da caixa.

Figura 11
Bloco de características adicionais da simulação.

Para finalizar o código, temos de informar ao LAMMPS qual será o comportamento dinâmico do sistema em análise. Desta forma, serão executadas (i) duas etapas de equilibração termodinâmica (cada uma com um tamanho distinto da caixa) e, após esta equilibração, (ii) mais passos de simulação para a tomada de resultados, na etapa chamada de produção. Todos esses comandos estão ilustrados na Figura 12. Nas linhas de código 1 e 2, encontram-se as informações sobre a dinâmica do sistema. Como escrito previamente, estamos realizando a DM à temperatura constante e fazendo uso do termostato de Langevin. Nestas linhas, encontram-se os comandos para implementar tal termostato. Em seguida, nas linhas 4 a 7, encontram-se os comandos para imprimir – na tela onde será executada a simulação – informações importantes a cada intervalo de passos, como: a energia potencial, a energia cinética, a temperatura, a pressão, a densidade e os lados da caixa, além do comando run, o qual faz com que o LAMMPS execute a primeira etapa da simulação (de equilibração da fase densa líquida), antes de aumentarmos a caixa no eixo z. Os comandos para aumento da caixa, encontram-se nas linhas 9 e 10. Após esse procedimento, realizamos uma outra equilibração à temperatura desejada. Feito todos estes procedimentos, está na hora de calcularmos o perfil de densidade ao longo do eixo z. Nas linhas 14 a 17, encontram-se os comandos relacionados à computação dessa propriedade, além do comando para rodar $1\times {10}^6$ passos para a produção de resultados.

Figura 12
Bloco de características dinâmicas da simulação, equilibração, e tomada de resultados.

5. Observando a Coexistência e o Diagrama de Fases

As fotografias da simulação, que foram obtidas utilizando as trajetórias das simulações realizadas e a versão gratuita do pacote de visualização Ovito [³⁰], estão ilustradas na Figura 13. Como se pode observar, no painel superior, temos o sistema na menor temperatura simulada $T=0.5$ com a fase líquida densa centrada e apenas algumas partículas povoando o espaço vazio. Nota-se que, à medida que aumentamos a temperatura do sistema (ver painéis inferiores), estamos fornecendo maior energia cinética, e desta forma possibilitando que mais partículas saiam da fase líquida e comecem a povoar os espaços vazios. Este efeito gera uma interface entre as fases líquida e vapor. Com o objetivo de se observar o efeito da temperatura sobre o sistema a a formação da coexistência líquido-vapor, um vídeo desta simulação foi criado e pode ser encontrado clicando-se aqui.

Figura 13
Zoom dos instantâneos da simulação entre as distâncias $-70\sigma$ e $70\sigma$ na direção z. Do painel superior ao inferior, temos os instantâneos do sistema nas temperaturas $T=0.50$, $T=0.80$, $T=0.90$ e $T=1.00$.

Ao se calcular o perfil de densidade $\rho (z)$ ao longo do eixo z, conseguimos apreciar o efeito da temperatura sobre o sistema. Este perfil está ilustrado na Figura 14(a) para as diversas temperaturas simuladas. Como se pode observar, temos na menor temperatura $T=0.5$ a fase líquida com maior densidade (ver curva azul) preenchendo o espaço entre $z=-10$ e z = 10. À medida que aumentamos a temperatura, mais partículas passam a povoar os espaços vazios, então, a densidade da fase líquida começa a diminuir. Isto pode ser observado pelo alargamento das extremidades da curva de densidade e diminuição do pico à medida que aumentamos a temperatura do sistema. Por exemplo, ao observar a curva para temperatura $T=1.0$, temos um pico bem menor em ralação a $T=0.5$ e as extremidades maiores e não nulas. Tal efeito de alargamento indica o aumento da densidade de partículas nos espaços vazios, ou seja, onde temos a fase vapor. Pode-se observar melhor esse efeito tomando-se apenas uma das interfaces, na região positiva, como ilustrado na Figura 14 (b).

Figura 14
(a) Perfis de densidade $\rho (z)$ ao longo do eixo z para as diversas temperaturas simuladas. (b) O mesmo gráfico, mas apenas tomando a parte positiva de z. Ambos os gráficos representam resultados obtidos diretamente da simulação para as diversas temperaturas simuladas. Para a obtenção da densidade da fase líquida ${\rho }_l$, densidade da fase vapor ${\rho }_v$, posição da interface z₀ e espessura da interface d recorre-se ao ajuste de curva via a equação (8). Esses resultados estão ilustrados na Tabela 2.

Para obter o diagrama de fases deste sistema, ajustamos os dados ilustrados na Figura 14 (b) com a função sigmoide, dada pela equação 8, para cada temperatura simulada, de forma a, se obter os valores da densidade da fase líquida ${\rho }_l$, da densidade da fase vapor ${\rho }_v$, da posição da interface z₀ e, finalmente, da espessura dessa interface d. Os resultados obtidos através desse ajuste estão ilustrados na Tabela 2. Em seguida, para se obter o valor da densidade crítica, tomamos a diferença entre ${\rho }_l$ e ${\rho }_v$ em função de cada temperatura T simulada e usando o expoente crítico experimental $\beta =0.35$, ajustamos esses dados usando a equação (2). Assim, obtivemos o valor da temperatura crítica do sistema, dada por $T_c=1.171$. Por outro lado, tomando a média aritmética entre ambas as densidades líquida e vapor em função de cada temperatura simulada, mas agora usando T_c - T onde T_c é a temperatura crítica obtida anteriormente, ajustamos esses valores através da equação (3), e foi possível obter o valor da densidade crítica ${\rho }_c=0.277$ do sistema. Ao compararmos os valores obtidos com o valor indexado pela literatura científica [⁴⁰], temos uma diferença entre nossos valores encontrados e os dados pela literatura de apenas $10.75\%$ para a temperatura crítica T_c e $10.35\%$ para a densidade crítica ${\rho }_c$. Tendo em mente que nosso sistema é finito, e, assim sendo, essas discrepâncias podem estar relacionadas com efeitos de tamanho finito. A influência de tamanho finito é uma boa sugestão para ser explorada por um professor de graduação em cursos de Mecânica Estatística, pois só se obtém o valor exato dessas quantidades no limite termodinâmico, ou seja, quando o volume da caixa tender ao infinito. Voltando aos valores obtidos para a temperatura e densidade críticas, T_c e ${\rho }_c$, respectivamente, o ponto dado por esses valores, no diagrama de fases, significa o último valor de temperatura e densidade para os quais ainda podemos encontrar a coexistência entre as fases líquida e vapor. Acima deles, o sistema torna-se indistinguível.

Tabela 2
Resultados obtidos após o ajuste de curva dado pela equação 8 para cada temperatura e seus valores de densidade da fase líquida ${\rho }_l$, densidade da fase vapor ${\rho }_v$, posição da interface z₀ e espessura da interface d.

T	${\rho }_l$	${\rho }_v$	z0	d
0.50	0.8822	0.0000	8.7153	0.6627
0.55	0.8594	0.0002	8.8549	0.7132
0.60	0.8364	0.0006	9.0403	0.8002
0.65	0.8124	0.0013	9.2125	0.9144
0.70	0.7868	0.0027	9.2152	1.0555
0.75	0.7599	0.0049	9.1400	1.1832
0.80	0.7313	0.0077	9.1182	1.3332
0.85	0.7006	0.0122	8.8688	1.5974
0.90	0.6669	0.0165	8.7127	2.0479
1.00	0.6080	0.0299	7.6881	4.1969

Com todos esses resultados em mãos, podemos montar o diagrama de fases $T\times \rho$, o qual se encontra ilustrado na Figura 15. Podemos observar a região de coexistência líquido-vapor abaixo da curva vermelha, onde temos as letras indicativas de cada fase (V+L). Não somente, ilustramos através dos pontos em azul, a média aritmética das densidades das fases líquida e vapor, segundo a lei dos diâmetros retilíneos dada pela equação 3. Além da sua extrapolação para se obter o ponto crítico (${\rho }_c$, T_c). É possível notar que, para valores de densidade e temperatura acima desse ponto, não há mais coexistência, pois a partir dele as fases se tornam indistinguíveis.

Figura 15
Diagrama de fase $T\times \rho$ obtido para esse sistema, cuja densidade global é $\rho =0.076$. As letras V e L correspondem as fases vapor e líquida, respectivamente. V+L a região de coexistência de ambas as fases (área vermelha). Além de apresentar os pontos críticos para a densidade ${\rho }_c$ e temperatura T_c. A curva azul não possui significado físico, mas representa a média aritmética das densidades das fases líquida e vapor, segundo a lei dos diâmetros retilíneos dada pela equação 3.

6. Considerações Finais

Neste estudo, apresentamos uma abordagem acessível (e economicamente viável) para a construção de diagramas de fases por meio de simulações de Dinâmica Molecular aplicadas a um sistema de moléculas interagentes de acordo com o potencial de Lennard-Jones. Além disso, apresentamos o software de simulação LAMMPS, disponível gratuitamente na internet e compatível com sistemas Linux, Windows e Mac.

A metodologia apresentada permite obter, através da análise de perfis de densidade para cada temperatura simulada, informações cruciais sobre o sistema, tais como temperatura crítica, densidade crítica, localização da interface entre as fases líquida e de vapor e a espessura dessa interface. Nossos resultados para o ponto crítico estão em consonância com os valores reportados na literatura científica, validando a precisão do nosso método. Isso nos permitiu identificar a região de coexistência entre as fases líquida e de vapor no diagrama de fase, e visualizar como ocorre essa coexistência.

Com este estudo e tendo em mente o constante avanço da tecnologia, esperamos incentivar professores de Física em todos os níveis de ensino a aproveitar, de forma gratuita, uma ferramenta valiosa para o ensino das transformações termodinâmicas. Isso proporciona um laboratório instrutivo para a construção de diagramas de fases, onde é possível explorar tanto os detalhes da construção desses diagramas, quanto a teoria das transições de fase e fenômenos críticos, adaptados a cada nível de ensino. A simplicidade do modelo permite que as simulações sejam realizadas em computadores pessoais comuns, ampliando o acesso e beneficiando uma ampla gama de estudantes. Obviamente, com mais acesso a recursos computacionais, é possível visualizar o comportamento de modelos mais complexos, como moléculas poliatômicas e mesmo macromoléculas. Isso, por sua vez, contribui para a disseminação do conhecimento e o aprimoramento do ensino da termodinâmica.

Agradecimentos

Sem financiamento público esta pesquisa seria impossível. Os autores são gratos ao Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, proc. 405479/2023-9, 441728/2023-5, e 304958/2022-0), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES, código de financiamento 0001), Fundação de Amparo à Pesquisa do Estado do Rio Grande do Sul (FAPERGS, TO 21/2551-0002024-5), pelo apoio financeiro.

Referências

• ^[1]
  H.B. Callen, Thermodynamics and an introduction to thermostatistics (Wiley, New York, 1985), 2 ed.
• ^[2]
  M.J. Oliveira, Termodinâmica (Livraria da Física, São Paulo, 2005).
• ^[3]
  J.P. Ansermet e S.D. Brechet, Principles of Thermodynamics (Cambridge University Press, Cambridge, 2019).
• ^[4]
  A.R.E. Oliveira, A History of the Work Concept: From Physics to Economics (Springer, Dordrecht, 2013).
• ^[5]
  W. Greiner, D. Rischke, L. Neise e H. Stöcker, Thermodynamics and Statistical Mechanics (Springer, New York, 2012).
• ^[6]
  I. Müller, A History of Thermodynamics: The Doctrine of Energy and Entropy (Springer, Berlin, 2007).
• ^[7]
  D. Kondepudi e I. Prigogine, Modern Thermodynamics: From Heat Engines to Dissipative Structures (CourseSmart, Wiley, 2014).
• ^[8]
  Y.A. Chang, Metallurgical and Materials Transactions A 37, 273 (2006).
• ^[9]
  M.L.L. Farias e M.A.A. Barbosa, Revista Brasileira de Ensino de Física 39, e4402 (2017).
• ^[10]
  N.O. Smith, Journal of Chemical Education 35, 125 (1958).
• ^[11]
  K. Doymus, Journal of Chemical Education 84, 1857 (2007).
• ^[12]
  B. Alder e T. Wainwright, Journal Of Chemical Physics 31, 459 (1959).
• ^[13]
  N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller e E. Teller, The Journal Of Chemical Physics 21, 1087 (1953).
• ^[14]
  D.C. Rapaport, The Art of Molecular Dynamics Simulation (Cambridge University Press, Cambridge 2004), 2 ed.
• ^[15]
  A. Pfennig, Journal of Foreign Language Education and Technology 3, 1 (2018).
• ^[16]
  A.M. Almeida, N.A Pimenta, K.S. Nascimento, L.C Gontijo e S.P. Eiras, Conexão Ciência (Online) 9, 16 (2015).
• ^[17]
  J.M. Boblick, Science Education 54, 77 (1970).
• ^[18]
  J. Richards, W. Barowy e D. Levin, Journal of Science Education and Technology 1, 67 (1992).
• ^[19]
  F. Esquembre, Computer Physics Communications 147, 13 (2002).
• ^[20]
  C.E. Wieman, W.K. Adams e K.K. Perkins, Science 322, 682 (2008).
• ^[21]
  M.S. Pfefferová, Journal of Technology and Information Education 7, 81 (2015).
• ^[22]
  Y. Li, L. Ma e Y. Shi, em: Education and Educational Technology, editado por Y. Wang (Springer, Berlin, 2012).
• ^[23]
  J. Abiasen e G. Reyes, International Journal of Asian Education 2, 480 (2021).
• ^[24]
  F. Almasri, Education and Information Technologies 27, 7161 (2022).
• ^[25]
  H.J. Banda e J. Nzabahimana, Phys. Rev. Phys. Educ. Res. 17, 023108 (2021).
• ^[26]
  J.E Jones, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 106, 463 (1924).
• ^[27]
  E. Guggenheim, The Journal Of Chemical Physics 13, 253 (1945).
• ^[28]
  D. Frenkel e B. Smit, Understanding molecular simulation: from algorithms to applications (Elsevier, San Diego, 2001), v. 1.
• ^[29]
  M.P. Allen e D.J. Tildesley, Computer Simulation of Liquids (Oxford University Press, New York, 2017), 2 ed.
• ^[30]
  A. Stukowski, Modelling and Simulation in Materials Science and Engineering 18, 015012 (2010).
• ^[31]
  J.E. Lennard-Jones, Proceedings of the Physical Society 43, 461 (1931).
• ^[32]
  H. Watanabe, N. Ito e C.K. Hu, The Journal of Chemical Physics 136, 204102 (2012).
• ^[33]
  F.J Blas, L.G. MacDowell, E. Miguel e G. Jackson, The Journal of Chemical Physics 129, 144703 (2008).
• ^[34]
  K.S. Silmore, M.P. Howard e A.Z. Panagiotopoulos, Molecular Physics 115, 320 (2017).
• ^[35]
  M. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation (Oxford Graduate Texts, Oxford, 2010).
• ^[36]
  H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola e J.R. Haak, The Journal of Chemical Physics 81, 3684 (1984).
• ^[37]
  S. Nosé, The Journal of Chemical Physics 81, 511 (1984).
• ^[38]
  W.G. Hoover, Physical Review A 31, 1695 (1985).
• ^[39]
  S. Plimpton, Journal of Computational Physics 117, 1 (1995).
• ^[40]
  D.M. Heyes, CMST 21, 169 (2015).

• 1
  Sugerimos ao leitor interessado, buscar a Fig. 8.1 da Ref. [²] para uma melhor compreensão desses patamares na região de coexistência das isotermas no espaço p-V
• 2
  Na literatura especializada, é comum nos referirmos à energia potencial de interação entre as partículas $U(r)$ como sendo apenas o potencial de interação entre as partículas.

Datas de Publicação

• Publicação nesta coleção
  16 Ago 2024
• Data do Fascículo
  2024

Histórico

• Recebido
  06 Nov 2023
• Revisado
  08 Jul 2024
• Aceito
  08 Jul 2024