From 72f54c716774cc957d6bee262eceaac29785b6f0 Mon Sep 17 00:00:00 2001
From: Johanna <johanna@zijderveld.de>
Date: Tue, 7 May 2024 11:12:51 +0200
Subject: [PATCH] copy over everything from example.ipynb file
---
docs/source/hubbard_1d.md | 123 ++++++++++++++++++++++++++++++++++++++
1 file changed, 123 insertions(+)
create mode 100644 docs/source/hubbard_1d.md
diff --git a/docs/source/hubbard_1d.md b/docs/source/hubbard_1d.md
new file mode 100644
index 0000000..90ebb17
--- /dev/null
+++ b/docs/source/hubbard_1d.md
@@ -0,0 +1,123 @@
+---
+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
+
+```{code-cell} ipython3
+import numpy as np
+import matplotlib.pyplot as plt
+import numpy as np
+import matplotlib.pyplot as plt
+import pymf
+```
+
+To simulate infinite systems, we provide the corresponding tight-binding model.
+
+We exemplify this construction by computing the ground state of an infinite spinful chain with onsite interactions.
+
+Because the ground state is an antiferromagnet, so we must build a two-atom cell. We name the two sublattices, $A$ and $B$. The Hamiltonian in is:
+$$
+H_0 = \sum_i c_{i, B}^{\dagger}c_{i, A} + c_{i, A}^{\dagger}c_{i+1, B} + h.c.
+$$
+We write down the spinful by simply taking $H_0(k) \otimes \mathbb{1}$.
+
+
+To build the tight-binding model, we need to generate a Hamiltonian on a k-point and the corresponding hopping vectors to generate a guess. We then verify the spectrum and see that the bands indeed consistent of two bands due to the Brillouin zone folding.
+
+```{code-cell} ipython3
+# Set number of k-points
+nk = 100
+ks = np.linspace(0, 2*np.pi, nk, endpoint=False)
+hamiltonians_0 = transforms.tb_to_khamvector(h_0, nk, 1, ks=ks)
+
+# Perform diagonalization
+vals, vecs = np.linalg.eigh(hamiltonians_0)
+# Plot data
+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()
+```
+
+Here, in the workflow to find the ground state, we use a helper function to build the initial guess. because we don't need a dense k-point grid in the self-consistent loop, we compute the spectrum later on a denser k-point grid.
+
+Finally, we compute the eigen0alues for a set of Ualues of $U$. For this case, since the interaction is onsite only, the interaction matrix is simply
+$$
+H_{int} =
+\left(\begin{array}{cccc}
+ U & U & 0 & 0\\
+ U & U & 0 & 0\\
+ 0 & 0 & U & U\\
+ 0 & 0 & U & U
+\end{array}\right)~.
+$$
+
+```{code-cell} ipython3
+def compute_phase_diagram(
+ Us,
+ nk,
+ nk_dense,
+ filling=2,
+):
+ gap = []
+ vals = []
+ for U in tqdm(Us):
+ # onsite interactions
+ h_int = {
+ (0,): U * np.kron(np.ones((2, 2)), np.eye(2)),
+ }
+ guess = utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
+ full_model = Model(h_0, h_int, filling)
+ mf_sol = solver(full_model, guess, nk=nk)
+ hkfunc = transforms.tb_to_kfunc(add_tb(h_0, mf_sol))
+ ks_dense = np.linspace(0, 2 * np.pi, nk_dense, endpoint=False)
+ hkarray = np.array([hkfunc(kx) for kx in ks_dense])
+ _vals = np.linalg.eigvalsh(hkarray)
+ _gap = (utils.compute_gap(add_tb(h_0, mf_sol), fermi_energy=0, nk=nk_dense))
+ gap.append(_gap)
+ vals.append(_vals)
+ return np.asarray(gap, dtype=float), np.asarray(vals)
+
+import xarray as xr
+
+ds = xr.Dataset(
+ data_vars=dict(vals=(["Us", "ks", "n"], vals), gap=(["Us"], gap)),
+ coords=dict(
+ Us=Us,
+ ks=np.linspace(0, 2 * np.pi, nk_dense),
+ n=np.arange(vals.shape[-1])
+ ),
+)
+
+# Interaction strengths
+Us = np.linspace(0.5, 10, 20, endpoint=True)
+nk, nk_dense = 40, 100
+gap, vals = compute_phase_diagram(Us=Us, nk=nk, nk_dense=nk_dense)
+
+ds.vals.plot.scatter(x="ks", hue="Us", ec=None, s=5)
+plt.axhline(0, ls="--", 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 Hartree-Fock dispersion should follow (see [these notes](https://www.cond-mat.de/events/correl11/manuscript/Lechermann.pdf))
+$$
+\epsilon_{HF}^{\sigma}(\mathbf{k}) = \epsilon(\mathbf{k}) + U \left(\frac{n}{2} + \sigma m\right)
+$$
+where $m=(\langle n_{i\uparrow} \rangle - \langle n_{i\downarrow} \rangle) / 2$ is the magnetization per atom and $n = \sum_i \langle n_i \rangle$ is the total number of atoms per cell. Thus, for the antiferromagnetic groundstate, $m=1/2$ and $n=2$. The gap thus should be $\Delta=U$. And we can confirm it indeed follows the expected trend.
--
GitLab