縮約ガウス関数にもとづく分子積分の計算法

本田 宏明, 小原 繁


Return

1 はじめに

非経験的量子化学計算に基づく分子系電子構造の理論研究を行なう際には、分子積分の高速計算が依然として重要かつ大きな問題である。とりわけ、遷移金属化合物や生体内重要分子の様な大型分子系においては分子積分の計算が最も労力を必要とする計算ステップとなる。最近求められている高い精度の量子化学計算を行なうには、より大きな基底関数系を使用することが必要であり、それには、原子軌道の展開項数を増すばかりでなく、角運動量の高い軌道も加える必要がある。この様に大きな基底関数系を用いる量子化学計算を実行するには分子積分の求値、とりわけ、電子反発積分(ERIs)の求値を高速に完了することが不可欠であり、これを可能にすることが大型分子系の広範な理論研究を促進することになる。
Boys らがガウス関数を基底関数とし分子積分を計算して以来、非経験的量子化学計算においてガウス関数が基底関数に使用されて来ており、分子積分の表式と計算法について多くの研究が行われてきた[1 - 21]。M.Honda, Sato, Obara らにより、Hypergaussian と名付づけられた Fourier核を含む一般化されたガウス関数に対し、原始関数を基底とする分子積分の漸化表式が提案された。この表式により、通常の電子反発積分のみならず、基底関数にexp(ikr) が含まれる種々の1、 2電子積分について導関数や2次導関数の積分表式の導出まで適用可能となった。
また最近、Ten-noにより縮約ガウス関数についての電子反発積分やその1次微分に対し漸化表式が提案されている[17]。これに対し、Hypergaussian を縮約関数形に拡張することによる縮約ガウス関数に対する分子積分の一般漸化表式が H.Honda, Yamaki, Obara らにより発表されており[2]、この手法は1986 に報告された Obara による種々の分子積分の表式[15]を一般化し、更に縮約基底にまで拡張した新しい Obara による分子積分の方法であるといえる。この中で提案されている縮約ガウス関数を基底関数に用いた分子積分の一般表式には次の3つの特長がある;
  1. 種々の分子積分に適用できる表式である。
  2. 効率良い算法として知られている漸化計算法を使用できる。
  3. 基底関数の縮約を漸化計算の開始積分(初期積分)についてのみ行なうだけで良い。
これらの特長から、種々の分子積分を共通の漸化計算手法を使用し計算することが可能となる。
次節では、一般の縮約ガウス関数についての漸化表式の導出を概観する。3節では具体的に電子反発積分の表式と計算アルゴリズムを概説した後、具体的な積分について詳述し、まとめとする。

2 縮約ガウス関数に対する一般漸化表式

縮約ガウス関数に対する分子積分の一般漸化表式をの導出のためのステップ(アイデア)は4つあり、以下 2.1~2.4 節に詳細する。

2. 1 Hypergaussian の利用

電子反発積分 <pmpn,ss>等を表現するために、一般的な原始基底関数と演算子を持つ2電子積分 <FaFb|F12|FcFd>に対し、演算子や関数を具体的な表式へと簡約することを可能とする簡約演算子 を作用させる。

理想的には <FaFb|F12|FcFd> を十分一般的に定義することで種々の積分表式が導出できるはずであり、そのための関数 F として以下の Hypergaussian を使用する。

2. 1. 1 Hypergaussian の定義

簡約演算子の作用により、種々の量子化学計算の場に現われる具体的な基底関数や演算子へと変換可能なHypergaussian FH

と、定義する[12]。ここで R = (Rx, Ry, Rz) は関数中心、r = (rx, ry, rz) は電子座標であり、f(r - R;n, z) については

と、軌道指数 z, 非負数整数組 n を持つガウス関数を表す。式(2)における (l)は関数中心 R に対する微分演算子を表す。

ここで非負整数組 l = (lx, ly, lz) はx, y, z それぞれの方向について何回微分が行われるかを示す。
式(2)における k = (kx, ky, kz)は3つの実数組を表しており、電子座標 r とは独立であると仮定する。

