UKAD MODELOWANIE PROCESW FIZYKOCHEMICZNYCH 1 czym si nie

UKŁAD <-> MODELOWANIE PROCESÓW FIZYKOCHEMICZNYCH 1) czym się nie będziemy zajmować 2) I czym zajmować się będziemy • geometria i granice układu • co „siedzi” w układzie • jakie „są relacje” układu z otoczeniem • czy mamy do czynienia ze stanem równowagi, czy też nie; jeśli nie, to czy proces jest np. stacjonarny • w jaki sposób rzeczywisty układ, zawierający ~1023 elementów „zmieścić” w komputerze • jaką zastosować metodę modelowania Geometria i granice układu: jak najbliżej rzeczywistości (np. układy objętościowe, układy niejednorodne (w polach zewnętrznych), pory, membrany, stopy, transport przez „rurki”, itp. )

W układzie „siedzą” cząsteczki. Cząsteczki owe, ich „zachowanie się” mamy modelować z x 1, y 1, z 1+kąty Eulera 1 y x 2, y 2, z 2+kąty Eulera 2 x przykład cząsteczek wody na czerwono: samemu uzupełnić wiadomości

oscylacja wiązań rotacja cząsteczka wody posiada też moment dipolowy i moment kwadrupolowy, tzn. + dipol -

polimer, zbudowany z grup atomów, fragment łańcucha w powiększeniu obok ILOŚĆ ZMIENNYCH (współrzędnych) niezbędnych do opisu ruchu cząsteczki zależy od przyjętego modelu. Jeśli cząsteczkę polimeru „zastąpimy” kulką, to wystarczą nam 3 współrzędne. Jeśli każdy segment polimeru potraktujemy jako kulkę i przyjmiemy, że polimer jest „giętki”, tzn. segmenty mogą się względem siebie poruszać, to liczba współrzędnych wyniesie 3*M, gdzie M jest ilością segmentów. Jeśli zaś polimer będzie „sztywny”, tzn. segmenty nie będą zmieniały swych wzajemnych położeń - to liczba współrzędnych =6 (jak dla ciała sztywnego o dowolnej geometrii)

Mechanika klasyczna; efekty kwantowe uwzględniamy jedynie na poziomie jednocząsteczkowym, tzn. w jednocząsteczkowych sumach stanów! W istocie w rozważanych później procesach efekty kwantowe prowadzą do zaniedbywalnie małych poprawek. Wyjątkiem są bardzo lekkie pierwiastki: wodór, hel i częściowo neon w niskich temperaturach. Nawet przy modelowaniu reakcji chemicznych nie będziemy brać pod uwagę efektów kwantowych. Dokładniej: potencjały oddziaływań używane w obliczeniach są otrzymywane w wyniku obliczeń kwantowo-mechanicznych. Jednak równania ruchu są klasyczne (równania Newtona), a nie równania Schroedingera z czasem - „pomijamy” więc zasady nieoznaczoności. W przypadku ruchów elektronów tak postapic nie wolno, ale w przypadku ruchów atomów, wiele tysięcy razy cięższych od elektronów, nie prowadzi do błedów

Podsumowując: „to co siedzi w układzie”, opis cząsteczek zależy od nas - od tego, co chcemy badać. Im bardziej „mikroskopowy” opis cząsteczek - tym trudniejsze i bardziej czasochłonne obliczenia. Z drugiej strony - im dokładniejszy model cząsteczek, tym poprawniejsze wyniki. Model jest zawsze kompromisem pomiędzy dokładnością a możliwością lub/i czasem obliczeń 1

„relacje” układu z otoczeniem, stan równowagi czy też „steady state” (przemian turbulentych badać nie będziemy) 1) Stan równowagi - zespoły statystyczne • mikrokanoniczny • wielki kanoniczny • izobaryczno-izotermiczny • ugólniony (Gibbsa), itp. 2) Stan stacjonarny - w jaki sposób stan taki się ustala czyli jak ustala się np. stała srednia liczba cząsteczek („przypływ i odpływ cząsteczek”) oraz jak zachodzi dysypacja energii. Przykład: stacjonarny przepływ płynu przez rurkę, kapilarę, por (slip-no slip) Przykład: materiały granularne (pseudo-ciecze; fluidyzacja)

