{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sympy\n",
    "\n",
    "import qsymm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & vt \\left(k_{x} - i k_{y}\\right) & vs \\left(k_{x} - i k_{y}\\right) & 0\\\\vt \\left(k_{x} + i k_{y}\\right) & 0 & 0 & vs \\left(k_{x} - i k_{y}\\right)\\\\vs \\left(k_{x} + i k_{y}\\right) & 0 & 0 & vt \\left(k_{x} - i k_{y}\\right)\\\\0 & vs \\left(k_{x} + i k_{y}\\right) & vt \\left(k_{x} + i k_{y}\\right) & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[               0, vt*(k_x - I*k_y), vs*(k_x - I*k_y),                0],\n",
       "[vt*(k_x + I*k_y),                0,                0, vs*(k_x - I*k_y)],\n",
       "[vs*(k_x + I*k_y),                0,                0, vt*(k_x - I*k_y)],\n",
       "[               0, vs*(k_x + I*k_y), vt*(k_x + I*k_y),                0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ham_kek_y = \"\"\"\n",
    "    vs * kron((k_x * sigma_x + k_y * sigma_y), eye(2))\n",
    "    + vt * kron(eye(2), (k_x * sigma_x + k_y * sigma_y))\n",
    "\"\"\"\n",
    "display(qsymm.sympify(ham_kek_y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_kek_y= qsymm.Model(ham_kek_y, momenta=['k_x', 'k_y'])\n",
    "candidates = qsymm.groups.hexagonal()\n",
    "sgy, cg = qsymm.symmetries(H_kek_y, candidates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_subgroups = qsymm.groups.generate_subgroups(sgy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset({1, R(π) C})"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Terms proportional to the identity matrix don't open a gap.\n",
    "identity_terms = [qsymm.Model({qsymm.sympify(1): np.eye(4)})]\n",
    "\n",
    "dd_sg = []\n",
    "for subgroup, gens in all_subgroups.items():\n",
    "    family = qsymm.continuum_hamiltonian(gens, 2, 0)\n",
    "    family = qsymm.hamiltonian_generator.subtract_family(family, identity_terms)\n",
    "    if family == []:\n",
    "        dd_sg.append(subgroup)\n",
    "\n",
    "dd_sg.sort(key=len)\n",
    "\n",
    "dd_sg[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset({1, R(2π/3), R(-2π/3), C, R(2π/3) C, R(-2π/3) C})"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rdd_sg = []\n",
    "for subgroup, gens in all_subgroups.items():\n",
    "    invs = {\n",
    "        qsymm.inversion(2) * g\n",
    "        for g in [\n",
    "            qsymm.identity(2),\n",
    "            qsymm.time_reversal(2),\n",
    "            qsymm.particle_hole(2),\n",
    "            qsymm.chiral(2),\n",
    "        ]\n",
    "    }\n",
    "    # Skip subgroups that have inversion\n",
    "    if subgroup & invs:\n",
    "        continue\n",
    "    family = qsymm.continuum_hamiltonian(gens, 2, 0)\n",
    "    family = qsymm.hamiltonian_generator.subtract_family(family, identity_terms)\n",
    "    if family == []:\n",
    "        rdd_sg.append(subgroup)\n",
    "\n",
    "rdd_sg.sort(key=len)\n",
    "\n",
    "rdd_sg[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "any(g.antisymmetry for sg in rdd_sg for g in sg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$e^{-i\\phi L}H(\\mathbf{k})e^{i\\phi L} = H(R_{\\phi}\\mathbf{k}) \\\\R_{{\\phi}} = R\\left(\\phi\\right)\\\\L = \\begin{bmatrix}1.+0.j&0.+0.j&0.+0.j&0.+0.j\\\\0.+0.j&0.+0.j&0.+0.j&0.+0.j\\\\0.+0.j&0.+0.j&0.+0.j&0.+0.j\\\\ 0.+0.j& 0.+0.j& 0.+0.j&-1.+0.j\\\\\\end{bmatrix} \\\\$"
      ],
      "text/plain": [
       "[\n",
       "exp(-i ϕ L)⋅H(k)⋅exp(i ϕ L) = H(R_ϕ⋅k)R_ϕ = R(ϕ)\n",
       "\n",
       "L = [[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
       "     [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
       "     [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
       "     [ 0.+0.j  0.+0.j  0.+0.j -1.+0.j]]]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_, cg = qsymm.symmetries(H_kek_y, candidates=[], continuous_rotations=True, prettify=True)\n",
    "cg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & vs \\left(k_{x} - i k_{y}\\right) & 0 & - \\Delta\\\\vs \\left(k_{x} + i k_{y}\\right) & 0 & \\Delta & 0\\\\0 & \\Delta & 0 & vs \\left(k_{x} + i k_{y}\\right)\\\\- \\Delta & 0 & vs \\left(k_{x} - i k_{y}\\right) & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[               0, vs*(k_x - I*k_y),                0,           -Delta],\n",
       "[vs*(k_x + I*k_y),                0,            Delta,                0],\n",
       "[               0,            Delta,                0, vs*(k_x + I*k_y)],\n",
       "[          -Delta,                0, vs*(k_x - I*k_y),                0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ham_kek_o= \"\"\"\n",
    "    vs * (\n",
    "        k_x * kron(eye(2), sigma_x)\n",
    "        + k_y * kron(sigma_z, sigma_y)\n",
    "    )\n",
    "    + Delta * kron(sigma_y, sigma_y)\n",
    "\"\"\"\n",
    "display(qsymm.sympify(ham_kek_o))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\left[ H(\\mathbf{{k}}), L \\right] = 0 \\\\L = \\begin{bmatrix} 0. +0.j&-0. +0.j&-0. +0.j& 0.5+0.j\\\\-0. -0.j& 0. +0.j& 0.5-0.j&-0. -0.j\\\\-0. -0.j& 0.5+0.j& 0. +0.j&-0. +0.j\\\\ 0.5-0.j&-0. +0.j&-0. -0.j& 0. +0.j\\\\\\end{bmatrix} \\\\\\\\e^{-i\\phi L}H(\\mathbf{k})e^{i\\phi L} = H(R_{\\phi}\\mathbf{k}) \\\\R_{{\\phi}} = R\\left(\\phi\\right)\\\\L = \\begin{bmatrix} 0.5+0.j&-0. -0.j&-0. -0.j&-0. +0.j\\\\-0. +0.j&-0.5-0.j& 0. -0.j&-0. +0.j\\\\-0. +0.j& 0. +0.j&-0.5+0.j&-0. +0.j\\\\-0. -0.j&-0. -0.j&-0. -0.j& 0.5+0.j\\\\\\end{bmatrix} \\\\$"
      ],
      "text/plain": [
       "[\n",
       "[H(k), L] = 0\n",
       "L = [[ 0. +0.j -0. +0.j -0. +0.j  0.5+0.j]\n",
       "     [-0. -0.j  0. +0.j  0.5-0.j -0. -0.j]\n",
       "     [-0. -0.j  0.5+0.j  0. +0.j -0. +0.j]\n",
       "     [ 0.5-0.j -0. +0.j -0. -0.j  0. +0.j]],\n",
       " \n",
       "exp(-i ϕ L)⋅H(k)⋅exp(i ϕ L) = H(R_ϕ⋅k)R_ϕ = R(ϕ)\n",
       "\n",
       "L = [[ 0.5+0.j -0. -0.j -0. -0.j -0. +0.j]\n",
       "     [-0. +0.j -0.5-0.j  0. -0.j -0. +0.j]\n",
       "     [-0. +0.j  0. +0.j -0.5+0.j -0. +0.j]\n",
       "     [-0. -0.j -0. -0.j -0. -0.j  0.5+0.j]]]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_kek_o= qsymm.Model(ham_kek_o, momenta=['k_x', 'k_y'])\n",
    "_, cg = qsymm.symmetries(H_kek_o, candidates=[], continuous_rotations=True, prettify=True)\n",
    "cg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "sites = ['A1', 'B1', 'A2', 'B2', 'A3', 'B3']\n",
    "norbs = [(site, 1) for site in sites]\n",
    "\n",
    "# Time reversal\n",
    "TR = qsymm.PointGroupElement(sympy.eye(2), True, False, np.eye(6))\n",
    "\n",
    "# Chiral symmetry\n",
    "C = qsymm.PointGroupElement(sympy.eye(2), False, True, np.kron(np.eye(3), ([[1, 0], [0, -1]])))\n",
    "\n",
    "# Atom A rotates into B, B into A.\n",
    "sphi = 2*sympy.pi/6\n",
    "RC6 = sympy.Matrix([[sympy.cos(sphi), -sympy.sin(sphi)],\n",
    "                  [sympy.sin(sphi), sympy.cos(sphi)]])\n",
    "permC6 = {'A1': 'B1', 'B1': 'A2', 'A2': 'B2', 'B2': 'A3', 'A3': 'B3', 'B3': 'A1'}\n",
    "C6 = qsymm.groups.symmetry_from_permutation(RC6, permC6, norbs)\n",
    "\n",
    "RMx = sympy.Matrix([[-1, 0], [0, 1]])\n",
    "permMx = {'A1': 'A3', 'B1': 'B2', 'A2': 'A2', 'B2': 'B1', 'A3': 'A1', 'B3': 'B3'}\n",
    "Mx = qsymm.groups.symmetry_from_permutation(RMx, permMx, norbs)\n",
    "\n",
    "symmetries = [C, TR, C6, Mx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Only need the unique hoppings, one weak and one strong\n",
    "hopping_vectors = [\n",
    "    ('A1', 'B1', np.array([0, -1])), # around ring\n",
    "    ('A2', 'B3', np.array([0, -1])), # next UC\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & e^{- i k_{y}} & 0 & 0 & 0 & e^{\\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}}\\\\e^{i k_{y}} & 0 & e^{- \\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}} & 0 & 0 & 0\\\\0 & e^{\\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} & 0 & e^{\\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}} & 0 & 0\\\\0 & 0 & e^{- \\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} & 0 & e^{i k_{y}} & 0\\\\0 & 0 & 0 & e^{- i k_{y}} & 0 & 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}} & 0 & 0 & 0 & e^{- \\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}} & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[                                 0,                       e**(-I*k_y),                                   0,                                  0,                                   0, e**(I*k_y/2)*e**(-sqrt(3)*I*k_x/2)],\n",
       "[                        e**(I*k_y),                                 0, e**(-I*k_y/2)*e**(-sqrt(3)*I*k_x/2),                                  0,                                   0,                                  0],\n",
       "[                                 0, e**(I*k_y/2)*e**(sqrt(3)*I*k_x/2),                                   0, e**(I*k_y/2)*e**(-sqrt(3)*I*k_x/2),                                   0,                                  0],\n",
       "[                                 0,                                 0,  e**(-I*k_y/2)*e**(sqrt(3)*I*k_x/2),                                  0,                          e**(I*k_y),                                  0],\n",
       "[                                 0,                                 0,                                   0,                        e**(-I*k_y),                                   0,  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),                                 0,                                   0,                                  0, e**(-I*k_y/2)*e**(-sqrt(3)*I*k_x/2),                                  0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & e^{\\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} & 0 & 0\\\\0 & 0 & 0 & 0 & e^{- \\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} & 0\\\\0 & 0 & 0 & 0 & 0 & e^{- i k_{y}}\\\\e^{- \\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}} & 0 & 0 & 0 & 0 & 0\\\\0 & e^{\\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}} & 0 & 0 & 0 & 0\\\\0 & 0 & e^{i k_{y}} & 0 & 0 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[                                  0,                                  0,          0, e**(I*k_y/2)*e**(sqrt(3)*I*k_x/2),                                  0,           0],\n",
       "[                                  0,                                  0,          0,                                 0, e**(-I*k_y/2)*e**(sqrt(3)*I*k_x/2),           0],\n",
       "[                                  0,                                  0,          0,                                 0,                                  0, e**(-I*k_y)],\n",
       "[e**(-I*k_y/2)*e**(-sqrt(3)*I*k_x/2),                                  0,          0,                                 0,                                  0,           0],\n",
       "[                                  0, e**(I*k_y/2)*e**(-sqrt(3)*I*k_x/2),          0,                                 0,                                  0,           0],\n",
       "[                                  0,                                  0, e**(I*k_y),                                 0,                                  0,           0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# construct Hamiltonian terms to symmetry-allowed terms\n",
    "family = qsymm.bloch_family(hopping_vectors, symmetries, norbs=norbs)\n",
    "qsymm.display_family(family)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-4., -1., -1.,  1.,  1.,  4.])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Gamma_terms = [term.subs(dict(k_x=0, k_y=0)) for term in family]\n",
    "evals, evecs = np.linalg.eigh(\n",
    "    np.sum([\n",
    "        list(term.values())[0] * coeff\n",
    "        for term, coeff in zip(Gamma_terms, [1, 2])\n",
    "    ],\n",
    "    axis=0)\n",
    ")\n",
    "\n",
    "evals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.33333333, -0.33333333,  0.33333333, -0.33333333])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# projector to low energy subspace\n",
    "proj = evecs[:, 1:5]\n",
    "UC3 = (C6 * C6).U\n",
    "\n",
    "# projected symmetry operator\n",
    "pC3 = proj.T.conj() @ UC3 @ proj\n",
    "\n",
    "assert np.allclose(pC3 @ pC3.T.conj(), np.eye(4))\n",
    "\n",
    "evals2, evecs2 = np.linalg.eig(pC3)\n",
    "\n",
    "# Check angles of the eigenvalues\n",
    "np.angle(evals2) / (2 * np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "sites = ['A1', 'B1', 'A2', 'B2', 'A3', 'B3']\n",
    "norbs = [(site, 1) for site in sites]\n",
    "\n",
    "# Time reversal\n",
    "TR = qsymm.time_reversal(2, U=np.eye(6))\n",
    "\n",
    "# Chiral symmetry\n",
    "C = qsymm.chiral(2, U=np.kron(np.eye(3), ([[1, 0], [0, -1]])))\n",
    "\n",
    "# Atom A rotates into B, B into A.\n",
    "sphi = 2*sympy.pi/3\n",
    "RC3 = sympy.Matrix([\n",
    "    [sympy.cos(sphi), -sympy.sin(sphi)],\n",
    "    [sympy.sin(sphi), sympy.cos(sphi)],\n",
    "])\n",
    "permC3 = {'A1': 'A1', 'B1': 'B3', 'A2': 'A2', 'B2': 'B1', 'A3': 'A3', 'B3': 'B2'}\n",
    "C3 = qsymm.groups.symmetry_from_permutation(RC3, permC3, norbs)\n",
    "\n",
    "RMx = sympy.Matrix([[-1, 0], [0, 1]])\n",
    "permMx = {'A1': 'A1', 'B1': 'B1', 'A2': 'A3', 'B2': 'B3', 'A3': 'A2', 'B3': 'B2'}\n",
    "Mx = qsymm.groups.symmetry_from_permutation(RMx, permMx, norbs)\n",
    "\n",
    "symmetries = [C, TR, C3, Mx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Only need the unique hoppings, one weak and one strong\n",
    "hopping_vectors = [\n",
    "    ('A1', 'B1', np.array([0, -1])), # Y\n",
    "    ('A2', 'B3', np.array([0, -1])), # other\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & e^{- i k_{y}} & 0 & e^{\\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} & 0 & e^{\\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}}\\\\e^{i k_{y}} & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0\\\\e^{- \\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}} & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0\\\\e^{- \\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} & 0 & 0 & 0 & 0 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[                                  0, e**(-I*k_y), 0, e**(I*k_y/2)*e**(sqrt(3)*I*k_x/2), 0, e**(I*k_y/2)*e**(-sqrt(3)*I*k_x/2)],\n",
       "[                         e**(I*k_y),           0, 0,                                 0, 0,                                  0],\n",
       "[                                  0,           0, 0,                                 0, 0,                                  0],\n",
       "[e**(-I*k_y/2)*e**(-sqrt(3)*I*k_x/2),           0, 0,                                 0, 0,                                  0],\n",
       "[                                  0,           0, 0,                                 0, 0,                                  0],\n",
       "[ e**(-I*k_y/2)*e**(sqrt(3)*I*k_x/2),           0, 0,                                 0, 0,                                  0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & e^{- \\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}} & 0 & e^{- \\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} & 0\\\\0 & e^{\\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} & 0 & e^{\\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}} & 0 & e^{- i k_{y}}\\\\0 & 0 & e^{- \\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}} & 0 & e^{i k_{y}} & 0\\\\0 & e^{\\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}} & 0 & e^{- i k_{y}} & 0 & e^{\\frac{i k_{y}}{2}} e^{\\frac{\\sqrt{3} i k_{x}}{2}}\\\\0 & 0 & e^{i k_{y}} & 0 & e^{- \\frac{i k_{y}}{2}} e^{- \\frac{\\sqrt{3} i k_{x}}{2}} & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0,                                  0,                                   0,                                  0,                                   0,                                 0],\n",
       "[0,                                  0, e**(-I*k_y/2)*e**(-sqrt(3)*I*k_x/2),                                  0,  e**(-I*k_y/2)*e**(sqrt(3)*I*k_x/2),                                 0],\n",
       "[0,  e**(I*k_y/2)*e**(sqrt(3)*I*k_x/2),                                   0, e**(I*k_y/2)*e**(-sqrt(3)*I*k_x/2),                                   0,                       e**(-I*k_y)],\n",
       "[0,                                  0,  e**(-I*k_y/2)*e**(sqrt(3)*I*k_x/2),                                  0,                          e**(I*k_y),                                 0],\n",
       "[0, e**(I*k_y/2)*e**(-sqrt(3)*I*k_x/2),                                   0,                        e**(-I*k_y),                                   0, e**(I*k_y/2)*e**(sqrt(3)*I*k_x/2)],\n",
       "[0,                                  0,                          e**(I*k_y),                                  0, e**(-I*k_y/2)*e**(-sqrt(3)*I*k_x/2),                                 0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# construct Hamiltonian terms to symmetry-allowed terms\n",
    "family = qsymm.bloch_family(hopping_vectors, symmetries, norbs=norbs)\n",
    "qsymm.display_family(family)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-5.19615242e+00, -8.69969719e-18,  0.00000000e+00,  0.00000000e+00,\n",
       "        0.00000000e+00,  5.19615242e+00])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Gamma_terms = [term.subs(dict(k_x=0, k_y=0)) for term in family]\n",
    "evals, evecs = np.linalg.eigh(\n",
    "    np.sum([\n",
    "        list(term.values())[0] * coeff\n",
    "        for term, coeff in zip(Gamma_terms, [1, 2])\n",
    "    ],\n",
    "    axis=0)\n",
    ")\n",
    "\n",
    "evals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.        ,  0.        ,  0.33333333, -0.33333333])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# projector to low energy subspace\n",
    "proj = evecs[:, 1:5]\n",
    "UC3 = C3.U\n",
    "\n",
    "# projected symmetry operator\n",
    "pC3 = proj.T.conj() @ UC3 @ proj\n",
    "\n",
    "assert np.allclose(pC3 @ pC3.T.conj(), np.eye(4))\n",
    "\n",
    "evals2, evecs2 = np.linalg.eig(pC3)\n",
    "\n",
    "# Check angles of the eigenvalues\n",
    "np.angle(evals2) / (2 * np.pi)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}