Ab initio вычисление молекулы водорода

from hfscf import *

def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=False): # ядерные(R) и 
                                                                               # электронные(r) 
                                                                               # степени свободы
    R = [0, r]
    if vbs: print("cоздаем матрицу перекрываний S") # создаем матрицу перекрывания S
    s_scf = S(R, b1, b2, b)
    if vbs: print("cоздаем гамильтониан Н")    # создаем матрицу Н - гамильтониан 
    h_scf = H(R=R, Z=Z, b1=b1, b2=b2, b=b)
        

    if vbs: print("диагонализируем матрицу S и получаем транспонированную матрицу Х")
    # диагонализируем матрицу S
    X = diagon(m=s_scf)
    # получаем транспонированную матрицу Х
    Xa = X.getH()
    
    # генерируем матрицу расчетной плотности P
    if vbs: print("генерируем матрицу расчетной плотности P")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)  
    
    # в цикле улучшаем матрицу P 
    if vbs: print("начинаем подсчет SCF")
    for iteracion in range(50):
        # генерируем матрицу Фока
        # F = H + G
        if vbs: print("генерируем матрицу Фока")
        g_scf = G(r=r, p=p_scf, b1=b1, b2=b2, b=b)
        f_scf = h_scf + g_scf   
        
        # генерируем матрицу F'
        # F' = X_adj * F * X
        if vbs: print("смена базиса F")
        f_tra = Xa * f_scf * X
        
        # диагонализируем построенную матрицу F' и генерируем матрицу C'
        if vbs: print("диагонализируем построенную матрицу F' и генерируем матрицу C'")
        c_tra = diagon2(m=f_tra)
        
        # генерируем матрицу коэффициентов С, умножая ее на Х
        # C = X * C'
        if vbs: print("**) Cгенерируем матрицу коэффициентов С")
        c_scf = X * c_tra
        
        # создаем новую матрицу расчетной плотности P из матрицы С
        if vbs: print("**) создаем новую матрицу расчетной плотности 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-->самосогласованное поле не сошлось:(")
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp,
            "EE":energy_electron, "ES":energy_sum}

Вызываем полученную функцию, которая в себе содержит вызовы функций, строющих графики и подсчитывающих электронную и полную энергии

SCF()

Цикл номер 1. завершен.

Цикл номер 2. завершен.

Цикл номер 3. завершен.

Цикл номер 4. завершен.

Цикл номер 5. завершен.

Цикл номер 6. завершен.

Цикл номер 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': 0.683433570256971}

Вынимаем из otput посчитанные значения электронной энергии молекулы и полной энергии молекулы

Электронная энергия молекулы равна -1.7974485548087507

Полная энергия молекулы равна 0.683433570256971