SciELO - Scientific Electronic Library Online

 
vol.35 issue10Acid-base reactions: concept, representation and generalization from the energy involved in transformationsChemical defense: history, warfare agent classification and action of neurotoxic agents author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

Share


Química Nova

Print version ISSN 0100-4042

Quím. Nova vol.35 no.10 São Paulo  2012

https://doi.org/10.1590/S0100-40422012001000032 

EDUCAÇÃO

 

Resolvendo a equação de Schrödinger utilizando procedimentos numéricos fundamentais

 

Solving the Schrödinger equation using fundamental numerical procedures

 

 

Rogério Custodio*; Maurício Ramalho Custodio; Eduardo José Creatto

Instituto de Química, Universidade Estadual de Campinas, CP 6154, 13083-970 Campinas - São Paulo, Brasil

 

 


ABSTRACT

A combination of the variational principle, expectation value and Quantum Monte Carlo method is used to solve the Schrödinger equation for some simple systems. The results are accurate and the simplicity of this version of the Variational Quantum Monte Carlo method provides a powerful tool to teach alternative procedures and fundamental concepts in quantum chemistry courses. Some numerical procedures are described in order to control accuracy and computational efficiency. The method was applied to the ground state energies and a first attempt to obtain excited states is described.

Keywords: Variational Monte Carlo; expectation value; fundamental and excited states.


 

 

INTRODUÇÃO

O método Monte Carlo corresponde a uma das técnicas numéricas mais poderosas ao alcance de cientistas, engenheiros e não especialistas. Essa popularidade está associada à facilidade, eficiência e precisão com que o método pode ser empregado. Seu espectro de aplicações é vasto e uma busca na literatura ou em bases de dados poderá dar um exemplo dos milhares de artigos e livros publicados em diferentes níveis de complexidade. Para uma abordagem didática do método o leitor pode consultar as refs. 1-8.

A aplicação do método Monte Carlo para resolver a equação de Schrödinger dependente ou independente do tempo recebeu a denominação genérica de Monte Carlo Quântico (MCQ) e uma abordagem simplificada, explorando as potencialidades do método para a solução de problemas de estrutura eletrônica, foi publicada em artigo recente.6 O texto concentrou-se na descrição de duas das mais utilizadas aproximações do MCQ: o Monte Carlo Variacional (MCV) e o Monte Carlo de Difusão (MCD), além de ilustrar com aplicações de altíssima precisão as características próprias de cada método, ressaltando vantagens e desvantagens. Tanto a leitura do artigo quanto de trabalhos especializados da literatura usualmente explicita o melhor desempenho em termos de exatidão do MCD em relação ao MCV. O MCV é visto como um método mais limitado e restrito à flexibilidade da função de onda utilizada, não possibilitando, como ocorre com o MCD, produzir resultados com precisão além da função de onda utilizada como guia da simulação.

Essa perspectiva vem sendo modificada com avanços recentes. Atualmente, é possível utilizar funções de onda altamente flexíveis e métodos de otimização eficientes com o MCV.9-11 Porém, o MCV pode ir além do simples ajuste de parâmetros em funções de onda flexíveis. Esta nova alternativa foi sugerida apenas como uma curiosidade de aplicação didática do princípio variacional e uso do método Monte Carlo na ref. 4. Realizando testes com o MCV alternativo (MCVA) verificou-se que o mesmo é computacionalmente simples, mas ineficiente quando aplicado em sistemas mais complexos, porém, permite uma compreensão clara de conceitos frequentemente utilizados em disciplinas de química quântica e mal compreendidos por alunos. O desenvolvimento numérico e o esforço computacional para atingir a precisão desejada levam ao questionamento da conveniência de procedimentos considerados como mandatórios para a resolução da equação de Schrödinger.

O objetivo deste artigo é apresentar esta nova versão do MCV como instrumento didático, explorando sua aplicabilidade em sistemas simples de forma não ortodoxa, permitindo uma visualização das necessidades matemáticas de alguns métodos usualmente abordados em disciplinas de química quântica.

 

PROCEDIMENTO NUMÉRICO

