VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ BRNO UNIVERSITY OF TECHNOLOGY
FAKULTA STROJNÍHO INŽENÝRSTVÍ ÚSTAV FYZIKÁLNÍHO INŽENÝRSTVÍ FACULTY OF MECHANICAL ENGINEERING INSTITUTE OF PHYSICAL ENGINEERING
PRVOPRINCIPIÁLNÍ STUDIUM STABILITY KRYSTALŮ PEVNÝCH LÁTEK FIRST-PRINCIPLES STUDY OF STABILITY OF SOLID CRYSTALS
BAKALÁŘSKÁ PRÁCE BACHELOR´S THESIS
AUTOR PRÁCE
ONDŘEJ PLESKOT
AUTHOR
VEDOUCÍ PRÁCE SUPERVISOR
BRNO 2014
doc. Mgr. MIROSLAV ČERNÝ, Ph.D.
Abstrakt Práce se zabývá modelováním elektronové struktury v krystalech pevných látek. Konkrétně jsou prvoprincipiálním kódem VASP spočteny hustoty valenčních elektronů a hustoty stavů tří krystalů, reprezentující tři typy vazeb. Jedná se o vazbu kovovou v krystalu hliníku, o vazbu kovalentní v krystalu diamantu a o vazbu iontovou v krystalu soli kamenné. Z vypočtených hodnot napětí a energie jsou následně určeny některé makroskopické parametry krystalů, jako jsou rovnovážný mřížkový parametr, modul objemové pružnosti a mez pevnosti při izotropním tahu.
Abstract This work deals with a modeling of an electronic structure of solid crystals. Specifically, densities of valence electrons, and the density of states are calculated using first-principles code VASP for three crystals representing three different types of bonding. These are the metal bonding in a crystal of aluminum, the covalent bonding in a diamond crystal, and an ionic bonding in the crystal of rock salt. From calculated values of stress and crystal energy are then determined some macroscopic parameters of crystals, such as the equilibrium lattice parameter, the bulk modulus and the theoretical strength under isotropic tension.
1
Čestné prohlášení Prohlašuji tímto, že jsem překládanou práci zpracoval samostatně na základě uvedené literatury, konzultací s odborníky a vlastních poznatků. V Brně dne
.........….............. Ondřej Pleskot
2
Poděkování Chtěl bych velmi poděkovat doc. Mgr. Miroslavu Černému, Ph.D. za jeho ochotu mi vždy dobře poradit, zejména v době seznamování se s programem VASP.
3
Obsah 1.
Úvod....................................................................................................................... 5
2.
Ab initio metody .................................................................................................... 6
3.
4.
5.
2.1.
Adiabatická aproximace ................................................................................. 6
2.2.
Teorie funkcionálu hustoty (DFT) .................................................................. 6
2.3.
Pseudopotenciály ............................................................................................ 7
2.4.
Výměnná a korelační energie ......................................................................... 7
VASP ..................................................................................................................... 8 3.1.
Historie ........................................................................................................... 8
3.2.
Základní vstupní soubory ............................................................................... 8
3.3.
Soubory používané programem VASP ........................................................ 10
3.4.
Postup práce s programem VASP ................................................................ 11
Výsledky výpočtů ................................................................................................ 12 4.1.
Vyhodnocování výsledků ............................................................................. 12
4.2.
Hliník ............................................................................................................ 13
4.3.
Diamant ........................................................................................................ 16
4.4.
NaCl .............................................................................................................. 20
Závěr .................................................................................................................... 24
Bibliografie ..................................................................................................................... 25
4
1. Úvod S vývojem výpočetní techniky se také rozšiřují možnosti využití počítačového modelování. Čím je výpočetní technika výkonnější a levnější, tím se více vyplácí začlenit modelování do nejrůznějších odvětví vědy a průmyslu. Jedním z těchto odvětví je také materiálový výzkum a fyzika pevných látek, kde se využívá spojitost mezi vlastnostmi materiálů a jejich elektronovou strukturou. To vedlo rozvoji ab initio výpočtů, které nejsou ovlivněny chybami měření při získávání vstupních dat pro výpočty, protože počítají mechanické vlastnosti jen na základě znalosti elektronové struktury. Tyto výpočty byly velmi náročné, a proto pro jejich rozvoj byl nutný vývoj výpočetních metod, ale zároveň i rozvoj výpočetní techniky. Přestože šel vývoj informačních technologií velmi rychle kupředu, tak se stejně muselo do modelů začlenit velké množství aproximací, což vedlo k rozvoji teorie funkcionálu hustoty, za kterou získal v roce 1998 prof. W. Kohn Nobelovu cenu [1]. To vše vedlo k možnosti provádět ab initio výpočty s velkou přesností, která dosáhla až atomového měřítka, kterého se experimentálně dosáhnout nedá. Pro tyto výpočty vzniklo mnoho programů, které se liší jak použitými metodami výpočtů, tak i použitými aproximacemi. Jedním z takovýchto programů je i program VASP, neboli Vienna Ab-initio Simulation Package. Jedná se o simulační program založený na teorii funkcionálu hustoty. Cílem práce bylo vymodelovat vybrané krystaly pevných látek a zvolenou ab initio metodou spočítat jejich elektronovou strukturu a vlastnosti. Konkrétně bylo mým úkolem seznámit se s programem VASP a pak s jeho pomocí spočítat pro tři vybrané krystaly zastupující tři základní vazby v pevných látkách, jejich rovnovážné mřížkové parametry, objemové moduly pružnosti, hustoty náboje v daných rovinách a meze pevnosti při izotropním tahu. Vybral jsem si krystal hliníku, jako typického zástupce krystalů s kovovou vazbou, krystal diamantu, jako typického zástupce krystalů s kovalentní vazbou a nakonec krystal soli kamenné, jako typického zástupce krystalů s iontovou vazbou. Práce dále obsahuje stručný popis programu VASP.
5
2. Ab initio metody Ab initio znamená prvoprincipiální, neboli vycházející z prvních principů, tedy z obecné charakteristiky látky. Účelem těchto metod je spočítat elektronovou strukturu ze znalosti krystalové struktury a protonových čísel. To nám následně umožní spočítat mechanické a magnetické vlastnosti, míru stability různých atomových uspořádání a rovnovážné pozice atomů. Výhodou těchto metod je, že nejsou zatíženy chybami spojenými se získáváním vstupních dat. První pokusy o takovéto výpočty probíhali již ve 30. letech. V 50. letech [2] začaly pokusy o výpočty pásové struktury, ale problémem byla velká náročnost těchto výpočtů. K výraznému urychlení došlo až s nástupem počítačů. Společně s vývojem počítačů šel i vývoj těchto metod, kdy s objevem nových výpočetních metod se umožnilo začít s výpočty složitějších systémů a to vedlo k pochopení nových souvislostí mezi elektronovou strukturou a vlastnostmi krystalů. Klíčovým momentem ve vývoji a využívání těchto metod přišel v devadesátých letech, kdy vývoj počítačů umožnil přenést řadu výpočtů z velkých superpočítačů na osobní počítače. Cílem výpočtů je spočítat stacionární stav elektronů pohybujících se v poli atomových jader, abychom mohli poté spočítat pohyb jader, kde jako potenciální energii bereme energii základního stavu elektronů. Ovšem v pevných látkách běžných rozměrů je velké množství interagujících částic a je nereálné řešit Schrodingerovu rovnici pro všechny částice. Proto je nutné zavést celou řadu aproximací.
2.1. Adiabatická aproximace Adiabatická aproximace, neboli Bornova-Oppenheimerova aproximace uvažuje oddělený pohyb elektronů a jader. Protože jsou jádra o tři až čtyři řády těžší než elektrony, pohybují se oproti nim výrazně pomaleji. Můžeme proto předpokládat, že se elektrony pohybují tak rychle, že jsou v základním stavu v poli jader. Pomalý pohyb jader pak představuje pro systém elektronů adiabatickou poruchu. To umožňuje studovat prvně pohyb systému elektronů v poli nehybných jader a až potom řešit problém pohybu jader v potenciálovém poli elektronů. Nevýhodou této aproximace je nemožnost popisu jevů souvisejících se vzájemnou interakcí jader a elektronů jako například rezistivitu či supravodivost.
2.2. Teorie funkcionálu hustoty (DFT) Teorie funkcionálu hustoty umožňuje určit všechny fyzikální vlastnosti základního stavu pevných látek, jako jsou celková energie, atd. na základě znalosti elektronové hustoty. Původně byl tento formalismus vyvinut pro popis chování klasických kapalin, ale Thomas a Fermi použili tyto statistické úvahy k aproximaci elektronového příspěvku [3]. Thomasova-Fermiho aproximace předpokládá, že se potenciál působící na elektrony mění dostatečně pomalu, 6
takže kinetická energie je rovna kinetické energii homogenního plynu volných elektronů o stejné hustotě, jaká se jeví lokálně. Celková energie elektronového systému je pak dána pomocí jednočásticové hustoty. Tato teorie měla dobré výsledky pro vysokou hustotu, ale pro normální hustoty byla nevyhovující. Například nepředpovídala vazbu mezi atomy (energie neměla minimum). V roce 1964 Hohenberg a Kohn [4] formulovali dva teorémy, jež umožnili, aby jednočásticová hustota postačovala k popisu základního stavu. Existenční teorém říká, že ze známé elektronové hustoty základního stavu lze určit (až na konstantu) vnější potenciál vyvolaný souborem jader či iontů. Tento potenciál následně určuje mnohočásticový hamiltonián a tedy znalost jednočásticové hustoty určuje úplný hamiltonián. Jakmile je znám hamiltonián, mohou být implicitně určeny všechny vlastnosti základního stavu systému. Variační princip se týká minimalizace celkové energie mnohoelektronového sytému změnami elektronové hustoty základního stavu. Vyhneme se tak nutnosti řešit komplikovanou Schrödingerovu rovnici a postačíme si pouze se změnami jednočásticové hustoty. Tyto dva teorémy nám tak umožní zjistit vlastnosti základního stavu jen za pomocí jednočásticové hustoty.
2.3. Pseudopotenciály Kdybychom chtěli pomocí jednoho bázového souboru specifikovat všechny elektronové stavy systému, pak bychom došli k závěru, že počet rovinných vln, které jsou k tomuto vyčerpávajícímu popisu potřeba, přesahuje výpočetní možnosti. Z tohoto důvodu se při výpočtech používá pseudopotenciálů, které mají za úkol radikálně snížit počet rovinných vln potřebných k výpočtu a zároveň zachovat vysokou přesnost výpočtu. Vnitřní elektrony, jejichž vlnové funkce jsou lokalizovány blízko jádra atomu, jsou velmi málo citlivé na změny v okolí atomu, proto je lze z řešení kvantově mechanických rovnic vyjmout a odpovídající část průběhu Coulombovského potenciálu lze nahradit pomocí pseudopotenciálu. Zároveň je pseudopotenciál nastaven tak, aby neměnil vlnové funkce elektronu ve vnějším elektronovém obalu, čímž je zachována přesnost výpočtu.
2.4. Výměnná a korelační energie Tvar funkcionálu výměnné a korelační energie není znám analyticky, a proto je nutné jej aproximovat. Dvě často využívané aproximace jsou LDA a GGA. LDA (local density approximation) je jednodušší a je založená na tom, že výměnná a korelační energie v každém bodě prostoru nahrazena energií odpovídající elektronovému plynu stejné hustoty. GGA (generalized gradient approximation) je zahrnuje rovněž korekce uvažující i gradient elektronové hustoty. Tato aproimace má více různých parametrizací, jako je například model PBE (Perdew-Burke-Ernzerhof) [5], který využívám při výpočtech programem VASP. 7
3. VASP VASP neboli Vienna Ab-initio Simulation Package je program pro simulace za použití ab-initio výpočtů a pseudopotenciálů. Vytvořili ho především Georg Kresse, Martijn Marsman a Jürgen Hafner z Vídeňské univerzity. Historie, princip fungování a jednotlivé soubory programu VASP jsou velmi důkladně popsané v manuálu [6] [7].
3.1. Historie Práce na programu VASP začali v září roku 1991 a vycházeli z verze programu CASTEP z roku 1989 od Mikea Paynea. V říjnu roku 1992 byly do kódu zahrnuty pseudopotenciály a self-konzistentní smyčka, což umožnilo efektivně používat VASP pro výpočty elektronové struktury kovů. V lednu roku 1993 se ke skupině připojil J. Furthmüller, který se zabýval kódy pro práci s hustotou nábojů, souborem INCAR a podobně. V únoru roku 1995 J. Furthmüller Vídeň opouští. Časem dostal VASP své finální jméno, stal se stabilním a univerzálním nástrojem k ab initio výpočtům. V září roku 1996 se začalo pracovat na paralelizaci kódu a J. M. Holender začlenil do VASPu komunikační jádro programu CETEP. Toto zkopírování zdrojového kódu způsobilo nemalé potíže. Většinu práce na paralelizaci provedl Georg Kresse a ukončil ji okolo ledna 1997. Následně pak kolem července roku 1998 bylo komunikační jádro kompletně přepsáno a tím byly odstraněny všechny části kódu pocházející z programu CETEP. V té době také probíhaly práce na začlenění metody PAW [8] (projector augmented wave), které byly dokončeny v prosinci roku 1999. Roku 2004 začal vývoj nové verze programu zahrnující Hartree a Fockovu aproximaci a teorii lineární odezvy.
3.2. Základní vstupní soubory Pro spuštění programu VASP jsou potřeba 4 základní soubory, jedná se o soubory INCAR, POTCAR, POSCAR a KPOINTS. Soubor INCAR je hlavní vstupní soubor programu VASP a určuje „co dělat a jak to dělat“. Může obsahovat vysoký počet parametrů, ale může být i prázdný. V takovém případě program používá vlastní výchozí hodnoty parametrů. Tento soubor pomáhá optimalizovat výpočty, ale už i výchozí parametry jsou pro běžné výpočty dostačující. Příklad souboru INCAR pro výpočet uhlíku v diamantové struktuře: C v diamantove strukture PREC = Accurate ENCUT = 520 EDIFF = 1E-05
První řádek je komentář a může obsahovat libovolný text. Další řádky obsahují informace o tom, které výchozí hodnoty se mají změnit a jak. 8
Nejčastěji se upravují hodnoty parametrů PREC – ovlivňuje přesnost s jakou se má počítat, ENCUT – mezní energie (přepisuje hodnotu ENMAX v souboru POTCAR) a EDIFF – konvergenční kritérium: když nastane situace, že je změna energie během jednoho kroku menší jak hodnota EDIFF, pak se výpočty ukončí. Soubor POTCAR obsahuje pseudopotenciál pro každý druh prvku v bázi. Pokud obsahuje báze více různých prvků (například NaCl obsahuje dva různé prvky v bázi – Na a Cl), tak je nutné sloučit soubory POTCAR těchto prvků do jednoho. Dále soubor obsahuje informace o atomech daných prvků (jejich hmotnost, valenci, pseudopotenciály, výchozí energetickou mez (ENMAX a ENMIN řádek) atd.), proto není nutné specifikovat valenci a hmotnost v INCAR souboru. Soubory POTCAR obsahující pseudopotenciály všech prvků periodické soustavy jsou dodané výrobcem programu VASP spolu s kódem. Soubor POSCAR obsahuje informace o primitivních translačních vektorech a o pozicích iontů. V tomto souboru je popsána primitivní buňka, jejímž opakováním vznikne model nekonečného krystalu. Příklad souboru POSCAR pro výpočet uhlíku v diamantové struktuře: C v diamantove strukture, fcc 3.567 0.00000000 0.50000000 0.50000000 0.50000000 0.00000000 0.50000000 0.50000000 0.50000000 0.00000000 2 Direct 0 0 0 0.25 0.25 0.25 První řádek je komentář. Druhý řádek udává velikost škálovací konstanty v angströmech. Jako škálovací konstanta byla v tomto příkladu použita velikost mřížkového parametru. Další tři řádky popisují translační vektory v násobcích škálovací konstanty. Na dalším řádku je určen počet atomů v bázi. Obsahuje-li báze atomy více prvků, udávají se jejich počty postupně. Pořadí musí být stejné jako v souboru POTCAR. Další řádek pomocí klíčového slova (přičemž důležité je jen první písmeno) určuje, zda dále uvedené hodnoty budou v násobcích škálovacího parametru (Direct), nebo v souřadnicích v Kartézském systému souřadnic (Cartesian). Další řádky popisují polohy bázových atomů. Soubor KPOINTS určuje rozdělení k-bodů v Brillouinově zóně. Zde je možnost nechat vytvořit síť k-bodů automaticky, nebo zadat body explicitně. Při explicitním zadání jsou opět dvě možnosti, a to buď přesně popsat souřadnice k-bodů, podobně jako souřadnice atomů v souboru POSCAR, nebo generováním řetězců k-bodů spojujících specifické body Brillouinovi zóny. Poslední způsob je vhodný pro výpočty pásových struktur.
9
3.3. Soubory používané programem VASP VASP užívá relativně velké množství vstupních a výstupních souborů. INCAR - vstupní STOPCAR - vstupní POTCAR - vstupní KPOINTS - vstupní IBZKPT - výstupní POSCAR - vstupní CONTCAR - výstupní CHGCAR - vstupní i výstupní CHG - výstupní WAVECAR - vstupní i výstupní TMPCAR - vstupní i výstupní EIGENVAL - výstupní DOSCAR - výstupní PROCAR - výstupní OSZICAR - výstupní PCDAT - výstupní XDATCAR - výstupní LOCPOT - výstupní ELFCAR - výstupní PROOUT – výstupní OUTCAR - výstupní O hlavních vstupních souborech bylo psáno výše, nyní se zaměřím na výstupní soubory, které používám. Jedná se o soubory CHGCAR, DOSCAR a OUTCAR. Soubor CHGCAR obsahuje údaje o elektronové struktuře. Hlavička souboru obsahuje informace ze souboru POSCAR. Po hlavičce souboru následuje pět sloupců s hodnotami hustoty valenčních elektronů. Soubor DOSCAR obsahuje informace o hustotě stavů. Příklad začátku souboru DOSCAR pro výpočet uhlíku v diamantové struktuře: C v diamantove strukture 15,29850416 -14,15472340 301 -14,155 0,0000E+00 0,0000E+00 -14,057 0,0000E+00 0,0000E+00 -13,958 0,0000E+00 0,0000E+00 -13,860 0,0000E+00 0,0000E+00 -13,762 0,0000E+00 0,0000E+00 ...
3,32699437
Hlavička souboru obsahuje komentář, informace v jakém rozsahu energií a pro kolik energií byla hustota stavů počítána a pak hodnotu energie odpovídající Fermiho mezi. Dále soubor pokračuje až do konce (v nejjednodušším případě jako je například tento) ve třech sloupcích, kde první 10
sloupec udává energii, druhý udává hustotu stavů a třetí integrovanou hustotu stavů. Soubor OUTCAR je hlavní výstupní soubor, který obsahuje shrnutí nejdůležitějších vstupních dat a informace o průběhu a výsledcích výpočtů. Obsahuje mimo jiné mřížkové parametry, tenzory napětí, celkovou energii a magnetický moment. Na konci souboru jsou také informace o množství paměti potřebné k výpočtům a také, jak dlouho výpočty trvaly.
3.4. Postup práce s programem VASP Prvně musím vytvořit vstupní soubory pro výpočet, jedná se o soubory KPOINTS, INCAR, POSCAR a POTCAR. Pro všechny výpočty používám mřížku o rozměrech 21x21x21 K-bodů definovanou v souboru KPOINTS. U souboru INCAR nastavuji jen hodnotu ENCUT, kterou nastavuji na hodnotu 130% největší hodnoty ENMAX ze všech prvků v bázi. Hodnoty ENMAX prvků jsou uvedeny v databázi pseudopotenciálů dodávané s programem VASP. Nakonec vytvořím soubor POTCAR pomocí databáze pseudopotenciálů. Po vytvoření vstupních souborů je spuštěn program VASP a vzniklé výstupní soubory, jako jsou OUTCAR, DOSCAR, CHGCAR, … jsou uloženy k dalšímu zpracování. Následně lze tento postup zopakovat pro změněnou hodnotu objemu. Celý postup lze samozřejmě zautomatizovat. Z uložených souborů OUTCAR jsou poté načteny složky tenzoru napětí a hodnota celkové energie, které jsou klíčové pro následné zpracovávání výsledků.
11
4. Výsledky výpočtů Při izotropním zatížení se mění jen objem, tvar (symetrie) krystalu zůstává stejná a proto pro popis deformace stačí jen relativní objem v. Relativní objem v je dán vztahem (1) kde V je aktuální objem a V0 je objem krystalu v základním stavu.
4.1. Vyhodnocování výsledků 4.1.1. Základní stav Mým úkolem je počítat vlastnosti systémů v základním stavu. Energii základního stavu krystalu E0 lze nalézt jako minimum závislosti celkové energie na objemu, popřípadě na mřížkovém parametru. Po té, co základní stav najdu, mohu znovu spustit program VASP, nyní jen pro odpovídající velikost škálovacího parametru. Takto získám například soubory CHGCAR a DOSCAR pro základní stav krystalu. 4.1.2. Modul objemové pružnosti Pro výpočet modulu objemové pružnosti B používám dva vzorce (2) (3) kde σ je mechanické napětí, v je relativní objem a E je celková energie systému. Ačkoli jsou oba vztahy v principu shodné, neboť platí , je v jejich praktickém použití rozdíl. Program VASP umožňuje počítat nejen celkovou energii E krystalu, ale i tenzor napětí . Rozdíl v hodnotách B, spočteném ze vztahů (2) a (3) musí být dostatečně malý. Výraznější rozdíl by signalizoval nedostatečnou přesnost výpočtů. 4.1.3. Hustota náboje Data obsahující informace o hustotě valenčních elekronů se nacházejí v souborech CHGCAR a CHG. Soubor CHG obsahuje v podstatě stejná data, jen zaokrouhlená, což zmenšuje velikost souboru. Pro vykreslení hustoty elektronů 12
však obvykle zcela postačí. Tato data je schopen zpracovat například program VESTA [9] 4.1.4. Hustota stavů Data obsahující informace o hustotě stavů se nacházejí v souboru DOSCAR. Tato data je možné vykreslit například programem Excel. 4.1.5. Mez izotropní pevnosti v tahu Pevnost v tahu stanovuji tak, že zvětšuji velikost objemu krystalu (při zachování jeho symetrie) a hledám maximální hodnotu mechanického napětí. Následné snižování napětí při zvětšování objemu je způsobeno zpřetrháním vazeb, tedy překročením meze pevnosti.
4.2. Hliník Hodnotu ENCUT pro výpočet hliníku nastavuji na 312 eV. Hliník krystalizuje v kubické plošně centrované mřížce (fcc) a jeho rovnovážný mřížkový parametr je 4,05 Å [10]. Minimum energie bylo nalezeno pro mřížkový parametr o velikosti 4,05 Å (obr. 1). Velikost mřížkového parametru tedy odpovídá experimentální hodnotě.
Al 16 14
E - E0 (meV)
12 10 8 6 4 2 0 3,95
4,00
4,05
4,10
4,15
mřížkový parametr (Å)
Obrázek 1: Závislost celkové energie na velikosti mřížkového parametru fcc krystalu Al.
13
Vypočtená velikost modulu objemové pružnosti ze závislosti velikosti napětí v tahu na velikosti relativního objemu podle vzorce (2) je asi 74 GPa (obr. 2). Hodnoty napětí byly získány z napěťového tenzoru, který program VASP vypisuje do souboru OUTCAR. Ze závislosti velikosti celkové energie na velikosti relativního objemu vychází podle vzorce (3) hodnota asi 77 GPa (obr. 3). Datovými body byl proložen kubický polynom, jehož koeficienty jsou uvedeny v obrázku. Experimentální hodnota je 76 GPa [10]. Velikost modulu objemové pružnosti vypočtená ze změny relativního objemu se od experimentální hodnoty liší o 2,6% a vypočtená ze změny celkové energie o 1,3% experimentální hodnoty.
Al 3 y = 74,443x - 74,36 R² = 0,9977
napětí v tahu (GPa)
2 1 0 0,96
0,97
0,98
0,99
1
1,01
1,02
1,03
1,04
-1 -2 -3
relativní objem
Obrázek 2: Závislost velikosti napětí v tahu na velikosti relativního objemu.
E - E0 (meV)
Al 9 8 7 6 5 4 3 2 1 0
y = 3867,2x2 - 7745,3x + 3878,1 R² = 0,998
0,95
0,96
0,97
0,98
0,99
1
1,01
1,02
1,03
1,04
1,05
relativní objem
Obrázek 3: Závislost celkové energie na velikosti relativního objemu.
14
1,06
Zobrazil jsem, pomocí programu VESTA, hustotu náboje valenčních elektronů v rovině (111) primitivní buňky fcc krystalu hliníku (obr. 4a). Tuto rovinu jsem pak zvýraznil ještě ve 2D (obr. 4b), aby bylo dobře vidět rozdělení valenčních elektronů. Je vidět, že hustota valenčních elektronů je sféricky symetrická okolo atomů. Největší hustota je ve vzdálenosti 0,3 násobku mřížkové konstanty od středu atomových sfér. Blíže k jádru hustota rychle klesá k nule, neboť s elektrony vnitřních slupek se nepočítá a zároveň s rostoucí vzdáleností od jader také klesá, ale pomaleji, neboť tam se stále mohou pohybovat volné elektrony.
Obrázek 4: a) Hustota valenčních elektronů v rovině (111) v primitivní buňce fcc krystalu Al. b) Daná rovina zvýrazněna ve 2D. Obr. 5 zobrazuje hustotu stavů a integrovanou hustotu stavů pro rovnovážný stav krystalu hliníku. Nulová energie zde odpovídá Fermiho mezi.
Al 1,0 0,9 0,8 hustota stavů
0,7 0,6 Hustota stavů
0,5 0,4
Integrovaná hustota stavů / 10
0,3 0,2 0,1 0,0
-8
-6
-4
-2
0
2
energie (eV)
15
4
6
Obrázek 5: Hustota stavů a integrovaná hustota stavů pro fcc krystal Al. V oblasti od -8 eV do -4 eV je hustota stavů velmi podobná hustotě stavů √
pro plyn volných elektronů [11]: ( ) , kde m je hmotnost elektronů a E je energie, pro kterou má být N(E) určena, h je Plankova konstanta. Integrovaná hustota stavů nabývá na Fermiho mezi hodnoty 3, což odpovídá počtu valenčních elektronů. Při výpočtu meze izotropní pevnosti v tahu mi maximální napětí vyšlo 11,1 GPa, což se liší o 0,4% od hodnoty 11,15 GPa, vypočtené Ogatou a kol. [12]. V obrázku 6 je zobrazena závislost izotropního napětí na velikosti relativního objemu. Jak je patrné, tak největší napětí je při relativním objemu 1,48, který menší o 1,8% než 1,507 [12].
Al mechanické napětí (GPa)
11,2
11,1
11,1 11,0 10,9 10,8 10,7 10,6 1,30
1,35
1,40
1,45
1,50
1,55
1,60
1,65
1,70
relativní objem
Obrázek 6: Závislost mechanického napětí na velikosti relativního objemu pro fcc krystal Al.
4.3. Diamant Hodnotu ENCUT pro výpočet hliníku nastavuji na 520 eV. Uhlík v diamantové struktuře krystalizuje v kubické plošně centrované mřížce (fcc) se dvěma atomy v bázi o souřadnicích (0 0 0) a (0,25 0,25 0,25) a jeho rovnovážný mřížkový parametr je 3,57 Å [10]. Minimum energie bylo nalezeno pro mřížkový parametr o velikosti 3,57 Å (obr. 8). Velikost mřížkového parametru tedy odpovídá experimentální hodnotě.
16
C 80 70
E - E0 (meV)
60 50 40 30 20 10 0 -10 3,45
3,50
3,55 3,60 mřížkový parametr (Å)
3,65
3,70
Obrázek 8: Závislost celkové energie na velikosti mřížkového parametru diamantu. Vypočtená velikost modulu objemové pružnosti ze závislosti velikosti napětí v tahu na velikosti relativního objemu podle vzorce (2) je asi 444 GPa (obr. 9). Ze závislosti velikosti celkové energie na velikosti relativního objemu vychází podle vzorce (3) hodnota asi 424 GPa (obr. 10). Experimentální hodnota je 432 GPa [10]. Velikost modulu objemové pružnosti vypočtená ze změny relativního objemu se od experimentální hodnoty liší o 2,8% a vypočtená ze změny celkové energie o 1,9% experimentální hodnoty.
C 20 y = 443,62x - 442,55 R² = 0,9983
napětí v tahu (GPa)
15 10 5 0 -5
0,96
0,97
0,98
0,99
1
1,01
1,02
1,03
1,04
-10 -15
relativní objem
Obrázek 9: Závislost velikosti napětí v tahu na velikosti relativního objemu.
17
C 25
E - E0 (meV)
20
y = 15026x2 - 30242x + 15216 R² = 0,9987
15 10 5 0 0,96 -5
0,97
0,98
0,99
1
1,01
1,02
1,03
1,04
1,05
relativní objem
Obrázek 10: Závislost celkové energie na velikosti relativního objemu. Zobrazil jsem, pomocí programu VESTA, hustotu náboje valenčních elektronů v rovině (112) primitivní buňky uhlíku v diamantové struktuře (obr. 11a). Tuto rovinu jsem pak zvýraznil ještě ve 2D (obr. 11b), aby bylo dobře vidět rozdělení valenčních elektronů. Je patrné, že největší hustota valenčních elektronů je mezi nejbližšími sousedy, kde je také nejsilnější vazba a naopak mimo spojnice nejbližší sousedů hustota velmi rychle klesá. Dále jsem vyznačil isoplochy pro stav s minimální energií (obr. 12a) a pro stav s maximálním napětím (obr. 12b). U stavu s maximálním napětím je objem krystalu 1,6 násobkem objemu krystalu s minimální energií a proto pro možnost porovnávání odpovídajících hladin zobrazuji u stavu s maximálním napětím izoplochy hladiny o 1,6 násobek menší. Z obrázků je dobře patrné snížení hustoty valenčních elektronů dané hladiny mezi nejbližšími atomy, které bylo způsobeno zpřetrháním vazeb.
Obrázek 11: a) Hustota valenčních elektronů v rovině (112) v primitivní buňce. b): Daná rovina zvýrazněna ve 2D.
18
Obrázek 12: a) Isoplochy hustoty valenčních elektronů v krystalu diamantu hladiny pro stav s minimální energií. b) Isoplochy hladiny pro stav s maximálním napětím. Kde je Bohrův poloměr. Obr. 13 zobrazuje hustotu stavů a integrovanou hustotu stavů pro rovnovážný stav krystalu diamantu. Z grafu hustoty náboje je dobře vidět velikost zakázaného pásu 3 eV, což je výrazně méně jak experimentální hodnota 5,4 eV [10]. Tato chyba je známý problém ab initio metod, kdy v rámci aproximací dochází k velkému podhodnocování velikosti zakázaného pásu. Z grafu integrované hustoty je zase dobře vidět, že interakcí se účastní 8 elektronů, neboli 4 valenční elektrony každého z obou atomů uhlíku tvořících krystalovou bázi.
C 2,0 1,8 1,6 hustota stavů
1,4 1,2
Hustota stavů
1,0 0,8
Integrovaná hustota stavů / 10
0,6 0,4 0,2 0,0
-5
-3
-1
1
3
5
energie (eV)
Obrázek 13: Hustota stavů a integrovaná hustota stavů pro krystal C v diamantové struktuře. Při výpočtu meze izotropní pevnosti v tahu mi maximální napětí vyšlo 79,9 GPa, které se liší o 9,8% od hodnoty 88,54 GPa [12]. Takto velký rozdíl je způsoben použitím jiných aproximací pro výpočet, kdy autoři článku [12] použili aproximace US-LDA (ultra měkké pseudopotenciály a aproximace lokální 19
hustoty), kdežto já používám PAW-GGA (aproximaci zobecněného gradientu). V obrázku 14 je zobrazena závislost izotropního napětí na velikosti relativního objemu. Jak je patrné, tak největší napětí je při relativním objemu 1,60, který je menší o 1,4% od hodnoty 1,623 [12].
C mechanické napětí (GPa)
80,0 79,9
79,8 79,6 79,4 79,2 79,0 78,8 1,45
1,50
1,55
1,60
1,65
1,70
relativní objem
Obrázek 14: Závislost mechanického napětí na velikosti relativního objemu krystalu diamantu.
4.4. NaCl Hodnotu ENCUT pro výpočet hliníku nastavuji na 340 eV. NaCl krystalizuje v kubické plošně centrované mřížce (fcc) se dvěma atomy v bázi o souřadnicích (0 0 0) a (0, 5 0,5 0,5). Tato struktura bývá rovněž označována jako B1. Rovnovážný mřížkový parametr je 5,64 Å [10]. Minimum energie bylo nalezeno pro mřížkový parametr o velikosti 5,65 Å (obr. 15). Velikost mřížkového parametru se liší od experimentální hodnoty o 0,2%.
20
NaCl 25
E - E0 (meV)
20 15 10 5 0 5,50
5,55
5,60
5,65
5,70
5,75
5,80
5,85
mřížkový parametr (Å)
Obrázek 15: Závislost celkové energie na velikosti mřížkového parametru. Vypočtená velikost modulu objemové pružnosti ze závislosti velikosti napětí v tahu na velikosti relativního objemu podle vzorce (2) je asi 24 GPa (obr. 16). Ze závislosti velikosti celkové energie na velikosti relativního objemu vychází podle vzorce (3) hodnota asi 25 GPa (obr. 17). Oba výsledky se velice dobře shodují s experimentální hodnotou 24 GPa [10].
NaCl 1
napětí v tahu (GPa)
0,8
y = 24,447x - 24,397 R² = 0,998
0,6 0,4 0,2 0 -0,2 0,96
0,97
0,98
0,99
1
1,01
1,02
1,03
1,04
-0,4 -0,6 -0,8
relativní objem
Obrázek 16: Závislost velikosti napětí v tahu na velikosti relativního objemu.
21
NaCl 3,5 3,0 y = 3445,9x2 - 6893,1x + 3447,2 R² = 0,9992
E - E0 (meV)
2,5 2,0 1,5 1,0 0,5 0,0 0,96
0,97
0,98
0,99
1
1,01
1,02
1,03
1,04
relativní objem
Obrázek 17: Závislost celkové energie na velikosti relativního objemu. Zobrazil jsem, pomocí programu VESTA, hustotu náboje vodivostních elektronů v rovině (110) primitivní buňky krystalu NaCl (obr. 18a). Tuto rovinu jsem pak zvýraznil ještě ve 2D (obr. 18b), aby bylo dobře vidět rozdělení náboje. To, že není vidět žádná hustota náboje v okolí atomů sodíku je způsobeno tím, že jsou sodíky a chlóry spojené iontovou vazbou, kde se valenční elektron sodíku přiblíží k chlóru a v podstatě se tváří, jako by mu patřil.
Obrázek 18: a) Hustota valenčních elektronů v rovině (110) v primitivní buňce krystalu NaCl. b) Daná rovina zvýrazněna ve 2D. Pozice atomů Na jsou znázorněny hnědě a atom Cl černě. Obr 19 zobrazuje hustotu stavů a integrovanou hustotu stavů pro rovnovážný stav krystalu NaCl. Z grafu hustoty náboje je dobře vidět velikost zakázaného pásu 6 eV, což je opět výrazně méně než experimentální hodnota 10 eV [13]. Z grafu integrované hustoty je zase dobře vidět, že interakcí se účastní 8 elektronů, neboli 7 valenčních elektronů chlóru a 1 valenční elektron 22
sodíku. Pík na křivce hustoty stavů pro energii cca -12,5 eV odpovídá 3s slupce atomu chlóru.
NaCl 16 14 hustota stavů
12 10 Hustota stavů
8 6
Integrovaná hustota stavů
4 2 0 -15
-10
-5
0
5
10
15
energie (eV)
Obrázek 19: Hustota stavů a integrovaná hustota stavů pro krystal NaCl. Při výpočtu meze izotropní pevnosti v tahu mi maximální napětí vyšlo 3,95 GPa. Tato hodnota je o 21,9% nižší než hodnota 5,06 GPa publikována Ogatou a kol. [12], vypočítaná s použitím aproximace US-LDA, kdežto já používám PAW-GGA. Můj výsledek je nižší o 8,1% oproti hodnotě 4,3 GPa, která byla spočtena pomocí Boruova-Mayerova empirického potenciálu [14]. Opět je rozdíl způsoben použitím jiných aproximací. V obrázku 20 je zobrazena závislost izotropního napětí na velikosti relativního objemu. Jak je patrné, tak největší napětí je při relativním objemu 1,56, který je větší o 2,2% než dříve publikovaná hodnota 1,527 [12].
mechanické napětí (GPa)
NaCl 3,96
3,95
3,95 3,94 3,93 3,92 3,91 1,40
1,45
1,50
1,55
1,60
1,65
1,70
relativní objem
Obrázek 20: Závislost mechanického napětí na velikosti relativního prodloužení pro krystal NaCl. 23
5. Závěr V rámci této bakalářské práce jsem se seznámil s programem VASP a pomocí něho jsem určil rovnovážné mřížkové parametry, objemové moduly pružnosti, hustoty náboje v daných rovinách a meze pevnosti při izotropním tahu u tří krystalů pevných látek. Jednalo se o krystal hliníku, který je typickým zástupcem krystalu s kovovou vazbou, o krystal diamantu, který je typickým zástupcem krystalu s kovalentní vazbou a o krystal soli kamenné, který je typickým zástupcem krystalu s iontovou vazbou. Hodnoty mřížkových parametrů a modulů objemové pružnosti se velmi dobře shodovaly s eperimentálními hodnotami uvedenými v použité literatuře. Po zobrazení hustot náboje byly dobře viditelné typické vlastnosti rozdílných typů vazeb. Z grafů hustoty stavů jsem určil velikost zakázaných pásů, které byly výrazně menší jak hodnoty v použité literatuře. Zároveň je vidět, že u kovové vazby zakázaný pás není. Nakonec jsem určil mez izotropní pevnosti v tahu a odpovídající mezní deformaci krystalů, které se většinou shodují v rámci 10% s hodnotami v použité literatuře, což bylo způsobeno použitím jiných aproximací při výpočtech.
24
Bibliografie [1] M. Friák, T. Hickel, F. Körmann, A. Udyansky, A. Dick, J. von Pezold, D. Ma, O. Kim, W.A. Counts, M. Šob, T. Gebhardt, D. Music, J. Schneider, D. Raabe, J. Neugebauer, Determining the Elasticity of Materials Employing Quantummechanical Approaches, 2010. [2] Henry Fritz Schaefer III, A history of ab initio computational quantum chemistry: 1950–1960, 1988. [3] Parr, Robert G.; Weitao, Yang, Density-Functional Theory of Atoms and Molecules, New York: Oxford University Press Inc, 1989. [4] Z. Havlas, Metody a aplikace teoretické organické chemie, Ústav organické chemie a biochemie Akademie věd České republiky, 1997. [5] John P. Perdew, Kieron Burke, Matthias Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77, 3865, 1996. [6] Georg Kresse, Martijn Marsman, and Jürgen Furthmüller, VASP the GUIDE, Wien: Computational Materials Physics, Faculty of Physics, Universität Wien, 2012. [7] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169, 1996. [8] G. Kresse, D. Joubert, Phys. Rev. B 59, 1758, 1999. [9] K. Momma, F. Izumi, VESTA: a Three-Dimensional Visualization System for Electronic and Structural Analysis, 2013. [10] C. Kittel, Introduction to solid state physics, New York: Wiley, 2005. [11] David Halliday, Robert Resnick, Jearl Walker, Fyzika, Prometheus, 2001. [12] Shigenobu Ogata, Ju Li, Naoto Hirosaki, Yoji Shibutani, Sidney Yip, Ideal shear strain of metals and ceramics, 2004. [13] Kratochvíl Bohumil, Švorčík Václav, Vojtěch Dalibor, Úvod do studia materiálů, Praha: Vysoká škola chemicko-technologická v Praze, 2005. [14] N. H. MacMillan, The theoretical strength of solids, Journal of Materials Science, 1972.
25
