{
"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": [
"# Eigenvalue solvers\n",
"\n",
"Solvers were first introduced in the [Band structure](http://docs.pybinding.site/page/tutorial/bands.html) section and then used throughout the tutorial to present the results of the various models we constructed. This section will take a more detailed look at the concrete [`lapack()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.lapack) and [`arpack()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.arpack) eigenvalue solvers and their common [`Solver`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver) interface.\n",
"\n",
"[Download this page as a Jupyter notebook](http://docs.pybinding.site/page/_downloads/solvers.ipynb)\n",
"\n",
"## LAPACK\n",
"\n",
"The [`Solver`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver) class establishes the interface of a solver within pybinding, but it does not contain a concrete diagonalization routine. For this reason we never instantiate the plain solver, only its implementations such as [`solver.lapack()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.lapack).\n",
"\n",
"The LAPACK implementation works on dense matrices which makes it well suited only for small systems. However, a great advantage of this solver is that it always solves for all eigenvalues and eigenvectors of a Hamiltonian matrix. This makes it perfect for calculating the entire band structure of the bulk or nanoribbons, as has been shown several times in this tutorial.\n",
"\n",
"Internally, this solver uses the [`scipy.linalg.eigh()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eigh.html#scipy.linalg.eigh) function for dense Hermitian matrices. See the [`solver.lapack()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.lapack) API reference for more details.\n",
"\n",
"## ARPACK\n",
"\n",
"The [`solver.arpack()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.arpack) implementation works on sparse matrices which makes it suitable for large systems. However, only a small subset of the total eigenvalues and eigenvectors can be calculated. This tutorial already contains a few examples where the ARPACK solver was used, and one more is presented below.\n",
"\n",
"Internally, the [`scipy.sparse.linalg.eigsh()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html#scipy.sparse.linalg.eigsh) function is used to solve large sparse Hermitian matrices. The first argument to [`solver.arpack()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.arpack) must be the pybinding [`Model`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Model.html#pybinding.Model), but the following arguments are the same as [`eigsh()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html#scipy.sparse.linalg.eigsh), so the solver routine can be tweaked as desired. Rather than reproduce the full list of options here, we refer you to the scipy [`eigsh()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html#scipy.sparse.linalg.eigsh) reference documentation. Here, we will focus on the specific features of solvers within pybinding.\n",
"\n",
"## Solver interface\n",
"\n",
"No matter which concrete solver is used, they all share a common [`Solver`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver) interface. The two primary properties are [`eigenvalues`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver.eigenvalues) and [`eigenvectors`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver.eigenvectors). These are the raw results of the exact diagonalization of the Hamiltonian matrix."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from pybinding.repository import graphene\n",
"model = pb.Model(graphene.monolayer())\n",
"model.hamiltonian.todense()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"solver = pb.solver.lapack(model)\n",
"solver.eigenvalues"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"solver.eigenvectors"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The properties contain just the raw data. However, [`Solver`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver) also offers a few convenient calculation methods. We'll demonstrate these on a simple rectangular graphene system."
]
},
{
"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(x=3, y=1.2)\n",
")\n",
"model.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, we'll take a look at the [`calc_eigenvalues()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver.calc_eigenvalues) method. While its job is essentially the same as the [`eigenvalues`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver.eigenvalues) property, there is one key difference: the property returns a raw array, while the method returns an [`Eigenvalues`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Eigenvalues.html#pybinding.Eigenvalues) result object. These objects have convenient functions built in and they know how to plot their data:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"solver = pb.solver.arpack(model, k=20) # for the 20 lowest energy eigenvalues\n",
"eigenvalues = solver.calc_eigenvalues()\n",
"eigenvalues.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The basic plot just shows the state number and energy of each eigenstate, but we can also do something more interesting. If we pass a position argument to [`calc_eigenvalues()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver.calc_eigenvalues) it will calculate the probability density $|\\Psi(\\vec{r})|^2$ at that position for each eigenstate and we can view the result using [`Eigenvalues.plot_heatmap()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.Eigenvalues.html#pybinding.Eigenvalues.plot_heatmap):"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"eigenvalues = solver.calc_eigenvalues(map_probability_at=[0.1, 0.6]) # position in [nm]\n",
"eigenvalues.plot_heatmap(show_indices=True)\n",
"pb.pltutils.colorbar()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case we are interested in the probability density at `[x, y] = [0.1, 0.6]`, i.e. a lattice site at the top zigzag edge of our system. Note that the given position does not need to be precise: the probability will be computed for the site closest to the given coordinates. From the figure we can see that the probability at the edge is highest for the two zero-energy states: numbers 9 and 10. We can take a look at the spatial map of state 9 using the [`calc_probability()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver.calc_probability) method:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"probability_map = solver.calc_probability(9)\n",
"probability_map.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The result object in this case is a [`StructureMap`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.StructureMap.html#pybinding.StructureMap) with the probability density $|\\Psi(\\vec{r})|^2$ as its data attribute. As expected, the most prominent states are at the zigzag edges of the system.\n",
"\n",
"An alternative way to get a spatial map of the system is via the local density of states (LDOS). The [`calc_spatial_ldos()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver.calc_spatial_ldos) method makes this easy. The LDOS map is requested for a specific energy value instead of a state number and it considers multiple states within a Gaussian function with the specified broadening:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ldos_map = solver.calc_spatial_ldos(energy=0, broadening=0.05) # [eV]\n",
"ldos_map.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The total density of states can be calculated with [`calc_dos()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver.calc_dos):"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dos = solver.calc_dos(energies=np.linspace(-1, 1, 200), broadening=0.05) # [eV]\n",
"dos.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our example system is quite small so the DOS does not resemble bulk graphene. The zero-energy peak stands out as the signature of the zigzag edge states.\n",
"\n",
"For periodic systems, the wave vector can be controlled using [`Solver.set_wave_vector()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver.set_wave_vector). This allows us to compute the eigenvalues at various points in k-space. For example:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from math import pi\n",
"\n",
"model = pb.Model(\n",
" graphene.monolayer(),\n",
" pb.translational_symmetry()\n",
")\n",
"solver = pb.solver.lapack(model)\n",
"\n",
"kx_lim = pi / graphene.a\n",
"kx_path = np.linspace(-kx_lim, kx_lim, 100)\n",
"ky_outer = 0\n",
"ky_inner = 2*pi / (3*graphene.a_cc)\n",
"\n",
"outer_bands = []\n",
"for kx in kx_path:\n",
" solver.set_wave_vector([kx, ky_outer])\n",
" outer_bands.append(solver.eigenvalues)\n",
"\n",
"inner_bands = []\n",
"for kx in kx_path:\n",
" solver.set_wave_vector([kx, ky_inner])\n",
" inner_bands.append(solver.eigenvalues)\n",
"\n",
"for bands in [outer_bands, inner_bands]:\n",
" result = pb.results.Bands(kx_path, bands)\n",
" result.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This example shows the basic principle of iterating over a path in k-space in order to calculate the band structure. However, this is made much easier with the [`Solver.calc_bands()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver.calc_bands) method. This was already covered in the [Band structure](http://docs.pybinding.site/page/tutorial/bands.html) section and will not be repeated here. But keep in mind that this calculation does not need to be done manually, [`Solver.calc_bands()`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#pybinding.solver.Solver.calc_bands) is the preferred way.\n",
"\n",
"## Further reading\n",
"\n",
"Take a look at the [`solver`](http://docs.pybinding.site/page/tutorial/../_api/pybinding.solver.html#module-pybinding.solver) and `results` reference pages for more detailed information. More solver examples are available throughout this tutorial."
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}