Ab initio вычисление молекулы водорода

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

# импортируем функцию Rony J. Letona
from hfscf import *

# Импортируем функцию расчёта орбиталей на основе подхода Рутана Холла
# ссылки на Теория (ref 7) : https://chemistlibrary.files.wordpress.com/2015/02/modern-quantum-chemistry.pdf


def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=True): 
    R = [0, r]
    # здесь R - ядерные, r - электронные степени свободы
    if vbs: print("\n1). Cоздаем матрицу перекрываний S")
    s_scf = S(R, b1, b2, b) # понятно из параметров функции
    if vbs: print("\n2). Cоздаем гамильтониан Н")    # матрицу Н
    h_scf = H(R=R, Z=Z, b1=b1, b2=b2, b=b)
        
    if vbs: print("\n3).Диагонализируем матрицу S и получаем транспонированную матрицу Х")
    # сначала работаем с матрицей S, а потом транспонируем Х
    X = diagon(m=s_scf)
    Xa = X.getH()
    
    if vbs: print("\n4). Генерируем матрицу расчетной плотности P") # p. 148
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)  
    
    # в цикле улучшаем матрицу P 
    if vbs: print("\n5). Начинаем подсчет SCF")
    for iteracion in range(50): # за 50 попыток должно сойтись
        # F = H + G
        if vbs: print("\n6). Генерируем матрицу Фока") # p. 141 eq. (3.154)
        g_scf = G(r=r, p=p_scf, b1=b1, b2=b2, b=b)
        f_scf = h_scf + g_scf   
        
        # генерируем матрицу F'
        # F' = X_adj * F * X
        if vbs: print("\n7). Смена базиса F")
        f_tra = Xa * f_scf * X
        
        # диагонализируем построенную матрицу F' и генерируем матрицу C' (ошибка в коде задания- у нас уже F', а не F)
        if vbs: print("\n8). Диагонализируем построенную матрицу F' и генерируем матрицу C'")
        c_tra = diagon2(m=f_tra)
        
        # построеним матрицу орбитальных коэффициентов C = X * C'
        if vbs: print("\n9). Cгенерируем матрицу коэффициентов С")
        c_scf = X * c_tra
        
        # создаем новую матрицу расчетной плотности P из матрицы С
        if vbs: print("\n10). создаем новую матрицу расчетной плотности P")
        p_temp = P(C=c_scf)
        
        
        # рассчитайтываем электронную энергию молекулы и полную энергию молекулы
        energy_electron = ener_elec(p_temp, h_scf, f_scf)
        energy_sum = ener_tot(r, Z, 0)        
        
        # строим графики 1D и 2D
        orbital(c_tra, r, b1, b2, b)
        orbital2D(c_scf, X, f_scf, r, b1, b2, b)
        
        
        print("\n Цикл " + str(iteracion + 1) + ". OK.\n") # выводим номер цикла
        
        # проверяем сходимость
        if np.linalg.norm(p_temp - p_scf) < 1E-4:  # p. 148
            
            print("\n\n-->самосогласованное поле сошлось")
            return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp, 
                    "EE":energy_electron, "ES":energy_sum}
        else:
            p_scf = p_temp
    print("\n\n-->самосогласованное поле не сошлось")
    print("\n с параметрами:\n")
    
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp,
            "EE":energy_electron, "ES":energy_sum}

SCF()

1). Cоздаем матрицу перекрываний S

2). Cоздаем гамильтониан Н

3).Диагонализируем матрицу S и получаем транспонированную матрицу Х

4). Генерируем матрицу расчетной плотности P

5). Начинаем подсчет SCF

6). Генерируем матрицу Фока

7). Смена базиса F

8). Диагонализируем построенную матрицу F' и генерируем матрицу C'

9). Cгенерируем матрицу коэффициентов С

10). создаем новую матрицу расчетной плотности P

 Цикл 1. OK.


6). Генерируем матрицу Фока

7). Смена базиса F

8). Диагонализируем построенную матрицу F' и генерируем матрицу C'

9). Cгенерируем матрицу коэффициентов С

10). создаем новую матрицу расчетной плотности P

 Цикл 2. OK.


6). Генерируем матрицу Фока

7). Смена базиса F

8). Диагонализируем построенную матрицу F' и генерируем матрицу C'

9). Cгенерируем матрицу коэффициентов С

10). создаем новую матрицу расчетной плотности P

 Цикл 3. OK.


6). Генерируем матрицу Фока

7). Смена базиса F

8). Диагонализируем построенную матрицу F' и генерируем матрицу C'

9). Cгенерируем матрицу коэффициентов С

10). создаем новую матрицу расчетной плотности P

 Цикл 4. OK.


6). Генерируем матрицу Фока

7). Смена базиса F

8). Диагонализируем построенную матрицу F' и генерируем матрицу C'

9). Cгенерируем матрицу коэффициентов С

10). создаем новую матрицу расчетной плотности P

 Цикл 5. OK.


6). Генерируем матрицу Фока

7). Смена базиса F

8). Диагонализируем построенную матрицу F' и генерируем матрицу C'

9). Cгенерируем матрицу коэффициентов С

10). создаем новую матрицу расчетной плотности P

 Цикл 6. OK.


6). Генерируем матрицу Фока

7). Смена базиса F

8). Диагонализируем построенную матрицу F' и генерируем матрицу C'

9). Cгенерируем матрицу коэффициентов С

10). создаем новую матрицу расчетной плотности P

 Цикл 7. OK.



-->самосогласованное поле сошлось

{'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 ]]),
 'EE': -1.7974485548087507,
 'ES': 0.683433570256971}

Электронная энергия финальной молекулы -1.797. Полная энергия молекулы 0.683