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

hfscf.py - файл с фунциями для расчёта основных матриц

In [2]:
import hfscf
import numpy as np

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

In [13]:
def SCF (r=1.4632, Z=[1, 1], b1=hfscf.GTO["H"], b2=hfscf.GTO["H"], b=3 , vbs=False):
    
    # r - электронные степени свободы, R - ядерные
    R = [0, r]
    
    # Создание матрицы перекрываний S
    if vbs: print ("*) Generating overlap matrix S.")
    s_scf = hfscf.S(R, b1, b2, b)
    
    # Создаём матрицу H (Гамильтониан)
    if vbs: print ("\n*) Generating Hamiltonian H.")
    h_scf = hfscf.H(R, Z, b1, b2, b)
        

    # Xa - сопряженная транспонированная матрица
    if vbs: print ("\n*) Diagonalizing matrix S and finding diagonal matrix X.")
    # Диагонализируем матрицу S
    X = hfscf.diagon(s_scf)
    # Находим матрицу X
    Xa = X.getH()
    
    # Рассчитываем матрицу расчетной плотности P
    if vbs: print("\n*) Creating density matrix P.")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)
    
    # Начинаем итерации
    if vbs: print("\n*) Starting with the SCF.")
        
    for iteration in range(50):
        
        # Конструируем матрицу Фока
        # F = H + G
        if vbs: print("\n**) Generating the Fock matrix: calculating integral of two electrons.")
        g_scf = hfscf.G(r, p_scf, b1, b2, b)
        f_scf = h_scf + g_scf   
        
        # Генерируем матрицу F'
        # F' = X_adj * F * X
        if vbs: print("**) Changing the base of F.")
        f_tra = Xa * f_scf * X
        
        # Диагонализируем F' и генерируем C'
        if vbs: print("**) Diagonalizing F' and generating C'.")
        c_tra = hfscf.diagon2(f_tra)
        
        # Конструируем матрицу C
        # C = X * C'
        if vbs: print("**) Constructing coefficient matrix C.")
        c_scf = X * c_tra
        
        # Конструируем матрицу Р на основе матрицы С
        if vbs: print("**) Recalculating density matrix P.")
        p_temp = hfscf.P(c_scf)
        
        print("\nFinished the " + str(iteration + 1) + ". iteration.\n")

        # Проверяем сходимость
        if np.linalg.norm(p_temp - p_scf) < 1E-4:
            print("\n\n-->The self-consistent field IF has converged!")
            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-->The self-consistent field has NOT converged!\nReview assumptions.")
    return {"S" : s_scf, "H" : h_scf , "X" :  X , "F" : f_scf ,"C" : c_scf, "P" : p_temp}
In [14]:
result = SCF()

Finished the 1. iteration.


Finished the 2. iteration.


Finished the 3. iteration.


Finished the 4. iteration.


Finished the 5. iteration.


Finished the 6. iteration.


Finished the 7. iteration.



-->The self-consistent field IF has converged!
In [15]:
result
Out[15]:
{'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 [16]:
print('Энергии орбиталей: ' + 
      str(hfscf.ener_orbs(result['X'], result['F'])[0]) + 
      ' и ' + 
      str(hfscf.ener_orbs(result['X'], result['F'])[1]))
print('Электронная энергия молекулы: ' + str(hfscf.ener_elec(result['P'], result['H'], result['F'])))
print('Полная энергия молекулы: ' + str(hfscf.ener_tot()))

Энергии орбиталей: -0.5646153594583119 и 0.6368673980935621
Электронная энергия молекулы: -1.7974485548087507
Полная энергия молекулы: 0.683433570256971

1D плот орбиталей

In [17]:
hfscf.orbital(result['C'])

2D плот орбиталей

In [18]:
hfscf.orbital2D(result['C'], result['X'], result['F'])