Simulação Monte Carlo com Repesagem Aplicada ao Calor Específico de Sólidos

Simulação Monte Carlo com Repesagem Aplicada ao Calor Específico de Sólidos

Roberto da Silva^* e J. R. Drugowich de Felício^†

Departamento de Física e Matemática

Faculdade de Filosofia, Ciências e Letras de Ribeirão Preto - USP

Av. Bandeirantes, 3900, 14040-901, Ribeirão Preto, São Paulo, SP, Brasil

Recebido em 9 de abril, 2002. Aceito em 12 de junho, 2002.

Apresenta-se uma estimativa numérica para a contribuição das vibrações da rede para o calor específico de um sólido cristalino. O método Monte Carlo é utilizado tanto em sua forma original quanto em sua versão mais moderna, que emprega a técnica da repesagem para extrair mais informações de cada simulação. Mostra-se como obter curvas contínuas para grandezas como a energia e o calor específico em função da temperatura a partir de simulação feita em apenas uma temperatura. Os resultados analíticos são incluídos para facilitar a comparação com os resultados da simulação.

A numerical estimate for the contribution of the lattice to the specific heat of a crystalline solid is presented. The Monte Carlo method is used such in its original form as in its modern version which explores the reweighting technique to extract more information from each simulation. One reveals how to get continuous plots for quantities like energy and specific heat, in function of the temperature, performing just one simulation at a convenient temperature. For the sake of comparison we present also the pertinent analytical results.

I Introdução

Nos últimos anos, avançou bastante o conhecimento a respeito de sistemas de muitas partículas que não podem ser descritos de forma analítica. O progresso decorre, em grande parte, da rapidez com que os computadores realizam os cálculos necessários para uma simulação do tipo Monte Carlo, por exemplo. Aqui, quando falamos em evolução, estamos nos referindo ao progresso obtido na eletrônica (''hardware''). Houve, porém, um enorme avanço também na forma de extrair resultados das simulações. Nesse caso, a evolução ocorreu na programação (''software''). Resumidamente, um conjunto de operações que se repetia uma centena de vezes passou a ser feita de uma única vez, seus resultados dando origem a uma série de outros apenas por um artifício de ''repesagem'', que será explicado na seqüência.

O objetivo principal desse artigo é introduzir essa nova maneira de fazer simulações, partindo de um sistema simples como é o cristal de Einstein (a primeira aproximação ao calor específico de sólidos que prevê o limite correto quando T ® 0). Para esse sistema, todos os resultados exatos são exibidos na próxima seção, podendo ser acompanhados por alunos de um curso introdutório de mecânica estatística. Na seção III, explicamos a técnica de simulação conhecida como algorítmo de Metropolis ou Monte Carlo com amostragem por importância. São apresentadas estimativas para a energia e o calor específico em algumas temperaturas e comparadas com os resultados da seção II. Na seção IV, a técnica do histograma com repesagem é apresentada e os resultados da simulação, feita em apenas uma temperatura, são utilizados para obter curvas contínuas da energia e do calor específico em um intervalo da temperatura. Alguns comentários são feitos para ilustrar como essa abordagem tem evoluído nos últimos anos.

II O modelo

Considerar um sólido cristalino como um conjunto de osciladores independentes, cada um deles com energia kT, de acordo com o teorema da equipartição, conduz a uma capacidade térmica igual a 3N₀k = 3R cal/mol, independente da temperatura. Esse é o resultado da aproximação de Dulong-Petit que vale para sólidos à temperatura ambiente. No entanto, quando a temperatura é reduzida, os resultados experimentais diferem dessa constante e tendem a zero quando T ® 0. Para corrigir a distorção da curva teórica, Einstein propôs, em 1907 [1], a quantização da energia dos átomos que vibram em torno dos sítios da rede. Seguindo uma idéia usada antes por Planck (1900) para tratar a radiação do corpo negro, ele imaginou que apenas alguns valores seriam permitidos para a energia de cada oscilador, o que ele escreveu como

onde n é um inteiro não negativo, h é a constante de Planck e n a freqüência. Com essa hipótese para a energia e tratando as vibrações como sendo independentes - todas com a mesma freqüência - Einstein chegou à expressão

para a contribuição da rede com o calor específico do sólido. Para calcular essa grandeza, basta lembrar que, em contato com um banho térmico (ensemble canônico), a probabilidade para um desses osciladores ocupar um estado com energia e decresce exponencialmente com a energia (fator de Boltzmann). Assim, com

e e dada pela equação (1), teremos para a energia média de cada oscilador a expressão

