Практикум 5

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

--2021-05-27 22:10:48--  http://kodomo.fbb.msu.ru/~golovin/qm/hfscf.py
Распознаётся kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)… 93.180.63.127
Подключение к kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:80... соединение установлено.
HTTP-запрос отправлен. Ожидание ответа… 302 Found
Адрес: https://kodomo.fbb.msu.ru/~golovin/qm/hfscf.py [переход]
--2021-05-27 22:10:49--  https://kodomo.fbb.msu.ru/~golovin/qm/hfscf.py
Подключение к kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:443... соединение установлено.
HTTP-запрос отправлен. Ожидание ответа… 200 OK
Длина: 13490 (13K) [text/x-python]
Сохранение в: «hfscf.py»

hfscf.py            100%[===================>]  13.17K  --.-KB/s    за 0.005s  

2021-05-27 22:10:49 (2.45 MB/s) - «hfscf.py» сохранён [13490/13490]

import numpy as np
import seaborn as sns
import matplotlib.pyplot as pyplot
import hfscf

Рассмотрим и откомментируем предложенную функцию

def SCF (r = 1.4632,
         Z=[1,1],
         b1 = hfscf.GTO["H"],
         b2 = hfscf.GTO["H"],
         b = 3,
         vbs=False):
    R = [0, r]
    #так как моделинуем молекулу водорода, то задаем радиус
    #задаем заряд молекулы
    if vbs: print("*) Generando matriz de traslape S.")
    #создаем матрицу перекрывания орбиталей
    s_scf = hfscf.S(R,
              b1,
              b2,
              b)
    if vbs: print("\n*) Generando hamiltoniano H.")
    #Рассчитываем одноэлектронный гамильтониан
    h_scf = hfscf.H(R,
                    Z,
                    b1,
                    b2,
                    b)
    # Diagonalizar matriz S y hallar matriz X
    #Диагонализируем матрицу S  - получаем матрицу трансформации Х
    #Получаем сопряженно транспонированную к Х Ха
    if vbs: print("\n*) Diagonalizando matriz S y hallando matriz diagonal X.")
    X = hfscf.diagon(m=s_scf)
    Xa = X.getH()
    
    # Estimar matriz de densidad P
    #Делаем предварительную оценку матрицы плотности
    if vbs: print("\n*) Creando matriz de densidad P.")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)  # Referencia (7) p. 148
    
    # Comenzar proceso iterativo
    # Начинаем итеративный процесс SCF, (1) рассчитываем F (фокиан) 
    #(2) по матрице плотности, строим F′=S**(−1/2)FS**(−1/2),
    #(3) Решаем |F’- EI|=0 для поиска собственных значений E и С’диаганолизации F’
    #(4)Рассчитываем орбитальные коэффициентыC=S**(−1/2)C′
    #(5)Считаем новую матрицу плотности из матрицы С
    #останавливаемся когда плотность прекращает сильно меняться
    if vbs: print("\n*) Comenzando con el SCF.")
    for iteracion in range(50):
        # Construir matriz de Fock F
        # F = H + G
        #Считаем фокиан - F (шаг 1)
        #фокиан = гамильтониан + двуэлектронный вклад(отталкивание)
        if vbs: print("\n**) Generando la matriz de Fock: calculando \
integrales de dos electrones.")
        g_scf = hfscf.G(r, p_scf, b1, b2, b)
        f_scf = h_scf + g_scf   # Referencia (7) p. 141 eq. (3.154)
        
        #Строим матрицу F' (шаг 2)
        # Construir matriz F'
        # F' = X_adj * F * X
        if vbs: print("**) Cambiando la base de F.")
        f_tra = Xa * f_scf * X
        
        #Диагонализируя матрицу F получаем матрицу С (шаг 3)
        # Diagonalizar matriz F y constuir matriz C'
        if vbs: print("**) Diagonalizando F' y generando C'.")
        c_tra = hfscf.diagon2(m=f_tra)
        
        #Строим матрицу C (шаг 4)
        # Матрица С - матрица коэффицентов перед базисными функциями
        # Construir matriz C
        # C = X * C'
        if vbs: print("**) Construyendo matriz de coeficientes C.")
        c_scf = X * c_tra
        
        #Строим матрицу плотности по матрице коэффицентов (шаг 5)
        # Construir matriz P a partir de matriz C
        if vbs: print("**) Recalculando matriz de densidad P.")
        p_temp = hfscf.P(C=c_scf)
        
        print("\nConcluida la " + str(iteracion + 1) + ". iteracion.\n")
        
        #Проверяем, сошлась ли наша новая матрица плотности с предыдущей
        # Revisar convergencia
        if np.linalg.norm(p_temp - p_scf) < 1E-4: # Referencia (7) p. 148
            print("\n\n-->El campo autoconsistente SI ha convergido!")
            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-->El campo autoconsistente NO ha convergido!\nRevisar supuestos.")
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}

Кто есть кто?

• S - матрица перекрывания, содержит интегралы перекрывания базисных функций.
• H - одноэлектронный гамильтониан системы. Содержит одноэлектронные интегралы.
• X - матрица перехода
• F - фокиан системы. Равен сумме матриц H - гамильтониан и G - двуэлектронное оталкивание
• C - матрица коэффицентов с которыми комбинируются базисные функции для полуения молекулярных орбиталей

Посчитаем орбитали

scf_mtxs = SCF(vbs=True)

*) Generando matriz de traslape S.

*) Generando hamiltoniano H.

*) Diagonalizando matriz S y hallando matriz diagonal X.

*) Creando matriz de densidad P.

*) Comenzando con el SCF.

**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 1. iteracion.


**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 2. iteracion.


**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 3. iteracion.


**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 4. iteracion.


**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 5. iteracion.


**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 6. iteracion.


**) Generando la matriz de Fock: calculando integrales de dos electrones.
**) Cambiando la base de F.
**) Diagonalizando F' y generando C'.
**) Construyendo matriz de coeficientes C.
**) Recalculando matriz de densidad P.

Concluida la 7. iteracion.



-->El campo autoconsistente SI ha convergido!

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

orbital_en = hfscf.ener_orbs(scf_mtxs['X'], scf_mtxs['F'])
elec_en = hfscf.ener_elec(scf_mtxs['P'], scf_mtxs['H'], scf_mtxs['F'])
total_en = hfscf.ener_tot(elec=elec_en)

print(f"Энергия связывающей орбитали:\t{orbital_en[0]}")
print(f"Энергия разрыхляющей орбитали:\t{orbital_en[1]}")
print(f"Электронная энергия:\t{elec_en}")
print(f"Полная энергия:\t\t{total_en}")

Энергия связывающей орбитали:	-0.5646153594583119
Энергия разрыхляющей орбитали:	0.6368673980935621
Электронная энергия:	-1.7974485548087507
Полная энергия:		-1.1140149845517797

Построим плот орбиталей

hfscf.orbital(scf_mtxs['C'])

hfscf.orbital2D(scf_mtxs['C'], scf_mtxs['X'], scf_mtxs['F'])

Две орбитали: связывающая (1) и разрыхляющая (2)