Transcript Szimulacio

4 ev Mernoki Informatika
BBTE Matematika es Informatika
Tudomanyos kutatasban ket terulet letezett kb
az 1600as evektol (Galilei kiserletei): ez a ket
terulet a kiserlet es az elmelet
Az elmelet matematikai leirast ad a
jelensegekre – kepleteket ir fel amelyek
osszefuggeseket allapitanak meg – ennek
prediktiv erteke van
A kiserlet mereseket vegez, ugy, hogy kozben
egy vagy tobb parametert ellenorzott modon
megvaltoztat – es ugy hogy egy vagy tobb
kulso tenyezo hatasat kikuszoboli vagy
minimizalja
Az elmelet es a kiserlet szoros kapcsolatban
allnak egymassal: az elmelet altal elore
kiszamolt hatasokat ellenorizni lehet
A kiserlet soran eszlelt uj jelensegekre
magyarazatot kell keresni elmeletileg
Az elmelet igy egyre jobb matematikai
modelleket gyart amelyek egyre jobban leirjak
azt ami a “valosagban” azaz a kiserletekben
tortennek
Az elso kiserletet Galileo Galilei vegezte a pisai
ferde toronybol (1589): kihajitott egy fem es
egy fa golyot amelyek azonos meretuek voltak
Arisztotelesz szerint a nehezebb targy
hamarabb kellene leessen (tehat a vasgolyo
hamarabb kellene leessen mint a fagolyo),
viszont egyszerre erkeztek meg
Arisztotelesz a surlodas miatt gondolta azt
amit gondolt: mindenhol jelen van a surlodas
korulottunk. Viszont az egyszeru eset az
hogyha nincs surlodas, ebben az esetben egy
test mozogni fog amig egy ero nem hat ra hogy
megallitsa. Ehhez kell hozzaadni a surlodast
hogy leirjuk azt az esetet amikor van surlodas.
Kiserlet > atirta az elmeletet
A Galilei kiserletet el lehet vegezni egy
vakumcsoben is – ez az amikor kiiktatjuk az
egyik tenyezot a kiserletbol (a levegovel valo
surlodast).
A matematikai leiras a leheto legegyszerubb
kell legyen (Occam borotvaja) es csak akkor
teszunk hozza ha szukseges hozzatenni ahhoz
hogy leirjunk valamit
Szuperpozicio elve: kulon-kulon le lehet irni a
hatasokat es ezek egyuttes hatasa a hatasok
osszege
f(x1+x2) = f(x1) + f(x2)
Vagyis a rendszerunk linearis es felbonthato
komponensekre es megertheto a
komponensek megertese alapjan
Pl. egy test eseten ahhoz hogy megertsuk a
viselkedeset, felbontjuk egyszeru
komponensekre, megertjuk ezen egyszeru
komponensek viselkedeset es osszegezzuk ezt
a viselkedest hogy megertsuk a nagymeretu
(makroszkopikus) test viselkedeset
Mikroszkopikus leiras > Makroszkopikus
viselkedes
Mikroszkopikus leiras > Makroszkopikus
viselkedes
Termodinamika: leirunk egy atomot vagy
molekulat, ebbol a leirasbol kiszamolhatjuk azt
hogy nagyon sok (10^26) atom vagy molekula
hogyan viselkedik
Mikor tudjuk ezt megcsinalni matmatikailag?
Ha az atomok fuggetlenek egymastol es eleg
egyetlen atomra megoldani a problemat
(peldaul idealis gaz modellje eseten)
Ha fuggetlenek, es nem viselkedik minden
atom ugyanugy, de ismerjuk az eloszlast
> Statisztikus Fizika
Pelda: idealis gaz modellje
Szukseg van arra hogy leirjam egyetlen
molekula mozgasat: egy tomeggel rendelkezo
kis pontnak veszem aminek a
mozgasegyenletet a Newton torvenye adja
meg f=ma. Idealis gaz eseten eltekintek a
gravitaciotol es a molekulak kozotti
kolcsonhatastol: csak az utkozesek szamitanak
Idealis gaz modellje
Az utkozesek kozul is csak a falakkal valo
utkozeseket tekintem. Ebben az esetben a
reszecskek impulzusa megvaltozik
p=mv az impulzus
Δp=2mvx az impulzusvaltozas a falnal
Δt = 2L/vx idonkent erkezik meg ujbol ide
F=Δp/Δt az ero amivel hat a falra
F=2mvx2/2L = mvx2/L
p=F/A = F/L2 a nyomas ezen a falon
p = mvx2/L3
Most terek at arra hogy hogyan osszegzem ezt
minden molekulara: atlagban mekkora a vx
sebesseguk ezt kell figyelembe vegyem
p = N m<vx2>/L3
(L3=V)
Ahol a <vx>-el az atlagos sebesseget jelolom
Ugyanolyan a sebesseg x,y,z iranyba
<vx2> = 1/3 <v2>
pV = νNA m 1/3 <v2> ahol hasznalom az
Avogadro-szamot (molekulak szama/mol)
pV = νNA m 1/3 <v2>
Ossze kell hozni az atlagsebesseget egy
makroszkopikus mennyiseggel – ehhez kellene
tudjam hogy oszlik el a molekulak mozgasi
energiaja (statisztikus fizika)
m <v2>/2 = 3/2kT minden szabadsagi fokra
ugyanannyi energia jut ami aranyos a
homerseklettel
pV = νNA k T ami mar csak makroszkopikus
pV = νNA k T ; NA k=R
pV = νRT altalanos gaztorveny, nagyon hasznos
mindenfele folyamat leirasara (izochor, izobar,
izoterm, politrop stb allapotvaltozasokra), le
lehet vele irni jol a hoerogepek mukodeset stb.
Az eredeti modell egyszerusegebol kovetkezik
hogy nem mindent tud leirni pl. amikor
cseppfolyositjak a gazakat
Ezt a modellt meg ki lehet boviteni arra hogy
leirja a cseppfolyositast, ekkor figyelembe kell
venni
 Azt hogy a molekulaknak van sajat terfogata
 Azt hogy van kozottuk egy vonzoero
Ha az atomok nagyon szorosan fuggnek
egymastol (pl egy kristalyracsban vannak) de ez
a kolcsonhatas ugyanugy nez ki mindenhol es a
szerkezet periodikus (ekkor ugyanis szinten
leirhatok egyetlen atomot, illetve az
eloszlasokat ismerve ki tudom szamolni a
rendszer viselkedeset)
Elmeleti leiras <-> kiserleti meresek
Mikor lehet gond az hogy kiszamoljam a
makroszkopikus viselkedest a mikroszkopikus
leiras alapjan?
Ha nem fuggetlenek egymastol az
atomok/molekulak: a kolcsonhatas nem
elhanyagolhato viszont nem is annyira eros
mint egy szilard krisztalyban ahol ismet
ugyanolyan (periodikusan)
Folyadekokra nem tudtak egy tesztelheto
modellt kesziteni.
Az elso konyv amit szamitogepes szimulaciorol
irtak, az Allen and Tildesley "Computer
Simulation of Liquids”
Az egyenleteket nem analitikusan hanem
numerikusan oldjuk meg.
Mit jelent az egyenletek numerikus megoldasa?
Azt jelenti hogy kiszamoljuk a reszecskek
palyajat az ido fuggvenyeben, tehat minden
reszecskere megkapjuk az x(t) y(t) es z(t)
egyenleteket. Mivel numerikus a megoldas,
nem magat a fuggvenyt kapjuk meg hanem a
fuggveny ertekeit egy sorozatban
x(t0), x(t0+dt), x(t0+2dt), x(t0+3dt), …
Ahhoz hogy ezt a megoldast megkapjuk,
szuksegunk van az egyetlen reszecsket leiro
mozgasegyenletre (equation of motion).
Esetunkben ez a mozgasegyenlet a klasszikus
Newton 2. torvenye lesz:
f = ma ahol a=d2x/dt2
d2x/dt2 = 1/m f ezt a differencialegyenletet kell
megoldani minden reszecskere
d2x1/dt2 = 1/m f1(x1,x2,x3,…)
d2x2/dt2 = 1/m f2(x1,x2,x3,…)
…
ez a csatolt differencialegyenlet-rendszerunk
van.Ez azert csatolt mert az erok
kiszamolasaban nem eleg tudni csak az illeto
reszecske poziciojat, hanem tudni kell minden
reszecske poziciojat – tehat az egyenleteket
egyszerre kell megoldani
1.
2.
3.
lepes: megadni minden reszecske poziciojat
es sebesseget (hatarfeltetelek a
differencialegyenlethez)
lepes: kiszamolni minden erot ami a
reszecskekre hat: ezt megtehetjuk mert
ismerjuk az osszes poziciot es sebesseget
egy kis lepest megtenni minden
reszecskevel – egy rovid dt ideig ugy
tekintjuk hogy az erok allandok
A 3. lepes utan ujbol megvannak a reszecskek
pozicioi es sebessegei, igy a 2. lepest ujra
megtehetjuk. A 2. es a 3.as lepest egymas utan
nagyon sokszor vegrehajtjuk (millioszor is) es
igy fokozatosan megkapjuk a reszecskek
palyajat, az x,y,z szamsorozatokat
Szimulaciokban nem ritka eset az hogy a dt
pikoszekundum nagysagrendu mig a teljes ido
nano vagy mikroszekundum nagysagrendu
Mi ez a kolcsonhatasi ero amit kiszamolunk a
reszecskek kozott? f1(x1,x2,x3,…)-el jeloltuk
ez a rendszerben a reszecske es minden mas
reszecske helyzetetol fugg igy altalanosan
felirva. Fel szoktak bontani kulonbozo
osszegekre ezt az erot, egy reszecske, ket
reszecske, harom reszecske stb. kolcsonhatasra
f1(x1,x2,x3,…) = f1(x1) + Σ f1(x1,xi) + Σ f1(x1,xi,xj)
+ Σ f1(x1,xi,xj,xk)+ …
az a tag ahol csak a reszecske koordinataja
szerepel, az a kulso eroket jelenti, ezek
fugghetnek attol hogy a reszecske hol van a
rendszerben (pl aramlas egy csoben: azt a
reszecsket amely kozepen van nagyobb erovel
huzza az aramlas mint ami a fal kozeleben van)
Σ f1(x1,xi)
ez a paronkenti kolcsonhatas (pairwise
interaction) – ennek van a legnagyobb szerepe
a rendszerben es ennek a szamolasa a legtobb
munka a molekularis dinamika szimulacio
futasa soran
Σ f1(x1,xi,xj) ez a harom-test kolcsonhatas, ezt
sokszor elhanyagoljak de vannak esetek amikor
fontos (pl elasztikus gombok vagy diszkek
eseten ez a tag nagyon fontos)
Σ f1(x1,xi)
hogy hat kolcson ket atom (molekula)?
Vegyuk a semleges atomok estetet. Ebben az
esetben a kolcsonhatast jol leirja a Van der
Waals kolcsonhatas: nagyon tavol egymastol az
atomok nem hatnak kolcson, egy adott
tavolsagon belul vonzzak egymast, mig hogyha
tul kozel kerulnek egymashoz akkor taszitjak
egymast
Ahhoz hogy ezt a kolcsonhatast (illetve
barmilyen atomok kozotti kolcsonhatast)
pontosan kiszamoljuk, meg kellene oldani az
atommagokra es az elektronokra is az ilelto
mozgasegyenleteket. Ezek a
mozgasegyenletek mar nem a klasszikus
egyenletek, hanem a Schrödinger- egyenletet
kellene megoldani: Hψ=Eψ. Ezekbol az
egyenletekbol megkapjuk a ψ hullamfuggvenyt
A szamitasokat eleg elvegezni csak az
elektronokra, mivel az atommagok tomege
sokkal nagyobb az elektronok tomegenel, ezert
mozgasuk is sokkal lassubb es igy kulon lehet
valasztani a ket mozgast: az elektronok gyors
es az atommagok lassu mozgasat (BornOppenheimer kozelites). Igy a megoldasban
elsosorban az elektronokra kiszamolt
hullamfuggvenyek erdekelnek minket.
A ψ hullamfuggveny magnitudoja (ψψ*) adja
meg az elektronok valoszinusegsuruseget a
terben. Ezekkel lehet pontosan kiszamolni azt
hogy milyen kolcsonhatas van kulonbozo
atomok kozott. Ha a szamolasokban ezt mindig
elvegezzuk azt az alapoktol indulo rendszernek
nevezik vagyis ab initio szamitasnak es ab initio
molekularis dinamikanak. Ez nagyon
szamitasigenyes, viszont pontos megoldas.
A gyakorlatban sokszor ezeket a szamitasokat
elvegzik egyszer es az igy kapott kolcsonhatast
hasznaljak fel a tovabbiakban anelkul hogy ujra
kiszamolnak az elektronokra a megoldasokat.
Ez a kvazi-klasszikus kozelites.
Ezek a kolcsonhatasi potencialok ki vannak
szamolva nagyon sok esetre elore es csak
felhasznaljak oket (AMBER, CHARMM, stb)
Nezzunk meg egy ilyen potencialt az egyik elso
potencialt amit hasznaltak: Lennard-Jones
vagy 6-12-es potencial
Ez ket semleges atom kolcsonhatasat irja le.
Ezek kozott Van der Waals erok hatnak.
Nagyon tavol levo atomok nem hatnak kolcson.
Noha az atomok semlegesek, az elektronfelho
kozpontja es az atommag kozpontja nem
mindig esik egybe, ilyenkor letrejon egy
elektromos dipolus. Ez a dipolus fluktual es
maga korul egy elektromos teret hoz letre. Ha
ezek a fluktualo dipolusok eleg kozel vannak
egymashoz akkor befolyasoljak egymast ugy,
hogy vonzoero jon letre a ket atom (vagy
molekula) kozott.
Ha az atomok nagyon kozel kerulnek
egymashoz akkor az elektronfelhok
atfedodnek es olyankor mar taszitas lep fel.
A kolcsonhatasi erot le lehet vezetni a
kolcsonhatasi potencialbol f(r)= - dU(r)/dr
A LJ potencialt azert nevezik 6-12es
potencialnak mert egy 1/r6 es egy 1/r12 tag is
van benne
Lennard-Jones (LJ) potencial:
r6 a vonzo
r12 a taszito ero
A LJ potencialt inkabb mar csak bemutato cellal
hasznaljak, az anyagok szimulalasara sokkal
jobb potencialok leteznek. Pl az embedded
atom potencial, kulonbozo elore kiszamolt
potencialok, polarizalhato potencialok, stb.
Bemutatunk meg egy egyszeru potencialt amit
az elso laborfeladatban fogunk hasznalni, ez a
toltott reszecskekre vonatkozik.
Coulomb potencial:
f = 1/(4πε0) q1q2/r2
amennyiben kolloidokrol beszelunk egy polaros
folyadekban, a folyadek arnyekolast vezet be
es ez az arnyekolas miatt a kolcsonhatas
hamarabb esik le (exponencialis tag):
f = 1/(4πε0) q1q2/r2 e-ar
Egyes potencialokat levaghatunk, ez azt jelenti
hogy egy adott rcutoff tavolsagon tul 0-nak
tekintjuk az eroket: igy egy csomo szamolast
megsporolunk.
Nem minden potencial vaghato le: 2
dimenzioban a LJ es a screened Coulomb
levaghato, a Coulomb nem vaghato le.
Mekkora rendszert tudunk szimulálni?
Egy millió atomot? Száz millió atomot?
Ezek mind nagyon kis méretek ahhoz képest
hogy mekkora egy makroszkopikus rendszer. A
makroszkopikus rendszer ugyanis ~1026
részecsét tartalmaz egyetlen kmolban
Nem vagyunk képesek ennyi részecskét
szimulálni: valamit ki kell találjunk
Ha a rendszer belsejében lezajló folyamatokra
vagyunk kíváncsiak (in the bulk), akkor úgy
tekinthetjük, hogy a rendszer felosztható sok
kis részre amelyek nagy többségben majdnem
ugyanúgy viselkednek (a felület nem, de nem
az érdekel most). Ezt a majdnem ugyanúgy-ot
kicseréljük arra hogy pontosan ugyanúgy. Ezzel
azt mondjuk hogy elég egy kis darabot
szimulálni amit saját maga másolatai vesznek
körül. Ez a periodikus határfeltétel (PBC)
PBC: Periodic Boundary Conditions
csak egyetlen kis
részt szimulálunk
–a teljes rendszert
leírjuk ezzel
PBC: Periodic Boundary Conditions
Két dolgot jelent ez a szimulációnkban:
 amikor egy részecske mozgása során kikerül
a dobozból, akkor az ellenkező oldalon bejön
a dobozba (PBC fold back).
 amikor egy részecske szomszédait nézzük, az
összes szomszédos dobozzal is számolni kell
Az összes szomszédos dobozzal levő
kölcsönhatás számolása:
 ha levágható a potenciál, akkor a dobozt
olyan nagynak vesszük hogy csak egyetlen
tükörképpel hasson kölcsön a részecske
 ha nem vágható le a potenciál akkor egy
végtelen sorozatot kell kiszámolni - helyes
sorrendberakással, csoportosítással ez
konvergál (Ewald summation)