《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

第一原理電子状態計算の新規手法に関する一提案：
記号的-数値的計算法
A proposal to first principles electronic structure calculation:
Symbolic-Numeric method
菊地章仁
Akihito Kikuchi
CANON INC.
30-2, Shimomaruko 3-chome, Ohta-Ku, Tokyo 146-8501-Japan
E-mail:kikuchi.akihito@canon.co.jp
要約
本研究では，記号的―数値的計算手法を利用した第一原理電子状態計算を試みる。ハートリ
ー・フォック・ロータン方程式を用いた分子軌道計算において，基礎方程式を連立多項式方程式系と見なして近似を行った後，数式処理によって同一の零点を持つ連立多項式方程式系を作成し，その解を数値的に求めることによって自己無撞着電子状態計算を実行する。
こうすることで，電子状態計算.・最適化計算・逆問題計算を，統一的に，すべて順問題として解くことができる。その際，試行錯誤的な反復計算は不要である。
Abstract
This study proposes an approach toward the first principles electronic structure
calculation with the aid of symbolic-numeric solving.

The symbolic computation

enables us to express the Hartree-Fock-Roothaan equation in an analytic form andapproximate it as a set of polynomial equations. By use of the Gröbner basis technique,the polynomial equations are transformed into other ones which have identical roots.
The converted equations take more convenient forms which will simplify the numerical
procedure, from which we can derive necessary physical properties in order, in an “a lacarte” way. This method enables us to solve the electronic structure calculation, theoptimization of any kind, or the inverse problem as a forward problem in a unified way,in which there is no need for iterative self-consistent procedures with trials and errors.
1.背景
第一原理的電子状態計算の基礎方程式は，電子－原子核系のエネルギー汎関数の極小条件に基づいて導出されるハートリー・フォック方程式あるいはコーン・シャム方程式である。これらの方程式はシュレーディンガー方程式の型に表される[1-3]。

⎛ 1
⎞
Za
ρ (r ' )⎜⎜ − Δ + ∑+ ∫ dr '+ V exc ⎟⎟φ i (r ) = Eiφ i (r )| r − r '|a | r − Ra |
⎝ 2
⎠

(1)

左辺括弧内の第二項は，電荷 Za をもつ原子核のポテンシャルである。第三項は，電荷分
布 ρ が作るクーロンポテンシャルである。第四項は，多電子系に働く量子力学的な相互作

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

用であり，交換相関ポテンシャルと呼ばれる。この方程式は，波動関数を何らかの基底関数系によって展開することによって数値的に解かれる。波動関数を原子局在基底の一次結合として表すと，次のような行列の固有値問題が得られる（Hartree-Fock-Roothaan 方程式，以下 HFR 方程式と略記）。

)
H ({R}, {α }, {Q}, {C})C ({R}, {α }, {Q}) =Sˆ ({R}, {α }, {Q})C ({R}, {α }, {Q})e({R}, {α }, {Q})(2)

ここで， Ĥ はハミルトニアン行列， Ŝ は重なり積分行列， C は波動関数（原子局在基底の一次結合係数）， e は固有値である。また，{R}は原子核の位置，{α}は基底関数の空間的な広がりを指定する係数(軌道指数)，{Q}は量子数であり，原子基底はこれらを用いて表現される。
HFR 方程式を解析式によって直接表現すると，これは， {R}，{α}，{Q}を変数にもつ多変数関数であり，超越関数を含む複雑な分子軌道方程式になる。これは一電子積分，二電子積分を構成する際の基底関数として，通常はガウシアン型(Gaussian Type Orbital, GTO)
またはスレーター型軌道(Slater Type orbital, STO)を用いるためである。 解析式型基底は，数値計算上は効率的であるが，解析表現が複雑なため数学的な解析は困難を伴う。しかしながら，この方程式に含まれる変数について多項式近似を行い，簡易な表式を得ることが可能である(ワイエルシュトラスの近似定理)。この目的では，テイラー展開を行えば十分である。分子軌道方程式は原子座標，軌道指数，量子数を用いて表される高次代数方程式の連立系となる。その具体的な表現は数式処理を用いることで計算可能である[4-11]。なお分子軌道法における HFR 方程式の多項式近似という概念は，安井によって開発されたものである[6-9]。安井は，多項式近似した HFR 方程式を「分子軌道代数方程式」と名付けている。
このことによって，分子積分の多項式近似の上に分子軌道の多変数理論を構築することができる。分子の性質は多変数関数として表現され，通常各変数間の連関を明確に追及することが可能になる。
通常の方法の電子状態計算は材料物性計算における順問題である。すなわち，物質の構造データの仮定値を入力し，電子状態計算を行い，全エネルギーが極小になるように構造を最適化し，最安定構造のもとで，種々の電子物性を求める。従来法には，次の課題がある。第一原理分子動力学法において，原子核は古典的なニュートン方程式に従うものであり，一方，波動関数は量子力学に従うと見なされる。その基礎となる理論が，いわゆる断熱近似である。断熱近似を適用することにより，波動関数と原子核とは，別々のモデルで取り扱われる。したがって，従来法では，緩和法による波動関数最適化と，原子配置の最適化という二つの計算フェーズを交互に繰り返しながら計算を進める。この方法は，数値的に安定である，とされる。しかし，波動関数と原子核配置の最適値を同時に求めるという見地からすると，回り道の多い，緩和問題の解法として非効率な手段になることがある。
例えば，原子核と波動関数の自由度を分離して扱っているため，原子核と波動関数のダイ
ナミクスが強く結合するような物理現象は，従来法では扱いにくいのである。これに対して，逆問題は，望ましい電子物性を与えるような物質構造を探求するということになる。
従来法を用いて逆問題を検討する際に，まず物質の構造を仮定し，これに対する電子物性を求め，計算結果が望ましい値に近づくように仮定を少しずつ変化させながら繰り返し計《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

算を行う。すなわち，仮定した構造を少しずつ変え，順問題を繰り返しとくことによって，試行錯誤的に逆問題の解へいたる，という道筋しかもたない。従来法の電子状態計算において，固有値計算・自己無撞着場計算・原子核系の構造緩和計算の順という一定の順序にしたがって，複数の緩和計算が入れ子構造になっている。未知変数は，いかなる場合でも，内側から外側のループまで，決まった順番にしたがって決定しなければならない。それが，従来法における逆問題の解法が，順問題の試行錯誤的な繰り返しになる理由である。
2. 記号的-数値的第一原理分子動力学・分子軌道法
われわれは，この状況を鑑み，次の手法を提案する。
「ハートリー・フォック・ロータン方程式を連立多変数多項式系として近似表現し，さらに，この連立多項式系を数式処理によって『より解きやすい形』に変更し，その根を記号
的-数値的に求めることによって固有状態を求める。」
本研究における課題は，分子軌道方程式を多項式近似した後に，どのように方程式の解を求め，有意義な情報を導くかということである。一般に，有限個数の解を持つ連立多項式系を解こうとする場合には，浮動小数点近似に基づき解は数値的に計算される。しかしながら，数値誤差等の影響によって，純粋に数値計算的な手法はしばしば不安定になり，
その不安定性が何処に向かうのか予測が困難になることが多い。そこで，記号計算(数式処理)と数値計算を組み合わせた解析方法(Symbolic-numeric solving)が提唱されている。
Symbolic-numeric solving では，記号計算によって，前処理として，記号計算を連立多項式系に適用し，同一の根を持ち，かつ安定かつ容易に数値計算が可能な連立方程式を作成する。さらに，この時点で，解の特徴，例えば，解が存在するか，どういう幾何学的構造を持つかなどを判断する。その後，数値計算に移行するのである。数学的背景は[12-15]を参照のこと。計算化学分野における記号計算の適用事例は，[16]にレビューがある。
本手法は，HFR 方程式を解くために，Symbolic-numeric solving の幾つかの手法を利用する。連立多変数多項式系として近似した HFR 方程式は，同一の根を持つ方程式系に変換される。一つの方法として，三角化のアルゴリズム(参考文献[12,13,14,15])を適用する。このアルゴリズムにおいて次のような変換が行われる。
出発点となる方程式系(The starting equations) , … ,f 1 ( x1 , x 2 ,..., x n ) = 0(3)f 2 ( x1 , x 2 ,..., x n ) = 0
M
f n ( x1 , x 2 ,..., x n ) = 0→,…,

