自由電子模型法

十進BASICインタプリタ(Windows, UNIX, Mac)はダウンロード可能.Macの最新OSでも実行可能(2021年3月).Pythonは,macOS Catalina の場合,バージョン2.7と3.7がプリインストールされている.

プログラミングしながらFMOを学習する

有機化学は憶える学問と言う人がいる.確かに, 憶えなければならないことが多い.ところが, 語学に象徴されるように, ほとんどの学問分野で最小限のことは憶えないと次の進展が望めない. したがっ て, 有機化学を記憶に頼る学問と決めつけるのは少々早計のようである. ここではBASICやPythonプログラ ミング(時代の要請を考慮し追記)を利用して, 分子軌道の概略を知り, 計算により算出された数値を使って反応のメカニズムを 知り, さらに生成物を予測することを試みてみよう. 注 ) プログラミングの部分は飛ばしていただいても結構です.

環化付加反応について

1928年,DielsとAlderは次式に示すような反応を発見し論文として発表した.当時,反応はプラス の部分とマイナスの部分が結合するという概念でとらえられていたが,この反応は中性物質同士の反 応であり,+電荷と-電荷を用いて説明することはできない(例 R+ + X- → RX).しかも,ただ加熱 するだけで炭素と炭素が二カ所で同時に結合するという大変不思議な反応であり,分子の骨核を一挙 に構築する反応として1950年ノーベル化学賞を受賞した.

十数年後,本反応はまったく新しい理論(1965)によって説明されることとなった.1981年我が国初 のノーベル化学賞に輝いた

福井謙一博士の「フロンティア軌道理論 [Frontier Molecular Orbital (FMO) Theory]」

である.

フロンティア 軌道 (FMO) 法

エチレンは12個の価電子を有している. 10個は骨核を作るために単結合に使用され,残った2個の電子は二重結合に使われる. 1個の軌道に入ることができる電子の数は2個だから, エネルギー の低い軌道から順番に電子を入れていくと下から数えて6個目の軌道でおしまいになる. この軌道は 最高被占軌道 (HOMO) と呼ばれ, 電子を与える軌道として最も活性である.

一方,電子の入っていない空の軌道のうち一番低い軌道を最低空軌道 (LUMO) と呼ぶ. 電子を受け 入れる軌道として機能する.これらの軌道を総称してフロンティア軌道(FMO) とよぶ.

今回のような付加反応では, 反応に関与するのは二重結合だけであるから, σ軌道のなすxy平面から垂直に突き出したπ軌道だけを取り出して議論してもよい. そうするとエチレンの反応は図のような2個のπ電子だけ を考慮すればよいことになる.

ブタジエンやヘキサトリエン等の共役二重結合を有する分子でも同様に考えることができる.

注)軌道の重なりの型と結合の種類については後述する.

エチレンの分子軌道

分子軌道は分子を構成する原子軌道の線形結合によって作られる.すなわち, 2個のp軌道からπ 軌道ができる際の式は次式で表される(数学的には2次方程式の解として得られる).

ψ1=c11+c22

ψ2=c11-c22

前式はエチレンのHOMO,後式はエチレンのLUMOに当たる. HOMOに入っている電子2個がエネルギーをもらってもっと高い軌道に遷移する際,受け入れる空の軌道と考えればよい.

自由電子模型法

新しく出来た分子軌道に電子がどのように分布しているかは簡単な波の式で表わすことができる. 定常波と呼ばれ, 波の両端は「節」になっている. Sinカー ブの180度までの半周期の「弧」を考えればよい. 1番目の軌道は弧が1個(節0個),2番目の軌道は弧(中心に節1個)が2個である.弧が1.5個など中途半端な波は存在しない.

作図で求めたエチレンの軌道

詳細はブタジエンの分子軌道を参照のこと. 波の高さ(振幅)は二乗した値の合計が1になるように決めればよい. 振幅を二乗した値は電子の存在確率すなわち電子密度を表す. 1個の電子の電子密度は全空間にわたって積分した値は1であることからも容易に理解できる. したがってエチレンの軌道係数は0.707である. 安定な軌道はψ1=0.707φ1+0.707φ2, 不安定な軌道はψ2=0.707φ1-0.707φ2で表される.

ブタジエンの分子軌道

上記の手続きは一般式で表すことができる.

ψn (x) = (2/L)1/2sin(n*pi*x/L)

n=1, 2, 3, 4 ----

ブタジエン (CH2=CH-CH=CH2)で説明しよう.

まず, 4個の炭素原子を任意の距離で等間隔に置く.

起点 C1-----C2-----C3-----C4 終点

|------|------|------|------|------|

|<--------------- L -------------->|

