Punkty krytyczne na hiperpowierzchni energii potencjalnej: minima, maksima i punkty siodłowe

Literatura:

1. A. Gołębiewski: Elementy mechaniki i chemii kwantowej, PWN, rozdz. 7. cz. 7.1 - 7.3.

2. W. Krysicki, L. Włodarski: Analiza matematyczna w zadaniach, część I, PWN, rozdz. 10, część II, rozdz. 2, § 1.9, § 1.12, § 1.16.

3. L. Maurin, M. Mączyński, T. Traczyk, Matematyka, podręcznik dla studentów wydziałów chemicznych, PWN, rozdz. 9.

Przydatny będzie program gnuplot, przy pomocy którgo można rysować wykresy funkcji. Program jest darmowy i można go zainstalować zarówno pod systmem Windows jak i Linux.

Wstęp

Hiperpowierzchnie energii potencjalnej.

Zachowanie cząsteczki lub układu cząsteczek określa zmiana jej energii potencjalnej w zależności od geometrii. Z powodu dużej różnicy mas jąder atomowych i elektronów, położenia jąder atomowych określają geometrię cząsteczki. Elektrony można traktować jako obiekty rozmyte, których stan dostosowuje się natychmiast do zmienionej geometrii jąder; jest to treścią przybliżenia adiabatycznego i dalej przybliżenia Borna-Oppenheimera (poz. 1 literatury). Energię w funkcji współrzędnych jąder nazywamy hiperopowierzchnią energii potencjalnej. Te współrzędne mogą być współrzędnymi kartezjańskimi lub wewnętrznymi. Analizę hiperpowierzchni energii potencjalnej przeprowadza się na ogół tylko względem tych współrzędnych, które ulegają zmianie w toku analizowanego procesu. Przykładowo, w badaniu powierzchni energii potencjalnej reakcji izomeryzacji cyjanowodoru (H-CN) do izocyjanowodoru (CN-H) wystarczy użyć jako zmiennych odległości H$\cdots$C i H$\cdots$N. Odległość C$\cdots$N odpowiada zawsze długości wiązania potrójnego między atomem węgla i atomem tlenu i jest w toku całego procesu efektywnie stała.

Hiperpowierzchnie energii potencjalnej odpowiadają energetyce reakcji chemicznych, przemian konformacyjnych (jeżeli rozpatrujemy jedną cząsteczkę i nie zmieniają się wiązania chemiczne) lub oddziaływaniom między cząsteczkami lub atomami/jonami zamkniętopowłokowymi. W każdym przypadku interesuje nas znalezienie punktów stabilnych, odpowiadających minimom energii potencjalnej, oraz stanów przejściowych, leżących na najmniej kosztownej energetycznie drodze od jednego minimum do drugiego (np. od substratów do produktów reakcji lub od jednej konformacji do drugiej). Wspólną ich cechą jest to, że są to tzw. punkty krytyczne lub punkty stacjonarne na hiperpowierzchni energii potencjalnej, w których na układ nie działają żadne siły. Punkty te znajduje się oraz określa ich charakter analizując pochodne i drugie pochodne energii jako funkcji współrzędnych.

Pochodna funkcji jednej zmiennej oraz jej interpretacja geometryczna.

Rozważmy na początek funkcję jednej zmiennej, np. znany z kursu chemii organicznej wykres energii butanu względem kąta obrotu wokół wiązania C-C. Pochodna funkcji jest definiowana jako granica ciągu ilorazów przyrostu wartości funkcji po zwiększeniu zmiennej przez przyrost zmiennej. Te ilorazy nazywają się ilorazami różnicowymi.

$$f^\prime(x) = \frac{\mathbf{d}f}{\mathbf{d}x}=\lim\limits_{h\righ... ...(x+h)-f(x)}{h} = \lim\limits_{\Delta x\rightarrow 0}\frac{\Delta f}{\Delta x}$$

Poniższy rysunek przedstawia ilutrację definicji pochodnej oraz jej interpretację geometryczną. Prowadzimy ciąg siecznych przez punkty $(x,f(x))$ $(x+h,f(x+h))$, zmniejszając stopniowo $h$. Sieczne zbiegają się do stycznej do wykresu funkcji w punkcie $(x,f(x))$. Ponieważ każda sieczna wraz z odcinkami $\overline{(x,x+h)}$ i $\overline{(f(x),f(x+h))}$ tworzy trójkąt prostokątny, iloraz różnicowy jest tangensem kąta tworzonego przez sieczną z odcinkiem równoległym do osi $x$ a przez to tangensem kąta nachylenia siecznej do osi $x$. Zatem pochodna funkcji w punkcie $x$ jest tangensem kąta nachylenia albo po prostu nachyleniem stycznej do wykresu funkcji w tym punkcie do osi $x$.

Pochodna funkcji nie zawsze istnieje. Jeżeli funkcja jest nieciągła, np. funkcja znaku, $[x]$, która przyjmuje wartość 1 dla $x>0$, 0 dla $x=0$ i -1 dla $x<0$, pochodna nie istnieje w punkcie $x=0$. Jeżeli funkcja ma ostrze” (np. $f(x)=\vert x\vert$) to pochodna w tym punkcie również nie istnieje (pochodna bezwzględnej wartości x wynosi $[x]$ z wyłączeniem $x=0$, w którym jest nieokreślona). Funkcję ciągłą, której pochodna jest zawsze określona i ciągła nazywamy funkcją klasy $C^1$.

Wybrane wzory na obliczanie pochodnych są podane w dodatku Użyteczne wzory”. Jest ważne, że reguły obliczania pochodnej zwane regułami różniczkowania umożliwiają obliczenia pochodnej każdej funkcji danej wzorem jawnym lub niejawnym (ten drugi przypadek zachodzi np. jeżeli obliczamy energię rozwiązując równania Hartree-Focka-Roothaana dla danej cząsteczki).

Różniczkując pochodną otrzymujemy drugą pochodną, kolejne różniczkowanie daje trzecią, czwartą i wyższe pochodne. Mając wszystkie pochodne funkcji w danym punkcie można obliczyć jej wartość w punktach sąsiednich, jest to rozwinięcie w szereg Taylora, które będzie często wykorzystywane w toku tego kursu:

