Warsztaty modelowania w nanofizyce

Podstawy obliczeń ab initio

Przemysław Piekarz

Zakład Komputerowych Badań Materiałów

Instytut Fizyki Jądrowej PAN, Kraków

Metody ab initio

Metody obliczeniowe oparte na prawach mechaniki kwantowej, stosowane do wyznaczania i badania podstawowych własności materiałowych takich jak:

• parametry sieci krystalicznej
• struktura elektronowa
• własności magnetyczne
• dynamika sieci
• stałe elastyczne
• własności optyczne

Oddziaływania w krysztale

Wiązania międzyatomowe, dzięki którym powstają układy molekularne, ciała stałe i nanostruktury, są efektem wzajemnych oddziaływań elektrostatycznych w układzie elektronów i jader atomowych. Funkcja falowa elektronów i jąder atomowych $\Psi(\vec{r}_i,\vec{R}_j)$ spełnia równanie Schrodingera

$ H\Psi(\vec{r}_i,\vec{R}_j) = E\Psi(\vec{r}_i,\vec{R}_j), $

gdzie $H$ jest Hamiltonianem układu, a $E$ jego energią.

$ H=-\frac{\hbar^2}{2m}\sum_{i}\nabla_i^2+V_{ext}(\vec{r})+\frac{1}{2}\sum_{i\neq j}\frac{e^2}{|\vec{r}_i-\vec{r}_j|} $

Można rozwiązać ściśle tylko dla szczególnych przypadków (np. atom wodoru).

Teoria funkcjonału gęstości (DFT)

P. Hohenberg and W. Kohn, Phys. Rev. 136, 864 (1964)

W teorii funkcjonału gęstości podstawową wielkością jest gęstość elektronowa określona w każdym punkcie kryształu $n(r)$.

Energia całkowita układu jest funkcjonałem gęstości:

$E[n]=T[n]+\int dr V_{ext}(\vec{r}) n(r) + E_H[n] + E_{xc}[n]$

Teoria DFT opiera się na dwóch twierdzeniach Hohenberga-Kohna.

Twierdzenia Hohenberga-Kohna

Twierdzenie I

Zewnętrzny potencjał układu oddziałujących elektronów $V_{ext}(\vec{r})$ jest jednoznacznie określony przez gęstość elektronową w stanie podstawowym $n_0(\vec{r})$.

Twierdzenie II

Dla ustalonego potencjału zewnętrznego $V_{ext}(\vec{r})$, funkcjonał energii $E[n]$ osiąga minimalną wartość $E_0$ dla gęstości elektronowej w stanie podstawowym $n_0(\vec{r})$.

$E[n_0] = E_0 = min$

Równanie Kohna-Shama

W. Kohn and L. J. Sham, Phys. Rev. 140, 1133 (1965)

$(-\frac{\hbar^2\nabla_i^2}{2m}+V_{KS})\psi_i(\vec{r})=\varepsilon_i\psi_i(\vec{r})$

$V_{KS}=V_{ext}+V_H+V_{xc}=V_{ext}+\int dr' \frac{n(r')}{|\vec{r}-\vec{r}'|}+\frac{\delta E_{xc}}{\delta n}$

$n(\vec{r})=\sum_{i,\sigma} |\psi_i^{\sigma}(\vec{r})|^2$

Energia wymienno-korelacyjna $E_{xc}$ nie jest znana dokładnie i wyliczana jest w przybliżeniu LDA lub GGA.

Przybliżenie lokalnej gęstości (LDA)

W przybliżeniu LDA energia wymienno-korelacyjna w każdym punkcie przestrzeni $\vec{r}$, gdzie gestość wynosi $n(\vec{r})$, równa jest energii wymienno-korelacyjnej jednorodnego gazu elektronowego o tej samej gęstości, $n=n(\vec{r})$:

$E_{xc}[n]=\int dr n(\vec{r}) \epsilon_{xc}(n),$

$\epsilon_{xc}(n)=\epsilon_x(n)+\epsilon_c(n)$

Gestość energii wymiennej wyliczana jest ze wzoru: $\epsilon_x(n)=-\frac{3(3\pi)^\frac{1}{3}}{4\pi}n^\frac{4}{3}$

Gęstość energii korelacyjnej $\epsilon_c(n)$ wyznaczana jest metodą kwantowego Monte Carlo:

D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980)

Przybliżenie uogólnionych gradientów (GGA)

W przybliżeniu uogólnionych gradientów (generalized gradient approximation - GGA), energia wymienno-korelacyjna jest funkcjonałem gęstości elektronowej i jej gradientów:

$ E_{xc}[n^{\uparrow},n^{\downarrow}]=\int dr f(n^{\uparrow},n^{\downarrow},\nabla n^{\uparrow},\nabla n^{\downarrow}) $

Funkcjonał GGA-PBE

J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996)

Pseudopotencjał

Przybliżony potencjał działający na elektrony walencyjne w obszarze rdzenia atomowego nazywany pseudopotencjałem. Pseudofunkcja falowa w obszarze międzywęzłowym jest równa dokładnej funcji falowej, a w obszarze rdzenia jest wolno-zmienną funkcją bez oscylacji i miejsc zerowych.

M. C. Payne et al., Reviews of Modern Physics 64, 1045 (1992)

Rodzaje pseudopotencjałów

• Pseudopotencjały zachowujące normę $\int_0^r dr r^2 |\psi(r)|^2=\int_0^r dr r^2 |\tilde{\psi}(r)|^2$

• Pseudopotencjały ultramiękkie (ultrasoft)

• PAW (projector augmented-wave)

Sieć krystaliczna

Przykład: tytanian baru BaTiO$_3$ - sieć tetragonalna

Parametry: stałe sieci $a$ i $c$, położenia atomów Ti i O

Nanostruktury

Przykład: monowarstwa Fe na powierzchni BaTiO$_3(001)$

Parametry: stała sieci $a$ i położenia atomów w każdej warstwie

Schemat iteracyjnej procedury optymalizacyjnej

Parametry procedury optymalizacyjnej

• Maksymalna energia kinetyczna $E_{cut}=\frac{\hbar^2 G_{cut}^2}{2m}$ definiuje ilość funkcji bazowych (fal płaskich):

  $\psi_\vec{k}(\vec{r})=\sum_\vec{G} c_{i,\vec{k}+\vec{G}}e^{ir(\vec{k}+\vec{G})}$

  $|\vec{k}+\vec{G}|<\vec{G}_{cut}$

• Siatka punktów w przestrzeni odwrotnej $k_x \times k_y \times k_z$ (np. $4\times4\times4$)

• Dokładność obliczeń $\Delta E = E_{i+1}- E_{i} < 10^{-7}$ eV

Energia całkowita

$E_{tot}(a_0)=min$

Stan równowagi

Stan podstawowy otrzymany w wyniku procedury optymalizacyjnej powinien spełniać warunki stanu równowagi badanego materiału.

• Energia całkowita przyjmuje minimum: $E_{tot}$ = min

• Całkowita siła działająca na każdy atom zeruje się: $\vec{F}_i=-\frac{\partial E_{tot}}{\partial{\vec{R}_i}}=0$

• Makroskopowe naprężenie w układzie równe jest naprężeniu wywołanym ciśnieniem zewnętrznym $\sigma_{\alpha\beta}=-\frac{1}{V}\frac{\partial E}{\partial u_{\alpha\beta}}$