はじめに
分子を対象とした量子化学計算プログラムを自作するため、以前の投稿では、核間距離1.4Bohrの水素分子（H2）の各種分子積分の値をハードコードして、SCFのドライバ部分を記述した。
しかし、これだと当然のことながら、何の実用にならない。そこで、今回は、基底関数を導入する。なお、分子積分そのものは、数式が複雑な上に、私自身自分で導出をしたことがないため、参考文献を上げておく。
今回のプログラムは、ここのステップ2においてある。
短縮ガウス型基底関数の導入
量子化学計算で使われる基底には、平面波、数値基底、局在基底などといったものがあるが、ここでは、分子の計算においてはもっとも一般的な、短縮ガウス型基底関数を用いる。これは、ガウス基底を線形結合したものである。
構成要素であるプリミティブガウス基底の表式は以下のようである。
\phi^{PGF} = \frac1{N}x^ly^mz^n\exp(-\alpha|r-R|^2 )
式中のrは任意の座標、Rは核座標、αは軌道係数（軌道の広がり方を表す）で、これは各原子ごとに最適化されている。ここで、注意が必要なのは、ESMLなどで入手できる規定関数のデータは、軌道の広がりを示すα(Orbital Exponent)の値と、このあと出てくる係数Cのみである。しかし、実際には、短縮ガウス関数を組み上げる前に、それぞれのプリミティブガウス関数ごとに、全空間で積分した時に１になるように規格化因子を求める必要がある（上の式のNのこと）。
短縮ガウス型基底関数は、これらの線形結合であるから、
\phi^{CGF} = \frac1{\Omega}\sum_i^{NContraction}C_i \phi^{PGF}(\alpha_i)
となる。
この短縮ガウス関数に関しても規格化が必要となる。要するに、自分自身との重なり積分が１になるように、規格化因子を求める必要があるということだ。
以上の情報を詰め込むために、短縮ガウス基底関数を表すクラスを下記のように作った。初期化する時に、内部で規格化因子を求めるようにしてある。
class PrimitiveGTO:
def __init__(self, lmn, origin, exponent):
self.lmn = lmn
self.exponent = exponent
self.origin = np.zeros(3)
for i in range(3):
self.origin[i] = origin[i]
self.nc = self.calculate_normalization_factor() # nc: Normalization Constant
def calculate_normalization_factor(self):
(l,m,n) = self.lmn
alpha = self.exponent
n = math.sqrt(
(math.pow(2,2*(l+m+n)+1.5)*math.pow(alpha,l+m+n+1.5)) \
/ (factorial2(2*l-1)*factorial2(2*m-1)*factorial2(2*n-1)*math.pow(math.pi,1.5)) )
return n
class ContractedGTO:
def __init__(self, lmn, origin, exponent_list, coefficient_list):
if len(coefficient_list) != len(exponent_list):
raise
self.pgto_list = []
n_contraction = len(exponent_list)
for i in range(n_contraction):
self.pgto_list.append( (coefficient_list[i], PrimitiveGTO(lmn, origin, exponent_list[i])) )
self.nc = self.calculate_normlization_factor() # nc: Normalization Constant
def calculate_normlization_factor(self):
s = 0.
for (coeff1, pgto1) in self.pgto_list:
for (coeff2, pgto2) in self.pgto_list:
s += coeff1 * coeff2 * overlap_PGTO(pgto1, pgto2)
return 1./math.sqrt(s)
def __len__(self):
return len(self.pgto_list)
def __getitem__(self, i):
return self.pgto_list[i]
#積分の評価
積分については短縮ガウスの展開について手短に述べる。プリミティブガウス関数同士の積の評価方法については述べない（というか私自身解析的に手で式を展開してないのでわかってない）。参考文献を読んでほしい。もしくは、プログラムを見てほしい。
HFを解くに当たっては、重なり積分(S)、運動エネルギー積分(T)、核引力積分($V_1,V_2,V_3...$)、電子間相互作用積分(G)の４種類の積分を求める必要があるが、もっとも簡単なのが重なり積分である（それでも結構大変）。これを例にとる。
重なり行列の行列要素の定義は、
S_{ij} = \int \phi_i^{CGF}(x)\phi_j^{CGF}(x) dx
これを、実際に求めるには、プログラムの中で短縮を展開してやる必要がある。プリミティブガウス関数の積分は別として、下記のような関数になるだろう。
def calc_S(basis_array):
dim = len(basis_array)
S = np.zeros( (dim, dim) )
for i in range(dim):
for j in range(dim):
S[i,j] = overlap_CGTO(basis_array[i], basis_array[j])
return S
def overlap_CGTO(cgto1, cgto2):
val = 0.
for i in range(len(cgto1)):
(coeff1, pgto1) = cgto1[i]
for j in range(len(cgto2)):
(coeff2, pgto2) = cgto2[j]
val += coeff1 * coeff2 * overlap_PGTO(pgto1, pgto2)
val = val * cgto1.nc * cgto2.nc
return val
ほぼ同様の形なので、下記に核引力積分、クーロン積分の呼び出し部分を書いておく。
# 核引力積分の行列
def calc_K(basis_array, atom_list):
dim = len(basis_array)
V = np.zeros( (dim,dim) )
for atom in atom_list:
_v = np.zeros( (dim, dim) )
for i in range(dim):
for j in range(dim):
_v[i,j] = nuc_attr_CGTO(basis_array[i], basis_array[j], atom)
V += _v
return V
def calc_G(D,basis_array):
# 電子間相互作用の行列(Gマトリックス）
# Szabo., pp141(3.154)
dim = len(basis_array)
G = np.zeros((dim,dim))
for u in range(dim):
for v in range(dim):
temp = 0.
for p in range(dim):
for q in range(dim):
doubleJ = elec_repulsion_CGTO(basis_array[u],basis_array[v],basis_array[p],basis_array[q])
K = 0.5 * elec_repulsion_CGTO(basis_array[u],basis_array[q],basis_array[p],basis_array[v])
temp += D[p,q] * (doubleJ - K)
G[u,v] = temp
return G
# 核引力積分の行列要素
def nuc_attr_CGTO(cgto1, cgto2, nuclear):
val = 0.
for i in range(len(cgto1)):
(coeff1, pgto1) = cgto1[i]
for j in range(len(cgto2)):
(coeff2, pgto2) = cgto2[j]
val += coeff1 * coeff2 * nuc_attr_PGTO(pgto1, pgto2, nuclear.pos)
val = -val * cgto1.nc * cgto2.nc * nuclear.atomic_number
return val
# 二電子積分
def elec_repulsion_CGTO(cgto1, cgto2, cgto3, cgto4):
val = 0.
for i in range(len(cgto1)):
(coeff1, pgto1) = cgto1[i]
for j in range(len(cgto2)):
(coeff2, pgto2) = cgto2[j]
for k in range(len(cgto3)):
(coeff3, pgto3) = cgto3[k]
for l in range(len(cgto4)):
(coeff4, pgto4) = cgto4[l]
val += coeff1*coeff2*coeff3*coeff4*elec_repulsion_PGTO(pgto1, pgto2, pgto3, pgto4)
val = val * cgto1.nc * cgto2.nc * cgto3.nc * cgto4.nc
return val
要するに、短縮ガウス基底をforループで展開してやり、係数と積分値を掛け合わせて、足し込んでいく。そして最後に、短縮基底ガウス関数自体の規格化因子を掛け合わせる。
多くの積分評価に関する論文で研究がされているのは、overlap_PGTO(1,2)とか、kinetic_PGTO(1,2)のようなプリミティブ関数の展開方法である。そして、本当に難しいのはこの残った部分なのだけど、それに関しては、各種論文やここにおいたプログラム(今回はstep2)を見てほしい。
これらの各積分が評価できるようになれば（文献や、Gaussianで数値をダンピングした結果が一致していれば）、あとは、前回作ったルーチンに渡してやれば自動的にエネルギーとオービタルを求めることができるようになる。
それにしても、やはりpythonで書いたことに加えて、式の対称性を利用した計算回数の最適化を全くしていないため、非常に遅い。水素分子（基底関数はSTO-3GのS軌道２つだけ）でもエネルギーの計算に0.76秒くらいかかっている。
参考文献
全般
Szabo, Attila, and Neil S. Ostlund. "Modern quantum chemistry: introduction to advanced electronic theory." NY: McGraw-Hill (1989).
積分のこと
- Taketa, Hiroshi, Sigeru Huzinaga, and Kiyosi O-ohata. "Gaussian-Expansion Methods for Molecular Integrals." Journal of the Physical Society of Japan 21.11 (1966): 2313-2324.
- Trygve Helgaker, Poul Jorgensen, Jeppe Olsen "Molecular Electronic-Structure Theory"
1は、ガウス関数の展開公式に関する古い論文（Errataがあるとかいう噂を聞いたが真偽不明）。これはGaussianプログラムからは引用されてなかったという昔話を以前聞いた。
2は、波動関数理論(HF, CI, Coupled Clusterなど)の実装に関する大変詳しいテキスト。DFTに関しては全く扱っていない。積分の評価方法も複数書いてあり、参考文献1の代わりにもなる。核引力積分と２電子積分を展開した時に出てくるboys関数の評価方法は、参考文献１にも書いていないのだが（今では簡単には入手できなさそうな文献を引用している）、その評価方法もきちんと言及してある。
その他、もうすこし現代的な積分の評価方法に関しては、こんなページを見つけた。というか、このリンク先のページほうがよく書かれており、本記事の存在意義はこの時点でもうすでにないと言っていい。