$$f(x+h)=f(x)+f^\prime(x)h+\frac{1}{2}f^{\prime\prime}(x)h^2 +\frac{1}{3!}f^{\prime\prime\prime}(x)h^3+\ldots+\frac{1}{n!}f^{(n)}(x)h^n+\ldots$$

$$= \sum\limits_{n=0}^\infty\frac{1}{n!}f^{(n)}(x)h^n$$

Zauważmy, że zerowa” pochodna funkcji, $f^{(0)}(x)$ jest oryginalną funkcją.

Minima i maksima funkcji jednej zmiennej.

Minimum jest definiowane jako taki punkt, w którym wartość funkcji jest najmniejsza w pewnym otoczeniu a maksimum jako taki, dla którego wartość funkcji jest największa w pewnym otoczeniu. Są to minima i maksima lokalne. W zagadnieniach fizycznych minima i maksima są zwykle minimami i maksimami właściwymi co oznacza, że wykres funkcji nie ma w tym miejscu ostrza” (przykładem funkcji z ostrzem” jest wartość bezwzględna, $f(x)=\vert x\vert$). Poniższy rysunek przedstawia wykres funkcji z minimum i maksimum właściwym (po lewej stronie) oraz minimum i maksimum niewłaściwym (po prawej stronie).

W minimum lub maksimum właściwym styczna do wykresu funkcji jest równoległa do osi $x$ co oznacza, że $f^\prime(x)=0$ (korzystamy z interpretacji pochodnej jako tangensa kąta nachylenia stycznej do wykresu funkcji do osi $x$). Zatem poszukiwanie minimum lub maksimum właściwego sprowadza się do znalezienia miejsc zerowych pierwszej pochodnej, co jest zilustrowane na poniższym rysunku.

Z powyższego rysunku również widać, że pierwsza pochodna jest funkcją rosnącą (zmienia znak z -” na +”) w minimum a malejącą (zmienia znak z +” na -”) w maksimum. Stąd bezpośrednio wynika, że do rozstrzygnięcia czy dany punkt jest minimum czy maksimum wystarczy znak drugiej pochodnej: dodatni wskazuje ma minimum a ujemny na maksimum. Do tego wniosku można dojść analizując pierwsze człony rozwinięcia funkcji w otoczeniu punktu krytycznego ($x^*$). Ponieważ w punkcie krytycznym $f^\prime(x^*)=0$, w dostatecznie mały otoczeniu punktu krytycznego mamy:

$$f(x)\approx f(x^*)+\frac{1}{2}f^{\prime\prime}(x^*)(x-x^*)^2$$

Ponieważ w wyrażeniu występuje kwadrat odchylenia od wartości $x^*$, zatem jeżeli $f^{\prime\prime}(x^*)>0$ to wychylenie z punktu $x^*$ doprowadzi do zwiększenia wartości funkcji (minimum) a jeżeli $f^{\prime\prime}(x^*)<0$ do zmniejszenia jej wartości (maksiumum).

Warunki konieczne i wystarczające minimum i maksimum właściwego można zatem zapisać następująco:

$$f^\prime(x)=0,\ \ f^{\prime\prime}(x)>0:\ \ \ \mathrm{minimum}$$

$$f^\prime(x)=0,\ \ f^{\prime\prime}(x)<0:\ \ \ \mathrm{maksimum}$$

Do zastanowienia:

Jak można określić charakter punktu krytycznego jeżeli $f^{\prime\prime}(x^*)=0$? Czy zawsze miejsce zerowe pochodnej funkcji odpowiada jej minimum lub maksimum?

Pochodne funkcji wielu zmiennych.

Pojęcie pochodnych można łatwo uogólnić na funkcje wielu zmiennych. Rozważmy funkcję dwóch zmiennych $f(x,y)$ (np. energię układu w przemianie cyjanowodoru w izocyjanowodór w funkcji odległości C$\cdots$H i N$\cdots$H). Wykresem tej funkcji jest powierzchnia w kartezjańskim układzie współrzędnych $(x,y,z)$. Pochodną cząstkową względem zmiennej $x$ obliczamy traktując $y$ jako stałą a pochodną względem $y$ traktując $x$ jako stałą. Dla odróżnienia od pochodnych zwyczajnych nie stosujemy notacji $f^\prime$” (w takim przypadku nie byłoby wiadomo jakiej zmiennej dotyczy różniczkowanie) ani $\frac{\mathrm{d}f}{\mathrm{d}x}$ i $\frac{\mathrm{d}f}{\mathrm{d}y}$ ale odpowiednio $\frac{\partial f}{\partial x}$ i $\frac{\partial f}{\partial y}$. Przykładowo, dla funkcji $f(x,y)=x^2+y-y^2$ mamy:

$$\frac{\partial f}{\partial x}=2x$$

$$\frac{\partial f}{\partial y}=1-2y$$

Geometrycznie moża to sobie wyobrazić tak, że robimy przekrój trójwymiarowego wykresu funkcji odpowiednio płaszczyzną równoległą do płaszczyzny $yz$ i przechodzącą przez punkt o ustalonej wartości $y$ lub płaszczyzną równoległą do płaszczyzny $xz$ i przechodzącą przez punkt o ustalonej wartości $x$. Na tych płaszczyznach mamy już dwuwymiarowe wykresy funkcji i odpowiednie pochodne cząstkowe możemy interpretować jako nachylenia stycznych do wykresów.

Analogicznie definiujemy drugie pochodne funkcji. Dla funkcji dwóch zmiennych występują 4 drugie pochodne, ponieważ np. pierwszą pochodną względem $x$ można zróżniczkować albo względem $x$ albo względem $y$. Notacja podana jest poniżej.

$$\frac{\partial^2 f}{\partial x^2} = \frac{\partial}{\partial x}\left(\frac{\partial f}{\partial x}\right)$$

$$\frac{\partial^2 f}{\partial y\partial x} = \frac{\partial}{\partial y}\left(\frac{\partial f}{\partial x}\right)$$

$$\frac{\partial^2 f}{\partial x\partial y} = \frac{\partial}{\partial x}\left(\frac{\partial f}{\partial y}\right)$$

$$\frac{\partial^2 f}{\partial y^2} = \frac{\partial}{\partial y}\left(\frac{\partial f}{\partial y}\right)$$

