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

Задание: https://kodomo.fbb.msu.ru/wiki/2016/8/MolSim/task5

Установим зависимости

 conda install -c anaconda sympy
 conda install -c anaconda numpy

Скачаем функцию для расчёта основных матриц hfscf.py

  wget http://kodomo.fbb.msu.ru/~golovin/qm/hfscf.py
  mv hfscf.py term8/pr5
In [2]:
import hfscf
import numpy as np
In [12]:
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]
    
    # Матрица перекрывания
    if vbs: print ("*) Generating overlap matrix S.")
    s_scf = hfscf.S(R, b1, b2, b)
    
    # Гамильтониан
    if vbs: print ("\n*) Generating Hamiltonian H.")
    h_scf = hfscf.H(R, Z, b1, b2, b)
        
    # Диагонализируем матрицу S и находим матрицу X
    if vbs: print ("\n*) Diagonalizing matrix S and finding diagonal matrix X.")
    X = hfscf.diagon(s_scf)
    Xa = X.getH()
    
    # Оцениваем плотность матрицы P
    if vbs: print("\n*) Creating density matrix P.")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)
    
    # Так как F' зависит о решения себя же, минимизацию энергии Хартри-Фока нужно проводить итеративно
    
    # Итеративно оптимизируем матрицу P
    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
        
        # Из C получаем матрицу P
        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 [13]:
r = 1.4632
Z = (1, 1)
b1 = hfscf.GTO["H"]
b2 = hfscf.GTO["H"]
b = 3

b1, b2
Out[13]:
({'a': {2: [1.309756377, 0.233135974],
   3: [3.42525091, 0.62391373, 0.1688554],
   6: [35.52322122,
    6.513143725,
    1.822142904,
    0.625955266,
    0.243076747,
    0.100112428]},
  'c': {2: [0.430128498, 0.678913531],
   3: [0.15432897, 0.53532814, 0.44463454],
   6: [0.00916359628,
    0.04936149294,
    0.1685383049,
    0.3705627997,
    0.4164915298,
    0.1303340841]}},
 {'a': {2: [1.309756377, 0.233135974],
   3: [3.42525091, 0.62391373, 0.1688554],
   6: [35.52322122,
    6.513143725,
    1.822142904,
    0.625955266,
    0.243076747,
    0.100112428]},
  'c': {2: [0.430128498, 0.678913531],
   3: [0.15432897, 0.53532814, 0.44463454],
   6: [0.00916359628,
    0.04936149294,
    0.1685383049,
    0.3705627997,
    0.4164915298,
    0.1303340841]}})
In [14]:
result = SCF(r, Z, b1, b2, b, vbs=True)
result

*) Generating overlap matrix S.

*) Generating Hamiltonian H.

*) Diagonalizing matrix S and finding diagonal matrix X.

*) Creating density matrix P.

*) Starting with the SCF.

**) Generating the Fock matrix: calculating integral of two electrons.
**) Changing the base of F.
**) Diagonalizing F' and generating C'.
**) Constructing coefficient matrix C.
**) Recalculating density matrix P.

Finished the 1. iteration.


**) Generating the Fock matrix: calculating integral of two electrons.
**) Changing the base of F.
**) Diagonalizing F' and generating C'.
**) Constructing coefficient matrix C.
**) Recalculating density matrix P.

Finished the 2. iteration.


**) Generating the Fock matrix: calculating integral of two electrons.
**) Changing the base of F.
**) Diagonalizing F' and generating C'.
**) Constructing coefficient matrix C.
**) Recalculating density matrix P.

Finished the 3. iteration.


**) Generating the Fock matrix: calculating integral of two electrons.
**) Changing the base of F.
**) Diagonalizing F' and generating C'.
**) Constructing coefficient matrix C.
**) Recalculating density matrix P.

Finished the 4. iteration.


**) Generating the Fock matrix: calculating integral of two electrons.
**) Changing the base of F.
**) Diagonalizing F' and generating C'.
**) Constructing coefficient matrix C.
**) Recalculating density matrix P.

Finished the 5. iteration.


**) Generating the Fock matrix: calculating integral of two electrons.
**) Changing the base of F.
**) Diagonalizing F' and generating C'.
**) Constructing coefficient matrix C.
**) Recalculating density matrix P.

Finished the 6. iteration.


**) Generating the Fock matrix: calculating integral of two electrons.
**) Changing the base of F.
**) Diagonalizing F' and generating C'.
**) Constructing coefficient matrix C.
**) Recalculating density matrix P.

Finished the 7. iteration.



-->The self-consistent field IF has converged!
Out[14]:
{'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 [15]:
hfscf.orbital(result['C'], r, b1, b2, b) #Строим график 1D

1D изображение орбиталей.

In [16]:
hfscf.orbital2D(result['C'], result['X'], result['F'], r, b1, b2, b, 0.02) #Строим график 2D

2D изображение орбиталей.

In [17]:
#Рассчитайтываем электронную энергию молекулы и полную энергию молекулы
orbital_energy = hfscf.ener_orbs(result['X'], result['F'])
electronic_energy = hfscf.ener_elec(result['P'], result['H'], result['F'])
total_energy = hfscf.ener_tot(r, Z, 0)

print(f'Orbital energies:\t{orbital_energy[0]:0.5E} and {orbital_energy[1]:0.5E}')
print(f'Electronic energy:\t{electronic_energy:0.5E}')
print(f'Total energy:\t{total_energy:0.5E}')

Orbital energies:	-5.64615E-01 and 6.36867E-01
Electronic energy:	-1.79745E+00
Total energy:	6.83434E-01

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

In [ ]: