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

Ниже задана функция расчёта орбиталей на основе подхода Рутана Холла

S - матрица перекрывания; p. 412 eq. (A.9)

F - матрица Фока;

X - диагонализированная матрица S;

H - гамильтониан; p. 141 eq. (3.153)

C - матрица коэффициентов.

In [22]:
from hfscf import *
import numpy as np
from IPython.display import display,Image
In [14]:
def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=False):
    R = [0, r]
    #r -- электронные степени свободы, R -- ядерные степени свободы
    if vbs: print("*) Generating overlap matrix S") #Генерируем матрицу перекрывания S
    s_scf = S(R, b1, b2, b)
    if vbs: print("\n*) Generating Hamiltonian H.") #Генерируем гамильтониан H
    h_scf = H(R, Z, b1, b2, b)
        
    # Диагонализировать матрицу S и найти матрицу X
    if vbs: print("\n*) Diagonalizing matrix S and finding diagonal matrix X.")#Диагонализируем матрицу S и ищем матрицу X
    X = diagon(m=s_scf)
    Xa = X.getH()
    
    # Cоздать матрицу электронной плотности P
    if vbs: print("\n*) Creating density matrix P.")#Cоздаем матрицу электронной плотности P
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)  # Referencia (7) p. 148
    
    # Запуск итеративных процессов
    if vbs: print("\n*) Starting with the SCF.") #Начинаем считать SCF.
    for iteracion in range(50):
        # Строим матрицу Фока F
        # F = H + G
        if vbs: print("\n**) Generating the Fock matrix: calculating \
integral of two electrons.") #Генерируем матрицу Фока: рассчитываем интегралы двух электронов.
        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("**) Changing the base of F.") # Меняем базис F
        f_tra = Xa * f_scf * X
        
        # диагонализировать матрицу F и построить C'
        if vbs: print("**) Diagonalizing F' and generating C'.") #диагонализируем F' и генерируем C'
        c_tra = diagon2(m=f_tra)
        
        # Строим матрицу C
        # C = X * C'
        if vbs: print("**) Constructing coefficient matrix C.") #Строим матрицу коэффициентов C
        c_scf = X * c_tra
        
        # Строим матрицу P на основе матрицы C
        if vbs: print("**) Recalculating density matrix P.") #Перерасчет матрицы плотности P
        p_temp = P(C=c_scf)
        
        # Calculating orbital, electronic and total energy.

        orbital_energy = hfscf.ener_orbs(X, f_scf)
        electronic_energy = hfscf.ener_elec(p_temp, h_scf, f_scf)
        total_energy = hfscf.ener_tot(r, Z, 0)
        
        # Plot orbitals.
        #hfscf.orbital(c_tra, r, b1, b2, b)
        #hfscf.orbital2D(c_scf, X, f_scf, r, b1, b2, b, 0.02)
        #Очень хотелось каждый раз показывать картиночки с орбиталями, но они выводились в отдельное окно и скрипт дальше не работал
        
        
        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,"OE":orbital_energy,"EE":electronic_energy, "TE":total_energy}

        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,"OE":orbital_energy,"EE":electronic_energy, "TE":total_energy}
In [15]:
result = SCF()

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


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


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


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


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


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


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



-->Самосогласованное поле сходится

Считаем электронную и полную энергию молекулы

In [16]:
print("Электронная энергия молекулы " + str(result.get("EE")))
print("Общая энергия молекулы " + str(result.get("TE")))

Электронная энергия молекулы -1.79744855481
Общая энергия молекулы 0.683433570257
In [28]:
#1D плот орбиталей, хз как нормально вставить картинку в notebook -- она выводится в отдельное окно :|
#orbital(result["C"])
Image("figure_1.png")
Out[28]:
In [27]:
#2D плот орбиталей
#orbital2D(result["C"],result["X"],result["F"])
Image("figure_2Dorb.png")
Out[27]: