Klasyczna dynamika molekularna w analizie

Klasyczna dynamika molekularna w analizie

strukturalnej materiałów

Dynamika molekularna (MD) służy do obliczania trajektorii w przestrzeni fazowej zbioru molekuł, z których każda oddzielnie podlega klasycznym równaniom ruchu.

Dynamika molekularna jest techniką symulacji komputerowych, w której ewolucja czasowa układu oddziaływujących ze sobą atomów, jest wynikiem rozwiązania klasycznych równań ruchu. Dla każdego z N rozpatrywanych atomów spełniona jest druga zasada dynamiki Newtona.

Kinematyczne równanie ruchu to pewna zależność (bądź układ

zależności), określająca położenie ciała w przestrzeni w funkcji

czasu.

W inercjalnym układzie odniesienia jeśli siły działające na ciało nie równoważą się (czyli wypadkowa sił F jest różna od zera), to ciało porusza się z przyspieszeniem wprost proporcjonalnym do siły

wypadkowej, a odwrotnie proporcjonalnym do masy ciała.

II zasada dynamiki Newtona

Rozwijamy funkcje i w szereg Taylora wokół punktu t: ^{r }

 ^{t   t}  ^{r }

 ^{t   t} 

 ^{t   t}    ^{ r t   t V}   ^{t }

^{ t}

^a   ^{t  ...}

r 







 ^{t   t}    ^{ r t   t V}   ^{t }

^{ t}

^a   ^{t  ...}

r 







Dodając stronami:

  ²    

  ^O   ^t

m t t F

t t

r t

r t

t r

i i i

i

i

         

 





Odejmując stronami:

       

2 O t

t

t t

r t

t t r

V

 











   



Algorytm Verleta (1967)

1. Działanie algorytmu Verleta może się rozpocząć, jeśli dane są dwa pierwsze położenia cząstek R

i R

.

2. Błąd jest większy dla prędkości niż dla położeń.

3. Do wyznaczenia położeń prędkości nie są potrzebne.

4. Duża dokładność w określeniu położenia.

5. Stabilność.

6. Prędkość wyznaczana jest jeden krok po położeniu.

( Rozwiązaniem problemu jest: algorytm velocity Verlet )

 ^{t t}  ^V   ^{t t}  ^a    ^{t a t t}  

V  _i     _i    _i   _i  

2 1

Algorytm Verleta (1967)

Energia potencjalna a siła

Jeżeli znany jest rozkład przestrzenny energii potencjalnej danego układu ciał, to można wyznaczyć siłę działającą na to ciało (układ ciał) obliczając gradient energii potencjalnej.

Zachowane:

• liczba cząstek – (N)

• objętość pudła – (V)

• energia całkowita – (E)

W rzeczywistości:

• brak zachowania E

• badania w stałej temperaturze – T=const

• badania pod stałym ciśnieniem – p=const

Schemat NVE – zespół mikrokanoniczny

Ograniczenie stałej temperatury wprowadza się przeskalowywując prędkości w czasie trwania symulacji, aby zapewnić właściwą relację pomiędzy całkowitą energią kinetyczną a pożądaną temperaturą.

Takie przeskalowanie jest równoważne wprowadzeniu pewnej siły (analog tarcia):

i

i

p

f  







 > 0 – proces ochładzania

 < 0 – proces podgrzewania

Schemat NVT – zespół kanoniczny termostat Nose-Hoovera

 

 



 

 

 N

i

B i

i

V Gk T

M

m

1 



M – „masa” termostatu

G – całkowita liczba stopni swobody układu ( wzrasta z 3N-3 do 3N-2 )

Schemat NPT– zespół izobaryczno-izotermiczny ( makrokanoniczny ) – schemat Andersena

Schemat ten nie wynika z żadnych zasad czy pojęć fizycznych, ale jest jedynie pewną metodą modelowania wpływu otoczenia na realistyczne układy fizyczne.

Wprowadza się dodatkową energię:



 



 



  

2

2 P

dt d E M

M – masa efektywna związana z rozszerzaniem się i kurczeniem objętości

 pod wpływem ciśnienia P

przy założeniu, że układ jest w kontakcie

z otoczeniem utrzymującym stałe ciśnienie.

Konstrukcja pola sił – CaZrO

ij ij ij ij

ij

ij r r

C r r

0 j i ij 6

4 q exp q

