Nowoczesne Metody Symulacyjne w Technice

Nowoczesne Metody Symulacyjne w Technice 6.2 Symulacje w nanotechnologii – symulacje wybranych zagadnień ĆWICZENIE: Symulacja zachowania zestalonego ...
2 downloads 5 Views 325KB Size
Nowoczesne Metody Symulacyjne w Technice

6.2 Symulacje w nanotechnologii – symulacje wybranych zagadnień ĆWICZENIE: Symulacja zachowania zestalonego argonu w temperaturze 30 K Wstęp i przygotowanie środowiska pracy Na dzisiejszym laboratorium zapoznamy się od strony praktycznej z programem DL_POLY, (http://www.cse.scitech.ac.uk/ccg/software/DL_POLY). DL_POLY jest popularnym narzędziem do przeprowadzania symulacji metodą dynamiki molekularnej (molecular-dynamics, MD), której podstawy omówiliśmy na wykładzie. DL_POLY jest programem darmowym, do którego dostęp można uzyskać po uprzedniej rejestracji (my jesteśmy już zarejestrowani). Po zarejestrowaniu otrzymuje się kod źródłowy programu (w języku FORTRAN-90), który można skompilować zarówno do pracy na komputerach szeregowych (jednoprocesorowych, "biurkowych") jak i równoległych (wieloprocesorowych, "klastrach obliczeniowych"). Do pracy równoległej program DL_POLY potrzebuje środowiska MPI. Z uwagi na to, że będziemy wykonywali niewielkie symulacje, wystarczy nam wersja szeregowa programu, a pracować będziemy wspólnie na mieszczącym się w s. 107A niewielkim, czteroprocesorowym serwerze jaca. Dla ułatwienia dysponujemy już wersją skompilowaną programu DL_POLY, gotową do pracy. Chwilami będziemy zaglądać do podręcznika użytkownika DL_POLY. Znajdziesz go pod adresem http://tiny.pl/lj7x. Upewnij się, że stanowisko przy którym siedzisz pracuje pod kontrolą systemu Linux – jeśli komputer uruchomił się pod Windows, uruchom go ponownie, wybierając podczas rozruchu system Linux. Środowisko linuksowe (konkretnie okna X11) będą nam potrzebne gdy dojdzie do wizualizacji otrzymanych przez nas wyników. Jeśli nie pracowałeś wcześniej w systemie Linux – nie przejmuj się, w dalszej części ćwiczenia nie będzie potrzebna jego znajomość. Ponieważ będziemy korzystać jednocześnie z jednego komputera, każdy z nas dysponuje osobnym kontem. Konta te mają identyfikatory mesuser1, mesuser2, ..., mesuser20. Aby każdy wybrał inne konto, przyjmijmy system, w którym numer stanowiska determinuje numer konta. Hasło dostępu jest takie same dla wszystkich kont – poda je prowadzący. Aby dostać się na komputer jaca, musisz najpierw zalogować się na maszynie lokalnej (tej, przed którą siedzisz). W tym celu podaj swój login i hasło. Pierwsze logowanie może zająć dłuższą chwilę, więc bądź cierpliwy. Po zalogowaniu się uruchom program terminala (ikona czarnego monitora na pasku skrótów na dole ekranu). Wydając polecenia w oknie terminala zalogujemy się na maszynę zdalną. Komputer jaca nie jest dostępny bezpośrednio z poziomu komputerów laboratoryjnych – w pierwszej kolejności musimy się zalogować na wydziałowy serwer, olimp. Zaloguj się na komputer olimp pisząc w oknie terminala ssh -X olimp Upewnij się, że użyłeś "X", a nie "x". Podaj to samo hasło. Jeśli wszystko poszło dobrze, powinieneś zobaczyć na ekranie powitanie komputera olimp, a znak zachęty konsoli zamieni się na "@olimp", co oznacza, że dostałeś się na komputer olimp. Teraz zaloguj się na komputer jaca, pisząc ssh -X jaca

Jeśli wszystko poszło dobrze, powinieneś zobaczyć na ekranie powitanie komputera jaca, a znak zachęty konsoli znowu zmieni się, tym razem na "@jaca", co potwierdza, że znajdujesz się na komputerze jaca. Program DL_POLY jest już zainstalowany, do uruchomienia go wystarczy wydać (później, teraz za wcześnie!) polecenie DLPOLY Większość programów do symulacji metodą MD znacząco różni się interfejsem od większości programów inżynierskich. Z uwagi na to, że obliczenia wykonują się na ogół bardzo długo (tygodniami) i są przeprowadzane zdalnie (nie na maszynie, na której wydajemy polecenia, a często na komputerze w innym mieście a nawet na innym kontynencie) programy MD na ogół pozbawione są zupełnie (!) interfejsu graficznego. Zazwyczaj przyjmują polecenia zapisane w postaci tekstowej w plikach, generując odpowiedzi również w postaci plików. Dla osób przyzwyczajonych do intuicyjnych interfejsów windowsowych, gdzie większość poleceń wydaje się za pomocą myszy podejście to może wydawać się nieco przedpotopowe. Musimy jednak pogodzić się z nim i wygenerować pliki potrzebne do przeprowadzenia symulacji. Aby ułatwić nam pracę, niektóre etapy zostały uproszczone do maksimum. Niestety, chwilami zmuszeni będziemy do myślenia i obliczenia kilku rzeczy na boku (na kartce, w arkuszu kalkulacyjnym albo przy użyciu Mathematiki) zanim rzucimy się do wykonywania symulacji. Przygotowanie danych do symulacji W ramach naszego ćwiczenia wykonamy symulację zachowania kryształu argonu w temperaturze 30 K. W tak niskich temperaturach argon występuje w stanie stałym, w postaci kryształu fcc (tj. o siatce kubicznej, centrowanej powierzchniowo), jak na rysunku poniżej.

Rys. 1. Fragment kryształu fcc. Argon nie jest szczególnie interesującym obiektem badań, ma wszak ważną zaletę – oddziaływania międzyatomowe dla pierwiastków gazów szlachetnych mają prostą naturę i dają dobrze się opisać za pomocą nieskomplikowanej funkcji potencjału. Dlatego właśnie warto pierwsze kroki w świecie symulacji molekularnych stawiać badając argon – bo jest prosty. Mimo to, większość umiejętności, które zdobędziemy w trakcie tego ćwiczenia przenosi się bezpośrednio na bardziej skomplikowane układy. Przyjrzymy im się następnym razem!

Aby obliczenia nie trwały zbyt długo (dzielimy się wszyscy jednym komputerem), będziemy badać niewielką próbkę – 500 atomów, czyli układ 5x5x5 czteroatomowych komórek prymitywnych. Dzięki zastosowaniu periodycznych warunków brzegowych powielimy ten układ w nieskończoność wzdłuż wszystkich osi kartezjańskich – w ten sposób będziemy symulowali zachowanie nie pojedynczego 500atomowego kryształu, a masowego (ang.: bulk) argonu – czyli pozbawionego powierzchni swobodnych, nieskończenie rozległego kryształu. DEFINIUJEMY POTENCJAŁ Oddziaływanie międzyatomowe w gazach szlachetnych dobrze opisuje prosty, dwuciałowy potencjał Lennarda-Jonesa, którego wykres znajdziesz poniżej.

energia oddziaływania

rij – odległość między atomami i-tym i j-tym. σ – parametr odległości, zależny od pierwiastka. Dla argonu wynosi 3.4 Å. ε – parametr energii, zależny od pierwiastka. Dla argonu wynosi 1.005 kJ/mol.

odległość między atomami Rys. 2. Potencjał Lennarda-Jonesa (linia niebieska) jako dobre przybliżenie potencjału obserwowanego empirycznie (linia czerwona). Źródło: Wikipedia. Spróbujmy pokusić się o zinterpretowanie wykresu potencjału. Opisuje on w jaki sposób zmienia się energia potencjalna pary atomów (dimeru) argonu, w zależności od odległości między nimi. Pamiętaj również, że gradient potencjału, ze zmienionym znakiem, opisuje siłę z jaką oddziaływują atomy, tj. siła działająca na atom i-ty wyraża się przez

 

dU  F i  r =−∇ U r =− r dr

(1)

Ponieważ operujemy tu tylko dwoma atomami, kierunek działania siły pokrywa się z prostą łączącą atomy, a wartość siły jest po prostu pochodną potencjału (ze zmienionym znakiem). Zatem wzrost funkcji potencjału oznacza działanie siły przyciągającej, malenie funkcji potencjału – siłę odpychającą. W ekstremach funkcji potencjału nie działa żadna siła.

Wróćmy do wykresu. Asymptotyczne dążenie do zera dla dużych odległości oznacza, że w miarę oddalania się od siebie atomy przestają ze sobą oddziaływać (ponieważ funkcja przestaje się nachylać). Gdy dwa atomy dostatecznie zbliżą się do siebie (spójrz np. na odległość 7.0 Å), odczuwają niewielką siłę przyciągającą (na skutek tzw. oddziaływań dyspersyjnych). Siła ta rośnie w miarę, jak atomy zbliżają się do siebie, o czym wnioskujemy z rosnącej stromości zbocza potencjału. Gdy atomy zbliżą się za bardzo, sile tej przedstawia się szybko rosnąca siła odpychająca (będąca skutkiem nadmiernego nakładania się chmur elektronowych obu atomów). Siły te równoważą się w minimum potencjału – tu energia potencjalna jest minimalna, jest to więc najbardziej "korzystna" (równowagowa) odległość pary atomów. Skoro mamy już pewne pojęcie o potencjale L-J, przygotujmy program DL_POLY do symulacji z jego udziałem. Informacji o potencjale (oraz o symulowanych molekułach) program szuka w pliku o nazwie FIELD. Zajrzyj do podręcznika użytkownika DL_POLY (sekcja 4.1.3) i przyjrzyj się jak wygląda format tego pliku. Korzystając z polecenia touch FIELD, utwórz pusty plik FIELD, uruchom edytor pisząc mcedit FIELD i wpisz do pliku co następuje: Linia 1: Komentarz, który możesz wykorzystać dowolnie. Linia 2: units kJ (informuje DL_POLY, że energie wyrażamy w kJ/mol). Linia 3: molecules 1 (mówi ile rodzajów substancji mamy; u nas jest tylko jeden – argon). Linia 4: nazwa substancji, np. "argon". W kolejnych liniach podajemy ile mamy molekuł i co się składa na molekułę. Przyjmiemy najprostszą konwencję, w której na molekułę będzie składał się tylko jeden atom argonu. Mamy zatem 500 jednooatomowych molekuł. Liczbę molekuł zadajemy dyrektywą nummols. Liczbę atomów w molekule – dyrektywą atoms. W następnej kolejności opisujemy atomy składające się na molekułę (w naszym przypadku – tylko jeden atom). Opis ten jest postaci:

(bez nawiasów kątowych). Nazwę (symbol chemiczny) argonu pewnie znasz. Masę atomową argonu (w jednostkach masy atomowej) wyszukaj w sieci. Ładunek obojętnego argonu, oczywiście wynosi zero. W kolejnej linii wpisz finish, co oznacza, że skończyłeś opisywać molekuły. Nadszedł czas na opisanie potencjałów. Chcemy tylko jeden potencjał, który należy do klasy dwuatomowych potencjałów krótkozasięgowych (van der Waalsa) (sekcja 2.3.1). Informujemy o tym DL_POLY pisząc w kolejnej linii vdw 1. W następnych liniach dla wszystkich możliwych par atomów piszemy linie postaci: ...

Skoro w naszym układzie występuje tylko argon, mamy do czynienia tylko z jedną możliwą parą: argonargon. Zatem nazwy (symbole chemiczne) atomów są łatwe do odgadnięcia. Nazwę potencjału sprawdź w sekcji 2.3.1. Następnie podaj parametry potencjału (pierwszym jest ε, drugim σ). Ostatnim słowem kluczowym w pliku FIELD musi być (w oddzielnej linii) close. Brawo, udało się nam opisać molekuły i oddziaływania występujące w symulacji. DEFINIUJEMY POŁOŻENIA ATOMÓW W kolejnym kroku musimy zadać początkową konfigurację układu. Do nas bowiem należy określenie, gdzie mają się znajdować atomy w chwili t = 0. Jeśli chcemy, możemy zadać też prędkości początkowe – jeśli tego nie zrobimy, zostaną one wylosowane z odpowiedniego rozkładu tak, aby pasowały do zadanej temperatury (odpowiada nam to). W programie DL_POLY położenia wszystkich atomów zadajemy explicite, podając współrzędne kartezjańskie każdego z atomów w pliku o nazwie CONFIG. Oczywiście nie będziemy sami obliczać położeń 500 atomów. Posłużymy się programem, który dla nas generuje położenia atomów rozłożonych na siatce fcc, zauważ jednak, że nie znamy stałej sieci (a) dla kryształu argonu, tj. nie wiemy jak odległe są od siebie atomy w krysztale. Wiemy natomiast, z eksperymentu, że w okolicy 30 K ma on gęstość 1.70 g/cm3. Wykonaj poniższe obliczenia, zwracając szczególną uwagę na jednostki – dzięki temu odgadniemy stałą sieci argonu.

d

a

a

Rys. 3. Fragment sieci fcc. a – stała sieci, d – odległość między najbliższymi sąsiadami. Ile kg waży 500 atomów argonu? (1 jednostka masy atomowej to 1.6605 × 10−27 kg ) .............................................................................................................................................................................................. .............................................................................................................................................................................................. .............................................................................................................................................................................................. .............................................................................................................................................................................................. Przy podanej na poprzedniej stronie gęstości jaką objętość zajmuje fragment kryształu zawierający 500 atomów? .............................................................................................................................................................................................. .............................................................................................................................................................................................. .............................................................................................................................................................................................. Jeśli pudło symulacyjne, w którym ma się mieścić wspomniany fragment, ma kształt sześcianu, jaki będzie jego bok? .............................................................................................................................................................................................. .............................................................................................................................................................................................. Jeśli na każdy bok przypada 5 stałych sieci, jaka jest stała sieci (a)? Jaka jest odległość między najbliższymi sąsiadami (d)? Porównaj tę odległość z odległością, dla której potencjał osiąga minimum (oszacuj z wykresu albo, lepiej, znajdź zero pochodnej). Powinieneś otrzymać zgodność z błędem nie większym niż 0.5%. Jeśli nie (spodziewamy się wartości d = ok. 3.80 Å), powtórz obliczenia. .............................................................................................................................................................................................. .............................................................................................................................................................................................. .............................................................................................................................................................................................. .............................................................................................................................................................................................. .............................................................................................................................................................................................. .............................................................................................................................................................................................. ..............................................................................................................................................................................................

Znamy już stałą sieci dla kryształu argonu w 30 K. Skorzystaj z przygotowanego dla ułatwienia programu genfcc (wpisz po prostu jego nazwę w terminalu), który generuje położenia atomów, po zadaniu rozmiarów układu, stałej sieci, masy i symbolu pierwiastka. W wyniku powinieneś otrzymać plik CONFIG, gotowy do użycia w DL_POLY. Rzuć okiem na jego zawartość (mcedit CONFIG), powinieneś zidentyfikować "na oko", które linie jakie mają znaczenie. Przyjrzyjmy się teraz, czy układ atomów zdefiniowany w pliku odpowiada naszym oczekiwaniom. Łatwy sposób na wizualizację położeń atomów zapewnia program rasmol. Ponieważ jednak oczekuje on plików wejściowych w nieco innym formacie, skorzystamy z niego pośrednio. Napisz show CONFIG

... aby uruchomić skrypt, który automatycznie pokaże Ci na ekranie wygląd stworzonej przez Ciebie konfiguracji początkowej. Powinieneś zobaczyć na ekranie okno programu rasmol, podobne do tego na następnej stronie:

Rys. 4. Okno programu rasmol, którego będziemy używać do wizualizacji symulowanego układu Spróbuj poruszać myszą przy wciśniętych różnych kombinacjach lewego i prawego przycisku myszy oraz shift i alt, żeby oswoić się z działaniem programu. Spróbuj dojść do tego, jak przybliżać i oddalać podgląd, jak go obracać i przesuwać. Uwaga – wciśnięcie klawisza ctrl wraz z ruchem myszą przełącza program w tryb oglądania przekrojów przez układ. Jeśli niechcący włączysz ten tryb (poznasz go po tym, że będziesz obserwował kulki atomowe "przekrojone"), możesz go opuścić przełączając opcję Options/Slab Mode. OPISUJEMY PRZEBIEG SYMULACJI Ostatni plik potrzebny do przeprowadzenia symulacji za pomocą programu DL_POLY nazywa się CONTROL, składa się na niego szereg dyrektyw, które opisują szczegóły przeprowadzenia symulacji (np. temperaturę, liczbę kroków, długość jednego kroku i wiele innych drobiazgów). Jeśli chcesz, w sekcji 4.1.1 podręcznika użytkownika znajdziesz komplet informacji o pliku CONTROL – listę wszystkich dyrektyw, wraz z opisem ich działania i parametrami, które przyjmują, przykładowy plik CONTROL oraz listę dyrektyw, które muszą w nim obowiązkowo wystąpić. Format pliku CONTROL jest następujący: Za wyjątkiem pierwszej linii (która zawiera komentarz), każda niepusta linia zawiera dyrektywę, po której mogą następować parametry. Ostatnią dyrektywą musi być finish. Linie puste i zaczynające się od znaku komentarza (#) są ignorowane. Utwórz pusty plik CONTROL, po czym wpisz do niego dyrektywy obowiązkowe. Będą to: • temp – jej parametr specyfikuje temperaturę (w K), w której przeprowadzamy symulację, • timestep – specyfikuje krok czasowy, z jakim całkowane są równania ruchu. Bezpieczną wartością jest 2 fs (w podręczniku sprawdź jakiej jednostki spodziewa się DL_POLY i dokonaj ewentualnego przeliczenia). • steps – to liczba kroków symulacji, które chcemy przeprowadzić. Weźmy na początek 10000 kroków.

• •

• •

– powie DL_POLY, że nie interesują nas oddziaływania kulombowskie między atomami (bo mamy obojętne atomy argonu). cutoff – promień odcięcia potencjału. Być może pamiętasz z wykładu, że żeby uniknąć obliczania oddziaływań każdy-z-każdym, zakładamy że potencjał zeruje się tożsamościowo powyżej pewnej odległości (np. powyżej 8 Å, patrząc na wykres). Przyjęcie zbyt dużej wartości spowolni obliczenia, natomiast przyjęcie za małej – zaważy na wynikach (zaczniemy ignorować niezaniedbywalne oddziaływania). Przyjmij bezpieczną wartość 10 Å. Sprawdź w dokumentacji jakiej jednostki spodziewa się program i dokonaj ewentualnego przeliczenia. delr – mało ważny parametr1, wpisz tutaj 10% promienia odcięcia, czyli 1.0. finish – jako ostatnia. no electrostatics

Normalnie należałoby wybrać jeszcze zespół termodynamiczny, w którym pracujemy (ensemble). Jeśli nie podamy żadnego, DL_POLY wybiera domyślnie zespół mikrokanoniczny (NVE). Oznacza to, że stałymi w naszej symulacji będą: liczba cząstek (N), objętość pudła (V) i całkowita energia (E). Niech tak zostanie. Dodajmy jeszcze kilka dyrektyw nieobowiązkowych (uwaga, koniecznie przed dyrektywą finish). • job time – specyfikuje maksymalny czas, po którym program powinien się zakończyć. Przy długich symulacjach istotne jest, aby zakończyć program w sposób kontrolowany, zanim zostanie on na siłę przerwany z zewnątrz (przez administratora, system kolejkowy czy system operacyjny). Dzięki temu DL_POLY może zapisać wszystkie informacje potrzebne do zrestartowania symulacji w późniejszym czasie, zanim jego działanie zostanie przerwane. Co dziwne, w tej wersji DL_POLY jeśli nie zadamy tej dyrektywy, czas działania ustawia się na 0 sekund, co skutkuje wykonaniem tylko jednego kroku, dlatego zadaj tu jakąś sensowną wartość (np. 600 sekund). • rdf – specyfikuje co ile kroków zbierane są statystyki potrzebne do wygenerowania radialnej funkcji rozkładu (o której mówiono na wykładzie). Ponieważ jest to szybka procedura, śmiało możemy zadać wielkość rzędu kilku kroków. Aby otrzymać zapis radialnej funkcji rozkładu do pliku, dodaj jeszcze dyrektywę print rdf (nie bierze parametrów). • stats – określa co ile kroków zbierane są statystyki dotyczące układu. Przyjmij, że odbywać się to będzie co 10 kroków. • traj – podanie tej dyrektywy wymusza zapis do pliku HISTORY pewnych informacji o układzie (położeń, prędkości, przyspieszeń) co określoną liczbę kroków. Podając "traj 0 100 1" życzymy sobie zapis położeń i prędkości, co 100 kroków począwszy od kroku 0. • equil i scale – pamiętasz zapewne z wykładu, że w pierwszej fazie symulacji układ powinien się zrównoważyć. W tym czasie sterujemy temperaturą układu za pomocą skalowania prędkości. Ważne jest też, aby w czasie równoważenia nie zbierać statystyk. Za pomocą dyrektywy equil zadajemy liczbę kroków poświęconych na równoważenie układu, przyjmijmy 2000 kroków. Dyrektywa scale mówi co ile kroków dokonujemy skalowania prędkości, przyjmijmy, że co 50. Mamy już przygotowane wszystkie pliki potrzebne do przeprowadzenia symulacji! Przeprowadzenie symulacji i analiza wyników Symulację startujesz wydając polecenie DLPOLY. Jeśli wszystko poszło dobrze, wyniki powinieneś dostać po ok. minucie (być może trochę dłużej, jeśli inni wykonują obliczenia w tym samym czasie). Jeśli zaraz po uruchomieniu programu (w ciągu kilku sekund albo i krócej) kończy się on, coś poszło nie tak. Zajrzyj do pliku OUTPUT – to plik, w którym DL_POLY przedstawia podstawowe informacje o symulacji, również komunikaty o błędach (jeśli te wystąpiły, stosowny komunikat pojawi się na końcu pliku). Jeśli 1 Opisuje on grubość tzw. powłoki Verleta. Wielkość ta jest blisko związana z promieniem odcięcia i mówi jak często trzeba uaktualniać listę sąsiadów dla każdego atomu. Gdy dwa lub więcej atomów przesunie się o więcej niż delr/2 względem położenia przy ostatnim uaktualnieniu, listę uaktualnia się na nowo. Dzięki takiemu podejściu nie trzeba uaktualniać listy sąsiadów (szkoda czasu) w każdym kroku.

błąd jest oczywisty – napraw go i uruchom DL_POLY ponownie. Jeśli nieoczywisty – zawołaj prowadzącego. Przejrzyj pobieżnie plik OUTPUT, jeśli zamiast liczb pojawiają się w nim nieskończoności (INF) lub nieliczby (NaN) – coś poszło ewidentnie źle – zawołaj prowadzącego. Oprócz pliku OUTPUT powinieneś dostać pliki: RDFDAT, HISTORY, STATIS, REVCON i REVIVE – to wyniki Twojej symulacji. Zajmiemy się nimi za chwilę. Plik REVCON zawiera końcową konfigurację układu (kolejno: położenia, prędkości i przyspieszenia każdego atomu) – rzuć okiem, czy wartości te wydają się sensowne (brak nieskończoności, cząstki nie poruszają się szybciej niż światło, etc.). Przypuszczamy, że nasza symulacja wykonała się pomyślnie. Czas, na obejrzenie wyników! Jak zachowują się poszczególne atomy symulowanego przez nas argonu – zobaczymy za chwilę. Na razie przyjrzyjmy się istotnym dla naszej symulacji wielkościom fizycznym. Zacznij od rzutu oka na wygenerowane przez DL_POLY pliki wyjściowe z wynikami symulacji. Będą to: OUTPUT, STATIS, HISTORY, RDFDAT, REVCON i REVIVE. OUTPUT – w pliku tym zawarte są podstawowe infomacje o przebiegu symulacji, w postaci czytelnej dla człowieka. Wydruk zaczyna się od podsumowania parametrów symulacji (możesz pobieżnie sprawdzić, czy widzisz tam oczekiwane wartości – jeśli nie, prawdopodobnie pomyliłeś się gdzieś w pliku CONTROL). Dalej widzimy próbkę konfiguracji początkowej – początkowe położenia i prędkości niektórych atomów, nie przyda nam się to. Kolejną, najobszerniejszą część pliku zajmuje wydruk istotnych dla symulacji wielkości fizycznych, zapisywanych (u nas) co 100 kroków. Ciekawsze wielkości w tabeli to: step – numer kroku, time (ps) – czas symulacji w ps, cpu (s) – czas wykonywania programu, eng_tot – całkowita energia układu, temp_tot – temperatura układu, eng_cfg – energia potencjalna układu, eng_cou – energia kulombowska układu (u nas zerowa, bo operujemy obojętnymi atomami), volume – objętość układu (u nas stała, bo wybraliśmy zespół NVE), press – ciśnienie panujące w układzie (wyrażone w kbar, czyli, na oko, tysiącach atmosfer). Oprócz samych wielkości tabela zawiera też ich bieżące średnie ze 100 kroków. Na razie rzuć okiem tylko na temperaturę – oczekujemy wartości bliskich 30 K – upewnij się, że takie dostałeś. Po tabeli mamy jeszcze kilka infomacji. Jedną z cenniejszych jest przybliżony współczynnik dyfuzji (D) – pozwala on wnioskować o stanie skupienia symulowanego układu. Dla zestalonego argonu spodziewamy się wartości bardzo bliskiej zeru. Zakłada się, że przy D powyżej ok. 10-5 cm2/s mamy do czynienia z cieczą – sprawdź, czy Twój argon jest jeszcze zestalony, czy też stopił się albo i zamienił w parę, potem porównaj 30K z temperaturą topnienia i wrzenia argonu i sprawdź, czy wyniki są zgodne z Twoimi oczekiwaniami. Wnioski: ................................................................................................................................................. ................................................................................................................................................................. ................................................................................................................................................................. Dalej w pliku OUTPUT mamy jeszcze kilka informacji, które nie będą dla nas istotne. STATIS – w pliku tym zawarte są statystyki układu – średnie istotnych wielkości liczonych co pewną liczbę kroków, którą zadaliśmy w pliku CONTROL. Tu kryje się większość wyników symulacji. Jak widzisz, plik ten składa się z bloków – każdy blok zawiera informacje o 37 wielkościach fizycznych, a nagłówek bloku opisuje do którego kroku symulacji odnosi się blok (10, 20, 30, ...) i jaki czas odpowiada temu krokowi (w ps). Aby zobaczyć cokolwiek ciekawego, będziemy musieli wydobyć z tego pliku istotne informacje, bo na pierwszy rzut oka niewiele tu widać. Jeśli chcesz wiedzieć, które z 37 liczb odpowiadają którym wielkością fizycznym – zajrzyj do podręcznika użytkownika DL_POLY. HISTORY – to plik z trajektorią układu. Znajdziesz w nim "zrzut" położeń i prędkości (na ogół również przyspieszeń, ale te wyłączyliśmy aby oszczędzić miejsce dyskowe) wszystkich atomów co 100 kroków (tyle chcieliśmy). Format tego pliku jest podobny jak pliku CONFIG, tylko zawiera on więcej "klatek" (konfiguracji). Na razie nie będziemy analizować tego pliku, przyjdzie na to czas później.

RDFDAT – to plik, w którym zapisana jest radialna funkcja rozkładu naszego układu, uśredniona po wszystkich krokach, w których były zbierane statystki (czyli od pierwszego kroku po równoważeniu, do końca symulacji). Przypomnij sobie z wykładu jaki jest fizyczny sens radialnej funkcji rozkładu. REVCON – to odpowiednik pliku CONFIG, z tą tylko różnicą, że zawiera konfigurację końcową układu (wraz z prędkościami i przyspieszeniami cząstek). Przyda nam się, gdy będziemy chcieli obejrzeć do jakiego stanu układ dotarł na końcu symulacji oraz gdy postanowimy kontynuować symulację (wtedy stanie się plikiem konfiguracji początkowej, CONFIG). Jeśli chcesz obejrzeć układ atomów składający się na konfigurację końcową, napisz show REVCON

REVIVE – zawiera komplet infomacji potrzebnych do kontynuacji symulacji, jeśli zajdzie taka potrzeba. Informacje te zapisane są w trybie binarnym, więc "na oko" nic z tego pliku nie zrozumiemy, ważne że DL_POLY wie co się tam kryje. TEMPERATURA Zacznijmy naszą analizę od przyjrzenia się, jak zmieniała się temperatura w naszej symulacji. Informacje o temperaturze układu (obliczonej na podstawie energii kinetycznej cząstek) zawarte są w plikach OUTPUT i STATIS. Skorzystamy z gotowego programu extract, który wydobywa te informacje i zapisuje w oddzielnym pliku, który łatwo możemy później wykreślić. Napisz po prostu extract temperature

...dostając w wyniku plik o nazwie "temperature". Sprawdź, co zawiera. W pierwszej kolumnie powinieneś zobaczy rosnące numery kroków, w drugiej kolumnie – bieżącą temperaturę, która powinna być bliska 30 K. Utwórz wykres temperatury korzystając z przygotowanego przez prowadzącego skryptu plot, pisząc po prostu wykresl temperature

... po chwili powinieneś dostać wykres przedstawiający jak zmieniała się temperatura w czasie. Spróbuj odpowiedzieć na poniższe pytania (przywołaj wiedzę z wykładu). Dlaczego temperatura nie jest stała, mimo że zażyczyliśmy sobie dokładnie 30 K? ................................................................................................................................................................ ................................................................................................................................................................ Jaka jest (jakościowa) różnica między oscylacjami temperatury w krokach 1-2000 a tymi w krokach 2001-10000? Jak ma się to do skalowania prędkości, które ma miejsce w trakcie równoważenia układu? ................................................................................................................................................................ ................................................................................................................................................................ ................................................................................................................................................................ CAŁKOWITA ENERGIA W kolejnym kroku przyjrzymy się całkowitej energii układu. Gdy przeprowadzamy symulację w zespole NVE (tak jak tę), energia całkowita jest całką ruchu (tj. nie zmienia się w trakcie symulacji). Sprawdź, czy dla naszej symulacji rzeczywiście ma to miejsce. Skorzystaj ponownie z programu extract i skryptu wykresl, podając zamiast "temperature" słowo "engtot" (od total energy). Spójrz na wykres i spróbuj odpowiedzieć na poniższe pytania.

Dlaczego energia całkowita w ogóle nie jest stała podczas równoważenia układu? ................................................................................................................................................................. ................................................................................................................................................................. ................................................................................................................................................................. Po skończonym równoważeniu energia powinna być stała, nigdy nie jest jednak doskonale zachowana, dlaczego? ................................................................................................................................................................. ................................................................................................................................................................. ................................................................................................................................................................. Miarą tego, jak dobrze energia jest zachowana jest wartość |ΔE/E|, czyli zmiana energii między dwoma krokami podzielona przez energię w danym kroku. Oczywiście interesują nas kroki po okresie równoważenia. Zajrzyj do pliku engtot, poza okres równoważenia i popatrz "na oko" jakiego rzędu są wartości |ΔE/E|. Przyzwoite wartości |ΔE/E| są rzędu 10-7-10-6, ale my, na podstawie naszych danych, możemy zbadać tylko jak zmienia się energia na przestrzeni 10 kroków, nie z kroku na krok, więc zadowoli nas nawet 10-5. Jakiego rzędu są |ΔE/E| w naszej symulacji? ................................................................................................................................................................ ................................................................................................................................................................ Jak możemy poprawić zachowanie energii w symulacji, jeśli okaże się że skoki energii są dla nas nieakceptowalnie duże? ................................................................................................................................................................ ................................................................................................................................................................ Jeśli energia zmienia się monotonicznie (wykazuje trend rosnący bądź malejący), musi dziać się coś złego (jeśli to zespół NVE) – albo energia pojawia się w układzie znikąd, albo znika. Upewnij się, że tak nie jest, obliczając różnicę energii na końcu symulacji i w pierwszym kroku po równoważeniu. Ile wynosi ta różnica? E2010-E10000 = ............................................................................................................................................ RADIALNA FUNKCJA ROZKŁADU Teraz przyjrzymy się radialnej funkcji rozkładu (radial distribution function, pair correlation function, g(r)). Przypomnij sobie z wykładu co opisuje ta funkcja, jakie są jej cechy charakterystyczne, ewentualnie spytaj prowadzącego albo poczytaj na: http://tiny.pl/cs36 i http://tiny.pl/cs3v. Program DL_POLY oblicza radialną funkcję rozkładu dla nas i zapisuje ją w pliku RDFDAT. Pierwsza kolumna pliku zawiera kolejne wartości r, a druga odpowiadające im wartości g(r). Wykreślić te dane możesz pisząc wykresl rdf

Piki w g(r) odpowiadają odległościom, dla których obserwujemy zwiększenie lokalnej gęstości układu. Wyznacz te odległości z dokładnością 0.1 Å (patrząc do pliku RDFDAT albo odczytując je z wykresu (posłuż się myszką, odczytywanie "na oko" nie będzie szczególnie precyzyjne). Dopasuj te odległości do charakterystycznych odległości w krysztale, poprzyj swoje rozumowanie obliczeniami. Pik #1: r1 = ........... Å. Pik #2: r2 = ........... Å. Pik #3: r3 = ........... Å. Pik #4: r4 = ........... Å. Pik #5: r5 = ........... Å.

Na poniższych fragmentach kryształu zaznacz, gdzie Ci wygodnie, odległości charakterystyczne r1..r5.

Skoro mamy do czynienia z kryształem, gdzie wszystkie atomy zajmują ściśle określone położenia (węzły) to dlaczego obserwowane piki są rozmyte? ................................................................................................................................................................. ................................................................................................................................................................. ANIMACJA ZACHOWANIA UKLADU Na koniec przyjrzymy się w jaki sposób poruszają się atomy symulowanego przez nas argonu. Oczywiście nie będziemy przyglądać się wszystkim 10000 krokom, zresztą mamy zapisane konfiguracje układu tylko z co setnego kroku. Położenia wszystkich atomów (oraz prędkości, ponieważ tak sobie zażyczyliśmy dyrektywą traj) w krokach 100, 200, ..., 20000 znajdziemy w pliku HISTORY. Oczywiście nie będziemy ręcznie generowali obrazków wszystkich 100 punktów symulacji, dla których mamy położenia atomów. Skorzystaj z gotowego skryptu generuj_film

... który po kilku minutach wygeneruje film .mpg ukazujący przebieg animacji. Odtwórz film pisząc mplayer anim.mpg

Patrząc na film koncentruj się na atomach znajdujących się bliżej środka. Atomy na zewnątrz pozornie skaczą jak oszalałe, ale to tylko złudzenie powodowane periodycznymi warunkami brzegowymi (atomy pozornie przeskakują z jednego brzegu pudła na drugi). Czy obserwowane zachowanie atomów argonu odpowiada Twoim oczekiwaniom?