日曜化学(3)：分子軌道法と可視化（Python/matplotlib）

いよいよ分子軌道を計算してみたいと思います。

今回の記事の内容を理解するとエチレンやブタジエンやベンゼンなどの分子軌道が計算でき、それをPythonのプログラムで可視化できるようになります。

これまで３回に渡って書いてきた「日曜化学シリーズ」の記事ですが、今回がまさに集大成となっています。

過去の記事を前提にお話しますので、まだの方はシリーズの過去記事をご覧になってください。
tsujimotter.hatenablog.com

（番外編の日曜化学(2.5)は読まなくても、今回の内容については大丈夫です。）

前回までの記事で計算したのは、水素様原子という原子核が１つ・電子が１つのものでした。

そうなると、原子核が２つ以上で電子が１つの状況（つまり分子）を計算したくなると思います。

上記の状況はポテンシャルによって表すことができますので、ハミルトニアンに反映させればシュレーディンガー方程式を立式すること自体は可能です。（これはあとでやりたいと思います。）

ところがどっこい、この状況のシュレーディンガー方程式を解こうと思うと、もはや厳密には解けなくなってしまうのです。

ここで「量子力学はこんなものなのか」とがっかりしないでください。近似的にであれば、実用的には十分解くことができるというのが今回のお話です。

そのための方法が分子軌道法です。

今回の記事では、分子軌道法の基本的な原理の紹介と、エチレンやベンゼンなどのいくつかの分子について、電子軌道を可視化してみたいと思います！

これまで通り、可視化に用いたPythonプログラムも紹介するので、ぜひ遊んでみてください。とても楽しいです！

目次：

1. 分子軌道法の考え方（概要）

2. 重ね合わせの原理

3. 分子軌道法の計算

4. 水素分子のσ結合

5. エチレン分子のπ結合

6. ヒュッケル法とブタジエンのπ電子共役系

7. ベンゼンのπ共有系

8. 分子軌道の可視化のためのPythonプログラム

9. おわりに

最終回はこちら！

参考文献

実はこの話、私が大学１年に入った頃に習った内容なのですが、当時はさっぱりわからなかったのを記憶しています。なんとなくよくわからない前提を紹介されて、そこから突然「安定性軌道」とか「反安定性軌道」といった、分子軌道法の様々な概念が紹介されていきました。私はいったい何に立脚して考えれば良いのか、よくわかりませんでした。

それもそのはず、本当は量子力学を知らないといけなかったのです。量子力学の基本的なことを勉強したあとで改めて読んでみると「なんだこんなことだったのか」と。

やっぱり基本的な前提があやふやな状態では、理解にはたどり着けないということですね。

ここまで読んできた方であれば十分ついていける内容です。かなりの長文記事ですが、丁寧な解説を心掛けたつもりですので、よろしければ最後まで読んでいただければ幸いです。

1. 分子軌道法の考え方（概要）

まずは大まかに分子軌道法の考え方について紹介します。

最終的な目標は、２つ以上の原子核で構成された分子に属する電子を考えて、その真の波動関数 $\Psi_{\text{真}}$ を求めることです。

このために分子の状況を表すハミルトニアン $\hat{H}_{\text{分子}}$ を考えます。ハミルトニアン $\hat{H}_{\text{分子}}$ に対応するシュレーディンガー方程式

$\hat{H}_{\text{分子}} \Psi = E \Psi$

を立式できて、これは固有値問題になっています。これを解くことができれば、分子内の電子が取りうる真のエネルギー $E_{\text{真}}$ が求まり、対応する真の波動関数 $\Psi_{\text{真}}$ が得られるはずです。

ところがです。このような状況のシュレーディンガー方程式は、複雑すぎて厳密には解けないことが知られています。さぁ困った！！

ここで必要になるのは、解ける問題にどうにかして帰着できないかという考え方です。我々が解くことができるのは、水素様原子（原子核１個・電子１個）の場合だけです。

たとえば２つの原子核の状況を考えて、原子核2がとーっても遠い位置にあるという状況を考えましょう。

このような状況では、電子にとっての状況は、原子核1の周りだけを回っているのとそう変わりません。つまり、原子核1についての、水素様原子の解を考えればよいわけです。電子が原子核2の周囲にいる場合も同様です。

そんなわけで、 $R$ が十分大きい場合においては、原子核1の周りでは電子が（原子核1についての）水素様原子の波動関数 $\chi_1$ に近似でき、原子核2の周りでは同じく（原子核2についての）水素様原子の波動関数 $\chi_2$ に近似できそうです。

もちろん、 $R$ が小さい場合にはそうはいきません。しかしながら、そうはいっても上記の状況と大きく変わることはないだろうという「仮定」を置くのです。なかなか大胆な仮定ですね。

ここで、二つの既知の波動関数 $\chi_1, \chi_2$ の線形結合

$c_1 \chi_1 + c_2 \chi_2$

を考えます。波動関数の線形結合のことを量子力学では「重ね合わせ」と呼びますが、文字通り電子の波が重なっているようなイメージです。

分子における真の波動関数 $\Psi_{真}$ は、このような重ね合わせによって近似できると考えるのです。もう少し正確にいえば、真の波動関数の近似

$\Psi_{真} \fallingdotseq c_1 \chi_1 + c_2 \chi_2$

になるような $(c_1, c_2)$ を探す問題だと捉えなおすわけです！　見事な発想の転換ですね！！

ここで「変分原理」と呼ばれる性質を活用します。 $(c_1, c_2)$ に対応する波動関数 $\Psi$ のエネルギーを $E(c_1, c_2)$ とすると、真のエネルギー $E_{\text{真}}$ に対して常に

$E_{\text{真}} \leq E(c_1, c_2)$

が成り立つというものです。

したがって、エネルギー $E(c_1, c_2)$ が最小になるような $(c_1, c_2)$ を求めれば、それが真のエネルギーの近似値になっていると考えることができそうです！　このようにして得られた軌道を「分子軌道」といいます。

このような方法により近似的に分子軌道を求める手法を、量子力学では「変分法」と呼ぶそうです。

以上が分子軌道法のアイデアです。いやー、アイデアだけでもめちゃめちゃ面白いですね！！

分子軌道法を具体的に理解するためには、まずは「重ね合わせの原理」について理解する必要があります。重ね合わせの原理と、重ね合わせの状況におけるエネルギーの計算方法について、次節で説明してから分子軌道法の具体的な計算にいきたいと思います。

2. 重ね合わせの原理

量子力学においては、電子などのミクロな粒子は「異なる量子状態の重ね合わせ」として存在するのだという話を聞いたことがある人はいるかと思います。

第１回の記事から紹介しているシュレーディンガー方程式を考えます。

$\hat{H} \Psi = E \Psi \tag{1}$

この方程式は、ハミルトニアン $\hat{H}$ についての固有値 $E$ に対する固有方程式になっていることは以前話しました。固有値に対応して固有空間が定まり、 $\Psi$ はその空間の元であるということでした。

２つの固有値 $E = E_1, \; E_2$ に対応する解として波動関数 $\Psi_1$ と $\Psi_2$ を考えます。これらの線形結合を

$\Psi = c_1 \Psi_1 + c_2 \Psi_2$

とすると、 $\Psi$ は $\Psi_1$ という状態と $\Psi_2$ という状態が「重なり合った状態」であると考えられます。この重ね合わせもまたシュレーディンガー方程式 $(1)$ の解であり、電子の取りうる状態なのです。

・・・といいたいところなのですが、これはそもそもシュレーディンガー方程式 $(1)$ の解にはなっていません。つまり、グレーの部分が間違いです。

どういうことか、以下の枠内で説明します。

$\Psi_1, \Psi_2$ は固有値 $E_1, E_2$ に対応する解なので

$\hat{H} \Psi_1 = E_1 \Psi_1$ 　や　 $\hat{H} \Psi_2 = E_2 \Psi_2$

は成り立ちます。

一方で、これらを足し合わせた $\Psi_1 + \Psi_2$ を考えたときに

$\hat{H} (\Psi_1 + \Psi_2) = E_1 (\Psi_1 + \Psi_2)$ 　や　 $\hat{H} (\Psi_1 + \Psi_2) = E_2 (\Psi_1 + \Psi_2)$

は、一般には成り立たないわけですね。（ $E_1 = E_2$ であれば成立します。）

では、重ね合わせの原理はどのように考えたら良いのか。式 $(1)$ ではなく「時間に依存するシュレーディンガー方程式」を考えるべきだったのです。

式 $(1)$ のシュレーディンガー方程式は、正確に言えば 時間に依存しないシュレーディンガー方程式となります。

一方で次の式 $(2)$ は時間に依存するシュレーディンガー方程式です。

$\displaystyle i\hbar \frac{\partial }{\partial t} \Psi = \hat{H} \Psi \tag{2}$

波動関数 $\Psi$ としては、前回は座標 ${\bf r}$ に対する関数として考えましたが、こちらは時刻 $t$ も含めた $({\bf r}, t)$ の関数として $\Psi({\bf r}, t)$ を考えることにするのです。だから「時間に依存する」というわけです。

実は、本来のシュレーディンガー方程式は式 $(2)$ の方なんですね。

では式 $(1)$ の方は何だったかというと、これは位置 ${\bf r}$ と時刻 $t$ で変数分離して得られる方程式だったのです。

$\Psi({\bf r}, t) = \phi_1({\bf r}) \, \phi_2(t)$

とおいて式 $(2)$ に代入すると

$\displaystyle i\hbar \phi_1({\bf r}) \frac{\partial }{\partial t} \phi_2(t) = \phi_2(t) \hat{H} \phi_1({\bf r})$

となりますが、これを $\phi_1({\bf r})\, \phi_2(t)$ で割ると

$\displaystyle \frac{i\hbar}{ \phi_2(t)} \frac{\partial }{\partial t} \phi_2(t) = \frac{1}{\phi_1({\bf r})} \hat{H} \phi_1({\bf r})$

となり、左辺は位置 ${\bf r}$ だけの関数、右辺は時刻 $t$ だけの関数となりますね。よって、定数関数です。

この定数を $E$ とおくと、これは物理的にはエネルギーという意味を持ちますが、

$\displaystyle \hat{H} \phi_1({\bf r}) = E \phi_1({\bf r}) \tag{3}$

$\displaystyle i\hbar \frac{\partial }{\partial t} \phi_2(t) = E \phi_2(t) \tag{4}$

となります。式 $(3)$ が時間に依存しないシュレーディンガー方程式というわけですね。

式 $(4)$ は簡単に解けて

$\displaystyle \phi_2(t) = A e^{-i\frac{E}{\hbar}t}$

となります。

結局、係数を $\phi_1$ の方に押し付けると

$\Psi({\bf r}, t) = \phi_1({\bf r}) \, e^{-i\frac{E}{\hbar}t}$

は、時間に依存するシュレーディンガー方程式 $(2)$ の解となります。

時間に依存する因子は、絶対値をとると $1$ になるので確率分布には影響しません。だから、前回までは特に問題にならなかったわけですね。

時間に依存するシュレーディンガー方程式 $(2)$ の解 $\Psi_1({\bf r}, t)$ および $\Psi_2({\bf r}, t)$ が得られたとすると、当然

$\displaystyle i\hbar \frac{\partial }{\partial t} \Psi_1({\bf r}, t) = \hat{H} \Psi_1({\bf r}, t)$

$\displaystyle i\hbar \frac{\partial }{\partial t} \Psi_2({\bf r}, t) = \hat{H} \Psi_2({\bf r}, t)$

が成立します。

ここで、両者の重ね合わせ $\Psi({\bf r}, t) = c_1 \Psi_1({\bf r}, t) + c_2 \Psi_2({\bf r}, t)$ を考えると、微分演算子やハミルトニアンの線形性より

$\displaystyle i\hbar \frac{\partial }{\partial t} (c_1 \Psi_1({\bf r}, t) + c_2 \Psi_2({\bf r}, t) ) = \hat{H} (c_1 \Psi_1({\bf r}, t) + \Psi_2({\bf r}, t) )$

が成立するわけですね。すなわち、 $\Psi({\bf r}, t) = c_1 \Psi_1({\bf r}, t) + c_2 \Psi_2({\bf r}, t)$ も時間に依存するシュレーディンガー方程式の解になるわけです。

これこそが重ね合わせの原理です。

時間に依存しないシュレーディンガー方程式 $(2)$ は、エネルギー $E$ を陽に含まないので、異なるエネルギーを持つ波動関数の重ね合わせを考えられるというのがミソです。

さて、エネルギー $E_1, E_2$ を持つ２つの波動関数 $\Psi_1, \; \Psi_2$ の重ね合わせとして

$\Psi = c_1 \Psi_1 + c_2 \Psi_2$

を考えたときに、 $\Psi$ のエネルギーはいったいどのように表せるのでしょうか。

$\Psi_1$ のエネルギーは $E_1$ で、 $\Psi_2$ のエネルギーは $E_2$ であり、 $\Psi_1$ と $\Psi_2$ のどちらが観測されるかどうかは確率的に決まります。したがって、エネルギーの期待値を考えれば良いでしょう。

