【統計力学も活用】水の自由エネルギーを量子化学計算ソフトで調べよう

自由エネルギーを量子化学計算ソフトで調べよう

分子の反応性や安定性を調べるためにはギブズエネルギーなどの自由エネルギーの情報は欠かせません。

じつは量子化学計算パッケージを使うことで、分子軌道だけではなく、自由エネルギーを簡便に調べることができます。

このページではフリーソフトであるGAMESSの計算結果を使って、分子の自由エネルギーがどのようにして得られるかを紹介します。

このページで紹介すること
  • 自由エネルギーへは電子エネルギーの寄与が最も大きい
  • 自由エネルギーへの補正値は統計力学の知見を使うと得られる
  • 並進運動は理想気体を仮定する
  • 回転運動は対称性の考慮が重要
  • 振動運動は全ての振動モードを取り扱う
目次

自由エネルギーの寄与

このページは、量子化学計算パッケージのGAMESSを使った水分子の計算結果をもとに作成しています。計算ログは下からダウンロードできます。(量子化学計算条件の詳細はログを確認願います。)

統計力学を使えば自由エネルギーがわかる

まず重要な知識ですが、自由エネルギーには電子エネルギーがほとんどの寄与を占めます。

例えば基底状態の水分子に対して量子化学計算を行うと、電子状態のエネルギー値として$-76.025\text{Hartree}=-199.60\text{GJ/mol}$が得られます。電子状態のエネルギーはとても大きい数値なので、自由エネルギーの値もほとんどこの値に近くなります。

   -----------------
   ENERGY COMPONENTS
   -----------------
   (中略)
  TOTAL ENERGY =     -76.0257177698

電子状態のエネルギーが大きいため、化学反応熱が大きくなります。

自由エネルギーはこの電子エネルギーに補正を加えることで得ることができます。

ギブズエネルギー$G$と内部エネルギー$U$には、$G=U+PV-TS$の関係があります。このうち、GAMESSでは、$PV$については理想気体を仮定することで、$G=U+RT-TS$によってギブズエネルギーを求めます。

GAMESSのデフォルト設定では、温度$298.15\text{K}$、大気圧の理想気体として自由エネルギーを計算しています。

つまり、内部エネルギー$U$とエントロピー$S$を計算できればギブズエネルギー$G$を計算できるね。

では、どのようにして$U$と$S$を得るのでしょうか?

答えは統計力学にあります。

統計力学に基づいた分配関数$Z$を計算できれば、内部エネルギーは$U=-\frac{\partial}{\partial\beta}\ln{Z}$、エントロピーは$S=\frac{U}{T}+k_B\ln{Z}$を用いて計算ができるのです。

               Q               LN Q
 ELEC.     1.00000E+00       0.000000
 TRANS.    3.00431E+06      14.915558
 ROT.      4.02909E+01       3.696126
 VIB.      1.00029E+00       0.000289
 TOT.      1.21081E+08      18.611973

これはGAMESSによる水分子の計算結果です。Qが分配関数$Z$で、LN Qが$\ln{Z}$に対応しています。

通常、電子エネルギーでは熱励起はほとんどないので$U_\text{elec}$と$S_\text{elec}$は無視できます

自由エネルギーは並進・回転・振動に分けて考える

では、どのようにして分子の分配関数$Z$を計算するのでしょうか?

「複雑なものは分けて考えよう」ということで、並進運動・回転運動・振動運動のそれぞれの寄与を合算することで、分配関数や自由エネルギーを計算することができます。

例えば$U_\text{total}=U_\text{trans}+U_\text{rot}+U_\text{vib}$ということだね。

統計力学を利用すると分子の並進、回転、振動の分配関数やエントロピーは次の式を利用することで計算できます。

並進・回転・振動の分配関数とエネルギー、エントロピー

GAMESSでも各寄与ごとに足しあわされて内部エネルギーやギブズエネルギーなどが計算されています。(下は水分子の例)

              E         H         G         CV        CP        S
           KJ/MOL    KJ/MOL    KJ/MOL   J/MOL-K   J/MOL-K   J/MOL-K
 ELEC.      0.000     0.000     0.000     0.000     0.000     0.000
 TRANS.     3.718     6.197   -36.975    12.472    20.786   144.800
 ROT.       3.718     3.718    -9.162    12.472    12.472    43.203
 VIB.      59.288    59.288    59.281     0.160     0.160     0.022
 TOTAL     66.725    69.204    13.144    25.103    33.417   188.025

例えば、$G_{\text{trans}}=U_{\text{trans}}+RT-TS_{\text{trans}}$$=3.718+\frac{8.3144}{1000}*298.15-298.15*\frac{144.8}{1000}$$=-36.9748$と計算されています。

以下では、熱力学量に対する並進、回転、振動による寄与についてそれぞれ解説していきます。

並進運動による寄与

背景

並進運動はエネルギーとしては大きな値は持ちませんが、エントロピーには大きな寄与があります。これは並進運動はほとんど連続したエネルギー準位を持つことに寄ります。

質量$M(\text{kg})$の分子($N_A\text{個}$)による並進運動の分配関数は$Z_{\text{trans}}=\frac{V^{N_A}}{N_A!}\Bigl(\frac{2\pi{M}k_BT}{h^2}\Bigl)^{\frac{3}{2}{N_A}}$と書くことができます。

$V$は大気圧下での理想気体の値を採用するので、$V=\frac{P}{R}=4.06261\times10^{-26}(\text{m}^3)$です

水でもベンゼンでも$V=4.06261\times10^{-26}(\text{m}^3)$とするんだね

内部エネルギーについては、理想気体として取り扱うので、$U=\frac{3}{2}RT$で計算できます。

エントロピーは分配関数を利用することで、$S=\frac{U}{T}+k_B\ln{Z}$から、$S=R\Bigl(\frac{3}{2}\ln(\frac{2\pi{M}k_BT}{h^2})+\ln{V}+\frac{5}{2}\Bigl)$とできます。これはザックール・テトローデ方程式と呼ばれる式でもあり、分子の重さだけに依存します。

実際の計算例

GAMESSによる水分子の計算結果を見てみましょう。

              E         H         G         CV        CP        S
           KJ/MOL    KJ/MOL    KJ/MOL   J/MOL-K   J/MOL-K   J/MOL-K
 ELEC.      0.000     0.000     0.000     0.000     0.000     0.000
 TRANS.     3.718     6.197   -36.975    12.472    20.786   144.800
 ROT.       3.718     3.718    -9.162    12.472    12.472    43.203
 VIB.      59.288    59.288    59.281     0.160     0.160     0.022
 TOTAL     66.725    69.204    13.144    25.103    33.417   188.025

まず、並進は「TRANS.」で表されており、$E=3.718(\text{KJ/MOL})$となっていますが、これは$U_{\text{trans}}=\frac{3}{2}RT$に$R=\frac{8.31}{1000}$、$T=298.15$を代入した数値と一致します。

また、エントロピーについても、$S=R\Bigl(\frac{3}{2}\ln(\frac{2\pi{M}k_BT}{h^2})+\ln{V}+\frac{5}{2}\Bigl)$に$M=18\times{M_0}$、$V=\frac{P}{R}=4.06261\times10^{-26}(\text{m}^3)$などを代入した数値と一致していることがわかります。

$M_0$は原子質量単位で$1.66054×10^{-27}(\text{kg})$だね

単位換算がややこしければSI単位系で常に考えましょう

この得られた$U$と$S$を$H=U+RT$、$G=H-TS$へ代入することでエンタルピー$H$とギブズエネルギー$G$が得られています。

回転運動による寄与

背景

3つの慣性モーメント$I_A,I_B,I_C(\text{kg}\cdot\text{m}^2)$を持つ分子の回転による内部エネルギーとエントロピーの寄与を考えましょう。

