How to actually find a Hartree-Fock ground state? I am interested in finding the Hartree-Fock ground state for a system of interacting fermions (with totally local scattering, so a delta-function interaction potential). I have read through some textbooks, and I am comfortable that extremizing (hopefully minimizing) the energy of the variational guess means that we must fulfil the Hartree-Fock equation:
$$ \hat{f}|\chi_a\rangle=\epsilon_a|\chi_a\rangle $$
where $\hat{f}$ is the Fock operator (on the Hilbert space formed by the set of single-particle wavefunctions, not the full Fock space) and the $\epsilon_a$, $|\chi_a\rangle$ are the eigenvalues and eigenfunctions that solve this problem (hopefully the ground state if we have a good initial guess in our iterative method).
It is from here that I am quite confused. Most books/resources I can find handwave that "we find the solution by iterative methods", but I want to actually implement this for a simple system. I have heard of a method in which you take a guess of the wavefunctions, construct the mean-field, use a numerical method to solve for the eigenfunctions of the Hartree-Fock equation assuming a fixed field, and then repeat until the field is self-consistent.


*

*Why would this method necessarily converge?  

*What sort of numerical method would I use to solve the Hartree-Fock equation for a fixed field (equivalent to solving the time-independent Schrödinger equation, I guess)? I've never done anything like this before.
 A: Implementing such a method is not a trivial task. Just FYI, there are multiple programs purpose-built to perform Hartree-Fock calculations. 
In terms of a good introduction to the theory of Hartree-Fock calculations, I found this pdf extremely helpful.
First of all, you have given the expression for the psuedo-exact (psuedo, because we have made the Born-Oppenheimer approximation and assumed a Slater determinant form to get here) 1-electron molecular orbitals $|\chi_a\rangle$, which are not known.
Formally, that's not a problem. We can just use some complete set of states $\{|\phi_\mu\rangle\}_{\nu=1...\infty}$ as a basis, and expand the exact (unknown) states by $|\chi_a\rangle = \sum_{\mu=1}^\infty C_{\mu a}|\phi_k\rangle$ in unknown coefficients $C_{\mu a}$.
Naturally, the requirement of infinitely many basis orbitals is not possible on a computer, so we must instead use a sufficiently large basis that the truncation doesn't matter to within our accuracy guidelines.
Following the analysis in the pdf, there are ultimately 4 matrices that we need.
$F_{\mu \nu}(C_{\alpha k}) = \langle \phi_\mu | \hat{f} | \phi_\nu \rangle$, the Fock matrix (The brackets are included here to represent dependence on $C$!)
$S_{\mu \nu}= \langle \phi_\mu | \phi_\nu \rangle $, the overlap integral of all of your basis functions (This does not change throughout the iteration.)
$C_{\alpha k}$, the coefficients that specify the expansion of the $k$th eigenfunction.
$\mathbf{\varepsilon}$, a diagonal matrix with entries corresponding to the 1-electron eigen-energies of $|\chi_a\rangle$ under Fock operator $\hat{f}$.
The choice of a basis "reduces" the earlier equation to
$$ \mathbf{F}(\mathbf{C}) \mathbf{C} = \mathbf{S} \mathbf{C} \mathbf{\varepsilon} $$
This type of equation is known as a generalised eigenvalue equation, and can be solved by most numerical linear algebra packages (e.g. Eigen, SciPy). Or rather, it would be if $F$ did not depend on $C$.
There are many, many algorithms for dealing with this, but a common approach is the "iteration" many authors refer to. In this setup, we start off with a guess for the form of $C$, get a better value of C and re-substitute and solve again.
$$ \mathbf{F}(\mathbf{C}^{(n)}) \mathbf{C}^{(n+1)} = \mathbf{S} \mathbf{C}^{(n+1)} \mathbf{\varepsilon} $$
where the superscripts on the $C$ matrices refer to the iteration number.

Why would the method necessarily converge

In general, it doesn't. We essentially hope that our starting guess was sufficiently close to the true solution that the algorithm can find the minimum. This is why it is of critical importance that the basis set is carefully chosen to be 'fairly close' to the true wavefunction from the get-go - A natural first choice is atomic hydrogenic wavefunctions, but these are computationally expensive to integrate and better results can be achieved with a larger collection of simpler functions. Such basis sets are tabulated for computational chemists at sites such as basis set exchange.

What sort of numerical method would I use to solve the Hartree-Fock equation?

Ultimately, it comes down to following this procedure:


*

*Choose a set of basis orbitals.

*Compute $\mathbf{S}$.

*Establish a guess for $C^{(0)}$.

*Solve $ \mathbf{F}(\mathbf{C}^{(n)}) \mathbf{C}^{(n+1)} = \mathbf{S} \mathbf{C}^{(n+1)} \mathbf{\varepsilon} $

*Repeat 4. until convergence.


A note on that convergence - generally, one tests to see if the new energies differ from the preceding energies by less than some cutoff, and so establish that the algorithm isn't "moving" very much in state space.
As for why it works at all - Conceivably, for sufficiently restricted conditions the step 4) above may define a contraction mapping $\mathcal{F} : M_{N\times M}(\mathbb{C}) \to M_{N\times M}(\mathbb{C}), \mathbf{C} \in M_{N\times M}(\mathbb{C})$, and so have convergence from the Banach fixed point theorem. For a general basis though, these requirements are not satisfied.
A: The existing answer by catalogue_number is perfectly adequate. I'll just point out that the simple fixed-point iterative scheme given by his points 1.-5. is not really used in practice in this "vanilla" form. Quantum chemistry programs use so-called convergence acceleration algorithms which 1) make convergence faster (fewer iterations are required) and 2) make the algorithms more robust (convergence to the minimum is achieved also for poor starting guesses). For example see:
Alejandro J Garza, Gustavo E. Scuseria, Comparison of self-consistent field convergence acceleration techniques, Journal of Chemical Physics 137, 054110 (2012)
https://www.semanticscholar.org/paper/Comparison-of-self-consistent-field-convergence-Garza-Scuseria/2aebb6f9a09cd83296af1530299a1a44b7a0a2dd
I can also point out that other approaches have been suggested for solving the Hartree-Fock (or Hartree-Fock-Roothaan) equations. For example Malbouisson and Vianna introduced in 1990 an algebraic method [J. Chim. Phys., Vol. 87 (1990), pp. 2017–2025, An algebraic method for solving Hartree-Fock-Roothaan equations. For an application see, e.g., Malbouisson, Sobrinho, de Andrade, Multiple Hartree-Fock Solutions of Systems Constituted with First Line Atoms: BH and FH molecules using the Double Zeta Base,  http://ojs.rpqsenai.org.br/index.php/rpq_n1/article/view/285 ]
A different numerical approach has been used by Shiozaki and Hirata [Grid-based numerical Hartree-Fock solutions of polyatomic molecules, PHYSICAL REVIEW A 76, 040503(R), 2007]
I remember seeing yet another approach which was based on Thouless theorem, but I couldn't find the exact reference.
