Fizika | Tanulmányok, esszék » Csabai István - Molekula dinamika

Molekula dinamika Számítógépes szimulációk 1n4i11/1 Csabai István ELTE Komplex Rendszerek Fizikája Tanszék 5.102 Email: csabaiΘcomplex.eltehu 2009 tavasz Nagy számú részecske szimulációja • A molekula dinamika (Molecular dynamics, (MD)) "részecskék" mikroszkopikus dinamikájának követésével foglalkozik. • Makroszkopikus rendszerben a részecskék száma 6 × 1023 nagyságrend¶. A jelenlegi számítógépekkel nem tudjuk ennyinek a mikroszkopikus dinamikáját követni, de lehetséges elég sok részecskét szimulálni, amelyekkel már termodinamikus tulajdonságok vizsgálhatóak. • Cél, hogy minél több részecskét, minél pontosabban szimuláljunk, ehhez közelítéseket vezetünk be. A részecskék között páronkénti kölcsönhatások vannak, tehát N részecske esetén elvileg N 2 er®hatást kell szimulációs lépésenként kiértékelni. • Ha az er® elég gyorsan lecseng (pl. a molekulák közt tipikusan fellép®

Van der Waals er®), akkor a távoli er®k nagyobb térrészekre kiátlagolhatóak, nagyon távol pedig elhanyagolhatóak, így a módszer jelent®sen gyorsítható. Az ún hosszú hatótávolságú er®knél (pl. a gravitáció) az elhanyagolás nem lehetséges, más közelítéseket kell alkalmazni, ezeket a molekula dinamikától megkülönböztetend®, N-test szimulációnak (N-body) nevezik. Molekula dinamika algoritmus vázlata • • • Legyen N (nagy számú) atomunk, melyeknek ismerjük kezd® sebességüket és pozíciójukat A Newton törvények értelmében a részecskék közt páronként számoljuk ki az er®t, és léptessük a részecskéket Figyeljük, hogy a rendszer a kezdeti tranziens után (nem egyensúlybeli kezdeti feltételek miatt) termodinamikai egyensúlyba került-e. Az egynsúly esetében teljesül az ekvipartíció:   hKi = • 1 mv 2 2 3 = kT 2 ahol k = 1.38 × 10−23 J/K; a Boltzmann állandó A szimuláció során hKi mutatja, hogy

egyensúlyban van-e a rendszer. Ha a rendszer elérte az egyensúlyt, mérhet®ek a termodinamikai változók: h®mérséklet (T ), nyomás (p), h®kapacitás (CV ), stb. Közelítések • Elvileg az atomok bels® szerkezetét is gyelembe vev® kvantummechanikát kellene használni. De pl T = 300K argon gázban az egy részecskére jutó mozgási energia 3 kT ∼ = 6.2 × 10−21 J = 0039eV 2 sokkal kisebb mint az elektronok gerjesztési energiája (11.6eV) Valamint az m = 669 × 10−26 kg tömeg¶ atomok de Brolglie hullámhossza 2π~ 2π~ ∼ =√ = 2.2 × 10−11 m, p 3mkT jóval kisebb mint az atom eektív mérete (r0 = 3.4 × 10−10 m), így az atomok közötti tipikus távolságnál is. A nemesgázoknak (pl az argon) zártak az elektronhéjai, így szimulációjuk során a fentiek gyelembe vételével, a klasszikus részecske közelítés elfogadható eredményre vezet. • Megviszgálva az id®skálákat, probléma lehet a termodinamikai egynsúly elérése

is. A szimulációs lépéshossznak kisebbnek kell lenni, mint az atomok közötti átlagos távolság megtételéhez szükséges id®, az atomok átlagos sebességével, ami tipikusan: r0 p 3kT /m ∼ = 7.9 × 10−13 s, azaz kevesebb mint 1 picosecundum. Ha N = 103 atomra minden lépésben N (N − 1)/2 er®t kell kiszámolni, akkor egy gyors gépen is órákig tarthat akár néhány nanosecundum szimulációja! Ezért fontos a kezdeti feltételeknek a termodinamikai egyensúlyhoz minél közelebbi beállítása. A Verlet algoritmus • • A bolygómozgásnál már megismert módon most is dierenciálegyenlet rendszert kell integrálni. Természetes választás lenne, az Euler-, vagy az annál pontosabb Runge-Kutta módszer. A legelterjedtebb MD dierenciálegyenlet léptet® algoritmus mégis inkább a Verlet algoritmus. • Jelölje R(t) = (r1 , . , rN ) a részecskék koordinátáit, hasonlóan V(t) a sebességeket, A(t) pedig a gyorsulásokat. A Verlet

