Prawo rozkładu Boltzmanna

Wstęp teoretyczny

Układy zawierające miriady komponentów, a takimi są zbiorowiska atomów lub cząsteczek związków chemicznych, nie poddają się opisowi szczegółowemu. Dlatego narzędziem odpowiednim do ich opisu jest mechanika statystyczna, która umożliwia badania zachowania całego układu na podstawie cech poszczególnych rodzajów obiektów, z których on się składa. W naszym przypadku umożliwia powiązanie właściwości pojedynczych molekuł (w szczególności hiperpowierzchni energii potencjalnej) z tym co się dzieje w laboratorium mokrym. Kluczowe jest tutaj określenie prawdopodobieństwa, że atom, cząsteczka lub cały układ znajdzie się w danym stanie (np. że cząsteczka przyjmie określoną konformację). Mówi o tym prawo rozkładu Boltzmanna.

Prawo rozkładu Boltzmanna określa prawdopodobieństwo znalezienia się układu w określonym stanie w zależności od energii tego stanu oraz temperatury otoczenia. Jeżeli układ występuje w stanach, którym nadamy numery 1,2,$\ldots$, $Z$ ($Z$ może być nieskończonością), to prawdopodobieństwo stanu $i$ o całkowitej energii $E_i$ jest proporcjonalne do

$$P_i \sim \exp\left(-\frac{E_i}{k_BT}\right)$$

gdzie $T$ jest temperaturą bezwzględną a $k_B=1,3806\times 10^{-23}$ jest stałą Boltzmanna, której jednostką jest J/K. W powyższym równaniu energia jest wyrażona w joulach. Jeżeli chcemy energię wyrazić w joulach/mol, zamiast stałej Boltzmanna do wzoru wstawiamy uniwersalną stałą gazową $R=8,3145$ J/(mol$\times$K).

$$P_i \sim \exp\left(-\frac{E_i}{RT}\right)$$

Stała gazowa jest iloczynem stałej Boltzmanna i liczby Avogadro ( $N_A=6,02214\times 10^{23}$).

Aby obliczyć bezwzględne prawdopodobieństwo, trzeba powyższe wyrażenie znormalizować tak, aby prawdopodobieństwo wystąpienia któregokolwiek ze stanów układu było równe jedności. Będziemy zawsze zakładać, że wystąpienie danego stanu jest zdarzeniem niezależnym od wystąpienia innych stanów. Jak na razie nie jest znane doświadczenie które przeczyłoby temu założeniu. W takim razie prawdopodobieństwo wystąpienia jednego z kilku stanów jest sumą prawdopodobieństw wystąpienia poszczególnych stanów.

Znormalizowane (bezwzględne) prawdopodobieństwo wystąpienia stanu o numerze $i$ wynosi

$$P_i = \frac{\exp\left(-\frac{E_i}{RT}\right)}{Q}$$

gdzie wielkość

$$Q=\sum\limits_{j=1}^Z\exp\left(-\frac{E_j}{k_BT}\right)$$

jest nazywana sumą statystyczną.

Bardzo często wiele stanów ma taką samą energię. Przykładem w mechanice klasycznej jest punkt materialny, na który nie działają żadne siły. Jego energia jest równa energii kinetycznej. Stan tego punktu jest opisany poprzez podanie położenia i prędkości. Energia kinetyczna zależy jedynie od kwadratu długości wektora prędkości i jest zatem niezależna od miejsca, w którym znajduje się punkt i od kierunku, w którym się porusza. Istnieje zatem nieskończenie wiele stanów tego punktu, którym odpowiada ta sama energia. Przykładem z chemii jest cząsteczka tlenu w podstawowym stanie trypletowym. Istnieją dokładnie 3 realizacje tego stanu i w nieobecności zewnętrznych pól każdemu odpowiada identyczna energia. Jeżeli przez $g_i$ oznaczymy liczbę stanów o tej samej energii (inaczej: degenerację stanu $i$) a zamiast indywidualnych stanów będziemy rozpatrywać grupy stanów o danej energii (każdą taką grupę będziemy dla uproszczenia nazywać stanem) to wyrażenia na prawdopodobieństwo i sumę statystyczną przyjmą postać

$$P_i = \frac{g_i\exp\left(-\frac{E_i}{k_BT}\right)}{Q}$$

