18. Hartree-Fock-Roothan (HF)

Warning

Si está utilizando Google Colab o la ejecución en línea, debe de ejecutar al inicio el siguiente código

!pip install pyscf

Tenemos que resolver:

\[ \textbf{F} \textbf{C} = \textbf{S}\textbf{C}\varepsilon \]

• Paso 1. Especificar molécula: Coordenadas de los núcleos \(\{R_A\}\), Carga de los núcleos \(\{Z_A\}\), Número de electrones \((N)\), y funciones base \(\{\phi_i\}\)

• Paso 2. Calcular \(S\), \(H\), \((ij|kl)\).

• Paso 3. Proponer una matriz \(C\).

• Paso 4. Calcular \(P\), \(J\) y \(K\).

• Paso 5. Calcular \(F=H+2J-K\)

• Paso 6. Resolver \(FC=SC\varepsilon\)

• Paso 7. \(E_{elec} = \sum_\mu \sum_\nu P_{\nu \mu} (H_{\mu \nu} + F_{\mu \nu})\)

• Paso 8. ¿\(E_i=E_{i-1}\)?, Sí: acabé. No: volver a paso 4.

• Paso 9. Calcular energía nuclear y sumarla a la energía electrónica.

# Descomentar estas líneas si está en modo online

#!pip install pyscf

import pyscf
import numpy as np
from scipy.linalg import eigh
import sympy as sp
sp.init_printing()

Paso 1. Especificar molécula: Coordenadas de los núcleos \(\{R_A\}\), Carga de los núcleos \(\{Z_A\}\), Número de electrones \((N)\), y funciones base \(\{\phi_i\}\).

H2 = pyscf.gto.Mole(atom = """
    H 0.0000  0.0000 0.0000
    H 0.0000  0.0000 0.7414 
    """,basis = "STO-3G")
H2 = H2.build()

Paso 2. Calcular \(S\), \(H\), \((ij|kl)\).

\[ S_{ij} = \int \psi_i^*(r) \psi_j(r) dr \]

\[ H_{ij} = \int \psi_i^*(r) \hat{H} \psi_j(r) dr \]

\[ H_{ij} = T_{ij} + V_{ij} \]

\[ \psi_i(r) = \sum_\mu a_\mu \phi_\mu(r) \]

\[ (\mu \nu|\sigma \lambda) = \int \int \phi_\mu^*(r_1) \phi_\nu(r_1) \frac{1}{r_{12}} \phi_\sigma^*(r_2) \phi_\lambda(r_2) dr_1 dr_2 \]

S = H2.intor('int1e_ovlp')
print("----------------Matriz S----------------")
sp.pprint(sp.Matrix(S))

T = H2.intor('int1e_kin')
V = H2.intor('int1e_nuc')
H=T+V
print("----------------Matriz H----------------")
sp.pprint(sp.Matrix(H))

I = H2.intor('int2e')
#print("Integrales (ij|kl):")
#print(I)

nbf = S.shape[0]
ndocc = int(H2.nelectron/2)
#print(nbf)
#print(ndocc)

----------------Matriz S----------------
⎡       1.0         0.658957120274098⎤
⎢                                    ⎥
⎣0.658957120274098         1.0       ⎦
----------------Matriz H----------------
⎡-1.12005114184512   -0.957732221404315⎤
⎢                                      ⎥
⎣-0.957732221404315  -1.12005114184512 ⎦

Paso 3. Proponer una matriz C.

C = np.zeros((nbf,nbf))
print(C)

[[0. 0.]
 [0. 0.]]

• Paso 4. Calcular \(P\), \(J\) y \(K\).

\[ P_{\mu \nu} = \sum_a^{N/2} C_{\mu a} C_{\nu a}^* \]

\[ J_{\mu \nu} = \sum_{\lambda \sigma} P_{\lambda \sigma} (\mu \nu | \sigma \lambda) \]

\[ K_{\mu \nu} = \sum_{\lambda \sigma} P_{\lambda \sigma} (\mu \lambda | \sigma \nu) \]

• Paso 5. Calcular \(F=H+2J-K\)

• Paso 6. Resolver \(FC=SC\varepsilon\)

• Paso 7. \(E_{elec} = \sum_\mu \sum_\nu P_{\nu \mu}(H_{\mu \nu} + F_{\mu \nu})\)

• Paso 8. ¿\(E_i=E_{i-1}\)?, Sí: acabé. No: volver a paso 4.

E_old = -1.0
converged=False

while(not converged):
    print("---------------------------------ITERACION---------------------------------")  

    print("------------Matriz C Entrada------------")
    sp.pprint(sp.Matrix(C))    
    
    #Paso 4. Calcular P, J y K
    P = np.zeros((nbf,nbf))
    for i in range(nbf):
        for j in range(nbf):
            for k in range(ndocc):
                    P[i][j] = P[i][j] + C[i][k]*C[j][k]
    print("----------------Matriz P----------------")
    sp.pprint(sp.Matrix(P))
                    
    J = np.zeros((nbf,nbf))
    for i in range(nbf):
        for j in range(nbf):
            for k in range(nbf):
                for l in range(nbf):
                    J[i][j] = J[i][j] + P[k][l]*I[i][j][l][k]
    print("----------------Matriz J----------------")
    sp.pprint(sp.Matrix(J))
                    
    K = np.zeros((nbf,nbf))
    for i in range(nbf):
        for j in range(nbf):
            for k in range(nbf):
                for l in range(nbf):
                    K[i][j] = K[i][j] + P[k][l]*I[i][l][k][j]         
    print("----------------Matriz K----------------")
    sp.pprint(sp.Matrix(K))                    

    #Paso 5. Calcular F = H + 2J - K
    F = H + 2*J - K
    print("----------------Matriz F----------------")
    sp.pprint(sp.Matrix(F))
    
    #Paso 6. Resolver FC=SCE
    E,C = eigh(F, S, eigvals_only=False)
    print("-------------Matriz C Salida-------------")
    sp.pprint(sp.Matrix(C))
    
    #Paso 7. Calcular E=sum_i sum_j P_ji (H_ij+F_ij)
    E_elec = 0.0
    for i in range(nbf):
        for j in range(nbf):
            E_elec = E_elec + P[j][i]*(H[i][j] + F[i][j])

    
    print("Energia Electronica: ",E_elec)

    #Paso 8. ¿$E_i=E_{i-1}$?, Si: acabe. No: volver a paso 4.
    if(abs(E_old - E_elec)< 0.0000001):
        converged = True
    else:
        E_old = E_elec

