Minimalizacja energii

W dziedzinie chemii obliczeniowej minimalizacja energii (zwana także optymalizacją energii , minimalizacją geometrii lub optymalizacją geometrii ) to proces znajdowania układu w przestrzeni zbioru atomów, w którym, zgodnie z pewnym obliczeniowym modelem wiązania chemicznego, sieć między -siła atomowa działająca na każdy atom jest akceptowalnie bliska zeru, a pozycja na powierzchni energii potencjalnej (PES) jest punktem stacjonarnym (opisanym dalej). Zbiorem atomów może być pojedyncza cząsteczka , jon , faza skondensowana , stan przejściowy , a nawet zbiór któregokolwiek z nich. Modelem obliczeniowym wiązań chemicznych może być na przykład mechanika kwantowa.

Na przykład, optymalizując geometrię cząsteczki wody , dąży się do uzyskania długości wiązań wodór-tlen i kąta wiązania wodór-tlen-wodór, które minimalizują siły, które w innym przypadku przyciągałyby atomy do siebie lub odpychały je od siebie.

Motywacją do przeprowadzenia optymalizacji geometrii jest fizyczne znaczenie otrzymanej struktury: zoptymalizowane struktury często odpowiadają substancji występującej w przyrodzie, a geometria takiej struktury może być wykorzystana w różnych badaniach eksperymentalnych i teoretycznych w dziedzinach budowy chemicznej , termodynamiki , kinetyki chemicznej , spektroskopii i innych.

Zazwyczaj, choć nie zawsze, proces ma na celu znalezienie geometrii określonego układu atomów, który reprezentuje lokalne lub globalne minimum energii. Zamiast poszukiwania globalnego minimum energetycznego, pożądana może być optymalizacja do stanu przejściowego , czyli punktu siodłowego na powierzchni energii potencjalnej. Ponadto niektóre współrzędne (takie jak długość wiązania chemicznego) mogą zostać ustalone podczas optymalizacji.

Geometria molekularna i interpretacja matematyczna

Geometria zbioru atomów może być opisana wektorem pozycji atomów. Może to być zbiór współrzędnych kartezjańskich atomów lub, w przypadku cząsteczek, mogą to być tak zwane współrzędne wewnętrzne utworzone z zestawu długości wiązań, kątów wiązań i kątów dwuściennych.

Mając zbiór atomów i wektor $r$ , opisujący położenie atomów, można wprowadzić pojęcie energii jako funkcji położeń $E (r)$ . Optymalizacja geometrii jest więc matematycznym problemem optymalizacyjnym , w którym należy znaleźć wartość $r$ , dla której $E (r)$ jest lokalnym minimum , czyli pochodną energii względem położenia atomów, $\partial mi / \partial r$ , jest wektorem zerowym i macierzą drugiej pochodnej układu ${\ Displaystyle {\ rozpocząć {pmatrix} {\ Frac {\ częściowe ^ {2} E}{\partial r_{i}\,\partial r_{j}}}\end{pmatrix}}_{ij}} , znana również$ jako macierz Hessego , która opisuje krzywiznę PES w $r$ , ma wszystkie dodatnie wartości własne (jest dodatnio określone ).

Szczególnym przypadkiem optymalizacji geometrii jest poszukiwanie geometrii stanu przejściowego ; jest to omówione poniżej.

Model obliczeniowy, który zapewnia przybliżoną wartość $E (r)$ mógłby opierać się na mechanice kwantowej (z wykorzystaniem teorii funkcjonału gęstości lub metod półempirycznych ), polach siłowych lub ich kombinacji w przypadku QM/MM . Korzystając z tego modelu obliczeniowego i wstępnego odgadnięcia (lub ansatz ) prawidłowej geometrii, stosuje się iteracyjną procedurę optymalizacji, na przykład:

oblicz siłę działającą na każdy atom (czyli $-\partial E /\partial r$ )
jeśli siła jest mniejsza niż pewien próg, zakończ
w przeciwnym razie przesuń atomy o obliczony krok $∆ r$ , który według przewidywań zmniejszy siłę
powtórzyć od początku

Praktyczne aspekty optymalizacji

Jak opisano powyżej, niektóre metody, takie jak mechanika kwantowa, mogą być użyte do obliczenia energii $E (r)$ , gradientu PES, to znaczy pochodnej energii względem położenia atomów, $\partial E / \partial r$ i druga macierz pochodna układu, $\partial\partial E /\partial r i \partial r j$ , znana również jako macierz Hessego , która opisuje krzywiznę PES w $r$ .

