Modelowanie procesw wzrostu krysztaw w skali atomowej

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: 88 80 244

e-mail: stach@unipress.waw.pl, mike@unipress.waw.pl

Zbigniew Żytkiewicz Instytut Fizyki PAN

02-668 Warszawa, Al. Lotników 32/46 E-mail: zytkie@ifpan.edu.pl

Wykład – 2 godz./tydzień – wtorek 9.15 – 11.00 Interdyscyplinarne Centrum Modelowania UW

Budynek Wydziału Geologii UW – sala 3075 http://www.icm.edu.pl/web/guest/edukacja

Modelowanie procesów wzrostu – obraz mikro

• Modelowanie - metody

• Modelowanie - zagadnienia i procesy fizyczne

Stanisław Krukowski

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

Metoda Monte Carlo – określenie przestrzeni zdarzeń

Aksjomatyczna definicja prawdopodobieństwa

• p(∅∅∅∅) = 0, p(E) = 1 • 0 ≤≤≤≤ p(A) ≤≤≤ 1≤

• p(A∪∪B) = p(A) + p(B)∪∪ zdarzenia się wykluczają tzn. gdy A ∩∩∩∩ B = ∅∅∅∅

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:

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ą

{ }

A

R

:

X

⇒

• Jednocześnie na algebrze zdarzeń losowych {A} określone jest prawdopodobieństwo, tzn. funkcja rzeczywista spełniająca warunki definicji aksjomatycznej.

{ }

A

[ ]

0

,

1

:

p

⇒

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

( )

x

0

P

≥

• unormowany

( )

x

1

P

i i

=

∑

• wyznaczany z pomiarów – np. rozkład energii o szerokości ∆E

( )

N

n

E

P

_k

=

k

gdzie n_k – liczba cząstek o energii w przedziale E_k,E_k+∆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.

p

( )

ξ

_i

= f( ) d

∫

ξ

ξ

f ( )

ξ

ξ

ξ

=

[0,1]

0

[0,1]

1

∈

∉







gdzie

Należy znaleźć procedurę generacji liczb losowych o rozkładzie prawdopodobieństwa danym funkcją rozkładu (próbkowanie)

(

)

p x

_k

=

f(x) dx

x

∈

∈

∫

Ω

Ω

Warunek zgodności - warunek równości prawdopodobieństw:

p x

(

>

y

)

= p( > )

ξ

ξ

_o

∫

∫

ξ ∞ −

ξ

0 x

d

=

dy

)

y

