Quatro alternativas para resolver a equação de Schrödinger para o átomo de hidrogênio

Four alternatives to solve Schrödinger equation

hydrogen atom; Schrödinger equation; numerical methods

Educação

QUATRO ALTERNATIVAS PARA RESOLVER A EQUAÇÃO DE SCHRÖDINGER PARA O ÁTOMO DE HIDROGÊNIO

Rogério Custodio, José Roberto dos Santos Politi^#, Maximiliano Segala, Roberto Luiz Andrade Haiduke

Departamento de Físico-Química, Instituto de Química, Universidade Estadual de Campinas, CP 6154, 13083-970 Campinas - SP

Márcio Cyrillo

Departamento de Física Aplicada, Instituto de Física, Universidade Estadual de Campinas, CP 6165, 13083-970 Campinas - SP

^# Endereço permanente: Instituto de Química, Universidade de Brasília, Campus Darcy Ribeiro - Asa Norte, CP 04478, 79910-900 Brasília ¾ DF.

Recebido em 25/4/01; aceito em 13/6/01

FOUR ALTERNATIVES TO SOLVE SCHRÖDINGER EQUATION. Quantum chemistry describes the hydrogen atom as one of the few systems that permits an exact solution of the Schrödinger equation. Students tend to consider that little can be learned from the hydrogen atom and forget that it can be used as a standard to test numerical procedures used to calculate properties of multielectronic systems. In this paper, four different numerical procedures are described in order to solve the Schrödinger equation for the hydrogen atom. The basic motivation is to identify new insights and methods that can be obtained from the application of powerful numerical techniques in a well-known system.

Keywords: hydrogen atom; Schrödinger equation; numerical methods.

INTRODUÇÃO

Nos quatro primeiros exemplares da revista Química Nova nos deparamos com uma série de artigos introdutórios sobre mecânica quântica¹⁴: o primeiro abordando principalmente aspectos relacionados com os primórdios da mecânica quântica¹, o segundo demonstrando rigorosamente o tratamento para o átomo de hidrogênio sem incluir efeitos relativísticos através da equação de Schrödinger², o terceiro com aplicações da mecânica quântica para sistemas simples³ e finalmente, o quarto que proporciona uma apresentação da teoria de perturbação⁴.

Com exceção do último tema, a abordagem para alunos iniciantes de química transmite a impressão de que estes são os limites confiáveis da mecânica quântica, uma vez que todos os problemas tratados apresentam soluções analíticas. O fato de apresentarem soluções analíticas também tende a sugerir erroneamente que o assunto foi esgotado dentro do nível de teoria apresentado e que, portanto, não deve haver mais nada que possa ser feito para extrair maiores informações a partir do modelo proposto.

Curiosamente, o Prof.Peixoto, autor destes quatro artigos, apresenta no segundo artigo a seguinte citação:

Chega mais perto e contempla as palavras.

Cada uma tem mil faces secretas sob a face neutra

E te pergunta, sem interesse pela resposta,

Pobre ou terrível, que lhe deres:

Trouxeste a chave?

"Procura da Poesia"

Carlos Drummond de Andrade

A sutileza desta citação associada com a apresentação de demonstrações com soluções analíticas deveriam deixar no leitor a impressão de que, embora estejamos observando um caminho rigoroso para representar e compreender o sistema em estudo, deve ser possível encontrar novos aspectos através de caminhos alternativos. Este caminho alternativo é o dos métodos numéricos e é praticamente desconhecido do aluno de química, apesar do seu poder de aplicabilidade em problemas de elevada complexidade e sem o qual não teríamos acesso a importantes aspectos teóricos, tais como o modelo orbital, propriedades de líquidos, etc.

Para demonstrar o poder que métodos numéricos podem introduzir em nosso conhecimento, procuramos neste artigo resolver a equação de Schrödinger para o átomo de hidrogênio de quatro maneiras diferentes. A primeira delas corresponde a uma rápida revisão da solução analítica e as outras três serão soluções aproximadas através de métodos numéricos, os quais nos levarão a perceber que, embora a solução da equação de Schrödinger seja muito bem conhecida para o átomo de hidrogênio, estaremos explorando novos aspectos sobre este e outros sistemas mais complexos, tais como, a possibilidade de uso de transformadas integrais para representar orbitais atômicos ou o uso de números aleatórios para a resolver a equação de Schrödinger (o método de Monte Carlo Quântico).

A Solução Analítica

Como mencionado acima, uma revisão da solução rigorosa da equação de Schrödinger para o átomo de hidrogênio pode ser encontrado na referência². Neste capítulo faremos uma rápida revisão de alguns dos aspectos mais importantes desta demonstração.

Em mecânica quântica a primeira coisa a ser preparada para modelarmos um sistema é a equação de Schrödinger. Levando-se em consideração que estamos interessados em uma descrição de um sistema em um estado estacionário, escrevemos a equação de Schrödinger independente do tempo:

sendo o operador hamiltoniano, Y a função de onda que descreve todo o sistema e E a energia de um dos estados desse sistema, que freqüentemente corresponde àquela do estado fundamental, ou seja, ao de menor energia.

O operador hamiltoniano, por sua vez, é descrito como a soma dos operadores de energia cinética e potencial :

Resolver a equação de Schrödinger corresponde a um trabalho sistemático de determinação da função de onda Y e da energia do sistema E. Conhecendo-se a função de onda pode-se determinar diversas propriedades do sistema. O controle inicial sobre a equação de Schrödinger está na construção do operador hamiltoniano e na determinação das condições de contorno que caracterizam o sistema de interesse. Desta forma, nosso primeiro passo consiste na construção dos operadores de energia potencial e cinética e, consequentemente, do operador hamiltoniano (eq.2).

O átomo de hidrogênio corresponde a um sistema de dois corpos com cargas opostas, um próton e um elétron. Estas duas cargas deslocam-se no espaço e estão sujeitas a ação de um potencial atrativo coulombiano expresso em coordenadas cartesianas por:

sendo q_p e q_eas cargas do próton e do elétron, respectivamente, e a distância entre as mesmas sendo representada por

em que podemos identificar as coordenadas do elétron (x_e,y_e,z _e) e as do próton (X, Y, Z). O operador de energia potencial⁵, corresponde a própria definição clássica desse potencial, .

O operador de energia cinética para este sistema é definido como:

ou utilizando o laplaciano (

²) para uma notação mais compacta:

sendo o primeiro termo à direita da igualdade o operador de energia cinética do elétron com massa m e o segundo termo o de energia cinética do próton com massa M. Na eq.4 ou 5 a constante , ou seja, equivale à constante de Planck dividida por 2p.

Definindo-se os operadores de energia cinética (eq.5) e potencial (eq.3), podemos estabelecer a equação de Schrödinger como:

e resolver a equação diferencial resultante para o átomo de hidrogênio. O tratamento matemático não é trivial e usualmente procuramos transformar o problema de duas partículas (comumente chamado de problema de dois corpos) em dois problemas de uma partícula (ou problema de um corpo). A manipulação da eq.6 pode nos levar a duas equações diferenciais que descreverão o movimento eletrônico ao redor do núcleo e o movimento de translação do átomo como um todo². Uma vez que estamos interessados particularmente no problema da distribuição eletrônica, vamos focalizar nossa atenção na equação diferencial que controla este evento.

Neste caso, o próton é colocado permanentemente no centro do sistema de coordenadas, (X=0,Y=0,Z=0), como se possuísse uma massa infinita, e o elétron se deslocasse ao redor do mesmo, representado por coordenadas cartesianas em relação ao próton iguais a (x,y,z), como se possuísse uma massa reduzida (µ) igual a:

