!wget http://kodomo.fbb.msu.ru/~golovin/qm/hfscf.py

import numpy as np
from hfscf import *

Здесь и далее используем обозначения:

• S -- матрица перекрывания
• F -- матрица Фока
• X -- матрица S, приведенная к диагональному виду
• H -- гамильтониан
• C -- матрица коэффициентов

def SCF(r=1.4632, Z=[1, 1], b1=GTO["H"], b2=GTO["H"], b=3, vbs=False):
    R = [0, r]
    if vbs:
        print("*) Calculate overlap S.")
    s_scf = S(R, b1, b2, b)
    if vbs:
        print("\n*) Calculate Hamiltonian H.")
    h_scf = H(R, Z, b1, b2, b)

    if vbs:
        print("\n*) Diagonalize matrix S to get X.")
    X = diagon(m=s_scf)
    Xa = X.getH()  # Сопряженная матрица

    if vbs:
        print("\n*) Create density matrix P.")
    p_scf = np.matrix(
        [[0, 0], [0, 0]], dtype=np.float64
    )  # Случайно инициализированная матрица плотности

    # Итеративно добиваемся схождения матрицы
    if vbs:
        print("\n*) Calculate SCF.")
    for iteracion in range(50):
        if vbs:
            print("\n**) Calculate Fock matrix")
        g_scf = G(r, p_scf, b1, b2, b)
        f_scf = h_scf + g_scf  # F = H + G

        if vbs:
            print("**) Transform F into F'.")
        f_tra = Xa * f_scf * X  # F' = X_adj * F * X

        if vbs:
            print("**) Diagonalize F, generate C.")
        c_tra = diagon2(m=f_tra)

        if vbs:
            print("**) Calculate coefficient C matrix.")
        c_scf = X * c_tra  # Рассчитать матрицу коэффициентов С = X * C'

        if vbs:
            print("**) Recalculate P matrix.")
        p_temp = P(C=c_scf)  # Пересчитать матрицу P с учетом полученных коэффициентов С

        print("\nIteration #" + str(iteracion + 1) + " is finished.\n")

        if np.linalg.norm(p_temp - p_scf) < 1e-4:  # Проверка сходимости
            print("\n\n-->Finished, it CONVERGES!")
            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-->Finished, no convergence has been achieved")
    return {"S": s_scf, "H": h_scf, "X": X, "F": f_scf, "C": c_scf, "P": p_temp}  # Вывести хоть что-то, пусть и не сошлось

molecule = SCF(vbs=True)
electron_energy = ener_elec(molecule["P"], molecule["H"], molecule["F"])
full_energy = ener_tot(elec=electron_energy)

# Электронная и полная энергия молекулы
print(round(electron_energy, 5), round(full_energy, 5))

*) Calculate overlap S.

*) Calculate Hamiltonian H.

*) Diagonalize matrix S to get X.

*) Create density matrix P.

*) Calculate SCF.

**) Calculate Fock matrix
**) Transform F into F'.
**) Diagonalize F, generate C.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.

Iteration #1 is finished.


**) Calculate Fock matrix
**) Transform F into F'.
**) Diagonalize F, generate C.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.

Iteration #2 is finished.


**) Calculate Fock matrix
**) Transform F into F'.
**) Diagonalize F, generate C.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.

Iteration #3 is finished.


**) Calculate Fock matrix
**) Transform F into F'.
**) Diagonalize F, generate C.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.

Iteration #4 is finished.


**) Calculate Fock matrix
**) Transform F into F'.
**) Diagonalize F, generate C.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.

Iteration #5 is finished.


**) Calculate Fock matrix
**) Transform F into F'.
**) Diagonalize F, generate C.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.

Iteration #6 is finished.


**) Calculate Fock matrix
**) Transform F into F'.
**) Diagonalize F, generate C.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.

Iteration #7 is finished.



-->Finished, it CONVERGES!
-1.79745 -1.11401

orbital(molecule["C"])

orbital2D(molecule["C"],molecule["X"],molecule["F"])