$$Q=\sum\limits_{j=1}^Zg_j\exp\left(-\frac{E_j}{k_BT}\right)$$

Jeżeli stanowi układu nie może być przyporządkowana liczba ale jest on opisywany zmienną lub wektorem zmiennych ciągłych to prawdopodobieństwo przechodzi w gęstość prawdopodobieństwa a suma statystyczna przechodzi w całkę. Najczęściej stan układu jest opisywany przez położenia i pędy. Wtedy wzór na gęstość prawdopodobieństwa przyjmuje postać

$$P(\mathbf{x},\mathbf{p})d^n\mathbf{x}d^n\mathbf{p} = \frac{\exp\left(-\frac{E(\mathbf{x},\mathbf{p})}{k_BT}\right)d^n\mathbf{x}d^n\mathbf{p}}{Q}$$

$$Q=\int\limits_{-\infty}^\infty\cdots\int\limits_{-\infty}^\infty \exp\left(-\frac{E(\mathbf{x},\mathbf{p})}{k_BT}\right)d^n\mathbf{x}d^n\mathbf{p}$$

Jeżeli pogrupujemy stany o tej samej energii, możemy zdefiniować prawdopodobieństwo jako funkcję jedynie energii

$$P(E)dE = \frac{\Omega(E)\exp\left(-\frac{E}{k_BT}\right)dE}{Q}$$

$$Q=\int\limits_{-\infty}^\infty \Omega(E)\exp\left(-\frac{E}{k_BT}\right)dE$$

gdzie $\Omega(E)$ jest gęstością stanów.

Często interesuje nas nie tyle bezwzględne prawdopodobieństwo danego stanu ale jego prawdopodobieństwo względne w stosunku do innego stanu. Jeżeli te stany oznaczymy przez I i II, to wzór przyjmuje postać

$$\frac{P_{II}}{P_I}=\frac{g_{II}\exp\left(-\frac{E_{II}}{k_BT}\rig... ...I}-E_I}{k_BT}\right)=\frac{g_{II}}{g_I}\exp\left(-\frac{\Delta E}{k_BT}\right)$$

Zatem żeby obliczyć stosunek prawdopododobieństw wystarczy różnica energii między stanami a nie trzeba znać bezwzględnych wartości energii. Należy zwrócić uwagę, że jeżeli zdefiniujemy

$$\Delta S = k_B\ln \frac{g_{II}}{g_I}$$

$$\Delta F = \Delta E - T\Delta S$$

jako odpowiednio różnicę entropii i energii swobodnej między stanami I i II to wyrażenie na stosunek prawdopodobieństw można zapisać następującym równaniem:

$$\frac{P_{II}}{P_I}=\exp\left(-\frac{\Delta F}{k_BT}\right)=K_{I\rightarrow II}$$

które jest wyrażeniem na zależność stałej równowagi reakcji przejścia formy I w formę II od temperatury jeżeli stan I i stan II odpowiadają formie I i II jakiegoś układu (np. dwóm formom tautomerycznym) a stosunek prawdopodobieństw wystąpienia każdej z form zidentyfikujemy ze stosunkiem ich ułamków molowych. Z tych rozważań również wynika, że entropia stanu jest proporcjonalna do logarytmu z jego degeneracji czyli liczby możliwości, na które stan może być zrealizowany.

Przykłady

W analizowanych przykładach oraz w zadaniach, z wyjątkiem przykładu 4 oraz zadań 5 i 6, będziemy operować energią przeliczoną na 1 mol. Zatem w równaniach będzie, poza wymienionymi wyjątkami, występować stała gazowa $R$ a nie stała Boltzmanna $k_B$.

Przykład 1

Oszacować ułamek molowy enolu acetyloacetonu w temperaturze pokojowej ($T=298$ K) jeżeli różnica energii formy enolowej i ketonowej wynosi 1,2 kcal/mol.

Rozwiązanie.

Przede wszystkim należy też zwrócić uwagę na fakt, iż ułamek molowy oznaczany literą $x$ jest liczbowo równy prawdopodobieństwu wystąpienia danej formy. Każda z form występuje tylko raz, więc degeneracja obu stanów wynosi 1. Zatem, zgodnie z rozważaniami podanymi w części teoretycznej, stosunek ułamka molowego formy enolowej do ułamka molowego formy ketonowej będzie równy