„sztuczne” powiększanie układu. Elementy układu (cząsteczki) oddziałują ze sobą. Oddziaływanie ma określony zasięg. Nawet w przypadku oddziaływań elektrostatycznych (proporcjonalnych do 1/odległość) istnieje taka odległośc, na której energia oddziaływań staje się praktycznie równa zero. Zatem od pewnej odległości elementy te (cząsteczki) stają się niezależne (w sensie statystycznym). R R R - zasięg oddziaływania R R

tzw. periodyczne warunki brzegowe

periodyczne warunki brzegowe dla danego układu są związane z jego geometrią. Dla przykładu, jeśli rozważamy płyn (płyn=(gaz. or. ciecz); płyn=(gaz||ciecz)) w porze szczelinowym (co to znaczy? ), to wówczas warunki periodyczne jedynie w płaszczyźnie (powiedzmy OXY); dla układu objętościowego - warunki periodyczne w trzech wymiarach. Z warunków periodycznych wynika, że gdy cząsteczka opuszcza układ z jednej strony, to musi wejść do niego z drugiej

Innymi słowy, jeśli układ (na pł. OXY) jest prostokątem: 0 <=x<=XL (x. ge. 0). and. (x. le. XL) 0<=y<=YL (y. ge. 0). and. (y. le. YL) i na początku ruchu cząsteczka była w punkcie x, spełaniającym powyższe warunki, a na końcu ruchu znalazła się w punkcie x_koniec> XL to w istocie pojawiła się w punkcie x_koniec 1=(x_koniec-XL)! Analogicznie, jeśli x_koniec <0 to w istocie x_koniec 1=XL+x_koniec (podobnie w kierunku OY)

XL YL x 1, y 1 -poczatkowe xk. yk -końcowe przesx, przesy - przesunięcie (0, 0) OXL=2. 0/XL OYL=2. 0/YL xk= x 1+przesx yk= y 1 +przesy xk= xk -XL*aint(OXL*xk -1. ) yk= yk -YL*aint(OYL*yk -1. ) fortran OXL=2. 0/XL; OYL=2. 0/YL; xk= x 1+przesx; yk= y 1 +przesy; xk= xk -XL*floor(OXL*xk -1. 0); yk= yk -YL*floor(OYL*yk -1. 0), C

Z warunkami periodycznymi wiąże się zasada minimalnej odległości oddziaływań (minimum image) x 2’, y 2’ dx’=x 1 -x 2’ dy’=y 1 -y 2’ odl’=sqrt(dx’*dx’+dy’*dy’) x 1, y 1 x 2, y 2 x 2’’, y 2’’ dx’’=x 1 -x 2’’ dy’’=y 1 -y 2’’ odl=sqrt(dx’’*dx’’+dy’’*dy’’) dx=x 1 -x 2 dy=y 1 -y 2 odl=sqrt(dx*dx+dy*dy) ze wszystkich - wybieramy najbliżej położoną

c** XL, YL - wymiary ukladu xl 2=XL/2. 0 yl 2=YL/2. 0 dx=x 1 -x 2 dx=xl 2 -abs(dx)) dy=y 1 -y 2 dy=yl 2 -abs(dy)) odl 1=sqrt(dx*dx+dy*dy) /* XL, YL - wymiary ukladu*/ xl 2=XL/2. 0 ; yl 2=YL/2. 0; dx=x 1 -x 2; dx=xl 2 -fabs(dx)); dy=y 1 -y 2; dy=yl 2 -fabs(dy)); odl 1=sqrt(dx*dx+dy*dy); FORTRAN C

Metoda modelowania: 1) dynamika molekularna: zalety i wady 2) metoda Monte Carlo - równoważność (za wyjątkiem układów nie będących ergodycznymi) średnich czasowych i średnich w zespole statystycznym 3) dynamika „brownowska”; układy ze stopniami swobody „działającymi” na różnej skali czasowej Koncentrować będziemy się na metodach 1) i 2)

W metodzie Monte Carlo możemy modelować układ w przestrzeni ciągłej (współrzędne cząsteczek są liczbami rzeczywistymi), lub nieciągłej. W ostatnim przypadku mówimy o tzw. UKŁADACH SIATKOWYCH, a współrzędne cząsteczek są wówczas liczbami naturalnymi „numerami” węzłów sieci. Istnieje też „siatkowy analog” dynamiki molekularnej tzw. dynamika boltzmannowska.