onde b = 1/kT,T é a temperatura e k a constante de Boltzmann. O produto de b por hn que aparece nas exponenciais costuma ser representado por q_E/T , onde q_E = hv/k é conhecida como a ''temperatura'' de Einstein.

Observando que

e que o somatório entre parênteses é uma progressão geométrica de razão q = e^-bhn, concluímos que

e, portanto,

De maneira análoga, pode-se obter o valor médio de e² e, finalmente, através da seguinte expressão

chegar à expressão (2), lembrando que o fator 3 é proveniente dos três graus de liberdade para os átomos da rede.

O cálculo da energia média do conjunto de osciladores pode ser feito lançando-se mão da expressão

onde a soma não é feita sobre os estados quânticos mas sim sobre os níveis de energia. Como a cada nível de energia estão associados diferentes estados dos osciladores, é preciso multiplicar o fator de Boltzmann pela degenerescência W do nível de energia E = n_ie_i, cujo valor exato é [2]

onde K = E/hn.

É importante observar que W(E) é uma função crescente da energia, enquanto o fator de Boltzmann é uma função que decresce exponencialmente com E. O produto de W(E) pelo fator exp(-bE), daqui por diante chamada P(E), tem valores próximos de zero quando a energia é baixa (W é pequeno) e também quando a energia é alta (a exponencial vai assintoticamente a zero para E ® ¥). O produto só tem valores apreciáveis em uma faixa de valores de energia que engloba a energia média. Essa faixa é tanto menor quanto maior for o número de osciladores. A Fig. 1 mostra o comportamento da probabilidade P(E) de encontrar o conjunto de 3N osciladores com energia igual a E, quando o sistema está em contato com um banho térmico à temperatura q_E/T = bhn = 1/3.

III O algoritmo de Metropolis

Para obter numericamente os resultados da seção anterior, vamos primeiramente fazer uma simulação Monte Carlo tradicional [3]. Essa técnica tem sido utilizada em uma gama de situações que inclui minimização de funções (método que ficou conhecido como ''recozimento simulado'' [4]), o cálculo de integrais multidimensionais [5], transições de fase estruturais em sólidos [6], além de outras aplicações [7]. Iniciamos a simulação com 3N osciladores (300 já é um número razoável). A energia de cada oscilador é igual a um número inteiro de quantas de energia hn que expressaremos por e_n/hv = n onde n Î IN (a menos da energia de ponto zero, hn/2). No início da simulação, as energias (em unidades de hn) são escolhidas aleatoriamente entre 0 e um teto n_MAX que estipulamos. Verificamos que 300 é um bom número para esse limite. A partir desse ponto, vamos evoluir o conjunto de 3N osciladores durante 2000 passos que servirão para termalizar o sistema. Como fazer isso? Bem, um novo nível é sorteado para cada oscilador e, na seqüência, comparado com o nível anterior. Se o nível sorteado for menor que o anterior (a energia decresce), o novo nível é aceito com probabilidade 1. Em caso contrário (De > 0), um número aleatório entre 0 e 1 é sorteado e comparado com a razão entre os fatores de Boltzmann dos níveis velho e novo (probabilidade de transição do velho para o novo nível). Se o número aleatório for menor que exp(-bDe), o novo estado é aceito. Se for maior, o oscilador continua com a energia anterior. Um passo de Monte Carlo é concluído quando todos os osciladores tiverem sido visitados uma vez (se eles forem escolhidos de forma determinística) ou quando 3N osciladores tiverem sido visitados aleatoriamente (caso em que cada oscilador também será visitado uma vez por passo, na média). Esse processo imita a troca de energia entre o sistema e um banho térmico à temperatura T e, de acordo com o princípio do balanceamento detalhado (veja [8]), depois de várias iterações desse tipo os estados serão visitados com probabilidade proporcional ao fator de Boltzmann exp(-be). Esse algoritmo recebe o nome de Metropolis, em homenagem ao primeiro autor de um trabalho publicado em 1953 [9], cuja originalidade estava justamente em propor a amostragem por importância. Em poucas palavras, trata-se de um método que despende pouco tempo visitando estados de baixa probabilidade, como é o caso dos estados cujas energias estejam fora da região onde ocorre o pico da Fig. 1.

Terminados os 2000 passos de Monte Carlo de termalização, o que significa que cada átomo foi visitado e teve a possibilidade de trocar o seu estado 2000 vezes, começamos a medir a energia média áEñ e também o seu segundo momento, áE²ñ . Para isso, 18000 passos serão executados, seguindo o mesmo procedimento anterior (resumido no quadro 1).