order) { g i }の辞書式順序グレブナー基底(Gröbner basis with lexicographic monomial《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）(4)

g1 ( x1 ) = 0
M
g 2 _ 1 ( x1 , x 2 ) = 0
M
g 2 _ m ( 2 ) ( x1 , x 2 ) = 0
g 3 _ 1 ( x1 , x 2 , x3 ) = 0
M
g n _ 1 ( x1 ,K, x n ) = 0
M
g n _ m ( n ) ( x1 ,K, x n ) = 0→ 三角化された多項式系(Triangular sets of polynomials){ }t1 ( x1 ) = 0(5)t 2 ( x1 , x 2 ) = 0
Mt n ( x1 , x 2 ,..., x n ) = 0
文献[12,13]に示されたアルゴリズムでは，出発点の方程式系 , … ,からまず辞書式順序のグレブナー基底{ g i }をつくる。このグレブナー基底も方程式系 , … ,と同等の根を持ち，数値的に「より解きやすい」形式になっている。すなわち未知変数の少ない多項式からは
じまり，未知変数の多い多項式にいたるように，多項式系が整理されている。しかし，多項式の数が出発点の方程式系よりも一般に大きくなる。この段階で，根の探索を実施する
こともできる。ここで，グレブナー基底に対し，三角化のアルゴリズムを適用すると，(5)に示すように，それぞれ未知変数x ，未知変数x , x ，…，未知変数x , x , … , x をもつ n 個の多項式から成る三角多項式系t , t , … , t を得る。なお，出発点の方程式系 , … ,複数の根を有する。そのため，方程式系 , … ,は，通常，の全零点を与えるために，三角多項式系

t , t , … , t が複数必要になることがあるが，参考文献の三角化アルゴリズムでは，必要な三角多項式系すべてを生成することが可能である。三角化された方程式系は，グレブナー基
底系よりも，方程式のサイズが小さく，解きやすいものである。こうして得られた三角化
方程式系t , t , … , t を用い，未知数x , x , … , x をひとつずつ順に求めていくことができる。
すなわち，HFR 方程式を多項式系とみなし，数式処理によって三角化することで，未知数である固有値，波動関数，各種パラメータの値を，順を追って決定することができるようになる。この際，数値計算に必要なのは，一変数に対するニュートン法あるいはその類似の方法だけである。
さらに，この種の方程式を解くために Symbolic-numeric solving の別の方法も利用できる。この解法の基礎は，Stickelberger の定理である。HFR 方程式を変数X , X , … , X で表される多項式の集合 , … ,
R

k

,…,

と見る。この多項式の集合は，数学的な見地では，可換環上の零次元イデアル I を構成し，イデアル I の零点（イデアル I を連立方程式と見たときの根）は剰余環 A＝R/I に対応する。 kは，多項式環の係数体であり，われわ
れの場合，有理数体もしくは実数体としてよい。A はk上有限次元のベクトル空間であり，このベクトル空間の基底は，X , X , … , X の単項式によって表される。したがって，A=R/I

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

において，各基底に X , X , … , Xを演算した結果も，A の基底の線形結合になる。すなわ

ち，X , X , … , X の積演算は，基底間の線形変換行列mh

X , X , … , X )として表されるのである。剰余環 A＝R/I の基底および変換行列m は，イデアル I のグレブナー基底を利用して求めることができる。Stickelberger の定理は，「行列m の固有ベクトル
I の零点ξ

ξ ,…,ξ

とイデアルが一対一の対応をする」ということを主張するものである。その対応はm ·である。ここで，固有ベクトル·(6)は，すべてのm について共通になる。（文献 14，

p101-130, ”From Enumerative Geometry to Solving Systems of Polynomial Equations”,
Frank Sottile を参照のこと。）零点 ξξ ,…,ξ

を求めるための数値計算は次のように行う。どれか一つの X に対して，上の固有値方程式を数値計算し， を求める。そうして得
られた
ξ

を他のm

ξ ,…,ξ

ji に掛け，固有値ξ ji を求めることによって，イデアル I の零点が計算される。

ここで，エネルギー汎関数，波動関数の規格化条件，HFR 方程式は波動関数（LCAO の係数）と固有値については多項式である。ここで，これらの解析式は，分子積分を用いて構成されるが，分子積分は，そのパラメータについては一般に解析関数で表され，多項式とはならない。[文献１,3-7]。そこで，分子積分は，そのパラメータに関する近似多項式に
よって置き換えるものとする。そうすることによって，波動関数と固有値ばかりでなく分子積分のパラメータまで含めた連立多変数多項式系，すなわち分子軌道代数方程式を作ることができる。ハートリー・フォック・ロータン方程式に他の拘束条件が加わる場合も，拘束条件の近似多項式を作り，方程式系に追加すればよい。方程式系の係数（実数）は有理数化しておけば，記号的に方程式系を変形する際に任意精度計算を行うことによって精度の低下を避けることができる。
この方法の一つの利点は次のようなものである。従来法における入力データは分子構造であり，出力は電子状態である。しかし，本法では，入力データは，分子構造に限定されない。HFR 方程式および付随の拘束条件に現れる任意の変数を選び入力値とすることができる。解こうとする問題が適切に設定されていれば，該入力値に対して，他の未知変数がどのような値をとるのかを計算できる。問題が適切に設定されているかどうか，すなわち
連立多変数多項式系に解があるかどうかは，代数学のイデアル理論に基づき，連立多変数
多項式系のグレブナー基底の零値を与える解集合が存在するか否かという条件によって判定する。
本手法のフローチャートを順に解説する。
１．エネルギー汎関数の解析式を求める。
変数として，波動関数（LCAO[Linear Combination of Atomic Orbitals]の係数），分子積分のパラメータ（原子座標，原子基底の減衰指数），その他の物理量を用いる。
２．拘束条件を解析式化する。汎関数および拘束条件中の係数（実数）は有理数化しておけば，記号的に方程式系を変形する際に精度の低下を避けることができる。
３．エネルギー汎関数と拘束条件は，固有値および波動関数（LCAO の係数）の多項式である。しかし，一般には，これらの式は，分子積分のパラメータの多項式にならない。
そこで，分子積分は，そのパラメータに関する近似多項式に置き換えるものとする。

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

