The foundation of using Hartree-Fock method to get the self-consistence wave function and potential is solving the Hartree-Fock equation, which is kind of differential equation.

The idea of the Hartree-Fock method is

1. Assume the trial wave function $\psi^{(0)}_\mu(i)$
2. put the trial wave function into the Hartree-Fock equation $\displaystyle F(\psi_{\nu}^{(0)}(j))$
3. Solve the equation and obtain a new wave function $\psi^{(1)}_\mu(i)$
4. Go to step 2. until converge.

In this post, we will discuss at least 2 method on the step 3.

One method is solve the Hartree-Fock equation using numerical method, for example, Rungu-Kutta method, or Numerov’s method. Those method can also obtain the energy when impose some conditions, for example, the behaviour of the wave function at long range.

The 2nd method is basis expansion. We assume the wave function is a superposition of some basis,

$\displaystyle \psi_\mu(i) = \sum_{\alpha}^{N_\alpha} a_{\mu\alpha} \phi_\alpha(i)$

Substitute into the Hartree equation (Hartree-Fock equation is similar),

$\displaystyle F \left(\sum_{j=i+1}^{N}\psi_\nu(j) \right) \psi_\mu(i) = H(i) \psi_\mu(i) + \sum_{j=i+1}^{N}\langle \psi_\nu(j) | G(i,j) | \psi_\nu(j) \rangle \psi_\mu(i)$

$\displaystyle \sum_{\alpha} a_{\mu\alpha} F |\phi_{\alpha}(i) \rangle \\ = \sum_{\alpha} a_{\mu\alpha} \left( H(i) |\phi_{\alpha}(i) \rangle + \sum_{k,l} \sum_{j=i+1}^N a_{\nu k}^{*} a_{\nu l}\langle \phi_{\nu k}(j) | G(i,j) | \phi_{\nu l} (j) \rangle |\phi_{\mu \alpha}(i) \rangle \right) \\ = \sum_{\alpha} a_{\mu\alpha} \left( H(i) + G(i) \right) |\phi_{\alpha}(i) \rangle \\ = \sum_{\alpha} a_{\mu \alpha} E |\phi_{\alpha} \rangle$

Acting $\langle \phi_{\mu \beta}(i)|$

$\displaystyle \sum_{\alpha} a_{\mu \alpha} \langle \phi_{\mu \beta} (i) | F | \phi_{\mu \alpha} (i) \rangle = \sum_{\alpha} F_{\beta \alpha} a_{\mu \alpha} = E a_{\mu \beta}$

Thus, we can solve the coefficient $a_{\mu \alpha}$, which is the eigen vector. Then, new trial functions with the new coefficients are used.

In the 2nd method, we can see the integral

$H_{\beta \alpha} = \langle \phi_\beta(i) | H(i)| \phi_\alpha (i) \rangle$

$G_{\beta \alpha}^{k l} = \langle \phi_\beta(i) \phi_k(j) | G(i,j) | \phi_l(j) \phi_\alpha(i) \rangle$

can be pre-calculated. The calculation becomes iterating the coefficient $a_{\mu \alpha}$.

In method 2, the basis does not change, but only the coefficients. The advantage, also the disadvantage of this method is that, all possible integrals are pre calculated, even though some of them will never be used. Once the integrals are calculated, it can be used for solving many similar problems. But the calculation time for the integrals may be long and consume memory.

In previous post, I did not pre-calculate the integral, I simply put the new trial function and integrate again.

This method is very similar to the method 2. However, the basis is the wave function itself. After one calculation, a new basis, which is the wave function, will be used. And the old basis will be replaced. In method 2, only the coefficient is replaced.

This method may somehow faster then method 2, because only nesscary integrals are calculated. In Mathematica, method 2 is generally slower, especially for larger basis.