Zadania do wykonania: wykonanie zadań jest obowiązkowe, niewykonanie 1/3 zadań oznacza brak zaliczenia wykładu. Termin wykonania zadań: 1 tydzień od podania na wykładzie; po tym terminie zadania przyjmowane nie będą. Zadania należy przesyłać pocztą elektroniczną na adres: czmpf@hermes. umcs. lublin. pl podając jako temat (TEMAT: BEZ POLSKICH ZNAKÓW!!!) numer zadania (lub zadań, jeśli więcej niż 1), imię i nazwisko Za błędy popełnione przy wysyłaniu poczty całkowitą odpowiedzialność ponosi wysyłający, żadne reklamacje przyjmowane nie będą. Pierwszy etap automatycznego sprawdzania polegać będzie na eliminacji rozwiązań zgodnych ze sobą w więcej niż w 40%, zgodnie z kryteriami programu PLAGIAT, stosowanymi w UMCS. Ponieważ etap ten jest automatyczny, nie istnieje możliwość reklamacji i odrzucone będą wszystkie prace nie spełniające powyższego kryterium.

Zadanie 1 odszukaj i podaj (powołując się na źródło) przynajmniej 5 przykładów rozmiarów cząsteczek, traktowanych jak kulki. Odszukaj i podaj (powołując się na źródło) przynajmniej 3 przykłady budowy cząsteczek wieloatomowych (średnice atomów lub grup funkcyjnych, odległości pomiędzy nimi, kąty między wiazaniami, itp. ) Odszukaj i podaj (powołując się na źródło) przynajmniej 3 przykłady danych na temat budowy ciał stałych (rodzaj sieci, wielkość komórki elementarnej, itp. ) Zadanie 2 Jaka jest struktura (jak są ułożone) kulki na płaszczyźnie, aby ich gęstość dwuwymiarowa była największa? Ile wynosi ta gęstość. Ile cząsteczek {azotu lub wody lub argonu lub tlenu lub rtęci lub metanu lub ksenonu}* znajdzie się na 1 cm 2 powierzchni, jeśli tworzą one warstwę najgęściej upakowaną? *jedna cząsteczka do wyboru

Zadanie 3 Co to jest moment dipolowy, podaj przykłady i wartości liczbowe momentów dipolowych dla 3 cząsteczek. Dla jakich wzajemnych orientacji energia oddziaływania 2 dipoli jest największa a kiedy najmniejsza? Co to jest moment kwadrupolowy? Podaj 2 przykłady. Dla jakich orientacji a)na płaszczyźnie b)w przestrzeni trójwymiarowej energia oddziaływania 2 kwadupoli jest największa? Zadanie 4 Wielkości zredukowane. Podaj przykłady zdefiniowania wielkości zredukowanych w termodynamice. Co możesz powiedzieć na temat zasady stanów odpowiadających sobie i jaki jest jej związek z definiowaniem wielkości zredukowanych? Zadanie 5. Podaj definicje kątów Eulera. Jak znając np. współrzędne środka masy oraz kąty Eulera ze środka masy do dowolnego punktu ciała sztywnego, wyliczyć współrzędne kartezjańskie tego punktu?

ODDZIAŁYWANIA MIĘDZYCZĄSTECZKOWE. Energia potencjalna (będziemy mówić potencjał) oddziaływań pary cząsteczek w określonej konfiguracji R 1 R 2 to różnica energii kwantowego stanu podstawowego dla tej konfiguracji oraz energii stanu, w którym 2 cząsteczki są nieskończenie od siebie odległe Obliczenia kwantowe -> aproksymacje, w termodynamice statystycznej musimy stosować aproksymacje (czas obliczeń!) Przedstawimy kilka najczęściej stosowanych aproksymacji

ADDYTYWNOŚĆ (po parach) Jeśli uij jest potencjałem międzycząsteczkowym pomiędzy cząsteczkami i oraz j, to energia oddziaływania N cząsteczek jest równa: taka sama zasada dotyczy obliczania energii oddziaływania 2 cząsteczek, zbudowanych z wielu atomów (grup atomów - grup funkcyjnych). Energia 2 takich cząsteczek jest sumą energii wszystkich oddziałujących centrów (grup). Centra energetyczne nie muszą być realnymi atomami. Np. energia dipol-dipol może być wyznaczona jako suma oddziaływań kulombowskich pomiędzy 4 (niekoniecznie całkowitymi) ładunkami:

q- - + r 1 + q+ q+ - q- q-q+/r 1 w modelu cząsteczli wody ładunek + może być umieszczony pomiedzy dwoma wodorami, w miejscu gdzie „fizycznie” nie ma żadnego atomu + -

