Számítógépes szimulációk a statisztikus fizikában – gyakorlat

2011/12. őszi szemeszter

Tartalomjegyzék

Bevezető

   Kedves Hallgatók!

   A korábbi évek tapasztalatai alapján a Számítógépes szimulációk a statisztikus fizikában c. tárgyon a féléves nagy házi feladat nehézséget szokott okozni a hallgatóság egy részének. Ezek a nehézségek azt jelzik, hogy a hallgatóság C, C++, FORTRAN programozástudása és tapasztalata szerényebb a feltételezettnél. Elhatároztuk, hogy idén tartunk egy technikai jellegű gyakorlatot Kertész Tanár Úr előadása mellé.

   A gyakorlat célja felfrissíteni a programozástudást, példákon keresztül az órán szereplő szimulációs technikákat gyakorolni, élőben bemutatni, kipróbálni. Aki jár, annak könnyebb lesz megírni a nagyházit.

   A gyakorlat fakultatív. A vizsgajegybe a gyakorlaton nyújtott teljesítmény sem pozitív, sem negatív irányban nem számít bele.

   Ez részünkről is egy első próbálkozás. Bármilyen ötletet, segítséget, építő kritikát szívesen fogadunk. Siker esetén várható, hogy a következő években is tartozni fog gyakorlat a tárgyhoz.

    Üdvözlettel,
    a (leendő) Gyakorlatvezetők.

Általános információk

Cygwin howto


1. hét

0. házi feladatok: Szoftverek telepítése

„Beadási határidő:” a 2. heti gyakorlat.

A következő programokat telepítse a saját számítógépére (az összes ajánlott program ingyenes):

(Linux alól az összes program telepíthető csomagkezelőn keresztül.)
  1. (hf01.zip) Fordítsa le az itt található kódot és készítse el doxygen segítségével a kód dokumentációját html formátumban.
  2. (hf02.zip) Töltse be a mellékelt fájlt gnuplotba, a kapott ábrát és a tex fájlt felhasználva pedig készítsen LaTeX segítségével egy pdf-et.
  3. A saját számítógépén (amelyen dolgozni fog a gyakorlatokon – ha nem saját számítógépen fog dolgozni, ez a házi feladat tárgytalan) hozza működőképes állapotba a BME WLAN kapcsolatot! Infó: https://net.bme.hu/wlan/. MAC cím regisztráció BME-hallgatóknak: https://accadmin.hszk.bme.hu/

2. hét

1. gyakorlat (szeptember 12.): C-ismétlés

  1. Bevezető előadás: pres1.pdf
  2. C-frissítő

    Példaprogram Rövid leírás, jegyzet
    macro1.c A legegyszerűbb makró: érték elnevezése.

    Példaprogram: kiírjuk α és sin(α) értékét 30 fokonként. Fordítási parancssor:
    gcc macro1.c -lm -Wall -o macro1.x
    ahol -lm jelenti a matematikai könytárral való linkelést (ha van a programban #include <math.h> és használunk is belőle függvéynt (pl. sin), akkor kell); -Wall = Warning all; -o macro1.x: a leendő futtatható program neve.

    macro2.c Paraméteres makró

    • A makrók paraméterezhetők. Úgy használandó, mint egy függvény, de nem az! A makrót a preprocesszor dolgozza fel. MAKRÓ = SZÖVEGHELYETTESÍTÉS!
    • Pro: a makró olyan gyors, mintha a meghívás helyére gépeltük volna a kódját. (A prepocesszor ténylegesen ezt csinálja.) Azaz gyorsabb, mint a függvényhívás. Rövid, gyakran használt műveletsorokat érdemes makróként definiálni.
    • Kontra: ha ugyanazzal az értékkel hívjuk meg sokszor a műveletsort, akkor érdemesebb csak 1x kiszámítani és tárolni. Mivel a makró csak szöveghelyettesítés, nincs típusellenőrzés, borulhat a precedencia (ha figylmetlenek vagyunk) és a mellékhatások miatt vicces dolgokat művelhet.
    • A makróparamétereket és magát a makrót is mindig zárójelekkel védjük. Ellenkező esetben ez történik: macro2-BUG.c.
    • Mellékhatásgubancra példa: macro3-DEFECT.c. Íratlan szabály: makróparaméterként ne adjunk át mellékhatással rendelkező kifejezést! (De leginkább függvényparaméterként se!)

    dice.c Túlindexelős példaprogram Példa: írjunk egy programot, ami 10 darab (pszeudo-)dobókockadobás eredményét egy tömbben tárolja és ki is írja.

    Tanulságok:

    • A tömbméret megadásának preferált módja: makróval.
    • Ha a véletlenszám-generátort ugyanarról a számról (seed = mag) indítom, ugyanazt a véletlen sorozatot fogja előállítani. Hasznos lehet hibakeresésnél.
    • A véletlenszám-generátor indítható a rendszeridővel, de ez csak 1 s felbontású. Nem ajánlott.
    • Hogyan kell 1-től 6-ig terjedő véletlen egészet sorsolni? Éles feladatban a rand() nem ajánlott.
    • Gyakorlat: egy figyelmetlenség miatt hiba kelthető, ha a tömbméretet átírjuk 9-re. Mi történik? Miért? Mi volt a figyelmetlenség?

      Végtelen ciklus keletkezik. Ebben a sorban van a figyelmetlenség:
        for (i=0; i<10; i++) {
      Így hangzik helyesen:
        for (i=0; i<MAXDICE; i++) {
      Ha a tömbméret 9, (ekkor az érvényes indexek: 0, 1, ..., 8) akkor dice[9] már nem a tömbben van, hanem i-t jelenti. A ciklusváltozót írjuk át minden lépésben 1 és 6 közti számra, így nehezen fogja elérni a kilépési feltételt.

    functions-memory.c Példaprogram: számítsuk ki 1 000 000 (pszeudo-)kockadobás alapján a dobás várható értékét és szórását!

    Hosszú lenne a teljes magyarázatot begépelni... Frissítsétek fel a következő témaköröket:

    • saját függvény írása, deklaráció, definíció;
    • függvényparaméterek: bemeneti értékek, kimeneti változók, tömb, visszatérési érték;
    • dinamikus memóriafoglalás, felszabadítás;
    • hibakezelés.

    Linkajánló:

    pow.c Több forrásfájl használata; parancsosri paraméterek; I/O átirányítás

    Példaprogram: olvassunk a standard bemenetről két valós számot tartalmazó sorokat a fájlvége jelig. Pl.: x1 y1 \ x2 y2 \ ... A standard kimenetre írjuk ki az x1 y1e \ x2 y2e \ ... adatokat, ahol e értékét a parancsorból olvassuk be a -e kapcsoló után!

    Példa: (fájlvége jel begépelése: Ctrl+D)
    ./pow.x -e 2.0
    1 3
    1.000000 9.000000

    A program a hf01.zip-ben található CommandLineParsing.h-n és CommandLineParsing.c-n alapul. Tegyük fel, hogy ezeket a fájlokat a hf01 könytárba másoltuk. A használathoz a forrásfájlunkban el kell helyezni a következő sort:
    #include "hf01/CommandLineParsing.h"
    A fordítás parancssora:
    gcc pow.c hf01/CommandLineParsing.c -lm -Wall -o pow.x

    Figyeljétek meg, hogy oldottam meg a bemenet soronkénti beolvasását a fájlvége jelig.

    Linux alatt (és Windows alatt is, csak azt nem próbáltam most ki) így irányíthatók át a standard ki-, be- és hibamenetek:
    ./pow.x -e 2.0 < identity.dat > square.dat 2> error-messages-go-here.txt
    Egy program kimenete közvetlenül egy másik program bemenetére irányítható a | jellel. (cat utasítás: kiírja a fájl tartalmát) Pl.:
    cat identity.dat | ./pow.x -e 2.0 | ./pow.x -e 0.5

    A Gnuplot tudja ábrázolni egy program kimenetét. Ehhez az ábrázolandó „fájl” „nevét” egy < jellel kell kezdeni:
    gnuplot
    plot '< ./pow.x -e 2.5 < identity.dat'
    plot '< ./pow.x -e 2.5 < identity.dat' with line linewidth 5
    quit


3. hét

1. adag kis házi feladat: C-ismétlés

Megoldások

  1. (ötperces – ) sinonehalfpi-BUG.c Írtam egy rövid programot sin( 1/2 pi ) kiszámítására, de nem a várt eredményt adja. Miért nem? Hogyan kell helyesen megírni?
  2. (ötperces – makró) macro4-BUG.c Miért nem a várt módon viselkedik a program? Hogyan kell helyesen megírni? Megjegyzés: csak akkor írjunk ilyet, ha tényleg szükségünk van rá!
  3. (rémálom stílusú forráskód) rondaprogram0.c Mit csinál ez a program? Írd át javítóbarát stílusúvá a forráskódot!
  4. (egész estés program – )

    Töltsd le a választott feladatnak megfelelő adatfájlt:

    1. worldstat.dat (csak az adatok)
    2. latexdoc-b.tex (teljes tex-fájl, a táblázat helyén a nyers adatokkal)

    Az adatfájlt be kell olvasni, és feldolgozva kiírni. Elegáns megoldás a standard bemeteről olvasás és standard kimenetre írás; futtatáskor I/O-átirányítás használata. (Kevésbé elegáns, de működőképes megoldás: megnyitni a be- és kimeneti fájlokat is, mint igazi fájlokat.) Az adatfájl oszlopai: a (int), b (double), c (double). Pl.:
    1038 08 15

    A bemenet egy-egy sorából készítsük el egy LaTeX-táblázat egy-egy sorát úgy, hogy az első oszlop (a) egész érték legyen, a második (b) és a harmadik oszlop (c) pedig normálalakban legyen. Legyen továbbá egy negyedik oszlop is a LaTeX-táblázatban, amely a c / b · 10–6 értéket tartalmazza fixpontos alakban. Azaz a fenti sor az átdolgozás után:
    1038 & $8.00\cdot 10^{0}$ & $1.50\cdot 10^{1}$ & 0.000001875 \\

    Választható feladatok:

    1. (még egyszerűbb) a program a tiszta adatfájlt kapja meg, és biztosra veheti, hogy soronként int double double adatokat talál benne. A feldolgozott fájt kézzel kell bemásolni a LaTeX-dokumentum (latexdoc-a.tex) megfelelő helyére, majd lefordíthatjuk.
    2. (egyszerű) a program minden sor esetén vizsgálja meg, hogy sikerült-e az int double double adatokat beolvasni. Ha igen, a kimenetre a LaTeX-formázott adatokat írja ki. Ha nem, változatlanul írja ki azt a sort, amit beolvasott. (Másoljon.) Így a teljes (nyers adatokat tartalmazó) LaTeX-dokumentumot megetethetjük a programunkkal, majd rögtön fordíthatjuk is. (Természetesen, ha bárhol máshol a dokumentumban be tudja olvasni az int double double adatokat, akkor ott is konvertálni fog, ha kell, ha nem. Emiatt ez a program inkább csak demonstrációs célokat szolgál.)
      Geek feladat: I/O-átirányítás segítségével egyetlen paranccsal oldjuk meg a nyers dokumentum transzformálását és a LaTeX-fordítást.


4. hét

2. gyakorlat (szeptember 26.): Véletlen számok, GSL

  1. Nem szokványos eloszlások – tábla-krétás óra
    1. Véletlen szám generálása adott eloszlásfüggvény szerint
    2. Példa: Pareto-eloszlás
    3. Corrected Pareto-eloszlás
    4. Rejection sampling
    5. Példa: rejection sampling alkalmazása a corrected Pareto-eloszlásra
    6. Empirikus eloszlásfüggvény
  2. Példaprogram a fentiekre + GSL használatára

    Példaprogram Rövid leírás, jegyzet
    gsl_example.c Nulladik GSL-es mintaprogram
    gyak2_rejection_sampling.c

    Doxygenes dokumentáció.

    Fordítási parancssor:
    gcc gyak2_rejection_sampling.c -lm -lgsl -lgslcblas -Wall -o gyak2_rejection_sampling.x
    ahol -lgsl és -lgslcblas jelenti a GSL-lel való linkelést. A GSL használatához a megfelelő header-fájlt be kell építeni a forráskódba: #include <gsl/gsl_rng.h>

    GNU Scientific Library: http://www.gnu.org/software/gsl (Reference Manual: +1 klikk.)


    pareto-test.jpeg A Pareto-eloszlás és a corrected Pareto-eloszlás empirikus eloszlásfüggévnye log–log skálán. (Azaz 1 – G(x).) Mintaméret: N = 10 000. A folytonos vonal az egzakt függvény.

    Az órán elfelejtettem mondani, hogy a példaprogramban fix seed-del (123456) indítottuk a véletlenszám-generátort. Ez arra lehet jó, hogy egy esetleges hibát reprodukálni tudjunk. Egy jó seed-elési módszert lásd az alábbi downloadseed.c-ben.

  3. 2 dimenziós Ising-modell

    Példaprogram Rövid leírás, jegyzet
    downloadseed.c
    downloadseed.h
    ising2d.c
    2 dimenziós Ising-modell

    Dokumentáció

    Fordítási parancssor:
    gcc ising2d.c downloadseed.c -lgsl -lgslcblas -Wall -o ising2d.x
    Opcionális fordítási paraméter (ha a processzor támogatja az SSE ill. SSE2 utasításkészletet): -msse ill. -msse2. (Bár én nem tapasztaltam mérhető gyorsulást... Nektek mi a tapasztalatotok?)


    ising2d-magn.jpeg, ising2d-susc.jpeg Mágnesezettség és szuszceptibilitás a hőmérséklet függvényében kis rendszerekre. (Amiket az órán számoltunk.) Figyeljétek meg, hogy a rendszerméret növelésével nem csak a szuszceptibilitás csúcsának nagysága növekszik, hanem a maximumhelye is eltolódik! Az egzakt kritikus hőmérséklet: TC = ≈ 2,269 185.1


    Az alábbi könyv nem csak hivatkozás, hanem ajánlott irodalom is. (3. fejezet: The Ising model and the Metropolis algorithm)
    1 Newman, Barkema, Monte Carlo Methods in Statistical Physics, p. 82, Oxford University Press, New York, (1999)


5. hét

2. adag kis házi feladat: Véletlen számok, GSL

Megoldások

  1. Írjuk meg az alábbi függvényeket, melyekkel multiplicative congruential (linear congruential) algoritmuson alapuló véletlenszám-generátort valósítunk meg.
    my_rng* my_rng_alloc(int a, int m, int c); // Helyfoglalás a generátornak és inicializálás. (nem seed!) Jelölések a jegyzet szerint. Hiba esetén adjon vissza NULL-pointert!
    void my_rng_set(my_rng* r, int seed); // seed beállítása
    int my_rng_get(my_rng* r); // véletlen egész szám visszaadása a generátor természetes értékkészletéből, melynek határait az int my_rng_min(my_rng* r) és int my_rng_max(my_rng* r) függvényekkel lehet lekérdezni
    int my_rng_min(my_rng* r);
    int my_rng_max(my_rng* r);
    double my_rng_uniform(my_rng* r); // UNI[0;1) eloszlású véletlen szám visszaadása
    void my_rng_free(my_rng* r); // a lefoglalt memória felszabadítása.
    A my_rng legyen egy saját típus, melyben tároljuk az adott generátor paramétereit és belső állapotát. Írjuk meg a függvényeket használó main függvényt, amely a standard kimetre ír egymillió UNI[0,1) eloszlású véletlen számot. (Soronként 1 számot.) (Tipp: a main függvénytől különálló .c és .h fájlokat használjunk a fenti függvények megvalósításához!)
  2. Írjunk programot, ami a standard bemenetről olvas double típusú számokat, a lista hosszát nem adjuk meg előre. Számítsuk ki az adatsorra az első néhány (mondjuk m = 20-ig) momentumot és a

    (auto)korrelációs függvényt t = 0..1000 értékekre. Ne használjunk feleslegesen sok memóriát!
    A bemenet formátuma:
    0.4290764009
    0.3599348921
    ...
    0.3685425046
    0.4074869850
    A kimeneten az adatokat két csoportra bontjuk és számozzuk gnuplot ábrázolás céljaira:
    #Momentumok j=0..m
    0     +1.00000
    1     +0.49766
    ...
     
    #Korrelaciok i=0..k
    0     +0.99516
    1     +0.00492
    ...
    Futtassuk le a tesztet legalább egymillió elemű véletlenszám-listákra, amiket a következő algoritmusokkal állítunk elő: Értelmezzük az eredményeket, melyek a „jó” véletlenszámok? Megbukik-e valamelyik generátor a teszten? (Emlékezzünk: az UNI[0,1) valószínűségi változó momentumai jól ismertek, a korrelációk „elfogadható” értékeire nehezebb jó becslést adni, tekintsük a random.org-ról származó adatokat „jó” véletlenszámnak.)
  3. Bónusz feladat: Adjunk meg olyan valószínűségi változót, (az eloszlássűrűség-függvényt, vagy egy algoritmust) amelynek értékkészlete a [0;1] intervallum vagy részhalmaza; a nulladik, első és második momentuma megegyezik az egyenletes eloszlás momentumaival, azaz 1, 1/2 és 1/3. A magasabb momentumok viszont különböznek az egyenletes eloszlásétól.


6. hét

3. gyakorlat (október 10.)

  1. Fürtszámlálás
  2. Importance sampling
  3. Maxwell-Boltzmann

7. hét

3. adag kis házi feladat: Kérdés

Megoldás nélküli házi feladat :) – Sajnos tényleg... :(


Félévközi kérdőív eredmények


8. hét

4. gyakorlat (október 24.): Konzultáció


9. hét

4. adag kis házi feladat


10. hét

5. gyakorlat (november 7.): Molekuladinamika Verlet-algoritmussal. Példa: 2D rugós kristálymodell

  1. Táblás bevezető

    Vegyünk egy L × L-es négyzetrácsot, periodikus határfeltétellel. ([L] = darab.) A rácsállandó a = 1, rácspontokban lévő atomok tömege m = 1. A szomszédok közé k rugóállandójú rugókat kapcsolunk, melyek elő vannak feszítve: nyújtatlan hosszuk l0. Az (i,j)-edik atom kitérése a rácsponthoz viszonyítva: (xi,j(t), yi,j(t)).

    Felírtuk a mozgásegyenleteket lineáris közelítésben. A megoldást lásd szilárdtestfizika gyakorlaton.

    Példaként felírtuk az (i,j)-edik és az (i-1,j)-edik atomok kölcsönhatásából származó mozgásegyenletet az (i,j)-edik atomra. (A programkódban // left-tel jelölt blokk.) A többi is teljesen hasonlóan jön ki. (Remélem, nem számoltam el… Ha van kompaktabb de még átlátható más ötletetek, elküldhetitek.)

    A differenciálegyenlet-rendszer kezdőfeltétele egy k hullámszámvektorú, ω körfrekvenciájú síkhullám:

    xi,j(t) = u sin( kx · x + ky · y – ω · t )
    yi,j(t) = v sin( kx · x + ky · y – ω · t )
    i,j(t) = – u · ω · cos( kx · x + ky · y – ω · t )
    i,j(t) = – v · ω · cos( kx · x + ky · y – ω · t )

    ahol az (u,v) vektor a hullám polarizációvektora.

  2. Verlet-algoritmus

    Lásd: Ea9.pdf: 7., 8. oldal.

  3. Példaprogram: MD szimuláció Verlet algoritmussal (verlet-crystal.c)


    k = (pi/10, pi/5) hullámszám esetén.


    k = (pi, 0) hullámszám esetén.

  4. Önálló feladat

    Írd át a fenti programot úgy, hogy az origóban lévő atom tömege 10 legyen! Legyen a hullámszám vektor k = (pi,0).
    (verlet-crystal.ver01.m.c)


11. hét

5. adag kis házi feladat: GSL FFT. Példa: 2D rugós kristálymodell

Megoldások

Ez egy kicsit melósabb, utánaolvasósabb házi feladat, de érdemes végigcsinálni, mert tanulságos és szép. :) Motivációként itt az én ábrám az 1. feladat eredményeként:

Megjegyzés: Nem hinném, hogy „élesben” bárki ezzel a módszerrel számolna diszperziós relációt, viszont több hasznos fogást lehet rajta gyakorolni. (Molekuladinamikai szimuláció, FFT, kétváltozós függvényábrázolás, adatfeldolgozás.)

  1. Olvass utána, hogyan kell használni a GSL FFT-jét! GSL Manual 16.0 – 16.3. fejezete. http://www.gnu.org/s/gsl/manual/html_node/Fast-Fourier-Transforms.html Egészítsd ki (és/vagy módosítsd) az órai programot úgy, hogy az origóban lévő atom x kitérését, mint az idő függvényét számítsa ki és transzformálja át frekvenciatérbe. Ha ezt a műveletet a Brillouin-zóna különböző k-pontjaiban elvégezzük, akkor a modell fonondiszperziója feltérképezhető. Az eredményt ábrázold olyan ábrán, amelynek a vízszintes tengelye a hullámszám, függőleges tengelye a körfrekvencia, és a fonon amplitúdóját (a Fourier-együttható abszolút értékét) a színnel ábrázolod. Ehhez ajánlott olvasmány: http://t16web.lanl.gov/Kawano/gnuplot/plotpm3d-e.html#6.7 Ábrázold a diszperziós relációt a Brillouin-zóna Γ → X → M → Γ vonalain! Tipp: állítsd be a színskálát úgy, hogy végig szépen kirajzolódjon a diszperziós ág! Ekkor persze a Γ pont közelében az értékek „túlcsordulnak” a skálán. A Γ pontban nem érdemes számolni.

    Ajánlott paraméterek: előfeszítés: l0 = 0.5; polarizáció: longitudinális. Figyelj a k-térbeli felbontásra: olyan k-pontokat használj, ahol a kezdőállapot megfelel a periodikus határfeltételnek! Figyelj a frekvenciafelbontásra: elegendően nagy felbontást használj, de csak annyira, hogy emberi léptékű maradjon a futásidő. :) Tipp: a Fourier-transzformációhoz nem kell Δt időközönként mintavételezni az x(t) függvényt.

  2. Az 1. feladatban elkészített programot futtasd le úgy is, hogy minden második atom tömege 1, minden másodiké 14/12, sakktáblaszerű elrendezésben! Hol vannak az optikai fononok? Megjelennek az eredményben?

  3. Legyen minden atom tömege 50%-os valószínűséggel 1, 50%-os valószínűséggel 14/12! (Véletlen ötvözet.) Most hogy néz ki a diszperziós reláció? Tipp: tanulságos 2 különböző realizációt végigfuttatni és összehasonlítani a kapott ábrákat.

  4. (Bónusz.) Elég sok paraméterrel lehet játszani, ha valaki talál szép és/vagy tanulságos beállítást, ossza meg velünk a levelezőlistán!


12. hét


13. hét


14. hét

6. gyakorlat (december 5.): Konzultáció


Év végi kérdőív eredmények