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

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 -- electron DoF, R -- nuclear DoF.
    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)
        
    # Диагонализируем матрицу S и находим диагональную матрицу X.
    if vbs: print ("*) Диагонализируем матрицу S и находим диагональную матрицу X.")
    X = hfscf.diagon(s_scf)
    Xa = X.getH()
    
    # Считаем матрицу плотности P.
    if vbs: print("*) Считаем матрицу плотности P.")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)
    
    # Начинаем итерации (Расчет F' зависит от собственного решения).
    if vbs: print("*) Основная часть SCF.")
    
    # 50 итераций.
    for iteration in range(50):
        
        # Строим матрицу Фока.
        # F = H + G
        if vbs: print("**) Строим матрицу Фока: вычисляем интеграл для 2 электронов.")
        g_scf = hfscf.G(r, p_scf, b1, b2, b)
        f_scf = h_scf + g_scf   
        
        # Строим матрицу F'.
        # F' = X_adj * F * X
        if vbs: print("**) Меняем базис матрицы F.")
        f_tra = Xa * f_scf * X
        
        # Диагонализируем F' and генерируем C'.
        if vbs: print("**) Диагонализируем F' and генерируем C'.")
        c_tra = hfscf.diagon2(f_tra)
        
        # Строим C.
        # 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' and генерируем C'.
**) Строим матрицу коэффициентов C.
**) Перерасчет матрицы плотности P.
Итерация 1. закончена.
**) Строим матрицу Фока: вычисляем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем F' and генерируем C'.
**) Строим матрицу коэффициентов C.
**) Перерасчет матрицы плотности P.
Итерация 2. закончена.
**) Строим матрицу Фока: вычисляем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем F' and генерируем C'.
**) Строим матрицу коэффициентов C.
**) Перерасчет матрицы плотности P.
Итерация 3. закончена.
**) Строим матрицу Фока: вычисляем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем F' and генерируем C'.
**) Строим матрицу коэффициентов C.
**) Перерасчет матрицы плотности P.
Итерация 4. закончена.
**) Строим матрицу Фока: вычисляем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем F' and генерируем C'.
**) Строим матрицу коэффициентов C.
**) Перерасчет матрицы плотности P.
Итерация 5. закончена.
**) Строим матрицу Фока: вычисляем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем F' and генерируем C'.
**) Строим матрицу коэффициентов C.
**) Перерасчет матрицы плотности P.
Итерация 6. закончена.
**) Строим матрицу Фока: вычисляем интеграл для 2 электронов.
**) Меняем базис матрицы F.
**) Диагонализируем F' and генерируем C'.
**) Строим матрицу коэффициентов C.
**) Перерасчет матрицы плотности P.
Итерация 7. закончена.

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

В модуле есть функции orbital, orbital2D, ener_orbs, ener_elec, ener_tot:

hfscf.orbital(scf['C'], r=1.4632, b1=hfscf.GTO["H"], b2=hfscf.GTO["H"], b=3)

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

array([-0.56461536,  0.6368674 ])

elec_e

-1.7974485548087507

total_e

0.683433570256971