Практикум 5

Для начала допишем функцию для расчёта орбиталей на основе подхода Рутана Холла (попутно поясняя каждый шаг и описывая используемые матрицы, также мы используем скрипт hsfcf.py с прописанными вспомогательными функциями).

import hsfcf
import numpy as np

def SCF (r = 1.4632, Z=[1,1], b1 = hsfcf.GTO["H"], b2 = hsfcf.GTO["H"], b = 3):
# Задаем центры атомов
    R = [0, r]
    
# Генерация матрицы перекрывания S
    s_scf = hsfcf.S(R, b1, b2, b)
    
# Генерация гамильтониана H для оператора Фока (core-Hamiltonian operator)
    h_scf = hsfcf.H(R, Z, b1, b2, b)
        
# Диагонализация матрицы S и поиск диагональной матрицы X (X = S^(-1/2))
# Xa - сопряженно-транспонированная матрица
    X = hsfcf.diagon(m=s_scf)
    Xa = X.getH()
    
# Создание матрицы плотности P
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)
    
# Начинаем проверку сходимости SCF (итеративно заполняем матрицу F, диагонализируем, проверяем разницу)
    for iteracion in range(50):
        
        # Расчет двух электронной части матрицы Фока (G)
        # Генерация матрицы Фока как:
        # F = H + G
        g_scf = hsfcf.G(r, p_scf, b1, b2, b)
        f_scf = h_scf + g_scf
        
        # Создание матрицы F'
        # F' = X_adj * F * X
        f_tra = Xa * f_scf * X
        
        # Диагонализация F' и генерация C'
        c_tra = hsfcf.diagon2(m=f_tra)
        
        # C = X * C'
        # Построение матрицы коэффициентов C
        c_scf = X * c_tra
        
        # Из матрицы C строим матрицу Р
        # Пересчет матрицы плотности P
        p_temp = hsfcf.P(C=c_scf)
        
        print("Завершил " + str(iteracion + 1) + " итерацию.\n")
        
        # Проверка сходимости
        if np.linalg.norm(p_temp - p_scf) < 1E-4:
            print("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("SCF не сошлось.")
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}

Применяем функцию (с параметрами по умолчанию)

orbitals = SCF()

Завершил 1 итерацию.

Завершил 2 итерацию.

Завершил 3 итерацию.

Завершил 4 итерацию.

Завершил 5 итерацию.

Завершил 6 итерацию.

Завершил 7 итерацию.

SCF сошлось.

На выходе - посчитанные матрицы

orbitals

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

Нарисуем орбитали в 1 измерении

hsfcf.orbital(orbitals['C'])

И в 2 измерениях

hsfcf.orbital2D(orbitals['C'], orbitals['X'], orbitals['F'])

Видны 2 орбитали (связывающая слева и разрыхляющая справа

orben = hsfcf.ener_orbs(orbitals['X'], orbitals['F'])
elecen = hsfcf.ener_elec(orbitals['P'], orbitals['H'], orbitals['F'])
toten = hsfcf.ener_tot()
print('Энергии орбиталей: ' + str(orben[0]) + ' и ' + str(orben[1]))
print('Электронная энергия молекулы: ' + str(elecen))
print('Полная энергия молекулы: ' + str(toten))

Энергии орбиталей: -0.5646153594583118 и 0.6368673980935621
Электронная энергия молекулы: -1.7974485548087507
Полная энергия молекулы: 0.683433570256971