例えば，分子積分のパラメータに関して，ある基準点（基準数値）を設定し，その基準点の周囲でテイラー展開することで，近似多項式を作ることができる。
４．多項式近似されたエネルギー汎関数と拘束条件について，エネルギー汎関数の極小条件に基づき，各変数に関する微分の解析式をつくる。これが解くべき多項式方程式系
である。
５．作業４で作られた多項式方程式系を数式処理し，同一の根を持つ連立多項式系に変換する。その際，多項式方程式系は，一旦グレブナー基底に変換し，ここで，解の存在可能性を検討する。解が孤立点の集合である場合，すなわち，解が数値として求められる場合は，グレブナー基底を三角化する。あるいは，グレブナー基底を利用し，Stickelberger の定理に基づき根の探索を固有値問題の形態に置き換える。
６．作業５で得られた方程式系を，数値的に解くことによって，この方程式系の根を求める。そうすることで，電子状態およびその他の未知変数の値が求められる。
3. 計算例
計算例を示す。以下，単位は原子単位系（atomic units）とする。まず，自己無撞着電子状態計算が可能であることを見る。対象として，水素分子を扱う。以下の解説では，最も単純な分子である水素原子の計算を示すが，本方法の適用対象は，二原子系・二電子系に限定されるものではないことに注意する。水素分子を取り上げた理由は，この系は，単純ではあるが，物質中に働く電子－電子間，電子－核間，核－核間の量子論的な相互作用を全て含むものであって，一般の多原子・多電子系の雛型と見なすことができるからである。水素分子の計算に必要な分子積分(一中心積分，二中心積分)は STO 基底を用いて作成した。エネルギー汎関数は，水素 A，B 其々の一中心積分，二中心積分，また波動関数の係数 a，b，c，d の解析式であり，核間距離 R については指数関数を含む。二中心積分の一例
を(7)に示す。
これは 1s 軌道間の二電子反発積分であり，クーロン型と呼ばれるものである。
その定義はdrdr ′

RA

′

RA
|

′|

RB

′

RB

であり，[1s(A)1s(A)||1s(B)1s(B)]と表すことができる。ここで 1s(A),1s(B)は原子 A,B に中心を持つ 1s 軌道の意である。式
中，za, zb, zc, zd は四つの 1s 軌道の軌道指数，R は核間距離，E R

はexp

za

zb R の意である。二電子反発積分として，この他に，Exchange 型[1s(A)1s(B)||1s(A)1s(B)]やhybrid 型[1s(A)1s(B)||1s(A)1s(B)]型の積分もある。一般には，分子積分は，この式よりもさらに複雑な超越関数を含む式となる。STO 基底を用いる場合は，指数関数や指数関数積分などが現れ，Exchange 型反発積分に至っては無限級数和となる。 しかし，如何なる複
雑な解析式であっても，これを有限次数の多項式で近似することで，数式処理が容易になる。しかも，STO 基底は，原子核近傍と遠距離の局在軌道の物理的な性質を GTO 軌道に
較べて精確に表現することが可能であり，原子座標に関するテイラー展開によって分子方
程式を多項式化するという目的に対して有利である。この理由から，本研究では STO 基底を用いた。ただし，以下に述べる処方は GTO 基底に対しても，また，AM1(Austin Model
1)[17]，PM3(Parameterized Model number 3)[18]，強結合モデル等の半経験的手法に対しても適用可能である。
[1s(A)1s(A)|1s(B)1s(B)]=(7)

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

本方法の順問題的第一原理分子電子動力学法に対する適用事例として，構造最適化（原子間距離の最適化）と UHF(Unrestricted Hartree-Fock)電子状態計算を同時に実行する。
ここで，水素 A，B の位置をR A およびR B とする。核間距離をrxA|x

R A | ，xB|xRA

R B とする。また，

R B |という記法を用いる。UHF 計算を行うので，上向き・下向きのスピンに対する試行波動関数を

φup ( x) = (a exp(− x A ) + b exp(− x B ) ) / π(8)φdown ( x) = (c exp(− x A ) + d exp(− xB ) ) / π

(9)

および

とする。それぞれに対応する固有値を ev，ew とする。なお，基底（軌道指数 1）を用いた
計算値は実験値と良く一致するものではない。実測を再現するためには，軌道指数を適切な値に設定する必要がある。軌道指数 1 の基底を用いるのは，単に，数式処理のコストを削減する目的である。
エネルギー汎関数は，解析式を数式処理によって作成したのち，R0=7/5(a.u.) の位置で，核間距離 r に関し四次のテイラー展開を行い，多項式化した。汎関数は，通常の分子軌道法の作法に従い[1]，電子－原子系のエネルギーに波動関数の規格直交条件(ノルム１)を課したものである。ラグランジュ乗数が固有値になる。多項式化した解析式において，実数係数は，絶対値 1/1000 以下のものを切り捨て，有理数で近似した。定数倍のような操作は行っておらず，単位は Hartree unit である。この式を(10)に示す。(10)中の記号は次の通りである。(a,b)は上向きスピン電子の LCAO 係数，(c,d)は下向きスピン電子の LCAO 係数，evは上向きスピン電子のエネルギー固有値，ew は下向きスピン電子のエネルギー固有値，rは核間距離である。式(10)の右辺各項は，すべてエネルギーの次元を持ち，したがって数値
係数は多変数単項式の次数に応じてそれぞれ異なる次元を持つことになる。
Ω＝(3571 - 1580*a^2 - 3075*a*b - 1580*b^2 - 1580*c^2 + 625*a^2*c^2 +1243*a*b*c^2 + 620*b^2*c^2 - 3075*c*d + 1243*a^2*c*d + 2506*a*b*c*d+ 1243*b^2*c*d - 1580*d^2 + 620*a^2*d^2 + 1243*a*b*d^2 +625*b^2*d^2 + 1000*ev - 1000*a^2*ev - 1986*a*b*ev - 1000*b^2*ev +1000*ew - 1000*c^2*ew - 1986*c*d*ew - 1000*d^2*ew - 5102*r +332*a^2*r + 284*a*b*r + 332*b^2*r + 332*c^2*r + 43*a*b*c^2*r +20*b^2*c^2*r+

284*c*d*r+

43*a^2*c*d*r+

80*a*b*c*d*r

+

(10)

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

43*b^2*c*d*r+

332*d^2*r+

20*a^2*d^2*r+

43*a*b*d^2*r-

63*a*b*ev*r - 63*c*d*ew*r + 3644*r^2 + 75*a^2*r^2 + 724*a*b*r^2 +75*b^2*r^2 + 75*c^2*r^2 - 401*a*b*c^2*r^2 - 124*b^2*c^2*r^2 +724*c*d*r^2 - 401*a^2*c*d*r^2 - 1372*a*b*c*d*r^2 - 401*b^2*c*d*r^2 +75*d^2*r^2 - 124*a^2*d^2*r^2 - 401*a*b*d^2*r^2 + 458*a*b*ev*r^2 +458*c*d*ew*r^2 - 1301*r^3 - 69*a^2*r^3 - 303*a*b*r^3 - 69*b^2*r^3 69*c^2*r^3 + 146*a*b*c^2*r^3 + 42*b^2*c^2*r^3 - 303*c*d*r^3 +146*a^2*c*d*r^3 + 618*a*b*c*d*r^3 + 146*b^2*c*d*r^3 - 69*d^2*r^3 +42*a^2*d^2*r^3 + 146*a*b*d^2*r^3 - 139*a*b*ev*r^3 - 139*c*d*ew*r^3+ 185*r^4 + 12*a^2*r^4 + 39*a*b*r^4 + 12*b^2*r^4 + 12*c^2*r^4 17*a*b*c^2*r^4 - 4*b^2*c^2*r^4 + 39*c*d*r^4 - 17*a^2*c*d*r^4 86*a*b*c*d*r^4 - 17*b^2*c*d*r^4 + 12*d^2*r^4 - 4*a^2*d^2*r^4 17*a*b*d^2*r^4 + 13*a*b*ev*r^4 + 13*c*d*ew*r^4)/1000