algoritmus egy lépése: Rn+1 = 2Rn − Rn−1 + τ 2 An + O(τ 4 ) Vn = • Rn+1 − Rn−1 + O(τ 2 ) 2τ El®nyök: • Gyorsabb mint a Runge-Kutta, egy lépés során, csak egyszer kell kiszámolni a gyorsulásokat • Majdnem olyan pontos, mint a 4-ed rend¶ Runge-Kutta, ahol a hiba O(τ 5 ) rend¶ • Jól meg®rzi az energiát, ami nagyon fontos a sokrészecske dinamikánál • Id®tükrözésre nem változik, ami fontos a részletes egyensúly (ergodicitás) feltétele miatt Rn+1 = 2Rn − Rn−1 + τ 2 An + O(τ 4 ) Rn+1 − Rn−1 Vn = + O(τ 2 ) 2τ • Hátrányok: • Három pontos rekurzió, azaz két el®z® lépést használ, nem indítható egy kezdeti feltételb®l • A sebesség hibája O(τ 2 ), viszont ha az er® nem függ a sebességt®l, akkor ez nem rontja a pozíciók pontosságát. • A sebesség és a pozíció nem egyszerre frissít®dik, a sebesség "le van maradva" • A problémák egy részén segít a velocity-Verlet

algoritmus: Rn+1 = Rn + τ Vn + τ2 An + O(τ 3 ) 2 τ Vn+1 = Vn + (An+1 + An ) + O(τ 3 ) 2 • Ez R-ben csak O(τ 3 ) rendben pontos, de nem használ két el®z® lépést, egyszerre frissíti a koordinátákat és sebességeket, valamint a sebesség is pontosabb. A 2 dimenziós Lennard-Jones rendszer • • Modelezzük N darab "2 dimenziós argon atom" viselkedését egy négyzetben. Azért, hogy a határokon ne kelljen speciális feltételekkel foglalkozni (fallal ütközések, nyílt rendszernél részecskeszám megmaradás), alkalmazzunk periodikus határfeltételt: azonosítsuk a szemközti oldalakat (tórusz toplógia), így egy "határtalan" rendszert kapunk. A részecskék közötti kölcsönhatást a van der Waals közelítésben a Lennard-Jones potenciál írja le: V (r) = 4V0 •   r0 12  r0 6 − , r r ahol r az atomok középpontja közötti távolság. Argonra V0 = 1.65 × 10−21 J azaz V0 /kB = 1198K és r0 = 3.41 × 10−10 m

A potenciál gradienséb®l számolható er®: ~ F(r) = −∇(r) =     24V0 r0 12  r0 6 2 − r r2 r r Forrás: Ahttp://www.doksihu Lennard-Jones potenciál r0 = 1 és V0 = 1 esetén: • • • • A potenciál r = r0 -nál 0 és r = 21/6 r0 -ben veszi fel minimumát, ami −V0 At er® er®sen taszító kis távolságokon a zárt elektronhéjak Pauli-féle kizárási "taszítása" miatt, nagy távolságokon pedig gyengén vonzó, az indukált elektromos dipóler®k (van der Waals kölcsönhatás) miatt A rendszer energiája: E= N X 1 i=1 • 2 mvi2 + X V (rij ), hiji ahol rij az hiji atompár közötti távolság. Egyensúlyban a sebességek eloszlásának a Maxwell-Boltzmann eloszlást kell követni, ami 2 dimenzióban: P(v)dv =  m  mv2 e− 2kT vdv. kT 2D Lennard-Jones szimuláció • A szimuláció során a valódi SI értékekkel numerikus szempontból nem praktikus számolni, ezért a változókat átskálázzuk (hasonló mint a

