In [13]:
from hfscf import *
import numpy as np

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

In [14]:
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(R, b1, b2, b)
    if vbs: print("\n*)Генерация гамильтониана H.")
    h_scf = H(R, Z, b1, b2, b)        
    # диагонализировать матрицу 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 = 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 и строим матрицу C '
        if vbs: print("**) Диагонализировать матрицу F' и строим 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)
        # Рассчитываем орбитальную, электронную и полную энергию молекулы, строим графики
        orbital_energy = ener_orbs(X, f_scf)
        electronic_energy = ener_elec(p_temp, h_scf, f_scf)
        total_energy = ener_tot(r, Z, 0)
        orbital(c_tra, r, b1, b2, b)
        orbital2D(c_scf, X, f_scf, r, b1, b2, b)
        print("\nЗавершена " + str(iteracion + 1) + " итерация.\n")
        
        # Revisar convergencia
        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-->Самосогласованное поле не сошлось. Проверьте стартовое приближение")
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}
In [15]:
SCF()

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


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


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


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


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


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


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



-->Самосогласованное поле сошлось!
Out[15]:
{'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 ]])}