From 4a66cc49fd56a983f763624d66d6423520a58680 Mon Sep 17 00:00:00 2001
From: Johanna <johanna@zijderveld.de>
Date: Tue, 7 May 2024 11:59:51 +0200
Subject: [PATCH] fix imports and remove unlooked at part
---
docs/source/hubbard_1d.md | 94 +++++++++++----------------------------
1 file changed, 27 insertions(+), 67 deletions(-)
diff --git a/docs/source/hubbard_1d.md b/docs/source/hubbard_1d.md
index 90ebb17..48cd977 100644
--- a/docs/source/hubbard_1d.md
+++ b/docs/source/hubbard_1d.md
@@ -13,47 +13,49 @@ kernelspec:
# 1d Hubbard
+To simulate how the package can be used to simulate infinite systems, we show how to use it with a tight-binding model in 1 dimension.
+We exemplify this by computing the ground state of an infinite spinful chain with onsite interactions.
+First, the basic imports are done.
+
```{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:
+After this, we start by constructing the non-interacting Hamiltonian. As we expect the ground state to be an antiferromagnet, we build a two-atom cell. We name the two sublattices, $A$ and $B$. The Hamiltonian is then:
$$
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}$.
+We write down the spinful part by simply taking $H_0(k) \otimes \mathbb{1}$.
+To translat ethis Hamiltonian into a tight-binding model, we specify the hopping vectors together with the hopping amplitudes. We ensure that the Hamiltonian is hermitian by letting the hopping amplitudes from $A$ to $B$ be the complex conjugate of the hopping amplitudes from $B$ to $A$.
-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
+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()}
+```
+
+We verify this tight-binding model by plotting the band structure and observing the two bands due to the Brillouin zone folding. In order to do this we transform the tight-binding model into a Hamiltonian on a k-point grid using the `tb_to_khamvector` and then diagonalize it.
```{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)
+hamiltonians_0 = pymf.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.
+After confirming that the non-interacting part is correct, we can set up the interacting Hamiltonian. We define the interaction similarly as a tight-binding dictionary. To keep the physics simple, we let the interaction be onsite only, which gives us the following interaction matrix:
-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}
@@ -64,60 +66,18 @@ H_{int} =
\end{array}\right)~.
$$
+This is simply constructed by writing:
+
```{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()
+h_int = {(0,): U * np.kron(np.ones((2, 2)), np.eye(2)),}
+```
+
+In order to find a meanfield solution, we combine the non interacting Hamiltonian with the interaction Hamiltonian and the relevant filling into a `Model` object. We then generate a starting guess for the meanfield solution and solve the model using the `solver` function. It is important to note that the guess will influence the possible solutions which the `solver` can find in the meanfield procedure. The `generate_guess` function generates a random Hermitian tight-binding dictionary, with the keys provided as hopping vectors and the values of the size as specified.
+```{code-cell} ipython3
+full_model = pymf.Model(h_0, h_int, filling)
+guess = pymf.generate_guess(frozenset(h_int), ndof=4)
+mf_sol = pymf.solver(full_model, guess, nk=nk)
```
-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.
+The `solver` function returns only the meanfield correction to the non-interacting Hamiltonian. To get the full Hamiltonian, we add the meanfield correction to the non-interacting Hamiltonian. To take a look at whether the result is correct, we first do the meanfield computation for a wider range of $U$ values and then plot the gap as a function of $U$.
--
GitLab