# 7. Defects and strain¶

This section will introduce `@site_state_modifier`

and
`@site_position_modifier`

which can be used to add defects and
strain to the model. These modifiers are applied to the structure of the system before the
Hamiltonian matrix is created.

`Download this page as a Jupyter notebook`

## 7.1. Vacancies¶

A `@site_state_modifier`

can be used to create vacancies in a crystal
lattice. The definition is very similar to the onsite and hopping modifiers explained in the
previous section.

```
def vacancy(position, radius):
@pb.site_state_modifier
def modifier(state, x, y):
x0, y0 = position
state[(x-x0)**2 + (y-y0)**2 < radius**2] = False
return state
return modifier
```

The `state`

argument indicates the current boolean state of a lattice site. Only valid sites
(`True`

state) will be included in the final Hamiltonian matrix. Therefore, setting the state of
sites within a small radius to `False`

will exclude them from the final system. The `x`

and `y`

arguments are lattice site positions. As with the other modifiers, the arguments are optional
(`z`

is not needed for this example) but the full signature of the site state modifier can be
found on its `API reference page`

.

This is actually very similar to the way a `FreeformShape`

works. In fact, it is possible
to create defects by defining them directly in the shape. However, such an approach would not be
very flexible since we would need to create an entire new shape in order to change either the
vacancy type or the shape itself. By defining the vacancy as a modifier, we can simply compose
it with any existing shapes:

```
from pybinding.repository import graphene
model = pb.Model(
graphene.monolayer(),
pb.rectangle(2),
vacancy(position=[0, 0], radius=0.1)
)
model.plot()
```

The resulting 2-atom vacancy is visible in the center of the system. The two vacant sites are
completely removed from the final Hamiltonian matrix. If we were to inspect the number of rows
and columns by looking up `model.hamiltonian.shape`

, we would see that the size of the matrix is
reduced by 2.

Any number of modifiers can be included in the model and they will compose as expected. We can take advantage of this and create four different vacancies, with 1 to 4 missing atoms:

```
model = pb.Model(
graphene.monolayer(),
pb.rectangle(2),
vacancy(position=[-0.50, 0.50], radius=0.1),
vacancy(position=[ 0.50, 0.45], radius=0.15),
vacancy(position=[-0.45, -0.45], radius=0.15),
vacancy(position=[ 0.50, -0.50], radius=0.2),
)
model.plot()
```

## 7.2. Layer defect¶

The site state modifier also has access to sublattice information. This can be used, for example,
with bilayer graphene to remove a single layer in a specific area. We’ll use the bilayer lattice
that’s included in the Material Repository. The `graphene.bilayer()`

lattice is laid out so that sublattices A1 and B1 belong to the top layer, while A2 and B2 are on
the bottom.

```
def scrape_top_layer(position, radius):
"""Remove the top layer of graphene in the area specified by position and radius"""
@pb.site_state_modifier
def modifier(state, x, y, sub_id):
x0, y0 = position
is_within_radius = (x-x0)**2 + (y-y0)**2 < radius**2
is_top_layer = np.logical_or(sub_id == 'A1', sub_id == 'B1')
final_condition = np.logical_and(is_within_radius, is_top_layer)
state[final_condition] = False
return state
return modifier
model = pb.Model(
graphene.bilayer(),
pb.rectangle(2),
scrape_top_layer(position=[0, 0], radius=0.5)
)
model.plot()
```

The central monolayer area is nicely visible in the figure. We can actually create the same
structure in a different way: by considering the `z`

position of the lattice site to distinguish
the layers. An alternative modifier definition is given below. It would generate the same figure.
Which method is more convenient is up to the user.

```
def scrape_top_layer_alt(position, radius):
"""Alternative definition of `scrape_top_layer`"""
@pb.site_state_modifier
def modifier(state, x, y, z):
x0, y0 = position
is_within_radius = (x-x0)**2 + (y-y0)**2 < radius**2
is_top_layer = (z == 0)
final_condition = np.logical_and(is_within_radius, is_top_layer)
state[final_condition] = False
return state
return modifier
```

Note

As with the onsite and hopping modifiers, all the arguments are given as numpy arrays.
Therefore, we must use the array-specific `np.logical_or()`

/
`np.logical_and()`

functions instead of the plain `or`

/`and`

keywords.

## 7.3. Strain¶

A `@site_position_modifier`

can be used to model the lattice site
displacement caused by strain. Let’s start with a simple triangular system:

```
from math import pi
model = pb.Model(
graphene.monolayer(),
pb.regular_polygon(num_sides=3, radius=2, angle=pi),
)
model.plot()
```

We’re going to apply strain in three directions, as if we are pulling outward on the vertices of
the triangle. The displacement function for this kind of strain is given below. The `c`

parameter
lets us control the intensity of the strain.

```
def triaxial_displacement(c):
@pb.site_position_modifier
def displacement(x, y, z):
ux = 2*c * x*y
uy = c * (x**2 - y**2)
return x + ux, y + uy, z
return displacement
```

The modifier function takes the `x`

, `y`

, `z`

coordinates as arguments. The displacement `ux`

, `uy`

