Uma abordagem simplificada do método Monte Carlo Quântico: da solução de integrais ao problema da distribuição eletrônica

UMA ABORDAGEM SIMPLIFICADA DO MÉTODO MONTE CARLO QUÂNTICO: DA SOLUÇÃO DE INTEGRAIS AO PROBLEMA DA DISTRIBUIÇÃO ELETRÔNICA Wagner Fernando Delfino Angelotti, André Luiz da Fonseca, Guilherme Barreto Torres e Rogério Custodio* Instituto de Química, Universidade Estadual de Campinas, CP 6154, 13084-971 Campinas – SP, Brasil Recebido em 24/4/07; aceito em 10/8/07; publicado na web em 19/12/07

Educação

Quim. Nova, Vol. 31, No. 2, 433-444, 2008

A SIMPLIFIED APPROACH TO THE QUANTUM MONTE CARLO METHOD: FROM THE SOLUTION OF INTEGRALS TO THE ELECTRONIC DISTRIBUTION PROBLEM. The paper presents an introductory and general discussion on the quantum Monte Carlo methods, some fundamental algorithms, concepts and applicability. In order to introduce the quantum Monte Carlo method, preliminary concepts associated with Monte Carlo techniques are discussed. Keywords: stochastic methods; Quantum Monte Carlo; electronic structure.

INTRODUÇÃO Métodos estocásticos têm sido freqüentemente utilizados para a solução dos mais diversos problemas de natureza microscópica ou macroscópica da matéria, dentre os quais o método de Monte Carlo é usualmente lembrado. Esses métodos promovem as simulações de átomos, moléculas, aglomerados etc. empregando números aleatórios e têm sido frequentemente utilizados no estudo de líquidos, soluções, sistemas amorfos etc.1. A aplicabilidade desses métodos em sistemas tão complexos, embora permita a determinação de importantes propriedades, não é apropriada para a obtenção de grandezas dinâmicas dos sistemas e, nesses casos, é conveniente o uso de métodos de dinâmica molecular, ou seja, resolver as equações de movimento, como as Equações de Newton. Uma importante contribuição no uso do método Monte Carlo ocorreu aproximadamente no início da década de 80 do século passado, quando essa metodologia passou a ser empregada na solução da Equação de Schrödinger para o estudo de propriedades eletrônicas e vibracionais de átomos e moléculas2. Esta aplicação específica, que passou a ser denominada de Monte Carlo Quântico, tem crescido significativamente nos últimos anos, sendo destacada atualmente como um dos métodos mais promissores para resolver a Equação de Schrödinger em termos de nível de precisão3,4. Embora a descrição dos procedimentos gerais para sua utilização seja freqüentemente mencionada na literatura5, verificam-se problemas de natureza formal e de estabilidade computacional que necessitam de uma avaliação mais profunda. A simplicidade de seus algoritmos e a precisão possível de ser alcançada sugere que um investimento nesta direção possibilitará um domínio importante no cálculo de propriedades atômicas e moleculares e fará parte dos métodos de uso convencional. Em termos educacionais, questões fundamentais são levantadas e tornam a noção estatística da interpretação das propriedades eletrônicas e vibracionais da matéria como algo natural. Neste sentido, torna-se desejável que os métodos introduzidos em disciplinas fundamentais de química quântica sejam complementados com a visão e utilização do método Monte Carlo Quântico na solução de problemas tão simples como átomos com poucos elétrons, efeitos de correlação eletrônica etc. ou tão complexos quanto a estrutura de sólidos ou líquidos. *e-mail: [email protected]

Desta forma, estabelecemos como objetivo deste trabalho apresentar uma visão sucinta e atual do método Monte Carlo Quântico, sua aplicabilidade, formas de solução mais comuns e conceitos associados a esta metodologia. A divisão do trabalho será feita na seguinte seqüência: na próxima seção será dada uma breve introdução ao método Monte Carlo, seguida das considerações sobre o formalismo e aplicação do Monte Carlo em sistemas quânticos. Dentre estas considerações, será feita uma abordagem sobre as duas versões mais utilizadas do Monte Carlo Quântico: o Monte Carlo Quântico Variacional6,7 e o Monte Carlo Quântico de Difusão8-13. Discussões sobre a eficiência e limitações do método serão analisadas. Por último, vários algoritmos computacionais serão apresentados, dentre eles alguns exemplos para um melhor entendimento da metodologia aqui exposta. O MÉTODO MONTE CARLO O método de Monte Carlo tem sido assim denominado em homenagem ao caráter aleatório proveniente dos jogos de roleta de Monte Carlo no Principado de Mônaco. Existem registros isolados de sua utilização na segunda metade do século XIX, quando foram realizadas experiências empregando informações casuísticas14. Porém, seu nome e, principalmente, o estabelecimento de um desenvolvimento sistemático do método data de 1944, durante a Segunda Guerra Mundial, época em que foi utilizado como ferramenta de pesquisa para o desenvolvimento da bomba atômica. A simplicidade de seus algoritmos e eficiência na obtenção de resultados em condições extremamente difíceis justifica sua utilização em diversas áreas do conhecimento, como economia, física, química, medicina, entre outras. Atualmente, a denominação “método de Monte Carlo” tornouse uma expressão geral associada ao uso de números aleatórios e estatística de probabilidade. Para que uma simulação de Monte Carlo esteja presente em um estudo basta que este faça uso de números aleatórios na verificação de algum problema15,16. Ao estimar a probabilidade de ocorrência de um evento, pode-se simular um número independente de amostras do evento e computar a proporção de vezes em que o mesmo ocorre. Um exemplo muito simples pode ser dado considerando-se a Figura 1. Na Figura 1a encontramos um quadrado de aresta a e uma elipse. A área do quadrado pode ser facilmente determinada como Área=a2. Embora a expressão para a área da elipse seja trivial, podemos determiná-la utilizando o mé-

434

Angelotti et al.

todo de Monte Carlo como um experimento numérico. Para isso, distribuímos aleatória e homogeneamente um conjunto arbitrário de pontos dentro da área do quadrado (Figura 1b). Vamos considerar que tenham sido distribuídos N pontos. Contamos quantos pontos estão dentro da elipse (ne), quantos estão fora (no) e, conseqüentemente, N= ne+ no. A probabilidade de que um determinado ponto tenha caído dentro da área da elipse é, portanto, ne/N. A área da elipse será determinada como sendo a probabilidade de encontrarmos pontos dentro da elipse multiplicados pela área do quadrado, ou seja, Área(elipse)=a2. ne/N. A estimativa da área da elipse possui um erro em sua determinação que pode ser minimizado com o aumento do número de pontos distribuídos dentro do quadrado. Quando o número de pontos tende a infinito, a área determinada tenderá ao valor exato.

Quim. Nova

número considerável de vezes para minimizar o erro na estimativa. Fundamentos Embora as determinações acima sejam práticas e muito utilizadas, formalmente o método Monte Carlo está freqüentemente associado à solução de integrais. O método é utilizado na solução de equações diferenciais pela conversão destas em equações integrais. A associação do método Monte Carlo com integrais é feita de maneira também muito simples e intuitiva e sua vantagem está na possibilidade de reduzir sistemas com grande número de dimensões através da determinação de uma média. Em cálculo numérico uma integral definida para uma dimensão é escrita como: (1)

Figura 1. (a) Elipse em um quadrado de aresta a e (b) distribuição de pontos aleatórios na superfície do quadrado

