A statisztikus fizikai szimulációk alapjai és a Monte Carlo módszer '12

Statisztikus fizikai szimulációk alapjai

Legyen egy sokaságunk, aminek az energiája $E=E(\vec{r})$, ahol $\vec{r}$ az összes szabadsái fokot tartalmazza (mechanikai rendszerre $\vec{r} = (\vec{q}, \vec{p})$). Ha a rendszer követi a Boltzmann-statisztikát, akkor a rendszer $[\vec{r}, \vec{r}+d \vec{r}]$ intervallumban való megtalálási valószínűsége:

$P(\vec{r})d\vec{r} = \frac{e^{-\beta E(\vec{r})}}{Z}d\vec{r}$

Egy A mennyiség átlagát a sokaságon a következőképp számoljuk (PS = phase space = fázistér):

$\langle A \rangle = \int_{PS} A(\vec{r}) \frac{e^{-\beta E(\vec{r})}}{Z}d\vec{r}$

Az integrált csak akkor tudjuk egzaktul kiszámolni, ha a fázistér minden pontjára kiértékeljük a függvényt, és úgy átlagolunk. A Monte-Carlo integrálással ezt a problémát oldhatjuk meg. Az $\int d\vec{r}$ helyére $\frac{1}{N} \sum$ írunk, így egy véges összeget kell elvégeznünk. A pontok számát (N) az határozza meg, hogy mekkora pontossággal akarjuk az integrált meghatározni. Az integrál dimenziójától függetlenül a hiba $1/\sqrt{N}$ szerint csökken, vagyis 10x nagyobb pontossághoz 100x annyi pontot kell kiszámolni.

Molekuladinamika

A molekuladinamika a részecskék mikroszkopikus dinamikájának követésével foglalkozik. Valódi rendszerben 10²³ nagyságrendű részecske van, ezt a mai számítógépekkel szimulálni lehetetlen. Azonban ennél kevesebb részecskét is elég szimulálnunk ahhoz, hogy a termodinamikai tulajdonságokat vizsgálhassuk. A szimulációkban a párkölcsönhatásokat vesszük figyelembe, ám lehetséges közelítéseket tenni. A részecskék közt a Van der Waals erő hat, ami elég gyorsan lecseng, így távoli részecskék közt elhanyagolható (ezen a pontot különbözik a molekuladinamika és a gravitációs N-test szimuláció, ahol az 1/r²-es erő miatt nem hanyagolható el a kölcsönhatás).

A szimulációban van tehát N darab részecske, melyekre a Newton-törvények alapján kiszámítjuk a rájuk ható erőt, majd léptetjük a rendszert. Mivel a kezdőállapotokat általában nem a termodinamikai egyensúlyból indítjuk, meg kell várni, hogy a rendszer beálljon abba. Hogy mikor állt be, azt az ekvipartíció segítségével mutathatjuk ki:

$\left \langle \frac{1}{2}mv^2 \right \rangle = \frac{3}{2}kT$

Ha a rendszer elérte az egyensúlyi állapotát, mérhetővé válnak a termodinamikai változók (hőmérséklet, nyomás, hőkapacitás, stb.)

Fizikai mennyiségek mérése a szimulációkban

Pillanatnyi hőmérséklet

Bár azt mondtuk, hogy a rendszert az egyensúlyi állapotában vizsgáljuk, a nem egyensúlyi állapotban is bevezethetünk egy pillanatnyi hőmérsékletnek nevezett mennyiséget. Termikus egyensúlyban igaz az ekvipartíció:

$3(N-1) \cdot \frac{1}{2}k_BT = \left \langle \frac{m}{2} \sum_{i=1}^N v_i^2 \right \rangle$,

ahol 3(N-1) a szabadsági fokok száma (a TKP 3 koordinátája van levonva). A nem egyensúlyi helyzetben, bár nem igaz az ekvipartíció, ezt a képletet alkalmazhatjuk a hőmérséklet meghatározására.

Teljes energia

$E = \frac{m}{2} \sum_{i=1}^N v_i^2 + \sum_{i \neq j} U(|r_i - r_j|)$

Hőkapacitás

Fluktuáció-disszipáció tétel alapján:

$C_V = \frac{1}{k_BT^2}[\langle E^2 \rangle - \langle E \rangle^2]$

Nyomás

A nyomást a viriál-tétel segítségével fejezhetjük ki:

$pV = Nk_BT + \frac{1}{3} \left \langle \sum_{i<j} r_{ij} \cdot F_{ij} \right \rangle$

Verlet-algoritmus

Molekuladinamikai szimulációkban legtöbbször a Verlet-algoritmust használják a diffegyenletek megoldására. A Runge-Kuttával szemben több előnye van:

• gyorsabb, mert egy lépésben csak egyszer kell a gyorsulásokat számolni
• majdnem olyan pontos, mint a RK4 ($O(\tau^4) \text{vs.} O(\tau^5)$)
• jól megőrzi az energiát
• időtükrözésre nem változik (ez a részleges egyensúly feltétele miatt fontos)

A Verlet-algoritmus egy lépése (R(t) a koordináták, V(t) a sebességek, A(t) a gyorsulások):

$R_{n+1} = 2R_n - R_{n-1} + \tau^2 A_n + O(\tau^4)\,$

$V_n = \frac{R_{n+1} - R_{n-1}}{2\tau} + O(\tau^2)$

Hátrányai:

• két előző lépést használ, így nem indítható 1 kezdeti feltételből
• a sebesség és a pozíció nem egyszerre frissítődik, a sebesség "le van maradva"

Megoldás: velocity-Verlet algoritmus:

$R_{n+1} = R_n + \tau V_n + \frac{\tau^2}{2} A_n + O(\tau^3)$

$V_{n+1} = V_n + \frac{\tau}{2} (A_{n+1} + A_n) + O(\tau^3)$

• R-ben már csak $O(\tau^3)$ pontosságú
• ez már nem 2 előző lépést használ
• koordináták és sebességek egyszerre frissülnek
• sebesség is $O(\tau^3)$ pontosságú

Gravitációs N-test szimuláció

Gravitációs N-test szimulációkor, a kölcsönhatásnak az említett $r^{2}$-es lecsengése miatt, nem lehet eltekinteni a hosszútávú kölcsönhatásoktól.

$F_{grav}=G\frac{m_{1}m_{2}}{r^{2}}$

A sok, nem elhanyagolható kölcsönhatás miatt a direkt-szimuláció hamar nehézségekbe ütközik. Ezért érdemes bevezetni olyan finomításokat, mint például a léptetési idő variálása: minden részecske saját időlépéssel rendelkezik, így a különböző dinamikájú részecskéket nem kell minden időpillanatban léptetni.

Emellett létezik két alapvető, elterjedt algoritmus, melyek különböző módon egyszerűsítik a rendszert:

• Fa-algoritmusok (Tree methods):

A Tree method lényege, hogy a teret köbös cellákra osztják (octree): első lépésben 8 kockára osztják a teret, majd a megadott mélységig a cellákat tovább osztják 8 részre (az első lépéshez hasonlóan). Ha a felbontás indexelése következetes, akkor a hierarchia miatt minden cella egyértelműen címezhető (pl. "138" az első cella 3. cellájának 8. cellája). A felbontás után csak a szomszédos/közeli cellák részecskéire számolják individuálisan a kölcsönhatást. A távolabbi cellákat egy részecskeként veszik figyelembe, melyet a cella össztömege reprezentál annak tömegközéppontjában.

(Hasonló tér-felosztást lehet elérni szimulációk során a KD-tree binári felbontással. A GADGET azonban pl az octree-t használja.)

• Részecske-háló (Particle mesh method - PM method):

A PM method-ban a teret diszkretizálják egy ráccsal, és kiszámolják a gravitációs-potenciált a Poisson-egyenlet segítségével a rácspontokra:

$\nabla^2\Phi=4\pi G{\rho},\,$

ahol G a Newton-konstans, $\rho$ pedig a sűrűség (a rácspontban lévő részecskék száma). Ez az egyenlet könnyen kezelhető FFT után a frekvencia térben:

$\hat{\Phi}=4\pi G\frac{\hat{\rho}}{k^2},$