まず、内部エネルギーについては、並進運動と同様にエネルギー等分配の法則から、$U_{\text{rot}}=\frac{3}{2}RT$で与えることができると考えます。

一方、エントロピーについては分配関数として$Z_{\text{rot}}=\frac{1}{\sigma}\Bigl[\pi^{\frac{1}{2}}\Bigl(\frac{2I_Ak_BT}{\hbar^2}\Bigl)^{\frac{1}{2}}$$\Bigl(\frac{2I_Bk_BT}{\hbar^2}\Bigl)^{\frac{1}{2}}\Bigl(\frac{2I_Ck_BT}{\hbar^2}\Bigl)^{\frac{1}{2}}\Bigl]^{N_A}$を利用します。$\sigma$は分子の回転対称数です。

例えば水分子だと$\sigma=2$です。

慣性モーメント$I_A(\text{kg}\cdot\text{m}^2)$の代わりに回転定数$A=\frac{\hbar}{4\pi{I_A}}(\text{Hz})$を導入して分配関数を書き換えると$Z_{\text{rot}}=\frac{1}{\sigma}\Bigl(\frac{\pi}{ABC}\Bigl)^{\frac{N_A}{2}}\Bigl(\frac{k_BT}{h}\Bigl)^{\frac{3N_A}{2}}$と分配関数を簡素に書くことができます。

これを利用すると、エントロピーは$S_{\text{rot}}=R\Bigl[\frac{3}{2}+\frac{1}{2}\ln\frac{\pi}{ABC}+\frac{3}{2}\ln\frac{k_BT}{h}-\ln\sigma\Bigl]$によって計算することができます。

実際の計算例

GAMESSによる水分子の計算結果を見てみましょう。

              E         H         G         CV        CP        S
           KJ/MOL    KJ/MOL    KJ/MOL   J/MOL-K   J/MOL-K   J/MOL-K
 ELEC.      0.000     0.000     0.000     0.000     0.000     0.000
 TRANS.     3.718     6.197   -36.975    12.472    20.786   144.800
 ROT.       3.718     3.718    -9.162    12.472    12.472    43.203
 VIB.      59.288    59.288    59.281     0.160     0.160     0.022
 TOTAL     66.725    69.204    13.144    25.103    33.417   188.025

回転は「ROT.」で表されており、$E=3.718(\text{KJ/MOL})$となっていますが、これは並進運動と同様に$U_{\text{rot}}=\frac{3}{2}RT$に$R=8.31$、$T=298.15$を代入した数値と一致します。

次にエントロピーについてですが、$S_{\text{rot}}=R\Bigl[\frac{3}{2}+\frac{1}{2}\ln\frac{\pi}{ABC}+\frac{3}{2}\ln\frac{k_BT}{h}-\ln\sigma\Bigl]$に回転定数$A,B,C$を代入する必要がありますが、これもGAMESSで出力されています。

 THE ROTATIONAL CONSTANTS ARE (IN GHZ)
    988.63138   406.32750   287.97130

つまり、$A=988*10^9,B=406*10^9,$$C=287*10^9(\text{Hz})$を代入し、水分子の対称性から$\sigma=2$であることに注意することでエントロピーが得られ、これはGAMESSの出力結果$S_{\text{rot}}=43.2(\text{J/mol}\cdot\text{K})$と一致します。

振動運動による寄与

背景

内部エネルギーに対する振動モードによる寄与は並進・回転に比べるとかなり大きくなります。それはゼロ点振動によって書く振動モードが必ず一定以上のエネルギーを持つためです。

振動運動のエネルギー準位幅は並進・回転に比べて大きいため、内部エネルギーについてエネルギー等分配の法則を利用して$U=\frac{\bigcirc}{\blacksquare}RT$のような形では与えることはできません。そのため、分配関数を利用して計算する必要があります。

