【わずか200行】pythonで量子化学のプログラミングをしよう

2024年10月10日2025年3月25日

量子化学計算は高度なシミュレーションソフト(GaussianやGAMESSなど)が必要と思っていませんか？

実は簡単な分子のハートリー・フォック計算（RHF法）であれば、わずか200行ほどプログラミングするだけで量子化学のシミュレーションができます。

最低限のpython知識だけで作成したRHF法のプログラミングコードを利用して量子化学を体験してみましょう。

HeH+をpythonで量子化学プログラミング
量子化学計算を行うためのインプット
必要となる電子積分の計算
初期分子軌道の準備
SCF計算開始
まとめ

HeH+をpythonで量子化学プログラミング

最低限のpython知識しかなくても、簡単な分子であれば200行程度プログラミングするだけで量子化学計算ができます。（必要となるpythonの知識は配列、ループと条件分岐くらいかと思います。）

以下のプログラミングコードを実行すると、最も簡単な異核二原子分子であるHeH $^{+}$ の量子化学計算（RHF）が実施され、HOMOとLUMOの分子軌道係数(coff_matrix)を得ることができます。

プログラミングコードはこちら

「python_RHF.py」などの名前で全コードを保存し、実行してみてください

「No module named ‘numpy’」とエラーが出たときは、numpyをインストールしてください。

import math
import numpy as np
import numpy.linalg as LA
np.set_printoptions(precision=6,floatmode='fixed')

# 原子数、基底関数数、電子数を入力
n_atm = 2           # HeH+の原子数は2つ
n_sto = 3           # STO-3Gは3つのガウス軌道の和で1つの原子軌道を表す
n_basis = 2         # Heの1s軌道、Hの1s軌道で基底関数数は2つとする
n_elec = 2          # HeH+の電子数は2
n_dim = 3           # x,y,zの3次元
dis_HeH = 1.4       # He-H+の結合長(bohr)

# 配列宣言 
s_matrix = np.zeros((n_basis,n_basis))
t_matrix = np.zeros((n_basis,n_basis))
v_matrix = np.zeros((n_basis,n_basis)) 
e_matrix = np.zeros((n_basis,n_basis,n_basis,n_basis))
m_matrix = np.zeros((n_basis,n_basis)) 
pc=np.zeros((n_dim))
pp=np.zeros((n_dim))
pq=np.zeros((n_dim))
PI = 3.1415926

# 原子核の電荷を入力
CHG = [2,1]

# 原子座標を入力
C = [[0.0000 ,0.0000,0.0000],
     [dis_HeH,0.0000,0.0000]]

# STO-3Gの係数を入力
EXX = [[ 6.362421394 ,1.1589229990,0.3136497915 ],
       [ 3.425250914 ,0.6239137298,0.1688554040 ]]

CXX = [[ 0.4406303066,0.4261584492,0.1328149745],
       [ 0.2769343551,0.2678388516,0.0834736711]]

# STO-3G基底関数の中心座標を入力
CO = [[0.0000 ,0.0000,0.0000],
      [dis_HeH,0.0000,0.0000]]

# 重なり行列を計算 s_t[n][m] = (n|m) 
for n in range(n_basis):
    for m in range(n_basis):
      d2 = (CO[n][0]-CO[m][0])**2+(CO[n][1]-CO[m][1])**2+(CO[n][2]-CO[m][2])**2
      for i in range(n_sto):
          for j in range(n_sto):
              zeta = EXX[n][i]+EXX[m][j]
              gzai = EXX[n][i]*EXX[m][j]/zeta
              sss  = math.exp(-gzai*d2)*((PI/zeta)**1.5)
              cxij = CXX[n][i]*CXX[m][j]
              s_matrix[n,m] +=  cxij * sss
print("s_matrix")                    
print(s_matrix)                    

# 運動エネルギー積分を計算 t_t[n][m] = (n|-1/2 nabla |m) 
for n in range(n_basis):
    for m in range(n_basis):
      d2 = (CO[n][0]-CO[m][0])**2+(CO[n][1]-CO[m][1])**2+(CO[n][2]-CO[m][2])**2
      for i in range(n_sto):
          for j in range(n_sto):
              zeta = EXX[n][i]+EXX[m][j]
              gzai = EXX[n][i]*EXX[m][j]/zeta
              sss  = math.exp(-gzai*d2)*((PI/zeta)**1.5)
              tss  = gzai*(3-2*gzai*d2)*sss
              cxij = CXX[n][i]*CXX[m][j]
              t_matrix[n,m] +=  cxij * tss