Pochodne $\frac{\partial^2 f}{\partial y\partial x}$ i $\frac{\partial^2 f}{\partial y\partial x}$ są nazywane pochodnymi mieszanymi. Na mocy twierdzenia Schwartza, jeżeli pochodne mieszane są ciągłe to są sobie równe:

$$\frac{\partial^2 f}{\partial x\partial y}=\frac{\partial^2 f}{\partial y\partial x}$$

Uogólnienie pochodnych cząstkowych na pochodne wyższych rzędów, jak również na pochodne funkcji o liczbie zmiennych większych niż dwie jest oczywiste. Pierwsze pochodne względem kolejnych zmiennych tworzą wektor gradientu funkcji; wektor ten wskazuje kierunek największego wzrostu funkcji w danym punkcie. Wektor gradientu oznaczamy $\mathbf{\nabla}f$ (czytaj: nabla-ef) lub $\mathbf{grad}f$. Dla funkcji $n$ zmiennych, oznaczanych jako $x_1, x_2,\ldots x_n$ mamy:

$$\mathbf{\nabla}f=\mathbf{grad}f = \left(\begin{matrix}\frac{\part... ...{\partial x_2}\\ \vdots \\ \frac{\partial f}{\partial x_n}\end{matrix}\right)$$

Z kolei drugie pochodne tworzą macierz zwaną hesjanem funkcji:

$$\mathbf{H}=\left(\begin{matrix} \frac{\partial^2f}{\partial x_1^... ...x_n\partial x_2}&\ldots&\frac{\partial^2f}{\partial x_n^2}\end{matrix}\right)$$

Punkty krytyczne funkcji wielu zmiennych.

W porównaniu z funkcją jednej zmiennej, w przypadku funkcji wielu zmiennych dochodzi jeszcze jeden rodzaj punktów krytycznych: pukty siodłowe. W przypadku funkcji dwóch zmiennych punkt siodłowy jest takim punktem, że funkcja w jednym z kierunków posiada minimum a w kierunku do niego prostopadłym posiada maksimum. Inaczej mówiąc, w jednym z kierunków, w odpowiednio małym otoczeniu punktu siodłowego wartość funkcji jest większa a w kierunku do niego prostopadłym mniejsza niż w tym punkcie. Geometrycznie otoczenie punktu siodłowego rzeczywiście przypomina siodło. Jeżeli ktoś woli wycieczki w góry to może sobie punkt siodłowy wyobrazić jako przełęcz. W przypadku funkcji więcej niż dwóch zmiennych (np. $n$ zmiennych) mówimy ogólnie o punktach siodłowych rzędu $m$ ($1\leq m<n$); w takim przypadku funkcja posiada maksimum w $m$ prostopadłych do siebie kierunkach a w pozostałych $n-m$ kierunkach posiada minimum. W przypadku minimum wartość funkcji jest w odpowiednio małym otoczeniu zawsze większa niż wartość funkcji w minimum a w przypadku maksiumum jest znacznie mniejsza.

Poniższy rysunek ilustruje minima, maksima i punt siodłowy funkcji dwóch zmiennych. Można zauważyć, że przejście przez punkt siodłowy od jednego minimum do drugiego wymaga pokonania najmniejszej bariery. Dlatego stany przejściowe w reakcjach chemicznych czy też przemianach konformacyjnych odpowiadają punktom siodłowym pierwszego rzędu. Procesy te nie biegną przez punkty siodłowe rzędów wyższych niż 1, ponieważ idąc z takiego punktu siodłowego można zawsze znaleźć punkt siodłowy pierwszego rzędu, który łączy stan początkowy ze stanem końcowym a odpowiada mniejszej niż punkt siodłowy wyższego rzędu barierze energetycznej. Dlatego w analizie hiperpowierzchni energii potencjalnej układów molekularnych znajdujemy jedynie minima (które odpowiadają substratom, produktom i produktom przejściowym w badaniu reakcji chemicznych lub stabilnym konformacjom w analizie konformacyjnej) i punkty siodłowe pierwszego rzędu (które odpowiadają stanom przejściowym). Wyjątkiem są sytuacje, w których energia jest funkcją jednej zmiennej; wtedy stany przejściowe odpowiadają maksimom.

Z rysunku można zauważyć, że w dowolnym punkcie krytycznym płaszczyzna styczna do wykresu jest równoległa do płaszczyzny $xy$ a to oznacza, że obie pochodne cząstowe funkcji są równe zeru. Podobnie będzie w przypadku funkcji więcej niż dwóch zmiennych. Zatem warunek konieczny aby dany punkt był punktem krytycznym ma następującą postać:

$$\frac{\partial f(x_1^*,x_2^*,\ldots x_n^*)}{\partial x_1}=0,\ \... ...l x_2}=0,\ldots, \frac{\partial f(x_1^*,x_2^*,\ldots x_n^*)}{\partial x_n}=0$$

Ponieważ w punkcie krytycznym druga pochodna nie jest, jak w przypadku funkcji jednej zmiennej, liczbą lecz macierzą, badanie charakteru punktu krytycznego jest dla funkcji wielu zmiennych bardziej skomplikowane. Okazuje się, że w otoczeniu punktu krytycznego można zawsze wybrać taki układ współrzędnych, że drugie pochodne mieszane są w nim równe zeru i hesjan staje się macierzą diagonalną. Rozważmy przypadek funkcji dwóch zmiennych i napiszmy przybliżone wyrażenie na energię w otoczeniu punktu krytycznego $(x^*,y^*)$ biorąc, podobnie jak w przypadku funcji jednej zmiennej, rozwinięcie w szereg Taylora do drugiego rzędu włącznie. Pamiętamy, że pierwsze pochodne w punkcie krytycznym są równe zeru oraz, że drugie pochodne mieszane są sobie równe.

