In [3]:
from hfscf import *
import hfscf

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

Завершена 1. итерация.
Завершена 2. итерация.
Завершена 3. итерация.
Завершена 4. итерация.
Завершена 5. итерация.
Завершена 6. итерация.
Завершена 7. итерация.

-->SCF сошлось

S - матрица перекрывания, интегралы перекрывания базисных функций

F - матрица Фока (фокиан), сумма гамильтониана H и двухэлектронного отталкивания G

X - диагональная матрица перехода

H - одноэлектронный гамильтониан для оператора Фока

C - матрица коэффициентов для комбинации базисных функций и получения молекулярных орбиталей

Сгенерированны орбитали:

In [5]:
o
Out[5]:
{'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 ]])}

D1 график¶

In [6]:
hfscf.orbital(o['C'])

D2 график¶

In [8]:
hfscf.orbital2D(o['C'], o['X'], o['F'])

Связывающая (1) и разрыхляющая (2) орбитали.

Рассчет энергии молекулы¶

In [10]:
orb = hfscf.ener_orbs(o['X'], o['F'])
el_e = hfscf.ener_elec(o['P'], o['H'], o['F'])
tot_e = hfscf.ener_tot(elec=el_e)

print('Энергии связывающей орбиталей:\t\t\t' + str(orb[0]))
print('Энергии разрыхляющей орбиталей:\t\t\t' + str(orb[1]))
print('Электронная энергия молекулы:\t\t\t' + str(el_e))
print('Полная энергия молекулы:\t\t\t' + str(tot_e))

Энергии связывающей орбиталей:			-0.5646153594583119
Энергии разрыхляющей орбиталей:			0.6368673980935621
Электронная энергия молекулы:			-1.7974485548087507
Полная энергия молекулы:			-1.1140149845517797