devido ao potencial de interação com o j-ésimo elétron. O i-ésimo elétron deverá se mover em um potencial médio Vi(~r) = Z d3r′ 1 |~r − ~r′| X j(6=i) |ϕj(~r′)|2− X m Z |~r − ~Rm| . (2.19)
Para encontrar o orbital ϕi(~r) do i-ésimo elétron devemos resolver a equação de Schrödinger
para uma partícula simples −∇ 2 2m + Vi(~r) ϕi(~r) = εiϕi(r). (2.20)
Acima temos um conjunto de equações que devem ser solucionadas numericamente de forma iterativa. As equações 2.19 e 2.20 são conhecidas como equações de Hartree e podem ser derivadas do princípio variacional usando 2.17 como função tentativa. A função de estado descrita na equação 2.17 não é antissimétrica, de forma que não satisfaz o princípio da exclusão de Pauli.
Para obter funções de estado antissimétricas utilizamos a seguinte expressão:
Ψ(~r1, ~r2, ..., ~rN) = 1 √ N ! ψ1(~r1) ψ2(~r1) . . . ψN(~r1) ψ1(~r2) ψ2(~r2) . . . ψN(~r2) ... ... ... ψ1(~rN) ψ2(~rN) . . . ψN(~rN) . (2.21)
O termo ~ri denota as coordenadas espaciais e de spin da partícula, ψi(~rj) é chamado de spin-
orbital e é constituído pelo produto do orbital espacial φi(~rj) com as funções de spin αi(σj)
correspondentes ao spin “up” ou “down”. Os spin-orbitais formam um conjunto ortonormal, de modo que
Z d~rψ∗
i(~r)ψj(~r) = δij. (2.22)
A equação 2.21 é conhecida como determinante de Slater e o fator 1/√N ! é devido a con- dição de normalização da função de estado. Esta forma de obter a função de estado de um sistema constituído de férmions foi sugerida por Heisenberg e Dirac, e aplicada ao problema multieletrônico por Slater [92]. A utilização do determinante de Slater como função tentativa
no método variacional, leva a um novo conjunto de equações conhecidas como equações de Hartree-Fock: −∇ 2 2m + Vi(~r) ψi(~r) − X j ψj(~r) Z d~r′ψ ∗ i(~r′)ψj(~r′) |~r − ~r′| = εiψi(r). (2.23)
O potencial obtido na aproximação de Hartree é expresso através do termo Vi(~r). O terceiro
termo à esquerda da igualdade acima é chamado de potencial de troca e correlação e é uma consequência da antissimetria da função de onda. Tal termo indica a correlação entre elétrons de spins paralelos. A implementação deste método utiliza uma técnica devida a Roothaan, onde expande-se os orbitais espaciais como uma combinação de um conjunto de funções de base de elétron, ou seja, uma combinação linear de orbitais atômicos (LCAO- Linear Combination of Atomic Orbitals).
O determinante de Slater leva a melhor solução para o problema de N-elétrons, mas não à solução exata. Uma parte dos efeitos de correlação entre os elétrons é perdido devido a variação do potencial colombiano exato visto por um elétron quando os outros se movimentam. Isto porque, o determinante de Slater leva em conta o fator de troca e correlação entre dois elétrons paralelos, mas o potencial coulombiano exato sentido por um elétron varia quando os outros se deslocam. Como consequência disto, temos um erro de correlação que é solucionado usando nessa teoria a função de estado exata dada por
|Φi = C0|Ψ0i + X ra Car|Ψ r ai + X a<b,r<s Cabrs|Ψ rs abi X a<b<c,r<s<t Cabcrst|Ψ rst abci + . . . , (2.24) onde |Ψr
ai, |Ψrsabi, |Ψrstabci, ..., são determinantes de Slater baseados no determinante usado
na aproximação de Hartree-Fock |Ψ0i, porém substituindo um ou mais orbitais ocupados
(índices a, b, c, ...) por orbitais virtuais desocupados (índices r, s, t,...). Dessa forma, Ψr a
representa um determinante para um estado excitado, onde o elétron que ocupava o orbital Ψa transitou para o orbital virtual Ψr, e assim sucessivamente. Se fosse possível obter a
energia total exata não relativística do sistema através do método variacional e da função de onda exata, seria possível obter a energia de correlação exata através da diferença entre a energia total exata não relativística e a energia total obtida pelo método de Hartree-Fock. Porém, apenas parte dessa energia de correlação é determinada, já que o cálculo exato não é possível, de forma que se utiliza o método da interação de configurações. Neste método, a solução para o problema de N-elétrons é obtida truncando a série 2.24 e usando-a como função tentativa no método variacional. O problema de tal procedimento é o alto custo
computacional, o que deixa o método viável apenas para problemas relativamente simples como moléculas com poucos átomos.
2.1.3
Teoria do Funcional da Densidade
A teoria do funcional da densidade (DFT- Density Functional Theory) baseia-se no fato de que as propriedades do estado fundamental de um sistema podem ser determinadas conhecendo-se a densidade eletrônica no estado fundamental. Isto leva a uma vantagem óbvia, em um sistema com N elétrons, por exemplo, sua função de onda contém 3N variáveis (4N variáveis se levarmos em conta o spin), enquanto sua densidade eletrônica depende de 3 variáveis, ou 4 considerando spin. Em suma, a equação de Schrödinger para N elétrons pode ser escrita como uma equação da densidade eletrônica que possui um número bem menor de variáveis. A DFT foi fundamentada por Hohenberg e Konh em 1964 [93]. No ano seguinte, Kohn e Sham mostraram uma importante aplicabilidade da teoria [94]. Por essas contribuições Walter Kohn veio a ser agraciado com o prêmio Nobel de Química em 1998 [95]. Ultimamente, a teoria do funcional da densidade tem se mostrado como o método mais eficiente para o cálculo de propriedades estruturais e eletrônicas do estado fundamental [90]. Teoremas de Hohenberg-Kohn
Dois teoremas propostos por Hohenberg e Kohn fornecem as bases teóricas para a DFT [93]. Vamos considerar um sistema com N elétrons, com ~ri = (xi, yi, zi) sendo a coor-
denada do i-ésimo elétron. Primeiro Teorema:
O potencial externo v(~r) sentido pelos elétrons é um funcional único da densi- dade eletrônica ρ(~r).
Se ψ é um estado fundamental do sistema, descrito por um hamiltoniano ˆH sub- metido a um potencial externo v(~r), onde ˆH = ˆT + ˆU + ˆV que expressa a energia cinética, somada a energia de interação eletrônica e energia potencial. Suponhamos também a exis- tência de um outro potencial externo v′(~r), resultando em ˆH′ = ˆT + ˆU + ˆV′ e um outro
estado fundamental ψ′. Vamos considerar que os dois potenciais geram a mesma densidade
eletrônica ρ(~r). Pelo princípio variacional temos
E = hψ| ˆT + ˆU + ˆV |ψi < hψ′| ˆT + ˆU + ˆV |ψ′i e (2.25) E′ = hψ′ | ˆT + ˆU + ˆV′ |ψ′ i < hψ| ˆT + ˆU + ˆV′ |ψi. (2.26) Então, podemos escrever
hψ| ˆH|ψi < hψ′ | ˆH|ψ′ i = hψ′ | ˆH′ |ψ′ i + hψ′ | ˆV − ˆV′ |ψ′ i ⇒ E < E′+ hψ′| ˆV − ˆV′|ψ′i. (2.27) Para hψ′| ˆH′|ψ′i, temos E′ < E + hψ| ˆV′− ˆV |ψi, (2.28) como ρ(~r) = hψ| PN
i=1δ(~r − ~ri)|ψi e ˆV =PNi=1v(~ri),
hψ| ˆV |ψi = N X i=1 Z d3r Z d3r1. . . Z d3riv(~r)δ(~r − ~ri) Z d3ri + 1 . . . Z d3rNψ∗ψ = Z ρ(~r)v(~r)d3r. (2.29) Com isso, as expressões 2.27 e 2.28 tornam-se, respectivamente,
E < E′+ Z ρ(~r)[v(~r) − v′(~r)]d3r e E′ < E + Z ρ(~r)[v′(~r) − v(~r)]d3r. Da soma das duas relações acima temos
E + E′ < E′+ E (2.30) que gera um absurdo. Na demostração, assumimos a mesma densidade eletrônica para po- tenciais externos distintos e dessa forma obtemos tal absurdo. Isto ocorre porque as funções de onda escolhidas foram diferentes. Podemos contornar esse problema e concluir que a uni- cidade de ρ(~r) exige que ψ = ψ′. Consequentemente, um observável físico é um funcional
Segundo Teorema:
A energia do estado fundamental é mínima para a densidade exata ρ(~r) e pode ser escrita como:
E[ρ] = hψ| ˆT + ˆU + ˆV |ψi (2.31) onde ψ é a função de onda no estado fundamental
Considere um estado, ψ não necessariamente fundamental que gera uma densidade eletrônica ρ. E seja ρ0 a densidade eletrônica do estado fundamental, proveniente de uma
função de onda do estado fundamental ψ0 temos que
se ρ 6= ρ0 ⇒ ψ 6= ψ0, ou seja, E > E0 e (2.32)
se ρ = ρ0 ⇒ ψ = ψ0, ou seja, E = E0. (2.33)
Considerando a equação 2.31 podemos escrever:
E[ρ] = hψ| ˆT + ˆU |ψi + hψ| ˆV |ψi (2.34) ou
E[ρ] = F [ρ] + hψ| ˆV |ψi, (2.35) onde F [ρ] é um funcional universal válido para qualquer sistema coulombiano e hψ| ˆV |ψi depende do sistema que estamos tratando. De forma análoga, para o estado fundamental temos
E[ρ0] = F [ρ0] + hψ0| ˆV |ψ0i. (2.36)
Aplicando o teorema variacional
E[ψ0] ≤ E[ψ] (2.37)
hψ0| ˆT + ˆU + ˆV |ψ0i ≤ hψ| ˆT + ˆU + ˆV |ψi (2.38)
F [ρ0] + hψ0| ˆV |ψ0i ≤ F [ρ] + hψ| ˆV |ψi (2.39)
E[ρ0] ≤ E[ρ]. (2.40)
Esses dois teoremas garantem a possibilidade de obter o estado fundamental de um sistema sem que se conheça diretamente sua função de onda.
Equações de Kohn-Sham
Dos teoremas de Hohenberg-Kohn sabemos que é possível calcular as propriedades de um sistema através de sua densidade eletrônica, mas não sabemos como calcular tais propri- edades. Para encontrar as equações que permitem realizar cálculos utilizando a densidade eletrônica considere a energia como o funcional
E[ρ] = Z
v(r)ρ(r)d3r + F [ρ], (2.41) onde F [ρ] é o funcional universal discutido anteriormente. Como as interações coulombianas são de longo alcance podemos separar a parte coulombiana clássica do funcional universal,
F [ρ] = 1 2 Z Z ρ(~r)ρ(~r′) |~r − ~r′| d 3rd3r′+ T 0[ρ] + Exc[ρ], (2.42)
onde T0[ρ] é a energia cinética e Exc[ρ] a energia de troca e correlação. Dessa forma a energia
total é dada por E[ρ] = Z v(~r)ρ(~r)d3r +1 2 Z Z ρ(~r)ρ(~r′) |~r − ~r′| d 3rd3r′ + T0[ρ] + Exc[ρ]. (2.43)
Considerando que a densidade eletrônica é dada por ρ(~r) = PN
i=1|ψi(~r)|2, onde N é o número
de elétrons; a energia cinética é T0[ρ] = −12 PiR ψi∗∇2ψid3r; e a ortonormalidade da função
de onda R ψi(~r)ψj(~r)d3~r = δij. Podemos definir o funcional da função de onda
Ω[ρ] = E[ρ] −X
ij
ǫij
Z
ψi(~r)ψjrd3~r. (2.44)
Aqui ǫij são os multiplicadores de Lagrange que asseguram a ortonormalidade das funções
de onda. Minimizando funcional Ω[ρ] com relação a ψ∗
i, chegamos a equação de Kohn-Sham
−12∇2+ vKS(~r) ψi = ǫiψi, i = 1, 2, ..., N, (2.45) onde vKS(~r) = v(~r) + Z ρ(~r) |~r − ~r′|+ δExc δρ . (2.46) Como ρ(~r) = N X i=1 |ψi(~r)|2, (2.47)
temos uma relação de dependência mútua entre vKS(~r) e ρ(~r). A solução para esse problema
vem através de cálculos autoconsistentes, onde o esquema mais adotado é mostrado na Figura 2.2.
O termo δExc
δρ é conhecido como potencial de troca e correlação cuja forma é des-
conhecida. Isto é de importância central na precisão dos cálculos da teoria do funcional da densidade e tem sido realizado inúmeros trabalhos para encontrar formas mais precisas de
δExc
δρ [96]. Dessa forma, um potencial de troca e correlação aproximado deve ser escolhido
de acordo com a precisão, ou a necessidade do problema envolvido. Existem várias apro- ximações para esse termo, as mais utilizadas são a aproximação da densidade local (LDA- Local Density-Functional Approximation) e a aproximação de gradiente generalizado (GGA- Generalized gradient Approximation).
Para aproximação da densidade local, assume-se que a densidade pode ser tratada localmente como um gás de elétrons uniforme, isto é, a energia de troca e correlação em cada ponto do sistema é a mesma que a de um gás de elétrons uniforme com a mesma densidade. Esta aproximação foi originalmente introduzida por Kohn e Sham e é válida para sistemas em que a densidade varia lentamente. Usando esta aproximação a energia de troca correlação para uma densidade é dada por
ELDA xc =
Z
ρ(~r)ǫxc[ρ(~r)]d~r, (2.48)
onde ǫxc é a correlação de troca de energia por partícula de um gás de elétrons homogêneo
de densidade ρ.
No entanto, para sistemas reais sabemos que a densidade de carga não é homogênea e a aproximação LDA consequentemente não é boa para sistemas em que ρ(~r) sofre uma forte variação. Com isso, uma maneira de melhorar a aproximação LDA é considerar o gradiente da densidade de carga ∇ρ(~r). Isso é feito com o seguinte funcional
EGGA xc [ρ] =
Z
ρ(~r)ǫxc[ρ(~r)]Fxc(ρ, ∇ρ)d~r, (2.49)
onde Fxc(ρ, ∇ρ) dá a correção sobre o exchange local. Essa aproximação é conhecida como
GGA. Existem várias propostas para o funcional EGGA
xc atualmente [97] as mais utilizadas são
baseadas nos trabalhos de Perdew-Burke-Ernerhof [98], de Lee-Yang-Parr [99], de Perdew e Wang [100], de Perdew [101] e Becke [102].
Figura 2.2: Ciclo de autoconsistência.
2.1.4
Pseudopotenciais
Além das aproximações e teoremas descritos anteriormente, a complexidade dos cál- culos de estrutura eletrônica de moléculas e sólidos pode ser reduzida significativamente através da teoria dos pseudopotenciais. Neste tratamento não é necessário levar em conta todos os elétrons do sistema, pois é razoável considerar que os estados eletrônicos desses dividam-se em estados do caroço e estados de valência. Os estados do caroço são aqueles que estão fortemente ligados aos núcleos e permanecem praticamente inalterados quando são postos em diferentes ambientes químicos. Já os estados de valência são aqueles responsáveis pelas ligações químicas. Dessa forma, podemos considerar somente os graus de liberdade dos elétrons de valência como uma aproximação válida para os cálculos de estrutura eletrônica de moléculas e sólidos.
diferentes ambientes químicos. Dessa forma, o pseudopotencial deve representar bem o com- portamento de elétrons fora da região de caroço. Para isso é necessário obedecer as seguintes condições:
• Os autovalores reais e pseudo devem ser iguais, e as pseudofunções de onda devem ser iguais as funções de onda reais para distância maiores do que um determinado raio de corte rcorte. Em suma,
ǫpsl = ǫreal
l e (2.50)
Ψpsl = Ψreal
l . (2.51)
Isso para r > rcorte. Para que seja contínua as derivadas de Ψ ps l e Ψ
real
l devem ser iguais
no raios de corte.
• A carga contida na esfera de raio de corte rcorte é igual para as duas funções de onda
(conservação de norma). Isso garante que o potencial eletrostático seja o mesmo fora do raio de corte.
• A derivada logarítmica das funções reais e suas derivadas em relação a energia devem concordar para r maior que o raio de corte. Esse fato garante que as propriedades de espalhamento sejam respondidas também pelo pseudopotencial.
Considerando-se uma blindagem eletrônica com simetria esférica devemos resolver a equação radial de Kohn-Sham: −12 d 2 dr2 + l(l + 1)) 2r2 + V [n, r] rRl = ǫlrRl. (2.52)
Invertendo a equação acima temos, Vlps = ǫl− l(l + 1) 2r2 + 1 2rΨl d2 dr2[rΨ P S l ]. (2.53)
Devemos retirar a blindagem dos elétrons de valência, já que essa depende do ambiente químico em que o pseudopotencial está. Subtraindo então o potencial de Hartree e o potencial de troca e correlação obtemos o pseudopotencial iônico:
Vionps = Vlps− VHps− Vps
xc. (2.54)
independente de l, e uma parte semi-local, de curto alcance e dependente de l de maneira que na forma de operador, o pseudopotencial iônico pode ser escrito como
Vionps = V ps ion,local+ Vsl X l Vsem|lihl|. (2.55)
Para raios grandes Vps
ion,local comporta-se como o potencial eletrostático com a carga equiva-
lente a carga de valência. A parte semi-local é truncada em um valor escolhido de l de forma a reproduzir o espalhamento atômico. A forma sugerida por Kleinman-Bylander para esse potencial semi-local é Vkb = |V slΨpsl ihVslΨpsl | hΨpsl |Vsl|Ψ ps l i , (2.56) onde Ψps
l é a pseudofunção de onda que contém o momento angular no qual foi calculado o
pseudopotencial.
2.2
SIESTA
Neste estudo, realizamos cálculos de estrutura eletrônica como implementado no có- digo SIESTA [103] (Spanish Initiative for Electronic Simulations with Thousands of Atoms) que é um dos programas computacionais mais conhecidos que usam a Teoria do Funcional da Densidade. O SIESTA trabalha com bases localizadas formadas por uma combinação linear de orbitais atômicos de alcance finito. Essas bases permitem calcular elementos de matriz do hamiltoniano de Kohn-Sham com o custo computacional que escala linearmente com o tamanho do sistema. Isso reduz o custo computacional trazendo uma efetiva redução no tempo de cálculo e fornece resultados com uma boa precisão.
2.2.1
Orbitais Atômicos
Para encontrar os orbitais atômicos numéricos utilizados no SIESTA resolve-se a equação de Schrödinger radial para pseudo-átomos isolados em uma simetria radial com a mesma aproximação para sólidos e moléculas. As funções de base a priori localizadas devem
ser determinadas utilizando condições de contorno, ou multiplicando-se os orbitais do átomo livre por uma dada função de corte. Obtidos os orbitais localizados, que necessariamente são nulos em uma região externa, a partir do raio de corte. Três condições devem ser obedecidas nesse tipo de base:
1. o número de orbitais por átomo;
2. o alcance dos raios de corte dos orbitais;
3. a forma de confinamento do orbitais atômicos numéricos.
Podemos utilizar orbitais atômicos numéricos com base simples SZ (single-ζ), ou com bases mais completas (double, multiple−ζ). Adicionando uma função chamada de polarização (P) podemos dar conta da flexibilidade angular. A base SZ também chamada de base mínima, possui somente uma função radial por momento angular, apenas para os estados ocupados na valência do átomo isolado. Com isso esta base requer um pequeno número de funções para descrever os elétrons de um átomo. Tal base é útil na na descrição qualitativa das ligações químicas e traz uma boa descrição para a banda de valência. Tendo necessidade de uma maior flexibilidade na parte angular adiciona-se uma segunda função por momento angular. Esta segunda função é conhecida como double − ζ. Este segundo orbital numérico deve reproduzir a função de onda a partir de um determinado raio e deve ser suave na origem, com rl(a − br2) onde os parâmetros a e b se ajustam de forma que esta função e sua derivada
sejam contínuas em um determinado raio re. Tal raio é tornado fixo de forma que a assíntota
do orbital tenha o valor de sua norma determinado. Dessa forma, gera-se o mesmo espaço de Hilbert tomando uma segunda função como sendo a diferença entre a função de onda original e essa nova função suave. Essa nova função é estritamente localizada em um raio re e reduz
o custo computacional. Esse mesmo procedimento é utilizado para calcular multiple − ζ, apenas com a escolha de valores para o raio re.
A vantagem de se usar orbitais atômicos estritamente localizados é que eles se anulam acima de um determinado raio de corte. O problema é encontrar uma maneira sistemática de definir todos os raios das funções de base, uma vez que tanto a exatidão como a eficiência computacional dependem deles. O modelo usual no qual todos os raios são definidos em uma só função é a correção de energia (energy shift), isto é, um incremento de energia sofrido pelo orbital quando está confinado. Tal processo aumenta a curvatura do orbital e sua energia cinética. Limitando-se todos os raios de maneira que este incremento seja o mesmo para todos os orbitais, gera-se uma base que evita a transferência de carga.
As funções de base devem se adaptar a forma do pseudopotencial na região próxima ao núcleo. Isso é obtido utilizando como base as soluções do hamiltoniano de Kohn-Sham para o pseudopotencial correspondente ao átomo livre. A forma dos orbitais para os raios maiores que o raio de corte depende da maneira pela qual se produz o confinamento. O potencial confinante em sua forma mais usual é importante porque evita problemas de descontinuidade abruptos, anulando-se a região do caroço, sendo contínuo assim em todas suas derivadas, a partir de um raio de corte interno ri e divergindo em rc, assegurando assim uma localização
suave.
Em nossos cálculos de estrutura eletrônica, utilizamos a teoria do funcional da densidade [94], com tratamento de troca e correlação dado pela aproximação de gradiente generalizado GGA parametrizado por Perdew, Burke e Ernzerhof [98] como implementado no código SIESTA [104]. Uma combinação linear de pseudo-orbitais atômicos com base double zeta e polarização de funções é usada [105, 106]. A interação entre nucleo iônico e elétrons de valência é dada pelo pseudopotencial de norma conservada de Troullier-Martins [107], na forma fatorizada de Kleinman-Bylander [108]. O critério de convergência dos cálculos autoconsistentes é de 10−4 eV para a energia total e densidade eletrônica. A otimização
completa de todas as posições atômicas é realizada com uma força atômica inferior a 0,1 eV/Å.