#はじめに
制限Hartree-Fockの実装（Part1) の続きです。
前回は抽象的な話が多かったので、今回は実際にネオン原子で考えていきます。
Python3.6でのコードは一番下にまとめてあります。
#ネオン原子について制限Hartree-Fockを用いる
一度制限Hartree-Fockの流れをまとめておきます。
(1) 核電荷、核座標を決める。
(2) 基底関数系を決める。
(3) 重なり行列Sを求める。
(4) 係数行列Cの初期値を設定（実際には密度行列Pの初期値を決める）。
(5) 重なり行列Sから変換行列Xを求める。
(6) Roothaanの行列方程式を解く($FC=SC\epsilon$)
(6-1) Fock行列Fを決定する。
(6-2) Fock行列を変換する($F'=X^{†}FX$)
(6-3) F'を対角化してC'と $\epsilon$を得る
(6-4) Cを得る($C=XC'$)(実際はここから密度行列Pを計算する)
(7) 設定した係数行列C(密度行列P)が新しく得た係数行列C(密度行列P)と十分に一致しているか確認する。十分なかったら新しく得た係数行列C(密度行列P)を元に(6)の行程を繰り返す。

ここから実際にネオン原子について解いていきます。
#(1)核電荷、核座標を決める。
ネオン原子なので各電荷は$Z_{A}=10$、核座標は$R_{A}=(0,0,0)$
#(2) 基底関数系を決める
今回はSTO-3Gという基底関数系を決めます。STO-3Gは基底関数がガウス型関数の線形結合であらわされ、軌道指数 $\xi_i$と係数 $d_{ij}$ が原子ごとに決まっています。
[1]新しい量子化学と[2]Self-Consistent Molecular Orbital Methods.IV.より軌道指数 $\xi_i$と係数 $d_{ij}$　を使ったネオン原子の基底関数は以下のようになります。

\phi_{1s}=0.444635\phi_{s}(10.2053)+0.535328\phi_{s}(37.7081)+0.154329\phi_{s}(207.0155)\\
\phi_{2s}=0.700115\phi_{s}(0.6232)+0.399513\phi_{s}(1.9162)-0.0999672\phi_{s}(8.2463)\\
\phi_{2px}=0.391957\phi_{px}(0.6232)+0.607684\phi_{px}(1.9162)+0.155916\phi_{px}(8.2463)\\
\phi_{2py}=0.391957\phi_{py}(0.6232)+0.607684\phi_{py}(1.9162)+0.155916\phi_{py}(8.2463)\\
\phi_{2pz}=0.391957\phi_{pz}(0.6232)+0.607684\phi_{pz}(1.9162)+0.155916\phi_{pz}(8.2463)\\
\Biggl(\phi_s(\alpha)=\biggl(\frac{8\alpha^3}{\pi^3}\biggr)^{\frac{1}{4}}e^{-\alpha r^2}\Biggr)\\
\Biggl(\phi_{px}(\alpha)=\biggl(\frac{128\alpha^5}{\pi^3}\biggr)^{\frac{1}{4}}xe^{-\alpha r^2}\Biggr)

#(3) 重なり行列Sを求める。
重なり行列Sは以下のようになります。

S=
\left(
\begin{array}{ccccc}
1.00000143& 0.2427812 & 0 & 0 & 0  \\
 0.2427812 & 0.99999979 & 0 & 0 & 0 \\
0 &  0 & 0.99999979 & 0 & 0\\
0 &  0 & 0 & 0.99999979 & 0\\
0 &  0 & 0 & 0 & 0.99999979\\
\end{array}\right)\\
S_{ij}=\int dr_1\phi^*_{i}(r_1)\phi_{j}(r_1)\\
以下\phi_{1s}=\phi_{1} \quad \phi_{2s}=\phi_{2} \quad \phi_{2px}=\phi_{3} \quad \phi_{2py}=\phi_{4} \quad \phi_{2pz}=\phi_{5}とする\\

#(4) 係数行列Cの初期値を設定（実際には密度行列Pの初期値を決める）。
密度行列Pを以下のように定めます。

P=
\left(
\begin{array}{ccccc}
0 & 0 & 0 & 0 & 0  \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0\\
\end{array}\right)\\

また密度行列Pと係数行列Cには以下の関係があります。

P_{ij}=2\sum_{a=1}^{5}C_{ia}C^*_{ja}

#(5) 重なり行列Sから変換行列Xを求める。
変換行列Xを正準直交化で求めます。
重なり行列Sの固有値を並べた行ベクトルをsとします。
重なり行列Sの固有ベクトルを対応する固有値を並べた順に左から並べ5×5の正方行列をUとします。（つまり対角行列）
そうすると

X=Us^{-\frac{1}{2}}\\
X_{ij}=U_{ij}/s^{\frac{1}{2}}_j

このXがSを単位行列に変換する行列です。

X^{†}SX=1 \quad (1は単位行列)

具体的には

X=
\left(
\begin{array}{ccccc}
0.63429062 & -0.81259361 & 0 & 0 & 0  \\
0.63428848 & 0.81259635 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0\\
0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 1\\
\end{array}\right)\\

#(6) Roothaanの行列方程式を解く
##(6-1) Fock行列Fを決定する
Fock行列は以下の二つの部分に分けることができ前半の部分は繰り返しの過程のなかで変化しません。

\begin{eqnarray*}
 F_{ij}&=&\int dr_1\phi*_{i}(r_1)h(1)\phi_{j}(r_1)+\sum_{a=1}^{5}\int dr_1\phi*_{i}(r_1)[2J_a(1)-K_a(1)]\phi_{j}(r_1)\\
 &=&H^{core}_{ij}+\sum_{a=1}^{5}2(ij|aa)-(ia|aj)
 \end{eqnarray*}\\
 (ij|kl)=\iint dr_1dr_2\phi^*_i(r_1)\phi_j(r_1)r^{-1}_{12}\phi^*_k(r_2)\phi_l(r_2)と定義される。

クーロン演算子、交換演算子は上で定義されている物を空間軌道に変えただけのものです。
ここで交換演算子の半分の項が消えているのはスピン軌道を空間軌道に直したときの計算の過程でスピン関数の直交性によって消えたからです。（これが交換演算子の表現行列は常に正なのでこの項はエネルギーを下げる役割を持ちスピンスピン相互作用につながる）
さらにインデックスaは求めたい軌道について取っているので、展開行列Cと基底関数で展開することができ以下のようになります。

\begin{eqnarray*}
F_{ij}&=&H^{core}_{ij}+\sum_{a=1}^{5}\sum_{k=1,l=1}^{5,5}C_{ka}C^*_{la}[2(ij|lk)-(ik|kl)]\\
&=&H^{core}_{ij}+\sum_{k=1,l=1}^{5,5}P_{kl}[(ij|lk)-\frac{1}{2}(ik|kl)]\\
&=&H^{core}_{ij}+G_{ij}
\end{eqnarray*}

実際に計算すると(Pの初期条件より $H_{core}=F$となるのは明らかです。)

H_{core}=F=
\left(
\begin{array}{ccccc}
-49.4245 & -11.8586 & 0 & 0 & 0  \\
-11.8586 & -13.1590 & 0 & 0 & 0 \\
0 & 0 & -3.0404 & 0 & 0\\
0 & 0 & 0 & -3.0404 & 0\\
0 & 0 & 0 & 0 & -3.0404\\
\end{array}\right)\\

##(6-2)Fock行列を変換する

F'=X^{†}FX=\left(
\begin{array}{ccccc}
-34.72091 & 18.69193 & 0 & 0 & 0  \\
18.69193 & -25.66367 & 0 & 0 & 0 \\
0 & 0 & -3.04049 & 0 & 0\\
0 & 0 & 0 & -3.04049 & 0\\
0 & 0 & 0 & 0 & -3.04049\\
\end{array}\right)\\

##(6-3) F'を対角化してC'を得る
F'を対角化してC'と $\epsilon$を得ます。

C'=\left(
\begin{array}{ccccc}
-0.785959 & 0.618278 & 0 & 0 & 0  \\
0.618278 & -0.785959 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0\\
0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 1\\
\end{array}\right)\\
\epsilon=\left(
\begin{array}{c}
-49.425000 & -10.959588 & -3.040491 & -3.040491 & -3.040491
\end{array}\right)\\