Isto é um procedimento de Monte Carlo. Certamente que a área da elipse poderia ter sido determinada de forma muito mais precisa, rápida e direta utilizando-se alternativas. Neste caso, o método de Monte Carlo poderia ser utilizado com o intuito de avaliar sua precisão numérica. Porém, se uma figura muito mais complicada estivesse inserida no quadrado ou se tivéssemos várias elipses, como mostrado na Figura 2, e estivéssemos interessados em determinar não só a área de cada elipse, mas também a área de possíveis sobreposições das mesmas, o método de Monte Carlo tornar-se-ia uma alternativa muito atrativa a ser considerada. Para resolver este problema significativamente mais complicado, aplicaríamos a mesma distribuição de pontos na superfície do quadrado e seriam contados os pontos que caíssem nas regiões de interesse. A área de cada região específica seria a área do quadrado multiplicada pela probabilidade de encontrar pontos naquela região, ou seja, a2.nx/N. Este exemplo permite caracterizar o Monte Carlo como uma técnica simples desde que tenhamos um mecanismo confiável de gerar números aleatórios e possamos repetir o experimento um

Figura 2. Várias elipses distribuídas em uma superfície quadrada com aresta a

A aproximação na Equação 1 pode ser justificada considerando que a coordenada x entre a e b pode ser subdividida em intervalos Δx e em cada intervalo se escolhe um valor de xi que é utilizado para determinar o valor de f(xi). Conhecendo-se f(xi) e Δx pode-se determinar a área do retângulo formado por essas duas quantidades e obter sua área. A soma de todas as áreas finitas corresponde à solução aproximada da integral. Sabemos que quanto menor o valor de Δx, melhor será a estimativa da integral e, conseqüentemente, maior o número de intervalos necessários para dividir o mesmo intervalo entre a e b. O intervalo Δx pode ser determinado considerando-se os limites de integração a e b e o número de divisões desejado k através de Δx= (b-a)/k. Substituindo-se esta forma de definir Δx na Equação 1, temos: (2) Observando-se a Equação 2 podemos identificar a solução da integral como sendo o resultado da determinação de uma média aritmética, ou seja:

(3) – sendo f a média de f(xi) dos pontos amostrados. A qualidade de uma média pode ser estimada por grandezas como a variância, dada por σ2 (IN) =< I2N >-< IN >2. A Equação 3 foi obtida considerandose que Δx é constante e os valores de xi foram definidos para cada intervalo Δx. Podemos questionar a necessidade de “intervalos constantes” ou se podemos determinar o valor da média utilizando uma amostragem de valores “aleatórios” de x. A resposta é evidente e integrais como a Equação 1 ou qualquer outro tipo de integral podem ser resolvidas pelo método de Monte Carlo por meio de uma amostragem uniforme, como vimos no capítulo anterior, de pontos {xi} escolhidos aleatoriamente no intervalo [a,b]. Para sistemas com maior número de dimensões a mesma técnica é aplicada. A única diferença está no fato de que cada coordenada deve ser amostrada dentro dos respectivos limites de integração. Matematicamente podemos escrever a integral multidimensional: (4)

Vol. 31, No. 2

Uma abordagem simplificada do método Monte Carlo quântico

ou: (5) ou ainda em notação mais compacta como: (6) Na Equação 5 cada somatório corresponde a amostragem em cada uma das coordenadas. Para integrais definidas de baixa dimensão técnicas de resolução, como a regra de Simpson, quadratura Gaussiana e outras, são mais eficientes que o método de Monte Carlo. Para sistemas multidimensionais (Equação 4) o esforço computacional cresce rapidamente com o tamanho do sistema, sendo a utilização do método Monte Carlo a mais indicada, uma vez que o erro associado ao cálculo independe da dimensão do sistema5. A implementação da solução das Equações 3 ou 6 através do método Monte Carlo é simplesmente a determinação dos valores da função f(x) com coordenadas obtidas através de números aleatórios. Como exemplo, apresentamos a determinação da integral: (7) Esta integral definida entre os limites 0 e 1 apresenta solução exata igual à π/4, ou seja, 0,7853982. Devemos utilizar um gerador de números aleatórios que produza valores de x entre estes limites. A primeira dúvida nesta etapa certamente será: quantos valores de x devem ser empregados? Não há um número específico. Esta decisão deve ser tomada durante a avaliação da média da função f(x), que neste exemplo é igual a 1/(1+x2). Para realizar uma escolha criteriosa do número de pontos, devemos prever também qual o nível de precisão desejado. Portanto, a maneira prática de obtermos a resposta desejada para as questões acima é gerar progressivamente diferentes valores de x e acompanhar o valor da média acumulando com os novos valores de f(x). A Figura 3 mostra o valor da média acumulada com o aumento do número de pontos para a realização da integração. A linha tracejada corresponde ao valor exato. Verifica-se que à medida que o número de pontos aumenta, o valor da média tende a uma constante que corresponde ao valor da integral e, não por coincidência,

Figura 3. Média acumulada de uma integração numérica utilizando o método Monte Carlo

435

tende ao valor exato. Na prática dificilmente conhecemos o valor exato da integral e, portanto, a solução aceita como correta corresponde àquela em que a média atingiu um valor constante. Neste exemplo em particular, conjuntos de 100 pontos foram acrescentados à média durante 100.000 vezes, o que corresponde a 10.000.000 pontos, produzindo uma média que atingiu o valor igual a 0,78509, ou seja, 0,00031 abaixo do valor exato. Para aumentarmos a precisão deste resultado é necessário aumentar o número de pontos na determinação da média. Este resultado surpreende em dois sentidos. O primeiro em demonstrar que o método de Monte Carlo é factível. O segundo em demonstrar que, neste caso, a obtenção do resultado numérico da integral foi consideravelmente ineficiente. Amostragem por preferência Em casos de integrais complicadas, o método poderia ser uma das poucas alternativas e possivelmente seria utilizado. Entretanto, existem técnicas poderosas que permitem acelerar a convergência do método Monte Carlo para valores precisos com um esforço computacional significativamente inferior. Uma das idéias mais simples está baseada no método denominado de amostragem preferencial (“Importance Sampling”)17. As amostragens anteriores foram feitas com coordenadas distribuídas uniformemente, ou seja, apresentavam a mesma probabilidade. O método de amostragem preferencial considera que serão escolhidos pontos que contribuem mais significativamente para o valor esperado da integral. Consideremos que os pontos selecionados apresentam uma densidade de probabilidade g(x). A Equação 1 pode ser escrita utilizando-se esta distribuição como: (8) Embora o resultado final seja o mesmo, a média será determinada pela equação:

(9) Comparando-se a Equação 9 e a Equação 3 verificamos algumas diferenças importantes. Notamos que a média está sendo determinada não por f apenas, como na Equação 3, mas pela razão f / g. Embora o resultado final seja o mesmo, a determinação dos valores de x está sendo feita de forma diferente. Na Equação 3 consideramos que todos os pontos gerados apresentavam exatamente a mesma probabilidade. Na Equação 9 esta escolha é feita através da função de distribuição g. Esta informação está expressa na Equação 9 pelo subscrito g no canto inferior e à direita da Equação. A Equação 3 pode ser considerada um caso particular da Equação 9. Dentro de um intervalo a e b, se todos os valores de x forem igualmente prováveis, a probabilidade de encontrar um valor qualquer será: 1/(b-a). Portanto, a densidade de probabilidade será: g(x)=1/(b-a). Substituindo-se esta função na Equação 9, obtemos a Equação 3. Os limites de integração estão incorporados na função densidade de probabilidade e não precisamos multiplicar a média de f por (b-a). A semelhança entre as Equações 9 e 3 desaparece se algumas regiões forem mais prováveis que outras. Nestes casos, a função densidade de probabilidade dependerá da posição que está sendo sorteada. A técnica de amostragem preferencial impõe algumas restrições de uso. Uma delas é que a densidade de probabilidade g(x)

436

Angelotti et al.

seja normalizada, ou seja: (10) A normalização garante que a função densidade de probabilidade proporcionará uma probabilidade finita em todo intervalo de interesse. Outra exigência é que g(x) seja tão próxima de f(x) quanto possível, mas que seja fácil de calcular, ou seja, é desejável que a razão entre f(x)/g(x) seja aproximadamente constante. Finalmente, uma vez que a função g(x) representa uma densidade de probabilidade, ela deve ser sempre positiva. Uma importante conseqüência da utilização deste método é que a variância será minimizada em relação à mesma determinação obtida através da distribuição uniforme. A obtenção de variáveis em uma distribuição específica é utilizada para distribuições usuais, como a distribuição gaussiana, exponencial etc.18. Porém, uma das formas mais populares de mapear regiões significativas a partir de uma distribuição, com um nível de simplificação considerável, é o algoritmo de Metropolis19. O algoritmo de Metropolis O método de Metropolis19 é considerado um dos dez algoritmos mais importantes do século passado e tem possibilitado a realização de simulações baseando-se na técnica de aceitar ou rejeitar alterações na modificação dos pontos distribuídos no espaço. O mapeamento das variáveis nos exemplos anteriores foi feito considerando-se apenas que novos pontos “sorteados” deveriam ser acrescentados através de uma distribuição específica. No algoritmo de Metropolis esta noção é alterada significativamente. Usualmente nos deparamos com denominações decorrentes do desenvolvimento formal do algoritmo de Metropolis. Por exemplo, a seqüência de alterações nas configurações ou coordenadas é controlada por um processo de Markov 20 . Em um processo Markoviano as condições finais de um processo não dependem das condições prévias. Uma cadeia de Markov pode ser especificada pela escolha de valores de probabilidades de transição, ou seja, a chance de um conjunto ou uma configuração a ser alterada em b, P(ra → rb). Podemos considerar cada cadeia de Markov no sistema movendo uma partícula fictícia, chamada “walker” (ou “psip”), em torno das diferentes configurações. Todos os “walkers” movem-se passo a passo de acordo com as probabilidades de transição P(ra → rb). Metropolis percebeu que a distribuição de “walkers” ficaria abaixo da distribuição exigida, g(r), contanto que as probabilidades de transição obedecessem à Equação: (11) Esta igualdade é denominada de balanceamento detalhado e é freqüentemente invocada no método Monte Carlo. Demonstrações rigorosas podem ser encontradas na literatura21 sobre a formalização do método de Metropolis. Computacionalmente podemos sistematizar o método de Metropolis para a solução de integrais utilizando amostragem preferencial através do seguinte algoritmo: 1) Criar um conjunto de pontos ou configurações {r1,r2,...,rM} arbitrários e calcular o valor de g(r1,r2,...,rM). Caso g(r) seja uma função de distribuição de uma dimensão, calcula-se g(r1), g(r2), ..., g(rM). 2) Alterar o valor das coordenadas empregando um mecanismo qualquer gerando {r’1,r’2,...,r’M} e calcular o(s) valor(es) de g para as coordenadas modificadas. 3) Efetuar a razão entre a densidade de probabilidade com as coordenadas novas por aquela obtida com as coordenadas

Quim. Nova

antigas, g(r’1,r’2,...,r’M)/g(r1,r2,...,rM)=P. Em casos de uma dimensão, calcular g(r’1)/g(r1)=Pr1→r’1, g(r’2)/g(r2)=Pr2→r’2, ..., g(r’M)/ g(rM)=PrM→r’M. 4) Comparar a probabilidade de transição P (individual ou coletiva) com um número aleatório com distribuição uniforme entre 0 e 1. Se P for maior que o número aleatório, então, as novas coordenadas devem aceitas e substituem as antigas. Se a razão P for menor que o número aleatório, as novas coordenadas devem ser descartadas e as coordenadas antigas devem ser preservadas. 5) Após realizar o teste, calcular o valor de f(r1,r2,...,rM)/ g(r1,r2,...,rM) ou, em caso de uma dimensão, f(r1)/g(r1), f(r2)/g(r2),..., f(rM)/g(rM) para determinar o valor médio e outras grandezas estatísticas. 6) Se a média não tiver alcançado um valor constante, voltar à etapa 2 e repetir o processo. Para ilustrar a aplicação do algoritmo acima, vamos aplicá-lo na solução da integral definida anteriormente pela Equação 7 utilizando uma função de distribuição que corresponde à reta de mínimos quadrados mostrada na Figura 4. A reta de mínimos quadrados normalizada, que será nossa função de distribuição, para esta função é igual a: g(x)=2x(1,051127-0,53865.x)/(2x1,051127-0,53865). Com o algoritmo de Metropolis, a média definida pela Equação 9 e o mesmo número de pontos utilizando anteriormente (100 pontos) modificados durante 10.000 passos, produziu uma média igual a 0,78542, ou seja, um valor 0,00009 acima do valor exato. Além da melhora significativa no resultado, notamos uma melhora estatística significativa no que diz respeito ao desvio padrão em ambos os cálculos. Sem amostragem preferencial o desvio padrão atingiu um valor igual a ±0,161, enquanto que com amostragem preferencial o resultado é cerca de 10 vezes inferior, ou seja, ±0,018.

Figura 4. Gráfico da função 1/(1 + x2) vs x (linha sólida) e reta nãonormalizada obtida por mínimos quadrados (linha tracejada)

Um detalhe importante não comentado anteriormente no uso do algoritmo de Metropolis diz respeito à taxa de aceitação. Esta quantidade é definida como o número de modificações aceitas dividida pelo número de total de modificações. Na alteração das coordenadas devemos incluir um mecanismo que controle os movimentos aceitos. Na literatura considera-se que a taxa de aceitação para obtermos valores adequados através do Monte Carlo é de 0,5, ou seja, 50% das modificações das coordenadas devem ser aceitas em média. No exemplo acima, utilizamos a seguinte expressão para alterar as coordenadas: r’a = rb + (0.5 – rand).δ

(12)

em que rand é um número aleatório entre 0 e 1 obtido de uma

Vol. 31, No. 2

Uma abordagem simplificada do método Monte Carlo quântico

distribuição uniforme e δ é um número que deve ser ajustado para produzir aceitação em torno de 50% quando submetida ao algoritmo de Metropolis. Usualmente testes iniciais são realizados com a finalidade de ajustar-se este parâmetro antes de acumular os valores da média. O método Monte Carlo Quântico O método de Monte Carlo recebe a denominação particular de Monte Carlo Quântico (MCQ) quando empregado na resolução da Equação de Schrödinger. Existem diferentes alternativas para realizar cálculos MCQ. As duas mais utilizadas, e que serão exploradas neste trabalho são: Monte Carlo Quântico Variacional (MCQV)22-26 e Monte Carlo Quântico de Difusão (MCQD)27-30. Embora o objetivo principal seja resolver a Equação de Schrödinger, que é uma equação diferencial, devemos lembrar que o Monte Carlo é uma metodologia empregada para a resolução de integrais. Logo, precisaremos converter a Equação de Schrödinger em uma equação integral. O método Monte Carlo Quântico Variacional

437

(17) Podemos amostrar o espaço de configurações representado pela função de distribuição |ΨT (q)|2 utilizando o algoritmo de Metropolis. Realizamos modificação das coordenadas q → q’ no espaço de configurações de 3N dimensões movendo-se uma ou mais partículas. A partícula i pode ser movida da posição qi para a posição qi’= qi+δqi, em que δqi deverá ser maior ou menor que zero e escolhida a partir de uma distribuição gaussiana. A taxa de aceitação do movimento é determinada através da Equação: (18) Fazendo uso do teorema do valor médio, a média ponderada da energia torna-se: (19)

O MCQV é o mais simples, dentre os métodos MCQ, para o cálculo de propriedades atômicas ou moleculares. O MCQV está baseado no método variacional31 para determinar energias, embora possa ser aplicado para qualquer propriedade definida por um operador quântico. Para encontrarmos o valor esperado da energia de um estado representado por uma função de onda de N partículas, minimizamos a integral do valor médio de 3N dimensões: (13) sendo q o vetor das 3N coordenadas eletrônicas e de N coordenadas de spin, isto é, q={r1ξ1, r2ξ2,..., rNξN}, em que r1={x1,y1,z1}, ..., rN ={xN,yN,zN}. O elevado custo computacional para se calcular uma integral multidimensional através de um método de quadratura convencional e a precisão diretamente relacionada ao número de dimensões do sistema recomenda o uso da aproximação de Monte Carlo. O MCQV usa diretamente o algoritmo de Metropolis na Equação 13. Uma função de onda teste ΨT(q) é usada como uma aproximação para a função de onda de muitos corpos verdadeira. Assim, um limite superior para a energia E0 do estado fundamental é obtido:

sendo M o número de pontos utilizados. O subscrito |Ψ(q)|2 indica que esta média é obtida através da amostragem da população de pontos distribuídos segundo a função: (20) A função ρ(q) é a densidade de probabilidade. Desta forma, o cálculo do valor esperado da energia é transformado de uma integral para uma média aritmética ponderada. O erro nesse cálculo é dado pelo desvio-padrão calculado convencionalmente como:

(21) Se considerarmos o caso em que a função de onda teste Ψ(q) é ^ autofunção de H , a energia local EL(q) é uma constante em qualquer ponto do espaço de configurações, isto é, (22)

(14) ^

Dividindo e multiplicando o lado esquerdo do operador H por ΨT obtemos:

(15) ^

sendo EL (q) = H ΨT (q) / ΨT (q) definido como a energia local do sistema. Utilizando o conceito de amostragem preferencial, podemos relacionar:

Como conseqüência, o valor médio da energia apresentará um desvio padrão nulo, indicando que quanto mais próxima Ψ(q) estiver da função de onda exata, mais precisa será a estimativa para o estado desejado. Portanto, costuma-se afirmar que o erro associado aos resultados é causado unicamente pela imprecisão da função de onda, sendo o desvio-padrão aplicado como um avaliador da qualidade da função de onda. O átomo de hélio será usado como exemplo do método MCQV. Uma função de onda tentativa incluindo efeito de correlação eletrônica com um parâmetro variacional (ζ) pode ser escrita como: (23)

(16) com a Equação 15 reescrevendo-a como:

sendo r1 e r2 as distâncias espaciais dos elétrons 1 e 2 em relação ao núcleo, r12 é a distância intereletrônica e ξ1 e ξ2 as coordenadas de spin. O operador Hamiltoniano em unidades atômicas para o átomo de He pode ser definido como:

438

Angelotti et al.

Quim. Nova

(24) Os dois primeiros termos à direita da Equação 24 são operadores de energia cinética dos elétrons 1 e 2, seguidos dos operadores de atração nuclear e o operador de repulsão intereletrônica. Utilizando-se a definição de energia local: (25) obtemos a expressão para a energia local dada por: (26) em que os termos r^12, r^1 e r^2 são vetores unitários para as respectivas distâncias espaciais. Devemos chamar a atenção para o fato de que a função de onda escolhida (Equação 23) foi fatorada em termos de uma parte espacial e uma parte de spin e o operador Hamiltoniano não depende das coordenadas de spin. Neste caso, a parte de spin será cancelada na Equação 25 e a energia local será um número real independente das funções de spin. Para funções de onda mais complexas que envolvam determinantes de Slater, alternativas numéricas são adotadas32 e serão discutidas adiante. Com a função de onda definida pela Equação 23 utilizamos o algoritmo de Metropolis:

(27) Com a Equação 27 definindo as modificações das coordenadas espaciais e com a Equação 26 para determinar a energia local, a energia média do sistema será definida pela Equação 19. Outras propriedades podem ser determinadas da mesma maneira. As alterações das coordenadas espaciais dos dois elétrons podem ser feitas ao mesmo tempo ou um elétron de cada vez. Usualmente, a segunda opção é utilizada, visto que a chance de alterar a configuração eletrônica é maior em relação às alterações simultâneas de todos os elétrons. Uma alternativa também utilizada no Monte Carlo Quântico Variacional é a de considerar não apenas a energia local de uma configuração de cada vez, mas trabalhar com um conjunto de configurações, modificar as coordenadas dos elétrons em todas as configurações e determinar a energia média final como a média das médias, ou seja:

(28) A Figura 5 corresponde ao gráfico de simulações envolvendo a modificação das coordenadas de 500 “walkers” ( N’ = 500) por 10000 passos (M’ = 10000) com diferentes valores de ζ. As coordenadas de cada elétron foram modificadas uma de cada vez e a energia local de cada configuração só foi calculada após aplicação do critério de Metropolis. A energia média com menor valor foi obtida para ζ=0,63 e apresentou um valor igual a -2,8665±0,00012 kcal/mol. Diferente dos métodos quânticos convencionais, o método MCQV está sujeito a uma incerteza que é calculada durante a si-

Figura 5. Comportamento da energia local acumulada para diferentes simulações para o átomo de He no estado fundamental. Cada simulação foi realizada com 500 walkers e 10000 passos, utilizando função de onda incluindo efeito de correlação explícita com diferentes parâmetros ζ

mulação. O desvio-padrão é obtido da amostragem realizada e do número de passos utilizado. Para que o desvio-padrão tenda a zero, devemos aumentar o número de passos e quando o número de passos tende a infinito, a determinação da integral tende ao valor exato e o desvio-padrão tende a zero. A grande vantagem do método Monte Carlo Variacional está no fato de que para se determinar o valor da energia ou qualquer outra propriedade empregando o teorema do valor médio basta utilizar o algoritmo de Metropolis e calcular o valor da propriedade local em um número grande de configurações. A desvantagem está na necessidade de termos que definir a função de onda tentativa e testar diferentes alternativas que nos levem à representação de menor energia possível. Para evitar o desenvolvimento de inúmeras funções de onda e o tedioso trabalho de testá-las, podemos lançar mão de uma alternativa mais poderosa, o Monte Carlo Quântico de Difusão. O método Monte Carlo Quântico de Difusão O MCQD é, provavelmente, o método de Monte Carlo Quântico mais utilizado em cálculos de estrutura eletrônica. Esse método, em princípio, produz resultados exatos, mas exige uma série de aproximações para que sua utilização prática seja possível. O MCQD consiste na solução da Equação de Schrödinger através do movimento browniano de partículas no espaço sujeitas a um processo probabilístico de criação e destruição dessas mesmas partículas. Intuitivamente, este processo foi idealizado a partir da associação fenomenológica da Equação de Schrödinger com uma Equação de difusão. A Equação de Schrödinger dependente do tempo para um conjunto de N partículas com coordenadas q é escrita como: (29) Comparando-se esta equação com uma equação de difusão do tipo: (30) em que C corresponde à concentração de um material sendo obser-

