#Вычисление параметров молекулы водорода
#Опишите :
#Добавьте коментарии в код функции
#Опишите объекты S,F,X,H,C

#Рассчитайте:
#Элеткронную энергию молекулы
#Полную энергию молекулы

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

import hfscf as a

import numpy as np

# S - Матрица перекрытий. Описывает наложение базисных функций друг на друга.
# F - Матрица Фока. Описывает поле, действующее на электрон со стороны остальных электронов.
# X - Матрица трансформации. Перевод в ортононрмированный базис.
# H - Core Hamiltionian. Гамильтониан описывает энергию системы.
# C - Матрица коэффициентов. Помогает расчитать молекулярные орбитали из базисных функций (атомных орбиталей).

def SCF (r = 1.4632, Z=[1,1], b1 = a.GTO["H"], b2 = a.GTO["H"], b = 3, vbs=False):
    R = [0, r]
    if vbs: print("*) Создание матрицы перехода S.")
    s_scf = a.S(R, b1, b2, b)
    if vbs: print("\n*) Создание гамильтониана H.")
    h_scf = a.H(R, Z, b1, b2, b)
        
    # Диагонализируем матрицу S и найдём матрицу X
    if vbs: print("\n*) Диагонализация матрицы S и нахождение диагональной матрицы X.")
    X = a.diagon(m=s_scf)
    Xa = X.getH()
    
    # Оценим матрицу плотности P
    if vbs: print("\n*) Создадим матрицу плотности P.")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)
    
    # Начнём итерацию
    if vbs: print("\n*) Вычислим SCF (Self-Consistent Field).")
    for iteracion in range(50):
        # Сделаем матрицу Фока - F
        # F = H + G
        if vbs: print("\n**) Создание матрицы Фока: вычисление \
двухэлектронных интегралов.")
        g_scf = a.G(r, p_scf, b1, b2, 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 = a.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 = a.P(C=c_scf)
        
        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}
        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}

c = SCF()

Итерация 1. завершена.


Итерация 2. завершена.


Итерация 3. завершена.


Итерация 4. завершена.


Итерация 5. завершена.


Итерация 6. завершена.


Итерация 7. завершена.



-->Есть сходимость.

print(c)

{'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 ]])}

#Посчитаем энергии
en1 = a.ener_orbs(c['X'], c['F'])
en2 = a.ener_elec(c['P'], c['H'], c['F'])
full_en = a.ener_tot(r=1.4632, Z=[1,1], elec=0)
print(f'Орбитальная энергия: {en1}')
print(f'Электронная энергия: {en2}')
print(f'Полная энергия молекулы: {full_en}')

Орбитальная энергия: [-0.56461536  0.6368674 ]
Электронная энергия: -1.7974485548087507
Полная энергия молекулы: 0.683433570256971

a.orbital(c['C'], r=1.4632, b1=a.GTO["H"], b2=a.GTO["H"], b=3) #1D плот орбиталей

a.orbital2D(c['C'], c['X'], c['F'], r=1.4632, b1=a.GTO["H"], b2=a.GTO["H"], b=3, delta=0.02) #2D плот орбиталей