Obliczenia energii całkowitej

Baza funkcji falowych

Według twierdzenia Blocha, funkcja falowa elektronu mo˙ze by´c wyra˙zona jako iloczyn fali płaskiej e^ikr o wektorze falowym k oraz pewnej funkcji okresowej uk(r):

ψ_k(r) = e^ikru_k(r) . (2.18) Okres funkcji u_kzwi ˛azany jest z periodyczno´sci ˛a sieci krystalicznej, tzn. u_k(r) = u_k(r + R_l), gdzie Rl jest wektorem tej sieci. Zgodnie z tym, po rozwini˛eciu funkcji okresowej w szereg

Fouriera, funkcj˛e falow ˛a elektronu mo˙zna przedstawi´c w postaci fal płaskich: ψ_jk(r) =^X

G

c_jk(G)e^i(k+G)r, (2.19)

w której G jest wektorem sieci odwrotnej, a współczynniki cjk s ˛a transformatami fourierow-skimi rozwini˛ecia funkcji okresowej.

Taki sposób przedstawienia funkcji falowej znacznie ułatwia wybór funkcji bazowych, które powinny charakteryzowa´c si˛e prostot ˛a i zupełno´sci ˛a. Pozwala na dobór odpowiedniej licz-by funkcji do osi ˛agania zadanego poziomu zbie˙zno´sci i w razie potrzeby rozszerzy´c baz˛e. Poza tym baza fal płaskich ułatwia obliczanie składników macierzy operatora energii kinetycznej oraz nie zale˙zy od poło˙ze´n rdzeni atomów, daj ˛ac wi˛ecej swobody w opisie całego układu. Z drugiej strony, aby uzyska´c jak najlepszy opis funkcji falowej elektronu trzeba zastosowa´c bardzo du˙z ˛a liczb˛e fal płaskich. Zazwyczaj na kształt funkcji falowej elektronu wpływaj ˛a fale płaskie o niskich energiach kinetycznych, natomiast te o wy˙zszych energiach maj ˛a na ni ˛a zni-komy wpływ. Powoduje to, ˙ze w obliczeniach mo˙zna pomin ˛a´c fale rozwini˛ecia powy˙zej pewnej granicznej energii, zwanej energi ˛a odci˛ecia E_cut. Pozwala to na kontrol˛e wielko´sci bazy oraz, w razie potrzeby, dostosowanie jej do konkretnych potrzeb.

Superkomórki

W przypadku obni˙zonej symetrii kryształów, wynikaj ˛acej z obecno´sci powierzchni, de-fektów lub dodatkowych atomów, wyznaczanie funkcji falowych staje si˛e utrudnione. Wygod-nie jest wtedy wykorzysta´c metod˛e superkomórek. Powi˛ekszona komórka krystaliczna umo˙z-liwia zamodelowanie np. defektu, który zostaje otoczony układem regularnej sieci. Gdy taka komórka jest odpowiednio du˙za, po jej powieleniu defekty znajduj ˛ace si˛e w s ˛asiaduj ˛acych su-perkomórkach nie powinny wpływa´c na siebie. Podobna sytuacja ma miejsce w przypadku powierzchni kryształu, która w płaszczy´znie ma układ periodyczny, natomiast w kierunku pro-stopadłym nast˛epuje gwałtowna zmiana jej symetrii. Rozwi ˛azaniem tego jest zamodelowanie superkomórki składaj ˛acej si˛e z wielowarstwowej płytki, przedzielonej obszarem pró˙zni w kie-runku prostopadłym do powierzchni (rys. 2.1). Odpowiednio dobrana grubo´s´c obszaru pró˙zni sprawia, ˙ze po powieleniu superkomórki w przestrzeni, powierzchnie w kolejnych powtarzaj ˛ a-cych si˛e komórkach nie oddziaływuj ˛a na siebie.

obszar pró ni

Rysunek 2.1: Superkomórka modeluj ˛aca powierzchni˛e (001) kryształu bcc oraz jej periodyczne powtó-rzenia.

Metoda punktów specjalnych

Baza fal płaskich pozwala na znaczne uproszczenie równa´n Kohna-Shama. Potencjał efektywny v_ef(r) mo˙zna rozwin ˛a´c w szereg Fouriera:

v_ef(r) =^X

G