この近似の精度についてまず議論する。多項式近似適用前および後の汎関数に，値
a=b=c=d=1，ev=ew=0 を代入し，核間距離 r についてプロットしたものを図 2 に示す。両者が概ねよい一致を示すのは r が 1~2 atomic unit の範囲であって，r がこの有効範囲を超えた場合，多項式近似は不適切であり，別の展開中心点 R0 を使用する必要がある。テイラー展開の次数を上げることによっても，近似が有効な r の範囲を広げることができる。

|Ω[４次テイラー展開]-Ω[厳密]|/|Ω[厳密]|
20
18
16
14

%

12
10
8
6
4
2
0
1

1.2

1.4

1.6

1.8

2

2.2

2.4

核間距離r(a.u.)

図 2 水素分子のエネルギー汎関数のテイラー展開と厳密式とのずれ。
ここで，水素分子の対称性を利用し，波動関数を対称・反対称波動関数の線形結合

φup ( x) = t (exp(− x A ) + exp(− x B )) / π + s(exp(− x A ) − exp(− x B )) / π

φdown ( x) = u (exp(− x A ) + exp(− x B )) / π + v(exp(− x A ) − exp(− x B )) / π

(11)
(12)

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）とする。これは，

a = t + s, b = t − s, c = u + v, d = u − v(13)という変換である。
そうすると，HFR 方程式は(14)に示す連立多項式方程式になる。 (t,s) は上向きスピン
電子の LCAO 係数，(u,v) は下向きスピン電子の LCAO 係数，ev は上向きスピン電子のエ
ネルギー固有値，ew は下向きスピン電子のエネルギー固有値，r は核間距離である。
式(14)は，式(10)を偏微分した後，有理数係数の分母を払って整数係数化したものである。
したがって，これらの式の物理量的次元は不定であり，式(10)の偏微分のものとは厳密には一致しない。以降の議論では，式(14)の根が計算できれば十分であり，定数を掛ける操作は
結果には影響しない。多項式系を整理する目的で，ブッフバーガーのアルゴリズムというガウス消去法に類似の操作を適用する[15]。その際，出発点の多項式系から，多項式対を選び，其々に特定の単項式を掛け，差を取ることによって，二者に共通の単項式を消去した
新しい多項式を生成する。新たに得られた多項式を出発点の多項式系に追加し，線形独立の観点に基づき不要な多項式を除去するという作業を反復することで，グレブナー基底に到達する。多項式間の加減算は，アルゴリズム上，物理量の次元を揃えた上で実行される。
方程式系の解自体は，物理量的次元を終始保つのである。ブッフバーガーのアルゴリズムの反復回数は事例ごとに異なる。計算の結果生成される多項式の物理量的次元は，数式処理の演算回数に応じて変化する不確定なものとなり，当初のエネルギー汎関数の微分の次元とは別のものとなる。今後例示する中間的な多項式の物理量的次元も，式(14)と同様に不定であると考えてよい。無論，逐一計算過程を追跡すれば物理量的次元は判る。

①

∂Ω ∂Ω
+
=0→
∂a ∂b

32*s*u*v*r^4-336*s*u*v*r^3+992*s*u*v*r^2-160*s*u*v*r+40*s*u*v-324*
t*u^2*r^4+2572*t*u^2*r^3-6448*t*u^2*r^2+584*t*u^2*r+19936*t*u^2+
156*t*v^2*r^4-1068*t*v^2*r^3+2248*t*v^2*r^2-80*t*v^2*r-32*t*v^2+26*t*ev*r^4-278*t*ev*r^3+916*t*ev*r^2-126*t*ev*r-7972*t*ev+126*t*r^4-8
82*t*r^3+1748*t*r^2+1896*t*r-12470*t=0

②

∂Ω ∂Ω
−
=0
∂a ∂b

→156*s*u^2*r^4-1068*s*u^2*r^3+2248*s*u^2*r^2-80*s*u^2*r-32*s*u^252*s*v^2*r^4+236*s*v^2*r^3-32*s*v^2*r^2-104*s*v^2*r+48*s*v^2-26*s*
ev*r^4+278*s*ev*r^3-916*s*ev*r^2+126*s*ev*r-28*s*ev-30*s*r^4+330*s*r^3-1148*s*r^2+760*s*r-170*s+32*t*u*v*r^4-336*t*u*v*r^3+992*t*u*v*r^2-160*t*u*v*r+40*t*u*v=0

③

∂Ω ∂Ω
+
=0
∂c ∂d

→156*s^2*u*r^4-1068*s^2*u*r^3+2248*s^2*u*r^2-80*s^2*u*r-32*s^2*u+32**t*v*r^4-336*s*t*v*r^3+992*s*t*v*r^2-160*s*t*v*r+40*s*t*v-324*t^
2*u*r^4+2572*t^2*u*r^3-6448*t^2*u*r^2+584*t^2*u*r+19936*t^2*u+26(14)

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）*u*ew*r^4-278*u*ew*r^3+916*u*ew*r^2-126*u*ew*r-7972*u*ew+126*u*r^4-882*u*r^3+1748*u*r^2+1896*u*r-12470*u=0

④

∂Ω ∂Ω
−
=0
∂c ∂d

→-52*s^2*v*r^4+236*s^2*v*r^3-32*s^2*v*r^2-104*s^2*v*r+48*s^2*v+32*s*t*u*r^4-336*s*t*u*r^3+992*s*t*u*r^2-160*s*t*u*r+40*s*t*u+156*t^2*v*r^4-1068*t^2*v*r^3+2248*t^2*v*r^2-80*t^2*v*r-32*t^2*v-26*v*ew*r^
4+278*v*ew*r^3-916*v*ew*r^2+126*v*ew*r-28*v*ew-30*v*r^4+330*v*r^
3-1148*v*r^2+760*v*r-170*v=0⑤∂Ω= φ up φ up − 1 = 0∂ (ev)→-13*s^2*r^4+139*s^2*r^3-458*s^2*r^2+63*s^2*r-14*s^2+13*t^2*r^4-1
39*t^2*r^3+458*t^2*r^2-63*t^2*r-3986*t^2+1000=0⑥∂Ω= φ down φ down − 1 = 0∂ (ew)→13*u^2*r^4-139*u^2*r^3+458*u^2*r^2-63*u^2*r-3986*u^2-13*v^2*r^
4+139*v^2*r^3-458*v^2*r^2+63*v^2*r-14*v^2+1000=0

