{
"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": [
"# Finite size\n",
"\n",
"This section introduces the concept of shapes with classes [`Polygon`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Polygon.html#pybinding.Polygon) and [`FreeformShape`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.FreeformShape.html#pybinding.FreeformShape) which are used to model systems of finite size. The sparse eigensolver [`arpack()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.arpack) is also introduced as a good tool for exactly solving larger Hamiltonian matrices.\n",
"\n",
"[Download this page as a Jupyter notebook](http://docs.pybinding.site/page/_downloads/finite.ipynb)\n",
"\n",
"## Primitive\n",
"\n",
"The simplest finite-sized system is just the unit cell of the crystal lattice."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from pybinding.repository import graphene\n",
"\n",
"model = pb.Model(graphene.monolayer())\n",
"model.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The unit cell can also be replicated a number of times to create a bigger system."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" pb.primitive(a1=5, a2=3)\n",
")\n",
"model.plot()\n",
"model.lattice.plot_vectors(position=[0.6, -0.25])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The [`primitive()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.primitive.html#pybinding.primitive) parameter tells the model to replicate the unit cell 5 times in the $a_1$ vector direction and 3 times in the $a_2$ direction. However, to model realistic systems we need proper shapes.\n",
"\n",
"## Polygon\n",
"\n",
"The easiest way to create a 2D shape is with the [`Polygon`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Polygon.html#pybinding.Polygon) class. For example, a simple rectangle:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def rectangle(width, height):\n",
" x0 = width / 2\n",
" y0 = height / 2\n",
" return pb.Polygon([[x0, y0], [x0, -y0], [-x0, -y0], [-x0, y0]])\n",
"\n",
"shape = rectangle(1.6, 1.2)\n",
"shape.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A [`Polygon`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Polygon.html#pybinding.Polygon) is initialized with a list of vertices which should be given in clockwise or counterclockwise order. When added to a [`Model`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Model.html#pybinding.Model) the lattice will expand to fill the shape."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" rectangle(width=1.6, height=1.2)\n",
")\n",
"model.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To help visualize the shape and the expanded lattice, the polygon outline can be plotted on top of the system by calling both plot methods one after another."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def trapezoid(a, b, h):\n",
" return pb.Polygon([[-a/2, 0], [-b/2, h], [b/2, h], [a/2, 0]])\n",
"\n",
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" trapezoid(a=3.2, b=1.4, h=1.5)\n",
")\n",
"model.plot()\n",
"model.shape.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In general, a shape does not depend on a specific material, so it can be easily reused. Here, we shall switch to a `graphene.bilayer()` lattice, but we'll keep the same `trapezoid` shape as defined earlier:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = pb.Model(\n",
" graphene.bilayer(),\n",
" trapezoid(a=3.2, b=1.4, h=1.5)\n",
")\n",
"model.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Freeform shape\n",
"\n",
"Unlike a [`Polygon`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Polygon.html#pybinding.Polygon) which is defined by a list of vertices, a [`FreeformShape`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.FreeformShape.html#pybinding.FreeformShape) is defined by a `contains` function which determines if a lattice site is inside the desired shape."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def circle(radius):\n",
" def contains(x, y, z):\n",
" return np.sqrt(x**2 + y**2) < radius\n",
" return pb.FreeformShape(contains, width=[2*radius, 2*radius])\n",
"\n",
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" circle(radius=2.5)\n",
")\n",
"model.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `width` parameter of [`FreeformShape`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.FreeformShape.html#pybinding.FreeformShape) specifies the bounding box width. Only sites inside the bounding box will be considered for the shape. It's like carving a sculpture from a block of stone. The bounding box can be thought of as the stone block, while the `contains` function is the carving tool that can give the fine detail of the shape.\n",
"\n",
"As with [`Polygon`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Polygon.html#pybinding.Polygon), we can visualize the shape with the [`FreeformShape.plot()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.FreeformShape.html#pybinding.FreeformShape.plot) method."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def ring(inner_radius, outer_radius):\n",
" def contains(x, y, z):\n",
" r = np.sqrt(x**2 + y**2)\n",
" return np.logical_and(inner_radius < r, r < outer_radius)\n",
" return pb.FreeformShape(contains, width=[2*outer_radius, 2*outer_radius])\n",
"\n",
"shape = ring(inner_radius=1.4, outer_radius=2)\n",
"shape.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The shaded area indicates the shape as determined by the `contains` function. Creating a model will cause the lattice to fill in the shape."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" ring(inner_radius=1.4, outer_radius=2)\n",
")\n",
"model.plot()\n",
"model.shape.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that the `ring` example uses `np.logical_and` instead of the plain `and` keyword. This is because the `x, y, z` positions are not given as scalar numbers but as numpy arrays. Array comparisons return boolean arrays:\n",
"\n",
"```python\n",
">>> x = np.array([7, 2, 3, 5, 1])\n",
">>> x < 5\n",
"[False, True, True, False, True]\n",
">>> 2 < x and x < 5\n",
"ValueError: ...\n",
">>> np.logical_and(2 < x, x < 5)\n",
"[False, False, True, False, False]\n",
"```\n",
"\n",
"The `and` keyword can only operate on scalar values, but `np.logical_and` can consider arrays. Likewise, `math.sqrt` does not work with arrays, but `np.sqrt` does.\n",
"\n",
"## Composite shape\n",
"\n",
"Complicated system geometry can also be produced by composing multiple simple shapes. The following example gives a quick taste of how it works. For a full overview of this functionality, see the [Composite shapes](http://docs.pybinding.site/page/tutorial/../advanced/shapes.html) section."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Simple shapes\n",
"rectangle = pb.rectangle(x=6, y=1)\n",
"hexagon = pb.regular_polygon(num_sides=6, radius=1.92, angle=np.pi/6)\n",
"circle = pb.circle(radius=0.6)\n",
"\n",
"# Compose them naturally\n",
"shape = rectangle + hexagon - circle\n",
"\n",
"model = pb.Model(graphene.monolayer(), shape)\n",
"model.shape.plot()\n",
"model.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Spatial LDOS\n",
"\n",
"Now that we have a ring structure, we can exactly diagonalize its `model.hamiltonian` using a [`Solver`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver). We previously used the [`lapack()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.lapack) solver to find all the eigenvalues and eigenvectors, but this is not efficient for larger systems. The sparse [`arpack()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.arpack) solver can calculate a targeted subset of the eigenvalues, which is usually desired and much faster. In this case, we are interested only in the 20 lowest energy states."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" ring(inner_radius=1.4, outer_radius=2)\n",
")\n",
"solver = pb.solver.arpack(model, k=20) # only the 20 lowest eigenstates\n",
"\n",
"ldos = solver.calc_spatial_ldos(energy=0, broadening=0.05) # eV\n",
"ldos.plot(site_radius=(0.03, 0.12))\n",
"pb.pltutils.colorbar(label=\"LDOS\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The convenient [`Solver.calc_spatial_ldos()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver.calc_spatial_ldos) method calculates the local density of states (LDOS) at every site for the given energy with a Gaussian broadening. The returned object is a [`StructureMap`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.StructureMap.html#pybinding.StructureMap) which holds the LDOS data. The [`StructureMap.plot()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.StructureMap.html#pybinding.StructureMap.plot) method will produce a figure similar to [`Model.plot()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Model.html#pybinding.Model.plot), but with a colormap indicating the LDOS value at each lattice site. In addition, the `site_radius` argument specifies a range of sizes which will cause the low intensity sites to appear as small circles while high intensity ones become large. The states with a high LDOS are clearly visible on the outer and inner edges of the graphene ring structure.\n",
"\n",
"## Further reading\n",
"\n",
"For more finite-sized systems check out the [examples section](http://docs.pybinding.site/page/tutorial/../examples/finite/index.html).\n",
"\n",
"## Example\n",
"\n",
"[Donwload source code](http://docs.pybinding.site/page/_downloads/finite_example.py)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"Model a graphene ring structure and calculate the local density of states\"\"\"\n",
"import pybinding as pb\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from pybinding.repository import graphene\n",
"\n",
"pb.pltutils.use_style()\n",
"\n",
"\n",
"def ring(inner_radius, outer_radius):\n",
" \"\"\"A simple ring shape\"\"\"\n",
" def contains(x, y, z):\n",
" r = np.sqrt(x**2 + y**2)\n",
" return np.logical_and(inner_radius < r, r < outer_radius)\n",
"\n",
" return pb.FreeformShape(contains, width=[2 * outer_radius, 2 * outer_radius])\n",
"\n",
"\n",
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" ring(inner_radius=1.4, outer_radius=2) # length in nanometers\n",
")\n",
"\n",
"model.plot()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# only solve for the 20 lowest energy eigenvalues\n",
"solver = pb.solver.arpack(model, k=20)\n",
"ldos = solver.calc_spatial_ldos(energy=0, broadening=0.05) # LDOS around 0 eV\n",
"ldos.plot(site_radius=(0.03, 0.12))\n",
"pb.pltutils.colorbar(label=\"LDOS\")\n",
"plt.show()"
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}