こんにちは。
今回は「ヒュッケル法」を python で計算する、という題でやっていきます。
おそらくこの記事を読んでいる方はヒュッケル法をある程度知っていてその計算を機械にやらせたい、という方だと思いますので原理の説明はあまり詳しくはしないことにします。
知らない方も、結論だけでも一応実行はできると思うのでなんとなく見てみてください。
目次
ヒュッケル法とは
概要
気づいている方もいるかもしれませんが、筆者は理系の中でも化学寄りの人間です。
中でも有機化学が好きで、これを学んでいると必ず出会うであろうものが「ヒュッケル法」です。
ヒュッケル法は、共役系の炭化水素について分子軌道関数を計算する方法で、分子軌道のエネルギー準位や電子密度が求められ、反応がどこで起こりやすいかなどがわかったりします。
（分子軌道はざっくりいうと「電子があるところ」。自分もはっきりと言語化できるほど理解ができておらず、「orbital」という表現くらい曖昧な感じなのですが…）
分子軌道という概念はディールス-アルダー反応の機構との関連を知ったときにイメージが掴めたような気がしたので、ぜひ調べてみてください。
「拡張ヒュッケル法」というものもあるようですが、今回は一番簡単であろう、いわゆる「単純ヒュッケル法」について書いていきます。
原理
原理の細かい説明は他のサイトを参照してください。
例えばこちら。
実際の手順としては、
$$ \psi = \sum_{r = 1}^{n} c_{r} \phi_{r} $$
と表されるとした分子軌道関数 $\psi$ を考えます。
永年行列式
$$ \begin{vmatrix} H_{11} - E S_{11} & H_{12} - E S_{12} & \dots & H_{1n} - E S_{1n} \\ H_{21} - E S_{21} & H_{22} - E S_{22} & \dots & H_{2n} - E S_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ H_{n1} - E S_{n1} & H_{n2} - E S_{n2} & \dots & H_{nn} - E S_{nn} \end{vmatrix} = 0 $$
をつくり、これを解いて分子軌道のエネルギー $E$ を求めます。
ここで、ヒュッケル法においては、
$S_{ij}$ は $i = j$ で 1、$i \neq j$ で 0
$H_{ij}$ は $i = j$ で $\alpha$、$i \neq j$ で結合があれば $\beta$、無ければ 0
と単純化して計算する点が肝です。
さらに、上の連立方程式と、規格化条件
から係数 $c_{r}$ を計算することで、分子軌道関数 $\psi$ がその分子を構成する原子 (水素以外) の数だけ求められます。
永年方程式を考えてみるとわかるかと思いますが、6炭素（ベンゼンなど）でも割と大変で、簡単に分子軌道の計算ができる方ではあるんですが、そうはいっても手計算でやるには時間がかかります。
そこをコンピュータにお任せして、簡単に計算できるようにしようというのが今回のねらいです。
手計算でやってみる
まずは流れを理解するためにも、最も単純な共役系であるエチレンについて計算してみます。
エチレンは「H2C=CH2」の形で表される分子で、便宜上、左の炭素を1、右の炭素を2と番号付けしておきます。
各原子の原子軌道関数 $\phi$ の線形結合から、分子軌道関数 $\psi$ は、
$$ \psi = c_{1} \phi_{1} + c_{2} \phi_{2} $$
と表される、ということにします。
また、炭素1と炭素2の間には結合があるため、永年行列式は、
$$ \begin{vmatrix} \alpha - E & \beta \\ \beta & \alpha - E \end{vmatrix} = 0 $$
となります。
一般的には、両辺を $\beta$ で割り、$\displaystyle \frac{ \alpha - E }{ \beta } = X$ と置いた、
$$ \begin{vmatrix} X & 1 \\ 1 & X \end{vmatrix} = 0 $$
という方程式を計算することが多いです。
（また、永年行列式をつくる際には、対角部分を $X$、それ以外の部分は結合があれば 1、なければ 0 とすればよいことになります。）
これを解くと、
$$ X^{2} - 1 = 0 \Leftrightarrow X^{2} = 1 \Leftrightarrow X = \pm 1 $$
となることから、
$$ \frac{ \alpha - E }{ \beta } = \pm 1 \Leftrightarrow E = \alpha \mp \beta $$
と、分子軌道2つのエネルギーが求められました。
さらに、永年行列式から $X c_{1} + c_{2} = 0$ と $c_{1} + X c_{2} = 0$ という方程式が得られ、これと規格化条件 $c_{1}^{2} + c_{2}^{2}= 1$ から、$X$ の値によって場合分けをしつつ、$c_{1}$ と $c_{2}$ の値を計算します。
すると、
$X = -1$ のとき、$\displaystyle c_{1} = c_{2} = \frac{ \sqrt{2} }{ 2 }$、
$X = 1$ のとき、$\displaystyle c_{1} = \frac{ \sqrt{2} }{ 2 }, c_{2} = - \frac{ \sqrt{2} }{ 2 }$、
と計算できるので、求められた分子軌道関数 $\psi$ は、
$$ \begin{align} \psi_{1} &= \frac{ \sqrt{2} }{ 2 } \phi_{1} + \frac{ \sqrt{2} }{ 2 } \phi_{2} \\ \psi_{2} &= \frac{ \sqrt{2} }{ 2 } \phi_{1} - \frac{ \sqrt{2} }{ 2 } \phi_{2} \end{align} $$
の2つになります。
ここから、炭素原子2つからエチレン分子をつくったとしたとき、
というエネルギー図が描け、炭素原子が単体でいるより、エチレン分子となった方が $2 \beta$ 分だけ安定化することができる、といえます。
（こんなようなものは、どこかで見たことがある方もいるかもしれません。）
$\psi_{1}$ や $\psi_{2}$ の右隣にあるのが、求めた分子軌道関数から分子軌道の形を模式的に表したものです。
プログラム化
方針
さて、ヒュッケル法がなんとなくわかったところで、プログラムにしてみます。
結局は、例えばエチレンでは、
$$ \begin{vmatrix} X & 1 \\ 1 & X \end{vmatrix} = 0 $$
というような永年方程式を解きたいので、これをどうするかが大きな問題となります。
2行2列や3行3列ならサラスの方法やなんかが使えますし、より大きな行列式でも余因子展開などの方法があります。
ですが高次の方程式となるので、解くのが大変です。
（二次方程式なら解の公式がありますが、それ以上は無いか、あっても複雑です。）
また、永年行列式が解けてエネルギー $E$ が求められたとしても、続いて係数 $c$ を計算するためにまた連立方程式を解かなければなりません。
というわけで、別の方法で計算します。
エチレンにおける永年行列式は
$$ \left( X \left( \begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix} \right) - \left( \begin{matrix} 0 & -1 \\ -1 & 0 \end{matrix} \right) \right) \left( \begin{matrix} c_{1} \\ c_{2} \end{matrix} \right) = 0 $$
と変形できるので、$X$ を固有値、係数 $c$ は $\left( \begin{matrix} c_{1} \\ c_{2} \end{matrix} \right)$ と表せて固有ベクトルとみなせます。
（固有値、固有ベクトルの話は自分もちゃんとは理解できていないので詳しくは説明できません…）
そしてなんと、pythonには「numpy」というライブラリがあり、そこには固有値、固有ベクトルを計算してくれる関数があるので、これを使えば数行で完結します。
ライブラリのインストール
ライブラリを使うので一応インストール方法も書いておきます。
pip install numpy
これを実行して待てば完了です。
呼び出すときは、
import numpy as np # npの部分は任意に変更可能 (このあとの呼び名を決められる)
と書いておけばOKです。
固有値・固有ベクトルの計算
numpyに存在する関数 eigh() は、引数に行列を入れるとその固有値や固有ベクトルを計算してくれます。
先程変形した永年行列式から、$\displaystyle \left( \begin{matrix} 0 & -1 \\ -1 & 0 \end{matrix} \right)$ が引数に入れるべき行列であり、元の永年行列式でいうところの、対角部分を 0 にして全体を-1倍したものになります。
プログラム的には、
import numpy as np # ライブラリの読み込み A = np.zeros((2, 2)) # 2行2列の行列を作成 A[0][1] = -1 # 炭素0と炭素1が結合している A += A.T # Aを転置したものを足す # print(A) X, c = np.linalg.eigh(A) # 固有値・固有ベクトルの計算 print(f"X = {X}") print(f"c = {c}")
これだけです。簡単。
実行すると、X = [-1. 1.]、c = [[-0.70710678 -0.70710678] [-0.70710678 0.70710678]] が出力されました。
先程の計算と合っていそうです。
ここで注意なのですが、固有ベクトルは (c[0][0], c[0][1]) や (c[1][0], c[1][1]) ではなく、((c[0][0], c[1][0])) や ((c[0][1], c[1][1])) です。
要は縦に並んだものをセットとして見るのが正解ですよ、ということです。
（c.T で転置すると見やすくなると思います。）
以上から、分子軌道関数とそのエネルギーは
$$ \begin{align} \psi_{1} &= 0.707 \phi_{1} + 0.707 \phi_{2} ~~~~~ &E_{1} = \alpha + \beta \\ \psi_{2} &= 0.707 \phi_{1} - 0.707 \phi_{2} ~~~~~ &E_{2} = \alpha - \beta \end{align} $$
と計算できました。
$\displaystyle \frac{ \sqrt{2} }{ 2 } = 0.707$ なので合っています。
試運転
せっかくなのでブタジエンについても計算してみます。
import numpy as np A = np.zeros((4, 4)) A[0][1] = -1 A[1][2] = -1 A[2][3] = -1 A += A.T print(A) X, c = np.linalg.eigh(A) print(f"X = {X}") print(f"c = {c}")
から、分子軌道関数とそのエネルギーは、
$$ \begin{align} \psi_{1} = 0.372 \phi_{1} + 0.602 \phi_{2} + 0.602 \phi_{3} + 0.372 \phi_{4} ~~~~~ &E_{1} = \alpha + 1.618 \beta \\ \psi_{2} = 0.602 \phi_{1} + 0.372 \phi_{2} - 0.372 \phi_{3} - 0.602 \phi_{4} ~~~~~ &E_{2} = \alpha + 0.618 \beta \\ \psi_{3} = 0.602 \phi_{1} - 0.372 \phi_{2} - 0.372 \phi_{3} + 0.602 \phi_{4} ~~~~~ &E_{3} = \alpha - 0.618 \beta \\ \psi_{4} = 0.372 \phi_{1} - 0.602 \phi_{2} + 0.602 \phi_{3} - 0.372 \phi_{4} ~~~~~ &E_{4} = \alpha - 1.618 \beta \end{align} $$
と計算され、正しい計算結果となりました。
π電子は4個あるのでブタジエン分子となることによって $4.472 \beta$ だけ安定化できます。
また、同じπ電子数のエチレン分子2個分と比べると $0.472 \beta$ だけエネルギーが低く、これが共役系が広がることによる安定化であるといえます。
その他、ベンゼンとシクロブタジエンの安定性の違いや、ナフタレンやアントラセンの反応しやすい点、共役系の長さによる色の違いなど、いろいろと面白いことがわかるのですが、それはまた別の機会に。
追記（2022/12/10）
ベンゼンなど、環状の分子を計算しようとするとずれることがわかりました。
ベンゼンの分子軌道とそのエネルギーは、
$$ \begin{align} \psi_{1} = 0.408 \phi_{1} + 0.408 \phi_{2} + 0.408 \phi_{3} + 0.408 \phi_{4} + 0.408 \phi_{5} + 0.408 \phi_{6} ~~~ &E_{1} = \alpha + 2 \beta \\ \psi_{2} = 0.577 \phi_{1} + 0.289 \phi_{2} - 0.289 \phi_{3} - 0.577 \phi_{4} - 0.289 \phi_{5} + 0.289 \phi_{6} ~~~ &E_{2} = \alpha + \beta \\ \psi_{3} = 0.000 \phi_{1} - 0.500 \phi_{2} - 0.500 \phi_{3} - 0.000 \phi_{4} + 0.500 \phi_{5} + 0.500 \phi_{6} ~~~ &E_{3} = \alpha + \beta \\ \psi_{4} = 0.000 \phi_{1} - 0.500 \phi_{2} + 0.500 \phi_{3} - 0.000 \phi_{4} - 0.500 \phi_{5} + 0.500 \phi_{6} ~~~ &E_{4} = \alpha - \beta \\ \psi_{5} = 0.577 \phi_{1} - 0.289 \phi_{2} - 0.289 \phi_{3} + 0.577 \phi_{4} - 0.289 \phi_{5} - 0.289 \phi_{6} ~~~ &E_{5} = \alpha - \beta \\ \psi_{6} = 0.408 \phi_{1} - 0.408 \phi_{2} + 0.408 \phi_{3} - 0.408 \phi_{4} + 0.408 \phi_{5} - 0.408 \phi_{6} ~~~ &E_{6} = \alpha - 2 \beta \\ \end{align} $$
と表されます。
しかし、numpy.linalg.eigh() を使って計算すると、係数 c について、
[[4.08248290e-01 -5.75580380e-01 4.51725513e-02 5.77350269e-01 -6.72121529e-18 -4.08248290e-01]
[ 4.08248290e-01 -3.26910767e-01 -4.75880955e-01 -2.88675135e-01 -5.00000000e-01 4.08248290e-01]
[ 4.08248290e-01 2.48669613e-01 -5.21053507e-01 -2.88675135e-01 5.00000000e-01 -4.08248290e-01]
[ 4.08248290e-01 5.75580380e-01 -4.51725513e-02 5.77350269e-01 -1.09414667e-18 4.08248290e-01]
[ 4.08248290e-01 3.26910767e-01 4.75880955e-01 -2.88675135e-01 -5.00000000e-01 -4.08248290e-01]
[ 4.08248290e-01 -2.48669613e-01 5.21053507e-01 -2.88675135e-01 5.00000000e-01 4.08248290e-01]]
と出力され、合いません。
X の値は合うので、固有ベクトルの計算方法の問題だと思いますが、何が悪いのかわかりません…
さて、固有値・固有ベクトルの計算は、numpy だけでしかできないわけではありません。
「scipy」というライブラリもあり、numpy と同じように scipy.linalg.eigh() という関数があります。
試しにと思いやってみると、
from numpy import zeros from scipy import linalg n = 6 A = zeros((n, n)) for i in range(0, n-1): A[i][i+1] = -1 A[0][n-1] = -1 A += A.T print(A) X, c = linalg.eigh(A) print(f"X = {X}") print(f"c = {c}")
[[-4.08248290e-01 5.77350269e-01 0.00000000e+00 0.00000000e+00 5.77350269e-01 4.08248290e-01]
[-4.08248290e-01 2.88675135e-01 -5.00000000e-01 -5.00000000e-01 -2.88675135e-01 -4.08248290e-01]
[-4.08248290e-01 -2.88675135e-01 -5.00000000e-01 5.00000000e-01 -2.88675135e-01 4.08248290e-01]
[-4.08248290e-01 -5.77350269e-01 -1.11022302e-16 -1.11022302e-16 5.77350269e-01 -4.08248290e-01]
[-4.08248290e-01 -2.88675135e-01 5.00000000e-01 -5.00000000e-01 -2.88675135e-01 4.08248290e-01]
[-4.08248290e-01 2.88675135e-01 5.00000000e-01 5.00000000e-01 -2.88675135e-01 -4.08248290e-01]]
合いました。
何故かはわかりません。
エチレンやブタジエンについても計算結果が合ったので、固有値・固有ベクトルの計算は scipy で行った方が良いと思われます。
まとめ
以上で、ヒュッケル法をpythonで計算できるようになりました。
今回はライブラリをガンガン使って簡単にやってみましたが、そのうち計算部分も自分で書いて、それらしく描画もできるようにして、最終的にはprocessingでも実装できるようにしたいと思っています。
また、拡張ヒュッケル法も実装できれば良いなと思っています。
追記（2022/12/12）
分子軌道の描画編も記事にしました。
以上