##(6-4) Cを得る(実際はここから密度行列Pを計算する)
C'からCを得ます。

C=XC'=\left(
\begin{array}{ccccc}
-1.00093 & 0.24649767 & 0 & 0 & 0  \\
0.003885 & -1.030834 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0\\
0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 1\\
\end{array}\right)\\

#(7)収束したか確認する
ここで収束したか確認します。今回はエネルギーEについて考え、新しいE’と前のEの絶対値の差が$10^{-6}$ 以下だったら収束したことにします。
エネルギーEと密度行列Pは以下のように関係しています。

E=\frac{1}{2}\sum_{i=1}^{5}\sum_{j=1}^{5}P_{ij}(H^{core}_{ij}+F_{ij})

(6-4)で得られたCからPを、PからGを、GからFを計算して上のエネルギーEの式に入れると以下のようになります。

P=\left(
\begin{array}{ccccc}
2.003743 & -7.778163×10^{-3} & 0 & 0 & 0  \\
-7.778163×10^{-3} & 3.019339×10^{-5} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0\\
\end{array}\right)\\
G=\left(
\begin{array}{ccccc}
6.025817 & 0.963170 & 0 & 0 & 0  \\
0.963170 & 2.68422 & 0 & 0 & 0 \\
0 & 0 & -7.90917 & 0 & 0\\
0 & 0 & 0 & -7.90917 & 0\\
0 & 0 & 0 & 0 & -7.90917\\
\end{array}\right)\\
F=\left(
\begin{array}{ccccc}
-43.39870 & -10.89551 & 0 & 0 & 0  \\
-10.89551 & -10.47480 & 0 & 0 & 0 \\
0 & 0 & -10.94966 & 0 & 0\\
0 & 0 & 0 & -10.94966 & 0\\
0 & 0 & 0 & 0 & -10.94966\\
\end{array}\right)\\
E=-92.82035578404\\