Potencjał Lennard-Jonesa parametry: energii i wielkości, sigma to tzw. średnica Lennard. Jonesa. Ponieważ we wszystkich wyrażeniach termodynamicznych epsilon występuje jako epsilon/k. B (stała Boltzmanna), to zazwyczaj podajemy wartość tego parametru jako epsilon/k. B. Wówczas jednostką tej wielkości jest TEMPERATURA. Potencjał Lennard-Jonesa jest dwuparametrowy; wielkości sigma i epsilon mogą być więc użyte jako jednostki długości i energii - a więc do definiowania wielkości zredukowanych, np. gestosc_zredukowana=gestosc*sigma**3 temperatura_zredukowana=k. B*temperatura/epsilon cisnienie_zredukowane=cisnienie*sigma**3/epsilon energia_wewnetrzna_zredukowana=energia_wewnetrzna/epsilon, itd.

u( r)/epsilon force sigma rcut u(rcut)=ucut r/sigma

Funkcja wyznaczająca potencjał LJ w jednostkach double lj(r, rmin, emax) zredukowanych double r, rmin, emax; real*8 function lj(r, rmin, emax) { double r, r 3, r 6; implicit none { if (r<=rmin); real*8 r, r 6 return(emax); if (r. le. rmin) then } lj=emax else { r 3=1. 0/r; else r 3=r 3*r 3; r 6=(1. d 0/r)**6 r 6=r 3*r 3; lj=4*r 6*(r 6 -1. d 0) r 3=4*r 6*(r 6 -1. 0); endif } return(r 3); return } end Uwaga: dlaczego rmin oraz emax? Jak je określić?

oddziaływanie z powierzchnią - na bazie potencjału LJ OZ x 1, z 1 r dx*dy OXY

całkujemy na całą płaszczyznę, od minus do plus nieskończoności Przechodzimy na współrzędne walcowe R 2=(x 1 -x)2+(y 1 -y)2; Rd. Rd =dxdy Potencjał Lennard-Jonesa (10, 4), opisujący oddziaływania z płaszczyzną

Inne modele potencjałów oddziaływa ń pary cząsteczek sferycznych (PARY CZĄSTECZEK, A NIE CZĄSTECZKI Z POWIERZCHNIA) 2) Potencjał Yukawy 3) potencjał sztywnych kul (BYŁO!) Kiedy stosujemy te potencjały?

Jakie zadania mogą być z tego materiału na kolokwium? Np. (a) przy wyznaczaniu potencjału oddziaływań z powierzchnią przeprowadzić analogiczne całkowanie dla potencjału Yukawy, napisać odpowiedni program, wykonać wykres. Całkowanie można przeprowadzić numerycznie (b) wyznaczyć numerycznie potencjał oddziaływania cząsteczki ze sferyczna cząsteczka koloidalną (całkowanie po płaszczyżnie, ale po powierzchni kuli) c ) jak wyżej, ale „w srodku” kuli, tj. wewnątrz poru sferycznego lub cylindrycznego d) wyznaczyć analog potencjału LJ(10, 4), jeśli uwzględnimy, że oddziaływanie zachodzi również z elementami wewnątrz ciała stałego e) wyznaczyć numerycznie drugi współczynnik wirialny dla potencjału LJ(12, 6) f) wyznaczyć numerycznie wartość stałej Henry’ego dla potencjału LJ(10, 4).

Zadania do wykonania 2 tygodnie od 24 lutego Zadanie 6 Jaką postać ma potencjał opisujacy oddziaływania pomiędzy dwoma idealnymi (punktowymi) dipolami? Napisz funkcje obliczającą ten potencjał. Zadanie 7 Jaką postać ma potencjał opisujacy oddziaływania pomiędzy dwoma idealnymi kwadrupolami ? Napisz funkcje obliczającą ten potencjał. Zadanie 8 Oddziaływania pomiędzy dwoma cząsteczkami koloidalnymi w roztworze atermalnym można opisać przy pomocy potencjału gaussowskiego (tzw. Gaussian overlap potential) lub Asekury Oozawy. Jaką postać ma jeden lub drugi potencjał. Napisz odpowiednią funkcję. lub - do wyboru

DYNAMIKA MOLEKULARNA - UKŁADY OBJĘTOŚCIOWE (bez pól zewnętrznych) Notacja: matematyczna i „komputerowa” N cząstek j rj ri=(xi, yi, zi), uij=u(rij) ri=(xi, yi, zi); uij=u(rij) i u(rij) ri r j’ rij=sqrt(dx 2+dy 2+dz 2) dx=min_image(xi-xj) dy=min_image(yi-yj) dz=min_image(zi-zj)

