import hfscf as scf

C:\Users\huvsa\hfscf.py:121: SyntaxWarning: "\P" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\P"? A raw string is also an option.
  plt.plot(dom, psi1, "r--", label="$\Psi_1$")
C:\Users\huvsa\hfscf.py:122: SyntaxWarning: "\P" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\P"? A raw string is also an option.
  plt.plot(dom, psi2, "b--", label="$\Psi_2$")
C:\Users\huvsa\hfscf.py:123: SyntaxWarning: "\P" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\P"? A raw string is also an option.
  plt.plot(dom, psi3, "g-", label="$\Psi_{1}^{2}$")
C:\Users\huvsa\hfscf.py:124: SyntaxWarning: "\P" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\P"? A raw string is also an option.
  plt.plot(dom, psi4, "k-", label="$\Psi_{2}^{2}$")

from sympy import *
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

Задание №1 Рассчёт матриц S,F,X,H,C, а также электронной и полной энергии молекулы водорода в атомных единицах.

def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=False):        #На вход функции мы подаём геометрию молекулы и базисные функции
    R = [0, r]
    if vbs: print("*) Строим матрицу перекрывания.")
    s_scf = scf.S(R, b1, b2, b) 
    if vbs: print("\n*) Строим гамильтониан H.")
    h_scf = scf.H(R, Z, b1, b2, b)
        
    # Diagonalizar matriz S y hallar matriz X
    if vbs: print("\n*) Диагонализируем матрицу S и находим диагональную матрицу X.")
    X = scf.diagon(m=s_scf)
    Xa = X.getH()
    
    # Estimar matriz de densidad P
    if vbs: print("\n*) Строим матрицу плотности P.")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)  # Referencia (7) p. 148
    
    # Comenzar proceso iterativo
    if vbs: print("\n*) Начинаем расчёт SCF.")
    for iteracion in range(50):
        # Construir matriz de Fock F
        # F = H + G
        if vbs: print("\n**) Вычисляем двухэлектронные интегралы и строим матрицу Фока.")
        g_scf = scf.G(r, p_scf)
        f_scf = h_scf + g_scf   # Referencia (7) p. 141 eq. (3.154)
        
        # Construir matriz F'
        # F' = X_adj * F * X
        f_tra = Xa * f_scf * X
        
        # Diagonalizar matriz F y constuir matriz C'
        if vbs: print("**) Диагонализируем матрицу Фока и строим матрицу C.")
        c_tra = scf.diagon2(m=f_tra)
        
        # Construir matriz C
        # C = X * C'
        if vbs: print("**) Строим матрицу коэффициентов C.")
        c_scf = X * c_tra
        
        # Construir matriz P a partir de matriz C
        if vbs: print("**) Перерасчитываем матрицу плотности P.")
        p_temp = scf.P(C=c_scf)
        
        print("\nЗавершаю " + str(iteracion + 1) + ". итерацию.\n")
        
        # Revisar convergencia
        if np.linalg.norm(p_temp - p_scf) < 1E-4: # Referencia (7) p. 148
            print("\n\n-->Самосогласованное поле сошлось.")
            eE=scf.ener_elec(p_temp, h_scf, f_scf)
            Etot=scf.ener_tot(r = 1.4632, Z=[1,1], elec=eE)
            return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp, "Electronic_energy":str(float(eE))+' a.u', "Total_energy":str(float(Etot))+' a.u'}
        else:
            p_scf = p_temp
    print("\n\n-->Самосогласованное поле НЕ сошлось.")
    return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}

SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=False)

Завершаю 1. итерацию.


Завершаю 2. итерацию.


Завершаю 3. итерацию.


Завершаю 4. итерацию.


Завершаю 5. итерацию.


Завершаю 6. итерацию.


Завершаю 7. итерацию.



-->Самосогласованное поле сошлось.

{'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 ]]),
 'Electronic_energy': '-1.7974485548087507 a.u',
 'Total_energy': '-1.1140149845517797 a.u'}

Задание № 2, теперь построим 1D и 2D плоты получённых молекулярных орбиталей.

scf.orbital(np.matrix([[ 0.55258114, -1.17442421],
         [ 0.55258013,  1.17442468]]), r=1.4632, b1 = GTO["H"], b2 = GTO["H"], b = 3)

scf.orbital2D (np.matrix([[ 0.55258114, -1.17442421],
         [ 0.55258013,  1.17442468]]), np.matrix([[ 0.55258063,  1.17442445],
         [ 0.55258063, -1.17442445]]), np.matrix([[-0.34684107, -0.57771139],
         [-0.57771139, -0.34684028]]), r = 1.4632, b1 = GTO["H"], b2 = GTO["H"], b = 3, delta=0.02)