ここで、 $\Psi_1$ のエネルギーについて考えると、 $\Psi_1$ は $E_1$ に対する固有関数なので

$\displaystyle \hat{H} \Psi_1 = E_1 \Psi_1$

が成り立っています。これに左から $\Psi_1$ の複素共役 $\Psi_1^*$ を掛けて空間全体で積分すると

$\displaystyle \int \Psi_1^* \hat{H} \Psi_1\, d{\bf r} = \int \Psi_1^* E_1 \Psi_1\, d{\bf r}$

となりますが、右辺は $E_1$ を外に出すことができて

$\displaystyle \int \Psi_1^* \hat{H} \Psi_1\, d{\bf r} = E_1 \int \Psi_1^* \Psi_1\, d{\bf r}$

となります。よって

$\displaystyle E_1 = \frac{\int \Psi_1^* \hat{H} \Psi_1\, d{\bf r}} {\int \Psi_1^* \Psi_1}\, d{\bf r}$

となりエネルギーが求まるわけですね。

これのアナロジーで、 $\Psi_1$ と $\Psi_2$ の重ね合わせである $\Psi = c_1 \Psi_1 + c_2 \Psi_2$ においても、次の式 $(5)$ によって、エネルギーの期待値 $E$ が求まると考えましょう。

$\displaystyle E = \frac{\int \Psi^* \hat{H} \Psi\, d{\bf r}} {\int \Psi^* \Psi\, d{\bf r}} \tag{5}$

実際、 $\Psi_1, \, \Psi_2$ を相異なるエネルギー（固有値） $E_1, \, E_2$ に対応する波動関数だとしましょう。ここで、ハミルトニアン $\hat{H}$ はエルミート演算子と呼ばれるものになっており、この性質により相異なる固有値に対応する固有ベクトルは直交します。すなわち

$\displaystyle \int \Psi_1^* \Psi_2 \, d{\bf r} = 0$

ということです。

この事実と、波動関数 $\Psi_1, \Psi_2$ の規格化条件

$\displaystyle \int \Psi_1^* \Psi_1 \, d{\bf r} = 1$

$\displaystyle \int \Psi_2^* \Psi_2 \, d{\bf r} = 1$

を使って $E$ をより具体的に計算することができます。以下はなかなかすごい計算ですが、丁寧に1行1行読んでいけば理解できるかと思います：

$\newcommand{\HLNO}{\unicode[\sans-serif,STIXGeneral]{x306E}}\displaystyle \begin{align} E &= \frac{\int (c_1^* \Psi_1^* + c_2^* \Psi_2^*) \hat{H} (c_1 \Psi_1 + c_2 \Psi_2)\, d{\bf r}} {\int (c_1^* \Psi_1^* + c_2^* \Psi_2^*)(c_1 \Psi_1 + c_2 \Psi_2)\, d{\bf r}} \\ &= \frac{\int (c_1^* \Psi_1^* + c_2^* \Psi_2^*) (c_1 \hat{H}\Psi_1 + c_2 \hat{H}\Psi_2)\, d{\bf r}} {\int (c_1^* \Psi_1^* + c_2^* \Psi_2^*)(c_1 \Psi_1 + c_2 \Psi_2)\, d{\bf r}} & (\,\because \; \hat{H}\,\text{を適用}\,) \\ &= \frac{\int (c_1^* \Psi_1^* + c_2^* \Psi_2^*) (c_1 E_1\Psi_1 + c_2 E_2\Psi_2)\, d{\bf r}} {\int (c_1^* \Psi_1^* + c_2^* \Psi_2^*)(c_1 \Psi_1 + c_2 \Psi_2)\, d{\bf r}} & (\,\because \; \hat{H}\Psi_i = E_i \Psi_i \; \text{より}\,) \\ &= \frac{\int (|c_1|^2 E_1 \Psi_1^*\Psi_1 + |c_2|^2 E_2 \Psi_2^*\Psi_2 + c_1^* c_2 E_2 \Psi_1^* \Psi_2 + c_2^* c_1 E_1 \Psi_2^* \Psi_1) \, d{\bf r}} {\int (|c_1|^2 \Psi_1^* \Psi_1 + |c_2|^2 \Psi_2^* \Psi_2 + c_1^* c_2 \Psi_1^* \Psi_2 + c_2^* c_1 \Psi_2^* \Psi_1)\, d{\bf r}} & (\,\because \; \text{式を展開}\,) \\ &= \frac{ |c_1|^2 E_1\, \int \Psi_1^*\Psi_1\, d{\bf r} + |c_2|^2 E_2\, \int \Psi_2^*\Psi_2\, d{\bf r} + c_1^* c_2 E_2\,\int \Psi_1^* \Psi_2\, d{\bf r} + c_2^* c_1 E_1\, \int \Psi_2^* \Psi_1 \, d{\bf r}} {|c_1|^2 \, \int \Psi_1^* \Psi_1\, d{\bf r} + |c_2|^2 \,\int \Psi_2^* \Psi_2\, d{\bf r} + c_1^* c_2 \, \int \Psi_1^* \Psi_2\, d{\bf r} + c_2^* c_1 \, \int \Psi_2^* \Psi_1\, d{\bf r}} & (\,\because \; \text{定数を積分}\HLNO\text{外に}\,) \\ &= \frac{ |c_1|^2 E_1 + |c_2|^2 E_2} {|c_1|^2 + |c_2|^2} & (\,\because \; \text{積分の直交性と規格化}\,) \\ &= \frac{|c_1|^2}{|c_1|^2 + |c_2|^2} E_1 + \frac{|c_2|^2}{|c_1|^2 + |c_2|^2} E_2 \\ \end{align}$

得られた式を観察しましょう。

係数 $\frac{|c_1|^2}{|c_1|^2 + |c_2|^2}$ の部分は重ね合わせの係数の絶対値の二乗をとって、その比率を考えたものとなっています。これは、 $\Psi_1$ が観測される確率だと解釈できますね。同様に、 $\frac{|c_2|^2}{|c_1|^2 + |c_2|^2}$ も $\Psi_2$ が観測される確率だと解釈できます。

すると、 $E$ はまさにエネルギーの期待値だと思うことができますね！

これにて準備が整いました。

3. 分子軌道法の計算

それでは、分子軌道法の具体的な計算を実行したいと思います。

ここから計算したいのは、水素分子のように原子核が２つある分子における、電子１つの波動関数です。

ハミルトニアン

$\displaystyle \hat{H}_{\text{分子}} = -\frac{\hbar}{2m_e}\left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}\right) - \frac{A}{r_1} - \frac{B}{r_2} + \frac{C}{R} \tag{6}$

を考えます。これはポテンシャルとして、各原子核から電子へのクーロン力と、原子核同士のクーロン力をそれぞれ考えたものとなっています。このハミルトニアンを使って、シュレーディンガー方程式

$\displaystyle i\hbar \frac{\partial}{\partial t} \Psi = \hat{H}_{\text{分子}} \Psi$

を考えるわけですね。

これは簡単に解けないので、知っている状況に帰着しようということになります。

もし仮に原子核2がなかったとすると、電子には原子核1の寄与しかありませんので、水素原子のシュレーディンガー方程式を解くことができます。

原子核 $1$ に対する水素原子の軌道に対応する波動関数を $\chi_1$ としましょう。水素原子のときにやったように、軌道は複数取りうるわけなので、どれか１つを選びます。
（あとで、具体的な軌道について議論することになります。）

また同様に原子核2についてのシュレーディンガー方程式の解をやはり $\chi_2$ とおくことにしましょう。

ここで、原子核の距離が十分大きいとすると、つまり $R$ が十分大きいと考えると $\chi_1, \; \chi_2$ は分子におけるシュレーディンガー方程式を近似的に満たすと考えられます。 $R$ を小さくしていくと、相互の原子核の影響を受けることになり、解の様子は $\chi_1, \chi_2$ から離れていきます。

分子軌道法では、分子のシュレーディンガー方程式の真の解 $\Psi_{\text{真}}$ は、 $\chi_1$ と $\chi_2$ の重ね合わせによって近似できると考えます。つまり、線形結合

$\Psi = c_1 \chi_1 + c_2 \chi_2$

を考えます。これが $\Psi_{\text{真}}$ の近似になるような $(c_1, c_2)$ の組を探せばよい、というのがアイデアです！

前節で一生懸命説明したように、波動関数の重ね合わせ $\Psi = c_1 \chi_1 + c_2 \chi_2$ についてのエネルギーの期待値は

$\displaystyle E = \frac{\int \Psi^* \hat{H}_{\text{分子}} \Psi\, d{\bf r}} {\int \Psi^* \Psi\, d{\bf r}} \tag{5再掲}$

で表せるのでした。これは $(c_1, c_2)$ についての 2変数関数 $E = E(c_1, c_2)$ となっています。

変分原理によれば、分子の取りうる真の軌道のエネルギー $E_{\text{真}}$ は、必ず

$E_{\text{真}} \leq E(c_1, c_2)$

を満たすはずなので、 $E(c_1, c_2)$ の極小値を決定すればそれが真のエネルギーの近似値と思えるという寸法です！

さぁ、面白くなってきました！！

具体的な計算を実行するために、エネルギーの式 $(5)$ に $\Psi = c_1 \chi_1 + c_2 \chi_2$ を代入して計算しましょう。

$\displaystyle \begin{align} E &= \frac{\int (c_1 \chi_1^* + c_2 \chi_2^*) \hat{H}_{\text{分子}} (c_1 \chi_1 + c_2 \chi_2)\, d{\bf r}} {\int (c_1 \chi_1^* + c_2 \chi_2^*)(c_1 \chi_1 + c_2 \chi_2)\, d{\bf r}} \\ &= \frac{\int (c_1 \chi_1^* + c_2 \chi_2^*) (c_1 \hat{H}_{\text{分子}}\chi_1 + c_2 \hat{H}_{\text{分子}}\chi_2)\, d{\bf r}} {\int (c_1 \chi_1^* + c_2 \chi_2^*)(c_1 \chi_1 + c_2 \chi_2)\, d{\bf r}} & (\,\because \; \hat{H}_{\text{分子}}\,\text{を適用}\,) \\ &= \frac{\int (c_1^2 \chi_1^*\hat{H}_{\text{分子}}\chi_1 + c_2^2 \chi_2^*\hat{H}_{\text{分子}}\chi_2 + c_1 c_2 \chi_1^* \hat{H}_{\text{分子}}\chi_2 + c_1 c_2 \chi_2^* \hat{H}_{\text{分子}}\chi_1) \, d{\bf r}} {\int (c_1^2 \chi_1^* \chi_1 + c_2^2 \chi_2^* \chi_2 + c_1 c_2 \chi_1^* \chi_2 + c_1 c_2 \chi_2^* \chi_1)\, d{\bf r}} & (\,\because \; \text{式を展開する}\,) \\ &= \frac{ c_1^2 \, \int \chi_1^*\hat{H}_{\text{分子}}\chi_1\, d{\bf r} + c_2^2 \, \int \chi_2^*\hat{H}_{\text{分子}}\chi_2\, d{\bf r} + c_1 c_2 \,\int \chi_1^* \hat{H}_{\text{分子}}\chi_2\, d{\bf r} + c_1 c_2 \, \int \chi_2^* \hat{H}_{\text{分子}}\chi_1 \, d{\bf r}} {c_1^2 \, \int \chi_1^* \chi_1\, d{\bf r} + c_2^2 \,\int \chi_2^* \chi_2\, d{\bf r} + c_1 c_2 \, \int \chi_1^* \chi_2\, d{\bf r} + c_1 c_2 \, \int \chi_2^* \chi_1\, d{\bf r}} & (\,\because \; \text{係数を積分}\HLNO\text{外に}\,) \end{align}$

分母と分子に積分がたくさん出てきているので、それぞれ次のように置きましょう。

$\displaystyle H_{11} = \int \chi_1^*\hat{H}_{\text{分子}}\chi_1\, d{\bf r}$
$\displaystyle H_{12} = \int \chi_1^*\hat{H}_{\text{分子}}\chi_2\, d{\bf r}$
$\displaystyle H_{21} = \int \chi_2^*\hat{H}_{\text{分子}}\chi_1\, d{\bf r}$
$\displaystyle H_{22} = \int \chi_2^*\hat{H}_{\text{分子}}\chi_2\, d{\bf r}$
$\displaystyle S_{11} = \int \chi_1^* \chi_1\, d{\bf r} = 1$ （ $\because$ $\chi_1$ は規格化されている）
$\displaystyle S_{12} = \int \chi_1^* \chi_2\, d{\bf r}$
$\displaystyle S_{21} = \int \chi_2^* \chi_1\, d{\bf r}$
$\displaystyle S_{22} = \int \chi_2^* \chi_2\, d{\bf r} = 1$ （ $\because$ $\chi_2$ は規格化されている）

これらの積分は名前がついていまして、何度も出てくるので説明しておきましょう。

★ $i = j$ のとき：

$H_{ii}$ ：クーロン積分
$S_{ii}$ ：規格化条件により、積分値は $1$ となる

★ $i \neq j$ のとき：

$H_{ij}$ ：共鳴積分
$S_{ij}$ ：重なり積分

これらを使うとエネルギーの期待値は

$\displaystyle E = \frac{ c_1^2 \, H_{11} + c_2^2 \, H_{22} + c_1 c_2 (H_{12} + H_{21})} { c_1^2 + c_2^2 + c_1 c_2 (S_{12} + S_{21}) }$

と表すことができます。

さらに、以下の枠内の理屈により $H_{12} = H_{21}$ や $S_{12} = S_{21}$ という性質もあるので、もう少し簡潔に

$\displaystyle E = \frac{ c_1^2 \, H_{11} + c_2^2 \, H_{22} + 2c_1 c_2 H_{12} } { c_1^2 + c_2^2 + 2c_1 c_2 S_{12} } \tag{7}$

とできます。

まず、ハミルトニアンはエルミート演算子という性質を持った演算子であり、このことから

$H_{12}^* = H_{21}$

の性質を持ちます。 $H_{12}$ も $H_{21}$ 実数なので、結果的に

$H_{12} = H_{21} \tag{8}$

であることが帰結されます。また

$S_{12} = S_{21} \tag{9}$

も同様です。

さて、これで波動関数 $\Psi = c_1 \chi_1 + c_2 \chi_2$ に対応するエネルギー（の期待値）を求めることができました。これは $c_1, c_2$ の関数だと思うことができます。

変数 $(c_1, c_2)$ に対する $E(c_1, c_2)$ の最小化問題を解いたときに、得られる解 $(c_1, c_2)$ に対応する波動関数が、分子の真の軌道の近似になっているというわけですね。

こういう方法のことを変分法というそうですが、要するに最小化問題だと思って差し支えないと思います。

最小化問題を解くための簡便な方法は、偏微分することです。すなわち、エネルギーを $c_1$ と $c_2$ でそれぞれ偏微分して

$\begin{cases} \displaystyle \frac{\partial E}{\partial c_1} = 0 \\ \displaystyle \frac{\partial E}{\partial c_2} = 0 \end{cases} \tag{10}$

を満たすような $c_1, c_2$ を求めればよいのです。

では、実際に偏微分を計算したいと思いますが、少し工夫します。エネルギーの式 $(7)$ の分母を両辺にかけて

$\displaystyle (c_1^2 + c_2^2 + 2c_1 c_2 S_{12})\, E = c_1^2 \, H_{11} + c_2^2 \, H_{22} + 2c_1 c_2 H_{12} \tag{11}$

とします。この全体を $c_1$ で偏微分すると

$\displaystyle (2c_1 + 2c_2 S_{12}) E + (c_1^2 + c_2^2 + 2c_1 c_2 S_{12})\frac{\partial E}{\partial c_1} = 2c_1 H_{11} + 2c_2 H_{12}$

となりますが、 $\frac{\partial E}{\partial c_1} =$ の形にすると

$\displaystyle \frac{\partial E}{\partial c_1} = \frac{2\left\{ (H_{11} - E)c_1 + (H_{12} - E S_{12})c_2 \right\}}{c_1^2 + c_2^2 + 2c_1 c_2 S_{12}} \tag{12}$

となります。

$\frac{\partial E}{\partial c_1} = 0$ であることは、式 $(12)$ の分子が $0$ になることと同値ですから

$(H_{11} - E)c_1 + (H_{12} - E S_{12})c_2 = 0$

が得られます。

同様に、 $c_2$ で偏微分すると

$(H_{12} - E S_{12})c_1 + (H_{22} - E )c_2 = 0$

も得られます。

すなわち、式 $(10)$ と同値な条件として $c_1, c_2$ を変数とする連立方程式

$\begin{cases} (H_{11} - E)c_1 + (H_{12} - E S_{12})c_2 = 0 \\ (H_{12} - E S_{12})c_1 + (H_{22} - E )c_2 = 0 \end{cases} \tag{13}$

が得られるわけです。

行列で表すと

$\displaystyle \begin{pmatrix} H_{11} - E & H_{12} - E S_{12} \\ H_{12} - E S_{12} & H_{22} - E \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \tag{14}$

ということですね。

さて、式 $(14)$ が $(c_1, c_2) = (0, 0)$ ではない解を持つための必要十分条件は

$\displaystyle \operatorname{det}\begin{pmatrix} H_{11} - E & H_{12} - E S_{12} \\ H_{12} - E S_{12} & H_{22} - E \end{pmatrix} = 0 \tag{15}$

であることです。

これは完全に線形代数の問題なので、この同値性が飲み込める方はそのまま進んでいただいて、気になる方は以下の証明を読んでください。

ここで使っているのは、線形代数の以下の定理です。この定理は上記で使っている事実の、ちょうど対偶になっています。

定理

$n$ 次正方行列 $A$ に対して、以下は同値：

$\operatorname{det} A \neq 0 \;\; \Longleftrightarrow \;\; A{\bf x} = {\bf 0} \, \text{が}\, {\bf x} \neq {\bf 0} \, \text{なる解を持たない}$

（証明）

$\Longrightarrow$ を示す：

仮定より $A$ は正則だから、逆行列 $A^{-1}$ が存在する。これを $A{\bf x} = {\bf 0}$ に左から掛けると ${\bf x} = {\bf 0}$ となり、 ${\bf x} \neq {\bf 0}$ なる解を持たないことが言えた。

$\Longleftarrow$ を示す：

仮定の条件を言い換えると「 $A{\bf x} = {\bf 0} \; \Longrightarrow \; {\bf x} = {\bf 0}$ 」ということである。行列 $A$ を $A = ({\bf a}_1, \ldots, {\bf a}_n)$ とベクトル表記し、 ${\bf x} = (x_1, \ldots, x_n)$ とすると、仮定の条件は

$x_1 {\bf a}_1 + \cdots + x_n {\bf a}_n = {\bf 0} \; \Longrightarrow \; {\bf x} = {\bf 0}$

ということである。これは、ベクトル ${\bf a}_1, \ldots, {\bf a}_n$ が一次独立であることを意味する。

ここで、一次独立なベクトルを並べた行列 $A$ は正則であることを示す。（これが言えれば、 $\operatorname{det}A = 0$ がしたがう。）

まず行列 $A$ に対する連立方程式 $A{\bf u} = {\bf b}$ を考える。これを行基本変形によって簡約化すると $A^r {\bf u} = {\bf b}^r$ が得られる。また、一次独立性より $A^r = I$ は単位行列より、 ${\bf u} = {\bf b}^r$ として解が得られる。したがって、任意の ${\bf b}$ について一意に解が定まる。

ここで、基底ベクトル ${\bf e}_i = (0, \ldots, 1, \ldots, 0)$ （ $i$ 番目が $1$ ）に対して ${\bf b} = {\bf e}_i$ の解を ${\bf u} = {\bf u}_i$ とおく。すると行列

$U = ({\bf u}_1, \ldots, {\bf u}_n)$

が定義できる。これを $A$ に左からかけると

$\begin{align} AU &= A({\bf u}_1, \ldots, {\bf u}_n) \\ &= (A{\bf u}_1, \ldots, A{\bf u}_n) \\ &= ({\bf e}_1, \ldots, {\bf e}_n) \\ &= I \end{align}$

であり、単位行列になる。右からかけても同じく単位行列になるので、 $U$ は $A$ の逆行列である。したがって $A$ は正則行列であり、 $\operatorname{det} A \neq 0$ が言えた。

（証明終わり）

あらためて、結論をいうと式 $(10)$

$\begin{cases} \displaystyle \frac{\partial E}{\partial c_1} = 0 \\ \displaystyle \frac{\partial E}{\partial c_2} = 0 \end{cases} \tag{10再掲}$

が $(c_1, c_2) \neq (0, 0)$ なる解を持つことと、式 $(15)$

$\displaystyle \operatorname{det}\begin{pmatrix} H_{11} - E & H_{12} - E S_{12} \\ H_{12} - E S_{12} & H_{22} - E \end{pmatrix} = 0 \tag{15再掲}$

は同値な条件ということになります。

式 $(15)$ は $E$ についての方程式となっていますが、この式を永年方程式というそうです。ちょっと、名前がかっこいいですね！

永年方程式は

（行列の $(i, j)$ -成分が $H_{ij} - E S_{ij}$ であるような対称的な行列の行列式） $= 0$

という形の方程式です。

現在は２つの原子核について考えているので $2\times 2$ 行列が出ていますが、一般の $n$ 個の原子核の場合においては $n\times n$ の同じ形の行列式が現れます。

さて、永年方程式 $(15)$ を解くために左辺の行列式を展開すると

$(H_{11} - E)(H_{22} - E) - (H_{12} - E S_{12})^2 = 0 \tag{16}$

となります。

これは $E$ についての２次方程式になっていますので、２次方程式を解けば高々２つのエネルギー $E = E_1, \, E_2$ が得られます。これらに対応する解が分子軌道というわけですね。

式 $(16)$ のままでは計算しづらいですが、もしも

$H_{11} = H_{22} \tag{17}$

であったと仮定しましょう。これは物理的には「原子軌道 $\chi_1, \; \chi_2$ のエネルギーがそれぞれ等しい」ということを意味します。

実際、 $H_{11}$ の定義を思い出すと

$\displaystyle H_{11} = \int \chi_1^* \hat{H}_{\text{分子}} \chi_1 \, d{\bf r}$

でしたが、これは $\chi_1$ のエネルギーを求める式

$\displaystyle H_{11} = \frac{\int \chi_1^* \hat{H}_{\text{分子}} \chi_1 \, d{\bf r}}{\int \chi_1^* \chi_1 \, d{\bf r}}$

において、分母を $1$ としたものになっています。 $\chi_1$ は正規化されているので、分母は $1$ となります。

よって、 $H_{11}$ は原子軌道 $\chi_1$ のエネルギーそのものだと思うことができますね。

$H_{22}$ も同様に原子軌道 $\chi_2$ のエネルギーだと思うことできます。 $H_{11} = H_{22}$ はこれらのエネルギーが等しいことを表します。

つまり、① 同じ元素による２原子核を考えていて、また、② $\chi_1$ と $\chi_2$ として同じ原子軌道を選択する、という仮定を置いている状況では、 $H_{11} = H_{22}$ は自然に成立するわけですね。

このような条件では $H_{11} = H_{22}$ より、式 $(16)$ は

$(H_{11} - E)^2 - (H_{12} - E S_{12})^2 = 0 \tag{18}$

と簡単化されます。

これは「二乗マイナス二乗＝０」の形になっていますので

$H_{11} - E + H_{12} - E S_{12} = 0$ 　または　 $H_{11} - E - H_{12} + E S_{12} = 0$

ということになります。つまり

$\displaystyle E = \frac{H_{11} + H_{12}}{1 + S_{12}}$ 　または　 $\displaystyle E = \frac{H_{11} - H_{12}}{1 - S_{12}}$

が解となります。

２つの解にそれぞれ $E_1 = \frac{H_{11} + H_{12}}{1 + S_{12}}$ 、 $E_2 = \frac{H_{11} - H_{12}}{1 - S_{12}}$ と名前をつけておきましょう。

ここで、 $H_{11} = H_{22} < 0$ （エネルギーが負であるから）であり、 $0 < S_{12} < 1$ （内積なので）であることを考慮に入れると

$E_1 < E_2$

が成り立ちます。

つまり、取りうる分子軌道は２つあって、それぞれ $E_1$ と $E_2$ というエネルギーを持ち、 $E_1$ の方がよりエネルギーが低い、ということになります。この $E_1$ に対応する軌道を結合性軌道といい、もうひとつの $E_2$ に対応する軌道を反結合性軌道といいます。

分子軌道のエネルギーを図に表すと次のようになります：

原子核１の原子軌道 $\chi_1$ と原子核２の原子軌道 $\chi_2$ が重なり合って、２つの分子軌道（結合性軌道と反結合性軌道）が生まれます。

結合性軌道の方は、エネルギーが元の原子軌道のものよりも低く、電子はより低いエネルギーの軌道を求めて結合性軌道に入ります。こうして安定な結合が形成されるというわけです。

おお、これが冒頭で言っていた「結合性軌道」と「反結合性軌道」だったのですね！　大学一年のときに分からなかった謎概念がようやく明らかになりました！

たしかに、エネルギーの極小値を求めると、必然的に２つの軌道が出てきますね！

結合性軌道と反結合性軌道の各エネルギー $E_1, E_2$ から、実際に軌道の係数 $c_1, c_2$ を求めてみましょう。

永年方程式の１個前の連立方程式

$\displaystyle \begin{pmatrix} H_{11} - E & H_{12} - E S_{12} \\ H_{12} - E S_{12} & H_{22} - E \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \tag{14再掲}$

に、 $E = E_1$ や $E = E_2$ を代入します。

