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

In [9]:
import hfscf
from hfscf import GTO
import numpy as np

Вам дан файл с фунциями для расчёта основных матриц hfscf.py
В следующей ячейке описана функция расчёта орбиталей на основе подхода Рутана Холла.

  • Добавьте коментарии в код функции
  • Опишите объекты S,F,X,H,C
  • Рассчитайте электронную энрегию молекулы
  • Рассчитайте полную эннергию молекулы
  • Постройте 1D и 2D плот орбиталей
In [7]:
def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=False):
    R = [0, r] 
    #Генерируем матрицу перекрывания S
    #S показывает, насколько перекрываются рассчитанные атомные орбитали
    if vbs: print("*) Генерирую матрицу перекрывания S.")
    s_scf = hfscf.S(R, b1, b2, b)
    if vbs: print("\n*) Генерирую гамильтониан H.")
    h_scf = hfscf.H(R, Z, b1, b2, b)
    
    #Диагонализируем матрицу S
    if vbs: print("\n*) Диагонализирую матрицу S и получаю диагональную матрицу X.")
    X = hfscf.diagon(m=s_scf)
    #Xa - сопряжённо-траснпонированная к X (транспонирование + замена каждого элемента на комплексно-сопряжённый)
    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.")
    #В ходе каждой итерации пытаемся приблизиться к таким значениям коэффициентов C, чтобы энергия молекулы была
    #минимальной. Всякий раз пересчитываем матрицу Фока (~энергия, хотим минимизировать), матрицу 
    #коэффициентов перед одноатомными орбиталями C и электронную плотность. 
    #Если в какой-то момент электронная плотность перестаёт сильно меняться, считаем, что сошлись.
    for iteracion in range(50):
        # Конструируем матрицу Фока F
        # F = H + G
        if vbs: print("\n**) Конструируем матрицу Фока F")
        g_scf = hfscf.G(r, p_scf, b1, b2, b) #двухэлектронная часть матрицы F
        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 = hfscf.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 = hfscf.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}

S - матрица перекрывания орбиталей. От того, насколько они перекрываются, зависит энергия системы.
F - матрица Фока. Некоторое приближение гамильтониана (в расчётах берём усреднённое взаимодействие электрона со всеми остальными электронами вместо того, чтобы учитывать взаимодействие каждой пары электронов отдельно)
X - диагонализированная матрица Фока
H - гамильтониан (~энергия)
C - матрица коэффициентов перед одноатомными орбиталями (складывая волновые функции для одноатомных орбиталей с этими коэффициентами получаем волновую функцию для молекулы)

In [10]:
my_orbitals = SCF()

Заканчиваю 1 итерацию.


Заканчиваю 2 итерацию.


Заканчиваю 3 итерацию.


Заканчиваю 4 итерацию.


Заканчиваю 5 итерацию.


Заканчиваю 6 итерацию.


Заканчиваю 7 итерацию.



-->Ура, сошлось!
In [11]:
my_orbitals
Out[11]:
{'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 ]])}
In [12]:
#1D графики волновых функций орбиталей атомов в молекуле водорода и квадратов волновых функций
hfscf.orbital(my_orbitals['C'])
In [13]:
#Соответствующие 2D графики
hfscf.orbital2D(my_orbitals['C'], my_orbitals['X'], my_orbitals['F'])
In [18]:
#Можем посчитать энергии
electron_energy = hfscf.ener_elec(my_orbitals['P'], my_orbitals['H'], my_orbitals['F'])
total_energy = hfscf.ener_tot()
In [16]:
#Электронная энергия молекулы
electron_energy
Out[16]:
-1.7974485548087507
In [17]:
#Полная энергия молекулы
total_energy
Out[17]:
0.683433570256971
In [ ]: