Examples

Everything related to constructing examples that motivate and teach how to use our package. This includes:

• coming up with new interesting physical systems
• cleaning up the existing ones

and likely more.

docs/source/graphene_example.md

@@ -15,7 +15,7 @@ kernelspec:
 
 ## The physics
 
-This tutorial serves as a simple example of using the meanfield algorithm in two dimensions in combination with using Kwant. We will consider a simple tight-binding model of graphene with a Hubbard interaction. The graphene system is first created using Kwant. For the basics of creating graphene with Kwant we refer to [this](https://kwant-project.org/doc/1/tutorial/graphene) tutorial.
+This tutorial serves as a simple example of using the mean-field algorithm in two dimensions in combination with using Kwant. We will consider a simple tight-binding model of graphene with a Hubbard interaction. The graphene system is first created using Kwant. For the basics of creating graphene with Kwant we refer to [this](https://kwant-project.org/doc/1/tutorial/graphene) tutorial.
 
 We begin with the basic imports
 
@@ -24,13 +24,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 import numpy as np
 import matplotlib.pyplot as plt
-import pymf.model as model
-import pymf.solvers as solvers
-import pymf.mf as mf
-import pymf.tb as tb
-import pymf.observables as observables
-import pymf.kwant_helper.kwant_examples as kwant_examples
-import pymf.kwant_helper.utils as kwant_utils
+import pymf
 ```
 
 ##  Preparing the model
@@ -38,6 +32,9 @@ import pymf.kwant_helper.utils as kwant_utils
 We first translate this model from a Kwant system to a tight-binding dictionary. In the tight-binding dictionary the keys denote the hoppings while the values are the hopping amplitudes.
 
 ```{code-cell} ipython3
+import pymf.kwant_helper.kwant_examples as kwant_examples
+import pymf.kwant_helper.utils as kwant_utils
+
 # Create translationally-invariant `kwant.Builder`
 graphene_builder, int_builder = kwant_examples.graphene_extended_hubbard()
 h_0 = kwant_utils.builder_to_tb(graphene_builder)
@@ -54,18 +51,18 @@ U=1
 V=0.1
 params = dict(U=U, V=V)
 h_int = kwant_utils.builder_to_tb(int_builder, params)
-_model = model.Model(h_0, h_int, filling=2)
+_model = pymf.Model(h_0, h_int, filling=2)
 ```
 
-To start the meanfield calculation we also need a starting guess. We will use our random guess generator for this. It creates a random Hermitian hopping dictionary based on the hopping keys provided and the number of degrees of freedom specified. As we don't expect the mean-field solution to contain terms more than the hoppings from the interacting part, we can use the hopping keys from the interacting part. We will use the same number of degrees as freedom as both the non-interacting and interacting part, so that they match.
+To start the mean-field calculation we also need a starting guess. We will use our random guess generator for this. It creates a random Hermitian hopping dictionary based on the hopping keys provided and the number of degrees of freedom specified. As we don't expect the mean-field solution to contain terms more than the hoppings from the interacting part, we can use the hopping keys from the interacting part. We will use the same number of degrees as freedom as both the non-interacting and interacting part, so that they match.
 
 ```{code-cell} ipython3
-guess = tb.utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
-mf_sol = solvers.solver(_model, guess, nk=18, optimizer_kwargs={'M':0})
-full_sol = tb.tb.add_tb(h_0, mf_sol)
+guess = pymf.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
+mf_sol = pymf.solver(_model, guess, nk=18, optimizer_kwargs={'M':0})
+full_sol = pymf.add_tb(h_0, mf_sol)
 ```
 
-After we have defined the guess, we feed it together with the model into the meanfield solver. The meanfield solver will return a hopping dictionary with the meanfield approximation. We can then add this solution to the non-interacting part to get the full solution. In order to get the solution, we specified the number of k-points to be used in the calculation. This refers to the k-grid used in the Brillouin zone for the density matrix.
+After we have defined the guess, we feed it together with the model into the mean-field solver. The mean-field solver will return a hopping dictionary with the mean-field approximation. We can then add this solution to the non-interacting part to get the full solution. In order to get the solution, we specified the number of k-points to be used in the calculation. This refers to the k-grid used in the Brillouin zone for the density matrix.
 
 ## Creating a phase diagram of the gap
 
@@ -73,7 +70,7 @@ We can now create a phase diagram of the gap of the interacting solution. In ord
 
 ```{code-cell} ipython3
 def compute_gap(h, fermi_energy=0, nk=100):
-    kham = tb.transforms.tb_to_khamvector(h, nk, ks=None)
+    kham = pymf.tb_to_khamvector(h, nk, ks=None)
     vals = np.linalg.eigvalsh(kham)
 
     emax = np.max(vals[vals <= fermi_energy])
@@ -87,10 +84,10 @@ Now that we can calculate the gap, we create a phase diagram of the gap as a fun
 def gap_and_mf_sol(U, V, int_builder, h_0):
     params = dict(U=U, V=V)
     h_int = kwant_utils.builder_to_tb(int_builder, params)
-    _model = model.Model(h_0, h_int, filling=2)
-    guess = tb.utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
-    mf_sol = solvers.solver(_model, guess, nk=18, optimizer_kwargs={'M':0})
-    gap = compute_gap(tb.tb.add_tb(h_0, mf_sol), fermi_energy=0, nk=300)
+    _model = pymf.Model(h_0, h_int, filling=2)
+    guess = pymf.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
+    mf_sol = pymf.solver(_model, guess, nk=18, optimizer_kwargs={'M':0})
+    gap = compute_gap(pymf.add_tb(h_0, mf_sol), fermi_energy=0, nk=300)
     return gap, mf_sol
 ```
 
@@ -112,8 +109,8 @@ We chose to initialize a new guess for each $U$ value, but not for each $V$ valu
 We can now compute the phase diagram and then plot it
 
 ```{code-cell} ipython3
-Us = np.linspace(0, 4, 5)
-Vs = np.linspace(0, 1.5, 5)
+Us = np.linspace(0, 4, 10)
+Vs = np.linspace(0, 1.5, 10)
 gap, mf_sols = compute_phase_diagram(Us, Vs, int_builder, h_0)
 plt.imshow(gap.T, extent=(Us[0], Us[-1], Vs[0], Vs[-1]), origin='lower', aspect='auto')
 plt.colorbar()
@@ -144,8 +141,8 @@ We choose a point in the phase diagram where we expect there to be a CDW phase a
 ```{code-cell} ipython3
 cdw_list = []
 for mf_sol in mf_sols.flatten():
-    rho, _ = mf.construct_density_matrix(tb.tb.add_tb(h_0, mf_sol), filling=2, nk=40)
-    expectation_value = observables.expectation_value(rho, cdw_order_parameter)
+    rho, _ = pymf.construct_density_matrix(pymf.add_tb(h_0, mf_sol), filling=2, nk=40)
+    expectation_value = pymf.expectation_value(rho, cdw_order_parameter)
     cdw_list.append(expectation_value)
 ```
 
@@ -181,10 +178,10 @@ Then, similar to what we did in the CDW phase, we calculate the expectation valu
 ```{code-cell} ipython3
 sdw_list = []
 for mf_sol in mf_sols.flatten():
-    rho, _ = mf.construct_density_matrix(tb.tb.add_tb(h_0, mf_sol), filling=2, nk=40)
+    rho, _ = pymf.construct_density_matrix(pymf.add_tb(h_0, mf_sol), filling=2, nk=40)
     expectation_values = []
     for order_parameter in order_parameter_list:
-        expectation_value = observables.expectation_value(rho, order_parameter)
+        expectation_value = pymf.expectation_value(rho, order_parameter)
         expectation_values.append(expectation_value)
 
     sdw_list.append(np.sum(np.array(expectation_values)**2))
@@ -199,3 +196,17 @@ plt.ylabel('U')
 plt.title('Spin Density Wave Order Parameter')
 plt.show()
 ```
+
+## Full phase diagram
+
+Finally, we can combine the gap, CDW and SDW phase diagrams into one plot. We naively do this by plotting the order parameter of CDW minus the order parameter of SDW. Furthermore, we normalize the gap such that it is between $0$ and $1$ and can thus be used for the transparency.
+
+```{code-cell} ipython3
+import matplotlib.ticker as mticker
+normalized_gap = gap/np.max(gap)
+plt.imshow(np.abs(cdw_list.T.real)-np.abs(sdw_list.T.real), extent=(Us[0], Us[-1], Vs[0], Vs[-1]), origin='lower', aspect='auto', cmap="coolwarm", alpha=normalized_gap.T)
+plt.colorbar(ticks=[-1.75, 0, 1.75], format=mticker.FixedFormatter(['SDW', '0', 'CDW']), label='Order parameter', extend='both')
+plt.xlabel('V')
+plt.ylabel('U')
+plt.show()
+```