Tight-Binding Model_SrCaCu2O4.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ccd756f3-44ee-4f2f-8f4c-ba06c91a888c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: could not identify MPI environment!\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Starting serial run at: 2023-08-11 12:38:05.129417\n"
     ]
    }
   ],
   "source": [
    "# importing necessary dependencies\n",
    "from itertools import product as itp\n",
    "from pythtb import *\n",
    "from triqs.lattice.tight_binding import TBLattice\n",
    "import sympy as sp\n",
    "from sympy import *\n",
    "import warnings\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1cbeb0cd",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def TB_to_sympy(TBL, analytical = True, precision = 6):\n",
    "    r\"\"\"\n",
    "    returns the analytical form of the momentum space hamiltonian of the tight-binding model \n",
    "    from a tight-binding lattice object by utilizing Fourier series\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    TBL: triqs TBLattice object\n",
    "        triqs tight binding object\n",
    "    analytical: boolean, default = True\n",
    "        whether to return the Hamiltonian in analytical (true) or numerical (false) form.\n",
    "    precision: integer, default = 6\n",
    "        specifies the number of digits in the floating point amplitudes. The default value is 6 but the user\n",
    "        can decrease it to help recognize similar hopping amplitudes, particularly for symmetrical hoppings\n",
    "        across the crystal lattice\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    Hk: NumPy array\n",
    "        the Hamiltonian of the tight-binding model in momentum space. It can be output in either numerical\n",
    "        form (Hk_numerical) or reduced analytical form (Hk) based on the user's choice. The default output\n",
    "        is the reduced analytical form. The numerical form depends solely on the k-space vector components\n",
    "        while the analytical form takes into account both the k-space vector components and the lattice\n",
    "        vectors\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    import sympy as sp\n",
    "    \n",
    "    # imaginary number\n",
    "    I = sp.I\n",
    "\n",
    "    # matrix from axis directions in momentum space\n",
    "    kx, ky, kz = sp.symbols(\"kx ky kz\", real = True)\n",
    "    k_space_matrix = sp.Matrix([kx, ky, kz])\n",
    "\n",
    "    # symbolic dot product representation between lattice unit vectors\n",
    "    # and momentum space matrix\n",
    "    a1k, a2k, a3k = sp.symbols(\"a1k a2k a3k\", real = True)\n",
    "    lattice = sp.Matrix([a1k, a2k, a3k])\n",
    "\n",
    "    # units contains the displacement vectors\n",
    "    # hops contains details about hopping of electrons such as orbital\n",
    "    # and hopping amplitude\n",
    "    if TBL.units.shape == (2, 2):\n",
    "        TBL_units = np.eye(3)\n",
    "        TBL_units[:2, :2] = TBL.units\n",
    "        TBL_hops = {key + (0,): val for key, val in TBL.hoppings.items()}\n",
    "    elif TBL.units.shape == (3,3):\n",
    "        TBL_units = TBL.units\n",
    "        TBL_hops = TBL.hoppings\n",
    "    # raises error for when the dimensions of the tb object is neither 2D nor 3D\n",
    "    else:\n",
    "        raise ValueError(\"This format of the tight-binding model is not implemented for this function.\")\n",
    "   \n",
    "    # number of orbitals involved in the unit cell\n",
    "    num_orb = TBL.n_orbitals\n",
    "\n",
    "    # maximum hopping distances of electrons in each direction\n",
    "    max_x, max_y, max_z = list(np.max(np.array(list(TBL_hops.keys())), axis = 0))\n",
    "\n",
    "    # number of cells involved in the hopping of electrons in each direction\n",
    "    num_cells_x, num_cells_y, num_cells_z = [2 * max_coord + 1 for max_coord in [max_x, max_y, max_z]]\n",
    "    \n",
    "    # real-space Hamiltonian\n",
    "    Hrij = np.zeros((num_cells_x, num_cells_y, num_cells_z, num_orb, num_orb), dtype = sp.exp)\n",
    "\n",
    "    # looping through hopping parameters of electrons involved in inter-orbital hoppings\n",
    "    for key, hopping in TBL_hops.items():\n",
    "        rx, ry, rz = key\n",
    "        hopping = np.around(hopping, precision)\n",
    "        Hrij[rx + max_x, ry + max_y, rz + max_z] = hopping\n",
    "\n",
    "    # basis of the exponential term in calculation of Hk\n",
    "    Hexp = np.empty_like(Hrij, dtype = sp.exp)\n",
    "\n",
    "    # perform Fourier transform\n",
    "    for xi, yi, zi in itp(range(num_cells_x), range(num_cells_y), range(num_cells_z)):\n",
    "        coefficients = np.array([xi - max_x, yi - max_y, zi - max_z])\n",
    "        r = lattice.dot(coefficients)\n",
    "        eikr = sp.exp(-I * r)\n",
    "        Hexp[xi, yi, zi, :, :] = eikr\n",
    "\n",
    "    # summation over all real space axes\n",
    "    Hk = np.sum(Hrij * Hexp, axis = (0, 1, 2))\n",
    "    \n",
    "    # rewriting exponential terms in Hamiltonian expression in terms of cosine\n",
    "    for i, j in itp(range(num_orb), repeat = 2):\n",
    "        Hk[i, j] = Hk[i, j].rewrite(sp.cos)\n",
    "\n",
    "    def _has_complex_exponential_sympy(matrix):\n",
    "        \"\"\"\n",
    "        Checks if a NumPy array containing SymPy elements has a complex exponential element.\n",
    "\n",
    "        Args:\n",
    "            matrix (NumPy array): The input NumPy array containing SymPy elements\n",
    "        \n",
    "        Returns:\n",
    "            bool: True if the matrix array contains a complex exponential element, False otherwise.\n",
    "\n",
    "        \"\"\"\n",
    "\n",
    "        for sublist in matrix:\n",
    "            for element in sublist:\n",
    "                if element.is_complex and element.has(sp.exp):\n",
    "                    return True\n",
    "        return False\n",
    "    \n",
    "    def _is_hermitian_sympy(matrix):\n",
    "        \"\"\"\n",
    "        Checks if a NumPy array containing SymPy elements is hermitian\n",
    "\n",
    "        Args:\n",
    "            matrix (NumPy array): The input NumPy array containing SymPy elements\n",
    "        \n",
    "        Returns:\n",
    "            bool: True if the matrix is a hermitian, False otherwise\n",
    "\n",
    "        \"\"\"\n",
    "        \n",
    "        n = matrix.shape[0]\n",
    "        for i in range(n):\n",
    "            for j in range(n):\n",
    "                if matrix[i, j] != matrix[j, i].conjugate():\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    # performing the check on the analytical Hamiltonian\n",
    "    if not _is_hermitian_sympy(Hk): return warnings.warn(\"The resulting Hamiltonian is not hermitian.\")\n",
    "    if _has_complex_exponential_sympy(Hk): return warnings.warn(\"\"\"Your expression has a complex exponential. \n",
    "                                                                    Choosing a different unit cell could make \n",
    "                                                                    your Hamiltonian expression real.\"\"\")\n",
    "    \n",
    "    if analytical: return Hk\n",
    "\n",
    "    # dealing with the numerical Hamiltonian\n",
    "\n",
    "    # convert to SymPy matrix to use substitutions method available in SymPy\n",
    "    Hk_numerical = sp.Matrix(Hk)\n",
    "\n",
    "    # obtaining individual displacement vectors\n",
    "    a1, a2, a3 = np.around(TBL_units, precision)\n",
    "\n",
    "    # numerical dot products between unit vectors\n",
    "    # and momentum space matrix\n",
    "    a1k_n = a1.dot(k_space_matrix)[0]\n",
    "    a2k_n = a2.dot(k_space_matrix)[0]\n",
    "    a3k_n = a3.dot(k_space_matrix)[0]\n",
    "    \n",
    "    # performing numerical dot product substitutions\n",
    "    Hk_numerical = Hk_numerical.subs(a1k, a1k_n)\n",
    "    Hk_numerical = Hk_numerical.subs(a2k, a2k_n)\n",
    "    Hk_numerical = Hk_numerical.subs(a3k, a3k_n)\n",
    "\n",
    "    # converting numerical Hamiltonian to NumPy array\n",
    "    Hk_numerical = np.array(Hk_numerical)\n",
    "\n",
    "    if not _is_hermitian_sympy(Hk_numerical): return warnings.warn(\"The resulting Hamiltonian is not hermitian.\")\n",
    "    if _has_complex_exponential_sympy(Hk_numerical): return warnings.warn(\"\"\"Your expression has a complex exponential. \n",
    "                                                                            Choosing a different unit cell could make \n",
    "                                                                            your Hamiltonian expression real.\"\"\")\n",
    "    return Hk_numerical\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "96ab7f9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from triqs.lattice.utils import TB_from_pythTB\n",
    "w90_input = w90('SrCaCu2O4 Files', 'sccuo2')\n",
    "fermi_ev = 9.6082\n",
    "w90_model = w90_input.model(zero_energy = fermi_ev, min_hopping_norm = 0.05, max_distance = None)\n",
    "w90_triqs = TB_from_pythTB(w90_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f024d970",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The analytical expression \n",
      " [[-1.002*cos(a3k) - 0.124*cos(2*a3k) + 0.018\n",
      "  -0.056*I*sin(a1k) - 0.056*I*sin(a1k - a3k) + 0.056*cos(a1k) + 0.056*cos(a1k - a3k)]\n",
      " [0.056*I*sin(a1k) + 0.056*I*sin(a1k - a3k) + 0.056*cos(a1k) + 0.056*cos(a1k - a3k)\n",
      "  -0.904*cos(a3k) - 0.128*cos(2*a3k) - 0.13]]\n",
      "The numerical expression \n",
      " [[-1.002*cos(3.939*kz) - 0.124*cos(7.878*kz) + 0.018\n",
      "  -0.056*I*sin(1.762*kx - 7.996*ky) + 0.056*I*sin(-1.762*kx + 7.996*ky + 3.939*kz) + 0.056*cos(1.762*kx - 7.996*ky) + 0.056*cos(-1.762*kx + 7.996*ky + 3.939*kz)]\n",
      " [0.056*I*sin(1.762*kx - 7.996*ky) - 0.056*I*sin(-1.762*kx + 7.996*ky + 3.939*kz) + 0.056*cos(1.762*kx - 7.996*ky) + 0.056*cos(-1.762*kx + 7.996*ky + 3.939*kz)\n",
      "  -0.904*cos(3.939*kz) - 0.128*cos(7.878*kz) - 0.13]]\n"
     ]
    }
   ],
   "source": [
    "# calling the functions in different scenarios\n",
    "print(\"The analytical expression \\n\", TB_to_sympy(w90_triqs, analytical = True, precision = 3))\n",
    "print(\"The numerical expression \\n\", TB_to_sympy(w90_triqs, analytical = False, precision = 3))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "TRIQS",
   "language": "python",
   "name": "triqs"
  },
  "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}