Merge branch 'examples' into 'main'

Examples

See merge request qt/kwant-scf!7

docs/source/algorithm.md → 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-consistency loop

docs/source/mf_notes.md → 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

docs/source/index.md

@@ -18,6 +18,8 @@ kernelspec:
 :maxdepth: 1
 :caption: Tutorials
 
+tutorial/hubbard_1d.md
+tutorial/graphene_example.md
 ```
 
 ```{toctree}
@@ -25,19 +27,23 @@ kernelspec:
 :maxdepth: 1
 :caption: Documentation
 
-mf_notes.md
-algorithm.md
+documentation/mf_notes.md
+documentation/algorithm.md
 documentation/pymf.md
 ```
 
 ## What is pymf?
 
+
 ## Why pymf?
 
+
 ## How does pymf work?
 
 ## What does pymf not do yet?
 
+* Superconductivity
+
 ## Installation
 
 ```bash

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 `pymf` to solve a simple 1D Hubbard model with onsite interactions.
+In this tutorial, we will apply `pymf` 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 pymf
+from pymf.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}`~pymf.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}`~pymf.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}`~pymf.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}`~pymf.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 = pymf.Model(h_0, h_int, filling=2)
+int_keys = frozenset(h_int)
+ndof = len(list(h_0.values())[0])
+guess = pymf.generate_guess(int_keys, ndof)
+mf_sol = pymf.solver(model, guess, nk=18)
+h_full = pymf.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}`~pymf.mf.construct_density_matrix` function.
+We then feed it into {autolink}`~pymf.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, _ = pymf.construct_density_matrix(h_full, filling=filling, nk=40)
+rho_0, _ = pymf.construct_density_matrix(h_0, filling=filling, nk=40)
+
+cdw_order_parameter = pymf.expectation_value(rho, cdw_operator)
+cdw_order_parameter_0 = pymf.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 = pymf.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 = pymf.Model(h_0, h_int, filling=filling)
+        guess = pymf.generate_guess(int_keys, ndof)
+        mf_sol = pymf.solver(model, guess, nk=18)
+        mf_sols.append(mf_sol)
+
+        gap = compute_gap(pymf.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, _ = pymf.construct_density_matrix(pymf.add_tb(h_0, mf_sol), filling=filling, nk=40)
+
+    # Compute CDW order parameter
+    cdw_list.append(np.abs(pymf.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(pymf.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 pymf
+```
+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}`~pymf.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 = pymf.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}`~pymf.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 = pymf.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}`~pymf.tb.utils.generate_guess` 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}`~pymf.tb.utils.generate_guess`.
+Finally, to solve the model, we use the {autolink}`~pymf.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 = 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 {autolink}`~pymf.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 = pymf.add_tb(h_0, mf_sol)
+
+hamiltonians = pymf.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 = pymf.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
+    full_model = pymf.Model(h_0, h_int, filling)
+    mf_sol = pymf.solver(full_model, guess, nk=nk)
+    return pymf.add_tb(h_0, mf_sol)
+
+
+def compute_gap_and_vals(full_sol, nk_dense, fermi_energy=0):
+    h_kgrid = pymf.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), vals
+
+
+def compute_phase_diagram(
+    Us,
+    nk,
+    nk_dense,
+):
+    gaps = []
+    vals = []
+    for U in Us:
+        full_sol = compute_sol(U, h_0, nk)
+        gap, _vals = compute_gap_and_vals(full_sol, nk_dense)
+        gaps.append(gap)
+        vals.append(_vals)
+
+    return np.asarray(gaps, dtype=float), np.asarray(vals)
+
+Us = np.linspace(0, 4, 40, endpoint=True)
+gap, vals = compute_phase_diagram(Us=Us, nk=20, nk_dense=100)
+
+plt.plot(Us, gap, 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