A diferença é que agora, ao final de cada passo, a energia do sistema será acumulada em uma variável ETOT, que servirá para estimar a energia média. Também o quadrado da energia será armazenado em uma outra variável EQUAD, que servirá para estimar a energia quadrática média, após o término dos 18000 passos. Como os estados pelos quais o sistema passa obedecem à distribuição de Boltzmann (o princípio do balanceamento detalhado garante isso), o cálculo da energia média pode ser feito como simples média aritmética das energias obtidas ao final de cada passo. O mesmo se aplica ao cálculo da energia quadrática média. Dessa forma, a conhecida fórmula para o cálculo da média da n-ésima potência da energia

no ensemble canônico, é substituída por uma simples média aritmética quando usamos o algoritmo de Metropolis (simulação por amostragem de importância):

onde N_MC é o número de passos de Monte Carlo válidos da simulação (número total menos o número de passos utilizados para a termalização).

Na tabela 1, apresentamos os valores para a energia média e para o calor específico para dois valores de temperatura, utilizando o método da amostragem por importância. O tempo de cada simulação foi de 5 segundos em um Pentium III de 933 MHz e 256 Mb de memória RAM, usando-se a linguagem Fortran [10]. Os erros foram estimados executando-se o programa com 5 diferentes sementes, o que significa que foram obtidos 5 pontos para o cálculo da média e do desvio-padrão da média para cada temperatura. Os valores são comparados com resultados exatos discutidos na seção anterior.

A dificuldade com esse procedimento é que para levantar uma curva de calor específico, por exemplo, há necessidade de repetir a simulacão para uma série de valores da temperatura. Na próxima seção, nós vamos mostrar como é possível aperfeiçoar a técnica e evitar a repetição do processo. Para isso, será preciso gerar um histograma que registre o número de vezes que cada energia foi visitada, durante os N_MC passos da simulação.

IV Histograma e a técnica de repesagem

Nos últimos anos, ganhou força uma abordagem que visa extrair mais informação das simulações. A idéia [11], já proposta nos anos 60 ([12],[13]) , consiste em separar o efeito de temperatura (que é um parâmetro da simulação) dos efeitos entrópicos, característicos do sistema e independentes da temperatura. Explicando melhor, quando se realiza uma simulação amostrando os estados por importância (usando o algoritmo de Metropolis, por exemplo), o número de vezes que cada estado aparece obedece à distribuição de Boltzmann. Porém o número de vezes que uma dada energia é encontrada será proporcional ao produto do fator de Boltzmann pela degenerescência daquele nível de energia W(E). Isso acontece porque a cada nível de energia estão associados diversos estados quânticos do sistema compatíveis com ela e, quanto maior esse número, maior será a chance daquela energia aparecer. Dessa forma, o que se obtém em uma simulação com amostragem por importância, na qual se constrói o histograma H_b(E), é o número de vezes que o sistema apresenta energia igual a E em N_MCpassos de Monte Carlo. Esse número é proporcional à degenerescência e contém o fator de Boltzmann, de forma que :

Deve-se ver, portanto, o termo H_b(E)/N_MC como uma estimativa para a probabilidade de encontrar o sistema com uma dada energia E, que no ensemble canônico é usualmente escrita como :

Na Fig. 2, esboçamos o histograma H(E) obtido em bhn = 1/3. O histograma foi obtido com N_MC = 20000 passos de Monte Carlo, dispensando 2000 passos para termalização, e se compara bem com o resultado exato. Podemos notar que a curva tem a mesma forma da Fig. 1 (linha cheia) e apresenta um pico na região da energia média (nesse exemplo, áEñ = 300(e^1/3-1)^-1erg » 758 erg).

Na seqüência, o que se faz é ''repesar'' o histograma para obtê-lo em uma outra temperatura T' = 1/kb', sem precisar repetir a simulação. Ressalte-se que esse procedimento tem limitações e em alguns casos é recomendável obter o histograma em três ou quatro temperaturas para depois chegar à melhor estimativa para W(E). Voltaremos a essa discussão no final desta seção, mas, por enquanto, vamos apenas mostrar como proceder para estender os resultados de um histograma simples para outros valores de temperatura. Para isso, basta seguir o raciocínio de Ferrenberg e Swendsen [11], observando-se que para um b^'¹ b, tendo em vista a expressão (13), podemos escrever

de forma que grandezas como a energia média e seus momentos possam ser calculadas nessa outra temperatura, utilizando-se simplesmente a equação

