FVMO法による1s型Gauss型関数系のみを用いた分子軌道法の可能性

笹金 光徳


Return

1 はじめに

Gauss型関数(Gaussian type function = GTF)は,単独では中心から遠距離のなだらかな裾野やs軌道におけるr = 0の尖端(cusp)をうまく表現できないが,重ね合わせることによって,原子の軌道を程よく表現でき,なおかつ積分計算が容易であることから,原子・分子の電子状態を計算する分子軌道法のプログラムで専ら利用されるようになった.また,藤永の水素原子に対する研究によって,各方位量子数lの軌道を記述する場合,主量子数nが最低のGTF (n = l + 1),すなわち,節のないGTFの重ね合わせにより,n > l + 1である原子軌道も効率よく表現できることが報告されている[1]. たとえば,水素原子の2s軌道をGTFの重ね合わせで表現するときに,1s型GTF (∝ exp[-ar2])を用いた場合のほうが2s型GTF(∝ r exp[-ar2])を用いるより合理的であることが示されている.現行の分子軌道法プログラムでは,基底関数は軌道指数(orbital exponent)の異なるいくつかGTFまたはその線形結合から構成されるが,こういった経緯から,s軌道はすべて1s型,p軌道はすべて2p型,d軌道はすべて3d型に統一されている.なお,ここで言う1s型,2p型,3d型とは,指数因子の外に等方的なrが現れない関数系であることを意味しているのであって,決して空間的な広がりを示唆しているわけではない.
また,分子の計算において,分子軌道における原子軌道からの歪みを正確に表現するためには,分子を構成する原子が本来持っている最上位の方位量子数より高次の方位量子数で軌道指数があまり小さくない関数を分極関数(polarization function)として加えることが必要である.また,分極関数の利用は,電子相関を考慮したさらに高精度な計算においても必須である.
一方,基底関数を構成するGTFの軌道指数は,通常はHartree-Fock-Roothaan-Hall法による原子の基底状態の計算から定められており,分子を構成するときに全体として起こる若干の軌道収縮や,電気陰性度の差に起因する分子内原子間の電子の移動に伴う若干の電子雲の収縮や膨張を記述するために,スケール因子(scaling factor)を軌道指数の調節に利用できるようなっている.しかし,陰イオンの計算や励起状態の計算となるとスケール因子での調節では不十分であり,電子雲の膨張を記述できる軌道指数の小さな関数を分散関数(diffuse function)として加えるのが常套手段である.たとえば,扱う現象が価電子励起(valence-electron excitation)の領域であれば,各方位量子数について,1~2個程度の分散関数を加えればかなり精度が改善される.しかるに,リドベルグ励起(Rydberg excitation)や内殻励起(core-electron excitation)を扱う場合に,分散関数で対処するだけでは,不十分で,軌道指数の何らかの最適化が望まれる.それと類似して,電磁場存在下での分子の電子状態の研究,言い換えれば,(超)分極率[(hyper)polarizability]や(超)磁化率[(hyper)magnetizability]を扱う研究では,比較的低エネルギー領域のリドベルグ励起や内殻励起の寄与が無視できないので,通常の分子軌道法を用いる場合,信頼性の高い結果を得るためには,分極関数と分散関数をふんだんに加える必要がある.
以上で述べたように,電子相関を考慮しない原子の基底状態の計算である程度最適化して得られた軌道指数をもつGTFより構成される基底関数を用いる通常の分子軌道法の計算では,分極関数,スケール因子,分散関数の存在を無視することはできない.しかも,その選び方にある程度の職人芸的経験則が必要であり,既存の分子軌道法プログラムの多くの利用者が自らの意志で的確な選択をすることは困難である.その結果,多くの場合にはプログラムに内蔵されている標準的な値をそのまま利用することになるが,すべての系においてそれで問題がないとは言い難い.
1997年に森が提案し,共同研究者と発展させたFVMO法[2 - 4]は,本質的のこういったすべてのことに対して職人芸的経験則を排除し,本当の意味の「非経験的分子軌道法」の実現を目指す試みであると解釈できる.すなわち,分子や原子団イオンのさまざまな状態において原子軌道からの歪みが,軌道指数と軌道中心の位置の双方を最適化することによって本質的に解決できるのである.筆者は,森らの開発したFVMO法による分子軌道法のプログラムであるGAMERAにROHFのパートやC2vの対称性を考慮したコードを付け加えた.

2 基本的な考え方