2. 1. 2 Hypergaussian の簡約

Hypergaussian FH に対し、簡約演算子 を作用させることでガウス関数と電子反発演算子を実際に導出することが可能である。
原始ガウス関数については

とすることで得られる。px へ簡約する場合は更に

とすれば良い。
更に電子反発演算子についても、

と、Hypergaussian とその簡約操作により表現できる。

2. 2 漸化関係の導出

前節の議論から、一般の2電子積分は5つの Hypergaussianの積を被積分関数として持つ2電子積分により表現可能といえる。ここで、

と被積分関数を とする。FHz が 0 でない限り、関数中の exp( - z (r - R)2) 部分の寄与により

となるが、 についても同様である。そのため

なる性質を持つ。ここで k = 1, 2 は電子番号、m = x, y, z である。この k = 1, 2 の2種類の微分・積分により2つの等式が得られるが、連立方程式とみたて解くことにより nl等に対する漸化表式を導くことが可能となる[12]。また、この一般の2電子積分についての漸化表式について、式(1)と同様 による簡約を行なうことにより、目的となる個別の積分の具体的漸化表式を得ることが可能となる。

2. 3 補助縮約被積分関数$\CG^{AC}$による一般の縮約分子積分

上記の議論を縮約ガウス関数の分子積分に拡張するためには被積分関数を補助縮約被積分関数へと以下のように拡張する必要がある。

1,2 は電子番号を表わしており、rzQ については以下で定義されている。

については以下で定義される。

ここで、ab の関数に対する ijについての和が縮約和に対応しており、ua ub は縮約係数となる。 の定義についても同様である。 により、一般の縮約ガウス関数に対する一般の2電子積分は縮約 Hypergaussian FCH により

と表わすことができる。またFaCHについては nala, ga, qa により特徴付けられていることをふまえ、nanb等を電子番号でまとめることにより、以下のように記述する。

この式は一般の原始ガウス関数の積分に係数をかけて縮約和をとった形となっている。このことから、原始ガウス関数の積分の漸化表式に対して縮約和をとる事で、一般の縮約ガウス関数の分子積分の漸化関係を代数的に求めることができ、以下の最終表式を得る。

ここで、Nm(n) = nm ,

であり、Xj'm(k)Yj'm(k) は以下の式をまとめたものである。

この表式に対して簡約演算子 を適用することにより種々の縮約ガウス関数に対する種々の演算子の2電子積分の漸化式を得ることが可能となる。通常の縮約ガウス関数の電子反発積分については3節で詳しく述べる。

2. 4 1電子積分への簡約

縮約ガウス関数についての一般の1電子積分は2電子積分に関する表式(4)を以下の Dirac のデルタ関数

により電子座標を簡約することで得られる。さらには簡約演算子 を適用することにより具体的な1電子積分漸化表式を得ることが出来る。Tables 1, - 3p 関数までの重なり積分や運動エネルギー積分、運動エネルギー積分、核引力積分を示す。

Table 1. Recurrence expressions for overlap integrals over s and p contracted Cartesian Gaussian functions

Table 2. Recurrence expressions for kinetic energy integrals over s and p contracted Cartesian Gaussian functions

Table 3. Recurrence expressions for nuclear attractive integrals over s and p contracted Cartesian Gaussian functions

3 電子反発積分の表式とアルゴリズム

3. 1 一般表式

前節までの議論によって導かれた電子反発積分 <abcd>の表式をTable 4(a)に掲げる。この積分の中のパラメータaは縮約軌道 F( r1, a) の軌道量子数を表わす[2]。電子反発積分はこの軌道量子数の増減に関する漸化関係を満足している。この漸化関係は通常垂直関係と呼ばれている。この垂直関係を数式で表現したものをTable 5(a)に掲載する。ここで、補助電子反発積分<abcd | m, ga gb Gab gc gd Gcd >[Table 4(b)]を使用して表記した[2]。

