Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Modelowanie molekularne - metody Monte Carlo Marcin Buchowiecki 24 lutego 2010 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Spis treści 1 Wprowadzenie 2 Pola siłowe 3 Klasyczna metoda Monte Carlo 4 Kwantowe metody Monte Carlo Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wykład 1. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Monte Carlo Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Monte Carlo Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Stanisław Ulam Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Stanisław Ulam Pomysł ten, nazwany pozniej metodą Monte Carlo, wpadł mi do głowy, kiedy podczas choroby stawiałem pasjanse. Zauważyłem, że znacznie praktyczniejszym sposobem oceniania prawdopodobieństwa ułożenia pasjansa jest wykładanie kart, czyli eksperymentowanie z tym procesem i po prostu zapisywanie procentu wygranych, niż próba obliczenia wszystkich możliwosci kombinatorycznych, których liczba rośnie wykładniczo Jest to zaskakujace z intelektualnego punktu widzenia, i choć może nie całkiem upokarzające, to jednak zmusza do skromnosci i pokazuje granice tradycyjnego, racjonalnego rozumowania. Jesli problem jest wystarczająąco złożony, próbowanie jest lepszym sposobem niż badanie wszystkich łańcuchów możliwości. Stanisław M. Ulam, Przygody matematyka. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Stanisław Ulam Pomysł polegał na wypróbowaniu tysięcy takich możliwosci z przypadkowym wybieraniem zdarzenia określającego los neutronu na każdym etapie procesu, przy użyciu liczb losowych (...) Po zbadaniu możliwych przebiegów procesu jedynie w kilku tysiącach przypadków będziemy mieli dobrą próbkę i przybliżoną odpowiedz na pytanie. Wszystko czego potrzeba, to metoda tworzenia takich przykładowych przebiegów. Tak sie złożyło, że właśnie powstały maszyny liczące... Wydaje mi się że nazwa Monte Carlo bardzo przyczyniła się do jej popularyzacji. Metoda została tak nazwana z powodu roli przypadku: generowania liczb losowych, które decyduja o przebiegu gry. Stanisław M. Ulam, Przygody matematyka. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Stanisław Ulam Cechą metody Monte Carlo jest to, że nigdy nie daje ona dokładnej odpowiedzi; wnioski z niej pokazuja raczej, że odpowiedz jest zawarta w pewnym przedziale błędu z takim a takim prawdopodobieństwem. ...zadanie polega na obliczeniu objętości obszaru określonego za pomoca pewnej liczby równan lub nierówności w wielowymiarowej przestrzeni. Zamiast korzystać z klasycznej metody polegającej na przybliżeniu tego obszaru siecią punktów czy komórek, co wymaga użycia miliardów pojedynczych elementów, można po prostu wybrać losowo kilka tysięcy punktów i za pomocą takiej próbki wyrobić sobie pojęcie o poszukiwanej wartości. Stanisław M. Ulam, Przygody matematyka. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Całkowanie metodą Monte Carlo Wybierzmy N przypadkowych punktów o rozkładzie jednorodnym w wielowymiarowym obszarze V (próbkowanie proste, random sampling), wtedy:
2 f - f 2 fdV H" V f ą V , N N N
1 1 2 2 f = f (xi), f = f (xi). N N i=1 i=1 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Całkowanie metodą Monte Carlo Odchylenie standardowe:
N E((X - E (X ))2) (f (xi) - f )2 i=1 = = = N N
N f 2 N f (xi ) N f 2 (xi ) - 2 f + i=1 i=1 i=1 N N N = N
N f 2 (xi ) 2 - 2 f 2 + f 2 f - f 2 i=1 N = = . N N Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Całkowanie metodą Monte Carlo " Dokładność rośnie jak N. Dla funkcji stałej - wynik dokładny. Redukcja wariancji - zamiana zmiennej tak, aby f była w przybliżeniu stała i próbkowanie nie za dużego obszaru. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Całkowanie metodą Monte Carlo Próbkowanie ważone (importance sampling) - jeżeli na dużym obszarze funkcja jest bliska zeru (jak w przypadku funkcji wykładniczej), to próbkowanie proste w wiekszości punktów będzie dawało znikomy wkład do wartości całki. Próbkując według rozkładu niejednorodnego w(x) można całkę zapisać następująco:
1 1 f (x) I = f (x)dx = w(x)dx. w(x) 0 0 Zakładając w(x) = u (x), gdzie u jest funkcją nieujemną i niemalejacą oraz u(0) = 0 i u(1) = 1, to
1 f [x(u)] I = du. w[x(u)] 0 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Całkowanie metodą Monte Carlo Po wygenerowaniu zbioru {ui}N według rozkładu jednorodnego: i=1 N
1 f [x(ui)] I H" , N w[x(ui)] i=1
(f /w)2 - f /w 2 , N zatem aby zredukowac błąd należy wybrać taką w(x) aby f /w była możliwie stała. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Szersze spojrzenie - symulacje molekularne Eksperyment komputerowy (computer experiment) - symulacja komputerowa odgrywa role eksperymentu zaprojektowanego do testowania teorii; można porównać wyniki symulacji z wynikami teorii - jeżeli wystepuje niezgodność to teoria jest zła. Uwaga! Symulacja komputerowa daje wgląd w zjawisko ale nie daje zrozumienia (tylko liczby). W niektórych dziedzinach są dobre teorie i nie potrzeba symulacji ale np. w teorii gęstych płynów jest bardzo mało danych eksperymentalnych. Często przewiduje się własności materiałów, np. łatwiej jest zmierzyć eksperymentalnie punkt zamarzania wody w warunkach pokojowych ale pod dużym ciśnieniem i temperaturą już nie. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Szersze spojrzenie - symulacje molekularne Przybliżenie Oppenheimera - równanie Schrodingera rodzielamy na zagadnienie ruchu elektronów w polu nieruchomych jąder (teoria struktury elektronowej w chemii kwantowej) i zagadnienie ruchu jąder w potencjale wyznaczonym przez energię elektronową. Pole silowe - przybliżenie potencjału energii elektronowej (kształtu hiperpowierzchni). Modelowanie molekularne Mechanika molekularna Dynamika molekularna Dynamika Monte Carlo Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Modelowanie molekularne Modelowanie molekularne - ogólny termin określający teoretyczne i obliczeniowe metody modelowania struktury i własności molekuł; metodami tymi bada się małe układy chemiczne oraz duże molekuly (białka) a także materiały (polimery itp.). Bada się: strukturę, dynamikę i własności termodynamiczne. Związki aktywne biologicznie - zwijanie się białek i ich stabilność, katalizę enzymatyczną, zmiany konformacji zwiazane z aktywnością biologiczną, rozpoznawanie molekularne, DNA, kompleksy membranowe. Najmniejszymi częściami układu są atomy, nie bierze się pod uwagę ich struktury (w przeciwieństwie do chemii kwantowej, która rozważa elektrony). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Mechanika molekularna Mechanika molekularna (MD) - metoda przewidywania stabilnej konfiguracji lub konformacji poprzez minimalizację V (R), czyli elektronowej energii molekuły lub układu molekuł. Lokalna MM - polega na szukaniu minimum lokalnego hiperpowierzchni energii potencjalnej np. metoda gradientowa (grad(V ) = 0; hessian); układ nie może poruszać się pod gorę potencjału. Giętka MM - wiązania traktowane są jak nierozrywalne sprężyny. Sztywna MM - długości wiązań i kąty pomiędzy nimi są ustalone (wartości doświadczalne); mniej zmiennych. MM z rozrywalnymi wiązaniami - oscylatory Morse a (-: obliczenie eksponenty jest droższe; 3 parametry a nie dwa; +: zwykle wiązania się nie rozrywają). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Mechanika molekularna Konformacje cykloheksanu - krzesłowa (1) i łodziowa (4): Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Globalna mechanika molekularna Globalna MM - niezależna od punktu początkowego; powinna odwiedzić każdy obszar hiperpowierzchni (problem ergodyczności). Schepens - liczba znalezionych konformerów jest proporcjonalna do czasu prowadzenia badań. Przykład: liczba minimów rośnie wykładniczo z ilością atomów w klastrze Lennarda-Jonnesa. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Globalna mechanika molekularna Analiza konformacyjna - znajdowanie konformacji w danych warunkach; w temperaturze <" 300K przeważają konformacje o niższej energii a zwykle najważniejsza jest konformacja odpowiadajaca minimum globalnemu. Uwaga! Rzeczywista konformacja nie zawsze jest minimum globalnym funkcji potencjału- lepszym kryterium jest minimum enegrii swobodnej E - TS (uwzględnia entropię a więc ilość dostępnych stanów przy danej energii) - w szerokiej studni potencjału jest więcej stanów oscylacyjnych niż w wąskiej
1 k ("E = h, = ). 2Ą Czasem molekuła może wskutek przeszkód sterycznych lub czynnika zewnętrznego zostać złapana w lokalne minimum (tzw. kinetyczne np. diament, fullereny), które jest różne od globalnego (termodynamicznego). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika molekularna Dynamika molekularna (MD) - bada ewolucję układu w czasie
przez rozwiązywanie równania Newtona mR = -"V (R); przez prędkości początkowe atomów uwzględnia się temperaturę. Trajektoria układu - zależność położeń atomów od czasu. Dynamika kwantowo-klasyczna - do opisu reakcji chemicznych np. układ enzym + rozpuszczalnik można podzielic na 3 regiony: I (centrum aktywne; równanie Schrdingera zależne od czasu), II (reszta enzymu; klasyczna MD), III (rozpuszczalnik; ośrodek ciągły o pewnej przenikalności dielektrycznej). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Krótki kurs strukturalnej chemii organicznej Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Krótki kurs strukturalnej chemii organicznej Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Krótki kurs strukturalnej chemii organicznej Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Krótki kurs strukturalnej chemii organicznej Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Krótki kurs strukturalnej chemii organicznej Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Krótki kurs strukturalnej chemii organicznej Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Krótki kurs strukturalnej chemii organicznej Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Krótki kurs strukturalnej chemii organicznej Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wykład 2. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Pola siłowe Pole siłowe (force field) - wyrażenie opisujące energię molekuły jako funkcję położenia jej jąder; jest to pole skalarne V (R). F = -"V - stąd nazwa; dim(R) = 3N W połowie XXw. ustalono podstawowe zależności geometryczne w molekułach ! wyjaśnianie budowy skomplikowanych molekuł i przewidywanie geometrii nowych. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Pola siłowe Pole siłowe metodami chemii kwantowej - należy rozwiązać równanie Schrdingera dla każdego położenia jąder Ri znajdując V (R) = {V (Ri)}. Przykład: Jeżeli N = 100, dla każdego wymiaru 100 punktów (wartości potencjału), założmy że obliczenie każdego punktu trwa 1 minutę, to czas obliczenia pola siłowego wynosi: 1003"1001min = 10600min. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Pola siłowe Wiązania - każde wiązanie ma długość r0 dla której jego wkład do energii będzie najmniejszy, w pobliżu minumum 1 można przyjąć kXY (r - r0) - wyrażeniu temu przypisuje się 2 uniwersalną ważność dla wiazań X - Y danego typu niezależnie od otoczenia tego wiązania (tj. to samo r0 i kXY ). Poprawki anharmoniczne - można uwzględnić zmniejszającą się siłę potrzebną do rozciągania wiązania. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Pola siłowe Kąty - mając dane optymalne długości A - B i B - C i zmieniając kąt A - B - C zmieniamy też odległość AC - atomy te nie są związane wiązaniem chemicznym, więc zmiana energii będzie mniejsza niż w przypadku zmiany długości wiązania. Istnieją optymalne kąty między danymi wiązaniami - odchylenia od nich także uwzględniamy członem 1 harmonicznym: kXYZ (ą - ą0)2 2 Przykład: kC-C -C = 0.0099molkcal podczas gdy stopien kcal kC -C = 317 . mol Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Pola siłowe Oddziaływania van der Waalsa - atomy oddalone od siebie bardziej niż 2 wiązania także oddziaływują. Promień van der Waalsa - jeżeli rF = odległość równowagi przy zbliżaniu do siebie dwóch atomów F i analogicznie rCl, to gdy zbliżamy do siebie H - F i Cl - H, to odległość równowagowa H" rF + rCl. Atomy zwykle nie zbliżają się na odległość większą niż suma ich promieni van der Waalsa (enegria gwałtownie rośnie) a odalane będą się przyciągać proporcjonalnie do r-6 ! jeżeli w odległości równowagowej energia wynosi -, to oddziaływanie opisuje wzór Lennarda-Jonesa (LJ 12-6): 12 6 r0 r0 VLJ(X , Y ) = - 2 . r r Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Pola siłowe Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Pola siłowe Oddziaływania elektrostatyczne - q1q2/r. Kąty torsyjne - w przypadku związanych atomów X - Y - Z - W rotacja Y - Z o kąt wprowadza dodatkową energię AXYZW (1 - cos n), gdzie n - krotność wystąpienia bariery przy pełnym obrocie np. etan n=3 (vs. etylen, n=2). Człony mieszane (sprzegające) np. wiązanie - kąt (stała siłowa zmiany kąta powinna zależeć od aktualnych długości wiązan). Człony uwzględniające przenikalność dielektryczną ośrodka - gdy uwzględniamy ośrodek bez jego struktury molekularnej, np. polaryzacja ośrodka przez molekułe. Inne np. wiązania wodorowe, polaryzowalność (na ładunek atomu wpływają oddziaływania z innymi atomami). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Pola siłowe Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Pola siłowe
1 1 V = kXY (r - r0)2 + kXYZ (ą - ą0)2 + VLJ(X , Y )+ 2 2 X -Y X -Y -Z XY
q1q2 + + AXYZW (1 - cos n) r tors XY Dodatkowa zaleta - pole siłowe daje funkcję V (R) w postaci wzoru dla dowolnego R. Zwykle dla tego samego zagadnienia jest wiele pól siłowych - często rezultaty bardzo podobne ale nie zawsze (nawet z tej samej geometrii początkowej). Parametrów z jednego pola siłowego nie powinno się używac w innym, gdyż są one optymalizowane dla danej postaci funkcyjnej, która w innym polu może być inna. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Pola siłowe Funkcje i parametry w polach siłowych mogą pochodzić z eksperymentów (spektroskopia, rtg, entalpia parowania) lub obliczeń kwantowomechanicznych. Pola siłowe zjednoczonego atomu (united-atom FF) - parametry nie są dane dla każdego atomu ale grupy typu CH3 traktuje się jako pojeyncze centrum oddziaływania. Gruboziarniste pola siłowe (coarse-grained FF) - nawet większe części molekuły traktują jako pojedynczy punkt np. w symulacji białek. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Problemy z polami siłowymi Siły van der Waalsa są silnie zależne od środowiska (wpływ rozpuszczalnika) - McLachlan (1963) - przyciaganie w rozpuszczalniku jest słabsze i zakłada że różne atomy przyciągają się słabiej niż identyczne; teoria ta dopuszcza oddziaływania czysto odpychające jak w ciekłym helu. Klasyczne pola siłowe opierają się na regułach kombinatorycznych (combinatorial rules) - energia oddziaływania różnych atomów jest średnią energii oddziaływań takich samych atomów np. E(C...N) = [E(C...C) + E (N...N)]/2. Struktura białek - uczestnicy CASP (Critical Assesment of Techniques for Protein Structure Prediction) nie próbują uniknąć tego, że minimalizacja energii lub MD prowadzi do struktury różnej od wyznaczonej doświadczalnie - zamiast tego rozwinięto inne wielkości do oceny przebiegu symulacji. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Popularne pola siłowe MM2 - stworzone do analizy konformacyjnej małych związków organicznych; posiada duży zbiór parametrów, który jest ciągle ulepszany i dostosowywany do różnych klas związków organicznych (MM3, MM4). CFF - stworzone do ujednolicenia obliczeń energii, struktury i wibracji czasteczek i krysztalow molekularnych; program CFF posłużył do rozwinięcia wielu innych. ECEPP - do modelowania peptydow i protein; dla uproszczenia używa ustalonych geometrii aminokwasów, więc zmiennymi są tylko kąty torsyjne. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Popularne pola siłowe AMBER, CHARMM, GROMOS - rozwinięte do MD makromolekuł, ale są powszechnie używane do minimalizacji energii - współrzędne wszystkich atomów są zmiennymi. EVB (Empirical Valence Bond) - reaktywne pole siłowe (reactive force field) - prawdopodobnie najlepsze do modelowania reakcji chemicznych w różnych środowiskach; pozwala obliczać energie swobodne aktywacji. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Modele wody Model wody - pole siłowe dla wody; ważne ze względu na specjalne własności wody i jej ważność jako rozpuszczalnika. Klasyfikujemy je ze względu na liczbę centrów definiujących model, sztywność - giętkość, efekty polaryzacyjne. Najprostsze modele są sztywne np. TIP3P. na a na b
kqiqj A B Eab = + - , 12 rij rOO B r00 i j k=332.1kcal/mol, q - ładunki cząstkowe, które mogą być rozmieszczone na atomach lub innych centrach (wolne pary elektronów). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Modele wody Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Modele wody TIP3P (model 3-centrowy) : rOH, (HOH), A, B, q(O) < 0, q(H) > 0. TIP5P (model 5-centrowy) : rOH, (HOH), r(OL), (LOL), A, B, q(L), q(H). TIP5P przewiduje takie własności jak geometria dimeru wody, radialny rozkład dyfrakcji neutronów, temperaturę maksymalnej gęstości wody. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Modele wody Model SPC/E - zawiera poprawkę polaryzacyjną
1 ( - 0)2 Epol = = 1.25kcal/mol, 2 ąi i - moment dipolowy spolaryzowanej cząsteczki H2O, 0 - moment dipolowy cząsteczki izolowanej, ąi - izotropowa stała polaryzowalności. Modele gruboziarniste - jedno- lub dwucentrowe. Czas obliczeń - proporcjonalny do ilości centrów w modelu np. dla 3-centrowego modelu jest 9 odległości dla każdej pary molekuł wody a dla 5-punktowego 4 4qq + 1OO = 17. Modele sztywne - wymagają dodatkowego kosztu na utrzymanie założonej struktury (istnieją odpowiednie algorytmy) ale za to można często zwiekszyć krok czasowy (time step). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo UFF - Dokładniejsze omówienie wybranego pola siłowego UFF (Universal Force Field) - zaprojektowane do opisu wszystkich atomow, więc obejmuje i zwiazki nieorganiczne. Parametry bazują na trzech zmiennych: rodzaju atomu, hybrydyzacji i wiązaniach. Typ atomu - 126 typów opisanych symbolem pierwiastka, geometrią/hybrydyzacją i jednym dodatkowym parametrem; 1 = linia, 2 = trójkat, R = aromatyczny, 3 = tetraedr, 4 = kwadrat, 5 = bipiramida trygonalna, 6 = oktaedr; np. N 3, Rh6; dodatkowy parametr - P 3 q fosfor tworzący wiazanie koordynacyjne np. (Ph3P)2PtCl2. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo UFF Postać funkcyjna E = ER + E + EĆ + E + EvdW + Eel - ER - wiązania; oscylatory harmoniczne lub Morse a; długość wiązania składa się z promieni atomów tworzących wiązanie (promienie brane są z ustalonych związków np. CR z benzenu), poprawki rzędu wiązania rBO = -0.1332(rI + rJ) ln (n) i poprawki elektroujemności " " rEN = rI rJ( I + J)2/(I rI + JrJ); " 2 ZI"ZJ stałe siłowe ER = E0 - FrIJ - G ! kIJ = (" ER )0 = rIJ "R2 " " 664.12ZI ZJ [kcal/(mol )], Z" - efektywne ładunki (w 3 rIJ jednostkach atomowych) dopasowane do pewnego zbioru danych. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo UFF - E - kąty między wiązaniami. Dla kątów (ogolnie) - rozwiniecie Fouriera m
Eł = K Cn cos (nł) [kcal/(mol mol2)] n=0 m
E = KIJK Cn cos (n), n=0 Cn - spełniają warunki brzegowe np. aby funkcja miała minimum przy danym 0. " 2 ZI"ZK E = E0 - F - ! KIJK = (" ER )0 rIK "R2 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo UFF - EĆ - kąty torsyjne. m
EĆ = KIJKL Cn cos (nĆIJKL), n=0 KIJKL, Cn - okreśłone są przez barierę rotacji VĆ, periodyczność potencjału i kąt równowagowy. - E - człony inwersyjne Dla atomu I związanego z trzema innymi atomami J, K , L używa się jedno- lub dwuczłonowego rozwinięcia Fouriera: E = KIJKL(C0 + C1 cos IJKL + C2 cos 2IJKL) [kcal/mol], IJKL = (oś IL, płaszczyzna IJK ); brane są pod uwagę wszystkie osie (IL, IJ, IK). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo UFF cos 2 : min = 90ć%, max = 0ć% (PH3) cos : min = 0ć%, max = 180ć% (etylen) Kombinacja liniowa powyższych przpadków pozwala opisać przypadki pośrednie. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Inwersja Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Inwersja Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo UFF - EvdW - oddziaływania van der Waalsa LJ 6-12: 12 6 xIJ xIJ EvdW = DIJ - 2 , r r DIJ [kcal/mol], xIJ [] - otrzymywane z reguł kombinatorycznych (xIJ - średnia arytmetyczna, dla kryształów - geometryczna; DIJ - geometryczna). - Eel - oddziaływania elektrostatyczne - Nie ma oddziaływań van der Waalsa i elektrostatycznych 1, 2 i 1.3. - Dokładność - błąd poniżej 0.1 i 5 - 10ć%. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wykład 3 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Metropolis Monte Carlo Metropolis Monte Carlo generuje proces Markowa konfiguracji/stanów próbnych o rozkładzie prawdopodobieństwa danym przez rozkład normalny. Ą(o n) - prawdopodobieństwo przejścia ze starej konfiguracji o do nowej n (macierz przejścia, transition matrix); macierz Ą musi spełniać warunek: po osiagnięciu rozkładu równowagowego konfiguracji układ musi w nim pozostać (można wyobrazić sobie symulację wielu kopii układu zachodzących jednocześnie). Oznacza to, że w stanie równowagi prawdopodobieństwo opuszczenia przez układ konfiguracji o musi być równe sumie prawdopodobieństw przejscia układu ze wszyskich innych konfiguracji n do o. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Metropolis Monte Carlo Przyjmuje się silniejszy warunek, że rowność taka musi zachodzić pomiędzy dwoma dowolnymi konfiguracjami układu o oraz n (równowaga szczegółowa, detailed balance): N (o)Ą(o n) = N (n)Ą(n o), gdzie N jest rozkładem prawdopodobieństwa znajdowania się układu w danej konfiguracji. Wybór macierzy przejścia spełniającej powyższy warunek nie jest jednoznaczny. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Metropolis Monte Carlo Macierz przejścia Ą można przedstawić w poniższej formie: Ą(o n) = ą(o n)acc(o n), gdzie ą(o n) określa prawdopodobieństwo danego ruchu próbnego (trial move) a acc(o n) jest prawdopodobieństwem zaakceptowania danego ruchu. Najczęściej wybieramy ą symetryczną (ą(o n) = ą(n o)): N (o)acc(o n) = N (n)acc(n o), acc(o n) N (n) = = exp{-[U(n) - U(o)]}. acc(n o) N (o) Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Metropolis Monte Carlo Ponownie, wybór acc nie jest jednoznaczny, Metropolis przyjął następującą postać (i z dotychczasowej praktyki wydaje się że jest to najlepszy wybór - daje bardzo efektywne próbkowanie przestrzeni konfiguracji układu):
N (n)/N (o), N (n) < N (o) acc(o n) = 1, N (n) N (o). Prawdopodobieństwo przejścia jest zatem dane wzorami: Ą(o n) = ą(o n) N (n) N (o) = ą(o n)[N (n)/N (o)] N (n) < N (o)
Ą(o o) = 1 - Ą(o n). n =o Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Metropolis Monte Carlo Decyzja czy ruch próbny zostanie zaakceptowany czy odrzucony: jeżeli wygenerowany został ruch próbny o n taki, że U(n) > U(o) to akceptujemy go z prawdopodobieństwem równym acc(o n) = exp{-[U(n) - U(o)]} < 1. Następnie należy wygenerować liczbę losową o rozkładzie jednorodnym rand " [0, 1]. Prawdopodobieństwo, że rand < acc(o n) wynosi acc(o n), zatem dany ruch wykonujemy, jeżeli rand < acc(o n), w przeciwnym wypadku jest on odrzucany. Ważne jest żeby Ą(o n) bylo ergodyczne, czyli pozwalało osiągnąć każdy punkt przestrzeni konfiguracyjnej z dowolnego innego punktu w skończonej liczbie kroków MC. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Algorytm metody Monte Carlo 1 Wybierz przypadkowo cząstkę (np. atom) i oblicz jej energię U(rN). 2 Zmnień jej położenie o przypadkową wartość r = r + " i oblicz nową energię U(r N) 3 Zaakceptuj wykonany ruch rN r N z prawdopodobieństwem acc(o n) = min(1, exp{-[U(n) - U(o)]}). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Liczby losowe z przedziału (0, 1) można odwzorować na dowolny przedział za pomocą funkcji liniowej "x"(0,1) y = a + (b - a)x ! y " (a, b) Liczby te zwykle uzyskujemy za pomocą algorytmów rekeurencyjnych - są to więc liczby pseudolosowe. Podstawową operacją wykorzystywaną w prostych generatorach jest modulo y = x mod m ! "n x = nm + y. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Dobre własności statystyczne - rozkład prawdopodobienstwa generowanych liczb powinien być możliwie bliski założonego. Test jednorodności rozkładu: Niech dana będzie d-wymiarowa regularna sieć o rozmiarach L L, gdzie L jest liczbą węzłów w danym kierunku. Jeżeli zapełniać taką sieć wybierając przypadkowo węzły (liczby pseudolosowe z przedziału (0, 1) można odwzorować na liczby całkowite wzorem [rL] + 1), to dobry generator powinien zapewnić możliwość osiągnięcia wszystkich punktów sieci. Jeżeli generujemy NLd liczb pseudolosowych, to średnio każdy węzeł sieci powinien być wylosowany N razy. Natomiast ułamek pustych węzłów sieci powinien maleć proporcjonalnie do exp(-N). Należy stosować możliwie duże sieci; zwykle L " (20, 100) oraz N " (10, 20). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Periodyczność generowanego ciągu liczb. Generator Lehmera: Ii+1 = (aIi + c)mod m, aby z liczb całkowitych otrzymać liczby rzeczywiste z przedziału (0, 1) oblicza się iloraz ri = Ii/m, gdzie m jest zwykle największą możliwą liczbą całkowitą dostepną dla danego procesora. Okres ciągu tak generowanych liczb jest nie większy niż m i zależy od a, c i I0. Przykład: m = 16, a = 3, c = 1, I0 = 2 ! 2, 7, 6, 3, 10, 15, 14, 11, 2, 7, ... Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Najlepiej, gdy m jest liczbą pierwszą np. m = 17 2, 7, 5, 16, 15, 12, 3, 10, 14, 9, 11, 0, 1, 4, 13, 6, 2, 7, ... Okres generatora musi być znacznie większy niż ilość liczb, które będą potrzebne do symulacji; zwykle powinno być to ponad 1010. Istotne znaczenie ma dobór parametrów. Jeżeli słowo używane przez procesor jest n-bitowe, to największa liczba całkowita to (jeden bit określa znak liczby) mmax = 2n-1 - 1. Przykład: dla procesora 32 bitowego jest to 231 - 1 = 2147483647 <" 109. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Maksymalny okres występuje dla c = 0 ale wtedy wlasności
statystyczne są zle, zatem przyjmuje się c = 0 i odpowiednio dobiera się a. Korelacja między kolejnymi liczbami generowanego ciągu jest proporcjonalna do odwrotności a ale z drugiej strony powinno " wybierać się a znacznie mniejsze od 2n-1. W praktyce zwykle stosuje się a = 16807, " a = 16807 < 232-1 H" 46341. Uwaga! Początkowa liczba I0 (ziarno generatora, seed) powinna być duża; w przeciwnym razie generator należy rozgrzac generując krótki ciąg liczb, które muszą być odrzucone przed rozpoczęciem właściwej symulacji. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Efekt Marsaglii - rozkłady generatorów liniowych (opartych o modulo) układają się na regularnych hiperpłaszczyznach wewnątrz kostki [0, 1]k. (U1, U2, ..., Uk), (U2, U3, ..., Uk+1), ... (U1, U2, ..., Uk), (Uk+1, Uk+2, ..., U2k), ... Ogólnie: Xn = (a1Xn-1 + a2Xn-2 + ... + akXn-k + c) modm, okres maksymalny - mk - 1. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Testy G. Marsaglii DIEHARD - pomogły wyeliminować wiele złych generatorów http://www.stat.fsu.edu/pub/diehard Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Generator Tauswortha - duża wydajność, okresy rzedu 1050 oraz bardzo dobra jednorodność rozkładu. Wykorzystuje operację XOR (exclusive OR; różne bity = 1, jednakowe = 0). Algorytm: 1 Generacja zbioru podstawowego M całkowitych liczb pseudolosowych przy pomocy generatora modulo. 2 Rozgrzanie generatora. 3 Generacja całkowitej liczby pseudolosowej Ik: Ik = XOR(Ik-q, Ik-p), Ik-q, Ik-p - liczby ze zbioru M; p i q - odpowiednio dobrane. 4 Obliczenie rzeczywistej liczby pseudolosowej Ik/Imax. 5 Modyfikacja zbioru podstawowego przez podstawienie Ik w miejsce jednego z dotychczasowych elementow. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Takie generatory nazywane są generatorami opartymi na rejestrach przesuwanych. bn = (a1bn-1 + ... + akbn-k) mod2, a1, ..., ak " {0, 1} (a + b) mod2 = a xor b Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Generator RANLUX (M.Lscher, 1993) - spełnia wszystkie znane testy statystyczne! Pozwala uzyskac wyniki o duzej precyzji. Oparty jest o generator RCARRY uzupełniony o algorytm odrzucania pewnch sekwencji liczb - usuwane są krótkozasięgowe korelacje; użycie wykładników Lapunowa i entropii Kolmogorowa. Okres: 5.2 10171. Xn = Xn-s Xn-r modm (r, s " N, r > s; r = 24, s = 10, m = 224) Inicjacja: X1, ..., Xr " (0, m), c - nie mogą to być dowolne liczby. xn-s yn-r modm =
xn-s - yn-r - cn-1 + m, cn = 1, gdy x - y - cn-1 < 0 = xn-s - yn-r - cn-1, cn = 0, gdy x - y - cn-1 0. cn - bity niosące . Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Korelacje są krótkozasięgowe i usówa się je nastepująco: używa się r liczb, odrzuca się p - r następnych itd. p r - parametr kontrolujący ilość odrzucanych liczb. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Zufall - oparty na ciągu Fibonacciego, generuje liczby z zakresu [0, 1); t = u(n - 273) + u(n - 607), u(n) = t - int(t). Fortran - call random seed; można generowac kilka liczb na raz. C - srand() - ziarno; rand() " (0, RandMax) )" N, gdzie RandMax > 32767. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Metoda Boxa-Mullera - generuje liczby pseudolosowe o rozkładzie normalnym: 1 2 " x <" N(0, 1) = e-x /2. 2Ą Zdefiniujmy funkcję
1 2 U(R) a" dxdy e-(x +y2)/2, 2Ą x2+y2 R2
2Ą R R2/2 1 2 2 U(R) = d dr re-r /2 = du e-u = 1 - e-R /2. 2Ą 0 0 0 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych lim U(R) = 0, lim U(R) = 1. R0 R " Dla liczby p " [0, 1]:
U(R) = p ! R = -2 log(1 - p). Po wprowadzeniu zmiennych s = 1 - p " [0, 1] i t " [0, 1] ponizsze x posiada rozkład normalny:
-2 log(s) cos(2Ąt) x = -2 log(s) sin(2Ąt). Uogolniając zmienna losowa o rozkładzie x <" N(, 2) dana jest wyrażeniami:
+ -2 log(s) cos(2Ąt) x = + -2 log(s) sin(2Ąt). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Generowanie liczb pseudolosowych Metoda Boxa-Mullera jest przypadkiem ogólnej metody odwracania dystrybuanty: Jeżeli U jest liczbą losową z rozkładu równomiernego na [0, 1], a F dystrybuanta (funkcja sciśle rosnąca, F (- ") = 0, F ( ") = 0), to zmienna losowa -1 X = F (U), ma rozkład o dystrybuancie F . Zalety: dokładna; prosta; szybka dla niektórych rozkładów. Wady: dystrybuanta powinna byc znana analitycznie i odwracalna; gdy całkujemy numerycznie - wolniej, mniej dokładnie, niestabilnosci numeryczne; trudna do zastosowanie dla rozkładow wielowymiarowych. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wykład 4 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Szukanie minimum globalnego - metoda przeskakiwania basenów oddziaływania Globalna optymalizacja: Problem komiwojażera. Znajdywanie optymalnej scieżki przesyłu danych w procesorach i internecie. Przewidywanie aktywnej struktury biomolekul. Klastery Lennarda-Jonesa (147-atomowy klaster posiada 1060!). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Szukanie minimum globalnego - metoda przeskakiwania basenów oddziaływania Alternatywne metody poszukiwania minimum globalnego: Algorytmy genetyczne (genetic algorithms; naśladowanie ewolucji genetycznej - mutacje). Metody deformacji hiperpowierzchni (hypersurface deformation methods; transformowanie powierzchni energii potencjalnej tak, aby została ona wygładzona; po minimalizacji - transformacja odwrotna). Symulowane schładzanie (simulated annealing; w wyższej temperaturze powierzchnia energii swobodnej jest prostsza) Dyfuzyjna metoda MC (wykorzystuje tunelowanie kwantowe do pokonania barier potencjału - znajduje się funkcję falową stanu podstawowego, która powinna lokalizować się z 0; metoda skaluje się jak 2N). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Szukanie minimum globalnego - metoda przeskakiwania basenów oddziaływania Metoda przeskakiwania basenów oddziaływania (basin-hopping): Następująca transformacja nie zmienia minimum globalnego ani minimow lokalnych: Ę(X) = min E (X), X " R3N, min - oznacza minimalizację lokalną zaczynając od X (np. metoda sprzeżonych gradientów). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Szukanie minimum globalnego - metoda przeskakiwania basenów oddziaływania Do funkcji Ę(X) stosujemy metodę Monte Carlo: Minimalizacja (kryterium zbieżności - średnia kwadratowa (root-mean-square) gradientu i "Ei są mniejsze niż pewna wartość). Każda współrzędną przemieszczamy o liczbę z przedziału [-1, 1] " step (step dobiera się tak, żeby ok. 0.5 ruchów było akceptowane; dzięki transformacji krok może być większy niż bez niej). Można przeprowadzić kilka symulacji zaczynając od różnych (przypadkowych) wartości początkowych. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Szukanie minimum globalnego - metoda przeskakiwania basenów oddziaływania Metoda ta pozwala zidentyfikowac wszystkie znane minima do LJ110. Wszystkie najtrudniejsze przypadki zostały zlokalizowane przez ogólny ( niestronniczy ) algorytm. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Mechanika statystyczna Mechanika statystyczna - pośredniczy między wielkościami mikroskopowymi a makroskopowymi; nie wszyskie wielkości fizyczne mogą być zmierzone w symulacji (odwrotnie: także wielkości mierzone w symulacji zwykle nie odpowiadają wielkościom mierzonym eksperymentalnie). Przykład: MC i MD pozwalają mierzyć chwilowe położenia atomów jednak doświadczenie nie daje takich informacji ale wielkości uśrednione zarówno w po atomach (molekułach) i czasie pomiaru. Większość symulacji zakłada, że dynamikę atomów i molekuł opisuje mechanika klasyczna - obliczenia są proste i najczęściej przybliżenie to jest usprawiedliwione. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Mechanika statystyczna &!(E, V , N) - liczba stanów układu N cząstek o energii E w objetości V ; układ opisany wielkościami (E, V , N) można znalezć z jednakowym prawdopodobieństwem w każdym z tych stanów. Równowaga termodynamiczna - ln &! ma wartość maksymalną.
"S Temperatura: 1/T = , S(N, V , E ) = kB ln &!(N, V , E). "E V ,N
" ln &!(E ,V ,N) = 1/(kBT ) = "E N,V Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Mechanika statystyczna Rozkład Boltzmanna: exp(-Ei/kBT ) Pi = exp(-Ej/kBT ) j Energia średnia (wartość oczekiwana energii):
1 "Q " ln Q E = EiPi = - = - , Q "(1/kBT ) "(1/kBT ) i Q - funkcja podziału. Energia swobodna Helmholtza:
F = -kBT ln Q = -kBT ln exp(-Ei/kBT ) . i Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Mechanika statystyczna Wartość średnia (oczekiwana) wiełkości A układu w równowadze termodynamicznej (thermal average):
exp(-Ei/kBT ) i|A|i i A = , exp(-Ej/kBT ) j równania tego nie można użyć, gdyż trzeba by wyznaczyć wszystkie stany wielocząstkowego układu kwantowego. Jednak korzystając z relacji exp(-Ei/kBT ) = i| exp(-H/kBT )|i ,
i| exp(-H/kBT )A|i Tr[exp(-H/kBT )A] i A = = . j| exp(-H/kBT )|j Tr exp(-H/kBT ) j Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Mechanika statystyczna exp(-K) exp(-U) = exp{(-[K + U + O([K, U]))}, [K, O] <" , Tr exp(-H) H" Tr[exp(-K) exp(-U)]. Jeżeli |k i |r są odpowiednio wektorami własnymi operatora pędu i polożenia:
Tr exp(-H) = r| exp(-U)|r r|k k| exp(-K)|k k|r . r,k Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Mechanika statystyczna r| exp(-U)|r = exp[-U(rN)],
N
2 k| exp(-K)|k = exp - pi /(2mi) , pi = ki, i=1 N r|k k|r = 1/V ,
1 2 Tr exp(-H) H" dpNdrN exp - pi /(2mi) + U(rN) hdNN! i a" Qclassical, hdN - objetość w przestrzeni fazowej, 1/N! - dla uwzględnienia nierozróżnialności identycznych cząstek. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Mechanika statystyczna
A = 2 dpNdrN exp{-[ pi /(2mi) + U(rN)]} i Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Ergodyczność Średnia po zespole statystycznym - średnia po wszystkich możliwych stanach kwantowych; zwykle tak nie myśli się o wartościach średnich ale oblicza się średnie z serii pomiarów w pewnym przedziale czasowym. Jeżeli mierzymy pewną wielkość x(r) (np. i(r)), to zakladając że dla wystarczająco długiego czasu symulacji średnia po czasie nie zależy od warunków początkowych (uwaga! to nie zawsze jest prawda) to średnia po czasie dana jest wzorem:
t 1 x(r) = lim x(r; t )dt . t " t 0 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Ergodyczność
limt " 1 t x(r; t )dt warunki poczatkowe t 0 x(r) = = warunki poczatkowe
t x(r; t ) 1 warunki poczatkowe = lim dt t " t warunki poczatkowe 0
t 1 x(r; rN(0), pN(0), t )drNdpN E lim dt t " t &!(N, V , E) 0
t 1 = lim x(r; rN(0), pN(0), t ) NVE dt . t " t 0 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Ergodyczność Uśrednianie po czasie można opuścić, ponieważ średnia po zespole statystycznym (NVE) nie zależy od czasu; uśrednianie po warunkach początkowych jest równoważne uśrednianiu po współrzędnych przestrzennych wyewoulowanych w czasie; zatem: x(r) = x(r) NVE . Powyższe równanie nazywamy hipotezą ergodyczną ; nie jest spelniona np. w układach metastabilnych (szkła) czy prawie harmoniczne ciała stałe. L = podejście MD; P = podejście MC Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo Zajmiemy się symulacjami w zespole statystycznym (N, V , E ).
dpNdrNA(pN, rN) exp{-H(pN, rN)}
A = . dpNdrN exp{-H(pN, rN)} Całki tego rodzaju nie można obliczyć zwykłymi metodami: jeżeli system jest DN-wymiarowy, to przy m punktach dla każdego wymiaru liczba punktów, w których trzeba obliczyć funkcję podcalkową wynosi mDN; przyjmując DN = 300 i m = 5, to 5300 = 51.43"210 = 10210. Dodatkowy problem - exp jest szybko zmnieniającą się funkcją która jest bliska zero na dużym obszarze, niestety nie można zastosować techniki próbkowania ważonego, gdyż nie znamy odpowiedniego rozkładu.
drNA(pN, rN) exp{-U(rN)}
A = . drN exp{-U(rN)} Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo Oznaczając funkcje rozkładu przez Z , gęstość prawdopodobieństwa znalezienia układu w sąsiedztwie rN wynosi exp[-U(rN)] N (rN) a" 0. Z Jeżeli generować punkty w przestrzeni według tego rozkładu, to średnio liczba punktów przypadających na otoczenie punktu rN jest równa ni = LN (rN), gdzie L jest całkowita liczba wygenerowanych punktów. N
1 N A H" niA(ri ). L i=1 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo Próbkowanie proste - można otrzymać całkowitą powierzchnię rzeki. Próbkowanie ważone - nie można jej otrzymać (co jest odpowiednikiem tego, że nie znamy Z ). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo Dlaczego ponownie liczyc stare konfiguracje? Prawdopodobienstwo przejścia musi być znormalizowane:
Ą(o n) = 1 ! Ą(o o) = 1 - Ą(o n), n n =o stąd wynika, że prawdopodobieństwo pozostania w starym stanie może być niezerowe. Przykład: Układ dwustanowy {E1, E2} - jeżeli nie będziemy liczyć starych konfiguracji to średnia zawsze będzie wynosić E = (E1 + E2)/2 (niezależnie od temperatury). Ponadto, błędy są większe. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wyklad 5 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - szczegóły techniczne Warunki brzegowe Symulacje MC mają dostarczać informacji o własnościach makroskopowej próbki. Większość symulacji bada strukturalne i termodynamiczne własności układów kilkuset-kilku tysięcy cząstek - nie jest to próbka makroskopowa. Dla tak małych układów nie można zakładać, że wybór warunków brzegowych nie wpływa na własności symulowanego układu. W trójwymiarowym układzie N cząstek ułamek molekuł znajdujących się na powierzchni układu jest proporcjonalna do N-1/3 np. dla sześciennego kryształu 1000 cząstek jest to 49% a dla 106 6%. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - szczegóły techniczne Periodyczne warunki brzegowe - układ wygląda jak nieskończony; objętość zawierająca N cząstek jest traktowana jak komórka nieskończonej sieci. Każda cząstka oddziaływuje z wszystkimi innymi w sieci (także swoimi własnymi obrazami). Całkowita energia N cząstek w dowolnej sześciennej komórce L L L:
1
Utot = u(|rij + nL|), 2 i,j,n gdzie n " Z3 i dla n = 0 wykluczane są wyrazy takie, że i = j. Pojawia się problem - powyższa suma jest nieskończona. Najczęściej oddziaływania są krótkozasięgowe a więc można dokonać obcięcia powyżej pewnego rc. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - szczegóły techniczne Periodyczne warunki brzegowe są zazwyczaj bardzo skuteczne ale należy uważać na możliwość dodatkowych korelacji, które mogą nie być obecne w rzeczywistym układzie. Fluktuacje mogą mieć długość co najwyżej = L. Przykładowo fluktuacje o dużej długości fali są ważne w pobliżu przejścia fazowego - w tym wypadku należy spodziewać się problemów przy stosowaniu warunków periodycznych. Uwaga! Wybranej komórce elementarnej nie należy przypisywać żadnego znaczenia - jej początek może być wybrany zupełnie dowolnie i nie zmieni to żadnej właśnosci układu. Ale kształt i orientacja jest ustalona. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - szczegóły techniczne Obcinanie oddziaływań Wybierając obcięcie rc dostatecznie duże błąd może być dowolnie mały. Przy periodycznych warunkach brzegowych, jeżeli rc < L/2 wystarczy uzwględnic oddziaływanie danej cząstki i tylko z najbliższymi obrazami każdej cząstki j. Potencjał nie jest dokładnie równy zero dla r rc, co powoduje błąd systematyczny. Jeżeli oddziaływanie maleje wystarczająco szybko (szybciej niż r-3; ważny przypadek oddziaływań kulombowskich nie spełnia tego warunku) można wprowadzić poprawkę:
"
N Utot = uc(rij) + dr u(r)4Ąr2, 2 rc igdzie uc - obcięta energia potencjalna, - średnia gęstość. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - szczegóły techniczne Mimo że wartość potencjału maleje z odległością miedzycząsteczkową r, to ilość cząstek biorących dających wkład do oddziaływania rośnie: liczba cząstek w odległości r od danego atomu rośnie asymptotycznie jak rd-1 (d - wymiar przestrzeni układu). Przykład - poprawka dla potencjału Lennarda-Jonesa (średnia energia na jeden atom):
" " 1 1 u = dr u(r)(r)4Ąr2 = 4Ą dr u(r)r2 2 2 rc 12 rc 6 " 1 = 16Ą dr r2 - 2 r r rc 9 3 8 = Ą3 - . 3 rc rc Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - szczegóły techniczne Proste obcięcie
u(r) r rc utrunc(r) = 0 r > rc Nieciągłość będzie dawać dodatkowy (niezaniedbywalny) wkład do ciśnienia; można to skorygować. Obcięcie i przesunięcie
u(r) - u(rc) r rc utr-sh(r) = 0 r > rc Nie ma nieciągłości; energia i ciśnienie takiego układu będzie inne niż w poprzednim przypadku i bez obcięcia ale można skorygowac ten efekt. Ten rodzaj obcięcia musi być używany ostrożnie do oddzialywań anizotropowych. Konwencja minimalnego obrazu - obliczane jest oddziaływanie z najbliższym obrazem każdej z cząstek. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - szczegóły techniczne Inicjalizacja - początkowe położenia cząstek mogą być dość dowolne (każde rozsądne jest akceptowalne) bo własności równowagowe układu nie powinny zależeć od warunków początkowych. Wiemy jak wygląda cząsteczka/i, które symulujemy. Ciało stałe - znamy strukturę krystaliczną. Uniezależnienie się od warunków początkowych daje tak zwane rozgrzanie układu. Jeżeli wyniki zależa od warunków początkowych, to: 1 Najprawdopodobniej nie został osiągnięty stan równowagi (próbkowanie przestrzeni konfiguracji może być nieprawidłowe). 2 Układ jest nieergodyczny. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - szczegóły techniczne Jednostki zredukowane - wybieramy dogodną jednostkę energii, długości i masy a wszystkie inne jednostki wyrażamy przez te wybrane (temperaturę, gęstość, ciśnienie). Przykład: uLJ(r) = f (r/), naturalny wybór to długość - , energia - oraz masa - m.
Wtedy np. jednostką temperatury jest /kB, czasu m/, ciśnienia - /3, gęstości - 1/3. Zredukowany potencjał u" a" u/ jest funkcją bezwymiarową zredukowanej odległości r" a" r/: 12 6 1 1 " uLJ(r") = 4 - . r" r" Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - szczegóły techniczne Korzyści z jednostek zredukowanych: Nieskończenie wiele kombinacji jednostek rzeczywistych odpowiada tym samym jednostkom zredukowanym np. Ar, T = 60K, = 840kg/m3 i Xe, T = 112K, = 1617kg/m3 " odpowiadają jednostkom zredukowanym " = 0.5, T = 0.5. Dla Ar parametry potencjału LJ: = 3.405 10-10m oraz = 1.6537 10-21J. kB 1.3806 10-23J/K " T = T = 60K = 0.5 1.6537 10-21J 840kg/m3 " = 3 = (1.6537 10-21J)3 0.03994kg/mol/6.02 10-23mol-1 = 0.5 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - szczegóły techniczne Prawie wszystkie wielkości są rzędu jedności (mniej wiecej: 10-3 - 103) - pozwala to wykrywać błedy (pojawienie się liczby o nietypowej wielkości) oraz zmniejsza błędy numeryczne (nie ma ryzyka niedomiaru (underflow) czy nadmiaru (overflow)). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - Równowaga kontra równowaga szczegółowa Jeżeli wybierać cząstkę i jej przemieszczane pojedynczo i następnie akceptować-odrzucać je, to wykonanie ruchu odwrotnego do wykonanego poprzednio jest takie same - warunek równowagi szczegółowej jest zachowany. Ale jeżeli przemieszczać wszystkie cząstki kolejno a nastepnie dokonywać akceptacji to prawdopodobieństwo ruchu odwrotnego do ruchu danej cząstki jest zerowe. W takim wypadku równowaga szczegółowa jest pogwałcona ale równowaga nadal jest zachowana i próbkowanie jest prawidłowe (zwykle). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - Równowaga kontra równowaga szczegółowa Większość algorytmów, które nie spełniają warunku równowagi szczegółowej jest błędna; szczególnie należy uważać na algorytmy łączące różne ruchy próbne. Prawdopodobnie nie jest znany przypadek, gdy algorytm zachowujący jedynie równowagę jest znacząco szybszy od zachowującego równowagę szczegółową. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - ruchy próbne Sposób generowania ruchów próbnych jest to problem wyboru macierzy ą. Optymalny schemat próbkowania - taki, który daje najmniejszy błąd statystyczny danej wielkości w danym czasie. Zatem nie jest to jednoznaczne pojęcie dla danej symulacji Bardziej praktyczną definicją będzie - suma kwadratów wszystkich zaakceptowanych przemieszczeń podzelona przez czas obliczeń. Uzasadnienie: błąd obserwabli jest odwrotnie proporcjonalny do liczby nieskorelowanych konfiguracji odwiedzonych w danym czasie obliczeń; a liczba ta jest miarą odległości w przestrzeni konfiguracyjnej. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - ruchy próbne Ruchy translacyjne - przemieszczamy jedynie środek masy przez dodanie liczby z przedziału [-"/2, "/2]: xi = xi + "(rand - 0.5) yi = yi + "(rand - 0.5) zi = zi + "(rand - 0.5) Ruch odwrotny jest równie prawdopodobny, więc ą jest symetryczna. Problemy: wybór wielkości ruchu " i czy przemieszczać cząstki jednocześnie czy po jednej (przypadkowo je wybierając - dla zachowania symetrii). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - ruchy próbne W przypadku układów takich jak symulacje fazy skondensowanej lepiej używać ruchów jednocząstkowych. Rozważmy N cząstek oddziaływujących zgodnie z potencjałem U(rN). Ruch próbny będzie zwykle odrzucany jeżeli "U kBT ale krok powinien być duży z nie za małym prawdopodobieństwem zaakceptowania, zatem: "x "! kBT . W przypadku ruchu jednocząstkowego:
"U 1 "2U "U = "rią + "rią"ri + ... "rią 2 "rią"ri = 0 + f (U)"ri2 + O("4). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - ruchy próbne Kreska - średnia po ruchach próbnych; - średnia po zespole statystycznym. "U H" kBT ! "ri2 H" kBT /f (U). Przy przemieszczaniu N cząstek pojedynczo, to najwięcej czasu zabiera obliczenie zmiany energii; przy zastosowaniu odpowiednich algorytmów czas ten wynosi tCPU = nN, gdzie n jest średnia liczbą cząstek oddziałującą z dana cząstką. Suma średnich kwadratowych przemieszczeń będzie proporcjonalna do N"r2 <" NkBT /f (U), więc średnie przemieszczenie na jednostkę czasu będzie proporcjonalne do kBT /(nf (U)) Jeżeli przemieszczać wszystkie cząstki jednocześnie - tCPU nie zmieni się ale suma średnich kwadratowych przemieszczeń będzie mniejsza o czynnik 1/N. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - ruchy próbne Uwaga! Powyższe rozważania zakładają, że każdy ruch kolektywny składa się z niezależnych przemieszczeń; istnieją jednak metody wykonywania efektywniejszych ruchów kolektywnych jeżeli przemieszczenia są zależne od siebie. Wielkoćś przemieszczenia " - wiele zródel doradza by prawdopodobieństwo zaakceptowania ruchu wynosilo 50% (metoda próbkowania może determinować inna wartość). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - ruchy próbne wykres Ądelta r2ż (delta) Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - ruchy próbne Kryterium efektywności nie pozwala jednak stwierdzić czy próbkowanie jest ergodyczne. Kryterium ergodyczności (Mountain, Thirumalai). ej(t) - średnia po czasie energii cząstki j w przedziale czasu t:
t 1 ej(t) = dt ej(t ). t 0 Średnia energia cząstki wynosi N
1 (t) = ej(t). N j=1 Wariancja N
1 2 E (t) = [ej(t) - (t)]2. N j=1 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - ruchy próbne Jeżeli wszystkie cząstki próbkują całą przestrzeń konfiguracji, to 2 E (t) 0, gdy t " 2 2 E (t)/E (0) E /t, gdzie E jest czasem charakterystycznym otrzymania nieskorelowanych próbek. Dobrą metodą optymalizacji efektywności MC jest minimalizacja iloczynu E i tCPU na ruch próbny. Dla potncjału LJ okazało się, że prawdopodobieństwo zaakceptowania ruchu 20% jest dwukrotnie efektywniejsze niż 50%. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - ruchy próbne Ruchy zmieniające orientację Molekuły sztywne liniowe. Jeżeli orientacja cząsteczki i dana jest wektorem jednostkowym i, to wektor nowej orientacji dany jest wzorem: t = , t = łĆ + i, v i t gdzie Ć jest wektorem jednostkowym wygenerowanym losowo. v Translacje i rotacje można wykonywać jednocześnie lub nie. W drugim przypadku dobrze jest wybierać typ ruchu próbnego w sposób przypadkowy. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - ruchy próbne Molekuły sztywne nieliniowe. Rotacja ciała sztywnego może być opisana kwaternionem 2 2 2 2 Q a" (q0, q1, q2, q3), q0, q1, q2, q3 = 1. Rotacja wektora sztywno przymocowanego do molekuły opisana jest macierzą ł ł 2 2 2 2 q0 + q1 - q2 - q3 2(q1q2 - q0q3) 2(q1q3 + q0q2) ł ł 2 2 2 2 R = 2(q1q2 + q0q3) q0 - q1 + q2 - q3 2(q2q3 - q0q1) . ł łł 2 2 2 2 2(q1q3 - q0q2) 2(q2q3 + q0q1) q0 - q1 - q2 + q3 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - ruchy próbne Molekuły giętkie. Ważne jest czy niektóre stopnie swobody zostały zamrożone. Jeżeli nie - wykonuje się normalne ruchy próbne poszczególnych atomów w cząsteczce. Jeżeli są sztywne więzy - powinno się robić małe ruchy próbne dla tych stopni swobody ale najlepiej jest użyc metody MD. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - metoda średnich blokowych Dane z MC nie są niezależne ale skorelowane. Średnia blokowa (block average) - średnia po pewnym przedziale czasu
tB 1 B = dtA(t). tB 0 Podczas symulacji można gromadzić średnie blokowe dla danej długości bloku tB a następnie obliczać je dla długości n tB przez zwykłe uśrednianie. Wariancja dla danej długości bloku: nB
1 (B) = (B - A )2. nB b=1 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - metoda średnich blokowych c Jeżeli tB tA, to c tA 2(B) H" A2 - A 2 . tB c Czas korelacji tA jest jednak nieznany; wielkość 2(B) P(tB) = tB A2 - A 2 c c w granicy dużego tB (tB tA) zmierza do tA; z wykresu P c od tB (lub 1/P od 1/tB) wyznacza się tA a więc i oszacowanie błędu . Jeżeli P(tB) jest silnie zmienna w granicy długiego czasu, to symulacja jest za krótka. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - metoda średnich blokowych Metoda Flyvbjerga-Petersena: niech {Ai}L będa kolejnymi i=1 próbkami wielkości A (zebranymi po ustaleniu się równowagi). Gdyby nie korelacja można by użyć zwykłych estymatorów L
1 A H" = Ai, L i=1 L
1 2(A) = A2 - A H" [Ai - ]2. L i=1 Użycie w tych wzorach coraz większych bloków prowadzi do usunięcia korelacji między blokami. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - metoda średnich blokowych Dane po transformacji A = 0.5(A2i-1 + A2i) i L = 0.5L nadal dają tą samą średnią. Nowa wariancja wynosi: L
1 2(A ) = A 2 - A H" [A - ]2. L i i=1 Powtarzamy tą operacę aż 2(A ) H" const. L - 1 Podobnie można obliczyć wariancję tej wielkości:
2(A ) 24(A ) 2(A) H" ą . L - 1 (L - 1)3 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo - metoda średnich blokowych Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo Średnie biegnące n
A n = Ai, i=1 gdzie n " (1, L). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dynamika Monte Carlo Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wyklad 6 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC W mechanice kwantowej propagator dany jest nastepującym wyrażeniem:
x iSL(x x ) K(x, x ) = DŻ exp , SL(x x ) = L(Ż Ż x x, )dt,
x 0 które w postaci zdyskretyzowanej (discretized path integral representation, DPI) [1,2] przyjmuje postać ł łł P/2 DN
mjP ł ł K(x, x ) = lim dx2 . . . dxP P " 2Ąi j=1 ł łł P DN P 2
iP i ł exp mj xjn+1 - xjn - V (xn)ł . 2 P n=1 j=1 n=1 1. R.P. Feynman and A.R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill, New York, 1965 2. R.Q. Topper, Adaptive Path Integral Monte Carlo Methods for Accurate Computation of Molecular Thermodynamic Properties, in Adv. Chem. Phys. 105, 117 (1999) Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Funkcja podziału w reprezentacji położeniowej:
" Q = dxne-$n = dx x|n n|x = dx x||x Ć Ć x n n
P
= dx1 dx2 dxP xn|e-$/P|xn+1 , xP+1 = x1, n=1 skąd wynika następująca postać operatora gęstości: = [e-$/P]P. Ć Przybliżenie wysokotemperaturowe/prymitywne - jeżeli /P 1, to ze wzoru Bakera-Campbella-Haudorffa [1] wynika: Ć Ć Ć e-$/P = e-V /(2P)e-T /Pe-V /(2P) + O((/P)3) Ć ! xn|e-$/P|xn+1 = e-V (xn)/(2P) xn|e-T /P|xn+1 e-V (xn)/(2P). 1. J. A. Barker, J. Chem. Phys. 70, 2914 (1979) Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Rozwinięcie |xn w reprezentacji pędowej daje z pomocą wzoru Lie-Trottera [1,2] DPI Q = lim QP , P " ł łł P/2 DN
P 1 ł exp - mj xjn+1 - xjn + V (xn)łł . ół 2 22 P n=1 j=1 n=1 Porównując z K(x, x) - Q jest całką po trajektoriach po orbicie periodycznej z czasem zespolonym " = -i . 1. B.J. Berne and D. Thirumilai, Annu. Rev. Phys. Chem. 37, 401 (1986) 2. H. Kleinert, Path integrals in Quantum Mechanics, Statistics, Polymer Physics and Financial Markets, World Scientific Publishing Co. Pte. Ltd. 2004 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC DN DN DN
P ip(xn-xn+1) DN dp Ć2/(2mj P j=1 n=1 = e e- pj ) (2Ą )DN ł łł P/2 P DN DN
mjP mjP(xn - xn+1)2 ł ł = exp 2Ą 2 2 2 j=1 n=1 j=1 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Potencjał efektywny: Ś(x1, ..., xP) = ŚK (x1, ..., xP) + ŚV (x1, ..., xP). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Interpretacja potencjału efektywnego (duże i małe P) Naturalny krok czasowy (wielkość przemieszczenia, time step)
2 2mP 1 x2 Krok kinetyczny = (exp(- x2) = exp(- )). 2mP 2 2 2 2 Krok potencjalny Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Rodzaje ruchów: Ruchy całego łańcucha (whole chain moves) Jednoelementowe ruchy typu klasycznego (classical single slice moves) Wieloelementowe ruchy typu klasycznego (classical multi slice moves) Jednoelementowe ruchy typu cząstki swobodnej (free particle single slice moves) Wieloelementowe ruchy typu cząstki swobodnej (free particle multi slice moves ;staging algorithm) Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Ruchy wieloelementowe są konieczne: Gdy wielkość łańcucha rośnie zastosowanie tylko ruchów jednoelemetowych spowoduje, że układ bedzie poruszał się coraz wolniej przez przestrzeń konfiguracyjną. Rozważając równanie dyfuzji Smoluchowskiego dochodzi się do wniosku, że dyfuzja środka masy c jest następująca: 1 A (c)2 = (r)2 . P2 P3 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Jednoelementowe ruchy typu klasycznego: Wybierz losowo jeden element łańcucha. Wygeneruj nowe położenie zgodnie z dystrybucja jednorodną: new xk = xk + [2 " rand() + 1] " ". Akceptacja ruchu acc(o n) = min(1, exp{-[Ś(n) - Ś(o)]}). Powtórzyć poprzednie kroki P razy (cykl - sweep). Oblicz estymatory żądanych wielkości. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Jednoelementowe ruchy typu cząstki swobodnej: Wybierz losowo jeden element łańcucha. Wygeneruj nowe położenie zgodnie z dystrybucja jednorodną:
Akceptacja ruchu N(n)ą(n o) exp{-Ś(n)} exp{-ŚK (o)} acc(o n) = = N(o)ą(o n) exp{-Ś(o)} exp{-ŚK (n)} exp{-ŚV (n)} exp[-(/P)V (xnew )] K = = . exp{-ŚV (o)} exp[-(/P)V (xK )] Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Powtórzyć poprzednie kroki P razy (cykl - sweep). Oblicz estymatory żądanych wielkości. Dla dużego P potencjał V /P staje się bardzo mały oraz xk-1 H" xx+1 H" xk, więc prawie wszystkie ruchy próbne są akceptowane. Nowe xk nie zależy od starego. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Wieloelementowe ruchy typu cząstki swobodnej: Wybierz losowo fragment łańcucha x1, x2, ..., xj, xj+1. Wygeneruj nowe położenie zgodnie z rozkładem: ą(o n) = exp{-ŚK (x1, x2, ..., xj, xj+1, ..., xP)} new new new = exp{-mP/2 2[... + (x1 - x2 )2 + (x2 - x3 )2 + . + (xjnew - xj+1)2 + ...]} Powyższy wielowymiarowy rozkład normalny można zdiagonalizować i przedstawić następująco j
2 k-1 gdzie k = = 2k-12. k mP k Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Akceptacja ruchu N(n)ą(n o) exp{-Ś(n)} exp{-ŚK (o)} acc(o n) = = N(o)ą(o n) exp{-Ś(o)} exp{-ŚK (n)} new exp{-ŚV (x1, x2 , ..., xjnew , ..., xP)} exp{-ŚV (n)} = = exp{-ŚV (o)} exp{-ŚV (x1, x2, ..., xj, ..., xP)} exp{-(/P)[V (xnew ) + V (xnew ) + ... + V (xnew )]} 2 3 j = . exp{-(/P)[V (x2) + V (x3) + ... + V (xj)]} Długość segmentu j wybiera się tak aby prawdopodobieństwo zaakceptowania ruchu wynosiło 40 - 50%. Powtórzyć poprzednie kroki P razy i obliczyć estymatory. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Staging transformation: mP exp{- [(x1 - x2)2 + (x2 - x3)2 + ... + (xj - xj+1)2]} 2 2 1 n + 1 na + b 1 (a - x)2 + (x - b)2 = [x - ]2 + (a - b)2 ! n n n + 1 n + 1 n = 1: 2 x1 + x3 1 1 2 (x1-x2)2+(x2-x3)2 = [x2- ]2+ (x1-x3)2 = 2u2+ (x1-x3)2 1 2 2 2 x1 + x3 " " u2 = x2 - x2 , x2 = . 2 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC n = 2: 1 3 x1 + 2x4 1 3 1 2 (x1-x3)2+(x3-x4)2 = [x3- ]2+ (x1-x4)2 = u3+ (x1-x4)2 2 2 3 3 2 3 x1 + 2x4 " " u3 = x3 - x3 , x3 = . 3 ... n = j 1 j x1 + (j - 1)xj+1 (x1 - xj)2 + (xj - xj+1)2 = [xj - ]2 j - 1 j - 1 j 1 j 1 2 + (x1 - xj+1)2 = uj + (x1 - xj+1)2 j j - 1 j x1 + (j - 1)xj+1 uj = xj - xj", xj" = . j Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Podsumowując: j
mP exp{- (xk - xk+1)2} 2 2 k=1 j
mP mP k 2 = exp{- (x1 - xj+1)2 - uk} 2 2j 2 2 k - 1 k=2 j
1 uk = A exp{- ( )2}, 2 k k=2 gdzie 2 k - 1 k = [ ]1/2, mP k x1 + (k - 1)xk+1 " " uk = xk - xk , xk = (k = 2, ..., j). k Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Przekształcenie odwrotne: x1 + (k - 1)xk+1 xk = uk + k Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Estymatory energii: Estymator termodynamiczny - bezpośrednie różniczkowanie funkcji rozdziału:
P N P
2 1 P n+1 n Qr =C [-DNP/2] dr1. . . drPexp - mi ri -ri + V (rn) , 22 P n=1 i=1 n=1
P N P
" log Qr DNP P 1 n+1 n E = - H" + - mi(ri - ri )2 + V (rn) . " 2 22 P n=1 i=1 n=1 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Estymator wirialny: " Zamiana zmiennych x a" r/ :
P N P
2 P Qr = C dx1 . . . dxP exp - mj xn+1 - xn - V (1/2xn) . i i 2 P n=1 i=1 n=1
P
" log Qr 1 d{V [1/2xn]} E = - " P d n=1
P
1 d{V [( + ")1/2-1/2rn]( + ")}
P d" n=1 1. J. Vanicek, W.H. Miller, J. Chem. Phys. 127, 114309 (2007) Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Estymator wirialny z usuniętym ruchem środka masy: Zmiana współrzędnych ys a" rs - r1, s = 2, ..., P i następnie " podstawienie x a" y/ daje
Qr = C [-DNP/2] dr1 dy2 . . . dyP
N P-1
2 P 2 s+1 s P exp{-[ mi (yi )2 + yi - yi + (yi )2 22 i=1 s=2
P
1 + V (r1) + V (r1 + ys) ]}, P s=2 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC
Qr = C [-DN/2] dr1 dx2 . . . dxP
N P-1
2 P exp{- mi (x2)2 + xs+1 - xs + (xP)2 i i i i 2 i=1 s=2
P
+ V (r1) + V (r1 + 1/2xs) ]}, P s=2
P
" log Qr DN 1 d{V [r1 + 1/2xn]} E = - + " 2 P d n=1 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Postępując analogicznie dla kolejnych elementów łancucha i uśredniając
P
" log Qr DN 1 d{V [rC + 1/2xn]} E = - + " 2 P d n=1
P
DN 1 d{V [rC + ( + ")1/2-1/2(rn - rC )]( + ")} + , 2 P d" n=1 P s C gdzie ri = ri /P. s=1 Procedura ta usuwa ruch środka masy redukując błąd statystyczny. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Przykład oscylator harmoniczny:
1 1 E = + 2 exp - 1 1 1 lim = =! E ( 1) 0 exp () - 1 1 E = ( 1) 2
2
Krok kinetyczny: = 2mP 1 " Krok potencjalny: exp(-1m2x2) = exp(-1 x2 ) ! = 2 2 2 m Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Zastosowania PIMC: PIMC jest metodą obliczania wielkości równowagowych w zespole kanonicznym. Układy zawierające lekkie jądra w niskich temperaturach wykazują efekty kwantowe w zakresie od reżimu kwaziklasycznego (np. stały Ar) do silnie kwantowego (np. nadciekły He) przejścia fazowe z efektami izotopowymi np. stały H2; transfer elektronów i protonów (częste w reakcjach biochemicznych) dynamika klastrów wody; 3 kondensacja Bose go-Einsteina (symulacja kondensacji He - duży sukces PIMC) Solwatowane elektrony; modyfikacja struktury rozpuszczalnika; zmiana mobilności elektronu; fotochemia; przejścia metal-dielektryk obserwowane w amoniakalnych roztworach metali alkalicznych. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Kwantowa dynamika Monte Carlo - PIMC Klastry; ważne z powodu umiejscowienia pomiędzy molekułami a ciałami stałymi np. wykazują zachowania analogiczne do przejść fazowych; pokazane zostało obniżenie temperatury przejść fazowych dla H2, D2, Ne. Adsorpcja; jedno z pierwszych zagadnień to adsorpcja helu na graficie - symulacje potwierdziły eksperymentalne diagramy fazowe; adsorpcja w strukturach porowatych - zeolity. Reakcje chemiczne (stałe szybkości reakcji chemicznych i ogólniej procesów kwantowych). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wyklad 7 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wariacyjna metoda Monte Carlo (VMC) i dyfuzyjna metoda Monte Carlo (DMC) Wariacyjna i dyfuzyjna metoda Monte Carlo - są to metody oparte na funkcji falowej; są one metodami chemii kwantowej. Służą one jako wzorce dla innych metod. Uzupełniają one mniej wymagającą metodę DFT dając dokładniejsze rezultaty i głębsze zrozumienie fizyki korelacji elektronowej. Algorytmy tych metod są wewnętrznie równoległe i pozwalają badać układy z nawet ponad 1000 elektronów. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wariacyjna metoda Monte Carlo (VMC) i dyfuzyjna metoda Monte Carlo (DMC) 4 Metoda VMC - po raz pierwszy zastosowana do ciekłego He (McMillan, 1965). Metoda DMC - w wersji ze stałymi węzłami (fixed-node) po raz pierwszy zastosowana do gazu elektronowego (Ceperley, Alder, 1980). Techniki minimalizacji wariancji do optymalizacji próbnych funkcji falowych - Umrigar, Wilson, Wilkins, 1988. Czas potrzebny do wyznaczenia energii z daną dokładnością w obydwu metodach skaluje się jak N3. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wariacyjna metoda Monte Carlo (VMC) i dyfuzyjna metoda Monte Carlo (DMC) Dla małych układów DMC osiąga dokładność chemiczną tzn. 1kcal/mol H" 0.04eV /molekua. Metoda ta korzystnie skaluje się z wielkością układu - dokładność nie spada zbyt szybko z N; aby osiągnąć dużą dokładność potrzebne są dokładne funkcje falowe ale stochastyczna ich generacja powoduje, że rezultaty symulacji są wolne od błędów skończonej bazy. W metodzie VMC jakość wyników jest całkowicie zdeterminowana przez jakość wybranej próbnej funkcji falowej. Obydwie metody najlepiej nadają się do obliczeń energii ze względu na korzystną własność zerowej wariancji - gdy funkcja próbna zbliża się do dokładnej (stanu podstawowego lub wzbudzonego) fluktuacje statystyczne energii zbliżają się do zera. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wariacyjna metoda Monte Carlo (VMC) i dyfuzyjna metoda Monte Carlo (DMC) Metody te nie są tak dobre dla stanów wzbudzonych ale mimo to są stosowane z sukcesem do obliczania wielu własności stanów wzbudzonych atomów, molekół i ciał stałych. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wariacyjna metoda Monte Carlo Opiera się na zasadzie wariacyjnej; T - funkcja próbna. T , "T - muszą być ciągłe tam, gdzie potencjał ma wartość skończoną.
Całki " T , " $T oraz " $2T muszą istnieć T T T (ostatnia zapewnia, że wariancja energii będzie skończona). Zasada wariacyjna daje ograniczenie górne na energię stanu podstawowego:
" (R)$T (R)dR T
EV = E0. " (R)T (R)dR T Aby zastosować algorytm Metropolisa przedstawiamy EV w postaci
|T (R)|2[T (R)-1$T (R)]dR
EV = E0. |T (R)|2dR Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wariacyjna metoda Monte Carlo Z powyższego wzoru wynika, że próbkujemy przestrzeń konfiguracyjną zgodnie z prawdopodobieństwem
|T (R)|2/ |T (R)|2dR. Energia lokalna w każdym z tych punktów jest dana przez EL = T (R)-1$T (R). Ostatecznie średnia energia z symulacji wynosi M
1 EV H" EL(Rm). M m=1 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Wariacyjna metoda Monte Carlo Ruchy próbne otrzymuje się zwykle z rozkładu Gaussa o początku w aktualnej pozycji danej cząstki. Inne wielkości niż energia także mogą być wyrażone jako 3N-wymiarowe całki. Można obniżyć wariancję z jaką oblicza się daną wielkość modyfikując estymator tak, aby wartość oczekiwana pozostała niezmieniona. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dyfuzyjna metoda Monte Carlo Jest to metoda rozwiązywania równania Schrdingera z urojonym czasem -"tŚ(R, t) = ($ - ET )Ś(R, t), gdzie t jest rzeczywistą zmienną mierzącą czas urojony. Zasada wariacyjna dla funkcji o ustalonej powierzchni węzłowej - jeżeli energia DMC dla danego obszaru ograniczonego powierzchnią węzłową (identyfikowaną zmienną ą ą) wynosi E0 to ą E0 E0, gdzie E0 jest energią (ścisłą) stanu podstawowego. Jeżeli powierzchnia węzłowa funkcji próbnej T (R) jest taka sama jak dokładnego stanu podstawowego, to zachodzi równość powyższych energii. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dyfuzyjna metoda Monte Carlo Twierdzenie kafelkowe (tiling theorem) - stany podstawowe wszytkich obszarów należą do tej samej klasy (zdefiniowanej przez równoważność ze względu na permutacje PR). Definiując funkcję f (R, t) = Ś(R, t)(R) oraz prędkość dryftu vD(R) = " ln |T (R)| = T (R)-1"T (R), z równanie Schrdingera z urojonym czasem przyjmuje postać 1 -"tf (R, t) = - "2f (R, t)+"[vD(R)f (R, t)]+[EL(R)-ET ]f (R, t). 2 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dyfuzyjna metoda Monte Carlo Odpowiednie równanie całkowe jest następujące
f (R, t + ) = G (R ! R , )f (R , t)dR , gdzie zmodyfikowana funkcja Greena
G(R ! R , ) a" T (R)G(R ! R , )T (R )-1. W przybliżeniu krótkiego czasu
G(R ! R , ) H" Gd(R ! R , )Gb(R ! R , ),
[R - R - vD(R )]2 Gd(R ! R , ) = (2Ą)-3N/2 exp - , 2 Gb(R ! R , ) = exp{-[EL(R) + EL(R ) - 2ET ]/2}. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dyfuzyjna metoda Monte Carlo Prędkość dryftu vD(R) powoduje przemieszczanie się walkerów w kierunku wzrastającego |T |. Także zawsze, gdy walker zbliża się do powierzchni węzłowej vD rośnie i oddala go od tej powierzchni. Funkcja Greena jest przybliżona zatem czasem następuje przekraczenie powierzchni węzłowej. Przebieg symulacji poprawia akceptacja z użyciem poniższego prawdopodobieństwa
Gd(R ! R , )Gb(R ! R, )T (R)2 paccept(R ! R ) = min 1, Gd(R ! R, )Gb(R ! R , )T (R )2
Gd(R ! R , )T (R)2 = min 1, . Gd(R ! R, )T (R )2 Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Dyfuzyjna metoda Monte Carlo Rezultatem metody DMC w opisanej wersji jest rozkład opisany dystrybucją f (R, t). Energia obliczana jest z użyciem tak zwanego estymatora mieszanego
0|$|T f (R, )EL(R)dR 1
ED = = lim H" EL(Rm). " 0|T f (R, )dR M m Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Algorytm DMC 1 Wygenerowanie zbiorów położeń początkowych zgodnie z pewną dystrubucją początkową (zwykle z |T |2 używając VMC). Obliczyć ich energie lokalne. 2 Obliczyć prędkość dryftu vD każdego walkera . 3 Jeżeli jest krokiem czasowym nowe położenie oblicza się według wzoru R = R + + vD(R ), gdzie jest 3D-wymiarowym wektorem liczb losowych o rozkładzie normalnym (o średniej zero i wariancji ). 4 Sprawdzić, czy walker przekroczył powierzchnię węzłową (przez sprawdzenie znaku funkcji próbnej) - gdy tak się stało wraca się do poprzedniego położenia. 5 Akceptacja ruchu z prawdopodobieństwem paccept(R ! R ). Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Algorytm DMC 1 Obliczenie liczby kopii, które kontynuują symulację Mnew = INT ( + exp{-[EL(R) + EL(R ) - 2ET ]/2}), gdzie - liczba losowa o rozkładzie jednorodnym na przedziale [0, 1]. 2 Zachować wielkości, którymi się interesujemy (np. średnie EL ze zbioru walkerów ). 3 Po początkowym rozgrzaniu symulacji kroki 2 - 7 powtarzać do uzyskania wystarczająco małego błędu. Wartość ET jest co pewien czas dostosowywana tak, aby utrzymać średnią ilość walkerów według wzoru ET ! ET - CE ln(Mact/Mave), gdzie CE jest dodatnią stałą kontrolującą szybkość osiągania zakładanej liczby Mave. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Funkcje próbne metod VMC i DMC Obliczenie funkcji próbnej jest najbardziej wymagającym etapem metod VMC i DMC. Zatem muszą być one zarówno dokładne i łatwe do obliczenia. Używa się funkcji typu Slatera-Jastrowa (pojedynczy wyznacznik Slatera pomnożony przez całkowicie symetryczny czynnik korekcyjny Jastrowa, który optymalizuje się). (X) = eJ(X)D(X), X = (x1, x2, ..., xN), xi = {ri, i}. Wyznacznik Slatera otrzymywany jest zwykle z metod LDA lub HF. Jednowyznacznikowa postać funkcji jest zwykle wystarczająca. Wprowadzenie Pola siłowe Klasyczna metoda Monte Carlo Kwantowe metody Monte Carlo Funkcje próbne metod VMC i DMC Czynnik Jastrowa zawiera zwykle człony jedno- i dwu-ciałowe N N N
1 J(X) = (xi) - u(xi, xj), 2 i=1 i=1 j=1,j=i
u opisuje korelację międzyelektronową a - korelację elektron-jądro. Zwykle ze zmiennych usuwa się spin i zapisuje funkcję próbną z iloczynem wyznaczników z orbitalami o spinie w dół i w górę (R) = eJ(R)Dę!(r1, ..., rNę!)D!(rNę!+1, ..., rN). Korzyść - z obliczeniowego punktu widzenia lepiej jest mieć dwa mniejsze wyznaczniki niż jeden duży.