Algorytm optymalizacji może wykorzystywać niektóre lub wszystkie $E (r)$ , $\partial E /\partial r$ i $\partial\partial E /\partial r i \partial r j$ , aby spróbować zminimalizować siły i teoretycznie może to być dowolna metoda, taka jak opadanie gradientu, sprzężenie gradientu lub metody Newtona, ale w praktyce algorytmy wykorzystujące wiedzę o krzywiźnie PES, czyli macierzy Hessego, okazują się lepsze. Jednak w przypadku większości systemów o znaczeniu praktycznym obliczenie macierzy drugiej pochodnej może być zbyt kosztowne i jest szacowane na podstawie kolejnych wartości gradientu, co jest typowe w optymalizacji quasi- niutonowskiej .

Wybór układu współrzędnych może mieć kluczowe znaczenie dla przeprowadzenia udanej optymalizacji. Na przykład współrzędne kartezjańskie są zbędne, ponieważ nieliniowa cząsteczka z $N$ ma $3 N -6$ wibracyjnych stopni swobody , podczas gdy zbiór współrzędnych kartezjańskich ma $3 N$ wymiarów. Ponadto współrzędne kartezjańskie są silnie skorelowane, to znaczy macierz Hesji ma wiele wyrazów niediagonalnych, które nie są bliskie zeru. Może to prowadzić do problemów numerycznych w optymalizacji, ponieważ np. trudno jest uzyskać dobre przybliżenie do macierzy Hessego, a jej dokładne obliczenie jest zbyt kosztowne obliczeniowo. Jednak w przypadku, gdy energia jest wyrażana za pomocą standardowych pól sił, opracowano wydajne obliczeniowo metody umożliwiające analityczne wyprowadzenie macierzy Hessego we współrzędnych kartezjańskich, przy jednoczesnym zachowaniu złożoności obliczeniowej tego samego rzędu, co w przypadku obliczeń gradientowych. Współrzędne wewnętrzne są zwykle mniej skorelowane, ale są trudniejsze do skonfigurowania i może być trudne opisanie niektórych systemów, takich jak te z symetrią lub dużymi fazami skondensowanymi. Wiele nowoczesnych pakietów oprogramowania do chemii obliczeniowej zawiera automatyczne procedury automatycznego generowania rozsądnych układów współrzędnych do optymalizacji.

Stopień ograniczenia swobody

Niektóre stopnie swobody można wyeliminować z optymalizacji, na przykład pozycje atomów lub długości i kąty wiązań mogą mieć ustalone wartości. Czasami są one określane jako zamrożone stopnie swobody.

Rysunek 1 przedstawia optymalizację geometrii atomów w nanorurce węglowej w obecności zewnętrznego pola elektrostatycznego. W tej optymalizacji pozycje atomów po lewej stronie są zamrożone. Ich interakcje z innymi atomami w układzie są nadal obliczane, ale zapobiega się zmianie położenia atomów podczas optymalizacji.

Rysunek 1: Ugięcia elektrostatyczne nanorurki węglowej w polu elektrycznym .

Optymalizacja stanu przejściowego

stanu przejściowego można określić, wyszukując punkty siodłowe na PES interesujących nas związków chemicznych. Punkt siodłowy pierwszego rzędu to pozycja na PES odpowiadająca minimum we wszystkich kierunkach z wyjątkiem jednego; punkt siodłowy drugiego rzędu to minimum we wszystkich kierunkach z wyjątkiem dwóch i tak dalej. Zdefiniowany matematycznie n- tego rzędu charakteryzuje się następująco: $\partial E /\partial r = 0$ a macierz Hessego $\partial\partial E /\partial r i \partial r j$ , ma dokładnie n ujemnych wartości własnych.

Algorytmy lokalizowania geometrii stanu przejściowego dzielą się na dwie główne kategorie: metody lokalne i metody półglobalne. Metody lokalne są odpowiednie, gdy punkt startowy optymalizacji jest bardzo zbliżony do rzeczywistego stanu przejściowego ( bardzo blisko zostanie zdefiniowany wkrótce), a metody półglobalne znajdują zastosowanie, gdy poszukuje się lokalizacji stanu przejściowego przy bardzo niewielkiej wiedzy a priori o jego geometria. Niektóre metody, takie jak metoda Dimera (patrz poniżej), należą do obu kategorii.

Wyszukiwania lokalne

Tak zwana optymalizacja lokalna wymaga wstępnego odgadnięcia stanu przejściowego, który jest bardzo zbliżony do rzeczywistego stanu przejściowego. Bardzo blisko zwykle oznacza, że początkowe przypuszczenie musi mieć odpowiednią macierz Hessian z jedną ujemną wartością własną lub ujemna wartość własna odpowiadająca współrzędnej reakcji musi być większa niż inne ujemne wartości własne. Ponadto wektor własny o najbardziej ujemnej wartości własnej musi odpowiadać współrzędnej reakcji, to znaczy musi reprezentować transformację geometryczną związaną z procesem, którego stan przejściowy jest poszukiwany.

