Skip to content
Snippets Groups Projects

Examples

Merged Examples
examples into main
Merged Kostas Vilkelis requested to merge examples into main
Compare
and
Show latest version
4 files
+ 188
− 61
Compare changes
  • Side-by-side
  • Inline
Files
4
docs/source/graphene_example.md
+ 81
− 60
@@ -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])
@@ -83,39 +80,37 @@ def compute_gap(h, fermi_energy=0, nk=100):
Now that we can calculate the gap, we create a phase diagram of the gap as a function of the Hubbard interaction strength $U$ and the nearest neighbor interaction strength $V$. We vary the onsite Hubbard interactio $U$ strength from $0$ to $2$ and the nearest neighbor interaction strength $V$ from $0$ to $1.5$.
```{code-cell} ipython3
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 = 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
```
```{code-cell} ipython3
def compute_phase_diagram(Us, Vs, int_builder, h_0):
gap = []
mf_sols = []
for U in Us:
gap_U = []
mf_sols_U = []
for V in Vs:
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)
converged=False
while not converged:
guess = tb.utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
try:
mf_sol = solvers.solver(_model, guess, nk=18, optimizer_kwargs={'M':0})
converged=True
except:
converged=False
gap_U.append(compute_gap(tb.tb.add_tb(h_0, mf_sol), fermi_energy=0, nk=300))
mf_sols_U.append(mf_sol)
guess = None
gap.append(gap_U)
mf_sols.append(mf_sols_U)
return np.asarray(gap, dtype=float), np.asarray(mf_sols)
gaps = []
mf_sols = []
for U in Us:
for V in Vs:
gap, mf_sol = gap_and_mf_sol(U, V, int_builder, h_0)
gaps.append(gap)
mf_sols.append(mf_sol)
gaps = np.asarray(gaps, dtype=float).reshape((len(Us), len(Vs)))
mf_sols = np.asarray(mf_sols).reshape((len(Us), len(Vs)))
return gaps, mf_sols
```
We chose to initialize a new guess for each $U$ value, but not for each $V$ value. Instead, for consecutive $V$ values we use the previous mean-field solution as a guess. We do this because the mean-field solution is expected to be smooth in the interaction strength and therefore by using an inspired guess we can speed up the calculation.
We can now compute the phase diagram and plot it.
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()
@@ -141,37 +136,63 @@ cdw_order_parameter = {}
cdw_order_parameter[(0,0)] = np.kron(sz, np.eye(2))
```
We choose a point in the phase diagram where we expect there to be a CDW phase and calculate the expectation value with the CDW order parameter. In order to do this we first construct the density matrix from the mean field solution.
We choose a point in the phase diagram where we expect there to be a CDW phase and calculate the expectation value with the CDW order parameter. In order to do this we first construct the density matrix from the mean field solution. We perform this calculation over the complete phase diagram where we calculated the gap earlier:
```{code-cell} ipython3
params = dict(U=0, V=2)
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})
full_sol = tb.tb.add_tb(h_0, mf_sol)
rho, _ = mf.construct_density_matrix(full_sol, filling=2, nk=40)
expectation_value = observables.expectation_value(rho, cdw_order_parameter)
print(expectation_value)
cdw_list = []
for mf_sol in mf_sols.flatten():
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)
```
We can also perform the same calculation over the complete phase diagram where we calculated the gap earlier:
```{code-cell} ipython3
cdw_list = np.asarray(cdw_list).reshape(mf_sols.shape)
plt.imshow(np.abs(cdw_list.T.real), extent=(Us[0], Us[-1], Vs[0], Vs[-1]), origin='lower', aspect='auto')
plt.colorbar()
plt.xlabel('V')
plt.ylabel('U')
plt.title('Charge Density Wave Order Parameter')
plt.show()
```
### Spin density wave
To check the other phase we expect in the graphene phase diagram, we construct a spin density wave order parameter. In our chosen graphene system the spin density wave has $SU(2)$ symmetry. This means that we need to sum over the pauli matrices when constructing this order parameter. We can construct the order parameter as:
```{code-cell} ipython3
expectation_value_list = []
sx = np.array([[0, 1], [1, 0]])
sy = np.array([[0, -1j], [1j, 0]])
s_list = [sx, sy, sz]
order_parameter_list = []
for s in s_list:
order_parameter = {}
order_parameter[(0,0)] = np.kron(sz, s)
order_parameter_list.append(order_parameter)
```
Then, similar to what we did in the CDW phase, we calculate the expectation value of the order parameter with the density matrix of the mean field solution over the complete phase diagram. The main subtlety here is that we need to sum over expectation value of all SDW direction order parameters defined in order to get the total spin density wave order parameter.
```{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)
expectation_value = observables.expectation_value(rho, cdw_order_parameter)
expectation_value_list.append(expectation_value)
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 = pymf.expectation_value(rho, order_parameter)
expectation_values.append(expectation_value)
sdw_list.append(np.sum(np.array(expectation_values)**2))
```
```{code-cell} ipython3
expectation_value_list = np.asarray(expectation_value_list).reshape(mf_sols.shape)
plt.imshow(np.abs(expectation_value_list.T.real), extent=(Us[0], Us[-1], Vs[0], Vs[-1]), origin='lower', aspect='auto')
sdw_list = np.asarray(sdw_list).reshape(mf_sols.shape)
plt.imshow(np.abs(sdw_list.T.real), extent=(Us[0], Us[-1], Vs[0], Vs[-1]), origin='lower', aspect='auto')
plt.colorbar()
plt.xlabel('V')
plt.ylabel('U')
plt.title('Charge Density Wave Order Parameter')
plt.title('Spin Density Wave Order Parameter')
plt.show()
```