Master3’s blog

LaTeXやExcelVBAなどの作例集

LaTeX作例17(3.5.3 STO-3Gを使ったHeH$^+$のSCF計算)

  • 量子化学に関する本を引用し、僕が書いたLaTeXの作例を紹介します
  • ポイントとしては、行列をarrayコマンドを用いて記述し、right, left, centerを表す{rr}{ll}{cc}絵を用いて行列要素を右揃え、左揃え、中央ぞろえに書き分けているところです。
  • プリアンブルは全部コピペして使ってるので、かなり余計なものも混ざってます。すいません
  • パッケージは基本的にデフォルトで入ってるやつが使われていると思います(たぶん)
  • ページ番号は原典と異なります
  • 『新しい量子化学―電子構造の理論入門』

    出版社 ‏ :  東京大学出版会 (1987/7/1)
  • 発売日 ‏ :  1987/7/1
  • 言語 ‏ :  日本語
  • 単行本 ‏ :  303ページ
  • ISBN-10 ‏ :  4130621114
  • ISBN-13 :  978-4130621113
  • [http://:title]

    3.5.3(later).tex - Google ドライブ

  • \documentclass{jsarticle}

    \usepackage{mathrsfs}

    \usepackage[dvipdfmx]{graphicx}

    \usepackage{parskip}

    \usepackage{indentfirst}

    \usepackage{amsmath,amssymb}

    \usepackage{braket}

    \usepackage{otf}

     

    \usepackage{calligra}

    \usepackage{calrsfs}

    \usepackage{mathrsfs}

     

    \usepackage{bm}

    \usepackage{okumacro}

    \begin{document}

     

    \subsection*{3.5.3 STO-3Gを使ったHeH$^+$のSCF計算(後編)}

     

    \parindent=1zw

     

    \hrulefill

     

    前に述べた他の2つの規格直交化の手続では、重なり行列を対角化することが必要である。2$\times$2行列の対角化は、1.1.6節で述べた。2次元の重なり行列に対しては、固有値は簡単に求まり、$s_1=1+S_{12}=1.4508$と$s_2=1-S_{12}=0.5492$である。対角化を行うユニタリー行列は

    $$\bm{U}=\left(

    \begin{array}{rr}

    [2]^{-1/2}&[2]^{-1/2}\\

    \;[2]^{-1/2}&-[2]^{-1/2}

    \end{array}

    \right)\eqno(3.259)$$

    となる。対称直交変換および正準直交変換を導くには、行列

    $$\bm{s}^{-1/2}=\left(

    \begin{array}{cc}

    s_1^{-1/2}&0\\

    0&s_2^{-1/2}

    \end{array}

    \right)=\left(

    \begin{array}{ll}

    0.8302&0.0\\

    0.0&1.3493

    \end{array}

    \right)\eqno(3.260)$$

    を使う。対称直交化では変換行列

    $$\bm{X}_{\text{symmetric}}=\bm{S}^{-1/2}=\bm{Us}^{-1/2}\bm{U}^\dagger=\left(

    \begin{array}{rr}

    1.0898&-0.2596\\

    -0.2596&1.0898

    \end{array}

    \right)\eqno(3.261)$$

    を用いる。一方、正準直交化では変換行列

    $$\bm{X}_{\text{cononical}}=\bm{Us}^{-1/2}=\left(

    \begin{array}{rr}

    0.5871&0.9541\\

    0.5871&-0.9541

    \end{array}

    \right)\eqno(3.262)$$

    を用いる。これら3つの直交化の間の関係は図3.7に示されている。もとの基底関数同士の間の角度$\theta$は、$\cos{\theta}=S_{12}$によって与えられる。Schmidtの直交化は1番目の基底関数をそのままにしておいて、それに直交する2番目の関数をつくる。対称直交化は、ベクトルの間の角度が90°になるまで対称的に開くことによって、もとの基底関数にできるだけ似ているような2個の新しい関数をつくる。正準直交化は、もとのベクトルの間を2等分するような1つのベクトルと、それに直交する2番目のベクトルをつくり出す。正準直交化を使うことにすると、変換された基底関数は

    $$\phi_1'=0.5871\phi_1+0.5871\phi_2\eqno(3.263)$$

    $$\phi_2'=0.9541\phi_1-0.9541\phi_2\eqno(3.264)$$

    となる。

     

    これで、SCFの反復手続をはじめる準備が整った。まず、密度行列の初期値の設定を行わなければならない。その際$\bm{G}$として要素がすべてゼロである行列を使うと便利である。つまり、すべての電子間相互作用を無視して($\bm{G}$をゼロ行列に置いて)得られる核$-$1電子ハミルトニアンをFock行列の最初の初期設定として使う。

    $$\bm{F}\simeq\bm{H}^{\rm core}=\left(

    \begin{array}{cc}

    -2.6527&-1.3472\\

    -1.3472&-1.7318

    \end{array}

    \right)\eqno(3.265)$$

    これが最も容易な初期設定の方法ではあるが、もっと複雑な状況においては出発点として良くない場合もある。計算のつぎの段階は、Fock行列を正準直交化された基底関数を使った表現に変換することである。

    $$\bm{F}'=\bm{X}^\dagger\bm{FX}=\left(

    \begin{array}{cc}

    -2.4397&-0.5158\\

    -0.5158&-1.5387

    \end{array}

    \right)\eqno(3.266)$$

    この行列を対角化する、つまり

    $$\bm{F}'\bm{C}'=\bm{C}'\bm{\varepsilon}\eqno(3.267)$$

    を解くと、係数のユニタリー行列

    $$\bm{C}'=\left(

    \begin{array}{rr}

    0.9104&0.4136\\

    0.4136&-0.9104

    \end{array}

    \right)\eqno(3.268)$$

    と2個の固有値

    $$\bm{\varepsilon}=\left(

    \begin{array}{ll}

    -2.6741&\;\;\;0.0\\

    \;\;\;0.0&-1.3043

    \end{array}

    \right)\eqno(3.269)$$

    が得られる。したがって、もとの基底関数についての係数はつぎのようになる。

    $$\bm{C}=\bm{XC}'=\left(

    \begin{array}{rr}

    0.9291&-0.6259\\

    0.1398&1.1115

    \end{array}

    \right)\eqno(3.270)$$

    式(3.270)と式(3.269)は、HeH$^{++}$の軌道と軌道エネルギーを与えていて、これをHeH$^+$の軌道と軌道エネルギーに対する初期設定として用いている。エネルギーの低い方の分子軌道$\psi_1$は、主として$\phi_1$からなり(係数は0.9291)、ほんのわずかしか$\phi_2$が混じっていない(係数は0.1398)ことに注意しよう。これは電子間反発を考えていないので、電子はより大きな核電荷をもったHe原子核の近くに集中する傾向があるということである。反復計算を開始して電子間反発を加えていくと、その効果は電子のHe原子核のまわりへの集中の程度を緩和して、電子をH原子核の方にやや“拡散”させる。この結果、電子間反発が減少する。

     

    式(3.270)から、つぎの反復段階に用いるHeH$^+$の密度行列に対する設定は

    $$\bm{P}=\left(

    \begin{array}{cc}

    1.7266&0.2599\\

    0.2599&0.0391

    \end{array}

    \right)\eqno(3.271)$$

    となる。$\bm{P}$の対角要素を見ると、電子密度のほとんどが、Heの近くに存在することが(あくまで定性的にだが)わかる。電子密度解析を行えば、このことがもっとはっきりとわかることであろう。式(3.271)の密度行列は、HeH$^{++}$のものではなくて(HeH$^{++}$のものとは2倍違っている)、核電荷の場の中の2個の相互作用をしていない電子の密度行列になっている。

     

    $\bm{P}$から、$\bm{G}$の初期値

    $$\bm{G}=\left(

    \begin{array}{cc}

    1.2623&0.3740\\

    0.3740&0.9890

    \end{array}

    \right)\eqno(3.272)$$

    をつくり、新しいFock行列

    $$\bm{F}=\bm{H}^{\rm core}+\bm{G}=\left(

    \begin{array}{cc}

    -1.3904&-0.9732\\

    -0.9732&-0.7429

    \end{array}

    \right)\eqno(3.273)$$

    が得られる。式(3.277)の行列$\bm{G}$の要素が正の値であることからわかるように電子間反発は正であるので、新しいFock行列の要素の値はもとの核-1電子ハミルトニアンによる設定値(式(3.265))に比べて、相当に絶対値が小さくなっている。さて、こうして得られたFock行列の固有値問題を解き、$\bm{C}$と$\bm{P}$の新しい設定値を得て、設定値とそれから得られる結果が一致するまで手続全体を繰り返す。付録Bには、この手続を行うプログラム、およびいま扱っている最小基底関数STO-3GによるHeH$^+$の計算出力が出ている。ここで行った説明に従って、その出力をたどってみる事を奨める。

     

    表3.5には、密度行列の要素と対応する全電子エネルギーが、反復回数ごとに与えられている。反復計算が進むにつれて、Hのまわりに電子が増加しHeのまわりの電子が減少していることがわかる。各反復段階における変分原理に従うエネルギーの値を得るためには、公式

    $$E_0=\frac{1}{2}\sum_\mu\sum_\nu P_{\nu\mu}(H_{\mu\nu}^{\rm core}+F_{\mu\nu})\eqno(3.274)$$

    において、$\bm{F}$をつくるのに用いたのと同じ密度行列$\bm{P}$を必ず使って計算しなければならない。したがって、エネルギーは新しい$\bm{F}$ができた直後に計算され、新しい$\bm{P}$ができた後に計算されるのではない。表3.5を見ると、エネルギーは上から単調に収束していることがわかる。エネルギーは変分原理に従う量なので、エネルギーの相対的誤差は、波動関数や密度行列における相対的誤差に比べて小さい。

     

    最終的に得られた波動関数と軌道エネルギーは

    \begin{flalign*}

    &&\bm{C}&=\left(

    \begin{array}{rr}

    0.8019&-0.7823\\

    0.3368&1.0684

    \end{array}

    \right)&\text(3.275)\\\\

    &&\bm{\varepsilon}&=\left(

    \begin{array}{ll}

    -1.5975&\;\;\;0.0\\

    \;\;\;0.0&-0.0617

    \end{array}

    \right)&\text(3.276)

    \end{flalign*}

    である。

     

    エネルギーの低い方の軌道、すなわち占有軌道$\psi_1$は、2つの係数の符号が同じであることから明らかなように結合性軌道である。これは、反復計算後もなお主としてHeの関数$\phi_1$からなっている。仮想軌道$\psi_2$は、符号の異なる2つの係数をもつ反結合性軌道である。こちらには、$\psi_1$と直交するためにHの関数$\psi_2$がより大きな重みをもっている。さて、Koopmansの定理を使えば、イオン化ポテンシャルと電子親和力を予測することができる。予測されるイオン化ポテンシャルは、正電荷を帯びた分子種の場合はいつもそうだが、大きな値をもつ(1.5975 a.u.=43.5 eV)。電子親和力は正(0.0617 a.u.=1.7 eV)であるので、HeH$^+$は1個の電子をとらえて安定化すると予測される。このことはHeHが安定な分子であるということとは違う。実際、HeH$^+$の解離した系(He + H$^+$)は大きな電子親和力をもち、これとHeH$^+$の解離エネルギーの和は、HeH$^+$の電子親和力よりずっと大きい。つまりHeHは不安定な分子である。

     

    Mullikenの電子密度解析は、$\bm{PS}$の対角要素から得られる。これによると、1.53個の電子が$\phi_1$に、残りの0.47個の電子が$\phi_2$に属することになる。正味の電荷でいえば+0.47がHeに+0.53がHに属していることになる。したがって、形式電荷+1は、2個の原子にほぼ均等に分かれていると予測される。L\"{o}wdinの電子密度解析はプライムのついた密度行列から得られる。このプライムのついた密度行列は、規格直交な基底関数$\{\phi_\mu'\}$についての密度行列である。Hの軌道に属する電子の数は

    $$P_{22}'=(\bm{S}^{1/2}\bm{PS}^{1/2})_{22}=2(\bm{S}^{1/2}\bm{C})_{12}^2=0.5273\eqno(3.277)$$

    である。これは、HeとH上におのおの+0.53と+0.47の正味の電荷があることを意味し、Mullikenの場合とほぼ同様の電荷分布が予測される。しかし、両者ではHeとHで正味の電荷の大小関係が逆になっていることは注意に値する。

     

    HeH$^+$の全エネルギーは、電子エネルギーに核間反発2/Rを加えることによって得られ、$-2.860662$ a.u.である。すでに計算した積分を使ってH、He$^+$およびHeのエネルギーも計算することができる。この基底関数によるH原子のエネルギーは、H$_2$の計算の場合と同じで、$T_{22}+V^2_{22}=-0.4666$ a.u.である。また同じ基底関数による1電子原子He$^+$のエネルギーは同様に$T_{11}+V^1_{11}=-1.9755$ a.u.である。He原子は$\phi_1$に2個の電子があり、そのエネルギーは運動エネルギーとHe原子核からの核引力エネルギーに加えて2つの電子の電子間反発からの寄与を含んでいる。したがって、He原子のエネルギーは$2(T_{11}+V^1_{11})+(\phi_1\phi_1|\phi_1\phi_1)$である。これらのエネルギーを表の正確な値と比較するために表3.6に示しておいた。

     

    表3.6の結果から、2つの解離過程に対する解離エネルギーが計算される。

    \begin{flalign*}

    &&{\rm HeH^+\to He+H^+}\hspace{10mm}&\Delta E=0.2168\;{\rm a.u.}&\text(3.278)\\

    &&{\rm HeH^+\to He^++H}\hspace{10mm}&\Delta E=0.4168\;{\rm a.u.}&\text(3.279)

    \end{flalign*}

    この結果、HeH$^+$は開殻系He$^++$Hに解離するのではなく、閉殻系He$+$H$^+$に解離することを正しく予測できた。しかし解離エネルギー0.2168 a.u.=5.90 eVは、正しい値2.04 eVに比べてかなり大きい。この原因は主として、私たちが用いたHeの軌道指数が2.0925で、これはHeH$^+$分子に対しては適当だが、解離生成物であるHe原子に対する最良の値1.6875に比べて大きすぎるためである。つまりHeのエネルギーの計算値が、HeH$^+$のエネルギーに比べて相対的に高すぎてしまったのである。

     

    図3.8には、私たちの用いた軌道指数による全体のポテンシャル曲線とWolniewiczの実質的に正確な結果が示されている。ここで計算されたSTO-3Gによる平衡結合長は1.3782 a.u.である。この値は、ポテンシャル曲線のくぼみの深さがひどく大きすぎるにもかかわらず、正確な結果に大体一致している。解離生成物が閉殻系であるので、H$_2$のときに直面したような困難(図3.5を参照)はなく、計算による解離過程の振舞は正しい。H$_2$に対して行ったように、この解離を解析的に考えてみよう。結合が伸びていくと、$\psi_1$における$\phi_1$の係数は大きくなり$\phi_2$の係数は小さくなっていく。したがって電子はますますHe原子核のまわりに集中することになる。極限では、$\psi_1$は$\phi_1$になってしまう。つまり

    $$\bm{C}_{R\to\infty}=\left(

    \begin{array}{cc}

    1.0&0.0\\

    0.0&1.0

    \end{array}

    \right)\eqno(3.280)$$

    である。対応する密度行列は

    $$\bm{P}_{R\to\infty}=\left(

    \begin{array}{cc}

    2.0&0.0\\

    0.0&2.0

    \end{array}

    \right)\eqno(3.281)$$

    となる。$R\to\infty$のとき全電子エネルギーは、すべての2中心積分をゼロに置くことによって計算できる。このときに残る積分は、$T_{11},\;T_{22},\;V^1_{11},\;V^2_{22},\;(\phi_1\phi_1|\phi_1\phi_1)$および$(\phi_2\phi_2|\phi_2\phi_2)$だけである。

     

    \hrulefill

     

    ここで述べたHeH$^+$の計算では、軌道指数に標準的な値を用いた。より良い計算は、各核間距離において最適化した軌道指数を用いて行うことができる。この手続を行ったとすると、Heの軌道指数$\zeta_1$は$R$が大きくなるにつれて原子に対する最良の値1.6875に向かって減少していく。Hの軌道指数の値は、$R$が大きいところではゼロへと減少していくであろう(解離生成物においては2個の電子はともにHeの上にあるので、$R$が大きいときHの基底関数が軌道に寄与できる唯一の可能性は、非常に広がってHe原子核の領域まで伸びていくことである)。この問題では2個の軌道指数しかないので、軌道指数を最適化するのは比較的容易である。たとえば、軌道指数を1個ずつ順々に最適化することを繰り返していくという方法がある。軌道指数がたくさんある一般的な問題では、極小点を見いだすためには非常に多くの局所的な極小点を有する多次元空間を探査しなければならない。こうした多くの軌道指数を扱う問題においては、ふつう軌道指数の最適化は行わない。

     

    H$_2$とHeH$^+$のモデル系について制限つき閉殻Hartree-Fock計算を説明したので、これから多原子分子のより実際的な計算結果を例示していこう。この例示の目的は、現在行われている計算を概説しようというのではなく、制限つき閉殻Hartree-Fock型のすべての計算の背景となっている考え方の要点を示し、そうした計算が実験値に比べてどのくらい良いか(あるいは、どのくらい悪いか)についての感じと直感的判断力を読者につかんでもらうことである。私たちの計算で扱う分子は、H$_2$、N$_2$およびCO、そして一連の10電子分子、CH$_4$, NH$_3$, H$_2$O, FHとする。これらの分子のおのおのに対して、いくつかの明確に定義された基底関数系を用いる。しかしながら、それらの計算結果を記述する前に、多原子系の基底関数の全般に関する問題を議論し、ここでの計算で使う特定の基底関数を説明しておく必要がある。

     

    \end{document}