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):
    
    # 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 = 1.4632
Z = (1, 1)
b1 = hfscf.GTO["H"]
b2 = hfscf.GTO["H"]
b = 3

result = 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!

result

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

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

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

Красиво: видим связывающую и разрыхляющую орбитали.

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