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

In [24]:
from hfscf import *
In [29]:
# подсчёт самосогласованного поля = selc-consiestent field = SCF

def SCF (r = 1.4632,     # расстояние между атомами
         Z=[1,1],        # заряды ядер
         b1 = GTO["H"],  # первый атом
         b2 = GTO["H"],  # второй атом
         b = 3,          # число гауссианов в приближении
         vbs=False       # verbose
        ):
    R = [0, r]
    if vbs: print("*) Создание матрицы перекрывания S.")
    s_scf =  S(R, b1, b2, b)
    if vbs: print("\n*) Создание гамильтониана H.")
    h_scf = H(R, Z, b1, b2, b)
        
    # Диагонализируем матрицу 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
        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' = X_adj * F * X
        if vbs: print("**) Создадим матрицу F'.")
        f_tra = Xa * f_scf * X
        
        # Диагонализируем матрицу F и построим матрицу C'
        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")
        
        # Проверим сходимость
        if np.linalg.norm(p_temp - p_scf) < 1E-4: # Referencia (7) p. 148
            print("\n\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-->НЕТ, самосогласованное поле не сошлось!\nПересмотрите исходные предположения.")
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}
In [15]:
result = SCF(vbs=True)

*) Создание матрицы перекрывания S.

*) Создание гамильтониана H.

*) диагонализируем матрицу S и найдём матрицу X.

*) Создание матрицы плотности P.

*) Начнём считать SCF.

**) Создание матрицы Фока: расчёт интегралов электронов.
**) Посчитаем базис F.
**)Диагонализируем матрицу F и построим матрицу C'.
**) Построим матрицу коэффициентов C.
**) Пересчитаем матрицу плотности P.

Завершена 1 итерация.


**) Создание матрицы Фока: расчёт интегралов электронов.
**) Посчитаем базис F.
**)Диагонализируем матрицу F и построим матрицу C'.
**) Построим матрицу коэффициентов C.
**) Пересчитаем матрицу плотности P.

Завершена 2 итерация.


**) Создание матрицы Фока: расчёт интегралов электронов.
**) Посчитаем базис F.
**)Диагонализируем матрицу F и построим матрицу C'.
**) Построим матрицу коэффициентов C.
**) Пересчитаем матрицу плотности P.

Завершена 3 итерация.


**) Создание матрицы Фока: расчёт интегралов электронов.
**) Посчитаем базис F.
**)Диагонализируем матрицу F и построим матрицу C'.
**) Построим матрицу коэффициентов C.
**) Пересчитаем матрицу плотности P.

Завершена 4 итерация.


**) Создание матрицы Фока: расчёт интегралов электронов.
**) Посчитаем базис F.
**)Диагонализируем матрицу F и построим матрицу C'.
**) Построим матрицу коэффициентов C.
**) Пересчитаем матрицу плотности P.

Завершена 5 итерация.


**) Создание матрицы Фока: расчёт интегралов электронов.
**) Посчитаем базис F.
**)Диагонализируем матрицу F и построим матрицу C'.
**) Построим матрицу коэффициентов C.
**) Пересчитаем матрицу плотности P.

Завершена 6 итерация.


**) Создание матрицы Фока: расчёт интегралов электронов.
**) Посчитаем базис F.
**)Диагонализируем матрицу F и построим матрицу C'.
**) Построим матрицу коэффициентов C.
**) Пересчитаем матрицу плотности P.

Завершена 7 итерация.



-->ДА, самосогласованное поле сошлось!
In [16]:
result
Out[16]:
{'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 ]])}

Где

  • S — матрица перекрывания спиновых орбиталей
  • H — матрица гамильтониана
  • X — диагонализированная матрица S
  • F — матрица Фока
  • C — матрица коэффициентов (умножение её на атомные орбитали даст молекулярные орбитали)
  • P — матрица плотности

Все матрицы имеют размер 2х2, потому что рассматриваются два электрона (соответственно, исходно имеются две волновые функции атомных орбиталей, и матрицы воздействующих на них операторов и должны быть 2x2).

In [17]:
orbital(result['C'])
In [25]:
orbital2D(result['C'], result['X'], result['F'])

На графиках видны две орбитали: "связывающая" $\psi_1$ (с отрицательной потенциальной энергией и положительной плотностью между атомами) и "разрыхляющая" $\psi_2$ (с положительной потенциальной энергией и нулевой плотностью посередине между атомами).

In [26]:
print(f"Общая электронная энергия: {(ener_elec(result['P'], result['H'], result['F'])):.2f} атомных единиц (a.u.)")

Общая электронная энергия: -1.80 атомных единиц (a.u.)
In [27]:
print(f"Общая энергия (с учётом взаимодействия ядер): {(ener_tot(elec=ener_elec(result['P'], result['H'], result['F']))):.2f} атомных единиц (a.u.)")

Общая энергия (с учётом взаимодействия ядер): -1.11 атомных единиц (a.u.)