1 Obliczenia Przesunięć Chemicznych NMR W Trybie Równoległym Krzysztof Woliński Wydział Chemii UMCS 20-031 Lublin, Pl. M.C. Skłodowskiej 3

2 Plan prezentacji (1) podstawy teoretyczne (2) transformacja kalibracyjna i metoda GIAO (3) tensor magnetycznego ekranowania (4) aspekty numeryczne (całki, CPHF) (5) obliczenia równoległe (model Master-Slave) (6) przykłady obliczeń dla dużych układów

3 Podstawy teoretyczne Zewnętrzne pole magnetyczne B powoduje powstanie w molekule dodatkowego pola magnetycznego. Jego wielkość i orientacja zależą od lokalnej struktury elektronowej. W pozycji jądra N to indukowane pole określone jest przez Gdziejest magnetyczną stałą ekranowania jądra N.

4 Oddziaływanie indukowanego pola z momentem magnetycznym jądra daje wkład do energii układu Stąd tensor magnetycznej stałej ekranowania określony jest jako druga pochodna energii układu w zewnętrznym polu względem tego pola i momentu magnetycznego jądra Podstawą obliczeń jest rachunek zaburzeń (podwójnych)

5 Zewnętrzne pole magnetyczne reprezentowane jest w hamiltonianie układu poprzez potencjał wektorowy A gdzie sumowanie obejmuje wszystkie elektrony. Na potencjał A dla jednego elektronu składa się potencjał pola zewnętrznego oraz potencjał wynikający z momentów magnetycznych jąder G oznacza dowolny punkt odniesienia (gauge origin).

6 Podstawienie potencjału A do hamiltonianu układu pozwala wyizolować człony perturbacyjne liniowe i kwadratowe ze względu na składowe pola oraz momentu magnetycznego jąder Wszystkie operatory zaburzeniowe są sumą odpowiednich operatorów jednoelektronowych.

7 Ponadto konieczna jest oczywiście znajomość funkcji falowej dla problemu niezaburzonego oraz poprawki I-go rzędu względem zewnętrznego pola magnetycznego potrzebne są następujące operatory jednoelektronowe Do wyznaczenia tensora magnetycznego ekranowania jądra

8 Funkcje falowe 0-go i I-go rzędu są rozwiązaniami, odpowiednio równania Schrödingera dla układu izolowanego oraz równania perturbacyjnego I-go rzędu Rozwiązanie powyższych równań wymaga wprowadzenia szeregu przybliżeń: - ruch elektronów w polu nieruchomych jąder - przybliżenie jednoelektronowe (metoda Hartree-Focka) - przybliżenie analityczne orbitali HF (Roothaan)

9 Ostatnie przybliżenie wprowadza pojęcie bazy funkcyjnej. Poszukiwane orbitale HF przedstawiane są w postaci liniowych kombinacji skończonej liczby znanych funkcji i=1,n (zakładamy 2n-elektronowy układ zamknięto-powłokowy) Tak przedstawione orbitale molekularne opisywać mogą zarówno układ izolowany jak również układ zaburzony. W tym ostatnim przypadku współczynniki rozwinięcia C oraz w ogólności funkcje bazy χ zależą od wielkości zaburzenia.

10 Powyższe przybliżenia pozwalają ostatecznie zapisać równania 0-go i I-go rzędu w zwartej formie macierzowej Równania te rozwiązywane są iteracyjnie ponieważ macierz Focka F zależy od poszukiwanych współczynników C.

11

12 Równanie perturbacyjne I-go rzędu można przekształcić do postaci gdzie sumowanie przebiega po orbitalach zajętych i oraz wirtualnych a. Z powodu sprzężonego charakteru ( F(D) ) równanie to nosi nazwę Coupled Perturbed Hartree-Fock (CPHF). Znajomość rozwiązań równań 0-go i I-go rzędu pozwala wyznaczyć interesującą nas poprawkę II-go rzędu do energii układu tzn. tensor magnetycznej stałej ekranowania jądra Gdzie Tr oznacza ślad macierzy.

13 Transformacja kalibracyjna i metoda GIAO Zewnętrzne pole magnetyczne pojawia się w hamiltonianie układu poprzez potencjał wektorowy, zależny od pewnego punktu kalibracji G Jednakże to nie potencjał pola, a jego natężenie określa stan układu w zewnętrznym polu. To natężenie dane jako jest niezależne od wyboru punktu odniesienia (kalibracji) G