A ) (

V  ^_^{ } 







  

Potencjał Buckinghama:

1. Struktura krystaliczna

Ca (0,0,0) Zr (0.5,0.5,0.5)

O (0.5,0.5,0) Pm3m a=4.020 Å

2 3 4 5 6 7 8 9

-0.5 0.0 0.5 1.0 1.5

V Ca-O [a.u.]

r [Å]

2. Stałe elastyczne

C₁₁ E

C₁₂ G

C₄₄ B

2 2

V V E B 

 

Krzywizna potencjału w okolicy minimum

Moduł ściśliwości:

Konstrukcja pola sił – CaZrO

Zastosowanie metod ab intio np. WIEN2k

1. Ładunek jonów

Jon formal oblicz.

Ca +2 +1.61

Zr +4 +2.56

O -2 -1.39

2. Stałe elastyczne

[GPa] [GPa]

c₁₁=360 E=270 C₁₂=76 B=171 C₄₄=87 G=109

Równanie stanu Birch–Murnaghan:

420 440 460 480 500

-9004.280 -9004.275 -9004.270 -9004.265 -9004.260 -9004.255 -9004.250 -9004.245 -9004.240

E [a.u.]

V [ų]

Konstrukcja pola sił – CaZrO

Jony A_ij [eV] ρ_ij [1/Å] Cij [A⁶]

Ca-O 7488.39 0.297 0 Zr-O 33211.0 0.189 0

O-O 31665.9 0.149 28.46

ij ij

ij ij

ij r r

C r r

0 j i ij 6

4 q exp q

A ) (

V  ^_^{ } 







  

Potencjał Buckinghama:

2 3 4 5 6 7 8 9

-1.0 -0.5 0.0 0.5 1.0 1.5

V [a.u.]

r [Å]

Ca-O

Zr-O

Konstrukcja pola sił – CaZrO

• ok. 11000 at. odpowiadające stechiometrii

• at. rozmieszczone tak

jak w komórce elementar.

•periodyczne warunki brzegowe

• def. współcz. potenc.

• krok czasowy: 2 fs(10^-15)

• temp. 300 K

npt

• T_p=300 K

• T_k=5000 K

• 20 000 it.

(1·10^-11 s )

npt

• T_p=5000 K

• T_k=5000 K

• 20 000 it.

(1·10^-11 s )

npt

• T_p=5000 K

• T_k=2300 K

• 20 000 it.

(1·10^-11 s ) npt

• T_p=2300 K

• T_k=2300 K

• 20 000 it.

(1·10^-11 s )

npt

• T_p=2300 K

• T_k=300 K

• 20 000 it.

(1·10^-11 s )

npt

• T_p=300 K

• T_k=300 K

• 20 000 it.

(1·10^-11 s )

Konstrukcja pola sił – CaZrO

Konstrukcja pola sił – CaZrO

Źle dobrany punkt obcięcia potencjału

2 3 4 5 6 7 8 9

-1.0 -0.5 0.0 0.5 1.0 1.5

V [a.u.]

r [Å]

Ca-O

Zr-O

Konstrukcja pola sił

ij ij

ij ij

ij r r

C r r

0 j i ij 6

4 q exp q

A ) (

V  ^_^{ } 







  

Symulacja szkła z układu P

O

-Fe

O

-Na

O

1. Propozycja składu szkła (1-x)(0.6P

O

-0.4Fe

O

)-xNa

O; x=0,0.1,0.2,…..

2. W trakcie syntezy 20-30 % jonów Fe

ulega redukcji do Fe

3. Nowy skład: (1-x)(0.55P

O

-0.30Fe

O

-0.15FeO)-xNa

O

Cel:

Określenie wpływu Na

O na właściwości strukturalne

szkła 60P

O

-40Fe

O

Symulacja szkła z układu P

O

-Fe

O

-Na

O

dla x=0 0.55P

O

0.30Fe

O

0.15FeO

110 at. P 60 at. Fe3+

15 at Fe2+

380 at O

Suma 565 atomów

Ułożenie przypadkowe lub nie

=2.91 g/cm

Z eksperymentu:

• Fe

/Fe

• gęstość

a=19.0 Å

periodyczne warunki brzegowe

Symulacja szkła z układu P

O

-Fe

O

-Na

O

Symulacja szkła z układu P

O

-Fe

O

-Na

O

timestep 1 fs = 10^-15 s

1 2 3 4 5 6 7 0

2 4 6 8 10 12 14 16 18 20 22 24

szkło 4 Cs - O

Al - O

B - O

Si - O

r [Å]

PDF [a.u.] ^{Na - O}

Ca - O

1 2 3 4 5 6 7

0.0 0.5 1.0 1.5 2.0 2.5 3.0

szkło 4

r [Å]

RDF [a.u.]

Symulacja szkła z układu P

O

-Fe

O

-Na

O

Symulacja szkła z układu P

O

-Fe

O

-Na

O

Symulacja szkła z układu P

O

-Fe

O

-Na

O

Symulacja szkła z układu P

O

-Fe

O

-Na

O

3.0x10⁶ 3.3x10⁶ 3.6x10⁶ 3.9x10⁶ 2

4 6 8 10 12 14 16 18

msd Na [Å2 ]

krok czasowy T = 1500 K

   

               















i

i i

n

i

i i

n

i

i

i

t x y t y z t z

n x MSD

1

2 1

2 1

2

0 0

1 0

t D MSD

 6

b t

a

MSD   

6 6

6

a t

b D a

t







1.5x10^-4 2.0x10^-4 2.5x10^-4 3.0x10^-4 3.5x10^-4 4.0x10^-4 -14

-12 -10 -8

ln DNa

1/T [1/K]



 



 

 RT

D E

D ₀ exp

Klastrowanie alkaliów (MAE)

Uszkodzenie radiacyjne

Ścinanie kruchego materiału

Pękanie

Krystalizacja - aluminium

Uderzenie kulki diamentowej na warstwę grafenu