---------------------------------ITERACION---------------------------------
------------Matriz C Entrada------------
⎡0  0⎤
⎢    ⎥
⎣0  0⎦
----------------Matriz P----------------
⎡0  0⎤
⎢    ⎥
⎣0  0⎦
----------------Matriz J----------------
⎡0  0⎤
⎢    ⎥
⎣0  0⎦
----------------Matriz K----------------
⎡0  0⎤
⎢    ⎥
⎣0  0⎦
----------------Matriz F----------------
⎡-1.12005114184512   -0.957732221404315⎤
⎢                                      ⎥
⎣-0.957732221404315  -1.12005114184512 ⎦
-------------Matriz C Salida-------------
⎡-0.548993777186558  1.2108225729832 ⎤
⎢                                    ⎥
⎣-0.548993777186557  -1.2108225729832⎦
Energia Electronica:  0.0
---------------------------------ITERACION---------------------------------
------------Matriz C Entrada------------
⎡-0.548993777186558  1.2108225729832 ⎤
⎢                                    ⎥
⎣-0.548993777186557  -1.2108225729832⎦
----------------Matriz P----------------
⎡0.301394167389564  0.301394167389564⎤
⎢                                    ⎥
⎣0.301394167389564  0.301394167389563⎦
----------------Matriz J----------------
⎡0.672609505851563  0.446338435641006⎤
⎢                                    ⎥
⎣0.446338435641006  0.672609505851562⎦
----------------Matriz K----------------
⎡0.590387599353547  0.528560342139022⎤
⎢                                    ⎥
⎣0.528560342139022  0.590387599353547⎦
----------------Matriz F----------------
⎡-0.365219729495538  -0.593615692261324⎤
⎢                                      ⎥
⎣-0.593615692261324  -0.365219729495538⎦
-------------Matriz C Salida-------------
⎡-0.548993777186558  -1.2108225729832⎤
⎢                                    ⎥
⎣-0.548993777186557  1.2108225729832 ⎦
Energia Electronica:  -1.830438380772959
---------------------------------ITERACION---------------------------------
------------Matriz C Entrada------------
⎡-0.548993777186558  -1.2108225729832⎤
⎢                                    ⎥
⎣-0.548993777186557  1.2108225729832 ⎦
----------------Matriz P----------------
⎡0.301394167389564  0.301394167389564⎤
⎢                                    ⎥
⎣0.301394167389564  0.301394167389563⎦
----------------Matriz J----------------
⎡0.672609505851563  0.446338435641007⎤
⎢                                    ⎥
⎣0.446338435641007  0.672609505851563⎦
----------------Matriz K----------------
⎡0.590387599353547  0.528560342139022⎤
⎢                                    ⎥
⎣0.528560342139022  0.590387599353547⎦
----------------Matriz F----------------
⎡-0.365219729495538  -0.593615692261324⎤
⎢                                      ⎥
⎣-0.593615692261323  -0.365219729495537⎦
-------------Matriz C Salida-------------
⎡-0.548993777186558  -1.2108225729832⎤
⎢                                    ⎥
⎣-0.548993777186558  1.2108225729832 ⎦
Energia Electronica:  -1.8304383807729592

Paso 9. Calcular energía nuclear y sumarla a la energía electrónica.

\[ E_{\rm nuc} = \sum_{A>B}\frac{Z_AZ_B}{r_{AB}} \]

\[ E_{Tot} = E_{elec} + E_{nuc} \]

E_nuc = H2.energy_nuc()
print("Energia nuclear: ", E_nuc)
E_T = E_elec + E_nuc
print("Energia Total: ", E_T)

Energia nuclear:  0.7137539936876182
Energia Total:  -1.116684387085341

18.1. Método Simple. PySCF

Geometría del agua

Átomo	X (A)	Y (A)	Z (A)
H	0.0000	0.0000	0.0000
H	0.0000	0.0000	0.7414

from pyscf import scf

H2 = pyscf.gto.Mole(atom = """
    H 0.0000  0.0000 0.0000
    H 0.0000  0.0000 0.7414 
    """,basis = "STO-3G")
H2 = H2.build()

rhf = scf.RHF(H2)
rhf.kernel()

converged SCF energy = -1.11668438708534

18.2. Referencias

• D.R. Hartree, The Wave Mechanics of an Atom with a Non-Coulomb Central Field. Part I. Theory and Methods, Math. Proc. Cambridge Philos. Soc. 24, 89 (1928).

• J.C. Slater, Note on Hartree’s method, Phys. Rev. 35, 210 (1930).

• V. Fock, Näherungsmethode zur Lösung des quantenmechanischen Mehrkörperproblems, Zeitschrift für Physik 61, 126 (1930).

• C.C.J. Roothaan, New Developments in Molecular Orbital Theory, Rev. Mod. Phys. 23, 69 (1951).

• Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; Dover Publications: Mineola, N.Y, 1996.