un=suma(N, uij), un jest funkcją położeń wszystkich N cząsteczek Siła działająca na i-tą cząsteczkę jest sumą sił od wszystkich pozostałych N-1 cząsteczek. Równania ruchu -klasyczna dynamika Newtona (przestrzeń trójwymiarowa, czasteczki sferyczne): siła=masa*przyspieszenie siła_na cząsteczkę_i=masa_i*przyspieszenie_i=(druga_pochodna_drogi_i_po_czasie) siła_i= - (gradient_potencjału_działającego_na_i) gradient=( / x, / y, / z) -fxi=mi*d 2 xidt 2 -fyi=mi*d 2 yidt 2 -fzi=mi*d 2 zidt 2

Dla potencjału LJ(12, 6) dx=x(i)-x(j) dy=y(i)-y(j) dz=z(i)-z(j) dx=dx-XL*nint(dx/XL) dy=dy-YL*nint(dy/YL) dz=dz-ZL*nint(dz/ZL) rsq=dx**2+dz**2+dy**2 robol=(sigma**6/rsq**4) robol 1=(sigma**12/rsq**7) robol=6. d 0*robol-12. d 0*robol 1 fxij=4. d 0*epsilon*robol*dz nint, sumowanie, ij<> ij

Algorytm Stroemera-Verleta. Inne algorytmy: M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, Claredon, Oxford 1982 (szukaj w sieci hasła: ccp 5, na stronie odszukać program library) D. J. Evans, G. P. Morriss, Statistical Mechanics of Nonequilibrium Liquids, Academic Press, London 1990

Schemat blokowy symulacji: 1. wczytaj dane początkowe (wielkość układu, N, parametry potencjału, wielkość kroku czasowego delta_t, jak długo mają trwać symulacje, KONFIGURACJĘ POCZĄTKOWĄ) 2. Dla każdej cząsteczki oblicz siły fxij, fyij, fzij 3. Dla każdej cząsteczki rozwiąż numerycznie równania Newtona, obliczając jej nowe położenie 4. Oblicz wielkości termodynamiczne, dynamiczne oraz charakteryzujące strukturę, które chcesz obliczać i których wartość średnią (uśrednioną po krokach symulacji) będziesz drukować powtarzaj 5. wydrukuj wyniki (wielkości średnie, konfigurację końcową)

A CO Z TEMPERATURĄ? Zasada ekwipartycji energii: w stanie równowagi wartość średnia energii kinetycznej wynosi (1/2)Nk. T na każdy stopień swobody, a więc: a co to prędkość vxi, vyi, vzi? do przecież pochodna xi po czasie, czyli vxi=(xi(t+delta_t)-xi (t))/delta_t

Problemy do przypomnienia/uzupełnienia: numeryczne rozwiązywanie równań różniczkowych rzędu drugiego, zwłaszcza metoda „predictor-corrector” oraz metoda „leap frogg” - źródło: odszukać na sieci numerical recipes, pobrać pliki PDF (języki FORTRAN lub C) DYNAMIKA MOLEKULARNA - ZESPÓŁ MIKROKANONICZNY (a co to jest zespół mikrokanoniczny? )

Omówimy dokładniej program pozwalający na obliczanie trajektorii ruchu cząsteczek metoda dynamiki molekularnej Program napisany jest przy pomocy meta-języka. Program główny wygląda następująco: program md call init t=0 do while(t<=tmax) call force call integrate t=t+deltat call sample enddo call results end program główny inicjujemy dane, pkt 1 zaczynamy, czas=0, krok czasowy=delta początek pętli po czasie obliczamy siły całkujemy równania ruchu czas=czas_stary+delta obliczamy średnie, które nas interesują koniec pętli czasowej wypisujemy wyniki

We wszystkich blokach programowych niżej rozpatruje się jedynie JEDNĄ wspólrzedna - x; Jeśli układ jest dwuwymiarowy - trzeba odpowiednio dodać y; dla układu trójwymiarowego - również wsp. z. Oczywistym jest, że odległości liczymy zawsze jako długości odpowiednich wektorów r 2=dx*dx - układ jednowymiarowy r 2=dx*dx+dy*dy - układ dwuwymiarowy r 2=dz*dz+dy*dy+dz*dz - układ trójwymiarowy

