import hsfcf
import numpy as np

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

    R = [0, r] # Задаем центры атомов
    s_scf = hsfcf.S(R, b1, b2, b) # Генерация матрицы перекрывания S
    h_scf = hsfcf.H(R, Z, b1, b2, b) # Генерация гамильтониана H для оператора Фока (core-Hamiltonian operator)
        

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

s=SCF()

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

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

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

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

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

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

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

SCF сошлось! Ура!

hsfcf.orbital(s['C']) #  орбитали в 1 измерении

hsfcf.orbital2D(s['C'], s['X'], s['F']) #и в двух

total = hsfcf.ener_tot()
print('Полная энергия молекулы: ' + str(total))


orbit = hsfcf.ener_orbs(s['X'], s['F'])
print('Энергии орбиталей: ' + str(orbit[0]) + ' и ' + str(orbit[1]))

electr = hsfcf.ener_elec(s['P'], s['H'], s['F'])
print('Электронная энергия молекулы: ' + str(electr))

Полная энергия молекулы: 0.683433570256971
Энергии орбиталей: -0.5646153594583119 и 0.636867398093562
Электронная энергия молекулы: -1.7974485548087507