UKŁAD <-> MODELOWANIE PROCESÓW FIZYKOCHEMICZNYCH 1) czym się nie będziemy zajmować

1. 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.)

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

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

4. 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)

5. 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

6. 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

7. „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 • kanoniczny • 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)

8. „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ń elektro- statycznych (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 - zasięg oddziaływania R R R

9. tzw. periodyczne warunki brzegowe

10. 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

11. 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_koniec1=(x_koniec-XL)! Analogicznie, jeśli x_koniec <0 to w istocie x_koniec1=XL+x_koniec (podobnie w kierunku OY)

12. XL x1,y1 -poczatkowe YL xk.yk -końcowe przesx,przesy - przesunięcie (0,0) OXL=2.0/XL OYL=2.0/YL xk= x1+przesx yk= y1 +przesy xk= xk -XL*aint(OXL*xk -1.) yk= yk -YL*aint(OYL*yk -1.) OXL=2.0/XL; OYL=2.0/YL; xk= x1+przesx; yk= y1 +przesy; xk= xk -XL*floor(OXL*xk -1.0); yk= yk -YL*floor(OYL*yk -1.0), fortran C

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

14. c** XL, YL - wymiary ukladu xl2=XL/2.0 yl2=YL/2.0 dx=x1-x2 FORTRAN dx=xl2-abs(xl2-abs(dx)) dy=y1-y2 dy=yl2-abs(yl2-abs(dy)) odl1=sqrt(dx*dx+dy*dy) /* XL, YL - wymiary ukladu*/ xl2=XL/2.0 ; yl2=YL/2.0; dx=x1-x2; C dx=xl2-fabs(xl2-fabs(dx)); dy=y1-y2; dy=yl2-fabs(yl2-fabs(dy)); odl1=sqrt(dx*dx+dy*dy);

15. 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)

16. 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.

17. 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ą odpowie- dzialność 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.

18. 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 cm2 powierzchni, jeśli tworzą one warstwę najgęściej upakowaną? *jedna cząsteczka do wyboru

19. 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ójwymiaro- wej 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?

20. ODDZIAŁYWANIA MIĘDZYCZĄSTECZKOWE. Energia potencjalna (będziemy mówić potencjał) oddziaływań pary cząsteczek w określonej konfiguracji R2 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 R1 Obliczenia kwantowe -> aproksymacje, w termodynamice statystycznej musimy stosować aproksymacje (czas obliczeń!) Przedstawimy kilka najczęściej stosowanych aproksymacji

21. 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:

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

23. 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/kB (stała Boltzmanna), to zazwyczaj podajemy wartość tego parametru jako epsilon/kB . 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=kB*temperatura/epsilon cisnienie_zredukowane=cisnienie*sigma**3/epsilon energia_wewnetrzna_zredukowana=energia_wewnetrzna/epsilon, itd.

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

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

26. oddziaływanie z powierzchnią - na bazie potencjału LJ OZ x1,x1,z1 r z1 dx*dy OXY

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

28. 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?

29. 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 nie 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).

30. Zadania do wykonania 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 2 tygodnie od 24 lutego

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

32. 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 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*d2xidt2 -fyi=mi*d2yidt2 -fzi=mi*d2zidt2

33. 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) robol1=(sigma**12/rsq**7) robol=6.d0*robol-12.d0*robol1 fxij=4.d0*epsilon*robol*dz Dla potencjału LJ(12,6) nint, sumowanie, ij<-> ij

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

35. 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ć 5. wydrukuj wyniki (wielkości średnie, konfigurację końcową) powtarzaj

36. A CO Z TEMPERATURĄ? Zasada ekwipartycji energii: w stanie równowagi wartość średnia energii kinetycznej wynosi (1/2)NkT 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

37. 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?)

38. 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 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 program md call init t=0 do while(t<=tmax) call force call integrate t=t+deltat call sample enddo call results end

39. 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 r2=dx*dx - układ jednowymiarowy r2=dx*dx+dy*dy - układ dwuwymiarowy r2=dz*dz+dy*dy+dz*dz - układ trójwymiarowy

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

41. 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

42. 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 endfor jeśli na tych bokach są to na tych bokach być nie mogą periodyczne war. brzegowe! a

43. cząsteczka następnej warstwy (wzdłuż osi OZ, prostopadłej do płaszczyzny rys. 2’ 1’ 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)) 2 1 3 a sieć heksagonalna a - jednostka A jak wyglądają tu warunki periodyczne??

44. 1 2 1. wstawiamymy 1 2. wstawiamy 2, aby r12>a 3. wstawiamy 3, aby r21,r32>a 4. wstawiamy 4, aby r41,r42,r42>a. itd.. Ale: sprawdzać należy stosując periodyczne warunki brzegowe! 3 4 Przypadkowe wstawianie

45. 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

46. 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/rcut12-1.0/rcut6), 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 uproszcze- nia będziemy go pomijać. W dokładnych obliczeniach należy go uwzględniać

47. procedure force energia=0 for i=1,n_czastek f(i)=0 endfor for i=1,n_czastek for j=i+1,n_czastek dx=x(i)-x(j) dx=dx-XL*nint(dx/XL) r2=dx*dx if (r2<rcut*rcut) then r2i=1.0/r2 r6i=r2i*r2i*r2i ff=48.0*r2i*r6i*(r6i-0.5) f(i)=f(i)+ff*dx f(j)=f(j)-ff*dx energia=energia+4*r6i*(r6i-1)-ucut end if endfor i endfor j end procedure zerujemy energię i siły f dla n_cza stek 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 to jest „odległość do kwadratu” sumujemy energię (minus ucut!)

48. 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

49. 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)

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