電場の加わった系に対する有限場摂動法(finite field method)を適用して,FVMO法で分極率および超分極率を算出する場合に,分極関数や分散関数を通常の分子軌道法よりも相当節約できることがわかっている.ここでは,論点を明確にするために,原子の計算を例にとって,FVMO法の特長を浮かび上がらせる.
立川らは,1999年にProceedings of International Workshop on Molecular Design of Photonic Materials (MDPM99)で,Hartree-FockレベルのHe原子の分極率が,1s型GTFだけで高精度に得られることを発表しているし,Full CIレベルのHe原子の分極率が(7s3p2d1f)でかなり正確に得られることを示している[5].また,筆者は,Hartree-FockレベルのBe原子の分極率aおよび第2超分極率gの計算をFVMO法で行い,International Conference of Computational Methods in Sciences and Engineering 2007 (ICCMSC 2007)において,(16s7p)でほぼ,Hartree-Fock極限値と同等の結果が得られることを発表した[6].
分極率aを正確に算出するためには,基底状態からの双極子遷移が許容となる主要な励起状態が忠実に記述できることが必要である.同様に,第2超分極率gを正確に算出するためには,基底状態からの2光子吸収で得られる状態,言い換えれば,aの算出に必要な主要な励起状態からの双極子遷移が許容となる多くの状態が記述できることが必要である.つまり,HeやBeの基底状態ではs軌道にしか電子が収容されていないので,aの算出にはp軌道への遷移がすべての許容遷移となる.一方,gの算出では,s軌道からp軌道への遷移の他に,p軌道からd軌道への遷移,p軌道からs軌道への遷移の寄与が必要となる.
つまり,軌道中心が固定されている通常の分子軌道法の場合には,Hartree-Fockレベルであっても,HeやBeの正確なaの算出には,2p, 3p, 4p, 5p ...が記述可能な分散性分極関数を含む必要がある.同様に,正確なgの算出には,3d, 4d, 5d, 6d ...が記述可能な分散性分極関数が必須なのである.それに対し,FVMO法を用いるとこのような分極関数を相当節約できるという事実が意味しているものは,1s型GTFの軌道中心をずらすことによって擬似的にp軌道の描像が表現可能であるということであり,2p型GTFの軌道中心をずらすことによって擬似的にd軌道の描像が得られているということである.また,1s型GTFの軌道中心をずらすことによって擬似的にd軌道の描像もある程度表せるのかもしれないという推測もできる.
そこで,本研究では,FVMO法プログラムGAMERAを用い,1s型GTFの軌道指数と軌道中心を最適化することによって,H原子の2p軌道や3d軌道をどの程度うまく表せるかどうかを試みることにした.
今回の試みは,上述のような外部電場への応答時のFVMO法の優位性に関する解析がきっかけとなっているが,p軌道やd軌道を1s型GTFで表現しようという試み自体は新しいことではない.1956年にPreussが提案したGaussian Lobe法[7]は,1960年代にWhittenらによって進展し[8 - 10],1s型GTFだけで基底関数系を形成するという試みがなされた.彼らの研究では,1s型GTFからp軌道やd軌道を定めるとき,エネルギーを最低にするよう軌道指数と軌道中心を直接最適化しているわけではなく,Clementiらが公表した高精度STO展開[11, 12]との重なりを最大化するよういくつかの制限下でフィッティングを行っているに過ぎない.また,定められた1s型GTFだけによる分子軌道計算を勧めているが,現在Gaussian Lobe法が標準的な手法として定着しているとは言い難い.これは,彼らが意図した「1s型GTFだけを使うことによる積分計算とアルゴリズムの簡便化とそれに伴うプログラムの高速化」が,支持されなかったことを示唆している.
これに対して,本研究では,変分原理に基づいて,軌道指数と軌道中心の最適化を行っている.さらに,この最適化は,Gaussian Lobe型の基底関数を定めることが最終目的ではなく,FVMO法の発展性を確認するための試行であるという点で,過去の研究とは全く意義が異なっている.
Table 1は,H原子の1s軌道を1s型GTFの線形結合で,2p軌道を2p型GTFの線形結合で,3d軌道を3d型GTFの線形結合でというように通常の方法で表したときの最適値を,展開項数を増やしながらまとめたものである.本表の軌道指数の最適化を含むすべての計算は,GAMERAで行った.得られたエネルギー値の計算精度については,小数点以下14桁まで正確である.その根拠は,すべての計算において,Virial比( -V/T )が14桁目まで,2.0に一致するからである.