|<------x---->|

次に両端の炭素原子から炭素ー炭素の距離分だけ延長した点同士の距離をLとする.起点から任意の 炭素原子までの距離をxとする.nは弧の数である.

実際に計算してみよう. 炭素ー炭素距離は任意だから1と置くとL=5である.

弧1個の場合 n=1

ψ1 (x) = (2/5)1/2sin(1*pi*x/5)となり

x 1 2 3 4

ψ1 (x) 0.372 0.601 0.602 0.372

である.同様に各軌道の各炭素上の振幅(係数)が計算できる.

n=4 0.372 -0.602 0.601 -0.372 (空軌道)

n=3 0.602 -0.372 -0.372 0.601 (空軌道)

n=2 0.601 0.372 -0.372 -0.602 電子2個(被占軌道)

n=1 0.372 0.601 0.602 0.372 電子2個(被占軌道)

一連の手続きは「自由電子模型」と呼ばれている有名な方法である.

軌道の形と位相について

二重結合のp軌道の形は球を上下に団子状にした形で表される. 水素などのs軌道は球型,さらに d軌道は鉢巻型などs, p, d軌道固有の形をしているので注意が必要である.

水素の分子軌道

エチレンの軌道と異なる点はs-軌道のため球状になっただけ

p-軌道とd-軌道の形

軌道は波の形をしているため, 波が正の方に振れる場合と負の方に振れる場合がある. これを単純に表すため上下の団子の色付けをして区別する.

左上図のような形をした軌道に電子が入った場合, A,C,Eの各点では電子の存在確率はゼロである. B, Dの点で電子密度は最大になるが, 位相が異なるため2個のp軌道の性質(反応性)はまったく異なる. 同じ位相を持つ軌道は結合するが, 位相が異なる場合は反発する.

右上の図はエチレン2分子が熱反応 (HOMOとLUMOが相互作用する) によりシクロブタンを生成しない理由を示している. 一方では安定化が起きるが,他方では不安定になり, 全体としては安定化は起きない.

このことは電子が波の性質を有していることから容易に理解できる. すなわち, 2個の波が重なるとき, お互いの波の位相が合えば強め合って振幅の大きな波になるが, 位相が異なる波が重なると打ち消し合い波は減衰する.

BASIC言語によるプログラミング

ここまで説明したら, 賢明なる諸君のことだから, 情報処理学などで習ったBASIC言語を思い出すだろう (一昔前の薬学部学生向けの記述). さっそく自由電子模型のプログラムを作ってみ よう.下記のプログラムは三角関数は「度」で考え,途中でラジアンに変換している.

十進BASICでは,角の大きさの単位は標準ではラジアンであるが, OPTION ANGLE DEGREES をプログラムのはじめに書くことで角度の単位を度 (degrees) に変えることができる. 60行は不要であり, 70行は wh = sin(a)に変更する. 実際の実行では文番号を書く必要はない.

10 input n 共役分子の炭素の数を入力

20 for i = n to 1 step -1 弧の数の多い方から計算 120行とループ

30 print i; 軌道の順を表示

40 for j = 1 to n 各原子に沿って計算処理 100行とループ

50 a = i*j*180/(n+1) 50-80行は係数の計算

60 d = 3.141592*a/180 度をラジアンへ換算

70 wh = sin(d)

80 coeff = sqr(2/(n+1))*wh

90 print using "###.###";coeff; 結果の表示

100 next j

110 print

120 next i

OPTION ANGLE DEGREES(十進BASICの場合)

INPUT n

for i = n to 1 step -1

print i;

for j = 1 to n

LET a = i*j*180/(n+1)

LET wh = SIN(a)

LET coeff = sqr(2/(n+1))*wh

print using "###.###":coeff;

next j

print

next i

END

>run

? 4

4 0.372 -0.602 0.601 -0.372

3 0.602 -0.372 -0.372 0.601 ---- LUMO

2 0.601 0.372 -0.372 -0.602 ---- HOMO

1 0.372 0.601 0.602 0.372

計算結果を図示すると分子軌道図ができあがる. ブタジエンのHOMOの係数は両端で大きく膨らんでおり,これらの炭素原子の反応性が高いことが分かる. ま た, 軌道は下から, 対称 (S),非対称 (A)の繰り返しになっていることおよび波の「節」が0,1,2,3個の順に増えていることに気づいてほしい.

十進BASICをMacで実行するには → 現在, OSに関係なく利用可能(2021)

1)VirtualBoxをインストールしてWindows XP OS等で実行する.

2)Intel Macで実行できる版(Windows版より不安定)を利用する.

Python言語によるプログラミング(最近追記)

最近はパソコン上でPythonが利用できる.BASIC言語でプログラミングしたFEM法の基本部分をPython言語に書き換えてみた.Python特有の文法の習得が必要である.


#FEM Method

import math

x=input('Number of pi electrons?')

n=int(x)

for i in range(n,0,-1):

for j in range(1,n+1):

a = i*j*180/(n+1)

wh = math.sin(math.radians(a))

c = math.sqrt(2/(n+1))*wh

print('{:<8,.3f}'.format(c), end=' ')

print('')


実行結果

Python 3.7.0 (v3.7.0:1bf9cc5093, Jun 26 2018, 23:26:24)

[Clang 6.0 (clang-600.0.57)] on darwin

Type "copyright", "credits" or "license()" for more information.

>>>

==================== RESTART: /Users/KH/pythonPRG/FEM.py ====================

Number of pi electrons?4

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

>>>

熱反応におけるFMO相互作用

熱反応はHOMO-LUMO相互作用で起こり, 光反応はHOMO-HOMOあるいはLUMO LUMO相互作用で起こる.

熱反応がHOMO-HOMO相互作用で起こるとすると,反応前後のエネルギー収支は下図のように描ける.

青矢印の安定化と赤矢印の安不定化が打ち消し合って, 反応系全体としては安定化しないので反応は起こらない.

HOMO-LUMO相互作用では両方の相互作用とも安定化に寄与する.

分子Aと分子Bが相互作用すると,2種類のHOMO-LUMO相互作用が可能であるが, 軌道準位間隔小さい方が主相互作用となる.

HOMO(A)とLUMO(B)の軌道準位間隔は, HOMO(B)とLUMO(A)の軌道準位間隔より小さいので, 電子の安定化(青)も大きい. その理由については, 別項(反応性を予測する計算式)を参照してほしい.

反応におけるFMO相互作用

光反応では, HOMOの電子1個が励起されてLUMOに移る. ここでは電子1個しか存在しない半占軌道をHOMO', LUMO`と記したが, SOMOと書く場合もある. 励起状態の分子A*と分子Aが相互作用すると, ペアリングしなかった二個の電子はそれぞれ安定化, 不安定化で相殺する. ペアリングした電子の分だけ反応系は安定化する.

反応性を予測する計算式

分子同士が接近し反応するには, 反応系全体が安定化する必要である. どの程度安定化エネルギーを獲得するかは次式で計算することができる.

上式はA分子とB分子が一点で反応する際の安定化エネルギーである. Cは係数, βは原子の種類に特有の定数(共鳴積分)である. Eの絶対値が大きくなるためには分母が小さく, 分子が大きくなればよい. 荒っぽい言い方だが, 分子より分母が大切であり, 一方の分子のHOMOと他方のLUMOが近く なれば反応性が大きくなる. 図の実線の相互作用の方が点線の相互作用より有利である. 上式は実線で示した安定化エネルギーを表している. 点線の相互作用も同様に算出できるが, フロンティア軌道論では一般にHOMO-LUMO間隔の小さい方に注目して反応を解析する. 反応点が複数の場合は安定化エネルギーの総和で判断する.

軌道エネルギーの見積もり

分子軌道のエネルギー準位を正確に求めるには高度の分子軌道計算を駆使しなければならないので, ここでは一般的な見積もり方を紹介しておこう.

○電子供与基が置換するとHOMOも LUMOも上昇する.

○電子吸引基が置換するとHOMOも LUMOも下降する.

○共役基が置換するとHOMO-LUMOが接近する.

実例で理解してほしい. 右図から明らかなようにブタジエンとエチレンの反応よりもアクリル酸との反応の方が起こりやすいことが分かる.

ブタジエンとアクロレインのDIels-Alder反応解析

ブタジエンの計算結果を使ってシクロペンタジエンとアクロレインのDiels-Alder反応 (4π+2π反応)を予測してみよう. アク ロレインはエチレンに電子吸引基であるカルボニル基が付いているからHOMO, LUMOとも低下して いる. そのため, シクロペンタジエンのHOMOとアクロレインのLUMOが相互作用することが予想できる. 双方の軌道はブタジエンのHOMO, LUMOで近似してよい.

係数を構造式に書き込んでみると3点(右図では4点)で位相が合い, endo付加が起こることが判る.

exo付加の遷移状態では結合が生成する2カ所だけで位相が合う.

endo付加を有利にする現象は二次軌道効果とよばれている.

ヘキサトリエンの軌道と環化付加反応

プログラムを実行し, ヘキサトリエンのπ電子数6を入力する.

>run

? 6

6 0.232 -0.418 0.521 -0.521 0.418 -0.232 A

5 0.418 -0.521 0.232 0.232 -0.521 0.418 S

4 0.521 -0.232 -0.418 0.418 0.232 -0.521 ---- LUMO A

3 0.521 0.232 -0.418 -0.418 0.232 0.521 ---- HOMO S

2 0.418 0.521 0.232 -0.232 -0.521 -0.418 A

1 0.232 0.418 0.521 0.521 0.418 0.232 S

ヘキサトリエンは6個のπ電子を有しているから, 実行結果の下から3行目の数値がHOMOである. その上の軌道がLUMOである.

例 シクロペンタジエンと無水マレイン酸の4π+2π付加反応

4カ所で位相が合っていることおよび係数の大きい所で反応が起きていることを確認してほしい. 結合する原子以外で位相が合っている相互作用は上述した「二次軌道効果」であり, sndo付加体が生成する.

p軌道の描き方

亜鈴型のp軌道を上から見ると白丸に見える. 下から見ると黒丸に見える.

単に円として描くことがあるので, 内容から判断する必要がある.

共役中員環化合物との反応

ヘキサトリエンの軌道はシクロヘプタトリエンやトロポンなどの反応にも適用可能である.

下記のように単純化して考える.p軌道は円で描いた.

シクロペンタジエンとの6π+4π付加反応では, exo付加体のみが得られる.

Endo付加では右側の構造の付加体が得られるはずである.

しかし, endo付加体は得られない.その理由として二次軌道効果が効かないためと説明されている,

×印で示した2個の相互作用はまったく位相が合わない. そのため反発力が働くと考えられる. 実験事実はフロンティア軌道法による予測と一致する.

フロンティア軌道論で説明できる反応や化学的性質

フロンティア軌道論は電子論とは異なり,試薬の攻撃する立体的な方向や生成物の予想を可能にするなど有用な理論である. 特に下記の反応予測などで威力を発揮する.

・芳香族の置換反応

・求核置換反応

・共役化合物の安定性(芳香族性)

・周辺環状反応(付加反応,シグマトロピー反応,キレトロピー反応,電子環状反応)

・光反応と熱反応の予測

熱反応および光反応の予測

熱反応はHOMO-LUMOで起る. 下記の2+2反応や4+4反応は加熱しても進行しないことが分かる.

実際には紫外線照射によって起る.

共役化合物の安定性(芳香族性)

分子を二つの部分に分け, 一方をHOMO, 他方をLUMOとして接点における位相関係を調べる.

合わない点がある場合, 電子が分子の周辺に均一に分布することができないため不安定である.

ベンゼンはブタジエンのHOMOとエチレンのLUMOに分割し, 接点の位相を調べると両接点で位相が合うため安定である.

シクロペンタジエノンの場合は, ブタジエンHOMOとエチレンLUMOの接点(赤表示)で位相が合わないため不安定である.

構造化学(平面性)

フロンティア軌道相互作用が分子の構造を決めることさえある.

ベンゼンはブタジエンのHOMOとエチレンのLUMOに分割すると両接点で位相が合うため, 電子がスムースに流れるため安定である.

一方シクロオクタテトラエンの場合, シクロヘキサトリエンHOMOとエチレンLUMOに分割すると片方の接点で位相が合わないため, 電子が周辺部位に局在化するため不安定である.位相が合わないのを避けるため, タブ(桶型)構造をとる.

分子内Diels-Alder反応を起こし, 6員環と4員環が縮環した化合物へ変化する.

シクロヘキサトリエンも同様な反応挙動を示し, 三員環と縮環した化合物を与える.

トロポンはシクロヘキサトリエンのHOMOの両端がエチレンLUMOの一端と位相が合い安定である. 単結晶X線解析により平面であることが知られている.

トロポンは平面構造

求核置換反応(SN2反応)

ハロゲン化合物の求核置換反応でYがバックサイトから攻撃すると位相が合う.

前面からの攻撃では, 一方とは位相が合うが 他方とは位相が合わない(波線)..

参考資料

軌道の接近方向

2個のp軌道が接近して相互作用して結合ができる際, 接近の仕方には方向性がある.

最も安定化 が得られる接近の仕方は枠に示した2通りであり. それぞれ二重結合と単結合を形成する.

四角に囲った軌道はπ軌道とσ軌道である.

軌道の重なり型と結合の種類

上から H-H H-C C-C C=Cである. それぞれに結合性軌道 HOMO, 反結合性軌道 LUMOがある.

芳香族の置換反応(ニトロ化)

ニトロ化の起こる位置の予想. NO2+イオンの攻撃はHOMO係数の大きい所に起こる.

◯や●はs軌道ではない. p軌道を上から見た図である.

電子環状反応機構

単結合のσ軌道が二重結合のπ軌道と相互作用することにより,単結合が開裂し,ジエンを与える反応.その際,単結合の回転方向は位相が合うように決められる.

3,3-シグマトロピー反応機構

π軌道とσ軌道の組み合わせでブタジエンのHOMOにあたる軌道を考える. その軌道とエチレンのLUMOの分子内相互作用で見事に説明できる(少し考える必要がある). 分子内環化付加反応と考えることもできる. 詳細は別稿「軌道の混合」を参照.

注)ブタジエンのHOMOは2個のエチレンのπ軌道が位相が合わないように接する形である(上図中程). CH2=CH-CH2-CH2のHOMOはエチレンのHOMOとC-C単結合のHOMOの組合せを考えればよい(上図右側).

電荷移動錯体の構造(キノンーハイドロキノンの分子複合体)

ハイドロキノンとその酸化体は安定な分子化合物を形成する.

配座決定

前者は不対電子とHOMOとの反発,後者は反結合性により遠ざかるように分子が配座を変える.

コンピュータなしで分子軌道を作図する方法

コンピュータがあれば分子の軌道係数を算出したり, 計算結果をグラフィック表示して眺めることにより, 一見して反応を理解することができる.しかし,コンピュータがなければできないことではな い. ここでは軌道混合の法則をもとに簡単な軌道を作図する方法を説明する.

軌道混合則

1.エネルギー順位の異なる軌道が混合する際, エネルギーの高い軌道が低い軌道に混じる際は,位相が合うように少し混じる.

2.エネルギー順位の異なる軌道が混合する際, エネルギーの低い軌道が高い軌道に混じる際は,位相が合わないように少し混じる

エチレンの軌道

一般に,相互作用する2個の軌道エネルギーが接近していればいるほど, 大きく分裂した2個の分子軌道を作 る.

その際.エネルギー的に低い軌道は位相が合うよう に結合する. もう一方の軌道は位相が合わないように結合する.2個のp軌道からエチレンの分子軌道ができる のに相当する.

カルボニルの軌道

炭素と酸素の場合は,酸素の方が炭素より軌道エネル ギーが低いので, 分子軌道ができる際には1:1で混じることはない. エネルギー差に反比例して混合がおこる. カ ルボニルのHOMOは酸素の方に電子が偏っており, LUMOは炭素の係数が大きく,HCNの付加反応等を説明 することができる.

2個のエチレンからブタジエンを作図する

今度は2個のエチレンの軌道からブタジエンの 軌道を作ってみよう. 考え方はエチレンと全く同 じである.

1.近い軌道同士 (HOMO同士, LUMO同士) からそれぞれ1対の結合性軌道および反結合性軌道をつくる.

2.次に,離れた軌道の混合を考慮して係数の大 きさを修正する.すなわち, HOMO同士から一対の軌道を合成する際には相手のLUMO軌道の影響を位相があうように少し加算する. 一方,LUMO 同士から一対の軌道を作る際は,相手のHOMO軌 道を反結合的に少し加算する.

3.全ての軌道にこの考え方を適用すると, 計算で求めた軌道とそっくりの分子軌道ができあがる. 次図を参考にして理解してほしい.

係数の大きさが修正される様子

最初にエチレンのHOMO, LUMOからブタジエンの軌道(1, 2, 3, 4)を作る.

ブタジエンのHOMOはエチレンのHOMOが位相が合わないように繋がる(軌道2).

2の軌道の緑枠の部分とLUMO(赤枠)が位相が合う(青線)ように相互作用する.その結果, 青枠1と青枠2を重ね合わせた青枠3が生成する.

LUMOも同様に考える.HOMOが干渉する場合は位相が合わない様に相互作用する(赤線).

斜めの長い矢印が少し干渉する様子

干渉によって生成した軌道は短い矢印(↑,↓)で示す.

重要な一般則(2個のHOMOから新しい分子のHOMOを作る)

HOMOを作る場合

接点で位相が合わない様 (out of phase) に接続する

LUMOを作る場合

接点で位相が合う様 (in phase) に接続する

例 ヘキサトリエンのHOMO, LUMOをブタジエンとエチレンのHOMO, LUMOから作る.

注 p軌道を円で表現

別稿(軌道の混合)にも詳述しているので,参考にしてほしい.

プログラムの利用

電子密度と結合次数

自由電子模型のプログラムに少し手を加えて,電子密度と結合次数を計算するプログラムを作ってみ よう.計算で求めた係数は一旦配列変数に格納し,取り出して2乗したり掛け合わせる等の演算処理 を行う.

電子密度

各軌道における各原子上の電子密度は係数の2乗を2倍(電子数)したものである.電子の入ってい る軌道について合計したものは一般に言うπ電子密度である.

10 input n

15 dim c(4,4)

20 for i = 4 to 1 step -1

30 print i;

40 for j = 1 to n

50 a = i*j*180/(n+1)

60 d = 3.141592*a/180

70 wh = sin(d)

80 c(i,j) = sqr(2/(n+1))*wh

90 print using "###.###";c(i,j);

100 next j

110 print

120 next i

125 print

130 for i = 4 to 1 step -1

132 print i;

135 for j = 1 to n

140 ed = 2*c(i,j)*c(i,j)

150 print using "###.###";ed;

160 next j

165 print

170 next i

>run<

? 4

4 0.372 -0.602 0.601 -0.372

3 0.602 -0.372 -0.372 0.601

2 0.601 0.372 -0.372 -0.602

1 0.372 0.601 0.602 0.372

4 0.276 0.724 0.724 0.276 ---- 空軌道(電子密度計算には考慮しない)

3 0.724 0.276 0.276 0.724 ---- 空軌道(電子密度計算には考慮しない)

2 0.724 0.276 0.276 0.724 ---- 2個の電子が詰まっている(電子密度計算に考慮する)

1 0.276 0.724 0.724 0.276 ---- 2個の電子が詰まっている(電子密度計算に考慮する)

↓ ↓

1.0 1.0 ---- π電子密度

結合次数

各軌道の結合次数は隣り合う係数を掛け合わせたものを2倍したもでである.位相が合わない場合は 反結合性軌道を反映して結合次数は負になる.電子の入っている軌道につき合計した値に1を加えた ものを全結合次数と言われている.結合次数は結合距離に換算することができる.ブタジエンの2, 3位の結合が1,2位にくらべて弱いことを半定量的に確認できる.

10 input n

15 dim c(4,4)

20 for i = 4 to 1 step -1

30 print i;

40 for j = 1 to n

50 a = i*j*180/(n+1)

60 d = 3.141592*a/180

70 wh = sin(d)

80 c(i,j) = sqr(2/(n+1))*wh

90 print using "###.###";c(i,j);

100 next j

110 print

120 next i

125 print

130 for i = 4 to 1 step -1132 print i;

135 for j = 1 to n-1

140 bo = 2*c(i,j)*c(i,j+1)

150 print using "###.###";bo;

160 next j

165 print

170 next i

>run

? 4

4 0.372 -0.602 0.601 -0.372

3 0.602 -0.372 -0.372 0.601

2 0.601 0.372 -0.372 -0.602

1 0.372 0.601 0.602 0.372

4 -0.447 -0.724 -0.447 ---- 空軌道(結合次数計算には考慮しない)

3 -0.447 0.276 -0.447 ---- 空軌道(結合次数計算には考慮しない)

2 0.447 -0.276 0.447 ---- 2個の電子が詰まっている(結合次数計算に考慮する)

1 0.447 0.724 0.447 ---- 2個の電子が詰まっている(結合次数計算に考慮する)

↓ ↓

0.894 0.448 ---- π結合次数

↓ ↓

1.894 1.448 ---- 全結合次数

Chipmunk BASIC 利用

「軌道係数のグラフィック表示」

上のような簡単なプログラムを利用して

「自由電子模型法で共役ポリエンの軌道係数」を求めるプログラムを以下に示す

塗りつぶしはサークルを半径を変えて描かせている.

5' free electron model

10 input n

15 dim c(n,n)

18 print"Eigenvalues"

20 for i = n to 1 step -1

30 print i;

40 for j = 1 to n

50 a = i*j*180/(n+1)

60 d = 3.141592*a/180

70 wh = sin(d)

80 c(i,j) = sqr(2/(n+1))*wh

90 print using "###.###";c(i,j);

100 next j

110 print

120 next i

125 print

128 print"Electron densities"

130 for i = n to 1 step -1

132 print i;

135 for j = 1 to n

140 ed = 2*c(i,j)*c(i,j)

150 print using "###.###";ed;

160 next j

165 print

170 next i

200 graphics 0

210 for i = n to 1 step -1:'input a$

220 for j = 1 to n step 1

230 c = abs(c(i,j)*30)

235 graphics moveto j*30,300-i*30

236 if c(i,j)<0 then

238 graphics color 0,0,100

250 for r=0 to c:graphics circle r:next r

255 else

258 graphics color 100,0,0

260 for r=0 to c:graphics circle r:next r

270 endif

280 next j

285 graphics moveto j*30,300-i*30+5:graphics drawtext str$(i)

290 next i

300 end

>run

? 4 ブタジエンの計算

Eigenvalues