⑦

∂Ω
=0
∂r

→312*s^2*u^2*r^3-1602*s^2*u^2*r^2+2248*s^2*u^2*r-40*s^2*u^2-104*s^2*v^2*r^3+354*s^2*v^2*r^2-32*s^2*v^2*r-52*s^2*v^2-52*s^2*ev*r^
3+417*s^2*ev*r^2-916*s^2*ev*r+63*s^2*ev-60*s^2*r^3+495*s^2*r^2-11
48*s^2*r+380*s^2+128*s*t*u*v*r^3-1008*s*t*u*v*r^2+1984*s*t*u*v*r-1
60*s*t*u*v-648*t^2*u^2*r^3+3858*t^2*u^2*r^2-6448*t^2*u^2*r+292*t^2*u^2+312*t^2*v^2*r^3-1602*t^2*v^2*r^2+2248*t^2*v^2*r-40*t^2*v^2+
52*t^2*ev*r^3-417*t^2*ev*r^2+916*t^2*ev*r-63*t^2*ev+252*t^2*r^3-13
23*t^2*r^2+1748*t^2*r+948*t^2+52*u^2*ew*r^3-417*u^2*ew*r^2+916*u^2*ew*r-63*u^2*ew+252*u^2*r^3-1323*u^2*r^2+1748*u^2*r+948*u^252*v^2*ew*r^3+417*v^2*ew*r^2-916*v^2*ew*r+63*v^2*ew-60*v^2*r^3+495*v^2*r^2-1148*v^2*r+380*v^2+740*r^3-3903*r^2+7288*r-5102=0
まず，本法で実際に第一原理計算が可能なことを示すため，R=7/5 を仮定して解く。方程式(14)において，

∂Ω
=0
∂r

(15)

5r − 7 = 0

(16)

という式を
と置き換える。この方程式系から生成された辞書式順序グレブナー基底は(17)である。単項

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

式順序は，s<t<u<v<ev<ew<r の順とする。今後の計算ではこの単項式順序を使う。そのため，グレブナー基底(17)の各要素において，変数は，これとは逆に，r, ew, ev, v, u, t, s の順序で現れる。HFR 方程式の段階では曖昧であった変数間の相互関連が，順次抽出されてくることに注意しよう。
J[1]=r-1.4
J[2]=0.000000000000000044923679950280179834153752545837120050869
800531719*ew^6+0.0000000000000000837848476849489179918230616638
79813241022968255965*ew^5+0.00000000000000005522431406589942780
2789765894852689587104815167352*ew^4+0.000000000000000014691161
157734078517931578598749996418029499346849*ew^3+0.0000000000000
000011979351185242681597299140874114661346717404726834*ew^2-0.00
000000000000000000839869963562374452466573406585698958493547765
49449*ew-0.00000000000000000000037230409053155950263983522183691
557632339096670694
J[3]=0.000000000139628765357057157466220713819530977198921180700
71*ev+0.00000032173656936235047932064238004296019881075131225526*ew^5+0.0000003329165247773600520718203430909808162869659481933
1*ew^4+0.000000065891103236483021260956200510555934469425614980
251*ew^3-0.00000001450219479632473049589549258544954048787083075
49*ew^2-0.000000002589663947030314527682380814926503763029302261
6944*ew+0.00000000005081919751322588124435323091396968621359149
4118136
J[4]=0.0035337421249281116351414330877752843125543561935616*v*ew^4+0.0021801710635532512041192523955501421176113052868825*v*ew^
3+0.00024683412076704120861849546037854544501562730166012*v*ew^
2-0.0000014554651684258252511501850703694711773679621730295*v*ew-0.000000075203364912461984107689153015151456273923670860722*v
J[5]=0.000000000000002469182263912597327197080785436650088580523
0838289*v^2-0.00000000000471902173197923393226385801209311776724
87391351614*ew^5-0.000000000007423372586143701255637248508913491
2758173742859754*ew^4-0.0000000000037044087563918384281802854504
031102187085751576865*ew^3-0.00000000000055056937948011551503818
116670012263028105621886299*ew^2+0.0000000000000028924710974820
572210152554673921687498485750622516*ew-0.0000000000000048519950
027595085878960476168529106170995782184575
J[6]=0.00082444529533914000480121749252034857614275087163634*u*ew^4+0.0015402471286109903529960234596133438365658562434456*u*ew^3+0.0010186155101412674041097562139745242960275920262626*u*ew^
2+0.00027330147151793449108809755546973584926704784843595*u*ew+
0.000023152649811839442946188162579039171063101525686031*u
J[7]=0.15500652085373657278272155048840425269081666861422*u*v*ew(17)

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）^2+0.096124486745923798453025086330083252761630136665135*u*v*ew+0.011178119900970398824611330746954877639528183070683*u*v
J[8]=0.000000000000002469182263912597327197080785436650088580523
0838289*u^2+0.00000000000066185419621956626251845869357753738900
545315696192*ew^5+0.0000000000010411459355962387028179380989969
086880336301691908*ew^4+0.0000000000005195522757006096209849875
8121689908883959384515193*ew^3+0.000000000000077218685315543857
629354030269499939648408016279623*ew^2-0.00000000000000040567605
788697774347415485721964229236559337475379*ew-0.0000000000000000
23368741490554450106147082788613043426219204168902
J[9]=0.061970918133163216567183870035153503888001513514567*t*ew^4+0.07786967298515361066019891920879105862334446468288*t*ew^3+0.0
29529287545743808709620447245061926950863619513321*t*ew^2+0.0032
178375021467349512602370986204354827641359493411*t*ew+0.00004345
6853031823070001289851734608283302642795735406*t
J[10]=0.05187413667297784999868866991024835274577850308566*t*v*ew^3+0.032981416568891625365664312121206773685246825760142*t*v*ew^
2+0.0042447701517134049943862889216523608128280851170331*t*v*ew+
0.000058600461674134415884220932579569266784799169692295*t*v
J[11]=0.021791703354265359880262646937291988076821313685722*t*u*ew^3+0.027041037783307743984286616783943211212082455989442*t*u*ew^2+0.0099601988616932931997530597003570949767046400427524*t*u*e
w+0.00097550648450924773054739454162691886820402290190827*t*u
J[12]=0.00000000000000246918226391259732719708078543665008858052
30838289*t^2-0.0000000000019897101060519652107395125478208499820
923810270269*ew^5-0.00000000000170255700897566555354942878105283
79598119659773602*ew^4-0.000000000000023483592200304199348742150
925412718812994059865964*ew^3+0.0000000000001967372781280848501
5748182293664380763085511692487*ew^2+0.000000000000019786078457
023833052592026181105821026866390535766*ew-0.0000000000000004421
9079290739404492269360637167466142801279295227
J[13]=0.056084833175126678686494811325106385425518369783604*s*ew^2+0.034127439269225675895016122448171947840127214030889*s*ew-0.
00066282300441698232705078773127382486945402101120119*s-0.040893
885020249465973553191938478590309960079399638*t*u*v*ew-0.0290103
01971762382682320743416650072000710482565746*t*u*v
J[14]=0.021791703354265359880262646937291988076821313685722*s*v*e
w-0.00041052862271488237533112188857801258113775891824163*s*v-t*u*ew^2-0.91523243331780979308516147067869265402585316644918*t*u*e
w-0.18279835739609284049315718860436264109118081404686*t*u
J[15]=0.082689621743132526006561521010847989515268982874621*s*u*e
w+0.05187413667297784999868866991024835274577850308566*s*u+t*v*e

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