なおエネルギーの単位は原子単位系であることに注意してください。
これを1回目として得られたPから再び(6)の行列方程式を解いて2回目の操作を行います。もし１回目のエネルギーと2回目のエネルギーの絶対値の差が $10^{-6}$ 以下だったら終了となり、$10^{-6}$ より大きかったら3回目の計算を同様に行います。
以下で各段階のエネルギーEと収束したときのC,$\epsilon$,P,E,Fを挙げます。

\begin{eqnarray*}
&1 回目:&-92.8203557840434\\
&2 回目:&-92.83520766588477\\
&3 回目:&-92.83522387033041\\
&4 回目:&-92.83522388775538(収束)\\
\end{eqnarray*}\\
C=\left(
\begin{array}{ccccc}
-0.997182 & 0.26126531 & 0 & 0 & 0  \\
-0.011351 & -1.030779 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0\\
0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 1\\
\end{array}\right)\\
\epsilon=\left(
\begin{array}{c}
-43.418609 & -8.223764 & -10.944756 & -10.944756 & -10.944756
\end{array}\right)\\
P=\left(
\begin{array}{ccccc}
1.988746 & 2.263897×10^{-2} & 0 & 0 & 0  \\
2.263897×10^{-2} & 2.577116×10^{-4} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0\\
\end{array}\right)\\
F=\left(
\begin{array}{ccccc}
-43.414422 & -10.91614 & 0 & 0 & 0  \\
-10.91614 & -10.484558 & 0 & 0 & 0 \\
0 & 0 & -10.94474 & 0 & 0\\
0 & 0 & 0 & -10.94474 & 0\\
0 & 0 & 0 & 0 & -10.94474\\
\end{array}\right)\\

#結果・考察

