{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import pybinding as pb\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"pb.pltutils.use_style()\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fields and effects\n",
"\n",
"This section will introduce [`@onsite_energy_modifier`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.onsite_energy_modifier.html#pybinding.onsite_energy_modifier) and [`@hopping_energy_modifier`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.hopping_energy_modifier.html#pybinding.hopping_energy_modifier) which can be used to add various fields to the model. These functions can apply user-defined modifications to the Hamiltonian matrix which is why we shall refer to them as modifier functions.\n",
"\n",
"[Download this page as a Jupyter notebook](http://docs.pybinding.site/page/_downloads/fields.ipynb)\n",
"\n",
"## Electric potential\n",
"\n",
"We can define a simple potential function like the following:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"@pb.onsite_energy_modifier\n",
"def potential(x, y):\n",
" return np.sin(x)**2 + np.cos(y)**2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here `potential` is just a regular Python function, but we attached a pretty `@` decorator to it. The [`@onsite_energy_modifier`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.onsite_energy_modifier.html#pybinding.onsite_energy_modifier) decorator gives an ordinary function a few extra properties which we'll talk about later. For now, just keep in mind that this is required to mark a function as a modifier for use with pybinding models. The `x` and `y` arguments are lattice site positions and the return value is the desired potential. Note the use of `np.sin` instead of `math.sin`. The `x` and `y` coordinates are `numpy` arrays, not individual numbers. This is true for all modifier arguments in pybinding. When you write modifier functions, make sure to always use `numpy` operations which work with arrays, unlike regular `math`.\n",
"\n",
"> Modifier arguments are passed as arrays for performance. Working with individual numbers would require calling the potential function individually for each lattice site which would be extremely slow. Arrays are much faster.\n",
"\n",
"To use the potential function, just place it in a [`Model`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Model.html#pybinding.Model) parameter list:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from pybinding.repository import graphene\n",
"\n",
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" pb.rectangle(12),\n",
" potential\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To visualize the potential, there's the handy [`Model.onsite_map`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Model.html#pybinding.Model.onsite_map) property which is a [`StructureMap`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.StructureMap.html#pybinding.StructureMap) of the onsite energy of the Hamiltonian matrix."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model.onsite_map.plot_contourf()\n",
"pb.pltutils.colorbar(label=\"U (eV)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The figure shows a 2D colormap representation of our wavy potential in a square system. The [`StructureMap.plot_contourf()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.StructureMap.html#pybinding.StructureMap.plot_contourf) method we just called is implemented in terms of matplotlib's `contourf` function with some slight adjustments for convenience.\n",
"\n",
"To make the potential more flexible, it's a good idea to enclose it in an outer function, just like this:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def wavy(a, b):\n",
" @pb.onsite_energy_modifier\n",
" def potential(x, y):\n",
" return np.sin(a * x)**2 + np.cos(b * y)**2\n",
" return potential\n",
"\n",
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" pb.regular_polygon(num_sides=6, radius=8),\n",
" wavy(a=0.6, b=0.9)\n",
")\n",
"model.onsite_map.plot_contourf()\n",
"pb.pltutils.colorbar(label=\"U (eV)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that we are using a system with hexagonal shape this time (via [`regular_polygon()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.regular_polygon.html#pybinding.regular_polygon)). The potential is only plotted inside the area of the actual system.\n",
"\n",
"We can make one more improvement to our `wavy` function. We'll add an `energy` argument:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def wavy2(a, b):\n",
" @pb.onsite_energy_modifier\n",
" def potential(energy, x, y):\n",
" v = np.sin(a * x)**2 + np.cos(b * y)**2\n",
" return energy + v\n",
" return potential"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `energy` argument contains the existing onsite energy in the system before the new potential function is applied. By adding to the existing energy, instead of just setting it, we can compose multiple functions. For example, let's combine the improved `wavy2` with a linear potential."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def linear(k):\n",
" @pb.onsite_energy_modifier\n",
" def potential(energy, x):\n",
" return energy + k*x\n",
" return potential\n",
"\n",
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" pb.regular_polygon(num_sides=6, radius=8),\n",
" wavy2(a=0.6, b=0.9),\n",
" linear(k=0.2)\n",
")\n",
"model.onsite_map.plot_contourf()\n",
"pb.pltutils.colorbar(label=\"U (eV)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see a similar wavy pattern as before, but the magnitude increases linearly along the x-axis because of the contribution of the `linear` potential.\n",
"\n",
"## About the decorator\n",
"\n",
"Now that you have a general idea of how to add and compose electric potentials in a model, we should talk about the role of the [`@onsite_energy_modifier`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.onsite_energy_modifier.html#pybinding.onsite_energy_modifier). The full signature of a potential function looks like this:\n",
"\n",
"```python\n",
"@pb.onsite_energy_modifier\n",
"def potential(energy, x, y, z, sub_id):\n",
" return ... # some function of the arguments\n",
"```\n",
"\n",
"This function uses all of the possible arguments of an onsite energy modifier: `energy`, `x`, `y`, `z` and `sub_id`. We have already explained the first three. The `z` argument is, obviously, the z-axis coordinate of the lattice sites. The `sub_id` argument tells us which sublattice a site belongs to. Its usage will be explained below.\n",
"\n",
"As we have seen before, we don't actually need to define a function to take all the arguments. They are optional. The `@` decorator will recognize a function which takes any of these arguments and it will adapt it for use in a pybinding model. Previously, the `linear` function accepted only the `energy` and `x` arguments, but `wavy` also included the `y` argument. The order of arguments is not important, only their names are. Therefore, this is also a valid modifier:\n",
"\n",
"```python\n",
"@pb.onsite_energy_modifier\n",
"def potential(x, y, energy, sub_id):\n",
" return ... # some function\n",
"```\n",
"\n",
"But the argument names must be exact: a typo or an extra unknown argument will result in an error. The decorator checks this at definition time and decides if the given function is a valid modifier or not, so any errors will be caught early.\n",
"\n",
"## Opening a band gap\n",
"\n",
"The last thing to explain about [`@onsite_energy_modifier`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.onsite_energy_modifier.html#pybinding.onsite_energy_modifier) is the use of the `sub_id` argument. It tells us which sublattice a site belongs to. If you remember from early on in the tutorial, [in the process of specifying a lattice](http://docs.pybinding.site/page/tutorial/lattice.html), we gave each sublattice a unique name. This name can be used to filter out sites of a specific sublattice. For example, let's add mass to electrons in graphene:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def mass_term(delta):\n",
" \"\"\"Break sublattice symmetry with opposite A and B onsite energy\"\"\"\n",
" @pb.onsite_energy_modifier\n",
" def potential(energy, sub_id):\n",
" energy[sub_id == 'A'] += delta\n",
" energy[sub_id == 'B'] -= delta\n",
" return energy\n",
" return potential"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that we don't need `x`, `y` or `z` arguments because this will be applied everywhere evenly. The `mass_term` function will add an energy `delta` to all sites on sublattice `A` and subtract `delta` from all `B` sites. Note that we are indexing the `energy` array with a condition on the `sub_id` array of the same length. This is a standard `numpy` indexing technique which you should be familiar with.\n",
"\n",
"The simplest way to demonstrate our new `mass_term` is with a graphene nanoribbon. First, let's just remind ourselves what a pristine zigzag nanoribbon looks like:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" pb.rectangle(1.2),\n",
" pb.translational_symmetry(a1=True, a2=False)\n",
")\n",
"model.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And let's see its band structure:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from math import pi, sqrt\n",
"\n",
"solver = pb.solver.lapack(model)\n",
"a = graphene.a_cc * sqrt(3)\n",
"bands = solver.calc_bands(-pi/a, pi/a)\n",
"bands.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that the bands touch at zero energy: there is not band gap. Now, let's include the mass term and compute the band structure again."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" pb.rectangle(1.2),\n",
" pb.translational_symmetry(a1=True, a2=False),\n",
" mass_term(delta=2.5) # eV\n",
")\n",
"solver = pb.solver.lapack(model)\n",
"bands = solver.calc_bands(-pi/a, pi/a)\n",
"bands.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We set a very high `delta` value of 2.5 eV for illustration purposes. Indeed, a band gap of 5 eV (`delta * 2`) is quite clearly visible in the band structure.\n",
"\n",
"## PN junction\n",
"\n",
"While we're working with a nanoribbon, let's add a PN junction along its main axis."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def pn_junction(y0, v1, v2):\n",
" @pb.onsite_energy_modifier\n",
" def potential(energy, y):\n",
" energy[y < y0] += v1\n",
" energy[y >= y0] += v2\n",
" return energy\n",
" return potential"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `y0` argument is the position of the junction, while `v1` and `v2` are the values of the potential (in eV) before and after the junction. Let's add it to the nanoribbon:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" pb.rectangle(1.2),\n",
" pb.translational_symmetry(a1=True, a2=False),\n",
" pn_junction(y0=0, v1=-5, v2=5)\n",
")\n",
"model.onsite_map.plot(cmap=\"coolwarm\", site_radius=0.04)\n",
"pb.pltutils.colorbar(label=\"U (eV)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Remember that the [`Model.onsite_map`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Model.html#pybinding.Model.onsite_map) property is a [`StructureMap`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.StructureMap.html#pybinding.StructureMap), which has several plotting methods. A contour plot would not look at all good for such a small nanoribbon, but the method [`StructureMap.plot()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.StructureMap.html#pybinding.StructureMap.plot) is perfect. As before, the ribbon has infinite length along the x-axis and the transparent sites represent the periodic boundaries. The PN junction splits the ribbon in half along its main axis.\n",
"\n",
"We can compute and plot the band structure:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"solver = pb.solver.lapack(model)\n",
"bands = solver.calc_bands(-pi/a, pi/a)\n",
"bands.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, let's create a square potential well. We could define a new modifier function, as before. But lets take a different approach and create the well by composing two PN junctions."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" pb.rectangle(1.2),\n",
" pb.translational_symmetry(a1=True, a2=False),\n",
" pn_junction(y0=-0.2, v1=5, v2=0),\n",
" pn_junction(y0=0.2, v1=0, v2=5)\n",
")\n",
"model.onsite_map.plot(cmap=\"coolwarm\", site_radius=0.04)\n",
"pb.pltutils.colorbar(label=\"U (eV)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It works as expected. This can sometimes be a nice and quick way to extend a model. The square well affects the band structure by breaking electron-hole symmetry:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"solver = pb.solver.lapack(model)\n",
"bands = solver.calc_bands(-pi/a, pi/a)\n",
"bands.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Magnetic field\n",
"\n",
"To model a magnetic field, we need to apply the Peierls substitution:\n",
"\n",
"$$\n",
"t_{nm} \\rightarrow t_{nm} \\text{e}^{i\\frac{2\\pi}{\\Phi_0} \\int_n^m \\vec{A}_{nm} \\cdot d\\vec{l}}\n",
"$$\n",
"\n",
"Here $t_{nm}$ is the hopping energy between two sites, $\\Phi_0 = h/e$ is the magnetic quantum, $h$ is the Planck constant and $\\vec{A}_{nm}$ is the magnetic vector potential along the path between sites $n$ and $m$. We want the magnetic field to be perpendicular to the graphene plane, so we can take the gauge $\\vec{A}(x,y,z) = (By, 0, 0)$.\n",
"\n",
"This can all be expressed with a [`@hopping_energy_modifier`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.hopping_energy_modifier.html#pybinding.hopping_energy_modifier):"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from pybinding.constants import phi0\n",
"\n",
"def constant_magnetic_field(B):\n",
" @pb.hopping_energy_modifier\n",
" def function(energy, x1, y1, x2, y2):\n",
" # the midpoint between two sites\n",
" y = 0.5 * (y1 + y2)\n",
" # scale from nanometers to meters\n",
" y *= 1e-9\n",
"\n",
" # vector potential along the x-axis\n",
" A_x = B * y\n",
"\n",
" # integral of (A * dl) from position 1 to position 2\n",
" peierls = A_x * (x1 - x2)\n",
" # scale from nanometers to meters (because of x1 and x2)\n",
" peierls *= 1e-9\n",
"\n",
" # the Peierls substitution\n",
" return energy * np.exp(1j * 2*pi/phi0 * peierls)\n",
" return function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `energy` argument is the existing hopping energy between two sites at coordinates (`x1`, `y1`) and (`x2`, `y2`). The function computes and returns the Peierls substitution as given by the equation above.\n",
"\n",
"The full signature of a [`@hopping_energy_modifier`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.hopping_energy_modifier.html#pybinding.hopping_energy_modifier) is actually:\n",
"\n",
"```python\n",
"@pb.hopping_energy_modifier\n",
"def function(energy, x1, y1, z1, x2, y2, z2, hop_id):\n",
" return ... # some function of the arguments\n",
"```\n",
"\n",
"The `hop_id` argument tells us which type of hopping it is. Hopping types can be specifically named during the creation of a lattice. This can be used to apply functions only to specific hoppings. However, as with all the modifier arguments, it's optional, so we only take what we need.\n",
"\n",
"To test out our `constant_magnetic_field`, we'll calculate the local density of states (LDOS), where we expect to see peaks corresponding to Landau levels. The computation method used here is explained in detail in the [Kernel polynomial method](http://docs.pybinding.site/page/tutorial/kpm.html) section of the tutorial."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" pb.rectangle(30),\n",
" constant_magnetic_field(B=200) # Tesla\n",
")\n",
"kpm = pb.kpm(model)\n",
"\n",
"ldos = kpm.calc_ldos(energy=np.linspace(-1, 1, 500), broadening=0.015, position=[0, 0])\n",
"ldos.plot()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The values of the magnetic field is exaggerated here (200 Tesla), but that is done to keep the computation time low for the tutorial (less than 0.5 seconds for this LDOS calculation).\n",
"\n",
"## Further reading\n",
"\n",
"Take a look at the [Modifiers](http://docs.pybinding.site/page/tutorial/../api.html#modifiers-api) API reference for more information.\n",
"\n",
"## Example\n",
"\n",
"[Download source code](http://docs.pybinding.site/page/_downloads/fields_example.py)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"PN junction and broken sublattice symmetry in a graphene nanoribbon\"\"\"\n",
"import pybinding as pb\n",
"import matplotlib.pyplot as plt\n",
"from pybinding.repository import graphene\n",
"from math import pi, sqrt\n",
"\n",
"pb.pltutils.use_style()\n",
"\n",
"\n",
"def mass_term(delta):\n",
" \"\"\"Break sublattice symmetry with opposite A and B onsite energy\"\"\"\n",
" @pb.onsite_energy_modifier\n",
" def potential(energy, sub_id):\n",
" energy[sub_id == 'A'] += delta\n",
" energy[sub_id == 'B'] -= delta\n",
" return energy\n",
"\n",
" return potential\n",
"\n",
"\n",
"def pn_juction(y0, v1, v2):\n",
" \"\"\"PN junction potential\n",
"\n",
" The `y0` argument is the position of the junction, while `v1` and `v2`\n",
" are the values of the potential (in eV) before and after the junction.\n",
" \"\"\"\n",
" @pb.onsite_energy_modifier\n",
" def potential(energy, y):\n",
" energy[y < y0] += v1\n",
" energy[y >= y0] += v2\n",
" return energy\n",
"\n",
" return potential\n",
"\n",
"\n",
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" pb.rectangle(1.2), # width in nanometers\n",
" pb.translational_symmetry(a1=True, a2=False),\n",
" mass_term(delta=2.5), # eV\n",
" pn_juction(y0=0, v1=-2.5, v2=2.5) # y0 in [nm] and v1, v2 in [eV]\n",
")\n",
"model.plot()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot the potential: note that pn_junction cancels out delta on some sites\n",
"model.onsite_map.plot(cmap=\"coolwarm\", site_radius=0.04)\n",
"pb.pltutils.colorbar(label=\"U (eV)\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# compute the bands\n",
"solver = pb.solver.lapack(model)\n",
"a = graphene.a_cc * sqrt(3) # nanoribbon unit cell length\n",
"bands = solver.calc_bands(-pi/a, pi/a)\n",
"bands.plot()\n",
"plt.show()"
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}