Ćwiczenie 1:

1

Akademia Górniczo–Hutnicza

Wydział Inżynierii Materiałowej i Ceramiki

Katedra Chemii Krzemianów i Związków Wielkocząsteczkowych

Instrukcja do ćwiczeń laboratoryjnych

Kierunek studiów: Technologia chemiczna Specjalność: Analityka i kontrola jakości Laboratorium: Analiza strukturalna materiałów

Ćwiczenie 1:

Wykorzystanie metod dynamiki molekularnej do

modelowania struktur szkieł boro-krzemianowych i procesów

dyfuzji.

2 2018/2019

1. Cel ćwiczenia

Celem ćwiczenia jest wykonanie symulacji klasycznej dynamiki molekularnej szkieł na bazie SiO2-B2O3-NaO w celu zbadania ich struktury oraz procesu dyfuzji w nich zachodzącej.

2. Wprowadzenie teoretyczne / obowiązujący zakres materiału

a.

Szkła boro-krzemianowe ze względu na dobra zdolność szkłotwórczą oraz dużą łatwość w przyjmowaniu różnego rodzaju odpadów do swojej struktury, były pierwszymi szkłami zastosowanymi do immobilizacji odpadów toksycznych i radioaktywnych.

Głównymi składnikami tych szkieł są: SiO2, Al2O3, B2O3 oraz Na2O. Atomy Si, Al i B posiadają koordynację 4 względem tlenu i tworzą z nim tetraedry. Tetraedry te łączą się narożami przez tzw. tleny mostkowe BO, powoduje to powstanie ciągłej sieci przestrzennej szkła.

Alkalia, metale przejściowe, ciężkie jony znajdujące się w lukach sieci nazywamy modyfikatorami, powodują one zrywanie połączeń między tetraedrami i w efekcie powstanie tlenów niemostkowych NBO. Skutkiem tego jest częściowe osłabienie i depolimeryzację struktury szkła. Modyfikatory posiadają zazwyczaj wyższe liczby koordynacyjne niż 5. Glin również może posiadać liczbę koordynacyjną równą 6, i wtedy pełni on funkcję modyfikatora. Bor natomiast występuje również w koordynacji 3.

3 Rys.1. Szkło boro-krzemianowe 47,2 SiO2 14,9 B2O3 4,4 Al2O3 33,5 Na2O mol%

b. Dynamika molekularna

Jest to jedna z technik komputerowych, w której bada się ewolucję czasową układu oddziaływujących ze sobą atomów. Polega ona na rozwiązywaniu klasycznych równań ruchu. Na każdy atom i z N rozpatrywanych działa siła, która w idealnym przypadku jest sumą wszystkich sił działających na dany atom w układzie:

(1)

-siła będąca sumą wektorową wszystkich sił działających na i-ty atom -masa i-tego atomu

-przyspieszenie działające na i-ty atom

Wektor przemieszczenia i-tego atomu można wyliczyć wtedy z zależności:

(2)

Gdzie:

-wektor przemieszczenia i-tego atomu -czas

 4 1

2

3 4

5 6

7

1-Teraedr [SiO4] 2-Teraedr [AlO₄]^- 3-Teraedr [BO4]^- 4- B o koordynacji 3

5-Atom Na 6-tlen mostkowy BO

7-tlen

niemostkowy NBO

4 Główny problem symulacji komputerowych to dobranie odpowiedniego modelu matematycznego opisującego rozważany problem. W przypadku dynamiki molekularnej jest to właściwy opis sił działających na poszczególny atom układu. W celu wyznaczenia siły wykorzystuje się funkcję , która zależy od położenia atomu i reprezentuje energię potencjalną danego atomu.

(3)

Oznacza to, że jeśli znamy potencjał opisujący interesujące nas odziaływania i początkowy stan układu możemy wyznaczyć położenie atomów w dowolnej chwili czasowej.

Przykładowym potencjałem opisującym odziaływania między jonami lub atomami (traktowane modelowo jako ładunki) może potencjał Buckingham’a o następującym wzorze:

(4)

Gdzie:

-to ładunki jonów w parze i-j

, , – to parametry dla danej pary atomów wyznaczone eksperymentalnie lub pochodzące z obliczeń teoretycznych.

c. Wprowadzenie do LAMMPS’a

Lammps jest jednym z programów umożliwiających symulacje z wykorzystaniem klasycznej dynamiki molekularnej MD. Jest on w pełni darmowy. Program jest dostępny na stronie: http://www.lammps.sandia.gov.

W programie tym operuje się z poziomu wiesza polecania systemu Windows. Aby wywołać wiersz opleceń wystarczy wpisać w polu wyszukiwania w menu Start cmd.