Table 4. The electron repulsion integrals(ERIs), the auxiliary ERIs,and the auxiliary initial ERIs

Table 5. The vertical and horizontal relations of ERIs

漸化計算を開始するには初期積分が必要である。この表式をTable 4(c)に掲載した。なお、電子反発積分は水平関係[Table 5(b)]と呼ばれる関係も満足している。この関係を使用すると電子反発積分<abcd>を計算する際に、軌道量子数 bd が零になった電子反発積分 <a0b0>を計算すれば良いことになる。これらの関係を利用して効率的な計算を行な方法を電子反発積分 <ppss> ( = <1a1b0c0d>)の場合について詳述する。d関数以上の角運動量を含む積分についても、水平関係と垂直関係の漸化計算を経ることで同様な計算が可能である。

3. 2 <ppss>の計算

電子反発積分<ppss>を計算する上で前述の水平関係と垂直関係の利用方法と初期積分の求値法を記す。

3. 2. 1 水平関係

電子反発積分 <ppss>( = <1100>において軌道量子数d0であるが、b1であり零にはなっていない。Table 5(a)に掲載した垂直関係はbdがともに零になっている時の表式であるので、この表式を使用するためにはbを零にする必要がある。このために水平関係[Table 5(b)]を使用する;

右辺の積分<dmnsss>と<pmsss>の計算には垂直関係を使用すればよい。

3. 2. 2 垂直関係

前記の積分<dmnsss>と<pmsss>を垂直関係[Table 5(a)]により軌道量子数の値を減らして行き、右辺に現われる総ての軌道量子数が零(つまり、初期積分)になるまでこの操作を繰り返す。例えば、<dmnsss>では

になって、初期積分以外の<pmsss|0 011 000>などが右辺に現われるので、これらの積分をさらに垂直関係を利用して初期積分になるまで垂直関係を利用する;

一方、(11)式の右辺第二項の積分<pmsss>では

になり、右辺の積分は総て初期積分になる。

3. 2. 3 初期積分の計算

総ての軌道量子数を零にしたならば、目的積分の計算に必要な総ての初期積分をリストアップしたことになる。<ppss>の場合の初期積分は15種類ある[Table 6(a)]。これらをTable 4(c)に掲げた四重の和を含む表式に基づいて計算する。この時、下記の工夫により演算数を減らす;パラメータgbの値(積分表記中の縦線|の右側第三番目の数値)は15種類の初期積分において、0,1,2の三通りの値になるが、どの値になるかは最外側の総和の計算の時にzpbの何乗を乗じるかにより決まるのであり、第三内側の総和計算( )においてはgb = 0と見なしておいてもよい。こうすると、第三内側までの総和計算( )においてはTable 6(b)に掲載した13種類の初期積分を計算するだけでよいことになる。 同様に考えると、第二内側までの総和計算( )において15種類の初期積分の内、gb = ga = Gab = 0にすることにより得られるTable 6(c)に掲載した6種類の初期積分を計算すればよく、さらに、最内側の総和計算( )においてはgb = ga = Gab = gc = 0にすることにより得られるTable 6(d)に掲載した6種類の初期積分を計算すればよいことになる。

Table 6. All the necessary initial ERIs for the calculation of<ppss>

3. 2. 4 計算手順

実際の計算は上記と逆順に行なう。Table 6(d)の6種類の初期積分の求値から計算を開始してTable 6(a)の15種類の初期積分までを計算し、次に、垂直関係[(14, 12, 13)式]により(11)式の右辺に現われる積分<dmnsss>と<pmsss>を求値して、最後に、水平関係[(11)式]により目的積分である<pmpnss>を求める。

3. 3 ループ構造

電子反発積分全体を計算するためのループ構造をFigure 1に描く。基底関数の中から4つの縮約関数を選び出し(a, b, c, dに関する4重ループ)、初期積分を求値した後、垂直関係と水平関係により電子反発積分の値を得る。初期積分の求値では原始関数についての総和計算(pb, pa, pd, pcの4重ループ)を行なう。この総和計算において、Kpapb、(zpa + zpb)-1のように関数abの対だけで値が決まる定数とKpcpd、(zpc + zpd)-1のように関数cdの対だけで値が決まる定数を計算する[Table 4(c)参照]。この計算ステップをFigure 1では「(C)」と表現してある。これらの対情報計算には指数計算や逆数計算を含まれているが、Figure 1のループ構造では、c d対情報をpapbが変わる毎に再計算することになり、余分に計算時間を消費する。
これを改善する方法の一つが、対情報をFigure 1の「(B)」で計算・保持し、pb, pa, pd, pc の4重ループ内では保持した値を参照するだけにする方法である。保持するデータの量は(10 × nc nd × 8)バイト程度である。ncndは通常3程度であるから、1KB以下であり保存領域としては充分に小さいものである。


Figure 1. The loop structure for the calculation of ERIs

なお、この方法も関数abが変わる毎に再計算することになる(余分に消費する計算時間は微少)。完全に再計算しない方法はFigure 1の「(A)」で対情報を計算・保持するものである。しかし、この方法では保持データの総量が約 (10 × nc nd × 8) × N2バイトになる。基底関数の総数Nが通常は5000程度であることを考えると、積分カットオフを考慮してもデータ総量は 200MB 以上となり、現在の通常の計算機においてもかなりのメモリ資源を消費してしまうことになる。また、Figure 1 の "Recursive caluculation" が pb, pa, pd, pc 縮約ループの外にあることからも分る通り、今回開発した計算法では縮約計算を初期積分計算の際に終えてしまい、漸化計算は縮約数に依存せず高速に実行可能であるという利点がある。

4 計算時間の比較

一酸化炭素分子の電子反発積分の計算時間を4種類の基底関数を使用して測定しTable 7にまとめた。基底関数 DZV, TZV については p 関数まで、 DZVP, TZVP についてはd関数までの積分を含んでいる。「GAMESS」は量子化学計算プログラムGAMESS[13](2005-06-27版)の Rys多項式と呼ばれる区分求積法からの電子反発積分計算時間であり、Figure 1の「(B)」において対情報を計算している。GAMESSプログラムについては Intel Fortran Compiler9.0(-O3 最適化)により、本プログラムについては Intel C++ Compiler9.0(-O3 最適化) によりコンパイルをし、Pentium 4(2.4GHz, 2GB Mem, 512KB L2) コンピュータ上にて測定を行なった。GAMESSよりも1.5〜2.0倍程度の速さとなっており、4種類の基底関数のいずれにおいても量子化学計算プログラムとして広く使用されているGAMESSよりも速く、特に大きなメモリ領域を必要としない方法になっている。

Table 7. Execution time and ratio of the evaluation of ERIs of CO molecule (10-2 seconds)

5 まとめ

Honda, Yamaki, Obaraが提案した、縮約ガウス関数を基底関数に関する分子積分の一般漸化表式を概観し、電子反発積分を効率良く高速計算するアルゴリズムを提案した。また、電子反発積分<ppss>について、このアルゴリズムを詳細に記した。一酸化炭素分子の電子反発積分についてGAMESSと計算時間を比較し、1.5~2.0倍の高速化が達成されることを確認した。更には対情報を計算する位置に関して、電子反発積分計算のループ構造を明示し、シェル4重ループの内側すぐの位置の場合に余分な計算時間が少なく、また、必要メモリ量も少ないことを明らかにした。
また、このアルゴリズムは電子反発積分以外の分子積分にも適用できるため、多様な理論化学研究の場に現われる種々の分子積分の高速求値法としても役立つものと期待される。

本研究の一部は科学技術振興調整費総合研究「科学技術計算専用ロジック組込み型プラットフォーム・アーキテクチャに関する研究」(代表 村上和彰 九州大学教授)によるものである。

参考文献

[ 1] H. F. King, M. Dupuis, Numerical integration using Rys polynomials, J. Comput. Phys., 21, 144-165 (1976).
[ 2] H. Honda, T. Yamaki, S. Obara, Molecular integrals evaluated over contracted gaussian functions by using auxiliary contracted hyper-gaussian functions, J. Chem. Phys., 117, 1457-1469 (2002).
[ 3] H. Taketa, S. Huzinaga, K. O-ohara, Gaussian-expansion methods for molecular integrals, J. Phys. Soc. Jpn., 21, 2313-2324 (1966).
[ 4] J. A. Pople, W. J. Hehre, Computation of electron repulsion integrals involving contracted gaussian basis functions, J. Comput. Phys., 27, 161-168 (1978).
[ 5] J. Rys, M. Dupuis, H. F. King, Computation of electron repulsion integrals using the rys quadrature method, J. Comput. Chem., 4, 154-157 (1983).
[ 6] K. Ishida, Rigorous algorithm for the electron repulsion integral over the generally contracted solid harmonic gaussian-type orbitals, J. Chem. Phys., 113, 7818-7829 (2000).
[ 7] K. Ishida, Accompanying coordinate expansion formulas derived with the solid harmonic gradient, J. Comput. Chem., 23, 378-393 (2002).
[ 8] L. E. McMurchie, E. R. Davidson, One- and two-electron integrals over cartesian gaussian functions, J. Comput. Phys., 26, 218-231 (1978).
[ 9] M. Dupuis, A. Marquez, The rys quadrature revisited: A novel formulation for the efficient computation of electron repulsion integrals over gaussian functions, J. Chem. Phys., 114, 2067-2078 (2001).
[10] M. Dupuis, J. Rys, H. F. King, Evaluation of molecular integrals over gaussian basis functions, J. Chem. Phys., 65, 111-116 (1976).
[11] M. Head-Gordon, J. A. Pople, A method for two-electron gaussian integral and integral derivative evaluation using recurrence relations, J. Chem. Phys., 89, 5777-5786 (1988).
[12] M. Honda, K. Sato, S. Obara, Formulation of molecular integrals over gaussian functions treatable by both the laplace and fourier transforms of spatial operators by using derivative of fourier-kernel multiplied gaussian, J. Chem. Phys., 94, 3790-3804 (1991).
[13] S. F. Boys, Electronic wave functions. i. a general method of calculation for the stationary states of any molecular system, Proc. Roy. Soc. London, A200, 542-554 (1950).
[14] S. Obara, A new recursion formula to evaluate molecular integras over gaussian functions, Proceedings of the 6th Joint EPS-APS International Conference on Physics Computing, 1994, 495-502.
[15] S. Obara, A. Saika, Efficient recursive computation of molecular integrals over cartesian gaussian functions, J. Chem. Phys., 84, 3963-3974 (1986).
[16] S. Obara, A. Saika, General recurrence formulas for molecular integrals over cartesian gaussian functions, J. Chem. Phys., 89, 1540-1559 (1988).
[17] S. Ten-no, An efficient algorithm for electron repulsion integrals over contracted gaussian-type functions, Chem. Phys. Lett., 211, 259-264 (1993).
[18] T. Yamaki, S. Obara, An efficient method of computating molecular integrals over contracted gaussian basis functions, Regional meeting of Hokkaido, Chemical Society of Japan, 2001, E16.
[19] V. R. Saunders, Science research council, Technical Report No. DL/SRF, Daresbury Laboratory, England, Publ., p.157.
[20] V. R. Saunders, Computational Techniques in Quantum Chemistry and Molecular Physics, Reidel, Dordrecht (1975).
[21] V. R. Saunders, Molecular Integrals for Gaussian Type Functions, Methods in Computational Physics, Reidel, Dordrecht, 1983.
[22] M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, J. A. Montgomery, General atomic and molecular electronic structure system, J. Comput. Chem., 14, 1347-1363 (1993).


Return