Convém ressaltar que essa técnica foi utilizada com sucesso em uma coleção de problemas de mecânica estatística [14] [15]. Na Fig. 3, esboçamos o calor específico e a energia média dos átomos, em função da temperatura, utilizando o método do histograma. Nesse gráfico, a simulação foi executada em bhn = 1.0 e a repesagem foi feita para as temperaturas T/q_E = (bhn)^-1 no intervalo [0.1,1.0]. Podemos constatar que, para altas temperaturas, a energia varia linearmente com a temperatura, o que implica um calor específico constante (lei de Dulong-Petit), como realmente se espera nesse limite. A curva obtida com a repesagem coincide com a curva teórica em todo o intervalo e apresenta, na região em que T ® 0, a mesma tendência que o resultado analítico. Vale ressaltar que a concordância entre as curvas não será boa se o histograma for construído em uma temperatura mais baixa, pois nesse caso níveis mais altos de energia dificilmente serão visitados. Para determinar objetivamente o intervalo de temperaturas em que a repesagem é válida [16], é preciso examinar a variação da energia média com a temperatura e exigir que os extremos do intervalo estejam associados a energias que estejam dentro de E ± 2s onde s é o desvio quadrático médio da distribuição P(E).

O passo seguinte é extrair apenas as informações que não dependam da temperatura. Essa é a tendência dos trabalhos mais recentes nos quais a busca é pela densidade de estados, com a qual todas as grandezas podem ser calculadas, seja qual for a temperatura. No presente caso, basta multiplicar o histograma H_b(E) por exp(bE) para eliminar o fator de Boltzmann (único que depende da temperatura) e ficar apenas com a densidade de estados W(E), fator entrópico que é independente da temperatura. Na Fig. 4, esboçamos um gráfico de W em função de E, obtido dessa maneira. Através da expressão (10), é possível também obter a curva teórica para W(E) e compará-la com a estimativa numérica. Note-se que a concordância é razoável em todo o intervalo de energia que foi explorado nessa simulação.

Uma outra variante que pode ser explorada nesse caso refere-se à dependência exclusiva com a temperatura. Para isso, devemos gerar um histograma da energia individual dos osciladores. Procedendo assim, não estaremos levando em conta a degenerescência dos níveis e mediremos, portanto, apenas a probabilidade de encontrar um oscilador com energia e. O novo histograma corresponde simplesmente ao fator de Boltzmann (conforme equação 3). O resultado dessa investigação, já dividido por uma constante adequada, está mostrado na Fig. 5, onde também está representada a curva teórica exp(-be) com bhn = 1/10.

Finalmente, alguns comentários sobre a validade dessa técnica e as inovações dos últimos anos são pertinentes. Há uma limitação da técnica do histograma que é imposta pela ausência de passagens pelos estados em que a probabilidade é muito baixa (asas da distribuição). É por essa razão que o intervalo em que fazemos a repesagem não pode abranger temperaturas que conduzam a energias médias fora da região onde P_b(E) seja apreciável. Também é preciso observar que o próprio intervalo em que a distribuição de probabilidades é apreciável vai ficando menor, quando o número de partículas (N) do sistema aumenta. Sendo assim, a extensão dos resultados para o caso de sistemas maiores deve necessariamente ser feita dentro de um intervalo menor de temperatura. Essa limitação fez com que outras abordagens fossem propostas para calcular a densidade de estados. Uma delas foi introduzida por Ferrenberg e Swendsen [17],[18]. Trata-se da técnica do histograma múltiplo, construído iterativamente a partir de simulações feitas em uma coleção de temperaturas. Uma extensão inteligente desse método foi proposta por Berg e Neuhaus [19] e ficou conhecida como ensemble multicanônico. Ela consiste em ''experimentar'' cada intervalo de energia com um banho que privilegie o aparecimento de estados daquela faixa. Com isso, a probabilidade de encontrar os estados daquele intervalo torna-se apreciável e de mesma ordem de grandeza que os estados de outros intervalos. É como se substituíssemos a distribuição da Fig. 1 por uma outra que fosse um platô. Naturalmente que será necessário retornar à verdadeira distribuição de probabilidades, artificialmente modificada pelos diversos banhos que compunham o ensemble multicanônico. Essa técnica tem sido utilizada com relativo sucesso mas é de difícil implementação [20]. Uma outra abordagem foi proposta por P. M. C. de Oliveira, T. J. Penna e H. J. Herrmann [21] e ficou conhecida como ''broad-histogram''. Essa nova abordagem e suas generalizações têm passado por diversos testes para verificar se ela obedece ao balanceamento detalhado. Tudo indica que, além de um interessante estudo de casos no tocante a métodos de amostragem estatística, ela será uma ferramenta de pesquisa muito útil ( [8],[22], [23]).

V Conclusões

