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

from hfscf import *

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

def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=False):
    R = [0, r]
    if vbs: print("*) Генерируем матрицу S.")
    s_scf = S()
    if vbs: print("\n*) Генерируем гамильтониан H.")
    h_scf = H()
        
    # Диагонализируем матрицу S и находим диагональную матрицу X
    if vbs: print("\n*) Диагонализируем матрицу S и находим диагональную матрицу X.")
    X = diagon(m=s_scf)
    Xa = X.getH()
    
    # Считаем матрицу плотности P
    if vbs: print("\n*) Считаем матрицу плотности P.")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)  # Referencia (7) p. 148
    
    # Начинаем итерации
    if vbs: print("\n*) Основная часть SCF.")
    for iteracion in range(50):
        # Строим матрицу Фока F
        # F = H + G
        if vbs: print("\n**) Строим матрицу Фока: считаем интеграл для 2 электронов.")
        g_scf = G(p = p_scf)
        f_scf = h_scf + g_scf   # Referencia (7) p. 141 eq. (3.154)
        
        # Строим марицу F'
        # F' = X_adj * F * X
        if vbs: print("**) Меняем базис матрицы F.")
        f_tra = Xa * f_scf * X
        
        # Диагонализируем матрицу F' and генерируем матрицу C'
        if vbs: print("**) Диагонализируем матрицу F' and генерируем матрицу C'.")
        c_tra = diagon2(m=f_tra)
        
        # Строим матрицу C
        # C = X * C'
        if vbs: print("**) Строим матрицу коэффициентов C.")
        c_scf = X * c_tra
        
        # Строим матрицу P из матрицы C
        if vbs: print("**) Пересчитываем матрицы плотности P.")
        p_temp = P(C=c_scf)
        
        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}
        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}

result = SCF(vbs = True)

*) Генерируем матрицу S.

*) Генерируем гамильтониан H.

*) Диагонализируем матрицу S и находим диагональную матрицу X.

*) Считаем матрицу плотности P.

*) Основная часть SCF.

**) Строим матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем матрицу F' and генерируем матрицу C'.
**) Строим матрицу коэффициентов C.
**) Пересчитываем матрицы плотности P.

Итерация 1. закончена.


**) Строим матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем матрицу F' and генерируем матрицу C'.
**) Строим матрицу коэффициентов C.
**) Пересчитываем матрицы плотности P.

Итерация 2. закончена.


**) Строим матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем матрицу F' and генерируем матрицу C'.
**) Строим матрицу коэффициентов C.
**) Пересчитываем матрицы плотности P.

Итерация 3. закончена.


**) Строим матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем матрицу F' and генерируем матрицу C'.
**) Строим матрицу коэффициентов C.
**) Пересчитываем матрицы плотности P.

Итерация 4. закончена.


**) Строим матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем матрицу F' and генерируем матрицу C'.
**) Строим матрицу коэффициентов C.
**) Пересчитываем матрицы плотности P.

Итерация 5. закончена.


**) Строим матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем матрицу F' and генерируем матрицу C'.
**) Строим матрицу коэффициентов C.
**) Пересчитываем матрицы плотности P.

Итерация 6. закончена.


**) Строим матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем матрицу F' and генерируем матрицу C'.
**) Строим матрицу коэффициентов C.
**) Пересчитываем матрицы плотности P.

Итерация 7. закончена.



-->Поле сходится!

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

#Электронная энергия
ener_elec(p, h, f)

-1.7974485548087507

#Полная энергия
ener_tot(elec = ener_elec(p, h, f))

-1.1140149845517797

#Построим 1D plot

orbital(c)

#Построим 2D plot

orbital2D(c, x, f)