5. Spanish for beginners.¶
import hfscf
import numpy as np
5.1 Some theory.¶
We know how to solve the electronic problem for the simplest atom, hydrogen, which has only one electron. We can pretend that the electrons do not interact with each other in the two electron system. So,
$$\psi(r_1, r_2) = \psi(r_1)\psi(r_2).$$ In general form: $$\psi_{HP}(r_1, r_2 \dots r_N) = \phi(r_1)\phi(r_2) \dots \phi(r_N).$$
But we know, that fermions have not only coordinates in space, but intrinsic spin coordinate ($\alpha$ or $\beta$). That is why we can substitute spatial orbital, $\phi(r_i)$ by $\chi(x_i)$, where $\chi(x) = \phi(r)\alpha.$ However this function fails to satisfy the antisymmetry principle, which states that a wavefunction describing fermions should be antisymmetric with respect to the interchange of any set of space-spin coordinates. Finally we can describe wavefunction for two electrones like this:
$$\psi(x_1, x_2) = \frac{1}{\sqrt{2}}[\chi(x_1)\chi(x_2) - \chi(x_2)\chi(x_1)]$$ In general form: $$\psi = \frac{1}{\sqrt{N}}\begin{vmatrix} \chi_1(x_1) & \chi_2(x_1) & \dots & \chi_N(x_1) \\ \chi_1(x_2) & \chi_2(x_2) & \dots & \chi_N(x_2) \\ \vdots & \vdots & \ddots & \vdots \\ \chi_1(x_N) & \chi_2(x_N) & \dots & \chi_N(x_N) \\ \end{vmatrix}$$
A determinant of spin orbitals is called a Slater determinant. An interesting consequence of this functional form is that the electrons are all indistinguishable. Each electron is associated with every orbital. That simply means that each electron moves independently of all the others except that it feels the Coulomb repulsion due to the average positions of all electrons. It also experiences an “exchange” interaction due to antisymmetrization. Now we can express electronic Hamiltonian: $$H_{el} = \sum_i{h(i)} + \sum_{i < j}{v(i, j)} + V_{NN},$$ where $h(i) = -\frac{1}{2}\nabla_i^2 - \sum_A{\frac{Z_A}{r_{iA}}}$ -- one-electron operator, $v(i,j) = \frac{1}{r_{ij}}$ -- two-electron operator, and $V_{NN}$ is just a constant for the fixed nuclear coordinates ${R}$. We can express average electronic energy represented by Hamiltonian in the particular state $\psi(r_i)$:
$$E_{el} = (\psi|H_{el}|\psi), \text{it is a bra-ket notation}.$$
Due to variational theorem obtained energy is always an upper bound to the true energy. So, we need to minimize obtained energy with respect to changes in the orbitals $\chi_i = \chi_i + \delta\chi_i$ . We also assume that the orbitals are orthonormal, and we want to ensure that our variational procedure leaves them orthonormal. Specific functional $L$ is employed. Setting the first variation $\delta L = 0$, we arrive at Hartree-Fock equations. They can be solved numerically or in the space spanned by a set of basis functions (Hartree-Fock-Roothan equations). The Hartree-Fock-Roothaan equations can be written in matrix form as:
$$FC=SC_{\epsilon},$$
where $F$ is a Fock matrix with changed base due to Hartree-Fock-Roothaan transformation, $C$ is a coefficient matrix that is needed for accounting electron-electron interactions, $S$ is an overlap matrix, and $C_{\epsilon}$ is a diagonal matrix of orbital energies.
# STO -- Slater-type orbitals equal to atomic orbitals
# GTO -- Gaussian-type orbitals -- saving large amount of computation time -- equal to STO
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]
# S matrix -- orbital overlap matrix -- square matrix that describes inter-relationships
# of an atomic basis set (set of functions that is used to represent
# the electronic wave function in the Hartree–Fock method.
if vbs: print ("*) Generating overlap matrix S.")
s_scf = hfscf.S(R, b1, b2, b)
# H matrix -- Hamiltonian for Fock operator.
# H = Ke + Pe, where Ke -- kinetic energy of electron, Pe -- potential energy of electron.
if vbs: print ("\n*) Generating Hamiltonian H.")
h_scf = hfscf.H(R, Z, b1, b2, b)
# Diagonalize matrix S and find matrix X.
# Diagonal matrix S = X == orthogonal vectors in basis set.
# Xa -- X_adjoint -- hermitian adjoint from the complex Hilbert space.
# Adjoints are useful for bra-ket notation in quantum mechanics and
# are used as bras.
if vbs: print ("\n*) Diagonalizing matrix S and finding diagonal matrix X.")
X = hfscf.diagon(s_scf)
Xa = X.getH()
# Estimate 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
# H -- Hamiltonian
# G = p*(J-0.5*K), where p -- p_psf, J -- Coulomb operator, K -- exchange operator
# J gives the Coulomb interaction of an electron in partiacular spin orbital
# with the average charge distribution of the other electrons.
# K defines the electron "exchange" energy. It arises from
# the antisymmetry requirement of the wavefunction.
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
# Since this step we introduce a basis set that transforms
# the Hartree-Fock equations into the Hartree-Fock-Roothaan equations.
# The Roothaan equations are a representation of the Hartree–Fock equation
# in a non orthonormal basis set which can be of GTO or STO.
# Build matrix F'.
# F' = X_adj * F * X
# As we transform HF equations into the HF-Roothaan equations
# we need to change the base of matrix F == F'.
if vbs: print("**) Changing the base of F.")
f_tra = Xa * f_scf * X
# Diagonalize matrix F' and construct matrix C'.
# Diagonilizing matrix F' is a necessary step because of (1) going to an
# orthogonal basis and therefore (2) eliminating overlap matrix S.
# C' = diagonal F'
if vbs: print("**) Diagonalizing F' and generating C'.")
c_tra = hfscf.diagon2(f_tra)
# Build matrix C.
# C = X * C'
# !!do not understand why do we need to multiply diagonal overlap matrix S
# and diagonal F' matrix!!
# i think that we should recalculate density matrix P using diagonal F' matrix
# (i thought that diagonilizing matrix F' is our intention not to use overlap matrix S).
# Product of multiplying diagonal overlap matrix S and diagonal matrix F' is a
# coefficient matrix C.
if vbs: print("**) Constructing coefficient matrix C.")
c_scf = X * c_tra
# Build matrix P from matrix C.
# Recalculate matrix P in HF-Roothaan equations system.
if vbs: print("**) Recalculating density matrix P.")
p_temp = hfscf.P(c_scf)
# Calculating orbital, electronic and total energy.
orbital_energy = hfscf.ener_orbs(X, f_scf)
electronic_energy = hfscf.ener_elec(p_temp, h_scf, f_scf)
total_energy = hfscf.ener_tot(r, Z, 0)
# Plot orbitals.
hfscf.orbital(c_tra, r, b1, b2, b)
hfscf.orbital2D(c_scf, X, f_scf, r, b1, b2, b, 0.02)
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,"OE":orbital_energy,
"EE":electronic_energy, "TE":total_energy}
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,"OE":orbital_energy,
"EE":electronic_energy, "TE":total_energy}
5.2 Displaying the result.¶
result = SCF()
print("Orbital energy equals to " + str(result.get("OE")[0]) + " and " + str(result.get("OE")[1]))
print("Electronic energy equals to " + str(result.get("EE")))
print("Total energy equals to " + str(result.get("TE")))