Uwagi wstępne

Gwałtowny rozwój techniki komputerowej dokonujący się w czasie ostatnich dziesięcio-leci umożliwił wykorzystanie mocy obliczeniowej komputerów do modelowania dynamiki procesów zachodzących w układach fizycznych. Celem modelowania komputerowego syste-mów złożonych z atosyste-mów i cząsteczek jest obliczenie makroskopowego zachowania układu mając do dyspozycji wiedzę na temat mikroskopowych oddziaływań pomiędzy atomami tworzącymi modelowany system. Opis makroskopowej ewolucji układu osiąga się poprzez obliczenie zmiany położeń atomów układu w czasie, otrzymując w rezultacie tzw. tra-jektorie ruchu atomów. Możliwość wizualizacji trajektorii ruchu cząstek w połączeniu z względnie szybkim czasem ich otrzymania umożliwia [ALL92]:

— zrozumienie oraz interpretację wyników doświadczalnych, — szacowanie wyników eksperymentu,

— wspomaganie badań doświadczalnych poprzez dostarczanie danych niemożliwych lub trudnych do zrealizowania w laboratorium.

Ostatni z powyżej przedstawionych punktów należy rozumieć jako możliwość bardzo pre-cyzyjnego zaplanowania oraz zaimplementowania warunków eksperymentu. W kontekście analizy powierzchni techniką SIMS/SNMS będzie to: brak zanieczyszczeń oraz defektów powierzchni, możliwość ustalenia dowolnej energii kinetycznej wiązki rozpylającej, a także rodzaju pocisków wchodzących w jej skład. Daje to zatem możliwość wyjścia poza ogra-niczenia warunków eksperymentu wynikające z posiadanej aparatury pomiarowej. Mode-lowanie komputerowe jest stosowane także do weryfikacji modeli teoretycznych [ALL92]. Przy pomocy modelowania komputerowego można zobaczyć czy układ rzeczywiście zacho-wuje się tak, jak to przewiduje teoria analityczna. Zgodność wyników stanowi pozytywną weryfikację modelu teoretycznego.

Modelowanie komputerowe posiada też szereg ograniczeń. Po pierwsze wynik symulacji komputerowej jest uzyskiwany dla bardzo szczegółowo zdefiniowanego stanu początkowe-go układu. Z tepoczątkowe-go powodu zmiana wartości nawet jednepoczątkowe-go z parametrów określających warunki symulacji (np. energii pocisku) powoduje konieczność wykonania powtórnych ob-liczeń. Tym samym, aby uzyskane wyniki miały przeglądowy charakter, trzeba poświęcić

wiele czasu, aby obliczyć ewolucje wielu różniących się nieznacznie układów. Problemem może być też szybkość obliczeń oraz precyzja obliczeń, a także konieczność składowa-nia ogromnych ilości danych. Problem czasu obliczeń jest związany z niedoskonałością komputerów. Nawet najszybszy komputer nie może prowadzić obliczeń z nieskończoną prędkością. Tym samym niemożliwe jest modelowanie ewolucji układu oraz śledzenie tra-jektorii cząstek przez tak długi czas, jak to ma miejsce w eksperymencie, gdzie rejestracja cząstek następuje po czasie rzędu mikrosekund. W symulacjach prowadzonych metodą dynamiki molekularnej trzeba ograniczyć się najczęściej do czasów rzędu pikosekund. Metody Monte Carlo pozwalają śledzić układ do czasów rzędu nanosekund. Drugim po-wodem ograniczenia czasowego trwania obliczeń jest kwestia błędów zaokrągleń. Niektóre parametry atomów układu są opisane przy pomocy liczb zmiennoprzecinkowych. Wy-konując obliczenia komputer dokonuje konwersji liczb zmiennoprzecinkowych zapisanych w formacie dziesiętnym do postaci binarnej, która jest przez niego obsługiwana. Każda liczba zmiennoprzecinkowa po konwersji do systemu binarnego ma postać nieskończonego ciągu zerojedynkowego. Ponieważ obliczenia są prowadzone ze skończoną precyzją, w pew-nym momencie ciąg zerojedynkowy, reprezentujący położenie, prędkość czy przyspieszenie danego atomu, ulega obcięciu wprowadzając niedokładności w opisie stanu układu. Im dłuższy czas obliczeń, tym większe są rozbieżności pomiędzy reprezentacją stanu układu, a jego stanem „rzeczywistym”. Niedoskonałość mocy obliczeniowej komputerów wpływa także na konieczność modelowania układów znacznie mniejszych niż te badane w praw-dziwym eksperymencie. Obecnie maksymalna liczba atomów w modelowanym układzie może wynosić kilka milionów. Dla porównania liczba atomów w jednym molu substancji to około 6.02 · 1023, czyli dużo, dużo więcej. Ograniczenie maksymalnej liczby atomów w układzie poważnie ogranicza skalę zjawisk, które mogą być poddane modelowaniu.