w^2+0.45088158483424862414602508768019029048144708971225*t*v*ew+0.0068176823481038061179522681529012450074191772534997*t*v
J[16]=0.061970918133163216567183870035153503888001513514567*s*u*v+7.54843402550665153673201812730961379367248454713989*t*ew^3+9.2
977596625794358858261233928027031703300204455612*t*ew^2+2.933336
92208165449675224101101090746851954196539922*t*ew+0.043698233453
379147442339681621264612101282386252875*t
J[17]=0.15500652085373657278272155048840425269081666861422*s*t-u*v*ew-0.3100659450211982179038598636924866298895502223605*u*v
J[18]=0.00000000000000246918226391259732719708078543665008858052
30838289*s^2+0.0000000000141866370031792066049786475011937123318
71881329982*ew^5+0.00000000001213923485139370932693489804708248
0718477074921104*ew^4+0.000000000000167438058973051827655402125
10939603762998421893102*ew^3-0.000000000001402737183325013576178
8614236060225048185772213358*ew^2-0.0000000000001410747786486290
6344428789190998451814760012308984*ew-0.000000000000001865792945
3315597187256856909000195939743454954377
グレブナー基底(17)を三角化すると方程式(18)となる。ここには五つの三角化された連立方程式が含まれる。一つの連立方程式は，七つの式を含み，変数は，r, ew, ev, v, u, t, s の順序で追加されていく。
(18)

[1]:
_[1]=r-1.4_[2]=0.082689621743132526006561521010847989515268982874621*ew+0.05187413667297784999868866991024835274577850308566_[3]=0.00000000013962876535705715746622071381953097719892118070071*ev+0.0000000000021872894505215127104699976329395744001395648217513_[4]=v_[5]=0.0000000000000024691822639125973271970807854366500885805230838289*u^2-0.00000000000000070387276012025253831506975421079541060390364383305_[6]=t_[7]=0.0000000000000024691822639125973271970807854366500885805230838289*s^2-0.0000000000000050186141759442129979160381879133091558873472043949[2]:_[1]=r-1.4_[2]=ew-0.018838757853893431115706795065601087293579573640499_[3]=ev-0.018838757853893431115706795065601087293579573640499_[4]=v^2-2.03250049593012100694952569568426974696181825868366_[5]=u_[6]=t_[7]=0.0000000000000024691822639125973271970807854366500885805230838289*s^2

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）-0.0000000000000050186141759442129979160381879133091558873472043949[3]:_[1]=r-1.4_[2]=ew^2+0.620131890042396435807719727384973259779100444721*ew+0.072113868754708844917382916961852739039135803532834_[3]=ev-ew_[4]=v^2-14.80995023428197133951905921507729315089177369253934*ew-5.77765333933354558148702071632256029046854951290411_[5]=u^2+2.07713129226284416025022121618132789954954241021854*ew+0.52526673231711754540325634713814260954516716824085_[6]=t^2+2.07713129226284416025022121618132789954954241021854*ew+0.52526673231711754540325634713814260954516716824085_[7]=s+8.55670765807561749977207690461050758297371934231746*t*u*v*ew+6.24834985328317561951025512690840116358160281553919*t*u*v[4]:_[1]=r-1.4_[2]=ew+0.62075494398358690094167914213207711792332315380224_[3]=ev+0.62075494398358690094167914213207711792332315380224_[4]=v_[5]=u^2-0.28506310384917288370291954789447830487339321330884_[6]=t^2-0.28506310384917288370291954789447830487339321330884_[7]=s[5]:_[1]=r-1.4_[2]=ew+0.015665034671961754863974784720813826678603082481373_[3]=ev+0.6273355178007699660610156274085545071738486487756_[4]=v^2-2.03250049593012100694952569568426974696181825868366_[5]=u_[6]=t^2-0.28506310384917288370291954789447830487339321330884_[7]=s式(17)(18)は，計算機代数ソフト Singular[11]の計算結果をそのまま示したものであり，記法は Singular の出力に従う。(17),(18)の記法の不統一は，それぞれ別々の作者がプログ
ラムの組み込み関数を作成しているためである。
式(17),(18)の係数は，非常に桁数の大きい少数を含む。これは，グレブナー基底を生成す
るアルゴリズムに固有の問題であって「係数膨張」と呼ばれる。ブッフバーガーのアルゴ
リズムに基づき，出発点の多項式系に対し，加減乗除を繰り返し，多項式の追加・除去を行うことで，グレブナー基底に到達する。その際に生成される中間的表現において，長大な多項式が生じ，その数値係数の一部が著しく小さく，また，別の一部が著しく大きい値《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

をとる場合がある[19]。この係数の不揃は最終結果にまで残るのである。求解精度を確保するために，任意精度の数値計算を行い，桁数の大きな数値係数を使用する必要がある。
式(18)の解は，複素数解も含むが，物理的に意味のある実数解のみを表 1 に示す。上向き・下向きスピンを持つ二つの電子が対称・反対称波動関数の二準位のどちらかに位置する，四つの組合せに対する波動関数が得られた。すなわち，基底状態と励起状態の両方が得られたことになる。
解１
電子１～対称軌道
電子２～対称軌道
基底状態

解２
解３
解４
電子１～反対称軌道 電子１～対称軌道
電子１～反対称軌道
電子２～対称軌道
電子２～反対称軌道 電子２～反対称軌道
励起状態
励起状態
励起状態

