link code to each algo step; minor adjustemnts to theory and algo

docs/source/documentation/algorithm.md

@@ -12,30 +12,30 @@ kernelspec:
 ---
 # Algorithm overview
 
-## Self-consistency loop
+## Self-consistent mean-field loop
 
 In order to calculate the mean-field interaction in {eq}`mf_infinite`, we require the ground-state density matrix $\rho_{mn}(R)$.
-However, the density matrix in {eq}`density` is a functional of the mean-field interaction $\hat{V}_{\text{MF}}$ itself.
+However, {eq}`density` is a functional of the mean-field interaction $\hat{V}_{\text{MF}}$ itself.
 Therefore, we need to solve for both self-consistently.
 
-We define a single iteration of a self-consistency loop:
+A single iteration of this self-consistency loop is a function that computes a new mean-field correction from a given one:
 
 $$
-\text{SCF}(\hat{V}_{\text{init, MF}}) \to \hat{V}_{\text{new, MF}},
+\text{MF}(\hat{V}_{\text{init, MF}}) \to \hat{V}_{\text{new, MF}},
 $$
 
-such that it performs the following operations given an initial mean-field interaction $\hat{V}_{\text{init, MF}}$:
-
+which is defined in {autolink}`~pymf.model.Model.mfield` method.
+It performs the following steps:
 1. Calculate the total Hamiltonian $\hat{H}(R) = \hat{H_0}(R) + \hat{V}_{\text{init, MF}}(R)$ in real-space.
-2. Fourier transform the total Hamiltonian to the momentum space $\hat{H}(R) \to \hat{H}(k)$.
-3. Calculate the ground-state density matrix $\rho_{mn}(k)$ in momentum space.
-    1. Diagonalize the Hamiltonian $\hat{H}(k)$ to obtain the eigenvalues and eigenvectors.
-    2. Calculate the fermi level $\mu$ given the desired filling of the unit cell.
-    3. Calculate the density matrix $\rho_{mn}(k)$ using the eigenvectors and the fermi level $\mu$ (currently we do not consider thermal effects so $\beta \to \infty$).
-4. Inverse Fourier transform the density matrix to real-space $\rho_{mn}(k) \to \rho_{mn}(R)$.
-5. Calculate the new mean-field interaction $\hat{V}_{\text{new, MF}}(R)$ via {eq}`mf_infinite`.
+2. ({autolink}`~pymf.mf.construct_density_matrix`) Compute the ground-state density matrix $\rho_{mn}(R)$:
+    1. ({autolink}`~pymf.tb.transforms.tb_to_kgrid`) Fourier transform the total Hamiltonian to momentum space $\hat{H}(R) \to \hat{H}(k)$.
+    2. ({autolink}`numpy.linalg.eigh`) Diagonalize the Hamiltonian $\hat{H}(R)$ to obtain the eigenvalues and eigenvectors.
+    3. ({autolink}`~pymf.mf.fermi_on_kgrid`) Calculate the fermi level given the desired filling of the unit cell.
+    4.  ({autolink}`~pymf.mf.construct_density_matrix_kgrid`) Calculate the density matrix $\rho_{mn}(k)$ using the eigenvectors and the fermi level.
+    5. ({autolink}`~pymf.tb.transforms.kgrid_to_tb`) Inverse Fourier transform the density matrix to real-space $\rho_{mn}(k) \to \rho_{mn}(R)$.
+3. ({autolink}`~pymf.mf.meanfield`) Calculate the new mean-field correction $\hat{V}_{\text{new, MF}}(R)$ using {eq}`mf_infinite`.
 
-## Self-consistency criterion
+## Self-consistency criteria
 
 To define the self-consistency condition, we first introduce an invertible function $f$ that uniquely maps $\hat{V}_{\text{MF}}$ to a real-valued vector which minimally parameterizes it:
 
@@ -43,13 +43,15 @@ $$
 f : \hat{V}_{\text{MF}} \to f(\hat{V}_{\text{MF}}) \in \mathbb{R}^N.
 $$
 
+In the code, $f$ corresponds to the {autolink}`~pymf.params.rparams.tb_to_rparams` function (inverse is {autolink}`~pymf.params.rparams.rparams_to_tb`).
 Currently, $f$ parameterizes the mean-field interaction by taking only the upper triangular elements of the matrix $V_{\text{MF}, nm}(R)$ (the lower triangular part is redundant due to the Hermiticity of the Hamiltonian) and splitting it into a real and imaginary parts to form a real-valued vector.
 
