#set text(14pt)

= SELF CONSISTENT FIELD CALCULATION 
\

== Step 1 (Initial Requirements)
To solve the many body problem (to find the ground state energy of our given molecule), first we need three parameters.

1. *Distance (r) :* we need the distance for each particle to each particle in order to calculate their potential. Mostly it require $r_(i j)$ for distance between $i^(t h)$ electron to $j^(t h)$ electron and $r_(i A)$ for the distance between $i^(t h)$ electron to $A^(t h)$ nucleus.

2. *Atomic Number (Z) :* The potential also depend on the $Z_A$ for an $A^(t h)$ atom with interact with $i^(t h)$ electron.

3. *N & K :* $N$ is the total number of electron and $K$ is the total number of atom for accounting nucleus.

== Step 2 (Wavefunction $Phi$)
- we take the initial wavefunction $Phi$, which is used to calculate the integral values. The wavefunction can be anything like gaussian wavefunction ($A e^(-alpha r^2)$), slatter type orbital wavefunction ($A e^(-alpha r)$) or anything but we need to find parameter $alpha$ and coefficient $A$. 
- The wavefunction can also take single or linear combination of something. For example: the three set of gaussian wavefuntion can take to approximate slatter type orbitals. for that three coefficient and parameter are need to find, that is known as contracted gaussian wave function.
- The coeffiecient and alpha can be taken from internet or the coefficient can be solved by normalising the wave function and the alpha parameter is solved by using hellmann feymann method used in variational principle.

== Step 3 (Calculate Integrals)

we need to calculate four integrals for He+H molecule.

1.  Overlap integral *S*

  $ S_("int") = (pi/(A+B))^(1.5) exp((-A*B*"Rab2")/(A+B)) $

2. Kinetic Energy integral *T*

  $ T_("int") = (A*B)/(A+B) (3- (2*A*B*"Rab2")/(A+B)) (pi/(A+B))^(1.5) exp((-A*B*"Rab2")/(A+B)) $

3. Potential Energy integral *V*

  $ V_("int") = - Z_c * 2 pi/(A+B) * op("F0") ((A+B)* "Rcp2") * exp((-A*B*"Rab2")/(A+B)) $

4. Electron-Electron potential energy *Two_E*

  $ op(T w o E)_("int") = 2.0*(pi**2.5)/((A+B)*(C+D)*sqrt(A+B+C+D))* \ op("F0")(((A+B)*(C+D)*"Rpq2")/(A+B+C+D))* \ exp((-A*B*"Rab2")/(A+B)-(C*D*"Rcd2")/(C+D)) $

== Step 4 (Find $H$, $S$, $X$, $TT$)
For the He+H molecule.
1. *Core Hamiltonian*
  $ H = mat(("T11"+"V11A"+"V11B"),("T12"+"V12A"+"V12B");("T12"+"V12A"+"V12B"),("T22"+"V22A"+"V22B")) $

2. *Overlap Matrix*
  $ S = mat(1,"S12";"S12",1) $

3. *$S^(-1/2)$*
  $ X = mat(1/sqrt(2 (1+"S12")),1/sqrt(2 (1-"S12"));1/sqrt(2 (1+"S12")),-1/sqrt(2 (1-"S12"))) $

4. *Coulumb and Exchange term*
  $ TT = mat("V1111","V2111","V2111","V2211";"V2111","V2121","V2121","V2221";"V2111","V2121","V2121","V2221";"V2211","V2221","V2221","V2222") $
== Step 5 (Guess Initial Density Matrix)
The initial density matrix of electron considered in this code $rho = mat(0,0;0,0)$

== Step 6 ( Calculate Fock Matrix to calculate energy) 
\

- $G = mat(0,0;0,0)$

- Update $G$
  $ G = sum_(i,j)^n sum_(k,l)^n G(i,j) + rho(k,l)* TT(i,j,k,l) - 1/2 TT(i,l,k,j)  $

- $F = H+G$
- Energy = $sum_i 1/2 rho(H+F)$

== Step 7 (Diagonalise Fock matrix to find new density matrix)
\ 
- $F' = S^(1/2) dot F dot S^(-1/2)$

- $|F' - E I|=0$
  - From that we can get eigen value $E$ and eigen vector $C'$

- $C = X dot C'$

- New density matrix 
  $ rho' = sum_(i,j)^n sum_k^1 rho'(i,j) + 2 * C(i,k)* C(j,k) $

== Step 8 (Check energy is converging)
\
$delta = sqrt(sum_i (rho' - rho)^2)/4$

- If delta is less than our appromiation limit ($i.e,1e-11$) we can stop the iteration then,
  - Energy = Energy + nuclear potential terms
  - Density matrix = $rho'$
  - Mullikan population = $rho' dot S$
  - New coeffiecient = C

- Else go back to step 6 with updated density matrix $rho'$