# импортируем функцию Rony J. Letona
from hfscf import *

def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=True):
    R = [0, r]
    if vbs: print("1) Генерируем матрицу перекрываний S")
    s_scf = S(R, b1, b2, b)
    if vbs: print("\n 2) Генерируем гамильтониан H")
    h_scf = H(R, Z, b1, b2, b)
        
    # Диагонализируем матрицу S и находим матрицу X
    if vbs: print("\n 3) Диагонализируем матрицу S и находим матрицу X")
    X = diagon(m=s_scf)
    Xa = X.getH()
    
    # Рассчитываем матрицу расчетной плотности P
    if vbs: print("\n 4) Создаем матрицу расчетной плотности P")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)  # Referencia (7) p. 148
    
    # Начинаем итерационный процесс
    if vbs: print("\n 5) Начинаем с функции SCF")
    for iteracion in range(50):
        # Конструируем матрицу Фока 
        # F = H + G
        if vbs: print("\n 6) Генерируем матрицу Фока: рассчитываем интегралы для 2х электронов")
        g_scf = G(r, p_scf, b1, b2, b)
        f_scf = h_scf + g_scf   
        
        # Генерируем матрицу F'
        # F' = X_adj * F * X
        if vbs: print("n\ 7) Смена базиса F")
        f_tra = Xa * f_scf * X
        
        # Диагонализируем F' и генерируем C'
        if vbs: print("n\ 8) Диагонализируем F' и генерируем C'")
        c_tra = diagon2(m=f_tra)
        
        # Конструируем матрицу C
        # C = X * C'
        if vbs: print("n\ 9) Конструируем матрицу коэффициентов C")
        c_scf = X * c_tra
        
        # Конструируем матрицу Р на основе матрицы С
        if vbs: print("n\ 10) Пересчитываем матрицу расчетной плотности P")
        p_temp = P(C=c_scf)
        
         # рассчитайтываем электронную энергию молекулы и полную энергию молекулы
        energy_electron = ener_elec(p_temp, h_scf, f_scf)
        energy_sum = ener_tot(r, Z, 0)        
        
        # строим графики 1D и 2D
        orbital(c_tra, r, b1, b2, b)
        orbital2D(c_scf, X, f_scf, r, b1, b2, b)
        
        print("\n Номер цикла " + str(iteracion + 1) + "\n")
        
        # Проверяем сходимость
        if np.linalg.norm(p_temp - p_scf) < 1E-4: 
            print("\n\n-->Расчитанное поле сходится")
            return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp, "EE":energy_electron, "ES":energy_sum}
        else:
            p_scf = p_temp
    print("\n\n-->Расчитанное поле НЕ сходится!\n")
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp, "EE":energy_electron, "ES":energy_sum}

orbitals = SCF()

1) Генерируем матрицу перекрываний S

 2) Генерируем гамильтониан H

 3) Диагонализируем матрицу S и находим матрицу X

 4) Создаем матрицу расчетной плотности P

 5) Начинаем с функции SCF

 6) Генерируем матрицу Фока: рассчитываем интегралы для 2х электронов
n\ 7) Смена базиса F
n\ 8) Диагонализируем F' и генерируем C'
n\ 9) Конструируем матрицу коэффициентов C
n\ 10) Пересчитываем матрицу расчетной плотности P

 Номер цикла 1


 6) Генерируем матрицу Фока: рассчитываем интегралы для 2х электронов
n\ 7) Смена базиса F
n\ 8) Диагонализируем F' и генерируем C'
n\ 9) Конструируем матрицу коэффициентов C
n\ 10) Пересчитываем матрицу расчетной плотности P

 Номер цикла 2


 6) Генерируем матрицу Фока: рассчитываем интегралы для 2х электронов
n\ 7) Смена базиса F
n\ 8) Диагонализируем F' и генерируем C'
n\ 9) Конструируем матрицу коэффициентов C
n\ 10) Пересчитываем матрицу расчетной плотности P

 Номер цикла 3


 6) Генерируем матрицу Фока: рассчитываем интегралы для 2х электронов
n\ 7) Смена базиса F
n\ 8) Диагонализируем F' и генерируем C'
n\ 9) Конструируем матрицу коэффициентов C
n\ 10) Пересчитываем матрицу расчетной плотности P

 Номер цикла 4


 6) Генерируем матрицу Фока: рассчитываем интегралы для 2х электронов
n\ 7) Смена базиса F
n\ 8) Диагонализируем F' и генерируем C'
n\ 9) Конструируем матрицу коэффициентов C
n\ 10) Пересчитываем матрицу расчетной плотности P

 Номер цикла 5


 6) Генерируем матрицу Фока: рассчитываем интегралы для 2х электронов
n\ 7) Смена базиса F
n\ 8) Диагонализируем F' и генерируем C'
n\ 9) Конструируем матрицу коэффициентов C
n\ 10) Пересчитываем матрицу расчетной плотности P

 Номер цикла 6


 6) Генерируем матрицу Фока: рассчитываем интегралы для 2х электронов
n\ 7) Смена базиса F
n\ 8) Диагонализируем F' и генерируем C'
n\ 9) Конструируем матрицу коэффициентов C
n\ 10) Пересчитываем матрицу расчетной плотности P

 Номер цикла 7



-->Расчитанное поле сходится

import hfscf
orben = hfscf.ener_orbs(orbitals['X'], orbitals['F'])
elecen = hfscf.ener_elec(orbitals['P'], orbitals['H'], orbitals['F'])
toten = hfscf.ener_tot()
print('Энергии орбиталей: ' + str(orben[0]) + ' и ' + str(orben[1]))
print('Электронная энергия молекулы: ' + str(elecen))
print('Полная энергия молекулы: ' + str(toten))

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