[Back to index]
[Back to Teaching materials]
UWAGA: zgodnie z zapowiedzią – w dniach 23 i 30 listopada zajęcia nie odbędą się; spotykamy się 6-go grudnia (09:45).
Zapoznamy się z jedną z implementacji teorii funkcjonału gęstości (ang. density functional theory, DFT) – program SIESTA. W szczególności: jak działają programy obliczeniowe, co można zrobić z pomocą DFT i jakie sa jej słabe i mocne strony. W ramach prostych ćwiczeń wykonamy obliczenia dla układów o zredukowanej wymiarowości: pojedynczy atom/cząsteczka.
A.) Za pomocą dowolnego edytora tekstu utworzyć plik wsadowy do obliczeń cząsteczki wody wykorzystując poniższy szablon (h2o.fdf):
SystemName Water molecule # Nazwa własna SystemLabel h2o # Symbol (bez spacji, np. wzór sumaryczny) NumberOfAtoms 3 # Liczba atomów w układzie NumberOfSpecies 2 # Liczba pierwiastków %block ChemicalSpeciesLabel 1 8 O # L. porządkowa, l. atomowa, symbol pierwiastka 2 1 H # L. porządkowa, l. atomowa, symbol pierwiastka %endblock ChemicalSpeciesLabel MaxSCFIterations 100 # Maksymalna liczba kroków cyklu samozgodnego LatticeConstant 10 Ang # Stała sieci krystalicznej [Ang] %block LatticeVectors # Wektory sieci 1.000 0.000 0.000 # X Y Z (skalowane przez stałą sieci) 0.000 1.000 0.000 0.000 0.000 1.000 %endblock LatticeVectors AtomicCoordinatesFormat Ang # Jednostka, w której wyrazone są współrzędne %block AtomicCoordinatesAndAtomicSpecies # Wspórzędne atomów w układzie 0.000 0.000 0.000 1 # X, Y, Z, liczba porządkowa pierwiastka 0.757 0.586 0.000 2 -0.757 0.586 0.000 2 %endblock AtomicCoordinatesAndAtomicSpecies
i zapisać w bieżącym katalogu roboczym przeznaczonym na wykonanie obliczeń.
B.) Do przeprowadzenia obliczeń niezbędny jest jeszcze jeden element: pseudopotencjał – należy pobrać go dla każdego z pierwiastków (wodór, tlen) według instrukcji podanej na zajęciach i zapisać w tym samym katalogu, w którym znajduje się główny plik wsadowy.
C.) Uruchomić obliczenia zapisując standardowe wyjście w pliku o nazwie out:
$ siesta < h2o.fdf | tee outlub
$ siesta < h2o.fdf > out
D.) Rozpisać konfiguracje elektronowe wszystkich atomów wchodzących w skład liczonego układu. Ile jest sumarycznie elektronów, a ile z nich uwzględniono w obliczeniach?
E.) Ile iteracji obliczeń zostało wykonanych?
F.) Ile wynosi energia całkowita wynikająca z rozwiązania równania Kohna-Shama [eV]?
G.) Zbadać wartości własne w pliku *.EIG (od eigenvalues). Ile ich jest? Czym są i co oznaczają?
H.) Obejrzeć strukturę liczonego układu wykonując poniższe polecenia:
$ xv2xsfJako label podajemy symbol (prefix) układu, tzn. h2o, a następnie:
$ xcrysden --xsf h2o.XSF
I.) Sformułować wnioski: czy otrzymane wyniki są poprawne? Jak to sprawdzić?
Przeprowadzone do tej pory obliczenia pozwalały na wyznaczenie energii całkowitej oraz właściwości elektronowych zadanego układu. Taki rachunek wymaga jednak dokładnej znajomości struktury liczonego materiału, tzn. współrzędnych tworzących go atomów, a w przypadku kryształów także stałych sieci krystalicznej. Czasami te informacje nie są znane lub są znane jedynie w przybliżeniu – można wówczas wykonać relaksację strukturalną, czyli wyznaczyć położenia atomów w układzie na podstawie obliczeń energii całkowitej: optymalna struktura = najniższa energia całkowita. W przypadku skomplikowanych (wieloatomowych) układów to zagadnienie jest nietrywialne, ale w poniższym ćwiczeniu wykonamy relaksację prostej (jednowymiarowej) cząsteczki wodoru (H2) demonstrując koncepcyjnie zasadę działania.
Do przeprowadzenia obliczeń możemy wykorzystać bezpośrednio programu SIESTA, jednak dla lepszego zilustrowania problemu zaczniemy od mniej efektywnej, ale bardziej obrazowej metody. Wykonamy szereg obliczeń energii całkowitej (E) cząsteczki wodoru różniących się jedynie długością wiązania (d), a następnie wykreślimy wartość E w funkcji d. Minimum otrzymanej funkcji wyznaczy optymalną długość wiązania (w którym energia układu jest najniższa).
Manualne wykonanie n niezależnych obliczeń dla różnych długości wiązania jest możliwe, choć nieefektywne. Żeby zautomatyzować ten proces skorzystamy z funkcjonalności powłoki Bash – pętli for o ogólnej konstrukcji:
for i in 1 2 3 4 5 do echo $i donelub
for i in `seq 1 1 5`
do
echo $i
done
Poniżej kompletny skrypt, którego analizę przeprowadzimy na zajęciach.
#!/bin/bash export LC_NUMERIC=C if [ -e energie.txt ] then rm energie.txt fi for i in `seq 0.4 0.1 2.0` do cat > h2.fdf << EOF SystemName Hydrogen molecule SystemLabel h2 NumberOfAtoms 2 NumberOfSpecies 1 %block ChemicalSpeciesLabel 1 1 H %endblock ChemicalSpeciesLabel MaxSCFIterations 100 LatticeConstant 10 Ang %block LatticeVectors 1.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 1.000 %endblock LatticeVectors AtomicCoordinatesFormat Ang %block AtomicCoordinatesAndAtomicSpecies 0.000 0.000 0.000 1 $i 0.000 0.000 1 %endblock AtomicCoordinatesAndAtomicSpecies EOF echo -n $i >> energie.txt siesta < h2.fdf | grep "E_KS(eV) =" | cut -d'=' -f2 >> energie.txt done
Po upewnieniu się, że posiadamy w bieżącym katalogu pseudopotencjał dla wodoru, możemy powyższy plik uruchomić jako skrypt nadając mu wcześniej prawo wykonywalności (patrz także: instrukcja):
$ chmod +x skrypt.sh $ ./skrypt.sh
Kiedy powłoka zwróci znak zachęty, możemy przeanalizować wyniki zawarte w pliku energie.txt, który powinien teraz zawierać dwie kolumny danych (x, y): długość wiązania i wartość energii całkowitej wyrażonej w elektronowoltach. Do wykreślenia wyznaczonej funkcji posłużymy się programem gnuplot:
$ gnuplot > plot "energie.txt" w lp lw 2 pt 5 ps 1.5
Z uzyskanego wykresu należy odczytać optymalną długość wiązania, dla której energia całkowita układu jest najniższa.
Nie jest to jedyna metoda wyznaczenia geometrii układu. SIESTA pozwala na przeprowadzenie automatycznej relaksacji z wykorzystaniem zaimplementowanych algorytmów optymalizacyjnych. W tym celu nie musimy już uruchamiać szeregu obliczeń, a tylko pojedynczą instancję programu. Możemy wykorzystać szablon dla cząsteczki wody (patrz punkt 2) modyfikując odpowiednio nazwę, liczbę atomów (i pierwiastków), ich współrzędne i inne. Następnie uzupełniamy nasz plik wsadowy (z rozszerzeniem .fdf) o następujący fragment:
MD.TypeOfRun CG # Rodzaj relaksacji (CG = Conjugate Gradients) MD.MaxForceTol 0.04 eV/Ang # Tolerancja na siły MD.NumCGsteps 1000 # Maksymalna liczba kroków relaksacjiZadeklarowaliśmy niniejszym użycie metody sprzężonych gradientów (CG) z kryterium zbieżności na poziomie 0.04 eV/Ang, a także maksymalną liczbę kroków relaksacji strukturalnej (1000). Po wykonaniu obliczeń możemy sprawdzić wyznaczoną strukturę odszukując w głównym pliku wynikowym (out) nowe współrzędne atomów lub rysując model dotychczasową metodą:
$ xv2xsfJako label podajemy symbol (prefix) układu, tzn. h2, a następnie:
$ xcrysden --xsf h2.XSF
Po uruchomieniu programu korzystamy z opcji Distance, żeby wyznaczyć odległość między wybranymi atomami. Porównujemy uzyskany wynik z długością wiązania odczytaną z wykresu E(d). Skąd się biorą ewentualne różnice i jak można poprawić jakość wyników?
Według instrukcji podanych na zajęciach: przeanalizujemy zagadnienie zbieżności i przeprowadzimy jej testy celem doboru parametrów obliczeń.
Struktura elektronowa jednowymiarowego łańcucha atomów wodoru – szablon pliku wsadowego (crystal.fdf):
SystemName 1D H chain # Nazwa własna SystemLabel crystal # Symbol (bez spacji, np. wzór sumaryczny) NumberOfAtoms 1 # Liczba atomów w układzie NumberOfSpecies 1 # Liczba pierwiastków SpinPolarized .false. %block ChemicalSpeciesLabel 1 1 H # L. porządkowa, l. atomowa, symbol pierwiastka %endblock ChemicalSpeciesLabel MaxSCFIterations 100 # Maksymalna liczba kroków cyklu samozgodnego LatticeConstant 10 Ang # Stała sieci krystalicznej [Ang] %block LatticeVectors # Wektory sieci 0.500 0.000 0.000 # X Y Z (skalowane przez stałą sieci) 0.000 0.500 0.000 0.000 0.000 0.500 %endblock LatticeVectors %block ProjectedDensityOfStates -20.00 20.00 0.200 1000 eV %endblock ProjectedDensityOfStates %block BandLines 1 0.000 0.000 0.000 50 0.500 0.000 0.000 %endblock BandLines BandLinesScale ReciprocalLatticeVectors AtomicCoordinatesFormat Ang # Jednostka, w której wyrazone są współrzędne %block AtomicCoordinatesAndAtomicSpecies # Wspórzędne atomów w układzie 0.0 0.0 0.0 1 # X, Y, Z, liczba porządkowa pierwiastka %endblock AtomicCoordinatesAndAtomicSpecies
Niezbędny będzie pseudopotencjał dla wodoru!
Prosze o przypomnienie wiadomości: sieć krystaliczna (rzeczywista i odwrotna), twierdzenie Blocha, teoria pasmowa, struktura elektronowa, masa efektywna. W tej części zajmiemy się układami periodycznymi (również tymi o zredukowanej wymiarowości).
Celem ćwiczenia jest zapoznanie się ze strukturą atomową i elektronową grafenu. Geometria komórki elementarnej jest określona wartościami parametrów [LatticeConstant], [LatticeVectors] oraz [AtomicCoordinatesAndAtomicSpecies]. Poniżej znajduje się zawartość bazowego pliku wsadowego, który zawiera podstawowy model układu w postaci komórki elementarnej 1x1.
SystemName Grafen SystemLabel grafen NumberOfAtoms 2 NumberOfSpecies 1 SpinPolarized .true. %block ChemicalSpeciesLabel 1 6 C %endblock ChemicalSpeciesLabel MaxSCFIterations 1000 LatticeConstant 1 Ang %block kgrid_Monkhorst_Pack 3 0 0 0 3 0 0 0 1 %endblock kgrid_Monkhorst_Pack %block LatticeVectors 2.463999822 0.000000000 0.000000000 1.231999911 2.133886441 0.000000000 0.000000000 0.000000000 15.00000000 %endblock LatticeVectors %block ProjectedDensityOfStates -20.00 20.00 0.200 1000 eV %endblock ProjectedDensityOfStates %block BandLines 1 0.00000 0.00000 0.00000 50 0.66667 0.33333 0.00000 30 0.50000 0.50000 0.00000 %endblock BandLines BandLinesScale ReciprocalLatticeVectors AtomicCoordinatesFormat Ang %block AtomicCoordinatesAndAtomicSpecies 0.000000000 0.0000000000 7.5000000000 1 1.232012231 0.7112883674 7.5000000000 1 %endblock AtomicCoordinatesAndAtomicSpecies
Pseudopotencjał dla węgla dostępny jest tutaj.
Po wykonaniu obliczeń można zwizualizować układ dotychczasową metodą – wykonując uprzednio skrypt xv2xsf (podając grafen jako prefix), a następnie:
xcrysden --xsf grafen.XSF
Analiza struktury pasmowej będzie wymagała przekonwertowania pliku *.bands do formatu kompatybilnego z programem gnuplot:
gnubands < grafen.bands > pasma
Niezbędne będzie odszukanie wartości energii Fermiego w głównym pliku wyjściowym (można skorzystać z grep). Wykres funkcji całkowitej gęstości stanów elektronowych można wykonać korzystając bezpośrednio z pliku grafen.DOS (proszę zwrócić uwagę na jego format – obliczenia są spinowo spolaryzowane!).
Plik wsadowy zawiera nowy blok danych – kgrid_Monkhorst_Pack – pozwalający zdefiniować sieć punktów k w przestrzeni odwrotnej. W ramach ćwiczenia proszę sprawdzić zbieżność obliczeń (energia, struktura pasmowa, a w szczególności DOS) względm liczby punktów k w kierunkach x i y.
Jaki powinien być kształt DOS w bliskim otoczeniu poziomu Fermiego? Uzyskany (parcjalny/lokalny) DOS można traktować jako pierwsze przybliżenie krzywej dI/dV spektroskopii tunelowej (ang. scanning tunneling spectroscopy – STS). Proszę porównać wynik z danymi odnalezionymi w dostępnej literaturze i zacytować w sprawozdaniu odpowiednie źródło.
Copyright (c) 2017, 2018 M. Hermanowicz