print("t_matrix")                    
print(t_matrix)                    

# 核引力積分を計算 v_t[n][m] = (n|1/r|m) 
for n in range(n_basis):
    for m in range(n_basis):
        for l in range(n_atm):
            for ndim in range(n_dim):
                pc[ndim] =C[l][ndim]
            for i in range(n_sto):
                for j in range(n_sto):
                    zeta1 = EXX[n][i]+EXX[m][j]
                    gzai1 = EXX[n][i]*EXX[m][j]/zeta1
                    cxij  = CXX[n][i]*CXX[m][j]
                    for ndim in range(n_dim):
                        pp[ndim] =(EXX[n][i]*CO[n][ndim]+EXX[m][j]*CO[m][ndim])/zeta1
                    d2_ab = (CO[n][0]-CO[m][0])**2+(CO[n][1]-CO[m][1])**2+(CO[n][2]-CO[m][2])**2
                    d2_pc = (pc[0]-pp[0])**2+(pc[1]-pp[1])**2+(pc[2]-pp[2])**2
                    if(d2_pc<0.0001):
                        f0 = 1
                    else:
                        f0 = 0.5*math.sqrt(PI/(zeta1*d2_pc))*math.erf(math.sqrt(zeta1*d2_pc))
                    vss = 2*PI/zeta1*math.exp(-gzai1*d2_ab)*f0
                    v_matrix[n,m] +=  cxij * vss * (-CHG[l])
print("v_matrix")
print(v_matrix)                    

# 2電子積分を計算 e_t[l][m][n][o] = (lm|no) 
for l in range(n_basis):
    for m in range(n_basis):
        for n in range(n_basis):
            for o in range(n_basis):
                d2_ab = (CO[l][0]-CO[m][0])**2+(CO[l][1]-CO[m][1])**2+(CO[l][2]-CO[m][2])**2
                d2_cd = (CO[n][0]-CO[o][0])**2+(CO[n][1]-CO[o][1])**2+(CO[n][2]-CO[o][2])**2
                for h in range(n_sto):
                    for i in range(n_sto):
                        for j in range(n_sto):
                            for k in range(n_sto):
                                zeta1 = EXX[l][h]+EXX[m][i]
                                gzai1 = EXX[l][h]*EXX[m][i]/zeta1
                                zeta2 = EXX[n][j]+EXX[o][k]
                                gzai2 = EXX[n][j]*EXX[o][k]/zeta2
                                rho   = zeta1*zeta2/(zeta1+zeta2)
                                cxijkl= CXX[l][h]*CXX[m][i]*CXX[n][j]*CXX[o][k]
                             
                                for ndim in range(n_dim):
                                    pp[ndim] =(EXX[l][h]*CO[l][ndim]+EXX[m][i]*CO[m][ndim])/zeta1
                                    pq[ndim] =(EXX[n][j]*CO[n][ndim]+EXX[o][k]*CO[o][ndim])/zeta2 

                                d2_pq = (pp[0]-pq[0])**2+(pp[1]-pq[1])**2+(pp[2]-pq[2])**2
                                z_ab = math.sqrt(2)*(PI**1.25)/zeta1*math.exp(-gzai1*d2_ab)
                                z_cd = math.sqrt(2)*(PI**1.25)/zeta2*math.exp(-gzai2*d2_cd)

                                if(d2_pq<0.0001):
                                    f0 = 1
                                else:
                                    f0 = 0.5*math.sqrt(PI/(rho*d2_pq))*math.erf(math.sqrt(rho*d2_pq))

                                e_ssss = 1.0/math.sqrt(zeta1+zeta2)*z_ab*z_cd*f0
                                e_matrix[l,m,n,o] +=  cxijkl * e_ssss
print()                                
print("START SCF CALCULATION")
# S行列を対角化  U' * S * U = M 
w_matrix, u_matrix = LA.eig(s_matrix)
u_d_matrix = u_matrix.T

# 変換行列XとX'をつくる、  X = U * M^(-0.5)  
for i in range(n_basis):
    m_matrix[i][i] = 1/math.sqrt(w_matrix[i])

