# Strained graphene superlattice#

We showcase the interface between `meanfi`

and `Kwant`

with a strained graphene supercell. In this tutorial, we qualitatively reproduce the results from [1].

We first create the atomistic model in `Kwant`

. The complete source code of this example can be found in `strained_graphene_kwant.py`

. To reduce the computational cost, we perform the calculations with a \(16 \times 16\) supercell whereas in 1 the calculations were performed with a \(25 \times 25\) supercell. Thus, the agreement throughout the tutorial is only qualitative.

## Show code cell source

```
import kwant
import matplotlib.pyplot as plt
import meanfi
import numpy as np
from meanfi.kwant_helper import utils
from scripts.pauli import s0, sx, sy, sz
from scripts.strained_graphene_kwant import create_system, plot_bands
sigmas = [sx, sy, sz]
```

We verify the band structure of the Kwant model along a high-symmetry k-path.

```
h0_builder, lat, k_path = create_system(n=16)
```

## Show code cell source

```
fsyst = kwant.wraparound.wraparound(h0_builder).finalized()
params = {"xi": 7}
def hk(k):
return fsyst.hamiltonian_submatrix(
params={**params, **dict(k_x=k[0], k_y=k[1])}, sparse=False
)
plot_bands(hk, k_path)
```

We now use the Kwant model to create the interacting Hamiltonian. Following [1], we consider only onsite interactions.

```
def func_hop(site1, site2):
return 0 * np.ones((2, 2))
def func_onsite(site, U):
return U * np.ones((2, 2))
int_builder = utils.build_interacting_syst(
h0_builder,
lat,
func_onsite,
func_hop,
max_neighbor=0,
)
```

After we have created the interacting system we can use MeanFi again for getting the solution. We turn both the non-interacting and interacting systems into tight binding dictionaries using the kwant utils. Then we combine them into a mean-field model.

```
from meanfi.kwant_helper import utils as utils
h0, data = utils.builder_to_tb(h0_builder, params={"xi": 7}, return_data=True)
params_int = dict(U=0.8)
ndof = [*h0.values()][0].shape[0]
filling = ndof // 2
h_int = utils.builder_to_tb(int_builder, params_int)
mf_model = meanfi.Model(h0, h_int, filling=filling)
```

Now getting the solution by providing a guess and the mean-field model to the solver. To accelerate the convergence, we use an antiferromagnetic guess.

```
def func_hop(site1, site2):
return 0 * np.ones((2, 2))
def func_onsite(site):
if site.family == lat.sublattices[0]:
return sz
else:
return -sz
guess_builder = utils.build_interacting_syst(
h0_builder,
lat,
func_onsite,
func_hop,
max_neighbor=0,
)
guess = utils.builder_to_tb(guess_builder)
```

Due to the large supercell, lwe use a coarse k-point grid. Furthermore, to speed up the calculation we provide additional `optimizer_kwargs`

.

```
mf_sol = meanfi.solver(
mf_model,
guess,
nk=2,
optimizer_kwargs={
"M": 10,
"f_tol": 1e-4,
"maxiter": 100,
"line_search": "armijo",
},
)
```

We now verify that the mean-field solution results in a gapped phase.

```
mf_ham = meanfi.add_tb(h0, mf_sol)
hk_mf = meanfi.tb_to_kfunc(mf_ham)
plot_bands(hk_mf, k_path)
```

Now we turn the mean-field corrections into a `kwant.Builder`

such that we can visualize observables using Kwant’s functionalities. We provide the mean-field solution `mf_sol`

as well as the sites and periods of the bulk system to the `tb_to_builder`

function.

```
mf_sol_builder = utils.tb_to_builder(mf_sol, data["sites"], data["periods"])
```

We now plot the magnetization of the system. First, we define the magnetization direction. We do extracting the exchange field at it site and finding the direction where all the spins are pointing to.

```
exchange_field = np.array([value[1] for value in list(mf_sol_builder.site_value_pairs())])
magnetization_direction = np.zeros(3)
for i, sigma in enumerate(sigmas):
for site_field in exchange_field:
magnetization_direction[i] += np.abs(np.trace(sigma @ site_field))
magnetization_direction /= np.linalg.norm(magnetization_direction)
reference_magnetization = np.transpose(sigmas, (1, 2, 0)) @ magnetization_direction
```

Now we plot the magnetization of the solution on the sites of the system. We do this by creating functions which calculate the magnetization and its magnitude with respect to the reference magnetization direction.

```
def magnetisation(site):
matrix = mf_sol_builder.H[site][1]
return np.trace(reference_magnetization @ matrix).real
def abs_magnetisation(site):
matrix = mf_sol_builder.H[site][1]
projected_magnetization = []
for sigma in sigmas[1:]:
projected_magnetization.append(np.trace(sigma @ matrix))
return np.sqrt(np.sum(np.array(projected_magnetization) ** 2).real)
def systemPlotter(syst, onsite, ax, cmap):
"""
Plots the system with the onsite potential given by the function onsite.
"""
sites = [*syst.sites()]
density = [onsite(site) for site in sites]
def size(site):
return 0.3 * np.abs(onsite(site)) / np.max(np.abs(density))
kwant.plot(
syst, site_color=onsite, ax=ax, cmap=cmap, site_size=size, show=False, unit=1
)
return np.min(density), np.max(density)
```

## Show code cell source

```
import matplotlib
from mpl_toolkits.axes_grid1 import make_axes_locatable
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
titles = ["Magnetisation", "Magnetisation magnitude"]
cmaps = ["coolwarm", "viridis"]
onsites = [magnetisation, abs_magnetisation]
for i, ax in enumerate(axs):
vmin, vmax = systemPlotter(mf_sol_builder, onsites[i], ax=ax, cmap=cmaps[i])
ax.axis("off")
ax.set_title(titles[i])
ax.set_aspect("equal")
fig.colorbar(ax.collections[0], ax=ax, shrink=0.5)
fig.show()
```

```
/home/docs/checkouts/readthedocs.org/user_builds/meanfi/conda/latest/lib/python3.11/site-packages/kwant/_plotter.py:77: RuntimeWarning: plotly is not available, if other engines are unavailable, only iterator-providing functions will work
warnings.warn("plotly is not available, if other engines are unavailable,"
```