Table 1. Accurate electronic energies (in a.u.) of H atom by optimized GTF expansions (calculated by GAMERA).
N1s state (2S)2p state (2P)3d state (2D)
1-0.42441318157839-0.11317684842090-0.05173798784956
2-0.48581271661627-0.12328871335863-0.05507968023553
3-0.49697925270505-0.12472760095647-0.05549257374609
4-0.49927840571435-0.12495181044060-0.05554644466053
5-0.49980983223189-0.12499056856056-0.05555410424618
6-0.49994557039665-0.12499798064052-0.05555530248720
7-0.49998329778916-0.12499953246659-0.05555550769972
8-0.49999456139075-0.12499988413336
9-0.49999813603794-0.12499996952343
10-0.49999933197686
11-0.49999975096518
12-0.49999990384689
Exact-0.5-0.125-0.05555555555556

3 H原子の2p軌道を1s型GTFで表す試み

動径方向の関数曲線の形を別にすれば,s軌道とp軌道の大きな違いは,もちろん対称性である.つまり,p軌道を1s型GTF

で表す場合,同じ軌道指数の2つのGTFを組みにして,空間的に反対称化する必要がある.本研究では,便宜上pxを1s型GTFで表すことにするが,k番目の擬似的2px展開関数 は,規格化因子は省略して次のように表される.

この式の軌道指数akと原子が置かれている原点からの距離xkを各展開項についてFVMO法で最適化した.(2)式の対称性の解は,ROHF計算の2番目の解として得られるが,C2v用のコードを活用して,対称性が崩れることを防止しながら,規約表現B1の最低解として収束させた.
Table 2は,1s型GTFの組みを用いて最適化されたエネルギーである.Table 1と比較すると同数の2p型GTFから得られる解よりもわずかに高いことが分かるが,その差は10-8 a.u.より小さいので実用上問題にはならない.これに対して,気になるのはVirial比が完全に2にならない点である.また,展開項数Nが7と8の場合には,いくつかの極小点が見つかっており,追いつめきれていない.ただし,どの解もエネルギー値はほぼ同じである.

Table 2. Optimized electronic energies (in a.u.) for 2p state (2P) of H atom by antisymmetric pairs of 1s-type GTF expansions (calculated by GAMERA).
NEnergyVirial ratio ( -V/T )
1-0.11317684690152.0000000179
2-0.12328871222152.0000000123
3-0.12472759999392.0000000103
4-0.12495180957052.0000000094
5-0.12499056775272.000000008
6-0.12499797988062.000000008
7-0.1249995317455(2.00000001)
8-0.124999883442(2.00000002)

Table 3は最適化された中心からの距離と軌道指数である.軌道指数については,Table 1に対応する2p型GTFによる展開の最適値も比較対照のため載せてあるが,比較すると,どの値も極めて近いがごくわずかに1s型GTFの組みのほうが大きくなっていることが確認できる.この原因を考えるために(2)式の右辺を展開して整理すると,この擬似的2pxは,4fxxx, 6hxxxxx, ...をわずかに含むことが確認できる.そして,これらの高次の関数は2p関数よりも外に広がっている.すなわち,これらの影響を緩和するために,軌道指数が2p型GTFによる展開の場合よりもわずかに大きくなっていると予測される.また,この高次の関数の影響がVirial比のずれを引き起こしている可能性があると思われるが,厳密な解析は困難に思われる.

Table 3. Optimized distance (in a.u.) from the center of H atom in antisymmetric pairs of 1s-type GTF expansions and the optimized orbital exponents in comparison with the optimized orbital exponents of standard 2p-type GTF expansions.
Nkantisymmetric pairs of 1s-type GTFs2p-type GTFs
Distance xkExponent akExponent ak
110.0663920.04527435170.0452707393683614
210.0337180.1392872680.13927845725851
20.0783710.03239494120.032392364985181
310.0189670.3370889500.33707270642986
20.0466400.07983970520.079834178056225
30.0867930.02468717580.024685343267111
410.011540.73235200.7323244890509
20.027940.173879380.1738701467249
30.056780.055629330.05562538530627
40.091240.0201667330.02016538916815
510.007781.485691.485642875022
20.017520.3526570.3526432251194
30.03580.11386610.1138598537210
40.06470.04282040.04281745880412
50.0930.017239550.01723853797304
610.005062.86542.86536159637
20.01170.679670.67965980673
30.02310.220190.220179601662
40.04290.0836760.083672161150
50.07100.0349650.034962864559
60.09080.01518430.0151835760134

一方,原点からの距離の最適化が意味しているのは,上述した4fxxx, 6hxxxxx等の高次の関数の寄与の重みを変えて最適化することだと思われるが,展開項数Nが大きくなるに従って,エネルギーの変化に対する影響が小さくなるようであり,N >=7では,初期値に依存して,一義に定まらないように思われる.また,Table 3Nが5と6の場合を比較しても,最小の軌道指数aNに対するxNNが4と5のときの傾向と逆になっていることがわかる.
以上のように,まだいくつか不明確な点もあるが,「2p軌道を1s型GTFで表す」という試みは概ね成功したと考えてよい.