実際、 $E = E_1 \; \left(\,= \frac{H_{11} + H_{12}}{1 + S_{12}} \,\right)$ を代入すると、1行目の式は

$\displaystyle \left(H_{11} - \left(\frac{H_{11} + H_{12}}{1 + S_{12}}\right)\right)c_1 + \left(H_{12} - \left(\frac{H_{11} + H_{12}}{1 + S_{12}}\right) S_{12}\right) c_2 = 0$

であり、計算すると

$c_1 = c_2 \tag{19}$

が得られます。
（式 $(14)$ の2行目からも同じ条件が得られます。）

また、 $\Psi = c_1 \chi_1 + c_2 \chi_2$ が正規化されているためには

$\displaystyle \int \Psi^* \Psi \, d{\bf r} = c_1^2 + c_2^2 + 2c_1 c_2 S_{12} = 1 \tag{20}$

である必要があります。

式 $(20)$ に式 $(19)$ を代入して、 $c_2$ を消去すると $c_1^2(2 + 2S_{12}) = 1$ となり

$\displaystyle c_1 = \frac{1}{\sqrt{2 + 2S_{12}}}$

が得られます。 $c_1 = c_2$ より

$\displaystyle \Psi = \frac{1}{\sqrt{2 + 2S_{12}}}\chi_1 + \frac{1}{\sqrt{2 + 2S_{12}}}\chi_2 \tag{21}$

が得られるということになります。これが結合性軌道の波動関数です。

同様に、 $E = E_2 \; \left(\,= \frac{H_{11} - H_{12}}{1 - S_{12}} \,\right)$ からは、 $c_1 = -c_2$ が得られ、反結合性軌道の波動関数

$\displaystyle \Psi^* = \frac{1}{\sqrt{2 - 2S_{12}}}\chi_1 - \frac{1}{\sqrt{2 - 2S_{12}}}\chi_2 \tag{22}$

が得られます。反結合性軌道の方は、波動関数の記号にアスタリスクをつけて表すのが一般的です。

さて、結合性軌道 $\Psi$ の方は、原子軌道の符号が一致するように重なるので、強め合うように結合が形勢されるというイメージです。一方で、反結合性軌道 $\Psi^*$ の方は、原子軌道の符号がちょうど反転していますので、逆に打ち消しあうようになるわけですね。

この辺のイメージは、実際の分子を想定した計算を実行することで確認してみましょう。

4. 水素分子のσ結合

より具体的な軌道を想定して分子軌道を計算し、実際に可視化を行うところまでやってみましょう。まず、最初に考えるのは水素分子です。

水素原子 ${}_1\text{H}$ は、1s軌道に1つの電子が入っている非常にシンプルな構造をしています。水素分子 $\text{H}_2$ は、各原子から $1$ つずつ電子が供給されて、2つの電子による共有結合が生じています。

結合の方向が結合軸（原子核同士を結んだ直線方向）に沿った結合になっており、このような結合のことを σ結合といいます。

このときの分子軌道を記述するためには、２つの水素原子の1s軌道に対応する波動関数 $\chi_1 = {\rm s}^{(1)}$ と $\chi_2 = {\rm s}^{(2)}$ を考えて、この線形結合

$c_1 \chi_1 + c_2 \chi_2$

を考えれば良いでしょう。

今、２原子間の結合を考えており、それぞれ同じ軌道（1s軌道）同士を考えています。したがって、対称性から $H_{11} = H_{22}$ と仮定してよいことになります。

このとき、取りうるエネルギーは

$\displaystyle E_1 = \frac{H_{11} + H_{12}}{1 + S_{12}}$ 　または　 $\displaystyle E_2 = \frac{H_{11} - H_{12}}{1 - S_{12}}$

であり、 $E_1 < E_2$ の関係にあります。

$E_1, E_2$ に対応する分子軌道は２つあり、それぞれ結合性軌道と反結合性軌道というのでした。今回は、σ結合を考えているので、結合性軌道を $\sigma_{1{\rm s}}$ 、反結合性軌道を $\sigma_{1{\rm s}}^*$ と表すことにします。

具体的には、この軌道は式 $(19), (20)$ により次のように表されるのでした：

$\displaystyle \sigma_{1{\rm s}} = \frac{1}{\sqrt{2 + 2S_{12}}}\chi_1 + \frac{1}{\sqrt{2 + 2S_{12}}}\chi_2$

$\displaystyle \sigma_{1{\rm s}}^* = \frac{1}{\sqrt{2 - 2S_{12}}}\chi_1 - \frac{1}{\sqrt{2 - 2S_{12}}}\chi_2$

あと必要な情報としては $H_{11}, \; H_{12}$ と $S_{12}$ です。

$H_{11}$ はクーロン積分、 $H_{12}$ は共鳴積分という名前がついているのでした。これらはエネルギーを計算するのに使う量なので、分子軌道の波動関数を計算するのには取り急ぎ必要ありません。

$S_{12}$ は重なり積分というのでした。これは、おおよそ「波動関数の重なり度合い」を表す量で、２つの原子軌道の波動関数がどれぐらい重なっているかを表します。（あまり重なっていなければ値は小さく、重なっていれば値は大きいということですね。）

各分子軌道の波動関数を計算するには $S_{12}$ が必要ですが、これには積分

$\displaystyle S_{12} = \int \chi_1^* \chi_2 \, d{\bf r}$

を計算すれば良いことになります。

ただ、この計算は面倒なので、ここでは教科書から数値をお借りします。参考文献の教科書によると、水素原子の1s軌道同士の重なり積分の値は $0.6$ 程度になるそうです。これはかなり大きな値だそうで、第２周期以上の元素においては1s軌道の重なりは0.2〜0.3程度になるそうです。水素原子はだいぶ例外的なのですね。
（たしかに、 $Z = 1$ なので半径は大きいため、重なりは大きくなりそうですね。）

そこで、 $S_{12} = 0.6$ としましょう。すると

$\displaystyle \sigma_{1{\rm s}} = \frac{1}{\sqrt{3.2}}\chi_1 + \frac{1}{\sqrt{3.2}}\chi_2 \tag{23}$

$\displaystyle \sigma_{1{\rm s}}^* = \frac{1}{\sqrt{0.8}}\chi_1 - \frac{1}{\sqrt{0.8}}\chi_2 \tag{24}$

となります。

実際、上記の結合性軌道と反結合性軌道の波動関数について、前回やったように「確率が90%になる領域」について可視化すると次のようになります。（プログラムはあとで紹介します。）

図：水素分子（1s軌道）の結合性軌道 $\sigma_{1{\rm s}}$

図：水素分子（1s軌道）の反結合性軌道 $\sigma_{1{\rm s}}^*$

前回までの記事同様、赤色は波動関数が正の部分、青色が負の部分をそれぞれ表しています。

図から、結合性軌道の方は、２つの1s軌道がくっついて結合を形成しているように見えます。一方、反結合性軌道の方は、位相が逆になっていることにより打ち消され、原子核と原子核の間の部分の確率密度が小さくなっているように見えます。

実際、反結合性軌道を $y$ 軸方向から見ると、こうなっています。

水素原子の1s軌道が近づいていくと、その重なりにより２つの安定な分子軌道が生じます。電子は各軌道にやはり２つずつ入ることができ、水素分子の場合は２つの電子が、より安定な結合性軌道に入ることになります。

結合性軌道では分子間の結合を形成します。何かの拍子に電子が不安定な反結合性軌道に入ると、上記のように電子の分布は結合部で打ち消しあい、結果的に結合が切れることになります。

なお、上記の可視化においては、２つの水素原子の位置を $(x, y, z) = (1, 0, 0)$ と $(x, y, z) = (-1, 0, 0)$ においた想定で計算しています。（つまり、原始間距離はボーア半径の2倍。）

このときの原始軌道 $\chi_1, \chi_2$ を計算するには、水素原子が原点に置かれているときの1s軌道の波動関数は量子数 $(n, l, m) = (1, 0, 0)$ の波動関数なので $\Psi_{(1, 0, 0)}(x, y, z)$ として

$\chi_1(x, y, z) = \Psi_{1,0,0}(x - 1, y, z)$

$\chi_2(x, y, z) = \Psi_{1,0,0}(x + 1, y, z)$

のように計算すればよいですね。

5. エチレン分子のπ結合

前節はσ結合の例を紹介しました。今度は π結合について考えたいと思います。π結合とは、σ結合と異なり、結合方向と直角の方向を向いた原子軌道同士の結合のことです。

具体的には、エチレン分子 $\text{CH}_2$ について考えてみましょう。

炭素原子には全部で6つの電子があり、内訳はK殻（1s軌道）に2個、L殻（2s軌道と2p軌道）に4個となっています。つまり、価電子（最外殻電子）は４つです。

エチレンの炭素同士は、炭素原子核1からは価電子のうち2個、もう一方の炭素原子核2からも価電子のうち2個、それぞれ2個ずつの電子を供給しあって結合しています。このような結合を二重結合というのでした。

この二重結合の内訳ですが、1本目が原子核の結合方向を向いた原子同士による結合（つまりσ結合）と、結合方向と垂直な方向を向いた原子軌道による結合（つまりπ結合）となっています。

π結合のイメージ

ここではこのπ結合にフォーカスしましょう。p軌道の方向を $z$ 軸とすると、π結合を構成する原子軌道は原子核１の ${\rm p}_z$ 軌道と、原子核２の ${\rm p}_z$ 軌道ということになります。

${\rm p}_z$ 軌道の形状（赤い部分は波動関数が正の部分、青い部分は波動関数が負の部分）

$\chi_1 = {\rm p}_z^{(1)}, \;\; \chi_2 = {\rm p}_z^{(2)}$ として、水素分子のときと同様に計算を行うと、結合性軌道 $\pi_{2{\rm p}}$ ・反結合性軌道 $\pi_{2{\rm p}}^*$ の波動関数はそれぞれ

$\displaystyle \pi_{2{\rm p}} = \frac{1}{\sqrt{2+2S_{12}}}\chi_1 + \frac{1}{\sqrt{2+2S_{12}}}\chi_2$

$\displaystyle \pi_{2{\rm p}}^* = \frac{1}{\sqrt{2-2S_{12}}}\chi_1 - \frac{1}{\sqrt{2-2S_{12}}}\chi_2$

となります。π結合なので、波動関数の記号は $\pi$ の記号を用います。

対応するエネルギー $E_1, E_2$ は

$\displaystyle E_1 = \frac{H_{11} + H_{12}}{1 + S_{12}}$
$\displaystyle E_2 = \frac{H_{11} - H_{12}}{1 - S_{12}}$

となるのでした。

ここでやはり、 $H_{11}, \; H_{12}, \; S_{12}$ の値が必要になります。

$H_{11}, \; H_{12}$ は計算するか、実験的に調べるしかありません。ひとまず定数として

$H_{11} = \alpha, \;\; H_{12} = \beta$

と置いておきましょう。（ $\alpha < 0, \;\; \beta < 0$ となります。）

π結合における重なり積分は、あとで述べるヒュッケル法という考え方によって $S_{12} = 0$ としてよいことが分かります。

したがって、結合性軌道 $\pi_{2{\rm p}}$ と反結合性軌道 $\pi_{2{\rm p}}^*$ は

$\displaystyle \pi_{2{\rm p}} = \frac{1}{\sqrt{2}}\chi_1 + \frac{1}{\sqrt{2}}\chi_2 \tag{26}$

$\displaystyle \pi_{2{\rm p}}^* = \frac{1}{\sqrt{2}}\chi_1 - \frac{1}{\sqrt{2}}\chi_2 \tag{27}$

となり、対応するエネルギーは

$\displaystyle E_1 = \alpha + \beta \tag{28}$

$\displaystyle E_2 = \alpha - \beta \tag{29}$

と簡単に求められます。

それでは、これらの分子軌道の様子を可視化してみましょう。

原子核の距離は、エチレンの炭素間の結合距離は 133.9 pm だそうなので、距離が $d = 133.9/52.9$ ボーア半径になるように炭素原子核を配置します。

すなわち

$\chi_1(x, y, z) = {\rm p}_z^{(1)}(x, y, z) = \Psi_{2,1,0}(x - d/2, \, y, \, z)$

$\chi_2(x, y, z) = {\rm p}_z^{(2)}(x, y, z) = \Psi_{2,1,0}(x + d/2, \, y, \, z)$

ということです。

また、炭素電子の2p軌道なので、有効核電荷は

$Z_{\text{eff}} = Z - \sigma = 6 - (0.35 \times 3 + 0.85\times 2) \tag{30}$

となります。

# 記事の一番最後に完全版のプログラムを掲載します。
# エチレン分子のπ_2p軌道
def ethylene_molecule(x, y, z, id):
    Z = 6        # 炭素原子を想定しているので Z = 6
    Z_eff = Z - (0.35*3 + 0.85*2)  # 有効核電荷（スレータ―則により計算）
    n = 2
    l = 1
    m = 0
    distance = 133.9/52.9  # 水素原子のボーア半径 52.9 pm, エチレン分子の原子間距離 133.9 pm
    if id == 0:
        # π_2p軌道（結合性軌道）
        return (f(x-distance/2, y, z, n, l, m, Z_eff)/np.sqrt(2.0)) + (f(x+distance/2, y, z, n, l, m, Z_eff)/np.sqrt(2.0)) 
    else:
        # π_2p*軌道（反結合性軌道）
        return (f(x-distance/2, y, z, n, l, m, Z_eff)/np.sqrt(2.0)) - (f(x+distance/2, y, z, n, l, m, Z_eff)/np.sqrt(2.0))

先ほどと同様のプログラム（最後にまとめて紹介します）で可視化したものが次の図になります。

左側が結合性軌道 $\pi_{1{\rm p}}$ で、右側が反結合性軌道 $\pi_{1{\rm p}}^*$ です。（波動関数が正の部分を赤色、負の部分を青色で表示しています。）

結合性軌道の方は、位相が揃っているので、分子に対してz軸方向上側と下側にそれぞれ電子雲ができていて、結合を形成しているように見えますね。

一方、反結合性軌道の方は、位相が反転している関係で分子間には「ちょうど電子が存在しない部分」ができているように見えます。

エネルギーについて見てみると、以下のようになります。

水素分子のときと同じように、原子軌道（ ${\rm p}_z$ 軌道）のエネルギーに対して、上下に２つの軌道が生まれます。エネルギーの低い方が結合性軌道で、高い方が反結合性軌道です。

また、 $S_{12} = 0$ と近似したことによって、エネルギーが $E_{1} = \alpha + \beta, \; E_{2} = \alpha - \beta$ と対称的になったことが重要なポイントです。（ $\alpha$ を中心に、 $\pm |\beta|$ だけ上下しています。）

エチレンの場合は、炭素原子の ${\rm p}_z$ 軌道電子が各1つずつ分子軌道に入るわけですが、2個なのでどちらも結合性軌道に入ります。どちらも結合性軌道に入るということで、安定した結合が形成されるようです。

電子２つはどちらも結合性軌道 $\pi_{2{\rm p}}$ に入る

6. ヒュッケル法とブタジエンのπ電子共役系

前節では分子として２炭素原子核によって構成されたものを考えてきましたが、ここからは原子核の数を増やしてみたいと思います。

例としてブタジエン $\text{C}_4 \text{H}_6$ を考えたいと思うのですが、炭素だけに着目すると構造式は次のようになります。

ちなみに、今回考えるブタジエンは、正確に言えば1-3ブタジエンと呼ぶべきものです。

他にも、二重結合の位置が異なる1-2ブタジエンという構造異性体がありますが、今回は扱いません。

4つの炭素原子を左から1,2,3,4と名前をつけることにします。

二重結合が1-2間, 3-4間にありますが、π結合の電子は必ずしも1-2間, 3-4間にのみあるわけではありません。実際は、1,2,3,4の炭素全体にまたがって存在しているのです。こういう状況を非局在化といいます。

つまり、「局」所的に「在」るわけではない、というわけですね。

（ブタジエンのように）二重結合・単結合が交互に並んだ結合において、π電子が非局在化しているような系を一般に π電子共役系といいます。

ブタジエンにおいては、分子を構成する4つの炭素原子にはそれぞれ ${\rm p}_z$ 軌道 $\chi_1, \chi_2, \chi_3, \chi_4$ があります。π電子が非局在化している分子軌道に対しては、これらの重ね合わせ（つまり線形結合）

$\displaystyle \Psi = c_1 \chi_1 + c_2 \chi_2 + c_3 \chi_3 + c_4 \chi_4 \tag{31}$

を考えることになるわけです。

あとは、この軌道に対するエネルギーの期待値

$\displaystyle E = \frac{\int \Psi^* \hat{H}_{\text{分子}} \Psi \, d{\bf r}}{\int \Psi^* \Psi \, d{\bf r}} \tag{32}$

を計算し、この極値を求めるとそれが分子軌道になるわけですね。

実際

$\displaystyle \frac{\partial E}{\partial c_1} = 0, \;\; \frac{\partial E}{\partial c_2} = 0, \;\; \frac{\partial E}{\partial c_3} = 0, \;\; \frac{\partial E}{\partial c_4} = 0 \tag{33}$

なる解を求めようとすると、 $4$ 変数の連立方程式が得られます。

この連立方程式が $(c_1, c_2, c_3, c_4) \neq (0, 0, 0, 0)$ なる解を持つことと同値な条件として、 $(\text{行列式}) = 0$ という形の条件が得られます。具体的には

$\operatorname{det} \begin{pmatrix} H_{11} - E S_{11} & H_{12} - E S_{12} & H_{13} - E S_{13} & H_{14} - E S_{14} \\ H_{21} - E S_{21} & H_{22} - E S_{22} & H_{23} - E S_{23} & H_{24} - E S_{24} \\ H_{31} - E S_{31} & H_{32} - E S_{32} & H_{33} - E S_{33} & H_{34} - E S_{34} \\ H_{41} - E S_{41} & H_{42} - E S_{42} & H_{43} - E S_{43} & H_{44} - E S_{44} \end{pmatrix} = 0 \tag{34}$

となりますが、これを永年方程式というわけですね。 $H_{ij}, \; S_{ij}$ は定数なので、これは $E$ についての４次方程式となります。

実際に、解くにあたっては、 $H_{ij}, \; S_{ij}$ をそれぞれ計算するわけですが、ちゃんと積分計算していくのはなかなか大変です。手計算では無理そうです。

そこで思い切った近似が必要になります。ここで使えるのがヒュッケル法です。

まず、ヒュッケル法はブタジエンのような π電子共役系で適用できる近似手法です。このような状況で、次のような近似を行います：

ヒュッケル法

$\displaystyle S_{ij} = \int \chi_i^* \chi_j \, d{\bf r} = \begin{cases} 1 & (\; i = j \; ) \\ 0 & (\; i \neq j \;) \end{cases} \tag{35}$

$\displaystyle H_{ij} = \int \chi_i^* \hat{H}_{\text{分子}} \chi_j \, d{\bf r} = \begin{cases} \alpha & (\; i = j \; ) \\ \beta & (\; i \neq j \;\; \text{かつ} \;\; (i, j) \; \text{間に}\,\sigma\,{結合あり} \;) \\ 0 & (\; i \neq j \;\; \text{かつ} \;\; (i, j) \; \text{間に}\,\sigma\,{結合なし} \;) \end{cases} \tag{36}$

ただし、 $\alpha < 0, \;\; \beta < 0$ は定数。

重なり積分 $S_{ij}$ の方はわかりやすいですね。

$S_{ii}$ は $\chi_i$ の全空間での確率なので規格化条件より $1$ であり、 $i\neq j$ の場合は $S_{ij} = 0$ と近似するというわけです。（波動関数同士の「重なり」がほとんどないと解釈できそうです。）

実は、エチレンのケースでも $S_{12} = 0$ としていたわけですが、これはヒュッケル法を反映していたわけですね。

一方の $H_{ij}$ は少し難しいルールですね。 $H_{ii} = \alpha$ の方は、原子核はすべて対等であるということで、定数としているわけです。（ $H_{ii}$ はクーロン積分というのでした。）

$i\neq j$ の場合は、原子核 $i$ と $j$ の間にσ結合があるかないかで値が決まります。今考えているのは、π電子なわけですが、σ結合があるということはその原子核間で結合があるということですね。その場合は $H_{ij} = \beta$ （定数）となります。そうでない場合は、要するに直接つながっていない原子核同士ということなので、その場合は $H_{ij} = 0$ とするわけです。（ $i\neq j$ の $H_{ij}$ を共鳴積分というのでした。）

たとえば、ブタジエンの場合は $1-2$ の間にはσ結合があるので $H_{12} = \beta$ です。一方、 $1-3$ や $1-4$ の間にはσ結合はないので、 $H_{13} = H_{14} = 0$ というわけです。

さて、ヒュッケル法を用いて永年方程式 $(34)$ を計算すると、次のようになります：

$\operatorname{det} \begin{pmatrix} \alpha - E & \beta & 0 & 0 \\ \beta & \alpha - E & \beta & 0 \\ 0 & \beta & \alpha - E & \beta \\ 0 & 0 & \beta & \alpha - E \end{pmatrix} = 0 \tag{37}$

ずいぶんと簡単な、そして対称的な形の行列式になりましたね。

ある１列をある定数 $C$ について $C$ 倍した行列の行列式は、元の行列式の $C$ 倍になります。これを使って、上記の行列の各列をすべて $1/\beta$ 倍すると

$\operatorname{det} \begin{pmatrix} \frac{\alpha - E}{\beta} & 1 & 0 & 0 \\ 1 & \frac{\alpha - E}{\beta} & 1 & 0 \\ 0 & 1 & \frac{\alpha - E}{\beta} & 1 \\ 0 & 0 & 1 & \frac{\alpha - E}{\beta} \end{pmatrix} = 0 \tag{38}$

となり、 $x = \frac{\alpha - E}{\beta}$ とおくと

$\operatorname{det} \begin{pmatrix} x & 1 & 0 & 0 \\ 1 & x & 1 & 0 \\ 0 & 1 & x & 1 \\ 0 & 0 & 1 & x \end{pmatrix} = 0 \tag{39}$

となります。

あとは式 $(39)$ の行列式を計算すれば良いですね！　まさに線形代数の世界です！

なお、計算過程を思い出すと、この行列式は

対角成分が $x$
$(i, j)$ の間にσ結合があれば $1$ 、そうでなければ $0$

というルールによって生成されます。そのため、機械的に導くことができます。

行列式 $(39)$ ですが、手計算はそれなりに大変です。

まず、 $4\times 4$ の行列式をそのまま計算するのは大変なので、第１列を使って余因子展開します。

$\begin{align} &\operatorname{det} \begin{pmatrix} x & 1 & 0 & 0 \\ 1 & x & 1 & 0 \\ 0 & 1 & x & 1 \\ 0 & 0 & 1 & x \end{pmatrix} \\ &= x\cdot \operatorname{det} \begin{pmatrix} x & 1 & 0 \\ 1 & x & 1 \\ 0 & 1 & x \end{pmatrix} - 1\cdot \operatorname{det} \begin{pmatrix} 1 & 0 & 0 \\ 1 & x & 1 \\ 0 & 1 & x \end{pmatrix} \\ &= x(x^3 - x - x) - (x^2 - 1) \\ &= x^4 - 2x^2 - x^2 + 1 \\ &= x^4 - 2x^2 + 1 - x^2 \\ &= (x^2 - 1)^2 - x^2 \\ &= (x^2 + x - 1) (x^2 - x - 1) \end{align}$

これが $0$ に一致するので

$(x^2 + x - 1) (x^2 - x - 1) = 0$

であり、よって

$\displaystyle x = \frac{-1 - \sqrt{5}}{2}, \; \frac{1 - \sqrt{5}}{2}, \; \frac{-1 + \sqrt{5}}{2}, \; \frac{1 + \sqrt{5}}{2} \tag{40}$

の４つが解となります。

$x = \frac{\alpha - E}{\beta}$ だったことを思い出すと

$\displaystyle \frac{\alpha - E}{\beta} = \frac{-1 - \sqrt{5}}{2}, \; \frac{1 - \sqrt{5}}{2}, \; \frac{-1 + \sqrt{5}}{2}, \; \frac{1 + \sqrt{5}}{2}$

であり、よって

$\displaystyle E = \alpha + \beta \frac{1 + \sqrt{5}}{2}, \;\; \alpha + \beta \frac{-1 + \sqrt{5}}{2}, \;\; \alpha + \beta \frac{1 - \sqrt{5}}{2}, \;\; \alpha + \beta \frac{-1 - \sqrt{5}}{2} \tag{41}$

となります。

つまり、分子軌道のエネルギーとしては

$\displaystyle \begin{cases} E_1 = \alpha + \beta \,\frac{1 + \sqrt{5}}{2} \\ E_2 = \alpha + \beta \,\frac{-1 + \sqrt{5}}{2} \\ E_3 = \alpha + \beta \,\frac{1 - \sqrt{5}}{2} \\ E_4 = \alpha + \beta \,\frac{-1 - \sqrt{5}}{2} \end{cases} \tag{42}$

の４つが取り得て、 $\alpha < 0, \; \beta < 0$ を考慮に入れると $E_1 < E_2 < E_3 < E_4$ となることがわかります。

つまり、エチレンの場合は結合性軌道と反結合性軌道の２つがあったわけですが、ブタジエンにおいては４つの分子軌道に分かれるというわけですね！

これらの４つに対応する軌道は、連立方程式

$\begin{pmatrix} x & 1 & 0 & 0 \\ 1 & x & 1 & 0 \\ 0 & 1 & x & 1 \\ 0 & 0 & 1 & x \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \\ c_3 \\ c_4 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \tag{43}$

に対して $E = E_1, E_2, E_3, E_4$ に対応する $x = \frac{\alpha - E}{\beta}$ を代入して計算できます。