A equação de Schrödinger para este sistema é escrita simplesmente como:

sendo o laplaciano é uma função apenas das coordenadas do elétron com massa µ e a função potencial reduz-se a:

De maneira mais explícita, a eq.8 será escrita como:

Embora tenhamos uma equação diferencial mais simples do que a da eq.6, para resolvê-la necessitamos ainda de alguma manipulação matemática. Uma das primeiras técnicas empregadas para resolver a eq.10 consiste na transformação das coordenadas cartesianas (x,y,z) em coordenadas esféricas polares (r,q,f). A relação existente entre estes dois sistemas de coordenadas pode ser verificado na Figura 1 e pelas expressões apresentadas na Tabela 1. Neste novo sistema de coordenadas, o operador de energia potencial converte-se em:

Enquanto que o operador de energia cinética pode ser representado como:

O que leva a equação de Schrödinger a ser escrita como:

A aparência mais complicada da eq.13, em relação a eq.10, é compensada pela possibilidade de realizarmos uma separação de variáveis, se admitirmos que a função de onda Y(r,q,f) pode ser fatorada como:

Caso a função de onda possa ser escrita como o produto de três funções que dependam exclusivamente de uma única variável, como na eq.14, a eq.13 poderá ser separada em três equações diferenciais ordinárias², que são:

sendo b e m_l constantes arbitrárias a serem determinadas.

Por inspeção, verificamos que uma possível solução para a eq.15 seria a de uma função matemática que derivada duas vezes em relação a variável f produzisse a própria função multiplicada por uma constante. Algumas funções matemáticas que apresentam esta propriedade poderiam ser sugeridas e podemos considerar que uma forma conveniente de representar a função de onda F é:

A função apresentada na eq.18 é complexa e levando-se em conta que a função de onda deve ser finita, unívoca e contínua, podemos verificar que esta representará o comportamento de uma função de onda aceitável ou bem comportada se m_l=0, ±1, ±2,...

A eq.16, por outro lado, apresenta uma grande semelhança com equações diferenciais bem conhecidas na literatura, denominadas de equações associadas de Legendre⁷. Para que a semelhança entre a eq.16 e as equações de Legendre seja completa devemos admitir que b = l(l + 1), sendo l=0, 1, 2, ..., |m_l|. Desta forma, por analogia com as equações de Legendre podemos dizer que as funções de onda dependentes de q apresentam a seguinte forma:

sendo funções conhecidas como polinômios associados de Legendre⁷. As restrições para que as funções Q existam estão vinculadas aos números inteiros l e m_l, ou seja, só existirão funções Q que apresentem valores de l e m_l que satisfaçam as seguintes restrições:

O produto das funções de onda Q e F corresponde a uma nova função que surge frequentemente em diferentes problemas físicos e é denominada de harmônicos esféricos, Finalmente, para a eq.17, que depende exclusivamente da variável r, a solução é conhecida e corresponde a equações diferenciais, que na literatura são conhecidas como equações associadas de Laguerre⁷. As soluções das equações de Laguerre possuem funções características que são representadas matematicamente por:

Nesta expressão,

sendo Z o número atômico, , n são números inteiros iguais a 1, 2, ..., e funções polinomiais conhecidas como polinômios associados de Laguerre⁷. Da mesma forma que ocorreu com a função Q, as únicas funções radiais R aceitáveis são aquelas cujos valores de l não são superiores a (n-1).

Desta forma, os estados eletrônicos possíveis para o elétron em um átomo de hidrogênio serão descritos por funções de onda obtidas através do produto das eqs.18, 19 e 21, que resultam em:

sendo N_n,l,mlum termo incluindo todas as constantes apresentadas nas eqs.18, 19 e 21 e x = 2Z/na₀.

Estes estados eletrônicos são caracterizados por energias orbitais que são determinadas através da comparação da eq.17 com as equações associadas de Laguerre. Sua forma geral corresponde a:

em que verificamos que a diferenciação da energia dos diferentes estados eletrônicos do átomo de hidrogênio ocorre apenas em termos de n². Considerando-se que para um determinado valor de n, podemos ter diferentes valores de l até um máximo de (n-1) e para cada valor de l, valores de m_l variando de 0, ±1, ±2, ..., ±l, percebemos que diversos estados deste átomo são degenerados. As propriedades desses estados eletrônicos são frequentemente discutidas qualitativamente em diferentes cursos de graduação e a demonstração da solução da equação de Schrödinger para o átomo de hidrogênio é considerada muitas vezes apenas como uma forma de justificar o surgimento e uso de regras de classificação de orbitais atômicos e uma maneira de caracterizar as propriedades da distribuição eletrônica em cada um destes estados eletrônicos.

Entretanto, uma vez que conheçamos a solução analítica de um problema, podemos posteriormente utilizá-la como uma referência para testar métodos alternativos de cálculo para sistemas mais complicados. Três alternativas serão abordadas a seguir e demonstrarão que em algumas circunstâncias poderemos obter novas informações, que eventualmente poderão ser utilizadas em diferentes sistemas.

O Método de Noumerov

A idéia fundamental por trás de uma solução numérica para uma equação diferencial é produzir uma tabela de valores das variáveis envolvidas no problema e condições iniciais que caracterizem restrições ao sistema, também denominadas de condições de contorno (como por exemplo, sugerir que uma função tenda a zero quando sua variável tender ao infinito: F(r ® ¥) = 0) e, através de alguma técnica iterativa, encontrar o valor da função para toda a tabela de valores das variáveis.

A equação diferencial para o átomo de hidrogênio pode ser resolvida numericamente e a solução para cada um dos estados do átomo de hidrogênio pode ser obtida de uma forma aproximada, mas com uma precisão tão grande quanto necessária.

Para que isto seja feito, aproveitamos a sugestão de separação de variáveis e desmembramos a equação de Schrödinger em suas três equações diferenciais dadas pelas eqs.15, 16 e 17. Observando-se estas três expressões, verificamos que as eqs.15 e 16 não dependem da função potencial e poderiam ser encontradas em problemas completamente diferentes daquele que estamos interessados. Para ilustrar a solução numérica que caracterize particularmente o átomo de hidrogênio, ou seja, que dependa do potencial, vamos nos concentrar apenas na solução numérica da equação diferencial que depende da variável radial r (eq.17).

Uma forma consideravelmente mais simples para a eq.17 pode ser obtida se substituirmos a função R_i, que corresponde a função que desejamos determinar, por uma outra função P_i, que está relacionada com R_i pela seguinte expressão:

O uso da eq.25 reduz a eq.17 a:

sendo:

A expressão radial para o átomo de hidrogênio (eq.17) tornou-se, com esta transformação, em uma equação diferencial ordinária de segunda ordem (eq.26). Resolvendo-se esta equação diferencial, determinaremos a função P_m(r) e a partir da eq.25 determinaremos a função R_m(r). Para que a solução numérica da eq.26 seja encontrada, precisamos inicialmente de uma expressão simplificada para a derivada segunda de uma função. A representação desta aproximação define o método numérico em si e também o seu desempenho em termos de estabilidade numérica e de convergência. Há uma vasta literatura básica tratando os casos mais utilizados^8-10.

Para criarmos uma expressão aproximada para a derivada segunda, podemos utilizar a expansão em série de Taylor de uma função. Considerando-se uma função arbitrária P(r), a série de Taylor nos permite estimar o valor desta função em um ponto r+h ou r-h, onde h define um intervalo pequeno na vizinhança de um ponto r a partir das equações:

Se efetuarmos a soma das eqs.28 e 29 eliminamos as derivadas primeira, terceira, quinta, etc. e preservamos as derivadas segunda, quarta, sexta, etc, obtendo uma expressão que poderá ser escrita como:

Uma maneira de determinar a derivada segunda aproximadamente seria desprezar todas as derivadas presentes na eq.30 diferentes da derivada segunda. Assim, poderíamos considerar a seguinte aproximação:

que nos diz que a derivada segunda de P(r) pode ser estimada com o conhecimento do valor desta função em dois pontos distintos e vizinhos a r.

Se observarmos a eq.30, verificaremos que ao isolarmos a derivada segunda (eq.31) a derivada quarta será multiplicada por um termo h² e as outras derivadas superiores nesta expressão serão multiplicadas por valores de h elevados a potências maiores. Desta forma, se o valor de h for muito pequeno (próximo de zero), h² será um número muito menor e consequentemente os outros valores de h elevados a potências maiores produzirão números ainda mais insignificantes. Baseados nesta informação, podemos dizer que o maior erro que deveremos encontrar na aproximação utilizada para produzir a eq.31 será proporcional a h².

Expansões em série de Taylor como aquelas apresentadas nas eqs.28 e 29, podem ser obtidas para intervalos 2h em vez de h. A combinação das expansões em série de potência a partir de somas ou subtrações podem produzir representações aproximadas para a derivada segunda com erros da ordem de h⁴. Como exemplo desta possibilidade, a equação abaixo foi obtida empregando-se esta técnica:

A eq.32 corresponde a uma representação mais precisa para a derivada segunda. A determinação do valor da derivada segunda através desta expressão depende do conhecimento da função P em cinco pontos diferentes {r+2h,r+h,r,r-h,r-2h}.

Embora imaginemos que a escolha pela eq.32 em relação a eq.31 seja óbvia para uma determinação mais precisa de P, devemos levar em consideração aspectos que dependem do processo de utilização destas aproximações. Um dos inconvenientes da eq.32 está na sua instabilidade numérica, que será elucidada abaixo, ou seja, no mecanismo de obtenção das funções P.

A busca por expressões que possuíssem os menores erros possíveis e que se mantivessem estáveis durante uma simulação numérica levou, entre outros, ao desenvolvimento do método de Noumerov^11-13. Apenas para ilustrar a popularidade deste método, Douglas Hartree, um dos pioneiros no cálculo de propriedades de átomos multieletrônicos, optou por este método durante o processo de solução das equações diferenciais que levam seu nome¹⁴. Este método é aplicável apenas às equações diferenciais de segunda ordem do tipo da eq.26. Assim como definimos uma função P(r) arbitrária relacionada a função radial R(r), Noumerov definiu uma nova função Q(r) da seguinte forma:

que também pode ser expandida em série de Taylor como nas eqs.28 e 29, resultando em:

que originará, de maneira análoga à eq.30, a seguinte expressão:

onde vemos que o termo em quarta ordem desapareceu, sendo o erro local agora proporcional a h⁶. Desprezando-se as derivadas de ordem superior a 2 teremos como aproximação para a derivada segunda a expressão:

Como podemos ver, a obtenção da derivada segunda pode ser realizada a partir das eqs.31, 32, 36 e de uma série de outras expressões que podem ser desenvolvidas com a preocupação de minimizar o erro na determinação da derivada segunda e também de fornecer soluções estáveis.

Tendo conhecimento das possíveis maneiras de representarmos a derivada segunda, podemos verificar como deveremos proceder para obter as funções radiais a partir da aproximação escolhida. Qualquer uma das expressões acima (eqs.31, 32 ou 36) poderá ser substituída na eq.26 e posteriormente definimos uma sistemática para determinarmos valores de Q ou P em função da coordenada r. Considerando a importância do método de Noumerov e sua maior complexidade em relação as outras duas aproximações apresentadas, nos concentraremos neste processo, o que permitirá ao leitor ter conhecimento suficiente para a solução das outras alternativas.

Inicialmente rearranjamos a eq.36 da seguinte maneira:

A derivada segunda de P será então substituída pela igualdade apresentada na eq.26 produzindo:

Por definição a função Q(r) deve ser escrita de acordo com a eq.33. Entretanto, a eq.33 pode ser obtida como função de P se substituirmos a eq.26 na eq.33, o que nos possibilita escrever Q como:

Esta definição de Q possibilita substituir a função P na eq.38 produzindo a seguinte expressão:

A eq.40 corresponde a expressão fundamental para a obtenção das soluções da equação diferencial (eq.26). O procedimento de obtenção da função Q, P e R segue os seguintes passos:

1. Definimos um conjunto arbitrário de valores de r igualmente espaçados pela distância h;

2. Definimos um valor arbitrário, mas que acreditemos esteja próximo do valor correto, para a energia do orbital E_m que desejamos descrever através da função R(r);

3. Calculamos os valores de g(r) através da eq.27 e tabelamos os valores ;

4. Construímos três outras colunas nesta tabela, uma para acomodar os valores de Q, outra para os valores de P e outra para os valores de R. Os valores de P serão obtidos a partir dos valores de Q pela eq.39 e os de R serão obtidos através de P pela eq.26;

5. Para determinarmos os valores de Q, partimos de dois valores iniciais de Q(r) e Q(r-h) e empregando-se a eq.40, encontramos o valor para o ponto Q(r+h). Uma vez determinado o valor de Q(r+h) podemos utilizar a eq.40 substituindo-se agora o valor de Q(r+h) no lugar de Q(r) e o de Q(r) no lugar de Q(r-h), determinando o valor de Q(r+2h). Este processo é então repetido para a determinação dos valores de Q(r+3h), Q(r+4h), etc. e assim sucessivamente, sempre se utilizando os dois valores de Q anteriores para determinação do valor de Q na coordenada de interesse;

6. Ao completarmos nossa tabela de valores das funções radiais devemos analisar a sua validade através das condições que normalmente esperamos para funções de onda bem comportadas: regularidade em todo intervalo, continuidade para a função e sua primeira derivada, etc. O número de nós das funções obtidas permite verificar e ordenar os estados de acordo com suas energias (autovalores). A função de onda correta deve ter um comportamento convergente nos limites da coordenada r. Caso a função que obtivemos apresente um comportamento divergente, isto indica que não obtivemos um auto-estado do sistema que estamos estudando e deveremos reiniciar todo o processo escolhendo-se novos valores para os parâmetros iniciais, ou seja, Q(r), Q(r-h) e E_m.

Este procedimento esconde em sua simplicidade uma série de dificuldades para sua implementação prática que, devemos ressaltar, pode ser realizada utilizando-se o mais simples dos computadores ou mesmo uma calculadora.

A primeira e mais complexa pergunta que podemos nos fazer sobre esta simulação é: como encontrar os dois primeiros pontos, Q(r) e Q(r-h), para gerarmos todos os valores de Q, P e R para a série de pontos r previamente escolhidos? Esta pergunta não possui resposta simples e muitas vezes devemos optar pela experimentação, ou seja, a repetição exaustiva do algoritmo acima com diversos valores iniciais até que obtenhamos as funções desejadas. Este é o maior preço a ser pago pela simplicidade e velocidade do método (Hartree conseguiu excelentes resultados¹⁴ para funções de onda atômicas com, em média, apenas 20 iterações do método de Noumerov).