Vol. 31, No. 2

Uma abordagem simplificada do método Monte Carlo quântico

vado durante o processo de difusão, nota-se uma semelhança considerável entre as duas expressões, sugerindo que, em princípio, a Equação de Schrödinger possa ser resolvida com os mesmos recursos numéricos que a Equação de difusão. Embora a Equação de Schrödinger apresente um componente imaginário em sua formulação, podemos mascarar esta diferença considerando a transformação do tempo real (t) para tempo imaginário (τ), ou seja: τ = i.t. Vale ressaltar também que os estados que se pretende obter por meio da Equação 29 são estados estacionários, ou seja, aqueles para os quais ∂Ψ(q,t)/∂t = 0. Neste caso, introduzimos um parâmetro arbitrário ET para deslocar a energia de referência a um valor arbitrário. Assim, realizando a transformação da coordenada temporal, introduzindo a energia de referência arbitrária ET e considerando que a Equação de Schrödinger está sendo empregada para o tratamento de elétrons (férmions), a Equação 29 passa a ser escrita, em unidades atômicas, como: (31) Para resolver a Equação de Schrödinger como uma Equação de difusão, ela é separada e avaliada sob dois aspectos distintos. O primeiro deles considera a Equação 30 sem o segundo termo à direita da igualdade:

439

durante a simulação, corresponde à energia do sistema, assim como, o valor de outras propriedades também será definido pelo valor médio da respectiva propriedade. Embora seja possível realizarmos simulações empregando estas duas estratégias, dinâmica browniana e criação e aniquilação de partículas definidas por uma cinética de primeira ordem, estas simulações são muito sensíveis às condições iniciais e o controle do processo necessita de ajustes. Funções potenciais coulombicas tendem a mais ou menos infinito quando as distâncias são iguais a zero. Experimentos computacionais empregando alternativas numéricas foram testados até que se definiu um formalismo matemático rigoroso justificando o método. O método MCQD encontrou sua fundamentação teórica no tratamento de equações diferenciais através do uso de funções de Green33,34. Formalmente, devemos levar em consideração que o método de Monte Carlo não pode ser aplicado diretamente às equações diferenciais. Portanto, antes de estabelecer uma simulação computacional para o processo de difusão, procura-se uma representação integral da Equação de Schrödinger. Utilizando-se o inverso do ^ Hamiltoniano, H -1, e aplicando-se do lado direito e esquerdo da Equação de Schrödinger independente do tempo, tem-se: (36) ou simplesmente, (37)

(32) ^

que corresponde à representação matemática do movimento browniano. A solução para esta equação diferencial corresponde a uma função de distribuição gaussiana do tipo:

Como o hamiltoniano é um operador diferencial, H -1 deve ser um operador integral e: (38)

(33) em que δτ representa um intervalo de tempo e D a constante de difusão. Para simular movimentos brownianos basta distribuir um conjunto de partículas com coordenadas arbitrárias no espaço e deslocar cada partícula aleatoriamente, obedecendo a uma distribuição gaussiana, como a Equação 33, com média zero e variância 2Dδt. A segunda parte da equação sem o termo browniano corresponde a uma equação diferencial ordinária de primeira ordem. Em química, sua associação está vinculada a uma cinética de reação de primeira ordem, ou seja: (34)

^

em que a integral do operador G (q’,q) corresponde a H -1. Para processos dependentes do tempo, como é o caso da Equação 31, a Equação 38 passa a ser escrita como: (39) Esta Equação determina que em um tempo τ a chance de alterar-se uma configuração q para uma configuração q’ se dá através das características da função G (q’,q; τ). Não é possível definir uma forma analítica para a função G (q’,q; τ). Para resolver a Equação 39 utilizamos uma aproximação válida para o limite em que τ tende a zero. Neste caso, a função de Green pode ser fatorada em um termo de difusão (Gdif) e um termo cinético (GB), ou seja:

Esta equação possui uma solução analítica bem conhecida: (40) (35) Em um processo cinético algumas substâncias são criadas e destruídas e para realizarmos uma simulação de tal processo, devemos associar a Equação de Schrödinger (Equação 31) ao processo de criação e aniquilação de partículas. A probabilidade que controla esse processo pode ser definida em termos da exponencial na Equação 35 e envolve, no caso da Equação de Schrödinger, a função potencial (V(q)-ET). Combinando-se a simulação do movimento browniano com a cinética de primeira ordem podemos estabelecer uma situação de equilíbrio entre os dois processos. Nesta situação de equilíbrio a energia média das configurações que aparecem e são preservadas

Estes dois termos estão diretamente associados com as Equações 33 e 35, respectivamente, o que nos permite definir o termo de difusão por: (41) e o termo cinético por: (42)

440

Angelotti et al.

Na Equação 42 a energia ET foi incorporada e corresponde à energia de referência. As expressões apresentadas para Gdif e GB permitem desenvolver um algoritmo para o Monte Carlo de Difusão capaz de calcular os valores médios de propriedades do estado fundamental. Contudo, esse algoritmo torna-se ineficiente em muitos casos, sobretudo pela divergência do termo {V (q) + V (q’)}/ 2 nas regiões em que há singularidade do potencial. Nessas condições, há variação significativa do número de partículas do sistema, o que confere grande incerteza estatística nos valores médios calculados. A forma de reduzir essas flutuações é por meio da técnica de amostragem preferencial, já vista acima, baseada no formalismo de Fokker-Planck5,35. Assim, combinando-se os melhores aspectos do formalismo de Fokker-Planck com o algoritmo de Metropolis desenvolve-se um método de Monte Carlo Quântico relativamente estável. Nesse caso, utiliza-se uma função de onda tentativa Φ como guia para a amostragem estatisticamente mais importante. Assim, define-se uma função densidade de probabilidade f (q, τ), tal que:

Quim. Nova

do sistema E0, eliminando o termo correspondente ao potencial. Desta forma, o número de partículas do ensemble permanece constante. Para uma função aproximada, não há o cancelamento entre as singularidades da energia cinética e potencial, acarretando em flutuação do número de partículas do ensemble para compensar as imperfeições da função tentativa. As partículas são eliminadas nas regiões em que Φ (q) > Ψ0 (q) e são criadas nas regiões que Φ (q) < Ψ0 (q). Assim, a função Φ (q) tende a se aproximar da função de onda exata à medida que a simulação prossegue. O termo de energia cinética é identificado como: (48) Considerando que a força F(q) permaneça constante durante todo o movimento, a correspondente função de Green para a parte cinética é: (49)

(43) sendo Ψ (q, τ) a função de onda dependente do tempo para o sistema. No decorrer da simulação, essa função deve se aproximar da função de onda exata do estado fundamental. Substituindo-se essa função densidade f (q, τ) na Equação de Schrödinger dependente do tempo (Equação 31) tem-se:

Essa Equação é apenas aproximada e não corresponde exatamente à solução da Equação de Fokker-Planck. O erro associado ao uso dessa Equação é proporcional a τ2, logo, quanto menor o valor de τ, menor o erro devido ao emprego dessa função aproximada. O deslocamento compatível com a Equação 49 é: (50)

(44) sendo EL a energia local e F um vetor de campo denominado força quântica (ou vetor velocidade) de Fokker-Planck. Essas funções são definidas pelas equações abaixo: (45)