ahol $\vec{k}$ az együtt-mozgó hullámszám-vektor (kalap jelöli a Fourier-transzformáltat). Ezek után $\vec{k}$-val szorozva és inverz Fourier-transzformálva megkapható a gravitációs-tér. Mivel a rács mérete korlátozza ezt a módszert, ezért a kis-skálás erők számításához általában más módszerekkel kombinálva alkalmazzák ezt az eljárást. (Néha adaptív rácsozást is alkalmaznak, ami nagyobb sűrűségnél nagyobb rács-felbontást jelent.)

(Ezen módszerek természetesen más N-test szimulációs területeket is használatosak lehetnek, mint pl gázok, folyadékok szimulációjakor...)

Meg lehet említeni a GADGET-2 szimulációs programot, melyben található egy hibrid módszer is: a Tree és PM kombinációja (TreePM)

A Metropolis algoritmus

A Metropolis algoritmus egy eljárás, amivel bonyolult valószínűségi eloszlások mintavételezhetőek. Másképp megfogalmazva szeretnénk előállítani egy $P$ eloszlást, a kérdés az, hogy milyen átmeneti valószínűségekkel $(W)$ vezérelt folyamat (egészen pontosan Markov-folyamat) vezet ilyenre? A módszer kihasználja a részletes egyensúly elvét, ez esetünkben azt jelenti, hogy:

$P_i W_{ij} = P_j W_{ji}\,$ teljesül minden $i, j$-re.

Ez még kevés feltétel az átmeneti valószínűségekre, ezért a Metropolisz megkötés a következő:

$W_{ij} = \mathrm{min}\left[1, \frac{P_j}{P_i} \right]$

Bár általánosan tetszőleges eloszlás előállítható, fizikában általában valamilyen egyensúlyi eloszlást szeretnénk előállítani. Például, ha termális egyensúlyt szeretnénk, akkor $P \propto \exp\left( -\beta E \right)$, és a Metropolisz megkötésben szereplő hányados is exponenciális alakban írható. A $P$-t előállító eljárás a következő:

1. Induljunk ki egy A konfigurációból.
2. Állítsuk elő a megváltoztatott B konfigurációt.
3. Számoljuk ki erre a két konfigurációra a megkötésben szereplő hányadost.
4. Generáljunk egy 0-1 intervallumon vett egyenletes eloszlású véletlen számot: X.
5. Ha $W_{AB} > X$ akkor a B konfigurációt fogadjuk el, ellenkező esetben marad az A.
6. Ismételjük 1-től.

A Monte-Carlo módszer

A Monte-Carlo módszernek nevezzük az olyan eljárásokat, amelyek a problémákat random számok és valószínűségek felhasználásával oldják meg. Az eljárás során ismétlődően kiértékelünk egy determinisztikus modellt, random számokat használva inputnak. Akkor használják, ha a feladat nagyon összetett, nemlineáris, bonyolultak a határfeltételei, vagy nagyon sok paramétertől függ, illetve sok dimenziós.

Használata:

1. Állítsuk föl a modellt: y = f(x₁, x₂, ..., x_q)
2. Generáljunk random számokat inputnak: x_i1, x_i2, ..., x_iq
3. Értékeljük ki a modellt, az eredményt tároljuk el y_i-ben
4. Ismételjük a 2. és 3. lépéseket n-szer
5. Elemezzük az eredményeket hisztogram, összesítő statisztikák, stb. segítségével

Klasszikus példa a $\pi$ kiszámítása:

Rajzoljunk egy r sugarú kört és egy köré írt négyzetet (2r oldalú). Ha véletlenszerűen mintavételezzük a négyzet területét, akkor a körben lévő pontok száma úgy aránylik a körben és négyzetben lévő pontok számához (tehát az összeshez), mint a síkidomok terülte. Így $\pi$ közelítőleg felírható egy hányadossal:

$\pi=4\frac{\#k\ddot{o}r\ddot{o}n\_bel\ddot{u}li\_pontok}{\#\ddot{o}sszes\_pont}=4\frac{r^{2}\pi}{(2r)^{2}}$

Itt is egyszerűen belátható, hogy a több mintavételezési pont pontosabb értéket eredményez.