$N$個の振動モード$\nu_1,\cdots,\nu_N(\text{Hz})$がある場合、分配関数は$Z_{\text{vib}}=\Pi_{i=1}^N\sum_{n=0}^\infty\exp[-\beta{h\nu_i}(n+\frac{1}{2})]$$=\Pi_{i=1}^N\frac{1}{2\sinh\frac{\beta{h\nu_i}}{2}}$で与えられ、これを利用すると内部エネルギーは$U_{\text{vib}}=\sum_ih\nu_i[\frac{1}{2}+\frac{1}{e^{\beta{h}\nu_i}-1}]$とできます。

ここで、振動数$\nu(\text{Hz})$ではなく波数$k=\frac{\nu}{c}(\text{m})$を利用すると、$U_{\text{vib}}=RT\sum_i\frac{\theta_i}{T}[\frac{1}{2}+\frac{1}{e^{\theta_i/T}-1}]$とも書けます。($\frac{\theta}{T}=\frac{hkc}{k_BT}$とおいています)

エントロピーについては、$S_{\text{vib}}=R\sum_i\Bigl[\frac{\theta_i}{T}\frac{1}{e^{\theta_i/T}-1}-\ln(1-e^{-\theta_i/T})\Bigl]$で計算できます。

実際の計算例

GAMESSによる水分子の計算結果を見てみましょう。

              E         H         G         CV        CP        S
           KJ/MOL    KJ/MOL    KJ/MOL   J/MOL-K   J/MOL-K   J/MOL-K
 ELEC.      0.000     0.000     0.000     0.000     0.000     0.000
 TRANS.     3.718     6.197   -36.975    12.472    20.786   144.800
 ROT.       3.718     3.718    -9.162    12.472    12.472    43.203
 VIB.      59.288    59.288    59.281     0.160     0.160     0.022
 TOTAL     66.725    69.204    13.144    25.103    33.417   188.025

振動は「VIB.」で表されており、$E_{\text{vib}}=59.288(\text{KJ/MOL})$と、並進、回転とは異なる数値です。これを得るためにはすべての振動モードを知る必要があります。

$N$個の原子で構成される非直線分子には$3N-6$個の振動モードが発生します。

水分子だと$3N-6=3$個だね

水の振動モードはGAMESSによってすべて出力されており、次の通りです。

  MODE FREQ(CM**-1)  SYMMETRY  RED. MASS  IR INTENS.
           (中略)
     7    1688.567    A1       1.086338    2.259042
     8    4033.275    A1       1.041776    0.471534
     9    4189.330    B1       1.087619    2.144479

振動モードの波数$k$に対して$1688,4033,4189(\text{cm}^{-1}$)を$U_{\text{vib}}=RT\sum_i\frac{\theta_i}{T}[\frac{1}{2}+\frac{1}{e^{\theta_i/T}-1}]$($\frac{\theta}{T}=\frac{hkc}{k_BT}$)へ代入すると、これもGAMESSの結果$U_{\text{vib}}=59(\text{KJ/MOL})$と一致します。

最後に振動のエントロピーについても$S_{\text{vib}}=R\sum_i\Bigl[\frac{\theta_i}{T}\frac{1}{e^{\theta_i/T}-1}-\ln(1-e^{-\theta_i/T})\Bigl]$を各振動モードで足し合わせるとGAMESSの結果$S_{\text{vib}}=0.022(\text{J/K}\cdot\text{MOL})$と一致することがわかります。

まとめ

このページではフリーソフトであるGAMESSの計算結果を参考にしながら、分子の自由エネルギーがどのようにして得られるかを紹介しました。

分子の反応性や安定性を調べるためにはギブズエネルギーなどの自由エネルギーの情報は欠かせません。

この自由エネルギーは量子化学計算パッケージを使うことで簡便に調べることができます。また、これらの自由エネルギーの計算結果が統計力学の知見を使っていることを知っていただければと思います。

目次