realtivitáselméletben szokásos c = 1 mértékegységrendszer.) Úgy választjuk meg a skálát, hogy az argon atom tömege legyen a tömegegység m = 1, r0 = 1 és V0 = 1. A rendszer karakterisztikus idejét s τ0 = • mr02 = 2.17 × 10−12 s, V0 ami pár picosecundum, szintén skálázzuk át, τ0 = 1. Vegyünk például N = 16 atomot, az id®lépés pedig legyen a karakterisztikus id® századrésze, τ = 0.01 Az egy részecskére ható er® kiszámolásához, össze kell adni az összes többi részecske hatását: Fi = X j =1 j 6= i Fj hat i−re • Az algoritmus • • • fontos új része a periodikus határfeltétel kezelése. Mivel az elemi L3 térfogatban lév® részecskék periodikusan ismétl®dnek, meg kell oldani, hogy ne legyen végtelen sok kölcsönhatás. A részecske saját képeinek hatása kioltja egymást, de más részecskék többszörösen több irányból hathatnának, ami olyan korrelációkat hozna be ami nem lenne zikai. A

legközelebbi kép nem mindig az aktuális cellában van, lehet annak valamelyik els® szomszédjában is. (Belátható, hogy távolabb nem lehet.) Például, lehet, hogy rij = |ri − ri | nagyobb mint egy szomszédos cellában lév® rij 0 = |ri − r0 i | Hogy egy adott részecskére minden más részecske csak egyszeresen hasson, válasszuk mindig a legközelebbi megjelenését (legnagyobb er®hatást), a fenti 1D példán i-hez a j 0 -t a j helyett. Gondoskodni kell a kódban arról is, hogy ha egy részecske pl. jobb oldalon elhagyja a négyzetet, akkor bal oldalon visszakerüljön, stb. a kezd®feltételek beállítása se. Ha nem realisztikus elrendezést és sebességeket állítunk be, akkor sok id®be telik míg a kezdeti tranziens után egyensúlyba kerül a rendszer. • Nem nyilvánvaló • Kondenzált állapotban az atomok általában lapcentrált köbös rács pontjaiban helyezkednek el, ez minimalizálja a potenciális energiát. Az egyszer¶ változatban e

helyett egyszer¶ köbös rácsot használunk. • A sebességeknek a Maxwell-Boltzmann eloszlást kellene követni, e helyett egyel®re egy adott intervallumban eloszló egyenletes véletlen értékekr®l indítjuk a sebességeket. Pillanatnyi h®mérséklet A Lennard-Jones er® konzervatív, rendszer teljes energiája megmarad. Termikus egyensúlyban teljesül az ekvipartíció tétele, amikor is a kinetikus energiából származtatható a h®mérséklet: 1 3(N − 1) × kB T = 2 * N mX 2 vi 2 + . i=1 h. i jelöli a termikus sokaság átlagot 3(N − 1), a bels® szabadsági fokok száma, a tömegközéppont mozgása nem járul hozzá a termikus energiához, ha az átlag sebesség v̄ 6= 0, akkor igazából vi2 helyett (vi − v̄)2 szerepel a kifejezésben. Nem-egyensúlyi helyzetben is lemérhetjük ezt a mennyiséget, ekkor a nem-egyensúlyi (nem igaz az ekvipartíció) jelleg megkülönböztetéseként pillanatnyi h®mérsékletnek nevezhetjük. Programok • •

A md.cpp program a legegyszer¶bb 3 dimenziós Lennard-Jones szimuláció. Nem használ periodikus határfeltételeket, a velocity-Verlet algoritmussal léptet, és a pillanatnyi h®mérsékletet egy fájlba menti ki. A md-gl.cpp program a velocity-Verlet algoritmust használja periodikus határfeltételek mellett. Folyamatosan méri a sebességhisztogramot, és összeveti a Maxwell-Boltzmann eloszlással. Az egérgombok segítségével indítható/állítható a szimuláció, illetve növelhet®/csökkenthet® a h®mérséklet. Fejlesztések a kódon • • A 3 dimenziós kódban is alkalmazzuk a periodikus határfeltételt! Alacsony h®mérsékleten, amikor a potenciális energia minimuma körül van a rendszer, a rendszer a lapcentrált köbös állapotban éri el az energiaminimumot. Az elemi cella 4 atomot tartalmaz (piros): • Akkor tudunk jól kitölteni egy kockát, ha az atomok száma N = 4M 3 , M = 1, 2, 3, . , pl 32 = 4 × 23 , 108 = 4 × 33 , 256, 500, 864, .

• Az 1 méret¶ elemi cellában lév® atomok pozíciói (0, 0, 0) (0.5, 05, 0) (05, 0, 05) (0, 05, 05) de hogy ne a határon legyenek toljuk el (0.5, 05, 05)-tel, így a kódban: // 4 atomic positions in double xCell[4] = {0.25, double yCell[4] = {0.25, double zCell[4] = {0.25, fcc unit cell 0.75, 075, 025}; 0.75, 025, 075}; 0.25, 075, 075}; Kezdeti sebességek Maxwell-Boltzmann eloszlás szerint • T h®mérsékleten a Maxwell-Boltzmann eloszlás  P(v) = m 2πkB T 3/2 e − 2 +v 2 +v 2 ) m(vx y z 2kB T . √ • • vagyis minden sebességkomponens 0 várható érték¶, T -vel arányos szórású Gauss eloszlással írható le. A kódban használt a Numerical Recipes-b®l átvett gasdev() függvény 0 várható érték¶, 1 szórású véletlen számokat állít el® a gépi egyenletes eloszlásokból a Box-Muller algoritmus segítségével. Az így generált kezd®sebességek átlaga 0-hoz közeli lesz, de nem egzaktul 0, a tömegközéppont sebességével

korrigálni kell a kezd®értékeket: vi := vi − vCM , ahol a tömegközépponti sebesség: P N vCM = Pi=1 N mvi i=1 m • Ezek után az 1-re normált szórású sebességeket a kívánt sebességre kell átskálázni hogy a kívánt T h®mérsékletet kapjuk: vi λvi s 2(N − 1)kB T λ= PN 2 i=1 mvi • A md2.cpp program tartalmazza a fenti javításokat, gyeli a pillanatnyi h®mérsékletet, és amennyiben nem felel meg a megjelölt értéknek, újra átskálázza a sebességeket. Gyorsítási lehet®ségek • • • A program legid®igényesebb része a páronkénti er®k kiszámítása, N (N − 1)/2 pár, azaz O(N 2 )-tel skálázik a futási id®. L. Verlet egy cikkében (Phys Rev 159, 98 (1967)) két gyorsítást is bevezetett: A potenciál levágása: Mivel a Lennard-Jones er® rövid hatótávú, nagysága r > r0 után gyorsan csökken, így bevezethet® egy levágási távolság rcutof f , amin túl praktikusan 0-nak tekinthet®. Így, ha az rcutof

f távolságon belül K részecske van és N = n × K , akkor a skálázás tesz®leges n-re O(n × K 2 ) rend¶. Máshogy megfogalmazva ha a s¶r¶séget állandóan tartjuk akkor n-ben nem négyzetes, hanem lineáris a skálázás. Nyilván akkor éri ez meg, ha rcutof f << L (L a rendszer mérete). A levágás érvénysítéséhez tudni kell, melyek azok a szomszédok, amelyekre rij = |ri − rj | < rcutof f . Mivel minden részecske mozog, ennek kiértékelése is O(N 2 ) feladatnak t¶nik. Ha azonban gyelembe vesszük, hogy egy-egy lépésben minden részecske keveset mozdul el, és választunk egy megfelel® L >> rmax > rcutof f távolságot, amelynél jobban néhány lépésen belül nem távolodnak el a részecskék. Így elég azon belül nyilvántartani a szomszédsági listát, és csak néha frissíteni a kölcsönható párok listáját. Verlet rcutof f = 25r0 , rmax = 32r0 értékeket javasolta, és 10 lépésenként frissítette a szomszédsági

