Monte Carlo simulace Monte Carlo ● Stochastické metody používající (pseudo)náhodná čísla ● Typické využití při výpočtech integrálů (vícerozměrných) ● Využití v biologii, chemii, matematice, finančnictví, počítačové grafice ● Hledáme střední hodnotu veličiny, která je výsledkem náhodného jevu. Monte Carlo - historie ● John von Neumann, Stanislaw Marcin Ulam ● Použití v projektu Manhattan – vývoj atomové bomby ● Studium chování neutronů v různých prostředích – Značné množství experimentálních dat – Nemožnost teoretického ani praktického řešení – Z experimentů věděli, že každý stý neutron je pohlcen Monte Carlo - historie ● Inspirace ruletou – Vybrali jedno políčko rulety, které symbolizovalo pohlcení – Roztočení rulety symbolizovalo pohyb neutronu – Pokud kulička dopadla na označené políčko, byl neutron pohlcen, u ostatních políček neutron pokračoval ve své dráze – Klasické hraní rulety bylo ale neefektivní – zapojení počítače ENIAC Monte Carlo - historie ● Nepočítačovým příkladem MC integrace je Buffonova jehla (1777) ● Zadání úlohy zní takto: Na podlaze leží list papíru, na kterém jsou nakresleny rovnoběžné linky stejně vzdálenými. Na papír házíme jehlu, jejíž délka je stejná jako vzdálenost dvou linek. Počítáme pravděpodobnost, že jehla po dopadu protne některou z linek. Buffonova jehla Buffonova jehla ● Vypočtená pravděpodobnost je rovna l / πd ● Provedením velkého počtu experimentů lze odhadnout hodnotu čísla π ● Při 5000 hodech vypočítal v roce 1850 Volf číslo π = 3,1596 Matematické zákony pro MC ● Zákon velkých čísel – aritmetický průměr n náhodných veličin se stejnou střední hodnotou se s rostoucím n za určitých předpokladů blíží této střední hodnotě. ● Centrální limitní věta – rozdělení výběrového průměru se po vhodné normalizaci blíží k normálnímu rozdělení. Monte Carlo – obecný postup ● Nejčastějším postupem MC metody je modelování náhodné veličiny X takové, že její střední hodnota E(X) je rovna hledané hodnotě a E(X) = a ● Provedeme-li n nezávislých realizací X1 , … , Xn náhodné veličiny X, můžeme odhadnout hodnotu a pomocí aritmetického průměru a= 1 n (X1+ X2+...+ Xn) Monte Carlo – obecný postup Otázky při použití MC metody 1) Jak vybrat nejvhodnější náhodnou veličinu pro určení hledané hodnoty? 2)Jak najít hodnoty použité veličiny? Nejprve nalezneme hodnoty v intervalu (0,1), které transformujeme funkcí na hledané hodnoty a s těmi pak pracujeme. Monte Carlo – obecný postup 1) Generování náhodných čísel yi s rovnoměrným rozdělením na intervalu (0, 1). 2) Transformace náhodných čísel yi na náhodná čísla zi se složitějším rozdělením. 3) Pomocí náhodných čísel zi se buď přímo počítají odhady charakteristik náhodné veličiny X nebo se počítají pomocí vhodného algoritmu hodnoty xi a odhady charakteristik náhodné veličiny. 4) Získané výsledky se zpracují statistickými metodami. Příklad ∫ 0 1 sin ( 1 x )dxVypočítejte integrál n = 1000000000 SUM = 0 FOR i = 1 TO n SUM = SUM + sin(1/RAND(0,1)) END PRINT SUM / n0 Výsledek: 0,504086, doba výpočtu 46 s, Generování náhodných čísel Podmínkou úspěchu je možnost náhodný proces mnohonásobně opakovat a k tomu je třeba mít k dispozici dostatečný počet náhodných čísel s nejrůznějšími distribučními funkcemi. Pro každou distribuční funkci není zapotřebí tvořit speciální generátor náhodných čísel, který odpovídá dané distribuční funkci. V praxi si vystačíme s generátorem náhodných čísel s rovnoměrným rozdělením na intervalu (0,1). Náhodné veličiny s jiným typem rozdělení lze potom získat vhodnou transformací. Generátor náhodných čísel ● Generátor náhodných čísel je výpočetní nebo fyzické zařízení, které generuje řadu náhodných čísel, tj. čísel která postrádají jakýkoliv vzor (sekvenci, předvídatelnost). ● Generátory náhodných čísel fungují tak, že nejprve vygenerují posloupnost náhodných čísel s rovnoměrným rozdělením (primární generátor) a poté je z nich transformací vytvořena posloupnost s požadovaným rozdělením. Generátory náhodných čísel ● Fyzikální (hrací kostka, mince). Častěji se ale používají generátory založené na kombinaci radioaktivního zářiče a detektoru (vyzařuje částice v náhodném množství). ● Pseudonáhodné generátory náhodných čísel Generátor pseudonáhodných čísel je efektivní deterministický program, který generuje posloupnost čísel, statistickými testy pokud možno nerozlišitelnou od náhodné. Pseudonáhodné generátory Vstupními daty pro pseudonáhodné generátory jsou skutečně (téměř) náhodné (pokud probíhá útok postranním kanálem navíc záměrně ovlivňované), leč krátké, posloupnosti zvané random seed, které jednoznačně určují další běh programu (generátoru). V důsledku determinističnosti těchto programů jsou na počítači s ohraničenou pamětí nevyhnutelně periodické, tedy po určité době (periodě) se generovaná posloupnost začne opakovat. Ta však může být velmi dlouhá, tudíž nedetekovatelná. Typy generátorů ● Mersenne Twister je generátor pseudonáhodných čísel vyvitutý Makoto Matsumotem a Takuji Nishimurou. Tento generátor nevyužívá aritmetických operací ve smyslu sčítání, odečítání, násobení či dělení. Využívá pouze bitového posunu a funkcí and, or a xor. Periodatohoto generátoru je 219937 - 1. Tento generátor je díky použitým operacím velice rychlý. Úspěšně prošel řadou testů statistické náhodnosti a je hojně využíván. Typy generátorů ● Lineární kongruentní generátor patří mezi nejjednodušší generátory. Používá matematické operace. ● Například lineární kongruentní generátor jazyka C++ má periodou 2147483647. Typy generátorů ● Subtract with carry je zpožděný Fibonacciho generátor. ● Teorie z hlediska matematiky je u tohoto generátoru zatím nejasná. Vlastnosti jsou odvozeny pouze ze statistických výsledků, ve kterých vychází lépe než Lineární kongruentní generátor. Tento generátor je jeden ze tří standardních generátorů knihovny jazyka C++. Rozšíření Monte Carlo metody V praxi můžeme narazit na problém, kdy náhodné generování stavů v mnoha případech povede k nesmyslným stavům a to v převládajícím množství. Potom bychom i při velkém počtu provedených experimentů nedošli ke správným výsledkům. Příkladem může být rozmístění n tuhých koulí (atomů) do zadaného objemu. Pokud dojde k překryvu alespoň dvou koulí, je takový stav nemožný. Možné stavy jsou pouze v případě, že se žádné dvě koule nebudou překrývat. Proto je třeba klasickou MC „vylepšit“. V takovém případě vyjdeme z vhodného stavu a změnou konfigurace tohoto stavu dojdeme k další konfiguraci. Markovův řetězec Markovův řetězec popisuje obvykle diskrétní náhodný (stochastický či pravděpodobnostní) proces, pro který platí, že pravděpodobnosti přechodu do následujícího stavu závisejí pouze na současném stavu, ne na předchozích stavech. Tato tzv. markovovská vlastnost dovoluje proces znázornit stavovým diagramem, kde z každého stavu (uzlu grafu) vycházejí hrany možných přechodů do dalšího stavu s připsanou pravděpodobností. Markovův řetězec ● Pravděpodobnost stavu v bodě 2 vychází z pravděpodobnosti v bodě 1. Přechod ze stavu 1 do stavu 2 je dán maticí přechodu W. – Můžeme tedy psát: π2 = π1 . W ● Krok 3 závisí na kroku 2 a již nezávisí na kroku 1. ● Matematicky: π3 = π2 . W = π1 . W .W Příklad Máme systém ve dvou stavech A, B. Pravděpodobnost, že systém A přejde do systému B je 90 %, pravděpodobnost, že systém zůstane v tomto stavu je tedy 10 %. Pravděpodobnost, že systém přejde z konfigurace B do konfigurace A je 30 %, pravděpodobnost, že zůstane ve stejné konfiguraci je 70 %. Metropolisův algoritmus V reálném případu ale řešíme problém s více než dvěma stavy. Proto musíme použít ještě sofistikovanější metody. Metropolisův algoritmus vychází z metod Markovových řetězců a hodí se právě pro analýzu systémů s více než dvěma možnými stavy. Metropolisův algoritmus ● Celý algoritmus můžeme popsat pomocí 4 následujících kroků – Vybereme libovolnou částici – Náhodně změníme její polohu – Vypočteme změnu potenciální energie systému – Pokud je nový systém energeticky vyšší, zůstaneme v původní konfiguraci. Jinak přejdeme do nové konfigurace. ● Uvedené kroky opakujeme mnohokrát. Simulované žíhání Tato metoda je úpravou Metropolisova algoritmu. Je určena k hledání globálního minima funkce. Stejně jako u Metropolisova algoritmu pro výpočet energie figuruje teplota systému. Ta je v případě Metropolise konstatní po celou dobu simulace. V případě simulovaného žíhání je na počátku teplota nastavena na vysokou hodnotu, což umožní přijetí téměř každé konfigurace systému, neboť je velmi málo konfigurací zamítnuto. Postupně dochází ke snižování teploty a tím jsou vybírány konfigurace s nižší energií. Na konci simulace je už přijata pouze jedna energeticky nejvýhodnější simulace. Aplikace MC v chemii Monte Carlo metoda je snadno implementovatelná na atomární systém, protože zde bereme v úvahu translační stupně volnosti. (Atomy bereme jako rigidní koule.) Algoritmus je snadno implementovatelný a výsledky získáme pro relativně malý počet simulačních kroků (desítky tisíc). Praktické problémy vyvstanou při aplikaci MC na molekulární systémy, zvláště pak na flexibilní molekuly s mnoha stupni volnosti. Zde je třeba zahrnout vnitřní stupně volnosti molekul, které musíme měnit. To vede v mnoha případech ke konfiguracím s velmi vysokou energií, které bývají zamítnuty. Musíme tedy provádět mnohonásobně větší počet MC simulací. Simulace rigidních molekul ● V případě rigidních nesférických molekul musíme brát v úvahu jak jejich orientaci, tak jejich vzájemnou pozici v simulovaném prostoru. V praxi se to provádí tak, že molekulu v prostoru posuneme a také s ní otočíme v průběhu jednoho MC simulačního kroku. Posun je popsaný pomocí změny polohy těžiště molekuly. Pro popis otočení molekuly máme několik možností. Připomeňme si matice posunu, zrcadlení, rotace, o kterých jsme se bavili v předchozích přednáškách. ● Například můžeme zvolit kombinaci rotací kolem souřadných os. Zahrnutí rotace společně s translací vede k většímu počtu zamítnutých konfigurací a je tedy třeba provádět mnohem větší počet MC simulací, abychom splnili požadovaná kritéria modelu. MC simulace flexibilních molekul MC simulace flexibilních molekul jsou, při použití klasické MC, možné pouze pro malé molekuly. Pro simulaci velkých flexibilních systémů musíme využít určitá „vylepšení metody“, jako například zafixování určitých rigidních částí molekuly. Nejjednodušší cesta změn poloh atomů zde není použitelná, proto je vhodné použít pro popis molekuly tzv. vnitřních souřadnic. Zde opět si znovu připomeneme, že i malá změna torzního úhlu uprostřed molekuly může vést k rozsáhlým změnám souřadnic atomů vzdálených od atomů rotující torze. Tato změna často vede k překryvu atomů v molekule a prudkému vzrůstu energie systému, což vede k zamítnutí změny. Zde platí přímá úměrnost mezi počtem stupňů volnosti systému a počtem Monte Carlo kroků. MC simulace polymerů Polymery jsou makromolekuly, které se skládají z opakujících se jednotek (monomerů), chemicky navzájem vázaných do molekuly. Některé polymery jsou tvořeny řetězci jednoho monomeru (polyethylen, polystyren). Jiné polymery jsou tvořeny směsí různých monomerů, příkladem jsou například bílkoviny, které jsou tvořeny dvaceti základními aminokyselinami. Spojení různých jednotek vede k rozdílným vlastnostem v dané části polymeru. Každá jednotka monomeru má své specifické konformační vlastnosti, což ovlivňuje celkové chování molekuly. Zmapování chování polymeru za různých podmínek je důležité pro praktické využití. Protože se jedná o simulace velkých systémů s mnoha stupni volnosti, nemůžeme pro generování použít klasické MC metody. Pro výpočet energií stavů jsou použity empirické modely založené na zjednodušených aproximacích, narůstající výpočetní výkon počítačů nám sice umožňuje studovat stále větší systémy, ale stále je třeba využívat další „smartness“. MC simulace polymerů ● Model mřížkové stavby polymeru – Do jednotlivých bodů na mřížce (2D, 3D)umísťujeme monomerní jednotky s popsanými vlastnostmi a hodnotíme, zda konfiguraci jednotek přijmeme nebo zamítneme. Jednou ze základních podmínek přijetí je, že generovaný řetězec se nesmí na mřížce sám protnout. Tento způsob je velice rychlý a efektivní a dokáže vygenerovat a ohodnotit velké množství stavů pro velké systémy. MC simulace polymerů ● Kontinuální model polymeru – staví polymer postupně z jednotek (můžeme si jej představit jako navlékání korálků na šňůru). Přidané jednotky (korálky) interagují se svými navázanými sousedy, ale také se svým okolím (nevazebné interakce). Zde ale opět narážíme na značnou různorodost stavů díky počtu volně rotovatelných vazeb, a proto je nutná aplikace dalších algoritmů, které umožní řešení problému. MC simulace Narůstající výpočetní výkon současných počítačů nám umožňuje provádět simulace velkých systémů za přítomnosti molekul solventu i za přítomnosti iontů. Při těchto simulacích se běžně používá kombinace molekulové dynamiky s Monte Carlo metodou. Ta je použita pro generování a optimalizaci pozic molekul solventu a iontů. Při této optimalizaci je studovaná molekula polymeru držena v tzv. zamrzlé pozici. Volba metody, použité aproximace, použití různých algoritmů vždy záleží na tom, jakou vlastnost molekul chceme studovat.