係数行列Cより1つ目の軌道は水素様原子の主に1sを,さらに1%程度の2sを含んでいる。
同様に2つ目の軌道は水素様原子の主に1sを,さらに25%程度の2sを含んでいる。
3つ目、4つ目、5つ目の軌道はそれぞれ水素様原子の2px、2py、2pzである。
3つ目、4つ目、5つ目の軌道は縮退している。
2つ目の2s軌道に相当する軌道が3つ目、4つ目、5つ目の軌道のエネルギーより低くなっており実際とは異なっている。

#励起エネルギー
つぎに基底関数に3s軌道も含ませ励起状態を考えます。この時の軌道指数は $\xi=2.1$ とします。密度行列Pと係数行列Cの関係で電子が入っている部分でのみ総和をとることに気を付けると
$$
E_0=-92.90533001571535(基底状態)\
E_1=-92.83767222130514(励起状態)
$$
このエネルギー差に対する光の波長は
$$
E=h\nu=\frac{hc}{\lambda}を用いて\
\lambda=673.24(nm)\
$$
ネオンランプの波長は580~700nmなので比較的よい結果が得られたと思います。
今後Hartree-Fockを改良する方法として主に二つあります。

基底関数系をより大きな関数系にする。（関数系を無限大まで大きくしたものをHartree-Fock極限と言います。)
基底状態のスレーター行列式と励起状態のスレーター行列の一部を基底としてSchrodinger方程式をとく。(Post-Hartree-Fockの一つ)

#Python3.6のコード

Gauss.py

import numpy as np
import math
import numpy.linalg as LA
from Basis import *
L = 3


class Gauss:
    def __init__(self, a, R, shape, dire, C=1):
        self.a = a
        self.R = R
        self.C = C
        self.shape = shape
        self.dire = dire

    def func(self, r):
        return np.exp(-a * LA.norm(r - r.R) ** 2)

    def __mul__(self, other):
        K = np.exp(-a * b / (a + b) * LA.norm(self.R - other.R) ** 2)
        Rp = (a * self.R + b * other.R) / (a + b)
        p = a + b
        return Gauss(p, Rp, K)


def sI(GA, GB):
    SA, SB = GA.shape, GB.shape
    DA, DB = GA.dire, GB.dire
    a, Ra, Ca, b, Rb, Cb = GA.a, GA.R, GA.C, GB.a, GB.R, GB.C
    if SA == "s" and SB == "s":
        x = (math.pi / (a + b)) ** (3 / 2)
        y = np.exp(-a * b / (a + b) * LA.norm(Ra - Rb) ** 2)
        return Ca * Cb * x*y
    if SA == SB and DA == DB:
        x = (1 / (2 * (a + b))) * (math.pi / (a + b)) ** (3 / 2)
        return Ca * Cb * x
    else:
        return 0


def Suv(i, j, d, CGF):
    num = 0
    for p in range(L):
        for q in range(L):
            num += d[p][i] * d[q][j] * sI(CGF[i][p], CGF[j][q])
    return num


def S(d, CGF):
    n = CGF.shape[0]
    M = np.zeros((n, n))
    for u in range(n):
        for v in range(n):
            M[u][v] = Suv(u, v, d, CGF)
    return M


def tI(GA, GB):
    SA, SB = GA.shape, GB.shape
    DA, DB = GA.dire, GB.dire
    a, Ra, Ca, b, Rb, Cb = GA.a, GA.R, GA.C, GB.a, GB.R, GB.C
    if SA == "s" and SB == "s":
        x = a * b / (a + b)
        y = 3 - 2 * a * b / (a + b) * LA.norm(Ra - Rb) ** 2
        z = (math.pi / (a + b)) ** (3 / 2)
        w = np.exp(-a * b / (a + b) * LA.norm(Ra - Rb) ** 2)
        return Ca * Cb * x * y * z * w
    if SA == SB and DA == DB:
        x=a*b/(a+b)
        y=5/(2*(a+b))
        z = (math.pi / (a + b)) ** (3 / 2)
        return Ca*Cb*x*y*z
    else:
        return 0


def Tuv(i, j, d, CGF):
    num = 0
    for p in range(L):
        for q in range(L):
            num += d[p][i] * d[q][j] * tI(CGF[i][p], CGF[j][q])
    return num


def T(d, CGF):
    n = CGF.shape[0]
    M = np.zeros((n, n))
    for u in range(n):
        for v in range(n):
            M[u][v] = Tuv(u, v, d, CGF)
    return M


def F0(x):
    if x == 0:
        return 1
    return 1 / 2 * (math.pi / x) ** (1 / 2) * math.erf(x ** (1 / 2))


def vI(GA, GB, Z, Rc):
    SA, SB = GA.shape, GB.shape
    DA, DB = GA.dire, GB.dire
    a, Ra, Ca, b, Rb, Cb = GA.a, GA.R, GA.C, GB.a, GB.R, GB.C
    Rp = (a * Ra + b * Rb) / (a + b)
    if SA == "s" and SB == "s":
        x = -2 * math.pi / (a + b) * Z
        y = np.exp(-a * b / (a + b) * LA.norm(Ra - Rb) ** 2)
        z = F0((a + b) * LA.norm(Rp - Rc) ** 2)
        return Ca * Cb * x * y * z
    if SA==SB and DA==DB:
        x=-math.pi*Z
        y=3*(a+b)**2
        return Ca*Cb*x/y
    else:
        return 0


def Vuv(i, j, d, CGF, Z, Rc):
    num = 0
    for p in range(L):
        for q in range(L):
            num += d[p][i] * d[q][j] * vI(CGF[i][p], CGF[j][q], Z, Rc)
    return num


def V_i(d, CGF, Z, Rc):
    n = CGF.shape[0]
    M = np.zeros((n, n))
    for u in range(n):
        for v in range(n):
            M[u][v] = Vuv(u, v, d, CGF, Z, Rc)
    return M


def H_core(d, CGF, Zlist, Rclist):
    M = T(d, CGF)
    for i in range(Zlist.size):
        M += V_i(d, CGF, Zlist[i], Rclist[i])
    return M


def twoI(GA, GB, GC, GD):
    j=change([GA.shape,GA.dire],[GB.shape,GB.dire],[GC.shape,GC.dire],[GD.shape,GD.dire])
    g=[]
    for k in j[0]:
        g.append([GA,GB,GC,GD][k])
    GA,GB,GC,GD=g[0],g[1],g[2],g[3]
    SA, SB = GA.shape, GB.shape
    a, Ra, Ca, b, Rb, Cb = GA.a, GA.R, GA.C, GB.a, GB.R, GB.C
    c, Rc, Cc, d, Rd, Cd = GC.a, GC.R, GC.C, GD.a, GD.R, GD.C
    Rp = (a * Ra + b * Rb) / (a + b)
    Rq = (c * Rc + d * Rd) / (c + d)
    p=a+b
    q=c+d
    C=Ca*Cb*Cc*Cd
    if j[1]==0:
       return 0
    if j[1]==1:
        x = 2 * math.pi ** (5 / 2)
        y = (a + b) * (c + d) * (a + b + c + d) ** (1 / 2)
        z = np.exp(-a * b / (a + b) * LA.norm(Ra - Rb) ** 2 - c * d / (c + d) * LA.norm(Rc - Rd) ** 2)
        w = F0((a + b) * (c + d) / (a + b + c + d) * LA.norm(Rp - Rq) ** 2)
        return Ca * Cb * Cc * Cd * x / y * z * w
    if j[1]==2:
        x=math.pi**(5/2)/2
        y=1/(6*p**2*(p+q)**(3/2))
        z=4/(p**2*q*(p+q)**(1/2))
        return C*x*(y-z)
    if j[1]==3:
        x=math.pi**2
        y=6*p*q*(p+q)**(3/2)
        return C*x/y
    if j[1]==4:
        x=2*(p*q/(p+q))**(1/2)*math.pi**(1/2)*(math.pi/p)**(3/2)*(math.pi/p)**(3/2)
        y=1/(4*p*q)
        z=3/(20*(p+q)**2)
        w=(-1/12)*(1/(p*(p+q))+1/(q*(p+q)))
        return C*x*(y+z+w)
    if j[1]==5:
        x = 2 * (p * q / (p + q)) ** (1 / 2) * math.pi ** (1 / 2) * (math.pi / p) ** (3 / 2) * (math.pi / p) ** (3 / 2)
        y = 1 / (4 * p * q)
        z = 1 / (20 * (p + q) ** 2)
        w = (-1 / 12) * (1 / (p * (p + q)) + 1 / (q * (p + q)))
        return C * x * (y + z + w)
    if j[1]==6:
        x = 2 * (p * q / (p + q)) ** (1 / 2) * math.pi ** (1 / 2) * (math.pi / p) ** (3 / 2) * (math.pi / p) ** (3 / 2)
        y= 1 / (20 * (p + q) ** 2)
        return C*x*y
    if j[1]==7:
        x = 2 * (p * q / (p + q)) ** (1 / 2) * math.pi ** (1 / 2) * (math.pi / p) ** (3 / 2) * (math.pi / p) ** (3 / 2)
        y = 1 / (4 * p * q)
        z = -3 / (64 * (p + q) ** 2)
        w = (-1 / 12) * (1 / (p * (p + q)) + 1 / (q * (p + q)))
        return C * x * (y + z + w)


