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

import hfscf
import numpy as np

def SCF (r=1.4632, Z=[1, 1], b1=hfscf.GTO["H"], b2=hfscf.GTO["H"], b=3 , vbs=False):
    # SCF is an orbital computing function based on the Rutan Hall approach
    
    # r -- electronic degrees of freedom, R -- nuclear degrees of freedom.
    R = [0, r]
    
    # Orbital overlap matrix S.
    if vbs: print ("*) Generating overlap matrix S.")
    s_scf = hfscf.S(R, b1, b2, b)
    
    # Hamiltonian for Fock operator, H matrix.
    if vbs: print ("\n*) Generating Hamiltonian H.")
    h_scf = hfscf.H(R, Z, b1, b2, b)
        
    # Diagonalizing S matrix and finding X matrix.
    # Xa (X_adjoint) hermitian adjoint from the complex Hilbert space.
    if vbs: print ("\n*) Diagonalizing matrix S and finding diagonal matrix X.")
    X = hfscf.diagon(s_scf)
    Xa = X.getH()
    
    # Estimating density matrix P.
    if vbs: print("\n*) Creating density matrix P.")
    p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)
    
    # Since F' depends on it’s own solution (through the orbitals), the process
    # of minimizing Hartree-Fock energy must be done iteratively.
    
    # Begin iterative process.
    if vbs: print("\n*) Starting with the SCF.")
        
    for iteration in range(50):
        
        # Build Fock F matrix.
        # 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   
        
        # Build matrix F'.
        # F' = X_adj * F * X
        if vbs: print("**) Changing the base of F.")
        f_tra = Xa * f_scf * X
        
        # Diagonalize matrix F' and construct matrix C'.
        if vbs: print("**) Diagonalizing F' and generating C'.")
        c_tra = hfscf.diagon2(f_tra)
        
        # Build matrix C.
        # C = X * C'
        if vbs: print("**) Constructing coefficient matrix C.")
        c_scf = X * c_tra
        
        # Build matrix P from matrix C.
        if vbs: print("**) Recalculating density matrix P.")
        p_temp = hfscf.P(c_scf)
        
        print("\nFinished the " + str(iteration + 1) + ". iteration.\n")

        # Review convergence
        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}

r, Z, b = 1.4632, (1, 1), 3
b1 = hfscf.GTO["H"]
b2 = hfscf.GTO["H"]

h = SCF(r, Z, b1, b2, b, vbs=True)

*) 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!

print(h)

{'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 ]])}

orbital_energy = hfscf.ener_orbs(h['X'], h['F'])
electronic_energy = hfscf.ener_elec(h['P'], h['H'], h['F'])
total_energy = hfscf.ener_tot(r, Z, 0)

print(f'Энергии орбиталей:\t{orbital_energy[0]:0.5E} and {orbital_energy[1]:0.5E}')
print(f'Электронная энергия молекулы:\t{electronic_energy:0.5E}')
print(f'Полная энергия молекулы:\t{total_energy:0.5E}')

Энергии орбиталей:	-5.64615E-01 and 6.36867E-01
Электронная энергия молекулы:	-1.79745E+00
Полная энергия молекулы:	6.83434E-01

hfscf.orbital(h['C'], r, b1, b2, b)

hfscf.orbital2D(h['C'], h['X'], h['F'], r, b1, b2, b, 0.02)