In [4]:
%matplotlib inline

import numpy as np
import hfscf as hf

np.set_printoptions(precision=6, suppress=True)

Что означают матрицы S, H, F, X, C¶

S — матрица перекрывания базисных функций.
Она показывает, насколько базисные функции "похожи" друг на друга.

H — одноэлектронный гамильтониан.
Это сумма кинетической энергии электронов и их притяжения к ядрам.

F — матрица Фока.
Это основной оператор в методе Хартри–Фока: [ F = H + G ] где (G) учитывает электрон-электронное взаимодействие.

X — матрица преобразования базиса.
Она нужна, чтобы перейти от неортогонального базиса к ортогональному и удобно диагонализовать задачу.

C — матрица коэффициентов молекулярных орбиталей.
Её столбцы показывают, из каких атомных орбиталей собираются молекулярные орбитали.

Функция SCF¶

Ниже записана функция самосогласованного поля (SCF). Она:

  1. строит матрицы S и H,
  2. получает матрицу X,
  3. начинает с нулевой матрицы плотности P,
  4. на каждой итерации строит F,
  5. пересчитывает коэффициенты C,
  6. обновляет матрицу плотности,
  7. проверяет сходимость.
In [5]:
def SCF(r=1.4632, Z=[1, 1], b1=hf.GTO["H"], b2=hf.GTO["H"], b=3, max_iter=50, tol=1e-4, vbs=True):
    # Координаты двух ядер на оси:
    # первое ядро ставим в 0, второе в точку r
    R = [0, r]

    # S показывает перекрывание атомных базисных функций
    S = hf.S(R=R, b1=b1, b2=b2, b=b)

    # Одноэлектронный гамильтониан H
    # H = T + V
    # T — кинетическая энергия
    # V — притяжение электронов к ядрам
    H = hf.H(R=R, Z=Z, b1=b1, b2=b2, b=b)

    # получили Х
    # X нужна для перехода к ортогональному базису
    X = hf.diagon(S)

    # Сопряжённо-транспонированная матрица X
    Xa = X.getH()

    # Начинаем с нулевой матрицы
    P = np.matrix([[0.0, 0.0],
                   [0.0, 0.0]], dtype=np.float64)

    # Флаг сходимости
    converged = False

    # Запускаем SCF-итерации
    for iteration in range(1, max_iter + 1):

        G = hf.G(r=r, p=P, b1=b1, b2=b2, b=b)
        F = H + G
        F_prime = Xa @ F @ X
        C_prime = hf.diagon2(F_prime)
        C = X @ C_prime
        P_new = hf.P(C)
        
        # Если новая и старая матрицы плотности почти одинаковы,
        # значит SCF сошёлся
        delta = np.linalg.norm(P_new - P)

        if vbs:
            print(f"Итерация {iteration}: ||P_new - P_old|| = {delta:.8f}")

        if delta < tol:
            P = P_new
            converged = True
            break

        # Если не сошлось — идём на следующую итерацию
        P = P_new

    # считаем энергию после сходимости 
    E_elec = hf.ener_elec(P, H, F)      # электронная энергия
    E_tot = hf.ener_tot(r=r, Z=Z, elec=E_elec)  # полная энергия = электронная + отталкивание ядер

    return {
        "S": S,
        "H": H,
        "X": X,
        "F": F,
        "C": C,
        "P": P,
        "E_elec": E_elec,
        "E_tot": E_tot,
        "converged": converged,
        "iterations": iteration
    }
In [6]:
result = SCF()

Итерация 1: ||P_new - P_old|| = 3.36923627
Итерация 2: ||P_new - P_old|| = 3.20813622
Итерация 3: ||P_new - P_old|| = 0.24080679
Итерация 4: ||P_new - P_old|| = 0.02165210
Итерация 5: ||P_new - P_old|| = 0.00195266
Итерация 6: ||P_new - P_old|| = 0.00017610
Итерация 7: ||P_new - P_old|| = 0.00001588
In [9]:
print("SCF сошелся?", result["converged"])
print("Число итераций:", result["iterations"])