加えて $\Psi = c_1 \chi_1 + c_2 \chi_2 + c_3 \chi_3 + c_4 \chi_4$ の規格化の条件

$c_1^2 + c_2^2 + c_3^2 + c_4^2 = 1 \tag{44}$

を適用すると、 $E_i$ に対応する $c_1^{(i)}, c_2^{(i)}, c_3^{(i)}, c_4^{(i)}$ がそれぞれの $i = 1, 2, 3, 4$ で求まります。

計算は飛ばしますが、対応する分子軌道の波動関数は

$\begin{cases} \Psi_1 = 0.372\, \chi_1 + 0.602\, \chi_2 + 0.602\, \chi_3 + 0.372\, \chi_4 \\ \Psi_2 = 0.602\, \chi_1 + 0.372\, \chi_2 - 0.372\, \chi_3 - 0.602\, \chi_4 \\ \Psi_3 = 0.602\, \chi_1 - 0.372\, \chi_2 - 0.372\, \chi_3 + 0.602\, \chi_4 \\ \Psi_4 = -0.372\, \chi_1 + 0.602\, \chi_2 - 0.602\, \chi_3 + 0.372\, \chi_4 \end{cases} \tag{45}$

となります。

ようやく分子軌道が計算できたので、この波動関数をプログラムで表現し可視化してみましょう。

ブタジエンの座標を図のようにおきます。

これを踏まえて、ブタジエンの波動関数を計算する関数を以下のように定義します：

# 記事の一番最後に完全版のプログラムを掲載します。
# ブタジエンのπ電子共役系
# http://www.nishino-labo.jp/pdf/enshu01_03_2012.pdf （原子間距離の参考）
def butadiene_molecule(x, y, z, id):
    Z = 6        # 炭素原子を想定しているので Z = 6
    Z_eff = Z - (0.35*3 + 0.85*2)  # 有効核電荷（スレータ―則により計算）
    n = 2
    l = 1
    m = 0
    distance_1 = 147.0/52.9  # 水素原子のボーア半径 52.9 pm, C2-C3の原子間距離 147.0 pm
    distance_2 = 134.0/52.9  # 水素原子のボーア半径 52.9 pm, C1-C2, C3-C4の原子間距離 134.0 pm
    
    # 各炭素原子の中心座標
    center = [
        [0.5 * distance_1 * np.cos(7.0*np.pi/6.0), 0.5 * distance_1 * np.sin(7.0*np.pi/6.0) - distance_2],
        [0.5 * distance_1 * np.cos(7.0*np.pi/6.0), 0.5 * distance_1 * np.sin(7.0*np.pi/6.0)],
        [0.5 * distance_1 * np.cos(np.pi/6.0), 0.5 * distance_1 * np.sin(np.pi/6.0)],
        [0.5 * distance_1 * np.cos(np.pi/6.0), 0.5 * distance_1 * np.sin(np.pi/6.0) + distance_2],
    ]
    
    # 係数
    c = [
        [0.372, 0.602, 0.602, 0.372],
        [0.602, 0.372, -0.372, -0.602],
        [0.602, -0.372, -0.372, 0.602],
        [-0.372, 0.602, -0.602, 0.372],
    ]
    
    res = 0.0
    for i in range(4):
        x0 = x - center[i][0]
        y0 = y - center[i][1]
        res += f(x0, y0, z, n, l, m, Z_eff) * c[id][i]

これを可視化したのが次の図になります：

左上が $\Psi_1$ 、右上が $\Psi_2$ 、左下が $\Psi_3$ 、右下が $\Psi_4$ です。

エネルギー図に表すと次のようになります。

さて、ここでπ電子の行き先について考えましょう。

ブタジエンの各炭素原子の ${\rm p}_z$ 軌道から全部で４つの電子が生じますが、これらがエネルギーが低い軌道から各２つずつ入っていきます。したがって、 $\Psi_1$ （エネルギーは $E_1$ ）に２つの電子が入り、 $\Psi_2$ （エネルギーは $E_2$ ）に２つの電子が入ります。

したがって、電子は最大で $E_2$ のところまでしか入らないというわけですね。 $E_3, \; E_4$ には空の軌道ということになります。

ここで一つ重要な用語を紹介したいのですが、電子によって占有されている最もエネルギーが高い軌道のことを HOMO（Highest Occupied Molecular Orbital）といい、電子が空であるような最もエネルギーが低い軌道のことを LUMO（Lowest Unoccupied Molecular Orbital）といいます。

なお、読み方ですが、HOMO（ホモ）・LUMO（ルモ）と呼ぶそうです。

この用語を使うと、ブタジエンのπ電子共役系においては、HOMOが $\Psi_2$ であり、LUMOは $\Psi_3$ ということになります。

HOMO-LUMOを考えるモチベーションとしては、HOMOの電子の反応しやすさが挙げられます。

HOMO-LUMOを起点として、種々の化学反応を説明するフロンティア軌道理論というものがあるそうです。今回は触れることはできませんが、とにかくそういう考え方があるのですね。

光の吸収という観点でもHOMO-LUMOを考えるのは重要です。

日曜化学(1)の記事では、こんなことを書いていました。水素原子に振動数 $\nu$ の光が入ってきたとき、電子がより高いエネルギーの軌道（励起状態）に遷移することがあります。遷移する条件はエネルギー差が $h\nu$ にちょうど一致することです。これは分子軌道でも同じことです。

すなわち、ちょうどエネルギー差が $\Delta E = E_3 - E_2$ に一致する光が入射すると、HOMOの電子が光を吸収してLUMOに励起するというわけです。

エチレン（二重結合１つ）とブタジエン（二重結合２つ）を考えましたが、同じ要領で二重結合／単結合のセットを増やしていくと、分子軌道の数が $2, 4, 6, 8, \ldots$ と増えていきます。

図の右に行けば行くほど、HOMOとLOMOの間のエネルギー差はどんどん小さくなっていきます。よりエネルギーの小さい（つまり、長い波長）の光を吸収するようになるということに他なりません。

こんな風に、分子軌道のエネルギーというのは、吸収する光の波長と大変関係が深いということですね。吸収する光の波長が可視光であれば、物質の色にも影響を与えることになります。

たとえば、β-カロテンという物質を考えましょう。ニンジンの中に入っている橙色の色素成分です。これはβ-カロテンの分子が補色である青色の光を吸収するからなのですが、この青色を吸収する仕組みが分子軌道からわかるのです。

β-カロテンの分子は、以下のような構造をしています：

二重結合が全部で11個あります。ブタジエンの議論の延長で、炭素22個によるπ電子共有系を考えると、HOMO/LUMOのエネルギー

$E_{\text{HOMO}} = \alpha + 0.136 \, \beta$

$E_{\text{LUMO}} = \alpha - 0.136 \, \beta$

であることが計算できます。したがって、エネルギー差は

$\Delta E = E_{\text{LUMO}} - E_{\text{HOMO}} = -0.272 \, \beta$

となります。これはエチレンのHOMO-LUMO間のエネルギー差の $1/7$ 以下です。より低い振動数の（つまり長波長の）光を吸収できるようになるというわけですね。これによって可視光の波長を吸収することができるようになるわけです。

こんな風に分子軌道法によって、我々が見ているものの「色」を分子軌道から説明できる可能性があるのです。これは面白くてたまらないですね！

これまで私自身のモチベーションをはぐらかしてきましたが、そろそろ白状すると、私の目的は色の理解です。

次回（最終回）の記事では、これまでの量子化学の理論を踏まえて、「とある身近な色」の仕組みについて考えていきたいと思います。

7. ベンゼンのπ共有系

最後の例としてベンゼンを扱いたいと思います。

ベンゼン $\text{C}_6\text{H}_6$ は二重結合と単結合が交互に並んで、円環状の構造を持つ分子です。

ベンゼン分子も二重結合を含むので、π電子が関係します。

二重結合と単結合が交互になっているからといって、π電子が各二重結合に局在化しているわけではありません。ブタジエンで考えたように、π電子が分子全体に広がって分布するようになります。こういう系をπ電子共役系というのでした。

ベンゼン分子のπ電子共役系の分子軌道を計算してみましょう。図の上にある炭素原子を $1$ として、そこから時計回りに $1, 2, 3, 4, 5, 6$ という名前をつけましょう。

対応する ${\rm p}_z$ 軌道の原子軌道を $\chi_1, \chi_2, \chi_3, \chi_4, \chi_5, \chi_6$ として、その重ね合わせとして波動関数

$\Psi = c_1 \chi_1 + c_2 \chi_2 + c_3 \chi_3 + c_4 \chi_4 + c_5 \chi_5 + c_6 \chi_6$

を考えます。対応するエネルギー（の期待値）は計算できて６変数 $c_1, c_2, c_3, c_4, c_5, c_6$ の関数になるわけですが、この極値を求めると分子軌道が得られるというわけですね。

ここまで丁寧に計算してきましたので、今回は途中の議論を省略して考えたいと思います。

まず、ヒュッケル法を用いて永年方程式を計算すると、次のようになります。

$\operatorname{det} \begin{pmatrix} \alpha - E & \beta & 0 & 0 & 0 & \beta \\ \beta & \alpha - E & \beta & 0 & 0 & 0 \\ 0 & \beta & \alpha - E & \beta & 0 & 0 \\ 0 & 0 & \beta & \alpha - E & \beta & 0 \\ 0 & 0 & 0 & \beta & \alpha - E & \beta \\ \beta & 0 & 0 & 0 & \beta & \alpha - E \end{pmatrix} = 0 \tag{46}$

ベンゼンは環構造になっているので、 $1$ 番の炭素と $6$ 番の炭素がつながっていることに注意しましょう。（行列の1-6成分と6-1成分に $\beta$ があることに注意。）

すべての列を $\beta$ で割って $x = \frac{\alpha - E}{\beta}$ とおくと

$\operatorname{det} \begin{pmatrix} x & 1 & 0 & 0 & 0 & 1 \\ 1 & x & 1 & 0 & 0 & 0 \\ 0 & 1 & x & 1 & 0 & 0 \\ 0 & 0 & 1 & x & 1 & 0 \\ 0 & 0 & 0 & 1 & x & 1 \\ 1 & 0 & 0 & 0 & 1 & x \end{pmatrix} = 0 \tag{47}$

となります。あとはこれを解けば分子軌道のエネルギーが得られます。

この計算は大変そうだなと思っていたら、木村巌先生（ @iwaokimura ）からSagemathを使えば計算できると教えていただいたので、これを使っていきます。

$x$ を固有値だと思うと、上の行列式は固有多項式（characteristic polynomial）だと思うことができます。これを使って、Sagemathで上記の行列から $xI$ を引いた行列 $A$

$A = \begin{pmatrix} 0 & 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 1 \\ 1 & 0 & 0 & 0 & 1 & 0 \end{pmatrix} \tag{48}$

について、固有多項式を計算します。

A = Matrix(QQ, [[0, 1, 0, 0, 0, 1], [1, 0, 1, 0, 0, 0], [0, 1, 0, 1, 0, 0], [0, 0, 1, 0, 1, 0], [0, 0, 0, 1, 0, 1], [1, 0, 0, 0, 1, 0]]); A
p = B.characteristic_polynomial(); p
factor(p)

実行すると、次が得られます。

[0 1 0 0 0 1]
[1 0 1 0 0 0]
[0 1 0 1 0 0]
[0 0 1 0 1 0]
[0 0 0 1 0 1]
[1 0 0 0 1 0]
x^6 - 6*x^4 + 9*x^2 - 4
(x - 2) * (x + 2) * (x - 1)^2 * (x + 1)^2

つまり

$\begin{align} &\operatorname{det} \begin{pmatrix} x & 1 & 0 & 0 & 0 & 1 \\ 1 & x & 1 & 0 & 0 & 0 \\ 0 & 1 & x & 1 & 0 & 0 \\ 0 & 0 & 1 & x & 1 & 0 \\ 0 & 0 & 0 & 1 & x & 1 \\ 1 & 0 & 0 & 0 & 1 & x \end{pmatrix} \\ &= x^6 - 6x^4 + 9x^2 - 4 \\ &= (x - 2) (x + 2) (x - 1)^2 (x + 1)^2 \end{align}$

ということですね。あっという間に計算できてしまいました。

$x = \frac{\alpha - E}{\beta}$ を思い出して、これが $\text{det} = 0$ なる解を求めると、次の４つが解となります。（重複度も含めると６つ。）

$E = E_1 = \alpha + 2\beta$
$E = E_2 = \alpha + \beta$ （重複度：2）
$E = E_3 = \alpha - \beta$ （重複度：2）
$E = E_4 = \alpha - 2\beta$

エネルギーの大小関係は、 $\alpha < 0, \; \beta < 0$ を考慮して $E_1 < E_2 < E_3 < E_4$ となります。

対応する波動関数も、固有ベクトルを計算することで求めることができます。

A.eigenvectors_left()

Sagemathでこれを実行すると、こうなります（見やすいように少し整形しています）：