-With this function, we define the self-consistency criterion as a fixed-point problem:
+With this, we define the self-consistency criterion as a fixed-point problem:
 
 $$
-f(\text{SCF}(\hat{V}_{\text{MF}})) = f(\hat{V}_{\text{MF}}).
+f(\text{MF}(\hat{V}_{\text{MF}})) = f(\hat{V}_{\text{MF}}).
 $$
 
-To solve this fixed-point problem, we utilize a root-finding function `scipy.optimize.anderson` which uses the Anderson mixing method to find the fixed-point solution.
-However, our implementation also allows to use other custom fixed-point solvers by either providing it to `solver` or by re-defining the `solver` function.
+Instead of solving the fixed point problem, we rewrite is as the difference of the two successive self-consistent mean-field iterations in {autolink}`~pymf.solvers.cost`.
+That re-defines the problem into a root-finding problem which is more consistent with available numerical solvers such as {autolink}`~scipy.optimize.anderson`.
+That is exactly what we do in the {autolink}`~pymf.solvers.solver` function, although we also provide the option to use a custom optimizer.

docs/source/documentation/mf_notes.md

@@ -31,24 +31,25 @@ Such a task is often infeasible for large systems and one needs to resort to app
 ## Mean-field approximaton
 
 The first-order perturbative approximation to the interacting Hamiltonian is the Hartree-Fock approximation also known as the mean-field approximation.
-The mean-field approximates the quartic term $\hat{V}$ in {eq}`hamiltonian` as a sum of bilinear terms weighted by the expectation values the remaining operators:
+The mean-field approximates the quartic term $\hat{V}$ in {eq}`hamiltonian` as a sum of bilinear terms weighted by the expectation values of the remaining operators:
 :::{math}
 :label: mf_approx
 \hat{V} \approx \hat{V}_{\text{MF}} \equiv \sum_{ij} v_{ij} \left[
-\braket{c_i^\dagger c_i} c_j^\dagger c_j - \braket{c_i^\dagger c_j} c_j^\dagger c_i \right]
+\braket{c_i^\dagger c_i} c_j^\dagger c_j - \braket{c_i^\dagger c_j} c_j^\dagger c_i \right],
 :::
-we neglect the superconducting pairing and constant offset terms.
-The expectation value terms  $\langle c_i^\dagger c_j \rangle$ are due to the ground-state density matrix and therefore act as an effective field acting on the system.
-The ground-state density matrix reads:
+where we neglect the constant offset terms and the superconducting pairing (for now).
+The expectation value terms  $\langle c_i^\dagger c_j \rangle$ are due to the ground state density matrix and act as an effective field on the system.
+The ground state density matrix reads:
 :::{math}
 :label: density
 \rho_{ij} \equiv \braket{c_i^\dagger c_j } = \text{Tr}\left(e^{-\beta \left(\hat{H_0} + \hat{V}_{\text{MF}} - \mu \hat{N} \right)} c_i^\dagger c_j\right),
 :::
 where $\beta = 1/ (k_B T)$ is the inverse temperature, $\mu$ is the chemical potential, and $\hat{N} = \sum_i c_i^\dagger c_i$ is the number operator.
+Currently, we neglect thermal effects and so $\beta \to \infty$.
 
 ## Finite tight-binding grid
 
-To simplify the mean-field Hamiltonian, we assume a finite, normalised orthogonal tight-binding grid defined by the single-particle basis states:
+To simplify the mean-field Hamiltonian, we assume a finite, normalised, orthogonal tight-binding grid defined by the single-particle basis states:
 
 $$
 \ket{n} = c^\dagger_n\ket{\text{vac}}
@@ -71,7 +72,7 @@ $$
 n \to n, R_n.
 $$
 
-Because of the translationaly invariance, the physical properties of the system are independent of the absolute unit cell position $R_n$ and rather depend on the relative position between the two unit cells $R_{nm} = R_n - R_m$:
+Because of the translational invariance, the physical properties of the system are independent of the absolute unit cell position $R_n$, but rather depend on the relative position between the two unit cells $R_{nm} = R_n - R_m$:
 
 $$
 \rho_{mn} \to \rho_{mn}(R_{mn}).