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

from hfscf import *
import numpy as np

def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 6, vbs=True):
  
    # b=6 -- аппроксимация по 6 гауссианам каждого атома
    R = [0, r] # расстояние между ядрами
    
    # расчитываем интеграл перекрывания
    if vbs: 
        print("*) Генерация матрицы S")
        
    s_scf = S(R,b1,b2,b)
    
    
    # расчитываем одноэлектронный гамильтониан
    if vbs: 
        print("\n*) Генерация гамильтониана H.")
        
    h_scf = H(R,Z,b1,b1,b)
    
        
    # Diagonalizar matriz S y hallar matriz X
    if vbs: 
        print("\n*) Diagonalizando matriz S y hallando matriz diagonal X.")
    X = diagon(m=s_scf)
    Xa = X.getH()
    
  
    # Расчитываем матрицу электронной плотности
    if vbs: 
        print("\n*) Создаю матрицу электронно плотности P.")
        
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)  # Referencia (7) p. 148
    
    #Повторям вычисления до сходимости или до 50 интераций
    
    if vbs: 
        print("\n*) Comenzando con el SCF.")
        
    for iteracion in range(50):
        # Конструирование матрицы Фока
        # F = H + G
        if vbs: 
            print("\n**) Generando la matriz de Fock: calculando integrales de dos electrones.")
        g_scf = G(r,p_scf,b1,b2,b)
        f_scf = h_scf + g_scf   # Referencia (7) p. 141 eq. (3.154)
        
        #Конструируем фокиан из одноэлектронного гамильтониана и обменного и кулоновского операторов
        
        # F' = X_adj * F * X
        if vbs: print("**) Cambiando la base de F.")
        f_tra = Xa * f_scf * X
        
        # Diagonalizar matriz F y constuir matriz C'
        if vbs: print("**) Diagonalizando F' y generando C'.")
        c_tra = diagon2(m=f_tra)
        
        # C = X * C'
        if vbs: 
            print("**) Конструирую матрицу коэффициентов C.")
        c_scf = X * c_tra
        
        #Расчитываем матрицу коэффициентов из уравения Хартри-Фока
        
        if vbs: 
            print("**) Recalculando matriz de densidad P.")
        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}

HF = SCF()

*) Generando matriz de traslape S.

*) Generando hamiltoniano H.

*) Diagonalizando matriz S y hallando matriz diagonal X.

*) Creando matriz de densidad P.

*) Comenzando con el SCF.

**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 1. iteracion.


**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 2. iteracion.


**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 3. iteracion.


**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 4. iteracion.


**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 5. iteracion.


**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 6. iteracion.


**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 7. iteracion.



-->El campo autoconsistente SI ha convergido!

S -- матрица перекрывания атомных орбиталей;

H -- одноэлектронный гамильтониан;

X -- ортогонализированная матрица S;

F -- фокиан;

P -- матрица электронной плотсности;

C -- матрица коэффициентов

c = HF['C']
x = HF['X']
f = HF['F']
p = HF['P']
h = HF['H']

Молекулярные орбитали в 2D и 3D

orbital(c)

orbital2D(c,x,f)

Орбиталь 1 -- связывающая, ее энергия отрицательная, орбиталь 2 -- разрыхляющая, ее энергия положительная.

Энергия электронов:

ener_elec(p, h, f)

-1.8139359654437917

Полная энергия:

ener_tot(elec=ener_elec(p,h,f))

-1.1305023951868205