(46) A inclusão da amostragem preferencial modifica as funções de Green utilizadas na estratégia de Monte Carlo. A Equação 44 guarda semelhanças com a Equação 31, o que permite sua resolução pela mesma metodologia. A parte cinética é simulada pelo processo de difusão e o termo correspondente ao potencial por meio da criação e destruição de partículas. Identifica-se o termo cinético da Equação 44 como sendo os dois primeiros à direita do sinal de igualdade e o termo correspondente ao potencial como o último termo à direita. A função de Green para a parte potencial (Equação 43) deve ser adaptada, substituindo-se simplesmente V(q) e V(q’) por EL. A função de Green para esta nova representação assume a forma: (47) O processo de criação e destruição de partículas é melhor definido pela Equação 47 que pela função de Green sem amostragem preferencial, pois a energia local é mais suave, minimizando as possíveis singularidades apresentadas pela respectiva função potencial. Quando Φ (q) = Ψ0 (q), sendo Ψ0 (q) a função de onda exata do sistema, a energia de referência ET torna-se igual à energia exata

sendo η um número aleatório de uma distribuição gaussiana com média zero e variância igual a um. Uma vez que o operador de energia cinética não é hermitiano, faz-se necessário utilizar o balanceamento detalhado, associado às Equações 47 e 49. Se utilizássemos diretamente o algoritmo de Metropolis, deveríamos utilizar apenas a razão entre as densidades de probabilidades da função de onda guia. Porém, como o operador de Green não é hermitiano, é necessário introduzir-se uma correção que permita sugerir que os caminhos de ida e volta de uma partícula ou configuração sejam equivalentes. Assim, a probabilidade de transição do algoritmo de Metropolis é modificada de acordo com a equação: (51) sendo (52) Esta estratégia utiliza uma função de Green como combinação de duas outras: a função de Green da Equação de Fokker-Planck (Equação 49) e a função de Green do processo de criação e destruição de partículas (Equação 47). Com todas estas adaptações, determina-se a energia média do sistema E através de uma média da energia local EL (Equação 45):

(53) No MCQD são obtidas configurações (N) que correspondem à distribuição f por meio do algoritmo de Metropolis e a energia é

Vol. 31, No. 2

Uma abordagem simplificada do método Monte Carlo quântico

calculada como uma média da energia local: (54) Devemos lembrar que o número de configurações muda durante a simulação e a média deve considerar essa mudança de população, tornando-se uma média ponderada. O algoritmo do método MCQD pode ser encontrado na referência 10. Como exemplo de aplicação deste algoritmo consideramos a molécula de H2 na geometria de equilíbrio (R=1,4011 u.a.) e estado fundamental sem função de correlação explícita. Neste cálculo o orbital 1σg da molécula foi obtido através de cálculo HartreeFock empregando-se em cada átomo combinação linear de funções de Slater 1s triplo zeta e uma função de polarização 2p orientada ao longo do eixo da ligação. A composição do orbital molecular foi otimizada no ambiente e corresponde a: 1σ g = [-0,02605.1s a(0,65982) + 0,24506.1s a(0,93256) + 0,32054.1s a (1,33002) – 0,02805.2pz a(1,82453) – 0,02605.1s b(0,65982) + 0,24506.1s b(0,93256) + 0,32054.1s b(1,33002) + 0,02805.2pz b (1,82453)]. Os valores entre parênteses correspondem aos expoentes das respectivas funções de Slater e também foram otimizados na molécula em nível Hartree-Fock. O sub-índice a ou b indica se as funções foram colocadas sobre o hidrogênio a ou b. Foram realizadas várias simulações com 1.000 configurações cada por 2.000.000 de passos e diferentes valores de τ. A Figura 6 apresenta os valores das diferentes energias médias obtidas em cada simulação. O valor da energia média extrapolado para τ=0 foi igual a -1,17440±2x10-5 u.a. O valor considerado exato para esta molécula no estado fundamental é igual a -1,17447 u.a.5 A inclusão de função de correlação explícita neste caso permitiria apenas melhorar a precisão numérica em termos da quinta ou sexta casa decimal. Este exemplo mostra ainda que uma boa função de onda guia, mesmo construída em nível Hartree-Fock, produz resultado com parte significativa de efeitos de correlação eletrônica.

441

pre um resultado melhor que o Variacional, uma vez que a função de onda serve apenas como um guia para a simulação. Para ser utilizada de forma eficiente, principalmente no MCQD, a função de onda aproximada deve considerar alguns aspectos gerais observados em funções de onda exatas36. São conhecidas várias propriedades das soluções analíticas da Equação de Schrödinger para sistemas simples e que devem ser obedecidas, para que a função de onda se aproxime da função exata, por exemplo: 1) a função deve ser contínua, unívoca, finita e a densidade eletrônica deve apresentar um comportamento assintótico apropriado com o aumento das distâncias nucleares da molécula e 2) a função de onda exata deve ser anti-simétrica com a permutação de qualquer par de elétrons. Além dessas propriedades globais, existem outras de caráter local. A energia local deve ter um valor finito para qualquer conjunto de posições espaciais dos elétrons. Portanto, as singularidades que surgem no potencial pela sobreposição de cargas devem ser compensadas pelo termo correspondente da energia cinética local. Essa compensação resulta numa descontinuidade na derivada primeira da função de onda quando duas partículas estão sobrepostas, conhecida como cúspide37. Para que a condição de cúspide seja suprimida, no caso da possibilidade de colisão núcleo-elétron, verifica-se que uma função de onda bem comportada satisfaz a relação: (55) sendo ri a distância do elétron i ao núcleo e Z é a carga nuclear. Se a função de onda não é anulada quando ri=0, a condição de cúspide nuclear é satisfeita se essa função exibe um comportamento exponencial em ri. As funções de base de Slater apresentam esta característica. Exemplificando para o átomo de He, a condição de cúspide nuclear é satisfeita multiplicando-se uma função de onda arbitrária por um termo exponencial: (56) sendo que a1 e b1 deverão ser iguais a –Z, enquanto que a2 e b2 são parâmetros que podem ser ajustados variacionalmente. Para sistemas moleculares utiliza-se o mesmo procedimento. O comportamento eletrônico em moléculas é determinado através de funções de onda moleculares obtidas a partir de funções monoeletrônicas, que são usualmente expansões em funções de base centradas nos núcleos atômicos. Assumindo que o comportamento dos elétrons próximos dos núcleos é completamente determinado pela forma analítica da função de base, a condição de cúspide nuclear impõe a mesma restrição à forma dessas funções de base nas regiões próximas dos núcleos5,38-40. Uma extensão para a condição de cúspide elétron-elétron também pode ser feita através das condições:

Figura 6. Convergência da energia média em função de τ para simulação da molécula de H2 na geometria de equilíbrio utilizando MCQD

A função de onda tentativa Embora a apresentação do Monte Carlo Quântico feita acima pareça relativamente simples, tanto em sua formulação Variacional quanto na Difusão, devemos ter em mente que a escolha da função de onda representa um passo importante em qualquer das abordagens. No MCQV serão obtidos resultados correspondentes aos da função de onda escolhida e MCQD a simulação proporcionará sem-

(57)

Dessa forma, a condição de cúspide eletrônico exige que a função de onda seja contínua, mas que tenha a derivada primeira em relação as coordenadas eletrônicas descontínuas quando os elétrons estiverem sobrepostos. Caso a função de onda e suas respectivas funções de base não satisfaçam as condições de cúspide eletrôni-

442

Angelotti et al.