O método MCV considera o ajuste de uma função de onda usando números aleatórios e minimizando a energia do estado fundamental de um sistema. A essência do princípio variacional está em garantir que qualquer função de onda tentativa Ψ sempre produzirá uma energia para o estado fundamental, que será sempre maior ou igual ao valor exato. A determinação da energia de um sistema quântico é feita através da expressão do teorema do valor médio definido como:

 

sendo Ĥ o operador hamiltoniano e Ψ uma função de onda tentativa.

Usualmente, quando se entra em contato com o princípio variacional pergunta-se o que exatamente pode ser utilizado como uma função de onda tentativa. Através dos postulados da mecânica quântica verifica-se que a única exigência para se elaborar uma função de onda é de que ela seja bem comportada, ou seja, finita, unívoca e contínua, além da exigência de que a primeira e segunda derivadas também obedeçam a tais critérios.12 Desta forma, a resposta para o que se pode usar como função de onda na Equação1 é literalmente qualquer função que preserve o caráter bem comportado. O que normalmente não se considera é a aplicação do princípio variacional para gerar a função de onda em si. Até onde sabemos, esta alternativa só foi indicada como exemplo didático por Giordano e Nakanishi.4 O método permite resolver a equação de Schrödinger numericamente, sem necessidade de recorrer às técnicas de separação de variáveis ou transformação de coordenadas ou modelo de partículas independentes, embora o leitor iniciante deva compreender que o uso desses recursos, além de facilitar a descrição do sistema, torna a solução analítica ou numérica mais simples e eficiente.

Para descrever o processo, a Equação1 será reescrita em uma dimensão como:

 

sendo que o operador hamiltoniano foi substituído pelo operador de energia cinética e o operador de energia potencial, respectivamente. Numericamente, a derivada segunda da energia cinética pode ser substituída por uma aproximação facilmente deduzida de séries de Taylor:13,14

  

O cálculo da derivada segunda depende do valor da função de onda em uma coordenada x e em duas posições vizinhas separadas por ±x. Substituindo-se a Equação3 na Equação 2, têm-se:

 

A Equação 4 pode ser ainda aproximada empregando-se, por exemplo, um processo de integração trapezoidal, ou seja, a área total a ser integrada é subdividida em vários trapézios e suas áreas são determinadas e somadas definindo aproximadamente a área total.13,14 Em termos práticos, substitui-se o elemento infinitesimal dx por um elemento finito x' e determina-se a área associada a cada elemento finito. Este elemento finito pode ser igual ou diferente daquele usado na derivada segunda dada pela Equação 4. Por conveniência, considera-se que os valores de x são igualmente espaçados e que os elementos finitos das integrais e derivada segunda são todos iguais: x'=x. Desta forma, o cálculo do valor médio da energia pode ser estimado para qualquer função de onda considerando-se a Equação 5:

 

Para se determinar E é necessário definir quais são as coordenadas {x} a serem utilizadas e quais são os valores da função de onda nos respectivos valores de x. No limite quando x0, os somatórios na Equação5 tendem aos valores exatos. Assim, a escolha de x será definida pela precisão numérica a ser atingida e a disponibilidade e eficiência computacional.

A construção da função de onda tentativa é definida pelo princípio variacional. Estabelece-se um valor arbitrário para cada Ψ(xi) e calcula-se a energia através da Equação 5. Sorteia-se uma determinada coordenada xi e modifica-se o valor da função de onda nessa coordenada por um número aleatório como, por exemplo: Ψ(xi) = Ψ(xi) + (0,5 - rand).d, sendo rand um número aleatório com distribuição uniforme entre 0 e 1 e δ um parâmetro que define o intervalo arbitrário de variação da função de onda. Se δ = 1, a função de onda oscilará entre os valores de ± 0,5. A energia para a função de onda modificada é calculada (Equação 5) e se o valor for menor, aceita-se a modificação na função de onda. Se a energia aumentar, preserva-se o valor original da função de onda. Este processo é repetido até que a energia do sistema atinja um limite mínimo desejado.

É conveniente mencionar que existe o método de diferenças finitas, utilizado em cálculos numéricos da equação de Schrödinger de sistemas simples.13,14 Neste método, resolve-se a equação diferencial e não a equação integral. A vantagem do presente método é que está fundamentado no teorema variacional e, em princípio, poderia ser aplicado para qualquer sistema. No método de diferenças finitas a solução é obtida quando o "chute" dado para a energia do sistema produz uma função de onda bem comportada.

 

