Simulated Structures: Computational Modeling in Modern Materials Science

A física do estado sólido moderna apoia-se sobre um paradoxo fundamental. As leis da mecânica quântica que governam os elétrons e os núcleos atômicos são conhecidas com extrema precisão desde as primeiras décadas do século XX. No entanto, aplicar essas leis para prever as propriedades de um cristal real contendo 10²³ partículas é uma tarefa impraticável do ponto de vista matemático. A equação de Schrödinger dependente do tempo para um sistema de múltiplos corpos contendo $N$ núcleos e $Z N$ elétrons exige a resolução de uma função de onda complexa de $3(N + Z N)$ variáveis espaciais, além das coordenadas de spin. Mesmo para um pequeno aglomerado de poucas dezenas de átomos de ferro, a quantidade de memória necessária para armazenar a função de onda excede a capacidade de qualquer supercomputador existente ou projetável. O desafio central da física da matéria condensada e da química quântica não é encontrar novas leis fundamentais, mas desenvolver aproximações criteriosas que tornem essas equações solváveis sem perder a fidelidade física necessária para prever propriedades de materiais reais.

A primeira simplificação estrutural para esse problema é a aproximação de Born-Oppenheimer. Fundamentada na diferença substancial de massa entre os núcleos atômicos e os elétrons — a massa do próton é aproximadamente 1836 vezes maior que a do elétron —, essa abordagem considera os núcleos atômicos como estáticos do ponto de vista dos elétrons. Como a dinâmica eletrônica ocorre em uma escala de tempo da ordem de attossegundos, enquanto o movimento nuclear acontece em femtossegundos, os elétrons respondem quase instantaneamente a qualquer mudança na posição dos núcleos. Isso permite separar a função de onda total do sistema em uma parte nuclear e uma parte eletrônica. Matematicamente, a equação de Schrödinger eletrônica assume a seguinte forma:

$$\hat{H}_{eletr} \Psi(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N) = E_{eletr} \Psi(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N)$$

Onde o hamiltoniano eletrônico $\hat{H}_{eletr}$ é composto por três termos distintos:

$$\hat{H}_{eletr} = -\frac{\hbar^2}{2m_e} \sum_{i=1}^{N} \nabla_i^2 - \sum_{i=1}^{N} \sum_{A=1}^{M} \frac{Z_A e^2}{4\pi\epsilon_0 |\mathbf{r}_i - \mathbf{R}_A|} + \frac{1}{2} \sum_{i=1}^{N} \sum_{j \neq i}^{N} \frac{e^2}{4\pi\epsilon_0 |\mathbf{r}_i - \mathbf{r}_j|}$$

O primeiro termo representa a energia cinética de todos os elétrons do sistema. O segundo termo descreve a atração eletrostática atrativa entre os elétrons nas posições $\mathbf{r}_i$ e os núcleos fixos nas posições $\mathbf{R}_A$ com carga atômica $Z_A$. Este termo é frequentemente denotado como o potencial externo $V_{ext}(\mathbf{r})$. O terceiro termo corresponde à repulsão Coulombiana entre os próprios elétrons. É justamente este último termo — a interação elétron-elétron — que impede a separação direta das equações em equações de um único elétron. O movimento de cada elétron está acoplado de forma instantânea ao movimento de todos os outros, um fenômeno conhecido como correlação eletrônica.

Apesar de simplificar a física do sistema, a aproximação de Born-Oppenheimer mantém o problema da dimensão excessiva. A função de onda eletrônica de muitos corpos $\Psi(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N)$ continua sendo um objeto matemático intrincado e impossível de ser computado diretamente para materiais complexos. Durante décadas, estratégias como o método de Hartree-Fock tentaram resolver o problema aproximando $\Psi$ como um único determinante de Slater de orbitais monoeletrônicos. Esse ajuste inclui a troca quântica exata exigida pelo princípio de exclusão de Pauli, mas ignora completamente a correlação eletrônica instantânea, o que gera erros substanciais em estimativas de energias de ligação, comprimentos de rede e propriedades eletrônicas.

Os Teoremas de Hohenberg-Kohn e a Densidade Eletrônica

Uma mudança de paradigma ocorreu em 1964, quando Pierre Hohenberg e Walter Kohn publicaram dois teoremas matemáticos que redefiniram a física da matéria condensada. A ideia central da Teoria do Funcional da Densidade (DFT) é substituir a função de onda de muitos corpos $\Psi$, que depende de $3N$ variáveis espaciais, pela densidade eletrônica do estado fundamental $n(\mathbf{r})$, um escalar que depende unicamente de três coordenadas espaciais $(x, y, z)$. A densidade eletrônica representa a probabilidade de encontrar qualquer um dos $N$ elétrons em um determinado volume no espaço e é definida como:

$$n(\mathbf{r}) = N \int |\Psi(\mathbf{r}, \mathbf{r}_2, \dots, \mathbf{r}_N)|^2 d\mathbf{r}_2 \dots d\mathbf{r}_N$$

A redução de $3N$ variáveis para apenas três coordenadas espaciais gera um ganho computacional extraordinário. No entanto, para que essa substituição faça sentido do ponto de vista físico, é necessário provar rigorosamente que a densidade contém todas as informações do sistema de muitos corpos. Isso foi demonstrado através dos dois teoremas de Hohenberg-Kohn.

Primeiro Teorema de Hohenberg-Kohn: O potencial externo $V_{ext}(\mathbf{r})$ atuando sobre um sistema de elétrons é determinado de forma única, a menos de uma constante aditiva, pela densidade eletrônica do estado fundamental $n(\mathbf{r})$.

A prova deste teorema é realizada por contradição (reductio ad absurdum). Suponha que existam dois potenciais externos distintos, $V_{ext}^{(1)}(\mathbf{r})$ e $V_{ext}^{(2)}(\mathbf{r})$, que difiram por mais do que uma constante, mas que resultem exatamente na mesma densidade eletrônica $n(\mathbf{r})$ no estado fundamental. Esses dois potenciais dariam origem a dois hamiltonianos distintos, $\hat{H}^{(1)}$ e $\hat{H}^{(2)}$, com duas funções de onda do estado fundamental diferentes, $\Psi^{(1)}$ e $\Psi^{(2)}$. Usando o princípio variacional para a energia do estado fundamental:

$$E^{(1)} < \langle \Psi^{(2)} | \hat{H}^{(1)} | \Psi^{(2)} \rangle = \langle \Psi^{(2)} | \hat{H}^{(2)} | \Psi^{(2)} \rangle + \langle \Psi^{(2)} | \hat{H}^{(1)} - \hat{H}^{(2)} | \Psi^{(2)} \rangle$$

Isso leva à relação:

$$E^{(1)} < E^{(2)} + \int n(\mathbf{r}) \left[ V_{ext}^{(1)}(\mathbf{r}) - V_{ext}^{(2)}(\mathbf{r}) \right] d\mathbf{r}$$

Invertendo os rótulos 1 e 2, obtém-se uma expressão simétrica:

$$E^{(2)} < E^{(1)} + \int n(\mathbf{r}) \left[ V_{ext}^{(2)}(\mathbf{r}) - V_{ext}^{(1)}(\mathbf{r}) \right] d\mathbf{r}$$

Ao somar as duas desigualdades, chega-se ao resultado absurdo $E^{(1)} + E^{(2)} < E^{(1)} + E^{(2)}$. Portanto, não podem existir dois potenciais externos distintos gerando a mesma densidade eletrônica do estado fundamental. Como $n(\mathbf{r})$ determina $V_{ext}(\mathbf{r})$ e o número total de elétrons $N = \int n(\mathbf{r}) d\mathbf{r}$, a densidade fixa o hamiltoniano completo e, por consequência, todas as propriedades quânticas do estado fundamental do sistema.

Segundo Teorema de Hohenberg-Kohn: Existe um funcional universal para a energia $E[n]$ em termos da densidade eletrônica $n(\mathbf{r})$, válido para qualquer potencial externo $V_{ext}(\mathbf{r})$. A energia exata do estado fundamental do sistema é o valor mínimo absoluto deste funcional, e a densidade que minimiza o funcional é a densidade real do estado fundamental.

A energia total pode ser escrita como:

$$E[n] = \int n(\mathbf{r}) V_{ext}(\mathbf{r}) d\mathbf{r} + F_{HK}[n]$$

Onde $F_{HK}[n]$ é o funcional universal de Hohenberg-Kohn, que contém a energia cinética dos elétrons e a energia de interação elétron-elétron:

$$F_{HK}[n] = T[n] + V_{ee}[n]$$

A beleza teórica dos teoremas de Hohenberg-Kohn é incontestável, mas eles apresentam uma limitação prática importante: não revelam a forma matemática exata de $F_{HK}[n]$. Sem conhecer essa expressão, o formalismo permanece inviável para cálculos numéricos do dia a dia.

O Esquema Mapeado de Kohn-Sham

