自由電子模型法

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

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

環化付加反応について

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

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

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

である．

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

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

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

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

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

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

エチレンの分子軌道

分子軌道は分子を構成する原子軌道の線形結合によって作られる．すなわち， ２個のｐ軌道からπ 軌道ができる際の式は次式で表される（数学的には２次方程式の解として得られる）．

ψ1=c1*φ1+c2*φ2

ψ2=c1*φ1-c2*φ2

前式はエチレンのＨＯＭＯ，後式はエチレンのＬＵＭＯに当たる． ＨＯＭＯに入っている電子２個がエネルギーをもらってもっと高い軌道に遷移する際，受け入れる空の軌道と考えればよい．

自由電子模型法

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

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

詳細はブタジエンの分子軌道を参照のこと． 波の高さ（振幅）は二乗した値の合計が１になるように決めればよい． 振幅を二乗した値は電子の存在確率すなわち電子密度を表す． １個の電子の電子密度は全空間にわたって積分した値は１であることからも容易に理解できる． したがってエチレンの軌道係数は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)で説明しよう．

まず， ４個の炭素原子を任意の距離で等間隔に置く．

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

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

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

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

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

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

弧１個の場合 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 電子２個（被占軌道）

n=1  0.372   0.601  0.602   0.372 電子２個（被占軌道）

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

軌道の形と位相について

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

水素の分子軌道

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

p-軌道とd-軌道の形

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

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

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

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

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

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

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

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

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

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

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

反応性を予測する計算式

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

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

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

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

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

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

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

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

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

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

係数を構造式に書き込んでみると３点（右図では４点）で位相が合い， endo付加が起こることが判る．

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

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

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

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

>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

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

例 シクロペンタジエンと無水マレイン酸の４π＋２π付加反応

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

p軌道の描き方

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

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

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

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

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

シクロペンタジエンとの６π＋４π付加反応では， exo付加体のみが得られる．

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

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

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

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

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

・芳香族の置換反応

・求核置換反応

・共役化合物の安定性（芳香族性）

・周辺環状反応（付加反応，シグマトロピー反応，キレトロピー反応，電子環状反応）

・光反応と熱反応の予測

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

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

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

共役化合物の安定性（芳香族性）

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

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

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

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

構造化学（平面性）

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

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

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

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

求核置換反応（SN2反応）

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

前面からの攻撃では， 一方とは位相が合うが 他方とは位相が合わない（波線）．．

参考資料

軌道の接近方向

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

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

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

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

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

芳香族の置換反応（ニトロ化）

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

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

電子環状反応機構

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

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

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

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

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

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

配座決定

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

プログラムの利用

電子密度と結合次数

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

電子密度

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

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 ---- ２個の電子が詰まっている（電子密度計算に考慮する）    

1   0.276  0.724  0.724  0.276 ---- ２個の電子が詰まっている（電子密度計算に考慮する） 

   ↓  ↓ 

   1.0 1.0 ---- π電子密度 

結合次数

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

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 ---- ２個の電子が詰まっている（結合次数計算に考慮する） 

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

   ↓  ↓ 

  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の利用

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環境のインタプリタ，平成１１年度熊大総合情報処理センターで採用) であり，グラフィック部分の手直しが必要である．

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

＜実行例＞

＜グラフィック処理＞