{
"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": [
"# Multi-orbital models\n",
"\n",
"In pybinding, if an onsite or hopping energy term is defined as a matrix (instead of a scalar), we refer to the resulting model as multi-orbital. The elements of the matrix term may correspond to different spins, electrons and holes, or any other degrees of freedom. These can have different physical meaning depending on the intend of the model. Because we're talking in generic terms here, we'll use orbital as a blanket term to refer to any degree of freedom, i.e. matrix element of an onsite or hopping term.\n",
"\n",
"This section describes how these models can be defined and how the presence of multiple orbitals affects modifier functions and the results obtained from solvers. In general, it is as simple as replacing a scalar value with a matrix while all of the principals described in the [Tutorial](http://docs.pybinding.site/page/advanced/../tutorial/index.html) still apply.\n",
"\n",
"[Download this page as a Jupyter notebook](http://docs.pybinding.site/page/_downloads/multiorbital.ipynb)\n",
"\n",
"## Onsite and hopping matrices\n",
"\n",
"Starting from the very beginning, the orbital count of a site is determined by the shape of the onsite energy matrix. Let's take a look at a few possibilities:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"lat = pb.Lattice([1, 0], [0, 1])\n",
"lat.add_sublattices(\n",
" (\"A\", [0.0, 0.0], 0.5), # single-orbital: scalar\n",
" (\"B\", [0.0, 0.2], [[1.5, 2j], # two-orbital: 2x2 Hermitian matrix\n",
" [-2j, 1.5]]),\n",
" (\"C\", [0.3, 0.1], np.zeros(2)), # two-orbital: zero onsite term\n",
" (\"D\", [0.1, 0.0], [[4, 0, 0], # three-orbital: only diagonal\n",
" [0, 5, 0],\n",
" [0, 0, 6]]),\n",
" (\"E\", [0.2, 0.2], [4, 5, 6]) # three-orbital: only diagonal, terse notation\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The onsite term is required to be a square Hermitian matrix. If a 1D array is given instead of a matrix, it will be interpreted as the main diagonal of a square matrix (see sublattices D and E which have identical onsite term specified with different notations).\n",
"\n",
"As seen above, sublattices don't need to all have the same orbital count. The only thing to keep in mind is that the hopping matrix which connect a pair of sublattice sites must have the appropriate shape: the number of rows must match the orbital count of the source sublattice and the number of columns must match the destination sublattice."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"lat.add_hoppings(\n",
" ([0, 1], \"A\", \"A\", 1.2), # scalar\n",
" ([0, 1], \"B\", \"B\", [[1, 2], # 2x2\n",
" [3, 4]]),\n",
" ([0, 0], \"B\", \"C\", [[2j, 0], # 2x2\n",
" [1j, 0]]),\n",
" ([0, 0], \"A\", \"D\", [[1, 2, 3]]), # 1x3\n",
" ([0, 1], \"D\", \"A\", [[7], # 3x1\n",
" [8],\n",
" [9]]),\n",
" ([0, 0], \"B\", \"D\", [[1j, 0, 0], # 2x3\n",
" [2, 0, 3j]])\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If a matrix of the wrong shape is given, an informative error is raised:\n",
"\n",
"```python\n",
">>> lat.add_one_hopping([0, 0], \"A\", \"B\", 0.6)\n",
"RuntimeError: Hopping size mismatch: from 'A' (1) to 'B' (2) with matrix (1, 1)\n",
">>> lat.add_one_hopping([0, 1], \"D\", \"D\", [[1, 2, 3],\n",
"... [4, 5, 6]])\n",
"RuntimeError: Hopping size mismatch: from 'D' (3) to 'D' (3) with matrix (2, 3)\n",
"```\n",
"\n",
"After the [`Lattice`](http://docs.pybinding.site/page/advanced/../_api/pybinding.Lattice.html#pybinding.Lattice) is complete, a [`Model`](http://docs.pybinding.site/page/advanced/../_api/pybinding.Model.html#pybinding.Model) can be built as usual:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = pb.Model(lat, pb.primitive(2, 2))\n",
"model.system.num_sites"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model.hamiltonian.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Sites refer to physical locations so their total count corresponds to the number of sublattices (A to E) multiplied by the number of times the unit cell is repeated. The Hamiltonian matrix is larger than `num_sites` due to the extra orbitals.\n",
"\n",
"## Effect on modifier functions\n",
"\n",
"The [`@onsite_energy_modifier`](http://docs.pybinding.site/page/advanced/../_api/pybinding.onsite_energy_modifier.html#pybinding.onsite_energy_modifier) and [`@hopping_energy_modifier`](http://docs.pybinding.site/page/advanced/../_api/pybinding.hopping_energy_modifier.html#pybinding.hopping_energy_modifier) functions work equally well for single- and multi-orbital models. In case of the latter, the `energy` argument of the modifiers will have a shape matching the onsite/hopping matrix term."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"@pb.onsite_energy_modifier\n",
"def potential(energy, x):\n",
" \"\"\"Linear onsite potential as a function of x for a 2-orbital model\"\"\"\n",
" return energy + np.eye(2) * x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note the [`np.eye(2)`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.eye.html#numpy.eye) in the code above. The number 2 matches the 2-orbital structure of a specific model. Without this, `energy + x` would also add the value to the off-diagonal elements of the onsite matrix which is not desirable in this case.\n",
"\n",
"The modifier defined above will only work for 2-orbital models. In general, we might want to create modifiers which work with any n-orbital model or with a mixed number of orbitals. For this we can use the `sub_id` modifier argument and its `.eye` attribute which supplies the correct matrix shape for any sublattice:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"@pb.onsite_energy_modifier\n",
"def potential(energy, x, sub_id):\n",
" \"\"\"Same as above, but works for any n-orbital model\"\"\"\n",
" return energy + sub_id.eye * x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Even more generally, if we wish to apply completely different functions to the various sublattices, the `sub_id` argument can be used to create different branches in the modifier:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"@pb.onsite_energy_modifier\n",
"def potential(energy, x, sub_id):\n",
" \"\"\"Applies different functions to different sublattices\"\"\"\n",
" if sub_id == \"A\":\n",
" return energy + x # we know sublattice A is single-orbital\n",
" elif sub_id == \"D\":\n",
" energy[x > 0] += sub_id.eye * x # the notation can be mixed with numpy indexing\n",
" return energy # apply only to sites where x > 0\n",
" elif sub_id == \"B\":\n",
" sigma_y = np.array([[0, -1j],\n",
" [1j, 0]])\n",
" return energy + sigma_y * 1.3 - np.eye(2) * 0.6 # add multiple 2x2 matrices\n",
" else:\n",
" return energy # leave the other sublattices unchanged"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This branching behavior is only supported by the `sub_id` and `hop_id` arguments. Do not try to create branches like this using any of the other modifier arguments:\n",
"\n",
"```python\n",
"\"\"\"Creating a position-dependent potential\"\"\"\n",
"# This is an error with anything except sub_id or hop_id\n",
"if x > 0:\n",
" return energy + 1\n",
"else:\n",
" return energy - 1\n",
"```\n",
"\n",
"```python\n",
"# Use this notation instead\n",
"energy[x > 0] += 1\n",
"energy[x <= 0] -= 1\n",
"```\n",
"\n",
"On the other hand, `sub_id` and `hop_id` can be used with either of these variants with just a single caveat:\n",
"\n",
"```python\n",
"\"\"\"Sublattice-dependent potential\"\"\"\n",
"# This always works with sub_id and hop_id\n",
"if sub_id == \"A\":\n",
" return energy + 1\n",
"else:\n",
" return energy - 1\n",
"```\n",
"\n",
"```python\n",
"# This only works when all sublattices have the same number of orbitals,\n",
"# but it will raise an error for mixed orbital counts.\n",
"energy[sub_id == \"A\"] += 1\n",
"energy[sub_id == \"B\"] -= 1\n",
"```\n",
"\n",
"## Local properties and plotting\n",
"\n",
"When examining the local properties of a multi-orbital model, it is important to make the distinction between system indices which correspond to sites (unique positions) and Hamiltonian indices which correspond to the onsite or hopping terms in the Hamiltonian.\n",
"\n",
"As shown in one of the previous examples, the number of sites in a system does not have to be equal to the size of the Hamiltonian matrix (`hamiltonian.shape[0] >= num_sites`). This affects how the system and Hamiltonian are indexed. System indices are always scalars and point to a single site position. For single-orbital models there is a 1:1 correspondence between system and Hamiltonian indices. However, for multi-orbital models the Hamiltonian indices are 1D arrays with a size corresponding to the number of orbitals on the target site."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = pb.Model(lat, pb.primitive(2, 2))\n",
"sys_idx = model.system.find_nearest(position=[0, 0], sublattice=\"D\")\n",
"sys_idx # <-- Points to a site on sublattice D which is closest to the target position."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model.system.x[sys_idx]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model.system.y[sys_idx]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ham_idx = model.system.to_hamiltonian_indices(sys_idx)\n",
"ham_idx # <-- Array of integers which can be used to index the Hamiltonian matrix."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ham = model.hamiltonian.todense()\n",
"ham[np.ix_(ham_idx, ham_idx)] # Returns the onsite hopping term of sublattice D."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Functions which compute various local properties take into account the presence of multiple orbitals on a single site. For example, when calculating the local density of states, one of the input parameters is the target site position. By default, the resulting LDOS is calculated as the sum of all orbitals but this is optional as shown in the following example:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"Calculate the LDOS in the center of a MoS2 quantum dot\"\"\"\n",
"from pybinding.repository import group6_tmd\n",
"\n",
"model = pb.Model(group6_tmd.monolayer_3band(\"MoS2\"),\n",
" pb.regular_polygon(6, 20))\n",
"\n",
"kpm = pb.kpm(model)\n",
"energy = np.linspace(-1, 3.8, 500)\n",
"broadening = 0.05\n",
"position = [0, 0]\n",
"\n",
"plt.figure(figsize=(6.7, 2.3))\n",
"\n",
"plt.subplot(121, title=\"Reduced -- sum of all orbitals\")\n",
"ldos = kpm.calc_ldos(energy, broadening, position)\n",
"ldos.plot(color=\"C1\")\n",
"\n",
"plt.subplot(122, title=\"Individual orbitals\")\n",
"ldos = kpm.calc_ldos(energy, broadening, position, reduce=False)\n",
"ldos.plot()"
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}