Aby program obliczył nasz plik wsadowy wchodzimy do folderu, gdzie się znajduje, wpisując polecenie w wierszy poleceń cd ścieżka do folderu np. cd Desktop\Laboratoria materiały reaktorowe. Aby program rozpoczął obliczenia możemy wykorzystać następujące komendy:

a. Lmp_serial.exe -in „nazwa pliku”

Jest to najprostsza komenda, niestety program wykonuje operacje jednowątkowo.

„nazwa pliku” – nasza nazwa pliku wsadowego

b. Lmp_serial -in „nazawa pliku” -pk omp „X” -sf omp

Wykorzystujemy tą komodę, gdy chcemy korzystać z więcej niż jednego rdzenia naszego procesora. Nie wymaga ona dodatkowego oprogramowania. „X”- to liczba wątków, które

5 może wykonywać na raz nasz procesor np. Lmp_serial -in in-3 -pk omp 4 -sf omp - program wczytuje plik o nazwie in-3 (nazwa musi być wpisana razem z rozszerzeniem jeśli plik takowe posiada np. nazawa.txt) i wykonuje nad nim obliczenia na 4 wątkach lub rdzeniach procesora.

Rys.2. Wiersz polecania, przykłady wpisanych komend.

Aby wykonać obliczenia potrzebny jest plik wsadowy, w którym zawarte są wszystkie informacje niezbędne do przeprowadzenia symulacji nad danym układem atomów. Można go napisać za pomocą najprostszych edytorów tekstu np. notatnik lub wordpad.

Przykład pliku wsadowego wraz z objaśnieniem komend:

#SZKLO 0.472 SiO2 14.9 B2O 4.4 Al2O3 33.5 Na2O POTENCJAL BMH - komentarze zaczynamy od znaku „#” obowiązuje on do końca linijki.

units metal -definiuje styl jednostek, jakie używa program np. „metal” (w tym stylu przykładowo energia jest w (eV), a odległość w (Å)).

atom_style charge -Definiujemy styl atomów „charge” atomy posiadają oprócz masy, ładunki.

dimension 3 -definiuje ilość wymiarów układu, w tym przypadku układ jest trójwymiarowy boundary p p p -warunki na brzegach układu, „p” oznacza warunki periodyczne. Tutaj w kierunku x y z zadane są warunki periodyczne.

read_data Poczatek.txt -wczytuje plik, w którym znajdują się informacje o początkowych położeniach atomów (w zależności od stylu również ich ładunkach, masie i typie).

restart 1000 szklo-1.restart -zapisuje plik restartu co 1000 kroków czasowych

group Si type 1 -nazywamy grupy atomów np. „Si” (nasza własna nazwa) przypisujemy im typ np. „1” (też nasz własny) musi być zgodny z typem użytym w pliku zawierającym informacje o początkowym położeniu atomów.

6 group B type 2

group Al type 3

group Na type 4

group O type 5

#PARAMETRY POTENCJALU

pair_style buck/coul/long 12 -Wybór używanego potencjału, w tym przypadku wybrano potencjał Buckingham’a, „12” oznacza zasięg, gdzie interakcje między cząstkami są obliczane bezpośrednio (dotyczy potencjału i odziaływań kulombowskich) kspace_style pppm 10e-6 -oznacza pożądany względny błąd w obliczeniach wartości sił, dotyczy on kulombowskich odziaływań dalekiego zasięgu dalej niż „12”, obliczane są one za pomocą metody pppm.

pair_coeff * * 0 1 0 -przypisuje wszystkim parą atomów zadane

parametry potencjału

pair_coeff 1 5 50306.24729 0.161 0.0 -dla pary atomów złożonej z typu 1 i 5 nadawane są odpowiednie parametry potencjału

velocity all create 10000 123456 -wszystkim atomom nadawane są prędkości tak aby ich średnia energia kinetyczna była tożsama z temperaturą 10000K liczba „123456” jest niezbędna do generatora liczb losowych i może być dowolna

#ALGORYTM SYMULACJI

dump id all xyz 1000 szklo.xyz -tworzy plik o rozszerzeniu xyz, który zawiera informację o położeniu „all” wszystkich atomów, „id” to nazwa polecenia (może być dowolna), „1000” oznacza, że zapisuje w tym pliku informacje co 1000 kroków czasowych. Plik ten można wykorzystać do wizualizacji układu np. za pomocą programu OVITO (jest on darmowy dostępny na stronie: http://wwwovito.org).

thermo 100 -termodynamiczne informacje o układzie np. temperatura, energia, objętość zapisywane są w pliku log.lammps oraz pokazywane w wierszu polecania co 100 kroków czasowych.

timestep 0.001 -długość kroku czasowego, w tym przypadku 0,001ps (zależy od stylu jednostek) fix 1 all nvt temp 10000 1000 1 -komenda fix pozwala zadać przebieg algorytmu

symulacji. „1” -nazwa fix’a (nasza nazwa), „all” -dla wszystkich grup atomów, „nvt”- algorytm, gdzie objętość układu jest stała, zmienia się ciśnienie. „10000” -to temperatura początkowa, „1000” -temperatura końcowa, „1” -czas relaksacji temperatury w układzie w jednostkach czasu.

run 20000 -określa ilość kroków czasowych w jakiej są wykonywane fix’y, które znajdują się przed tą komendą.

Dokładna instrukcja do programu wraz z listą komend znajduje się na stronie programu, w zakładce Manual.

7 Rys.3. Wiersz polecania, w czasie pracy programu pokazywane są parametry

termodynamiczne układu co zadaną ilość kroków czasowych.

3. Przebieg ćwiczenia

a. Wykonanie symulacji dla szkła o składzie: 50 SiO2 (50-x) B2O3 x Na2O mol%

o Symulację należy wykonać według algorytmu w Tabela 1.

Układ na początku symulacji składa się z X atomów umieszczonych w sześciennym pudełku o periodycznych warunkach brzegowych i krawędzi długości wyliczonej z gęstości podanej przez prowadzącego. Długość kroku czasowego w symulacji zadana jest na 1fs.

Tabela 1. Algorytm symulacji

kroki czasowe Temperatura [K] Stałe ciśnienie [bar]

Od 0 do 10000 10000 Brak (Stała objętość)

Od 1001 do 20000 Spada od 10000 do 1200 Od 20001 do 30000 1200

Od 30001 do 31000 Spada od 1200 do 1000 1 Od 31001 do 41000 1000

Od 41001 do 61000 Spada od 1000 do 300 Od 61001 do 71000 300

Odziaływania między atomami w symulacji opisywane są przez potencjał Buckingham’a. Parametry potencjału podane są w Tabeli 2.

8 Tabela 2. Parametry potencjału

Para atomów Ładunek efektywny (e)

(eV) (Å) (eV Å⁶)

Si-O 1,890 50306,24 0,161 46,30

B-O 1,4175 15176,81 0,150 9,08

Al-O 1,4175 18538,47 0,172 34,58

Na-O 0,4725 120304,15 0,170 0,00

O-O -0,945 9022,82 0,265 85,09

o Z uzyskanych wyników symulacji należy sporządzić wykresy funkcji RDF i liczby koordynacyjnej w zależności od odległości dla podanych par atomów.

Funkcja RDF (funkcja rozkładu radialnego) jest to prawdopodobieństwo znalezienia atomu w odległości z przedziału (r, r + dr) od danego atomu odniesienia.

Wykresy powinny być tak wykonane, aby położenie pierwszego maksimum funkcji RDF na osi x było dobrze widoczne. Odległość tego maksimum określa średnią długość wiązania w parze atomów. W przypadku wykresu LK w funkcji odległości powinno być dobrze widoczne pierwsze przegięcie funkcji (schodek).

o Dla pierwszego maksimum funkcji RDF odczytać długość wiązania pary atomów.

o Następnie odczytać promień odcięcia, jest to położenie na osi x minimum funkcji RDF zaraz po pierwszym maksimum. Dla odczytanego promienia odcięcia należy odczytać wartość LK (liczby koordynacyjnej) z jej wykresy dla danej pary.

Odczytane liczby z wykresów zamieścić w tabeli 3.

Tabela 3. Wyniki

Para atomów Długość wiązania [Å] Promień odcięcia [Å] LK [-]

Np. Si-O

b. Dyfuzja Na w szkle.

Dla wcześniej wysymulowanego szkła wyznaczyć średnie kwadratowe przemieszczenie (MSD) dla sodu w temperaturach zadanych przez prowadzącego.

(5)

Program Lammps pozwala na wyznaczenie MSD w funkcji czasu (kroku czasowego):

(6) Wtedy współczynnik a=6D. Należy przeliczyć współczynnik dyfuzji na [cm²/s].

9 Następnie obliczyć energię aktywacji z zależności Arrheniusa:

(7) R – stała gazowa

Wykorzystać w tym celu tak samo jak w przypadku współczynnika dyfuzji regresję liniową.

(8)

Na osi y wykresu powinna znajdować się wartość lnD, na osi x: 1/T.

4. Opracowanie wyników

Sprawozdanie powinno zawierać: wykresy, tabelę 3, wykresy zawierające proste regresji liniowej dla MSD i współczynnika aktywacji, oraz obliczone wartości , D i krótkie wnioski.