[
   (2, [(1, 1, 1, 1, 1, 1)], 1), 
   (-2, [(1, -1, 1, -1, 1, -1)], 1), 
   (1, [(1, 0, -1, -1, 0, 1), (0, 1, 1, 0, -1, -1)], 2), 
   (-1, [(1, 0, -1, 1, 0, -1), (0, 1, -1, 0, 1, -1)], 2)
]

出力の形式としては、以下のようになっています：

[
    (固有値1, [(固有ベクトル11, 固有ベクトル12, ..., 固有ベクトル1n),  固有値1の重複度), 
    (固有値2, [(固有ベクトル21, 固有ベクトル22, ..., 固有ベクトル2m),  固有値2の重複度), 
    ...
]

固有値 $x = 2$ （エネルギーは $E_1 = \alpha - 2\beta$ ）に対応する固有ベクトルは

$(c_1, c_2, c_3, c_4, c_5, c_6) = (1, 1, 1, 1, 1, 1)$

ということになります。

正確にいうと、定数 $C$ をかけた $(C, C, C, C, C, C)$ が、この固有値に対応する固有ベクトル全体となります。言い換えると、これがエネルギー $E_1$ の波動関数全体（の空間）ということになります。これを

$c_1^2 + c_2^2 + c_3^2 + c_4^2 + c_5^2 + c_6^2 = 1 \tag{49}$

によって規格化すると、分子軌道

$\displaystyle \Psi_1 = \frac{1}{\sqrt{6}} \chi_1 + \frac{1}{\sqrt{6}} \chi_2 + \frac{1}{\sqrt{6}} \chi_3 + \frac{1}{\sqrt{6}} \chi_4 + \frac{1}{\sqrt{6}} \chi_5 + \frac{1}{\sqrt{6}} \chi_6 \tag{50}$

が得られます。

ここで注意したいのは、重複度が２以上のときです。重複度 $n \geq 2$ の固有値に対応する固有ベクトルの次元は $n$ 次元になるということを考慮する必要があるからです。（重複度が2であれば2次元。）

たとえば、固有値 $x = 1$ （エネルギーは $E_2 = \alpha - \beta$ ）に対応する固有ベクトルの次元は $2$ となるわけですが、２次元なので基底が２個となります。実際

$v_1 = (1, 0, -1, -1, 0, 1)$
$v_2 = (0, 1, 1, 0, -1, -1)$

が基底となります。固有ベクトル全体の空間は

$C_1 v_1 + C_2 v_2$

となります。

ところで、今たまたま上のような基底 $\langle v_1, v_2\rangle$ をとっていますが、基底の取り方は特に決まっているわけではありません。つまり、基底を取り直しても何ら問題ないということです。

ここではひとまず

$\displaystyle v_1 + \frac{1}{2}v_2, \;\;\; v_1$

を基底として選ぶことにしましょう。（特に意味はないのですが、参考文献の教科書の波動関数に合わせるためです。）

規格化すると、対応する波動関数は

$\begin{cases} \Psi_2 = \frac{1}{\sqrt{3}} \chi_1 + \frac{1}{2\sqrt{3}} \chi_2 - \frac{1}{2\sqrt{3}} \chi_3 - \frac{1}{\sqrt{3}} \chi_4 - \frac{1}{2\sqrt{3}} \chi_5 + \frac{1}{2\sqrt{3}} \chi_6 \\ \Psi_3 = \frac{1}{2} \chi_2 - \frac{1}{2} \chi_3 - \frac{1}{2} \chi_5 - \frac{1}{2} \chi_6 \end{cases} \tag{51}$

となります。

こんな要領で、各固有値に対応する６つの分子軌道が次のように得られます：

$\begin{cases} \Psi_1 = \frac{1}{\sqrt{6}} \chi_1 + \frac{1}{\sqrt{6}} \chi_2 + \frac{1}{\sqrt{6}} \chi_3 + \frac{1}{\sqrt{6}} \chi_4 + \frac{1}{\sqrt{6}} \chi_5 + \frac{1}{\sqrt{6}} \chi_6 \\ \Psi_2 = \frac{1}{\sqrt{3}} \chi_1 + \frac{1}{2\sqrt{3}} \chi_2 - \frac{1}{2\sqrt{3}} \chi_3 - \frac{1}{\sqrt{3}} \chi_4 - \frac{1}{2\sqrt{3}} \chi_5 + \frac{1}{2\sqrt{3}} \chi_6 \\ \Psi_3 = \frac{1}{2} \chi_2 + \frac{1}{2} \chi_3 - \frac{1}{2} \chi_5 - \frac{1}{2} \chi_6 \\ \Psi_4 = \frac{1}{\sqrt{3}} \chi_1 - \frac{1}{2\sqrt{3}} \chi_2 - \frac{1}{2\sqrt{3}} \chi_3 + \frac{1}{\sqrt{3}} \chi_4 - \frac{1}{2\sqrt{3}} \chi_5 - \frac{1}{2\sqrt{3}} \chi_6 \\ \Psi_5 = \frac{1}{2} \chi_2 - \frac{1}{2} \chi_3 + \frac{1}{2} \chi_5 - \frac{1}{2} \chi_6 \\ \Psi_6 = \frac{1}{\sqrt{6}} \chi_1 - \frac{1}{\sqrt{6}} \chi_2 + \frac{1}{\sqrt{6}} \chi_3 - \frac{1}{\sqrt{6}} \chi_4 + \frac{1}{\sqrt{6}} \chi_5 - \frac{1}{\sqrt{6}} \chi_6 \end{cases} \tag{52}$

さて、これにてベンゼンのπ電子共役系の分子軌道がすべて求まりました。いよいよ分子軌道を可視化したいと思います！

まず、ベンゼンの炭素原子 $i = 1, 2, 3, 4, 5, 6$ の座標は

$\displaystyle (x, y) = \left(\sin \frac{2\pi (i-1)}{6}, \;\; \cos \frac{2\pi (i-1)}{6}\right) \tag{53}$

のように設定します。

これを踏まえて分子軌道の波動関数を以下のように定義します。

# 記事の一番最後に完全版のプログラムを掲載します。
# ベンゼン分子のπ電子共役系
def benzene_molecule(x, y, z, id):
    Z = 6        # 炭素原子を想定しているので Z = 6
    Z_eff = Z - (0.35*3 + 0.85*2)  # 有効核電荷（スレータ―則により計算）
    n = 2
    l = 1
    m = 0
    distance = 139/52.9  # 水素原子のボーア半径 52.9 pm, ベンゼン分子の炭素間距離 139 pm
    
    # 各炭素原子の中心座標
    center = []
    for i in range(6):
        cx = distance * np.sin(2*i*np.pi/6)
        cy = distance * np.cos(2*i*np.pi/6)
        center.append([cx,cy])
    
    # 係数
    c = [
        [1/np.sqrt(6), 1/np.sqrt(6), 1/np.sqrt(6), 1/np.sqrt(6), 1/np.sqrt(6), 1/np.sqrt(6)],
        [1/np.sqrt(3), 1/(2.0*np.sqrt(3)), -1/(2.0*np.sqrt(3)), -1/np.sqrt(3), -1/(2.0*np.sqrt(3)), 1/(2.0*np.sqrt(3))],
        [0.0, 1/2.0, 1/2.0, 0.0, -1/2.0, -1/2.0],
        [1/np.sqrt(3), -1/(2.0*np.sqrt(3)), -1/(2.0*np.sqrt(3)), 1/np.sqrt(3), -1/(2.0*np.sqrt(3)), -1/(2.0*np.sqrt(3))],
        [0.0, 1/2.0, -1/2.0, 0.0, 1/2.0, -1/2.0],
        [1/np.sqrt(6), -1/np.sqrt(6), 1/np.sqrt(6), -1/np.sqrt(6), 1/np.sqrt(6), -1/np.sqrt(6)],
    ]
    
    res = 0.0
    for i in range(6):
        x0 = x - center[i][0]
        y0 = y - center[i][1]
        res += f(x0, y0, z, n, l, m, Z_eff) * c[id][i]
    
    return res

これを可視化したものが次の６つの図となります。

ベンゼンのエネルギー図を書くと次のようになります。

炭素原子の ${\rm p}_z$ 軌道から計 $6$ つの電子が供給されるので、これが分子軌道にエネルギーが低い方から入っていきます。したがって、 $E_2$ までの軌道が充填されHOMOとなり、 $E_3$ がLUMOとなるわけですね。

8. 分子軌道の可視化のためのPythonプログラム

最後に、これまで可視化に使ってきたプログラムを紹介して終わりにします。

日曜化学シリーズを通して、Pythonのmatplotlibを使った可視化を行っています。したがって、以下のライブラリを事前にインストールする必要があります。

numpy
matplotlib

また、可視化方法としてmatplotlibではなく、Plotlyというライブラリを使ったバージョンを使いたい方は、Plotlyもインストールしておいてください。
（今回の記事の図は、すべてPlotlyを使用しています。）

可視化においては、第２回（日曜化学(2)：3次元空間における電子雲の計算（Python/matplotlib））の記事内で紹介した方法３を用いています。

得られた分子軌道の波動関数 $\Psi(x, y, z)$ に対して、確率密度関数 $|\Psi(x, y, z)|^2$ の値が高い順に点を表示していき、確率が 90% になった段階でやめる、という方法を使って可視化しています。

ただし、分子軌道の計算にあたって方法３と少し変えた部分があり、その点を注意しておきます。

ヒュッケル法を適用するときに、 $S_{ij} = 0$ と近似していました。この影響で「分子軌道の確率密度関数 $|\Psi(x, y, z)|^2$ を全空間で積分したときの値が $1$ に一致しない」という問題が発生します。

分子軌道 $\Psi = c_1 \chi_1 + c_2 \chi_2$ を全空間で積分すると

$\displaystyle \int |\Psi(x, y, z)|^2 dxdydz = c_1^2 + c_2^2 + 2c_1 c_2 S_{12} = 1 \tag{54}$

となります。

$S_{12}$ は本当は $0$ ではないですが、 $S_{12} = 0$ として、 $c_1^2 + c_2^2 = 1$ となるように係数 $c_1, c_2$ を決めるのがヒュッケル法でした。ところが、こうしてしまうと $c_1^2 + c_2^2$ の時点で $1$ になってしまうので、確率全体としては $2c_1 c_2 S_{12}$ の分だけずれが生じるわけですね。そのため、全空間の確率が $1$ を下回ったり、 $1$ を超えてしまったりするわけです。

これは前回紹介した可視化方法「確率90%の領域を描画する」においてクリティカルな問題となります。全空間での確率の総和が90%を下回ってしまう場合が起きてしまうからです。（たとえば、エチレンの $\pi_{2{\rm p}}^*$ 軌道の絶対値の二乗を $-5 \leq x, y, z \leq 5$ の範囲で積分すると、 $0.729$ 程度の積分値になります。）

この状況下でアルゴリズムをそのまま実行した結果、全領域を描画してしまうことになります。結果的に全域が真っ赤（真っ青）になってしまうのです。

それでは困りますね。今回は、とりあえずの対策として「あらかじめ全空間の積分を計算しておき、その90%の領域を描画する」というようにプログラムを修正しました。一応、それっぽい描画にはなるので、ひとまずの解決としていますが、私の中で納得はいっていないので後々修正するかもしれません。

# 「90%サンプリング」による分子軌道の可視化
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as axes3d

import math


# 動径波動関数（ボーア半径を a_0 = 1 に正規化）
def radial_wave_function(r, n, l, Z_eff):
    rho = 2.0 * Z_eff * r / n
    if n == 1 and l == 0:
        return (Z_eff ** 1.5) * np.exp(-rho/2.0) * 2.0
    elif n == 2 and l == 0:
        return (Z_eff ** 1.5) * np.exp(-rho/2.0) * (2.0 - rho) / (2.0*np.sqrt(2.0))
    elif n == 2 and l == 1:
        return (Z_eff ** 1.5) * np.exp(-rho/2.0) * rho / (2.0*np.sqrt(6.0))
    elif n == 3 and l == 0:
        return (Z_eff ** 1.5) * np.exp(-rho/2.0) * (6.0 - 6.0*rho + rho**2) / (9.0*np.sqrt(3.0))
    elif n == 3 and l == 1:
        return (Z_eff ** 1.5) * np.exp(-rho/2.0) * (4.0 - rho) * rho / (9.0*np.sqrt(6.0))
    elif n == 3 and l == 2:
        return (Z_eff ** 1.5) * np.exp(-rho/2.0) * (rho**2) / (9.0*np.sqrt(30.0))


# 球面調和関数（ただし、実関数表示したもの）
def spherical_harmonics(theta, phi, l, m):
    if l == 0:
        if m == 0:
            # l=0, m=0
            return np.sqrt(1.0/(4*np.pi))
    if l == 1:
        if m == 0:
            # l=1, m=0
            return np.sqrt(3.0/(4.0*np.pi)) * np.cos(theta)
        if m == 1:
            # l=1, m=+1
            return np.sqrt(3.0/(4.0*np.pi)) * np.sin(theta) * np.cos(phi)
        if m == -1:
            # l=1, m=+1
            return np.sqrt(3.0/(4.0*np.pi)) * np.sin(theta) * np.sin(phi)
    if l == 2:
        if m == 0:
            return np.sqrt(5.0/(16.0*np.pi)) * (3.0*(np.cos(theta)**2) - 1.0)
        if m == 1:
            return np.sqrt(15.0/(4.0*np.pi)) * np.cos(theta) * np.sin(theta) * np.cos(phi)
        if m == -1:
            return np.sqrt(15.0/(4.0*np.pi)) * np.cos(theta) * np.sin(theta) * np.sin(phi)
        if m == 2:
            return np.sqrt(15.0/(16.0*np.pi)) * (np.sin(theta)**2) * np.cos(2*phi)
        if m == -2:
            return np.sqrt(15.0/(16.0*np.pi)) * (np.sin(theta)**2) * np.sin(2*phi)


        
# 水素原子の波動関数 Ψ_{n,l,m}(x, y, z)
def f(x, y, z, n, l, m, Z_eff):
    # 座標系を (x,y,z) -> (r,θ,φ) に変換
    r = np.sqrt(x*x + y*y + z*z)
    
    theta = 0.0
    phi = 0.0
    if r > 0:
        theta = np.arccos( z / r )
    if y == 0:
        if x < 0:
            phi = np.pi
    elif x*x+y*y > 0:
        phi = np.sign(y) * np.arccos(x/np.sqrt(x*x+y*y))
        
    # 動径波動関数と球面調和関数の積を計算して出力
    return (radial_wave_function(r, n, l, Z_eff) * spherical_harmonics(theta, phi, l, m))


    
# 水素分子のσ_1s軌道
def hydrogen_molecule(x, y, z, id):
    Z = 1        # 炭素原子を想定しているので Z = 1
    n = 1
    l = 0
    m = 0
    distance = 2.0  # 水素分子の原子核間距離
    if id == 0:
        # σ_1s軌道（結合性軌道）
        return (f(x-distance/2, y, z, n, l, m, Z)/np.sqrt(3.2)) + (f(x+distance/2, y, z, n, l, m, Z)/np.sqrt(3.2))
    else:
        # σ_1s*軌道（反結合性軌道）
        return (f(x-distance/2, y, z, n, l, m, Z)/np.sqrt(0.8)) - (f(x+distance/2, y, z, n, l, m, Z)/np.sqrt(0.8))



# エチレン分子のπ_2p軌道
def ethylene_molecule(x, y, z, id):
    Z = 6        # 炭素原子を想定しているので Z = 6
    Z_eff = Z - (0.35*3 + 0.85*2)  # 有効核電荷（スレータ―則により計算）
    n = 2
    l = 1
    m = 0
    distance = 133.9/52.9  # 水素原子のボーア半径 52.9 pm, エチレン分子の原子間距離 133.9 pm
    if id == 0:
        # π_2p軌道（結合性軌道）
        return (f(x-distance/2, y, z, n, l, m, Z_eff)/np.sqrt(2.0)) + (f(x+distance/2, y, z, n, l, m, Z_eff)/np.sqrt(2.0)) 
    else:
        # π_2p*軌道（反結合性軌道）
        return (f(x-distance/2, y, z, n, l, m, Z_eff)/np.sqrt(2.0)) - (f(x+distance/2, y, z, n, l, m, Z_eff)/np.sqrt(2.0)) 


    
# ブタジエンのπ電子共役系
# http://www.nishino-labo.jp/pdf/enshu01_03_2012.pdf （原子間距離の参考）
def butadiene_molecule(x, y, z, id):
    Z = 6        # 炭素原子を想定しているので Z = 6
    Z_eff = Z - (0.35*3 + 0.85*2)  # 有効核電荷（スレータ―則により計算）
    n = 2
    l = 1
    m = 0
    distance_1 = 147.0/52.9  # 水素原子のボーア半径 52.9 pm, C2-C3の原子間距離 147.0 pm
    distance_2 = 134.0/52.9  # 水素原子のボーア半径 52.9 pm, C1-C2, C3-C4の原子間距離 134.0 pm
    
    # 各炭素原子の中心座標
    center = [
        [0.5 * distance_1 * np.cos(7.0*np.pi/6.0), 0.5 * distance_1 * np.sin(7.0*np.pi/6.0) - distance_2],
        [0.5 * distance_1 * np.cos(7.0*np.pi/6.0), 0.5 * distance_1 * np.sin(7.0*np.pi/6.0)],
        [0.5 * distance_1 * np.cos(np.pi/6.0), 0.5 * distance_1 * np.sin(np.pi/6.0)],
        [0.5 * distance_1 * np.cos(np.pi/6.0), 0.5 * distance_1 * np.sin(np.pi/6.0) + distance_2],
    ]
    
    # 係数
    c = [
        [0.372, 0.602, 0.602, 0.372],
        [0.602, 0.372, -0.372, -0.602],
        [0.602, -0.372, -0.372, 0.602],
        [-0.372, 0.602, -0.602, 0.372],
    ]
    
    res = 0.0
    for i in range(4):
        x0 = x - center[i][0]
        y0 = y - center[i][1]
        res += f(x0, y0, z, n, l, m, Z_eff) * c[id][i]
        
    return res

# ベンゼン分子のπ電子共役系
def benzene_molecule(x, y, z, id):
    Z = 6        # 炭素原子を想定しているので Z = 6
    Z_eff = Z - (0.35*3 + 0.85*2)  # 有効核電荷（スレータ―則により計算）
    n = 2
    l = 1
    m = 0
    distance = 139/52.9  # 水素原子のボーア半径 52.9 pm, ベンゼン分子の炭素間距離 139 pm
    
    # 各炭素原子の中心座標
    center = []
    for i in range(6):
        cx = distance * np.sin(2*i*np.pi/6)
        cy = distance * np.cos(2*i*np.pi/6)
        center.append([cx,cy])
    
    # 係数
    c = [
        [1/np.sqrt(6), 1/np.sqrt(6), 1/np.sqrt(6), 1/np.sqrt(6), 1/np.sqrt(6), 1/np.sqrt(6)],
        [1/np.sqrt(3), 1/(2.0*np.sqrt(3)), -1/(2.0*np.sqrt(3)), -1/np.sqrt(3), -1/(2.0*np.sqrt(3)), 1/(2.0*np.sqrt(3))],
        [0.0, 1/2.0, 1/2.0, 0.0, -1/2.0, -1/2.0],
        [1/np.sqrt(3), -1/(2.0*np.sqrt(3)), -1/(2.0*np.sqrt(3)), 1/np.sqrt(3), -1/(2.0*np.sqrt(3)), -1/(2.0*np.sqrt(3))],
        [0.0, 1/2.0, -1/2.0, 0.0, 1/2.0, -1/2.0],
        [1/np.sqrt(6), -1/np.sqrt(6), 1/np.sqrt(6), -1/np.sqrt(6), 1/np.sqrt(6), -1/np.sqrt(6)],
    ]
    
    res = 0.0
    for i in range(6):
        x0 = x - center[i][0]
        y0 = y - center[i][1]
        res += f(x0, y0, z, n, l, m, Z_eff) * c[id][i]
    
    return res





# サンプリングの設定
N = 50                    # サンプリング個数
x_range = 5               # 座標の範囲 (-x_range, x_range)
delta = 2.0*x_range / N   # サンプリング間隔
#print(delta)


# 3次元空間を格子状にサンプリング（間隔 delta, N^3 個）
data_list = []
for i in range(N):
    for j in range(N):
        for k in range(N):
            x = i*delta - x_range
            y = j*delta - x_range
            z = k*delta - x_range
            
            # 座標 (x,y,z) における波動関数 Ψ(x,y,z) を計算
            # psi = hydrogen_molecule(x,y,z, 0)   # id = 0 or 1
            # psi = ethylene_molecule(x,y,z, 1)   # id = 0 or 1
            # psi = butadiene_molecule(x,y,z, 0)   # id = 0, 1, 2, 3
            psi = benzene_molecule(x,y,z, 0)   # id = 0, 1, 2, 3, 4, 5
            
            # 波動関数の2乗によって電子の確率密度関数 |Ψ(x,y,z)|^2 を計算
            pdf = psi**2
            
            # 確率密度(|Ψ(x,y,z)|^2, p.d.f), x, y, z, 波動関数の生の値 Ψ(x,y,z) の順にデータを格納
            data_list.append([pdf, x, y, z, psi])


# 確率密度が高い順にサンプリングデータをソートする
sorted_list = sorted(data_list, key=lambda data: data[0], reverse=True)  # 確率密度(p.d.f)でソート



# 波動関数が正のデータを格納する変数
x1_list = []
y1_list = []
z1_list = []

# 波動関数が負のデータを格納する変数
x2_list = []
y2_list = []
z2_list = []



# 現在考えている全体の空間（体積: (x_range)^3）における積分値を計算しておく（あとで使う）
def integral_in_all_spaces():
    integral = 0.0
    i = 0
    while i < len(sorted_list):
        pdf = sorted_list[i][0]
        psi = sorted_list[i][4]

        # |Ψ(x,y,z)|^2 ΔV を加算
        integral += (pdf * delta**3)
        i += 1
    print("グリッド数: {}, 全空間の積分: {}".format(i, integral) )  
    return integral

all_integral = integral_in_all_spaces()



# |Ψ|^2 ΔV の総和( volume )が 0.9 を超えるまで足し合わせる（体積が0.9を超えたら終了）
integral = 0.0
i = 0
while integral < 0.9 * all_integral: # 全空間の積分値に対して 90 % になる領域を描画するように変更
    pdf = sorted_list[i][0]
    x = sorted_list[i][1]
    y = sorted_list[i][2]
    z = sorted_list[i][3]
    psi = sorted_list[i][4]
    
    # |Ψ(x,y,z)|^2 ΔV を加算
    integral += (pdf * delta**3)
    
    if psi > 0:
        x1_list.append(x)
        y1_list.append(y)
        z1_list.append(z)
    else:
        x2_list.append(x)
        y2_list.append(y)
        z2_list.append(z)   
    i += 1



# figureを生成する
X1 = np.array(x1_list)
Y1 = np.array(y1_list)
Z1 = np.array(z1_list)

X2 = np.array(x2_list)
Y2 = np.array(y2_list)
Z2 = np.array(z2_list)



#'''
# matplotlibによる可視化
fig = plt.figure(figsize=(8.0, 8.0))
ax = fig.add_subplot(1,1,1, projection='3d')

ax.set_xlim3d(-x_range, x_range)
ax.set_ylim3d(-x_range, x_range)
ax.set_zlim3d(-x_range, x_range)

# 軸の設定
ax.set_xlabel("x", size = 14)
ax.set_ylabel("y", size = 14)
ax.set_zlabel("z", size = 14)

ax.plot(X2,Y2,Z2,color='b',marker="o",linestyle='None')
ax.plot(X1,Y1,Z1,color='r',marker="o",linestyle='None')

plt.show()
#'''



'''
# Plotlyによる可視化
import plotly.graph_objects as go

fig = go.Figure(data=[go.Scatter3d(
    x = X1,
    y = Y1,
    z = Z1,
    mode='markers',
    marker = dict(
        size = 2,
        color = 'red'
    )
), go.Scatter3d(
    x = X2,
    y = Y2,
    z = Z2,
    mode='markers',
    marker = dict(
        size = 2,
        color = 'blue'
    )
)])

fig.update_layout(
    scene = dict(
        camera=dict(eye=dict(x=1.7, y=-1.7, z=1.7)),  #eye=dict(x=1.25, y=1.25, z=1.25)), 
        xaxis = dict(range=[-x_range, x_range]),
        yaxis = dict(range=[-x_range, x_range]),
        zaxis = dict(range=[-x_range, x_range]),
        aspectmode='cube', #this string can be 'data', 'cube', 'auto', 'manual'
        #a custom aspectratio is defined as follows:
        aspectratio=dict(x=1, y=1, z=0.95)
    ),
    width=800,
    height=700,
)

fig.show()
'''

プログラム内の次の箇所を変えると（コメントアウトを変更すると）、表示する分子軌道の種類を変えることができます。（デフォルトはベンゼンの $\Psi_1$ 軌道。）

            # 座標 (x,y,z) における波動関数 Ψ(x,y,z) を計算
            # psi = hydrogen_molecule(x,y,z, 0)   # id = 0 or 1
            # psi = ethylene_molecule(x,y,z, 1)   # id = 0 or 1
            # psi = butadiene_molecule(x,y,z, 0)   # id = 0, 1, 2, 3
            psi = benzene_molecule(x,y,z, 0)   # id = 0, 1, 2, 3, 4, 5

また、デフォルトでは、matplotlibで描画するようになっていますが、matplotlibの箇所をコメントアウトして、コメントアウトされているplotlyの箇所を戻せばplotlyで可視化できます。（plotlyを使用する際には、plotlyのインストールが別途必要です。）