Biorąc pod uwagę powyższe warunki wstępne, lokalny algorytm optymalizacji może następnie poruszać się „w górę” wzdłuż wektora własnego o najbardziej ujemnej wartości własnej i „w dół” wzdłuż wszystkich innych stopni swobody, używając czegoś podobnego do metody quasi-Newtona.

Metoda Dimera

Metodę dimerów można wykorzystać do znalezienia możliwych stanów przejściowych bez znajomości ostatecznej struktury lub do udoskonalenia dobrego odgadnięcia struktury przejściowej. „Dimer” jest tworzony przez dwa obrazy bardzo blisko siebie na PES. Metoda polega na przesunięciu dimeru w górę z pozycji początkowej podczas obracania dimeru w celu znalezienia kierunku najniższej krzywizny (ostatecznie ujemnej).

Aktywacja Techniki Relaksacyjnej (ART)

Technika Relaksacji Aktywacji (ART) jest również otwartą metodą znajdowania nowych stanów przejściowych lub udoskonalania znanych punktów siodłowych w PES. Metoda podąża w kierunku najniższej ujemnej krzywizny (obliczonej za pomocą algorytmu Lanczosa ) na PES, aby osiągnąć punkt siodłowy, rozluźniając się w prostopadłej hiperpłaszczyźnie pomiędzy każdym „skokiem” (aktywacją) w tym kierunku.

Metody łańcucha stanów

Metody łańcuchów stanów można wykorzystać do znalezienia przybliżonej geometrii stanu przejściowego na podstawie geometrii reagenta i produktu. Wygenerowana przybliżona geometria może następnie służyć jako punkt wyjścia do udoskonalenia poprzez wyszukiwanie lokalne, które zostało opisane powyżej.

Metody łańcuchów stanów wykorzystują szereg wektorów, czyli punktów na PES, łączących reagent i produkt interesującej nas reakcji, $r reagenta$ i $r produktu$ , dyskretyzując w ten sposób ścieżkę reakcji. Bardzo często punkty te nazywane są kulkami ze względu na analogię do zestawu kulek połączonych sznurkami lub sprężynami, które łączą reagent i produkty. Seria kulek jest często początkowo tworzona przez interpolację między $r reagentem$ a $r produktem$ , na przykład dla serii kulek $N + 1$ , kulka $i$ może być dana przez

${\ Displaystyle \ mathbf {r} _ {i} = {\ Frac {i} {N }}\mathbf {r} _{\mathrm {iloczyn}}+\left(1-{\frac {i}{N}}\right)\mathbf {r} _{\mathrm {reagent} }}$

gdzie $ja \in 0, 1, ..., N$ . Każda z kulek $r i$ ma energię $E (r i)$ i siły $-\partial E /\partial r i$ są one traktowane w procesie optymalizacji z ograniczeniami, który ma na celu uzyskanie możliwie najdokładniejszego odwzorowania ścieżki reakcji. Aby to osiągnąć, należy zastosować ograniczenia dotyczące odstępów, tak aby każda kulka $ri nie$ została po prostu zoptymalizowana do geometrii reagenta i produktu.

$ri osiągane$ przez rzutowanie składowych siły na każdą stopkę lub alternatywnie ruch każdej stopki podczas optymalizacji, które są styczne do ścieżki reakcji. Na przykład, jeśli dla wygody zdefiniowano, że $g i = \partial E /\partial r i$ , to gradient energii na każdej kulce pomniejszony o składową gradientu energii, która jest styczna do ścieżki reakcji, jest określony wzorem

${\ Displaystyle \ mathbf {g} _ {i} ^ {\ perp} = \ mathbf {g } _{i}-\mathbf {\tau} _{i}(\mathbf {\tau} _{i}\cdot \mathbf {g} _{i})=\left(I-\mathbf {\tau } _{i}\mathbf {\tau} _{i}^{T}\right)\mathbf {g} _{i}}$

gdzie $I$ jest macierzą tożsamości, a $τ i$ jest wektorem jednostkowym reprezentującym styczną ścieżkę reakcji w $r i$ . Dzięki projekcji składowych gradientu energii lub kroku optymalizacji, które są równoległe do ścieżki reakcji, algorytm optymalizacji znacznie zmniejsza tendencję każdej kulki do bezpośredniej optymalizacji do minimum.

Tranzyt synchroniczny

Najprostszą metodą łańcucha stanów jest metoda liniowego tranzytu synchronicznego (LST). Działa poprzez pobieranie interpolowanych punktów między geometrią reagentów i produktów i wybieranie tego o najwyższej energii do późniejszego udoskonalania poprzez wyszukiwanie lokalne. Metoda kwadratowego tranzytu synchronicznego (QST) rozszerza LST, umożliwiając paraboliczną ścieżkę reakcji, z optymalizacją punktu o najwyższej energii prostopadle do paraboli.