x_matrix = u_matrix @ m_matrix

# ハミルトニアン行列hをつくり、左右からX'とXをかけてa行列を作る a = X' * h * X  
h_matrix = t_matrix + v_matrix
a_matrix = x_matrix.T @ h_matrix @ x_matrix

# a行列を対角化する固有ベクトルwをつくる w' * a * w = eps 固有値の小さい順番で並び変える
eps_matrix,w_matrix = LA.eig(a_matrix)
index = np.argsort(eps_matrix)[::1]
w_matrix = w_matrix[:,index]

# 変換行列xと固有ベクトルwから初期分子軌道係数coffが得られる coff = x * w
coff_matrix = x_matrix @ w_matrix

# SCF計算開始、電子エネルギーeneが収束するまで行う
ene_save = 0
for n_scf in range(10):
    p_matrix = np.zeros((n_basis,n_basis))
    for n in range(n_basis):
        for m in range(n_basis):
            for i in range( int(n_elec/2) ):
                p_matrix[n][m] +=  2* coff_matrix[n][i] * coff_matrix[m][i] 
    
    # フォック行列hをつくる
    f_matrix =  t_matrix + v_matrix
    for n in range(n_basis):
        for m in range(n_basis):
            for i in range(n_basis):
                for j in range(n_basis):
                    f_matrix[n][m] +=  p_matrix[i][j]*(e_matrix[n][m][i][j]-0.5*e_matrix[n][i][m][j])

    # 電子エネルギーeneを計算し、 ene = 0.5 * p * ( h + f )　出力する
    ene = 0
    for n in range(n_basis):
        for m in range(n_basis):
            ene +=  0.5*p_matrix[n][m]*(h_matrix[m][n]+f_matrix[m][n])
    print("scf_roop",n_scf+1,"elec_ene",round(ene,7))

    # 収束判定を行う
    if(abs(ene - ene_save))<10 ** (-7):
        print("SCF CONVERSION")
        print()
        break

    # 次回roopの時に収束判定を行えるようにするため、エネルギ値をene_saveとして保存する
    ene_save = ene        

    # フォック行列hに、左右からX'とXをかけてa行列を作る a = X' * h * X  
    a_matrix = x_matrix.T @ f_matrix @ x_matrix
    
    # a行列を対角化する固有ベクトルwをつくる w' * a * w = eps 固有値の小さい順番で並び変える
    eps_matrix,w_matrix = LA.eig(a_matrix)
    index = np.argsort(eps_matrix)[::1]
    w_matrix = w_matrix[:,index]

    # 変換行列xと固有ベクトルwから分子軌道係数coffが得られる coff = x * w
    coff_matrix = x_matrix @ w_matrix

# 得られた分子軌道を出力する
print("coff_matrix")
print(coff_matrix)

例えば、coff_matrixの1列目の値 ${- 0.872167, - 0.202795}$ はHOMO軌道におけるHe 1s, H 1sの軌道係数を表しています。

次のセクションから、このプログラミングコードを11ステップに分けて順番に解説していきます

量子化学計算を行うためのインプット

量子化学計算に限らず、何かしらシミュレーションをする上では、前提となるインプット条件が必要になります。STEP1～STEP3まででpythonでHeH $^{+}$ を量子化学計算するためのインプットを記述します。

STEP

原子数、基底関数数、電子数を入力

今回はHeH $^{+}$ の量子化学計算を行うので、HeとHの電荷（ $2, 1$ ）と結合長1.4bohr(0.7408Å)を設定します。

# 原子数、基底関数数、電子数を入力
n_atm = 2           # HeH+の原子数は2つ
n_sto = 3           # STO-3Gは3つのガウス軌道の和で1つの原子軌道を表す
n_basis = 2         # Heの1s軌道、Hの1s軌道で基底関数数は2つとする
n_elec = 2          # HeH+の電子数は2
n_dim = 3           # x,y,zの3次元
dis_HeH = 1.4       # He-H+の結合長(bohr)

STEP

配列宣言

量子化学計算では行列計算が必ず必要なので、行列を表すための配列を宣言します。RHF計算で必要な行列の次元の多くは基底関数の数と同じになります。

# 配列宣言 
s_matrix = np.zeros((n_basis,n_basis))
t_matrix = np.zeros((n_basis,n_basis)) 
（以下略）

