こんにちは。
今回は、「python でヒュッケル法の計算をする」の第2回です。
前回、ヒュッケル法をプログラムにし、分子軌道関数の計算ができるようになりました。
ですがそれを図にする方法はまだ触れていませんでした。
ということで、今回は分子軌道の可視化を行ってみます。
目次
概要
「分子軌道の可視化」とは言いましたが、電子にははっきりとした軌道 (orbit) があるわけではなく、この辺りに居る確率が高いぞというところを軌道 (orbital) と呼んでいるため、実際は電子の存在確率を可視化することになります。
（ぼんやりしているので、「電子雲」なんて呼ばれたりもします。）
実はこのブログの初期の方で、シュレーディンガー方程式から計算される原子軌道を描画する、ということをやっていました。
このときは軌道関数そのものを図にしているので、存在確率の描画ではありません。
したがって、別の方法で図にしなければなりません。
さて、ヒュッケル法では共役系のπ電子のみを考えているので、炭素原子はsp2混成であり、p軌道が大きく関わってきます。
前回の記事では、分子の骨格を描いて各炭素原子の位置に係数 $c$ の符号や大きさに対応したp軌道を配置していく、という方法で分子軌道の概形を図示しました。
例えばベンゼンなら、
$$ \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} $$
このように図示できます。
ここで、図中のp軌道の白や黒は各原子関数の正負と対応しています。
しかし、実際の分子軌道関数は原子軌道関数の和で表されると考えているため、正と正、負と負の間は強め合い、正と負の間は弱め合うようになり、純粋なp軌道の集まりではなくなります。
そんな図をこれから描こうというのが今回の目標です。
（純粋なp軌道の集まりと考えたものが上の図に相当します。）
（百聞は一見に如かずということで、先に完成図を示しておきます。）
方針
先程、分子軌道の描画は電子の存在確率の描画であることを言いました。
ではどうやって存在確率を計算するのかということになりますが、これは分子軌道関数の2乗を積分することで求められます。
つまり、適当な範囲での積分ならその範囲における電子の存在確率になりますし、全範囲でなら 1 となります。
（前回触れた、規格化条件が全範囲の場合に相当します。）
今回は、一定の範囲内での存在確率という形ではなく、ある一点での存在確率 (確率密度？) として図とすることとします。
（また、その計算で得られる数値にちゃんとした意味は無いと考えてください。）
こちらなどでは一定範囲での存在確率の形で表されています。
というわけで今回行う描画方法をざっくり説明すると、
描画範囲内からランダムに一点選ぶ
そこでの電子の存在確率を $\psi ^{2}$ から計算する
一定以上の確率がある点のみ描画する
となります。
ここで、分子軌道関数 $\psi$ の計算には原子軌道関数 $\phi$ が必要になります。
先程p軌道が重要になると書きましたが、その内のひとつである、
$$ \phi = \frac{ 1 }{ 4 \sqrt{ 2} } \left( \frac{ Z }{ a_{0} } \right) ^{ \frac{3}{2} } \sigma \mathrm{e} ^{ -\frac{\sigma}{2} } \mathrm{cos} \theta , ~~~ \sigma = \frac{ Zr }{ a_{0} } $$
という関数を使います。
ただし、正確性は求めていないので係数は考えずに
$$ \phi =r \mathrm{e} ^{ -\frac{r}{2} } \mathrm{cos} \theta $$
と簡略化することにします。
では、実際にプログラムにしていきます。
プログラム
分子軌道関数の計算
まずは前回つくった、ヒュッケル法の計算部分のコードを挙げておきます。
（何故わざわざ scipy を使っているかなど、詳細は前回の記事を参照してください。）
import numpy as np from scipy import linalg n_carbon = 4 # 炭素原子の数 (4ではブタジエンとなる) A = np.zeros((n_carbon, n_carbon)) # n行n列の零行列を生成 # 永年行列式の作成 for i in range(0, n_carbon-1): A[i][i+1] = -1 A += A.T # 転置したものを足す print(A) X, c = linalg.eigh(A) # 固有値・固有ベクトルの計算と保存 c = c.T # 係数を保存した行列を転置 # 出力 print(f"X = {X}") print(f"c = {c}")
c.T を用いて転置することで、c[0] を分子軌道関数 $\psi _{1}$ の係数、c[1] を $\psi _{2}$ の係数、…と見やすくできました。
実行すると、係数として
[[ 0.37174803 0.60150096 0.60150096 0.37174803]
[-0.60150096 -0.37174803 0.37174803 0.60150096]
[ 0.60150096 -0.37174803 -0.37174803 0.60150096]
[ 0.37174803 -0.60150096 0.60150096 -0.37174803]]
が出力され、ブタジエンの分子軌道関数は
$$ \begin{align} \psi_{1} = 0.372 \phi_{1} + 0.602 \phi_{2} + 0.602 \phi_{3} + 0.372 \phi_{4} \\ \psi_{2} = 0.602 \phi_{1} + 0.372 \phi_{2} - 0.372 \phi_{3} - 0.602 \phi_{4} \\ \psi_{3} = 0.602 \phi_{1} - 0.372 \phi_{2} - 0.372 \phi_{3} + 0.602 \phi_{4} \\ \psi_{4} = 0.372 \phi_{1} - 0.602 \phi_{2} + 0.602 \phi_{3} - 0.372 \phi_{4} \end{align} $$
と求められました。
炭素原子の座標の計算
続いて、炭素原子の座標を計算します。
今はブタジエンを想定しているので、
（実際には二重結合と単結合で結合長は違いますが、簡単のため同じ dist という長さだとして考えることにします。）
構造最適化をすればいちいち計算しなくても最も安定な配座での炭素の座標が得られるのでしょうが、とりあえず今はこの形から手で計算することにします。
（ブタジエンならこの形が最安定だと思います。）
この図から炭素原子の座標を計算し、ベクトル化すると、
bond_length = 500 pos_carbon = [] pos_carbon.append(np.array([-bond_length/2-bond_length*np.cos(np.pi/3), -bond_length*np.sin(np.pi/3), 0])) pos_carbon.append(np.array([-bond_length/2, 0, 0])) pos_carbon.append(np.array([bond_length/2, 0, 0])) pos_carbon.append(np.array([bond_length/2+bond_length*np.cos(np.pi/3), bond_length*np.sin(np.pi/3), 0]))
都合上、左（もしくは右）の炭素から順番に配列 pos_carbon に入れてください。
存在確率 (確率密度) の計算
では、電子の存在確率を計算していきます。
まずは描画範囲からランダムに点を選びます。
length = 1000 #選ぶ範囲の一辺の長さ n_point = 10000 # 選ぶ点の個数 rng = np.random.default_rng() x = rng.uniform(-length, length, n_point) y = rng.uniform(-length, length, n_point) z = rng.uniform(-length, length, n_point)
続いて、分子軌道関数に代入して値を計算します。
psi = [] num_psi = 0 # 描画したい分子軌道関数の番号 pos = [] for i in range(n_point): pos.append(np.array([x[i], y[i], z[i]])) # 選んだ点をベクトルに変換 dist = [] theta = [] phi = 0 for j in range(n_carbon): dist.append(np.linalg.norm(pos[i]-pos_carbon[j])) # 各炭素原子と選んだ点との間の距離 theta.append(np.arccos(pos[i][2]/dist[j])) # 各炭素原子から見たときの選んだ点の方向 phi += c[num_psi][j]*dist[j]*np.exp(-dist[j]/100)*np.cos(theta[j]) # 原子軌道関数を足していく psi.append(phi) # 分子軌道関数
そして、これが必要か定かではありませんが、一応全ての $\psi ^{2}$ の和に対する割合として存在確率を計算することにします。
# psi^2の和 s = 0 for i in range(n_point): s += psi[i]*psi[i] # ある点での存在確率 for i in range(n_point): prob = psi[i]*psi[i]/s
描画のために、$\psi$ の値の正負によってその点での座標を別々の配列に入れていきます。
x_p = [] y_p = [] z_p = [] x_n = [] y_n = [] z_n = [] for i in range(n_point): prob = psi[i]*psi[i]/s # 存在確率 # 確率が低い点は描画しない if prob < 0.00001: continue # psiの正負によって分ける if psi[i] > 0: x_p.append(pos[i][0]) y_p.append(pos[i][1]) z_p.append(pos[i][2]) if psi[i] < 0: x_n.append(pos[i][0]) y_n.append(pos[i][1]) z_n.append(pos[i][2])
if prob < 0.00001 の 0.00001 はテキトーなのでお好みで変えてみてください。
分子軌道の描画
各点での存在確率が計算や $\psi$ の値による選り分けができたので、これを描画します。
ここでは matplotlib というライブラリを使います。
（インストールは pip install matplotlib でできます。）
from matplotlib import pyplot as plt # 3次元でのプロット fig = plt.figure(figsize=(10, 10)) ax = fig.add_subplot(projection='3d') # psiが正の点は赤く、負の点は青く描く ax.plot(x_p, y_p, z_p, color='r', marker="o", linestyle='None') ax.plot(x_n, y_n, z_n, color='b', marker="o", linestyle='None')
どこに炭素原子があるかも描いておきます。
# 各炭素原子の座標をx,y,zで分けて配列に入れる x_carbon = [] y_carbon = [] z_carbon = [] for i in range(n_carbon): x_carbon.append(pos_carbon[i][0]) y_carbon.append(pos_carbon[i][1]) z_carbon.append(pos_carbon[i][2]) ax.plot(x_carbon, x_carbon, x_carbon, color='k', marker="o") # 黒で描画
最後に、できた図を表示するために、
plt.show()
をつけます。
以上で、プログラムが完成しました。
全てまとめると、
import numpy as np from scipy import linalg from matplotlib import pyplot as plt n_carbon = 4 # 炭素原子の数 (4ではブタジエンとなる) A = np.zeros((n_carbon, n_carbon)) # n行n列の零行列を生成 # 永年行列式の作成 for i in range(0, n_carbon-1): A[i][i+1] = -1 A += A.T # 転置したものを足す # print(A) X, c = linalg.eigh(A) # 固有値・固有ベクトルの計算と保存 c = c.T # 係数を保存した行列を転置 # 出力 # print(f"X = {X}") # print(f"c = {c}") # 炭素原子の座標 bond_length = 500 pos_carbon = [] pos_carbon.append(np.array([-bond_length/2-bond_length*np.cos(np.pi/3), -bond_length*np.sin(np.pi/3), 0])) pos_carbon.append(np.array([-bond_length/2, 0, 0])) pos_carbon.append(np.array([bond_length/2, 0, 0])) pos_carbon.append(np.array([bond_length/2+bond_length*np.cos(np.pi/3), bond_length*np.sin(np.pi/3), 0])) # ランダムに選ぶ length = 1000 #選ぶ範囲の一辺の長さ n_point = 10000 # 選ぶ点の個数 rng = np.random.default_rng() x = rng.uniform(-length, length, n_point) y = rng.uniform(-length, length, n_point) z = rng.uniform(-length, length, n_point) # 分子軌道関数の計算 psi = [] num_psi = 0 # 描画したい分子軌道関数の番号 pos = [] for i in range(n_point): pos.append(np.array([x[i], y[i], z[i]])) # 選んだ点をベクトルに変換 dist = [] theta = [] phi = 0 for j in range(n_carbon): dist.append(np.linalg.norm(pos[i]-pos_carbon[j])) # 各炭素原子と選んだ点との間の距離 theta.append(np.arccos(pos[i][2]/dist[j])) # 各炭素原子から見たときの選んだ点の方向 phi += c[num_psi][j]*dist[j]*np.exp(-dist[j]/100)*np.cos(theta[j]) # 原子軌道関数を足していく psi.append(phi) # 分子軌道関数 # psi^2の和 s = 0 for i in range(n_point): s += psi[i]*psi[i] # 描画のための準備 x_p = [] y_p = [] z_p = [] x_n = [] y_n = [] z_n = [] for i in range(n_point): prob = psi[i]*psi[i]/s # 存在確率 # 確率が低い点は描画しない if prob < 0.00001: continue # psiの正負によって分ける if psi[i] > 0: x_p.append(pos[i][0]) y_p.append(pos[i][1]) z_p.append(pos[i][2]) if psi[i] < 0: x_n.append(pos[i][0]) y_n.append(pos[i][1]) z_n.append(pos[i][2]) # 3次元でのプロット fig = plt.figure(figsize=(10, 10)) ax = fig.add_subplot(projection='3d') # psiが正の点は赤く、負の点は青く描く ax.plot(x_p, y_p, z_p, color='r', marker="o", linestyle='None') ax.plot(x_n, y_n, z_n, color='b', marker="o", linestyle='None') # 各炭素原子の座標をx,y,zで分けて配列に入れる x_carbon = [] y_carbon = [] z_carbon = [] for i in range(n_carbon): x_carbon.append(pos_carbon[i][0]) y_carbon.append(pos_carbon[i][1]) z_carbon.append(pos_carbon[i][2]) ax.plot(x_carbon, y_carbon, z_carbon, color='k', marker="o") # 黒で描画 plt.show()
これを実行すると、
（ここでは num_psi = 0 なので、$\psi _{1}$ が描画されています。）
num_psi を 1~3 まで変えて実行してみると、
それっぽいです。
（ただし仕様上 plot は上に重ねられていくため、赤と青の前後が入れ替わってしまっている部分もあります。）
欲を言えばちゃんと面として描画したいですが、これはこれで「電子雲」らしくてよいのではないでしょうか。
gifとして保存
さて、一応描画はできましたが、もっといろんな角度から見たいです。
視点を変える方法として、
ax.view_init(elev=20, azim=0)
を追加して、elev や azim の値を変えていく方法がありますが、どうせなら動かしたい。
そんなときは、gif にすれば良い。
詳しくはこちらを参照してください。
大して難しいことはしておらず、視点を変えた画像を何枚もつくってくっつけているだけです。
というわけで、一番初めに載せたような gif のできあがりです。
グリグリ動かすのも一応方法はあるらしいですが、自分の環境ではうまくいきませんでした。
（動かせたので追記しました。(2023/10/07)）
追記（2022/12/19）
描画する点をランダムに取っていましたが、格子状に取るようにしてみました。
隙間が埋まって、よりそれらしくなったのではないでしょうか。
プログラム上は、「ランダムに選ぶ」の部分を以下のコードに置き換えればOKです。
# 格子状に区切る length = 1000 l = 50 x = [] y = [] z = [] for i in range(-length, length, l): for j in range(-length, length, l): for k in range(-length, length, l): x.append(i) y.append(j) z.append(k) n_point = len(x)
うまくいかないときは、l や描画する prob の範囲を変えてみてください。
追記（2023/10/07）
約1年ぶり、今更感は否めないですが、追記です。
図を出力してもグリグリ動かせなく、gif にするしかなかったわけですが、環境を変えたらできたというお話です。
これまでは「JupyterLab」というソフトウェアを使っていました。
というのも、python は他のソフトを入れない限りコマンドプロンプトで動かすものだと思っていたので、開発環境をどこかから拾ってこないとやりづらそうだと考えていたからです。
しかし、python をダウンロードすると「IDLE」 なるものも一緒に入ってきて、実はそれが開発環境として普通に使える、ということに最近気づきました。
試しにとこの記事で作ったコードを動かしてみると、普通にグリグリ動かせました。
まとめ
以上で、ヒュッケル法の計算結果を python で描画することができるようになりました。
前回も書きましたが、同じことを processing でもできるようにしたり、拡張ヒュッケル法も実装したり、ということも考えています。
お楽しみに。