Popychana elastyczna opaska

W metodzie Nudged Elastic Band (NEB) kulki wzdłuż ścieżki reakcji symulują siły sprężystości oprócz sił chemicznych $-\partial E /\partial r i$ , aby optymalizator utrzymywał ograniczenie odstępu. Konkretnie, siła $f i$ działająca na każdy punkt i jest wyrażona wzorem

${\ Displaystyle \ mathbf {f} _ {i} = \ mathbf {f} _ {i} ^ {\ równolegle} - \ mathbf {g} _ {i} ^ {\ sprawca }}$

Gdzie

${\ Displaystyle \ mathbf {f} _ {i} ^ {\ równoległy} =k\left[\left(\left(\mathbf {r} _{i+1}-\mathbf {r} _{i}\right)-\left(\mathbf {r} _{i}-\ mathbf {r} _{i-1}\right)\right)\cdot \tau _{i}\right]\tau _{i}}$

jest siłą sprężystości równoległą do ścieżki w każdym punkcie $r i$ ( k jest stałą sprężystości, a $τ i$ , jak poprzednio, jest wektorem jednostkowym reprezentującym ścieżkę reakcji styczną w $r i$ ).

W tradycyjnej implementacji punkt o najwyższej energii jest używany do późniejszego uściślenia w wyszukiwaniu lokalnym. Istnieje wiele odmian metody NEB, w tym obraz wspinaczkowy NEB, w którym punkt o najwyższej energii jest wypychany w górę podczas procedury optymalizacji, aby (miejmy nadzieję) uzyskać geometrię jeszcze bliższą geometrii stanu przejściowego . Istnieją również rozszerzenia obejmujące regresję procesu Gaussa w celu zmniejszenia liczby ocen. W przypadku układów o geometrii nieeuklidesowej (R^2), takich jak układy magnetyczne, metoda jest modyfikowana do geodezyjnego podsuwania elastycznego pasma.

Metoda łańcuchowa

Metoda łańcuchowa wykorzystuje splajny łączące punkty, $r i$ , do pomiaru i wymuszenia ograniczeń odległości między punktami oraz do obliczenia stycznej w każdym punkcie. W każdym kroku procedury optymalizacyjnej punkty mogą być przesuwane zgodnie z działającą na nie siłą prostopadłą do ścieżki, a następnie, jeśli warunek równej odległości między punktami nie jest już spełniony, punkty można ponownie rozmieścić za pomocą splajnu reprezentacja ścieżki do generowania nowych wektorów z wymaganymi odstępami.

Odmiany metody sznurkowej obejmują metodę rosnącej struny, w której odgadywanie ścieżki jest narastane od punktów końcowych (tj. Reagenta i produktów) w miarę postępu optymalizacji.

Porównanie z innymi technikami

Optymalizacja geometrii zasadniczo różni się od symulacji dynamiki molekularnej . Ten ostatni symuluje ruch cząsteczek w czasie, pod wpływem temperatury, sił chemicznych, prędkości początkowych, ruchów Browna rozpuszczalnika i tak dalej, poprzez zastosowanie praw ruchu Newtona . Oznacza to, że obliczane trajektorie atomów mają pewne znaczenie fizyczne. Z kolei optymalizacja geometrii nie tworzy „trajektorii” o jakimkolwiek znaczeniu fizycznym – dotyczy minimalizacji sił działających na każdy atom w zbiorze atomów, a ścieżka, przez którą to osiąga, nie ma znaczenia. Różne algorytmy optymalizacji mogą dać ten sam wynik dla minimalnej struktury energetycznej, ale dochodzą do tego inną drogą.

Zobacz też

Graf złożony z ograniczeń
Cięcia grafów w widzeniu komputerowym – aparat do rozwiązywania problemów widzenia komputerowego , które można sformułować w kategoriach minimalizacji energii
Zasady energetyczne w mechanice budowli

Linki zewnętrzne

Przepisy liczbowe w Fortranie 77

Dodatkowe referencje

Payne i in., „Techniki iteracyjnej minimalizacji do obliczeń energii całkowitej ab initio: dynamika molekularna i gradienty koniugatów”, Recenzje współczesnej fizyki 64 (4), s. 1045–1097. (1992) (streszczenie)
Stich i in., „ minimalizacja gradientu sprzężenia funkcjonału energii: nowa metoda obliczania struktury elektronowej ”, Physical Review B 39 (8), s. 4997–5004, (1989)
Chadi, „ Podejście minimalizacji energii do geometrii atomowej powierzchni półprzewodnikowych ”, Physical Review Letters 41 (15), s. 1062–1065 (1978)