Michal Hermanowicz's web site | Last update: 2017-11-09 | A Lynx-friendly site!

[Back to index]
[Back to Materiały dydaktyczne / Teaching resources]

Symulacje komputerowe z pierwszych zasad (SKIZ) – FT / stopień I / sem. 5

  1. Informacje ogólne
    1. Organizacja
    2. Literatura – dostępna w bibliotece PP – i zasoby elektroniczne
    3. Zaliczenie
    4. Sprawozdanie
  2. Materiały do zajęć
    1. Informacje wstępne
    2. Zapoznanie się ze środowiskiem pracy
    3. Implementacja teorii funkcjonału gęstości – program SIESTA
    4. Relaksacja strukturalna
    5. Zbieżność
    6. Ciała stałe i teoria pasmowa
    7. Projekt

I. Informacje ogólne

0. Organizacja

Wszelkie informacje organizacyjne zostaną podane na pierwszych zajęciach.

1. Literatura – dostępna w bibliotece PP – i zasoby elektroniczne

2. Zaliczenie

Ocena końcowa będzie w pierwszym przybliżeniu średnią poniższych zmodyfikowaną o aktywność (lub jej brak) na zajęciach (obie oceny musza być pozytywne): * Termin i zakres materiału zostaną ustalone na zajęciach.

3. Sprawozdanie

Sprawozdanie powinno zostać wykonane w formacie PDF za pomocą dowolnego programu (sugerowane jest jednak wolne oprogramowanie, np. LaTeX, Libre Office lub inne). Plik proszę wysłać za pomocą poczty elektronicznej na adres podany na zajęciach (w ustalonym terminie). Raport powinien zawierać:
  1. wstęp – co jest liczone i w jakim celu. Proszę scharakteryzować badany układ – czym jest, jak jest zbudowany i jakie ma właściwości fizyczne/chemiczne. Warto jest cytować źródła literaturowe (Scopus, Google Scholar) i odszukać, co wiadomo na temat badanego układu;
  2. metoda – jaką metodą wykonano obliczenia. Proszę pamiętać, że oprogramowanie, którego używamy to tylko jedna z wielu implementacji metody. Innymi słowy – fizyka jest ważniejsza niż nazwy programów. Słowa kluczowe, które charakteryzują metodę: teoria funkcjonału gęstości (ang. density functional theory, DFT), przybliżenie lokalnej gęstości (ang. local density approximation, LDA), baza fal płaskich, pseudopotencjał;
  3. wyniki – co uzyskano w wyniku obliczeń. Proszę zawrzeć w raporcie użyte parametry, scharakteryzować użytą komórkę elementarną (ile atomów zawiera, jakie pierwiastki, ile elektronów walencyjnych zawiera hamiltonian, etc.). W sprawozdaniu proszę umieścić uzyskane wykresy wraz z opisem tego, co przedstawiają i (jeżeli to możliwe) porównać uzyskane wyniki z danymi literaturowymi;
  4. podsumowanie – wnioski z wykonanej pracy.

[Powrót na górę strony]

II. Materiały do zajęć

0. Informacje wstępne

Zajęcia odbywają się w pracowni komputerowej (sala 602) – regulamin pracowni/zasady BHP zostaną omówione na pierwszym spotkaniu i są wywieszone w sali. Pracujemy w ładowanym domyślnie systemie operacyjnym (Debian GNU/Linux) – do zajęć przydatna będzie instrukcja jego użytkowania dostępna TUTAJ (warto mieć ją otwartą w trakcie zajęć).

1. Zapoznanie się ze środowiskiem pracy

Według instrukcji na zajęciach. Zagadnienia: terminal systemu operacyjnego GNU/Linux, system plików, elementy programowania w językach interpretowanych (Bash), edytor tekstu i inne.

2. Implementacja teorii funkcjonału gęstości – program SIESTA

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 out
lub
$ 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:

$ xv2xsf
Jako 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ć?

[Powrót na górę strony]

3. Relaksacja strukturalna

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
done
lub
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 relaksacji
Zadeklarowaliś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ą:
$ xv2xsf
Jako 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?

[Powrót na górę strony]

4. Zbieżność

Według instrukcji podanych na zajęciach: przeanalizujemy zagadnienie zbieżności i przeprowadzimy jej testy celem doboru parametrów obliczeń.

[Powrót na górę strony]

5. Ciała stałe i teoria pasmowa

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).

[Powrót na górę strony]

6. Projekt

Według instrukcji podanych na zajęciach: indywidualne projekty polegające na zbadaniu zadanego układu (w oparciu o omówione wcześniej zagadnienia).

Copyright (c) 2017 M. Hermanowicz