procedure init bardzo ważne, inaczej układ symv=0 by dryfował, sumv 2 - skalowanie sumv 2=0 prędkości do temperatury for i=1, n_czastek x(i)=polożenie_na_sieci_i położenia na sieci przypadkowe prędkości v(i)=(random-0. 5) ZACHOWANIE PĘDU! sumv=sumv+v(i) TEMPERATURA sumv 2=sumv 2+v(i)*v(i) RÓWNOWAGOWANIE endfor sumv=sumv/n_czastek sumv 2=sumv 2/n_czastek fs=sqrt(3*temp/sumv 2) for i=1, n_czastek pętla skalująca prędkosci, zero v(i)=(v(i)-sumv)*fs wypadkowego pędu, układ ma endfor zadana temperaturę end_procedure

UWAGA jeśli stosować będziemy algorytm Verleta, to w miejscu strzałki należy dodać xm(i)=x(i)-v(i)*deltat, gdzie deltat jest krokiem czasu do obliczenia położeń w czasie t+deltat algrorytm ten wymaga znajomosci położeń w czasie t oraz w czasie t-deltat. Te ostatnie „przybliżamy”, zakładając ruch jednostajny w czasie deltat, czyli: polożenie_stare=polożenie_aktualne-predkość*deltat

SIEĆ, układ 2 -wymiarowy nsqr=sqrt(n_czastek) for i=1, nsqr x(i)=(i-1)*a for j=1, nsqr y(j)=(j-1)*a endfor jeśli na tych bokach są a to na tych bokach być nie mogą periodyczne war. brzegowe!

cząsteczka następnej warstwy (wzdłuż osi OZ, prostopadłej do płaszczyzny rys. 1’ 1 2’ 2 a sieć heksagonalna a - jednostka A jak wyglądają tu warunki periodyczne? ? 3 1 - (0, 0) 2 - (1, 0) 3 - (3, 0) itd. . 2 -warstwa: 1’ -(0. 5, sqrt(3)/2) 2’ - (1. 5, sqrt(3)/2) idt. . 3 -wartstwa 1’’ - (0, sqrt(3)) 2’’- (1, sqrt(3))

1 2 1. wstawiamymy 1 2. wstawiamy 2, aby r 12>a 3 4 3. wstawiamy 3, aby r 21, r 32>a 4. wstawiamy 4, aby r 41, r 42>a. itd. . Ale: sprawdzać należy Przypadkowe wstawianie stosując periodyczne warunki brzegowe!

Zadanie 9 Co to jest rozkład prędkości Maxwella-Boltzmanna? Dlaczego w procedurze inicjującej nie losujemy prędkości zgodnie z tym rozkładem? Zadanie 10 Napisz procedurę inicjującą dla układu dwu lub trójwymiarowego, losując położenia cząsteczek na sieci kwadratowej (sześciennej) lub najgęstszego upakowania lub przypadkowo. Język: fortran lub c

Procedura force Jeśli algorytm symulacji działa poprawnie, to cząsteczki nie mogą „za bardzo” zbliżyć się do siebie, nie ma więc potrzeby wprowadzania parametru „rmin”. Natomiast zazwyczaj układ jest większy niż (zasięg korelacji)*2, dlatego wprowadzamy parametr rcut wartość potencjału LJ(12, 6) w pkt. rcut to ucut=4*(1. 0/rcut 12 -1. 0/rcut 6), PAMIETAJ: epsilon i sigma potencjału LJ(12, 6) to nasze JEDNOSTKI W istocie potencjał LJ(12, 6) w jest NIECIĄGŁY w rcut, bo do r>rcut u(r )=0, ale ucut jest różne od zera. To oznacza, że również jego pochodna du( r)/dr na DODATKOWY „wkład” w r=rcut. Jeśli funkcja w danym punkcie ma „schodek”, to pochodna w tym punkcie ma człon z funkcją delta Diraca! Ten człon jest mały (tym mniejszy, im rcut większe) i dla uproszczenia będziemy go pomijać. W dokładnych obliczeniach należy go uwzględniać

procedure force energia=0 for i=1, n_czastek f(i)=0 endfor i=1, n_czastek for j=i+1, n_czastek dx=x(i)-x(j) dx=dx-XL*nint(dx/XL) r 2=dx*dx if (r 2<rcut*rcut) then r 2 i=1. 0/r 2 r 6 i=r 2 i*r 2 i ff=48. 0*r 2 i*r 6 i*(r 6 i-0. 5) f(i)=f(i)+ff*dx f(j)=f(j)-ff*dx energia=energia+4*r 6 i*(r 6 i-1)-ucut end if endfor i endfor j end procedure zerujemy energię i siły f dla n_cza stek to jest „odległość do kwadratu” jeśli cząstka i działa na j, to j działa na i taka sama siła, ale z przeciwnym zwrotem periodyczne war. brzegowe, nint - to najbliższ całkowita sumujemy siły od wszystkich cąstecek sumujemy energię (minus ucut!)

Zadanie 10 Układ cząsteczek oddziałujących potencjałem LJ(12, 6) jest w kontakcie z powierzchnią (płaszczyzna OXY). Powierzchnia ta jest nieruchoma i oddziałuje na cząsteczki układu potencjałem v(z) opisanym funkcja LJ(10, 4). Napisz równania na siły wynikające z tych oddziaływań. Zmodyfikuj tak procedurę force, aby poprawnie obliczała siły. UWAGA: w takim układzie NIE MA PERIODYCZNYCH WARUNKÓW BRZEGOWYCH W KIERUNKU Z; aby być konsekwentnym, należy również umieścić symetrycznie taką samą powierzchnie w z=ZL i uczynić ZL dużo większym niż pozostałe wymiary układu: XL i YL

Zadanie 11 Właściwości funkcji delta Diraca: definicja, całkowanie tej funkcji, obliczanie pochodnej tej funkcji. Funkcja Diraca to analog delty Kroneckera w przestrzeni liczb rzeczywistych Źródło: np. uzupełnienia matematyczne w Mechanice Kwantowej Davydova (Dawidowa, jest po polsku)

Całkowanie równań ruchu procedure integrate sumv - sprawdzenie „dryfowania” sumv=0 zasada zachowania pędu, sumv 2=0 for i=1, n_czastek do obliczenia temperatury xx=2*x(i)-xm(i)*deltat+deltat* f(i) vi=(xx-xm(i))/(2*deltat) równania ruchu całkowane zgodsumv=sumv+vi sumv 2=sumv 2+vi*vi nie z równaniami na następnym xm(i)=x(i) slajdzie x(i)=xx endfor temp=sumv 2/(2*n_czastek) etot=(energia+sumv 2)/(2*n_czastek) end procedure całkowita energia, energia pot. (energia) liczona w force

calkowanie równań ruchu: x(t+teltat)=x(t)+v(t)*deltat+f(t)*deltat/(2*masa)+O(deltat^3) x(t-deltat)=x(t)vv(t)*deltat+f(t)*deltat/(2*masa)+O(deltat^3) dodając stronami: x(t+deltat)+x(t-deltat)=2*x(t)+f(t)*deltat/masa lub x(t+deltat)= - x(t-deltat)=2*x(t)+f(t)*deltat/masa podobnie postępując z prędkością otrzymamy: v(t)=(t(t+deltat)-t(t-deltat))/(2*deltat) powyższe równania zastosowano w procedurze integrate

RÓWNOWAGOWANIE: ENEGRGIA lub TEMPERATURA kroki czasu podczas równowagowania musimy „uzgodnić” siły z prędkościami. Zaczynamy symulacje od przypadkowych prędkosci, ich rozkład nie jest zgodny z rozkładem MB. Równowagowanie - to osiągniecie tej zgodności. Dlatego też podczas równowagowania należy (w procedurze integrate) dopisać linie z procedury inicjalizacji procedura inicjalizacji

Zadanie 12 Napisz (albo w meta-języku, albo w dowolnym języku programowania) analog procedury integrate dla równowagowania układu.

Uwagi podsumowujące: Najprostszy układ: całkowicie izolowany, bez pól zewnętrznych, siły krótkiego zasięgu, cząsteczki sferyczne. W przypadku cząsteczek niesferycznych - dodatkowe równania ruchu, np. dla rotacji, Średnia energia (kinetyczna) rotacji - związek z temperaturą. Równowagowanie termiczne - rotacja i translacja „razem” czy oddzielnie? Reakcje chemiczne, np. asocjacja, tworzenie wiązań typu wiązania wodorowego

Uzupełnienia (samodzielnie, nieobowiązkowo) Inne metody całkowania równań ruchu (inne wersje procedury inegrate) Równania Liouville’a, zagadnienie odwracalności w czasie danego algorytmu. Niestabilność Lapunowa.

Co chcemy liczyć - czyli procedura sample 1. Wielkości termodynamiczne (T, energia, potencjał chemiczny, ciśnienie (albo ugólniej: składowe tensora ciśnienia (naprężeń), a więc i napiecie powierzchniowe, napięcie liniowe), swobodna energia Helmholtza, ciepło własciwe. . . 2. Wielkości dynamiczne, np. współczynnik dyfuzji, czas „życia” powstających kompleksów. . . 3. Wielkości strukturalne: funkcje korelacji (gęstość lokalna, radialna funkcja rozkładu, liczby koordynacyjne, niezmienniki (inwarianty) strukturalne + histogramy

FUNKCJA RADIALNA

g(r) 1 odległość

prawdopodobieństwo znalezienia innej cząsteczki na odległości r 12 od pierwszej cząsteczki (jakakolwiek pierwsza, jakakolwiek „inna”) - dwucząsteczkowa funkcja rozkładu, n 2(x 1, y, z 1, x 2, y 2, z 2) Periodyczne warunki brzegowe Układ objętościowy -> n 2(x 1, y, z 1, x 2, y 2, z 2=n 2(r 12) Układ idealny (brak oddziaływań - cząsteczki nieskorelowane) n 2(r 12)=n 2 g 2(r 12), n-gęstość, g 2(r 12) - funkcja radialna Zatem dla układu bez oddziaływań -> g 2(r 12)=1 jeśli r 12 -> nieskończoność, to brak oddziaływań, czyli g 2(r 12)->1 Zatem jeśli obowiązuje zasada addytywności po parach, to termodynamikę określa funkcja radialna, np. średnia energia wewnętrzna=energia kinetyczna+energia potencjalna

Energia kinetyczna - z zasady ekwipartycji, tj. ze średniej temperatury Energia potencjalna (w 3 wymiarach): (patrz zadanie 13) gdzie potencjał u(r 12) pod całka jest np. potencjałem LJ(12, 6); rmax jest zasięgiem oddziaływań. Wykażemy dalej, że również ciśnienie oblicza się z funkcji radialnej

Dalej: liczby koordynacji oblicza się jako całki z funkcji radialnej. Jeśli pierwsza „strefa” koordynacji to zasięg minimum oddziaływania LJ(12, 6) (r=(21/6)*sigma) to liczba ta wynosi Funkcja radialna określa nam prawdopodobieństwo znalezienia pary (jakichkolwiek) cząsteczek na odległości r

procedure radial if(iradial) then {ileradial=ileradial+1 for i=1, n_czastek-1 for j=i+1, n_czastek dx=x(i)-x(j) dx=dx-XL*nint(dx/XL) r=sqrt(dz*dz) if(r>BOX/2) then { ir=int(r/step) g(ir)=g(ir)+1 }} end procedure Zmienne: procedura wywoływana być powinna po kroku równowag. wtedy iradial=prawda, przeciwnie iradial=fałsz. step to krok z jakim tablicujemy funkcje radialną, g(ir), zaś mox to min(XL, Yl, ZL) zazwyczaj step jest rzędu od 0. 01 do 0. 05 średnic cząsteczek. Tablica g(r ) powinna być tak zadeklarowa na, aby 0<=r<=BOX/2 Do wyliczenia średnich wartości g(ir) na końcu symulacji musimy wiedzieć ile razy obliczaliśmy funkcję radialną!, dlatego też zmienna ileradial powinna być wyzerowana na początku programu głównego

Wydruk funkcji radialnej Oczywiście, w 2 wymiarach gestosc=n_czastel/(XL*YL) gestosc=n_czastek/(XL*YL*ZL) for i=0, int(BOIX/step) oblicz_volume oblicza r=(i+0. 5)*step oblicz_volume objętość (lub pole dla 2 wym. ) g(i)=g(i)/(ileradial*gestosc*volume pomiedzy sferami i+1 a i * n_czastek) volume=(4/3)*pi*((r(i+1))3 -(ri)3) print r, g(i) lub endfor volume=2*pi*((r(i+1))2 -(ri)2) ri=i*step

Zadanie 13 Zakładając, że funkcja radialna dana jest wyrażeniem g 2(r )=exp[-u( r)/k. B T}, gdzie T jest temperaturą, k. B - stała Boltzmanna a u( r) - potencjałem LJ, napisz program wyznaczający średnią wartość energii potencjalnej z równania Obliczenia przeprowadź w 2 lub 3 wymiarach Zadanie 14. Wykaż, że jeżeli wyrażenie na swobodną energię Helmholtza (lub równanie stanu) dane są równaniem wirialnym, kończącym się na 2 współczynniku wirialnym, to jest to jednoznaczne z założeniem, że funkcja radialna dana jaest wyrażeniem g 2(r )=exp[-u( r)/k. B T}