A PARTÍCULA NA CAIXA 1D E 2D E A SEPARAÇÃO DE VARIÁVEIS

Serão considerados neste item dois exemplos bem conhecidos: a partícula na caixa unidimensional (1D) e a partícula na caixa bidimensional (2D) com extremidades com potencial infinito. A partícula na caixa 1D será definida considerando-se uma caixa com dimensões de 1,0 u.a. de comprimento. As coordenadas inicial e final dos limites da caixa foram: xi = -0,5 u.a. e xf = 0,5 u.a. Nestas extremidades as condições de contorno determinam que Ψ(x) = 0,0 para x -0,5 u.a. e x +0,5 u.a. Para valores de x dentro da caixa, o potencial é nulo e a função de onda e a energia do sistema devem ser determinadas. As soluções analíticas para este modelo são dadas pelas seguintes expressões: e , sendo nx um número inteiro maior do que zero, a o comprimento da caixa, m a massa da partícula e h a constante de Planck. Assim, o valor exato para a energia de um elétron confinado em um fio com comprimento de 1 u.a. para nx = 1 é igual a Enx = 4,93480 u.a.

Uma vez que o potencial é nulo dentro da caixa, a Equação 5 em unidades atômicas é simplificada para:

 

A Figura 1 apresenta um diagrama da função de onda inicial, Ψinicial, e a função de onda final, Ψfinal, para uma simulação empregando 20 valores de x e minimizando-se a Equação 6. Os valores da função de onda nas extremidades foram mantidos em zero, satisfazendo as condições de contorno, durante toda simulação. A função de onda inicial foi arbitrariamente definida como tendo valores iguais a 1,00 para 14 pontos e zero para 4 pontos (Figura 1), sendo normalizada durante toda simulação pela Equação 7:

 

 

A simulação continuou por um número de passos suficiente para atingir um nível de precisão em energia de 10-5 u.a. Na Tabela 1 são apresentados valores de energias para simulações com diferentes números de pontos para a partícula na caixa 1D (segunda coluna) e para a partícula na caixa 2D (terceira coluna), cujos resultados serão discutidos posteriormente. A primeira e óbvia conclusão é que a tendência ao valor exato ocorre com o aumento do número de pontos utilizado. Pode-se perceber que, embora a precisão alcançada seja excelente, para se atingir o resultado exato, é necessário um número significativo de pontos.

 

 

A razão entre o número de modificações na função de onda que proporcionam uma energia menor dividido pelo número total de modificações define o que se conhece como taxa de aceitação. Em cálculos Monte Carlo,6 a simulação é controlada para fornecer uma taxa de aceitação em torno de 0,5. No presente método, a taxa de aceitação é ajustada a partir do valor de δ utilizado para modificar Ψ(x). Entretanto, à medida que a simulação transcorre, a função de onda se aproxima da solução exata e se o valor de δ for fixo, a taxa de aceitação tende a zero. Nossos testes indicam que há necessidade de se modificar δ quando a taxa apresentar valores inferiores a 0,01 ou 1%. Para uma maior eficiência, esta redução na ordem de grandeza em δ pode ser feita dividindo-se este parâmetro por um número menor do que 10, cada vez que a taxa de aceitação atinge um valor inferior a 1% a cada 10.000 passos. Também é conveniente repetir-se o cálculo para um processo que já tenha alcançado taxa de aceitação menor do que 1%, uma vez que a amostragem do espaço é aleatória e pode ocorrer do espaço amostrado ter maior ou menor taxa de aceitação.

Usualmente, em disciplinas introdutórias de quântica após o modelo da partícula na caixa 1D é apresentado o tratamento para a partícula na caixa 2D ou 3D para ilustrar o uso do recurso de separação de variáveis. Com o MCVA é possível resolver a equação de Schrödinger para a partícula na caixa 2D sem utilizar a separação de variáveis. A função de onda 2D será dependente das coordenadas (x,y), ou seja, Ψ = Ψ(x,y). A partícula em uma caixa 2D estará livre dentro de uma área definida no intervalo entre -0,5 u.a. x 0,5 u.a. e -0,5 u.a. y 0,5 u.a. Fora desses limites o potencial é infinito e a função de onda será nula. O operador hamiltoniano para a partícula na caixa 2D é exclusivamente o operador de energia cinética. Substituindo-se o operador de energia cinética para duas dimensões na Equação 1 e rearranjando, têm-se:

A Equação 8 pode ser aproximada por somatórios em cada coordenada e a derivada segunda sobre a função de onda 2D pode ser descrita pela Equação 3, lembrando-se que se tem agora derivadas parciais, que estão sendo aplicadas em uma coordenada de cada vez, o que leva à Equação 9:

 

Nesta equação, x e y são constantes. Para se utilizar o MCVA deve-se sortear pares de coordenadas (xi,yi), modificar o respectivo valor da função de onda e calcular o valor de energia média com a Equação 9. A aplicação desse procedimento foi realizada com uma função de onda inicial semelhante àquela apresentada na Figura1, mas agora considerando uma superfície constante no plano (x,y) e o valor igual a zero nos limites e fora da caixa para satisfazer as condições de contorno. Os resultados para simulações com diferentes números de pontos encontram-se na coluna 2D da Tabela 1. Pode-se verificar que quanto maior o número de pontos, mais próximo está o resultado da solução exata. Além desse aspecto óbvio já exemplificado na caixa 1D, pode-se constatar que as energias para a caixa 2D são exatamente iguais a duas vezes as respectivas energias 1D. Se a separação de variáveis fosse considerada, esse seria o resultado correto. Portanto, a avaliação numérica demonstra que o processo de separação de variáveis, além de permitir a solução analítica das equações diferenciais para este sistema, é rigorosamente correto e numericamente poderia ser realizado em um tempo computacional muito menor.

Existiria alguma vantagem em se resolver o sistema 2D sem lançar mão do método de separação de variáveis, mesmo que seja uma solução numérica? Para este sistema em particular a resposta é: não! Nada substitui o bom senso ou as soluções analíticas. Mas, o objetivo de uma solução numérica não é necessariamente produzir números, mas informações. O esforço computacional utilizado para resolver a partícula na caixa 2D é significativamente maior do que para 1D. O tempo gasto para a resolução MCVA da caixa 2D demonstra que o método é factível, mas converge muito lentamente para um número grande de pontos. Para contornar este problema e reduzir o tempo de CPU, determinou-se a função de onda ótima para uma malha bruta, por exemplo, 10 x 10 e utilizou-se esta função de onda para determinar valores aproximados de funções de onda para uma malha mais fina, por exemplo, 12 x 12. A partir da otimização desta malha, interpolou-se valores para a malha 14 x 14 e assim sucessivamente. Em nossos cálculos considerou-se a interpolação linear através do método de diferença dividida de Newton13 utilizando a Equação 10:

 

sendo xi+1, xi, yi+1 e yi valores conhecidos da malha bruta e xi+1 > x > xi e yi+1 > y > yi, os valores das coordenadas onde se pretende estimar o valor de Ψ(x,y).

O tempo de cálculo para a caixa 2D, se o sistema for construído de forma semelhante ao exemplo da caixa 1D, deverá ser da ordem de t2, sendo t o tempo para completar o cálculo para a caixa 1D, ou seja, o aumento no número de dimensões acarreta necessariamente em um aumento significativo no tempo de CPU.

O procedimento a ser empregado para caixas com potencial interno é exatamente o mesmo. A única diferença para esses casos é que no uso das Equações 5 ou 8 tem-se que calcular também o termo potencial.

 

OSCILADORES HARMÔNICOS E ANARMÔNICOS

Na aplicação do método para o oscilador harmônico ou anarmônico nos deparamos com a necessidade de definir os limites de integração de ± no caso de uso de coordenadas cartesianas. Essas integrais são denominadas de integrais impróprias e existem diferentes abordagens para resolvê-las numericamente. A abordagem tecnicamente mais correta corresponde a reescrever as integrais da Equação 1, dividindo-se cada uma em duas ou mais componentes. Para uma integral unidimensional cujos limites de integração estão entre ± pode-se dividi-la arbitrariamente em três integrais:

 