4 H原子の3d軌道を1s型GTFで表す試み

前章と同様に,1s型GTFの組合せで,H原子の3d軌道を表そうと試みた.通常の分子軌道法で使われているデカルト座標型のd関数の中で,3sの寄与を含まないのは,dxydxzdyzの3つであり,ここではdxyを選んだが,次のような4つの1s型GTF重ね合わせで,3 dxy軌道と同じ対称性の関数になることがわかる.

N = 1の場合について最適化を試みたが,残念なことにうまく収束しなかった.その理由は,(3)式の対称性の解がROHF計算の4番目の解となるが,最適化の行程で対称性が崩れるのを防止するためにはGAMERAに最低限でも点群C4を組み込む必要があるためだと判断できる.ただし,前章での議論から,N = 1のときの軌道指数a1Table 1に対応する3d型GTFでの3d軌道の最適化で得られる軌道指数に近いことが予測される.また,中心からの距離に対応するパラメータrkは,最適化に際してそれほど敏感でないことも予想される.そこで,一点計算を繰り返してみた.まず,軌道指数を,N = 1のときの3d型GTFでの3d軌道の最適値(0.01478228224273023)に固定し,rk = 0.3で計算を行ったところ,電子エネルギーとして-0.05173800592638 a.u.が得られた.Table 1との比較により,この時点で3d型GTFのときの最適値より低い解が得られていることが分かる.さらに,rkを小さくするとエネルギーは下降するが,関数間の重なり積分がどんどん大きくなり,計算精度が失われていく現象が見られ,絞り込むことができない.すなわち変分原理を崩してしまう.
以上により,現時点のプログラムではこれ以上の解析は困難であり,3d軌道が1s型GTFである程度表現できそうだという示唆が得られたことで,議論を留めたい.

5 考察

3章と4章の報告によって,厳密性を欠くものの,対称性を考慮すれば,1s型GTFだけでp軌道やd軌道の描像がある程度表現できそうなことが確かめられた.これは,FVMO法を用いれば,1s型の基底関数だけで分子の計算ができることを示唆している.ただし,このような分子軌道プログラムの開発が実用上どの程度可能か,優位性はあるか,問題点はないか,といったことを十分吟味する必要がある.Gaussian Lobe法が定着しなかったという過去の事実は,1s型GTFだけによるプログラムに,少なくとも現行までの時点では優位性がさほどなかったことを示唆している.しかしながら,大規模系への拡張やデータベース技術との連携によって優位性が確保できるのではないかと筆者は考えている.また,こういった点の解析については,若い世代にも大いに期待したい.

本研究が,森和英氏の存在なくして成り立たないことはもちろんであるが,彼がこの世に出現させたFVMO法を発展させた方々,GAMERAに改善を加えた方々の存在ももちろん無視できない.その中でも,全体を管理し,積分計算部分を改良し精度を高める等多くの貢献をした横浜市立大学の立川仁典教授と最適化ルーチンの改善に多大なる貢献をした労働安全衛生総合研究所の大塚輝人氏に深く感謝します.さらに,筆者らの研究を支え続けている高千穂大学情報メディアセンター長の鈴木一成教授に対し,心からお礼を申し上げます.

参考文献

[ 1] S. Huzinaga, J. Chem. Phys., 42, 1293 (1965).
[ 2] K. Mori and C. Ohe, Mem. Kokushikan Univ. CIS, 18, 62 (1997).
[ 3] K. Taneda and K. Mori, Chem. Phys. Letters, 298, 293 (1998).
[ 4] K. Mori, Mem. Kokushikan Univ. CIS, 20, 50 (1997).
[ 5] M. Tachikawa, K. Sasagane, and Y. Osamura, Nonlinear Optics, 26, 43 (2000).
[ 6] K. Sasagane, American Institute of Physics Selected Papers from ICNAAM 2007 and ICCMSE 2007, ed. by T. E. Simos, G. Maroulis, G. Psihoyios, C. Tsitouras, Springer, Dordrecht (2009), pp. 40-43.
[ 7] H. Preuss, Z. Naturforch., 11, 823 (1956).
[ 8] J. L. Whitten, J. Chem. Phys., 39, 349 (1963).
[ 9] J. L. Whitten, J. Chem. Phys., 44, 359 (1966).
[10] J. D. Petke, J. L. Whitten, and A. W. Douglas, J. Chem. Phys., 51, 256 (1969).
[11] E. Clememti, C. C. J. Roothaan, and M. Yoshimine, Phys. Rev., 127, 1618 (1962).
[12] E. Clememti, IBM J. Res. Develop. Suppl., 9, 2 (1965).


Return