14 Tak więc, zmieniając potencjał poprzez wybór innego punktu kalibracji G (na nieskończenie wiele sposobów) otrzymujemy zawsze takie samo natężenie pola B. Jednakże każda zmiana kalibracji G zmienia hamiltonian układu zaburzonego gdyż są w nim człony od tego punktu zależne. Stąd też różne wielkości obliczane z takim hamiltonianem będą zależeć od wyboru punktu kalibracji co jest niefizyczne. (człony te znikają tylko dla ścisłych rozwiązań) Jest to tzw. problem transformacji kalibracyjnej.

15 W przypadku użycia tradycyjnych, niezależnych od zaburzenia baz funkcyjnych (takich jak w SCF) tensor magnetycznej stałej ekranowania jądra zależy od wyboru punktu kalibracji G zarówno w swojej części dia- jak i paramagnetycznej. W części diamagnetycznej punkt G pojawia się w H N (11), reprezentacji macierzowej operatora h N (11) (w danej bazie funkcyjnej)

16 W części paramagnetycznej od punktu kalibracji G zależy macierz gęstości I-go rzędu. Ta zależność pojawia się poprzez operator h (10) który jest częścią macierzy Focka I-go rzędu Człony zależne od G znikają tylko w przypadku zupełnej bazy funkcyjnej (granica Hartree-Focka) co jest w praktyce nieosiągalne. Skończenie wymiarowe bazy funkcyjne używane w obliczeniach kwantowo chemicznych powodują naruszenie niezmienniczości kalibracyjnej. Pierwszą tego manifestacją jest zależność obliczanych wielkości (stała ekranowania, podatność magnetyczna) od położenia molekuły w układzie współrzędnych.

17 Istnieje kilka metod usunięcia niefizycznej zależności magnetycznej stałej ekranowania od wyboru punktu kalibracji. Najbardziej popularną jest metoda Gauge-Including Atomic Orbitals (GIAO). W metodzie GIAO baza funkcyjna zależy od zaburzenia tj. od zewnętrznego pola magnetycznego gdzie A B (R  ) jest wartością potencjału zewnętrznego pola w pozycji centrum orbitalu atomowego   Zależność fukcji GIAO od zewnętrznego pola powoduje pojawienie się dodatkowych członów w wyrażeniach na stałą ekranowania.

18 Te dodatkowe człony zawierają pierwsze pochodne funkcji GIAO względem pola Oraz Następuje pełne kasowanie członów zawierających G i ostatecznie zarówno część dia- jak i paramagnetyczna stają się niezależne od punktu kalibracji.

19 Tensor magnetycznego ekranowania Spektroskopia NMR nierozerwalnie łączy ekranowanie pola magnetycznego z jądrami atomowymi. Jaką rolę odgrywają jądra z niezerowym spinem w eksperymencie NMR ? O efekcie ekranowania a tym samym o wielkości i orientacji lokalnego pola magnetycznego decyduje struktura elektronowa układu. Zjawisko to ma miejsce w każdym punkcie molekuły, a nie tylko w pozycji jej jąder. Jądra z niezerowym spinem są w eksperymencie NMR sondą pozwalającą zmierzyć lokalne pola magnetyczne.

20 Czy można określić magnetyczą stałą ekranowania w dowolnym punkcie w molekule, a nie tylko w pozycji jąder ? w pomiarze NMR - NIE na drodze obliczeń - TAK Tensor magnetycznego ekranowania w dowolnym punkcie R molekuły obliczamy dokładnie tak samo jak w pozycji jakiegoś jądra. W punkcie, w którym nie ma jądra z niezerowym spinem nie można zmierzyć w NMR lokalnego pola magnetycznego bo brakuje probierza. Hipotetycznie taką sondą może być neutron (ładunek zero, spin ½) umieszczony (jakoś) w danym punkcie. Może więc: NMR = Neutron Magnetic Resonance ?

21 Tensor magnetycznego ekranowania w dowolnym punkcie R molekuły można rozpatrywać jako ciągłą funkcję współrzędnych punktu co prowadzi do pojęcia trójwymiarowej powierzchni magnetycznego ekranowania (magnetic shielding surface) Analiza powierzchni ekranowania daje dużo więcej informacji niż znajomość jądrowych stałych ekranowania. Dla przykładu, stałe ekranowania 13C i 15N w HCN i HNC są zupełnie różne : HCN - 13C: 81.27 ppm, 15N: -33.78 ppm HNC - 13C: 10.87 ppm, 15N: 102.82 ppm Jednakże powierzchnie ekranowania są w obu przypadkach prawie takie same

