Практикум 5

Оценим молекулярные орбитали молекулы водорода методом методом Хартри-Фока-Рутана: считая МО=ЛКАО, найдем коэфициенты в линейной комбинации.

from pr5.hfscf import *
import numpy as np

def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=False):
    R = [0, r]
    # S -- матрица перекрывания, описывает взаимодействие в базисном наборе АО
    if vbs: print("*) Рассчитать матрицу перекрывания S.")
    s_scf = S(R, b1, b2, b)
    # H -- матрица гамильтониана в ХФ-приближении
    if vbs: print("\n*) Рассчитать гамильтониан H.")
    h_scf = H(R, Z, b1, b2, b)
        
    # Найти матрицу X (диагонализировать матрицу 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)
        # F -- гамильтониан молекулы, переписанный как функция координат одного электрона в приближаниях ХФ
        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 на базис атомных орбиталей: F' = X_adj * F * X)
        if vbs: print("**) Замена базиса F ")
        f_tra = Xa * f_scf * X
        
        # Построить матрицу C' (диагонализировать матрицу F')'
        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)
        
        print("\nВыполнена " + str(iteracion + 1) + ". итерация.\n")
        
        #Рассчитать энергию полученных орбиталей
        energy_orb = ener_orbs(X, f_scf)
        energy_ele = ener_elec(p_temp, h_scf, f_scf)
        energy_tot = ener_tot(r, Z, 0)
        
        # Построить изображения орбиталей
        orbital(c_tra, r, b1, b2, b)
        orbital2D(c_scf, X, f_scf, r, b1, b2, b)
        
        # Проверить сходимость
        if np.linalg.norm(p_temp - p_scf) < 1E-4: # Referencia (7) p. 148
            print("\n\n-->Схождение самосогласованного поля!")
            print('Электронная энергия: {0}, полная энергия: {1}'.format(energy_ele, energy_tot))
            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}

SCF()

Выполнена 1. итерация.

Выполнена 2. итерация.

Выполнена 3. итерация.

Выполнена 4. итерация.

Выполнена 5. итерация.

Выполнена 6. итерация.

Выполнена 7. итерация.

-->Схождение самосогласованного поля!
Электронная энергия: -1.7974485548087507, полная энергия: 0.683433570256971

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

Отметим, что после 2 итерации вид решения качественно не меняется