def TWO(u, v, w, x, CGF,d):
    num = 0
    for p in range(L):
        for q in range(L):
            for r in range(L):
                for s in range(L):
                    num += d[p][u] * d[q][v] * d[r][w] * d[s][x] * twoI(CGF[u][p], CGF[v][q], CGF[w][r], CGF[x][s])
    return num
def change(Ashapes,Bshapes,Cshapes,Dshapes):
    Ashape,Bshape,Cshape,Dshape=Ashapes[0],Bshapes[0],Cshapes[0],Dshapes[0]
    shapes=[Ashape,Bshape,Cshape,Dshape]
    Adire,Bdire,Cdire,Ddire=Ashapes[1],Bshapes[1],Cshapes[1],Dshapes[1]
    dires=[Adire,Bdire,Cdire,Ddire]
    shapelist=["s","p"]
    direlist=[0,"x","y","z"]
    if shapes.count("s")==4:
        return [(0,1,2,3),1]
    if shapes.count("s")==3:
        p=shapes.index("p")
        return [(0,1,2,3),0]
    if shapes.count("s")==2:
        if dires.count("x")==2 or dires.count("y")==2 or dires.count("z")==2:
            if Ashape==Bshape=="p":
                return [(0,1,2,3),2]
            elif Cshape==Dshape=="p":
                return [(2,3,0,1),2]
            elif Ashape==Cshape=="p":
                return [(0,1,2,3),3]
            elif Bshape==Dshape=="p":
                return [(1,0,3,2),3]
            elif Ashape==Dshape=="p":
                return [(0,1,3,2),3]
            elif Bshape==Cshape=="p":
                return [(1,0,2,3),3]
        else:
            return [(0,1,2,3),0]
    if shapes.count("s")==1:
        return [(0,1,2,3),0]
    if shapes.count("s")==0:
        if dires.count("x")*dires.count("y")*dires.count("z")>=2:
            return [(0,1,2,3),0]
        elif dires.count("x")==4 or dires.count("y")==4 or dires.count("z")==4:
            return [(0,1,2,3),4]
        elif dires.count("x") == 3 or dires.count("y") == 3 or dires.count("z") == 3:
            j=[0,1,2,3]
            if dires.count("x")==3:
                if dires.count("y")==1:
                    i=dires.index("y")
                    j.remove(i)
                    j.append(i)
                    return [j,5]
                elif dires.count("z")==1:
                    i=dires.index("z")
                    j.remove(i)
                    j.append(i)
                    return [j,5]
            elif dires.count("y")==3:
                if dires.count("x")==1:
                    i=dires.index("x")
                    j.remove(i)
                    j.append(i)
                    return [j,5]
                elif dires.count("z")==1:
                    i=dires.index("z")
                    j.remove(i)
                    j.append(i)
                    return [j,5]
            elif dires.count("z")==3:
                if dires.count("x")==1:
                    i=dires.index("x")
                    j.remove(i)
                    j.append(i)
                    return [j,5]
                elif dires.count("y")==1:
                    i=dires.index("y")
                    j.remove(i)
                    j.append(i)
                    return [j,5]
        elif dires.count("x")==2 or dires.count("y")==2 or dires.count("z")==2:
            if Adire==Bdire:
                return [(0,1,2,3),7]
            elif Adire==Cdire:
                return [(0,1,2,3),6]
            elif Adire==Ddire:
                return [(0,1,3,2),6]

