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_scf = S(R, b1, b2, b)
    if vbs: print("\n*) Построение гамильтониана")
    h_scf = H(R, Z, b1, b2, b)
        
    # Диагонализация S и получение матриц X и Xa (присоединенной матрицы X)
    if vbs: print("\n*) Диагонализация матрицы перекрывания")
    X = diagon(m=s_scf)
    Xa = X.getH()
    
    # Создание матрицы P
    if vbs: print("\n*) Построение исходной матрицы плотности")
    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 = H + G
        if vbs: print("\n**) Вычисление интегралов электронов, построение матрицы Фока")
        g_scf = G(r, p_scf, b1, b2, b)
        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'
        if vbs: print("**) Диагонализация матрицы F' и получение C'")
        c_tra = diagon2(m=f_tra)
        
        # Расчет С
        # C = X * C'
        if vbs: print("**) Вычисление матрицы коэффициентов")
        c_scf = X * c_tra
        
        # Расчет P из матрицы коэффициентов C
        if vbs: print("**) Пересчет матрицы плотности")
        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--> Сходимость SCF достигнута")
            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--> Нет сходимости SCF\nПересмотрите предположения")
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}

h2_scf = SCF(vbs=True)

*) Построение матрицы перекрывания

*) Построение гамильтониана

*) Диагонализация матрицы перекрывания

*) Построение исходной матрицы плотности

*) Старт итераций SCF

**) Вычисление интегралов электронов, построение матрицы Фока
**) Построение матрицы F'
**) Диагонализация матрицы F' и получение C'
**) Вычисление матрицы коэффициентов
**) Пересчет матрицы плотности

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


**) Вычисление интегралов электронов, построение матрицы Фока
**) Построение матрицы F'
**) Диагонализация матрицы F' и получение C'
**) Вычисление матрицы коэффициентов
**) Пересчет матрицы плотности

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


**) Вычисление интегралов электронов, построение матрицы Фока
**) Построение матрицы F'
**) Диагонализация матрицы F' и получение C'
**) Вычисление матрицы коэффициентов
**) Пересчет матрицы плотности

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


**) Вычисление интегралов электронов, построение матрицы Фока
**) Построение матрицы F'
**) Диагонализация матрицы F' и получение C'
**) Вычисление матрицы коэффициентов
**) Пересчет матрицы плотности

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


**) Вычисление интегралов электронов, построение матрицы Фока
**) Построение матрицы F'
**) Диагонализация матрицы F' и получение C'
**) Вычисление матрицы коэффициентов
**) Пересчет матрицы плотности

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


**) Вычисление интегралов электронов, построение матрицы Фока
**) Построение матрицы F'
**) Диагонализация матрицы F' и получение C'
**) Вычисление матрицы коэффициентов
**) Пересчет матрицы плотности

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


**) Вычисление интегралов электронов, построение матрицы Фока
**) Построение матрицы F'
**) Диагонализация матрицы F' и получение C'
**) Вычисление матрицы коэффициентов
**) Пересчет матрицы плотности

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


--> Сходимость SCF достигнута

h2_scf

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

elec = ener_elec(h2_scf['P'], h2_scf['H'], h2_scf['F'])
total = ener_tot(elec=elec)
print('Электронная энергия молекулы: {:.2f} a. u. ({:.2f} ккал/моль), общая энергия: {:.2f} a. u. ({:.2f} ккал/моль)'.format(elec, elec * 627.51, total, total * 627.51))

Электронная энергия молекулы: -1.80 a. u. (-1127.92 ккал/моль), общая энергия: -1.11 a. u. (-699.06 ккал/моль)

orbital(h2_scf['C'])

orbital2D(h2_scf['C'], h2_scf['X'], h2_scf['F'])

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

Расчет SCF¶

Визуализация орбиталей¶