$$f(x,y)\approx f(x^*,y^*)+\frac{1}{2}\biggl[ \frac{\partial f(x^*,... ...x^2}(x-x^*)^2+ 2\frac{\partial f(x^*,y^*)}{\partial x\partial y}(x-x^*)(y-y^*)+$$

$$\frac{\partial f(x^*,y^*)}{\partial y^2}(y-y^*)^2\biggr]= f(x,^*,... ...}{2}(x-x^*,y-y^*)\mathbf{H}\left(\begin{matrix}x-x^*\\ y-y^*\end{matrix}\right)$$

Wyrażenie typu $\mathbf{z}^T\mathbf{Hz}$, gdzie $\mathbf{z}$ jest dowolnym wektorem niezerowym, nazywa się formą kwadratową macierzy $\mathbf{H}$. Można zatem powiedzieć, że energia w dostatecznie małym otoczeniu punktu krytycznego, liczona w odniesieniu do energii w punkcie krytycznym, jest w przybliżeniu formą kwadratową hesjanu.

Przesuwamy teraz początek układu współrzędnych do punktu $(x^*,y^*)$ oraz obracamy osie układu tak, aby w przekształconym układzie współrzędnych znikły mieszane pochodne cząstkowe (na wykresie warstwicowym funkcji są to kierunki prostopadłe do jej warstwic). Oznaczmy nowe współrzędne przez $\xi$ oraz $\eta$. Oczywiście, osie $\xi$ i $\eta$ są do siebie prostopadłe tak jak osie oryginalnego układu współrzędnych. Wtedy przybliżone wyrażenie na wartość funkcji w otoczeniu punktu krytycznego (który w nowym układzie współrzędnych ma współrzędne (0,0)) ma postać:

$$f(\xi,\eta) = f(0,0) + \frac{1}{2} \left[\frac{\partial^2 f(0,0)}{\partial\xi^2}\xi^2+\frac{\partial^2 f(0,0)}{\partial\eta^2}\eta^2\right]$$

Prawa strona zawiera wyrażenia zależące tylko od kwadratu $\xi$ albo $\eta$. Zatem jeżeli obie drugie pochodne w nowym układzie współrzędnych będą większe od zera to punkt będzie minimum (wartość funkcji po oddaleniu się od punktu krytycznego w dowolnym kierunku będzie się zwiększać), jeżeli obie będą mniejsze od zera to maksimum (wartość funkcji po oddaleniu się od punktu krytycznego w dowolnym kierunku będzie się zmniejszać) a jeżeli jedna będzie większa a druga mniejsza od zera to punkt będzie punktem siodłowym (wartość funkcji rośnie wzdłuż jednego kierunku a maleje wzdłuż kierunku doń prostopadłego). Przykład współrzędnych ($\xi,\eta$) dla jednego z minimów oraz punktu siodłowego na ścieżce reakcji izomeryzacji cynajowodoru do izocyjanowodoru jest pokazany na poniższym rysunku.

Łatwo zauważyć, że znajdowanie współrzędnych $\xi$ i $\eta$ jest równoważne diagonalizacji hesjanu. Diagonalizacja polega bowiem na przekształceniu macierzy, przy pomocy transformacji podobieństwa, w macierz diagonalną. W przypadku funkcji dwóch zmiennych mamy:

$$\mathbf{V}^T\mathbf{H}\mathbf{V}= \left(\begin{matrix}\lambda_1&0\\ 0&\lambda_2\end{matrix}\right)$$

przy czym kolumny macierzy $\mathbf{V}$ tworzą wektory ortonormalne:

$$\mathbf{v_1}\circ\mathbf{v_1}=v_{11}^2+v_{21}^2=1$$

$$\mathbf{v_1}\circ\mathbf{v_2}=v_{11}v_{12}+v_{21}v_{22}=0$$

$$\mathbf{v_2}\circ\mathbf{v_2}=v_{12}^2+v_{22}^2=1$$

Kolumny macierzy wektorów własnych $\mathbf{V}$ wyznaczają kierunki osi $\xi$ i $\eta$, natomiast wartości własne są drugimi pochodnymi funkcji w tych kierunkach w punkcie krytycznym:

$$\mathbf{V}=\left(\begin{matrix}\xi_x&\eta_x\\ \xi_y&\eta_y\end{matrix}\right)$$

$$\lambda_1=\frac{\partial^2 f}{\partial\xi^2},\ \lambda_2=\frac{\partial^2f}{\partial\eta^2}$$

Aby zbadać charakter punktu krytycznego potrzebne są tylko wartości własne hesjanu. Znajduje się je rozwiązując równanie wiekowe dla wyznacznika macierzy:

$$\det(\mathbf{H}-\lambda\mathbf{I}) = \left\vert\begin{matrix} \f... ...}&\ldots&\frac{\partial^2 f}{\partial x_n^2}-\lambda\end{matrix}\right\vert=0$$

gdzie $\mathbf{I}$ jest macierzą jednostkową.

$$\mathbf{I}=\left(\begin{matrix}1&0&\ldots&0\\ 0&1&\ldots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\ldots&1\end{matrix}\right)$$

Równanie to ma $n$ pierwiastków rzeczywistych, co wynika wprost z faktu, że hesjan jest macierzą hermitowską.

Podsumowując, punkty krytyczne na hiperpowierzchni energii potencjalnej można scharakteryzować jak w poniższej tabelce.

Rodzaj	Charakterystyka	Znaczenie
Miminum	$\lambda_1>0,\lambda_2>0,\ldots\lambda_n>0$	Substrat, produkt, produkt przejściowy, stabilna konformacja
Punkt siodłowy 1-go rzędu	$\lambda_1<0,\lambda_2>0,\ldots\lambda_n>0$	Stan przejściowy
Punkt siodłowy $m$-tego rzędu ($m>1$)	$\lambda_1<0,\lambda_2<0,\ldots\lambda_m<0$, $\lambda_{m+1}>0\ldots\lambda_n>0$	Bez znaczenia
Maksimum	$\lambda_1<0,\lambda_2<0,\ldots\lambda_n<0$	Bez znaczenia dla $n>1$, stan przejściowy dla $n=1$

Przykłady

Przykład 1

Znaleźć długość wiązania odpowiadającą minimum energii cząsteczli azotu, jeżeli energia potencjalna cząsteczki w zależności od długości wiązania ($d$) jest dana następującym przybliżonym wzorem:

$$E(d) = 815(d-1,09\ \mathrm{$\AA$})^2$$

Rozwiązanie:

W pierwszym kroku obliczamy pochodną energii względem $d$.