Sendo a função f(x) convergente, a integral central com limites entre a e b representa a componente mais importante. As outras duas integrais são denominadas de "caudas" da integral. A escolha dos limites a e b deve ser feita de tal maneira que o efeito das caudas seja pequeno. A maneira mais simples para se calcular os efeitos de cauda é resolver estas integrais em espaçamento não constante. A integral central poderá ser dividida em espaçamentos constantes e as integrais de cauda terão os respectivos valores de x definidos por uma transformação de coordenadas como, por exemplo, x = 1/s. Nesta transformação os limites para, por exemplo, a primeira cauda passam de - a a para valores finitos de 0 a 1/a. A escolha de um intervalo s constante permite a determinação de valores de si e, consequentemente, dos valores de xi correspondentes. O espaçamento em x deixa de ser constante, mas a discretização dada pela Equação 5 pode ser utilizada, desde que se considere para cada ponto xi um xi correspondente. A derivada segunda no operador de energia cinética também deve ser modificada, sendo definida neste trabalho pelo método de diferença dividida de Newton,13 ou seja:

 

Quando xi+2 - xi+1 = xi+1 - xi tem-se que xi+2 - xi = 2(xi+1 - xi) que, substituído na Equação 12 leva à Equação 3.

Para ilustrar o uso desta técnica de integração no MCVA, a equação de Schrödinger será resolvida numericamente para o oscilador harmônico e o potencial de Morse: V(x) = 0,5x2 e , respectivamente.

As integrações para o oscilador harmônico foram feitas considerando-se diferentes malhas com diferentes limites de integração. Foram escolhidos três limites de integração internos (a e b): ±1 u.a., ±2 u.a. e ±5 u.a. No limite interno entre ±1 u.a. o efeito de cauda não é desprezível e no limite de ±5 u.a. é praticamente nulo. Cada integral de cauda foi dividida em malhas de 10 pontos e a função de onda nos dois últimos pontos quando o limite tende a infinito foi mantida igual a zero durante todo o processo de integração. Na Tabela 2 encontram-se valores destas simulações.

 

 

O primeiro aspecto que chama a atenção é que no limite interno de ±1 u.a. com 10 pontos de integração, o valor da energia difere em apenas 0,00052 u.a. do valor exato. O mesmo cálculo com o mesmo número de pontos, mas com um limite interno de ±2 u.a. apresenta um desvio de 0,00656 u.a. em relação ao valor exato e para o limite interno de ±5 u.a. o desvio vai para 0,03733 u.a. Esse aumento no desvio da energia para os limites de integração maiores são consequência do maior espaçamento em x. Para os limites de integração de ±1 u.a., ±2 u.a. e ±5 u.a. os intervalos x são de 0,22222 u.a., 0,44444 u.a. e 1,11111 u.a., respectivamente. Em outras palavras, a densidade de pontos para o limite de ±1 u.a. é significativamente maior e, portanto, descreve a região integrada muito melhor do que ocorre com os limites de ±2 u.a. e ±5 u.a. Para que as integrais envolvendo os limites maiores apresentem a mesma densidade de pontos, o resultado para 10 pontos com o limite de integração de ±1 u.a. deve ser comparado com resultados para 18 e 45 pontos nos limites de integração de ±2 u.a. e ±5 u.a., respectivamente. Se for feita esta comparação, verifica-se ainda que o resultado para o limite de ±1 u.a. se encontra mais próximo do valor exato do que os outros dois cálculos envolvendo limites de integração maiores. A resposta para esta tendência é uma combinação de: a) precisão da região descrita, b) precisão nas integrais de cauda e c) possível cancelamento de erros dos três cálculos. Se o comportamento das três simulações for analisado à medida que o número de pontos aumenta, verifica-se que as simulações com maior limite de integração se aproximam sistematicamente do valor exato. Para as outras duas integrais com limites de convergência menores, os resultados não necessariamente convergem para o resultado exato, porque as integrais de cauda apresentam um papel não desprezível. A maior ou menor aproximação em relação ao resultado exato depende se a discretização representa de maneira mais apropriada ou não a região de cauda. Isto pode ser feito aumentando-se o número de pontos nessa região ou escolhendo-se um método de integração numérico mais elaborado. De maneira geral, pode-se perceber que a escolha conveniente dos limites de integração pode levar a excelentes resultados se a malha na integral central for densa e se as integrações das integrais de cauda descreverem apropriadamente as regiões, tendendo ao limite apropriado.

