merge from main

.gitlab-ci.yml

@@ -5,7 +5,7 @@ before_script:
   - export MAMBA_ROOT_PREFIX=micromamba
   - micromamba create -yf docs/environment.yml
   - eval "$(micromamba shell hook --shell bash)"
-  - micromamba activate pymf-docs
+  - micromamba activate meanfi-docs
 
 run tests:
   script:
@@ -65,3 +65,35 @@ run pre-commit:
     paths:
       - .pre-commit-cache
       - micromamba
+
+prepare zips:
+  image: gitlab.kwant-project.org:5005/qt/research-docker
+  before_script: []
+  when: manual
+  script:
+    - zip -r zenodo.zip *
+  artifacts:
+    paths:
+      - zenodo.zip
+
+publish to test pypi:
+  needs:
+    - run tests
+    - run docs
+  rules:
+    - if: '$CI_COMMIT_TAG =~ /^v[0-9]+\.[0-9]+\.[0-9]+.*\+test$/' # vX.Y.Z.post1+test
+  script:
+    - micromamba install -c conda-forge -y hatch hatch-vcs
+    - hatch build
+    - hatch publish -u __token__ -a $PYPI_TEST_TOKEN -r test
+
+publish to pypi:
+  needs:
+    - run tests
+    - run docs
+  rules:
+    - if: '$CI_COMMIT_TAG =~ /^v[0-9]+\.[0-9]+\.[0-9]+[^+]*$/' # No +test
+  script:
+    - micromamba install -c conda-forge -y hatch hatch-vcs
+    - hatch build
+    - hatch publish -u __token__ -a $PYPI_TOKEN

.gitmodules

-[submodule "data"]
-	url = ../kwant-scf-data
-	path = data

.mailmap

+Kostas Vilkelis <kostasvilkelis@gmail.com>
+R. Johanna Zijderveld <johanna@zijderveld.de>
+Anton R. Akhmerov <meanfi@antonakhmerov.org>
+Antonio L.R. Manesco <am@antoniomanesco.org>

.readthedocs.yaml

+version: 2
+
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "mambaforge-4.10"
+
+conda:
+  environment: docs/environment.yml
+
+sphinx:
+   configuration: docs/source/conf.py

AUTHORS.md