STEP

原子核の位置とSTO-3Gの係数を入力

分子軌道に基づく量子化学計算では基底関数の設定が欠かせません。今回はSTO-3Gと呼ばれる基底関数を設定しています。これを使うと、例えばHeの1s軌道は $χ_{He1s} = 0.44 e^{- 6.36 r^{2}}$ $+ 0.42 e^{- 1.15 r^{2}} + 0.13 e^{- 0.31 r^{2}}$ と表現できます。（ $r$ は原子からの距離）

# 原子核の電荷を入力
CHG = [2,1]

# 原子座標を入力
C = [[0.0000 ,0.0000,0.0000],
     [dis_HeH,0.0000,0.0000]]

# STO-3Gの係数を入力
EXX = [[ 6.362421394 ,1.1589229990,0.3136497915 ],
       [ 3.425250914 ,0.6239137298,0.1688554040 ]]

CXX = [[ 0.4406303066,0.4261584492,0.1328149745],
       [ 0.2769343551,0.2678388516,0.0834736711]]

# STO-3G基底関数の中心座標を入力
CO = [[0.0000 ,0.0000,0.0000],
      [dis_HeH,0.0000,0.0000]]

これまでですでに40行超えてるけど、大丈夫？

大丈夫です。最終的に200行（コメント、空行含む）で作成できます。

必要となる電子積分の計算

SCF計算を行う前に、重なり行列 $S$ 、運動エネルギー積分 $T$ 、核引力積分 $V$ 、2電子積分 $E$ を計算する必要があります。STEP4～STEP7でこれらの行列要素を計算していきます。

このページで紹介しているプログラミングコードの内、およそ半分程度（90行）がこの行列計算の部分です。

STEP

重なり行列の計算

2つのs軌道 $χ_{A} = \exp (- a_{i} (r - r_{A})^{2})$ と $χ_{B} = \exp (- b_{j} (r - r_{B})^{2})$ どうしの重なり積分は $ζ_{i j} = a_{i} + b_{j}$ 、 $ξ_{i j} = a_{i} * b_{j} / ζ_{i j}$ と2つの軌道の距離 $d 2_{i j}$ を利用すると $s s s = \exp (- ξ_{i j} * d 2_{i j}) * (\frac{π}{ζ_{i j}})^{1.5}$ で計算できます。

これを利用して、pythonで重なり行列のコードを作成すると次の通りになります。

# 重なり行列を計算 s_matrix[n][m] = (n|m) 
for n in range(n_basis):
    for m in range(n_basis):
      d2 = (CO[n][0]-CO[m][0])**2+(CO[n][1]-CO[m][1])**2+(CO[n][2]-CO[m][2])**2
      for i in range(n_sto):
          for j in range(n_sto):
              zeta = EXX[n][i]+EXX[m][j]
              gzai = EXX[n][i]*EXX[m][j]/zeta
              sss  = math.exp(-gzai*d2)*((PI/zeta)**1.5)
              cxij = CXX[n][i]*CXX[m][j]
              s_matrix[n,m] +=  cxij * sss
print("s_matrix")                    
print(s_matrix)

STEP

運動エネルギー積分計算

2つのs軌道どうしの運動エネルギー積分は $t s s = ξ_{i j} (3 - 2 ξ_{i j} * d 2_{i j}) * s s s$ で計算できます。

これを利用して、pythonで運動量エネルギー積分のコードを作成すると次の通りになります。

# 運動エネルギー積分を計算 t_matrix[n][m] = (n|-1/2 nabla |m) 
for n in range(n_basis):
    for m in range(n_basis):
      d2 = (CO[n][0]-CO[m][0])**2+(CO[n][1]-CO[m][1])**2+(CO[n][2]-CO[m][2])**2
      for i in range(n_sto):
          for j in range(n_sto):
              zeta = EXX[n][i]+EXX[m][j]
              gzai = EXX[n][i]*EXX[m][j]/zeta
              sss  = math.exp(-gzai*d2)*((PI/zeta)**1.5)
              tss  = gzai*(3-2*gzai*d2)*sss
              cxij = CXX[n][i]*CXX[m][j]
              t_matrix[n,m] +=  cxij * tss
print("t_matrix")                    
print(t_matrix)

STEP

核引力積分計算