v_ef(G)e^iGr. (2.20) Wstawiaj ˛ac (2.19) i (2.20) do równa´n Kohna-Shama (2.13) oraz całkuj ˛ac po obj˛eto´sci komórki elementarnej, otrzymuje si˛e niesko´nczony układ równa´n liniowych na współczynniki c_jk:

X k⁰  1 2|k + G|²δ_GG0 + vef(G, G⁰)  | {z } H_k+G,k+G0 cjk(G⁰) = εjkcjk(G) , (2.21) gdzie H_k+G,k+G0 s ˛a elementami macierzy Hamiltona. Diagonalizacja macierzy pozwala ob-liczy´c współczynniki c_jk. Dobre przybli˙zenie mo˙zna uzyska´c, znajduj ˛ac rozwi ˛azania jedynie dla warto´sci wektora falowego k mieszcz ˛acych si˛e w pierwszej strefie Brillouina. Otrzymane współczynniki c_jkpozwalaj ˛a nast˛epnie, za pomoc ˛a wzoru (2.19), obliczy´c g˛esto´s´c elektronow ˛a. Metoda ta wi ˛a˙ze si˛e z wykorzystaniem du˙zej liczby fal płaskich na ka˙zdy atom. Powoduje to ogromny wzrost ich liczby, w miar˛e zwi˛ekszania rozmiaru superkomórki.

Wyliczenie g˛esto´sci elektronowej wi ˛a˙ze si˛e z obliczeniem ´sredniej wielko´sci funkcji n w przestrzeni k po obj˛eto´sci strefy Brillouina (Ω_BZ):

¯ n = ¹ Ω_BZ Z BZ n(k)dk , (2.22)

Wyliczenie takiej całki jest bardzo skomplikowane, gdy˙z wymagana jest znajomo´s´c funkcji n dla wielu punktów k. Metod ˛a, która znacznie poprawia wydajno´s´c oblicze´n, poprzez znaczn ˛a redukcj˛e punktów k, jest metoda punktów specjalnych.

W metodzie tej, zamiast wykonywa´c całkowanie po wszystkich punktach k w strefie Bruillouina, przeprowadza si˛e sumowanie po pewnych wybranych punktach specjalnych k, któ-rym zostaje przypisana odpowiednia waga:

1 ΩBZ Z BZ n(k)dk −→ ^X k w_kn(k) . (2.23) Punkty specjalne s ˛a wybierane w taki sposób, aby otrzyma´c optymaln ˛a zbie˙zno´s´c funkcji. Jed-na z metod generowania punktów specjalnych została podaJed-na przez Monkhorsta i Packa [77]. Poniewa˙z czas oblicze´n zwi ˛azany jest z wyborem liczby punktów k, stosuje si˛e mo˙zliwie jak najmniejsz ˛a ich liczb˛e, która pozwala osi ˛agn ˛a´c wymagan ˛a zbie˙zno´s´c. Szybko´s´c oblicze´n zale˙zy tak˙ze od typów atomów. W przypadku metali, gdzie pasma elektronowe przebiegaj ˛a przez po-wierzchni˛e Fermiego, uzyskanie zbie˙zno´sci jest wydłu˙zone. Popraw˛e mo˙zna uzyska´c stosuj ˛ac rozmycie powierzchni Fermiego [78].

Energia całkowita układu

Zgodnie z równaniem Schrödingera (2.5), wzór na energi˛e całkowit ˛a przyjmuje posta´c:

Etotal = Ekin+ EH + Eej + Exc+ Ejj, (2.24)

gdzie E_kin jest energi ˛a kinetyczn ˛a nieoddziałuj ˛acych elektronów, E_H – energi ˛a Hartree’ego, która opisuje oddziaływania kulombowskie elektronów, Eej – energi ˛a oddziaływania elektron-jon, E_xc – energi ˛a wymienno-korelacyjn ˛a, natomiast E_jj – energi ˛a Madelunga, opisuj ˛ac ˛a od-działywania jon-jon. Jedynym problemem pojawiaj ˛acym si˛e w obliczeniach energii całkowitej jest znajdowanie funkcjonału energii wymienno-korelacyjnej, które wymaga stosowania przy-bli˙ze´n.

Do obliczenia energii całkowitej wykorzystuje si˛e metod˛e samouzgodnion ˛a. Wymaga ona okre´slenia pocz ˛atkowej g˛esto´sci elektronowej, która słu˙zy do obliczenia potencjałów Har-tree’go oraz wymienno-korelacyjnego. Nast˛epnie, poprzez diagonalizacj˛e macierzy H_k+G,k+G0, rozwi ˛azywane s ˛a równania Kohna-Shama. Otrzymane stany własne słu˙z ˛a do policzenia nowej g˛esto´sci elektronowej. Proces ten jest nast˛epnie powtarzany do momentu otrzymania satysfak-cjonuj ˛acej zbie˙zno´sci.

