Pr5: Вычисление параметров молекулы водорода

import hfscf
import numpy as np

def SCF (r=1.4632, Z=[1, 1], b1 = hfscf.GTO["H"], b2 = hfscf.GTO["H"], b = 3, vbs=False):
    R = [0, r]
    
    if vbs: print ("*) Генерируем матрицу перекрывания S.")
    s_scf = hfscf.S(R, b1, b2, b)
    
    if vbs: print ("*) Генерируем гамильтониан H.")
    h_scf = hfscf.H(R, Z, b1, b2, b)
        
    if vbs: print ("*) Проводим диагонализацию матрицы S и считаем диагональную матрицу X.")
    X = hfscf.diagon(s_scf)
    Xa = X.getH()
    
    if vbs: print("*) Считаем матрицу плотности P.")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)
    
    # начало итеративного процесса
    if vbs: print("*) Обрещение к функции SCF.")
    for iteration in range(50):
        if vbs: print("**) Генерируем матрицу Фока: считаем интеграл для 2 электронов.")
        g_scf = hfscf.G(r, p_scf, b1, b2, b)
        f_scf = h_scf + g_scf   
        
        # F' = X_adj * F * X
        if vbs: print("**) Меняем базис матрицы F.")
        f_tra = Xa * f_scf * X
        
        if vbs: print("**) Проводим диагонализацию F' и строим C'.")
        c_tra = hfscf.diagon2(f_tra)
        
        # C = X * C'
        if vbs: print("**) Строим матрицы коэффициентов C.")
        c_scf = X * c_tra
        
        # Вычисляем P из C.
        if vbs: print("**) Заново считаем матрицу плотности P.")
        p_temp = hfscf.P(c_scf)
        
        print("Итерация " + str(iteration + 1) + ". завершилась.")

        # Проверка сходимости
        if np.linalg.norm(p_temp - p_scf) < 1E-4:
            print("\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}

scf = SCF(vbs=True)

*) Генерируем матрицу перекрывания S.
*) Генерируем гамильтониан H.
*) Проводим диагонализацию матрицы S и считаем диагональную матрицу X.
*) Считаем матрицу плотности P.
*) Обрещение к функции SCF.
**) Генерируем матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Проводим диагонализацию F' и строим C'.
**) Строим матрицы коэффициентов C.
**) Заново считаем матрицу плотности P.
Итерация 1. завершилась.
**) Генерируем матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Проводим диагонализацию F' и строим C'.
**) Строим матрицы коэффициентов C.
**) Заново считаем матрицу плотности P.
Итерация 2. завершилась.
**) Генерируем матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Проводим диагонализацию F' и строим C'.
**) Строим матрицы коэффициентов C.
**) Заново считаем матрицу плотности P.
Итерация 3. завершилась.
**) Генерируем матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Проводим диагонализацию F' и строим C'.
**) Строим матрицы коэффициентов C.
**) Заново считаем матрицу плотности P.
Итерация 4. завершилась.
**) Генерируем матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Проводим диагонализацию F' и строим C'.
**) Строим матрицы коэффициентов C.
**) Заново считаем матрицу плотности P.
Итерация 5. завершилась.
**) Генерируем матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Проводим диагонализацию F' и строим C'.
**) Строим матрицы коэффициентов C.
**) Заново считаем матрицу плотности P.
Итерация 6. завершилась.
**) Генерируем матрицу Фока: считаем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Проводим диагонализацию F' и строим C'.
**) Строим матрицы коэффициентов C.
**) Заново считаем матрицу плотности P.
Итерация 7. завершилась.

-->Поле сходится.

# 1D-график орбиталей
hfscf.orbital(scf['C'], r=1.4632, b1=hfscf.GTO["H"], b2=hfscf.GTO["H"], b=3)

# 2D-график орбиталей
hfscf.orbital2D(scf['C'], scf['X'], scf['F'], r=1.4632, b1=hfscf.GTO["H"], b2=hfscf.GTO["H"], b=3, delta=0.02)

orbs_e = hfscf.ener_orbs(scf['X'], scf['F'])  # орбитальная энергия
elec_e = hfscf.ener_elec(scf['P'], scf['H'], scf['F'])  # электронная энергия
total_e = hfscf.ener_tot()  # полная энергия

orbs_e, elec_e, total_e

(array([-0.56461536,  0.6368674 ]), -1.7974485548087507, 0.683433570256971)