Моделирование молекулы водорода

from hfscf import *
import numpy as np

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

def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=False):
    R = [0, r]
    if vbs: print("*) Calculate overlap matrix S.")
    s_scf = S(R,b1,b2,b) # Рассчитать матрицу перекрывания S
    if vbs: print("\n*) Calculate hamiltonian H.")
    h_scf = H(R, Z, b1, b2, b) # Рассчитать матрицу гамильтониана ядра H
              
    if vbs: print("\n*) Diagonalize matrix S into X.")
    X = diagon(m=s_scf) # Диагонализировать матрицу перекрывания
    Xa = X.getH() # Рассчитать сопряженную матрицу
    
    if vbs: print("\n*) Create density matrix P.")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)  # Создать матрицу плотности Р, наполненную случайными числами
    
    if vbs: print("\n*) Calculate SCF.")
    for iteration in range(50): #Циклически добиваемcя схождения матрицы

        if vbs: print("**) Calculate Fock matrix.")
        g_scf = G(r, p_scf, b1, b2, b) # Рассчитать матрицу одноэлектронных потенциалов Фока
        f_scf = h_scf + g_scf   # Рассчитать матрицу Фока F = H + G
        
        if vbs: print("**) Calculate F'.")
        f_tra = Xa * f_scf * X # Рассчитать матрицу F' = X_adj * F * X
        
        if vbs: print("**) Diagonalize F' and calculate C'.")
        c_tra = diagon2(m=f_tra) # Рассчитать С' как диагонализацию F'
        
        if vbs: print("**) Calculate coefficient C matrix.")
        c_scf = X * c_tra # Рассчитать матрицу коэффициентов С = X * C'
        
        if vbs: print("**) Recalculate P matrix.")
        p_temp = P(C=c_scf) # Пересчитать матрицу P с учетом полученных коэффициентов С
        
        print("Finished iteration № " + str(iteration + 1))
        
        if np.linalg.norm(p_temp - p_scf) < 1E-4: #Проверка сходимости матрицы
            print("\n-->Finished, convergence achieved!")
            return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp} # Если Р сошлась, прервать цикл
            skip
        else:
            p_scf = p_temp
    print("\n-->Finished, convergence not achieved")
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp} #В случае отсутствия сходимости выдть итоговые параметры

Расчитаем все параметры для молекулы водорода, выведем электронную и общую энерии атома водорода

mol=SCF(vbs=True)
E_el=ener_elec(mol['P'],mol['H'],mol['F'])
E_tot=ener_tot(elec=E_el)
print(E_el,"\n",E_tot)

*) Calculate overlap matrix S.

*) Calculate hamiltonian H.

*) Diagonalize matrix S into X.

*) Create density matrix P.

*) Calculate SCF.
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
Finished iteration № 1
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
Finished iteration № 2
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
Finished iteration № 3
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
Finished iteration № 4
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
Finished iteration № 5
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
Finished iteration № 6
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
Finished iteration № 7

-->Finished, convergence achieved!
-1.7974485548087504 
 -1.1140149845517793

Одномерное и двумерное отображение орбиталей

orbital(mol["C"])

orbital2D(mol["C"],mol["X"],mol["F"])