位置 $(C_{x}, C_{y}, C_{z})$ にある原子核Cによる2つのs軌道の核引力積分は、2つのs軌道の内分点 $(P_{x}, P_{y}, P_{z})$ との距離 $d 2_{p c}$ を利用することで、 $v s s = \frac{2 π}{ζ} \exp (- ξ * d 2_{i j}) * f 0$ で計算できます。ただし、 $f 0 = \int_{0}^{1} \exp (- ζ_{i j} * d 2_{p c} t^{2}) d t$ $= (0.5 * \sqrt{\frac{π}{ζ_{i j} * d 2_{p c}}} erf (ζ_{i j} * d 2_{p c})$ です。

これを利用すると核引力積分のコードを作成できます。

コードはこちら

# 核引力積分を計算 v_matrix[n][m] = (n|1/r|m) 
for n in range(n_basis):
    for m in range(n_basis):
        for l in range(n_atm):
            for ndim in range(n_dim):
                pc[ndim] =C[l][ndim]
            for i in range(n_sto):
                for j in range(n_sto):
                    zeta1 = EXX[n][i]+EXX[m][j]
                    gzai1 = EXX[n][i]*EXX[m][j]/zeta1
                    cxij  = CXX[n][i]*CXX[m][j]
                    for ndim in range(n_dim):
                        pp[ndim] =(EXX[n][i]*CO[n][ndim]+EXX[m][j]*CO[m][ndim])/zeta1
                    d2_ab = (CO[n][0]-CO[m][0])**2+(CO[n][1]-CO[m][1])**2+(CO[n][2]-CO[m][2])**2
                    d2_pc = (pc[0]-pp[0])**2+(pc[1]-pp[1])**2+(pc[2]-pp[2])**2
                    if(d2_pc<0.0001):
                        f0 = 1
                    else:
                        f0 = 0.5*math.sqrt(PI/(zeta1*d2_pc))*math.erf(math.sqrt(zeta1*d2_pc))
                    vss = 2*PI/zeta1*math.exp(-gzai1*d2_ab)*f0
                    v_matrix[n,m] +=  cxij * vss * (-CHG[l])
print("v_matrix")
print(v_matrix)

STEP

2電子積分を計算

4つのs軌道 $χ_{A} = \exp (- a_{i} (r - r_{A})^{2})$ 、 $χ_{B} = \exp (- b_{j} (r - r_{B})^{2})$ 、 $χ_{C} = \exp (- c_{k} (r - r_{C})^{2})$ と $χ_{D} = \exp (- d_{l} (r - r_{D})^{2})$ による2電子積分は $ζ_{1} = a_{i} + b_{j}$ 、 $ζ_{2} = c_{k} + d_{l}$ 、 $ξ_{1} = \frac{a_{i} * b_{j}}{ζ_{1}}$ 、 $ξ_{2} = \frac{c_{k} * d_{l}}{ζ_{2}}$ 、軌道間の距離 $d 2_{A B}$ 、 $d 2_{C D}$ と、軌道を内分点PとQを利用すると、 $e s s s s = \frac{2 π^{5 / 2}}{ζ_{1} * ζ_{2}}$ $\exp (- ξ_{1} * d 2_{a b}) \exp (- ξ_{2} * d 2_{c d}) * f 0$ で表せます。ただし、 $ρ = \frac{ζ_{1} * ζ_{2}}{ζ_{1} + ζ_{2}}$ 、 $f 0 = \int_{0}^{1} \exp (- ρ * d 2_{p q} t^{2}) d t$ です。

初期分子軌道の準備

今の段階で $S$ 、 $V$ 、 $T$ など必要な行列計算が完了しており、 $F C = S C ϵ$ を解けば分子軌道 $C$ を得ることができますが、 $S$ があることによって行列の固有値方程式になっていません。

そのため、STEP8で $F C = S C ϵ$ を固有値方程式 $A C^{'} = ϵ C^{'}$ の形に変換するための変換行列 $X$ を重なり行列 $S$ から作成し、これを利用してSTEP9で初期分子軌道を生成させます。

STEP

S行列から変換行列Xをつくる

変換行列 $X$ は重なり行列 $S$ を $M$ へ対角化する行列 $U$ を使って、 $X = U * M^{- \frac{1}{2}}$ で作成します。

# S行列を対角化  U' * S * U = M 
w_matrix, u_matrix = LA.eig(s_matrix)

# 変換行列XとX'をつくる、  X = U * M^(-0.5)  
for i in range(n_basis):
    m_matrix[i][i] = 1/math.sqrt(w_matrix[i])

x_matrix = u_matrix @ m_matrix

STEP

初期分子軌道を得る

ここまででほぼ下準備は完了です。

変換行列 $X$ を使って、ハミルトニアン行列 $H = T + V$ をA行列（ $A = X^{†} * H * X$ ）に変換し、A行列の固有ベクトル $W$ を使うことで分子軌道を $C = X * W$ で得ることができます。

# ハミルトニアン行列hをつくり、左右からX'とXをかけてa行列を作る a = X' * h * X  
h_matrix = t_matrix + v_matrix
a_matrix = x_matrix.T @ h_matrix @ x_matrix

# a行列を対角化する固有ベクトルwをつくる w' * a * w = eps 固有値の小さい順番で並び変える
eps_matrix,w_matrix = LA.eig(a_matrix)
index = np.argsort(eps_matrix)[::1]
w_matrix = w_matrix[:,index]

# 変換行列xと固有ベクトルwから分子軌道係数coffが得られる coff = x * w
coff_matrix = x_matrix @ w_matrix

SCF計算開始

下準備に長らく時間がかかりましたが、初期分子軌道を得たので、STEP10～STEP11でSCF（繰り返し計算）を行います。

SCF部分は残りの50行ほどです。

STEP

SCF計算開始

繰り返し計算ではすでに得られている分子軌道 $C$ をもとに密度行列 $P$ を得て、電子間反発が加味されたフォック行列 $F$ を計算します。

この密度行列 $P$ とフォック行列 $F$ から電子状態のエネルギー $E$ を計算します。

# SCF計算開始、電子エネルギーeneが収束するまで行う
ene_save = 0
for n_scf in range(10):
    p_matrix = np.zeros((n_basis,n_basis))
    for n in range(n_basis):
        for m in range(n_basis):
            for i in range( int(n_elec/2) ):
                p_matrix[n][m] +=  2* coff_matrix[n][i] * coff_matrix[m][i] 
    
    # フォック行列hをつくる
    f_matrix =  t_matrix + v_matrix
    for n in range(n_basis):
        for m in range(n_basis):
            for i in range(n_basis):
                for j in range(n_basis):
                    f_matrix[n][m] +=  p_matrix[i][j]*(e_matrix[n][m][i][j]-0.5*e_matrix[n][i][m][j])

    # 電子エネルギーeneを計算し、 ene = 0.5 * p * ( h + f )　出力する
    ene = 0
    for n in range(n_basis):
        for m in range(n_basis):
            ene +=  0.5*p_matrix[n][m]*(h_matrix[m][n]+f_matrix[m][n])
    print("scf_roop",n_scf+1,"elec_ene",round(ene,7))

STEP

収束判定

SCF計算ではいつ繰り返しを終了するかを判定しなければなりません。

今回は電子状態のエネルギー $E$ が繰り返し計算前後で値が一定値以下になれば計算を終了することにして、得られた分子軌道を出力するようにしています。

電子状態のエネルギーが一定になっていなければフォック行列 $F$ について、 $F C = S C ϵ$ を解いて、新しい分子軌道 $C$ を得て、計算ループを継続します。

# 収束判定を行う
    if(abs(ene - ene_save))<10 ** (-7):
        print("SCF CONVERSION")
        print()
        break

    # 次回roopの時に収束判定を行えるようにするため、エネルギ値をene_saveとして保存する
    ene_save = ene        

    # フォック行列hに、左右からX'とXをかけてa行列を作る a = X' * h * X  
    a_matrix = x_matrix.T @ f_matrix @ x_matrix
    
    # a行列を対角化する固有ベクトルwをつくる w' * a * w = eps 固有値の小さい順番で並び変える
    eps_matrix,w_matrix = LA.eig(a_matrix)
    index = np.argsort(eps_matrix)[::1]
    w_matrix = w_matrix[:,index]

    # 変換行列xと固有ベクトルwから分子軌道係数coffが得られる coff = x * w
    coff_matrix = x_matrix @ w_matrix

# 得られた分子軌道を出力する
print("coff_matrix")
print(coff_matrix)