Fizyka, technologia oraz modelowanie
wzrostu kryształów
Stanisław Krukowski i Michał Leszczyński
Instytut Wysokich Ciśnień PAN
01-142 Warszawa, ul Sokołowska 29/37
tel: [phone]
e-mail: [e-mail], [e-mail]
Zbigniew Żytkiewicz
Instytut Fizyki PAN
02-668 Warszawa, Al. Lotników 32/46
E-mail: [e-mail]
Wykład  2 godz./tydzień  piątek 13.15  15.00
Interdyscyplinarne Centrum Modelowania UW
Budynek Wydziału Geologii UW  sala 3089
http://www.unipress.waw.pl/~stach/wyklad_ptwk_2007/
Modelowanie procesów wzrostu  obraz mikro
Stanisław Krukowski
Modelowanie procesów wzrostu  dwa obrazy
" Modelowanie - metody
" Modelowanie - zagadnienia i procesy fizyczne
Modelowanie procesów wzrostu  zagadnienia
" Modelowanie w skali makroskopowej
- obrazowanie procesów transportu podczas wzrostu
kryształów (masy, energii pędu)
- wyznaczanie naprężeń w strukturach niejednorodnych
- wyznaczanie własności elektrycznych układów
elektronicznych
- wyznaczanie własności optycznych układów elektronicznych
" Modelowanie w skali atomowej
- wyznaczanie struktury i morfologii kryształów
- wyznaczanie charakterystycznych własności energetycznych
dla stanów równowagowych
- wyznaczanie charakterystycznych własności energetycznych
dla procesów kinetycznych
Modelowanie procesów wzrostu  metody
" Modelowanie w skali makroskopowej
- metoda skończonej różnicy
- metoda skończonej objętości
- metoda elementu skończonego
" Modelowanie w skali atomowej
- metoda Monte Carlo
- metoda dynamiki molekularnej
- metody ab initio - DFT
Metoda Monte Carlo
" Określenie przestrzeni zdarzeń elementarnych
" Definicja zmiennej losowej
" Wyznaczenie rozkładu prawdopodobieństwa zmiennej
losowej
" Próbkowanie rozkładu
" Wyznaczenie własności fizycznych
Metoda Monte Carlo  określenie przestrzeni zdarzeń
Aksjomatyczna definicja prawdopodobieństwa
Prawdopodobieństwo - miara p, zdefiniowana na zbiorze (algebrze)
zdarzeń losowych {Ai, i = 1 2 3 n}, spełniająca następujące zależności:
" p("
") = 0, p(E) = 1
"
"
" 0 d" p(A) d" 1
d" d"
d" d"
d" d"
zdarzenia się wykluczają tzn. gdy A )" B = "
)" "
" p(A*" )" "
*"B) = p(A) + p(B) )" "
*"
*"
Prawdopodobieństwo zdarzenia warunkowego, tzn. zajścia zdarzenia A pod
warunkiem że zaszło zdarzenie B jest równe:
p(A )"B)
p(A|B) =
p(B)
Zmienna losowa
" Załóżmy, że algebra zdarzeń losowych A jest odwzorowywana na
funkcję o wartościach rzeczywistych. Funkcję taką nazywamy
zmienną losową
X :{A}�! R
" Jednocześnie na algebrze zdarzeń losowych {A} określone jest
prawdopodobieństwo, tzn. funkcja rzeczywista spełniająca
warunki definicji aksjomatycznej.
p :{A}�! [0,1]
Każdej wartości zmiennej losowej X można przyporządkować
prawdopodobieństwo odpowiadające sumie prawdopodobieństw
wykluczających się zdarzeń dla których zmienna losowa przyjmuje wartość X.
Zależność prawdopodobieństwa p od wartości zmiennej losowej nazywamy
rozkładem zmiennej losowej.
Własności rozkładów zmiennej losowej
" dodatnio określony
P(x)e" 0
" unormowany
"P(x )= 1
i
i
" wyznaczany z pomiarów  np. rozkład energii o szerokości "E
nk
P(Ek )=
N
gdzie nk  liczba cząstek o energii w przedziale Ek,Ek+"E , N  całkowita
liczba cząstek.
" Rozkłady ciągłe i dyskretne
Rozkład zmiennej losowej nazywamy dyskretnym, gdy jego zakres
jest skończony lub policzalny.
Zależności pomiędzy
różnymi rozkładami
prawdopodobieństw
Próbkowanie  sampling
Dany jest generator liczb losowych o rozkładzie jednorodnym: tzn.
1 � "[0,1]
ńł
f(�) = �ł0 �
p(� ) = f(�) d�
gdzie
i
+"
"[0,1]
ół
Należy znalez
ć procedurę generacji liczb losowych o rozkładzie
prawdopodobieństwa danym funkcją rozkładu (próbkowanie)
p x " &! = f(x) dx
( )
k
+"
x " &!
Warunek zgodności - warunek równości prawdopodobieństw:
�
x x
p(x > y) = p(� > � )
� = f(y) dy
o
+"f (y)dy = +"d �
+"
- "
-" 0
" Rozkład Lorentza
x
Ą
1 1 1 1
x = tan�ł 2� - 1
�ł ( )�ł
�ł
dy = �
f(x) =
+" 2
2
�ł łł
2
Ą 1 + y
Ą 1 + x - "
Techniki odrzucania
f(x)
Załóżmy że mamy funkcję rozkładu
jednej zmiennej określoną na
przedziale [0,1] za pomocą
skomplikowanej zależności:
0 1 x
Próbkowanie tego rozkładu za pomocą technik odrzucania:
" generacja liczby losowej �1 z przedziału [0,1]
�
�
�
" generacja liczby losowej �2 z przedziału [0,1]
�
�
�
" sprawdzamy warunek
� < f(� ) �! x = �
2 1 1
" jeżeli warunek nie jest spełniony, wykonujemy procedurę generacji
liczby od początku
Wada metody - dla pewnych rozkładów prawdopodobieństw jest ona
stosunkowo mało efektywna.
Algorytm próbkowania Metropolisa
Algorytm ten został skonstruowany do wykonywania symulacji
własności równowagowych cieczy (Metropolis et al , J. Chem. Phys.
21 (1953) 1087
Równowagowy rozkład prawdopodobieństwa dla układu o
temperaturze T opisany jest za pomocą zespołu kanonicznego Gibbsa:
E p , q
�ł ( )�ł 1 V(qi)�ł
1
�ł
i i
f(qi) = exp�ł-
f p , q = exp�ł-
( ) �ł
�ł
i i
Z kT
Z �ł kT łł
�ł łł
Własności rozkładu:
- rozkład określony na przestrzeni bardzo wielu zmiennych
- rozkład posiadający bardzo ostre maksimum w pobliżu stanu
równowagi makroskopowej.
Własności te powodują że standardowe metody jego próbkowania są w
zasadzie niemożliwe do wykonania.
Algorytm próbkowania Metropolisa - procedura
" Wybieramy zespół zmiennych początkowych q(t)
" Znajdujemy stan q'(t+1)
" Oszacujemy prawdopodobieństwo znalezienia się w stanie q'(t+1)
" Używamy zasady równowagi szczegółowej
p(q|q') f(q) E(q) - E(q')�ł
Q = = = exp�ł-
�ł �ł
�ł łł
p(q'|q) f(q') k T
" Jeżeli Q > 1 przejście zachodzi
" Jeżeli Q < 1 generujemy liczbę losową �, gdy � < Q przejście
� �
� �
� �
również zachodzi, jeżeli � > Q to próbę odrzucamy.
�
�
�
Algorytm Metropolisa jest więc rodzajem błądzenia przypadkowego w
przestrzeni mikrostanów.
Funkcja korelacji gęstości  faza gazowa i ciekła
Faza gazowa
Faza ciekła
Funkcja korelacji gęstości  faza stała
Niższa gęstość
Wyższa gęstość
Metoda dynamiki molekularnej - numeryczna
analiza zachowania układów dynamicznych
"H
"H
pi = -
&
Formalizm Hamiltona
qi =
&
"qi
"pi
2
pi
gdzie
H(p,q) =
+ V({qi})
"
2mi
i
Dla sił zachowawczych tzn. dla potencjału niezależnego od prędkości
można wprowadzić zależność pomiędzy pędem i prędkością vi = pi/m:
"
&
qi = vi
mivi = Fi
&
Fi = - V({qi})
"qi
Można sformułować w postaci jednego równania:
Fi({q})
&&i
q = fi({q})=
mi
Obliczanie potencjału sił  mechanika kwantowa
Problem  rozwiązanie równania Schr�dingera na funkcję falową �(R,r):
�
�
�
"�(R, r)
ih
= H �(R, r)
"t
R oraz r  współrzędne położeń jąder atomowych i elektronów.
Przybliżenie Borna-Oppenheimera  zaniedbanie ruchu jąder przy
rozwiązywaniu równania Schr�dingera.
"�R (r)
ih
= H(r;R) �R (r)
�(R, r) = f(R)�R(r)
"t
Warunkiem stosowania metody DM jest aby długość fali de Broglie a
ruchu jąder była znacznie mniejsza niż średnia odległość pomiędzy atomami
1
3
1
h V
�ł �ł
3
= � << v =
�ł �ł
2ĄmkBT N
�ł łł
Potencjał Lennarda-Jonesa: gazy szlachetne
ńł
�ł�ł qij �ł-12 �ł qij �ł-6 łł
" Potencjał Lennarda-
�ł
�ł �ł
śł q < rc
4��ł�ł �ł - �ł �ł
�ł �ł
(r )
V qij =
�ł �
�ł śł
Jonesa �ł łł
�ł�ł � łł �ł
�ł
0
q > rc
ół
1,00
0,75
0,50
0,25
0,00
-0,25
-0,50
-0,75
-1,00
0,8 1,0 1,2 1,4 1,6 1,8 2,0 2,2 2,4 2,6 2,8 3,0
r/�
V/
�
Anizotropowy potencjał azotu N2
ŃB
Z
0.0 0E+000
Ś
�B
-5.00 E-022
-1.00 E-021
Y
0 0 0 0 0 0
R
-1.50 E-021
0 0 90 0 0 0
0 0 90 0 90 0
0 0 90 90 90 0
ŃA
0 0 45 0 45 0
-2.00 E-021
0 0 45 0 135 0
�A
Ś
-2.50 E-021
0.30 0 .35 0.40 0.45 0.50 0.55 0. 60
X
R [nm]
Hartree-Fock approx. (van der Avoird et
al. - 1984)
E [J]
Anizotropowy potencjał azotu N2  CCSD(T)
4.00E-021
�A=0
�A=90
�B=0
�B=0
4.00E-021
ĆB=0
ĆB=0
2.00E-021
2.00E-021
0.00E+000
0.00E+000
-2.00E-021
0.40 0.45 0.50 0.55 0.60 0.65
0.35 0.40 0.45 0.50 0.55
R (nm)
R(nm)
�A=90
�A=90
4.00E-021 �B=90
�B=90
ĆB=0
1.50E-020
ĆB=90
2.00E-021
1.00E-020
5.00E-021
0.00E+000
0.00E+000
-2.00E-021
0.30 0.35 0.40 0.45 0.50 0.55
0.35 0.40 0.45 0.50 0.55 0.60
R(nm)
R(nm)
(P. Strąk - 2006)
E(J)
E(J)
E(J)
E(J)
Obliczanie sił  najbardziej pracochłonny element metody
dynamiki molekularnej
Potencjał oddziaływania
V(q)= V(R) = �R (r) H(r;R)�R (r)
Substancje  potencjały oddziaływań
1. Kryształy i ciecze pierwiastków gazów szlachetnych -
oddziaływania dwucząstkowe Lennarda-Jonesa (krótkozasięgowe)
2. Kryształy i ciecze jonowe  oddziaływania dwucząstkowe
kulombowskie (długozasięgowe)
3. Kryształy półprzewodników  oddziaływania trójcząstkowe
(krótkozasięgowe)
4. Kryształy i ciecze metali prostych  oddziaływania kolektywne
(krótkozasięgowe)
Metody całkowania równań ruchu
Rozwiązanie problemu  zastąpienie ewolucji ciągłej ewolucją o
skończonym przedziale czasowym. Równanie różniczkowe
&&
q = f(q)
całkujemy w skończonej różnicy czasów oznaczając wynik w n-tym kroku
czasowym przez xn.
Dwa typy metod całkowania równań ruchu:
1. Metody otwarte  metody predyktora
Oznacza ona że wartość xn+1 jest wyznaczona wyłącznie przez wartości
w chwilach poprzednich tzn. xn, xn+1 , itd
2. Metody zamknięte  metody predyktor+korektor
Oznacza ona że wynik yn+1 jest wyznaczony poprzez zależność od xn,
xn-1, tzn. jak poprzednio, lecz następnie wartość tę koryguje się
poprzez użycie tej samej zależności dla xn+1
Metody otwarte: metoda Eulera
Metoda Eulera polega na rozwinięciu w szereg Taylora względem
czasu dla położeń q oraz prędkości v
1 1
&& &&n
xn+1 = xn + xnh + x h2 = xn + vnh + hfn
2 2
vn+1 = vn + hvn = vn + hfn
gdzie fn = f(xn)
&
Zachowanie metody Eulera jest bardzo złe  jest to metoda względnie
mało stabilna
Metody zamknięte: zmodyfikowana metoda Eulera
Predyktor
1
1
*
yn+1 = xn + hvn + h2fn fn = (fn + f(yn+1))
2
2
Korektor
1 1
*
xn+1 = xn + hvn + h2fn = yn+1 + [f(yn+1)- f(yn )]
2 4
*
fn+1 = f xn+1
( )
vn+1 = vn + hfn
Wadą zmodyfikowanej metody Eulera jest obliczanie siły dwukrotnie
w każdym kroku czasowym co wymaga dużych mocy obliczeniowych
 jest to jednak metoda bardziej stabilna
Metody otwarte: metoda Verleta
Metoda Verleta polega na użyciu rozwinięcia względem czasu w dwu
krokach czasowych
"x
x' a"
n
1 1
gdzie
"t
x = x - hx' + h2x'' - h3x'''
n-1 n n n n
2 6
"2x
x'' a"
n
"t2
1 1
xn+1 = xn + hx' + h2x'' + h3x'''
n n n "3x
2 6
x''' a"
n
"t3
Po dodaniu stronami otrzymujemy niezależny od prędkości algorytm
xn+1 = 2xn - xn-1 + h2f (xn )
Algorytm dla prędkości
xn+1 - xn-1
( )
vn =
2h
Własności metody Verleta - prosty, szybki algorytm, dokładny do
rzędu O(h4)  powodują że jest to metoda bardzo często stosowana.
Oscylator harmoniczny:
porównanie różnych metod
całkowania równań ruchu -
położenie
Względne odchylenie od
rozwiązania dokładnego
Steps per cycle  na ile kroków
został podzielony jeden okres
drgań oscylatora - i
T
h =
i
Oscylator harmoniczny:
porównanie różnych metod
całkowania równań ruchu 
dryft energii
Dopasowanie metoda
najmniejszych kwadratów
energii  cześć proporcjonalną
do długości kroku nazywamy
dryftem
Steps per cycle  na ile kroków
został podzielony jeden okres
drgań oscylatora - i
T
h =
i
Potencjały zależne od kierunku  hybrydyzacja sp3
Krzem  potencjał oddziaływania Stillingera-Webera:
" Oddziaływania dwucząstkowe
1
ńłA -p -q
łł
rij�ł
�ł exp�ł
( )
�łr - aśł r < a
f2 r = �ł Br - r
v2 rij = � f2�ł �ł ( )
( ) �ł
�ł �ł
�ł łł
�
�ł
0 r > a
ół
Założono że głębokość potencjału jest określona przez parametr, natomiast
- jest wyznaczona przez żądanie aby f2(21/6) = 0.
" Oddziaływania trójcząstkowe
f3(ri, rj,rk)= h(rij,rik,�jik)
�ł ri rj rk �ł
v3 ri , rj , rk = �f3 , ,
�ł �ł
( )
( ) (
+ h rji, rjk,�ijk + h rki,rkj,�ikj )
� � �
�ł łł
2
ńł
�ł łł
ł ł 1
�ł
 exp�ł + cos�ł�jik + rij < a i rik < a
�ł
�ł �ł
śł
h(rij,rik ,�jik)=
3
�ł
�ł łł
�łr - a rik - a śł
ij
�ł �ł
�ł
0 rij > a lub rik > a
ół
Potencjały Stillingera-Webera  wybór parametrów (Si)
Krzem  potencjał oddziaływania Stillingera-Webera:
" A = 7.049 556 277
" B = 0.602 224 5584
" p = 4
" q = 0
" a = 1.80
"  = 21.0



" ł = 1.20
ł
ł
ł
Potencjał dwucząstkowy w funkcji odległości znika dla r = 1.8.
Faza ciekła i stała krzemu  stabilność
" Energia sieci w funkcji
" Średnia energia potencjalna
gęstości
w funkcji temperatury
Potencjał Stillingera-Webera dobrze odtwarza stabilność energetyczną Si.
Metody mechaniki kwantowej  ab initio 
metoda funkcjonału gęstość (DFT)
Procedura obliczeniowa:
" rozwiązanie równań mechaniki kwantowej w formalizmie
funkcjonału gęstości (DFT) albo w formalizmie funkcjonału
gęstości w modelu ciasnego wiązania (DFTB)
" obliczenie sił wynikających z rozwiązań kwantowo-mechanicznych
 twierdzenie Hellmana-Feynmana
" wykonanie całkowania równań ruchu
" obliczenie wielkości uśrednionych
Używamy formalizmu kanonicznego Hamiltona  współrzędnych i
pędów kanonicznych  kwantyzacja: przyporządkowanie
zmiennym kanonicznym operatorów położenia i pędu.
Hamiltonian układu - przypadek bezspinowy (najbardziej
uproszczony!)
Hamiltonian układu wielu cząstek H  jest operatorem działającym na
współrzędne elektronowe r oraz jonowe R. Dla układów atomów, o
niezbyt wysokich liczbach atomowych można zaniedbać efekty
relatywistyczne. Wówczas w najprostszym przypadku, hamiltonian
układu można przedstawić w następującej postaci :
r
)
r
h2 h2
H(R, r) = - "R - "r
" "
ą i
2Mą 2mi
ą i
ZąZ�e2
Ząe2 e2
+ +
r r -
r
r r
" " r "
ri - rj
Rą - R� ri - Rą
ą<� i,ą i< j
Funkcja falowa układu, zależy od współrzędnych jąder  R, oraz
elektronów  r
r
� = �(r,Rą) i =1..N ą =1..N'
ri
Przybliżenie adiabatyczne i Borna-Oppenheimera: położenia jąder są
traktowane jako parametr (ruch jąder: mechanika klasyczna).
Energia układu
E(R) = T + EN-e + Ee-e + EN-N
" Energia kinetyczna: T
h2 r h2
�(r r r ) - "r �(r r ,..rN ) Ćm(r ) "r Ćm(r )
r1, r2,..rN r1, r2 ri ri
"
i i
2m 2m
i
T = = -
r
"
�(r r ,..rN ) �(r r r ) Ćm(r ) Ćm(r )
r1, r2 r1, r2,..rN ri ri
i
" Energia oddziaływania jąder: EN-N
ZąZ�e2
EN-N =
r r
"
Rą - R�
ą`"�
Energia układu - cd
E(R) = T + EN-e + Ee-e + EN-N
" Energia oddziaływania jądro - elektron : EN-e
Ząe2 Ząe2
�(r r r ) - �(r r r ) Ćm(r ) Ćm(r )
r1, r2,..rN r1, r2,..rN ri ri
"
Rą - ri Rą - ri
i,ą
EN-e = = -
"
�(r r r ) �(r r r ) Ćm(r ) Ćm(r )
r1, r2,..rN r1, r2,..rN ri ri
i,ą
" Energia oddziaływania elektron - elektron
e2
�(r r r �(r r r
r1, r2,..rN) r1, r2,..rN)
"
ri - ri
i`" j
Ee-e = = 2J - K
�(r r r �(r r r
r1, r2,..rN) r1, r2,..rN)
gdzie J oraz K odpowiadają energii odpychania ładunków oraz energii
korelacji i wymiany.
Równanie Kohna-Shama
" Energię układu, znajdującego się w polu zewnętrznym Vext(r) można
przedstawić jako funkcjonał gęstości E[�], w stanie równowagi jest
określona przez warunek minimum, na przestrzeni funkcji falowych
spełniających warunek normalizacji, tzn.:
�E[�]
= �j
�Ć*
i
gdzie �j  jest mnożnikiem Lagrange a wynikającym z zasady normalizacji:
Ći Ćj = �ij
" Otrzymujemy szereg sprzężonych, nieliniowych równań na
funkcję Dj
ńł �ł
h2
" + Ve-e - Ve-N + �xc + Vext żłĆi = �iĆi
�ł-
2m
ół �ł
�j
Równanie Kohna-Shama w bazie
" Równania Kohna-Shama są więc przekształcone na formalne
równania na współczynniki Cij:
"H (C)Cjk = �iSijCjk
ij
j
gdzie:
h2
Hij = �i - " + Ve-e - Ve-N + �xc + Vext �j
Sij = �i �j
2m
Ponieważ Hij zależą od C- należy rozwiązywać równania iteracyjnie.
Najprostsza metoda kolejnych podstawień (SS  succesive substitutions)
polega na linearyzacji poprzez podstawienie C z poprzednich kroków.
Procedura iteracyjnego rozwiązywania równań
I. Wybierz wyjściowy zespół współczynników Cij
II. Utwórz wyjściowy zbiór orbitali molekularnych
III. Oblicz rozkład gęstości
IV. Oblicz nieliniowe wyrazy w hamiltonianie
V. Utwórz Hij
VI. Rozwiąż równanie Kohna-Shama na współczynniki Cij
VII. Oblicz nowy zbiór orbitali molekularnych
VIII. Oblicz rozkład gęstości
IX. Jeżeli rozkład gęstości nie zmienił się poza pewna granice  koniec
iteracji
X. Jeżeli rozkład gęstości się zmienił  wracamy do punktu IV
Warunek zbieżności zależy od fizycznych własności układu. Dla układów
o silnie niejednorodnych rozkładach rozkład może być obliczony przy
mniejszej liczbie iteracji, natomiast dla układów metali należy
wykonywać znacznie więcej iteracji.
Oddziaływanie N2 z Ga(l)
N2
Energia bariery na rozpad
6
N2 molecule
horizontally and
cluster of 19 atoms
5.8 eV
In
Ga
4
Al
4.8 eV
2
3.2 eV
0
1 2 3 4
Distance from surface(A)
Ga
Energia dysocjacji swobodnej S. Krukowski and Z. Romanowski
Obliczenia kwantowe, Dmol, DFT
cząsteczki N2 -- 9.8 eV/cząsteczkę
Excess energy [eV]
Dysocjacja N2 na powierzchni Ga
Dysocjacja N2 na powierzchni Ga
3,5
1.0A
3,0
2,5
2,0
1.6A
1,5
2.6A
1,0
0,8 1,2 1,6 2,0 2,4 2,8 3,2 3,6 4,0
d [A]
S. Krukowski and Z. Romanowski
QM DFT
N - N spacing [A]
Wzrost azotku galu na powierzchni GaN (0001): PA MBE
Struktura powierzchni GaN (0001)
Gęstość ładunku elektronowego dla
- diagram fazowy
atomu N na powierzchni GaN (0001)
Powierzchnia energii dla atomu N na
powierzchni GaN (0001) pokrytej In
Energia bariery na skok atomu N:
- powierzchnia czysta  1.3 eV
- powierzchnia pokryta In  05 eV
J. Neugebauer, T. Zywietz, M. Scheffler, J.E. Northrup,
H. Chen, R.M. Feenstra, PRL 90 (2003) 056101
Podziękowania dla ICM UW
Podziękowania dla ICM UW
" Za możliwość korzystania z oprogramowania
komercyjnego Fidap (ANSYS Inc.) oraz Abaqus
(Dassault SystŁmes)
" Wykorzystanie komputerów w ramach grantu
obliczeniowego G15-9