Molekulové modelování a simulace c Jiří Kolafa ([email protected]), 17. listopadu 2015, Ústav Fyzikální chemie, VŠCHT Praha Tento text obsahuje m.j. upravené části skript Skripta Fyzikální chemie – bakalářský a magisterský kurz, J. Novák a kol., VŠCHT Praha 2008, 1. vydání, ISBN 978-80-7080-675-3, a Úvod do molekulárních simulací – Metody Monte Carlo a molekulární dynamiky, I. Nezbeda, J. Kolafa, M. Kotrla, Univerzita Karlova, Praha 2002 Chyby a nedostatky hlaste prosím autorovi. Obsah 1 Úvod 3 2 Ideální plyn v kinetické teorii 4 2.1 Stavová rovnice ideálního plynu ......................... 4 2.2 Ekvipartiční princip ................................ 6 2.3 Energetická rovnice ................................ 6 3 Statistická termodynamika 7 3.1 Pravděpodobnost ................................. 7 3.1.1 Mikrokanonický soubor .......................... 7 3.1.2 Kanonický soubor ............................. 9 3.2 Termodynamika v kanonickém souboru ..................... 13 3.2.1 Helmholtzova energie ........................... 14 3.2.2 Ideální plyn ................................ 15 3.2.3 Klasická Helmholtzova energie ...................... 17 3.3 Pravděpodobnostní rozložení molekulárních rychlostí ............. 18 4 Síly mezi molekulami 19 4.1 Jednoatomové molekuly ............................. 20 4.1.1 Dva neutrální atomy ........................... 20 4.1.2 Více atomů ................................ 22 4.2 Ionty ........................................ 22 4.3 Polarizovatelnost ................................. 22 4.4 Víceatomové molekuly .............................. 23 4.4.1 Nevazebné síly .............................. 23 4.4.2 Vazebné síly ................................ 25 4.5 Vnější síly ..................................... 26 5 Klasické mřížkové modely 27 5.1 Isingův model ................................... 27 5.2 Mřížkový plyn ................................... 27 5.3 Binární slitina ................................... 28 5.4 Model polymeru .................................. 28 1
Tento text obsahuje m.j. upravené části skript Skripta Fyzikální chemie – bakalářský amagisterský kurz, J. Novák a kol., VŠCHT Praha 2008, 1. vydání, ISBN 978-80-7080-675-3,
a Úvod do molekulárních simulací – Metody Monte Carlo a molekulární dynamiky,I. Nezbeda, J. Kolafa, M. Kotrla, Univerzita Karlova, Praha 2002
Chemie, vzešlá ze středověké alchymie a nejrůznějších empirických poznatků o látkách, seještě v osmnáctém století zabývala pouze „sléváním vodičekÿ. Přelom v chápání složeníhmoty nastal v roce 1803, kdy J. Dalton publikoval moderní verzi atomové hypotézy. Tentomodel vedl nejen k rozvoji samotné chemie (Mendělejevova periodická tabulka na konci 19.století), ale stal se i základem Boltzmannovy kinetické teorie a statistické termodynamiky,která úspěšně vysvětlila (nejen) chování plynů. Všechny důkazy o existenci atomů a molekulvšak zůstávaly nepřímé – atomy nikdo neviděl1, a tak ještě začátkem 20. stol. mohl E. Machpovažovat atomovou teorii za pouhou neopodstatněnou hypotézu.
Na začátku 20. stol. se však podařilo rozložit atom na menší částečky – jádra a elektrony– a později jádra na protony a neutrony a ty dále na kvarky a gluony (to je však již doménoufyziky a vysokých energií). Tyto experimentální poznatky (spolu s problémem záření černéhotělesa či fotoefektu) vedly k formulování kvantové teorie. Její aplikace na systémy složenéz elektronů a jader, tzv. kvantová chemie, umí vysvětlit existenci, stabilitu a vlastnostiatomů a molekul, průběh chemických reakcí, spektroskopické vlastnosti a mnoho dalšíchjevů. Shoda teorie s experimentem dosahuje u jednoduchých systémů (které umíme přesněspočítat) neuvěřitelných dvanácti desetinných míst2.
Vraťme se ale k fyzikální chemii. V úvodním kurzu jste se seznámili s klasickou termo-dynamikou. Principy této teorie byly formulovány axiomaticky – jako čtyři termodynamickézákony nebo „větyÿ (nultá až třetí) zobecňující experimentální zkušenosti. Pro odvození ter-modynamických vztahů nebylo podstatné, zda jsou zkoumané systémy složeny elektronů ajader nebo z atomů a molekul. Toto jednak odpovídá historickému vývoji, protože základytermodynamiky byly položeny dříve, než byla všeobecně přijata atomová hypotéza (o kvan-tové teorii ani nemluvě), jednak nám umožňuje použít termodynamiku i na systémy, které sez atomů neskládají (záření, neutronové hvězdy), i když toto vás jako chemiky trápit nebude.
Takže situace je jednoduchá. Máme kvantovou chemii, která je schopná (v principu) po-psat celou chemii. Vezmeme tedy systém, který nás zajímá, spočítáme, kolik jakých jadertam máme a kolik elektronů, nasypeme to do počítače, stiskneme tlačítko a přečteme výsle-dek. Tomu se říká výpočet ab initio, protože používáme pouze hodnoty základních přírodníchkonstant (Planckova konstanta, elementární náboj, hmotnosti částic). Pokud nám (třeba prointerpretaci spekter či chemických reakcí v plynné fázi) stačí málo atomů, jde to dobře. Po-tíž nastává, když potřebujeme studovat velké systémy (DNA, kapaliny), protože výpočetnínáročnost kvantových metod prudce vzrůstá s počtem elektronů a jader. Nicméně dnes jižumíme nasypat do počítače 100 jader kyslíku, dvakrát tolik protonů a 1000 elektronů a do-stat malinkou krychličku kapalné vody a spočítat an initio její hustotu, permitivitu apod.s přesností, která sice není závratná, ale ani směšná.
Jestliže chceme studovat systémy složené z ohromného množství atomů nebo molekul,musíme udělat mezikrok, zbavit se elektronů a jader a vytvořit model molekuly. Tenmůže být extrémně jednoduchý (hmotné body v případě ideálního plynu), jednoduchý (tvrdákulička pro vysvětlení chování plynu) či složitý (model proteinu). Jednou z možností tvorbymodelu jsou kvantové výpočty jedné či několika málo molekul. Alternativou (jedinou možnou
1To se změnilo až nedávno. Zařízení jako „Atomic force microscopeÿ umí detekovat a zobrazit jednotlivéatomy.2Toto není samozřejmé. Např. teorie, která by z 92 protonů a 143 neutronů spočítala dostatečně přesně
řekněme poločas rozpadu 235U, není známa.
3
jádraelektrony
makroskopickévlastnosti
molekulárnímodel
-kvantová mechanika
(kvantové simulace)@@@@@@R
kvantová
mechanika
statis
tická
mec
hani
ka
kine
tická
teor
ie
simul
ace
Obr. 1: Dělba práce mezi kvantovou chemií a statistickou termodynamikou
v předkvantové éře) je matematické vyjádření apriorních představ o velikosti a tvaru molekula silách mezi nimi.
I když máme model, není vždy jednoduché získat chování látky. Podle přístupů a použi-tých aproximací používáme následující postupy
• kinetická teorie plynů vychází z představy molekul pohybujících se prostorem a občasse srážejících,
• statistická termodynamika zkoumá chování systémů mnoha částic statisticky na zá-kladě pojmu pravděpodobnosti,
• molekulární simulace na základě molekulárního modelu (a výsledků obou výše uvede-ných metod) „uvaříÿ v počítači model látky složený z mnoha molekul.
Mezi kvantovou chemií a těmito metodami tedy dochází k dělbě práce (viz obr. 1).Na základě představy hmoty složené z atomů je možné jednoduše odvodit mnoho vztahů,
které jsme zatím brali pouze jako experimentální poznatky, např. stavovou rovnici ideálníhoplynu či závislost (přesněji nezávislost) vnitřní energie ideálního plynu na objemu. Počítatlze viskozitu plynů i rychlost chemických reakcí.
2 Ideální plyn v kinetické teorii
Kinetická teorie plynů vychází z představy atomů nebo molekul, které se pohybují prosto-rem. Pohyb je neuspořádaný (chaotický), a proto je možné jej popisovat pouze statisticky.Kinetická energie je mírou teploty plynu. Zde se budeme zabývat pouze ideálním plynem.Jeho molekuly jsou tak malé, že se (téměř) nesrážejí.
2.1 Stavová rovnice ideálního plynu
Tlak plynu je způsoben nárazy molekul plynu na stěnu. Při odvození vztahu pro výpočettlaku budeme uvažovat N molekul plynu o hmotnosti m uzavřených v krychli o straně L,kterou umístíme do počátku souřadnicového systému tak, aby její stěny byly kolmé na osy.Rychlost libovolné i-té molekuly je vektor ~vi = (vi,x, vi,y, vi,z), její hmotnost mi.
4
-
6
*
HHH
HHH
HHY
*
x
y
L
Obr. 2: Pohyb molekuly v krychli o straně L.
Sledujme nyní molekulu i, viz obr. 2, a její nárazy do stěny krychle v x = 0. K tomunám stačí složka rychlosti vi,x ve změru osy x. Molekula překoná vzdálenost L mezi kolmýmistěnami za čas L/vi,x. O srážkách molekul se stěnou předpokládáme, že jsou dokonale pružné,tedy že úhel dopadu se rovná úhlu odrazu a mění se pouze směr rychlosti, nikoli její velikost.Po odrazu je tedy vi,x stejně veliké s opačným znaménkem. Složky vi,y a vi,z se nemění.Molekula podruhé do stejné stěny tedy narazí za čas t = 2L/vi,x.
Síla je rovna změně hybnosti za jednotku času. Hybnost molekuly i je ~P = m~v. Protožese molekula odrazí opačným směrem, je změna hybnosti při nárazu ve směru osy x rovna∆Px = 2mivi,x. Průměrná síla způsobená nárazy jedné molekuly je tedy
Fi,x =∆Pxt
=2mivi,x2L/vi,x
=miv
2i,x
L
Tlak je síla ode všech N molekul dělená plochou, tedy
p =
∑Ni=1 FiL2
=
∑Ni=1miv
2i,x
L3
Kinetická energie jedné molekuly je
12miv
2i =
12mi(v
2i,x + v2
i,y + v2i,z)
Protože předpokládáme, že se molekuly pohybují náhodně ve všech směrech, jsou příspěvkyod všech tří souřadnic stejné. Kinetická energie plynu (všech N molekul) je tedy
Ekin =12
N∑i=1
miv2i =
32
N∑i=1
miv2i,x
a tlak můžeme vyjádřit jako
p =
∑Ni=1miv
2i,x
L3=
23Ekin
V
Jinak napsáno
pV =23Ekin (1)
5
Ale kde je teplota? Jak tento vztah souvisí se stavovou rovnicí ideálního plynu, kterouznáte ve tvaru pV = nRT? Potíž je v tom, že mechanika sama o sobě nezná pojem teploty.Teplota je pojem termodynamiky. Srovnáním (1) s pV = nRT dostaneme tvrzení, že teplotaje mírou kinetické energie molekul.
2.2 Ekvipartiční princip
Rozeberme tento vztah podrobněji. Zaveďme nejprve pojem počet mechanických stupňůvolnosti f . To je počet mechanických proměnných (souřadnic), kterými je soustava popsána.V našem případě je f = 3N , protože každá molekula je popsána trojrozměrným vektorem.Pak
pV = nRT =N
NART =
f
3kBT =
23Ekin
kde kB = R/NA = 1,3806503·10−23 J K−1 je Boltzmannova konstanta. Výraz Ekin je složencelkem z f členů tvaru 1
2miv2i,k, kde k je x, y nebo z. V průměru tedy na každý stupeň
volnosti připadá energieEkin
f=
12kBT (2)
To je speciální případ tzv. ekvipartičního principu. Vrátíme-li se zpět k molárním jed-notkám (f = 3NA), máme pro izochorickou molární tepelnou kapacitu jednoatomového plynu
CVm =
(∂Um
∂T
)V
=
(∂Ekin
∂T
)V
=32R (3)
což dobře souhlasí s experimentem pro vzácné plyny. U lineárních molekul jako N2 přistupujíještě rotace, které jsou popsány dvěma stupni volnosti, tedy CVm = 5
2R.
Ve skutečnosti jsou dva atomy N popsány celkem šesti stupni volnosti, každý atom třemi. Vzhledem k sílevazby lze však ukázat, že stupeň volnosti odpovídající vibracím nejen že nelze popisovat klasickou mechanikou, aleza běžných teplot budou prakticky všechny molekuly v základním stavu. Proto ho neuvažujeme. Právě vzhledemk významnosti kvantových efektů je použitelnost ekvipartičního teorému omezena na nejjednodušší molekuly (aomezený obor teplot).
Kdybychom však vibrační stupně volnosti (např. pro těžší atomy a vyšší teploty) mohli popisovat klasicky,dostali bychom kBT na každý vibrační stupeň volnosti. To je proto, že člen 12kBT se počítá za každý kvadratickýčlen v energii. V připadě translací a rotací máme jeden kvadratický člen, např. 12mv
2x, na stupeň volnosti v kinetické
energii, v případě vibrací máme navíc harmonický (kvadratický) potenciál (44).
2.3 Energetická rovnice
Pro úplný termodynamický popis látky nám stavová rovnice nestačí. Jednotlivé molekulyideálního plynu navzájem neinteragují, a proto je jedno, jak jsou daleko od sebe; vnitřníenergie ideálního plynu proto nezávisí na objemu (ani tlaku). Závisí však (pro molekulys vnitřními stupni volnosti) na teplotě, Uid. plyn = U(T ); Pro jednoatomové molekuly bezvnitřních stupňů volnosti dokonce víme, jak: z (1) dostaneme
U = Ekin =32pV =
32nRT
6
Z obou rovnic (stavové a energetické) za použití druhého termodynamického zákona pak(přes Carnotův cyklus) dostaneme, že teplota T zavedená výše je tzv. absolutní termo-dynamickou teplotou, tedy po dělení elementárního vyměněného tepla dQ touto teplotoudostaneme úplný diferenciál stavové funkce (totiž entropie): dS = dQ/T .
3 Statistická termodynamika
Z kvantové teorie víte, že systém (jeden atom, jedna molekula, 6·1023 molekul. . . ) se nacházív mnoha kvantových stavech – vlastních stavech Hamiltoniánu systému. Ve statistické ter-modynamice se jim říká také mikrostavy. Budeme je označovat symbolicky ψ. Pro velkésystémy je těchto stavů sice neuvěřitelně mnoho, je jich však vždy konečně mnoho3, takžemůžeme bez problémů psát výrazy jako
∑ψ, tj. součet přes všechny stavy. Komplementárním
popisem je přístup klasické mechaniky, kde mikrostav systému sestávajícího, dejme tomu,z N stejných atomů je zcela popsán, jestliže známe polohy ~ri a rychlosti ~vi všech atomů. Po-lohy jsou vektory v trojrozměrném prostoru, ~ri = (xi,yi,zi); místo rychlostí budeme používat(z důvodů, které přesahují rozsah této přednášky) hybnosti ~pi = mi~ri, kde mi je hmotnostčástice číslo i. Vzniká zde ovšem problém – místo sumy
∑ψ zde máme integrály. Pro zjedno-
dušení zápisu budeme proto při odvození základních vztahů používat jazyk kvantové teorie,integrálům ale neuniknete, vrátíme se k nim později.
Základní myšlenkou statistické termodynamiky je, že makroskopické vlastnosti hmoty(řekněme hustota whisky) jsou výsledkem pohybu a vzájemné interakce obrovského množstvímolekul (vody a ethanolu). Kromě dostatečně velkého množství molekul je musíme sledovatpo dostatečně dlouhou dobu. Např. o tlaku plynu lze hovořit pouze tehdy, jestliže nárazůmolekul na stěnu nasbíráme dostatečné množství. Tlak je tedy časovou střední hodnotouokamžitého tlaku – jednotlivých nárazů. Tomuto makroskopickému zprůměrovaného projevuvšech mikrostavů říkáme makrostav.
Jednotlivé mikrostavy se ve vyvíjejícím se systému (tzv. trajektorii) vyskytují s určitoupravděpodobností. Budeme se snažit tuto pravděpodobnost určit. Množina všech mikrostavůdaného systému včetně přiřazených pravděpodobností se nazývá (statistický) soubor. Mů-žeme si ho představit i tak, že necháme daný systém vyvíjet po velmi dlouhou dobu a kon-figurace si budeme zapisovat. Časem dostaneme obrovskou množinu (soubor) mikrostavů.
3.1 Pravděpodobnost
3.1.1 Mikrokanonický soubor
Nejjednodušším souborem je mikrokanonický soubor. Je to vlastně izolovaný systém, je-hož objem je konstantní a nevyměňuje si s okolím žádnou energii. Celková energie tohotosystému je proto konstantní; v jazyce kvantové teorie to znamená, že stavy jsou degene-rované. Z přednášek o kvantové teorii víte, že systému v pevném objemu V lze přiřaditvlastní kvantové stavy ψ tak, že každému stavu (přesněji: mikrostavu) odpovídá hodnotavlastní energie E(ψ). Je-li v systému mnoho částic, je počet těchto stavů těžko představitelnéobrovské číslo, což nám však nemůže zabránit ve výpočtech. Kromě toho jsou stavy častodegenerované, tedy jedné hodnotě energie vyhovuje mnoho stavů ψ.
3Za určitých podmínek spočetně mnoho, takže řady jsou nekonečné, to však další popis neovlivní
7
Uvažujme nyní izolovaný systém, nedochází tedy k výměně ani práce (V = konst; jinédruhy práce než objemovou nebudeme uvažovat) ani energie ve formě tepla. Energie systémuE je konstantou a je rovna (až na případnou aditivní konstantu) vnitřní energii systému.Stále však může existovat mnoho stavů s touto energií. Základním předpokladem statistickétermodynamiky je demokracie: jestliže nemáme žádný důvod, proč nějaký stav preferovat,jsou všechny stavy stejně pravděpodobné4. Všechny stavy mají proto stejnou pravděpodob-nost výskytu
π(ψ) =1W
pro všechny stavy ψ takové, že E(ψ) = E
kde π(ψ) označuje pravděpodobnost stavu ψ a W je počet stavů. Makrostav, tedy to, copozorujeme v přírodě, si můžeme představit tak, že tyto mikrostavy rychle přecházejí jedenna druhý. V dlouhém čase je pravděpodobnost výskytu všech stavů stejná.
Z kvantové teorie víte, co je to hodnota pozorovatelné (veličiny) X ve stavu ψ, X(ψ)(nejjednodušším příkladem je energie E , jiným třeba výška molekuly nad hladinou moře).Makroskopická hodnota libovolné veličiny X (tedy to, co nakonec měříme v laboratoři jakovýsledek složený ze všech možných stavů), je dána aritmetickým průměrem přes všechnymožné stavy, tedy
〈X〉 =
∑ψX(ψ)
W(4)
kde symbolem 〈X〉 je označena střední hodnota a X(ψ) je hodnota veličiny X v kvantovémstavu ψ; W =
∑ψ 1 je počet stavů.
Příklad. „Kvantovýÿ systém hrací kostka se po hodu vyskytuje v jednom s W = 6 stavů, kterémůžeme symbolicky zapsat jako 1,2,3,4,5,6. Předpokládejme, že hrajeme velmi primitivní hru:hráč vsadí (dá bankéři) 1 korunu a hodí; padne-li 6, bere 5 Kč, jinak nic. Pro výhru tedy platíX(1) = X(2) = X(3) = X(4) = X(5) = −1 (nic nevyhraje, odevzdal ale 1 Kč) a X(6) = 4(vyhraje 5, ale 1 vsadil). Průměrná výhra je
〈X〉 =
∑ψX(ψ)
W=−1− 1− 1− 1− 1 + 4
6= −1
6
Hráč tedy v průměru prohraje 1/6 Kč v jednom kole. Samozřejmě může náhodou vyhrát. Má-livšak majitel ve svém kasinu 6 hracích stolů a na každém se sehraje 1000 takových her za večer, lzecelkem jistě počítat s inkasem okolo 1000 Kč.
Kvantová mechanika již obsahuje pojem pravděpodobnosti, takže výše uvedený předpoklad je přirozený.Pokud však používáme klasickou mechaniku, máme místo vlastní funkce ψ tzv. konfiguraci molekul, tedy prokaždou molekulu musíme znát polohu všech jejích atomů a rovněž všechny rychlosti. Jestliže toto známe, umímev principu předpovědět vývoj systému, tedy polohy i rychlosti všech atomů v libovolném budoucím čase. Výšeuvedený předpoklad pak znamená, že každou konfigurací (o dané energii) systém v nekonečně dlouhém časeprojde stejněkrát. Tento předpoklad se nazývá ergodická hypotéza. I když byla pro některé jednoduché systémyvyvrácena, pro složité systémy, které se vyvíjejí chaoticky, vyhovuje dobře.
4Pokud to neodporuje nějakému zákonu zachování, třeba hybnosti
8
3.1.2 Kanonický soubor
Předpokládejme nyní, že dovolíme systému vyměňovat teplo s okolím, tedy ponoříme ho dotermostatu (který si můžeme představit jako jiný mnohem větší systém). Energie systémujiž nebude přesně konstantní, ale bude fluktuovat okolo jisté střední hodnoty. (Pokud jevám tato představa cizí, představte si, že náš systém je tvořen jen jednou molekulou, třebasvelkou. Ta bude občas narážet do molekul termostatu a jednou bude mít energii vyšší, jednounižší.) Poznamenejme ještě, že příčina fluktuace energie našeho systému je v tom, že systémpřechází z jedné energetické hladiny na jinou, protože za konstantního objemu jsou vlastnífunkce ψ i energetické hladiny neměnné. Nemůžeme však již předpokládat, že všechny stavy(s různou energií) jsou stejně pravděpodobné, ale naopak jejich zastoupení je dáno jistoufunkcí π(E), která je však již jen funkcí energie, neboť stavy se stejnou energií mají stejnoupravděpodobnost.
Pokusme se nyní odvodit tvar funkce π(E). K tomu předpokládejme, že ve styku s ter-mostatem máme dva (nikoliv nutně stejné) systémy, 1 a 2. Pravděpodobnost, že systém 1 máenergii E1, je π(E1), pravděpodobnost, že systém 2 má energii E2, je π(E2). Nyní systémyspojíme do jednoho, ale tak, aby se minimálně ovlivňovaly. Energie spojeného systému budesoučtem energií obou subsystémů, E1+2 = E1 +E2. Pravděpodobnost, že systém 1 naleznemeve stavu s energií E1 a zároveň systém 2 nalezneme ve stavu s energií E2, je dána součinemobou pravděpodobností,
π(E1+2) = π(E1 + E2) = π(E1)π(E2)
Po zlogaritmování máme
lnπ(E1 + E2) = lnπ(E1) + lnπ(E2) (5)
Pokud má tato rovnice platit pro jakékoliv hodnoty E1 a E2, musí být lnπ(E) lineární funkcíE, tedy
lnπ(Ei) = αi − βEi =⇒ π(Ei) = eαi−βEi (6)
kde αi a β jsou konstanty. Konstanta αi je přitom aditivní, tedy α1+2 = α1 + α2, a zá-visí jak na velikosti systému tak na (zatím neznámé) teplotě. Je dána normalizací, součetpravděpodobností se totiž musí rovnat jedničce,
1 =∑ψ
π(ψ) =∑ψ
eα−βE(ψ) = eα∑ψ
e−βE(ψ) =⇒ e−α =∑ψ
e−βE(ψ)
kde index i u α jsme již pro jednouchost vynechali (vztah platí pro obecný systém). Fyzikálnívýznam α odvodíme až v odd. 3.2.
Konstanta β je daná vlastnostmi termostatu. Veličina, která je stejná v termostatu iv systému, je podle nultého zákona (nulté věty) empirická teplota, tedy β bude souvisets teplotu. Veličina β nezávisí na systému – teplota je u všech systémů v rovnováze stejná.Faktor eα−βE(ψ) resp. e−βE(ψ) se nazývá Boltzmannova pravděpodobnost resp. váha čifaktor.
Vztah (4) pro střední hodnotu musíme zobecnit. V (4) byla pravděpodobnost stavuπ(ψ) = 1/W , nyní je dána funkcí energie (6), π(ψ) = π(E(ψ)). Každý stav započtemes pravděpodobností, se kterou se vyskytuje,
〈X〉 =∑ψ
X(ψ)π(E(ψ)) =∑ψ
X(ψ)eα−βE(ψ) =
∑ψX(ψ)e−βE(ψ)∑
ψ e−βE(ψ)(7)
9
Např. vnitřní energie je nyní dána střední hodnotou
U = 〈E〉 (8)
Možná je na místě poněkud ukecaná vsuvka ilustrující princip střední hodnoty. Řekněme, že studujeme vlivnadmořské výšky na obsah hemoglobinu v krvi. Zajistíme si vzorek obyvatel, které žijí ve výšce h = 0 m nad mořem(dejme tomu přesněji v rozsahu výšek 0 až 1 m) a zjistíme, že mají v krvi určité množství hemoglobinu, kteréoznačíme X(0), a stejně pro další výšky. Známe tedy funkci X(h), která udává průměrné množství hemoglobinuv krvi (tato funkce s výškou roste). Jaké je ale průměrné množství hemoglobinu v krvi náhodně vybraného člověka?Jistě nemůžeme průměrovat jen funkci X(h), protože na pobřeží žije víc lidí než ve velehorách. Označme početlidí žijící ve výšce h symbolem N(h). Tímto číslem budeme vážit zjištěné množství hemoglobinu. Střední obsahhemoglobinu v krvi je
〈X〉 =∑hX(h)N(h)∑
hN(h)
Toto je analogií pravé strany rov. (7). Všimněte si také, že∑hN(h) je celkový počet lidí. Ve spojitém případě
(výšky nedělíme po metru) přejdou sumy na integrály:
〈X〉 =
∫X(h)N(h)dh∫N(h)dh
Dále si můžeme definovat pravděpodobnost, že člověk žije ve výšce h, vztahem π(h) = N(h)/∑hN(h). Pak
ovšem〈X〉 =
∑h
X(h)π(h)
případně
〈X〉 =∫X(h)π(h)dh
Zbývá určit konstantu β. Protože vztahy platí pro libovolné systémy, zvolíme si ta-kový, který umíme spočítat. Nejjednodušším systémem je asi jednoatomový ideální plyn(složený z hmotných bodů). Vypočtěme střední hodnotu energie jedné jeho molekuly. Je-li~v = (vx,vy,vz) rychlost molekuly (vektor) a m její hmotnost, je její energie E = 1
2m~v2 =
12m(v2
x + v2y + v2
z). Podle (7) máme (pro spojité hodnoty ~v přejde součet na integrál)
〈E〉 =
∑ψ E(ψ)π(E(ψ))∑
ψ π(E(ψ))=
∫12m~v
2π(12m~v
2) d~v∫π(1
2m~v2) d~v
(9)
kde∫
d~v =∫∞−∞
∫∞−∞
∫∞−∞ dvxdvydvz (integruje se přes všechny možné hodnoty rychlosti). Po
poněkud nepříjemném výpočtu dostaneme
〈E〉 =32
1β
(10)
Tento výpočet bohužel vyžaduje trochu matematiky. Snad si ještě vzpomenete na vzorec (Gaussův integrál)∫ ∞−∞
exp(−Ax2) =√π/A (11)
Užitečným trikem nejen ve statistické termodynamice je derivace podle parametru. Vztah (11) zderivujeme podleparametru A,
∂
∂A
∫ ∞−∞
exp(−Ax2) =∂
∂A
√π/A∫ ∞
−∞(−x2) exp(−Ax2) = −1
2√π/A3/2
10
Tento vztah spolu s (11) použijeme v (9) a dostaneme výsledek (10).
Energie jednoho molu je NA krát větší než energie jedné molekuly, U = NA32
1β, což se
má rovnat U = 32RT podle rovnice (3). Z toho dostaneme hledaný vztah pro β
β =1kBT
kdekB = R/NA = 1,38065·10−23 J K−1
je tzv. Boltzmannova konstanta. Boltzmannův faktor je tedy
e−E(ψ)/(kBT ) (12)
Ve výše uvedených vztazích je E skutečná energie daného souboru molekul (vyjádřená v J). Budeme-li, jakje v chemii zvykem, vyjadřovat energii na jeden mol, tedy v J mol−1, bude β = 1/(RT ) a příslušné výrazy budoumít povědomý tvar exp[−energie/(RT )].
Zkusme se podívat, co z toho, co znáte z fyzikální chemie, se dá vysvětlit na základěBoltzmannovy pravděpodobnosti.
• Aby mohla proběhnout chemická reakce, musí reagující látky mít určitou minimálníenergii, tzv. aktivační energii E∗. Podíl molekul, které získají tuto energii a zreagují zajednotku času, je úměrný Boltzmannovu faktoru exp
(−E∗RT
). Pro rychlostní konstantu
pak dostaneme Arrheniův vztah, k = A exp(−E∗RT
), kde A je (prakticky) konstanta.
• Energie potřebná k přenesení molekuly (v „chemickýchÿ jednotkách jednoho molu)z kapaliny do páry je ∆výpHm, pravděpodobnost nalezení molekuly v páře je pak
úměrná exp(−∆výpHRT
). Pravděpodobnost nalezení molekuly je úměrná hustotě a ta
je v ideálním plynu úměrná tlaku. Výsledek souhlasí s integrovaným tvarem Clausi-ovy-Clapeyronovy rovnice
p = p0 exp
[−∆výpH
R
(1T− 1T0
)]= const · exp
(−∆výpH
RT
)který byl odvozen za předpokladu, že ∆výpHm nezávisí na teplotě.
Další aplikace je obzvlášť instruktivní.
Příklad. Vypočtěte tlak vzduchu ve výšce h = 8850 m za teploty 0 C, je-li na hladině mořenormální tlak. Předpokládejte, že zemská atmosféra je v rovnováze, tedy zanedbejte vítr, vlivslunečního záření atd.Řešení. Energie molekuly o rychlosti v je dána vztahem
E = mgh+12m~v2
kde g je tíhové zrychlení. Druhý člen je (v průměru) pro všechny molekuly stejný, protože teplotav rovnováze je nezávislá na výšce, a nemusíme jej proto uvažovat (tj. členy jej obsahující se vykrátí).Hustota (počet molekul v objemu) je úměrná pravděpodobnosti π(E) = exp(−βmgh). Tlak je
11
úměrný pro ideální plyn hustotě, tedy této pravděpodobnosti, a hledaná závislost tlaku na výšce(tzv. barometrická rovnice) je
p = pst exp(−βmgh) = pst exp
(−mghkBT
)= pst exp
(−Mgh
RT
)(13)
Dosadíme střední molární hmotnost vzduchu M = 29 g mol−1 a máme
p = 101,325 exp
(−29·10−3 · 9,81 · 8850
8,314 · 273
)= 33,4 kPa
Barometrickou rovnici můžeme také získat ze stavové rovnice ideálního plynu a podmínky, žetlak v určité výšce je dán tíhou veškerého vzduchu nad daným místem. Uvažujme pokles tlakudp = p(h + dh) − p(h) při změně výšky o dh. Hydrostatický tlak sloupce vzduchu o výšce dh je−dh g% (tlak s výškou klesá, proto záporné znaménko). Hustota je podle stavové rovnice ideálníhoplynu rovna % = Mp/(RT ), a proto
dp = −dh g% = −dhgMp
RT
po separaci proměnných a integraci s podmínkou p(0) = pst∫ p
pst
dpp
= −∫ h
0dh
gM
RT
dostaneme opět barometrickou rovnici (13).Nad vzorcem je užitečné se zamyslet. Nechť je molekula v určité výšce. Pod vlivem nárazů
ostatních molekul a tíhové síly se může pohybovat nahoru i dolů; pravděpodobnost, že se budepohybovat nahoru je vždy menší než že se bude pohybovat dolů. Máme-li molekulu u hladinymoře, dolů se již ovšem utéci nemůže. Po dosazení h = 80 m snadno zjistíme, že ve výšce 80 m jepravděpodobost výskytu molekuly o 1 % menší než u hladiny (tedy je 0.99), stejně tak ve výšce 160m je pravděpodobost výskytu molekuly o 1 % menší než ve výšce 80 m, celkem tedy 0.992 hodnotyu hladiny, protože na každém místě se může nezávisle „rozhodnoutÿ, zda poputuje nahoru či dolů,přičemž preference směru dolů je v každém místě stejná. Jinými slovy, pravděpodobost výskytumolekuly ve výšce 160 m je součinem pravděpodobností, že popoleze dvakrát o 80 m, zatímco jejíenergie je dvojnásobná.
Představme si nyní, že máme molekuly dvě; jsou obě v ideálním plynu, takže se neovlivňují.Pravděpodobnost, že nalezneme obě dvě ve výšce 80 m, je rovna druhé mocnině pravděpodobnosti,že tam nalezneme jednu, totiž 0.992 násobku pravděpodobnosti molekul u hladiny moře. Faktor0.992 je stejný jako pro jednu molekulu u hladiny moře a druhou ve výšce 160 m – protože v oboupřípadech má spojený systém obou molekul stejnou energii, a pravděpodobnost závisí jen na energii.
Při aplikaci vztahu (7) na velké systémy je užitečné si uvědomit ještě jednu skutečnost.I když ve vzorci máme funkci e−E(ψ)/(kBT ) klesající s teplotou5, zastoupení stavu o dané energii(pravděpodobnost nalezení energie E u systému ve styku s termostatem) není klesající funkcíE, ale má ostré maximum (pík). To je proto, že počet stavů rychle roste s energií, zpravidlajako vyšší mocnina energie (srov. např. (26)). Zastoupení stavu s danou energií u velkéhosystému je dáno jako součin rychle rostoucího počtu stavů a rychle klesající Boltzmannovypravděpodobnosti, viz obr. 3. Šířka píku se zmenšuje s velikostí systému (úměrně N−1/2, kdeN je počet molekul).
5Uvědomte si, že E(ψ) měříme vždy od tzv. základního stavu, tedy stavu s nejnižší energií, a proto je E(ψ)vždy kladné. Kdyby stav s nejnižší energií neexistoval, energie by se stále snižovala a systém by explodoval.
12
1 2 3 4 5E
0
5
10
15
20
π(E)
Obr. 3: Pravděpodobnost stavu s danou energií ( ). Křivka : Boltzmannův faktor (12) pro1/(kBT ) = 10, křivka : počet stavů úměrný E30
3.2 Termodynamika v kanonickém souboru
Vztahy předchozího oddílu nám umožňují pracovat s pravděpodobnostmi stavů a počítatstřední hodnoty veličin, které lze vyjádřit jako funkce konfigurace. Takovými funkcemi jsounapříklad tlak a energie (a tedy i entalpie). Neumožňují však spočítat entropii a odvozenéveličiny (Helmholtzova a Gibbsova energie), které tvoří jádro chemické termodynamiky.
Abychom odvodili korespondenci s termodynamikou, napišme vnitřní energii, která jestřední hodnotou energie
U =∑ψ
π(ψ)E(ψ)
Malá změna této veličiny je
dU =∑ψ
π(ψ) · dE(ψ) +∑ψ
dπ(ψ) · E(ψ) (14)
Člen dE(ψ) znamená, že se změnila energetická hladina (energie stavu ψ), člen dπ(ψ) zna-mená, že se změnila pravděpodobnost výskytu stavu ψ. Z termodynamiky víme, že se tatozměna má rovnat
dU = −p dV + TdS (15)
Rovnají se oba vztahy?První člen v (14) odpovídá vlivu změny energetických hladin E(ψ). Aby se tyto hladiny
změnily, musíme změnit vnější podmínky, například objem. Představme si proto, že mámeve válci s pístem systém ve stavu ψ a tento píst posuneme o malé dx. Tím se změní vlastnífunkce ψ a následně se energie změní o dE(ψ). Tato změna energie se rovná mechanické práci,tedy až na znaménko součinu síly F a dráhy dx, tedy dE(ψ) = −Fdx = −F/A · d(Ax) =−p(ψ) dV , je to tedy práce objemová (A je plocha pístu). Veličina p(ψ) je tlak danéhostavu, skutečný tlak dostaneme jako střední hodnotu (vážený průměr) přes všechny stavy,tedy p =
∑ψ π(ψ)p(ψ). První člen v (14) je tedy objemová práce −p dV .
Druhý člen v (14) naopak zahrnuje vliv změny pravděpodobnostního rozložení π(ψ) zakonstantního objemu za neměnných stavů ψ i hladin energie E(ψ). Bude tedy odpovídat
13
výměně tepla; např. při dodání tepla vzroste zastoupení stavů s vyšší energií na úkor stavůs energií nižší. Tento člen by tedy měl odpovídat členu TdS ve vztahu (15). Abychom hle-danou korespondenci dostali, vyjádřeme nejprve E(ψ) pomocí pravděpodobnosti z (6),
E(ψ) =1β
[α− lnπ(ψ)]
a dosaďme do druhého členu v (14),∑ψ
dπ(ψ)E(ψ) =∑ψ
dπ(ψ)1β
[α− lnπ(ψ)] = − 1β
∑ψ
dπ(ψ) · lnπ(ψ)
Při úpravě jsme použili vztah ∑ψ
dπ(ψ) = 0 (16)
který snadno odvodíme diferencováním normovací podmínky∑
ψ π(ψ) = 1 (vyjádřeno slovy,abyste se nebáli derivací: když π(ψ) vzroste, dπ(ψ) je kladné, když π(ψ) klesne, je dπ(ψ)záporné; protože součet je jedna, musejí některé dπ(ψ) vzrůst a jiné klesnout). Rovnici (16)použijeme nyní ještě jednou, abychom ukázali, že∑
ψ
dπ(ψ) · lnπ(ψ) = d∑ψ
π(ψ) lnπ(ψ)
Druhý sčítanec v (14) je tedy
∑ψ
dπ(ψ) · E(ψ) = −kBT d
[∑ψ
π(ψ) lnπ(ψ)
]
Porovnáním s TdS dostaneme entropii jako funkci pravděpodobnosti stavů,
S = −kB
∑ψ
π(ψ) lnπ(ψ) (17)
Vrátíme-li se zpátky k izolovanému systému, je π(ψ) = 1/W pro E = E(ψ) a π(ψ) = 0 proE 6= E(ψ) a máme
S = kB lnW (18)
kde W je počet stavů. Toto je slavná Boltzmannova rovnice pro entropii. Entropie je tedytím větší, čím více způsoby lze daný stav realizovat, přičemž závislost je logaritmická. Tosouvisí s aditivitou entropie: Uvažujeme-li opět spojené systémy 1+2, platí W1+2 = W1W2,a proto S1+2 = S1 + S2.
3.2.1 Helmholtzova energie
Zbývá určit fyzikální význam parametru α. Z (17) po dosazení za lnπ z (6) dostaneme
S = −kB
∑ψ
π(ψ)[α− βE(ψ)] = −kBα +U
T(19)
14
a tedy
α =U − TSkBT
=F
kBT
kde F je vám dobře známá Helmholtzova energie. Z toho zároveň plyne důležitý vztah
F = −kBT ln
[∑ψ
e−βE(ψ)
](20)
kde suma v hranatých závorkách je stejná jako ve jmenovateli (7), nazývá se kanonickápartiční funkce nebo statistická suma a budeme ji značit Z. Z (20) lze alespoň v prin-cipu spočítat libovolnou termodynamickou veličinu (je to opakování Fyzikální chemie probakaláře)
p = −∂F∂V
(21)
S = −∂F∂T
(22)
U = F − TS =∂βF
∂β(23)
H = U + pV (24)
G = F + pV (25)
3.2.2 Ideální plyn
Vypočteme tlak a molární entropii jednoatomového ideálního plynu za teploty T v objemu V .Podle výsledků, které jste získali na přednáškách z kvantové teorie, je energie jedné molekulyv nádobě tvaru kvádru o stranách a, b a c s pevnými stěnami rovna
E =h2
8m
(n2x
a2+n2y
b2+n2z
c2
)(26)
Kanonická partiční funkce je
Z1 =∞∑
nx=1
∞∑ny=1
∞∑nz=1
exp(−βE)
Je-li nádoba dost velká, můžeme nahradit sumaci integrací
Z1 =∫ ∞
0
∫ ∞0
∫ ∞0
exp(−βE) dnxdnydnz
a použít opět vzorce (11). Dostaneme
Z1 =V
Λ3
kde V = abc je objem a
Λ =h
(2πmkBT )1/2(27)
15
je tzv. de Broglieova tepelná vlnová délka. Molekuly ideálního plynu navzájem neinteragují,a proto E =
∑Ni=1 Ei, kde se sčítá přes všechny molekuly v systému. Kdyby byly všechny
molekuly různé (ale se stejnou hmotností m), dostali bychom
Z =N∏i=1
Zi = ZN1
Tak tomu ale není. V kvantovém popisu, ze kterého jsme vyšli, jsou stejné atomy nerozliši-telné6. Partiční funkci je proto nutno vydělit počtem způsobů, kterými lze molekuly označit(očíslovat). Nejsme totiž v principu schopni rozlišit stav, kdy je molekula vody č. 1 v le-vém horním rohu nádoby a molekula vody č. 2 v pravém spodním rohu nádoby od stavu,kdy jsou tyto molekuly opačně očíslovány. Počet očíslování N stejných molekul se rovnáN ! = 1 · 2 · 3 · · ·N . Správný výsledek je proto
Z =1N !ZN
1
a podle (20) je Helmholtzova energie rovna
F = −kBT lnZ = −kBT
(N ln
V
Λ3− lnN !
)= −kBT
(N ln
V
Λ3−N lnN +N
)Pro odhad lnN ! jsme použili aproximaci spolu s integrací per partes (′ značí derivaci)
lnN ! = ln 1 + ln 2 + · · ·+ lnN =N∑i=1
ln i
≈∫ N
1lnx dx =
∫ N
1x′ lnx dx = [x lnx]N1 −
∫ N
1x(lnx)′ dx = N lnN −N
Tlak je dán vztahem
p = −∂F∂V
=NkBT
V=nRT
V
což je stavová rovnice ideálního plynu. Výpočet entropie je poněkud zdlouhavý, protože nateplotě závisí i Λ. Vyjde tzv. Sackurova–Tetrodeova formule
S = −∂F∂T
= kBN
(ln
V
NΛ3+
52
)a po přepočtu na mol
Sm = R
(ln
V
NAΛ3+
52
)= konst +R lnV + CV lnT
kde CV = 32R, což je správná hodnota pro jednoatomový plyn.
6Z hlediska klasické teorie není tento předpoklad vůbec samozřejmý. Jeho vynechání však vede k chybnýmhodnotám entropie, např. entropie vzroste po „smícháníÿ dvou stejných vzorků stejného plynu (tzv. Gibbsůvparadox). Vy víte, že entropie vzroste pouze po smíchání různých plynů.
16
3.2.3 Klasická Helmholtzova energie
Velké systémy zpravidla nepopisujeme kvantově (byť s aproximacemi) jako ve výše uvede-ném příkladu, ale klasicky. Předpokládejme, že systém se skládá z N atomů, jejichž polohypopíšeme N vektory ~ri, i = 1, 2 . . . N . Místo rychlostí ~vi, použijeme (z důvodů, které vyžadujíhlubší znalosti mechaniky) hybnosti, ~pi = m~vi. Energie systému atomů je E = Epot + Ekin,kde o potenciální energii Epot (často se označuje i U) se více dovíte v následující kapitolea kinetická energie je Ekin =
∑Ni=1 p
2/(2m). V tzv. semiklasické limitě (provedené ve výšeuvedeném příkladu pro ideální plyn) přejde (20) na
F = −kBT ln
[1
N !h3N
∫e−βEd~r1 . . . d~rN d~p1 . . . d~pN
](28)
= −kBT ln
[1
N !Λ3N
∫e−βEpotd~r1 . . . d~rN
](29)
Tento vztah budeme často používat. Až na faktor 1/(N !h3N) je to přímočará aplikace (20)na klasický soubor molekul. Člen 1/N ! vyplývá z nerozlišitelnosti molekul a jeho zanedbánívede ke špatné hodnotě entropie. Člen 1/h3N vyplývá z kvantové teorie; bez něj se F liší okonstantu – je vzhledem ke špatnému standardnímu stavu. Výraz
QN =∫
e−βEpotd~r1 . . . d~rN
se nazývá konfigurační integrál; je to integrál z Boltzmannova faktoru přes všechny konfi-gurace systému. Ještě uveďme výraz pro střední hodnotu veličiny závislé pouze na polohách(konfiguraci). Podle (7) platí
〈X〉 =
∫X(~r1, . . . ,~rN)e−βEpotd~r1 . . . d~rN∫
e−βEpotd~r1 . . . d~rN(30)
Příklad. Aplikujte rovnici (29) na jednoatomový ideální plyn a vypočtěte tlak a chemický poten-ciál.Řešení. Pro ideální plyn platí Epot = 0. Integrál v (29) je proto∫
d~r1 . . . d~rN =∫V
d~r1 ×∫V
d~r2 × · · · ×∫V
d~rN = V N
Dále použijeme Stirlingovu aproximaci lnN ! ≈ N lnN −N . Vyjde
F = NkBT
[lnN − 1 + ln
(Λ3
V
)]Z toho snadno vypočteme tlak
p = −(∂F
∂V
)T
= −NkBT∂ ln(1/V )
∂V=NkBT
V=nRT
V(31)
což je stavová rovnice ideálního plynu, ježto N = NAn a R = NAkB. Chemický potenciál jednoa-tomového ideálního plynu spočteme nejsnáze z rovnice pro otevřený systém
dF = −SdT − pdV + µdN
17
Snadno zderivujeme
µ =
(∂F
∂N
)T ,V
= kBT
[lnN − 1 + ln
(Λ3
V
)]+NkBT
[1N
]= kBT ln
(NΛ3
V
)= kBT ln
(pΛ3
kBT
)kde v poslední úpravě jsme použili stavovou rovnici ideálního plynu, pV = NkBT .
3.3 Pravděpodobnostní rozložení molekulárních rychlostí
Rychlost molekul plynu se vlivem vzájemných srážek neustále mění, takže stanovení skutečnérychlosti vybrané molekuly v daném okamžiku není možné. Ve velkém souboru molekullze však pomocí statistických zákonitostí získat informaci o zastoupení rychlostí, jimiž semolekuly pohybují. Pravděpodobnost, že molekulu 1 nalezneme
• v krychličce o velikosti dx1dy1dz1 se souřadnicemi v intervalu (x1, x1 + dx1), (y1, y1 +dy1) a (z1, z1 + dz1) a zároveň
• s rychlostmi v intervalu (v1x, v1x + dvx), (v1y, v1y + dvy), (v1z, v1z + dvz),
• a stejně pro molekuly 2, 3, . . . , N ,
je úměrná Boltmannovu faktoru (12)
exp
(−Epot + Ekin
kBT
)= exp
(−Epot
kBT
)exp
(−12m1v
21x
kBT
)exp
(−1
2m1v21y
kBT
)exp
(−12m1v
21z
kBT
)· · ·
Pravděpodobnost je vyjádřena jako součin, protože jednotlivé složky kinetické energie 12m1v
21x,
12m1v
21y atd. jsou nezávislé navzájem i na potenciální energii. Proto pravděpodobnost, že x-
-ová složka rychlosti částice 1 je v intervalu (v1x, v1x + dvx), je úměrná
exp
(−12m1v
21x
kBT
)bez ohledu na polohy a rychlosti ostatních částic a složky rychlosti částice 1 v osách y a z.Stejný vztah platí pro všechny kartézské složky rychlosti a všechny molekuly. Rovnici jezvykem normalizovat, tedy znásobit konstantou tak, aby normalizované rozložení7 fx(vx)dvxudávalo přímo pravděpodobnost, že rychlost ve směru x je v intervalu (vx, vx + dvx). Pakplatí normalizační podmínka ∫ +∞
−∞fx(vx)dvx = 1
tj. pravděpodobnost, že rychlost je jakákoliv (kdekoliv v intervalu (−∞,+∞)), je jedna.Takto normalizovaná funkce se nazývá hustota pravděpodobnosti. Podle Gaussova inte-grálu (11) snadno najdete, že
fx(vx) =√
m1
2πkBTexp
(−12m1v
21x
kBT
)7Také se používají termíny (pravděpodobnostní) rozdělení či distribuce.
18
-1000 -500 0 500 1000
vx/ms-1
0
0.0005
0.001
0.0015
fx(vx)/sm-1
600 K
300 K
0 500 1000 1500 2000
v/ms-1
0
0.0002
0.0004
0.0006
0.0008
0.001
0.0012
f(v)/sm-1600 K
300 K
Obr. 4: Maxwellovo–Boltzmannovo rozložení rychlostí pro argon za různých teplot. Vlevo: hus-tota pravděpodobnosti nalezení složky vx rychlosti, vpravo: hustota pravděpodobnosti nalezenírychlosti v
Ti, kdo dobře poslouchali přednášky z matematické statistiky, si jistě všimli, že fx(vx) nenínic jiného než Gaussovo rozložení se směrodatnou odchylkou σ = kBT/m.
Funkce fx(vx) je pro dvě teploty vynesena na obr. 4 vlevo. Se zvyšující teplotou sedistribuce rozšiřuje, neboť stoupá zastoupení vyšších rychlostí.
Rozdělení na tři kartézské složky se někdy hodí, častěji nás ale zajímá, jaká je pravdě-podobnost, že rychlost v = |~v| je v intervalu (v, v+ dv). Protože objem slupky (v, v+ dv) je4πv2dv, je tato pravděpodobnost
f(v) = 4π
(m1
2πkBT
)3/2
exp
(−12m1v
2
kBT
)v2 (32)
Závislost f(v) je na obr. 4 vpravo. Vrchol funkce f(v), tj. nejpravděpodobnější rychlost,se s růstem teploty posunuje k vyšším hodnotám. Plochy pod jednotlivými křivkami jsoustejné, totiž jedna, neboť distrubuce jsou normalizované.
Výše uvedené vztahy se nazývají Maxwellovo nebo také Maxwellovo–Boltzmannovorozložení. Při odvození jsme nikde nepoužili předpoklad, že látka je v plynném stavu.Vztahy by tedy měly platit i pro kapalinu nebo pevnou látku. Omezením je však použitíklasické mechaniky, které je tím méně oprávněné, čím lehčí je atom nebo molekula (kapalnéhélium), čím silnější je interakce (krystal) a čím nižší je teplota.
4 Síly mezi molekulami
Klasický spojitý systém je soubor N částic, jež spolu interagují jistým mezičásticovým poten-ciálem neboli silovým polem. Na úrovni klasických spojitých systémů lze úspěšně studovatmnoho zajímavých aplikací na různých prostorových škálách od mikroskopické počínaje (ka-paliny, biologické makromolekuly) přes mesoskopickou (např. granulární materiály) až pomakroskopickou (třeba kulová hvězdokupa). Máme-li model, můžeme studovat látku meto-dami „molekulárního modelováníÿ, kterým rozumíme nejen simulace, ale obecně jakékolivstudie tvaru či konformace molekul. Často vystačíme „jenÿ s minimem potenciální energie
19
(tzv. optimalizace struktury). Typickým úkolem je interpretace rozptylových a jiných (NMR)experimentů, tzv. refinement, kdy se struktura stanovená s velkými experimentálními chy-bami doladí tak, aby odpovídaly délky vazeb, atomy se nepřekrývaly atd.
4.1 Jednoatomové molekuly
Uvažujme pro jednoduchost nejprve síly mezi sféricky symetrickými objekty – atomy neboionty. Vyjadřujeme je pomocí potenciálu neboli interakční (potenciální) energie, která jepouze funkcí vzdálenosti, r, obou částic, u(r).
4.1.1 Dva neutrální atomy
Dva atomy vzácného plynu (např. argonu) velmi daleko od sebe na sebe prakticky nepůsobí.Jejich interakční energie je nulová, limr→∞ u(r) = 0.
Jsou-li atomy blízko u sebe, musí se odpuzovat. Z kvantové mechaniky lze získat proenergii odpuzování (repulze) aproximaci
urep(r) ∼ e−αr , α > 0 (33)
Typické hodnoty parametru α jsou řádu 1011 m−1 = 1/(10 pm).Jsou-li molekuly dál od sebe, přitahují se. Pro neutrální molekuly platí aproximace
udisp = −Cr6
(34)
Kvalitativní vysvětlení této disperzní neboli Londonovy energie lze získat pomocí mecha-nismu fluktuující dipól–indukovaný dipól.
Elektrony okolo jádra jsou v neustálém pohybu. I když v průměru je rozložení elektronů okolo jádra řekněmeargonu sféricky symetrické, v daném okamžiku to tak není. V prvním přiblížení lze okamžité rozložení nábojenahradit elektrickým dipólem. Za malý okamžik má dipól jinou velikost a orientaci. Dipól okolo sebe budí elektricképole, jehož intenzita ubývá se vzdáleností jako 1/r3. Jiný atom umístěný do tohoto elektrického poli se polarizuje– vznikne indukovaný dipól o velikosti, která je úměrná intenzitě pole, tj. 1/r3. Tento dipól interaguje s původním,energie je úměrná součinu intenzity pole a velikosti dipólu a podle le Chatelierova–Braunova principu je záporná,protože indukovaný dipól působí proti změně, která jej vyvolala. Energie tedy ubývá se vzdáleností jako 1/r3 ×1/r3 = 1/r6. Poznamenejme ještě, že pro velmi velké vzdálenosti (nad stovky nm) se z důvodu konečné rychlostišíření elektrického pole nestihnou oba dipóly synchronizovat a energie klesá ještě rychleji – jako 1/r7.
Spojením přitažlivé a odpudivé části dostaneme potenciál, kterým se aproximují reálnéinterakce. Můžeme ho zapsat např. takto
u(r) = Ae−Br − C/r6 (35)
Nazývá se exp-6, Buckinghamův nebo též Busingův. Exponenciála ubývá pro dostatečněvelké vzdálenosti rychleji než C/r6, a proto potenciál ve velkých vzdálenostech správně vy-stihuje skutečnou interakci, −C/r6. Při zmenšování vzdálenosti převládne odpuzování danéčlenem Ae−Br. Bohužel při ještě menších vzdálenostech opět převládne záporný člen −C/r6
(jsme v oblasti vzdáleností, kdy předpoklady odvození disperzní interakce neplatí), což jetechnická závada, kterou je nutno za některých okolností opravit, jinak simulace může hava-rovat.
V praxi se často používá jednodušší tvar odpudivé části potenciálu,
urep = const r−m (36)
kde m > 6. Spojením rov. (36) s prvním členem rov. (34) vznikne tzv. m–n neboli Mieůvpotenciál,
uMie(r) = 4ε[(σr
)m−(σr
)n]. (37)
Volba m = 12 a n = 6 dává jednoduchý a přesto poměrně realistický Lennard-Jonesův(LJ) potenciál,
uLJ = 4ε
[(σr
)12−(σr
)6]
= ε
[(σvdW
r
)12− 2
(σvdW
r
)6], (38)
jehož průběh je ukázán na obr. 5. V (38) značí σvdW/2 van der Waalsův poloměr, σ =2−1/6σvdW se někdy nazývá kolizní průměr a ε je hloubka potenciálové jámy.
Protože se ukazuje, že strukturní vlastnosti molekulárních systémů neutrálních částic jsou určeny předevšímkrátkodosahovými silami, studují se systémy tuhých částic kopírujících tvar a velikost molekul a představující taknejhrubší aproximaci pro urep. Typickým příkladem je potenciál tuhých koulí (hard sphere, HS),
uHS(r) =
∞ pro r < σ ,
0 pro r > σ .(39)
Teoreticky výhodným mezistupněm mezi realistickými spojitými potenciály a potenciálem tuhých koulí je tzv.model pravoúhlé potenciálové jámy (square-well),
uSW(r) =
∞ pro r < σ ,
−ε pro σ < r < λσ ,
0 pro r > λσ ,
(40)
kde bezrozměrný parametr λ určuje dosah přitažlivé části potenciálu. V limitě λ→ 1, ε→∞ dostaneme lepkavétuhé koule (sticky hard spheres), v limitě λ→∞, ε→ 0 přitažlivé pozadí (Kacův potenciál).
21
4.1.2 Více atomů
Reálné systémy se skládají z mnoha atomů. Můžeme být v pokušení všechny párové interakcejednoduše sečíst
U(~rN) ≈∑i<j
u(~ri, ~rj) (41)
kde u(~ri, ~rj) = u(|~rj−~rj|) = u(~rij) je párový příspěvek. Symbolem∑
i<j označujeme zkráceněsoučet přes všechny páry, přesněji (pro N atomů)
∑i<j
=N∑j=2
j−1∑i=1
Bohužel, tato tzv. aproximace párové aditivity neplatí přesně. Budeme-li zkoumat trojiceatomů argonu, zjistíme, že odchylka
představuje okrouhle 10 % celé interakce. Zahrnutí těchto tzv. tříčásticových (obecně více-částicových) interakcí je však složité (s výjimkou speciálního typu, tzv. polarizovatelnosti,viz dále), a proto se pravidla postupuje jinak. Párový potenciál nahradíme efektivní verzí,která pro typické konfigurace (např. v kapalině) již v průměrném smyslu všechny mnohočás-ticové příspěvky obsahuje. Takový potenciál pak funguje dobře pro kapalinu, naopak všakpopisuje hůře volný pár částic a potažmo plynnou fázi.
4.2 Ionty
Máme-li částice elektricky nabité, např. ionty s náboji q1 a q2, musíme k interakční energii(např. Lennard-Jonesově) přidat ještě coulombickou interakci. V soustavě SI platí
uCoul(r) =1
4πε0
q1q2
r, (42)
kde ε0 je permitivita vakua. Tento potenciál ubývá s rostoucí vzdáleností velmi pomalu, cožz hlediska simulací (ale i teorie) vede k řadě komplikací.
4.3 Polarizovatelnost
Coulombický potenciál (42) na jednotlivých interakčních centrech nevystihuje celou elektric-kou interakci molekul. Působí-li na atom (molekulu) i elektrické pole o intenzitě ~Ei, dojdek deformaci elektronových obalů a ke vzniku (v prvním přiblížení) indukovaného dipólu ovelikosti
~µi = αSI · ~Ei (43)
Pro sféricky symetrický atom je vektor µi rovnoběžný s ~Ei a tedy polarizovatelnost αSI jeskalár, obecně ale pro nesymetrickou molekulu není ~µi rovnoběžné s ~Ei a αSI je tenzor.
Jednotkou polarizovatelnosti v SI je C m2 V−1 = C2m2 J−1. Obvykle se ale uvádí veličina α = αSI/(4πε0),které se správně v SI říká objem polarizovatelnosti, vzhledem k lenosti a konzervativnosti vědců ale nejčastějiuslyšíte jen „polarizovatelnostÿ. Je totožná s polarizovatelností v soustavě CGS (která je ve škole zakázaná,
22
Obr. 6: Oligopeptid v silovém poli CHARMM22 (vlevo), jeho parciální náboje (v e, uprostřed) amodel vody TIP4P (záporný náboj je posunut z kyslíku do bodu M).
ač teoretiky stále hojně používaná). Její rozměr je m3, nejčastěji se ale uvádí v A3 (1 A = 10−10m). Objempolarizovatelnosti je řádově roven desetině objemu molekuly. Objem polarizovatelnosti v A3 převedeme na SIznásobením faktorem 1.11265 · 10−40.
Protože indukované dipóly budí opět elektrostatické pole, je výsledný potenciál soustavypolarizovatelných molekul dán implicitně (nutno řešit rovnici pro self-konzistentní pole) anení ani párově aditivní, což prodražuje jakékoliv výpočty. Proto se kromě speciálních simu-lačních metod, které jsou již mimo záběr těchto skript, častěji používají efektivní potenciály.
4.4 Víceatomové molekuly
Svět jednoatomových molekul by byl velmi chudý. Zpravidla potřebujeme studovat molekulyod malých (např. voda) po velké (např. proteiny). Základním stavebním prvkem modelů jeatom či přesněji řečeno atomové jádro; někdy (zvláště u starších modelů) se alifatické skupinyjako -CH3 nahrazují jedním „sjednocenýmÿ neboli „rozšířenýmÿ atomem (anglicky „unitedÿnebo „extendedÿ atom – viz obr. 6), výjimečně se naopak přidává pomocné centrum mimoatomy (model vody TIP4P). Počet typů atomů je však větší než chemických prvků, protoženapř. atom uhlíku ve skupině >C=O má jiné vlastnosti než v benzenu. Atomy také zpravidlanesou parciální náboje, viz obr. 6.
Síly mezi atomy se dělí na vazebné (bonded) a nevazebné (nonbonded). Nevazebnésíly působí mezi atomy na různých molekulách a také mezi atomy jedné molekuly, kterénejsou ovlivněny chemickými vazbami. Atomy oddělené více než třemi vazbami v řetězci(např. koncové vodíky v propanu CH3-CH2-CH3) již interagují jen nevazebně, atomy spojenévazbou (symbolicky 1–2) či oddělené dvěma vazbami (1–3) interagují jen vazebnými silami,atomy oddělené třemi vazbami v řetězci (1–4) jsou mezní případ, ve kterém se zpravidlakombinují oba druhy sil.
4.4.1 Nevazebné síly
I mezi atomy schovanými v molekulách působí stejné síly jako mezi volnými atomy: odpudivé,Londonovy přitažlivé a elektrické.
23
Nenabitá část Odpudivé a Londonovy přitažlivé síly se obvykle spojují do Lennard--Jonesova potenciálu (38) nebo exp-6.
Parametry Lennard-Jonesova potenciálu se tabelují jako parametry pro jednotlivé druhyatomů. Rozlišují se přitom atomy s různými chemickými vlastnostmi, např. C v karbonylovéskupině >CO má jiné parametry než C v methylu -CH3. Studujeme-li látky složené z vícedruhů atomů, řekněmě oxid uhličitý CO2, potřebujeme kromě interakce C. .C a O. .O ještěkřížovou interakci C. .O. Tyto parametry počítáme z vhodných kombinačních pravidel.Pro Lennard-Jonesův potenciál se často používá Lorentzovo–Berthelotovo pravidlo
σij =σi + σj
2, εij =
√εiεj,
tj. atomové poloměry jsou aditivní a energie jsou dány geometrickým průměrem. V novějšíchsilových polích se zpravidla používá geometrický průměr i pro σ.
Coulombovy síly Elektrony v molekule jsou rozloženy mnohem složitěji než u sférickysymetrických iontů. Toto rozložení nábojů se snažíme vystihnout pomocí parciálních nábojůumístěných (obvykle) na atomových jádrech, méně často pomocí elektrických dipólů (čivyšších multipólů) umístěných tamtéž. Náboje interagují Coulombovým potenciálem (42).
Parciální náboje můžeme získat z vhodného kvantově chemického softwaru. Standardnímetodou je CHELPG (CHarges from Electrostatic Potentials using a Grid based method):okolo molekuly rozmístíme testovací body a hledáme takové rozmístění parciálních nábojů najádrech, aby pole v testovacích bodech co nejlépe odpovídalo kvantově-chemickému výpočtu.
Abychom ušetřili výpočetní čas, síly pro vzdálenost atom–atom větší než určitá vzdá-lenost se zanedbávají. Protože však jen zanedbání (useknutí, cutoff ) potenciálu by vedloke skoku v potenciálu a tedy k nekonečné síle, což vadí především v molekulové dynamice,potenciál se buď posune o konstantu, aby byl spojitý, nebo se vyhladí nebo se kombinují obapostupy. Obvyklé vzdálenosti, od kterých se zanedbávají disperzní síly, jsou 8–15 A nebo2.5–4 σ. Chyba způsobená zanedbáním se může do značné míry omezit přičtením korekce,kterou zpravidla počítáme za předpokladu, že ve větší vzdálenosti než je cutoff jsou ostatníatomy rozmístěny rovnoměrně.
Výše uvedený trik lze použít v nouzi i pro elektrostatické síly, je však nutno pečlivě kom-binovat useknutí, posun potenciálu a vyhlazení, protože elektrostatické síly jsou silné. O něcolepší metodou je nahrazení polárních skupin (např. >C=O) od určité vzdálenosti bodovýmdipólem, který se pak již lépe spojitě „usekáváÿ; silové pole je totiž obvykle navrženo poskupinách, které jsou elektricky neutrální. Nyní ale nejčastěji používají speciální metody,Ewaldova sumace a metoda reakčního pole, které řeší problém elektrostatické interakce po-mocí dobře definovaných aproximací.
Ewaldova sumace je založena na sečtení interakcí v periodických okrajových podmínkáchpřes všechny periodické obrazy nábojů v simulační buňce. To ale matematicky znamená sečísttrojrozměrnou nekonečnou řadu, která ještě ke všemu pomalu (a jen relativně) konverguje.Matematickým trikem se tato řada převede na dvě řady, které však konvergují mnohemrychleji.
Tento trik se dá interpretovat tak, že k bodové náboje odstíníme Gaussovým rozložením stejného nábojeopačného znaménka. V určité vzdálenosti (která se volí menší než polovina hrany simulační buňky) je náboj prak-ticky odstíněn. Interakce stíněných nábojů je dána místo funkce qiqj/(4πε0rij) funkcí qiqjerfc(αrij)/(4πε0rij),
24
1
2 3
4
φ
12
3 4φ
Obr. 7: Vlastní (vlevo) a nevlastní (vpravo) torze
která rychle ubývá se vzdáleností; α je konstanta. Zbývá započíst interakci Gaussovsky rozmytých nábojů. To sepočítá v k-prostoru pomocí Fourierovy transformace nábojového rozložení. Pro větší systémy se přitom dosáhnevyšší efektivity využitím mřížky (particle mesh) a algoritmu tzv. rychlé Fourierovy transformace.
Metoda reakčního pole je vhodná jen pro dipolární systémy. Interakci vybraného di-pólu s ostatními dipóly počítáme plně jen do určitého „cutoffuÿ. Všechny vzdálenější dipólynahradíme spojitým dielektrikem, které se elektrickým polem vybraného dipólu polarizuje.
4.4.2 Vazebné síly
Vazby, pokud nejsou v simulaci pevné, se obvykle aproximují harmonickým potenciálem
Ubond(r) =Kbond
2(r − r0)2 (44)
kde r je vzdálenost atomů, r0 její rovnovážná hodnota a Kbond silová konstanta (tuhostvazby).
Podobně se řeší vazebné úhly, tedy síly 1-3:
Uangle(α) =Kangle
2(α− α0)2 . (45)
Síly 1–4 jsou dvojí. Vlastní torze (proper torsions) tvaru
Udih(φ) =Kdih
2cos(nφ) , (46)
kde φ je tzv. dihedrální úhel (viz obr. 7) a n je číslo, např. pro torzní potenciál C–C–C–Cv butanu je n = 3, případně se sčítá několik členů pro n = 0, 1, 2, 3.
Nevlastní torze (improper torsions) slouží ke „zpevněníÿ skupin s sp2 hybridizací jakoH2C=O, kde uhlík leží v rovině s ostatními atomy, případně k zajištění tetraedrické geometrietří vazeb okolo skupiny CH, jestliže je tato reprezentována jedním interakčním centrem (vizCH1E obr. 6):
Uimprop(φ) =Kimprop
2(φ− φ0)2 (47)
Popis interakce mezi molekulami vychází z toho, že molekula je složena z atomů. Tytoatomy nejsou v molekule „ztracenyÿ, ale atom jedné molekuly může interagovat s atomy jinémolekuly. Molekula se tedy skládá z interakčních center a mezimolekulární potenciál lze psátve tvaru
u(1, 2) =∑a∈1
∑b∈2
uab(|~r2,b − ~r1,a|) (48)
25
kde symbolické proměnné (1,2) označují soubor souřadnic molekul 1 a 2 a ~ri,a značí polohovývektor interakčního centra a na molekule i. Za potenciál uab(r), který závisí pouze na vzdále-nosti mezi centry, se volí některý z jednoduchých párových potenciálů; obvyklá je kombinacecoulombické (42) a neutrální (např. (38)) interakce.
Interakční centra nejčastěji odpovídají přímo jednotlivým atomovým jádrům v molekule,v závislosti na úrovni aproximace interakcí mohou však odpovídat i celým skupinám atomů(„unitedÿ nebo „extendedÿ atom), případně naopak mohou být v molekule pomocná cen-tra mimo jádra. Pro malé molekuly se obvykle rozložení těchto center považuje za fixní amolekula nebude tedy mít žádné vnitřní stupně volnosti (je „tuháÿ). Potom jednu molekulupopisujeme kromě polohy referenčního bodu („centraÿ, obvykle těžiště) ještě orientací, cožje pro lineární molekuly vektor osy molekuly (či ekvivalentně dva úhly ve sférických sou-řadnicích) a pro obecné molekuly buď orientační matice nebo tři Eulerovy úhly, případněkvaternion.
4.5 Vnější síly
Působí-li na systém vnější pole (gravitační, elektromagnetické), umíme obvykle napsat téžpříslušnou interakční energii. Stěny nádoby jsou však také vnějším polem a podobně lzechápat i kapalinu uzavřenou v póru, molekuly adsorbované na povrchu tuhé látky, apod.Potřebujeme proto i potenciálový model materiálu, s kterým náš systém intraguje. Nejpřes-nější ale také nejnáročnější je složit stěnu či pór z jednotlivých atomů. Tak se modelujípředevším různé pórézní materiály, např. zeolity. Je-li danou stěnou krystal, musíme rozlišitrůzné krystalové roviny. Někdy stačí jen jedna vrstva atomů, ale i tak může jejich početsnadno překročit tisícovku. Na druhém konci zjednodušování je, podobně jako model tuhýchkoulí pro atomy, model tuhé bezstrukturní stěny. Je-li rovina tuhé stěny v kartézských sou-řadnicích definována rovnicí z = 0, je potenciální energie částice o souřadnicích ~r = (x, y, z)dána vztahem
Utuhá stěna(~r) =
∞, pro z < 0,
0 pro z ≥ 0(49)
Mezistupněm vhodným pro realističtější potenciály je model měkké bezstrukturní stěny.Představme si, že v poloprostoru z > 0 jsou rovnoměrně rozmístěny částice s číselnou husto-tou ρ. Atomy vně stěny nechť interagují s částicemi stěny potenciálem u(r). Odhlédneme-liod jednotlivých atomů a budeme-li látku stěny považovat za kontinuum, dostaneme efektivnípotenciál částice-stěna po integraci přes poloprostor z > 0:
Uměkká stěna(~r) =∫z′>0
u(~r + ~r′)d~r′ =∫ ∞−∞
dx′∫ ∞−∞
dy′∫ ∞
0dz′u(~r + ~r′) (50)
Speciálně pro Lennard-Jonesův potenciál u = uLJ, viz (38), který je typu 12–6, dostanemepo integraci potenciál typu 9–3:
ULJ−stěna(~r) = περ
[445
(σz
)9− 2
3
(σz
)3]
(51)
Integrál v (50) vypočteme tak, že nejprve integrujeme přes roviny z′ = const v polárních souřadnicích r′, φ′
Obr. 8: Schéma Isingova modelu v rovině (vlevo) a mřížkového modelu polymeru (vpravo)
5 Klasické mřížkové modely
Mřížkové modely se používají v mnoha oblastech fyziky i chemie, a proto je obtížné podatvyčerpávající definici. Podstatné je, že prostorové souřadnice jsou diskrétní a jejich počet jekonečný. Obvykle se vychází z nějaké pravidelné (krystalové) mřížky, např. kubické. Každémuvrcholu mřížky označenému indexem i je přiřazena proměnná si, případně n-tice proměnných(vektor). Ty mohou nabývat hodnot z jisté množiny, která může být diskrétní i spojitá. Modelje definován interakční (též konfigurační) energií U(si) (někdy se nazývá Hamiltonián).
5.1 Isingův model
Původní motivací pro vytvoření tohoto modelu bylo studium feromagnetismu, má však řaduaplikací i v dalších oblastech. S každým vrcholem mřížky je spojena skalární proměnná si,která se nazývá spin a která nabývá dvou hodnot +1 nebo −1 (viz obr. 8).
Konfigurační energie je dána vztahem
U = −J∑<i,j>
sisj + h∑i
si (52)
kde J je síla interakce (pro feromagnet je J > 0) a h je vnější magnetické pole. První suma jepřes všechny dvojice nejbližších sousedů <i, j> (hrany mřížky) a druhá suma přes všechnyvrcholy. Isingův model je exaktně řešitelný v jedné a pro h = 0 i ve dvou dimenzích, vefyzikálně zajímavých třech dimenzích musíme použít buď aproximativní řešení nebo simulace.
5.2 Mřížkový plyn
Je to nejhrubší model pro systém N klasicky interagujících částic. Celý uvažovaný prostor, aťuž plochu či objem, si rozdělíme na malé buňky a předpokládáme, že částice leží ve středechbuněk. Skalární proměnná ni spojená s vrcholem i má význam počtu částic v buňce. Odpu-divé síly mezi částicemi popíšeme tak, že budeme uvažovat nejvýše jednu částici v buňce, tedyni ∈ 0, 1 (víc se jich tam nevejde), přitažlivou interakci pak můžeme omezit jen na sou-sedící buňky obsazené částicemi. Uvažujeme-li mřížkový plyn v grandkanonickém souboru,tedy s proměnným počtem částic a chemickým potenciálem µ, je jeho interakční energie
27
nízká teplota rychle ochlazený systém vysoká teplota kritický bod
Obr. 9: Typické konfigurace dvoudimenzionálního Isingova modelu
rovnaU = ε
∑<i,j>
ninj − µ∑i
ni (53)
kde ε je energie páru částic. Ze substituce ni = (1 + si)/2 vyplývá, že (53) je, až na definicikonstant, ekvivalentní Isingově energii (52).
5.3 Binární slitina
Model slitiny, jejíž dva atomy jsou natolik podobné, že se mohou zaměňovat v krystalovémřížce, lze popsat energií
U = −∑<i,j>
εkikj +∑i
µki , ki ∈ ,
kde ε , , ε , a ε , jsou energie interakce sousedních atomů a µ , µ jejich chemicképotenciály. Model je ekvivalentní mřížkovému plynu (ni = 0 ∼ ki = , ni = 1 ∼ ki = ).
Existují různá zobecnění Isingova modelu. Spinová proměnná si může nabývat více neždvou diskrétních hodnot (Pottsovy modely), nebo může být i spojitá, např. omezená nasféru s2
i = 1 (spojitý Heisenbergův model). Významnou aplikací jsou kvantové teorie pole,kde obor hodnot spinových proměnných je grupa transformací, např. SU(2), SO(3) aj.
5.4 Model polymeru
Jiným příkladem klasického mřížkového modelu je model polymeru, obr. 8. Zde se obvyklepracuje s hranami mřížky, které buď jsou nebo nejsou obsazeny článkem polymerního řetězce.Podmínkou je, že se řetězec nesmí protínat. Pokud články již nijak neinteragují a neuvažu-jeme rozvětvené řetězce, je model ekvivalentní tzv. náhodné procházce bez protínání (randomself-avoiding walk).
6 Molekulová dynamika
V metodě molekulové dynamiky (MD) sledujeme vývoj systému složeného z atomů/molekulv „reálnémÿ čase. Podle typu modelu můžeme rozlišit různé metody
28
• Tuhé koule a další modely s nespojitým potenciálem. Tyto částice se pohybují prosto-rem rovnoměrně přímočaře, dokud se nesrazí. V okamžiku nárazu změní směr (podlezákonů pružného rázu: celková energie i hybnost se nemění).
Algoritmus takové simulace je řízen událostmi. Vytvoříme si tabulku všech možnýchsrážek koulí, případně dalších možných událostí (měření v rovnoměrných časovýchintervalech ap.). Vybereme událost, ke které dojde nejdříve a provedeme ji; zpravidlaje pak nutné přepočítat některé další události (molekuly změnily směr).
• Klasická molekulová dynamika se spojitým potenciálem. Těmito metodami se budemev dalším textu zabývat podrobněji.
• Dynamika s náhodnými silami. Některé stupně volnosti (rozpouštědlo, nahrazení většískupiny atomů jedním „coarse-grainedÿ atomem) nahradíme Gaussovsky náhodnými„šťouchanciÿ. Systém jinak simulujeme stejně jako v klasické MD. Aby šťouchancesystém nezahřívaly, přidáváme tření (molekuly se rovnoměrně zpomalují); uděláme-lito správně, dostaneme systém při zadané teplotě T (Langevinův termostat). Postup setéž nazývá Brownowská dynamika. Ztrácíme tím ale zcela hydrodynamické chování(proudění kapaliny; např. nelze stanovit viskozitu) a nezachovává se hybnost.
Složitější varianta zvaná disipativní částicová dynamika šťouchá náhodně symetrickyvždy do páru částic. Zachovává se tak hybnost a lze studovat i hydrodynamické jevy.Tyto metody jsou vhodné pro velmi velké systémy, kdy si nemůžeme dovolit atomárnírozlišení.
• Metoda dráhového integrálu (path integral) umožňuje správně popsat kvantové cho-vání jader. Např. jsou správně popsány nulové kmity lehkých atomů (vodíku). Pů-vodní systém je však stále popsán potenciální energií (silovým polem) na Bornově––Oppenheimerově úrovni. Časová náročnost značně stoupá.
• Kvantové simulace typu Car–Parrinello (a odvozené) nepotřebují meziatomový poten-ciál, protože integrují v čase vývoj vlnové funkce. I přes množství aproximací jsou velmičasově náročné.
6.1 Klasická molekulová dynamika
Uvažujme atomární systém popsaný vektory poloh atomů ~ri, i = 1, . . . , N ; těchto 3N číselzapíšeme úsporně jako ~rN . Interakce jsou popsány mezimolekulovým potenciálem nebolisilovým polem, U(~rN).
Abychom mohli studovat vývoj systému, potřebujeme síly. Ty jsou dány gradientempotenciálu
~fi = −∂U(~rN)∂~ri
i = 1, . . . , N (54)
To je celkem N vektorů, tj. celkem 3N čísel.
Párové síly jsou dány součtem párových příspěvků přes všechny dvojice molekul:
U =∑i<j
u(rij)
29
Síla na částici i je pak
~fi =N∑j=1j 6=i
~fji ≡ −N∑j=1j 6=i
du(rji)drji
~rjirji
kde používáme značení ~rij = ~rj − ~ri, rij = |~rij |. Síla na atom i je tedy rovna součtu párových sil ode všechostatních atomů. Zároveň platí ~fji = −~fij (Newtonův zákon akce a reakce).
Podle druhého Newtonova zákona se částice, na kterou působí síla ~fi, se pohybuje sezrychlením
d2~r
dt2≡ ~ri =
~fimi
, i = 1, . . . , N
Abychom tuto soustavu obyčejných diferenciálních rovnic druhého řádu mohli řešit, potře-bujeme znát počáteční podmínky. To jsou polohy ~ri(0) a rychlosti ~ri(0) v počátečním čase(bez újmy na obecnosti pokládáme t = 0).
V metodě molekulové dynamiky pro spojité potenciály řešíme tuto rovnici metodou ko-nečných diferencí. Výsledkem je tzv. trajektorie, tj. posloupnost konfigurací ~ri(t) a rychlostí~ri(t) pro diskrétní časy t = kh, kde k je celé číslo a h integrační krok (proměnný kroknebudeme uvažovat).
Numerické matematika vyvinula celou řadu integračních metod:
• Rungeovy–Kuttovy metody umožňují změnu integračního kroku a jsou vyššího řádu,tj. velmi přesné. Nejsou však časově reverzibilní (viz dále) a potřebují několik výpočtůpravé strany (tj. sil) na integrační krok, takže nejsou efektivní. Prakticky se pro MDnepoužívají.
• Metody typu prediktor–korektor vycházejí ze znalosti trajektorie v několika předcho-zích časech (historie). Z této historie predikují následující hodnotu, která je upřesněnakorekcí, jež je již založena na výpočtu pravé strany. V MD se někdy používá variantazvaná Gearova metoda.
• Symplektické metody jsou nejen časově reverzibilní, ale pro dynamické systémy za-jišťují, že celková energie8 je aproximována s určitou chybou, která v čase neroste akterou lze zkrácením integačního kroku zmenšit. Takové chování je pro MD ideální.Nejjednodušší i nejčastěji používanou symplektickou metodou je Verletova metoda (aekvivalentní varianty).
6.1.1 Verletova metoda
Z Taylorova rozvoje vyplývá následující vzorec pro výpočet druhé derivace
~ri(t) =~ri(t− h)− 2~ri(t) + ~ri(t+ h)
h2+O(h2)
8Klasická mechanika zavádí tzv. Hamiltonův formalismus, ve kterém je systém popsán polohami ~rN a(sdruženými) hybnostmi ~pN a vyvíjí se podle tzv. Hamiltonových rovnic, které jsou ekvivalentní Newtono-vým pohybovým rovnicím. Podle těchto rovnic se zachovává objem elementu fázového prostoru d~rNd~pN .Symplektické metody aproximují toto zachování s určitou omezenou chybou, která v čase neroste.
30
Obr. 10: Metoda leap-frog
Symbol O(h2) znamená, že chyba je řádově h2. Po reorganizaci a zanedbání chyby máme
~ri(t+ h) = 2~ri(t)− ~ri(t− h) + h2~fi(t)mi
(55)
To je již návod pro výpočet poloh částic v čase t + h, známe-li jejich polohy v časech t at − h a síly v čase t. Není to však řešení počáteční úlohy, jak jsme si stanovili výše. Pokudznáme v čase t = 0 polohy a rychlosti, musíme si ještě dopočítat odpovídající polohu v časet = −h, např. Taylorovým rozvojem
~ri(−h) = ~ri(0)− h~ri(0) +h2
2
~fi(0)mi
+O(h3)
nebo s menší přesností bez členu se silami. Ze znalosti poloh v časech t = 0 a t = −h pakvypočteme polohy v čase t = h, pak t = 2h, atd.
Ve vzorci však bohužel nevystupují rychlosti, které potřebujeme např. pro výpočet ki-netické energie a potažmo teploty. Můžeme si je (po provedení Verletova kroku) dopočítattakto
~ri(t) =~ri(t+ h)− ~ri(t− h)
2hExistují však i jiné možnosti.
Verletovu metodu často uvidíte v některém jiném ekvivalentním tvaru. Trajektorie je totožná, někdy jsou všakmírně (tj. až na chybu O(h2)) odlišné rychlosti. Nejnázornější je tzv. metoda leap-frog. Protože rychlost je dráha(tj. změna polohy) za jednotku času (totiž h), můžeme pro průměrnou rychlost v intervalu (t, t + h) napsat
~v(t,t+h) =~r(t + h)− ~r(t)
h
Pokud se tato rychlost nemění příliš prudce, je přibližně rovna rychlosti v polovině tohoto intervalu,
~v(t + h/2) =~r(t + h)− ~r(t)
h
Podobně zrychlení je změna rychlosti za jednotku času, tedy
~a(t) =~v(t + h/2)− ~v(t− h/2)
h
31
což je ovšem rovno ~f/m. Algoritmus jednoho kroku metody leap-frog je proto dán aplikací dvou dosazení
~v(t + h/2) = ~v(t− h/2) + h~f
m~r(t + h) = ~r(t) + ~v(t + h/2)h
Rychlosti zde známe v polovičních časech. Rychlost v čase t dostaneme např. jako průměr
~ri(t) =~v(t + h/2) + ~v(t− h/2)
2
z čehož lze snadno odvodit ekvivalenci s Verletovou metodou.
Rovnice (55) se nezmění, jestliže zaměníme t→ −t a h→ −h. Verletova metoda je tedyčasově reverzibilní. To mj. znamená, že celková energie nemá žádný trend (drift), tj. sys-tematicky neroste ani neklesá. To neznamená, že je celková energie úplně přesně konstantní(jak by měla správně být): náhodně kolísá s přesností danou výrazem O(h2), který jsmepři odvození zanedbali9. Pokud aplikujeme Verletovu metodu na oběh planety okolo Slunce,bude obíhat v podstatě správně po elipse, jejíž parametry se nemění (energie je konstantní),ale tato elipsa se bude poněkud stáčet, což neodpovídá realitě. Blíže chemickým aplikacímje integrace harmonického oscilátoru. Energie je opět pěkně konstantní, ale frekvence se odpřesné poněkud liší. V reálných simulacích nám mírná změna frekvence nevadí, ale změnaenergie by nám vadila. Proto je Verletova metoda (a její klony) tak populární.
Nevýhodou Verletovy metody je, že rychlost v čase t známe až po provedení kroku.Nemůžeme proto (bez dalších úprav) integrovat rovnice s pravou stranou obsahující rychlosti;tohoto tvaru jsou např. rovnice popisující rotaci tuhého tělesa či Noséův–Hooverův termostat.
Existují symplektické metody vyššího řádu. Větší význam má však jiná kategorie me-tod na podobném principu: metody několikanásobného kroku (multiple timestep methods).Pokud je část silového pole odpovědná za rychlé pohyby (např. vibrace vazeb) rychlá i k vý-počtu, integrujeme tyto rychlé pohyby s krátkým integračním krokem. Jednou za několikrychlých kroků spočteme všechny síly (i ty, které se mění pomalu, ale jejichž výpočet jenáročnější) a provedeme „superkrokÿ.
6.1.2 Gearovy metody
Tyto metody vycházejí ze znalosti historie v několika předchozích časech, ~r(t)N , ~r(t − h)N ,~r(t−2h)N , . . . Na základě těchto hodnot predikujeme pomocí polynomu10 hodnotu ~rp(t+h)N .
Hodnoty ~r(t)N , ~r(t − h)N , ~r(t − 2h)N jsou obecně nepřesné a chyby se po provedenípredikce propagují do dalších kroků; pokud bychom jenom predikovali a nepoužili vůbecpravou stranu (síly), chyby by exponenciálně rychle rostly. Stejně jako jsme predikovali ~r(t+h)N , můžeme predikovat i druhou derivaci ~rp(t+h)N . Tu ovšem můžeme nezávisle vypočítatze sil, které spočteme v predikovaných polohách: ~r(t+h)N = f(~rp(t+h)N)/m. Tušíme, že totovypočtené zrychlení bude asi lepší. Chyba, které bychom se dopustili, kdybychom zrychlenínevypočetli, je rovna rozdílu E = ~r(t+ h)N − ~rp(t+ h)N . Trik vedoucí k fungující integračnímetodě je dán tím, že chybu E (jiné kritérium nepřesnosti kroku nemáme!) ponásobímeurčitými koeficienty (naleznete je ve speciální literatuře) a přičteme k historii tak, abychom
9Jak jsme se již zmínili, je to dokonce ještě lepší, chyba je omezená, protože metoda je symplektická10Zpravidla se používá jiné, avšak zcela ekvivalentní vyjádření – místo historie známe kromě ~r(t)N ještě
derivace ~r(t)N , ~r(t)N , atd.; Taylorův rozvoj v bodech t, t− h, . . . je právě výše uvedená historie.
32
Obr. 11: Vývoj celkové energie v závislosti na metodě integrace a délce kroku h = 0.005 ps ( ),h = 0.01 ps ( ) a h = 0.02 ps ( ); některé křivky jsou pro přehlednost posunuty o vhodnýnásobek 100 K. Výsledky simulace 216 atomů Lennard-Jonesova modelu argonu při teplotě přibližně150 K a hustotě 1344 kg m−3.
33
dostali „správnéÿ zrychlení (to je krok zvaný korektor). To lze udělat mnoha způsoby. Volímetakový způsob, aby chyby, které se i po korekci propagují z kroku do kroku, nerostly, alenaopak po několika krocích vymřely.
Metoda využívající pouze tři body historie je „náhodouÿ ekvivalentní Verletově metodě.Metody vyššího řádu jsou přesnější (např. při integraci pohybu planety nedávají takovýposun perihelia, frekvence harmonického oscilátoru je přesnější), bohužel však nejsou časověreverzibilní, a proto celková energie buď systematicky klesá nebo roste. Výhodou metody jevšak, že ji lze bez problému použít i pro pravou stranu používající rychlosti.
6.2 Volba integrátoru a integrační krok
Obecně lze vzhledem k jednoduchosti a dobrému zachování energie doporučit Verletovumetodu. Pokud nevadí drift energie (např. proto, že stejně používáme termostat) a naopakmáme složitější systém s pravou stranou obsahující rychlosti, lze volit Gearovu integraci.
V obou případech je nutno zvolit vhodný integrační krok. Příliš krátký krok je neefektivní– nejpomalejší část simulace, výpočet sil, se provádí zbytečně často. Příliš dlouhý krok vedek velkým chybám, které mohou vést i ke krachu simulace (dojde k překryvu částic, síla vzrostenad únosnou mez a v příštím kroku se částice posune opět do oblasti překryvu nebo někamdo velké vzdálenosti). Verletova metoda je přitom mnohem odolnější. Základním kritériemvhodné délky kroku je přitom právě přesnost zachování celkové energie. Čím lehčí atom, tímse pohybuje rychleji a tím musí být krok kratší. Pro systémy obsahující vodíky se volí krokokolo 1 fs, pro modely bez vodíků (např. kapalný argon) stačí několikanásobek.
6.3 Teplota
Výše uvedená simulace dává mikrokanonický soubor s konstantní celkovou energií. (Kromětoho se mohou zachovávat i další integrály pohybu, např. hybnost.) Teplotu stanovíme z ki-netické energie pomocí ekvipartičního principu
Tkin =Ekin12kBf
(56)
kde f je počet stupňů volnosti. Ten je roven 3N minus počet zachovávajících se veličin(např. pro kapalinu v periodických okrajových podmínkách je f = 3N − 4: odčítáme 1 zazachování energie a 3 za zachování hybnosti; moment hybnosti se v periodických okrajovýchpodmínkách nezachovává).
Teplota Tkin fluktuuje v průběhu simulace. Čím větší systém, tím jsou fluktuace menší(typická velikost odchylky od průměru je úměrná N−1/2). Teplotu pak spočteme jako časovoustřední hodnotu.
Příklad. Kolik integrálů pohybu (stupňů volnosti) se zachovává přisimulaci tekutiny v nekonečně dlouhém válcovém póru, je-li systém pe-riodický ve směru osy válce?
Řešení. Systém je translačně invariantní ve směru osy válce, bude se tedy zachovávat hybnost(v simulaci ji nastavíme na 0). Dále je systém rotačne invariantní kolem téže osy, bude se tedyzachovávat moment hybnosti okolo této osy (v simulaci ho nastavíme na 0). Dále se zachovávácelková energie. Pro výpočet teploty použijeme f = 3N − 3.
34
6.4 Molekulová dynamika za konstantní teploty
Simulace za konstantní energie (mikrokanonické) jsou nepraktické. Abychom měli systémza konstantní teploty T , musíme přidat termostat. Metody můžeme rozdělit na ty, kterédávají skutečný kanonický soubor, a metody přibližné. Pro přesné kanonické metody platí,že pravděpodobnost nalezení konfigurace s celkovou energií Etot je úměrná Boltzmannovufaktoru exp(−Etot/kBT ). Z toho mj. plyne, jak fluktuuje skutečná (kinetická) teplota okoloteploty termostatu T 11. Přibližné metody dávají sice v průměru správnou teplotu, 〈Tkin〉 = T ,ale fluktuace jsou jiné.
Berendsenova metoda Nejznámější přibližnou metodou je metoda přeškálování rych-lostí. Protože kinetická teplota je úměrná dvojmoci rychlostí, můžeme po dokončení krokuznásobit všechny rychlosti faktorem (T/Tkin)1/2, tj.
~vi,new = ~vi(T/Tkin)1/2
Kinetická teplota spočtená po znásobení je rovna T . Při dalším vývoji by se ovšem opětodchýlila, a tak je je nutno násobit po každém kroku. To však (zvláště pro malé systémy)narušuje pohybové rovnice. Proto se provede podobný krok správným směrem – je-li skutečnáteplota příliš vysoká, rychlosti se sníží, ale jen o trošku
~vi,new = ~vi(T/Tkin)q, q < 1/2
Po mnoha těchto krocích teplota začne fluktuovat okolo nastavené. Metodu lze přepsat12 dodiferenciálního tvaru, který se nazývá Berendsenova metoda
~ri =~fimi
− q
h
Tkin − TT
(57)
Člen s q/h je vlastně tření. Je-li teplota příliš vysoká, částice brzdí, a naopak.Berendsenova metoda je schopna rychle ochladit nerovnovážný systém, trpí však jedním
neduhem. Budeme-li simulovat klastr (malou kapku) ve vakuu (bez toho, abychom uměledrželi moment hybnosti na nule), dostaneme časem nikoliv kapku dané teploty, ale rychlerotující krystalek o nulové teplotě (flying icecube) – kinetická energie sice formálně odpovídádané teplotě, ale je schovaná jen v jednom pohybu. V periodických okrajových podmínkáchje tento artefakt méně významný, ale stejně se doporučuje udržovat q co nejmenší.
Člen q/h lze také zapsat jako 1/(2τ), kde τ je typický korelační čas. Berensenův termostat se chová jakoreálný termostat, kterému také trvá nějakou dobu, než se teplota zkumavky ustálí. Abychom vztah pro τ odvodili,zanedbejme síly v (57), znásobme (57) výrazem mi~ri a sečtěme přes všech N částic. Dostaneme
12Tkin = − q
h
Tkin − TT
Tkin
Po převodu na ∆T12
∆T = − qh
∆TTTkin ≈ −
q
h∆T
ježto T ≈ Tkin. Řešení je ∆T = exp(−t/τ), τ = h/(2q), tj. teplota Tkin se přibližuje nastavené T exponenciálněs časovou konstantou (korelačním časem) τ . Vzhledem k zanedbání sil platí toto jen pro ideální plyn, v reálnémsystému je to však stále dobrá řádová aproximace (platí tím hůř, čím víc se liší tepelná kapacita od ideální).
11Z těchto fluktuací lze např. spočítat tepelnou kapacitu systému.12Spočtěte rozdíl ~vi,new−~vi během jednoho kroku h, vyjádřete ho pomocí ∆T = Tkin−T a předpokládejte,
že ∆T T .
35
Tabulka 1: Srovnání termostatůNosé–Hoover⊕ kanonický⊕ velmi kvalitní⊕ vhodný i pro malé systémy
(po rozšíření na N-H řetězec)
oscilace, decoupling(nutno pečlivě nastavit τ)
horší pro start pohybové rovnice s rychlostí
Berendsen⊕ jednoduchý⊕ exponenciální relaxace
(tj. vhodný i pro start)
flying icecube nekanonický velmi špatný pro malé systémy
Andersenova metoda Přesný kanonický soubor dávají metody založené na přiřazenírychlosti z Maxwellova–Boltzmannova rozdělení rychlostí. V Andersenově metodě si (občas)zvolíme náhodně částici, o její původní rychlost se nezajímáme a vybereme novou rychlostz Maxwellova–Boltzmannova rozdělení,
π(xi) =1
σ√
2πexp
(− x2
i
2σ2
), σ2 =
⟨x2i
⟩=kBT
mi
To lze interpretovat, jako kdyby se dotčená částice vykoupala v termostatu; v pravděpodob-nostním smyslu má teplotu termostatu. Můžeme také toto přiřadit nové rychlosti jednou začas pro všechny částice najednou (Maxwellova metoda).
I když tato metoda dává správný kanonický soubor a vzorkuje i hybnost (tedy f = 3Nbez ohledu na okrajové podmínky), má jednu vadu – nezachovává kinetiku. Veličiny jakodifuzivita nebo viskozita nelze takto počítat. Také při výpočtu počtu stupňů volnosti simusíme uvědomit, že integrály pohybu se nezachovávají.
Variantou s malými „šťouchanciÿ je již zmíněný Langevinův termostat.Všechny tyto metody relaxují k nastavené teplotě exponenciálně – jako reálný termostat.
Nosé–Hoover Často používaný kanonický termostat, který zároveň nemění kinetické vlast-nosti. Hlavní myšlenkou je přidání další dynamické proměnné (stupně volnosti) s k pohy-bovým rovnicím. Tato proměnná má rovněž „rychlostÿ s, a proto i kinetickou energii Ms
2 s2
i potenciální energii, která se rovná −fkBT ln s (T je teplota a f počet stupňů volnosti).Pohybové rovnice vyjádřené (z důvodů, které zde nelze vysvětlit) pomocí logaritmu ξ = ln sjsou
~ri =~fimi
− ~riξ
ξ =
(Tkin
T− 1
)τ−2
kde τ je časová konstanta termostatu je
τ =
√Ms
fkBT
36
0 1 2
t/ps
0
500
1000
T/K
Berendsen
0 1 2
t/ps
0
500
1000
T/K
Berendsen CM
0 1 2
t/ps
0
500
1000
T/K
Berendsen IN
0 1 2
t/ps
0
500
1000
T/K
Maxwell CM
0 1 2
t/ps
0
500
1000
T/K
Nose
0 1 2
t/ps
0
500
1000
T/K
Andersen CM
Obr. 12: Srovnání termostatů na příkladu simulace SPC/E vody startované z „kubického krystaluÿ(viz vpravo). CM znamená, že se termostatují jen translace (pohyb těžiště), IN znamená, že setermostatují jen rotace. — celková kinetická teplota, — jen z rychlosti těžiště, — jen z rotací
Celková energie tohoto rozšířeného systému se zachovává, nás však ve výsledku proměnnás (resp. ξ) nezajímá, a proto vlastně sledujeme střední hodnotu přes všechny hodnoty pro-měnné s. Lze ukázat, že (pro ergodický systém) vznikne kanonický soubor, přičemž správněkanonicky jsou i veličiny obsahující rychlosti.
Rozdíl oproti ostatním metodám je, že (po zjednodušení) je výsledná rovnice druhéhořádu (jako harmonický oscilátor). Proto při startu z velmi nerovnovážné konfigurace docházík oscilacím a teplota nerelaxuje dost rychle k nastavené. To, jak metoda pracuje, je citlivé nastanovení parametru (času oscilace) τ . Je-li příliš dlouhý, neinteraguje pomocná proměnná s(ξ) se zbytkem systému, a oscilace se dlouho drží. Výhodnější je mít τ co nejkratší, ale tak,abychom nemuseli zkracovat integrační krok.
Existuje i kanonická varianta Berendsenova termostatu
6.5 Dynamika s vazbami
Pokud při integraci pohybových rovnic požadujeme, aby se zachovávaly nějaké veličiny (např.vzdálenosti atomů, tj. délky chemických vazeb), mluvíme o dynamice s vazbami (constraintdynamics. Případ, kdy fixujeme délky vazeb, je asi nejčastější. Alternativou jsou vazby po-psané (nejčastěji) harmonickým potenciálem. Důvodem k fixování vazeb je několik:
• Vazby vibrují rychle, potřebovali bychom tedy příliš krátký integrační krok. (Tentoproblém lze nicméně řešit metodami několikanásobného kroku.)
• Vazby vibrují rychle, a proto je přenos energie mezi vibračními a ostatními stupnivolnosti pomalý. (Lze napravit např. Andersonovým termostatem, tím však ztrácímečasový vývoj.)
37
~r(t− h)
~r(t)~r(t+ h)
~rSHAKE(t+ h)
λ~r(t)
Obr. 13: Algoritmus SHAKE pro matematické kyvadlo
• Rychle vibrující jsou v realitě kvantovány (a často jsou v nulovém bodu), takže inte-grace klasickou mechanikou stejně nemá smysl.
Představme si, že máme klasický systém s vibrujícími vazbami a vazebnými úhly a žezvyšujeme tuhost harmonických potenciálů (pružin). Tím dochází jednak k zvyšování frek-vencí, jednak k zmenšování amplitudy vibrací. Avšak limita pro nekonečně velké tuhostinení rovna systému s pevnými vazbami! Uvažujme například řetězec C–C–C–C (united--atom model butanu) s tím, že dihedrální potenciál je nulový. Potom (v limitě nekonečněvelké tuhosti) je rozdělení dihedrálních úhlů uniformní, molekulu najdeme se stejnou prav-děpodobností v konformaci syn i anti. Nikoli však, pokud zafixujeme jak vazby, tak úhly,zde je větší pravděpodobnost syn konformace. Hlavní potíž je zde s fixováním úhlů; dá seukázat, že u většiny chemických systémů dává fixace vazeb (s ohebnými úhly) výsledky máloodlišné od limity pro nekonečné tuhosti. Ostatně, fixovat úhly není dobré z fyzikálních dů-vodů, molekula je pak nerealisticky rigidní; výjimkou jsou tuhé molekuly např. vody (fixacíobou vazeb i úhlu dostaneme tuhé těleso).
6.5.1 SHAKE
Nejběžnější metodou pro integraci pohybových rovnic s vazbami je algoritmus SHAKE. Jezaložen na Verletově integraci. Probereme si ho nejprve na příkladu matematického kyvadladélce l (hmotný bod na neroztažitelné niti), viz obr. 13.
Představme si, že použijeme rovnici (55) pro jeden krok integrace bez ohledu na přítom-nost vazeb, ale s tím, že v předchozích krocích byly vazby splněny; toto označíme indexemVerlet. V čase t + h vazby splněny nebudou. Splněny by byly, kdybychom od vnější síly ~fpůsobící na závaží odečetli odstředivou sílu. Pak bychom mohli vlákno přestřihnout, protožesíla působící podél vlákna by byla přesně nahrazena odstředivou silou:
~r(t+ h) = ~rVerlet(t+ h)− h2
m~fc(t) ≡ 2~r(t)− ~r(t− h) + h2
~f(t)− ~fc(t)m
(58)
Fiktivní sílu ~fc(t) neznáme, známe však její směr, který je rovnoběžný s vláknem, tj. ~r(t),což zapíšeme takto:
h2 ~fc(t)m
= λ~r(t) (59)
38
kde λ je neznámé. Spočteme ho z podmínky, že v čase t+h budou mít vazby správné délky.Pro naše kyvadlo tedy
|~r(t+ h)| = |~r(t)| = l (60)
z čehož po umocnění na druhou vypočteme, zanedbávajíce členy vyššího řádu,
λ =|~rVerlet(t+ h)|2 − |~r(t)|2
2~rVerlet(t+ h) · ~r(t)(61)
Pro dva atomy i a j spojené vazbou modifikujeme výše uvedené vzorce tak, aby odstředivésíly působící na oba atomy měly stejnou velikost a opačná znaménka, tedy opravný vektorλ~r(t) se rozdělí v inverzním poměru hmotností:
~ri(t+ h) = ~rVerlet,i(t+ h) + λ1/mi
1/mi + 1/mj
~rij , (62)
~rj(t+ h) = ~rVerlet,j(t+ h)− λ 1/mj
1/mi + 1/mj
~rij , (63)
kde
λ =|~rVerlet,ij(t+ h)|2 − |~rij(t)|2
2~rVerlet,ij(t+ h) · ~rij(t)(64)
a ~rVerlet,i je dáno vzorcem (55).Výše uvedený postup zachovává těžiště i hybnost soustavy atomů, což je podstatné pro
bezproblémovou aplikaci na složité molekuly. Pro ně totiž postupujeme iteračně. Procházímevšemi vazbami v cyklu a aplikujeme výše uvedenou korekci. Jelikož vazby jsou obecně spo-jené, korigováním jedné vazby můžeme způsobit chybu vazbě jiné. Postup proto opakujemetak dlouho, až bude chyba všech vazeb menší než nějaká předem daná přesnost; proto takénevadilo, že jsme v (61) zanedbali členy vyššího řádu (stejně můžeme nahradit jmenovatel(61) třeba |~r(t)|2 = l2).
Podobně jako v lineárních iteračních metodách můžeme konvergenci urychlit zavedením relaxace: λ vypočtenépomocí (61) znásobíme jistým číslem q. Ukazuje se, že vhodná hodnota je asi q = 1.3, což souvisí s typickouvelikostí vazebných úhlů.
Algoritmus SHAKE můžeme použít i pro MD rigidních molekul; pro osově symetrickémolekuly stačí jedna vazba, pro nesymetrické jsou nutny tři do trojúhelníka. Je nevhodné(ale při dostatečné pečlivosti ne nemožné) mít úlohu přeurčenu zavedením více vazeb, prosimulaci methanu je tedy vhodné zvolit např. čtyři vodíky jako referenční a polohu uhlíkuz nich počítat a rovněž přepočítat síly podle zákonů platných pro tuhá tělesa.
Algoritmus SHAKE je založen na Verletově metodě a má proto jednu podstatnou výhodu:je časově reverzibilní s přesností, kterou dosáhneme při iteracích pro délky vazeb.
7 Monte Carlo
Smyslem počítačových simulací je generovat konfigurace systému mnoha částic a tyto kon-figurace pak použít ke stanovení různých termodynamických či strukturních veličin, což ob-vykle představuje výpočet nějaké střední hodnoty. Metoda MC generuje konfigurace právěs ohledem na efektivní výpočet středních hodnot. Název metody pak pochází z toho, že narozdíl od deterministické MD používá generátor náhodných čísel, tj. počítačový kód, kterýprodukuje náhodná čísla s danými statistickými vlastnostmi.
39
Obr. 14: K výpočtu pravděpodobnosti protnutí linky náhodně hozenou jehlou. Poloha jehly jepopsaná vzdáleností od linky z a úhlem θ, k protnutí dojde v oblasti parametrů vyznačené naobrázku vpravo
7.1 Monte Carlo integrace
Často potřebujeme spočítat mnohodimenzionální integrál. Můžeme využít metodu „náhod-ného stříleníÿ, též známou jako naivní Monte Carlo.
Nepočítačovým příkladem Monte Carlo integrace je Buffonova jehla. Mějme linkovaný papír s linkami vzdá-lenými d. Pravděpodobnost, že náhodně hozená jehla délky l, l ≤ d, protne linku, je 2l/πd.
Tento vzorec odvodíme následovně. Podle obr. 14 je pravděpodobnost protnutí rovna poměru obsahu modréoblasti k celkové ploše oblasti (d/2)× (π/2),
π =1d/2
∫ d/2
0dz
1π/2
∫ π/2
0
(z <
l
2cos θ
)dθ =
1d/2
1π/2
∫ π/2
0
l
2cos θdθ =
2lπd
kde výraz (a < b) dává 1, pokud nerovnost platí, jinak 0. Nechť házíme ncelkem krát a z toho jehla protne linkunprotne krát. Číslo π vypočteme podle vzorce
π ≈ 2lπd
, kde π =nprotnencelkem
Dále nás zajímá, s jakou přesností jsme π stanovili. Dá se ukázat, že standardní chyba odhadu π je
δπ ≈√
π(1− π)n− 1
Obecně při výpočtu integrálu funkce f(x1, . . . , xD) přes oblast Ω vD-rozměrném prostoruplatí: ∫
Ωf(x1, . . . , xD)dx1 . . . dxD ≈
|Ω|n
n∑k=1
f(x(k)1 , . . . , x
(k)D ) (65)
kde (x(k)1 , . . . , x
(k)D ) značí k-tý náhodný bod v oblasti Ω, jejíž D-objem je |Ω|. Metoda náhod-
ného střílení není příliš přesná, je však jednoduchá a v mnoha případech složitých mnoho-násobných integrálů jediná možná.
Zkusme spočítat integrál∫ 10 sin(1/x) dx. Integrand vypadá velmi ošklivě, takže běžné metody (jako Simpso-
novo pravidlo) by bez dalších úprav byly nepřesné. Výpočet metodou MC je velmi jednoduchý, v snadno pocho-pitelném „pseudokóduÿ např.:
40
n := 1000000000sum := 0FOR i:=1 TO n
sum := sum + sin(1/u0,1)ENDPRINT sum/n 0 1
x
-1
0
1
sin
(1/x
)
kde u0,1 je náhodné číslo z intervalu (0, 1) (v různých programovacích jazycích dostupné např. jako rand() neborandom(1) apod.). Použili jsme vzorec (65) s Ω = (0, 1), takže |Ω| = 1.
Výše uvedený kód implementovaný v C mi dal 0.504086; snadno lze dopočítat i standardní chybu tohotovýsledku, 0.000020. Výpočet trval 46 s na 1 procesoru Intel Core DUO E8600, 3.33 GHz. Přesná hodnota je13
sin(1)− Ci(1) = 0.50406706.
7.2 Importance sampling
V případě výpočtu střední hodnoty systému mnoha částic, rov. (30), však tato metoda zcelaselhává z jednoho prostého důvodu: náhodný výběr nedělá rozdíl mezi konfiguracemi, kterémají velkou pravděpodobnost výskytu a tudíž podstatně přispívají k hodnotě 〈X〉, a málopravděpodobnými či přímo nemožnými. Nejjednodušeji se tato skutečnost demonstruje napřípadu systému N tuhých koulí. Generovat náhodně konfigurace tohoto systému znamenánáhodně volit polohy N koulí v objemu V (náhodně vybírat konfiguraci ~rN,(k) v (30)) avypočítat vždy součin funkce X s Boltzmannovým faktorem exp(−βU). Pro libovolnou kon-figuraci však s velikou pravděpodobností (pro vysoké hustoty téměř s jistotou) dojde k tomu,že alespoň dvě koule se budou protínat, a tedy Boltzmannův faktor bude nula. Pravděpodob-nost získání možné konfigurace je téměř nula a výraz (30) nebude vůbec definován (dělenínuly nulou).
Řešení, které se nabízí k odstranění tohoto problému, je toto: při výpočtu střední hod-noty nebudeme uvažovat libovolné konfigurace, ale přednostně ty, které podstatně přispívajík hodnotě integrálu (angl. importance sampling). Otázkou pak je
1. jak toto realizovat a
2. dostaneme-li správný odhad 〈X〉.
Pokud se týká problému 1., je samozřejmě obtížné vytvořit možnou (dostatečně pravdě-podobnou) konfiguraci pouhým vkládáním částic do prázdného prostoru (viz výše uvedenýpříklad). Máme-li však nějakou (dostatečně pravděpodobnou) konfiguraci, pak by již nemělbýt takový problém z této konfigurace vytvořit jinou dostatečně pravděpodobnou konfiguraci.Intuitivně lze tedy navrhnout následující schéma, které pro jednoduchost vysvětlíme opět nasystému tuhých koulí. V tomto případě je totiž Boltzmannův faktor buď nula, nebo jedna, atedy každá konfigurace, ve které nedochází k překryvu žádných koulí, má stejnou pravděpo-dobnost výskytu, zatímco konfigurace s překryvem se nikdy nevyskytnou. Předpokládejmetedy, že se nám podařilo nějakým způsobem vytvořit možnou konfiguraci. Změníme-li nynípolohu jedné (nebo i více) koulí tak, že opět nedojde k překryvu, dostaneme další možnoukonfiguraci, a tak můžeme pokračovat a vytvořit posloupnost konfigurací, která bude našívybranou posloupností v kanonickém souboru.
13v systému Maple: integrate(sin(1/x),x=0..1); evalf(%);v systému Mathematica: Integrate[Sin[1/x],x,0,1]//N
41
Jak uvidíme dále, lze výše uvedený princip rozšířit a najít vhodnou vybranou posloupnostkonfigurací k výpočtu střední hodnoty i ve složitějším případě, kdy Boltzmannův faktornabývá více hodnot než jen dvou. Rigorózně se tento problém řeší pomocí Markovovýchřetězců, což nakonec zajistí i správnost výsledku (bod 2.). Zkusme to však nejprve zcelaintuitivně a napišme algoritmus zvaný Metropolisova metoda
1. Vybereme částici, i (např. náhodně)
2. Zkusíme s ní náhodně hýbnout, např.:
xzki = xi + u(−d,d) ,
yzki = yi + u(−d,d) ,
zzki = zi + u(−d,d)
u(−d,d) je náhodné číslo rovnoměrně rozdělené v intervalu (−d, d). To znamená, žepravděpodobnost opačného pohybu je stejná
3. Spočteme změnu potenciální energie, ∆U = U zk − U
4. – Je-li ∆U ≤ 0, změnu přijmeme a pokračujeme s novou konfigurací– Je-li ∆U > 0, změnu přijmeme s pravděpodobností exp(−β∆U); to znamená, žes pravděpodobností 1− exp(−β∆U) pokračujeme se starou konfigurací
Tento MC krok opakujeme mnohokrát.
„Přijmout s pravděpodobností πÿ se realizuje následujícím algoritmem:IF u(0,1) < π
THEN přijmoutELSE odmítnout
Pokud vám tento algoritmus dělá potíže, uvědomte si, že funguje u obou limitních případů: pro π = 0 nenípodmínka u(0,1) < π nikdy splněna (všechna náhodná čísla z intervalu (0, 1) jsou větší než π = 0; pro π = 1 jepomínka splněna vždy (všechna náhodná čísla z intervalu (0, 1) jsou menší než π = 1).
Proč dostaneme správný kanonický soubor? Nechť při pohybu z konfigurace A na B dojdeke zvýšení energie, U(B) − U(A) = ∆U > 0. Pak daný zkušební pohyb (novou konfiguraciB) přijmeme s pravděpodobností exp(−β∆U) < 1. Se stejnou pravděpodobností, se kteroujsme se pokusili A změnit na B, se někdy pokusíme změnit B na A (protože v kroku 2je pravděpodobností opačného pohybu stejná). Takovou změnu ale vždy přijmeme, protoženyní U(A) − U(B) = ∆U < 0. Změnu A → B provádíme s pravděpodobností exp(−β∆U),změnu B→ A vždy, tedy A→ B provádíme exp(−β∆U) krát častěji než B→ A. Proto propoměr pravděpodobností platí
π(B)π(A)
= exp(−β∆U) =exp[−βU(B)]exp[−βU(A)]
bez ohledu na to, zda je ∆U kladné či záporné (A a B vystupují v úvahách zcela symetricky).Stejná úvaha platí pro všechny uvažované páry změn konfigurací, všechny poměry jsou dánypoměrem Boltzmannových pravděpodobností, tedy pravděpodobnost nalezení jakékoliv do-sažitelné konfigurace musí být úměrná Boltzmannovu faktoru.
42
7.3 Trocha teorie: náhodné veličiny
Zkusme nyní to samé formulovat trochu vědečtěji.Symbolem S označme náhodnou veličinu. Ta nabývá hodnoty z jisté (pro jednoduchost
konečné) množiny stavů (pro nás: konfigurací) Ai, i = 1, . . . ,M , s jistou pravděpodobnostiπ(Ai) = πi, přičemž platí
∑i πi = 1. Náhodnou veličinou může být např. hod hrací kostkou.
Pak M = 6 a πi = 1/6 (je-li kostka zcela pravidelná a jestliže nefixlujeme).Markovův řetězec je posloupnost náhodných veličin S(k), k = 1, . . . ,∞ taková, že
událost, kterou pozorujeme v „časeÿ k + 1, závisí na tom, jaká událost byla pozorovánav čase k. Konkrétně, jestliže se v čase k vyskytne událost Ai s pravděpodobností π(k)
i (tj.náhodná veličina S(k) nabývá hodnoty Ai s pravděpodobností π
(k)i ), pak v čase k + 1 se
událost Aj vyskytne s pravděpodobností, pro kterou platí
π(k+1)j =
M∑i=1
π(k)i Wi→j (66)
nebo zapsáno vektorově pro π = (π1, . . . ,πM) (řádkový vektor)
π(k+1) = π(k) ·W (67)
Veličina W se nazývá matice přechodu. Prvky Wi→j ≡ W (Ai → Aj) mají fyzikálnívýznam pravděpodobnosti přechodu ze stavu Ai do stavu Aj (jsou tedy nezáporné). Maticepřechodu musí splňovat normovací podmínku
M∑j=1
Wi→j = 1 pro všechna i (68)
Tato podmínka znamená, že z konfigurace Ai vznikne (s pravděpodobností jedna) jednaz konfigurací Aj, j = 1, . . .M .
Markovovým řetězcem může být třeba posloupnost hodů (jednou) kostkou. Můžeme zdepřipustit, že při jednom hodu je číslo, které padne, poněkud ovlivněno předchozí polohou(nedostatečně kostkou v dlani zatřepeme), není však (přímo) ovlivněnou polohou před dvěmahody.
Máme-li Markovův řetězec, můžeme si položit řadu otázek. Například, vyjdu-li ze stavuAi, jaká je pravděpodobnost toho, že po k krocích dojdu do stavu Aj? Jak bude závisetvýskyt stavu Aj na k? Tyto otázky osvětlíme nejlépe na příkladu.
Mám v kanceláři počítač. Ten mi vždy funguje, ale horší je to se sítí, ta někdy funguje, někdy ne. Dlouhodobýmpozorováním jsem zjistil, že
1. funguje-li síť dnes, je 90% pravděpodobnost, že bude fungovat i zítra(tj. spadne s pravděpodobností 10 %);
2. nefunguje-li síť dnes, pak s pravděpodobností 70 % nebude fungovat ani zítra(atjťáci ji spraví s pravděpodobností 30 %) .
Stav sítě tedy nabývá dvou hodnot, A1 = funguje a A2 = nefunguje. Pravděpodobnosti v čase k lze popsatdvourozměrným vektorem
π(k) = (π(k)1 ,π(k)2 ).
Matice přechodu je v tomto případě
W =(
0.9 0.10.3 0.7
)
43
Jestliže tedy včera byl stav sítě popsán vektorem π(1), pak pro dnešní stav platí:
π(2) = π(1) ·W
Tedy např. je-li π(1) = (1, 0), pak π(2) = (0.9, 0.1), a je-li π(1) = (0, 1), pak π(2) = (0.3, 0.7). Takto mohupokračovat a druhý den mám rozložení
π(3) = π(2) ·W = π(1) ·W2
a tedy
π(3) =
(0.84, 0.16), je-li π(1) = (1, 0),
(0.48, 0.52), je-li π(1) = (0, 1)
Pokračuji-li dále, zjistím, že rozložení π(k) pro velká k nebude už vůbec záviset na π(1) a dostanu tzv. limitnírozložení
limn→∞
π(n) = π = (0.75.0.25) (69)
tedy průměrná pravděpodobnost fungování je 75%. Opravdu je tomu tak ale vždycky? Zkusme něco jiného. Byljsem na dovolené a nevím, zda včera síť fungovala nebo ne. Jeden kolega říká, že síť nefungovala, druhý tvrdí, žeano. Pravděpodobnosti pro a proti jsou tedy asi stejné a počáteční stav proto vezmeme ve tvaru
π(1) = (0.5, 0.5)
Jaká je tedy pravděpodobnost, že síť bude fungovat, až dnes přijdu do práce?
π(2) = π(1) ·W = (0.6, 0.4)
a tedy pravděpodobnost je 60%. Budu-li nyní pokračovat ve výpočtu pravděpodobnosti dál, zjistím, že po několikadnech dostanu opět rozložení (69).
Ve statistické fyzice nás zajímá měření veličin. Zde může jako příklad veličiny (pozorovatelné) sloužit výdělek:jestliže síť funguje, vydělám X(funguje) = 2000 Kč za den, jestliže nefunguje, nemohu pracovat a beru pouzeX(nefunguje) = 500 Kč. Průměrný výdělek je dán střední hodnotou
Matice přechodu v našem jednoduchém příkladu má tu vlastnost, že systém ztrácí „pa-měťÿ, tj. po jisté době pravděpodobnost πj nalezení jevu Aj vůbec nezávisí na tom, z jakéhorozložení jsme vyšli. A jak tento příklad souvisí s původním problémem vybrané posloup-nosti, s níž budeme procházet konfigurační prostor? Jednoduše. Mějme posloupnost stavůčili nyní konfigurací A(k)nk=1 vybranou z Markovova řetězce S(1),S(2), . . . s limitním roz-ložením
πj =exp(−βUj)∑Mk=1 exp(−βUk)
≡ exp(−βUj)Q
(70)
kde jsme označili Uj = U(Aj). Pak střední hodnota veličiny X podél tohoto řetězce,
X =1n
n∑k=1
X(k) (71)
kde X(k) = X(A(k)), bude pro rostoucí n konvergovat k souborové střední hodnotě danévztahem
〈X〉 =M∑j=1
πjX(Aj) ≡M∑j=1
πjXj . (72)
Zbývá určit podmínky, za kterých π(k) = π(1) ·Wk konverguje k limitnímu rozložení π.Jestliže:
44
1. všechny stavy jsou dosažitelné z libovolného stavu v konečném čase s nenulovou prav-děpodobností a
2. žádný stav není periodický (stav Ai je periodický, jestliže existuje perioda m taková,že je-li π(k)
i = 0, pak π(k+m)i = 0 a je-li π(k)
i 6= 0, pak π(k+m)i 6= 0),
pak se množina stavů nazývá ergodická a pro libovolné počáteční rozložení pravděpodob-ností π(1) existuje limita π = limk→∞ π(k). Rozložení pravděpodobnosti π je tedy řešenímrovnice
π ·W = π (73)
a toto řešení je jediné.
Jinými slovy, vektor stavů π je vlastním levým vektorem stochastické maticeW. Lze ukázat, že všechna dalšívlastní čísla jsou v absolutní hodnotě menší než 1.
7.3.1 Určení matice přechodu
Potřebujeme tedy zkonstruovat posloupnost (Markovův řetězec) konfigurací tak, aby se prav-děpodobnost výskytu jednotlivých konfigurací rovnala Boltzmannově váze (70), která takbude představovat limitní rozložení jisté, zatím neznámé matice přechodu. Toto je právěopačný problém než ten, který se obvykle řeší v teorii Markovových řetězců, tj. k danématici přechodu nalézt limitní rozdělení.
Pro určení matice přechodu máme celkem tři podmínky:
Wi→j ≥ 0 pro všechna i, j = 1, . . . ,M (74)M∑j=1
Wi→j = 1 pro všechna i = 1, . . . ,M (75)
π ·W = π . (76)
Poslední rovnice vyjadřuje podmínku tzv. detailní rovnováhy. Toto je nutná podmínka za-jišťující vlastnosti limitního (stacionárního) rozložení pravděpodobnosti. Umíme ji splnit,požadujeme-li silnější podmínku, tzv. podmínku mikroskopické reverzibility. Stačí, aby
πiWi→j = πjWj→i (77)
a podmínka detailní rovnováhy je splněna.
Platnost podmínky (76) za předpokladu platnosti (77) snadno ukážeme, aplikujeme-li operátor∑Mi=1 na obě
strany rovnice (77) a dosadíme-li jedničku za∑Mi=1Wj→i (viz rov. (75)) na pravé straně.
Rovnice (74)–(76) neurčují matici přechodu jednoznačně. Soustava (75) a (76) totiž před-stavuje celkem 2M rovnic pro M ×M neznámých.
Využijme nyní toho, že πi = exp(−βUi)/Q je limitní rozložení. Po dosazení do (77)dostaneme buď
Wj→i
Wi→j=
πi
πj
= exp[−β(Ui − Uj)] (78)
nebo Wj→i = Wi→j = 0. Jak je vidět, podařilo se nám zbavit, alespoň v poměru pravděpo-dobností, neznámé funkce Q; v obecném případě (jsou-li πi jiné než Boltzmannovy pravdě-podobnosti) tento výsledek znamená, že nám stačí zadat jen relativní pravděpodobnosti anemusíme se starat o součet čísel πi.
45
Metropolisova metoda. Snadno se nyní přesvědčíme, že matice přechodu definovanávztahy
Wi→j =
αi→j pro i 6= j a πj ≥ πi,
αi→jπj
πi
pro i 6= j a πj < πi,
1−∑k, k 6=i
Wi→k pro i = j
(79)
kde αi→j je libovolná symetrická stochastická (tj. platí pro ni (74) a (75)) matice, splňujepodmínku mikroskopické reverzibility. Toto je matice navržená Metropolisem a používanádodnes (a odpovídající algoritmu na str. 42). Poznamenejme ještě, že první dva řádky v (79)lze zapsat kompaktně jako
Wi→j = αi→j min
1,
πj
πi
pro i 6= j (80)
Algoritmus jsme již uvedli. Zopakujme si ho nyní včetně počítačového „pseudokóduÿ:
1. Zvolíme částici, kterou se bude hýbat, mřížkový bod, . . .
4. Konfiguraci přijmeme (A(k+1) = Azk) s pravděpodobností ppřij,v opačném případě odmítneme:
Varianta 1 Varianta 2 Varianta 3u := u(0,1) u := u(0,1) IF ∆U < 0IF u < min1, e−β∆U IF u < e−β∆U THEN A(k+1) := Azk
THEN A(k+1) := Azk THEN A(k+1) := Azk ELSEELSE A(k+1) := A(k) ELSE A(k+1) := A(k) u := u(0,1)
IF u < e−β∆U
THEN A(k+1) := Azk
ELSE A(k+1) := A(k)
5. k := k + 1 a opakujeme od začátku
Metoda tepelné lázně. Jinou možností, používanou především pro mřížkové modelys krátkodosahovým potenciálem, je tzv. metoda tepelné lázně (heat-bath method). Vybe-reme si mřížkový bod (příp. i skupinu) a pro jeho dané okolí vypočteme Boltzmannovyfaktory všech jeho možných stavů. Sečteme, faktory vydělíme součtem (částečnou statistic-kou sumou), čímž dostaneme Boltzmannovy pravděpodobnosti. Nový spin vybereme z tétomnožiny s Boltzmannovou pravděpodobností. Metodu lze interpretovat tak, že daný spin„ponoříme do lázněÿ s danou teplotou. Metoda tepelné lázně je vhodnou alternativou, jestližeMetropolisova metoda se stane neefektivní pro příliš mnoho odmítnutých konfigurací. Vý-hodná je tehdy, pokud lze všechny pravděpodobnosti pro všechna okolí předem tabelovat.
Více matematicky: Předpokládejme, že umíme spočítat částečnou statistickou sumu pro jistou podmnožinuCpart celého konfiguračního prostoru C, který předpokládáme diskrétní. Matice přechodu metody tepelné lázně je
Wi→j = exp(−βUj) /∑
Ak∈Cpart
exp(−βUk) pro Ai, Aj ∈ Cpart (81)
46
Abychom byli konkrétní, zvolme si za Cpart množinu všech konfigurací lišících se jen hodnotou proměnné na jednommřížkovém bodu b. Tuto proměnnou (spin) označíme sb. Pak lze napsat
W (sb → s′b) = exp[−βU(s′b)] /∑s∗b
exp[−βU(s∗b)] (82)
kde závislost veličin na ostatních spinech, které se v daném kroku nemění, již nevyznačujeme. Konfigurace spinu sbvybíráme tedy s Boltzmannovou váhou danou energiemi konfigurací a teplotou, což si lze představit jako ponořeníspinu do termostatu či „tepelné lázněÿ o dané teplotě T . Spin tak „zapomeneÿ na předchozí stav a Wi→j nezávisína Ai.
7.3.2 Zkušební změna konfigurace
Zbývá odpovědět na otázku, co je matice αi→j vystupující v maticích přechodu (79). Tatomatice udává podmíněnou pravděpodobnost generování zkušební konfigurace Aj z konfigu-race Ai. Pro klasický spojitý systém budeme mít místo matice funkci α(~rN → ~r′N).
Ve standardní Metropolisově metodě je tato matice symetrická. Máme-li tedy konfigu-raci Ai, pak konfiguraci Aj musíme vygenerovat s pravděpodobností stejnou, jakou bychomz konfigurace Aj generovali konfiguraci Ai. Podobně je tomu ve spojitém případě, kde prav-děpodobnost je nahrazena hustotou pravděpodobnosti; nesmíme jen zapomenout na to, žetato hustota pravděpodobnosti je definovaná vzhledem ke kartézským proměnným ~rN a můžese změnit, jestliže bychom chtěli přejít k jiným souřadnicím (např. sférickým). Dále musímemít při návrhu zkušebního posunutí na paměti ergodičnost, tedy aby systém mohl (pomocíco nejmenšího počtu kroků) přecházet z jedné části konfiguračního prostoru do jiné. Spe-ciálně ve spojitých modelech je nutno dát pozor na to, aby systém snadno překonal velképotenciálové bariéry, např. u vnitřních stupňů volnosti. Podobně u některých mřížkovýchmodelů může nastat situace, kdy není možné jednu konfiguraci přeměnit na druhou pouzepostupnými záměnami jednotlivých spinů; pak je nutné měnit v jednom kroku celou skupinuspinů najednou.
Uveďme několik příkladů.
Diskrétní mřížový systém. U simulace Isingova feromagnetu máme dva stavy, + a−. Zkušební změna konfigurace je jednoduchá: vybereme (náhodně) spin a změníme jehohodnotu na opačnou. Pak aplikujeme Metropolisovo kritérium, zda tuto změnu přijememe.14
Atomární systém. Zkušební translace může být dána krokem 2 na str. 42. Asi efektivnějšíje mít ve trojrozměrném prostoru toto zkušební posunutí izotropní, např. uvnitř koule opoloměru d. Podstatné je, aby k opačnému pohybu došlo se stejnou pravděpodobností.
Směsi Mikroreverzibilitu musíme zachovat ve všech krocích návrhu algoritmu. Kdybychomnapříklad simulovali ternární směs složenou z molekul A, B a C a střídali pohyby v pořadí
[: (náhodná molekula typu A) (náhodná molekula typu B) (náhodná molekula typu C) :]
([: :] značí repetici), nebylo by to správné, protože opačné pořadí C, B, A není stejné.
14Pro tento jednodudchý systém existují mnohem efektivnější metody, v nichž se překlápí celé oblastimřížky najednou
47
0.005
0.010
0.015
0.020
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
δP
d
♦
♦
♦♦
♦♦ ♦
♦
♦
♦
♦
0.005
0.010
0.015
0.020
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
χ
♦
♦
♦♦
♦♦♦
♦
♦
♦
♦
Obr. 15: Závislost standardní chyby měření viriálového tlaku δP na maximální délce zkušebníhoposunutí d (vlevo) a na zlomku přijatých konfigurací χ (vpravo). Každý bod byl získán z 107 MCkroků (krok = pokus pohnout jednou částicí) na systému 64 LJ atomů (σ = 1, ε = 1) při T = 1.2,ρ = 0.8.
7.3.3 Zlomek přijetí a nastavení parametrů
Důležitou charakteristikou MC simulací je tzv. zlomek přijetí (acceptance ratio), čili poměr
χ =počet přijatých konfigurací
počet všech generovaných konfigurací(83)
Uvažujeme-li simulace klasických spojitých systémů, bude tento poměr zřejmě záviset navelikosti zkušebního posunutí d v rovnici xzk
i = xi + u(−d,d) či v jiné podobné. Bude-li např.d příliš malé, bude změna energie v jednom kroku malá, Boltzmannův faktor exp(−β∆U)blízký jedničce, téměř každá konfigurace bude přijata a χ bude blízké jedničce, ale efektivitasimulace bude nízká, protože konfigurace se od sebe málo liší. Bude-li d velké, bude změnaenergie v jednom kroku velká a většinou kladná, Boltzmannův faktor exp(−β∆U) blízký nulea téměř každá konfigurace bude odmítnuta, což je zřejmě také neefektivní. Existuje protojistá optimální hodnota zlomku přijatých konfigurací, která mnohdy leží okolo χ = 1/3.
Ve speciálních případech může být optimální χ jiné. Příkladem jsou řídké systémy, zvláštěblízko kritického bodu: zde je efektivnější mít délku zkušebního posunutí srovnatelnou s veli-kostí simulovaného systému i za cenu velmi malého zlomku přijatých konfigurací (třeba 1 %)než malou změnu konfigurace v jednom kroku.
7.4 Rosenbluthovo vzorkování
Simulace polymerů (MC i MD) jsou zvláště obtížné, protože pohyby dlouhých a často pro-pletených řetězců jsou navzájem bržděny a dlouho trvá, než se změní konformace. To jefyzikální vlastnost, kterou lze těžko vylepšit v molekulové dynamice, jež simuluje reálnou
48
p =1
4× 3× 3× 3p =
14× 3× 3× 2
Obr. 16: Budujeme lineární řetězec na mřížce jako náhodnou procházku bez protínání tak, žev okamžiku přidání článku (červené tyčinky) volíme náhodně pokračování. Řetězec vpravo mávětší pravděpodobnost, protože v posledním kroku jsme volili jen ze dvou pokračování, zatímcovlevo ze tří
dynamiku s fyzikálně dlouhými relaxačními časy. Možnosti metod MC jak efektivně vzorko-vat konfigurační prostor jsou větší. Zjevně nebudou příliš efektivní lokální pohyby jednohočlánku řetězce. Je však možné rotovat najednou několika články; za příklad nechť sloužítzv „crankshaft moveÿ: \/\ → /\/. Ještě efektivnější je (u lineárních řetězců) pohyb zvaný„reptationÿ (housenkový pohyb): článek z jednoho konce přemístíme na druhý konec.
Jinou metodou je zapomenout na to, co jsme se dosud o metodách MC naučili, a snažitse vytvořit konfiguraci polymeru od začátku. I když to může být obtížné, výhodou je, žekonfigurace jsou zcela nekorelované.
Abychom si vysvětlili princip, uvažujme nejprve polymer v dobrém rozpouštědle15. Tentopolymer lze popsat modelem náhodné procházky bez protínání. Navíc pro jednoduchostpoložme polymer na mřížku, viz obr. 16. Budujme řetězec od malého čtverečku. Na začátkumáme 4 možnosti, jak vést první vazbu. Pro druhou vazbu máme ale jen tři možnosti, protoženemůže zpátky. To samé platí pro třetí vazbu. Ale čtvrté vazby se liší u konfigurace vlevoa vpravo: zatímco vlevo jsme si vybrali jednu ze tří možností, vpravo máme jen dvě (vazbadoleva by vedla k překryvu s počátečním článkem). Konfiguraci vpravo proto dostanemečastěji, než bychom měli.
Řešením je znásobit váhu takto získané konfigurace číslem, které se rovná součinu všechtěchto možností (tj. jmenovatelem zlomků). Tomuto součinu R se říká Rosenbluthova váhakonfigurace:
R =N∏i=1
Ri
kde součin je přes všechny články postupně budovaného řetězce. U náhodné procházky jsouRi celá čísla – počty pokračování. Zobecnění na libovolnou interakční energii je přímočaré,místo nuly (nelze pokračovat, protože by se atomy překrývaly) a jedničky (lze pokračovat)máme Boltzmannův faktor změny potenciální energie po přidání jednoho článku,
Ri =k∑l=1
exp[−βUi(l)]
15V dobrém rozpouštědle se řetězec polymeru rozvine, protože jeho články se přitahují přednostně k mo-lekulám rozpouštědla, resp. se navzájem odpuzují. Naopak ve špatném rozpouštědle se svine, protože jevýhodnější, aby sousedily dva články než článek a rozpouštědlo. Theta-rozpouštědlo (θ-rozpouštědlo) je„něco meziÿ – článek interaguje stejně s molekulou rozpouštědla jako s jiným článkem.
49
kde součet je přes všechna pokračování řetězce v daném kroku i a Ui(l) značí energii danéhopokračování (u 2D systému z obr. 16 bylo k = 4). Další článek z k vybíráme s pravděpodob-ností úměrnou Boltzmannově faktoru exp[−βUi(l)]/Ri.
Princip lze zobecnit i na reálný řetězec (tj. není na mřížce) tak, že místo k pokračovánína mřížce (ve směrech ±x, ±y . . . ) generujeme určitý počet k náhodných směrů. Je uži-tečné si uvědomit, že pro k = 1 generujeme zcela náhodný řetězec, jehož váha je pak dánaBoltzmannovým faktorem, protože
R =N∏i=1
Ri =N∏i=1
exp[−βUi(1)] = exp
[−β
N∑i=1
Ui(1)
]kde Ui(1) je změna energie po přidání článku a součet není nic jiného než celková potenciálníenergie. Dostáváme tedy „naivníÿ Monte Carlo integraci (která asi nebude příliš efektivní).Naopak pro velmi velké k vybíráme další článek ve shodě s Boltzmannovou pravděpodobností(metoda tepelné lázně) a generovaný řetězec vzorkuje kanonické rozdělení pravděpodobnostímnohem lépe. I přes to je metoda tím méně efektivní, čím delší řetězec tvoříme; typicky mávětšina řetězců malé váhy (nejsou signifikantní) a málo řetězců váhy velké (ty nás zajímají).
8 Měření veličin v simulacích
Veličiny stanovované v simulacích můžeme zhruba rozdělit na
1. Makroskopické veličiny:
(a) Mechanické veličiny (energie, tlak atd.),
(b) Entropické veličiny (chemický potenciál).
2. Strukturní (mikroskopické) veličiny.
3. Další pomocné veličiny vypovídající o vývoji a stavu systému (parametry uspořádání,integrály pohybu v MD).
4. Transportní vlastnosti (kinetické veličiny, např. difuzivita a viskozita).
Simulace produkuje posloupnost konfigurací čili trajektorii. Převážná většina veličin X„měřenýchÿ v průběhu simulace je dána aritmetickým průměrem (označíme ho X) hod-not stanovených z jednotlivých konfigurací, Xi
16. Tímto průměrem aproximujeme skutečnoustřední hodnotu, 〈X〉.
〈X〉 ≈ X =1M
M∑i=1
Xi
K tomu, aby X bylo dobrým odhadem 〈X〉, musí být splněny dvě podmínky:
• Systém musí být v rovnováze, tj. na grafu Xi nesmí být patrný žádný trend. Před zapo-četím měření (produktivního běhu) musíme tedy dostatečně dlouho simulovat (míchat,zrovnovážňovat) a sledovat vývoj systému, viz obr. 17.
• Měříme dostatečně dlouho.16Obvykle nepoužíváme všechny konfigurace, ale provedeme několik MC či MD kroků, aby se konfigurace
dostatečně změnila, než provedeme měření. Např. v MD může být časový krok 1 fs, ale vzorkujeme po 10 fs,v MC mezi dvěma měřeními hýbneme každou částicí alespoň jednou.
50
Obr. 17: Před započetím „měřeníÿ středních hodnot musíme systém zrovnovážnit. Závislost tlakuna čase simulace pro model vody startovaný z pravidelného rozmístění molekul vody
8.1 Odhad chyby aritmetického průměru
Kromě střední hodnoty nás ovšem zajímá, jak přesný je výsledek. Přesnost vyjadřujemepomocí standardní chyby δX, která je definována vztahem
δX2 = 〈(X − 〈X〉)2〉
Tomuto vzorečku často studenti nerozumějí. V simulaci či reálném měření totiž máme jen jednu sérii měřeníXi, i = 1, . . . ,M , a jeden aritmetický průměr X, který ovšem nevyjde přesně jako správná střední hodnota 〈X〉.Kdybychom však provedli tu samou simulaci mnohokrát, dostali bychom mnoho výsledků X, které by se ovšemponěkud lišily; jejich aritmetický průměr by se blížil k 〈X〉 mnohem lépe. Protože však X je náhodná veličina,má určité rozdělení, totiž mnohdy Gaussovo; zjevně také platí 〈X〉 = 〈X〉. Pokud stanovíme δX a rozdělení jeGaussovo, víme, že s pravděpodobností 68 % bude 〈X〉 v intervalu hodnot X ± δX a s pravděpodobností 95 %v intervalu X ± 2δX. Potíž je v tom, že odhad δX musíme založit na jedné posloupnosti hodnot Xi, protože nicjiného nemáme.
Pokud by hodnoty Xi byly nezávislé, nebyl by s odhadem standardní chyby δX problém;následující vzoreček jistě znáte
δX2 ≈∑M
i=1(Xi −X)2
M(M − 1)(84)
Tak je tomu bohužel jen málokdy, protože data jsou korelovaná, viz obr. 18. To je proto, žei jednotlivé konfigurace jsou korelované. Vzorec (84) lze upravit na
δX2 ≈∑M
i=1(Xi −X)2
M(M − 1)× (1 + 2τ)
kde τ je korelační čas (též délka). Lze pro něj odvodit vztah
τ =∞∑k=1
ck, ck =〈∆Xi∆Xi+k〉〈∆X2〉
kde ∆Xi = Xi − 〈X〉, což při výpočtu aproximujeme ∆Xi ≈ Xi − X. Číslo ck se nazýváautokorelační koeficient. Lidově vyjádřeno, korelační čas vyjadřuje čas (počet kroků), zakterý se hodnoty veličiny X „podstatně změníÿ.
51
nekorelovaná data
korelovaná data (korelační koeficient = 0.9)
korelovaná data (korelační koeficient = 0.99)
korelovaná data (různé typy korelací)
Obr. 18: Různě korelovaná data. Korelační čas je vyznačen vodorovou červenou úsečkou
Obr. 19: Odhad chyby průměru korelovaných dat. Vlevo: zprůměrujeme dostatečně dlouhé bloky,které již lze považovat za nekorelované. Vpravo: nakreslíme kumulativní klouzavý průměr (odzačátku simulace), standardní chyba je zhruba 0.6 krát rozdíl maximum − minimum z druhépoloviny tohoto klouzavého průměru. (Obrázky jsou ze stejných dat, ale osy y mají různá měřítka)
52
Chceme-li tedy odhadnout chyby aritmetického průměru z korelovaných dat, musíme buďznát odhad τ nebo nějak korelace odstranit. Asi nejčastěji používaná je bloková metoda:posloupnost dat se rozdělí na bloky, které jsou delší než τ , takže je můžeme považovat zanekorelované. Vypočteme průměry přes bloky a z těchto blokových průměrů pak odhadnemechyby podle vzorce (84). Druhou dosti přibližnou leč v praxi vyhovující metodou je metodazaložená na kumulativním klouzavém průměru od začátku simulace (cumulative running(moving) average). Vynášíme graf veličiny
Xn =1n
n∑i=1
Xi
v závislosti na n. Tomuto grafu se někdy říká konvergenční profil, protože X = Xn alimn→∞Xn = 〈X〉, tj. hodnoty konvergují ke střední hodnotě. Stanovíme minimum Xmin
a maximum Xmax z druhé poloviny grafu (tj. pro n ∈ M/2, . . . ,M). Standardní chyba
je přibližně rovna 0.6(Xmax − Xmin). Výsledek je zhruba stejně kvalitní, jako kdybychompoužili deset bloků.
Chybu zapisujeme ve tvaru 12.3 ± 4 nebo v jednotkách posledního místa jako 12.3(4).Ve fyzice se tím rozumí standardní chyba, tj. výsledek je v intervalu (11.9, 12.7) s prav-děpodobností 68 %. V biologii a mezi inženýry je však obvyklejší uvádět chybu na hladiněvýznamnosti 95 %, která je rovna dvojnásobku standardní chyby17. Pokud tedy uvádímeveličinu s chybou (a to bychom měli!), musíme také dodat, jakou chybu máme na mysli.
8.2 Mechanické veličiny
8.2.1 Energie a teplota
Již jsme se setkali s kinetickou energií, z níž se počítá (v klasických simulacích) teplota.Vnitřní energie je průměrem z kinetické i potenciální energie,
U = 〈Ekin〉+ 〈Epot〉
Protože kinetická část je jaksi triviální (je stejná i pro ideální plyn), často se udává jen〈Epot〉; říká se jí také reziduální vnitřní energie.
8.2.2 Tlak
Pokud simulujeme za konstantního objemu, počítáme i další základní mechanickou veličinu,totiž tlak. Použijeme vztah (31) vlevo, nespokojíme se však s ideálním plynem a dosadíme zaF celý vzorec (29), jen jsme změnili symbol pro potenciální energii (místo Epot píšeme U).
βP =1
QNV T
∂QNV T
∂V, QNV T =
1N !Λ3N
∫V N
exp [−βU(~r1, . . . , ~rN)] d~r1 . . . d~rN (85)
Potíž výpočtu je v tom, že derivujeme přes objem, který nemáme v integrandu, ale v mezíchintegrace. Pomůžeme si trikem: zavedeme bezrozměrné (přeškálované) souřadnice ~ξi:
~ri = V 1/3~ξi
17Pro hnidopichy více desetinných míst: 1.959964 násobek
53
Integrál (85) je celkem 3N rozměrný, a proto Jacobián transformace je V N (pro jeden vektorplatí d~ri = V d~ξi), a proto
QNV T =∫
1Nexp
[−βU(V 1/3~ξ1, V
1/3~ξN)]V Nd~ξ1 · · · d~ξN
Každý vektor ~ξi se integruje přes jednotkovou periodickou krychli. Za integrálem je součin:po derivaci členu V N dostaneme (po dělení QNV T podle (85)) ideální část stejně jako v (31),u druhého členu zatím jen zderivujeme exponenciálu (a vydělíme QNV T , což vede ke středníhodnotě):
P =N
VkT −
⟨(∂U(V 1/3~ξN)
∂V
)~ξ1,...,~ξN
⟩Toto již je užitečný vzorec, jedině musíme derivaci podle V počítat v simulacích nume-
ricky, např.: (∂U
∂V
)~ξ1,...,~ξN
≈ U(V + ∆V )− U(V −∆V )2∆V
Metoda se nazývá virtuální změna objemu, protože objem měníme jen pro účely výpočtuderivace a se změněným objemem nesimulujeme. Provádíme tedy normální simulaci a občaszměníme objem o +∆V a −∆V , kde ∆V je dostatečně malé. Zároveň se změnou objemupřeškálujeme všechny souřadnice stejným faktorem [(V ±∆V )/V ]1/3 (tj. necháme ξi nezmě-něné), tedy celou konfigurace „nafouknemeÿ nebo „stlačímeÿ, spočteme potenciální energiiU a pak derivaci.
Derivaci můžeme spočítat i analyticky, protože V se vyskytuje ve všech N vektorovýchargumentech U , použijeme vzorec pro derivaci složené funkce mnoha proměnných a dosta-neme součet:(
∂U(V 1/3~ξN)∂V
)~ξ1,...,~ξN
=N∑i=1
13V −2/3~ξi ·
∂U
∂~ri=
13V
N∑i=1
~ri ·∂U
∂~ri= − 1
3V
N∑i=1
~ri · ~fi
kde ~fi je síla na částici i. Součet vpravo se nazývá viriál sil.
Ve skutečnosti tento obecný vzorec nelze jednoduše použít v periodických okrajových podmínkách. Za před-pokladu párové aditivity a dosahu potenciálu menšího než je polovina boxu můžeme odvodit
P = ρkT − 13V
∑i<j
〈riju′(rij)〉 = ρkT +1
3V
∑i<j
〈rijfij〉
kde f je párová síla.
8.3 Entropické veličiny
Veličiny, které v sobě obsahují entropii, což jsou zvláště Helmholtzova energie a chemickýpotenciál, nelze spočítat jako jednoduchou střední hodnotu. Používá se několik přístupů:
• metoda termodynamické integrace, kdy podobně jako v termodynamice integrujeme(mechanické veličiny) přes teplotu a objem (příp. tlak či další veličiny jako „zapojovacíparametrÿ);
54
• metoda vkládání částice, kdy z pravděpodobnosti vložení spočteme chemický potenciál,
• výpočet reverzibilní práce integrací síly,
• metoda lokální hustoty.
8.3.1 Metoda termodynamické integrace
Chceme-li stanovit (např.) Helmholtzovu energii, musíme ji znát v nějakém vhodném bodě,např. pro ideální plyn nebo ve stavu krystalu (ideální Einsteinův krystal). Potom si vzpo-meneme na kurz fyzikální chemie, kde jsme odvodili vzorce(
∂F
∂V
)T
= −P ∂(βF )∂β
= E
Tyto vztahy budeme integrovat numericky. Musíme proto pseudoměřit v mnoha bodech(hodnotách p a T 18) a pak integrovat např. lichoběžníkovým (jsou-li body blízko u sebe)nebo Simpsonovým pravidlem.
Potenciální energie se vyskytuje ve statisticko-termodynamických vztazích násobená in-verzní teplotou β = 1/kBT . Pokud tedy měníme teplotu, je to to samé, jako bychom potenciálznásobili faktorem rovným inverznímu poměru teplot. Myšlenku můžeme rozšířit zavedenímtzv. zapojovacího parametru. Chceme spočítat rozdíl Helmholtzových energií systémů po-psaných potenciálními energiemi U1 a U0; pro názornost si můžete třeba představit ion a jehonenabitou variantu, z rozdílu pak lze spočítat třeba aktivitní koeicient iontu. Obě energiespojíme
U(λ) = U0 + λ(U1 − U0)
Derivace konfiguračního integrálu je
∂Q
∂λ=∫
∂
∂λe−βUd~rN = −β
∫∂U(λ)∂λ
e−βU(λ)d~rN
a rozdíl Helmholtzových energií
F (1) = F (0) +∫ 1
0
⟨∂U(λ)∂λ
⟩λ
dλ = F (0) +∫ 1
0〈U1 − U0〉λdλ
kde simulujeme několik systémů lišících se hodnotou λ (např. ion, ion s polovičním nábojema ion bez náboje).
8.3.2 Widomova metoda vkládání částice
Chemický potenciál je definován jako změna Gibbsovy energie při vratném přidání částice(v chemii molu částic) k systému za konstantní teploty a tlaku,
µi =
(∂G
∂Ni
)T,p,Nj ,j 6=i
18Existují rafinované metody, kdy lze provést sérii simulací pro různé např. teploty a z dat pak spočítathodnoty veličin při libovolné teplotě uvnitř intervalu (multiple histogram reweighting), dále lze sérii simulacíprovádět paralelně s tím, že mezi jednotlivými simulacem občas „přeskakujemeÿ (parallel tempering).
55
Přesněji slovem „přidáníÿ rozumíme přesun ze standardního stavu, dostaneme tedy chemickýpotenciál vzhledem ke standardnímu stavu. Protože změna Gibbsovy energie je rovna práci(jiné než objemové), vyjadřuje chemický potenciál jakýsi energetický obsah dané částice,její schopnost konat práci (přičemž je k dispozici dost tepla z termostatu a dost energiez barostatu, vše převáděno vratně). Tuto myšlenku můžeme převést do formy použitelnév simulacích několika způsoby.
Uvažujme otevřený systém a pro jednoduchost zápisu nikoliv směs, ale čistou látku.Protože se lépe pracuje s konstantním objemem, napíšeme vztah pro změnu Helmholtzovyenergie v otevřeném systému. Platí
dF = −SdT − PdV + µdN
Chemický potenciál je pak dán podobně,
βµ =
(∂(βF )∂N
)V,T
= −(∂ lnZNV T
∂N
)V,T
Počet částic N ale není spojitá veličina, takže nejmenší diference je rovna jedné,
∂ lnZNV T∂N
=lnZN+1,V T − lnZNV T
(N + 1)−N= ln
(ZN+1,V T
ZNV T
)e−βµ =
ZN+1,V T
ZNV T=
1(N + 1)Λ3
QN+1
QN
≈ 1NΛ3
QN+1
QN
Po porovnání s chemickým potenciálem jednoatomového ideálního plynu (vnitřní stupněvolnosti zde pro jednoduchost neuvažujeme), µid = kBT ln(NΛ3/V ), dostaneme
exp(−βµres) = exp[−β(µ− µid)] =1V
QN+1
QN
kde µres = µ − µid je tzv. reziduální chemický potenciál, tj. vztažený k ideálnímu plynu zadané teploty a objemu.
Rozepišme nyní energii systému N + 1 částic jako
UN+1 = UN + Ψ(N)
kde Ψ je energie interakce (N + 1)-ní částice s N částicemi. Konfigurační integrál systémuN + 1 částic je pak
QN+1 =∫
exp(−βUN+1) d~r1 . . . d~rN+1
=∫
exp(−βUN − βΨ) d~r1 . . . d~rN+1
= QN
∫〈e−βΨ〉N d~rN+1
Pod integrálem máme střední hodnotu počítanou v systému N částic, který simulujeme. Po-slední integrál přes d~rN+1 spočteme metodou Monte Carlo (náhodné střílení); objem oblastije∫
1d~rN+1 = V , a proto
exp(−βµres) =1V
∫〈e−βΨ〉Nd~rN+1 = (〈e−βΨ〉N)N+1 (86)
56
Obr. 20:Vlevo: chemický potenciál lze stanovit integrací střední síly naměřené pro částici v různýchpolohách. Vpravo: necháme-li na částice působit vnější pole, rozmístí se tak, aby součet chemického+ vnějšího potenciálu byl konstantní
kde symbol (X)N+1 značí střední hodnotu veličiny X přes polohy (N + 1)-ní částice v ob-jemu V .
Celý postup je tedy založen na normální NV T simulaci (MC nebo MD). Do systémuobčas vložíme (N + 1)-ní částici a spočítáme, o kolik se změnila energie (totiž o Ψ). Vloženáčástice je virtuální (fiktivní, ghost), tedy není součástí simulovaného systému. Z mnohatěchto virtuálních vložení částice spočteme střední hodnotu Boltzmannova faktoru 〈e−βΨ〉Na podle vzorce (86) pak reziduální chemický potenciál.
Metoda selhává pro příliš velké rozpuštěnce v hustých systémech. Existuje varianta po-stupného vkládání. Nejprve vložíme malou částici a stanovíme její chemický potenciál. Paksimulujeme systém s touto malou částicí (již není virtuální, ale reálná). Částici se pokou-šíme zvětšit (přidat náboje ap.) a opět stanovíme střední hodnotu Boltzmannova faktoru;logaritmus této střední hodnoty je příspěvek k chemickému potenciálu.
8.3.3 Reverzibilní práce integrací střední síly
Z termodynamiky víme, že chemický potenciál je roven vratné práci. Pokud budeme částicidržet (molekulu např. za těžiště) pevně na místě, můžeme v průběhu simulace stanovitstřední sílu, která na částici v daném systému působí. Integrací této síly po dráze získámechemický potenciál resp. jeho změnu mezi dvěma stavy
∆µi = −∫ ~ri(2)
~ri(1)〈~fi〉 · d~ri
kde ~fi je síla působící na částici i. Nutno provést sérii simulací pro různé polohy částice.
8.3.4 Metoda lokální hustoty/koncentrace
Nechť na rozpuštěnce i působí vnější potenciál U exti (~r) (např. „gravitaceÿ). V rovnováze je
součet chemického a vnějšího potenciálu konstantní
µi(~r) + U exti (~r) = const
čili∆µi = µi(~r2)− µi(~r1) = −[U ext
i (~r2)− U exti (~r1)]
57
Při výpočtu v simulaci počkáme na ustavení rovnováhy a potom naměříme histogram hustotči koncentrací částice i. Je-li např. potenciál U(~r1) tak vysoký, že koncentrace v této polozeje velmi nízká a aktivitní koeficient jednička, snadno vypočteme chemický potenciál rozpuš-těnce. Ze simulace pak získáme chemický potenciál (a tedy i aktivitní koeficient) v oblastivysokých koncentrací.
8.4 Strukturní veličiny
Základní strukturní veličinou ve fyzice tekutin je radiální distribuční funkce (RDF, též párovákorelační nebo distribuční funkce). Vyjadřuje pravděpodobnost nalezení částice ve vzdále-nosti r od jiné částice. Tato pravděpodobnost je normalizovaná tak, že pro nekorelovanéčástice (idelní plyn) je hodnota RDF rovna jedničce (a proto také pro kapaliny je RDF proobě částice daleko od sebe rovna 1). Radiální distribuční funkci dostaneme inverzní Fouriero-vou transformací rozptylových dat. Nejčastěji se používá rozptyl neutronů nebo rentgenovazáření.
Mějme N molekul ideálního plynu v objemu V . Jaká je pravděpodobost, že naleznuv objemu dV částici (tj. kteroukoliv z N částic)? Je rovna NdV/V . Jaká je pravděpodobost,že naleznu částici v objemu dV ve vzdálenosti r od vybrané částice (červeně na obr. 21vlevo)? Je rovna (N − 1)dV/V , protože červenou částici už neuvažujeme a protože částiceideálního plynu se neovlivňují; a ovšem můžeme pro N velké napsat jen NdV/V . Jinak tomuale je v husté kapalině. Pravděpodobnost nalezení částice závisí na vzdálenosti (např. velmiblízko nebude už žádná, protože molekuly se nepřekrývají). Tento vliv vyjádříme faktoremg(r). Pro velké vzdálenosti v kapalině platí, že g(r) → 1 (stejně jako v ideálním plynu),protože molekuly již o sobě nevědí.
Radiální distribuční funkce (RDF) jednoduchých tekutin (se sféricky symetrickými mo-lekulami) odrážejí slupkovou strukturu, viz obr. 22 vlevo. Kolem každé molekuly je mírněnepravidelná „slupkaÿ molekul, v podstatě ve vzdálenosti odpovídající minimu potenciálu(kromě velkých tlaků). Této slupce odpovídá první vrchol na RDF. Slupka prvních sousedůponěkud brání molekulám v přístupu, a proto následuje minimum. Druhé slupce odpovídájiž méně významné maximum. Naproti tomu struktura vody je ovlivněna tetraedrickýmuspořádáním vodíkových vazeb, a proto druhé maximum je blíže, viz obr. 22 vpravo.
Obr. 21: Radiální distribuční funkce vyjadřuje, kolikrát je v průměru víc částic ve vzdálenosti rokolo vybrané částice než by tomu bylo v ideálním plynu
58
Obr. 22: Vlevo: Radiální distribuční kapalného argonu, vpravo: distribuční funkce modelu kapalnévody
Matematicky je RDF izotropní a homogenní tekutiny dána vztahem
g(r) ≡ g(r12) =N(N − 1)ρ2QNVT
∫. . .
∫exp[−βU(~r1, ~r2, . . . , ~rN )] d~r3 . . . d~rN
Stejně jako Boltzmannův faktor exp[−βU(~r1, ~r2, . . . , ~rN )]/QNV T udává pravděpodobnost, že částice 1. . .N na-lezneme v daných polohách, po integraci přes d~r3 . . . d~rN dostaneme pravděpodobnost, že částice 1 a 2 majípolohy ~r1, ~r2 pro všechny možné polohy ostatních částic. Protože kapalina je homogenní a izotropní (na rozdíltřeba od krystalu), můžeme napsat g(~r1, ~r2) = g(r12) jen jako funkci vzdálenosti. Faktor N(N − 1)/ρ2 zajišťuje,že g(r) = 1− 1/N pro ideální plyn.
Radiální distribuční funkci v simulacích počítáme tak, že projdeme všechny páry částica započteme jedničku do přihrádky histogramu odpovídající určitému rozsahu vzdáleností,dejme tomu [i∆r, (i+ 1)∆r), kde ∆r je malé.
Často je názornější veličina zvaná kumulativní koordinační číslo (kumulativní RDF). Jedefinována vztahem
N(r) = ρ
∫ r
0g(r′)4πr′2dr′
Integruje se přes element objemu dV = 4πr′2dr′ (povrch koule 4πr′2 krát dr′). Počet částicideálního plynu v tomto objemu by byl dV N/V , v kapalině pak g(r)dV N/V . Veličina N(r)vyjadřuje proto průměrný počet molekul do vzdálenosti r. Často se integruje do vzdálenostiodpovídající prvnímu minimu na RDF, výsledné koordinační číslo je rovno průměrnémupočtu molekul v první slupce. Např. u vody jsou v první slupce necelé čtyři molekuly, alev průměru nalezneme jen dva vodíky okolo kyslíku, viz obr. 23.
8.5 Transportní vlastnosti
Transportem v nerovnovážné termodynamice rozumíme tok způsobený jistou termodyna-mickou silou (která se dá vyjádřit jako gradient chemického potenciálu). Např. tok teplaje způsoben rozdílem teplot, termodynamickou silou je gradient teploty (který lze převéstna gradient chemického potenciálu molekul za různých teplot); elektrický proud19 je způso-ben elektrickým polem (gradientem elektrochemického potenciálu), tok hybnosti gradientem
19Způsobený iontovou vodivostí, protože pohyb elektronů v kovech nelze klasickými simulačními metodamistudovat
59
0 1 2 3 4 5 6
r/A
0
1
2
3
4
5
6
g(r
) n
eb
o N
(r)
NHO(r) (H okolo O)
gHO(r)
NOO(r)
gOO(r)
voda
Obr. 23: Simulační distribuční funkce modelu kapalné vody a kumulativní koordinační čísla. Na-lezneme minimum na gOO(r) (přibližně 3.2 A) a odečteme hodnotu NOO (cca. 3.6). Pouze vodíkypatřící k jiné molekule jsou započteny do gHO(r), a proto je koordinační číslo těsně pod 2, jinak bybylo také skoro 4
rychlosti (je příčinou viskozity). Tyto jevy nějak souvisejí s rychlostí částic, a proto to jsouz hlediska simulací kinetické jevy. Obecně je tok vyjádřen vektorem, který vyjadřuje, kolikveličiny proteče jednotkovou plochou (kolmo) za jednotku času, a pro malé toky předpoklá-dáme lineární úměrnost:
~J = −ξ ~F
kde F je termodynamická síla a ξ příslušný transportní koeficient.Tyto veličiny nelze zapsat jako střední hodnotu funkce závislé jen na polohách (jako
v (30)), ale potřebujeme i rychlosti. Proto k jejich výpočtu potřebujeme molekulovou dy-namiku. Můžeme použít „normálníÿ rovnovážnou simulaci a vzorec (obsahující rychlosti),tento postup se označuje zkratkou EMD (equilibrium molecular dynamics). Tím se budemezabývat v dalším textu.
Můžeme také použít „fyzikálníÿ přístup: příslušnou sílu zkrátka zapneme a sledujeme tok;např. zapneme elektrické pole a měříme proud. Tento postup se označuje zkratkou NEMD(nonequilibrium molecular dynamics). Výhodou je koncepřní jednoduchost. Nevýhodou je, žesystém není v rovnováze, a proto v systému dochází k disipaci energie a zahřívá se. Pro většísíly není odezva lineární, pro menší síly je odezva malá a těžko rozlišitelná od šumu. Protomusíme provést několik simulací pro různé hodnoty termodynamické síly a extrapolovatk nule.
Postupy EMD si probereme na příkladu difuze.
8.5.1 První Fickův zákon
Představte si, že nalijete do sklenice sirup a opatrně převrstvíte vodou. I když vodu sesirupem nezamícháte, časem dostane homogenní směs – sladkou obarvenou vodu. Příčinouje difuze, čili pohyb molekul z místa s nižší koncentraci do místa s vyšší koncentrací.
Rychlost difuze popíšeme vektorem ~Ji, který vyjadřuje množství látky, které protečejednotkovou plochou (kolmo ke směru vektoru) za jednotku času. Pokud množství látky
60
vyjádříme v molech, bude jednotkou toku v SI mol m−2 s−1 a koncentrace mol m−3; prohmotnostní vyjádření jsou jednotky kg m−2 s−1 a kg m−3.
Příčinou difuzního toku je gradient koncentrace, který zapisujeme různými symboly
~∇ci = grad ci =
(∂
∂x,∂
∂y,∂
∂z
)ci =
(∂ci∂x
,∂ci∂y
,∂ci∂z
)(∇ se nazývá „nablaÿ). Gradient koncentrace je vektor, který míří ve směru zvyšování kon-centrace a jehož velikost je úměrná tomu, jak rychle se koncentrace se vzdáleností zvyšuje.
Tok ~Ji látky i je v nejjednodušším případě úměrný gradientu koncentrace,
~Ji = −Di~∇ci
Veličina D se nazývá koeficient difuze (též difuzivita). Její jednotka je m2 s−1, ať používámejednotky založené na látkovém množství nebo hmotnosti.
Příklad. Trubice tvaru U délky l = 20 cm a průřezu A = 0.3 cm2 má na obou koncích fritu. Jedenkonec je ponořen v Coca-Cole (11 hm.% cukru) a druhý v čisté vodě. Kolik cukru prodifunduje zaden? Dsacharoza(25 C) = 5.2·10−6 cm2 s−1.
Řešení. 110 g cukru v litru odpovídá hmotnostní koncentraci cw = 110 g dm−3 = 110 kg m−3. Gra-dient koncentrace je grad cw = cw/l = 550 kg m−4. Difuzivitu převedeme na SI:D = 5.2·10−6 cm2 s−1
= 5.2·10−10 m2 s−1. Tok je tedy J = D grad cw = 2.86·10−7 kg m−2 s−1. Po znásobení plochouA = 0.3 cm2 a časem 1 den vyjde množství cukru jako
m = JAt = 2.86·10−7 kg m−2 s−1 × 0.3·10−4 m2 × 24× 602 s = 7.4·10−7 kg = 0.74 mg
8.5.2 Einsteinova–Smoluchowského rovnice
Zkusme se na empirický první Fickův zákon podívat z hlediska termodynamické síly a mo-lekul. Tok látky je dán střední rychlostí molekul ~vi:
~Ji = ~vici
Termodynamická síla je gradientem chemického potenciálu20:
~Fi = −~∇(µiNA
)= −kBT
ci~∇ci
kde jsme použili vztah pro ideální roztok, µi = µi + RT ln(ci/cst), a ~∇ ln ci = (1/ci)~∇ci(derivace složené funkce).
Pohybuje-li se molekula rychlostí ~vi, působí na ní síla odporu prostředí (přibližně) úměrnárychlosti:
~F třeníi = −λi~vi
20Pamatuj: rozdíl chemických potenciálů = reverzibilní práce, kterou částice vykoná při pohybu z jednohomísta do druhého
61
Obr. 24: Ploškou dydz proteče do krychličky zespoda Jx(x)dydzlátky za jednotku času, horem odteče Jx(x + dx)dydz. RozdílJx(x+ dx)dydz − Jx(x)dydz = (∂Jx/∂x)dxdydz
kde λi je koeficient tření. V ustáleném stavu jsou obě síly v rovnováze,
~F třeníi + Fi = 0 tj. λi~vi = λi
~Jici
= −kBT
ci~∇ci
Odvodili jsme první Fickův zákon, jen konstantu nazýváme jinak. Porovnáním s Fickovýmzákonem ve tvaru ~Ji = −Di
~∇ci dostaneme hledaný vztah mezi (makroskopickou) difuzivitoua (mikroskopickým) koeficientem tření molekuly v prostředí, tzv. Einsteinovu rovnici (téžEinsteinovu–Smoluchowského):
Di =kBT
λi(87)
Pro velké kulovité molekuly o poloměru Ri v kapalině o viskozitě η platí Stokesův vzorecpro koeficient tření, λi = 6πηRi, a proto ~Fi = 6πηRi~vi, a Einsteinova rovnice přejde naEinsteinovu–Stokesovu rovnici
Di =kBT
6πηRi
Pro malé molekuly platí jen řádově, protože kapalina není v tomto rozlišení kontinuum amolekula nemusí být kulatá, přesnější je tato rovnice pro větší kulovité koloidní částice.
Příklad. Odhadněte velikost molekuly sacharozy. Viskozita vody je0.891·10−3 m−1 kg s−1 při 25 C.Řešení. Obrácením Einsteinovy–Stokesovy rovnice dostaneme
Uvažujme krychličku dx × dy × dz. Tok ploškou dydz v poloze x (viz obr. 24) je malinkojiný než v poloze x+ dx, to samé platí pro ostatní dva směry. Po sečtení dostaneme změnulátkového množství za jednotku času (rovnici kontinuity)
dndt
=
(∂Jx∂x
)dxdydz+
(∂Jy∂y
)dxdydz+
(∂Jz∂z
)dxdydz = −D
[∂2c
∂x2+∂2c
∂y2+∂2c
∂z2
]dxdydz
Protože n/(dxdydz) = c, dostaneme tzv. druhý Fickův zákon
∂ci∂t
= −Di∆ci (88)
kde ∆ = ~∇ · ~∇ = ∂2
∂x2+ ∂2
∂y2+ ∂2
∂z2se nazývá Laplaceův operátor. Matematici znají tento
typ rovnice pod termínem rovnice pro vedení tepla (protože teplo se šíří podle stejnérovnice).
62
0
x
0
200
y
0
x
0y
-3 -2 -1 0 1 2 3
x
0.1
0.2
0.3
0.4
0.5 t=1
t=2
t=3
t=4
t=5
c(x
,t)
Obr. 25: Brownův pohyb. Vlevo: tři trajektorie začínající v bodě (0, 0); koncové body jsou vyzna-čeny zeleným kolečkem. Uprostřed: koncové body 10 000 trajektorií; Vpravo: závislost koncentracev ose x na čase (c(0, 0) =∞) resp. pravděpodobnosti výskytu (částice je v počátku pro t = 0)
8.5.4 Einsteinův vztah
Sledujme nyní trajektorii jedné částice během určitého času, viz obr. 25. Představme si, žeten samý experiment provedeme mnohokrát a vždy zaznamenáme koncový bod. Dostaneme„mrakÿ bodů. Můžeme se ptát, jak jsou tyto body rozloženy, neboli jaká je hustota pravdě-podobnosti, že se za daný čas částice dostane do daného místa. Hustota pravděpodobnostije ale (skoro) to samé co koncentrace, platí tedy pro ni rovnice vedení tepla (88) s tím, žev počátečním čase je částice se 100% pravděpodobností v počátku; je tedy popsána „funkcíÿ(přesněji distribucí), jejíž integrál je jednička, je nekonečno v počátku a nula jinde. Takové„funkciÿ se říká Diracova delta-funkce. Snadno se přímým derivováním ověří, že
1D: c(x, t) = (4πDt)−1/2 exp
(− x2
4Dt
)řeší (88) v jedné dimenzi a
3D: c(~r, t) = (4πDt)−3/2 exp
(− r2
4Dt
)ve třech dimenzích. Tyto funkce jsou normalizované (pro každé t); např. v 1D:∫ ∞
−∞c(x, t) dx = 1
K ověření potřebujete znát vzoreček (Gaussův integrál)∫∞−∞ exp(−x2) dx = π1/2, po substi-
tuci∫∞−∞ exp(−Ax2) dx = (π/A)1/2. Pak také snadno vypočtete
1D: 〈x2〉 =∫x2c(x, t)dx = 2Dt, 3D: 〈r2〉1/2 =
√6Dt (89)
Tento výsledek je významný. Znamená, o kolik se v průměru (lépe řečeno „kvadratickém prů-měruÿ) posune difundující částice za čas t. Za čtyřnásobek času nalezneme částici v průměrudvakrát tak daleko.
63
0 2 4 6 8 10
t/ps
0
0.02
0.04
0.06
0.08
0.1
MS
D(t
)/n
m2
Cl-Na+
0 2 4 6 8 10
t/ps
0
1
2
3
MS
D(t
)/a
rb.u
.
whole box
Kohlrausch Na+ + Cl-
Obr. 26: Vlevo: střední kvadratické posunutí MSD pro oba ionty v tavenině NaCl; Vpravo: veličinaMSCD pro určení vodivosti a součet příspěvků na základě Kohlrauschova zákona (vodivosti získámz difuzivit a jako nezávislé sečtu)
Tento výsledek lze ověřit i pomocí modelu náhodné pocházky. Předpokládejme jen jednu dimenzi (osu x).Nechť se za jistý časový krok h částice posune náhodně (v důsledku náhodných nárazů ostatních molekul) o∆x = +(2D)1/2 (s pravděpodobností 1/2) a o ∆x = −(2D)1/2 (s pravděpodobností 1/2). Z tzv. centrálnílimitní věty známé z matematické statistiky plyne, že po mnoha krocích dostanu to samé Gaussovo rozloženíc(x, t) jako výše. Pro uvedené pravděpodobnosti pohybu doleva a doprava dostaneme, že za čas 2h se částiceposune o −2∆x s pravděpodobností 1/4, o 0 s pravděpodobností 1/2 (může se vrátit zprava i zleva) a o 2∆xs pravděpodobností 1/4. Obecně ji v čase nh najdeme v poloze (n− 2k)∆x, k ∈ 0, . . . , n, s pravděpodobností(nk
)2−n. Po troše matematiky dostaneme v limitě n→∞ opět Gaussovo rozložení.
Vzorec (89) je vhodný k přímému výpočtu difuzivity v simulacích. Provedeme ještě prů-měrování přes všechny částice (nechť je jich N stejného druhu) a spočteme veličinu
MSD(t) =1
6N
N∑i=1
[~r(0)− ~r(t)]2 (90)
případně ji ještě zprůměrujeme přes několik simulací (či bloků vybraných z jedné simulace).Veličinu vyneseme do grafu, obr. 26; difuzivita je směrnice lineární části křivky (po vynechánípřechodových jevů na začátku).
Analogický vztah existuje pro iontovou vodivost. Místo jednotlivých poloh částic ale musíme uvažovat polohyvážené náboji,
MSCD(t) =16
∑i
qi[~ri(t)− ~ri(0)]
2(91)
Směrnice křivky je rovna kBTV κ, kde κ je měrná vodivost. Zde je součet přes ionty unvitř druhé mocniny, protoževodivost je funkcí celé simulační buňky, nikoliv jedné částice. Musíme tedy průměrovat přes několik simulací(bloků), protože bychom neměli dostatečnou statistiku.
64
8.5.5 Greenovy–Kubovy vzorce
Existuje ještě jeden přístup, tzv. Greenovy–Kubovy vzorce. Zpravidla se odvozují na základěpředstavy malé poruchy, která ovlivňuje systém vyvíjející se v čase. Např. pro vysvětlenídifuze se uvažuje malá síla působící na částici a zapnutá v určitém čase, pro vodivost zapnuelektrické pole ap. Matematickým zpracováním v rámci tzv. teorie lineární odezvy (kteráleží mimo rozsah těchto skript) pro difuzi vyjde
D =13
∫ ∞0〈~ri(0) · ~ri(t)〉dt (92)
Dá se ukázat, že vzorce (92) a (90) jsou ekvivalentní.
Zkusme to pro jednu souřadnici. Integrál totiž snadno vypočteme:
D = limt→∞
∫ t
0〈x(0)x(t′)〉dt′ = lim
t→∞
[〈x(0)x(t′)〉
]t0
Z důvodu časové symetrie pohybových rovnic můžeme provést transformaci t → −t (znaménko změní takéderivace: x(0)→ −x(0))[
Věc v integrálu (92) se po normalizaci nazývá rychlostní autokorelační funkce a je zají-mavá sama o sobě:
cv(t) =〈x1(0)x1(t)〉〈x1(0)x1(0)〉
Má následující vlastnosti:
• Je sudá (symetrická v čase), cv(t) = cv(−t). To vyplývá ze symetrie pohybových rovnic.
• Platí cv(0) = 0 (tečna v počátku je rovnoběžná s osou x). To je přímým důsledkempředchozího tvrzení a existence derivací.
• Pro krátké časy je blízká 1, protože částice za krátký čas nestačí změnit (působenímostatních částic) směr letu.
• Po jisté době překmitne do záporných hodnot. To je důsledkem odrazu částice odostatních. Odraz však není dokonalý, protože ostatní částice se náhodně pohybují, aproto hodnota nedosahuje −1.
• Platí limt→∞ cv(t) = 0, tj. po dlouhé době částice zapomene, jakým směrem letělav čase t = 0.
65
0 0.5 1 1.5
t/ps
0
0.5
1
cv(t)
Obr. 27: Rychlostní autokorelační funkce kapalného argonu: teplota 150 K, hustota 1344 kg m−3,teplota 120 K, hustota 1680 kg m−3. Výsledky trajektorie dlouhé 100 ps pro 216 Lennard-Jonesovýchčástic
![147.33.74.135147.33.74.135/knihy/uid_isbn-978-80-7080-747-7/978-80-7080-747-7.pdf · +HMQD%RKGDQ ,QIRUPDþQtWHUPRG\QDPLND, VWHG 9\VRNiãNRODFKHPLFNR WHFKQRORJLFNiY3UD]H 3UDKD ,6%1](https://static.dokumenty.site/doc/80x56/5c776aeb09d3f2c43b8c0642/14733741351473374135knihyuidisbn-978-80-7080-747-7978-80-7080-747-7pdf.jpg)