$$\frac{x_{enol}}{x_{keton}} = \exp\left(-\frac{E_{enol}-E_{keton}}{RT}\right) = \exp\left(-\frac{\Delta E}{RT}\right)$$

Po wyrażeniu $\Delta E$ w J/mol, $\Delta E=1,2\times 4184=5023,2$ J/mol oraz wstawieniu $R=8,3145$ J/(mol$\times$K) i $T=298$ K otrzymujmey $x_{enol}/x_{keton}=0,1317$. Ponieważ suma ułamków molowych formy enolowej i formy ketonowej jest równa 1

$$x_{keton} = 1 - x_{enol}$$

można łatwo wyliczyć zawartość formy enolowej.

$$\frac{x_{enol}}{1 - x_{enol}} = 0,1317$$

$$x_{enol}=\frac{0,1317}{1+0,1317} = 0,1163 \approx 11,6~\%$$

Przykład 2

W czasie 1 ns symulacji dynamiki molekularnej w temperaturze T=298 K wiązania kationu wapniowego przez kalbindynę (białko wiążące wapń) stwierdzono, że kation wapniowy przebywał w miejscu wiążącym przez 368 ps. Oszacować zmianę energii swobodnej po związaniu kationu wapniowego przez kalbindynę.

Rozwiązanie.

Wiązanie kationów wapniowych (zielone kulki) przez kalbindynę jest zilustrowane na poniższym rysunku. Kationy wapnia są wiązane przez grupy karboksylowe łańcuchów bocznych kwasu glutaminowego i asparaginowego znajdujące się w pętlach.

Powyższe zadanie można rozwiązać w analogiczny sposób jak poprzedni przykład. Definiujemy różnicę energii między stanem niezwiązanym ($nzw$) i związanym ($zw$).

$$\Delta E = E_{zw} - E_{nzw}$$

Przede wszystkim zauważmy, że zamiast energii” występuje tutaj energia swobodna”. To oznacza, że wielkość ta zawiera w sobie możliwą degenerację stanów oraz inne efekty związane z entropią takie jak przeorganizowanie struktury rozpuszczalnika i białka po związaniu. Tym niemniej w prawie rozkładu Boltzmanna energia swobodna pełni taką rolę jak energia.

Zakładamy ergodyczność symulacji, czyli możliwość odwiedzania obu stanów w toku całej symulacji proporcjonalnie do ich prawdopodobieństw. Wtedy zaobserwowanie w czasie 1 ns (1000 ps) formy związanej przez łącznie 368 ps jest równoważne zaobserwowaniu w danej chwili czasu 368 cząsteczek białka związanych z kalbindyną, w symulacji prowadzonej dla układu zawierającego 1000 cząsteczek kalbindyny i 1000 jonów wapniowych. Dlatego stosunek czasów obserwacji obu form można przyrównać do stosunku prawdopodobieństw tych stanów. Mamy $t_{zw}$ = 368 ps, $t_{nzw}$ = 1000-368=632 ps a zatem

$$\frac{P_{zw}}{P_{nzw}} = \frac{t_{zw}}{t_{nzw}} = \frac{368}{632} = \exp\left(-\frac{\Delta F}{RT}\right)$$

Stąd różnica energii swobodnej wynosi

$$\Delta F = -RT\ln\left(\frac{t_{zw}}{t_{nzw}}\right) = -8,3145\ti... ...,5823 = 1340\ \frac{\mathrm{J}}{\mathrm{mol}} = 0,32\ \mathrm{\frac{kcal}{mol}}$$

Jak widać, obliczona zmiana energi swobodnej jest dodatnia czyli bardziej prawdopodobna jest sytuacja kiedy układ jest zdysocjowany. Ten fakt wydaje się nieoczywisty ponieważ ujemnie naładowane zdeprotonowane w warunkach fizjologicznych grupy karboksylowe łańcuchów bocznych glutaminy i asparaginy znajdujących się w pętlach wiążących kalbindyny oddziałują silniej z kationami wapniowymi niż cząsteczki wody. Decydujący jest czynnik entropowy. Cząsteczek wody jest po pierwsze więcej a po drugie związanie kationu wapniowego z kalbindyną powoduje zmniejszenie liczby możliwych stanów układu, czyli spadek entropii, a co za tym idzie, wzrost energii swobodnej.