(

f

ξ

= f(y) dy

-x ∞

∫

• Rozkład Lorentza f x( ) = 1 1 1 + x2 π 1 π ξ 1 1 + y2 dy = -x ∞

∫

x = tan π 2 - 1

(

ξ

)

2      

Techniki odrzucania

x f(x)

0 1

Załóżmy że mamy funkcję rozkładu jednej zmiennej określoną na

przedziale [0,1] za pomocą skomplikowanej zależności:

Próbkowanie tego rozkładu za pomocą technik odrzucania:

• generacja liczby losowej ξξξξ₁ z przedziału [0,1]

• generacja liczby losowej ξξξξ₂ z przedziału [0,1]

• sprawdzamy warunek

• jeżeli warunek nie jest spełniony, wykonujemy procedurę generacji liczby od początku

ξ

₂

< f(

ξ

₁

) ⇒

x =

ξ

₁

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:

(

)

(

)

f p q p q kT i i i i , = 1 , Z exp - E      

( )

( )

      kT q V -exp Z 1 = q f _i i 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

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

Q =

=

f(q)

f(q' )

= exp

-E(q) - E(q' )

k T

p q q

p q q

( | ' )

( '| )













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

Funkcja korelacji gęstości – faza stała

Niższa gęstość

Metoda dynamiki molekularnej - numeryczna

analiza zachowania układów dynamicznych

&

p

_i

= -

H

q

_i

∂

∂

Formalizm Hamiltona _q&

i = H p_i ∂ ∂

{ }

(

_i

)

i i 2 i _{+ V q} 2m p = ) q , p ( H

∑

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ą v_i = p_i/m:

gdzie i i

=

v

q&

m v

_i

&

_i

_{= F}

_i

(

{ }

)

i i i V q q -= F ∂ ∂

Można sformułować w postaci jednego równania:

{ }

( )

( )

{ }

i i i i

m

q

F

q

f

q

&

=

=

&

Obliczanie potencjału sił – mechanika kwantowa

Problem – rozwiązanie równania Schrödingera na funkcję falową ΨΨΨ(R,r):Ψ

(

)

r) (R, H = t r , R i Ψ ∂ Ψ ∂ h

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

f

r)

(R,

=

ϕ

_R

Ψ

( )

=

H

(

r;

R

)

(r)

t

r

i

R

ϕ

_R

∂

∂ϕ

h

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

3 1 3 1 B N V v T mk 2 h       = << Λ = π

Potencjał Lennarda-Jonesa: gazy szlachetne

• Potencjał Lennarda-Jonesa

( )

                   σ −       σ ε = − − c c 6 ij 12 ij ij r > q 0 r < q q q 4 q V r 0,8 1,0 1,2 1,4 1,6 1,8 2,0 2,2 2,4 2,6 2,8 3,0 -1,00 -0,75 -0,50 -0,25 0,00 0,25 0,50 0,75 1,00 V /ε r/σ

Anizotropowy potencjał azotu N

₂ X Y Z ϑ_A ϕ_A Θ Φ ϑ B ϕ B R 0.30 0 .35 0.40 0.45 0.50 0.55 0. 60 -2.50 E-021 -2.00 E-021 -1.50 E-021 -1.00 E-021 -5.00 E-022 0.0 0E+000 E [ J ] R [nm] 0 0 0 0 0 0 0 0 90 0 0 0 0 0 90 0 9 0 0 0 0 90 90 9 0 0 0 0 45 0 4 5 0 0 0 45 0 13 5 0

Hartree-Fock approx. (van der Avoird et al. - 1984)

Anizotropowy potencjał azotu N

₂

– CCSD(T)

(P. Strąk - 2006) 0.40 0.45 0.50 0.55 0.60 0.65 0.00E+000 2.00E-021 4.00E-021 θ_A=0 θ_B=0 φ_B=0 E (J ) R (nm) 0.35 0.40 0.45 0.50 0.55 -2.00E-021 0.00E+000 2.00E-021 4.00E-021 θ A=90 θ B=0 φ B=0 E (J ) R(nm) 0.35 0.40 0.45 0.50 0.55 0.60 -2.00E-021 0.00E+000 2.00E-021 4.00E-021 θ_A=90 θ_B=90 φ B=0 E (J ) R(nm) 0.30 0.35 0.40 0.45 0.50 0.55 0.00E+000 5.00E-021 1.00E-020 1.50E-020 θ_A=90 θ_B=90 φ_B=90 E (J ) R(nm)

Obliczanie sił – najbardziej pracochłonny element metody

dynamiki molekularnej

Potencjał oddziaływania

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)

( )

q

V

( )

R

=

(r)

H

(

r;

R

)

(r)

V

=

ϕ

_R

ϕ

_R

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

Dwa typy metod całkowania równań ruchu:

1. Metody otwarte – metody predyktora

Oznacza ona że wartość x_n+1 jest wyznaczona wyłącznie przez wartości

w chwilach poprzednich tzn. x_n, x_n+1 , itd

2. Metody zamknięte – metody predyktor+korektor

Oznacza ona że wynik y_n+1 jest wyznaczony poprzez zależność od x_n,

x_n-1, tzn. jak poprzednio, lecz następnie wartość tę koryguje się

Metody otwarte: metoda Eulera

n n n 2 n n n 1 n

hf

2

1

+

h

v

+

x

h

x

2

1

+

h

x

+

x

=

x

₊

&

&

&

&

=

n n n n 1 n

v

h

v

v

hf

v

₊

=

+

&

=

+

Metoda Eulera polega na rozwinięciu w szereg Taylora względem czasu dla położeń q oraz prędkości v

Zachowanie metody Eulera jest bardzo złe – jest to metoda względnie mało stabilna

Metody zamknięte: zmodyfikowana metoda Eulera

Predyktor

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

y

_n+1

= x + hv +

1

h f

2 _n

2

n n

(

n

(

n 1

)

)

* n f f y 2 1 f = + ₊ Korektor

(

) ( )

[

n 1 n

]

1 + n * n 2 n n 1 n

f

y

f

y

4

1

y

=

f

h

2

1

+

hv

+

x

=

x

₊

+

₊

−

v

_{n 1+}

=

v

_n

+

hf

_n*

f

n+1

= f x

(

n+1

)

Metody otwarte: metoda Verleta

Metoda Verleta polega na użyciu rozwinięcia względem czasu w dwu krokach czasowych

Własności metody Verleta - prosty, szybki algorytm, dokładny do

rzędu O(h4_{) – powodują że jest to metoda bardzo często stosowana.}

gdzie x_n+1 = x_n + hx_n + 1h x2 _n + h x3 _n 2 1 6 ' '' ''' x_n−1 = x_n − hx_n + 1 h x2 _n − h x3 _n 2 1 6 ' '' '''

x

_n+1

=

2

x

_n

−

x

_n−1

+

h f x

2

(

_n

)

3 3 '' ' n 2 2 '' n ' n

t

x

x

t

x

x

t

x

x

∂

∂

≡

∂

∂

≡

∂

∂

≡

Po dodaniu stronami otrzymujemy niezależny od prędkości algorytm

Algorytm dla prędkości

(

)

v

x

x

h

n n n

=

+1

−

−1

2

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

i

T

h

=

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

i

T

h

=

Potencjały zależne od kierunku – hybrydyzacja sp

3

• Oddziaływania dwucząstkowe

Krzem – potencjał oddziaływania Stillingera-Webera:

( )

(

)

f r r _{r a} r q 2 0 = A Br exp 1 < a r > a -p − −         −

( )

v r_{2 ij} = rij      ε σ f₂

Założono że głębokość potencjału jest określona przez parametr, natomiast - jest wyznaczona przez żądanie aby f₂(21/6_{) = 0.}

• Oddziaływania trójcząstkowe

(

)

v₃ r r r_i, _j, _k = f₃ri , rj , rk      ε σ σ σ

(

) (

)

(

ji jk ijk

) (

ki kj ikj

)

jik ik ij k j i 3

,

r

,

r

h

,

r

,

r

h

,

r

,

r

h

r

,

r

,

r

f

θ

+

θ

+

θ

=

(

)

     > > < <       + θ         − γ + − γ λ = θ a r lub a r 0 a r i a r 3 1 cos a r a r exp , r , r h ik ij ik ij 2 jik ik ij jik ik ij

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 • γγγγ = 1.20 • λλλλ = 21.0

Faza ciekła i stała krzemu – stabilność

• Średnia energia potencjalna

w funkcji temperatury

• Energia sieci w funkcji gęstości

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

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 :

∑

∑

∑

∑

∑

< α _α α β < α _{α β} β α α α

−

+

−

−

−

+

∆

−

∆

−

=

α j i _{i j} 2 , i _i 2 2 i r i 2 R 2

r

r

e

R

r

e

Z

R

R

e

Z

Z

m

2

M

2

)

r

,

R

(

H

i

r

r

r

r

r

r

h

h

r

r

)

Funkcja falowa układu, zależy od współrzędnych jąder – R, oraz elektronów – r

(

r

_i

,

R

)

i

=

1

..

N

α

=

1

..

N

'

Ψ

=

Ψ

_α

r

r

Przybliżenie adiabatyczne i Borna-Oppenheimera: położenia jąder są traktowane jako parametr (ruch jąder: mechanika klasyczna).

Energia układu

• Energia kinetyczna: T

( )

R T E_{N e} E_{e e} E_{N N} E = + ₋ + ₋ + ₋

(

)

(

)

(

) (

)

( )

( )

( )

( )

∑

∑

φ φ φ ∆ φ − = Ψ Ψ Ψ ∆ − Ψ = i m i m i i m r 2 i m N 2 1 N 2 1 N 2 1 i r 2 N 2 1 r r r m 2 r r ,.. r , r r ,.. r , r r ,.. r , r m 2 r ,.. r , r T i i r r r h r r r r r r r r r r h r r r

• Energia oddziaływania jąder: E_N-N

∑

β ≠ α _{α β} β α −

−

=

R

R

e

Z

Z

E

2 N N

r

r

Energia układu - cd

• Energia oddziaływania jądro - elektron : E_N-e

( )

R T E_{N e} E_{e e} E_{N N}

E = + ₋ + ₋ + ₋

• Energia oddziaływania elektron - elektron

(

)

(

)

(

) (

)

( )

( )

( )

( )

∑

∑

α α α α α α − φ φ φ − φ − = Ψ Ψ Ψ − − Ψ = , i m i m i i m i 2 i m N 2 1 N 2 1 N 2 1 , i i 2 N 2 1 e N r r r r R e Z r r ,.. r , r r ,.. r , r r ,.. r , r r R e Z r ,.. r , r E r r r r r r r r r r r r r r r r

(

)

(

)

(

r ,r ,..r

) (

r ,r ,..r

)

2J K r ,.. r , r r r e r ,.. r , r E N 2 1 N 2 1 N 2 1 j i i i 2 N 2 1 e e = − Ψ Ψ Ψ − Ψ =

∑

≠ − r r r r r r r r r r r r

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 V_ext(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.:

[ ]

j * i

E

ε

=

δφ

ρ

δ

• Otrzymujemy szereg sprzężonych, nieliniowych równań na

funkcję ф_j

ij j

i

φ

=

δ

φ

gdzie ε_j – jest mnożnikiem Lagrange’a wynikającym z zasady normalizacji:

i i i ext xc N e e e 2

V

V

V

m

2

_



φ

=

ε

φ









+

ε

+

−

+

∆

−

h

_{− −}

Równanie Kohna-Shama w bazie

• Równania Kohna-Shama są więc przekształcone na formalne równania na współczynniki C_ij:

( )

_{ij jk} j i jk ij

C

C

S

C

H

∑

=

ε

gdzie:

Ponieważ H_ij 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.

j ext xc N e e e 2 i ij

V

V

V

m

2

H

=

χ

−

h

∆

+

₋

−

₋

+

ε

+

χ

j i ij

S

=

χ

χ

j

χ

Procedura iteracyjnego rozwiązywania równań

I. Wybierz wyjściowy zespół współczynników C_ij 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 H_ij

VI. Rozwiąż równanie Kohna-Shama na współczynniki C_ij 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

Oddziaływanie N

₂

z Ga(l)

S. Krukowski and Z. Romanowski Obliczenia kwantowe, Dmol, DFT

Ga

N

₂

Energia bariery na rozpad

1 2 3 4 0 2 4 6 3.2 eV 4.8 eV 5.8 eV N2 molecule horizontally and cluster of 19 atoms In Ga Al E xc e s s e n e rg y [ e V ]

Distance from surface(A)

Energia dysocjacji swobodnej cząsteczki N₂ -- 9.8 eV/cząsteczkę

Dysocjacja N

Dysocjacja N

₂₂

na powierzchni Ga

na powierzchni Ga

S. Krukowski and Z. Romanowski QM DFT 0,8 1,2 1,6 2,0 2,4 2,8 3,2 3,6 4,0 1,0 1,5 2,0 2,5 3,0 3,5 2.6A 1.6A 1.0A N N s p a c in g [ A ] d [A]

Wzrost azotku galu na powierzchni GaN (0001): PA MBE

J. Neugebauer, T. Zywietz, M. Scheffler, J.E. Northrup, H. Chen, R.M. Feenstra, PRL 90 (2003) 056101

Gęstość ładunku elektronowego dla atomu N na powierzchni GaN (0001) Struktura powierzchni GaN (0001)

- diagram fazowy

Powierzchnia energii dla atomu N na powierzchni GaN (0001) pokrytej In

Energia bariery na skok atomu N: - powierzchnia czysta – 1.3 eV

Podzi

Podzi

ę

_ę

kowania dla ICM UW

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