# 分子軌道法 ー 計算代数幾何学からのアプローチ
背景
文献[KIKUCHI2013]においてAkihito Kikuchiは、独特の分子軌道法を展開している。
この方法の特徴は次の通り。
- 一電子・二電子積分を原子間距離の解析関数を以て表現する。
- 波動関数（LCAO係数）も同様に未定変数とする。
- エネルギー汎関数(H)、波動関数の規格直交条件(S)は、原子間距離およびLCAO係数の解析関数を以て表現する
- 波動関数の規格直交条件(S)のもとでエネルギー汎関数(H)を最小化する。ラグランジュ未定定数(e)が一電子軌道エネルギーとなる。 目的関数は(H)-(e)(S)である。
- 目的関数を、原子間距離変数Rの特定の値の周囲でテイラー展開する。よって目的関数は他変数の多項式である。
- 目的関数Fをその引数（LCAO係数、核間距離R）で偏微分することで、その根が目的関数の極小を与える連立多項式を作る。
- この連立多項式を、「多項式環のイデアル」とみなし、グレブナ基底・準素イデアル分解と呼ばれる計算代数幾何学の手法を適用し、数値的に「解きやすく」、また、電子構造がおのずと開示されるような、複数の連立多項式に分解する。分解後の小さい多項式系を数値的に解く。
[KIKUCHI2013]では一電子・二電子積分はSTO(SLater-Type Orbital)を使って計算している。しかし、STOを用いて二電子積分の解析式をどう計算するかは、多原子系では今なお解決できていない問題である。
そこで、この記事では、STOの代わりに、量子化学計算の標準的な方法であるGTO(Gaussian-Type Orbital)を使って[KIKUCHI2013]の方針を追うことを目的にする。
[KIKUCHI2013] @article{KIKUCHI2013, author="Kikuchi,Akihito", title="An approach to first principles electronic structure computation by symbolic-numeric computation", journal="QScience Connect",
volume="2013:14",year=2013, note="http://dx.doi.org/10.5339/connect.2013.14"}
モデル
HeH+分子を考える。Szabo and Ostlund [szabo2012modern]では、この分子の計算を目的とした古風なFortranプログラムを示している。
本記事の研究では、このFortranプログラムを出発点とする。まずこのプログラムをpythonの数値計算プログラムに書き換え、ついで、pythonの現代的な記号処理ライブラリ(sympy)を用いて、関連するすべての電子積分の解析式・テーラー展開を実行し、多項式型のエネルギー汎関数および波動関数の正規直交条件を構成した。ついで、記号処理による偏微分によって、解くべき連立多項式を作成した。この連立多項式を、計算代数システム「Singular」で処理することで数値解を求める。
[szabo2012modern]@book{szabo2012modern,
title={Modern quantum chemistry: introduction to advanced electronic structure theory},
author={Szabo, Attila and Ostlund, Neil S},
year={2012},
publisher={Courier Corporation}
}
#!/usr/bin/env python
import sympy as sympy
from sympy import symbols, Function,expand,core,sqrt
import numpy as np
from sympy import symbols, Function,expand,sqrt
import numpy as np
import copy
S12,T11,T12,T22,V11A,V12A,V22A,V11B,V12B,V22B,V1111,V2111,V2121,V2211,V2221,V2222=[0 for _ in range(16)]
S=[[0 for _ in range(2)] for _ in range(2)]
X=[[0 for _ in range(2)] for _ in range(2)]
NM=[[0 for _ in range(2)] for _ in range(2)]
XT=[[0 for _ in range(2)] for _ in range(2)]
H=[[0 for _ in range(2)] for _ in range(2)]
F=[[0 for _ in range(2)] for _ in range(2)]
G=[[0 for _ in range(2)] for _ in range(2)]
C=[[0 for _ in range(2)] for _ in range(2)]
FPRIME=[[0 for _ in range(2)] for _ in range(2)]
CPRIME=[[0 for _ in range(2)] for _ in range(2)]
P=[[0 for _ in range(2)] for _ in range(2)]
ODLPP=[[0 for _ in range(2)] for _ in range(2)]
E=[[0 for _ in range(2)] for _ in range(2)]
TT=[[[[0 for _ in range(2)] for _ in range(2)] for _ in range(2)] for _ in range(2)]
IOP=2
#
# CHOOSE STO-NG
#
N=3
R=1.4632
#
# Gaussian exponents:
# He:ZETA1; H;ZETA2
#
ZETA1=2.0925
ZETA2=1.24
#
#
#
ZA=2.0
ZB=1.0
#
# ATOMIC DISTANCE
#
R=symbols("R",positive=True)
# WAVEFUNCTION (x,y)
x, y = symbols("x y")
def F0(ARG):
#
# F0 FUNCTION
#
PI = np.pi
if type(ARG)==float and ARG < 1.0e-6:
return 1 -ARG/3.
if type(ARG)==sympy.core.numbers.Zero and ARG < 1.0e-6:
return 1 -ARG/3.
else:
#print("F0:ARG",ARG)
return sqrt(PI/ARG)*sympy.erf(sqrt(ARG))/2
def getS(A,B,RAB2):
#
# Calculates the overlap matrix
#
PI=np.pi
#print(A,B,RAB2)
return (PI/(A+B))**1.5*sympy.exp(-A*B*RAB2/(A+B))
def getT(A,B,RAB2):
#
# Calculates the kinetic energy.
#
PI=np.pi
return A*B/(A+B)*(3.0-2.0*A*B*RAB2/(A+B))*(PI/(A+B))**1.5*sympy.exp(-A*B*RAB2/(A+B))
def V(A,B,RAB2,RCP2,ZC):
#
# CALCULATES UN-NORMALIZED NUCLEAR ATTRACTION INTEGRALS
#
PI=np.pi
V=2.0*PI/(A+B)*F0((A+B)*RCP2)*sympy.exp(-A*B*RAB2/(A+B))
return -V*ZC
def TWOE(A,B,C,D,RAB2,RCD2,RPQ2):
#
# CALCULATES TWO-ELECTRON INTEGRALS FOR UN-NORMALIZED PRIMITIVES
# A,B,C,D ARE THE EXPONENTS ALPHA, BETA, ETC.
# RAB2 EQUALS SQUARED DISTANCE BETWEEN CENTER A AND CENTER B, ETC.
#
PI=np.pi
return 2.0*(PI**2.5)/((A+B)*(C+D)*np.sqrt(A+B+C+D))*F0((A+B)*(C+D)*RPQ2/(A+B+C+D))*sympy.exp(-A*B*RAB2/(A+B)-C*D*RCD2/(C+D))
def INTGRL_SYMBOL(IOP,N,R,ZETA1,ZETA2,ZA,ZB):
#
# CALCULATES ALL THE BASIC INTEGRALS NEEDED FOR SCF CALCULATION
#
global S12,T11,T12,T22,V11A,V12A,V22A,V11B,V12B,V22B,V1111,V2111,V2121,V2211,V2221,V2222
COEF=[[1.0,0.0,0.0],[0.678914,0.430129,0.0],[0.444635,0.535328,0.154329]]
EXPON=[[0.270950,0.0,0.0],[0.151623,0.851819,0.0],[0.109818,0.405771,2.22766]]
R2=R*R
N=3
A1=[0 for I in range(N)]
D1=[0 for I in range(N)]
A2=[0 for I in range(N)]
D2=[0 for I in range(N)]
PI=np.pi
for i in range(N):
#print(i,EXPON[N-1][i],COEF[N-1][i],ZETA1,ZETA2)
A1[i]=EXPON[N-1][i]*ZETA1**2
D1[i]=COEF[N-1][i]*((2.0*A1[i]/PI)**0.75)
A2[i]=EXPON[N-1][i]*(ZETA2**2)
D2[i]=COEF[N-1][i]*((2.0*A2[i]/PI)**0.75)
#print("A1",A1,"D1",D1,"A2",A2,"D2",D2)
S12=0.0
T11=0.0
T12=0.0
T22=0.0
V11A=0.0
V12A=0.0
V22A=0.0
V11B=0.0
V12B=0.0
V22B=0.0
V1111=0.0
V2111=0.0
V2121=0.0
V2211=0.0
V2221=0.0
V2222=0.0
for I in range(N):
for J in range(N):
RAP=A2[J]*R/(A1[I]+A2[J])
RAP2=RAP**2
RBP2=(R-RAP)**2
S12+=getS(A1[I],A2[J],R2)*D1[I]*D2[J]
T11=T11+getT(A1[I],A1[J],0.0)*D1[I]*D1[J]
#print(I,J,getT(A1[I],A1[J],0.0),D1[I],D1[J],T11)
T12=T12+getT(A1[I],A2[J],R2)*D1[I]*D2[J]
T22=T22+getT(A2[I],A2[J],0.0)*D2[I]*D2[J]
V11A=V11A+V(A1[I],A1[J],0.0,0.0,ZA)*D1[I]*D1[J]
V12A=V12A+V(A1[I],A2[J],R2,RAP2,ZA)*D1[I]*D2[J]
V22A=V22A+V(A2[I],A2[J],0.0,R2,ZA)*D2[I]*D2[J]
V11B=V11B+V(A1[I],A1[J],0.0,R2,ZB)*D1[I]*D1[J]
V12B=V12B+V(A1[I],A2[J],R2,RBP2,ZB)*D1[I]*D2[J]
V22B=V22B+V(A2[I],A2[J],0.0,0.0,ZB)*D2[I]*D2[J]
for I in range(N):
for J in range(N):
for K in range(N):
for L in range(N):
RAP=A2[I]*R/(A2[I]+A1[J])
RBP=R-RAP
RAQ=A2[K]*R/(A2[K]+A1[L])
RBQ=R-RAQ
RPQ=RAP-RAQ
RAP2=RAP*RAP
RBP2=RBP*RBP
RAQ2=RAQ*RAQ
RBQ2=RBQ*RBQ
RPQ2=RPQ*RPQ
V1111=V1111+TWOE(A1[I],A1[J],A1[K],A1[L],0.0,0.0,0.0)*D1[I]*D1[J]*D1[K]*D1[L]
V2111=V2111+TWOE(A2[I],A1[J],A1[K],A1[L],R2,0.0,RAP2)*D2[I]*D1[J]*D1[K]*D1[L]
V2121=V2121+TWOE(A2[I],A1[J],A2[K],A1[L],R2,R2,RPQ2)*D2[I]*D1[J]*D2[K]*D1[L]
V2211=V2211+TWOE(A2[I],A2[J],A1[K],A1[L],0.0,0.0,R2)*D2[I]*D2[J]*D1[K]*D1[L]
V2221=V2221+TWOE(A2[I],A2[J],A2[K],A1[L],0.0,R2,RBQ2)*D2[I]*D2[J]*D2[K]*D1[L]
V2222=V2222+TWOE(A2[I],A2[J],A2[K],A2[L],0.0,0.0,0.0)*D2[I]*D2[J]*D2[K]*D2[L]
def COLLECT(IOP,N,R,ZETA1,ZETA2,ZA,ZB):
#
# THIS TAKES THE BASIC INTEGRALS FROM COMMON AND ASSEMBLES THE
# RELEVENT MATRICES, THAT IS S,H,X,XT, AND TWO-ELECTRON INTEGRALS
#
global H,S,X,TT,XT
H[0][0]=T11+V11A+V11B
H[0][1]=T12+V12A+V12B
H[1][0]=H[0][1]
H[1][1]=T22+V22A+V22B
S[0][0]=1.
S[0][1]=S12
S[1][0]=S[0][1]
S[1][1]=1.
TT[0][0][0][0]=V1111
TT[1][0][0][0]=V2111
TT[0][1][0][0]=V2111
TT[0][0][1][0]=V2111
TT[0][0][0][1]=V2111
TT[1][0][1][0]=V2121
TT[0][1][1][0]=V2121
TT[1][0][0][1]=V2121
TT[0][1][0][1]=V2121
TT[1][1][0][0]=V2211
TT[0][0][1][1]=V2211
TT[1][1][1][0]=V2221
TT[1][1][0][1]=V2221
TT[1][0][1][1]=V2221
TT[0][1][1][1]=V2221
TT[1][1][1][1]=V2222
return
def MATOUT(A,IM,IN,M,N,LABEL):
print(LABEL)
print(A)
return
def MULT(A,B,C,IM,M):
for i in range(2):
for j in range(2):
C[i][j]=0
for k in range(2):
C[i][j]+=A[i][k]*B[k][j]
return
def FORMG():
#
# Calculates G matrix from density matrix and two-electron integrals
#
global G,P,TT
for I in range(2):
for J in range(2):
G[I][J]=0
for K in range(2):
for L in range(2):
#print(I,J,K,L,P[K][L],TT[I][J][K][L],TT[I][L][J][K],G[I][J])
G[I][J]+=P[K][L]*(TT[I][J][K][L]-0.5*TT[I][L][J][K])
#print("G",G)
return
def SCF_SYMBOL(IOP,N,R,ZETA1,ZETA2,ZA,ZB):
#
# PREPARES THE ANALYTIC FORMULA OF THE TOTAL ENERGY.
#
global G,C,FPRIME,CPRIME
PI=np.pi
CRIT=1.0e-4
MAXIT=25
ITER=0
for I in range(2):
for J in range(2):
P[I][J]=0.
P[0][0]=x*x*2
P[0][1]=x*y*2
P[1][0]=x*y*2
P[1][1]=y*y*2
MAXIT=1
for ITER in range(MAXIT):
FORMG()
for i in range(2):
for j in range(2):
F[i][j]=H[i][j]+G[i][j]
EN=0
for i in range(2):
for j in range(2):
EN+=0.5*P[i][j]*(H[i][j]+F[i][j])
ENT=EN+ZA*ZB/R
return EN,ENT
#
# COMPUTES THE ELECTRON ENERGY (EN) and THE TOTAL ENERGY(ENT)
#
INTGRL_SYMBOL(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
COLLECT(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
EN,ENT=SCF_SYMBOL(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
from sympy import Symbol, cos, series
#
print("CHECK THE TOTAL ENERGY: BY THIS SUBSTITUTION WE'LL GET THE VALUE CLOSE TO -2.86...")
#
print(ENT.subs([(x,-0.80191751),(y,-0.33680049),(R,1.463)]))
# TAYLOR EXPANSION OF THE TOTAL ENERGY AROUNR R0=1.5
ENTS=series(expand(ENT),R,1.5)
ENTS2=expand(ENTS.removeO())
print("CHECK THE TOTAL ENERGY: BY THIS SUBSTITUTION WE'LL GET THE VALUE CLOSE TO -2.86...")
print(sympy.N(ENTS2.subs([(x,-0.80191751),(y,-0.33680049),(R,1.463)])))
#
# THE NORMALIZARION CONDITION OF THE WAVE FUNCTION
#
W=[0,0]
W[0]=S[0][0]*x+S[0][1]*y
W[1]=S[1][0]*x+S[1][1]*y
NORMC=expand(x*W[0]+y*W[1])
#
# THE NORMALIZATION CONDITION AFTER TAYLOR EXPANSION
#
NORMCS=expand(series(NORMC,R,1.5).removeO())
print("CHECK THE NOTMALIZATION: BY THIS SUBSTITUTION WE'LL GET THE VALUE CLOSE TO 1...")
print(NORMCS.subs([(x,-0.80191751),(y,-0.33680049),(R,1.463)]))
#
# TURN THE POLYNOMIAL ENERGY FUNCTIONAL INTO A SYMPY-POLYNOMIAL
#
ENTS3=sympy.poly(sympy.N(ENTS2))
#
# GET THE OBJECTIVE FUNTION WITH INTEGER COEFFICIENTS
#
def getINT(x,N):
#
# ROWND DOWN x*N AFTER THE DECIMAL POINT and get AN INTEGER M
# THEN X is approximated by M/N.
#
return (int(np.floor(x*N)))
AT=ENTS3.terms()
getF=0
for tm in AT:
p0,p1,p2=tm[0]
cf=tm[1]
#print(p0,p1,p2,cf)
getF+=x**p0*y**p1*R**p2*getINT(cf,10000)
AT=sympy.poly(NORMCS-1).terms()
getNS=0
for tm in AT:
p0,p1,p2=tm[0]
cf=tm[1]
#print(p0,p1,p2,cf)
getNS+=x**p0*y**p1*R**p2*getINT(cf,10000)
print((NORMCS-1).subs([(x,-0.80191751),(y,-0.33680049),(R,1.463)]))
print(ENTS3.subs([(x,-0.80191751),(y,-0.33680049),(R,1.463)]))
print(getNS.subs([(x,-0.80191751),(y,-0.33680049),(R,1.463)]))
e=symbols("e")
#
# GET THE OBJECTIVE FUNCTION
#
OF=getF-2*e*getNS
#
# PREPARES A Singular SCRIPT.
#
# If you are a user of Singular, you simply do as follows.
#
# >> Singular < SCRIPT.txt
#
stringQ='LIB "solve.lib";option(redSB);\n'
stringQ+='ring r=0,(x,y,e,R),dp;\n'+'poly OBJ='+str(OF)+';\n'
stringQ+='list diffs;\n'
stringQ+='for(int i=1;i<=nvars(r); i=i+1){diffs=insert(diffs,diff(OBJ,var(i)));}\n'
stringQ+='poly fR=100*R-146;\n'
stringQ+='ideal I=fR;\n'
stringQ+='for(int i=1;i<=nvars(r)-1; i=i+1){I=I+diff(OBJ,var(i));}\n'
stringQ+='ideal SI=std(I);\n'
stringQ+='ring s=0,(x,y,e,R),lp;\n'
stringQ+='setring s;\n'
stringQ+='ideal j=fglm(r,SI);\n'
stringQ+='def R=triang_solve(j,10);\n'
stringQ+='setring R;rlist;'
text_file = open("SCRIPT.txt", "w")
#write string to file
on = text_file.write(stringQ)
#close file
text_file.close()
stringQ
#
# GET THE DIFFERENTIALS OF THE OBJECTIVE FUNCTION.
#
DS=[]
for t in [x,y,e,R]:
DS.append(sympy.diff(OF,t))
index=0
for i in DS:
print("poly f"+str(index)+"=",i,";")
index+=1
今回は、簡単なコメントを入れたプログラムを示す。このプログラムは、実行終了直前にサブプログラムを生成し、これは「Singular」のスクリプトになっている。
このプログラムではHeH+のRHF計算を行っている。
数値計算の結果は、概ねSabo-Ostlundの結果を再現するものである。ただし、Sabo-Ostlundでは得ることのない、励起状態を基底状態ともに、一つの計算で取得することができる。
以下に示すのは、Singularのスクリプトである。
1．目的関数をpoly OBJに定義
2. (x,y)はHe,Hの1s軌道に対するLCAOの係数
3. eはエネルギー固有値
4. Rは核間距離
5. d(OBJ)/dx, d(OBJ)/dy, d(OBJ)/de, d(OBJ)/dRを計算。
6. イデアルIを生成元（d(OBJ)/dx, d(OBJ)/dy, d(OBJ)/de, 100*R-146)で定義。これは
d(OBJ)/dR=0の代わりに、Rを1.46に固定するという処理。
7. Iのグレブナー基底を作る。
8. グレブナー基底を用いて多項式連立方程式を数値的に解く。この際、準素イデアル分解を利用している。
LIB "solve.lib";option(redSB);
ring r=0,(x,y,e,R),dp;
poly OBJ=-782*R**4*x**3*y - 1034*R**4*x**2*y**2 - 1036*R**4*x**2 - 1037*R**4*x*y**3 + 1799*R**4*x*y + 187*R**4*y**2 + 2633*R**4 + 4977*R**3*x**3*y + 3982*R**3*x**2*y**2 + 8767*R**3*x**2 + 6383*R**3*x*y**3 - 10588*R**3*x*y - 303*R**3*y**2 - 19754*R**3 - 6472*R**2*x**3*y + 5751*R**2*x**2*y**2 - 30430*R**2*x**2 - 9875*R**2*x*y**3 + 6621*R**2*x*y - 5877*R**2*y**2 + 59259*R**2 - 19598*R*x**3*y - 47149*R*x**2*y**2 + 54517*R*x**2 - 12230*R*x*y**3 + 75627*R*x*y + 29906*R*y**2 - 88889*R - 2*e*(-419*R**4*x*y + 3017*R**3*x*y - 6319*R**2*x*y - 2791*R*x*y + 10000*x**2 + 19093*x*y + 10000*y**2 - 10000) + 13071*x**4 + 48009*x**3*y + 68137*x**2*y**2 - 90397*x**2 + 36259*x*y**3 - 153801*x*y + 7746*y**4 - 65726*y**2 + 66666;
list diffs;
for(int i=1;i<=nvars(r); i=i+1){diffs=insert(diffs,diff(OBJ,var(i)));}
poly fR=100*R-146;
ideal I=fR;
for(int i=1;i<=nvars(r)-1; i=i+1){I=I+diff(OBJ,var(i));}
ideal SI=std(I);
ring s=0,(x,y,e,R),lp;
setring s;
ideal j=fglm(r,SI);
def R=triang_solve(j,10);
setring R;rlist;
いくつかの解を得る。複素数の解は除外し、実数解は以下の４つ。
[1][2]が基底状態、[3][4]が励起状態である。
[1]:
[1]:
x:-0.8014629666
[2]:
y:-0.3370604896
[3]:
e:-1.599681024
[4]:
R:1.46
[2]:
[1]:
x:0.8014629666
[2]:
y:0.3370604896
[3]:
e:-1.599681024
[4]:
R:1.46
[3]:
[1]:
x:0.6039377034
[2]:
y:-1.11522499
[3]:
e:-0.5378807636
[4]:
R:1.46
[4]:
[1]:
x:-0.6039377034
[2]:
y:1.11522499
[3]:
e:-0.5378807636
[4]:
R:1.46
結論
この記事では 第一原理量子化学計算（分子軌道計算）を計算代数幾何学の手法に基づいて実行する一例を示した。従来の逐次近似（自己無同着）計算法を用いる必要は全くない。基本的な道具は、解析式の多項式近似と、グレブナー基底の生成、これを用いた数値計算である。