Niedokładności w opisie stanu układu są również związane z wyborem konkretnej metody prowadzenia obliczeń. Typowe problemy, z którymi trzeba się zmierzyć, zostaną przedstawione na przykładzie dynamiki molekularnej, ponieważ program obliczeniowy, z którego korzystano, został napisany w myśl założeń tej metody.

Zagadnienie ewolucji układu fizycznego może być modelowane przy zastosowaniu dwóch podejść: deterministycznego i stochastycznego [HEE97]. Metody deterministyczne polega-ją na numerycznym rozwiązywaniu równań ruchu. Zaczynapolega-jąc od dobrze zdefiniowanych punktów początkowych w przestrzeni fazowej, otrzymuje się nowe położenia i pędy

opisu-jące stan cząstek w modelowanym układzie. Metody stochastyczne wykorzystują rachunek prawdopodobieństwa.

Dynamika molekularna jest przedstawicielem metod deterministycznych. Pozwala ona na wyznaczenie trajektorii w przestrzeni fazowej zbioru N atomów, z których każdy pod-lega klasycznym równaniom ruchu Netwona [ALL92]:

d2~ri(t) dt2 = ¹ mi X i<j ~F_ij, i, j = 1, ..., N. (16)

We wzorze 16, ~Fijjest siłą oddziaływania atomów i oraz j, mijest masą i-tego atomu, a~ri(t) jego wektorem położenia. Aby znaleźć siłę oddziaływania pomiędzy atomami, konieczna jest znajomość potencjału ich oddziaływania:

~Fij = −∇Vij(r). (17)

Już w tym momencie można zauważyć, że założony newtonowski charakter ruchu cząstek wyklucza możliwość modelowania przy pomocy dynamiki molekularnej efektów kwantowych. Wszystkie cząstki wchodzące w skład układu są neutralne bez możliwości wzbudzenia czy jonizacji. Zderzenia pomiędzy nimi zachodzą na zasadzie balistycznych kolizji zachowujących całkowitą energię, pęd oraz liczbę cząstek.

Istnieje wiele algorytmów całkujących równanie 16. W dalszej części tego rozdziału zo-stanie zaprezentowany jeden z nich - metoda predictor–corrector rzędu piątego [ALL92], która została zaimplementowana do programu liczącego, z którego korzystano podczas prowadzenia obliczeń. Jednak idea każdego z algorytmów jest podobna. Mianowicie, przy pomocy znanych wielkości opisujących stan układu w chwili „t”, zostają obliczone nowe wartości położeń, prędkości oraz przyspieszeń atomów w czasie „t+δt”. Wprowadzenie kroku czasowego δt oznacza, że problem ciągły zostaje zastąpiony problemem dyskret-nym. Najbardziej bezpośrednim sposobem dyskretyzacji równania 16 jest zastosowanie rozwinięcia w szereg Taylora, np. dla dowolnej zmiennej K [RAL75], [HEE97]:

K(t + δt) = K(t) +ⁿ⁻¹X i=1 (δt)i i! diK(t) dti + Rn, (18) gdzie diK(t)

dti oznacza i-tą pochodną po czasie zmiennej K. Dyskretyzacja problemu powo-duje powstanie błędu aproksymacji Rn, związanego z obcięciem szeregu zdefiniowanego

równaniem 18 na wyrazie rzędu n. Im więcej uwzględni się członów rozwinięcia (im większy rząd metody), tym większa będzie dokładność obliczeń. Równocześnie wydłuży się jednak czas obliczeń.

Bardzo ważne jest, aby właściwie dobrać wartość kroku czasowego. Wybór zbyt du-żego kroku czasowego spowoduje pominięcie, w procesie modelowania ewolucji systemu, wielu stanów pośrednich układu, które istotnie mogłyby wpłynąć na stany następne. Tym samym otrzymany wynik końcowy nie byłby zbyt wiarygodny. Z kolei wybór zbyt małego kroku czasowego sprawi, że moc komputera zostanie użyta do obliczenia wielu praktycznie identycznych stanów pośrednich układu, a otrzymanie każdego z nich zabierze pewien czas. Idealny wybór wartości kroku czasowego to taki, który umożliwia jak najszybsze prowa-dzenie obliczeń, nie powodując przy tym utraty dokładności obliczeń. Należy pamiętać również, że modelowane zjawisko może być charakteryzowane przez więcej niż jedną skalę czasową. Wtedy należy rozważyć wprowadzenie kroku czasowego o zmiennej wartości, automatycznie dopasowującej się do aktualnych potrzeb. Praktyczna wskazówka dotyczą-ca wyboru kroku czasowego stwierdza, że wybór jego wartości powinien być taki, aby fluktuacje wartości energii całkowitej nie przekroczyły kilku procent fluktuacji wartości energii potencjalnej [HEE97]. Krok czasowy rzędu 10−14s jest wystarczającą wartością wg [HEE97] dla układów, których atomy oddziałują potencjałem Lennarda-Jonesa [AZI86].

Poza wyborem rzędu metody oraz wartości kroku czasowego jest jeszcze kilka kwestii, które mogą wpłynąć na szybkość obliczeń nie pogarszając znacząco ich jakości. Jedną z nich jest nieuwzględnianie oddziaływania pomiędzy istotnie oddalonymi od siebie atoma-mi. Okazuje się, że uwzględnianie w obliczeniach oddziaływań zachodzących pomiędzy bardzo odległymi od siebie atomami jest najbardziej czasochłonną częścią każdego kro-ku maszynowego symulacji. Przez krok maszynowy można rozumieć jeden cykl obliczeń „przenoszący” układ ze stanu opisanego w czasie „t”, do stanu w czasie „t+δt”. Dla od-działywań dwuciałowych w celu obliczenia pełnej siły (uwzględniając przyczynki do siły od wszystkich atomów w układzie) działającej na dany atom, trzeba wykonać ∼ N2 ope-racji, ponieważ tyle oddziaływań dwuciałowych trzeba obliczyć w układzie składającym się z N atomów. Większość z nich jest zupełnie niepotrzebna. Dlatego korzystne byłoby pominięcie oddziaływań pomiędzy atomami odległymi o taką odległość, dla której poten-cjał ich oddziaływania jest bliski zeru. Definiuje się odległość obcięcia potenpoten-cjału RCUT (z ang. „cut off”), powyżej której uznaje się, że atomy nie oddziałują. Aby nie sprawdzać za każdym razem, które atomy oddziałują ze sobą, a które nie, dla każdego atomu sporządza

się listę jego najbliższych sąsiadów [ALL92]. Tylko oddziaływania z atomami z tej listy są brane do wyliczenia całkowitej siły działającej na dany atom. Koszt aktualizacji listy dalej jest rzędu N2, ale nie wykonuje się jej w każdym kroku.