grogu/test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.14.3'"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tqdm import tqdm\n",
    "from sys import getsizeof\n",
    "from timeit import default_timer as timer\n",
    "\n",
    "import numpy as np\n",
    "import sisl\n",
    "import sisl.viz\n",
    "from src.grogu_magn.useful import *\n",
    "from mpi4py import MPI\n",
    "import pickle\n",
    "from numpy.linalg import inv\n",
    "import warnings\n",
    "\n",
    "\"\"\"    \n",
    "# Some input parsing\n",
    "parser = argparse.ArgumentParser()\n",
    "parser.add_argument('--kset'    , dest = 'kset'   , default  = 2         , type=int  , help = 'k-space resolution of Jij calculation')\n",
    "parser.add_argument('--kdirs'   , dest = 'kdirs'  , default  = 'xyz'                 , help = 'Definition of k-space dimensionality')\n",
    "parser.add_argument('--eset'    , dest = 'eset'   , default  = 42        , type=int  , help = 'Number of energy points on the contour')\n",
    "parser.add_argument('--eset-p'  , dest = 'esetp'  , default  = 10        , type=int  , help = 'Parameter tuning the distribution on the contour')\n",
    "parser.add_argument('--input'   , dest = 'infile' , required = True                  , help = 'Input file name')\n",
    "parser.add_argument('--output'  , dest = 'outfile', required = True                  , help = 'Output file name')\n",
    "parser.add_argument('--Ebot'    , dest = 'Ebot'   , default  = -20.0     , type=float, help = 'Bottom energy of the contour')\n",
    "parser.add_argument('--npairs'  , dest = 'npairs' , default  = 1         , type=int  , help = 'Number of unitcell pairs in each direction for Jij calculation')\n",
    "parser.add_argument('--adirs'   , dest = 'adirs'  , default  = False                 , help = 'Definition of pair directions')\n",
    "parser.add_argument('--use-tqdm', dest = 'usetqdm', default  = 'not'                 , help = 'Use tqdm for progressbars or not')\n",
    "parser.add_argument('--cutoff'  , dest = 'cutoff' , default  = 100.0     , type=float, help = 'Real space cutoff for pair generation in Angs')\n",
    "parser.add_argument('--pairfile', dest = 'pairfile', default  = False    ,             help = 'File to read pair information')\n",
    "args = parser.parse_args()\n",
    "\"\"\"\n",
    "# runtime information\n",
    "times = dict()\n",
    "times[\"start_time\"] = timer()\n",
    "########################\n",
    "# it works if data is in downloads folder\n",
    "########################\n",
    "sisl.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-12.806878959999999"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dat = sisl.io.siesta.eigSileSiesta(\n",
    "    \"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.EIG\"\n",
    ")\n",
    "siesta_eigs = dat.read_data()\n",
    "siesta_eigs.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-5.82448514"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ef = dat.read_fermi_level()\n",
    "Ef"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================================================================================================================================================================\n",
      "Input file: \n",
      "/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\n",
      "Output file: \n",
      "./Fe3GeTe2.pickle\n",
      "Number of nodes in the parallel cluster:  1\n",
      "================================================================================================================================================================\n",
      "Cell [Ang]: \n",
      "[[ 3.79100000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [-1.89550000e+00  3.28310231e+00  0.00000000e+00]\n",
      " [ 1.25954923e-15  2.18160327e-15  2.05700000e+01]]\n",
      "================================================================================================================================================================\n",
      "DFT axis: \n",
      "[0 0 1]\n",
      "Quantization axis and perpendicular rotation directions:\n",
      "[1 0 0]  --»  [array([0, 1, 0]), array([0, 0, 1])]\n",
      "[0 1 0]  --»  [array([1, 0, 0]), array([0, 0, 1])]\n",
      "[0 0 1]  --»  [array([1, 0, 0]), array([0, 1, 0])]\n",
      "================================================================================================================================================================\n",
      "Parameters for the contour integral:\n",
      "Number of k points:  10\n",
      "k point directions:  xy\n",
      "Ebot:  -13\n",
      "Eset:  600\n",
      "Esetp:  10000\n",
      "================================================================================================================================================================\n",
      "Setup done. Elapsed time: 3484.176097333 s\n",
      "================================================================================================================================================================\n"
     ]
    }
   ],
   "source": [
    "################################################################################\n",
    "#################################### INPUT #####################################\n",
    "################################################################################\n",
    "path = (\n",
    "    \"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\"\n",
    ")\n",
    "outfile = \"./Fe3GeTe2\"\n",
    "\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 for this\n",
    "ref_xcf_orientations = [\n",
    "    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",
    "magnetic_entities = [\n",
    "    dict(atom=3, l=2),\n",
    "    dict(atom=4),\n",
    "    dict(atom=5, l=2),\n",
    "]\n",
    "pairs = [\n",
    "    dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])),\n",
    "    dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])),\n",
    "    dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])),\n",
    "    dict(ai=0, aj=2, Ruc=np.array([-1, -1, 0])),\n",
    "    dict(ai=1, aj=2, Ruc=np.array([-1, -1, 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",
    "# Brilloun zone sampling and Green function contour integral\n",
    "kset = 10\n",
    "kdirs = \"xy\"\n",
    "ebot = -13\n",
    "eset = 600\n",
    "esetp = 10000\n",
    "################################################################################\n",
    "#################################### INPUT #####################################\n",
    "################################################################################\n",
    "\n",
    "# MPI parameters\n",
    "comm = MPI.COMM_WORLD\n",
    "size = comm.Get_size()\n",
    "rank = comm.Get_rank()\n",
    "root_node = 0\n",
    "\n",
    "# rename outfile\n",
    "if not outfile.endswith(\".pickle\"):\n",
    "    outfile += \".pickle\"\n",
    "\n",
    "simulation_parameters = dict(\n",
    "    path=path,\n",
    "    outpath=outfile,\n",
    "    scf_xcf_orientation=scf_xcf_orientation,\n",
    "    ref_xcf_orientations=ref_xcf_orientations,\n",
    "    kset=kset,\n",
    "    kdirs=kdirs,\n",
    "    ebot=ebot,\n",
    "    eset=eset,\n",
    "    esetp=esetp,\n",
    "    parallel_size=size,\n",
    ")\n",
    "\n",
    "# digestion of the input\n",
    "# read sile\n",
    "fdf = sisl.get_sile(path)\n",
    "# read in hamiltonian\n",
    "dh = fdf.read_hamiltonian()\n",
    "simulation_parameters[\"cell\"] = fdf.read_geometry().cell\n",
    "\n",
    "# unit cell index\n",
    "uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0])\n",
    "\n",
    "if rank == root_node:\n",
    "    print_parameters(simulation_parameters)\n",
    "    times[\"setup_time\"] = timer()\n",
    "    print(f\"Setup done. Elapsed time: {times['setup_time']} s\")\n",
    "    print(\n",
    "        \"================================================================================================================================================================\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xyz[-3:]: red, green, blue\n",
      "2.745163300331324\n",
      "2.5835033632437767\n",
      "2.583501767937866\n",
      "2.583541444641373\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABMUAAAGyCAYAAADksXO/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAADgcklEQVR4nOy9eZwlRZnv/YuIzDxLrV29VrO0jaDQgM02zbA06CjQqDCojDPOVdH3qiOv3utcxg1HgXZmRMbl4o7eOzNc5M47Mt6LwqityOggwjQqimLbCk7TslTvXetZMjMi3j8iM09mnjynzqk6VaeW5/v5BBHxxJqnm66s34l4Hqa11iAIgiAIgiAIgiAIgiCIZQTv9gYIgiAIgiAIgiAIgiAIYr4hUYwgCIIgCIIgCIIgCIJYdpAoRhAEQRAEQRAEQRAEQSw7SBQjCIIgCIIgCIIgCIIglh0kihEEQRAEQRAEQRAEQRDLDhLFCIIgCIIgCIIgCIIgiGUHiWIEQRAEQRAEQRAEQRDEsoNEMYIgCIIgCIIgCIIgCGLZQaIYQRAEQRAEQRAEQRAEsewgUYwgCIIgCIIgCIIgCIJYdsypKPbAAw/giiuuwPr168EYw9e+9rVEu9YaN9xwA4aHh1EoFPCyl70MTzzxxFxuiSAIgiAIgugA9J5HEARBEMRiZ05FsampKWzevBmf+9znMtv/9m//Fp/+9Kdx2223YefOnejp6cFll12GSqUyl9siCIIgCIIgZgm95xEEQRAEsdhhWms9LwsxhrvvvhtXXXUVAPPt4fr16/EXf/EXePe73w0AGBsbw9q1a3H77bfjT/7kT+ZjWwRBEARBEMQsofc8giAIgiAWI1a3Ft6zZw/27duHl73sZZFtYGAA5557Lh5++OGGL0vVahXVajWqK6Vw5MgRrFy5EoyxOd83QRAEQRCLH601JiYmsH79enBOLlY7Db3nEQRBEATRLdp5z+uaKLZv3z4AwNq1axP2tWvXRm1Z3Hzzzdi+ffuc7o0gCIIgiOXB008/jWOPPbbb21hy0HseQRAEQRDdppX3vK6JYjPl+uuvx3XXXRfVx8bGcPzxx2PPnj3o6+vr4s46h+d5+N73voeXvOQlsG2729uZV5bzswPL+/np2ZfnswPL+/kX9bMrCetLFwKT+zCT8zv+a/8R7vA5XXv+iYkJbNy4ccm8OywVGr3nPf300+jv7+/4euf89X2oeAo73rUVxw4VOz4/MUOUApQHSA+QLqBkraylyWXQrvygXAEq40BlFCiP1vLyWL0Ncnb723QV8Iefnd0cBEEQxJwxPj6O4447rqX3vK6JYuvWrQMA7N+/H8PDw5F9//79OOOMMxqOy+VyyOVydfahoaE5eVnqBp7noVgsYuXKlYvvl6RZspyfHVjez0/PvjyfHVjez7+on33PDwBvP5BrVxJjQP96YPPl8KTq2vOH69GVvLmh0+95/f39c/KeJ3I94Fyir78f/f1LUBTTOhCUqkZI8t1AUAqSX62JTdKNCUxhORCfIoEqq81PiVThfH6yLN3kPM3m1LMUrZqRA9BMymccKK4CelYB+QHAygN2EbDzgF0AVmwEznuHKRMEQRALmlbe87omim3cuBHr1q3D/fffH70cjY+PY+fOnbj22mu7tS2CIAiCIFphcv8MBgUvJts+CnABSNXRLRELh2X1nqdkTGQKBadAbIpEp2pN8PGryBamYgJVelxazGoqcFWT62BeYmrNPdwGhAMIy+TcBkSQ4mXhANxKlZ36vvlBoGe1Eb96VtdSYQVAfgYJgiCWDXMqik1OTuLJJ5+M6nv27MHPfvYzDA0N4fjjj8ef//mf46//+q9x0kknYePGjfjQhz6E9evXR5GLCIIgCIJYoPSunb5Pmv71RhDbdGXn90PMO0v2PU9JYOogMLHPpMl9tXK8Xh6tXedbTIgcYOUCkShXE49Cwalt8cmZuTAVlcM5svo5RkSnU50EQRDEHDCnotiPf/xjvOQlL4nqoY+Ia665Brfffjve+973YmpqCm9729swOjqKCy+8EDt27EA+n5/LbREEQRAEMVs2nG9ErvERZJ9EYUDfMPCq24zA0LvWjOFivndKzBFL4T1PQ+Ml/KdYde/tQOWAOQE5uR/QszjFyO2Y6OTUhKeEEOUAVihExVLa1taY+DoN5uQWiUsEQRAEEWNORbEXv/jF0LrxkW3GGD784Q/jwx/+8FxuAwAgpYTneXO+TifwPA+WZaFSqUDKhf/to23bEIJ+ySEIglhWcAFsuwW4640w1yLjP++DX7ovvwU44eIubI6YDxbSe95MuRCP4Tb7E7D2pEQwxoGeNUDfWiPu9q4F+taZ1BvkxaFssYpEJ4IgCIJYNCy66JPtorXGvn37MDo62u2ttIzWGuvWrcPTTz+9aBwADw4OYt26dYtmvwRBEEQH2HQl8No7gB3vA8afq9npmiSxSPhLfjssplA68ZUo/t4baiJYcZW5PkgQBEEQxJJmyf+0DwWxNWvWoFgsLgrRRimFyclJ9Pb2gi9wR59aa5RKJRw4cAAAEhGmCIIgiGXApiuBk18B7H3IXDuja5LEImII4wCAsd9/L4onbu7ybgiCIAiCmG+WtCgmpYwEsZUrV3Z7Oy2jlILrusjn8wteFAOAQsGEpD5w4ADWrFlDVykJgiCWG1wAG7d2excE0TYsvPbLFv77FkEQBEEQnWdJvwGEPsSKxWKXd7L0CT/jxeK3jSAIgiAI4jAGAAD9P/408Mu7gWcfBaYOA018pREEQRAEsXRY0ifFQhbDlcnFDn3GBEEQBEEsNj6nX4OPs8+i51d3Ab+6q9bg9AKDxwODG4L8eGDFBqBndeBYP1+L8mjlA0f7OWOjdyKCIAiCWDQsC1GMIAiCIIgOoCT5DiOWFN/EhWCei5tOPYie0rPA6O+AyX2AOwkc2GVSu4iYQNZIOLNyQTloE07K3qxPG/1JoCMIgiCIppAotsh405vehNHRUXzta1/raF+CIAiCaMquexpEmbyFokwSixYGhn+WL8Y7L3sxelb2GKNXBsaeAY7uBUb3GqFsdK+pl48C0gX8qkmyaupxpGuSOzH/D5QmFOiEHZQdE1VTOCmbDXA7w24l+4Rl3sDeVn871YcEdoIgCGL+IVGsBaTSeGTPERyYqGBNXx5bNg5B8O588/apT30KmvxcEARBEPPJrnuAu94IIPXzZ3zE2F97BwljxKIkdP+QeLWyC8Cqk0xqBa0DoawC+K4RyuKiWVQO+1Qz+gf1RNltME81Y53YmDihQLcYYDwQx1JiWVRuILhZOaD/GGBoI7DiecCKjea6q7C7/UQEQRDEIoBEsWnY8fgItt+7CyNjlcg2PJDHjVdswrbThud9PwMDA/O+JkEQBLGMSF+RPO5cc0IsLYgBgY0BO94PnPwKOulBLDrCrzhn9XUjY7Xri90mLdCFQpnyawKZ9GLlBvbM/l6qHGtXDex1/YNceWaNxN6V2S8qmY/WFowDA8fWRLJcH7D5dcC602Y/N0EQBLGkIFGsCTseH8G1dz5a96K0b6yCa+98FF94/VlzJox9/etfx8c//nE8+eSTKBaLOPPMM/H1r38d73jHOxJXIr/61a9i+/btdf16enrmZF8EQRDEEibrimRxJVA63GSQBsafNULaxq1zvkWC6CS2xYEq8KuRcTxvZXHxBw5aSALddCiVIaZl5CottKUEN3cKGHsaOLIHOPoUcPhJM2b0dybtecCs9+N/AD7wLPlZIwiCIBKQKNYAqTS237ur2ffi2H7vLlyyaV3Hr1KOjIzgLW95C2655Ra8+tWvxsTEBH7wgx/UXZscGRnB6173Ovzt3/4tXvWqVzXsRxAEQRDT0uiKZFNBLMbk/o5viSDmmledeQz+7sE9+H//96Mo2ALrBvJY25/D2v481vXnTT5g8rX9Oazpy8OxeLe3vfBQKnmFM3FSLX0lNH5ttNLkumiztoxyWNcqe48XvIsEMYIgCKIOEsUa8MieI4krk2k0gJGxCh7ZcwTnPX9lR9ceGRmB7/t41atehec973kAgNNPP71hv1e/+tXYsGFDw34
      "text/plain": [
       "<Figure size 1500x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABNoAAAHACAYAAAB0/gUQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJqUlEQVR4nO3de5xVdb0//tcAMnhjytABPJPXEq9AmoTZV7FJUvNIvzKzjhBftcvRfhmlSV6wUilNIxXlpJmpecR7nfSHEUf0dKQ8IpyjpR4vIHgZ1KPOACrIzP79MYfRkYszsGb2zPB8+lgPnM9aa+/3Z8Nen/15zdprVZRKpVIAAAAAgI3Sq9wFAAAAAEBPIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAoQJ9yF9AVNTU15fnnn8/WW2+dioqKcpcD0O2VSqUsXbo0gwcPTq9efsdjnAEolnFmTcYagGK1dawRtK3F888/n5qamnKXAdDjLF68OH/3d39X7jLKzjgD0DGMM28z1gB0jPcaawRta7H11lsnaX7x+vfvX+ZqALq/hoaG1NTUtBxfN3XGGYBiGWfWZKwBKFZbxxpB21qsPrW6f//+BiWAAvnqSjPjDEDHMM68zVgD0DHea6xxAQMAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKCtCEuWJL/61drXPfZY8tvfdm490N0sW5b84hfJ8ccnX/1qcsstyVtvlbsq6DqampJLLkleey259trkhBOSE09Mrr8+efXV5NJLm7cBAHqcFSuaPwY0Nq657tVXmz9GA+u2alVy223J177WPOWcNi1ZurTjnq9Pxz30JmLZsuSQQ5K//S155ZXkO995e91jjyUHH5y8/HLy+98nn/502cqELuuee5IxY5qPdL17N7ddeWWyww7J3Xcnu+1W1vKgS/jmN5PLL09OPTVZuTLp87/D91VXJX37Nrc9/XTys5+Vt04AoFClUvKlLzWHBA8+2Hx+x+qPzK++mnzqU8ncuc2/izvttLKWCl3Sk08mo0c3f1Re/RH6V79q/lh9663JoYcW/5zOaNtYW26ZfP7zzf//3e8mF13U/P+rQ7YlS5I990z2269sJUKX9cQTyeGHNwfWpVLzrxpWrWpe9+yzzSH2smXlrRG6gkMOaf5z5crmP9/5XlndtnobAKDHqKhoDtp6906uuy4ZP775zLZ3hmwDBjR/pAZae/315o/IzzzT/PPqj9ClUvO6I49sPmeqaIK2jVVRkZxzTnL22c0/f/e7zV/nWR2y7bNPMmtW89EPaO3nP28+0q3tK2+NjckLLyS/+U3n1wVdzX/9V9JrPUN2RUXy8MOdVw8A0Gk+97lk+vS3w7bPfjb55CffDtnuuSfZa69yVwldz/TpyeLFa//adVNT8/Lznxf/vIK2Irw7bLvqKiEbtMXNN799Vs663Hpr59QCXdnNN6//GmylUvM2AECPtDpsS5J/+Zdk3rxkm22EbLA+t922/t9Vr1qV3HRT8c8raCtKRUVy7LGt2444QsgG6/PGG+tfXyoly5d3Ti3QlbXlfeC9AgA92ruvErHDDsnuu5enFugOli9/7/uFvflm8c8raCvK6muyvdPkyW9fsw1Y0z77rP9XDH36JMOGdVo50GUNH/72lY/Xpnfv5m0AgB5p9TXZ3mnevLev2QasaejQt2+AsDa9enXMGaGCtiK888YH++yTvPRS62u2Cdtg7U4+ef2/Yli1qvkezLCpO/749X+Kbmxs3gYA6HHefeODhx9ObrllzRskAK197Wvrv1JRU1PzlLRograNtWxZ85Uo331NtnffIOGWW8paJnRJX/hC8sUvNn/1uqLi7fbVZ7mdd17z+wo2dTfeuPb2d75v/vmfO6cWAKDTlErN12d7940P3n2DhLPOKnel0PUMGZJccEHz/7/zi1Srp5//z/+T/MM/FP+8graNtdVWzWHAsGGtb3zwzhskjB6dfOYz5awSuqZevZLrr08uuyzZZZe32z/60eabIHz/++WrDbqSM85ofo9MnJjsuefb7Xvv3fw+2XXX5nUAQI9SUdE83dxllzVvfLA6bNtzz445Kwd6glNPTX772+RjH3u7baedmu82etNN6786y4aqKJVKpeIftntraGhIVVVV6uvr079//7bt9NZbyWabrX3dqlXr/2Iw0PzruoaG5iPdVluVuxoKtkHH1R5sg16Pd44lDQ3Nn7y33nrNdQCbIOPMmrwmPcv6hnofA6Btli9vfr/079/6iyFt1dbjqrdjUdYVsiWOetAWFRVJVVW5q4Cu651jybsHduMMAPRo6xvqfQyAttlyy855Hl8dBQAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAKUNWi77777cuSRR2bw4MGpqKjIHXfc0Wr9bbfdlkMPPTQf+MAHUlFRkfnz57/nY15zzTWpqKhotfTr169jOgBAl2acAQAAOlNZg7bly5dn6NChmTp16jrXH3jggfnJT37Srsft379/XnjhhZblmWeeKaJcALoZ4wwAANCZ+pTzyQ877LAcdthh61x/3HHHJUkWLlzYrsetqKjIwIEDN6Y0AHoA4wwAANCZeuQ12pYtW5YddtghNTU1Oeqoo/LXv/51vduvWLEiDQ0NrRYAWBfjDAAAsDY9LmjbbbfdcvXVV+e3v/1trr/++jQ1NeWAAw7Is88+u859Jk+enKqqqpalpqamEysGoDsxzgAAAOvS44K2kSNHZuzYsRk2bFgOOuig3Hbbbdl2223zT//0T+vcZ+LEiamvr29ZFi9e3IkVA9CdGGcAAIB1Kes12jrDZpttluHDh+fJJ59c5zaVlZWprKzsxKoA6CmMMwAAwGo97oy2d2tsbMzDDz+cQYMGlbsUAHog4wwAALBaWc9oW7ZsWaszABYsWJD58+dnm222yQc/+MG88sorWbRoUZ5//vkkyeOPP54kGThwYMvd3saOHZvtt98+kydPTpL88Ic/zMc+9rHsuuuuee2113LhhRfmmWeeyQknnNDJvQOg3IwzAABAZypr0Pbggw9m1KhRLT9PmDAhSTJu3Lhcc801+d3vfpfx48e3rP/iF7+YJJk0aVLOOeecJMmiRYvSq9fbJ+a9+uqrOfHEE1NXV5f3v//92XfffXP//fdnjz326IQeAdCVGGcAAIDOVFEqlUrlLqKraWhoSFVVVerr69O/f/9ylwPQ7Tmutub1AChWVz+u3nfffbnwwgszd+7cvPDCC7n99tszZsyY9e4ze/bsTJgwIX/9619TU1OTM888M1/5ylfa/Jxd/TUB6G7aelzt8ddoAwAAKKfly5dn6NChmTp1apu2X7BgQY444oiMGjUq8+fPzymnnJITTjghd999dwdXCsDG6vF3HQUAACinww47LIcddlibt582bVp22mmnXHTRRUmS3XffPX/605/ys5/9LKNHj+6oMgEogDPaAAAAupA5c+aktra2Vdvo0aMzZ86cde6zYsWKNDQ0tFoA6HyCNgAAgC6krq4u1dXVrdqqq6vT0NCQN954Y637TJ48OVVVVS1LTU1NZ5QKwLsI2gAAALq5iRMnpr6+vmVZvHhxuUsC2CS5RhsAAEAXMnDgwCxZsqRV25IlS9K/f/9svvnma92nsrIylZWVnVEeAOvhjDYAAIAuZOTIkZk1a1artpkzZ2bkyJFlqgiAthK0AQAAdKBly5Zl/vz5mT9/fpJkwYIFmT9/fhYtWpSk+WufY8eObdn+61//ep5++umcdtppeeyxx3L55Zfnpptuyre//e1ylA9AOwjaAAAAOtCDDz6Y4cOHZ/jw4UmSCRMmZPjw4Tn77LOTJC+88EJL6JYkO+20U+68887MnDkzQ4cOzUUXXZSrrroqo0ePLkv9ALSda7QBAAB0oIMPPjilUmmd66+55pq17jNv3rwOrAqAjuC
      "text/plain": [
       "<Figure size 1500x500 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(15, 5))\n",
    "plt.subplot(121)\n",
    "plt.plot(np.sort(dh.eigh()), marker=\"o\", linestyle=\" \", label=\"sisl\")\n",
    "plt.plot(np.sort(siesta_eigs[0, 0]), marker=\"o\", linestyle=\" \", label=\"siesta\")\n",
    "plt.ylim(None, 10)\n",
    "plt.xlim(None, 75)\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.subplot(122)\n",
    "DOS = sisl.physics.electron.DOS(np.linspace(-15, 85, 1000), dh.eigh())\n",
    "plt.plot(DOS, np.linspace(-15, 85, 1000))\n",
    "DOS = sisl.physics.electron.DOS(np.linspace(-15, 85, 1000), siesta_eigs[0, 0])\n",
    "plt.plot(DOS, np.linspace(-15, 85, 1000))\n",
    "plt.ylim(None, 10)\n",
    "\n",
    "coords = dh.xyz[-3:]\n",
    "\n",
    "shift = np.array([-1, 0, 0]) @ simulation_parameters[\"cell\"]\n",
    "\n",
    "\n",
    "plt.figure(figsize=(15, 5))\n",
    "plt.subplot(131)\n",
    "plt.scatter(coords[:, 0], coords[:, 2], color=[\"r\", \"g\", \"b\"])\n",
    "plt.scatter(\n",
    "    (coords + shift)[:, 0], (coords + shift)[:, 2], color=[\"r\", \"g\", \"b\"], marker=\"x\"\n",
    ")\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"z\")\n",
    "plt.subplot(132)\n",
    "plt.scatter(coords[:, 1], coords[:, 2], color=[\"r\", \"g\", \"b\"])\n",
    "plt.scatter(\n",
    "    (coords + shift)[:, 1], (coords + shift)[:, 2], color=[\"r\", \"g\", \"b\"], marker=\"x\"\n",
    ")\n",
    "plt.xlabel(\"y\")\n",
    "plt.ylabel(\"z\")\n",
    "plt.subplot(133)\n",
    "plt.scatter(coords[:, 0], coords[:, 1], color=[\"r\", \"g\", \"b\"])\n",
    "plt.scatter(\n",
    "    (coords + shift)[:, 0], (coords + shift)[:, 1], color=[\"r\", \"g\", \"b\"], marker=\"x\"\n",
    ")\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")\n",
    "print(\"xyz[-3:]: red, green, blue\")\n",
    "\n",
    "print(np.linalg.norm(coords[0] - coords[1]))\n",
    "print(np.linalg.norm(coords[0] - coords[2]))\n",
    "print(np.linalg.norm(coords[2] - coords[1]))\n",
    "print(np.linalg.norm(coords[0] - (coords + shift)[2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hamiltonian and exchange field rotated. Elapsed time: 3484.986699166 s\n",
      "================================================================================================================================================================\n"
     ]
    }
   ],
   "source": [
    "NO = dh.no  # shorthand for number of orbitals in the unit cell\n",
    "\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 = (\n",
    "    dh.tocsr(dh.S_idx)\n",
    "    .toarray()\n",
    "    .reshape(NO, dh.n_s, NO)\n",
    "    .transpose(0, 2, 1)\n",
    "    .astype(\"complex128\")\n",
    ")\n",
    "\n",
    "\n",
    "# Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation\n",
    "U = np.vstack(\n",
    "    [np.kron(np.eye(NO, dtype=int), [1, 0]), np.kron(np.eye(NO, dtype=int), [0, 1])]\n",
    ")\n",
    "# This is the permutation that transforms ud1ud2 to u12d12\n",
    "# That is this transforms FROM SPIN BOX to ORBITAL BOX => U\n",
    "# the inverse transformation is U.T u12d12 to ud1ud2\n",
    "# That is FROM ORBITAL BOX to SPIN BOX => U.T\n",
    "\n",
    "# From now on everything is in SPIN BOX!!\n",
    "hh, ss = np.array(\n",
    "    [\n",
    "        U.T @ np.block([[h11[:, :, i], h12[:, :, i]], [h21[:, :, i], h22[:, :, i]]]) @ U\n",
    "        for i in range(dh.lattice.nsc.prod())\n",
    "    ]\n",
    "), np.array(\n",
    "    [\n",
    "        U.T\n",
    "        @ np.block([[sov[:, :, i], sov[:, :, i] * 0], [sov[:, :, i] * 0, sov[:, :, i]]])\n",
    "        @ U\n",
    "        for i in range(dh.lattice.nsc.prod())\n",
    "    ]\n",
    ")\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())]  # equation 77\n",
    "XCF = np.array(\n",
    "    [\n",
    "        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",
    "    ]\n",
    ")  # equation 77\n",
    "\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(\n",
    "        f\"Exchange field has non negligible scalar part. Largest value is {max_xcfs}\"\n",
    "    )\n",
    "\n",
    "if rank == root_node:\n",
    "    times[\"H_and_XCF_time\"] = timer()\n",
    "    print(\n",
    "        f\"Hamiltonian and exchange field rotated. Elapsed time: {times['H_and_XCF_time']} s\"\n",
    "    )\n",
    "    print(\n",
    "        \"================================================================================================================================================================\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Site and pair dictionaries created. Elapsed time: 3485.081068083 s\n",
      "================================================================================================================================================================\n"
     ]
    }
   ],
   "source": [
    "# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions\n",
    "for mag_ent in magnetic_entities:\n",
    "    parsed = parse_magnetic_entity(dh, **mag_ent)  # parse orbital indexes\n",
    "    mag_ent[\"orbital_indeces\"] = parsed\n",
    "    mag_ent[\"spin_box_indeces\"] = blow_up_orbindx(parsed)  # calculate spin box indexes\n",
    "    # if orbital is not set use all\n",
    "    if \"l\" not in mag_ent.keys():\n",
    "        mag_ent[\"l\"] = \"all\"\n",
    "    if isinstance(mag_ent[\"atom\"], int):\n",
    "        mag_ent[\"tags\"] = [\n",
    "            f\"[{mag_ent['atom']}]{dh.atoms[mag_ent['atom']].tag}({mag_ent['l']})\"\n",
    "        ]\n",
    "        mag_ent[\"xyz\"] = [dh.xyz[mag_ent[\"atom\"]]]\n",
    "    if isinstance(mag_ent[\"atom\"], list):\n",
    "        mag_ent[\"tags\"] = []\n",
    "        mag_ent[\"xyz\"] = []\n",
    "        # iterate over atoms\n",
    "        for atom_idx in mag_ent[\"atom\"]:\n",
    "            mag_ent[\"tags\"].append(\n",
    "                f\"[{atom_idx}]{dh.atoms[atom_idx].tag}({mag_ent['l']})\"\n",
    "            )\n",
    "            mag_ent[\"xyz\"].append(dh.xyz[atom_idx])\n",
    "\n",
    "    # calculate size for Greens function generation\n",
    "    spin_box_shape = len(mag_ent[\"spin_box_indeces\"])\n",
    "\n",
    "    mag_ent[\"energies\"] = []  # we will store the second order energy derivations here\n",
    "\n",
    "    # These will be the perturbed potentials from eq. 100\n",
    "    mag_ent[\"Vu1\"] = []  # so they are independent in memory\n",
    "    mag_ent[\"Vu2\"] = []\n",
    "\n",
    "    mag_ent[\"Gii\"] = []  # Greens function\n",
    "    mag_ent[\"Gii_tmp\"] = []  # Greens function for parallelization\n",
    "    for i in ref_xcf_orientations:\n",
    "        # Rotations for every quantization axis\n",
    "        mag_ent[\"Vu1\"].append([])\n",
    "        mag_ent[\"Vu2\"].append([])\n",
    "        # Greens functions for every quantization axis\n",
    "        mag_ent[\"Gii\"].append(\n",
    "            np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n",
    "        )\n",
    "        mag_ent[\"Gii_tmp\"].append(\n",
    "            np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n",
    "        )\n",
    "\n",
    "# for every site we have to store 2x3 Greens function (and the associated _tmp-s)\n",
    "# in the 3 reference directions, because G_ij and G_ji are both needed\n",
    "for pair in pairs:\n",
    "    # calculate distance\n",
    "    xyz_ai = magnetic_entities[pair[\"ai\"]][\"xyz\"]\n",
    "    xyz_aj = magnetic_entities[pair[\"aj\"]][\"xyz\"]\n",
    "    xyz_aj = xyz_aj + pair[\"Ruc\"] @ simulation_parameters[\"cell\"]\n",
    "    pair[\"dist\"] = np.linalg.norm(xyz_ai - xyz_aj)\n",
    "\n",
    "    # calculate size for Greens function generation\n",
    "    spin_box_shape_i = len(magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"])\n",
    "    spin_box_shape_j = len(magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"])\n",
    "    pair[\"tags\"] = []\n",
    "    for mag_ent in [magnetic_entities[pair[\"ai\"]], magnetic_entities[pair[\"aj\"]]]:\n",
    "        tag = \"\"\n",
    "        # get atoms of magnetic entity\n",
    "        atoms_idx = mag_ent[\"atom\"]\n",
    "        orbitals = mag_ent[\"l\"]\n",
    "\n",
    "        # if magnetic entity contains one atoms\n",
    "        if isinstance(atoms_idx, int):\n",
    "            tag += f\"[{atoms_idx}]{dh.atoms[atoms_idx].tag}({orbitals})\"\n",
    "\n",
    "        # if magnetic entity contains more than one atoms\n",
    "        if isinstance(atoms_idx, list):\n",
    "            # iterate over atoms\n",
    "            atom_group = \"{\"\n",
    "            for atom_idx in atoms_idx:\n",
    "                atom_group += f\"[{atom_idx}]{dh.atoms[atom_idx].tag}({orbitals})--\"\n",
    "            # end {} of the atoms in the magnetic entity\n",
    "            tag += atom_group[:-2] + \"}\"\n",
    "        pair[\"tags\"].append(tag)\n",
    "    pair[\"energies\"] = []  # we will store the second order energy derivations here\n",
    "\n",
    "    pair[\"Gij\"] = []  # Greens function\n",
    "    pair[\"Gji\"] = []\n",
    "    pair[\"Gij_tmp\"] = []  # Greens function for parallelization\n",
    "    pair[\"Gji_tmp\"] = []\n",
    "    for i in ref_xcf_orientations:\n",
    "        # Greens functions for every quantization axis\n",
    "        pair[\"Gij\"].append(\n",
    "            np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype=\"complex128\")\n",
    "        )\n",
    "        pair[\"Gij_tmp\"].append(\n",
    "            np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype=\"complex128\")\n",
    "        )\n",
    "        pair[\"Gji\"].append(\n",
    "            np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n",
    "        )\n",
    "        pair[\"Gji_tmp\"].append(\n",
    "            np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n",
    "        )\n",
    "\n",
    "if rank == root_node:\n",
    "    times[\"site_and_pair_dictionaries_time\"] = timer()\n",
    "    print(\n",
    "        f\"Site and pair dictionaries created. Elapsed time: {times['site_and_pair_dictionaries_time']} s\"\n",
    "    )\n",
    "    print(\n",
    "        \"================================================================================================================================================================\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "k loop:   0%|          | 0/100 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k set created. Elapsed time: 3485.092548083 s\n",
      "================================================================================================================================================================\n"
     ]
    }
   ],
   "source": [
    "kset = make_kset(dirs=kdirs, NUMK=kset)  # generate k space sampling\n",
    "wkset = np.ones(len(kset)) / len(kset)  # generate weights for k points\n",
    "kpcs = np.array_split(kset, size)  # split the k points based on MPI size\n",
    "kpcs[root_node] = tqdm(kpcs[root_node], desc=\"k loop\")\n",
    "\n",
    "if rank == root_node:\n",
    "    times[\"k_set_time\"] = timer()\n",
    "    print(f\"k set created. Elapsed time: {times['k_set_time']} s\")\n",
    "    print(\n",
    "        \"================================================================================================================================================================\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rotations done perpendicular to quantization axis. Elapsed time: 3485.347085125 s\n",
      "================================================================================================================================================================\n"
     ]
    }
   ],
   "source": [
    "# this will contain the three hamiltonians in the reference directions needed to calculate the energy variations upon rotation\n",
    "hamiltonians = []\n",
    "\n",
    "# iterate over the reference directions (quantization axes)\n",
    "for i, orient in enumerate(ref_xcf_orientations):\n",
    "    # obtain rotated exchange field\n",
    "    R = RotMa2b(scf_xcf_orientation, orient[\"o\"])\n",
    "    rot_XCF = np.einsum(\"ij,jklm->iklm\", R, XCF)\n",
    "    rot_H_XCF = sum(\n",
    "        [np.kron(rot_XCF[i], tau) for i, tau in enumerate([tau_x, tau_y, tau_z])]\n",
    "    )\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 = (\n",
    "        hTRS + rot_H_XCF\n",
    "    )  # equation 76 #######################################################################################\n",
    "\n",
    "    hamiltonians.append(\n",
    "        dict(orient=orient[\"o\"], H=rot_H)\n",
    "    )  # store orientation and rotated Hamiltonian\n",
    "\n",
    "    # these are the rotations (for now) perpendicular to the quantization axis\n",
    "    for u in orient[\"vw\"]:\n",
    "        Tu = np.kron(np.eye(NO, dtype=int), tau_u(u))  # section 2.H\n",
    "\n",
    "        Vu1 = 1j / 2 * commutator(rot_H_XCF_uc, Tu)  # equation 100\n",
    "        Vu2 = 1 / 8 * commutator(commutator(Tu, rot_H_XCF_uc), Tu)  # equation 100\n",
    "\n",
    "        for mag_ent in magnetic_entities:\n",
    "            idx = mag_ent[\"spin_box_indeces\"]\n",
    "            # fill up the perturbed potentials (for now) based on the on-site projections\n",
    "            mag_ent[\"Vu1\"][i].append(Vu1[:, idx][idx, :])\n",
    "            mag_ent[\"Vu2\"][i].append(Vu2[:, idx][idx, :])\n",
    "\n",
    "if rank == root_node:\n",
    "    times[\"reference_rotations_time\"] = timer()\n",
    "    print(\n",
    "        f\"Rotations done perpendicular to quantization axis. Elapsed time: {times['reference_rotations_time']} s\"\n",
    "    )\n",
    "    print(\n",
    "        \"================================================================================================================================================================\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x16b694a30>"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAGkCAYAAABtmxHBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACytElEQVR4nO29eZwU1dU+/swwMCAKiMqmIKgg4gIEBFETMfILLokgGqNBg0viG/ctRn0TNJooLvkad42JUZO4JCZiXk2iURTUsIgg7gIqKgoDKgKyDbP0749bt6r63Od23WlmGHo4z+fjR6a7uqruqaru85zznHPKcrlcDgqFQqFQbOEob+4TUCgUCoUiBPqDpVAoFIqSgP5gKRQKhaIkoD9YCoVCoSgJ6A+WQqFQKEoC+oOlUCgUipKA/mApFAqFoiSgP1gKhUKhKAnoD5ZCoVAoSgL6g6VQKBSKkkDJ/mDdcccd6N27N9q2bYvhw4fj5Zdfbu5TAgBMmjQJ+++/P7bbbjt06dIFY8eOxfz58/O22bBhA84++2zssMMO2HbbbXHsscdi2bJlzXTG+bjuuutQVlaGCy64IH5tSzvfTz/9FCeddBJ22GEHtGvXDvvuuy9eeeWV+P1cLocrrrgC3bt3R7t27TBq1CgsXLiw2c63rq4OEydORJ8+fdCuXTvsvvvu+OUvf4l0V7TmPOcXXngB3/nOd9CjRw+UlZXh8ccfz3s/5NxWrFiB8ePHo0OHDujUqRNOP/10rFmzZrOfb01NDS699FLsu+++aN++PXr06IEf/OAHWLJkSbOdb9Y5S/z4xz9GWVkZbr755mY755Dzfeedd3D00UejY8eOaN++Pfbff398/PHH8ftN8b1Rkj9Yf/nLX3DRRRfhyiuvxNy5czFw4ECMHj0ay5cvb+5Tw7Rp03D22Wdj5syZeOaZZ1BTU4NvfetbWLt2bbzNhRdeiCeeeAKPPvoopk2bhiVLlmDcuHHNeNYGs2fPxm9/+1vst99+ea9vSef75Zdf4qCDDkLr1q3x73//G2+//Tb+3//7f9h+++3jbW644QbceuutuPvuuzFr1iy0b98eo0ePxoYNG5rlnK+//nrcdddduP322/HOO+/g+uuvxw033IDbbrttizjntWvXYuDAgbjjjjvo+yHnNn78eLz11lt45pln8OSTT+KFF17AGWecsdnPd926dZg7dy4mTpyIuXPn4rHHHsP8+fNx9NFH5223Oc8365zTmDx5MmbOnIkePXo4720pNgaA999/HwcffDD69++PqVOn4vXXX8fEiRPRtm3beJsm+d7IlSCGDRuWO/vss+O/6+rqcj169MhNmjSpGc+KY/ny5TkAuWnTpuVyuVxu5cqVudatW+ceffTReJt33nknByA3Y8aM5jrN3FdffZXr27dv7plnnskdcsghufPPPz+Xy21553vppZfmDj74YO/79fX1uW7duuVuvPHG+LWVK1fmKisrcw8//PDmOEUHRx11VO60007Le23cuHG58ePH53K5LeucAeQmT54c/x1ybm+//XYOQG727NnxNv/+979zZWVluU8//XSzni/Dyy+/nAOQ++ijj3K5XPOeby7nP+dPPvkkt/POO+fefPPN3K677pr7zW9+E7+3pdn4e9/7Xu6kk07yfqapvjdKjmFt3LgRc+bMwahRo+LXysvLMWrUKMyYMaMZz4xj1apVAIDOnTsDAObMmYOampq88+/fvz969erVrOd/9tln46ijjso7L2DLO9//+7//w9ChQ/Hd734XXbp0weDBg/G73/0ufn/RokWoqqrKO9+OHTti+PDhzWbfAw88EFOmTMGCBQsAAK+99hpeeuklHHHEEVvsOVuEnNuMGTPQqVMnDB06NN5m1KhRKC8vx6xZszb7OUusWrUKZWVl6NSpE4At83zr6+tx8skn45JLLsHee+/tvL8lnXN9fT3++c9/ol+/fhg9ejS6dOmC4cOH54UNm+p7o+R+sD7//HPU1dWha9euea937doVVVVVzXRWHPX19bjgggtw0EEHYZ999gEAVFVVoU2bNvHDY9Gc5//II49g7ty5mDRpkvPelna+H3zwAe666y707dsXTz/9NM4880ycd955eOCBB+Lztee3JZwvAFx22WU44YQT0L9/f7Ru3RqDBw/GBRdcgPHjxwPYMs/ZIuTcqqqq0KVLl7z3Kyoq0Llz52Y//w0bNuDSSy/FiSeeiA4dOgDYMs/3+uuvR0VFBc477zz6/pZ0zsuXL8eaNWtw3XXX4fDDD8d//vMfHHPMMRg3bhymTZsWn29TfG9UbMqJKwrj7LPPxptvvomXXnqpuU/Fi8WLF+P888/HM888kxd/3lJRX1+PoUOH4tprrwUADB48GG+++SbuvvtuTJgwoZnPjuOvf/0rHnzwQTz00EPYe++9MW/ePFxwwQXo0aPHFnvOLQE1NTU4/vjjkcvlcNdddzX36XgxZ84c3HLLLZg7dy7Kysqa+3QyUV9fDwAYM2YMLrzwQgDAoEGDMH36dNx999045JBDmuzYJcewdtxxR7Rq1cpRmyxbtgzdunVrprNycc455+DJJ5/E888/j1122SV+vVu3bti4cSNWrlyZt31znf+cOXOwfPlyfO1rX0NFRQUqKiowbdo03HrrraioqEDXrl23qPPt3r07BgwYkPfaXnvtFauT7DltSffHJZdcErOsfffdFyeffDIuvPDCmNFuiedsEXJu3bp1cwRPtbW1WLFiRbOdv/2x+uijj/DMM8/E7ArY8s73xRdfxPLly9GrV6/4Gfzoo49w8cUXo3fv3lvcOe+4446oqKjIfA6b4nuj5H6w2rRpgyFDhmDKlCnxa/X19ZgyZQpGjBjRjGdmkMvlcM4552Dy5Ml47rnn0KdPn7z3hwwZgtatW+ed//z58/Hxxx83y/kfdthheOONNzBv3rz4v6FDh2L8+PHxv7ek8z3ooIOcMoEFCxZg1113BQD06dMH3bp1yzvf1atXY9asWc12f6xbtw7l5fmPWqtWrWJPdUs8Z4uQcxsxYgRWrlyJOXPmxNs899xzqK+vx/Dhwzf7Odsfq4ULF+LZZ5/FDjvskPf+lna+J598Ml5//fW8Z7BHjx645JJL8PTTT29x59ymTRvsv//+BZ/DJvueK1qu0Yx45JFHcpWVlbn7778/9/bbb+fOOOOMXKdOnXJVVVXNfWq5M888M9exY8fc1KlTc0uXLo3/W7duXbzNj3/841yvXr1yzz33XO6VV17JjRgxIjdixIhmPOt8pFWCudyWdb4vv/xyrqKiInfNNdfkFi5cmHvwwQdz22yzTe7Pf/5zvM11112X69SpU+4f//hH7vXXX8+NGTMm16dPn9z69eub5ZwnTJiQ23nnnXNPPvlkbtGiRbnHHnsst+OOO+Z++tOfbhHn/NVXX+VeffXV3KuvvpoDkLvppptyr776aqyqCzm3ww8/PDd48ODcrFmzci+99FKub9++uRNPPHGzn+/GjRtzRx99dG6XXXbJzZs3L+8ZrK6ubpbzzTpnBqkS3NznnHW+jz32WK5169a5e+65J7dw4cLcbbfdlmvVqlXuxRdfjPfRFN8bJfmDlcvlcrfddluuV69euTZt2uSGDRuWmzlzZnOfUi6XMxJQ9t99990Xb7N+/frcWWedldt+++1z22yzTe6YY47JLV26tPlOWkD+YG1p5/vEE0/k9tlnn1xlZWWuf//+uXvuuSfv/fr6+tzEiRNzXbt2zVVWVuYOO+yw3Pz585vpbHO51atX584///xcr169cm3bts3ttttuuZ/97Gd5X6DNec7PP/88vWcnTJgQfG5ffPFF7sQTT8xtu+22uQ4dOuROPfXU3FdffbXZz3fRokXeZ/D5559vlvPNOmcG9oO1pdjY4t57783tscceubZt2+YGDhyYe/zxx/P20RTfG2W5XKrcXqFQKBSKLRQll8NSKBQKxdYJ/cFSKBQKRUlAf7AUCoVCURLQHyyFQqFQlAT0B0uhUCgUJQH9wVIoFApFSaBkf7Cqq6vxi1/8AtXV1c19KsEotXPW821alNr5AqV3znq+TY/
      "text/plain": [
       "<Figure size 480x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAGkCAYAAABtmxHBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACxeUlEQVR4nO19eZgV1bX96gEaBAFB6QYFQQUBByCgiPIiKr/gkNgIRlFicHiaOCPGqEnQqFE0vmeMxjHPqEmcYiIaTWKiIDjQIoI4gQwyCEqDiszQY/3+OHWqTu2zz63Tl26ae9nr+/yk69Zwzq6qe8/ae+29C4IgCCAQCAQCwS6OwuYegEAgEAgEPpAfLIFAIBDkBOQHSyAQCAQ5AfnBEggEAkFOQH6wBAKBQJATkB8sgUAgEOQE5AdLIBAIBDkB+cESCAQCQU5AfrAEAoFAkBOQHyyBQCAQ5ARy9gfrvvvuQ48ePdCqVSsMGTIE77zzTnMPCQAwefJkHHHEEdhzzz3RuXNnjBo1CgsXLkzss337dlx66aXo1KkT2rZtizFjxmDNmjXNNOIkbr/9dhQUFGDChAnRtl1tvJ9//jl+8IMfoFOnTmjdujUOO+wwvPvuu9HnQRDghhtuQJcuXdC6dWuMGDECixcvbrbx1tXVYdKkSejZsydat26NAw88ELfccgvMqmjNOebXX38d3/ve99C1a1cUFBTg+eefT3zuM7Z169Zh3LhxaNeuHTp06IALLrgAmzdv3unjrampwbXXXovDDjsMbdq0QdeuXfHDH/4QX3zxRbONN23MFD/+8Y9RUFCAu+++u9nG7DPeBQsW4NRTT0X79u3Rpk0bHHHEEfjss8+iz5vieyMnf7CeeeYZTJw4ETfeeCPmzp2L/v37Y+TIkVi7dm1zDw0zZszApZdeirfffhuvvPIKampq8J3vfAdbtmyJ9rnqqqvw4osv4tlnn8WMGTPwxRdfYPTo0c04aoXZs2fjoYcewuGHH57YviuN95tvvsExxxyDFi1a4F//+hfmz5+P//3f/8Vee+0V7fPrX/8a99xzDx588EHMmjULbdq0wciRI7F9+/ZmGfMdd9yBBx54AL/73e+wYMEC3HHHHfj1r3+Ne++9d5cY85YtW9C/f3/cd9997Oc+Yxs3bhw+/vhjvPLKK3jppZfw+uuv46KLLtrp4926dSvmzp2LSZMmYe7cuXjuueewcOFCnHrqqYn9duZ408ZsYsqUKXj77bfRtWtX67NdxcYA8Omnn2LYsGHo06cPpk+fjg8++ACTJk1Cq1aton2a5HsjyEEceeSRwaWXXhr9XVdXF3Tt2jWYPHlyM46Kx9q1awMAwYwZM4IgCIL169cHLVq0CJ599tlonwULFgQAgoqKiuYaZrBp06agV69ewSuvvBIce+yxwZVXXhkEwa433muvvTYYNmyY8/P6+vqgrKwsuPPOO6Nt69evD0pKSoKnnnpqZwzRwimnnBKcf/75iW2jR48Oxo0bFwTBrjVmAMGUKVOiv33GNn/+/ABAMHv27Giff/3rX0FBQUHw+eef79TxcnjnnXcCAMGKFSuCIGje8QaBe8yrVq0K9t133+Cjjz4K9t9//+A3v/lN9NmuZuMzzzwz+MEPfuA8pqm+N3KOYVVXV2POnDkYMWJEtK2wsBAjRoxARUVFM46Mx4YNGwAAHTt2BADMmTMHNTU1ifH36dMH3bt3b9bxX3rppTjllFMS4wJ2vfH+/e9/x+DBg/H9738fnTt3xsCBA/H73/8++nzZsmWorKxMjLd9+/YYMmRIs9n36KOPxtSpU7Fo0SIAwPvvv48333wTJ5100i47Zg2fsVVUVKBDhw4YPHhwtM+IESNQWFiIWbNm7fQxU2zYsAEFBQXo0KEDgF1zvPX19TjnnHNwzTXX4JBDDrE+35XGXF9fj3/84x/o3bs3Ro4cic6dO2PIkCEJt2FTfW/k3A/WV199hbq6OpSWlia2l5aWorKysplGxaO+vh4TJkzAMcccg0MPPRQAUFlZiZYtW0Yvj0Zzjv/pp5/G3LlzMXnyZOuzXW28S5cuxQMPPIBevXrh3//+Ny6++GJcccUVePzxx6Px6vHtCuMFgOuuuw5jx45Fnz590KJFCwwcOBATJkzAuHHjAOyaY9bwGVtlZSU6d+6c+Ly4uBgdO3Zs9vFv374d1157Lc466yy0a9cOwK453jvuuAPFxcW44oor2M93pTGvXbsWmzdvxu23344TTzwR//nPf3Daaadh9OjRmDFjRjTepvjeKN6RgQsy49JLL8VHH32EN998s7mH4sTKlStx5ZVX4pVXXkn4n3dV1NfXY/DgwbjtttsAAAMHDsRHH32EBx98EOPHj2/m0fH4y1/+gieeeAJPPvkkDjnkEMybNw8TJkxA165dd9kx5wNqampwxhlnIAgCPPDAA809HCfmzJmD3/72t5g7dy4KCgqaezipqK+vBwCUl5fjqquuAgAMGDAAM2fOxIMPPohjjz22ya6dcwxr7733RlFRkaU2WbNmDcrKypppVDYuu+wyvPTSS3jttdew3377RdvLyspQXV2N9evXJ/ZvrvHPmTMHa9euxbe+9S0UFxejuLgYM2bMwD333IPi4mKUlpbuUuPt0qUL+vXrl9jWt2/fSJ2kx7QrPR/XXHNNxLIOO+wwnHPOObjqqqsiRrsrjlnDZ2xlZWWW4Km2thbr1q1rtvHrH6sVK1bglVdeidgVsOuN94033sDatWvRvXv36B1csWIFrr76avTo0WOXG/Pee++N4uLi1PewKb43cu4Hq2XLlhg0aBCmTp0abauvr8fUqVMxdOjQZhyZQhAEuOyyyzBlyhRMmzYNPXv2THw+aNAgtGjRIjH+hQsX4rPPPmuW8Z9wwgn48MMPMW/evOi/wYMHY9y4cdG/d6XxHnPMMVaawKJFi7D//vsDAHr27ImysrLEeDdu3IhZs2Y12/OxdetWFBYmX7WioqJopborjlnDZ2xDhw7F+vXrMWfOnGifadOmob6+HkOGDNnpY9Y/VosXL8arr76KTp06JT7f1cZ7zjnn4IMPPki8g127dsU111yDf//737vcmFu2bIkjjjgi43vYZN9zWcs1mhFPP/10UFJSEjz22GPB/Pnzg4suuijo0KFDUFlZ2dxDCy6++OKgffv2wfTp04PVq1dH/23dujXa58c//nHQvXv3YNq0acG7774bDB06NBg6dGgzjjoJUyUYBLvWeN95552guLg4uPXWW4PFixcHTzzxRLDHHnsEf/7zn6N9br/99qBDhw7BCy+8EHzwwQdBeXl50LNnz2Dbtm3NMubx48cH++67b/DSSy8Fy5YtC5577rlg7733Dn7605/uEmPetGlT8N577wXvvfdeACC46667gvfeey9S1fmM7cQTTwwGDhwYzJo1K3jzzTeDXr16BWedddZOH291dXVw6qmnBvvtt18wb968xDtYVVXVLONNGzMHqhLc2WNOG+9zzz0XtGjRInj44YeDxYsXB/fee29QVFQUvPHGG9E5muJ7Iyd/sIIgCO69996ge/fuQcuWLYMjjzwyePvtt5t7SEEQKAko99+jjz4a7bNt27bgkksuCfbaa69gjz32CE477bRg9erVzTdoAvqDtauN98UXXwwOPfTQoKSkJOjTp0/w8MMPJz6vr68PJk2aFJSWlgYlJSXBCSecECxcuLCZRhsEGzduDK688sqge/fuQatWrYIDDjgg+PnPf574Am3OMb/22mvsMzt+/HjvsX399dfBWWedFbRt2zZo165dcN555wWbNm3a6eNdtmyZ8x187bXXmmW8aWPmwP1g7So21njkkUeCgw46KGjVqlXQv3//4Pnnn0+coym+NwqCwEi3FwgEAoFgF0XOxbAEAoFAsHtCfrAEAoFAkBOQHyyBQCAQ5ATkB0sgEAgEOQH5wRIIBAJBTkB+sAQCgUCQE8jZH6yqqir88pe/RFVVVXMPxRu5NmYZb9Mi18YL5N6YZbxNj5055pzNw9q4cSPat2+PDRs2JOqE7crItTHLeJsWuTZeIPfGLONteuzMMTcrw9pV29wLBAKBYNd
      "text/plain": [
       "<Figure size 480x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "HK, SK = hsk(hamiltonians[2][\"H\"], ss, dh.sc_off, np.array([0.3, 0, 0]))\n",
    "\n",
    "myhk = dh.Hk(k=np.array([0.3, 0, 0]))\n",
    "\n",
    "abs(HK - myhk.toarray()).max()\n",
    "\n",
    "plt.matshow(abs(HK))\n",
    "plt.matshow(abs(myhk.toarray()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting matrix inversions\n",
      "Total number of k points: 100\n",
      "Number of energy samples per k point: 600\n",
      "Total number of directions: 3\n",
      "Total number of matrix inversions: 180000\n",
      "The shape of the Hamiltonian and the Greens function is 84x84=7056\n",
      "Memory taken by a single Hamiltonian is: 0.015625 KB\n",
      "Expected memory usage per matrix inversion: 0.5 KB\n",
      "Expected memory usage per k point for parallel inversion: 900.0 KB\n",
      "Expected memory usage on root node: 87.890625 MB\n",
      "================================================================================================================================================================\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "k loop: 100%|██████████| 100/100 [1:39:50<00:00, 59.90s/it]  "
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculated Greens functions. Elapsed time: 3878.870911625 s\n",
      "================================================================================================================================================================\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "if rank == root_node:\n",
    "    print(\"Starting matrix inversions\")\n",
    "    print(f\"Total number of k points: {kset.shape[0]}\")\n",
    "    print(f\"Number of energy samples per k point: {eset}\")\n",
    "    print(f\"Total number of directions: {len(hamiltonians)}\")\n",
    "    print(\n",
    "        f\"Total number of matrix inversions: {kset.shape[0] * len(hamiltonians) * eset}\"\n",
    "    )\n",
    "    print(f\"The shape of the Hamiltonian and the Greens function is {NO}x{NO}={NO*NO}\")\n",
    "    # https://stackoverflow.com/questions/70746660/how-to-predict-memory-requirement-for-np-linalg-inv\n",
    "    # memory is O(64 n**2) for complex matrices\n",
    "    memory_size = getsizeof(hamiltonians[0][\"H\"].base) / 1024\n",
    "    print(\n",
    "        f\"Memory taken by a single Hamiltonian is: {getsizeof(hamiltonians[0]['H'].base) / 1024} KB\"\n",
    "    )\n",
    "    print(f\"Expected memory usage per matrix inversion: {memory_size * 32} KB\")\n",
    "    print(\n",
    "        f\"Expected memory usage per k point for parallel inversion: {memory_size * len(hamiltonians) * eset * 32} KB\"\n",
    "    )\n",
    "    print(\n",
    "        f\"Expected memory usage on root node: {len(np.array_split(kset, size)[0]) * memory_size * len(hamiltonians) * eset * 32 / 1024} MB\"\n",
    "    )\n",
    "    print(\n",
    "        \"================================================================================================================================================================\"\n",
    "    )\n",
    "\n",
    "comm.Barrier()\n",
    "# ----------------------------------------------------------------------\n",
    "\n",
    "# make energy contour\n",
    "# we are working in eV now  !\n",
    "# and sisl shifts E_F to 0 !\n",
    "cont = make_contour(emin=ebot, enum=eset, p=esetp)\n",
    "eran = cont.ze\n",
    "\n",
    "# ----------------------------------------------------------------------\n",
    "# sampling the integrand on the contour and the BZ\n",
    "for k in kpcs[rank]:\n",
    "    wk = wkset[rank]  # weight of k point in BZ integral\n",
    "    # iterate over reference directions\n",
    "    for i, hamiltonian_orientation in enumerate(hamiltonians):\n",
    "        # calculate Greens function\n",
    "        H = hamiltonian_orientation[\"H\"]\n",
    "        HK, SK = hsk(H, ss, dh.sc_off, k)\n",
    "        # HK = HK - Ef * SK\n",
    "        # Gk = inv(SK * eran.reshape(eset, 1, 1) - HK)\n",
    "\n",
    "        # solve Greens function sequentially for the energies, because of memory bound\n",
    "        Gk = np.zeros(shape=(eset, HK.shape[0], HK.shape[1]), dtype=\"complex128\")\n",
    "        for j in range(eset):\n",
    "            Gk[j] = inv(SK * eran[j] - HK)\n",
    "\n",
    "        # store the Greens function slice of the magnetic entities (for now) based on the on-site projections\n",
    "        for mag_ent in magnetic_entities:\n",
    "            mag_ent[\"Gii_tmp\"][i] += (\n",
    "                Gk[:, mag_ent[\"spin_box_indeces\"], :][:, :, mag_ent[\"spin_box_indeces\"]]\n",
    "                * wk\n",
    "            )\n",
    "\n",
    "        for pair in pairs:\n",
    "            # add phase shift based on the cell difference\n",
    "            phase = np.exp(1j * 2 * np.pi * k @ pair[\"Ruc\"].T)\n",
    "\n",
    "            # get the pair orbital sizes from the magnetic entities\n",
    "            ai = magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"]\n",
    "            aj = magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]\n",
    "\n",
    "            # store the Greens function slice of the magnetic entities (for now) based on the on-site projections\n",
    "            pair[\"Gij_tmp\"][i] += Gk[:, ai][..., aj] * phase * wk\n",
    "            pair[\"Gji_tmp\"][i] += Gk[:, aj][..., ai] / phase * wk\n",
    "\n",
    "# summ reduce partial results of mpi nodes\n",
    "for i in range(len(hamiltonians)):\n",
    "    for mag_ent in magnetic_entities:\n",
    "        comm.Reduce(mag_ent[\"Gii_tmp\"][i], mag_ent[\"Gii\"][i], root=root_node)\n",
    "\n",
    "    for pair in pairs:\n",
    "        comm.Reduce(pair[\"Gij_tmp\"][i], pair[\"Gij\"][i], root=root_node)\n",
    "        comm.Reduce(pair[\"Gji_tmp\"][i], pair[\"Gji\"][i], root=root_node)\n",
    "\n",
    "if rank == root_node:\n",
    "    times[\"green_function_inversion_time\"] = timer()\n",
    "    print(\n",
    "        f\"Calculated Greens functions. Elapsed time: {times['green_function_inversion_time']} s\"\n",
    "    )\n",
    "    print(\n",
    "        \"================================================================================================================================================================\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Magnetic entities integrated.\n",
      "Pairs integrated.\n",
      "Magnetic parameters calculated.\n",
      "##################################################################### GROGU OUTPUT #############################################################################\n",
      "================================================================================================================================================================\n",
      "Input file: \n",
      "/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\n",
      "Output file: \n",
      "./Fe3GeTe2.pickle\n",
      "Number of nodes in the parallel cluster:  1\n",
      "================================================================================================================================================================\n",
      "Cell [Ang]: \n",
      "[[ 3.79100000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [-1.89550000e+00  3.28310231e+00  0.00000000e+00]\n",
      " [ 1.25954923e-15  2.18160327e-15  2.05700000e+01]]\n",
      "================================================================================================================================================================\n",
      "DFT axis: \n",
      "[0 0 1]\n",
      "Quantization axis and perpendicular rotation directions:\n",
      "[1 0 0]  --»  [array([0, 1, 0]), array([0, 0, 1])]\n",
      "[0 1 0]  --»  [array([1, 0, 0]), array([0, 0, 1])]\n",
      "[0 0 1]  --»  [array([1, 0, 0]), array([0, 1, 0])]\n",
      "================================================================================================================================================================\n",
      "Parameters for the contour integral:\n",
      "Number of k points:  10\n",
      "k point directions:  xy\n",
      "Ebot:  -13\n",
      "Eset:  600\n",
      "Esetp:  10000\n",
      "================================================================================================================================================================\n",
      "Atomic information: \n",
      "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
      "[atom index]Element(orbitals)        x [Ang]        y [Ang]        z [Ang]        Sx        Sy        Sz        Q        Lx        Ly        Lz        Jx        Jy        Jz\n",
      "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
      "[3]Fe(2)        -7.339158738013707e-06        4.149278510690423e-06        11.657585837928032\n",
      "\n",
      "[4]Fe(all)        -7.326987662162937e-06        4.158274523275774e-06        8.912422537596708\n",
      "\n",
      "[5]Fe(2)        1.8954667088117545        1.0943913231921656        10.285002698393109\n",
      "\n",
      "================================================================================================================================================================\n",
      "Exchange [meV]\n",
      "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
      "Magnetic entity1        Magnetic entity2       [i  j  k]       d [Ang]\n",
      "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
      "[3]Fe(2)        [4]Fe(all)        [0 0 0]        d [Ang] 2.745163300331324\n",
      "Isotropic:  -38.08908138560059\n",
      "DMI:  [-3.50315141e-01  2.63885818e-04 -4.29362132e-06]\n",
      "Symmetric-anisotropy:  [ 1.17647885e+00 -1.96468230e-05 -7.74651208e-07 -1.96468230e-05\n",
      "  1.26142216e+00  3.57154955e-04 -7.74651208e-07  3.57154955e-04\n",
      " -2.43790101e+00]\n",
      "J:  [-3.69126025e+01 -1.96468230e-05 -7.74651208e-07 -1.96468230e-05\n",
      " -3.68276592e+01  3.57154955e-04 -7.74651208e-07  3.57154955e-04\n",
      " -4.05269824e+01]\n",
      "Energies for debugging: \n",
      "array([[-4.04406378e-02, -3.50672296e-04,  3.49957987e-04,\n",
      "        -4.01805196e-02],\n",
      "       [-4.06133270e-02, -2.63111167e-07,  2.64660469e-07,\n",
      "        -4.03504319e-02],\n",
      "       [-3.34747989e-02,  1.53532017e-08,  2.39404444e-08,\n",
      "        -3.34747732e-02]])\n",
      "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
      "array([-0.04035043, -0.0334748 , -0.04044064])\n",
      "Test J_xx = E(y,z) = E(z,y)\n",
      "-0.04035043189424629 -0.03347477318614506\n",
      "\n",
      "[3]Fe(2)        [5]Fe(2)        [0 0 0]        d [Ang] 2.5835033632437767\n",
      "Isotropic:  -65.6147625726105\n",
      "DMI:  [ 3.55875702e+00 -6.14766359e+00  2.13990280e-05]\n",
      "Symmetric-anisotropy:  [-0.08599962  0.03698877 -0.10152953  0.03698877 -0.29087488 -0.04979049\n",
      " -0.10152953 -0.04979049  0.3768745 ]\n",
      "J:  [-6.57007622e+01  3.69887710e-02 -1.01529530e-01  3.69887710e-02\n",
      " -6.59056375e+01 -4.97904945e-02 -1.01529530e-01 -4.97904945e-02\n",
      " -6.52378881e+01]\n",
      "Energies for debugging: \n",
      "array([[-6.53648847e-02,  3.60854752e-03, -3.50896653e-03,\n",
      "        -6.58426415e-02],\n",
      "       [-6.51108914e-02,  6.24919312e-03, -6.04613406e-03,\n",
      "        -6.54759423e-02],\n",
      "       [-6.59686334e-02, -3.69673720e-05, -3.70101700e-05,\n",
      "        -6.59255821e-02]])\n",
      "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
      "array([-0.06547594, -0.06596863, -0.06536488])\n",
      "Test J_xx = E(y,z) = E(z,y)\n",
      "-0.06547594226467608 -0.06592558212190268\n",
      "\n",
      "[4]Fe(all)        [5]Fe(2)        [0 0 0]        d [Ang] 2.583501767937866\n",
      "Isotropic:  -63.461479457470716\n",
      "DMI:  [-3.46188832e+00  5.96672549e+00  2.26187324e-05]\n",
      "Symmetric-anisotropy:  [-0.0829036   0.03178786  0.09341113  0.03178786 -0.2854846   0.04272782\n",
      "  0.09341113  0.04272782  0.3683882 ]\n",
      "J:  [-6.35443831e+01  3.17878581e-02  9.34111270e-02  3.17878581e-02\n",
      " -6.37469641e+01  4.27278222e-02  9.34111270e-02  4.27278222e-02\n",
      " -6.30930913e+01]\n",
      "Energies for debugging: \n",
      "array([[-6.32236137e-02, -3.50461614e-03,  3.41916050e-03,\n",
      "        -6.36877197e-02],\n",
      "       [-6.29625688e-02, -6.06013662e-03,  5.87331436e-03,\n",
      "        -6.33196114e-02],\n",
      "       [-6.38062084e-02, -3.17652393e-05, -3.18104768e-05,\n",
      "        -6.37691547e-02]])\n",
      "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
      "array([-0.06331961, -0.06380621, -0.06322361])\n",
      "Test J_xx = E(y,z) = E(z,y)\n",
      "-0.06331961140706378 -0.06376915470527415\n",
      "\n",
      "[3]Fe(2)        [5]Fe(2)        [-1 -1  0]        d [Ang] 2.5834973202859075\n",
      "Isotropic:  -65.62000802135631\n",
      "DMI:  [-7.07503443e+00  3.04878078e-04  3.59246447e-06]\n",
      "Symmetric-anisotropy:  [-3.95641341e-01 -1.37462093e-04  1.36999229e-04 -1.37462093e-04\n",
      "  1.87904813e-02  1.22453077e-01  1.36999229e-04  1.22453077e-01\n",
      "  3.76850860e-01]\n",
      "J:  [-6.60156494e+01 -1.37462093e-04  1.36999229e-04 -1.37462093e-04\n",
      " -6.56012175e+01  1.22453077e-01  1.36999229e-04  1.22453077e-01\n",
      " -6.52431572e+01]\n",
      "Energies for debugging: \n",
      "array([[-6.49844951e-02, -7.19748751e-03,  6.95258136e-03,\n",
      "        -6.52968607e-02],\n",
      "       [-6.55018192e-02, -4.41877308e-07,  1.67878849e-07,\n",
      "        -6.60401128e-02],\n",
      "       [-6.59055744e-02,  1.41054558e-07,  1.33869629e-07,\n",
      "        -6.59911860e-02]])\n",
      "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
      "array([-0.06604011, -0.06590557, -0.0649845 ])\n",
      "Test J_xx = E(y,z) = E(z,y)\n",
      "-0.06604011276766036 -0.06599118595746928\n",
      "\n",
      "[4]Fe(all)        [5]Fe(2)        [-1 -1  0]        d [Ang] 2.583495745338251\n",
      "Isotropic:  -63.4637006802454\n",
      "DMI:  [ 6.84888899e+00 -7.88878725e-04  1.52778682e-05]\n",
      "Symmetric-anisotropy:  [-3.84273968e-01 -1.41648052e-04  1.12499451e-04 -1.41648052e-04\n",
      "  1.63124769e-02 -1.19287820e-01  1.12499451e-04 -1.19287820e-01\n",
      "  3.67961491e-01]\n",
      "J:  [-6.38479746e+01 -1.41648052e-04  1.12499451e-04 -1.41648052e-04\n",
      " -6.34473882e+01 -1.19287820e-01  1.12499451e-04 -1.19287820e-01\n",
      " -6.30957392e+01]\n",
      "Energies for debugging: \n",
      "array([[-6.28360612e-02,  6.96817681e-03, -6.72960117e-03,\n",
      "        -6.31426268e-02],\n",
      "       [-6.33554172e-02,  6.76379274e-07, -9.01378176e-07,\n",
      "        -6.38701956e-02],\n",
      "       [-6.37521496e-02,  1.56925920e-07,  1.26370184e-07,\n",
      "        -6.38257537e-02]])\n",
      "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
      "array([-0.0638702 , -0.06375215, -0.06283606])\n",
      "Test J_xx = E(y,z) = E(z,y)\n",
      "-0.06387019555954467 -0.06382575373745485\n",
      "\n",
      "[3]Fe(2)        [5]Fe(2)        [-1  0  0]        d [Ang] 2.583541444641373\n",
      "Isotropic:  -65.60128795737425\n",
      "DMI:  [ 3.57594200e+00  6.14751321e+00 -2.99860089e-05]\n",
      "Symmetric-anisotropy:  [-0.08550091 -0.03685709  0.10144213 -0.03685709 -0.29158846 -0.04351401\n",
      "  0.10144213 -0.04351401  0.37708937]\n",
      "J:  [-6.56867889e+01 -3.68570877e-02  1.01442127e-01 -3.68570877e-02\n",
      " -6.58928764e+01 -4.35140091e-02  1.01442127e-01 -4.35140091e-02\n",
      " -6.52241986e+01]\n",
      "Energies for debugging: \n",
      "array([[-6.53521011e-02,  3.61945601e-03, -3.53242799e-03,\n",
      "        -6.58308241e-02],\n",
      "       [-6.50962960e-02, -6.24895534e-03,  6.04607109e-03,\n",
      "        -6.54612308e-02],\n",
      "       [-6.59549288e-02,  3.68271017e-05,  3.68870737e-05,\n",
      "        -6.59123469e-02]])\n",
      "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
      "array([-0.06546123, -0.06595493, -0.0653521 ])\n",
      "Test J_xx = E(y,z) = E(z,y)\n",
      "-0.06546123082429328 -0.06591234690586369\n",
      "\n",
      "[4]Fe(all)        [5]Fe(2)        [-1  0  0]        d [Ang] 2.5835398672184064\n",
      "Isotropic:  -63.445681369296956\n",
      "DMI:  [-3.44476779e+00 -5.96603346e+00 -2.37584314e-05]\n",
      "Symmetric-anisotropy:  [-0.08466531 -0.031648   -0.0935334  -0.031648   -0.28385999  0.04835705\n",
      " -0.0935334   0.04835705  0.3685253 ]\n",
      "J:  [-6.35303467e+01 -3.16479979e-02 -9.35333990e-02 -3.16479979e-02\n",
      " -6.37295414e+01  4.83570464e-02 -9.35333990e-02  4.83570464e-02\n",
      " -6.30771561e+01]\n",
      "Energies for debugging: \n",
      "array([[-6.32065968e-02, -3.49312484e-03,  3.39641074e-03,\n",
      "        -6.36666376e-02],\n",
      "       [-6.29477154e-02,  6.05956686e-03, -5.87250006e-03,\n",
      "        -6.33048155e-02],\n",
      "       [-6.37924451e-02,  3.16242395e-05,  3.16717563e-05,\n",
      "        -6.37558779e-02]])\n",
      "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
      "array([-0.06330482, -0.06379245, -0.0632066 ])\n",
      "Test J_xx = E(y,z) = E(z,y)\n",
      "-0.06330481549381357 -0.06375587786178244\n",
      "\n",
      "================================================================================================================================================================\n",
      "Runtime information: \n",
      "Total runtime: 395.5379132080002 s\n",
      "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
      "Initial setup: 0.14075858300020627 s\n",
      "Hamiltonian conversion and XC field extraction: 0.811 s\n",
      "Pair and site datastructure creatrions: 0.094 s\n",
      "k set cration and distribution: 0.011 s\n",
      "Rotating XC potential: 0.255 s\n",
      "Greens function inversion: 393.524 s\n",
      "Calculate energies and magnetic components: 0.702 s\n"
     ]
    }
   ],
   "source": [
    "if rank == root_node:\n",
    "    # iterate over the magnetic entities\n",
    "    for tracker, mag_ent in enumerate(magnetic_entities):\n",
    "        # iterate over the quantization axes\n",
    "        for i, Gii in enumerate(mag_ent[\"Gii\"]):\n",
    "            storage = []\n",
    "            # iterate over the first and second order local perturbations\n",
    "            for Vu1, Vu2 in zip(mag_ent[\"Vu1\"][i], mag_ent[\"Vu2\"][i]):\n",
    "                # The Szunyogh-Lichtenstein formula\n",
    "                traced = np.trace((Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii), axis1=1, axis2=2)\n",
    "                # evaluation of the contour integral\n",
    "                storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))\n",
    "\n",
    "            # fill up the magnetic entities dictionary with the energies\n",
    "            magnetic_entities[tracker][\"energies\"].append(storage)\n",
    "        # convert to np array\n",
    "        magnetic_entities[tracker][\"energies\"] = np.array(\n",
    "            magnetic_entities[tracker][\"energies\"]\n",
    "        )\n",
    "    print(\"Magnetic entities integrated.\")\n",
    "\n",
    "    # iterate over the pairs\n",
    "    for tracker, pair in enumerate(pairs):\n",
    "        # iterate over the quantization axes\n",
    "        for i, (Gij, Gji) in enumerate(zip(pair[\"Gij\"], pair[\"Gji\"])):\n",
    "            site_i = magnetic_entities[pair[\"ai\"]]\n",
    "            site_j = magnetic_entities[pair[\"aj\"]]\n",
    "\n",
    "            storage = []\n",
    "            # iterate over the first order local perturbations in all possible orientations for the two sites\n",
    "            for Vui in site_i[\"Vu1\"][i]:\n",
    "                for Vuj in site_j[\"Vu1\"][i]:\n",
    "                    # The Szunyogh-Lichtenstein formula\n",
    "                    traced = np.trace((Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2)\n",
    "                    # evaluation of the contour integral\n",
    "                    storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))\n",
    "            # fill up the pairs dictionary with the energies\n",
    "            pairs[tracker][\"energies\"].append(storage)\n",
    "        # convert to np array\n",
    "        pairs[tracker][\"energies\"] = np.array(pairs[tracker][\"energies\"])\n",
    "\n",
    "    print(\"Pairs integrated.\")\n",
    "\n",
    "    # calculate magnetic parameters\n",
    "    for pair in pairs:\n",
    "        J_iso, J_S, D, J = calculate_exchange_tensor(pair)\n",
    "        pair[\"J_iso\"] = J_iso * sisl.unit_convert(\"eV\", \"meV\")\n",
    "        pair[\"J_S\"] = J_S * sisl.unit_convert(\"eV\", \"meV\")\n",
    "        pair[\"D\"] = D * sisl.unit_convert(\"eV\", \"meV\")\n",
    "        pair[\"J\"] = J * sisl.unit_convert(\"eV\", \"meV\")\n",
    "\n",
    "    print(\"Magnetic parameters calculated.\")\n",
    "\n",
    "    times[\"end_time\"] = timer()\n",
    "    print(\n",
    "        \"##################################################################### GROGU OUTPUT #############################################################################\"\n",
    "    )\n",
    "\n",
    "    print_parameters(simulation_parameters)\n",
    "    print_atoms_and_pairs(magnetic_entities, pairs)\n",
    "    print_runtime_information(times)\n",
    "\n",
    "    # remove clutter from magnetic entities and pair information\n",
    "    for pair in pairs:\n",
    "        del pair[\"Gij\"]\n",
    "        del pair[\"Gij_tmp\"]\n",
    "        del pair[\"Gji\"]\n",
    "        del pair[\"Gji_tmp\"]\n",
    "    for mag_ent in magnetic_entities:\n",
    "        del mag_ent[\"Gii\"]\n",
    "        del mag_ent[\"Gii_tmp\"]\n",
    "        del mag_ent[\"Vu1\"]\n",
    "        del mag_ent[\"Vu2\"]\n",
    "    # create output dictionary with all the relevant data\n",
    "    results = dict(\n",
    "        parameters=simulation_parameters,\n",
    "        magnetic_entities=magnetic_entities,\n",
    "        pairs=pairs,\n",
    "        runtime=times,\n",
    "    )\n",
    "    # save dictionary\n",
    "    with open(outfile, \"wb\") as output_file:\n",
    "        pickle.dump(results, output_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (3105939143.py, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  Cell \u001b[0;32mIn[111], line 1\u001b[0;36m\u001b[0m\n\u001b[0;31m    ========================================\u001b[0m\n\u001b[0m    ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "========================================\n",
    " \n",
    "Atom Angstrom\n",
    "# Label,        x           y           z          Sx           Sy           Sz         #Q           Lx           Ly           Lz           Jx           Jy           Jz\n",
    "--------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
    "Te1         1.8955      1.0943      13.1698     -0.0000     0.0000      -0.1543    # 5.9345       -0.0000      0.0000       -0.0537      -0.0000      0.0000       -0.2080   \n",
    "Te2         1.8955      1.0943      7.4002      0.0000      -0.0000     -0.1543    # 5.9345       0.0000       -0.0000      -0.0537      0.0000       -0.0000      -0.2080   \n",
    "Ge3         -0.0000     2.1887      10.2850     0.0000      0.0000      -0.1605    # 3.1927       -0.0000      0.0000       0.0012       0.0000       0.0000       -0.1593   \n",
    "Fe4         -0.0000     0.0000      11.6576     0.0001      -0.0001     2.0466     # 8.3044       0.0000       -0.0000      0.1606       0.0001       -0.0001      2.2072    \n",
    "Fe5         -0.0000     0.0000      8.9124      -0.0001     0.0001      2.0466     # 8.3044       -0.0000      0.0000       0.1606       -0.0001      0.0001       2.2072    \n",
    "Fe6         1.8955      1.0944      10.2850     0.0000      0.0000      1.5824     # 8.3296       -0.0000      -0.0000      0.0520       -0.0000      0.0000       1.6344    \n",
    "==================================================================================================================================\n",
    " \n",
    "Exchange meV\n",
    "--------------------------------------------------------------------------------\n",
    "# at1     at2   i  j  k    #    d (Ang)\n",
    "--------------------------------------------------------------------------------\n",
    "Fe4     Fe5     0  0  0    #    2.7452\n",
    "Isotropic -82.0854\n",
    "DMI 0.12557 -0.00082199  6.9668e-08\n",
    "Symmetric-anisotropy -0.60237    -0.83842 -0.00032278 -1.2166e-05 -3.3923e-05\n",
    "--------------------------------------------------------------------------------\n",
    "Fe4     Fe6     0  0  0    #    2.5835\n",
    "Isotropic -41.9627\n",
    "DMI 1.1205     -1.9532   0.0018386\n",
    "Symmetric-anisotropy 0.26007 -0.00013243     0.12977   -0.069979   -0.042066\n",
    "--------------------------------------------------------------------------------\n"
   ]
  }
 ],
 "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}