A solução prática para o impasse do funcional universal foi apresentada por Walter Kohn e Lu Jeu Sham em 1965. Eles propuseram mapear o sistema real de elétrons interagentes — que é extremamente difícil de resolver — em um sistema hipotético de elétrons não interagentes que se movimentam sob a influência de um potencial efetivo local, $V_{eff}(\mathbf{r})$, escolhido de tal forma que a densidade eletrônica resultante seja idêntica à do sistema real interagente.

Para esse sistema fictício sem interação direta, a energia cinética pode ser calculada exatamente a partir de orbitais monoeletrônicos $\psi_i(\mathbf{r})$, denominados orbitais de Kohn-Sham:

$$T_s[n] = -\frac{\hbar^2}{2m_e} \sum_{i=1}^{N} \int \psi_i^*(\mathbf{r}) \nabla^2 \psi_i(\mathbf{r}) d\mathbf{r}$$

A densidade eletrônica total do sistema é dada pela soma do módulo ao quadrado desses orbitais ocupados:

$$n(\mathbf{r}) = \sum_{i=1}^{N} |\psi_i(\mathbf{r})|^2$$

Kohn e Sham reescreveram o funcional de energia total dividindo-o em parcelas conhecidas e isolando a parte desconhecida em um único termo:

$$E[n] = T_s[n] + \int n(\mathbf{r}) V_{ext}(\mathbf{r}) d\mathbf{r} + E_H[n] + E_{xc}[n]$$

Onde $E_H[n]$ é a energia de Hartree, que descreve a repulsão eletrostática clássica da distribuição de carga contínua $n(\mathbf{r})$ consigo mesma:

$$E_H[n] = \frac{e^2}{8\pi\epsilon_0} \iint \frac{n(\mathbf{r}) n(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r} d\mathbf{r}'$$

O termo $E_{xc}[n]$ é o funcional de troca e correlação. Ele atua como uma lixeira conceitual para a física desconhecida do sistema, contendo:

• A diferença entre a energia cinética real de muitos corpos e a energia cinética fictícia de elétrons não interagentes ($T[n] - T_s[n]$).
• As correções quânticas de troca resultantes do princípio de exclusão de Pauli.
• As correções de correlação eletrônica derivadas de forças instantâneas repulsivas intereletrônicas.
• A correção do termo de autointeração espúrio presente na energia de Hartree.

Ao aplicar o princípio variacional minimizando a energia em relação aos orbitais de Kohn-Sham sob a restrição de ortonormalidade ($\langle \psi_i | \psi_j \rangle = \delta_{ij}$), obtém-se o conjunto de equações mono-eletrônicas de Kohn-Sham:

$$\left( -\frac{\hbar^2}{2m_e} \nabla^2 + V_{eff}(\mathbf{r}) \right) \psi_i(\mathbf{r}) = \epsilon_i \psi_i(\mathbf{r})$$

O potencial efetivo $V_{eff}(\mathbf{r})$ é dado por:

$$V_{eff}(\mathbf{r}) = V_{ext}(\mathbf{r}) + V_H(\mathbf{r}) + V_{xc}(\mathbf{r})$$

Onde o potencial de Hartree é $V_H(\mathbf{r}) = \frac{e^2}{4\pi\epsilon_0} \int \frac{n(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r}'$ e o potencial de troca e correlação é definido como a derivada funcional:

$$V_{xc}(\mathbf{r}) = \frac{\delta E_{xc}[n]}{\delta n(\mathbf{r})}$$

Estas equações devem ser resolvidas de maneira autoconsistente. Como o potencial efetivo depende da densidade eletrônica $n(\mathbf{r})$, e a própria densidade depende dos orbitais $\psi_i(\mathbf{r})$ obtidos pela resolução da equação de Kohn-Sham, adota-se um ciclo iterativo de cálculo até alcançar a convergência numérica.

Aproximações do Funcional de Troca e Correlação

Embora a formulação de Kohn-Sham seja exata, a forma analítica de $E_{xc}[n]$ permanece desconhecida. Toda a utilidade prática da DFT depende do uso de aproximações eficientes para esse funcional. A evolução dessas aproximações é representada metaforicamente como a "Escada de Jacó" dos funcionais, proposta por John Perdew, em que cada degrau adiciona novos parâmetros físicos para melhorar a precisão.

1. Aproximação da Densidade Local (LDA)

O primeiro degrau da escala assume que o funcional de troca e correlação em cada ponto do espaço é idêntico ao de um gás de elétrons homogêneo de mesma densidade. A energia de troca e correlação é expressa por:

$$E_{xc}^{LDA}[n] = \int n(\mathbf{r}) \epsilon_{xc}^{unif}(n(\mathbf{r})) d\mathbf{r}$$

Onde $\epsilon_{xc}^{unif}(n)$ é a energia de troca e correlação por partícula de um gás de elétrons uniforme. A parte de troca é conhecida analiticamente por equações teóricas de Paul Dirac, enquanto a parte de correlação foi calculada com altíssima precisão via simulações de Monte Carlo Quântico por Ceperley e Alder em 1980. A LDA funciona surpreendentemente bem em sistemas onde a densidade varia suavemente, como em metais simples (alumínio, sódio). No entanto, tende a superestimar as energias de ligação (gerando supercoesão) e a subestimar comprimentos de ligação e constantes de rede em moléculas e isolantes.

2. Aproximação de Gradiente Generalizado (GGA)

Para tratar sistemas com variações espaciais mais bruscas na densidade, o segundo degrau introduz a dependência em relação ao gradiente local da densidade eletrônica, $\nabla n(\mathbf{r})$:

$$E_{xc}^{GGA}[n] = \int n(\mathbf{r}) \epsilon_{xc}^{GGA}(n(\mathbf{r}), |\nabla n(\mathbf{r})|) d\mathbf{r}$$

Os funcionais GGA mais difundidos no estudo de materiais são o PBE (Perdew-Burke-Ernzerhof) e o BLYP (Becke-Lee-Yang-Parr). O PBE corrige os excessos da LDA, relaxando as ligações, o que resulta em parâmetros de rede maiores e melhor concordância com energias de coesão em cristais e moléculas orgânicas. No entanto, o PBE pode subestimar ligeiramente a coesão em redes metálicas densas, motivo pelo qual foi desenvolvido o funcional PBEsol, ajustado para sólidos cristalinos.

3. Funcionais Híbridos e Meta-GGA

Degraus superiores adicionam a laplaciana da densidade $\nabla^2 n(\mathbf{r})$ ou a densidade de energia cinética dos orbitais (meta-GGA, como o funcional SCAN). Subindo ainda mais na hierarquia, surgem os funcionais híbridos, que misturam uma fração da energia de troca exata calculada pelo método de Hartree-Fock com a troca e correlação da GGA. O exemplo mais notável na ciência das superfícies e semicondutores é o HSC06 (Heyd-Scuseria-Ernzerhof) e o B3LYP na química teórica. Os funcionais híbridos corrigem o erro de autointeração de funcionais semi-locais e fornecem estimativas de bandgaps em semicondutores muito superiores às da LDA e GGA.

Conjuntos de Bases, Teorema de Bloch e Pseudopotenciais

Para transformar as equações diferenciais de Kohn-Sham em equações matriciais passíveis de solução computacional, é necessário expandir os orbitais eletrônicos $\psi_i(\mathbf{r})$ em um conjunto de funções de base conhecidas. A escolha desse conjunto varia conforme a geometria do sistema estudado.

Em sistemas periódicos tridimensionais, como cristais de semicondutores, óxidos ou ligas metálicas, a simetria de translação da rede facilita a aplicação do Teorema de Bloch. Este teorema estabelece que os orbitais eletrônicos em um potencial periódico $V(\mathbf{r} + \mathbf{R}) = V(\mathbf{r})$ podem ser escritos como o produto de uma onda plana por uma função $u_{n,\mathbf{k}}(\mathbf{r})$ com a mesma periodicidade da rede cristalina:

$$\psi_{n,\mathbf{k}}(\mathbf{r}) = e^{i\mathbf{k}\cdot\mathbf{r}} u_{n,\mathbf{k}}(\mathbf{r})$$

Onde $n$ é o índice da banda de energia, $\mathbf{k}$ é o vetor de onda situado dentro da primeira Zona de Brillouin no espaço recíproco, e $\mathbf{R}$ representa qualquer vetor de translação da rede direta. Como $u_{n,\mathbf{k}}(\mathbf{r})$ é periódica na rede cristalina, ela pode ser expandida em uma série de Fourier em termos dos vetores da rede recíproca $\mathbf{G}$:

$$u_{n,\mathbf{k}}(\mathbf{r}) = \sum_{\mathbf{G}} c_{n,\mathbf{k}+\mathbf{G}} e^{i\mathbf{G}\cdot\mathbf{r}}$$

Substituindo isso na equação de Bloch, o orbital é representado diretamente por uma soma de ondas planas:

$$\psi_{n,\mathbf{k}}(\mathbf{r}) = \sum_{\mathbf{G}} c_{n,\mathbf{k}+\mathbf{G}} e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}}$$