4 0.372 -0.602 0.602 -0.372

3 0.602 -0.372 -0.372 0.602 ---- LUMO

2 0.602 0.372 -0.372 -0.602 ---- HOMO

1 0.372 0.602 0.602 0.372

Electron densities

4 0.276 0.724 0.724 0.276

3 0.724 0.276 0.276 0.724

2 0.724 0.276 0.276 0.724

1 0.276 0.724 0.724 0.276

>run

? 6 ヘキサトリエンの計算

Eigenvalues

6 0.232 -0.418 0.521 -0.521 0.418 -0.232

5 0.418 -0.521 0.232 0.232 -0.521 0.418

4 0.521 -0.232 -0.418 0.418 0.232 -0.521 ---- LUMO

3 0.521 0.232 -0.418 -0.418 0.232 0.521 ---- HOMO

2 0.418 0.521 0.232 -0.232 -0.521 -0.418

1 0.232 0.418 0.521 0.521 0.418 0.232

Electron densities

6 0.108 0.349 0.543 0.543 0.349 0.108

5 0.349 0.543 0.108 0.108 0.543 0.349

4 0.543 0.108 0.349 0.349 0.108 0.543

3 0.543 0.108 0.349 0.349 0.108 0.543

2 0.349 0.543 0.108 0.108 0.543 0.349

1 0.108 0.349 0.543 0.543 0.349 0.108

十進BASICの利用

OPTION ANGLE DEGREES

input n

SET WINDOW -1,2*n+2,-1,2*n+2

draw grid(1,1)

FOR i = 1 TO n

for j = 0 to n+1 step 0.1

LET a = i*j*180/(n+1)

LET wh = sin(a)

LET c = sqr(2/(n+1))*wh

PLOT LINES: j,i*2+c;

SET LINE COLOR "RED"

IF c<0 THEN

SET LINE COLOR "BLUE"

END IF

IF mod(j,1.0)=0 THEN

FOR k = 0 TO 9 STEP 1

PLOT LINES: j+0.01*k,i*2;j+0.01*k,i*2+c

next k

PLOT LINES

SET LINE COLOR "BLUE"

END IF

NEXT j

PLOT LINES: 0,i*2

next i

END

OPTION ANGLE DEGREES

input n

SET WINDOW -1,2*n+2,-1,2*n+2

draw grid(1,1)

FOR i = 1 TO n

for j = 0 to n+1

LET a = i*j*180/(n+1)

LET wh = sin(a)

LET c = sqr(2/(n+1))*wh

LET r= ABS(C/2)

IF c<0 THEN

SET AREA STYLE "HOLLOW"

DRAW disk WITH SCALE(r)*SHIFT(j,I*2-r)

SET AREA STYLE "SOLID"

DRAW disk WITH SCALE(r)*SHIFT(j,I*2+r)

else

SET AREA STYLE "SOLID"

DRAW disk WITH SCALE(r)*SHIFT(j,I*2-r)

SET AREA STYLE "HOLLOW"

DRAW disk WITH SCALE(r)*SHIFT(j,I*2+r)

END if

NEXT j

next i

END

Pythonプログラミング(追記)

プログラミング言語教育環境の変化を考慮して,最近追記した.

#FEM Method

import math

import matplotlib.pyplot as plt

fig = plt.figure()

ax = fig.add_axes([0, 0, 1.0, 1.0])

n=int(input('Number of pi electrons?:'))

print('Your input data = ', n)

for i in range(n,0,-1):

for j in range(1,n+1):

w = math.sin(i*j*3.141592/(n+1))

c = math.sqrt(2/(n+1))*w

# 中心(x,y)で半径rの円を描画

x=j/(n+1)

y=i/(n+1)

r=abs(c)/8

if c<0:

circle = plt.Circle((x,y+r/2),r/2,fc="#0000FF")

ax.add_patch(circle)

circle = plt.Circle((x,y-r/2),r/2,fc="#0000FF",fill=False)

ax.add_patch(circle)

else:

circle = plt.Circle((x,y+r/2),r/2,fc="#0000FF",fill=False)

ax.add_patch(circle)

circle = plt.Circle((x,y-r/2),r/2,fc="#0000FF")

ax.add_patch(circle)

print('{:<8,.3f}'.format(c), end=' ')

print('')

plt.show()

ヘキサトリエンの実行結果

係数は別ウインドウに表示される.

Scratchプログラミング

小学校で教えているScratch言語を使ってFEMプログラムを作ってみた.結果はリスト(配列)に入れ,スクロールして見る方式である.

参考資料

FBASICによるプログラム例(サポート終了)

上記のプログラムを統合し新たな機能を加えたプログラムを作ってみよう.

