4. Obliczenia struktury pasmowej
4.1. Podstawy teoretyczne obliczeń DFT
Ten podrozdział został opracowany na podstawie wprowadzenia do obliczeń w pakiecie WIEN2k [Cottenier 2002] i instrukcji użytkownika pakietu WIEN2k [Blacha 2012].
Ciało stałe to związany układ masywnych i naładowanych dodatnio jąder atomowych oraz wielokrotniemniej masywnych i ujemnie naładowanych elektronów. Dla układu N jąder atomowych liczba cząstek oddziałowujących elektromagnetycznie wynosi N+ZN, gdzie Z to ładunek jądra atomowego. Hamiltonian takiego układu wielu oddziaływujących cząsteczek ma postać: ̂ ∑ ⃗⃗ ∑ ⃗⃗ ∑ | ⃗ | ∑ | | ∑ | ⃗ ⃗ | (4.1)
gdzie to masa jądra atomowego w położeniu ⃗ , to masa elektronu w położeniu . Pierwszy i drugi składnik hamiltonianu to operatory energii kinetycznej jąder atomowych i elektronów. Trzy kolejne składniki opisują oddziaływanie kulombowskie między elektronami i jądrami, elektronami i elektronami (elektron-elektron), jądrami – jądrami.
Rozwiązania równania Schrödingera dla hamiltonianu (4.1) nie są znane, jednak po wprowadzaniu pewnych uproszczeń możliwe jest znalezienie rozwiązań przybliżonych. Pierwsze przybliżenie jest skutkiem dużej różnicy mas elektronów i jąder atomowych (trzy rzędy wielkości i więcej). Z punktu widzenia elektronów masywne jądra atomowe nie poruszają się, zatem elektrony zawsze są w stanie równowagi termodynamicznej z powolnymi jądrami atomowymi. W pierwszym przybliżeniu, nazywanym przybliżeniem Borna-Oppenheimera (B-O) [Born 1927], zakłada się, że jądra atomowe to naładowane cząstki, których położenie jest niezmienne w czasie. Konsekwencją przybliżenia B-O jest uproszczenie hamiltonianu (4.1) do sumy trzech składników:
̂ ̂ ̂ ̂ (4.2)
gdzie ̂ to operator energii kinetycznej elektronów, ̂ to operator oddziaływania elektron-elektron a ̂ operator opisujący potencjał elektryczny od jąder atomowych (ładunek jądra atomowego Z i położenie jądra).
Kolejnym krokiem do znalezienia rozwiązań jest zastosowanie dwóch twierdzeń Hohenberga i Kohna [Hohenberg 1964].
Twierdzenie I. Istnieje odpowiedniość między gęstością elektronową w stanie podstawowym układu wielo-elektronowego (atomu, cząsteczki, ciała stałego) ( ) i potencjałem , z czego bezpośrednio wynika, że wartość oczekiwana każdego operatora ̂ dla układu w stanie podstawowym jest funkcjonałem gęstości elektronowej w stanie podstawowym
⟨ | ̂| ⟩ [ ] (4.3)
Twierdzenie II. Dla obserwabli ̂ będącej hamiltonianem ̂, funkcjonał energii całkowitej w stanie podstawowym [ ] [ ] ma postać:
102
[ ] ⟨ | ̂ ̂| ⟩⏟
[ ]
⟨ | ̂ | ⟩ [ ] ∫ ( ) ( ) (4.4) gdzie [ ] to funkcjonał gęstości elektronowej Hohenberga-Kohna, który jest uniwersalny dla każdego układu wieloelektronowego. [ ] osiąga minimum równe energii całkowitej w stanie podstawowym, dla gęstości elektronowej w stanie podstawowym odpowiadającej
.
Z twierdzenia I wynika, że rozkład gęstości elektronowej dla układu w stanie podstawowym zawiera dokładnie tyle samo informacji o układzie, co jego funkcja falowa. Na podstawie twierdzenia II wiadomo, że funkcjonał [ ] jest uniwersalny dla każdego układu wieloelektronowego, ponieważ nie zawiera on żadnych informacji o jądrach atomowych (np. położenia jąder atomowych), a znajomość rozkładu gęstości elektronowej dla układu w stanie podstawowym pozwala na obliczenie wkładu zewnętrznego potencjału do energii całkowitej układu. Ponadto według twierdzenia II dla nieskończonego zbioru możliwych gęstości elektronowych, istnieje gęstość elektronowa, która minimalizuje funkcjonał energii i jest gęstością elektronową układu w stanie podstawowym odpowiadającą zewnętrznemu potencjałowi ( ).
Na podstawie pracy [Kohn 1965] funkcjonał gęstości elektronowej Hohenberga-Kohna [ ] można przedstawić, jako sumę trzech funkcjonałów:
[ ] [ ] [ ] ( ⏟ [ ] [ ])
[ ]
, (4.5)
gdzie to funkcjonał energii kinetycznej nieoddziaływującego gazu elektronowego, to funkcjonał Hartreego, to funkcjonał wymienny, to funkcjonał korelacyjny, to funkcjonał energii korelacyjno-wymiennej. Zatem funkcjonał energii całkowitej dla układu w stanie podstawowym [ ] ma postać:
[ ] [ ] [ ] [ ] [ ] . (4.6)
Powyższe równanie można interpretować, jako funkcjonał energii nieoddziaływującego gazu elektronów, znajdującego się w dwóch zewnętrznych potencjałach: od jąder atomowych i potencjale związanym z efektami korelacyjno-wymiennymi. Hamiltonian dla takiego gazu został nazwany Hamiltonianem Kohna-Shama. Ma on postać:
̂ ̂ ̂ ̂ [ ] ̂ ⃗⃗ ∫| | ( ) , (4.7) gdzie potencjał korelacyjno-wymienny dany jest przez pochodną funkcjonału:
̂ [ ]
. (4.8)
Według teorii Kohna Shama [Kohn 1965] gęstość elektronowa dla stanu podstawowego N-elektronowego układu wynosi:
( ) ∑ ( ) ( ) (4.9) gdzie ( ) to jednocząstkowe funkcje falowe.
Zbiór N jednocząstkowych funkcji falowych ( ) jest najniższym-energetycznie rozwiązaniem równania Kohna-Shama (K-S):
103 Jednocząstkowe funkcji falowe ( ) nie opisują elektronów, tylko kwasi-cząstki, które nie mają fizycznego odpowiednika. Całkowita gęstość wszystkich kwasi cząsteczek jest równa gęstości elektronowej układu. Równanie Kohna Shama (4.10) jest odpowiednikiem równania Schrödingera dla układu nieoddziaływujących kwasi cząstek. Jest ono dużo łatwiejsze do rozwiązania niż równanie Schrödingera dla układu wielu oddziaływujących elektronów. Rozwiązanie równania K-S pozwala znaleźć gęstość elektronową dla układu w stanie podstawowym. Operatory i zależą od gęstości elektronowej ( ), która sama zależy od jednocząstkowych funkcji falowych , które z kolei są poszukiwanym rozwiązaniem równania K-S. Zatem mamy tutaj do czynienia z problemem samouzgodnionym, ponieważ rozwiązania określają postać rozwiązywanego równania K-S ( i w zależą od postaci ). Aby rozwiązać równanie samouzgodnione stosuje się procedurę przedstawioną na Rys. 4.1. W pierwszym kroku przyjmuje się arbitralnie pewną początkową gęstość elektronową i na jej podstawie konstruuje się hamiltonian , (w tym wyrażenia na i ). Następnie rozwiązuje się równanie K-S otrzymując, jako rozwiązania wartości własne i jednocząstkowe funkcje falowe , na podstawie, których tworzy się gęstość elektronową . Obliczona gęstość elektronowa służy do stworzenia nowego (poprawionego) hamiltonianu i rozwiązania poprawionego równania K-S, którego rozwiązaniem są inne (poprawione) funkcje falowe , na podstawie, których zostanie utworzona nowa gęstość elektronowa , która zostanie wykorzystana w kolejnej iteracji. Powyższą procedurę przeprowadza się do momentu, gdy obliczona gęstość elektronowa w n-tym kroku stanie się zgodna z gęstością obliczoną w kroku wcześniejszym (n-1). W pakiecie WIEN2k za kryterium zgodności (zbieżności) przyjęto zgodność gęstości elektronowych dla trzech kolejnych kroków iteracji [Blacha 2012].
Rys. 4.1. Przebieg obliczeń samouzgodnionych
stosowany przy rozwiązywaniu równania Kohna-Shama. Obliczenia kończą się w momencie, gdy obliczona gęstość jest równa w pewnym przybliżeniu gęstości [Cottenier 2002].
104 Drugim przybliżeniem, jakie należy zastosować, aby rozwiązać problem wieloelektronowego układu oddziaływujących cząstek jest założenie postaci funkcjonału wymienno-korelacyjnego. Najprostszym przybliżeniem jest przybliżenie zwane LSDA (ang. Local Spin Density Approximation), w którym zakłada się funkcjonał korelacyjno-wymienny w postaci:
∫ ( ) ( ( )) , (4.11)
gdzie ( ) to funkcja korelacyjno-wymienna jednorodnego gazu elektronowego, której postać znana jest z obliczeń numerycznych.
W przybliżeniu LSDA zakłada się, że energia korelacyjno-wymienna odpowiadająca rozkładowi gęstości elektronowej ( ) może zostać znaleziona przez podzielenie objętości materiału na nieskończenie małe objętości o stałej gęstości elektronowej. Każda z tych objętości wnosi wkład do całkowitej energii korelacyjno-wymiennej równy wkładowi energii korelacyjno-wymiennej od identycznej objętości wypełnionej jednorodnym gazem elektronowym, którego gęstość ładunkowa jest taka sama jak gęstość ładunkowa materiału w tej objętości.
Rozwinięciem LSDA jest przybliżenie uogólnionego gradientu GGA (ang. Generalized Gradient Approximation), w którym dodatkowo uwzględnia się wkład do całkowitej energii korelacyjno-wymiennej pochodzący od gradientu gęstości ładunkowej. W przybliżeniu GGA funkcjonał korelacyjno-wymienny ma postać:
∫ ( ) ( ( ) ( ( ))) . (4.12)
W odróżnieniu od przybliżenia LSDA, przybliżenie GGA pozwala na wprowadzenie różnych parametrów opisujących gradient gęstości ładunkowej. Możliwość parametryzacji zaowocowała powstaniem licznych odmian przybliżenia GGA. Sposób parametryzacji i wartości parametrów dobiera się na podstawie analizy dużej liczby danych doświadczalnych. Powstały także odmiany GGA, które starają się odtworzyć funkcjonał korelacyjno-wymienny z pierwszych zasad i są wolne od parametrów [Cottenier 2002,
Perdew 1996]. Przy przeprowadzeniu obliczeń struktury elektronowo-pasmowej LaPbMnFeO wybrano bezparametrowe przybliżenie PBEsol-GGA [Perdew 2008].
Ostatecznie po omówionych przybliżeniach otrzymuje się zestaw n równań (odpowiednik do równania Schrödingera), którego rozwiązaniem jest zbiór n funkcji falowych :
( ⃗⃗ ∫| | ( ) )
⏟
̂
( ) ( ) , (4.13)
gdzie ̂ to jednocząstkowy hamiltonian, to funkcja falowa n-tej kwasi-cząstki a to operator wymienno-korelacyjny w przybliżeniu LSDA, GGA lub innym.
Trzecie przybliżenie dotyczy postaci jaką mają rozwiązania równania (4.13). Przyjmuje się, że każda z funkcji falowych , będąca rozwiązaniem (4.13) może być przedstawiona jako kombinacja liniowa (p- licznej) bazy funkcji (b – oznacza inną bazę) i zbioru współczynników według poniższego równania:
105 Baza funkcji może być bazą nieskończoną. Jednak liczbę i rodzaj funkcji dobiera się tak, aby jak najdokładniej opisywała funkcje falowe i zarazem była możliwe mała. Dokładność i czas obliczeń, a zatem stopień przybliżenia, rośnie wraz z p.
Podsumowując. Korzystając z wymienionych założeń został wyprowadzony odpowiednik równania Schrödingera (4.13), który opisuje właściwości układu wielo-elektronowego. Rozwiązaniem układu n- równań (4.13) są odpowiednie rozkłady gęstości elektronowych dla układu w stanie podstawowym pod warunkami: założenia racjonalnych postaci operatorów wymienno-korelacyjnego oraz potencjału jądra atomowego i wybraniu bazy funkcji własnych. Wszystkie właściwości układu wielu oddziałujących elektronów można obliczyć o ile znany jest rozkład gęstości elektronowych.