Outro problema é a localização dos autovalores, ou seja, das energias dos orbitais, E_m, que pode ser contornado se possuirmos algum conhecimento prévio sobre o sistema que estamos estudando, por exemplo, dados experimentais sobre energias de ionização.

Finalmente temos o sério problema de propagação de erros. As eqs. 30, 32 e 35 mostram que o erro depende do estado estudado. Por exemplo, para a eq.35, temos que o erro total será proporcional à integração do erro local em todo o intervalo, sendo dado por:

Este comportamento implica no aumento de erro para o cálculo de estados excitados de alta ordem¹⁴. Entretanto, considerando-se a simplicidade do mecanismo numérico de solução da equação de Schrödinger, uma série de técnicas podem ser encontradas na literatura que possibilitem a sua solução minimizando uma parte significativa de erros provenientes das aproximações utilizadas como, por exemplo, a correção de Cooley ¹¹.

O Método da Combinação Linear das Funções de Base

Considerando-se as dificuldades de se resolver numericamente a equação de Schrödinger para sistemas moleculares e ainda, levando-se em consideração a possibilidade de uma solução aproximada empregando-se o método de Hartree-Fock¹⁴ para sistemas multieletrônicos, C.C.J.Roothaan, em 1951, estabeleceu uma associação formal do método de Hartree-Fock com um outro modelo matemático denominado de combinação linear dos orbitais atômicos (LCAO/Linear Combination of Atomic Orbitals)¹⁵.

Como vimos anteriormente, a solução da equação de Schrödinger consiste na determinação da função de onda e da energia do estado do sistema que estamos estudando. Para um sistema como o átomo de hidrogênio, estaremos procurando a representação dos orbitais atômicos e as suas respectivas energias. O método LCAO, que deve ser denominado de forma mais correta de método de combinação linear de funções de base, consiste em substituir parte da informação matemática que desconhecemos, por outras que possamos ter algum controle e que nos permita um método de determinação das propriedades que estamos interessados. No caso do átomo de hidrogênio, o método LCAO consiste em substituir os orbitais hidrogenóides, f_i, por uma combinação linear de funções de base matemáticas c_j. Matematicamente, estamos considerando que:

sendo c_jm coeficientes da expansão do m-ésimo orbital hidrogenóide, que correspondem a uma medida de quanto cada uma das funções de base c_j colaboram na constituição desse orbital.

A eq.41 nos mostra que trocamos uma função desconhecida, Y_H,m, por três outras informações também desconhecidas. Em primeiro lugar, a eq.41 não nos impõe qualquer restrição sobre qual deve ser o tipo de função de base c_j a ser utilizada; em segundo, não é definido o número k de funções que deverá ser empregado para caracterizar corretamente um orbital hidrogenóide e em terceiro, nenhuma informação é oferecida com relação a como serão determinados os coeficientes de combinação linear c_jm.

Embora tenhamos inicialmente a impressão de que o problema da determinação das funções orbitais ficou mais complicado, em geral, informações complementares nos auxiliam na definição de respostas que podem proporcionar um nível considerável de precisão na determinação da solução da equação de Schrödinger.

Inicialmente, o controle sobre quais são as funções de base que eventualmente podem ser utilizadas para representar um orbital hidrogenóide, pode ser definido pelo primeiro postulado da mecânica quântica. Neste postulado são apresentadas condições para que uma função de onda possa ser considerada como aceitável^1,5. Desta forma, poderemos escolher qualquer expressão matemática como função de base, desde que ao construirmos a função orbital, a função de onda seja uma função bem comportada.

O número destas funções de base que devem ser combinadas para produzir uma representação exata da função orbital deve tender a infinito. Na prática nunca utilizamos um conjunto infinito de funções e no caso do método LCAO, a precisão dos resultados estará ligada fundamentalmente ao número de funções que empregaremos em nossos cálculos.

Finalmente, os valores dos coeficientes de combinação linear poderão ser determinados a partir do princípio variacional^15,16. Através deste critério, os valores dos coeficientes de combinação linear ou de outros parâmetros que eventualmente possam ser incorporados nas funções de base, podem ser ajustados de tal forma que a energia eletrônica do sistema seja mínima. Segundo o método variacional, a melhor função de onda tentativa nunca apresentará uma energia menor do que a função de onda exata. Desta forma, se calcularmos a energia de um estado m qualquer com uma função de onda tentativa Y_m, independente da maneira como construímos esta função bem comportada, a menor energia que poderemos obter será na melhor das hipóteses o resultado exato.

A definição do tipo de função de base, do número de funções a serem utilizadas e a determinação dos coeficientes de combinação linear constituem a essência do método LCAO. Enquanto os dois primeiros fatores dependem de uma escolha arbitrária, a determinação dos coeficientes de combinação linear apresenta uma sistemática muito bem conhecida e que será explorada com maiores detalhes a seguir.

Para determinarmos os valores dos coeficientes de combinação linear, devemos antes de mais nada encontrar uma forma para determinar a energia do sistema. Considerando-se que temos uma função de onda tentativa, é muito pouco provável que esta função de onda seja auto-função do operador hamiltoniano (satisfazendo a eq.1) e, desta forma, a única maneira de determinarmos a energia do sistema é através do teorema do valor médio (3^o postulado da mecânica quântica) definido pela expressão:

sendo t a representação das coordenadas generalizadas.

Substituindo-se a eq.41 na eq.42 e rearranjando, encontramos:

que pode ser escrita em uma notação mais compacta como:

sendo

O mínimo de energia em relação aos coeficientes é encontrado fazendo-se com que a derivada primeira da energia em relação a cada coeficiente seja igual a zero:

Derivando-se a eq.44 em relação, por exemplo, ao coeficiente c_j e rearranjando-se, obtemos a seguinte expressão:

Para cada um dos coeficientes de combinação linear teremos uma equação semelhante à eq.46. A única diferença estará na função de base que estará sendo utilizada para calcular as integrais H_ij e S_ij. De maneira geral, teremos a seguinte situação, em uma forma mais explicita, depois de realizarmos todas as derivadas de E_m em relação a cada c:

A primeira linha foi obtida derivando-se a eq.44 em relação a c₁, a segunda linha foi obtida em relação a c₂ e assim sucessivamente.

A série de equações representadas acima pode ser escrita na forma matricial como:

que normalmente é encontrada na literatura em uma notação mais compacta como:

Esta expressão é freqüentemente denominada em mecânica quântica de equação secular e sua solução é muito bem conhecida em cálculo numérico. Os passos para sua resolução correspondem a aplicar algum tipo de transformação que faça com que a matriz S seja convertida em uma matriz identidade, 1. Usualmente esta transformação consiste em encontrar uma matriz A capaz de fazer com que a seguinte relação seja verdadeira: A⁺SA = 1, sendo A⁺ a adjunta de A. Esta transformação corresponde a fazer com que todas as funções de base empregadas na combinação linear sejam ortogonais entre si, quando A⁺=A^-1.

Ao ortogonalizarmos as funções de base, a eq.51 adquire a seguinte expressão:

que pode ser re-escrita como:

Resolver a eq.51 corresponde a determinar os valores de C' e E_m. A eq.51 apresenta uma solução trivial que não acrescentará nenhuma informação física razoável ao nosso sistema, ou seja, considerar que todos os coeficientes sejam iguais a zero. Uma segunda alternativa corresponde a fazer com que:

Em outras palavras, uma solução não-trivial corresponde a fazer com que o determinante da eq.52 seja igual a zero. A solução para este determinante faz com que tenhamos condições de determinar o valor de E_m. Para obtermos os valores dos coeficientes c, substituímos o valor de E_m adequado na eq.51 e determinamos os valores de C'. Para obtermos os valores de C, empregamos a matriz de ortogonalização através da expressão: AC' = C.