co, devemos multiplicar a função de onda por um termo exponencial semelhante ao da Equação 56, mas satisfazendo as condições de cúspide da Equação 57. As funções que devem ser incorporadas à função de onda aproximada não são únicas: existem várias alternativas na literatura38-40. Além das questões relacionadas ao cúspide eletrônico, que estão intimamente relacionadas com a correlação eletrônica, o Monte Carlo Quântico considera uma aproximação para a função de onda que merece atenção. As expressões nas definições da densidade de probabilidade (Equações 16, 18, 52) ou propriedades locais (Equações 22, 45, 46) só podem ser utilizadas se a função de onda empregada permitir o cancelamento das funções de spin entre o numerador e o denominador. A função de onda típica em simulações MCQ que permite o cancelamento das coordenadas de spin é descrita como um produto entre determinantes de Slater de funções de spins opostos2,32,41, ou seja:

Quim. Nova

utilizadas funções monoeletrônicas que individualmente obedeçam a essa mesma exigência. As funções de base de Slater apresentam este comportamento e, por isso, são as mais empregadas no MCQ. Ao contrário, as funções gaussianas não apresentam um perfil adequado nas regiões próximas dos núcleos. Por outro lado, a condição de cúspide eletrônica (Equação 57) não pode ser satisfeita por funções de onda compostas por um único determinante porque não apresentam dependência com a distância intereletrônica (rij). Usualmente, multiplica-se a função de onda da Equação 58 por funções de correlação explícitas42, ou seja: Ψcorr = e-U

(60)

sendo que U pode ser uma função de correlação do tipo BoysHandy43: (61) em que I são índices nucleares e i e j índices eletrônicos; ou uma função do tipo Pade-Jastrow, por exemplo, com 2 parâmetros para átomos ou moléculas (62)

(58) sendo φi a parte espacial da função de onda monoeletrônica. Usualmente, na função Ψ(q) = Ψα (qα) Ψβ (qβ), incorpora-se uma função de correlação explícita Ψcorr, ou seja, Ψ(q) = Ψα (qα) Ψβ (qβ)Ψcorr, que será descrita adiante. Por enquanto, descreveremos apenas a função de onda dada pela Equação 58. Essa função de onda permite fatorar as funções de spin das funções espaciais, possibilitando o cancelamento das funções de spin no cálculo das propriedades locais; é finita, normalizável e permite o cálculo de propriedades locais sem grandes dificuldades5. Essa última característica pode ser verificada observando-se a Equação 46 que requer a determinação da derivada da função de onda com relação às coordenadas espaciais dos elétrons. Essa derivada será simplesmente a soma dos determinantes, cujas colunas foram seqüencialmente substituídas por suas respectivas derivadas. Para a derivada da função de onda em relação a um único elétron α com coordenadas espaciais r1, por exemplo, tem-se que:

(59)

Além dessa característica, a função de onda fatorada pode satisfazer a condição de cúspide nuclear se na sua construção forem

Os parâmetros CkI da Equação 61 e a e b da Equação 62 são ajustados para satisfazerem as condições de cúspide e minimizar a energia do sistema. A inclusão dessas funções de correlação explícita corrige os problemas de cúspide nuclear e eletrônico, aumentando a precisão dos resultados e tornam as simulações mais estáveis. Vantagens e desvantagens do uso do MCQ Algumas características importantes podem ser comentadas em relação a este método, tanto para o MCQV quanto para o MCQD: 1) o problema do conjunto de base não é muito discutido no MCQD uma vez que o uso de funções de base de Slater permite a inclusão de conjuntos com qualidade superior em relação às funções gaussianas, com um número de primitivas significativamente menor. Devemos considerar que erros devido ao uso de um conjunto de bases finito são bem menores, já que a função de onda multieletrônica é utilizada para expandir a função de onda guia exigida para a amostragem preferencial. 2) Seus algoritmos são mais simples em relação a métodos convencionais de estrutura eletrônica alcançando excelente precisão na determinação de propriedades. 3) Os algoritmos para o MCQ são facilmente adaptados para computadores paralelos e escalonam linearmente com o número de processadores. Não há congestionamento na memória do computador ou necessidade de grande capacidade de disco para sistemas grandes. 4) Funções de onda multieletrônicas com dependência explícita das distâncias entre as partículas (rij) podem ser utilizadas sem dificuldade. 5) Propriedades de estados fundamentais, alguns estados excitados, barreiras de reação química e outras propriedades podem ser determinadas dentro de um sistema unificado. 6) O uso de pseudopotenciais em um cálculo produz escalamento computacional da ordem de Z3,5, sendo Z a carga nuclear. 7) Boa convergência para a maioria dos sistemas estudados. O MCQ também sofre com algumas limitações que fazem com que seu uso seja considerado. Dentre as desvantagens podemos citar: 1) demanda um custo computacional razoável para calcular derivadas de primeira e segunda ordem da energia local total com

Vol. 31, No. 2

Uma abordagem simplificada do método Monte Carlo quântico

respeito às posições atômicas e, conseqüentemente, estimar forças interatômicas e constantes de força, ou seja, inabilidade na otimização de geometria e simulação dinâmica. 2) Os resultados são obtidos com uma barra de erro estatística que decai somente com o inverso da raiz quadrada do tempo computacional, ou seja, do número de passos da simulação e número de configurações consideradas. 3) Para sólidos, a necessidade de cálculos de propriedades do “bulk” produz erros sistemáticos de tamanho finito. 4) Apesar de uma sensível melhora no cálculo de propriedades de estados excitados, estas informações ainda são raras quando comparadas às propriedades calculadas por outros métodos correla-cionados. 5) A precisão das simulações é dependente das posições dos nós da função de onda tentativa e, durante as simulações, a possibilidade de uma partícula assumir as coordenadas de um nó da função de onda produzirá efeitos indesejáveis na simulação e devem ser controlados. Dos artigos de revisão divulgados na literatura destacamos um, de M. D. Towler44. Este artigo propõe doze questões referentes aos problemas do MCQ e é um excelente referencial para se avaliar se o método atenderá, ou não, às expectativas dos usuários. O custo computacional do MCQ pode ser comparado com outros métodos através da Tabela 1. Os métodos de estrutura eletrônica convencionais são muito precisos para cálculos de estrutura eletrônica, porém, funcionam bem somente para sistemas contendo poucos elétrons. Para sistemas maiores, o MCQ torna-se muito competitivo em termos de precisão e tempo computacional para cálculos de estrutura eletrônica. Na Tabela 1 podemos verificar que os cálculos MCQ realizados com todos os elétrons para sistemas grandes se tornam muito caros computacionalmente, principalmente quando comparado com a TFD. Quando o cálculo MCQ é feito com elétrons de valência utilizando pseudopotenciais, este tempo computacional cai quase que pela metade e, devido a sua alta precisão, pode-se tornar uma alternativa mais eficaz e precisa que a TFD e, principalmente, em relação às demais metodologias apresentadas na Tabela 1. Tabela 1. Comparação do escalonamento computacional aproximado para alguns métodos de cálculos de estrutura eletrônica Método HF DFT (LDA) MP2 MP3 MP4 CISD CCSD MCQa MCQb

Escalamento N3 N3 N5 N6 N7 N6 N6 5,5 N – N6,5 N3,5

a

O método MCQ considerando todos os elétrons. b O método MCQ considerando pseudopotenciais e apenas elétrons de valência. COMENTÁRIOS FINAIS Disciplinas iniciais de mecânica quântica na área de química seguem usualmente a sistemática de introduzir os postulados, as soluções exatas da Equação de Schrödinger para sistemas simples e, eventualmente, a apresentação do método Hartree-Fock. Nas disciplinas mais avançadas são discutidas questões associadas a modelos utilizados para a correção dos efeitos de correlação eletrônica. A impressão final em qualquer uma das fases é de que são necessários conhecimentos matemáticos e computacionais extrema-

