F
AKULTETA ZA KEMIJO IN KEMIJSKO TEHNOLOGIJODIPLOMSKO DELO
Jaka Prelog
Ljubljana, 2021
F
AKULTETA ZA KEMIJO IN KEMIJSKO TEHNOLOGIJOUNIVERZITETNI ŠTUDIJSKI PROGRAM 1. STOPNJE KEMIJA
Monte Carlo simulacije mehkih delcev
DIPLOMSKO DELO
Jaka Prelog
M
ENTOR: prof. dr. Tomaž Urbič
Ljubljana, 2021
IZJAVA O AVTORSTVU
diplomskega dela
Spodaj podpisani Jaka Prelog sem avtor diplomskega dela z naslovom:
Monte Carlo simulacije mehkih delcev (angl. Monte Carlo simulations of soft particles).
S svojim podpisom zagotavljam, da:
je diplomsko delo rezultat mojega raziskovalnega dela pod mentorstvom prof. dr. Tomaža Urbiča;
sem poskrbel, da so dela in mnenja drugih avtorjev, ki jih uporabljam v predloženem diplomskem delu, navedena oziroma citirana v skladu z navodili;
se zavedam, da je plagiatorstvo, v katerem so tuje misli oziroma ideje predstavljene kot moje lastne, kaznivo po zakonu (Zakon o avtorskih in sorodnih pravicah – uradno prečiščeno besedilo (ZASP-UPB3) (Ur. list RS, št. 16/2007);
sem poskrbel za slovnično in oblikovno korektnost diplomskega dela;
je elektronska oblika diplomskega dela identična tiskani obliki diplomskega dela.
V Ljubljani, 30. avgusta 2021 Podpis avtorja:
Zahvala
»How can giving be better than receiving? You easily forget what you receive but you do not forget when you give. That is why the joy even of giving thanks lasts a long time.«
(J. M. Seok)
V tej zahvali se želim spomniti vseh, ki ste na kakršenkoli način pripomogli k nastanku diplomskega dela ali mi drugače pomagali v času mojega študija. Diplomsko delo sem opravljal na Katedri za fizikalno kemijo Fakultete za kemijo in kemijsko tehnologijo, Univerze v Ljubljani, pod mentorstvom prof. dr. Tomaža Urbiča, zato gre najprej iskrena zahvala prav Vam. Sprejeli ste me z odprtimi rokami, me med nastajanjem diplomskega dela spodbujali, usmerjali, strokovno vodili in mi bili vselej pripravljeni pomagati.
Nenazadnje se vam moram zahvaliti za vso potrpežljivost in razlago nerazumljivega.
Zahvaljujem se tudi prof. dr. Juriju Reščiču za strokovni pregled diplomske naloge in napotke za njeno izboljšavo. Za slovnični in pravopisni pregled se zahvaljujem Mojci Karnet, univ. dipl. prof. slov., za pregled kratkega povzetka v angleščini pa Mini Mazouzi, univ. dipl. prof. angl. in nem.
Posebna zahvala gre tudi mojim staršem, mami in očetu, ker sta mi omogočila študij v Ljubljani. Hvala vama za vso podporo, spodbude in nesebično (finančno) pomoč v času študija!
Hvala tudi mojemu dragemu sostanovalcu in dobremu prijatelju Klemnu za vse računalniške nasvete. Bil si mi v veliko pomoč. Prav tako hvala vsem, ki ste me kakorkoli podpirali in mi stali ob strani na moji izobraževalni poti.
Povzetek: Z izrazom »mehka snov« mislimo na snovi, kot so koloidi, tekoči kristali, surfaktanti, polimeri itd. Lastnosti nekaterih od teh sistemov je mogoče modelirati s pomočjo Gaussovega potenciala. V tem diplomskem delu bomo predstavili rezultate, dobljene po modelu Gaussove sredice (angl. the Gaussian-core model, GCM). V ta namen smo izvedli obsežne Monte Carlo simulacije v kanoničnem (izohorno- izotermnem) – 𝑁𝑉𝑇 in izotermno-izobarnem – 𝑁𝑝𝑇 ansamblu ter opravili analizo termodinamičnih in strukturnih lastnosti za posamezni ansambel.
Ključne besede: računalniška simulacija Monte Carlo, mehki delci, model Gaussove sredice, ansambel, termodinamične in strukturne lastnosti
Monte Carlo simulations of soft particles
Abstract: The term »soft substance« refers to substances such as colloids, liquid crystals, surfactants, polymers, etc. The properties of some of these systems can be modelled using the Gaussian potential. In this work, we will present the results obtained according to the Gaussian-core model (GCM). For this purpose, we performed extensive Monte Carlo simulations in the canonical (isochoric-isothermal) – 𝑁𝑉𝑇 and isothermal-isobaric – 𝑁𝑝𝑇 ensemble and made an analysis of thermodynamic and structural properties for each ensemble.
Keywords: Monte Carlo simulations, soft particles, the Gaussian-core model, ensemble, thermodynamic and structural properties
Kazalo
1 Uvod ... 1
1.1 Lennard-Jonesov potencial ... 2
1.2 Model Gaussove sredice ... 2
1.3 Simulacija polimerov ... 3
2 Namen dela ... 5
3 Teoretične osnove ... 7
3.1 Postulata statistične termodinamike ... 8
3.2 Ansambli ... 9
3.2.1 Termodinamika kanoničnega ansambla ... 11
3.2.2 Termodinamika izotermno-izobarnega ansambla ... 13
3.3 Zasnova računalniške simulacije Monte Carlo ... 14
3.3.1 Generiranje začetne konfiguracije ... 14
3.3.2 Problem makroskopskega števila delcev ... 15
3.3.3 Uravnoteženje sistema ... 17
3.3.4 Produkcijski del računalniške simulacije Monte Carlo ... 18
3.4 Metropolisov algoritem ... 19
3.5 Termodinamične in strukturne lastnosti ... 23
3.5.1 Reducirane količine ... 23
3.5.2 Izračun termodinamičnih količin ... 23
3.5.3 Radialna porazdelitvena funkcija ... 25
4 Rezultati ... 27
4.1 Tehnične podrobnosti simulacije ... 27
4.2 Strukturne lastnosti mehkih delcev ... 28
4.3 Posnetki sistemov ... 32
4.4 Termodinamične lastnosti ... 35
5 Zaključek ... 39
6 Literatura ... 41
Seznam uporabljenih kratic in simbolov
Oznaka Pomen
𝑘𝐵 Boltzamnnova konstanta 1,38∙10‒23 JK‒1
ℎ Planckova konstanta 6,626∙10‒34 Js
𝜎 enota za dolžino 𝜀 enota za energijo 𝑁 število delcev 𝑝𝑁
⃑⃑⃑⃑⃑ set gibalnih količin 𝑁 delcev 𝑟𝑁
⃑⃑⃑⃑ set prostorskih koordinat 𝑁 delcev
< > ansambelsko povprečje 𝐸 celotna energija sistema 𝑈 notranja energija sistema 𝐴 Helmholtzova prosta energija 𝐺 Gibbsova prosta entalpija
𝑆 entropija
𝜇 kemijski potencial
𝑀 število generiranih konfiguracij
𝑄𝑁𝑉𝑇 statistična vsota kanoničnega ansambla
∆𝑁𝑝𝑇 statistična vsota izotermno-izobarnega ansambla 𝑟𝐶 mejni radij
ℋ Hamiltonian
𝑔(𝑟) radialna (parska) porazdelitvena funkcija 𝜌 številčna gostota
𝑃𝑖 verjetnost za kvantno stanje 𝑖
𝑍𝑁𝑉𝑇 konfiguracijski integral v kanoničnem ansamblu 𝑟𝑖𝑗 razdalja med delcema 𝑖 in 𝑗
𝐿 dolžina roba simulacijske celice 𝑉∗ reducirana prostornina
𝑇∗ reducirana temperatura 𝑝∗ reduciran tlak
𝐶𝑉 toplotna kapaciteta pri konstantni prostornini 𝐶𝑝 toplotna kapaciteta pri konstantnem tlaku
𝐻 entalpija
𝜅 izotermna stisljivost 𝛼 razteznostni koeficient 𝜏 čas opazovanja sistema
𝑁(𝑠) verjetnostna gostota, da sistem najdemo v stari konfiguraciji 𝑁(𝑛) verjetnostna gostota, da sistem najdemo v novi konfiguraciji 𝜈(𝑠) potencialna energija starega stanja
𝜈(𝑛) potencialna energija novega stanja 𝜉 naključno število
1
1 Uvod
»What do we mean by soft matter? Americans prefer to call it 'complex fluids.' This is a rather ugly name, which tends to discourage young students.«
(Pierre-Gilles de Gennes) Urejenost delcev v neki snovi je pogosto rezultat odbojnih in privlačnih sil. Vse sile so v osnovi elektrostatske in temeljijo na Coulombovem zakonu, vendar izraz »elektrostatske«
navadno uporabljamo v ožjem smislu, ko govorimo o silah med nabitimi delci. Če so sile močne, je urejenost tridimenzionalna (npr. v kristalih); če so šibke, pa ostane le lokalna urejenost, ki se s povečevanjem razdalje med delci izgubi. Za nekatere snovi se danes pogosto uporabljata izraza »trda« oz. »mehka« snov (angl. hard matter oz. soft matter).
Kriterij, s katerim ločujemo med trdo in mehko snovjo, je energija vezi med delci. Če je energija primerljiva s termično, 𝑘𝐵𝑇, kjer je vrednost temperature 𝑇 blizu sobne, snov proglasimo za mehko; v kolikor je energija večja od termične, pa za trdo. Izraz »mehka snov« temelji na makroskopskih lastnosti teh snovi, kamor štejemo koloide, tekoče kristale, površinsko aktivne snovi (surfaktante) in polimere. To so snovi, ki lahko pod določenimi pogoji tečejo, to pa pomeni, da je urejenost v njih šibka [1, 2].
Sile, ki vodijo do urejanja delcev v sistemu, so lahko odbojne ali privlačne. Odbojne sile so kratkosežne in pogosto sorazmerne z 𝑥−2 (Coulombove sile med ioni)1, privlačne pa daljnosežne in pogosto odvisne od 𝑥−6. Ker so vse sile gradienti potencialne energije, nas zanima, kako je potencialna energija odvisna od razdalje [2].
Odbojne in privlačne sile med delci pa ne odločajo zgolj o tem, ali bo sistem termodinamsko stabilen ali ne, ampak tudi o njegovi strukturi na mikroskopskem in makroskopskem nivoju. Struktura sistema je tista, ki narekuje njegove termodinamske lastnosti. Če med delci v sistemu prevladujejo močne odbojne interakcije, je tak sistem termodinamsko stabilen. Delci se v tem sistemu uredijo v kristalinične strukture. Če prevladujejo močne privlačne interakcije, je sistem termodinamsko in kinetično nestabilen. Delci tvorijo prepletene mrežne fraktalne strukture. V praksi je zelo pomembna situacija, ko so med delci privlačne in odbojne interakcije po velikosti primerljive. V nasprotju s tem, kar bi pričakovali, se izkaže, da lahko delci v tem primeru tvorijo zelo goste trdne sklade. Praksa pokaže, da kratkosežni »mehki« odboj med delci in šibek privlak na velikih razdaljah vodita do zelo efektivnega zlaganja delcev. Eden izmed poglavitnih ciljev takšnih študij je napovedovanje in zagotavljanje želenih struktur v preučevanih sistemih. Ravnotežna termodinamika sicer pove, kaj je možno in kaj ni, vendar pa ima v praksi na koncu zadnjo besedo kinetika [2].
1 Taka zveza (oz. odboj) navadno velja med dvema izoliranima ionoma ali molekulama v vakuumu. Če se delca približata na razdaljo, kjer nastopijo Paulijeve odbojne sile, je tak odboj bolj strm.
Monte Carlo simulacije mehkih delcev Jaka Prelog
2
1.1 Lennard-Jonesov potencial
Matematični model, ki opisuje potencialno energijo interakcij med dvema nevtralnima delcema glede na njuno razdaljo, je Lennard-Jonesov potencial. Pravzaprav gre za parski potencial, ki razmeroma dobro obravnava karakteristike sistemov, kjer med delci prevladujejo van der Waalsove interakcije. Enačba (1.1) je uporabna tam, kjer se srečamo s hkratnim delovanjem privlačnih in odbojnih van der Waalsovih sil med delci [2]:
𝑈𝐿𝐽(𝑟) = 4𝜀 [(𝜎 𝑟)
12
− (𝜎 𝑟)
6] (1.1)
kjer je 𝑟 razdalja med delcema, parametra 𝜀 in 𝜎 pa opisujeta globino potencialne jame (𝜀) in razdaljo (𝜎), kjer je je potencial enak 0 [3]. Ta dva parametra uporabimo za formulacijo reduciranih količin.
Prvi člen v enačbi (1.1) predstavlja kratkosežne odbojne interakcije in pada z dvanajsto potenco2. Odbojne sile tega tipa ustrezajo interakcijam, ki so posledica izključenega volumna (in ne elektrostatskih sil). Pojav je znan kot Paulijevo izključitveno načelo (ali Paulijeva prepoved), zaradi katerega delcev ne gre poljubno stiskati skupaj, saj bi v nasprotnem primeru prišlo do prekrivanja elektronskih ovojnic. To je tudi eden od razlogov, zakaj ne pademo skozi tla. Drugi člen v enačbi (1.1) pa predstavlja daljnosežne privlačne interakcije in pada s šesto potenco razdalje med delci. Lennard-Jonesov potencial poznamo tudi pod imenom 12–6 potenčni zakon [2].
1.2 Model Gaussove sredice
Tipični primer sistema delcev mehke snovi je suspenzija mezoskopskih delcev (tj. delcev, katerih velikost je reda od 1 μm do 1 nm). Dejstvo, da je tak sistem večji od sistema atomskih delcev, nam lahko olajša eksperimentalne raziskave, saj je mogoče z ustreznimi spremembami sistem mezoskopskih delcev poljubno prilagajati. Ker pa so delci mehkih snovi navadno kompleksni agregati, zgrajeni iz več tisoč atomov ali molekul, je pri raziskavah nemogoče upoštevati vse sile, ki delujejo nanje [4].
Zaradi zapletene notranje strukture mezoskopskih agregatov lahko ti delci prodirajo, se prekrivajo ali stisnejo. Te posebne značilnosti vodijo do nepričakovanih učinkov na njihove strukturne lastnosti, pa tudi na fazno obnašanje, in jih je potrebno pri izračunih upoštevati [4].
Take sisteme je mogoče modelirati s pomočjo Gaussovega potenciala z uporabo modela Gaussove sredice (angl. the Gaussian-core model, GCM).
2 Eksponent pri kratkosežnih odbojnih interakcijah je v območju 9–15; matematično najprimernejši eksponent je 12. Predznak teh interakcij je vedno pozitiven.
3
V pričujočem diplomskem delu bomo za simulacije strukturnih in termodinamičnih lastnosti modelnega sistema mehkih delcev uporabili model Gaussove sredice, ki ga opiše enačba (1.2) [5]:
𝑈𝐺 = 𝜀 𝑒−(𝜎𝑟)
2 (1.2)
pri tem je 𝜀 energijski parameter (pove nam velikost vrha Gaussove funkcije) in 𝜎 velikostni parameter, ki pove, kako hitro Gaussova funkcija pada z razdaljo [5, 9].
Tako imenovani model Gaussove sredice (GCM) je prvotno predstavil Stillinger [8], ki je zapisal, da se delci po modelu GCM obnašajo kot trde krogle z vedno večjim premerom, ko temperatura in gostota limitirata proti 0. Takšen potencial je kljub dejstvu, da je omejen pri popolnem prekrivanju med delci, vseeno sprejemljiv kot učinkovit potencial [6]. Med drugim lahko opazimo celo vrsto nenavadnih pojavov in lastnosti, kot je recimo nenormalna odvisnost korelacijskih funkcij od gostote ali pojav, kjer lahko z nenehnim spreminjanjem spremenljivke stanja (navadno temperature, tlaka ali koncentracije) določena faza med enim prehodom izgine [7].
Model Gaussove sredice, predstavljen v enačbi (1.2), si zasluži pozornost iz vsaj treh razlogov [8]:
1. Model lahko prikaže vedenje relativno »mehkih« interakcij, podobnih tistim pri resničnih snoveh, ki tvorijo plastične kristale (z nizko talilno entropijo in veliko trdnostjo trdne faze);
2. Pri modelu lahko ovrednotimo izraze v asimptotah visokotemperaturnih vrst proste energije;
3. S pomočjo tega modela razmeroma enostavno ustvarjamo natančne zbirke podatkov z računalniškimi simulacijami [8].
1.3 Simulacija polimerov
Model Gaussove sredice nima odbojnih trdnih jeder, ki so navadno prisotna v modelih fluidov, zato je uporaben za simulacijo in modeliranje sistemov, sestavljenih iz velikih mehkih molekul, kot so recimo polimeri [9].
Polimeri so makromolekule, ki so med seboj povezane iz posameznih fragmentov. Pri simulacijah polimerov si moramo pomagati z različnimi (približnimi) modeli, saj je polimere zaradi velikega števila prostostnih stopenj relativno težko simulirati. Primer takega modela je mrežni model, pri katerem prostostne stopnje, kot so recimo dolžine
Monte Carlo simulacije mehkih delcev Jaka Prelog
4
vezi, in kote zamrznemo, ostanejo pa interagirajoče enote (efektivni monomeri), ki lahko zavzamejo le mesta na tridimenzionalni mreži [10].
Obravnavanje polimerov kot »mehkih« delcev je že stara ideja. Polimerne tuljave v raztopini, zvezdasti polimeri in mikrogeli so tipični primeri mehkih delcev. Zaradi mehkobe so ti sistemi vsestransko uporabni za tehnološke aplikacije, pa tudi za izdelavo novih pametnih materialov [9].
5
2 Namen dela
Namen in cilj diplomskega dela je pregledati in pojasniti nekatere vidike in značilnosti uporabe modela Gaussove sredice za simuliranje mehkih delcev z računalniškimi simulacijami Monte Carlo. V ta namen smo izvedli obsežne Monte Carlo simulacije v kanoničnem (izohorno-izotermnem) – 𝑁𝑉𝑇 in izotermno-izobarnem – 𝑁𝑝𝑇 ansamblu.
Želeli smo ugotoviti, kako spreminjanje temperature in gostote vpliva na strukturo kanoničnega sistema oziroma kako spreminjanje temperature in tlaka vpliva na strukturo izotermno-izobarnega sistema.
S pomočjo računalniških simulacij Monte Carlo smo prav tako želeli pridobiti informacije o termodinamičnih lastnostih mehkih delcev, kot so razteznostni (toplotni) koeficient (𝛼), toplotna kapaciteta pri konstantnem tlaku (𝐶𝑝), gostota sistema (𝜌) in koeficient izotermne stisljivosti (𝜅), in sicer pri spreminjanju tlaka in temperature v izotermo- izobarnem ansamblu, ter lastnostih, kot so tlak (𝑝), povprečna energija na delec (𝑈/𝑁) in toplotna kapaciteta pri konstantni prostornini (𝐶𝑉), pri spremembah temperature in gostote sistema v kanoničnem ansamblu.
7
3 Teoretične osnove
Statistična termodinamika predstavlja povezavo med makroskopskimi lastnostmi sistema v termodinamičnem ravnotežju ter mikroskopskim obnašanjem in interakcijami delcev, ki se pojavljajo znotraj preiskovanega sistema. Z uporabo statistične termodinamike lahko iz velikega števila stanj, ki jih sistem v daljšem času zavzame, razložimo makroskopske pojave in izračunamo merljive makroskopske lastnosti sistema iz karakteristik posameznih molekul. Pozornost je usmerjena na statistično ravnotežje. To ne pomeni, da so se delci prenehali gibati (mehansko ravnotežje), marveč da se ansambel ne spreminja. Statistična termodinamika nam torej priskrbi matematične zveze med merljivimi količinami za sisteme v ravnotežju, ne pove pa ničesar o trenutnih vrednostih;
pozna zgolj povprečja. V primeru, da bi poznali položaje in (vektorje) hitrosti vseh delcev v nekem trenutku, bi lahko z znanjem klasične in kvantne mehanike predvideli obnašanje sistema v poljubnem času [11]. Toda to ni mogoče iz dveh razlogov. Prvi razlog je, da hkrati ne moremo poznati dejanskih položajev in (vektorjev) hitrosti delcev – Heisenbergov princip nedoločljivosti. W. Heisenberg je namreč pokazal, da ne moremo npr. hkrati izmeriti natančne hitrosti in položaja delca v določeni smeri. To obliko principa nedoločljivosti zapišemo kot ∆𝑥 ∙ ∆𝑝𝑥 ≥ ℎ, kjer je ℎ Planckova konstanta, ∆𝑥 napaka pri meritvi koordinate 𝑥 in ∆𝑝𝑥 napaka pri meritvi komponente gibalne količine 𝑝𝑥 [12]. Drugi razlog pa je ta, da imamo preprosto preveliko število delcev in število stanj za posamezne delce, da bi lahko iz znanih stanj določili makroskopske lastnosti posameznih molekul [11]. To so torej zadostni razlogi, da se poslužimo statističnih metod, ki omogočajo napoved najbolj verjetnega obnašanja velikega števila delcev (ansambla) [11].
Omenili smo, da s statistično termodinamiko ne moremo poznati trenutnih vrednosti, zato razne fizikalne količine podamo kot časovna povprečja. Za primer vzemimo mero za povprečno hitrost delcev – temperaturo (𝑇). Merilnik temperature bo kot rezultat pokazal povprečno vrednost, ki se za celotni sistem zdi konstantna, v resnici pa temperatura majhnega števila delcev sistema niha (fluktuira) [11]. Časovno povprečje mikroskopskih konfiguracij sistema zapišemo kot:
〈𝐴〉𝑡 = lim
𝜏→∞
1
𝜏∫𝜏 𝐴 {𝑝⃑⃑⃑⃑⃑ (𝑡), 𝑟𝑁 ⃑⃑⃑⃑ (𝑡)} 𝑑𝑡𝑁
𝑡=0 (3.1)
𝑝𝑁
⃑⃑⃑⃑⃑ (𝑡) ≡ (𝑝⃑⃑⃑ , 𝑝1 ⃑⃑⃑⃑ , … , 𝑝2 ⃑⃑⃑⃑ ); 𝑝𝑁 ⃑⃑⃑ = (𝑝𝑖 𝑥, 𝑝𝑦, 𝑝𝑧)
𝑖
𝑟𝑁
⃑⃑⃑⃑ (𝑡) ≡ (𝑟⃑⃑⃑ , 𝑟1 ⃑⃑⃑ , … , 𝑟2 ⃑⃑⃑⃑ ); 𝑟𝑁 ⃑⃑ = (𝑟𝑖 𝑥, 𝑟𝑦, 𝑟𝑧)
𝑖
kjer je 𝜏 čas opazovanja sistema, 𝑝⃑⃑⃑⃑⃑ 𝑁 set gibalnih količin vseh 𝑁 delcev v sistemu in 𝑟⃑⃑⃑⃑ 𝑁 set vseh prostorskih koordinat [10].
Monte Carlo simulacije mehkih delcev Jaka Prelog
8
Časovno povprečje 〈𝐴〉𝑡 lahko izračunamo s pomočjo računalniške simulacije molekulske dinamike. Slednja se od simulacij Monte Carlo razlikuje po načinu generiranja reprezentativnih konfiguracij sistema. V primeru računalniške simulacije molekulske dinamike poskušamo ponazoriti gibanje in interagiranje molekul v nekem časovnem obdobju. Delcem določimo njihovo premikanje s silami, ki bodo delovale nanje kot to določajo klasične Newtonove enačbe gibanja. Ta metoda je v celoti deterministična in ima časovno resolucijo, vendar pa numeričen problem limite v neskončnosti ostane [10, 13]. V primeru računalniške simulacije Monte Carlo pa gre v bistvu za stohastične procese, ki temeljijo na naključnem premikanju naključnih delcev. Simulacije Monte Carlo nimajo časovne resolucije, saj dolžine simulacije ne merimo v številu časovnih korakov, kot to počnemo v računalniški simulaciji molekulske dinamike, marveč v ciklih, iz katerih pa ne moremo neposredno izračunati časovnega povprečja 〈𝐴〉𝑡. To nam omogoči šele 1. postulat statistične termodinamike, s pomočjo katerega lahko iz ansambelskega povprečja dobimo časovno povprečje [13].
3.1 Postulata statistične termodinamike
Prvi postulat statistične termodinamike je v literaturi znan pod imenom ergodijska hipoteza. Pove nam, da je v primeru, če je sistem ergodijski, časovno povprečje poljubne količine enako ansambelskemu povprečju te količine [11, 14]:
〈𝐴〉𝑡 = 〈𝐴〉𝑒 (3.2)
V tem primeru ansambelsko povprečje neke količine 𝐴 lahko zapišemo kot vsoto vseh stanj 𝑖, ki jih dani sistem lahko zavzame, pomnoženih z verjetnostjo 𝑃𝑖 za to stanje [13]:
〈𝐴〉𝑒 = ∑ 𝐴(𝑟 )
𝑖
∙ 𝑃𝑖 (3.3)
Časovno povprečje neke količine v sistemu je pravzaprav posledica prehajanj tega sistema preko različnih kvantnih stanj. Za izračun povprečja potrebujemo verjetnosti pojavljanja kvantnih stanj 𝑖 ter vrednost količine 𝐴 v posameznih kvantnih stanjih. Da ergodijska hipoteza drži, kažejo rezultati simulacij, ki so v primeru računalniške simulacije Monte Carlo ali simulacije molekulske dinamike povsem enaki [11, 14].
Drugi postulat statistične termodinamike pravi tako: Kadar velja 𝐸𝑗 = 𝐸𝑘, potem velja tudi 𝑃𝑗 = 𝑃𝑘, kar poznamo kot osnovno načelo enakih apriornih verjetnosti (angl.
principle of equal a priori probabilities). Z drugimi besedami, če imata dve mikroskopski stanji enako energijo, sta tudi enako verjetni. Pri nekih konstantnih pogojih je verjetnost za posamezno kvantno stanje odvisna le od njegove energije [11, 15]. Naj se še tako trudimo, ta dva postulata ne moremo matematično dokazati, lahko pa na podlagi eksperimentov potrdimo, da ne poznamo primera, ko bi bila prekršena [15].
9
3.2 Ansambli
Ansambel je skupek med seboj neodvisnih sistemov. V vsakem ansamblu imamo nekaj neodvisnih, konstantnih količin, ki so potrebne za opis makroskopskega termodinamskega stanja vseh sistemov v ansamblu, in količine, ki se spreminjajo (fluktuirajo) [11]. Sistemi v ansamblu se razlikujejo na molekularnem nivoju (različna klasična ali kvantna stanja), njihovo termodinamično stanje pa je enako [16].
Replike (slike) sistema se razlikujejo po razporeditvi delcev oz. gibalnih količinah.
Ansambelsko povprečje po vseh stanjih zapišemo z enačbo (3.4)3:
〈𝐴〉𝑒 = ∬ 𝐴 (𝑝⃑⃑⃑⃑⃑ , 𝑟𝑁 ⃑⃑⃑⃑ ) 𝜌 (𝑝𝑁 ⃑⃑⃑⃑⃑ , 𝑟𝑁 ⃑⃑⃑⃑ ) 𝑑𝑝𝑁 ⃑⃑⃑⃑⃑ 𝑑𝑟𝑁 ⃑⃑⃑⃑ 𝑁 (3.4) Integracija teče po vseh koordinatah položajev in vseh gibalnih količinah delcev. V primeru kanoničnega (izohorno-izotermnega) sistema, kjer so temperatura (𝑇), prostornina (𝑉) in število delcev (𝑁) konstantni, verjetnostno gostoto podamo z Boltzmannovo porazdelitvijo:
𝜌 (𝑝⃑⃑⃑⃑⃑ , 𝑟𝑁 ⃑⃑⃑⃑ ) =𝑁 𝑒−
ℋ(𝑝⃑⃑⃑⃑⃑⃑ ,𝑟𝑁⃑⃑⃑⃑⃑ )𝑁 𝑘𝐵𝑇
𝑄𝑁𝑉𝑇
(3.5)
𝑄𝑁𝑉𝑇 = 1 𝑁!
1
ℎ3𝑁∬ 𝑒−
ℋ(𝑝⃑⃑⃑⃑⃑⃑ ,𝑟𝑁⃑⃑⃑⃑⃑ )𝑁
𝑘𝐵𝑇 𝑑𝑝⃑⃑⃑⃑⃑ 𝑑𝑟𝑁 ⃑⃑⃑⃑ 𝑁 (3.6) kjer je 𝑘𝐵 Boltzmannova konstanta, ℋ Hamiltonian sistem, 𝑄𝑁𝑉𝑇 pa statistična vsota (particijska funkcija) pri konstantni temperaturi (𝑇), prostornini (𝑉) in številu delcev (𝑁).
Hamiltonian predstavlja vsoto kinetične energije, ki je odvisna od gibalnih količin delcev, in potencialne energije, ki je odvisna od koordinat položajev delcev [10, 12]:
ℋ(𝑝 , 𝑟 ) = 𝐾(𝑝 ) + 𝑈(𝑟 ) (3.7) Za opis makroskopskega stanja zadostuje le nekaj spremenljivk, ki stanje v celoti opišejo in jih lahko neodvisno spreminjamo. Glede na željene konstantne pogoje ločimo več vrst ansamblov, ki so zbrani v preglednici (3.1). Pravi ansambel mora ustrezati načinu priprave in karakterizaciji sistema – z drugimi besedami, ansambel mora odražati znanje o tem sistemu.
3 Enačba (3.4) je zaradi preglednejšega zapisa poenostavljena (6𝑁-kratni integral je zapisan kot dvojni integral). Zavedati se moramo, da upravljamo s šestimi definiranimi koordinatami, in sicer s 3𝑁 koordinatami gibalne količine in 3𝑁 koordinatami položajev.
Monte Carlo simulacije mehkih delcev Jaka Prelog
10
1. Mikrokanonični ansambel (𝑁, 𝑉, 𝐸) =konst.
Izolirani oz. mikrokanonični ansambel je najpreprostejši možni sistem, v katerem so vsa mikroskopska stanja enako verjetna (imajo enako, predpisano energijo).
Ansambel je zaprt in izoliran (z okolico ne izmenjuje snovi, niti energije). Vsi sistemi v ansamblu imajo enako prostornino (𝑉), število delcev (𝑁) in energijo (𝐸) [11, 15].
2. Kanonični ansambel (𝑁, 𝑉, 𝑇) =konst.
Gre za zaprt sistem, ki z okolico izmenjuje energijo v obliki toplote. Posamezni sistemi ansambla so med seboj neodvisni. Kanonični ansambel vsebuje stanja z različno energijo, a enako sestavo. Različna stanja v ansamblu imajo različne verjetnosti, odvisno od njihove celotne energije [11].
3. Velekanonični ansambel (𝜇, 𝑉, 𝑇) =konst.
Ta ansambel ima konstantno prostornino (𝑉), temperaturo (𝑇) in kemijski potencial (𝜇). Od zgornjih dveh se razlikuje v tem, da je odprt, torej z okolico lahko izmenjuje snov in energijo v obliki toplote. Posamezni sistemi v ansamblu lahko imajo različno število molekul in različno energijo, zato imamu tu dodaten opravek s fluktuacijami in s številom delcev. Različni sistemi lahko imajo različno število delcev, to pa pomeni, da so različni tudi razpoložljivi energijski nivoji.
Različna stanja v ansamblu imajo različne verjetnosti, odvisno od njihove celotne energije in skupnega števila delcev. Tak ansambel je primeren za študij adsorpcije [11, 15].
4. Izotermno-izobarni ansambel (𝑁, 𝑝, 𝑇) =konst.
V termodinamski limiti so vsi ansambli termodinamsko ekvivalentni, zato si lahko zamislimo še druge vrste ansamblov. Edini pogoj pri tem je, da izberemo spremenljivke, ki v celoti opišejo stanje sistema, pri tem pa mora biti vsaj ena spremenljivka odvisna od velikosti, saj v nasprotnem primeru ne moremo določiti velikosti sistema. Ker v laboratoriju često raziskujemo sisteme pri konstantnem tlaku, lahko konstruiramo ansambel, s toplotno prevodnimi in gibljivimi stenami.
Ker ima tak sistem možnost spreminjanja prostornine, se vzdržuje konstantni tlak.
Tak ansambel je izotermno-izobarni, v katerem imajo različna stanja, različne verjetnosti, odvisno od njihove celotne energije in prostornine [11].
Za potrebe te diplomske naloge smo simulacije izvajali v kanoničnem (izohorno- izotermnem) – 𝑁𝑉𝑇 ansamblu in izotermno-izobarnem – 𝑁𝑝𝑇 ansamblu. V naslednjih poglavjih bomo značilnosti teh dveh ansamblov pogledali pobližje.
11
Preglednica 3.1: Najpogostejši ansambli v statistični termodinamiki. Prirejeno po: [17].
Ansambel Konstantne količine
Fluktuirajoče (spremenljive) količine
Korespondenčni termodinamski sistem
Mikrokanonični
(𝑁𝑉𝐸) Število delcev (𝑁) Prostornina (𝑉) Notranja energija (𝑈)
Energija
kvantnega stanja vsakega delca (𝜀𝑖)
Izolirani
Kanonični
(𝑁𝑉𝑇) Število delcev (𝑁) Prostornina (𝑉) Temperatura (𝑇)
Notranja energija
(𝑈𝑖) Zaprti
Velekanonični
(𝜇𝑉𝑇) Kemijski potencial (𝜇) Prostornina (𝑉)
Temperatura (𝑇)
Število delcev (𝑁𝑖) Notranja energija (𝑈𝑖)
Odprti
Izotermno-izobarni
(𝑁𝑝𝑇) Število delcev (𝑁)
Tlak (𝑝)
Temperatura (𝑇)
Prostornina (𝑉𝑖) Notranja energija (𝑈𝑖)
/
3.2.1 Termodinamika kanoničnega ansambla
Zanima nas, katere zveze veljajo v kanoničnem sistemu, saj bomo na podlagi teh lahko določili termodinamične lastnosti sistema. Če poznamo particijsko funkcijo za kanonični ansambel 𝑄𝑁𝑉𝑇, lahko termodinamične lastnosti izračunamo iz lastnosti molekul, ki so v sistemu. Omenili smo, da je za kanonični sistem porazdelitvena gostota enaka Boltzmannovi porazdelitvi, zato verjetnost za 𝑖-to stanje (𝑃𝑖) izračunamo po enačbi (3.8):
𝑃𝑖 = 𝑒−
𝐸𝑖(𝑁,𝑉) 𝑘𝐵𝑇
𝑄𝑁𝑉𝑇
(3.8)
kjer je 𝑘𝐵 Boltzmannova konstanta, 𝑇 absolutna temperatura, 𝐸𝑖(𝑁, 𝑉) energija stanja 𝑖 pri konstantnem številu delcev 𝑁 in konstantni prostornini 𝑉 ter 𝑄𝑁𝑉𝑇 statistična vsota, ki normira verjetnost:
𝑄𝑁𝑉𝑇 = ∑ 𝑒−
𝐸𝑖(𝑁,𝑉) 𝑘𝐵𝑇 𝑖
(3.9)
Monte Carlo simulacije mehkih delcev Jaka Prelog
12
Povprečja nato računamo po enačbi (3.3) in dobimo:
〈𝐴𝑖〉𝑒 = ∑ 𝐴𝑖𝑒−
𝐸𝑖(𝑁,𝑉) 𝑘𝐵𝑇 𝑖
∑ 𝑒−
𝐸𝑖(𝑁,𝑉) 𝑘𝐵𝑇 𝑖
(3.10)
Tak način računanja povprečij je možen zgolj za mehanične količine (npr. energija), ki so definirane za vsako kvantno stanje 𝑖. Preostale termodinamske količine lahko izračunamo s pomočjo Helmholtzove proste energije, preko enačb klasične termodinamike:
𝐴𝑁𝑉𝑇 = −𝑘𝐵𝑇 ln(𝑄𝑁𝑉𝑇) (3.11) Izkaže se, da je Helmholtzova prosta energija neposredno sorazmerna statistični vsoti za kanonični sistem 𝑄𝑁𝑉𝑇. Enačba (3.11) predstavlja najpomembnejšo povezavo, saj lahko iz te enačbe in znanih termodinamičnih zvez izpeljemo enačbe med 𝑄𝑁𝑉𝑇 in preostalimi termodinamskimi količinami [11]:
1. Entropija:
𝑆𝑁𝑉𝑇= 𝑘𝐵ln(𝑄𝑁𝑉𝑇) + 𝑘𝐵𝑇 (𝜕 ln(𝑄𝑁𝑉𝑇)
𝜕𝑇 )
𝑉,𝑁
(3.12) 2. Energija:
𝐸𝑁𝑉𝑇 = 𝑘𝐵𝑇2(𝜕 ln(𝑄𝑁𝑉𝑇)
𝜕𝑇 )
𝑉,𝑁
(3.13) 3. Tlak4:
𝑝𝑁𝑉𝑇 = 𝑘𝐵𝑇 (𝜕 ln(𝑄𝑁𝑉𝑇)
𝜕𝑉 )
𝑇,𝑁
(3.14) 4. Kemijski potencial5:
𝜇𝑁𝑉𝑇 = −𝑘𝐵𝑇 (𝜕 ln(𝑄𝑁𝑉𝑇)
𝜕𝑁 )
𝑇,𝑉
(3.15)
4 Enačbo (3.14) uporabimo v primeru, če opazujemo eno kvantno stanje, zato enačbe ne smemo zamešati s povprečnim tlakom, ki je odvod proste energije na volumen (upoštevamo energijski in entropijski prispevek).
5 Enačba (3.15) je enačba kemijskega potenciala za enokomponentni sistem.
13
3.2.2 Termodinamika izotermno-izobarnega ansambla
Izotermno-izobarni ansambel je zaprt sistem pri konstantnem tlaku (𝑝) in temperaturi (𝑇) ter z okolico lahko izmenjuje energijo v obliki dela. Podobno kot pri kanoničnem ansamblu lahko izpeljemo verjetnost za 𝑖-to stanje (𝑃𝑉,𝑖), ki pa je v tem primeru odvisna od energije (𝐸𝑉,𝑖) in prostornine (𝑉):
𝑃𝑉,𝑖 = 𝑒−
𝐸𝑉,𝑖(𝑁,𝑉) 𝑘𝐵𝑇 −𝑝𝑉
𝑘𝐵𝑇
∆𝑁𝑝𝑇
(3.16)
∆𝑁𝑝𝑇= ∑ ∑ 𝑒−
𝐸𝑉,𝑖(𝑁,𝑉) 𝑘𝐵𝑇 −𝑝𝑉
𝑘𝐵𝑇 𝑖
𝑉
(3.17)
kjer je 𝑘𝐵 Boltzmannova konstanta, 𝑇 absolutna temperatura, 𝐸𝑉,𝑖(𝑁, 𝑉) energija stanja 𝑖 pri konstantnem številu delcev 𝑁 in konstantni prostornini 𝑉 ter ∆𝑁𝑝𝑇 statistična vsota.
Povprečja nato računamo po enačbi (3.3). Preostale termodinamske količine lahko izračunamo s pomočjo Gibbsove proste entalpije preko enačb klasične termodinamike:
𝐺𝑁𝑝𝑇 = −𝑘𝐵𝑇 ln(∆𝑁𝑝𝑇) (3.18) Izkaže se, da je Gibbsova prosta entalpija neposredno sorazmerna statistični vsoti izotermno-izobarnega ansambla ∆𝑁𝑝𝑇. Omenili bomo še nekaj enačb, ki v celoti opišejo izotermno-izobarni sistem [11, 15]:
1. Entropija:
𝑆𝑁𝑝𝑇= − (𝜕𝐺𝑁𝑝𝑇
𝜕𝑇 )
𝑁,𝑝
(3.19) 2. Prostornina:
𝑉𝑁𝑝𝑇 = (𝜕𝐺𝑁𝑝𝑇
𝜕𝑝 )
𝑁,𝑇
(3.20) 3. Kemijski potencial:
𝜇𝑁𝑝𝑇 = (𝜕𝐺𝑁𝑝𝑇
𝜕𝑁 )
𝑝,𝑇
(3.21)
Monte Carlo simulacije mehkih delcev Jaka Prelog
14
3.3 Zasnova računalniške simulacije Monte Carlo
Simulacijo navadno začnemo z generiranjem začetne konfiguracije sistema. Nato nadaljujemo z uravnoteženjem sistema (ekvilibracijo) in nazadnje izvedemo še produkcijski del simulacije, v katerem zajemamo karakteristike posameznih konfiguracij in računamo ustrezna termodinamična povprečja [10].
3.3.1 Generiranje začetne konfiguracije
Ko govorimo o generiranju začetne konfiguracije, si je v začetku potrebno izbrati dober model za sistem, ki ga proučujemo. Pred začetkom simulacije, moramo poznati še [10]:
1. karakteristike sistema – npr. temperaturo, tlak, koncentracijo komponent …;
2. ansambel – na podlagi tega, pri katerih pogojih bomo preučevali sistem, izberemo ustrezni ansambel (kanonični, velekanonični, mikrokanonični, izotermno- izobarni, Gibbsov …), saj je od tega odvisen izračun termodinamičnih lastnosti;
3. velikost osnovne simulacijske celice – za velikost osnovne simulacijske celice se odločimo na podlagi daljnosežnih interakcij med delci. Navadno velikost celice izberemo tako, da prvi delec, ki je postavljen na sredino celice, ne bo več čutil drugega delca, ki je zunaj celice. Na razdalji polovice roba celice je torej potencial med dvema delcema približno nič. Prav tako moramo še upoštevati, da sta število delcev v simulacijski celici in velikost simulacijske celice med seboj odvisni, in da želimo imeti v simulacijski celici izbrano koncentracijo delcev (vsaj reda velikosti 100);
4. obliko osnovne simulacijske celice – najpogosteje je to kocka (za simulacije v tridimenzionalnem sistemu) oz. kvadrat (za simulacije v dvodimenzionalnem sistemu). Koordinate delcev izberemo naključno (delce naključno razporedimo po prostoru), učinkovitejše pa je, če delce izberemo v konfiguraciji, ki je čim bližja ravnotežni. Za začetne lege delcev najpogosteje izberemo kar ploskovno centrirano kocko. Pomagamo si lahko s podatki rentgenske praškovne difrakcije ali jedrske magnetne resonance – NMR. Da se izognemo nestabilnostim, je pred samo simulacijo dobro eliminirati visokoenergijske napetosti v konfiguraciji, kar lahko storimo med simulacijo z energijsko minimizacijo sistema. Pazimo le, da energijsko minimizacijo storimo pred uravnoteženjem sistema.
Generator naključnih števil
Ključno vlogo v računalniški simulaciji Monte Carlo ima generator naključnih števil.
Pomembno je, da ta deluje pravilno. V ta namen imamo različne teste, s katerimi lahko preverimo delovanje generatorja (npr. statistični χ2 test). Zavedati se moramo, da v resnici naključna števila generiramo z matematičnimi algoritmi, torej v osnovi ta niso naključna;
15
rečemo, da so psevdonaključna. Psevdonaključna števila s pridom uporabljamo v simulacijah, saj so praktična zaradi svoje obnovljivosti in hitrega tvorjenja. Števila si sledijo po točno določenem algoritmu, niz pa se po preteku določenega časa ponovi.
Generator se zažene iz poljubnega začetnega stanja s pomočjo semena, tj. naključne začetne vrednosti. Tako bo zaporedje pri enakem semenu zmeraj enako. Psevdonaključna števila morajo zadoščati naključni statistiki: biti morajo enakomerno razporejena, neodvisna eden od drugega, imeti morajo ustrezno ponovljivost; vsota dveh zaporednih psevdonaključnih števil pa mora biti z enako verjetnostjo soda kot liha. Merila, ki jih izpolnjuje dober generator naključnih števil, so [10]:
1. učinkovitost – generator je hiter;
2. uniformnost – vsa števila so enako verjetna;
3. neodvisnost – iz trenutnega števila ne gre napovedati naslednjega števila;
4. dolga perioda niza – zaporedje se ponovi po preteku določenega časa in po velikem številu generiranih psevdonaključnih števil.
3.3.2 Problem makroskopskega števila delcev
Kljub visoki zmogljivosti računalnikov še vedno s težavo modeliramo sisteme z velikim številom delcev. Povsem nemogoče pa je v računalniški simulaciji ponazoriti sistem z makroskopskim številom delcev, kot je recimo sistem z Avogadrovim številom delcev reda velikosti 1023. Takšne naloge ne zmore noben računalnik, saj smo omejeni s hitrostjo izvajanja programa, zmogljivostjo procesorjev in velikostjo pomnilnika. Posledično so delci v simulacijski celici izpostavljeni površinskim pojavom. To pomeni, da na delce, ki so bližje površini simulacijske celice, delujejo sile, ki so drugačne od tistih, ki delujejo na preostale delce v notranjosti. Problem pride toliko bolj do izraza pri manjši prostornini sistema, saj je takrat večji delež delcev izpostavljen vplivom površine. Za ponazoritev vpliva površine vzemimo 1 L vode pri sobni temperaturi, kjer je približno 3,3∙1025 molekul vode. Delež molekul, ki so v stiku s površino, je približno 6∙10–7, medtem ko so v simulacijski celici pri povsem enakih pogojih, vendar zgolj s 1000 molekulami vode, v stiku s površino pravzaprav vse [10, 15]. Obstajajo načini, ki lahko ta problem zmanjšajo, ne pa tudi povsem odpravijo. V nadaljevanju bomo predstavili nekaj najpogostejših načinov.
Periodični robni pogoji
V kolikor z računalniško simulacijo proučujemo vpliv površine, potem v tej točki s simulacijo nimamo težav, v nasprotnem primeru si pomagamo s periodičnimi robnimi pogoji (angl. Periodic Boundary Condition, PBC). Pri tem simulacijsko celico neskončno krat repliciramo v vseh smereh po prostoru, kot je to nazorno prikazano na sliki (3.1).
Monte Carlo simulacije mehkih delcev Jaka Prelog
16
S tem zagotovimo, da bodo na delec delovale enake sile, kot bi v realnem makroskopskem sistemu [10, 13, 15]. Če se delec denimo v izhodiščni simulacijski celici premakne, se hkrati z njim premaknejo tudi vsi kloni tega delca v preostalih periodičnih simulacijskih celicah. V primeru, da delec zapusti celico, njegovo mesto zavzame replika iz sosednje celice. Tako zaradi periodičnih robnih pogojev število delcev v vsaki simulacijski celici ostane nespremenjeno, sočasno pa eliminiramo površinske efekte [10].
Približek najbližje slike
Navadno periodične robne pogoje uporabljamo v povezavi s približkom najbližje slike (angl. minimum image convention). Tak način je najbolj primeren za simulacije Monte Carlo. Pogosto namreč moramo izračunati energijo celotnega sistema. Ker pa v praksi to ni mogoče zaradi neskončnega sistema, kjer bi morali upoštevati interakcije vseh delcev z vsemi replikami, s približkom najbližje slike, upoštevamo zgolj tiste interakcije replike delca (od vseh možnih), ki je izbranemu delcu najbližje. Vse ostale interakcije lahko zanemarimo le v primeru, če so dovolj kratkosežne. Pri računanju interakcij med dvema delcema v kartezičnih koordinatah, kjer velja:
∆𝑥 = |𝑥𝑖 − 𝑥𝑗|,
(3.22)
∆𝑦 = |𝑦𝑖− 𝑦𝑗| in
∆𝑧 = |𝑧𝑖 − 𝑧𝑗|,
razdaljo 𝑟𝑖𝑗 med delcema 𝑖 in 𝑗 izračunamo po enačbi:
𝑟𝑖𝑗 = √(∆𝑥2+ ∆𝑦2+ ∆𝑧2 (3.23) V kolikor je recimo razlika ∆𝑥 med koordinatama delcev večja od polovice dimenzije simulacijske celice 𝐿𝑥
2, v izračunu 𝑟𝑖𝑗, uporabimo 𝐿𝑥
2 − ∆𝑥:
∆𝑥𝑛𝑜𝑣 = {
∆𝑥 č𝑒 ∆𝑥 ≤𝐿𝑥 2 ,
∆𝑥 −𝐿𝑥
2 𝑠𝑖𝑐𝑒𝑟.
(3.24)
Podobno storimo tudi za vse preostale koordinate [20]. Potencial je navadno sferično simetričen, zato je približek najbližje slike možno uporabiti skupaj s sferičnim
»rezanjem« potenciala (angl. r cutoff).
Na določeni razdalji, kjer je 𝑟𝐶 ≤ 𝐿
2 in 𝐿 predstavlja dolžino roba osnovne celice, upoštevamo zgolj tiste interakcije z delci, katerih oddaljenost od izhodiščnega delca je manjša od 𝑟𝐶. Rezanje potenciala nam lahko povzroči težave pri računalniški simulaciji, saj v točki 𝑟𝐶 potencialna energija ni več zvezna. Kot alternativo lahko eventualno uporabimo »preklopno funkcijo« polinoma, s katerim pomnožimo potencial, s čimer odpravimo nezveznosti [10].
17
Slika 3.1: Shematski prikaz periodičnih robnih pogojev skupaj z uporabo približka najbližje slike v dveh dimenzijah6. Prirejeno po [19].
3.3.3 Uravnoteženje sistema
Pri računalniški molekulski simulaciji Monte Carlo je mogoče vzorčiti le ravnotežno porazdelitev stanj. Namen uravnoteženja sistema je premakniti začetno konfiguracijo sistema v termodinamično ravnotežje, kjer se sistem ne bo sistematično spreminjal. To pomeni, da če v sistem vnesemo majhno motnjo, se bo ta vrnil nazaj v ravnotežno stanje.
Za uravnoteženje sistema navadno uporabimo Metropolisovo metodo, ki je opisana v poglavju 3.5. Če želimo določiti, kdaj sistem doseže ravnotežno stanje, med simulacijo nujno spremljamo različne termodinamične in strukturne lastnosti sistema in to odvisno od izbranih vzdrževanih konstantnih pogojev, ki smo jih izbrali pred začetkom simulacije.
Tako lahko metodo Monte Carlo uporabimo za generiranje velikega števila naključnih stanj sistema in iz teh po teoriji statistične termodinamike izračunamo ansambelska povprečja. Za vsako sprejeto stanje izračunamo zgoraj navedene termodinamične lastnosti, ki pa jih podamo kot povprečje. Povprečno notranjo energijo (𝑈), prostornino (𝑉) in število delcev (𝑁) izračunamo po enačbi (3.25) kot povprečje po konfiguracijah, generiranih tekom simulacije:
6 V treh dimenzijah imamo namesto kvadratov kocke, ki so po prostoru razporejene še po tretji dimenziji.
Monte Carlo simulacije mehkih delcev Jaka Prelog
18
〈𝐴〉 = 1
𝑀∑ 𝐴 (𝑝⃑⃑⃑⃑⃑ , 𝑟𝑁 ⃑⃑⃑⃑ )𝑁
𝑀
𝑖=1
(3.25)
kjer je 𝑀 število generiranih konfiguracij (oz. ciklov) v simulaciji.
Da med ansambelskim povprečjem 〈𝐴〉 in povprečjem simulacij res velja enakost, mora biti število tvorjenih konfiguracij (število ciklov) dovolj veliko (𝑀 → ∞). Pri simulaciji Monte Carlo so povprečne količine odvisne zgolj od razporeditve delcev, zato njihovih gibalnih količin 𝑝⃑⃑⃑⃑⃑ 𝑁 sploh ni potrebno računati [10, 19]. Enačba (3.25) se poenostavi:
〈𝐴〉 = 1
𝑀∑ 𝐴 (𝑟⃑⃑⃑⃑ )𝑁
𝑀
𝑖=1
(3.26)
3.3.4 Produkcijski del računalniške simulacije Monte Carlo
Ko sistem uravnotežimo in se preiskovane količine ustalijo okrog povprečne vrednosti, rečemo, da nihajo oz. fluktuirajo okrog povprečne vrednosti, v produkcijskem delu simulacije izračunamo termodinamične, dinamične in strukturne lastnosti izbranega sistema. Če želimo zagotoviti pravilne in natančne rezultate, ki jih bo možno primerjati z eksperimentalnimi, moramo pustiti, da simulacija teče dovolj dolgo. Za oceno zanesljivosti simulacije moramo izračune večkrat ponoviti, zato simulacijo navadno organiziramo v serijah, kjer mora vsaka serija vsebovati dovolj veliko število konfiguracij, da zadosti pogoju termodinamičnega ravnotežja. Začetna konfiguracija nove serije je enaka končni konfiguraciji predhodne serije. V splošnem je število generiranih konfiguracij znotraj posamezne serije (razporeditev, hitrost in orientacija delcev) odvisno od sistema, ki ga preučujemo, zato mora biti to število dovolj veliko, da izračunane lastnosti sistema ne bodo odvisne od števila generiranih konfiguracij. Dolžino računalniške simulacije Monte Carlo merimo v ciklih. V vsakem ciklu izvedemo premik vseh delcev (𝑁 premikov) v naši simulacijski celici. Če molekula ni sferična, izvedemo še 𝑁 rotacij, v primeru konstantnega tlaka še spremembo prostornine oz. v primeru odprtega sistema še spremembo števila delcev. V vsaki seriji izvedemo določeno število ciklov. Tekom vsake serije izračunamo termodinamična povprečja količin, kot so npr.
energija, prostornina, število delcev (enačba 3.26), tlak (enačba 3.48), toplotna kapaciteta (enačba 3.52) in temperatura. Rezultate iz posameznih ciklov povprečimo, prav tako povprečimo še rezultate iz vsake serije. Povprečje znotraj vsake serije tako predstavlja eno »meritev« in končen rezultat naše preiskovane količine [10].
19
3.4 Metropolisov algoritem
Statistično vsoto v kanoničnem ansamblu opiše enačba (3.6). Če želimo Hamiltonian (ℋ) v enačbi (3.6) zapisati ločeno kot seštevek kinetične in potencialne energije, dobimo:
𝑄𝑁𝑉𝑇 = 1 𝑁!
1 ℎ3𝑁∫ 𝑒−
|𝑝⃑⃑⃑⃑⃑⃑ |𝑁2
2𝑚𝑘𝐵𝑇𝑑𝑝⃑⃑⃑⃑⃑ ∫ 𝑒𝑁 −
𝜈(𝑟⃑⃑⃑⃑⃑ )𝑁
𝑘𝐵𝑇 𝑑𝑟⃑⃑⃑⃑ 𝑁 (3.27) Prispevka potencialne in kinetične energije med seboj nista odvisna, zato ju lahko ločimo.
Tako iz 6𝑁-kratnega integrala iz enačbe (3.6) dobimo produkt dveh 3𝑁-kratnih integralov. Prispevek kinetične energije lahko rešimo analitično, saj so gibanja posameznih delcev neodvisna v vseh smereh prostora, neodvisne pa so tudi posamezne smeri gibanja določenega delca, zato lahko 3𝑁-kratni integral razstavimo na 3𝑁 neodvisne integrale. Prispevek potencialne energije pa ni odvisen od hitrosti delcev, marveč od koordinat, tj. konfiguracije sistema, kjer so vse koordinate delcev med seboj odvisne. Tako nam ostane konfiguracijski integral 𝑍𝑁𝑉𝑇, ki pa ga ne moremo izračunati analitično [11]:
𝑍𝑁𝑉𝑇= ∫ 𝑒−
𝜈(𝑟⃑⃑⃑⃑⃑ )𝑁
𝑘𝐵𝑇 𝑑𝑟⃑⃑⃑⃑ 𝑁 (3.28)
Λ = ℎ
√2𝜋𝑚𝑘𝐵𝑇
7(3.29)
𝑄𝑁𝑉𝑇 = 𝑍𝑁𝑉𝑇
𝑁! Λ3𝑁 = 1
𝑁! Λ3𝑁∫ 𝑒−
𝜈(𝑟⃑⃑⃑⃑⃑ )𝑁
𝑘𝐵𝑇 𝑑𝑟⃑⃑⃑⃑ = 𝑄𝑁 𝑁𝑉𝑇𝑖𝑑 𝑄𝑁𝑉𝑇𝑒𝑥 (3.30) Opazimo, da so termodinamične lastnosti sestavljene iz idealnega dela 𝑄𝑁𝑉𝑇𝑖𝑑 in presežnega dela 𝑄𝑁𝑉𝑇𝑒𝑥 , kjer je presežni del posledica interakcij med delci. Podobno kot particijsko funkcijo, lahko tudi ostale termodinamske funkcije razdelimo na idealni in presežni del [11]. Enačba (3.31) ponazarja to lastnost na primeru Helmholtzove proste energije za kanonični ansambel:
𝐴𝑁𝑉𝑇 = −𝑘𝐵𝑇 ln(𝑄𝑁𝑉𝑇) = −𝑘𝐵𝑇(ln(𝑄𝑁𝑉𝑇𝑖𝑑 ) + ln(𝑄𝑁𝑉𝑇𝑒𝑥 )) = 𝐴𝑖𝑑𝑁𝑉𝑇+ 𝐴𝑁𝑉𝑇𝑒𝑥 (3.31) Presežno (potencialno) energijo izračunamo po enačbi (3.32) [10]:
〈𝜈 (𝑟⃑⃑⃑⃑ )〉 =𝑁 ∫ 𝜈 (𝑟⃑⃑⃑⃑ ) 𝑒𝑁 −
𝜈(𝑟⃑⃑⃑⃑⃑ )𝑁 𝑘𝐵𝑇 𝑑𝑟⃑⃑⃑⃑ 𝑁 𝑍𝑁𝑉𝑇
(3.32)
7 Enačbo (3.29) dobimo z integracijo po gibalni količini.
Monte Carlo simulacije mehkih delcev Jaka Prelog
20
Ker tega 3𝑁-kratnega integrala ne moremo analitično rešiti, se poslužimo uporabe numeričnih metod. Kot alternativo izberemo integracijo Monte Carlo, kjer za vsako naključno generirano konfiguracijo izračunamo vrednost funkcije v števcu enačbe (3.32) in imenovalcu (𝑍𝑁𝑉𝑇) ter jih seštevamo [4]:
〈𝜈 (𝑟⃑⃑⃑⃑ )〉 =𝑁 ∑ 𝜈𝑖(𝑟⃑⃑⃑⃑ ) 𝑒𝑁 −
𝜈𝑖(𝑟⃑⃑⃑⃑⃑ )𝑁 𝑘𝐵𝑇 𝑁𝑝𝑜𝑠𝑘𝑢𝑠𝑜𝑣
𝑖=1
∑ 𝑒−
𝜈𝑖(𝑟⃑⃑⃑⃑⃑ )𝑁 𝑘𝐵𝑇 𝑁𝑝𝑜𝑠𝑘𝑢𝑠𝑜𝑣
𝑖=1
(3.33)
Žal pa takšen račun ni uspešen, saj vsoti v končnem številu poskusov ne bosta konvergirali [10]. Popolnoma naključno namreč upoštevamo zelo malo verjetne konfiguracije, ki k povprečju ne prispevajo bistveno. Ta problem rešimo tako, da vzorčimo le vplivna (verjetnejša) stanja konfiguracijskega prostora, ki bistveno prispevajo k povprečnim vrednostim. Metodo je leta 1953 predlagal N. Metropolis skupaj s sodelavci [18]. Metropolisov algoritem za generiranje konfiguracij sistema uporablja Boltzamnnovo verjetnost, ki jih nato v vsoti enakovredno upoštevamo [10]. V algoritmu tvorimo t. i. Markovsko verigo stanj (angl. Markov Chain Monte Carlo methods, MCMC), kjer je novo stanje sistema zmeraj odvisno le od trenutnega stanja in ne od prejšnjih [21]. Omenili smo, da v vsakem ciklu izvedemo premik vseh delcev (𝑁 premikov). Najmanjšo enoto Markovske verige, kjer naključno izbrani delec premaknemo ali zarotiramo za naključen korak, imenujmo premik. Vsak premik je sestavljen iz treh faz [13, 15]:
1. Izbira: Izberemo delec 𝑖 s konfiguracijo 𝑟⃑⃑⃑⃑ 𝑖𝑁, ki ga bomo premaknili. Ta izbira ne vpliva na rezultat, uporabimo pa jo samo na začetku. Imenujmo jo stara konfiguracija (𝑠).
2. Premik: Tvorimo naključen premik delca 𝑖 in novo konfiguracijo označimo z (𝑛).
Potem izračunamo novo energijo sistema 𝜈(𝑛) in jo primerjamo s staro energijo sistema 𝜈(𝑠). Po enačbi (3.34) izračunamo spremembo celotne energije sistema, ki jo povzroči premik:
∆𝜈 = 𝜈(𝑛) − 𝜈(𝑠) (3.34)
Odločitev: Odločiti se moramo, ali novo stanje sprejmemo ali ne. Če se energija celotnega sistema zmanjša, potem novo stanje vedno sprejmemo (stanje z nižjo energijo je bolj ugodno), v nasprotnem primeru pa računalnik generira novo naključno število 𝜉 med 0 in 1, ki ga primerja z Boltzmannovim faktorjem (𝑓), kot to opisuje enačba (3.35):
𝑓 = {1, ∆𝜈 ≤ 0 𝑒−
∆𝜈
𝑘𝐵𝑇, ∆𝜈 > 0 (3.35)
21
Če je 𝜉 < 𝑓, novo konfiguracijo delca sprejmemo, v nasprotnem pa obdržimo staro konfiguracijo delca. Celoten algoritem je prikazan na sliki (3.2). Postopek ponavljamo, dokler se povprečje merjene termodinamične količine ne ustali.
Slika 3.2: Shematski prikaz algoritma računanja z Metropolisovo Monte Carlo metodo.
Prirejeno po [14].
Povprečno energijo, kakor tudi druge termodinamične lastnosti, izračunamo kot aritmetično povprečje po vseh konfiguracijah po enačbi (3.36) [10]:
〈𝜈 (𝑟⃑⃑⃑⃑ )〉 =𝑁 1
𝑁𝑝𝑜𝑠𝑘𝑢𝑠𝑜𝑣 ∑ 𝜈𝑖
𝑁𝑝𝑜𝑠𝑘𝑢𝑠𝑜𝑣
𝑖=1
(𝑟⃑⃑⃑⃑ ) 𝑁 (3.36)
kjer upoštevamo energije vseh sprejetih konfiguracij (štejemo tudi stara stanja, v primeru, da sistem ostane v starem stanju).
Pri vzorčenju vplivnega dela konfiguracijskega prostora Metropolisov algoritem išče premike, ki stanje najbolj izboljšajo, izogiba pa se lokalnim ekstremom. Prehod iz starega v novo stanje je dan s prehodno verjetnostjo 𝑃(𝑠 → 𝑛). Prehodna verjetnost mora v Markovski verigi izpolnjevati določene pogoje [13, 15]:
𝑟𝑠𝑡𝑎𝑟
⃑⃑⃑⃑⃑⃑⃑⃑ = 𝑟⃑⃑⃑⃑ 𝑖𝑛 𝜈(𝑠) = 𝜈 (𝑟𝑖𝑁 ⃑⃑⃑⃑ ) 𝑖𝑁
𝑟𝑛𝑜𝑣
⃑⃑⃑⃑⃑⃑⃑ = 𝑟⃑⃑⃑⃑ + 𝛥 𝑖𝑛 𝜈(𝑛) = 𝜈(𝑟𝑖𝑁 ⃑⃑⃑⃑⃑⃑⃑ ) 𝑛𝑜𝑣
𝑓 =𝑒−
𝜈(𝑛) 𝑘𝐵𝑇
𝑒−
𝜈(𝑠) 𝑘𝐵𝑇
Trenutna konfiguracija
Testna konfiguracija
Izračun Boltzmannovega faktorja
Monte Carlo simulacije mehkih delcev Jaka Prelog
22
1. Stacionarnost procesa (angl. detailed balance condition) in princip mikroskopske reverzibilnosti: Za sistem v ravnotežju mora veljati, da je premik iz starega v novo stanje enako verjeten kot premik iz novega v staro stanje, pri čemer sistem mora ostati v ravnotežju, ko ga doseže:
𝑁(𝑠)𝑃(𝑠 → 𝑛) = 𝑁(𝑛)𝑃(𝑛 → 𝑠) (3.37) kjer je 𝑁(𝑠) verjetnost stare konfiguracije in 𝑁(𝑛) verjetnost nove konfiguracije.
Zvezo med verjetnostjo in energijo opisuje enačba (3.38):
𝑃(𝑠 → 𝑛)
𝑃(𝑛 → 𝑠)=𝑁(𝑛) 𝑁(𝑠) = 𝑒−
1
𝑘𝐵𝑇(𝜈(𝑛)−𝜈(𝑠))
(3.38) 2. Ergodijski pogoj: V končnem številu premikov mora biti dosegljiva vsaka točka
faznega prostora, kar s simboli zapišemo kot 𝑃(𝑖 → 𝑗) > 0.
Nove konfiguracije določimo tako, da trenutni koordinati naključno izbranega delca prištejemo maksimalni premik 𝑑𝑟𝑚𝑎𝑥, ki ga pomnožimo z naključnim številom 𝜉 med 0 in 1. Tako delec z enako verjetnostjo premikamo v katerokoli smer (pozitivno ali negativno). Pazimo le, da maksimalni premik določimo tako, da bo vsaj polovica vseh premikov sprejetih. Če je premik premajhen, se sistemi med seboj ne bodo signifikantno razlikovali, saj bo večina premikov sprejetih; v kolikor pa izberemo prevelik premik, bi v sistemu prihajalo do prekrivanj med delci in premiki večinoma sploh ne bi bili sprejeti [13, 19].
𝑥𝑖𝑛𝑜𝑣 = 𝑥𝑖𝑠𝑡𝑎𝑟+ (𝜉 −1
2) 𝑑𝑟𝑚𝑎𝑥
(3.39) 𝑦𝑖𝑛𝑜𝑣 = 𝑦𝑖𝑠𝑡𝑎𝑟 + (𝜉 −1
2) 𝑑𝑟𝑚𝑎𝑥 𝑧𝑖𝑛𝑜𝑣 = 𝑧𝑖𝑠𝑡𝑎𝑟+ (𝜉 −1
2) 𝑑𝑟𝑚𝑎𝑥
Če želimo v sistemu namesto konstantne prostornine konstanten tlak (izotermno-izobarni ansambel), moramo omogočiti širjenje oz. krčenje simulacijske celice (spremembo prostornine) za 𝛿𝑉𝑚𝑎𝑥 [10]:
𝑉𝑛𝑜𝑣 = 𝑉𝑠𝑡𝑎𝑟 + 𝜉𝛿𝑉𝑚𝑎𝑥 (3.40) kjer je 𝜉 naključno število med −1 in 1. S spremembo prostornine simulacijske celice se spremenijo tudi vse koordinate delcev. Da se v sistemu ohranja konstantni tlak, vse koordinate pomnožimo s faktorjem (𝑉𝑛𝑜𝑣⁄𝑉𝑠𝑡𝑎𝑟)13, novo konfiguracijo pa primerjamo s faktorjem (𝑓𝑉):
𝑓𝑉 = 𝑒−
(𝜈(𝑛)−𝜈(𝑠))+𝑝(𝑉𝑛𝑜𝑣−𝑉𝑠𝑡𝑎𝑟)−𝑁𝑘𝐵𝑇 ln𝑉𝑛𝑜𝑣 𝑉𝑠𝑡𝑎𝑟
𝑘𝐵𝑇 (3.41)
23
3.5 Termodinamične in strukturne lastnosti
3.5.1 Reducirane količine
Med računalniško simulacijo pogosto naletimo na prevelike ali premajhne rezultate, ki za shranjevanje v pomnilniku po nepotrebnem zavzamejo prostor, zato se pogosto odločimo, da fizikalne količine pretvorimo v reducirane fizikalne količine (v nadaljevanju jih bomo označili z »∗«). Prednost reduciranih količin se izkaže predvsem v praksi, saj rezultate lažje ocenimo in si tako ustvarimo dober občutek o obnašanju sistemov. Izberemo si dva nova parametra za energijo (𝜀) in dolžino (𝜎) ter z njuno pomočjo formuliramo naslednje količine [13]:
1. Reducirana notranja energija:
𝑈∗ =𝑈
𝜀 (3.42)
2. Reduciran tlak:
𝑝∗ =𝑝𝜎3
𝜀 (3.43)
3. Reducirana gostota:
𝜌∗ = 𝜌𝜎3 (3.44)
4. Reducirana prostornina:
𝑉∗ = 𝑉
𝜎3 (3.45)
5. Reducirana temperatura:
𝑇∗ =𝑘𝐵𝑇
𝜀 (3.46)
kjer je 𝑘𝐵 Boltzmannova konstanta.
3.5.2 Izračun termodinamičnih količin
V tem poglavju bomo predstavili enačbe, ki smo jih uporabili pri računalniški simulaciji Monte Carlo za izračun želenih termodinamičnih količin.
1. Notranja energija:
Povprečno notranjo energijo (𝑈) izračunamo kot povprečje po vseh konfiguracijah, generiranih tekom simulacije:
𝑈 = 〈𝐸〉 = 1 𝑀∑ 𝐸𝑖
𝑀
𝑖=1
(3.47)