$$E^\prime(d) = (815d^2-2\times 815\times 1,09d-815\times 1,09^2)^\prime$$

$$=2\times815d-2\times 815\times 1,09=2\times 815(d-1,09)$$

Jeżeli zastosować wzór na pochodną funkcji złożonej to można uniknąć podnoszenia wyrażenia w nawiasie do kwadratu i późniejszego zwijania wyrażenia na pochodną:

$$E(g(d))=815g(d)^2, g(d)=d-1,09$$

$$E^\prime(g(d)) = E^\prime(g)g^\prime(d) = 2\times 815(g(d))\times 1 = 2\times 815(d-1,09)$$

W następnym kroku szukamy miejsca zerowego pochodnej. W takim punkcie na układ nie działają żadne siły czyli energia może mieć minimum, maksimum lub punkt przegięcia:

$$E^{\prime}(d^*)=2\times 815(d^*-1,09)=0 \Rightarrow d^*=1,09$$

Pozostaje ustalić czy punkt ten stanowi minimum. W tym celu liczymy drugą pochodną energii i sprawdzamy jej znak:

$$E^{\prime\prime}(d^*)=(2\times 815(d^*-1,09))^\prime = 2\times 815 > 0$$

Zatem znaleziony punkt stanowi minimum energii.

Wykresy energii i jej pochodnej w funkcji długości wiązania są podane na poniższym rysunku.

Przykład 2

Dany jest układ składający się z protonu położonego pomiędzy dwoma takimi samymi atomami, będącymi akceptorem protonu, gdzie istnieje możliwość przeniesienia protonu od jednego akceptora do drugiego:

$$\large\mathrm{A-H\cdots A \rightleftharpoons A\cdots H-A}$$

Zakładając, że początek układu współrzędnych znajduje się w połowie odległości między atomami akceptora X i że proton porusza się tylko na osi X$\cdots$X, krzywa energii potencjalnej jest dana wzorem:

$$E(x) = a(x^2 - b^2)^2 - c$$

