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

from hfscf import *

Функция расчёта орбиталей на основе подхода Рутана Холла:

# S - матрица перекрывания
# F - матрица Фока
# X - диагонализированная матрица S
# H - гамильтониан ядра
# C - матрица коэффициентов

def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3):
    R = [0, r]
    #генерируем матрица перекрывания S
    s_scf = S(R=R, b1=b1, b2=b2, b=b)       # Referencia (7) p. 412 eq. (A.9)
    #генерируем гамильтониан ядра H 
    h_scf = H(R=R, Z=Z, b1=b1, b2=b2, b=b)  # Referencia (7) p. 141 eq. (3.153)
        
    #диагонализируем матрицу S
    X = diagon(m=s_scf)
    #находим сопряженную транспонированную матрицу
    Xa = X.getH()
    
    #создаем матрицу плотности P
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)  # Referencia (7) p. 148
    
    #пробуем итеративно улучшить матрицу плотности
    for iteracion in range(50):
        # строим матрицу Фока F
        # F = H + G
        g_scf = G(r=r, p=p_scf, b1=b1, b2=b2, b=b)
        f_scf = h_scf + g_scf     # Referencia (7) p. 141 eq. (3.154)
        
        # строим матрицу F' 
        # F' = X_adj * F * X
        f_tra = Xa * f_scf * X
        
        # диагонализируем матрицу F' и строим матрицу C'
        c_tra = diagon2(m=f_tra)
        
        # строим матрицу орбитальных коэффициентов C
        # C = X * C'
        c_scf = X * c_tra
        
        # считаем новую матрицу плотности P из матрицы C
        p_temp = P(C=c_scf)
        
        #print("\nConcluida la " + str(iteracion + 1) + ". iteracion.\n")
         
        # проверяем сходимость
        if np.linalg.norm(p_temp - p_scf) < 1E-4: # Referencia (7) p. 148
            #print("\n\n-->El campo autoconsistente SI ha convergido!")
            return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}
        else:
            p_scf = p_temp
    print("\n\n-->El campo autoconsistente NO ha convergido!\nRevisar supuestos.")
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}

r = 1.4632
Z=[1,1]
b1 = GTO["H"]
b2 = GTO["H"]
b = 3

res = SCF(r, Z, b1, b2, b)

Рассчитаем электронную энергию молекулы (E_elec) и полную энергию молекулы (E_tot)

E_elec = ener_elec(res['P'], res['H'], res['F'])
E_tot = ener_tot(r, Z, E_elec)
print(E_elec)
print(E_tot)

-1.7974485548087507
-1.1140149845517797

Построим 1D и 2D плот орбиталей

orbital(res['C'], r, b1=b1, b2=b2, b=b)
orbital2D(res['C'], res['X'], res['F'], r, b1=b1, b2=b2, b=b)