22

23

24 Rozpatrywanie  R jako ciągłej funkcji R pozwala na jej całkowanie. Jeśli zrobimy to po całej przestrzeni R 3 to otrzymamy zintegrowany tensor magnetycznego ekranowania (integrated magnetic shielding) opisujący ekranujące własności całej molekuły. Okazuje się, że ta wielkość jest bezpośrednio związana z podatnością magnetyczną  gdzie  0 jest przenikalnością magnetyczną próżni. Podobne relacje łączą również inne własności magnetyczne:

25 Aspekty numeryczne obliczeń magnetycznych stałych ekranowania Obliczenia magnetycznych stałych ekranowania jąder są pod względem numerycznym dość złożone. Są bardziej skomplikowane niż obliczenia gradientu (sił działających na atomy w molekule) ale dużo prostsze od obliczeń hessianu (stałych siłowych). Obliczenia magnetycznych stałych ekranowania zajmują ok. 1.5-2.0 razy więcej czasu niż rozwiązanie równania Schrödingera dla izolowanej molekuły (SCF). Dla porównania, obliczenie gradientu jest ok. 3.0 razy szybsze niż SCF, a hessianu nawet 10-15 razy wolniejsze.

26 Co trzeba policzyć i ile to kosztuje Należy obliczyć szereg całek 1- i 2-elektronowych oraz ich pierwszych pochodnych względem pola magnetycznego oraz rozwiązać równania CPHF. Krok obliczeniowy koszt 1. 1-el. całki : S (10),T (10),H (10) praktycznie zerowy 2. 1-el. Całki : V (10) minimalny 3. 1-el. Całki H N (11 ), H N (01) mały ale rośnie z N 4. 2-el. Całki (ij|kl) (1) duży 5. 2-el. Całki (ij|kl) (0) i CPHF duży (iteracje)

27 Najbardziej czasochłonnym etapem we wszystkich obliczeniach kwantowo-chemicznych jest wyznaczanie różnego rodzaju całek 2-elektronowych w danej bazie funkcyjnej oraz operacje z nimi związane (konstrukcja macierzy typu Focka, transformacje). Całki 2-elektronowe mają w ogólności następującą postać Całki 2-el. są trudne do liczenia i jest ich bardzo dużo ~ m 4. W dużych molekułach znaczna ich część znika i nie musi być wyznaczana. W praktyce zostaje ~ m 2.5 (m 2 ) ale to wciąż b.dużo.

28 Obliczone całki 2-el. są następnie kontraktowane z elementami macierzy gęstości tworząc macierz typu Focka W ogólności jedna całka 2-el daje wkład do 6-ciu elementów G. W obliczeniach magnetycznych stałych ekranowania mamy : m 4 (m 3 ) całek (ij|kl) (0) oraz 3 macierze G(D (1),g (0) ) ( n razy) 6*m 4 (m 3 ) całek (ij|kl) (1) oraz 3 macierze G(D (0),g (1) ) ( 1 raz ) Te obliczenia są bardzo intensywne numerycznie i czasochłonne. Jak realizować je effektywnie ? W trybie równoległym tzn. z wykorzystaniem N procesorów

29 Obliczenia równoległe Obliczenia równoległe można realizować na kilka sposobów. Ważnym kryterium jest tu liczba procesorów, które mają być lub mogą być w obliczeniach użyte. Możliwe jest użycie nawet tysięcy procesorów (Massively Parallel Computing) (np. NWchem).Takie podejście wymaga specjalnych maszyn, które są oczywiście bardzo drogie (miliony $). Alternatywne i nieporównywalnie tańsze podejście to klaster kilku-kilkudziesięciu (4-32,64, a nawet 256) procesorów PC.

30 W przypadku klasteru z ze stosunkowo niewielką liczbą procesorów odpowiednim modelem organizacji paralelizacji jest model Master –slave W modelu tym na jednym z procesorów (Master) wykonywany jest główny program, w którym następuje rozdzielenie zadań między pozostałe procesory (slaves). Komunikacja między procesorami realizowana jest za pomocą pakietu Parallel Virtual Machine (PVM). W przypadku obliczeń NMR komunikacja taka ma miejsce na drodze master-slave-master Każdy slave wykonuje własny podprogram (taki sam) realizujący zlecone zadania. W naszym przypadku są to obliczenia różnych całek 2-el. oraz konstrukcja odpowiednich macierzy typu Focka.

