{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"OMP_NUM_THREADS\"] = \"4\" # export OMP_NUM_THREADS=4\n",
    "os.environ[\"OPENBLAS_NUM_THREADS\"] = \"4\" # export OPENBLAS_NUM_THREADS=4 \n",
    "os.environ[\"MKL_NUM_THREADS\"] = \"4\" # export MKL_NUM_THREADS=6\n",
    "os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"4\" # export VECLIB_MAXIMUM_THREADS=4\n",
    "os.environ[\"NUMEXPR_NUM_THREADS\"] = \"4\" # export NUMEXPR_NUM_THREADS=6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import eigh\n",
    "#import matplotlib.pyplot as plt\n",
    "import warnings\n",
    "import sisl\n",
    "from grogu.useful import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.15.1'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.set_printoptions(linewidth=2*75)\n",
    "sisl.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this cell mimicks an input file\n",
    "fdf = sisl.get_sile('./Fe3GeTe2/Fe3GeTe2.fdf')\n",
    "# this information needs to be given at the input!!\n",
    "scf_xcf_orientation=np.array([0,0,1])                #z\n",
    "# list of reference directions for around which we calculate the derivatives\n",
    "# o is the quantization axis, v and w are two axes perpendicular to it\n",
    "# at this moment the user has to supply o,v,w on the input. \n",
    "# we can have some default sfor this\n",
    "ref_xcf_orientations=[dict(o=np.array([1,0,0]),vw=[np.array([0,1,0]),np.array([0,0,1])]),\n",
    "                      dict(o=np.array([0,1,0]),vw=[np.array([1,0,0]),np.array([0,0,1])]),\n",
    "                      dict(o=np.array([0,0,1]),vw=[np.array([1,0,0]),np.array([0,1,0])]),]\n",
    "\n",
    "# human readable definition of magnetic entities\n",
    "magnetic_entities=[dict(atom=3    ,l=2),\n",
    "                   dict(atom=4    ,l=2),\n",
    "                   dict(atom=5    ,l=2),\n",
    "                   dict(atom=[3,4],)]\n",
    "\n",
    "# pair information\n",
    "pairs=[dict(ai=0,aj=1,Ruc=np.array([0,0,0])),     # isotropic should be -82 meV\n",
    "       dict(ai=0,aj=2,Ruc=np.array([0,0,0])),     # these should all be around -41.9 in the isotropic part\n",
    "       dict(ai=1,aj=2,Ruc=np.array([0,0,0])),\n",
    "       dict(ai=0,aj=2,Ruc=np.array([-1,0,0])),\n",
    "       dict(ai=1,aj=2,Ruc=np.array([-1,0,0]))]\n",
    "\n",
    "# Bz sampling\n",
    "nk=np.array([10,10,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[5], line 3\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# digestion of the input\u001b[39;00m\n\u001b[1;32m      2\u001b[0m \u001b[38;5;66;03m# read in geometry and hamiltonian\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m geo \u001b[38;5;241m=\u001b[39m \u001b[43mfdf\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mread_geometry\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m      4\u001b[0m dh \u001b[38;5;241m=\u001b[39m fdf\u001b[38;5;241m.\u001b[39mread_hamiltonian()\n\u001b[1;32m      6\u001b[0m uc_in_sc_idx\u001b[38;5;241m=\u001b[39mdh\u001b[38;5;241m.\u001b[39mlattice\u001b[38;5;241m.\u001b[39msc_index([\u001b[38;5;241m0\u001b[39m,\u001b[38;5;241m0\u001b[39m,\u001b[38;5;241m0\u001b[39m])\n",
      "File \u001b[0;32m~/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.12/site-packages/sisl/io/siesta/fdf.py:1440\u001b[0m, in \u001b[0;36mfdfSileSiesta.read_geometry\u001b[0;34m(self, output, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1434\u001b[0m order \u001b[38;5;241m=\u001b[39m _parse_order(\n\u001b[1;32m   1435\u001b[0m     kwargs\u001b[38;5;241m.\u001b[39mpop(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124morder\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m),\n\u001b[1;32m   1436\u001b[0m     {\u001b[38;5;28;01mTrue\u001b[39;00m: [\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mXV\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mnc\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mTSHS\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSTRUCT\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mHSX\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfdf\u001b[39m\u001b[38;5;124m\"\u001b[39m], \u001b[38;5;28;01mFalse\u001b[39;00m: \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfdf\u001b[39m\u001b[38;5;124m\"\u001b[39m},\n\u001b[1;32m   1437\u001b[0m     output,\n\u001b[1;32m   1438\u001b[0m )\n\u001b[1;32m   1439\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m f \u001b[38;5;129;01min\u001b[39;00m order:\n\u001b[0;32m-> 1440\u001b[0m     v \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mgetattr\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43mf\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43m_r_geometry_\u001b[39;49m\u001b[38;5;132;43;01m{\u001b[39;49;00m\u001b[43mf\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlower\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;132;43;01m}\u001b[39;49;00m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1441\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m v \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m   1442\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m v\n",
      "File \u001b[0;32m~/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.12/site-packages/sisl/io/siesta/fdf.py:1511\u001b[0m, in \u001b[0;36mfdfSileSiesta._r_geometry_fdf\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1506\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_r_geometry_fdf\u001b[39m(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[1;32m   1507\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Returns Geometry object from the FDF file\u001b[39;00m\n\u001b[1;32m   1508\u001b[0m \n\u001b[1;32m   1509\u001b[0m \u001b[38;5;124;03m    NOTE: Interaction range of the Atoms are currently not read.\u001b[39;00m\n\u001b[1;32m   1510\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m-> 1511\u001b[0m     lattice \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mread_lattice\u001b[49m\u001b[43m(\u001b[49m\u001b[43morder\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mfdf\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1513\u001b[0m     \u001b[38;5;66;03m# No fractional coordinates\u001b[39;00m\n\u001b[1;32m   1514\u001b[0m     is_frac \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mFalse\u001b[39;00m\n",
      "File \u001b[0;32m~/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.12/site-packages/sisl/io/siesta/fdf.py:774\u001b[0m, in \u001b[0;36mfdfSileSiesta.read_lattice\u001b[0;34m(self, output, *args, **kwargs)\u001b[0m\n\u001b[1;32m    768\u001b[0m order \u001b[38;5;241m=\u001b[39m _parse_order(\n\u001b[1;32m    769\u001b[0m     kwargs\u001b[38;5;241m.\u001b[39mpop(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124morder\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m),\n\u001b[1;32m    770\u001b[0m     {\u001b[38;5;28;01mTrue\u001b[39;00m: [\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mXV\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mnc\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mTSHS\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSTRUCT\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mHSX\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfdf\u001b[39m\u001b[38;5;124m\"\u001b[39m], \u001b[38;5;28;01mFalse\u001b[39;00m: \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfdf\u001b[39m\u001b[38;5;124m\"\u001b[39m},\n\u001b[1;32m    771\u001b[0m     output,\n\u001b[1;32m    772\u001b[0m )\n\u001b[1;32m    773\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m f \u001b[38;5;129;01min\u001b[39;00m order:\n\u001b[0;32m--> 774\u001b[0m     v \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mgetattr\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43mf\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43m_r_lattice_\u001b[39;49m\u001b[38;5;132;43;01m{\u001b[39;49;00m\u001b[43mf\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlower\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;132;43;01m}\u001b[39;49;00m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    775\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m v \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m    776\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m v\n",
      "File \u001b[0;32m~/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.12/site-packages/sisl/io/siesta/fdf.py:781\u001b[0m, in \u001b[0;36mfdfSileSiesta._r_lattice_fdf\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    779\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_r_lattice_fdf\u001b[39m(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[1;32m    780\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Returns `Lattice` object from the FDF file\"\"\"\u001b[39;00m\n\u001b[0;32m--> 781\u001b[0m     s \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mLatticeConstant\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43munit\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mAng\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m    782\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m s \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m    783\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m SileError(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mCould not find LatticeConstant in file\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "File \u001b[0;32m~/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.12/site-packages/sisl/io/sile.py:707\u001b[0m, in \u001b[0;36msile_fh_open.<locals>._wrapper.<locals>.pre_open\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    705\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mhasattr\u001b[39m(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfh\u001b[39m\u001b[38;5;124m\"\u001b[39m):\n\u001b[1;32m    706\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m func(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n\u001b[0;32m--> 707\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28mself\u001b[39m:\n\u001b[1;32m    708\u001b[0m     reset(\u001b[38;5;28mself\u001b[39m)\n\u001b[1;32m    709\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m func(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n",
      "File \u001b[0;32m~/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.12/site-packages/sisl/io/sile.py:1107\u001b[0m, in \u001b[0;36mSile.__enter__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m   1105\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__enter__\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m   1106\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Opens the output file and returns it self\"\"\"\u001b[39;00m\n\u001b[0;32m-> 1107\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_open\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1108\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\n",
      "File \u001b[0;32m~/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.12/site-packages/sisl/io/sile.py:1098\u001b[0m, in \u001b[0;36mSile._open\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m   1096\u001b[0m             \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfh \u001b[38;5;241m=\u001b[39m gzip\u001b[38;5;241m.\u001b[39mopen(\u001b[38;5;28mstr\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfile), mode\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_mode)\n\u001b[1;32m   1097\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m-> 1098\u001b[0m         \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfh \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfile\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_mode\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1100\u001b[0m \u001b[38;5;66;03m# the file should restart the file-read (as per instructed)\u001b[39;00m\n\u001b[1;32m   1101\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_line \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0\u001b[39m\n",
      "File \u001b[0;32m/opt/homebrew/Cellar/python@3.12/3.12.7_1/Frameworks/Python.framework/Versions/3.12/lib/python3.12/pathlib.py:1013\u001b[0m, in \u001b[0;36mPath.open\u001b[0;34m(self, mode, buffering, encoding, errors, newline)\u001b[0m\n\u001b[1;32m   1011\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mb\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m mode:\n\u001b[1;32m   1012\u001b[0m     encoding \u001b[38;5;241m=\u001b[39m io\u001b[38;5;241m.\u001b[39mtext_encoding(encoding)\n\u001b[0;32m-> 1013\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mio\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mbuffering\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mencoding\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43merrors\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnewline\u001b[49m\u001b[43m)\u001b[49m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf'"
     ]
    }
   ],
   "source": [
    "# digestion of the input\n",
    "# read in geometry and hamiltonian\n",
    "geo = fdf.read_geometry()\n",
    "dh = fdf.read_hamiltonian()\n",
    "\n",
    "uc_in_sc_idx=dh.lattice.sc_index([0,0,0])\n",
    "\n",
    "# get indecies for orbitals of the magnetic entities\n",
    "for i,d in enumerate(magnetic_entities):\n",
    "    parsed = parse_magnetic_entity(dh,**d)\n",
    "    magnetic_entities[i]['orbital_indeces'] = parsed\n",
    "    magnetic_entities[i]['spin_box_indeces'] = blow_up_orbindx(parsed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WE WILL NOT NEED THIS!!\n",
    "eigfile=sisl.io.siesta.eigSileSiesta('./Fe3GeTe2/Fe3GeTe2.EIG')\n",
    "EF=eigfile.read_fermi_level()\n",
    "EF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "NO=dh.no # shorthand for number of orbitals in the unit cell\n",
    "# preprocessing Hamiltonian and overlap matrix elements\n",
    "h11  = dh.tocsr(dh.M11r)\n",
    "h11 += dh.tocsr(dh.M11i)*1.0j\n",
    "h11  = h11.toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128')\n",
    "\n",
    "h22  = dh.tocsr(dh.M22r)\n",
    "h22 += dh.tocsr(dh.M22i)*1.0j\n",
    "h22  = h22.toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128')\n",
    "\n",
    "h12  = dh.tocsr(dh.M12r)\n",
    "h12 += dh.tocsr(dh.M12i)*1.0j\n",
    "h12  = h12.toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128')\n",
    "\n",
    "h21  = dh.tocsr(dh.M21r)\n",
    "h21 += dh.tocsr(dh.M21i)*1.0j\n",
    "h21  = h21.toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128')\n",
    "\n",
    "sov = dh.tocsr(dh.S_idx).toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128')\n",
    "\n",
    "# Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation\n",
    "\n",
    "U=np.vstack([np.kron(np.eye(NO,dtype=int),[1,0]),np.kron(np.eye(NO,dtype=int),[0,1])])\n",
    "# This is the permutation that transforms ud1ud2 to u12d12\n",
    "# That is this transforms FROM SPIN BOX to ORBITAL BOX => U\n",
    "\n",
    "# the inverse transformation is U.T u12d12 to ud1ud2\n",
    "# That is FROM ORBITAL BOX to SPIN BOX => U.T\n",
    "\n",
    "# this is a test\n",
    "# u12d12=np.array([np.kron([1,1],np.arange(NO)),np.kron([0,1],np.ones(NO))],dtype=int).T\n",
    "# ud1ud2=np.array([np.kron(np.arange(NO),[1,1]),np.kron(np.ones(NO),[0,1])],dtype=int).T\n",
    "# np.allclose(u12d12,U@ud1ud2)\n",
    "\n",
    "hh,ss = \\\n",
    "np.array([U.T@np.block([[h11[:,:,i],h12[:,:,i]],\n",
    "                        [h21[:,:,i],h22[:,:,i]]])@U for i in range(dh.lattice.nsc.prod())]),      \\\n",
    "np.array([U.T@np.block([[sov[:,:,i]  ,sov[:,:,i]*0],\n",
    "                        [sov[:,:,i]*0,sov[:,:,i]  ]])@U for i in range(dh.lattice.nsc.prod())])\n",
    "\n",
    "\n",
    "# symmetrizing Hamiltonian and overlap matrix to make them hermitian \n",
    "for i in range(dh.lattice.sc_off.shape[0]):\n",
    "    j = dh.lattice.sc_index(-dh.lattice.sc_off[i])\n",
    "    h1,h1d=hh[i],hh[j]\n",
    "    hh[i],hh[j]=(h1+h1d.T.conj())/2,(h1d+h1.T.conj())/2\n",
    "    s1,s1d=ss[i],ss[j]\n",
    "    ss[i],ss[j]=(s1+s1d.T.conj())/2,(s1d+s1.T.conj())/2    \n",
    "\n",
    "# identifying TRS and TRB parts of the Hamiltonian\n",
    "TAUY = np.kron(np.eye(NO),tau_y)\n",
    "hTR  = np.array([TAUY@hh[i].conj()@TAUY for i in range(dh.lattice.nsc.prod())])\n",
    "hTRS = (hh+hTR)/2\n",
    "hTRB = (hh-hTR)/2\n",
    "\n",
    "# extracting the exchange field\n",
    "traced=[spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())]\n",
    "\n",
    "XCF = np.array([np.array([f['x'] for f in traced]),\n",
    "                np.array([f['y'] for f in traced]),\n",
    "                np.array([f['z'] for f in traced])])\n",
    "# Check if exchange field has scalar part\n",
    "max_xcfs=abs(np.array(np.array([f['c'] for f in traced]))).max()\n",
    "if  max_xcfs > 1e-12:\n",
    "    warnings.warn(f\"Exchange field has non negligeble scalr part. Largest value is {max_xcfs}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# collecting information for all proposed rotations, this might generate a lot of  unnceessary data\n",
    "# information beeing collected: \n",
    "# - matrices rotating the direction of the exchange field to the proposed reference direction\n",
    "# - rotated xc field \n",
    "# - single commutators of the unit cell localized exchange field for the two variations perpendicular to the reference direction\n",
    "# - double commutators of the unit cell localized exchange field for the two variations perpendicular to the reference direction\n",
    "\n",
    "\n",
    "# Transformation matrices for rotating the exchange field\n",
    "for i,orient in enumerate(ref_xcf_orientations):\n",
    "\n",
    "    # obtain rotated exchange field    \n",
    "    R=RotMa2b(scf_xcf_orientation,orient['o'])\n",
    "    \n",
    "    rot_XCF = np.einsum('ij,jklm->iklm',R, XCF)\n",
    "    rot_H_XCF = sum([np.kron(rot_XCF[i],tau) for i,tau in enumerate([tau_x,tau_y,tau_z])])\n",
    "    rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx]\n",
    "    \n",
    "    # obtain total Hamiltonian with the rotated exchange field\n",
    "    rot_H = hTRS+rot_H_XCF\n",
    "    \n",
    "    # Update magnetic entities with commutator data for each \n",
    "    for mag_ent in magnetic_entities:\n",
    "        mag_ent['Commutator_Data'] = []\n",
    "\n",
    "    for u in orient['vw']:\n",
    "        \n",
    "        Tu = np.kron(np.eye(NO,dtype=int),tau_u(u))\n",
    "\n",
    "        Vu1 = 1j/2*comm(rot_H_XCF_uc,Tu)\n",
    "        Vu2 = 1/8*comm(comm(Tu,rot_H_XCF_uc),Tu)\n",
    "\n",
    "        for mag_ent in magnetic_entities:\n",
    "            idx = mag_ent['spin_box_indeces']\n",
    "            mag_ent['Commutator_Data'].append(dict(Vu1=Vu1[idx,:][:,idx],\n",
    "                                                   Vu2=Vu2[idx,:][:,idx]))\n",
    "\n",
    "\n",
    "# TODO STORE ABOVE INFO!!!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k  = np.array([0.11441,0.234432,0.0])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This goes inside the for loop for k \n",
    "k  = np.array([0,0.1,0.0])\n",
    "\n",
    "\n",
    "phases = np.exp(-2j*np.pi*k@dh.lattice.sc_off.T)\n",
    "\n",
    "np.exp(-1j * np.dot(np.dot(np.dot(dh.rcell, k), dh.cell),dh.lattice.sc_off.T))\n",
    "ss.shape,rot_H_XCF.shape,hTRS.shape,phases.shape\n",
    "\n",
    "HK=np.einsum('abc,a->bc',hTRS+rot_H_XCF,phases)\n",
    "SK=np.einsum('abc,a->bc',ss,phases)\n",
    "\n",
    "vals,vecs=eigh(HK,SK)\n",
    "\n",
    "#This goes inside the for loop for ze\n",
    "ze = -.1+0.000j\n",
    "GK=vecs@np.diag(1/(ze-vals))@vecs.conj().T\n",
    "\n",
    "#phase=np.exp(1j*np.dot(np.dot(k,dh.rcell),offset))\n",
    "Gii=GK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]\n",
    "Gjj=GK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]\n",
    "Gij=GK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]*np.exp(-k@offset*2*np.pi)\n",
    "Gji=GK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]*np.exp( k@offset*2*np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.allclose(np.linalg.inv(ze*SK-HK),GK), abs(np.linalg.inv(ze*SK-HK)-GK).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hsk(dh,k=(0,0,0)):\n",
    "    '''\n",
    "    One way to speed up Hk and Sk generation\n",
    "    '''\n",
    "    k = np.asarray(k, np.float64) # this two conversion lines\n",
    "    k.shape = (-1,)               # are from the sisl source\n",
    "\n",
    "    # this generates the list of phases\n",
    "    phases = np.exp(-1j * np.dot(np.dot(np.dot(dh.rcell, k), dh.cell),\n",
    "                                 dh.lattice.sc_off.T))\n",
    "\n",
    "    HK11 = np.einsum('abc,c->ab',dh.h11,phases)\n",
    "    HK12 = np.einsum('abc,c->ab',dh.h12,phases)\n",
    "    HK21 = np.einsum('abc,c->ab',dh.h21,phases)\n",
    "    HK22 = np.einsum('abc,c->ab',dh.h22,phases)\n",
    "\n",
    "    SK  = np.einsum('abc,c->ab',dh.sov,phases)\n",
    "\n",
    "    return HK11,HK12,HK21,HK22,SK\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "    \n",
    "\n",
    "rot_xcf_mats=[]\n",
    "for i,rot in enumerate(rots):\n",
    "    R=RotM(*rot)             # rotation matrix to go from scf to reference direction\n",
    "    R@scf_xcf_orientation    # exchange field orientation \n",
    "\n",
    "\n",
    "rot_xcf_mats = [RotM(*rot) for rot in rots] \n",
    "# directions to ratate around \n",
    "dirs_xcf = [R@scf_xcf_orientation for R in rot_xcf_mats]\n",
    "dirs_deriv=[]\n",
    "for o in dirs_xcf:\n",
    "    x,y,z = np.eye(3)\n",
    "    R = RotMa2b(z,o)\n",
    "    dirs_deriv.append([R@x,R@y])\n",
    "\n",
    "\n",
    "# geting the exchange field part of the reference Hamiltonians for the rotated exchange field\n",
    "dh.hxc_rots=[]\n",
    "for mtx in rot_xcf_mats:\n",
    "    rotated_xcf = np.einsum('ij,jklm->iklm',mtx, XCF)\n",
    "    dh.hxc_rots.append(sum([np.kron(rotated_xcf[i],tau) \n",
    "           for i,tau in enumerate([tau_x,tau_y,tau_z]) ]))\n",
    "\n",
    "\n",
    "\n",
    "cummutators=[]\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i=0\n",
    "\n",
    "dh.hxc_rots[i][uc_in_sc_idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dh.hTRS.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.kron(np.random.rand(3,3),tau_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dh.hh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sanity check for the exchange field orientations\n",
    "xcf_test=np.array([[abs(xcfield[i].c).max(),\n",
    "                    abs(xcfield[i].x).max(),\n",
    "                    abs(xcfield[i].y).max(),\n",
    "                    abs(xcfield[i].z).max()] for i in range(dh.lattice.nsc.prod())])\n",
    "# Check how the exchange field components behave\n",
    "plt.plot(np.linalg.norm(dh.lattice.sc_off,axis=1),xcf_test[:,0],'o',label='scalar')\n",
    "plt.plot(np.linalg.norm(dh.lattice.sc_off,axis=1),xcf_test[:,1],'o',label='x')\n",
    "plt.plot(np.linalg.norm(dh.lattice.sc_off,axis=1),xcf_test[:,2],'o',label='y')\n",
    "plt.plot(np.linalg.norm(dh.lattice.sc_off,axis=1),xcf_test[:,3],'o',label='z')\n",
    "plt.yscale('log')\n",
    "plt.ylim(1e-22,1e1)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i=41\n",
    "DAT=np.vstack([np.hstack([dh.h11[:,:,i]+EF*dh.sov[:,:,i],dh.h12[:,:,i]]),\n",
    "               np.hstack([dh.h21[:,:,i],                 dh.h22[:,:,i]]+EF*dh.sov[:,:,i])])\n",
    "dh.lattice.sc_off[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i,v in enumerate(dh.lattice.sc_off):\n",
    "    if np.allclose(v,np.array([1,0,0])):\n",
    "        print(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.pcolor((U.T@DAT@U)[:6,:6].imag)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dh.h11  = dh.tocsr(dh.M11r)\n",
    "dh.h11 += dh.tocsr(dh.M11i)*1.0j\n",
    "dh.h11  = dh.h11.toarray().astype('complex128')\n",
    "dh.h22  = dh.tocsr(dh.M22r)\n",
    "dh.h22 += dh.tocsr(dh.M22i)*1.0j\n",
    "dh.h22  = dh.h22.toarray().astype('complex128')\n",
    "\n",
    "dh.h12  = dh.tocsr(dh.M12r)\n",
    "dh.h12 += dh.tocsr(dh.M12i)*1.0j\n",
    "dh.h12  = dh.h12.toarray().astype('complex128')\n",
    "dh.h21  = dh.tocsr(dh.M21r)\n",
    "dh.h21 += dh.tocsr(dh.M21i)*1.0j\n",
    "dh.h21  = dh.h21.toarray().reshape(dh.no,dh.n_s,dh.no).transpose(0,2,1).astype('complex128')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dh.tocsr(dh.S_idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hsk(dh,k=(0,0,0)):\n",
    "    '''\n",
    "    One way to speed up Hk and Sk generation\n",
    "    '''\n",
    "    k = np.asarray(k, np.float64) # this two conversion lines\n",
    "    k.shape = (-1,)               # are from the sisl source\n",
    "\n",
    "    # this generates the list of phases\n",
    "    phases = np.exp(-1j * np.dot(np.dot(np.dot(dh.rcell, k), dh.cell),\n",
    "                                 dh.lattice.sc_off.T))\n",
    "\n",
    "    HK11 = np.einsum('abc,c->ab',dh.h11,phases)\n",
    "    HK12 = np.einsum('abc,c->ab',dh.h12,phases)\n",
    "    HK21 = np.einsum('abc,c->ab',dh.h21,phases)\n",
    "    HK22 = np.einsum('abc,c->ab',dh.h22,phases)\n",
    "\n",
    "    SK  = np.einsum('abc,c->ab',dh.sov,phases)\n",
    "\n",
    "    return HK11,HK12,HK21,HK22,SK\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myHk=(lambda x:np.vstack([np.hstack([x[0],x[1]]),\n",
    "                          np.hstack([x[2],x[3]])]))(hsk(dh))\n",
    "mySk=(lambda x:np.vstack([np.hstack([x[-1],0*x[-1]]),\n",
    "                          np.hstack([0*x[-1],x[-1]])]))(hsk(dh))\n",
    "\n",
    "myHkshift=myHk+EF*mySk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H0=np.vstack( [np.hstack([dh.h11[:,:,0],dh.h12[:,:,0]]),\n",
    "               np.hstack([dh.h21[:,:,0],dh.h22[:,:,0]])])\n",
    "\n",
    "H0shifted=np.vstack( [np.hstack([dh.h11[:,:,0]+EF*dh.sov[:,:,0],dh.h12[:,:,0]]),\n",
    "                      np.hstack([dh.h21[:,:,0],                 dh.h22[:,:,0]+EF*dh.sov[:,:,0]])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(U.T@H0shifted@U)[:3,:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.pcolor((U.T@H0shifted@U)[:5,:5])\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.eigvalsh(myHk+EF*mySk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dh.Hk().toarray()[:3,:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.matshow(np.abs(dh.Hk().toarray()-myHk))\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(np.linalg.eigvalsh(dh.Hk().toarray())-np.linalg.eigvalsh(myHk))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "magnetic_entities = [parse_magnetic_entity(dh,atom=[0,1],l=2),parse_magnetic_entity(dh,atom=2,l=2)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO\n",
    "\n",
    "- Definition of magnetic entities:\n",
    "    * Through simple sequence o forbitals in the unit cell\n",
    "    * Through atom specification\n",
    "    * Through atom and orbital specification\n",
    "- Separation of TR and TRB components of the Hamiltonian, Identification of the exchange field. \n",
    "- Definition of commutator expressions, old projection matrix elements\n",
    "- Efficient calculation of Green's functions\n",
    "- Calculation of energy and momentum resolved derivatives\n",
    "- Parallel BZ and serial energy integral"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}