Basis.py

from Gauss import *
def convert(a, xil):
    for i in xil:
        xi=i[0]
        for j in i[1]:
            a[j] = a[j]*xi**2
    return a

Hartree-FockNe.py

Ra = np.array([0, 0, 0])
d = np.array([[0.444635, 0.700115, 0.391957, 0.391957, 0.391957],
              [0.535328, 0.399513, 0.607684, 0.607684, 0.607684],
              [0.154329, -0.0999672, 0.155916, 0.155916, 0.155916]])#1s,2s,,2p,2p,2p

xi1 = 9.64
xi2 = 2.88
xi3 = 2.1
a_1 = np.array([[0.109818, 0.405771, 2.22766],
              [0.0751386, 0.231031, 0.994203],
              [0.0751386, 0.231031, 0.994203],
              [0.0751386, 0.231031, 0.994203],
              [0.0751386, 0.231031, 0.994203]])  #上から1s,2s,2p,2p,2p,3s
a = convert(a_1, [[xi1, [0]], [xi2, [1, 2, 3, 4]]])
K1 = (2*a/math.pi)**(3/4)
K2 = (128*a**5/math.pi**3)**(1/4)

CGF = np.array([[Gauss(a[0][0], Ra, "s", 0, K1[0][0]), Gauss(a[0][1], Ra, "s", 0, K1[0][1]), Gauss(a[0][2], Ra, "s", 0,K1[0][2])],
              [Gauss(a[1][0], Ra, "s", 0, K1[1][0]), Gauss(a[1][1], Ra, "s", 0, K1[1][1]), Gauss(a[1][2], Ra, "s", 0,K1[1][2])],
              [Gauss(a[2][0], Ra, "p", "x", K2[2][0]), Gauss(a[2][1], Ra, "p", "x", K2[2][1]), Gauss(a[2][2], Ra, "p", "x", K2[2][2])],
              [Gauss(a[3][0], Ra, "p", "y", K2[3][0]), Gauss(a[3][1], Ra, "p", "y", K2[3][1]), Gauss(a[3][2], Ra, "p", "y", K2[3][2])],
              [Gauss(a[4][0], Ra, "p", "z", K2[4][0]), Gauss(a[4][1], Ra, "p", "z", K2[4][1]), Gauss(a[4][2], Ra, "p", "z", K2[4][2])]])
Zlist = np.array([10])
Rclist = np.array([Ra])
def makeTwo(d, CGF):
    n = CGF.shape[0]
    Twodict = dict()
    for u in range(n):
        vdict = dict()
        for v in range(n):
            wdict = dict()
            for w in range(n):
                xdict = dict()
                for x in range(n):
                    xdict[x] = TWO(u, v, w, x, CGF, d)
                wdict[w] = xdict
            vdict[v] = wdict
        Twodict[u] = vdict
    return Twodict


def koyuu(S):
    koyuu, U = LA.eig(S)
    koyuu_ = np.sort(koyuu)
    U_= np.zeros((U.shape[0], U.shape[0]))
    for i in range(koyuu.shape[0]):
        I = koyuu_[i]
        k = list(koyuu).index(I)
        koyuu[k] = 0
        for j in range(U.shape[0]):
            U_[j][i] = U[j][k]
    return (koyuu_, U_)

def X_syn(S):
    s, U = koyuu(S)
    s_12 = np.zeros((U.shape[0], U.shape[0]))
    for i in range(U.shape[0]):
        s_12[i][i] = s[i] ** (-1 / 2)
    X = U @ s_12 @ U.T
    return X