Powy˙zszy sposób otrzymywania energii całkowitej daje bardzo dokładne warto´sci ener-gii całkowitej. Niestety, ze wzgl˛edu na rosn ˛ac ˛a liczb˛e wymaganych fal płaskich wraz ze zwi˛ek-szaniem si˛e układu (około 100 fal na atom), metoda ta jest wydajna jedynie dla niewielkich układów. Jednym ze sposobów na zwi˛ekszenie wydajno´sci jest ewolucja macierzy Hamiltona podczas wyliczania stanów własnych [79]. Znacznie redukuje to czas potrzebny na osi ˛agni˛ecie odpowiedniego poziomu zbie˙zno´sci.

Rozdział 3

Szczegóły oblicze ´n

Prezentowane w niniejszej pracy wyniki oblicze´n zostały otrzymane przy u˙zyciu pro-gramu VASP (Vienna ab initio simulation package) [80, 81], którego podstaw ˛a jest teoria funk-cjonału g˛esto´sci. Obliczenia te wi ˛a˙z ˛a si˛e z zastosowaniem pewnych przybli˙ze´n oraz dobraniem odpowiednich parametrów. Wymagana w obliczeniach energia wymienno-korelacyjna została opisana za pomoc ˛a przybli˙zenia GGA [70, 82] w spinowo-spolaryzowanej wersji PBE [71]. Oddziaływania elektronowo-jonowe opisane zostały potencjałem PAW [76]. Jako stany walen-cyjne przyj˛ete zostały stany 3d⁷4s¹ dla atomów ˙zelaza, 3s² dla magnezu oraz 2s²2p⁴ dla tlenu. Obsadzenia stanów realizowane były za pomoc ˛a metody Methfessela-Paxtona pierwszego rz˛e-du [83] z rozmyciem powierzchni Fermiego równym 0,2 eV.

3.1 Parametry zastosowane w obliczeniach

Warto´sci parametrów takich jak energia odci˛ecia bazy fal płaskich, liczba punktów spe-cjalnych k, czy grubo´s´c płytki modeluj ˛acej powierzchni˛e ˙zelaza zostały dobrane w wyniku prze-prowadzenia serii testów. Zwi˛ekszanie tych parametrów pozwala na uzyskanie dokładniejszych wyników, lecz wi ˛a˙ze si˛e tak˙ze ze zwi˛ekszeniem kosztu obliczeniowego oraz wydłu˙zeniem cza-su oblicze´n. Dlatego te˙z odpowiednie dobranie tych parametrów jest konieczne, aby zachowa´c rozs ˛adek pomi˛edzy dokładno´sci ˛a oblicze´n oraz czasem potrzebnym do ich uzyskania.

Warto´s´c energii odci˛ecia bazy fal płaskich została ustalona na podstawie serii testów przeprowadzonych dla kryształu ˙zelaza bcc. Podczas nich energia była zwi˛ekszana od 268 eV (warto´sci minimalnej dla atomów Fe) do 700 eV. Rysunki 3.1 i 3.2 przedstawiaj ˛a zmiany energii całkowitej kryształu Fe w przeliczeniu na jeden atom (Etotal) oraz stałej sieci kryształu ˙zelaza

300 400 500 600 700 E cut [eV] -8.318 -8.316 -8.314 -8.312 -8.310 -8.308 -8.306 -8.304 E ^total [eV / atom]

Rysunek 3.1: Zale˙zno´s´c energii całkowitej (E_total) kryształu ˙zelaza bcc od energii odci˛ecia (Ecut).

300 400 500 600 700 E_cut [eV] 2.810 2.815 2.820 2.825 2.830 2.835 a ^Fe [Å]

