Ab initio вычисление молекулы водорода

Дан файл с фунциями для расчёта основных матриц hfscf.py, откуда импортирую функции:

from hfscf import *

Задаю функцию расчёта орбиталей на основе подхода Рутана Холла:

(*) Объекты: S - матрица перекрывания, F - матрица Фока, X - диагонализированная матрица S, H - матрица гамильтониана ядра, C - матрица коэффициентов.

def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=False):
    R = [0, r]
    if vbs: print("*) Генерирую матрицу перекрывания S.")
    s_scf = S(R,b1,b2,b)
    if vbs: print("\n*) Генерирую матрицу гамильтониана ядра H.")
    h_scf = H(R,Z,b1,b2,b)
        
    # Диагонализирую матрицу S и рассчитываю сопряженную матрицу X
    if vbs: print("\n*) Диагонализирую матрицу S и рассчитываю сопряженную матрицу X")
    X = diagon(m=s_scf)
    Xa = X.getH()
    
    # Создаю матрицу плотности P
    if vbs: print("\n*) Создаю матрицу плотности P")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)  
    
    # С помощью итераций "улучшаю" матрицу P
    if vbs: print("\n*) Запуск SCF")
    for iteraction in range(50):
        # Создаю матрицу Фока F
        # F = H + G
        if vbs: print("\n**) Генерирую матрицу Фока")
        g_scf = G(r,p_scf,b1,b2,b)
        f_scf = h_scf + g_scf  
        
        # Строю матрицу F'
        # F' = X_adj * F * X
        if vbs: print("**) Строю матрицу F'")
        f_tra = Xa * f_scf * X
        
        # Диагонализирую матрицу F' и строю матрицу C'
        if vbs: print("**) Диагонализирую матрицу F' и строю матрицу C'")
        c_tra = diagon2(m=f_tra)
        
        # Строю матрицу C
        # C = X * C'
        if vbs: print("**) Строю матрицу коэффициентов C")
        c_scf = X * c_tra
        
        # Пересчитаю матрицу P на основе матрицы C
        if vbs: print("**) Пересчитаю матрицу плотности P")
        p_temp = P(C=c_scf)
        
        print("\n Произведено итераций: " + str(iteraction + 1))
        
        # Проверяю сходимость
        if np.linalg.norm(p_temp - p_scf) < 1E-4: 
            print("--> Сходится")
            return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}
        else:
            p_scf = p_temp
    print("--> Не сходится")
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}

Рассчитаю электронную энергию молекулы и полную энергию молекулы:

orb = SCF(vbs=True)

*) Генерирую матрицу перекрывания S.

*) Генерирую матрицу гамильтониана ядра H.

*) Диагонализирую матрицу S и рассчитываю сопряженную матрицу X

*) Создаю матрицу плотности P

*) Запуск SCF

**) Генерирую матрицу Фока
**) Строю матрицу F'
**) Диагонализирую матрицу F' и строю матрицу C'
**) Строю матрицу коэффициентов C
**) Пересчитаю матрицу плотности P

 Произведено итераций: 1

**) Генерирую матрицу Фока
**) Строю матрицу F'
**) Диагонализирую матрицу F' и строю матрицу C'
**) Строю матрицу коэффициентов C
**) Пересчитаю матрицу плотности P

 Произведено итераций: 2

**) Генерирую матрицу Фока
**) Строю матрицу F'
**) Диагонализирую матрицу F' и строю матрицу C'
**) Строю матрицу коэффициентов C
**) Пересчитаю матрицу плотности P

 Произведено итераций: 3

**) Генерирую матрицу Фока
**) Строю матрицу F'
**) Диагонализирую матрицу F' и строю матрицу C'
**) Строю матрицу коэффициентов C
**) Пересчитаю матрицу плотности P

 Произведено итераций: 4

**) Генерирую матрицу Фока
**) Строю матрицу F'
**) Диагонализирую матрицу F' и строю матрицу C'
**) Строю матрицу коэффициентов C
**) Пересчитаю матрицу плотности P

 Произведено итераций: 5

**) Генерирую матрицу Фока
**) Строю матрицу F'
**) Диагонализирую матрицу F' и строю матрицу C'
**) Строю матрицу коэффициентов C
**) Пересчитаю матрицу плотности P

 Произведено итераций: 6

**) Генерирую матрицу Фока
**) Строю матрицу F'
**) Диагонализирую матрицу F' и строю матрицу C'
**) Строю матрицу коэффициентов C
**) Пересчитаю матрицу плотности P

 Произведено итераций: 7
--> Сходится

E_electron = ener_elec(orb['P'], orb['H'], orb['F'])
print(E_electron)
E_total = ener_tot(r=1.4632, Z=[1,1], elec=E_electron)
print(E_total)

-1.7974485548087504
-1.1140149845517793

Строю 1D и 2D плот орбиталей:

orbital (orb["C"], r = 1.4632, b1 = GTO["H"], b2 = GTO["H"], b = 3)

orbital2D (orb["C"], orb["X"], orb["F"], r = 1.4632, b1 = GTO["H"], b2 = GTO["H"], b = 3, delta=0.02)