def X_can(S):
    #s, U = koyuu(S)
    s,U=LA.eig(S)
    s_12 = np.zeros((U.shape[0], U.shape[0]))
    for i in range(U.shape[0]):
        s_12[i][i] = s[i]**(-1/2)
    X = U@s_12
    return X


def calc_G(P, CGF):
    n = CGF.shape[0]
    def Guv(u, v, P):
        Guv = 0
        for i in range(n):
            for k in range(n):
                Guv += P[i, k]*(Twodict[u][v][k][i]-Twodict[u][i][k][v]*(1/2))
        return Guv
    M=np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            M[i][j] = Guv(i, j, P)
    return M


def calc_P(C, elecnum):
    n = CGF.shape[0]
    def Puv(u, v, C, elecnum):
        Puv = 0
        for a in range(int(elecnum/2)):
            Puv += C[u, a]*C[v, a]
        Puv *= 2
        return Puv
    M = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            M[i][j] = Puv(i, j, C, elecnum)
    return M

def calc_P_(C, elecnum):
    n = CGF.shape[0]
    def Puv_(u, v, C, elecnum):
        Puv_ = 0
        for a in range(int(elecnum/2)):
            Puv_ += C[u, a]*C[v, a]
        Puv_ *= 2
        Puv_ -= C[u, 2]*C[v, 2]
        Puv_ += C[u, 5]*C[v, 5]
        return Puv_
    M = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            M[i][j] = Puv_(i, j, C, elecnum)
    return M


def E_0(P, F):
    E_0 = 0
    for u in range(CGF.shape[0]-1):
        for v in range(CGF.shape[0]-1):
            E_0 += P[v, u]*(H[u, v]+F[u, v])
    E_0 /= 2
    return E_0


def E_1(P,F):
    E_1 = 0
    for u in range(CGF.shape[0]):
        for v in range(CGF.shape[0]):
            if u == 4:
                pass
            elif v == 4:
                pass
            else:
                E_1 += 1/2*P[v, u]*(H[u, v]+F[u, v])
    E_1 = E_1
    #E_1 /= 2
    return E_1


def HF(F, P, elecnum, n=1, E = 0):
    E_ = E
    F_prime = X_can1.T@F@X_can1
    e, C_prime = LA.eig(F_prime)
    #e, C_prime = koyuu(F_prime)
    C = X_can1@C_prime
    print("C:",E_0(P, F))
    P = calc_P(C, elecnum)
    G = calc_G(P, CGF)
    F = H+G
    E = E_0(P, F)
    print(n, "回目:", E, P)
    if n == 1:
        print("C:", C)
        print("e:", e)
        print("P:", P)
        print("G:", G)
        print("F:", F)
        print("E:", E)
        HF(F, P, elecnum, n+1, E)
    else:
        if abs(E_-E)<=1*10**(-6):
            print("C:", C)
            print("e:", e)
            print("P:", P)
            print("F:", F)
            print("E_0:",E)
            print("E_1:",E_1(P,F))
            print("波長は：", 45.55 / (E_1(P, F) - E), "nm")
            return 0
        else:
            HF(F, P, elecnum, n+1, E)
P = np.zeros((CGF.shape[0], CGF.shape[0]))
Twodict = makeTwo(d, CGF)
S1 = S(d, CGF)

X_can1=X_can(S1)
H = H_core(d, CGF, Zlist, Rclist)
F = H+calc_G(P, CGF)
HF(F, P, 2)

#参考文献
[1] A.ザボ, N.S.オストランド『新しい量子化学―電子構造の理論入門〈上〉』東京大学出版会（1987）
[2] W.J.HEHRE, R.DITCHFIELD, R.F.STEWART, and J.A.POPLE「Self-Consistent Molecular Orbital Methods.IV Use of Gaussian Expansions of Slater-Type Orbitals.Extension to Second-Row Molecles」THE JOURNAL OF CHEMICAL PHYSICS VOLUME 52,NUMBER 5, 1 MARCH 1970