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

Импортировать функции из файла hfscf.py

from hfscf import *

Задать функцию расчета орбиталей

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 matrix S.")
    s_scf = S(R,b1,b2,b) # Рассчитать матрицу перекрывания S
    if vbs: print("\n*) Calculate hamiltonian H.")
    h_scf = H(R, Z, b1, b2, b) # Рассчитать матрицу гамильтониана ядра H
              
    if vbs: print("\n*) Diagonalize matrix S into 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*) Comenzando con el SCF.")
    for iteration in range(50): # Это цикл, в котором матрица плотности должна сходиться

        if vbs: print("**) Calculate Fock matrix.")
        g_scf = G(r, p_scf, b1, b2, b) # Рассчитать матрицу одноэлектронных потенциалов Фока
        f_scf = h_scf + g_scf   # Рассчитать матрицу Фока F = H + G
        
        if vbs: print("**) Calculate F'.")
        f_tra = Xa * f_scf * X # Рассчитать матрицу F' = X_adj * F * X
        
        if vbs: print("**) Diagonalize F' and calculate C'.")
        c_tra = diagon2(m=f_tra) # Рассчитать С' как диагонализацию F'
        
        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(str(iteration + 1) + "-st iteration done.")
        
        if np.linalg.norm(p_temp - p_scf) < 1E-4: # Проверка сходимости Р, то есть, сильно ли она поменялась с 
                                                  # предыдущей итерации
            print("\n-->Finished, convergence achieved!")
            return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp} # Если Р сошлась, прервать цикд
            skip
        else:
            p_scf = p_temp
    print("\n-->Finished, convergence not achieved")
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp} # P не сошлась, выдать что есть

Рассчитать все матрицы для атома водорода (параметры по умолчанию)

mol=SCF(vbs=True)

*) Calculate overlap matrix S.

*) Calculate hamiltonian H.

*) Diagonalize matrix S into X.

*) Create density matrix P.

*) Comenzando con el SCF.
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
1-st iteration done.
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
2-st iteration done.
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
3-st iteration done.
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
4-st iteration done.
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
5-st iteration done.
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
6-st iteration done.
**) Calculate Fock matrix.
**) Calculate F'.
**) Diagonalize F' and calculate C'.
**) Calculate coefficient C matrix.
**) Recalculate P matrix.
7-st iteration done.

-->Finished, convergence achieved!

Рассчитать электронную и общую энергии атома водорода

E_el=ener_elec(mol['P'],mol['H'],mol['F'])
E_tot=ener_tot(elec=E_el)
print(E_el,E_tot)

-1.7974485548087507 -1.1140149845517797

Построить 1D и 2D плот орбиталей

orbital(mol["C"])

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