bcc (aFe) w zale˙zno´sci od energii odci˛ecia (Ecut). Zauwa˙zy´c mo˙zna wyra´zne stabilizowanie si˛e E_total od warto´sci E_cut równej 400 eV. Z kolei warto´s´c stałej sieci, pocz ˛awszy od E_cut wyno-sz ˛acej 300 eV, oscyluje w okolicy 2,830 Å, natomiast od energii równej 500 eV, przyjmuje stał ˛a warto´s´c 2,832 Å. W zwi ˛azku z powy˙zszymi rezultatami, warto´s´c energii odci˛ecia zastosowa-na w obliczeniach wła´sciwo´sci kryształu oraz powierzchni ˙zelaza, a tak˙ze adsorpcji cz ˛asteczek MgO oraz wzrostu FeO na powierzchni Fe(001) została ustalona na 500 eV. Ze wzgl˛edu na wy-dajno´s´c, obliczenia adsorpcji atomów O i Mg na powierzchni Fe(001) zostały przeprowadzone dla E_cut=400 eV.

Siatka punktów specjalnych k została wygenerowana metod ˛a Monkhorsta-Packa [77]. Rozmiar siatki równie˙z został dobrany poprzez testowe obliczenia kryształu ˙zelaza. Na rysun-ku 3.3 wykre´slona została energia całkowita kryształu w funkcji rozmiaru siatki punktów k. Na wykresie zauwa˙zy´c mo˙zna malej ˛ace oscylacje E_total wraz ze wzrostem liczby punktów k. Po-zwoliło to ustali´c siatk˛e punktów specjalnych wykorzystan ˛a w dalszych obliczeniach w zale˙zno-´sci od rozmiaru komórki zastosowanej w obliczeniach. Dla oblicze´n litego kryształu ˙zelaza, siat-ka punktów k została ustalona na 12×12×12. Taki dobór siatki pozwolił potem na odpowiednie przeskalowanie liczby punktów k dla wi˛ekszych komórek powierzchniowych. W obliczeniach powierzchni Fe(001) w komórce powierzchniowej 1×1 zastosowano siatk˛e 12×12×1. Do ob-licze´n adsorpcji zastosowana została siatka 6×6×1 dla komórki 2×2 Fe(001) oraz 4×4×1 dla komórki 3×3 Fe(001).

Powierzchnia ˙zelaza Fe(001) została zamodelowana za pomoc ˛a płytki składaj ˛acej si˛e

4 6 8 10 12 14 16 18 Liczba punktów k -8.325 -8.320 -8.315 -8.310 -8.305 -8.300 -8.295 E ^total [eV / atom]

Rysunek 3.3: Zmiany energii całkowitej (Etotal) kryształu ˙zelaza bcc w funkcji rozmiaru siatki punktów specjalnych k.

z warstw ˙zelaza oraz grubej warstwy pró˙zni, wynosz ˛acej 17 Å. Wybór grubo´sci płytki został poprzedzony kolejn ˛a seri ˛a testów, w których liczba warstw Fe była zwi˛ekszana z 5 do 16 mo-nowarstw. Płytka relaksowana była obustronnie. Pozycje atomów w komórce były optymalizo-wane do momentu, w którym siły działaj ˛ace na ka˙zdy atom były mniejsze ni˙z 0,01 eV/Å. Na rysunku 3.4 przedstawione zostały zmiany relaksacji powierzchni Fe(001) w zale˙zno´sci od gru-bo´sci płytki. Relaksacje zostały zdefiniowane jako zmiany odległo´sci pomi˛edzy i-t ˛a oraz j-t ˛a warstw ˛a ˙zelaza w płytce (dij) w stosunku do odległo´sci mi˛edzypłaszczyznowej w krysztale Fe (d) według nast˛epuj ˛acego wzoru:

∆_ij = ^dij − d

d ^{. (3.1)}

Z rysunku zauwa˙zy´c mo˙zna stabilizacj˛e relaksacji rozpoczynaj ˛ac ˛a si˛e dla płytek o grubo´sciach 9–10 monowarstw. Ze wzgl˛edu na symetri˛e, dalsze obliczenia były wykonywane dla płytki o grubo´sci 9 monowarstw ˙zelaza. Pozwala to na uzyskanie dokładnych wła´sciwo´sci fizycznych powierzchni Fe(001) oraz jest wystarczaj ˛aca do oblicze´n adsorpcji.

Cz˛e´s´c oblicze´n została przeprowadzona dla powierzchni Fe(001) relaksowanej jedno-stronnie. W takim przypadku zastosowana została poprawka dipolowa, której celem jest usu-ni˛ecie asymetrii pola elektrycznego powstaj ˛acej z drugiej strony płytki.

6 8 10 12 14 16 Liczba warstw -3 -2 -1 0 1 2 3 4 ∆ [%] ∆₁₂ ∆₂₃ ∆₃₄