7. Narzędzia aplikacyjne metod DFT_________
7. Narzędzia aplikacyjne metod DFT
Przedstawione w poprzednich rozdziałach podstawy metod DFT wskazują iż
formalizm ten w wydajny i stosunkowo dokładny sposób pozwala wyznaczyć strukturę
elektronową badanych układów katalitycznych. Jednak uzbieżnienie struktury elektronowej
realizowane poprzez schemat Kohna Shama dla stworzonego przez nas modelu samo
w sobie nie niesie jeszcze pełnej informacji przydatnej dla chemika. W rozdziale tym
prezentowany jest przeglÄ…d metod opartych na DFT wykorzystywanych w ramach pracy do
interpretacji wyników doświadczalnych.
7.1 Analiza Powierzchni Energii Potencjalnej
Powierzchnia energii potencjalnej (PES, ang. Potential Energy Surface) zdefiniowana
dla metod chemii kwantowej jest funkcją położeń jąder atomowych układu w przybliżeniu
adiabatycznym. Dokładna znajomość tej powierzchni pozwala wyznaczyć stabilne geometrie
badanych układów (substraty i produkty reakcji chemicznych) oraz geometrie stanów
przejściowych stowarzyszonych z badaną reakcją w oparciu o analizę hesjanu14 czyli
macierzy drugich pochodnych energii potencjalnej. Wartości własne definiują z jakim
punktem na PES mamy do czynienia15 (Rys. 7.1):
" stabilne geometrie rozpoznajemy w punktach PES, dla których wszystkie wartości
własne są dodatnie
" stany przejściowe zdefiniowane jako punkty siodłowe pierwszego rzędu mają jedną
wartość ujemną (oznacza to urojoną wartość jednego z drgań normalnych
odpowiedzialnego za przejście przez stan przejściowy)
" mniej interesujące dla chemika maksima (wszystkie wartości na diagonali ujemne)
i siodła wyższych rzędów (rząd siodła = liczba wartości ujemnych)
"V
16
W tak określonych punktach krytycznych gradient potencjału "V = jest równy zeru,
"Xi
a zatem siła ( "V ) działająca na atomy także jest równa zero.
14
formalna nazwa to macierz Hessego
15
należy pamiętać, że sześć wartości własnych (odpowiadających rotacjom i translacjom układu) jest równych 0
16
Xi oznaczają współrzędne kartezjańskie, w których opisujemy konfigurację N jąder układu (i = 1,2,....3N)
72
__________Metodyka Badawcza
Rys. 7.1: Schemat wielowymiarowej powierzchni energii potencjalnej.
Oczywistym jest fakt, iż ze względu na znaczny wymiar przestrzeni w której rozpięta
jest funkcja PES dla większych układów, wyznaczenie jej kształtu dla całego zakresu zmian
położeń jąder atomowych jest praktycznie niemożliwe. Dlatego w obliczeniach stosuje się
odpowiednie metody poszukiwania stabilnych geometrii substratów i produktów oraz
aplikacje mające za zadanie odnalezienie na powierzchni PES punktów związanych ze
stanami przejściowymi badanych reakcji.
Optymalizacja geometrii
Zadaniem algorytmów optymalizacyjnych jest znalezienie minimum PES po starcie
z punktu poczÄ…tkowego zdefiniowanego arbitralnie (geometria startowa). IstniejÄ… trzy
podstawowe klasy wykorzystywanych algorytmów optymalizacyjnych z których każda
charakteryzuje siÄ™ wadami i zaletami.
Metoda największego spadku (ang. steepest descents, SD) wykorzystuje iteracyjny algorytm
poszukiwania liniowego przemieszczajÄ…c siÄ™ w kierunku przeciwnym do wskazywanego
przez wektor gradientu. Gdy w rezultacie cięcia powierzchni w tak zdefiniowanym kierunku
wyznaczone zostaje pierwsze przybliżenie minimum, procedura startuje z niego w kierunku
prostopadłym do poprzedniego. Zaletą metody jest prostota oraz niewielkie wymagania co do
pamięci komputerowej  przechowywany jest bowiem tylko wektor gradientu. Wady metody
to przed wszystkim wolnozbieżność w pobliżu minimum17 oraz tendencja do wpadania
17
Każdy kolejny krok korelowany jest z poprzednim w prymitywny sposób (prostopadłość geometryczna) co
prowadzi do utraty części informacji o kształcie funkcji z poprzedniego kroku
73
7. Narzędzia aplikacyjne metod DFT_________
w oscylacje w przypadku przeszukiwania podłużnych dolin potencjału [67]. W praktyce
algorytm stosowany jest wyłącznie do wstępnej optymalizacji układów dalekich od minimum.
Metody sprzężonych gradientów (ang. Conjugate Gradient, CG) w pierwszym kroku
optymalizacji wykorzystujÄ… algorytm SD ale kolejne kierunki przeszukiwania (wektory d)
budowane sÄ… w oparciu o poprzedni kierunek oraz wyznaczony w nowym punkcie wektor
gradientu (g), zgodnie z wzorem:
di = -gi + ²idi-1 (7.1)
Implemenatcje tej metody różniÄ… siÄ™ wÅ‚aÅ›nie doborem staÅ‚ej ²i a najszerzej stosowane to
metoda Polaka-Ribiere'a [68] i Fletchera-Reevesa [69] dla których stałe te mają poniżej
podaną postać:
t
gtgi g (g - g )
FR PR
i i i i -1
²i = ² =
(7.2) i (7.3)
i
t
gt gi-1 g g
i-1 i -1 i -1
Metody te sprawdzają się najlepiej w przypadku powierzchni łatwych do przybliżenia
funkcjami kwadratowymi. W innych obszarach często w toku optymalizacji muszą być
restartowane (wyzerowanie parametru ²i). Mimo to wykazujÄ… one znacznie lepszÄ… zbieżność
od metod SD, nie zwiększając znacznie wymagań pamięciowych (program magazynuje dwa
wektory: gradientu wykonywanego kroku i kierunek poprzedniego etapu).
Metody Newtona-Raphsona zapewniają dokładniejsze wyznaczenie kierunku poszukiwania
minimum poprzez uwzględnieniu informacji różniczkowej drugiego rzędu. Model otoczenia
bieżącego punktu ma wtedy postać kwadratową:
1
f (x) H" f (x0) + gt (x - x0) + (x - x0)t H(x - x0)
(7.4)
2
gdzie H jest symetryczną i dodatnio określoną macierzą Hessego aproksymowanej funkcji zaś
g wektorem gradientu w bieżącym punkcie. Z przyrównania pochodnej modelu względem (x -
xo) do 0 otrzymuje siÄ™ kierunek d poszukiwania minimum (kierunek Newtona):
d = -H-1g (7.5)
Zaletą metod jest możliwość poszukiwania innych niż minima stanów stacjonarnych
(algorytm szuka najbliższego stanu stacjonarnego) oraz dobra zbieżność w pobliżu punktu
stacjonarnego. Wadą metody są przede wszystkim problemy z dokładnym wyznaczeniem
74
__________Metodyka Badawcza
Hesjanu dlatego opracowano metody z iteracyjnÄ… aktualizacjÄ… tej macierzy18 nazywane
metodami quasi-Newtonowskimi19 (lub metodami zmiennej metryki). Najpopularniejsze
formuły aktualizacji to BFGS20 [70,71]:
"g"gt H"x"xtH
"H = -
(7.6)
BFGS
"gt"x "xtH"
I DFP21 [72] aktualizująca bezpośrednio odwrotność hesjanu. Także przechowywanie
i diagonalizacja H staje się problemem dla większych układów. W związku faktem, iż
większość pozadiagonalnych wyrażeń jest bliska zeru hesjan przybliża się obciętą formą
zachowując tylko podmacierze o małych wymiarach (np. 3x3) zbudowane na elementach
diagonalnych. Diagonalizacja takiej serii macierzy jest znacznie tańsza niż hesjanu o dużym
wymiarze.
Poszukiwanie stanów przejściowych
Poprawne określenie geometrii stanów przejściowych charakteryzujących reakcję
powierzchniowÄ… jest warunkiem koniecznym dla zrozumienia zachodzÄ…cego procesu
katalitycznego. Algorytmy matematyczne pozwalające lokalizować stany przejściowe na PES
są zazwyczaj bardzo złożone obliczeniowo, a ich stosowalność niestety nie jest uniwersalna.
Możemy wyróżnić tutaj dwie podstawowe grupy  metody interpolacyjne i metody lokalne.
Metody interpolacyjne są najczęściej stosowane w chemii kwantowej i bazują na interpolacji
pomiędzy strukturami reagentów i produktów. Ich główną wadą jest fakt, iż nie poszukują
wprost punktu w którym gradienty się zerują lecz opierają sie tylko na porównaniu energii
kolejnych obrazów układu generowanych przez interpolację. Zaletą jest możliwość ich użycia
z wykorzystaniem startowych struktur dalekich od TS oraz tani koszt obliczeń (obliczamy
tylko gradient). Przykładem takich metod są metody NEB oraz LST/QST.
NEB [73,74]  (ang. Nudged Elastic Band) jest wieloobrazową metodą określania
ścieżki o minimalnej energii pomiędzy określonymi reagentami i produktami. Polega ona na
optymalizacji serii obrazów wygenerowanych na podstawie interpolacji pomiędzy strukturą
produktów i reagentów. Geometria obrazów (punktów PES) jest optymalizowana z więzami,
18
w ogólnym przypadku startowy hesjan jest początkowo macierzą jednostkową, a uzupełnia się go w każdym
kroku o drugie pochodne energii obliczanych w kierunkach poszukiwania - pozostałe kierunki pozostają
arbitralne
19
metody z bezpośrednim wyznaczaniem hesjanu są nazywane metodami Newtona
20
od nazwisk twórców: Broyden, Fletcher, Goldfarb, Shanno
21
od nazwisk twórców: Davidion, Fletcher, Powell
75
7. Narzędzia aplikacyjne metod DFT_________
które powodują iż optymalizacja prowadzona jest względem: (A) sił działających na atomy
prostopadle do kierunku współrzędnej reakcji, (B) narzuconych sił utrzymujących geometrię
struktur zgodną z ścieżkę energii minimalnej (MEP, Minimal Energy Pathway) [67]. Po
optymalizacji22 struktura o najwyższej energii jest dobrym przybliżeniem stanu
przejściowego. Metoda CI NEB [75] (amg. Climbing, wspinanie) jest wyjściem poza
klasyczną metodę poprzez dodanie dla obrazu o najwyższej energii etapu, w którym układ
maksymalizuje swoją energię dzięki odwróceniu kierunku działającego nań potencjału
(wspinanie do właściwego siodła). Kolejnym rozwinięciem metody jest podejście
zaproponowane przez Maragakisa poprawiające wydajność i dokładność metody w
odnajdywaniu geometrii TS [76].
LST/QST (ang. Linear/Quadratic Synchronous Transit): to grupa metod po raz
pierwszy wprowadzonych przez Halgrena i Lipscomba [77]. Startowym etapem tej metody
jest znalezienie pierwszego przybliżenia stanu przejściowego (TS1) realizowane dzięki
procedurze (LST). Interpolacja liniowa pomiędzy reagentami i substratami daje serię
geometrii, z których ta związana z najwyższą energią definiowana jest jako TS1. Z tego
punktu na PES odbywa się minimalizacja energii w kierunku sprzężonym ze współrzędną
reakcji (CG). Następnie dzięki interpolacji kwadratowej (QST)23 poszukiwany jest drugie
przybliżenie geometrii stanu przejściowego (TS2) i wykonywana kolejna minimalizacja
metodą CG. Kroki QST i CG są powtarzane iteracyjnie do osiągnięcia zbieżności [78].
Metody lokalne wymagają wstępnego zlokalizowania TS. Klasycznym przykładem
tego typu metod jest metoda Newtona-Raphsona wykorzystujÄ…ca wprost Hesjan do
odnalezienia geometrii stanu przejściowego poprzez śledzenie jego ujemnej składowej (ang.
eigenvector following). Metoda maksymalizuje energię układu wzdłuż tego kierunku
minimalizując ją w kierunkach pozostałych. Oczywistą wada jest konieczność wyznaczenia
hesjanu i koszty zwiÄ…zane z jego przechowywaniem i operacjami na nim wykonywanymi.
Analiza wibracyjna
Hiperpowierzchnia energii potencjalnej ma dla większych molekuł bardzo złożoną
postać. Funkcje tą charakteryzuje wiele punktów odpowiadających stanom stacjonarnym, a jej
krzywizna w ich otoczeniu niesie informację o widmie oscylacyjnym układu w danej
22
jedna z wad metody wynika z konieczności optymalizacji wszystkich struktur definiujących pierwsze
przybliżenie MEP, mimo że znajomość ich energii i geometrii nie wnosi nic do geometrii szukanego TS-u
23
trzy punkty definiujÄ…ce parabolÄ™ to substrat, produkt i TS1
76
__________Metodyka Badawcza
konfiguracji położeń jądrowych. Zakładając harmoniczne przybliżenie małych wychyleń,
możliwe jest wyznaczenie drgań normalnych24 układu definiujących jego strukturę
oscylacyjną za pomocą analizy zmian energii układu w momencie wychylania kolejnych
atomów z położeń równowagi25. Założenie harmoniczności potencjału sprowadza się do
aproksymacji dna studni potencjału każdego atomu parabolą, a zastosowanie tego podejścia
pozwala w rozwinięciu potencjału V(R) w szereg Taylora26 zaniedbać wyższe człony
rozwinięcia, sprowadzając wyrażenie na potencjał od wychylenia do postaci [4]:
ëÅ‚ öÅ‚
ëÅ‚ öÅ‚
V (R0 + x) = V (R0 ) +
(7.7)
"ìÅ‚ "V ÷Å‚ xi + 1 "ìÅ‚ "2V ÷Å‚ xi x j
ìÅ‚ ÷Å‚
÷Å‚
"xi
ij ij
íÅ‚ Å‚Å‚0 2 ìÅ‚ "xi"x j
íÅ‚ Å‚Å‚0
Dla powyższego problemu rozwiązanie równań Newtona pozwala wyznaczyć częstość
oraz amplitudę każdego z drgań. Analiza wibracyjna przeprowadzana jest w praktyce przez
kody obliczeniowe standardowo implementowane w pakietach do obliczeń kwantowo-
chemicznych.
Znajomość postaci drgań normalnych układu pozwala wygenerować teoretyczne
widmo oscylacyjne badanego modelu, które może zostać skonfrontowane z eksperymentem.
Pozwala to nie tylko na potwierdzenie geometrii stanu przejściowego (problem omówiony w
rozdziale 0) lecz także poprzez zastosowanie formalizmu termodynamiki atomowej
umożliwia obliczenie wpływu temperatury na entalpię swobodną badanego układu. Ostatni
problem zostanie szerzej omówiony w rozdziale 7.3.
24
drganie normalne to drganie harmoniczne o zdefiniowanej częstości angażujące wszystkie atomy układu -
wychylenia atomów charakteryzują się dla wybranego drgania identyczną fazą
25
każdy z atomów wychylany jest w trzech kierunkach w "prawo" i w "lewo"
26
dzięki założeniu małych wychyleń
77
7. Narzędzia aplikacyjne metod DFT_________
7.2 Struktura elektronowa
Dla zdefiniowanych metodami DFT geometrii reagentów charakterystycznych dla
badanego układu chemicznego wyznaczona struktura elektronowa (poziomy elektronowe,
przestrzenny rozkład gęstości elektronowej i potencjału elektrostatycznego, funkcja gęstości
stanów) niesie informację, którą wykorzystać możemy do interpretacji wyników. Analiza
kształtu orbitali dostarcza informacji o ich wiążącym/antywiążącym charakterze, zaś ich
energia pozwala określić np. wstępne deskryptory reaktywności chemicznej  potencjał
chemiczny i twardość chemiczną [79]. Dzięki zastosowaniu analizy populacyjnej możliwe
jest określenie atomowych ładunków cząstkowych, rzędów wiązań oraz innych istotnych
właściwości molekularnych takich jak momenty dipolowe czy parametry spektroskopowe.
Funkcja gęstości stanów wygenerowana dla modelu może zostać połączona z widmem
elektronowym dla badanego materiału. Analiza populacyjna i funkcja gęstości stanów zostaną
omówione w dalszej części rozdziału.
Analiza populacyjna - ładunki cząstkowe i rzędy wiązań
Dla wyznaczenia ładunków cząstkowych analiza populacyjna może opierać się na
podziale gęstości elektronowej na dwa zasadnicze sposoby [67]:
" w przestrzeni Hilberta (podział pomiędzy orbitale atomowe)
" w przestrzeni fizycznej (podział pomiędzy baseny27 lub przyczynki atomowe)
Przykładem analizy pierwszego typu jest analiza Mullikena [80,81] która jest jedną
z wielu metod bazujących na sumowaniu wkładów od poszczególnych orbitali atomowych,
w których ładunek qA na atomie A można określić jako:
x
qA = Z (S ) Pjk(S1-x)
A""" ij ik
(7.8)
i"A j k
gdzie ZA jest Å‚adunkiem jÄ…dra danego atomu, i, j, k sÄ… indeksami orbitali atomowych (jedna
z sum przebiega po orbitalach atomu A), P jest macierzą rzędów wiązań i ładunków, S jest
macierzą całek nakładania, a x jest pewną arbitralną liczbą z zakresu 0 do 1. Wadami tak
sformułowanej analizy są przede wszystkim ukryta definicja atomu w cząsteczce określonego
27
Basenem atomowym nazywamy obszar przestrzeni w układzie rzeczywistym ograniczony powierzchniami
stacjonarnymi zdefiniowanymi poprzez gradient gÄ™stoÅ›ci ( "Á = 0 ). WewnÄ…trz basenu gradient jest różny od
zera.
78
__________Metodyka Badawcza
przez lokalizowane na nim orbitale i zwiÄ…zane z tym problemy z funkcjami dyfuzyjnymi
i polaryzacyjnymi. Niewątpliwe zalety to prostota i stosowalność dla dowolnego typu funkcji
falowej w przybliżeniu LCAO. Do omawianej grupy analiz należy także analiza Löwdina
[82,83], która dzięki zastosowaniu bazy ortogonalizowanej jest stabilniejsza od analizy
Mullikena przy zwiększaniu liczby funkcji bazy.
Drugi podział charakteryzuje naturalną dla metod DFT analizę populacyjną Hirshfelda
[84], opierajÄ…cÄ… siÄ™ na caÅ‚kowitej gÄ™stoÅ›ci elektronowej28 ukÅ‚adu molekularnego Á0 oraz
atomowych gÄ™stoÅ›ciach elektronowych (Áa dla atomu a). W podejÅ›ciu tym wprowadza siÄ™
tzw. gÄ™stość odksztaÅ‚cenia Ád (ang. deformation density) zdefiniowanÄ… nastÄ™pujÄ…co:
Ád = Á0 (r) - (r - Ra ))
"(Áa
(7.9)
a
gdzie: Á0 jest gÄ™stoÅ›ciÄ… elektronowÄ… w czÄ…steczce
Áa(r - Ra) jest gÄ™stoÅ›ciÄ… elektronowÄ… izolowanego atomu umieszczonego w Ra
Za pomocą gęstości odkształcenia można obliczyć ładunki atomowe za pomocą wyrażenia:
qa = (r)Å"Wa (r)dr
(7.10)
d
+"Á
gdzie Wa(r) jest wkładem ze strony danego atomu do całkowitej gęstości w danym punkcie
(wagą statystyczną), zdefiniowaną następująco:
Áa (r - Ra )
Wa (r) =
(7.11)
(r - Rb )
"Áb
b
Zaletami analizy Hirshfelda są precyzyjna i rozsądna fizycznie definicja ładunków,
sÅ‚absza zależność wyników od bazy niż w metodzie Mullikena, oraz fakt ½
reprezentowalności gęstości. Wadą jest konieczność całkowania numerycznego.
Do grupy metod dzielących gęstość w przestrzeni rzeczywistej należą także analiza
Badera [85] i Vornoia [86]. Obie metody przypisują gęstość do basenów atomowych, przy
czym ta pierwsza proponuje bardziej realny wybór basenów atomowych w oparciu
o topologiczne cechy pola gęstości elektronowej [4].
28
w praktycznych realizacjach gęstość elektronowa konstruowana jest z orbitali, jednakże jest to tylko kwestia
implementacji
79
7. Narzędzia aplikacyjne metod DFT_________
Rzędy wiązań (indeksy wartościowości) są znacznie trudniejsze do zdefiniowania od
ładunków przede wszystkim z braku klasycznego odpowiednika dla bardziej
skomplikowanych przypadków, lecz ich przydatność (ocena efektów wiążących między
parami atomów w stanach przejściowych, przewidywanie reaktywności atomu w cząsteczce)
są motorem dla opracowywania metod ich określania. Możemy wyróżnić dwa sposoby
podejścia do problemu:
" podejście bezwzględne oparte na strukturze elektronowej samej molekuły,
" podejście różnicowe wykorzystujące analizę zmian struktury elektronowej podczas
formowania czÄ…steczki.
Szeroko stosowanym przedstawicielem pierwszego z podejść jest analiza Meyera
w której rząd wiązania pomiędzy atomami A-B definuje się jako [87,88]:
BOA-B = [(PÄ… S) (PÄ… S) + (P² S) (P² S) ]
"" µ½ ½µ µ½ ½µ
(7.12)
µ"A½"B
gdzie gdzie PÄ… i P² sÄ… macierzami gÄ™stoÅ›ci dla spinu Ä… i ².
Zaletą formalizmu Meyera jest zgodność wyznaczonych BO dla modelowych rzędów
wiązań w prostych modelach29 oraz prostota i popularność samej metody. Wadą jest opis
tylko kowalencyjnej składowej oddziaływania międzyatomowego.
Opierające się na pojęciu promolekuły30 podejście różnicowe wyraża rząd wiązania
przez zmiany w macierzy gęstości związane z tworzeniem molekuły z promolekuły.
Opracowana metoda Nalewajskiego-Mrozka [89,90,91] dla obliczeni rzędu wiązania wymaga
zdefiniowania indeksów kowalencyjnych i jonowych31 pochodzących od oddziaływań jedno
i dwuatomowych. Metoda poprawnie wyznacza rzędy wiązań nawet o znacznej jonowości
a dodatkowo dostarcza nie tylko sumarycznego rzędu wiązania ale i poszczególnych wkładów
(jonowy  kowalencyjny, jednoatomowy  dwuatomowy). Wadą jest silna zależność od
arbitralnie definiowanej premolekuły.
29
rząd wiązania kowalencyjnego wynosi tyle ile mamy par wiążących (Lewis and Kossel), dla cząsteczek 2
atomowych: rząd = (Nb - Na)/2 gdzie Nb liczba elektronów na orbitalach wiążących, Na  na antywiążących []
30
układ złożony z izolowanych fragmentów przesuniętych w miejsca ich położeń w cząsteczce (nazwa ang.
Separated Atoms Limit, SAL)
31
indeksy kowalencyjne wywodzą się z części wymiennej zaś jonowe z kulombowskiej
80
__________Metodyka Badawcza
Funkcja gęstości stanów
Funkcja gęstości stanów (z ang. density of states, DOS) dla danego pasma n wyraża
liczbę stanów przypadającą na jednostkę energii i zdefiniowana jest wzorem:
r
dk
(7.13)
Nn(E) = ´ (E - En (k))
+" 3
4Ä„
gdzie En(k) opisuje dyspersję danego pasm. Całkowanie we wzorze (7.13) odbywa się po
strefie Brillouine'a. Całkowita funkcja gęstości stanów jest otrzymywana dzięki całkowaniu
po wszystkich pasmach, i daje całkowitą liczbę elektronów w modelu. W układach
spolaryzowanych spinowo możliwe jest wprowadzenie niezależnych funkcji DoS dla
elektronów Ä… i ²32 [92]. Funkcja DoS jest bardzo użytecznÄ… koncepcjÄ… pozwalajÄ…cÄ… na Å‚atwÄ…
wizualną analizę struktury elektronowej badanego układu. Szerokość pasma przewodzenia,
wartość przerwy energetycznej czy położenie charakterystycznych poziomów lokalnych
dostępne jest wprost z wykresu funkcji DoS, co pozwalają na jakościową interpretację
i korelację danych obliczeniowych z danymi doświadczalnymi.
7.3 Termodynamika atomistyczna
Teoretyczne podejście do katalizy heterogenicznej ma w pierwszym przybliżeniu
poprawnie opisać centra aktywne odpowiedzialne za proces katalityczny na poziomie
molekularnym a następnie wyjaśnić sposób oddziaływania (adsorpcji, aktywacji) reagentów
z powierzchnią katalizatorów. Określenie struktury stanów przejściowych pomiędzy
reagentami i produktami modelowanego procesu to kolejny krok prowadzÄ…cy do zrozumienia
zachodzących na katalizatorze procesów. Związana z tym problemem eksploracja statycznej
powierzchni energii potencjalnej (opisana w rozdziale 7.1) pozwala na zdefiniowanie
geometrii i energii stanów stacjonarnych w zerowej temperaturze i bez uwzględnienia
ciśnienia parcjalnego reagentów. Elektronowa energia (Eelec) uzyskana z obliczeń DFT nie
niesie więc informacji o wpływie temperatury na energetykę układu oraz zaniedbuje
niezależną od temperatury energię drgań zerowych (ZPE, z ang. Zero Point Energy).
Naturalnym postępowaniem jest więc próba uwzględnienia w wynikach parametrów
termodynamicznych (p,T). Tego typu modelowanie realizować można na bazie
termodynamiki statystycznej, wykorzystując pojęcie cząsteczkowej funkcji rozdziału [93].
Podejście to pozwala na uzyskanie wglądu w kształt powierzchni entalpii swobodnej (G),
32
suma takich funkcji odtwarza całkowity DoS zaś różnica definiuje spinową funkcję gęstości stanów
81
7. Narzędzia aplikacyjne metod DFT_________
która w jednoznaczny sposób definiuje stabilność układu przypisanego danemu punktowi na
analizowanej powierzchni (w warunkach izobarycznych) [94] oraz jest wprost powiÄ…zana z
kinetyką reakcji poprzez teorię stanu przejściowego Eyringa, Evansa i Polanyiego [95,96].
Dodatkowo oprócz przewidywania stabilności i kierunków przemian w układzie możliwe
staje się także generowanie diagramów termodynamicznych określających sposób adsorpcji
powierzchniowej w funkcji temperatury i ciśnienia.
Termodynamika atomistyczna  podstawowe wzory
Najważniejszym pojęciem termodynamiki statystycznej jest funkcja rozdziału (suma
stanów) q, która wywodząc się z rozkładu Boltzmanna podaje średnią liczbę stanów
dostępnych termicznie dla cząsteczki w temperaturze jaką ma układ [93]. Dla ogólnego
przypadku funkcja ta wyraża się wzorem33:
-µ kT
i
q =
"e
(7.14)
i
Gdzie:
k -stała Boltzmanna
µ -energia stanu liczona wzglÄ™dem stanu najniższego
Energię cząsteczkową można rozłożyć na następujące udziały: energię oscylacji,
translacji, rotacji i przejść elektronowych. W przypadku gdy energia jest sumą niezależnych
składowych, funkcja rozdziału jest iloczynem funkcji rozdziału odpowiadających tym
przyczynkom. Możemy więc zdefiniować osobno translacyjną, rotacyjną, oscylacyjną i
elektronową funkcję rozdziału cząsteczki (odpowiednio: qtrans, qrot, qosc, qel), a wypadkową
funkcję przedstawić jako ich iloczyn. Znaczenie funkcji rozdziału polega na tym, że zawiera
ona całą informacje potrzebną do określenia właściwości termodynamicznych układu
niezależnych cząstek. Opierając się na klasycznych wzorach [93], przy założeniu przybliżenia
gazu doskonałego, poprzez funkcje rozdziału możemy zdefiniować przyczynki termiczne
pozwalające wyliczać interesujące z punktu modelowania w katalizie heterogenicznej funkcje
stanu [97]. Wzory definiujące odpowiednie funkcje rozdziału oraz sposób otrzymania
wkładów termodynamicznych do energii i entropii przedstawiono w Aneksie I. Zdefiniowane
tak wkłady pozwalają wyznaczyć potrzebne w obliczeniach entalpii swobodnej funkcje stanu
układu co omówione zostanie w oparciu o wyrażenia zawarte w tabeli 7.1. Energię
33
poniższy wzór stosuje się dla przypadku stanów niezdegenerowanych
82
__________Metodyka Badawcza
wewnętrzną (U)34 definiujemy jako sumę energii elektronowej oraz termodynamicznych
wkładów do energii. Entalpię jednego mola gazu doskonałego (H) otrzymujemy dodając do
energii wewnętrznej człon pV = RT. Entropia układu w zadanej temperaturze (S) to suma
wkładów od poszczególnych rodzajów ruchu zaś entalpia swobodna może ostatecznie zostać
zdefiniowana dzięki znajomości powyżej opisanych funkcji termodynamicznych za pomocą
ostatniego równania tabeli 7.1.
Tabela 7.1: Funkcje stanu układu dostępne z obliczeń termodynamiki atomistycznej
Energię wewnętrzną U U(T) = ZPE + Eosc(0 T), + Erot + Etrans
Entalpię H H(T) = U(T) + RT (założenie: pV = RT)
EntropiÄ™ S S(T) = Sosc + Srot + Strans
EntalpiÄ™ swobodnÄ… G G(T) = H(T) - TS(T)
Tak zdefiniowana entalpia swobodna zawiera więc w sobie energię elektronową
(Eelec), energię drgań zerowych (ZPE) jako przyczynki niezależne od temperatury a także
termicznie zależne wkłady do entropii i energii. Tak zdefiniowana wartość jest standartowo
obliczana poprzez programy obliczeniowe chemii kwantowej (miedzy innymi Mopac,
Gaussian, DMol3) i pozwala wprost definiować względną stabilność modelowanych układów
z uwzględnieniem warunków termodynamicznych, oraz szacować wartości stałych szybkości
zachodzenia rozważanych reakcji, zgodnie z termodynamiczną postacią podstawowego
równania teorii stanu przejściowego [95,98]:
#
ëÅ‚ öÅ‚
kBT - "G
ìÅ‚ ÷Å‚
k = expìÅ‚ ÷Å‚
(7.15)
h RT
íÅ‚ Å‚Å‚
gdzie "G# to entalpia swobodna aktywacji stanowiąca różnicę w entalpii swobodnej
pomiędzy reagentami a układem w TS i określająca ich stosunek liczbowy w stanie
stacjonarnym [99].
34
pogrubiona czcionka dla symboli termodynamiczny ułatwia odróżnienie ich od zdefiniowanych w Aneksie
składowych (Sosc + Srot + Strans)
83
7. Narzędzia aplikacyjne metod DFT_________
Termodynamiczne diagramy adsorpcji
Jedynie w wypadku porównywania układów o tym samym składzie atomowym ich
względna stabilność może być wyrażona wprost jako różnica stacjonarnych energii
elektronowych czy też na poziomie termodynamiki atomistycznej jako różnica odpowiednich
funkcji termodynamicznych. Dla układów które nie są zgodne względem składu chemicznego
takie porównanie zawodzi [100]. Uszeregowanie stabilności takich układów jest możliwe na
bazie założeń termodynamiki atomistycznej przy wprowadzeniu rezerwuaru atomów które
dopełniają stechiometrię modelowanego układu, prowadząc do sytuacji w której za każdym
razem porównywane są modele o identycznym składzie chemicznym.
Definicja "magazynu" atomów nie jest jednoznaczna i możemy jako odnośnika użyć
wyników uzyskanych dla atomów izolowanych lub połączonych w cząsteczki, będących
w fazie gazowej lub budujących ciało stałe. Na szczęście takie ogólne rozważania zostają
ułatwione w przypadku analizy konkretnego problemu. Dla modelowania stabilności
roztworów stałych porównanie w większości przypadków będzie polegać na odniesieniu do
ciałostałowych układów zaś w przypadku analizowanej w niniejszej pracy adsorpcji
molekularnej zachodzącej na powierzchni katalizatora rezerwuarem atomów będzie faza
gazowa o odpowiednim składzie cząsteczkowym pozostająca w równowadze z modelowaną
powierzchnią [45, 46,101], w określonej temperaturze i ciśnieniu reagentów (Rys. 7.2).
Rys. 7.2: Schematyczne przedstawienie powierzchni (model slabowy) w kontakcie z modelem referencyjnym
(cząsteczki w fazie gazowej). Oba modele mogą wymieniać atomy w procesach adsorpcji i desorpcji, jednak
całkowita liczba atomów danego typu i pozostaje stała (Ni(tot) = const.)
84
__________Metodyka Badawcza
Przy takich założeniach Nitot niech oznacza całkowitą liczbę cząsteczek adsorbenta i,
Niads liczbÄ™ czÄ…steczek zaadsorbowanych zaÅ› Nigas liczbÄ™ czÄ…steczk w fazie referencyjnej fazie
gazowej. Założenie iż układ znajduje się w równowadze pozwala traktować sumę entalpii
swobodnych pochodzÄ…cych od modelu slabowego (Gslab) i rezerwuaru czÄ…stek w fazie
gazowej jako minimum funkcji G(T,p,N) [100]:
ads
slab gas
G (T , p,{N })+ Gigas (T , p, N )= min
i " i
(7.16)
i
Zapisując entalpię swobodną fazy referencyjnej (gazowej) za pomocą potencjałów
chemicznych:
gas
Gigas (T , p, Nigas )= (Nitot - Niads )µ (T , p)
(7.17)
i
możemy w związku założeniem stanu równowagi35 zapisać powierzchniową entalpię
swobodną układu jako [100]:
1 ëÅ‚
ads
G(T, p)= ìÅ‚Gslab(T, p,{Niads})- µ(T, p)öÅ‚
÷Å‚
(7.18)
"Ni
A
íÅ‚ i Å‚Å‚
gdzie A jest powierzchniÄ… modelu slabowego normalizujÄ…cÄ… wynik.
Wyznaczone minimum powyższej funkcji dla założonych T i p zapewnia iż
zlokalizowany zostaje taki podział atomów pomiędzy fazy referencyjną i powierzchniową
(Niads), że mamy do czynienia ze stanem równowagi układu. Dla większych układów takie
podejście jest niepraktyczne ze względu na kosztowne obliczenia wymagające wyznaczenia
energii elektronowej (Eel) i wkładów termicznych (Eosc(0 T), Erot, Etrans Sosc, Srot, Strans
,
omówione w poprzednim podrozdziale) a następnie znalezienia zależności G od Niads.
Podejściem stosowanym w praktyce jest więc obliczenie dla serii układów o różnej strukturze
(pokrycie reagentami, rodzaj eksponowanych płaszczyzn) zależności od potencjału
chemicznego układu (a więc pośrednio od temperatury i ciśnienia parcjalnego reagentów).
Dla zdefiniowanych zmiennych (µ lub T i p) możemy wówczas przewidzieć, który
z zakładanych modeli adsorpcji charakteryzuje się najniższą entalpią swobodną, a więc jest
najbardziej stabilny w danych warunkach. Oczywistym wynikającym z powyższych rozważań
ograniczeniem metody jest problem adekwatności zaproponowanych modeli adsorpcji.
Obliczenia prowadzą do wybrania układu najstabilniejszego spośród wszystkich
zdefiniowanych, jednak istnieje ryzyko że układ najbardziej stabilny z przyjętych do obliczeń
nie odpowiada rzeczywistości. Niemniej liczne prace badawcze oparte na takim podejściu
35
w stanie równowagi potencjały chemiczne fazy referencyjnej i modelowanej powierzchni są identyczne
85
7. Narzędzia aplikacyjne metod DFT_________
"nie wprost" dowodzą iż przy zachowaniu systematyczności i wykorzystaniu intuicji
chemicznej możliwe jest uwzględnienie takich struktur które pozostają w dobrej zgodności
z eksperymentem [94,101].
W praktyce oblicza się nie bezwzględną entalpię swobodną układu lecz jedynie jej
zmianÄ™ zwiÄ…zanÄ… z modelowanym procesem:
"rG = G(powierzchnia pokryta) - (G(powierzchnia czysta) + G(faza gazowa)) (7.19)
Zakładając, że różnicę pomiędzy entalpiami swobodnymi ciał stałych możemy przybliżyć
jako różnicę ich energii elektronowej (Eel dostępna wprost z obliczeń DFT):
ort 36
"rG = "Eel + Niads"(ZPE + TS + Eort )
(7.20)
gdzie:
"Eel = Eel  ( Eel + Eel )
(powierzchnia pokryta) (powierzchnia czysta) (faza gazowa)
zaś drugi człon równania (7.20) odnosi się do zmiany energii Gibbsa układu w wyniku
adsorbcji Niads cząsteczek z której wyłączono część związaną z energią elektronową (człon
pierwszy) i w dalszej części bÄ™dziemy go oznaczać w skrócie symbolem "µ.
Korzystając z definicji potencjału chemicznego cząsteczki [93] i możemy zapisać:
"µi = "µi0(T) + RT‡ln(pi/po) (7.21)
gdzie standardowy potencjaÅ‚ chemiczny molekuÅ‚ w fazie gazowej ("µ0(T)) dany jest wzorem
wykorzystującym poszczególne przyczynki omówione w poprzednim rozdziale:
"µ0(T) = ZPE + Eosc(0 T) + Erot + Etrans + RT - T(Sosc + Srot + Strans) (7.22)
Zauważmy że wykorzystanie wzoru (7.21) prowadzi do liniowej zależności pomiędzy
"µ i ln(pi/po). Zależność temperaturowa "µ0(T) jest bardziej skomplikowana i zależy od
postaci cząsteczkowej funkcji rozdziału [97]. Dodatkowo w tak zdefiniowanym podejściu nie
musimy znać wprost energii powierzchniowej badanego układu w przybliżeniu energii
stacjonarnej. Różnice entalpii swobodnej ("rG) wyznaczane są względem modelu odniesienia
który dla problemu adsorpcji wybierany jest zazwyczaj jako powierzchnia niepokryta
reagentami.
36
Sort = Sosc + Srot + Strans, oraz analogicznie Eort = Eosc + Erot + Etrans
86
__________Metodyka Badawcza
7.4 Energie powierzchni i morfologia krystalitów
Powierzchnie ciała stałego są z definicji niekorzystne energetycznie, gdyż tworzenie
powierzchni wymaga wykonanie pewnej pracy związanej ze zrywaniem wiązań
międzyatomowych. Dla powierzchni o różnym wskaz
niku Millera istnieje bliski i intuicyjny
związek miedzy swobodną energią powierzchni a jej strukturą. W ogólności najbardziej
stabilnymi (korzystnymi energetycznie) powierzchniami sÄ… te charakteryzujÄ…ce siÄ™ wysokimi
gęstościami atomów powierzchniowych i wysokimi liczbami koordynacyjnymi terminalnych
atomów.
Anizotropia powierzchniowej energii swobodnej jest ściśle związana
z równowagowym kształtem krystalitów [102] co wykazał Gibbs [103]. Suma energii
powierzchniowych krystalitu o pewnej objętości musi być minimalna zgodnie z wyrażeniem:
(7.23)
+"Å‚dA = min
gdzie ł oznacza energią powierzchniową zaś A to powierzchnia. Kształt kryształu może być
w prosty sposób określony dzięki zastosowaniu konstrukcji Wulffa [104]. Warunek
minimalnej energii powierzchniowej kryształu można określić korzystając z zależności:
Å‚
hkl
(7.24)
hhkl = cons.
gdzie łhkl to energia swobodna powierzchni o kierunku oznaczonym hkl zaś hhkl to odległość
od środka kryształu do powierzchni hkl.
Konstrukcja Wulffa definiuje kształt kryształu jako wewnętrzny obrys wielościanu
utworzonego na bazie powierzchni prostopadłych do wektorów zaczepionych w środku
układu, wskazujących kierunek narastania powierzchni, których długość związana jest
z energią powierzchniową. Schemat generowania kształtu Wulffa dla przypadku
dwuwymiarowego przedstawiono na rysunku (Rys. 7.3). Ostatecznie znajomość symetrii
układu oraz energii powierzchniowych dla różnych terminacja kryształu pozwala określić
jego morfologiÄ™ (Rys. 7.4).
87
7. Narzędzia aplikacyjne metod DFT_________
Rys. 7.3: Konstrukcja Wulffa dla przypadku dwuwymiarowego. Obwiednia w kolejnych ćwiartkach układu
współrzędnych została wygenerowana zgodnie z symetrią układu.
Rys. 7.4: Kształt trójwymiarowych kryształów (symetria kubiczna Oh1, Fm3m) określony na podstawie
konstrukcji Wulffa. a) energia powierzchni 100 znacznie niższa od pozostałych b) energia powierzchni 111
niższa niż 100
Dla określania energii powierzchniowych katalizatora wykorzystać można obliczenia
ab initio prowadzone w modelu periodycznym (omówione w rozdziale 12). Dodatkowo
rozważenie adsorpcji cząsteczek w połączeniu z podejściem termodynamiki atomistycznej
(rozdział 7.3) na rozważanych powierzchniach pozwala przewidzieć także zmiany morfologii
materiału związane ze zmianami parametrów termodynamicznych (temperatura, ciśnienia
parcjalne reagentów) które wpływają na swobodna energię powierzchni [100,101, 46].
88