443

mente sofisticados para se atingir um nível de precisão aceitável mesmo para sistemas simples. Neste sentido, o método MCQ pode ser considerado como uma alternativa aos métodos convencionais com aplicações imediatas em qualquer nível de dificuldade desejado, permitindo a obtenção de resultados precisos a partir de funções de onda simples. Efeitos de correlação eletrônica, por exemplo, passam a ser um termo possível de ser determinado com um mínimo de esforço computacional e novas perspectivas de conceitos físicos e matemáticos são abertas ao iniciante, como cúspide eletrônico e nuclear, a noção estatística dos modelos etc. Não existem métodos que resolvam todos os problemas. Mas, é desejável que os problemas sejam atacados pelos métodos mais precisos em qualquer escala de investigação. Como discutido acima, o método Monte Carlo Quântico apresenta deficiências que têm feito parte do dia-a-dia das atividades de vários pesquisadores ao redor do mundo. Não será um método absoluto, mas continuará sendo uma alternativa extremamente atraente por muitos anos para diversas situações onde os métodos convencionais só poderão ser empregados com recursos computacionais extremamente sofisticados e dispendiosos. AGRADECIMENTOS Ao apoio financeiro recebido da FAPESP e CNPq na forma de auxílio à pesquisa, bolsa de doutorado (FAPESP/W. F. D. Angelotti) e de iniciação científica (Pibic - CNPq/G. B. Torres). REFERÊNCIAS

1. Landau, D. P.; Binder, K.; A Guide to Monte Carlo Simulations in Statistical Physics, University Press: Cambridge, 2000. 2. Ceperley, D. M.; Alder, B. J.; Phys. Rev. Lett. 1980, 45, 566. 3. Anderson, J. B.; J. Chem. Phys. 1992, 96, 3702. 4. Anderson, J. B.; Traynor, C. A.; Boghosian, B. M.; J. Chem. Phys. 1991, 95, 7418. 5. Hammond, B. L.; Lester Jr., W. A.; Reynolds, P. J.; Monte Carlo Methods in Ab Initio Quantum Chemistry, World Scientific: Singapore, 1994. 6. Moskowitz, J. W.; Kalos, M. H.; Int. J. Quantum Chem. 1981, 20, 1107. 7. Yoshida, T.; Iguchi, K.; J. Chem. Phys. 1989, 91, 4249. 8. Bianchi, R.; Cremaschi, P.; Morosi, G.; Puppi, C.; Chem. Phys. Lett. 1988, 148, 86. 9. Manten, S.; Luchow, A.; J. Chem. Phys. 2001, 115, 5632. 10. Reynolds, P. J.; Ceperley, D. M.; Alder, B. J.; Lester Jr., W. A.; J. Chem. Phys. 1982, 77, 5593. 11. Caffarel, M.; Claverie, P.; J. Chem. Phys. 1988, 88, 1088. 12. Caffarel, M.; Claverie, P.; J. Chem. Phys. 1988, 88, 1100. 13. Ceperley, D. M.; Alder, B. J.; J. Chem. Phys. 1984, 81, 5833. 14. Pllana, S.; History of Monte Carlo Method, 2002, http:// www.geocities.com/CollegePark/Quad/2435/, acessada em Março 2007. 15. Machline, C.; Motta, I. de S.; Schoeps, W.; Weil, K. E.; Manual de Administração da Produção, Fundação Getúlio Vargas: Rio de Janeiro, 1970, vol. 2. 16. Woller, J.; The Basics of Monte Carlo Simulations, Physical Chemistry Lab, Spring: University of Nebraska – Lincoln, 1996. 17. Kalos, M. H.; Whitlock, P. A.; Monte Carlo Methods – Volume I: Basics, John Wiley & Sons, 1986. 18. Press, W. H.; Teukolsky, S. A.; Vetterling, W.T.; Flannery, B. P.; Numerical Recipes in FORTRAN 77: The Art of Scientific Computing, 2 nd ed., Cambridge University Press: Cambridge, 1997. 19. Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E.; J. Chem. Phys. 1953, 21, 1087. 20. Allen, M. P.; Tildesley, D. J.; Computer Simulation of Liquids, Oxford Science Publications: Oxford, 1990. 21. Feller, W.; An Introduction to Probability Theory and its Applications, 3rd ed., Wiley: New York, 1968, vol. 1. 22. Bahnsen, R.; Eckstein, H.; Schattke, W.; Fitzer, N.; Redmer, R.; Phys. Rev. B: Condens. Matter Mater Phys. 2001, 63, 235415. 23. Bressanini, D.; Reynolds, P. J.; Adv. Chem. Phys. 1998, 105, 37. 24. Eckstein, H.; Schattke, W.; Physica A 1995, 216, 151.

444 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36.

Angelotti et al. Pratt, L. R.; Phys. Rev. A: At., Mol., Opt. Phys. 1989, 40, 6077. Schmidt, K. E.; Moskowitz, J. W.; J. Chem. Phys. 1990, 93, 4172. East, A. L. L.; Rothstein, S. M.; Vrbik, J.; J. Chem. Phys. 1988, 89, 4880. Reynolds, P. J.; Tobochnik, J.; Gould, H.; Computers in Physics 1990, 4, 662. Anderson, J. B.; J. Chem. Phys. 1975, 63, 1499. Kosztin, I.; Faber, B.; Schulten, K.; Am. J. Phys. 1996, 64, 633. Schiff, L. I.; Quantum Mechanics, 3rd ed., MacGraw Hill, 1987. Foulkes, W. M. C.; Mitas, L.; Needs, R. J.; Rajagopal, G.; Rev. Mod. Phys. 2001, 73, 33. Barton, G.; Elements of Green’s Functions and Propagation – Potentials, Diffusion and Waves, Claredon: Oxford, 1995. Scherer, C.; Métodos Computacionais da Física, Ed. Livraria da Física: São Paulo, 2005. Grimm, R. C.; Storer, R. G.; J. Comp. Phys. 1971, 41, 1247. Szabo, A.; Ostlund, N. S.; Modern Quantum Chemistry, Dover, 1996.

Quim. Nova

37. Kato, T.; Comm. Pure Appl. Math. 1957, 10, 151. 38. Umrigar, C . J.; Wilson, K. G.; Wilkins, J. W.; Phys. Rev. Lett. 1988, 60, 1719. 39. Moskowitz, J. W.; Schmidt, K. E.; J. Chem. Phys. 1992, 97, 3382. 40. Huang, S.-Y.; Sun, Z.; Lester Jr., W. A.; J. Chem. Phys. 1990, 92, 597. 41. Wigner, E.; Phys. Rev. 1934, 46, 1002. 42. Jastrow, R.; Phys. Rev. 1955, 98, 1479. 43. Boys, S. F.; Handy, N. C.; Proc. R. Soc. London Series A 1969, 310, 63. 44. Towler, M. D.; Psi-k Newsletter 2003, 60, 166; http://psi-k.dl.ac.uk/ index.html?newsletters, acessada em Março 2007. 45. Politi, J. R. S.; Custodio, R.; J. Chem. Phys. 2003, 118, 4781. 46. Bajdich, M.; Mitas, L.; Drogny, G.; Wagner, L. K.; Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 075131. 47. Troyer, M.; Wiese, U-J.; Phys. Rev. Lett. 2005, 94, 170201.