is computed and the modified coordinates are returned. The `z`

argument is returned unchanged but
we still need it here because the modifier is expected to always return all three.

```
model = pb.Model(
graphene.monolayer(),
pb.regular_polygon(num_sides=3, radius=2, angle=pi),
triaxial_displacement(c=0.15)
)
model.plot()
```

As seen in the figure, the displacement has been applied to the lattice sites and the new position data is saved in the system. However, the hopping energies have not been modified yet. Every hopping element of the Hamiltonian matrix is equal to the hopping energy of pristine graphene:

```
>>> np.all(model.hamiltonian.data == -2.8)
True
```

We now need to use the new position data to modify the hopping energy according to the relation
\(t = t_0 e^{-\beta (\frac{d}{a_{cc}} - 1)}\), where \(t_0\) is the original unstrained
hopping energy, \(\beta\) controls the strength of the strain-induced hopping modulation,
\(d\) is the strained distance between two atoms and \(a_{cc}\) is the unstrained
carbon-carbon distance. This can be implemented using a
`@hopping_energy_modifier`

:

```
@pb.hopping_energy_modifier
def strained_hopping(energy, x1, y1, z1, x2, y2, z2):
d = np.sqrt((x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2)
beta = 3.37
w = d / graphene.a_cc - 1
return energy * np.exp(-beta*w)
```

The structural modifiers (site state and position) are always automatically applied to the model
before energy modifiers (onsite and hopping). Thus, our `strain_hopping`

modifier will get the new
displaced coordinates as its arguments, from which it will calculate the strained hopping energy.

```
model = pb.Model(
graphene.monolayer(),
pb.regular_polygon(num_sides=3, radius=2, angle=pi),
triaxial_displacement(c=0.15),
strained_hopping
)
```

Including the hopping modifier along with the displacement will yield position dependent hopping energy, thus the elements of the Hamiltonian will no longer be all equal:

```
>>> np.all(model.hamiltonian.data == -2.8)
False
```

However, it isn’t convenient to keep track of the displacement and strained hoppings separately. Instead, we can package them together in one function which is going to return both modifiers:

```
def triaxial_strain(c, beta=3.37):
"""Produce both the displacement and hopping energy modifier"""
@pb.site_position_modifier
def displacement(x, y, z):
ux = 2*c * x*y
uy = c * (x**2 - y**2)
return x + ux, y + uy, z
@pb.hopping_energy_modifier
def strained_hopping(energy, x1, y1, z1, x2, y2, z2):
l = np.sqrt((x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2)
w = l / graphene.a_cc - 1
return energy * np.exp(-beta*w)
return displacement, strained_hopping
```

The `triaxial_strain`

function now has everything we need. We’ll apply it to a slightly larger
system so that we can clearly calculate the local density of states (LDOS). For more information
about this computation method see the Kernel polynomial method section. Right now, it’s enough to know that
we will calculate the LDOS at the center of the strained system, separately for sublattices
A and B.

```
model = pb.Model(
graphene.monolayer(),
pb.regular_polygon(num_sides=3, radius=40, angle=pi),
triaxial_strain(c=0.0025)
)
kpm = pb.kpm(model)
for sub_name in ['A', 'B']:
ldos = kpm.calc_ldos(energy=np.linspace(-1, 1, 500), broadening=0.03,
position=[0, 0], sublattice=sub_name)
ldos.plot(label=sub_name, ls="--" if sub_name == "B" else "-")
pb.pltutils.legend()
```

Strain in graphene has an effect similar to a magnetic field. That’s why we see Landau-level-like features in the LDOS. Note that the zero-energy peak has double intensity on one sublattice but zero on the other: this is a unique feature of the strain-induced pseudo-magnetic field.

## 7.5. Example¶

```
"""Strain a triangular system by pulling on its vertices"""
import pybinding as pb
import numpy as np
import matplotlib.pyplot as plt
from pybinding.repository import graphene
from math import pi
pb.pltutils.use_style()
def triaxial_strain(c):
"""Strain-induced displacement and hopping energy modification"""
@pb.site_position_modifier
def displacement(x, y, z):
ux = 2*c * x*y
uy = c * (x**2 - y**2)
return x + ux, y + uy, z
@pb.hopping_energy_modifier
def strained_hopping(energy, x1, y1, z1, x2, y2, z2):
l = np.sqrt((x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2)
w = l / graphene.a_cc - 1
return energy * np.exp(-3.37 * w)
return displacement, strained_hopping
model = pb.Model(
graphene.monolayer(),
pb.regular_polygon(num_sides=3, radius=2, angle=pi),
triaxial_strain(c=0.1)
)
model.plot()
plt.show()
```

```
plt.figure(figsize=(7, 2.5))
grid = plt.GridSpec(nrows=1, ncols=2)
for block, energy in zip(grid, [0, 0.25]):
plt.subplot(block)
plt.title("E = {} eV".format(energy))
solver = pb.solver.arpack(model, k=30, sigma=energy)
ldos_map = solver.calc_spatial_ldos(energy=energy, broadening=0.03)
ldos_map.plot()
pb.pltutils.colorbar(label="LDOS")
plt.show()
```