A vantagem das ondas planas é sua ortogonalidade e simplicidade computacional, o que facilita o cálculo de derivadas operacionais (como a energia cinética) via Transformadas Rápidas de Fourier (FFT). Contudo, tentar descrever as oscilações ultrarrápidas das funções de onda dos elétrons mais internos (elétrons de caroço) próximos ao núcleo atômico exigiria um número absurdamente grande de ondas planas na série de Fourier.

Para solucionar esse entrave, adota-se a aproximação de Pseudopotenciais. A maioria das propriedades físicas e químicas dos materiais — como condutividade, módulo de compressão, absorção óptica e reatividade — depende quase exclusivamente dos elétrons de valência. Os elétrons de caroço estão fortemente ligados e permanecem praticamente inalterados durante as reações químicas ou mudanças de fase.

O conceito de pseudopotencial consiste em substituir o potencial Coulombiano divergente do núcleo $-Z e^2/r$ e os elétrons de caroço por um potencial efetivo suave que atua apenas sobre os elétrons de valência. Além disso, a função de onda real com nós ortogonais perto do núcleo é substituída por uma "pseudofunção de onda" suave dentro de um raio de corte $r_c$, coincidindo exatamente com a função de onda real para distâncias $r > r_c$.

Existem diferentes abordagens de pseudopotenciais:

• Pseudopotenciais Norm-Conserving: Exigem que a carga total contida dentro da esfera de raio $r_c$ seja idêntica para a função de onda real e para a pseudofunção de onda. Garante reprodutibilidade em diferentes ambientes químicos, mas exige um número elevado de ondas planas para elementos com orbitais $2p$ ou $3d$ muito localizados (como oxigênio e metais de transição).
• Pseudopotenciais Ultrasoft (Vanderbilt): Relaxam a condição de conservação de norma, permitindo o uso de pseudofunções de onda substancialmente mais suaves. Isso reduz dramaticamente o número de ondas planas necessárias no cálculo, diminuindo o custo computacional.
• Método Projector Augmented-Wave (PAW): Desenvolvido por Peter Blöchl, o método PAW combina a flexibilidade formal do uso de pseudopotenciais ultrasoft com o rigor de ter acesso completo às funções de onda de todos os elétrons perto do núcleo, reconstruindo analyticamente as oscilações através de operadores de projeção. O método PAW tornou-se o padrão na maioria dos pacotes de simulação de materiais, como VASP, Quantum ESPRESSO e ABINIT.

Parâmetros Práticos e Critérios de Convergência

Realizar uma simulação de DFT confiável exige a configuração criteriosa de parâmetros numéricos para garantir que os resultados físicos sejam independentes do tamanho da base matemática e da discretização do espaço. Dois parâmetros são fundamentais nesse processo: a energia de corte e a malha de amostragem no espaço recíproco.

Energia de Corte (Cutoff Energy)

A expansão formal da função de onda em ondas planas $e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}}$ exige, teoricamente, um número infinito de vetores da rede recíproca $\mathbf{G}$. Na prática, truncamos a série limitando a energia cinética máxima das ondas planas incluídas no conjunto de base:

$$\frac{\hbar^2 |\mathbf{k}+\mathbf{G}|^2}{2m_e} \leq E_{cut}$$

Essa energia limite é denominada energia de corte (ecut). Quanto maior o valor de $E_{cut}$, maior a base de ondas planas, maior a precisão dos cálculos e maior o tempo computacional necessário. Para determinar o valor adequado de $E_{cut}$, realiza-se um **estudo de convergência**: executa-se o mesmo cálculo variando $E_{cut}$ progressivamente até que a variação na energia total por átomo seja inferior a um limite preestabelecido (por exemplo, 1 meV/átomo).

Amostragem do Espaço Recíproco ($k$-points)

Cálculos de grandezas como a densidade eletrônica total $n(\mathbf{r})$ ou a energia de Hartree exigem a integração de funções sobre toda a primeira Zona de Brillouin (BZ) no espaço recíproco:

$$n(\mathbf{r}) = \frac{1}{\Omega_{BZ}} \sum_{n} \int_{BZ} f(\epsilon_{n,\mathbf{k}}) |\psi_{n,\mathbf{k}}(\mathbf{r})|^2 d\mathbf{k}$$

A integração contínua é aproximada por uma soma discreta sobre um conjunto finito de pontos $\mathbf{k}$. A estratégia padrão mais empregada é o esquema de **Monkhorst-Pack**, que gera uma grade homogênea e uniforme de pontos no espaço recíproco (ex.: $6 \times 6 \times 6$ ou $12 \times 12 \times 12$).