追加する新機能

・結合次数から結合距離を算出する

・係数の二乗は電子密度を表すことに着目して円による表示を行う

グラフィック処理コマンドはBASIC言語の種類に依存するので書き換えが必要である.ここではこれまで情報処理学実習で使用してきた富士通の FBASIC(コンパイラ)を使用した.現在薬学部パソコンにインストールされているのは仮称十進BASIC (Windows95環境のインタプリタ,平成11年度熊大総合情報処理センターで採用) であり,グラフィック部分の手直しが必要である.


10 ' fem method PSIn(x)=SQR(2/a)*SIN(n*PI*x/a) 1989 CENTER JISSHU

20 CLS

30 INPUT "number of electrons";N

40 DIM C(N,N)

50 FOR I=1 TO N:FOR J=1 TO N:C(I,J)=0:NEXT J:NEXT I

60 FOR I=N TO 1 STEP -1

70 ' PRINT I;

80 FOR J=1 TO N

90 A=I*J*1/(N+1)

100 D=3.14159*A

110 WH=SIN(D)

120 C(I,J)=SQR(2/(N+1))*WH

130 ' PRINT USING"###.####";COEFF;

140 NEXT J

150 'PRINT

160 NEXT I

170 '------ COEFFICIENT ------

180 PRINT:PRINT"------ COEFFICIENT ------"

190 PRINT"Atom ";:FOR I=1 TO N:PRINT USING"& & ";"C"+MID$(STR$(I),2,2);:NEXT I:PRINT

200 FOR I=N TO 1 STEP -1:PRINT"Orb ";I;

210 FOR J=1 TO N

220 PRINT USING"###.####";C(I,J);

230 NEXT J

240 PRINT

250 NEXT I

260 LINE INPUT A$

270 '------ ELECTRON DENSITY ------

280 PRINT:PRINT"------ ELECTRON DENSITY ------"

290 PRINT" ";:FOR I=1 TO N:PRINT USING"& & ";"C"+MID$(STR$(I),2,2);:NEXT I:PRINT

300 FOR J=1 TO N

310 ED=0

320 FOR I=1 TO N/2

330 D=C(I,J)^2

340 ED=ED+D*2

350 NEXT I

360 PRINT USING"###.####";ED;

370 NEXT J

380 LINE INPUT A$

390 '------ BOND ORDER ------

400 PRINT:PRINT"------ BOND ORDER ------"

410 PRINT" ";:FOR I=1 TO N-1:PRINT USING"& &";"C"+MID$(STR$(I),2,2)+"-"+"C"+MID$(STR$(I+1),2,2);:NEXT I:PRINT

420 FOR J=1 TO N-1

430 BO=0

440 FOR I=1 TO N/2

450 B=C(I,J)*C(I,J+1)

460 BO=BO+B*2

470 NEXT I

480 PRINT USING"###.####";BO+1;

490 NEXT J

500 LINE INPUT A$

510 PRINT:PRINT"------ BOND LENGTH ------"

520 PRINT" ";:FOR I=1 TO N-1:PRINT USING"& &";"C"+MID$(STR$(I),2,2)+"-"+"C"+MID$(STR$(I+1),2,2);:NEXT I:PRINT

530 FOR J=1 TO N-1

540 BO=0

550 FOR I=1 TO N/2

560 B=C(I,J)*C(I,J+1)

570 BO=BO+B*2

580 L=1.515-.177/(1+1.05*(1/BO-1))

590 NEXT I

600 PRINT USING"###.####";L;

610 NEXT J

620 LINE INPUT A$:CLS

630 '------ DRAWING MO ------

640 SYMBOL(150,30),"FEM MOLECULAR ORBITAL",2,2

650 FOR I=N TO 1 STEP -1

660 FOR J=1 TO N

670 DX=400/(N+2):X=150+J*DX

680 DY=360/(N+2):Y=360-I*DY

690 R=C(I,J):K=DY*.8

700 IF ABS(C(I,J))<.001 THEN C(I,J)=0!

710 IF J=1 THEN LINE(X-DX,Y)-(X+N*DX,Y),PSET

720 IF R<0 THEN 740

730 CIRCLE(X,Y),R^2*K,0,,,,F:CIRCLE(X,Y),R^2*K ,7:GOTO 760

740 CIRCLE(X,Y),R^2*K ,7,,,,F

750 ' PRINT C(I,J);

760 SYMBOL(X-20,Y+DY/2),LEFT$(STR$(C(I,J)),5),1,1

770 ' LINE(X-DX,Y)-(X+DX,Y),PSET

780 NEXT J

790 PRINT

800 next I

805 input dummy

810 end

<実行例>

<グラフィック処理>