Przykład 3

Podczas reakcji nitrowania chlorobenzenu grupa nitrowa zostaje podstawiona w położeniach orto, meta i para tworząc mieszaninę, której zawartość procentowa została przedstawiona poniżej. Wiedząc, że reakcja odbywała się w temperaturze 300 K, wyznaczyć względne różnice energii potencjalnej dla otrzymanych produktów. Degeneracja w pozycji orto i meta wynosi 2, natomiast w pozycji para 1.

Rozwiązanie.

Z podanych zawartości procentowych produktów można wyliczyć ich ułamki molowe otrzymując $x_o$ = 0,30, $x_m$ = 0,01, $x_p$ = 0,69. Przyjmując formę para jako odniesienie, możemy następnie zastosować prawo rozkładu Boltzmanna do obliczenia różnicy energii pomiędzy formą orto i para ( $\Delta E_{o-p}=E_o-E_p$) oraz meta i para ( $\Delta E_{m-p}=E_m-E_p$).

$$\frac{x_o}{x_p} = \frac{g_o}{g_p}\exp\left(-\frac{\Delta E_{o-p}}{RT}\right)$$

$$\Delta E_{o-p}=-RT\ln\frac{g_px_o}{g_ox_p} = -8,3145\times 298\ti... ...{2\times 0,69} = 3806\ \mathrm{\frac{J}{mol}} = 0,91\ \mathrm{\frac{kcal}{mol}}$$

$$\Delta E_{m-p}=-RT\ln\frac{g_px_m}{g_mx_p} = -8,3145\times 298\ti... ...2\times 0,69} = 12290\ \mathrm{\frac{J}{mol}} = 2,94\ \mathrm{\frac{kcal}{mol}}$$

Stąd można również obliczyć różnicę energii pomiędzy formą orto i formą meta, $\Delta E_{m-o}=\Delta E_{m-p}-\Delta E_{o-p} = 2,94-0,91 = 2,01$ kcal/mol.

Przykład 4

Wyprowadzić prawo rozkładu Maxwella-Boltzmanna, określające rozkład gęstości prawdopodobieństwa cząsteczek gazu doskonałego w zależności od całkowitej prędkości $v$ (inaczej: od długości wektora prędkości).

Rozwiązanie.

Energia pojedynczej cząstki gazu doskonałego jest energią kinetyczną. Jeżeli prędkości w kierunku $x$, $y$ i $z$ wynoszą odpowiednio $v_x$, $v_y$ i $v_z$ a masa cząstki wynosi $m$ to jej energia kinetyczna dana jest wzorem

$$E_k = \frac{1}{2}m\left(v_x^2+v_y^2+v_z^2\right) = \frac{1}{2}mv^2$$

$$v = \sqrt{v_x^2+v_y^2+v_z^2}$$

Z powyższego wzoru wynika, że energia kinetyczna zależy tylko od całkowitej prędkości $v$ a nie zależy od kierunku wektora prędkości. Należy zatem określić ile stanów o określonej całkowitej prędkości odpowiada przedziałowi od $v$ do $v+dv$. Ilość tych stanów jest proporcjonalna do objętości sferycznej skorupy o promieniu $v$ i grubości $dv$ (taką skorupę można z kolei porównać skórki pomarańczy). Ilustracją tego jest poniższy rysunek, na którym czarne wektory oznaczają wektory prędkości o długości nie większej niż $v$ a czerwone o długości nie większej niż $v+dv$; ich ilości" wynoszą odpowiednio $N(v)$ i $N(v+dv)$ a różnica $dN=N(v+dv)-N(v)$ jest właśnie ilością wektorów prędkości o długości zawartej od $v$ do $v+dv$.

Uwaga! Formalnie zarówno $N(v)$ jak i $N(v+dv)$ są nieskończone więc takie rozważania nie są do końca ścisłe. Dokładniejsze rozważania polegają na użyciu jakobianu przekształcenia ze wpółrzędnych kartezjańskich do sferycznych i wycałkowaniu po kącie zenitalnym $\theta$ i azymutalnym $\phi$. Odpowiednie wyprowadzenia można znaleźć w rozdziale 5.4 II tomu podręcznika W. Krysicki, L. Włodarski Analiza matematyczna w zadaniach”.