SCF сошелся? True
Число итераций: 7
In [10]:
print("S =")
print(result["S"])
print()

print("H =")
print(result["H"])
print()

print("X =")
print(result["X"])
print()

print("F =")
print(result["F"])
print()

print("C =")
print(result["C"])
print()

print("P =")
print(result["P"])

S =
[[1.      0.63749]
 [0.63749 1.     ]]

H =
[[-1.099204 -0.919548]
 [-0.919548 -1.099204]]

X =
[[ 0.552581  1.174424]
 [ 0.552581 -1.174424]]

F =
[[-0.346841 -0.577711]
 [-0.577711 -0.34684 ]]

C =
[[ 0.552581 -1.174424]
 [ 0.55258   1.174425]]

P =
[[0.610692 0.610691]
 [0.610691 0.61069 ]]

Энергии¶

Электронная энергия считается по матрицам P, H и F.

Полная энергия молекулы получается как: Etot = Eelec + Enuc,

где Enuc — отталкивание двух ядер.

In [16]:
print(f"Электронная энергия = {float(result['E_elec']):.6f} Ha")
print(f"Полная энергия     = {float(result['E_tot']):.6f} Ha")

Электронная энергия = -1.797449 Ha
Полная энергия     = -1.114015 Ha

Энергии отдельных орбиталей¶

In [13]:
orbital_energies = hf.ener_orbs(result["X"], result["F"])

print("Энергии орбиталей:")
print(orbital_energies)

Энергии орбиталей:
[-0.564615  0.636867]

1D график орбиталей¶

На этом графике будут показаны:

  • две волновые функции,
  • квадраты этих функций, то есть плотности вероятности.
In [14]:
hf.orbital(result["C"], r=1.4632, b1=hf.GTO["H"], b2=hf.GTO["H"], b=3)
No description has been provided for this image

На первом, 1D-графике показаны две орбитали. Красная кривая — это первая орбиталь. Она не обращается в ноль между двумя ядрами. Это значит, что электронная плотность есть и между атомами. А если плотность есть между атомами, это помогает удерживать их вместе. Поэтому такая орбиталь — связывающая.

Синяя кривая — вторая орбиталь. У неё в середине есть переход через ноль. Это значит, что между атомами появляется область, где орбиталь исчезает. Такая орбиталь уже не помогает связывать атомы, а наоборот, соответствует разрыхлению связи. Поэтому её называют разрыхляющей.

Линии квадратов — это уже не сама волновая функция, а плотность вероятности. То есть они показывают, где электрон скорее всего находится.

Зелёная кривая показывает, что для первой орбитали электронная плотность есть и около ядер, и между ними. Это ещё раз говорит, что орбиталь связывающая.

Чёрная кривая показывает два максимума по бокам, но в центре провал почти до нуля. Значит, между ядрами электрон там почти не находится. Это характерно для разрыхляющей орбитали.

2D график орбиталей¶

Построим двумерную карту плотности для двух молекулярных орбиталей.

In [15]:
hf.orbital2D(result["C"], result["X"], result["F"], r=1.4632, b1=hf.GTO["H"], b2=hf.GTO["H"], b=3, delta=0.03)
No description has been provided for this image

Для Orbital 1 яркие области около двух ядер как бы соединены между собой. Это значит, что плотность распределена между атомами, и такая орбиталь способствует образованию связи.

Для Orbital 2 яркие области разделены тёмной полосой посередине. Эта тёмная область — место, где плотность мала. Значит, между атомами нет хорошего «электронного мостика», и такая орбиталь уже разрыхляет связь.

По энергиям смысл такой же:

первая орбиталь имеет более низкую энергию (−0.564615), значит она выгоднее и заполняется первой;

вторая орбиталь имеет более высокую энергию (0.636867), поэтому она менее выгодна.

Вывод: первая орбиталь — связывающая: она помогает удерживать два атома водорода вместе. Вторая орбиталь — разрыхляющая. Если бы электроны сидели на ней, связь ослаблялась бы. Графики получились, по моему мнению, правильные и хорошо соответствуют молекуле H₂.