Da mesma forma que para a partícula na caixa, o oscilador harmônico pode ser resolvido analiticamente. Efeitos de anarmonicidade normalmente são considerados definindo-se potenciais sofisticados. A resolução da equação de Schrödinger com esses potenciais não é trivial e é usualmente efetuada numericamente. Para testar a eficiência do método MCVA, considerou-se a resolução da equação de Schrödinger utilizando-se o potencial de Morse. A Tabela 3 apresenta o resultado da energia do ponto zero estimado experimentalmente para uma série de moléculas diatômicas e resultados calculados com o MCVA. A Tabela 3 também apresenta os parâmetros utilizados para representar o potencial de Morse para cada diatômica. Os limites de integração da integral central foram escolhidos para que as integrais de cauda tivessem valores desprezíveis. Entretanto, mantiveram-se as caudas das integrais com 10 pontos de integração e a integral de cauda à esquerda foi avaliada sem necessidade de transformação de coordenadas, uma vez que o intervalo é finito (de 0 ao limite inferior da integral central). As simulações foram feitas com 200 pontos e a convergência na energia com uma precisão da ordem de 10-6 u.a. Os resultados demonstram uma precisão considerável entre os valores calculados e os estimados experimentalmente. O maior desvio encontrado foi da ordem de 3 cm-1 (~1 x 10-5 u.a.), que pode ser considerado extremamente preciso em relação a outros cálculos teóricos. O nível de dificuldade foi limitado apenas ao tempo de CPU para que a energia atingisse o limite desejado.

Embora estes resultados sejam estimulantes em termos de simplicidade e nível de precisão, o cálculo de espectros infravermelho para diatômicas depende da estimativa de estados excitados. O método discutido até o presente foi sugerido apenas para o estado fundamental. Pode-se verificar que as exigências para o cálculo de estados excitados podem ser consideradas, incluindo-se um mínimo de conceitos fundamentais discutidos em disciplinas de quântica.

Estados excitados

A única exigência adicional necessária para calcular estados excitados é a ortogonalidade entre esses estados. Para introduzir o tema de maneira simplificada, será apresentada a possibilidade de determinação de um único estado excitado. Considere que o estado fundamental será descrito por uma função de onda Ψa e o primeiro estado excitado, por uma função de onda Ψb. Considere ainda que estas duas funções de onda não são ortogonais, mas são normalizadas. Em outras palavras, a integral de recobrimento entre esses dois estados não é nula: Sab = Ψ*aΨbdx 0.

Existem diferentes maneiras de se ortogonalizar uma função de onda. Neste caso, será utilizado o método de ortogonalização de Schmidt.16,17 Considere que a função do estado fundamental e excitado ortogonais são ΨA e ΨB, respectivamente. O método de Schmidt considera que ΨA = Ψa e que ΨB = N(Ψb + λΨa), sendo N um fator de normalização e l um fator que determinará a condição de ortogonalidade entre os dois estados. A definição destes dois parâmetros é feita através das integrais de normalização (∫Ψ*BΨBdx = 1 = ∫Ψ*AΨAdx) e ortogonalização (∫Ψ*AΨBdx = 0), levando .

Para determinar a função de onda ΨB e a energia desse estado excitado, é necessário que se conheça a função de onda do estado fundamental ou, então, que ambas sejam geradas na mesma simulação. A energia do estado A deve ser obtida pela Equação 5 e a energia do estado B será definida por:

A última integral à direita corresponde à energia do estado A, a primeira integral à direita da igualdade é a energia da função de onda não ortogonal do estado B que está sendo modificada durante a simulação, para produzir o mínimo de energia desse estado, e a integral intermediária corresponde à energia de interação das funções de onda não ortogonais. Estas três integrais serão discretizadas como mostrado anteriormente. As funções de onda Ψa e Ψb devem ser sempre normalizadas (Equação 7) antes do cálculo das respectivas energias. Isto sendo feito, a discretização dada pela Equação 5 não necessita das integrais nos denominadores. A integral de recobrimento Sab é discretizada da mesma forma como todas as outras integrais. Finalmente, uma vez programada a Equação 13, basta aplicar o princípio variacional modificando aleatoriamente a função de onda Ψb para a determinação da energia e função de onda do primeiro estado excitado.

A Figura 2a ilustra as funções de onda Ψa e Ψb iniciais antes da normalização para a determinação simultânea do estado fundamental e primeiro estado excitado do oscilador harmônico e a Figura 2b mostra as funções ortogonais e otimizadas para os dois estados. Nota-se que a função de onda inicial do estado B foi propositalmente construída com um nó para facilitar a convergência da simulação. Funções de onda construídas com números aleatórios podem ser consideradas, mas consumirão um tempo de CPU superior para convergir. Nesta simulação simultânea do estado fundamental e do primeiro estado excitado à respectiva energia de cada estado com limites de integração de ±5 u.a., malhas de 100 pontos com caudas contendo mais 10 pontos proporcionaram energias iguais a 0,49968 e 1,49855 u.a. O valor exato para estes dois estados corresponde a 0,50 e 1,50 u.a. Como nas aplicações anteriores, este exemplo deixa claro que a combinação do princípio variacional e do método de Monte Carlo em uma versão de malha fixa permite resolver a equação de Schrödinger não apenas para o estado fundamental como também para estados excitados. Para considerar outros estados excitados, basta apenas que as respectivas funções de onda mantenham a ortogonalidade em relação a todos os estados já descritos.

 


 

Os exemplos apresentados neste artigo proporcionam um método numérico alternativo, não ortodoxo, possibilitando uma visão da relevância de alguns conceitos fundamentais em química quântica e a conveniência de outros recursos usualmente apresentados nesta disciplina como mandatórios para a resolução da equação de Schrödinger. As perspectivas de aplicação do método em sistemas mais complexos proporcionarão uma visão alternativa de problemas fundamentais em disciplinas similares. Em artigo publicado anteriormente foram apresentadas e discutidas quatro alternativas diferentes para resolver a equação de Schrödinger para o átomo de hidrogênio.18 Com o presente método é possível antecipar que existe, portanto, uma quinta maneira para resolver a equação de Schrödinger para esse sistema. Nesta quinta alternativa, o tratamento numérico levará naturalmente à necessidade de se resolver um problema que usualmente não é percebido nas soluções anteriores. Por meio da necessidade de se definir uma malha de pontos no processo de integração numérica, verifica-se que dependendo da escolha da malha, surge naturalmente um problema denominado de correlação elétron-núcleo, conhecido também pelo nome de cúspide nuclear.19 Da mesma forma, a utilização do mesmo procedimento para sistemas com um, dois ou mais elétrons suscitará o questionamento das necessidades e conveniências de procedimentos, tais como, a aproximação de Born-Oppenheimer, o modelo de partículas independentes de Hartree ou Hartree-Fock,17 ou como o cálculo de propriedades atômicas ou moleculares, como as forças de Hellmann-Feynman,20,21 pode ser afetado com as malhas utilizadas, embora as energias dos sistemas sejam razoavelmente bem descritas, etc. Trabalhos nestas questões fundamentais estão em andamento. Gostaríamos de ressaltar que grande parte dos cálculos realizados neste artigo foi efetuada por dois alunos de iniciação científica utilizando laptops convencionais em linguagens de programação de escolha pessoal. Interessados em exemplo de programa podem contatar o autor para correspondência.

 

CONCLUSÃO

A combinação do princípio variacional com o teorema do valor médio e o método de Monte Carlo corresponde a uma ferramenta didática e acadêmica extremamente poderosa, no sentido de possibilitar uma avaliação criteriosa da relevância de técnicas usualmente consideradas imperativas na solução da equação de Schrödinger. Neste artigo esta combinação de métodos, que denominamos de MCVA, foi aplicada em sistemas simples, tais como, partícula na caixa 1D e 2D e oscilador harmônico e anarmônico representado pelo potencial de Morse. Os três primeiros sistemas foram considerados para apresentar de forma didática a simplicidade do método, algumas dificuldades e alternativas numéricas e precisão possível de ser alcançada. O oscilador anarmônico foi apresentado como um exemplo de aplicação do método em um sistema usualmente considerado como de solução não trivial. Com a solução para o estado fundamental, o artigo esboçou as condições necessárias para o cálculo de estados excitados empregando os mesmos princípios e acrescentando a condição de ortogonalidade entre os estados. O mesmo alto nível de precisão foi demonstrado para o cálculo do primeiro estado excitado para o exemplo do oscilador harmônico.

Além da precisão, a apresentação do método tem o objetivo de suscitar o questionamento por parte de alunos e usuários de métodos computacionais sobre alguns aspectos frequentemente utilizados na resolução de problemas em química quântica.

O novo Monte Carlo variacional, como todo método, possui limitações que serão discutidas e aprofundadas futuramente no tratamento de sistemas eletrônicos. De maneira geral, pela facilidade de uso para sistemas simples e caráter não ortodoxo na resolução da equação de Schrödinger, a combinação dos procedimentos é extremamente atrativa e reforça fundamentos e métodos usados em trabalhos de química computacional.

 

AGRADECIMENTOS

Ao apoio financeiro das agências de fomento FAPESP, CNPq e FAEPEX - UNICAMP. A disponibilidade computacional do Centro Nacional de Processamento de Alto Desempenho - CENAPAD, Campinas, também está auxiliando no desenvolvimento deste projeto, ao qual somos gratos.

 

REFERÊNCIAS

1. DeVries, P. L.; A First Course in Computational Physics, John Wiley & Sons: New York, 1994.         [ Links ]

2. Thijssen, J. M.; Computational Physics, Cambridge University Press: Cambridge, 2003.         [ Links ]

3. Gould, H.; Tobochnik, J.; Christian, W.; An Introduction to Computer Simulation Methods, 3rd ed., Pearson - Addison Wesley: New Jersey, 2007.         [ Links ]

4. Giordano, N. J.; Nakanishi, H.; Computational Physics, 2nd ed., Pearson - Addison Wesley: New Jersey, 2007.         [ Links ]

5. Morgon, N. H.; Coutinho, K., eds.; Métodos de Química Teórica e Modelagem Molecular, Livraria da Física Ed.: São Paulo, 2007.         [ Links ]

6. Angelotti, W. F. D.; Fonseca, A. L.; Torres, G. B.; Custodio, R.; Quim. Nova 2008, 31, 433.         [ Links ]

7. Novo, J. B. M.; Dias, L. C.; Quim. Nova 2011, 34, 707.         [ Links ]

8. Lopez-Castillo, A.; de Souza Filho, J. C.; Quim. Nova 2007, 30, 1759.         [ Links ]

9. Bouabça, T.; Braïda, B.; Caffarel, M.; J. Chem. Phys. 2010, 133, 044111.         [ Links ]

10. Toulouse, J.; Umrigar, C. J.; J. Chem. Phys. 2008, 128, 174101.         [ Links ]

11. Umrigar, C. J.; Toulouse, J.; Filippi, C.; Sorella, S.; Hennig, R. G.; Phys. Rev. Lett. 2007, 98, 110201.         [ Links ]

12. Peixoto, E. M. A.; Quim. Nova 1978, 2, 10.         [ Links ]

13. Carnahan, B.; Luther, H. A.; Wilkes, J. O.; Applied Numerical Methods, John Wiley & Sons, Inc.: New York, 1969.         [ Links ]

14. Acton, F. S.; Numerical Methods that Usually Work, Mathematical Association of America: Washington, 1990.         [ Links ]

15. Huber, K. P.; Herzberg, G.; Molecular Spectra and Molecular Structure. IV. Constants of diatomic molecules, van Nostrand Reinhold: New York, 1979.         [ Links ]

16. Custodio, R.; Quim. Nova 1987, 9, 89.         [ Links ]

17. Szabo, A.; Ostlund, N. S.; Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Dover: New York, 1996.         [ Links ]

18. Custodio, R.; Politi, J. R. S.; Segala, M.; Haiduke, R. L. A.; Cyrillo, M.; Quim. Nova 2002, 25, 159.         [ Links ]

19. Hättig, C.; Klopper, W.; Köhn, A.; Tew, D. P.; Chem. Rev. 2012, 112, 4.         [ Links ]

20. Deb, B. M.; Rev. Mod. Phys. 1973, 45, 22.         [ Links ]

21. Mohallem, J. R.; Custodio, R.; Chacham, H.; Vianna, R.O.; Int. J. Quantum Chem. Sup. 1992, 26, 311.         [ Links ]

 

 

Recebido em 27/2/12
Aceito em 31/5/12
Publicado na web em 31/8/12

 

 

* e-mail: roger@iqm.unicamp.br

Creative Commons License All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License