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

Загрузим функцию расчёта орбиталей на основе подхода Рутана Холла. Данная функция включает в себя: S - матрицу перекрывания; F - матрицу Фока; X - диагонализированную матрица S; H - гамильтониан ядра; C - матрицу коэффициентов

import hfscf
import numpy as np

def SCF (r = 1.4632, Z = [1, 1], b1 = hfscf.GTO["H"], b2 = hfscf.GTO["H"], b = 3 , vbs = False):
    
    # r - электронные степени свободы, R - ядерные степени свободы
    R = [0, r]
    
    if vbs: print("Генерируем матрицу перекрываний S.")
    s_scf = hfscf.S(R, b1, b2, b)
    
    # H matrix - Гамильтониан для оператора Фрока
    # H = Ke + Pe, где Ke - кинетическая энергия электрона, Pe - потенциальная энергия электрона
    if vbs: print ("\n*) Генерируем Гамильтониан H.")
    h_scf = hfscf.H(R, Z, b1, b2, b)
        
    if vbs: print("Диагонализация матрицы S и поиск диагональной матрицы X.")
    X = hfscf.diagon(s_scf)
    Xa = X.getH()
    
    if vbs: print("Создаем матрицу плотности P.")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)
    
    # Since F' depends on it’s own solution (through the orbitals), the process
    # of minimizing Hartree-Fock energy must be done iteratively.
    
    # Запуск итеративных процессов ( так как F' зависит от своего собственного решения, процесс 
    # минимизации энергии Харти-Фока нужно делать итеративно)
    if vbs: print("Начинаем считать SCF.")  
    for iteration in range(50):
        # Строим матрицу Фока F 
        # F = H + G
        if vbs: print("Генерация матрицы Фока: вычисление интегралов двух электронов.")
        g_scf = hfscf.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
        
        if vbs: print("Диагонализация F' и генерация C'.")
        c_tra = hfscf.diagon2(f_tra)
        
        #C = X * C'
        if vbs: print("Построение матрицы коэффициентов C.")
        c_scf = X * c_tra
        
        # Строим матрицу P на основе матрицы C
        if vbs: print("Расчет матрицы плотности P.")
        p_temp = hfscf.P(c_scf)
        
        # Вычисляем орбитальную, электронную и полную энергию молекул 

        orbital_energy = hfscf.ener_orbs(X, f_scf)
        electronic_energy = hfscf.ener_elec(p_temp, h_scf, f_scf)
        total_energy = hfscf.ener_tot(r, Z, 0)
        
         # Построим графики
        hfscf.orbital(c_tra, r, b1, b2, b)
        hfscf.orbital2D(c_scf, X, f_scf, r, b1, b2, b, 0.02)
        
        print("\n Итерация " + str(iteration + 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,"OE":orbital_energy,
                   "EE":electronic_energy, "TE":total_energy}
        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,"OE":orbital_energy,
            "EE":electronic_energy, "TE":total_energy}

result = SCF()

 Итерация 1. завершена.

 Итерация 2. завершена.

 Итерация 3. завершена.

 Итерация 4. завершена.

 Итерация 5. завершена.

 Итерация 6. завершена.

 Итерация 7. завершена.



-->Нашлось самосогласованное поле!

print("Орбитальная энергия равна " + str(result.get("OE")[0]) + " and " + str(result.get("OE")[1]))
print("Электроная энергия равна " + str(result.get("EE")))
print("Полная энергия равна " + str(result.get("TE")))

Орбитальная энергия равна -0.564615359458 and 0.636867398094
Электроная энергия равна -1.79744855481
Полная энергия равна 0.683433570257