1P053 分子間相互作用研究のための原子の高精度 RHF 計算（北教大釧路）1【序】近年のコンピュータの発達に伴ない分子軌道計算が大型分子系へ適用できる様になってきて、比較的小型の酵素について分子軌道計算が行われ報告される様になってきた。大型分子系の分子間相互作用を量子化学計算を用いて研究していく上で発生すると懸念される問題点として、過度の BSSE(Basis Set SuperpositionError) の発生、これを改善するためにより豊富な基底関数系を使用するとこんどは基底関数系の線形従属性の上昇、および、膨大な棄却電子反発積分の累積による過度な引力、などがある。これらを解決する基本的な方策の一つはより高い精度で分子間相互作用の量子化学計算を行うことであり、そのためには、高精度で原子の量子化学計算を行っておくことが必要である。本研究室で開発した、種々の数学関数と誤差関数の高精度ライブラリを使用して高精度 RHF 計算プログラムを開発したので報告する。

2【原子核モデルと核引力積分】

小原繁を使用し、この動径関数に規格化していないガウス関数を用いる：£¤l2Rnl (R) = r exp −ζn r

なお、規格化定数は次式になる。
"rNnl =Uρ (R) =

半径 R0 の球内に核荷電が均一に分布していると考えるモデルで、R0 は原子核半径の二乗平均平方根 hR2 ir
R0 =