s（電子１係数）0.00000-1.42566
0.00000t（電子１係数）-0.53391
0.00000-0.53391u（電子２係数）-0.53391-0.53391
0.00000v(電子２係数）0.00000
0.00000-1.42566ev（電子１固有値）-0.62075-0.01567-0.62734ew（電子２固有値）-0.62075-0.62734-0.01567r（核間距離）1.40000
1.40000
1.40000
Etot（全エネルギー）-1.09624-0.49115-0.49115
表1 方程式系(18)の実数解。電子１（上向スピン），電子２（下向スピン）である。

基底状態のみを得るために，方程式(13)に，式

s=v=0

(19)

を付け加える。（これは，この事例のみに適用可能な手段である。一般に，基底状態は，全固有値の和が最小になるという条件を満たす解として得ることができる。基底状態を特定するためには，固有値のみをまず計算すればよい。このとき，辞書式順序のグレブナー基
底またはこれを三角化した方程式系を作る際に，固有値のみを含む方程式群をまず作成する。そして，固有値のみの方程式を解けばよい。）こうすることで方程式系を簡単なものに置き換え解くことができる。このとき，核間距離 r 及び波動関数は同時に最適化される(カー・パリネロ型解法)。r に関する方程式を(20)，実数解を表 2 に示す。
[1]:
8942144364*r^23-435341589039*r^22+9813157241157*r^21-134458128500631*r^20+1251986164962728*r^19-8584760758387395*r^18+48176522279858253*r^17-254992901607817871*r^16+1360184656773665254*r^15-668541270541318423
5*r^14+26848712421674517351*r^13-82265960807423324641*r^12+1853704803
18135661708*r^11-295651827763150999108*r^10+307426892321213994312*r^9-148683667595876075980*r^8-97338526988608612178*r^7+24577251883657979
1529*r^6-200002425723099153061*r^5+47298638179277635737*r^4+460063481
88804187952*r^3-41646082527529600720*r^2+13118922400543578496*r-18697
47053688110592=0[2]:10313892*r^15-376866027*r^14+6245669754*r^13-61144647973*r^12+3877646
99571*r^11-1646957525797*r^10+4691411679124*r^9-8760215434992*r^8+102
81598671237*r^7-7316755042677*r^6+3784010771997*r^5-2194016637700*r^4299532295668*r^3+1482785614608*r^2-746000940352*r+36100845312=0(20)-1.42566
0.00000
0.00000-1.42566
0.018840.01884
1.40000
0.15503

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）[3]:
10313892*r^15-376866027*r^14+6245669754*r^13-61144647973*r^12+3785967
95571*r^11-1405917509797*r^10+1917133055124*r^9+9044372533008*r^8-551
18065080763*r^7+103464030245323*r^6+92432281739997*r^5-77079701000570
0*r^4+1063674493728332*r^3+652030557238608*r^2-2854269358708352*r+195
4998898509312=0

r

ev

t

実数解１

-1.812

-6.6

0.846

実数解２

1.652

-0.578

0.545

実数解３

6.010

-17.5850.983

表 2 方程式系(20)の実数解
解は，正負の実数値および複素数値を含む。意味のある解（r＞0）は，二つあるが，テイラー展開の有効範囲内にある解は，r～1.6 のみである。この解と実測値(r~1.4)のずれは，第一に四次テイラー展開と数値係数の有理化の粗さに由来する数値誤差によるものである，
また，軌道指数を最適化していない試行関数を用いていることにもよる。
次に本法を用いた逆問題の解法例を示す。水素分子の占有，非占有状態の固有エネルギ
ー差がある特定の値を持つ状況を想定し，その値の出るような核間距離を求めることを考える。この実施例は，例えば，格子定数を変化させて望ましいバンドギャップを求めるような「逆問題」の雛型である。今回は RHF(Restricted Hartree-Fock) 計算を行う。占有・非占有状態の波動関数を

φocc ( x) = s(exp(− x A ) + exp(− x B )) / π
φunocc ( x) = t (exp(− x A ) − exp(− xB )) / π

(21)
(22)

とし，固有値を eocc，eunocc とした。この場合，占有状態を決める方程式は，上にあげた UHF
順問題計算と同様の方式で得ることができる。さらに非占有状態を決める方程式と占有・非占有状態間の直交条件が追加される。このとき必要な方程式系は(23)に示す。式詳細は省略する。

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

(23)

∂Ω
=0
∂s
∂Ω
=0
∂t
∂Ω
= φocc φocc − 1 = 0∂ (eocc )∂Ω= φunocc φunocc − 1 = 0∂ (eunocc )φocc φunocc = 0
eocc − eunocc = E gap
例として，Egap=eunocc-eocc=0.9 を与える r を求めてみよう。実数解は，r＝－1.103，0.307，
1.643，3.958 である。テイラー展開の有効範囲にある正の実数解は r=1.643 のみである。rを連続的に変化させた場合の占有・非占有準位の固有値を図 3 に示すが，この図と比べて適正な結果が得られたことがわかる。
0.8
0.6
0.4
0.2
固
有
値

0
-0.2

占有準位
1

1.2

1.4

1.6

1.8

2

非占有準位

-0.4
-0.6
-0.8

核間距離(R)

図3水素分子の占有・非占有状態の固有値の核間距離依存性
では，この場合，構造は安定であるか。無論，力場計算を行えば，簡単に判定ができる。
一方，Symbolic-Numeric 計算の立場からは，次のような判定を行うことができる。方程式
系(23)に汎関数の核間距離 r に関する極値条件

∂Ω
=0
∂r

(24)

を追加する。このときのグレブナー基底の計算結果は，多項式の集合と見たとき，ただ一つの要素(定数多項式)のみを含む{1}となる。
（今回用いた計算機代数ソフト Singular のグレ
ブナー基底では，定数多項式は１となるように規格化される。）グレブナー基底の零点は，元来の多項式系の根を与えるが，そもそも多項式とみたとき{1}は零にならない。したがって，今回考える問題は解を持たないと結論できる。

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

次に，これとは別種の数値解法である，Stickelberger の定理を利用した解法の実施例を示す。例題は，水素分子の基底状態の構造最適化であり，RHF 計算を実行する。試行関数を(25)

φocc (r ) = t (exp(−rA ) + exp(−rB )) / πとする。固有値を ev，核間距離を r とする。前述の実施例と同様に，ハートリー・フォック方程式は，変数 t，ev，r で表される多項式となる。この多項式の集合は，数学的には多項式環 A = R[t, ev, r] における零次元イデアル I を構成し，イデアル I の零点は剰余環 A/Iに対応する。R は実数体 R である。A は R 上有限次元のベクトル空間であり，このベクトル空間の基底は，t，ev，r の単項式によって表わされる。これを(26)に示す。
b[1]=t*ev*r^3b[2]=t*r^4b[3]=t^3*ev

b[4]=t*ev^3b[5]=t^3*rb[6]=t*ev^2*r b[7]=t*ev*r^2b[8]=t*r^3b[9]=ev*r^3b[10]=r^4b[11]=t^3b[12]=t^2*ev

b[13]=t*ev^2 b[14]=ev^3b[15]=t^2*rb[16]=t*ev*r

b[17]=ev^2*rb[18]=t*r^2

b[19]=ev*r^2b[20]=r^3b[21]=t^2b[22]=t*ev

b[23]=ev^2b[24]=t*rb[25]=ev*rb[26]=r^2b[27]=tb[28]=evb[29]=rb[30]=1(26)

ベクトル空間 A/I の基底に対する変数ｔの積演算を表す変換行列 mt を(27)に示す。
(27)

変数 t，ev の積演算を表す変換行列も同様に計算することができる。変換行列 mt の固有ベ
クトルv を求め，関係式

mt ⋅ vξ = ξ t ⋅ vξm ev ⋅ v ξ = ξ ev ⋅ v ξ

(28)
(29)

から t，ev の値を固有値 ξt，ξev として求めることができる。このとき，実数解は，表 2 に示した多項式の三角化解法によって得られたものと同等である。

《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

4. まとめ以下，本手法の効果を順に挙げる。
(利点 1) 本手法では，基礎方程式系であるハートリー・フォック・ロータン方程式及び拘束条件・最適化条件を多項式系で表すため，電子状態計算と付加条件とがシームレスに結合される。従来法では，固有値計算のループ，自己無撞着計算のループ，最適化計算のループが入れ子になっており複雑である。一方，本手法においては，これらの入れ子になったループが，連立多変数多項式系の根の探索という単一のフローに統合されており，数値計算時に見通しのよいものであり，手短に必要な情報に到達することが可能である。従来法で行われたような，原子核・波動関数の自由度に関する独立の緩和計算を交互に反復する必要がない。本法では原子核と波動関数のダイナミクスが強く結合するような物理現象
も，断熱近似の域を超えて，取り扱うことができるであろう。
(利点 2)
本手法においては基礎方程式系である HFR 方程式及び拘束条件・最適化条件を多項式系で表すもう一つ利点は次の通りである。多項式系は，未知変数間の関係式を与えるものである。したがって，多項式系によって記述される未知変数を求めたい場合，計算の入力・出力に，それぞれ，適切な変数を振り分けることができる。すなわち，入力値は構造であり，出力は電子状態である，という従来法の逆問題の枠組みに限定されない。したがって，従来法でいう順問題・最適化問題・逆問題の区別が解消され，すべて順問題として統一的に解くことができる。逆問題を扱う場合には，問題の適切性を議論しなくてはならない。本方法では，基礎方程式系をグレブナー基底に変換した段階で逆問題の適切性，すなわち解の存在の検証が可能である。代数学のイデアル理論に基づき，グレブナー基底の零値を与える解集合が存在するか否かが判る。解集合が存在する場合，それは孤立点か，あるいは，一次元以上の次元を持つ集合であるかがわかる。解集合が孤立点の場合，複素数あるいは実数の範囲で，幾つ解があるかも判定できる。
(利点３)本手法を自己無撞着場計算に適用する利点は次のようなものである。三角化された多項式系を一変数ずつ数値的に求めることによって計算できる。収束性に問題の発生することの多い，多次元行列に対する固有値計算も，平均場の逐次近似計算も不要である。
(利点４)本手法によって，最適化計算する場合の利点として，次の点も挙げることができる。例として，原子構造の最適化を考える。変数消去の順番を工夫し，原子座標のみの多項式を得ることができる。この多項式の根が安定構造である。すなわち，SCF(Self-Consistent Field,自己無撞着場)計算を行わずに，最適な原子配置を求めることができるのである。他のパラメータを最適化する際も同様の利点がある。
本研究は「分子積分の多項式近似」「分子軌道代数方程式」という概念が現実の系の電子状態計算に適用可能であるのみならず，既存手法では不可能あるいは困難が大きい問題も視野に収めうるものを示すものである。ただし，現状では，必ずしも精度の良い計算ではない。この理由の一つは，計算機性能の制約である。記号処理の計算負荷が大きく，このために，近似多項式の次数を小さくし，また，数値係数を低精度で有理数近似を行わざるを得ない。さらに，記号処理上の原理的な問題もある。式(14)に見るように，出発点とした多項式連立方程式系の数値係数の大きさはほぼ揃っている。しかし，これを数式処理し，グレブナー基底を生成すると，数値係数の大きさは極めて不揃いなものになる。この種の，係数が極端に膨れ上がる事態は，計算コストを増大させるばかりでなく数値解上の精度悪化の原因となるが，現時点では有効な対策は少ない。しかしながら，今後の計算機性能の《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

向上と数式処理技術の改良によって，本手法における十分な精度の達成と大規模計算への適用が可能になるであろう。
謝辞
本研究において，分子積分の解析式およびテイラー展開式は数式処理ソフト
Mathematica を用いて作成した。分子軌道代数方程式の記号的―数値的求解は，計算機代数システム Singular を用いて実施した。
筆者は，分子軌道代数方程式および分子積分の多項式近似法の概念についてご議論を頂いた同僚の安井潤博士に感謝の意を表します。また，同氏に，STO 分子積分生成用数式処理プログラムを提供頂いたことを感謝します。
参考文献
[1]

藤永茂「分子軌道法」岩波書店， 1980 年。

[2]

金森 順次郎， 川村 清 ， 米沢 富美子 ， 寺倉 清之 ，
「岩波講座 現代の物理学 7
固体－構造と物性」

[3]

岩波書店

1994。

Dominik Marx and Jürg Hutter, “Ab Initio Molecular Dynamics: Theory and
Implementation”, published in “Modern Methods and Algorithms of Quantum
Chemistry, Proceedings, Second Edition”, J. Grotendorst (Ed.), John von
Neumann Institute for Computing, Jülich, NIC Series, Vol. 3, ISBN3-00-005834-6, pp. 329 - 477, 2000.

[4]

A. D. McLean and M. Yoshimine, “Tables of Linear Molecule Wavefunctions”,
IBM J. Res. Dev. Supplement, November, 1967.

[5]

Jun Yasui and A. Saika, J. Chem. Phys 76, 468(1982).

[6]

安井潤，「量子材料計算の基礎としての分子積分」，DVXα研究協会会報，
Vo1.22,No.1&2，7‐11(2009).

[7]

安井潤，”Polynomial Expressions of Molecular Integrals Functionals overSlater-Type-Orbitals and Its application to the Extension ofHartree-Fock-Roothaan Equation”, DVXα研究協会会報，Vol.23,No.1&2，54‐59(2010).

[8]

安井潤，”Algebraic Molecular Orbital Equation”, DVXα研究協会会報，
Vol.23,No.1&2，47‐54(2011).

[9]

J. Yasui, to be published in Progress in Theoretical chemistry and Physics(Springer) (2012).

[10] 数式処理ソフト Mathematica, http://www.wolfram.com.
[11] 計算機代数システム Singular, W. Decker, G.-M. Greuel, G. Pfister, H.
Schoenemann, http://www.singular.uni-kl.de/.
[12] H. M. Möller, “On decomposing systems of polynomial equations with finitely
many solutions”, Appl. Algebra Eng. Commun. Comput. 4, 217 - 230, (1993).
[13] D. Lazard, “Solving zero-dimensional algebraic systems”, J. Symb. Comp. 13, 117 132, (1992).
[14] Frank Sottile , in “Computations in algebraic geometry with Macaulay 2” editedby David Eisenbud, Daniel R. Grayson, Michael E. Stillman, and Bernd《寄稿論文》

物性研究・電子版 Vol, 1, No. 3, 013101 （2012年11月号）

Sturmfels, published by Springer-Verlag in September 25, 2001, as
number 8 in the series “Algorithms and Computations in Mathematics”,
ISBN 3-540-42230-7, p101-130.
[15] 丸山正樹「共立叢書・現代数学の潮流 グレブナー基底とその応用」

共立出版，

2002 年。
[16] M. P. Barnett, J. F. Capitani, J. von zur Gathen and J. Gerhard, “Symbolic
calculation in chemistry: selected examples”, Int. J. Quant. Chem. 100 (2) 80 - 104,
2004.
[17] M. J. S.Dewar, E.G. Zoebisch, E.F.Healy, J.J.P. Stewart, J. Am. Chem. Soc. 107,
3902 (1985).
[18] J.J.P. Stewart, J. Comput. Chem. 10 (1989).
[19] M.Brickenstein,
http://www.mathematik.uni-kl.de/~zca/Reports_on_ca/35/paper_35_full.ps.gz.