Zatem degeneracja energii przy dla prędkości $v$ jest równa powierzchni skorupy pomnożonej przez jej grubość.

$$g(v) = 4\pi v^2 dv$$

Stąd prawdopodobieństwo znalezienia cząstki w przedziale prędkości od $v$ do $v+dv$ jest równe

$$P(v)dv = \frac{1}{Q}4\pi v^2\exp\left(-\frac{mv^2}{2k_BT}\right) dv$$

gdzie $Q$ jest sumą statystyczną albo inaczej czynnikiem normalizującym prawdopodobieństwo. Pozostaje wyliczyć $Q$.

$$Q = 4\pi\int\limits_{0}^{\infty}v^2\exp\left(-\frac{mv^2}{2k_BT}\right)dv$$

Aby obliczyć powyższą całkę zauważmy, że całka

$$\int\limits_{0}^{\infty}x^2\exp\left(-ax^2\right)dx$$

gdzie w naszym przypadku

$$a = \frac{m}{2k_BT}$$

jest pochodną całki z funkcji Gaussa (funkcji błędu) względem parametru $a$ ze znakiem minus.

$$\frac{d}{da}\int\limits_{0}^{\infty}\exp\left(-ax^2\right)dx= -\int\limits_{0}^{\infty}x^2\exp\left(-ax^2\right)dx$$

Jeżeli całkowanie odbywa się w granicach $\pm\infty$ to całka występująca po prawej stronie nazwa się całką Poissona, całką Laplace'a lub całką Gaussa. Wyraża się ona wzorem podanym poniżej, którego wyprowadzenie można znaleźć w W. Krysicki, L. Włodarski, Analiza matematyczna w zadaniach”, tom II, rozdz. 4.5. W naszym przypadku całkujemy tylko po połowie dziedziny a ponieważ funkcja podcałkowa jest symetryczna, całka jest równa połowie całki Gaussa.

$$\int\limits_{0}^{\infty}x^2\exp\left(-ax^2\right)dx = \frac{1}{2... ..._{-\infty}^{\infty}x^2\exp\left(-ax^2\right)dx=\frac{1}{2}\sqrt{\frac{\pi}{a}}$$

Zatem

$$-\frac{d}{da}\int\limits_{0}^{\infty}\exp\left(-ax^2\right)dx = -\frac{d}{da}\frac{1}{2}\sqrt{\frac{\pi}{a}} = \frac{1}{4}\sqrt{\frac{\pi}{a^3}}$$

co daje następujące wyrażenie na sumę statystyczną i na rozkład gęstości prawdopodobieństwa

$$Q = \sqrt{\left(\frac{2\pi k_BT}{m}\right)^3}$$

$$P(v)dv = 4\pi\left(\frac{m}{2\pi k_BT}\right)^{\frac{3}{2}}v^2\exp\left(-\frac{mv^2}{2k_BT}\right)dv$$

Wykresy gęstości rozkładu gęstości prawdopodobieństwa od całkowitej prędkości dla gazu złożonego z atomów argonu w różnych temperaturach są pokazane na poniższym rysunku.

Zadania

Zadanie 1

Pewna pochodna antrachinonu (AQ) interkaluje (wiąże się) z makromolekułami DNA, z których każda posiada 10 miejsc wiążących typu I i 40 miejsc wiążących typu II. Inaczej mówiąc, degeneracja stanu związanego I wynosi $g_{I}$ = 10 a stanu związanego II wynosi $g_{II}$ = 40. Wiedząc, że różnica pomiędzy energią wiązania AQ z miejscem II a miejscem I wynosi $\Delta E$ = $E_{AQ(II)}$ - $E_{AQ(I)}$ = 0,597 kcal/mol obliczyć ułamek molowy AQ związanego w miejscu wiążącym II w temperature $t=20^\circ C$.

Odpowiedź:

Ułamek molowy $x_{AQ(II)} = 0,59$

Zadanie 2

Pewna pochodna akrydyny (AK) interkaluje (wiąże się) z makromolekułami DNA, z których każda posiada 10 miejsc wiążących typu I i 30 miejsc wiążących typu II. Inaczej mówiąc, degeneracja stanu związanego I wynosi $g_I$ = 10 a stanu związanego II wynosi $g_{II}$ = 30. Wiedząc, że ułamek molowy AK związanej w miejscu wiążącym I w temperature $t=30^\circ$ C wynosi $x_{AK(I)} = 0,80$ obliczyć różnicę energii wiązania AK przez miejsca wiążące typu I i II.

Odpowiedź:

Różnica energii $\Delta E$ = 1,50 kcal/mol.

Zadanie 3

2-chloroetanol w dimetylosulfotlenku (DMSO) występuje w konformacji naprzemianległej ($t$) oraz skośnej ($g$), przy czym geometryczne prawdopodobieństwo konformacji skośnej jest dwa razy większe niż konformacji naprzemianległej; inaczej mówiąc degeneracje stanów $g$ i $t$ wynoszą odpowiednio $g_g$ = 2, $g_t$ = 1. W temperaturze $t=25^\circ$C stosunek ułamka molowego konformacji skośnej do konformacji naprzemianległej wynosi $x_{g}$/$x_{t}$ = 10.

a) Obliczyć różnicę energii pomiędzy konformacją skośną a naprzemianległa ($\Delta E$ = $E_{t}$ - $E_{g}$).

Odpowiedź:

Różnica energii $\Delta E$ = 0,95 kcal/mol.

b) Wyprowadzić wzór na zależność ułamka molowego konformacji naprzemianległej od temperatury.

Odpowiedź:

$$x_t(T) = \frac{\frac{1}{2}\exp\left(-\frac{\Delta E}{RT}\right)}{... ...ht)}=\frac{\exp\left(-\frac{480}{T}\right)}{2+\exp\left(-\frac{480}{T}\right)}$$

Zadanie 4

Ułamek molowy niskospinowego kompleksu Ni(II)MePy(BZ)$_2$)($t$Buacac) (degeneracja $g$ = 1) wynosi $x_{n}$ = 0,2 w temperaturze T=300 K. Obliczyć różnicę energii między kompleksem wysokospinowym (degeneracja $g$ = 3) a kompleksem niskospinowym, $\Delta E$ = $E_{w}$ - $E_{n}$.

Odpowiedź:

Różnica energii $\Delta E$ = -0,17 kcal/mol.

Zadanie 5

Jaka musi być temperatura aby 1% tlenu cząsteczkowego znalazł się w stanie singletowym? Różnica energii między stanem singletowym ( $^1\Delta_g$) a podstawowym stanem trypletowym ($^3\sum_g$) wynosi $\Delta E = 21$ kcal/mol a degeneracje tych stanów wynoszą odpowiednio $g_0$ = 3 (stan podstawowy trypletowy) i $g_1$ = 1 (stan wzbudzony singletowy).

Odpowiedź:

T = 3022 K.

Zadanie 6

Zmodyfikować wyprowadzenie rozkładu Maxwella-Boltzmanna z przykładu 4 tak aby gęstość prawdopodobieństwa była wyrażona w funkcji nie prędkości ale (a) całkowitego pędu (długości wektora pędu) cząstki i (b) energii kinetycznej cząstki gazu doskonałego. Uwaga! Należy zwrócić uwagę jak będzie się wtedy wyrażała gęstość stanów o danym pędzie lub energii i odpowiednio zmodyfikować równanie, co również będzie wymagało modyfikacji lub ponownego wyprowadzenia wzoru na sumę statystyczną (Q). Pęd oraz energia kinetyczna w funkcji pędu wyrażają się następującymi wzorami:

$$\mathbf{p}=m\mathbf{v}$$

$$p = m\sqrt{v_x^2+v_y^2+v_z^2}$$

$$E_k = \frac{p_x^2+p_y^2+p_z^2}{2m}$$

Zadanie 7*

Wyprowadzić wzór na gęstość prawdopodobieństwa znalezienia się oscylatora harmonicznego charakteryzującego się masą zredukowaną $\mu$ oraz stałą siłową $k$ w zależności od współrzędnej normalnej $x$ i pędu $p$. Energia oscylatora wyraża się wzorem

$$E = E_k+E_p = \frac{p^2}{2\mu}+\frac{1}{2}kx^2$$