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

from hfscf import *

функция ниже дает ссылку сюда: https://chemistlibrary.files.wordpress.com/2015/02/modern-quantum-chemistry.pdf

#Определим расчет орбиталей на базе подхода Рутхана-Хола
def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=True):
    R = [0, r] #r - электронные степени свободы, 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 (Xa - сопряженно-транспонированная матрица)
    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

    # Начало итерационного процесса улучшения матрицы P
    if vbs: print("\n 5) Обратимся к функции SCF.")
    for iteracion in range(50):
        # Конструирование матрицы Фока 
        # F = H + G
        if vbs: print("\n 6) Генерирование матрицы Фока: вычисление интегралов для двух электронов.")
        g_scf = G(r, p_scf, b1, b2, b)
        f_scf = h_scf + g_scf # Referencia (7) p. 141 eq. (3.154)
 
        # Генерирование матрицы 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)
        
         # Рассчет электронной энергии молекулы и полной энергии молекулы
        E_elec = ener_elec(p_temp, h_scf, f_scf)
        E_tot = ener_tot(r, Z, E_elec) 
        print('Электронная энергия молекулы ', E_elec)
        print('Полная энергия молекулы ', E_tot)
        
        # Построение графиков 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: # Referencia (7) p. 148
            print("\n\n-->Расчитанное поле сходится!")
            return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp, "EE":E_elec, "ES":E_tot}
        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":E_elec, "ES":E_tot}

SCF()

1) Генерирование матрицы перекрывания S.

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

 3) Диагонализирование матрицы S и поиск матрицы X.

 4) Рассчет матрицы расчетной плотности P.

 5) Обратимся к функции SCF.

 6) Генерирование матрицы Фока: вычисление интегралов для двух электронов.

 7) Смена базиса F.

 8) Диагонализирование F' и генерирование C'.

 9) Конструирование матрицы коэффициентов C.

 10) Перерасчет матрицы расчетной плотности P.
Электронная энергия молекулы  -1.7284206318164699
Полная энергия молекулы  -1.0449870615594987

 Завершение цикла 1.


 6) Генерирование матрицы Фока: вычисление интегралов для двух электронов.

 7) Смена базиса F.

 8) Диагонализирование F' и генерирование C'.

 9) Конструирование матрицы коэффициентов C.

 10) Перерасчет матрицы расчетной плотности P.
Электронная энергия молекулы  -1.5654299820991575
Полная энергия молекулы  -0.8819964118421865

 Завершение цикла 2.


 6) Генерирование матрицы Фока: вычисление интегралов для двух электронов.

 7) Смена базиса F.

 8) Диагонализирование F' и генерирование C'.

 9) Конструирование матрицы коэффициентов C.

 10) Перерасчет матрицы расчетной плотности P.
Электронная энергия молекулы  -1.795075334555484
Полная энергия молекулы  -1.111641764298513

 Завершение цикла 3.


 6) Генерирование матрицы Фока: вычисление интегралов для двух электронов.

 7) Смена базиса F.

 8) Диагонализирование F' и генерирование C'.

 9) Конструирование матрицы коэффициентов C.

 10) Перерасчет матрицы расчетной плотности P.
Электронная энергия молекулы  -1.7974292055137375
Полная энергия молекулы  -1.1139956352567664

 Завершение цикла 4.


 6) Генерирование матрицы Фока: вычисление интегралов для двух электронов.

 7) Смена базиса F.

 8) Диагонализирование F' и генерирование C'.

 9) Конструирование матрицы коэффициентов C.

 10) Перерасчет матрицы расчетной плотности P.
Электронная энергия молекулы  -1.7974483974402236
Полная энергия молекулы  -1.1140148271832526

 Завершение цикла 5.


 6) Генерирование матрицы Фока: вычисление интегралов для двух электронов.

 7) Смена базиса F.

 8) Диагонализирование F' и генерирование C'.

 9) Конструирование матрицы коэффициентов C.

 10) Перерасчет матрицы расчетной плотности P.
Электронная энергия молекулы  -1.7974485535391342
Полная энергия молекулы  -1.1140149832821633

 Завершение цикла 6.


 6) Генерирование матрицы Фока: вычисление интегралов для двух электронов.

 7) Смена базиса F.

 8) Диагонализирование F' и генерирование C'.

 9) Конструирование матрицы коэффициентов C.

 10) Перерасчет матрицы расчетной плотности P.
Электронная энергия молекулы  -1.7974485548087507
Полная энергия молекулы  -1.1140149845517797

 Завершение цикла 7.



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

{'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 ]]),
 'EE': -1.7974485548087507,
 'ES': -1.1140149845517797}