From 314cdf8026efbf0ebe5be80f8772695b5f623af3 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Thu, 7 Jan 2021 14:31:55 +0100 Subject: [PATCH 1/4] start implementation of period finding algorithm --- example_notebooks/bloch_generator.ipynb | 1538 +++++++++++++++++++++++ 1 file changed, 1538 insertions(+) create mode 100644 example_notebooks/bloch_generator.ipynb diff --git a/example_notebooks/bloch_generator.ipynb b/example_notebooks/bloch_generator.ipynb new file mode 100644 index 0000000..2c1d7aa --- /dev/null +++ b/example_notebooks/bloch_generator.ipynb @@ -0,0 +1,1538 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import sympy\n", + "from collections import OrderedDict" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import qsymm" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Graphene" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We use the Hamiltonian generator to reproduce the spinless nearest neighbour tight binding Hamiltonian for graphene.\n", + "\n", + "The generators of the symmetry group are time-reversal symmetry, sublattice (or chiral) symmetry, and threefold rotation symmetry." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Time reversal\n", + "TR = qsymm.time_reversal(2, U=np.eye(2))\n", + "\n", + "# Chiral symmetry\n", + "C = qsymm.chiral(2, U=np.array([[1, 0], [0, -1]]))\n", + "\n", + "# Atom A rotates into A, B into B, use exact sympy representation\n", + "sphi = 2*sympy.pi/3\n", + "RC3 = sympy.Matrix([[sympy.cos(sphi), -sympy.sin(sphi)],\n", + " [sympy.sin(sphi), sympy.cos(sphi)]])\n", + "C3 = qsymm.PointGroupElement(RC3, False, False, np.eye(2))\n", + "\n", + "symmetries = [C, TR, C3]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are two carbon atoms per unit cell (A and B) with one orbital each. The lattice is triangular, and only includes hoppings between nearest neighbour atoms. This restricts hoppings to only those between atoms of different types, such that each atom couples to three neighbouring atoms. Using the symmetrization strategy to generate the Hamiltonian, it is sufficient to specify hoppings to one such neighbour along with the symmetry generators, and we take the vector $(1, 0)$ to connect this neighbouring pair of atoms." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "norbs = OrderedDict([('A', 1), ('B', 1)]) # A and B atom per unit cell, one orbital each\n", + "hopping_vectors = [('A', 'B', [0, 1])] # Hopping between neighbouring A and B atoms" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Generate the Hamiltonian." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "family = qsymm.bloch_family(hopping_vectors, symmetries, norbs)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & e^{i k_{y}} + e^{- \\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} + e^{- \\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}}\\\\e^{\\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} + e^{\\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}} + e^{- i k_{y}} & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, e**(I*k_y) + e**(-I*k_y/2)*e**(sqrt(3)*I*k_x/2) + e**(-I*k_y/2)*e**(-sqrt(3)*I*k_x/2)],\n", + "[e**(I*k_y/2)*e**(sqrt(3)*I*k_x/2) + e**(I*k_y/2)*e**(-sqrt(3)*I*k_x/2) + e**(-I*k_y), 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "qsymm.display_family(family)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Scale the bond length in terms of the graphene lattice constant, and have the function return a list of BlochModel objects. For this we can use a more convenient definition of the rotation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "C3 = qsymm.rotation(1/3, U=np.eye(2))\n", + "symmetries = [C, TR, C3]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "norbs = OrderedDict([('A', 1), ('B', 1)]) # A and B atom per unit cell, one orbital each\n", + "hopping_vectors = [('A', 'B', [0, 1/np.sqrt(3)])] # Hopping between neighbouring A and B atoms" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "family = qsymm.bloch_family(hopping_vectors, symmetries, norbs, bloch_model=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "qsymm.display_family(family)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Three-orbital tight binding model for monolayer MX$_2$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We use the Hamiltonian generator to reproduce the tight binding model for monolayer MX$_2$ published in Phys. Rev. B 88, 085433 (2013).\n", + "\n", + "The generators of the symmetry group of the tight binding model are time reversal symmetry, mirror symmetry and threefold rotation symmetry." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "# Time reversal\n", + "TR = qsymm.time_reversal(2, np.eye(3))\n", + "\n", + "# Mirror symmetry \n", + "Mx = qsymm.mirror([1, 0], np.diag([1, -1, 1]))\n", + "\n", + "# Threefold rotation on d_z^2, d_xy, d_x^2-y^2 states.\n", + "C3U = np.array([[1, 0, 0],\n", + " [0, -0.5, -np.sqrt(3)/2],\n", + " [0, np.sqrt(3)/2, -0.5]])\n", + "# Could also use the predefined representation of rotations on d-orbitals\n", + "Ld = qsymm.groups.L_matrices(3, 2)\n", + "C3U2 = qsymm.groups.spin_rotation(2 * np.pi * np.array([0, 0, 1/3]), Ld)\n", + "# Restrict to d_z^2, d_xy, d_x^2-y^2 states\n", + "mask = np.array([1, 2 ,0])\n", + "C3U2 = C3U2[mask][:, mask]\n", + "assert np.allclose(C3U, C3U2)\n", + "\n", + "C3 = qsymm.rotation(1/3, U=C3U)\n", + "\n", + "symmetries = [TR, Mx, C3]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we specify the hoppings to include. The tight binding model has a triangular lattice, three orbitals per M atom, and nearest neighbour hopping." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "# One site per unit cell (M atom), with three orbitals\n", + "norbs = OrderedDict([('a', 3)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Each atom has six nearest neighbour atoms at a distance of one primitive lattice vector. Since we use the symmetrization strategy to generate the Hamiltonian, it is sufficient to specify a hopping to one nearest neighbour atom along with the symmetry generators. We take the primitive vector connecting the pair of atoms to be $(1, 0)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "# Hopping to a neighbouring atom one primitive lattice vector away\n", + "hopping_vectors = [('a', 'a', [1, 0])]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Generate the tight binding Hamiltonian." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "family = qsymm.bloch_family(hopping_vectors, symmetries, norbs, bloch_model=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Hamiltonian family should include 8 linearly independent components, including the onsite terms." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "len(family)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "qsymm.display_family(family, nsimplify=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 4-site model for monolayer WTe$_2$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We use the Hamiltonian generator to reproduce the tight binding model for monolayer WTe$_2$ published in Phys. Rev. X 6, 041069 (2016).\n", + "\n", + "The generators of the symmetry group of the tight binding model are time reversal symmetry, glide reflection and inversion symmetry." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "# Define 4 sites with one orbital each\n", + "sites = ['Ad', 'Ap', 'Bd', 'Bp']\n", + "norbs = OrderedDict([(site, 1) for site in sites])\n", + "\n", + "# Define symbolic coordinates for orbitals\n", + "rAp = qsymm.sympify('[x_Ap, y_Ap]')\n", + "rAd = qsymm.sympify('[x_Ad, y_Ad]')\n", + "rBp = qsymm.sympify('[x_Bp, y_Bp]')\n", + "rBd = qsymm.sympify('[x_Bd, y_Bd]')\n", + "\n", + "# Define hoppings to include\n", + "hopping_vectors = [('Bd', 'Bd', np.array([1, 0])),\n", + " ('Ap', 'Ap', np.array([1, 0])),\n", + " ('Bd', 'Ap', rAp - rBd),\n", + " ('Ap', 'Bp', rBp - rAp),\n", + " ('Ad', 'Bd', rBd - rAd), \n", + " ]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "# Inversion\n", + "perm_inv = {'Ad': 'Bd', 'Ap': 'Bp', 'Bd': 'Ad', 'Bp': 'Ap'}\n", + "onsite_inv = {site: (1 if site in ['Ad', 'Bd'] else -1) * np.eye(1) for site in sites}\n", + "inversion = qsymm.groups.symmetry_from_permutation(-np.eye(2), perm_inv, norbs, onsite_inv)\n", + "\n", + "# Glide\n", + "perm_glide = {site: site for site in sites}\n", + "onsite_glide = {site: (1 if site in ['Ad', 'Bd'] else -1) * np.eye(1) for site in sites}\n", + "glide = qsymm.groups.symmetry_from_permutation(np.array([[-1, 0],[0, 1]]), perm_glide, norbs, onsite_glide)\n", + "\n", + "# TR\n", + "time_reversal = qsymm.time_reversal(2, np.eye(4))\n", + "\n", + "gens = {glide, inversion, time_reversal}\n", + "sg = qsymm.groups.generate_group(gens)\n", + "len(sg)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Generate the tight binding Hamiltonian." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "family = qsymm.bloch_family(hopping_vectors, gens, norbs=norbs)\n", + "print(len(family), [len(fam) for fam in family])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "qsymm.display_family(family)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Square lattice with 4 sites in the UC" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Make model with square lattice that has 4 sites in the unit cell related by 4-fold rotation. Sites have spin-1/2 and we add time reversal and particle-hole symmetry." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "hopping_vectors = [('a', 'b', np.array([1, 0])),\n", + " ('b', 'a', np.array([1, 0])),\n", + " ('c', 'd', np.array([1, 0])),\n", + " ('d', 'c', np.array([1, 0])),\n", + " ('a', 'c', np.array([0, 1])),\n", + " ('c', 'a', np.array([0, 1])),\n", + " ('b', 'd', np.array([0, 1])),\n", + " ('d', 'b', np.array([0, 1])), \n", + " ]\n", + "\n", + "# Define spin-1/2 operators\n", + "S = qsymm.groups.spin_matrices(1/2)\n", + "# Define real space rotation generators in 2D\n", + "L = qsymm.groups.L_matrices(d=2)\n", + "\n", + "sites = ['a', 'b', 'c', 'd']\n", + "norbs = OrderedDict([(site, 2) for site in sites])\n", + "\n", + "perm_C4 = {'a': 'b', 'b': 'd', 'd': 'c', 'c': 'a'}\n", + "onsite_C4 = {site: qsymm.groups.spin_rotation(2*np.pi * np.array([0, 0, 1/4]), S) for site in sites}\n", + "C4 = qsymm.groups.symmetry_from_permutation(\n", + " qsymm.groups.spin_rotation(2*np.pi * np.array([1/4]), L, roundint=True), \n", + " perm_C4, norbs, onsite_C4)\n", + "\n", + "# Fermionic TR\n", + "time_reversal = qsymm.time_reversal(2, \n", + " np.kron(np.eye(4), qsymm.groups.spin_rotation(2*np.pi * np.array([0, 1/2, 0]), S)))\n", + "\n", + "# define strange PH symmetry\n", + "particle_hole = qsymm.particle_hole(2, np.eye(8))\n", + "\n", + "gens = {C4, time_reversal, particle_hole}\n", + "\n", + "sg = qsymm.groups.generate_group(gens)\n", + "len(sg)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "family = qsymm.bloch_family(hopping_vectors, gens, norbs=norbs)\n", + "print(len(family), [len(fam) for fam in family])\n", + "qsymm.display_family(family)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### C4T symmetric model" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "8" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "hopping_vectors = [('a', 'a', np.array([1, 0])), \n", + " ]\n", + "\n", + "# Define spin-1/2 operators\n", + "S = qsymm.groups.spin_matrices(1/2)\n", + "# Define real space rotation generators in 2D\n", + "L = qsymm.groups.L_matrices(d=2)\n", + "\n", + "sites = ['a']\n", + "norbs = OrderedDict([(site, 4) for site in sites])\n", + "\n", + "perm_C4 = {'a': 'a'}\n", + "onsite_C4 = {site: np.kron(np.eye(2),\n", + " qsymm.groups.spin_rotation(2*np.pi * np.array([0, 0, 1/4]), S)\n", + " )\n", + " for site in sites}\n", + "C4 = qsymm.groups.symmetry_from_permutation(\n", + " qsymm.groups.spin_rotation(2*np.pi * np.array([1/4]), L, roundint=True), \n", + " perm_C4, norbs, onsite_C4)\n", + "\n", + "# Fermionic TR\n", + "T = qsymm.time_reversal(2, \n", + " np.kron(np.eye(2),\n", + " qsymm.groups.spin_rotation(2*np.pi * np.array([0, 1/2, 0]), S),\n", + " )\n", + " )\n", + "\n", + "C4T = C4 * T\n", + "\n", + "sg = qsymm.groups.generate_group(gens)\n", + "len(sg)" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "{1, R(π), R(π/2), R(-π/2), T, R(π) T, R(π/2) T, R(-π/2) T}" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sg" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "14 [1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]\n" + ] + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0 & 0\\\\0 & 1 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[1, 0, 0, 0],\n", + "[0, 1, 0, 0],\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 1 & 0\\\\0 & 0 & 0 & 1\\\\1 & 0 & 0 & 0\\\\0 & 1 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 0, 1, 0],\n", + "[0, 0, 0, 1],\n", + "[1, 0, 0, 0],\n", + "[0, 1, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 1 & 0\\\\0 & 0 & 0 & 1\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, 0],\n", + "[0, 0, 1, 0],\n", + "[0, 0, 0, 1]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & i & 0\\\\0 & 0 & 0 & - i\\\\- i & 0 & 0 & 0\\\\0 & i & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, I, 0],\n", + "[ 0, 0, 0, -I],\n", + "[-I, 0, 0, 0],\n", + "[ 0, I, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0 & 0\\\\0 & e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x), 0, 0, 0],\n", + "[ 0, e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x), 0, 0],\n", + "[ 0, 0, 0, 0],\n", + "[ 0, 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0\\\\i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x), 0, 0],\n", + "[I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0, 0, 0],\n", + "[ 0, 0, 0, 0],\n", + "[ 0, 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}} & 0\\\\0 & 0 & 0 & e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}}\\\\e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0 & 0\\\\0 & e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x), 0],\n", + "[ 0, 0, 0, e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x)],\n", + "[e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x), 0, 0, 0],\n", + "[ 0, e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x), 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}}\\\\0 & 0 & i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0\\\\0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0\\\\i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x)],\n", + "[ 0, 0, I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0],\n", + "[ 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x), 0, 0],\n", + "[I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}} & 0\\\\0 & 0 & 0 & e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}}\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, 0],\n", + "[0, 0, e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x), 0],\n", + "[0, 0, 0, e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x)]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}}\\\\0 & 0 & i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x)],\n", + "[0, 0, I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & - e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\\\e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}} & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, -e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x), 0, 0],\n", + "[e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x), 0, 0, 0],\n", + "[ 0, 0, 0, 0],\n", + "[ 0, 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & i e^{i k_{x}} + i e^{i k_{y}} + i e^{- i k_{y}} + i e^{- i k_{x}} & 0\\\\0 & 0 & 0 & - i e^{i k_{x}} - i e^{i k_{y}} - i e^{- i k_{y}} - i e^{- i k_{x}}\\\\- i e^{i k_{x}} - i e^{i k_{y}} - i e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0 & 0\\\\0 & i e^{i k_{x}} + i e^{i k_{y}} + i e^{- i k_{y}} + i e^{- i k_{x}} & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, I*e**(I*k_x) + I*e**(I*k_y) + I*e**(-I*k_y) + I*e**(-I*k_x), 0],\n", + "[ 0, 0, 0, -I*e**(I*k_x) - I*e**(I*k_y) - I*e**(-I*k_y) - I*e**(-I*k_x)],\n", + "[-I*e**(I*k_x) - I*e**(I*k_y) - I*e**(-I*k_y) - I*e**(-I*k_x), 0, 0, 0],\n", + "[ 0, I*e**(I*k_x) + I*e**(I*k_y) + I*e**(-I*k_y) + I*e**(-I*k_x), 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & - e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}}\\\\0 & 0 & e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}} & 0\\\\0 & - e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\\\e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}} & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, 0, -e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x)],\n", + "[ 0, 0, e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x), 0],\n", + "[ 0, -e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x), 0, 0],\n", + "[e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x), 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & - e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}}\\\\0 & 0 & e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}} & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, -e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x)],\n", + "[0, 0, e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x), 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "familyTR = qsymm.bloch_family(hopping_vectors, {T, C4}, norbs=norbs)\n", + "print(len(familyTR), [len(fam) for fam in familyTR])\n", + "qsymm.display_family(familyTR)" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "20 [1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]\n" + ] + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0 & 0\\\\0 & 1 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[1, 0, 0, 0],\n", + "[0, 1, 0, 0],\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 1 & 0\\\\0 & 0 & 0 & 1\\\\1 & 0 & 0 & 0\\\\0 & 1 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 0, 1, 0],\n", + "[0, 0, 0, 1],\n", + "[1, 0, 0, 0],\n", + "[0, 1, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 1 & 0\\\\0 & 0 & 0 & 1\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, 0],\n", + "[0, 0, 1, 0],\n", + "[0, 0, 0, 1]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & i & 0\\\\0 & 0 & 0 & - i\\\\- i & 0 & 0 & 0\\\\0 & i & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, I, 0],\n", + "[ 0, 0, 0, -I],\n", + "[-I, 0, 0, 0],\n", + "[ 0, I, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}e^{i k_{y}} + e^{- i k_{y}} & 0 & 0 & 0\\\\0 & e^{i k_{x}} + e^{- i k_{x}} & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[e**(I*k_y) + e**(-I*k_y), 0, 0, 0],\n", + "[ 0, e**(I*k_x) + e**(-I*k_x), 0, 0],\n", + "[ 0, 0, 0, 0],\n", + "[ 0, 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0\\\\i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x), 0, 0],\n", + "[I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0, 0, 0],\n", + "[ 0, 0, 0, 0],\n", + "[ 0, 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & e^{i k_{y}} + e^{- i k_{y}} & 0\\\\0 & 0 & 0 & e^{i k_{x}} + e^{- i k_{x}}\\\\e^{i k_{y}} + e^{- i k_{y}} & 0 & 0 & 0\\\\0 & e^{i k_{x}} + e^{- i k_{x}} & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, e**(I*k_y) + e**(-I*k_y), 0],\n", + "[ 0, 0, 0, e**(I*k_x) + e**(-I*k_x)],\n", + "[e**(I*k_y) + e**(-I*k_y), 0, 0, 0],\n", + "[ 0, e**(I*k_x) + e**(-I*k_x), 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & e^{i k_{y}} - e^{- i k_{y}}\\\\0 & 0 & i e^{i k_{x}} - i e^{- i k_{x}} & 0\\\\0 & i e^{i k_{x}} - i e^{- i k_{x}} & 0 & 0\\\\- e^{i k_{y}} + e^{- i k_{y}} & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, 0, e**(I*k_y) - e**(-I*k_y)],\n", + "[ 0, 0, I*e**(I*k_x) - I*e**(-I*k_x), 0],\n", + "[ 0, I*e**(I*k_x) - I*e**(-I*k_x), 0, 0],\n", + "[-e**(I*k_y) + e**(-I*k_y), 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}e^{i k_{x}} + e^{- i k_{x}} & 0 & 0 & 0\\\\0 & e^{i k_{y}} + e^{- i k_{y}} & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[e**(I*k_x) + e**(-I*k_x), 0, 0, 0],\n", + "[ 0, e**(I*k_y) + e**(-I*k_y), 0, 0],\n", + "[ 0, 0, 0, 0],\n", + "[ 0, 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & - i e^{i k_{x}} + i e^{- i k_{x}}\\\\0 & 0 & e^{i k_{y}} - e^{- i k_{y}} & 0\\\\0 & - e^{i k_{y}} + e^{- i k_{y}} & 0 & 0\\\\- i e^{i k_{x}} + i e^{- i k_{x}} & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, 0, -I*e**(I*k_x) + I*e**(-I*k_x)],\n", + "[ 0, 0, e**(I*k_y) - e**(-I*k_y), 0],\n", + "[ 0, -e**(I*k_y) + e**(-I*k_y), 0, 0],\n", + "[-I*e**(I*k_x) + I*e**(-I*k_x), 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & e^{i k_{x}} + e^{- i k_{x}} & 0\\\\0 & 0 & 0 & e^{i k_{y}} + e^{- i k_{y}}\\\\e^{i k_{x}} + e^{- i k_{x}} & 0 & 0 & 0\\\\0 & e^{i k_{y}} + e^{- i k_{y}} & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, e**(I*k_x) + e**(-I*k_x), 0],\n", + "[ 0, 0, 0, e**(I*k_y) + e**(-I*k_y)],\n", + "[e**(I*k_x) + e**(-I*k_x), 0, 0, 0],\n", + "[ 0, e**(I*k_y) + e**(-I*k_y), 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & e^{i k_{y}} + e^{- i k_{y}} & 0\\\\0 & 0 & 0 & e^{i k_{x}} + e^{- i k_{x}}\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, 0],\n", + "[0, 0, e**(I*k_y) + e**(-I*k_y), 0],\n", + "[0, 0, 0, e**(I*k_x) + e**(-I*k_x)]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}}\\\\0 & 0 & i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x)],\n", + "[0, 0, I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & e^{i k_{x}} + e^{- i k_{x}} & 0\\\\0 & 0 & 0 & e^{i k_{y}} + e^{- i k_{y}}\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, 0],\n", + "[0, 0, e**(I*k_x) + e**(-I*k_x), 0],\n", + "[0, 0, 0, e**(I*k_y) + e**(-I*k_y)]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & - e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\\\e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}} & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, -e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x), 0, 0],\n", + "[e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x), 0, 0, 0],\n", + "[ 0, 0, 0, 0],\n", + "[ 0, 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & i e^{i k_{y}} + i e^{- i k_{y}} & 0\\\\0 & 0 & 0 & - i e^{i k_{x}} - i e^{- i k_{x}}\\\\- i e^{i k_{y}} - i e^{- i k_{y}} & 0 & 0 & 0\\\\0 & i e^{i k_{x}} + i e^{- i k_{x}} & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, I*e**(I*k_y) + I*e**(-I*k_y), 0],\n", + "[ 0, 0, 0, -I*e**(I*k_x) - I*e**(-I*k_x)],\n", + "[-I*e**(I*k_y) - I*e**(-I*k_y), 0, 0, 0],\n", + "[ 0, I*e**(I*k_x) + I*e**(-I*k_x), 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & i e^{i k_{y}} - i e^{- i k_{y}}\\\\0 & 0 & e^{i k_{x}} - e^{- i k_{x}} & 0\\\\0 & - e^{i k_{x}} + e^{- i k_{x}} & 0 & 0\\\\i e^{i k_{y}} - i e^{- i k_{y}} & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, 0, I*e**(I*k_y) - I*e**(-I*k_y)],\n", + "[ 0, 0, e**(I*k_x) - e**(-I*k_x), 0],\n", + "[ 0, -e**(I*k_x) + e**(-I*k_x), 0, 0],\n", + "[I*e**(I*k_y) - I*e**(-I*k_y), 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & - e^{i k_{x}} + e^{- i k_{x}}\\\\0 & 0 & i e^{i k_{y}} - i e^{- i k_{y}} & 0\\\\0 & i e^{i k_{y}} - i e^{- i k_{y}} & 0 & 0\\\\e^{i k_{x}} - e^{- i k_{x}} & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, 0, -e**(I*k_x) + e**(-I*k_x)],\n", + "[ 0, 0, I*e**(I*k_y) - I*e**(-I*k_y), 0],\n", + "[ 0, I*e**(I*k_y) - I*e**(-I*k_y), 0, 0],\n", + "[e**(I*k_x) - e**(-I*k_x), 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & - i e^{i k_{x}} - i e^{- i k_{x}} & 0\\\\0 & 0 & 0 & i e^{i k_{y}} + i e^{- i k_{y}}\\\\i e^{i k_{x}} + i e^{- i k_{x}} & 0 & 0 & 0\\\\0 & - i e^{i k_{y}} - i e^{- i k_{y}} & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, -I*e**(I*k_x) - I*e**(-I*k_x), 0],\n", + "[ 0, 0, 0, I*e**(I*k_y) + I*e**(-I*k_y)],\n", + "[I*e**(I*k_x) + I*e**(-I*k_x), 0, 0, 0],\n", + "[ 0, -I*e**(I*k_y) - I*e**(-I*k_y), 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & - e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}}\\\\0 & 0 & e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}} & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, -e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x)],\n", + "[0, 0, e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x), 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "familyC4T = qsymm.bloch_family(hopping_vectors, {C4T}, norbs=norbs)\n", + "print(len(familyC4T), [len(fam) for fam in familyC4T])\n", + "qsymm.display_family(familyC4T)" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "6 [4, 4, 4, 4, 4, 4]\n" + ] + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}- e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} - e^{- i k_{x}} & 0 & 0 & 0\\\\0 & e^{i k_{x}} - e^{i k_{y}} - e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[-e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) - e**(-I*k_x), 0, 0, 0],\n", + "[ 0, e**(I*k_x) - e**(I*k_y) - e**(-I*k_y) + e**(-I*k_x), 0, 0],\n", + "[ 0, 0, 0, 0],\n", + "[ 0, 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & - e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} - e^{- i k_{x}} & 0\\\\0 & 0 & 0 & e^{i k_{x}} - e^{i k_{y}} - e^{- i k_{y}} + e^{- i k_{x}}\\\\- e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} - e^{- i k_{x}} & 0 & 0 & 0\\\\0 & e^{i k_{x}} - e^{i k_{y}} - e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, -e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) - e**(-I*k_x), 0],\n", + "[ 0, 0, 0, e**(I*k_x) - e**(I*k_y) - e**(-I*k_y) + e**(-I*k_x)],\n", + "[-e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) - e**(-I*k_x), 0, 0, 0],\n", + "[ 0, e**(I*k_x) - e**(I*k_y) - e**(-I*k_y) + e**(-I*k_x), 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & - i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} + i e^{- i k_{x}}\\\\0 & 0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}} & 0\\\\0 & i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0\\\\- i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} + i e^{- i k_{x}} & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, 0, -I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) + I*e**(-I*k_x)],\n", + "[ 0, 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x), 0],\n", + "[ 0, I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0, 0],\n", + "[-I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) + I*e**(-I*k_x), 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & - e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} - e^{- i k_{x}} & 0\\\\0 & 0 & 0 & e^{i k_{x}} - e^{i k_{y}} - e^{- i k_{y}} + e^{- i k_{x}}\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 0, 0, 0],\n", + "[0, 0, 0, 0],\n", + "[0, 0, -e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) - e**(-I*k_x), 0],\n", + "[0, 0, 0, e**(I*k_x) - e**(I*k_y) - e**(-I*k_y) + e**(-I*k_x)]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & - i e^{i k_{x}} + i e^{i k_{y}} + i e^{- i k_{y}} - i e^{- i k_{x}} & 0\\\\0 & 0 & 0 & - i e^{i k_{x}} + i e^{i k_{y}} + i e^{- i k_{y}} - i e^{- i k_{x}}\\\\i e^{i k_{x}} - i e^{i k_{y}} - i e^{- i k_{y}} + i e^{- i k_{x}} & 0 & 0 & 0\\\\0 & i e^{i k_{x}} - i e^{i k_{y}} - i e^{- i k_{y}} + i e^{- i k_{x}} & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, -I*e**(I*k_x) + I*e**(I*k_y) + I*e**(-I*k_y) - I*e**(-I*k_x), 0],\n", + "[ 0, 0, 0, -I*e**(I*k_x) + I*e**(I*k_y) + I*e**(-I*k_y) - I*e**(-I*k_x)],\n", + "[I*e**(I*k_x) - I*e**(I*k_y) - I*e**(-I*k_y) + I*e**(-I*k_x), 0, 0, 0],\n", + "[ 0, I*e**(I*k_x) - I*e**(I*k_y) - I*e**(-I*k_y) + I*e**(-I*k_x), 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}}\\\\0 & 0 & e^{i k_{x}} - i e^{i k_{y}} + i e^{- i k_{y}} - e^{- i k_{x}} & 0\\\\0 & - e^{i k_{x}} - i e^{i k_{y}} + i e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\\\- e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, 0, 0, e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x)],\n", + "[ 0, 0, e**(I*k_x) - I*e**(I*k_y) + I*e**(-I*k_y) - e**(-I*k_x), 0],\n", + "[ 0, -e**(I*k_x) - I*e**(I*k_y) + I*e**(-I*k_y) + e**(-I*k_x), 0, 0],\n", + "[-e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x), 0, 0, 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "family_new = qsymm.hamiltonian_generator.subtract_family(familyC4T, familyTR, prettify=True)\n", + "print(len(family_new), [len(fam) for fam in family_new])\n", + "qsymm.display_family(family_new)" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(1-0j)" + ] + }, + "execution_count": 37, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from scipy import linalg as la\n", + "la.det(np.kron(np.eye(4), 2 * S[1]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Find periods from BlochModel" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "ham = qsymm.hamiltonian_generator.hamiltonian_from_family(family, tosympy=False)" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{(array([0., 1.]), c0):\n", + "[[0.+0.j 1.+0.j]\n", + " [0.+0.j 0.+0.j]],\n", + "\n", + "(array([-0.8660254, 0.5 ]), c0):\n", + "[[0.+0.j 0.+0.j]\n", + " [1.+0.j 0.+0.j]],\n", + "\n", + "(array([ 0., -1.]), c0):\n", + "[[0.+0.j 0.+0.j]\n", + " [1.+0.j 0.+0.j]],\n", + "\n", + "(array([-0.8660254, -0.5 ]), c0):\n", + "[[0.+0.j 1.+0.j]\n", + " [0.+0.j 0.+0.j]],\n", + "\n", + "(array([0.8660254, 0.5 ]), c0):\n", + "[[0.+0.j 0.+0.j]\n", + " [1.+0.j 0.+0.j]],\n", + "\n", + "(array([ 0.8660254, -0.5 ]), c0):\n", + "[[0.+0.j 1.+0.j]\n", + " [0.+0.j 0.+0.j]],\n", + "\n", + "}" + ] + }, + "execution_count": 49, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ham = qsymm.BlochModel(ham)\n", + "ham" + ] + }, + { + "cell_type": "code", + "execution_count": 96, + "metadata": {}, + "outputs": [], + "source": [ + "from collections import defaultdict\n", + "from qsymm.model import BlochModel, BlochCoeff\n", + "from itertools import product\n", + "\n", + "one = sympy.sympify(1)\n", + "\n", + "def model_periods(model):\n", + " \"\"\"Find the periods and positions of orbitals\n", + " for the lattice model where it is is periodic.\n", + " It needs to be a bloch Hamiltonian written in the\n", + " real space convention (i.e. hopping phases correspond\n", + " to the real space sepatarion of orbitals)\"\"\"\n", + " # convert to BlochModel\n", + " model = BlochModel(model)\n", + " n = model.shape[0]\n", + " d = len(model.momenta)\n", + "\n", + " # dictionary of relative positions for each pair of\n", + " # orbitals, allow multiple relative positions\n", + " rel_pos_nn = defaultdict(lambda: set())\n", + " # populate it with all the hoppings in the model\n", + " for (hop, _), val in model.items():\n", + " for a, b in zip(*np.nonzero(val)):\n", + " rel_pos_nn[a, b].add(BlochCoeff(hop, one))\n", + "\n", + " # keep adding further neighbor relative positions\n", + " # until we find at least d distances to every site\n", + " rel_pos = defaultdict(lambda: set(), {ind: set.copy() for ind, set in rel_pos_nn.items()})\n", + " while any([len(rel_pos[a, b]) < d for a, b in product(range(n), repeat=2)]):\n", + " for ((a1, b1), hops1), ((a2, b2), hops2) in product(rel_pos.items(), rel_pos_nn.items()):\n", + " if b1 == a2:\n", + " for (hop1, _), (hop2, _) in product(hops1, hops2):\n", + " rel_pos[a1, b2].add(BlochCoeff(hop1 + hop2, one))\n", + " return rel_pos" + ] + }, + { + "cell_type": "code", + "execution_count": 98, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 128 ms, sys: 0 ns, total: 128 ms\n", + "Wall time: 121 ms\n" + ] + } + ], + "source": [ + "%%time\n", + "rel_pos = model_periods(ham)\n", + "rel_pos" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[False, False, False, False]" + ] + }, + "execution_count": 78, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "[rel_pos[a, b] == set() for a, b in product(range(2), repeat=2)]" + ] + }, + { + "cell_type": "code", + "execution_count": 79, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "dict_items([((0, 1), {(array([0., 1.]), 1), (array([-0.8660254, -0.5 ]), 1), (array([ 0.8660254, -0.5 ]), 1)}), ((1, 0), {(array([-0.8660254, 0.5 ]), 1), (array([ 0., -1.]), 1), (array([0.8660254, 0.5 ]), 1)}), ((0, 0), {(array([-0.8660254, 1.5 ]), 1), (array([-1.73205081, 0. ]), 1), (array([-0.8660254, -1.5 ]), 1), (array([0., 0.]), 1), (array([0.8660254, 1.5 ]), 1), (array([ 0.8660254, -1.5 ]), 1), (array([1.73205081, 0. ]), 1)}), ((1, 1), {(array([-0.8660254, 1.5 ]), 1), (array([-0.8660254, -1.5 ]), 1), (array([ 0.8660254, -1.5 ]), 1), (array([-1.73205081, 0. ]), 1), (array([0., 0.]), 1), (array([0.8660254, 1.5 ]), 1), (array([1.73205081, 0. ]), 1)})])" + ] + }, + "execution_count": 79, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rel_pos.items()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:.conda-kwant_conda]", + "language": "python", + "name": "conda-env-.conda-kwant_conda-py" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} -- GitLab From 0b2535ee4c18d14f32960e969a27fff46c210819 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Fri, 8 Jan 2021 12:53:20 +0100 Subject: [PATCH 2/4] continue implementation of period finding algorithm --- example_notebooks/bloch_generator.ipynb | 302 ++++++++++++++++++------ 1 file changed, 230 insertions(+), 72 deletions(-) diff --git a/example_notebooks/bloch_generator.ipynb b/example_notebooks/bloch_generator.ipynb index 2c1d7aa..fcbcbb0 100644 --- a/example_notebooks/bloch_generator.ipynb +++ b/example_notebooks/bloch_generator.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -13,7 +13,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -38,7 +38,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -66,7 +66,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -83,7 +83,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -92,7 +92,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "metadata": {}, "outputs": [ { @@ -123,13 +123,8 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, + "execution_count": 104, + "metadata": {}, "outputs": [], "source": [ "C3 = qsymm.rotation(1/3, U=np.eye(2))\n", @@ -138,13 +133,8 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, + "execution_count": 105, + "metadata": {}, "outputs": [], "source": [ "norbs = OrderedDict([('A', 1), ('B', 1)]) # A and B atom per unit cell, one orbital each\n", @@ -153,13 +143,8 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, + "execution_count": 106, + "metadata": {}, "outputs": [], "source": [ "family = qsymm.bloch_family(hopping_vectors, symmetries, norbs, bloch_model=True)" @@ -167,14 +152,24 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true + "execution_count": 107, + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & e^{\\frac{i k_{x}}{2}} e^{- \\frac{\\sqrt{3} i k_{y}}{6}} + e^{\\frac{\\sqrt{3} i k_{y}}{3}} + e^{- \\frac{i k_{x}}{2}} e^{- \\frac{\\sqrt{3} i k_{y}}{6}}\\\\e^{\\frac{i k_{x}}{2}} e^{\\frac{\\sqrt{3} i k_{y}}{6}} + e^{- \\frac{\\sqrt{3} i k_{y}}{3}} + e^{- \\frac{i k_{x}}{2}} e^{\\frac{\\sqrt{3} i k_{y}}{6}} & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, e**(I*k_x/2)*e**(-sqrt(3)*I*k_y/6) + e**(sqrt(3)*I*k_y/3) + e**(-I*k_x/2)*e**(-sqrt(3)*I*k_y/6)],\n", + "[e**(I*k_x/2)*e**(sqrt(3)*I*k_y/6) + e**(-sqrt(3)*I*k_y/3) + e**(-I*k_x/2)*e**(sqrt(3)*I*k_y/6), 0]])" + ] + }, + "metadata": {}, + "output_type": "display_data" } - }, - "outputs": [], + ], "source": [ "qsymm.display_family(family)" ] @@ -1352,7 +1347,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 108, "metadata": {}, "outputs": [], "source": [ @@ -1361,40 +1356,40 @@ }, { "cell_type": "code", - "execution_count": 49, + "execution_count": 109, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "{(array([0., 1.]), c0):\n", - "[[0.+0.j 1.+0.j]\n", - " [0.+0.j 0.+0.j]],\n", - "\n", - "(array([-0.8660254, 0.5 ]), c0):\n", + "{(array([0.5 , 0.28867513]), c0):\n", "[[0.+0.j 0.+0.j]\n", " [1.+0.j 0.+0.j]],\n", "\n", - "(array([ 0., -1.]), c0):\n", - "[[0.+0.j 0.+0.j]\n", - " [1.+0.j 0.+0.j]],\n", + "(array([ 0.5 , -0.28867513]), c0):\n", + "[[0.+0.j 1.+0.j]\n", + " [0.+0.j 0.+0.j]],\n", + "\n", + "(array([-0.5 , -0.28867513]), c0):\n", + "[[0.+0.j 1.+0.j]\n", + " [0.+0.j 0.+0.j]],\n", "\n", - "(array([-0.8660254, -0.5 ]), c0):\n", + "(array([-0. , 0.57735027]), c0):\n", "[[0.+0.j 1.+0.j]\n", " [0.+0.j 0.+0.j]],\n", "\n", - "(array([0.8660254, 0.5 ]), c0):\n", + "(array([-0.5 , 0.28867513]), c0):\n", "[[0.+0.j 0.+0.j]\n", " [1.+0.j 0.+0.j]],\n", "\n", - "(array([ 0.8660254, -0.5 ]), c0):\n", - "[[0.+0.j 1.+0.j]\n", - " [0.+0.j 0.+0.j]],\n", + "(array([-0. , -0.57735027]), c0):\n", + "[[0.+0.j 0.+0.j]\n", + " [1.+0.j 0.+0.j]],\n", "\n", "}" ] }, - "execution_count": 49, + "execution_count": 109, "metadata": {}, "output_type": "execute_result" } @@ -1406,13 +1401,18 @@ }, { "cell_type": "code", - "execution_count": 96, + "execution_count": 148, "metadata": {}, "outputs": [], "source": [ "from collections import defaultdict\n", "from qsymm.model import BlochModel, BlochCoeff\n", + "from qsymm.kwant_linalg_lll import lll\n", "from itertools import product\n", + "from functools import reduce\n", + "from math import gcd\n", + "# from https://github.com/lan496/hsnf\n", + "from hsnf import smith_normal_form, row_style_hermite_normal_form\n", "\n", "one = sympy.sympify(1)\n", "\n", @@ -1435,75 +1435,233 @@ " for a, b in zip(*np.nonzero(val)):\n", " rel_pos_nn[a, b].add(BlochCoeff(hop, one))\n", "\n", - " # keep adding further neighbor relative positions\n", - " # until we find at least d distances to every site\n", + " # Keep adding further neighbor relative positions up to n steps,\n", + " # this guarantees we find all linearly independent (over integers)\n", + " # lattice vectors.\n", + " # This may take long if there are many orbitals, perhaps there is\n", + " # a better stopping condition.\n", " rel_pos = defaultdict(lambda: set(), {ind: set.copy() for ind, set in rel_pos_nn.items()})\n", - " while any([len(rel_pos[a, b]) < d for a, b in product(range(n), repeat=2)]):\n", + " for _ in range(n-1):\n", " for ((a1, b1), hops1), ((a2, b2), hops2) in product(rel_pos.items(), rel_pos_nn.items()):\n", " if b1 == a2:\n", - " for (hop1, _), (hop2, _) in product(hops1, hops2):\n", - " rel_pos[a1, b2].add(BlochCoeff(hop1 + hop2, one))\n", - " return rel_pos" + " for (hop1, hop2) in product(hops1, hops2):\n", + " rel_pos[a1, b2].add(hop1 * hop2)\n", + "\n", + " # Lattice vectors are those connecting the same site type\n", + " vecs = np.array([vec for ((a, b), hops) in rel_pos.items() if a == b\n", + " for vec, _ in hops])\n", + " # Find primitive vectors\n", + " prim_vecs = primitive_lattice_vecs(vecs)\n", + " return prim_vecs\n", + "\n", + "\n", + "def make_int(R):\n", + " # If close to an integer array convert to integer array, else\n", + " # return None\n", + " R_int = np.array(np.round(R), int)\n", + " if qsymm.linalg.allclose(R, R_int):\n", + " return R_int\n", + " else:\n", + " return None\n", + "\n", + "\n", + "def primitive_lattice_vecs(vecs, method='hnf'):\n", + " \"\"\"Find a set of primitive lattice vectors spanning the\n", + " same lattice as vecs.\"\"\"\n", + " vecs = np.asarray(vecs)\n", + " d = vecs.shape[1]\n", + " # pick out a set of small linearly independent vectors\n", + " ### TODO this likely works even if vecs don't span the full space\n", + " norms = np.linalg.norm(vecs, axis=1)\n", + " vecs = vecs[np.argsort(norms)]\n", + " bas = []\n", + " for v in vecs:\n", + " new_bas = bas + [v]\n", + " rank = np.linalg.matrix_rank(new_bas)\n", + " if rank == len(new_bas):\n", + " bas = new_bas\n", + " if rank == d:\n", + " break\n", + " else:\n", + " raise ValueError('The vectors do not span the full space.')\n", + " bas = np.vstack(bas)\n", + " # rewrite all vectors in this new basis\n", + " vecs = np.linalg.solve(bas.T, vecs.T).T\n", + " # find the common denominator to make all integer\n", + " denom = lcd(vecs)\n", + " vecs = make_int(vecs * denom)\n", + " bas = bas / denom\n", + " # find primitive lattice vectors\n", + " if method == 'snf':\n", + " D, L, R = smith_normal_form(vecs)\n", + " prim_vecs = D @ np.linalg.inv(R)\n", + " elif method == 'hnf':\n", + " prim_vecs, U = row_style_hermite_normal_form(vecs)\n", + " else:\n", + " raise ValueError(\"Method must be either 'hnf' or 'snf'.\")\n", + " prim_vecs = prim_vecs[:d]\n", + " # convert back to original basis\n", + " prim_vecs = prim_vecs @ bas\n", + " return lll(prim_vecs)[0]\n", + "\n", + "\n", + "def lcd(a, tol=1e-9):\n", + " \"\"\"Calculate approximate least common denominator for\n", + " the entries of floating point array `a`.\"\"\"\n", + " a = np.asarray(a).flatten()\n", + " a = a % 1\n", + " ### This would be nice to vectorize\n", + " denoms = np.array([farey(x)[1] for x in a])\n", + " return reduce(lambda a, b: a * b // gcd(a, b), denoms, 1)\n", + " \n", + "\n", + "def farey(x, tol=1e-9):\n", + " \"\"\"Farey sequence rational approximation of `0 <= x <= 1`\n", + " with tolerance `tol`. Adapted from\n", + " https://www.johndcook.com/blog/2010/10/20/best-rational-approximation/\n", + " \"\"\"\n", + " N = int(tol**(-1/2) / np.sqrt(5))\n", + " a, b = 0, 1\n", + " c, d = 1, 1\n", + " while (b <= N and d <= N):\n", + " # Add these checks to avoid additional loops\n", + " if abs(x - a/b) <= tol:\n", + " return a, b\n", + " elif abs(x - c/d) <= tol:\n", + " return c, d\n", + " mediant = (a+c)/(b+d)\n", + " if abs(x - mediant) <= tol:\n", + " return a+c, b+d\n", + " elif x > mediant:\n", + " a, b = a+c, b+d\n", + " else:\n", + " c, d = a+c, b+d\n", + "\n", + " if (b > N):\n", + " return c, d\n", + " else:\n", + " return a, b" ] }, { "cell_type": "code", - "execution_count": 98, + "execution_count": 149, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "CPU times: user 128 ms, sys: 0 ns, total: 128 ms\n", - "Wall time: 121 ms\n" + "55.8 ms ± 4.06 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" ] + }, + { + "data": { + "text/plain": [ + "array([[ 1.00000000e+00, -2.77555756e-16],\n", + " [ 5.00000000e-01, 8.66025404e-01]])" + ] + }, + "execution_count": 149, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "%%time\n", - "rel_pos = model_periods(ham)\n", - "rel_pos" + "# %%time\n", + "%timeit lat_vecs = model_periods(ham)\n", + "lat_vecs" + ] + }, + { + "cell_type": "code", + "execution_count": 114, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext line_profiler" ] }, { "cell_type": "code", - "execution_count": 78, + "execution_count": 150, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "[False, False, False, False]" + "Timer unit: 1e-06 s\n", + "\n", + "Total time: 0.000181 s\n", + "File: \n", + "Function: farey at line 111\n", + "\n", + "Line # Hits Time Per Hit % Time Line Contents\n", + "==============================================================\n", + " 111 def farey(x, tol=1e-9):\n", + " 112 \"\"\"Farey sequence rational approximation of `0 <= x <= 1`\n", + " 113 with tolerance `tol`. Adapted from\n", + " 114 https://www.johndcook.com/blog/2010/10/20/best-rational-approximation/\n", + " 115 \"\"\"\n", + " 116 28 81.0 2.9 44.8 N = int(tol**(-1/2) / np.sqrt(5))\n", + " 117 28 13.0 0.5 7.2 a, b = 0, 1\n", + " 118 28 11.0 0.4 6.1 c, d = 1, 1\n", + " 119 28 13.0 0.5 7.2 while (b <= N and d <= N):\n", + " 120 # Add these checks to avoid additional loops\n", + " 121 28 40.0 1.4 22.1 if abs(x - a/b) < tol:\n", + " 122 24 16.0 0.7 8.8 return a, b\n", + " 123 4 4.0 1.0 2.2 elif abs(x - c/d) < tol:\n", + " 124 4 3.0 0.8 1.7 return c, d\n", + " 125 mediant = (a+c)/(b+d)\n", + " 126 if x == mediant:\n", + " 127 return a+c, b+d\n", + " 128 elif x > mediant:\n", + " 129 a, b = a+c, b+d\n", + " 130 else:\n", + " 131 c, d = a+c, b+d\n", + " 132 \n", + " 133 if (b > N):\n", + " 134 return c, d\n", + " 135 else:\n", + " 136 return a, b" ] }, - "execution_count": 78, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "[rel_pos[a, b] == set() for a, b in product(range(2), repeat=2)]" + "%lprun -f farey model_periods(ham)" ] }, { "cell_type": "code", - "execution_count": 79, + "execution_count": 88, "metadata": {}, "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 142 ms, sys: 24 µs, total: 142 ms\n", + "Wall time: 138 ms\n" + ] + }, { "data": { "text/plain": [ - "dict_items([((0, 1), {(array([0., 1.]), 1), (array([-0.8660254, -0.5 ]), 1), (array([ 0.8660254, -0.5 ]), 1)}), ((1, 0), {(array([-0.8660254, 0.5 ]), 1), (array([ 0., -1.]), 1), (array([0.8660254, 0.5 ]), 1)}), ((0, 0), {(array([-0.8660254, 1.5 ]), 1), (array([-1.73205081, 0. ]), 1), (array([-0.8660254, -1.5 ]), 1), (array([0., 0.]), 1), (array([0.8660254, 1.5 ]), 1), (array([ 0.8660254, -1.5 ]), 1), (array([1.73205081, 0. ]), 1)}), ((1, 1), {(array([-0.8660254, 1.5 ]), 1), (array([-0.8660254, -1.5 ]), 1), (array([ 0.8660254, -1.5 ]), 1), (array([-1.73205081, 0. ]), 1), (array([0., 0.]), 1), (array([0.8660254, 1.5 ]), 1), (array([1.73205081, 0. ]), 1)})])" + "array([[1., 0.],\n", + " [0., 1.]])" ] }, - "execution_count": 79, + "execution_count": 88, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "rel_pos.items()" + "%%time\n", + "primitive_lattice_vecs([[11, 1], [1, 0], [0, 0]], method='hnf')" ] }, { -- GitLab From a983184cd152f63f54289cd7fabed8bae2d357f4 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Fri, 8 Jan 2021 17:36:23 +0100 Subject: [PATCH 3/4] finish implementation of automatic candidate finding --- example_notebooks/bloch_generator.ipynb | 1696 ----------------------- qsymm/linalg.py | 40 + qsymm/symmetry_finder.py | 153 +- qsymm/tests/test_symmetry_finder.py | 7 +- 4 files changed, 192 insertions(+), 1704 deletions(-) delete mode 100644 example_notebooks/bloch_generator.ipynb diff --git a/example_notebooks/bloch_generator.ipynb b/example_notebooks/bloch_generator.ipynb deleted file mode 100644 index fcbcbb0..0000000 --- a/example_notebooks/bloch_generator.ipynb +++ /dev/null @@ -1,1696 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import sympy\n", - "from collections import OrderedDict" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "import qsymm" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Graphene" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We use the Hamiltonian generator to reproduce the spinless nearest neighbour tight binding Hamiltonian for graphene.\n", - "\n", - "The generators of the symmetry group are time-reversal symmetry, sublattice (or chiral) symmetry, and threefold rotation symmetry." - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "# Time reversal\n", - "TR = qsymm.time_reversal(2, U=np.eye(2))\n", - "\n", - "# Chiral symmetry\n", - "C = qsymm.chiral(2, U=np.array([[1, 0], [0, -1]]))\n", - "\n", - "# Atom A rotates into A, B into B, use exact sympy representation\n", - "sphi = 2*sympy.pi/3\n", - "RC3 = sympy.Matrix([[sympy.cos(sphi), -sympy.sin(sphi)],\n", - " [sympy.sin(sphi), sympy.cos(sphi)]])\n", - "C3 = qsymm.PointGroupElement(RC3, False, False, np.eye(2))\n", - "\n", - "symmetries = [C, TR, C3]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "There are two carbon atoms per unit cell (A and B) with one orbital each. The lattice is triangular, and only includes hoppings between nearest neighbour atoms. This restricts hoppings to only those between atoms of different types, such that each atom couples to three neighbouring atoms. Using the symmetrization strategy to generate the Hamiltonian, it is sufficient to specify hoppings to one such neighbour along with the symmetry generators, and we take the vector $(1, 0)$ to connect this neighbouring pair of atoms." - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "norbs = OrderedDict([('A', 1), ('B', 1)]) # A and B atom per unit cell, one orbital each\n", - "hopping_vectors = [('A', 'B', [0, 1])] # Hopping between neighbouring A and B atoms" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Generate the Hamiltonian." - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "family = qsymm.bloch_family(hopping_vectors, symmetries, norbs)" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [ - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & e^{i k_{y}} + e^{- \\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} + e^{- \\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}}\\\\e^{\\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} + e^{\\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}} + e^{- i k_{y}} & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, e**(I*k_y) + e**(-I*k_y/2)*e**(sqrt(3)*I*k_x/2) + e**(-I*k_y/2)*e**(-sqrt(3)*I*k_x/2)],\n", - "[e**(I*k_y/2)*e**(sqrt(3)*I*k_x/2) + e**(I*k_y/2)*e**(-sqrt(3)*I*k_x/2) + e**(-I*k_y), 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "qsymm.display_family(family)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Scale the bond length in terms of the graphene lattice constant, and have the function return a list of BlochModel objects. For this we can use a more convenient definition of the rotation" - ] - }, - { - "cell_type": "code", - "execution_count": 104, - "metadata": {}, - "outputs": [], - "source": [ - "C3 = qsymm.rotation(1/3, U=np.eye(2))\n", - "symmetries = [C, TR, C3]" - ] - }, - { - "cell_type": "code", - "execution_count": 105, - "metadata": {}, - "outputs": [], - "source": [ - "norbs = OrderedDict([('A', 1), ('B', 1)]) # A and B atom per unit cell, one orbital each\n", - "hopping_vectors = [('A', 'B', [0, 1/np.sqrt(3)])] # Hopping between neighbouring A and B atoms" - ] - }, - { - "cell_type": "code", - "execution_count": 106, - "metadata": {}, - "outputs": [], - "source": [ - "family = qsymm.bloch_family(hopping_vectors, symmetries, norbs, bloch_model=True)" - ] - }, - { - "cell_type": "code", - "execution_count": 107, - "metadata": {}, - "outputs": [ - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & e^{\\frac{i k_{x}}{2}} e^{- \\frac{\\sqrt{3} i k_{y}}{6}} + e^{\\frac{\\sqrt{3} i k_{y}}{3}} + e^{- \\frac{i k_{x}}{2}} e^{- \\frac{\\sqrt{3} i k_{y}}{6}}\\\\e^{\\frac{i k_{x}}{2}} e^{\\frac{\\sqrt{3} i k_{y}}{6}} + e^{- \\frac{\\sqrt{3} i k_{y}}{3}} + e^{- \\frac{i k_{x}}{2}} e^{\\frac{\\sqrt{3} i k_{y}}{6}} & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, e**(I*k_x/2)*e**(-sqrt(3)*I*k_y/6) + e**(sqrt(3)*I*k_y/3) + e**(-I*k_x/2)*e**(-sqrt(3)*I*k_y/6)],\n", - "[e**(I*k_x/2)*e**(sqrt(3)*I*k_y/6) + e**(-sqrt(3)*I*k_y/3) + e**(-I*k_x/2)*e**(sqrt(3)*I*k_y/6), 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "qsymm.display_family(family)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Three-orbital tight binding model for monolayer MX$_2$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We use the Hamiltonian generator to reproduce the tight binding model for monolayer MX$_2$ published in Phys. Rev. B 88, 085433 (2013).\n", - "\n", - "The generators of the symmetry group of the tight binding model are time reversal symmetry, mirror symmetry and threefold rotation symmetry." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [], - "source": [ - "# Time reversal\n", - "TR = qsymm.time_reversal(2, np.eye(3))\n", - "\n", - "# Mirror symmetry \n", - "Mx = qsymm.mirror([1, 0], np.diag([1, -1, 1]))\n", - "\n", - "# Threefold rotation on d_z^2, d_xy, d_x^2-y^2 states.\n", - "C3U = np.array([[1, 0, 0],\n", - " [0, -0.5, -np.sqrt(3)/2],\n", - " [0, np.sqrt(3)/2, -0.5]])\n", - "# Could also use the predefined representation of rotations on d-orbitals\n", - "Ld = qsymm.groups.L_matrices(3, 2)\n", - "C3U2 = qsymm.groups.spin_rotation(2 * np.pi * np.array([0, 0, 1/3]), Ld)\n", - "# Restrict to d_z^2, d_xy, d_x^2-y^2 states\n", - "mask = np.array([1, 2 ,0])\n", - "C3U2 = C3U2[mask][:, mask]\n", - "assert np.allclose(C3U, C3U2)\n", - "\n", - "C3 = qsymm.rotation(1/3, U=C3U)\n", - "\n", - "symmetries = [TR, Mx, C3]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Next, we specify the hoppings to include. The tight binding model has a triangular lattice, three orbitals per M atom, and nearest neighbour hopping." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [], - "source": [ - "# One site per unit cell (M atom), with three orbitals\n", - "norbs = OrderedDict([('a', 3)])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Each atom has six nearest neighbour atoms at a distance of one primitive lattice vector. Since we use the symmetrization strategy to generate the Hamiltonian, it is sufficient to specify a hopping to one nearest neighbour atom along with the symmetry generators. We take the primitive vector connecting the pair of atoms to be $(1, 0)$." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [], - "source": [ - "# Hopping to a neighbouring atom one primitive lattice vector away\n", - "hopping_vectors = [('a', 'a', [1, 0])]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Generate the tight binding Hamiltonian." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [], - "source": [ - "family = qsymm.bloch_family(hopping_vectors, symmetries, norbs, bloch_model=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The Hamiltonian family should include 8 linearly independent components, including the onsite terms." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [], - "source": [ - "len(family)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [], - "source": [ - "qsymm.display_family(family, nsimplify=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# 4-site model for monolayer WTe$_2$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We use the Hamiltonian generator to reproduce the tight binding model for monolayer WTe$_2$ published in Phys. Rev. X 6, 041069 (2016).\n", - "\n", - "The generators of the symmetry group of the tight binding model are time reversal symmetry, glide reflection and inversion symmetry." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [], - "source": [ - "# Define 4 sites with one orbital each\n", - "sites = ['Ad', 'Ap', 'Bd', 'Bp']\n", - "norbs = OrderedDict([(site, 1) for site in sites])\n", - "\n", - "# Define symbolic coordinates for orbitals\n", - "rAp = qsymm.sympify('[x_Ap, y_Ap]')\n", - "rAd = qsymm.sympify('[x_Ad, y_Ad]')\n", - "rBp = qsymm.sympify('[x_Bp, y_Bp]')\n", - "rBd = qsymm.sympify('[x_Bd, y_Bd]')\n", - "\n", - "# Define hoppings to include\n", - "hopping_vectors = [('Bd', 'Bd', np.array([1, 0])),\n", - " ('Ap', 'Ap', np.array([1, 0])),\n", - " ('Bd', 'Ap', rAp - rBd),\n", - " ('Ap', 'Bp', rBp - rAp),\n", - " ('Ad', 'Bd', rBd - rAd), \n", - " ]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [], - "source": [ - "# Inversion\n", - "perm_inv = {'Ad': 'Bd', 'Ap': 'Bp', 'Bd': 'Ad', 'Bp': 'Ap'}\n", - "onsite_inv = {site: (1 if site in ['Ad', 'Bd'] else -1) * np.eye(1) for site in sites}\n", - "inversion = qsymm.groups.symmetry_from_permutation(-np.eye(2), perm_inv, norbs, onsite_inv)\n", - "\n", - "# Glide\n", - "perm_glide = {site: site for site in sites}\n", - "onsite_glide = {site: (1 if site in ['Ad', 'Bd'] else -1) * np.eye(1) for site in sites}\n", - "glide = qsymm.groups.symmetry_from_permutation(np.array([[-1, 0],[0, 1]]), perm_glide, norbs, onsite_glide)\n", - "\n", - "# TR\n", - "time_reversal = qsymm.time_reversal(2, np.eye(4))\n", - "\n", - "gens = {glide, inversion, time_reversal}\n", - "sg = qsymm.groups.generate_group(gens)\n", - "len(sg)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Generate the tight binding Hamiltonian." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [], - "source": [ - "family = qsymm.bloch_family(hopping_vectors, gens, norbs=norbs)\n", - "print(len(family), [len(fam) for fam in family])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [], - "source": [ - "qsymm.display_family(family)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Square lattice with 4 sites in the UC" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Make model with square lattice that has 4 sites in the unit cell related by 4-fold rotation. Sites have spin-1/2 and we add time reversal and particle-hole symmetry." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [], - "source": [ - "hopping_vectors = [('a', 'b', np.array([1, 0])),\n", - " ('b', 'a', np.array([1, 0])),\n", - " ('c', 'd', np.array([1, 0])),\n", - " ('d', 'c', np.array([1, 0])),\n", - " ('a', 'c', np.array([0, 1])),\n", - " ('c', 'a', np.array([0, 1])),\n", - " ('b', 'd', np.array([0, 1])),\n", - " ('d', 'b', np.array([0, 1])), \n", - " ]\n", - "\n", - "# Define spin-1/2 operators\n", - "S = qsymm.groups.spin_matrices(1/2)\n", - "# Define real space rotation generators in 2D\n", - "L = qsymm.groups.L_matrices(d=2)\n", - "\n", - "sites = ['a', 'b', 'c', 'd']\n", - "norbs = OrderedDict([(site, 2) for site in sites])\n", - "\n", - "perm_C4 = {'a': 'b', 'b': 'd', 'd': 'c', 'c': 'a'}\n", - "onsite_C4 = {site: qsymm.groups.spin_rotation(2*np.pi * np.array([0, 0, 1/4]), S) for site in sites}\n", - "C4 = qsymm.groups.symmetry_from_permutation(\n", - " qsymm.groups.spin_rotation(2*np.pi * np.array([1/4]), L, roundint=True), \n", - " perm_C4, norbs, onsite_C4)\n", - "\n", - "# Fermionic TR\n", - "time_reversal = qsymm.time_reversal(2, \n", - " np.kron(np.eye(4), qsymm.groups.spin_rotation(2*np.pi * np.array([0, 1/2, 0]), S)))\n", - "\n", - "# define strange PH symmetry\n", - "particle_hole = qsymm.particle_hole(2, np.eye(8))\n", - "\n", - "gens = {C4, time_reversal, particle_hole}\n", - "\n", - "sg = qsymm.groups.generate_group(gens)\n", - "len(sg)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [], - "source": [ - "family = qsymm.bloch_family(hopping_vectors, gens, norbs=norbs)\n", - "print(len(family), [len(fam) for fam in family])\n", - "qsymm.display_family(family)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### C4T symmetric model" - ] - }, - { - "cell_type": "code", - "execution_count": 27, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [ - { - "data": { - "text/plain": [ - "8" - ] - }, - "execution_count": 27, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "hopping_vectors = [('a', 'a', np.array([1, 0])), \n", - " ]\n", - "\n", - "# Define spin-1/2 operators\n", - "S = qsymm.groups.spin_matrices(1/2)\n", - "# Define real space rotation generators in 2D\n", - "L = qsymm.groups.L_matrices(d=2)\n", - "\n", - "sites = ['a']\n", - "norbs = OrderedDict([(site, 4) for site in sites])\n", - "\n", - "perm_C4 = {'a': 'a'}\n", - "onsite_C4 = {site: np.kron(np.eye(2),\n", - " qsymm.groups.spin_rotation(2*np.pi * np.array([0, 0, 1/4]), S)\n", - " )\n", - " for site in sites}\n", - "C4 = qsymm.groups.symmetry_from_permutation(\n", - " qsymm.groups.spin_rotation(2*np.pi * np.array([1/4]), L, roundint=True), \n", - " perm_C4, norbs, onsite_C4)\n", - "\n", - "# Fermionic TR\n", - "T = qsymm.time_reversal(2, \n", - " np.kron(np.eye(2),\n", - " qsymm.groups.spin_rotation(2*np.pi * np.array([0, 1/2, 0]), S),\n", - " )\n", - " )\n", - "\n", - "C4T = C4 * T\n", - "\n", - "sg = qsymm.groups.generate_group(gens)\n", - "len(sg)" - ] - }, - { - "cell_type": "code", - "execution_count": 28, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [ - { - "data": { - "text/plain": [ - "{1, R(π), R(π/2), R(-π/2), T, R(π) T, R(π/2) T, R(-π/2) T}" - ] - }, - "execution_count": 28, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "sg" - ] - }, - { - "cell_type": "code", - "execution_count": 30, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "14 [1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]\n" - ] - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0 & 0\\\\0 & 1 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[1, 0, 0, 0],\n", - "[0, 1, 0, 0],\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 1 & 0\\\\0 & 0 & 0 & 1\\\\1 & 0 & 0 & 0\\\\0 & 1 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[0, 0, 1, 0],\n", - "[0, 0, 0, 1],\n", - "[1, 0, 0, 0],\n", - "[0, 1, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 1 & 0\\\\0 & 0 & 0 & 1\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, 0],\n", - "[0, 0, 1, 0],\n", - "[0, 0, 0, 1]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & i & 0\\\\0 & 0 & 0 & - i\\\\- i & 0 & 0 & 0\\\\0 & i & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, I, 0],\n", - "[ 0, 0, 0, -I],\n", - "[-I, 0, 0, 0],\n", - "[ 0, I, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0 & 0\\\\0 & e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x), 0, 0, 0],\n", - "[ 0, e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x), 0, 0],\n", - "[ 0, 0, 0, 0],\n", - "[ 0, 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0\\\\i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x), 0, 0],\n", - "[I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0, 0, 0],\n", - "[ 0, 0, 0, 0],\n", - "[ 0, 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}} & 0\\\\0 & 0 & 0 & e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}}\\\\e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0 & 0\\\\0 & e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x), 0],\n", - "[ 0, 0, 0, e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x)],\n", - "[e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x), 0, 0, 0],\n", - "[ 0, e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x), 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}}\\\\0 & 0 & i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0\\\\0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0\\\\i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x)],\n", - "[ 0, 0, I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0],\n", - "[ 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x), 0, 0],\n", - "[I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}} & 0\\\\0 & 0 & 0 & e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} + e^{- i k_{x}}\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, 0],\n", - "[0, 0, e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x), 0],\n", - "[0, 0, 0, e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) + e**(-I*k_x)]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}}\\\\0 & 0 & i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x)],\n", - "[0, 0, I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & - e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\\\e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}} & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, -e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x), 0, 0],\n", - "[e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x), 0, 0, 0],\n", - "[ 0, 0, 0, 0],\n", - "[ 0, 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & i e^{i k_{x}} + i e^{i k_{y}} + i e^{- i k_{y}} + i e^{- i k_{x}} & 0\\\\0 & 0 & 0 & - i e^{i k_{x}} - i e^{i k_{y}} - i e^{- i k_{y}} - i e^{- i k_{x}}\\\\- i e^{i k_{x}} - i e^{i k_{y}} - i e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0 & 0\\\\0 & i e^{i k_{x}} + i e^{i k_{y}} + i e^{- i k_{y}} + i e^{- i k_{x}} & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, I*e**(I*k_x) + I*e**(I*k_y) + I*e**(-I*k_y) + I*e**(-I*k_x), 0],\n", - "[ 0, 0, 0, -I*e**(I*k_x) - I*e**(I*k_y) - I*e**(-I*k_y) - I*e**(-I*k_x)],\n", - "[-I*e**(I*k_x) - I*e**(I*k_y) - I*e**(-I*k_y) - I*e**(-I*k_x), 0, 0, 0],\n", - "[ 0, I*e**(I*k_x) + I*e**(I*k_y) + I*e**(-I*k_y) + I*e**(-I*k_x), 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & - e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}}\\\\0 & 0 & e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}} & 0\\\\0 & - e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\\\e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}} & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, 0, -e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x)],\n", - "[ 0, 0, e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x), 0],\n", - "[ 0, -e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x), 0, 0],\n", - "[e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x), 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & - e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}}\\\\0 & 0 & e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}} & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, -e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x)],\n", - "[0, 0, e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x), 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "familyTR = qsymm.bloch_family(hopping_vectors, {T, C4}, norbs=norbs)\n", - "print(len(familyTR), [len(fam) for fam in familyTR])\n", - "qsymm.display_family(familyTR)" - ] - }, - { - "cell_type": "code", - "execution_count": 31, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "20 [1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]\n" - ] - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0 & 0\\\\0 & 1 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[1, 0, 0, 0],\n", - "[0, 1, 0, 0],\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 1 & 0\\\\0 & 0 & 0 & 1\\\\1 & 0 & 0 & 0\\\\0 & 1 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[0, 0, 1, 0],\n", - "[0, 0, 0, 1],\n", - "[1, 0, 0, 0],\n", - "[0, 1, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 1 & 0\\\\0 & 0 & 0 & 1\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, 0],\n", - "[0, 0, 1, 0],\n", - "[0, 0, 0, 1]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & i & 0\\\\0 & 0 & 0 & - i\\\\- i & 0 & 0 & 0\\\\0 & i & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, I, 0],\n", - "[ 0, 0, 0, -I],\n", - "[-I, 0, 0, 0],\n", - "[ 0, I, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}e^{i k_{y}} + e^{- i k_{y}} & 0 & 0 & 0\\\\0 & e^{i k_{x}} + e^{- i k_{x}} & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[e**(I*k_y) + e**(-I*k_y), 0, 0, 0],\n", - "[ 0, e**(I*k_x) + e**(-I*k_x), 0, 0],\n", - "[ 0, 0, 0, 0],\n", - "[ 0, 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0\\\\i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x), 0, 0],\n", - "[I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0, 0, 0],\n", - "[ 0, 0, 0, 0],\n", - "[ 0, 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & e^{i k_{y}} + e^{- i k_{y}} & 0\\\\0 & 0 & 0 & e^{i k_{x}} + e^{- i k_{x}}\\\\e^{i k_{y}} + e^{- i k_{y}} & 0 & 0 & 0\\\\0 & e^{i k_{x}} + e^{- i k_{x}} & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, e**(I*k_y) + e**(-I*k_y), 0],\n", - "[ 0, 0, 0, e**(I*k_x) + e**(-I*k_x)],\n", - "[e**(I*k_y) + e**(-I*k_y), 0, 0, 0],\n", - "[ 0, e**(I*k_x) + e**(-I*k_x), 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & e^{i k_{y}} - e^{- i k_{y}}\\\\0 & 0 & i e^{i k_{x}} - i e^{- i k_{x}} & 0\\\\0 & i e^{i k_{x}} - i e^{- i k_{x}} & 0 & 0\\\\- e^{i k_{y}} + e^{- i k_{y}} & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, 0, e**(I*k_y) - e**(-I*k_y)],\n", - "[ 0, 0, I*e**(I*k_x) - I*e**(-I*k_x), 0],\n", - "[ 0, I*e**(I*k_x) - I*e**(-I*k_x), 0, 0],\n", - "[-e**(I*k_y) + e**(-I*k_y), 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}e^{i k_{x}} + e^{- i k_{x}} & 0 & 0 & 0\\\\0 & e^{i k_{y}} + e^{- i k_{y}} & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[e**(I*k_x) + e**(-I*k_x), 0, 0, 0],\n", - "[ 0, e**(I*k_y) + e**(-I*k_y), 0, 0],\n", - "[ 0, 0, 0, 0],\n", - "[ 0, 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & - i e^{i k_{x}} + i e^{- i k_{x}}\\\\0 & 0 & e^{i k_{y}} - e^{- i k_{y}} & 0\\\\0 & - e^{i k_{y}} + e^{- i k_{y}} & 0 & 0\\\\- i e^{i k_{x}} + i e^{- i k_{x}} & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, 0, -I*e**(I*k_x) + I*e**(-I*k_x)],\n", - "[ 0, 0, e**(I*k_y) - e**(-I*k_y), 0],\n", - "[ 0, -e**(I*k_y) + e**(-I*k_y), 0, 0],\n", - "[-I*e**(I*k_x) + I*e**(-I*k_x), 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & e^{i k_{x}} + e^{- i k_{x}} & 0\\\\0 & 0 & 0 & e^{i k_{y}} + e^{- i k_{y}}\\\\e^{i k_{x}} + e^{- i k_{x}} & 0 & 0 & 0\\\\0 & e^{i k_{y}} + e^{- i k_{y}} & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, e**(I*k_x) + e**(-I*k_x), 0],\n", - "[ 0, 0, 0, e**(I*k_y) + e**(-I*k_y)],\n", - "[e**(I*k_x) + e**(-I*k_x), 0, 0, 0],\n", - "[ 0, e**(I*k_y) + e**(-I*k_y), 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & e^{i k_{y}} + e^{- i k_{y}} & 0\\\\0 & 0 & 0 & e^{i k_{x}} + e^{- i k_{x}}\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, 0],\n", - "[0, 0, e**(I*k_y) + e**(-I*k_y), 0],\n", - "[0, 0, 0, e**(I*k_x) + e**(-I*k_x)]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}}\\\\0 & 0 & i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x)],\n", - "[0, 0, I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & e^{i k_{x}} + e^{- i k_{x}} & 0\\\\0 & 0 & 0 & e^{i k_{y}} + e^{- i k_{y}}\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, 0],\n", - "[0, 0, e**(I*k_x) + e**(-I*k_x), 0],\n", - "[0, 0, 0, e**(I*k_y) + e**(-I*k_y)]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & - e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\\\e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}} & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, -e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x), 0, 0],\n", - "[e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x), 0, 0, 0],\n", - "[ 0, 0, 0, 0],\n", - "[ 0, 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & i e^{i k_{y}} + i e^{- i k_{y}} & 0\\\\0 & 0 & 0 & - i e^{i k_{x}} - i e^{- i k_{x}}\\\\- i e^{i k_{y}} - i e^{- i k_{y}} & 0 & 0 & 0\\\\0 & i e^{i k_{x}} + i e^{- i k_{x}} & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, I*e**(I*k_y) + I*e**(-I*k_y), 0],\n", - "[ 0, 0, 0, -I*e**(I*k_x) - I*e**(-I*k_x)],\n", - "[-I*e**(I*k_y) - I*e**(-I*k_y), 0, 0, 0],\n", - "[ 0, I*e**(I*k_x) + I*e**(-I*k_x), 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & i e^{i k_{y}} - i e^{- i k_{y}}\\\\0 & 0 & e^{i k_{x}} - e^{- i k_{x}} & 0\\\\0 & - e^{i k_{x}} + e^{- i k_{x}} & 0 & 0\\\\i e^{i k_{y}} - i e^{- i k_{y}} & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, 0, I*e**(I*k_y) - I*e**(-I*k_y)],\n", - "[ 0, 0, e**(I*k_x) - e**(-I*k_x), 0],\n", - "[ 0, -e**(I*k_x) + e**(-I*k_x), 0, 0],\n", - "[I*e**(I*k_y) - I*e**(-I*k_y), 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & - e^{i k_{x}} + e^{- i k_{x}}\\\\0 & 0 & i e^{i k_{y}} - i e^{- i k_{y}} & 0\\\\0 & i e^{i k_{y}} - i e^{- i k_{y}} & 0 & 0\\\\e^{i k_{x}} - e^{- i k_{x}} & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, 0, -e**(I*k_x) + e**(-I*k_x)],\n", - "[ 0, 0, I*e**(I*k_y) - I*e**(-I*k_y), 0],\n", - "[ 0, I*e**(I*k_y) - I*e**(-I*k_y), 0, 0],\n", - "[e**(I*k_x) - e**(-I*k_x), 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & - i e^{i k_{x}} - i e^{- i k_{x}} & 0\\\\0 & 0 & 0 & i e^{i k_{y}} + i e^{- i k_{y}}\\\\i e^{i k_{x}} + i e^{- i k_{x}} & 0 & 0 & 0\\\\0 & - i e^{i k_{y}} - i e^{- i k_{y}} & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, -I*e**(I*k_x) - I*e**(-I*k_x), 0],\n", - "[ 0, 0, 0, I*e**(I*k_y) + I*e**(-I*k_y)],\n", - "[I*e**(I*k_x) + I*e**(-I*k_x), 0, 0, 0],\n", - "[ 0, -I*e**(I*k_y) - I*e**(-I*k_y), 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & - e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}}\\\\0 & 0 & e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}} & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, -e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x)],\n", - "[0, 0, e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x), 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "familyC4T = qsymm.bloch_family(hopping_vectors, {C4T}, norbs=norbs)\n", - "print(len(familyC4T), [len(fam) for fam in familyC4T])\n", - "qsymm.display_family(familyC4T)" - ] - }, - { - "cell_type": "code", - "execution_count": 32, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "6 [4, 4, 4, 4, 4, 4]\n" - ] - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}- e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} - e^{- i k_{x}} & 0 & 0 & 0\\\\0 & e^{i k_{x}} - e^{i k_{y}} - e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[-e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) - e**(-I*k_x), 0, 0, 0],\n", - "[ 0, e**(I*k_x) - e**(I*k_y) - e**(-I*k_y) + e**(-I*k_x), 0, 0],\n", - "[ 0, 0, 0, 0],\n", - "[ 0, 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & - e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} - e^{- i k_{x}} & 0\\\\0 & 0 & 0 & e^{i k_{x}} - e^{i k_{y}} - e^{- i k_{y}} + e^{- i k_{x}}\\\\- e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} - e^{- i k_{x}} & 0 & 0 & 0\\\\0 & e^{i k_{x}} - e^{i k_{y}} - e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, -e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) - e**(-I*k_x), 0],\n", - "[ 0, 0, 0, e**(I*k_x) - e**(I*k_y) - e**(-I*k_y) + e**(-I*k_x)],\n", - "[-e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) - e**(-I*k_x), 0, 0, 0],\n", - "[ 0, e**(I*k_x) - e**(I*k_y) - e**(-I*k_y) + e**(-I*k_x), 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & - i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} + i e^{- i k_{x}}\\\\0 & 0 & i e^{i k_{x}} + e^{i k_{y}} - e^{- i k_{y}} - i e^{- i k_{x}} & 0\\\\0 & i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} - i e^{- i k_{x}} & 0 & 0\\\\- i e^{i k_{x}} - e^{i k_{y}} + e^{- i k_{y}} + i e^{- i k_{x}} & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, 0, -I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) + I*e**(-I*k_x)],\n", - "[ 0, 0, I*e**(I*k_x) + e**(I*k_y) - e**(-I*k_y) - I*e**(-I*k_x), 0],\n", - "[ 0, I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) - I*e**(-I*k_x), 0, 0],\n", - "[-I*e**(I*k_x) - e**(I*k_y) + e**(-I*k_y) + I*e**(-I*k_x), 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & - e^{i k_{x}} + e^{i k_{y}} + e^{- i k_{y}} - e^{- i k_{x}} & 0\\\\0 & 0 & 0 & e^{i k_{x}} - e^{i k_{y}} - e^{- i k_{y}} + e^{- i k_{x}}\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[0, 0, 0, 0],\n", - "[0, 0, 0, 0],\n", - "[0, 0, -e**(I*k_x) + e**(I*k_y) + e**(-I*k_y) - e**(-I*k_x), 0],\n", - "[0, 0, 0, e**(I*k_x) - e**(I*k_y) - e**(-I*k_y) + e**(-I*k_x)]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & - i e^{i k_{x}} + i e^{i k_{y}} + i e^{- i k_{y}} - i e^{- i k_{x}} & 0\\\\0 & 0 & 0 & - i e^{i k_{x}} + i e^{i k_{y}} + i e^{- i k_{y}} - i e^{- i k_{x}}\\\\i e^{i k_{x}} - i e^{i k_{y}} - i e^{- i k_{y}} + i e^{- i k_{x}} & 0 & 0 & 0\\\\0 & i e^{i k_{x}} - i e^{i k_{y}} - i e^{- i k_{y}} + i e^{- i k_{x}} & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, -I*e**(I*k_x) + I*e**(I*k_y) + I*e**(-I*k_y) - I*e**(-I*k_x), 0],\n", - "[ 0, 0, 0, -I*e**(I*k_x) + I*e**(I*k_y) + I*e**(-I*k_y) - I*e**(-I*k_x)],\n", - "[I*e**(I*k_x) - I*e**(I*k_y) - I*e**(-I*k_y) + I*e**(-I*k_x), 0, 0, 0],\n", - "[ 0, I*e**(I*k_x) - I*e**(I*k_y) - I*e**(-I*k_y) + I*e**(-I*k_x), 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/latex": [ - "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} - e^{- i k_{x}}\\\\0 & 0 & e^{i k_{x}} - i e^{i k_{y}} + i e^{- i k_{y}} - e^{- i k_{x}} & 0\\\\0 & - e^{i k_{x}} - i e^{i k_{y}} + i e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0\\\\- e^{i k_{x}} + i e^{i k_{y}} - i e^{- i k_{y}} + e^{- i k_{x}} & 0 & 0 & 0\\end{matrix}\\right]$" - ], - "text/plain": [ - "Matrix([\n", - "[ 0, 0, 0, e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) - e**(-I*k_x)],\n", - "[ 0, 0, e**(I*k_x) - I*e**(I*k_y) + I*e**(-I*k_y) - e**(-I*k_x), 0],\n", - "[ 0, -e**(I*k_x) - I*e**(I*k_y) + I*e**(-I*k_y) + e**(-I*k_x), 0, 0],\n", - "[-e**(I*k_x) + I*e**(I*k_y) - I*e**(-I*k_y) + e**(-I*k_x), 0, 0, 0]])" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "family_new = qsymm.hamiltonian_generator.subtract_family(familyC4T, familyTR, prettify=True)\n", - "print(len(family_new), [len(fam) for fam in family_new])\n", - "qsymm.display_family(family_new)" - ] - }, - { - "cell_type": "code", - "execution_count": 37, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [ - { - "data": { - "text/plain": [ - "(1-0j)" - ] - }, - "execution_count": 37, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "from scipy import linalg as la\n", - "la.det(np.kron(np.eye(4), 2 * S[1]))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Find periods from BlochModel" - ] - }, - { - "cell_type": "code", - "execution_count": 108, - "metadata": {}, - "outputs": [], - "source": [ - "ham = qsymm.hamiltonian_generator.hamiltonian_from_family(family, tosympy=False)" - ] - }, - { - "cell_type": "code", - "execution_count": 109, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "{(array([0.5 , 0.28867513]), c0):\n", - "[[0.+0.j 0.+0.j]\n", - " [1.+0.j 0.+0.j]],\n", - "\n", - "(array([ 0.5 , -0.28867513]), c0):\n", - "[[0.+0.j 1.+0.j]\n", - " [0.+0.j 0.+0.j]],\n", - "\n", - "(array([-0.5 , -0.28867513]), c0):\n", - "[[0.+0.j 1.+0.j]\n", - " [0.+0.j 0.+0.j]],\n", - "\n", - "(array([-0. , 0.57735027]), c0):\n", - "[[0.+0.j 1.+0.j]\n", - " [0.+0.j 0.+0.j]],\n", - "\n", - "(array([-0.5 , 0.28867513]), c0):\n", - "[[0.+0.j 0.+0.j]\n", - " [1.+0.j 0.+0.j]],\n", - "\n", - "(array([-0. , -0.57735027]), c0):\n", - "[[0.+0.j 0.+0.j]\n", - " [1.+0.j 0.+0.j]],\n", - "\n", - "}" - ] - }, - "execution_count": 109, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "ham = qsymm.BlochModel(ham)\n", - "ham" - ] - }, - { - "cell_type": "code", - "execution_count": 148, - "metadata": {}, - "outputs": [], - "source": [ - "from collections import defaultdict\n", - "from qsymm.model import BlochModel, BlochCoeff\n", - "from qsymm.kwant_linalg_lll import lll\n", - "from itertools import product\n", - "from functools import reduce\n", - "from math import gcd\n", - "# from https://github.com/lan496/hsnf\n", - "from hsnf import smith_normal_form, row_style_hermite_normal_form\n", - "\n", - "one = sympy.sympify(1)\n", - "\n", - "def model_periods(model):\n", - " \"\"\"Find the periods and positions of orbitals\n", - " for the lattice model where it is is periodic.\n", - " It needs to be a bloch Hamiltonian written in the\n", - " real space convention (i.e. hopping phases correspond\n", - " to the real space sepatarion of orbitals)\"\"\"\n", - " # convert to BlochModel\n", - " model = BlochModel(model)\n", - " n = model.shape[0]\n", - " d = len(model.momenta)\n", - "\n", - " # dictionary of relative positions for each pair of\n", - " # orbitals, allow multiple relative positions\n", - " rel_pos_nn = defaultdict(lambda: set())\n", - " # populate it with all the hoppings in the model\n", - " for (hop, _), val in model.items():\n", - " for a, b in zip(*np.nonzero(val)):\n", - " rel_pos_nn[a, b].add(BlochCoeff(hop, one))\n", - "\n", - " # Keep adding further neighbor relative positions up to n steps,\n", - " # this guarantees we find all linearly independent (over integers)\n", - " # lattice vectors.\n", - " # This may take long if there are many orbitals, perhaps there is\n", - " # a better stopping condition.\n", - " rel_pos = defaultdict(lambda: set(), {ind: set.copy() for ind, set in rel_pos_nn.items()})\n", - " for _ in range(n-1):\n", - " for ((a1, b1), hops1), ((a2, b2), hops2) in product(rel_pos.items(), rel_pos_nn.items()):\n", - " if b1 == a2:\n", - " for (hop1, hop2) in product(hops1, hops2):\n", - " rel_pos[a1, b2].add(hop1 * hop2)\n", - "\n", - " # Lattice vectors are those connecting the same site type\n", - " vecs = np.array([vec for ((a, b), hops) in rel_pos.items() if a == b\n", - " for vec, _ in hops])\n", - " # Find primitive vectors\n", - " prim_vecs = primitive_lattice_vecs(vecs)\n", - " return prim_vecs\n", - "\n", - "\n", - "def make_int(R):\n", - " # If close to an integer array convert to integer array, else\n", - " # return None\n", - " R_int = np.array(np.round(R), int)\n", - " if qsymm.linalg.allclose(R, R_int):\n", - " return R_int\n", - " else:\n", - " return None\n", - "\n", - "\n", - "def primitive_lattice_vecs(vecs, method='hnf'):\n", - " \"\"\"Find a set of primitive lattice vectors spanning the\n", - " same lattice as vecs.\"\"\"\n", - " vecs = np.asarray(vecs)\n", - " d = vecs.shape[1]\n", - " # pick out a set of small linearly independent vectors\n", - " ### TODO this likely works even if vecs don't span the full space\n", - " norms = np.linalg.norm(vecs, axis=1)\n", - " vecs = vecs[np.argsort(norms)]\n", - " bas = []\n", - " for v in vecs:\n", - " new_bas = bas + [v]\n", - " rank = np.linalg.matrix_rank(new_bas)\n", - " if rank == len(new_bas):\n", - " bas = new_bas\n", - " if rank == d:\n", - " break\n", - " else:\n", - " raise ValueError('The vectors do not span the full space.')\n", - " bas = np.vstack(bas)\n", - " # rewrite all vectors in this new basis\n", - " vecs = np.linalg.solve(bas.T, vecs.T).T\n", - " # find the common denominator to make all integer\n", - " denom = lcd(vecs)\n", - " vecs = make_int(vecs * denom)\n", - " bas = bas / denom\n", - " # find primitive lattice vectors\n", - " if method == 'snf':\n", - " D, L, R = smith_normal_form(vecs)\n", - " prim_vecs = D @ np.linalg.inv(R)\n", - " elif method == 'hnf':\n", - " prim_vecs, U = row_style_hermite_normal_form(vecs)\n", - " else:\n", - " raise ValueError(\"Method must be either 'hnf' or 'snf'.\")\n", - " prim_vecs = prim_vecs[:d]\n", - " # convert back to original basis\n", - " prim_vecs = prim_vecs @ bas\n", - " return lll(prim_vecs)[0]\n", - "\n", - "\n", - "def lcd(a, tol=1e-9):\n", - " \"\"\"Calculate approximate least common denominator for\n", - " the entries of floating point array `a`.\"\"\"\n", - " a = np.asarray(a).flatten()\n", - " a = a % 1\n", - " ### This would be nice to vectorize\n", - " denoms = np.array([farey(x)[1] for x in a])\n", - " return reduce(lambda a, b: a * b // gcd(a, b), denoms, 1)\n", - " \n", - "\n", - "def farey(x, tol=1e-9):\n", - " \"\"\"Farey sequence rational approximation of `0 <= x <= 1`\n", - " with tolerance `tol`. Adapted from\n", - " https://www.johndcook.com/blog/2010/10/20/best-rational-approximation/\n", - " \"\"\"\n", - " N = int(tol**(-1/2) / np.sqrt(5))\n", - " a, b = 0, 1\n", - " c, d = 1, 1\n", - " while (b <= N and d <= N):\n", - " # Add these checks to avoid additional loops\n", - " if abs(x - a/b) <= tol:\n", - " return a, b\n", - " elif abs(x - c/d) <= tol:\n", - " return c, d\n", - " mediant = (a+c)/(b+d)\n", - " if abs(x - mediant) <= tol:\n", - " return a+c, b+d\n", - " elif x > mediant:\n", - " a, b = a+c, b+d\n", - " else:\n", - " c, d = a+c, b+d\n", - "\n", - " if (b > N):\n", - " return c, d\n", - " else:\n", - " return a, b" - ] - }, - { - "cell_type": "code", - "execution_count": 149, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "55.8 ms ± 4.06 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" - ] - }, - { - "data": { - "text/plain": [ - "array([[ 1.00000000e+00, -2.77555756e-16],\n", - " [ 5.00000000e-01, 8.66025404e-01]])" - ] - }, - "execution_count": 149, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# %%time\n", - "%timeit lat_vecs = model_periods(ham)\n", - "lat_vecs" - ] - }, - { - "cell_type": "code", - "execution_count": 114, - "metadata": {}, - "outputs": [], - "source": [ - "%load_ext line_profiler" - ] - }, - { - "cell_type": "code", - "execution_count": 150, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "Timer unit: 1e-06 s\n", - "\n", - "Total time: 0.000181 s\n", - "File: \n", - "Function: farey at line 111\n", - "\n", - "Line # Hits Time Per Hit % Time Line Contents\n", - "==============================================================\n", - " 111 def farey(x, tol=1e-9):\n", - " 112 \"\"\"Farey sequence rational approximation of `0 <= x <= 1`\n", - " 113 with tolerance `tol`. Adapted from\n", - " 114 https://www.johndcook.com/blog/2010/10/20/best-rational-approximation/\n", - " 115 \"\"\"\n", - " 116 28 81.0 2.9 44.8 N = int(tol**(-1/2) / np.sqrt(5))\n", - " 117 28 13.0 0.5 7.2 a, b = 0, 1\n", - " 118 28 11.0 0.4 6.1 c, d = 1, 1\n", - " 119 28 13.0 0.5 7.2 while (b <= N and d <= N):\n", - " 120 # Add these checks to avoid additional loops\n", - " 121 28 40.0 1.4 22.1 if abs(x - a/b) < tol:\n", - " 122 24 16.0 0.7 8.8 return a, b\n", - " 123 4 4.0 1.0 2.2 elif abs(x - c/d) < tol:\n", - " 124 4 3.0 0.8 1.7 return c, d\n", - " 125 mediant = (a+c)/(b+d)\n", - " 126 if x == mediant:\n", - " 127 return a+c, b+d\n", - " 128 elif x > mediant:\n", - " 129 a, b = a+c, b+d\n", - " 130 else:\n", - " 131 c, d = a+c, b+d\n", - " 132 \n", - " 133 if (b > N):\n", - " 134 return c, d\n", - " 135 else:\n", - " 136 return a, b" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "%lprun -f farey model_periods(ham)" - ] - }, - { - "cell_type": "code", - "execution_count": 88, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "CPU times: user 142 ms, sys: 24 µs, total: 142 ms\n", - "Wall time: 138 ms\n" - ] - }, - { - "data": { - "text/plain": [ - "array([[1., 0.],\n", - " [0., 1.]])" - ] - }, - "execution_count": 88, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "%%time\n", - "primitive_lattice_vecs([[11, 1], [1, 0], [0, 0]], method='hnf')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python [conda env:.conda-kwant_conda]", - "language": "python", - "name": "conda-env-.conda-kwant_conda-py" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.8" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/qsymm/linalg.py b/qsymm/linalg.py index 05a4a74..c75d80a 100644 --- a/qsymm/linalg.py +++ b/qsymm/linalg.py @@ -4,6 +4,8 @@ import scipy.linalg as la import scipy.sparse.linalg as sla import scipy from numbers import Number +from math import gcd +from functools import reduce from scipy.optimize import minimize from scipy.spatial.distance import cdist from scipy.sparse.csgraph import connected_components @@ -610,3 +612,41 @@ def _inv_int(A): if _A != A or abs(la.det(A)) != 1: raise ValueError('Input needs to be an invertible integer matrix') return ta.array(np.round(la.inv(_A)), int) + + +def lcd(a, tol=1e-9): + """Calculate approximate least common denominator for + the entries of floating point array `a`.""" + a = np.asarray(a).flatten() + a = a % 1 + denoms = np.array([farey(x)[1] for x in a]) + return reduce(lambda a, b: a * b // gcd(a, b), denoms, 1) + + +def farey(x, tol=1e-5): + """Farey sequence rational approximation of `0 <= x <= 1` + with tolerance `tol`. Adapted from + https://www.johndcook.com/blog/2010/10/20/best-rational-approximation/ + """ + # Would be nice to vectorize this function + N = int(1/tol) + a, b = 0, 1 + c, d = 1, 1 + while (b <= N and d <= N): + # Add these checks to avoid additional loops + if abs(x - a/b) <= tol: + return a, b + elif abs(x - c/d) <= tol: + return c, d + mediant = (a+c)/(b+d) + if abs(x - mediant) <= tol: + return a+c, b+d + elif x > mediant: + a, b = a+c, b+d + else: + c, d = a+c, b+d + + if (b > N): + return c, d + else: + return a, b diff --git a/qsymm/symmetry_finder.py b/qsymm/symmetry_finder.py index 6056201..e4fac0b 100644 --- a/qsymm/symmetry_finder.py +++ b/qsymm/symmetry_finder.py @@ -5,11 +5,13 @@ import itertools as it import numpy as np import scipy.linalg as la import scipy.sparse +# from https://github.com/lan496/hsnf +from hsnf import smith_normal_form, row_style_hermite_normal_form from .linalg import matrix_basis, nullspace, split_list, simult_diag, commutator, \ prop_to_id, sparse_basis, mtm, family_to_vectors, solve_mat_eqn, \ - allclose -from .model import Model, BlochModel, BlochCoeff + allclose, lcd +from .model import Model, BlochModel, BlochCoeff, One from .groups import PointGroupElement, ContinuousGroupGenerator, generate_group, \ set_multiply, L_matrices, spin_rotation, time_reversal, \ particle_hole, chiral, inversion, rotation, mirror, PrettyList @@ -31,7 +33,7 @@ def symmetries(model, candidates=None, continuous_rotations=False, Model or BlochModel which represents family of Hamiltonians. Every symbolic prefactor is treated as a free parameter, and model.momenta as independent momentum variables. - candidates : iterable of PointGroupElements or None + candidates : iterable of PointGroupElements or None or 'auto' Set of candidate PointGroupElements used for finding discrete symmetries. Must have .U attribute set to None. It is advised that candidates forms a group, as combinations of @@ -42,6 +44,9 @@ def symmetries(model, candidates=None, continuous_rotations=False, Models and BlochModels. If None (default), only discrete onsite symmetries (time reversal, particle-hole, chiral) are found. + If 'auto' the candidates are inferred from the periodicity of the + Hamiltonian. Only works for Bloch Hamiltonians (Models that are + convertible to BlochModel). continuous_rotations : bool (default False) Whether to search for continuous rotation symmetries. generators : bool (default False) @@ -90,6 +95,13 @@ def symmetries(model, candidates=None, continuous_rotations=False, # Make discrete onsite symmetries dim = len(model.momenta) candidates = generate_group({time_reversal(dim), particle_hole(dim), chiral(dim)}) + elif candidates == 'auto': + try: + model = BlochModel(model) + except: + raise ValueError("candidates='auto' only works if model is convertible to BlochModel.") + periods = model_periods(model) + candidates = bravais_point_group(periods, tr=True, ph=True) if candidates: disc_sym, _ = discrete_symmetries(model, set(candidates), Ps=Ps, @@ -991,7 +1003,6 @@ def bravais_group_3d(neighbors, num_eq, sets_eq, verbose=False): def equals(vectors): # group equivalent vectors based on length and angles - one = kwant_continuum.sympify('1') vector_sets = defaultdict(list) # Take abs because every vector has opposite pair overlaps = np.abs(vectors @ vectors.T) @@ -1000,7 +1011,7 @@ def equals(vectors): length = np.array([overlaps[i, i]]) # Symmetry equivalent vectors must have the same signature signature = np.concatenate([length, sorted(overlaps[i]), sorted(angles[i])]) - key = BlochCoeff(signature, one) + key = BlochCoeff(signature, One()) vector_sets[key].append(vector) vector_sets = sorted( @@ -1063,12 +1074,140 @@ def twofold_axis(vectors, neighbors): def check_bravais_symmetry(neighbors, group): - one = kwant_continuum.sympify('1') neighbors = np.vstack([neighbors, -neighbors]) - neighbors = {BlochCoeff(vec, one) for vec in neighbors} + neighbors = {BlochCoeff(vec, One()) for vec in neighbors} for g in group: r_neighbors = {BlochCoeff(g.R @ hop, coeff) for (hop, coeff) in neighbors} if not neighbors == r_neighbors: return False else: return True + + +def model_periods(model, return_coords=False, norbs=None): + """Find the periods and positions of orbitals + for the lattice model where it is is periodic. + It needs to be a bloch Hamiltonian written in the + real space convention (i.e. hopping phases correspond + to the real space sepatarion of orbitals). + + Parameters: + ----------- + model : qsymm.Model + Lattice Hamiltonian (needs to be convertible to BlochModel). + + return_coords: boolean, default False + Whether to return the orbital coordinates as well, relative + to the first orbital. + + norbs : iterable of integers, default None + Sequence of number of consecutive orbitals belonging to the + same site. If `None`, all orbitals are treated as separate + sites. + + Returns: + -------- + prim_vecs : ndarray + Primitive lattice vectors of the translation symmetries as + row vectors. + + atom_coords : ndarray + Site coordinates as row vectors. + """ + # convert to BlochModel + model = BlochModel(model) + d = len(model.momenta) + + if norbs is None: + norbs = [1] * model.shape[0] + + N, ranges = 0, [] + for n in norbs: + ranges.append(slice(N, N + n)) + N += n + n = len(ranges) + + # dictionary of relative positions for each pair of + # orbitals, allow multiple relative positions + rel_pos_nn = defaultdict(lambda: set()) + # populate it with all the hoppings in the model + for (hop, _), val in model.items(): + for a, b in it.product(range(n), repeat=2): + if np.count_nonzero(val[ranges[a], ranges[b]]) > 0: + rel_pos_nn[a, b].add(BlochCoeff(hop, One())) + + # Keep adding further neighbor relative positions up to n steps, + # this guarantees we find all linearly independent (over integers) + # lattice vectors. + # This may take long if there are many orbitals, perhaps there is + # a better stopping condition. + rel_pos = defaultdict(lambda: set(), {ind: set.copy() for ind, set in rel_pos_nn.items()}) + for i in range(n-1): + for ((a1, b1), hops1), ((a2, b2), hops2) in it.product(rel_pos.items(), rel_pos_nn.items()): + if b1 == a2: + for (hop1, hop2) in it.product(hops1, hops2): + rel_pos[a1, b2].add(hop1 * hop2) + + # Lattice vectors are those connecting the same site type + vecs = np.array([vec for ((a, b), hops) in rel_pos.items() if a == b + for vec, _ in hops]) + # Find primitive vectors + prim_vecs = primitive_lattice_vecs(vecs) + + if not return_coords: + return prim_vecs + + # Find coordinates relative to first orbital + atom_coords = [] + for a in range(n): + coords = np.array([vec for vec, _ in rel_pos[0, a]]) + if len(coords) == 0: + raise ValueError("The TB model does not form a connected graph. Try specifying `norbs`.") + pos = coords[0] + pos -= cvp(pos, prim_vecs)[0] @ prim_vecs + atom_coords.append(pos) + atom_coords = np.array(atom_coords) + + return prim_vecs, atom_coords + + +def primitive_lattice_vecs(vecs, method='hnf'): + """Find a set of primitive lattice vectors spanning the + same lattice as vecs.""" + vecs = np.asarray(vecs) + d = vecs.shape[1] + # pick out a set of small linearly independent vectors + ### TODO this likely works even if vecs don't span the full space + norms = np.linalg.norm(vecs, axis=1) + vecs = vecs[np.argsort(norms)] + bas = [] + for v in vecs: + new_bas = bas + [v] + rank = np.linalg.matrix_rank(new_bas) + if rank == len(new_bas): + bas = new_bas + if rank == d: + break + else: + raise ValueError('The vectors do not span the full space.') + bas = np.vstack(bas) + # rewrite all vectors in this new basis + vecs = np.linalg.solve(bas.T, vecs.T).T + # find the common denominator to make all integer + denom = lcd(vecs) + bas = bas / denom + vecs = vecs * denom + vecs_int = np.array(np.round(vecs), int) + assert allclose(vecs, vecs_int) + # find primitive lattice vectors + if method == 'snf': + D, L, R = smith_normal_form(vecs_int) + prim_vecs = D @ np.linalg.inv(R) + elif method == 'hnf': + prim_vecs, U = row_style_hermite_normal_form(vecs_int) + else: + raise ValueError("Method must be either 'hnf' or 'snf'.") + prim_vecs = prim_vecs[:d] + # convert back to original basis + prim_vecs = prim_vecs @ bas + return lll(prim_vecs)[0] diff --git a/qsymm/tests/test_symmetry_finder.py b/qsymm/tests/test_symmetry_finder.py index 36aa2a0..46f5489 100644 --- a/qsymm/tests/test_symmetry_finder.py +++ b/qsymm/tests/test_symmetry_finder.py @@ -7,6 +7,7 @@ from ..symmetry_finder import * from ..symmetry_finder import _reduced_model, _reduce_hamiltonian, bravais_point_group from ..linalg import * from .. import kwant_continuum +from ..groups import * sigma = np.array([[[1, 0], [0, 1]], [[0, 1], [ 1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]) @@ -533,7 +534,7 @@ def test_bloch(): False, False, None) TR = time_reversal(realspace_dim=2) PH = particle_hole(realspace_dim=2) - gens_hex_2D ={Mx, C6, TR, PH} + gens_hex_2D = {Mx, C6, TR, PH} hex_group_2D = generate_group(gens_hex_2D) assert len(hex_group_2D) == 48 @@ -544,6 +545,10 @@ def test_bloch(): assert [P.shape for P in Ps] == [(1, 1, 1)] assert len(sg) == 24 assert sg == generate_group({Mx, C6, TR}) + # Test automatic candidates + pg, _ = symmetries(H6, candidates='auto') + assert len(pg) == 24 + assert set(pg) == hexagonal(sympy_R=False, tr=True, ph=False) # extend model to add SOC ham62 = 'eye(2) * (' + ham6 + ') +' -- GitLab From 0b4345198b8e49d3f9e551d6ceb68eb9c05b3a69 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Mon, 11 Jan 2021 01:49:52 +0100 Subject: [PATCH 4/4] optimize period finding --- qsymm/symmetry_finder.py | 50 +++++++++++++++++++++-------- qsymm/tests/test_symmetry_finder.py | 4 +++ 2 files changed, 40 insertions(+), 14 deletions(-) diff --git a/qsymm/symmetry_finder.py b/qsymm/symmetry_finder.py index e4fac0b..bf61a9a 100644 --- a/qsymm/symmetry_finder.py +++ b/qsymm/symmetry_finder.py @@ -1139,20 +1139,42 @@ def model_periods(model, return_coords=False, norbs=None): # Keep adding further neighbor relative positions up to n steps, # this guarantees we find all linearly independent (over integers) # lattice vectors. - # This may take long if there are many orbitals, perhaps there is - # a better stopping condition. - rel_pos = defaultdict(lambda: set(), {ind: set.copy() for ind, set in rel_pos_nn.items()}) - for i in range(n-1): - for ((a1, b1), hops1), ((a2, b2), hops2) in it.product(rel_pos.items(), rel_pos_nn.items()): - if b1 == a2: - for (hop1, hop2) in it.product(hops1, hops2): - rel_pos[a1, b2].add(hop1 * hop2) - - # Lattice vectors are those connecting the same site type - vecs = np.array([vec for ((a, b), hops) in rel_pos.items() if a == b - for vec, _ in hops]) - # Find primitive vectors - prim_vecs = primitive_lattice_vecs(vecs) + rel_pos = defaultdict(lambda: set(), + {(a, a): {BlochCoeff(np.zeros((d,)), One())} for a in range(n)}) + rel_pos_next = defaultdict(lambda: set()) + for i in range(n): + for a, b in it.product(range(n), repeat=2): + rel_pos_next[a, b] = {hop1 * hop2 + for c in range(n) + for (hop1, hop2) in it.product(rel_pos[a, c], rel_pos_nn[c, b]) + } + for a, b in it.product(range(n), repeat=2): + rel_pos[a, b] |= rel_pos_next[a, b] + + # Lattice vectors are those connecting the same site type + vecs = np.array([vec for ((a, b), hops) in rel_pos.items() if a == b + for vec, _ in hops]) + + # Find primitive vectors + try: + prim_vecs = primitive_lattice_vecs(vecs) + except ValueError: + # there weren't enough vectors, add more + continue + + # Take realtive positions modulo lattice vectors + for a, b in it.product(range(n), repeat=2): + rel_pos[a, b] = {BlochCoeff((np.linalg.solve(prim_vecs.T, vec) % 1) @ prim_vecs, + One()) + for vec, _ in rel_pos[a, b]} + rel_pos_next[a, b] = {BlochCoeff((np.linalg.solve(prim_vecs.T, vec) % 1) @ prim_vecs, + One()) + for vec, _ in rel_pos_next[a, b]} + + if all(len(rel_pos_next[a, b] - rel_pos[a, b]) == 0 + for a, b in it.product(range(n), repeat=2)): + # There were no new vectors modulo lattice translations, stop + break if not return_coords: return prim_vecs diff --git a/qsymm/tests/test_symmetry_finder.py b/qsymm/tests/test_symmetry_finder.py index 46f5489..b45e95b 100644 --- a/qsymm/tests/test_symmetry_finder.py +++ b/qsymm/tests/test_symmetry_finder.py @@ -559,6 +559,10 @@ def test_bloch(): assert [P.shape for P in Ps] == [(1, 2, 2)] assert len(sg) == 24 assert sg == generate_group({Mx, C6, TR}) + # Test automatic candidates + pg, _ = symmetries(H62, candidates='auto') + assert len(pg) == 24 + assert set(pg) == hexagonal(sympy_R=False, tr=True, ph=False) # Add degeneracy ham63 = 'kron(eye(2), ' + ham62 + ')' -- GitLab