Вычисление параметров молекулы водорода

import hfscf
import numpy as np

def SCF (r = 1.4632, Z=[1,1], b1 = hfscf.GTO["H"], b2 = hfscf.GTO["H"], b = 3):

    R = [0, r] # Задаем центры атомов
    
    s_scf = hfscf.S(R, b1, b2, b) # Создаем матрицу перекрывания S
    h_scf = hfscf.H(R, Z, b1, b2, b) # Создаем гамильтониан H для оператора Фока
        
    X = hfscf.diagon(m=s_scf) # Переводим матрицу S в диагональную и находим диагональную матрицу Х
    Xa = X.getH() # Получаем транспонированную матрицу Xa
    
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64) #  Создаем матрицу плотности Р
    
    for iteracion in range(50): # Проверяем сходимость SCF
        g_scf = hfscf.G(r, p_scf, b1, b2, b) # Расчет двухэлектронной части матрицы Фока 
        f_scf = h_scf + g_scf # Генерация матрицы Фока
        f_tra = Xa * f_scf * X  # Создание матрицы F' = X_adj * F * X
        c_tra = hfscf.diagon2(m = f_tra) # Диагонализируем матрицу F' и получаем транспонированную матрицу C'
        c_scf = X * c_tra # Построение матрицы коэффициентов C C = X * C'
        p_temp = hfscf.P(C = c_scf) # Из матрицы C строим матрицу Р и пересчитываем матрицы плотности P
        
        print("Итерация " + str(iteracion + 1) + " окончена.\n")
        
        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}

o = SCF()

Итерация 1 окончена.

Итерация 2 окончена.

Итерация 3 окончена.

Итерация 4 окончена.

Итерация 5 окончена.

Итерация 6 окончена.

Итерация 7 окончена.

Сошлось

o # орбитали

{'S': matrix([[0.99999999, 0.63749012],
         [0.63749012, 0.99999999]]), 'H': matrix([[-1.09920375, -0.91954841],
         [-0.91954841, -1.09920375]]), 'X': matrix([[ 0.55258063,  1.17442445],
         [ 0.55258063, -1.17442445]]), 'F': matrix([[-0.34684107, -0.57771139],
         [-0.57771139, -0.34684028]]), 'C': matrix([[ 0.55258114, -1.17442421],
         [ 0.55258013,  1.17442468]]), 'P': matrix([[0.61069182, 0.61069071],
         [0.61069071, 0.6106896 ]])}

Представим в 1D

hfscf.orbital(o['C'])

Представим в 2D

hfscf.orbital2D(o['C'], o['X'], o['F'])