Uma outra maneira usualmente empregada em programas de mecânica quântica corresponde a empregar métodos que diagonalizem a matriz H' na eq.50. A matriz de diagonalização de H' será a matriz C'. Todos estes métodos são muito bem estabelecidos e empregados rotineiramente em cálculos mecânico-quânticos^7-9.

A obtenção dos coeficientes de combinação linear em cálculos quânticos de sistemas multieletrônicos pode proporcionar informações extremamente relevantes, principalmente quando tratamos sistemas moleculares. Estes coeficientes podem ser relacionados à distribuição de carga em diferentes regiões da molécula. Conhecido como análise populacional de Mulliken¹⁷, este método pode ser considerado como um dos mais simples e populares para determinação da densidade eletrônica sobre os átomos ou ligações químicas em uma molécula. Entretanto, os coeficientes de combinação linear podem esconder informações que normalmente passam despercebidas, mesmo por especialistas na área. Para explorarmos esta possibilidade, vamos retornar ao átomo de hidrogênio e verificar que novas informações podem ser obtidas a partir de cálculos convencionais.

Para a construção das funções do átomo de hidrogênio poderemos empregar a eq.41 e para isto precisamos definir uma expressão matemática para as funções de base. A solução analítica da equação de Schrödinger para este átomo no estado fundamental proporciona uma energia eletrônica igual a 0,5u.a. com uma função orbital 1s que pode ser escrita em coordenadas polares esféricas como:

sendo a constante de normalização

Uma das alternativas para representar a função 1s aproximadamente através da eq.41 seria escolher, por exemplo, funções de base do tipo gaussianas, ou seja:

sendo a constante de normalização N' = (2a/p)^¾. A função dada na eq.54 não é autofunção do operador hamiltoniano e se utilizada para representar a energia do átomo de hidrogênio individualmente proporcionará um resultado desastroso. Podemos verificar que esta função apresenta a mesma simetria do orbital Y_1s, ou seja, apresenta a mesma dependência angular da função dada na eq.53 e, portanto, corresponde a uma função de base adequada para a obtenção da função Y_1s. Desta forma, a função de onda para o átomo de hidrogênio pode ser escrita como:

ou seja, combinamos diferentes funções gaussianas para produzir a função Y_1s. Como sabemos, precisamos determinar o número k de funções de base e os coeficientes de combinação linear c na eq.55. Entretanto, observando-se a eq.55 verificamos que desconhecemos uma outra informação, os expoentes das funções gaussianas, a_i. Estes expoentes podem ser determinados variacionalmente empregando-se o mesmo critério adotado para os coeficientes de combinação linear. Uma dificuldade prática na otimização destes expoentes está no fato de que os mesmos não são parâmetros lineares e, portanto, necessitam de métodos consideravelmente mais caros computacionalmente para serem obtidos. Uma alternativa simples para a determinação de todos os k expoentes das gaussianas consiste em estabelecer uma regra através de expressões matemática que forneçam o valor aproximado dos expoentes, em relação aqueles que seriam obtidos se os mesmos fossem calculados através da minimização de energia em relação a todos os expoentes da série que está sendo utilizada. Um exemplo de tal série aproximada pode ser:

A eq.56 mostra que todos os k expoentes a serem empregados com as gaussianas podem ser obtidos se determinarmos os valores de a e b. Novamente, um critério que pode ser empregado para a determinação de a e b pode ser o ajuste variacional. A vantagem de se empregar a eq.56 em relação à otimização de todos os expoentes é óbvia. Na otimização de todos os expoentes deveríamos realizar a aplicação de técnica apropriada para a obtenção de todos os k expoentes ótimos de todas as gaussianas, enquanto que a utilização da eq.56 restringe a otimização a apenas 2 parâmetros para a obtenção de todos os k expoentes da série. Certamente que utilizar a eq.56 deverá proporcionar algum prejuízo numérico ao cálculo da energia e da própria função de onda da molécula, porém isto poderá ser compensado empregando-se um número maior de funções gaussianas para representar a função Y_1s.

Considerando-se que o problema da determinação dos expoentes pode ser resolvido através da otimização dos parâmetros a e b, seria interessante verificarmos a dependência do cálculo da energia do átomo de hidrogênio no estado fundamental com o número de funções gaussianas. Na tabela 2 são apresentados os valores da energia para o átomo de hidrogênio, calculados com diferentes números de funções gaussianas com expoentes otimizados empregando-se a eq.56 e obtendo-se os coeficientes de combinação linear através dos métodos discutidos acima (por exemplo, eq.52). Os resultados mostram que o uso de um número reduzido de funções gaussianas compromete seriamente a qualidade das energias calculadas. Também podemos perceber que ao aumentarmos o número de funções gaussianas, a energia eletrônica do hidrogênio converge para o valor exato de 0,5u.a.. Esta informação reforça o fato de que uma representação adequada, empregando-se a combinação linear de funções de base proporciona o resultado exato se o número de funções tende ao infinito.

A possibilidade de desenvolver a série de expoentes para as funções gaussianas empregando séries geométricas como a da eq.56 nos faz pensar se não existe algum tipo de função matemática que forneça a tendência dos coeficientes de combinação linear do respectivo orbital molecular. A resposta para esta questão pode ser determinada de diferentes maneiras. Optamos pela via dos testes computacionais e para isto devemos observar atentamente a eq.55. Nesta equação, considerando que estamos utilizando funções gaussianas de mesma simetria e que cada função gaussiana é distinguida de todas as outras pelo valor do expoente a. Podemos dizer que, uma vez que cada coeficiente está associado a cada uma das gaussianas, estes também devem ser uma função dos respectivos expoentes, ou seja, c_i=c_i(a_i ). Desta forma, seria interessante verificar se podemos observar alguma correlação matemática entre os coeficientes de combinação linear com os respectivos expoentes.

Das diversas possibilidades de correlacionar estas duas quantidades, a Figura 2 nos mostra uma que sugere um comportamento matemático característico dos coeficientes em relação ao logaritmo neperiano dos respectivos expoentes. Nesta figura foram incluídos coeficientes obtidos com diferentes números de gaussianas e podemos verificar uma tendência bem definida no comportamento dos coeficientes à medida que o número de gaussianas aumenta. A expansão da função 1s realizada com 3 gaussianas (3G) apresenta uma tendência quase linear em relação ao lna. Para a expansão envolvendo 5 funções gaussianas (5G) verificamos um desvio da linearidade com uma possível tendência dos coeficientes na região de lna < 0 para valores apontando para a abcissa. As funções de base seguintes (6G, 8G, 10G, 15G e 20G) apresentam o comportamento de uma função contínua e que tende a zero nos limites de lna. A Figura 2 mostra que a amplitude da possível função contínua diminui à medida que o número de gaussianas é aumentado na expansão do orbital Y_1s.

O comportamento dos coeficientes de combinação linear frente aos respectivos expoentes e a representação da função Y_1s pela eq.55 sugere que provavelmente uma maneira exata de representar as funções do átomo de hidrogênio seja substituir o somatório da eq.55 por uma integral se fizermos com que o número de gaussianas tenda ao infinito. Assim, uma nova representação para Y_1s seria:

Nesta equação, os coeficientes de combinação linear c_i foram substituídos por uma função contínua f_i(a), denominada de função peso. A função de base continua sendo uma função gaussiana, embora possamos empregar outros tipos de funções que levem a função de onda a ser bem comportada, e o limite de integração é definido de maneira a abranger todo o espaço possível de expoentes a. Pode-se realizar uma transformação de coordenadas de tal forma que a integral definida na eq.57 seja re-escrita como:

Tanto a eq.57 quanto a eq.58 podem ser empregadas para representar a função orbital exatamente. Estas formas de representar as funções orbitais são familiarmente identificadas como transformadas integrais^8,9 e foram sugeridas independentemente em dois modelos matemáticos equivalentes: a) o método de bases even tempered, por Raffenetty¹⁸ e Feller e Ruedenberg¹⁹ e b) o Método da Coordenada Geradora por Griffin, Hill e Wheeler^20,21, sendo o modelo proposto por estes últimos autores um caso mais genérico do que o propostos pelos autores das bases even tempered. As condições para termos uma representação adequada de uma função orbital através da eq.57 ou eq.58 é de que a função peso f_i(a) seja uma função contínua, suave e com comportamento assintótico tendendo a zero nas extremidades dos limites de integração, exatamente o tipo de tendência que observamos na Figura 2.

Neste ponto, valeria a pena estabelecer uma correlação quantitativa entre a função peso f_i(a), que representa exatamente a função orbital, com os coeficientes de combinação linear c_i. A diferença óbvia entre as duas quantidades é que uma representa informações provenientes de um espaço contínuo, enquanto a outra apresenta a mesma representação em um espaço discreto. Desta forma, para interconvertermos ambas as quantidades poderemos, por exemplo, discretizar a eq.57 ou eq.58, ou seja, identificá-las com uma integral que possa ser resolvida numericamente.

Tomando-se a eq.58 para o processo de integração numérica, pode-se integrá-la aproximadamente subdividindo-se o eixo lna em pequenos intervalos idênticos com um intervalo igual a Dlna. Quanto menor o intervalo em Dlna, mais precisa será a determinação da integral numericamente. Considerando-se esta subdivisão de lna, a eq.58 pode ser escrita como:

Comparando as eq.59 com a eq.55 verifica-se que a relação entre os coeficientes da expansão e a função peso é dada por:

o que nos permite determinar a função peso discretizada através da equação:

A eq.61 sugere que o gráfico apresentado na Figura 2 seja elaborado novamente através de um gráfico de f_1s(a) vs ln(a) ao invés de c_i vs ln(a). Uma vez que os conjuntos de base foram elaborados empregando-se a eq.56, a determinação do Dlna pode ser obtida através do ln(b) para cada conjunto de gaussianas empregado nos diferentes cálculos. Assim, dividindo-se cada um dos coeficientes de combinação linear do cálculo com três gaussianas pelo respectivo ln(b), e procedendo-se da mesma forma para todos os outros cálculos realizados com os respectivos ln(b), pode-se traçar o gráfico apresentado na Figura 3.

A Figura 3 mostra que o comportamento de f_i é característico do orbital que está sendo analisado e que independe do número de funções com os quais o cálculo foi realizado, uma vez que todos os pontos de qualquer um dos cálculos estão exatamente sobre a mesma curva. Este comportamento regular e específico não é exclusivo do orbital atômico 1s do hidrogênio. É possível obtermos através dos coeficientes de combinação linear informações sobre a verdadeira representação do orbital atômico ou molecular de sistemas multieletrônicos de maneira semelhante a realizada para o átomo de hidrogênio. O conhecimento de que esta função existe para cada orbital atômico ou molecular permite avaliar a qualidade de funções de base nesses sistemas e o comportamento assintótico das funções peso auxiliam no desenvolvimento ou modificações em conjuntos de funções de base finitos^22-32.

Devemos chamar a atenção para o fato de que, embora o uso de funções gaussianas possibilitem a descrição de funções orbitais através de uma transformada integral, isto não necessariamente ocorrerá caso sejam utilizados outros tipos de funções de base.

O Método Monte Carlo Quântico

O último modelo a ser discutido neste artigo para resolver a equação de Schrödinger para o átomo de hidrogênio está baseado em uma categoria de métodos denominados de estocásticos³³, mais particularmente, no método de Monte Carlo³⁴. Esta denominação deve-se ao fato das simulações de Monte Carlo utilizarem seqüências de números aleatórios em seu processo de cálculo. Sua popularidade está associada com a simplicidade das técnicas numéricas envolvidas e, consequentemente, os campos de aplicação são extremamente variados, podendo-se citar por exemplo: desenvolvimento de reatores nucleares, cromodinâmica quântica, radioterapia, fluxo de tráfego, evolução estelar, exploração de petróleo, etc³⁵.

Nosso interesse particular neste trabalho, corresponde a empregar esta técnica para resolver a equação de Schrödinger. Nesta categoria de aplicações, o método de Monte Carlo recebe a designação particular de método Monte Carlo Quântico. De maneira geral, existem duas linhas distintas baseadas nesta possibilidade: a) o método Monte Carlo Quântico Variacional e b) o método Monte Carlo Quântico de Difusão. Independente da maneira como o método é empregado nestes dois casos, ele normalmente estará procurando resolver uma equação integral e não uma diferencial.

Mais especificamente, o método Monte Carlo Variacional procura determinar o valor médio de qualquer propriedade atômica ou molecular, conhecendo-se uma função de onda arbitrária que possa descrever o sistema, integrando-se a equação:

A característica deste método está na maneira como a integral da eq.62 é resolvida.

De nosso conhecimento de cálculo, a integral definida de uma função g(x) pode ser obtida como:

se admitirmos que o intervalo Dx é constante, seu valor poderá ser determinado como:

ou seja, o valor aproximado da integral dada pela eq.63 é determinado como o valor médio da função g(x) multiplicado pela diferença entre os limites de integração. Neste caso particular, consideramos que o espaçamento entre os valores de x são constantes. Porém, a média de qualquer propriedade pode ser definida com qualquer conjunto de valores de x. O método de Monte Carlo considera que a variável x é determinada aleatoriamente. Em outras palavras, escolhe-se um número N de valores de x aleatórios e determina-se o valor de g(x) para cada variável. A integral da função g(x) será a média de todos os valores calculados dentro do intervalo definido por b e a.

Desta forma, a eq.62 pode ser resolvida de maneira semelhante se tivermos conhecimento da função de onda que descreve o sistema. Uma vez que a função de onda corresponde a algo que pretendemos determinar, o método de Monte Carlo deve ser considerado variacional no sentido de que poderemos criar qualquer tipo de função de onda que imaginarmos. Na melhor das hipóteses estaremos empatando com a função de onda exata do sistema. Se o operador na eq.62 for o hamiltoniano, esta hipótese será confirmada pelo valor da energia calculada para cada função de onda tentativa que desenvolvermos.

Um outro aspecto que deve ser comentado é que, na prática, o método de amostragem para a determinação do valor aproximado da integral é definido por sistemas mais sofisticados do que o sugerido no exemplo de integração da função g(x). Usualmente, as simulações de Monte Carlo estão associadas a um mecanismo bastante eficiente de amostragem denominado de método de Metrópolis³⁶.

O método de Monte Carlo de Difusão, por sua vez, baseia-se na grande semelhança existente entre a equação de Schrödinger dependente do tempo e a equação que descreve um processo de difusão. A equação de Schrödinger em sua forma dependente do tempo para um sistema unidimensional é dada por:

sendo m a massa reduzida da partícula e V(x) o potencial a que a partícula está submetida.

Por outro lado, um processo de difusão pode ser representado pela seguinte equação diferencial:

sendo C a concentração da substância sofrendo difusão e D e k são duas constantes, das quais a primeira é conhecida como coeficiente de difusão.

Comparando-se a eq.65 com a eq.66 verifica-se uma semelhança considerável e podemos imaginar que, em princípio, a equação de Schrödinger poderia ser resolvida através de simulações que sejam caracterizadas pela eq.66. Uma das poucas diferenças entre as duas equações é que a equação de Schrödinger (eq.65) apresenta uma componente imaginária. No sentido de aumentar a semelhança entre as duas expressões podemos definir uma unidade de tempo imaginária . Após transformarmos a coordenada temporal da eq.65 para esta nova escala de tempo obtemos:

A simulação da eq.67 ou 65 deve levar em consideração dois fatores e para termos uma idéia mais clara dos mesmos vamos analisar a eq.66 por partes.

Se eliminarmos o termo ¾kC da eq.66, a equação de difusão estará caracterizando essencialmente o movimento aleatório de partículas, comumente denominado de movimento browniano:

Na equação de Schrödinger, este termo corresponde a determinação da energia cinética do sistema. A solução analítica para esta equação diferencial é bem conhecida e corresponde a seguinte expressão³⁶:

Por outro lado, se desconsiderarmos o primeiro termo a direita da eq.66 teremos uma expressão que caracteriza um processo cinético de primeira ordem:

Em uma cinética de primeira ordem, dependendo do sinal e valor de k deveremos verificar um processo em que a concentração da substância aumenta ou diminui, ou seja, deveremos observar a criação ou aniquilação de componentes do sistema. Também podemos perceber que este termo que caracteriza uma cinética de primeira ordem está diretamente relacionado com a energia potencial na equação de Schrödinger (eq.67). A solução analítica para a eq.70 é conhecida e pode ser escrita como:

Além das duas informações acima, 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 estabelecermos uma simulação computacional para o processo de difusão, devemos procurar uma representação integral de equação de Schrödinger para que o método de Monte Carlo possa ser empregado de maneira satisfatória.

Uma forma compreensível de tal processo será apresentada em termos da equação de Schrödinger independente do tempo. Se desenvolvermos o inverso do operador hamiltoniano

^-1 e aplicarmos do lado direito e esquerdo da eq.1, teremos:

ou simplesmente:

Como o operador é um operador diferencial, o operador ^-1 deve ser um operador integral. Formalmente, a eq.73 pode ser então escrita como:

onde a integral do operador G(y,x) corresponde a representação do operador . Um maior aprofundamento nos aspectos matemáticos deste tipo de procedimento estão incorporados no campo das funções de Green³⁷. Para processos dependentes do tempo, como é o caso da eq.67, a eq.74 passa a ser escrita como:

Esta equação nos diz que em um determinado tempo t a chance de passarmos de uma configuração x para uma configuração y se dá através das características da função G(y,x;t).

Não podemos definir uma forma analítica para a função G(y,x;t) para sistemas complicados. Uma alternativa para resolvermos a eq.75 consiste em utilizar uma aproximação válida somente para o limite em que t tende a zero. Neste caso, a função de Green pode ser fatorada em um termo de difusão (G_dif) e um termo cinético (G_B), ou seja:

Estes dois termos estão diretamente associados as eqs.69 e 71, respectivamente, o que nos permite definir o termo de difusão por:

e o termo cinético por:

Na eq.78 a energia E_Tcorresponde ao valor médio da energia do sistema em estudo.

O desenvolvimento das eqs.77 e 78, associados ao uso de algoritmos eficientes para mapear todo o espaço das coordenadas eletrônicas permite simular o comportamento quântico do sistema. Uma vez que as informações que descrevem a natureza quântica de qualquer sistema microscópico é estatística, uma simulação da distribuição eletrônica deve ser definida através de valores médios de qualquer propriedade de interesse. Desta forma, ao procurarmos descrever o comportamento do elétron no espaço, deveríamos definir as coordenadas desse elétron em todas as possíveis posições do espaço (infinitas posições no espaço) ou então acompanhar a evolução do elétron em um tempo infinito. Quando os valores médios de qualquer propriedade eletrônica definida, ou em termos de um tempo infinito, ou em termos de um conjunto infinito de coordenadas forem os mesmos dizemos que temos uma integral ergódica. Este é o caso das integrais resolvidas através do método de Monte Carlo Quântico. Na prática, porém, ao invés de escolhermos qualquer uma das duas condições anteriores, definimos um número finito de posições no espaço e acompanhamos a evolução temporal das partículas nessas posições por um tempo suficiente para que a energia média do sistema permaneça constante. Matematicamente estamos dizendo que se escolhermos um conjunto de coordenadas arbitrárias no espaço para nosso elétron, {x}, poderemos determinar o valor médio da energia do sistema em relação a colocação de nosso elétron em cada uma das posições previamente definidas. Isto corresponderá a determinação do valor médio da energia em um tempo t qualquer. Assim, submetendo-se as posições do elétron a algum tipo de deslocamento, que depende essencialmente do tempo t, e observando-se o valor médio da energia em diferentes unidades de tempo, poderemos verificar o comportamento temporal da média espacial.

A extensão desta média temporal deverá ser escolhida de tal forma a produzir um valor estatístico adequado. Este valor estatístico poderá ser analisado, por exemplo, a partir da variância associada aos valores calculados ou através da análise de histogramas das diferentes energias obtidas em diferentes intervalos de tempo.

Um possível algoritmo de aplicação do Monte Carlo de Difusão ao átomo de hidrogênio encontrado na literatura pode ser resumido nas seguintes etapas³⁸:

1) Gera-se um número inicial de configurações de partículas (diferentes conjuntos de posições no espaço), que normalmente é da ordem de centenas. Na literatura estas configurações são denominadas de psips ou walkers. Estas configurações iniciais não podem ser muito diferentes da configuração de equilíbrio do átomo de hidrogênio e devem ser obtidas a partir de algum cálculo prévio. Normalmente, isto é feito definindo-se os pontos necessários para produzir resultados confiáveis de uma simulação de Monte Carlo Variacional de uma função de onda tentativa.

2) Escolhe-se um valor para E_T que deve ser uma aproximação ao valor da energia do átomo de hidrogênio.

3) Calcula-se o potencial eletrostático de uma configuração inicial (V(x)). Este potencial é simplesmente o potencial de Coulomb entre o psip, que é uma representação do elétron, e o núcleo do átomo de hidrogênio (eq.3). Este valor de energia potencial deve ser acumulado;

4) Dada esta configuração, o elétron é movimentado de forma que sua nova coordenada é gerada conforme:

onde c é um número aleatório pertencente à uma distribuição gaussiana com variância 2Dt. Em unidades atômicas o coeficiente de difusão D associado a equação de Schrödinger é igual a 1/2. O tempo t apresentado nesta expressão corresponde a um intervalo de tempo e não ao tempo total da simulação. Portanto, 2Dt=constante.

5) Calcula-se o potencial eletrostático da nova configuração (V(y)).

6) A função G_B está associada a criação e destruição de configurações eletrônicas. Desta forma, se uma configuração eletrônica obtida a partir do deslocamento eletrônico representar uma diminuição de energia potencial para o sistema, devemos multiplicar o número de configurações nesta região do espaço, caso contrário, deveremos destruir esta configuração eletrônica. Portanto, a partir de G_B obtemos o número de cópias desta configuração que serão geradas, M_B, de maneira que:

onde x é um número aleatório com distribuição uniforme. Caso M_Bseja igual a zero, esta nova configuração é rejeitada e mantém-se a configuração antiga;