listát. (A frissítés gyakorisága és rmax közötti összefüggés a sebességek valószín¶ségeloszlásából becsülhet®). Verlet így 10-szeres sebességnövekedést ért el a pontosság jelent®s romlása nélkül. A levágás kisebb hibákat behoz a szimulációba és elrontja az energia megmaradást. A hiba korrigálható a potenciál módosításával: • Szomszéd-lista: • • Ucorr (r) = U (r) − d U (rcutof f )(r − rcutof f ) dr (Nincs a kiadott kódokban implementálva.) Gyorsított programok • • A md3.cpp program tartalmazza a fenti Verlet-féle közelítésesket, gyorsítást. A md3-gl.cpp program a fenti kód 3 dimenziós OpenGL változata. Termodinamikai mennyiségek • Teljes energia (mozgási és potenciális energiák összege): N E= mX 2 X vi + U (|ri − rj |) 2 i=1 • i6=j H®kapacitás (a uktuáció-disszipáció tétel alapján):  CV = ∂E ∂T  = V i 1 h 2 2 E − hEi kB T 2 Az átlagolás elvileg sokaságátlag

(pl. sok különböz® véletlen futtatás eredménye), de egyensúlyban ezt jól közelíti az id®átlag. A h®kapacitást a h®mérséklet-uktuációkból is becsülhetjük: T 2   3 3N kB 2 − hT i = N (kB T ) 1 − 2 2CV 2 • Véges zárt rendszernél a nyomást mérhetnénk a doboz falánál történt impulzusváltozásokból is, de a Viriál-tétel alapján a következ® kifejezés is használható: 1 P V = N kB T + 3 * X + rij · Fij i<j A Kompresszibilitási faktor kifejezi, mennyire vagyunk távol az ideális gáz állapottól: 1 PV =1− Z= N kB T 3N kB T * X + rij · Fij i<j Nagy s¶r¶ségen a taszító potenciál miatt Z > 1, kisebb s¶r¶ségeken pedig a vonzó van der Waals er®k miatt Z < 1. Szabad részecskékre, ideális gázra lenne Z = 1. Feladatok 1 2 3 Értsük meg az md.cpp programot és ábrázoljuk kimenetét! Értelmezzük a kapott görbét! Vizsgáljuk meg az md2.cpp programban a javításokat. A program futtatási

eredményét elemezve vizsgáljuk meg, hogy mennyit segített a kezdeti feltételek pontosabb beállítása, és mennyi id® alatt áll be a kívánt h®mérséklet! Értsük meg a md3.cpp program hogyan kezeli a levágásokat és a szomszédsági listát, és annak frissítését. Kísérletezzünk más rcutof f , rmax és updateInterval paraméterekkel, nézzük meg hogyan változik a sebesség és a pontosság az eredeti md2.cpp programhoz képest. Írjuk meg a teljes energiát, nyomást, h®kapacitást, kompresszibilitást számoló függvényeket! Próbáljuk meg a fázisátalakulási h®mérsékletet meghatározni a rendparaméterek (h®kapacitás, kompresszibilitás) gyelésével. Vessük össze az értékeket a valódi argon olvadás és forráspontjával! Hogyan befolyásolja N értéke a tapasztaltakat? Szorgalmi feladatok: Lapozz! Szorgalmi feladatok: 4 5 Alakítsuk át a programot úgy, hogy a periodikus határfeltétel helyett zárt, kemény fala legyen, ahonnan

visszapattannak a részecskék! (Egyszer¶ megoldás a sebességek és pozíciók tükrözése, kicsit bonyolultabb a falat is tömör Lennard-Jones potenciállal leírható anyagnak venni.) Mérjük a nyomást a falakra kifejtett er® segítségével, és ellen®rizzük mennyire teljesül a pV = N kB T törvény! Értsük meg a md3-gl.cpp programban a 3D graka használatát! Szerkesszük össze a kepler.cpp programmal úgy, hogy a gravitációs 3-test (vagy soktest) probléma dinamikáját ábrázolja. Vajon alkalmazható itt is egy levágási hossz a kód gyorsítására?