Skip to content
Snippets Groups Projects

Examples

Merged Examples
examples into main
Merged Kostas Vilkelis requested to merge examples into main
Viewing commit 4a66cc49
Prev Next
Show latest version
1 file
+ 27
− 67
Compare changes
  • Side-by-side
  • Inline
docs/source/hubbard_1d.md
+ 27
− 67
@@ -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$.