A densidade de pontos $\mathbf{k}$ necessária varia de acordo com o comportamento eletrônico do material:

• Isolantes e Semicondutores: Apresentam bandas de energia contínuas separadas por um gap de energia. A integração é numericamente bem-comportada, exigindo grades de pontos $\mathbf{k}$ relativamente esparsas (como $4 \times 4 \times 4$).
• Metais: Possuem superfícies de Fermi complexas, onde as bandas cruzam o nível de Fermi. Isso gera descontinuidades abruptas na função de ocupação $f(\epsilon_{n,\mathbf{k}})$, demandando grades de pontos $\mathbf{k}$ densas (como $16 \times 16 \times 16$ ou maiores) para obter convergência adequada da energia.

Para contornar o problema da superfície de Fermi em metais e evitar instabilidades numéricas no ciclo de autoconsistência, utilizam-se técnicas de **smearing** (alongamento). Essas técnicas substituem o degrau de ocupação por uma função suave que simula uma temperatura eletrônica fictícia. Os esquemas mais populares são o smearing de Fermi-Dirac, o esquema de Methfessel-Paxton e o método dos tetraedros ajustado por Blöchl.

Estudo de Caso Prático: Propriedades do Silício Cristalino

Para aplicar o formalismo teórico da DFT em um problema prático da ciência dos materiais, vamos detalhar a determinação do parâmetro de rede de equilíbrio, do módulo de compressão $B_0$ e da estrutura de bandas eletrônicas do silício cristalino na estrutura diamante (grupo espacial $Fd\bar{3}m$).

A célula unitária do silício cristalino na estrutura diamante contém 8 átomos de silício (ou 2 átomos na célula primitiva FCC). Cada átomo coordena-se tetraedricamente com quatro vizinhos imediatos por meio de ligações covalentes $sp^3$.

1. Cálculo do Parâmetro de Rede e Módulo de Compressão

O primeiro passo de uma caracterização computacional é o relaxamento estrutural para encontrar a geometria de menor energia. Variamos o volume da célula unitária primitiva do silício de forma controlada ao redor do volume esperado e calculamos a energia total $E(V)$ para cada volume $V$. Nesses cálculos, os parâmetros de convergência adotados foram um cutoff de ondas planas $E_{cut} = 400 \text{ eV}$ e uma grade de Monkhorst-Pack de $8 \times 8 \times 8$ para a amostragem de pontos $\mathbf{k}$.

Os dados computados de energia total em função do volume são ajustados por uma equação de estado fenomenológica. Uma das formulações mais robustas para sólidos sob compressão isotrópica é a **Equação de Estado de Birch-Murnaghan** de 3ª ordem:

$$E(V) = E_0 + \frac{9 V_0 B_0}{16} \left[ \left( \left( \frac{V_0}{V} \right)^{\frac{2}{3}} - 1 \right)^3 B_0' + \left( \left( \frac{V_0}{V} \right)^{\frac{2}{3}} - 1 \right)^2 \left( 6 - 4 \left( \frac{V_0}{V} \right)^{\frac{2}{3}} \right) \right]$$

Nesta expressão, as variáveis ajustáveis são:

• $E_0$: Energia total no estado fundamental minimizado.
• $V_0$: Volume de equilíbrio do cristal.
• $B_0$: Módulo de compressão isotrópico no volume de equilíbrio ($B_0 = -V \left( \frac{\partial P}{\partial V} \right)_T$).
• $B_0'$: Derivada do módulo de compressão em relação à pressão ($\frac{\partial B_0}{\partial P}$).

A tabela a seguir apresenta os dados brutos calculados para a energia do silício em função do volume por átomo, utilizando o funcional PBE-GGA:

Volume / átomo ($\text{\AA}^3$)	Parâmetro de Rede $a$ ($\text{\AA}$)	Energia Total / átomo ($\text{eV}$)
18.15	5.260	-107.8210
18.88	5.330	-107.8895
19.63	5.400	-107.9254
20.02	5.435	-107.9312
20.40	5.470	-107.9288
21.19	5.540	-107.9015
22.01	5.610	-107.8480

Ao realizar o ajuste não linear dos dados da tabela acima pela Equação de Birch-Murnaghan, obtém-se o volume de equilíbrio $V_0 = 20.04 \text{ \AA}^3/\text{átomo}$, que corresponde a um parâmetro de rede cúbico de $a_0 = (8 \times V_0)^{1/3} = 5.433 \text{ \AA}$. O módulo de compressão resultante é $B_0 = 88.2 \text{ GPa}$.

2. Cálculo da Estrutura de Bandas Eletrônicas

Uma vez determinada a estrutura cristalina de menor energia (volume de equilíbrio $V_0$), realiza-se o cálculo da estrutura de bandas eletrônicas. Para isso, resolvem-se as equações de Kohn-Sham ao longo de um caminho contínuo de alta simetria na primeira Zona de Brillouin da rede FCC. Os pontos críticos padrão no espaço recíproco são:

• $\Gamma = (0, 0, 0)$: Centro da Zona de Brillouin.
• $X = \left(\frac{1}{2}, 0, \frac{1}{2}\right)$: Centro das faces quadradas no espaço recíproco.
• $L = \left(\frac{1}{2}, \frac{1}{2}, \frac{1}{2}\right)$: Centro das faces hexagonais.
• $W = \left(\frac{1}{2}, \frac{1}{4}, \frac{3}{4}\right)$: Vértices da BZ.

O gráfico de dispersão $E_n(\mathbf{k})$ revela o comportamento semicondutor do silício. O topo da banda de valência (VBM) está localizado exatamente no centro da zona de Brillouin, no ponto $\Gamma$, e é triplamente degenerado (desconsiderando o acoplamento spin-órbita). O mínimo da banda de condução (CBM) não ocorre no ponto $\Gamma$, mas ao longo do segmento de alta simetria $\Gamma \rightarrow X$, localizado a aproximadamente $85\%$ da distância em direção à borda do ponto $X$. Essa descontinuidade espacial entre o máximo da valência e o mínimo da condução caracteriza o silício como um semicondutor de bandgap indireto.

Análise Comparativa e Validação Experimental

A validação dos métodos computacionais fundamenta-se na comparação direta entre valores previstos por simulação e medidas experimentais obtidas por técnicas analíticas avançadas, como a difração de raios X (XRD) e o espalhamento inelástico de nêutrons.

A tabela abaixo compara os parâmetros de rede ($a_0$), o módulo de compressão ($B_0$) e o bandgap fundamental ($E_g$) calculados por diferentes funcionais de troca e correlação com os dados experimentais obtidos a 300 K:

Método / Funcional	$a_0$ ($\text{\AA}$)	Erro $a_0$ (%)	$B_0$ ($\text{GPa}$)	$E_g$ ($\text{eV}$)	Tipo de Bandgap
LDA	5.381	-0.92%	96.5	0.48	Indireto ($\Gamma \rightarrow X$)
GGA (PBE)	5.433	+0.04%	88.2	0.62	Indireto ($\Gamma \rightarrow X$)
Híbrido (HSE06)	5.435	+0.07%	90.1	1.11	Indireto ($\Gamma \rightarrow X$)
Experimental (300K)	5.431	-	97.88	1.12	Indireto ($\Gamma \rightarrow X$)

A análise comparativa revela tendências sistemáticas bem conhecidas no meio acadêmico:

• A LDA subestima o parâmetro de rede de equilíbrio em aproximadamente $0.9\%$, fenômeno conhecido como supercoesão. Por prever ligações atômicas mais curtas e rígidas, a LDA sobrestima o módulo de compressão ($96.5 \text{ GPa}$).
• O funcional GGA (PBE) produz excelente concordância para a geometria cristalina do silício, errando o parâmetro de rede por apenas $+0.04\%$, um desvio praticamente dentro da incerteza térmica experimental. O módulo de compressão é ligeiramente subestimado ($88.2 \text{ GPa}$).
• O Problema do Bandgap: Tanto a LDA ($0.48 \text{ eV}$) quanto o PBE ($0.62 \text{ eV}$) subestimam severamente o valor do bandgap real do silício ($1.12 \text{ eV}$). Essa deficiência não se deve à falta de convergência numérica, mas é uma limitação formal das equações de Kohn-Sham com funcionais semilocais, que não incluem a descontinuidade do potencial de troca e correlação ao adicionar um elétron extra ao sistema.
• O funcional híbrido HSE06 recupera a física correta dos autoestados não ocupados ao introduzir troca exata de Hartree-Fock a curtas distâncias, prevendo um bandgap de $1.11 \text{ eV}$, em consonância com as medidas experimentais.

Limitações da DFT e o Tratamento de Sistemas Correlacionados

Apesar de seu sucesso estrondoso na descrição de metais simples, semicondutores covalentes e compostos orgânicos, a Teoria do Funcional da Densidade padrão (LDA e GGA) falha criticamente quando aplicada a um grupo amplo de materiais avançados conhecidos como **sistemas fortemente correlacionados**.

Essa classe abrange compostos contendo elementos de transição com orbitais $3d$ ou $4f$ parcialmente preenchidos, tais como óxidos de níquel ($\text{NiO}$), óxidos de ferro ($\text{FeO}$), cupratos supercondutores e terras raras. Nesses sistemas, os elétrons estão fortemente localizados ao redor dos sítios atômicos e as repulsões Coulombianas intereletrônicas diretas dominam a dinâmica do material. A DFT semilocal espalha a densidade eletrônica desses orbitais localizados para delocalizá-los erroneamente pela rede (erro de delocalização resultante da autointeração espúria).

O caso clássico é o do $\text{FeO}$ e do $\text{NiO}$. Experimentos mostram que ambos são isolantes Mott-Hubbard antiferromagnéticos com bandgaps expressivos. No entanto, cálculos realizados via LDA ou PBE prevêem incorretamente que o $\text{NiO}$ e o $\text{FeO}$ são condutores metálicos. Esse descompasso inviabiliza a modelagem de baterias de íon-lítio à base de catodos de fosfato de ferro-lítio ($\text{LiFePO}_4$) ou óxidos de cobalto ($\text{LiCoO}_2$).

A Correção DFT+U (Modelo de Hubbard)

Para contornar essa deficiência sem elevar o custo computacional ao nível da química quântica de muitos corpos, Anisimov e colaboradores desenvolveram o método **DFT+U**. Esta abordagem combina a descrição de DFT padrão para elétrons delocalizados de valência ($s$ e $p$) com um termo corretivo extra, derivado do Modelo de Hubbard, aplicado estritamente aos orbitais fortemente localizados ($d$ ou $f$):

$$E_{DFT+U} = E_{DFT} + E_{Hubbard} - E_{double-counting}$$

Onde o termo corretivo em uma formulação simplificada (proposta por Dudarev) depende de um único parâmetro efetivo $U_{eff} = U - J$:

$$E_{Hubbard} = \frac{U_{eff}}{2} \sum_{\sigma} \text{Tr} \left( \mathbf{n}^{\sigma} ( \mathbf{1} - \mathbf{n}^{\sigma} ) \right)$$

Aqui, $\mathbf{n}^{\sigma}$ é a matriz de ocupação atômica dos orbitais localizados para o spin $\sigma$, $U$ representa a repulsão Coulombiana no mesmo sítio atômico, e $J$ é o parâmetro de troca de Hund. A penalidade energética introduzida pelo termo $+U$ força a ocupação dos orbitais a assumir valores inteiros (completamente cheios ou completamente vazios), penalizando ocupações fracionárias artificiais. Isso abre um gap de energia característico de isolantes de Mott e corrige as posições relativas dos níveis eletrônicos d nos espectros de fotoemissão.

Outra limitação crítica dos funcionais semilocais é a incapacidade de descrever forças de dispersão a longa distância, como as interações de van der Waals (vdW). As forças de vdW derivam de flutuações quânticas não locais de dipolo-dipolo elétrico instantâneo. Como a LDA e a GGA consideram apenas a densidade ou o gradiente de densidade local, elas não capturam essas flutuações, falhando em prever o espaçamento interplanar em materiais bidimensionais (como grafeno, dissulfeto de molibdênio e nitreto de boro) e a cristalização de redes moleculares e polímeros.

Para corrigir essa falha, desenvolveram-se abordagens empíricas e ab-initio para integrar o efeito vdW na DFT:

• Métodos de Correção Empírica (Grimme - DFT-D2, DFT-D3, DFT-D4): Adicionam um termo de energia de dispersão do tipo $C_6/R^6$ à energia total calculada pela DFT, onde os coeficientes $C_6$ dependem do elemento químico e do seu ambiente de coordenação local. É um método rápido que atinge excelente precisão com baixo custo computacional.
• Funcionais Não Locais de vdW (vdW-DF): Incorporam diretamente no funcional de troca e correlação um termo de integração dupla não local sobre a densidade em dois pontos distantes no espaço, $n(\mathbf{r})$ e $n(\mathbf{r}')$, capturando a resposta de polarizabilidade sem uso de parâmetros empíricos.

Com a adição contínua de correções físicas como o DFT+U, funcionais híbridos e correções de dispersão de van der Waals, a Teoria do Funcional da Densidade consolidou-se como o pilar da simulação de materiais. Ela oferece uma ponte quantitativa confiável entre os postulados fundamentais da mecânica quântica e o comportamento termodinâmico, estrutural e eletrônico de sistemas reais, atuando como o ponto de partida essencial para a engenharia de materiais auxiliada por computador.