5 2hR i3が hR2 i と質量数 A の経験式
p1/3hR2 i = (0.836A

π(8)(2l)!!Z(2ζ)l+1(9)Z exp(−η/2) X (2l)!!R02l+2η(2l − 2i)!!η il

hn0 l|A均 |nli

= −i=0

ZR02l+2−(3η − 2l − 3)Fl+1 (η/2)
2η2l+2ZR0−exp(−η/2)2η(10)になり、ガウス分布モデルでは(2l + 1)!!Z[2ξ(1 + γ)]l+1 (1 + γ)1/2(2l)!!Z−[2(ζn0 + ζn ]l+1 (1 + γ)3/2= −×µ¶jlX(2j + 1)!!
γ(2j)!!

(11)

1+γになる。ここで、ζ 、η 、γ 、および、Fm (T ) は、それぞれ、≡ ζn0 + ζn ,≡ 2(ζn0 + ζn )R02 ,ζ 0 + ζn
≡ n,
Z ξζ
ηγ(2)

(12)

1Fm (T )dt t2m exp(−T t2 )≡0+ 0.570) fmに一致するように決める。
もう一つはガウス分布モデルである：³ ´3/2£¤ξG2ρ (R) = eZ.になり、均一球モデルでは

hn0 l|Aガ |nli

(1)

(R > R0 )

#1/2

2ζnπ

hn0 l|A点 |nli = −j=0

(R ≤ R0 )(7)

点電荷原子核モデルの核引力積分 hn0 l|A点 |nli は原子核のモデルには従来の点電荷モデルの他に有限の大きさを持つ原子核のモデルがいくつかある。[1] 本
研究ではその中の二つを使用する。一つ目は、均一球モデルである：(3eZ ,4πR030.

2(4ζn )l+1(2l + 1)!!

,

exp −ξR

(3)

3.

(4)

原子核の中心の周りに核荷電がガウス分布するモデル
で、拡がりの度合を表す ξ は均一球モデルと同様に hR2 iの期待値
3hR i =2ξ2(5)が質量数 A との経験式 (3) を満足するように決める。
次にこれらのモデルにおける核引力積分の表式を記す。まず、基底関数 ωnlm は実の動径関数 Rnl (r) に球面調和関数 Ylm (θ, ϕ) を乗じたもの
ωnlm = Rnl (r) Ylm (θ, ϕ).

である。

(6)

【RHF 収束法】

指数関数型直交行列 Us により現在の規格直交化軌道
φs を改良したもの φs+1 へ変換することにすると

s+1 = s Us ,Us = exp(As ),exp(As ) = 1 + As +µAs =0xs−xs01¶2A2s + ...,(13).

回転角行列 xs は基底関数の総数 N の自乗のオーダーの要素数を持っているが、以後、これを全要素数を一次元に配置した列ベクトルとして扱う。回転角列ベク
トル xs に関するエネルギーの一次微分は Fock 行列の非対角ブロックの要素に密接に関連していて、例えば、
閉殼系 RHF の場合、

∂E= 4F空, 占 + O(xs )∂xs(14)になる。[2,3] したがって、xs = 0 でのエネルギーの一次微分 gs は¯∂E ¯gs ≡= 4F空, 占(15)¯∂xs xs =0となって厳密に計算できることになる。つまり、回転角
xs と一次微分 gs のいずれの厳密値も SCF 計算の繰り返し毎に得られることになる。
一方、一般に、エネルギー E が未知の点 x0 の近傍において x の二次関数になっている時
1(x − x0 )T B(x − x0 )2

E = E0 +(16)x に関する一次微分 g がg≡∂E= B(x − x0 )∂x(17)x0 = x − B0g004【水素原子の RHF 解】ここに記した方法を取え入れた高精度 RHF 計算プロ
グラム QuCCIE を開発し水素原子の RHF 計算を行った。計算に使用した浮動小数点数は 15 倍精度を有し 10進数にすると 139 桁余りの精度を持っている。基底関数の軌道因子 ζn は 1.25 のベキ表現にしておき、
ζn = 1.25n ,になるから極値点 x0 は−1

この等式から分かるように、定数 Pj 、Qj と列ベクトルxj 、χj 、および、先頭の B−10 （例えば対角行列に選ぶ）だけを保持していれば SCF 計算が可能であり、N 4 オーダーの要素数を持つ B−1s そのものを保持しなくてもよい。これは Almlof らの方法に相当するが、彼らの方法
では２定数 Pj 、Qj の代わりに１つの列ベクトル ∆j を用いた３列ベクトル xj 、χj 、および、∆j を保持する方法になっている。[5,6](18)000になり、２点 x と x での一次微分 g と g には次の関係が成立する(x00 − x0 ) = B−1 (g00 − g0 ).

(19)

これらの関係式 (18) と (19) を SCF 法で使用するパラメータで表記すると
xs+1 = −B−1
厳密 gs ,

(20)

xs = B−1
厳密 (gs+1 − gs )(21)とになる。(21) 式は近似 Hessian 行列 B−1s を改良して−1Bs+1 を得るために使用し、(20) 式は新らたに得られたB−1s+1 を使用して新規の回転角 xs+1 を計算するために使用する。つまり、χs を

s ≡ B−1s (gs+1 − gs )(22)と定義すると、Bs が B厳密 とわずかに異なるので χs はxs とわずかに異なるが、これらのパラメータを使って−1B−1s に補正項を加えて新規の Bs+1 を得る。
−1TTT
B−1s+1 = Bs + Ps xs xs + Qs (xs s + s xs )n = M, M − 1, M − 2, ..., M − 200. (28)

項数を 201 項 (n = M, M − 1, ..., M − 200) に固定してM 値を最適化した。
点電荷モデルとガウス分布モデルの各原子核モデルにおける最適 M 値、RHF エネルギー、および、原子核中心での電荷密度 ρ(0) を表 1 に掲げる。
点電荷モデルの計算は水素原子 (点電荷モデル) の解
析解をガウス関数展開で近似したものになる。計算された RHF エネルギーが厳密値と 28 桁以上も一致していることから、この展開の精度がかなり高いことが判かる。ガウス分布原子核モデルの RHF エネルギーの計算値も同程度の精度を持つものと期待する。
点電荷モデルでの原子核上の電荷密度がガウス分布モデルのものよりも大きな値になっている。これを反映して点電荷モデルの M 値 (M = 179) がガウス分布モデルのもの (M = 134) よりも大きくなり原子核近傍の記述のために多くの基底関数が使用されている。
表 1: 点電荷モデルとガウス分布モデルの各原子核モデ(23)

ルにおける最適 M 値、RHF エネルギー、および、原
ここで、定数 Ps と Qs は次式が成立するように決め、 子核中心での電荷密度 ρ(0)xs = B−1s+1 (gs+1 − gs )次のようになるµT ∆1+ s s

1xTs ∆sxTs ∆s
Qs = − T1 ,xs ∆s∆s ≡ gs+1 − gs .

Ps =(24)¶,(25)

こうして得た B−1s+1 を用いて新規の回転角 xs+1 を計算する。
xs+1 = −B−1s+1 gs .

(26)

ここに記した RHF 収束法は、途中の定義が多少異なるが、擬二次法である Broyden-Fletcher-GoldfarbShanno (BFGS) 法に該当する。[4] なお、B−1s+1 は (23)
式により再帰的に表現されているので xs+1 を次式で計算することができる。
xs+1= −B−1
0 gs−sX£Pj xj xTj gs + Qj (xj Tj gs + j xTj gs )点電荷原子核モデル水素原子
M点 = 179
E点 = −.49999 99999 99999 99999 99999 99942 a.u.
ρ点 (0) = .31830 98859 20809 88259 71148 94837 a.u.
ガウス分布原子核モデル水素原子
Mガ = 134
Eガ = −.49999 99995 29394 49245 93399 50434 a.u.
ρガ (0) = .31829 42994 02107 84496 64495 29498 a.u.

¤j=0(27)

1. L. Visscher and K. G. Dyall, Atomic Data and Nuclear Data Tables, 67, 207-224 (1997).
2. C. C. J. Roothaan and P. S. Bagus, Methods in Computational Physics, Vol. II, Academic Press (1963).
3. G. L. Malli and J. P. Olive, Doctor thesis, Universityof Chicago.
4. R. Fletcher, Practical methods of optimization 2nded., Wiley, 1987.
5. T. H. Fischer, and J. E. Almlof, Journal of Physical
Chemistry, Vol. 96, 1992, p. 9768.
6. G. Chaban, M. W. Schmidt, and M. S. Gordon, Theoretical Chemistry Accounts, Vol. 97, 1997, pp. 88-95.