[計算Tip] | 簡略化χSTO | ガウス型χGTO | 基底関数χCGTO | 積分法 | 固有値 解法 | 拡張Huckel | 収斂判定 | 物理量 | 軌道法 拡張 | [描画Tip] | 軌道表面 | 立体表現 | 量子化学計算勉強室ホーム |
CD付きで本になりました。 「はじめての量子化学計算<基礎と可視化>」(工学社) |
Schroedingerの方程式の導出・方程式解など基本的なことは、「国立科学博物館で学ぶ物理学 : シュレーディンガー方程式/原子の電子雲」を参照してください。
: 分子コア積分
: 交換積分(多中心積分)
: 密度行列
: 重なり積分
(4) HF-SCF固有値問題の解法 (永年方程式の行列表記 : F=(Frs)、C=(Cir)、S=(Srs)、E=diag(εi))
<距 離> 1au = 0.5292 Å (a0=ε0h2/(πme2)=0.5292 Å) <エネルギー> 1au = 627.17 Kcal/mol
原子軌道は、水素類似原子の解であるスレータ型原子軌道を出発点とするが、距離raの低次数部分は
主量子数nの下位の軌道との組み合わせにより表され、固有値問題を解く時の重なり積分Sの変換(基底関数の直交化
操作)でその不都合が解消されるので、raの最高次数部分で代表する。
即ち、3s軌道を例にとれば
規格化条件の積分を、実距離raの座標での積分ではなく、無次元距離 r の座標での積分で行うと
[U] 分子軌道法HF-SCFを解くためのTip集
(Tip-1) 簡略化スレータ型原子軌道χSTO
ではなく、以下の式とする。但し、係数は規格化条件 ∫χ2dv=1 を満足するように決める。
詳細は「簡略化スレータ型原子軌道の具体的表現」(pdf)参照
(Tip-2) ガウス型軌道χGTO
スレータ型軌道は多中心積分が不可能なため以下に示すガウス型軌道を採用する。 係数Nは、規格化条件 ∫χ2dv=1 より決まる。
但し r=ra/a0
但し (ex) 9!! = 9・7・5・3・1
詳細は「ガウス型軌道の具体的表現」(pdf)参照
(Tip-3) 基底関数系
スレータ型原子軌道を、以下に示すようにガウス型軌道の線形結合で近似(展開)する。 ガウス型軌道を用いる典型的基底関数 STO-nG (n=3の時、STO-3Gという)なら
これらの短縮(展開)係数di、短縮軌道指数αi は、規格化条件 ∫χCGTOχCGTOdv=1 を満足しなければならない。すなわちfor s 軌道
for p軌道
for d軌道
これらdi、αi は、この条件化で以下の式の値を最小とするものとして求められる。但し、χSTOは簡略型を使用する。
2s,2pに対してはdiは異なるがαi は同一(手計算時代の手計算負荷軽減の名残か?)として求める。αi は、スケール因子ζ=1に対応するものとして得られる。
詳細は「基底関数系の短縮係数・短縮軌道指数の具体的数値」(pdf)参照
手っ取り早く全データがほしい方は「EMSL Baias Set Library」から入手
(Tip-4) 積分Frs,Srs,Hrs
HF-SCFの計算手順は簡単であるが、最も困難なのがこの積分であり、この方法が殆どの教科書には書かれていない。
任意の演算子をH(1/2▽2、ZA/rAμ、1/rμν)とすると
となるので、重なり積分、分子コア積分(2つに分けて、運動エネルギー関連Hrs1、クーロン関連積分Hrs2 )および多中心積分を、 s軌道同士の積分を例にとり、以下のとおり定義する。 ここに現れる(SaSb)、[SaSb]、<SaSb>、(SaSb|ScSd)は、次のように解析的に解ける。但し、F0(x)=erf(x)/x である。 (誤差関数erf(x)=∫xexp(-u2)du ): p=α+β、s=αβ/(α+β)、Rab=(a〜b)距離
: Rs=(p〜q)距離
: q=γ+δ、t=γδ/(γ+δ)、Rcd=(c〜d)距離、Rpq=(p〜q)距離
s軌道以外の軌道との積分(ex(Sa Xb))、s軌道以外同士の積分(ex (XaYb))は、これらの微分で求めることができるが、膨大な式となるため すべてを書き下してプログラム化することは不可能である。そこで、プログラム上では、微分を1回するたびに微分形 の係数や指数等々を記憶させて、最終微分形を作ることでこの問題を解決することができる。
他の積分も同様であるが、特に多中心積分(rs|tu)は、以下の性質 (rs|tu)=(sr|tu)=(rs|ut)=(tu|rs) etc を有するので、プログラム作成においてはそれらの対称性を巧く利用して計算回数を減少させる工夫が肝要である。
計算法概要詳細は「重なり積分・分子コア積分・交換積分の計算法概要」(pdf)参照
微分ルールの体系的詳細は「重なり積分・分子コア積分・交換積分の微分法整理」(pdf)参照
具体的な式詳細例は「重なり積分・分子コア積分・交換積分の具体的表現」(pdf)参照
(Tip-5) 固有値問題の解法
一般固有値問題FC=SCE(Sが単位行列ではない)を、通常の固有値問題FC=CEの形に変形してから問題を解き、分子軌道の係数Cir(固有ベクトル) および各軌道のエネルギーεi(固有値)を計算する。
Sの固有値Λ=diag(λi)、固有ベクトルX=(xij)は、SX=XΛの解であるが、 X'=XΛ-1/2 ここに Λ-1/2=diag(λi-1/2) F'=tX'FX' C'=X'-1C とおくと、FC=SCE が F'C'=C'Eとなり、EとC'が得られる。求める答えは C=X'C'
通常の固有値問題SX=XΛやFC'=C'Eの解法には種々の方法があるが、収斂速度は速くないが固有値と固有ベクトルが同時に求められ且つ論理が簡単 であることから、Jacob法を適用する。収斂判定(詳細はpdfを参照)は単純に言えば一般的に10^(-6)とするが、場合によっては収斂しないこともあり、 そのときは一時的に10^(-5)とした後10^(-6)で収斂することもある。
固有値問題の解法の詳細は「固有値問題の解法(Jacob法)」(pdf)参照
(Tip-6) 分子軌道の係数Cirの初期値(拡張Huckel法)
分子軌道の係数Cirの初期値はすべて0から出発してもかなりの問題は計算可能であるが、時々おかしな答えとなる時も あるので、適切な初期値を使うようにする事が必要である。この初期値としては、拡張Huckel法の解がよく使用される。
拡張Huckel法では、原子価電子のみのスレータ型軌道を用いて解こうとするものであるが、スレータ型軌道による交換 積分(rs|tu )の計算が困難なため、以下の近似のもとに計算する。
即ち、Frsの対角要素は原子価状態イオン化ポテンシャルIpを用い、非対角要素は経験的なWolfberg-Helholzの近似式を用い, Frr=-Ip Frs=(1/2)KSrs(Frr+Fss) 但し 一般的に K=1.75 こうすると、1回の固有値計算でFC=SCEの解C=(Cir)を得ることができる。ここでも、重なり積分の計算が結構複雑であり、またイオン化ポテンシャルIpや、スレータ則(内殻省略)によるスレータ型 軌道の原子の電荷Z'、主量子数n'の取り扱いが簡単とはいえない。
拡張Huckel法での重なり積分の概要は「拡張Huckel法とその計算法」(pdf)参照
Ip,Z',n'の詳細は「イオン化ポテンシャル・スレータ則」(pdf)参照
拡張Huckel法での重なり積分の詳細(4p軌道までの計算)は「拡張Huckel法とその計算法の補足」(pdf)参照
拡張Huckel法での重なり積分の具体的表現(4s軌道までの計算)は「拡張Huckel法における重なり積分具体的表現」(pdf)参照
(Tip-7) HF-SCF収斂判定基準
計算フローチャートではCir=Cir'と記したが、実際的には、密度行列 P = (Ptu) を用いて判定する。
(1) Ptuの計算 Tip-5で求めた分子軌道のエネルギーεiの小さい順に並び替えて、一番小さいεiの分子軌道を1番、 その次を2番と番号も入れ替え、i=1〜N(分子軌道数)まで加算(一般的に「被占有軌道すべてについて 加算する」という。)してPtuを求める。初期値も0行列即ち P = (0)でもよいが、Tip-6で述べた密度行列 P(Ptu)を使用するのが好ましい。
(2) 収斂判定 上のようにして得られたPtuと前回の計算で得られたP'tuを用い、通常以下の判定基準とする。(後述する振動解析のための数値微分のときは(-5)乗くらい必要)
(Tip-8) その他の物理量
(1)全電子エネルギー/分子エネルギー
全電子エネルギーは、拡張Huckel法の場合は各軌道のエネルギーεiの単純和でよいが、HF-SCF法では以下のとおり。: 全電子エネルギー: 核反発エネルギー: 分子エネルギー
(2)電荷(Mullkikenのpopulation analysis)
(3)双極子モーメント
詳細は「分子エネルギー、電荷、双極子モーメント」(pdf)参照
(Tip-9) 分子軌道法の拡張
(1) 局在化分子軌道
(2) 分子構造解析
(3) 基準振動解析
(4) 赤外線吸収強度/IR解析
詳細は「分子軌道法の拡張」(pdf)参照
核座標微分は「原子核座標での微分」(pdf)参照
(Tip-10) 分子軌道表面の決定
分子軌道は無限に広がったものとして計算しているため、理屈では全空間に電子が存在し表面がない。そこで、各分子軌道 の電子の確率密度(実数表現の波動関数の2乗)を数値積分して、たとえば95%の電子存在確率を与える最外層の波動関数 一定値の面を分子軌道の面とする。我々の計算ではおおむね以下の値となる。 電子存在確率密度=0.001 波動関数=(0.001)の平方根 一般的には、波動関数の値を利用者が選択するようにしたものが多いようである。
(Tip-11) 分子軌道表面の立体表現
どんなコンピュータグラフィクスの本にも書かれているので不必要かもしれないが、老婆心ながらポイントを紹介しておこう。
<反射強度計算> I=I0{cos2(θ)+0.5cos(i)} i:入射角、θ:視線と反射光の角 これを入射光RGB(I0)それぞれに適用し(白色光ならすべてのRGBは同一強度)、 反射面の本来の色RGBに加える。 <軌道表面探索と立体化> @ 見たい角度に図を回転する座標変換を決定する。 A 見たい角度に回転した座標での(x,y)面の各点で、 視点側と反対側の、想定軌道より外側のz座標を決める。 B それら2点の座標を元の座標に戻し、戻した2点間で、 たとえば指定の電子存在確率密度を与えるzの値(視点に近い側の値)を探索する。 C こうして得た元の座標での点(x,y,z)での法線を計算する。 D その点ならびに法線を再び座標変換し、座標変換した法線を基準に反射光の方向を求める。 E 反射光の方向ならびに視線との角度を求めて、上に述べる反射強度計算をする。
使用した書籍を以下に紹介しておくので、説明不足又は意味不明などはこれらの書籍をみてほしい。
(1) (三訂)量子化学入門(上・下) 米澤他 化学同人 1983 第3版 : 入門から応用まで全般 (2) 電子構造の理論入門 新しい量子化学(上・下) ザボ他 東京大学出版会 2003 6版 : 特に多中心積分が参考 (3) 分子軌道法 藤永 岩波書店 1988 第2刷 : 基底関数 (4) 分子科学講座2 分子理論と分子計算 樋口編 共立出版 1986 1版 : 紹介していない開殻制限HF法(ROHF)の研究用 (5) 実験化学講座(3) 基本操作V 日本化学会 丸善 1991 4版 : 入門全般 (6) 実験化学講座(12) 計算化学 日本化学会 丸善 2003 5版 : 振動解析法他、新解析法
以上を読破してもp軌道以上を含むプログラムが作成できるところまでは書かれてはいないのが現実です。