+# MeanFi authors
+
+## Current meanfi maintainers
+- Kostas Vilkelis
+- R. Johanna Zijderveld
+- Anton R. Akhmerov
+- Antonio L.R. Manesco
+
+## Other contributors
+- Isidora Araya Day
+- José L. Lado
+
+## Funding
+
+The project was developed in [Delft University of
+Technology](https://www.tudelft.nl/en/).
+We acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement No. 828948
+(AndQC) and [Netherlands
+Organisation for Scientific Research (NWO/OCW)](https://www.nwo.nl/) 016.Vidi.189.180 grant.

LICENSE

+BSD 2-Clause License
+
+Copyright (c) 2024, MeanFi authors
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

README.md

-# Project Name
+# `MeanFi`
 
-## Research Goal
+## What is `MeanFi`?
 
-## Research Plan
+`MeanFi` is a Python package that performs self-consistent Hartree-Fock calculations on tight-binding models.
+It aims to find the groundstate of a Hamiltonian with density-density interactions
 
-## Working on this project
-Configure the project by running
+$$
+\hat{H} = \hat{H_0} + \hat{V} = \sum_{ij} h_{ij} c^\dagger_{i} c_{j} + \frac{1}{2} \sum_{ij} v_{ij} \hat{n}_i \hat{n}_j,
+$$
 
-    ./setup
+and computes the mean-field correction $\hat{V}_{\text{MF}}$ which approximates the interaction term:
+
+$$
+\hat{V} \approx \hat{V}_{\text{MF}} = \sum_{ij} \tilde{v}_{ij} c^\dagger_{i} c_{j}.
+$$
+
+For more details, refer to the [theory overview](docs/source/documentation/mf_notes.md) and [algorithm description](docs/source/documentation/algorithm.md).
+
+## How to use `MeanFi`?
+
+The calculation of a mean-field Hamiltonian is a simple 3-step process:
+
+1. **Define**
+
+    To specify the interacting problem, use a `Model` object which collects:
+    - Non-interacting Hamiltonian as a tight-binding dictionary.
+    - Interaction Hamiltonian as a tight-binding dictionary.
+    - Particle filling number in the unit cell.
+2. **Guess**
+
+    Construct a starting guess for the mean-field correction.
+
+3. **Solve**
+
+    Solve for the mean-field correction using the `solver` function and add it to the non-interacting part to obtain the total mean-field Hamiltonian.
+
+```python
+import meanfi
+
+#Define
+h_0 = {(0,) : onsite, (1,) : hopping, (-1,) : hopping.T.conj()}
+h_int = {(0,) : onsite_interaction}
+model = meanfi.Model(h_0, h_int, filling=2)
+
+#Guess
+guess = meanfi.guess_tb(guess_hopping_keys, ndof)
+
+#Solve
+mf_correction = meanfi.solver(model, guess)
+h_mf = meanfi.add_tb(h_0, mf_correction)
+```
+
+For more details and examples on how to use the package, we refer to the [tutorials](docs/source/tutorial/hubbard_1d.md).
+
+## Why `MeanFi`?
+
+Here is why you should use `MeanFi`:
+
+* Simple
+
+    The workflow is straightforward.
+    Interface with `Kwant` allows easy creation of complicated tight-binding systems and interactions.
+
+* Extensible
+
+    `MeanFi`'s code is structured to be easy to understand, modify and extend.
+
+* Optimized numerical workflow
+
+    Introduces minimal overhead to the calculation of the mean-field Hamiltonian.
+
+
+## What `MeanFi` doesn't do (yet)
+
+Here are some features that are not yet implemented but are planned for future releases:
+
+- **Superconductive order parameters**. Mean-field Hamiltonians do not include pairing terms.
+- **General interactions**. We allow only density-density interactions (e.g. Coulomb) which can be described by a second-order tensor.
+- **Temperature effects**. Density matrix calculations are done at zero temperature.
+
+## Installation
+
+```
+pip install meanfi
+```
+
+## Citing `MeanFi`
+
+If you have used `MeanFi` for work that has led to a scientific publication, please cite us as:
+
+```bibtex
+@misc{meanfi,
+  author = {Vilkelis, Kostas and Zijderveld,  R. Johanna and Akhmerov, Anton R. and Manesco, Antonio L.R.},
+  doi = {10.5281/zenodo.11149850},
+  month = {5},
+  title = {MeanFi},
+  year = {2024}
+}
+```

docs/environment.yml

-name: pymf-docs
+name: meanfi-docs
 
 channels:
   - conda-forge

docs/source/AUTHORS.md

+../../AUTHORS.md
\ No newline at end of file

docs/source/_static/css/custom.css

+html[data-theme="light"] {
+    --pst-color-primary: rgb(4, 87, 13);
+    --pst-color-secondary: rgb(66, 96, 150);
+    --pst-color-border-tippy: rgb(0, 0, 0);
+    --pst-color-inline-code-links: rgb(4, 87, 13);
+    --pst-color-link-hover: var(--pst-color-secondary)
+}
+
+html[data-theme="dark"] {
+    --pst-color-primary: rgb(176, 217, 162);
+    --pst-color-secondary: rgb(125, 166, 244);
+    --pst-color-border-tippy: white;
+    --pst-color-inline-code-links: var(--pst-color-primary);
+    --pst-color-link-hover: var(--pst-color-secondary)
+}
+
+.tippy-box {
+    background-color:var(--pst-color-surface);
+    color:var(--pst-color-text-base);
+    border: 1px solid var(--pst-color-border-tippy);
+}

docs/source/conf.py

@@ -13,9 +13,9 @@
 import os
 import sys
 
-import pymf  # noqa: F401
+import meanfi  # noqa: F401
 
-package_path = os.path.abspath("../pymf")
+package_path = os.path.abspath("../meanfi")
 # Suppress superfluous frozen modules warning.
 os.environ["PYDEVD_DISABLE_FILE_VALIDATION"] = "1"
 sys.path.insert(0, package_path)
@@ -23,14 +23,14 @@ sys.path.insert(0, package_path)
 
 # -- Project information -----------------------------------------------------
 
-project = "pymf"
-copyright = "2024, pymf developers"
-author = "pymf developers"
-gitlab_url = "https://gitlab.kwant-project.org/qt/pymf"
+project = "MeanFi"
+copyright = "2024, MeanFi developers"
+author = "MeanFi developers"
+gitlab_url = "https://gitlab.kwant-project.org/qt/meanfi"
 
 # The full version, including alpha/beta/rc tags
-release = pymf.__version__
-major, minor = pymf.__version_tuple__[:2]
+release = meanfi.__version__
+major, minor = meanfi.__version_tuple__[:2]
 version = f"{major}.{minor}"
 
 # -- General configuration ---------------------------------------------------
@@ -71,6 +71,8 @@ intersphinx_mapping = {
 
 default_role = "autolink"
 
+latex_elements = {"extrapackages": r"\usepackage{braket}"}
+
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ["_templates"]
 
@@ -86,9 +88,10 @@ autoclass_content = "both"
 # a list of builtin themes.
 #
 html_theme = "sphinx_book_theme"
+html_title = "MeanFi"
 
 html_theme_options = {
-    "repository_url": "https://gitlab.kwant-project.org/qt/kwant-scf",
+    "repository_url": "https://gitlab.kwant-project.org/qt/meanfi",
     "use_repository_button": True,
     "use_issues_button": True,
     "use_edit_page_button": True,
@@ -112,5 +115,5 @@ html_theme_options = {
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-# html_static_path = ["_static"]
-# html_css_files = ["local.css"]
+html_static_path = ["_static"]
+html_css_files = ["css/custom.css"]

docs/source/documentation/algorithm.md

+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.4
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+# Algorithm overview
+
+## Self-consistent mean-field loop
+
+To calculate the mean-field interaction in {eq}`mf_infinite`, we require the ground-state density matrix $\rho_{mn}(R)$.
+However, {eq}`density` is a function of the mean-field interaction $\hat{V}_{\text{MF}}$ itself.
+Therefore, we need to solve for both self-consistently.
+
+A single iteration of this self-consistency loop is a function that computes a new mean-field correction from a given one:
+
+$$
+\text{MF}(\hat{V}_{\text{init, MF}}) \to \hat{V}_{\text{new, MF}},
+$$
+
+which is defined in {autolink}`~meanfi.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. ({autolink}`~meanfi.mf.density_matrix`) Compute the ground-state density matrix $\rho_{mn}(R)$:
+    1. ({autolink}`~meanfi.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}`~meanfi.mf.fermi_on_kgrid`) Calculate the fermi level given the desired filling of the unit cell.
+    4.  ({autolink}`~meanfi.mf.density_matrix_kgrid`) Calculate the density matrix $\rho_{mn}(k)$ using the eigenvectors and the fermi level.
+    5. ({autolink}`~meanfi.tb.transforms.kgrid_to_tb`) Inverse Fourier transform the density matrix to real-space $\rho_{mn}(k) \to \rho_{mn}(R)$.
+3. ({autolink}`~meanfi.mf.meanfield`) Calculate the new mean-field correction $\hat{V}_{\text{new, MF}}(R)$ using {eq}`mf_infinite`.
+
+## 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:
+
+$$
+f : \hat{V}_{\text{MF}} \to f(\hat{V}_{\text{MF}}) \in \mathbb{R}^N.
+$$
+
+In the code, $f$ corresponds to the {autolink}`~meanfi.params.rparams.tb_to_rparams` function (inverse is {autolink}`~meanfi.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 real and imaginary parts to form a real-valued vector.
+
+With this, we define the self-consistency criterion as a fixed-point problem:
+
+$$
+f(\text{MF}(\hat{V}_{\text{MF}})) = f(\hat{V}_{\text{MF}}).
+$$
+
+Instead of solving the fixed point problem, we rewrite it as the difference of the two successive self-consistent mean-field iterations in {autolink}`~meanfi.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}`~meanfi.solvers.solver` function, although we also provide the option to use a custom optimizer.

docs/source/documentation/meanfi.md

+# Package reference
+
+## Interactive problem definition
+
+To define the interactive problem, we use the following class:
+
+```{eval-rst}
+.. autoclass:: meanfi.model.Model
+   :members: mfield
+```
+
+## Mean-field and density matrix
+
+```{eval-rst}
+.. automodule:: meanfi.mf
+   :members: meanfield, density_matrix, density_matrix_kgrid, fermi_on_kgrid
+   :show-inheritance:
+```
+
+## Observables
+
+```{eval-rst}
+.. automodule:: meanfi.observables
+   :members: expectation_value
+   :show-inheritance:
+```
+
+## Solvers
+
+```{eval-rst}
+.. automodule:: meanfi.solvers
+   :members: solver, cost
+   :show-inheritance:
+```
+
+## Tight-binding dictionary
+
+### Manipulation
+
+```{eval-rst}
+.. automodule:: meanfi.tb.tb
+   :members: add_tb, scale_tb
+   :show-inheritance:
+```
+
+### Brillouin zone transformations
+
+```{eval-rst}
+.. automodule:: meanfi.tb.transforms
+   :members:
+   :show-inheritance:
+```
+
+### Parametrisation
+
+```{eval-rst}
+.. automodule:: meanfi.params.rparams
+   :members:
+   :show-inheritance:
+```
+
+### Utility functions
+
+```{eval-rst}
+.. automodule:: meanfi.tb.utils
+   :members:
+   :show-inheritance:
+```
+
+## `kwant` interface
+
+```{eval-rst}
+.. automodule:: meanfi.kwant_helper.utils
+   :members:
+   :show-inheritance:
+```

docs/source/documentation/mf_notes.md

+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.4
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+# Theory overview
+
+## Interacting problems
+
+In physics, one often encounters problems where a system of multiple particles interacts with each other.
+In this package, we consider a general electronic system with density-density interparticle interaction:
+
+:::{math}
+:label: hamiltonian
+\hat{H} = \hat{H_0} + \hat{V} = \sum_{ij} h_{ij} c^\dagger_{i} c_{j} + \frac{1}{2} \sum_{ij} v_{ij} c_i^\dagger c_j^\dagger c_j c_i
+:::
+
+where $c_i^\dagger$ and $c_i$ are the creation and annihilation operators respectively for fermion in state $i$.
+The first term $\hat{H_0}$ is the non-interacting Hamiltonian which by itself is straightforward to solve on a single-particle basis by direct diagonalizations made easy through packages such as [kwant](https://kwant-project.org/).
+The second term $\hat{V}$ is the density-density interaction term between two particles, for example Coulomb interaction.
+To solve the interacting problem exactly, one needs to diagonalize the full Hamiltonian $\hat{H}$ in the many-particle basis which grows exponentially with the number of particles.
+Such a task is often infeasible for large systems and one needs to resort to approximations.
+
+## Mean-field approximation
+
+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 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],
+:::
+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 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:
+
+$$
+\ket{n} = c^\dagger_n\ket{\text{vac}}
+$$
+
+where $\ket{\text{vac}}$ is the vacuum state.
+We project our mean-field interaction in {eq}`mf_approx` onto the tight-binding grid:
+
+:::{math}
+:label: mf_finite
+V_{\text{MF}, nm} = \braket{n | \hat{V}_{\text{MF}} | m} =  \sum_{i} \rho_{ii} v_{in} \delta_{nm} - \rho_{mn} v_{mn},
+:::
+where $\delta_{nm}$ is the Kronecker delta function.
+
+## Infinite tight-binding grid
+
+In the limit of a translationally invariant system, the index $n$ that labels the basis states partitions into two independent variables: the unit cell internal degrees of freedom (spin, orbital, sublattice, etc.) and the position of the unit cell $R_n$:
+
+$$
+n \to n, R_n.
+$$
+
+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}).
+$$
+
+That allows us to re-write the mean-field interaction in {eq}`mf_finite` as:
+
+:::{math}
+:label: mf_infinite
+V_{\text{MF}, nm} (R) =  \sum_{i} \rho_{ii} (0) v_{in} (0) \delta_{nm} \delta(R) - \rho_{mn}(R) v_{mn}(R),
+:::
+
+where now indices $i, n, m$ label the internal degrees of freedom of the unit cell and $R$ is the relative position between the two unit cells in terms of the lattice vectors.

docs/source/index.md

@@ -11,29 +11,26 @@ kernelspec:
   name: python3
 ---
 
-# pymf
-
 ```{toctree}
 :hidden:
 :maxdepth: 1
 :caption: Tutorials
 
-mf_notes.md
+tutorial/hubbard_1d.md
+tutorial/graphene_example.md
 ```
 
-## What is pymf?
-
-## Why pymf?
-
-## How does pymf work?
-
-## What does pymf not do yet?
-
-## Installation
+```{toctree}
+:hidden:
+:maxdepth: 1
+:caption: Documentation
 
-```bash
-pip install .
+documentation/mf_notes.md
+documentation/algorithm.md
+documentation/meanfi.md
+AUTHORS.md
 ```
-## Citing
 
-## Contributing
+```{include} ../../README.md
+:relative-docs: docs/source/
+```

docs/source/mf_notes.md

-# Self-consistent mean field algorithm
-
-## Mean-field approximation
-
-The full hamiltonian is:
-
-$$
-\hat{H} = \hat{H_0} + \hat{V} = \sum_{ij} h_{ij} c^\dagger_{i} c_{j} + \frac{1}{2} \sum_{ijkl} v_{ijkl} c_i^\dagger c_j^\dagger c_l c_k
-$$
-
-We assume the dominant part of the ground state wavefunction comes from $\hat{H}_0$. Let's assume $b_i$ operators diagonalize the unperturbed hamiltonian:
-
-$$
-c_i^\dagger = \sum_{k} U_{ik} b_k^\dagger,
-$$
-
-such that the unperturbed groundstate wavefunction is:
-
-$$
-| 0 \rangle = \Pi_{E_i \leq \mu } b_i^\dagger |\textrm{vac}\rangle.
-$$
-
-Based on this definition, we define the normal ordering operator $:ABC...:$ such that it fulfills:
-
-$$
-:ABC...: | 0 \rangle = 0
-$$
-
-which practically means it orders $b_i$ operators based on whether its above or below the Fermi level $\mu$.
-
-Under this definition of normal ordering, we define the Wick's expansion of the interaction term:
-
-$$
-\begin{multline}
-c_i^\dagger c_j^\dagger c_l c_k = :c_i^\dagger c_j^\dagger c_l c_k: \\+  \overline{c_i^\dagger c_j^\dagger} :c_l c_k: + \overline{c_i^\dagger c_k} :c_j^\dagger c_l: - \overline{c_i^\dagger c_l} :c_j^\dagger c_k: + \overline{c_l c_k} :c_i^\dagger c_j^\dagger: - \overline{c_j^\dagger c_k} :c_i^\dagger c_l: + \overline{c_j^\dagger c_l} :c_i^\dagger c_k: \\
-\overline{c_i^\dagger c_j^\dagger} \overline{c_l c_k} - \overline{c_i^\dagger c_l} \overline{c_j^\dagger c_k} + \overline{c_i^\dagger c_k} \overline{c_j^\dagger c_l}
-\end{multline}
-$$
-
-where the overline defines Wick's contraction:
-
-$$
-\overline{AB} = AB - :AB:.
-$$
-
-The expectation value of the interaction with respect to the $| 0 \rangle$ is:
-
-$$
-\langle 0 | c_i^\dagger c_j^\dagger c_l c_k | 0 \rangle = \langle 0 | \overline{c_i^\dagger c_j^\dagger} \overline{c_l c_k} - \overline{c_i^\dagger c_l} \overline{c_j^\dagger c_k} + \overline{c_i^\dagger c_k} \overline{c_j^\dagger c_l}  | 0 \rangle
-$$
-
-where we can forget about all the normal ordered states since those give zero acting on the unperturbed groundstate. To evaluate this further, we utilize the mean-field approximation:
-
-$$
-A B \approx \langle A \rangle B + A \langle B \rangle - \langle A \rangle \langle B \rangle
-$$
-
-onto the contractions such that we get:
-
-$$
-\langle  c_i^\dagger c_j^\dagger c_l c_k \rangle \approx \langle c_i^\dagger c_j^\dagger \rangle \langle  c_l c_k \rangle + \langle c_i^\dagger c_k \rangle  \langle c_j^\dagger c_l \rangle - \langle c_i^\dagger c_l \rangle \langle c_j^\dagger c_k \rangle
-$$
-
-note $\langle A B \rangle \approx \langle A \rangle  \langle B \rangle$ assuming mean-field.
-
-To consider excitations from the groundstate, we make use of the mean-field approximation defined above:
-
-$$
-\begin{multline}
-c_i^\dagger c_j^\dagger c_l c_k \approx \\
-\langle c_i^\dagger c_j^\dagger \rangle c_l c_k + \langle c_i^\dagger c_k \rangle c_j^\dagger c_l - \langle c_i^\dagger c_l \rangle c_j^\dagger c_k + \langle c_l c_k \rangle c_i^\dagger c_j^\dagger - \langle c_j^\dagger c_k \rangle c_i^\dagger c_l + \langle c_j^\dagger c_l \rangle c_i^\dagger c_k + \\
-\langle c_i^\dagger c_j^\dagger \rangle \langle  c_l c_k \rangle + \langle c_i^\dagger c_k \rangle  \langle c_j^\dagger c_l \rangle - \langle c_i^\dagger c_l \rangle \langle c_j^\dagger c_k \rangle
-\end{multline}
-$$
-
-Where we made use of the following operations:
-
-$$
-:c_i^\dagger c_j^\dagger c_l c_k: \approx 0
-$$
-
-$$
-\overline{c_i^\dagger c_k} \overline{c_j^\dagger c_l} \approx \langle \overline{c_i^\dagger c_k} \rangle \overline{c_j^\dagger c_l} + \overline{c_i^\dagger c_k} \langle \overline{c_j^\dagger c_l} \rangle - \langle \overline{c_i^\dagger c_k} \rangle \langle \overline{c_j^\dagger c_i} \rangle =  \langle c_i^\dagger c_k \rangle \overline{c_j^\dagger c_l} + \overline{c_i^\dagger c_k} \langle c_j^\dagger c_l \rangle - \langle c_i^\dagger c_k \rangle \langle c_j^\dagger c_l \rangle
-$$
-
-$$
-\overline{c_i^\dagger c_k} :c_j^\dagger c_l: \approx \langle \overline{c_i^\dagger c_k} \rangle  :c_j^\dagger c_l: +  \overline{c_i^\dagger c_k} \langle :c_j^\dagger c_l: \rangle - \langle \overline{c_i^\dagger c_k} \rangle \langle :c_j^\dagger c_l: \rangle = \langle \overline{c_i^\dagger c_k} \rangle  :c_j^\dagger c_l:
-$$
-
-
-$$
-\langle \overline{c_i^\dagger c_k} \rangle  = \langle c_i^\dagger c_k - :c_i^\dagger c_k: \rangle  = \langle c_i^\dagger c_k \rangle
-$$
-
-
-Without any superconducting terms, the form simplifies to:
-
-$$
-\begin{multline}
-c_i^\dagger c_j^\dagger c_l c_k \approx
-\langle c_i^\dagger c_k \rangle c_j^\dagger c_l - \langle c_i^\dagger c_l \rangle c_j^\dagger c_k - \langle c_j^\dagger c_k \rangle c_i^\dagger c_l + \langle c_j^\dagger c_l \rangle c_i^\dagger c_k + \\
-\langle c_i^\dagger c_k \rangle  \langle c_j^\dagger c_l \rangle - \langle c_i^\dagger c_l \rangle \langle c_j^\dagger c_k \rangle
-\end{multline}
-$$
-
-## Finite size
-
-### Coulomb interaction
-
-We simplify the interaction term through the MF approximation to get:
-
-$$
-V = \frac{1}{2}\sum_{ijkl} v_{ijkl} c_i^{\dagger} c_j^{\dagger} c_l c_k
-\approx
-\frac12 \sum_{ijkl} v_{ijkl} \left[ \langle c_i^{\dagger} c_k \rangle c_j^{\dagger} c_l - \langle c_j^{\dagger} c_k \rangle c_i^{\dagger} c_l - \langle c_i^{\dagger} c_l \rangle c_j^{\dagger} c_k + \langle c_j^{\dagger} c_l \rangle c_i^{\dagger} c_k \right]
-$$
-(assuming no superconductivity)
-
-and an additional constant part:
-
-$$
-V_0 =  \frac{1}{2} \sum_{ijkl} v_{ijkl} \left(\langle c_j^{\dagger} c_l \rangle \langle c_i^{\dagger} c_k \rangle - \langle c_j^{\dagger} c_k \rangle \langle c_i^{\dagger} c_l \rangle \right).
-$$
-
-The interaction reads:
-
-$$
-v_{ijkl} = \iint w_{i}^*(r) w_{j}^*(r') V(r, r') w_{k}(r) w_l(r') dr dr' = \\
-\iint  V(|r - r'|) w_{i}^*(r)w_{k}(r) w_{j}^*(r')  w_l(r') dr dr'
-$$
-
-whereas $w_i$ is a wannier function on site i (and corresponding dof). Whenever one interchanges $i \to j, k \to l$, the Coulomb term is preserved $v_{ijkl} = v_{jilk}$
-
-To make things more understandable, we are also going to explicitly split up position and spin indices: $i \to i \times \sigma$. In this notation, the Coulomb integral reads:
-
-$$
-v_{ijkl}^{\sigma_i \sigma_j \sigma_k \sigma_l} =
-\iint V(|r - r'|) w_{i\times\sigma_i}^{*} (r)w_{k \times \sigma_k}(r) w_{j \times \sigma_j}^{*}(r')  w_{l\times \sigma_l}(r') dr dr' \delta_{\sigma_i \sigma_k} \delta_{\sigma_{j} \sigma_l}
-$$
-
-On a fine tight-binding model, we have:
-
-$$
-v_{ijkl}^{\sigma_i \sigma_j \sigma_k \sigma_l} = v_{ij} \delta_{ik} \delta_{jl} \delta_{\sigma_i \sigma_k} \delta_{\sigma_j \sigma_l}
-$$
-
-where $v_{ij} = V(r_i, r_j)$.
-
-We shall re-define $i$ index to absorb spin:
-
-$$
-\delta_{ik} \times \delta_{\sigma_{i} \sigma_{k}} \to \delta_{ik}
-$$
-
-in this notation the above reads:
-
-$$
-v_{ijkl} = v_{ij} \delta_{ik} \delta_{jl}
-$$
-
-The mean-field terms are:
-
-$$
-\langle c_i^{\dagger} c_j\rangle = \langle \Psi_F|c_i^{\dagger} c_j | \Psi_F \rangle
-$$
-
-whereas $|\Psi_F \rangle = \Pi_{i=0}^{N_F} b_i^\dagger |0\rangle$. To make sense of things, we need to transform between $c_i$ basis (position + internal dof basis) into the $b_i$ basis (eigenfunction of a given mean-field Hamiltonian):
-
-$$
-c_i^\dagger = \sum_{k} U_{ik} b_k^\dagger
-$$
-
-where
-
-$$
-U_{ik} = \langle{i|\psi_k} \rangle.
-$$
-
-That gives us:
-
-
-$$
-c_i^{\dagger} c_j = \sum_{k, l} U_{ik} U_{lj}^* b_k^\dagger b_{l}
-$$
-
-and its expectation value gives us the mean-field ... field $F_{ij}$:
-
-$$
-F_{ij} = \langle c_i^{\dagger} c_j\rangle =  \sum_{k, l} U_{ik} U_{lj}^* \langle \Psi_F| b_k^\dagger b_{l}| \Psi_F \rangle =  \sum_{k} U_{ik} U_{kj}^{*}
-$$
-
-whereas I assumed the indices for wavefunctions $k,l$ are ordered in terms of increasing eigenvalue. We pop that into the definition of the mean-field correction $V$:
-
-
-$$
-\begin{multline}
-V_{nm} = \frac12 \sum_{ijkl} v_{ijkl} \langle n| \left[ \langle c_i^{\dagger} c_k \rangle c_j^{\dagger} c_l - \langle c_j^{\dagger} c_k \rangle c_i^{\dagger} c_l - \langle c_i^{\dagger} c_l \rangle c_j^{\dagger} c_k + \langle c_j^{\dagger} c_l \rangle c_i^{\dagger} c_k \right] |m\rangle = \\
- \frac12 \sum_{ijkl} v_{ijkl} \left[ +\delta_{jn}\delta_{lm} F_{ik} -\delta_{in}\delta_{lm} F_{jk} -\delta_{jn}\delta_{km} F_{il} + \delta_{in}\delta_{km} F_{jl} \right] = \\
-\frac12 \left[ \sum_{ik} v_{inkm} F_{ik} - \sum_{jk} v_{njkm} F_{jk} - \sum_{il} v_{inml} F_{il} + \sum_{jl} v_{njml} F_{jl} \right] = \\
--\sum_{ij} F_{ij} \left(v_{inmj} - v_{injm} \right)
-\end{multline}
-$$
-
-where I used the $v_{ijkl} = v_{jilk}$ symmetry from Coulomb.
-
-For a tight-binding grid, the mean-field correction simplifies to:
-
-$$
-\begin{multline}
-V_{nm} = - \sum_{ij} F_{ij} \left(v_{inmj} - v_{injm} \right) = \\
--\sum_{ij}F_{ij} v_{in} \delta_{im} \delta_{nj} + \sum_{ij}F_{ij} v_{in} \delta_{ij} \delta_{nm} = \\
--F_{mn} v_{mn} + \sum_{i} F_{ii} v_{in} \delta_{nm}
-\end{multline}
-$$
-
-the first term is the exchange interaction whereas the second one is the direct interaction.
-
-Similarly, the constant offset term reads:
-
-$$
-\begin{multline}
-V_0 = \frac{1}{2} \sum_{ijkl} v_{ijkl} \left(F_{jl} F_{ik} - F_{jk} F_{il} \right) = \\
-\frac{1}{2} \sum_{ijkl} v_{ij} \delta_{ik} \delta_{jl} \left(F_{jl} F_{ik} - F_{jk} F_{il}\right) \\
-= \frac{1}{2} \sum_{ij} v_{ij} \left(F_{ii} F_{jj} - F_{ji} F_{ij}\right)
-\end{multline}
-$$
-
-where we identify the first term as the exchange (mixes indices) and the right one as the direct (diagonal in indices).
-
-## Translational Invariance
-
-The above assumed a finite tight-binding model - all $nm$-indices contain the position of all atoms (among other dof). In this section tho we want to consider an infinite system with translational invariance.
-
-To begin with we deconstruct a general matrix $O_{nm}$ into the cell degrees of freedom ($nm$) and the position of the the cell itself ($ij$):
-
-$$
-O_{nm} \to O^{ij}_{nm}
-$$
-
-and we will Fourier transform the upper indices into k-space:
-
-$$
-O_{mn}(k) = \sum_{ij} O_{nm}^{ij} e^{-i k (R_i-R_j)}
-$$
-
-where I assumed $O$ (and thus all operators I will consider here) is local and thus diagonal in k-space.
-
-Now lets rewrite our main result in the previous section using our new notation:
-
-$$
-V_{nm}^{ij} =-F_{mn}^{ij} v_{mn}^{ij} + \sum_{r,p} F_{pp}^{rr} v_{pn}^{rj} \delta_{nm} \delta^{ij}
-$$
-
-Lets first consider the second (direct) term. Lets express the corresponding $F$ term in k-space:
-
-$$
-F_{pp}^{rr} = \int e^{i k (R_r-R_r)} F_{pp}(k) dk = \int F_{pp}(k) dk
-$$
-
-Notice that in the final expression, there is no $rr$ dependence and thus this term is cell-periodic. Therefore, we shall redefine it as cell electron density $\rho$:
-$$
-F_{pp}^0 = F_{pp}(R = 0) = \int F_{pp}(k) dk
-$$
-
-Now since $\rho$ has no $r$ dependence, we can proceed with the sum:
-
-$$
-\sum_{r} v_{pn}^{rj} = \int v_{pn}(k) e^{ik R_j} \sum_{r} e^{-i k R_r} dk = \int v_{pn}(k) e^{ik R_j} \delta_{k, 0} dk = v_{pn}(0)
-$$
-
-We are finally ready to Fourier transform the main result. Invoking convolution theorem and the results above gives us:
-
-$$
-V_{nm}(k) = \sum_{p} F_{pp}^0 v_{pn}(0) \delta_{nm} -F_{mn}(k) \circledast v_{mn}(k) = V_n^D - F_{mn}(k) \circledast v_{mn}(k)
-$$
-
-which does make sense. The first term (direct) is a potential term coming from the mean-field and the second term (exchange) is purely responsible for the hopping.
-
-The constant offset is:
-$$
-V_0 = \frac{1}{2} \sum_{r,s} \rho_r v_{rs}(0) \rho_s- \\ \frac{1}{2} tr\left[\int_{BZ} \left(F \circledast v\right)(k) F(k) dk \right]
-$$
-
-## Short-Range interactions
-
-In the case of short-range interactions, it is much more convenient to go back to real space to both store objects and perform the operations. In real space the mean-field part of the Hamiltonian reads:
-
-$$
-V_{nm}(\mathbf{R}) = V_n^D \delta(\mathbf{R}) - F_{mn}(\mathbf{R}) v_{mn}(\mathbf{R})
-$$
-
-(the first term might need some prefactor from Fourier transformation)
-
-where $\mathbf{R}$ is a sequence of integers representing real-space unit cell indices.
-
-### Proposed Algorithm
-Given an initial Hamiltonian $H_0 (R)$ and the interaction term $v(R)$ in real-space, the mean-field algorithm is the following:
-
-0. Start with a mean-field guess in real-space: $V(R)$.
-1. Fourier transform tight-binding model and the mean-field in real space to a given k-grid: $H_0(R) + V(R) \to H_0(k) + V(k)$
-2. Diagonalize and evaluate the density matrix: $H_0(k) + V(k) \to F(k)$
-3. Fourier transform the density matrix back to real-space: $F(k) \to F(R)$
-4. Evaluate the new mean-field Hamiltonian $V(R)$ according to the equation above.
-5. Evaluate self-consistency metric (could be based either on $V(R/k)$ or $F(R/k)$). Based on that, either stop or go back to 1 with a modified $V(R)$ starting guess.

docs/source/tutorial/graphene_example.md

+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.4
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
+# Interacting graphene
+
+In the previous tutorial, we showed how to use `MeanFi` to solve a simple 1D Hubbard model with onsite interactions.
+In this tutorial, we will apply `MeanFi` to more complex system: graphene with onsite $U$ and nearest-neighbour $V$ interactions.
+The system is more complicated in every aspect: the lattice structure, dimension of the problem, complexity of the interactions.
+And yet, the workflow is the same as in the previous tutorial and remains simple and straightforward.
+
+## Building the system with `kwant`
+
+### Non-interacting part
+
+As in the previous tutorial, we could construct a tight-binding dictionary of graphene by hand, but instead it is much easier to use [`kwant`](https://kwant-project.org/) to build the system.
+For a more detailed explanation on `kwant` see the [tutorial](https://kwant-project.org/doc/1/tutorial/graphene).
+
+
+```{code-cell} ipython3
+import numpy as np
+import matplotlib.pyplot as plt
+import kwant
+
+import meanfi
+from meanfi.kwant_helper import utils
+
+s0 = np.identity(2)
+sx = np.array([[0, 1], [1, 0]])
+sy = np.array([[0, -1j], [1j, 0]])
+sz = np.diag([1, -1])
+
+# Create graphene lattice
+graphene = kwant.lattice.general([(1, 0), (1 / 2, np.sqrt(3) / 2)],
+                                 [(0, 0), (0, 1 / np.sqrt(3))], norbs=2)
+a, b = graphene.sublattices
+
+# Create bulk system
+bulk_graphene = kwant.Builder(kwant.TranslationalSymmetry(*graphene.prim_vecs))
+# Set onsite energy to zero
+bulk_graphene[a.shape((lambda pos: True), (0, 0))] = 0 * s0
+bulk_graphene[b.shape((lambda pos: True), (0, 0))] = 0 * s0
+# Add hoppings between sublattices
+bulk_graphene[graphene.neighbors(1)] = s0
+```
+
+The `bulk_graphene` object is a `kwant.Builder` object that represents the non-interacting graphene system.
+To convert it to a tight-binding dictionary, we use the {autolink}`~meanfi.kwant_helper.utils.builder_to_tb` function:
+
+```{code-cell} ipython3
+h_0 = utils.builder_to_tb(bulk_graphene)
+```
+
+### Interacting part
+
+We utilize `kwant` to build the interaction tight-binding dictionary as well.
+To define the interactions, we need to specify two functions:
+* `onsite_int(site)`: returns the onsite interaction matrix.
+* `nn_int(site1, site2)`: returns the interaction matrix between `site1` and `site2`.
+
+We feed these functions to the {autolink}`~meanfi.kwant_helper.utils.build_interacting_syst` function, which constructs the `kwant.Builder` object encoding the interactions.
+All we need to do is to convert this object to a tight-binding dictionary using the {autolink}`~meanfi.kwant_helper.utils.builder_to_tb` function.
+
+```{code-cell} ipython3
+def onsite_int(site, U):
+    return U * sx
+
+def nn_int(site1, site2, V):
+    return V * np.ones((2, 2))
+
+builder_int = utils.build_interacting_syst(
+    builder=bulk_graphene,
+    lattice=graphene,
+    func_onsite=onsite_int,
+    func_hop=nn_int,
+    max_neighbor=1
+)
+params = dict(U=0.2, V=1.2)
+h_int = utils.builder_to_tb(builder_int, params)
+```
+
+Because `nn_int` function returns the same interaction matrix for all site pairs, we set `max_neighbor=1` to ensure that the interaction only extends to nearest-neighbours and is zero for longer distances.
+
+## Computing expectation values
+
+As before, we construct {autolink}`~meanfi.model.Model` object to represent the full system to be solved via the mean-field approximation.
+We then generate a random guess for the mean-field solution and solve the system:
+
+```{code-cell} ipython3
+filling = 2
+model = meanfi.Model(h_0, h_int, filling=2)
+int_keys = frozenset(h_int)
+ndof = len(list(h_0.values())[0])
+guess = meanfi.guess_tb(int_keys, ndof)
+mf_sol = meanfi.solver(model, guess, nk=18)
+h_full = meanfi.add_tb(h_0, mf_sol)
+```
+
+To investigate the effects of interaction on systems with more than one degree of freedom, it is more useful to consider the expectation values of various operators which serve as order parameters.
+For example, we can compute the charge density wave (CDW) order parameter which is defined as the difference in the charge density between the two sublattices.
+
+To calculate operator expectation values, we first need to construct the density matrix via the {autolink}`~meanfi.mf.density_matrix` function.
+We then feed it into {autolink}`~meanfi.observables.expectation_value` function together with the operator we want to measure.
+In this case, we compute the CDW order parameter by measuring the expectation value of the $\sigma_z$ operator acting on the graphene sublattice degree of freedom.
+```{code-cell} ipython3
+cdw_operator = {(0, 0): np.kron(sz, np.eye(2))}
+
+rho, _ = meanfi.density_matrix(h_full, filling=filling, nk=40)
+rho_0, _ = meanfi.density_matrix(h_0, filling=filling, nk=40)
+
+cdw_order_parameter = meanfi.expectation_value(rho, cdw_operator)
+cdw_order_parameter_0 = meanfi.expectation_value(rho_0, cdw_operator)
+
+print(f"CDW order parameter for interacting system: {np.round(np.abs(cdw_order_parameter), 2)}")
+print(f"CDW order parameter for non-interacting system: {np.round(np.abs(cdw_order_parameter_0), 2)}")
+```
+
+We see that the CDW order parameter is non-zero only for the interacting system, indicating the presence of a CDW phase.
+
+## Graphene phase diagram
+
+In the remaining part of this tutorial, we will utilize all the tools we have developed so far to create a phase diagram for the graphene system.
+
+To identify phase changes, it is convenient to track the gap of the system as a function of $U$ and $V$.
+To that end, we first create a function that calculates the gap of the system given the tight-binding dictionary and the Fermi energy.
+
+```{code-cell} ipython3
+def compute_gap(h, fermi_energy=0, nk=100):
+    kham = meanfi.tb_to_kgrid(h, nk)
+    vals = np.linalg.eigvalsh(kham)
+
+    emax = np.max(vals[vals <= fermi_energy])
+    emin = np.min(vals[vals > fermi_energy])
+    return np.abs(emin - emax)
+```
+
+And proceed to compute the gap and the mean-field correction for a range of $U$ and $V$ values:
+
+```{code-cell} ipython3
+Us = np.linspace(0, 4, 10)
+Vs = np.linspace(0, 1.5, 10)
+
+gaps = []
+mf_sols = []
+for U in Us:
+    for V in Vs:
+        params = dict(U=U, V=V)
+        h_int = utils.builder_to_tb(builder_int, params)
+
+        model = meanfi.Model(h_0, h_int, filling=filling)
+        guess = meanfi.guess_tb(int_keys, ndof)
+        mf_sol = meanfi.solver(model, guess, nk=18)
+        mf_sols.append(mf_sol)
+
+        gap = compute_gap(meanfi.add_tb(h_0, mf_sol), fermi_energy=0, nk=100)
+        gaps.append(gap)
+gaps = np.asarray(gaps, dtype=float).reshape((len(Us), len(Vs)))
+mf_sols = np.asarray(mf_sols).reshape((len(Us), len(Vs)))
+
+plt.imshow(gaps.T, extent=(Us[0], Us[-1], Vs[0], Vs[-1]), origin='lower', aspect='auto')
+plt.colorbar()
+plt.xlabel('V')
+plt.ylabel('U')
+plt.title('Gap')
+plt.show()
+```
+
+This phase diagram has gap openings at the same places as shown in the [literature](https://arxiv.org/abs/1204.4531).
+
+We can now use the stored results in `mf_sols` to fully map out the phase diagram with order parameters.
+On top of the charge density wave (CDW), we also expect a spin density wave (SDW) in a different region of the phase diagram.
+We construct the SDW order parameter with the same steps as before, but now we need to sum over the expectation values of the three Pauli matrices to account for the $SU(2)$ spin-rotation symmetry.
+
+```{code-cell} ipython3
+s_list = [sx, sy, sz]
+cdw_list = []
+sdw_list = []
+for mf_sol in mf_sols.flatten():
+    rho, _ = meanfi.density_matrix(meanfi.add_tb(h_0, mf_sol), filling=filling, nk=40)
+
+    # Compute CDW order parameter
+    cdw_list.append(np.abs(meanfi.expectation_value(rho, cdw_operator))**2)
+
+    # Compute SDW order parameter
+    sdw_value = 0
+    for s_i in s_list:
+      sdw_operator_i = {(0, 0) : np.kron(sz, s_i)}
+      sdw_value += np.abs(meanfi.expectation_value(rho, sdw_operator_i))**2
+    sdw_list.append(sdw_value)
+
+cdw_list = np.asarray(cdw_list).reshape(mf_sols.shape)
+sdw_list = np.asarray(sdw_list).reshape(mf_sols.shape)
+```
+
+Finally, we can combine the gap, CDW and SDW order parameters into one plot.
+We naively do this by plotting the difference between CDW and SDW order parameters and indicate the gap with the transparency.
+
+```{code-cell} ipython3
+import matplotlib.ticker as mticker
+normalized_gap = gaps/np.max(gaps)
+plt.imshow((cdw_list - sdw_list).T, extent=(Us[0], Us[-1], Vs[0], Vs[-1]), origin='lower', aspect='auto', cmap="coolwarm", alpha=normalized_gap.T, vmin=-2.6, vmax=2.6)
+plt.colorbar(ticks=[-2.6, 0, 2.6], format=mticker.FixedFormatter(['SDW', '0', 'CDW']), label='Order parameter', extend='both')
+plt.xlabel('V')
+plt.ylabel('U')
+plt.show()
+```

docs/source/tutorial/hubbard_1d.md

+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.4
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+# 1D Hubbard model
+
+## Background physics
+
+To show the basic functionality of the package, we consider a simple interacting electronic system: a 1D chain of sites that allow nearest-neighbor tunneling with strength $t$ and on-site repulsion $U$ between two electrons if they are on the same site.
+Such a model is known as the 1D [Hubbard model](https://en.wikipedia.org/wiki/Hubbard_model) and is useful for understanding the onset of insulating phases in interacting metals.
+
+To begin, we first consider the second quantized form of the non-interacting Hamiltonian.
+Because we expect the interacting ground state to be antiferromagnetic, we build a two-atom cell and name the two sublattices $A$ and $B$.
+These sublattices are identical to each other in the non-interacting case $U=0$.
+The non-interacting Hamiltonian reads:
+
+$$
+\hat{H_0} = - t \sum_\sigma \sum_i \left(c_{i, B, \sigma}^{\dagger}c_{i, A, \sigma} + c_{i, A, \sigma}^{\dagger}c_{i+1, B, \sigma} + \textrm{h.c}\right).
+$$
+
+where $\textrm{h.c}$ is the hermitian conjugate, $\sigma$ denotes spin ($\uparrow$ or $\downarrow$) and $c_{i, A, \sigma}^{\dagger}$ creates an electron with spin $\sigma$ in unit cell $i$ of sublattice $A$.
+Next up, is the interacting part of the Hamiltonian:
+
+$$
+\hat{V} = U \sum_i \left(n_{i, A, \uparrow} n_{i, A, \downarrow} + n_{i, B, \uparrow} n_{i, B, \downarrow}\right).
+$$
+
+where $n_{i, A, \sigma} = c_{i, A, \sigma}^{\dagger}c_{i, A, \sigma}$ is the number operator for sublattice $A$ and spin $\sigma$.
+The total Hamiltonian is then $\hat{H} = \hat{H_0} + \hat{V}$.
+With the model defined, we can now proceed to input the Hamiltonian into the package and solve it using the mean-field approximation.
+
+## Problem definition
+
+### Non-interacting Hamiltonian
+
+First, let's get the basic imports out of the way.
+
+```{code-cell} ipython3
+import numpy as np
+import matplotlib.pyplot as plt
+import meanfi
+```
+Now let us translate the non-interacting Hamiltonian $\hat{H_0}$ defined above into the basic input format for the package: a **tight-binding dictionary**.
+The tight-binding dictionary is a python dictionary where the keys are tuples of integers representing the hopping vectors and the values are the hopping matrices.
+For example, a key `(0,)` represents the onsite term in one dimension and a key `(1,)` represents the hopping a single unit cell to the right.
+In two dimensions a key `(0,0)` would represent the onsite term and `(1,0)` would represent hopping to the right in the direction of the first reciprocal lattice vector.
+In the case of our 1D Hubbard model, we only have an onsite term and hopping a single unit cell to the left and right.
+Thus our non-interacting Hamiltonian  becomes:
+
+```{code-cell} ipython3
+hopp = np.kron(np.array([[0, 1], [0, 0]]), np.eye(2))
+h_0 = {(0,): hopp + hopp.T.conj(), (1,): hopp, (-1,): hopp.T.conj()}
+```
+Here `hopp` is the hopping matrix which we define as a kronecker product between sublattice and spin degrees of freedom: `np.array([[0, 1], [0, 0]])` corresponds to the hopping between sublattices and `np.eye(2)` leaves the spin degrees of freedom unchanged.
+In the corresponding tight-binding dictionary `h_0`, the key `(0,)` contains hopping within the unit cell and the keys `(1,)` and `(-1,)` correspond to the hopping between the unit cells to the right and left respectively.
+
+To verify the validity of `h_0`, we evaluate it in the reciprocal space using the {autolink}`~meanfi.tb.transforms.tb_to_kgrid`, then diagonalize it and plot the band structure:
+
+```{code-cell} ipython3
+nk = 50 # number of k-points
+ks = np.linspace(0, 2*np.pi, nk, endpoint=False)
+hamiltonians_0 = meanfi.tb_to_kgrid(h_0, nk)
+
+vals, vecs = np.linalg.eigh(hamiltonians_0)
+plt.plot(ks, vals, c="k")
+plt.xticks([0, np.pi, 2 * np.pi], ["$0$", "$\pi$", "$2\pi$"])
+plt.xlim(0, 2 * np.pi)
+plt.ylabel("$E - E_F$")
+plt.xlabel("$k / a$")
+plt.show()
+```
+
+which seems metallic as expected.
+
+### Interaction Hamiltonian
+
+We now proceed to define the interaction Hamiltonian $\hat{V}$.
+To achieve this, we utilize the same tight-binding dictionary format as before.
+Because the interaction Hamiltonian is on-site, it must be defined only for the key `(0,)` and only for electrons on the same sublattice with opposite spins.
+Based on the kronecker product structure we defined earlier, the interaction Hamiltonian is:
+
+```{code-cell} ipython3
+U = 2
+s_x = np.array([[0, 1], [1, 0]])
+h_int = {(0,): U * np.kron(np.eye(2), s_x),}
+```
+Here `s_x` is the Pauli matrix acting on the spin degrees of freedom, which ensures that the interaction is only between electrons with opposite spins whereas the `np.eye(2)` ensures that the interaction is only between electrons on the same sublattice.
+
+### Putting it all together
+
+To combine the non-interacting and interaction Hamiltonians, we use the {autolink}`~meanfi.model.Model` class.
+In addition to the Hamiltonians, we also need to specify the filling of the system --- the number of electrons per unit cell.
+```{code-cell} ipython3
+filling = 2
+full_model = meanfi.Model(h_0, h_int, filling)
+```
+
+The object `full_model` now contains all the information needed to solve the mean-field problem.
+
+## Solving the mean-field problem
+
+To find a mean-field solution, we first require a starting guess.
+In cases where the non-interacting Hamiltonian is highly degenerate, there exists several possible mean-field solutions, many of which are local and not global minima of the energy landscape.
+Therefore, the choice of the initial guess can significantly affect the final solution depending on the energy landscape.
+Here the problem is simple enough that we can generate a random guess for the mean-field solution through the {autolink}`~meanfi.tb.utils.guess_tb` function.
+It creates a random Hermitian tight-binding dictionary based on the hopping keys provided and the number of degrees of freedom within the unit cell.
+Because the mean-field solution cannot contain hoppings longer than the interaction itself, we use `h_int` keys as an input to {autolink}`~meanfi.tb.utils.guess_tb`.
+Finally, to solve the model, we use the {autolink}`~meanfi.solvers.solver` function which by default employes a root-finding [algorithm](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.anderson.html) to find a self-consistent mean-field solution.
+
+```{code-cell} ipython3
+filling = 2
+full_model = meanfi.Model(h_0, h_int, filling)
+guess = meanfi.guess_tb(frozenset(h_int), ndof=4)
+mf_sol = meanfi.solver(full_model, guess, nk=nk)
+```
+
+The {autolink}`~meanfi.solvers.solver` function returns only the mean-field correction to the non-interacting Hamiltonian in the same tight-binding dictionary format.
+To get the full Hamiltonian, we add the mean-field correction to the non-interacting Hamiltonian and plot the band structure just as before:
+
+```{code-cell} ipython3
+h_mf = meanfi.add_tb(h_0, mf_sol)
+
+hamiltonians = meanfi.tb_to_kgrid(h_mf, nk)
+vals, vecs = np.linalg.eigh(hamiltonians)
+plt.plot(ks, vals, c="k")
+plt.xticks([0, np.pi, 2 * np.pi], ["$0$", "$\pi$", "$2\pi$"])
+plt.xlim(0, 2 * np.pi)
+plt.ylabel("$E - E_F$")
+plt.xlabel("$k / a$")
+plt.show()
+```
+
+the band structure now shows a gap at the Fermi level, indicating that the system is in an insulating phase!
+
+
+We can go further and compute the gap for a wider range of $U$ values:
+
+```{code-cell} ipython3
+def compute_sol(U, h_0, nk, filling=2):
+    h_int = {
+        (0,): U * np.kron(np.eye(2), np.ones((2, 2))),
+    }
+    guess = meanfi.guess_tb(frozenset(h_int), len(list(h_0.values())[0]))
+    full_model = meanfi.Model(h_0, h_int, filling)
+    mf_sol = meanfi.solver(full_model, guess, nk=nk)
+    return meanfi.add_tb(h_0, mf_sol)
+
+
+def compute_gap(full_sol, nk_dense, fermi_energy=0):
+    h_kgrid = meanfi.tb_to_kgrid(full_sol, nk_dense)
+    vals = np.linalg.eigvalsh(h_kgrid)
+
+    emax = np.max(vals[vals <= fermi_energy])
+    emin = np.min(vals[vals > fermi_energy])
+    return np.abs(emin - emax)
+
+
+def compute_phase_diagram(
+    Us,
+    nk,
+    nk_dense,
+):
+    gaps = []
+    for U in Us:
+        full_sol = compute_sol(U, h_0, nk)
+        gaps.append(compute_gap(full_sol, nk_dense))
+
+    return np.asarray(gaps, dtype=float)
+
+Us = np.linspace(0, 4, 40, endpoint=True)
+gaps = compute_phase_diagram(Us=Us, nk=20, nk_dense=100)
+
+plt.plot(Us, gaps, c="k")
+plt.xlabel("$U / t$")
+plt.ylabel("$\Delta{E}/t$")
+plt.show()
+```
+
+We see that at around $U=1$ the gap opens up and the system transitions from a metal to an insulator.  In order to more accurately determine the size of the gap, we chose to use a denser k-grid for the diagonalization of the mean-field solution.

examples/codes

-../codes/
\ No newline at end of file