31 Realizacja obliczeń równoległych odbywa się następująco : Master : 1. Wysyła do wszystkich slaves podstawowe dane (geometria,baza, żądane dokładności itd..) 2. Wysyła do wszystkich slaves odpowiednie macierze gęstości 3. Tworzy listę całek (bloków całek), które mają być obliczane (tysiące-setki tysięcy bloków różnego typu i wielkości) 4. Wysyła pierwsze N bloków całek do N kolejnych slaves 5. Czeka na sygnał od slave ”skończyłem” i wysyła następne bloki do zgłaszających się slaves (first-come first-served) Jeśli bloki się już wyczerpały i nie ma więcej zadań to slave dostaje sygnał ”koniec pracy” co jest dla niego również poleceniem ”odeślij swoją część macierzy Focka” 6. Odbiera od slaves części macierzy Focka i je sumuje

32 Slave : 1. Odbiera od master podstawowe dane 2. Odbiera od master podstawowe macierze gęstości D 3. Tworzy listę całek (bloków całek), które mają być obliczane 4. Odbiera od master numer swojego pierwszego bloku 5. Oblicza całki i kontraktuje je z odpowiednimi elementami macierzy D tworząc lokalną część macierzy Focka F 6. Wysyła do master sygnał ”skończyłem” 7. Odbiera od master numer kolejnego bloku i powtarza 5-6 albo ”koniec pracy” i wysyła do master swoją część macierzy F Slaves rozpoczynają pracę w tym samym czasie i pożądane jest aby ją jednocześnie zakończyły. W tym celu początkowo do slaves wysyłane są duże bloki całek, a w miarę postępu obliczeń coraz mniejsze.

33 Przykłady obliczeń dla dużych układów Wszystkie prezentowane tu obliczenia wykonane zostały przy użyciu pakietu PQS (Parallel Quantum Solutions) na klasterze 11-CPU (PC) : 1 procesor AMD Athlon XP 2800+ (2105 MHz, 2.0GB ) 6 procesorów AMD Athlon XP 2800+ (1802 MHz, 1.5 GB) 4 procesory AMD Athlon XP 3000+ (2105 MHz, 2.0GB ) Procesory połączone są 100Mb/s kablem Ethernet

34 ----------------------------------------------------------------------------------------------------- Czas obliczeń (min) magnetycznych stałych ekranowania jąder Baza funkcyjna 6-311G-dp, klaster 11 cpu ----------------------------------------------------------------------------------------------------- Molekuła wzór atom baza SCF giao cphf h1n NMR TOT ----------------------------------------------------------------------------------------------------- C 70 70 1260 45 25 27 0 55 99 C 54 H 18 72 1080 20 14 23 0 39 49 C 96 H 24 120 1872 70 54 86 2 146 217 Taxol C 47 H 51 NO 14 113 1422 92 74 63 7 146 238 Chlorophyll C 55 H 72 N 4 O 5 Mg 137 1610 100 64 66 11 144 245 Valinomycin C 54 H 90 N 6 O 18 168 1944 160 120 105 19 250 411 Gramicidin A C 99 H 140 N 20 O 17 276 3288 600 467 643 88 1230 1830 ----------------------------------------------------------------------------------------------------- Średnia efektywność obliczeń równoległych : SCF - 9.03 NMR - 9.21 TOT - 9.10 -----------------------------------------------------------------------------------------------------

35 Czas obliczeń (min) magnetycznych stałych ekranowania jąder Baza funkcyjna 6-311G-dp, klaster 11 cpu ----------------------------------------------------------------------------------------------------- Molekuła DFT atom baza SCF giao cphf h1n NMR TOT ----------------------------------------------------------------------------------------------------- Taxol B3LYP 113 1422 126 80 50 7 143 273 WAH 128 83 - 7 96 224 Chlorophyll B3LYP 137 1610 160 69 49 11 135 295 WAH 178 70 - 11 87 265 Valinomycin B3LYP 168 1944 218 133 92 19 257 475 WAH 242 137 - 19 169 411 ----------------------------------------------------------------------------------------------------- Średnia efektywność obliczeń równoległych : B3LYP : SCF - 9.19 NMR - 9.95 TOT - 9.58 WAH : SCF - 9.06 NMR - 10.5 TOT - 9.55 -----------------------------------------------------------------------------------------------------