Utilizamos neste texto simulações de Monte Carlo para calcular a energia da rede e a sua contribuição para o calor específico de um sólido isolante, de acordo com a aproximação de Einstein. Trabalhando primeiramente com o método da amostragem por importância (algoritmo de Metropolis), fomos capazes de encontrar valores tanto para a energia quanto para o calor específico, para alguns valores da temperatura. Os resultados encontrados numericamente comparam-se muito bem com os resultados exatos.

Em seguida, introduzimos o método da repesagem para obter curvas contínuas para o calor específico e para a energia em função da temperatura. Elas estão em completo acordo com os resultados analíticos apresentados na seção II. Vale a pena ressaltar que a extensão pela técnica da repesagem no presente caso deve ser feita para temperaturas menores do que a da simulação (valores de b maiores). Isso porque, fazendo-se a simulação com uma temperatura mais alta, os diversos níveis de energia podem ser ocupados durante a simulação. Verificamos pela curva da energia um crescimento linear com a temperatura, no limite de altas temperaturas. É esse comportamento que leva a um calor específico constante nesse mesmo limite, conforme estabelece a Lei de Dulong-Petit. Mais interessante, no entanto, foi observar a variação do calor específico com a temperatura, conforme preconiza a teoria de Einstein e garante a experiência. Os alunos saberão explorar esse modelo de brinquedo para testar as idéias fundamentais da mecânica estatística. Ele tem a vantagem de não ter limitação para os níveis de energia, como ocorre com sistemas de spins.

Agradecimentos

Agradecemos à FAPESP, pelo apoio financeiro. JRDF agradece também a Jorge Chahine, Nelson A. Alves, Nestor Caticha e Valter Líbero pelas discussões esclarecedoras sobre o assunto.

References

[1] A. Einstein, Ann. Physik, 22, 186 (1907).

[2] S. R. A. Salinas, Introdução à Física Estatística, Edusp (1997).

[3] V. L. Líbero, Revista Brasileira de Ensino de Física, 22, 346 (2000).

[4] S. Kirkpatrick, C. D. Gelatt, M. P. Vecchi, Science 220, 671 (1983), J. Stat. Phys. 34, 975 (1984).

[5] H. Gould, J. Tobochnik, An Introduction to Computer Simulation Methods, Addison Wesley (1996).

[6] G. Meissner, K. Binder, Phys. Rev. B 12, 3948(1975).

[7] M. E. J. Newman, G. T. Barkema, Monte Carlo Simulation in Statistical Physics, Clarendin Press-Oxford (1999).

[8] D. P. Landau e K. Binder, A guide to Monte Carlo simulations in Statistical Physics, Cambridge University Press (2000).

[9] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller e E. Teller, J. Chem. Phys. 21, 1087 (1953).

[10] Os programas relativos e este artigo podem ser obtidos com os autores.

[11] A. M. Ferrenberg e R. H. Swendsen, Phys. Rev. Lett. 61, 2635 (1988).

[12] D. A. Chesnut e Z. W. Salsburg, J. Chem. Phys. 38, 2861 (1963).

[13] G. M. Torrie e J. P. Valleau, J. Comp. Phys. 23, 187 (1977).

[14] R. H. Swendsen,A. M. Ferrenberg, J. S. Wang, Top. Apl. Phys, 71, 75 (1995).

[15] N. A. Caticha, J. Chahine e J. R. Drugowich de Felício, Phys. Rev. A 40, 7431 (1989); Phys. Rev. B 43, 1173 (1991).

[16] N. A. Alves, Phys. Rev. D 46, 3678 (1992).

[17] A. M. Ferrenberg e R. H. Swendsen, Phys. Rev. Lett. 63, 1196 (1989).

[18] J. R. Drugowich de Felício e V. L. Líbero, Am. J. Phys. 64, 1281 (1996).

[19] B. A. Berg, T. Neuhaus, Phys. Rev. Lett. 68 9 (1992).

[20] N. A. Alves, J. R. Drugowich de Felício, U. H. E. Hansmann, J. Phys. A 33, 7489 (2000).

[21] P. M. C de Oliveira, T. J. Penna e H. J. Herrmann, Braz. J. Phys. 26, 677 (1996).

[22] A. R. Lima, P.M.C. de Oliveira e T. J. Penna, Solid State Communications 114, 447 (2000), J. Stat. Phys. 99, 691 (2000).

[23] P. M. C. de Oliveira, Braz. J. Phys. 30, 766 (2000) e 195 (2000).

^† [e-mail]

Datas de Publicação

• Publicação nesta coleção
  01 Dez 2006
• Data do Fascículo
  Jun 2002

Histórico

• Aceito
  12 Jun 2002
• Recebido
  09 Abr 2002