7) Tomamos uma nova configuração e repete-se os passos 3 a 6 até que o número total de configurações seja atingido;

8) A partir das novas configurações, repete-se o procedimento até que o tempo total de simulação seja alcançado. Durante a simulação e periodicamente após um certo número de passos a energia tentativa E_T deve ser corrigida por:

em que representa a média da energia acumulada no passo 3;

9) Por fim, a energia do átomo de hidrogênio é dada por , que pode ser acompanhada durante a simulação.

Outras propriedades também podem ser calculadas de maneira muito simples pelo método de Monte Carlo de Difusão, bastando que seu valor seja acumulado durante a simulação e, ao término da mesma, obtemos uma média dos valores encontrados.

Devemos chamar a atenção de que estas simulações baseadas no método Monte Carlo de Difusão estão normalmente associadas ao estado fundamental do sistema. A justificativa para isto pode ser obtida considerando-se que a solução para a equação de Schrödinger dependente do tempo, que também pode ser escrita como:

pode ser representada como um somatório de um número infinito de estados Y_k, que são autofunções de , e que, portanto, levam a função de onda do sistema a ser representada como:

Podemos perceber, que para um tempo suficientemente longo, t'®¥, a função de onda tenderá a sua configuração de menor energia, Y₀.

Da mesma forma que discutido nos capítulos anteriores, o processo descrito aqui pode ser aplicado em sistemas multieletrônicos diversos. A simulação para tais sistemas pode apresentar problemas de ordem prática (convergência dos resultados) que são frequentemente minimizadas por técnicas de amostragem do espaço de configurações consideravelmente mais sofisticadas do que a discutida acima.

COMENTÁRIOS FINAIS

Como mencionado no início deste trabalho, o principal objetivo deste artigo foi o de sugerir uma reflexão profunda, principalmente para alunos iniciantes, sobre aspectos que consideramos muitas vezes como definitivos. Um exemplo tão simples e corriqueiro como o átomo de hidrogênio pode se transformar em um veículo extremamente poderoso para explorarmos novos modelos matemáticos, que eventualmente poderão nos transportar para a solução de problemas significativamente mais complexos. Os temas apresentados não se esgotam ao texto apresentado e uma sondagem cuidadosa na literatura especializada poderá demonstrar muitos outros aspectos tão curiosos quanto os discutidos neste artigo.

REFERÊNCIAS

• 1. Peixoto, E.M.A.; Quim. Nova 1978, 1, 5.
• 2. Peixoto, E.M.A.; Quim. Nova 1978, 2, 10.
• 3. Peixoto, E.M.A.; Quim. Nova 1978, 3, 13.
• 4. Peixoto, E.M.A.; Quim. Nova 1978, 4, 21.
• 5. Uma descrição sucinta dos postulados da mecânica quântica pode ser encontrada em http://www.chemkeys.com/bra/md/tdqq_4/pdmq_2/pdmq_2.htm
• 6. Ryzik, I.M.; Gradstein, I.S.; Tables of Integrals, Series and Products, Academic Press, New York, 1972.
• 7. Cunha, C.; Métodos numéricos para as engenharias e ciências aplicadas, Editora da Unicamp, São Paulo, 1993.
• 8. Atkinson, K. E.; "An Introduction to Numerical Analysis.", John Wiley & Sons, New York, 1978.
• 9. Hamming, R. W.; "Numerical methods for scientists and engineers.", McGraw-Hill, New York, 1973.
• 10. Blatt, J.M.; J. Comp. Phys 1967, 1, 382.
• 11. Blukis, U.; Howell, J.M.; J. Chem. Educ. 1983,60, 207.
• 12. Kubach, C.; J. Chem. Educ. 1983, 60, 212.
• 13. Hartree, D.R.; "The calculation of Atomic Structure", John Wiley & Sons, New York, 1957.
• 14. Roothaan, C.C.J.; Rev. Mod. Phys. 1951, 23, 69.
• 15. Szabo, A.; Ostlund, N.S.; Modern Quantum Chemistry, Dover Publications, New York, 1996.
• 16. Custodio, R.; Morgon, N.H.; de Andrade, J.C.; http://www.chemkeys.com/bra/md/tdqq_4/metlca_1/metlca_1.htm
• 17. Mulliken, R.S.; J.Chem.Phys. 1955, 23, 1833, 1841, 2338, 2343.
• 18. Raffenetty, R.C.; J. Chem. Phys. 1973, 58; 4452.
• 19. Feller, D.F.; Ruedenberg, K.; Theoret. Chim. Acta 1979, 20, 231.
• 20. Hill, D.L.; Wheeler, J.A.; Phys. Rev. 1953, 89, 1102.
• 21. Griffin, J.J.; Wheeler, J.A.; Phys. Rev. 1957, 108, 311.
• 22. Mohallem, J.R.; Dreizler, R.M.; Trsic, M.; Int. J. Quantum Chem. Symp. 1986, 20, 45.
• 23. da Costa, H.F.M.; da Silva, A.B.F.; Trsic, M.; Simas, A.M.; Aquino, A.J.A.; J.Mol. Struct. (THEOCHEM) 1990, 210, 63.
• 24. da Silva, A.B.F.; da Costa, H.F.M.; Trsic, M.; Mohallem, J.R.; Molec.Phys. 1989, 68, 433.
• 25. da Costa, H.F.M.; Simas, A.M.; Smith Jr., V.H.; Trsic, M.; Chem.Phys.Lett. 1992, 192, 195.
• 26. Custodio, R.; Giordan, M.; Morgon, N.H.; Goddard, J.D.; Int.J.Quantum Chem. 1992, 42, 411.
• 27. Custodio, R.; Goddard, J.D.; Giordan, M.; Morgon, N.H.; Can.J.Chem. 1992, 70, 580.
• 28. Custodio, R.; Goddard, J.D.; J.Mol.Struct.(THEOCHEM) 1992, 277, 263.
• 29. Vianna, R.O.; Custodio, R.; Chacham, H.; Mohallen, J.R.; Int.J.Quantum Chem.Symp. 1992, 26, 311.
• 30. Custodio, R.; Goddard, J.D.; J.Mol.Struct. (THEOCHEM) 1993, 281, 75.
• 31. Custodio, R.; Davis, W.M.; Goddard, J.D.; J.Mol.Struct. (THEOCHEM) 1994, 315, 163.
• 32. Morgon, N.H.; Custodio, R.; Riveros, J.M.; Chem.Phys.Lett. 1995, 235, 436.
• 33. Heermann, D.W.; "Computer Simulation Methods in Theoretical Physics", Springer-Verlag, 2^nd ed., New.York, 1990.
• 34. Kalos, M.H.; Whitlock, P.A.; "Monte Carlo Methods", John Wiley & Sons, 1^aed., vol 1, New York, 1986.
• 35. "Introduction to Monte Carlo Methods", Computacional Science Education Project 1995.
• 36. Hammond, B.L.; Lester Jr, W.A.; Reynolds, P.J.; Monte Carlo Methods in Ab Initio Quantum Chemistry, World Scientific, Singapore, 1994.
• 37. G.Barton; "Elements of Green's Functions and Propagation", Oxford University Press, 3^a ed, Oxford, 1995.
• 38. Lester Jr, W.A.; Recent Advances in Quantum Monte Carlo Methods, Chapter 6, World Scientific, Singapore, 1997.

Datas de Publicação

• Publicação nesta coleção
  31 Jul 2002
• Data do Fascículo
  Fev 2002

Histórico

• Aceito
  13 Jun 2001
• Recebido
  25 Abr 2001