gdzie $a$, $b$ i $c$ są parametrami (przy czym $a>0$ i $b\geq 0$ a $x$ jest współrzędną określającą położenie protonu między atomami X, znaleźć położenia wszystkich minimów i maksimów oraz odpowiadające im wartości energii potencjalnej, jak również barierę przejścia między dwoma stanami w procesie przeniesienia protonu, jako funkcje parametrów $a$, $b$ i $c$.

Rozwiązanie:

Najpierw szkicujemy badany układ z zaznaczeniem zmiennej funkcji energii (współrzędnej reakcji):

Korzystając ze wzoru na pochodną funkcji złożonej, obliczamy pochodną i znajdujemy jej miejsca zerowe:

$$E^\prime(x)=2a(x^2-b^2)\times 2x=0$$

$$x^2-b^2=0\ \mathrm{lub}\ x=0$$

$$x_1 = -b,\ x_2=0,\ x_3=b$$

Następnie, korzystając z wzoru na pochodną iloczynu, obliczamy drugą pochodną i wyliczamy jej wartości w punktach $x_1$, $x_2$ i $x_3$.

$$E^{\prime\prime}(x)=2a\times 2x\times 2x+2a(x^2-b^2)\times 2=4a(3x^2-b^2)$$

$$E^{\prime\prime}(-b) = 4a(3(-b)^2-b^2) = 8ab^2 > 0 \ (a>0)$$

$$E^{\prime\prime}(0) = 4a(0-b^2) = -4ab^2 < 0$$

$$E^{\prime\prime}(b) = 4a(3b^2-b^2) = 8ab^2 > 0$$

Zatem punkty $x_1$ i $x_3$ są minimami energii a punkt $x_2$ stanowi maksimum.

Obliczamy wartości energii w tych trzech punkach:

$$E(x_1)=E(x_3)=-c$$

$$E(x_2)=ab^4-c$$

Bariera energetyczna jest różnicą energii w maksimum i w minimum:

$$\Delta E^* = E(x_2)-E(x_1)=E(x_3)-E(x_1)=ab^4$$

Bariera bardzo szybko maleje z parametrem określającym położenia minimów ($b$). Gdy $b=0$, bariera znika i na krzywej energii pozostaje tylko jedno minimum, odpowiadające położeniu protonu dokładnie pomiędzy atomami akceptora. Taka sytuacja występuje w układach z bardzo silnym wiązaniem wodorowym, np. kwaśnych fluorkach albo kwaśnych maleinianach.

Krzywe energii dla $a=1$, $c=0$ oraz trzech różnych wartości $b$ ilustruje rysunek.

Przykład 3

Znaleźć położenia oraz wartości kąta $\omega$ (zdefiniowanego na poniższym rysunku) wszystkich minimów i maksimów energii torsyjnej w przedziale $[-\pi,\pi]$ dla obrotu wokół wiązania C-N grupy peptydowej w cząsteczce N-metyloacetamidu (CH$_3$-CONH-NHCH$_3$) oraz obliczyć barierę przejścia między strukturami odpowiadającymi minimom energii. Jakie konfiguracje cząsteczki odpowiadają tym minimom?

Energia potencjalna w funkcji kąta $\omega$ jest dana następującym równaniem:

$$E(\omega) = 1,2[1+\cos\omega] + 10[1-\cos(2\omega)]$$

Rozwiązanie:

Korzystając z wzorów na pochodną funkcji cosinus oraz na pochodną funkcji złożonej, obliczamy pochodną energii:

$$E^\prime(\omega) = -1,2\sin\omega + 2\times 10\sin(2\omega)$$

$$=-1,2\sin\omega+20\sin(2\omega) = \sin\omega(40\cos\omega-1,2)$$

W ostatnim przekształceniu skorzystaliśmy ze wzoru $\sin(2\omega)=2\sin\omega\cos\omega$.

Pochodna ma miejsca zerowe gdy $\sin\omega=0$ lub $\cos\omega=1,2/40=0,03$. Zatem na krzywej energii mamy w przedziale $[-\pi,\pi]$ (lub $[-180^\circ,180^\circ]$) następujące punkty krytyczne:

$$\omega_1=-\pi\ (-180^\circ),\ \omega_2=-1,5408\ (-88,28^\circ), \omega_3=0,\ \omega_4=1,5408\ (88,28^\circ),\ \omega_5=\pi\ (180^\circ)$$

Aby określić charaktery tych punktów, obliczamy drugą pochodną:

$$E^{\prime\prime}(\omega) = -1,2\cos\omega+2\times 20\cos(2\omega)$$

$$=-1,2\cos\omega + 40(2\cos^2\omega - 1)$$

Skorzystaliśmy z wzoru redukcyjnego $\cos(2\omega)=2\cos^2\omega-1$.

$$E^{\prime\prime}(-\pi)=E^{\prime\prime}(\pi)=-1,2(-1)+40(2\times(-1)^2-1)=41,2\ \Rightarrow\ \mathrm{minimum}$$

$$E^{\prime\prime}(-1,5408)=E^{\prime\prime}(1,5408)=-1,2\times 0,03+40(2\times 0,03^2-1) =$$

$$-39,964 < 0\ \Rightarrow\ \mathrm{maksimum}$$

$$E^{\prime\prime}(0)=-1,2+40(2-1)=38,8>0\ \Rightarrow\ \mathrm{minimum}$$

W drugiej linii podstawiliśmy $\cos\omega_2=\cos\omega_4=0,03$.

W tabelce są zestawione poszczególne punkty krytyczne, ich charakter oraz wartości energii a wykres energii jest podany poniżej. Minimum, dla którego $\omega=0$ jest konformacją ^2.1 cis a to z $\omega=\pi\ (180^\circ)$ konformacją trans. Różnica energii między konformacją trans i cis wynosi $\Delta E = E(0) - E(\pi)=2,4$ kcal/mol a bariera przejścia od konformacji cis do konformacji trans wynosi $\Delta E^*=E(1,5408)-E(0)=20,018$ kcal/mol (aby obliczyć barierę trzeba skorzystać z wzoru redukcyjnego na $\cos(2\omega)$ i wstawić do wzoru wartość cosinusa w maksimum).

$\omega [^\circ]$	$-180$	$-88,28$	0	$88,28$	$180$
Charakter	min.	maks.	min.	maks.	min.
$E(\omega)$[kcal/mol]	0	$20,018$	$2,4$	$20,018$	0

Przykład 4

Fragment powierzchni energii przeniesienia protonu od donora $A$ do akceptora $B$ w pewnym układzie można opisać następującym równaniem, w którym $x=d_{\mathrm{A\cdots H}}$ a $y=d_{\mathrm{B\cdots H}}$:

$$E(x,y) = 2(x-y)^2-(x+y-3)^2$$

Znaleźć punkt lub punkty krytyczne tej funkcji i określić ich charakter (minium, maksimum lub punkt siodłowy).

Rozwiązanie:

W pierwszym kroku wyliczamy pochodne cząstkowe funkcji względem poszczególnych zmiennych:

$$\frac{\partial E}{\partial x} = 2\times 2(x-y)-2(x+y-3) = 2x-6y+6$$

$$\frac{\partial E}{\partial y} =-2\times 2(x-y)-2(x+y-3) =-6x+2y+6$$

Punktowi krytycznemu odpowiada taka para $(x^*,y^*)$, że obie pochodne cząstkowe przyjmują wartość zero. Znalezienie tego punktu wymaga roziązania układu dwóch równań liniowych z dwiema niewiadomymi:

$$2x-6y=-6$$

$$-6x+2y=-6$$

Rozwiązaniem jest $x^*=1,5$ $\AA$, $y^*=1,5$ $\AA$.

Charakter punktu krytycznego określamy poprzez obliczenie wartości własnych macierzy drugich pochodnych (hesjanu) energii w punkcie krytycznym:

$$\frac{\partial^2 E}{\partial x^2} = \frac{\partial{(2x-6y+6)}}{\partial x}=2$$

$$\frac{\partial^2 E}{\partial y\partial x} = \frac{\partial{(2x-6y+6)}}{\partial y}=-6$$

$$\frac{\partial^2 E}{\partial x\partial y} = \frac{\partial{(-6x+2y+6)}}{\partial x}=-6$$

$$\frac{\partial^2 E}{\partial y^2} = \frac{\partial{(-6x+2y+6)}}{\partial y}=2$$

Zwróćmy uwagę, że wszystkie drugie pochodne są liczbami. Aby wyliczyć wartości własne, układamy wyznacznik wiekowy hesjanu:

$$\det(\mathbf{H}-\lambda\mathbf{I}) = \left\vert\begin{matrix}2-\lambda&6\\ 6&2-\lambda\end{matrix}\right\vert=0$$

Po rozwinięciu wyznacznika dostajemy równanie kwadrartowe na $\lambda$

$$(2-\lambda)^2-36=0\ \Rightarrow\ 2-\lambda=\pm 6$$

Stąd $\lambda_1=-4$, $\lambda_2=8$. Punkt krytyczny jest punktem siodłowym (jedna wartość własna hesjanu jest ujemna).

Powierzchnia energii potencjalnej jest przedstawiona na poniższym rysunku.

Przykład 5$^*$

Powierzchnia energii konformacyjnej pewnej cząsteczki jest opisywana poprzez dwa kąty torsyjne, oznaczone jako $\tau_1$ i $\tau_2$. Znaleźć minima, punkty siodłowe i maksima na tej powierzchni zakładając, że oba kąty torsyjne leżą w przedziale $[-\pi,\pi]$ ( $[-180^\circ,180^\circ]$).

$$E(\tau_1,\tau_2)=\cos(\tau_1+\tau_2)-\cos(\tau_1-\tau_2)$$

Tak, jak w przykładzie 4, obliczamy pochodne cząstkowe względem każdej ze zmiennych i znajdujemy ich miejsca zerowe:

$$\frac{\partial E}{\partial\tau_1} = -\sin(\tau_1+\tau_2)+\sin(\tau_1-\tau_2)=0$$

$$\frac{\partial E}{\partial\tau_2} = -\sin(\tau_1+\tau_2)-\sin(\tau_1-\tau_2)=0$$

Dodając a następnie odejmując powyższe równania stronami dostajemy:

$$-2\sin(\tau_1+\tau_2)=0\ \Rightarrow\ \tau_1+\tau_2=k\pi, k\in C$$

$$2\sin(\tau_1-\tau_2)=0\ \Rightarrow\ \tau_1-\tau_2=m\pi, m\in C$$

Dodając i odejmując stronami prawe strony powyższych równań i dzieląc wyniki obustronnie przez 2 dostajemy:

$$\tau_1=\frac{k+m}{2}$$

$$\tau_2=\frac{k-m}{2}$$

Liczby całkowite $k$ i $m$ należy dobrać tak, aby zarówno $\tau_1$, jak i $\tau_2$ mieściły się w przedziale $[-\pi,\pi]$. Z powyższego wzoru wynika wprost, że $\tau_1$ i $\tau_2$ mogą być albo jednocześnie wielokrotnościami $\pi\ (180^\circ)$ albo $\pi/2\ (90^\circ)$, natomiast kombinacje takie jak np. $\tau_1=\pi$, $\tau_2=\pi/2$ nie są punktami krytycznymi na powierzchni energii.

Charaktery poszczególnych punktów krytycznych są zestawione w poniższej tabelce. Dla przykładu poniżej dalej podane określenie charakteru punktów o współrzędnych odpowiednio $(\pi/2,\pi/2)$, $(\pi/2,-\pi/2)$ oraz $(0,0)$; obliczenia w pozostałych przypadkach są podobne.

$\tau_1$$\tau_2$	$-\pi$	$-\pi/2$	0	$\pi/2$	$\pi$
$-\pi$	siodło		siodło		siodło
$-\pi/2$		min		Maks
0	siodło		siodło		siodło
$\pi/2$		Maks		min
$\pi$	siodło		siodło		siodło

Obliczamy drugie pochodne:

$$\frac{\partial^2 E}{\partial\tau_1^2}=-\cos(\tau_1+\tau_2)+\cos(\... ...ial^2 E}{\partial\tau_2\partial\tau_1}=-\cos(\tau_1+\tau_2)-\cos(\tau_1-\tau_2)$$

$$\frac{\partial^2 E}{\partial\tau_1\partial\tau_2}=-\cos(\tau_1+\t... ... \frac{\partial^2 E}{\partial\tau_2^2}=-\cos(\tau_1+\tau_2)+\cos(\tau_1-\tau_2)$$

Dla $\tau_1=\pi/2$, $\tau_2=\pi/2$:

$$\mathbf{H}=\left(\begin{matrix}-\cos\pi+\cos 0&-\cos\pi-\cos 0\\ ... ...\cos 0\end{matrix}\right)= \left(\begin{matrix}2&0\\ 0&2\end{matrix}\right)$$

Ponieważ hesjan ma tylko niezerowe elementy diagonalne, każdy z nich jest jego wartością własną. Obie wartości własne są dodatnie a zatem punkt $(\pi/2,\pi/2)$ stanowi minimum energii. Ponieważ cosinus jest funkcją parzystą, jego siostrzanym punktem jest $(-\pi/2,\pi/2)$.

Dla $\tau_1=\pi/2$, $\tau_2=-\pi/2$:

$$\mathbf{H}=\left(\begin{matrix}-\cos 0+\cos\pi&-\cos 0-\cos\pi\\ ... ...s\pi\end{matrix}\right)= \left(\begin{matrix}-2&0\\ 0&-2\end{matrix}\right)$$

Ten punkt stanowi maksimum (oba elementy diagonalne, które są jednocześnie wartościami własnymi, są ujemne). Jego siostrzanym punktem jest $(-\pi,\pi)$.

Dla $\tau_1=0$, $\tau_2=0$:

$$\mathbf{H}=\left(\begin{matrix}-\cos 0+\cos 0&-\cos 0-\cos 0\\ ... ...os 0\end{matrix}\right)= \left(\begin{matrix}0&-2\\ -2&0\end{matrix}\right)$$

Równanie wiekowe ma postać:

$$\det(\mathbf{H}-\lambda\mathbf{I})=\left\vert\begin{matrix}-\lamb... ...end{matrix}\right\vert\ \Rightarrow\ \lambda^2=4,\ \lambda_1=-2,\ \lambda_2=2$$

Ten punkt jest punktem siodłowym. Jego siostrzanymi punktami są wszystkie punkty takie, że każdy z kątów wynosi $-\pi$, 0 lub $\pi$.

Poniższy rysunek ilustruje badaną powierzchnię energii oraz znalezione na niej punkty krytyczne.

Przykład 6

W poniższej tabeli są zestawione energie oraz wartości własne czterech punktów krytycznych pewnej reakcji chemicznej, którą można opisać przy pomocy trzech współrzędnych. Jeden z tych punktów stanowią substraty a jeden produkty.

Punkt	1	2	3	4
Energia [kcal/mol]	-210	-260	-200	-185
$\lambda_1$ [kcal/mol/$\AA$$^2$]	100	155	-150	-100
$\lambda_2$ [kcal/mol/$\AA$$^2$]	125	295	55	-31
$\lambda_3$ [kcal/mol/$\AA$$^2$]	300	489	176	52

Na podstawie tych danych:

1. Określić charakter punktów krytycznych i ocenić, które z nich są punktami na ścieżce reakcji.

2. Wiedząc, że reakcja jest egzotermiczna, określić, który punkt odpowiada substratom, który produktom, a który stanowi przejściowemu.

3. Obliczyć energię reakcji oraz barierę energetyczną reakcji.

Rozwiązanie:

1. Punkty 1 i 2 są minimami energii (wszystkie wartości własne hesjanu są dodatnie), punkt 3 jest punktem siodłowym pierwszego rzędu (jedna wartość własna ujemna) a punkt 4 jest punktem siodłowym drugiego rzędu (dwie wartości własne ujemne). Punkt 4 nie jest punktem na ścieżce reakcji ponieważ reakcje biegną tylko przez punkty siodłowe pierwszego rzędu.

2. Spośród minimów niższa energia odpowiada punktowi 2, zatem punkt 1 odpowiada substratom a punkt 2 produktom reakcji. Punkt 3 odpowiada stanowi przejściowemu reakcji.

3. Efekt energetyczny reakcji wynosi $\Delta E=E_2-E_1=-260-(-210)=-50$ kcal/mol (odejmujemy energię substratów od energii produktów) a bariera energetyczna reakcji wynosi $\Delta E^*=-200-(-210) = 10$ kcal/mol (odejmujemy energię substratów od energii stanu przejściowego).

Zadania

Uwaga!

W poniższych zadaniach odległości i długości wiązań są dane w angströmach, kąty płaskie w stopniach a energia jest wyrażona w kcal/mol.

Zadanie 1

Znaleźć wartość kąta H-O-H ($\alpha$) w cząsteczce wody, która odpowiada minimum energii cząsteczki jeżeli energia jest w przybliżeniu dana poniższym równaniem:

$$E(\alpha)=47(\alpha-104,5)^2$$

Odpowiedź:

$\alpha^*=104,5^\circ$

Zadanie 2

Znaleźć odległość $r^*$ odpowiadającą minimum energii układu dwóch atomów helu, których energia jest dana poniższym równaniem. Sprawdzić, czy znaleziony punkt odpowiada minimum energii. Obliczyć energię w minimum.

$$E(r)=0,02\left[\left(\frac{2,3}{r}\right)^{12}-2\left(\frac{2,3}{r}\right)^6\right]$$

Odpowiedź:

$r^*=2,3$ $\AA$, $E^*=-0.02$ kcal/mol.

Zadanie 3

Znaleźć położenie minimum energii w funkcji stałych $C$ i $D$ potencjału 12-10, używanego do obliczania energii oddziaływania proton-akceptor protonu w wiązaniach wodorowych. Energia jest dana poniższym równaniem. Sprawdzić, czy znaleziony punkt jest minimum energii.

$$E(r)=\frac{C}{r^{12}}-\frac{D}{r^{10}}$$

Odpowiedź:

$r^*=\sqrt{\frac{12C}{10D}}$.

Zadanie 4$*$

Dobrać tak parametr $k$ w funkcji wykładników $n$ i $m$ ($n>m$) w potencjale $n-m$, aby minimum energii wypadało dokładnie przy wartości $r^\circ$. Jaka będzie wtedy wartość energii w minimum?

$$E(r)=\epsilon\left[\left(\frac{r^\circ}{r}\right)^n-k\left(\frac{r^\circ}{r}\right)^m\right]$$

Odpowiedź:

$k=n/m$, $E^*=-\epsilon$.

Zadanie 5

Znaleźć wszystkie minima i maksima w przedziale $[-\pi,\pi]$ energii w funkcji kąta obrotu $\tau$ wokół wiązania S-S w disiarczku dimetylu CH$_3$-S-S-CH$_3$, jeżeli energia jest dana poniższym równaniem:

$$E(\tau)=3,5[1+\cos(2\tau)]$$

Określić barierę energetyczną obrotu wokół tego wiązania oraz naszkicować przybliżony wykres energii w funkcji kąta $\tau$.

Odpowiedź:

Minima: $\tau=\pm\pi/2$, maksima: $\tau=-\pi, 0, \pi$, $\Delta E^*=7$ kcal/mol.

Zadanie 6

Znaleźć wszystkie minima i maksima w przedziale $[-\pi,\pi]$ energii jako funkcji kąta obrotu wokół wiązania C-C w cząsteczce etanu.

$$E(\tau)=1,3[1+\cos(3\tau)]$$

Odpowiedź: Minima: $\tau=-\pi,-\pi/3,\pi/3,\pi$, maksima: $\tau=-2\pi/3, 0, 2\pi/3$.

Zadanie 7

Znaleźć punkty krytyczne na wycinku hiperpowierzchni energii potencjalnej przeniesienia protonu od donora $A$ do akceptora $B$, danej jako funkcja odległości $A\cdots H$ ($x$) oraz $B\cdots H$ ($y$), danej następującym wzorem:

$$E(x,y) = 2(x-y)^2+(x+y-3)^2$$

Odpowiedź: $(x^*=1,5; y^*=1,5)$, minimum.

Zadanie 8

Energię potencjalną cząsteczki nadtlenku benzoilu Ph-C(=O)-O-O-(O=C)-Ph w zależności od kąta walencyjnego O-O-C $(\theta)$ oraz kąta torsyjnego obrotu wokół wiązania O-O $(\tau)$ można opisać równaniem:

$$E(\theta,\tau)=50 (\theta - 120^\circ)^2 + 2,5[1 + \cos(2\tau)]$$

Znaleźć punkty krytyczne oraz określić ich typ (minimum, maksimum, punkt siodłowy) na tej hiperpowierzchni energii potencjalnej, traktując $\theta$ i $\tau$ jako zmienne. Kąt $\theta$ zmienia się w przedziale $[0^\circ,180^\circ]$ a kąt $\tau$ w przedziale $[-180^\circ,180^\circ]$.

Odpowiedź:

Minima występują dla $(\theta=120^\circ,\tau=\pm 90^\circ)$ a punkty siodłowe dla $(\theta=120^\circ,\tau=-180^\circ,0^\circ, 180^\circ)$. Powierzchnia energii nie posiada maksimów.

Zadanie 9

Które zaznaczone (symbolami lub liczbami) punkty na załączonej powierzchni energii potencjalnej blokowanej alaniny są minimami a które punktami siodłowymi?

Zadanie 10

Na powierzchni energii potencjalnej odpowiadającej pewnej przemianie konformacyjnej, która może być opisana trzema parametrami geometrycznymi, mającymi sens kątów dwuściennych, znaleziono 4 punkty krytyczne, dla których wartości energii ($E$, kcal/mol) oraz wartości własne hesjanu energii [kcal/(mol$\times$$\AA$$^2$)] są następujące:

Punkt 1: $E=0, \lambda_1=10, \lambda_2=20, \lambda_3=30$

Punkt 2: $E=15, \lambda_1=-5, \lambda_2=-4, \lambda_3=-4$

Punkt 3: $E=5, \lambda_1=-5, \lambda_2=5, \lambda_3=10$

Punkt 4: $E=-8, \lambda_1=10, \lambda_2=10, \lambda_3=20.$

Określić charakter (minimum, maksimum, punkt siodłowy oraz rząd punktu siodłowego) każdego z tych punktów krytycznych. Który czy które z nich odpowiadają punktom na ścieżce przemiany konformacyjnej? Odpowiedzi uzasadnić.